無人機多光譜影像的馬鈴薯地上生物量估算

劉 楊，孫 乾，黃 玨，馮海寬，王嬌嬌，楊貴軍

1. 農業部農業遙感機理與定量遙感重點實驗室，北京農業信息技術研究中心，北京 100097 2. 山東科技大學測繪科學與工程學院，山東 青島 266590 3.南京農業大學國家信息農業工程技術中心，江蘇 南京 210095 4.國家農業信息化工程技術研究中心，北京 100097

引 言

地上生物量(above-ground biomass，AGB)是作物活力和凈初級生產力的基礎，能夠表征作物光合作用強度和評估作物營養狀況以及長勢情況[1]。故快速無損地獲取AGB信息可監測作物長勢變化，這對于指導農業生產管理和精準施肥至關重要。傳統上，AGB測量是采用破壞性取樣方法，需要人工收割作物、 稱重和記錄，這使得大面積、 長時間的監測作物變為困難[2]。

隨著遙感平臺的多樣性和影像空間分辨率的提高，遙感技術被當今認為是監測和估算大面積作物AGB最有效的技術手段[3]。其中低空無人機遙感平臺具有機動靈活，能夠提供更高的時間、 空間和光譜分辨率的遙感數據等優勢，成為了精準農業中定量觀測使用最為頻繁的技術工具[4]。當前，無人機平臺搭載的傳感器主要有3種(數碼、 多光譜和高光譜)。數碼相機雖然價格低廉，但是波段數目較少，難以解譯作物冠層光譜中包含的復雜信息[5]。高光譜傳感器在高精度表征光譜響應能力方面雖表現突出，但價格昂貴，數據后續處理較為復雜[6]。而多光譜傳感器由于經濟上更適宜，且包含監測農學參數信息重要的紅邊波段，因此在農業定量遙感中得到廣泛關注[7]。國內外學者基于無人機多光譜遙感技術開展了作物理化參數監測的大量研究，如Han等[4]使用6種植被指數結合不同的回歸技術，評估了玉米不同生育期AGB的估算效果； 肖武等[8]利用22種植被指數結合BP神經網絡有效地估算采煤沉陷區上方玉米作物的AGB； Qi等[7]使用8種植被指數通過人工神經網絡估算花生不同生育期的葉面積指數； 孫詩睿等[9]運用10種植被指數結合隨機森林估算冬小麥葉面積指數； Zheng[10]和Brinkhoff等[11]分別通過19種和6種植被指數估算冬小麥和水稻氮素含量，結果發現RDVI和NDRE植被指數估算效果最佳； 陳鵬等[12]和奚雪等[13]分別使用植被指數估算馬鈴薯和冬小麥葉綠素含量。

以上基于無人機多光譜遙感技術的研究，主要通過單個或多個植被指數估算作物不同生育期的理化參數，沒有引入新的模型變量去參與建模。然而，隨著生育期的推進，植被指數可能會達到飽和，當理化參數值較大時，會存在低估的現象。為了更全面地估算作物多生育的理化參數，一些學者通過多類別變量或者回歸技術提高了估算模型精度，較好地實現作物長勢監測[14-19]。但還未見有報道關于從無人機多光譜影像中提取高頻信息(high frequency information，HFI)和作物株高(Hdsm)并結合植被指數估算馬鈴薯不同生育期的AGB。因此，本研究基于光譜信息以及融合結構信息采用偏最小二乘回歸(partial least squares regression，PLSR)和嶺回歸(ridge regression，RR)方法構建馬鈴薯不同生育期的AGB估算模型，探究融合新的模型因子對AGB估算結果的影響，從而為提高馬鈴薯AGB監測效果提供新的方法。

1 實驗部分

1.1 研究位置與試驗設計

