基于無人機多光譜影像的夏玉米SPAD估算模型研究

馮 浩 楊禎婷 陳 浩 吳莉鴻 李 成 王乃江

(1.西北農林科技大學水利與建筑工程學院， 陜西楊凌 712100；2.中國科學院水利部水土保持研究所， 陜西楊凌 712100；3.西北農林科技大學中國旱區節水農業研究院， 陜西楊凌 712100；4.中國電建集團西北勘測設計研究院有限公司， 西安 710065)

0 引言

近年來，隨著精準農業的發展，農業遙感技術應用的廣度和深度不斷擴展，不僅要求能獲得傳統的產量、種植面積，還要能夠及時獲取作物的健康狀況，進而指導田間管理[1-2]。葉綠素是植物進行光合作用的主要色素，是表征植物營養脅迫、光合作用能力的重要因子，是進行作物長勢監測的重要指標之一[3-4]。葉綠素相對含量(Soil and plant analyzer development，SPAD)與葉綠素含量密切相關，可以反映植物葉片葉綠素含量，測量SPAD是獲取葉綠素含量的一種無損、便捷的方法[5-6]。

目前遙感數據源主要有地物光譜儀、無人機以及衛星。地物光譜儀主要進行點尺度數據的獲取，難以進行大范圍的觀測，同時數據采集會對農田造成一定程度的破壞。衛星遙感數據時間和空間分辨率均較大，所以衛星遙感更適合區域尺度產量、種植面積和作物長勢等的監測。相比之下，無人機具有成本低、靈活性高、機動性強、時空分辨率高、操作方便等特點，進行田塊尺度作物長勢監測有無可比擬的優勢[7]。

在利用遙感數據估算SPAD的研究中，主要有3種方法：經驗模型、物理模型以及經驗-物理耦合模型[8]，其中，以有關經驗模型的研究居多[9-14]。喬浪等[12]從無人機RGB圖像中提取了10種顏色特征和6種紋理特征，采用BP神經網絡建立玉米冠層葉綠素含量檢測模型，模型的決定系數為0.72。牛慶林等[13]將可見光植被指數與多光譜植被指數結合，利用逐步回歸和隨機森林回歸方法估算SPAD值，模型的決定系數均能達到0.88以上。

經驗模型是建立基于冠層影像得到的光譜指數、植被指數等與實測SPAD之間的數學統計關系，參數少，模型簡單，應用簡便，也能夠取得精度較高的建模結果，但缺乏物理意義，不同時間地點建模結果相差大，難以移用。近年來，有研究者嘗試將氣象數據和植被指數結合，建立分層線性模型進行產量等的估算，取得了一些成果。LI等[15]結合高光譜和氣象數據，利用分層線性模型建立了小麥產量和品質的年際可擴展預測模型，結果表明，分層線性模型是一種提高冬小麥產量和籽粒蛋白質含量預測穩定性的方法。ZHU等[16]以歸一化植被指數和氣象數據作為輸入參數，建立分層線性模型估算2016—2019年吉林省玉米產量，結果顯示，分層線性模型優于線性回歸和多元線性回歸方法。XU等[17]將植被指數與歐洲中尺度天氣預報中心氣象數據使用分層線性模型耦合，構建河北省4個小麥產區籽粒蛋白質含量預測模型，研究結果表明，使用HLM方法(R2=0.524)預測的籽粒蛋白質含量比線性模型(R2=0.286)顯示出更高的準確性。

以上研究表明，加入氣象數據可以有效提高產量和籽粒蛋白質含量的模型精度，提升經驗模型在時空上的通用性。然而，以上研究都是與一元線性或多元線性模型做對比，融合了氣象數據的分層線性模型是否比其他模型精度高需要進一步探究。同時，對于能否提升SPAD估算模型的精度也是未知的。因此，本研究將基于無人機多光譜影像提取植被指數，以篩選后的植被指數作為模型輸入變量，采用偏最小二乘法、隨機森林回歸和分層線性模型分別構建2020年夏玉米拔節期、抽雄期、灌漿期以及全生育期的SPAD估算模型，進行模型評價，選出最優估算模型，以期為SPAD估算提供參考。

