基于Sentinel-2多光譜數據的枸杞葉面積指數反演

余剛馮克鵬

(1.寧夏大學土木與水利工程學院，寧夏 銀川 750001；2.寧夏節水灌溉與水資源調控工程技術研究中心,寧夏 銀川 750021；3.旱區現代農業水資源高效利用教育部工程研究中心，寧夏 銀川 750021)

引言

葉面積指數(Leaf area index，LAI)是模擬陸地生態過程、生物地球化學循環以及水熱循環的重要參數[1]，在農業、林業、生態、氣候變化等領域已經得到了廣泛的應用[2]。LAI能夠反映植被長勢個體和群體特征，是農作物長勢監測的一個關鍵生態參數[3]，因此快速、準確、可靠地獲取農作物LAI具有重要意義。

隨著遙感技術的飛速發展，過去的幾十年時間里，許多新發射的中分辨率對地觀測衛星投入了研究應用，如1999年12月18日、2002年5月4日發射的MODIS衛星，2013年2月11日發射的Landsat-8衛星，2013年4月26日發射的高分1號衛星，2015年6月23日、2017年3月7日由歐空局發射的Sentinel-2A和Sentinel-2B多光譜遙感衛星。Sentinel-2多光譜遙感衛星搭載包括空間分辨率為10m的3個可見光波段以及1個近紅外波段，空間分辨率為20m的特有3個紅邊參數波段和短波紅外波段，空間分辨率為60m的近紅外波段等，共計13個波段。Sentinel-2A和Sentinel-2B雙軌衛星同時運行，時間分辨率縮短至5d，為遙感衛星數據在農業遙感中的應用提供了更多的可能性[4-7]。

Sentinel-2衛星所搭載的紅邊參數波段在農業遙感中具有很高的有效性[8]。Frampton等[9]基于模擬的Sentinel-2數據提出了用于估算作物葉綠素和葉面積指數LAI的IRECI(inverted red edge chlorophyll index)和 S2REP(Sentinel-2 red-edge position)紅邊參數植被指數，Korhonen等[10]開展的森林葉面積指數估算，以及蘇偉等[11]利用Sentinel-2數據進行的玉米葉面積指數估算，這些研究基于Sentinel-2遙感影像所構建的的各種估算模型，都證實了Sentinel-2多光譜數據用于LAI估算的可行性，但研究中所構建的LAI估算模型均是針對特定研究區域和研究對象，無法直接用于寧夏地區枸杞LAI估算。同時在這些研究中，所構建的LAI估算模型大多是物理模型和經驗模型，無法直接應用在新的數據源上。

機器學習算法通過相關數據擬合一個靈活的模型，學習輸入變量與輸出變量之間的關系，將模型參數調整為在一個獨立驗證集中預測誤差最小，使模型獲得最好的泛化能力，其可以同時輸入多個與LAI相關性高的植被指數，并具有很好的非線性擬合能力，可提高模型的穩定性及估測精度[12]，是農作物LAI估算的重要手段。近年來，已有眾多學者使用機器學習算法進行LAI遙感估算[13-15]。梁棟等[16]基于支持向量機回歸算法，將不同植被指數與波段反射率作為輸入參數，在冬小麥多個生育期反演LAI具有較好的適用性。

本研究通過大面積田間實驗實測LAI，構建基于Sentinel-2多光譜遙感影像提取的波段反射率和計算的植被指數，構建3種不同數據輸入模型，運用機器學習算法進行估算，對比分析反演精度，研究結果將進一步豐富枸杞LAI估算方法，提高估算精度。

1 材料與方法

1.1 研究區概況

研究區位于寧夏回族自治區，地理坐標為E104°17′～107°39′，N35°14′～39°23′。海拔1100～1200m，平均氣溫-9℃，年降水量在150～600mm，屬中溫帶干旱氣候帶，干旱少雨，光照條件充足，蒸發強烈。

圖1 研究區位置

1.2 采樣區設置

研究區內枸杞種植廣泛，本研究選取包括2年期、4年期和多年期枸杞植株種植園，共6個采樣區，分布于寧夏中北部。枸杞采取溝壟相間的種植方式，在每個采樣區內選取面積較大且相對均質的枸杞地塊，每個地塊隨機設置1個10m×10m樣方。在全部采樣區共設置了80個樣方，并采用GPS記錄樣方經緯度坐標。采區設置及枸杞植株物理量如表1所示。

表1 樣區設置及枸杞植株物理量

1.3 枸杞LAI測定

