一类具时滞森林病虫害传染病模型

摘 要：控制森林病虫害有利于增加森林碳汇量，对实现我国“双碳”战略具有重要意义。为了研究森林病虫害对林区的影响，建立了一个含饱和发生率的森林病虫害传播时滞微分方程模型，研究了平衡点的存在性与稳定性，以及平衡点处Hopf分支的存在性，同时利用多尺度方法推导出系统在Hopf分支临界点附近的规范型。结合全国森林清查报告等实际数据，对模型中的参数进行分析与预测，得出各参数的取值范围。选取合理参数进行数值模拟，并研究了部分参数变化对模型的影响。最后根据模拟结果，结合黑龙江省林区森林病虫害的防控现状给出结论和建议。

关键词：森林病虫害；时滞微分方程；稳定性；Hopf分支；多尺度方法；规范型

DOI：10.15938/j.jhust.2024.03.018

中图分类号： O175.1

文献标志码： A

文章编号： 1007-2683（2024）03-0149-10

A Kind of Forest Disease and Pest Infection Model with Time Delay

FEI Xingyan， DING Yuting

（College of Science， Northeast Forestry University， Harbin 150000， China）

Abstract：Controlling forest diseases and insect pests is beneficial to increase forest carbon sink. It is of great significance to realize China′s “double carbon” strategy. In order to study the impact of forest diseases and insect pests on the forest， the authors establish a delayed differential equation model of forest diseases and insect pests infection with saturation incidence， study the existence and stability of the equilibrium and the existence of Hopf bifurcation near the equilibrium. At the same time， the authors use the multiple time scales method to derive the normal form of the system near the critical point of Hopf bifurcation. Based on the data such as the forest inventory report of China， we analyze and predict the parameters of the model， obtain the range of the parameters， select reasonable parameters for numerical simulation， and study the influence of the changes of some parameters. Finally， according to the simulation results， we provide some relevant conclusions and suggestions based on the current situation of forest disease and insect pest prevention and control in forest areas of Heilongjiang Province.

Keywords：forest diseases and insect pests; delayed differential equation; stability; Hopf bifurcation; multiple time scales method; normal form

0 引 言

在气候异常变化的背景下，实现全球“碳中和”已经成为国际共识，我国也提出了“在2060年前实现碳中和”这一战略目标［1］。森林生态系统包含了生物圈中大部分的碳，陆地植被与大气之间有90%以上的碳交换由森林完成，且森林面积及覆盖类型的变化将会引起全球碳循环的变动［2］。然而我国是世界上受森林病虫害影响最为严重的国家之一［3］，以食叶害虫为例，叶片受到虫害后，叶片面积减少导致光合作用减弱，从而影响森林固碳量。有统计表明：每年植食性昆虫对森林生态系统的净初级生产力破坏达到5%［4］。因此，我国森林病虫害问题对实现“双碳”战略目标的影响不容忽视。

目前，学者们对森林病虫害的探讨主要在监测识别、防治对策、损失评估等领域［5-9］。在林区中，虽然树木彼此独立，但是林区中的杂草、枝梗、媒介昆虫、土壤等均可以作为传播媒介，导致病虫害在林区中蔓延。例如，柑橘的黄龙病通过媒介昆虫橘木虱传播［10］，松树的松材线虫病通过媒介昆虫松褐天牛传播［11］……从而有很多学者基于传染病模型建立了微分方程模型刻画林区内部健康树木与染病树木面积的连续变化的情况［12-13］。由于林区面积较大，实际情况中染病树木变化率呈非线性变化，因此以饱和接触率代替线性增长率更具合理性［14-15］。此外，考虑到对染病树木采取防治措施后存在一定的恢复时间，因此也将考虑树木从染病状态恢复至健康状态的时间延迟。

综合以上分析，本文建立一类含饱和发生率的时滞微分方程模型，分析该模型的平衡点的存在性、稳定性，并讨论Hopf分支的产生条件及其约化在中心流形上的规范型。

1 数学建模

在没有外界干预的情况下，林区面积的增长情况符合logistic模型的形式。森林病虫害发生后，会通过媒介昆虫传播至健康树木，导致健康树木面积的扩张速率减小，而染病树木面积加速增长。考虑到一棵染病树木的影响能力有限，此处用饱和发生率来刻画病虫害的感染情况。森林病虫害经过治理后，有一部分树木会被治愈，小部分树木会由于受害过于严重而失去固碳能力。

基于现有的流行病学数学模型［16］，构建如下含饱和发生率的时滞微分方程模型。x（t），y（t）分别表示t时刻该林区健康树木和染病树木的面积。考虑到树木的自然增长与死亡，建立如下模型：

