马云强 杨 斌 赵 宁 张新民

（1.西南林业大学生物多样性保护学院 云南省森林灾害预警与控制重点实验室，云南 昆明 650233；2.西南林业大学生命科学学院，云南 昆明 650233）

云南松（Pinus yunnanensis），常绿乔木，在云南省分布较广，在云南的亚热带高原，从南到北，从东到西，海拔700～3 200 m 都有大面积分布，约占云南林地面积的1/5[1]。云南松分布区为亚热带高原季风类型，属于中亚热带热量水平，年平均气温在13～19℃之间。作为西南地区的乡土树种和荒山绿化造林先锋树种，云南松具有重要的经济价值和生态服务功能，在保持水土、涵养水源、改善环境等方面起着重要作用。云南切梢小蠹（Tomicus yunnanensis）是云南松的重要虫害，可通过柱梢和柱干快速侵染活立木，使寄主大批死亡。自1980 年以来，云南省已经暴发了两次云南切梢小蠹的灾害，累计危害云南松林面积高达百万公顷，使云南松成片林木枯死的面积约9.3 万hm2，虫害暴发后很难控制，严重威胁生态安全，对经济造成巨大损失[2−3]。因此，对云南松受云南切梢小蠹危害进行精细化监测，及时对危害程度进行准确的判断，做到“早发现、早防治”是非常有必要的。

高光谱、激光雷达（LiDAR）技术已经成为农业、林业、草原病虫害监测的重要手段。高光谱数据包含丰富的病虫害危害光谱信息，但空间立体结构信息不足；LiDAR 数据空间分辨率高、立体结构信息丰富，但病虫害危害特征信息不明显。利用机载LiDAR 的数据可以实现单株冠层的精准划分，解决使用高光谱进行病虫害危害程度的提取时出现的“同物异谱、同谱异物”问题，提高分类诊断精度。Anderson 等[4]研究了波形LiDAR 数据与高光谱数据同时应用森林调查中，2 种数据融合会比单独使用一类数据提高精度约5%～8%。Jones 等[5]对高光谱和机载LiDAR 的立体空间融合进行了研究，对两类数据提取的冠层、高程值进行像素级融合，从而提高分类精度，最高可提高11%。赵旦[6]研究了单木提取和高光谱融合进行森林分类的方法，实现了复杂林分下的森林分类。针对云南切梢小蠹的研究多集中于生理生化等方面，而对于利用遥感数据监测云南松受云南切梢小蠹的危害程度分级研究多使用多光谱数据或单一的卫星遥感数据，因此，结合高光谱和机载LiDAR 数据对云南松受云南切梢小蠹危害程度进行分类诊断是非常有意义的。

本研究将以无人机采集的石林县黑龙潭片区云南松林高光谱和LiDAR 数据为基础，通过调查120 株云南松的冠层枯稍率，分析样本冠层枯稍率与光谱特征，通过Matlab 软件，构建分类诊断模型，基于高光谱和机载LiDAR 数据融合提取的云南松单株冠层光谱，对研究区内云南松受云南切梢小蠹危害程度进行分类诊断。以期为云南松受云南切梢小蠹危害的监测、预警和防治提供理论依据和参考。

1 研究区概况

研究区位于云南省昆明市石林县城东部的黑龙潭水库旁（图1），地处东经103°19′53.454″～103°20′17.780″，北纬24°45′58.536″～24°46′19.000″。石林县是林业有害生物国家级中心测报点，云南切梢小蠹是主要测报对象之一，全县乔木林面积为54 864.2 hm2，其中云南松作为该县的优势树种之一，面积约为26 049.5 hm2，占全县乔木面积的47.48%，年均云南切梢小蠹发生面积上千公顷[7]。在石林县黑龙潭水库林区的研究区内，通过查阅石林县森防站记录数据，选取了以云南切梢小蠹危害为主的林区作为研究区域，面积约23.33 hm2。

