基于Sentinel-2A的棉花種植面積提取及產量預測

王匯涵，張 澤，康孝巖，林 皎，印彩霞，馬露露， 黃長平，呂 新※

（1.石河子大學農學院/新疆生產建設兵團綠洲生態農業重點實驗室，石河子 832000； 2.中國科學院空天信息創新研究院，北京 100094；3.塔里木大學農學院，阿拉爾 843300）

0 引 言

棉花是重要的經濟作物。新疆作為中國最大的棉花種植基地，2021年棉花產量512.9萬t，占全國棉花產量的89.5%，在棉花產業中發揮著舉足輕重的作用。棉花產量是指導田間管理的指標之一。及時、準確地預測棉花產量在棉田經營管理、農業決策制定等方面具有重要的價值和意義。

遙感技術具有快速、宏觀、動態、成本低等優點，能夠克服傳統作物估產方法的局限性，已成為產量估測的新技術。目前，已有學者對植被指數和作物產量的模型進行了大量的研究，基于時間序列MODIS_NDVI的農作物產量預測研究較為廣泛。趙文亮等根據 NOAA/AVHRR NDVI時間序列得出作物生長期內的歸一化植被指數（Normalized Difference Vegetation Index, NDVI）累計值或最大值，通過線性擬合，建立作物產量遙感預測模型。Ji等利用MOD09Q1數據產品計算NDVI時間序列，用于作物大規模產量預測。由于地物類型復雜多樣，MODIS數據受空間分辨率的限制，對預測效果具有負面影響，低空間分辨率仍然是遙感技術預測產量應用中的一個問題。劉煥軍等利用30 m空間分辨率的Landsat_5_TM與Landsat_7_ETM影像構建非多天合成的NDVI時間序列數據建立棉花估產模型。魯新新等以Landsat 8遙感影像為數據源，確定NDVI與棉花產量相關性最優的算法模型。這些研究證明了傳統植被指數與產量之間的關系，避免了MODIS數據低空間分辨率、多天合成產生混合像元的問題。但也有研究表明，單一植被指數NDVI在較大的葉面積指數（Leaf Area Index，LAI）下易出現飽和現象，并受土壤背景影響，與棉花產量相關性較低（=0.47），Zhao等認為，皮棉產量與棉花生長早期測得的冠層反射率指數相關性最好（為0.56～0.89）。朱婉雪等基于冬小麥不同生育期影像數據和植被指數，構建了多個估產模型，最終確定最佳估產時期和植被指數為抽穗灌漿期和增強型植被指數（Enhanced Vegetation Index，EVI），決定系數達到0.70。

綜上所述，以往研究主要通過建立單一植被指數與棉花產量的線性或非線性回歸方程，比較分析其相關性，選用相關性最高的某一植被指數與棉花產量的回歸方程作為棉花產量預測模型，而棉花產量形成是一個連續變化的動態過程，受冠層結構、生物量和葉綠素含量等因素的影響。單一植被指數的應用簡化了棉花產量預測的復雜性，增加了模型的誤差和不穩定性。針對該問題，本研究以新疆莫索灣墾區為研究區，基于Google Earth Engine（GEE）平臺，獲取棉花3個生育時期的Sentinel-2A影像，并進行預處理，重采樣至10 m空間分辨率，計算14種植被指數，探討不同指數預測棉花產量的能力，確定棉花產量預測的最佳生育時期，建立最優的棉花產量預測模型，實現莫索灣墾區棉花產量預測制圖，以期為棉花精準管理提供決策信息。

1 研究區與數據來源

1.1 研究區概況

研究區位于新疆兵團第八師莫索灣墾區，地處天山北麓中段，準噶爾盆地南部，東經84°58′～86°24′，北緯43°26′～45°20′，平均海拔300～500 m，屬典型的溫帶大陸性氣候，冬季長而嚴寒，夏季短而炎熱，年平均降水量180～270 mm，年平均氣溫7.5～8.2 ℃，光熱資源豐富，日照2 318～2 732 h，無霜期147～191 d。研究區域內種植少量小麥、苜蓿等作物，以棉花為主栽作物。

1.2 數據來源

以2020年Sentinel-2A影像數據為主要數據源，包含13個光譜波段，具體波段參數如表1所示。在GEE平臺中，利用JavaScript編程計算云信息波段QA60的Bit10和Bit11的值，設置二者的值均為0得到云掩膜，使用云掩膜去除影像中的云信息，得到無云影像，并重采樣至10 m分辨率。

