基于高光谱特征的雅氏落叶松尺蠖虫口密度估算

白力嘎, 黄晓君,3,*, Ganbat DASHZEBEGD, Mungunkhuyag ARIUNAAD,Tsagaantsooj NANZADD, Altanchimeg DORJSUREN, 包 刚, 佟斯琴,包玉海, 银 山, Enkhnasan DAVAADORJ

(1. 内蒙古师范大学地理科学学院, 呼和浩特 010022; 2. 内蒙古自治区遥感与地理信息系统重点实验室, 呼和浩特 010022;3. 内蒙古自治区蒙古高原灾害与生态安全重点实验室, 呼和浩特 010022; 4. 蒙古国科学院地理与地质研究所, 蒙古国乌兰巴托 15170;5. 蒙古国科学院综合实验生物学研究所, 蒙古国乌兰巴托 13330)

森林在生态系统中扮演着至关重要的角色。它不仅在生物圈能量维持和生物多样性保持中发挥重要作用，而且对自然界矿元素和水循环中亦起到不可替代的作用(李双成和杨勤业, 2000)。森林害虫暴发会使林木大量失叶，生理功能逐步衰退，甚至死亡。这导致森林生态功能急剧减弱，而造成森林生态系统严重破坏，并威胁生态安全。可见，监测森林虫害对害虫防控及森林生态安全保护具有实际意义。虫口密度能够直接表征森林害虫的发生程度，它可作为森林虫害监测的重要指标。从而开展森林害虫虫口密度估算研究具有重要意义。

雅氏落叶松尺蠖Erannisjacobsoni是蒙古高原典型落叶松Larixsibirica害虫。其危害主要集中在每年6月份的幼虫期，在这期间其幼虫暴食针叶，落叶松针叶从顶部开始脱落，直至全部损失；同时，落叶松针叶水分含量和叶色素含量等生理指标急剧下降，生理功能衰退，光合作用减弱，导致森林生态功能严重衰退，甚至引起生态环境问题。多年来，雅氏落叶松尺蠖灾害在蒙古国落叶松林频发，严重破坏了落叶松林生态系统。更值得关注的是该虫在蒙古国具有向东或东南方向传播态势，而我国大兴安岭林区恰好达到此虫中度适生条件，存在入侵风险(黄晓君等, 2018)，并且一旦入侵极容易形成优势种群，将严重破坏我国森林生态系统，造成巨大经济损失。可见，监测该森林虫害对森林生态安全保护具有重要意义。在此背景下，雅氏落叶松尺蠖害虫虫口密度的及时并大面积监测势在必行，以期防控该虫灾害。

因虫口数量指标可直接描述害虫发生的严重程度，虫口数量相关调查与虫口密度研究始终受到人们广泛关注(Cayuelaetal., 2014)。主要有以下方面研究：一是虫口数量指标调查与计算研究和虫口数量指标与其他因素之间的关系研究，如：Valmorbida等(2018)采用人工计数法统计蛴螬害虫虫口密度，分析了害虫地理分布特征；Jyoti等(2020)采用棋盘式取样法，统计西藏飞蝗Locustamigratoriatibetensis不同年龄组分类的数量与昆虫总数，并研究了该虫在草地和农田中的空间分布格局；Saghafipour等(2020)利用诱蝇器诱捕白粉虫，并统计其数量，分析了气候变化与白粉虫虫口暴发的关系；Kutuk等(2010)在0.25 m2的若干样方内，采用人工计数法统计了麦扁盾蝽Eurygasterintegriceps越冬的和新生代的虫口密度，研究了田间越冬成虫密度、卵和成虫寄生率对该虫新生代虫口密度的影响。二是虫口密度预测研究和虫口密度与危害程度的关系研究，如：许向利等(2012)采用棋盘式10点取样法调查麦田地下害虫虫口密度，并探析地下害虫对麦苗危害程度与虫口密度的关系；Nyeko等(2010)研究了成年瘿姬小蜂Leptocybeinvasa虫口密度的时空变化并评估了桉树上的瘿姬小蜂侵染情况；Night等(2011)通过调查卢旺达木薯害虫数量，分析了害虫发生率和严重程度；Debuse等(2019)调查种植园蛀干害虫危害数量，探究了害虫发生率和严重程度的相关性。