图1 研究区概况图Fig.1 The map of study area

2 材料与方法

2.1 数据来源及处理

Liu 等[8]对云南松地面全树枯梢率和冠层枯梢率进行了相关性分析发现两者之间相关系数达到0.995，呈极显著正相关。构建一元线性回归方程，利用检验样本检验，方程的预测精度为0.992，均方根误差仅有0.033，证明了通过冠层枯梢率反演地面全树枯梢率，进而判断云南松单株危害程度是可行的。因此，本研究基于冠层对云南松受云南切梢小蠹危害程度进行分类诊断。

2.1.1 样本数据

2019 年11 月18 到30 日，根据国家林业局发布的《林业有害生物发生及成灾标准》（LY/T 1681—2006）[9]（表1），对研究区内云南松单株样本枯稍率进行调查，通过人工地面目测数出和记录点的位置、枯梢数量和健康梢数量，采用式（1）计算枝梢被害率。分别在林区采集健康（枝梢被害率10%）、轻度危害（10%≤枝梢被害率≤20%）、中度危害（21%≤枝梢被害率≤50%）、重度危害（枝梢被害率≥51%）的云南松样本各30 个，共计120 个样本。

表1 机载高光谱S185 成像仪详细技术参数Table 1 Detailed technical parameters of airborne hyperspectral S185 imager

2.1.2 高光谱数据

2019 年11 月17 日11：00—14：00 利用DJI M600 Pro 搭载UHD S185 高光谱成像仪进行高光谱数据采集（图2），在无人机起飞前，进行全白和全黑标定（图2a），得到反射率为100%和0%的数据。飞行高度为100 m，平均飞行海拔高度1832 m，飞行速度为5 m/s，重叠率为60%，采集范围为450～950 nm 的125 个波段，空间分辨率为0.1 m，采集区域与机载LiDAR 为同一区域，飞行8 个架次。设备详细参数见（表1）。

图2 机载高光谱采集设备Fig.2 Airborne hyperspectral acquisition equipment

对采集的高光谱数据进行预处理，包括辐射定标、几何校正、地形校正、大气校正四个步骤。辐射定标是厂商提供完成的，几何校正、地形校正和大气校正均在软件ENVI5.5 中完成，校正完成后利用软件Agisoft PhotoScan 进行影像拼接。根据样本的坐标位置在影像上确定每个冠层的位置和大小（图3），每个冠层绘制1 个ROI（图4a），提取整个ROI 的平均反射率作为单株云南松的光谱反射率，并利用Savitzky−Golay 方法对原始光谱反射率进行平滑降噪处理（图4b）。

图3 样本与高光谱影像叠加Fig.3 Sample canopy and image overlay

图4 样本冠层光谱提取Fig.4 Sample canopy spectrum extraction

2.1.3 机载LiDAR 数据

2019 年11 月16 日，利用DJI M600 Pro 搭载蜂鸟无人机LiDAR 系统（图5）对研究区进行了机载LiDAR 数据采集。根据中国林学会发布的《无人机遥感监测异常变色木操作规程》（T/CSF 002—2018）[10]，采集飞行高度为60 m，飞行海拔平均在1 792 m，航带旁向重叠度在60%，影像旁向重叠度要求在30%，航向重叠度60%。飞行了14 个航带，飞行速度为6 m/s，采集飞行面积约23.33 hm2，坐标系统采用WGS84，UTM 3 度带投影，采集的数据包含格式为*.T04 的POS 数据以及格式为*.isf 的原始点云文件，都储存在设备中的存储盘内，其他的数据还包括设备参数以及天线杆臂值用来辅助解算。设备详细参数见（表2）。

表2 Genius LiDAR 系统主要技术参数Table 2 Genius LiDAR system main technical parameters

图5 机载LiDAR 采集设备Fig.5 Airborne LiDAR acquisition equipment

