基于無人機多光譜數(shù)據(jù)的水稻LAI反演與應用

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

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

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

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

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

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

1 研究區(qū)域與數(shù)據(jù)

1.1 研究區(qū)域概況

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

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

華南試驗區(qū)位于海南省陵水黎族自治縣多品種雜交水稻試驗基地(18°31′47.10″ N，110°03′34.90″ E)，屬典型的熱帶島嶼季風型氣候，全年高溫，干濕季分明，土壤肥沃(圖1(b))。試驗區(qū)總面積為518.24 m2，分成了12個小塊，種植了珞優(yōu)9348和豐兩優(yōu)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 數(shù)據(jù)獲取

1.2.1無人機影像數(shù)據(jù)

試驗采用大疆公司生產(chǎn)的多軸八旋翼無人機，搭載了Mini MCA 12通道數(shù)碼相機(Mini-MCA 12, Terracam Inc., Chatsworth, CA, USA)，相機每個通道1 280像素×1 024像素，各通道波長及波段寬度見表1。無人機影像數(shù)據(jù)采集盡量選擇在無云無風、中午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葉面積指數(shù)測量

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

1.3 數(shù)據(jù)處理

1.3.1幾何校正

Mini MCA多光譜相機為推掃式相機，在實驗室內(nèi)對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)

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

1.3.2輻射校正

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

Ri=DNi×Gaini+Offseti

(3)

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

1.3.3水稻分類

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

2 試驗方法

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

2.1 植被指數(shù)經(jīng)驗模型

植被指數(shù)提取出植被的光譜反射特征，與葉面積指數(shù)表現(xiàn)出強相關性?；谥脖恢笖?shù)的經(jīng)驗模型需要確定3個關鍵要素：光譜特征參數(shù)、經(jīng)驗關系數(shù)和用于模型參數(shù)計算的葉面積指數(shù)[15]。本研究中使用2種回歸模型：

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

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

選取8種經(jīng)典的植被指數(shù)參與模型反演，如表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)的核心思想是貝葉斯估計，嶺回歸則是在貝葉斯線性回歸的損失函數(shù)中加入了一個范數(shù)的懲罰項。貝葉斯嶺回歸不僅可以解決極大似然估計找那個存在的過擬合問題，而且對樣本數(shù)據(jù)的利用率是100%，同時對于樣本中存在多重共線性的情況保持較好的穩(wěn)定性。

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

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

圖2 技術路線Fig.2 Technology roadmap

表2 植被指數(shù)計算公式Table 2 Formula of vegetation index

2.3 模型驗證及評價標準

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

3 結果與分析

3.1 鄂州試驗區(qū)建模精確度

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

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

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

3.2 海南試驗區(qū)驗證精確度

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

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

圖4 不同模型隨波段數(shù)變化的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

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

3.3 模型應用及評估

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

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

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

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

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

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

4 討論與結論

4.1 討論

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

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

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

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

4.2 結論

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