剛竹毒蛾脅迫下毛竹葉綠素與葉光譜特征的關(guān)系及其變化

胡新宇, 許章華, 3, 5, 6*, 黃旭影, 8， 張藝偉， 陳秋霞，王 琳， 劉 輝， 劉智才

1. 福州大學(xué)環(huán)境與安全工程學(xué)院， 福建 福州 350108 2. 福建省資源環(huán)境監(jiān)測與可持續(xù)經(jīng)營利用重點(diǎn)實(shí)驗(yàn)室， 福建 三明 365004 3. 空間數(shù)據(jù)挖掘與信息共享教育部重點(diǎn)實(shí)驗(yàn)室, 福建 福州 350108 4. 福州中谷海創(chuàng)科技發(fā)展有限公司， 福建 福州 350108 5. 福州大學(xué)先進(jìn)制造學(xué)院， 福建 泉州 362200 6. 福州大學(xué)信息與通信工程博士后科研流動站， 福建 福州 350108 7. 福建農(nóng)林大學(xué)公共管理學(xué)院， 福建 福州 350002 8. 南京大學(xué)國際地球系統(tǒng)科學(xué)研究所， 江蘇 南京 210023

引 言

據(jù)估算， 全球竹林總面積為30 538.35×103hm2[1]。 全國第九次森林資源清查顯示， 我國有竹林面積逾640萬hm2， 竹子種類和竹林面積穩(wěn)居全球首位。 毛竹具有根系密集、 生長周期短、 經(jīng)濟(jì)與生態(tài)價(jià)值高等特點(diǎn)， 是我國重要經(jīng)濟(jì)竹種， 竹產(chǎn)業(yè)年產(chǎn)值已超2 000億元， 成為許多地區(qū)踐行習(xí)近平總書記“綠水青山就是金山銀山”理念、 助力鄉(xiāng)村振興的重要資源保障。 然而， 蟲害一直威脅著竹林健康， 并呈現(xiàn)種類多、 致?lián)p嚴(yán)重的態(tài)勢， 導(dǎo)致林分質(zhì)量下降， 竹林價(jià)值受損。 其中， 剛竹毒蛾是我國竹林最主要的一種食葉性害蟲， 具有群發(fā)性、 周期性、 危害極嚴(yán)重等特征， 廣泛分布于福建、 江西、 浙江等多個(gè)省區(qū)。 葉綠素含量是綠色植物生長狀態(tài)的重要評價(jià)指標(biāo)， 是評價(jià)植物光合作用能力、 監(jiān)測生長環(huán)境優(yōu)劣、 判斷植物生理營養(yǎng)是否充足的重要依據(jù)。 葉綠素含量監(jiān)測主要可采用分光光度計(jì)、 熒光或活體葉綠素儀等進(jìn)行測定[2]， 但這些方法效率偏低并可能影響植物的正常生長， 同時(shí)也難以全面掌握林分的整體狀態(tài)。 高光譜技術(shù)為葉綠素的定量化監(jiān)測提供了一種簡便有效且破壞度較小的途徑， 相比于多光譜， 其對光譜信息的細(xì)化也為深入探索植被葉綠素變化規(guī)律提供了重要基礎(chǔ)。 Sonobe等[3]發(fā)現(xiàn)， 采用預(yù)處理技術(shù)與機(jī)器學(xué)習(xí)結(jié)合的方式對兩種山葵葉綠素含量估測獲得較好效果； 韓浩坤等[4]在RVI， PSNDb， GNDVI750的基礎(chǔ)上建立了糜子冠層葉綠素含量的高光譜反射率模型； Zhang等[5]發(fā)現(xiàn)反射率一階微分值構(gòu)建的多元回歸方程以及修正的綠色歸一化植被指數(shù)對雷竹葉片葉綠素的擬合效果較好。

蟲害侵染是植被葉綠素等理化參數(shù)發(fā)生改變的重要誘因， 反之葉綠素也可作為植被健康狀況的主要判定依據(jù)之一。 在農(nóng)林病蟲害領(lǐng)域， Souza等[6]等通過衡量葉綠素a在除草劑脅迫下的熒光特征值分析玉米基因是否受草地夜蛾危害的影響； 有報(bào)道采用葉綠素與含水量特征， 將松材線蟲害分為5種等級， 并采用光譜特征參數(shù)構(gòu)建了光譜特征模型； 王景旭等[7]通過實(shí)驗(yàn)發(fā)現(xiàn)隨著切梢小蠹的侵入時(shí)間延長， 葉綠素含量呈現(xiàn)先平緩后急劇的下降趨勢； 有研究選取包含相對葉綠素含量在內(nèi)的理化參數(shù)作為因子， 對剛竹毒蛾危害進(jìn)行了檢測； 竇志國等[8]分析了健康和蟲害蘆葦葉片高光譜反射率與葉綠素含量間的相關(guān)關(guān)系， 建立了葉綠素含量紅邊位置和全波段高光譜反演估算模型； 周曉等[9]以二齡稻縱卷葉螟幼蟲為試驗(yàn)材料， 分析不同投蟲量條件下、 不同生育期水稻冠層的原始光譜、 三邊參數(shù)和SPAD(soil and plant analyzer development)值的變化規(guī)律， 建立了水稻葉綠素相對含量的回歸估算模型。

