不同分辨率無人機數(shù)碼影像的馬鈴薯地上生物量估算研究

劉 楊，馮海寬，孫 乾，楊福芹，楊貴軍

1.農(nóng)業(yè)部農(nóng)業(yè)遙感機理與定量遙感重點實驗室，北京農(nóng)業(yè)信息技術(shù)研究中心，北京 100097 2.山東科技大學測繪科學與工程學院，山東 青島 266590 3.國家農(nóng)業(yè)信息化工程技術(shù)研究中心，北京 100097 4.北京市農(nóng)業(yè)物聯(lián)網(wǎng)工程技術(shù)研究中心，北京 100097 5.河南工程學院土木工程學院，河南 鄭州 451191

引 言

地上生物量(above ground biomass，AGB)是作物生長發(fā)育過程中進行光合作用的產(chǎn)物，能夠反映作物長勢狀況和營養(yǎng)狀態(tài)，與作物產(chǎn)量有著密切的聯(lián)系[1]。因此，快速精準地監(jiān)測AGB信息對指導農(nóng)田精準管理和挖掘生產(chǎn)潛力具有重要的意義[2]。傳統(tǒng)的AGB測定主要以田間取樣調(diào)查為主，雖能夠達到一定的估算精度，但破壞性強、非動態(tài)性和花費成本高，且局限于小面積[3-4]。近年來，將遙感技術(shù)用于反演作物的生化參數(shù)，檢測生態(tài)環(huán)境信息等方面得到廣泛關注[5]。相較于其他遙感技術(shù)，無人機遙感具有運營成本較低、操作更為簡單、可云下獲取高分辨率影像等優(yōu)勢，更適合田塊尺度觀測，是現(xiàn)代化農(nóng)業(yè)向定量精細化邁入的重要技術(shù)手段[6-7]。

目前，基于無人機遙感平臺搭載的多種傳感器(如高光譜和激光雷達)相繼用于作物AGB獲取，但因傳感器價格昂貴，制約了無人機技術(shù)在作物AGB監(jiān)測中的應用。數(shù)碼相機傳感器因價格低廉、空間分辨率較高且獲取數(shù)據(jù)后續(xù)處理較為簡單等特點，作為一種低成本獲取作物AGB的方式，已經(jīng)得到國內(nèi)外許多學者的關注。崔日鮮等[2]基于無人機數(shù)碼相機獲取冬小麥影像，通過構(gòu)建植被指數(shù)利用BP神經(jīng)網(wǎng)絡和多元線性回歸估算冬小麥AGB，顯示神經(jīng)網(wǎng)絡AGB估算模型優(yōu)于多元線性回歸模型。陶惠林等[8]基于無人機數(shù)碼影像，采用逐步回歸、偏最小二乘和隨機森林3種方法以植被指數(shù)、植被指數(shù)加實測株高和植被指數(shù)加從作物表面模型(crop surface model，CSM)中獲得的估測株高3種變量估算冬小麥AGB，結(jié)果表明每種方法中加入估測株高能夠明顯提高模型的精度。Zhang等[9]基于無人機數(shù)碼影像，采用偏最小二乘、BP神經(jīng)網(wǎng)絡、支持向量機和隨機森林4種方法以5種光學植被指數(shù)構(gòu)建冬小麥AGB估算模型，結(jié)果表明偏最小二乘法相比其他方法構(gòu)建的模型效果更優(yōu)。以上成果表明，無人機數(shù)碼影像能夠很好的對作物AGB進行有效監(jiān)測，但存在一定的不足，僅利用數(shù)碼影像的光譜信息計算植被指數(shù)進行參數(shù)估算，模型的精度容易出現(xiàn)飽和現(xiàn)象。為了解決這一問題，眾多學者通過引入作物株高和紋理特征分別結(jié)合光譜特征來估算作物參數(shù)，達到了理想的效果。相較于作物株高而言，紋理信息更具多面性，能夠反映圖像的灰度屬性和空間位置關系，將其與光譜信息結(jié)合能夠更好地減弱僅利用植被指數(shù)反演作物參數(shù)出現(xiàn)的低估現(xiàn)象，提高估算模型的適用性[10]。當前，利用數(shù)碼影像以二者融合為變量估算作物AGB的研究較少，主要以高光譜和多光譜為主。劉暢等[11]利用無人機高光譜影像，比較分析光譜特征、紋理特征和二者結(jié)合對冬小麥AGB估算精度的影響，結(jié)果表明，光譜特征和紋理特征相結(jié)合能夠解決AGB估算出現(xiàn)的過早飽和問題，而且模型估算效果優(yōu)于單一光譜特征或單一紋理特征。Zheng等[12]基于無人機多光譜影像，以光譜信息和光譜信息結(jié)合紋理特征2種類型變量對水稻AGB進行估算，結(jié)果表明以后者為輸入變量能夠更好的估測AGB。Liu等[13]基于無人機多光譜影像，采用偏最小二乘和隨機森林2種方法以光譜特征、光譜特征+紋理特征估算油菜AGB，經(jīng)比較可知2種方法以融合指標構(gòu)建的模型效果較優(yōu)。紋理特征的提取十分依靠于獲得的影像分辨率，以上成果并未考慮到不同飛行高度下所獲得的影像特征對小麥、油菜和水稻AGB估算精度的影響。因此，探討不同分辨率的無人機數(shù)碼影像對作物AGB監(jiān)測結(jié)果的影響，具有科學指導意義。

