新冠肺炎疫情城市个体仿真模型及评估系统

开云电竞pg软件测绘学院教授

张健钦

1 模型意义

2019年12月以来，湖北省武汉市部分医院陆续发现了多例有华南海鲜市场暴露史的不明原因肺炎病例，后证实为2019新型冠状病毒（SARS-COV-2）感染引起的急性呼吸道传染病COVID-19。随着肺炎疫情的蔓延，有必要对疫情的趋势分析，后续各种复工等措施造成的风险进行模拟和评估。同时，模拟病毒在城市中的传播，也能够为以后类似的病毒传染病的防治提供一定的参考。

2 传染病模型研究概述

目前，对于传染病模型的研究主要有以下几类：

第一类是动力学模型。此类模型主要通过各类微分方程进行构建。常见的模型有SIR模型、SEIR模型等。许多学者在这类模型的建立与应用领域进行了详细的研究。例如陈田木(2015)建立有、无干预措施的易感者-潜伏期者-染病者/隐性感染者-移出者-传播介质模型, 利用长沙市戊肝发病数据、人口学数据、干预措施的特点和模型假设，通过Matlab7进行仿真计算求解；陈慧林(2017)等基于已有的SEIR模型,加入隐性感染者分箱得到SLICAR模型,并以2009年春季的甲型H1N1和1918年秋季的Spanish流感数据为实例分析；钟剑明(2019)等过SEIAR模型对流感暴发疫情隔离和停课效果进行分析。这类模型适合群体模型的数值模拟，能够很好地反映的疫情的发展趋势。

第二类主要是统计学模型或机器学习模型。此类模型主要借助各类机器学习的算法进行研究与预测。常见的模型包括回归模型、ARIMA、LSTM等。这些模型在传染病疫情预测上有着广泛的应用。例如刘如春(2014)使用Logistic模型对长沙市某校甲型H1N1流感暴发疫情进行拟合,确定模型的参数,并进行拟合优度检验,探讨模型；曾四清(2019)利用Joinpoint回归模型对2014年广东省登革热大流行疫情为例进行分析,探讨该方法用于传染病大流行疫情时间趋势变化特征分析及预测方面的适用性；黄国(2019)构建时间序列ARIMA乘积季节性模型,预测江门市手足口病发病趋势,探讨该模型在预测手足口病发病率中的应用。从这些例子可以看出，这一类的模型的优点在于能够针对已有的疫情真实数据进行拟合，进而对接下来疫情演变进行预测。

第三类是智能体仿真模型以及网络模型。这一类模型主要通过agent仿真或者网络和进行模拟分析，能够反映疫情的扩散细节，结合可视化技术可以对政策评估和后期疫情演变预测提供一定的参考。因此许多的学者也通过此类模型对传染病进行模拟与预测，如赵翌僮（2017）采用基于Agent的仿真建模方法,生成北京暴发H1N1流感的仿真数据，基于4组分布形状与规模不同的数据,对2种空间聚类算法的疫情监测结果进行分析；瞿毅臻（2008）用Repast仿真建模平台对SARS的传播过程进行模拟,从而研究SARS传播宏观层面的规律,并根据实验结果提出了疫情防控的策略和措施；黄甜（2015）基于复杂网络理论构建了手足口病疫情预测预警模型，并探讨了在云南省手足口病疫情防控中的应用。

综上所述，动力学、统计学或机器学习模型可以在数值上做出合理的预测，适合建立群体模型，但不能反映疫情演变的空间特征与相关的细节。智能体仿真模型或网络模型可以反映模型计算过程的细节，能够反映患者的空间位置变化，适合建立个体模型，但是智能体模拟过程中需要进行大量的并行计算，因此运算量有上限限制。在城市疫情模型研究中，群体模型和个体模型都能够反映疫情的演变特征和政策调控的效果，但个体模型的结果更适合用于可视化与反映疫情扩散细节，决策者与公众可以通过空间数据的可视化，直观了解疫情变化和评估风险。因此采用个体模型进行疫情扩散的模拟与复工复产政策的评估。

3 研究区域与数据

3.1 研究区域概况

本次研究选取的地区是北京市，北京市位于华北平原北部，地势西北高东南低，气候为温带大陆性季风气候，是全国的政治中心、文化中心、国际交往中心、科技创新中心，常驻人口2153.6万人。北京是中国铁路的枢纽之一，拥有2座大型的民用机场。市内交通发达，人员出行量大。