dxdt=rx1－xN－βxy1+αx+my（t－τ）

dydt=βxy1+αx－my（t－τ）－by（1）

其中：r为健康林区面积的自然增长率；N为健康林区的环境容纳量；βxy/（1+αx）为森林病虫害传播的饱和接触率；β为有效接触率；α为饱和参数；m为林区森林病虫害的治理率；b为森林病虫害对该林区树木的致死率；τ为染病树木的恢复时间。假设模型中所有系数均为正数。

在接下来的章节中主要分析模型 （1） 的动力学性质。

2 模型分析

模型 （1） 始终存在无病虫害平衡点E0=（N，0）。当β>α（m+b）时，系统 （1） 存在病虫害平衡点E*=（x*，y*），其中

x*=m+bβ－α（m+b）， y*=rbx*（1－x*N）

下面主要讨论病虫害平衡点E*附近的解的动力学性质。

模型 （1） 在病虫害平衡点E*处的特征方程为

λ2+λ（－r+a+c－m）－m（－r+a+c）+cn+m（λ－r+a）e－λτ=0（2）

其中

a=2rNx*，c=βy*（1+αx*）2，n=βx*1+αx*

当τ=0时，特征方程（2）变为如下形式：

λ2+λ（－r+a+c）+c（n－m）=0（3）

给出如下假设

（H1）：－r+a+c>0，

当假设（H1）成立时，特征方程（3）有两个负实根，模型（1）的病虫害平衡点E*在τ=0时是局部渐近稳定的。

当τ>0时，由于m≠n，则有c（m－n）≠0，因此λ=0不是方程的根，模型（1）在E*处无不动点分支。假设λ=iω为方程（2）的根，代入特征方程 （2）并分离实虚部有

－ω2－［m（－r+a+c）－cn］=

－m（a－r）cosωτ－mωsinωτ

ω（－r+a+c－m）=

－mωcosωτ+m（a－r）sinωτ（4）

根据方程 （4） 可以推导出如下结果：

Qsin（ωτ）=ω（a－r）（－r+a+c－m）+ω3m（a－r）2+mω2+

ω［m（－r+a+c）－cn］m（a－r）2+mω2

Pcos（ωτ）=（a－r）［m（－r+a+c）－cn］m（a－r）2+mω2－

ω2（c－m）m（a－r）2+mω2（5）

设z=ω2，方程 （4） 化为

h（z）=z2+D1z+D2=0（6）

其中

D1=（a+c－r）2－2cn， D2=－cb［2m（a－r）－cb］

当D2<0即m∈cb2（a－r），+∞时，方程 （6） 有一个正根，记为z1，则有h′（z1）>0；当D1<0，D2>0，即m∈0，cb2（a－r），n∈（a+c－r）22c，+∞时，方程 （6） 有两个正根，记为z2<z3，则有h′（z2）<0，h′（z3）>0。

不失一般性，假设方程 （6） 有正根zk（k=1，2，3），则由 （5） 可得

τ（j）k=1ωk（arccosPk+2jπ），Qk≥0

1ωk（2π－arccosPk+2jπ），Qk<0

k=1，2，3，j=0，1，2，…（7）

其中

Qk=sin（ωkτ）=

ωk（a－r）（－r+a+c－m）+ω3k+ωk［m（－r+a+c）－cn］m（a－r）2+mω2k

Pk=cos（ωkτ）=

（a－r）［m（－r+a+c）－cn］－ω2k（c－m）m（a－r）2+mω2k

ωk=zk，k=1，2，3

由方程 （2） 可得Hopf分支横截条件如下：

Redλdτ－1λ=iωk=h′（zk）m2［ω2k+（a－r）2］≠0（8）

定理1 当假设 （H1） 成立，模型 （1） 平衡点E*的稳定性和Hopf分支的存在性有如下结论：

1）当m∈0，cb2（a－r），n∈0，（a+c－r）22c时，对于τ≥0，平衡点E*是局部渐近稳定的，模型 （1） 无Hopf分支。

2）当m∈cb2（a－r），+∞时，方程 （6） 有一个正根z1。当τ∈［0，τ（0）1），平衡点E*是局部渐近稳定的；当τ≥τ（0）1，平衡点E*是不稳定的；当τ=τ（j）1 （j=0，1，2，…），平衡点E*处产生Hopf分支。

