一类基于阶段结构的蚊群压制模型的动力学分析

Dynamic analysis of a mosquito population suppression model based on stage structure

摘要

围绕蚊媒种群的生物控制问题，构建了基于昆虫不育技术的阶段结构种群动力学模型，系统探讨了不育蚊子在常值连续释放与周期脉冲释放两种策略下的生态调控机制。探讨了影响种群动态的关键阈值条件，揭示了不同释放模式下种群灭绝与持续存活的规律。理论和数值结果表明，适当延长不育蚊子的脉冲释放周期，在确保防控效果的同时，可显著降低现场释放的频次，所得结果可为权衡防控效能与实施成本提供参考。

Abstract

The biological control of mosquito-borne populations was investigated by developing stage-structured population dynamics models based on the sterile insect technique. The ecological regulation mechanisms of constant continuous release and periodic impulsive release of sterile mosquitoes were systematically explored. The key threshold conditions influencing population dynamics were examined， revealing the patterns of population extinction and sustained survival under different release strategies. The theoretical results and numerical simulations indicated that appropriately extending the pulse release interval of sterile mosquitoes could significantly reduce the frequency of on-site releases while maintaining control efficacy， which can provide references for balancing control effectiveness and implementation costs in practical applications.

Graphical abstract

关键词

蚊群压制 / 释放阈值 / 平衡点 / 稳定性

Key words

mosquito population suppression / release threshold / equilibrium / stability

0 引言

在全球气候变化与城市化进程的双重驱动下，蚊媒传染病正加速突破原有地理限制，其传播范围呈现向中高纬度地区持续扩张的态势。在此背景下，由于登革热、寨卡病毒等主要蚊媒传染病的疫苗免疫屏障尚未成熟，因此通过媒介种群调控阻断传播链已成为当前公共卫生防控的关键途径^［1］。相较于传统物理清除和化学灭杀等被动式防控手段，昆虫不育技术（Sterile Insect Technique，SIT）作为第三代生物防控技术，凭借其生态安全性、持续防控效能及无化学残留等核心优势，已被世界卫生组织列为重点推广的媒介生物控制策略。该技术通过系统性释放辐射不育或遗传修饰的雄性蚊虫，破坏目标种群的生殖潜能，最终实现种群数量级衰减。

已有研究者通过构建微分方程模型，深入解析了不育个体与野生种群的动态互作关系^［2-4］。基于关键参数（如释放周期、性别比例、竞争系数）的系统优化，研究者已提出多种精准释放方案，成功实现将目标种群压制至流行病学安全阈值以下的干预目标^［5-7］。这种数学模型指导下的SIT技术，为应对新发蚊媒传染病的跨境传播提供了可量化的防控工具。

基于蚊子种群的发育生物学特征与行为生态学规律，文献［7］提出了基于全生命周期的高维动力学模型，文献［8-9］采用降维建模方法对其进行了简化和分析。特别地，文献［9］通过引入合理假设，创新性构建了如下具有生长阶段的三维系统：

