基于GF-1衛星數據監測灌區灌溉面積方法研究
——以東雷二期抽黃灌區為例

宋文龍，李 萌，3，路京選，盧奕竹，史楊軍，賀海川

（1.中國水利水電科學研究院，北京 100038；2.水利部防洪抗旱減災工程技術研究中心，北京 100038；3.首都師范大學 資源環境與旅游學院，北京 100048；4.渭南市東雷二期抽黃工程管理局，陜西 渭南 714000）

1 研究背景

準確監測灌溉面積對于預防干旱、優化灌區水資源利用和水資源調控配置具有重要意義。據統計，世界上約有農田1.5 億～1.7 億hm2，約占地球陸地總面積的17%［1-2］，農業用水量約占人類用水的92%［3］。作為農業大國，我國灌溉面積約占總耕地面積的50%，在全球氣候變暖、干旱等極端事件增加、經濟發展用水需求猛增情況下，水資源短缺形式日趨嚴峻。準確監測灌溉面積時空分布及其變化，是掌握灌溉用水量、灌溉用水去向的關鍵，對監測灌溉系統運轉、預防干旱、優化作物種植結構、提高灌溉用水效率、優化灌區水資源配置均具有重要意義。

國內外農業灌溉面積監測研究在決策樹算法［4-5］、時空螺旋曲線與變化矢量分析［6］、物候特征［7］、支持向量機［8］、自動農田算法［9］、光譜匹配［10］、k-means 和ISOCLASS 算法［11］、自動機器學習算法［12-13］、植被指數與干旱指數反演［14］等方法方面取得諸多進展。其中，光譜匹配技術（SMTs）能夠提取長時間序列多期遙感影像的光譜反射率或特征指數變化特征，在周期性農業灌溉面積監測方面具有獨特優勢。Thenkabai 最先應用這種方法開展土地類型識別與灌溉面積識別研究［10］。2006年，世界水資源管理研究所（IWMI）通過使用多種數據源組成的具有巨大數據量的數據集，根據降雨量與高程信息對全球范圍區域進行分類，并使用非監督分類對區域進行聚類，而后使用光譜匹配技術對全球范圍內的土地利用類型、灌溉面積分布進行了識別，發布了世界第一份1 km 尺度全球灌溉區域圖（GIAM）［15］。2012年，IWMI又針對非洲和亞洲進行了重點研究，繪制了250 m 空間分辨率的灌溉與雨養區分布圖［16］。2016年，Ambika 將這種方法應用到了印度，使用了250 m 空間分辨率的MODIS NDVI 數據（MOD13Q1），先通過匹配計算分類，再使用決策樹通過NDVI 閾值劃分灌溉區域，獲得了印度2000—2015年的灌溉區域分布圖［17］，相比IWMI 250 m 空間分辨率的灌溉產品提升了精度。Teluguntla 結合定量光譜匹配技術（QSMTs）與農田自動分類算法，使用連續14年的250 m 空間分辨率MODIS NDVI 產品數據對澳大利亞的灌溉信息與種植強度進行了提取［18］。

我國具有氣候多樣、農業種植結構復雜、田塊破碎、灌區信息化水平不高等國情，農業統計缺乏空間信息，可公開獲取的灌溉遙感產品缺乏且空間分辨率與精度難以滿足我國實際需求［19］，我國灌溉面積遙感監測對數據空間分辨率要求更高，灌溉面積遙感監測技術需要取得進一步突破。本文將針對我國國情，應用較高分辨率的國產衛星遙感數據，通過光譜匹配方法像元尺度應用，并引入自適應閾值算法，構建高分辨率灌溉面積遙感監測新方法，對我國西北干旱半干旱區典型渠灌灌區開展灌溉面積遙感監測研究。

2 研究方法與數據來源

