具有避難所和Ornstein-Uhlenbeck過程的隨機三種群捕食模型的持久與滅絕

摘要： 目前大多數的捕食模型都是用白噪聲描述環境的隨機性，針對真實情境下，模型中的參數可能滿足Ornstein-Uhlenbeck過程的問題，提出具有避難所和Ornstein-Uhlenbeck過程的隨機三種群捕食模型。利用伊藤公式、微分不等式及隨機分析理論對模型進行動力學行為分析。證明了模型全局正解的存在唯一性，分別得到了每個種群平均持久生存和滅絕的充分條件，最后利用python進行數值模擬，驗證了定理中所得到的結論，并進一步研究了避難所對種群數量的影響。

關鍵詞： Ornstein-Uhlenbeck過程；避難所；平均持久性；滅絕性

中圖分類號："" O175; TB3文獻標識碼： A

Persistence and Extinction of a Stochastic Three-population Predation

Model with Refuge and Ornstein-Uhlenbeck Process

SU Xiaoming， CHEN Bo

（School of Science， Shenyang University of Technology， Shenyang 110870， China）

Abstract：

Most current predation models describe the stochastic nature of the environment in terms of white noise， for the problem that the parameters in the model may satisfy the Ornstein-Uhlenbeck process in a real situation， a stochastic three-species predation model with refuge and Ornstein-Uhlenbeck process is proposed. The dynamic behavior of the model is analyzed by its formula， differential inequality and stochastic analysis theory. The existence and uniqueness of the global positive solution of the model are proved， sufficient conditions for the average persistence and extinction of each population were obtained separately， finally， python is used to carry out numerical simulation to verify the conclusions obtained in the theorem， and further study the impact of refuge on population size.

Keywords：

Ornstein-Uhlenbeck processes; refuge; average persistence; extinction

0 引言

捕食模型作為生態學和數學生態學中重要的種群模型之一，由于能準確捕捉到種群的動態行為變化而得到學者們的廣泛研究[13]。近年來，許多學者在捕食模型中考慮了獵物避難所[48]。在自然環境之中，獵物為了有更多的機會得以生存，會利用各種避難所來迷惑正在追捕的捕食者，比如：洞穴、喬木、茂密的植被等[910]。Yue[5]提出并研究了如下具有改進的Leslie-Gower、Holling-type II功能反應函數和獵物避難所的捕食模型：

dx（t）dt=xf1－b1x－a1－my1－mx+kdy（t）dt=yf2－ry1－mx+k（1）

其中，x（t），y（t）分別表示獵物種群和捕食者種群在t時刻的密度；參數f1，f2，b1，a，r，k，m都是正常數。f1，f2分別代表獵物x和捕食者y的增長率；b1表示獵物x的種內競爭系數；a度量了由于y的捕食導致的獵物x的最大單位減少率；k代表環境對獵物x和捕食者y的保護率；m∈［0，1］表示避難所系數，（1－m）x代表可供捕食的獵物數量；r和a具有相似的生物學意義。

文獻[5]考慮的是兩種群的模型，然而在現實環境中，存在著多個捕食者去捕食一個獵物的情況。因此，受文獻[11]的啟發，本文將模型（1）改進為兩個捕食者捕食一個獵物的模型：

dx（t）dt=xf1－b1x－a11－my11－mx+k1－a21－my21－mx+k2

dy1（t）dt=y1f2－r11－my11－mx+k1

dy2（t）dt=y2f3－r21－my21－mx+k2（2）

其中，x（t），y1（t），y2（t）分別代表獵物和兩個捕食者在t時刻的種群密度，其它參數的生物學意義與（1）相同。

另一方面，種群的動態行為不可避免地會受到環境波動的影響，學者們常用的方法是往模型中引入白噪聲，建立隨機捕食模型，而在實際情況中，模型中的參數可能滿足Ornstein-Uhlenbeck過程（均值回歸過程）[1216]。受文獻[1516]的啟發，本文假設模型（2）中的增長率fi （i=1，2，3）滿足均值回歸過程：

df～i（t）=βifi－f～i（t）dt+ηidBi（t），i=1，2，3（3）

其中，Bi（t） （i=1，2，3）為標準布朗運動。對式（3）應用文[17]中的伊藤公式：

f～i（t）=fi+fi0－fie－βit+ηi∫t0e－βit－sdBi（s）

=fi+fi0－fie－βit+σi（t）dBi（t）dt，i=1，2，3（4）

其中，fi0=f～i0，σi（t）=ηi2βi1－e－2βit，βigt;0為回歸速率，ηi為波動強度。

將式（4）代入模型（2）可得隨機捕食模型：

dx（t）=xf1+f10－f1e－β1t－b1x－a11－my11－mx+k1－a21－my21－mx+k2dt+xσ1（t）dB1（t）

dy1（t）=y1f2+f20－f2e－β2t－r11－my11－mx+k1dt+y1σ2（t）dB2（t）

dy2（t）=y2f3+f30－f3e－β3t－r21－my21－mx+k2dt+y2σ3（t）dB3（t）（5）

為討論方便，令

R3+=z∈R3|zigt;0，i=1，2，3，gi（t）=fi－η2i4βi+η2ie－2βit4βi

g-i=limt→+SymboleB@t－1∫t0gisds=fi－η2i4βii=1，2，3

x（0）=x0gt;0，y1（0）=y10gt;0，y2（0）=y20gt;0R+=0，+SymboleB@，X（t）=x（t），y1（t），y2（t）

1 模型正解的存在唯一性

定理1 對任意初值x0，y10，y20∈R3+，模型（5）存在唯一的全局正解x（t），y1（t）y2（t）∈R3+，且依概率1存在于R3+中。

證明：作變換

x（t）=eu（t），y1（t）=ev1（t），y2（t）=ev2（t）

則

u（t）=lnx（t），v1（t）=lny1（t），v2（t）=lny2（t）

對模型（5）應用伊藤公式得

