基于阶段式SIR-F模型的新冠肺炎疫情评估及预测

凡友荣 杨 涛 孔华锋

1(公安部第三研究所 上海 201204) 2(武汉商学院 湖北 武汉 430056)

0 引 言

根据世界卫生组织公布的数据，截至2020年5月27日，全球新冠肺炎确诊病例已达548万，超过10万确诊人数的国家包括美国、巴西、俄罗斯、英国、意大利、法国、德国、土耳其、印度、伊朗、秘鲁。其中美国、巴西目前的每日新增确诊数超过1万人；意大利、法国等国家的疫情趋缓。由于印度存在人口密集、医疗条件差等问题，世界卫生组织表示，人类对抗新冠疫情能否取得决定性胜利，未来很大程度上取决于印度控制该病毒的能力。针对全球各国不同的疫情状况，分别进行现状分析及防控措施的影响研究至关重要。

SIR模型、SEIR(Susceptible-Exposed-Infected-Recovered)模型是传染病动力学研究的经典模型，目前已有多项基于SIR、SEIR模型的新冠疫情研究。尹楠[1]基于SIR模型进行了有限区域内新冠病毒传播的模拟分析。梅文娟等[2]将极限学习机与SIR模型相结合构建了极限IR实时预测模型。盛华雄等[3]将SIR模型与Logistic模型相结合，对不同防控措施进行了模拟和对比。喻孜等[4]使用基于时变参数的 SIR 模型预测得到了疫情拐点和最大确诊数。Zhou等[5]基于SEIR模型计算得到了COVID-19的再生数。林俊锋[6]在SEIR模型的基础上增加了“隐形传播者”，对隐形传播者的数量进行了预测。范如国等[7]基于SEIR模型对3种不同潜伏期下的疫情进行了拐点预测。

本文使用的SIR-F模型在SIR模型的基础上增加了未确诊的感染者(具有传染性)，这一点与SEIR模型相同；本文模型还将“移除者”(Removed)细分为“治愈者”(Recovered)和“死亡者”(Fatal)，从而能够进一步分析治愈人数、死亡人数的变化趋势。此外，考虑到在数据拟合的实践中存在过拟合的现象，并且实际的疫情发展过程分为多个阶段，每个阶段的趋势特征各不相同，因此本文先进行阶段划分，再将疫情数据分阶段地进行拟合分析，从而提高分析结果的准确性。

1 技术路线

本文的技术路线如图1所示。首先计算各国新冠肺炎确诊数量的增长因子，根据该因子将所有国家分为3个增长类型：“停滞”式、“爆发”式、“波动”式；然后在“停滞”增长态势的国家中选择典型例子进行案例分析，由于意大利是全球重点疫情国家之一，且其防控工作已取得显著成效，因此选择意大利为例，对其疫情发展过程和防控措施进行深入分析；最后针对“爆发”式增长的国家，从中识别目前爆发态势最严峻的国家，并针对其进行疫情的深入分析，提出防控措施建议。

在案例分析中，首先使用Python调用Prophet的“突变点自动监测”功能对疫情数据进行阶段划分，可用于疫情数据突变点的自动识别。基于SIR-F模型进行分阶段的数据拟合、参数估计，对比不同阶段的参数，从而发现各阶段疫情的变化情况，然后结合防控措施的实施时间，研究防控措施对疫情发展产生的影响。对于“爆发”式增长的国家，在模型拟合的基础上进行疫情预测，并且基于不同的感染系数进行预测比对，从而模拟不同强度防控措施下的疫情发展趋势。

2 模型构建

2.1 增长态势识别

本文通过增长因子来评估 “新冠肺炎”确诊人数增长的态势，增长因子的计算公式如下：

(1)

式中：GF为增长因子；ΔCn为第n天确诊人数增加的数量。计算每个国家近7天的GF，进行增长情况的分类：

(2)

GC为增长态势所属类别，GC=1代表“爆发”的增长态势，GC=-1代表“停滞”的增长态势，GC=0代表“波动”的增长态势。

2.2 SIR-F模型

由于COVID-19疫情比较特殊，有一部分感染者在死亡前仍未正式确诊，将这类情况增加到经典的SIR模型中，构建SIR-F模型。SIR-F模型中从未感染者转变为感染者、治愈者、死亡者的过程如图 2 所示。

图2 SIR-F模型中感染者的转化过程

图2中，S代表未感染者，S*代表未确诊的感染者，I代表感染者，R代表治愈者，F代表死亡者，α1为S*的死亡系数，α2为I的死亡系数，β为感染系数，γ为I的治愈系数。累计确诊数量C=I+R+F。

根据上述模型设计，构建其动力学方程如下：

(3)

(4)

(5)

(6)

式中：N为总的人口数量，N=S+I+R+F；T为疫情开始的总时长。为求解方便，将上述函数方程进行无量纲处理，设(S,I,R,F)=N×(x,y,z,w)，(T,α1,α2,β,γ)=(τt,θ,τ-1k,τ-1ρ,τ-1σ)，得到以下微分方程组：

(7)