已有研究多著眼于改進(jìn)算法并利用光譜信息估測植被葉綠素含量， 或者僅探索蟲害發(fā)生前后葉綠素含量的變化。 盡管上述研究已較為豐富， 關(guān)于蟲害脅迫下葉光譜信息與葉綠素含量的變化機(jī)理卻鮮有涉及。 剛竹毒蛾作為威脅毛竹林健康、 制約竹產(chǎn)業(yè)高質(zhì)量與可持續(xù)發(fā)展的主要因素， 探索其蟲害發(fā)生發(fā)展與葉綠素、 葉光譜特征關(guān)系具有重要意義。 研究從地面非成像高光譜數(shù)據(jù)入手， 分析并篩選了4種高光譜指標(biāo)， 基于4種回歸模型構(gòu)建不同致?lián)p狀態(tài)下毛竹葉片葉綠素含量估測模型， 分析比較基于不同剛竹毒蛾危害等級下SPAD與毛竹葉光譜特征關(guān)系， 為進(jìn)一步深入探索剛竹毒蛾蟲害響應(yīng)與毛竹葉片理化參數(shù)變化研究奠定基礎(chǔ)。

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

1.1 野外調(diào)查與毛竹葉片樣本采集

試驗(yàn)區(qū)位于福建省南平市順昌縣， 地理坐標(biāo)為117°29′—118°14′E， 26°38′～27°12′N[圖1(a—c)]。 順昌縣為亞熱帶海洋性季風(fēng)氣候， 干濕季明顯， 森林資源豐富， 是我國南方的重點(diǎn)林區(qū)， 有“中國竹子之鄉(xiāng)”的美譽(yù)。 截至2018年底， 全縣林地面積達(dá)16.7萬hm2， 其中竹林面積為4.4萬hm2， 毛竹立竹量1.1億根。 長期以來， 剛竹毒蛾均是威脅該縣毛竹生長的主要食葉型害蟲， 嚴(yán)重制約竹林健康與竹林產(chǎn)業(yè)可持續(xù)發(fā)展。 據(jù)2019年順昌縣林業(yè)有害生物統(tǒng)計(jì)數(shù)據(jù)， 該縣當(dāng)年剛竹毒蛾發(fā)生面積為1.34萬hm2， 約占全縣竹林面積的1/3。

圖1 順昌縣地理區(qū)位圖(a)、 研究區(qū)遙感影像測量點(diǎn)2D影像(b)和3D影像(c)

剛竹毒蛾的年生代數(shù)因地而異， 福建省多為3代。 課題組于2020年7月—9月(剛竹毒蛾第二代)赴試驗(yàn)區(qū)， 以小班為單位開展現(xiàn)場踏查。 通過隨機(jī)采樣的方式， 于不同林樣地中采集不同竹齡及健康程度的毛竹葉片樣本。 不同于林分尺度依靠蟲株率、 蟲口數(shù)量/密度的方法， 本文葉片尺度危害等級的確定采用綜合判定法： (1)根據(jù)剛竹毒蛾的危害機(jī)制以及國家林業(yè)局發(fā)布的《林業(yè)有害生物發(fā)生及成災(zāi)標(biāo)準(zhǔn)》， 將單株失葉率(健康： 0%、 輕度危害： 0～25%、 中度危害： 25%～50%、 重度危害： >50%)及蟲口數(shù)量(健康： <10條、 輕度危害： 10～30條、 中度危害： 31～80、 重度危害： >80條)列入危害等級劃分的參考因子。 除此之外， 因毛竹林養(yǎng)分的制造、 積累、 分配、 消耗不平衡的周期變化， 使得毛竹出筍、 成竹年和換葉年交替進(jìn)行， 故毛竹林也分為大小年。 根據(jù)竹林歷年的出筍、 成竹平均數(shù)作為劃分大小年的臨界判別值， 高于或等于該值即為大年， 小于該值即為小年。 處于小年期的毛竹會顯示出落葉、 枯黃等特征， 小年階段毛竹處于換葉期， 營養(yǎng)用于竹鞭生長， 其外部形態(tài)與蟲害葉片類似但受蟲害影響較少， 故若不對其加以區(qū)分則會對毛竹葉片光譜信息收集及葉綠素含量測定工作造成極大干擾。 本研究劃分危害等級的葉片均于大年毛竹林中采集， 并額外采集了部分小年毛竹葉片作為參照； (2)以植物保護(hù)、 森林保護(hù)等學(xué)科背景的高校學(xué)者及長期從事森防檢疫工作的林業(yè)從業(yè)人員為對象， 采用專家咨詢法對危害等級進(jìn)行最終判定。

1.2 毛竹葉綠素與葉光譜數(shù)據(jù)的獲取

課題組于試驗(yàn)區(qū)使用日本Mioolta公司生產(chǎn)的美能達(dá)牌SPAD-502P葉綠素計(jì)對毛竹葉片葉綠素含量快速無損監(jiān)測， 通過測量葉片的吸收率， 得到的數(shù)據(jù)應(yīng)與葉片內(nèi)部結(jié)構(gòu)的葉綠素濃度成正比， 通常情況下葉片葉綠素濃度的相對值即為SPAD值。 每葉片于葉尖、 葉中、 葉基處測得3個(gè)SPAD值， 取其平均值作為該葉片的SPAD值。