du（t）=g1（t）+f10－f1e－β1t－b1eu（t）－a11－mev1（t）1－meu（t）+k1－a21－mev2（t）1－meu（t）+k2dt+σ1（t）dB1（t）

dv1（t）=g2（t）+f20－f2e－β2t－r11－mev1（t）1－meu（t）+k1dt+σ2（t）dB2（t）

dv2（t）=g3（t）+f30－f3e－β3t－r21－mev2（t）1－meu（t）+k2dt+σ3（t）dB3（t）（6）

滿足初始條件u0=lnx0，v10=lny10，v20=lny20。

顯然模型（6）的系數滿足局部Lipschitz條件，那么在t∈0，τe上存在唯一的局部解u（t），v1（t），v2（t），這里τe為爆破時刻。

由伊藤公式可得x（t），y1（t），y2（t）=eu（t），ev1（t），ev2（t）是模型（5）具有初始值x0，y10，y20∈R3+的唯一局部正解。

證明解是全局的，即證明τe=+SymboleB@。考慮如下隨機微分方程：

dΦ（t）=Φ（t）f1+f10－f1e－β1t－b1Φ（t）dt+Φ（t）σ1（t）dB1（t）（7）

dΨ1（t）=Ψ1（t）f2+f20－f2e－β2t－r11－mΨ1（t）1－mΦ（t）+k1dt+Ψ1（t）σ2（t）dB2（t）（8）

dψ1（t）=ψ1（t）f2+f20－f2e－β2 t－r11－mψ1（t）k1dt+ψ1（t）σ2（t）dB2（t）（9）

dΨ2（t）=Ψ2（t）f3+f30－f3e－β3 t－r21－mΨ2（t）1－mΦ（t）+k2dt+Ψ2（t）σ3（t）dB3（t）（10）

dψ2（t）=ψ2（t）f3+f30－f3e－β3 t－r21－mψ2（t）k2dt+ψ2（t）σ3（t）dB3（t）（11）

dφ（t）=φ（t）f1+f10－f1e－β1 t－b1φ（t）－a11－mΨ1（t）k1－a21－mΨ2（t）k2dt+

φ（t）σ1（t）dB1（t）（12）

其中，

Φ0=x0，Ψ10=y10，ψ10=y10

Ψ20=y20，ψ20=y20，φ0=x0

根據微分方程比較定理[18]，對于t∈0，τe有

φ（t）≤x（t）≤Φ（t），ψ1（t）≤y1（t）≤Ψ1（t），ψ2（t）≤y2（t）≤Ψ2（t）a.s.（13）

根據文獻[19]中的定理2.2，式（7）～（12）有解為

Φ（t）=exp∫t0 g1sds－f10－f1β1e－β1t－1+∫t0σ1sdB1sx－1（0）+b1∫t0exp∫s0g1τdτ－f10－f1β1e－β1s－1+∫s0σ1τdB1τds（14）

Ψ1（t）=exp∫t0 g2sds－f20－f2β2e－β2 t－1+∫t0σ2sdB2sy－11（0）+∫t0r11－m1－mΦs+k1exp∫s0 g2τdτ－f20－f2β2e－β2 s－1+∫s0σ2τdB2τds（15）

ψ1（t）=exp∫t0 g2sds－f20－f2β2e－β2 t－1+∫t0σ2sdB2sy－11（0）+r11－mk1∫t0exp∫s0g2τdτ－f20－f2β2e－β2s－1+∫s0σ2τdB2τds（16）

Ψ2（t）=exp∫t0g3sds－f30－f3β3e－β3t－1+∫t0σ3sdB3sy－12（0）+∫t0r21－m1－mΦs+k2exp∫s0g3τdτ－f30－f3β3e－β3s－1+∫s0σ3τdB3τds（17）

ψ2（t）=exp∫t0g3sds－f30－f3β3e－β3t－1+∫t0σ3sdB3sy－12（0）+r21－mk2∫t0exp∫s0g3τdτ－f30－f3β3e－β3s－1+∫s0σ3τdB3τds（18）

φ（t）=exp－∫t0G（s）－g1sds－f10－f1β1e－β1t－1+∫t0σ1sdB1sx－1（0）+b1∫t0exp－∫s0G（τ）－g1τdτ－f10－f1β1e－β1s－1+∫t0σ1τdB1τds（19）

其中，G（t）=a11－mΨ1（t）k1+a21－mΨ2（t）k2。

由于，φ（t），Φ（t），ψ1（t），Ψ1（t），ψ2（t），Ψ2（t）在0，+SymboleB@上全局存在，所以τe=+SymboleB@。

2 種群的持久與滅絕

定義1 令X（t）是模型（5）的解。1）若種群X滿足limt→+X（t）=0， a.s.，則種群X滅絕；2）若種群X滿足limt→+t－1∫t0X（s）ds=M， a.s.，則種群X平均持久生存；其中M為一個正常數。

引理1[20] 假設X（t）∈CΩ×R+，R0+，其中R0+：=a|agt;0，a∈R。

1）如果存在正常數λ0，T和λ≥0，滿足當t≥T時，

lnX（t）≤λT－λ0∫t0Xsds+∑ni=1βiBi（t）

則limt→+sup1t∫t0X（s）ds≤λλ0，a.s.。其中βi為常數，1≤i≤n。

2）如果存在正常數λ0，T和λ≥0，滿足當t≥T時，

lnX（t）≥λT－λ0∫t0Xsds+∑ni=1βiBi（t）

則limt→+sup1t∫t0X（s）ds≥λλ0，a.s.。其中βi為常數，1≤i≤n。

引理2 假設g-1gt;0，則有如下結論：1）如果g-2gt;0，那么limt→+t－1lny1（t）=0，a.s.；2）如果g-3gt;0，那么limt→+t－1lny2（t）=0，a.s.。

證明：不失一般性，本文僅證明g-2gt;0。令T充分大，當t≥T時有