3）当m∈0，cb2（a－r），n∈（a+c－r）22c，+∞时，方程 （6） 有两个正根z2和z3（假设z2<z3），有τ（0）2>τ（0）3，则m∈N，使0<τ（0）3<τ（0）2<τ（1）3<τ（1）2<…<τ（m－1）2<τ（m）3<τ（m+1）3。当τ∈［0，τ（0）3）∪∪ml=1（τ（l－1）2，τ（l）3），平衡点E*是局部渐近稳定的；当τ∈∪m－1l=0（τ（l）3，τ（l）2）∪（τ（m）3，+∞），平衡点E*是不稳定的；当τ=τ（j）k （k=2，3. j=0，1，2，…），平衡点E*处产生Hopf分支。

根据参数含义，参数m和n分别可以表示森林病虫害的治理率和传播率。该定理中m和n的取值范围与病虫害平衡点处分支现象的生物学意义相统一，将在4.2节中深入讨论。

3 Hopf分支的规范型

树木受到病虫害干扰或因病虫害致死后，森林的固碳能力下降，而受干扰后森林的更新过程会对森林结构与功能等多个方面产生影响［17］。因此关注染病树木的恢复时间是很必要的，选取染病树木的恢复时间τ作为分支参数。假设模型 （1） 的病虫害平衡点E*在τ=τc=τ（j）k处经历Hopf分支，此时特征方程 （2） 有一对纯虚根λ=±iω，其中τ（j）k由式（7）给出。Hopf分支的研究方法有多种［18-19］，下面利用多尺度方法推导模型 （1） 在病虫害平衡点E*附近的Hopf分支规范型。

首先，令tMT ExtraaA@t/τ，将病虫害平衡点E*转移到原点，即

=x－x*

=y－y*（9）

为了方便表达，仍然用x和y表示和，则方程 （2） 可写为

dxdt=τrx1－x+2x*N－β（x+x*）（y+y*）1+α（x+x*）+

my（t－1）+βx*y*1+αx*

dydt=τβ（x+x*）（y+y*）1+α（x+x*）－my（t－1）－

by－βx*y*1+αx*（10）

将式 （10） 中的分式在病虫害平衡点处进行泰勒展开，得到如下方程：

dxdt=τ［a1x+a2y+a3x2+a4xy+a5x3+a6x2y+

my（t－1）+…］

dydt=τ［b1x+b2y+b3x2－a4xy－a5x3－a6x2y－

my（t－1）+…］（11）

其中

a1=r－a－βy*（1+αx*）2， a2=－βx*1+αx*，

a3=αβy*（1+αx*）3－rN，a4=－β（1+αx*）2，

a5=αβy*（1+αx*）4， a6=αβ（1+αx*）3，

b1=βy*（1+αx*）2， b2=βx*1+αx*－b，

b3=－αβy*（1+αx*）3，

为了方便计算，将方程 （11） 记作

Z·（t）=τN1Z（t）+τN2Z（t－1）+τF（Z（t），Z（t－1））（12）

其中

Z（t）=［x（t），y（t）］T，Z（t－1）=

［x（t－1），y（t－1）］T

N1=a1a2b1b2，N2=0m0－m

设h为方程 （12） 的线性矩阵对应于特征值λ=iωτ的特征向量，h*为方程 （12） 的线性矩阵的伴随矩阵对应于特征值λ=－iωτ的特征向量，满足

〈h*，h〉=h*Th=1（13）

通过计算可得，

h=1，iω－a1a2+me－iωτ

h*=d－b1iω+a1，1T

d=（a1+iω）（a2+meiωτ）－b1（a2+meiωτ）－a21－2iωa1+ω2（14）

假设方程 （12） 有如下形式的解

Z（t）=Z（T0，T1，T2，…）=∑+∞k=1εkZk（T0，T1，T2，…）（15）

其中

Z（T0，T1，T2，…）=［x（T0，T1，T2，…），y（T0，T1，T2，…）］T

Zk（T0，T1，T2，…）=［xk（T0，T1，T2，…），yk（T0，T1，T2，…）］T

由于

ddt=T0+εT1+ε2T2+…=D0+εD1+ε2D2+…（16）

其中Di=Ti表示微分算子，则

Z·（t）=εD0Z1+ε2D1Z1+ε3D2Z1+ε2D0Z2+ε3D1Z2+ε3D0Z3+…（17）

将y（t+1）=y（T0－1，ε（T0－1），ε2（T0－1），…）在y（T0－1，T1，T2，…）处进行泰勒展开，得到

y（t－1）=εy1，τc+ε2y2，τc+ε3y3，τc－ε2D1y1，τc－ε3D2y1，τc－ε3D1y2，τc（18）

