基于Sentinel-2時序影像的地塊尺度灌溉耕地提取

王 琰，周亞男，汪順營

（河海大學水文水資源學院，南京 211100）

0 引 言

灌溉是提升糧食產量的最主要條件之一，灌溉農業的單產水平是雨養農業的2.5 倍[1]。此外，灌溉消耗了近80%的人類用水[2]，實際灌溉面積和灌溉空間范圍之間有著很強的不確定性，這也進一步導致了在估計灌溉實際用水時的不確定性。準確掌握灌溉耕地的數量和空間分布等信息，對國內國際的糧食安全、水資源管理甚至氣候變化領域都有重要的現實意義。然而，傳統的灌溉數據統計方法通常要耗費大量的時間和人力，不僅效率不高、精度低，且可獲取的范圍有限。隨著科技的發展，遙感、GIS 等更加高效的技術手段有著覆蓋范圍大、探測周期短、成本低等優點，已廣泛應用于農田地物識別研究中[2]。

在遙感提取灌溉耕地的方向上，學者們的研究大致有以下幾種方法：第一種是基于指數的方法，通過各類指數從土壤濕度、干旱程度等各方面判斷地塊是否進行灌溉。例如Deines 等人[8]開發了兩個新的組合指數AGI 和WGI，嘗試將濕度信息和綠度指數相結合來擴大灌溉狀態的差異，在測試中該指數與傳統指數相比具有更高的重要性，可將該指數應用于其他農業地區灌溉耕地的分類。但基于指數的方法需要依賴大量的現有統計數據用于數學模型分析與驗證，對于一些缺少相關資料的地區難以使用該方法建立指數，因此該方法能夠應用的范圍有限；第二種是基于機器學習的方法，機器學習作為一門源于人工智能和統計學的學科，是當前數據分析領域重點研究方向之一，分類問題及其算法是機器學習的一個重要分支。如Ketchum 等人[9]建立了各類土地覆蓋的地理空間數據庫，基于谷歌地球引擎使用Landsat 衛星圖像以及氣候、氣象、地形數據來訓練隨機森林分類器并預測土地類型，繪制了30 m 分辨率的灌溉地圖IrrMapper。但各類別的機器學習算法都有各自擅長的領域和難以克服的缺陷，沒有一種算法可以解決所有問題，此外，數據降維、特征選擇等因素會對分類算法的發展產生很大的影響，因此在實際應用中，應該結合實際進行比較和選擇適當的分類算法以達到更高的分類精度；第三種是空間分配的方法，Zhu 等人[10]估計出每個像素的灌溉潛力，并提出一種中國灌區的空間分配模型，將格網內灌溉潛力最高的像素識別為灌溉像素，再把所有灌溉像素進行組合以獲得整個研究區的灌溉地圖。這種空間分配的方法主要是通過構建一定的分配規則后將像素識別為灌溉，但其中建立分配規則的依據通常是一些影響灌溉耕地空間分布的特征量，由于這些特征量過于依賴樣本數據，所以在很大程度上影響了結果的準確性；還有一些其他方法，例如朱秀芳等人[11]利用降水、實際蒸散發和潛在蒸散發數據提出了雨養指示線的概念，以此來表征耕地受灌溉可能性的大小，為灌溉耕地制圖提供了指示意義強的灌溉特征參量。楊永民等人[12]使用水云模型實測土壤水分數據，利用散射系數的時序變化探測灌溉信號，提取實際灌溉面積。但這種使用非遙感數據的方法所制作的地圖尺度較小，在制作大范圍灌溉地圖方面有所欠缺。

綜上所述，學者們研究所用的數據大多是中低分辨率的多源衛星影像及氣候數據，制作的國內外灌溉產品多為大尺度。然而有研究表明，以地塊為單元的農業遙感分析更符合農業生產的應用，而傳統基于像元的遙感分析方法往往會受到椒鹽噪聲的干擾，使地塊的基本形態遭到破壞，因而分類精度有限，所以本文以地塊尺度展開研究是核心要點之一。此外，在耕地較為精細、破碎的農田區域，單時相衛星影像所做出的地圖往往達不到所需精度[13]，而時序遙感提供了大量的高質量、連續的地表信息，也方便進行多源數據融合和分析，所以一直是遙感研究的熱門領域[14]，故本文使用了2020年3-10月的Sentinel-2高分辨率時間序列影像，這在獲取連續的高質量研究區相關地表信息方面具有一定優勢。同時本文還運用了機器學習中的XGBoost模型進行灌溉耕地的訓練和預測，通過使用集成了機器學習、時序遙感、以地塊為尺度這3個特點的創新思路來探究該方法應用于灌溉耕地提取的潛力，并進行對比驗證以證明該方法的必要性。經驗證，本文的方法可實際應用于灌溉耕地的調查與監測，同時為該方向后續的研究提供了參考。

