基于無人機高光譜影像的玉米葉綠素含量估測

郭 松， 常慶瑞， 鄭智康， 蔣丹垚， 高一帆， 宋子怡， 姜時雨

(西北農林科技大學資源環境學院，陜西 楊凌 712100)

葉綠素作為植被葉片的主要色素，是植被光合作用的核心，其含量能反映植被的長勢,指導田間氮肥施用[1-2]。原子吸收光譜法、研磨法[3]、分光光度法[4]等傳統的植被葉綠素含量檢測方法不僅對樣品本身有著不可逆的損害，還耗時耗力；SPAD儀[5]、植物多酚葉綠素測量計[6]等新型葉綠素含量檢測儀器雖然能有效避免損害樣品，增加檢測效率，但是需要接觸測量，不適合大規模樣本和面尺度下的葉綠素含量檢測[7]。近年來，高光譜遙感技術因其具有數據量大、精度高的特點，在作物生理生化參數反演研究方面逐漸被學者們所重視。然而，地面高光譜遙感存在靈活性較低，獲取數據時人力、物力損耗較大的缺點；航天高光譜遙感對機場等基礎設施和專業人員的技術水平要求比較高；航空高光譜遙感面臨重訪周期長，受大氣影響較大，數據質量偏低的挑戰[8-9]。上述缺陷限制了高光譜遙感在作物生理生化參數反演方面的應用，而低空無人機高光譜遙感的出現，使得上述問題得到一定解決，低空無人機高光譜影像的數據質量高，受大氣影響小，較為便捷且學習成本低[10]。

目前，在基于低空無人機遙感平臺反演植被生理生化參數方面，國內外學者做了大量研究。陳鵬飛等[11]發現對高空間分辨率無人機影像剔除土壤背景和計算紋理信息，均能提升棉花氮素估測模型的精度，建模決定系數從0.33上升至0.57。陳曉凱等[12]基于隨機森林算法反演冬小麥葉面積指數(LAI)，并對無人機搭載UHD185傳感器獲取的冬小麥影像進行填圖，填圖驗證均方根誤差為0.27。陶惠林等[13]通過無人機高光譜影像數據計算不同植被指數，運用多種回歸方法建立冬小麥各生育期的長勢監測指標(GMI)估測模型，結果表明開花期的多元回歸模型精度最高，其歸一化均方根誤差為15.9%。Singhal等[14]結合核脊回歸算法和無人機影像光譜信息，反演姜黃葉片的葉綠素值，結果表明該機器學習算法具有良好的建模能力，最佳模型的決定系數達0.76。Colorado等[15]利用無人機影像提取光譜曲線，采取多元線性回歸、支持向量回歸以及神經網絡回歸估測水稻葉片氮素含量，發現機器學習算法(支持向量回歸以及神經網絡回歸)的建模效果優于多元線性回歸。綜上所述，利用低空無人機影像數據反演植被生理生化參數的研究已經較為成熟，配備的傳感器包括Rededge、UHD185等，使用的方法由傳統回歸算法向機器學習算法擴展，但是目前有關無人機搭載S185高光譜相機以及基于核極限學習機算法建模的報道較少。

本研究擬以關中地區抽雄期玉米為研究對象，使用無人機搭載S185高光譜相機獲取研究區的遙感影像，結合地面采樣點實測的葉綠素含量，構建傳統回歸模型和核極限學習機回歸模型，從中選取魯棒性、泛化性較好的玉米葉綠素含量估測模型并填圖，以期為關中地區的綠色農業發展提供技術、決策支持。

1 材料與方法

1.1 研究區概況和試驗設計

