CN112700885B - 一种基于卡尔曼滤波辨识新冠病毒传播模型参数的方法 - Google Patents

Abstract

本发明提供一种基于卡尔曼滤波理论辨识新型冠状病毒传播模型参数的方法，用于解决现有技术中存在的传染病模型不考虑潜伏期以及传染率不考虑时间变化影响的技术问题，实现步骤为：获取新型冠状肺炎疫情数据集，建立传染病传播模型，使用卡尔曼滤波算法补充所述传染病传播模型中的未知参数，再对识别出的传染率建模成随时间递减的时变函数，建立传染率变化模型，根据传染率变化模型更新所述传染病传播模型得到最优化传染病传播模型，使用所述最优化传染病传播模型预测传染病传播情况，本发明有效提高了新冠病毒传染模型参数预测的准确性。

Description

一种基于卡尔曼滤波辨识新冠病毒传播模型参数的方法

技术领域

本发明涉及新型冠状病毒传播模型预测技术领域，尤其涉及一种基于卡尔曼滤波理论辨识新型冠状病毒传播模型参数的方法。

背景技术

目前，新型冠状病毒传播模型大都是采用的仓室模型(compartmental model)。传染病仓室模型是传染动力学中最为常用的传染病模型，其中SIR模型是最基本最经典的模型。根据这次新冠病毒传播的特点，本专利提出了一种基于卡尔曼滤波理论辨识新型冠状病毒传播模型参数来预测传染病的演化情况。

已有技术存在忽略传播率随时间的变化情况，传播模型不考虑潜伏期，预测精度差，计算大，参数识别准确性低等缺点。

发明内容

本发明提出一种基于滤波理论辨识新型冠状病毒传播模型参数的方法，其特征在于，包括如下步骤：

步骤1获取新型冠状肺炎疫情的时间序列数据集；

步骤2建立传染病传播模型；

步骤3使用卡尔曼滤波算法补充所述传染病传播模型中的未知参数；

步骤4建立传染率变化模型；

步骤5根据传染率变化模型更新所述传染病传播模型得到最优化传染病传播模型；

步骤6使用所述最优化传染病传播模型预测传染病传播情况。

进一步地，步骤2所述的建立传染病传播模型，其实现步骤为：

步骤1获取发病地区的疫情数据集，即无症状感染者的数据、有症状感染者的数据、治愈者的数据以及死亡者的数据；

步骤2将所述发病地区人群分为六个类别，分别为：S代表易感者，E代表处于潜伏期的人群，A代表无症状感染者，I代表感染者(有症状的)，R代表治愈者，D代表死亡者；

步骤3传染病传播模型由六个常微分方程组成，描述了每个类别在每个阶段随时间的演化情况；

X(t)＝{S(t),E(t),A(t),I(t),R(t),D(t)}，代表各个类别每天的数据，代表各个类别每天的变化情况，也就是变化率，在这里变化率＝当天(t)人数-前一天(t-1)人数来计算。α表示易感者(S)被无症状患者(A)传染的概率，β表示易感者(S)被有症状患者(I)传染的概率，ε表示处在潜伏期(E)的人转化为无症状患者(A)的概率，θ表示处在潜伏期(E)的人转化为有症状患者(I)的概率，μ表示无症状患者(A)出现临床症状的概率，v表示无症状患者(A)恢复的概率，ρ表示有症状患者(I)恢复的概率，τ表示死亡率。

进一步地，步骤(3)所述的使用卡尔曼滤波算法补充所述传染病传播模型中的未知参数，其实现步骤为：

步骤1根据所述传染病传播模型得到未知参数：α，β，ε，θ，μ，v，ρ，τ

步骤2将所述传染病传播模型整理为：

其中X(t)_i＝{S(t),E(t),A(t),I(t),R(t),D(t)}，分别代表这六类类别，代表每类类别的变化率，f_ik(X₁,X₂,...X_n)表示第i个方程中的第k个函数，，即：i＝1时，i＝2时，/>i＝3时，/>i＝4时，/>i＝5时，/>i＝6时，f_6k＝I(t),k＝1。因为微分方程组里有重复的参数，为了简便计算，所以选取上述微分方程组中的(2),(4)和(5)表达式进行参数的估计，即利用i＝2,4,5时，/>的表达式进行参数估计，此时f_ik函数分别为：v_i表示均值为零的高斯白噪声，w_ik为当前微分方程所需要估计的系数；