\left\{\begin{array}{l}\frac{\mathrm{d}J}{\mathrm{d}t}=bA(\mathrm{1}-\frac{J}{K})\frac{A}{A+\gamma g}-(c+d_{\mathrm{0}})J,\\ \frac{\mathrm{d}A}{\mathrm{d}t}=vhJ-d_{\mathrm{1}}A,\\ \frac{\mathrm{d}g}{\mathrm{d}t}=\Lambda (t)-d_{\mathrm{2}}g,\end{array}\right.(1)

式中：J=J(t)、A=A(t)和g=g(t)分别表示环境中水生阶段野蚊、成年雌蚊和不育蚊子在t时刻的数量，b为野蚊产卵率，K为环境容纳量，c为孵化率，\gamma表征不育蚊子的交配竞争系数，为转化率，v为雌雄比例参数，d_i（i=\mathrm{0,1},\mathrm{2}）分别对应水生阶段野蚊、成年雌蚊和不育蚊子的死亡率，\Lambda (t)为不育蚊子的释放函数。

文献［6］针对系统（1）构建了时变控制策略的最优控制框架，但其研究局限于有限时间内的脉冲式干预。与之相反，本文将聚焦于该动力系统的长期渐近行为，旨在建立保障种群持续压制的释放强度临界准则。

系统（1）仅刻画了水生阶段野蚊的密度制约效应，未考虑成蚊阶段的资源竞争机制。事实上，成蚊产卵过程存在显著的空间竞争：当栖息地产卵点接近饱和时，成蚊间通过信息素介导产生生殖抑制现象。为此，本文在种群动态方程中引入负指数项{\mathrm{e}}^{-\eta A}量化雌蚊密度依赖的繁殖效率衰减，得到改进后模型如下：

\left\{\begin{array}{l}\frac{\mathrm{d}J}{\mathrm{d}t}=bA\left(\mathrm{1}-\frac{J}{K}\right)\frac{A}{A+\gamma g}{\mathrm{e}}^{-\eta A}-(c+d_{\mathrm{0}})J,\\ \frac{\mathrm{d}A}{\mathrm{d}t}=vhJ-d_{\mathrm{1}}A,\\ \frac{\mathrm{d}g}{\mathrm{d}t}=\Lambda (t)-d_{\mathrm{2}}g\mathrm{。}\end{array}\right.(2)

本文基于模型（2），解析不育蚊子的两种典型释放策略（常值连续释放、周期脉冲释放）对目标种群压制的动力学机制，确定维持野蚊种群灭绝的最低释放强度。

1 无干预情形下野蚊种群发展模型

当未实施不育蚊子干预时，野生蚊群的动态演化遵循自治系统：

\left\{\begin{array}{l}\frac{\mathrm{d}J}{\mathrm{d}t}=bA\left(\mathrm{1}-\frac{J}{K}\right){\mathrm{e}}^{-\eta A}-(c+d_{\mathrm{0}})J,\\ \frac{\mathrm{d}A}{\mathrm{d}t}=vhJ-d_{\mathrm{1}}A\mathrm{。}\end{array}\right.(3)

对于系统（3），很容易得到如下结论：

命题1 对于任意初始条件J(\mathrm{0})>\mathrm{0},A(\mathrm{0})>\mathrm{0}，系统（3）的解在可行域{\mathbf{R}}_{+}^{\mathrm{2}}内保持非负且全局有界。

定义区域

\mathrm{\Omega }=\left\{(J,A)\in {\mathbf{R}}_{+}^{\mathrm{2}}|\mathrm{0}\le J\le K,\mathrm{0}\le A\le \frac{vhK}{d_{\mathrm{1}}}\right\},

易知，该区域为系统（3）的全局正向不变集。

此外，系统（3）存在灭绝平衡点P_{\mathrm{0}}(\mathrm{0,0})。当正平衡点P^{*}(J^{*},A^{*})存在时，必有

{\mathrm{ℛ}}_{\mathrm{0}}\left(\mathrm{1}-\frac{J^{*}}{K}\right)={\mathrm{e}}^{\frac{vh\eta }{d_{\mathrm{1}}}{J}^{*}},A^{*}=\frac{vh}{d_{\mathrm{1}}}{J}^{*},

式中：{\mathrm{ℛ}}_{\mathrm{0}}=\frac{vhb}{d_{\mathrm{1}}(c+d_{\mathrm{0}})}。

显然，当且仅当{\mathrm{ℛ}}_{\mathrm{0}}>\mathrm{1}时，系统（3）的正平衡点P^{*}(J^{*},A^{*})存在，称{\mathrm{ℛ}}_{\mathrm{0}}为野蚊种群的基本子代数。

命题2 对于系统（3），当{\mathrm{ℛ}}_{\mathrm{0}}\le \mathrm{1}时，灭绝平衡点P_{\mathrm{0}}是全局渐近稳定的；当{\mathrm{ℛ}}_{\mathrm{0}}>\mathrm{1}时，正平衡点P^{*}是全局渐近稳定的。

证明 应用Dulac函数B(J,A)=\mathrm{1}，可证系统（3）在正象限{\mathbf{R}}_{+}^{\mathrm{2}}内不存在周期轨道。

首先分析灭绝平衡点P_{\mathrm{0}}(\mathrm{0,0})的雅可比矩阵：

{\mathit{𝒥}}_{P_{\mathrm{0}}}=\left(\begin{array}{cc}-(c+d_{\mathrm{0}})& b\\ vh& -d_{\mathrm{1}}\end{array}\right)\mathrm{。}

通过直接计算可得

由此推得：当{\mathrm{ℛ}}_{\mathrm{0}}\le \mathrm{1}时，行列式\mathrm{d}\mathrm{e}\mathrm{t}{\mathit{𝒥}}_{P_{\mathrm{0}}}>\mathrm{0}，P_{\mathrm{0}}为局部指数渐近稳定的结点；当{\mathrm{ℛ}}_{\mathrm{0}}>\mathrm{1}时，行列式，P_{\mathrm{0}}退化为鞍点。

由系统（3）不存在其他稳定的周期解与正平衡点可知，当{\mathrm{ℛ}}_{\mathrm{0}}\le \mathrm{1}时，灭绝平衡点P_{\mathrm{0}}是全局渐近稳定的。

当{\mathrm{ℛ}}_{\mathrm{0}}>\mathrm{1}时，对正平衡点P^{*}(J^{*},A^{*})进行局部稳定性分析，其雅可比矩阵为

{\mathit{𝒥}}_{P^{*}}=\left(\begin{array}{cc}-u_{\mathrm{1}}-(c+d_{\mathrm{0}})& u_{\mathrm{2}}\\ vh& -d_{\mathrm{1}}\end{array}\right),

式中：

\begin{array}{l}{u}_{\mathrm{1}}=b{A}^{*}/(K{\mathrm{ℛ}}_{\mathrm{0}}(\mathrm{1}-\frac{J^{*}}{K})),\\ u_{\mathrm{2}}=b(\mathrm{1}-\frac{J^{*}}{K})(\mathrm{1}-\eta A^{*})/{\mathrm{ℛ}}_{\mathrm{0}}\mathrm{。}\end{array}

通过计算可得

\mathrm{d}\mathrm{e}\mathrm{t}{\mathit{𝒥}}_{P^{*}}=d_{\mathrm{1}}{u}_{\mathrm{1}}+\frac{\mathrm{1}}{{\mathrm{ℛ}}_{\mathrm{0}}}vbh\eta A^{*}>\mathrm{0}\mathrm{。}

根据Routh-Hurwitz判据，P^{*}在此情形下为稳定的焦点或结点。进一步应用上述类似的方法可证其全局吸引性。证毕。

下面始终假设{\mathrm{ℛ}}_{\mathrm{0}}>\mathrm{1}。

2 常值连续释放策略

当实施不育蚊子的常值连续释放策略（\Lambda (t)\equiv {\Lambda }_{\mathrm{0}}）时，系统（2）退化为如下自治系统：

\left\{\begin{array}{l}\begin{array}{l}\frac{\mathrm{d}J}{\mathrm{d}t}=bA\left(\mathrm{1}-\frac{J}{K}\right)\frac{A}{A+\gamma g}{\mathrm{e}}^{-\eta A}-\\ (c+d_{\mathrm{0}})J,\end{array}\\ \frac{\mathrm{d}A}{\mathrm{d}t}=vhJ-d_{\mathrm{1}}A,\\ \frac{\mathrm{d}g}{\mathrm{d}t}={\Lambda }_{\mathrm{0}}-d_{\mathrm{2}}g\mathrm{。}\end{array}\right.(4)

不育蚊子的方程

\frac{\mathrm{d}g}{\mathrm{d}t}={\Lambda }_{\mathrm{0}}-d_{\mathrm{2}}g,

具有全局渐近稳定的平衡点g^{*}={\Lambda }_{\mathrm{0}}/d_{\mathrm{2}}，因此系统（4）可降维为极限系统：

\left\{\begin{array}{l}\begin{array}{l}\frac{\mathrm{d}J}{\mathrm{d}t}=bA\left(\mathrm{1}-\frac{J}{K}\right)\frac{A}{A+\gamma \frac{{\Lambda }_{\mathrm{0}}}{d_{\mathrm{2}}}}{\mathrm{e}}^{-\eta A}-\\ (c+d_{\mathrm{0}})J,\end{array}\\ \frac{\mathrm{d}A}{\mathrm{d}t}=vhJ-d_{\mathrm{1}}A\mathrm{。}\end{array}\right.(5)

易知，系统（5）的解非负有界，且

\mathrm{\Omega }=\{(J,A)|\mathrm{0}\le J\le K,\mathrm{0}\le A\le \frac{vhK}{d_{\mathrm{1}}}\}

是它的一个正向不变集。此外，灭绝平衡点P_{\mathrm{0}}(\mathrm{0,0})始终存在。

当系统（5）的正平衡点P^{*}(J^{*},A^{*})存在时，

\begin{array}{l}{J}^{*}=\frac{d_{\mathrm{1}}}{vh}{A}^{*},\hspace{1em}\\ {\mathrm{ℛ}}_{\mathrm{0}}\left(\mathrm{1}-\frac{d_{\mathrm{1}}}{Kvh}{A}^{*}\right)A^{*}/(A^{*}+\gamma \frac{{\Lambda }_{\mathrm{0}}}{d_{\mathrm{2}}})={\mathrm{e}}^{\eta A^{*}}\mathrm{。}\end{array}

为分析上述关于A^{*}的方程根的存在性，定义

f_{\mathrm{1}}(x)={\mathrm{ℛ}}_{\mathrm{0}}\frac{(\mathrm{1}-\alpha x)x}{x+\beta },\hspace{1em}{f}_{\mathrm{2}}(x)={\mathrm{e}}^{\eta x},

式中：\alpha =\frac{d_{\mathrm{1}}}{Kvh}，\beta =\gamma \frac{{\Lambda }_{\mathrm{0}}}{d_{\mathrm{2}}}。

通过直接计算，可得函数f_{\mathrm{1}}的性质：

f_{\mathrm{1}}(\mathrm{0})=\mathrm{0},\hspace{1em}\underset{x\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}{f}_{\mathrm{1}}(x)=-\mathrm{\infty },

\begin{array}{l}{f}_{\mathrm{1}}^{{}^{\mathrm{\text{'}}}}(x)={\mathrm{ℛ}}_{\mathrm{0}}\frac{\beta -\mathrm{2}\alpha \beta x-\alpha x^{\mathrm{2}}}{{(x+\beta )}^{\mathrm{2}}},\\ f_{\mathrm{1}}^{{}^{\mathrm{\text{'}}}}(\mathrm{0})=\frac{{\mathrm{ℛ}}_{\mathrm{0}}}{\beta }>\mathrm{0}\mathrm{。}\end{array}

定义临界值

x_{\beta }=-\beta +\frac{\sqrt[]{\alpha \beta (\alpha \beta +\mathrm{1})}}{\alpha },

则有f_{\mathrm{1}}^{{}^{\mathrm{\text{'}}}}(x_{\beta })=\mathrm{0}，函数f_{\mathrm{1}}(x)在区间(\mathrm{0},x_{\beta })单调递增，在(x_{\beta },+\mathrm{\infty })单调递减，于x_{\beta }处达到极大值（见图1）。

显然，存在临界值{\beta }^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}>\mathrm{0}，使得当\beta ={\beta }^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}时，曲线f_{\mathrm{1}}(x)与指数函数曲线f_{\mathrm{2}}(x)={\mathrm{e}}^{\eta x}在点x^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}\in (\mathrm{0},x_{\beta }]处相切，且满足

\left\{\begin{array}{l}{\mathrm{ℛ}}_{\mathrm{0}}\frac{(\mathrm{1}-\alpha x^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}})x^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}}{x^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}+{\beta }^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}}={\mathrm{e}}^{\eta x^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}},\\ {\mathrm{ℛ}}_{\mathrm{0}}\frac{{\beta }^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}-\mathrm{2}\alpha {\beta }^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}{x}^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}-\alpha (x^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}{)}^{\mathrm{2}}}{(x^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}+{\beta }^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}{)}^{\mathrm{2}}}=\eta {\mathrm{e}}^{\eta x^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}}\mathrm{。}\end{array}\right.

