1.   引　言

新冠病毒肆虐全球已经三年. 认识新冠疫情传播规律, 是预测和制定应对策略的基础, 是社会和政府共同关注的焦点, 也是科技工作者不懈努力要达成的目标. 2022年12月7日国务院联防联控机制综合组在坚持“动态清零”总方针下, 对防控措施进行了优化^{[1, 2]}, 此后, 一波新的规模巨大的疫情在各地陆续爆发并蔓延. 中国大陆拥有多个人口千万的大型城市, 由于人口众多而且集中, 大型城市面临更加严峻的挑战. 给出大城市疫情发展趋势和关键的时间节点, 成为迫切需要解决的问题^[3].

新冠病毒不断变异和扩散, 使得疫情具有长期性和地域差异性, 每次爆发也有新的特点. 新冠病毒疫情有复杂的时空特征^[4-6], 病例数时间序列呈现出依赖于管控政策的混沌特征^[7]. 人们采用回归模型^[8-16]、疾病传播模型^[3,17-21]、机器学习^[22-30] 等对一些国家或全球疫情进行了预测, 取得了一些有指导价值的结果. 比如, 采用长短记忆神经网络(LSTM)对病例数预测, 标准差可以达到$45.72$, 对疫情应对准备工作起到了指导作用^[25].

防控措施优化后, 无法获得完整准确的疫情统计数据, 采用大数据技术获得相关数据进行预测, 成为目前人们努力的方向, 例如百度APP上线的“百度疫情指数”^[31]试图给出实时的预测. 但是, 这些预测缺乏必要的理论基础. 本文综合考虑疾病传播的正反馈作用、疫情和人们出行行为之间的负反馈作用、以及阳性感染者转阴后重新被感染等因素, 给出了一个改进的虫口模型, 用以描述疫情传播动力学过程. 采用疫情初期少量疫情相关数据, 该模型实现了对特定时间和地域疫情的预测.

具体来说, 该模型很好地复现了上海市2022年3月1日到6月28日新冠病毒阳性感染者人数曲线. 利用百度指数^[32]查询了2022年12月1日到20日关键词“布洛芬”(ibuprofen) 的搜索量, 采用给出的模型对这一搜索量的走势进行了回归. 假设这一走势与每日新增被感染人数趋势一致, 对上海2022年12月爆发的疫情进行了预测. 结果表明, 如果病毒阳性感染者转阴后有$80\%$的人不再被感染, 这波疫情将在2023 年1月10 左右每日新增达到峰值, 在2023年3月7日恢复到2022年12月1日的水平. 如果病毒阳性感染者转阴后不再被病毒感染, 这波疫情将在2022年12月31日左右每日新增达到峰值, 每日新增在这一峰值前后各一周内处于持续高位, 在2023年2月10日恢复到2022年12月1日的水平.

2.   基于出行行为调整的疫情传播模型和材料

2.1   基于出行行为调整的疫情传播动力学模型

考虑一个离散映射模型:

其中, $t = 1, 2, \cdots$是迭代时间步. $X(t)$是第$t$时间步内人群中新增加的阳性感染者人数(确诊或者无症状感染者), $\lambda_0$是人们出行完全不受疫情影响下, 每一时间步内新增被感染者数目与人们接触次数的比率, 也就是病毒极限传播速度. 这一模型基于如下几点考虑:

1)在$t + 1$时间步内新增加阳性感染者人数$X(t + 1)$正比于$t$时间步内新增阳性感染人数$X(t)$, 也正比于人群中潜在的可感染人数$P(t)$.

2)阳性感染者转阴后, 重新被感染的概率明显降低. 用参数$\gamma$表示阳性感染者转阴后不再重复感染的比率, 在$t$时间步内潜在的可被感染人群总数就是$A(t)$ (见方程(4)). 第$t$个时间步结束之后, $A(t)$中$X(t)$变为阳性感染者, 剩下的$P(t) =A(t) - X(t)$就是第$t + 1$时间步内可感染人数(见方程(2)). 综合考虑事实1)和2), 可以得到(1)式.

3)随着每日新增阳性感染者人数的增加, 人们会降低出行次数. 用$\lambda(t)$表示人们出行的衰减情况, 其中$\alpha$ 描述衰减速率. 由于$X(t) < 1$, $\alpha$ 越小, 意味着人们出行越谨慎.