表1 Sentinel-2A波段信息 Table 1 Band information of Sentinel 2A

地物分類數據：數據來自Google Earth 高分辨率影像和野外調查。為了充分掌握研究區域內的地物類型及其典型樣本的空間位置，調查時利用Google Earth 高分辨率影像記錄代表性地塊的經緯度坐標及相應的植被類型。樣本數據分布均勻，覆蓋研究區域范圍，用于農田提取的樣本共計609，包括非農田樣本282，農田樣本327；用于棉花提取的樣本共計207，其中非棉花樣本74，棉花樣本133。

棉花實測產量數據：在棉花收獲期開展測產工作，在研究區內，每個團場選擇4～5個采樣站點，每個站點附近范圍內選取5個樣點，共計95個樣點進行模型構建，同時選擇21個樣點進行模型驗證（圖1）。按照五點采樣法，調查株數、鈴數，取一膜一米的棉絮稱干質量，計算平均鈴質量。

圖1 研究區域及樣點分布圖 Fig.1 Study area and distributions of sample points

植被指數是多個光譜帶和波長的轉換，可以解釋冠層結構和光合作用的變化。根據光譜特性，本文研究了產量預測中常用的14種光譜植被指數，包括6個冠層結構指數及8個葉綠素相關指數（表2）。

表2 文中選用的植被指數 Table 2 The selected vegetation index

2 數據處理及分析方法

利用GEE平臺對Sentinel-2A影像進行預處理。采用兩步分級制圖策略，利用分類算法剔除非農作物分布區域，合成棉花不同生育時期的影像，實現棉田與非棉田的劃分。利用順序向前算法（Sequential Forward Selection，SFS）對波段和植被指數進行特征篩選，基于偏最小二乘回歸模型（Partial Least Squares Regression，PLSR）構建不同生育時期的棉花產量預測模型，以確定棉花估產最佳時期，并將模型應用到棉花產量制圖中。

2.1 兩步分級棉花種植區域提取策略

本研究中實施兩步分級策略用于棉花種植區域提取。第一步，將研究區域劃分為5個主要的土地覆蓋類別（農田、荒漠、建筑、城市綠地和水體），基于10 m空間分辨率Sentinel-2A影像，選擇17個特征用于農田提取，包括12個光譜特征、3個植被指數特征、2個地形特征（Digital Elevation Model ，DEM）特征剔除非農作物分布區域，得到準確的農田掩膜。第二步，將農田劃分為棉花區和非棉花區。對研究區影像進行掩膜處理得到農田植被的總體分布，確定棉花5個生育時期，將各生育期的5幅合成圖像組合，形成包括12個光譜特征、3個植被指數特征、2個DEM特征的新圖像，進行棉花種植區提取。

2.2 分類方法

利用隨機森林（Radom Forest, RF）、支持向量機（Support Vector Machine, SVM）算法和決策樹（Classification and Regression Tree, CART）三種分類方法，結合研究區域及算法本身的特點，在樣本數據中隨機選取 70%作為訓練數據，剩下的30%作為驗證數據。分類性能的評估基于混淆矩陣（Confusion Matrix），用行列的矩陣形式表示。每一行代表參考數據中的實際類，每一列代表預測類，對角線元素表示正確分類的像素個數?；煜仃囍休^常使用的具體評價指標有總體精度（Overall Accuracy, OA）、制圖精度（Producer’s Accuracy, PA）、用戶精度（User’s Accuracy, UA）以及Kappa系數（Kappa Coefficient, KC）等，這些指標從不同的角度反映作物分類的精度。

2.3 特征選擇

由于棉花早期生長階段的田間空間變異性小，考慮土壤背景反射率的潛在負面影響，開花后期的棉花產量和植被指數之間的相關性更好，因此本研究選取棉花3個主要生育時期的遙感影像。7月11日－7月20日，7月21日－8月30日，8月31日－9月10日，分別對應棉花花期、鈴期、吐絮期。利用SFS方法對波段和植被指數進行篩選。

2.4 棉花產量預測模型