步骤3使用卡尔曼滤波算法预测w_ik，即对应为α，β，ε，θ，μ，v，ρ，τ。

进一步地，步骤(4)所述的建立传染率变化模型，其实现为：

传染率变化模型为时间相关的单调递减的函数，具体表示为：

α(t)＝α₀e^-ωt (8)

β(t)＝β₀e^-ωt (9)

其中，ω∈(0,1)，α₀，β₀是α和β的初始感染率。

进一步地，步骤(5)所述的根据传染率变化模型更新所述传染病传播模型得到最优化传染病传播模型，其实现为：

根据所述传染率变化模型，更新所述传染病传播模型，更新后的传播模型可表示为：

因为获取的疫情数据是按天变化的，所以上述微分方程组是以天作为单位时间进行连续变化的。假设E＝2*(A+I),S＝N-E-A-I-R-D,N是该发病地区的总人口。其中，X(t)＝{S(t),E(t),A(t),I(t),R(t),D(t)}，代表各个类别每天的数据，代表各个类别每天的变化情况，也就是变化率，在这里变化率＝当天(t)人数-前一天(t-1)人数来计算；α表示易感者(S)被无症状患者(A)传染的概率，β表示易感者(S)被有症状患者(I)传染的概率，ε表示处在潜伏期(E)的人转化为无症状患者(A)的概率，θ表示处在潜伏期(E)的人转化为有症状患者(I)的概率，μ表示无症状患者(A)出现临床症状的概率，v表示无症状患者(A)恢复的概率，ρ表示有症状患者(I)恢复的概率，τ表示死亡率。

本发明基于一种基于卡尔曼滤波理论辨识新型冠状病毒传播模型参数的方法，可以有效提高传染病传播预测精度，解决传播模型不考虑潜伏期，预测精度差，计算大，参数识别准确性低等问题。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作一简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明的整体流程图；

图2为本发明卡尔曼识别传染病传播模型参数算法流程图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

图1为本发明的整体流程图，本发明提出一种基于卡尔曼滤波理论辨识新型冠状病毒传播模型参数的方法，其特征在于，包括如下步骤：

(1)获取新型冠状肺炎疫情的时间序列数据集；

(2)建立传染病传播模型；

(3)使用卡尔曼滤波算法补充所述传染病传播模型中的未知参数；

(4)建立传染率变化模型；

(5)根据传染率变化模型更新所述传染病传播模型得到最优化传染病传播模型；

(6)使用所述最优化传染病传播模型预测传染病传播情况；

步骤(2)所述的建立传染病传播模型，其实现步骤为：

(2a)获取发病地区的每日疫情报告数据，即无症状感染者的数据、有症状感染者的数据、治愈者的数据以及死亡者的数据；

(2b)将所述发病地区人群分为六个类别，分别为：S代表易感者，E代表处于潜伏期的人群，A代表无症状感染者，I代表感染者(有症状的)，R代表治愈者，D代表死亡者；

(2c)传染病传播模型由六个常微分方程组成，描述了每个类别在每个阶段随时间的演化情况。

因为获取的疫情数据是按天变化的，所以上述微分方程组是以天作为单位时间进行连续变化的。假设E＝2*(A+I),S＝N-E-A-I-R-D,N是该发病地区的总人口。其中，X(t)＝{S(t),E(t),A(t),I(t),R(t),D(t)}，代表各个类别每天的数据，代表各个类别每天的变化情况，也就是变化率，在这里变化率＝当天(t)人数-前一天(t-1)人数来计算。α表示易感者(S)被无症状患者(A)传染的概率，β表示易感者(S)被有症状患者(I)传染的概率，ε表示处在潜伏期(E)的人转化为无症状患者(A)的概率，θ表示处在潜伏期(E)的人转化为有症状患者(I)的概率，μ表示无症状患者(A)出现临床症状的概率，v表示无症状患者(A)恢复的概率，ρ表示有症状患者(I)恢复的概率，τ表示死亡率。

步骤(3)所述的使用卡尔曼滤波算法补充所述传染病传播模型中的未知参数，其实现步骤为：

(3a)根据所述传染病传播模型得到未知参数：α，β，ε，θ，μ，v，ρ，τ

(3b)将所述传染病传播模型整理为：