4)关注上海、北京、广州等大型城市, 以便模型中涉及到的统计假设能很好地得到满足.

本模型是在虫口模型(logistic model)基础上建立的, 称为改进的虫口模型(modified logistic model). 需要强调的是, 本文是对中国大陆城市某一波疫情进行预测. 在上海2022年3月到6月的疫情中, 实行了严格的“动态清零”策略, 新冠病毒阳性感染者能够在很短的时间内被发现, 其与未感染者之间的接触传染途径也被有效阻断, 因此每日新增阳性感染者人数代表了可感染人群接触的病毒源的规模. 在上海2022年12月以来的这波疫情中, 人们则是以发热等生理症状和指标作为被感染的判据. 人们延续了3年来养成的对疫情高度重视、积极应对的习惯, 感染者一旦被确认, 其与未被感染者接触传播病毒的行为也会被本人及其家庭成员等严格控制. 而这次疫情的特点是, 每个家庭内一旦出现有症状的感染者, 其他成员也往往会在1—2 天内出现症状, 这可以认为是潜在传播源的传播时间. 每日新增人数也能代表可感染人群接触到的病毒源的规模. 因此, 在模型中每日新增正比于上一时间步新增, 而不是多个时间步内感染者人数之和. 模型需要根据每次疫情爆发和传播特点进行相应的修改.

2.2   数　据

为了检验模型的的正确性, 我们收集了上海市2022年3月1日到6月28日新冠病毒阳性感染者数据. 在这一时期, 随着疫情的爆发, 上海市围绕“动态清零”总方针实行了一系列严格的防控措施, 包括全民核酸筛查、建设大量方舱对阳性感染者和密接人员进行集中管理、对人群的活动范围进行严格限制等. 考虑到无症状感染者中只有很小的比例被确诊为新冠肺炎, 意味着这一部分人所处的病毒感染源可能会具有更强的传染能力, 因此把每日新增无症状和确诊人数之和作为每日新增新冠病毒阳性感染者人数, 忽略重复计数问题. 该数据由上海市人民政府办公厅在“上海发布”微信公众号每日向公众发布^[33], 如图1所示.

2022年12月初以来, 由于疫情防控措施的优化^[1,2], 多地迎来新一波疫情, 给社会运行和医疗资源带来前所未有的挑战. 同时, 原来建立起来的每日新增病毒阳性感染者追踪系统也不再有效运转, 这为疫情预测带来了特殊困难. 我们采用百度指数网站上提供的关键词搜索工具^[32], 收集了自2022年12月1日到20日每日“布洛芬”的搜索量. 上海2022年12月这波疫情爆发之前, 北京等地已经陆续爆发大规模的新冠疫情, “布洛芬”已经被普遍认为是缓解阳性感染者症状的有效药物之一; 并且人们根据自己的症状自行判断服用这一药物的时间和剂量, 不再遵从“遵照医嘱服药”原则, 因此, 一个合理的假设是“布洛芬”的搜索量的变化趋势能够很好地反映疫情发展趋势. 另外, 在广州、北京、上海等多地陆续爆发的疫情, 是在2022 年12月防控措施优化后很短的时间内集中爆发的, 另一个合理的假设是, 各地的疫情传播速度等特征应该一致. 对上海这波疫情的预测, 正是利用了这两个假设.

3.   结　果

