基于高光譜特征的松材線蟲嶺回歸估測模型研究

張素蘭 黃金龍 秦 林 李宏群

(1.長江師范學院大數據與智能工程學院， 重慶 408100； 2.昆士蘭大學園藝科學中心， 布里斯班 4072；3.長江師范學院電子信息工程學院， 重慶 408100； 4.長江師范學院三峽庫區環境監測與災害防治工程研究中心， 重慶 408100)

0 引言

松材線蟲病是由松材線蟲寄生在松樹上引起的毀滅性病害，1982年首次傳入中國，并被確定為我國主要的森林病蟲害[1-2]。松樹感染初期，感病枝梢針葉由綠變黃再到紅褐色，并由局部發展到整樹針葉出現萎蔫，直到全株枯萎死亡。松樹一旦感染松材線蟲病，整株松樹最快40 d左右即可枯死，成片松林從最初少數死樹到林相被毀只需5年左右[3]。因此，及時準確地掌握松材線蟲病發生的地理位置、面積、危害等級等，對于松材線蟲病的防治尤為重要[4]。目前監測技術手段仍是現場觀察，費時費力，且獲得的信息具有滯后性，使得松材線蟲病很難徹底控制和根除，嚴重破壞了森林生態系統[5]。

近年來，隨著信息技術的發展，高光譜遙感技術作為一種新的監測手段被廣泛運用于植物生長監測中[6-8]。當植物受到病蟲危害時，綠葉中細胞活性、含水率和葉綠素等生理指標都會變化，引起反射光譜特性上的差異，特別是紅色區和近紅外區的光譜特征的差異[9]。因此，基于高光譜遙感數據，分析病蟲害危害程度與原始光譜、光譜特征、植被指數等之間的關系，確定不同種類病蟲害監測的敏感波段，建立反演模型，可對森林病蟲害進行有效監測[10]。NASI等[11]利用光譜指數和K近鄰分類器估測云杉受皮甲蟲危害程度，對于判別健康和死亡兩種類型的總體準確率可達90%。CALDERON等[12]根據生化參數和病害指數，利用支持向量機(SVM)和線性判別分析(LDA)分類方法對橄欖樹黃萎病進行分類，其中LDA在初始和低嚴重程度上分類效果較優，分類準確率可達到71.4%和75.0%。

在松材線蟲病害高光譜遙感監測方面，國內外學者主要研究了敏感波段、植被指數及光譜特征檢測松材線蟲病害[13-16]。KIM等[13]基于現有10個高光譜植被指數，引入了綠邊與紅邊面積指數，判斷是否感染松材線蟲病。JU等[14]研究了健康馬尾松感染松材線蟲病的最佳光譜波段或波段組合，發現759 nm處的波譜的一階導數為區別馬尾松是否染病的最有效波段。黃明祥等[15]對受松材線蟲感染的馬尾松樹的時間序列光譜變化進行了分析，得出近紅外最大一階導數(FD)值、紅邊FD值和藍色邊緣FD值之和的比值可作為松材線蟲檢測的有效因素。張衡等[16]分析了波長593 nm處光譜反射率的一階導數光譜特征，結合葉綠素質量分數，判斷馬尾松是否感染了松萎蔫病，對肉眼可見感病特征前的旱期階段具有良好的監測效果。

研究表明，松樹光譜特征與其健康程度具有較大的相關性，但利用其相關性估測松林感染松材線蟲病等級的研究鮮有報道。本文以重慶永勝林場馬尾松為研究對象，利用70個采樣區的野外高光譜數據，分析松材線蟲危害下馬尾松綠針的光譜特征，構建最大反射率、反射率總和、綠峰反射率高度、紅谷吸引深度及其比率和歸一化值等14個特征參數。利用嶺跡分析篩選最優特征參數，進行嶺回歸建模反演感染松材線蟲病害等級，并對模型估算精度進行驗證，以期為實現松材線蟲害早期監測和防治提供支持。

