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

張素蘭 黃金龍 秦 林 李宏群

(1.長(zhǎng)江師范學(xué)院大數(shù)據(jù)與智能工程學(xué)院， 重慶 408100； 2.昆士蘭大學(xué)園藝科學(xué)中心， 布里斯班 4072；3.長(zhǎng)江師范學(xué)院電子信息工程學(xué)院， 重慶 408100； 4.長(zhǎng)江師范學(xué)院三峽庫(kù)區(qū)環(huán)境監(jiān)測(cè)與災(zāi)害防治工程研究中心， 重慶 408100)

0 引言

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

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

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

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

1 數(shù)據(jù)獲取

1.1 數(shù)據(jù)采集

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

1.2 數(shù)據(jù)處理

對(duì)70個(gè)光譜數(shù)據(jù)進(jìn)行歸一化處理，計(jì)算式為

(1)

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

Ri——樣本i光譜反射率

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

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

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

樣本號(hào)坐標(biāo)經(jīng)度/(°)緯度/(°)海拔/m光譜反射率ρ400…ρ1100蟲害等級(jí)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 光譜特征參數(shù)構(gòu)建

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

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

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

2.2 嶺回歸估計(jì)模型

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

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

簡(jiǎn)化模型和提高模型的魯棒性[24]。因此本研究采用嶺回歸分析優(yōu)選特征參數(shù)，并用于估測(cè)模型構(gòu)建。

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

Y=Xβ+ε

(2)

式中X——特征參數(shù)組成的觀測(cè)陣

Y——等級(jí)量化值組成的向量陣

β——回歸系數(shù)ε——隨機(jī)誤差

(3)

式中XT——觀測(cè)陣X的轉(zhuǎn)置矩陣

待定參數(shù)的嶺估計(jì)嶺回歸系數(shù)計(jì)算公式為

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

(4)

式中k——嶺跡參數(shù)

Ip——p階單位矩陣

p——建模樣本數(shù)

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

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

3 結(jié)果與分析

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

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

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

圖3 特征參數(shù)嶺跡圖Fig.3 Ridge traces for characteristic parameters

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

選取特征參數(shù)X4、X5、X9、X10、X12，由圖3可知，當(dāng)k=0.2時(shí)，上述5個(gè)特征參數(shù)趨于穩(wěn)定。因此取k=0.2進(jìn)行嶺回歸建模，得到基于5個(gè)特征參數(shù)的蟲害等級(jí)估測(cè)模型為

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

(5)

使用決定系數(shù)R2及均方根誤差RMSE評(píng)價(jià)估測(cè)模型預(yù)測(cè)能力。R2用于判定模型的穩(wěn)定程度，其值越接近1則估測(cè)模型的穩(wěn)定性越好；RMSE用于表征模型的準(zhǔn)確度，其值越小則估測(cè)模型精度越高。

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

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

將剩余21個(gè)樣本數(shù)據(jù)使用式(5)進(jìn)行估測(cè)，所得到的結(jié)果進(jìn)行驗(yàn)證，驗(yàn)證計(jì)算式為

(6)

式中Ei——樣本i的估測(cè)精度

Yri——樣本i的實(shí)測(cè)數(shù)據(jù)

Yei——樣本i的估測(cè)數(shù)據(jù)

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

表5 估測(cè)模型驗(yàn)證結(jié)果Tab.5 Validation results of estimation model

圖5 嶺回歸估測(cè)模型Y1驗(yàn)證Fig.5 Validation of ridge regression estimation model Y1

4 結(jié)論

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

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

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

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

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