PLSR是一種融合了主成分分析、回歸分析和典型相關分析等方法的多元數據統計分析模型，能夠同時考慮各特征之間的相關性以及與因變量之間的相關性，實現數據的降維、信息綜合與篩選，較好地解決特征共線性問題，有效剔除不重要的變量。本文利用SPXY（sample set partitioning based on joint x-y distance）算法同時計算樣本間距離對樣品集進行劃分，基于PLSR構建不同生育時期棉花產量預測模型，并采用獨立驗證數據對預測模型進行評估和驗證，最終確定棉花最佳估產模型。評價指標包括：決定系數（），均方根誤差（RMSE）和相對誤差RE（%），其值越大和 RMSE 值越小，說明模型的預測性能和準確性越高。均方根誤差（RMSE）和相對誤差RE（%）公式如下：

式中P、O和分別為預測值、觀測值和實測值的平均值，為樣本數量。當RE≤10%時，模型較優，10%＜RE≤20%模型適中，RE＞20%模型較差。

3 結果與分析

3.1 棉花種植區域提取

結合多種分類特征對莫索灣墾區棉花種植區域進行提取。利用驗證點數據對農田與非農田分類進行交叉驗證，結果如圖2a所示，其中，RF、SVM、CART的UA分別為0.96、0.97、0.91；PA分別為0.92、0.83、0.74；OA精度分別為0.94、0.90、0.83；KC分別為0.89、0.81、0.67。除UA相差較小以外，其他指標差距較大，RF算法分類效果明顯優于SVM和CART，更適用于農田與非農田的提取。

在農田掩膜基礎上，根據棉花不同生育時期影像，進行棉田與非棉田分類，分類結果如圖3b所示。計算分類結果的混淆矩陣，對不同分類方法的分類結果進行精度驗證，驗證結果如圖2b所示，RF算法的PA為0.95，SVM算法PA為0.93，CART的PA為0.90。RF與CART分類精度相差不大，SVM的精度最低，其OA為0.89，KC為0.78。總體而言，RF算法最佳，對2020年莫索灣墾區的棉花種植面積進行計算，并與實際種植面積進行比較，分類結果面積為60 400 hm，實際面積為64 866.7 hm，相對誤差為6.9%，與團場棉花的種植現狀相吻合。

圖2 基于不同算法的分類精度評價 Fig.2 Classification accuracy evaluation based on different algorithms

圖3 基于隨機森林算法的兩步分級棉花種植面積提取 Fig.3 Two-step grading cotton planting area extraction based on random forest algorithm

3.2 不同生育時期特征變量選擇

為進一步分析不同生育時期波段/植被指數對棉花估產的影響，對研究區棉花產量與不同生育時期波段/植被指數之間的相關性進行統計（圖4）?；ㄆ?、鈴期、吐絮期波段/植被指數與產量的相關性趨勢基本相同，其中B6波段與產量相關性最高，隨著生育時期的推進，相關系數逐漸增大，依次為0.37、0.47、0.53。其次是B7波段，相關系數依次為0.32、0.45、0.51。近紅波段僅次于紅邊波段，B8A波段與產量的相關系數分別為0.31、0.42、0.50?？梢?，Sentinel-2A特有的紅邊波段以及近紅外波段對棉花產量預測具有一定幫助。

圖4 不同特征變量（波段/植被指數）與產量的相關性分析 Fig.4 Correlation analysis between different characteristic variables (Band/Vegetation Index) and yield

相關分析結果能夠反映出單個波段/植被指數的預測性能，但并不能有效推斷出多波段/植被指數組合的預測能力。利用SFS方法對波段/植被指數進行特征選擇，結果如表3所示，就篩選順序而言，在三個生育時期中篩選的第1個特征都為紅邊波段（B6波段），進一步證明紅邊波段在棉花產量預測中的重要性。就數量而言，波段在花期對棉花產量較為敏感，篩選數量多于冠層結構指數和葉綠素相關指數，依次為6個、1個、4個。鈴期葉綠素相關植被指數篩選數量（6個）明顯多于波段（3個）和冠層結構指數（2個）。吐絮期冠層結構相關指數5個，其次波段4個，葉綠素相關指數2個。冠層結構相關指數隨著生育期的推進，對產量的敏感性增大，在棉花生育后期對產量預測模型的貢獻較多，葉綠素相關指數則呈現先增長后減小的趨勢，在鈴期貢獻最為突出。弱相關性To指數均被篩選為特征變量，相關系數依次為0.02、0.03、0.10，可見，弱相關性特征的預測性能未必低于獨立性弱的強相關性特征組合。對產量預測篩選特征波段時不僅僅考慮單一特征的性能，還要衡量特征之間的獨立性和冗余度。

