摘 要：控制森林病蟲(chóng)害有利于增加森林碳匯量，對(duì)實(shí)現(xiàn)我國(guó)“雙碳”戰(zhàn)略具有重要意義。為了研究森林病蟲(chóng)害對(duì)林區(qū)的影響，建立了一個(gè)含飽和發(fā)生率的森林病蟲(chóng)害傳播時(shí)滯微分方程模型，研究了平衡點(diǎn)的存在性與穩(wěn)定性，以及平衡點(diǎn)處Hopf分支的存在性，同時(shí)利用多尺度方法推導(dǎo)出系統(tǒng)在Hopf分支臨界點(diǎn)附近的規(guī)范型。結(jié)合全國(guó)森林清查報(bào)告等實(shí)際數(shù)據(jù)，對(duì)模型中的參數(shù)進(jìn)行分析與預(yù)測(cè)，得出各參數(shù)的取值范圍。選取合理參數(shù)進(jìn)行數(shù)值模擬，并研究了部分參數(shù)變化對(duì)模型的影響。最后根據(jù)模擬結(jié)果，結(jié)合黑龍江省林區(qū)森林病蟲(chóng)害的防控現(xiàn)狀給出結(jié)論和建議。

關(guān)鍵詞：森林病蟲(chóng)害；時(shí)滯微分方程；穩(wěn)定性；Hopf分支；多尺度方法；規(guī)范型

DOI：10.15938/j.jhust.2024.03.018

中圖分類號(hào)： O175.1

文獻(xiàn)標(biāo)志碼： A

文章編號(hào)： 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 引 言

在氣候異常變化的背景下，實(shí)現(xiàn)全球“碳中和”已經(jīng)成為國(guó)際共識(shí)，我國(guó)也提出了“在2060年前實(shí)現(xiàn)碳中和”這一戰(zhàn)略目標(biāo)［1］。森林生態(tài)系統(tǒng)包含了生物圈中大部分的碳，陸地植被與大氣之間有90%以上的碳交換由森林完成，且森林面積及覆蓋類型的變化將會(huì)引起全球碳循環(huán)的變動(dòng)［2］。然而我國(guó)是世界上受森林病蟲(chóng)害影響最為嚴(yán)重的國(guó)家之一［3］，以食葉害蟲(chóng)為例，葉片受到蟲(chóng)害后，葉片面積減少導(dǎo)致光合作用減弱，從而影響森林固碳量。有統(tǒng)計(jì)表明：每年植食性昆蟲(chóng)對(duì)森林生態(tài)系統(tǒng)的凈初級(jí)生產(chǎn)力破壞達(dá)到5%［4］。因此，我國(guó)森林病蟲(chóng)害問(wèn)題對(duì)實(shí)現(xiàn)“雙碳”戰(zhàn)略目標(biāo)的影響不容忽視。

目前，學(xué)者們對(duì)森林病蟲(chóng)害的探討主要在監(jiān)測(cè)識(shí)別、防治對(duì)策、損失評(píng)估等領(lǐng)域［5-9］。在林區(qū)中，雖然樹(shù)木彼此獨(dú)立，但是林區(qū)中的雜草、枝梗、媒介昆蟲(chóng)、土壤等均可以作為傳播媒介，導(dǎo)致病蟲(chóng)害在林區(qū)中蔓延。例如，柑橘的黃龍病通過(guò)媒介昆蟲(chóng)橘木虱傳播［10］，松樹(shù)的松材線蟲(chóng)病通過(guò)媒介昆蟲(chóng)松褐天牛傳播［11］……從而有很多學(xué)者基于傳染病模型建立了微分方程模型刻畫(huà)林區(qū)內(nèi)部健康樹(shù)木與染病樹(shù)木面積的連續(xù)變化的情況［12-13］。由于林區(qū)面積較大，實(shí)際情況中染病樹(shù)木變化率呈非線性變化，因此以飽和接觸率代替線性增長(zhǎng)率更具合理性［14-15］。此外，考慮到對(duì)染病樹(shù)木采取防治措施后存在一定的恢復(fù)時(shí)間，因此也將考慮樹(shù)木從染病狀態(tài)恢復(fù)至健康狀態(tài)的時(shí)間延遲。

