一类具有收获和恐惧效应的捕食生态模型动力学问题研究

Dynamical Analysis of a Predator-Prey Ecological Model with Harvesting and Fear Effect

摘要: 本文提出了一类具有线性收获和恐惧效应的Leslie-Gower捕食生态模型，推导出模型所有可能平衡点存在与稳定的临界条件，给出模型发生跨临界分岔、鞍结分岔、Hopf分岔的阈值条件，模拟出模型特定分岔动力学性态的演化进程，从分岔动力学的角度揭示种群生长共存模式及其内在驱动机制。希望这些结果有利于拓宽捕食生态模型复杂动力学问题研究框架。

Abstract: This paper proposed a Leslie-Gower predator-prey ecological model with linear harvesting and fear effect. The critical conditions for the existence and stability of all possible equilibrium points are derived. The threshold conditions for the occurrence of transcritical bifurcation, saddle-node bifurcation, and Hopf bifurcation are explored. The evolutionary processes of specific bifurcation dynamical behaviors are simulated. The population growth coexistence modes and their underlying driving mechanisms could be revealed from the perspective of bifurcation dynamics. In a word, it is my hope that these results will help broaden the research framework for complex dynamic problems of predator-prey ecological models.

1. 引言

随着环境问题加剧和人类活动对自然生态的破坏加剧，水体富营养化已成为全球性的水环境问题之一[1] [2]。在经济、高效评估与应对自然环境挑战的需求驱动下，理论生态学家与数学家通过构建生态数学模型，对生态模型的复杂作用机制展开量化解析。其中，捕食者–食饵互作关系因其对食物网稳定和能量流动的关键调控作用，始终是生态模型研究的核心命题。此类研究不仅揭示物种协同演化规律，还为生物入侵防控、生态恢复工程等实践提供一定理论支撑[3]。

在1920年代，Lotka [4]与Volterra [5]共同构建了一类捕食者–食饵模型理论体系。自此之后，众多研究者持续拓展该模型的生态学解释维度，后续研究通过引入多维度生态因子(包括环境污染暴露、捕捞强度调控、食饵种群庇护效应、种群恐惧响应、食物资源补给及疾病传播机制) [6]-[10]，逐步构建起解释种间互作的复杂模型分析框架，相关成果为现代生态捕食动力学研究奠定了重要理论基础。

同时，针对Lotka-Volterra捕食者–食饵模型，许多数学家进一步扩展了捕食者–食饵模型，得到了几类新的Leslie-Gower捕食者–食饵模型[11] [12]。同时，一类新的Leslie-Gower捕食者–食饵模型可以在文献[13]中发现，主要创新是食饵和捕食者都具有各自的生长速度上限

