以福建省南平市顺昌县为试验区, 该县位于闽西北部, 地理坐标为117° 29'— 118° 14'E, 26° 38'— 27° 12'N, 其主要受亚热带海洋性季风影响, 气候温和, 干湿季明显。 顺昌县的森林资源丰富, 是我国南方的重点林区, 有中国竹子之乡的美誉。 截至2018年底, 全县林地面积达16.7万hm², 其中竹林面积为4.4万hm², 毛竹林面积约3.21万hm²。 近年来, 顺昌县各类森林虫害事件频发, 诸如刚竹毒蛾、 马尾松毛虫等主要森林害虫的发生面积居高不下, 严重损害了顺昌县林业经济的可持续发展。

1.2.1 叶片光谱数据采集

于2019年8月20日— 23日赴实地开展调研工作(该时期处于刚竹毒蛾第一代幼虫末期, 其体型最大, 食叶量最多)。 选用美国ASD公司(Analytical Spectral Devices)生产的Field Spec3地物光谱仪及配套的植物光谱探测器(Plant probe)采集不同竹龄及健康程度的毛竹叶片(健康、 受刚竹毒蛾危害、 小年竹叶)光谱数据。 波长范围为350~2 500 nm, 采样间隔为1.4 nm(350~1 000 nm)和2 nm(1 001~2 500 nm), 重采样间隔为1 nm。 为确保光谱数据的准确性, 每隔20 min进行一次白板校正。 分别测定各叶片的叶尖、 叶中、 叶基部的光谱值, 每个部位取5次测定的平均值, 随后利用该仪器配套的光谱处理软件(View Spec Pro)剔除明显错误的光谱数据, 最后求取有效数据的算数均值, 从而得到该叶片的光谱反射率数据。

1.2.2 刚竹毒蛾虫害等级确定标准

由于现有的林业虫害等级判定标准均以小班、 株为单位^[6], 考虑到尺度差异, 采用株— 叶片的综合评判方式进行虫害等级划分。 首先以虫口数量(无危害: < 10条、 轻度危害: 10~30条、 中度危害: 31~80条、 重度危害: > 80条)和失叶率(无危害: 0%~10%、 轻度危害: 10%~25%、 中度危害: 25%~50%、 重度危害: > 50%)为分级依据, 于各虫害等级毛竹冠层的上、 中、 下部采集叶片; 随后依据各叶片的叶损量(即病斑、 缺刻面积之和占叶片原始面积的比例), 进一步确定叶片尺度下的刚竹毒蛾虫害等级(无危害: 0%~10%、 轻度危害: 10%~25%、 中度危害: 25%~50%、 重度危害: > 50%)。 对于缺刻程度较轻的叶片, 沿叶片边缘将其补齐后计算其原始面积; 对于缺刻程度较为严重的叶片, 则根据其残存部分的形状, 并以最近位置的完整叶片为参照, 补齐后计算该叶片的原始面积。 叶损量的计算公式为式(1)

LLC=ALLA×100%(1)

式中: LLC为叶损量; A_LL为病斑、 缺刻面积之和, A为叶片原始面积。

1.2.3 高光谱指数计算

植被指数可作为描述植被状态的重要参量, 选择较为常用的20种植被指数进行分析。 根据已有研究^{[7, 8]}, 在一般植被的光谱曲线当中, “ 绿峰” (因叶绿素对绿光的弱吸收而形成的小反射峰)、 “ 红谷” (因叶绿素对红光的强吸收而形成的反射率低谷)、 红边(植被在680~780 nm间反射率快速变化的区域)范围内的光谱信息可为其健康程度评估提供丰富的参考依据。 经过综合考虑, 尝试以“ 绿峰” — “ 红谷” 、 “ 红谷” — 红边范围内的光谱导数变化为切入点, 建立植被健康程度评估指数, 其构建思路为: (1)求取各叶片在550~780 nm波长范围内的光谱一阶导数; (2)运用遗传算法选择特征波段, 重复运行该算法10次, 取入选8次及以上的波段; (3)根据筛选结果, 以健康叶片在红光波段上的反射率拐点(667 nm)为分界点, 利用归一化差值公式建立指数。 经过筛选, 入选的波段分别为572, 576, 618, 652, 712, 717和777 nm, 由此产生12组待验指数。 所用到及构建的指数名称及公式参见表1。

表1

Table 1

表1(Table 1)