模型中需要确定的参数有$\lambda_0, \alpha, \gamma$和每个时间步长对应的实际天数$d_{\rm eff}$. 其中$\gamma$可依据疫情中的经验在一定范围内调节. 计算中取$\gamma = 0.8$和$1$, 也就是考虑80%的阳性感染者转阴后不再重新感染和一旦阳性感染者转阴后就不再感染两种情况. 实证数据和模型计算表明, 依据阳性感染者每日新增人数的变化, 疫情的发展可以划分为三个阶段. 在疫情爆发初期, 也就是阳性感染者成规模地出现, 但是其数量远没有达到峰值这段时间, 阳性感染者人数随着时间几何指数上升$\sim {\rm{exp}}(\beta_\uparrow t)$, 称为几何指数增长期; 在疫情末期, 也就是阳性感染者人数经过峰值并明显回落, 直到疫情结束, 阳性感染者人数随着时间几何指数衰减$\sim {\rm{exp}}(-\beta_\downarrow t)$, 称为几何指数衰退期. 在几何指数增长期结束和几何指数衰减期开始之前这段时间, 阳性感染者人数增长/衰减速度明显减小, 但是阳性感染者人数保持在高位, 并出现峰值, 称为高发期. 增长期和衰减期的速率会随着模型中每个时间步代表的现实中的天数同时改变, 但是其比值${\beta_\uparrow}/{\beta_\downarrow}$ 不变, 也与取定的初始感染者人数$X(1)$ 无关. 考虑到这一特点, 在重现已经有的疫情数据的时候, 我们在$[0, 4]$范围调试$\lambda_0$, 同时在$[a, 10]$范围调试$\alpha$, 得到与实证数据相同的增长/衰减速率比${\beta_\uparrow}/{\beta_\downarrow}$, 其中$a>0$. 再根据疫情实证数据的真实时间跨度和模型计算得到的疫情横跨的时间步长数, 得到每个时间步长对应的实际天数$d_{\rm eff}$. 在进行预测的时候, 采用相似的方法, 调整参数$\lambda_0$和$\alpha$得到疫情相关数据的几何指数增长趋势(本文是对布洛芬的搜索量进行回归), 然后采用得到的参数进行预测. $\beta_\uparrow, \beta_\downarrow$和$d_{\rm eff}$这三个参数表征了整个曲线的特征, 确定了峰值出现时间和大小、疫情持续时间等我们感兴趣的量.

图 1  复现中国上海市2022年3月1日至6月28日新冠病毒感染阳性病例数

Fig. 1.  Reproduction of the infected numbers for COVID-19 in the duration from March 1 to June 28, 2022 in Shanghai, China

下载: 全尺寸图片 幻灯片

在预测的时候, 还需要确定疫情的起始点. 考察对“布洛芬”的搜索量会发现, 在2022年12月1日之前搜索量围绕一个较低的水平线上下涨落, 之后搜索量开始按照几何指数增长. 上海市规模性疫情也出现在12月1日之后. 因此, 把“布洛芬”搜索量按照几何指数增长的起始点默认为上海疫情的起始点. 显然, 这两个起始点并不重合, 而是有一个小的差异. 真正的疫情起始点, 只能采用直接的阳性感染者人数数据才能得到. 我们的计算结果发现, 这一起始点并不影响我们感兴趣的几何指数增长期、高发期和几何指数衰退期跨越的时长, 以及阳性感染者人数峰值在疫情持续时间段中出现的位置及其大小(详见“结果”部分关于图5的讨论).

理论模型能够很好地复现2022年3月1日至6月28日上海新冠病毒阳性感染者人数曲线 (见图1). 当$\gamma = 0.8$时(黑色实线), $(\beta_\uparrow, \beta_\downarrow, d_{\rm eff}) = (0.2070, \, 0.1287,\, {105}/{190}) = (0.2070,\, 0.1287,\, 0.5526)$, $(\lambda_0, \alpha) = (1.13, 0.9)$. 当$\gamma = 1$时(红色实线), $(\beta_\uparrow, \beta_\downarrow, d_{\rm eff}) = (0.2128, 0.1738, \dfrac{105}{183}) = (0.2128, 0.1738, 0.5738)$, $(\lambda_0, \alpha) = (1.13, 0.8)$. 这意味着人们的出行处于较低水平, 并且随着政府发布的每日新增病毒阳性感染者人数增加, 会大幅度地降低出行水平.

图2给出了理论模型对2022年12月1日到20日“布洛芬”搜索量的回归结果(黑色实线). 假设这一回归得到的参数值反映了这波上海疫情的特征, 也就是疫情发展曲线和“布洛芬”搜索量回归曲线形状相同, 只需要在垂直方向移动. 当$\gamma = 0.8$和$1.0$时分别采用$d_{\rm eff} = 0.5526$和0.5738换算模型中单位时间步实际对应的时间, 得到对上海这波疫情的预测曲线. 模型参数值是$(\lambda_0, \alpha) = (1.17, 0.7)$, 也就是说, 虽然疫情防控措施得到了优化, 面对日益严峻的疫情, 人们延续了与严格管控期间相似的出行调整策略.

图 2  预测中国上海市2022年12月1日开始爆发的新冠疫病毒感染阳性病例数

Fig. 2.  Prediction of the infected numbers for COVID-19 starting from December 1 in Shanghai, China

下载: 全尺寸图片 幻灯片