2.1 研究區域概況東雷二期抽黃灌區（見圖1），位于東經109°10′—110°10′，北緯34°41′—35°00′，橫跨陜西省富平、大荔、蒲城三縣，屬黃土波狀臺原區，海拔介于385～600 m，地勢東南低，西北高。灌區屬溫帶大陸性季風氣候，春季溫暖干旱，夏季多暴雨，并有短暫大風，冬季嚴寒干燥，雨雪稀少，年平均降水量519～552 mm，干旱年份僅為360～390 mm，且降水時空分配不均衡，50%以上雨量集中在7—9月，年蒸發量約為1700～2000 mm。灌區灌溉方式為抽取黃河水進行渠道引水漫灌。設計灌溉面積8.5 萬hm2，分春灌、夏灌和冬灌。主要作物為小麥和玉米，占灌區種植總面積的90%以上。主要作物生長周期及灌溉時間詳見表1。

圖1 東雷二期抽黃灌區示意圖

表1 主要作物生長周期及灌溉時間

2.2 數據

2.2.1 衛星遙感數據 數據源為2018年GF-1 號16 m 多光譜數據，來源為中國資源衛星公共網站平臺（http：//www.cresda.com/CN/），詳細參數見表2。使用2018年研究區云覆蓋量低于15%的影像，并編譯成由92 個波段組成的數據集（共有23 景數據，每景4 個波段）。此數據集用于準確獲取不同類別的時序光譜特征信息［20-21］。

對遙感數據進行輻射定標、幾何校正、大氣校正等預處理，并計算NDVI 得到研究區各像元的NDVI 時間序列曲線。NDVI 能夠較好地反映植被生長狀態，進行歸一化處理后可以更好地限制數據范圍，有利于后續的光譜匹配計算，用于提取灌溉信息。

表2 GF-1 多光譜傳感器信息

2.2.2 實測樣點數據 像元尺度的光譜匹配算法需要端元光譜特征信息，即灌溉區域小麥和玉米的NDVI 時間序列曲線，因此在研究區開展野外調研與采樣工作，獲得研究區主要地物的真實地面樣點，用于提取灌溉面積及精度驗證。

結合歷史灌溉情況與衛星遙感影像，選擇具有典型代表意義的灌溉區域與非灌溉區域，在選取灌溉區域時應注意避免附近有建筑物與道路的干擾，避開農田的邊界處，以防在取點時樣本點中同時包含農田與其他干擾地物。最終在研究區內共選取225 個真實樣點數據（見圖2），主要記錄樣點位置、作物類型及灌溉情況等信息。

圖2 樣點分布

2.2.3 IWMI 灌溉產品 使用國際水管理研究所（International Water Management Institute，IWMI）的數據產品IWMI-Irrigated Area Map Asia（2000—2010）and Africa（2010）作為驗證參考數據（http：//waterdata.iwmi.org/）。該產品使用250m 空間分辨率MODIS NDVI 合成數據（MOD13Q1），通過ISODATA 算法進行非監督分類，之后依據各種類的NDVI 時間序列曲線確定土地利用類型，得出亞洲和非洲的灌溉區域分布圖。IWMI 灌溉產品將農業區域分為灌溉區和雨養區，并且根據作物種植強度將灌溉區具體劃分成單季作物灌溉區域、雙季作物灌溉區域和連續作物灌溉區域。

2.3 研究方法有別于IWMI 采用的聚類光譜匹配，本文針對高分辨率遙感影像，將光譜匹配方法應用于像元尺度，保證所有像元的時序光譜曲線參與匹配計算，減少聚類過程造成的誤差，并引入OTSU 自適應閾值算法自動確定灌溉面積提取閾值。OTSU 算法的引入可保證每個像元計算相似度信息的同時，自動確定相似度合理變化范圍，對于不同年份與不同數據情況有良好變動適應性，無需人為干預即可提取灌溉區域。基于像元尺度光譜匹配計算的灌溉面積遙感監測方法適用于高分辨率衛星遙感數據，滿足提取小地塊灌溉信息的需求，能夠提高灌溉面積監測結果的精度。

2.3.1 光譜匹配方法 光譜匹配法是一種量化端元光譜（樣本光譜）與目標光譜（待測光譜）相似度的方法。根據光譜匹配原理，光譜匹配計算方法可以分為統計算法、光譜波形特征算法、光譜編碼匹配算法以及特征空間算法。本文使用統計算法和光譜波形特征算法計算光譜相似度，通過三個指標對研究區主要作物端元光譜與目標光譜匹配程度進行定量分析［10］，具體算法如下：

