基于高光譜數(shù)據(jù)與針葉電阻率的松材線蟲病早期診斷研究

劉 飛,談家金*,謝谷艾,蘇 峻,葉建仁

1. 南京林業(yè)大學(xué)林學(xué)院,南方現(xiàn)代林業(yè)協(xié)同創(chuàng)新中心,江蘇 南京 210037 2. 江西省林業(yè)科學(xué)院林業(yè)有害生物防治研究所,江西 南昌 330013 3. 南京林業(yè)大學(xué)理學(xué)院,江蘇 南京 210037

引 言

松材線蟲(Bursaphelenchusxylophilus)是目前對(duì)我國(guó)森林危害最嚴(yán)重的林業(yè)檢疫性有害生物,由其引起的松材線蟲病是一種松樹毀滅性病害[1-3]。我國(guó)自1982年發(fā)現(xiàn)該病以來(lái),病情經(jīng)過(guò)41年不斷發(fā)展,目前全國(guó)共有19個(gè)省(自治區(qū)、直轄市)701個(gè)縣級(jí)行政區(qū)發(fā)生疫情,松材線蟲病的準(zhǔn)確診斷對(duì)該病的有效防控至關(guān)重要。該病現(xiàn)行的診斷技術(shù)包括林間癥狀診斷法、病原線蟲鑒定法、流膠法等,這些方法受到條件或技術(shù)的限制而不完善;為了能夠有效發(fā)現(xiàn)和控制松材線蟲病,其爆發(fā)前的早期診斷更是成為當(dāng)前研究的熱點(diǎn)[4]。如高瑞賀等發(fā)現(xiàn)健康馬尾松感染松材線蟲病后,其針葉的蒸騰速率和凈光合速率變化顯著,可作為松材線蟲病早期診斷的指標(biāo)[5];Gaspar等在接種松材線蟲后,發(fā)現(xiàn)松樹3-蒈烯揮發(fā)量相比健康松樹顯著提高,可用于松材線蟲病早期診斷的指標(biāo)[6]。

大量研究表明,植物在受生物脅迫或非生物脅迫過(guò)程中,其自身色素、含水量等生理生化指標(biāo)會(huì)發(fā)生一定的變化[7-8]。而植物對(duì)電磁輻射的吸收與反射特性會(huì)隨著這些生理生化參數(shù)的改變而改變,在遙感上表現(xiàn)為不同的光譜反射率。因而遙感在研究植物脅迫方面應(yīng)用廣泛,可見光區(qū)和近紅外區(qū)是指示植物色素含量變化和健康情況的關(guān)鍵波段[9-11]。健康馬尾松在感染松材線蟲病后,其葉片內(nèi)部結(jié)構(gòu)會(huì)被破壞,光合色素和水分含量也會(huì)發(fā)生變化[12];伴隨病害不斷的發(fā)展,馬尾松的光譜反射率呈現(xiàn)出較大變化[13]。松材線蟲病在病情發(fā)展的早期階段,感病松樹自身不會(huì)發(fā)生明顯癥狀,且葉綠素和水分含量變化不明顯,但當(dāng)肉眼可以觀察到針葉葉色明顯變色時(shí),感病松樹會(huì)在短期內(nèi)萎蔫死亡。目前研究中,無(wú)論是基于光譜本身的直接診斷方法,或是基于生理生化組分的間接診斷方法在病害的早期階段效果并不理想[14-15]。

為了實(shí)現(xiàn)對(duì)感病植株的早期診斷,可以考慮尋找一種基于光譜的新方法用于判斷植株是否感病,這將對(duì)松材線蟲病的防治具有重要的實(shí)踐意義與應(yīng)用價(jià)值。馬尾松感染松材線蟲病后除生理生化指標(biāo)發(fā)生變化外,其生物電也必然會(huì)發(fā)生一定的變化[16]。本研究擬從電學(xué)角度出發(fā),在研究感病松樹針葉電阻率變化的基礎(chǔ)上,將基于馬尾松針葉反射光譜數(shù)據(jù)建立針葉電阻率的預(yù)測(cè)模型,以實(shí)現(xiàn)在感病后的早期階段對(duì)馬尾松松材線蟲病的快速判別,為利用無(wú)人機(jī)、衛(wèi)星遙感對(duì)松林進(jìn)行松材線蟲病監(jiān)測(cè)提供理論依據(jù)。

