鞍山师范学院，数学与信息科学学院，辽宁 鞍山

收稿日期：2022年6月11日；录用日期：2022年7月6日；发布日期：2022年7月13日

摘要

本文考虑到杀虫剂的残留作用、喷洒前后对种群的作用方式的改变以及天敌资源的有限性，首先建立并系统地研究了资源有限下杀虫剂具有残留作用的固定时刻的害虫控制切换系统，分析关键参数对害虫灭绝阈值的影响；其次建立资源有限下依赖于害虫种群数量的状态害虫控制切换系统，通过数值模拟分析在规定的时间内影响杀虫剂的使用频率的因素。

关键词

害虫综合控制，局部渐近稳定性，全局吸引，资源有限，杀虫剂的残留作用

Mathematical Simulations of Integrated Pest Control with Pesticide Residues under Limited Resources

Lingbi Tang, Xiaohu Liu, Bing Liu^*, Jinyang Li

School of Mathematics and Information Science, Anshan Normal University, Anshan Liaoning

Received: Jun. 11^th, 2022; accepted: Jul. 6^th, 2022; published: Jul. 13^th, 2022

ABSTRACT

Considering the residual effect of pesticides, the change of the mode of action of the population before and after spraying pesticides and the limitation of natural enemy resources, firstly, a fixed-time pest control switching system with residual effect of pesticides under limited resources is established and systematically studied, and the influence of key parameters on the pest eradication threshold is analyzed. Secondly, a state pest control switching system that is dependent on the amount of pest populations under limited resources is established, and the factors that affect the frequency of pesticide use in the specified time are analyzed by numerical simulations.

Keywords:Integrated Pest Management, Local Asymptotic Stability, Global Attractiveness, Limited Resources, Residual Effects of Pesticides

Copyright © 2022 by author(s) and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

目前，害虫仍然是影响世界各国农业生产的主要危害之一。2019年草地贪夜蛾的入侵和2020年东非蝗灾，摧毁了数万公顷农田，对许多国家的农业生产造成了严重损害。已有很多学者利用脉冲动力系统来模拟固定时刻喷洒杀虫剂和释放天敌的过程，研究害虫的综合治理问题 [1] - [6]。但实际上，害虫控制的目的是控制害虫在经济危害水平之下，只有当害虫种群的密度达到经济阈值ET (Economic Threshold)时才实施害虫控制策略，可以利用状态依赖的脉冲动力系统来刻画，这方面的相关工作可以参阅文献 [7] [8] [9]。上述有关固定时刻和状态依赖的害虫控制模型都是假设杀虫剂瞬间、成比例杀死害虫和天敌的，并没有考虑到杀虫剂喷洒以后在相当长的时间里对害虫、天敌仍然持续起作用，在此背景下文献 [10] [11] [12] 引入杀虫剂作用函数，模拟杀虫剂对害虫的残留作用。上述文献均假设天敌是常数释放的，考虑到天敌资源的有限性，Li等在文献 [13] 中建立并研究了非线性脉冲控制的害虫控制模型，给出害虫灭绝周期解全局渐近稳定的条件。本文在上述工作的基础上，综合考虑杀虫剂的残留作用、喷洒前后对种群的作用方式的改变以及天敌资源的有限性，分别建立并研究固定时刻和状态依赖的害虫综合治理切换模型，分析主要参数对害虫控制和喷洒频率的影响。

2. 资源有限下固定时刻害虫综合治理模型的研究

2.1. 模型建立

文献 [2] 建立并研究了一个具有Ivlev型功能性反应函数等周期不同固定时刻喷洒杀虫剂和释放天敌的害虫综合治理模型，但是并没有考虑到杀虫剂的残留作用、喷洒前后对种群的作用方式的改变以及天敌资源有限性，基于此本文建立如下模型：