g-i－εt≤∫t0gisds≤g-i+εt，" expg-i－εt≥2expg-i－εT

由式（14）有

Φ（t）=exp∫t0g1sds－f10－f1β1e－β1t－1+∫t0σ1sdB1sx－1（0）+b1∫t0exp∫s0g1τdτ－f10－f1β1e－β1s－1+∫s0σ1τdB1τds≤exp∫t0g1sds－f10－f1β1e－β1t－1+∫t0σ1sdB1sb1∫t0exp∫s0g1τdτ－f10－f1β1e－β1s－1+∫s0σ1τdB1τds≤expg-1+εt－f10－f1β1e－β1t－1+∫t0σ1sdB1sb1expmin0≤v≤t－f10－f1e－β1v－1β1+∫v0σ1τdB1τ∫tTexpg-1－εsds=g-1－εexpg-1+εt－f10－f1β1e－β1t－1+∫t0σ1sdB1sb1expg-1－εt－expg-1－εTexpmin0≤v≤t－f10－f1e－β1v－1β1+∫v0σ1τdB1τ

≤2g-1－εexpg-1+εt－f10－f1β1e－β1t－1+∫t0σ1sdB1sb1expg-1－εtexpmin0≤v≤t－f10－f1e－β1v－1β1+∫v0σ1τdB1τ=2g-1－εb1e2εtK1（t）

其中，

K1（t）=exp－f10－f1β1e－β1t－1+∫t0σ1sdB1sexpmin0≤v≤t－f10－f1e－β1v－1β1+∫v0σ1τdB1τ

由K1（t）≥1可得

∫tTr11－m1－mΦs+k1exp∫s0g2（τ）dτ－f20－f2e－β2s－1β2+∫s0σ2（τ）dB2（τ）ds≥∫tTr11－m1－m2g-1－εb1e2εtK1（t）+k1expg-2－εs－f20－f2e－β2s－1β2+∫s0σ2（τ）dB2（τ）ds≥∫tTr11－m1－m2g-1－εb1+k1e2εtK1（t）expg-2－εs－f20－f2e－β2s－1β2+∫s0σ2（τ）dB2（τ）ds=b1r11－m1－m2g-1－ε+b1k1∫tTK－11（s）expg-2－3εs－f20－f2e－β2s－1β2+∫s0σ2（τ）dB2（τ）ds≥b1r11－m21－mg-1－ε+b1k11g-2－3εexpg-2－3εt－expg-2－3εTmin0≤v≤tK2（v）=K3（t）expg-2－3εt－expg-2－3εT

其中，

K2（t）=K－11（t）exp－f20－f2e－β2s－1β2+∫t0σ2（τ）dB2（τ）K3（t）=b1r11－m21－mg-1－ε+b1k11g-2－3εmin0≤v≤tK2（v）

所以由式（15）有

1Ψ1（t）≥exp－∫tTg2（s）ds+f20－f2e－β2t－T－1β2－∫tTσ2（s）dB2（s）×K3（t）expg-2－3εt－expg-2－3εT≥exp－∫tTg2（s）ds+f20－f2e－β2t－T－1β2－∫tTσ2（s）dB2（s）×12K3（t）expg-2－3εt≥K4（t）exp（－4εt）

其中，

K4（t）=12K3（t）exp∫T0g2（s）ds+f20－f2e－β2t－T－1β2－∫tTσ2（s）dB2（s）（20）

所以有

t－1lnΨ1（t）lt;－t－1lnK4（t）+4ε（21）

由limt→+SymboleB@t－1∫t0σi（s）dBi（s）=0（i=1，2，3），結合式（20），如果g-2gt;0，則有limt→+SymboleB@t－1lnK4（t）=0，a.s.

結合式（13）和式（21）有limt→+SymboleB@supt－1lny1（t）≤limt→+SymboleB@supt－1lnΨ1（t）≤0，a.s.。

對式（9），應用伊藤公式可得

dlnψ1（t）=g2（t）+f20－f2e－β2t－r11－mψ1（t）k1dt+σ2（t）dB2（t）

對上式兩邊先積分再除以t可得

t－1lnψ1（t）=t－1lnψ10+t－1∫t0g2sds－f20－f2e－β2t－1β2t－r11－mk1t∫t0ψ1sds+t－1∫t0σ2sdB2（s）（22）

對任意的εgt;0，存在Tgt;0，使得當t≥T時有

g-2－ε2≤t－1lnψ10+t－1∫t0g2（s）ds－f20－f2e－β2t－1β2t≤g-2+ε2

由式（22），對于t≥T時有

t－1lnψ1（t）≤g-2+ε－r11－mk1t∫t0ψ1sds+t－1∫t0σ2sdB2（s）（23）

t－1lnψ1（t）≥g-2－ε－r11－mk1t∫t0ψ1sds+t－1∫t0σ2sdB2（s）（24）

選取充分小的ε滿足g-2－εgt;0，再利用引理1中的1）和2）有

g-2－εk1r11－m≤limt→+SymboleB@inft－1∫t0ψ1sds≤limt→+SymboleB@supt－1∫t0ψ1sds≤g-2+εk1r11－m，a.s.

由ε的任意性可得

limt→+SymboleB@t－1∫t0ψ1sds=g-2k1r11－m，a.s.（25）

這說明了limt→+SymboleB@t－1lnψ10=0，a.s.。由式（13）有

limt→+SymboleB@inft－1lny1（t）≥limt→+SymboleB@t－1lnΨ1（t）≤0，a.s.（26）

定理2 對模型（5）：

1）如果g-1lt;0，g-2lt;0和g-3lt;0，那么獵物種群x和捕食者種群y1，y2均滅絕，即

limt→+SymboleB@x（t）=0，limt→+SymboleB@y1（t）=0，limt→+SymboleB@y2（t）=0，a.s.

2）如果g-1lt;0，g-2gt;0和g-3lt;0，那么獵物種群x和捕食者種群y2均滅絕，而捕食者種群y1平均持久生存，即