綜合以上分析，本文建立一類含飽和發(fā)生率的時(shí)滯微分方程模型，分析該模型的平衡點(diǎn)的存在性、穩(wěn)定性，并討論Hopf分支的產(chǎn)生條件及其約化在中心流形上的規(guī)范型。

1 數(shù)學(xué)建模

在沒(méi)有外界干預(yù)的情況下，林區(qū)面積的增長(zhǎng)情況符合logistic模型的形式。森林病蟲(chóng)害發(fā)生后，會(huì)通過(guò)媒介昆蟲(chóng)傳播至健康樹(shù)木，導(dǎo)致健康樹(shù)木面積的擴(kuò)張速率減小，而染病樹(shù)木面積加速增長(zhǎng)。考慮到一棵染病樹(shù)木的影響能力有限，此處用飽和發(fā)生率來(lái)刻畫(huà)病蟲(chóng)害的感染情況。森林病蟲(chóng)害經(jīng)過(guò)治理后，有一部分樹(shù)木會(huì)被治愈，小部分樹(shù)木會(huì)由于受害過(guò)于嚴(yán)重而失去固碳能力。

基于現(xiàn)有的流行病學(xué)數(shù)學(xué)模型［16］，構(gòu)建如下含飽和發(fā)生率的時(shí)滯微分方程模型。x（t），y（t）分別表示t時(shí)刻該林區(qū)健康樹(shù)木和染病樹(shù)木的面積?？紤]到樹(shù)木的自然增長(zhǎng)與死亡，建立如下模型：

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

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

其中：r為健康林區(qū)面積的自然增長(zhǎng)率；N為健康林區(qū)的環(huán)境容納量；βxy/（1+αx）為森林病蟲(chóng)害傳播的飽和接觸率；β為有效接觸率；α為飽和參數(shù)；m為林區(qū)森林病蟲(chóng)害的治理率；b為森林病蟲(chóng)害對(duì)該林區(qū)樹(shù)木的致死率；τ為染病樹(shù)木的恢復(fù)時(shí)間。假設(shè)模型中所有系數(shù)均為正數(shù)。

在接下來(lái)的章節(jié)中主要分析模型 （1） 的動(dòng)力學(xué)性質(zhì)。

2 模型分析

模型 （1） 始終存在無(wú)病蟲(chóng)害平衡點(diǎn)E0=（N，0）。當(dāng)β>α（m+b）時(shí)，系統(tǒng) （1） 存在病蟲(chóng)害平衡點(diǎn)E*=（x*，y*），其中

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

下面主要討論病蟲(chóng)害平衡點(diǎn)E*附近的解的動(dòng)力學(xué)性質(zhì)。

模型 （1） 在病蟲(chóng)害平衡點(diǎn)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*

當(dāng)τ=0時(shí)，特征方程（2）變?yōu)槿缦滦问剑?/p>

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

給出如下假設(shè)

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

當(dāng)假設(shè)（H1）成立時(shí)，特征方程（3）有兩個(gè)負(fù)實(shí)根，模型（1）的病蟲(chóng)害平衡點(diǎn)E*在τ=0時(shí)是局部漸近穩(wěn)定的。

當(dāng)τ>0時(shí)，由于m≠n，則有c（m－n）≠0，因此λ=0不是方程的根，模型（1）在E*處無(wú)不動(dòng)點(diǎn)分支。假設(shè)λ=iω為方程（2）的根，代入特征方程 （2）并分離實(shí)虛部有

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

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

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

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

根據(jù)方程 （4） 可以推導(dǎo)出如下結(jié)果：

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）

設(shè)z=ω2，方程 （4） 化為

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

其中

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

當(dāng)D2<0即m∈cb2（a－r），+∞時(shí)，方程 （6） 有一個(gè)正根，記為z1，則有h′（z1）>0；當(dāng)D1<0，D2>0，即m∈0，cb2（a－r），n∈（a+c－r）22c，+∞時(shí)，方程 （6） 有兩個(gè)正根，記為z2<z3，則有h′（z2）<0，h′（z3）>0。

