带有恐惧和脉冲控制的捕食系统动力学分析

Dynamics Analysis of Predation System with Fear and Impulsive Control

摘要: 长期以来，化学农药是小麦蚜虫防治的主要手段。然而，随着小麦种植面积的增大，化学农药被大量使用，害虫抗药性增强、环境污染和生态破坏等问题日益突出。因此，寻找一种更加环保、可持续的防治方法成为当务之急。近年来，生物防治作为一种绿色环保的防治措施逐渐受到关注。使用生物防治手段可以保障农业安全，保护生态环境，从而推动农业可持续发展，促进农业科技创新。结合上述情况，建立一个数学模型来描述小麦生态系统，寻找合适便捷的调控手段，在环境污染无法根除的情况下保证小麦蚜虫的防治和小麦种植的持续健康发展具有十分重要的现实意义。建立状态反馈脉冲微分方程描述小麦蚜虫的防治过程，用脉冲策略来描述喷洒农药杀死小麦蚜虫和人工投放食蚜蝇幼虫，通过对状态反馈脉冲模型中阶一周期解的研究，将所得结论作为小麦蚜虫防治的理论基础，对小麦种植业的发展有重大的指导意义。

Abstract: For a long time, chemical pesticides have been the main means of controlling wheat aphids. However, with the expansion of wheat planting areas, the extensive use of chemical pesticides has led to increasingly prominent problems such as enhanced pest resistance, environmental pollution, and ecological damage. Therefore, finding a more environmentally friendly and sustainable control method has become an urgent task. In recent years, biological control, as a green and environmentally friendly control measure, has gradually attracted attention. The application of biological control methods can ensure agricultural safety, protect the ecological environment, and promote the sustainable development of agriculture and agricultural technological innovation. In light of the above situation, establishing a mathematical model to describe the wheat ecosystem and seeking appropriate and convenient control measures to ensure the control of wheat aphids and the sustainable and healthy development of wheat planting in the context of inevitable environmental pollution is of great practical significance. A state feedback impulsive differential equation is established to describe the process of controlling wheat aphids, and the impulsive strategy is used to describe the spraying of pesticides to kill wheat aphids and the artificial release of aphid larvae. Through the study of the first-order periodic solution in the state feedback impulsive model, the obtained conclusions can serve as the theoretical basis for the control of wheat aphids and have significant guiding significance for the development of the wheat planting industry.

1. 介绍

探索自然界中物种间相互作用的微分方程模型一直是生态学和数学领域的研究热点之一。而其中最引人入胜且经典的模型之一就是捕食者–被捕食者模型，也称为predator-prey模型。该模型主要描述了两个物种之间的相互作用关系，其中一个物种称为捕食者，另一个物种称为被捕食者(猎物)。二者的种群规模变化受到很多因素的影响，例如繁殖行为、捕食行为、病毒感染、自然死亡等。该模型以二者之间的捕食关系为基础展开研究，研究过程中既揭示了捕食者如何影响猎物的种群规模变化，也描述了猎物如何影响捕食者的种群规模变化。通过研究该模型可以较好地解释自然界中的一些现象，为生物多样性保护和害虫防治等问题提供了重要指导。

最早的捕食模型是意大利数学家Alfred J. Lotka和美国数学家Vito Volterra分别在20世纪初提出的，故称为Lotka-Volterra模型。该模型直观地描述了捕食者和猎物之间的种群规模变化，可以帮助我们更好地理解物种间的相互作用关系。该模型如下所示：

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=x\left(r-ax-\beta y-d_1\right)\\ \frac{\text{d}y}{\text{d}t}=y\left(\alpha \beta x-by-d_2\right)\end{array}$ (1.1)

其中，$x\left(t\right)$ 为猎物在t时刻的种群规模，$y\left(t\right)$ 为t时刻捕食者的种群规模，r为猎物种群的内在生长率，a为猎物的种内竞争率，$\beta$ 为捕食者对猎物的攻击率，$d_1$ 和$d_2$ 分别为猎物和捕食者的自然死亡率，$\alpha$ 为猎物的转化率，b为捕食者的种内竞争率。