表1 用于评估毛竹叶片健康程度的高光谱指数 Table 1 Hyperspectral indices for estimating the health level of Moso bamboo leaves 指数名称 计算公式 参考文献 NDVI (R800-R670)/(R800+R670) [9] PSRI (R680-R500)/R750 [9] MCARI [(R700-R670)-0.2(R700-R550)]× (R700/R670) [9] SIPI (R800-R445)/(R800-R680) [10] GNDVI (R750-R550)/(R750+R550) [10] NDCI (R762-R527)/(R762+R527) [10] NPCI (R680-R430)/(R680+R430) [10] mSR705 (R750-R445)/(R705+R445) [10] VOG 1 R740/R720 [11] VOG 2 (R734-R747)/(R715+R726) [11] VOG 3 (R734-R747)/(R715+R720) [11] RVSI (R712+R752)/2-R732 [12] PRI (R531-R570)/(R531+R570) [13] NDWI (R857-R1 241)/(R857+R1 241) [14] LCI (R850-R710)/(R850-R680) [14] DSWI (R803+R549)/(R1 659+R681) [14] NDLI [log(1/R1754)-log(1/R1680)]/ [log(1/R1 754)+log(1/R1 680)] [15] NDNI [log(1/R1 510)-log(1/R1 680)]/ [log(1/R1 510)+log(1/R1 680)] [15] CIgreen (R800/R550)-1 [16] CIrededge (R800/R720)-1 [16] FDVI572, 712 (DR572-DR712)/(DR572+DR712) 本研究 FDVI572, 717 (DR572-DR717)/(DR572+DR717) 本研究 FDVI572, 777 (DR572-DR777)/(DR572+DR777) 本研究 FDVI576, 712 (DR576-DR712)/(DR576+DR712) 本研究 FDVI576, 717 (DR576-DR717)/(DR576+DR717) 本研究 FDVI576, 777 (DR576-DR777)/(DR576+DR777) 本研究 FDVI618, 712 (DR618-DR712)/(DR618+DR712) 本研究 FDVI618, 717 (DR618-DR717)/(DR618+DR717) 本研究 FDVI618, 777 (DR618-DR777)/(DR618+DR777) 本研究 FDVI652, 712 (DR652-DR712)/(DR652+DR712) 本研究 FDVI652, 717 (DR652-DR717)/(DR652+DR717) 本研究 FDVI652, 777 (DR652-DR777)/(DR652+DR777) 本研究 Note: Ri is the reflectance in i band; DRi is the first derivative reflectance in i band 注: Ri表示波段i处的光谱反射率; DRi表示波段i处的一阶导数光谱反射率

表1 用于评估毛竹叶片健康程度的高光谱指数 Table 1 Hyperspectral indices for estimating the health level of Moso bamboo leaves

XGBoost(eXtreme gradient boosting)是在GBDT(gradient boosting decision tree)的基础上改进的一种提升树模型, 其在当前的机器学习和数据挖掘领域中有着深远的影响。 该模型的原理是让新的基模型(CART回归树)不断地拟合先前模型的偏差, 直至模拟效果达到最佳。 相较于传统的GBDT, XGBoost将目标函数泰勒展开至二阶, 保留了目标函数的更多有效信息, 并添加了正则化项用以避免过拟合; 此外, XGBoost还特别设计了针对稀疏数据的算法, 无需对缺失特征进行填充预处理。 该模型的理论框架及具体说明可参见文献[17]。 研究表明, 经过适当的调参, XGBoost可得到比随机森林更好的分类结果^{[18, 19]}。 XGBoost模型的关键参数有:

(1)n_estimators: 生成分类树的个数, 即迭代次数;

(2)learning_rate: 每次迭代的步长, 即学习率;

(3)max_depth: 每棵树的最大深度, 用于避免过拟合;

(4)gamma: 节点分裂所需的最小损失函数下降值;

(5)subsample: 训练样本占整体样本的比例;

(6)objective: 目标函数, 本研究multi: softmax多分类器。

在建立模型的过程中, 为了得到更好的模拟效果, 需剔除贡献度较低的特征。 因此, 在保证每次模拟精度为当次最高的前提下, 根据模型给出的变量重要性排序, 依次剔去贡献度较低的特征指数, 直至模拟效果达到最佳。 此外, 为了避免实验的片面性, 对样本数据进行随机划分(重复6次), 得到6个试验集(样本量442个)和验证集(样本量110个), 分别设为实验1— 6, 最后根据各组的检验结果对模型的效果进行综合评价。 以上过程基于Python 3.6.5平台实现。