不失一般性，假設(shè)方程 （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 當(dāng)假設(shè) （H1） 成立，模型 （1） 平衡點(diǎn)E*的穩(wěn)定性和Hopf分支的存在性有如下結(jié)論：

1）當(dāng)m∈0，cb2（a－r），n∈0，（a+c－r）22c時(shí)，對(duì)于τ≥0，平衡點(diǎn)E*是局部漸近穩(wěn)定的，模型 （1） 無(wú)Hopf分支。

2）當(dāng)m∈cb2（a－r），+∞時(shí)，方程 （6） 有一個(gè)正根z1。當(dāng)τ∈［0，τ（0）1），平衡點(diǎn)E*是局部漸近穩(wěn)定的；當(dāng)τ≥τ（0）1，平衡點(diǎn)E*是不穩(wěn)定的；當(dāng)τ=τ（j）1 （j=0，1，2，…），平衡點(diǎn)E*處產(chǎn)生Hopf分支。

3）當(dāng)m∈0，cb2（a－r），n∈（a+c－r）22c，+∞時(shí)，方程 （6） 有兩個(gè)正根z2和z3（假設(shè)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。當(dāng)τ∈［0，τ（0）3）∪∪ml=1（τ（l－1）2，τ（l）3），平衡點(diǎn)E*是局部漸近穩(wěn)定的；當(dāng)τ∈∪m－1l=0（τ（l）3，τ（l）2）∪（τ（m）3，+∞），平衡點(diǎn)E*是不穩(wěn)定的；當(dāng)τ=τ（j）k （k=2，3. j=0，1，2，…），平衡點(diǎn)E*處產(chǎn)生Hopf分支。

根據(jù)參數(shù)含義，參數(shù)m和n分別可以表示森林病蟲(chóng)害的治理率和傳播率。該定理中m和n的取值范圍與病蟲(chóng)害平衡點(diǎn)處分支現(xiàn)象的生物學(xué)意義相統(tǒng)一，將在4.2節(jié)中深入討論。

3 Hopf分支的規(guī)范型

樹(shù)木受到病蟲(chóng)害干擾或因病蟲(chóng)害致死后，森林的固碳能力下降，而受干擾后森林的更新過(guò)程會(huì)對(duì)森林結(jié)構(gòu)與功能等多個(gè)方面產(chǎn)生影響［17］。因此關(guān)注染病樹(shù)木的恢復(fù)時(shí)間是很必要的，選取染病樹(shù)木的恢復(fù)時(shí)間τ作為分支參數(shù)。假設(shè)模型 （1） 的病蟲(chóng)害平衡點(diǎn)E*在τ=τc=τ（j）k處經(jīng)歷Hopf分支，此時(shí)特征方程 （2） 有一對(duì)純虛根λ=±iω，其中τ（j）k由式（7）給出。Hopf分支的研究方法有多種［18-19］，下面利用多尺度方法推導(dǎo)模型 （1） 在病蟲(chóng)害平衡點(diǎn)E*附近的Hopf分支規(guī)范型。

首先，令tMT ExtraaA@t/τ，將病蟲(chóng)害平衡點(diǎn)E*轉(zhuǎn)移到原點(diǎn)，即

=x－x*

=y－y*（9）

為了方便表達(dá)，仍然用x和y表示和，則方程 （2） 可寫(xiě)為

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） 中的分式在病蟲(chóng)害平衡點(diǎn)處進(jìn)行泰勒展開(kāi)，得到如下方程：

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，

為了方便計(jì)算，將方程 （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

設(shè)h為方程 （12） 的線性矩陣對(duì)應(yīng)于特征值λ=iωτ的特征向量，h*為方程 （12） 的線性矩陣的伴隨矩陣對(duì)應(yīng)于特征值λ=－iωτ的特征向量，滿足

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

通過(guò)計(jì)算可得，

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

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

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

假設(shè)方程 （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，…）處進(jìn)行泰勒展開(kāi)，得到

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，…。