通过消去指数项得

\begin{array}{l}{\beta }^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}-\mathrm{2}\alpha {\beta }^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}{x}^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}-\alpha (x^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}{)}^{\mathrm{2}}=\\ \eta x^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}(\mathrm{1}-\alpha x^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}})(x^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}+{\beta }^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}})\mathrm{。}\end{array}

进一步定义两个多项式函数：

\begin{array}{ll}{g}_{\mathrm{1}}(x)& =-\alpha x^{\mathrm{2}}-\mathrm{2}\alpha {\beta }^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}x+{\beta }^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}},\\ g_{\mathrm{2}}(x)& =\eta x(\mathrm{1}-\alpha x)(x+{\beta }^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}})\mathrm{。}\end{array}

分析可知，g_{\mathrm{1}}(x)=\mathrm{0}存在唯一正根，记为x_{{\beta }^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}}；g_{\mathrm{2}}(x)=\mathrm{0}存在3个根，分别为-{\beta }^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}\mathrm{、}\mathrm{0}\mathrm{、}\frac{\mathrm{1}}{\alpha }（见图2）。比较它们根的分布关系，有

\begin{array}{l}\frac{\mathrm{1}}{\alpha }-x_{{\beta }^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}}=\\ \frac{\sqrt[]{\mathrm{1}+\alpha {\beta }^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}}}{\alpha }(\sqrt[]{\mathrm{1}+\alpha {\beta }^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}}-\sqrt[]{\alpha {\beta }^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}})>\mathrm{0},\end{array}