葉片光譜數(shù)據(jù)利用ASD Field Spec4 HandHeld光譜儀(Analytical Spectral Device， US)及配套的植物光譜探測器(Unit 16558 plant probe)獲取， 該儀器的波長范圍為350～2 500 nm， 光譜分辨率為1.4 nm(350～1 000 nm)和1.1 nm(1 001～2 500 nm)， 重采樣間隔為1 nm， 共有2 151個(gè)波段。 為確保光譜數(shù)據(jù)的準(zhǔn)確性， 每隔20 min進(jìn)行一次白板校正。 每片葉子在黑色背景下測量其10次正面光譜， 在該儀器配套的光譜處理軟件(View Spec Pro)環(huán)境下剔除異常值后取光譜數(shù)據(jù)的平均值作為樣本實(shí)測光譜數(shù)據(jù)。 根據(jù)研究需要， 采用The Unscrambler X10.4軟件中的Savizky-Golay卷積平滑算法對原始光譜進(jìn)行預(yù)處理， 截選400～2 000 nm光譜數(shù)據(jù)予以統(tǒng)計(jì)。

1.3 毛竹葉光譜特征指標(biāo)計(jì)算與關(guān)系模型建立

1.3.1 特征光譜

高光譜數(shù)據(jù)波段數(shù)較多， 波段信息冗余且復(fù)雜。 為了更有效地識別篩選目標(biāo)地物， 節(jié)約計(jì)算成本， 采用連續(xù)投影算法(successive projection algorithm， SPA)提取特征光譜信息。 該算法是一種前向循環(huán)的特征變量選擇算法， 利用向量的投影分析， 選取含有最低冗余度和最小共線性的有效波長， 對光譜波長進(jìn)行優(yōu)選， 減少建模所需變量個(gè)數(shù)， 改善建模條件。 該算法步驟如下：

(1) 在光譜矩陣中任選一條光譜列向量作為起始向量；

(2) 分別計(jì)算剩余列向量在起始向量的正交平面上的投影向量；

(3) 挑選出最大的投影作為下一次投影的起始向量， 直到挑選變量個(gè)數(shù)達(dá)到最大所需個(gè)數(shù)；

(4) 將提取的所有波長組合進(jìn)行多元線性回歸， 最小的均方根誤差(RMSE)所對應(yīng)組合即為最優(yōu)的波長組合。

經(jīng)過上述步驟篩選出特征波長為469， 702和760 nm， RMSE為0.062 4。

有研究表明， 高光譜衍生指標(biāo)與葉綠素含量之間的估測模型是葉綠素變化檢測研究的有效途徑。 因此， 除原始特征光譜外， 本研究選取植被指數(shù)、 光譜微分、 紅邊參數(shù)3種高光譜衍生指標(biāo)以期更清晰全面地探究毛竹葉片葉綠素與葉光譜特征之間的關(guān)系。

1.3.2 植被指數(shù)

植被指數(shù)類型眾多， 部分研究顯示[10-11]， 葉光譜指數(shù)與SPAD之間的經(jīng)驗(yàn)?zāi)Ｐ涂捎行Ч浪闳~綠素含量。 根據(jù)應(yīng)用尺度的不同， 可以分為葉片尺度和冠層尺度兩類光譜指數(shù)。 葉片尺度上應(yīng)用的植被指數(shù)結(jié)構(gòu)相對簡單， 常用的類型主要有比值型(SR)、 差值型(D)、 歸一化差值型(ND)、 新二重差值型(DDn)改進(jìn)版本等。 除使用原始光譜構(gòu)建光譜指數(shù)外， 也可使用一階微分光譜構(gòu)建光譜指數(shù)(可消除基線漂移和平緩背景干擾的影響， 并且可以放大光譜曲線中的細(xì)微變化)[12]。 所選取20種植被指數(shù)作為特征因子， 既包括基于傳統(tǒng)多光譜數(shù)據(jù)所構(gòu)建的寬帶指數(shù)， 如適用于濃密植被的RVI與描述土壤-植被系統(tǒng)的SAVI等； 亦包括較寬帶指數(shù)更為靈敏的窄帶指數(shù)， 如對葉冠層微小變化敏感的NDVI705與對葉綠素濃度敏感度較高的VOG系列指數(shù)等[13-21]。 本工作選取典型高光譜植被指數(shù)(表1)應(yīng)用于毛竹葉片葉綠素含量的反演研究。

1.3.3 光譜微分

經(jīng)過微分處理的光譜數(shù)據(jù)不僅能消除基線漂移和背景干擾的影響， 還可以放大光譜曲線中的細(xì)微變化， 提供比原始光譜更高分辨率、 更高清晰度的光譜輪廓。 一階微分光譜可更好地表示原始光譜數(shù)據(jù)的變化率和極值點(diǎn)， 二階微分則突出了光譜曲線中切線斜率的變化速度。 計(jì)算公式如下

(1)

(2)

式(1)和式(2)中， FDRλi為波段i和波段i+1之間的一階微分光譜值； SDRλi為波段i和i+1之間的二階微分光譜值；λi為第i波段的波長值；Rλi+2，Rλi+1，Rλi分別代表λi+2，λi+1，λi波段處的原始光譜反射率； Δλ代表波段λi+1到λi之間的波長差值。 下文使用469FDR， 469SDR分別代表469 nm處一階微分光譜值與二階微分光譜值， 其他特征波長處同。