本工作利用馬鈴薯塊莖增長期不同飛行高度的無人機數(shù)碼影像，結(jié)合地上部實測生物量數(shù)據(jù)，分析不同分辨率影像下獲取的光譜特征和紋理特征對馬鈴薯植株地上部生物量估算結(jié)果的影響，為無人機遙感監(jiān)測AGB時，對飛行高度和影像分辨率的選擇提供技術(shù)參考。

1 實驗部分

1.1 研究區(qū)域與實驗設計

在北京市昌平區(qū)小湯山鎮(zhèn)國家精準農(nóng)業(yè)研究示范基地的馬鈴薯試驗田開展試驗，該處位于北緯40°10′34″，東經(jīng)116°26′39″，平均海拔為36 m，氣候類型為暖溫帶半濕潤大陸性季風氣候。試驗為小區(qū)隨機試驗設計，采用了中薯5(P1)和中薯3(P2)這2個均為早熟的馬鈴薯品種，試驗田共計48個小區(qū)，每個小區(qū)面積為6.5 m×5 m。試驗區(qū)共設密度試驗(N區(qū))、氮素試驗(S區(qū))、鉀肥試驗(K區(qū))3個試驗區(qū)，每種試驗3種重復，每種重復進行不同程度的密度、氮素和鉀肥處理，其中密度設置了3種水平(T1：60 000株·hm-2、T2：72 000株·hm-2、T3：84 000株·hm-2)，氮素設置了4種水平(N0：0 kg尿素·hm-2、N1：244.65 kg尿素·hm-2、N2：489.15 kg尿素·hm-2(正常處理，15 kg純氮)、N3：733.5 kg尿素·hm-2)，鉀肥設置了3種水平(K0：0 kg鉀肥·hm-2、K1：970.5 kg鉀肥·hm-2(N區(qū)和S區(qū)均為K1處理)、K2：1 941 kg鉀肥·hm-2)。為了更好的獲取試驗田的位置，在試驗小區(qū)周圍均勻布控11個GCP(k1—k11由0.3 m×0.3 m的木板和埋于地下的木樁組成,其上有黑白標志的塑料板，目的是準確確定木板的中心位置)并用差分GPS測定其三維空間位置，詳細的試驗設計見圖1所示。

圖1 馬鈴薯研究區(qū)域和試驗設計Fig.1 The experimental area and design for potato study

1.2 地面數(shù)據(jù)獲取

馬鈴薯AGB通過收獲法獲取，在每個小區(qū)中選取能夠代表小區(qū)長勢水平的3棵植株，進行取樣并迅速帶回實驗室將其莖葉分離，隨后用流水洗凈，105 ℃殺青，80 ℃烘干48 h以上，直到質(zhì)量恒定再進行稱量。將植株莖和葉的干質(zhì)量求和得到樣本干質(zhì)量，最終通過群體密度和樣本干質(zhì)量，得到每個小區(qū)的馬鈴薯地上生物量。計算公式如式(1)

