基于光学和雷达遥感信息的灌溉信号分析及灌溉面积提取方法研究—以华北平原灵寿县磁右灌区为例

杨永民，顾 涛，吴 迪，龙爱华，邱建秀，刘宏鑫

(1.中国水利水电科学研究院，北京 100038；2.水利部防洪抗旱减灾工程技术研究中心，北京 100038；3.中国灌溉排水发展中心，北京 100054；4.中山大学 地理科学与规划学院，广东 广州 510275)

1 研究背景

人类活动的水文效应一直是国内外关注的焦点之一，人类活动特别是灌溉用水是影响陆面水文循环最难刻画的部分[1-2]。目前，绝大部分的水文模型对灌溉事件的刻画仍然采用预设的土壤水分消退阈值法或灌溉制度进行参数化，存在较大的不确定性[1，3]。国内外已公布的全球和区域灌溉面积资料存在较大的不确定性，难以反映灌溉时间和次灌溉区域的分布[4-5]。因此，准确的提取区域实际灌溉面积是进一步深入评估人类取用水活动对区域陆地水循环影响研究的关键前提。灌溉面积是灌区有效实施用水总量控制和定额管理的重要用水参数之一，科学准确的灌溉面积调查数据是实施农业水资源管理、农业用水效率评估、农业节水、灌区现代化管理、流域水资源管理等应用领域的重要基础[6-10]。目前，区域灌溉面积主要通过调查统计方式获取，时效性和准确性较差，且调查和统计工作量较大。因此，传统灌溉面积获取方法已不能满足灌区现代化管理需求，亟需建立更为科学、合理的技术方法，尤其在基于遥感等多源信息的高分辨率实际灌溉面积快速提取方面亟待加强研究，为农业用水管理提供科学、客观依据。

遥感作为地面观测的重要手段，为区域实际灌溉面积的动态监测提供了重要的技术手段。现有的基于卫星遥感的灌溉面积调查方法主要包括以下几类：(1)植被指数阈值法，基于归一化植被指数或垂直植被指数方法使用阈值法进行灌溉探测[11-12]。目前大部分的研究工作主要集中在可见光和近红外谱段范围，由于光学遥感指数对土壤水分的敏感性较差，此类方法的不确定性较大[13]。(2)遥感地表温度法，基于地表温度和下垫面土壤水分之间的关系来分析识别灌溉活动[14-15]。但受天气因素和目前热红外卫星遥感分辨率的限制，实现连续、动态监测灌溉时空分布的难度较大[16]。(3)结合遥感地表蒸散和土壤水分的大范围灌溉探测法，此类方法空间分辨率较粗，主要应用于大范围的灌溉用水分析计算[17-19]。区域实际灌溉面积的动态监测对区域水资源管理等方面的应用至关重要。目前，传统的基于光学植被指数阈值法和遥感地表温度的灌溉面积提取方法对灌溉活动导致的土壤和植被含水量的变化响应较差，提取的灌溉范围及面积存在较大的不确定性。此外，传统的基于光学遥感的灌溉面积探测方法本质上依赖于土地利用分类，大部分的工作是在耕地提取的基础上进行的简单分析和判别，并不能有效地揭示灌溉事件的时空格局，也难以对不同水源的灌溉区域进行有效区分。

为此，本文开展基于雷达遥感信息的灌溉信号分析及次灌面积提取方法的研究，以期为农业用水管理等方面应用提供科学、客观依据。本文首先使用水云模型结合实测土壤水分数据进行正向模拟，以量化灌溉前后土壤水分变化导致C波段雷达后向散射系数的变化特征。然后，提出一种基于时序差值和局部阈值法的实际灌溉面积提取方法，使用高频的哨兵1号雷达遥感观测数据开展华北平原灵寿县磁右灌区灌溉面积的提取，并结合灌区实地调查资料对提取结果进行验证。

2 研究区域及数据

