基于無人機(jī)多光譜影像的番茄冠層SPAD預(yù)測研究

田軍倉 楊振峰 馮克鵬 丁新軍

(1.寧夏大學(xué)土木與水利工程學(xué)院， 銀川 750021； 2.寧夏節(jié)水灌溉與水資源調(diào)控工程技術(shù)研究中心， 銀川 750021)

0 引言

番茄是我國北方主要栽培蔬菜之一，其市場需求大、經(jīng)濟(jì)效益高[1-2]。葉綠素是植被營養(yǎng)脅迫、葉片氮素含量、植物病蟲害、農(nóng)作物估產(chǎn)的重要反映因子[3-4]，精確預(yù)測番茄冠層不同垂直位置葉綠素含量可為番茄生產(chǎn)提供指導(dǎo)。近年來，已有學(xué)者應(yīng)用無人機(jī)對田塊尺度的精準(zhǔn)農(nóng)業(yè)進(jìn)行了相關(guān)研究[5]。無人機(jī)[6]遙感監(jiān)測具有靈活性高、成本低、操作簡單、作業(yè)周期短、時(shí)效性強(qiáng)等特點(diǎn)[7]，在農(nóng)業(yè)領(lǐng)域，其應(yīng)用從傳統(tǒng)的作物種植面積、區(qū)域總產(chǎn)、品種單產(chǎn)向作物健康、作物品質(zhì)、土壤濕度、農(nóng)業(yè)污染控制等更深層方向發(fā)展[8-9]，無人機(jī)搭載各類光譜傳感器技術(shù)成為精準(zhǔn)農(nóng)業(yè)的重要技術(shù)支撐。研究表明，光譜反射及透射特性與葉片及冠層葉綠素含量具有明顯的相關(guān)性[10-13]，光譜特性成為預(yù)測農(nóng)作物葉綠素含量的重要因素。

植被凈光合速率、葉綠素?zé)晒鈪?shù)及葉面積指數(shù)因葉片位置不同而呈現(xiàn)規(guī)律性變化[14]。葉片特性受各種條件的影響，在冠層垂直方向上表現(xiàn)出差異性(如光合作用)[15]。隨著生長期的變化，營養(yǎng)物質(zhì)在冠層垂直方向上轉(zhuǎn)移、聚集，以維持代謝平衡[16]。TAWANDA等[17]對多種植物冠層不同垂直位置(莖下部、中部、上部)的葉片性狀及光譜特性進(jìn)行研究，結(jié)果表明，冠層垂直位置對單株性狀的葉片光譜特性有顯著影響。謝傳奇等[18]利用高光譜成像儀對灰霉病脅迫下番茄葉片葉綠素進(jìn)行檢測。GARA等[19]通過采集不同時(shí)期的樹冠葉片反射特性，研究了垂直樹冠剖面(冠層上部、下部)生物化學(xué)和形態(tài)在PROSPECT模型反演葉綠素含量時(shí)影響模型性能的因素。丁永軍等[20]利用光譜輻射儀采集玻璃溫室環(huán)境下番茄葉片光譜信息，提取敏感波段，建立葉綠素預(yù)測模型。SCHLEMMER等[21]基于遙感數(shù)據(jù)利用近紅外綠邊和紅邊光譜信息，分析了玉米葉片和冠層氮素與葉綠素之間的關(guān)系。陳鵬等[22]通過融合光譜與紋理信息篩選特征變量，對生育期馬鈴薯冠層葉綠素含量進(jìn)行預(yù)測建模。

在農(nóng)作物葉綠素遙感反演方法的研究中，基于物理模型的反演方法通過建立冠層光譜反射率與農(nóng)作物結(jié)構(gòu)參數(shù)及生化參數(shù)之間的輻射傳輸方程[23-25]來預(yù)測葉綠素含量，但因輸入?yún)?shù)組合的不同，得到的冠層反射率存在差異，此外，建模復(fù)雜程度高，其應(yīng)用受到限制[26]。以敏感波段及植被指數(shù)與葉綠素含量建立回歸關(guān)系的經(jīng)驗(yàn)統(tǒng)計(jì)模型[27-28]具有操作簡單、可行性強(qiáng)的特點(diǎn)，因此，被廣泛應(yīng)用于農(nóng)作物葉綠素定量遙感監(jiān)測中。