在对捕食模型的研究过程中，有专家学者发现，捕食者对猎物的影响不仅仅是直接捕食，还存在因捕食者存在猎物而产生的间接影响，我们将这种影响称为恐惧[1]-[4]。猎物因捕食者存在而产生恐惧，从而减少外出觅食的频率或缩小觅食范围，从而影响猎物种群的生长和繁殖。将恐惧因素加入到捕食模型中来考虑，使得捕食模型的研究更具有实际意义、更系统。该模型如下所示：

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=x\left(\frac{r}{1+kx}-ax-\beta y-d_1\right)\\ \frac{\text{d}y}{\text{d}t}=y\left(\alpha \beta x-by-d_2\right)\end{array}$ (1.2)

该模型中，k为猎物因捕食者存在而产生的恐惧系数。

在我国甘肃，大部分家庭都是以小麦种植为生，小麦作为中国人日常需要的主食之一，在我国的需求量日益上升。在人们开始大面积种植小麦来满足国人日益增长的粮食需求时，更大的问题也摆在了面前。在小麦种植过程中，小麦蚜虫是很常见的一种害虫，对小麦的产量和品质都产生了很大的影响[5]-[7]。小麦蚜虫在小麦的整个生育期均可发生，但以穗期危害最为严重。小麦蚜虫通过刺吸麦株茎叶和嫩穗的汁液，导致小麦生长受阻，严重时叶片发黄、枯死，影响小麦产量。其次，小麦蚜虫还能传播小麦黄矮病等病毒病，进一步加剧小麦种植的损失。

长期以来，化学农药是小麦蚜虫防治的主要手段。然而，随着小麦种植面积的增大，化学农药被大量使用，害虫抗药性增强、环境污染和生态破坏等问题日益突出[8] [9]。

因此，寻找一种更加环保、可持续的防治方法成为当务之急。近年来，生物防治作为一种绿色环保的防治措施逐渐受到关注[10]-[15]。利用天敌昆虫、真菌、细菌等生物防治手段来控制蚜虫种群数量，具有高效、低毒、对环境友好的特点。使用生物防治手段可以保障农业安全，保护生态环境，从而推动农业可持续发展，促进农业科技创新。

结合上述情况，建立一个数学模型来描述小麦生态系统，寻找合适便捷的调控手段，在环境污染无法根除的情况下保证小麦蚜虫的防治和小麦种植的持续健康发展具有十分重要的现实意义。建立状态反馈脉冲微分方程描述小麦蚜虫的防治过程，用脉冲策略来描述喷洒农药杀死小麦蚜虫和人工投放食蚜蝇幼虫，通过对状态反馈脉冲模型中阶一周期解的研究[14]-[16]，将所得结论作为小麦蚜虫防治的理论基础，对小麦种植业的发展有重大的指导意义。因此，本文考虑建立一个蚜虫–大灰后食蚜蝇模型，并通过施加状态脉冲的方法进行蚜虫的防治，制定一个控制策略来达到更好的害虫防治效果。根据实际情况，我们建立以下模型：

$\{\begin{array}{l}\begin{array}{l}\frac{\text{d}x}{\text{d}t}=x\left(\frac{r}{1+ky}-ax-\beta y-d_1\right)\\ \frac{\text{d}y}{\text{d}t}=y\left(\alpha \beta x-by-d_2\right)\end{array}\}x<h\\ \begin{array}{l}\Delta x=-px\\ \Delta y=q\end{array}\}x=h\end{array}$ (1.3)

该模型中，p表示只对蚜虫有杀伤力的农药的喷洒率，q表示食蚜蝇幼虫的常数投放量。在害虫防治过程中，h是一个施加脉冲的信号，也就是该系统中蚜虫的最大容纳量。当蚜虫的数量值未达到阈值h时，蚜虫和大灰后食蚜蝇按照系统(1.3)的前两个方程生长繁殖；当蚜虫的数量达到阈值h时，系统(1.3)的后两个方程将采取按比例喷洒农药和常数投放食蚜蝇幼虫措施。