综上所述，虫口密度获取手段主要以传统的人工调查为主，其有很多局限性，如费时费力、成本高以及不能满足用户快速大范围监测的需求，从而亟待有个快速、大范围监测虫口密度的方法框架来满足林业灾害防控部门的需求。随高光谱遥感快速发展，基于区域尺度的虫口密度快速估算将成为可能。随害虫虫口密度逐渐上升，树木失叶量随之增加，叶片光合作用减弱，蒸腾速率下降，内部生化组分和外观形态发生变化，而高光谱遥感波段窄、通道多，能够捕捉到森林植被内部生化组分和外观形态变化特征。这说明虫口密度与受害林木冠层光谱特征之间存在着密切关系，利用高光谱特征实现虫口密度估算是可能的。基于此，本研究通过蒙古国雅氏落叶松尺蠖暴发典型区的落叶松光谱实测数据和虫口密度调查数据，分析虫口密度与高光谱特征的关系，提取敏感光谱特征，结合回归算法建立虫口密度高光谱估算模型，揭示了高光谱特征对雅氏落叶松尺蠖虫口密度的估算潜力，为森林虫口密度大范围遥感监测提供了实验理论基础。

1 材料与方法

1.1 试验区概况

以蒙古国后杭爱省伊和塔米尔、青克尔、巴图青格乐以及肯特省宾德尔雅氏落叶松尺蠖暴发区作为4个试验区，共选取了110株不同程度受害的西伯利叶落叶松Larixsibirica样本树(如图1所示)，于2016年6月和2019年6月开展虫口密度和高广谱数据调查。

图1 蒙古国杭爱省和肯特省4个试验区样本树分布图

1.2 数据获取与计算

1.2.1虫口密度数据获取：样本树雅氏落叶松尺蠖虫口数量调查，设样本树外观形状为锥形(图2)，并计算了虫口密度(样本树总虫数与样本树投影面积的比值)。其具体计算方法如下：

图2 样本树模拟图

第1步,统计样本树树节n，测量相邻树节之间的高度h以及第1树节到树顶的高度H，然后测量第1树节平行轮生枝垂直于地面的投影长度R1，并计算其平均长度。第2步,通过人工计数法获得第1树节上的虫口数量，然后结合虫口数量削减系数，计算剩余各树节上的虫口数量。最后,通过各树节上的虫口数量计算虫口总数，并结合第1树节投影面积，计算当前样本树虫口密度，记为POPD，详情见文献(黄晓君, 2019)。计算结果显示，POPD最大值为302，最小值为10，平均值为131.36，标准差为72.56，其符合正态分布(图3)，说明POPD能够满足建模需求。

图3 正态性检验

1.2.2高光谱数据获取：针对不同程度受害的样本树，与虫口数量调查同步，利用SVC HR-1024(3.5 nm@350～1 000 nm， 9.5 nm@1 000～1 850 nm， 6.5 nm@1 850～2 500 nm)和ASD FieldSpec4(3 nm@350～1 000 nm， 10 nm@1 000～2 500 nm)地物光谱仪，测量了其冠层光谱反射率。实测工作在晴朗无云无风的天气条件下，于北京时间10∶30-14∶30期间进行的。另外，为避免地面植被对测量值的干扰，在目标针叶下方放置了黑布。在作业中，光谱仪探头垂直向下，视场角为25°，离测量目标约为20 cm。样本树垂直方向上分为上、中、下3个层次，在各层次上选择一支典型树枝进行测量20次，每次测量前后均使用校正白板对光谱仪进行校正(白板的反射率为1)，将20×3次的光谱反射率的平均值作为样本树冠层光谱反射率。

1.2.3微分光谱反射率计算：实测光谱的过程中受外界和仪器噪声干扰，会产生异常和重叠光谱反射率的现象。从而在剔除异常光谱和重叠光谱之后，取光谱反射率平均值，再借助Smooth函数对其进行平滑处理，然后通过ENVI软件的Spectral Math工具，对各样本树平滑光谱反射率求一阶导数，得到微分光谱反射率，记为DSR。光谱微分变换可消除大气效应,降低背景干扰,净化光谱数据，从而突出与植被本质变化有关的光谱吸收和反射特征信息(黄晓君, 2019)。