將τ視為分支參數(shù)，設(shè)τ=τc+εμ，其中τc=τ（j）k為Hopf分支臨界值，μ為擾動(dòng)參數(shù)，ε為無(wú)量綱尺度參數(shù)。將等式（15）～（18）代入方程（12），分別比較等號(hào)左右兩邊ε，ε2，ε3的系數(shù)，得到如下表達(dá)式：

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見(jiàn)式（14）。將式（22）代入等式（20）的右側(cè)，將會(huì)有eiωτcT0項(xiàng)的系數(shù)向量記為m1，通過(guò)可解條件〈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）

接下來(lái)求解方程 （20） 。由于μ是一個(gè)較小的擾動(dòng)參數(shù)，只考慮μ對(duì)低階項(xiàng)的影響，則方程 （20） 有如下形式的解：

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

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

其中c.c.表示前面項(xiàng)的共軛。將式（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）

接下來(lái)將表達(dá)式 （22） 和 （25） 代入方程（21） 的右側(cè)，并將含有eiωτcT0項(xiàng)的系數(shù)向量記為m2，通過(guò)可解條件〈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）

設(shè)GMT ExtraaA@G/ε，從而得到模型 （1） 截?cái)嗟饺蔚腍opf分支的規(guī)范型

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

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

將G=reiθ代入規(guī)范型 （28） 中，得到如下極坐標(biāo)下Hopf分支的規(guī)范型：

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

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

定理2 對(duì)于模型 （1） ，當(dāng)Re（M）μRe（H）<0時(shí)，模型 （1） 平衡點(diǎn)E*附近存在Hopf分支周期解。

1）當(dāng)Re（M）μ<0，分支周期解是不穩(wěn)定的；若Re（M）<0（Re（M）>0），則分支周期解是前向的（后向的）。

2）當(dāng)Re（M）μ>0，分支周期解是穩(wěn)定的；若Re（M）<0 （Re（M）>0），則分支周期解是后向的（前向的）。

4 數(shù)值模擬

黑龍江省有我國(guó)最大的國(guó)有林區(qū)和森林資源基地［20］，在增強(qiáng)森林碳匯方面具有巨大潛力。本文以黑龍江省的森林資源數(shù)據(jù)為例，分析并預(yù)測(cè)各參數(shù)的取值范圍，選擇合理參數(shù)進(jìn)行數(shù)值模擬，并分析模擬結(jié)果給出建議。

4.1 參數(shù)分析

1）健康樹(shù)木面積的自然增長(zhǎng)率r。

基于我國(guó)國(guó)家統(tǒng)計(jì)年鑒［21］、全國(guó)森林資源清查報(bào)告及黑龍江省森林所占面積的相關(guān)數(shù)據(jù)，計(jì)算得出1984年—2021年黑龍江省森林面積的年增長(zhǎng)率，如表1所示。

根據(jù)表1，取數(shù)據(jù)的最大值與最小值作為參數(shù)范圍的上線與下限，可以得出黑龍江省樹(shù)木面積的自然增長(zhǎng)率r∈（0.003 6，0.022 3）。

2）理想條件下該區(qū)域的森林最大容納量N。

結(jié)合近40年的8組黑龍江省森林面積數(shù)據(jù)，對(duì)于這類“小樣本”的不確定性系統(tǒng)，通過(guò)灰色預(yù)測(cè)GM（1，1）模型能夠?qū)崿F(xiàn)對(duì)系統(tǒng)運(yùn)行行為、演化規(guī)律的描述和有效監(jiān)控［22］。因此構(gòu)建GM（1，1）模型來(lái)預(yù)測(cè)該區(qū)域的最大森林面積，借助SPSS軟件進(jìn)行模型計(jì)算。

為了使模型通過(guò)級(jí)比檢驗(yàn)，取數(shù)據(jù)單位為萬(wàn)公頃，原始數(shù)據(jù)序列確定為

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）