在Python(版本3.7.6)环境下，使用超参数自动优化框架Optuna对上述微分方程进行求解，从而得到各个参数(θ,k,ρ,σ)的估计值。

2.3 防控措施的影响函数

各国政府通过疫情防控措施改变民众的行为，从而控制疫情的发展。本文引入感染系数β受民众行为影响的函数，如式(8)所示。设感染系数β的影响因子包括：每分钟内未感染者与感染者在外出过程中接触的次数c，每周内未感染者外出的天数g，唾液中病毒存在的可能性v，有效佩戴口罩的人数比例m，洗手降低的感染概率we，其他影响因素δ。

(8)

式(8)与SIR-F模型相结合，能够模拟居民的外出、佩戴口罩、洗手等行为对疫情发展的影响。

3 模拟仿真

3.1 数据来源

本文基于的全球疫情数据来源于约翰·霍普金斯大学的系统科学与工程中心(https://systems.jhu.edu/ )，包含了全球217个国家在110天内(2020年1月22日至2020年5月11日)的累计确诊人数、现存确诊人数、死亡人数、治愈人数，总数据量为23 804行。各国防控措施数据来源于牛津大学的流行病预测研究组织(http://epidemicforecasting.org/)，包含了96个国家在2019年12月18日至2020年4月15日期间采取的管控措施，总数据量为1 633行。

3.2 全球疫情的增长态势分析

将全球疫情数据导入式(1)和式(2)，计算得到96个国家的疫情增长态势所属类别GC，结果示例如表1所示。

表1 全球各国疫情增长态势的分类结果

处于“爆发”增长态势的包括印度、俄国、爱沙尼亚、巴西、肯尼亚等53个国家。将这类国家的疫情数据聚合，增长趋势如图3所示。

图3 “爆发”式增长国家的疫情数据变化趋势

处于“停滞”增长态势的包括意大利、日本、利比亚、安哥拉、尼加拉瓜、越南等28个国家。将这类国家的确诊数据聚合，增长趋势如图4所示。

图4 “停滞”式增长国家的疫情数据变化趋势

处于“波动”增长态势的包括加拿大、以色列、美国、英国等138个国家。将这类国家的确诊数据聚合，增长趋势如图5所示。

图5 “波动”式增长国家的疫情数据变化趋势

观察图3至图5可知，“爆发”式呈现迅速上升的趋势；“停滞”式的增长态势已基本得到了控制；“波动”式的未来发展方向尚不明确。“爆发”式和“波动”式增长的国家是全球疫情控制的重点，因此有必要研究“停滞”式增长的国家所采取的防控措施，为其他国家的疫情防控提供有效建议。

3.3 “停滞”式增长案例分析

1) 意大利疫情趋势分析。意大利目前的疫情态势如图6所示。根据式(1)计算得到意大利确诊人数的增长因子变化态势如图7所示，其增长因子GF<1的持续天数为18天。

图6 意大利疫情数据趋势图

图7 意大利增长因子变化趋势图

使用Prophet 对意大利的累计确诊数据进行突变点的识别，如图8所示。

图8 意大利累计确诊数据的突变点识别

由图8可知，存在3个突变点：“2020-03-17”、“2020-04-05”、“2020-04-25”，则意大利的疫情发展过程可以划分为4个阶段，如表2所示。每个阶段的变化规律有明显的变化，因此下文将分阶段地进行拟合分析。

表2 意大利疫情发展的4个阶段

2) SIR-F模型拟合及参数估计。将4个阶段的真实数据分别导入SIR-F模型进行拟合，拟合效果如图9所示。

图9 意大利4个阶段感染人数的拟合效果

对式(7)进行求解，计算得到各阶段SIR-F模型的参数估计值，如表3所示，由该表可知参数ρ下降了，即感染系数β下降了，并且有效再生数RT明显下降，这说明意大利的疫情控制已取得了一定成效。

表3 意大利4个阶段的参数估计值

3) 感染率降低的原因分析。意大利在3月采取了多项限制居民出行的管控措施，如表4所示。

表4 意大利3月采取的限制出行措施

考虑到政府的管控措施从开始执行到产生实际的疫情控制效果存在时间延迟，因此认为第1阶段为实施管控措施之前，第2、第3、第4阶段为实施管控措施之后。因此本文假设居民在第1阶段的出行情况与平时相同，设意大利民众平时每周出行5.6天，即g1=5.6。

根据2020年3月13日报道，由于出行限制、公共活动场所的关闭，使得居民相互之间的潜在接触次数降低了19%[8]，因此将采取管控措施之前的c取值为1，采取管控措施之后的c取值为0.81。

设出行限制措施只影响了式(8)中的g和c，其他值不变，得到：

(9)

将g1=5.6，c1=1，ci=0.81(i=2,3,4)，以及表3中各阶段的ρ值代入式(9)，计算得到g2=2.92，g3=1.06，g4=0.51。可见民众每周的出行天数已经下降到了0.51天。

上述分析过程从实际的疫情数据反推得到了实施管控后的每周出行天数，可知民众的每周出行天数减少了90.9%，使得感染系数降低了92.6%，证明了限制出行的管控措施对新冠病毒的扩散起到了明显的抑制作用；同时证明了本文的分析模型能够较灵敏地模拟现实中防疫管控下新冠肺炎的传播情况。