2.1 研究区概况本文以河北省灵寿县磁右灌区南部为研究区。灵寿县磁右灌区位于太行山中部东麓，南距省会石家庄市30 km，东与正定、行唐县城毗邻。灌区范围内有耕地面积33.17万亩，占灵寿县耕地面积的92%，总播种面积 53.68万亩，设计灌溉面积26万亩，有效灌溉面积25.6万亩。由于灌区渠系及建筑物年久失修，近几年横山岭水库和燕川水库实灌面积在8万亩左右，另有机井灌溉6万多亩，黄壁庄水库提水灌溉近万亩，其余均为旱作种植。本文选择其中的地形平坦区域为研究区(图1)，以减轻地形因素对雷达信号的影响。灌区种植作物以小麦、玉米为主。

图1 研究区概况

磁右灌区属暖温带大陆性季风气候，四季分明，雨热同季。多年平均气温12.8 ℃，最低气温-24.2 ℃，最高气温41.1 ℃。多年平均日照时数2661 h，多年平均无霜期212 d。多年平均降雨量为 618.3 mm，降雨多集中在6—9月，占全年的75%以上。降雨时空分布不均，冬春干旱少雨，夏秋洪涝多灾。

2.2 数据收集为分析磁右灌区的雷达后向散射系数的变化特征，本文收集了2021年3月17日、3月29日、4月10日、4月22日、5月4日、5月16日共计6期哨兵1号雷达卫星数据和2021年3月23日哨兵2号光学卫星遥感数据。对收集的雷达遥感数据进行斑点滤波、辐射定标、地理编码、后向散射系数时序滤波和影像裁剪等处理。为开展灌溉信号分析，在研究区选择4个片区，分别为青同片区、台南片区、大干斗片区和桃园片区(4个片区具体位置已在图1中标注)，各个片区的面积分别为1818亩、518亩、2749亩、6552亩。其中青同片区地块较为平整，台南和大干斗片区地势有一定的起伏，这三个片区都以渠灌和管灌为主。桃园片区以桃树和林果为主，主要采用滴灌方式。2021年春灌期间，分别于2021年3月23日和4月23日赴磁右灌区开展两次灌溉实地调研，并进行灌区灌溉地块和非灌溉地块的信息采集(如图1所示)，实地采集的灌溉地块信息将用于分析和评估遥感探测灌溉信号的能力和精度。为开展模型参数的率定，本文还收集了2019和2020年度河北省深州市灌溉试验站大田观测的土壤水分数据和同期的哨兵1号雷达卫星数据。其中，土壤水分观测频率为1 h，大田作物为小麦和玉米。

3 模型方法

3.1 水云模型水云模型是由Attema和Ulaby基于辐射传输方程提出针对植被覆盖区的微波散射模型[20]。水云模型在植被覆盖区域的微波辐射传输建模和土壤水分反演中有着广泛的应用[21]。水云模型假定植被层为影响微波辐射传输的水云，忽略植被与土壤表面的多次散射，模型中的变量为植被参数和土壤含水量。在给定入射角下，水云模型可表达为：

(1)

τ2=exp(-2BmVsecθ)

(2)

(3)

(4)

为优化和校准水云模型的4个关键参数(A、B、C和D)，本文基于收集的河北省深州市灌溉试验站2019年度大田观测的土壤水分数据和同期的哨兵1号雷达卫星时序数据，使用如下目标函数进行参数率定：

(5)

本文使用SciPy算法库中的非线性最小二乘优化算法进行水云模型4个关键参数的优化和校准。其中，水云模型的4个关键参数的初始值参考Dabrowska-Zielinska等[20]的文献设定，优化后的模型参数估计值见表1所示。

表1 水云模型参数估计值

小麦生育期水云模型参数优化校准前后模型模拟值和哨兵1号雷达卫星观测的后向散射系数的对比见图2所示，优化后的水云模型的模拟值和卫星观测值相关系数由优化前的0.58提升到0.88，优化后的水云模型和观测值拟合较好。本文进一步使用优化的水云模型对小麦返青期和拔节期的灌溉事件前后雷达后向散射系数进行模拟，分析量化灌溉事件后向散射系数的变化特征，为基于雷达的区域灌溉面积提取提供支撑。

图2 水云模型参数校准前、后模拟值与雷达卫星后向散射系数的对比