首先运用软件TerraSolid 对点云数据进行分类，然后获得数字高程模型（Digital Elevation Model，DEM)、数字地表模型（Digital Surface Model，DSM），对DSM 和DEM 作差值处理得到CHM，使用基于冠层控制的冠层高度模型优化算法对异常点进行处理，得到优化后的CHM（图6b），最后计算点云数据的高程值与DEM 的对应像元的差值得到高程归一化点云数据（图6c），具体处理流程见（图6a）。

图6 机载LiDAR 数据预处理Fig.6 Airborne LiDA data preprocessing

2.2 高光谱特征提取方法

在高光谱中消除冗余及干扰噪声等信息带来的误差，有效提取特征信息是建立分析模型的关键。本研究在软件MATLAB 2018 中，采用连续小波变换（CWT）对云南松冠层样本的光谱数据进行特征提取，具体流程见图7。CWT 能够对整条光谱信息曲线在连续波长和尺度上进行分解，有利于光谱信息中的精细部分进行解析。已有研究者将CWT 应用于水稻（Oryza sativa）稻瘟病的区分诊断[11]和冬小麦（Triticum aestivum）叶片含水量的回归模型[12]中，取得了很好的效果。

图7 连续小波变换特征提取流程Fig.7 Feature extraction process of CWT

2.3 云南松受云南切梢小蠹危害分类诊断方法

对基于CWT 提取的云南松冠层高光谱特征信息，在软件MATLAB 2018 中，采用BP 神经网络分类，分类诊断出健康、轻度危害、中度危害和重度危害的云南松。BP 神经网络在农业和病虫害高光谱分析中应用较多[13−14]，神经网络的节点激活函数通常采用S 型函数，见式（2）。

式中：p为Sigmoid 参数，用于对激活函数进行调整。

2.4 分类诊断评价

利用混淆矩阵（CM）统计分类结果，M表示一个混淆矩阵，假设样本有n个类别，则M∈Rn×n，M(i,j) 表示第j个类别被分成第i个类别的样本数量。根据总体分类精度（Overall Accuracy，OA）来评价分类结果，见式（3）。

式中：OAi表示类别i的分类精度，hi表示分类正确的数量，ki表示类别包含的样本总数。

3 结果与分析

3.1 冠层云南松受云南切梢小蠹危害光谱分析

云南松冠层受云南切梢小蠹不同危害程度的光谱反射率及一阶导数见图8，光谱反射率可以反映云南松冠层受云南切梢小蠹不同危害程度。在绿光波段（520～570 nm），不同危害程度光谱反射率云南松冠层大致呈现先上升后下降的变化趋势，即“绿峰”，健康、轻度危害、中度危害光谱反射率变化不明显，重度危害光谱反射率最低；在红光波段（640～700 nm），各个危害程度的云南松冠层光谱反射率大致呈现先下降后上升的变化趋势，即“红谷”，健康云南松冠层光谱反射率最低，重度危害云南松冠层“红谷”有消失的迹象，随着危害程度加重，“红谷”处反射率值有上升趋势；在近红外波段（720～940 nm），不同危害程度云南松冠层光谱反射率波动明显，随着云南切梢小蠹危害程度加重，云南松冠层光谱反射率显著降低。

图8 云南松冠层受云南切梢小蠹不同危害程度光谱反射率及光谱一阶导数平均曲线Fig.8 Spectral reflectance and first derivative average curve of P.yunnanensis canopy affected by T.yunnanensis in different degrees

光谱一阶导数可以反映云南松冠层受云南切梢小蠹不同危害程度。在“绿边”（510～540 nm），健康、轻度危害、中度危害、重度危害存在明显波峰，随着危害程度的加重，峰值逐渐降低；在550～630 nm，健康、轻度危害、中度危害有明显的波谷和波峰，重度危害波动不明显；在“红边”（680～750 nm）有明显的波峰，且峰值达到最高，随着危害程度的加重，峰高逐渐降低，同时红边位置向蓝光方向发生轻微偏移；在750～950 nm，各危害程度波动差异不明显，但大致可以看出重度危害的值大于其他危害程度。综上所述，可以通过光谱反射率来确定云南松受云南切梢小蠹危害程度。

