污染水体中捕食种群具有常数投放率的脉冲状态反馈控制模型

Impulsive State Feedback Control Model of Predator Population with Constant Release Rate in Polluted Water

摘要: 水体富营养化引起的破坏水生生态系统的现象时有发生，本文研究了污染水体中捕食种群具有常数投放率和脉冲状态反馈控制的食饵捕食系统。通过构造李雅普诺夫函数得到了系统在无脉冲效应的情况下全局渐近稳定的充分条件。利用半连续动力系统的几何理论和后继函数方法，得到了脉冲作用下系统具有阶一周期解，且得到了阶一周期解具有存在和唯一性。最后，通过数值模拟来说明我们的理论结果及生物学意义。

Abstract: The destruction of aquatic ecosystems caused by eutrophication of water bodies occurs from time to time. In this paper, a predator-prey system with constant release rate and impulsive state feedback control of predator population in polluted water is studied. By constructing Lyapunov function, a sufficient condition for the global asymptotic stability of the system without impulse effect is ob-tained. By using the geometric theory of semi-continuous dynamic system and the successor func-tion method, it is obtained that the system has order-1 periodic solution under the action of impulse, and the existence and uniqueness of order-1 periodic solution are obtained. Finally, we illustrate our theoretical results and biological significance by numerical simulation.

1. 引言

近年来随着城市的快速发展，工业废水、农业污染、城市生活污水等使得水资源的污染加重，由水体富营养化引起的破坏水生生态系统的现象时有发生。水体富营养化对浮游生物和鱼类都产生了很大的影响，不仅破坏了水生生态系统的平衡和稳定，也影响了水产养殖业的经济效益。因此，调节水生生态系统中各种群之间的相互作用关系，对于维持系统的平衡和稳定具有重要意义。

具有状态反馈脉冲效应的常微分方程被称为半连续动力系统 [1]。近年来，半连续动力系统在科学领域的应用越来越广泛 [2] - [12]。在调控生态系统的过程中，很多学者都通过脉冲状态反馈控制的方法来调节系统的平衡和稳定。例如，Zhang等人提出了一个考虑营养利用率的食饵捕食模型，并通过状态反馈脉冲控制来调节养殖水体的生态平衡 [4]。Guo等人提出并研究了脉冲放鱼和脉冲喷洒杀藻剂的藻鱼系统 [5]。Fu和Chen研究了具有两个状态依赖脉冲控制的水葫芦生态系统，讨论了阶一周期解和阶二周期解的存在性和稳定性 [6]。在此基础上，本文则提出了一个污染水体中具有脉冲状态反馈控制且捕食种群具有常数投放率的食饵捕食系统。

论文其余部分的组织如下。在第二节中，我们将介绍捕食种群具有常数投放率的自由发展模型和相应的状态反馈脉冲模型。在第三节中，对自由发展模型进行了定性分析，并分情况讨论了脉冲状态反馈控制模型的动力学性能。最后一节给出了一些数值模拟并说明了生物学意义。

2. 建立模型

2.1. 自由发展模型

首先我们考虑如下模型。

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=\frac{rx}{k+\omega x}-\lambda xy,\\ \frac{\text{d}y}{\text{d}t}=hxy-fy+u.\end{array}$ (2.1)

假设 $x,y$ 分别表示t时刻食饵(浮游动植物)和捕食者(鱼类)的种群密度；浮游动植物以 $f\left(x\right)=r/\left(k+\omega x\right)$ 的速率增长(其中 $r,k,\omega >0$ )， $k+\omega x$ 是食饵的环境容纳量，且食饵受到自身密度的制约，即随着食饵数量的增加食饵的增长速率逐渐减缓； $\lambda$ 为鱼类对浮游动植物的捕食率；h为鱼类捕食的浮游动植物吸收转化为可存活的鱼类后代的比率；f为水污染导致的鱼类的死亡率；u为鱼类种群的常数投放率。

2.2. 状态反馈脉冲模型