冠層不同位置的葉片特性存在差異，冠層整體葉綠素預(yù)測研究忽略了冠層不同垂直位置的葉綠素含量不同，使預(yù)測結(jié)果不能精確反映冠層垂直方向上葉綠素的變化情況，而在實(shí)際生產(chǎn)中，冠層不同位置的葉綠素精準(zhǔn)預(yù)測是指導(dǎo)田間管理的重要依據(jù)。因此，本文基于無人機(jī)多光譜技術(shù)[29]，獲取番茄種植基地遙感影像，建立番茄冠層上部(TU)、中部(TM)、下部(TL)位置及整個(gè)冠層(TC)葉綠素含量預(yù)測模型。結(jié)合以往葉綠素建模方法研究[30-32]，綜合考慮建模可行性及預(yù)測效果，選取偏最小二乘法、支持向量機(jī)及BP神經(jīng)網(wǎng)絡(luò)3種方法進(jìn)行預(yù)測建模。

1 材料與方法

1.1 研究區(qū)概況

研究區(qū)位于寧夏回族自治區(qū)賀蘭縣常信鄉(xiāng)五渠村(北緯38°38′11″,東經(jīng)106°16′49″)，屬中溫帶干旱氣候帶，干旱少雨，光照條件充足，蒸發(fā)強(qiáng)烈。研究區(qū)分為Ⅰ、Ⅱ區(qū)，番茄均采用溝壟相間的種植方式，實(shí)驗(yàn)品種為2個(gè)水平：Ⅰ區(qū)為英冠218；Ⅱ區(qū)為英冠128。番茄種植行距為80 cm，株距為40 cm，株高在1.6～2.5 m之間，種植時(shí)間為2019年5月，土壤類型均為壤土，灌水量、施肥類型及施肥量等田間管理方式均無特定實(shí)驗(yàn)處理，Ⅰ區(qū)種植行為東西走向，Ⅱ區(qū)種植行為南北走向(圖1)。

圖1 研究概況Fig.1 Study area location

1.2 無人機(jī)影像數(shù)據(jù)獲取

無人機(jī)影像數(shù)據(jù)獲取時(shí)間分別為2019年7月12日(坐果期)、2019年8月6日(結(jié)果初期)和2019年9月15日(結(jié)果晚期)，選擇天氣晴朗無風(fēng)、太陽光照強(qiáng)烈的12:00—14:00采集圖像。采用eBee無人機(jī)智能系統(tǒng)(SenseFly公司,瑞士)，搭載Sequoia五通道多光譜傳感器，參數(shù)如表1所示。飛行前用灰度板輻射定標(biāo)，定標(biāo)反射率為0.18，垂直地面飛行，設(shè)置飛行高度為120 m,航向重疊率80%，旁向重疊率70%，質(zhì)量控制點(diǎn)數(shù)量為6個(gè)，獲取正射影像圖。

表1 Sequoia多光譜傳感器參數(shù)Tab.1 Sequoia multispectral camera parameters

獲取的單幅影像數(shù)據(jù)通過Pix4D軟件的農(nóng)業(yè)多光譜模板進(jìn)行拼接；將帶有POS數(shù)據(jù)的Sequoia多光譜圖像剔除后加載到軟件中，采用UTM投影方式，輸出坐標(biāo)系為CGCS2000；導(dǎo)入預(yù)先布設(shè)的控制點(diǎn)坐標(biāo)信息進(jìn)行手動刺點(diǎn)；通過點(diǎn)云加密及紋理特征匹配，最終生成完整影像。

1.3 番茄葉綠素測定

研究表明，葉綠素含量與SPAD(Soil and plant analysis development)值具有良好的相關(guān)性，SPAD測量是葉綠素?zé)o損測量的有效方法[33]。在無人機(jī)影像獲取的同時(shí)，采用SPAD-502Plus型手持式葉綠素測量儀進(jìn)行同步測量,測量樣本點(diǎn)的株高H，最小值為160 cm，最大值為246 cm，平均值為198.1 cm，標(biāo)準(zhǔn)差為0.052 cm。