1.2.4光谱指数计算：(1)原始光谱指数(OSI)：考虑雅氏落叶松尺蠖数量变化会使林木针叶叶绿素、水分和外观形态发生变化的关系，并结合前人研究经验，选择了30个与受害过程林木生化组分和外观状态变化有密切关系的光谱指数(Chengetal., 2011; Heetal., 2018; 张卓然, 2018)。前人已构建的光谱指数设为原始光谱指数(OSI)，如表1所示。通过Matlab编程计算原始光谱指数。

表1 原始光谱指数(OSI)及其定义

(2)改进型光谱指数(MSI)：为进一步挖掘OSI的潜在POPD估算能力，对原始光谱指数的波长进行调整，以期实现光谱指数改进。原始光谱指数的改进方法是在OSI现有的波长两侧扩展20 nm，以1 nm为步长，对构建指数的所有波长进行排列组合，得到一组新的光谱指数集SIj。在此基础上，分别计算SIj与POPD的拟合优度(R2)，然后拟合优度最高值对应的光谱指数设为改进型光谱指数(MSI)。光谱指数集SIj获得的具体算法如下：

(1)

[Yi,Yi-1,…,Y1]=ndgrid(Xi,Xi-1,…,X1)

(2)

Zji=[Y1(:),Y2(:),…,Yi(:)]

(3)

SIj=[D1,D2,…,Dj]

(4)

式中X1,X2和Xi分别为构建光谱指数的第1， 2和i波段波长组成的行矩阵(相邻元素间隔，即波长间隔为1 nm)；Am1，Bn1和Ck1为最小波长值；Amm，Bnn和Ckk为最大波长值；ndgrid()为生成多维矩形网格函数；Y1,Y2, …,Yi为通过ndgrid()函数复制X1,X2, …,Xi等行矩阵并生成的i维矩形网格阵列(每格网值代表波长值)；Zji为构建光谱指数的各波段所有波长排列组合堆积的j行i列矩阵；SIj是以不同波长组合所构建的D1,D2, …,Di等j个光谱指数形成的列矩阵，即对于一个光谱指数的所有可能波长组合构建的光谱指数。利用上述方法获得了所有原始光谱指数的改进型光谱指数。

1.3 数据分析与建模

1.3.1光谱特征敏感性分析：从POPD数据来看，POPD变化有一定规律，其具体表现为随林木受害程度的增加，POPD先逐步上升，等到一定峰值后，呈下降态势。显然，POPD随受害程度的加剧呈现了多项式的倒“U”型曲线特征。基于此，选择样本训练集的OSI, MSI和DSR等高光谱特征，借助最小二乘法思想(最小化误差的平方和寻找数据的最佳函数匹配)，分别与相应POPD进行了多项式曲线拟合，并获得了拟合优度(R2)。当拟合的R2越高说明对应的光谱特征敏感性越强，反之亦然。

1.3.2敏感光谱特征提取：(1)敏感DSR选取：通过Matlab编程，在350～1 800 nm范围内提取DSR与POPD的拟合优度(R2)最高时对应的DSR，设为最敏感特征，并作为单个高光谱特征的估算模型输入变量。与此同时，利用Matlab的Findpeaks函数，以10 nm为步长，提取R2峰值，获得峰值对应的一组DSR，作为初始敏感特征组。在此基础上，借助连续投影算法(SPA)剔除初始敏感特征组内冗余信息，保留少数有意义的敏感DSR，并作为多个高光谱特征的估算模型输入变量。SPA是一种使建模变量共线性最小化的前向变量选择算法，它的优势在于提取初始数据中的少数几列数据就可以概括绝大部分光谱变量的信息，能够排除初始光谱矩阵中冗余的信息，减少数据量，最大程度地减少了信息重叠，可用于光谱特征的筛选，提高建模效率。Findpeaks函数和SPA算法详情见文献(黄晓君等, 2019)。

(2)敏感MSI选取：通过Matlab编程选取MSI与POPD的拟合优度(R2)最高时对应的MSI，确定为最敏感MSI，并作为单个高光谱特征的估算模型输入变量。与此同时，运用SPA算法，对由30种MSI组成的光谱指数集进行降维处理，选取少数的敏感MSI，作为多个高光谱特征的估算模型输入变量。

2 结果