1 實(shí)驗(yàn)部分

1.1 松材線蟲接種、病情觀察及樣品采集

實(shí)驗(yàn)地設(shè)在安徽省休寧縣商山鎮(zhèn)五方水庫(kù)的一處馬尾松林內(nèi)(118°21′63″E,29°69′23″N;156.7 m)。2021年8月17號(hào)采用注干法對(duì)8～9年生馬尾松接種松材線蟲,實(shí)驗(yàn)共設(shè)置處理和對(duì)照兩組,處理組和對(duì)照組的馬尾松分別設(shè)置6棵和4棵。在松樹主干距離地面1.5 m處,使用直徑6 mm的鉆頭(鉆頭與主干夾角約45°)鉆入木質(zhì)部4 cm,在鉆孔內(nèi)放入適量的無(wú)菌棉球;使用1 mL的槍頭(槍頭前端已被修剪為注射器針頭狀,增加接觸面積)插入鉆孔,用微量進(jìn)液器注入線蟲懸浮液,處理組接種線蟲量為30 000條·株-1(濃度30 000條·mL-1),對(duì)照組接等量無(wú)菌水。接種后每4～7天對(duì)馬尾松發(fā)病狀況進(jìn)行觀察并記錄感病指數(shù)。馬尾松病情分級(jí)與感病指數(shù)按談家金等的方法[17],略加修改后計(jì)算:0級(jí)正常,針葉綠色,代表值為0;1級(jí)為1/2以下針葉褪綠、1/4以下葉發(fā)黃,代表值為1;2級(jí)為1/2以上葉褪綠、1/4～3/4葉發(fā)黃,代表值為2;3級(jí)為3/4以上葉發(fā)黃、1/2以下葉變紅,代表值為3;4級(jí)為1/2以上葉變紅、植株瀕死或死亡,代表值為4。計(jì)算公式:感病指數(shù)=Σ(各病級(jí)株數(shù)×病級(jí)代表值)×100/(總株數(shù)×最高病級(jí)的代表值)。接種當(dāng)天以及接種后7、11、15、19和26 d,采集外觀表現(xiàn)健康的針葉進(jìn)行光譜和電阻的測(cè)定。樣品采集時(shí)間設(shè)置在上午10:00—12:00,選擇晴朗無(wú)風(fēng)無(wú)云或少云天氣,從馬尾松向陽(yáng)面3個(gè)位置(即冠層上、中、下3個(gè)位置各取一個(gè)點(diǎn))采集當(dāng)年生成熟針葉,將采集好的樣本針葉用密封袋封口,放置在冰盒內(nèi)帶回室內(nèi)進(jìn)行光譜和電阻的測(cè)定。

1.2 針葉光譜反射率測(cè)定

采用Ocean Optics USB2000+光譜儀對(duì)針葉光譜反射率進(jìn)行測(cè)定,該光譜儀由美國(guó)Ocean Insight公司生產(chǎn),其作用波段范圍為350～1 000 nm,光譜分辨率為1.5 nm,光譜采樣間隔為0.44 nm。光譜儀探頭視場(chǎng)角25.4°,探頭距離樣本表面距離8 cm,采用10 W鹵素?zé)魹楣庠?方位角為45°(光源與樣本表面的夾角)。每次測(cè)定前,使用白色面板進(jìn)行光譜校正。測(cè)定儀器各配件如圖1所示,開關(guān)電源(a)、探頭(b)、樣本槽(c)和光源(d)均放置在內(nèi)部密閉黑色的箱子內(nèi),如圖1(A)所示。測(cè)定過(guò)程中設(shè)定光譜平均次數(shù)為10次,積分時(shí)間100 ms。對(duì)針葉樣本進(jìn)行測(cè)定時(shí),使用黑色雙面膠將針葉中間部分整齊平鋪在黑色樣本板上,并剪去多余部分,針葉之間不能有縫隙,如圖1(B)所示。將冠層上、中、下3個(gè)位置光譜反射率平均值作為該植株的光譜反射率。