AGB=[(莖干重+葉干重)/取樣株數(shù)]×種植密度

(1)

1.3 無人機數(shù)碼影像獲取

采用大疆公司生產(chǎn)的精靈系列4A產(chǎn)品在馬鈴薯塊莖增長期進行無人機遙感數(shù)據(jù)獲取，無人機參數(shù)見表1所示。無人機飛行時間為中午12:00，此時太陽光照強度穩(wěn)定。人為設置了10，20，30，40和50 m共5種無人機飛行高度，相應的地面分辨率約為0.43，0.85，1.30，1.83和2.05 cm。

表1 無人機部分參數(shù)Table 1 Partial parameters of UAV

1.4 影像預處理

通過利用俄羅斯Agisoft LLC公司研發(fā)的基于運動三維重建(structure from motion，SFM)技術(shù)的PhotoScan 掃描軟件對獲取的不同高度的數(shù)碼影像進行拼接處理，得到馬鈴薯試驗區(qū)域的數(shù)字正射影像(digital orthophoto map，DOM)。首先將帶有經(jīng)緯度、高度、偏航角、俯仰角和旋轉(zhuǎn)角共6種POS信息的每張影像導入軟件中，還原影像拍攝時的姿態(tài)和空間位置；將GCP添加到軟件中，利用GCP的三維坐標信息優(yōu)化照片對齊，生成飛行區(qū)域的稀疏點云；利用稀疏點云構(gòu)建密集點云，進而生成飛行區(qū)域的空間格網(wǎng)并建立紋理信息；生成馬鈴薯試驗田的DOM。

1.5 影像指數(shù)選取

利用數(shù)碼影像指數(shù)可以較好的反映作物的AGB信息，因此，通過DOM，利用ENVI軟件提取每個試驗小區(qū)影像的平均紅波段值R、綠波段值G和藍波段值B。并將R，G和B三個波段的DN值進行歸一化處理，分別得到r，g和b。計算公式如下

r=R/(R+G+B)

(2)

g=G/(R+G+B)

(3)

b=B/(R+G+B)

(4)

根據(jù)文獻中選取的能夠很好反映AGB信息的15個數(shù)碼影像指數(shù)，加上定義的3個，總共18個影像指數(shù)用來估算馬鈴薯地上部生物量，具體如表2所示。

表2 數(shù)碼影像指數(shù)Table 2 Digital image indexes

1.6 紋理特征選取

紋理特征通過像素與空間鄰域的灰度分布變化來表現(xiàn)，能夠反映目標地物與周圍環(huán)境的空間位置關系，是遙感影像中一種重要的圖像信息。通過ENVI軟件，利用灰度共生矩陣提取紅、綠、藍3個波段的0°，45°，90°和135°4個方向的紋理特征，不同方向的紋理特征分別求平均值，得到每個波段的8個紋理特征，分別為對比度(contrast，con)、相關性(correlation，cor)、差異性(dissimilarity，dis)、熵(entropy，ent)、同質(zhì)性(homogenetity，hom)、二階矩(second moment，sm)、方差(variance，var)和均值(mean)。紅波段的紋理特征為con_R，cor_R，dis_R，ent_R，hom_R，sm_R，var_R和mean_R。綠波段的紋理特征為con_G，cor_G，dis_G，ent_G，hom_G，sm_G，var_G和mean_G。藍波段的紋理特征為con_B，cor_B，dis_B，ent_B，hom_B，sm_B，var_B和mean_B。

1.7 數(shù)據(jù)分析方法

1.7.1 多重共線性分析

采用方差膨脹因子(variance inflation factor,VIF)評價輸入變量之間存在的自相關性，當方差膨脹因子值較大時，輸入變量之間存在高度自相關，會造成數(shù)據(jù)冗余和計算困難，因此有必要對數(shù)據(jù)進行降維。一般當VIF<10可視為輸入變量之間不存在自相關性；當10≤VIF≤20時，認為輸入變量之間存在一定的自相關性；當VIF>20時，認為輸入變量之間自相關性較為嚴重。