1 估算模型構建

1.1 研究區概況

研究區位于陜西省楊凌區西北農林科技大學灌溉試驗站(34°20′N，108°24′E，海拔521 m)。該區域屬于暖溫帶半濕潤氣候區，多年平均氣溫為12.9℃，多年平均降水量為630 mm，降水集中在7—10月。試驗采用隨機區組設計，根據施氮量不同和是否追肥共設置7個處理：CK、N1-1、N1-2、N2-1、N2-2、N3-1、N3-2。施氮量分為4個水平，分別為CK(不施肥)、N1(105 kg/hm2)、N2(210 kg/hm2)、N3(315 kg/hm2)；N1-1、N2-1、N3-1為氮肥全部基施，N1-2、N2-2、N3-2為60%氮肥基施，剩余40%在拔節期追施。每個處理設置3個重復，共21個小區，每個小區面積為20 m2。供試玉米品種為秦龍14，生育期內不灌溉，施磷肥量90 kg/hm2。研究區位置及試驗設置如圖1所示。

圖1 研究區概況Fig.1 Location of study area and experimental design

1.2 無人機影像獲取與預處理

本研究采用大疆M600 Pro搭載RedEdge多光譜相機作為無人機多光譜影像采集系統，如圖2所示，多光譜相機的參數如表1所示。

圖2 無人機多光譜影像采集系統Fig.2 UAV for acquiring multispectral images

分別于2020年7月20日(拔節期)、2020年8月9日(抽雄期)、2020年8月26日(灌漿期)進行無人機影像采集。采集時，在晴朗無云的天氣條件進行，飛行時間選擇11：00—13：00，設置飛行高度為30 m，航向重疊度為80%，旁向重疊度為80%。無人機在起飛之前，對灰板進行拍照，用于之后的輻射定標，拍照時避免其受陰影干擾。

表1 多光譜相機參數及灰板對其中心波長的反射率Tab.1 Parameters of multispectral camera and reflectivity of calibrated reflectance panel to its center wavelength

篩選出包含試驗區的單幅影像，使用Agisoft LCC公司開發的PhotoScan軟件進行拼接，得到試驗區各波段正射影像。將拼接好的各波段影像導入ArcGIS，采用

(1)

式中Ri——地面目標第i波段反射率

Di——地面目標第i波段DN值

Rbi——灰板第i波段反射率

Dbi——灰板第i波段DN值

進行輻射定標，將圖像數字量化值(Digital number，DN)轉化為反射率，從而得到各波段反射率影像。

1.3 植被指數選取與計算

結合前人的研究成果，選取了8種植被指數用于估算夏玉米SPAD值。各個植被指數的公式及來源如表2所示。

表2 植被指數匯總Tab.2 Definitions of vegetation indices

基于各波段反射率影像計算得到植被指數。進一步，采用植被指數閾值法剔除土壤背景，通過目視解譯的方式確定每個植被指數中植被與土壤背景的閾值，經柵格計算器運算后剔除土壤背景。根據試驗小區位置和大小構建矢量文件，統計矢量范圍內植被指數的均值作為這一小區的植被指數。以上過程均在ArcGIS中完成。

1.4 田間數據獲取

在無人機影像采集的同時進行SPAD的測量。使用SPAD-502Plus便攜式葉綠素測定儀測量植株的穗位葉，在葉子中部避開葉脈隨機選取5個位置測定，取平均值作為該植株的SPAD值，每個小區選7株進行測量，其平均值作為該小區的SPAD值。不同生育期獲得的SPAD值樣本特征如表3所示。

表3 SPAD值樣本特征統計Tab.3 Descriptive statistics of SPAD values

