基于GF-1卫星数据和P-M方法的灌溉需水量空间格局研究

王志东，程 锐，宋 妮，王景雷，韩其晟，赵西宁

（1.西北农林科技大学 水利与建筑工程学院，陕西 杨凌 712100；2.西北农林科技大学旱区农业水土工程教育部重点实验室，陕西 杨凌 712100；3.中国水利学会，北京100053；4.中国农业科学院 农田灌溉研究所/农业农村部作物需水与调控重点开放实验室，河南 新乡 453002；5.中国科学院 水利部水土保持研究所，陕西 杨凌 712100)

0 引言

【研究意义】随着水资源短缺日益严重，农业粮食生产与农业用水之间的矛盾日益突出[1]，农业用水量高占比的干旱半干旱地区水资源供需矛盾更加明显[2]，因此合理规划农业水土资源是当前的研究热点[3-5]。作物种植结构和作物需水量与灌溉需水量是农业水土资源合理规划的基础资料，但由于宏观研究区域面积广、跨度大，种植结构的空间分布往往难以通过人工调查获取。此外，不同地区光热条件不同，作物需水量和灌水量存在差异，利用定点气象站数据代表特定区域计算其作物需水量时存在相应的误差。准确有效地提取研究区种植结构，分析区域尺度作物需水量和灌溉需水量的时空分布可以为作物种植区域适宜性分析、区域农业水资源时间调度和种植结构调整提供有力依据。【研究进展】随着3S 技术的发展，利用遥感卫星数据获取地物覆盖信息的精度和应用范围不断提高，遥感数据在农业领域的应用日益广泛[6-8]。作物种植结构是农业发展研究的重要指标，但由于下垫面作物类型繁杂，以往单一时相的遥感图像往往无法满足精度要求[9]。常用的Landsat8 卫星数据分辨率为30 m，分类精度往往不能满足研究需求，或者只能提取单一植被[10]。作物需水量和灌溉需水量是灌溉农业的关键性指标，也是农业水土资源规划的主要参数。以往宏观研究中常以单一点的数据来代表整个区域，忽略了各参数的空间差异[11]；在时间尺度上常以年灌溉定额来计算农业需水量[12]，不利于农业水资源的年内分配。【切入点】针对以上问题，本文利用分辨率为16 m 的GF-1 卫星（国产）数据构建作物生育期NDVI时间序列，结合使用1年多时相遥感卫星与实地调研数据实现区域种植结构的提取，利用ArcGIS 软件中的空间插值探究作物需水量与灌溉需水量的空间分布，为区域灌溉制度的制定和农业水土资源的有效配置提供参考。【拟解决的关键问题】基于作物物候期多时相NDVI提取作物种植结构，并使用ArcGIS 软件实现对区域作物需水量和灌溉需水量的时空分布分析。

1 材料与方法

1.1 研究区概况

阿克苏河流域位于我国新疆维吾尔自治区西南部，地处N 40°00'—42°27'，E 75°35'—82°00'之间，具体位置见图1（为保证研究地区完整性，本文将阿拉尔以南的塔里木河流域的部分区域一并纳入研究区，以下统称阿克苏河绿洲区）。流域内，绿洲区占地面积约13 000 km2，是新疆主要灌溉绿洲农业区和重要粮食、水果、干果生产基地之一，也是棉花主要生产区。流域内气候干燥，蒸发量大，降水稀少，晴天多，日照时间长，年均总日照时间2 570.9～2 996.8 h，热量资源丰富。流域年均降雨量64 mm，年均蒸发量1 890 mm，属典型纯灌溉农业区[13]。流域区域主要包括乌什、温宿、阿克苏、阿瓦提、柯坪及农一师阿拉尔市的部分团场。流域内主要粮食作物为冬小麦、夏玉米、春玉米和水稻等，主要经济作物为棉花，林果类作物主要为红枣、核桃、苹果、香梨和杏等。

图1 研究区示意Fig.1 The study area

1.2 数据来源帧

1.2.1 遥感影像数据

根据2020年各月的研究区内GF-1 遥感影像，筛选得到可用图像共16 景，遥感数据基本信息见表1。

表1 GF-1 影像列表Table 1 List of GF-1 images

1.2.2 实地数据调研