馬鈴薯試驗在北京市昌平區小湯山鎮國家精準農業研究示范基地進行，該區域位于北緯40°10′32″—40°10′36″，東經116°26′34″—116°26′41″，氣候類型為暖溫帶半濕潤大陸性季風氣候。試驗設計為小區完全隨機試驗，共設計密度試驗、 氮素試驗、 鉀肥試驗3個試驗區。試驗品種均為早熟的中薯5和中薯3，每個品種進行相同的控制試驗，每種試驗重復3次，每種重復進行不同程度的密度、 氮素和鉀肥處理，其中密度、 氮素和鉀肥具體處理詳情見文獻[14]。小區總計48個，每個小區面積為32.5 m2。為了更精確地獲取試驗田的數字表面模型，在試驗小區周圍均勻布控11個地面控制點(ground control point，GCP)，并用差分GPS測定其三維空間位置。

1.2 無人機多光譜數據獲取及處理

分別于天空晴朗、 無風無云的2019年4月20日(裸土期)、 2019年5月28(塊莖形成期)、 2019年6月10日(塊莖增長期)和2019年6月20日(淀粉積累期)進行無人機多光譜遙感數據采集。各次起飛地點固定和飛行航線基本保持一致，操作時間為12:00，飛行高度為30 m，航向重疊度為80%，旁向重疊度為85%。獲取影像數據前，首先采集傳感器自帶的光譜反射率校正板數據，用于多光譜影像像元亮度值的標定。采用八旋翼電動無人機搭載Parrot Sequoia農業遙感專用的4通道多光譜相機，其由光照傳感器和多光譜傳感器組成，能獲取1個1 600萬像素的RGB影像和4個120萬像素的單波段影像。多光譜相機內置的波段參數信息見表1所示。

表1 多光譜傳感器的波段參數Table 1 Band parameters of multispectral sensor

無人機多光譜影像獲取后需要預處理，處理前，首先對無人機采集的多光譜影像進行篩選，去除姿態角異常，成像存在問題的影像。將選取的多光譜影像和11個地面控制點導入到以運動結構為核心算法的Agisoft PhotoScan Professional軟件中，通過GCP的空間位置對多光譜影像作地形精糾正，各時期的校正誤差均小于2 cm?；诿芗c云生成試驗區域的數字正射影像(digital orthophoto map，DOM)和數字表面模型(digital surfacemap，DSM)。采用ENVI 5.1軟件通過校正板獲取各波段的反射率影像，并以TIF格式儲存。使用ArcGIS軟件基于劃分的小區繪制出最大面積矢量并對其編號，分別計算統計出所有感興趣區域在各波段的平均光譜反射率，將其作為各小區馬鈴薯冠層在該波段下的光譜反射率。

1.3 地面數據采集及處理

地面數據采集時，獲取了各小區馬鈴薯實測AGB和株高數據。馬鈴薯AGB和株高的測定方法見文獻[14]。

1.4 植被指數的選取

根據以往的研究成果，各生育期選取9種在AGB估算方面表現較好的多波段組合植被指數(multi-band combined vegetation indices，MVIs)和本研究提取的4個單波段植被指數(single-band vegetation indices，SVIs)用于構建馬鈴薯AGB估算模型，具體的模型參數見表2所示。

1.5 高頻信息獲取

二維離散小波變換的影像分解技術通過包含低通和高通的濾波器組對無人機多光譜影像做兩次濾波，分別得到水平方向(HL)、 垂直方向(LH)和對角線方向(LL)3種高頻信息和1種低頻信息，其中高頻信息反映的是影像的大致概貌和輪廓，與影像的真實信息接近。由于紅邊是植被反射率從紅光波段到近紅外波段快速升高變化的區域，是區分植被和地物的最顯著標志，經常被用來估算作物的理化參數[2]。因此，提取各生育期紅邊波段的3種高頻信息(REGHL，REGLH和REGLL)用來參與模型構建。

1.6 分析方法

馬鈴薯株高的提取，利用數據管理、 分析和可視化顯示為平臺的ArcGIS軟件，將馬鈴薯塊莖形成期、 塊莖增長期和淀粉積累期的DSM分別與裸土期的DSM進行差值運算，可以得到相應的作物高度模型，最后利用ROI工具提取出各生育期試驗小區的馬鈴薯平均高度Hdsm。

表2 模型參數Table 2 Model parameters