3.2 基于光学和雷达信息的实际灌溉面积提取方法雷达后向散射系数对土壤和植被含水量较为敏感，而灌溉事件会导致地块的土壤含水量和植被含水量显著增加，因此可基于雷达后向散射系数的时序变化对灌溉信号进行探测，对实际灌溉面积进行提取。对此，本文提出一种基于时序差值和局部阈值法的实际灌溉面积提取方法，主要包括以下步骤：

(1)基于全局阈值的分割计算。基于相邻时序的雷达卫星遥感数据计算雷达后向散射系数差值，基于全局阈值对潜在灌溉地块进行分割计算：

ΔVV=VVt2-VVt1

(6)

Irrig1=ΔVV>Th

(7)

式中：ΔVV为时段VV极化后向散射系数的差值；VVt2和VVt1分别为t2和t1时刻VV极化后向散射系数值；Th为全局分割阈值；Irrig1为基于全局阈值分割的潜在灌溉区域。全局分割阈值可基于3.1节水云模型对灌溉事件的量化模拟分析确定，本文中的全局分割阈值为1 dB。

(2)局部阈值分割的计算。为了抑制裸土地表粗糙度变化对时序后向散射系数产生的影响，引入植被指数对作物种植区域进行增强：

ΔEVV=ΔVV·NDVI

(8)

与大范围的降水导致的区域后向散射系数的变化特征不同，灌溉事件导致雷达后向散射系数呈现明显的局地变化特征，可使用局部阈值法对潜在灌溉地块进行分割计算。首先，将研究区划分为L×L的n个网格。由于灵寿县磁右灌区地块多为矩形的畦田，畦长为20～60 m，畦宽为2～10 m。本文结合研究区地块特征信息，选用畦长的30～50倍作为网格大小，在试验分析的基础上确定L为2 km。针对每个网格内的图像数据进行阈值分割计算：

Irrig2，i=ΔEVV>Thi

(9)

然后，合并n个网格计算结果得到区域的分割计算结果：

(10)

式中：Irrig2，i为第i个网格的潜在灌溉区；Thi为第i个网格内ΔEVV的均值；Irrig2为研究区基于局部阈值分割的潜在灌溉区。

(3)基于作物分类和灌溉制度的潜在灌溉区判别。受到雷达噪声等因素的影响，基于全局和局部阈值分割提取的潜在灌溉区仍存在一定的不确定性。对此，可进一步结合研究区的作物分类和灌溉制度对潜在的灌溉区进行判别：

Irrig3=Cropm

(11)

式中：Cropm为研究区常年同期灌溉的作物类型；Irrig3为基于灌溉制度的研究区潜在灌溉区。最后，对以上三个潜在灌溉区进行综合，获得最终的区域实际灌溉范围：

Irrig=Irrig1∩Irrig2∩Irrig3

(12)

4 结果与分析

4.1 基于水云模型的灌溉事件模拟分析本文使用深州试验站2020年小麦地块春季灌溉期观测的土壤含水量和叶面积指数数据，基于水云模型对灌溉前后的雷达后向散射系数进行正向模拟分析，如图3所示。由图可知，在3月15日的灌溉事件后，10 cm深度的土壤含水量由0.18 cm3/cm3增至0.31 cm3/cm3，土壤含水量上升明显。3月17日以后土壤含水量呈现明显的消退趋势，到4月5日降至0.23 cm3/cm3。期间叶面积指数呈现持续上升趋势，由3月初的1 m2/m2逐渐增到3 m2/m2左右。模拟显示雷达后向散射系数对灌溉事件导致的土壤含水量变化较为敏感，灌溉事件导致模拟的VV和VH极化的雷达后向散射系数分别增加了1.51 dB和1.14 dB，其中VV极化的响应更为显著。同时，与灌溉事件后土壤水分的消退相对应，3月17日以后模拟的VV和VH极化的雷达后向散射系数都呈现缓慢下降趋势。由水云模型的正向模拟可知雷达后向散射系数对灌溉事件所导致的土壤含水量的变化有较好的响应，VV极化对灌溉事件的响应比VH极化敏感。基于模拟分析可知，在华北平原小麦春灌时期，灌溉事件导致的后向散射系数变化超过1 dB，而哨兵1号雷达卫星的绝对辐射精度优于1 dB[22]，灌溉事件有望通过雷达卫星进行探测。