1.3.4 紅邊參數(shù)

紅邊參數(shù)可指示植物的健康狀況， 在病蟲害檢測方面有重要作用， 常見的有葉綠素吸收比值指數(shù)(chlorophyll absorption ratio index， CARI)、 紅邊斜率(red edge slope， RES)、 紅邊面積(red edge area， REA)、 紅邊差值植被指數(shù)(red edge difference vegetation index， REDVI)、 紅邊比值植被指數(shù)(red edge ratio index， RERVI)、 紅邊歸一化差值植被指數(shù)(red edge normalized difference vegetation index， RENDVI)， 計(jì)算公式分別為

b=ρ550-550a

即紅邊覆蓋680～780 nm范圍內(nèi)一階微分光譜的最大值。

表1 葉綠素相關(guān)植被指數(shù)

即紅邊覆蓋680～780 nm范圍內(nèi)一階微分光譜的總和。

REDVI=ρ780-ρ680

1.3.5 關(guān)系模型構(gòu)建原理

采用多算法構(gòu)建比較毛竹葉綠素與葉光譜特征的關(guān)系模型， 探究在傳統(tǒng)回歸模型與機(jī)器學(xué)習(xí)回歸模型下， 所選特征是否可有效模擬葉綠素含量與葉光譜特征關(guān)系。 綜上所述， 采用多元線性回歸、 嶺回歸、 隨機(jī)森林回歸、 XGBoost回歸4種模型對毛竹葉片SPAD進(jìn)行擬合。 通過對比4種模型的優(yōu)劣確定SPAD為最佳估測模型， 基于模型估測結(jié)果分析不同危害程度下毛竹葉片葉綠素含量與葉片光譜特征關(guān)系。

(1)多元線性回歸： 是多元回歸分析中最基礎(chǔ)、 最簡單的一種， 即一個(gè)樣本中有多個(gè)特征的線性回歸問題， 對于一個(gè)有n個(gè)特征的樣本i而言， 其回歸結(jié)果如式(3)

(3)

式(3)中，w被統(tǒng)稱為模型的參數(shù)， 其中w0為截距(intercept)，w1～wn為回歸系數(shù)(regression coeffcient)。

(2)嶺回歸： 又稱為吉洪諾夫正則化(Tikhonov regulation)， 是一種專用于共線性數(shù)據(jù)分析的有偏估計(jì)回歸方法， 是在多元線性回歸的損失函數(shù)上加上了正則項(xiàng)， 通過放棄最小二乘法的無偏性， 以損失部分信息、 降低精度為代價(jià)獲得回歸系數(shù)， 使結(jié)果更符合實(shí)際。 該算法不是為了提升模型表現(xiàn)， 而是為了修復(fù)漏洞所設(shè)計(jì)。

(3)隨機(jī)森林： 是一種典型的Bagging集成算法， 其所有基評估器都是決策樹， 回歸樹所集成的森林即為隨機(jī)森林回歸器(圖2)。 其優(yōu)勢在于對數(shù)據(jù)集的適應(yīng)能力強(qiáng)， 具有很好的抗噪性能和極強(qiáng)的擬合能力但是不會產(chǎn)生過擬合現(xiàn)象[22]。

圖2 隨機(jī)森林結(jié)構(gòu)圖

(4)XGBoost(extreme gradient boosting)： 是在GBDT(gradient boosting decision tree)的基礎(chǔ)上改進(jìn)的一種Boosting模型， 憑借其可擴(kuò)展性強(qiáng)， 運(yùn)行速度快等優(yōu)勢， 在機(jī)器學(xué)習(xí)和數(shù)據(jù)挖掘領(lǐng)域中有著深遠(yuǎn)的影響。 與隨機(jī)森林不同的是， XGBoost的決策樹建立過程是根據(jù)算法中線建立的決策樹情況而定， 而隨機(jī)森林中決策樹的建立是相互獨(dú)立的， 每顆子樹為一個(gè)弱分類器[23]。 具體目標(biāo)函數(shù)如式(4)和式(5)

(4)

(5)

式(4)和式(5)中，T為葉子節(jié)點(diǎn)的個(gè)數(shù)， ‖W‖為葉子節(jié)點(diǎn)向量的模，γ表示節(jié)點(diǎn)切分的難度，λ表示L2正則化系數(shù)。

本實(shí)驗(yàn)選取370片毛竹葉片樣本， 經(jīng)判定各危害等級葉片數(shù)量為無危害： 83、 輕度危害： 71、 中度危害： 70、 重度危害： 89、 小年： 57。 將其中70%作為訓(xùn)練集(樣本容量為259)， 30%作為測試集(樣本容量為111)代入4個(gè)模型之中進(jìn)行評估。

2 結(jié)果與討論

2.1 毛竹葉綠素含量變化分析

剛竹毒蛾作為典型食葉性害蟲威脅毛竹林健康， 主要通過侵蝕葉片影響毛竹葉綠素含量， 使植物葉片組織發(fā)生改變， 光合作用減少進(jìn)而改變毛竹葉片光譜吸收特征。 圖3顯示了全體葉片樣本、 不同危害等級及小年情況下毛竹葉片樣本SPAD的變化情況。 當(dāng)剛竹毒蛾蟲害初發(fā)時(shí)， 毛竹葉中的各種營養(yǎng)成分含量尚未出現(xiàn)明顯變化； 隨著失葉程度的不斷加重， 竹體的養(yǎng)分循環(huán)系統(tǒng)被完全破壞， 由此導(dǎo)致其中的養(yǎng)分元素含量產(chǎn)生顯著變化， 進(jìn)而出現(xiàn)病斑、 枯萎等癥狀， 此時(shí)毛竹葉片樣本SPAD值整體呈現(xiàn)下降的趨勢。