文章的安排如下：我们先对无脉冲动力系统(1.2)进行非负平衡点存在性和稳定性的分析，以及极限环不存在性的分析，然后对脉冲动力系统(1.3)进行阶一周期解的存在性、唯一性和稳定性分析，最后通过数值模拟来验证结果。

2. 无脉冲动力系统的分析

在这一部分，我们将讨论所有平衡点的类型以及正平衡点的存在条件，此外，还会证明系统的有界性和极限环的不存在性。

2.1. 解的正性和有界性

2.1.1. 解的正性

从生物学实际情况来考虑，系统(1.2)的初始条件必须满足$x\left(0\right)\ge 0,y\left(0\right)\ge 0$ 。

从系统(1.2)的第一个方程，我们得到

$x\left(t\right)=x\left(0\right)\exp \left[{\displaystyle {\int }_0^t\left(\frac{r}{1+ky\left(s\right)}-ax\left(s\right)-\beta y\left(s\right)-d_1\right)}\right]$ ,

同理，从系统(1.2)的第二个方程，我们可以得到

$y\left(t\right)=y\left(0\right)\exp \left[{\displaystyle {\int }_o^t\left(\alpha \beta x\left(s\right)-by\left(s\right)-d_2\right)}\right]$

于是，当$x\left(0\right)>0$ 和$y\left(0\right)>0$ 时，有$x\left(t\right)>0$ 和$y\left(t\right)>0$ 。因此，初始点在第一象限内得到的所有解都保持在第一象限内。

2.1.2. 解的有界性

我们由系统(1.2)的第一个方程可以得出$\frac{\text{d}x}{\text{d}t}\le x\left(r-ax\right)$ ，所以由文献[3]中的比较定理可以得出$\underset{t\to +\infty }{\lim }\sup x\left(t\right)\le \frac{r}{a}$ 。

令$W\left(t\right)=x\left(t\right)+y\left(t\right)$ ，则有

$\begin{array}{c}\frac{\text{d}w}{\text{d}t}=\frac{\text{d}x}{\text{d}t}+\frac{\text{d}y}{\text{d}t}\\ =x\left(\frac{r}{1+ky}-ax-\beta y-d_1\right)+y\left(\alpha \beta x-by-d_2\right)\\ \le rx-d_{1}x-d_{2}y\\ \le rx-MW\end{array}$

其中，$M=\min \left\{d_1,d_2\right\}$ 。

由比较定理得：$\underset{t\to +\infty }{\lim }\sup W\left(t\right)\le \frac{r^2}{Ma}$ 。

所以在区域$\Omega =\left\{\left(x,y\right)\in R_{+}^2,W\left(t\right)\le \frac{r^2}{Ma}+\epsilon ,\forall \epsilon >0\right\}$ 中解是有界的。

2.2. 平衡点的局部稳定性

定理1 当$a{d}_1+\alpha \beta \left(d_1-r\right)<0$ 时，系统(1.2)有唯一正平衡点；当$a{d}_1+\alpha \beta \left(d_1-r\right)>0$ 时，系统(1.2)无正平衡点。

证明 由$\frac{\text{d}x}{\text{d}t}=0$ 得到系统(1.2)的垂直等倾线为$x=0,x=\frac{r}{a\left(1+ky\right)}-\frac{\beta y}{a}-\frac{d_1}{a}$ ；由$\frac{\text{d}y}{\text{d}t}=0$ 得到系统(1.2)的水平等倾线为$y=0,x=\frac{by+d_2}{\alpha \beta }$ 。

于是，系统正平衡点满足方程