$\{\begin{array}{l}\begin{array}{l}\frac{\text{d}x\left(t\right)}{\text{d}t}=rx\left(t\right)\left(1-x\left(t\right)\right)-\left(1-{\text{e}}^{-ax\left(t\right)}\right)y\left(t\right),\\ \frac{\text{d}y\left(t\right)}{\text{d}t}=\left(k\left(1-{\text{e}}^{-ax\left(t\right)}\right)-D\right)y\left(t\right),\end{array}\}t\in \left[nT,\left(n+l\right)T\right),\\ \begin{array}{l}\frac{\text{d}x\left(t\right)}{\text{d}t}=rx\left(t\right)\left(1-x\left(t\right)\right)-\left(1-{\text{e}}^{-ax\left(t\right)}\right)y\left(t\right)-b_1\left(t\right)x\left(t\right),\\ \frac{\text{d}y\left(t\right)}{\text{d}t}=\left(k\left(1-{\text{e}}^{-ax\left(t\right)}\right)-D-b_2\left(t\right)\right)y\left(t\right),\end{array}\}t\in \left[\left(n+l\right)T,\left(n+1\right)T\right)\\ \begin{array}{l}x\left(t^{+}\right)=x\left(t\right),\\ y\left(t^{+}\right)=y\left(t\right)+\frac{\sigma }{1+\theta y\left(t\right)},\end{array}\}t=\left(n+1\right)T,n\in N=\left\{0,1,2,\cdots \right\},\end{array},$ (2.1)

其中 $x\left(t\right)$，$y\left(t\right)$ 分别表示在t时刻害虫和天敌的种群数量；a表示天敌对害虫的捕食效率，a表示天敌对害虫的捕食效率，函数 $1-{\text{e}}^{-ax}$ 是Ivlev型功能反应函数；T是释放天敌的周期； $l\left(0<l<1\right)$ 表示在每个周

期T内喷洒杀虫剂的时刻； $\frac{\sigma }{1+\theta y}$ 表示 $t=nT$ 时刻释放天敌的量， $\theta$ 为形状参数；这里

$b_i\left(t\right)=m_i{\text{e}}^{-{\delta }_i\left(t-\left(n+l\right)T\right)},\left(n+l\right)T\le t<\left(n+1\right)T,i=1,2,$

分别表示杀虫剂喷洒以后对害虫和天敌的作用， $m_i\left(i=1,2\right)$ 分别为杀虫剂对害虫和天敌作用的杀死率； ${\delta }_i\left(i=1,2\right)$ 分别为杀虫剂对害虫和天敌作用的衰减率。模型中其他参数生物意义同 [2]。

2.2. 系统(2.1)的动力学性质

对于子系统：