历史上北京曾经出现过多次传染病流行的情况，最早记载是在元朝至元年间，之后又多次传染病暴发的情况，主要出现过的疫情有天花、伤寒、鼠疫、白喉等。最近的疫情大规模流行是2003年的“非典”和2010年的甲型H1N1流感。本次新冠肺炎疫情从1月20日出现2例输入性病例开始影响北京。1月24日，北京启动一级应急响应，停止各类大型活动。截止3月11日，北京出现病例435例。疫情也主要从开始的湖北地区输入、本地扩散逐渐转为境外输入。

3.2 研究数据及处理

3.2.1 基础地理信息数据处理

基础地理数据主要用到了北京市建筑物矢量数据与土地利用数据、北京市街道人口数据。其中建筑物矢量数据包括建筑物层数、面积等信息。土地利用数据包括土地代码类型等信息。具体数据如下图所示。

图3-1 街道人口数据示例

图3-2 北京建筑轮廓数据

图3-3 北京市土地利用数据

数据处理大致思路为：利用ArcGIS空间连接工具赋予建筑物的类型，统计街道内的住宅区面积总和，进而计算建筑物内的人口密度，从而计算每个建筑物内的人口数量。流程图如下所示。

图3-4 基础地理信息数据处理流程

3.2.2 城市居民出行统计数据

城市交通出行特征统计数据来源于《2019年北京市交通发展年度报告》，本文主要参考了其中的出行目的与出行时间分布特征，用于确定智能体的出行规则。

出行目的：劳动人口（上下班）约占45%，上下学约占5.7%。

图3-5 北京市居民出行目的分布（来源于北京市交通发展年报）

出行时间分布：早高峰时段7:00-9:00出行量占全日出行总量的32.7%，晚高峰时段17:00-19:00出行量占全日出行总量的20.8%。

图3-6 北京市居民出行目的分布（来源于北京市交通发展年报）

此外我们引入了百度地图的迁徙大数据，用于分析城市居民的迁入迁出状况，为后面复工复产分析提供数据 支撑。

图3-7 迁入迁出数据

3.2.3 疫情数据

疫情数据截取自丁香园的疫情数据，主要包括日期、确诊、治愈、死亡、康复、疑似病例等数据。疫情数据如图所示。

图3-8 北京市疫情数据

4 模型设计

4.1 人群智能体规则

人群智能体的属性主要包括个人识别ID、出行群体、地址信息、状态。其中个人识别ID（pid）为智能体的唯一识别特征。出行群体记录了智能体所属出行群体类型（学生、上班族、其他）。地址类信息主要有区县、住址、时刻和位置。住址用生成该智能体的所在建筑物ID表示（home_id），时刻和位置则记录了智能体不同时刻和在该时刻所处建筑物的ID（pos）。状态（state）则记录了智能体不同健康状态，包括易感（S）、无症状感染或潜伏（E）、轻症（Im）、无症状（In）、进入医院隔离救治（H）、确诊（C）、重症（Is）、恢复（R）、死亡（D）。具体属性字段如表3-1所示。

表4-1 人群智能体属性列表

属性字段	数据类型	说明
pid	int32	个人识别id
home_id	int64	住址id
county	object	区县
group	object	出行群体
state	object	状态
xy	object	时刻和对应建筑物位置id
pos（i）	int64	i时刻所在建筑物id
infected_duration	int64	潜伏期
seviere_duration	int64	轻症转重症时长
hospital_duration	int64	治疗时长
hospitalization	bool	是否收治
I_duration	int64	病程
confirm	bool	是否确诊

人群智能体的规则主要有两部分内容：感染规则和出行规则。

首先是出行规则，出行规则根据出行群体的不同分为三组。根据不同组别将为每个智能体分配不同时刻所在对应建筑物的ID。其中，上学和上班族各个时长服从正态分布，其他人群出行特点服从均匀分布。

上学：出行的主要轨迹为家-学校-（其它）-家。出发时刻为早上7点前后，往来通勤时长约15到20分钟，在校时长约8到10个小时。

