基于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）年毛灌溉需水量隨河流流向從上游往下游呈減小趨勢。溫宿縣由于主要作物為果樹和水稻，灌溉需水量在全年都比較高；時間尺度上，灌溉需水量明顯冬季小于其他季節，夏季最大。同時在夏季除烏什縣部分地區外，全流域灌區都處在高需水狀態，因此在夏季應加大中下游徑流輸送量以滿足灌溉需求。對比潛在（凈）灌水量和實際（毛）灌水量可以看出，提高灌溉水利用系數對提高水資源利用效率至關重要。在保證渠道輸水效率的基礎上，應促進傳統灌溉農業向高效節水灌溉農業發展，推廣膜下滴灌、噴灌等高效節水的灌溉方式，以提灌溉水利用系數，最終實現水資源的高效利用。