圖3 不同狀態(tài)下毛竹葉片樣本SPAD變化趨勢

2.2 毛竹葉光譜特征變化分析

毛竹葉片的光譜曲線存在極為明顯的峰谷特征， 這與大部分綠色植被的光譜特性基本一致。 當(dāng)剛竹毒蛾危害發(fā)生時(shí)， 失葉狀態(tài)下的寄主會產(chǎn)生失綠、 缺水等病癥， 其光譜特征隨即產(chǎn)生顯著變化[圖4(a,b)]， 具體表現(xiàn)為： (1)健康竹葉在520～570 nm波長范圍內(nèi)因葉綠素吸收較少而形成了一個(gè)綠色波段的小反射峰， 稱之為“綠峰”。 (2)在640～700 nm波長范圍因葉綠素的強(qiáng)吸收形成了一個(gè)紅色波段吸收谷， 稱之為“紅谷”。 由于葉片的多孔薄壁細(xì)胞結(jié)構(gòu)特性對近紅外光強(qiáng)烈的反射， 形成光譜上的最高峰。 (3)當(dāng)毛竹遭受蟲害侵蝕后光譜發(fā)生明顯變化， 主要表現(xiàn)為可見光波段范圍內(nèi)的“綠峰”和“紅谷”逐漸消失， “紅邊”斜率減小， 近紅外波長反射率降低。 隨著蟲害程度加重， 這一現(xiàn)象更加明顯。 導(dǎo)致這一現(xiàn)象發(fā)生的主要原因是毛竹遭受蟲害侵蝕后， 竹葉葉綠素含量相對減少， 細(xì)胞結(jié)構(gòu)遭到破壞。 因此， 根據(jù)光譜間的差異， 進(jìn)行竹葉葉綠素含量的估算是可行的。

圖4 不同狀態(tài)下的毛竹葉片樣本(a)及其光譜信息(b)

相比之下， 小年葉片和健康葉片的光譜差異更為明顯。 在可見光—近紅外波段范圍內(nèi)， 小年葉片的反射率總體上高于健康及受害葉片； 至短波紅外波段區(qū)間， 其反射率雖仍高于健康葉片， 但與受害葉片的差異并不算太明顯， 雖然在1 450及1 940 nm波段上， 各健康狀態(tài)葉片依然會呈現(xiàn)出一定的規(guī)律， 但總體而言， 小年葉片的識別研究重點(diǎn)應(yīng)該落在可見光-近紅外范圍內(nèi)。

2.3 關(guān)系模型特征指標(biāo)篩選分析

2.3.1 基于全樣本的葉光譜特征指標(biāo)篩選分析

用于反演葉綠素的特征指標(biāo)眾多， 為減少計(jì)算量， 避免造成實(shí)驗(yàn)結(jié)果冗余、 復(fù)雜， 需首先對上述植被指數(shù)進(jìn)行篩選。 Pearson相關(guān)系數(shù)， 是一種線性相關(guān)系數(shù)， 通常用來度量兩個(gè)變量X和Y之間的相互關(guān)系， 其計(jì)算公式為[24]

(6)

式(6)中， Cov(X,Y)代表X與Y的協(xié)方差， Var(X)和Var(Y)分別表示X和Y的方差。 當(dāng)相關(guān)系數(shù)為1時(shí)，X與Y的關(guān)系為Y=aX+b，a>0； 當(dāng)相關(guān)性為-1時(shí)，X和Y的關(guān)系為Y=aX+b，a<0。 如果X和Y相互獨(dú)立， 相關(guān)系數(shù)為0(圖3、 圖4)。 本研究中將毛竹葉SPAD作為變量X， 將葉光譜特征指標(biāo)作為變量Y進(jìn)行Pearson相關(guān)分析。

圖5顯示， 植被指數(shù)對毛竹葉片SPAD的響應(yīng)最為顯著。 依照Pearson相關(guān)系數(shù)從大到小排列， 有5種植被指數(shù)的Pearson相關(guān)系數(shù)較為突出， 分別為： VOG2(-0.651**)，R515/R570(0.643**)， CIred(0.607**)， PRI(0.606**)及NDVI705(0.606**)， 故將上述5個(gè)植被指數(shù)作為SPAD估測毛竹葉片全體樣本SPAD模型的特征指標(biāo)。 而通過分析特征光譜與光譜微分， 可以發(fā)現(xiàn)在469和702 nm處原始光譜對SPAD的相關(guān)性較高， 而在760 nm處， 光譜微分則呈現(xiàn)較為明顯的相關(guān)性。 紅邊參數(shù)RES， REA， REDVI與SPAD的相關(guān)性較高而CARI， RERVI， RENDVI卻呈現(xiàn)較低的相關(guān)性。

圖5 全體葉片樣本葉光譜特征指標(biāo)Pearson相關(guān)分析

2.3.2 考慮不同危害等級的葉光譜特征指標(biāo)篩選分析