limt→+SymboleB@t－1∫t0y1sds=k1g-2r11－m，a.s.

3）如果g-1lt;0，g-2lt;0和g-3gt;0，那么獵物種群x和捕食者種群y1均滅絕，而捕食者種群y2平均持久生存，即

limt→+SymboleB@t－1∫t0y2sds=k2g-3r21－m，a.s.

4）如果g-1lt;0，g-2gt;0和g-3gt;0，那么獵物種群x滅絕，而捕食者種群y1，y2平均持久生存，即

limt→+SymboleB@t－1∫t0y1sds=k1g-2r11－m，limt→+SymboleB@t－1∫t0y2sds=k2g-3r21－m，a.s.（27）

5）如果g-1gt;0，g-2lt;0和g-3lt;0，那么捕食者種群y1，y2均滅絕，而獵物種群x平均持久生存，即

limt→+SymboleB@t－1∫t0xsds=g-1b1a.s.（28）

6）如果g-1gt;0，g-2gt;0和g-3lt;0，那么捕食者種群y2滅絕，此外，

（1）如果g-1lt;a1g-2r1，那么獵物種群x滅絕，而捕食者種群y1平均持久生存，即

limt→+SymboleB@t－1∫t0y1sds=k1g-2r11－m

（2）如果g-1gt;a1g-2r1，那么

limt→+SymboleB@t－1∫t0xsds=r1g-1－a1g-2r1b1，limt→+SymboleB@t－1∫t0y1s1－mxs+k1ds=g-2r1（1－m），a.s.

7）如果g-1gt;0，g-2lt;0和g-3gt;0，那么捕食者種群y1滅絕，此外，

（1）如果g-1lt;a2g-3r2，那么獵物種群x滅絕，而捕食者種群y2平均持久生存，即

limt→+SymboleB@t－1∫t0y2sds=k2g-3r21－m

（2）如果g-1gt;a2g-3r2，那么

limt→+SymboleB@t－1∫t0xsds=r2g-1－a2g-3r2b1，limt→+SymboleB@t－1∫t0y2s1－mxs+k2ds=g-3r2（1－m），a.s.

8）如果g-1gt;0，g-2gt;0和g-3gt;0，那么

（1）如果g-1lt;a1g-2r1+a2g-3r2，那么獵物種群x滅絕，而捕食者種群y1，y2平均持久生存，即

limt→+SymboleB@t－1∫t0y1sds=k1g-2r11－m，limt→+SymboleB@t－1∫t0y2sds=k2g-3r21－m

（2）如果g-1gt;a1g-2r1+a2g-3r2，那么

limt→+SymboleB@t－1∫t0xsds=r2g-1－a2g-3r1－a1r2g-2r1r2b1，a.s.，

limt→+SymboleB@t－1∫t0y1s1－mxs+k1ds=g-2r1（1－m），a.s.，limt→+SymboleB@t－1∫t0y2s1－mxs+k2ds=g-3r2（1－m），a.s.（29）

證明：對模型（5）應用伊藤公式可得

dlnx（t）=g1（t）+f10－f1e－β1t－b1x（t）－a11－my1（t）1－mx（t）+k1－a21－my2（t）1－mx（t）+k2dt+

σ1（t）dB1（t）dlny1（t）=g2（t）+f20－f2e－β2t－r11－my1（t）1－mx（t）+k1dt+σ2（t）dB2（t）dlny2（t）=g3（t）+f30－f3e－β3t－r21－my2（t）1－mx（t）+k2dt+σ3（t）dB3（t）

對上式進行積分可得

lnx（t）－lnx0=∫t0g1（s）ds－f10－f1e－β1t－1β1－b1∫t0x（s）ds－

a11－m∫t0y1s1－mxs+k1ds－a21－m∫t0y2s1－mxs+k2ds+∫t0σ1（s）dB1（s）（30）

lny1（t）－lny10=∫t0g2（s）ds－f20－f2e－β2t－1β2－r11－m∫t0y1s1－mxs+k1ds+∫t0σ2（s）dB2（s），（31）

lny2（t）－lny20=∫t0g3（s）ds－f30－f3e－β3t－1β3－r21－m∫t0y2s1－mxs+k2ds+∫t0σ3（s）dB3（s）（32）

1）對于足夠大的t，由式（30）有

t－1lnx（t）x0≤g-1+ε+t－1∫t0σ1（s）dB1（s）－f10－f1e－β1t－1β1t

由于limt→+SymboleB@t－1∫t0σ1（s）dB1（s）=0和g-1+εlt;0，那么limt→+SymboleB@x（t）=0，a.s.。

同理，如果g-2lt;0，那么limt→+SymboleB@y1（t）=0，a.s.；如果g-3lt;0，那么limt→+SymboleB@y2（t）=0，a.s.。

2）注意到g-1lt;0，g-3lt;0，由1）的證明可知limt→+SymboleB@x（t）=0，limt→+SymboleB@y2（t）=0， a.s.；對于足夠大的t，由式（31）可得

lny1（t）－lny10≤g-2+εt－r11－m1－mε+k1∫t0y1sds+∫t0σ2（s）dB2（s）（33）

lny1（t）－lny10≥g-2－εt－r11－mk1－1－mε∫t0y1sds+∫t0σ2（s）dB2（s）（34）

對式（33）、式（34）分別應用引理1中的1），2）可得到

limt→+SymboleB@supt－1∫t0y1sds≤g-2+ε1－mε+k1r11－m，a.s.

limt→+SymboleB@inft－1∫t0y1sds≥g-2－εk1－1－mεr11－m a.s.

由ε的任意性，可得到limt→+SymboleB@t－1∫t0y1sds=k1g-2r11－m，a.s.。

3）證明類似于2），所以省略。

4）由于g-1lt;0，從而由1）的結論可知limt→+SymboleB@x（t）=0，a.s.；式（27）的證明類似于2），所以省略。