依番茄植株特性，取上部(2H/3～H)、中部(H/3～2H/3)、下部(0～H/3)各5片新舊程度存在差異、均勻分布的葉片，使用SPAD-502Plus型手持式葉綠素測量儀，每個(gè)葉片測量3組SPAD值，每個(gè)部位共獲得5組數(shù)據(jù)，求平均值作為測量數(shù)據(jù)，對每個(gè)樣本點(diǎn)的TU、TM、TL測量數(shù)據(jù)求取平均值作為TC的SPAD測量值。用RTK工具記錄樣本點(diǎn)位置，樣本點(diǎn)均勻布置于研究區(qū)內(nèi)。共設(shè)70個(gè)實(shí)測樣本點(diǎn)，每個(gè)樣本點(diǎn)實(shí)測番茄TU、TM、TL、TC 4組數(shù)據(jù)，獲取開花坐果期、結(jié)果初期、結(jié)果晚期3個(gè)生育期共840組數(shù)據(jù)。整個(gè)實(shí)驗(yàn)期內(nèi)番茄SPAD值箱線圖如圖2所示。

圖2 番茄SPAD值箱線圖Fig.2 SPAD value box map of tomato

1.4 多光譜輻射定標(biāo)

eBee無人機(jī)在飛行準(zhǔn)備時(shí)，搭載Sequoia多光譜傳感器進(jìn)行室外輻射定標(biāo)。輻射定標(biāo)時(shí)，Sequoia多光譜傳感器與灰度板垂直距離為1.5～1.8 m，并確定無陰影遮擋灰度板，所獲取的單通道反射率圖像如圖3a所示。在Pix4D軟件合成多光譜影像時(shí)，取距灰度板外邊緣1/3圖像部分的數(shù)字量化值(DN),番茄多光譜DN值轉(zhuǎn)換為反射率的公式為

(1)

式中R1——目標(biāo)地物反射率

DN1——目標(biāo)地物數(shù)字量化值

DN2——灰度板數(shù)字量化值

R2——灰度板反射率，取0.18

圖3 單通道反射率圖像Fig.3 Single channel reflectivity images

1.5 反射率提取

文獻(xiàn)[34-35]表明，綠光、紅光、近紅外及紅邊波段均與葉綠素值具有明顯的相關(guān)性。在RTK測量樣點(diǎn)位置之前，設(shè)置RTK測量儀坐標(biāo)系為CGCS2000，標(biāo)高為2 m，布設(shè)RTK矩形測點(diǎn)，在矩形測點(diǎn)內(nèi)實(shí)測番茄SPAD值。將RTK記錄樣點(diǎn)坐標(biāo)輸出并在影像圖中顯示，如圖4所示。根據(jù)顯示位置，在ENVI 5.3支持下，通過RTK工具記錄樣點(diǎn)位置，繪制感興趣區(qū)域(Region of interest, RoI)，RoI尺寸為20像素×20像素，利用New Region of Interest工具的計(jì)算統(tǒng)計(jì)功能提取各波段RoI反射率數(shù)據(jù)。

1.6 植被指數(shù)選取

植被指數(shù)(Vegetation index，VI)從數(shù)學(xué)角度出發(fā)能更多地反映植被特征，從而提高葉綠素反演精度[23]。結(jié)合以往研究成果，本研究選取的植被指數(shù)如表2所示。

1.7 多光譜植被指數(shù)變化敏感度分析

光譜敏感性分析(Sensitivity analysis，SA)是植被理化參量進(jìn)行定量反演的基礎(chǔ)[42]。為進(jìn)一步研究所選取多光譜植被指數(shù)對番茄SPAD值變化的敏感度，假設(shè)多光譜植被指數(shù)為模型指標(biāo)參數(shù)，實(shí)測番茄SPAD值為不確定因素，光譜參數(shù)敏感度指數(shù)計(jì)算公式為

圖4 RoI區(qū)域的建立Fig.4 Creating of RoI regions

表2 多光譜植被指數(shù)Tab.2 Vegetation index

(2)

式中S——光譜參數(shù)敏感度指數(shù)

ΔSP——多光譜植被指數(shù)變化率

ΔSPAD——番茄實(shí)測SPAD值變化率

根據(jù)式(2)分別計(jì)算各植被指數(shù)在TU、TM、TL、TC位置上的敏感度指數(shù)S。S越接近1，說明光譜參量對番茄SPAD值的線性敏感程度越好；S越遠(yuǎn)離1，說明光譜參量對番茄SPAD值的非線性敏感程度越好。