此時(shí)模型通過(guò)級(jí)比檢驗(yàn)，說(shuō)明數(shù)據(jù)適合構(gòu)建GM（1，1）模型。接下來(lái)對(duì)數(shù)據(jù)進(jìn)行擬合預(yù)測(cè)，結(jié)果如表2所示。

選用向后10期的數(shù)據(jù)，即大約50年后森林面積可以達(dá)到3 276.142萬(wàn)公頃。下面對(duì)模型情況進(jìn)行檢驗(yàn)，檢驗(yàn)結(jié)果如表3所示。

模型相對(duì)誤差值的最大值0.029<0.1說(shuō)明模型擬合效果達(dá)到較高要求；模型級(jí)比偏差值的最大值0.042<0.1，說(shuō)明模型擬合效果達(dá)到較高要求。

根據(jù)2022黑龍江省統(tǒng)計(jì)年鑒［23］中土地狀況數(shù)據(jù)，黑龍江省土地總面積為4 706.9萬(wàn)公頃，預(yù)測(cè)值已經(jīng)占總面積的70%?？紤]到耕地、濕地等土地的面積對(duì)林地?cái)U(kuò)張的限制，認(rèn)為森林的最大容納量N∈（32，34）（單位：百萬(wàn)公頃）。

3）自然災(zāi)害對(duì)健康樹(shù)木的有效接觸率β。

肖風(fēng)勁等［24］基于中國(guó)林業(yè)統(tǒng)計(jì)資料統(tǒng)計(jì)了中國(guó)森林生態(tài)系統(tǒng)的主要風(fēng)險(xiǎn)源概率，其中森林病蟲(chóng)害的總概率是82.71%。即健康樹(shù)木有82.71%的概率感染森林病蟲(chóng)害，因此取β=0.827。

4）飽和參數(shù)α。

當(dāng)α=0，βxy1+αx變?yōu)榫€性型發(fā)生率βxy，而α>0時(shí)的情況更貼近實(shí)際。當(dāng)病蟲(chóng)害平衡點(diǎn)E*存在，說(shuō)明飽和參數(shù)滿足α<β。因此，飽和參數(shù)的范圍選取為α∈（0，β）=（0，0.827），此時(shí)森林病蟲(chóng)害的接觸感染率具備了飽和特性，且α的增大會(huì)導(dǎo)致接觸感染率相對(duì)降低。

5）林區(qū)對(duì)森林病蟲(chóng)害的治理率m。

通過(guò)國(guó)家統(tǒng)計(jì)年鑒［21］收集到黑龍江省林區(qū)歷年遭受森林病蟲(chóng)害的防治情況數(shù)據(jù)，如圖1所示。

取歷年數(shù)據(jù)的最大值與最小值作為參數(shù)范圍的上限和下限，得出參數(shù)范圍m∈（0.78，0.90）。

6）病蟲(chóng)害對(duì)樹(shù)木的致死率b。

本文認(rèn)為所有染病樹(shù)木的發(fā)展情況最終分為兩部分，一部分經(jīng)治理后康復(fù)，另一部分最終死亡，則治理率和死亡率存在關(guān)系：b=1－m，如圖2所示。因此b∈（0.10，0.22）。

4.2 數(shù)值模擬

根據(jù)4.1節(jié)的分析，選取多組數(shù)據(jù)的平均值，確定如下參數(shù)：

參數(shù)I：r=0.01; N=33; β=0.827; α=0.792; m=0.84; b=0.16;

參數(shù)I滿足條件β>α（m+b），因此模型 （1） 存在病蟲(chóng)害平衡點(diǎn)E*=（28.57，0.24）。

1）首先分析模型 （1） 在τ=0時(shí)的情況，即染病樹(shù)木的恢復(fù)時(shí)間是瞬時(shí)的。根據(jù)2021年黑龍江省森林面積情況，選取（21.62，0.46）為初值。參數(shù)I滿足假設(shè) （H1） ，平衡點(diǎn)E*局部漸近穩(wěn)定，如圖3所示。