2.1.1 不同健康程度毛竹叶片的光谱特征

毛竹叶片的光谱曲线存在极为明显的峰谷特征, 与大部分绿色植被的光谱特性基本一致。 当刚竹毒蛾危害发生时, 寄主会产生失绿、 缺水等病症, 其光谱特征随即产生显著变化(图1), 具体表现为: (1)在可见光范围内, 虫害胁迫引发的色素变化使得“ 绿峰” 及“ 红谷” 的反射率变化尤为明显, 表现为受害叶片在绿光波段的反射率随虫害等级的上升而降低, 红光波段的反射率则逐渐上升; (2)各虫害等级叶片在近红外波段上的光谱差异则更为显著, 由于刚竹毒蛾幼虫在取食的过程中破坏了叶片的组织结构, 导致光在叶片缺刻附近的散射程度被削弱, 同时在色素、 水分等因素的多重影响之下, 导致其在近红外波段的反射率随虫害等级的上升而不断降低; (3)至短波红外时, 失水使得受害叶片的反射率普遍高于健康叶片, 尤其在1 450及1 940 nm这两个典型的水汽吸收波段的差异尤为显著。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 不同虫害等级毛竹叶片光谱反射率曲线Fig.1 Spectrum of moso bamboo leaves in different health levels

小年叶片和健康叶片的光谱差异则更为明显。 在可见光— 近红外波段范围内, 小年叶片的反射率总体上高于健康及受害叶片; 至短波红外, 其反射率虽仍高于健康叶片, 但与受害叶片的差异则并不算太明显, 虽然在1 450及1 940 nm波段, 各健康状态叶片依然会呈现出一定的规律, 但总体而言, 小年叶片的识别研究重点应该落在可见光— 近红外范围内。

2.1.2 不同受害类型叶片的光谱特征

不同受害类型可分为缺刻型叶片、 红褐色病斑型叶片、 灰白色木质化病斑型叶片三类[图2(a, b, c)], 各类叶片光谱特征亦具有极为显著的差异。 对于每一种受害类型, 各选5片有代表性的叶片, 取其光谱均值进行分析。 由图3可以看出, 较之健康叶片, 缺刻型叶片的光谱并未表现出太大的差异; 相比之下, 红褐色病斑型叶片的光谱特征则产生了一定的变化, 其在可见光部分的反射率高于缺刻型叶片, 且红光及黄光波段的反射率出现了一定程度的升高, 但总体上依然保有植被的基本光谱特征; 灰白色病斑型叶片在可见光-近红外波长范围内的光谱峰谷特征已经消失, 且反射率总体上高于另外两种受害类型叶片(其近红外波段的反射率甚至高于健康叶片), 已基本失去了植被的基本光谱特征。 由此可见, 即便是相同虫害等级的叶片, 其光谱特征仍有可能会存有较大差异。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 不同受害类型的叶片 (a): 缺刻型叶片; (b): 红褐色病斑型叶片; (c): 灰白色木质化病斑型叶片Fig.2 Leaves with different damage types (a): Indentation-only leaves; (b): Red-brown disease spots leaves; (c): Gray-white disease spots leaves

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 不同类型病斑叶片的光谱特征Fig.3 Spectrum of moso bamboo leaves with different damage types

2.2.1 入选的特征变量及其重要性排序

经过反复调试及筛选, 最终确定模型的最佳输入变量为FDVI_{576, 717}, VOG 1, PRI, DSWI, NPCI, NDWI和RVSI。 总体而言, 这些特征变量呈现出了较为清晰的虫害响应规律, 不同受害等级的叶片之间存在着较大的差异性。 例如, 随着受害等级的加重, PRI和FDVI_{576, 717}呈现明显的下降趋势。 但从各类别叶片的特征变量值及其方差来看, 即便是相同虫害等级的叶片, 其值仍存在较大的不同。 根据实验1— 6的变量重要性排序可知, 7个特征变量的贡献度依次为PRI> FDVI_{576, 717}> NPCI> DSWI> VOG 1> RVSI> NDWI[图4(a— h)]。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 不同健康程度叶片的特征变量值及变量重要性排序 (a): FDVI576, 717; (b): VOG1; (c): RPI; (d): DSWI; (e): NPCI; (f): NDWI; (g): RVSI; (h): 各变量的重要性得分Fig.4 The individual characteristic variables values for leaves with different health levels and the feature importance score (a): FDVI576, 717; (b): VOG1; (c): RPI; (d): DSWI; (e): NPCI; (f): NDWI; (g): RVSI; (h): The importance score of each characteristic variable