1.8 反演建模

偏最小二乘法(Partial least squares，PLS)是WORD于1983年提出的主要用于研究因變量對于自變量回歸的多元統(tǒng)計(jì)方法[43]。

支持向量機(jī)(Support vector machine，SVR)是一種運(yùn)用統(tǒng)計(jì)學(xué)原理的機(jī)器算法[44]。本文采用徑向基(RBF)函數(shù)，通過交叉驗(yàn)證最終確定懲罰因子C為1.66，參數(shù)δ為0.5。

BP神經(jīng)網(wǎng)絡(luò)基于誤差逆向傳播算法，可實(shí)現(xiàn)正向傳播和誤差反向傳播過程[45]。文中采用輸入、隱含、輸出3層結(jié)構(gòu)，最終確定輸入層-隱含層-輸出層結(jié)構(gòu)為9-16-1，學(xué)習(xí)速率為0.01，迭代次數(shù)5 000次。

1.9 模型驗(yàn)證

為避免人為因素對選取模型和驗(yàn)證樣本的影響，本文在Matlab平臺支持下，應(yīng)用Divide函數(shù)分別將TU、TM、TL、TC 4種冠層位置實(shí)測樣本按照建模集∶驗(yàn)證集為3∶1的比例進(jìn)行多次重復(fù)分割，直至分割后的建模集和驗(yàn)證集樣本滿足時(shí)空分布均勻。最終確定建模集為158個(gè)樣本數(shù)據(jù)，驗(yàn)證集為51個(gè)樣本數(shù)據(jù)，固定建模集和驗(yàn)證集樣本數(shù)據(jù)，進(jìn)行分析建模，采用決定系數(shù)R2(Coefficient of determination)、均方根誤差(Root mean square error，RMSE)、平均相對誤差(Mean relative error，MRE)評價(jià)模型反演精度。

圖5 不同位置處番茄光譜特征參數(shù)與SPAD值相關(guān)系數(shù)矩陣Fig.5 Matrix diagram of correlation coefficient between spectral characteristic parameters and SPAD value

2 結(jié)果與分析

2.1 SPAD值與多光譜植被指數(shù)相關(guān)性分析

對表2所示的多光譜植被指數(shù)與番茄坐果期、結(jié)果初期、結(jié)果晚期在TU、TM、TL、TC 4種冠層位置上的SPAD值分別進(jìn)行相關(guān)性分析,結(jié)果如圖5所示。

由圖5a可知，TU位置上，REG、NIR、NDVI、OSAVI、CIred-edge、TVI與SPAD值的相關(guān)系數(shù)在0.60～0.68之間(P<0.01)，而RVI、CVI、MCARI與SPAD值相關(guān)系數(shù)在0.51～0.55之間(P<0.01)。由圖5b可知，TM位置上，OSAVI與SPAD值相關(guān)系數(shù)最高，為0.68，CIred-edge最低，為0.53(P<0.01)。由圖5c可知，TL位置上，REG、NIR、NDVI、TVI與SPAD值相關(guān)系數(shù)在0.51～0.53之間，OSAVI最高，為0.58(P<0.01)，MCARI、CIred-edge與SPAD值相關(guān)系數(shù)較低，分別為0.40和0.41(P<0.01)；由圖5d可知，TC位置上，除CIred-edge外，其他植被指數(shù)與SPAD值相關(guān)系數(shù)在0.62～0.72之間(P<0.01)，REG、OSAVI較高，分別為0.72和0.71。綜合比較，所選植被指數(shù)在TC位置上整體相關(guān)性較好，其次為TU和TM位置上，在 TL位置上相關(guān)程度最差；OSAVI、REG、NIR、TVI與SPAD值相關(guān)系數(shù)在4種冠層位置下最穩(wěn)定；CIred-edge在TU位置上與SPAD值相關(guān)性程度優(yōu)于TM、TL和TC位置，CVI、MCARI、RVI在TC位置與SPAD相關(guān)性程度較好。

2.2 多光譜植被指數(shù)敏感度分析

