一類具時滯森林病蟲害傳染病模型

摘 要：控制森林病蟲害有利于增加森林碳匯量，對實現我國“雙碳”戰略具有重要意義。為了研究森林病蟲害對林區的影響，建立了一個含飽和發生率的森林病蟲害傳播時滯微分方程模型，研究了平衡點的存在性與穩定性，以及平衡點處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.

（編輯：溫澤宇）