实地调研于2020年11月开展，利用手持GPS进行现场标点。共采集研究区内主要地物类型及作物样本点2 541 个，并随机分为分类样本和验证样本，样本数量较少的白萝卜、甘草、辣椒和青储玉米等并作其他作物，桃、葡萄和枸杞等并作其他果树，类别及样本数见表2。根据阿克苏绿洲区作物生育期的相关文献及实地调研结果，阿克苏绿洲区主要作物及果树物候期见表3。调研结果作为NDVI曲线构建结果分析的参考依据。

表2 实采样本数量Table 2 The number of actual samples

表3 研究区主要作物及果树物候期Table 3 Major crop and fruiter phenological period in study area

1.2.3 气象数据

气象数据来源于国家气象信息中心（http://data.cma.cn）发布的阿克苏河绿洲区内及附近4个气象站（沙雅、柯坪、阿拉尔和阿克苏）的逐日平均气温、最低、最高气温、降雨、风速、日照时间等。

1.3 研究方法

参考作物生育期构建研究区作物NDVI曲线，在ENVI 5.3 软件中利用实采样本点依照监督分类对研究区作物种植结构进行提取；根据研究区气象数据和ArcGIS 软件中的空间插值得出研究区作物需水量和灌溉需水量的空间分布。

1.3.1 遥感数据处理

利用ENVI 5.3软件对GF-1卫星遥感数据进行预处理，预处理内容包括正射校正、辐射定标、大气校正、拼接和裁剪等步骤。

1.3.2 构建全年归一化植被指数NDVI时间序列

利用遥感图像红光波段与近红外波段波长计算NDVI，即GF-1 号卫星的第3、第4 波段。计算式为：

式中：ρNIR为近红外波段波长；ρR为红外波段波长。

根据预处理后的GF-1 影像，选取样本点中各作物的纯净像元，采用ENVI 5.3 软件的统计工具构建研究区内主要作物生育期NDVI曲线（图2（a）），并构建各植被纯净像元最小NDVI曲线与水域、荒地和建设用地纯净像元最大NDVI曲线作对比（图2（b））。

图2 NDVI 时间序列Fig.2 NDVI time series

从图2（a）可以看出，作物类与果树类NDVI曲线有明显区别，除冬小麦-夏玉米外，作物类NDVI曲线大多在7、8月出现峰值，9、10月大幅下降，果树类NDVI曲线变化较为平缓。由图2（b）可知，8月各作物最小NDVI值远大于非植被类最大NDVI值，因此取8月水域、荒地和建设用地中NDVI最大值0.09 提取植被，并以此制作掩膜提取文件。

1.3.3 种植结构分类

利用以上植被掩膜数据对构建的NDVI时间序列进行掩膜提取，并进行作物分类识别。经过所有分类样本点分离度的统计，果树与作物之间分类度均大于1.8，各作物之间样本分离度均大于1.7。此外，由于研究区面积较大，遥感影像区域色差及多年生果树树龄存在差异，同一生育期NDVI存在较大差异，各果树之间样本分类度在1.65 左右。波段组合的NDVI时间序列图像，分别选择监督分类中的最大似然法、最小距离法、神经网络法、马氏距离法进行分类。

1.3.4ET0计算及IDW 反距离权重插值

1）Penman-Monteith 公式

本文选用P-M 公式计算绿洲区气象站点2020年逐日参考作物需水量ET0[14]。

式中：Δ 为饱和水汽压曲线斜率（kPa/℃）；γ为干湿表常数（kPa/℃）；U2为距地面2 m 处的风（m/s）；Rn为冠层表面净辐射（MJ/（m2·d））；G为土壤热通量（MJ/（m2·d））；T为日平均气温（℃）；ea为饱和水汽压（kPa）；ed为实际水汽压（kPa）。

2）反距离权重插值法

Inverse Distance Weighted(IDW)是GIS 中简便，常用的空间插值方法，以插值点和样本点间距进行加权平均，离插值点越近的样本权重越大。本文根据气象数据分别对参考作物需水量ET0，降雨量P逐月进行IDW 插值。

本文参考各县灌溉水利用系数，对灌溉水利用系数进行插值，使之更加符合实际情况。

1.3.5 作物灌溉需水量计算

1）作物需水量计算

FAO 中作物腾发量（需水量）可以通过ET0和Kc的乘积来计算。本文基于单作物系数法求解作物需水量，见式（3）。单作物系数法中将作物生育期分为生长初期、快速生长期、生长中期和生长后期4 个阶段[15]，各生育阶段Kc取值不同。本文以作物物候期为依据将Kc按分段单值平均法分配至各月（见表4）。根据前文分类结果，利用ArcGIS 重分类逐月定义空间不同作物Kc值，与流域绿洲区逐月ET0空间插值结果相乘得到流域作物逐月作物需水量空间分布。