由圖3可知，健康樹(shù)木的面積會(huì)先達(dá)到峰值，再緩慢下降至平衡狀態(tài)后保持穩(wěn)定；染病樹(shù)木的面積會(huì)先達(dá)到最小值，再緩慢上升至平衡狀態(tài)后保持穩(wěn)定。事實(shí)上τ=0是較為理想的狀態(tài)，樹(shù)木染病過(guò)后存在恢復(fù)期，接下來(lái)將重點(diǎn)對(duì)τ>0的情況進(jìn)行討論。

2）考慮到不同的樹(shù)種的恢復(fù)期不同，取τ=0.5模擬模型 （1） 病蟲(chóng)害平衡點(diǎn)E*附近的解的走向?；诤邶埥∩仲Y源清查報(bào)告，選取2021年的數(shù)據(jù)（x（t），y（t））=（21.62，0.46），t∈［－τ，0］作為初始函數(shù)。模擬結(jié)果如圖4所示。

由圖4所示，隨著時(shí)間的推移，x和y在病蟲(chóng)害平衡點(diǎn)E*=（28.57，0.24）處是局部漸近穩(wěn)定的。也就是說(shuō)，當(dāng)染病樹(shù)木恢復(fù)期在臨界值內(nèi)，健康樹(shù)木面積和染病樹(shù)木面積會(huì)發(fā)生一段增減后達(dá)到平衡狀態(tài)，森林病蟲(chóng)害的感染情況會(huì)最終趨于穩(wěn)定。此外，考慮到健康樹(shù)木、染病樹(shù)木的面積在短期內(nèi)的走向不同，建議調(diào)整病蟲(chóng)害的防治措施，可以將病蟲(chóng)害平衡點(diǎn)控制在更為理想的狀態(tài)。

接下來(lái)討論τ∈（0，τc）時(shí)，分別取τ1=0.3， τ2=0.6， τ3=0.9時(shí)模型 （1） 達(dá)到平衡狀態(tài)時(shí)的現(xiàn)象。模擬結(jié)果如圖5所示。

如圖5所示，隨著τ的增大，x和y在達(dá)平衡前的波動(dòng)程度也越來(lái)越大。這XBY7OzdncQcsYmRh8wPm4q7R3yXXpAoZsQbdf7uM/gs=說(shuō)明染病樹(shù)木的恢復(fù)期越長(zhǎng)，森林生態(tài)系統(tǒng)內(nèi)部不穩(wěn)定因素帶來(lái)的影響程度越大。因此要盡量縮小染病樹(shù)木的恢復(fù)期，以維持林區(qū)內(nèi)部的穩(wěn)定。

3）將參數(shù)I代入式（7），求得臨界時(shí)滯τc=1.15。?。▁（t），y（t））=（28，0.3），t∈［－τ，0］為初值函數(shù)。令τ=1.155>τc=1.15，進(jìn)行數(shù)值模擬，如圖6所示。

根據(jù)定理2，當(dāng)恢復(fù)時(shí)間τ超過(guò)臨界時(shí)滯τc且在臨界時(shí)滯附近，模型 （1） 的病蟲(chóng)害平衡點(diǎn)E*的穩(wěn)定性發(fā)生變化。當(dāng)染病樹(shù)木的恢復(fù)時(shí)間大于臨界τc時(shí)，由規(guī)范型式（28）可得Re（M）=0.006 014>0，Re（H）=－0.000 120<0，則Re（M）μRe（H）<0且Re（M）μ>0。由定理2可知，模型（1）在τc附近有前向的穩(wěn)定的周期解，如圖6所示。該現(xiàn)象反應(yīng)了染病樹(shù)木的面積會(huì)長(zhǎng)期在病蟲(chóng)害平衡點(diǎn)E*附近波動(dòng)，保證治理率m∈cb2（a－r），+∞時(shí)，健康樹(shù)木的面積和染病樹(shù)木的面積將會(huì)是此消彼長(zhǎng)的狀態(tài)，短期內(nèi)出現(xiàn)的感染面積增大、減小都為可控現(xiàn)象。此外，當(dāng)恢復(fù)時(shí)間τ遠(yuǎn)大于臨界時(shí)滯τc時(shí)，模型（1）病蟲(chóng)害平衡點(diǎn)E*不穩(wěn)定。因此要將樹(shù)木恢復(fù)時(shí)間盡可能控制在（0，τc）內(nèi)，以維持森林生態(tài)系統(tǒng)的穩(wěn)定。