對全體葉片樣本劃分危害等級與小年葉片， 再次進(jìn)行Pearson相關(guān)分析， 結(jié)果如圖6(a—e)所示。

圖6 不同危害等級下葉片樣本葉光譜特征指標(biāo)Pearson相關(guān)分析

不同致?lián)p狀態(tài)下的毛竹葉片， 呈現(xiàn)趨勢為健康狀態(tài)及小年葉片樣本相關(guān)性較佳的特征指標(biāo)較多， 而與輕度、 中度及重度危害葉片樣本相關(guān)性較佳的特征指標(biāo)較少， 具體表現(xiàn)為： 與健康狀態(tài)的毛竹葉片相關(guān)性較佳的植被指數(shù)有： CIred(0.732**)， VOG2(-0.743**)， ARVI(0.549**)，R515/R570(0.551**)， DVI(0.522**)； 與小年毛竹葉片相關(guān)性較佳的植被指數(shù)有： PRI(0.706**)， NDVI705(0.704**)， VOG1(0.696**)， CIred(0.65**)。 與輕度、 中度危害葉片樣本SPAD相關(guān)性較高的指標(biāo)主要集中于紅邊參數(shù)， 與輕度蟲害相關(guān)的指標(biāo)有： RENDVI(0.438**)， RERVI(0.42**)及REDVI(0.264*)； 與中度蟲害相關(guān)的指標(biāo)與輕度蟲害相同， 其Pearson相關(guān)系數(shù)分別為： RENDVI(0.354**)， RERVI(0.36**)及REDVI(0.197*)； 與重度蟲害的毛竹葉片相關(guān)性較佳的植被指數(shù)有： VOG2(-0.443**)， CIred(0.435**)， NDVI705(0.371**)， PRI(0.359**)。 據(jù)此將上述指標(biāo)代入模型進(jìn)行計(jì)算。

圖7 4種模型對毛竹葉片SPAD檢測結(jié)果(R2與RMSE分別為模型擬合度與均方根誤差， 下同)

2.4 葉綠素與葉光譜特征模型構(gòu)建與關(guān)系分析

2.4.1 基于全樣本的模型構(gòu)建及葉綠素、 葉光譜關(guān)系分析

將上述選取的特征指標(biāo)分別代入多元線性回歸、 嶺回歸、 隨機(jī)森林及XGBoost 4種回歸模型之中， 以特征指標(biāo)作為自變量， SPAD作為因變量對全部葉片樣本進(jìn)行葉綠素估測， 模型評價(jià)指標(biāo)為擬合度(R2)及均方根誤差(root mean square error， RMSE)[圖7(a—d)]。 其中嶺回歸正則化參數(shù)α=0.01； 隨機(jī)森林回歸部分重要參數(shù)的調(diào)參結(jié)果為： random_state=90， n_estimators=141， max_depth=12， min_samples_split=5， min_samples_leaf=2， max_features=9； XGBoost部分重要參數(shù)的調(diào)參結(jié)果為： random_state=420， obj=reg: linear， max_depth=1， eta=0.02， gamma=50， nfold=4， lambda=1， alpha=0， colsample_bytree=1， colsample_bylevel=1， colsample_bynode=1。

從上述模型估測結(jié)果中， 可以發(fā)現(xiàn)基于總?cè)~片樣本環(huán)境下多元線性回歸模型對SPAD的檢測結(jié)果最佳，R2為0.753 7， RMSE為3.015 0， 其擬合方程為式(7)

SPAD=9.223+0.542×NDVI705+53.734×PRI+15.35×CIred+20.565×R515/R570+73.1×VOG2

(7)

XGBoost對SPAD的擬合結(jié)果較差，R2為0.711 3， RMSE為3.636 1。 值得注意的是， 嶺回歸及隨機(jī)森林的RMSE均較高， 分別為7.456 0與7.822 6， 說明這兩種模型對于SPAD的檢測較為不穩(wěn)定。 綜合上述結(jié)果， 無法清晰地看出毛竹葉SPAD與葉光譜特征的關(guān)系， 故本研究對全體樣本進(jìn)行危害等級及小年葉片劃分， 以求探究在不同危害等級下， 葉光譜關(guān)系模型估測結(jié)果發(fā)生何種變化。

2.4.2 考慮不同危害等級的模型構(gòu)建及葉綠素、 葉光譜關(guān)系分析

根據(jù)2.1中的特征指標(biāo)篩選結(jié)果， 將CIred， VOG2， ARVI，R515/R570， DVI作為健康葉片的特征指標(biāo)， RENDVI， RERVI， REDVI作為輕度危害及中度危害葉片的特征指標(biāo)， VOG2， CIred， NDVI705， PRI作為重度危害葉片的特征指標(biāo)， PRI， NDVI705， VOG1， CIred作為小年葉片的特征指標(biāo)。 將上述特征指標(biāo)分別代入多元線性回歸模型、 嶺回歸模型、 隨機(jī)森林回歸模型及XGBoosting回歸模型， 結(jié)果分別如圖8(a—t)所示。

圖8 4種模型對不同危害等級下毛竹葉片SPAD檢測結(jié)果