對所選取建模的9個(gè)多光譜植被指數(shù)進(jìn)行SPAD值敏感度分析，如圖6所示。由圖可知，S變化范圍在0.44～1.69之間。OSAVI 在TU、TM、 TL、TC位置上的S變化幅度最小，說明OSAVI 對SPAD值變化線性敏感程度較好，NDVI次之；MCARI在TU、TM、 TL、TC位置上的S變化幅度最大，說明MCARI 對SPAD值變化非線性敏感程度較高，CIred-edge次之。TL位置上，各植被指數(shù)的S變化范圍均高于TU、TM、TC位置。綜合來看，各植被指數(shù)在各冠層位置上線性敏感度由優(yōu)到差的順序依次為OSAVI、NDVI、REG、NIR、TVI、RV、CVI、CIred-edge、MCARI。結(jié)合相關(guān)性分析，OSAVI、REG、NIR、TVI與SPAD值的相關(guān)系數(shù)在4種冠層位置最穩(wěn)定。因此，TU、TM、TC位置上，本研究選取OSAVI、NDVI、REG、NIR 4種植被指數(shù)作為 PLS建模的自變量，TL位置上，各植被指數(shù)的S變化幅度較大，通過減少模型輸入變量以提高線性回歸程度，選取OSAVI和NDVI兩種植被指數(shù)作為PLS建模的自變量。考慮到BP、SVR建模時(shí)為非線性映射，因此以全變量(9種植被指數(shù))作為自變量進(jìn)行建模。

圖6 多光譜植被指數(shù)敏感度指數(shù)變化范圍Fig.6 Sensitivity range of spectral parameters

2.3 番茄的SPAD值反演建模

基于植被指數(shù)的番茄冠層不同垂直位置的SPAD值預(yù)測建模及驗(yàn)證結(jié)果如表3所示。由表3可知，在TU位置上，SVR模型預(yù)測結(jié)果最優(yōu)，R2、RMSE、MRE分別為0.68、2.84和0.39。在TM位置上，SVR模型的R2比PLS和BP模型分別高0.11和0.04；RMSE比PLS模型降低了0.26，比BP模型高0.18；MRE介于兩者中間。TL位置上，SVR模型R2最高，BP模型RMSE最低。TC位置上，R2由大到小依次為SVR、PLS、BP模型，PLS模型RMSE最小，為2.22，BP模型MRE最小，為0.05。同一預(yù)測模型縱向比較，TU和TC位置上R2高于TM和TL位置。在TC、TU、TM、TL位置上PLS和BP模型R2依次減小。SVR模型在TU位置上的R2最高，在TM位置上MRE最小，為0.04。PLS模型，MRE值在TC位置上比TU、TM、TL位置分別低0.30、0.09和0.18，BP模型的MRE從小到大依次為TM、TU、TC、TL。

圖7為不同預(yù)測模型的驗(yàn)證結(jié)果。由圖可知，TU、TM、TL、TC位置上，SVR模型驗(yàn)證集R2均高于PLS、BP模型。在TU、TL位置上，BP模型RMSE值均高于PLS和SVR模型，在TM和TC位置上，PLS模型RMSE值最高。同一模型的驗(yàn)證集在不同位置比較結(jié)果為，TU位置上驗(yàn)證效果最優(yōu)，TL位置上驗(yàn)證效果最差。不同模型之間比較，基于全變量的SVR模型在建模集和驗(yàn)證集的擬合結(jié)果最優(yōu)，模型較穩(wěn)定。

2.4 基于最優(yōu)預(yù)測模型的番茄SPAD值可視化制圖