研究區位于陜西省楊陵區西北農林科技大學試驗農場(108.092°E，36.303°N)，地處關中平原中部，屬于溫帶季風氣候，四季分明，年均氣溫12.85 ℃，平均降水量630.35 mm，平均海拔480 m，試驗地土壤為塿土，水肥充足，適宜耕作。此次試驗開展于2021年6月20日，試驗共設置20個小區，小區面積30 m2(5 m×6 m)，每個小區選取3個采樣點并記錄其GPS坐標(圖1)，施用氮、磷、鉀3種肥料，其中氮肥和鉀肥均設5個水平，磷肥僅設1個水平，共10個處理，每個處理重復2次，氮肥處理N1～N5純氮施用量分別為0 kg/hm2、50 kg/hm2、100 kg/hm2、150 kg/hm2、200 kg/hm2，鉀肥處理K1～K5的K2O施用量分別為0 kg/hm2、144 kg/hm2、288 kg/hm2、432 kg/hm2、576 kg/hm2，磷肥處理P1的P2O5施用量為1 125 kg/hm2，所有肥料播種前全部施入，后期不追肥，田間管理與當地一般大田一致。10個處理分別為N1P1K2、N2P1K2、N3P1K2、N4P1K2、N5P1K2、N2K1、N2K2、N2K3、N2K4、N2K5。

圖1 研究區玉米試驗田位置Fig.1 Location of maize test field in the study area

1.2 測定項目和方法

1.2.1 無人機高光譜影像和RGB影像的獲取 無人機高光譜影像數據通過M600 Pro無人機(大疆創新科技有限公司產品)搭載S185高光譜相機(Cubert公司產品)獲取。S185高光譜相機采用畫幅式成像技術，結合了高光譜數據的精確性和快照成像的高速性，可在0.001 s內采集整個視場內的高光譜影像數據，影像光譜范圍為450～998 nm，包含138個波段，光譜采樣間隔4 nm。影像采集于2021年8月20日，此時玉米處于抽雄期，是其生長發育的關鍵時期[16]，無人機飛行時間為14∶00-16∶00，天氣晴朗無風，飛行前進行黑白板校正，設定飛行高度為50 m，飛行速度為5 m/s，航向與旁向重疊度均為80%，鏡頭垂直向下，影像地面分辨率為1.4 cm。由于獲取的無人機高光譜影像缺失地理參考，故使用自帶GPS儀的精靈4 RTK無人機(大疆創新科技有限公司產品)同步獲取研究區RGB影像，設定的飛行參數與M600 Pro無人機一致，用于高光譜影像的地理配準和幾何校正。

1.2.2 葉綠素含量測定 玉米葉片葉綠素含量使用植物多酚-葉綠素測量計Dualex scientific+TM(上海澤泉科技股份有限公司產品)測定，測定項目中的Chl值代表葉綠素含量[6]，單位為μg/cm2，測定時間為無人機影像采集后的1 h內，每個采樣點選擇3片健康、無病蟲害的冠層葉片，于葉尖、中部、基部3個位置分別測量3個Chl值，每個采樣點獲取27個Chl值，取平均值作為該采樣點的Chl值。

1.3 數據處理

在S185高光譜相機處理軟件Cubert Utils Touch中進行高光譜影像的波段融合，并輸出為單張Multiff格式的高光譜影像，使用Agisoft Photoscan軟件分別將單張Multiff影像與精靈4 RTK無人機獲取的RGB影像拼接為研究區的完整影像。在軟件ENVI 5.3中以精靈4 RTK無人機獲取的RGB影像作為基準影像，以無人機高光譜影像作為待校準影像，實施地理配準和幾何校正，使高光譜影像具有正確的地理坐標，并消除幾何變形，然后對高光譜影像進行平滑濾波以減小影像數據采集過程中的隨機噪聲，依據采樣點的GPS坐標，在高光譜影像上生成半徑為30 cm的圓形感興趣區(ROI)，通過ROI提取出各采樣點的反射光譜曲線。在前人研究結果[17]的基礎上，本次研究提取的光譜曲線范圍在450～850 nm，影像處理流程如圖2顯示。依據從小到大的原則對采集的各采樣點Chl值進行排序，通過分層抽樣法以3∶1的比例獲取建模樣本45個，驗證樣本15個。

1.4 模型建立

通過各采樣點的原始光譜以及一階導數光譜獲取建模參數。使用任意2波段組合的方式構建單因素建模參數，包括歸一化光譜指數(NDSI)、差值光譜指數(DSI)和比值光譜指數(RSI)[18]；基于連續投影算法(SPA)構建多因素建模參數。單因素模型的構建方式有線性回歸、指數回歸、冪回歸、對數回歸以及多項式回歸，多因素模型的構建方式包括多元線性回歸(MLR)以及核極限學習機回歸(KELMR)，對比各類模型的精度，從中篩選出反演效果最佳的模型。所有模型的建立以及高光譜遙感影像的填圖均在軟件Matlab的環境下實現。