表3 特征變量篩選結果 Table 3 Characteristic band screening results

3.3 不同生育時期棉花產量預測模型

本研究利用SPXY算法對樣品集進行劃分，以產量為變量，波段/植被指數為變量，將70%的樣本訓練，剩余30%為預測集樣本?；赑LSR將篩選出的特征變量與產量構建模型，結果如圖5。鈴期模型預測效果最佳（=0.62，RMSE=625.5 kg/hm，RE=8.87%），優于吐絮期（=0.51，RMSE=789.45 kg/hm，RE=11.06%）、花期（=0.48，RMSE=686.4 kg/hm，RE=9.86%）。預測模型的精度隨著生育時期的推進有了較大的提升。

圖5 基于偏最小二乘算法不同生育時期產量預測模型評價 Fig.5 Evaluation of yield prediction model for different growth periods based on PLSR algorithm

3.4 棉花產量預測模型測試及制圖

利用2020年棉花鈴期的遙感影像數據（10 m空間分辨率），基于構建的鈴期產量預測模型，選用獨立的21個樣本點測試模型性能。測試結果如圖6所示，=0.61，RMSE=425.25 kg/hm，RE=6.61%，與所建立模型的驗證結果接近（=0.62，RMSE=625.5 kg/hm，RE=8.87%），比較理想。

圖6 鈴期模型預測產量與實測產量的比較 Fig.6 Comparison between predicted yield and measured yield of Boll period model

此外，將已篩選的波段/植被指數作為輸入因子，逐像元計算棉花產量，對莫索灣墾區棉花產量分布制圖。從圖7可以看出，棉花產量分布在3 000～9 000 kg/hm范圍內，與實際情況相符，總體而言，棉花產量由北向南逐漸減少，高值區出現在150團，低值區則在147團、148團均有分布。

圖7 新疆莫索灣墾區棉花產量預測 Fig.7 Prediction of cotton yield in Xinjinag Mosuwan Reclamation Region

4 討 論

4.1 不同算法提取棉花種植區域精度驗證結果

本研究采取兩步分級提取策略實現了新疆莫索灣墾區2020年棉花種植面積提取。與李方杰等的研究方法相比，本文使用10 m空間分辨率Sentinel-2A影像數據，在影像合成上突出了棉花的特征，且在模型訓練時加入植被指數和地形特征，實現較高精度的棉花提取，該提取策略聚焦田間植被，篩選得到的特征針對性更強，理論上所得結果的適用性和推廣性更優且不受其他非農地物的影響，在棉花衛星遙感提取中具有重要意義。以往研究表明，機器學習算法與傳統算法相比在處理高維復雜數據時可以有效利用數據特征實現更高的分類精度。然而，使用多種機器學習算法進行對象提取的研究較少，在相同數據上比較多種算法有利于在棉花種植面積提取中找到最合適的算法。鑒于此，本研究憑借GEE存儲大量遙感數據的優勢，獲取了2020年覆蓋莫索灣墾區的Sentinel-2A數據影像，同時實現了去云、影像鑲嵌、裁剪等一系列預處理，有效避免了影像鑲嵌帶來的色差問題，保證了分類數據源的純潔度與質量，利用GEE中提供的分類器，對比RF、SVM、CART三種算法在棉花種植區域提取中的分類效果。從PA、UA、OA和KC檢驗結果來看，RF對農田與非農田分類的PA為0.92，UA為0.96，OA為0.94，KC為0.89，優于SVM和CART。棉田與非棉田分類中，RF與CART表現相似，RF的PA（0.95）略高于CART（0.90），由此可見RF的分類效果最佳，這與哈斯圖亞等的研究結果相符。CART容易出現過度擬合現象，而RF是一個決策樹的集合，只考慮特征的一個子集，最終結果取決于所有的樹，有更好的泛化能力和更少的過度擬合。該分類方法對農作物的識別能力較強，分類精度較高，能夠滿足精準農業遙感監測的需要。

4.2 波段/植被指數與產量的關系及不同生育時期預測效果