于是一定存在唯一临界点x^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}\in (\mathrm{0},x_{{\beta }^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}}]满足

g_{\mathrm{1}}(x^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}})=g_{\mathrm{2}}(x^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}})。

于是针对系统（5），有如下结果：

定理1 设基本子代数满足{\mathrm{ℛ}}_{\mathrm{0}}>\mathrm{1}，则存在临界释放强度{\Lambda }^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}=d_{\mathrm{2}}{\beta }^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}/\gamma >\mathrm{0}，使得系统（5）的平衡点分布呈现以下特性：

1） 亚临界释放（）：存在两个正平衡点P_{\mathrm{1}}^{*}(J_{\mathrm{1}}^{*},A_{\mathrm{1}}^{*})与P_{\mathrm{2}}^{*}(J_{\mathrm{2}}^{*},A_{\mathrm{2}}^{*})，满足，；

2） 临界释放（{\Lambda }_{\mathrm{0}}={\Lambda }^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}）：存在唯一退化正平衡点P^{*}(J^{*},A^{*})；

3） 超临界释放（{\Lambda }_{\mathrm{0}}>{\Lambda }^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}）：不存在正平衡点。

定理2 当{\mathrm{ℛ}}_{\mathrm{0}}>\mathrm{1}时，系统（5）的动力学行为具有如下阈值响应特性：

1） 种群灭绝状态：当{\Lambda }_{\mathrm{0}}>{\Lambda }^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}时，灭绝平衡点P_{\mathrm{0}}(\mathrm{0,0})全局渐近稳定；

2） 双稳态区域：当时，系统呈现双稳特性，其中P_{\mathrm{0}}和高密度平衡点P_{\mathrm{2}}^{*}是局部渐近稳定的，而低密度平衡点P_{\mathrm{1}}^{*}是不稳定的鞍点；

3） 临界状态：当{\Lambda }_{\mathrm{0}}={\Lambda }^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}时，P_{\mathrm{0}}是局部稳定的，而退化平衡点P^{*}是半稳定的。

证明 类似命题2的论证，系统（5）在正象限{\mathbf{R}}_{+}^{\mathrm{2}}内不存在闭轨线，其雅可比矩阵具有如下形式：

\mathit{𝒥}(J,A)=\left(\begin{array}{cc}-\frac{\mathrm{1}}{K}bA\cdot \frac{A}{a_{\mathrm{0}}}{\mathrm{e}}^{-\eta A}-(c+d_{\mathrm{0}})& {\mathit{𝒥}}_{\mathrm{12}}\\ vh& -d_{\mathrm{1}}\end{array}\right),

式中：

{\mathit{𝒥}}_{\mathrm{12}}=b\left(\mathrm{1}-\frac{J}{K}\right){\mathrm{e}}^{-\eta A}\left(\frac{\mathrm{2}A-\eta A^{\mathrm{2}}}{a_{\mathrm{0}}}-\frac{A^{\mathrm{2}}}{(a_{\mathrm{0}}{)}^{\mathrm{2}}}\right),