$\{\begin{array}{l}\dot{x}=rx\left(1-\frac{x}{K}\right)-f\left(x\right)y,\hfill \\ \dot{y}=sy\left(1-\frac{y}{hx}\right),\hfill \end{array}$ (1)

其中$x$ 和$y$ 分别代表食饵和捕食者的种群密度；$r$ 和$s$ 分别表示食饵和捕食者种群的内禀增长率；$hx$ 是捕食者的环境容纳量，其与食饵种群密度有关；$K$ 表示是环境对食饵的最大容纳量，并且它与食饵种群数量成正比；$f\left(x\right)$ 被称为Leslie-Gower项。显然，不同的功能反应函数$f\left(x\right)$ 和生态因素产生不同的影响，这可能导致多样的动力学结果。因此，许多研究学者通过Leslie-Gower捕食者–食饵模型引入了不同的功能反应函数和生态因素。

针对此模型中的捕食者对食饵密度的线性功能反应函数，1965年Holling [14]在大量实验和分析的基础上，提出了三类不同类型的功能反应函数，其中的一类为Holling II型功能反应函数，即$f\left(x\right)=\frac{qx}{x+a}$ ，在此基础上改进的捕食生态模型如下

$\{\begin{array}{l}\dot{x}=rx\left(1-\frac{x}{K}\right)-\frac{qxy}{x+a},\hfill \\ \dot{y}=sy\left(1-\frac{y}{hx}\right).\hfill \end{array}$ (2)

Allee在1931年观察了种群的波动，发现如果种群密度过低，出生率会下降，而死亡率会上升，这类现象被称为Allee效应。Allee效应对逻辑增长的影响可用$x-M$ 的形式来表示，其中$M$ 是Allee效应阈值。同时，Allee效应主要分为弱Allee效应和强Allee效应[15]-[17]，强Allee效应是指种群密度低于一定阈值时，种群增长率为负；弱Allee效应是指种群密度低于一定的阈值时，种群增长率为正。据文献[18]可知，对于强Allee效应$0<M<1$ ，对于弱Allee效应$-1<M<0$ 。

如果假设模型(2)中的食饵种群可以通过群体效应或者集群效应来提升存活几率，这种反捕食行为就是存活上的Allee效应，因此我们引入Allee效应，模型(2)可转化为如下模型

$\{\begin{array}{l}\dot{x}=rx\left(1-\frac{x}{K}\right)\left(x-M\right)-\frac{bxy}{c+x},\hfill \\ \dot{y}=sy\left(1-\frac{y}{hx}\right).\hfill \end{array}$ (3)

许多实验[19]-[26]表明，捕食者并不一定仅仅通过捕猎来影响食饵种群。有时，捕食者种群的存在会影响食饵种群行为和心理，从而改变食饵种群的数量。然而，许多学者只考虑了捕食者的直接捕食，忽视了捕食者本身存在会对食饵密度产生影响的事实。因此，有必要讨论捕食者对食饵带来的恐惧影响。Wang等人[25]提出了一类具有恐惧效应的捕食–食饵生态模型，其中恐惧函数可表示为$f\left(x\right)=\frac{1}{1+\alpha y}$ ，$\text{α}$ 表示食饵的恐惧程度，研究发现恐惧程度越高，生态模型就越稳定，但是当恐惧水平降低时，生态模型会产生多个极限环。此后，许多学者受到这一思想的影响，开始研究具有恐惧效应的捕食–食饵生态模型。Cong等人在论文[24]中发现恐惧效应会使模型从混沌状态转变为稳定状态。Wang等人[26]研究了一类改进的Leslie-Gower生态模型的特性，表明恐惧效应丰富了模型的动力学性态。

如果我们假设模型(3)中的食饵受到捕食者种群捕食时，食饵种群会产生一定的恐惧效应，进而降低食饵种群的繁殖能力，这就是恐惧生态学效应。同时，为了维持捕食者种群与食饵种群的持久生存，需要对捕食者种群进行一定的收获，这样不仅可以降低捕食者种群的数量，也可以为保护措施的实施提供一定的资金支持，有利于维护生态系统的持久生态平衡。因此，针对模型(2)，我们同时考虑食饵种群的Allee和恐惧效应，以及对捕食者种群进行线性收获，这既可以比较真实直观刻画捕食者和食饵的种群作用生态关系，又可以维持保护机制的持久循环和种群持久生长共存模式的稳定。因此，我们提出新的模型如下

$\{\begin{array}{l}\dot{x}=\frac{rx}{1+ay}\left(1-\frac{x}{K}\right)\left(x-M\right)-\frac{bxy}{c+x},\hfill \\ \dot{y}=sy\left(1-\frac{y}{hx}\right)-qmEy,\hfill \end{array}$ (4)

其中$q$ 为捕食者的捕获系数；$m\left(0<m<1\right)$ 为可供捕获的种群比例；$E$ 为捕捞强度；$r$ 和$s$ 分别表示食饵和捕食者种群的内禀增长率；$K$ 表示环境对食饵的最大容纳量；$hx$ 是捕食者的环境容纳量，其与食饵种群密度有关；$\alpha$ 表示食饵的恐惧程度；$c$ 为半饱和系数；$M\left(0<M<K\right)$ 表示Allee效应阈值。

经过无量纲化

$\overline{x}=\frac{x}{k},\overline{y}=\frac{b}{rkc}y,\overline{t}=rkt.$

原模型(4)变为

$\{\begin{array}{l}\dot{x}=\frac{x}{1+ky}\left(1-x\right)\left(x-m\right)-\frac{xy}{1+\alpha x},\hfill \\ \dot{y}=\delta y\left(\beta -\frac{y}{x}\right)-\lambda y,\hfill \end{array}$ (5)

其中$k=\frac{rKca}{b},m=\frac{M}{K},\alpha =\frac{k}{c},\delta =\frac{sc}{Kbh},\beta =\frac{bh}{rc},\lambda =\frac{qmE}{rK}$ 都为正数。

2. 模型的有界性

定理1：若初值$\left(x\left(0\right),y\left(0\right)\right)\in \Omega =\left\{\left(x,y\right)\in {ℝ}^2|x>0,y>0\right\}$ ，则模型(5)的解总是正的。

证明。模型(5)简化为：

$\{\begin{array}{l}\dot{x}=xF\left(x,y\right)\hfill \\ \dot{y}=yG\left(x,y\right)\hfill \end{array}$

其中

$F\left(x,y\right)=\frac{\left(1-x\right)\left(x-m\right)}{1+ky}-\frac{y}{1+\alpha x},G\left(x,y\right)=\delta \left(\beta -\frac{y}{x}\right)-\lambda$

由于$\dot{x}=xF\left(x,y\right)$ 等价于$\frac{\text{d}x}{x}=F\left(x,y\right)\text{d}t$ ，得到：

$x\left(t\right)=x\left(0\right)\exp {\displaystyle {\int }_0^t\left\{\frac{1}{1+ky\left(\xi \right)}\left(1-x\left(\xi \right)\right)\left(x\left(\xi \right)-m\right)-\frac{y\left(\xi \right)}{1+\alpha x\left(\xi \right)}\right\}\text{d}\xi }$

由于$\dot{y}=yG\left(x,y\right)$ 等价于$\frac{\text{d}y}{y}=F\left(x,y\right)\text{d}t$ ，得到

$y\left(t\right)=y\left(0\right)\exp {\displaystyle {\int }_0^t\left\{\delta \left(\beta -\frac{y\left(\xi \right)}{x\left(\xi \right)}\right)-\lambda \right\}\text{d}\xi }$

这证明$\left(x\left(0\right),y\left(0\right)\right)\in \Omega =\left\{\left(x,y\right)\in {ℝ}^2|x>0,y>0\right\}$ ，模型(5)的解总是正的。

定理2：若参数满足$\lambda <\delta \beta$ ，则模型(5)解可以被限制在一个正集合内： $\overline{\Omega }=\left\{\left(x\left(t\right),y\left(t\right)\right)\in {ℝ}_{+}^2|0<x\left(t\right)<1,0<y\left(t\right)<\beta -\frac{\lambda }{\delta }\right\}$ 。

证明：对模型(5)的第一个方程进行分析：

$\dot{x}<x\left(1-x\right)\left(x-m\right)$ .

分三类情况讨论解的有界性：

1) 当$x\left(t\right)>1$ 时，则动力学方程满足$\dot{x}<0$ ；

2) 当$0<x\left(t\right)<m$ 时，则动力学方程满足$\dot{x}<0$ ；

3) 当$m<x\left(t\right)<1$ 时，通过分析可得：$x\left(t\right)<1-\frac{1}{1+C_1{e}^{\left(1-m\right)t}}$ ，进而证明上极限满足：$\underset{t\to +\infty }{\lim \sup }y\left(t\right)\le 1$ 。

(1)和(2)两种情况可以通过对原方程进行变换分析进一步论证，构造上界估计：

$\dot{x}<x\left(1-x\right)\left(1-m\right)$

变量分离后得到可积分形式：

$\left(\frac{1}{x}+\frac{1}{1-x}\right)\text{d}x<\left(1-m\right)\text{d}t$

积分运算导出解析关系：

$\frac{x}{1-x}<C_1{\text{e}}^{\left(1-m\right)t}$

整理后得到关键约束：

$x\left(t\right)<1-\frac{1}{1+C_1{\text{e}}^{\left(1-m\right)t}}$

由此可证：在$\left(x\left(t\right),y\left(t\right)\right)\in {ℝ}_{+}^2$ 定义域内，模型始终满足$0<x\left(t\right)<1$ 。

进而对第二个方程进行捕食者种群动力学分析，将通过两种情况来确定状态变量$y\left(t\right)$ 的渐近性态：

1) 当$y\ge \beta -\frac{\lambda }{\delta }$ 时，存在$\dot{y}<0$ ，与定理1正解条件矛盾

2) 当$0<y<\beta -\frac{\lambda }{\delta }$ 时，通过微分不等式得：

$y\left(t\right)<\beta -\frac{\lambda }{\delta }-\frac{\beta -\frac{\lambda }{\delta }}{1+C_2{\text{e}}^{\delta \left(\beta -\frac{\lambda }{\delta }\right)t}}$

通过上极限分析可得：

$\underset{t\to +\infty }{\lim \sup }y\left(t\right)\le \beta -\frac{\lambda }{\delta }.$

(1)的非存在性可通过反证法直接证得。现论证(2)的动力学约束机制，对动力学方程进行变量分离处理：

$\left(\frac{1}{y}+\frac{1}{\beta -\frac{\lambda }{\delta }-y}\right)\text{d}y<\delta \left(\beta -\frac{\lambda }{\delta }\right)\text{d}t,$

通过双侧积分得到关键关系式：

$\frac{y}{\beta -\frac{\lambda }{\delta }-y}<C_2{\text{e}}^{\delta \left(\beta -\frac{\lambda }{\delta }\right)t}.$

整理后得到显式上界估计：

$y<\beta -\frac{\lambda }{\delta }-\frac{\beta -\frac{\lambda }{\delta }}{1+C_2{\text{e}}^{\delta \left(\beta -\frac{\lambda }{\delta }\right)t}}.$

渐近分析表明当$t\to +\infty$ 时，最终得：

$\underset{t\to +\infty }{\lim \sup }y\left(t\right)\le \beta -\frac{\lambda }{\delta }.$

结合对$x\left(t\right)$ 的约束条件，与当前对$y\left(t\right)$ 的估计，可得不变集：

$\Omega =\left\{\left(x\left(t\right),y\left(t\right)\right)\in {ℝ}_{+}^2|0<x\left(t\right)<1,0<y\left(t\right)<\beta -\frac{\lambda }{\delta }\right\}.$

证毕。

3. 平衡点的存在性

显然，边界平衡点$E_1\left(1,0\right)$ 和$E_2\left(m,0\right)$ 总是存在的。模型(5)的正平衡点满足以下方程：

$\{\begin{array}{l}\frac{\left(1-x\right)\left(x-m\right)}{1+ky}=\frac{y}{1+\alpha x},\hfill \\ \delta \left(\beta -\frac{y}{x}\right)=\lambda .\hfill \end{array}$ (6)

从模型(6)的第二个方程中，我们得到$y=\left(\beta -\lambda \delta \right)x$ 。为了满足正平衡点的存在性，要求$\beta -\frac{\lambda }{\delta }>0$ 。将$y=\left(\beta -\frac{\lambda }{\delta }\right)x$ 代入方程(6)的第一个方程中，我们得到：

$\alpha x^3+\left[1+{\left(\beta -\frac{\lambda }{\delta }\right)}^{2}k-\alpha m-\alpha \right]x^2+\left(\beta -\frac{\lambda }{\delta }+\alpha m-m-1\right)x+m=0$ .

令

$f\left(x\right)=\alpha x^3+\left[1+{\left(\beta -\frac{\lambda }{\delta }\right)}^{2}k-\alpha m-\alpha \right]x^2+\left(\beta -\frac{\lambda }{\delta }+\alpha m-m-1\right)x+m$ .

显然$f\left(0\right)=m>0$ ，因为$\alpha >0$ ，所以$\underset{x\to -\infty }{\lim }f\left(x\right)=-\infty$ ，则总是存在一个零点属于$\left(-\infty ,0\right)$ 。

对$f\left(x\right)$ 求导得到：

$h\left(x\right)=f^{\prime }\left(x\right)=3\alpha x^2+2\left[1+{\left(\beta -\frac{\lambda }{\delta }\right)}^{2}k-\alpha m-\alpha \right]x+\beta -\frac{\lambda }{\delta }+\alpha m-m-1$ .

显然，$3\alpha >0$ ，因此抛物线函数$h\left(x\right)$ 向上开口。

对于$h\left(x\right)$ ，我们令：

$A=\alpha ,B=2\left[1+{\left(\beta -\frac{\lambda }{\delta }\right)}^{2}k-\alpha m-\alpha \right],C=\beta -\frac{\lambda }{\delta }+\alpha m-m-1$ .

然后，对于$h\left(x\right)=0$ ，我们有判别式$\Delta =B^2-4AC$ 。

为了验证模型(5)具有两个正平衡点，对相关方程进行了数值模拟，参数取值分别为$k=0.88$ ，$m=0.3$ ，$\alpha =0.6$ ，$\delta =0.8$ ，$\beta =1$ ，$\lambda =0.66$ 。结果如图1所示，根据$f\left(x\right)$ 和$h\left(x\right)$ 的性质，我们可以确定在区间$\left[0,1\right]$ 内，模型(5)具有两个正平衡点。

定理3：对于正平衡点的个数，有：

• 当$\text{Δ}>0$ 时，$h\left(x\right)$ 有两个根，分别为$x_{h1}=\frac{-B-\sqrt{\Delta }}{2A}$ ，$x_{h2}=\frac{-B+\sqrt{\Delta }}{2A}$ 。

1) 当$x_{h2}\ge 0$ 时：

a) 如果$f\left(x_{h2}\right)>0$ ，则没有正平衡点。

Figure 1. Existence of equilibrium point

图1. 平衡点的存在性情况

b) 如果$f\left(x_{h2}\right)=0$ ，则有一个平衡点$E_{+}\left(x_{+},y_{+}\right)$ ，其中$x_{+}$ 是一个双正根。

c) 如果$f\left(x_{h2}\right)<0$ ，则存在两个正平衡点，$E_{+}^{*}\left(x_{+}^{*},y_{+}^{*}\right)$ 和$E_{+}^{**}\left(x_{+}^{**},y_{+}^{**}\right)$ ，其中$x_{+}^{*}<x_{+}^{**}$ 。$x_{+}^{*}$ 和$x_{+}^{**}$ 都表示单一的正根。

2) 当$x_{h2}<0$ 时，则没有正平衡点。

• 当$\text{Δ}\le 0$ 时，模型没有正平衡点。

证明：根据求根公式当$\text{Δ}>0$ 时，$h\left(x\right)$ 有两个根，分别为$x_{h1}=\frac{-B-\sqrt{\Delta }}{2A}$ ，$x_{h2}=\frac{-B+\sqrt{\Delta }}{2A}$ 。