由于水体受有机物污染后，会发生生物降解作用，在分解过程中消耗水中的氧气，使水中的氧气量降低，导致鱼类等会因缺氧而死亡。同时鱼类种群数量减少也会导致水体中浮游动植物大量繁殖；为达到一定的产鱼量，同时也为了避免水华的产生，需要常量投放一些鱼苗来补充受水污染影响而死亡的鱼类数量。但当鱼类种群数量过多时，水体的浮游动植物又不能满足鱼类对食物的需求，所以我们会通过状态反馈脉冲控制的方法捕捞成鱼，进而保持食饵捕食系统的稳定，同时也能产生经济效益。设阈值为m，则得到状态反馈脉冲系统为：

$\{\begin{array}{l}\begin{array}{l}\frac{\text{d}x}{\text{d}t}=\frac{rx}{k+\omega x}-\lambda xy,\\ \frac{\text{d}y}{\text{d}t}=hxy-fy+u,\end{array}\}y<m\\ \Delta y=-\beta y,y=m\end{array}$ (2.2)

当鱼类数量小于m时，食饵捕食系统按照系统(2.2)的前两个方程发展。当鱼类种群数量达到m时，开始以 $\beta \left(0<\beta <1\right)$ 为收获率收获成鱼。根据实际的生物意义，下面的内容我们只考虑系统(2.2)在区域 $\left\{\left(x,y\right)|x>0,y>0\right\}$ 的情况。

3. 模型分析

3.1. 自由发展模型(2.1)的定性分析

模型(2.1)满足如下方程：

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=\frac{rx}{k+\omega x}-\lambda xy=P\left(x,y\right),\\ \frac{\text{d}y}{\text{d}t}=hxy-fy+u=Q\left(x,y\right).\end{array}$ (3.1)

它的垂直等倾线为 $x=0$ 和 $y=r/\left[\lambda (k+\omega x)\right]$，水平等倾线为 $y=u/\left(f-hx\right)$。解方程组(3.1)可知，若 $rf-u\lambda k>0$，则系统(2.1)存在一个边界平衡点 $N_1\left(0,u/f\right)$ 和一个正平衡点 $E\left(x^{*},y^{*}\right)$，其中 $x^{*}=\left(rf-u\lambda k\right)/\left(rh+u\lambda \omega \right)$， $y^{*}=\left(rh+u\lambda \omega \right)/\left(\lambda \omega f+\lambda kh\right)$。接下来，我们假设 $rf-u\lambda k>0$

总是成立。

下面分析平衡点 $N_1\left(0,u/f\right)$ 和 $E\left(x^{*},y^{*}\right)$ 的稳定性。系统(2.1)的雅可比矩阵为：

$J=\left(\begin{array}{cc}\frac{rk}{{\left(k+\omega x\right)}^2}-\lambda y& -\lambda x\\ hy& hx-f\end{array}\right)$ .

系统(2.1)在 $N_1\left(0,u/f\right)$ 点的雅可比矩阵为：

$J\left(N_1\right)=\left(\begin{array}{cc}\frac{r}{k}-\frac{\lambda u}{f}& 0\\ \frac{hu}{f}& -f\end{array}\right)$ .

经计算可得 $Det\left(J\left(N_1\right)\right)=\left(u\lambda k-rf\right)/k<0$，即 $N_1$ 是一个鞍点。

系统(2.1)在 $E\left(x^{*},y^{*}\right)$ 点的雅可比矩阵为：

$J\left(E\right)=\left(\begin{array}{cc}\frac{-r\omega x^{*}}{{\left(k+\omega x^{*}\right)}^2}& -\lambda x^{*}\\ h{y}^{*}& -\frac{u}{y^{*}}\end{array}\right)$ .

经计算可得

$Det\left(J\left(E\right)\right)=\frac{r\omega u{x}^{*}}{{\left(k+\omega x^{*}\right)}^2{y}^{*}}+h\lambda x^{*}{y}^{*}>0$， $Tr\left(J\left(E\right)\right)=-\frac{r\omega x^{*}}{{\left(k+\omega x^{*}\right)}^2}-\frac{u}{y^{*}}<0$ .

显然，点E是一个局部渐进稳定的结点或焦点。

因此，得到结论如下。

定理1若 $rf-u\lambda k>0$，则模型(2.1)的边界平衡点 $N_1\left(0,u/f\right)$ 是一个鞍点，正平衡点 $E\left(x^{*},y^{*}\right)$ 是一个局部渐近稳定的结点或焦点。