a_{\mathrm{0}}=A+\gamma \frac{{\Lambda }_{\mathrm{0}}}{d_{\mathrm{2}}}\mathrm{。}

在P_{\mathrm{0}}(\mathrm{0,0})处，雅可比矩阵简化为

{\mathit{𝒥}}_{P_{\mathrm{0}}}=\left(\begin{array}{cc}-(c+d_{\mathrm{0}})& \mathrm{0}\\ vh& -d_{\mathrm{1}}\end{array}\right),

其特征值，，故P_{\mathrm{0}}为局部渐近稳定的结点。

当{\Lambda }_{\mathrm{0}}>{\Lambda }^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}时，系统（5）无其他稳定平衡点。结合Poincaré-Bendixson定理可得P_{\mathrm{0}}的全局吸引性。

当时，系统（5）存在两个正平衡点P_{\mathrm{1}}^{*}(J_{\mathrm{1}}^{*},A_{\mathrm{1}}^{*})与P_{\mathrm{2}}^{*}(J_{\mathrm{2}}^{*},A_{\mathrm{2}}^{*})。利用J_i^{*}=\frac{d_{\mathrm{1}}}{vh}{A}_i^{*}，雅可比行列式可重构为

\begin{array}{l}\mathrm{d}\mathrm{e}\mathrm{t}{\mathit{𝒥}}_{P_i^{*}}=d_{\mathrm{1}}{a}_i{\mathrm{e}}^{-\eta A_i^{*}}\left(\frac{\mathrm{1}}{K}b{A}_i^{*}+\frac{vhb}{d_{\mathrm{1}}}(\mathrm{1}-\frac{J_i^{*}}{K})\right)-\\ vhb\left(\mathrm{1}-\frac{J_i^{*}}{K}\right){\mathrm{e}}^{-\eta A_i^{*}}{a}_i\left(\mathrm{2}-\eta A_i^{*}-a_i\right)=\\ d_{\mathrm{1}}{a}_i{\mathrm{e}}^{-\eta A_i^{*}}\left(\frac{b}{K}{A}_i^{*}+\frac{vhb}{d_{\mathrm{1}}}(\mathrm{1}-\frac{d_{\mathrm{1}}}{Kvh}{A}_i^{*})-\right.\\ \left.\frac{vhb}{d_{\mathrm{1}}}(\mathrm{1}-\frac{d_{\mathrm{1}}}{Kvh}{A}_i^{*})\left(\mathrm{2}-\eta A_i^{*}-a_i\right)\right)=\\ (c+d_{\mathrm{0}})d_{\mathrm{1}}{A}_i^{*}{\mathrm{e}}^{-\eta A_i^{*}}\frac{\mathrm{d}}{\mathrm{d}A}(f_{\mathrm{1}}-f_{\mathrm{2}}){|}_{A=A_i^{*}},\end{array}

式中：

{\mathrm{e}}^{\eta A_i^{*}}={\mathrm{ℛ}}_{\mathrm{0}}\left(\mathrm{1}-\frac{d_{\mathrm{1}}}{Kvh}{A}_i^{*}\right)a_i,

\begin{array}{l}\frac{\mathrm{d}}{\mathrm{d}A}(f_{\mathrm{1}}-f_{\mathrm{2}}){|}_{A=A_i^{*}}=-{\mathrm{ℛ}}_{\mathrm{0}}\frac{\alpha A_i^{*}}{A_i^{*}+\beta }+\\ {\mathrm{ℛ}}_{\mathrm{0}}(\mathrm{1}-\alpha A_i^{*})\frac{\beta }{(A_i^{*}{+\beta )}^{\mathrm{2}}}-\eta {\mathrm{e}}^{\eta A_i^{*}},\end{array}

a_i=\frac{A_i^{*}}{A_i^{*}+\gamma g^{*}}\mathrm{。}

当时，对于低密度平衡点P_{\mathrm{1}}^{*}，有 \frac{\mathrm{d}}{\mathrm{d}A}(f_{\mathrm{1}}-f_{\mathrm{2}}){|}_{A=A_{\mathrm{1}}^{*}}>\mathrm{0}，即，故P_{\mathrm{1}}^{*}为鞍点。对于高密度平衡点P_{\mathrm{2}}^{*}，有，即\mathrm{d}\mathrm{e}\mathrm{t}{\mathit{𝒥}}_{P_{\mathrm{2}}^{*}}>\mathrm{0}，且，故P_{\mathrm{2}}^{*}为稳定的结点或焦点。

当{\Lambda }_{\mathrm{0}}={\Lambda }^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}时，上述两个正平衡点P_{\mathrm{1}}^{*}和P_{\mathrm{2}}^{*}合并为退化平衡点P^{*}，故P^{*}是半稳定的。证毕。

3 周期脉冲释放策略

当采用周期脉冲释放策略时，模型（2）转化为如下脉冲微分方程：

