一类具有随机环境扰动的捕食者–食饵生态模型动力学问题研究

Research on Dynamics of a Predator-Prey Ecological Model with Stochastic Environmental Disturbances

摘要: 本文提出了一类带有随机环境扰动的捕食者–食饵种群生态模型，对其特定动力学行为进行理论分析与数值模拟。数学理论工作集中讨论了模型解的全局正性、随机最终有界性，同时分析了种群灭绝、持久生存和系统稳态分布的充分条件。数值模拟结果发现，环境白噪声强度是影响捕食者和食饵种群存活与灭绝的关键因素，从概率动力学的角度揭示了种群生长共存模式及其内在驱动机制，为拓宽捕食生态模型复杂动力学问题的研究框架提供一定的理论支撑。

Abstract: In this paper, we proposed a predator-prey population ecological model with stochastic environ- mental disturbances, and conducted theoretical analysis and numerical simulation of its specific dynamic behavior. The mathematical theory work focused on the global positivity and stochastic ultimate boundedness of model solutions, and analyzed the sufficient conditions for population extinction, persistent survival, and system steady-state distribution. Numerical simulations revealed that the intensity of environmental white noise was a key factor affecting the survival and extinction of both predator and prey populations. From the perspective of probabilistic dynamics, this study uncovered the population growth coexistence mode and their intrinsic driving mechanisms, which could provide some theoretical support for expanding the research framework of complex dynamical problems in predator-prey ecological models.

1. 引言

自然种群生态系统不可避免地受到环境扰动的影响，这源于栖息地固有的随机性与不可预测性。温度波动、降雨量变异、台风等非生物因素，会直接作用于种群的生存与繁衍过程，进而影响种群生长动态。事实上，在环境多变性与随机过程的双重驱动下，自然种群的生物量通常围绕长期均值上下波动，而非稳定于某一固定平衡态[1]。大量实验证实，环境随机扰动对种群生态系统的稳定性具有显著调控作用[2] [3]。研究进一步表明，在持续变化的生境条件下，这类随机干扰甚至会对生态种群的生存产生决定性影响[4]。值得注意的是，环境波动的调控效应可通过生态模型参数得以量化。相关研究表明，环境波动可显著影响模型中的增长率、死亡率和竞争系数等关键参数[5]，而将这些模型参数视为服从一定分布的随机变量，更能够精准描述参数变化的不确定性，从而更贴合自然种群生态系统的实际动态[6]。因此，研究随机种群生态系统比确定性模型更具实际价值，这也使得随机种群模型成为生态学领域的研究热点，受到学者的广泛关注[7]-[9]。

在种群生态模型的构建中，向确定性模型引入随机性的方式具有多样性，其中两种核心手段因适配性强、贴合实际研究场景，被广泛应用于相关研究。第一种常用手段为参数扰动法，旨在通过修正模型参数来量化环境变异性对种群动态的影响，具体操作是将确定性模型中固定不变的参数(如种群增长率、死亡率)，替换为该参数的均值与随机误差项之和[10] [11]。从理论层面来看，根据中心极限定理，这类随机误差项会呈现正态分布特征，而该分布与自然环境中普遍存在的高斯白噪声具有高度适配性，因此可近似为高斯白噪声[12]，目前该参数扰动策略已在各类生物种群模型中获得广泛应用与验证[13]-[15]。另一种核心手段为方程扰动法，区别于参数层面的调整，其通过直接在种群生长方程中嵌入随机扰动项实现随机性引入，且扰动强度的大小与状态变量(如种群生物量)偏离其确定性平衡态的程度正相关[16]。该方法的优势在于，不仅能够助力研究人员探究确定性动力学系统的恢复能力，还可深入解析系统长期演化中的随机行为特征[17]，因其在生物系统(尤其是种群生态系统)研究中展现出显著的有效性与实用性，已成为学界的研究热点，受到广泛关注[18] [19]。两种方法互补适配，为揭示随机环境下种群动态规律提供了重要技术支撑。

2. 模型构建

本文聚焦环境波动对种群动态的调控作用，假设环境随机波动可直接影响捕食者与食饵种群的死亡率。为量化该调控效应，本文借鉴文献[13] [14]的随机建模思路，对文献[20]提出的确定性捕食者–食饵模型进行修正，将模型中两类种群的恒定死亡项$-m_1$ 、$-m_2$ 分别替换为$-m_1+{\sigma }_1\text{d}{B}_1\left(t\right)$ 与$-m_2+{\sigma }_2\text{d}{B}_2\left(t\right)$ ，以此嵌入环境随机扰动。其中，$\text{d}{B}_i\left(t\right)$ 为布朗运动$B_i\left(t\right)\left(i=1,2\right)$ 的微分，对应环境白噪声扰动项；${\sigma }_i^2\left(i=1,2\right)$ 表征噪声强度，反映环境波动的剧烈强度；布朗运动$B_i\left(t\right)\left(i=1,2\right)$ 相互独立，定义于完备概率空间$\left(\Omega ,\mathbb{F},\mathbb{P}\right)$ ，且关于右连续滤过${\left\{{\mathbb{F}}_t\right\}}_{t\in {ℝ}_{+}}$ 满足标准布朗运动的概率性质[21]。

基于以上分析，本节提出如下环境随机扰动下的捕食者–食饵动力学模型：

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=rx\left(t\right)\frac{1+ax\left(t\right)}{1+ax\left(t\right)+by\left(t\right)}-m_{1}x\left(t\right)-c{x}^2\left(t\right)-\frac{\alpha x\left(t\right)y\left(t\right)}{\beta +x^2\left(t\right)}+{\sigma }_{1}x\left(t\right)\text{d}{B}_1\left(t\right),\\ \frac{\text{d}y}{\text{d}t}=\frac{\lambda \alpha x\left(t\right)y\left(t\right)}{\beta +x^2\left(t\right)}-\frac{qey\left(t\right)}{c_{1}e+c_{2}y\left(t\right)}-m_{2}y\left(t\right)+{\sigma }_{2}y\left(t\right)\text{d}{B}_2\left(t\right).\end{array}$ (1.1)

其中，$r$ 是食饵种群的繁殖项；$b$ 是恐惧程度；$a$ 是恐惧遗留效应程度；$c$ 表示食饵种群的内在竞争系数；$m_1$ 是食饵种群的自然死亡率；$\alpha$ 表示捕食者种群对食饵种群的捕获率；$\beta$ 表示半饱和系数；$\lambda$ 表示捕食者种群转化食饵种群为自身生长的能量转化率；$q$ 表示捕食者种群可捕获性系数；$e$ 表示用于收获捕食者种群的捕获努力量；$c_1$ 表示用于衡量各类技术手段在捕获中所占的比重，$c_2$ 表示捕获率和处理率；$m_2$ 表示捕食者种群的自然死亡率。

3. 数学理论推导

3.1. 预备知识

记${ℝ}_{+}=\left[0,+\infty \right)$ ，${ℝ}_{+}^n=\left\{\left(x_1,\cdots ,x_n\right)\in {ℝ}^n:x_i>0,i=1,2,\cdots ,n\right\}$ ，且$\left|x\right|=\sqrt{{\displaystyle {\sum }_{i=1}^n{x}_i^2}}$ 。除非另有说明，我们始终假设$\left(\Omega ,{\mathbb{F}}_t,{\left\{{\mathbb{F}}_t\right\}}_{t\ge 0},\mathbb{P}\right)$ 是一个满足通常条件的完备概率空间(即滤过${\left\{{\mathbb{F}}_t\right\}}_{t\ge 0}$ 满足右连续性，且${\mathbb{F}}_0$ 包含所有$\mathbb{P}$ -零集)。

一般地，考虑如下形式的$n$ 维随机微分方程：

$\text{d}x\left(t\right)=f\left(x\left(t\right),t\right)\text{d}t+g\left(x\left(t\right),t\right)\text{d}B\left(t\right),t\in \left[t_0,T\right]$ (1.2)

初始值为$x\left(t_0\right)=x_0\in {ℝ}^n$ ，其中，$B\left(t\right)$ 是定义在完备概率空间$\left(\Omega ,\mathbb{F},\mathbb{P}\right)$ 上的$n$ 维标准布朗运动(也称维纳过程)。记$C^{2,1}\left({ℝ}^n\times {ℝ}_{+},ℝ\right)$ 为所有定义在${ℝ}^n\times {ℝ}_{+}$ 上的非负函数$V\left(x,t\right)$ 构成的函数族，且$V\left(x,t\right)$ 关于$x$ 二阶连续可微，关于$t$ 一阶连续可微。定义与方程(1.2)相关的微分算子$ℒ$ 如下[22]：

$ℒ=\frac{\partial }{\partial t}+{\displaystyle \sum }_{i=1}^n{f}_i\left(x,t\right)\frac{\partial }{\partial x_i}+\frac{1}{2}{\displaystyle \sum }_{i,j=1}^n{\left[g^{\text{T}}\left(x,t\right)g\left(x,t\right)\right]}_{ij}\frac{{\partial }^2}{\partial x_i\partial x_j}.$