根據地面樣點數據提取端元光譜C=（s1，s2，…，sn），從數據集中獲取目標光譜Xi，j=（b1，b2，…，bn）。

（1）形狀定量。光譜關聯相似度（Spectral Correlation Similarity，SCS）計算公式如下：

式中：Xi，j為目標光譜的NDVI 值；C 為端元光譜的NDVI 值；i，j 分別為第i 行、第j 列；n 為數據集層數；μi，j為目標光譜NDVI 均值；μs為端元光譜NDVI 均值；σX為目標光譜的NDVI 標準偏差，σC為端元光譜的NDVI 標準偏差。

SCS 是衡量端元光譜與目標光譜NDVI 時間序列曲線形狀相似度的指標。SCS 值越高，端元光譜與目標光譜的NDVI 時間序列曲線就越相似。當SCS 值為0 時，相似度最小；當SCS 為1 時，相似度最大，即兩種光譜完全一致。所有SCS 值在0～1 范圍之外的像元即可認為與端元光譜相似度差異過大，直接剔除。

（2）距離定量。歐氏距離（Euclidian Distance Similarity，EDS）是用于衡量端元光譜與目標光譜在光譜空間中距離的指標，在光譜空間中距離越近則越相似。EDS 計算公式如下：

通常歐氏距離越大則代表著兩種光譜差異性越大，反之則差異性越小。為了方便計算將上述公式獲得結果通過歸一化處理，將其標準范圍處理至0 到1。具體公式計算如下所示：

式中，EDSnormal值的范圍為在0～1，它從光譜特征空間距離上測量端元光譜與目標光譜NDVI 時間序列曲線間相關性的大小。M，m 分別代表最大、最小的歐氏距離。通過對EDS 的計算可以把SCS 運算結果中滿足條件像元點進行計算，理想狀態下，0 表示二者完全一致，1 表示二者完全不相關。此外，EDS 受像元數量的影響，一旦有新的像元參與計算，該式中的M 和m 即可能會發生變化，對不同研究對象有較好的變動適應性。同時歐氏距離的局限性在于不同年份氣候情況可能導致NDVI 光譜曲線會發生上下平移，進而會影響判斷結果的一致性。因此，僅使用光譜距離的特征描述往往還不足以較好的完成光譜匹配的工作，需要進一步對端元光譜與目標光譜的相似程度進行量化。

（3）形狀距離綜合定量。光譜相似值（Spectral Similarity Value，SSV）綜合了SCS（形狀定量）與EDS（距離定量）的特點，從NDVI 時間序列曲線形狀和光譜特征空間距離兩方面衡量了端元光譜與目標光譜間的形似度，其計算公式如下：

式中SSV 值越小，光譜之間越相似，其范圍通常介于0～1.414 之間。不同作物之間，NDVI 時間序列曲線的差異較大，SSV 值高；對于同一種作物，灌溉區域的NDVI 時間序列曲線比非灌溉區域具有更高的一致性，SSV 值更低。

2.3.2 OTSU 算法 以實地采樣獲取的灌溉區域玉米和小麥的NDVI 時間序列曲線分別作為端元光譜，利用上述衡量光譜相似度的指標計算研究區每個像元的SSV 值。同類型的光譜相似度越高，SSV值越低；不同類的光譜相似度越低，SSV 值越高。需要確定合理分割閾值，當SSV 值小于該閾值時與端元光譜識別為同一類別。因此，引入OTSU 自適應閾值算法計算SSV 分割閾值來判斷是否為小麥或玉米的灌溉區域，小于該閾值即為灌溉區域，從而識別研究區的灌溉面積空間分布情況，極大減少了人為干預造成的影響。

OTSU 算法是無參量的一種自適應閾值選取方法［22］，該方法可以給定一個臨時閾值，然后計算閾值兩側范圍的像元SSV 的方差值，當方差達到最大所對應的像元SSV 為最佳閾值，即目標區域與背景區域差異最大化，而后將相似度大于或等于該閾值的歸并成一類，小于該閾值的歸并成另一類，得到相應的二值化圖像。具體原理算法如下：