其中yj，τc=y（T0－1，T1，T2，…），j=1，2，3，…。

将τ视为分支参数，设τ=τc+εμ，其中τc=τ（j）k为Hopf分支临界值，μ为扰动参数，ε为无量纲尺度参数。将等式（15）～（18）代入方程（12），分别比较等号左右两边ε，ε2，ε3的系数，得到如下表达式：

D0x1－τc［a1x1+a2y1+my1，τc］=0

D0y1－τc［b1x1+b2y1－my1，τc］=0（19）

D0x2－τc［a1x2+a2y2+my2，τc］=

－D1x1+μ［a1x1+a2y1+my1，τc］〗+

τc［a3x21+a4x1y1－mD1y1，τc］

D0y2－τc［b1x2+b2y2－my2，τc］=

－D1y1+μ［b1x1+b2y1－my1，τc］+

τc［b3x21－a4x1y1+mD1y1，τc］（20）

D0x3－τc［a1x3+a2y3+my3，τc］=－D2x1－D1x2+

μ［a1x2+a2y2+my2，τc+a3x21+a4x1y1］+

τc［2a3x1x2+a4（x1y2+x2y1）+a5x31+a6x21y1］－

m［μD1y1，τc+τc（D2y1，τc+D1y2，τc）］

D0y3－τc［b1x3+b2y3－my3，τc］=－D2y1－D1y2+

μ［b1x2+b2y2－my2，τc+b3x21－a4x1y1］+

τc［2b3x1x2－a4（x1y2+x2y1）－a5x31－a6x21y1］+

m［μD1y1，τc+τc（D2y1，τc+D1y2，τc）］（21）

方程 （19） 的解具有如下形式

Z1=GheiωτcT0+G—h—e－iωτcT0（22）

其中，G=G（T1，T2，T3，…），h见式（14）。将式（22）代入等式（20）的右侧，将会有eiωτcT0项的系数向量记为m1，通过可解条件〈h*，m1〉=0，可以求解

GT1=MμG（23）

其中

M=b1（a1+a2h2+mh2e－iωτc）－b1（1+τcmh2e－iωτc）－（a1－iω）（1－τcmh2e－iωτc）

（a1－iω）（b1+b2h2－mh2e－iωτc）b1（1+τcmh2e－iωτc）－（a1－iω）（1－τcmh2e－iωτc）

h2=iω－a1a2+me－iωτc（24）

接下来求解方程 （20） 。由于μ是一个较小的扰动参数，只考虑μ对低阶项的影响，则方程 （20） 有如下形式的解：

x2=g1G2e2iωτcT0+c.c.+l1GG—

y2=g2G2e2iωτcT0+c.c.+l2GG—（25）

其中c.c.表示前面项的共轭。将式（25）代入方程（20） ，可以求得

g1=（a3+a1h2）（2iω－b2+me－2iwτc）+（2iω－a1）（2iω－b2+me－2iwτc）－b1（a2+me－2iwτc）

（b3－a4h2）（a2+me－2iwτc）（2iω－a1）（2iω－b2+me－2iwτc）－b1（a2+me－2iwτc）

g2=b1（a3+a1h2）+（b3－a4h2）（2iω－a1）（2iω－a1）（2iω－b2+me－2iwτc）－b1（a2+me－2iwτc）

l1=（2a3+a4h2+a4h2）（－b2+m）a1（b2－m）－b1（a2+m）－

（2b3－a4h2－a4h2）（－a2－m）a1（b2－m）－b1（a2+m）

l2=b1（2a3+a4h2+a4h2）－a1（2bFLEtqB52YvrwGxkjOqHutb/GNYVzMiyKv0YcCXLS3EI=3－a4h2－a4h2）a1（b2－m）－b1（a2+m）

接下来将表达式 （22） 和 （25） 代入方程（21） 的右侧，并将含有eiωτcT0项的系数向量记为m2，通过可解条件〈h*，m2〉=0，则可以求解

GT2=HG2G—（26）

其中

H=b1R2－（a1－iω）R4－b1R1+（a1－iω）R3

R1=－1－τcmh2e－iωτc

R2=τc［2a3（P1P+S1S）+a4（P1Ph2+S1Sh2－

P2P－S2S）+3a5+a6（h2+2h2）］

R3=－1+τcmh2e－iωτc

R4=τc［2b3（P1P+S1S）－a4（P1Ph2+S1Sh2+

P2P+S2S）－3a5－a6（h2+2h2）］（27）

设GMT ExtraaA@G/ε，从而得到模型 （1） 截断到三次的Hopf分支的规范型

