基于無人機多光譜數據的水稻LAI反演與應用

王 靖 彭 漪* 劉小娟 莫佳才 梁 婷

(1.武漢大學 遙感信息工程學院，武漢 430079； 2. 武漢大學 生命科學學院，武漢 430072)

水稻是我國重要糧食作物之一，2019—2020年我國水稻的消費量約為2億t，占全世界水稻消費總量的40.8%，水稻生產直接影響到我國的糧食供應與世界糧食市場的平穩[1]。葉面積指數(Leaf area index, LAI)是單位土地面積上植株單側葉片平鋪面積的總和，與作物的色素含量、碳循環、生物量、物候等生化參數密切相關，是表征作物生長狀況以及預測作物生物量的重要參數[2]。

傳統葉面積指數測量方法費時費力，便捷、高效的遙感技術助力精準農業。遙感技術依據傳感器搭載平臺可大致分為地面平臺和衛星平臺2 類。在地面平臺上，通過測量多光譜或高光譜數據計算紅邊參數、植被指數和其他光譜特征參數建立LAI反演的經驗或半經驗模型[3-5]仍是目前估算LAI的主流方法，地面平臺可獲取高空間分辨率和高時間分辨率的數據，但在測量大范圍的田塊時人工成本較高；在衛星平臺上，利用影像對LAI的估算也取得了一定的發展，如MODIS-LAI[6]和GLASS LAI[7]產品已經在大范圍得到應用。由于衛星LAI產品空間分辨率較低，無法為精準農業提供支撐。無人機兼具有高空間分辨率和高時間分辨率的特點，可以短時間內獲取較大范圍的高分辨光譜影像[8]。依據無人機影像數據反演LAI的植被指數法[9]、改進型光譜特征參數法預測LAI[10]以及使用無人機影像數據的紋理特征估計LAI[11]等方面均已有研究。

基于各種技術平臺的不同反演方法雖然可以取得較好的反演精確度，但大都只是將單個試驗區的數據分成建模集和驗證集來進行反演與精確度驗證，應用上受限于試驗區域的跨度與遙感平臺的限制，對于不同區域模型的可移植性缺乏驗證。機器學習本質上屬于基于先驗知識的統計模型，常用的植被指數經驗回歸反演LAI方法也屬于機器學習方法的一種[12]。由于機器學習模型僅學習已知樣本，并基于統計學知識對結果進行解釋，因此當環境脅迫因素發生改變，模型訓練容易過擬合或欠擬合，再加上不同模型機理的限制，會對不同模型的遷移性產生不同影響[13]。

鄂州與海南區域跨度大，兩地地理環境因素差異大，并且試驗區數據包含了不同品種的水稻。因此，為檢測機器學習模型反演LAI的可移植性，本研究擬以鄂州與海南兩地的水稻為研究對象，依據無人機平臺獲取的水稻不同生育期的多光譜影像數據，采用鄂州水稻影像光譜數據作為建模集，并采用海南試驗區的水稻影像數據作為驗證集，評估海南試驗區的水稻生長狀況，以期為水稻LAI機器學習估計模型的遷移性驗證提供依據，同時為海南地區水稻品種選擇提供參考。

1 研究區域與數據

1.1 研究區域概況

本研究選擇華中及華南地區2個不同地域、不同氣候類型的試驗區。

華中試驗區位于湖北省鄂州市試驗基地(30°22′22.27″ N, 114°45′7.03″ E)，地處長江中游南岸，屬亞熱帶季風氣候，夏季高溫多雨，適宜水稻生長。試驗區分為48個小塊，種植48種不同品種的水稻，總面積約為2 000 m2(圖1(a))。2019年5月11日育秧，6月9日移栽，移栽密度為22.5×106株/km2，施肥適量且均等。

華南試驗區位于海南省陵水黎族自治縣多品種雜交水稻試驗基地(18°31′47.10″ N，110°03′34.90″ E)，屬典型的熱帶島嶼季風型氣候，全年高溫，干濕季分明，土壤肥沃(圖1(b))。試驗區總面積為518.24 m2，分成了12個小塊，種植了珞優9348和豐兩優4號2 種水稻，并分別施用了4種不同的氮肥水平，N00、N08、N12和N16分別代表純氮0、12 000、18 000、240 000 kg/km2，氮肥分3次施用，其中基肥(50%)、幼穗分化期(25%)、抽穗期(25%)。磷肥、鉀肥均基肥施用(磷肥：P2O59 000 kg/km2；鉀肥：18 000 kg/km2)。2017年11月10日育秧，12月10日移栽，移栽密度為22.5×106株/km2。