2.1 高光谱特征对POPD的敏感性分析

2.1.1高光谱特征对不同等级POPD的响应：为分析高光谱特征对POPD的响应特征，根据落叶松健康和不同受害程度，对POPD分了4个等级，即健康(POPD： 0～20)、轻度(POPD： 20～120)、中度(POPD： 120～230)和重度(POPD： 230～最大值或最大值～最小值)，并绘制了POPD的4个等级对应的高光谱特征变化线(图4)。由图可知，POPD由健康→轻度→中度→重度等级变化过程中，DSR和MSI(归一化值)具有明显的层次变化响应。具体表现为：对DSR而言(图4: A)，随POPD等级的增加，DSR在502～545和678～761 nm波段内呈下降趋势，而在554～594, 1 111～1 169和1 284～1 344 nm波段内呈上升态势，其中在678～761 nm波段内DSR层次变化尤其明显，并出现了波峰，说明该波段中DSR对POPD变化的响应较显著。对MSI而言(图4: B)，随POPD等级的加重，MSI的ARI1, RGI和CTR2等光谱指数归一化值逐步上升，而其余MSI的光谱指数归一化值呈下降趋势，表明MSI随POPD等级变化具有明显的层次变化规律。这是因为POPD等级增加使落叶松冠层结构发生变化，进而影响林木冠层光谱反射率。可见，利用DSR和MSI估算POPD是可行的。

2.1.2基于拟合优度的敏感性分析: (1)微分光谱反射率：为进一步分析DSR对POPD的敏感性，将DSR逐波段(350～1 800 nm)与POPD进行多项式曲线拟合，得到了各波长对应的拟合优度(R2)，并与其线性拟合优度(R2)做了比较(图5)。结果显示，DSR与POPD的多项式曲线拟合效果较好，线性拟合效果较差。在505～531, 553～585, 691～755, 1 275～1 351, 1 459～1 625和1 726～1 775 nm波段DSR多项式曲线拟合的R2表现突出，尤其553～585 nm和691～755 nm中更明显，其R2均0.4(P<0.001)以上，其中DSR572(波长572 nm处的DSR)的拟合优度(R2=0.5821)最大(图5: A)。从DSR572拟合的POPD变化来看(图5: B)，多项式曲线拟合优度明显优于线性拟合，其R2提高了0.5151，说明DSR与POPD之间有着密切的多项式曲线关系。这是因为553～585 nm位于黄边波段，处于绿峰与红谷之间，其变化被绿峰和红谷波段反射率同时控制。当落叶松健康时，其针叶生化组分含量较高(如色素、水分、氮等)，绿峰和红谷波段光谱反射率差异最大，黄边的坡最陡；随POPD的增加，落叶松受害程度加剧，致使失叶量逐渐增加，针叶叶绿素含量、水分含量和氮含量等下降，使得绿峰和红谷波段反射率随之呈现下降和上升趋势，导致黄边的坡呈逐步平缓态势；等到落叶松重度受害阶段时，POPD从峰值开始下降，失叶量剧增，针叶生化组分含量直降，绿峰和红谷的波峰和波谷减退，甚至有消失现象，使得黄边最终变得平缓。691～755 nm为红边波段，处于红谷与近红外反射肩之间，其光谱反射率变化与针叶叶绿素含量和针叶细胞结构等因素有密切关系。当落叶松健康时，其针叶叶绿素含量较高，红谷和近红外反射肩波段反射率分别达最低和最高值，使红边坡度最大；当POPD逐渐上升时，落叶松受害程度增加，导致失叶量逐渐上升，针叶叶绿素含量下降，针叶细胞结构等发生变化，使得红谷深度减小和近红外反射肩降低，导致红边坡度呈减小态势；等POPD增加到峰值又开始下降时，针叶大量损失，叶绿素含量直降，受害落叶松冠层纹理发生变化，红谷与近红外反射肩波段反射率差异减小，使得红边呈平缓态势。上述波段光谱反射率变化会使微分光谱反射率具有明显的响应。因此，DSR可作为估算POPD的重要指标。

图5 DSR与POPD的线性和多项式曲线拟合效果