1 研究區和試驗數據

1.1 研究區概況

研究區是一個以南普拉特河為主要水源的農業灌溉集中區，有良好的灌溉數據基礎，它位于美國科羅拉多州東北部，這里氣候溫和干燥，年平均溫度在10～20 ℃左右，以丘陵為主，主要農作物有玉米、小麥等。這里的農田大部分依靠灌溉，灌溉系統主要通過地下水和河流水源來進行，常用的灌溉方法包括噴灌、滴灌、漫灌等，灌溉面積僅次于加利福尼亞、德克薩斯和愛達荷，居全國第四位（見圖1）。

圖1 研究區示意圖Fig.1 Schematic diagram of the study area

1.2 數據及預處理

1.2.1 Google地圖影像

研究選用Google 地圖影像來提取研究區的農田地塊。Google 地圖影像具有3個可見光波段、且其空間分辨率為1 m，能夠精細地描述研究區的農業種植場景。

1.2.2 Sentinel-2影像

本文的主要數據源是高分辨率多光譜成像數據哨兵2 號（Sentinel-2）衛星，地面分辨率分別為10 m、20 m 和60 m。本研究下載了2020 年3 月4 日至10 月25 日共48 幅L1C 級的遙感影像，使用歐空局（ESA）發布的插件Sen2cor 對這些影像進行大氣校正，生成L2A 級產品數據，并在ENVI 中對影像進行云掩膜。

1.2.3 灌溉樣本數據 CIL

研究選用2020 年的科羅拉多灌溉耕地數據集(Colorado Irrigated Lands, CIL)來訓練和驗證所提出的識別模型。灌溉耕地數據集以矢量多邊形的形式提供了灌溉耕地的空間分布，并記錄了作物類型、灌溉方式、灌溉面積等信息（見表1）。

灌溉耕地數據集為模型訓練和驗證提供了足夠的正樣本。首先將灌溉耕地數據集CIL空間連接到研究區耕地地塊專題圖FP 上，然后計算FP 上落入每個耕地地塊P 內的灌溉面積A；當A 大于地塊P 總面積的70%時，我們認為地塊P 為灌溉正樣本。接著，在耕地地塊專題圖FP 上，從非灌溉樣本地塊中隨機選取與灌溉正樣本相等數量的地塊，作為負樣本集。最后，研究共獲得2 000 個地塊灌溉樣本，并按照2∶2∶6 的比例，將其劃分為訓練集、驗證集和測試集。

2 研究方法

本文研究流程如圖2 所示，主要包括以下內容：利用Google地圖影像進行耕地地塊的提取與修正，并結合經過預處理的Sentinel-2 時序數據集進行空間映射及地塊特征處理，得到地塊時序特征數據集后，結合灌溉樣本，利用XGBoost模型訓練預測灌溉耕地，最后對分類展開分析與評價。

圖2 研究流程圖Fig.2 Research flow chart

2.1 耕地地塊提取

地塊提取參考了Wang 等人[20]提出的邊界語義融合深度卷積網絡，該方法的核心是將具有準確空間定位的淺層邊界特征與用于類別識別的深層語義特征相結合進行地塊識別，從而預測出整個研究區的耕地地塊并進行修正。圖3展示了模型預測出的耕地地塊結果，結果表明，即使在農田分布密集的地區，模型也可以精確提取出耕地的精細邊界，并且對河流、建筑等干擾地物有著良好的辨別能力[圖3（a）、圖3（b）]，這為本研究灌溉耕地提取提供了空間約束。

圖3 耕地預測結果Fig.3 Farmland prediction results

2.2 特征構建與組合

對48景Sentinel-2數據基于地塊單元進行分區統計，分別獲取每一期影像中b2～b8、b8a、b11～b12 共10 個波段的均值、中值、標準差及眾數作為特征變量，用于遙感分類的候選變量。為了比較采用不同光譜特征下的灌溉耕地識別精度，設計了11 種不同的組合方式如表2 所示，意在對比均值、中值、標準差、眾數這4類光譜特征及不同種類之間的相互組合在灌溉耕地識別上的不同效果，優選出其中精度最高的特征組合方式對研究區進行灌溉耕地識別。