1.7.2 主成分分析

主成分分析(principal components analysis，PCA)，又稱K-L變換，本質(zhì)上是對數(shù)據(jù)進行降維，即用很少的不相關變量來表達原始信息，使信息更簡單化，是當前應用最為廣泛的數(shù)據(jù)降維方法之一。現(xiàn)實中，為了能描述研究問題，往往選擇多個自變量，變量之間卻不完全獨立，導致所表達的信息重疊，為了簡化信息和提升運算效率，需要對數(shù)據(jù)進行降維處理。以光譜和紋理特征構(gòu)建AGB估算模型前，首先對光譜信息和紋理特征進行K-L變換，依據(jù)特征值大于1的原則，分別選取前n個主成分來建立估算模型。

1.7.3 多元線性回歸

多元線性回歸(multiple linear regression，MLR)可以同時使用多個自變量來預測因變量，要求自變量與因變量之間要有較好的相關性。分別以光譜特征、紋理特征以及二者結(jié)合作為自變量，馬鈴薯AGB作為因變量來綜合構(gòu)建模型。

1.8 精度評價指標

挑選2/3樣本數(shù)據(jù)(32個，重復1和重復2)作為建模集，1/3樣本數(shù)據(jù)(16個，重復3)作為驗證集以此來構(gòu)建馬鈴薯AGB估算模型。將采用決定系數(shù)(coefficient of determination,R2)、均方根誤差(root mean square error,RMSE)、標準均方根誤差(normalized root mean square error,NRMSE)來評價模型的精度。R2越接近于1，RMSE和NRMSE越低，其估測模型的精度就越高。

2 結(jié)果與討論

2.1 影像指數(shù)與AGB的相關性

將表2選取的影像指數(shù)與馬鈴薯AGB進行相關性分析，將相關系數(shù)絕對值按從大到小排列，選取不同飛行高度下前10個影像指數(shù)作為模型變量，其結(jié)果見表3。由表可知，不同飛行高度下，前10個影像指數(shù)與馬鈴薯AGB均達到0.01顯著水平(p<0.01)，相關系數(shù)絕對值都在0.606以上。隨著飛行高度的增加，VARI，GRVI，NDI，MGRVI，EXGR，g，r+b和EXG共8個影像指數(shù)與AGB的相關性均在前10以內(nèi)，相關性順序略有變動。當飛行高度為10和30 m時，影像指數(shù)EXR和RGBVI相關性都位于前10；當飛行高度為20和50 m時，影像指數(shù)RGBVI和g-b相關性也都位于前10；而當飛行高度為40 m時，影像指數(shù)g-b和EXR的相關性雖位于前10但都處于后兩位。因此可以看出，不同飛行高度下，與AGB相關性較好的前10個影像指數(shù)類別差異小，不同分辨率的數(shù)碼影像對光譜信息組成的影像指數(shù)與AGB的相關性排序差別不大。

表3 不同高度下影像指數(shù)和AGB的相關性Table 3 Correlation between image indexes of different altitudes and AGB

2.2 紋理特征與AGB的相關性

將提取的不同飛行高度下的影像紋理特征與AGB進行相關性分析，將相關系數(shù)絕對值按從大到小排列，選取前8個紋理特征作為模型變量，其結(jié)果見表4。由表可知，不同飛行高度下，前8個紋理特征與AGB均達到0.01顯著水平(p<0.01)，相關系數(shù)絕對值在0.366以上。隨著飛行高度的增加，篩選的前8個紋理特征類別差異較大，同種紋理特征在不同飛行高度下與AGB的相關性排序波動明顯。當飛行高度為50 m時，紋理特征與AGB相關系數(shù)相較于前4種高度比較小，可知紋理特征的提取非常依賴于數(shù)碼影像的分辨率。

表4 不同高度下紋理特征和AGB的相關性Table 4 Correlation between texture features of different altitudes and AGB

2.3 不同飛行高度下光譜特征估算馬鈴薯AGB