1.2 數據獲取

1.2.1無人機影像數據

試驗采用大疆公司生產的多軸八旋翼無人機，搭載了Mini MCA 12通道數碼相機(Mini-MCA 12, Terracam Inc., Chatsworth, CA, USA)，相機每個通道1 280像素×1 024像素，各通道波長及波段寬度見表1。無人機影像數據采集盡量選擇在無云無風、中午12點左右太陽高度角近似90°的時候進行，航高約為200 m。

湖北省鄂州雜交水稻試驗基地選定13個測量日期為：2019-06-26、2019-07-02、2019-07-06、2019-07-16、2019-07-21、2019-07-26、2019-08-01、2019-08-06、2019-08-11、2019-08-17、2019-08-21、2019-08-26、2019-09-02。海南省陵水雜交水稻試驗基地選定9個測量日期為：2018-01-29、2018-02-02、2018-02-07、2018-02-12、2018-02-20、2018-03-03、2018-03-11、2018-03-25、2018-04-01。

圖1 湖北鄂州(a)與海南陵水(b)試驗田Fig.1 Experimental field in Ezhou, Hubei Province (a) and Lingshui, Hainan Province (b)

表1 無人機傳感器12通道中心波長及波段寬度Table 1 Central wavelength and band width of 12 bands of UAV sensor

1.2.2葉面積指數測量

葉面積指數的測量使用LAI-2200植物冠層分析儀(LI-COR 2010)(LI-COR Inc., Lincoln, Nebraska, USA)對水稻進行非破壞性有效的葉面積值測量。LAI-2200配備魚眼光學傳感器，其中5 個同心環的天頂角中心分別為7°、22°、38°、52°和68°，通過測量上下輻射冠層，估算5個角度的冠層光截距和透射率，進而使用Beer-Lambert定律來計算葉面積指數，是室外試驗廣泛使用的葉面積測量儀器[14]。為了避免直射陽光引起的測量誤差，測量時間在06:30—09:30和16:30—19:30進行，與MCA無人機測量日期同步。

1.3 數據處理

1.3.1幾何校正

Mini MCA多光譜相機為推掃式相機，在實驗室內對MCA相機的鏡頭進行幾何校正，使得12個波段的多光譜影像位于同一坐標系，并消除影像的鏡頭畸變，用以下公式對無人機多光譜影像進行幾何校正：

Δx=(x-x0)(k1r2+k2r4+k3r6)+p1(r2+2(x-x0)2)+2p2(x-x0)(y-y0)+α(x-x0)+β(y-y0)

(1)

Δy=(y-y0)(k1r2+k2r4+k3r6)+p2(r2+2(y-y0)2)+2p1(x-x0)(y-y0)

(2)

利用以上公式幾何糾正單波段影像，其他鏡頭根據相似關系進行配準。

1.3.2輻射校正

此外，將數字灰度值(Digital number，DN)轉化成具有實際物理意義的地物光譜反射率ρ需要進行輻射校正。利用地面標準輻射定標版的已知反射率，與影像DN值建立線性回歸模型對影像進行輻射定標，6塊定標毯在可見光至近紅外波段范圍內的反射率值依次為：0.03、0.12、0.24、0.36、0.56和0.80，根據以下公式定標：

Ri=DNi×Gaini+Offseti

(3)

式中：Ri為第i波段的地表反射率；i=1,2,…， 12；DNi為傳感器在第i波段的DN值；Gaini為第i波段的增益系數；Offseti為第i波段的偏置值。

1.3.3水稻分類

圖像預處理完成后，將每個田塊中所有水稻像元的反射率進行平均計算，將平均值作為每個田塊LAI測量值所對應的光譜反射率，因此需要先將水稻像元從背景中提取出來。計算像元NDVI，并使用K均值聚類將像元分為植被與非植被2 類，從而提取田塊中的水稻像元。利用目視勾選田塊像元，檢驗分類精確度，均達到了92%以上。

2 試驗方法

本研究將光譜數據及光譜數據經過計算得到的植被指數與測量LAI建模，采用傳統植被指數經驗模型以及機器學習的方法，完整的技術路線圖如圖2 所示。