1 數據獲取

1.1 數據采集

重慶松材線蟲病從6月初開始發病， 7月發生的病樹最多，8月逐漸下降，本研究采樣時間為2017年6—8月，共設置70個采樣區。從各個采樣區內選取受害程度不一的典型松樹針葉進行測量。數據采集過程中使用FieldSepc4型野外光譜分析儀，作用波段為400～1 100 nm，光譜分辨率為1 nm。根據實際需求，測量過程中設定光譜平均次數為5次，暗電流平均次數為10次，白板平均次數為5次，共采集了70個馬尾松植株的有效光譜數據。使用FieldSepc4型野外光譜分析儀配套軟件RS3，優化FieldSpec4型儀器采集的原始計數、輻射亮度/輻射照度等相關數據，計算其光譜反射率。

1.2 數據處理

對70個光譜數據進行歸一化處理，計算式為

(1)

式中ρi——樣本i的歸一化光譜反射率

Ri——樣本i光譜反射率

Rmax、Rmin——樣本最大、最小光譜反射率

經采集樣本專業人員鑒定，按照松樹針葉受害程度，依次劃分為健康、輕度、中度、重度和枯死5個等級，如圖1所示[17]，并量化為蟲害等級1，2，…，5，從而得到樣本光譜數據及蟲害等級數據，如表1所示。

圖1 不同受害程度的松樹針葉Fig.1 Needle leaves suffering from different levels of disease

樣本號坐標經度/(°)緯度/(°)海拔/m光譜反射率ρ400…ρ1100蟲害等級129.8067107.23243860.0359…0.70122229.8069107.23384270.0293…0.65783329.8072107.23464330.0134…0.57683429.8074107.23524420.0364…0.85882???????7029.8064107.23444130.0227…0.95612

2 研究方法

2.1 光譜特征參數構建

光譜特征是物質在電磁波的作用下，在特定波長位置形成的反映物質成分和結構信息的光譜吸收和反射特征[18]。由于受松材線蟲害侵蝕，馬尾松針葉色素、含水率以及細胞構造發生劇烈變化，使得其光譜反射率曲線產生明顯的變化特征[19]。圖2為健康、輕度、中度、重度和枯死樣本所對應光譜反射率曲線[17]。由圖2可知，健康樣本在可見光550 nm有明顯的反射峰，在680 nm處有明顯的吸收谷，在680～780 nm波段(紅邊)反射率急劇上升，在780～1 100 nm波段(近紅外)，反射率相對平穩且有最大值。隨著蟲害程度的加深，樣本反射峰、吸收谷、紅邊反射率陡峭程度、近紅外最大反射率降低。特別當樣本處于枯死狀態時，反射峰和吸收谷明顯消失，紅邊及近紅外區反射率被逐漸拉平。這種光譜上的差異使得應用高光譜遙感技術估測病害程度成為可能。

圖2 健康與感病植株光譜反射率曲線Fig.2 Spectral reflectance curves for healthy and infected pines

根據文獻[13,15,18,20]對波段的劃分，選擇綠光區(490～560 nm)、黃光區(560～590 nm)、紅光區(620～680 nm)、紅邊(680～780 nm)、近紅外(780～1 100 nm)波段。計算上述波段范圍內光譜反射率最大值和反射率總和，500～670 nm波段范圍內綠峰反射高度和560～760 nm波段范圍內紅谷吸收深度、其比值和歸一化值，共14個特征參數，分別表示為Mg、My、Mr、Mre、Minr、Sg、Sy、Sr、Sre、Sinr、GH、RD、DH、D_H，如表2所示。由于從光譜特征上可區分枯死樣本，因此后續研究只針對健康、輕度、中度、重度4種類型。

2.2 嶺回歸估計模型