表2 不同特征類型組合方式Tab.2 Combination of different feature types

2.3 灌溉耕地提取模型

本研究使用的是XGBoost 算法，由Chen 等人[21]在2016 年提出，這是一種基于集成學習原理Boosting 的機器學習方法，除了有精度高、速度快的優點，XGBoost 算法還可以對影像缺失值進行處理[22]。它以CART 決策樹作為基分類器，充分利用了多核CPU 并行計算的優勢，大幅度提高了模型的運算速度和預測分類精度[23]，基模型為：

式中：xi是第i個樣本的預測值；fk(xi)是第k棵樹對數據集中第i個樣本的計算分數；F是所有樹的集合。

并定義XGBoost目標函數為：

式中：N為樣本的數量；為損失函數；為正則化項。

損失函數衡量模型與數據的吻合程度，正則化項衡量模型的復雜性。對目標函數的泰勒展開式進行整合、重組，轉化為與預測殘差相關的多項式，得到葉節點最優權重和目標值最優解分別為：

2.4 精度評價指標

在機器學習的分類任務中，常用的精度評價指標有很多，本研究中使用的指標共5 種（以下公式中，TP 指真正例（True Positive），TN 指真反例（True Negative），FP 指假正例（False Positive），FN 指假反例（False Negative）），相關指標定義詳見文獻[24]。

（1）總體精度（Overall Accuracy - OA）?？傮w精度衡量的是分類正確的比例，計算公式如下：

（2）F1分數（F1_Score）。F1是精確率和召回率的調和均值，同時考慮了查準率和查全率，F1值越大認為學習器的性能越好，計算公式如下：

（3）精確率（Precision）。精確率又叫查準率，是指被預測為正樣本的所有樣本中預測正確的占比，代表對正樣本結果的預測準確程度，計算公式如下：

（4）召回率（Recall）。召回率又叫查全率，是指在實際為正的樣本中被預測為正樣本的概率，計算公式如下：

（5）Kappa 系數（Kappa Coefficient）。Kappa 系數綜合考慮了分類器預測結果與實際情況之間的差異和隨機誤差，可以衡量分類器在各類別上的分類能力是否超過了隨機選擇的水平，計算公式如下：

式中：po也就是總體分類精度；pe是所有類別分別對應的實際與預測數量的乘積的總和除以樣本總數的平方。

3 結果與分析

3.1 精度分析

3.1.1 特征組合精度分析

本研究中，在地塊尺度下利用11 種特征組合進行灌溉耕地的識別，得到不同組合方式下的識別精度（表3）。結果表明，包含了均值、中值、標準差這三類光譜特征的組合8獲得了最高精度0.850 3，對灌溉耕地的識別最有利。由組合1、2、3可知單個變量的加入對提升分類精度影響不大，從組合11的精度表現可以看出眾數的參與會導致分類精度下降，在灌溉耕地的分類中應避免輸入該特征以降低被誤分類的可能?；诖?，在提取灌溉耕地時可將均值、中值、標準差輸入模型，以提高識別精度，準確地提取特征并分類，同時也保證了模型的泛化能力和穩定性。

表3 各組合不同指標下的分類精度Tab.3 Classification accuracy under different indicators of each combination

3.1.2 時相敏感性分析

從時間維度來探究灌溉耕地識別精度的變化規律，將不同長度時間序列遙感影像數據以月為單位輸入模型，數據長度等差增長，將完整的時間序列進行分割，得到了不同時間節點的灌溉耕地識別情況。

圖4 中精度總體呈逐漸上升趨勢，當加入4 月和5 月的特征時，精度增長幅度最大，此時正是農作物迅速生長的季節，后續繼續加入特征時增長幅度逐漸趨于平緩，是由于農作物于9、10月份完全成熟進入了收獲期。再結合研究區當地自然條件，5 月前和9 月后的天氣多為寒潮及大范圍降雪[24]，積雪的覆蓋不利于影像的識別，這也是導致精度不高的原因之一。綜上得出下列結論：對灌溉耕地提取較為敏感的時相集中在農作物生長季的中后期，該時相的遙感影像用于提取灌溉耕地可獲取較高的制圖精度。