设$V\left(x,t\right)\in C^{2,1}\left({ℝ}^n\times {ℝ}_{+},ℝ\right)$ ，则

$ℒV=V_t+V_{x}f\left(x,t\right)+\frac{1}{2}trace\left[g^{\text{T}}\left(x,t\right)V_{xx}g\left(x,t\right)\right],$

其中，

$V_t=\frac{\partial V}{\partial t},V_x=\left(\frac{\partial V}{\partial x_1},\frac{\partial V}{\partial x_2},\cdots ,\frac{\partial V}{\partial x_n}\right),V_{xx}={\left(\frac{{\partial }^2}{\partial x_i\partial x_j}\right)}_{n\times n}.$

根据伊藤公式，若$x\left(t\right)\in {ℝ}^n$ ，则

$\text{d}V=LV\left(x\left(t\right),t\right)\text{d}t+V_x\left(x\left(t\right),t\right)g\left(x\left(t\right),t\right)\text{d}B\left(t\right).$

接下来我们引入平稳分布的判别准则。在此之前，考虑如下随机微分方程：

$x\left(t\right)=f\left(x\left(t\right)\right)\text{d}t+{\displaystyle \sum }_{r=1}^k{\sigma }_r\left(x\right)\text{d}{B}_r\left(t\right),$ (1.3)

其中，$x\left(t\right)$ 是$n$ 维欧几里得空间${ℝ}^n$ 上的齐次马尔可夫过程。扩散矩阵为$D\left(x\right)=\left(d_{ij}\left(x\right)\right)$ ， $d_{ij}\left(x\right)={\displaystyle {\sum }_{r=1}^k{\sigma }_r^i\left(x\right){\sigma }_r^j\left(x\right)}$ 。因此，以下引理给出平稳分布的判别准则。

引理1 [23] 若存在有界开集$E\subset {ℝ}^n$ ，其边界$\Gamma$ 光滑且正则，满足以下条件：

(i) 对所有$x\in E$ ，扩散矩阵$D\left(x\right)$ 严格正定；

(ii) 存在非负$C^2$ 函数$V\left(x\right)$ 与正常数$M$ ，使得对任意$x\in {ℝ}^n/E$ ，有$ℒV\le -M$ 。

则方程(1.3)的解$x\left(t\right)$ 是一个具有平稳分布$\mu (\cdot )$ 的平稳马尔可夫过程，且对测度$\mu$ 可积的任意函数$g(\cdot )$ ，有

$\mathbb{P}\left(\underset{t\to \infty }{\lim }\frac{1}{t}{\displaystyle {\int }_0^{t}g\left(x\left(t\right)\right)\text{d}t}={\displaystyle {\int }_{{ℝ}^n}g\left(x\right)\mu \left(\text{d}x\right)}\right)=1.$

定义1 [24] 设$x\left(t\right)$ 为方程(1.2)的解，则

(1) 若$\lim {\text{sup}}_{t\to \infty }\frac{1}{t}{\displaystyle {\int }_0^{t}x\left(s\right)\text{d}s}<0$ ，称$x\left(t\right)$ 灭绝；

(2) 若$\lim {\text{sup}}_{t\to \infty }\frac{1}{t}{\displaystyle {\int }_0^{t}x\left(s\right)\text{d}s}=0$ ，称$x\left(t\right)$ 均值非持久；

(3) 若$\lim {\text{sup}}_{t\to \infty }\frac{1}{t}{\displaystyle {\int }_0^{t}x\left(s\right)\text{d}s}>0$ ，称$x\left(t\right)$ 均值弱持久；

(4) 若$\lim {\text{inf}}_{t\to \infty }\frac{1}{t}{\displaystyle {\int }_0^{t}x\left(s\right)\text{d}s}>0$ ，称$x\left(t\right)$ 均值强持久。

3.2. 动力学分析

在本节中，我们主要研究系统(1.1)全局正解的存在唯一性与随机最终有界性、捕食者–食饵种群的灭绝与均值持久性行为，以及解的稳态分布。

3.2.1. 系统(1.1)的全局正解的存在性和唯一性

在探究系统(1.1)的动力学行为之前，首要任务是明确随机解的非负性与全局存在性，这是开展后续研究的基础性结论。为此，本小节将采用李雅普诺夫分析法，讨论系统(1.1)全局正解的存在性问题。基于生态学意义，我们仅考虑系统(1.1)的非负解，并给出如下结论。

定理3.1 对任意给定的初始值$\left(x\left(0\right),y\left(0\right)\right)\in {ℝ}_{+}^2$ ，系统(1.1)在${ℝ}_{+}$ 上存在唯一解$\left(x\left(t\right),y\left(t\right)\right)$ ，且该正解以概率1保持在${ℝ}_{+}^2$ 内，即对所有$t>0$ ，几乎处处有$\left(x\left(t\right),y\left(t\right)\right)\in {ℝ}_{+}^2$ 。

证明：根据文献[25]中引理2.1的方法，显然可知，对任意给定的初始值$\left(x\left(0\right),y\left(0\right)\right)\in {ℝ}_{+}^2$ ，系统(3.1)的所有系数均满足局部利普希茨连续，因此该系统在区间$t\in \left[0,{\tau }_e\right)$ 上存在唯一的局部解$\left(x\left(t\right),y\left(t\right)\right)$ ，其中${\tau }_e$ 表示爆破时间。接下来，我们只需证明该解是全局解，即证明${\tau }_e=\infty$ (几乎处处)。取充分大的整数$n_0\ge 1$ ，使得$\left(x_0,y_0\right)\in \left[\frac{1}{n_0},n_0\right]$ 。对每个整数$n\ge n_0$ ，定义如下停时：

${\tau }_n=\inf \left\{t\in [0,{\tau }_e):\min \left\{\left(x\left(t\right),y\left(t\right)\right)\right\}\le \frac{1}{n}\text{or}\min \left\{\left(x\left(t\right),y\left(t\right)\right)\right\}\ge n\right\},$

并规定$\inf \varnothing =\infty$ ($\varnothing$ 表示空集)。显然，${\tau }_n$ 随$n\to \infty$ 单调递增。记${\tau }_{\infty }=\underset{n\to +\infty }{\lim }{\tau }_n$ ，则有${\tau }_{\infty }\le {\tau }_e$ (几乎处处)。若能证明${\tau }_{\infty }=\infty$ (几乎处处)，则可得${\tau }_e=\infty$ ，且对所有$t\ge 0$ ，几乎处处有$\left(x\left(t\right),y\left(t\right)\right)\in {ℝ}_{+}^2$ 。换句话说，完成证明只需验证${\tau }_{\infty }=\infty$ (几乎处处)。若上述结论不成立，则存在常数$T>0$ 与$ϵ\in \left(0,1\right)$ ，使得$\mathbb{P}\left\{{\tau }_{\infty }\le T\right\}>ϵ$ 。因此，存在整数$n_1\ge n_0$ ，使得对所有$n\ge n_1$ ，

$\mathbb{P}\left\{{\tau }_{\infty }\le T\right\}\ge ϵ,$ (1.4)

定义如下$C^2$ 函数$\overline{V}:{ℝ}_{+}^2\to {ℝ}_{+}$ ：

$\overline{V}\left(x,y\right)=\lambda \left(x-1-\log x\right)+\left(y-1-\log y\right).$

显然，函数$\overline{V}\left(x,y\right)$ 是非负的。对函数$\overline{V}\left(x,y\right)$ 应用伊藤公式，可得

$\text{d}\overline{V}\left(x,y\right)=ℒ\overline{V}\left(x,y\right)\text{d}t+\lambda {\sigma }_1\left(x-1\right)\text{d}{B}_1\left(t\right)+{\sigma }_2\left(y-1\right)\text{d}{B}_2\left(t\right),$ (1.5)

其中，微分算子$ℒ\overline{V}:{ℝ}_{+}^2\to {ℝ}_{+}$ 的定义为：

$\begin{array}{c}ℒ\overline{V}\left(x,y\right)=\left(x-1\right)\left(r\frac{1+ax}{1+ax+by}-cx-\frac{\alpha y}{\beta +x^2}-m_1\right)+\frac{{\sigma }_1^2}{2}\\ +\left(y-1\right)\left(\frac{\lambda \alpha x}{\beta +x^2}-\frac{qe}{c_{1}e+c_{2}y}-m_2\right)+\frac{{\sigma }_2^2}{2}\\ \le m_2+\frac{{\sigma }_2^2}{2}+\left(m_1-r+\frac{{\sigma }_1^2}{2}\right)-c{x}^2+\left(r+c-m_1\right)x+m_{2}y\\ \le K+m_{2}y,\end{array}$