上班族：出行轨迹为家-工作场所（办公用地、机关政府、工矿用地）-（其它）-家。出发时刻为早上7点-8点前后，往来通勤时长约为45分钟到1小时，工作时长约为8到10个小时，休闲或其他活动时间为2小时左右。

其他人群：出行轨迹为家-其它-家。出行时刻在早7点至晚7点。外出时间最多约6小时。

其次是感染规则。感染规则为：根据不同时刻和不同建筑物内的人员分布情况，采用小世界网络的方式，表示感染者所接触其他智能体。并在相应的小世界网络上完成各个状态的转换，状态转换的规则如下图所示。

图4-1 智能体转换规则

状态转换方法为，出现感染者时，感染者以及与感染者接触的智能体构成一个小世界网络。每个智能体会生成一个随机数，当该随机数小于某一状态的转换概率时，智能体向该状态进行转换。

4.2 环境智能体规则

环境智能体为建筑物、政府（社会）政策，其中建筑物智能体的属性包括建筑物所属ID，经纬度、楼层信息、建筑类型等。政府（社会）的属性为时间和对应政策（触发事件）。建筑物属性如表3-2所示。

表4-2 建筑物智能体属性

属性字段	数据类型	说明
gid	int64	建筑物id
floor	int64	楼层数
level	int64	建筑物类型代码
x	float64	经度
y	float64	纬度

环境智能体的规则有两方面：生成小世界网络和按时间调整政策。

在感染者所在的不同的建筑物类型内，根据不同的建筑物类型，生成对应场景的小世界网络。网络的节点数、网络数、重连概率根据不同场景的中小尺度模型测定。

另一个规则是根据不同时间段生成相应的政策，疫情的传播时间节点主要有3个：

（1）输入性病例出现：1月20日出现2例输入性病例。

（2）1月24日，春节假期开始，同时北京市启动一级应急响应，各个小区开始封闭，人员出入停滞。

（3）2月9日开始，北京各个企业开始复工，上班族恢复出行轨迹，学生和其他人群依旧居家。

详细规则设置如表4-3所示。

表4-3 不同时间阶段与对应复工原则

时间阶段	规则描述	参数滤定
1月18日-1月23日	1月21日广州市第一例输入病例 1月23日广州市启动重大突发事件一级响应，人的行为规则不变	考虑到病毒潜伏期的原因将模型运行的时间提前3天 出行规则参考武汉市个体模型
1月24日-2月10日	2月10日广州开始复工， 只有部分的人出行，且每天仅出行一次，目的地随机，一个小时之后回到家中	出行人数人为设定为出行率αMove ="0.001 设定仅潜伏期E和易感者S可以出门，其余状态如果不被医院收治则呆在家中。
2月10日-3月18日	开始复工，学生呆在家中，退休人员和上班族维持原来出行规则。总计模型模拟50天	考虑复工的实际情况设定规则 可设定复工率αWork 设定仅潜伏期E和易感者S可以出门，其余状态如果不被医院收治则呆在家中 35day之后强行收治所有非易感者，不具备传染性

5 模型实现

5.1 模型初始化

剔除完建筑物中异常值数据后，根据处理完的住宅数据和城市的出行结构模拟生成了不同的人群仿真数据，类别分成上班族、学生、其他人群，每一个人拥有ID、区县、出行时间地点、状态、群体、住址ID、接触过其他人等属性。流程图如下所示。

图5-1 模型初始化流程

生成相应的人群数据后，将人群数据中的时刻与位置数据填充至24小时，从而得到了一个人24小时中的每一个时刻的位置信息，从而完成人群智能体数据构建。填充方法为：获取人群智能体的指定时刻与所在位置，在每一个时刻区间的中间时刻中，填入区间左端点时刻的所在建筑物id。同时，从建筑物的矢量数据中获取建筑物的id、楼层数、建筑类型、经纬度，生成建筑物智能体。这两组数据用于接下来的网络构建与仿真计算。

初始化完成后，生成的人群智能体的出门时间分布如下图所示。

图5-2 人群智能体出门时间分布

5.2 模型运行结果分析

模型利用python实现，并利用pymc3库测定了模型中的主要参数（如感染率）。模型初步模拟了50天的疫情结果，结果如下图所示。

图5-3 模型运行结果

从模拟结果来看，模型较为贴近实况，增长较快的时间段是第10~20天，之后疫情增速放缓。