圖4 不同時序長度下總體精度變化圖Fig.4 Overall accuracy variation chart under different time series lengths

總體來看，灌溉耕地識別的最佳時間與作物的生長階段特征關系明顯，在成熟期作物出現較為獨特、明顯的特征后識別精度將會大幅度提升，而當作物早期特征不明顯時，應使用完整的生長序列提取灌溉耕地以保證識別精度。從3月份到5月份精度迅速增加，分類精度提升效果明顯，說明這一時期的數據含有較多信息量。繼續增加時間長度，可以看到精度仍不斷上升，證明時間序列長度的增加可以有效提升分類精度，作物在8月份之后可以達到84%以上的精度，是一個比較理想的水平。當使用完整時間序列長度的數據時，識別精度最高，可見隨著時間序列長度的增加，特征數量也隨之增加，識別精度能夠達到最高水平。

以上結果表明，與傳統的單時相影像相比，結構化的時間序列影像數據不僅可以降低天氣對識別的影響，同時還蘊含了更多的特征信息，有利于灌溉耕地的提取。

3.2 灌溉耕地提取結果

根據3.1.1 節的結論，將精度最高的特征組合8 放入模型生成研究區灌溉耕地空間分布圖，并給出該模型的混淆矩陣，對精度進行評價，同時還使用特征重要性這一評價標準對每個特征要素在模型預測中的作用進行了評估。

3.2.1 灌溉耕地空間分布

由圖5可知，南北兩端灌溉耕地的規模很小，多數分布于研究區的中西部地區及東北部。灌溉耕地靠近南普拉特河以及一些水庫，是由于美國西半部干旱缺水，所以在許多河干、支流上興建了大型的引水工程用于農業灌溉，靠近河流便于獲取水源進行輸水，依靠灌溉手段才能更好地發展農業。例如，在該州的艾爾伯特縣和沙拉摩亞縣，大量種植玉米、小麥和大豆等灌溉作物。模型提取出的灌溉耕地大量聚集在研究區的最左側邊緣，因為這里靠近科羅拉多河流域，在這片流域內有許多重要的灌溉水源，如科羅拉多河的支流莫阿布河、弗拉格斯塔夫河、漢尼拔河等，這些河流的水源較為充足，可以為灌溉作物提供足夠的水資源。

圖5 研究區灌溉耕地預測結果Fig.5 Prediction results of irrigated farmland in the research area

3.2.2 精度評價

此次分類模型各項表現如表4所示，Kappa系數達到0.69，顯示出了相對較高的一致性水平，總體精度達到85.03%，意味著模型能夠正確地分類大部分樣本，其中灌溉耕地分類精度為86.76%，非灌溉耕地分類精度為82.30%。綜上，該分類模型在整體上表現良好，并且在灌溉耕地和非灌溉耕地的分類任務上都有著有優異的表現，各項精度指標均優于Zhu 等人[10]的研究。

表4 灌溉耕地分類混淆矩陣Tab.4 Classification confusion matrix of irrigated farmland

3.2.3 特征及重要性評價

（1）最佳特征變量。優選出20 個表現最好的特征變量如圖6 所示，其中包含了10 個均值特征、7 個標準差特征、3 個中值特征，0個眾數特征，可以看出均值和標準差對提升整體分類精度的影響較大。在眾多特征變量中，5 月13 日影像b8波段的標準差、8 月21 日影像b1 和b3 波段的均值的特征重要性遠高于其他特征變量，在分類中起到重要作用。

圖6 組合8前20個最佳特征變量Fig.6 Combining the top 20 best feature variables of 8

（2）時序特征。即使是相同類型的特征變量在不同時相也會具有不同的敏感性，圖7中每一期影像的特征重要性都由均值、中值、標準差這三類變量共同組成，但其在灌溉耕地識別中的作用截然不同。在美國，玉米、大豆等主要農作物的物候期一般集中在4-11 月，從圖7 可以明顯看出有利于識別的遙感影像集中在5-9月，而該時間段正是農作物的生長期及成熟期，此結論與3.1.2 節相符。也有部分影像的重要性為0，可能是這幾期影像云量較大范圍較廣，覆蓋了整個研究區，使用這些重要性較低的特征變量進行制圖對分類精度的提升可能沒有作用甚至起到副作用，所以在進行灌溉制圖前有必要進行特征的選擇工作。