• 当$\Delta >0$ 时，根据求根公式$h\left(x\right)$ 有两个根，分别为$x_{h1}=\frac{-B-\sqrt{\Delta }}{2A}$ ，$x_{h2}=\frac{-B+\sqrt{\Delta }}{2A}$ 。

1) 当$x_{h2}\ge 0$ 时：因为$f\left(0\right)=m>0$ ，$h\left(x\right)$ 是$f\left(x\right)$ 的导数且当$x\in \left(-\infty ,x_{h1}\right)$ 时，$h\left(x\right)\ge 0$ ；当$x\in \left(x_{h1},x_{h2}\right)$ ，时$h\left(x\right)\le 0$ 。当$x\in \left(x_{h2},+\infty \right)$ 时，$h\left(x\right)\ge 0$ 。所以$x\in \left(-\infty ,x_{h1}\right)$ 时，$h\left(x\right)$ 单调递增。当$x\in \left(x_{h1},x_{h2}\right)$ ，$h\left(x\right)$ 单调递减。当$x\in \left(x_{h2},+\infty \right)$ 时，$h\left(x\right)$ 单调递增。所以有：

a) 如果$f\left(x_{h2}\right)>0$ ，则没有正平衡点。

b) 如果$f\left(x_{h2}\right)=0$ ，则有一个平衡点$E_{+}\left(x_{+},y_{+}\right)$ ，其中$E_{+}$ 是一个双正根。

c) 如果$f\left(x_{h2}\right)<0$ ，则存在两个正平衡点，$E_{+}^{*}\left(x_{+}^{*},y_{+}^{*}\right)$ 和$E_{+}^{\text{**}}\left(x_{+}^{\text{**}},y_{+}^{\text{**}}\right)$ ，其中$x_{+}^{*}<x_{+}^{**}$ 。$x_{+}^{*}$ 和$x_{+}^{**}$ 都表示单一的正根。

2) 同理，当$x_{h2}<0$ 时，则没有正平衡点。

• 同理当$\text{Δ}\le 0$ 时，模型没有平衡点。

定理4：对于定理3中存在两个正平衡点的情况，总有$h\left(x_{+}^{*}\right)<0$ ，$h\left(x_{+}^{**}\right)>0$ 。针对定理4，我们利用数值模拟进行的验证，参数取值和图1是一样的，其结果如图2所示。

Figure 2. The $h\left(x\right)$ value of positive equilibrium point

图2. 正平衡点的$h\left(x\right)$ 值

4. 平衡点的稳定性

现在，我们将讨论平衡点的稳定性。

模型(5)的雅可比矩阵为：

$J=\left[\begin{array}{cc}\frac{\left(1-x\right)\left(x-m\right)}{1+ky}-\frac{x\left(x-m\right)}{1+ky}+\frac{x\left(1-x\right)}{1+ky}-\frac{y}{1+\alpha x}+\frac{\alpha xy}{{\left(1+\alpha x\right)}^2}& -\frac{x\left(1-x\right)\left(x-m\right)k}{{\left(1+ky\right)}^2}-\frac{x}{1+\alpha x}\\ \frac{\delta y^2}{x^2}& \delta \left(\beta -\frac{y}{x}\right)-\frac{\delta y}{x}-\lambda \end{array}\right].$

因此，我们得出以下定理。

4.1. 边界平衡点$E_1$ 的稳定性

定理5：边界平衡点$E_1\left(1,0\right)$ 的稳定性讨论如下：

1) 如果$\lambda <\beta \delta$ ，则$E_1$ 是一个鞍点。

2) 如果$\lambda >\beta \delta$ ，则$E_1$ 是一个稳定结点。

3) 如果$\lambda =\beta \delta$ ，则$E_1$ 是一个吸引鞍结点。也就是说，$E_1$ 的一个适当小的邻域被两个分界线沿着$E_1$ 的上方和下方分成两部分，这两部分趋向于$E_1$ 。左半平面是一个抛物面区域，而右半平面是两个双曲区域。

证明：$E_1$ 的雅可比矩阵为：

$J_{E_1}=\left(\begin{array}{cc}m-1& -\frac{1}{\alpha +1}\\ 0& \beta -\frac{\lambda }{\delta }\end{array}\right)$ .

显然，$J_{E_1}$ 有两个特征值${\mu }_1=m-1$ 和${\mu }_2=\beta -\frac{\lambda }{\delta }$ 。因此：

1) 如果${\mu }_2=\beta -\frac{\lambda }{\delta }>0$ ，即，$\lambda <\beta \delta$ ，则$E_1$ 是一个双曲鞍点。

2) 如果${\mu }_2=\beta -\frac{\lambda }{\delta }<0$ ，即，$\lambda >\beta \delta$ ，则$E_1$ 是一个双曲稳定结点。

3) 如果${\mu }_2=\beta -\frac{\lambda }{\delta }=0$ ，即，$\lambda =\beta \delta$ ，则两个特征值为${\mu }_1=m-1$ 和${\mu }_2=0$ 。

为了研究当$\lambda =\beta \delta$ 时$E_1$ 的稳定性，我们通过$\left(X,Y\right)=\left(x-1,y\right)$ 将$E_1$ 转换到原点。同时，我们将模型(5)在原点附近展开到二阶，然后模型(5)变为：

$\{\begin{array}{l}\dot{X}=\left(m-1\right)X-\frac{1}{\alpha +1}Y+\left(m-2\right)X^2+a_{11}XY+P_1\left(X,Y\right),\hfill \\ \dot{Y}=-\delta Y^2+\delta X{Y}^2+Q_1\left(X,Y\right),\hfill \end{array}$ (7)

其中

$a_{11}=-\frac{1+{\left(\alpha +1\right)}^2\left(m-1\right)}{{\left(\alpha +1\right)}^2}$ .

$P_1\left(X,Y\right)$ 至少是$\left(X,Y\right)$ 的3阶函数，$Q_1\left(X,Y\right)$ 至少是$\left(X,Y\right)$ 的4阶函数。

通过变换：

$\left(\begin{array}{c}X\\ Y\end{array}\right)=\left(\begin{array}{cc}1& \frac{1}{1+\alpha }\\ 0& m-1\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)$ .

模型(7)变为标准形式：

$\{\begin{array}{l}\dot{x}=\left(m-1\right)x+\left(m-2\right)x^2+b_{11}xy+b_{02}{y}^2+P_2\left(x,y\right),\hfill \\ \dot{y}=-\delta \left(m-1\right)y^2+\delta \left(m-1\right)x{y}^2+\frac{\left(m-1\right)\delta }{\alpha +1}{y}^3+Q_2\left(x,y\right),\hfill \end{array}$ (8)

其中

$b_{11}=\frac{a_{11}\alpha m-a_{11}\alpha +m{a}_{11}-m+2m-4}{\alpha +1},$

$b_{02}=\frac{\alpha \delta m^2+a_{11}\alpha m-2\alpha m\delta +m^2\delta -a_{11}\alpha +m{a}_{11}+\alpha \delta -2m\delta -a_{11}+\delta +m-2}{{\left(\alpha +1\right)}^2},$

$P_2\left(x,y\right)$ 是至少是$\left(x,y\right)$ 的3阶函数，$Q_2\left(x,y\right)$ 至少是$\left(x,y\right)$ 的4阶函数。

引入一个新变量$\tau$ ，通过$\tau =\left(m-1\right)t$ ，我们得到：

$\{\begin{array}{l}\dot{x}=x+\frac{m-2}{m-1}{x}^2+\frac{b_{11}}{m-1}xy+\frac{b_{02}}{m-1}{y}^2+P_3\left(x,y\right),\hfill \\ \dot{y}=-\delta y^2+\delta x{y}^2+\frac{\delta }{\alpha +1}{y}^3+Q_3\left(x,y\right),\hfill \end{array}$ (9)

$P_3\left(x,y\right)$ 是至少是$\left(x,y\right)$ 的3阶函数，$Q_3\left(x,y\right)$ 是至少是$\left(x,y\right)$ 的4阶函数。

由于模型(9)的第二个方程中的$y^2$ 的系数为$-\delta <0$ ，所以平衡点$E_1\left(1,0\right)$ 是一个吸引鞍结点，这意味着$E_1$ 的一个适当小的区域被两个分界线沿着$E_1$ 的上方和下方分成两部分，这两部分都趋向于$E_1$ 。左半平面是一个抛物体扇形，右半平面是两个双曲扇形。

4.2. 平衡点$E_2$ 的稳定性

对平衡点$E_2$ 的研究与$E_1$ 的稳定性分析类似。模仿定理5的证明，我们可得到下面的定理。

定理6：边界平衡点$E_2\left(m,0\right)$ 的稳定性讨论如下：

1) 如果$\lambda <\beta \delta$ ，则$E_2$ 是一个结点。

2) 如果$\lambda >\beta \delta$ ，则$E_2$ 是一个鞍点。

3) 如果$\lambda =\beta \delta$ ，则$E_2$ 是一个排斥鞍结点。

4.3. 内部平衡点的稳定性

设$x^{*}$ 和$y^{*}$ 为内部平衡点的横纵坐标。则内部平衡点的雅可比矩阵为：