(2)光谱指数：为揭示MSI对POPD的敏感性，计算MSI与POPD的多项式曲线拟合优度，并与改进前光谱指数(OSI)的拟合优度进行比较，同时与其线性拟合结果做了对比(图6)。结果显示，MSI与OSI相比，其多项式曲线拟合的R2有提高，说明通过改进处理增强了光谱指数的敏感性。MSI的两种拟合相比，多项式曲线拟合效果有明显提升。具体表现为：MSI多项式曲线拟合效果最佳为TVI(R2=0.5386)，其次为TCARI(R2=0.5034)和MCARI(R2=0.5019)，这比其线性拟合的R2分别提高了0.4479, 0.4015和0.4001，比OSI多项式曲线拟合的R2分别提高了0.0567, 0.1557和0.0886(图6: A)。从TVI拟合的POPD变化来讲(图6: B)，多项式曲线拟合优度明显优于线性拟合优度，其曲线轨迹更接近于POPD的变化趋势。这是因为TVI, TCARI和MCARI的波长分布在绿峰、红谷和红边波段内，对针叶生化组分含量变化的敏感性较强，从而随落叶松POPD呈倒“U”型曲线变化，林木受害程度增加，导致失叶量上升，针叶叶绿素含量逐渐下降，绿峰波段光谱反射率亦呈下降趋势，而红谷波段光谱反射率呈上升态势，使得这些MSI与POPD产生密切非线性关系。上述3种MSI的红边波段的波长调整较明显，如TVI：750 nm→730 nm; TCARI和MCARI：700 nm→680 nm。显然，730 nm和680 nm与调整前波长相比更接近红边位置和红谷位置，提高针叶叶绿素含量变化的感知能力，增强了TVI, TCARI和MCARI对POPD的敏感性。可见，MSI比OSI敏感，更合适作为估算POPD的指标。

图6 OSI, MSI与POPD的线性和多项式曲线拟合效果

2.2 敏感高光谱特征提取结果与分析

基于多项式拟合优度的敏感高光谱特征提取结果显示，DSR和MSI的最敏感特征分别为DSR572和TVI，其拟合的R2分别为0.5821和0.5386(P<0.001)。DSR572处于黄边波段的针叶光谱反射率降低较快的波长(572 nm)上，受绿峰和红谷波段反射率变化的控制，进而与针叶生化组分有密切联系。它可描述针叶色素状态，能够感知落叶松健康状况。TVI是通过530, 650和730 nm等波长反射率构建的指数，受绿峰、红谷和红边变化的控制，与针叶生化组分有密切相关，是表征针叶叶绿素的重要指标。因此，利用DSR572和TVI估算雅氏落叶松尺蠖POPD是合理的。

从DSR和MSI敏感特征组选取结果显示，DSR敏感特征组为DSR572, DSR593, DSR617, DSR738, DSR1210, DSR1323和DSR1487；MSI敏感特征组为TVI, RGI, NDWI, OSAVI2和GM2。显然，基于SPA算法选取的DSR和MSI敏感特征组指标充分地捕捉了雅氏落叶松尺蠖危害下落叶松针叶光谱反射和吸收特征。DSR敏感特征组指标分布在黄边(572, 593和617 nm)、红边(738 nm)、近红外反射肩(1 210和1 323 nm)和近红外吸收谷(1 487 nm)等波段，其中黄边和红边波段分别位于红谷左侧和右侧，能够反映针叶叶绿素吸收特征(Frazieretal., 2014)；近红外反射肩主要受植被细胞结构和水分的影响，对针叶水分吸收特征和细胞受损引起的反射特征较敏感(Zhangetal., 2018)；近红外吸收谷主要受植被水分的影响，能够描述受害过程针叶水分吸收特征(Xiangetal., 2011)。MSI的敏感特征组指标分布在绿峰、红谷、红边和近红外反射肩等波段，其中TVI(530, 650和730 nm), RGI(570和670 nm)和GM2(718和730 nm)能够表征植被色素状态，对受害过程针叶叶绿素吸收特征较敏感(Sankaran and Ehsani, 2011)；OSAVI2(650和780 nm)可削减土壤背景干扰，能够反映针叶叶绿素吸收特征(Stagnarietal., 2018)；NDWI(840和1 260 nm)可描述植被叶片水分的状态，能够表征针叶水分吸收特征(Murphyetal., 2019)。随雅氏落叶松尺蠖POPD变化，受害落叶松针叶色素、水分和细胞结构发生变化，而DSR和MSI敏感特征组指标可以捕捉到针叶的这些变化特征。因此，通过SPA算法提取的上述两组高光谱特征，估算雅氏落叶松尺蠖POPD是可行的。

