基于动态隔离的工业互联网恶意软件传播模型研究

Research on malware propagation model of industrial internet based on dynamic isolation

摘要

为更真实地描述恶意软件在工业互联网中的传播动态行为，基于动力学系统理论，引入时延因素，构建可变隔离率的时延恶意软件传播模型。针对恶意软件感染早期阶段感染节点数量较小的特点，模型根据感染主机数量的变化动态调整隔离率；对模型进行数值分析和仿真实验，结合理论推导得出Hopf分岔的阈值。研究结果表明：时延的存在导致系统出现Hopf分岔现象，当时延小于阈值时，系统保持稳定；当时延超过阈值时，系统呈现不稳定状态。表明该模型能够更准确地预测恶意软件传播的动态过程。研究结论为工业互联网安全防护策略的制定提供理论依据。

Abstract

To more accurately describe the dynamic behavior of malware propagation in the industrial internet, this paper constructs a delay malware propagation model with a variable isolation rate, based on the theory of dynamical systems and incorporating delay factors. This model dynamically adjusts the isolation rate with the change in the number of infected hosts, addressing the characteristic of a small number of infected nodes in the early stages of malware infection. By performing numerical analysis and simulation experiments, and combining theoretical derivation, the Hopf bifurcation threshold is obtained. The research results show that the presence of delay leads to the occurrence of Hopf bifurcation. When the delay is below the threshold, the system remains stable; when the delay exceeds the threshold, the system becomes unstable. This proves that the model can more accurately predict the dynamic process of malware propagation. The research conclusions provide a theoretical basis for the formulation of security protection strategies in the industrial internet.

Graphical abstract

关键词

网络安全 / 工业互联网 / 恶意软件传播模型 / 动态隔离 / 时延系统 / Hopf分岔阈值

Key words

network security / industrial internet / malware propagation model / dynamic isolation / time-delay system / Hopf bifurcation threshold

0 引言

工业互联网作为工业数字化与智能化转型的重要基础^[1-4]，所面临的安全威胁日益严峻。网络攻击带来的生产中断和数据泄露等问题，对制造业构成巨大威胁^[5-8]。由于现有的恶意软件传播模型普遍存在固定隔离率、忽略时延、节点异质性考虑不足、网络拓扑结构过于简化等问题，如何利用传播动力学模型描述恶意软件在工业互联网中的传播过程，并结合工业互联网中的具体情况，设计具有针对性的防控策略，已经成为当前研究的重点问题。

针对现有恶意软件模型无法适用于所有工业中的细节等问题，学者们提出了SIS模型^[9-11]和SIR模型^[12-14]等用于研究。XIAO等^[15]研究恶意软件在网络理系统中的传播动力学特性及其临界点预测方法，为理解复杂系统中的恶意软件传播机制奠定了理论基础。REY等^[16]基于个体的数学模型，提出一种专门用于研究工业关键基础设施中的恶意软件传播行为。MAHBOUBI等^[17]对异构移动恶意软件进行了系统性分析，识别出影响传播过程的关键参数。上述研究虽然能够初步揭示恶意软件传播的基本特性，但未能充分考虑传播过程中参数的动态演化特征。工业互联网环境下的恶意软件传播往往伴随着参数的动态变化过程，单一的静态分析方法难以准确预测传播趋势。YAO等^[18]提出考虑时滞效应和可变感染率的蠕虫传播模型，探讨动态传播过程中变量的演化规律，通过引入可变感染率提升预测精度。YU等^[19]在此基础上进一步改进包含双感染率的SEI2RS模型，实现了对网络物理系统中双感染率恶意软件传播的建模，但该模型未考虑时延效应。YANG等^[20]提出一种双时延SIDQR模型，在文献[19]的基础上强调了时间延迟在传播过程中的作用，提升了动态传播模型的复杂度。综上可知，现有的模型虽然在一定程度上可以模拟恶意软件传播行为，但忽略了隔离率的动态变化在实际工业系统安全防护中的重要价值。

因此，本文基于动力学系统理论，结合时延因素和可变隔离率，构建一个更符合工业互联网场景的恶意软件传播模型。采用动力学系统理论描述恶意软件的传播过程，将节点状态转化为数学变量，分析其动态演化规律；引入时延因素，反映出工业互联网中存在的时间延迟情况；并通过数值分析和仿真实验，研究系统的稳定性条件，推导Hopf分岔阈值。结合理论推导和仿真数据实验，验证模型对恶意软件动态传播的预测能力。

1 模型设计

设工业互联网上主机的总数为N，SIQVD模型的状态转换见图1。

该模型中，易感主机在遭受恶意软件攻击后会转变为感染主机，模型采用动态隔离策略来控制恶意软件在易感主机间的传播。网络管理员可根据实时情况动态调整感染主机的隔离率，而被隔离的主机可以通过安装补丁等方式转变为恢复态。假设所有主机存在易感（S）、感染（I）、隔离（Q）、恢复（V）和延迟（D）5种状态。SIDQV模型为