图3 灌溉事件后向散射系数的变化

基于雷达后向散射系数的模拟结果可进一步分析灌溉事件探测的可行性。哨兵1号卫星在我国华北平原的重访周期为12 d。为此，基于相隔12 d的时序雷达后向散射系数差值(ΔVV12)分析灌溉事件探测的最佳时期。由图4可知，在灌溉前雷达后向散射系数基本保持稳定，ΔVV12在0附近。而在灌溉后的10 d内ΔVV12都超过1 dB。灌溉事件后，随着土壤含水量的消退，雷达后向散射系数从峰值开始降低。在灌溉事件发生的13 d之后，ΔVV12变为负值，并随着时间的推移不断增大。由以上分析可知，使用哨兵1号雷达时序后向散射系数的差值可以用于探测灌溉事件，特别是在灌溉事件发生后10 d内，相隔12 d的后向散射系数的差值显著。此外，灌溉土壤水分的消退特征也可被后向散射系数的时序差值所反映。

图4 雷达后向散射系数12 d时序差值

4.2 区域灌溉信号探测分析使用收集的哨兵1号雷达时序遥感数据进行区域灌溉信号的探测分析，结果如图4所示。图5(a)为区域的哨兵2号光学影像数据，用作参考和对比。3月17日至3月29灌区开展2021年度春季第一轮次的渠水灌溉。图5(b)为3月17日与3月29日的雷达后向散射系数时序差值，由图可见青同片区周边和台南片区的后向散射系数增加显著，这与片区灌溉活动导致土壤含水量的增加密切相关。图5(c)为3月29日与4月10的雷达后向散射系数时序差值。期间灌区已结束第一轮次渠水灌溉，青同等片区的雷达后向散射系数呈现明显的减小趋势，这与灌溉后土壤水分消退密切相关。但在青同片区以东和以南区域雷达后向散射系数增加明显，这是由于使用地下水的井灌区域仍在进行第一轮次的春季灌溉所致。

图5 基于时序雷达后向散射系数的区域灌溉信号探测分析

4月16日至4月26日灌区开展2021年度春季第二轮次的渠水灌溉。图5(d)为4月10日和4月22日的雷达后向散射系数时序差值图。由图可知，在青同、大干斗片区以及研究区西部地区雷达后向散射系数呈现明显的增加趋势，而在靠近河流的井灌区域雷达后向散射系数呈现明显的减小趋势。图5(e)为4月22日和5月4日的雷达后向散射系数的时序差值图，青同等渠水灌溉区域的雷达后向散射系数随着土壤水分的消退呈现减小趋势，而青同片区东南部井灌区域雷达后向散射系数呈现增加趋势。图5(f)为5月4日和5月16日的雷达后向散射系数的时序差值图，灌区大部分区域雷达后向散射系数都呈现增加趋势，这是由于5月14日灵寿县普降超过30 mm的降水所致，因降水丰沛，该灌区未开展春季第三轮次的灌溉。

基于以上分析可知，雷达后向散射系数随着灌溉事件土壤水分的显著增加和消退变化呈现明显的先增大后减小的时域特征。灵寿县磁右灌区的渠灌时间较为集中，每次灌溉持续时间在10 d左右，而井灌持续事件更长，由于井灌区和渠灌区灌溉时间的差异，可在雷达后向散射系数的增减变化的空间格局上得到反映。此外，与降水事件导致的大范围后向散射系数变化不同，灌溉事件导致雷达后向散射系数呈现明显的局地变化特征。

4.3 区域实际灌溉面积提取研究区作物分类是区域潜在灌溉区判别的基础，本文基于哨兵2号卫星数据，结合实地调研采集的训练样本，训练随机森林分类算法实现研究区作物分类提取。研究区的土地利用分类结果如图6(a)所示。使用分层抽样法提取110个评估样本点对分类结果进行评估，结果如表2所示，总体分类精度为90%，说明分类结果可靠，可进一步用于潜在灌溉区域的判别。研究区的主要作物以小麦为主，面积约为15.2万亩。