3.2 基于CWT 特征提取结果

首先利用小波基函数对原始光谱进行分解，小波基选用Daubechies（DBN）系列中的DB4 进行小波分解和重构的误差最小，所以对于健康、轻度、中度、重度的小波系数计算采用 DB4 作为小波基函数。云南松冠层样本的光谱反射率波段数n=125，分解尺度为m=10。得到小波系数矩阵后，将小波系数矩阵与枯梢率做相关系数检测，发现当分解尺度为8～10 时，得到的相关系数趋于稳定，所以提取最大分解尺度为7，如图9 所示，在1～7 尺度中相关性较高的小波系数较多，不能实现降维的目的。

图9 小波特征分析及相关性分析Fig.9 Characteristics analysis and correlation analysis of CWT

将120 个样本的125×7 个小波系数对应的变量与枯梢率进行多元线性回归和多重共线性分析，经过不断调整模型，调整R2逐渐增大，标准估计误差逐渐减小，在第22 次拟合后，得到了共线性比较弱且对枯梢率相关性比较强的16 个特征（表3）。此时的调整R2达到了0.991，说明云南松不同受害程度的枯梢率与16 个特征的相关性很大。从表中可以看出，这16 个特征的方差膨胀系数VIF 均不大于10，说明它们之间的共线性较弱。此外，t检验的概率P值（显著性Sig.）均小于0.05，说明这些特征对枯梢率来说是显著的。

表3 连续小波特征及共线性统计Table 3 CWT characteristics and collinearity statistics

3.3 基于BP 神经网络训练分类诊断模型

本研究将120 个云南松样本数据的70%作为训练集，30%作为验证集。以提取的16 个小波特征作为BP 神经网络分类算法的训练集，以正切Sigmoid 函数作为BP 神经网络算法输入层和隐含层的传输函数，输出层则采用线性传输函数，经过多次训练，最终确定，当隐含层为8 层时，训练分类精度达到最高，为94.05%（表4），验证精度为94.44%（表5）。

表4 基于BP 神经网络的分类训练结果Table 4 Classification training results based on BP neural network

表5 基于BP 神经网络的分类验证结果Table 5 Classification verification results based on BP neural network

3.4 基于机载LiDAR 识别与提取云南松单株树冠

单株分割是提取单株树冠结构的关键，首先利用局部最大值法对高程归一化点云数据进行探测，将探测到的单株的树顶点，即高程最大值点，作为单株的位置，然后采用归一化切割（Ncut）算法[15]，使用全局最大值替代局部最大值为先验条件，并作规则约束，实现单株的精确探测（图10），基于NCut（M，N）算法的分类结果见（图11）。

图10 树顶探测结果Fig.10 Detection results of treetop

图11 基于Ncut 的单株分割结果Fig.11 Segmentation results of single wood based on Ncut

3.5 基于高光谱和机载LiDAR 的单株云南松危害诊断

3.5.1 云南松单株区域高光谱提取

利用单株的分类结果提取出矢量数据，单株矢量数据包含了单株的位置和冠幅信息，利用ENVI5.5，将单株冠幅分类结果进行ROI 的自动生成，进行单株冠层矢量与高光谱信息的关联，计算单株冠层对应范围内的高光谱像元（图12）。高光谱采集的分辨率为0.2 m，每个单株冠幅内包含了超过100 个高光谱像元，因此针对每个单株所对应的高光谱像元，取其平均值作为单株对应的高光谱值，最终得到单株云南松冠层矢量数据，其中属性包含冠层范围内所有波段的高光谱均值。

图12 云南松单株区域高光谱提取Fig.12 Regional hyperspectral extraction of P.yunnanensis

3.5.2 云南松受云南切梢小蠹危害程度分类可视化