Cubert Utils Touch、Agisoft Photoscan、ENVI5.3為處理數據的軟件名稱。圖2 影像處理流程Fig.2 Processing flows of images

(1)

(2)

2 結果與分析

2.1 玉米葉片Chl值光譜特征分析

2.1.1 玉米葉片Chl值基本統計特征 供試玉米葉片Chl值統計結果(表1)表明，建模集和驗證集的各項統計特征均與總數據集較為接近，3類樣本的Chl值平均值和變異系數分別為47.75～48.20、6.61%～6.88%，說明使用分層抽樣法劃分的建模集和驗證集代表性較好[19]，各集群所反映的采樣點Chl值較一致。

表1 玉米葉片葉綠素含量(Chl值)統計特征

2.1.2 Chl值與光譜反射率的相關性 圖3顯示，原始光譜反射率與Chl值相關性曲線變化較為平緩，且所有波段反射率與Chl值均呈現負相關，通過0.01相關性檢驗的波段主要位于500～624 nm以及710～850 nm，在766 nm波長處相關性最佳，相關系數為-0.61。一階導數光譜反射率與Chl值相關性曲線的走勢波動性較大，通過0.01相關性檢驗的波段在504～654 nm零星分布，在686～766 nm集中分布，相關性最佳的波長為714 nm，相關系數為0.64。整體上看，在通過0.01相關性檢驗的波長數量上，原始光譜要大于一階導數光譜，但是在相關性大小方面，一階導數的最大相關系數絕對值要高于原始光譜。

圖3 不同類型光譜反射率與Chl值的相關性Fig.3 Correlation between different types of spectral reflectance and Chl values

2.1.3 Chl值與光譜指數的相關性 根據光譜指數計算公式(表2)，分析統計Chl值與任意2波段組合光譜指數的相關性，繪制相關系數等勢圖(圖4)，選擇相關性最好的波長組合方式確定敏感光譜指數。從不同類型光譜指數來看，僅差值光譜指數(DSI)存在連片的相關性較高的區域，主要分布在700～850 nm，而歸一化光譜指數(NDSI)、比值光譜指數(RSI)的相關系數偏低且分布較為零散。從相關性最佳來看，原始光譜和一階導數光譜中均為DSI與Chl值的相關性大于NDSI和RSI，相關系數分別為0.71、0.68，表明DSI是此次研究中的最優光譜指數。

表2 光譜指數計算方法

a:PS-RSI;b:PS-NDSI;c:PS-DSI;d:FDS-RSI;e:FDS-NDSI;f:FDS-DSI。PS：原始光譜；FDS：一階導數光譜；RSI：比值光譜指數；NDSI：歸一化光譜指數；DSI：差值光譜指數。圖4 不同類型光譜指數與Chl值的相關性等勢圖Fig.4 Equipotential diagram of correlation between different types of spectral indices and Chl value

2.2 單因素模型構建及精度評價

以表2中的光譜指數作為自變量，玉米葉片Chl值為因變量，構建單因素估算模型，應用線性回歸、2次曲線、指數回歸和對數回歸4種方法，篩選出各光譜指數的最佳模型(表3)。從不同類型光譜指數來看，DSI的建模效果最好；從不同光譜類型來看，原始光譜構建模型整體優于一階導數光譜。由此可得，單因素模型中的最優模型是基于原始光譜的DSI模型，其建模、驗證R2均為0.52，MRE分別為4.09%、4.36%，使用單因素估測玉米Chl值時，應優先考慮該模型。

表3 不同光譜類型的單因素模型

2.3 多因素模型構建及精度評價