圖1 針葉光譜反射率測(cè)定儀器各配件(A)黑色箱子內(nèi)部各組件:(a)開關(guān)電源;(b)探頭;(c)樣本槽;(d)光源;(B)樣品板與待測(cè)針葉樣本Fig.1 Accessories of needle spectral reflectance measurement instrument(A) Internal components of the black box:(a) Switching power supply;(b) The probe;(c) The sample slot;(d) The light source;(B) Sample plates and sample of needles to be tested

1.3 針葉電阻率測(cè)定

采用M4070 LCR測(cè)試儀測(cè)定針葉電阻,該儀器是由中國(guó)晶研儀器科技有限公司生產(chǎn),其量程范圍為0.00 Ω～10.000 MΩ,基本測(cè)量準(zhǔn)確度0.5%。使用3D打印機(jī)打印黑色薄板,將兩個(gè)電極插入對(duì)應(yīng)位置,在半球形凹槽內(nèi)放入導(dǎo)電膏(濟(jì)寧高新區(qū)金諾特醫(yī)用凝膠廠生產(chǎn)),方便針葉電阻的測(cè)量,如圖2所示。將針葉橫截面近似看成半橢圓形,剪取針葉中部4 cm,使用游標(biāo)卡尺(分辨率0.01 mm)對(duì)針葉的寬度和厚度進(jìn)行測(cè)量。將針葉與M4070 LCR測(cè)試儀連接后3～5 s的讀數(shù)作為其電阻值,從而計(jì)算出其電阻率。冠層向陽(yáng)面上、中、下3個(gè)位置各取5根針葉進(jìn)行電阻的測(cè)定,將5根針葉電阻率平均值(電阻率異常值的樣本除外)作為對(duì)應(yīng)冠層位置針葉電阻率,將上、中、下3個(gè)位置電阻率平均值作為該植株的針葉電阻率。

圖2 針葉電阻測(cè)定過(guò)程Fig.2 Needle resistance measurement process

1.4 數(shù)據(jù)分析

1.4.1 數(shù)據(jù)預(yù)處理

地面光譜與針葉電阻同步采樣共6次,將獲取的數(shù)據(jù)去除明顯異常值后,共獲得52組數(shù)據(jù)。其中,處理組含34組數(shù)據(jù),對(duì)照組含18組數(shù)據(jù)。使用Excel數(shù)據(jù)分析工具對(duì)光譜數(shù)據(jù)和電阻數(shù)據(jù)處理匯總;使用SPSS軟件對(duì)電阻率數(shù)據(jù)進(jìn)行分析,并在Origin 2018中繪圖。

在MATLAB R2021b中調(diào)用Savitzky-Golay函數(shù)對(duì)光譜數(shù)據(jù)進(jìn)行平滑處理,既可以消除數(shù)據(jù)中的噪聲部分,又能將數(shù)據(jù)中重要的信息保留,提高了數(shù)據(jù)的平滑性。將平滑后的光譜定義成原始光譜(OR)。光譜數(shù)據(jù)自身具有復(fù)雜性和高維數(shù)特點(diǎn),為了提取更多有用的信息,需要對(duì)光譜數(shù)據(jù)進(jìn)行一些復(fù)雜的數(shù)學(xué)變換,為預(yù)測(cè)模型的創(chuàng)建提供充足的信息源。目前通常的光譜變換包括:一階導(dǎo)數(shù)(FD)、二階導(dǎo)數(shù)(SD)、對(duì)數(shù)變換(LOG)、倒數(shù)變換(1/R)和連續(xù)統(tǒng)去除法(CR)。

