一类具有饱和发生率的SEIQV网络蠕虫模型的稳定性分析

Stability Analysis of a Class of SEIQV Network Worm Model with Saturated Incidence Rate

摘要

本文建立了一类考虑隔离和安装杀毒软件策略的具有饱和发生率的SEIQV网络蠕虫模型， 并利用下一代矩阵法得到了模型的基本再生数R_{\mathrm{0}}。以\mathrm{L}\mathrm{a}\mathrm{S}\mathrm{a}\mathrm{l}\mathrm{l}\mathrm{e}不变集原理证明了当R_{\mathrm{0}}\le \mathrm{1}时无蠕虫平衡点E^{\mathrm{0}}是全局渐近稳定的； 以\mathrm{H}\mathrm{u}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{z}判据证明了当R_{\mathrm{0}}>\mathrm{1}时正平衡点E^{\mathrm{*}}是局部渐近稳定的； 利用\mathrm{L}\mathrm{i}\mathrm{­}\mathrm{M}\mathrm{u}\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{n}\mathrm{e}\mathrm{y}几何方法证明了正平衡点E^{\mathrm{*}}的全局渐近稳定性。通过数值模拟验证了所得结论的正确性， 同时证明， 减少感染设备与其他设备的连接， 提高人们对杀毒软件的认识， 可有效控制蠕虫病毒的传播。

Abstract

A type of network worm model with saturated incidence rate was established， this model took into account isolation and antivirus software installation strategies. The basic reproduction number of the model was derived using the next generation matrix method； The LaSalle invariant set principle was utilized to prove that the worm free equilibrium point E^{\mathrm{0}} was globally asymptotically stable when R_{\mathrm{0}}\le \mathrm{1}. The Hurwitz criterion was employed to demonstrate that the positive equilibrium point E^{\mathrm{*}} was locally asymptotically stable when R_{\mathrm{0}}>\mathrm{1}.And the Li-Mulowney geometric method was used to prove the global asymptotic stability of the positive equilibrium point E^{\mathrm{*}}. The correctness of these conclusions was verified through numerical simulation， and it was found that reducing the connectivity between infected devices with other devices， as well as enhancing individuals’ awareness of antivirus software， could significantly curb the spread of worm viruses.

Graphical abstract

关键词

网络蠕虫 / 基本再生数 / 隔离 / 稳定性 / 平衡点

Key words

worm / basic regeneration number / isolation / stability / equilibrium point

0 引　言

网络蠕虫病毒是一种可以自我复制并感染其他计算机的计算机程序^［1］， 其具有生物病毒的一些特征， 可以像生物病毒在人与人之间传播一样从一台计算机传播到另一台计算机， 网络蠕虫的恶意传播对网络安全造成了严重威胁。

自从Kephart等^［2‐3］对蠕虫病毒的传播进行数学建模以来， 许多学者对不同的数学模型^［4‐11］进行了研究。Badshah等^［4］建立了一类具有非线性关联函数的SEIQRS模型， 得到了基本再生数和平衡点， 并证明了平衡点的稳定性。Xiao等^［5］研究了一类具有隔离的SEIQR蠕虫传播模型， 发现隔离（即切断感染设备与其它设备的联系）能有效减少感染设备的影响。Kim等^［6］建立了包含计算机和USB外接设备两个子框架的计算机病毒传播模型， 发现提高完全恢复率能够有效减少感染计算机和感染USB设备的数量。Gao等^［7］研究了一类VEIQS网络蠕虫模型， 发现安装杀毒软件能够有效减少感染设备的数量。Wang等^［8］建立了具有饱和发生率的计算机蠕虫模型， 研究了计算机蠕虫在利用USB传输数据过程中的动态行为。Madhusudanan等^［9］建立了具有holling型发生率的SVEIR计算机蠕虫模型， 证明了Hopf分岔的存在， 并研究了Hopf分岔的性质。Wang等^［10］建立了具有两个子框架的网络蠕虫模型， 研究了安装杀毒软件策略后计算机与可移动设备之间交互动力学的影响。王瑞玲等^［11］建立了具有饱和发生率的网络蠕虫模型， 发现隔离能有效减少感染设备的数量。

上述文献并未将饱和发生率和隔离以及安装杀毒软件策略结合起来， 且大多采用构造Lyapunov函数法来证明正平衡点的全局稳定性， 本文将其结合起来， 建立了一类具有饱和发生率的SEIQV网络蠕虫模型， 并采用Li-Mulowney几何方法^［12］证明了正平衡点的全局稳定性。

1 模型的建立

将计算机种群N\left(t\right)分为\mathrm{5}个状态： 易感状态S\left(t\right)， 潜伏状态E\left(t\right)， 感染状态I\left(t\right)， 隔离状态Q\left(t\right)， 安全状态V\left(t\right)。其中， S\left(t\right)， E\left(t\right)， I\left(t\right)， Q\left(t\right)， V\left(t\right)分别表示在t时刻处于易感状态、 潜伏状态、 感染状态、 隔离状态、 安全状态的设备的数量。 蠕虫传播过程如图 1 所示。

图 1 中， b为外部设备连入网络的速率； \mu为自然淘汰率； \omega为杀毒软件的覆盖率； \beta为感染设备对易感设备的传染率； m为杀毒软件的生效率； \sigma为从潜伏状态到感染状态的转移率； \gamma为从感染状态到免疫状态的转移率； n为感染状态计算机修复的成功率； \phi为从隔离状态到安全状态的转移率； q_{\mathrm{1}}， q_{\mathrm{2}}， q_{\mathrm{3}}分别为易感状态、 潜伏状态、 感染状态的隔离率。参数b， \mu， \omega， m， q_{\mathrm{1}}， q_{\mathrm{2}}， q_{\mathrm{3}}， \sigma， \gamma， \phi均为不大于1的正数， \beta， \alpha均大于等于0。假设不考虑计算机因病毒入侵所导致的淘汰率， 处于安全状态的计算机有永久免疫。

