松材线虫病害扩散预测模型的设计

姚 缘

上海市嘉定区林业站，上海 201822

借助计算机技术，林业工作人员尝试构建松材线虫病害扩散预测模型。一方面，根据松材线虫病害聚集分布特点，建立基于区域的病虫害扩散预测模型。与传统的数值测算方法相比，基于区域的预测模型对病虫害扩散的推演全部在二维平面上完成，通过对样条曲线进行拟合，得到病虫害区域参数，并结合能量函数预估一定时间内的病虫害扩散区域。基于森林病虫害影响力，使用回归方法预测病虫害在二维平面上的传播方向与距离，经过参数校正之后得到病虫害扩散预测轮廓。与元胞自动机方法相比，基于区域病虫害扩散预测模式能够更为准确地预估病虫害扩散机制，能够以较为直观的方式帮助林业从业者了解松材线虫病害扩散情况。

另一方面，根据松材线虫病害高传染率的特点，创建基于传染病动力学的病虫害扩散预测模型。该模型通过对比自然状态与施药状态下，松材线虫病害的扩散差异，以动态化的方式展示SEI模型（自然状态下松材线虫病害扩散模型）与SEIR（施药状态下松材线虫病害扩散模型）模拟扩散情况，为松材线虫病害防治提供数据参考。以下针对这2种病害扩散预测模型进行详细分析。

1 基于区域的病虫害扩散预测模型

1.1 图像拟合

构建预测模型之后，相关工作人员需要使用专业的信息提取设备获得目标区域的矢量数据，通过分批处理病虫害扩散区域图像对矢量图像进行拟合，进而得到病虫害区域坐标点集合。该过程主要分为3个步骤：（1）调节图像尺寸，根据具体的设计需求，将此次模型设计中的图像尺寸调节为W=500，H=500；（2）将矢量图像转变为二值图像，其中深色的区域代表未受到松材线虫病害影响，浅色区域代表已经受到了松材线虫病害影响；（3）提取二值图像上的病虫害扩散区域轮廓点坐标[1]。

提取轮廓点坐标之后，采用均匀间隔提取模式，选取轮廓指标点并将其拟合为曲线型值点，并根据实际需要，调整设计型值点的精度。需要注意的是，闭合曲线首尾需要重合，因此，设计人员假定病虫害轮廓提取点的数量为Ncdp-1个(此处的Ndtp-1数值为100)，则共需要提取Ndtp个型值点。为了保障采样点在图像呈均匀分布状态，工作人员要确定每一个型值点在二维平面上的间隔，其公式如下：

上述2个公式中，floor表示向下取整函数，Ncdp代表坐标点的数量，Nitv则表示型值点取值间隔，Nitv则表示双精度浮点型。先利用浮点数确定二维图像中的坐标值，在确定坐标序数的情况下对大型值数据进行四舍五入，进而得到间隔数据（图1）。

图1 Nitv与Nitv提取型值点的区别

分析图1可以发现，2种数值提取方式在轮廓点末端存在巨大的差异。第2种数值提取方式中，目标区域（Ω区）提取的型值点比第1种方式提取的型值点更加均匀。因此，此次模型设计决定使用第22种型值提取方式。

1.2 病虫害扩散模拟

1.2.1 确定病虫害扩散因素 在正式开始模拟之前，需要确定影响松材线虫病害扩散的主要因素。经过研究发现，一方面，病虫害的扩散与林地的坡度呈正相关关系，即坡度越大，病虫害发病概率越高；另一方面，林地的海拔与病虫害扩散呈负相关关系，即林地的海拔越高，松材线虫病害的发病概率越低。基于这种特点，相关工作人员将上述2种因素作为影响松材线虫病害扩散的特征。

1.2.2 计算病虫害无差别扩散轮廓 假设某林区发生松材线虫病害，设计人员将不受海拔与坡度影响的病虫害扩散轮廓与受海拔与坡度影响的病虫害扩散轮廓进行对比。在对比过程中，相关工作人员将只由于松材线虫繁殖而产生的病虫害扩散称为无差别扩散。在计算无差别扩散时，工作人员利用3次B样条闭合曲线面积的拟合，得到无差别扩散轮廓，利用面积值计算公式得到无差别扩散面积与质心数据，其公式如下：

式（3）中，m00代表封闭3次B样条曲线零阶距，即扩散面积。xi与yi为B样条曲线控制点的横纵坐标。在得到面积之后，利用1阶距m10，m01计算质心，其公式为：

式（4）中，xc表示型心的横坐标，yc表示型心的纵坐标。

在确定病虫害扩散轮廓面积和质心之后，可以开始对无差别扩散情况进行模拟。先计算2次病虫害扩散范围的面积倍数，假设St为在t时刻病虫害的面积，则t时刻与t+1时刻病虫害面积的倍数为通过计算得到t+1时刻病虫害无差别扩散轮廓Pnui(t+1)，其计算公式为：

式（5）、（6）中，Pnui(t+1)（0）与Pnui(t+1)（1）表示在t+1时扩散轮廓的横坐标与纵坐标，而Pt（0）与Pt（1）表示在t时病虫害实际扩散的横坐标与纵坐标[2]。

1.3 预测模拟计算

假设t+1时刻无差别扩散轮廓的型值 点 为Pnui(t+1):(x0,x1,x2,…,xNdtp－1)，考 虑坡度和海拔因素的病虫害扩散轮廓型值点为Pt+1:(x0,x1,x2,…,xNdtp－1)，对比二者之间的差异。