在估測模型構建過程中，特征參數篩選尤為關鍵。上述14個特征參數中，有些特征參數相關性較強，不適合全部用于建模。為了提高模型估計準確度，需刪除一些相關性較強的特征參數。常用的變量優選方法主要有：基于自適應重加權[21]、多元逐步線性回歸[22]和遺傳算法[23]等。嶺估計是一種改進的最小二乘法，其本質是一種放棄最小二乘法的無偏性與部分精確度，尋求效果稍差但更符合實際情況的回歸過程，其在特征提取方面效果較優，并可

表2 光譜特征參數定義Tab.2 Definition of spectral characteristic parameters

簡化模型和提高模型的魯棒性[24]。因此本研究采用嶺回歸分析優選特征參數，并用于估測模型構建。

多元線性回歸模型可表示為

Y=Xβ+ε

(2)

式中X——特征參數組成的觀測陣

Y——等級量化值組成的向量陣

β——回歸系數ε——隨機誤差

(3)

式中XT——觀測陣X的轉置矩陣

待定參數的嶺估計嶺回歸系數計算公式為

β(k)=(XTX+kIp)-1XTY

(4)

式中k——嶺跡參數

Ip——p階單位矩陣

p——建模樣本數

為了改善方陣XTX特征根趨于0的情況，式(4)中引入嶺跡參數k，使所求估計系數盡可能達到最優值。將估計系數β作為k的函數，當k在[0,1]范圍內變化時，以k為橫坐標，β(k)為縱坐標，其在平面直角坐標中所對應的圖形即為嶺跡。嶺跡可直觀地反映各特征參數對蟲害等級估測的作用和相互關系，借助嶺跡分析，可有效地選取影響蟲害等級的主要特征參數。

基于嶺跡分析篩選特征參數的原則有[24]：①隨著嶺跡參數k值增加，嶺回歸系數很不穩定的特征參數刪除。②隨著嶺跡參數k值增加，嶺回歸系數迅速趨于零的特征參數刪除。③其嶺回歸系數雖然穩定但其絕對值比較小的特征參數刪除。根據嶺跡分析，當剩余特征參數的嶺跡趨于基本穩定時，取嶺跡參數k的值，根據嶺跡參數k計算得特征參數的待定系數，進而可得蟲害等級估測模型(式(2))。

3 結果與分析

對于采樣的70個光譜數據，根據光譜曲線去除枯死樣本，剩余56個光譜數據。為了使建模數據集和測試數據集可以充分反映研究區松材線蟲害程度，將56個樣本數據按光譜反射率平均值從大到小進行排序，等間隔抽取35個作為建模數據集，21個作為測試數據集。對于用于建模的35個樣本提取其特征參數與蟲害等級，對特征參數組成的自變量X，計算其方陣XTX，然后對其進行主成分分析，分析其是否存在線性關系，計算XTX的特征值及特征向量，結果如表3所示。其中9個特征參數特征值較小，表明他們之間線性相關性較大，因此需剔除。但需要剔除哪些特征參數，還需借助嶺跡圖來進行判斷。將特征參數組成的自變量X進行中心標準化，利用Matlab進行嶺回歸分析，得到各特征參數嶺跡如圖3所示。

表3 方陣XTX特征值Tab.3 Eigenvalue of square matrix XTX

基于嶺跡分析篩選特征參數原理，最符合條件的特征參數有4個，分別為X4(Mre)、X5(Minr)、X9(Sre)、X10(Sinr)，其嶺跡變化平穩且不趨于0，其次是X12(RD)。分別計算健康、輕度、中度、重度4種樣本(圖2)特征參數X4、X5、X9、X10、X12，如表4所示。由表4可知，隨著蟲害程度的加深，特征參數X4、X5、X9、X10逐漸減小，與蟲害程度呈負相關；特征參數X12總體呈上升趨勢，與蟲害程度呈正相關，與嶺跡分析結果一致。由此可見，紅谷、紅邊和近紅外區域是反映蟲害程度的敏感區域，因此利用高光譜特征估算蟲害程度是可行的。