$\{\begin{array}{l}x=\frac{r}{a\left(1+ky\right)}-\frac{\beta y}{a}-\frac{d_1}{a}\\ x=\frac{by+d_1}{\alpha \beta }\end{array}$

即得到方程

$A_1{y}^2+A_{2}y+A_3=0$ (1.4)

其中，$A_1=kab+k\alpha {\beta }^2>0,A_2=ab+ka{d}_1+\alpha {\beta }^2\text{+}k\alpha \beta d_1>0,A_3=a{d}_1+\alpha \beta \left(d_1-r\right)$ 。

通过计算可以得到$\Delta =A_2^2-4A_1{A}_3>0$ ，即该方程(1.4)有两个根；在计算两根之积$y_1{y}_2=\frac{A_3}{A_1}=\frac{a{d}_1+\alpha \beta \left(d_1-r\right)}{kab+k\alpha {\beta }^2}$ 、两根之和$y_1+y_2=-\frac{A_2}{A_1}<0$ 时，我们可以得到以下结论：

当$a{d}_1+\alpha \beta \left(d_1-r\right)<0$ 时，系统(1.2)有唯一正平衡点；

当$a{d}_1+\alpha \beta \left(d_1-r\right)>0$ 时，系统(1.2)无正平衡点。

为了简便起见，本文我们只研究在条件1) $a{d}_1+\alpha \beta \left(d_1-r\right)<0$ 下的情况，此时存在唯一正平衡点$E_2\left(x^{*},y^{*}\right)$ (其中$y^{*}=\frac{-A_1+\sqrt{A_2^2-4A_1{A}_3}}{2A_1},x^{*}=\frac{b{y}^{*}+d_1}{\alpha \beta }$ )。显然，根据系统(1.2)，我们通过计算还可得到边界平衡点$E_0\left(0,0\right)$ 和$E_{\text{1}}\left(\frac{r-d_1}{a},0\right)$ 。

通过计算，我们可以得到平衡点$E_0\left(0,0\right)$ 的Jacobian矩阵为$J_{E_0}=\left(\begin{array}{cc}r-d_1& 0\\ 0& -d_2\end{array}\right)$ ，得到两个特征根为${\lambda }_{01}=r-d_1>0,{\lambda }_{02}=-d_2<0$ ，则$E_0$ 为一个鞍点；计算得到$E_{\text{1}}\left(\frac{r-d_1}{a},0\right)$ 的Jacobian矩阵为$J_{E_1}=\left(\begin{array}{cc}\left(r-d_1\right)\left(1-\frac{2\alpha }{\beta }\right)& \left(-kr-\beta \right)\frac{r-d_1}{\beta }\\ 0& \alpha \left(r-d_1\right)-d_2\end{array}\right)$ ，得到两个特征根为${\lambda }_{11}=\left(r-d_1\right)\left(1-\frac{2\alpha }{\beta }\right),{\lambda }_{12}=\alpha \left(r-d_1\right)-d_2$ 。

当满足条件(2) $\frac{\beta }{2}<\alpha <\frac{d_2}{r-d_1}$ 时，$E_{\text{1}}$ 为局部渐近稳定的；计算得到$E_2$ 的Jacobian矩阵为$J_{E_0}=\left(\begin{array}{cc}-a{x}^{*}& \frac{-kr{x}^{*}}{{\left(1+k{y}^{*}\right)}^2}-\beta y^{*}\\ \alpha \beta y^{*}& -b{y}^{*}\end{array}\right)$ ，得到$D=ab{x}^{*}{y}^{*}+\alpha \beta y^{*}\left[\frac{kr{x}^{*}}{{\left(1+k{y}^{*}\right)}^2}+\beta y^{*}\right]>0,T=-a{x}^{*}-b{y}^{*}<0$ ，则$E_2$ 是一个局部渐近稳定的结点或焦点。

定理2 系统(1.2)在区域A中不存在极限环。

证明 由前文可知，系统(1.2)的水平等倾线为$y=\frac{\alpha \beta x-d_2}{b}$ ，作一条直线$l_1:x=x_1=\frac{r-d_1}{a}$ 与其水平等倾线交于点$Q\left(\frac{r-d_1}{a},\frac{1}{b}\left(\frac{\alpha \beta \left(r-d_1\right)}{a}-d_2\right)\right)$ ，点Q作平行于x轴的直线$l_2:y=y_1=\frac{1}{b}\left(\frac{\alpha \beta \left(r-d_1\right)}{a}-d_2\right)$ ，因此得到一个由$OPQ{E}_1$ 组成的Bendixson环(如图1所示)。