其中，$K=m_2+\frac{{\sigma }_2^2}{2}+\left(m_1-r+\frac{{\sigma }_1^2}{2}\right)-\frac{c^2}{4\left(r+c-m_1\right)}$ 。

注意到对所有$y>0$ ，$y\le 2\left(x-1-\log x\right)+2\log 2\le 2\overline{V}\left(x,y\right)+2\log 2$ ，则

$ℒ\overline{V}\left(x,y\right)\le K+2m_2\log 2+2m_2\overline{V}\le \Upsilon \left(1+\overline{V}\right),$ (1.6)

其中，$\Upsilon =\max \left\{K+2m_2\log 2+2\log 2\right\}$ 。

根据(1.5)与(1.6)，可得

$ℒ\overline{V}\left(x,y\right)\le \Upsilon \left(1+\overline{V}\right)+\lambda {\sigma }_1\left(x-1\right)\text{d}{B}_1\left(t\right)+{\sigma }_2\left(y-1\right)\text{d}{B}_2\left(t\right)$ (1.7)

对(1.7)的两边从0到${\tau }_n\wedge T$ ，并取期望

$\mathbb{E}\overline{V}\left(x\left({\tau }_n\wedge T\right),y\left({\tau }_n\wedge T\right)\right)\le \overline{V}\left(x_0,y_0\right)+\Upsilon {\displaystyle {\int }_0^{{\tau }_n\wedge T}\mathbb{E}\overline{V}\left(x,y\right)\text{d}t}+\Upsilon T.$

由Gronwall不等式可得

$\mathbb{E}\overline{V}\left(x\left({\tau }_n\wedge T\right),y\left({\tau }_n\wedge T\right)\right)\le \left[\overline{V}\left(x_0,y_0\right)+\Upsilon T\right]\exp \left(\Upsilon T\right)$ (1.8)

根据(1.4)，对一切$n\ge n_1$ ，我们有$\mathbb{P}\left({\Omega }_n\right)\ge ϵ$ ，其中${\Omega }_n=\left\{{\tau }_n\le T\right\}$ 。因此，对每一个$\chi \in {\Omega }_n$ ，至少有一个$x\left({\tau }_n\wedge T\right)$ 或$y\left({\tau }_n\wedge T\right)$ 等于$n$ 或$\frac{1}{n}$ 。故$\overline{V}\left(x\left({\tau }_n\wedge T\right),y\left({\tau }_n\wedge T\right)\right)$ 不小于$\min \left\{n-1-\log n,\frac{1}{n}-1+\log n\right\}$ 。所以，

$\overline{V}\left(x\left({\tau }_n\wedge T\right),y\left({\tau }_n\wedge T\right)\right)\ge \left(n-1-\log n\right)\wedge \left(\frac{1}{n}-1+\log n\right).$

根据(1.8)，则

$\begin{array}{c}\left[\overline{V}\left(x_0,y_0\right)+\Upsilon T\right]\exp \left(\Upsilon T\right)\ge \mathbb{E}\left[I_{{\Omega }_n\left(\chi \right)}\overline{V}\left(x\left({\tau }_n,\chi \right),y\left({\tau }_n,\chi \right)\right)\right]\\ \ge ϵ\left[\left(n-1-\log n\right)\wedge \left(\frac{1}{n}-1+\log n\right)\right],\end{array}$

其中，$I_{{\Omega }_n\left(\chi \right)}$ 是${\Omega }_n$ 的特征函数。当$n\to +\infty$ 时，有

$\infty >\left[\overline{V}\left(x_0,y_0\right)+\Upsilon T\right]\exp \left(\Upsilon T\right)\ge ϵ\left[\left(n-1-\log n\right)\wedge \left(\frac{1}{n}-1+\log n\right)\right]=\infty ,$

这与假设相矛盾。因此，${\tau }_{\infty }=\infty$ (几乎处处)成立。

证毕。

3.2.2. 系统(1.1)解的随机最终有界性

定理3.2 对任意给定的初始值$\left(x\left(0\right),y\left(0\right)\right)\in {ℝ}_{+}^2$ ，则系统(1.1)的解$\left(x\left(t\right),y\left(t\right)\right)$ 满足

$\underset{t\to +\infty }{\lim }\sup \left[x\left(t\right)+\frac{1}{\lambda }y\left(t\right)\right]<+\infty ,a.s.$

此外，该解是随机最终有界的，即存在一对常数$\overline{x}$ ，$\overline{y}>0$ 使得$x\left(t\right)\le \overline{x}$ 和$y\left(t\right)\le \overline{y}$ ，其中$\overline{x}$ 与$\overline{y}$ 分别表示$x\left(t\right)$ 、$y\left(t\right)$ 的随机上界。

证明：定义如下$C^2$ 函数$V:{ℝ}_{+}^2\to {ℝ}_{+}$ ：

$V\left(t\right)=x\left(t\right)+\frac{1}{\lambda }y\left(t\right).$

根据系统(1.1)，$V\left(t\right)$ 导数为

$\begin{array}{c}\text{d}V\left(x,y\right)=\left[rx\frac{1+ax}{1+ax+by}-c{x}^2-m_{1}x-\frac{1}{\lambda }\left(\frac{qey}{c_{1}e+c_{2}y}+m_{2}y\right)\right]\text{d}t+{\sigma }_{1}x\text{d}{B}_1\left(t\right)+\frac{{\sigma }_2}{\lambda }y\text{d}{B}_2\left(t\right)\\ \le \left[\left(r+m_2\right)x-b{x}^2-m_{2}x-\frac{m_2}{\lambda }y\right]\text{d}t+{\sigma }_{1}x\text{d}{B}_1\left(t\right)+\frac{{\sigma }_2}{\lambda }y\text{d}{B}_2\left(t\right)\\ \le \left(\frac{{\left(r+m_2\right)}^2}{4b}-m_{2}V\right)\text{d}t+{\sigma }_{1}x\text{d}{B}_1\left(t\right)+\frac{{\sigma }_2}{\lambda }y\text{d}{B}_2\left(t\right).\end{array}$

设$U\left(t\right)$ 为如下随机微分方程的解：