2.3 基于高光谱特征的POPD估算精度分析

通过DSR和MSI的单高光谱特征(最敏感特征Uni-DSR/Uni-MSI)和多高光谱特征(敏感特征组Mul-DSR/Mul-MSI)，运用多项式回归(PR)和支持向量机回归(SVMR)算法，建立了雅氏落叶松尺蠖POPD估算模型，并对其进行了精度比较(表2)。

表2 POPD估算模型性能评价

为更加直观地揭示POPD估算模型性能优劣，采用模型训练集和验证集的估算值与实测值之间分别进行1∶1直线拟合(图7)。由图可知，在单高光谱特征的估算模型中，Uni-MSI-SVMR模型训练集与验证集的数据分布点相较其他估算模型更加均匀分布在1∶1直线周围，说明Uni-MSI-SVMR比其他模型性能相对较好。在多高光谱特征的估算模型中，Mul-DSR-SVMR模型训练集与验证集的数据分布点相较其他估算模型更均匀并集中分布在1∶1直线两侧，说明Mul-DSR-SVMR比其他模型精度较优且稳定性更好。SVMR与PR模型相比，其训练集与验证集的数据分布点表现出较好的拟合效果，而PR模型的数据点离散程度较为明显，说明高光谱特征结合SVMR算法能够构建性能较高的POPD估算模型。这是因为在具备敏感高光谱特征的前提下，SVMR能够避免过学习或欠学习的问题，其泛化能力较强，可应用于雅氏落叶松尺蠖POPD估算。

图7 训练集和验证集的估算值与实测值1∶1直线拟合

总之，地面非成像高光谱观测试验表明，微分光谱和光谱指数不仅在单株尺度上能够无损、快速监测雅氏落叶松尺蠖虫口密度，而且为后续大区域尺度上虫口密度监测提供了理论 实验依据。外界和背景的干扰会影响数据的客观性，而微分光谱和光谱指数可降低大气和背景干扰影响，在大面积虫口密度估算中也具有很大潜力。这为森林虫害遥感监测指出了新方向，将在今后的研究中对其进行探究和考证。

3 讨论

本研究从雅氏落叶松尺蠖暴发典型区选取不同程度受害的110株落叶松，通过调查与计算样本树害虫POPD和实测林木冠层光谱反射率，分析高光谱特征(微分光谱反射率和高光谱指数)对POPD的敏感性，借助SPA算法提取敏感高光谱特征，结合多项式回归(PR)和支持向量机回归(SVMR)算法，构建POPD估算模型并揭示了高光谱特征的估算潜力，得出以下结果和结论：

(1)多项式曲线拟合相较线性拟合，其更能挖掘光谱指数和微分光谱反射率(DSR)对POPD的敏感性。DSR的敏感性明显的波段主要在黄边和红边波段内，其中在572 nm处的敏感性最显著，DSR572拟合的R2达0.5821(P<0.001)。30种改进型光谱指数(MSI)中最敏感的指数为TVI，其R2达0.5386(P<0.001)，其次为TCARI(R2=0.5034)和MCARI(R2=0.5019)，这比原始光谱指数(OSI)拟合的R2分别提高了0.0567, 0.1557和0.0886。从而MSI和DSR可成为估算POPD的敏感指标。

(2)通过SPA算法提取的DSR敏感指标(DSR572, DSR593, DSR617, DSR738, DSR1210, DSR1323和DSR1487)和MSI敏感指数(TVI, RGI, NDWI, OSAVI2和GM2)充分捕捉了雅氏落叶松尺蠖危害下落叶松针叶叶绿素吸收特征、水分吸收特征以及针叶细胞受损引起的反射特征。可见，SPA是适用于POPD敏感光谱特征提取的一种有效方法。

(4)多个比单个高光谱特征的POPD估算精度更高，并且其SVMR模型性能始终优于PR模型。SVMR泛化能力较强能够避免过学习或欠学习的问题，可应用于雅氏落叶松尺蠖POPD估算。