圖7 時序特征重要性Fig.7 Importance of temporal features

（3）不同類型特征的重要性。從圖8可以看出，均值變量在分類中的重要性最高，對灌溉耕地的識別最有幫助；標準差的重要性略低于均值，但總體重要性得分也高；比起均值和標準差，中值這一類型的特征變量在參與分類任務時的作用較小，但對灌溉耕地和非灌溉耕地的分類也具有一定影響。

圖8 特征類別重要性Fig.8 Importance of feature categories

4 討 論

4.1 分類方法對比

本文利用高分辨率Google 影像，結合研究區耕地的形態特征，利用邊界語義融合深度卷積網絡提取當地的耕地邊界，該方法與常規模型進行地塊分割或提取的方法相比，對本文農業灌溉集中區的耕地提取更加精細，鮮少有錯提、漏提的現象，為后續耕地分類提供了良好的基礎。在使用時序遙感影像方面，與單時相遙感影像相比，本文充分考慮科羅拉多州地區作物種植類型及耕地類型復雜多樣的特點，通過綜合多類型、多時相的特征變量提升分類精度。在分類模型選擇方面，和常用的隨機森林方法相比，本文選擇的XGBoost分類方法已被證實有速度快、精度高等優點，且對于本實驗采用的光學數據源由于云的覆蓋導致的特征信息缺失的部分，XGBoost分類方法能對缺失值進行相關處理[22]，因此該方法更適用于本研究。

4.2 分類結果不確定性分析

由于研究區內作物種植類型及耕地類型的多樣性，使分類工作成為一個很大的挑戰，分類過程中引入的一些誤差會導致分類結果的不確定性，具體內容如下：

（1）地塊提取的工作仍有一定的進步空間，由于影像的質量的問題，還是存在一定的地塊錯分、漏提的現象，這也進一步導致了在耕地識別時發生錯誤的可能性。

（2）不同類型耕地上種植的作物可能由于生長期相近，導致有相似的特征信息，這不利于區分耕地類型。

（3）光學遙感影像易受天氣影響，當云霧或水汽等遮擋時會導致圖像質量下降或無法成像，導致在特征構建時有缺失現象，會在一定程度上影響模型分類精度。

4.3 問題與展望

由于Sentinel-2 易受天氣影響的原因，在后續的工作中，考慮增加Sentinel-1 等微波數據以多源數據融合的方式進一步展開研究，利用微波遙感滲透力強、全天候工作的優勢，以提高分類精度。此外還可構建物候、紋理、地形等特征解決“同物異譜”或“異物同譜”的現象，通過更多的特征組合方式建立更加精確的灌溉耕地提取模型。在算法方面，數據降維、特征選擇等因素會對分類算法產生很大的影響，可以結合實際對比擇優使用更加穩定、精確的模型，選擇適當的分類算法以達到更高的分類精度。

5 結 論

本文以美國科羅拉多州南普拉特河流域的農業灌溉集中區為研究區，綜合利用高分辨率遙感影像和Sentinel-2 時間序列影像開展地塊尺度的灌溉耕地提取，分析不同時相和不同特征組合情況下的識別精度，探究了XGBoost模型在時序遙感灌溉耕地提取中的應用潛力，得出以下結論：

（1）相比于傳統基于像元的遙感灌溉提取，本文基于地塊的提取方法不會破環地塊基本形態，在耕地情況破碎、復雜的區域也能夠實現灌溉耕地的精確提取，且長時間序列的影像蘊含更豐富的特征信息，為灌溉耕地的提取提供了有力的支撐，結合XGBoost機器學習模型，此次分類總體精度高達85.03%；

（2）在特征的構建與組合中，不同類型的光譜特征表現不同，精度表現最好的是均值、標準差、中值這三類光譜特征的組合，Kappa 系數達到0.69，其中灌溉耕地分類精度為86.76%，非灌溉耕地分類精度為82.30%；

（3）從時相敏感性分析結果可知，對灌溉耕地和非灌溉耕地的區分較為敏感的時相集中在農作物生長季的中后期，其分類精度隨著時序長度的增加而不斷提高。

本研究利用高分辨率影像的光譜特征進行灌溉耕地的識別，為灌溉制圖提供了研究思路，對準確掌握灌溉耕地的數量和空間分布等信息有重要的現實意義，將來有望廣泛應用于大區域乃至全球灌溉耕地空間產品研制。