$J^{*}=\left[\begin{array}{cc}\frac{\left(1-x^{*}\right)\left(x^{*}-m\right)}{1+k{y}^{*}}-\frac{x^{*}\left(x^{*}-m\right)}{1+k{y}^{*}}+\frac{x^{*}\left(1-x^{*}\right)}{1+k{y}^{*}}-\frac{y^{*}}{1+\alpha x^{*}}+\frac{\alpha x^{*}{y}^{*}}{{\left(1+\alpha x^{*}\right)}^2}& -\frac{k{x}^{*}\left(1-x^{*}\right)\left(x^{*}-m\right)}{{\left(1+k{y}^{*}\right)}^2}-\frac{x^{*}}{1+\alpha x^{*}}\\ \frac{\delta y^{*2}}{x^{*2}}& \delta \left(\beta -\frac{y^{*}}{x^{*}}\right)-\frac{\delta y^{*}}{x^{*}}-\lambda \end{array}\right].$

在这一节中，我们通过简化雅可比矩阵来讨论内部平衡点的稳定性。显然，内部平衡点满足方程(6)。

因此，我们得到：

$\begin{array}{c}{J}_{11}^{*}=\frac{\left(1-x^{*}\right)\left(x^{*}-m\right)}{1+k{y}^{*}}-\frac{x^{*}\left(x^{*}-m\right)}{1+k{y}^{*}}+\frac{x^{*}\left(1-x^{*}\right)}{1+k{y}^{*}}-\frac{y^{*}}{1+\alpha x^{*}}+\frac{\alpha x^{*}{y}^{*}}{{\left(1+\alpha x^{*}\right)}^2}\\ =x^{*}\left[\frac{1-2x^{*}+m}{1+k{y}^{*}}+\frac{\alpha y^{*}}{{\left(1+\alpha x^{*}\right)}^2}\right]\\ =x^{*}\left[\frac{1-2x^{*}+m}{1+k{y}^{*}}+\frac{\alpha \left(1-x^{*}\right)\left(x^{*}-m\right)}{\left(1+\alpha x^{*}\right)\left(1+k{y}^{*}\right)}\right]\\ =x^{*}\frac{-3\alpha x^{*2}-2\left(1-\alpha m-\alpha \right)x^{*}-\left(\alpha m-m-1\right)}{\left(1+k{y}^{*}\right)\left(1+\alpha x^{*}\right)}\\ =-x^{*}\frac{3\alpha x^{*2}+2\left(1-\alpha m-\alpha \right)x^{*}+\left(\alpha m-m-1\right)}{\left(1+k{y}^{*}\right)\left(1+\alpha x^{*}\right)}\end{array}$

$\begin{array}{c}{J}_{12}^{*}=-\frac{k{x}^{*}\left(1-x^{*}\right)\left(x^{*}-m\right)}{{\left(k{y}^{*}+1\right)}^2}-\frac{x^{*}}{1+\alpha x^{*}}\\ =-x^{*}\left[\frac{k{y}^{*}}{\left(1+k{y}^{*}\right)\left(1+\alpha x^{*}\right)}+\frac{1}{1+\alpha x^{*}}\right]\\ =-x^{*}\frac{2k{y}^{*}+1}{\left(1+k{y}^{*}\right)\left(1+\alpha x^{*}\right)},\end{array}$

$J_{21}^{*}=\frac{\delta y^2}{x^2}=\delta {\left(\beta -\frac{\lambda }{\delta }\right)}^2,$

$J_{22}^{*}=\lambda -\delta \beta .$

然后，我们得到内部平衡点的雅可比矩阵：

$J^{*}=\left[\begin{array}{cc}-x^{*}\frac{3\alpha x^{*2}+2\left(1-\alpha m-\alpha \right)x^{*}+\left(\alpha m-m-1\right)}{\left(1+k{y}^{*}\right)\left(1+\alpha x^{*}\right)}& -x^{*}\frac{2k{y}^{*}+1}{\left(1+k{y}^{*}\right)\left(1+\alpha x^{*}\right)}\\ \delta {\left(\beta -\frac{\lambda }{\delta }\right)}^2& \lambda -\delta \beta \end{array}\right].$

因此，该矩阵的行列式和迹分别为：

$Det\left(J^{*}\right)=\frac{x^{*}\delta \left(\beta -\frac{\lambda }{\delta }\right)}{\left(1+k{y}^{*}\right)\left(1+\alpha x^{*}\right)}\left[\left(\beta -\frac{\lambda }{\delta }\right)\left(1+2k{y}^{*}\right)+3\alpha x^{*2}+2\left(1-\alpha m-\alpha \right)x^{*}+\left(\alpha m-m-1\right)\right],$

$Tr\left(J^{*}\right)=-x^{*}\frac{3\alpha x^{*2}+2\left(1-\alpha m-\alpha \right)x^{*}+\left(\alpha m-m-1\right)}{\left(1+k{y}^{*}\right)\left(1+\alpha x^{*}\right)}-\delta \left(\beta -\frac{\lambda }{\delta }\right).$

$Det\left(J^{\text{*}}\right)$ 可以简化为：

$\begin{array}{c}Det\left(J^{*}\right)=\frac{x^{*}\delta \left(\beta -\frac{\lambda }{\delta }\right)}{\left(1+k{y}^{*}\right)\left(1+\alpha x^{*}\right)}\left[\left(\beta -\frac{\lambda }{\delta }\right)\left(1+2k{y}^{*}\right)+3\alpha x^{*2}+2\left(1-\alpha m-\alpha \right)x^{*}+\left(\alpha m-m-1\right)\right]\\ =\frac{x^{*}\delta \left(\beta -\frac{\lambda }{\delta }\right)}{\left(1+k{y}^{*}\right)\left(1+\alpha x^{*}\right)}\left[\left(\beta -\frac{\lambda }{\delta }\right)+2\left(\beta -\frac{\lambda }{\delta }\right)k{y}^{*}+3\alpha x^{*2}+2\left(1-\alpha m-\alpha \right)x^{*}+\left(\alpha m-m-1\right)\right]\\ =\frac{x^{*}\delta \left(\beta -\frac{\lambda }{\delta }\right)}{\left(1+k{y}^{*}\right)\left(1+\alpha x^{*}\right)}\left[\left(\beta -\frac{\lambda }{\delta }\right)+2{\left(\beta -\frac{\lambda }{\delta }\right)}^{2}k+3\alpha x^{*2}+2\left(1-\alpha m-\alpha \right)x^{*}+\left(\alpha m-m-1\right)\right]\\ =\frac{x^{*}\delta \left(\beta -\frac{\lambda }{\delta }\right)}{\left(1+k{y}^{*}\right)\left(1+\alpha x^{*}\right)}\left[3\alpha x^{*2}+2\left(1+{\left(\beta -\frac{\lambda }{\delta }\right)}^{2}k-\alpha m-\alpha \right)x^{*}+\left(\alpha m-m-1\right)\right]\\ =\frac{x^{*}\delta \left(\beta -\frac{\lambda }{\delta }\right)}{\left(1+k{y}^{*}\right)\left(1+\alpha x^{*}\right)}h\left(x\right).\end{array}$ (10)

下面讨论平衡点$E_{+}$ 的稳定性

首先，我们通过平移$\left(X,Y\right)=\left(x-x_{+},y-y_{+}\right)$ 将平衡点$E_{+}$ 转换为原点，然后模型(5)变为