\left\{\begin{array}{l}\left.\begin{array}{l}\frac{\mathrm{d}J}{\mathrm{d}t}=\frac{b{A}^{\mathrm{2}}{\mathrm{e}}^{-\eta A}}{A+\gamma g}\left(\mathrm{1}-\frac{J}{K}\right)-(c+d_{\mathrm{0}})J,\\ \frac{\mathrm{d}A}{\mathrm{d}t}=vhJ-d_{\mathrm{1}}A,\\ \frac{\mathrm{d}g}{\mathrm{d}t}=-d_{\mathrm{2}}g,\end{array}\right\}t\ne nT,\\ \left.\begin{array}{l}J(n{T}^{+})=J(nT),\\ A(n{T}^{+})=A(nT),\\ g(n{T}^{+})={\Lambda }_n+g(nT),\end{array}\right\}\hspace{1em}t=nT,n\in \mathbf{N},\end{array}\right.(6)

式中：T是不育蚊子的脉冲释放周期，{\Lambda }_n表示第n个周期的不育蚊子释放量。

不育蚊子种群的动态变化由下面子模型（7）描述：

\left\{\begin{array}{l}\frac{\mathrm{d}g(t)}{\mathrm{d}t}=-d_{\mathrm{2}}g(t),\hspace{1em}t\ne nT,\\ g(t^{+})=g(t)+{\Lambda }_n,\hspace{1em}t=nT\mathrm{。}\end{array}\right.(7)

当采用等量周期释放策略（{\Lambda }_n\equiv \Lambda）时，模型（7）存在全局渐近稳定的正周期解

\tilde{g}(t)=\frac{\Lambda {\mathrm{e}}^{-d_{\mathrm{2}}(t-nT)}}{\mathrm{1}-{\mathrm{e}}^{-d_{\mathrm{2}}T}},\hspace{1em}t\in (nT,(n+\mathrm{1})T],

且满足

\underset{t\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\left|g(t)-\tilde{g}(t)\right|=\mathrm{0}\mathrm{。}

由于聚焦于野生蚊群的渐近动力学行为，可以将不育蚊子数量固定为周期函数\tilde{g}(t)，于是得到如下非自治周期系统：

\left\{\begin{array}{l}\frac{\mathrm{d}J}{\mathrm{d}t}=\frac{b{A}^{\mathrm{2}}{\mathrm{e}}^{-\eta A}}{A+\gamma \tilde{g}(t)}(\mathrm{1}-\frac{J}{K})-(c+d_{\mathrm{0}})J,\\ \frac{\mathrm{d}A}{\mathrm{d}t}=vhJ-d_{\mathrm{1}}A\mathrm{。}\end{array}\right.(8)

易知，灭绝平衡点Q_{\mathrm{0}}(\mathrm{0,0})是系统（8）的唯一平衡点，需要确定不育蚊子的释放周期及释放量的取值，使之全局稳定。记

\begin{array}{l}{\Lambda }_p^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}\triangleq \frac{\mathrm{2}b(\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}(d_{\mathrm{2}}T)-\mathrm{1})}{\eta \mathrm{e}\gamma d_{\mathrm{2}}T}\times \\ \mathrm{m}\mathrm{i}\mathrm{n}\left\{\frac{\mathrm{1}}{d_{\mathrm{1}}},\frac{\mathrm{1}}{\mathrm{2}\mathrm{m}\mathrm{i}\mathrm{n}\{c+d_{\mathrm{0}},d_{\mathrm{1}}\}-vh}\right\}\mathrm{。}\end{array}

定理3 对于任意释放周期T>\mathrm{0}，当脉冲释放量满足

\Lambda >{\Lambda }_p^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}

时，系统（8）的灭绝平衡点P_{\mathrm{0}}(\mathrm{0,0})全局渐近稳定。

证明 简单起见，定义一个周期平均效应参数

\begin{array}{l}〈{\tilde{g}}^{-\mathrm{1}}(t)〉=\frac{\mathrm{1}}{T}{\int }_{\mathrm{0}}^T\frac{\mathrm{1}}{\tilde{g}(t)}\mathrm{d}t=\\ \frac{\mathrm{1}-{\mathrm{e}}^{-d_{\mathrm{2}}T}}{T\Lambda }{\int }_{\mathrm{0}}^T{\mathrm{e}}^{d_{\mathrm{2}}t}\mathrm{d}t=\frac{\mathrm{2}(\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}(d_{\mathrm{2}}T)-\mathrm{1})}{d_{\mathrm{2}}T\Lambda }\mathrm{。}\end{array}

利用不等式A{\mathrm{e}}^{-\eta A}\le \mathrm{1}/(\eta \mathrm{e}),对系统（8）进行上界估计可得

\begin{array}{l}\frac{\mathrm{d}(J+A)}{\mathrm{d}t}\le bA\frac{A{\mathrm{e}}^{-\eta A}}{A+\gamma \tilde{g}(t)}-\\ (c+d_{\mathrm{0}})J+vhJ-d_{\mathrm{1}}A\le \\ \left(\frac{b}{\eta \mathrm{e}\gamma \tilde{g}(t)}-d_{\mathrm{1}}\right)A-(c+d_{\mathrm{0}}-vh)J,\end{array}

于是有