4）參數(shù)α、r和m的變化對(duì)模型 （1） 病蟲(chóng)害平衡點(diǎn)E*的影響。

首先研究飽和系數(shù)α的變化。在3.1節(jié)中已經(jīng)討論了病蟲(chóng)害平衡點(diǎn)E*存在時(shí)的α的范圍，選取3組不同的飽和參數(shù)α=0.785，α=0.790，α=0.795，代入?yún)?shù)I進(jìn)行數(shù)值模擬，如圖7所示。

由圖7可知，飽和系數(shù)α越大，模型 （1） 穩(wěn)定后健康樹(shù)木所占的面積越大，染病樹(shù)木所占面積越小。飽和發(fā)生率中的α反映了抑制作用［25］，因此考慮用α的大小衡量健康樹(shù)木對(duì)病蟲(chóng)害的抵抗力。也就是說(shuō)，α越大，樹(shù)木的抵抗越強(qiáng)，健康樹(shù)木越多，染病樹(shù)木越少，從而更有益于我國(guó)實(shí)現(xiàn)碳中和。

接下來(lái)研究森林面積自然增長(zhǎng)率r的變化對(duì)模型 （1） 的影響，模擬結(jié)果如圖8所示。

森林面積的自然增長(zhǎng)率r的增大，會(huì)導(dǎo)致模型 （1） 到達(dá)平衡狀態(tài)所需時(shí)長(zhǎng)減小，但同時(shí)染病樹(shù)木在平衡狀態(tài)時(shí)的面積也相應(yīng)增大。這說(shuō)明森林面積的擴(kuò)張速度需要適度，速度過(guò)慢會(huì)導(dǎo)致達(dá)平衡的時(shí)間延長(zhǎng)，速度過(guò)快會(huì)導(dǎo)致受病蟲(chóng)害影響的樹(shù)木面積在所有樹(shù)木中的占比增大，因此，將增長(zhǎng)速率控制在適中的情況對(duì)森林生態(tài)系統(tǒng)的發(fā)展更有益。

下面研究林區(qū)森林病蟲(chóng)害的治理率m的變化對(duì)模型 （1） 的影響，模擬結(jié)果如圖9所示。

森林病蟲(chóng)害的治理率越高，則模型 （1） 達(dá)到平衡狀態(tài)后健康樹(shù)木的所占面積越大，染病樹(shù)木所占面積越小。且治理率越高，達(dá)平衡狀態(tài)的速度越快，森林生態(tài)系統(tǒng)內(nèi)部波動(dòng)也越小。因此建議提高治理率，有利于擴(kuò)大健康樹(shù)木的面積、控制森林病蟲(chóng)害的蔓延，從而提升該林區(qū)的固碳能力。

5 結(jié) 論

本文基于森林病蟲(chóng)害的傳播機(jī)理，建立了一個(gè)具時(shí)滯的含有飽和發(fā)生率的森林病蟲(chóng)害傳播微分方程模型，討論了平衡點(diǎn)存在的條件和穩(wěn)定性，以及Hopf分支存在的條件。同時(shí)，采用多時(shí)間尺度法計(jì)算Hopf分支的規(guī)范型，分析了Hopf分支周期解的穩(wěn)定性和分支方向。基于黑龍江省森林資源數(shù)據(jù)，確定了合理的參數(shù)范圍，并通過(guò)選取適當(dāng)參數(shù)、對(duì)照不同參數(shù)進(jìn)行了多次數(shù)值模擬。最后根據(jù)模擬結(jié)果，在第4章中給出了結(jié)論和建議。