当$\gamma = 0.8$, 也就是80%阳性感染者转阴后不再感染时(红色实线), 模型参数是$(\beta_\uparrow, \beta_\downarrow, d_{\rm eff})= (0.2700,\; 0.07589,\; 0.5526)$. 与上海3—6月份的疫情相比, 这波疫情爆发后上升速度明显更快, 而衰退速度明显变慢. 2022年12月1日到2023 年1月1日是几何指数增长阶段. 从2023年1月2日开始进入高发期, 在1月10日左右达到顶峰, 高发期持续到1月18号结束, 历时17天. 随后将进入几何指数衰退期, 到2023年3月7日恢复到2022年12月初的疫情水平. 当$\gamma = 1$, 也就是阳性感染者转阴后不再感染时(蓝色实线), 模型参数是$(\beta_\uparrow,\; \beta_\downarrow, \;d_{\rm eff}) = (0.2700, 0.1891, 0.5738)$. 与上海3—6月份的疫情相比, 这波疫情爆发后上升速度明显更快, 而衰退速度与之前相近. 2022年12月1日到2022年12月24日是指数增长期. 从2022年12月25日开始进入高发期, 在12月31日左右达到顶峰, 高发期持续到2023年1月8号结束, 历时15天. 随后进入几何指数衰退期, 到2023年2月10日恢复到2022年12月初的疫情水平.

需要指出的是, 由于2022年12月以来没有可靠的新冠病毒阳性感染者人数数据, 这里的每日新增只具有相对的意义, 并不代表实际预测人数. 在图2中, 出于视觉效果的考虑, 这些数据乘了一个适当的数, 使得曲线在垂直方向移动. 另外, 搜集到的“布洛芬”搜索量从12月1日到20日, 共20个数据. 开始的时候搜索量几何指数增长, 后边几天的数据明显偏离几何指数增长. 可能的原因是, 当服用“布洛芬”药物的方法在人群中大量普及之后, 人们可以靠相互交流(如家庭内部人员之间的交流)来正确服用药物, 其搜索量应该不再能体现新增人数. 因此, 我们只采用了几何指数增长部分确定模型参数.

图3和图4分别给出了$\gamma = 0.5$和$0.8$时, $\lambda_0$和$\alpha$ (人们出行行为水平及其随疫情变化)对疫情发展的影响. 可以发现$\lambda_0$和$\gamma$会决定几何指数增长期阳性感染者人数增加的速度($\beta_\uparrow$)、阳性感染者人数峰值的大小和出现的位置. 人们只有在每日新增阳性感染者人数很多的时候才考虑减少出行, 也就是当$\alpha$取小值的时候, 阳性感染者人数峰值较大而衰退期较短. 当人们的出行强烈地依赖每日新增阳性感染者人数, 也就是当$\alpha$取大值的时候, 阳性感染者人数峰值较小, 但是疫情会持续更长的时间. 一个有意思的结论是, 人们出行策略的变动($\alpha$)不影响每日新增阳性感染者人数峰值出现的时间.

图 3  出行行为参数对疫情演化的影响($\gamma = 0.5$)

Fig. 3.  Impact of travel behavior on epidemic process ($\gamma = 0.5$)

下载: 全尺寸图片 幻灯片

图 4  出行行为参数对疫情演化的影响($\gamma = 0.8$)

Fig. 4.  Impact of travel behavior on epidemic process ($\gamma = 0.8$)

下载: 全尺寸图片 幻灯片

在疫情爆发后, 地铁出行客流量也被广泛地用于对疫情发展状况的估计. 地铁系统人员密集、环境封闭、流动性强, 是新冠病毒扩散最畅通的渠道之一, 也是乘客被感染高风险场所. 一个合理的推论是, 地铁乘客将是疫情爆发后首先被感染的人群, 地铁客流量的锐减、低谷以及增加, 会比整个人群病毒阳性感染者人数的快速增加、峰值以及衰减等提前到来. 图5(a)为$\gamma = 1$时的阳性感染者人数预测曲线(与图2中的蓝色线相同, 采用线性坐标), 用虚线画出了每日新增阳性感染者人数峰值的${1}/{2}$位置. 图5(b)为地铁乘客流量曲线, 并把这一曲线向右平移了$8$天, 用虚线给出了疫情初期乘客人数最多$N_{\rm{max1}}$, 疫情中期乘客人数达到谷底$N_ {\rm{min}}$和疫情后期乘客人数逐渐增加达到新的峰值$N_ {\rm{max2}}$的位置, 以及乘客人数衰减一半$(N_{\rm{min}}\;+ N_{\rm{max1}})/2$和从谷底增加一半$(N_{\rm{min}}+N_{\rm{max2}})/2$的位置. 发现阳性感染者人数预测曲线和地铁客流量变化曲线的这些特征能够很好地符合. 图5(b)中红色实心圆为每日上海地铁客流量, 实黑线为去除掉周六周日和节假日(2023年1月1日和2日为元旦假期)的上海地铁客流量.