\begin{array}{l}\frac{\mathrm{d}(J+A)}{\mathrm{d}t}\le \mathrm{m}\mathrm{a}\mathrm{x}\{-(c+d_{\mathrm{0}}-\\ vh),\frac{b}{\eta \mathrm{e}\gamma \tilde{g}(t)}-d_{\mathrm{1}}\}(J+A),\end{array}

进一步，得到

\begin{array}{l}J((n+\mathrm{1})T)+A((n+\mathrm{1})T)\le \\ (J(nT)+A(nT))\times \\ \mathrm{m}\mathrm{a}\mathrm{x}\{{\mathrm{e}}^{-(c+d_{\mathrm{0}}-vh)T},{\mathrm{e}}^{(\frac{b}{\eta \mathrm{e}\gamma }〈{\tilde{g}}^{-\mathrm{1}}(t)〉-d_{\mathrm{1}})T}\}\mathrm{。}\end{array}

由于，故有c+d_{\mathrm{0}}-vh>\mathrm{0}。如果，那么

也即

\Lambda >\frac{\mathrm{2}b(\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}(d_{\mathrm{2}}T)-\mathrm{1})}{\eta \mathrm{e}\gamma d_{\mathrm{2}}T{d}_{\mathrm{1}}}\mathrm{。}

此时可推出\underset{n\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}(J(nT)+A(nT))=\mathrm{0}，这表明\underset{t\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}(J(t)+A(t))=\mathrm{0}，于是有J(t)+A(t)指数收敛至零，即系统（8）的任意解指数收敛至灭绝平衡点P_{\mathrm{0}}(\mathrm{0,0})。

此外，还可以构造正定函数：

V_{\mathrm{0}}(J,A)=\frac{\mathrm{1}}{\mathrm{2}}(J^{\mathrm{2}}+A^{\mathrm{2}}),

其沿着系统（8）的方向导数满足

\begin{array}{l}{\dot{V}}_{\mathrm{0}}=J\frac{\mathrm{d}J}{\mathrm{d}t}+A\frac{\mathrm{d}A}{\mathrm{d}t}=-(c+d_{\mathrm{0}})J^{\mathrm{2}}-d_{\mathrm{1}}{A}^{\mathrm{2}}+\\ bJA\left(\mathrm{1}-\frac{J}{K}\right)\frac{A{\mathrm{e}}^{-\eta A}}{A+\gamma \tilde{g}(t)}+vhJA\le \\ -\mathrm{m}\mathrm{i}\mathrm{n}\{c+d_{\mathrm{0}},d_{\mathrm{1}}\}(J^{\mathrm{2}}+A^{\mathrm{2}})+\\ JA\left(\frac{b}{\eta \mathrm{e}}\cdot \frac{\mathrm{1}}{A+\gamma \tilde{g}(t)}+vh\right)\le \\ -\mathrm{2}\mathrm{m}\mathrm{i}\mathrm{n}\{c+d_{\mathrm{0}},d_{\mathrm{1}}\}{V}_{\mathrm{0}}+\\ \left(\frac{b}{\eta \mathrm{e}\gamma }\cdot \frac{\mathrm{1}}{\tilde{g}(t)}+vh\right)V_{\mathrm{0}}\le \\ \left(\frac{b}{\eta \mathrm{e}\gamma \tilde{g}(t)}+vh-\mathrm{2}\mathrm{m}\mathrm{i}\mathrm{n}\{c+d_{\mathrm{0}},d_{\mathrm{1}}\}\right)V_{\mathrm{0}}\mathrm{。}\end{array}

若

则。这表明，当\Lambda >{\Lambda }_p^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}时，系统（8）的平衡点P_{\mathrm{0}}(\mathrm{0,0})是全局渐近稳定的。证毕。

4 数值模拟

相关动力学参数基于热带岛屿蚊媒控制实证研究数据^［8-9］，具体取值如下：b=10，v=0.5，γ=1，c=0.05，h=0.01，d₀=0.03，d₁=0.04，d₂=0.12。

环境容纳量参数K的测算基于岛屿生态系统承载力的评估方法^［8-9］。以74 hm²热带岛屿为研究样区，通过估测得野生雄蚊密度为69只/hm²。根据种群性别平衡假设（v=\mathrm{0.5}），稳态下成年雌蚊数量等于雄蚊数量，为A^{*}=\mathrm{69}\times \mathrm{74}=\mathrm{5106}。由此通过水生期的承载力方程求得

K=\frac{d_{\mathrm{1}}{A}^{*}}{vh}/(\mathrm{1}-\frac{\mathrm{1}}{{\mathrm{ℛ}}_{\mathrm{0}}}{\mathrm{e}}^{\eta A^{*}})\approx \mathrm{46}\mathrm{243}\mathrm{。}

它表征了最大环境承载能力。此外，雌蚊的资源竞争系数\eta =\mathrm{3.57}\times {\mathrm{10}}^{-\mathrm{4}}，这借鉴了文献［10］中针对白纹伊蚊的取值。

计算显示，野蚊种群的基本子代数

{\mathrm{ℛ}}_{\mathrm{0}}=\frac{vhb}{d_{\mathrm{1}}(c+d_{\mathrm{0}})}=\mathrm{15.625}>\mathrm{1}\mathrm{。}