當影像數據量級為L（G=0，1，…，L-1）時，初始閾值為h 將圖像分為目標區域T 與背景區域B 兩部分，pi代表SSV 為i 的像元出現的概率，ni代表SSV 為i 的像元數量，N 代表像元總數，T、B 兩部分概率為：

對應的T、B 兩部分的均值為：

則存在閾值h0即為最佳閾值。

3 結果分析

3.1 主要作物端元光譜通過野外實地采樣，獲得了玉米與小麥灌溉區域的端元光譜（見圖3）。圖3（a）為玉米灌溉區（6月至10月）9 個實測樣本點的NDVI 時間序列曲線，從中篩選5～6 條相似程度最高的曲線作為端元參考光譜，然后再使用端元光譜均值化進行光譜匹配計算。均值化可以保證所獲得的均值光譜在代表該種類光譜特征信息的同時減少數據波動對結果的影響，且進行種類光譜廓線匹配計算時縮小同類別差異性。同理，圖3（b）為小麥灌溉區（10月至6月）共10 個實測樣本點的NDVI 時間序列曲線，提取小麥端元光譜均值并計算。

圖3 玉米與小麥端元光譜

3.2 主要作物灌溉區域相似性參數空間分布根據像元尺度光譜匹配算法得到的主要作物灌溉區域識別相似性參數SCS、EDS 和SSV 分布如圖4所示，研究區內未能灌溉到的地方主要位于灌區最北（區域1）、中部（區域2）與右下位置（區域3）。根據實地調查了解，最北側主要是由于沒有渠系覆蓋；而蒲城附近（區域2）是由于靠近城鎮，在進行城鎮建設時灌溉渠系損壞導致灌溉情況較差；右下位置為大荔系統，以種植果樹為主，主要作物種植面積相對較少。因此提取結果與實際情況相符。

圖4 主要作物灌溉區域相似性參數空間分布

3.3 灌溉區域相似性分割閾值通過光譜匹配方法獲得了研究區主要作物小麥與玉米相似度量化空間信息。進而對獲得的量級分布直方圖使用OTSU 算法得到兩種作物對應的相似度分割閾值（見圖5），結果表明小麥灌溉區最佳分割閾值為0.3985，玉米灌溉區最佳分割閾值為0.4639。

圖5 OTSU 法計算直方圖

3.4 主要作物灌溉區域分布應用小麥和玉米灌溉區域相似性分割閾值，得到主要作物小麥和玉米的灌溉區域空間分布圖（見圖6）。根據作物種植強度情況，研究區內的灌溉區域可分為雙季輪作灌溉區（小麥和玉米輪作）、單季小麥灌溉區和單季玉米灌溉區，以雙季輪作灌溉區為主，單季作物灌溉區主要分布在城鎮周圍與渠系未灌溉區域，由于地塊破碎、種植結構復雜造成。多數非植被覆蓋區域（如城鎮、村莊、渠道、道路等）識別效果較好，但在一些小的建筑物周圍有小部分誤識別區域，是由于混合像元所導致。

對灌區各灌溉子系統（各子系統分布參見圖2）的主要糧食作物小麥和玉米的灌溉面積進行統計（見表3）。各灌溉子系統灌溉面積由大到小依次是流曲、孫鎮、興鎮、荊姚、劉集、蒲城和大荔。除大荔外，其他灌溉子系統均是小麥和玉米輪作灌溉區為主，占比大于50%，荊姚占比達到81.72%。大荔果樹種植面積較大，糧食作物灌溉主要是單季小麥和單季玉米，小麥和玉米輪作灌溉較少。

圖6 主要作物灌溉區域分布圖

表3 各灌溉子系統灌溉面積情況 （單位：hm2）