2.3.1 連續投影算法提取多因素建模參數 連續投影算法是一種特征選擇方法，通過連續投影的方式從原始光譜矩陣中選出預測響應變量的子集[20]，該算法不僅確保所篩選的光譜矢量共線性最小，還實現了對響應變量預測能力的最大化，其原理為第一次迭代時任選光譜矩陣一列賦予迭代向量，此后分別計算該迭代向量對剩余向量的投影并提取最大投影的波段，提取完成后開始第二次迭代循環，直至所篩選的波段數量等于給定數量，隨著提取波段數量的增加，每循環一次即進行一次篩選波段對響應變量的預測能力計算。以預測能力最佳為原則，分別將原始光譜和一階導數光譜中通過0.01相關性檢驗的波段反射率作為自變量，Chl值作為響應變量，提取多因素建模參數(表4)，提取出的原始光譜、一階導數光譜多因素建模參數分別為14個和8個，二者的多因素建模參數波長集中在506～766 nm，即綠光波段至近紅外波段。

表4 不同類型光譜多因素建模參數

2.3.2 多因素模型建立及精度比較 由于SPA所提取的建模參數在保證預測性最佳的同時還使得變量之間的共線性最低，故多因素模型的建立不再考慮建模參數的共線性問題，使用多元線性回歸(MLR)法、核極限學習機(KELM)算法構建多因素模型。

核極限學習機是單隱含層前饋神經網絡的一種，由Huang等[21]提出，其輸入層到隱含層的權值是隨機產生的，而隱含層到輸出層的權值是計算產生的，隨機產生的權值通過不斷迭代，能得到較好的預測結果，該算法在普通極限學習機的基礎上引入了“核”的概念，與普通神經網絡相比，提升了學習速率和泛化性，在大樣本下也能保證較快的計算速率，與隨機森林回歸相比，在樣本數量偏少時也能快速學習到數據的規律，并且預測的結果能超越給定的訓練集，但是在建模樣本向高維空間映射的過程中，由于其函數參數的生成具有一定的隨機性，往往需要迭代多次才能取得較好的結果。為了減小計算成本，本研究采用獅群優化算法(LSO)[22]對KELM迭代過程進行優化，使計算過程向著R2增大、MRE減小的方向收斂。采用徑向基函數(RBF)作為KELM的核函數，獅群數量(POP)以及迭代次數(N)分別設置為10、50，優化所得的正則化系數(C)、核函數參數(S)以及多因素建模結果(表5)顯示，與單因素模型相比，多因素模型精度有較大提升，各類多因素模型的建模R2、驗證R2均在0.57以上，高于單因素的R2最大值0.52。就建模方法而言，MLR算法的建模效果要低于KELM算法；從光譜類型來看，原始光譜的模型精度優于一階導數光譜；綜合比較可得，原始光譜的LSO-KELM模型是所有方法中的最優模型，其建模、驗證R2分別為0.73、0.70，相應MRE處于較低水平，為3.56%和3.53%。

表5 不同光譜類型下多因素回歸模型

2.4 模型精度比較

基于單因素、多因素模型中的最優模型，分別進行高光譜遙感影像Chl值反演估算，同時將反演得到的預測值與實地采樣得到的實測值進行擬合，實現模型精度比較，1∶1線設置為虛線。圖5顯示，PS-DSI模型和PS-LSO-KELM模型下的Chl值反演圖基本符合小區排列，但是外在表現有所差異。與PS-DSI模型相比，PS-LSO-KELM模型的反演圖不僅體現了玉米冠層Chl值的空間分布，還體現了玉米冠層的形態特征，Chl值的位置信息更加準確，蘊含的信息更多，可見在Chl值的遙感反演估算方面，機器學習算法優于傳統回歸算法。圖6顯示，PS-LSO-KELM模型的預測值比PS-DSI模型更加貼近擬合線，同時擬合線與1∶1線差距不大，故模型精度更高，其擬合方程R2為0.70，因此，PS-LSO-KELM模型是此次研究中的最優模型。

A：基于原始光譜的差值光譜指數模型；B：基于原始光譜的獅群算法優化的核極限學習機算法模型。Chl值：葉綠素含量。圖5 不同類型最優模型下玉米冠層Chl值空間分布圖Fig.5 Spatial distribution of Chl values in maize canopy under different types of optimal models

3 討 論