在无干预情形下，系统（3）存在全局渐近稳定的共存平衡点P^{*}（见图3）。当采用常值连续释放策略时，依据定理1中给出的临界释放阈值{\Lambda }^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}，可选取，此时，系统（5）呈现双稳态特征，灭绝平衡点P_{\mathrm{0}}与高密度平衡点P_{\mathrm{2}}^{*}都是局部稳定的平衡点（如图4所示）。

针对周期脉冲释放策略，基于每周2至3次的野外不育蚊释放实践数据，选取T=\mathrm{2}与T=\mathrm{4}两种释放周期，在周期性脉冲释放条件下开展对比分析。为建立统一的评估标准，当\mathrm{l}\mathrm{g}(A(t))数值趋近于-1阈值时，即判定野生蚊子种群已达到功能性灭绝状态。这一判定标准主要是为了简化不同释放周期效果的横向比较。

基于定理3的阈值解析方法计算得出，当释放周期分别设定为T=\mathrm{2}和T=\mathrm{4}时，对应的临界释放阈值分别为{\Lambda }_p^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}(\mathrm{2})=\mathrm{3.31}\times {\mathrm{10}}^{\mathrm{4}}和{\Lambda }_p^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}(\mathrm{4})=\mathrm{6.72}\times {\mathrm{10}}^{\mathrm{4}}。

图5和图6中数值模拟结果显示，两种释放方案在灭蚊时效性方面差异不显著：T=\mathrm{2}方案实现种群压制耗时298 d，累计释放不育蚊子\mathrm{4.94}\times {\mathrm{10}}^{\mathrm{6}}只；而T=\mathrm{4}方案仅延长了2 d，压制耗时300 d，总释放量为\mathrm{5.04}\times {\mathrm{10}}^{\mathrm{6}}只，二者控制效果基本相同。

值得关注的是，两种方案在操作成本上呈现显著差异：T=\mathrm{4}方案通过延长释放间隔，将释放事件的实施频次从149次锐减至75次，释放工作量降低约50%。这种周期优化策略在保证防控效能的同时，大幅缩减了人力与物流成本。

4 结束语

通过动力学建模与阈值分析，系统探讨了基于昆虫不育技术（SIT）的蚊媒种群压制策略问题。首先，构建了包含水生期密度制约、成蚊资源竞争与不育蚊子交配竞争的三维阶段结构模型，揭示了种群压制的双稳态特性与临界释放阈值的存在性。然后，证实了基本子代数{\mathrm{ℛ}}_{\mathrm{0}}对自然种群稳定性的核心作用，并推导出连续释放与脉冲释放策略的临界阈值{\Lambda }^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}与{\Lambda }_p^{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}(T)，为蚊子种群压制策略设计提供了量化依据。最后，通过数值模拟，验证了脉冲释放策略通过延长周期至4 d，可在维持等效控制效果（灭绝时间300 d）的同时，减少50%的释放频次，显著降低操作成本，这表明“低频高量”释放模式具有经济性优势，所得理论结果可为热带地区蚊媒控制工程提供“效率-成本”均衡优化的科学参考。

参考文献

[1] CYRANOSKI D.Sterile mosquitoes near take-off［J］.Nature，2008，453： 435.

[2] CHANG Lijie， YU Jianshe， ZHENG Bo. A periodic delay differential equation model for mosquito suppression based on Beverton-Holt-type of birth［J］. The Journal of Geometric Analysis，2025，35： 17.

[3] LUO Huifen， YANG Huichao， YU Jianshe，et al.The dynamics of Wolbachia spread in mosquito population under period-2 perturbation of environments［J］.Advances in Continuous and Discrete Models，2025，2025： 15.

[4] ZHU Zhongcai， FENG Xiaomei， HE Xue，et al. Mirrored dynamics of a wild mosquito population suppression model with Ricker-type survival probability and time delay［J］. Mathematical Biosciences and Engineering，2024，21（2）： 1884-1898.

[5] HUANG Mingzhan， SONG Xinyu， LI Jia. Modelling and analysis of impulsive release of sterile mosquitoes［J］. Journal of Biological Dynamics，2017，11（1）： 147-171.

[6] STRUGAREK M， BOSSIN H， DUMONT Y. On the use of the sterile insect release technique to reduce or eliminate mosquito populations［J］. Applied Mathematical Modelling，2019，68： 443-470.

[7] DUMONT Y， THUILLIEZ J. Human behaviors： A threat to mosquito control？［J］. Mathematical Biosciences，2016，281： 9-23.

[8] ALMEIDA L， DUPREZ M， PRIVAT Y，et al. Optimal control strategies for the sterile mosquitoes technique［J］. Journal of Differential Equations，2022，311： 229-266.

[9] ALMEIDA L， DUPREZ M G， PRIVAT Y，et al. Mosquito population control strategies for fighting against arboviruses［J］. Mathematical Biosciences and Engineering，2019，16（6）： 6274-6297.

[10] BLIMAN P A， CARDONA-SALGADO D， DUMONT Y，et al. Implementation of control strategies for sterile insect techniques［J］. Mathematical Biosciences，2019，314： 43-60.

基金资助

国家自然科学基金项目(12271466)

河南省自然科学基金项目(252300420346)

河南省科技攻关项目(242102110374)