摘要：本文建立了一類具有線性收獲率及Holling II型功能反應(yīng)函數(shù)的共享資源型捕食-食餌模型，分析了該系統(tǒng)平衡點的存在性與穩(wěn)定性。系統(tǒng)可能存在無食餌種群的平衡點，無捕食者種群的平衡點以及正平衡點，通過定性分析得到無食餌種群的平衡點在不同的條件下可以是穩(wěn)定的平衡點、鞍點或者鞍結(jié)點，而無捕食者種群的平衡點始終是鞍點，并給出了正平衡點穩(wěn)定的條件。最后，對系統(tǒng)進行數(shù)值模擬，發(fā)現(xiàn)收獲率的不同取值可以改變系統(tǒng)的穩(wěn)定狀態(tài)。因此，利用理論分析結(jié)果可得到在保證生態(tài)穩(wěn)定的前提下人類捕撈的最優(yōu)收獲。

關(guān)鍵詞：收獲率；捕食-食餌模型；Holling II型功能反應(yīng)函數(shù)；平衡點

中圖分類號： O175.1文獻標志碼：A文獻標識碼

Establishment and analysis of a shared resource predator-prey model with harvest

ZHANG "Xueting，LI "Guihua*，LUO "Hang

（School of Mathematics， North University of China， Taiyuan，Shanxi 030051，China）

Abstract： "In this paper， a shared resource predator-prey model with linear harvesting rate and Holling type II functional response function is established， and the existence and stability of the equilibria of the system are analyzed. There may exist the prey-free equilibrium， the predator-free equilibrium， and the positive equilibrium in the system. Through qualitative analysis， we give the conditions that the prey-free equilibrium is stable， saddle， or saddle-node one， respectively， and find that the predator-free equilibrium always be saddle one. Furthermore， we obtain the conditions that the positive equilibrium is stable. Finally， some numerical simulations are carried out. It was found that different values of harvesting rate can change the stable state of the system. Hence， the optimal harvest for human fishing can be obtained while ensuring ecological stability using theoretical analysis.

Key words： harvesting rate;predator-prey model;Holling type II functional response function;equilibrium

1引言

種群生態(tài)學是生態(tài)學的一個重要分支，起源于應(yīng)用昆蟲學、人口統(tǒng)計學和水產(chǎn)資源學3個學科。近年來，由于傳統(tǒng)魚類資源過度開發(fā)，人們對于從海洋中獲取新型食物的興趣也日益濃厚，造成了由于過度捕撈而導致物種滅絕的情況。例如分布在我國長江流域的白鱘，由于上世紀70年代人類的大肆捕撈，最終在2019年12月正式宣告滅絕。為了在保護生物多樣性的同時也滿足人類發(fā)展的需要，越來越多的學者開始研究數(shù)學生態(tài)模型。

早在1925年，研究海洋生態(tài)系統(tǒng)的數(shù)學家就開始了對掠肉魚和食用魚之間捕食者與食餌關(guān)系的研究，Alfred J. Lotka和Vito Volterra隨之建立了著名的Lotka-Volterra（L-V）模型，之后關(guān)于捕食者-食餌系統(tǒng)的研究便進入了蓬勃發(fā)展的階段。在眾多改進的L-V模型中，應(yīng)用非常廣泛的捕食模型還有Logistic模型、Leslie模型等。在這些經(jīng)典的模型中，環(huán)境承載能力都假設(shè)為一個固定的值，考慮到現(xiàn)實的可變環(huán)境，且捕食者和食餌往往會共享一些資源。在2013年Safuan等［1］提出了一種捕食者和食餌共享資源的三維生態(tài)模型，其中環(huán)境承載能力與生物資源成比例，捕食者和食餌種群都可以通過與生物資源相互作用直接改變可用資源的量，模型如下：

dxdt=r1x1-xpz-axy，

dydt=r2y1-yqz+bxy，

dzdt=z（c-dx-ey），（1）