G·=MμG+HG2G—（28）

其中，M和H分别由等式 （24） 和 （27） 给出。

将G=reiθ代入规范型 （28） 中，得到如下极坐标下Hopf分支的规范型：

r·=Re（M）μr+Re（H）r3

θ·=Im（M）μr+Im（H）r3（29）

定理2 对于模型 （1） ，当Re（M）μRe（H）<0时，模型 （1） 平衡点E*附近存在Hopf分支周期解。

1）当Re（M）μ<0，分支周期解是不稳定的；若Re（M）<0（Re（M）>0），则分支周期解是前向的（后向的）。

2）当Re（M）μ>0，分支周期解是稳定的；若Re（M）<0 （Re（M）>0），则分支周期解是后向的（前向的）。

4 数值模拟

黑龙江省有我国最大的国有林区和森林资源基地［20］，在增强森林碳汇方面具有巨大潜力。本文以黑龙江省的森林资源数据为例，分析并预测各参数的取值范围，选择合理参数进行数值模拟，并分析模拟结果给出建议。

4.1 参数分析

1）健康树木面积的自然增长率r。

基于我国国家统计年鉴［21］、全国森林资源清查报告及黑龙江省森林所占面积的相关数据，计算得出1984年—2021年黑龙江省森林面积的年增长率，如表1所示。

根据表1，取数据的最大值与最小值作为参数范围的上线与下限，可以得出黑龙江省树木面积的自然增长率r∈（0.003 6，0.022 3）。

2）理想条件下该区域的森林最大容纳量N。

结合近40年的8组黑龙江省森林面积数据，对于这类“小样本”的不确定性系统，通过灰色预测GM（1，1）模型能够实现对系统运行行为、演化规律的描述和有效监控［22］。因此构建GM（1，1）模型来预测该区域的最大森林面积，借助SPSS软件进行模型计算。

为了使模型通过级比检验，取数据单位为万公顷，原始数据序列确定为

X（0）=（1 561.52，1 616.20，1 760.31，1 797.50，1 926.97，1 962.13，1 990.46，2 162.30）（30）

此时模型通过级比检验，说明数据适合构建GM（1，1）模型。接下来对数据进行拟合预测，结果如表2所示。

选用向后10期的数据，即大约50年后森林面积可以达到3 276.142万公顷。下面对模型情况进行检验，检验结果如表3所示。

模型相对误差值的最大值0.029<0.1说明模型拟合效果达到较高要求；模型级比偏差值的最大值0.042<0.1，说明模型拟合效果达到较高要求。

根据2022黑龙江省统计年鉴［23］中土地状况数据，黑龙江省土地总面积为4 706.9万公顷，预测值已经占总面积的70%。考虑到耕地、湿地等土地的面积对林地扩张的限制，认为森林的最大容纳量N∈（32，34）（单位：百万公顷）。

3）自然灾害对健康树木的有效接触率β。

肖风劲等［24］基于中国林业统计资料统计了中国森林生态系统的主要风险源概率，其中森林病虫害的总概率是82.71%。即健康树木有82.71%的概率感染森林病虫害，因此取β=0.827。

4）饱和参数α。

当α=0，βxy1+αx变为线性型发生率βxy，而α>0时的情况更贴近实际。当病虫害平衡点E*存在，说明饱和参数满足α<β。因此，饱和参数的范围选取为α∈（0，β）=（0，0.827），此时森林病虫害的接触感染率具备了饱和特性，且α的增大会导致接触感染率相对降低。

5）林区对森林病虫害的治理率m。

通过国家统计年鉴［21］收集到黑龙江省林区历年遭受森林病虫害的防治情况数据，如图1所示。

取历年数据的最大值与最小值作为参数范围的上限和下限，得出参数范围m∈（0.78，0.90）。

6）病虫害对树木的致死率b。

本文认为所有染病树木的发展情况最终分为两部分，一部分经治理后康复，另一部分最终死亡，则治理率和死亡率存在关系：b=1－m，如图2所示。因此b∈（0.10，0.22）。

4.2 数值模拟

根据4.1节的分析，选取多组数据的平均值，确定如下参数：

参数I：r=0.01; N=33; β=0.827; α=0.792; m=0.84; b=0.16;

参数I满足条件β>α（m+b），因此模型 （1） 存在病虫害平衡点E*=（28.57，0.24）。

1）首先分析模型 （1） 在τ=0时的情况，即染病树木的恢复时间是瞬时的。根据2021年黑龙江省森林面积情况，选取（21.62，0.46）为初值。参数I满足假设 （H1） ，平衡点E*局部渐近稳定，如图3所示。