${\rm A}=\left\{\left(x,y\right)|0\le x\le \frac{r-d_1}{a},0\le y\le \frac{1}{b}\left(\frac{\alpha \beta \left(r-d_1\right)}{a}-d_2\right)\right\}$

Figure 1. Area A

图1. 区域A

对于$y>0$ ，我们有${\frac{\text{d}x}{\text{d}t}|}_{x=x_1}=\frac{r-d_1}{a}\left(\frac{r}{1+ky}-r-\beta y\right)<0$ ，于是系统的轨线总是从$l_1$ 的右边穿向左边；对于$0<x<\frac{r-d_1}{a}$ ，我们有${\frac{\text{d}y}{\text{d}t}|}_{y=y_1}=y_1\alpha \beta \left(x-\frac{r-d_1}{a}\right)<0$ ，于是系统的轨线总是从$l_2$ 的上方穿向下方。因此，系统(1.2)是一致有界的。

定义Dulac函数$D\left(x,y\right)=\frac{1}{xy}$ ，我们将得到

$\frac{\partial \left(DP\right)}{\partial x}+\frac{\partial \left(DQ\right)}{\partial y}=-\frac{a}{y}-\frac{b}{x}<0$

因此，根据Bendixson Dulac定理可得出，系统(1.2)在区域A内不存在极限环。

2.3. 全局稳定性

设定一个Lyapunov函数为：

$g\left(x,y\right)=\left[\left(x-x^{*}\right)-x^{*}\log \frac{x}{x^{*}}\right]+\left[\left(y-y^{*}\right)-y^{*}\log \frac{y}{y^{*}}\right]$

我们可以得到$g\left(x,y\right)>0,g\left(x^{*},y^{*}\right)=0$ ，于是根据系统(1.2)我们可以得到

$\begin{array}{c}\frac{\text{d}g}{\text{d}t}=\frac{x-x^{*}}{x}\frac{\text{d}x}{\text{d}t}+c\frac{y-y^{*}}{y}\frac{\text{d}y}{\text{d}t}\\ =\left(x-x^{*}\right)\left(\frac{r}{1+ky}-ax-\beta y-d_1\right)+c\left(y-y^{*}\right)\left(\alpha \beta x-by-d_2\right)\end{array}$ (1.5)

而由正平衡点的存在性我们可以得到

$\{\begin{array}{l}\frac{r}{1+k{y}^{*}}-a{x}^{*}-\beta y^{*}-d_1=0\\ \alpha \beta x^{*}-b{y}^{*}-d_2=0\end{array}$ (1.6)

结合式(1.5)和式(1.6)我们可以得到

$\begin{array}{c}\frac{\text{d}g}{\text{d}t}=\frac{x-x^{*}}{x}\frac{\text{d}x}{\text{d}t}+c\frac{y-y^{*}}{y}\frac{\text{d}y}{\text{d}t}\\ =-\frac{kr}{\left(1+ky\right)\left(1+k{y}^{*}\right)}\left(x-x^{*}\right)\left(y-y^{*}\right)-\left(\beta -c\alpha \beta \right)\left(x-x^{*}\right)\left(y-y^{*}\right)-a{\left(x-x^{*}\right)}^2-cb{\left(y-y^{*}\right)}^{\text{2}}\end{array}$

如果$c=\frac{1}{\alpha }$ 时，则$\frac{\text{d}g}{\text{d}t}<0$ ，所以正平衡点是全局渐近稳定的。

3. 状态反馈脉冲控制系统的分析

对于正平衡点$E_2$ 和系统中h的位置有两种情况：

1) $h>x^{*}$ ：此时蚜虫的数量未达到可控经济损失最大额时蚜虫的数量，无需人工干预；

2) $h\le x^{*}$ ：正平衡点随着时间的发展蚜虫的数量变多，所造成的经济损失超过了能承受的范围，此时需要人工干预来减少小麦的损失。

在这一节中，我们只考虑情况2)，即当蚜虫的种群密度达到了一定量时，我们通过喷洒农药杀死蚜虫和投放食蚜蝇幼虫的双重手段来控制蚜虫的数量，由此建立以下模型：