其中X(t)_i＝{S(t),E(t),A(t),I(t),R(t),D(t)}，分别代表这六类群体，代表每类群体的变化率，f_ik(X(t)₁,X(t)₂,...X(t)₇)表示第i个方程中的第k个函数，，即：i＝1时，/>i＝2时，/>i＝3时，/>i＝4时，/>i＝5时，/>i＝6时，f_6k＝I(t),k＝1。因为微分方程组里有重复的参数，为了简便计算，所以选取上述微分方程组中的(2),(4)和(5)表达式进行参数的估计，即利用i＝2,4,5时，/>的表达式进行参数估计，此时f_ik函数分别为：v_i表示均值为零的高斯白噪声，w_ik为当前微分方程所需要估计的系数；

(3c)使用卡尔曼滤波算法预测w_ik，即对应为α，β，ε，θ，μ，v，ρ，τ；

使用卡尔曼滤波算法预测w_ik，可见图2，下面的步骤中w_ik用w表示，其步骤为：

1.数据准备

假设一共观测到m天数据，将S、E、A、I、R和D六个类别的每日数据和各自的变化率整理好，形成一个m(行数)×12(列数)的矩阵x，在这个数组中前六列(1-6列)分别代表每类群体的每日数据；后六列(7-12列)各自对应自己的变化率。

2.初始化：

设定卡尔曼滤波的相关初始值。t＝0时刻的状态估计误差协方差权向量w系统噪声和量测噪声得相关矩阵Q和R，最大迭代次数N。

3.计算预测输出

y(t)＝w^T(t)X(t) (8)

X(t)就是指的f_ik(X(t)_1,X(t)₂,...X(t)₇)函数，因为只需要用微分方程组中(2),(4),(5)表达式来识别参数，X(t)也就是在i＝2,4,5时，f_ik的函数，分别为： 在i＝2时，对应矩阵x中：

X(t)是一个m*4的矩阵；i＝4时，对应矩阵x中，X(:,1)＝x(1:m,2)；X(:,2)＝x(1:m,3)；X(:,3)＝-x(1:m,4)；X(:,4)＝-x(1:m,4)；X(t)是一个m*4的矩阵，i＝5时，X(:,1)＝x(1:m,3)；X(:,2)＝x(1:m,4)；X(t)是一个m*2的矩阵。

4.预测误差输出

e(t)＝d(t)-y(t) (9)

d(t)就是就是{S(t),E(t),A(t),I(t),R(t),D(t)}六类群体中某一类的变化率，y(t)是卡尔曼滤波的预测值。

5.计算t时刻的卡尔曼增益

6.权向量更新

w(t+1)＝w(t)+K(t)e(t) (11)

7.更新状态估计误差协方差

P(t+1)＝P(t)-K(t)X^T(t)P(t)+Q (12)

8.t＝t+1，转到步骤3，迭代进行权向量更新，直至t＝N，输出辨识结果。

步骤(4)所述的建立传染率变化模型，其实现为：

传染率变化模型为时间相关的单调递减的函数，具体表示为：

α(t)＝α₀e^-ωt (13)

β(t)＝β₀e^-ωt (14)

其中，ω∈(0,1)，α₀，β₀是α和β的初始感染率，是卡尔曼滤波算法辨识的结果；

步骤(5)所述的根据传染率变化模型更新所述传染病传播模型得到最优化传染病传播模型，其实现为：

根据所述传染率变化模型，更新所述传染病传播模型，更新后的传播模型可表示为：

因为获取的疫情数据是按天变化的，所以上述微分方程组是以天作为单位时间进行连续变化的。假设E＝2*(A+I),S＝N-E-A-I-R-D,N是该发病地区的总人口。其中，X(t)＝{S(t),E(t),A(t),I(t),R(t),D(t)}，代表各个类别每天的数据，代表各个类别每天的变化情况，也就是变化率，在这里变化率＝当天(t)人数-前一天(t-1)人数来计算。α表示易感者(S)被无症状患者(A)传染的概率，β表示易感者(S)被有症状患者(I)传染的概率，ε表示处在潜伏期(E)的人转化为无症状患者(A)的概率，θ表示处在潜伏期(E)的人转化为有症状患者(I)的概率，μ表示无症状患者(A)出现临床症状的概率，v表示无症状患者(A)恢复的概率，ρ表示有症状患者(I)恢复的概率，τ表示死亡率。

步骤(6)所述使用所述最优化传染病传播模型预测传染病传播情况，其步骤为：

(6a)设置S、E、A、I、R、D的初始值以及预测时间T的长度；

(6b)调用matlab中的ode45函数求解微分方程组(4)；

(6c)得到一个T*6的输出数组Y，T代表行数，6代表列数，即{S，E，A，I，R，D}这六个类别。每一行代表的含义就是该天这六个类别的数量。