接下来讨论系统(2.1)在平衡点E处的全局渐近稳定性。

首先，构造一个Lyapunov函数：

$V\left(x,y\right)=\left[\left(x-x^{*}\right)-x^{*}\ln \left(\frac{x}{x^{*}}\right)\right]+\eta \left[\left(y-y^{*}\right)-y^{*}\ln \left(\frac{y}{y^{*}}\right)\right]$ (3.2)

此函数沿着时间的导数计算为：

$\begin{array}{c}\frac{\text{d}V}{\text{d}t}=\frac{x-x^{*}}{x}\frac{\text{d}x}{\text{d}t}+\eta \frac{y-y^{*}}{y}\frac{\text{d}y}{\text{d}t}\\ =\left(x-x^{*}\right)\left(\frac{r}{k+\omega x}-\lambda y\right)+\eta \left(y-y^{*}\right)\left(hx-f+\frac{u}{y}\right)\end{array}$ (3.3)

考虑平衡点方程：

$\{\begin{array}{l}\frac{r}{k+\omega x^{*}}-\lambda y^{*}=0\\ h{x}^{*}-f+\frac{u}{y^{*}}=0\end{array}$ (3.4)

将方程(3.4)带入(3.3)可得：

$\begin{array}{c}\frac{\text{d}V}{\text{d}t}=\left(x-x^{*}\right)\left(\frac{r}{k+\omega x}-\lambda y-\frac{r}{k+\omega x^{*}}+\lambda y^{*}\right)+\eta \left(y-y^{*}\right)\left(hx-f+\frac{u}{y}-h{x}^{*}+f-\frac{u}{y^{*}}\right)\\ =\left(x-x^{*}\right)\left[\frac{r\omega \left(x^{*}-x\right)}{\left(k+\omega x\right)\left(k+\omega x^{*}\right)}-\lambda \left(y-y^{*}\right)\right]+\eta \left(y-y^{*}\right)\left[h\left(x-x^{*}\right)+\frac{u\left(y^{*}-y\right)}{y{y}^{*}}\right]\\ =-\frac{r\omega {\left(x-x^{*}\right)}^2}{\left(k+\omega x\right)\left(k+\omega x^{*}\right)}+\left(h\eta -\lambda \right)\left(x-x^{*}\right)\left(y-y^{*}\right)-\frac{\eta u{\left(y-y^{*}\right)}^2}{y{y}^{*}}\end{array}$ (3.5)

当 $\eta =\lambda /h$ 时，有 $\frac{\text{d}V}{\text{d}t}=-\frac{r\omega {\left(x-x^{*}\right)}^2}{\left(k+\omega x\right)\left(k+\omega x^{*}\right)}-\frac{\eta u{\left(y-y^{*}\right)}^2}{y{y}^{*}}<0$，故 $\frac{\text{d}V}{\text{d}t}$ 在第一象限是负定的。

因此，我们得到如下结论。

定理2系统(2.1)在平衡点 $E\left(x^{*},y^{*}\right)$ 处是全局渐进稳定的。

3.2. 模型(2.2)的阶一周期解的存在唯一性

根据上述分析，可知系统(2.1)的平衡点E是一个全局渐进稳定的结点或焦点(如图1所示)。接下来我们将要讨论当平衡点E是焦点时，系统(2.2)是否存在阶1周期解以及阶1周期解的唯一性。

Figure 1. Phase diagram of system (2.1). 1) $E\left(x^{*},y^{*}\right)$ be a node point; 2) $E\left(x^{*},y^{*}\right)$ be a focus point

图1. 系统(2.1)的相图。1) $E\left(x^{*},y^{*}\right)$ 是一个结点；2) $E\left(x^{*},y^{*}\right)$ 是一个焦点

情况1. $m\le y^{*}$

当 $m\le y^{*}$ 时，直线 $y=m$ 落在平衡点 $E\left(x^{*},y^{*}\right)$ 的下方，我们有如下结论。

定理3当 $E\left(x^{*},y^{*}\right)$ 为焦点且 $m\le y^{*}$ 时，系统(2.2)存在一个唯一的阶1周期解。