$\{\begin{array}{l}\begin{array}{l}\frac{\text{d}x}{\text{d}t}=x\left(\frac{r}{1+ky}-ax-\beta y-d_1\right)\\ \frac{\text{d}y}{\text{d}t}=y\left(\alpha \beta x-by-d_2\right)\end{array}\}x<h,\\ \begin{array}{l}\Delta x=-px\\ \Delta y=q\end{array}\}x=h.\end{array}$ (1.3)

首先建立系统(1.3)的脉冲集$M:x=h$ 和相集$N:x=\left(1-p\right)h$ 。

定理3 系统(1.3)存在阶一周期解。

证明 $B\left(x_B,y_B\right)$ 是相集N上的一点，其满足$0<y_B<\epsilon$ ，$\epsilon$ 是充分小的正数。在系统(1.3)的作用下从B点出发的轨线与脉冲集M相交于$B_1\left(x_{B_1},y_{B_1}\right)$ ，系统(1.3)的脉冲函数将$B_1$ 映射到$B_2\left(\left(1-p\right)x_{B_1},y_{B_1}+q\right)$ (如图2所示)，可以得到

Figure 2. Existence of first-order periodic solutions in the system (1.3)

图2. 系统(1.3)中阶一周期解的存在性

$F\left(B\right)=y_{B_2}-y_B=y_{B_1}+q-\epsilon >0.$ (1.7)

我们在点$B\left(x_B,y_B\right)$ 的上方取一点$A\left(x_A,y_A\right)$ ，点A满足$y_A>y_{E_2}+q$ ，通过对每个区域的轨线分析得出从点A出发的轨线在系统的作用下会与脉冲集交于点$A_1\left(x_{A_1},y_{A_1}\right)$ ，此时$y_{A_1}<y_{E_2}$ ，之后在脉冲函数$\phi (x,y)$ 的作用下到达点$A_2\left(\left(1-p\right)x_{A_1},y_{A_1}+q\right)$ ，可以得到

$F\left(A\right)=y_{A_2}-y_A=y_{A_1}+q-y_A<0.$ (1.8)

根据式(1.7)和(1.8)可得，在点A和B之间一定至少存在一点C使得$F\left(C\right)=0$ 。因此，系统(1.3)存在阶一周期解(如图3所示)。

Figure 3. The first-order periodic solution of the system (1.3)

图3. 系统(1.3)的阶一周期解

定理4 若$h\left({\eta }_1\right)-h\left({\eta }_1^{+}\right)\ge aph$ ，系统(1.3)的阶一周期解是轨道渐近稳定的。

证明 设系统(1.3)的一个阶一周期解是$x=\xi \left(t\right),y=\eta \left(t\right)$ ，周期为T，设

$\begin{array}{l}C=\left(\xi \left(0\right),\eta \left(0\right)\right)=\left({\xi }_0,{\eta }_0\right),\\ C^{\prime }=\left(\xi \left(T\right),\eta \left(T\right)\right)=\left({\xi }_1,{\eta }_1\right)=\left(h,{\eta }_1\right),\\ C^{″}=\left(\xi \left(T^{+}\right),\eta \left(T^{+}\right)\right)=\left({\xi }_1^{+},{\eta }_1^{+}\right)=\left({\xi }_0,{\eta }_0\right)=\left(\left(1-p\right)h,{\eta }_1+q\right),\end{array}$

由系统(1.3)可知

$P\left(x,y\right)=\frac{rx}{1+ky}-a{x}^2-\beta xy-d_{1}x,Q\left(x,y\right)=\alpha \beta xy-b{y}^2-d_{2}y,$

$A\left(x,y\right)=-px,B\left(x,y\right)=q,\phi \left(x,y\right)=x-h.$

于是有

$\frac{\partial P}{\partial x}=\frac{r}{1+ky}-2ax-\beta y-d_1,\frac{\partial P}{\partial y}=\frac{-krx}{{\left(1+ky\right)}^2}-\beta x,$