1.5 模型構建方法

偏最小二乘法(Partial least squares，PLS)是一種通過最小化誤差的平方和尋找一組數據最佳匹配函數的方法。它將相關性分析、主成分分析和多元線性回歸等統計方法有效地結合起來，可以在自變量具有多重相關性和樣本點個數較少的條件下建立模型。

隨機森林回歸(Random forest regression，RF)是利用多個決策樹對樣本進行訓練并預測的一種機器學習算法，有很強的抗干擾能力和抗過擬合能力。同時，還具有訓練速度快、無需對輸入數據做處理等優點。

分層線性模型(Hierarchical linear model，HLM)是由不同層次的自變量解釋同一變量的一體化模型，是一種具有分層或嵌套數據結構的混合模型形式。在本研究中，采用包含兩層的分層線性模型，第1層模型是基于植被指數的SPAD估算模型，第2層模型中因變量是第1層模型的截距和斜率，氣象要素(降雨量、最高氣溫)作為自變量，具體形式為

SPADmj=ψ0j+ψ1jOSAVImj+ψ2jCIremj+emj

(2)

其中

(3)

式中ψ0j——第1層模型的截距

ψ1j、ψ2j——OSAVI、CIre的斜率

emj——隨機誤差

Tmaxj、Pj——最高氣溫、降雨量

γ00、γ10、γ20——第2層模型的截距

γ01、γ11、γ21——最高氣溫的斜率

γ02、γ12、γ22——降雨量的斜率

u0j、u1j、u2j——隨機誤差

1.6 模型精度評價

采用決定系數(Coefficient of determination，R2)、均方根誤差(Root mean square error，RMSE)和相對分析誤差(Relative prediction deviation，RPD)3個指標評價模型反演精度。其中，R2越大，RMSE越小，表明模型效果越好。而RPD對模型的評價分為3個等級：當RPD<1.4時，認為模型無法對樣本進行預測；1.4≤RPD<2時，認為模型效果一般，可以用來對樣本進行粗略預測；RPD≥2時，認為模型具有極好的預測能力。

2 結果與分析

2.1 SPAD與植被指數的相關性分析

計算得到拔節期、抽雄期和灌漿期不同試驗處理下各小區的8個植被指數，與相對應的田間實測SPAD組成數據集。分別對這3個生育期的SPAD與植被指數進行相關性分析，結果如圖3所示。

圖3 植被指數與SPAD的相關系數Fig.3 Correlation coefficients between vegetation indices and SPAD values

由圖3可知，在3個生育期，NDVI、GNDVI、RVI、OSAVI、MSR和CIre與SPAD均呈正相關關系，而MCARI、NRI與SPAD之間均是負相關關系。從圖3a中可以看出，在拔節期，NDVI、RVI、OSAVI、MSR、CIre的相關系數為0.59～0.77，在0.01水平上顯著；GNDVI、MCARI與SPAD有較低的相關性，相關系數分別為0.52和-0.49，在0.05水平上相關性顯著；NRI與SPAD的相關系數只有-0.03，關系不顯著。由圖3b可知，在抽雄期，OSAVI與SPAD的相關系數最高，為0.63；NRI最低，為-0.35。GNDVI與SPAD的相關性在0.05水平上顯著，NRI不顯著，其余在0.01水平上顯著。由圖3c可以看出，在灌漿期，NRI與SPAD的相關系數仍舊是最低的，為-0.35，相關關系不顯著；除此之外的植被指數與SPAD均具有較強的相關關系，在0.01水平上顯著，其中，CIre與SPAD的相關系數最高，為0.74。

綜合來看，NRI與SPAD的相關性較差；其余植被指數與SPAD呈現較強且穩定的相關關系，均能達到0.05水平上顯著。

2.2 反演建模