2.2.2 模型检测效果分析

模型的检测结果显示, 实验1— 6的检测精度分别为73.64%, 78.18%, 73.64%, 72.73%, 73.64%和74.55%, 总平均精度为74.39%(表2)。 细分各健康程度叶片的分级结果, 健康叶片的识别精度最高, 达到了94.55%; 其次为重度危害叶片, 识别精度为84.12%; 再次为轻度危害叶片, 识别精度为74.93%; 小年叶片的识别精度则为71.10%; 中度危害叶片的识别精度是最低的, 仅有33.48%。 总体而言, 模型对健康、 轻度危害、 重度危害及小年叶片的识别精度基本达到了预期, 但对中度危害叶片的识别效果仍有待提升。

表2

Table 2

表2(Table 2)

表2 模型检测结果 Table 2 Model detection results 实验 分级精度/% 单次实 验精度 /% 健康 轻度 危害 中度 危害 重度 危害 小年 1 82.76 75.00 52.94 83.33 50.00 73.64 2 96.43 93.75 25.00 88.57 60.00 78.18 3 100.00 64.71 27.27 87.50 80.00 73.64 4 96.67 88.24 21.74 74.19 88.89 72.73 5 91.43 63.16 42.11 80.77 72.73 73.64 6 100.00 64.71 31.82 90.32 75.00 74.55 平均分级精度/% 94.55 74.93 33.48 84.12 71.10 总平均精度/% 74.39

表2 模型检测结果 Table 2 Model detection results

当刚竹毒蛾危害发生时, 其幼虫食出的缺刻使得竹叶中的水分迅速流失, 导致叶绿素的合成速率出现下降。 随着失叶程度的不断加剧, 寄主用于光合作用的叶面积减少, 同时叶绿素、 水分含量的下降进一步削弱了寄主的光合作用能力。 由于植被吸收的能量通常比用于光合作用的能量要多, 其自身的各种光保护机制使植被能够释放存在潜在威胁的多余能量, 但当光合效率减弱到一定程度时, 光能吸收-释放的平衡将被打破, 过多的能量将导致致命的感光氧化。 此外, 光合效率的减弱进一步导致竹体内的水分无法被有效耗解, 由此引发恶性循环, 竹节内逐渐产生积水, 而竹冠则不断干枯。 因此, 选择可反映以上表征变化的指标是确定虫害胁迫程度的关键。 根据最终入选的7个特征变量, PRI指示了不同健康程度叶片的光能利用率变化; FDVI_{576, 717}及NPCI分别反映了各叶片间的色素含量差异; DSWI和NDWI代表了各叶片的水分含量情况; 而VOG 1及RVSI则突出了不同健康程度叶片的红边光谱差异。 本研究选取的特征变量是较为合理的。

根据模型的检测结果显示, 其总平均精度约在74%, 从传统分类的角度看, 此精度似乎并不算高, 分析其原因: (1)不同毛竹个体间的养分含量、 光谱特征本就存有差异, 而当研究尺度落在叶片上时, 这种差异会进一步被放大, 例如大年和小年竹叶或是同株毛竹上新叶与老叶的光谱差异均有可能大于健康叶片和受害叶片; (2)受害竹叶的光谱特征与其距离虫害发生时间的长短有关, 当虫害初发时, 竹叶中的各种养分含量尚未产生明显变化, 而随着失叶程度的不断加剧, 竹体的养分循环系统将被完全破坏, 由此导致其中的养分元素含量产生显著变化, 进而出现病斑、 枯萎等症状, 此时的叶片光谱特征与初受害时期便有着明显的差异, 这亦是导致中度危害叶片成为误分“ 重灾区” 的主要原因。 因此, 尽管本文选用了当前在机器学习领域中表现极为出色的XGBoost模型, 但得到的精度结果仍旧没有达到令人满意的程度, 不过考虑到叶片尺度研究所面临的多种不确定性因素, 该结果也显得相对可以接受。 另外, 本研究重点也并非在此, 而旨在探索不同健康状况叶片的光谱变化机理, 从而为后续冠层尺度的研究奠定理论基础。