2.1 植被指數經驗模型

植被指數提取出植被的光譜反射特征，與葉面積指數表現出強相關性。基于植被指數的經驗模型需要確定3個關鍵要素：光譜特征參數、經驗關系數和用于模型參數計算的葉面積指數[15]。本研究中使用2種回歸模型：

線性模型：y=ax+b非線性模型：y=alnx+b

式中：y為計算的植被指數；x為實測的LAI值；a和b分別為待求解的系數。

選取8種經典的植被指數參與模型反演，如表2 所示，其中ρgreen、ρred、ρrededge、ρNIR所用波段分別為550、670、720、800 nm。

2.2 機器學習方法

采用比較常見的用于多元向量擬合的機器學習方法，將不同波段組合的光譜反射率作為特征值，投入回歸算法進行訓練。使用的方法有：貝葉斯嶺(Bayesian ridge, BR)[22]、隨機森林(Random forest, RF)[23]和梯度提升回歸(Gradient boosting regression, GBR)[24]這3種常見的機器學習算法。

貝葉斯嶺回歸(BR)的核心思想是貝葉斯估計，嶺回歸則是在貝葉斯線性回歸的損失函數中加入了一個范數的懲罰項。貝葉斯嶺回歸不僅可以解決極大似然估計找那個存在的過擬合問題，而且對樣本數據的利用率是100%，同時對于樣本中存在多重共線性的情況保持較好的穩定性。

隨機森林(RF)是一種以決策樹為基分類器的集成學習算法。每棵子樹在訓練時從原始的樣本中隨機進行有放回的采樣，同時在構建決策樹時從總的特征中隨機的選擇一部分特征用于訓練。因此隨機森林中基分類器的多樣性不僅來自樣本的擾動，還來自特征的擾動，這就使得最終集成的泛化性能可以通過個體學習器之間的差異度的增加而進一步提升。

梯度提升回歸(GBR)算法的思想借鑒于梯度提升法，基本原理是根據當前模型損失函數的負梯度信息來訓練新加入的弱分類器，然后將訓練好的弱分類器以累加的形式結合到現有模型中。GBR算法不需要對數據進行縮放就可以表現的很好，而且適用于二元特征與連續特征同時存在的數據集。

圖2 技術路線Fig.2 Technology roadmap

表2 植被指數計算公式Table 2 Formula of vegetation index

2.3 模型驗證及評價標準

本研究中，鄂州試驗區共有587個有效樣本，其中抽穗前有381個樣本，滿足正態分布；海南試驗區共有215個有效樣本，皆為水稻抽穗前數據，其中珞優9348有107個樣本，豐兩優4號有108個樣本，皆滿足正態分布。本研究使用K折交叉驗證對模型反演精確度進行評估[25]，避免模型對樣本的偏重，直接估計泛化誤差。模型的精確度評價由決定系數(R2)和變異系數(CV)來表征。

3 結果與分析

3.1 鄂州試驗區建模精確度

K折交叉驗證選擇K=3，將鄂州試驗區的樣本分成3份，每份包含127個樣本點。選擇MTVI1、EVI2、OSAVI、MTCI、NDRE、CIgreen、NDVI和CIrededge這8種植被指數，分別用線性函數和非線性函數進行回歸，結果見圖3。由圖3可知，除了CIgreen、CIrededge，MTCI這3 種指數，其他指數的線性模型精確度都要低于非線性模型。非線性模型精確度最高的植被指數為EVI2(R2=0.66，CV=29.83%)；線性模型精確度最高的植被指數為MTCI(R2=0.64，CV=30.96%)。而NDVI(R2=0.25)和CIgreen(R2=0.34)兩種植被指數，線性和非線性模型的R2都小于0.4，因此后續海南試驗區驗證模型精確度不使用這兩種指數。

從1到12個波段依次增加波段數量，并列出所有可能的波段組合，然后投入機器學習算法。該窮舉法是為了找到機器學習最佳的波段組合，判斷最佳波段組合的依據是CV達到最小，結果見圖4。從圖中可以看到，4種機器學習算法在波段數量為1～4 的時候精確度迅速提升，但再隨著波段數的增加，大約在波段數>5的時候精確度不再出現明顯提升，在波段數>9的時候，RF、GBR和BR算法的擬合精確度反而出現下降的情況。