$\frac{\partial Q}{\partial x}=\alpha \beta y,\frac{\partial Q}{\partial y}=\alpha \beta x-2by-d_2,$

$\frac{\partial A}{\partial x}=-p,\frac{\partial A}{\partial y}=0,\frac{\partial B}{\partial x}=0,\frac{\partial B}{\partial y}=0,\frac{\partial \phi }{\partial x}=1,\frac{\partial \phi }{\partial y}=0.$

根据类Poincare准则，有

$\begin{array}{c}\Delta =\frac{P^{+}\left(\frac{\partial B}{\partial y}\frac{\partial \phi }{\partial x}-\frac{\partial B}{\partial x}\frac{\partial \phi }{\partial y}+\frac{\partial \phi }{\partial x}\right)+Q^{+}\left(\frac{\partial A}{\partial y}\frac{\partial \phi }{\partial x}-\frac{\partial A}{\partial x}\frac{\partial \phi }{\partial y}+\frac{\partial \phi }{\partial y}\right)}{P\frac{\partial \phi }{\partial x}+Q\frac{\partial \phi }{\partial y}}\\ =\frac{P_{+}}{P}\\ =\left(1-p\right)\frac{\frac{r}{1+k{\eta }_1^{+}}-a{\xi }_1^{+}-\beta {\eta }_1^{+}-d_1}{\frac{r}{1+k{\eta }_1}-a{\xi }_1-\beta {\eta }_1-d_1},\end{array}$

${\displaystyle \underset{0}{\overset{T}{\int }}\frac{\partial P}{\partial x}\left(\xi \left(t\right),\eta \left(t\right)\right)\text{d}t}=\ln \left(\frac{{\xi }_1}{{\xi }_0}\right)-{\displaystyle \underset{0}{\overset{T}{\int }}ax\text{d}t}=\ln \left(\frac{{\xi }_1}{\left(1-p\right){\xi }_1}\right)-{\displaystyle \underset{0}{\overset{T}{\int }}ax\text{d}t},$

${\displaystyle \underset{0}{\overset{T}{\int }}\frac{\partial Q}{\partial x}\left(\xi \left(t\right),\eta \left(t\right)\right)\text{d}t}=0.$

于是得出

$\begin{array}{c}\mu =\Delta \cdot \exp {\displaystyle \underset{0}{\overset{T}{\int }}\left(\frac{\partial P}{\partial x}\left(\xi \left(t\right),\eta \left(t\right)\right)+\frac{\partial Q}{\partial x}\left(\xi \left(t\right),\eta \left(t\right)\right)\right)\text{d}t}\\ \le \frac{P^{+}}{P}\\ \le \frac{\frac{r}{1+k{\eta }_1^{+}}-a{\xi }_1^{+}-\beta {\eta }_1^{+}-d_1}{\frac{r}{1+k{\eta }_1}-a{\xi }_1-\beta {\eta }_1-d_1}.\end{array}$

根据系统(1.3)的性质可以得出等式上下都为正，令

$h\left(y\right)=\frac{r}{1+ky}-\beta y,$

于是，当系统(1.3)满足$x=\frac{10}{1+0.4y}-4.5y-1.5$

$h\left({\eta }_1\right)-h\left({\eta }_1^{+}\right)\ge aph,$

可以得到$\Vert \mu \Vert \le 1$ ，系统(1.3)的阶一周期解是轨道渐近稳定的。

4. 数值模拟及结论

在本节中，我们选取参数$k=0.4$ ，$\alpha =0.8$ ，$\beta =0.9$ ，$d_1=0.3$ ，$a=0.2$ ，$b=1$ ，$d_2=1.5$ ，$r=2$ ，使用MATLAB进行数值模拟，来验证上述得到的理论结果。

首先在没有脉冲的条件下，选取初始点(10, 3)以及上述参数条件进行MATLAB仿真，图4是无脉冲系统(1.2)的时间序列图。虽然它是稳定的，但是无法阻止蚜虫对小麦系统的侵害，无法提高小麦种植的经济效益。

之后在脉冲控制的条件下，我们通过数值模拟验证了系统(1.3)阶一周期解的存在性和稳定性(图5~8)。