建立模型

\left\{\begin{array}{l}\frac{\mathrm{d}S}{\mathrm{d}t}=b-\frac{\beta SI}{\mathrm{1}+\alpha I}-\left(m\omega +\mu +q_{\mathrm{3}}\right)S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}=\frac{\beta SI}{\mathrm{1}+\alpha I}-\left(\sigma +\mu +q_{\mathrm{2}}\right)E+\left(\mathrm{1}-n\right)\gamma I,\\ \frac{\mathrm{d}I}{\mathrm{d}t}=\sigma E-\left(\mu +\gamma +q_{\mathrm{1}}\right)I,\\ \frac{\mathrm{d}Q}{\mathrm{d}t}=q_{\mathrm{3}}S+q_{\mathrm{2}}E+q_{\mathrm{1}}I-\left(\mu +\phi \right)Q,\\ \frac{\mathrm{d}V}{\mathrm{d}t}=m\omega S+\phi Q+n\gamma I-\mu V\mathrm{。}\end{array}\right.（1）

令N(t)=S(t)+E(t)+Q(t)+I(t)+R(t)， N\left(t\right)为t时刻移动设备的总数。

将系统（1）中的等式相加得

\frac{\mathrm{d}N\left(t\right)}{\mathrm{d}t}=b-\mu N\left(t\right),

从而可得，

N\left(t\right)=\frac{b}{\mu }+\left(N\left(\mathrm{0}\right)-\frac{b}{\mu }\right){\mathrm{e}}^{-\mu t},

故 \underset{t\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}N\left(t\right)=\frac{b}{\mu }\mathrm{。}

由此可得系统（1）的正向不变集为

\begin{array}{c}\Omega =\\ \left\{\left(S,E,I,Q,V\right)\in R_{+}^{\mathrm{5}}\left|\mathrm{0}\le S+E+I\right.+Q+V\le \frac{b}{\mu }\right\}\mathrm{。}\end{array}

2 基本再生数及平衡点的存在性

通过计算， 系统（1）的无蠕虫平衡点为E^{\mathrm{0}}\left(S_{\mathrm{0}},\mathrm{0,0},Q_{\mathrm{0}},V_{\mathrm{0}}\right)， 其中

S_{\mathrm{0}}=\frac{b}{m\omega +\mu +q_{\mathrm{3}}}， Q_{\mathrm{0}}=\frac{q_{\mathrm{3}}{S}_{\mathrm{0}}}{\mu +\phi }，

V_{\mathrm{0}}=\frac{m\omega S_{\mathrm{0}}+\phi Q_{\mathrm{0}}}{\mu }\mathrm{。}

利用下一代矩阵法^［13］求基本再生数R_{\mathrm{0}}， 则

\mathit{ℱ}=\left[\begin{array}{c}\frac{\beta SI}{\mathrm{1}+\alpha I}\\ \mathrm{0}\end{array}\right]，

\mathit{𝒱}=\left[\begin{array}{c}\left(\sigma +\mu +q_{\mathrm{2}}\right)E-\left(\mathrm{1}-n\right)\gamma I\\ -\sigma E+\left(\mu +\gamma +q_{\mathrm{1}}\right)I\end{array}\right]，

F=\left(\begin{array}{cc}\mathrm{0}& \beta S_{\mathrm{0}}\\ \mathrm{0}& \mathrm{0}\end{array}\right)，

\mathit{V}=\left[\begin{array}{cc}\sigma +\mu +q_{\mathrm{2}}& -\left(\mathrm{1}-n\right)\gamma \\ -\sigma & \mu +\gamma +q_{\mathrm{1}}\end{array}\right]，

{\mathit{V}}^{-\mathrm{1}}=\left[\begin{array}{cc}\frac{\mu +\gamma +q_{\mathrm{1}}}{A}& \frac{\left(\mathrm{1}-n\right)\gamma }{A}\\ \frac{\sigma }{A}& \frac{\sigma +\mu +q_{\mathrm{2}}}{A}\end{array}\right]，

令

\begin{array}{c}A=\left|\mathit{V}\right|=\\ \left(\sigma +\mu +q_{\mathrm{2}}\right)\left(\mu +\gamma +q_{\mathrm{1}}\right)-\sigma \left(\mathrm{1}-n\right)\gamma ，\end{array}

\mathit{F}{\mathit{V}}^{-\mathrm{1}}=\left(\begin{array}{cc}\beta S_{\mathrm{0}}\frac{\sigma }{A}& \beta S_{\mathrm{0}}\frac{\sigma +\mu +q_{\mathrm{2}}}{A}\\ \mathrm{0}& \mathrm{0}\end{array}\right)\mathrm{。}

基本再生数R₀是 FV^-1的谱半径

\begin{array}{c}{R}_{\mathrm{0}}=\rho \left(F{V}^{-\mathrm{1}}\right)=\beta S_{\mathrm{0}}\frac{\sigma }{A}=\\ \frac{\sigma \beta b}{\left(m\omega +\mu +q_{\mathrm{3}}\right)\left[\left(\sigma +\mu +q_{\mathrm{2}}\right)\left(\mu +\gamma +q_{\mathrm{1}}\right)-\sigma \left(\mathrm{1}-n\right)\gamma \right]}\mathrm{。}\end{array}

定理 1　当R_{\mathrm{0}}>\mathrm{1}时， 系统（1）在\Omega内存在唯一的正平衡点E^{\mathrm{*}}\left(S_{\mathrm{*}},E_{\mathrm{*}},I_{\mathrm{*}},Q_{\mathrm{*}},V_{\mathrm{*}}\right)。

证明　通过计算， 可得