其中x表示食餌種群的密度，y表示捕食者種群的密度，z表示生物資源的量。r1，r2分別表示食餌種群和捕食者種群的內(nèi)稟增長率，生物資源以恒定速率c增長，常數(shù)d和e分別表示食餌和捕食者對生物資源的消耗率，pz，qz分別代表了食餌種群和捕食者種群的環(huán)境承載能力，其中p，qgt;0，且p+q≤1。若pgt;q，意味著食餌可以獲得更高比例的生物資源。ax為功能反應(yīng)函數(shù)，是反映捕食者與食餌之間能量傳遞的量。作者通過理論分析計算了平衡點的存在性，通過數(shù)值模擬發(fā)現(xiàn)系統(tǒng)可能發(fā)生Hopf分支。1965年生物數(shù)學家Holling曾經(jīng)在試驗的基礎(chǔ)上，對不同物種提出了3種不同的功能反應(yīng)函數(shù)F（x）。其中2種分別為：（1）F（x）=ax（Holling I型），適合海藻類、浮游動物等低等生物群；（2）F（x）=ax1+wx（Holling II型），主要適用于無脊椎動物。Capone等［2］引入了Holling II型功能反應(yīng)函數(shù)，推廣了環(huán)境承載能力與生物資源成正比的三維模型。

隨著漁業(yè)的不斷發(fā)展，如何既使生物種群長久持續(xù)生存，又能使人們獲得可觀的經(jīng)濟利潤成為了一個廣泛關(guān)注的課題。為了幫助人們更好地掌握捕獲量，防止過度捕撈，在模型中引入和考慮種群收獲至關(guān)重要。目前，常數(shù)收獲項、線性收獲項是較為主流的收獲項，另外還有更復雜的Michaelis-Menten（M-M）型收獲項。常數(shù)收獲適用于已經(jīng)高度開發(fā)或容易受到過度捕撈的漁業(yè)資源；而比例收獲和M-M型收獲適用于較為穩(wěn)定或能夠適應(yīng)一定程度捕撈的漁業(yè)資源。Ang等［3］在具有可變環(huán)境承載力的模型（1）中，引入了對捕食者的非線性M-M型收獲，討論了模型的動力學性態(tài)和最優(yōu)收獲。Ganguli等［4］在模型（1）的基礎(chǔ)上，引入了對捕食者和食餌的線性收獲，分析了該模型的動力學行為并求解了最優(yōu)采收政策。

在本文中，考慮一種食餌種群為海洋生態(tài)系統(tǒng)中的無脊椎動物，選擇了更符合漁業(yè)中生物種群動態(tài)的Holling II型功能反應(yīng)函數(shù)，再考慮到人類的捕撈行為，加入對食餌種群的線性收獲。模型建立如下：

dxdt=r1x1-xpz-axy1+wx-hx，

dydt=r2y1-yqz+bxy1+wx，

dzdt=zc-dx-ey，（2）

其中a，b，w都是常數(shù)，h表示恒定比例的收獲率。下面對這個模型進行動力學性態(tài)分析。

2平衡點的存在性

為了簡化計算，對系統(tǒng)（2）進行無量綱化，令=aqxpr2，=ayr2，=aqzr2。為了方便，仍用x，y，z代替，，，則系統(tǒng)變化為

dxdt=r1x1-xz-r2xy1+αx-hx，

dydt=r2y1-yz+βxy1+αx，

dzdt=z（c-δx-σy），（3）

其中α=wpr2aq，β=bpr2aq，δ=dpr2aq，σ=er2a。

為了分析系統(tǒng)（3）平衡點的存在性，令系統(tǒng)（3）的右端為零

xr11-xz-r2y1+αx-h=0，（4）

yr21-yz+βx1+αx=0，（5）

z（c-δx-σy）=0.（6）

系統(tǒng)（3）存在邊界平衡點E1=0，cσ，cσ和E2=cδ，0，r1cδ（r1-h）（hlt;r1時存在）。下面分析正平衡點的存在性，由方程（6）得

x=c-σyδ，（7）

將（7）代入方程（5），得

z=r2（δ+α（c-σy））yr2δ+（r2α+β）（c-σy），（8）

將（7）、（8）代入方程（4），得

Ay2+By+C=0，（9）

其中