1.4.2 特征波段選擇

高光譜數(shù)據(jù)具有高維、復(fù)雜的特點(diǎn),為了篩選出與針葉電阻率相關(guān)性高的波段,采用隨機(jī)森林進(jìn)行特征波段篩選。隨機(jī)森林(random forest,RF)是一種非線性模型,以決策樹作為基學(xué)習(xí)器,使用裝袋(Bagging)原理來(lái)處理訓(xùn)練數(shù)據(jù)集[18]。由于其計(jì)算過(guò)程相對(duì)簡(jiǎn)單,對(duì)樣本和變量的選擇過(guò)程是隨機(jī)的,因而能夠在多種領(lǐng)域應(yīng)用并且均可以得到很好的分類效果[19]。這些特點(diǎn)使得RF特別適用于光譜數(shù)據(jù)的研究。本次實(shí)驗(yàn)使用RF進(jìn)行特征波段的篩選,RF中樹的數(shù)量nTrees設(shè)置300,輸入的特征數(shù)量10。

1.4.3 最小二乘支持向量機(jī)回歸模型

最小二乘支持向量機(jī)(least squares support vector machine,LSSVM)算法是基于支持向量機(jī)算法的一種改進(jìn),極大提高了求解問(wèn)題的速度和收斂精度[20],且不再受到樣本數(shù)量的影響,具有很強(qiáng)的泛化能力。

采用LSSVM算法對(duì)針葉光譜反射率與針葉電阻率進(jìn)行回歸分析,從而對(duì)馬尾松針葉電阻率進(jìn)行反演。以光譜變換下的特征波段為自變量,針葉電阻率為因變量,基于LSSVM算法構(gòu)建馬尾松針葉電阻率預(yù)測(cè)模型。為了使得模型測(cè)試的結(jié)果更客觀,80%的樣本選作建模集,20%的樣本用作驗(yàn)證集。主要用以下三個(gè)指標(biāo)評(píng)估馬尾松針葉電阻率預(yù)測(cè)模型的精度:決定系數(shù)(coefficient of determination,R2)、均方根誤差(root mean square error,RMSE)、平均絕對(duì)誤差(mean absolute error,MAE)。當(dāng)預(yù)測(cè)模型的R2趨近1、RMSE和MAE越小時(shí),表明其精度越高。

2 結(jié)果與討論

2.1 接種后馬尾松病情發(fā)展

接種后11 d,發(fā)現(xiàn)有1棵接種線蟲的馬尾松1個(gè)2級(jí)側(cè)枝針葉開始褪綠,此時(shí),病情指數(shù)為4;接種后15 d,該側(cè)枝針葉變?yōu)辄S褐色,此時(shí),病情指數(shù)為4;接種后19 d,該側(cè)枝針葉變成紅褐色,針葉已枯死,此時(shí),已有3棵接種線蟲的馬尾松針葉變色,病情指數(shù)為12.5。接種后26 d,已有4棵接種線蟲的馬尾松極少數(shù)針葉變色,病情指數(shù)為16.7。接種后67 d,接種線蟲的6棵馬尾松少數(shù)針葉均已變色,病情指數(shù)為25。

2.2 接種后馬尾松針葉電阻率的變化

接種實(shí)驗(yàn)結(jié)果顯示,馬尾松顯癥早期,針葉電阻率顯著升高。由圖3可知,接種后15 d,此時(shí)僅1棵接蟲馬尾松的1個(gè)側(cè)枝針葉變色,接蟲株和對(duì)照株的電阻率開始達(dá)到極顯著差異,這種差異一直維持到接種后26 d(p<0.01)。整個(gè)實(shí)驗(yàn)期間,相比接蟲株的針葉電阻率,對(duì)照株的針葉電阻率波動(dòng)較小。因此,從電阻率的大小可以初步判斷當(dāng)馬尾松針葉電阻率大于323 Ω·m時(shí),該松樹已經(jīng)感病。

