具有接种效应带时变时滞的SVEIR模型的稳定性分析

Stability Analysis of SVEIR Models with Inoculated Effect and Time-Varying De-lay

摘要: 考虑了一类受接种疫苗影响的带时变时滞的SVEIR疾病传播模型，针对不同取值下的病毒基本再生数R₀，分别采取Liapunov函数法和Routh-Hurwitz判别法对模型平衡点的存在性、稳定性、分岔性进行研究。最后对不同的R₀进行数值模拟，验证理论分析结果的正确性。

Abstract: In this paper, we consider a SVEIR disease transmission model with time-varying delay affected by vaccination, for the basic regeneration number of virus R₀ under different values, Liapunov func-tion method and Routh-Hurwitz criterion were used to study the existence, stability and bifurcation of model equilibrium points. Finally, the numerical simulation of different R₀ is carried out to verify the correctness of theoretical analysis results.

1. 引言

自2020年新型冠状病毒发生以来，已有众多学者致力于新冠肺炎疫情传播机理与特性的研究，并建立各种新冠肺炎病毒传播的动力学模型，通过模型对疫情的发展趋势做预测分析，详见参考文献 [1] [2] [3] [4]。还有许多学者针对经典的SEIR模型通过改变模型参数研究各类病毒传播模型的动力学性质，如参考文献 [5] [6] [7]。但是由于新冠肺炎的潜伏期较长，具有一定的时变滞后性，因此，控制时滞也是一种常见的模型引入方式。很多学者开始选择采用时变系统描述传染病的动力学特性，如文献 [8] 基于季节性变化的SIR模型研究了流感类传染病的疫苗接种和治疗策略。文献 [9] 将传染率定义为时间、易感染人口数、感染人口数和总人口数的非线性函数。文献 [10] 讨论了潜伏期时滞的时变SEIR模型的最优疫苗接种策略。文献 [11] 研究了基于时滞SIR模型关于H1N1型流感的最优疫苗接种策略，并考虑了疫苗起效的延迟。

本文在经典SEIR模型的基础上，通过引入潜伏期时滞因素、疫苗接种率因素及接种疫苗的免疫率因素，更加真实地描述传染病的传播规律。研究了一类具有疫苗接种影响效应的带时变时滞的SVEIR模型，其微分方程形式如式(1)所示。数值仿真结果表明，在满足各种约束的条件下，接种疫苗可以有效地抑制传染病的传播。