$\{\begin{array}{l}\text{d}U\left(t\right)=\left(\frac{{\left(r+m_2\right)}^2}{4b}-m_{2}V\right)\text{d}t+{\sigma }_{1}x\text{d}{B}_1\left(t\right)+\frac{{\sigma }_2}{\lambda }y\text{d}{B}_2\left(t\right),\\ U\left(0\right)=V\left(0\right).\end{array}$

显然，该方程的解$U\left(t\right)$ 有如下形式：

$U\left(t\right)=\frac{{\left(r+m_2\right)}^2}{4b{m}_2}+\left(U\left(0\right)-\frac{{\left(r+m_2\right)}^2}{4b{m}_2}\right)\exp \left(-m_{2}t\right)+M\left(t\right),$

其中

$M\left(t\right)={\sigma }_1{\displaystyle {\int }_0^t\exp \left[-m_2\left(t-s\right)\right]x\left(s\right)\text{d}{B}_1\left(s\right)}+\frac{{\sigma }_2}{\lambda }{\displaystyle {\int }_0^t\exp \left[-m_2\left(t-s\right)\right]y\left(s\right)\text{d}{B}_2\left(s\right)}$

是一个实值连续局部鞅，且满足$M\left(t\right)=0$ (几乎处处)。根据随机比较定理，可得$V\left(t\right)\le U\left(t\right)$ (几乎处处)。后续证明过程与文献[26]的引理2.3及文献[27]的定理3.2类似，因此我们有

$\underset{t\to +\infty }{\lim }\sup \left[x\left(t\right)+\frac{1}{\lambda }y\left(t\right)\right]<+\infty ,a.s.$

此外，由系统(1.1)的第一个方程，可得

$\text{d}x\le x\left(r-m_1-cx\right)\text{d}t+{\sigma }_{1}x\text{d}{B}_1\left(t\right).$

类似于文献[28]中的定理3.2的分析，我们有

$\underset{t\to +\infty }{\lim }\sup E\left[x\right]\le \frac{r-m_1}{c},a.s.$

即对任意的${\epsilon }_1>0$ ，总存在${\delta }_1>0$ 使得

$\underset{t\to +\infty }{\lim }\sup \mathbb{P}\left\{x>{\epsilon }_1\right\}\le {\delta }_1,a.s.$

其中${\delta }_1=\frac{r-m_1}{c}$ 。因此，对任意的${\epsilon }_1>0$ ，存在${\delta }_1>0$ 使得$\mathbb{P}\left\{x>{\epsilon }_1\right\}\le {\delta }_1$ 。即$x$ 是随机最终有界的，存在$\overline{x}>0$ 使得对所有的$t\in \left[0,t_e\right)$ ，有

$x\left(t\right)\le \overline{x},a.s.$

类似的，我们可以证得$y$ 是随机最终有界的，存在$\overline{y}>0$ 使得对所有的$t\in \left[0,t_e\right)$ ，有

$y\left(t\right)\le \overline{y},a.s.$

证毕。

注1 随机最终有界性表明，系统对剧烈且持续性的种群波动具有内在恢复力，从而保障食物网的稳定性与长期行为的可预测性。也就是说，即便面临环境变化或随机扰动，捕食者与食饵种群仍能得以生存。

注2 根据该定理，很容易得到对任意给定的初始值$\left(x\left(0\right),y\left(0\right)\right)\in {ℝ}_{+}^2$ ，系统(1.1)的解$\left(x\left(t\right),y\left(t\right)\right)$ 满足以下性质：

$\underset{t\to +\infty }{\text{lim}}\sup \frac{\ln x\left(t\right)}{t}\le 0,\underset{t\to +\infty }{\text{lim}}\sup \frac{\ln y\left(t\right)}{t}\le 0,a.s.$

3.2.3. 系统(1.1)种群灭绝与持久性

种群灭绝与持久性是兼具理论与现实意义的核心动力学行为。因此，厘清系统(1.1)中捕食者–猎物种群的长期存续趋势(灭绝或持久)至关重要。下文将给出两类种群均值灭绝与均值持久的判定依据。

定理3.3 对任意给定的初始值$\left(x\left(0\right),y\left(0\right)\right)\in {ℝ}_{+}^2$ ，$\left(x\left(t\right),y\left(t\right)\right)$ 是系统(1.1)的解，则我们有如下结果：

(i) 若$r-m_1-\frac{{\sigma }_1^2}{2}<0$ ，则食饵种群是灭绝的；

(ii) 若$r-m_1-\frac{{\sigma }_1^2}{2}=0$ ，则食饵种群是均值非持久的；

(iii) 若$\frac{r}{1+a\overline{x}+b\overline{y}}-\frac{\alpha \overline{y}}{\beta }-m_1-\frac{1}{2}{\sigma }_1^2>0$ ，则食饵种群是持久的。

证明：应用伊藤公式到系统(1.1)的第一个方程，可得

$\text{d}\ln x=\left(r\frac{1+ax}{1+ax+by}-m_1-cx-\frac{\alpha y}{\beta +x^2}-\frac{1}{2}{\sigma }_1^2\right)\text{d}t+{\sigma }_1\text{d}{B}_1\left(t\right)$ (1.9)

上式两边同时在$\left[0,t\right]$ 积分，并除以$t$ 可得

$\frac{1}{t}\frac{\ln x\left(t\right)}{\ln x\left(0\right)}\le r-m_1-\frac{1}{2}{\sigma }_1^2-c{\displaystyle {\int }_0^{t}x\left(s\right)\text{d}s}+\frac{M_1\left(t\right)}{t}$ (1.10)

其中，$M_1\left(t\right)={\displaystyle {\int }_0^t{\sigma }_1\left(s\right)\text{d}{B}_1\left(s\right)}$ 。根据局部鞅的大数定理[22]，可得

$\underset{t\to \infty }{\lim }\frac{M_1\left(t\right)}{t}=0,a.s.$

根据极限存在的性质，对任意的${\epsilon }_2>0$ ，存在常数$T_2>0$ 使得只要$t\ge T_2$ ，有

$\frac{M_1\left(t\right)}{t}\le {\epsilon }_2$ (1.11)

(i) 由(1.10)与(1.11)，根据文献[29]中引理4，可得

$\underset{t\to \infty }{\lim \sup }\frac{1}{t}{\displaystyle {\int }_0^{t}x\left(s\right)\text{d}s}\le \frac{r-m_1-\frac{1}{2}{\sigma }_1^2+{\epsilon }_2}{c},a.s.$ (1.12)

显然，$\underset{t\to \infty }{\lim }x\left(t\right)=0,a.s.$ 如果$r-m_1-\frac{{\sigma }_1^2}{2}<0$ ，即食饵种群是灭绝的。

(ii) 由条件$r-m_1-\frac{{\sigma }_1^2}{2}=0$ 及(1.12)可得

$\underset{t\to \infty }{\lim \sup }\frac{1}{t}{\displaystyle {\int }_0^{t}x\left(s\right)\text{d}s}\le \frac{{\epsilon }_2}{c},a.s.$

由${\epsilon }_2$ 的任意性与$x$ 的非负性，则

$\underset{t\to \infty }{\lim \sup }\frac{1}{t}{\displaystyle {\int }_0^{t}x\left(s\right)\text{d}s}=0,a.s.$

即食饵种群是均值非持久的。

(iii) 由(1.9)及定理3.2可得

$\text{d}\ln x\ge \left(\frac{r}{1+a\overline{x}+b\overline{y}}-cx-\frac{\alpha \overline{y}}{\beta }-m_1-\frac{1}{2}{\sigma }_1^2\right)\text{d}t+{\sigma }_1\text{d}{B}_1\left(t\right)$

因此

$\frac{\ln x\left(t\right)}{\ln x\left(0\right)}\ge \left(\frac{r}{1+a\overline{x}+b\overline{y}}-\frac{\alpha \overline{y}}{\beta }-m_1-\frac{1}{2}{\sigma }_1^2\right)\text{d}t-c{\displaystyle {\int }_0^{t}x\left(s\right)\text{d}s}+{\sigma }_1\text{d}{B}_1\left(t\right)$

显然，如果条件$\frac{r}{1+a\overline{x}+b\overline{y}}-\frac{\alpha \overline{y}}{\beta }-m_1-\frac{1}{2}{\sigma }_1^2>0$ ，根据文献[29]中引理4，可得

$\underset{t\to \infty }{\lim \inf }\frac{1}{t}{\displaystyle {\int }_0^{t}x\left(s\right)\text{d}s}\ge \frac{1}{c}\left(\frac{r}{1+a\overline{x}+b\overline{y}}-\frac{\alpha \overline{y}}{\beta }-m_1-\frac{1}{2}{\sigma }_1^2\right)>0.$

即食饵种群是持久的。

证毕。

定理3.4 对任意给定的初始值$\left(x\left(0\right),y\left(0\right)\right)\in {ℝ}_{+}^2$ ，$\left(x\left(t\right),y\left(t\right)\right)$ 是系统(1.1)的解，则我们有如下结果：

(I) 若$r-m_1-\frac{{\sigma }_1^2}{2}\le 0$ ，则捕食者种群是灭绝的；

(II) 若$r-m_1-\frac{{\sigma }_1^2}{2}>0$ ，则

(i) 若$A=\frac{\lambda \alpha \left(r-m_1-\frac{1}{2}{\sigma }_1^2\right)}{c\beta }-m_2-\frac{1}{2}{\sigma }_2^2>0$ ，则捕食者种群是灭绝的；

(ii) 若$B=\frac{\lambda \alpha }{\beta +{\overline{x}}^2}\left(\frac{r}{1+a\overline{x}+b\overline{y}}-m_1-\frac{1}{2}{\sigma }_1^2\right)-c\left(\frac{q}{c_1}+m_2+\frac{1}{2}{\sigma }_2^2 \right)>0$ ，则捕食者种群是持久的。

证明：(I) 由定理3.3结论可知，

$\underset{t\to \infty }{\lim \sup }\frac{1}{t}{\displaystyle {\int }_0^{t}x\left(s\right)\text{d}s}\le 0,a.s.$

应用伊藤公式到系统(1.1)的第二个方程，可得

$\text{d}\ln y=\left(\frac{\lambda \alpha x}{\beta +x^2}-\frac{qe}{c_{1}e+c_{2}y}-m_2-\frac{1}{2}{\sigma }_2^2\right)\text{d}t+{\sigma }_2\text{d}{B}_2\left(t\right)$ (1.13)

上式两边同时在$\left[0,t\right]$ 积分，并除以$t$ 可得

$\frac{1}{t}\frac{\ln y\left(t\right)}{\ln y\left(0\right)}\le -m_2-\frac{1}{2}{\sigma }_2^2+\frac{\lambda \alpha }{\beta }{\displaystyle {\int }_0^{t}x\left(s\right)\text{d}s}+\frac{M_2\left(t\right)}{t}$ (1.14)

其中，$M_2\left(t\right)={\displaystyle {\int }_0^t{\sigma }_2\left(s\right)\text{d}{B}_2\left(s\right)}$ 。根据局部鞅的大数定理[22]，可得

$\underset{t\to \infty }{\lim }\frac{M_2\left(t\right)}{t}=0a.s.$ (1.15)

因此

$\underset{t\to \infty }{\lim \sup }\frac{1}{t}{\displaystyle {\int }_0^{t}y\left(s\right)\text{d}s}\le -m_2-\frac{1}{2}{\sigma }_2^2<0$

即捕食者种群是灭绝的。

(II) 现在我们考虑情况$r-m_1-\frac{{\sigma }_1^2}{2}>0$ 。

(i) 由(1.10)及局部鞅的大数定理[22]，可得

$\underset{t\to \infty }{\lim \sup }\frac{1}{t}{\displaystyle {\int }_0^{t}x\left(s\right)\text{d}s}\le \frac{r-m_1-\frac{1}{2}{\sigma }_1^2}{c},a.s.$

结合(1.14)与(1.15)，可得

$\underset{t\to \infty }{\lim \sup }\frac{1}{t}{\displaystyle {\int }_0^{t}y\left(s\right)\text{d}s}\le \frac{\lambda \alpha \left(r-m_1-\frac{1}{2}{\sigma }_1^2\right)}{c\beta }-m_2-\frac{1}{2}{\sigma }_2^2\triangleq \text{Α}$

显然，如果$\text{Α}<0$ ，则捕食者种群是灭绝的。

(ii) 由(1.9)与(1.13)可得

$\frac{1}{t}\frac{\ln x\left(t\right)}{\ln x\left(0\right)}\ge \frac{r}{1+a\overline{x}+b\overline{y}}-m_1-\frac{1}{2}{\sigma }_1^2-c\frac{1}{t}{\displaystyle {\int }_0^{t}x\left(s\right)\text{d}s}-\frac{\alpha }{\beta }\frac{1}{t}{\displaystyle {\int }_0^{t}y\left(s\right)\text{d}s}+\frac{M_1\left(t\right)}{t}$ (1.16)

与

$\frac{1}{t}\frac{\ln y\left(t\right)}{\ln y\left(0\right)}\ge \frac{\lambda \alpha }{\beta +{\overline{x}}^2}\frac{1}{t}{\displaystyle {\int }_0^{t}x\left(s\right)}-\frac{q}{c_1}-m_2-\frac{1}{2}{\sigma }_2^2+\frac{M_2\left(t\right)}{t}$ (1.17)

则$\frac{\lambda \alpha }{\beta +{\overline{x}}^2}\times$ (1.16) $+c\times$ (1.17)可得

$\begin{array}{c}\frac{1}{t}\frac{\ln x\left(t\right)}{\ln x\left(0\right)}\times \frac{\lambda \alpha }{\beta +{\overline{x}}^2}+\frac{1}{t}\frac{\ln y\left(t\right)}{\ln y\left(0\right)}\times c\ge \left[\frac{\lambda \alpha }{\beta +{\overline{x}}^2}\left(\frac{r}{1+a\overline{x}+b\overline{y}}-m_1-\frac{1}{2}{\sigma }_1^2\right)-c\left(\frac{q}{c_1}+m_2+\frac{1}{2}{\sigma }_2^2\right)\right]\\ -\frac{\lambda {\alpha }^2}{\beta \left(\beta +{\overline{x}}^2\right)}\frac{1}{t}{\displaystyle {\int }_0^{t}y\left(s\right)\text{d}s}+\frac{\lambda \alpha }{\beta +{\overline{x}}^2}\frac{M_1\left(t\right)}{t}+c\frac{M_2\left(t\right)}{t}\end{array}$

由注2可得

$\begin{array}{c}\frac{1}{t}\frac{\ln y\left(t\right)}{\ln y\left(0\right)}\ge \frac{1}{c}\left[\frac{\lambda \alpha }{\beta +{\overline{x}}^2}\left(\frac{r}{1+a\overline{x}+b\overline{y}}-m_1-\frac{1}{2}{\sigma }_1^2\right)-c\left(\frac{q}{c_1}+m_2+\frac{1}{2}{\sigma }_2^2\right)\right]\\ -\frac{\lambda {\alpha }^2}{\beta \left(\beta +{\overline{x}}^2\right)}\frac{1}{t}{\displaystyle {\int }_0^{t}y\left(s\right)\text{d}s}+\frac{\lambda \alpha }{c\left(\beta +{\overline{x}}^2\right)}\frac{M_1\left(t\right)}{t}+\frac{M_2\left(t\right)}{t}\end{array}$

根据文献[29]中引理4，可得

$\underset{t\to \infty }{\lim \inf }\frac{1}{t}{\displaystyle {\int }_0^{t}y\left(s\right)\text{d}s}\ge \frac{\beta \left(\beta +{\overline{x}}^2\right)}{c\lambda {\alpha }^2}\left[\frac{\lambda \alpha }{\beta +{\overline{x}}^2}\left(\frac{r}{1+a\overline{x}+b\overline{y}}-m_1-\frac{1}{2}{\sigma }_1^2\right)-c\left(\frac{q}{c_1}+m_2+\frac{1}{2}{\sigma }_2^2\right)\right]$

显然，如果条件$B>0$ ，则有

$\underset{t\to \infty }{\lim \inf }\frac{1}{t}{\displaystyle {\int }_0^{t}y\left(s\right)\text{d}s}>0$

即捕食者种群是持久的。

证毕。

注3 定理3.3与定理3.4表明，系统(1.1)中食饵种群灭绝情况下，捕食者种群一定是灭绝的。但当食饵种群存在时，捕食者种群的生存与否取决于环境随机扰动强度。这一结论进一步说明了环境随机扰动对种群生态系统稳定性的重要影响。

3.2.4. 系统(1.1)稳态的存在性

遍历性是随机微分方程所建种群动力学模型的一项基本特征。具体而言，遍历性意味着系统存在唯一的平稳分布，这一性质为种群长期的周期性动态变化提供了生态学解释。本节将分析系统(1.1)的长期统计行为，并推导其平稳分布存在且唯一的充分条件。

定理3.5 若条件

$\Pi =\frac{\lambda \alpha }{\beta +{\overline{x}}^2}\left(\frac{r}{1+a\overline{x}+b\overline{y}}-m_1-\frac{{\sigma }_1^2}{2}\right)-c\left(\frac{q}{c_1}+m_2+\frac{{\sigma }_2^2}{2}\right)>0$

成立，则对任意给定的初始值$\left(x\left(0\right),y\left(0\right)\right)\in {ℝ}_{+}^2$ ，系统(3.1)具有唯一的遍历性稳态分布。

证明：要证明该定理成立，需逐个验证引理1中的两个条件成立。经简单计算可知，系统(1.1)的扩散矩阵$D=\text{diag}\left({\sigma }_1^2{x}^2,{\sigma }_2^2{y}^2\right)$ 在${ℝ}_{+}^2$ 的任意紧子集上都是正定的，这表明引理1的条件(i)是成立的。

现在，我们证明引理1的条件(ii)是成立的。定义如下$C^2$ 函数$f:{ℝ}_{+}^2\to ℝ$ ：

$f\left(x,y\right)=\frac{1}{p+1}{\left(x+\frac{1}{\lambda }y\right)}^{p+1}-u\left(\frac{e_1}{m+\alpha \eta A+\overline{P}}\ln x+c\ln y-\frac{c{e}_1}{d{\left(m+\alpha \eta A+\overline{P}\right)}^2}y\right)-x-y.$

其中，$p$ 是正常数，满足约束$0<p<\frac{m_2}{2^{p+2}\left[{\sigma }_1^2\vee {\sigma }_2^2\right]}$ ，而$u>0$ 的取值将在后续推导中确定。由于函数$f\left(x,y\right)$ 是连续函数，则函数存在唯一驻点$\left(x_{*},y_{*}\right)\in {ℝ}_{+}^2$ 且$f\left(x,y\right)$ 在该点取最小值。因此，构造如下非负$C^2$ 函数$V:{ℝ}_{+}^2\to ℝ$ ：

$\begin{array}{c}V\left(x,y\right)=f\left(x,y\right)-f\left(x_{*},y_{*}\right)\\ =\frac{1}{p+1}{\left(x+\frac{1}{\lambda }y\right)}^{p+1}-u\left(\frac{\lambda \alpha }{\beta +{\overline{x}}^2}\ln x+c\ln y-\frac{\lambda {\alpha }^2}{m_2{\left(\beta +{\overline{x}}^2\right)}^2}y\right)\\ -x-y-f\left(P^{*},Z^{*}\right)\\ =V_1\left(x,y\right)+V_2\left(x,y\right)+V_3\left(x,y\right)-f\left(x_{*},y_{*}\right).\end{array}$

应用伊藤公式到$V_1\left(x,y\right)$ ，可得

$\begin{array}{c}ℒV_1\left(x,y\right)\le {\left(x+\frac{1}{\lambda }y\right)}^p\left[\left(r-m_1\right)x-c{x}^2-\frac{1}{\lambda }{m}_{2}Z\right]+\frac{p}{2}{\left(x+\frac{1}{\lambda }y\right)}^{p-1}\left[{\sigma }_1^2{x}^2+{\sigma }_2^2{\left(\frac{1}{\lambda }\right)}^2{y}^2\right]\\ \le \left(r-m_1\right)x{\left(x+\frac{1}{\lambda }y\right)}^p-c{x}^{p+2}-{\left(\frac{1}{\lambda }\right)}^{p+1}{m}_2{y}^{p+1}+\frac{p}{2}{\left(x+\frac{1}{\lambda }y\right)}^{p+1}\left[{\sigma }_1^2\vee {\sigma }_2^2\right].\end{array}$

由于${\left|{\displaystyle {\sum }_{i=1}^n{q}_i}\right|}^n\le k^{n-1}{\displaystyle {\sum }_{i=1}^n{\left|q_i\right|}^n}$ ，则

$\begin{array}{c}ℒV_1\left(x,y\right)\le -\frac{c}{2}{x}^{p+2}-\frac{1}{2}{m}_2{\left(\frac{1}{\lambda }\right)}^{p+1}{y}^{p+1}-\frac{1}{4}{m}_2{\left(\frac{1}{\lambda }\right)}^{p+1}{y}^{p+1}\\ -\left(\frac{1}{4}{m}_2-p{2}^p\left[{\sigma }_1^2\vee {\sigma }_2^2\right]\right){\left(\frac{1}{\lambda }\right)}^{p+1}{y}^{p+1}+g\\ \le -\frac{c}{2}{x}^{p+2}-\frac{1}{2}{m}_2{\left(\frac{1}{\lambda }\right)}^{p+1}{y}^{p+1}+g,\end{array}$

其中，

$g={\text{sup}}_{\left(x,y\right)\in {ℝ}_{+}^2}\left\{-\frac{c}{2}{x}^{p+2}-\frac{1}{4}{m}_2{\left(\frac{1}{\lambda }\right)}^{p+1}{y}^{p+1}+\left(r-m_1\right)x{\left(x+\frac{1}{\lambda }y\right)}^p+p{2}^p{x}^{p+1}\left[{\sigma }_1^2\vee {\sigma }_2^2\right]\right\}<\infty .$

类似的，我们可以得到

$\begin{array}{c}ℒV_2\left(x,y\right)=-u\left[\frac{\lambda \alpha }{\beta +{\overline{x}}^2}\ln x+c\ln y-\frac{\lambda {\alpha }^2}{m_2{\left(\beta +{\overline{x}}^2\right)}^2}y\right]\\ =-u[\frac{\lambda \alpha }{\beta +{\overline{x}}^2}\left(r\frac{1+ax}{1+ax+by}-cx-\frac{\alpha y}{\beta +x^2}-m_1-\frac{{\sigma }_1^2}{2}\right)\\ +c\left(\frac{\lambda \alpha x}{\beta +x^2}-\frac{qe}{c_{1}e+c_{2}y}-m_2-\frac{{\sigma }_2^2}{2}\right)\\ -\frac{\lambda {\alpha }^2}{m_2{\left(\beta +{\overline{x}}^2\right)}^2}\left(\frac{\lambda \alpha xy}{\beta +x^2}-\frac{qey}{c_{1}e+c_{2}y}-m_{2}y\right)]\\ \le -u\Pi +\frac{u{\lambda }^2{\alpha }^3}{m_2{\beta }^3}xy+\frac{u\lambda q{\alpha }^2}{m_2{c}_1{\beta }^2}y.\end{array}$

和

$\begin{array}{c}ℒV_3\left(x,y\right)=-rx\frac{1+ax}{1+ax+by}+c{x}^2+\frac{\alpha xy}{\beta +x^2}+m_{1}x-\frac{\lambda \alpha xy}{\beta +x^2}+\frac{qey}{c_{1}e+c_{2}y}+m_{2}y\\ \le c{x}^2+m_{1}x+\frac{\alpha }{\beta }xy+\left(\frac{q}{c_1}+m_2\right)y.\end{array}$

因此

$\begin{array}{c}ℒV\left(x,y\right)=ℒV_1\left(x,y\right)+ℒV_2\left(x,y\right)+ℒV_3\left(x,y\right)\\ \le -u\Pi +c{x}^2+m_{1}x+\left(\frac{u{\lambda }^2{\alpha }^3}{m_2{\beta }^3}+\frac{\alpha }{\beta }\right)xy+\left(\frac{q}{c_1}+m_2+\frac{u\lambda q{\alpha }^2}{m_2{c}_1{\beta }^2}\right)y\\ -\frac{c}{2}{x}^{p+2}-\frac{1}{2}{m}_2{\left(\frac{1}{\lambda }\right)}^{p+1}{y}^{p+1}+g.\end{array}$

考虑如下紧子集$E$ ：

$E=\left\{\left(x,y\right)\in {ℝ}_{+}^2:\epsilon <x<\frac{1}{\epsilon },\epsilon <y<\frac{1}{\epsilon }\right\},$

其中，$\epsilon$ 为充分小的正常数，并满足如下条件：

$-\Pi u+\left(\frac{c}{m}+\frac{cu{e}_1^2}{d{m}^3}\right)\epsilon +b{\epsilon }^2+S_1\le -1,$ (1.18)

$-\Pi u+S_2-\frac{b}{4}{\epsilon }^{-p-2}\le -1,$ (1.19)

$-\Pi u+\left(\frac{cu{e}_1{e}_{2}A}{d{m}^3}+d+\theta +\frac{c}{m}+\frac{cu{e}_1^2}{d{m}^3}\right)\epsilon +S_3\le -1,$ (1.20)

与

$-\Pi u+S_4-\frac{1}{4}\left(d-\frac{e_2}{\alpha \eta }\right){\left(\frac{c}{e_1}\right)}^p{\epsilon }^{-p-1}\le -1,$ (1.21)

其中，

$S_1={\text{sup}}_{\left(x,y\right)\in {ℝ}_{+}^2}\left\{\left(\frac{q}{c_1}+m_2+\frac{u\lambda q{\alpha }^2}{m_2{c}_1{\beta }^2}\right)y+\left(\frac{u{\lambda }^2{\alpha }^3}{m_2{\beta }^3}+\frac{\alpha }{\beta }\right)y^2+g\right\},$

$\begin{array}{l}{S}_2={\text{sup}}_{\left(x,y\right)\in {ℝ}_{+}^2}\{\left[c+\frac{1}{2}\left(\frac{u{\lambda }^2{\alpha }^3}{m_2{\beta }^3}+\frac{\alpha }{\beta }\right)\right]x^2+m_{1}x+\left(\frac{q}{c_1}+m_2+\frac{u\lambda q{\alpha }^2}{m_2{c}_1{\beta }^2}\right)y\\ +\frac{1}{2}\left(\frac{u{\lambda }^2{\alpha }^3}{m_2{\beta }^3}+\frac{\alpha }{\beta }\right)y^2-\frac{c}{4}{x}^{p+2}-\frac{1}{2}{m}_2{\left(\frac{1}{\lambda }\right)}^{p+1}{y}^{p+1}+g\},\end{array}$

$S_3={\text{sup}}_{\left(x,y\right)\in {ℝ}_{+}^2}\left\{\left(c+\left(\frac{u{\lambda }^2{\alpha }^3}{m_2{\beta }^3}+\frac{\alpha }{\beta }\right)\epsilon \right)x^2+m_{1}x\right\},$

与

$\begin{array}{l}{S}_4={\text{sup}}_{\left(x,y\right)\in {ℝ}_{+}^2}\{\left[c+\frac{1}{2}\left(\frac{u{\lambda }^2{\alpha }^3}{m_2{\beta }^3}+\frac{\alpha }{\beta }\right)\right]x^2+m_{1}x+\left(\frac{q}{c_1}+m_2+\frac{u\lambda q{\alpha }^2}{m_2{c}_1{\beta }^2}\right)y\\ +\frac{1}{2}\left(\frac{u{\lambda }^2{\alpha }^3}{m_2{\beta }^3}+\frac{\alpha }{\beta }\right)y^2-\frac{c}{2}{x}^{p+2}-\frac{1}{4}{m}_2{\left(\frac{1}{\lambda }\right)}^{p+1}{y}^{p+1}\}.\end{array}$

则

${ℝ}_{+}^2/E=E_1\cup E_2\cup E_3\cup E_4,$

其中，

$E_1=\left\{\left(x,y\right)\in {ℝ}_{+}^2:0<x\le \epsilon \right\},E_2=\left\{\left(x,y\right)\in {ℝ}_{+}^2:x\ge \frac{1}{\epsilon }\right\},$

$E_3=\left\{\left(x,y\right)\in {ℝ}_{+}^2:0<y\le \epsilon \right\},E_4=\left\{\left(P,Z\right)\in {ℝ}_{+}^2:y\ge \frac{1}{\epsilon }\right\}.$

接下来我们证明对任意的$\left(x,y\right)\in {ℝ}_{+}^2/E$ ，$ℒV\left(x,y\right)$ 的负性。

情况(1)：若$\left(x,y\right)\in E_1$ ，则$xy\le \epsilon y\le \epsilon \left(1+y^2\right)$ ，从而

$\begin{array}{c}ℒV\left(x,y\right)\le -u\Pi +c{x}^2+m_{1}x+\left(\frac{u{\lambda }^2{\alpha }^3}{m_2{\beta }^3}+\frac{\alpha }{\beta }\right)xy+\left(\frac{q}{c_1}+m_2+\frac{u\lambda q{\alpha }^2}{m_2{c}_1{\beta }^2}\right)y\\ -\frac{c}{2}{x}^{p+2}-\frac{1}{2}{m}_2{\left(\frac{1}{\lambda }\right)}^{p+1}{y}^{p+1}+g\\ \le -u\Pi +c{\epsilon }^2+m_1\epsilon +\left(\frac{u{\lambda }^2{\alpha }^3}{m_2{\beta }^3}+\frac{\alpha }{\beta }\right)\epsilon \left(1+y^2\right)+\left(\frac{q}{c_1}+m_2+\frac{u\lambda q{\alpha }^2}{m_2{c}_1{\beta }^2}\right)y\\ -\frac{c}{2}{\epsilon }^{p+2}-\frac{1}{2}{m}_2{\left(\frac{1}{\lambda }\right)}^{p+1}{y}^{p+1}+g\\ \le -u\Pi +\left(\frac{u{\lambda }^2{\alpha }^3}{m_2{\beta }^3}+\frac{\alpha }{\beta }+m_1\right)\epsilon +c{\epsilon }^2+S_1.\end{array}$

根据(1.18)，对任意的$\left(x,y\right)\in E_1$ ，我们可得$ℒV\left(x,y\right)\le -1$ 。

情况(2)：若$\left(x,y\right)\in E_2$ ，有

$\begin{array}{c}ℒV\left(x,y\right)\le -u\Pi +c{x}^2+m_{1}x+\left(\frac{u{\lambda }^2{\alpha }^3}{m_2{\beta }^3}+\frac{\alpha }{\beta }\right)xy+\left(\frac{q}{c_1}+m_2+\frac{u\lambda q{\alpha }^2}{m_2{c}_1{\beta }^2}\right)y\\ -\frac{c}{2}{x}^{p+2}-\frac{1}{2}{m}_2{\left(\frac{1}{\lambda }\right)}^{p+1}{y}^{p+1}+g\\ \le -u\Pi +c{x}^2+m_{1}x+\left(\frac{u{\lambda }^2{\alpha }^3}{m_2{\beta }^3}+\frac{\alpha }{\beta }\right)\frac{x^2+y^2}{2} +\left(\frac{q}{c_1}+m_2+\frac{u\lambda q{\alpha }^2}{m_2{c}_1{\beta }^2}\right)y\\ -\frac{c}{4}{x}^{p+2}-\frac{c}{4}{x}^{p+2}-\frac{1}{2}{m}_2{\left(\frac{1}{\lambda }\right)}^{p+1}{y}^{p+1}+g\\ \le -u\Pi -\frac{c}{4}{\epsilon }^{p+2}+S_2.\end{array}$

根据(1.19)，对任意的$\left(x,y\right)\in E_2$ ，我们可得$ℒV\left(x,y\right)\le -1$ 。

情况(3)：若$\left(x,y\right)\in E_3$ ，则$xy\le x\epsilon \le \epsilon \left(1+x^2\right)$ ，从而

$\begin{array}{c}ℒV\left(x,y\right)\le -u\Pi +c{x}^2+m_{1}x+\left(\frac{u{\lambda }^2{\alpha }^3}{m_2{\beta }^3}+\frac{\alpha }{\beta }\right)xy+\left(\frac{q}{c_1}+m_2+\frac{u\lambda q{\alpha }^2}{m_2{c}_1{\beta }^2}\right)y\\ -\frac{c}{2}{x}^{p+2}-\frac{1}{2}{m}_2{\left(\frac{1}{\lambda }\right)}^{p+1}{y}^{p+1}+g\\ \le -u\Pi +c{x}^2+m_{1}x+\left(\frac{u{\lambda }^2{\alpha }^3}{m_2{\beta }^3}+\frac{\alpha }{\beta }\right)\epsilon \left(1+x^2\right)+\left(\frac{q}{c_1}+m_2+\frac{u\lambda q{\alpha }^2}{m_2{c}_1{\beta }^2}\right)\epsilon \\ -\frac{c}{2}{x}^{p+2}-\frac{1}{2}{m}_2{\left(\frac{1}{\lambda }\right)}^{p+1}{y}^{p+1}+g\\ \le -u\Pi +\left(\frac{q}{c_1}+m_2+\frac{u\lambda q{\alpha }^2}{m_2{c}_1{\beta }^2}+\frac{u{\lambda }^2{\alpha }^3}{m_2{\beta }^3}+\frac{\alpha }{\beta }\right)\epsilon +S_3\end{array}$

根据(1.20)，对任意的$\left(x,y\right)\in E_3$ ，我们可得$ℒV\left(x,y\right)\le -1$ 。

情况(4)：若$\left(x,y\right)\in E_4$ ，有

$\begin{array}{c}ℒV\left(x,y\right)\le -u\Pi +c{x}^2+m_{1}x+\left(\frac{u{\lambda }^2{\alpha }^3}{m_2{\beta }^3}+\frac{\alpha }{\beta }\right)xy+\left(\frac{q}{c_1}+m_2+\frac{u\lambda q{\alpha }^2}{m_2{c}_1{\beta }^2}\right)y\\ -\frac{c}{2}{x}^{p+2}-\frac{1}{2}{m}_2{\left(\frac{1}{\lambda }\right)}^{p+1}{y}^{p+1}+g\\ \le -u\Pi +c{x}^2+m_{1}x+\left(\frac{u{\lambda }^2{\alpha }^3}{m_2{\beta }^3}+\frac{\alpha }{\beta }\right)\frac{x^2+y^2}{2} +\left(\frac{q}{c_1}+m_2+\frac{u\lambda q{\alpha }^2}{m_2{c}_1{\beta }^2}\right)y\\ -\frac{c}{2}{x}^{p+2}-\frac{1}{4}{m}_2{\left(\frac{1}{\lambda }\right)}^{p+1}{y}^{p+1}-\frac{1}{4}{m}_2{\left(\frac{1}{\lambda }\right)}^{p+1}{y}^{p+1}+g\\ \le -u\Pi -\frac{1}{4}{m}_2{\left(\frac{1}{\lambda }\right)}^{p+1}{y}^{p+1}+S_4.\end{array}$

根据(1.21)，对任意的$\left(x,y\right)\in E_4$ ，我们可得$ℒV\left(x,y\right)\le -1$ 。

综上，我们可得对任意的$\left(x,y\right)\in {ℝ}_{+}^2/E$ ，都有$ℒV\left(x,y\right)\le -1$ 成立。即引理1中的条件2是成立的。因此，系统(1.1)存在遍历性稳态分布。

证毕。

4. 数值模拟分析

本小节将通过数值模拟，验证理论结果的有效性，并进一步探究随机系统(1.1)的复杂动力学行为。具体地，采用Milstein高阶方法[30]，探究环境白噪声对系统(1.1)的影响。在后续数值模拟中，除特别说明外，参数取值均如下：

$r=3,\alpha =0.8,\beta =2.5,\lambda =0.7,c=0.6,m_1=0.2,$

$m_2=0.1,q=0.5,e=0.1,c_1=0.6,c_2=0.8,a=0.2,b=0.3,$

并初始值设定为$\left(x\left(0\right),y\left(0\right)\right)=\left(2.5,3.5\right)$ 。

为探究环境白噪声对捕食者–食饵种群动力学行为的影响，本文首先通过控制变量，假设随机系统(1.1)不存在环境波动扰动。换言之，系统(1.1)退化为对应的确定性系统。在此条件下，确定性系统存在唯一的共存平衡点$E^{*}\left(0.5993,3.5191\right)$ ，且该平衡点具有局部渐近稳定性(见图1)，这表明无论初始状态如何，系统最终均会趋于该平衡状态。显然，这个数值模拟结果说明食饵种群和捕食者种群可以形成常稳态生长共存模式。

Figure 1. The dynamic behavior and phase diagram of predator and prey in deterministic system

图1. 确定性系统的捕食者与食饵种群动力学行为与相图

接下来，通过固定环境白噪声强度${\sigma }_2=0.01$ ，令${\sigma }_1$ 在$\left[0,2\right]$ 范围内变动，探究环境白噪声对随机系统(1.1)动力学行为的影响，分析白噪声强度对捕食者–食饵种群生存–灭绝特性的作用，其余参数均与图1保持一致。从图2(a)清晰地可以看到，随着${\sigma }_1$ 取值的增加，在I，II与III三个不同区域内，系统(1.1)中食饵种群的动力行为由平均持续生存转变为灭绝。事实上，根据定理3.3，当白噪声强度${\sigma }_1$ 达到临界阈值$\left({\sigma }_1\approx 2.3664\right)$ 时，系统(3.1)中食饵种群将趋于灭绝，这一理论结论与本文的数值分析结果相吻合。图2(b)展示了与图2(a)中时间区间$t\in \left[0,1000\right]$ 相对应，${\sigma }_1$ 分别取0.1，2.2，2.8三个不同数值时，系统(3.1)中食饵种群的动力学变化。由图2可知，固定环境白噪声强度取值比较小时，食饵种群密度在常稳态附件出现振荡，食饵种群密度整体幅度变化不大。但是当固定环境白噪声强度取值适中时，食饵种群密度峰度大幅度增大，这说明适度的环境白噪声有利于食饵种群的持久生存。然而，当固定环境白噪声强度取值超过关键值时，会导致食饵种群快速灭绝。

Figure 2. When white noise ${\sigma }_2=0.01$ , the impact of environmental white noise ${\sigma }_1\left(0\le {\sigma }_1\le 2\right)$ on the dynamics of the system (1.1)

图2. 固定参数${\sigma }_2=0.01$ 时，环境白噪声${\sigma }_1\left(0\le {\sigma }_1\le 2\right)$ 对系统(1.1)动力学的影响

取${\sigma }_1=0.01$ 时，不难验证定理3.3与定理3.4成立，因此系统(1.1)满足均值持久性，且存在唯一的平稳分布(见图3)。值得注意的是，从图3(a)和图3(b)中可以观察到，系统(1.1)的种群密度在其对应确定性系统稳态的小邻域内波动，这表明强度较弱的环境白噪声，长期来看无法改变捕食者与食饵种群的生存状态，会以均值稳态的生长共存模式存在。此外，捕食者与食饵的种群密度分别在区间(0.1, 1.6)和(2.7, 4.5)内呈正态分布，其中红色曲线代表捕食者与食饵种群的概率密度函数，矩形的宽度表示不同组距，矩形的高度则对应组频。有趣的是，捕食者与食饵种群密度分布的频率峰值分别出现在0.5993和3.5191处(见图3(c)、3(d))，这两个数值几乎以确定性系统的稳态平衡点为中心。换言之，在此条件下系统的稳定性表现更为突出，生长共存模式是以平衡点为中心的均值稳态。当系统处于中等噪声强度${\sigma }_1=0.1$ 时，通过对比图3(c)、图3(d)与图4(c)、图4(d)可以发现，食饵种群密度的分布呈现显著的右偏特征，分布范围更宽且峰值更低；而捕食者种群密度的分布形态基本保持不变，但同样出现分布范围拓宽、峰值降低的现象。这意味着白噪声强度的增加，会使捕食者种群的密度分布在高位区间内更为分散。此外，我们从图3(a)、图3(b)与图4(a)、图4(b)可以发现，围绕确定性稳态波动的捕食者与食饵种群密度波动幅度，会随白噪声强度的增大而加剧，表明环境白噪声会放大捕食者与食饵种群密度的随机振荡幅度，即白噪声强度越高，系统越难以维持稳定状态，这个研究结果说明随着白噪声强度的增加，食饵种群与捕食者种群所成的生态稳定性越来越弱。然而，进一步将白噪声强度增大至${\sigma }_1=2.8$ 时，易验证$r-m_1-\frac{{\sigma }_1^2}{2}=-2.24<0$ 。根据定理3.3(i)的结论，食饵种群将趋于灭绝，而捕食者种群由于专一依赖食饵作为唯一食物来源，在食饵消亡后失去能量供给，亦无法避免灭绝的命运。即便确定性系统的正平衡点是渐近稳定的，这一结果也可能破坏生态系统的平衡(见图5)。从生物学角度出发，高强度的环境波动作用于食饵种群死亡率，会直接导致其灭绝；而捕食者会因食物资源匮乏相继走向灭绝，这一过程与环境波动对捕食者种群的直接作用无关。值得指出的是，在确定性系统中，食饵与捕食者种群并不会同时灭绝(即平衡点(0,0)是不稳定的)，但随机系统却实现了这一结果，即环境白噪声让确定性系统中原本不稳定的零平衡点变为稳定状态，这说明高强度的环境白噪声对捕食生态系统稳定性的影响作用较大，完全导致生态稳定性消失，进而诱发食饵种群与捕食者种群趋于灭绝。

Figure 3. (a), (b): When white noise ${\sigma }_1=0.01,{\sigma }_2=0.01$ , the evolution of predator and prey population density for the system (1.1) and its corresponding deterministic system; (c), (d): Distribution of predator and prey population density in the system (1.1)

图3. (a)，(b)：白噪声强度为${\sigma }_1=0.01,{\sigma }_2=0.01$ 时，系统(1.1)及其对应确定性系统的捕食者与食饵种群密度随时间演化；(c)，(d)：系统(1.1)中的捕食者与食饵种群密度分布

Figure 4. (a), (b): When white noise ${\sigma }_1=0.1,{\sigma }_2=0.01$ , the evolution of predator and prey population density for the system (1.1) and its corresponding deterministic system; (c), (d): Distribution of predator and prey population density in the system (1.1)

图4. (a)，(b)：白噪声强度为${\sigma }_1=0.1,{\sigma }_2=0.01$ 时，系统(1.1)及其对应确定性系统的捕食者与食饵种群密度随时间演化；(c)，(d)：系统(1.1)中的捕食者与食饵种群密度分布

Figure 5. (a), (b): When white noise ${\sigma }_1=2.8,{\sigma }_2=0.01$ , the evolution of predator and prey population density for the system (1.1) and its corresponding deterministic system

图5. (a)，(b)：白噪声强度为${\sigma }_1=2.8,{\sigma }_2=0.01$ 时，系统(1.1)及其对应确定性系统的捕食者与食饵种群密度随时间演化

由此可见，捕食者–食饵种群的生存状态高度依赖于环境白噪声的强度，过高的噪声强度不利于捕食者–食饵种群的持久生存，这也意味着在现实场景中，合理调控环境白噪声的强度，有望实现对捕食者与食饵种群生物量的管控。需要强调的是，图2、图3、图4及图5的结果均证实：相较于对应的确定性系统，纳入环境白噪声的捕食者–食饵种群随机系统，能更精准地模拟和预测捕食者与食饵种群的生长动力学行为。研究白噪声强度${\sigma }_2$ 对随机系统(1.1)动力学的影响时，可观察到相似的作用规律，因此不再赘述。

5. 结论

捕食者与食饵种群间的相互作用是一种普遍存在的生态现象，在真实的种群生态系统中，这一作用过程往往受到诸多外在因素的共同影响。环境随机波动作为一类重要的外在因素，对种群动力学过程，尤其对种群的生存与灭绝，具有至关重要的影响[8] [9] [13]-[15] [31]-[33]。为更深入地阐释捕食者与食饵种群间的相互作用机制，本文考虑外部环境噪声的影响，构建了一类带有环境随机扰动的捕食者–食饵种群生态系统，并通过理论分析与数值模拟相结合的方式，探究了环境随机扰动对捕食者–食饵种群动力学的影响机制。本文首先借助李雅普诺夫方法，证明了系统(1.1)存在唯一的全局正解；并基于随机比较原理，验证了该解以概率1保持随机最终有界性。在此基础上，进一步推导了仅捕食者种群灭绝、食饵与捕食者种群均走向灭绝的充分条件，同时也得到了食饵与捕食者种群实现共存的判别条件。尤为关键的是，当系统满足共存条件时，通过构造合适的李雅普诺夫函数，本文证明了该系统存在唯一的遍历平稳分布。数值分析结果验证了本文的理论结论，同时进一步表明，系统的动力学行为与环境白噪声密切相关。

环境随机性让系统的设定更贴合实际，而白噪声的存在不利于捕食者与食饵种群的生存。如图2所示，白噪声强度可显著改变食饵种群的平均持续性与灭绝状态，这一结果与本文的理论分析结论相一致。事实上，当白噪声强度较弱时，其并不会改变捕食者与食饵种群的生存状态，这表明该生态系统具备一定的自我调节能力。值得注意的是，随着环境白噪声强度的增大，系统的随机特性愈发显著，表现出不同程度的随机波动(见图3、图4)。此外，若进一步增大白噪声强度，捕食者与食饵种群密度的振荡幅度会大幅加剧，进而导致捕食者与食饵种群均快速趋于灭绝(见图5)。出现这一结果的原因在于，在此强度的环境白噪声扰动下，捕食者与食饵种群均难以适应环境噪声的干扰。因此，环境白噪声强度可作为潜在的调控机制，影响捕食者与食饵种群的生存与灭绝。

综上，环境随机扰动能够调控捕食者–食饵系统的动力学行为，进而在维持捕食者与食饵共存、避免种群发生灾难性灭绝、维系生态系统稳定性及生物多样性方面发挥重要作用。然而，仍有一些颇具研究价值的问题值得进一步探讨。可将更多影响种群动力学的现实因素纳入本文所构建的模型，例如食饵避难效应[34] [35]、空间扩散效应[20] [36]、脉冲扰动[8] [13] [37]、时滞效应[38] [39]等，并探究这些因素如何影响复杂模型的动力学特征。本文将以上方向列为未来的研究重点。

基金项目

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

NOTES

^*通讯作者。