\left\{\begin{array}{l}{S}_{\mathrm{*}}=\frac{A\left(\mathrm{1}+\alpha I_{\mathrm{*}}\right)}{\sigma \beta },\\ I_{\mathrm{*}}=\frac{\sigma \beta b-\left(m\omega +\mu +q_{\mathrm{3}}\right)A}{A\left[\beta +\alpha \left(m\omega +\mu +q_{\mathrm{3}}\right)\right]},\\ E_{\mathrm{*}}=\frac{\mu +\gamma +q_{\mathrm{1}}}{\sigma }{I}_{\mathrm{*}},\\ Q_{\mathrm{*}}=\frac{q_{\mathrm{3}}{S}_{\mathrm{*}}+q_{\mathrm{2}}{E}_{\mathrm{*}}+q_{\mathrm{1}}{I}_{\mathrm{*}}}{\mu +\phi },\\ V_{\mathrm{*}}=\frac{m\omega S_{\mathrm{*}}+\phi Q_{\mathrm{*}}+n\gamma I_{\mathrm{*}}}{\mu V_{\mathrm{*}}},\end{array}\right.

I_{\mathrm{*}}=\frac{\sigma \beta b-\left(m\omega +\mu +q_{\mathrm{3}}\right)A}{A\left[\beta +\alpha \left(m\omega +\mu +q_{\mathrm{3}}\right)\right]}=

\frac{\sigma \beta b/\left[\left(m\omega +\mu +q_{\mathrm{3}}\right)A\right]-\mathrm{1}}{\left[\beta +\alpha \left(m\omega +\mu +q_{\mathrm{3}}\right)\right]/\left(m\omega +\mu +q_{\mathrm{3}}\right)}=

\frac{R_{\mathrm{0}}-\mathrm{1}}{\left[\beta +\alpha \left(m\omega +\mu +q_{\mathrm{3}}\right)\right]/\left(m\omega +\mu +q_{\mathrm{3}}\right)}，

A=\left|V\right|=

\left(\sigma +\mu +q_{\mathrm{2}}\right)\left(\mu +\gamma +q_{\mathrm{1}}\right)-\sigma \left(\mathrm{1}-n\right)\gamma \mathrm{。}

当R_{\mathrm{0}}>\mathrm{1}时， I_{\mathrm{*}}>\mathrm{0}， 系统（1）在\Omega内存在唯一的正平衡点E^{\mathrm{*}}\left(S_{\mathrm{*}},E_{\mathrm{*}},I_{\mathrm{*}},Q_{\mathrm{*}},V_{\mathrm{*}}\right)。

3 平衡点的稳定性

定理 2　当R_{\mathrm{0}}\le \mathrm{1}时， E^{\mathrm{0}}是局部渐近稳定的； 当R_{\mathrm{0}}>\mathrm{1}时， E^{\mathrm{0}}是不稳定的。

证明　系统（1）在无蠕虫平衡点E^{\mathrm{0}}处的Jacobian矩阵为

\mathit{J}\left(E^{\mathrm{0}}\right)=\left[\begin{array}{ccccc}-\left(m\omega +\mu +q_{\mathrm{3}}\right)& \mathrm{0}& -\frac{\beta b}{m\omega +\mu +q_{\mathrm{3}}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -\left(\sigma +\mu +q_{\mathrm{2}}\right)& \frac{\beta b}{m\omega +\mu +q_{\mathrm{3}}}+\left(\mathrm{1}-n\right)\gamma & \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \sigma & -\left(\mu +\gamma +q_{\mathrm{1}}\right)& \mathrm{0}& \mathrm{0}\\ q_{\mathrm{3}}& q_{\mathrm{2}}& q_{\mathrm{1}}& -\left(\mu +\phi \right)& \mathrm{0}\\ m\omega & \mathrm{0}& n\gamma & \phi & -\mu \end{array}\right]，

显然\mathit{J}\left(E^{\mathrm{0}}\right)有3个负特征根， 即

{\lambda }_{\mathrm{1}}=-\left(m\omega +\mu +q_{\mathrm{3}}\right)， {\lambda }_{\mathrm{2}}=-\left(\mu +\phi \right)， {\lambda }_{\mathrm{3}}=-\mu，

其它特征值为

{\lambda }^{\mathrm{2}}+\left(\mathrm{2}\mu +\sigma +\gamma +q_{\mathrm{1}}+q_{\mathrm{2}}\right)\lambda +

\left(\sigma +\mu +q_{\mathrm{2}}\right)\left(\mu +\gamma +q_{\mathrm{1}}\right)-

\sigma \left[\frac{\beta b}{m\omega +\mu +q_{\mathrm{3}}}+\left(\mathrm{1}-n\right)\gamma \right]=\mathrm{0}，

c=\left(\sigma +\mu +q_{\mathrm{2}}\right)\left(\mu +\gamma +q_{\mathrm{1}}\right)-

\sigma \left[\frac{\beta b}{m\omega +\mu +q_{\mathrm{3}}}+\left(\mathrm{1}-n\right)\gamma \right]\mathrm{。}

当R_{\mathrm{0}}\le \mathrm{1}时， c>\mathrm{0}， 则

{\lambda }_{\mathrm{4}}{\lambda }_{\mathrm{5}}=c>\mathrm{0}，

所以， {\lambda }_{\mathrm{4}}和{\lambda }_{\mathrm{5}}均为负数。\mathit{J}\left(E^{\mathrm{0}}\right)的特征值均为负数， 故无蠕虫平衡点E^{\mathrm{0}}局部渐近稳定。

当R_{\mathrm{0}}>\mathrm{1}时， ， 则， 故{\lambda }_{\mathrm{4}}和{\lambda }_{\mathrm{5}}一正一负， 无蠕虫平衡点E^{\mathrm{0}}不稳定。

定理 3　当R_{\mathrm{0}}\le \mathrm{1}时， E^{\mathrm{0}}全局渐近稳定。

证明　构造Lyapunov函数

H\left(E,I\right)=\sigma E+\left(\sigma +\mu +q_{\mathrm{2}}\right)I\mathrm{。}

计算H沿着系统（1）的导数， 可得

\frac{\mathrm{d}H}{\mathrm{d}t}=\sigma E^{\text{'}}+\left(\sigma +\mu +q_{\mathrm{2}}\right)I^{\text{'}}=

\sigma \frac{\beta SI}{\mathrm{1}+\alpha I}+\sigma \left(\mathrm{1}-n\right)\gamma I-\left(\sigma +\mu +q_{\mathrm{2}}\right)\left(\mu +\gamma +q_{\mathrm{1}}\right)I\le

\sigma \beta SI+\sigma \left(\mathrm{1}-n\right)\gamma I-\left(\sigma +\mu +q_{\mathrm{2}}\right)\left(\mu +\gamma +q_{\mathrm{1}}\right)I=\frac{\sigma \beta S_{\mathrm{0}}}{R_{\mathrm{0}}}\left(\frac{R_{\mathrm{0}}S}{S_{\mathrm{0}}}-\mathrm{1}\right)I，

当R_{\mathrm{0}}\le \mathrm{1}时， \frac{\mathrm{d}H}{\mathrm{d}t}\le \mathrm{0}， 并且\frac{\mathrm{d}H}{\mathrm{d}t}=\mathrm{0}当且仅当I=\mathrm{0}。由LaSalle不变集原理^［14］可知， 当R_{\mathrm{0}}\le \mathrm{1}时， E^{\mathrm{0}}在\Omega内是全局渐近稳定的。

定理 4　当R_{\mathrm{0}}>\mathrm{1}时， E^{\mathrm{*}}是局部渐近稳定的。

证明　系统（1）在E^{\mathrm{*}}处的Jacobian矩阵为

\mathit{J}\left(E^{\mathrm{*}}\right)=\left[\begin{array}{ccccc}-\frac{\beta I_{\mathrm{*}}}{\mathrm{1}+\alpha I_{\mathrm{*}}}-\left(m\omega +\mu +q_{\mathrm{3}}\right)& \mathrm{0}& -\frac{\beta S_{\mathrm{*}}}{{\left(\mathrm{1}+\alpha I_{\mathrm{*}}\right)}^{\mathrm{2}}}& \mathrm{0}& \mathrm{0}\\ \frac{\beta I_{\mathrm{*}}}{\mathrm{1}+\alpha I_{\mathrm{*}}}& -\left(\sigma +\mu +q_{\mathrm{2}}\right)& \frac{\beta S_{\mathrm{*}}}{{\left(\mathrm{1}+\alpha I_{\mathrm{*}}\right)}^{\mathrm{2}}}+\left(\mathrm{1}-n\right)\gamma & \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \sigma & -\left(\mu +\gamma +q_{\mathrm{1}}\right)& \mathrm{0}& \mathrm{0}\\ q_{\mathrm{3}}& q_{\mathrm{2}}& q_{\mathrm{1}}& -\left(\mu +\phi \right)& \mathrm{0}\\ m\omega & \mathrm{0}& n\gamma & \phi & -\mu \end{array}\right]\mathrm{。}

易知{\lambda }_{\mathrm{1}}=-\mu， {\lambda }_{\mathrm{2}}=-\left(\mu +\phi \right)是\mathit{J}\left(E^{\mathrm{*}}\right)的两个负特征根， 只需考虑{\lambda }^{\mathrm{3}}+B_{\mathrm{1}}{\lambda }^{\mathrm{2}}+B_{\mathrm{2}}\lambda +B_{\mathrm{3}}=\mathrm{0}的根， 其中

\begin{array}{c}{B}_{\mathrm{1}}=\frac{\beta I_{\mathrm{*}}}{\mathrm{1}+\alpha I_{\mathrm{*}}}+\mathrm{3}\mu +m\omega +\sigma +\gamma +\\ q_{\mathrm{1}}+q_{\mathrm{2}}+q_{\mathrm{3}}>\mathrm{0},\end{array}

B_{\mathrm{2}}=\left(\frac{\beta I_{\mathrm{*}}}{\mathrm{1}+\alpha I_{\mathrm{*}}}+m\omega +\mu +q_{\mathrm{3}}\right)\left(\mathrm{2}\mu +\sigma +\gamma +q_{\mathrm{1}}+q_{\mathrm{2}}\right)+A-\frac{A}{\mathrm{1}+\alpha I_{\mathrm{*}}}>\mathrm{0}，

\begin{array}{c}{B}_{\mathrm{3}}=\frac{\beta I_{\mathrm{*}}}{\mathrm{1}+\alpha I_{\mathrm{*}}}A+(m\omega +\mu +q_{\mathrm{3}})·\\ \left(A-\frac{A}{\mathrm{1}+\alpha I_{\mathrm{*}}}\right)>\mathrm{0}，\end{array}

B_{\mathrm{1}}{B}_{\mathrm{2}}-B_{\mathrm{3}}=

\frac{\beta I_{\mathrm{*}}}{\mathrm{1}+\alpha I_{\mathrm{*}}}\left(\frac{\beta I_{\mathrm{*}}}{\mathrm{1}+\alpha I_{\mathrm{*}}}+m\omega +\mu +q_{\mathrm{3}}\right)·

(\mathrm{2}\mu +\sigma +\gamma +q_{\mathrm{1}}+q_{\mathrm{2}})-\frac{\beta I_{\mathrm{*}}}{\mathrm{1}+\alpha I_{\mathrm{*}}}\frac{A}{\mathrm{1}+\alpha I_{\mathrm{*}}}+

(m\omega +\mu +q_{\mathrm{3}})\left(\frac{\beta I_{\mathrm{*}}}{\mathrm{1}+\alpha I_{\mathrm{*}}}+m\omega +\mu +q_{\mathrm{3}}\right)·

(\mathrm{2}\mu +\sigma +\gamma +q_{\mathrm{1}}+q_{\mathrm{2}})+

(\mathrm{2}\mu +\sigma +\gamma +q_{\mathrm{1}}+q_{\mathrm{2}})\left(\frac{\beta I_{\mathrm{*}}}{\mathrm{1}+\alpha I_{\mathrm{*}}}+m\omega +\mu +q_{\mathrm{3}}\right)·

(\mathrm{2}\mu +\sigma +\gamma +q_{\mathrm{1}}+q_{\mathrm{2}})+

(\mathrm{2}\mu +\sigma +\gamma +q_{\mathrm{1}}+q_{\mathrm{2}})\left(A-\frac{A}{\mathrm{1}+\alpha I_{\mathrm{*}}}\right)\mathrm{。}

由不等式

-\frac{\beta I_{\mathrm{*}}}{\mathrm{1}+\alpha I_{\mathrm{*}}}\frac{A}{\mathrm{1}+\alpha I_{\mathrm{*}}}+(\mathrm{2}\mu +\sigma +\gamma +q_{\mathrm{1}}+q_{\mathrm{2}})·

(\mathrm{2}\mu +\sigma +\gamma +q_{\mathrm{1}}+q_{\mathrm{2}})\frac{\beta I_{\mathrm{*}}}{\mathrm{1}+\alpha I_{\mathrm{*}}}=\frac{\beta I_{\mathrm{*}}}{\mathrm{1}+\alpha I_{\mathrm{*}}}·

\begin{array}{c}[(\mathrm{2}\mu +\sigma +\gamma +q_{\mathrm{1}}+q_{\mathrm{2}})(\mathrm{2}\mu +\sigma +\gamma +q_{\mathrm{1}}+q_{\mathrm{2}})-\\ ((\sigma +\mu +q_{\mathrm{2}})(\mu +\gamma +q_{\mathrm{1}})-\sigma (\mathrm{1}-n)\gamma )/\end{array}(\mathrm{1}+\alpha I_{\mathrm{*}})]>\mathrm{0},

可得B_{\mathrm{1}}{B}_{\mathrm{2}}-B_{\mathrm{3}}>\mathrm{0}。根据Hurwitz判据^［15］， H_{\mathrm{1}}=B_{\mathrm{1}}>\mathrm{0}， H_{\mathrm{2}}=B_{\mathrm{1}}{B}_{\mathrm{2}}-B_{\mathrm{3}}>\mathrm{0}， H_{\mathrm{3}}=B_{\mathrm{3}}{H}_{\mathrm{2}}>\mathrm{0}， \mathit{J}（E^{\mathrm{*}}）的其它3个特征根也是负的， 故E^{\mathrm{*}}是局部渐近稳定的。

引理 1^［16］　系统（1）在\Omega内是一致持续的， 等价于存在一个紧的吸引集K\subset \Omega。

定理 5　当R_{\mathrm{0}}>\mathrm{1}时， E^{\mathrm{*}}是全局渐近稳定的。

证明　用\mathrm{L}\mathrm{i}\mathrm{­}\mathrm{M}\mathrm{u}\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{n}\mathrm{e}\mathrm{y}几何方法^［12］来证明。据文献［17］的定理2.3， 当且仅当R_{\mathrm{0}}>\mathrm{1}时系统（1）在\Omega内一致持续， 故单连通集\Omega内存在一个紧的吸引集K。

考虑系统（1）的子系统

\left\{\begin{array}{l}\frac{\mathrm{d}S}{\mathrm{d}t}=b-\frac{\beta SI}{\mathrm{1}+\alpha I}-\left(m\omega +\mu +q_{\mathrm{3}}\right)S,\\ \frac{\mathrm{d}E}{\mathrm{d}t}=\frac{\beta SI}{\mathrm{1}+\alpha I}-\left(\sigma +\mu +q_{\mathrm{2}}\right)E+\left(\mathrm{1}-n\right)\gamma I,\\ \frac{\mathrm{d}I}{\mathrm{d}t}=\sigma E-\left(\mu +\gamma +q_{\mathrm{1}}\right)I,\end{array}\right.（2）

系统（2）在E^{\mathrm{*}}处的Jacobian矩阵为

\mathit{J}\left(E^{\mathrm{*}}\right)=\left[\begin{array}{ccc}-\frac{\beta SI}{\mathrm{1}+\alpha I}-\left(m\omega +\mu +q_{\mathrm{3}}\right)& \mathrm{0}& -\frac{\beta S}{{\left(\mathrm{1}+\alpha I\right)}^{\mathrm{2}}}\\ \frac{\beta SI}{\mathrm{1}+\alpha I}& -\left(\sigma +\mu +q_{\mathrm{2}}\right)& \frac{\beta S}{{\left(\mathrm{1}+\alpha I\right)}^{\mathrm{2}}}+\left(\mathrm{1}-n\right)\gamma \\ \mathrm{0}& \sigma & -\left(\mu +\gamma +q_{\mathrm{1}}\right)\end{array}\right]

它的第二加性复合矩阵^［18］为

{\mathit{J}}^{\left[\mathrm{2}\right]}\left(E^{\mathrm{*}}\right)=\left[\begin{array}{ccc}-\frac{\beta I}{\mathrm{1}+\alpha I}-M& \frac{\beta S}{{\left(\mathrm{1}+\alpha I\right)}^{\mathrm{2}}}+\left(\mathrm{1}-n\right)\gamma & \frac{\beta S}{{\left(\mathrm{1}+\alpha I\right)}^{\mathrm{2}}}\\ \sigma & -\frac{\beta I}{\mathrm{1}+\alpha I}-N& \mathrm{0}\\ \mathrm{0}& \frac{\beta I}{\mathrm{1}+\alpha I}& -K\end{array}\right]，

其中

M=\mathrm{2}\mu +\mathrm{m}\omega +\sigma +q_{\mathrm{2}}+q_{\mathrm{3}}，

N=\mathrm{2}\mu +m\omega +\gamma +q_{\mathrm{1}}+q_{\mathrm{3}}，

K=\mathrm{2}\mu +\sigma +\gamma +q_{\mathrm{1}}+q_{\mathrm{2}}，

选P(S,E,I)=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\mathrm{1},\frac{E}{I},\frac{E}{I}\right)， 则{\mathit{P}}_f{\mathit{P}}^{-\mathrm{1}}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\mathrm{0},\frac{E^{\text{'}}}{E}-\frac{I^{\text{'}}}{I},\frac{E^{\text{'}}}{E}-\frac{I^{\text{'}}}{I}\right)， 矩阵\mathit{B}={\mathit{P}}_f{\mathit{P}}^{-\mathrm{1}}+\mathit{P}{\mathit{J}}^{\left[\mathrm{2}\right]}{\mathit{P}}^{-\mathrm{1}}可写成分块矩阵

\mathit{B}=\left[\begin{array}{cc}{B}_{\mathrm{11}}& B_{\mathrm{12}}\\ B_{\mathrm{21}}& B_{\mathrm{22}}\end{array}\right]，{\mathit{B}}_{\mathrm{11}}=-\frac{\beta I}{\mathrm{1}+\alpha I}-M，

{\mathit{B}}_{\mathrm{12}}=\left(\frac{I}{E}\left(\frac{\beta S}{{(\mathrm{1}+\alpha I)}^{\mathrm{2}}}+(\mathrm{1}-n)\gamma \right),\frac{\beta SI}{E{(\mathrm{1}+\alpha I)}^{\mathrm{2}}}\right),

{\mathit{B}}_{\mathrm{21}}={\left(\frac{E}{I}\sigma ,\mathrm{0}\right)}^{\mathrm{T}}，

{\mathit{B}}_{\mathrm{22}}=\left[\begin{array}{cc}\frac{E^{\text{'}}}{E}-\frac{I^{\text{'}}}{I}-\frac{\beta I}{\mathrm{1}+\alpha I}-N& \mathrm{0}\\ \frac{\beta I}{\mathrm{1}+\alpha I}& \frac{E^{\text{'}}}{E}-\frac{I^{\text{'}}}{I}-K\end{array}\right]\mathrm{。}

令(u,v,w)为R^{\mathrm{3}}上的一个向量， 定义其范数‖\cdot ‖为

‖(u,v,w)‖=\mathrm{m}\mathrm{a}\mathrm{x}\left\{\left|u\right|,\left|v\right|+\left|w\right|\right\}，

关于范数的Lozinskil测度为\mu (\mathit{B})， \mu (\mathit{B})\le \mathrm{s}\mathrm{u}\mathrm{p}\left\{g_{\mathrm{1}},g_{\mathrm{2}}\right\}， 其中g_{\mathrm{1}}={\mu }_{\mathrm{1}}({\mathit{B}}_{\mathrm{11}})+\left|{\mathit{B}}_{\mathrm{12}}\right|， g_{\mathrm{12}}=\left|{\mathit{B}}_{\mathrm{21}}\right|+{\mu }_{\mathrm{1}}({\mathit{B}}_{\mathrm{22}})， \left|{\mathit{B}}_{\mathrm{12}}\right|， \left|{\mathit{B}}_{\mathrm{21}}\right|为l_{\mathrm{1}}矩阵范数， {\mu }_{\mathrm{1}}为l_{\mathrm{1}}范数的Lozinskil测度。

{\mu }_{\mathrm{1}}({\mathit{B}}_{\mathrm{11}})=-\frac{\beta I}{\mathrm{1}+\alpha I}-M， \left|{\mathit{B}}_{\mathrm{21}}\right|=\frac{\sigma E}{I}，

\left|{\mathit{B}}_{\mathrm{12}}\right|=\mathrm{m}\mathrm{a}\mathrm{x}\left\{\frac{I}{E}(\frac{\beta S}{{(\mathrm{1}+\alpha I)}^{\mathrm{2}}}+(\mathrm{1}-n)\gamma ),\right.

\begin{array}{c}\left.\frac{\beta SI}{E{(\mathrm{1}+\alpha I)}^{\mathrm{2}}}\right\}=\frac{I}{E}\left[\frac{\beta S}{{(\mathrm{1}+\alpha I)}^{\mathrm{2}}}+(\mathrm{1}-n)\gamma \right]=\\ \frac{\beta SI}{E{(\mathrm{1}+\alpha I)}^{\mathrm{2}}}+\frac{I(\mathrm{1}-n)\gamma }{E}，\end{array}

{\mathit{B}}_{\mathrm{22}}^{\text{'}}=\left[\begin{array}{cc}\frac{E^{\text{'}}}{E}-\frac{I^{\text{'}}}{I}-N& \mathrm{0}\\ \frac{\beta S}{\mathrm{1}+\alpha I}& \frac{E^{\text{'}}}{E}-\frac{I^{\text{'}}}{I}-K\end{array}\right]，

取{\mathit{B}}_{\mathrm{22}}^{\text{'}}的对角元素的最大值即得{\mu }_{\mathrm{1}}(B_{\mathrm{22}})。

{\mu }_{\mathrm{1}}\left({\mathit{B}}_{\mathrm{22}}\right)=\mathrm{m}\mathrm{a}\mathrm{x}\left\{\frac{E^{\text{'}}}{E}-\frac{I^{\text{'}}}{I}-N,\frac{E^{\text{'}}}{E}-\frac{I^{\text{'}}}{I}-K\right\}=

\frac{E^{\text{'}}}{E}-\frac{I^{\text{'}}}{I}-\mathrm{m}\mathrm{i}\mathrm{n}\left\{N,K\right\}，

g_{\mathrm{1}}=-\frac{\beta I}{\mathrm{1}+\alpha I}-M+\frac{\beta SI}{E{(\mathrm{1}+\alpha I)}^{\mathrm{2}}}+\frac{I}{E}(\mathrm{1}-n)\gamma ，

g_{\mathrm{2}}=\frac{\sigma E}{I}+\frac{E^{\text{'}}}{E}-\frac{I^{\text{'}}}{I}-\mathrm{m}\mathrm{i}\mathrm{n}\left\{N,K\right\}，

\frac{E^{\text{'}}}{E}=\frac{\beta SI}{E(\mathrm{1}+\alpha I)}-(\sigma +\mu +{\mathrm{q}}_{\mathrm{2}})+\frac{(\mathrm{1}-n)\gamma I}{E}，

\frac{I^{\text{'}}}{I}=\frac{\sigma E}{I}-(\mu +\gamma +q_{\mathrm{1}})，

\begin{array}{c}{g}_{\mathrm{1}}=-\frac{\beta SI}{\mathrm{1}+\alpha I}-M+\frac{\beta SI}{E(\mathrm{1}+\alpha I)}+\frac{I(\mathrm{1}-n)\gamma }{E}+\\ \left(\frac{\beta SI}{E{(\mathrm{1}+\alpha I)}^{\mathrm{2}}}-\frac{\beta SI}{E(\mathrm{1}+\alpha I)}\right)\le \\ -\frac{\beta SI}{\mathrm{1}+\alpha I}-(\mu +m\omega +q_{\mathrm{3}})+\frac{E^{\text{'}}}{E}，\end{array}

\begin{array}{c}{g}_{\mathrm{2}}=\frac{I^{\text{'}}}{I}+(\mu +\gamma +q_{\mathrm{1}})+\frac{E^{\text{'}}}{E}-\frac{I^{\text{'}}}{I}-\mathrm{m}\mathrm{i}\mathrm{n}\left\{N,K\right\}\le \\ \frac{E^{\text{'}}}{E}-\mu ，\end{array}

\begin{array}{c}\mu \left(\mathit{B}\right)\le \mathrm{s}\mathrm{u}\mathrm{p}\left\{g_{\mathrm{1}},g_{\mathrm{2}}\right\}\le \\ \mathrm{s}\mathrm{u}\mathrm{p}\left\{-\frac{\beta SI}{\mathrm{1}+\alpha I}-\mu +\frac{E^{\text{'}}}{E},\frac{E^{\text{'}}}{E}-\mu \right\}\le \frac{E^{\text{'}}}{E}-\mu \mathrm{。}\end{array}

从而

\frac{\mathrm{1}}{t}{\int }_{\mathrm{0}}^t\mu \left(\mathit{B}\right)ds\le \frac{\mathrm{1}}{t}{\int }_{\mathrm{0}}^t\left(\frac{E^{\text{'}}}{E}-\mu \right)\mathrm{d}s=\frac{\mathrm{1}}{t}\mathrm{l}\mathrm{n}\frac{E\left(t\right)}{E\left(\mathrm{0}\right)}-\mu ，

根据文献［12］中定理3.5， (S_{\mathrm{*}},E_{\mathrm{*}},I_{\mathrm{*}})是全局渐近稳定的。

由系统（1）的子系统

\left\{\begin{array}{l}\frac{\mathrm{d}Q}{\mathrm{d}t}=q_{\mathrm{3}}S+q_{\mathrm{2}}E+q_{\mathrm{1}}I-\left(\mu +\phi \right)Q,\\ \frac{\mathrm{d}V}{\mathrm{d}t}=m\omega S+n\gamma I+\phi Q-\mu V,\end{array}\right.

可得

\left\{\begin{array}{l}\frac{\mathrm{d}Q}{\mathrm{d}t}=q_{\mathrm{3}}{S}_{\mathrm{*}}+q_{\mathrm{2}}{E}_{\mathrm{*}}+q_{\mathrm{1}}{I}_{\mathrm{*}}-\left(\mu +\phi \right)Q_{\mathrm{*}},\\ \frac{\mathrm{d}V}{\mathrm{d}t}=m\omega S_{\mathrm{*}}+n\gamma I_{\mathrm{*}}+\phi Q_{\mathrm{*}}-\mu V_{\mathrm{*}},\end{array}\right.

Q(t)=Q(\mathrm{0}){\mathrm{e}}^{-(\mu +\phi )t}+\frac{q_{\mathrm{3}}{S}_{\mathrm{*}}+q_{\mathrm{2}}{E}_{\mathrm{*}}+q_{\mathrm{1}}{I}_{\mathrm{*}}}{\mu +\phi }，

V(t)=V(\mathrm{0}){\mathrm{e}}^{-\mu t}+\frac{m\omega S_{\mathrm{*}}+n\gamma I_{\mathrm{*}}+\phi Q_{\mathrm{*}}}{\mu }，

取极限， 可得

Q(t)\to Q_{\mathrm{*}}， V(t)\to V_{\mathrm{*}}， t\to \mathrm{\infty },

所以， 正平衡点E^{\mathrm{*}}是全局渐近稳定的。

4 数值模拟和敏感性分析

为验证以上结果， 用Matlab对系统（1）的结论进行数值模拟， 具体情况如表 1 所示。

图 2 模拟了无蠕虫平衡点E^{\mathrm{0}}的稳定性， 此时基本再生数， 可看出t趋于无穷时， 潜伏设备和感染设备均为零， 恶意蠕虫传播结束， E^{\mathrm{0}}是全局渐近稳定的， 验证了定理3。

图 3 中， 令参数\omega =\mathrm{0.10}， 其余参数和图 2 一样， 此时， 基本再生数R_{\mathrm{0}}=\mathrm{1.651}\mathrm{8}>\mathrm{1}， 可看出t趋于无穷时， 潜伏设备和感染设备都大于零， 恶意蠕虫会持续存在， E^{\mathrm{*}}是全局渐近稳定的， 验证了定理5。

基本再生数R_{\mathrm{0}}和其它参数的关系如图 4。其中， \sigma， \beta， b与R_{\mathrm{0}}呈正相关， m， \omega， \gamma， q_{\mathrm{1}}， q_{\mathrm{2}}， q_{\mathrm{3}}， 与R_{\mathrm{0}}呈负相关， 且m， \omega， q_{\mathrm{1}}的变化对R_{\mathrm{0}}更敏感。结果表明， 增大感染设备的隔离率q_{\mathrm{1}}， 增大杀毒软件的覆盖率\omega以及生效率m都可使R_{\mathrm{0}}减小， 即通过减少感染设备与其他设备的连接， 提高人们对杀毒软件的认识， 可有效控制蠕虫病毒传播。

5 结　论

本文研究了一类具有饱和发生率的SEIQV网络蠕虫传播模型， 考虑了安装杀毒软件和隔离对蠕虫传播的影响。计算了基本再生数R_{\mathrm{0}}， 并证明了正平衡点的存在性， 通过LaSalle不变集原理， 证明了当R_{\mathrm{0}}\le \mathrm{1}时无蠕虫平衡点E^{\mathrm{0}}是全局渐近稳定的， 蠕虫病毒传播结束； 用\mathrm{L}\mathrm{i}\mathrm{­}\mathrm{M}\mathrm{u}\mathrm{l}\mathrm{o}\mathrm{w}\mathrm{n}\mathrm{e}\mathrm{y}几何方法证明了当R_{\mathrm{0}}>\mathrm{1}时正平衡点E^{\mathrm{*}}是全局渐近稳定的， 蠕虫病毒将会在网络中传播。通过Matlab验证了结论的准确性， 并发现通过增大感染设备的隔离率q_{\mathrm{1}}， 增大杀毒软件的覆盖率\omega以及生效率m， 都可减少R_{\mathrm{0}}， 即通过减少感染设备与其他设备的连接， 提高人们对杀毒软件的认识， 可有效控制蠕虫病毒传播。

参考文献

[1] COHEN F. Computer viruses： theory and experiments［J］. Computers & Security， 1987， 6（1）： 22-35.

[2] KEPHART J O， WHITE S R. Directed-graph epidemiol-ogical models of computer viruses［C］//Proceedings of 1991 IEEE Computer Society Symposium on Research in Security and Privacy， 1991： 343-359.

[3] KEPHART J O， WHITE S R. Measuring and modeling computer virus prevalence［C］//Proceedings 1993 IEEE Computer Society Symposium on Research in Security and Privacy. IEEE， 1993： 2-15.

[4] BADSHAH Q. Global stability of SEIQRS computer virus propagation model with non-linear incidence function［J］. Applied Mathematics， 2015， 6（11）： 1926-1938.

[5] XIAO X， FU P， DOU C， et al. Design and analysis of SEIQR worm propagation model in mobile internet［J］. Communications in Nonlinear Science and Numerical Simulation， 2017， 43： 341-350.

[6] KIM K S， IBRAHIM M M， JUNG I H， et al. Mathematical analysis of the effectiveness of control strategies to prevent the autorun virus transmission propagation［J］. Applied Mathematics and Computation， 2020， 371： 124955.

[7] GAO Q， ZHUANG J. Stability analysis and control strategies for worm attack in mobile networks via a VEIQS propagation model［J］. Applied Mathematics and Computation， 2020， 368： 124584.

[8] WANG R， XUE Y， XUE K. Dynamic analysis and optimal control of worm propagation model with saturated incidence rate［J］. Frontiers in Physics， 2023， 10： 1098040.

[9] MADHUSUDANAN V， SRINIVAS M N， NWOKOYE C N H， et al. HOPF-bifurcation analysis of delayed computer virus model with holling type iii incidence function and treatment［J］. Scientific African， 2022， 15： e01125.

[10] WANG R， XUE Y. Stability analysis and optimal control of worm propagation model with saturated incidence rate［J］. Computers & Security， 2023， 125： 103063.

[11] 王瑞玲， 薛亚奎. 具有饱和发生率的网络蠕虫病毒模型分析［J］. 华中师范大学学报（自然科学版）， 2023， 57（3）： 347-353.

[12] WANG Ruiling， XUE Yakui. Analysis of network worm virus model with saturated incidence rate［J］. Journal of Central China Normal University（Natural Sciences）， 2023， 57（3）： 347-353. （in Chinese）

[13] LI M Y， MULDOWNEY J S.A geometric approach to global-stability problems［J］. SIAM Journal on Mathematical Analysis， 1996， 27（4）： 1070-1083.

[14] VAN DEN DRIESSCHE P， WATMOUGH J. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission［J］. Mathematical Biosciences， 2002， 180（1/2）： 29-48.

[15] LASALLE J P. The stability of dynamical systems［M］. Philadelphia： SIAM Press， 1976.

[16] 王拉娣.传染病动力学模型及控制策略研究［D］.上海： 上海大学， 2005.

[17] BUTLER G， WALTMAN P. Persistence in dynamical systems［J］. Journal of Differential Equations， 1986， 63（2）： 255-263.

[18] LIU X， YANG L. Stability analysis of an SEIQV epidemic model with saturated incidence rate［J］. Nonlinear Analysis： Real World Applications， 2012， 13（6）： 2671-2679.

[19] 任亚鑫， 薛亚奎.一类社交网络中受媒介影响的信息传播动力学模型分析与最优控制［J］.中北大学学报（自然科学版）， 2023， 44（1）： 1-7.

[20] REN Yaxin， XUE Yakui. Analysis and optimal control of a dynamic model for information propagation affected by medium in social networks［J］. Journal of North University of China（Natural Science Edition）， 2023， 44（1）： 1-7. （in Chinese）

基金资助

国家自然科学青年基金(11301491)

山西省自然科学青年基金资助项目(2018010221040)