$\{\begin{array}{l}{S}^{\prime }\left(t\right)=a-\left(d+b\left(t\right)\right)S\left(t\right)-\left({\alpha }_{1}E\left(t\right)+{\alpha }_{2}I\left(t\right)\right)S\left(t\right)\\ V^{\prime }\left(t\right)=b\left(t\right)S\left(t\right)-dV\left(t\right)-\left({\beta }_{1}E\left(t\right)+{\beta }_{2}I\left(t\right)\right)V\left(t\right)\\ E^{\prime }\left(t\right)=\left({\alpha }_{2}S\left(t-\tau \right)+{\beta }_{2}V\left(t-\tau \right)\right)I\left(t-\tau \right)-\left(d+k+\alpha \left({\gamma }^{\prime }\right)\right)E\left(t\right)\\ I^{\prime }\left(t\right)=kE\left(t\right)-\left(d+e+r\right)I\left(t\right)\\ R^{\prime }\left(t\right)=\alpha \left({\gamma }^{\prime }\right)E\left(t\right)+rI\left(t\right)-dR\left(t\right)+b\left(t\right)c\left(t\right)S\left(t\right)\end{array}$ (1)

其中， $S\left(t\right)$ ：为易感者； $V\left(t\right)$ ：接种疫苗者； $E\left(t\right)$ ：为潜伏者； $I\left(t\right)$ ：为感染者；

$R\left(t\right)$ ：为治愈者；a：人口常数输入率；d：自然死亡率；k：潜伏者的确诊率；

${\alpha }_1$ ：易感者与潜伏期接触时的传播率； ${\alpha }_2$ ：易感者与确诊患者接触时的传播率；

${\beta }_1$ ：接种者与潜伏期接触时的传播率； ${\beta }_2$ ：接种者与确诊患者接触时的传播率；

${\gamma }_1$ ：接种者人群潜伏期移出率； ${\gamma }_2$ ：未接种者人群潜伏期移出率； $\delta$ ：疾病传播率；

$\alpha \left({\gamma }^{\prime }\right)={\gamma }_1\delta +\left(1-\delta \right){\gamma }_2$ ：潜伏期人群受疫苗接种信息影响的移出率； $\tau$ ：平均潜伏期；

e：因病死亡率；r：染病者的治愈率； $b\left(t\right)$ ：易感者的疫苗接种率，满足 $0\le b\left(t\right)\le 1$ ；

$c\left(t\right)$ ：接种人群的成功免疫率，且 $c\left(t\right)=c_0+c_1\exp \left(-t/c_2\right)$，其中， $c_0,c_1,c_2$ 为非负参数，并满足： $c\left(0\right)=c_0+c_1$， $c\left(+\infty \right)=c_0$。

2. 模型的建立

由于系统(1)里 $R\left(t\right)$ 项是独立存在的，故只需研究系统(1)的子系统(2)。

$\{\begin{array}{l}{S}^{\prime }\left(t\right)=a-\left(d+b\left(t\right)\right)S\left(t\right)-\left({\alpha }_{1}E\left(t\right)+{\alpha }_{2}I\left(t\right)\right)S\left(t\right)\\ V^{\prime }\left(t\right)=b\left(t\right)S\left(t\right)-dV\left(t\right)-\left({\beta }_{1}E\left(t\right)+{\beta }_{2}I\left(t\right)\right)V\left(t\right)\\ E^{\prime }\left(t\right)=\left({\alpha }_{2}S\left(t-\tau \right)+{\beta }_{2}V\left(t-\tau \right)\right)I\left(t-\tau \right)-\left(d+k+\alpha \left({\gamma }^{\prime }\right)\right)E\left(t\right)\\ I^{\prime }\left(t\right)=kE\left(t\right)-\left(d+e+r\right)I\left(t\right)\end{array}$ (2)

由现实生活可知， $S\left(t\right),V\left(t\right),E\left(t\right),I\left(t\right)$ 都是正数，故系统(2)的初始状态空间为 $\Omega =\{\left(S\left(t\right),V\left(t\right),E\left(t\right),I\left(t\right)\right)|S\left(t\right)\ge 0,V\left(t\right)\ge 0,E\left(t\right)\ge 0,I\left(t\right)\ge 0\}$。

由于 $\Omega$ 具有正不变性，因此，只需考虑初始条件处于 $\Omega$ 内的解，并验证解在 $\Omega$ 内具有存在唯一性。

3. 平衡点的存在性

令系统(2)的微分方程右侧等于0，计算得无病平衡点 $E_0$ 和地方病平衡点 $E_{e1}^{*},E_{e2}^{*}$。

$E_0=\left(\frac{a}{d+b\left(t\right)},\frac{ab\left(t\right)}{d\left(d+b\left(t\right)\right)},0,0\right)$， $E_{e1}^{*}=\left(S_1^{*},V_1^{*},E_1^{*},I_1^{*}\right)$， $E_{e2}^{*}=\left(S_2^{*},V_2^{*},E_2^{*},I_2^{*}\right)$

令 $M=d+k+\alpha \left({\gamma }^{\prime }\right)$， $N=d+e+r$

则 $I_1^{*},I_2^{*}$ 是方程(3)的解。

$\begin{array}{l}MN\left({\alpha }_{1}N+{\alpha }_{2}k\right)\left({\beta }_{1}N+{\beta }_{2}k\right)I^2+\{MN\left[kd\left({\alpha }_{1}N+{\alpha }_{2}k\right)+k\left(d+b\left(t\right)\right)\left({\beta }_{1}N+{\beta }_{2}k\right)\right]\\ -ak\left({\alpha }_{1}N+{\alpha }_{2}k\right)\left({\beta }_{1}N+{\beta }_{2}k\right)\}I+MN{k}^{2}d\left(d+b\left(t\right)\right)-a{k}^2\left(d\left({\alpha }_{1}N+{\alpha }_{2}k\right)+b\left(t\right)\left({\beta }_{1}N+{\beta }_{2}k\right)\right)=0\end{array}$ (3)

且， $E^{*}=\frac{N}{k}{I}^{*}$， $S^{*}=\frac{a}{d+b\left(t\right)+{\alpha }_1{E}^{*}+{\alpha }_2{I}^{*}}$， $V^{*}=\frac{b\left(t\right)}{d+{\beta }_1{E}^{*}+{\beta }_2{I}^{*}}{S}^{*}$

若令 $u_1=\frac{{\alpha }_1}{k}N+{\alpha }_2$， $u_2=\frac{{\beta }_1}{k}N+{\beta }_2$

则(3)式可转变为(4)式

$MN{k}^2{u}_1{u}_2{I}^2+\left\{MN\left[k^{2}d{u}_1+k^2\left(d+b\left(t\right)\right)u_2\right]-a{k}^3{u}_1{u}_2\right\}I+MN{k}^{2}d\left(d+b\left(t\right)\right)-a{k}^3\left(d{u}_1+b\left(t\right)u_2\right)=0$ (4)

(4)式两边同时消去 $k^2$，进而可得(5)式

$MN{u}_1{u}_2{I}^2+\left\{MN\left[d{u}_1+\left(d+b\left(t\right)\right)u_2\right]-ak{u}_1{u}_2\right\}I+MNd\left(d+b\left(t\right)\right)-ak\left(d{u}_1+b\left(t\right)u_2\right)=0$ (5)

4. 再生系数

采用文献 [12] 中的下一代矩阵法(Next Generation Matrix)计算病毒基本再生数，仅考虑携带病毒潜伏者和感染者情形，作如下记号

$\begin{array}{l}\{\begin{array}{l}{F}_E=\left({\alpha }_{2}S\left(t-\tau \right)+{\beta }_{2}V\left(t-\tau \right)\right)I\left(t-\tau \right)\\ F_I=0\end{array}\\ \{\begin{array}{l}{V}_E=\left(d+k+\alpha \left({\gamma }^{\prime }\right)\right)E\left(t\right)\\ V_I=-kE\left(t\right)+\left(d+e+r\right)I\left(t\right)\end{array}\end{array}$ (6)

令

$F=\left(\begin{array}{cc}\frac{\partial F_E}{\partial E}& \frac{\partial F_E}{\partial I}\\ \frac{\partial F_I}{\partial E}& \frac{\partial F_I}{\partial I}\end{array}\right)$， $V=\left(\begin{array}{cc}\frac{\partial V_E}{\partial E}& \frac{\partial V_E}{\partial I}\\ \frac{\partial V_I}{\partial E}& \frac{\partial V_I}{\partial I}\end{array}\right)$， (7)

由(6)和(7)联立得， $F=\left(\begin{array}{cc}0& {\alpha }_{2}S\left(t-\tau \right)+{\beta }_{2}V\left(t-\tau \right)\\ 0& 0\end{array}\right)$， $V=\left(\begin{array}{cc}d+k+\alpha \left({\gamma }^{\prime }\right)& 0\\ -k& d+e+r\end{array}\right)=\left(\begin{array}{cc}M& 0\\ -k& N\end{array}\right)$

进而有，

$\begin{array}{c}F\cdot V^{-1}=\left(\begin{array}{cc}0& {\alpha }_{2}S\left(t-\tau \right)+{\beta }_{2}V\left(t-\tau \right)\\ 0& 0\end{array}\right)\cdot \frac{1}{MN}\left(\begin{array}{cc}N& 0\\ k& M\end{array}\right)\\ =\left(\begin{array}{cc}0& {\alpha }_{2}S\left(t-\tau \right)+{\beta }_{2}V\left(t-\tau \right)\\ 0& 0\end{array}\right)\cdot \left(\begin{array}{cc}\frac{1}{M}& 0\\ \frac{k}{MN}& \frac{1}{N}\end{array}\right)\\ =\left(\begin{array}{cc}\frac{k\left({\alpha }_{2}S\left(t-\tau \right)+{\beta }_{2}V\left(t-\tau \right)\right)}{MN}& \frac{{\alpha }_{2}S\left(t-\tau \right)+{\beta }_{2}V\left(t-\tau \right)}{N}\\ 0& 0\end{array}\right)\\ =\left(\begin{array}{cc}\frac{k\left({\alpha }_2\frac{a}{d+b\left(t\right)}+{\beta }_2\frac{a}{d+b\left(t\right)}\frac{b\left(t\right)}{d}\right)}{MN}& \frac{{\alpha }_2\frac{a}{d+b\left(t\right)}+{\beta }_2\frac{a}{d+b\left(t\right)}\frac{b\left(t\right)}{d}}{N}\\ 0& 0\end{array}\right)\end{array}$

所以系统(2)的再生系数谱半径为

$R_0=\rho \left(F\cdot V^{-1}\right)=\frac{k\left({\alpha }_2\frac{a}{d+b\left(t\right)}+{\beta }_2\frac{a}{d+b\left(t\right)}\cdot \frac{b\left(t\right)}{d}\right)}{MN}$。

5. 平衡点的稳定性分析

定理1：当 $R_0=1$ 时，系统(2)在无病平衡点 $E_0$ 处存在分岔。

证明：系统(2)在无病平衡点 $E_0$ 处的雅可比矩阵为

$J_{E_0}=\left(\begin{array}{cccc}-\left(d+b\left(t\right)\right)& 0& -\frac{a{\alpha }_1}{d+b\left(t\right)}& -\frac{a{\alpha }_2}{d+b\left(t\right)}\\ b\left(t\right)& -d& -\frac{a{\beta }_{1}b\left(t\right)}{d\left(d+b\left(t\right)\right)}& -\frac{a{\beta }_{2}b\left(t\right)}{d\left(d+b\left(t\right)\right)}\\ 0& 0& -M& \frac{a{\alpha }_2}{d+b\left(t\right)}+\frac{ab\left(t\right){\beta }_2}{d\left(d+b\left(t\right)\right)}\\ 0& 0& k& -N\end{array}\right)$

则可求得 $J_{E_0}$ 的特征根为

$\begin{array}{l}{\lambda }_1=-d\\ {\lambda }_2=-\left(d+b\left(t\right)\right)\\ {\lambda }_3=\frac{-\left(M+N\right)-\sqrt{{\left(M-N\right)}^2+4\frac{\left(ad{\alpha }_2+ab\left(t\right){\beta }_2\right)k}{d\left(d+b\left(t\right)\right)}}}{2}\\ {\lambda }_4=\frac{-\left(M+N\right)+\sqrt{{\left(M-N\right)}^2+4\frac{\left(ad{\alpha }_2+ab\left(t\right){\beta }_2\right)k}{d\left(d+b\left(t\right)\right)}}}{2}\end{array}$

当 $R_0=1$ 时，易知， $\frac{k\left(ad{\alpha }_2+ab\left(t\right){\beta }_2\right)}{d\left(d+b\left(t\right)\right)}=MN$

所以， ${\lambda }_4=0$

故， $E_0$ 为非双曲平衡点，则系统(2)在无病平衡点 $E_0$ 处存在分岔现象。

定理2：当 $R_0<1$ 时，系统(2)在无病平衡点 $E_0$ 处渐进稳定。

证明：当 $R_0<1$ 时， $\frac{k\left(ad{\alpha }_2+ab\left(t\right){\beta }_2\right)}{d\left(d+b\left(t\right)\right)}<MN$，由定理1可知，

${\lambda }_4=\frac{-\left(M+N\right)+\sqrt{{\left(M-N\right)}^2+4\frac{\left(ad{\alpha }_2+ab\left(t\right){\beta }_2\right)k}{d\left(d+b\left(t\right)\right)}}}{2}<\frac{-\left(M+N\right)+\sqrt{{\left(M-N\right)}^2+4MN}}{2}=0$

所以可得 ${\lambda }_1<0,{\lambda }_2<0,{\lambda }_3<0,{\lambda }_4<0$。

由Lyapnunov定理可知，系统(2)在无病平衡点 $E_0$ 处渐进稳定。

定理3：当 $R_0>1$ 时，系统(2)在无病平衡点 $E_0$ 不稳定。

证明：由定理1可知， ${\lambda }_1<0,{\lambda }_2<0,{\lambda }_3<0$

${\lambda }_4=\frac{-\left(M+N\right)+\sqrt{{\left(M-N\right)}^2+4\frac{\left(ad{\alpha }_2+ab\left(t\right){\beta }_2\right)k}{d\left(d+b\left(t\right)\right)}}}{2}>\frac{-\left(M+N\right)+\sqrt{{\left(M-N\right)}^2+4MN}}{2}=0$

故，系统(2)在无病平衡点 $E_0$ 不稳定。

定理4：当 $R_0>1$ 时，系统(2)在地方病平衡点 $E_e^{*}$ 局部渐进稳定。

证明：由于系统(2)在地方病平衡点 $E_{e1}^{*},E_{e2}^{*}$ 处的证明方法类似，故取 $E_{e1}^{*}=E_e^{*}$，可得到系统(2)在地方病平衡点 $E_e^{*}$ 处的雅可比矩阵为

$J_{E_e}^{*}=\left(\begin{array}{cccc}-d-b\left(t\right)-{\alpha }_1{E}^{*}-{\alpha }_2{I}^{*}& 0& -{\alpha }_1{S}^{*}& -{\alpha }_2{S}^{*}\\ b\left(t\right)& -d-{\beta }_1{E}^{*}-{\beta }_2{I}^{*}& -{\beta }_1{V}^{*}& -{\beta }_2{V}^{*}\\ {\alpha }_2{I}^{*}& {\beta }_2{I}^{*}& -M& {\alpha }_2{S}^{*}+{\beta }_2{V}^{*}\\ 0& 0& k& -N\end{array}\right)$

则可求得 $J_{E_e}^{*}$ 的特征方程为 $a_4{\lambda }^4+a_3{\lambda }^3+a_2{\lambda }^2+a_1\lambda +a_0=0$

其中，

$\begin{array}{l}{a}_4=1\\ a_3=M+N+2d+b\left(t\right)+\left({\alpha }_1+{\beta }_1\right)E^{*}+\left({\alpha }_2+{\beta }_2\right)I^{*}\\ a_2=MN+\left(M+N\right)\left(d+b\left(t\right)+{\alpha }_1{E}^{*}+{\alpha }_2{I}^{*}+d+{\beta }_1{E}^{*}+{\beta }_2{I}^{*}\right)+\left(d+{\beta }_1{E}^{*}+{\beta }_2{I}^{*}\right)\left(d+b\left(t\right)+{\alpha }_1{E}^{*}+{\alpha }_2{I}^{*}\right)\\ +{\alpha }_1{\alpha }_2{S}^{*}{I}^{*}+{\beta }_1{\beta }_2{V}^{*}{I}^{*}-k\left({\alpha }_2{S}^{*}+{\beta }_2{V}^{*}\right)\end{array}$

$\begin{array}{l}{a}_1=MN\left(d+b\left(t\right)+{\alpha }_1{E}^{*}+{\alpha }_2{I}^{*}+d+{\beta }_1{E}^{*}+{\beta }_2{I}^{*}\right)+\left(M+N\right)\left(d+{\beta }_1{E}^{*}+{\beta }_2{I}^{*}\right)\left(d+b\left(t\right)+{\alpha }_1{E}^{*}+{\alpha }_2{I}^{*}\right)\\ +\left(b\left(t\right){\alpha }_1{\beta }_2+k{\alpha }_2^2+{\alpha }_1{\alpha }_2\left(N+d+{\beta }_1{E}^{*}+{\beta }_2{I}^{*}\right)\right)S^{*}{I}^{*}+\left({\beta }_1{\beta }_2\left(N+d+b\left(t\right)+{\alpha }_1{E}^{*}+{\alpha }_2{I}^{*}\right)+k{\beta }_2^2\right)V^{*}{I}^{*}\\ -k\left({\alpha }_2{S}^{*}+{\beta }_2{V}^{*}\right)\left(d+b\left(t\right)+{\alpha }_1{E}^{*}+{\alpha }_2{I}^{*}+d+{\beta }_1{E}^{*}+{\beta }_2{I}^{*}\right)\\ a_0=\left(N{\alpha }_1+k{\alpha }_2\right)\left(b\left(t\right){\beta }_2+{\alpha }_2\left(d+{\beta }_1{E}^{*}+{\beta }_2{I}^{*}\right)\right)S^{*}{I}^{*}+\left(N{\beta }_1+k{\beta }_2\right){\beta }_1\left(d+b\left(t\right)+{\alpha }_1{E}^{*}+{\alpha }_2{I}^{*}\right)V^{*}{I}^{*}\\ +\left(d+{\beta }_1{E}^{*}+{\beta }_2{I}^{*}\right)\left(d+b\left(t\right)+{\alpha }_1{E}^{*}+{\alpha }_2{I}^{*}\right)\left(MN-k\left({\alpha }_2{S}^{*}+{\beta }_2{V}^{*}\right)\right)\end{array}$

由Routh-Hurwitz判定准则得

$\begin{array}{l}{\Delta }_1=a_4=1\\ {\Delta }_2=\left(\begin{array}{cc}{a}_3& a_1\\ a_4& a_2\end{array}\right)\\ {\Delta }_3=\left(\begin{array}{ccc}{a}_3& a_1& 0\\ a_4& a_2& a_0\\ 0& a_3& a_1\end{array}\right)\\ {\Delta }_4=\left(\begin{array}{cccc}{a}_3& a_1& 0& 0\\ a_4& a_2& a_0& 0\\ 0& a_3& a_1& 0\\ 0& a_4& a_2& a_0\end{array}\right)\end{array}$

易知， $a_0,a_1,a_2,a_3,a_4$ 均大于0，所以可得

$\begin{array}{l}{\Delta }_1>0\\ {\Delta }_2=a_2{a}_3-a_1{a}_4>0\\ {\Delta }_3=a_1{a}_2{a}_3-a_1^2{a}_4-a_0{a}_3^2>0\\ {\Delta }_4=a_0{\Delta }_3>0\end{array}$

由Routh-Hurwitz准则可知，系统(2)在地方病平衡点 $E_e^{*}$ 处局部渐进稳定。

6. 数值模拟实验与结论

通过数值模拟，验证理论分析结果的正确性。

1) 取 $a=0.01$， $d=0.03$， $b\left(t\right)=0.8$， $c\left(t\right)=0.5$， $k=0.9$， ${\alpha }_1=0.6$， ${\alpha }_2=0.85$， ${\beta }_1=0.1$， ${\beta }_2=0.1$， $e=0.7$， $r=0.9$，则 $R_0<1$ 满足所求。随着接种疫苗人数的增加，潜伏期人数越来越少，最终易感患者人数将趋向0，故无病平衡点 $E_0$ 全局渐近稳定，疫情最终会得到控制，如图1所示。

Figure 1. When $R_0<1$, the trends of disease-free equilibrium point

图1. $R_0<1$ 时无病平衡点的变化趋势

2) 取 $a=0.1$， $d=0.05$， $b\left(t\right)=0.3$， ${\alpha }_2=0.45$， ${\beta }_2=0.2$， $e=0.3$， $r=0.5$，则 $R_0>1$。

随着疫苗接种人数达到饱和时，接种人数会越来越少，接种疫苗的免疫率会随着接种疫苗的时间发生变化，此时潜伏期人数越来越多，易感患者人数也会在一定程度上增加，会出现大面积爆发的情况，地方病平衡点 $E_e^{*}$ 局部渐近稳定，疾病会演化为当地流行病。如图2所示。

Figure 2. When $R_0>1$, the trends in the equilibrium of endemic diseases

图2. $R_0>1$ 时地方病平衡点的变化趋势

3) 取 $a=0.02$， $d=0.003$， $b\left(t\right)=0.7$， $k=0.25$， ${\alpha }_1=0.15$， ${\alpha }_2=0.35$， ${\beta }_1=0.15$， ${\beta }_2=0.25$， ${\gamma }_1=0.7$， ${\gamma }_2=0.3$， $\delta =0.3$， $e=0.1$， $r=0.9$，则 $R_0=1$。经过一段时间后明显出现分岔现象，在分界处，潜伏期和确诊人数比例会发生显著变化。如图3所示。

Figure 3. When $R_0=1$, the trends in the equilibrium of endemic diseases

图3. $R_0=1$ 时地方病平衡点的变化趋势

本文重点讨论了受接种疫苗的比例及疫苗免疫成功率影响的SVEIR模型，通过以上讨论可知，当再生系数 $R_0<1$ 时，系统会存在一个无病平衡点，随着政府的管控，疾病会消失；当再生系数 $R_0>1$ 时，系统平衡点局部稳定，且疾病会转化为地方流行病；当再生系数 $R_0=1$ 时，系统会出现临界分岔现象，同时也印证了接种疫苗是疾病防控的关键措施之一。

基金项目

福建省教育教学教改项目(FBJG20170154)；漳州市自然科学基金项目(ZZ2018J26)；福建省高校产学合作项目(2018H6018)；教育部产学协同育人项目(JGH2019003和JGH2019023)。