表4 Kc 值按月分配Table 4 Allocate Kc values by month

2）逐月作物净灌溉需水量（NIR）

作物净灌溉需水量取作物需水量与月作物生育期内有效降雨量的差值，通过ArcGIS 软件空间分析工具逐月计算绿洲区作物净灌溉需水量空间分布：

式中：ETc为各月作物需水量；P为各月（生育期内）降雨量；0.52 为降水利用系数[15]。

3）全年毛灌溉需水量

毛灌溉需水量计算式为：

式中：W是全年全流域毛灌溉需水量；μ为灌溉水利用系数。

2 结果与分析

2.1 分类结果

由图2（a）NDVI曲线可以看出，单作物分类时，冬小麦在6月成熟收割时夏玉米尚未出苗，NDVI值明显降低，是区分冬小麦、夏玉米和其他作物的最好时期；春玉米于8月底收获后，9月的NDVI值明显低于棉花和水稻；水稻于9月进入黄熟期后，NDVI值低于棉花，是区分水稻和棉花的关键时期。单林果分类时，由于杏树种植密度低于其他果树类，整个生育期NDVI值均低于其他果树类，便于区分杏树类；而核桃、苹果和香梨的树冠较大，绿度值高于红枣，因此它们的NDVI值显著大于红枣，便于将核桃及苹果和香梨与红枣区分开。

采用不同监督分类方法（最大似然、最小距离、神经网络、马氏距离）对遥感地物进行识别，并利用验证样本求各分类方法总体精度和Kappa系数。结果如表5所示。可以看出，最大似然法分类精度明显高于其他3 类，因而选择最大似然法对阿克苏河流域绿洲区耕地结构进行分类，（图3），并建立精度混淆矩阵验证不同类别用户精度和制图精度（表6）。

表 5 监督分类方法比较Table 5 Comparison of supervised classification methods

图3 最大似然法分类结果Fig.3 Classification results by maximum likelihood method

表6 最大似然法验证及分类结果Table 6 Maximum likelihood method verification and classification results

分类结果见表6，耕地总面积为7 563.54 km2。由于2020年统计数据尚未公布，分析研究区2013—2018年耕地数据得知耕地总量变化不大，对比分类结果耕地总面积与2018年阿克苏地区年鉴耕地面积，耕地总体提取精度达99.98%（7 563.54/7 734.99）[13]。由表6 和图3 可以看出，研究区粮食作物面积占比约20.74%，主要分布在流域中上游，其中冬小麦夏玉米主要分布在乌什和温宿县，根据实地调研结果，乌什县存在部分冬小麦和核桃、红枣、杏树等的套种（标样本点时，新种植果树标为冬小麦，往年果树则标为果树）；水稻主要分布在温宿县和农一师六团水稻连；玉米主要分布在乌什、温宿、阿克苏等流域上游地区。经济作物棉花面积占比达到30.08%，主要分布在阿克苏市、阿拉尔市、农一师一团、二团、五团南部及阿瓦提南部等流域中下游。果树占比约30.42%，整个流域均有分布。其中红枣主要分布在流域下游的阿瓦提、阿拉尔；核桃主要在阿瓦提、阿克苏；苹果香梨主要分布在温宿和阿拉尔。

2.2 ET0 计算结果

利用P-M 公式计算各气象站全年逐日，逐月ET0变化曲线见图4，阿克苏流域全年参考作物需水量主要在5—8月，占全年ET0的60%左右。空间分布上，柯坪站由于海拔高，温度低，ET0低于其他3 站。

图4 绿洲区四气象站逐日逐月ET0Fig.4 The ET0 of four weather stations in Green Region day and month

2.3 作物需水量与灌溉需水量空间分布

如图5（a）所示，区域作物需水量ETc在空间上存在较大差异。温宿县和阿克苏市沿河地区是主要核桃种植区，年参考作物需水量ET0明显大于其他地区，而夏玉米主要种植区的乌什县ETc相对较小。同一种作物不同种植区域需水量也存在差异，对于区域最重要的经济作物棉花，需水量也存在较大差异，ETc范围在704.4～801.8 mm 之间，阿克苏市和温宿县附近的棉花种植区ETc明显高于其他地区。