構建AGB估算模型時，若輸入變量存在嚴重的共線性問題，則會降低模型的穩定性和準確性。因此，采用2種方法(PLSR和RR)估算馬鈴薯各生育期AGB。PLSR利用了多元線性回歸、 典型相關分析和主成分分析結合為一體，可以提供一種多對多的線性回歸建模方法，能夠消除自變量之間的相關性，用較少數據來估測因變量[14]。RR是一種專用于共線性數據分析的有偏估計回歸方法，實質上是一種改良的最小二乘估計法，通過放棄最小二乘法的無偏性，以損失部分信息、 降低精度為代價獲得回歸系數更為符合實際、 更可靠的回歸方法[2]。

1.7 精度評價

各生育期提取48組數據構成樣本數據集，選取32個用來構建模型，剩余16個用來驗證模型效果。為了評估不同模型的估算效果和穩定性，選用決定系數(coefficient of determination，R2)和標準均方根誤差(normalized root mean square error, NRMSE)作為精度評價指標。

2 結果與討論

2.1 馬鈴薯株高的提取

3個生育期共提取144個馬鈴薯株高，為了評估提取的精度，將地面采集的實測株高數據和Hdsm作比較分析，R2為0.87，NRMSE為14.34%，說明提取的株高Hdsm較為可靠。

2.2 植被指數、 高頻信息和Hdsm與AGB的相關性

將表2的13種植被指數、 紅邊波段的3種高頻信息和基于無人機多光譜遙感技術提取的Hdsm分別與馬鈴薯各生育期的AGB作相關性分析，得到結果見表3所示。由表3可知，各生育期提取的模型參數均與AGB達到0.01顯著水平，相關性均從塊莖形成期到淀粉積累期先升高后降低。整體上，各模型參數與AGB的相關性由高到低依次為多波段組合植被指數、 高頻信息、 單波段植被指數和Hdsm。塊莖形成期相關性最高的模型參數為RVI，相關性系數為0.750，塊莖增長期和淀粉積累期相關性最高的模型參數均為GNDVI，相關性系數分別為0.762和0.759。

表3 植被指數和Hdsm與馬鈴薯AGB的相關性系數Table 3 The correlation coefficients of vegetationindices and Hdsm with potato AGB

2.3 基于植被指數估算AGB

為了評估原始4個單波段植被指數(SVIs)和9個多波段組合植被指數(MVIs)估算AGB的能力，分別使用PLSR和RR方法構建各生育期AGB估算模型，其結果見表4所示。由表可知，各生育期使用2種方法基于SVIs和MVIs構建的模型效果均從塊莖形成期到淀粉積累期先好后變差。使用同種方法以MVIs構建的模型精度更高，穩定性更強，且均在塊莖增長期達到最佳估算效果(PLSR： 建模R2=0.65，NRMSE=17.48%； 驗證R2=0.68，NRMSE=16.71%。RR： 建模R2=0.62，NRMSE=18.42%； 驗證R2=0.65，NRMSE=17.75%)。另外，從3個生育期的估算模型精度和穩定性來看，使用PLSR方法估算AGB的效果要優于RR方法。

2.4 基于植被指數結合HFI或Hdsm估算AGB

為了探究融合新的模型因子(HFI或Hdsm)對AGB估算結果的影響，將各生育期提取的13種植被指數(vegetation indices，VIs)分別結合紅邊波段的3種高頻信息和Hdsm使用PLSR和RR方法建立各生育期的AGB估算模型，其結果見表5所示。從表5可以看出，各生育期基于植被指數結合HFI或Hdsm利用2種方法估算AGB的效果同表4的結果保持一致，也從塊莖形成期到淀粉積累期先好后變差。相較于僅以植被指數構建的AGB估算模型(表4)，結合HFI或Hdsm明顯提高了估算模型的擬合性和穩定性(表5)，其中植被指數融合紅邊波段的高頻信息效果較優，2種方法均在塊莖增長期達到最高的估算效果(PLSR： 建模R2=0.72，NRMSE=15.44%； 驗證R2=0.74，NRMSE=15.33%。RR： 建模R2=0.67，NRMSE=16.34%； 驗證R2=0.70，NRMSE=15.62%)。通過分析各生育期的表5建模和驗證結果，發現使用PLSR方法基于植被指數結合高頻信息或Hdsm估算AGB的效果也優于RR方法。