利用PLS、RF分別建立拔節期、抽雄期、灌漿期的SPAD反演模型，采用PLS、RF和HLM分別建立全生育期SPAD的估算模型，建模和驗證結果如表4所示。PLS、RF、HLM模型預測值和實測值的關系分別如圖4～6所示。

表4 不同估算模型的建模和驗證結果Tab.4 Calibration and validation results of different estimation models

圖4 基于偏最小二乘法的SPAD估算模型預測值和實測值的關系Fig.4 Relationships between predicted values of SPAD based on partial least squares method and measured values

圖5 基于隨機森林的SPAD估算模型預測值和實測值的關系Fig.5 Relationships between predicted values of SPAD based on random forest and measured values

圖6 基于分層線性模型的SPAD估算模型預測值和 實測值的關系Fig.6 Relationships between predicted values of SPAD based on hierarchical linear model and measured values

基于SPAD與植被指數的相關性分析結果，選擇NDVI、GNDVI、RVI、OSAVI、MCARI、MSR、CIre這7個植被指數進行PLS和RF建模。從表4中可以看出，同一方法對不同生育期SPAD的反演效果存在差異，同一生育期中不同方法的表現也不同。就生育期來看，在拔節期，RF模型的建模集R2為0.91，驗證集R2為0.81，均高于PLS建模集和驗證集的R2；RF模型建模集和驗證集的RMSE均低于PLS模型，分別為1.00和1.24；RF模型的RPD大于2，PLS模型的RPD為1.94，說明RF模型對SPAD具有極好的預測能力，而PLS模型效果一般，可以對SPAD進行粗略估計。在抽雄期，RF模型的建模集和驗證集R2分別比PLS模型高0.51和0.46，RMSE分別低2.12和1.97，RPD顯示RF模型具有極好的預測能力，而PLS模型無法對SPAD進行預測。在灌漿期，RF模型各個評價指標均表明RF模型精度比PLS模型高，其中，兩個模型在驗證集上表現出差不多的效果，兩個模型均為一般模型。就建模方法來看，和PLS模型相比，RF模型對各生育期的SPAD反演效果都有不同程度的提升，尤其是在抽雄期。PLS和RF模型精度均在拔節期為最高，PLS模型在抽雄期表現最差，驗證集R2只有0.35，不能用來反演SPAD，而RF模型在各生育期表現較為穩定，在灌漿期反演精度相對較低。

由于用以上7個植被指數進行HLM建模時，出現過擬合現象，建立的模型在驗證集上效果很差，所以基于SPAD與植被指數的相關性分析和RF模型的變量重要性分析，選擇OSAVI和CIre作為輸入變量建立HLM模型。在全生育期，不論是建模集還是驗證集，R2從大到小依次為RF、HLM、PLS，RMSE從小到大分別為RF、HLM、PLS，RPD從大到小為RF、HLM、PLS，說明相比之下RF模型對SPAD的反演效果最好、精度最高，但3個模型對SPAD的反演效果均一般，只能進行粗略估計。

2.3 基于最佳估算模型的SPAD可視化制圖

與分生育期模型相比，全生育期模型在各生育期的預測精度有所降低，所以將分生育期RF模型應用到各生育期，得到不同生育期試驗地各小區的SPAD空間分布圖，如圖7所示。

圖7 不同生育期SPAD空間分布圖Fig.7 SPAD spatial distribution maps in different growth stages

整體上呈現出從拔節期到抽雄期再到灌漿期SPAD持續增大的趨勢，這與實測結果一致。在不同生育期，CK處理的SPAD均最小，集中在48.5～52；不同施氮處理，N1處理的SPAD較低，而N2和N3處理的比較相近；相同施氮水平下的追肥與不追肥處理差異不明顯。圖7中信息與實際測量結果的規律一致，說明基于無人機多光譜影像建立經驗模型進行田塊尺度夏玉米SPAD估算是可行的。

3 討論