將經驗回歸模型中擬合精確度最好的EVI2非線性模型(CV=29.83%)以及MTCI線性模型(CV=30.96%)與其他機器學習算法一起進行比較，在投入相同的波段數下，機器學習模型的精確度高于植被指數回歸模型，其中在投入2、3波段數時GBR模型精確度最優，CV分別為28.46%和26.45%。

3.2 海南試驗區驗證精確度

將海南試驗區的數據輸入用鄂州數據集訓練好的EVI2、MTCI、MTVI1、CIrededge、NDRE、OSAVI等6 種植被指數的線性與非線性回歸模型，以及最佳2波段、3波段組合的BR、GBR、RF這3種種機器學習模型。共計18個模型。為了消除不同波段量綱的影響，使用CV作為評價標準，結果見表3。表3比較了鄂州試驗區的建模精確度和海南試驗區的驗證精確度：首先，OSAVI的線性與非線性回歸模型在海南試驗區CV>1 000%，已經完全不適用，一定程度上反映鄂州與海南兩地的土壤條件差異；其次，與線性模型與機器學習模型相比，非線性的植被指數模型的推廣性較差。除去MTVI1的非線性擬合模型推廣后依然較為穩定(鄂州區CV=31.96%，海南區CV=30.73%)，其余非線性植被指數模型在推廣至海南試驗區后，CV升高且超過了40%。

x代表LAI，y1代表線性模型的回歸結果，y2代表非線性模型的回歸結果。 x represents LAI, y1 represents the regression result of linear model, and y2 represents the regression result of nonlinear model.圖3 8種植被指數線性與非線性模型模型Fig.3 Linear model and nonlinear model of 8 vegetation index

將海南試驗區的估計值與真實值進行比較，結果見圖5。首先在試驗設計之初，測量日期2月12日—2月20日是水稻生長迅猛的時間，由于并未縮短測量間隔，因此造成了LAI 2～3的一段空缺，但這不影響模型的驗證。從圖5中可以發現： RF算法當LAI較高時，易出現偏低估計，且精確度點分布較散。BR算法精確度點分布則出現了偏離的情況，在LAI較低時估計偏低，在LAI較高時則估計偏高。對于植被指數的線性擬合模型來說，基本上在LAI較低的情況下都會出現估計偏低的情況，對于植被指數的非線性擬合模型，在LAI較高的情況下估計分散，且估計普遍偏高。

圖4 不同模型隨波段數變化的CVFig.4 Coefficient of variations (CV) of different models under different number of bands

表3 不同模型鄂州建模集精確度與海南驗證集精確度的比較Table 3 Comparison of the accuracy of different models Ezhou modeling set and Hainan verification set %

圖5 不同模型估計值與真實值得比較Fig.5 Comparison of the estimated values and mesured values of different models

綜上可知：各模型推廣至海南試驗區后精確度最高的是GBR二波段模型，海南試驗區的CV為26.58%(鄂州試驗區CV為28.91%)；植被指數EVI2的線性模型次之，海南試驗區的CV為27.90%(鄂州試驗區CV為33.78%)。

3.3 模型應用及評估

在區分4種施氮差異的基礎上，分別對2個品種驗證模型精確度，然后不區分品種綜合驗證模型精確度，顯示了EVI2線性經驗模型和GBR二波段模型在2種水稻品種間的精確度差異，得到表4。由表4中可知：對于同一氮水平，使用EVI2線性模型在2個不同品種間的差異最大為5.95%(N16水平下)；在品種間差異較小，而GBR模型在N16水平下，不同品種間的差異達到12.69%(N16水平下)，品種間差異較大。因此，盡管GBR模型的整體精確度(CV=26.58%)比EVI2經驗模型(CV=27.90%)高，但EVI2在品種間的差異更小，且單獨驗證不同品種的精確度與統一驗證所有品種間的差異不大，因此可以使用EVI2經驗模型推廣至海南地區的水稻LAI反演。

表4 EVI2經驗模型與GBR二波段模型在海南試驗區品種差異系數Table 4 Coefficient of variation between EVI2 empirical model and GBR two band model in Hainan experimental area %

EVI2線性模型推廣至海南試驗區后各田塊葉面積指數的分布情況見圖8。隨著水稻的生長，各田塊LAI逐漸增大，2018-02-12—2018-02-20直觀脹現水稻從分蘗期到拔節期呈現迅猛生長的態勢。