玉米冠層Chl值與S185高光譜相機的原始光譜主要在綠光波段和近紅外波段處通過0.01相關性檢驗，與一階導數光譜則主要在近紅外波段處通過0.01相關性檢驗，這與王婷婷等[6]等的研究結果基本一致。構造的光譜指數中，2種光譜下均為DSI的相關性最好，這與前人的研究結果[12,23-24]不一致，一方面是因為光譜指數本質上是一種消除背景噪聲、增強反演特征的算法，故該算法對儀器、作物種類以及光譜預處理方式等具有一定依賴性；另一方面是尺度效應，前人的光譜曲線多是地面遙感獲取，屬于葉片尺度，本研究的光譜曲線來自于低空無人機影像，屬于冠層尺度。這兩方面原因都會造成光譜指數對Chl值的敏感性差異。未來的研究應該增加關中地區冠層玉米葉片的采樣數量，同時還要分多個生育期進行試驗，基于S185高光譜相機確立不同生育期的最優光譜指數。

A：基于原始光譜的差值光譜指數模型；B：基于原始光譜的獅群算法優化的核極限學習機算法模型。x：Chl值的實測值；Y：Chl值的預測值；R2：決定系數。圖6 不同類型最優模型預測值與實測值分布Fig.6 Distribution of predicted and measured values of different types of optimal models

SPA算法屬于特征選擇的一種，與基于相關性的特征選擇算法(FCBF)相比，SPA所篩選的自變量數量可控，所建立模型相對簡單,與傳統的特征提取算法如主成分分析(PCA)等相比，SPA提取的自變量可解釋性高，所建立的模型具有一定的物理化學意義。

在處理非線性回歸問題時，機器學習算法往往比傳統回歸算法效果更佳[25]。原因在于機器學習算法是通過決策樹或模擬人類神經元工作的方式來學習數據本身的規律，相較于單一的傳統回歸算法，機器學習算法體現了群體的特點。核函數的選取以及參數的不同會影響到KELM對Chl值的擬合結果，KELM的核函數包括徑向基核、線性核、小波核等，由于不同核函數的側重點不同，也有學者提出雙核函數的想法，充實機器學習算法的實用性[26]。麻雀搜索算法、灰狼算法等群體智能算法能快速確定KELM的最佳參數，減小KELM的運行成本，本研究采用獅群算法優化KELM，以較少的迭代次數取得較好的玉米冠層Chl值擬合結果。

在本研究構建的所有模型中，多因素模型擬合效果優于單因素模型，機器學習算法優于傳統回歸算法，這與前人的研究結果[27-28]一致。實際應用時，通過無人機搭載S185高光譜相機獲取目標區域抽雄期玉米高光譜遙感影像，從中提取模型參數，應用于已經訓練好的KELM模型中即可得到玉米冠層葉片Chl值的空間分布圖。采樣點的Chl值高，說明玉米長勢較好，無需追肥；采樣點的Chl值低，說明玉米長勢較差，需要適量追肥。但是此次研究的區域局限在關中地區，研究時期限制在抽雄期，構建的模型是否適用于其他地區和其他生育期的玉米還有待測試。故未來應該在不同地區、不同生育期背景下開展低空無人機搭載S185高光譜相機反演玉米冠層Chl值的研究，在更多樣本的支持下，確立相應的敏感波段，以期未來可以通過無人機高光譜技術監測大區域的玉米Chl值。

4 結 論

本研究以陜西省關中地區抽雄期玉米為研究對象，基于玉米葉片Chl值和S185高光譜影像，構建反演玉米葉片Chl值的單因素模型、多因素模型以及獅群算法優化的核極限學習機模型并填圖，得到如下結論：第一，玉米抽雄期冠層原始光譜反演Chl值的敏感波段主要集中于綠光波段至近紅外波段，一階導數光譜反演Chl值的敏感波段主要集中于近紅外波段，二者的單因素最優建模參數均為DSI，相關系數分別為0.71、0.68。第二，連續投影算法降維效果較好，從原始光譜和一階導數光譜中分別提取了14個、8個建模參數，該算法未對光譜作修改，所提取的建模參數可解釋性較好。第三，從填圖效果以及預測值驗證結果來看，PS-LSO-KELM模型是此次研究中的最優模型，其建模R2和驗證R2分別達到了0.73、0.70，MRE分別為3.56%、3.53%。