证明：令直线 $y=\left(1-\beta \right)m$ 分别与水平等倾线和垂直等倾线交于 $A\left(x_A,\left(1-\beta \right)m\right)$， $B\left(x_B,\left(1-\beta \right)m\right)$

两点(如图2所示)。

Figure 2. The existence and uniqueness of order-1 periodic solution for $m\le y^{*}$

图2. $m\le y^{*}$ 时阶一周期解的存在和唯一性

首先，我们证明阶一周期解的存在性。从点A出发的轨迹交脉冲集M于点A₁，然后经过脉冲效应点A₁被映射到相集上的点A₂。根据轨迹趋势，点A₂一定落在点A的右侧，于是有 $x_{A_2}>x_A$，且点A₂是点A的后继点，则点A的后继函数为 $f\left(A\right)=x_{A_2}-x_A>0$。从点B出发的轨迹交脉冲集M于点B₁，然后经过脉冲效应点B₁被映射到相集上的点B₂。根据轨迹趋势，点B₂一定落在点B的左侧，于是有 $x_{B_2}<x_B$，且点B₂是点B的后继点，则点B的后继函数为 $f\left(B\right)=x_{B_2}-x_B<0$。根据解的存在唯一性定理，在A和B之间一定存在一个点Q使得 $f\left(Q\right)=0$。因此，从点Q开始的轨迹 $L_{Q{Q}_1}$ 和相应的脉冲线 $\overline{Q_{1}Q}$ 构成系统(2.2)的一个阶一周期解。

接下来，我们考虑阶一周期解的唯一性。在线段AB上任意取两个点 $I\left(x_I,\left(1-\beta \right)m\right)$， $J\left(x_J,\left(1-\beta \right)m\right)$，使点I，J满足 $x_A<x_I<x_J<x_B$。于是从点I出发的轨迹交脉冲集M于点I₁，然后经过脉冲效应被映射到点I₂，则点I₂是点I的后继点。根据轨迹趋势可知， $x_I<x_{I_2}$。同理，从点J出发的轨迹交脉冲集M于点J₁，然后经过脉冲效应被映射到点J₂，则点J₂是点J的后继点。根据轨迹趋势可知， $x_{J_2}<x_J$。由于 $x_I<x_J$，则 $x_I<x_{I_2}<x_{J_2}<x_J$。于是点I，J的后继函数满足：

$f\left(J\right)-f\left(I\right)=\left(x_{J_2}-x_J\right)-\left(x_{I_2}-x_I\right)=\left(x_{J_2}-x_{I_2}\right)-\left(x_J-x_I\right)<0$ .

这表明后继函数 $f\left(A\right)$ 在线段AB上是单调递减的。因此，在线段AB上仅存在一个点Q使得 $f\left(Q\right)=0$，从而系统(2.2)的阶1周期解是唯一的。

情况2. $0<\left(1-\beta \right)m<y^{*}<m$

当 $0<\left(1-\beta \right)m<y^{*}<m$ 时，直线 $y=m$ 在平衡点 $E\left(x^{*},y^{*}\right)$ 的上方，直线 $y=\left(1-\beta \right)m$ 在平衡点 $E\left(x^{*},y^{*}\right)$ 的下方，直线 $y=m$ 与水平等倾线 $\text{d}y/\text{d}t=0$ 相交，设其交点为 $C\left(x_C,m\right)$，并设通过点C的轨迹为L_C，设直线 $y=m$ 与轨线L_C的负延长线的交点为 $C_0\left(x_{C_0},m\right)$，且C₀是直线 $y=m$ 上距离点C最近的点，从C₀到C的轨线段记为 $L_{C_{0}C}$，记由轨线段 $L_{C_{0}C}$ 、直线 $y=m$ 和两坐标轴围成的区域为 $\Omega$ (如图3)，我们得到如下结论。

定理4若 $E\left(x^{*},y^{*}\right)$ 为焦点且 $0<\left(1-\beta \right)m<y^{*}<m$ 时，则存在一个参数 ${\beta }^{*}\in \left(0,1\right)$ 满足：

当 $0<\beta <{\beta }^{*}$ 时，系统(2.2)不存在阶一周期解，且所有的轨迹都将趋于平衡点E；