圖6 2種水稻在4種氮水平下不同日期的LAI估計值比較Fig.6 Comparison of LAI estimates of two rice varieties at different dates under four nitrogen levels

圖6是9個不同時期、2個水稻品種(珞優9348和豐兩優4號)、4種施氮條件下LAI的均值比較。海南土壤肥沃，土質本身含有豐富的氮元素，十分適宜水稻種植，在水稻分蘗期(2018-01-29—2018-02-12)，隨著氮肥量的增加，水稻葉面積指數反而出現下降的狀況；在水稻拔節期(2018-02-20—2018-03-25)，水稻葉片迅速生長，LAI迅速達到一個較高水平(LAI>8)；當水稻進入孕穗期(2018-03-25—2018-04-01)，葉鞘吸氮量降低，穗吸氮量將顯著提高，累積干物質，此時LAI將會降低。另外，在同氮水平下，珞優8348的LAI性狀預測表現明顯高于豐兩優4號。該差異在水稻分蘗期十分明顯(2018-01-29—2018-02-12)，但隨著水稻生長，差異減小，在拔節期(2018-02-20—2018-03-25)兩者差異較小，而隨著進入孕穗期(2018-04-01)，差異又迅速增大。N00、N08、N12和N16水平下，珞優9348的LAI值分別比豐兩優4號超出1.03、0.85、0.56和0.43。

圖7 EVI2線性模型反演海南試驗區水稻LAIFig.7 EVI2 linear model inversion of rice LAI in Hainan experimental area

4 討論與結論

4.1 討論

本研究使用無人機獲取的多光譜影像數據，通過植被指數的線性、非線性模型以及機器學習模型對水稻建立LAI反演模型，使用鄂州試驗區的數據建模后推廣至海南試驗區。

利用植被指數經驗模型反演LAI時，非線性模型優于線性模型，這是由于在水稻分蘗中后期，葉片覆蓋度增加越來越慢，導致冠層光譜在紅外波段的反射率增長愈加緩慢，因此LAI與植被指數之間呈現出非線性關系[21]。此外，NDVI與CIgreen的經驗模型精確度較低，是由于水稻分蘗后期葉片迅速遮蓋田地，紅波段反射率迅速下降，而近紅外波段反射率遠高于紅波段，造成NDVI迅速飽和；而在水稻拔節期，綠波段變化微小以及比值的構造使得CIgreen的樣本點非常的分散。使用機器學習模型反演LAI，由于過多的波段數存在冗余信息，不僅會增加算法復雜度，且對提升模型精確度沒有幫助。

在水稻葉片迅速生長階段，葉鞘吸氮量占水稻各器官比例最高。不同品種水稻的葉鞘的吸氮能力不同，基本水平保持在80%以上[26]。當LAI迅速達到一個較高水平(LAI>8)時，需要適量追加氮肥。本研究發現：以海南陵水試驗區的土質為基礎，珞優9348適宜添加的氮水平是N08，過量增加氮肥對水稻生長無明顯促進；豐兩優4號適宜添加的氮水平是N12，N16水平的過量氮肥則會使得水稻LAI下降。

珞優9348屬于氮高效品種，對于氮素的吸收利用效率高，抽穗期積累了大量的氮素為灌漿期提供良好基礎，因此是高產品種。從珞優9348的葉面積指數來看，全生育期的高LAI水平有助于水稻的光合作用，干物質總量增加。但珞優9348氮高效機理尚不清楚。

4.2 結論

GBR二波段模型(鄂州建模集CV=28.91%，海南驗證集CV=26.58%)與EVI的線性模型(鄂州建模集CV=33.78%，海南驗證集CV=27.90%)具有較好的推廣性；水稻種植施氮需要根據生長時期進行調整，以海南為例，分蘗期、拔節期水稻生長可以依賴土壤中的氮元素，分蘗期、孕穗期則需要追加氮肥；同等氮水平下，珞優9348是氮元素吸收利用更高效的品種。此外，水稻的葉面積指數雖然可以一定程度上反映出水稻的生長狀況，但并不能完全指示水稻產量狀況，例如有的水稻會出現徒長的情況，LAI很高但產量很低。因此，有必要進一步探究無人機多光譜數據與產量的關系，從無人機原始數據到水稻的抽穗期、LAI預測、產量預測等形成一條完整便捷的產品鏈，更好地助益農業生產。