由图3可知，健康树木的面积会先达到峰值，再缓慢下降至平衡状态后保持稳定；染病树木的面积会先达到最小值，再缓慢上升至平衡状态后保持稳定。事实上τ=0是较为理想的状态，树木染病过后存在恢复期，接下来将重点对τ>0的情况进行讨论。

2）考虑到不同的树种的恢复期不同，取τ=0.5模拟模型 （1） 病虫害平衡点E*附近的解的走向。基于黑龙江省森林资源清查报告，选取2021年的数据（x（t），y（t））=（21.62，0.46），t∈［－τ，0］作为初始函数。模拟结果如图4所示。

由图4所示，随着时间的推移，x和y在病虫害平衡点E*=（28.57，0.24）处是局部渐近稳定的。也就是说，当染病树木恢复期在临界值内，健康树木面积和染病树木面积会发生一段增减后达到平衡状态，森林病虫害的感染情况会最终趋于稳定。此外，考虑到健康树木、染病树木的面积在短期内的走向不同，建议调整病虫害的防治措施，可以将病虫害平衡点控制在更为理想的状态。

接下来讨论τ∈（0，τc）时，分别取τ1=0.3， τ2=0.6， τ3=0.9时模型 （1） 达到平衡状态时的现象。模拟结果如图5所示。

如图5所示，随着τ的增大，x和y在达平衡前的波动程度也越来越大。这XBY7OzdncQcsYmRh8wPm4q7R3yXXpAoZsQbdf7uM/gs=说明染病树木的恢复期越长，森林生态系统内部不稳定因素带来的影响程度越大。因此要尽量缩小染病树木的恢复期，以维持林区内部的稳定。

3）将参数I代入式（7），求得临界时滞τc=1.15。取（x（t），y（t））=（28，0.3），t∈［－τ，0］为初值函数。令τ=1.155>τc=1.15，进行数值模拟，如图6所示。

根据定理2，当恢复时间τ超过临界时滞τc且在临界时滞附近，模型 （1） 的病虫害平衡点E*的稳定性发生变化。当染病树木的恢复时间大于临界τc时，由规范型式（28）可得Re（M）=0.006 014>0，Re（H）=－0.000 120<0，则Re（M）μRe（H）<0且Re（M）μ>0。由定理2可知，模型（1）在τc附近有前向的稳定的周期解，如图6所示。该现象反应了染病树木的面积会长期在病虫害平衡点E*附近波动，保证治理率m∈cb2（a－r），+∞时，健康树木的面积和染病树木的面积将会是此消彼长的状态，短期内出现的感染面积增大、减小都为可控现象。此外，当恢复时间τ远大于临界时滞τc时，模型（1）病虫害平衡点E*不稳定。因此要将树木恢复时间尽可能控制在（0，τc）内，以维持森林生态系统的稳定。

4）参数α、r和m的变化对模型 （1） 病虫害平衡点E*的影响。

首先研究饱和系数α的变化。在3.1节中已经讨论了病虫害平衡点E*存在时的α的范围，选取3组不同的饱和参数α=0.785，α=0.790，α=0.795，代入参数I进行数值模拟，如图7所示。

由图7可知，饱和系数α越大，模型 （1） 稳定后健康树木所占的面积越大，染病树木所占面积越小。饱和发生率中的α反映了抑制作用［25］，因此考虑用α的大小衡量健康树木对病虫害的抵抗力。也就是说，α越大，树木的抵抗越强，健康树木越多，染病树木越少，从而更有益于我国实现碳中和。

接下来研究森林面积自然增长率r的变化对模型 （1） 的影响，模拟结果如图8所示。

森林面积的自然增长率r的增大，会导致模型 （1） 到达平衡状态所需时长减小，但同时染病树木在平衡状态时的面积也相应增大。这说明森林面积的扩张速度需要适度，速度过慢会导致达平衡的时间延长，速度过快会导致受病虫害影响的树木面积在所有树木中的占比增大，因此，将增长速率控制在适中的情况对森林生态系统的发展更有益。

下面研究林区森林病虫害的治理率m的变化对模型 （1） 的影响，模拟结果如图9所示。

森林病虫害的治理率越高，则模型 （1） 达到平衡状态后健康树木的所占面积越大，染病树木所占面积越小。且治理率越高，达平衡状态的速度越快，森林生态系统内部波动也越小。因此建议提高治理率，有利于扩大健康树木的面积、控制森林病虫害的蔓延，从而提升该林区的固碳能力。