除了上述管控措施，意大利的疫情防控和诊疗工作均参照了中国方案，这也是意大利疫情控制效果较好的重要原因。意大利从5月4日已进入抗疫和恢复经济的并行阶段，新闻显示62%的意大利人对于复工感到害怕，因此出行习惯的恢复还需要时间。若出行量逐渐增加，应从其他方面控制感染率，例如提高戴口罩、勤洗手习惯的普及率等，从而防止疫情再次反弹。

3.4 “爆发”式增长案例分析

印度的增长因子GF>1的持续天数为全球最长(66天)，其增长态势极其严峻，因此下文将印度作为“爆发”式增长的典型例子进行深入分析。

1) 印度疫情趋势分析。印度目前的疫情态势如图10、图11所示。

图10 印度疫情数据趋势图

图11 印度增长因子变化趋势图

对印度的累计确诊数据进行突变点的识别，如图12所示，存在2个突变点(“2020-04-03”、“2020-04-28”)，则印度的疫情发展过程可以划分为3个阶段，如表5所示。

图12 印度累计确诊数据的突变点识别

表5 印度疫情发展的3个阶段

2) SIR-F模型参数估计。与意大利案例的分析方法相同，将印度各阶段的疫情数据分别导入SIR-F模型，参数估计结果如表6所示，可以看出参数ρ和RT存在一定程度的下降，但第3阶段的值仍较高，说明目前印度的疫情扩散仍然较为迅速。

表6 印度3个阶段的参数估计值

3) 印度疫情预测。以第3阶段的参数值进行疫情预测，结果如图13所示，印度感染数量将在2020年10月16日达到最高峰，最高值为3.84亿人。根据纽约时报3月27日发布的一篇有关印度疫情的文章，作者Laxminarayan教授预测印度感染人数将达到3至5亿，这与本文的计算结果基本符合。

图13 目前管控措施下的印度疫情预测

若印度采取更强有力的防控措施，例如为印度贫民地区提供更多口罩、酒精等医疗物资，提高核酸检测的覆盖率，更精准地进行病患隔离等，将能有效降低感染系数。假设印度的ρ值从目前的0.041降低为0.026，再次进行疫情预测，结果如图14所示。在ρ为0.026的假设下，印度感染数量的最高值降为2.01亿人，最高值到来的时间推迟到了2021年3月7日。

图14 假设ρ值降至0.026的疫情预测

通过上述对比分析可知，若印度增强防控措施的力度，将感染系数在目前的基础上降低36.59%，感染人数的最高值即可降低47.66%，到达最高感染人数的日期将推迟142天。

4) 印度疫情防控的现状分析及建议。由上述模拟分析可知，印度的疫情防控将是一个漫长的过程，且感染人数十分巨大，这与印度目前的实际情况相符。虽然印度从3月25日实施了全国封锁，采取了较严格的管控措施，但是由于印度存在人口基数大且密度大、核酸检测率不高、医疗资源有限、无症状感染者多等问题，导致疫情防护工作很难快速取得成效。此外，从5月1日开始印度开通专列帮助农民工返乡，其中不少人在抵达目的地后被检测出新冠肺炎阳性，因此此次返乡活动提高了病毒传播的风险。以印度德里为例，在近1万人的确诊病例中，超过75%是无症状或轻症感染者，且目前整体的核酸检测率未达到0.2%，在这样的情况下，放开公共交通会有极大的风险。

综上，从控制疫情的角度考虑，建议印度在感染人数的峰值(2020年10月16日)到来前控制公共交通的恢复比例；出现社区传染的贫民窟是疫情防控的关键；必须继续执行严格的限行措施才能防止疫情更加恶化。此外，防疫消毒、核酸检测、精准隔离、医疗条件等均是阻止印度成为全球疫情“震中”地区的必要条件。

4 结 语

本文将SIR-F模型与疫情增长态势的分类、防控措施影响分析方法相结合，以真实疫情数据为基础，使用Python实现模拟仿真，实现了对疫情的趋势分析、重点疫情地区的识别、疫情阶段划分、防控措施对疫情的影响分析、疫情预测分析。以意大利为“停滞”增长态势的例子，以印度为“爆发”增长态势的例子，对其疫情的发展过程进行了分阶段的拟合分析，实验证明，分阶段的拟合方法能比较准确地刻画疫情数据随时间的变化规律，并且在逻辑上更符合客观事实。本文分析方法同样适用于其他国家或地区，可以为新冠肺炎等类似传染病疫情防控提供较大的参考和指导价值。

实验结果证明了意大利的出行限制措施对疫情已经起到了有效作用。印度的确诊数量虽然低于美国等多个国家，但其爆发式增长状况极其严峻，预测得到其感染人数最高值为3.84亿人；若印度增强防控力度，使得感染系数降低36.59%，则其感染人数最高值将降低47.66%。最后，结合印度的最新情况，本文建议印度应重点加强贫民窟的疫情防控力度，在2020年10月16日之前继续控制公共交通，并结合其他精准防控手段阻止其疫情继续恶化。