新冠肺炎全国疫情趋势预测——经典SIR模型的改进探索

研究概述

全国新型冠状病毒2019-nCOV自武汉爆发以来，过去两个多月在全世界迅猛传播，对人类社会经济与工作生活带来严重的影响。在此期间，广大医学专家与医护工作者身先士卒奋战在抗击疫情的前线；科研学者同时也从自己研究领域出发，采用多种算法与模型，力求挖掘病毒传播及疫情发展的时空规律。1月23日，武汉采取“封城”措施严阻疫情向各地传播。在此之前，湖北省累计确诊病例为444例(武汉425例)，其他省市病例数量相对较少；在此之后，疫情累计确认数量全国日渐增加，其中湖北省疫情发展趋势十分严重(Fig.1)，各省市医疗团队也及时采取有效防疫措施，启动公共卫生事件一级响应来控制病毒进一步传播。于此同时，广大群众更加关注疫情发展趋势、增长速率及数量“拐点”，本文借助全国病例数据(1月23日-2月23日)，通过建立并改进经典SIR传染病模型，模拟预测不同情境下未来新冠肺炎传播情况，给各省市疫情防控工作提供有力参考。

Fig.1 2019-nCOV新冠肺炎累计确诊病例(2020年1月23日-2月23日)全球发展趋势图.

人口迁徙影响下，SIR模型结构再加工

动力学模型是对新冠肺炎疫情进行数学预测的重要方法之一，它由微分方程组、以及合理参数与初始值设定组成。经典传染病动力学模型主要包括SIR模型、SEIR模型：前者将人群分成易感染人群(Susceptible)、感染人群(Infectious)、康复人群(Recovered)；后者加入病毒处与潜伏期人群(Exposed)，表示该阶段人群已被病毒感染但未呈现明显症状。在给定最初感染人数前提下，通过病毒感染率、潜伏期、康复率等一系列参数作用影响，实现不同时期不同阶段人群数量增减。

自新冠肺炎疫情爆发以来，广大学者用SIR模型、SEIR模型对病毒传播进行大量模拟，分别给出未来不同时期疫情发展走势，这对我们充分理解2019-nCOV病毒传播机理提供有效帮助。借助一定时间内官方公布疾病确认，康复与死亡人数，以及春运期间人口迁徙数据，我们尝试建立新冠肺炎传播模型，如图：

变量I：研究指出2019-nCOV病毒潜伏同样具备传染性，这里我们选择SIR模型进行模拟。变量R：这里特指该城市每天从疫情感染到住院隔离人群数量(住院隔离、康复、死亡)。变量S：因各地区最初感染人群很大概率近期有到过武汉，我们用S1与S2分别表示该城市的本地与武汉流入易感染人群。四个变量初值设为(S10, S20, I0, R0)，其中本地易感染人群S10= N0 (N0表示当地户籍总人口)，武汉流入易感染人群S20=Nw-I0 (Nw表示1月1日-1月24日累计武汉流入总人口)。 

蒙特卡洛模拟下，SIR模型参数再拟合

不同于经典SIR模型，β1与β2分别表示城市的本地户籍人群与春运期间武汉流入人群疫情日感染率，其中β1<β2；γ表示患者从感染病毒到入院被隔离，即广义上疫情日隔离速率。根据截至1月23日，武汉累计确诊肺炎病例425例的研究结果表面病毒在患者潜伏期平均5.2天，患者从病毒暴露到入院接受隔离平均9.1天，则广义病毒感染周期5.2+9.1=14.3天，即日隔离速率为1/14.3=0/07。

不同于常规函数拟合，对β1、β2的估计我们采用马尔科夫链蒙特卡洛(MCMC)算法得到。简单来说，MCMC是一种常用采样方法，即通过在概率空间中随机采样以近似兴趣参数(Parameter of Interest)的后验分布。因篇幅有限，这里简单提及算法所涉及贝叶斯统计思想。贝叶斯统计描述事物发生的概率受到先验概率假设和事件发生的可能性的影响。如果先验概率(Prior)和似然概率(Likelihood)分布如下，我们可以得到后验概率(Posterior)的分布由两者概率“组合”计算得到(Fig.2)。

Fig.2 贝叶斯统计先验概率Prior、似然概率Likelihood、后验概率Posterior示意图[2]. 

MCMC算法涉及到的“蒙特卡洛”、“马尔科夫链”及“验收-拒绝抽样”可参考该文献[2]。在95%置信区间，我们用1月25日-2月22日累计确诊病例数作为训练样本，借助PyMC3 python脚本模块[3]，求解得到各城市感染率β1=exp(-at)，β2=exp(-at)，其中参数作为各城市调节系数。几个典型城市参数拟合如下(Fig.3)。 

Fig.3 典型城市新冠肺炎累计确诊病例与MCMC拟合曲线散点图1月24日-2月22日(武汉，孝感，北京，成都，昆明，大同). 

模型验证

我们用1月24日-2月22日累计确诊病例作为训练样本得到参数β1、β2、γ表达式，再用2月23日新增肺炎病例数作为模型Eq1模拟结果的验证信息。结果显示(Fig.4)，中国大陆除武汉外334个城市疫情模拟结果与新增确诊病例数量大致吻合(R2=0.74)，而采用经典SIR模型的肺炎疫情预测结果相对较差(R2=0.14)。

