Study on early diagnosis model of pine wilt disease in infected forest area of southern Sichuan

在自贡市自流井区荣边镇和飞龙峡镇交界区域，即E：104°38′08″~104°39′45″，N：29°14′42″~29°16′13″之间，平均海拔440m。植被包括马尾松(Pinus massoniana Lamb.)、香樟(Cinnamomum camphora (L.) Presl)、大桉(Eucalyptus grandis Hill ex Maiden)、杉木(Cunninghamia lanceolata (Lamb.) Hook.)、麻栎(Quercus acutissima Carr.)、乌桕(Triadica sebifera (Linnaeus) Small)、慈竹(Bambusa emeiensis L. C. Chia & H. L. Fung)、毛竹(Phyllostachys edulis (Carriere) J. Houzeau)、油茶(Camellia oleifera Abel.)等。自贡地处中亚热带湿润季风气候，平均气温17.9℃，极端最高气温40℃，最低气温−2.8℃。平均相对湿度77%，夏、秋季节相对湿度超过80%。紫色土类占土壤总面积的50.08%，黄壤土类占13.73%，黑褐土类占0.71%。

采用北京安洲科技自主研发高光谱成像系统。包括：SOC710VP便携式高光谱成像仪、便携式电脑（采集软件SpecView）和98%反射率的硫酸钡白板。高光谱成像仪数据参数参考表1。

利用松墨天牛作为传播媒介的作用^[21]，5月中旬对实验样株进行预处理。选取疫木伐桩500 m范围内外观表现正常的马尾松，将松墨天牛高效持久型诱芯（呈康APF-1）挂于树干上。6月上旬，从预处理的样株中选取实验样株，选择标准：树体上有松墨天牛产卵刻槽，或用刀片剥离树皮，发现树体流胶明显减少。结合植株胸径和立地条件，共选择9株实验样株。此外，在远离实验样株且500 m范围内未发现疫木的孤立林班内选择健康马尾松作空白对照，实验样株分布见表2。

结合松墨天牛在我省的生物学特性（4月中下旬开始羽化，5月中旬为羽化高峰期），从6月中旬开始，对实验样株进行光谱影像采集。在天气晴朗光线充足的条件下，利用高枝剪从样株四个方位随机采摘松针，松针采集量控制在100枚以上，将松针表面简单清洁后架设光谱仪进行采集，影像采集前用硫酸钡白板进行校正。测量过程中，避免人员在目标区两侧走动。每隔一段时间进行一次光谱采集，直至马尾松表现出明显的受害症状为止。每个样品重复进行3次影像采集，取其中成像质量最高的一组进行研究。

选用ENVI 5.3遥感图像处理软件，通过生成增益系数和反射率转换，提取目标物光谱曲线^[22]。

增益系数生成：将采集的高光谱影像raw文件导入ENVI，针对硫酸钡白板范围绘制ROI，提取并导出平均DN值的txt文件，通过白板反射率计算出增益系数，提取波段对应的增益系数另存为txt文件备用。增益系数计算公式：增益系数=白板反射率/白板平均DN值/100。

反射率转换：将高光谱影像raw文件导入ENVI Classic，在Edit Header中编辑调用备用增益系数，通过工具箱中的预处理对增益系数进行应用，另存文件并加载其真彩色影像以查看光谱反射曲线。

通过ENVI打开保存文件，将可选择范围移动至目标物，对光谱曲线进行平滑处理并编辑图片，导出波段及对应反射率的txt数据文件。

获取不同波段及反射率数据，利用相应模型进行各种变换，提取能够反映马尾松感染松材线虫病不同感病阶段的敏感波段或波段指数^[8,23]。结合不同监测时间、无重复双因素分析方差值建立相应公式进行比较并提取敏感波段。选择在高光谱数据处理中具有良好效果的植被指数，如差值、比值、归一化和加强指数（见表3），利用不同植被指数及不同敏感波段组合，结合指数、对数、线性、二阶导数和幂指数等在辨别细微光谱差异上的优势，构建感病马尾松林间早期诊断模型^[6]。