5）因為g-2lt;0，g-3lt;0，所以由1）的結論可知limt→+SymboleB@y1（t）=0，limt→+SymboleB@y2（t）=0。式（28）的證明類似于2），所以省略。

6）因為g-3lt;0，所以由1）的結論知limt→+SymboleB@y2（t）=0。

（1）由式（30）×r1-式（31）×a1可得

t－1r1lnx（t）x0=t－1a1lny1（t）y10+t－1r1∫t0g1sds－t－1a1∫t0g2sds－f10－f1r1e－β1t－1β1t+

f20－f2a1e－β2t－1β2t－t－1b1r1∫t0xsds－t－1a21－mr1∫t0y2s1－mxs+k2ds+

t－1r1∫t0σ1（s）dB1（s）－t－1a1∫t0σ2（s）dB2（s）≤t－1a1lny1（t）y10+t－1r1∫t0g1sds－t－1a1∫t0g2sds－f10－f1r1e－β1t－1β1t+

f20－f2a1e－β2t－1β2t－t－1b1r1∫t0xsds+t－1r1∫t0σ1（s）dB1（s）－t－1a1∫t0σ2（s）dB2（s）（35）

由引理2，對任意的εgt;0，存在Tgt;0，使得當t≥T時，

a1tlny1（t）y10≤ε4，r1lnx0t≤ε4，－f10－f1r1e－β1t－1β1t+f20－f2a1e－β2t－1β2t≤ε4，r1t∫t0σ1（s）dB1（s）－a1t∫t0σ2（s）dB2（s）≤ε4，t－1r1∫t0g1sds－t－1a1∫t0g2sds≤r1g-1－a1g-2+a1ε+r1ε

將以上各不等式帶入式（35）可得

t－1r1lnx（t）≤a1+r1+1ε+r1g-1－a1g-2（36）

令ε充分小，使得0lt;εlt;a1g-2－r1g-1a1+r1+1，所以limt→+SymboleB@x（t）=0， a.s.；limt→+SymboleB@t－1∫t0y1sds=k1g-2r11－m的證明類似于2），所以省略，（1）的證明完畢。

（2）由式（31）可得

lny1（t）t－lny10t=t－1∫t0g2sds－f20－f2e－β2t－1β2t－

t－11－mr1∫t0y1s1－mxs+k1ds+t－1∫t0σ2（s）dB2（s）

由引理2和limt→+SymboleB@t－1∫t0σ2（s）dB2（s）=0可得

limt→+SymboleB@t－1∫t0y1s1－mxs+k1ds=g-2r11－m（37）

由式（30）可得

lnx（t）t=∫t0g1sdst－f10－f1e－β1t－1β1t+lnx0t－b1∫t0xsdst－

a11－mt∫t0y1s1－mxs+k1ds－a21－mt∫t0y2s1－mxs+k2ds+t－1∫t0σ1（s）dB1（s）（38）

所以對任意的εgt;0，存在Tgt;0，使得對t≥T，由limt→+SymboleB@y2（t）=0和式（37）有

－a1g-2r1－ε≤－f10－f1e－β1t－1β1t+lnx0t

－a11－mt∫t0y1s1－mxs+k1ds－

a21－mt∫t0y2s1－mxs+k2ds≤－a1g-2r1+ε

將不等式帶入式（38），對于t≥T可得

lnx（t）t≥g-1－a1g-2r1－ε－b1∫t0xsdst+t－1∫t0σ1（s）dB1（s）（39）

lnx（t）t≤g-1－a1g-2r1+ε－b1∫t0xsdst+t－1∫t0σ1（s）dB1（s）（40）

選取充分小的ε以滿足0lt;εlt;g-1－a1g-2r12，對式（39）、式（40）分別使用引理1中的1），2）可得：

g-1－εr1－a1g-2r1b1≤limt→SymboleB@inft－1∫t0xsds≤limt→SymboleB@supt－1∫t0xsds≤g-1+εr1－a1g-2r1b1，a.s.

由ε的任意性，有limt→+SymboleB@t－1∫t0xsds=r1g-1－a1g-2r1b1，a.s.，6）的證明完畢。

7）證明類似于6），所以省略。

8）證明

（1）由式（35）×r2-式（32）×a2r1可得

r1r2tlnx（t）x0=a1r2tlny1（t）y10+a2r1tlny2（t）y20+r1r2t∫t0g1（s）ds－a1r2t∫t0g2（s）ds－a2r1t∫t0g3（s）ds－f10－f1r1r2e－β1t－1β1t+a1e－β2t－1f20－f2r2β2t+f30－f3a2r1e－β3t－1β3t－t－1b1r1r2∫t0xsds+r1r2t∫t0σ1（s）dB1（s）－a1r2t∫t0σ2（s）dB2（s）－a2r1t∫t0σ3（s）dB3（s）（41）

由引理2可知，對任意的εgt;0，存在Tgt;0，使得當t≥T時有

a1r2tlny1（t）y10+a2r1tlny2（t）y20≤ε4，r1r2lnx0t≤ε4，－f10－f1r1r2e－β1t－1β1t+a1e－β2t－1f20－f2r2β2t+f30－f3a2r1e－β3t－1β3t≤ε4，r1r2t∫t0σ1（s）dB1（s）－a1r2t∫t0σ2（s）dB2（s）－a2r1t∫t0σ3（s）dB3（s）≤ε4，r1r2t∫t0g1（s）ds－a1r2t∫t0g2（s）ds－a2r1t∫t0g3（s）ds≤r1r2g-1－a1r2g-2－a2r1g-3+a2r1+r1r2+a1r2ε

將上述不等式帶入式（41）可得

lnx（t）r1r2t≤r1r2g-1－a1r2g-2－a2r1g-3+a2r1+r1r2+a1r2+1ε（42）

選取充分小的ε，使得0lt;εlt;a1r2g-2+a2r1g-3－r1r2g-1a2r1+r1r2+a1r2+1，由式（42）可知limt→+SymboleB@x（t）=0；