基于不同危害等級下對比上述模型的檢測結(jié)果， 主要有以下幾點(diǎn)特征： (1)分析模型R2， 多元線性回歸模型的檢測效果要優(yōu)于其他模型， 不同危害等級多元線性回歸的R2分別為健康葉片： 0.882 3、 輕度蟲害： 0.180 2、 中度蟲害： 0.360 4、 重度蟲害： 0.467 7、 小年葉片： 0.732 4； (2)分析模型RMSE， 多元線性回歸模型亦普遍較為穩(wěn)定， 不同危害等級多元線性回歸的RMSE分別為健康葉片： 1.638 8、 輕度蟲害： 3.335 4、 中度蟲害： 3.886 7、 重度蟲害： 2.601 8、 小年葉片： 2.375 4。 綜合判斷， 多元回歸分析可較好地估測小年及健康葉片SPAD， 對于輕度-中度-重度危害葉片的估測效果較差。

估測健康葉片的回歸方程為式(8)

SPAD健=-31.839+16.005×CIred+11.833×VOG2+49.892A×RVI+92.883×R515/R570+339.104×DVI

(8)

估測小年葉片的回歸方程為式(9)

SPAD小=40.049+9.036×PRI-21.091×NDVI705+27.81×VOG1+76.019×CIred

(9)

2.5 蟲害脅迫下毛竹葉SPAD關(guān)系模型的不確定性原因

為分析出現(xiàn)SPAD關(guān)系模型結(jié)果變化原因， 對毛竹葉片樣本進(jìn)行討論。 圖9顯示了不同危害等級下葉片樣本在3種高光譜特征指標(biāo)構(gòu)成的特征空間中的分布情況， 健康葉片與小年葉片樣本主要集中于右上視面與左上視面， 而輕度危害-中度危害-重度危害樣本則較為離散地分布于特征空間中， 這表明健康與小年?duì)顟B(tài)下的葉片葉光譜特征較為明顯， 而其他危害等級下的葉片在所選取的高光譜特征指標(biāo)中表現(xiàn)不突出， 故上述研究中3種危害等級葉片在Pearson相關(guān)分析中所表現(xiàn)出的相關(guān)性較差， 同時(shí)4種模型對其估測結(jié)果也較低。

圖9 全體葉片樣本及不同危害等級下葉片樣本葉光譜指標(biāo)特征空間

2.6 毛竹葉綠素估測模型結(jié)果分析

葉綠素是毛竹個(gè)體中參與光合作用的重要色素， 其通過捕捉光能轉(zhuǎn)變?yōu)榛瘜W(xué)能用于毛竹個(gè)體的生長和代謝。 葉綠素含量的高低可直接反映毛竹的生境狀況， 而蟲害脅迫下致使毛竹葉片生境及理化參數(shù)改變已成為阻礙估測毛竹葉片葉綠素含量的重要因素。 因此， 研究如何快速且準(zhǔn)確估測毛竹葉綠素含量及解析其與葉光譜特征關(guān)系變化具有重要意義。 本研究首先分析了不同危害等級毛竹葉片的劃分依據(jù)， 在此基礎(chǔ)上對全體葉片樣本及不同危害等級下的毛竹葉片進(jìn)行葉光譜特征指標(biāo)篩選， 再將篩選的指標(biāo)代入4種模型中。 由此發(fā)現(xiàn)， 多元線性回歸模型對全體毛竹葉片樣本SPAD具有較好的估測能力， 而針對不同危害等級下的葉片樣本， 模型的估測效果也隨之發(fā)生改變。 究其原因主要有以下兩方面：

首先， 無論特征指標(biāo)和葉片致?lián)p狀態(tài)如何變化， 嶺回歸和多元線性回歸的檢測結(jié)果總是十分接近， 并保留相同的變化趨勢。 事實(shí)上， 通過選取不同的正則化參數(shù)α， 發(fā)現(xiàn)當(dāng)α趨于更小值時(shí)， 嶺回歸的決定系數(shù)R2也更加趨近于多元線性回歸(圖10)。 這是因?yàn)閹X回歸是對多重共線性問題進(jìn)行回歸分析的一種統(tǒng)計(jì)方法， 是在多元線性回歸的損失函數(shù)上加上了正則項(xiàng)， 表達(dá)為ω的L2范式乘以正則化系數(shù)α， 其損失函數(shù)的完整表達(dá)式寫作式(10)

(10)

假設(shè)特征矩陣結(jié)構(gòu)為(m,n)， 系數(shù)ω的結(jié)構(gòu)是(1,n)， 則有

(XTX+αI)ω-XTy

(11)

式(11)中， RSS為殘差平方和，I為結(jié)構(gòu)為n×n的單位矩陣。 假設(shè)原本的特征矩陣中存在共線性， 則方陣XTX就會不滿秩。 對于存在“高度相關(guān)關(guān)系”的矩陣， 可以通過調(diào)大α控制參數(shù)向量ω的偏移， 即模型越不容易受到共線性的影響。 因此若一個(gè)數(shù)據(jù)集于嶺回歸中使用各種正則化參數(shù)取值下模型表現(xiàn)沒有明顯上升， 則說明數(shù)據(jù)沒有多重共線性， 反之亦然。

圖10 不同正則化參數(shù)下模型R2及RMSE變化

針對隨機(jī)森林回歸模型與XGBoost回歸模型， 盡管通過調(diào)參縮小了模型的過擬合問題， 但過擬合現(xiàn)象仍然存在。 本研究中所涉及的數(shù)據(jù)量較少且數(shù)值變化較為穩(wěn)定， 無法很好發(fā)揮出隨機(jī)森林與XGBoost模型的優(yōu)勢， 故在本研究中估測精度較多元回歸分析略低。