圖3 接種松材線蟲后馬尾松針葉電阻率的變化**表示在0.01水平差異顯著Fig.3 Changes in needle resistivity of Pinus massoniana inoculated with B. xylophilus** Shows significant correlation difference at 0.01 level

2.3 接種后馬尾松針葉光譜數(shù)據(jù)特征分析

由于部分波段區(qū)域出現(xiàn)大量噪聲,主要對(duì)490～980 nm范圍內(nèi)的馬尾松針葉反射光譜進(jìn)行分析。為了豐富基于光譜建立針葉電阻率模型的信息,使用5種數(shù)學(xué)變換方式對(duì)所有有效光譜數(shù)據(jù)進(jìn)行光譜預(yù)處理,僅對(duì)數(shù)變換沒(méi)有改變光譜曲線的趨勢(shì)[圖4(a—f)]。

圖4 不同預(yù)處理方法的光譜變化(a):OR;(b):FD;(c):SD;(d):LOG;(e):1/R;(f):CRFig.4 Spectral changes of different pretreatment methods(a):OR;(b):FD;(c):SD;(d):LOG;(e):1/R;(f):CR

由表1可知,通過(guò)隨機(jī)森林算法對(duì)各預(yù)處理方法進(jìn)行特征波段篩選后發(fā)現(xiàn)存在較大不同。原始光譜(OR)的特征波段集中在近紅外區(qū)域,一階導(dǎo)數(shù)和連續(xù)統(tǒng)去除法光譜的特征波段集中分布在綠光和近紅外區(qū)域,對(duì)數(shù)變換和倒數(shù)變換篩選出的波段主要分布在紅光和近紅外區(qū)域。綜合比較發(fā)現(xiàn),特征波段主要分布在近紅外區(qū)域,說(shuō)明病害發(fā)生后的早期針葉內(nèi)細(xì)胞結(jié)構(gòu)可能發(fā)生了變化,這與高瑞賀等研究的感染松材線蟲病引起馬尾松針葉光合作用相關(guān)機(jī)構(gòu)改變的結(jié)果相一致[5]。大量研究表明,在可見光和近紅外區(qū)域的反射光譜在病蟲害監(jiān)測(cè)中具有很高的應(yīng)用價(jià)值,它們通常與病害的發(fā)展過(guò)程密切相關(guān)[9-10]。

表1 基于隨機(jī)森林算法篩選的特征波段Table 1 Feature bands screening based on random forest algorithm

2.4 馬尾松針葉電阻率預(yù)測(cè)模型分析

基于原始光譜進(jìn)行五種光譜變換,使用隨機(jī)森林對(duì)原始光譜和各預(yù)處理后的光譜數(shù)據(jù)進(jìn)行篩選,從而得到與針葉電阻率相關(guān)特征波段,最終利用LSSVM算法建立的針葉電阻率預(yù)測(cè)模型。該模型的預(yù)測(cè)效果匯總于表2,發(fā)現(xiàn)不同預(yù)處理間預(yù)測(cè)效果的綜合表現(xiàn)具有差異,其中二階導(dǎo)數(shù)模型預(yù)測(cè)效果綜合表現(xiàn)最好,建模集與驗(yàn)證集的平均R2、RMSE和MAE分別為0.848、32.332和7.067。

表2 基于不同光譜變換篩選輸入變量的模型表現(xiàn)Table 2 Performance of model based on the different spectral transformation input variable selection