在确定Emin数值的基础上，病虫害扩散轮扣型值点之间最小距离Dist与对应点角度Angle可视为模型数据集中的Y值，计算对应点连接后形成的直线与水平方向射线之间的夹角，计算公式如下：

在得到θ的具体数据之后，工作人员就能够预测t+1时刻的病虫害扩散轮廓（表1）。

表1 松材线虫病害扩散预测流程

2 基于传染病动力学松材线虫病虫害扩散的预测模型

构建基于区域的松材线虫病虫害扩散预测模型，能够使相关工作人员以较为直观的方式观测病虫害扩散情况，但该模型无法准确描述病虫害扩散流行规律。针对这一问题，工作人员尝试利用传染病动力学理论构建SEI模型（自然状态下松材线虫病害扩散模型）和SEIR模型（施药状态下松材线虫病害扩散模型）。通过对比，帮助工作人员了解松材线虫病害扩散动态化过程。

2.1 建立模型

作为一种外来入侵型病原生物，松材线虫对林业的影响巨大，这种生物将松褐天牛作为主要宿主，其繁殖过程会堵塞松树内部的导管，导致松树枯萎。在正式构建模型之前，工作人员假定松褐天牛在林区之中均匀分布，并且松材线虫仅在病树中存在。假设松材线虫的传播方式与人类传染病的传播方式相同，分别设计SEI与SEIR模型[3]。

2.1.1 SEI模型 如果林业管理者并未发现林区中存在松材线虫，则在自然状态下，只需要五年左右的时间，松材线虫就能摧毁一整片松林。由于松材线虫为外来入侵物种，因此，中国本土松树并不具备针对松材线虫的抵抗能力。一旦松树被感染，如果不进行人工干预，松树几乎没有可能自行恢复健康。基于这种情况，在设计模型时，不考虑非人工干预状态下病树转变为移出者R的可能。

当松树枯萎后，松材线虫还可以继续在松树内存活，以松树内的真菌为食。同时，松材线虫的主要传播媒介松褐天牛习惯在枯萎的松树中产卵繁殖，因此，在设计SEI模型时仍然将枯萎的松树视为感染者I。

此外，松材线虫病有着很长的潜伏期，松树在发病初期其表面并无明显变化。通常情况下，松褐天牛并不会从健康的松树中飞出，因此可以将处于染病初期的松树视作潜伏者E。

基于上述背景，工作人员将自然状态下的松材线虫病害扩散模型划分为3个舱室，再利用SEI模型在不同的舱室中进行模拟（图2）。

图2 SEI模型

SEI模型中，S表示健康但不具备免疫力的松树；E表示处于染病初期但不具备传染性的松树；I表示已经受到感染的松树，其微分方式模型如下：

该公式中，系数β代表染病松树将松材线虫传染给其他松树的有效扩散率，其中，潜伏者E有一定的概率（μ）表现出病症，松材线虫病害的潜伏期即1/μ。

2.1.2 SEIR模型 目前，松材线虫防治方式主要包括树干注药、砍伐疫木以及诱捕器诱杀等方式，横向对比之下，树干注药的防治效果较为优秀。因此，此次设计以树干注药为基础建立SEIR模型。

在松树的树干中注入药液，利用松树的蒸腾作用将药液输送至松树的各个部位。需要注意的是，由于健康的松树内部含有松脂，不具备脂溶性特点的药液无法融入松脂，使得药液无法被输送至松树的每一个部位[4]。因此，健康的松树有一定的概率成为移出者R。

一些松树年代久远，具有很高的观赏价值和研究价值，因此相关工作人员研制出用于急救的药液，但此类药品的价格较高，无法普遍使用。因此，在建立SEIR模型时，排除松树因得到了急救药而转好的可能性，其微分方程模型如下：

该公式中，系数σ表示健康松树在注射药液之后获得免疫能力并转变为移出者R的概率。

2.2 可视化设计

为了能够让相关工作人员以较为直观的方式观察松材线虫病扩散过程，设计人员对上述模型进行可视化升级。利用计算机技术设计SEI与SEIR模型可视化界面（图3）。

图3 SEI与SEIR模型可视化界面

分析图3可以发现，此次模型设计使用的参数中，C（林区松树总数量）为29 753，点击“All Green”按钮，UI页面实时反馈的深色部分为健康松树分布位置，点击“SEIimulating”和“SEIRimulating”按钮，可以发现UI界面中出现黄色区域，代表受到松材线虫病害影响的松树分布范围。点击“Intitalization”按钮使得UI页面恢复初始状态[5]。

此时，工作人员可以根据实际需求设定迭代间隔时间，并点击橙色的“开始按钮”，UI页面会按照设定好的迭代时间，展示健康松树与受到松材线虫病害影响的松树分布变化情况，通过这种方式帮助工作人员了解该林业松材线虫病害扩散情况。

3 结束语

松材线虫病害对中国林业资源开发工作产生了很大的影响，如何解决松材线虫病害问题已经成为摆在林业从业者面前的一个难题。想要提高病虫害防治效率，就要对松材线虫病害扩散情况进行精准预测，并以此为基础，制定切实可行的预防措施。因此，相关工作人员尝试基于病害区域和传染病动力学建立松材线虫病害扩散预测模型，并对该模型进行可视化设计，方便林业从业人员对松材线虫病害的扩散进行动态观察，为更好地保护林业资源提供技术支撑。