在各月ETc的基础上，逐月对区域降雨量进行插值、叠加得到全年降雨量。如图5（b）所示，整个绿洲区降雨量偏少，从西北向东南呈减小趋势，尤其是阿拉尔部分地区全年降雨量仅在9 mm 左右，灌溉仍然是保证作物生长的主要水分来源。本文以各月作物需水量为基础，结合区域灌溉水利用系数空间插值结果（图6（a）），得出流域绿洲区毛灌溉水量的时空分布情况见图6（b）和图7。阿拉尔市农一师农业机械化程度高于绿洲区其他地区，灌溉系统配置也更加完善，因而灌溉水利用系数高于区域其他地区。区域全年毛灌溉需水量与全年作物需水量时空分布规律相似，春夏2 季占全年需水量和毛灌溉需水量的70%左右。阿拉尔虽然降雨略少，但由于灌溉水利用系数较高，棉花，红枣等灌溉需水量均低于其他同类作物种植区。

图5 作物需水量与全年降雨量空间分布Fig.5 Spatial distribution of crop water demand and annual rainfall

图6 灌溉水利用系数与全年毛灌溉量空间分布Fig.6 Spatial distribution of irrigation water utilization coefficient and annual gross irrigation amount

图7 毛灌溉需水量各季节空间分布Fig.7 Spatial distribution of irrigation water demand in each season

3 讨论

GF-1 卫星于2013年4月26日在酒泉卫星发射中心成功发射，具有同一地区的短时多次拍照的优势，其中WFV 数据分辨率达16 m。本文采用多时相GF-1号卫星数据构建NDVI曲线对阿克苏河绿洲区提取作物种植结构。虽然区域面积跨度广，不同时相卫星图像拼接存在色差，但由于样本点数量大、分布均匀，监督分类法精度基本满足分类要求。不同监督分类方法中，最大似然法应用范围广，精度明显高于其他分类方法[16-17]，这与以往区域尺度监督分类结果一致[18]。根据阿克苏地区和新疆建设兵团年鉴[19-20]，结合实地调研发现，阿合苏河绿洲区棉花主要分布在阿克苏市与农一师部分团场，冬小麦主要分布在乌什县和温宿县，水稻主要分布在温宿县，除红枣外的果树类主要分布在阿克苏、阿瓦提和温宿，与本文分类结果基本吻合。不同作物NDVI曲线规律存在显著差异，与以往研究结果相似[21-22]。本项研究分类精度较高主要是因为构建了作物生育期的NDVI时间序列。

ET0曲线显示，阿克苏河流域绿洲区主要作物需水时期在5—8月，柯坪站海拔较高，温度低，ET0明显低于其他3 站。采用空间插值法逐月对区域ET0进行插值，可以更准确地反映出区域ET0的空间变化规律，相比以往直接以站点数据代替某一区域[23-24]更符合实际。通过对有效降雨和灌溉水利用系数的空间插值得出绿洲区毛灌溉需水量四季的时空分布情况（图6（b）、图7），全年毛灌溉需水量在乌什、温宿等流域上游灌区明显较高。对于主要经济作物棉花，流域下游降雨量和灌溉水利用系数存在差异，即使同是棉花主要种植区的阿克苏市和阿拉尔市，毛灌溉需水量也存在差异，这与王梅等[25]研究结果相似。结合阿克苏地区节水背景，可考虑将部分阿克苏市附近地区的棉花种植区域南移，这也符合集中种植便于管理的要求。本可以为灌溉制度的制定和农业水资源利用的时空格局优化提供一定的依据。

4 结论

1）建立多时相NDVI序列，通过实地采取样本点与监督分类结合的方法对阿克苏河流域绿洲区域主要作物和果树进行解译，最大似然法分类结果最优，通过验证样本混淆矩阵验证总体分类精度达到93.08%，Kappa系数达到0.913，表明利用本文的方法可实现区域尺度作物种植结构的提取。

2）年毛灌溉需水量随河流流向从上游往下游呈减小趋势。温宿县由于主要作物为果树和水稻，灌溉需水量在全年都比较高；时间尺度上，灌溉需水量明显冬季小于其他季节，夏季最大。同时在夏季除乌什县部分地区外，全流域灌区都处在高需水状态，因此在夏季应加大中下游径流输送量以满足灌溉需求。对比潜在（净）灌水量和实际（毛）灌水量可以看出，提高灌溉水利用系数对提高水资源利用效率至关重要。在保证渠道输水效率的基础上，应促进传统灌溉农业向高效节水灌溉农业发展，推广膜下滴灌、喷灌等高效节水的灌溉方式，以提灌溉水利用系数，最终实现水资源的高效利用。