地面數據采集于2019年8月15—18日，使用LAI-2200C冠層分析儀測量樣方內每行枸杞LAI值，取樣方內所有行平均值為該樣方LAI測量結果。LAI-2200C植物冠層分析儀是利用“魚眼”光學傳感器測量植被冠層5個不同天頂角方向的散射天空輻射衰減，通過5個天頂角方位的透光度結果計算出葉面積指數、空隙比等樹冠結構參數。LAI測量時應盡量避免在直射光背景下測量，測量模式采用AAAABBBB。第1個A值采用白色散射帽，測量人員背對太陽，使陽光直射散射帽，保持探桿水平測量；第2個A值將探桿移動至測量人員頭部陰影下測量；第3個A值將白色散射帽移除，保持在測量者頭部陰影下測量；第4個A值使用45°散射帽，正常測量。B值在樣方內部采用兩壟之間對角線測量，每條對角線上采集4個B值，第1個B值取在壟上，第2個B值取在兩壟之間1/4處，第3個B值取在兩壟中間，第4個B值取離壟3/4處。實測LAI統計如表2所示。

表2 樣區實測LAI統計

1.4 Sentinel-2遙感影像

遙感影像使用Sentinel-2多光譜衛星影像，通過歐洲航空局的數據共享網站(https://scihub. copernicus.eu/dhus/#/home)下載，獲取與地面觀測時間同步的影像數據，空間分辨率10m，成像時間分別為2019年8月15日和18日。下載的遙感影像均為已經過輻射校正和幾何校正處理的Level-1C級大氣上層表觀反射率。遙感影像大氣校正采用SNAP-Sen2Cor軟件處理，經大氣校正后的遙感影像通過最近鄰插值法重采樣至10m后用于提取單波段反射率和計算植被指數。Sentinel-2多光譜波段信息見表3。

表3 Sentinel-2多光譜傳感器參數

1.5 植被指數

利用Sentinel-2多光譜衛星可見光及近紅外波段計算多種植被指數，主要包括土壤校正型植被指數、紅邊參數植被指數、葉綠素含量植被指數以及傳統的近紅外植被指數，共計21種植被指數[17-21]。各植被指數基于Sentinel-2數據各波段計算公式如表4所示。本研究所采用的植被指數使用歐洲航空局SNAP軟件計算得到。

表4 植被指數及基于Sentinel-2數據的計算方法

1.6 機器學習算法

多元線性回歸模型(Multivariable Linear Regression，MLR)通常用來描述變量y和x之間的隨機線性關系，用來研究一個應變量依賴多個自變量的變化關系，特點是形式簡單、易于建模[33]。

逐步多元線性回歸模型(Stepwise Multiple Linear Regression，SMLR)是將變量逐個引入模型，每引入一個解釋變量后都要進行F檢驗，并對已經選入的解釋變量逐個進行t檢驗。這是一個反復的過程，直到既沒有顯著的解釋變量選入回歸方程，也沒有不顯著的解釋變量從回歸方程中剔除為止，以保證最后所得到的解釋變量集是最優、最簡單的[34]。

支持向量機(Support Vector Machine，SVM)是一種二分類模型，其目的是尋找一個超平面對樣本進行分割，分割的原則是間隔最大化，最終轉換為一個凸二次規劃問題來求解[35]。

高斯過程回歸模型(Gaussian Process Regression，GPR)是使用高斯過程先驗對數據進行回歸分析的非參數模型[36]。GPR有著嚴格的統計學習理論基礎,對處理高維數、小樣木、非線性等復雜的問題具有很好的適應性，且泛化能力強[37,38]。

樹的集成模型(Ensembles of Trees)是使用基本模型創建組合模型的學習方法。主要包括梯度提升樹和隨機森林2種，這2個算法都使用決策樹作為基礎模型，兩者都是學習集成樹的算法，但是訓練過程不同[31]。

回歸樹模型(Regression Trees)由Breiman等人在1984年提出，是應用廣泛的決策樹學習方法。由特征選擇、樹的生成及剪枝組成，既可以用于分類也可以用于回歸。回歸樹算法流程類似于分類樹，分枝時窮舉每一個特征的每一個閾值，來尋找最優切分特征j和最優切分點s，衡量的方法是平方誤差最小化，分枝達到預設的終止條件(如葉子個數上限)就停止[32]。

1.7 評價指標

選用決定系數(R2)、均方根誤差(RMSE)進行分析檢驗，其中R2表示模擬值與實測值的擬合優度，其值越接近1，表示模型的參考價值越高；RMSE主要用于模型驗證，反映了模擬值與真實值的偏離程度，其值越小，模型精度越高。指標的計算公式：

2 結果與分析

2.1 遙感影像反射率特征

對Sentinel-2多光譜遙感影像各波段反射率在不同LAI值時的變化特征進行分析，LAI值選取本研究中實測LAI最大值(LAI=2.62)、最小值(LAI=0.22)和平均值(LAI=1.01)的樣方所對應的Sentienl-2多光譜反射率進行對比，如圖2所示。3種不同LAI值情況下，光譜曲線均在400～700nm的可見光區域冠層反射率均呈現先上升再下降趨勢，在藍光波長450nm和黃光波長650nm左右有1個吸收谷，在綠光波長550nm左右有1個反射峰，這是因為植物葉片中葉綠素對藍光和紅光輻射吸收較強，對綠光輻射反射較強。在700～850nm的近紅外區間內，反射率隨著波長增大而增大，此譜段急劇上升，具有陡而近于直線的形態，在紅邊波段850nm左右出現最大值，形成“紅邊”現象。主要原因是由葉片細胞結構造成的，在這一區域，葉片結構導致對光譜的反射和透射達到最大，吸收最小，所以出現反射峰。隨著LAI值增大，反射率數值上開始逐漸增大，光譜曲線的位置增高，這是由于隨著葉片大小、數量和層數的增加，光線在葉片的不同層面發生多次反射，光線除了在上層葉片形成反射以外，還可以透過上層葉片在下層葉片發生二次反射，導致植株總體反射率急劇上升。在850～2200nm區間，反射率隨著波長增大而減小，這是因為在這一波長范圍內，光線被植物葉片中的水分和二氧化碳吸收，導致反射率減小。

圖2 不同LAI值反射率變化規律圖

2.2 波段反射率及植被指數與LAI相關關系分析

Sentinel-2多光譜遙感影像的波段反射率及各類植被指數與枸杞LAI的皮爾遜相關性分析見表5。

表5 波段反射率及植被指數與實測LAI相關關系

對比各單波段反射率與LAI的相關性，可見光波段B1～B4和B11、B12波段與LAI呈負相關；紅邊參數波段B5～B7以及近紅外波段B8和B8a與LAI均在(P<0.01)水平呈極顯著相關且相關系數均大于0.6，其中近紅外波段B8相關性表現最佳，這與已有研究證實的LAI反演中紅邊參數波段及近紅外波段為敏感波段，以及LAI與這些波段的反射率呈極顯著線性相關關系的結論一致[39]。對比各類植被指數的表現發現，本研究所選植被指數與LAI均在(P<0.01)水平呈極顯著相關，其中紅邊參數植被指數IRECI與LAI相關性最好，具有最大的相關系數。

2.3 LAI估算模型構建

隨機將樣方分為訓練集和測試集，其中訓練集60組，測試集20組。基于Sentinel-2多光譜遙感影像提取的波段反射率和計算的植被指數，本研究共建立3種數據輸入模式。模式1以波段反射率為自變量，對應實測LAI值為因變量；模式2以植被指數為自變量，對應實測LAI值為因變量；模式3以波段反射率和植被指數為自變量，對應實測LAI值為因變量。將訓練集數據驅動MLR、SML、GPR、SVM、Rtree和Etree算法，進行多次計算，分別統計各種算法決定系數(R2)最高的10次計算結果，并求算術平均數。結果見表6和圖3。

從表6、圖3可知，訓練集數據中，在3種數據輸入模式中，GPR算法均表現了較強的預測能力，決定系數(R2)分別為0.803、0.668、0.699，均高于同模型中其它算法；對應均方根誤差(RMSE)為0.402、0.453、0.441，均低于同模型中其它算法。

圖3 訓練集中各輸入模式不同算法R2和RMSE箱線圖

表6 訓練集中各輸入模式不同算法R2和RMSE統計

將訓練好的模型帶入測試集中，統計3種數據輸入模式的決定系數(R2)和均方根誤差(RMSE)，結果如表7，箱線圖如4所示。可知基于最優的GPR算法，3種數據輸入模式中，基于波段反射率的模式1具有最好的反演精度，決定系數(R2)高于其它模式，均方根誤差(RMSE)低于其它模式。

表7 測試集中GPR算法各模式R2和RMSE統計

圖4 測試集中各輸入模式R2和RMSE箱線圖

3 結論

輸入模式的選擇對農作物生理參數遙感反演的精度影響較大，各種機器學習算法預測精度差距較大，建立不同數據輸入模式和定量評價不同算法反演精度對農作物生長遙感監測意義重大。本研究利用Sentinel-2多光譜遙感衛星數據及枸杞實測LAI數據，探討了Sentinel-2多光譜衛星數據用于枸杞LAI估算的特點和精度。主要得到以下結論。

分析基于Sentinel-2多光譜數據提取的波段反射率和實測LAI值，可以得出在400～700nm可見光區域冠層反射率呈現先上升再下降趨勢，在700～850nm的近紅外區間內，反射率隨著波長增大而增大，在紅邊波段850nm左右出現最大值，形成“紅邊”現象，在850～2200nm區間，反射率隨著波長增大而減小。

基于Sentinel-2多光譜衛星數據提取的波段反射率和以此計算的植被指數，可見光波段B1～B4和B9～B12波段與LAI呈負相關；紅邊參數波段B5～B7以及近紅外波段B8和B8a與LAI呈極顯著相關，其中近紅外波段B8相關性表現最佳。對比各類植被指數的表現發現，植被指數與LAI均呈極顯著相關，其中紅邊參數植被指數IRECI與LAI相關性最好，具有最大的相關系數。

基于Sentinel-2多光譜衛星數據提取的波段反射率和以此計算的植被指數，在劃分的訓練集中，3種數據輸入模式下，GPR算法均表現較強的預測能力，具有最高的決定系數(R2)和最低的均方根誤差(RMSE)。將訓練好的模型帶入測試集中，以波段反射率為數據輸入模式具有最好的反演精度。