用简单的数学模型，预测新型肺炎接下来将会如何发展

（Image by Gerd Altmann from Pixabay）

撰稿人：沙威 浙江大学

新型肺炎牵动全国人民的心，各方力量都在共同努力，最终打赢这场战役指日可待。中国科研人，在基因测序、疫苗研发、病毒溯源、模型预测等领域均做出了巨大努力，我们对这种新型病毒的科学认识也不断增强。

作为一个计算电磁学领域的研究者，唯一能做的是科普工作，谈谈新型肺炎的简单数值模型和参数辨识（反演）方法。于是，组织团队成员合作，2天内完成了简单的数学模型。感谢安徽大学谢国大（导师黄志祥，浙大交换学生）完成了SIR/SEIR模型数值求解，徐婷（浙江大学）完成简单早期确诊病例的指数拟合，黄成念（厦门大学，即将来浙大攻博）完成数据收集，袁帅（浙江大学）和我完成参数反演。

所有数理模型都有应用范围，本文谈到的模型，是传染病模型中最简单的一类。（1）模型假设是封闭系统（总人口数守恒），模型无法刻画交通流动等复杂网络问题。（2）模型采用的物理参数是时间不变的，而随着政府强力干预措施，参数实际上是时变的。（3）模型是宏观模型，无法刻画个体行为，如超级传播者等。但通过简单模型，可从动力学角度直观理解病毒传播的物理图景；并在现有官方数据基础上，对疾病传播做出科学预测。

假设感兴趣的人群总数是N，定义三类人群：易感者用S(t)表示，感染者用I(t)表示，康复者（包括死亡者）用R(t)表示。康复者有抗体，不会再被传染。无论任何时间，S(t)+I(t)+R(t)=N，即人数守恒。SIR 模型可以写成如下常微分方程组:

因为交叉乘积项IS存在，上述方程是非线性常微分方程组。把三个子方程相加，右侧等于0，即可得到人数守恒条件。β是一个极重要的参数，叫传染率。每天感染人数的增加，等于感染率（感染风险），感染者数量，易感者数量的乘积。这里因为交叉项数值很大，所以除以N，避免β太小。政府干预，口罩使用，自我在家隔离，都强烈影响β。γ是另一个重要参数，叫康复率，是一个感染者从感染到治愈的平均时间（平均患病时间）的倒数。再生数α=β/γ，是指一个患者在其患病期内平均传染的人数。该模型很像光学激光器中的速率方程，其三类人群可类比成三个原子能级，传染率、康复率类比于跃迁概率，人数守恒类比于粒子数守恒。

上述模型只要给定初值，就可用一般的数值方法求解。这里我们采用预测校正显式方法，以第一个子方程为例，其离散化形式为:

给定参数N=25000，β=0.32，γ=1/21，I(0)=3，R(0)=0，仿真结果如下：

可以看出，发病40天后，感染者数量达到顶峰，后期逐渐下降。如果用指数坐标，画出前40天感染者的预测人数图，你会发现在30天前，其增长方式都是指数速率的，后期增长速率变低，达到峰值。

通过这个模型，对新型肺炎的早期数据，即对公布的确诊人数可以做曲线拟合，进行科学预测。以前七天确诊人数（黑点），进行数据拟合，预测后三天的数据。发现红点真实数据和黑色预测曲线，吻合良好，但略微减少。说明政府干预和市民努力起到了作用，几天后指数预测模型将失效，确诊人数增加速率放缓。

3. SEIR 模型及修正 

很多传染病，具有潜伏期。因此，需要添加一类新型人群，即潜伏者E(t)。SEIR模型其对应的常微分方程组为:

这里，κ叫潜伏率，是平均潜伏时间的倒数。但是根据新型肺炎的报道，潜伏者也可以传染病毒给易感人群，因此上述经典模型，需要被修正如下，

这里，λ是潜伏者对易感人群的传染率。

4. 参数辨识与反演 

根据官方报道，目前可获得的数据，是确诊者人数和康复者（包括死亡者）人数，感染者人数可用确诊人数减去康复者（包括死亡者）人数。我们需要拟合6个参数 P={ S(0), E(0), β, λ, κ, γ}，即易感者初始值S(0)，潜伏者初始值E(0)，传染率β和λ，潜伏率κ，康复率γ。优化问题，可用最大误差最小化方法（不会造成某些拟合点误差大的问题），MATLAB优化工具箱求解。

这里，ID 和 RD 是实际公布的感染者和康复者（包括死亡者）数据。因为治愈人数太少，这里我们用相对误差作为判据。冯诺依曼说过：用四个参数我可以拟合出一头大象，而用五个参数我可以让它的鼻子晃。越多参数，让解的唯一性存疑。考虑到目前数据很少，从单纯研究和科普目的出发，只用简单SIR模型来拟合。取6天数据和7天数据，最大相对误差8.8%和10.4%，结果如下： 

可以看出，后面几天的误差明显高于前几天，说明人工干预效果显著，简单SIR模型不再适用。确诊人数悲观（6天数据）和乐观（7天数据）的数值分析如下：

本文所涉及的MATLAB程序和数据，仅供科研和教学使用。如有需要，请提供真实姓名和单位信息，并用单位邮箱发送至作者邮箱索取。邮箱: weisha@zju.edu.cn