在構(gòu)建AGB估算模型前，首先對每種飛行高度下的前10個影像指數(shù)進行共線性診斷，判斷輸入變量之間的共線性。依據(jù)各個飛行高度下的共線性診斷VIF>20，得知輸入變量之間存在嚴重共線性，因此有必要對前10個影像指數(shù)進行PCA處理。根據(jù)主成分特征值大于1的原則，確定不同飛行高度下的最佳主成分為1，1，1，1和2個?；赑CA-MLR方法，得到不同飛行高度下利用光譜特征估算馬鈴薯AGB的結(jié)果如表5所示。由表可知，在飛行高度10～50 m之間，隨著影像分辨率降低，光譜信息估算馬鈴薯AGB的效果整體上呈現(xiàn)下降趨勢變化，建模R2從0.69遞減到0.65，建模RMSE從226.07 kg·hm-2遞增到240.69 kg·hm-2，建模NRMSE從21.22%遞增到25.49%；模型驗證效果同樣呈下降趨勢，R2從0.73遞減到0.70，RMSE從199.21 kg·hm-2遞增到239.25 kg·hm-2，NRMSE從17.74%遞增到18.89%。

表5 不同高度下影像光譜信息的AGB估算結(jié)果Table 5 Prediction results of AGB from image spectral information with different altitudes

2.4 不同飛行高度下紋理特征估算馬鈴薯AGB

依據(jù)各個飛行高度下的共線性診斷VIF>20，得知輸入的紋理特征之間存在嚴重共線性，因此通過PCA降維處理，確定不同飛行高度下的最佳主成分為1，1，1，1和3個?；赑CA-MLR方法，得到不同飛行高度下利用紋理特征估算馬鈴薯AGB的結(jié)果如表6所示。由表可知，隨著飛行高度的增加，紋理特征估算馬鈴薯AGB的效果逐漸變差，建模R2從0.75遞減到0.68，RMSE從208.86 kg·hm-2遞增到230.69 kg·hm-2，NRMSE從17.22%遞增到21.44%；驗證R2從0.83遞減到0.73，RMSE從161.65 kg·hm-2遞增到201.26 kg·hm-2，NRMSE從16.39%遞增到18.99%。以上結(jié)果說明，飛行高度的增加會使影像分辨率降低，進而影響紋理特征的提取，造成馬鈴薯AGB的估算效果不斷變差。所有模型中，當飛行高度為10 m及對應的影像分辨率為0.43 cm時，估算效果最好。

表6 不同高度下影像紋理特征對生物量的反演結(jié)果Table 6 Prediction results of biomass from image texture features with different altitudes

2.5 不同飛行高度下光譜與紋理特征結(jié)合估算馬鈴薯AGB

將整合后的輸入變量進行PCA處理，確定不同飛行高度下的最佳主成分為1，1，1，2和3個?；赑CA-MLR方法，得到不同飛行高度下利用光譜信息+紋理特征估算馬鈴薯AGB的結(jié)果如表7所示。由表可知，隨著飛行高度的增加，光譜信息+紋理特征估算AGB效果逐漸變差，建模R2從0.82遞減到0.70，RMSE從183.05 kg·hm-2遞增到225.85 kg·hm-2，NRMSE從15.61%遞增到19.64%；同樣驗證R2從0.85遞減到0.75，RMSE從146.51 kg·hm-2遞增到184.29 kg·hm-2，NRMSE從14.36%遞增到18.11%。當無人機飛行高度為10 m及對應的地面分辨率為0.43 cm時，估算效果最好。此外，與表5和表6結(jié)果相比，同種高度下，通過融合指標構(gòu)建的AGB估算模型精度最高，效果最佳。

表7 不同高度下影像光譜信息+紋理特征對AGB的估算結(jié)果Table 7 Prediction results of AGB from image spectral information and texture features with different altitudes