limt→+SymboleB@t－1∫t0y1sds=k1g-2r11－m，limt→+SymboleB@t－1∫t0y2sds=k2g-3r21－m的證明類似于2），所以省略。

（2）對式（31），式（32）利用引理2和limt→+SymboleB@t－1∫t0σi（s）dBi（s）=0， i=1，2，3，有

limt→+SymboleB@1t∫t0y1s1－mxs+k1ds=g-2r11－m，limt→+SymboleB@1t∫t0y2s1－mxs+k2ds=g-3r21－m，a.s.

那么對任意的εgt;0，存在Tgt;0，使得當t≥T時有

－g-3a2r2－a1g-2r1－ε≤－f10－f1e－β1t－1β1t－a21－mt∫t0y2s1－mxs+k2ds－a11－mt∫t0y1s1－mxs+k1ds+lnx0t≤－g-3a2r2－a1g-2r1+ε（43）

將式（43）帶入式（38）可得

t－1lnx（t）≥g-1－g-3a2r2－a1g-2r1－ε－t－1b1∫t0xsds+t－1∫t0σ1（s）dB1（s）t－1lnx（t）≤g-1－g-3a2r2－a1g-2r1+ε－t－1b1∫t0xsds+t－1∫t0σ1（s）dB1（s）

選取充分小的ε以滿足0lt;εlt;a1r2g-2+a2r1g-3－r1r2g-12r1r2，由引理1可得

g-1－εr2－a2g-3r1－a1r2g-2b1r1r2≤limt→+SymboleB@inft－1∫t0xsds≤limt→+SymboleB@supt－1∫t0xsds≤g-1+εr2－a2g-3r1－a1r2g-2b1r1r2

由ε的任意性可得出limt→+SymboleB@t－1∫t0xsds=－a2g-3+r2g-1r1－a1r2g-2b1r1r2。

3 數值模擬

本節使用python進行數值仿真實驗，以進一步說明理論結果的正確性，利用文獻[21]中的Milstein方法對模型（5）進行數值求解，則模型（5）對應的離散化形式為

xi+1=xi+xif1+f10－f1e－β1kΔt－b1xi－a11－my1i1－mxi+k1－a21－my2i1－mxi+k2Δt+

xiη12β11－e－2β1kΔtΔtξi1+xiη211－e－2β1kΔt4β1ξi12－1Δt

y1i+1=y1i+y1if2+f20－f2e－β2kΔt－r11－my1i1－mxi+k1Δt+

y1iη22β21－e－2β2kΔtΔtξi2+y1iη221－e－2β2kΔt4β2ξi22－1Δt，y2i+1=y2i+y2if3+f30－f3e－β3kΔt－r21－my2i1－mxi+k2Δt+y2iη32β31－e－2β3kΔtΔtξi3+y2iη231－e－2β3kΔt4β3ξi32－1Δt

其中，ξi1，ξi2和ξi3i=1，2，3，…，n是具有分布N（0，1）的隨機變量且相互獨立，Δt代表時間步長。

本節參照文獻[11]，[16]對各參數進行取值，不失一般性，本文僅討論波動強度ηi對種群動態行為的影響。

1） 在圖1中，令x（0）=0.4，y1（0）=0.6，y2（0）=0.2，f1=0.5，f2=0.5，f3=0.3，f10=0.3，f20=0.25，f30=0.2，b1=0.25，a1=0.36，a2=0.4，r1=0.5，r2=0.47，k1=k2=1，β1=0.31，β2=0.35，β3=0.42，m=0.1。圖1a～d僅η2ii=1，2，3取值不同，其余參數取值一樣。

（1）取η21=0.7，η22=0.9，η23=0.6，通過計算可得g-1=－0.065lt;0，g-2=－0.143lt;0，g-3=－0.057lt;0，因此滿足定理21）中的參數條件，由理論結果可知，所有種群均滅絕。數值模擬結果見圖1a。

（2）取η21=0.7，η22=0.3，η23=0.6，通過計算可得g-1=－0.065lt;0，g-2=0.286gt;0，g-3=－0.057lt;0，因此滿足定理22）中的參數條件，由理論結果可知，獵物種群x和捕食者種群y2均滅絕，而捕食者種群y1平均持久生存，即

limt→+SymboleB@t－1∫t0y1sds=k1g-2r11－m=0.636

數值模擬結果如圖1b所示。

（3）取η21=0.3，η22=0.8，η23=0.6，通過計算可得g-1=0.258gt;0，g-2=－0.071lt;0，g-3=－0.057lt;0，

因此滿足定理25）中的參數條件，由理論結果可知，捕食者種群y1，y2均滅絕，而獵物種群x平均持久生存，即

limt→+SymboleB@t－1∫t0xsds=g-1b1=1.032

數值模擬結果如圖1c所示。

（4）取η21=0.7，η22=0.3，η23=0.3，通過計算可得g-1=－0.065lt;0，g-2=0.286gt;0，g-3=0.121gt;0，

因此滿足定理24）中的參數條件，由理論結果可知，獵物種群x滅絕，而捕食者種群y1，y2平均持久生存，即

limt→+SymboleB@t－1∫t0y1sds=k1g-2r11－m=0.636，limt→+SymboleB@t－1∫t0y2sds=k2g-3r21－m=0.286，a.s.

數值模擬結果如圖1d所示。

比較圖1a與圖1b，可以發現僅對波動強度η2的大小進行改變，而其他參數不變時，圖像發生了變化，當η2縮小，捕食者種群y1由原來的滅絕狀態變為平均持久生存狀態；比較圖1a與圖1c，當η1和η2縮小，其他參數不改變時，獵物種群x由原來的滅絕狀態變為平均持久生存狀態；比較圖1a與圖1d，當η2和η3縮小，其他參數也不改變時，捕食者種群y1，y2由原來的滅絕狀態變為平均持久生存狀態。