\left\{\begin{array}{l}\mathrm{d}S/\mathrm{d}t=-\beta I\left(t\right)S\left(t\right)+\delta V\left(t\right)\\ \mathrm{d}I/\mathrm{d}t=\beta I\left(t\right)S\left(t\right)-\gamma I\left(t\right)-\alpha I\left(t\right)\\ \mathrm{d}Q/\mathrm{d}t=\alpha I\left(t\right)-\theta Q\left(t\right)\\ \mathrm{d}V/\mathrm{d}t=\theta Q\left(t\right)+\gamma I\left(t-\tau \right)-\delta V\left(t\right)\\ \mathrm{d}D/\mathrm{d}t=\gamma I\left(t\right)-\gamma I\left(t-\tau \right)\end{array}\right.，（1）

式中：t为时间，s；S为易感者的台数；I为感染者的台数；Q为被隔离者的台数；V为接受疫苗的个体台数；D为延迟状态的个体台数；\alpha为感染者的被隔离率；\beta为感染率；\gamma为感染者的恢复率；\theta为被隔离者向接受疫苗者的转化率；\delta为接受疫苗者失去免疫的速率；\tau为时间延迟。

假设\alpha为t的函数，记为

\alpha \left(t\right)={\alpha }_{\mathrm{0}}{f}_{\mathrm{1}}\left(I\left(t\right)\right)，

式中：{\alpha }_{\mathrm{0}}为初始感染率；f_{\mathrm{1}}\left(\cdot \right)为感染宿主I的非线性函数，表达式为

f_{\mathrm{1}}\left(I\left(t\right)\right)=\frac{I}{\mathrm{1}+aI}，（2）

满足：①f_{\mathrm{1}}\left(\mathrm{0}\right)=\mathrm{0}；②f_{\mathrm{1}}\text{'}\left(I\right)>\mathrm{0}；③；④。

采用微分方程组表示SIDQV模型为

\left\{\begin{array}{l}\mathrm{d}S/\mathrm{d}t=-\beta I\left(t\right)S\left(t\right)+\delta V\left(t\right)\\ \mathrm{d}I/\mathrm{d}t=\beta I\left(t\right)S\left(t\right)-\gamma I\left(t\right)-{\alpha }_{\mathrm{0}}f\left(I\left(t\right)\right)\\ \mathrm{d}Q/\mathrm{d}t={\alpha }_{\mathrm{0}}f\left(I\left(t\right)\right)-\theta Q\left(t\right)\\ \mathrm{d}V/\mathrm{d}t=\theta Q\left(t\right)+\gamma I\left(t-\tau \right)-\delta V\left(t\right)\\ \mathrm{d}D/\mathrm{d}t=\gamma I\left(t\right)-\gamma I\left(t-\tau \right)\end{array}\right.，（3）

式中， f\left(I\left(t\right)\right)=f_{\mathrm{1}}\left(I\left(t\right)\right)I\left(t\right)。

2 平衡稳定性与分岔分析

为分析系统平衡点的稳定性条件，探讨在时延参数变化下系统是否出现分岔现象，揭示系统的内在动态特性，进行如下分析。

为验证系统的平衡稳定性，需证明式（3）有且只有1个正平衡点。基本再生数R_{\mathrm{0}}可以衡量恶意软件传播的临界条件，当R_{\mathrm{0}}>1时，恶意软件会持续传播，而当R_{\mathrm{0}}<1时，恶意软件最终会消失。

根据式（3）可得

R_{\mathrm{0}}=\frac{\beta N}{\left(\gamma +\alpha f\mathrm{\text{'}}\left(\mathrm{0}\right)\right)}。

定理1 当R_{\mathrm{0}}\ge \mathrm{1}时，式（3）有1个独特正平衡点为E^{\mathrm{*}}=\left(S^{\mathrm{*}},I^{\mathrm{*}},D^{\mathrm{*}},Q^{\mathrm{*}},V^{\mathrm{*}}\right)，其中，

S^{\mathrm{*}}=\frac{{\alpha }_{\mathrm{0}}f\left(I^{\mathrm{*}}\right)+\gamma I^{\mathrm{*}}}{\beta I^{\mathrm{*}}}，I^{\mathrm{*}}=I^{\mathrm{*}}，Q^{\mathrm{*}}=\frac{{\alpha }_{\mathrm{0}}f\left(I^{\mathrm{*}}\right)}{\theta }，V^{\mathrm{*}}=\frac{{\alpha }_{\mathrm{0}}f\left(I\right)+\gamma I^{\mathrm{*}}}{\delta }，D^{\mathrm{*}}=\gamma \tau I^{\mathrm{*}}。

证明： 当式（3）稳定时，满足

\left\{\begin{array}{l}-\beta I\left(t\right)S\left(t\right)+\delta V\left(t\right)=\mathrm{0}\\ \beta I\left(t\right)S\left(t\right)-\gamma I\left(t\right)-{\alpha }_{\mathrm{0}}f\left(I\left(t\right)\right)=\mathrm{0}\\ {\alpha }_{\mathrm{0}}f\left(I\left(t\right)\right)-\theta Q\left(t\right)=\mathrm{0}\\ \theta Q\left(t\right)+\gamma I\left(t-\tau \right)-\delta V\left(t\right)=\mathrm{0}\\ \gamma I\left(t\right)-\gamma I\left(t-\tau \right)=\mathrm{0}\end{array}\right.，（4）

令I\left(t\right)>\mathrm{0}，则有

S^{\mathrm{*}}=\frac{{\alpha }_{\mathrm{0}}f\left(I^{\mathrm{*}}\right)+\gamma I^{\mathrm{*}}}{\beta I^{\mathrm{*}}},Q^{\mathrm{*}}=\frac{{\alpha }_{\mathrm{0}}f\left(I^{\mathrm{*}}\right)}{\theta }，

V^{\mathrm{*}}=\frac{{\alpha }_{\mathrm{0}}f\left(I\right)+\gamma I^{\mathrm{*}}}{\delta },D^{\mathrm{*}}=\gamma \tau I^{\mathrm{*}}\mathrm{。}

由于式（4）中的主机总数为二进制，系统平衡条件（即所有主机台数的总和始终保持恒定）为

H\left(I\right)=\frac{{\alpha }_{\mathrm{0}}{f}_{\mathrm{1}}\left(I\left(t\right)\right)+\gamma I}{\beta I}+I+\frac{{\alpha }_{\mathrm{0}}{f}_{\mathrm{1}}\text{'}\left(I\left(t\right)\right)}{\theta }+

\frac{{\alpha }_{\mathrm{0}}{f}_{\mathrm{1}}\text{'}\left(I\left(t\right)\right)+\gamma I}{\delta }-N=\mathrm{0}\mathrm{。}（5）

对式（5）求导得

H\mathrm{\text{'}}\left(I\right)=\frac{{\alpha }_{\mathrm{0}}{f}_{\mathrm{1}}\text{'}\left(I\left(t\right)\right)}{\theta }+\mathrm{1}+\frac{{\alpha }_{\mathrm{0}}{f}_{\mathrm{1}}\text{'}\left(I\left(t\right)\right)+\gamma }{\delta }+

\frac{{\alpha }_{\mathrm{0}}{f}_{\mathrm{1}}\text{'}\left(I\left(t\right)\right)+\gamma }{\beta I}-\frac{\gamma I+{\alpha }_{\mathrm{0}}{f}_{\mathrm{1}}\text{'}\left(I\left(t\right)\right)}{\beta I^{\mathrm{2}}}。（6）

将式（2）以及各参数代入式（6）得

H\text{'}\left(I\right)=\frac{\beta \delta \theta I^{\mathrm{2}}+\beta \gamma \theta I^{\mathrm{2}}+\mathrm{10}I\beta \gamma \theta }{\beta \theta \delta \left(I^{\mathrm{2}}+\mathrm{10}I+\mathrm{25}\right)}+

\frac{\mathrm{25}\alpha \beta \delta +\mathrm{25}\alpha \beta \theta -\mathrm{5}\alpha \delta \theta }{\beta \theta \delta \left(I^{\mathrm{2}}+\mathrm{10}I+\mathrm{25}\right)}+

\frac{\mathrm{25}\beta \delta \theta +\mathrm{25}\beta \gamma \theta }{\beta \theta \delta \left(I^{\mathrm{2}}+\mathrm{10}I+\mathrm{25}\right)}\mathrm{。}（7）

由于参数均为正值，故H\mathrm{\text{'}}\left(I\right)>\mathrm{0}，函数单调递增。若存在H\left(I\right)=\mathrm{0}的正根，则，基于此可得

\begin{array}{l}H\left(\mathrm{0}\right)=\underset{I\to \mathrm{0}}{\mathrm{l}\mathrm{i}\mathrm{m}}H\left(I\right)=-N+\frac{{\alpha }_{\mathrm{0}}{f}_{\mathrm{1}}\text{'}\left(\mathrm{0}\right)}{\theta }+\\ {}^{}{}^{}{}^{}{}^{}{}^{}{}^{}{}^{}{}^{}{}^{}{}^{}{}^{}\frac{{\alpha }_{\mathrm{0}}{f}_{\mathrm{1}}\text{'}\left(\mathrm{0}\right)}{\delta }+\frac{{\alpha }_{\mathrm{0}}{f}_{\mathrm{1}}\text{'}\left(\mathrm{0}\right)}{\beta }+\frac{\gamma }{\delta }+I=\\ {}^{}{}^{}{}^{}{}^{}{}^{}{}^{}{}^{}{}^{}{}^{}{}^{}{}^{}\frac{I\beta \delta \theta -N\beta \delta \theta +\alpha \beta \delta f_{\mathrm{1}}\mathrm{\text{'}}\left(I\right)}{\beta \delta \theta }+\end{array}

{}^{}{}^{}{}^{}{}^{}{}^{}{}^{}{}^{}{}^{}{}^{}{}^{}{}^{}\frac{\alpha \beta \theta f_{\mathrm{1}}\text{'}\left(I\right)+\alpha \theta \delta f_{\mathrm{1}}\text{'}\left(I\right)+\beta \gamma \theta }{\beta \delta \theta }。（8）

将参数代入式（8）可得，表明系统存在1个正平衡点。证毕。

由于Q\left(t\right)=N-S\left(t\right)-I\left(t\right)-V\left(t\right)-D\left(t\right)，将式（3）简化为

\left\{\begin{array}{l}\mathrm{d}S/\mathrm{d}t=-\beta I\left(t\right)S\left(t\right)+\delta V\left(t\right)\\ \mathrm{d}I/\mathrm{d}t=\beta I\left(t\right)S\left(t\right)-\gamma I\left(t\right)-{\alpha }_{\mathrm{0}}f\left(I\left(t\right)\right)\\ \mathrm{d}V/\mathrm{d}t=\theta \left(N-S\left(t\right)-I\left(t\right)-V\left(t\right)-D\left(t\right)\right)+\\ {}^{}\gamma I\left(t-\tau \right)-\delta V\left(t\right)\\ \mathrm{d}D/\mathrm{d}t=\gamma I\left(t\right)-\gamma I\left(t-\tau \right)\end{array}\right.。（9）

求式（9）关于平衡点E^{\mathrm{*}}的雅可比矩阵为\mathit{J}\left(E^{\mathrm{*}}\right)=\left(\begin{array}{cccc}-\beta I& -\beta S& \delta & \mathrm{0}\\ \beta I& \beta S-\gamma -{\alpha }_{\mathrm{0}}f\text{'}\left(I\right)& \mathrm{0}& \mathrm{0}\\ -\theta & -\theta +\gamma {\mathrm{e}}^{-\lambda \tau }& -\theta -\delta & -\theta \\ \mathrm{0}& \gamma -\gamma {\mathrm{e}}^{-\lambda \tau }& \mathrm{0}& \mathrm{0}\end{array}\right)。（10）

令\mathit{J}\left(E^{\mathrm{*}}\right)=\mathrm{0}，得到其特征方程为

P\left(\lambda \right)+B\left(\lambda \right){\mathrm{e}}^{-\lambda \tau }=\mathrm{0}，（11）

式中，

\left\{\begin{array}{l}P\left(\lambda \right)={\lambda }^{\mathrm{4}}+p_{\mathrm{3}}{\lambda }^{\mathrm{3}}+p_{\mathrm{2}}{\lambda }^{\mathrm{2}}+p_{\mathrm{1}}\lambda +p_{\mathrm{0}}\\ B\left(\lambda \right)=b_{\mathrm{1}}\lambda +b_{\mathrm{0}}\end{array}\right.。（12）

由式（10）可得

p_{\mathrm{3}}=I\beta -S\beta +\alpha f\text{'}\left(I\right)+\delta +\gamma +\theta，

p_{\mathrm{2}}=I\alpha \beta f\text{'}\left(I\right)+I\beta \delta +I\beta \gamma +I\beta \theta -

S\beta \delta -S\beta \theta +\alpha \delta f\text{'}\left(I\right)+\alpha f\text{'}\left(I\right)\theta +

\delta \gamma +\delta \theta +\gamma \theta，

p_{\mathrm{1}}=I\alpha \beta \delta f\mathrm{\text{'}}\left(I\right)+I\alpha \beta f\text{'}\left(I\right)\theta +I\beta \delta \gamma +

{}^{}I\beta \delta \theta +I\beta \gamma \theta -S\beta \delta \theta +\alpha \delta f\text{'}\left(I\right)\theta +\delta \gamma \theta

p_{\mathrm{0}}=I\beta \delta \gamma \theta，b_{\mathrm{1}}=-I\beta \delta \gamma，b_{\mathrm{0}}=-I\beta \delta \gamma \theta。

定理2 如果模型满足条件\mathrm{A}\left(\mathrm{1}\right)：

①p_{\mathrm{3}}>\mathrm{0}；②p_{\mathrm{3}}{p}_{\mathrm{2}}-\left(p_{\mathrm{1}}+b_{\mathrm{1}}\right)>\mathrm{0}；③\left(p_{\mathrm{1}}+b_{\mathrm{1}}\right)>\mathrm{0}，

则正平衡点E^{\mathrm{*}}是稳定且局部范围内没有时延的。

证明：当\tau =\mathrm{0}时，式（9）可化简为

{\lambda }^{\mathrm{4}}+p_{\mathrm{3}}{\lambda }^{\mathrm{3}}+p_{\mathrm{2}}{\lambda }^{\mathrm{2}}+p_{\mathrm{1}}\lambda +p_{\mathrm{0}}+\left(b_{\mathrm{1}}\lambda +b_{\mathrm{0}}\right){\mathrm{e}}^{-\lambda \tau }=\mathrm{0}\mathrm{。}（13）

根据Routh-Hurwitz准则，式（13）所有根都有负实部。因此，正平衡点E^{*}是局部无时滞渐近稳定的。证毕。

令\lambda =\mathrm{i}\omega，代入式（11），分离实部和虚部，其中，实部为

{\omega }^{\mathrm{4}}-p_{\mathrm{2}}{\omega }^{\mathrm{2}}+p_{\mathrm{0}}+b_{\mathrm{0}}\mathrm{c}\mathrm{o}\mathrm{s}\left(\omega \tau \right)+b_{\mathrm{1}}\omega \mathrm{s}\mathrm{i}\mathrm{n}\left(\omega \tau \right)=\mathrm{0}，（14）

虚部为

-p_{\mathrm{3}}{\omega }^{\mathrm{3}}+p_{\mathrm{1}}\omega -b_{\mathrm{0}}\mathrm{s}\mathrm{i}\mathrm{n}\left(\omega \tau \right)+{}^{}{}^{}{b}_{\mathrm{1}}\omega \mathrm{c}\mathrm{o}\mathrm{s}\left(\omega \tau \right)=\mathrm{0}\mathrm{。}（15）

将式（14）和式（15）平方后相加得

{b_{\mathrm{0}}}^{\mathrm{2}}+{\left(b_{\mathrm{1}}\omega \right)}^{\mathrm{2}}={\left({\omega }^{\mathrm{4}}-p_{\mathrm{2}}{\omega }^{\mathrm{2}}+p_{\mathrm{0}}\right)}^{\mathrm{2}}+

{}^{}{}^{}{}^{}{}^{}{}^{}{}^{}{\left(-p_{\mathrm{3}}{\omega }^{\mathrm{3}}+p_{\mathrm{1}}\omega \right)}^{\mathrm{2}}\mathrm{。}（16）

化简式（16）得

{\omega }^{\mathrm{8}}+D_{\mathrm{3}}{\omega }^{\mathrm{6}}+D_{\mathrm{2}}{\omega }^{\mathrm{4}}+D_{\mathrm{1}}{\omega }^{\mathrm{2}}+D_{\mathrm{0}}=\mathrm{0}，（17）

式中：D_{\mathrm{3}}={p_{\mathrm{3}}}^{\mathrm{2}}-\mathrm{2}{p}_{\mathrm{2}}；D_{\mathrm{2}}=\mathrm{2}{p}_{\mathrm{0}}+{p_{\mathrm{2}}}^{\mathrm{2}}-\mathrm{2}{p}_{\mathrm{1}}{p}_{\mathrm{3}}；D_{\mathrm{1}}=p_{\mathrm{1}}-\mathrm{2}{p}_{\mathrm{0}}{p}_{\mathrm{2}}-{b_{\mathrm{1}}}^{\mathrm{2}}；D_{\mathrm{0}}={p_{\mathrm{0}}}^{\mathrm{2}}-{b_{\mathrm{0}}}^{\mathrm{2}}。

令{\omega }^{\mathrm{2}}=z，将式（17）简化为z的多项式，得

h\left(z\right)=z^{\mathrm{4}}+D_{\mathrm{3}}{z}^{\mathrm{3}}+D_{\mathrm{2}}{z}^{\mathrm{2}}+D_{\mathrm{1}}z+D_{\mathrm{0}}。（18）

定理3 假设式（18）满足\mathrm{A}\left(\mathrm{2}\right)：①D_{\mathrm{0}}\ge \mathrm{0}，\Delta \ge \mathrm{0}，z_{\mathrm{1}}>\mathrm{0}或h\left(z_{\mathrm{1}}\right)>\mathrm{0}；②D_{\mathrm{0}}\ge \mathrm{0}，；③不存在z^{\mathrm{*}}\in \left(z_{\mathrm{1}},z_{\mathrm{2}},z_{\mathrm{3}}\right)使z^{\mathrm{*}}>\mathrm{0}和h\left(z^{\mathrm{*}}\right)\le \mathrm{0}，则系统的正平衡点E^{\mathrm{*}}是绝对稳定的。也就是说E^{\mathrm{*}}对任何\tau \ge \mathrm{0}都是渐近稳定的。

定理3中\Delta为特征多项式（18）的判别式，用于判断根的实虚性质；z_{\mathrm{1}}、z_{\mathrm{2}}、z_{\mathrm{3}}表示特征方程的实数根。

证明：假设式（18）中的模型满足条件\mathrm{A}\left(\mathrm{2}\right)，确保式（11）的根满足稳定性条件。

由前文推导可知，式（18）至少有1个正实根{\omega }_{\mathrm{0}}，式（11）有1对纯虚根\pm \mathrm{i}{\omega }_{\mathrm{0}}。这意味着系统存在周期性振荡。将纯虚根\pm \mathrm{i}{\omega }_{\mathrm{0}}代入式（14）和式（15）进行联立，可得时延阈值{\tau }_k>\mathrm{0}，使系统满足稳定性条件。

{\tau }_k=\frac{\mathrm{1}}{{\omega }_{\mathrm{0}}}\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}\left[\frac{b_{\mathrm{0}}\left({\omega }^{\mathrm{4}}-p_{\mathrm{2}}{\omega }^{\mathrm{2}}+p_{\mathrm{0}}\right)}{{b_{\mathrm{0}}}^{\mathrm{2}}+{\left(b_{\mathrm{1}}\omega \right)}^{\mathrm{2}}}+\right.

\left.\frac{\left(p_{\mathrm{3}}{\omega }^{\mathrm{3}}-p_{\mathrm{1}}\omega \right)b_{\mathrm{1}}\omega }{{b_{\mathrm{0}}}^{\mathrm{2}}+{\left(b_{\mathrm{1}}\omega \right)}^{\mathrm{2}}}+\frac{\mathrm{2}k\mathrm{\pi }}{{\omega }_{\mathrm{0}}}\right]。（19）

令\lambda \left(\tau \right)=v\left(\tau \right)+\mathrm{i}\omega \left(\tau \right)为式（11）的根，当\tau ={\tau }_k时，满足v\left({\tau }_k\right)=\mathrm{0}，\omega \left({\tau }_k\right)={\omega }_{\mathrm{0}}。

由上述推导可知，式（11）所示特征方程在一定的范围内存在1对纯虚根，可以证明存在1个时延使系统保持一定的稳定性。证毕。

定理4 设h\mathrm{\text{'}}\left(z_{\mathrm{0}}\right)\ne \mathrm{0}。若\tau ={\tau }_{\mathrm{0}}，则\pm \mathrm{i}{\omega }_{\mathrm{0}}为式（11）所示特征方程的一对纯虚根，且满足条件\mathrm{A}\left(\mathrm{2}\right)，则

\frac{\mathrm{d}\mathrm{R}\mathrm{e}\lambda \left({\tau }_{\mathrm{0}}\right)}{\mathrm{d}\tau }>\mathrm{0}。（20）

因此，当时延\tau >{\tau }_k时，式（11）存在至少1个实部为正的特征根,此时系统不稳定。对式（11）两边求导，结合b_{\mathrm{1}}和b_{\mathrm{0}}的相关项，可得

\begin{array}{l}{\left(\frac{\mathrm{d}\lambda }{\mathrm{d}\tau }\right)}^{-\mathrm{1}}=\frac{\mathrm{4}{\lambda }^{\mathrm{3}}+\mathrm{3}{p}_{\mathrm{3}}{\lambda }^{\mathrm{2}}+\mathrm{2}{p}_{\mathrm{2}}\lambda }{\left(b_{\mathrm{1}}\lambda +b_{\mathrm{0}}\right)\lambda {\mathrm{e}}^{-\lambda \tau }}+\\ {}^{}{}^{}{}^{}{}^{}{}^{}{}^{}{}^{}{}^{}{}^{}{}^{}{}^{}{}^{}{}^{}\frac{p_{\mathrm{1}}+b_{\mathrm{1}}{\mathrm{e}}^{-\lambda \tau }-\left(b_{\mathrm{1}}\lambda +b_{\mathrm{0}}\right)\tau {\mathrm{e}}^{-\lambda \tau }}{\left(b_{\mathrm{1}}\lambda +b_{\mathrm{0}}\right)\lambda e^{-\lambda \tau }}=\\ {}^{}{}^{}{}^{}{}^{}{}^{}{}^{}{}^{}{}^{}{}^{}{}^{}{}^{}{}^{}{}^{}\frac{\left(\mathrm{4}{\lambda }^{\mathrm{3}}+\mathrm{3}{p}_{\mathrm{3}}{\lambda }^{\mathrm{2}}+\mathrm{2}{p}_{\mathrm{2}}\lambda +p_{\mathrm{1}}\right){\mathrm{e}}^{\lambda \tau }}{\left(b_{\mathrm{1}}\lambda +b_{\mathrm{0}}\right)\lambda }+\end{array}

\frac{b_{\mathrm{1}}}{\left(b_{\mathrm{1}}\lambda +b_{\mathrm{0}}\right)\lambda }-\frac{\tau }{\lambda }。（21）

用式（14）和式（15）简化式（11）的实部导数部分，可得

\begin{array}{l}\mathrm{s}\mathrm{g}\mathrm{n}{\left[\frac{\mathrm{d}\mathrm{R}\mathrm{e}\lambda }{\mathrm{d}\tau }\right]}_{\tau ={\tau }_k}=\mathrm{s}\mathrm{g}\mathrm{n}{\left[\mathrm{R}\mathrm{e}\left(\frac{\mathrm{d}\lambda }{\mathrm{d}\tau }\right)\right]}_{\lambda =\mathrm{i}{\omega }_{\mathrm{0}}}^{-\mathrm{1}}=\\ \mathrm{s}\mathrm{g}\mathrm{n}\frac{{{\omega }_{\mathrm{0}}}^{\mathrm{2}}}{\Gamma }\left(\mathrm{4}{{\omega }_{\mathrm{0}}}^{\mathrm{6}}+\mathrm{3}{D}_{\mathrm{3}}{{\omega }_{\mathrm{0}}}^{\mathrm{4}}+\mathrm{2}{D}_{\mathrm{2}}{{\omega }_{\mathrm{0}}}^{\mathrm{2}}+D_{\mathrm{1}}\right)=\end{array}

\mathrm{s}\mathrm{g}\mathrm{n}\frac{{{\omega }_{\mathrm{0}}}^{\mathrm{2}}}{\Gamma }\left\{h\text{'}\left({{\omega }_{\mathrm{0}}}^{\mathrm{2}}\right)\right\}=\mathrm{s}\mathrm{g}\mathrm{n}\left\{h\text{'}\left({{\omega }_{\mathrm{0}}}^{\mathrm{2}}\right)\right\}，（22）

式中，\Gamma ={\left(b_{\mathrm{0}}{\omega }_{\mathrm{0}}\right)}^{\mathrm{2}}+{\left(b_{\mathrm{1}}{{\omega }_{\mathrm{0}}}^{\mathrm{2}}\right)}^{\mathrm{2}}。

遵循假设\mathrm{A}\left(\mathrm{2}\right)和假设h\text{'}\left({{\omega }_{\mathrm{0}}}^{\mathrm{2}}\right)\ne \mathrm{0}，可得

\frac{\mathrm{d}\left(\mathrm{R}\mathrm{e}\lambda \right)}{\mathrm{d}\tau }{\left|\right.}_{\tau ={\tau }_k}>\mathrm{0}。（23）

根据Routh-Hurwitz准则，式（11）的根在虚轴上从左向右交叉。当一个参数连续变化，从小于某值到大于某值时，系统变得不稳定。因此，根据泛函微分方程的Hopf分岔定理，横向条件成立，满足Hopf分岔的条件。证毕。

定理5 假设满足条件\mathrm{A}\left(\mathrm{1}\right)和\mathrm{A}\left(\mathrm{2}\right)：

①当，系统是稳定的，当\tau >{\tau }_{\mathrm{0}}，系统是不稳定的；②当式（3）满足\mathrm{A}\left(\mathrm{2}\right)时，在\tau ={\tau }_{\mathrm{0}}时发生Hopf分岔。

由此可得：当时，系统稳定，可以通过实施相关的方法和策略来抑制恶意软件的传播；当\tau >{\tau }_{\mathrm{0}}时，系统不稳定，恶意软件的传播将难以得到有效的控制。

3 数值模拟与实验

为验证模型的有效性和理论分析的准确性，针对和R_{\mathrm{0}}>\mathrm{1}两种情形开展数值仿真实验，分析模型的动态行为。从定性与定量两个维度对理论推导结果进行验证，探讨各参数对恶意软件传播过程的影响，同时对比分析理论预期与实际仿真结果。

3.1 时的实验分析

当时，设主机总数为4×10⁵台。感染率β为0.83×10^-6，传染宿主的恢复率\gamma为0.25，被隔离主机的免疫率\theta为0.08，接受补丁的主机以\delta为0.17的速率丧失免疫力，感染主机的隔离率\alpha为0.1。受隔离策略影响的模型中，隔离率随实际感染主机数的变化函数为

f\left(I\left(t\right)\right)=I\left(t\right)/\left(\mathrm{1}+aI\left(t\right)\right)，

式中，a代表隔离率对受感染主机数量的敏感性。当a=0时，隔离率为常数。设初始感染主机为10⁵台，其他状态均为易感状态。为验证第二部分中推导的关于时延结论的正确性，进行如下数值实验。

当τ=20<τ₀、a=10^-5时，5种主机数量的变化情况见图2（a）。由图2（a）可知，感染主机的数量持续降低，当t=60 s时，感染主机的数量为0。可知，在时延较短时，防御策略可迅速响应感染情况并快速实施隔离策略，从而在传染过程中形成正反馈的控制效应，使系统具有良好的稳定性。表明当时，工业互联网中恶意软件的传播可通过及时动态的隔离机制得到遏制。

当τ=70＞τ₀、a=10^-5时，5种主机数量的变化情况见图2（b）。由图2（b）可见，5种主机数量变化波动较大，当时延较大时，防御措施的响应滞后，难以在感染初期及时隔离被感染的主机，导致恶意软件在工业互联网的局部区域反复传播。在实际工业互联网环境中，防御系统的时延若超出阈值范围时，防护措施不能及时覆盖所有受感染节点，出现局部感染暴发的风险，此时就需要引入更加智能化和自适应的防护策略，补偿时延所带来的不利影响。

通过仿真实验来验证理论分析和数值结果的一致性和正确性。其中，在τ=20<τ₀时，易感主机、感染主机、隔离主机和恢复主机数量的理论值与模拟值的对比结果见图3。由图3可见，4种主机的模拟值与理论值总体趋势一致，其中感染主机和恢复主机的数量曲线几乎完全重合，表明这2种主机在传播过程中的转换速率及比例符合理论预期。易感主机的数量在最低点处存在微小偏差，但随时间的增加逐渐减小为0，这表明模型参数在动态调控下能够准确捕捉易感主机状态变化的细节。此外，隔离主机的数量在峰值处略有偏差，但随着时间的推移逐渐减小为0。仿真实验结果验证了理论计算的准确性，表明该模型能够准确描述不同类型主机在感染过程中的状态转换规律，有效预测恶意软件在工业网络中的传播趋势。

3.2 R_{\mathrm{0}}>\mathrm{1}时的实验分析

当R_{\mathrm{0}}>\mathrm{1}时，初始感染主机是5台，其他主机状态均为易感状态。为验证第二部分关于时延结论的正确性，进行如下数值实验，其他参数设置与时相同。

（1） τ对系统稳定性影响

τ对5种主机数量的影响见图4，由图4（a）可见，当，a=10^-5时，在t=250 s处，各类主机都趋近于理论值，此时防御措施能迅速响应感染动态，系统中各状态变量的调整过程呈现平滑收敛的趋势，系统在平衡点E^*处为局部渐近稳定。此时，系统表现为线性稳定状态，易感主机、感染主机、隔离主机和恢复主机之间的转换遵循预期的动态规律，各状态之间的相互作用达到平衡。

由图4（b）可见，当\tau =\mathrm{70}>{\mathrm{\tau }}_{\mathrm{0}}、a=10^-5时，曲线波动较大，此时防御措施响应滞后，系统中恶意软件的传播过程呈周期性振荡。这种周期性振荡现象正是Hopf分岔所预测的动态不稳定状态的特征，说明系统由于延迟反馈效应在局部失去了稳定。此外，这种动态不稳定现象在实际工业互联网中意味着防护系统可能无法及时、有效地抑制恶意软件的传播，易引发大规模安全事件。

基于第2部分的理论推导和Hopf分岔理论分析，得到a=10^-5时的系统分岔图，见图5。通过计算雅可比矩阵的特征值，在\tau为58附近发现一对共轭纯虚根，此时系统的实部正好穿过0值点，出现Hopf分岔。当系统时延参数τ大于58时，系统将从稳定状态转变为周期性振荡状态，说明工业互联网中防御措施的时延效应对系统稳定性有显著影响。在实际应用中，如果防御系统的响应时延超过该临界值，恶意软件可能会以周期性爆发的方式扩散，从而增加预测和防御的难度。基于此，在实际网络中可以进一步预测工业互联网中的稳定性边界，为制定有效的网络安全抑制策略提供理论依据。

为进一步验证时延对工业互联网的影响，探讨不同时延下，同一坐标内感染主机台数的变化情况（见图6），τ分别取20、30、50、70。

由图6可知，当时延较小时，系统相对稳定。当时延达到阈值后，系统稳定性逐渐降低，开始出现振荡现象。这表明时延的增加可延长系统的响应时间，导致隔离策略难以快速实施，使恶意软件的持续传播。在实际应用中，当工业互联网中的时延较长时，防御系统必须具备较高的应急响应能力，例如快速调整隔离策略、增加系统的防护力度，以应对恶意软件的传播。

（2） 参数a对系统稳定性影响

不同参数a对被感染主机的数量的影响见图7，由图7（a）可知，a从10^-7增至10^-4时，感染主机的峰值数量由8×10⁴ 台上升至2.15×10⁵ 台。由图7（b）可知，a越大，感染主机数量变化曲线的振幅越大。这表明，较大的a值会使隔离措施过于激进，影响网络的正常运行和系统效率。但较小的a值会使系统对感染主机的响应不足。因此，在工业互联网的实际应用中，防御系统应根据实时的网络状况和感染态势，动态调整a值，以实现既能有效控制恶意软件传播，又能保证系统高效运行的需求。

（3） 感染率对系统稳定性影响

当τ<τ₀ 和τ>τ₀时，不同感染率对被感染主机台数的影响情况见图8。由图8可见，感染率越大，感染主机台数的峰值也随之增大，表明感染率直接影响了恶意软件的传播范围和扩散速度。当感染率较高时，恶意软件能够迅速传播至更多主机，从而使感染主机的数量迅速增加，且台数的峰值明显更高。

对比图8（a）和图8（b）可见，时延对恶意软件传播影响更显著。在时延较小时，系统能够有效抑制高感染率下的恶意软件传播，保持系统的稳定性。而在时延较大时，即使降低感染率，恶意软件的扩散仍不可控，导致系统不稳定。因此，综合考虑感染率和时延以及系统的响应能力，是确保工业互联网安全的核心逻辑。在实际应用中，应通过优化网络架构和提高防御系统响应能力，减少时延的影响，确保恶意软件的传播能够在初期得到有效控制。

（4） α对系统稳定性影响

当τ＜τ₀ 和τ＞τ₀时，不同α对被感染主机的数量的影响情况见图9。由图9（a）可见，α较大时，恶意软件传播受到有效抑制，工业互联网能够快速恢复安全状态；而当α较小时，感染范围扩大，感染主机数量增速加快且峰值显著提高。由图9（b）可见，在α较小时，恶意软件会持续存在并长期传播；随着隔离率提高，感染主机数量迅速减少，证明提高隔离率能有效遏制恶意软件扩散，但由于时延较大，系统仍不稳定。综合分析图9可知，隔离率与感染主机台数的峰值呈负相关关系，α越小，恶意软件传播越难控制，感染规模也越大。因此，在实际应用中，为防止恶意软件的蔓延，必须根据传播阶段的不同动态调整隔离率。在感染初期，系统需要保证一定的运行效率，可以维持较低的隔离率，以减少过早的资源消耗和对系统性能的影响；而当感染扩散加速时，必须及时提高隔离率，以控制感染范围，确保系统不陷入失控状态。

（5） \gamma对系统稳定性影响

不同\gamma对感染主机数量的影响见图10。由图10可见，当\gamma较高时，感染主机能快速转变为恢复状态，有效减缓恶意软件传播速度，促使系统趋于稳定；反之，当恢复率较低时，主机恢复速率显著降低，恶意软件传播周期延长，扩散风险加剧。低恢复率会延长感染主机的恢复时间，增加恶意软件的持续传播时间。这种现象不仅会增加恶意软件威胁的消除难度，还会显著加剧工业互联网的安全风险。

由图10（a）可知，时延小于阈值时，恢复率的优化能够更为有效地控制恶意软件的传播，减少潜在的安全威胁。而由图10（b）可知，时延过大时，防御系统对恶意软件的响应变得迟缓，恢复措施无法及时有效地隔离和清除感染主机，导致恶意软件的传播失控。在这种情况下，提高恢复率的效果被时延的影响所削弱，系统的稳定性无法得到有效保障。

因此，在实际应用中，除了提高恢复率外，还应优化网络的时延特性，减少时延对系统恢复的负面影响，从而构建更为高效和稳定的防御体系。

综上，在工业互联网中，通过合理选择参数值可有效控制恶意软件传播，进而制定更优化的防御策略。在实际应用中，可以根据网络的实际情况和攻击的不同阶段，合理调整这些参数，灵活应对恶意软件的传播趋势，避免系统陷入失控状态。

（6） 理论值与模拟值对比

为了更加准确地描述恶意软件的传播行为，本部分设计了仿真实验来验证理论分析和数值模拟结果的一致性和正确性：

4种主机在τ=20＜τ₀时的理论值与模拟值的对比情况见图11，由图11可知，4种主机的模拟值与理论值的总体趋势变化一致。证明数值实验可以准确反映恶意软件传播的动态过程，模型在不同阶段的表现与实际系统的运行情况相符。易感主机的数量在初期迅速下降，而感染主机数量在短时间内急剧攀升至峰值后持续衰减。这表明，防护机制未完全介入之前的传播初期，恶意软件快速传播并迅速感染大量主机。随着隔离机制的介入，传播范围得到抑制，最终促使工业互联网系统恢复稳定状态。这一规律验证了模型对恶意软件传播动力学特征的准确描述。

（7） 对比实验

为验证模型的有效性，通过数值仿真实验，在相同参数设置下对比SIDQV模型与经典SIDR模型（不考虑隔离率的变化），结果见图12，由图12可知，本文提出的SIDQV模型的感染峰值明显低于SIDR模型。说明在相同的隔离资源下，动态调整隔离率可以有效降低感染峰值和感染总量，控制恶意软件的传播。SIDQV模型通过实时调整隔离率应对网络中的感染主机数量，使得系统能够灵活响应感染的变化。在传播初期，系统可采用较低的隔离率，以保证网络的流畅性和资源的有效利用。随着感染主机数量的增加，系统则会自动提高隔离率，限制更多主机的感染。相比之下，SIDR模型始终依赖固定的隔离策略，导致感染主机的数量在较长时间内维持较高水平，无法如SIDQV模型那样快速有效地遏制恶意软件的扩散。因此，引入可变隔离率和时延因素，可使SIDQV模型能够提供比传统静态模型更精确的传播预测与更高效的防控措施。

4 结论

（1）可变隔离率的SIDQV模型在特定条件下具有一个正平衡点，且当时延超过阈值时，系统会失去稳定并出现Hopf分岔现象。

（2）在参数合理调控下，采用动态隔离策略可降低恶意软件的感染峰值，并减小其恶意软件的传播速度和范围，增强工业互联网的防护能力。因此，本文提出的模型和抑制策略能够为工业互联网的防护提供明确的参数指导，提高其对恶意软件传播的应对能力与稳定性。

（3）模型仍存在一定的局限性。首先，在处理超大规模网络时，模型可能会面临计算资源和时间方面的挑战。其次，模型参数调整的过程较为复杂，需要针对不同的网络环境进行细致调整。在未来的工作中，可以考虑更复杂的网络拓扑结构并引入更多的实际约束，以进一步提高模型的逼真度和工程实用性。

参考文献

[1] WANG C, JIANG C X, WANG J J,et al.Blockchain-aided network resource orchestration in intelligent Internet of Things[J].IEEE Internet of Things Journal,2023,10(7):6151-6163.

[2] KUN X, WANG Z L, ZHOU Z A,et al.Design of industrial Internet of Things system based on machine learning and artificial intelligence technology[J]. Journal of Intelligent & Fuzzy Systems, 2021, 40(2): 2601-2611.

[3] MENDONÇA R V, SILVA J C, ROSA R L,et al.A lightweight intelligent intrusion detection system for industrial Internet of Things using deep learning algorithms[J].Expert Systems,2022,39(5):e12917.

[4] KHAN I A, KESHK M, PI D C,et al.Enhancing IIoT networks protection: a robust security model for attack detection in Internet Industrial Control Systems[J].Ad Hoc Networks,2022,134:102930.

[5] TRAUTMAN L J.Industrial cyber vulnerabilities: lessons from stuxnet and the Internet of Things[J].SSRN Electronic Journal,2017,72:761-826.

[6] BINDRA A. Securing the power grid: protecting smart grids and connected power systems from cyberattacks[J]. IEEE Power Electronics Magazine, 2017, 4(3): 20-27.

[7] MAKRAKIS G M, KOLIAS C, KAMBOURAKIS G, et al. Industrial and critical infrastructure security: technical analysis of real-life security incidents[J]. IEEE Access, 2021, 9: 165295-165325.

[8] AKBARZADEH A.Dependency based risk analysis in Cyber-Physical Systems[D].Trondheim:Norwegian University of Science and Technology,2023:5-8.

[9] FLÓREZ J, MONTOYA G A, LOZANO-GARZÓN C.An epidemiological SIS malware spreading model based on Markov chains for IoT networks[C]//Proceedings of Eighth International Congress on Information and Communication Technology.Singapore: Springer Nature Singapore, 2023:663-672.

[10] JHUN B,JO M, KAHNG B.Simplicial SIS model in scale-free uniform hypergraph[J].Journal of Statistical Mechanics:Theory and Experiment,2019,12(12):123207.

[11] DADLANI A, KUMAR M S, KIM K,et al.Stability and immunization analysis of a malware spread model over scale-free networks[J].IEEE Communications Letters,2014,18(11):1907-1910.

[12] DEL REY A M, VARA R C, GONZÁLEZ S R.A computational propagation model for malware based on the SIR classic model[J]. Neurocomputing,2022,484:161-171.

[13] TAGHVAEI A, GEORGIOU T T, NORTON L,et al.Fractional SIR epidemiological models[J].Scientific Reports,2020,10:20882.

[14] TAYNITSKIY V, GUBAR E, ZHU Q Y.Optimal impulse control of bi-virus SIR epidemics with application to heterogeneous Internet of Things[C]//2017 Constructive Nonsmooth Analysis and Related Topics (Dedicated to the Memory of V.F.Demyanov).May 22-27,2017,St. Petersburg,Russia.IEEE,2017:1-4.

[15] XIAO M, CHEN S, ZHENG W X,et al.Tipping point prediction and mechanism analysis of malware spreading in cyber-physical systems[J]. Communications in Nonlinear Science and Numerical Simulation, 2023,122:107247.

[16] DEL REY A M, ENCINAS A H, VAQUERO J M,et al.A method for malware propagation in industrial critical infrastructures[J].Integrated Computer-Aided Engineering,2016,23(3):255-268.

[17] MAHBOUBI A, CAMTEPE S, MORARJI H.A study on formal methods to generalize heterogeneous mobile malware propagation and their impacts[J].IEEE Access,2017,5:27740-27756.

[18] YAO Y, FU Q, YANG W,et al.An epidemic model of computer worms with time delay and variable infection rate[J].Security and Communication Networks,2018,2018:9756982.

[19] YU Z H, GAO H X, WANG D,et al.SEI²RS malware propagation model considering two infection rates in cyber-physical systems[J].Physica A:Statistical Mechanics and its Applications,2022,597:127207.

[20] YANG W, FU Q, YAO Y.A malware propagation model with dual delay in the industrial control network[J].Complexity,2023,2023: 8823080.

基金资助

辽宁省自然基金科学项目(2022-MS-291┫辽宁省教育厅高校基本科研项目┣LJKMZ20220781)