本研究構(gòu)建AGB估算模型前，發(fā)現(xiàn)當飛行高度為50 m及對應影像分辨率為2.05 cm時，紋理特征與AGB的相關性較低，主要原因是新生葉片的長度和寬度小于2 cm，此分辨率下難以獲取完整的影像紋理特征。以光譜信息、紋理特征以及二者結(jié)合分別估算馬鈴薯AGB的結(jié)果表明：同種高度下光譜信息+紋理特征估算效果最優(yōu)，單一紋理特征次之，單一光譜信息最差。此研究結(jié)果與劉暢等[11]估算冬小麥AGB時，得出以光譜信息和紋理特征融合估算精度最高，單一紋理特征次之，單一光譜最差結(jié)論一致。與陳鵬等[15]估算馬鈴薯葉綠素含量時，得出結(jié)論略有不同，相同之處都是光譜信息+紋理特征融合指標估算精度高于單一光譜、紋理特征指標，不同在于，本工作得出光譜信息估算效果最差，而陳鵬等得出紋理特征估算效果最差。主要原因可能是二者估算的作物理化參數(shù)不同，其次是采用的傳感器和影像分辨率不同。在探究不同飛行高度下數(shù)碼影像對馬鈴薯AGB估算精度的影響時，本研究人為設置了共5種飛行高度并獲取相應地影像分辨率，增加了3種變量對AGB估算能力的分析，結(jié)果發(fā)現(xiàn)同種飛行高度下光譜信息+紋理特征估算AGB的精度相較于單一光譜和紋理特征分別提高了18.84%，11.76%，10.45%，5.97%，7.69%和9.33%，4.13%，2.78%，2.89%，2.94%，表明僅利用數(shù)碼影像的光譜信息估算馬鈴薯AGB，模型的精度容易出現(xiàn)飽和現(xiàn)象，而影像的紋理特征不受環(huán)境的影響，所得到的模型效果同種分辨率下優(yōu)于光譜信息，綜合考慮了光譜和紋理信息對AGB估算的貢獻性，從而結(jié)合二者的優(yōu)勢，在一定程度上彌補了單一指標估算AGB能力的不足，提高了模型適用性。本研究還發(fā)現(xiàn)無人機數(shù)碼影像分辨率在0.43～2.05 cm之間變化，光譜信息、紋理信息以及光譜+紋理信息會隨影像分辨率提高，其估算AGB的精度逐漸變好，這與賈丹等[10]的研究結(jié)果略有不同。賈丹等研究結(jié)果表明，隨著影像分辨率提高，光譜信息對氮素濃度反演結(jié)果影響不大，而基于紋理信息以及光譜+紋理信息的氮濃度反演精度隨分辨率增強而提高，主要原因其一可能是實驗設置的飛行高度不同，其二研究所用的傳感器和作物對象不同，賈丹等研究利用多光譜估算冬小麥氮素濃度，其三本文通過設計氮素和密度梯度試驗，會導致每個小區(qū)的植被覆蓋度發(fā)生變化，植株葉片的表現(xiàn)形態(tài)和冠層空間結(jié)構(gòu)會明顯不同，這會導致作物紋理特征和冠層光譜有所差別。綜上可知，對于同種傳感器來講，要想獲得高分辨率影像，必須降低無人機飛行高度，但是這樣會降低影像采集效率。因此，在對AGB進行監(jiān)測時，需要根據(jù)現(xiàn)實需求并考慮模型精度來選擇一定的飛行高度。本研究也存在一些缺陷，僅使用一年的馬鈴薯塊莖增長期的地上部生物量數(shù)據(jù)進行探討，得到的模型適用性還需要進一步驗證，今后應采集更多年份、生育期及不同作物的樣品數(shù)據(jù)，對影像分辨率如何影響AGB估算結(jié)果進行深入分析。

3 結(jié) 論

(1)隨著數(shù)碼影像分辨率降低，與馬鈴薯AGB相關的前10個植被指數(shù)類別和相關性排序差異不大；與馬鈴薯AGB相關的前8個紋理特征類別有所差異且相關性排序波動明顯，飛行高度為50 m時，影像的紋理特征與AGB相關性最低。

(2)同種分辨率影像下，光譜信息+紋理特征估算AGB的效果最優(yōu)，其次為單一紋理特征模型，而單一光譜模型表現(xiàn)最差。

(3)隨著數(shù)碼影像分辨率提高，光譜信息、紋理信息以及光譜+紋理信息估算AGB的精度逐漸變好。