A=r22δ2+r2ασ（r1-h）δ+r1σ2（r2α+β）

=-r2ασδ（h-h*），

B=-r2（r1-h）δ2-r2（r1σ+αc（r1-h））δ-

2cr1σ（r2α+β）

=（r2δ2+r2αcδ）h-h**，

C=r1c（r2（αc+δ）+βc），

h*=r22δ2+r1r2ασδ+r1σ2（r2α+β）r2ασδ，

h**=r1r2σδ+2cr1σ（r2α+β）+r1r2δ2+r1r2αcδr2δ2+r2αcδ.

由x，zgt;0得0lt;ylt;cσ，則方程（9）正根的存在性應(yīng)在區(qū)間0，cσ內(nèi)分析。

記：

Δ=B2-4AC，

顯然，hlt;h*，hlt;h**分別是Agt;0和Blt;0的充要條件。分析方程（9）正根的存在性，有以下結(jié)論：

（1）當hlt;min{h*，h**}時，若Δ≥0，方程（9）有2個正根；若Δlt;0，方程（9）無正根；

（2）當h**lt;hlt;h*時，方程（9）無正根；

（3）當hgt;h*時，方程（9）有唯一一個正根。

通過總結(jié)上面的分析，有以下定理：

定理1對于系統(tǒng)（3）

（1）系統(tǒng)始終存在無食餌種群的平衡點E1=0，cσ，cσ；

（2）當hlt;r1時， 系統(tǒng)存在無捕食者種群的平衡點E2=cδ，0，r1cδ（r1-h）；

（3）當hlt;min{h*，h**}時，若Δ≥0，系統(tǒng)存在2個正平衡點；當hgt;h*時，系統(tǒng)存在唯一一個正平衡點。則系統(tǒng)的正平衡點記為E*i=c-σy*iδ，y*i，r2（δ+α（c-σy*i））y*ir2δ+（r2α+β）（c-σy*i）。

3平衡點的穩(wěn)定性

系統(tǒng)（3）在平衡點E處的雅克比矩陣為

J|E=

r1-h-2r1xz-r2y（1+αx）2-r2x1+αxr1x2z2

βy（1+αx）2r2-2r2yz+βx1+αxr2y2z2

-δz-σzc-δx-σy.

下面分別討論每個平衡點的穩(wěn)定性。

對于平衡點E1，系統(tǒng)（3）在E1=0，cσ，cσ處的Jacobian矩陣為

JE1=r1-h-r2cσ00

βcσ-r2r2

-δcσ-c0.

顯然JE1存在1個特征值λ1=r1-h-r2cσ，其余的2個特征值λ2，λ3滿足方程：

λ2+r2λ+r2c=0.

容易判別Re（λ2），Re（λ3）lt;0。則平衡點E1的穩(wěn)定性由λ1決定。λ1lt;0，即hgt;r1-r2cσ，E1是穩(wěn)定的；λ1gt;0時，即hlt;r1-r2cσ，E1是不穩(wěn)定鞍點；λ1=0時，即h=r1-r2cσ，有1個零特征值，則需進一步分析。做線性變化（x^，y^，z^）=x，y-cσ，z-cσ把E1平移到原點。系統(tǒng)化為：

dx^dt=r1x^1-x^z^+cσ-r2x^y^+cσ1+αx^-hx^，

dy^dt=r2y^+cσ1-y^+cσz^+cσ+βx^y^+cσ1+αx^，

dz^dt=z^+cσc-δx^-σy^+cσ.（10）

將1z^+cσ和11+αx^在原點處泰勒展開，系統(tǒng)（10）可變?yōu)椋?/p>

dx^dt=a200x^2+a210x^2y^+a201x^2z^+…，

dy^dt=βcσx^-r2y^+r2z^+b200x^2+b110x^y^+b011y^z^+b020y^2+b210x^2y^+b021y^2z^+…，

dz^dt=-δcσx^-cy^-δx^z^-σy^z^，（11）

其中

a200=r2cασ+r1σc，a210=r2α，a201=r1σ2c2，

b200=-αβcσ，b110=β，b011=2r2σc，

b020=-r2σc，b210=-αβ，b021=r2σ2c2.