自2020年6月中旬开始，结合实验样株反映和天气状况，野外共进行5次（6月19日、7月19日、8月16日、9月8日和9月30日）高光谱影像采集，获取了47组141号高光谱数据。每幅影像光谱介于402.6 nm~1005.5 nm之间，可分为176个波段，每个波段3.7 nm。利用ENVI对获取的松针影像进行处理，获得了平滑的光谱反射曲线，见图1。

Figure 1.  Spectral image processing results of plants in different periods

通过一段时间的野外监测，发现3号样株于第2次（65 d）野外影像采集时出现枯黄现象，6号和7号样株于第3次（90 d）野外影像采集时前后出现明显枯死现象。对照株与感病植株光谱差异显著，主要表现为：随着感病程度深入，绿光区的反射峰逐渐减弱，直至病株在外观上表现出明显枯萎，光谱反射曲线上的反射峰或吸收谷特征消失，曲线斜率显著降低，见图2。

Figure 2.  Spectral reflection curves of infected and control strains

结合感病植株和对照植株光谱反射曲线，选择青光、绿光、红光和近红外光谱反射率变化明显的区域提取相应波段进行差异显著性分析。波段集中在472.6 nm~488.7 nm，550.8 nm~567.3 nm，668.5 nm~685.6 nm，751.6 nm~851.1 nm，其中前三个波段范围各选6个波段，近红外宽平反射范围选8个波段，各波段及其反射率数值见表4。

选择同一时间段对感病株与对照株进行无重复双因素分析，方差分析表明差异显著性。利用公式：

式中，T_N为监测总天数；T_n为具体监测天数；B_n为具体监测时间与监测总天数比值；时间越早，B_n值越大；C为野外监测总次数；F_c某次监测特定波段下光谱反射率方差值大小排序；R_c为某次监测特定波段的光谱反射率方差加权值，R值越大，则该相应波段越敏感。

通过对不同波段范围内不同时间光谱反射率方差进行加权平均，提取整个感病过程中的敏感波段。最终提取4个敏感波段，即488.7 nm、550.8 nm、682.2 nm和779.8 nm，见图3。

Figure 3.  Significant differences between the experimental group and the control group in different bands and different times

选用四种常用的植被指数，结合不同波段进行分析。在敏感波段中，仅779.8 nm波段属于近红外波段，可用于NDVI和EVI计算；488.7 nm和550.8 nm为蓝光波段；682.2 nm为红光波段。基于感病株数据，利用指数、对数、线性、二阶导数和幂指数等建立模型回归拟合，四种植被指数中DVI、RVI、EVI的幂指数效果最佳，DNVI的二阶导数效果最佳。通过不同波段组合对比，EVI以550.8 nm蓝光波段建立的模型回归拟合最佳，其R²为0.9754，不同植被指数回归拟合见表5。

建立早期诊断模型的重点在确保对照株稳定的前提下，显示出与感病株的显著差异。基于EVI建立的感病株与对照株病程反演，通过感病株二次回归拟合得到诊断模型。最佳回归拟合方程为指数型，其R²为0.9754，表明感病株感病后期病程加剧，与对照株形成显著差异；在该拟合中对照株整体稳定，其R²为0.0199，为几种趋势分析中的最佳模型。而选用DVI、RVI、NDVI建立的模型中，对照株稳定性较差，与感病株形成交叉曲线而造成一定干扰，难以对病程中植株感病做出定性判断，早期诊断模型见图4。

Figure 4.  Inversion model of plant susceptibility

通过感病指数K值可以发现：当K>6时，能确定植株感病已经到中后期，感病时间超过90 d；当K<1时，可知感病株处于感病早期（36 d以内），此时通过肉眼无法感知植株针叶颜色和长势的变化；当1<K<6时，此时植株表现存在一定差异，如树龄较小、长势较差的植株会出现明显的外观变化，而树龄较大或长势较好的植株不会出现肉眼可见的变化。

马尾松感染松材线虫病早期诊断模型：

式中，K为马尾松感病指数；EVI为加强植被指数；R_NIR为近红外波段779.8 nm的反射率；R_RED为红光波段682.2 nm反射率；R_BLUE为蓝光波段550.8 nm反射率。