本研究基于無人機獲取的多光譜影像提取8個植被指數，經篩選后作為模型輸入變量，利用PLS、RF分別建立拔節期、抽雄期和灌漿期的SPAD反演模型，采用PLS、RF和HLM進行全生育期SPAD建模。

CIre與SPAD的相關系數波動大，且在灌漿期明顯高于其他植被指數。CIre是紅邊波段和近紅外波段構建的植被指數。首先，HORLER等[26]研究發現波長700 nm左右的反射率與葉綠素有良好的相關性，這使得CIre與SPAD具有強相關性。其次，紅邊波段對葉綠素含量的變化非常敏感[27]，在植被指數中引入紅邊波段，可以增強其對葉綠素含量變化的敏感性[28]。此外，紅邊波段參與構建的植被指數多數容易受外界環境因素的影響[29-30]。因此，CIre與SPAD的相關性又具有不穩定性。

各個生育期的最優模型均是RF模型，主要原因是隨機森林算法自身的優越性[31]。它對數據的包容力強，并且隨機選取數據和待選特征，構建多個決策樹對樣本進行訓練，最終取每個決策樹預測值的平均值作為最終結果，因此，使得模型具有較高的精度以及很強的抗干擾能力[32-33]。RF模型在很多研究中都表現較好，姚雄等[34]研究表明RF模型為馬尾松林葉面積指數的最佳估測模型，程千等[35]利用無人機多時相植被指數估測冬小麥產量，RF模型估測精度最高，且穩定性更好。偏最小二乘法和隨機森林回歸都在拔節期達到最佳，原因可能是拔節期數據分布好，變異系數小，而到夏玉米生育中后期對照組和施肥的SPAD值差距拉大，變異系數增大。

HLM模型結合了夏玉米生育期的氣象數據(降雨量、最高氣溫)和多光譜數據，就本研究來看，HLM模型介于PLS模型和RF模型之間。相比于PLS模型，HLM模型對SPAD估算精度有所提升，而和RF模型相比，HLM模型精度稍差一些。這與已有研究結果相類似，HLM模型可以一定程度提升線性模型的估算精度和穩定性[15-17]。

分層線性模型為無人機遙感數據估算SPAD提供了新的方向，同時它也存在著巨大的潛力。本研究中，HLM模型融合了降雨量和最高氣溫這兩個氣象因子，還應進一步關注其他氣象參數對SPAD的影響，例如太陽輻射，從而提高SPAD估算精度。同時，如果考慮更多的輸入參數，HLM模型的復雜性將增加，所以如何平衡好模型的精度和復雜程度也是一個值得探究的問題。

4 結論

(1)在本研究選取的8個植被指數中，除NRI之外，NDVI、OSAVI、GNDVI、RVI、MCARI、MSR、CIre與SPAD均顯著相關。其中，OSAVI、NDVI與SPAD呈現出較強且穩定的相關性，在拔節期、抽雄期和灌漿期的相關系數都較高；灌漿期，CIre表現突出，相關系數達到0.74。

(2)各個生育期的最優模型均是RF模型，在拔節期、抽雄期、灌漿期和全生育期，驗證集R2分別為0.81、0.81、0.73、0.61，RMSE分別為1.24、2.32、3.13、3.20。和PLS模型相比，RF模型在抽雄期極大提升了估算精度，PLS模型驗證集的R2、RMSE、RPD分別為0.35、4.29、1.24，RF模型的R2、RMSE、RPD分別為0.81、2.32、2.29。

(3)PLS模型在各生育期表現不穩定，特別是在抽雄期，精度最低，而RF模型在各生育期表現較為穩定，PLS和RF模型精度均在拔節期為最高。與各生育期分別建立的估算模型相比，全生育期模型的估算精度有所降低。

(4)對于SPAD估算模型，將降雨量、最高氣溫這兩個氣象因子與植被指數耦合的HLM模型可以一定程度提升線性模型的估算精度和穩定性，但其精度低于RF模型。