2） 在圖2中，令x（0）=0.4，y1（0）=0.6，y2（0）=0.2，f1=0.6，f2=0.8，f3=0.3，f10=0.3，f20=0.25，f30=0.2，b1=0.25，a1=0.36，a2=0.4，r1=0.5，r2=0.47，k1=k2=1，β1=0.81，β2=0.65，β3=0.42，m=0.1。圖2a，b僅η2ii=1，2，3取值不同，其余參數取值一樣。

（1）取η21=0.5，η22=0.06，η23=0.7，通過計算可得g-1=0.446gt;0，g-2=0.777gt;0，g-3=－0.117lt;0，因此滿足定理2-6）-（1）中的參數條件，由理論結果可知，捕食者種群y2滅絕，由g-1lt;a1g-2r1=0.559知，獵物種群x滅絕，而捕食者種群y1平均持久生存，即

limt→+SymboleB@t－1∫t0y1sds=k1g-2r11－m=1.727

數值模擬結果如圖2a所示。

（2）取η21=0.3，η22=0.56，η23=0.7，通過計算可得g-1=0.507gt;0，g-2=0.585gt;0，g-3=－0.117lt;0，因此滿足定理26）（2）中的參數條件，由理論結果可知，捕食者種群y2滅絕，由g-1gt;a1g-2r1=0.421知，獵物種群x和捕食者種群y1平均持久生存，即

limt→+SymboleB@t－1∫t0xsds=r1g-1－a1g-2r1b1=0.343，limt→+SymboleB@t－1∫t0y1s1－mxs+k1ds=g-2r1（1－m）=1.3

數值模擬結果如圖2b所示。

比較圖2a和b，當縮小η1，放大η2，而其他參數不改變時，獵物種群x由原來的滅絕狀態轉變為平均持久生存狀態。

3） 在圖3中，令x（0）=0.4，y1（0）=0.6，y2（0）=0.2，f1=0.8，f2=0.25，f3=0.1，f10=0.4，f20=0.05，f30=0.05，b1=0.4，a1=0.8，a2=0.6，r1=0.4，r2=0.3，k1=k2=1，β1=0.61，β2=0.85，β3=0.89，m=0.1。

圖3a，b僅η2ii=1，2，3取值不同，其余參數取值一樣。

（1）取η21=0.5，η22=0.006，η23=0.009，通過計算可得g-1=0.595gt;0，g-2=0.248gt;0，g-3=0.097gt;0，因此滿足定理28）（1）中的參數條件，由g-1lt;a1g-2r1+a2g-3r2=0.690知，獵物種群x滅絕，而捕食者種群y1，y2平均持久生存，即

limt→+SymboleB@t－1∫t0y1sds=k1g-2r11－m=0.689，limt→+SymboleB@t－1∫t0y2sds=k2g-3r21－m=0.359

數值模擬結果如圖3a所示。

（2）取η21=0.000 2，η22=0.006，η23=0.009，通過計算可得g-1=0.80gt;0，g-2=0.248gt;0，g-3=0.097gt;0，因此滿足定理28）（1）中的參數條件，由g-1gt;a1g-2r1+a2g-3r2=0.690知，獵物種群x和捕食者種群y1，y2平均持久生存，即

limt→+SymboleB@t－1∫t0xsds=r2g-1－a2g-3r1－a1r2g-2r1r2b1=0.275limt→+SymboleB@t－1∫t0y1s1－mxs+k1ds=g-2r1（1－m）=0.689

limt→+SymboleB@t－1∫t0y2s1－mxs+k2ds=g-3r2（1－m）=0.359

數值模擬結果如圖3b所示。

比較圖3a和b，當對η1進行大幅度縮小時，獵物種群x由原來的滅絕狀態轉變為平均持久生存狀態，此時3個種群都為平均持久生存狀態。通過對圖1～圖3中的圖片分別進行比較，可以發現，當改變波動強度ηi的大小時，種群的持久生存和滅絕狀態也會發生改變。

4） 在圖4中，令x（0）=0.4，y1（0）=0.6，y2（0）=0.2，f1=0.3，f2=0.2，f3=0.3，f10=0.2，f20=0.05，f30=0.14，b1=0.6，a1=0.3，a2=0.3，r1=0.96，r2=0.98，k1=k2=1，β1=0.82，β2=0.87，β3=0.9，η21=0.03，η22=0.007，η23=0.009。圖4a～c僅避難所系數m取值不同，其余參數取值一樣。

選取避難所系數m=0，0.4，0.8，容易驗證選取的參數值均滿足定理28）（6）的參數條件，因此所有種群平均持久生存。由圖4可知，隨著獵物避難所的增加，獵物種群的數量有所上升，而捕食者種群的數量有一定程度的下降，說明一定數量的獵物避難所有利于獵物種群的生存，但不利于捕食者種群的生存。

4 結論

本文在模型（1）的基礎上建立了一個具有獵物避難所和Ornstein-Uhlenbeck過程的隨機三種群捕食模型。證明了模型（5）全局正解的存在唯一性，得到了模型（5）中種群平均持久生存和滅絕的充分條件，揭示了Ornstein-Uhlenbeck過程對種群平均持久生存和滅絕的影響：1）對g-i=fi－η2i4βii=1，2，3，定理2中的1）說明如果回歸速率βi相對較小，而波動強度ηi相對較大，那么無論是捕食者種群還是獵物種群都會滅絕。

2）定理2中的2）說明βi（i=1，3）相對較小，ηi（i=1，3）相對較大時，獵物種群x和捕食者種群y2均滅絕；當βi（i=2）相對較大，ηi（i=2）相對較小時，捕食者種群y1平均持久生存。

3）定理2中的3）說明βi（i=1，2）相對較小，ηi（i=1，2）相對較大時，獵物種群x和捕食者種群y1均滅絕；當βi（i=3）相對較大，ηi（i=3）相對較小時，捕食者種群y2平均持久生存。