進一步做如下變換使系統(tǒng)（11）對角化

x^

y^

z^=100

p1p2p3

p4-11x～

y～

z～，

其中

p1=-δσ，p2=-r2+r22-4r2c2c，

p3=r2+r22-4r2c2c，p4=-βc-r2δr2σ.

對角化之后的系統(tǒng)（11）為：

dx～dt=a200x～2+S1（x～，y～，z～），

dy～dt=-r2（2+p2）y～+S2（x～，y～，z～），

dz～dt=r2（2-p3）z～+S3（x～，y～，z～），

其中Sl（x～，y～，z～）是（x～，y～，z～）中的冪級數(shù)，S1（x～，y～，z～）表示冪大于等于3的所有項，S2（x～，y～，z～）、S3（x～，y～，z～）表示其冪大于等于2的所有項。由于a200=r2cασ+r1σcgt;0，因此E1為鞍結(jié)點。

綜上所述，有以下結(jié)論：

定理2對于系統(tǒng)（3），無食餌種群的平衡點E1的穩(wěn)定性如下：

（1） 當hgt;r1-r2cσ時，E1是穩(wěn)定的；

（2） 當hlt;r1-r2cσ時，E1是不穩(wěn)定鞍點；

（3）當h=r1-r2cσ時，E1是鞍結(jié)點。

對于平衡點E2，系統(tǒng)（3）在E2=cδ，0，r1cδ（r1-h）處的Jacobian矩陣為

JE2=-r1+h-r2cδ+αc（r1-h）2r1

0r2+βcδ+αc0

-r1cr1-h-σr1cδ（r1-h）0.

顯然JE2存在1個特征值λ1=r2+βcδ+αcgt;0，其余的2個特征值λ2，λ3滿足方程：

λ2+（r1-h）λ+4c（r1-h）=0.

由平衡點E2的存在性可知hlt;r1，可得Re（λ2），Re（λ3）lt;0，故E2為鞍點。

對于平衡點E*i，系統(tǒng)在E*i=c-σy*iδ，y*i，r2（δ+α（c-σy*i））y*ir2δ+（r2α+β）（c-σy*i）處的Jacobian矩陣為

JE*=-r1x*z*+r2αx*y*（1+αx*）2-r2x*1+αx*r1x*2z*2

βy*1+αx*-αβx*y*（1+αx*）2-r2y*z*r2y*2z*2

-δz*-σz*0.

其特征方程為

λ3+a2λ2+a1λ+a0=0，

其中

a2=r1x*z*-r2αx*y*（1+αx*）2+r2y*z*，

a1=r1x*z*-r2αx*y*（1+αx*）2r2y*z*+r1δx*2z*+

r2σy*2z*+r2βx*y*（1+αx*）3，

a0=r2x*y*z*2r1x*z*-r2y*（1+αx*）+

r1βx*2y*（1+αx*）2z*2+

r2y*2z*2r1x*z*-r2αx*y*（1+αx*）2.

根據(jù)Routh-Hurwitz判據(jù)，當a2，a0gt;0且a1a2-a0gt;0時，其所有特征根具有負實部，即E*i是穩(wěn)定的。

定理3對于系統(tǒng)（3），無捕食者種群的平衡點E2和正平衡點E*i的穩(wěn)定性如下：

（1） 無捕食者種群的平衡點E2始終為鞍點；

（2） 當a2，a0gt;0且a1a2-a0gt;0時，正平衡點E*i是穩(wěn)定的。

4數(shù)值模擬

本節(jié)通過Matlab模擬系統(tǒng)（3）的動力學性態(tài)。為了體現(xiàn)收獲率h對模型的影響，取不同的h分別作圖。

圖1中，當h=0.9時，hgt;r1-r2cσ，系統(tǒng)趨于無食餌種群的平衡點E1（圖1A）。當h=0.6時，hlt;r1-r2cσ，無食餌種群的平衡點E1不穩(wěn)定，系統(tǒng)趨于正平衡點E*（圖1B），驗證了定理2。此時hlt;h*=1.303，hlt;h**=1.356。取3組不同的初值做了關(guān)于y的時間序列圖。顯然，當h=0.9時，系統(tǒng)最終趨于平衡點E1穩(wěn)定（圖2A）。當h=0.6時，系統(tǒng)最終趨于正平衡點E*穩(wěn)定（圖2B），則系統(tǒng)（3）存在1個穩(wěn)定的正平衡點。