5 结 论

本文基于森林病虫害的传播机理，建立了一个具时滞的含有饱和发生率的森林病虫害传播微分方程模型，讨论了平衡点存在的条件和稳定性，以及Hopf分支存在的条件。同时，采用多时间尺度法计算Hopf分支的规范型，分析了Hopf分支周期解的稳定性和分支方向。基于黑龙江省森林资源数据，确定了合理的参数范围，并通过选取适当参数、对照不同参数进行了多次数值模拟。最后根据模拟结果，在第4章中给出了结论和建议。

数值模拟部分的数据选取主要来源于黑龙江省林区的相关数据，事实上，该模型对于其他省份同样适用，可以用于预测森林病虫害的发展趋势。此外，林区还会遭受火灾、气候灾害等风险的影响而导致固碳能力下降，且这些风险的发生具有不确定性，这些因素将在后续的研究中进行讨论。

参 考 文 献：

［1］ 邹才能， 薛华庆， 熊波， 等. “碳中和”的内涵、创新与愿景［J］.天然气工业，2021，41（8）：46.

ZOU Caineng， XUE Huaqing， XIONG Bo， et al. Connotation， Innovation and Vision of “Carbon Neutral”［J］. Natural Gas Industry，2021，41（8）：46.

［2］ 马晓哲， 王铮. 土地利用变化对区域碳源汇的影响研究进展［J］.生态学报，2015，35（17）：5898.

MA Xiaozhe， WANG Zheng. Progress in the Study on the Impact of Land-use Change on Regional Carbon Sources and Sinks［J］. Acta Ecological Sinica，2015，35（17）：5898.

［3］ 杨忠岐， 王小艺， 张翌楠， 等. 以生物防治为主的综合控制我国重大林木病虫害研究进展［J］.中国生物防治学报，2018，34（2）：163.

YANG Zhongqi， WANG Xiaoyi， ZHANG Yinan， et al. Research Advances of Chinese Major Forest Pests by Integrated Management Based on Biological Control［J］. Chinese Journal of Biological Control，2018，34（2）：163.

［4］ 李丹， 王文帆， 刘滨凡. 黑龙江林业碳汇的探究［J］.科技创新与生产力，2022（2）：59.

LI Dan， WANG Wenfan， LIU Binfan. Research on Forestry Carbon Sequestration in Heilongjiang Province［J］. Sci-tech Innovation and Productivity，2022（2）：59.

［5］ 郭仲伟， 吴朝阳， 汪箫悦. 基于卫星遥感数据的森林病虫害监测与评价［J］.地理研究，2019，38（4）：831.

GUO Zhongwei， WU Chaoyang， WANG Xiaoyue. Forest Insect-disease Monitoring and Estimation Based on Satellite Remote Sensing Data［J］. Geographical Research，2019，38（4）：831.

［6］ YANG Bo， CAO Weiqun， TIAN Chengming. Visual analysis of impact factors of forest pests and diseases［J］. Journal of V-isualization，2019，22（6）：1257.

［7］ LINNAKOSKI R， KASANEN R， DOUNAVI A， et al. Editorial： Forest Health Under Climate Change： Effects on Tree Resilience， and Pest and Pathogen Dynamics［J］. Frontiers in Plant Science，2019，10：1157.

［8］ 施汉钰， 刘瑰琦， 曹昊， 等. 间作东北铁线莲对病虫害生物防治的研究［J］.森林工程，2018，34（5）：15.

SHI Hanyu， LIU Guiqi， CAO Hao， et al. The Biological Control Effect of Intercropping Clematis manshurica on Diseases and Insect Pests［J］. Forest Engineering，2018，34（5）：15.

［9］ 孙丽萍， 谭少亨， 周宏威， 等. 基于YOLOv5的林业有害生物检测与识别［J］.森林工程，2022，38（5）：104.

SUN Liping， TAN Shaoheng， ZHOU Hongwei， et al. Forestry Pests Detection and Identification Based on YOLOv5［J］. Forest Engineering，2022，38（5）：104.

［10］王建平， 高淑京， 邹圣龙. 一类柑橘黄龙病动力学模型的研究［J］.赣南师范学院学报，2013，34（3）：60.

WANG Jianping， GAO Shujing， ZOU Shenglong. Study on Dynamic Model of Huanglongbing［J］. Journal of Gannan Normal University，2013，34（3）：60.

［11］张蒙， 师向云. 一类具有时滞的松材线虫非线性传播模型的稳定性分析和Hopf分支［J］.信阳师范学院学报（自然科学版），2017，30（1）：17.