表4 基于植被指數使用PLSR和RR估算馬鈴薯AGBTable 4 Estimation of potato AGB using partial least squares regression and ridge regression based on vegetation indices

表5 基于植被指數結合HFI或Hdsm使用PLSR和RR估算馬鈴薯AGBTable 5 Estimation of potato AGB using partial least squares regression and ridgeregression based on vegetation indices combined with HFI or Hdsm

2.5 基于植被指數結合HFI和Hdsm估算AGB

將植被指數結合HFI和Hdsm作為模型輸入參數，同樣使用PLSR和RR方法構建各生育期AGB估算模型，其結果見表6所示。從表6可以看出，整個生育期基于VIs+Hdsm+HFI使用2種方法估算AGB效果變化趨勢(先好后變差)與表4和表5結果相一致，其中基于融合所有特征為變量得到的估算結果最出色(表6)。2種方法也均在塊莖增長期達到最佳估算精度(PLSR： 建模R2=0.73，NRMSE=15.22%； 驗證R2=0.75，NRMSE=14.62%。RR： 建模R2=0.69，NRMSE=15.56%； 驗證R2=0.71，NRMSE=15.47%)。各生育期基于融合所有特征的PLSR-AGB模型的R2分別提高7.69%和2.94%，7.35%和1.38%，4.83%和3.17，RR-AGB模型的R2分別提高10%和6.45%，6.15%和2.99%，10.34%和4.9%(相較于表4)。分析表6建模和驗證結果可知，馬鈴薯3個生長期使用PLSR方法基于融合特征估算AGB效果也同樣優于RR方法。

表6 基于植被指數結合HFI和Hdsm使用PLSR和RR估算馬鈴薯AGBTable 6 Estimation of potato AGB using partial least squares regression and ridgeregression based on vegetation indices combined with HFI and Hdsm

通過馬鈴薯3個生育期的無人機多光譜影像，結合GCP生成了試驗田的DSM，提取了各生育期株高Hdsm，實測株高和Hdsm擬合的R2為0.87，證實了基于DSM提取的Hdsm效果較優，此結論與他人提取作物株高的研究基本相同[5]。但是，基于無人機多光譜影像提取的馬鈴薯株高相比于實測株高偏小，這是因為多光譜傳感器獲取的是馬鈴薯冠層數據，而地面實測的株高最高點位置在生成3D點云時被去除，導致基于DSM提取的Hdsm較??； 另外，獲取的多光譜影像中包含一些裸土像元，在提取Hdsm時也參與了運算，這同樣會使提取的株高偏低。

將4個單波段的植被指數、 9個多波段組合的植被指數、 紅邊波段的3種高頻信息和Hdsm分別與AGB作相關性分析，結果表明各生育期的模型參數與AGB的相關性均達到0.01顯著水平，這說明提取的各類參數都能夠反映作物的長勢情況[1]。3個生育期的相關性均從塊莖形成期到淀粉積累期先升高后降低，這與馬鈴薯作物自身的生長狀況相關。前期主要表現為馬鈴薯生殖器官發育，莖節和葉片不斷地增長，當進入塊莖增長期，馬鈴薯植株地上各部位的鮮重達到峰值，植被覆蓋度為整個生育期的最優時期，此階段提取的冠層光譜反射率能夠真實反映作物AGB的變化情況。而生長后期地上的同化物需要向地下塊莖轉移，基部葉片自下而上逐漸衰老變黃，馬鈴薯作物長勢變差，植被覆蓋度也明顯降低，此時提取的光譜信息并不是馬鈴薯植株冠層真實的反射率，使得上述3種光譜參數與AGB的相關性降低。