4）定理2中的4）說明βi（i=2，3）相對較大，ηi（i=2，3）相對較小時，捕食者種群y1，y2平均持久生存；當βi（i=1）相對較小，ηi（i=1）相對較大時，獵物種群x滅絕。

5）定理2中的5）說明βi（i=2，3）相對較小，ηi（i=2，3）相對較大時，捕食者種群y1，y2均滅絕；當βi（i=1）相對較大，ηi（i=1）相對較小時，獵物種群x平均持久生存。

6）定理2中的6）說明在g-i（i=1，2）適當大的條件下，（1）表明在r1足夠小，a1足夠大時，獵物種群x滅絕，而捕食者種群y1平均持久生存；（2）表明在r1足夠大，a1足夠小時，獵物種群x和捕食者種群y1平均持久生存。當βi（i=3）相對較小，ηi（i=3）相對較大時，捕食者種群y2滅絕。

7）定理2中的7）說明在g-i（i=1，3）適當大的條件下，（1）表明在r2足夠小，a2足夠大時，獵物種群x滅絕，而捕食者種群y2平均持久生存；（2）表明在r2足夠大，a2足夠小時，獵物種群x和捕食者種群y2平均持久生存。當βi（i=2）相對較小，ηi（i=2）相對較大時，捕食者種群y1滅絕。

8）定理2中的8）說明在g-i（i=1，2，3）適當大的條件下，（1）表明在r1，r2足夠小，a1，a2足夠大時，獵物種群x滅絕，而捕食者種群y1，y2平均持久生存；（2）表明在r1，r2足夠大，a1，a2足夠小時，獵物種群x和捕食者種群y1，y2平均持久生存。

最后利用數值仿真驗證了理論結果，并且也模擬了獵物避難所對種群動態行為的影響。

參考文獻：

[1]BERRYMAN A A. The orgins and evolution of predator-prey theory[J]. Ecology， 1992， 73（5）： 15301535.

[2]HAQUE M. A detailed study of the Beddington-DeAngelis predator-prey model[J]. Mathematical Biosciences， 2011， 234（1）： 116.

[3]XIE Y， WANG Z， MENG B， et al. Dynamical analysis for a fractional-order prey-predator model with Holling III type functional response and discontinuous harvest[J]. Applied Mathematics Letters， 2020， 106： 106342.

[4]LUO Y， ZHANG L， TENG Z， et al. Global stability for a nonautonomous reaction-diffusion predator-prey model with modified Leslie-Gower Holling-II schemes and a prey refuge[J]. Advances in Difference Equations， 2020， 2020（1）： 116.

[5]YUE Q. Dynamics of a modified Leslie-Gower predator-prey model with Holling-type II schemes and a prey refuge[J]. SpringerPlus， 2016， 5（1）： 112.

[6]ZHANG H， CAI Y， FU S， et al. Impact of the fear effect in a prey-predator model incorporating a prey refuge[J]. Applied Mathematics and Computation， 2019， 356： 328337.

[7]QI H， MENG X. Threshold behavior of a stochastic predator-prey system with prey refuge and fear effect[J]. Applied Mathematics Letters， 2021， 113： 106846.

[8]DAS S， MAHATO P， MAHATO S K. Disease control prey-predator model incorporating prey refuge under fuzzy uncertainty[J]. Modeling Earth Systems and Environment， 2021， 7（4）： 21492166.

[9]MUKHERJEE D. The effect of refuge and immigration in a predator-prey system in the presence of a competitor for the prey[J]. Nonlinear Analysis： Real World Applications， 2016， 31： 277287.

[10] MUKHERJEE D， MAJI C. Bifurcation analysis of a Holling type II predator-prey model with refuge[J]. Chinese Journal of Physics， 2020， 65： 153162.

[11] XU Y， LIU M， YANG Y. Analysis of a stochastic two-predators one-prey system with modified Leslie-Gower and Holling-typeⅡ schemes[J]. Journal of Applied Analysis amp; Computation， 2017， 7（2）： 713727.

[12] TIAN B， YANG L， CHEN X， et al. A generalized stochastic competitive system with Ornstein-Uhlenbeck process[J]. International Journal of Biomathematics， 2021， 14（1）： 2150001.

[13] TONG Z， LIU A. A censored Ornstein-Uhlenbeck process for rainfall modeling and derivatives pricing[J]. Physica A： Statistical Mechanics and its Applications， 2021， 566： 125619.

[14] WANG W， CAI Y， DING Z， et al. A stochastic differential equation SIS epidemic model incorporating Ornstein-Uhlenbeck process[J]. Physica A： Statistical Mechanics and its Applications， 2018， 509： 921936.

[15] ZHOU D， LIU M， LIU Z. Persistence and extinction of a stochastic predator-prey model with modified Leslie-Gower and Holling-type II schemes[J]. Advances in Difference Equations， 2020， 2020（1）： 115.

[16] GAO Y， YAO S. Persistence and extinction of a modified Leslie-Gower Holling-typeⅡ predator-prey stochastic model in polluted environments with impulsive toxicant input[J]. Mathematical Biosciences and Engineering， 2021， 18（4）： 48944918.

[17] MAO X. Stochastic Differential Equations and Applications[M]. Chichester（UK）：Horwood Publishing， 2007.

[18] IKEDA N， WATANABE S. Stochastic Differential Equations and Diffusion Processes[M]. Amsterdam： North-Holland： Elsevier， 2014.

[19] JIANG D， SHI N. A note on nonautonomous logistic equation with random perturbation[J]. Journal of Mathematical Analysis and Applications， 2005， 303（1）： 164172.

[20] LIU M， WANG K， WU Q. Survival analysis of stochastic competitive models in a polluted environment and stochastic competitive exclusion principle[J]. Bulletin of mathematical biology， 2011， 73（9）： 19692012.

[21] HIGHAM D J. An algorithmic introduction to numerical simulation of stochastic differential equations[J]. SIAM Review， 2001， 43（3）： 525546.

（責任編輯 李 進）