本发明基于一种基于卡尔曼滤波理论辨识新型冠状病毒传播模型参数的方法，可以有效提高传染病传播预测精度，解决传播模型不考虑潜伏期，预测精度差，计算大，参数识别准确性低等问题。

最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。

Claims (1)

1.一种基于卡尔曼滤波理论辨识新型冠状病毒传播模型参数的方法，其特征在于，包括如下步骤：

步骤1获取新型冠状肺炎疫情的时间序列数据集；

步骤2建立传染病传播模型；

所述的建立传染病传播模型，其实现步骤为：

步骤2-1获取发病地区的疫情数据集，即无症状感染者的数据、有症状感染者的数据、治愈者的数据以及死亡者的数据；

步骤2-2将所述发病地区人群分为六个类别，分别为：S代表易感者，E代表处于潜伏期的人群，A代表无症状感染者，I代表感染者，R代表治愈者，D代表死亡者；

步骤2-3传染病传播模型由六个常微分方程组成，描述了每个类别在每个阶段随时间的演化情况；

其中，X(t)＝{S(t),E(t),A(t),I(t),R(t),D(t)}，代表每个类别每天的数据，代表各个类别每天的变化情况，也就是变化率，在这里变化率＝当天人数-前一天人数来计算，α表示易感者S被无症状患者A传染的概率，β表示易感者S被有症状患者I传染的概率，ε表示处在潜伏期E的人转化为无症状患者A的概率，θ表示处在潜伏期E的人转化为有症状患者(I)的概率，μ表示无症状患者A出现临床症状的概率，v表示无症状患者A恢复的概率，ρ表示有症状患者I恢复的概率，τ表示死亡率；

步骤3使用卡尔曼滤波算法补充所述传染病传播模型中的未知参数；

步骤3所述的使用卡尔曼滤波算法补充所述传染病传播模型中的未知参数，其实现步骤为：

步骤3-1根据所述传染病传播模型得到未知参数：α，β，ε，θ，μ，v，ρ，τ；

步骤3-2将所述传染病传播模型整理为：

其中X(t)_i＝{S(t),E(t),A(t),I(t),R(t),D(t)},分别代表这六个类别，代表每个类别的变化率，f_ik(X(t)₁,X(t)₂,...X(t)₇)表示第i个方程中的第k个函数，即：i＝1时，时，/>时，/>时，时，/>时，f_6k＝I(t),k＝1；因为微分方程组里有重复的参数，为了简便计算，所以选取上述微分方程组中的(2),(4)和(5)表达式进行参数的估计，即利用i＝2,4,5时，/>的表达式进行参数估计，此时f_ik函数分别为：

v_i表示均值为零的高斯白噪声，w_ik为当前微分方程所需要估计的系数；

步骤3-3使用卡尔曼滤波算法预测w_ik，即对应为α，β，ε，θ，μ，v，ρ，τ；

步骤4建立传染率变化模型；

步骤4所述的建立传染率变化模型，其实现为：

传染率变化模型为时间相关的单调递减的函数，具体表示为：

α(t)＝α₀e^-ωt(8)

β(t)＝β₀e^-ωt(9)

其中，ω∈(0,1)，α₀，β₀分别是α和β的初始感染率；

步骤5根据传染率变化模型更新所述传染病传播模型得到最优化传染病传播模型；

步骤5所述的根据传染率变化模型更新所述传染病传播模型得到最优化传染病传播模型，其实现为：

根据所述传染率变化模型，更新所述传染病传播模型，更新后的传播模型可表示为：

因为获取的疫情数据是按天变化的，所以上述微分方程组是以天作为单位时间进行连续变化的；假设E＝2*(A+I),S＝N-E-A-I-R-D,N是该发病地区的总人口，其中，X(t)＝{S(t),E(t),A(t),I(t),R(t),D(t)}，代表各个类别每天的数据，代表各个类别每天的变化情况，也就是变化率，在这里变化率＝当天人数-前一天人数来计算；α表示易感者S被无症状患者A传染的概率，β表示易感者S被有症状患者I传染的概率，ε表示处在潜伏期E的人转化为无症状患者A的概率，θ表示处在潜伏期E的人转化为有症状患者I的概率，μ表示无症状患者A出现临床症状的概率，v表示无症状患者A恢复的概率，ρ表示有症状患者I恢复的概率，τ表示死亡率；

步骤6使用所述最优化传染病传播模型预测传染病传播情况。