接下來模擬了平衡點的三維空間圖。取3個不同的初始點（x，y，z）=0.5，0.2，0.3， （x，y，z）=（0.2，0.2，0.5），（x，y，z）=（0.3，0.1，0.6）。

當h=0.9時，系統(tǒng)趨于平衡點E1，當h=0.6時，系統(tǒng)趨于平衡點E*，如圖3所示。在圖4中，模擬了隨著收獲率h的變化，正平衡點E*的時間序列圖。圖像表明，通過減少收獲率h，食餌種群的數(shù)量會逐漸增大。因此，要想使食餌種群保持相對高的數(shù)量，應(yīng)當適當降低收獲率h。

5總結(jié)

本文考慮在實際的生態(tài)系統(tǒng)中，捕食者和食餌常常會共享資源。為了更好地理解捕食者和食餌之間共享資源的相互作用以及收獲的影響，建立了具有線性收獲的Holling II型功能反應(yīng)函數(shù)的捕食者和食餌共享資源的三維生態(tài)模型。對系統(tǒng)（2）進行無量綱化后，分析所得系統(tǒng)（3）其平衡點的存在性和穩(wěn)定性。然后取合適的參數(shù)模擬正平衡點E*的全局漸近穩(wěn)定性。最后，模擬收獲率h對系統(tǒng)（3）的影響，發(fā)現(xiàn)降低收獲率h到一定條件下，可使種群在更高數(shù)量下共存。因此，可以得到在保證生態(tài)穩(wěn)定的前提下人類捕撈的最優(yōu)收獲，即在帶有收獲項的捕食-食餌系統(tǒng)中，如果想讓食餌種群保持在一個穩(wěn)定的存活狀態(tài)，收獲率h不能過高，要滿足hlt;r1-r2cσ，否則，食餌種群就會滅絕。

參考文獻（References）

［1］SAFUAN H M， SIDHU H S， JOVANOSKI Z， et al. Impacts of a biotic resource enrichment on a predator-prey population［J］. Bulletin of Mathematical Biology， 2013， 75（10）： 1798-1812.

［2］CAPONE F， CARFORA M F， et al. On the dynamics of an intraguild predator-prey model［J］. Mathematics amp; Computers in Simulation， 2018， 149： 17-31.

［3］ANG T K， SAFUAN H N. Dynamical behaviors and optimal harvesting of an intraguild prey-predator fishery model with Michaelis-Menten type predator harvesting［J］. Biosystems， 2021， 202（2）： 1-10.

［4］GANGULI C， KAR T K， MONDAL P K. Optimal harvesting of a prey-predator model with variable carrying capacity［J］. International Journal of Biomathematics， 2017， 10（5）： 185-208.

［5］SHANG Z C， QIAO Y H， DUAN L J， et al. Stability and bifurcation analysis in a nonlinear harvested predator-prey model with simplified Holling type IV functional response［J］. International Journal of Bifurcation and Chaos in Applied Sciences and Engineering， 2020， 30（14）： 1-23.

［6］XIAO Z， LI Z， ZHU Z， et al. Hopf bifurcation and stability in a Beddington-DeAngelis predator-prey model with stage structure for predator and time delay incorporating prey refuge［J］. Open Mathematics， 2019， 17（1）： 141-159.

［7］DAWED M Y， KEBEDOW K G. Coexistence and harvesting optimal policy in three species food chain model with general Holling type functional response［J］. Natural Resource Modeling， 2021， 34（3）： 1-25.

［8］吳耀沖. 具有收獲率的時滯捕食者-食餌擴散系統(tǒng)的分支分析［D］. 廣東： 廣東工業(yè)大學， 2023.