$\{\begin{array}{l}\frac{\text{d}y\left(t\right)}{\text{d}t}=-Dy\left(t\right),t\in \left[nT,\left(n+l\right)T\right),\\ \frac{\text{d}y\left(t\right)}{\text{d}t}=\left(-D-b_2\left(t\right)\right)y\left(t\right),t\in \left[\left(n+l\right)T,\left(n+1\right)T\right),\\ y\left(t^{+}\right)=y\left(t\right)+\frac{\sigma }{1+\theta y\left(t\right)},t=\left(n+1\right)T,\end{array}$ (2.2)

类似于 [13] 中定理2.2的证明，可得

定理2.1 如果 $\theta \ne 0$，系统(2.2)存在一个全局渐近稳定的正T周期解 $y^{*}\left(t\right)$，并且对模型(2.2)的任何解 $y\left(t\right)$ 都有：当 $t\to \infty$ 时， $y\left(t\right)\to y^{*}\left(t\right)$，其中

$y^{*}\left(t\right)=\{\begin{array}{l}{y}_1^{*}\left(t\right),t\in \left[nT,\left(n+l\right)T\right),\\ y_2^{*}\left(t\right),t\in \left[\left(n+l\right)T,\left(n+1\right)T\right),\end{array}$ (2.3)

这里

$y_1^{*}\left(t\right)=y^{*}{\text{e}}^{-D\left(t-nT\right)},y_2^{*}\left(t\right)=y^{*}{\text{e}}^{-D\left(t-nT\right)-\frac{m_2}{{\delta }_2}\left(1-{\text{e}}^{-{\delta }_2\left(t-\left(n+l\right)T\right)}\right)},$

$y^{*}=\frac{-1+\sqrt{1+4\theta \sigma {\text{e}}^{-DT-\frac{m_2}{{\delta }_2}\left(1-{\text{e}}^{-{\delta }_2\left(1-l\right)T}\right)}{\left(1-{\text{e}}^{-DT-\frac{m_2}{{\delta }_2}\left(1-{\text{e}}^{-{\delta }_2\left(1-l\right)T}\right)}\right)}^{-1}}}{2\theta {\text{e}}^{-DT-\frac{m_2}{{\delta }_2}\left(1-{\text{e}}^{-{\delta }_2\left(1-l\right)T}\right)}}.$

因此，由定理2.1知系统(2.2)存在害虫根除周期解 $\left(0,y^{*}\left(t\right)\right)$。

令

$R_0\triangleq \frac{rT}{a{\displaystyle {\int }_0^{lT}{y}_1^{*}\left(t\right)\text{d}t}+{\displaystyle {\int }_{lT}^T\left(y_2^{*}\left(t\right)+b_1\left(t\right)\right)\text{d}t}},$

定理2.2 如果 $R_0<1$，系统(2.2)的害虫根除周期解 $\left(0,y^{*}\left(t\right)\right)$ 是局部渐近稳定的。

证明 设 $\left(x\left(t\right),y\left(t\right)\right)$ 是以 $\left(x\left(0\right),y\left(0\right)\right)$ 为初始值的系统(2.2)的解，则 $\left(0,y^{*}\left(t\right)\right)$ 的局部渐近稳定性可以通过对周期解附近小振幅摄动的解的动态行为来研究。令

$\{\begin{array}{l}\tilde{x}\left(t\right)=x\left(t\right),\\ \tilde{y}\left(t\right)=y\left(t\right)-y^{*}\left(t\right),\end{array}$

则 $\left(\tilde{x}\left(t\right),\tilde{y}\left(t\right)\right)$ 可以写成如下形式

$\left(\begin{array}{l}\tilde{x}\left(t\right)\\ \tilde{y}\left(t\right)\end{array}\right)=\Phi \left(t\right)\left(\begin{array}{l}\tilde{x}\left(0\right)\\ \tilde{y}\left(0\right)\end{array}\right),0\le t<T,$

其中， $\Phi \left(t\right)={\Phi }_1\left(t\right)\cdot {\Phi }_2\left(t\right)$，并且 ${\Phi }_1\left(T\right)$，${\Phi }_2\left(T\right)$ 满足

$\frac{\text{d}{\Phi }_1\left(t\right)}{\text{d}t}=\left(\begin{array}{cc}r-a{y}_1^{*}\left(t\right)& 0\\ ak{y}_1^{*}\left(t\right)& -D\end{array}\right){\Phi }_1\left(t\right),t\in \left[0,lT\right),$

$\frac{\text{d}{\Phi }_2\left(t\right)}{\text{d}t}=\left(\begin{array}{cc}r-a{y}_2^{*}\left(t\right)-b_1\left(t\right)& 0\\ ak{y}_2^{*}\left(t\right)& -D-b_2\left(t\right)\end{array}\right){\Phi }_2\left(t\right),t\in \left[lT,T\right),$

因此

${\Phi }_1\left(t\right)=\left(\begin{array}{cc}{\text{e}}^{{\displaystyle {\int }_0^t\left(r-a{y}_1^{*}\left(t\right)\right)\text{d}t}}& 0\\ {*}_1& {\text{e}}^{-Dt}\end{array}\right),$

${\Phi }_2\left(t\right)=\left(\begin{array}{cc}{\text{e}}^{{\displaystyle {\int }_{lT}^t\left(r-a{y}_2^{*}\left(t\right)-b_1\left(t\right)\right)\text{d}t}}& 0\\ {*}_2& {\text{e}}^{-{\displaystyle {\int }_{lT}^t\left(D+b_2\left(t\right)\right)\text{d}t}}\end{array}\right).$

${*}_1$ 与 ${*}_2$ 的表达式与后面的计算无关。在脉冲点 $t=T$ 处，系统(2.1)的线性化矩阵如下：

$\left(\begin{array}{l}\tilde{x}\left(T^{+}\right)\\ \tilde{y}\left(T^{+}\right)\end{array}\right)=\left(\begin{array}{cc}1& 0\\ 0& 1-\frac{\theta \sigma }{{\left(1+\theta y^{*}\left(T\right)\right)}^2}\end{array}\right)\left(\begin{array}{l}\tilde{x}\left(T\right)\\ \tilde{y}\left(T\right)\end{array}\right),$

从而得单值矩阵

$\begin{array}{c}M=\left(\begin{array}{cc}1& 0\\ 0& 1-\frac{\theta \sigma }{{\left(1+\theta y_p\left(T\right)\right)}^2}\end{array}\right){\Phi }_1\left(lT\right)\cdot {\Phi }_2\left(T\right)\\ =\left(\begin{array}{cc}1& 0\\ 0& 1-\frac{\theta \sigma }{{\left(1+\theta y_p\left(T\right)\right)}^2}\end{array}\right)\left(\begin{array}{cc}{\text{e}}^{{\displaystyle {\int }_0^{lT}\left(r-a{y}_1^{*}\left(t\right)\right)\text{d}t}}& 0\\ {*}_1& {\text{e}}^{-DlT}\end{array}\right)\left(\begin{array}{cc}{\text{e}}^{{\displaystyle {\int }_{lT}^T\left(r-a{y}_2^{*}\left(t\right)-b_1\left(t\right)\right)\text{d}t}}& 0\\ {*}_2& {\text{e}}^{-{\displaystyle {\int }_{lT}^T\left(D+b_2\left(t\right)\right)\text{d}t}}\end{array}\right),\end{array}$

经计算矩阵M的两个特征值分别为：

${\lambda }_1={\text{e}}^{{\displaystyle {\int }_0^{lT}\left(r-a{y}_1^{*}\left(t\right)\right)\text{d}t}+{\displaystyle {\int }_{lT}^T\left(r-a{y}_2^{*}\left(t\right)-b_1\left(t\right)\right)\text{d}t}},$

${\lambda }_2=\left(1-\frac{\theta \sigma }{{\left(1+\theta y^{*}\left(T\right)\right)}^2}\right){\text{e}}^{-{\displaystyle {\int }_{lT}^T\left(D+b_2\left(t\right)\right)\text{d}t}}<1,$

因此，由Floquet理论知，当 $R_0<1$ 时，系统(2.1)的害虫根除周期解 $\left(0,y^{*}\left(t\right)\right)$ 是局部渐近稳定的。

令

$R_1=\frac{M_{s}T}{{\displaystyle {\int }_0^{lT}{y}_1^{*}\left(t\right)\text{d}t}+{\displaystyle {\int }_{lT}^T{y}_2^{*}\left(t\right)\text{d}t}},$ 其中 $M_s=\underset{x\ge 0,t\ge 0}{\sup }\left\{\frac{x\left(t\right)r\left(1-x\left(t\right)\right)}{1-{\text{e}}^{-ax\left(t\right)}}\right\}.$

定理2.3 如果 $R_1<1$，系统(2.1)的害虫根除周期解 $\left(0,y^{*}\left(t\right)\right)$ 是全局吸引的。

证明 由 $R_1<1$，所以存在充分小的 $\epsilon >0$，使得

${\displaystyle {\int }_0^{lT}\left(M_s-\left(y_1^{*}\left(t\right)-\epsilon \right)\right)\text{d}t}+{\displaystyle {\int }_{lT}^T\left(M_s-\left(y_2^{*}\left(t\right)-\epsilon \right)\right)\text{d}t}<0$,

由系统(2.1)有

$\{\begin{array}{l}\frac{\text{d}y\left(t\right)}{\text{d}t}\ge -Dy\left(t\right),t\in \left[nT,\left(n+l\right)T\right),\\ \frac{\text{d}y\left(t\right)}{\text{d}t}\ge \left(-D-b_2\left(t\right)\right)y\left(t\right),t\in \left[\left(n+l\right)T,\left(n+1\right)T\right),\\ y\left(t^{+}\right)=y\left(t\right)+\frac{\sigma }{1+\theta y\left(t\right)},t=\left(n+1\right)T,\end{array}$ (2.4)

由定理2.1及脉冲微分方程比较定理知，当t充分大时， $y\left(t\right)>y^{*}\left(t\right)-\epsilon$，为方便起见，假设对所有 $t\ge 0$，$y\left(t\right)>y^{*}\left(t\right)-\epsilon$ 均成立，由系统(2.1)的第一个和第三个方程，不妨设

$p\left(x\right)=1-{\text{e}}^{-ax},g\left(x\right)=r\left(1-x\right),$

则有

$\frac{\stackrel{\dot{}}{x}}{p\left(x\right)}=\frac{xg\left(x\right)}{p\left(x\right)}-y\left(t\right),t\in \left[nT,\left(n+l\right)T\right),$

$\frac{\stackrel{\dot{}}{x}}{p\left(x\right)}=\frac{x\left(g\left(x\right)-b_1\left(t\right)\right)}{p\left(x\right)}-y\left(t\right),t\in \left[\left(n+l\right)T,\left(n+1\right)T\right),$

当x足够小时， $p\left(x\right)\approx p^{\prime }\left(0\right)\cdot x$，其中， $p^{\prime }\left(0\right)=a>0$，因此积分

${\displaystyle {\int }_{\delta }^{x\left(t\right)}\frac{\text{d}x}{p\left(x\right)}}\to -\infty \iff x\left(t\right)\to 0,$

令

$G\left(x\right)={\displaystyle {\int }_{\delta }^x\frac{\text{d}s}{p\left(s\right)}},$

其中 $\delta >0$，因此当 $t\to \infty$ 时， $G\left(x\right)\to -\infty$，此时，害虫种群趋于灭绝。在 $\left[nT,\left(n+1\right)T\right)$ 上，有

$\begin{array}{c}G\left(x\left(\left(n+1\right)T^{+}\right)\right)={\displaystyle {\int }_{\delta }^{x\left(\left(n+1\right)T^{+}\right)}\frac{\text{d}s}{p\left(s\right)}}={\displaystyle {\int }_{\delta }^{x\left(\left(n+1\right)T\right)}\frac{\text{d}s}{p\left(s\right)}}+{\displaystyle {\int }_{x\left(\left(n+1\right)T\right)}^{x\left(\left(n+1\right)T^{+}\right)}\frac{\text{d}s}{p\left(s\right)}}\\ =G\left(x\left(\left(n+1\right)T\right)\right)+{\displaystyle {\int }_{x\left(\left(n+1\right)T\right)}^{x\left(\left(n+1\right)T^{+}\right)}\frac{\text{d}s}{p\left(s\right)}}=G\left(x\left(n{T}^{+}\right)\right)+{\displaystyle {\int }_{x\left(nT\right)}^{x\left(\left(n+1\right)T\right)}\frac{\text{d}s}{p\left(s\right)}}+{\displaystyle {\int }_{x\left(\left(n+1\right)T\right)}^{x\left(\left(n+1\right)T^{+}\right)}\frac{\text{d}s}{p(\; s\; )}}\end{array}$

$\begin{array}{l}\le G\left(x\left(n{T}^{+}\right)\right)+{\displaystyle {\int }_{nT}^{\left(n+l\right)T}\left(\frac{xg\left(x\right)}{p\left(x\right)}-\left(y_1^{*}\left(t\right)-\epsilon \right)\right)\text{d}t}\\ +{\displaystyle {\int }_{\left(n+l\right)T}^{\left(n+1\right)T}\left(\frac{x\left(g\left(x\right)-b_1\left(t\right)\right)}{p\left(x\right)}-\left(y_2^{*}\left(t\right)-\epsilon \right)\right)\text{d}t}+{\displaystyle {\int }_{x\left(\left(n+1\right)T\right)}^{x\left(\left(n+1\right)T^{+}\right)}\frac{\text{d}s}{p\left(s\right)}}\\ \le G\left(x\left(n{T}^{+}\right)\right)+{\displaystyle {\int }_{nT}^{\left(n+l\right)T}\left(\frac{xg\left(x\right)}{p\left(x\right)}-\left(y_1^{*}\left(t\right)-\epsilon \right)\right)\text{d}t}\\ +{\displaystyle {\int }_{\left(n+l\right)T}^{\left(n+1\right)T}\left(\frac{x\left(g\left(x\right)-b_1\left(t\right)\right)}{p\left(x\right)}-\left(y_2^{*}\left(t\right)-\epsilon \right)\right)\text{d}t},\end{array}$

因此，对 $\forall t>0$，都存在一个正整数m，使得 $t\in \left[mT,\left(m+1\right)T\right)$ 时，

$G\left(x\left(t\right)\right)-G\left(x_0\right)\le m\left({\displaystyle {\int }_0^{lT}\left(M_s-\left(y_1^{*}\left(t\right)-\epsilon \right)\right)\text{d}t}+{\displaystyle {\int }_{lT}^T\left(M_s-\left(y_2^{*}\left(t\right)-\epsilon \right)\right)\text{d}t}\right)+{\displaystyle {\int }_{mT}^t\left(M_s-\left(y^{*}\left(t\right)-\epsilon \right)\right)\text{d}t},$

显然积分 ${\displaystyle {\int }_{mT}^t\left(M_s-\left(y^{*}\left(t\right)-\epsilon \right)\right)\text{d}t}$ 有上界。又 $t\to \infty ,m\to \infty$，如果

${\displaystyle {\int }_0^{lT}\left(M_s-\left(y_1^{*}\left(t\right)-\epsilon \right)\right)\text{d}t}+{\displaystyle {\int }_{lT}^T\left(M_s-\left(y_2^{*}\left(t\right)-\epsilon \right)\right)\text{d}t}<0,$

那么， $G\left(x\right)\to -\infty$，即当 $t\to \infty$ 时， $x\left(t\right)\to 0$。

下面我们证明 $t\to \infty$ 时， $y\left(t\right)\to y^{*}\left(t\right)$。由上面可知， $t\to \infty$ 时， $x\left(t\right)\to 0$。为方便起见，假设对 $t>0$，$p\left(x\right)<\epsilon$ 均成立，从而

$-Dy\left(t\right)\le \frac{\text{d}y\left(t\right)}{\text{d}t}\le \left(k\epsilon -D\right)y\left(t\right),t\in \left[nT,\left(n+l\right)T\right),$

$-\left(D+b_2\left(t\right)\right)y\left(t\right)\le \frac{\text{d}y\left(t\right)}{\text{d}t}\le \left(k\epsilon -D-b_2\left(t\right)\right)y\left(t\right),t\in \left[\left(n+l\right)T,\left(n+1\right)T\right),$

对上面两个不等式右端，考虑下面的微分方程

$\{\begin{array}{l}\frac{\text{d}z\left(t\right)}{\text{d}t}=\left(k\epsilon -D\right)y\left(t\right),t\in \left[nT,\left(n+l\right)T\right),\\ \frac{\text{d}z\left(t\right)}{\text{d}t}=\left(k\epsilon -D-b_2\left(t\right)\right)y\left(t\right),t\in \left[\left(n+l\right)T,\left(n+1\right)T\right),\\ z\left(t^{+}\right)=z\left(t\right)+\frac{\sigma }{1+\theta z\left(t\right)},t=\left(n+1\right)T,\\ z\left(0^{+}\right)=y\left(0^{+}\right),\end{array}$ (2.5)

用与定理2.1相同的方法，系统(2.5)存在全局渐近稳定的周期解：

$z^{*}\left(t\right)=\{\begin{array}{l}{z}_1^{*}\left(t\right),t\in \left[nT,\left(n+l\right)T\right),\\ z_2^{*}\left(t\right),t\in \left[\left(n+l\right)T,\left(n+1\right)T\right),\end{array}$

其中，

$z_1^{*}\left(t\right)=z^{*}{\text{e}}^{-\left(D-k\epsilon \right)\left(t-nT\right)},$

$z_2^{*}\left(t\right)=z^{*}{\text{e}}^{-\left(D-k\epsilon \right)\left(t-nT\right)-\frac{m_2}{{\delta }_2}\left(1-{\text{e}}^{-{\delta }_2\left(t-\left(n+l\right)T\right)}\right)},$

$z^{*}=\frac{-1+\sqrt{1+4\theta \sigma {\text{e}}^{-\left(D-k\epsilon \right)T-\frac{m_2}{{\delta }_2}\left(1-{\text{e}}^{-{\delta }_2\left(1-l\right)T}\right)}{\left(1-{\text{e}}^{-\left(D-k\epsilon \right)T-\frac{m_2}{{\delta }_2}\left(1-{\text{e}}^{-{\delta }_2\left(1-l\right)T}\right)}\right)}^{-1}}}{2\theta {\text{e}}^{-\left(D-k\epsilon \right)T-\frac{m_2}{{\delta }_2}\left(1-{\text{e}}^{-{\delta }_2\left(1-l\right)T}\right)}},$

因此，当 $t\to \infty$ 时， $y\left(t\right)\to y^{*}\left(t\right)$。所以，系统(2.2)的害虫根除周期解 $\left(0,y^{*}\left(t\right)\right)$ 是全局吸引的。

$y^{*}\left(t\right)-{\epsilon }_1\le y\left(t\right)\le z^{*}\left(t\right)+{\epsilon }_1,$

故当 $t>t_1$ 时，令 $\epsilon \to 0$，则有

$y^{*}\left(t\right)-{\epsilon }_1\le y\left(t\right)\le y^{*}\left(t\right)+{\epsilon }_1,$

因此，当 $t\to \infty$ 时， $y\left(t\right)\to y^{*}\left(t\right)$。所以，系统(2.2)的害虫根除周期解 $\left(0,y^{*}\left(t\right)\right)$ 是全局吸引的。

经过简单计算可以得到 $R_0<R_1$，因此有

定理2.4 如果 $R_1<1$，则系统(2.1)的害虫根除周期解 $\left(0,y^{*}\left(t\right)\right)$ 是全局渐近稳定的。

2.3. 主要参数对害虫根除阈值的影响

首先，考虑喷洒杀虫剂时刻和天敌释放量对害虫根除阈值 $R_0$ 的影响。图1(a)和图1(b)分别表明阈值 $R_0$ 关于T是单调增函数，关于天敌释放量 $\sigma$ 是单调减函数。这说明生物控制周期T的减小和天敌释放量的增加有利于害虫控制。其次，图2(a)和图2(b)分别给出 $R_0$ 关于l和T的等高线， $m_1$ 和 $m_2$ 的等高线，结果表明 $R_0$ 分别是关于l的单调增加函数， $m_1$ 的单调减函数， $m_2$ 是单调增函数，因此，越早喷洒杀虫剂，杀虫剂对害虫杀死率越大，杀虫剂对天敌杀死率越小越有利于害虫控制。

图1. (a) 生物控制周期T对阈值R₀的影响，参数取值为r = 0.8，k = 0.8，σ = 10，θ = 25，a = 0.5，D = 0.3，l = 0.7，δ₂ = 0.2，σ = 10，m₁ = 0.7，m₂ = 0.5，δ₁ = 0.5，T = 1；(b) 释放天敌数量σ对阈值R₀的影响，参数取值为r = 0.8，k = 0.8，σ = 10，θ = 5，a = 0.1，D = 0.15，l = 0.7，δ₂ = 0.5，σ = 10，m₁ = 0.5，m₂ = 0.1，δ₁ = 0.1，T = 2

图2. (a) R₀关于T和l的等高线，参数取值为a = 0.2，r = 0.8，k = 0.8，σ = 9.5，θ = 2，D = 0.2，δ₁ = 0.4，δ₂ = 0.1，T = 3，l = 0.1，m₁ = 0.4，m₂ = 0.2；(b) R₀关于m₁和m₂的等高线，参数取值为a = 0.4，r = 0.8，k = 0.6，σ = 5，θ = 14，D = 0.4，δ₁ = 0.5，δ₂ = 0.3，T = 1，l = 0.4，m₁ = 0.2，m₂ = 0.1

图1(a)同时表明阈值 $R_0$ 关于T是单调增函数，因此存在临界值 $T_{\max }$，当 $T<T_{\max }$ 时，害虫根除周期解是局部渐近稳定的(见图3(a))；当 $T>T_{\max }$ 时，害虫种群最终会爆发(见图3(b))。

图3. 系统(3.2)害虫种群和天敌种群的时间序列图，(a) T = 0.5，(b) T = 7，其他参数取值为a = 0.5，r = 0.6，k = 0.8，σ = 3，θ = 10，D = 0.3，δ₁ = 0.4，δ₂ = 0.6，l = 0.6，m₁ = 0.2，m₂ = 0.1，x₀ = 1，y₀ = 2

3. 资源有限下状态依赖害虫综合治理模型的研究

使用杀虫剂能够有效控制害虫的暴发，但是杀虫剂的频繁使用会使害虫产生抗药性，而且对害虫的天敌也产生危害。所以，只有当害虫的种群数量达到ET时，才喷洒杀虫剂。假设害虫种群数量达到ET的时刻分别为 ${\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_j,\cdots$ 且 ${\lambda }_1<{\lambda }_2<\cdots <{\lambda }_j<\cdots$，建立如下模型。

$\{\begin{array}{l}\begin{array}{l}\frac{\text{d}x\left(t\right)}{\text{d}t}=rx\left(t\right)\left(1-x\left(t\right)\right)-\left(1-{\text{e}}^{-ax\left(t\right)}\right)y\left(t\right),\\ \frac{\text{d}y\left(t\right)}{\text{d}t}=\left(k\left(1-{\text{e}}^{-ax\left(t\right)}\right)-D\right)y\left(t\right),\end{array}\}t\in \left[0,{\lambda }_1\right),\\ \begin{array}{l}\frac{\text{d}x\left(t\right)}{\text{d}t}=rx\left(t\right)\left(1-x\left(t\right)\right)-\left(1-{\text{e}}^{-ax\left(t\right)}\right)y\left(t\right)-b_1\left(t\right)x\left(t\right),\\ \frac{\text{d}y\left(t\right)}{\text{d}t}=\left(k\left(1-{\text{e}}^{-ax\left(t\right)}\right)-D-b_2\left(t\right)\right)y\left(t\right),\end{array}\}t>{\lambda }_1,t\ne {\lambda }_j,\\ \begin{array}{l}x\left(t^{+}\right)=x\left(t\right),\\ y\left(t^{+}\right)=y\left(t\right)+\frac{\sigma }{1+\theta y\left(t\right)},\end{array}\}t=\left(n+1\right)T,\\ b_i\left(t\right)=m_i{\text{e}}^{-{\delta }_i\left(t-{\lambda }_j\right)},i=1,2,{\lambda }_j\le t\le {\lambda }_{j+1}.\end{array}$ (3.1)

由于系统(3.1)的理论分析较为复杂，因此下面仅通过数值模拟分析影响杀虫剂喷洒频率的因素。设害虫和天敌的初始值分别为 $x\left(0\right)=x_0<ET$，$y\left(0\right)=y_0$。其他参数意义与系统(2.1)中的相同。

考虑 $t\in \left[0,300\right]$ 杀虫剂的喷洒频率，参数取值为 $a=0.5$，$r=0.6$，$k=0.8$，$\theta =10$，$D=0.3$，${\delta }_1=0.4$，${\delta }_2=0.6$，$l=0.6$，$m_1=0.2$，$m_2=0.1$，$T=4$，$\sigma =4$。在图4(a)中，取在图4(a)中，取 $\left(x\left(0\right),y\left(0\right)\right)=\left(0.3,0.4\right)$ 时，喷洒一次杀虫剂就可以控制害虫在ET之下。当害虫和天敌的初值分别为 $\left(0.6,0.4\right)$ 、 $\left(0.5,0.25\right)$ 和 $\left(0.6,0.3\right)$ 时，从图4(b)、图4(c)和图4(d)中可以看出，三种不同初值状态下杀虫剂分别需要喷洒2次、3次、4次才能成功控制害虫。固定 $\left(x\left(0\right),y\left(0\right)\right)=\left(0.6,0.3\right)$，当 $\sigma =4$ 变为 $\sigma =1.5$ 后，杀虫剂喷洒的次数由4次变为3次(见图4(e))。固定 $\left(x\left(0\right),y\left(0\right)\right)=\left(0.3,0.4\right)$，改变周期令 $T=3$，此时不需要喷洒杀虫剂就可以使害虫种群小于ET (见图4(f))。上述结果表明，杀虫剂的喷洒频率依赖于害虫和天敌种群的初值、天敌释放量和天敌释放周期。

图4. 害虫种群和天敌种群的初始密度、天敌释放量、释放天敌周期对杀虫剂使用次数的影响

4. 结论

本文首先建立了一个资源有限下杀虫剂具有残留的固定时刻害虫综合控制的切换模型。结果表明天敌资源有限情况下，害虫根除周期解局部渐近稳定性并不能保证其全局渐近稳定性。由于 $R_0<R_1$，所以只有当 $R_1<1$ 时，害虫灭绝周期解是全局渐近稳定的。同时数值模拟分析表明缩短释放天敌周期，每周期内越早喷洒杀虫剂，增加天敌的释放量，增强杀虫剂对害虫的杀死率，减小杀虫剂对天敌的作用将有利于害虫控制。其次，为了避免频繁使用杀虫剂，考虑只有当害虫种群数量达到ET时才喷洒杀虫剂，建立了相应的状态依赖的害虫综合控制切换模型。结果表明，在规定时间内如果控制害虫在经济阈值之下，杀虫剂的使用频率依赖于害虫和天敌的初值、天敌的释放量和天敌释放的周期。本文所得结论可以为农业生产部门提供借鉴。

基金项目

国家自然科学基金项目(12171004)，大学生创新创业项目(202110169054)。

文章引用

唐玲碧,刘晓虎,刘 兵,李金洋. 有限资源下杀虫剂具有残留作用的害虫综合控制的数学模拟研究
Mathematical Simulations of Integrated Pest Control with Pesticide Residues under Limited Resources[J]. 应用数学进展, 2022, 11(07): 4509-4518. https://doi.org/10.12677/AAM.2022.117477

参考文献