2.7 蟲害與毛竹葉片葉綠素及葉光譜特征關(guān)系分析

當(dāng)剛竹毒蛾危害發(fā)生時(shí)， 其幼蟲食出的缺刻使竹葉中的水分迅速流失， 導(dǎo)致葉綠素的合成速率下降。 隨著危害程度的不斷加劇， 殘余竹葉中的水分及葉綠素含量不斷降低， 進(jìn)而影響其光合效率。 由于植被吸收的能量通常比用于光合作用的能量更多， 其自身的各種光保護(hù)機(jī)制使植被能夠釋放存在潛在威脅的多余能量。 但當(dāng)光合效率減弱到一定程度時(shí)， 光能吸收-釋放的平衡被打破， 過多的能量將導(dǎo)致致命的感光氧化。 與此同時(shí)， 光合效率減弱進(jìn)一步導(dǎo)致竹體內(nèi)的水分無法被有效耗解， 由此引發(fā)惡性循環(huán)， 竹節(jié)內(nèi)逐漸產(chǎn)生積水， 而竹冠則不斷干枯。 因此， 選擇可反映以上表征變化的指標(biāo)是確定其蟲害脅迫程度的關(guān)鍵。 本研究從原始光譜、 植被指數(shù)、 紅邊參數(shù)及光譜微分4個(gè)光譜特征指標(biāo)入手， 從多角度指示毛竹SPAD在不同危害程度下葉光譜特征的變化差異， 由此發(fā)現(xiàn)， 隨著蟲害程度的加重， 葉光譜特征指標(biāo)的相關(guān)性呈現(xiàn)顯著下降趨勢， 致使模型對存在蟲害脅迫的葉片SPAD估測效果不佳。 可見， 當(dāng)毛竹葉SPAD與葉光譜特征的關(guān)系趨向紊亂時(shí)， 可預(yù)示蟲害發(fā)生。

3 結(jié) 論

基于實(shí)測毛竹葉片樣本SPAD及葉光譜數(shù)據(jù)， 采用Pearson相關(guān)法對不同危害等級下毛竹葉片SPAD的高光譜衍生指標(biāo)進(jìn)行篩選， 據(jù)此代入多元線性回歸、 嶺回歸、 隨機(jī)森林回歸及XGBoost 回歸4種模型中對毛竹葉片SPAD信息予以估測， 基于模型估測指標(biāo)分析葉綠素與葉光譜特征關(guān)系變化， 得出以下結(jié)論：

(1)隨著蟲害程度的加劇， SPAD普遍呈現(xiàn)下降的趨勢。

(2)蟲害脅迫下葉光譜特征發(fā)生明顯變化， 其可見光波段范圍內(nèi)的“綠峰”和紅谷逐漸消失， “紅邊”斜率減小， 近紅外波長反射率降低。

(3)基于總樣本， 根據(jù)Pearson相關(guān)法篩選出估測全體葉片樣本SPAD的最佳指標(biāo)為： VOG2(-0.651**)，R515/R570(0.643**)， CIred(0.607**)， PRI(0.606**)及NDVI705(0.606**)。 依據(jù)上述指標(biāo)擬合出的最佳模型為多元線性回歸模型，R2為0.753 7， RMSE為3.015 0， 其擬合方程為： SPAD=9.223+0.542NDVI705+53.734PRI+15.35CIred+20.565R515/R570+73.1VOG2。

(4)基于不同危害等級下的葉片樣本， 多元回歸模型是最佳估測模型， 隨著危害等級變化多元線性回歸的估測效果也隨之改變， 分別為， 健康葉片：R2=0.882 3， RMSE=1.638 8； 輕度蟲害葉片：R2=0.180 2， RMSE=3.335 4； 中度蟲害葉片：R2=0.360 4， RMSE=3.886 7； 重度蟲害葉片：R2=0.467 7， RMSE=2.601 8； 小年葉片：R2=0.732 4， RMSE=2.375 4。 可以看出， 模型對于健康葉片與小年葉片的SPAD估測效果較好； 對于輕度危害—中度危害—重度危害葉片的SPAD估測效果較差。

(5)通過構(gòu)建由不同葉光譜特征組成的特征空間， 分析不同危害等級下葉光譜特征關(guān)系不確定的原因， 為小年與健康葉片樣本SPAD較為集中地分布在特征空間中的不同視面， 而其他危害等級葉片樣本的SPAD則呈現(xiàn)較為離散地分布， 致使樣本特征不夠突出。

綜上所述， 通過本方法所選取的敏感性指標(biāo)并利用多元線性回歸模型進(jìn)行計(jì)算， 對毛竹葉片SPAD具有良好的估測效果， 而通過分析葉光譜特征與蟲害發(fā)生關(guān)系， 則可發(fā)現(xiàn)在不同危害等級下光譜特征表征也隨之改變， 同時(shí)基于不同危害等級亦可發(fā)現(xiàn)當(dāng)毛竹葉SPAD與葉光譜特征的關(guān)系趨向紊亂時(shí)， 則預(yù)示著蟲害的發(fā)生。 本研究可為毛竹健康監(jiān)測及其剛竹毒蛾危害響應(yīng)與預(yù)警研究提供理論支持。