數(shù)值模擬部分的數(shù)據(jù)選取主要來(lái)源于黑龍江省林區(qū)的相關(guān)數(shù)據(jù)，事實(shí)上，該模型對(duì)于其他省份同樣適用，可以用于預(yù)測(cè)森林病蟲(chóng)害的發(fā)展趨勢(shì)。此外，林區(qū)還會(huì)遭受火災(zāi)、氣候?yàn)?zāi)害等風(fēng)險(xiǎn)的影響而導(dǎo)致固碳能力下降，且這些風(fēng)險(xiǎn)的發(fā)生具有不確定性，這些因素將在后續(xù)的研究中進(jìn)行討論。

參 考 文 獻(xiàn)：

［1］ 鄒才能， 薛華慶， 熊波， 等. “碳中和”的內(nèi)涵、創(chuàng)新與愿景［J］.天然氣工業(yè)，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］ 馬曉哲， 王錚. 土地利用變化對(duì)區(qū)域碳源匯的影響研究進(jìn)展［J］.生態(tài)學(xué)報(bào)，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］ 楊忠岐， 王小藝， 張翌楠， 等. 以生物防治為主的綜合控制我國(guó)重大林木病蟲(chóng)害研究進(jìn)展［J］.中國(guó)生物防治學(xué)報(bào)，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］ 李丹， 王文帆， 劉濱凡. 黑龍江林業(yè)碳匯的探究［J］.科技創(chuàng)新與生產(chǎn)力，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］ 郭仲偉， 吳朝陽(yáng)， 汪簫悅. 基于衛(wèi)星遙感數(shù)據(jù)的森林病蟲(chóng)害監(jiān)測(cè)與評(píng)價(jià)［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］ 施漢鈺， 劉瑰琦， 曹昊， 等. 間作東北鐵線蓮對(duì)病蟲(chóng)害生物防治的研究［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的林業(yè)有害生物檢測(cè)與識(shí)別［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］王建平， 高淑京， 鄒圣龍. 一類柑橘黃龍病動(dòng)力學(xué)模型的研究［J］.贛南師范學(xué)院學(xué)報(bào)，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］張蒙， 師向云. 一類具有時(shí)滯的松材線蟲(chóng)非線性傳播模型的穩(wěn)定性分析和Hopf分支［J］.信陽(yáng)師范學(xué)院學(xué)報(bào)（自然科學(xué)版），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］楊婉君， 羅勇， 胡亦鄭. 一類森林病蟲(chóng)害傳染病模型研究［J］.溫州大學(xué)學(xué)報(bào)（自然科學(xué)版），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］許立濱， 李冬梅， 董在飛. 帶有飽和發(fā)病率的離散SIR傳染病模型的穩(wěn)定性及分支問(wèn)題［J］.哈爾濱理工大學(xué)學(xué)報(bào)，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］王晶囡， 逯蘭芬， 楊德中， 等. 具時(shí)滯腫瘤T細(xì)胞免疫系統(tǒng)Hopf分支及腫瘤抑制［J］.哈爾濱理工大學(xué)學(xué)報(bào)，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］丁宇婷， 蔣麗， 蘇琳琳. 時(shí)滯變頻調(diào)壓供水系統(tǒng)的穩(wěn)定性及Hopf分支性質(zhì)［J］.哈爾濱理工大學(xué)學(xué)報(bào)，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］李靜， 焦文敬. 黑龍江國(guó)有林區(qū)社會(huì)救助多元主體協(xié)作問(wèn)題研究——基于網(wǎng)絡(luò)治理理論［J］.林業(yè)經(jīng)濟(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］中國(guó)統(tǒng)計(jì)年鑒［EB/OL］.中國(guó)統(tǒng)計(jì)出版社，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］黑龍江統(tǒng)計(jì)年鑒［EB/OL］.中國(guó)統(tǒng)計(jì)出版社，2022. （2023-02-13）［2023-06-27］.

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

［24］肖風(fēng)勁， 歐陽(yáng)華， 程淑蘭， 等. 中國(guó)森林健康生態(tài)風(fēng)險(xiǎn)評(píng)價(jià)［J］.應(yīng)用生態(tài)學(xué)報(bào)，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傳染病傳播數(shù)學(xué)模型［J］.河南科技大學(xué)學(xué)報(bào)（自然科學(xué)版），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.

（編輯：溫澤宇）