根據(jù)表2中各植被指數(shù)計(jì)算公式，得到每個(gè)像素點(diǎn)的植被指數(shù)特征。根據(jù)2.3節(jié)可知，基于全變量構(gòu)建的SVR模型預(yù)測能力最優(yōu)。因此，將每個(gè)像素點(diǎn)的植被指數(shù)特征作為模型參量代入構(gòu)建的SVR預(yù)測模型中，從而得到每個(gè)像素點(diǎn)的SPAD值，將柵格分類顯示，得到番茄SPAD值可視化制圖，如圖8所示，從左到右依次為開花坐果期、結(jié)果初期、結(jié)果晚期。由圖可知，番茄冠層各部位葉片SPAD值在開花坐果期、結(jié)果初期、結(jié)果晚期總體呈現(xiàn)衰減趨勢。總體來看，與實(shí)測番茄SPAD值箱線圖(圖2)規(guī)律一致。在開花坐果期，番茄植株生長茂盛，營養(yǎng)物質(zhì)向上聚集，因此，上層葉片SPAD值是番茄縱向生長監(jiān)測的關(guān)鍵；在結(jié)果期，果實(shí)生產(chǎn)量大，且集中在中層部位，因此，中層SPAD值監(jiān)測是及時(shí)防治番茄結(jié)果期病蟲害的關(guān)鍵。通過番茄冠層各部位SPAD值可視化填圖，更清晰直觀地說明番茄冠層不同部位葉片SPAD值的差異，對更加及時(shí)準(zhǔn)確地掌握番茄植株?duì)I養(yǎng)狀況分布情況以及給灌水、施肥和病蟲害監(jiān)測等具有重要意義。

表3 不同估算模型的SPAD值建模與驗(yàn)證結(jié)果Tab.3 Modeling and verification results of different estimation models for SPAD value

2.5 討論

(1)本文研究結(jié)果顯示，無人機(jī)多光譜植被指數(shù)OSAVI、REG、NIR、NDVI與番茄冠層不同垂直位置上的SPAD值均具有良好的相關(guān)性及線性敏感度，并且隨著冠層位置降低，其與對應(yīng)的SPAD值相關(guān)性程度及線性敏感度也依次減小，原因可能是番茄株高較高，行距、株距較小，隨著冠層高度降低，葉片透光能力減弱，使得下層葉片對光譜信息的貢獻(xiàn)較小。后續(xù)研究可以考慮用傾斜攝影等方法實(shí)現(xiàn)光譜信息在垂直位置上的分層，進(jìn)而增強(qiáng)垂直地面獲取的光譜信息對冠層不同垂直位置葉片生理特性預(yù)測的精確性。

(2)SVR模型將低維的非線性關(guān)系在高維空間中線性映射，因而具有較好的魯棒性[46]。在番茄冠層不同位置上用SPAD值預(yù)測建模，SVR模型較BP神經(jīng)網(wǎng)絡(luò)和PLS模型的R2高，主要原因是SVR模型通過交叉驗(yàn)證得到的懲罰因子C和參數(shù)δ使得模型擬合情況最優(yōu)；BP模型受到隱含層不確定性和樣本數(shù)量(建模集158個(gè)樣本，預(yù)測集51個(gè)樣本)的影響，在一定程度上影響了BP模型學(xué)習(xí)及預(yù)測能力；對于多元回歸，即使添加一個(gè)不顯著的自變量，模型R2都會增大，產(chǎn)生過擬合現(xiàn)象[31]，本文通過敏感度分析篩選少數(shù)幾個(gè)PLS模型建模變量，建模變量的數(shù)量成為PLS模型的R2較低的主要原因。

3 結(jié)論

(1)開花坐果期，番茄冠層上層葉片的SPAD值高于中層和下層葉片，結(jié)果初期和結(jié)果晚期，番茄中層葉片的SPAD值高于上層和下層葉片。

(2)所選取的光譜植被指數(shù)與番茄SPAD值的相關(guān)性在TU和TC位置上優(yōu)于TM和TL位置，OSAVI、REG、NIR、TVI與SPAD值的相關(guān)性在不同冠層位置上均最高。

(3)光譜參數(shù)線性敏感度分析表明，各光譜植被指數(shù)在TU、TM、TL、TC冠層位置上線性敏感度由優(yōu)到差依次為OSAVI、NDVI、REG、NIR、TVI、RVI、CVI、CIred-edge、MCARI。

(4)同一模型不同冠層位置上的SPAD值預(yù)測的R2由大到小依次為TU、TM、 TL，表明利用冠層光譜信息預(yù)測冠層上、中、下層SPAD值的精確性依次降低。

(4)采用篩選變量的PLS模型、全變量的SVR模型及BP模型對番茄SPAD值進(jìn)行預(yù)測建模，結(jié)果表明：SVR模型結(jié)果最優(yōu)，在TU、TM、TL、TC位置上，R2分別為0.68、0.64、0.55、0.65，RMSE分別為2.84、2.75、3.51、2.26。