各預(yù)測(cè)模型間除了倒數(shù)變換和連續(xù)統(tǒng)去除法變換外,其他光譜變換后的預(yù)測(cè)模型比原始光譜預(yù)測(cè)模型驗(yàn)證集的R2均有所提升,RMS和MAE均有所降低。一階導(dǎo)數(shù)預(yù)測(cè)模型的改善效果較小,建模集與驗(yàn)證集的平均R2僅提升1.5%,但RMSE和MAE卻分別升高28%和1.7%。而二階導(dǎo)數(shù)模型預(yù)測(cè)效果提升較為明顯,R2提升了4.2%,RMSE和MAE分別降低2.5%和18.9%。對(duì)數(shù)模型預(yù)測(cè)效果也有所提升,R2提升了3.1%,但RMSE和MAE分別升高了57.3%和11.3%。綜合比較,基于二階導(dǎo)數(shù)選擇的特征變量在最小二乘支持向量機(jī)算法預(yù)測(cè)模型中綜合表現(xiàn)最優(yōu),特征波段為594.986、646.107、646.451、782.896、784.841、839.164、863.890、902.021、947.901、962.315;其次為L(zhǎng)OG、FD、1/R和CR,這表明對(duì)光譜預(yù)處理的必要性,處理后的光譜數(shù)據(jù)增強(qiáng)了樣本的光譜反應(yīng)特征,進(jìn)一步提升反演模型的精度。

圖5為基于二階導(dǎo)數(shù)(SD)變換的LSSVM建模方法反演結(jié)果散點(diǎn)圖,從圖5中可以清晰看出,該模型的實(shí)測(cè)值樣點(diǎn)和預(yù)測(cè)值樣點(diǎn)多數(shù)分布在y=x線附近,能夠較好地對(duì)馬尾松針葉電阻率進(jìn)行估測(cè)。因此,認(rèn)為SD-RF-LSSVM可作為馬尾松針葉電阻率大小的最優(yōu)估測(cè)模型。相較于傳統(tǒng)的林間癥狀診斷法、病原線蟲鑒定法以及流膠法[4],本文提出的基于SD-RF-LSSVM模型的光譜檢測(cè)方法不僅可以實(shí)現(xiàn)對(duì)該病害的快速無(wú)損診斷;此外,該方法能夠在病樹極少數(shù)針葉變色時(shí)就完成對(duì)病害的監(jiān)測(cè),相比之下,傳統(tǒng)的遙感監(jiān)測(cè)方法需要等到病樹針葉大部分已經(jīng)明顯變色才能發(fā)現(xiàn)病害[16]。因此,該方法可以在病樹出現(xiàn)極少數(shù)針葉變色時(shí)進(jìn)行早期預(yù)警和監(jiān)測(cè),及時(shí)發(fā)現(xiàn)病害的發(fā)生和變化情況,從而更好地制定病害防治措施。

圖5 基于二階導(dǎo)數(shù)數(shù)據(jù)變換的最小二乘支持向量機(jī)模型檢驗(yàn)結(jié)果Fig.5 Test results of LSSVM model based on the SD data transformation

3 結(jié) 論

多年來(lái),我國(guó)在松材線蟲病的致病機(jī)理及防治方法上取得許多重要進(jìn)展,然而在該病的早期診斷方面一直未有突破,而利用遙感技術(shù)對(duì)該病的早期監(jiān)測(cè)和預(yù)警長(zhǎng)期處在探索階段。本研究表明,使用針葉反射光譜反演針葉電阻率是可行的,且SD-RF-LSSVM建立的預(yù)測(cè)模型精度最高,可用于針葉電阻率的快速估測(cè),為實(shí)現(xiàn)基于遙感的松材線蟲病早期診斷與監(jiān)測(cè)提供了思路與方法。

實(shí)驗(yàn)僅對(duì)針葉葉片層面的光譜特征進(jìn)行分析,雖可為將來(lái)冠層尺度研究提供參考,但冠層尺度的松樹冠層監(jiān)測(cè)更復(fù)雜,光譜識(shí)別與分析受到冠層形狀、土壤、植被等因素影響。需要在冠層尺度比較健康和感病馬尾松的光譜特征差異,從而為使用遙感技術(shù)實(shí)現(xiàn)病害的快速有效識(shí)別提供方法。此外樣本數(shù)量較少,模型需要進(jìn)一步檢驗(yàn)和驗(yàn)證。同時(shí),環(huán)境條件對(duì)結(jié)果的影響需要進(jìn)一步研究,增強(qiáng)模型的穩(wěn)定性和適用性是今后研究的重點(diǎn)。