当 ${\beta }^{*}\le \beta <1$ 时，系统(2.2)存在一个唯一的阶一周期解。

证明：首先证明 ${\beta }^{*}\in \left(0,1\right)$ 的存在性。

我们知道直线 $y=m$ 与水平等倾线 $\text{d}y/\text{d}t=0$ 相交于点 $C\left(x_C,m\right)$，且通过点C的轨迹为L_C，过点C做一条垂直于x轴的直线 $x=x_C$，则直线 $x=x_C$ 与轨迹L_C相交，设其交点为 $C_1\left(x_C,h_{C_1}\right)$，其中C₁是直线 $x=x_C$ 上距离点C最近的点。我们再过点C₁做直线 $y=h_{C_1}$。由于函数 $y=\left(1-\beta \right)m$ 关于 $\beta$ 是单调的，则这里一定存在一个参数 ${\beta }^{*}\in \left(0,1\right)$ 使得 $h_{C_1}=\left(1-{\beta }^{*}\right)m$。于是，当 $\beta ={\beta }^{*}$ 时，从点C₁出发的轨迹交脉冲集M于点C，且经过脉冲效应后，点C被映射到点C₁。因此，轨迹 $L_{C_{1}C}$ 和线段 $\overline{C{C}_1}$ 构成了系统(2.2)的一个阶一周期解，其中轨线 $L_{C_{1}C}$ 表示点C₁和点C之间的部分。

当 $\beta <{\beta }^{*}$ 时，直线 $y=\left(1-\beta \right)m$ 落在直线 $y=\left(1-{\beta }^{*}\right)m$ 的上方，则从某一点出发的轨迹经过有限次脉冲效应将落在轨线 $L_{C_{1}C}$ 的左侧，并最终趋向于平衡点E。因此，当 $\beta <{\beta }^{*}$ 时，系统(2.2)不存在阶一周期解(如图3所示)。

当 ${\beta }^{*}<\beta <1$ 时，直线 $y=\left(1-\beta \right)m$ 落在直线 $y=\left(1-{\beta }^{*}\right)m$ 的下方。因为直线 $x=x_C$ 与直线 $y=\left(1-\beta \right)m$ 相交，我们将其交点表示为 $C_2\left(x_{C_2},\left(1-\beta \right)m\right)$ (如图4所示)。根据轨迹趋势，从点 $\left(x_0,y_0\right)\in \Omega$ 出发的轨迹一定与脉冲集M相交，且交点一定落在点C的右侧，则从点C₂出发的轨迹交脉冲集M于点 $C_3\left(x_{C_3},m\right)$，根据轨迹趋势，交点C₃一定落在点C的右侧。经过脉冲效应点C₃被映射到相集上的点 $C_4\left(x_{C_4},\left(1-\beta \right)m\right)$，则点C₄是点C₂的后继点，且 $x_{C_4}>x_{C_2}$，于是点C₂的后继函数为 $f\left(C_2\right)=x_{C_4}-x_{C_2}>0$。已知直线 $y=\left(1-\beta \right)m$ 与垂直等倾线 $\text{d}x/\text{d}t=0$ 交于点 $B\left(x_B,\left(1-\beta \right)m\right)$，则从点B出发的轨迹与脉冲集M交于点B₁，然后经过脉冲效应点B₁被映射到点B₂，由轨迹趋势可知点B₂一定落在点B的左侧，于是有 $x_{B_2}<x_B$，则点B的后继函数为 $f\left(B\right)=x_{B_2}-x_B<0$。

Figure 3. Existence of ${\beta }^{*}$ and trajectory tendency for $\beta <{\beta }^{*}$ , $m>y^{*}$

图3. ${\beta }^{*}$ 的存在性及 $\beta <{\beta }^{*}$， $m>y^{*}$ 时的轨线趋势

因此，由解的存在唯一性定理可知，在点C₂和B之间一定存在一个点Q，使得 $f\left(Q\right)=0$，即从点Q开始的轨迹 $L_{Q{Q}_1}$ 和相应的脉冲线 $\overline{Q_{1}Q}$ 构成系统(2.2)的一个阶1周期解。系统(2.2)的阶1周期解的唯一性的证明方法类似于定理3。