以SVIs(x1)，MVIs(x2)，VIs結合Hdsm(x3)，VIs結合HFI(x4)，以及VIs結合Hdsm和HFI(x5)為模型輸入變量，使用PLSR和RR方法估算馬鈴薯各生育期的AGB。結果發現，基于5種變量使用同樣的方法構建的模型效果變化趨勢，均從塊莖形成期到淀粉積累期先好后變差，這與模型參數和AGB的相關性變化趨勢相一致。各生育期以不同變量使用同樣的方法估算AGB的精度由高到低依次為x5>x4>x3>x2>x1，主要因為融合結構信息(HFI和Hdsm)解決了生育期效應引起的植被指數飽和性問題，提高了模型預測能力，這與劉楊等[14]得到的結論一致。相較于單波段植被指數(x1)估算模型，基于多波段組合植被指數(x2)構建的模型精度和穩定性都較優，主要因為通過多波段組合的植被指數能夠去除或者降低背景土壤對馬鈴薯冠層光譜信息的影響，增強了植被指數與AGB的敏感性(表3)，以此提高了AGB估算精度，如Bispo等[15]研究結果也認為多波段組合的植被指數估算AGB效果較優。當采用波段組合的植被指數估算AGB時，大多通過綠、 紅和近紅外波段的反射率經過波段運算得到寬波段參數[15]，忽略了紅光波段與近紅外區域的紅邊參數，紅邊是植被特有的光譜特征，在整個生育期內與作物參數敏感性較高，因此紅邊位置對于研究AGB的動態變化非常重要[18]。然而，僅僅通過植被指數估算作物不同生育期的AGB，隨著生育期的推進，植被指數會出現飽和的現象，這會造成估算AGB不準確[4]。因此，將提取紅邊波段的3種高頻信息(HFI)和作物株高(Hdsm)一起融入到植被指數中形成新的模型因子來估算各生育期的AGB，結果表明結合作物光譜信息和結構信息構建的模型精度最高，穩定性最強(表6)，這與Li等[1]和Yang等[3]研究作物理化參數結果一致，都表明融入作物結構信息能夠解決植被指數造成的低估現象。

為了減弱模型參數之間的自相關性，使用PLSR和RR方法構建各生育期AGB估算模型，探究了這2種方法估算AGB的效果。結果表明每種變量以PLSR方法構建的模型R2較大，NRMSE較小，說明此方法估算效果要優于RR方法，這與Tao等[17]和Yue等[19]研究冬小麥AGB結果一致，都證明PLSR方法估算效果較優。綜上，研究作物不同生育期AGB時，模型參數類別和建模方法的選取對構建AGB估算模型的精度影響不同。因此利用無人機多光譜遙感技術估算作物AGB時，模型參數的適用性和建模方法的篩選還需要進一步研究。此外，還需要獲取不同年限和不同地點的馬鈴薯冠層光譜數據對模型的外推性進行驗證，以便得到一個更具有代表性的AGB估算模型。

3 結 論

(1)基于無人機多光譜影像提取的Hdsm與馬鈴薯實測株高之間有較高的擬合性(R2=0.87，NRMSE=14.34%)，表明提取的Hdsm精度可靠。

(2)選取的4個單波段植被指數、 9個多波段組合的植被指數、 紅邊波段的3種高頻信息和提取的Hdsm與AGB相關性均達到0.01顯著水平，各種變量與AGB的相關性變化趨勢保持一致，均從塊莖形成期到淀粉積累期先升高后降低。塊莖形成期，相關系數最高的模型變量為RVI(r=0.750)，塊莖增長期和淀粉積累期，相關系數最高的模型變量均為GNDVI(r=0.762和0.759)。

(3)各生育期以5種變量使用相同的回歸方法估算馬鈴薯AGB，從塊莖形成期到淀粉積累期估算效果先好后變差，其中估算精度由高到低依次為x5>x4>x3>x2>x1。

(4)各生育期使用PLSR以不同變量估算AGB的效果要優于RR方法，其中在塊莖增長期使用PLSR結合x5變量估算馬鈴薯AGB，精度達到最高(建模R2=0.73，NRMSE=15.22%； 驗證R2=0.75，NRMSE=14.62%)。