3.5 精度驗證使用實測樣點數據對提取結果進行混淆矩陣精度驗證（見表4），使用高精度數據提取灌溉面積與種植強度結果的總體精度為88.27%，Kappa 系數為0.8308。對于分布最廣的雙季輪作灌溉區域的提取效果最佳，精度達到90.14%；非灌溉區域的提取精度為89.80%，其主要與單季小麥、單季玉米灌溉區發生混淆，與單小麥灌溉區誤差較大，種類錯分誤差為13.04%；單季玉米灌溉區提取精度為89.47%；單季小麥灌溉區的提取精度相對最差，為78.26%。結果分析發現主要是混合像元對提取結果影響最大，單季小麥與單季玉米灌溉區主要分布在城鎮附近，受人類活動影響，導致種植結構分布復雜且細小散碎，提取難度較大。

表4 GF 基于實測樣點的結果精度評價 （單位：%）

3.6 合理性分析圖7為國際水管理研究所（IWMI）發布的MODIS-250m 全球灌溉區域分布成果的研究區部分，是目前可公開獲得精度較高的灌溉區域分布數據，對其重采樣成16 m 分辨率作為參考數據。東雷二期研究區總面積為99 008.06 hm2，IWMI 研究區灌溉面積結果統計為82 745.37 hm2，使用GF-1 結果統計灌溉面積為81 571.58 hm2，二者計算灌溉面積結果差異為1.42%。對本研究獲得的灌溉區域與IWMI 灌溉區域數據進行空間比較，空間重疊率為87.29%。但二者在不同種植強度灌溉區類型及其分布方面提取結果存在較大差異，據灌區實際調研與本研究結果顯示，研究區主要包括雙季輪作（小麥和玉米輪作）灌溉區域與單季作物（小麥、玉米）灌溉區域，而IWMI 發布的成果僅分為雙季作物灌溉區域與三季/連續作物灌溉區域，本研究成果與實際情況相符，而IWMI 成果有誤。

如圖8所示，本研究成果與IWMI-Irrigated Area Map Asia（2000—2010）相比，在空間分辨率與地物識別方面具有顯著優勢，可以更有效區分出建筑物、道路、渠道與田塊信息；受數據分辨率限制，IWMI 結果無法獲取灌溉田塊分布細節，較多的建筑物與道路存在混淆錯分情況。 與IWMI 結果比較，本研究成果能更好地反映我國灌區實際灌溉情況，一方面是由于采用了更高空間分辨率的遙感數據，另一方面相較于基于聚類的光譜匹配計算方法，基于像元尺度的光譜匹配方法適用于較高空間分辨率的灌溉區域遙感監測研究。

圖7 IWMI 發布的灌溉區域分布

圖8 基于GF-1 的灌溉面積監測結果與IWMI 結果細節對比

4 結論

本研究結合野外調研采樣，應用16 m 較高空間分辨率GF-1 衛星影像，獲取了研究區主要作物的時間序列光譜特征信息，構建了像元尺度光譜匹配方法，并使用OTSU 算法自動獲得灌溉區域分割閾值，獲得主要作物種植強度及其灌溉區域分布，取得的主要結論如下：

（1）2018年東雷二期抽黃灌區灌溉面積為81 571.58 hm2，其中雙季輪作（小麥與玉米輪作）灌溉面積為40 335.88 hm2，單季小麥灌溉面積為15 276.94 hm2，單季玉米灌溉面積為14 059.14 hm2。各灌溉子系統灌溉面積由大到小排序依次是流曲、孫鎮、興鎮、荊姚、劉集、蒲城和大荔。

（2）使用地面調研采樣數據對灌溉面積提取結果進行驗證，總體精度為88.27%，Kappa 系數為0.8308。對雙季輪作（小麥與玉米輪作）灌溉區域與單季玉米灌溉區域提取效果較好，精度分別為90.14%和89.47%；單季小麥灌溉區域提取精度較差，精度為78.26%，易與其他類型地物發生混淆。

（3）本文提出的像元尺度的光譜匹配計算方法，更適用于基于高空間分辨率遙感數據的灌溉信息提取，通過引入自適應閾值算法，最大程度減少人為影響，保證了灌溉信息提取的較高精度。對比IWMI 結果，本研究成果可以更有效區分城鎮、村莊、道路、渠道、田塊等邊界信息，可有效識別小田塊灌溉分布，在作物種植強度及其灌溉區域分布方面更符合實際情況。