$\{\begin{array}{l}\dot{X}={a^{\prime }}_{10}X+{a^{\prime }}_{01}Y+{a^{\prime }}_{20}{X}^2+{a^{\prime }}_{11}XY+{a^{\prime }}_{02}{Y}^2+P_4\left(X,Y\right),\hfill \\ \dot{Y}={b^{\prime }}_{10}X+{b^{\prime }}_{01}Y+{b^{\prime }}_{20}{X}^2+{b^{\prime }}_{11}XY+{b^{\prime }}_{02}{Y}^2+Q_4\left(X,Y\right),\hfill \end{array}$ (11)

其中

$\begin{array}{l}{a^{\prime }}_{10}=-\frac{x_{+}\left[3\alpha x_{+}^2+2\left(1-\alpha m-\alpha \right)x_{+}+\left(\alpha m-m-1\right)\right]}{\left(k{y}_{+}+1\right)\left(\alpha y_{+}+1\right)},\\ {a^{\prime }}_{01}=-\frac{x_{+}\left(2k{y}_{+}+1\right)}{\left(k{y}_{+}+1\right)\left(\alpha x_{+}+1\right)},{a^{\prime }}_{20}=\frac{\left(m-3x_{+}+1\right)}{k{y}_{+}+1}-\frac{y_{+}\alpha }{{\left(\alpha x_{+}+1\right)}^3},\end{array}$

${a^{\prime }}_{11}=-\frac{\left(2x_{+}m-3x_{+}^2-m+2x_{+}\right)k}{{\left(k{y}_{+}+1\right)}^2}+\frac{1}{{\left(\alpha x_{+}+1\right)}^2},{a^{\prime }}_{02}=\frac{x_{+}{k}^2\left(x_{+}m-x_{+}^2-m+x_{+}\right)}{{\left(k{y}_{+}+1\right)}^3}$

${b^{\prime }}_{10}=\delta {\left(\beta -\frac{\lambda }{\delta }\right)}^2,{b^{\prime }}_{01}=\delta \left(\beta -\frac{\lambda }{\delta }\right),{b^{\prime }}_{20}=-\frac{\delta y_{+}^2}{x_{+}^3},{b^{\prime }}_{11}=\frac{2\delta y_{+}}{x_{+}^2},{b^{\prime }}_{02}=-\frac{\delta }{x_{+}}$

$P_4\left(x,y\right)$ 和$Q_4\left(x,y\right)$ 是至少为$\left(x,y\right)$ 的三阶无穷可微函数。

情况1：$\lambda \ne x_{+}\frac{3\alpha x_{+}^2+2\left(1-\alpha m-\alpha \right)x_{+}+\left(\alpha m-m-1\right)}{\left(1+k{y}_{+}\right)\left(1+\alpha x_{+}\right)}+\delta \beta$

在这种情况下，模型(11)的雅可比矩阵为：

$J_{E_{+}}=\left(\begin{array}{cc}{a^{\prime }}_{10}& {a^{\prime }}_{01}\\ \left(\beta -\frac{\lambda }{\delta }\right){b^{\prime }}_{01}& {b^{\prime }}_{01}\end{array}\right)$ .

由于$h\left(x_{+}\right)=0$ ，且${b^{\prime }}_{01}\ne 0$ ，因此

$\begin{array}{c}Det\left(J_{E_{+}}\right)={a^{\prime }}_{10}{b^{\prime }}_{01}-{a^{\prime }}_{01}\left(\beta -\frac{\lambda }{\delta }\right){b^{\prime }}_{01}\\ ={b^{\prime }}_{01}\left[{a^{\prime }}_{10}-{a^{\prime }}_{01}\left(\beta -\frac{\lambda }{\delta }\right)\right]\\ =\frac{x_{+}\delta \left(\beta -\frac{\lambda }{\delta }\right)}{\left(1+k{y}_{+}\right)\left(1+\alpha x_{+}\right)}h\left(x_{+}\right)=0.\end{array}$

$J_{E_{+}}$ 的特征值为${\mu }_1=0$ ，${\mu }_2={a^{\prime }}_{10}+{b^{\prime }}_{01}$ ，${a^{\prime }}_{01}=\frac{{a^{\prime }}_{10}}{\beta -\frac{\lambda }{\delta }}$ 。在这种情况下，模型(11)的雅可比矩阵为：

$J_{E_{+}}=\left(\begin{array}{cc}{a^{\prime }}_{10}& \frac{{a^{\prime }}_{10}}{\beta -\frac{\lambda }{\delta }}\\ \left(\beta -\frac{\lambda }{\delta }\right){b^{\prime }}_{01}& {b^{\prime }}_{01}\end{array}\right)$ .

接下来，我们通过以下变换将模型(11)简化为标准形式：

$\left(\begin{array}{l}X\hfill \\ Y\hfill \end{array}\right)=\left(\begin{array}{cc}{a^{\prime }}_{10}& -1\\ \delta {\left(\beta -\frac{\lambda }{\delta }\right)}^2& \beta -\frac{\lambda }{\delta }\end{array}\right)\left(\begin{array}{l}x\hfill \\ y\hfill \end{array}\right)$

模型(11)变为以下形式

$\{\begin{array}{l}\dot{x}={c^{\prime }}_{10}x+{c^{\prime }}_{20}{x}^2+{c^{\prime }}_{11}xy+{c^{\prime }}_{02}{y}^2+P_5\left(x,y\right),\hfill \\ \dot{y}={d^{\prime }}_{20}{x}^2+{d^{\prime }}_{11}xy+{d^{\prime }}_{02}{y}^2+Q_5\left(x,y\right),\hfill \end{array}$ (12)

其中

${c^{\prime }}_{10}=\delta p+{a^{\prime }}_{10},{c^{\prime }}_{20}=-\frac{{\delta }^2{p}^4{a^{\prime }}_{02}-{\delta }^2{p}^3{b^{\prime }}_{02}+{a^{\prime }}_{10}\delta p^2{a^{\prime }}_{11}+2\delta p^2{a^{\prime }}_{10}{b^{\prime }}_{02}-{a^{\prime }}_{10}^{2}p{b^{\prime }}_{02}+{a^{\prime }}_{10}^2{a^{\prime }}_{20}}{\delta p-{a^{\prime }}_{10}},$

${c^{\prime }}_{11}=-\frac{2p^3\delta {a^{\prime }}_{02}-\delta p^2{a^{\prime }}_{11}-4p^2\delta {b^{\prime }}_{02}+p{a^{\prime }}_{10}{a^{\prime }}_{11}+4{a^{\prime }}_{10}p{b^{\prime }}_{02}-2{a^{\prime }}_{10}{a^{\prime }}_{20}}{\delta p-{a^{\prime }}_{10}},{c^{\prime }}_{02}=-\frac{p^2{a^{\prime }}_{02}-p{a^{\prime }}_{11}-4{b^{\prime }}_{02}p+{a^{\prime }}_{20}}{\delta p-{a^{\prime }}_{10}},$

${d^{\prime }}_{20}=\frac{p\left({\delta }^3{p}^4{a^{\prime }}_{02}+{\delta }^2{p}^2{a^{\prime }}_{10}{a^{\prime }}_{11}-{\delta }^2{p}^2{a^{\prime }}_{10}{b^{\prime }}_{02}+2\delta p{a^{\prime }}_{10}^2{b^{\prime }}_{02}+\delta {a^{\prime }}_{10}^2{a^{\prime }}_{20}-{a^{\prime }}_{10}^3{b^{\prime }}_{02}\right)}{\delta p-{a^{\prime }}_{10}},$

${d^{\prime }}_{11}=\frac{p\left(2{\delta }^2{p}^3{a^{\prime }}_{02}-{\delta }^2{p}^2{a^{\prime }}_{11}+\delta p{a^{\prime }}_{10}{a^{\prime }}_{11}-4\delta p{a^{\prime }}_{10}{b^{\prime }}_{02}-2\delta {a^{\prime }}_{10}{a^{\prime }}_{20}+4{a^{\prime }}_{10}^2{b^{\prime }}_{02}\right)}{\delta p-{a^{\prime }}_{10}},$

${d^{\prime }}_{02}=\frac{p\left(\delta p^2{a^{\prime }}_{02}-\delta p{a^{\prime }}_{11}+\delta {a^{\prime }}_{20}-4{a^{\prime }}_{10}{b^{\prime }}_{02}\right)}{\delta p-{a^{\prime }}_{10}}.$

且$p=\beta -\frac{\lambda }{\delta }$ ，$P_5\left(x,y\right)$ 和$Q_5\left(x,y\right)$ 是至少为$\left(x,y\right)$ 的三阶无穷可微函数。

我们通过引入新变量$\tau$ ，使得$\tau =\left(\delta p+{a^{\prime }}_{10}\right)t$ ，为简便起见，我们仍然使用$t$ 而不是$\tau$ ，得到

$\{\begin{array}{l}\dot{x}=x+{e^{\prime }}_{20}{x}^2+{e^{\prime }}_{11}xy+{e^{\prime }}_{02}{y}^2+P_6\left(x,y\right),\hfill \\ \dot{y}={f^{\prime }}_{20}{x}^2+{f^{\prime }}_{11}xy+{f^{\prime }}_{02}{y}^2+Q_6\left(x,y\right),\hfill \end{array}$ (13)

其中${e^{\prime }}_{ij}=\frac{{c^{\prime }}_{ij}}{\delta p+{a^{\prime }}_{10}}$ ，${f^{\prime }}_{ij}=\frac{{d^{\prime }}_{ij}}{\delta p+{a^{\prime }}_{10}}$ (当$i+j=2$ 时)，$P_6\left(x,y\right)$ 和$Q_6\left(x,y\right)$ 是至少为$\left(x,y\right)$ 的三阶无穷可微函数。如果模型(13)的第二个方程中的$y^2$ 项系数${f^{\prime }}_{02}\ne 0$ ，即${d^{\prime }}_{02}\ne 0$ ，根据参考文献[27]可知，$E_{+}$ 是鞍结点。

通过简单计算，我们得到

${d^{\prime }}_{02}=\frac{p\left(\delta p^2{a^{\prime }}_{02}-\delta p{a^{\prime }}_{11}+\delta {a^{\prime }}_{20}-4{a^{\prime }}_{10}{b^{\prime }}_{02}\right)}{\delta p-{a^{\prime }}_{10}}=0.$

同时，根据方程(6)，我们得到$p=\frac{y_{+}}{x_{+}}$ 。也就是说，如果

$\lambda \ne \delta \left(\frac{{a^{\prime }}_{11}-\frac{x_{+}}{y_{+}}{a^{\prime }}_{20}+4\frac{x_{+}}{\delta y_{+}}{a^{\prime }}_{10}}{{a^{\prime }}_{20}}+\beta \right),$

$i.e.,\lambda \ne \frac{\delta {\left(k{y}_{+}+1\right)}^3}{k^2{x}_{+}\left(x_{+}-1\right)\left(m-x_{+}\right)}\left[-\frac{k\left(2m{x}_{+}-3x_{+}^2-m+2x_{+}\right)}{{\left(k{y}_{+}+1\right)}^2}-\frac{4{\alpha }^2{x}_{+}^2+6\alpha x_{+}+1}{{\left(\alpha x_{+}+1\right)}^3}-\frac{x_{+}\left(5m-11x_{+}+5\right)}{y_{+}\left(k{y}_{+}+1\right)}\right]+\delta \beta ,$ 则$E_{+}$ 是一个鞍结点。

情况2：$\lambda =x_{+}\frac{3\alpha x_{+}^2+2\left(1-\alpha m-\alpha \right)x_{+}+\left(\alpha m-m-1\right)}{\left(1+k{y}_{+}\right)\left(1+\alpha x_{+}\right)}+\delta \beta$ 。

在这种情况下，模型(11)的雅可比矩阵为：

$J_{E_{+}}=\left(\begin{array}{cc}{a^{\prime }}_{10}& \frac{{a^{\prime }}_{10}}{\beta -\frac{\lambda }{\delta }}\\ -\left(\beta -\frac{\lambda }{\delta }\right){a^{\prime }}_{10}& -{a^{\prime }}_{10}\end{array}\right)$ .

接下来，我们通过以下变换将模型(11)简化为标准形式：

$\left(\begin{array}{c}X\\ Y\end{array}\right)=\left(\begin{array}{cc}1& 0\\ -p& 1\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right).$

模型(11)变为以下形式

$\{\begin{array}{l}\dot{x}=c_{01}y+c_{20}{x}^2+c_{11}xy+c_{02}{y}^2+P_7\left(x,y\right),\hfill \\ \dot{y}=d_{20}{x}^2+d_{11}xy+d_{02}{y}^2+Q_7\left(x,y\right),\hfill \end{array}$ (14)

其中

$c_{01}=\frac{a_{10}}{p},c_{20}=p^2{a}_{02}-p{a}_{11}+a_{20},c_{11}=-2a_{02}p+a_{11},c_{02}=a_{02},$

$d_{20}=p^3{a}_{02}-p^2{a}_{11}+p^2{b}_{02}+p{a}_{20}-p{b}_{11}+b_{20},$

$d_{11}=-2p^2{a}_{02}+p{a}_{11}-2p{b}_{02}+b_{11},d_{02}=a_{02}p+b_{02}$

且$p=\beta -\frac{\lambda }{\delta }$ ，$P_7\left(x,y\right)$ 和$Q_7\left(x,y\right)$ 是至少为$\left(x,y\right)$ 的三阶无穷可微函数。

我们通过引入新变量$\tau =\frac{a_{10}}{p}t$ ，为了简便起见，我们仍然使用$t$ 而不是$\tau$ ，得到

$\{\begin{array}{l}\dot{x}=y+e_{20}{x}^2+e_{11}xy+e_{02}{y}^2+P_8\left(x,y\right),\hfill \\ \dot{y}=f_{20}{x}^2+f_{11}xy+f_{02}{y}^2+Q_8\left(x,y\right),\hfill \end{array}$ (15)

其中$e_{ij}=\frac{c_{ij}}{a_{10}}p$ ，$f_{ij}=\frac{d_{ij}}{a_{10}}p$ (当$i+j=2$ 时)。$P_8\left(x,y\right)$ 和$Q_8\left(x,y\right)$ 是至少为$\left(x,y\right)$ 的三阶无穷可微函数。

现在，我们根据文章[27]定理7.3的方法证明$E_{+}$ 是一个尖点。首先，我们做以下定义：

$W\left(x,y\right)\triangleq e_{20}{x}^2+e_{11}xy+e_{02}{y}^2+P_8\left(x,y\right),T\left(x,y\right)\triangleq f_{20}{x}^2+f_{11}xy+f_{02}{y}^2+Q_8\left(x,y\right).$

如果$e_{20}\ne 0$ ，根据$y+W\left(x,y\right)=0$ ，我们可以得到

$y=\varphi \left(x\right)=-e_{20}{x}^2+\cdots$

然后我们有

$\phi \left(x\right)\triangleq T\left(x,\varphi \left(x\right)\right)=f_{20}{x}^2+\cdots$

$\gamma \left(x\right)\triangleq \frac{\partial W}{\partial x}\left(x,\varphi \left(x\right)\right)+\frac{\partial T}{\partial x}\left(x,\varphi \left(x\right)\right)=\left(2e_{20}+f_{11}\right)x+\cdots$

因此，$k=2m$ ，$m=1$ ，$a_{2m}=f_{20}$ ，$N=1$ ，$B_N=2e_{20}+f_{11}$ ，$E_{+}$ 是一个尖点。

总所周知，模型(15)等价于以下模型：

$\{\begin{array}{l}\dot{x}=y,\hfill \\ \dot{y}=f_{20}{x}^2+\left(f_{11}+2e_{20}\right)xy+Q_9\left(x,y\right),\hfill \end{array}$

其中$Q_9\left(x,y\right)$ 是至少为$\left(x,y\right)$ 的三阶无穷可微函数。

$e_{20}=\frac{p}{a_{10}}\left[\frac{kp\left(3k{p}^2{x}^2+xm+px-x^2-m+x\right)}{{\left(ky+1\right)}^2\left(\alpha x+1\right)}+\frac{m-3x+1}{ky+1}+\frac{\left(px+y\right)\alpha +p}{{\left(\alpha x+1\right)}^3}\right]$

$f_{20}=\frac{p}{a_{10}}\left[\frac{k{p}^2\left(3k{p}^2{x}^2+xm+px-x^2-m+x\right)}{{\left(ky+1\right)}^2\left(\alpha x+1\right)}+\frac{p\left(m-3x+1\right)}{ky+1}+\frac{\left[\left(px+y\right)\alpha +p\right]p}{{\left(\alpha x+1\right)}^3}-\frac{4\delta p^2}{x}\right]$

$f_{11}=-\frac{p}{a_{10}}\left[\frac{kp\left(4k{p}^2{x}^2+xm+px-x^2-m+x\right)}{{\left(ky+1\right)}^2\left(\alpha x+1\right)}+\frac{p}{{\left(\alpha x+1\right)}^2}-\frac{4\delta p}{x}\right].$

此外，如果$f_{20}\ne 0$ 且$f_{11}+2e_{20}\ne 0$ ，即$\lambda \ne {\lambda }_1$ 且$\lambda \ne {\lambda }_2$ ，则$E_{+}$ 是一个具有2余维的尖点；如果$f_{20}=0$ 或$f_{11}+2e_{20}=0$ ，即$\lambda ={\lambda }_1$ 或$\lambda ={\lambda }_2$ ，则$E_{+}$ 是一个具有至少3余维的尖点，其中

${\lambda }_1=\frac{x}{4}\left(\frac{kp\left(2p^2{x}^{2}k+mx+px-x^2-m+x\right)}{{\left(ky+1\right)}^2\left(\alpha x+1\right)}+\frac{\left(px+2y\right)\alpha +p}{{\left(\alpha x+1\right)}^3}+\frac{2\left(m-3x+1\right)}{ky+1}\right)+\delta \beta ,$

${\lambda }_2=-\frac{x}{4}\left(\frac{kp\left(3k{p}^2{x}^2+xm+px-x^2-m+x\right)}{{\left(ky+1\right)}^2\left(\alpha x+1\right)}+\frac{m-3x+1}{ky+1}+\frac{\left(px+y\right)\alpha +p}{{\left(\alpha x+1\right)}^3}\right)+\delta \beta .$

5. 分岔分析

为了探讨收获项对模型(5)动力学性态的影响，我们选择$\lambda$ 作为分岔控制参数，探讨一些特定分岔动力学行为，如跨临界分岔、鞍结分岔和Hopf分岔。

5.1. 跨临界分岔

定理7：当$\lambda ={\lambda }_{TC}=\beta \delta$ 时，模型(5)经历跨临界分岔。

证明：如果$\lambda ={\lambda }_{TC}=\beta \delta$ ，则在$E_1$ 处的雅可比矩阵可以写成如下形式：

$J_{E_1}=\left(\begin{array}{cc}m-1& -\frac{1}{\alpha +1}\\ 0& 0\end{array}\right).$

然后，我们通过Sotomayor定理[28]检查跨临界分岔的条件是否成立。假设$J_{E_1}$ 和$J_{E_1}^{\text{T}}$ 对应零特征值的特征向量分别为$V$ 和$W$ ，分别为：

$V=\left(\begin{array}{c}{v}_1\\ v_2\end{array}\right)=\left(\begin{array}{c}\frac{1}{\alpha +1}\\ m-1\end{array}\right),W=\left(\begin{array}{c}{w}_1\\ w_2\end{array}\right)=\left(\begin{array}{c}0\\ 1\end{array}\right).$

经过简单计算，我们得到：

$F_{\lambda }\left(E_1;{\lambda }_{TC}\right)={\left(\begin{array}{c}0\\ -y\end{array}\right)}_{\left(E_1;{\lambda }_{TC}\right)}=\left(\begin{array}{c}0\\ 0\end{array}\right),$

$D{F}_{\lambda }\left(E_1;{\lambda }_{TC}\right)V=\left(\begin{array}{cc}0& 0\\ 0& -1\end{array}\right)\left(\begin{array}{c}\frac{1}{\alpha +1}\\ m-1\end{array}\right)=\left(\begin{array}{c}0\\ 1-m\end{array}\right),$

$\begin{array}{c}{D}^2{F}_{\lambda }\left(E_1;{\lambda }_{TC}\right)\left(V,V\right)={\left(\begin{array}{c}\frac{{\partial }^2{F}_1}{\partial x^2}{v}_1{v}_1+2\frac{{\partial }^2{F}_1}{\partial x\partial y}{v}_1{v}_2+\frac{{\partial }^2{F}_1}{\partial y^2}{v}_2{v}_2\\ \frac{{\partial }^2{F}_2}{\partial x^2}{v}_1{v}_1+2\frac{{\partial }^2{F}_2}{\partial x\partial y}{v}_1{v}_2+\frac{{\partial }^2{F}_2}{\partial y^2}{v}_2{v}_2\end{array}\right)}_{\left(E_1;{\lambda }_{TC}\right)}\\ =\left(\begin{array}{c}\frac{\left(2m-x\right)\alpha -2-2{\left(\alpha +1\right)}^2{\left(m-1\right)}^{2}k}{{\left(\alpha +1\right)}^3}\\ -2\delta {\left(m-1\right)}^2\end{array}\right).\end{array}$

因此，我们得到：

$W^{\text{T}}{F}_{\lambda }\left(E_1;{\lambda }_{TC}\right)=0,$

$W^{\text{T}}\left[D{F}_{\lambda }\left(E_1;{\lambda }_{TC}\right)V\right]=1-m\ne 0,$

$W^{\text{T}}\left[D^2{F}_{\lambda }\left(E_1;{\lambda }_{TC}\right)\left(V,V\right)\right]=-2\delta {\left(m-1\right)}^2\ne 0.$

因此，当$\lambda ={\lambda }_{TC}=\beta \delta$ 时，模型(5)将经历跨临界分岔。

5.2. 鞍结分岔

分岔面表明，当参数$\lambda$ 改变时，模型(5)将发生鞍结分岔，正平衡点的数量将从零个变为两个，下面给出定理。

定理8：当参数$\text{λ}$ 满足条件$\lambda ={\lambda }_{SN}$ 时，模型(5)经历鞍结分岔，鞍结分岔的阈值参数$\lambda ={\lambda }_{SN}$ 由$f\left(x_{h2}\right)$ 给出。

证明：通过Sotomayor定理[28]检查鞍结分岔的条件是否成立。首先，如果$\lambda ={\lambda }_{SN}$ ，则在$E_{+}$ 处的雅可比矩阵可以写成如下形式：

$J_{E_{+}}=\left(\begin{array}{cc}{a^{\prime }}_{01}p& {a^{\prime }}_{01}\\ \delta p^2& \delta p\end{array}\right),$

其中$p=\beta -\frac{\lambda }{\delta }$ ，${a^{\prime }}_{01}=-\frac{x_{+}\left(2k{y}_{+}+1\right)}{\left(k{y}_{+}+1\right)\left(\alpha x_{+}+1\right)}$ 。

显然，矩阵$J_{E_{+}}$ 的一个特征值为零。假设$J_{E_{+}}$ 对应零特征值的特征向量为$V$ 和$W$ ，则：

$V=\left(\begin{array}{c}{v}_1\\ v_2\end{array}\right)=\left(\begin{array}{c}1\\ -p\end{array}\right),W=\left(\begin{array}{c}{w}_1\\ w_2\end{array}\right)=\left(\begin{array}{c}{a^{\prime }}_{01}\\ \delta p\end{array}\right).$

因此，我们有：

$F_{\lambda }\left(E_{+};{\lambda }_{SN}\right)={\left(\begin{array}{c}0\\ -y\end{array}\right)}_{\left(E_{+};{\lambda }_{SN}\right)}=\left(\begin{array}{c}0\\ -y_{+}\end{array}\right),$

$\begin{array}{l}{D}^2{F}_{\lambda }\left(E_{+};{\lambda }_{SN}\right)\left(V,V\right)\\ ={\left(\begin{array}{c}\frac{{\partial }^2{F}_1}{\partial x^2}{v}_1{v}_1+2\frac{{\partial }^2{F}_1}{\partial x\partial y}{v}_1{v}_2+\frac{{\partial }^2{F}_1}{\partial y^2}{v}_2{v}_2\\ \frac{{\partial }^2{F}_2}{\partial x^2}{v}_1{v}_1+2\frac{{\partial }^2{F}_2}{\partial x\partial y}{v}_1{v}_2+\frac{{\partial }^2{F}_2}{\partial y^2}{v}_2{v}_2\end{array}\right)}_{\left(E_1;{\lambda }_{SN}\right)}\\ =\left(\begin{array}{c}\frac{2kp\left(3k{p}^2{x}^2+mx+px-x^2-m+x\right)}{{\left(k\left(\alpha x+1\right)y\left(\alpha x+1\right)+1\right)}^2}+\frac{\left(2px+y\right)\alpha +2p}{{\left(\alpha x+1\right)}^3}+\frac{m-3x+1}{ky+1}\\ -\frac{6\delta p^2}{x_{+}}\end{array}\right).\end{array}$

最后，我们得到：

$W^{\text{T}}{F}_{\lambda }\left(E_{+};{\lambda }_{SN}\right)=-\delta p{y}_{+}\ne 0,$

$\begin{array}{l}{W}^{\text{T}}\left[D^2{F}_{\lambda }\left(E_{+};{\lambda }_{SN}\right)\left(V,V\right)\right]\\ ={a^{\prime }}_{01}\left(\frac{2kp\left(3k{p}^2{x}^2+mx+px-x^2-m+x\right)}{{\left(k\left(\alpha x+1\right)y\left(\alpha x+1\right)+1\right)}^2}+\frac{\left(2px+y\right)\alpha +2p}{{\left(\alpha x+1\right)}^3}+\frac{m-3x+1}{ky+1}\right)-\frac{6{\delta }^2{p}^3}{x_{+}}\ne 0.\end{array}$

因此，当参数$\lambda$ 满足条件$\lambda ={\lambda }_{SN}$ 时，模型(5)经历鞍结分岔。

5.3. Hopf分岔

基于前面定理可知，当$E_{+}^{*}$ 存在时，它始终是鞍点，而$E_{+}^{**}$ 可能是汇点、源点或中心。考虑$\lambda$ 为Hopf分岔控制参数，Hopf分岔的阈值参数由$Tr\left(J_{E_{+}^{**}}\right)=0$ 给出，即$\lambda ={\lambda }_H$ 。同时，结合方程(10)和(4)，我们发现${Det\left(J_{E_{+}^{**}}\right)|}_{\lambda ={\lambda }_H}>0$ 。随着$\lambda$ 的变化并越过阈值${\lambda }_H$ ，$E_{+}^{**}$ 的稳定性发生变化。接下来我们将证明这一结论。

定理9：如果$f\left(x_{h2}\right)<0$ ，则模型(5)有两个正平衡点$E_{+}^{*}$ 和$E_{+}^{**}$ ，其中$E_{+}^{*}$ 是鞍点，且$x_{+}^{**}$ 是单一的正根，则模型(5)将在$\lambda ={\lambda }_H$ 附近经历Hopf分岔。

证明：我们已经得到当$\lambda ={\lambda }_H$ 时，${Det\left(J_{E_{+}^{**}}\right)|}_{\lambda ={\lambda }_H}>0$ 且$Tr\left(J_{E_{+}^{**}}\right)=0$ 。因此，我们只需要检查Hopf分岔的横向条件是否成立。即：

$\frac{\text{d}}{\text{d}\lambda }Tr\left(J_{E_{+}^{**}}\right)=1-{\frac{\text{d}}{\text{d}\lambda }\left[x_{+}^{**}\frac{3\alpha x_{+}^{**2}+2\left(1-\alpha m-\alpha \right)x_{+}^{**}+\left(\alpha m-m-1\right)}{\left(1+k{y}_{+}^{**}\right)\left(1+\alpha x_{+}^{**}\right)}\right]|}_{\lambda ={\lambda }_H}\ne 0$

通过我们的数值模拟，横向条件得到了满足，这意味着$E_{+}^{**}$ 在$\lambda ={\lambda }_H$ 处通过Hopf分岔失去其稳定性。

接下来，需要确定极限环的方向和稳定性，然后计算第一Lyapunov数。

通过将$E_{+}^{**}$ 变换为原点，我们将模型(5)转化为如下形式：

$\{\begin{array}{l}\dot{u}={\alpha }_{10}u+{\alpha }_{01}v+{\alpha }_{11}uv+a_{20}{u}^2+a_{02}{v}^2+{\alpha }_{30}{u}^3+{\alpha }_{21}{u}^{2}v+{\alpha }_{12}u{v}^2+{\alpha }_{03}{v}^3+P_{10}\left(u,v\right)\\ \dot{v}={\beta }_{10}u+{\beta }_{01}v+{\beta }_{11}uv+{\beta }_{20}{u}^2+{\beta }_{02}{v}^2+{\beta }_{30}{u}^3+{\beta }_{21}{u}^{2}v+{\beta }_{12}u{v}^2+{\beta }_{03}{v}^3+Q_{10}\left(u,v\right)\end{array}$

其中

${\alpha }_{10}=-\frac{x_{+}^{**}\left[3\alpha x_{+}^{**2}+2\left(1-\alpha m-\alpha \right)x_{+}^{**}+\left(\alpha m-m-1\right)\right]}{\left(k{y}_{+}^{**}+1\right)\left(\alpha y_{+}^{**}+1\right)},{\alpha }_{01}=-\frac{x_{+}^{**}\left(2k{y}_{+}^{**}+1\right)}{\left(k{y}_{+}^{**}+1\right)\left(\alpha x_{+}^{**}+1\right)},$

${\alpha }_{20}=\frac{\left(m-3x_{+}^{**}+1\right)}{k{y}_{+}^{**}+1}-\frac{y_{+}^{**}\alpha }{{\left(\alpha x_{+}^{**}+1\right)}^3},{\alpha }_{11}=-\frac{\left(2x_{+}^{**}m-3x_{+}^{**2}-m+2x_{+}^{**}\right)k}{{\left(k{y}_{+}^{**}+1\right)}^2}+\frac{1}{{\left(\alpha x_{+}^{**}+1\right)}^2},$

${\alpha }_{02}=\frac{x_{+}^{**}{k}^2\left(x_{+}^{**}m-x_{+}^{**2}-m+x_{+}^{**}\right)}{{\left(k{y}_{+}^{**}+1\right)}^3},{\alpha }_{30}=-\frac{y_{+}^{**}{\alpha }^2}{{\left(\alpha x_{+}^{**}+1\right)}^4}-\frac{1}{k{y}_{+}^{**}+1},$

${\alpha }_{21}=\frac{\alpha }{{\left(\alpha x_{+}^{**}+1\right)}^3}-\frac{k\left(m-3x_{+}^{**}+1\right)}{{\left(k{y}_{+}^{**}+1\right)}^2},{\alpha }_{12}=\frac{k^2\left(2x_{+}^{**}m-3x_{+}^{**2}-m+2x_{+}^{**}\right)}{{\left(k{y}_{+}^{**}+1\right)}^3},$

${\alpha }_{03}=-\frac{k^3{x}_{+}^{**}{y}_{+}^{**}}{\left(\alpha x_{+}^{**}+1\right){\left(k{y}_{+}^{**}+1\right)}^3},$

${\beta }_{10}=\frac{\delta y_{+}^{**2}}{x_{+}^{**2}},{\beta }_{01}=\frac{\delta y_{+}^{**}}{x_{+}^{**}},{\beta }_{20}=-\frac{\delta y_{+}^{**2}}{x_{+}^{**3}},{\beta }_{11}=\frac{2\delta y_{+}^{**}}{x_{+}^{**2}},{\beta }_{02}=-\frac{\delta }{x_{+}^{**}},$

${\beta }_{30}=\frac{\delta y_{+}^{**2}}{x_{+}^{**4}},{\beta }_{21}=-\frac{2\delta y_{+}^{**}}{x_{+}^{**3}},{\beta }_{12}=\frac{\delta }{x_{+}^{**2}},{\beta }_{03}=0,$

且$P_{10}\left(u,v\right)$ 和$Q_{10}\left(u,v\right)$ 是定义在$\left(u,v\right)$ 上至少四阶光滑函数。

第一Lyapunov数为

$\begin{array}{c}{l}_1=\frac{3\pi }{2{\alpha }_{01}det{\left(J_{E_{+}^{**}}\right)}^{\frac{3}{2}}}[-\frac{3}{k{y}_{+}^{**}+1}+\frac{\alpha y_{+}^{**}}{x_{+}^{**3}}+\frac{k\left(2m{x}_{+}^{**}-3x_{+}^{**2}-m+2x_{+}^{**}\right)\left(k-m+3x_{+}^{**}-1\right)}{{\left(k{y}_{+}^{**}+1\right)}^3}\\ -\frac{2k^2\delta y_{+}^{**}}{\left(\alpha x_{+}^{**}+1\right){\left(k{y}_{+}^{**}+1\right)}^2}+\frac{-m+3x_{+}^{**}-1}{{\left(\alpha x_{+}^{**}+1\right)}^2\left(k{y}_{+}^{**}+1\right)}+2y_{+}^{**2}\delta \frac{\left(m-3x_{+}^{**}+1\right){\left(\alpha x_{+}^{**}+1\right)}^3+\alpha y_{+}^{**}\left(k{y}_{+}^{**}+1\right)}{x_{+}^{**3}\left(k{y}_{+}^{**}+1\right){\left(\alpha x_{+}^{**}+1\right)}^3}\\ -\frac{k^2{x}_{+}^{**}{y}_{+}^{**}}{{\left(\alpha x_{+}^{**}+1\right)}^2{\left(k{y}_{+}^{**}+1\right)}^2}+\frac{\delta \left(2\delta x_{+}^{**2}{y}_{+}^{**}+2\delta y_{+}^{**3}+x_{+}^{**3}\right)}{x_{+}^{**5}}-\frac{k^3{x}_{+}^{**}{y}_{+}^{**}\left(2m{x}_{+}^{**}-3x_{+}^{**2}-m+2x_{+}^{**}\right)}{\left(\alpha x_{+}^{**}+1\right){\left(k{y}_{+}^{**}+1\right)}^4}\\ +\frac{\alpha k^2{x}_{+}^{**2}{y}_{+}^{**}-k{y}_{+}^{**}\alpha \left(2m{x}_{+}^{**}-3x_{+}^{**2}-m+2x_{+}^{**}\right)}{{\left(\alpha x_{+}^{**}+1\right)}^3{\left(k{y}_{+}^{**}+1\right)}^2}-\frac{3\alpha y_{+}^{**}\left(x_{+}^{**}{\alpha }^2+\alpha +\frac{1}{3}\right)}{{\left(\alpha x_{+}^{**}+1\right)}^5}].\end{array}$

基于上述理论，我们获得模型(5)发生跨临界分岔、鞍结分岔和Hopf分岔的阈值条件，为后续研究种群生长共存模式及其内在驱动机制提供了一定的理论基础。同时，也从理论上揭示了一些关键参数取值大小严重影响模型(5)的特定分岔动力学性态。

6. 数值模拟分析

为了验证理论推导结果的有效性和可行性，探索关键参数如何影响模型(5)的特定分岔动力学性态，对模型(5)进行相关动力学模拟。基于各类关键参数的生物学意义和以往发表论文文献，并结合相关定理给出的阈值条件，当我们选取参数$\lambda$ 作为控制参数，其它参数取值可为$k=0.88$ ，$m=0.3$ ，$\alpha =0.6$ ，$\delta =0.8$ ，$\beta =1$ 。依据定理7可知${\lambda }_{TC}=0.8$ 。由图3可知，当$\lambda =0.83>{\lambda }_{TC}=0.8$ 时，模型(5)拥有两个边界平衡点，且分别为鞍点与稳定的结点，这暗示了捕食者种群与食饵种群不能形成生长共存模式，捕食者种群最终走向灭绝。当$\lambda =0.77<{\lambda }_{TC}=0.8$ 时，模型(5)拥有两个边界平衡点和两个内平衡点，这暗示了捕食者种群与食饵种群能形成生长共存模式。同时，从图3也可知模型(5)发生一次跨临界分岔动力学行为，导致模型(5)中的种群可以生长共存模式。

Figure 3. The dynamic process of transcritical bifurcation in the model (5)

图3. 模型(5)的跨临界分岔动态过程

通过观察图4可知，模型(5)发生了一次鞍结分岔，可以产生一个稳定的内平衡点，这暗示了捕食者种群与食饵种群可以形成一个常稳态生长共存模式。因此，从图3和图4可知，由控制参数取值所诱发的模型(5)跨临界分岔与鞍结分岔是导致种群生长共存模式形成的内在驱动机制之一。

Figure 4. Dynamic process of saddle node bifurcation in the model (5)

图4. 模型(5)的鞍结分岔动态过程

当我们选取参数$\lambda$ 作为控制参数，其它参数取值为$k=0.2$ ，$m=0.3$ ，$\alpha =0.05$ ，$\delta =0.3$ ，$\beta =0.4$ ，依据定理9可知${\lambda }_H=0.0597360$ 。从图5可知，模型(5)会发生Hopf分岔，产生极限环。同时计算第一Lyapunov数，得到$l_1=17843.436399956\pi >0$ ，因此可得模型(5)具有一个不稳定的极限环，这暗示了捕食者种群与食饵种群可以形成一个不稳定的周期振荡生长共存模式。

Figure 5. The Hopf bifurcation dynamic process in the model (5)

图5. 模型(5)的Hopf分岔动态过程

基于数值模拟结果可知，随着控制参数取值的变化，模型(5)会发生跨临界分岔、鞍结分岔和Hopf分岔动力学行为，进而导致捕食者种群与食饵种群可以形成常稳态生长共存模式和不稳定的周期振荡生长共存模式，这也直接揭示了模型特定分岔动力学行为是种群生长共存模式产生的内在驱动机制。

7. 总结

本论文考虑食饵种群的Allee和恐惧效应，并对捕食者种群进行收获，提出了一类具有线性收获和恐惧效应的Leslie-Gower捕食生态模型，证明了模型(5)所有可能平衡点的存在性和稳定性，推导出模型(5)发生跨临界分岔、鞍结分岔和Hopf分岔的阈值条件，这些理论条件不仅为后续数值模拟提供参数取值的理论基础，而且直接指出了关键参数取值大小深刻影响模型的动力学性态和种群持久生存模式。

对模型(5)的特定分岔动力学性态进行数值模拟，从图3可知，模型(5)发生跨临界分岔动力学行为，这可以导致捕食者种群从逐渐灭绝状态转变为逐渐短暂生存状态。从图4可知，模型(5)发生鞍结分岔动力学行为，这迫使食饵种群与捕食者种群形成常数稳态生长共存模式。从图5可知，模型(5)发生Hopf分岔动力学行为，这驱动食饵种群与捕食者种群以周期振荡模式持久生存。这些数值模拟结果显示随着参数$\lambda$ 取值的变化，模型(5)的分岔动力学性态发生了本质的改变，进而驱使模型(5)所代表的捕食生态系统从不持续生存状态变为可持续生存状态，这不仅代表捕食生态系统具有更强的内在自我调节韧性，并揭示人为调控措施是维持捕食生态系统持久生存的外在调控机制。

总而言之，希望本论文的研究结果有利于拓展捕食生态模型动力学性态的研究体系。

基金项目

国家自然科学基金(31570364)、浙江省大学生科技创新活动计划(新苗人才计划) (2024R429A010和2025R444A017)、温州大学国家级大学生创新创业训练计划项目(JWXC2024081)。

NOTES

^*通讯作者。