本研究利用GEE平臺獲取并預處理棉花三個生育時期的Sentinel-2A遙感影像，構建14種植被指數，與產量進行相關性分析。與花期相比，棉花產量和生育后期（鈴期/吐絮期）的波段/植被指數之間相關性更高。可能由于現花期田間空間變異較小，土壤背景反射率存在負面影響，這一結果與前人研究一致。由于棉花吐絮后冠層結構變化較大，吐絮期的冠層結構指數與產量的相關性高于其他的生育時期，并且有明顯的差異性，而葉綠素相關植被指數則在棉花盛鈴期表現較好，冠層結構和葉綠素含量相關植被指數對作物生長過程中的不同方面（即形態和生理）都有反應，因此兩者都包括在最終模型中。波段/植被指數與棉花產量的相關分析結果能夠反映出單個波段/植被指數的反演性能，但并不能有效推斷出多波段/植被指數組合的反演能力。由于波段/植被指數之間有著不同程度的信息重疊、冗余以及相關性，從而獨立性強的弱相關特征組合的反演性能未必低于獨立性弱的強相關的特征組合。

由于Sentinel-2A特有的紅邊波段與產量有較高的響應，將與冠層結構和葉綠素含量相關的植被指數結合起來建立模型，在不同生育時期的模型對比中，鈴期的模型預測效果最好(=0.62，RMSE=625.5 kg/hm，RE=8.87%)，優于吐絮期（=0.51，RMSE=789.45 kg/hm，RE=11.06 %），花期的預測結果較差（=0.48，RMSE=686.4 kg/hm，RE=9.86%）。前人研究也發現隨著季節的推進，產量預測效果改善。本研究中的預測模型在鈴期時有所改善，主要歸因于特征篩選時鈴期葉綠素相關指數的貢獻更大。棉花吐絮期模型精度低于鈴期的原因可能在于棉花在進入吐絮期后，白色棉絮漸漸顯露出來，冠層結構發生變化，吐絮前后冠層反射率產生較大差異，植被指數受棉絮影響，導致低估?，F蕾營養生長與生殖生長并進，營養物質及環境資源競爭激烈，其中鈴期環境條件的影響最明顯，該生育期為棉花產量形成的關鍵期。

4.3 莫索灣墾區棉花產量預測模型測試及制圖

基于上述構建的產量預測模型，選用單獨的21個樣本點對模型性能進行測試，其結果較為理想（=0.61，RMSE=425.25 kg/hm，RE=6.61%），并實現了新疆莫索灣墾區棉花產量預測制圖。產量預測值分布在3 000～9 000 kg/hm的范圍內，對于每個像元而言，存在高估或低估屬于正?，F象。前人的研究也表明，由于作物管理投入的差異，也會存在同一區域內產量變異性大，在合理誤差范圍之中?？傮w而言，分布呈現出高產集中在北部區域，由北向南遞減，這可能與150團緊鄰沙漠有關，沙漠增溫效應可以使準噶爾盆地100 km內的區域溫度有不同程度的提高，距沙漠近的區域有較好的熱量資源，促進生理活動，增加產量。

5 結 論

1）本研究采用兩步分級提取策略，有效避免了非農田的影響，準確提取莫索灣墾區棉花種植區域。其中隨機森林算法表現最佳，農田與非農田分類的制圖精度PA和用戶精度UA分別達到了0.92和0.96，分類總體精度OA為0.94，Kappa系數KC為0.89，棉田與非棉田分類的PA和UA分別達到了0.95和0.87，分類OA為0.92，KC為0.83。分類結果面積為60 400 hm，實際面積為64 866.7 hm，相對誤差為6.9%，表明該分類方法對棉田分類精度較高，能夠滿足精準農業遙感監測的需要。

2）本研究基于SFS-PLSR算法構建了3個生育時期的棉花產量預測模型，利用實測數據對模型進行驗證，結果表明，模型在鈴期表現出了最佳產量預測能力（=0.62，RMSE=625.5 kg/hm，RE=8.87%），鈴期是新疆莫索灣墾區棉花產量預測的最佳生育時期。根據棉花種植區域圖實現棉花產量預測制圖，莫索灣墾區的棉花產量范圍為3 000～9 000 kg/hm。

本研究提出的兩步分級策略能夠有效識別棉花種植區域，綜合考慮波段/植被指數組合在產量預測中的表現，為新疆莫索灣墾區棉花產量預測提供了一種新的思路，可為合理水肥配置、精準種植、農作物生長過程監測提供數據支撐，促進精準農業快速發展。