在自然条件下完成了10株马尾松的5次监测，其中3号、6号和7号样株出现明显的死亡症状，采样镜检发现有松材线虫。林缘与林间植株光谱曲线差异不明显，林缘植株受松材线虫侵染的风险更大，幼树染病后表现症状时间更短。不同立地条件，植株光谱曲线有变化，即4号和5号样株虽没出现明显的枯黄现象，但光谱曲线与对照株有明显差异。结合现地进行分析，4号样株位于峭壁上，5号样株位于水渠边，这些因素使植株长势受影响，进而在光谱曲线上反映出差异，但是与感病株相比差异不显著。在提取敏感波段时，加入感病过程中的时间因子，增加早期监测比重。样品在采集过程中，易受周边环境的影响，比如光强、大气、人为活动等，在样本的采集时需规范化操作，注意对样品保鲜保湿，防止失水蒸发，避免后期数据分析的过程中出现误判。

结合已有研究进行讨论：

（1）实验条件：部分研究的实验株为人工接种^[13,24]，自接种日开始计算，时间覆盖自7月至9月。该实验结合松墨天牛生物学特性，利用诱芯引诱其完成松材线虫接种，时间覆盖6月至9月，即在秋季普查之前，尽量模拟林间自然条件。有研究在林间选择3个研究区域，在未知是否会发生松材线虫病的情况下进行^[8]，研究结果出现不同程度的波动，这可能与病害程度以及出现症状到死亡所经历的时间有关，即染病时间不确定。

（2）对不同区域的光谱曲线比较：感病植株整体变化趋势相同，即“绿峰”和近红外宽平反射区降低，红光区“红谷”上升，光谱曲线波动性降低。数据采集仪不同，大部分研究基于ASD高光谱仪^[8,13,14,24]，其分辨率3.5 nm，而SOC710VP光谱分辨率为1.3 nm。此外，华东地区马尾松“红边”反射峰值较小，略高于0.5^[24-25]，使得感病植株的峰值变化差异显著增加。由于不同光谱成像仪或不同区域植株长势的差异，有必要对当地植株进行光谱采集，且造成的反射率差异会对敏感度具有影响。

（3）敏感波段提取：在450~780 nm波长范围内反射光谱的曲线分形维数的变化是植被健康状况发生变化的综合描述，对马尾松松材线虫病害的发生具有一定的指示作用^[8]。700 nm至800 nm为“红边”波段，该区域反射率会显著变化，它是由于叶绿素的强吸收和近红外叶子内部散射导致的高反射所形成的^[26]。我们提取的敏感波段介于其中，一些研究选用的近红外波段为810 nm，实验中发现其敏感性相对较弱。在提取敏感波段过程中，选择不同时间感病株与对照株进行无重复双因素分析，而一些研究直接在可见光和近红外波段中选用某特定波段，没有对波段敏感性进行比较。

（4）反演模型构建：首先，大量研究表明植被的生物化学参数与光谱反射率之间的关系基本上是非线性的，由此选择神经网络等进行拟合回归，但是输入参数选取主要依靠主观经验^[13]，引入多因子进行分析可能造成较大误差，且林间操作繁琐，效率较低。其次，选用何种光谱指数进行建模很重要。部分研究选择RVI^[12]，在感病模型中，对照株与感病株在19 d前形成交汇，65 d至76 d比值降低，之后上升在下降，等于感病植株43d前值^[24]，感病植株稳定性较差，在判读时会造成干扰。此外，分形理论也是单纯对光谱曲线进行分析，但反射光谱曲线局部变异会对整个曲线带来影响，故其选择方差法计算反射光谱曲线的分形维数^[8]。

30~90 d为感病的早中期，不同植株表现存在一定差异。3号和6号样株处于幼龄期，在65 d前后表现出明显的枯死状，7号样株处于近成熟期，在9月后表现肉眼可见的感病状。此时为松材线虫病早期研究的关键期，可结合反演预测模型进行定性判断。但建立的诊断模型对于植株龄期没有进行分析，且研究中涉及的样株数量偏少，表现出症状的样株更少，今后可以利用激光雷达或林业资源数据探究不同龄期的感病诊断研究。此外，GaiaField-Pro系列高光谱成像仪质量轻便，便于无人机携载。当下部分固定翼无人机在饱满续航能力下可以飞行2.5 h，单次作业面积可以达到2600 hm²，利用该早期诊断模型与无人机结合的优势，为今后大面积开展松材线虫病林间早期监测与防治提供了可能。