表2 研究区作物分类精度评估

图6 磁右灌区作物种植结构及灌溉面积提取

基于本文提出的提取方法可对区域的实际灌溉面积进行提取，如图6所示。图6(b)为3月17日至3月29日小麦返青期第一次渠灌范围提取结果，灌溉面积约为9.20万亩。图6(c)为4月10日至4月22日小麦拔节期第二次渠灌范围提取结果，灌溉面积约为10.64万亩。使用实地调查采集的样本点数据进行验证(见表3所示)，灌溉地块的提取精度为76.9%，非灌溉地块的提取精度为76.2%，总体的提取精度为76.6%。磁右灌区井灌区的灌溉时间晚于渠灌时间，且灌溉周期较长，井灌区和渠灌区的雷达后向散射系数存在差异。使用两次灌溉周期中3月29日和4月10日的雷达数据，使用灌溉探测方法对井灌范围进行提取，结果如图6(d)所示。井灌面积约为3.9万亩。基于遥感提取结果显示井灌区域主要集中在研究区的东南部和靠近河谷的区域，这与实地调研情况及灵寿县灌溉区井灌区的范围基本一致。

表3 灌溉面积精度评估

基于提取的灌溉面积和片区灌溉用水计量可进一步分析灌溉片区的用水状况(见表4所示)，青同片区一水灌溉面积为1226.4亩，亩均用水量为90 m3，灌溉水深为134.9 mm，二水灌溉面积为1047.7亩，亩均用水量为108.9 m3，灌溉水深为163.3 mm。大干斗片区一水灌溉面积为939.7亩，亩均用水量为70.9 m3，灌溉水深为106.3 mm，二水灌溉面积为722亩，亩均用水量为98.9 m3，灌溉水深为148.3 mm。青同和大干斗片区两次灌溉的亩均用水量分别为98 m3和83 m3。青同片区的亩均用水量略大于大干斗片区，磁右灌区小麦的设计灌溉定额为65 m3，实际测算表明两个片区的实际亩均用水量都高于设计的灌溉定额。

表4 灌溉用水测算分析

5 结论与展望

实际灌溉面积是区域水土资源利用的关键指标，灌溉面积监测是农业水资源利用和流域水资源管理等应用的基础。针对目前基于区域灌溉探测识别方法存在的不足，本文使用水云模型结合实测土壤水分数据进行模拟，量化灌溉事件导致后向散射系数的变化，构建一种基于时序差值和局部阈值法的实际灌溉面积提取方法，使用高频次的哨兵1号雷达卫星遥感数据开展华北平原灵寿县磁右灌区灌溉面积的提取，并结合灌区实地调查资料对提取结果进行验证。本文结论如下：

(1)使用水云模型模拟显示春季小麦地块的灌溉事件可导致C波段的雷达后向散射系数变化超过1 dB，灌溉信号可被绝对辐射精度优于1 dB的雷达卫星探测。

(2)与大范围的降水事件导致的大范围后向散射系数变化不同，灌溉事件导致雷达后向散射系数呈现明显的局地变化特征，雷达后向散射系数的时序变化特征对磁右灌区的渠水灌溉事件有较好的响应。此外，由于井灌区和渠灌区存在明显的灌溉时间差异，可通过雷达后向散射系数的增减变化得到反映。但是，部分农事活动导致下垫面粗糙度的显著变化也会对雷达的后向散射系数产生显著影响。

(3)提出一种基于时序差值和局部阈值法的实际灌溉面积提取方法，基于高频次的雷达遥感信息可对灌区的灌溉事件进行探测。本文提出的方法可为区域次灌溉面积调查和流域水资源管理等应用提供技术支撑，结合实际灌溉用水数据和实际灌溉面积调查可调查分析区域的实际用水状况，进一步评估灌溉定额的合理性和实用性。但本方法基于全局和局部阈值分割方法，阈值的确定会影响灌溉区域的提取。本文采用水云模型模拟分析确定全局分割阈值，局部阈值分割仍然需要结合区域地块特征信息进行试验确定。本方法在平坦地区的效果较好，但应用于复杂下垫面的农田区域，仍然有待进一步的完善。