利用训练好的BP 神经网络分类诊断方法对研究区内23 hm2林区进行单株云南松分类，共提取出11 029 株云南松，其中健康10 142 株、轻度危害490 株、中度危害266 株、重度危害131株。选用采样的120 株样本对分类结果进行精度验证，如表6 所示。从表6 可以看出，重度危害的分类精度最高为100%，其次是健康、中度的分类精度为93.33% 和90.00%，轻度危害的分类精度最低为80.00%，总体分类精度为90.83%。

表6 基于CWT 的BP 神经网络算法分类结果及验证Table 6 The classification results and verification of BP neural network algorithm based on CWT

将所有分类结果通过唯一编号与单株云南松冠层矢量数据相关联。利用ArcGIS Pro 以分类结果字段为渲染字段，进行分类结果可视化（图13）。

图13 云南松受云南切梢小蠹危害程度分类可视化Fig.13 Classification and mapping of damage degree of T.yunnanensis to P.yunnanensis

4 结论与讨论

本研究在云南松切梢小蠹蛀害区域，采用无人机进行大面积的高光谱数据和LiDAR 数据的采集。以冠层枯梢率划分云南松受云南切梢小蠹危害程度为基础，分析云南松不同受害程度的冠层光谱特征并基于小波变换提取光谱特征，采用BP 神经网络分类算法训练单株云南松受云南切梢小蠹危害程度分类诊断模型。利用机载LiDAR 数据提取DEM、DSM、CHM 和高程归一化点云数据，采用局部最大值方法探测云南松单株树顶，通过NCut 分割方法进行云南松单株树冠提取，将提取的云南松单株冠层矢量数据与高光谱数据融合，获得整片林区云南松每株冠层高光谱数据。基于单株云南松冠层光谱数据，利用分类诊断模型对研究区内云南松受云南切梢小蠹危害程度进行分类诊断，一共提取出11 029 株云南松，其中健康10 142 株、轻度危害490 株、中度危害266 株、重度危害131 株。利用120 个样本对分类结果进行验证，得到总体分类精度为90.83%，最后将分类诊断结果进行可视化，可清晰的看到研究区内不同受害程度的云南松分布，实现了云南松受云南切梢小蠹危害程度的精细化监测，可为大面积的云南切梢小蠹实时精准防治提供依据。

植被光谱反射率及一阶导数特征可以较好地反映植被受病虫害危害情况，郭伟等[16]研究冬小麦全蚀病发现，不同危害等级的冬小麦光谱反射率在“绿边”出现“波峰”，在“红边”出现“波谷”特征，在近红外波段，随着冬小麦病害程度的加深，冬小麦光谱反射率显著降低。连玲等[17]在研究枣截形叶螨危害冠层光谱特征中发现，冠层光谱反射率有“绿峰”和“红谷”，在近红外波段，随着叶螨危害的加重，光谱反射率渐渐降低；光谱一阶曲线在“红边”存在“蓝移”现象。本研究发现，受云南切梢小蠹危害的云南松，在近红外波段，光谱反射率和光谱一阶曲线随着危害程度的加重逐渐降低，一阶微分曲线存在“红边”“蓝移”现象。这说明光谱反射率和光谱一阶曲线可以很好地反映云南松受云南切梢小蠹不同危害程度的情况。在森林单株提取研究中，赵旦[6]与王濮等[18]对比标记控制分水岭和Ncut分割方法发现，Ncut 的提取精度较高，本研究利用Ncut 分割方法进行云南松单株树冠提取，分割提取效果较好。在不同受害程度光谱特征提取及分类诊断方法研究中，本研究基于CWT 的BP 神经网络分类方法精度达到了90.83%，但是本研究的研究区是以云南切梢小蠹危害为主的林区，没有涉及到复杂林分的单株提取，缺少验证。因此，在下一步的研究中，需要通过高光谱及其他生理特征，研究枯梢的生理性和病理性的危害的区别，以达到更全面的监测和分类诊断。