图 5  上海阳性感染者预测曲线与地铁乘客数量曲线比较　(a)上海2022年12月到2023年2月新冠阳性每日新增预测; (b)地铁乘客人数, 曲线向右移动了8天

Fig. 5.  Comparison of the prediction of infected number per day with the number of metro passengers in Shanghai China: (a) Prediction of COVID-19 virus infected numbers in the duration from Dec. 2022 to Feb. 2023 in Shanghai, China; (b) metro passenger numbers, where the time has been shifted rightly 8 days

下载: 全尺寸图片 幻灯片

预测曲线和地铁出行流量曲线在时间上有$8$天的差异, 除了乘坐地铁人群比整个社会人群面临更高的感染风险外, 另一个可能的来源是布洛芬搜索量的几何指数增长开始时间与阳性感染者人数几何指数增长的起始点之间的差. 实际上, 只有疫情达到一定规模并且自行使用布洛芬药物缓解症状被普遍认可, 搜索量才会随着阳性感染者人数的增加而增加.

因此, 地铁客流量、对“布洛芬”药物的搜索量等间接数据, 能帮助我们取得疫情传播速度等特征, 并进一步确定疫情峰值及其到达时间、疫情从峰值衰减到正常的时间等重要的疫情传播特征. 但是这些间接数据不能给出疫情起始点信息, 代替不了对疫情监测得到的阳性感染者数据提供的信息.

4.   结论和建议

新冠病毒不断变异和传播, 意味着疫情多次爆发、每次爆发特点各异. 随着防控措施的优化, 原来建立起来的防控体系不再提供详实完整的疫情数据. 从疫情传播动力学机制出发, 采用有限的疫情相关数据, 对每一次疫情进行预测, 成为亟需解决的问题. 在虫口模型基础上, 考虑人们出行和疫情的相互影响, 给出了一个预测新冠病毒阳性感染者人数的离散动力学模型. 假设“布洛芬”搜索量的变化正确地反映了人群中新冠病毒阳性感染者人数的变化趋势, 对上海市2022年12月初爆发的新一波疫情的趋势和关键节点进行了预测. 主要结论如下:

1)本文发展起来的模型能够很好地重现新冠病毒阳性感染者历史数据.

2)对上海市2022年12月开始的这波疫情进行了预测. 如果考虑80%的病毒阳性感染者转阴后不再感染, 这次疫情与2022年3—6月的疫情相比, 阳性感染者人数在几何指数增长期增长速度更快, 在几何指数衰减期衰减速度更慢; 在2023年1月2日到18日处于高发期, 1月10日出现峰值. 如果考虑病毒阳性感染者转阴后不再感染, 阳性感染者人数在2022年12月1日到24日之间几何指数增长; 2022年12月25日到2023年1月8日处于高发期, 2023年1月2日达到峰值; 从2023年1月8日开始几何指数衰减, 2月10日将达到2022年12月初的水平.

3)人们出行过高地考虑每日新增阳性感染者人数的时候(小$\alpha$), 阳性感染者人数达到峰值的时间几乎不受影响, 峰值会明显地变小, 而峰值之后衰退速度明显减缓. 所以, 在抗疫资源和社会对疫情承受力允许的范围内, 应尽可能减小每日新增阳性感染者人数对出行策略的影响.

4) 2022年12月7日疫情防控措施优化后, 人们的出行行为并没有发生明显的改变. 这是政策制定部门需要深入分析的一个现象.

5)需要强调的是, 在疫情初期对模型中参数进行准确的估计, 是模型预测的关键. 上海疫情发生前“布洛芬”搜索量、疫情发生后地铁客流量等疫情相关的数据, 都能提供有效的信息, 但是这些间接信息各有特点, 并不能完全替代直接的病毒阳性感染者数据. 因此, 建议疫情相关部门严格按照统计学中抽样设计方法建立疫情监测体系, 为疫情预测提供直接数据.