［9］王菲， 呂堂紅， 周林華. 具時滯和線性收獲項的Leslie-Gower捕食-食餌系統(tǒng)Hopf分支［J］. 哈爾濱師范大學自然科學學報， 2022， 38（5）： 1-9.

WANG F， LV T H， ZHOU L H. Hopf bifurcation of Leslie-Gower predator-prey system with time delay and linear harvesting term［J］. Journal of Natural Sciences of Harbin Normal University， 2022， 38（5）： 1-9.

［10］CHAKRABORTY K， JANA S， KAR T K. Global dynamics and bifurcation in a stage structured prey-predator fishery model with harvesting［J］. Applied Mathematics and Computation， 2012， 218（18）： 9271-9290.

［11］XIAO M， CAO J. Hopf bifurcation and non-hyperbolic equilibrium in a ratio-dependent predator-prey model with linear harvesting rate： Analysis and computation［J］. Mathematical and Computer Modelling， 2009， 50（3/4）： 360-379.

［12］SAFUAN H M， MUSA S B. Food chain model with competition interaction in an environment of a biotic resource［C］// Advances in Industrial and Applied Mathematics. Proceedings of 23rd Malaysian National Symposium of Mathematical Sciences （SKSM23）. Malaysia， 2016， 1750（1）： 030015.

［13］ANG T K， SAFUAN H M. Harvesting in a toxicated intraguild predator-prey fishery model with variable carrying capacity［J］. Chaos， Solitons amp; Fractals， 2019， 126： 158-168.

［14］ABDULLAH S N， SAFUAN H M， NOR M E， et al. Harvesting effects on population model of competitive interaction with shared biotic resource［C］// Mathematical Sciences as the Core of Intellectual Excellence. Proceedings of the 25th National Symposium on Mathematical Sciences （SKSM25）. Malaysia， 2018， 1974： 020078.

［15］張志揚， 李忠. 一類具有線性捕獲和其他食物來源的捕食-食餌模型的穩(wěn)定性研究［J］. 延邊大學學報（自然科學版）， 2023， 49（1）： 21-25.

ZHANG Z Y， LI Z. Stability study of a predator-prey model with linear capture and other food sources［J］. Journal of Yanbian University （Natural Science Edition）， 2023， 49（1）： 21-25.

［16］張欣欣， 劉萍. 具有非線性捕獲量的捕食-食餌模型的動力學分析［J］. 哈爾濱師范大學自然科學學報， 2017， 33（2）： 9-12.

ZHANG X X， LIU P. Dynamic analysis of a predator-prey model with nonlinear capture quantity［J］. Journal of Natural Sciences of Harbin Normal University， 2017， 33（2）： 9-12.

［17］馮倩. 一類具有庇護所與非線性收獲項的捕食系統(tǒng)的穩(wěn)定性分析［D］. 甘肅： 蘭州交通大學， 2022.

［18］王甜. 具有收獲項的捕食-食餌模型的分支分析［D］. 安徽： 安徽大學， 2022.

［19］張蓬霞. 一類食餌具有收獲率的帶功能反應(yīng)的捕食模型定性分析［J］. 太原師范學院學報（自然科學版）， 2021，20（2）： 33-35.

ZHANG P X. Qualitative analysis of a predatory model with functional response and harvesting rate for a type of prey［J］. Journal of Taiyuan Normal University （Natural Science Edition）， 2021，20（2）： 33-35.

［20］TIAN Y， GAO Y， SUN K. Qualitative analysis of exponential power rate fishery model and complex dynamics guided by a discontinuous weighted fishing strategy［J］. Communications in Nonlinear Science and Numerical Simulation， 2023， 118（4）： 1-15.

（責任編輯：編輯郭蕓婕）

收稿日期：中文收稿日期2023-11-09

基金項目：山西省自然科學基金項目（201901D111179）

作者簡介：張雪婷（2000—），女，碩士研究生，專業(yè)方向為應(yīng)用數(shù)學。

*通信作者：李桂花（1973—），女，教授，從事應(yīng)用數(shù)學方向的研究，e-mail：liguihua@nuc.edu.cn。

DOI：10.13880/j.cnki.65-1174/n.2024.23.019

文章編號：1007-7383（2024）03-0383-07