Figure 4. The order-1 periodic solution for ${\beta }^{*}<\beta <1$ , $m>y^{*}$

图4. ${\beta }^{*}<\beta <1$， $m>y^{*}$ 时的阶一周期解

4. 数值分析和生物学意义

我们选择了 $r,k,\omega ,\lambda ,h,f,u,m,\beta$ 为参数，分情况对系统(2.2)的阶一周期解的存在性进行了数值模拟。令 $r=4,k=2,\omega =2,\lambda =0.8,h=0.65,f=0.3,u=0.2$，通过数值计算，得到点 $E\left(x^{*},y^{*}\right)$ 是一个焦点，其中 $x^{*}=0.3011,y^{*}=1.922$。

令 $m=1.5,\beta =0.3$，此时 $m\le y^{*}$，根据定理3，出现阶一周期解。如图5，我们可以看到系统(2.2)在点 $\left(0.9456,1.054\right)$ 处存在一个阶一周期解。

令 $m=2.21$，此时 $m>y^{*}$。如图6所示，当 ${\beta }^{*}=0.56$ 时在点 $\left(0.3631,0.9726\right)$ 处存在一个阶一周期解；当 $\beta =0.3<{\beta }^{*}$ 时，如图7所示，从初始点 $\left(1.389,1.57\right)$ 开始的轨迹趋向于平衡点 $E\left(x^{*},y^{*}\right)$。当 $\beta =0.75>{\beta }^{*}$ 时，如图8所示，此时在点 $\left(0.6961,0.5307\right)$ 处有一个阶一周期解。以上这些都表明数值理论结果与定理4是一致的。

Figure 5. Phase portrait and time series of system (2.2) for $m\le y^{*}$

图5. $m\le y^{*}$ 时系统(2.2)的相图和时间序列图

Figure 6. Phase portrait and time series of system (2.2) for $\beta ={\beta }^{*}=0.56$ , $m>y^{*}$

图6. $\beta ={\beta }^{*}=0.56$， $m>y^{*}$ 时系统(2.2)的相图和时间序列图

Figure 7. Phase portrait and time series of system (2.2) for $\beta =0.3<{\beta }^{*}$ , $m>y^{*}$

图7. $\beta =0.3<{\beta }^{*}$， $m>y^{*}$ 时系统(2.2)的相图和时间序列图

Figure 8. Phase portrait and time series of system (2.2) for $\beta =0.75>{\beta }^{*}$ , $m>y^{*}$

图8. $\beta =0.75>{\beta }^{*}$， $m>y^{*}$ 时系统(2.2)的相图和时间序列图

根据数值模拟，我们得到如下结论。对于自由发展系统(2.1)，存在一个稳定的正平衡点 $E\left(x^{*},y^{*}\right)$，其中 $x^{*}=\left(rf-u\lambda k\right)/\left(rh+u\lambda \omega \right)$， $y^{*}=\left(rh+u\lambda \omega \right)/\left(\lambda \omega f+\lambda kh\right)$。当食饵的常数投放率u增大时，平衡点 $E\left(x^{*},y^{*}\right)$ 会向左上方移动(如图9所示)。这说明常量投放鱼苗会更多的消耗水体的浮游动植物，有助于预防水华的发生；与此同时，常量投放一些鱼苗还会补充受水污染影响而死亡的鱼类数量，从而达到一定的产鱼量。

Figure 9. Relationship between equilibrium E and constant release rate u

图9. 平衡点E与常数投放率u之间的关系

以上数值模拟表明，具有脉冲状态反馈控制且捕食种群具有常数投放率的食饵捕食系统能够形成阶一周期解，且周期解是轨道渐进稳定的。这说明在一定的条件下，污染水体中的食饵和捕食种群数量能够呈周期性的震荡，且最终能够控制在适当的范围内。因此理论上只要给出合适的阈值水平，就可以通过捕食种群的常数投放率u来控制食饵捕食种群的数量，从而保持水生生态系统的平衡与稳定。

此外，如果考虑节约资源和控制投放成本，我们还提出了一个对捕食种群进行脉冲投放的生物经济模型，这将是下面的研究工作。

NOTES

^*通讯作者。