圖3 特征參數嶺跡圖Fig.3 Ridge traces for characteristic parameters

表4 不同蟲害程度樣本特征參數Tab.4 Values of characteristic parameters for samples in different healthy levels

選取特征參數X4、X5、X9、X10、X12，由圖3可知，當k=0.2時，上述5個特征參數趨于穩定。因此取k=0.2進行嶺回歸建模，得到基于5個特征參數的蟲害等級估測模型為

Y=4.554 4-0.261 2X4-0.011 8X5-0.069 3X9- 0.002 7X10+1.791 6X12

(5)

使用決定系數R2及均方根誤差RMSE評價估測模型預測能力。R2用于判定模型的穩定程度，其值越接近1則估測模型的穩定性越好；RMSE用于表征模型的準確度，其值越小則估測模型精度越高。

為了比較基于特征參數X4、X5、X9、X10和X12的嶺回歸估測模型穩定性和準確度，實驗中分別對14個特征參數進行二進制編碼G=[X1X2…X14]。當特征參數Xi用于建模時，將其編碼為1，否則為0。當每個特征參數分別編碼為0、1時， 共構建214-1個嶺回歸估測模型，其均方根誤差(RMSE)統計如圖4所示。由圖4可知，基于X4、X5、X9、X10、X12構建的估測模型R2=0.868 6，RMSE為0.273 5；基于14個特征參數的多元嶺回歸估測模型R2=0.859 5，RMSE為0.288 0，且其余估測模型的RMSE均大于0.273 5，表明嶺跡分析能夠優選特征參數。

圖4 不同特征參數對應的嶺回歸估測模型均方根 誤差(RMSE)Fig.4 Root-mean-square error (RMSE) of ridge regression models with different characteristic parameters

將剩余21個樣本數據使用式(5)進行估測，所得到的結果進行驗證，驗證計算式為

(6)

式中Ei——樣本i的估測精度

Yri——樣本i的實測數據

Yei——樣本i的估測數據

其結果如表5及圖5所示。由表5可知，模型Y1決定系數R2為0.829 4，均方根誤差RMSE為0.328 0，平均精度為87.15%，與模型Y2相比其決定系數與精度均較大，而均方根誤差較小。從圖5可知，模型Y1得到的估測數據與實測數據緊密分布在y=x周圍，估測值與實際值相差不大。因此估測模型嶺回歸估測模型Y1具有較好的估測能力。

表5 估測模型驗證結果Tab.5 Validation results of estimation model

圖5 嶺回歸估測模型Y1驗證Fig.5 Validation of ridge regression estimation model Y1

4 結論

(1)健康樣本在綠光區(490～560 nm)和紅光區(620～680 nm)有明顯的反射峰和吸收谷，在紅邊內(680～760 nm)反射率急劇上升，在近紅外區域內(780～1 100 nm)反射率相對平穩且有最大值。

隨著蟲害程度的加深，樣本反射峰、吸收谷、紅邊反射率陡峭程度、近紅外最大反射率降低。

(2)根據嶺跡圖可知，紅邊反射率最大值X4、近紅外區反射率最大值X5，紅邊反射率總和X9、近紅外區反射率總和X10和紅谷吸收深度X12嶺跡曲線變化穩定且不趨于零，符合嶺跡分析篩選標準，因此選取這個5個特征參數用于嶺回歸建模。

(3)當k=0.2時，特征參數X4、X5、X9、X10和X12嶺跡趨于穩定。計算嶺回歸系數β(k=0.2)，進行嶺回歸建模，得松材線蟲嶺回歸估測模型，其決定系數為0.868 6，均方根誤差為0.273 5。

(4)使用測試數據集對嶺回歸估測模型精度進行驗證，并將結果與基于14個光譜特征的嶺回歸估測模型結果進行對比，可得其決定系數為0.829 4，均方根誤差為0.328 0，平均精度為87.15%，均優于基于14個特征參數的嶺回歸估測模型。