Figure 4. x and y with respect to t time series diagram without impulses

图4. 无脉冲系统x和y关于t时间序列图

(a) (b)

(c)

Figure 5. Phase diagrams and time series diagrams of system (1.3) with the initial point being (0.5, 0.8), p = 0.5, and q = 0.3: (a) Time series diagram of x with respect to t; (b) Time series diagram of y with respect to t; (c) Phase diagram of the system (1.3)

图5. 系统(1.3)在初始点为(0.5, 0.8)、p = 0.5、q = 0.3条件下的相图和时间序列图：(a) x关于t的时间序列图；(b) y关于t的时间序列图；(c) 系统(1.3)的相图

(a) (b)

(c)

Figure 6. Phase diagrams and time series diagrams of system (1.3) with the initial point being (0.5, 0.8), p = 0.4, and q = 0.3: (a) Time series diagram of x with respect to t; (b) Time series diagram of y with respect to t; (c) Phase diagram of the system (1.3)

图6. 系统(1.3)在初始点为(0.5, 0.8)、p = 0.4、q = 0.3条件下的相图和时间序列图：(a) x关于t的时间序列图；(b) y关于t的时间序列图；(c) 系统(1.3)的相图

(a) (b)

(c)

Figure 7. Phase diagrams and time series diagrams of system (1.3) with the initial point being (0.5, 0.8), p = 0.5, and q = 0.2: (a) Time series diagram of x with respect to t; (b) Time series diagram of y with respect to t; (c) Phase diagram of the system (1.3)

图7. 系统(1.3)在初始点为(0.5, 0.8)、p = 0.5、q = 0.2条件下的相图和时间序列图：(a) x关于t的时间序列图；(b) y关于t的时间序列图；(c) 系统(1.3)的相图

(a) (b)

(c)

Figure 8. Phase diagrams and time series diagrams of system (1.3) with the initial point being (0.65, 0.8), p = 0.5, and q = 0.3: (a) Time series diagram of x with respect to t; (b) Time series diagram of y with respect to t; (c) Phase diagram of the system (1.3)

图8. 系统(1.3)在初始点为(0.65, 0.8)、p = 0.5、q = 0.3.条件下的相图和时间序列图：(a) x关于t 的时间序列图；(b) y关于t的时间序列图；(c) 系统(1.3)的相图

根据上述对四组数据的数值模拟结果，我们进行更深入的分析和讨论。首先在图像中可以清晰地看出，阶一周期解的存在性与参数p、q以及初值条件无关。这一发现意味着，无论在模型中如何调整这些参数的初始设定，阶一周期解都会稳定地存在。然而，尽管阶一周期解的存在性与这些参数无关，但其具体的位置与参数q密切相关。这表明，通过调整q值，可以在一定程度上影响周期解在相空间中的具体位置，从而为系统的调控提供了一定的自由度。这一发现对于理解和优化生物控制系统有重要的意义。

此外通过对这些数据的分析，我们发现在恐惧因子的影响下，蚜虫和大灰后食蚜蝇的种群数量都受到了一定影响，从而影响生物系统的动态平衡。当蚜虫的数量增加到了一定阈值，我们必须进行人工干预。通过喷洒农药来杀死蚜虫，可以迅速降低蚜虫的种群数量，缓解其对小麦的侵害；同时投放食蚜蝇幼虫可以增强天敌的捕食能力，进一步抑制蚜虫的数量增长，使两个物种在该系统中能够达到一个动态平衡，以此来减少小麦蚜虫对小麦的侵害，提高小麦种植基地的经济效益，促进祖国农业的发展。同时，根据实际意义，q值不宜过大，因为投放食蚜蝇幼虫不但需要时间成本来等待幼虫长大具备捕食能力，而且食蚜蝇幼虫长大为成虫也有一定的损失率。这意味着，如果我们过度依赖通过投放食蚜蝇幼虫来控制蚜虫的数量，会面临高成本、低效率的问题。因此，在选择q值时，需要综合实际效果和成本问题，找到一个最优的解决方案。

NOTES

^*通讯作者。