Fig.4 改进SIR模型(a)与经典SIR模型(b)对除武汉外全国2月23日新增(核酸检测)病例. 

模型预测

我们分别给出未来两种情景模式模拟结果：(a)情景1维持现有趋势，(b)情景2随着复工增加人与人接触。截至3月23日，中国大陆累计确诊病例在情景1预计91883例，在情景2预计110615例(95%置信区间)。可以看出大部分感染情况更多集中在华中地区与东部沿海地区，相对而言西部地区预测病例较少(Fig.5)。

Fig.5基于改进SIR模型预测3月23日中国大陆肺炎累计确诊病例: (a)情景1(b)情景2.  

新冠肺炎新增病例数量“拐点”预测，即对各地区疫情数量达到峰值的日期模拟对于延阻病毒的传播有一定的帮助。在情景1中，即现有防疫工作维持之下，预计在2月24日中国大陆除武汉外的城市中，88.6%已达到新增病例峰值，预计在3月17日中国大陆除武汉外的所有城市均已到达确诊病例峰值。在情景2中，即随着复工过程持续，人流回暖，多个城市病例“拐点”大致会延缓至2月底(Fig.6)。

各地新冠肺炎新增病例归零日期的预测也是此次研究工作之一。两种情景模拟结果均表明，武汉市新增病例数在3月23日还不能减少归零，即疾病还不能完全消除，接近武汉的城市新增病例数为零的日期也会推迟出现；其余城市的模拟结果，在情景1正常趋势下，93%的城市每日新增病例会在二月底基本保持为零，在情景2复工作用影响下，只有不到80%的城市每日新增病例会在二月底基本保持为零，并且湖北省仍有三座城市在3月23日后仍有新增病例的出现(Fig.6)。

Fig.6.每日新增肺炎病例数抵达峰值日期模拟：(a)情景1，(b)情景2；每日新增肺炎病例数等于零日期模拟：(c)情景1，(d)情景2.

结语

我们利用改进后SIR模型对中国大陆新冠肺炎发展趋势进行时空验证与模拟。通过2月23日结果对比，全国城市(武汉除外)模型拟合结果与当日新增病例大致相同(R2=0.74)，如果保持现有防疫措施，全国93%城市每日新增病例数量2月底可达到峰值；如果复工过程持续，2月底病例到达峰值的城市比例将不到80%。相比已有模型研究，该研究对经典SIR模型有如下情况改进：在模型结构方面，我们利用春运期间人口迁徙数据，在方程组中加入有武汉迁入的易感染人群部分，使得各城市的模拟结果充分体现出空间异质性；在模型参数方面，我们利用马尔科夫链蒙特卡洛模拟，分别给出各城市感染率拟合曲线，这对于准确把握各地区疫情发展走势提供了精细的参数化结果。另一方面，现有模型还有一定缺陷，例如武汉市2月13-15日突然新增临床病例数据未能很好在模型模拟显示。虽然相关部门已经严格限制各地交通来往，但随着复工展开，模型在一定程度上可以揭示出三月份新冠肺炎的发展规模，这对各部门防疫工作提供一定参考。

参考文献

[1] Wang Xia, Tang Sanyi, Chen Yong, Feng Xiaomei, Xiao Yanni, Xu Zongben. When will be the resumption of work in Wuhan and its surrounding areas during COVID-19 epidemic? A data-driven network modeling analysis. SCIENTIA SINICA Mathematica, 2020.2. 17.

[2] 简洁清晰解释马尔科夫链蒙特卡洛方法 AI火箭营 2019.3.16.

[3] Salvatier, J., Wiecki, T. V., Fonnesbeck, C., 2016. Probabilistic programming in Python using PyMC3. PeerJ Comput. Sci. https://doi.org/10.7717/peerj-cs.55.

[4] 丁香园·丁香医生新型冠状病毒感染 肺炎疫情实时动态

 http://ncov.dxy.cn/ncovh5/view/pneumonia，2020.

文章链接

相关论文近期发表于在线预印版杂志medRxiv，Xiaolin Zhu, Aiyin Zhang, Shuai Xu, Pengfei Jia, Xiaoyue Tan, Jiaqi Tian, Tao Wei, Zhenxian Quan, Jiali Yu, Spatially Explicit Modeling of 2019-nCoV Epidemic Trend based on Mobile Phone Data in Mainland China. medRxiv，2020. doi: https://doi.org/10.1101/2020.02.09.20021360

团队介绍

本研究主要由中国城市规划设计研究院学术信息中心(创新中心)联合香港理工大学土地测量和地理信息学系助理教授朱孝林研究团队完成。香港理工大学土地测量及地理资讯学系朱孝林老师的研究团队主要从事空间数据（卫星遥感数据，无人机数据，社会感知大数据）的处理和分析，主要研究方向包括遥感数据的自动化处理算法开发，地表定量信息反演，生态系统动态监测，以及突发事件时空建模与风险预测。参与该工作的成员博士生田家旗、徐帅、谭晓悦，硕士生张爱因。本研究具体细节可通过电子邮件联系xiaolin.zhu@polyu.edu.hk

文字撰写

贾鹏飞 中国城市规划设计研究院学术信息中心(创新中心)