ZHANG Meng， SHI Xiangyun. Stability Analysis and Hopf Bifurcation of the Nonlinear Model for Pine Wilt Disease with Time Delay［J］. Journal of Xinyang Normal University （Natural Science Edition），2017，30（1）：17.

［12］杨婉君， 罗勇， 胡亦郑. 一类森林病虫害传染病模型研究［J］.温州大学学报（自然科学版），2018，39（2）：31.

YANG Wanjun， LUO Yong， HU Yizheng. The Study on an Epidemic Model of Forest Insect Pests［J］. Journal of Wenzhou University （Natural Science Edition），2018，39（2）：31.

［13］WANG Dingjiang， ZHANG Yingying. Stability Analysis of the Forest Insect Pests Model With Time Delays［J］. Journal of Biomathematics，2013，28（2）：211.

［14］ANDERSON R M， MAY R M. Population Biology of Infectious Diseases： Part I［J］. Nature，1979，280（5721）：361.

［15］许立滨， 李冬梅， 董在飞. 带有饱和发病率的离散SIR传染病模型的稳定性及分支问题［J］.哈尔滨理工大学学报，2017，22（3）：117.

XU Libin， LI Dongmei， DONG Zaifei. The Stability and Bifurcation Problems of Discrete SIR Epidemic Model with Saturation Incidence［J］. Journal of Harbin University of Science and Technology，2017，22（3）：117.

［16］BRAUER F， CASTILLO-CHAVEZ C. Mathematical Models in Population Biology and Epidemiology［M］. New York： Spring： 275.

［17］ANDERSON-TEIXEIRA K J， MILLER A D， MOHAN J E， et al. Altered Dynamics of Forest Recovery Under a Changing Climate［J］. Global Change Biology，2013，19（7）：2001.

［18］王晶囡， 逯兰芬， 杨德中， 等. 具时滞肿瘤T细胞免疫系统Hopf分支及肿瘤抑制［J］.哈尔滨理工大学学报，2020，25（2）：130.

WANG Jingnan， LU Lanfen， YANG Dezhong， ZHANG Yanqiao. Hopf Bifurcation of Delayed Tumor-T Cell Immune System and Inhibition of Tumor Growth［J］. Journal of Harbin University of Science and Technology，2020，25（2）：130.

［19］丁宇婷， 蒋丽， 苏琳琳. 时滞变频调压供水系统的稳定性及Hopf分支性质［J］.哈尔滨理工大学学报，2018，23（2）：149.

DING Yuting， JIANG Li， SU Linlin. Stability and Hopf Bifurcation of Delayed Variable Frequency Water-Supply System［J］. Journal of Harbin University of Science and Technology，2018，23（2）：149.

［20］李静， 焦文敬. 黑龙江国有林区社会救助多元主体协作问题研究——基于网络治理理论［J］.林业经济，2019，41（5）：26.

LI Jing， JIAO Wenjing. Research on Multi-subject Cooperation in Social Assistance of Heilongjiang State-owned Forest Regions——Based on Network Governance Theory［J］. Forestry Economics，2019，41（5）：26.

［21］中国统计年鉴［EB/OL］.中国统计出版社，2022. （2023-02-15）［2023-06-26］.

China Statistics Yearbook［EB/OL］. China Statistics Press，2022. （2023-02-15）［2023-06-26］.

［22］MA Xin， WU Wenqing， ZHANG Yuanyuan. Improved GM（1，1） Model Based on Simpson Formula and Its Applications［J］. The Journal of Grey System，2019，31（4）：33.

［23］黑龙江统计年鉴［EB/OL］.中国统计出版社，2022. （2023-02-13）［2023-06-27］.

Heilongjiang Statistics Yearbook［EB/OL］.China Statistics Press，2022. （2023-02-13）［2023-06-27］.

［24］肖风劲， 欧阳华， 程淑兰， 等. 中国森林健康生态风险评价［J］.应用生态学报，2004（2）：349.

XIAO Fengjin，OUYANG Hua，CHENG Shulan，et al. Forest Health Ecological Risk Assessment in China［J］. Chinese Journal of Applied Ecology，2004（2）：349.

［25］王婷， 王辉， 胡志兴. 一类非线性SEIRS传染病传播数学模型［J］.河南科技大学学报（自然科学版），2017，38（2）：84.

WANG Ting， WANG Hui， HU Zhixing. A Kind of Nonlinear SEIRS Epidemic Spread Mathematic Model［J］. Journal of Henan University of Science and Technology （Natural Science），2017，38（2）：84.

（编辑：温泽宇）