流空间邻近关系约束下的流行病分布空间异常探测方法

石 岩，王 达，陈袁芳，陈炳蓉，赵冰冰，邓 敏

中南大学地球科学与信息物理学院，湖南 长沙 410083

流行病长期威胁人类生命健康安全。2019年12月武汉市出现的新型冠状病毒引发的具有高强传播力的肺炎(即COVID-19)疫情[1]，发达的公路、铁路与航空网络为人群跨区域流动提供了极大便利，进一步提高了以武汉为传播中心的疫情扩散与局部区域爆发的风险[2]。探测疫情扩散态势的空间异常区域并探索其时变特征，将有助于疫情时空扩散规律精确捕捉与动态风险有效评估，从而辅助政府部门制定兼具针对性与适应性应急管控决策[3]。

相关领域学者就流行病分布的地理空间异常探测开展了一系列研究工作。①基于数据统计的挖掘。通过对研究区域中病例绝对数量、发病率等疾病专题属性值进行有序统计与分析，发现传染病异常高爆发区域[4-5]。②基于地图可视化的挖掘。由于传统的数据统计方法难以反映传染病的时空分异特征，相关学者受地图思维启发，构建了流行病数据时空可视化分析模型，利用视觉语言形象、直观地揭示流行病时空分异规律。早在1854年，文献[6]将地图可视化方法首次用于疾病归因研究，以点状地图符号的方式标记伦敦市霍乱病例的空间分布，发现了霍乱疫情与公共井水污染之间的直接关联关系。文献[7]采用点状地图符号绘制我国县级尺度乙肝病例空间分布地图，通过视觉分析发现我国中部、南部以及东部的病例分布密度显著偏高。中科院地理所和中国科学预防医学科学院联合编制了《中华人民共和国鼠疫与环境图集》，这是我国首部系统、准确反映鼠疫空间分布的大型专题地图集[8]；针对COVID-19疫情，文献[9]针对省级行政单元，采用基于等级划分的方法绘制包含新增确诊病例数、累计确诊病例数等属性的专题地图，实现了疫情发展的动态可视化展示。③基于空间聚类的挖掘。流行病分布作为一类典型的地理现象，必然具有符合地理学第一定律的空间相关性[10]和第二定律的空间异质性[11]，而以上两类研究由于对数据缺乏时空约束下的统计分析，难以深层挖掘传染病的空间热点、离群等异常分布。为此，地理学领域学者提出基于空间聚类的传染病异常分布挖掘方法，将显著偏离全局/局部分布的空间簇或离群单元识别为异常区域或异常点。例如，空间扫描统计方法借助移动扫描窗口，以空间扩展分析的方式统计判别数据集中具有显著高密度病例分布的热点区域[12-13]；在流行病学中，莫兰指数(Moran’s I)、Getis G等空间相关性指数被广泛用于定量刻画疾病分布的全局与局部相关程度以及热点-冷点区域[14-15]；另外，文献[16]利用蚁群聚类算法探测传染病聚集区域。

流行病传播过程受主观(如防控措施、人群流动等)与客观(如地理环境、人群分布等)因素综合影响[17-18]，不同因素间协同作用的区域差异性导致疫情扩散态势呈现出显著的时空分异特征。然而，现有相关研究大多围绕地理学第一定律[10]，在欧氏空间邻近约束下探测疾病分布的空间异常区域[14-15]，这种策略仅能发现空间邻域范围内的疫情属性特征显著偏离区域，通常具有较强的可解释性(如与传播中心人群交互强度较大导致疫情态势异常严重)，难以深层识别由其他潜在因素引发的空间异常。实际上，区域间人群流动是流行病发展-前期快速传播的核心主导因素[19]。地理学第三定律指出：区域间地理环境配置越相近，则地理专题属性值越相似[20]。为此，本文将由传播中心至各空间单元的人群迁出强度视作一种给定的地理环境配置，那么空间区域之间由传播中心迁入的人群流量序列越相似，其疫情属性特征差异越小。基于以上分析，本文从人群跨区域流动形成的“流空间”视角出发，构建一种定量的流空间邻近关系来度量地理环境配置的相似度，提出一种流空间邻近约束下的流行病空间异常探测方法，以探测人群流动外的其他潜在相关地理环境配置导致的疫情态势空间异常分布区域，从而有效识别不同阶段疫情扩散潜在高风险区域，指导潜在风险诱因精准追踪，辅助区域疫情防控措施的针对性修正与改进。

1 基于流空间邻近约束的流行病空间异常探测方法

给定包含疫情专题属性时间序列的空间单元集合，本文将疫情分布空间异常定义为“任一时段内，疫情专题属性值在流空间邻域范围呈现显著偏离的空间单元”，空间异常在时间维的动态演变可以深层揭示疫情扩散的时空分异规律。针对与疫情发展态势相关的多维专题属性特征(主要包括累计/新增确诊病例数、确诊病例增长率等)，本文首先借助地理探测器[21]对疫情属性特征分布与传播中心人群流出强度进行关联分析，从中提取与人群流动呈显著相关关系的疫情专题属性；其次，基于传播中心人群流出强度的相似性度量计算空间单元间的流空间邻近度，自适应构建流空间权重矩阵；最后，融合疫情专题属性的流空间梯度变化量，提出改进的莫兰指数(Moran’s I)实现疫情专题属性全局空间分布特征判别与局部空间异常区域探测。本文总体研究策略如图1所示。

图1 本文总体研究策略Fig.1 The proposed strategy for detecting epidemic spatial anomalies in this study

1.1 疫情属性特征分布与人群流出强度关联分析

给定研究区域中疫情专题属性空间分布和传播中心人群流出序列， 采用地理探测器定量分析人群流出强度对疫情属性空间分布的解释程度(图2)。

图2 疫情属性特征分布与人群流出强度关联分析Fig.2 Illustration of association analysis between spatial distributions of epidemic attributes and crowd flows

(1) 记人群从传播中心流出至空间单元Vi的强度时间序列为Fi={fi1,fi2,…,fim,…,fit}，其中fim表示第m天传播中心流出至Vi的人群数量与总流出量之比，那么单元Vj(j≠i)与Vi的流空间距离可以表达为

FDij=

(1)

(2) 采用基于模块度的社团发现方法[22]进行空间单元分区。进而，基于q统计量将疫情属性空间分布与人群流出强度间的关联度表达为[21]

(2)

式中，l为人群流出强度分区后的区域数量；n和Nh分别为总单元数和区域h单元数；S和Sh分别表示整体研究区域和区域h疫情属性值方差。最后，采用F分布对q值进行显著性检验[21]，即

(3)

(4)

1.2 流空间权重矩阵自适应构建

空间权重矩阵是规则化定量表达空间单元对间邻近关系的有效工具之一[23-24]。本文将传播中心人群流出强度序列作为疫情态势的地理场景约束，借助高斯核函数将任意两个空间单元Vi与Vj(j≠i)间的流空间邻近度表达为

(5)

式中，带宽σ用于描述流空间中空间单元相互作用范围，对此本文采用交叉验证法自适应确定带宽值[24]

(6)

选择CV最小值对应的带宽为最优带宽。基于此，可以构建以下流空间对称权重矩阵

(7)

为进一步消除量纲影响，对矩阵WF中各元素进行行标准化可得到非对称标准化矩阵WF_S

(8)

1.3 疫情专题属性空间局部变化梯度计算

给定研究区域中第m天疫情专题属性空间分布数据集合Am={a1m,a2m,…,aim,…,anm}，其中aim表示空间单元Vi在第m天的疫情专题属性值，首先将单元Vi的流空间疫情专题属性一阶差异计算为

(9)

在此基础上，采用如式(10)所示的疫情专题属性空间局部变化梯度表征单元Vi的疫情专题属性二阶差异

(10)

式中，疫情属性二阶差异量Gradi越大，单元Vi越可能被判别为异常。例如，图3(b)、图3(c)分别为针对图3(a)专题属性空间分布数据计算得到的空间单元专题属性一阶差异和二阶差异分布，可以明显发现与专题属性一阶差异相比，专题属性二阶差异能够有效消除由于属性值量级不同而引起的空间异常误判，从而真正探测获得具有属性值局部显著突变的空间异常单元。

图3 局部空间异常简例Fig.3 An example of local spatial anomalies

1.4 疫情属性全局空间分布特征判别

从本文对空间异常的定义可以看出，其本质是在研究区域呈现空间正相关性条件下的一种局部空间负相关关系，即空间异常探测的前提在于需要保证数据具有显著全局空间正相关性。对此，本文在流空间邻近关系约束下，采用以下改进的全局Moran’s I指数进行空间相关性统计判别，表达为

(11)

(12)

式中，Iobs为全局Moran’s I观测值；E(Inull)和D(Inull)为零假设下全局Moran’s I均值与标准差，分别表达为

(13)

给定显著性水平α′，若Z>Zα′/2或Z

1.5 局部空间异常区域探测

若疫情专题属性值在流空间具有显著正相关性，融合2.2节构建的标准化流空间权重矩阵与2.3节定义的加权局部变化梯度提出一种改进的局部Moran’s I指数，从而实现在流空间邻近约束下探测疫情属性局部异常区域。具体而言，任一空间单元Vi的流空间局部Moran’s I指数可以表达为

(14)

式中，std(⋅)表示所有空间单元疫情属性局部变化梯度的标准差函数。基于此，在任一空间单元与其流空间邻域单元之间疫情专题属性值相互独立的零假设下，构造局部Moran’s I指数Z统计量，进行显著性检验

(15)

式中，E(Ii_null)和D(Ii_null)分别为零假设下空间单元Vi处局部Moran’s I的期望和方差，计算为

(16)

(17)

若Flow_TADi_m>0，则表明空间单元Vi疫情专题属性值显著异常偏高，反之为显著异常偏低。在此基础上，对于任一空间异常单元Vi，本文进一步定义“显著异常偏高时长”(Ti_H)和“显著异常偏低时长”(Ti_L)用以量化该单元判别为显著异常偏高/偏低的时序变化特征，表达为

(18)

(19)

式中，I(·)为值等于0或1的判别函数；M为研究数据总天数。

2 试验分析

2.1 试验区域与数据

本节采用我国新型冠状病毒肺炎(COVID-19)疫情态势空间分布时序数据，通过与在欧氏空间邻近约束下直接采用疫情专题属性特征的Moran’s I指数分析方法[14-15]和分级可视化方法[9]进行对比分析，以验证本文方法的有效性和优越性。我国各省市于2020年1月24日起开始逐日公开COVID-19疫情发展态势，2月15日后各地新增病例数趋于平稳下降，因此试验数据主要涵盖我国各地市卫生健康委员会在2020年1月23日—2月15日期间每日更新的累计确诊病例数、新增确诊病例数和确诊病例增长率等3类疫情专题属性的空间分布。疫情始发地武汉市自1月23日实施“封城”措施，考虑国家卫健委公布的COVID-19潜伏期最长为14 d的流行病学先验知识[25]，可以推知各地市1月23日的确诊病例最早在1月10日被感染。百度地图迁徙数据(http:∥qianxi.baidu.com/)于2020年1月初，每天对地级市尺度人群跨区域流动情况进行无出行方式差别的更新[26]，并且我国百度用户在各大搜索引擎用户中的渗透率高达90.9%[27]，另外结合现有相关研究可以有力说明百度迁徙数据在反映城市间人群流动情况的实时性和真实性[28]。为此，本文将采用百度地图迁徙数据记录的2010年1月10日—1月23日期间从武汉市流出至各地市的人群比例信息以量化武汉市与各地人群流动强度。

由COVID-19潜伏期的时间滞后影响可以大致推断，2020年1月10日—1月23日期间武汉市的流出人群为其他地市2月4日前新增确诊病例的主要来源。为此，针对新增确诊病例数和确诊病例增长率，仅分析两类属性在2020年1月24日—2月4日期间的空间分布与人群跨区域流动因素间的相关性。基于模块度的社团发现方法将研究区域划分为5个社区，模块度值为0.654，表1为3类疫情专题属性空间分布与人群流动强度间的q值及显著性检验结果，从中可以发现由于确诊病例增长率对累计确诊病例基数变化敏感，使得该变量与人群流动强度相关性较弱；2月4日前新增确诊病例数，以及研究时段的累计确诊病例数均与人群流动强度显著相关。基于以上分析，下面将重点针对累计确诊病例数和新增确诊病例数探测我国地市尺度疫情分布的空间异常区域。

表1 疫情专题属性——武汉人群流出强度间q值及显著性检验结果

2.2 疫情属性分布空间异常探测

首先对武汉迁出至全国各城市人口比例与各地疫情态势属性进行线性回归分析，如图4所示。可以发现整体上武汉迁出人口比例与新增、累计确诊病例数具有显著线性正相关关系，即新增、累计确诊病例数随武汉迁出至各地人口比例单调递增，从而验证了“空间区域之间由传播中心迁入的人群流量序列越相似，其疫情属性特征差异越小”这一论点。

图4 武汉迁出人口比例与疫情态势专题属性回归分析Fig.4 Regression analysis results between outflow population rate from Wuhan and confirmed case numbers in each city

表2和表3分别给出了欧氏空间和流空间两种分析视角下疫情累计和新增确诊病例数属性的全局Moran’s I值与显著性检验结果，两类疫情专题属性在两个空间中的全局Moran’s I值均大于0且在显著性水平α=0.05下显著，表明累计和新增确诊病例数在欧氏空间和流空间中均呈显著空间聚集分布。

表2 欧氏空间和流空间累计确诊病例数空间分布特征判别

表3 欧氏空间和流空间新增确诊病例数空间分布特征判别

在此基础上，考虑到疫情潜伏期影响，1月23日武汉“封城”前流出人口中的潜在病例将在2月4日前被陆续确诊，因此本文将2月4日确定为关键时间节点。进而，分别采用本文方法、欧氏空间局部Moran’s I指数分析法和可视化分析方法对疫情属性分布进行空间异常探测。通过分析探测结果可以发现：

(1) 采用中国疾病预防控制中心的等级划分方法可以绘制图5所示的2月4日累计确诊病例数和新增确诊病例数空间分布专题图。这种人为等级划分的策略仅能凸显如湖北省内各市、北京和重庆等专题数值较高的区域，难以从视觉层面直接从中发现隐含的疫情属性分布局部异常区域。

图5 2020年2月4日疫情专题属性空间分布可视化Fig.5 Spatial distributions of COVID-19 related attribute values on February 4, 2020

(2) 图6和图7表明在欧氏空间邻近约束下，累计和新增确诊病例数一直保持在湖北省内部分城市呈显著HH聚集，并随时间逐渐以武汉市为中心向孝感、黄冈、鄂州、咸宁、黄石、荆门、随州等邻近城市扩张。结合百度迁徙信息可知，上述城市来自武汉的人群流量显著大于未被识别为空间异常的城市。另外，由武汉迁入的人群比例越大，城市判别为累计和新增确诊病例数显著异常的时长值越大。这说明传统欧氏空间局部Moran’s I指数分析法仅能发现与武汉市人群流动强度较大导致的疫情属性空间异常区域。

图6 欧氏空间累计确诊病例数空间异常偏高区域分布Fig.6 Spatial distribution of regions with anomalous large cumulative confirmed case numbers constrained by adjacency in Euclidean space

图7 欧氏空间新增确诊病例数空间异常偏高区域分布Fig.7 Spatial distribution of regions with anomalous large newly confirmed case numbers constrained by adjacency in Euclidean space

(3) 图8(a)和图9所示流空间邻近约束的探测结果表明2月4日前疫情累计和新增确诊病例数的显著异常区域具有相似的空间分布，均主要集中在北京、上海、广州、深圳等超一线特大城市，以及重庆、成都、合肥、杭州等湖北邻近省会城市或直辖市。上述城市不仅吸引各地流动人群聚集，且城市内部人群出行频率较高，导致累计/确诊病例数异常偏高、持续时间长。

图8 流空间累计确诊病例数空间异常偏高区域分布Fig.8 Spatial distribution of regions with anomalous large cumulative confirmed case numbers constrained by adjacency in flow space

图9 流空间新增确诊病例数空间异常偏高区域分布Fig.9 Spatial distribution of regions with anomalous large newly confirmed case numbers constrained by adjacency in flow space

(4) 本文方法还可以准确识别哈尔滨等在疫情期间确诊病例数持续异常偏高的城市。哈尔滨新闻办曾报道当地连续多日出现聚集性疫情病例,截至2020年2月3日63例确诊病例中有58.7%的病例由聚餐等聚集活动被感染(https:∥www.sohu.com/a/370680380_120513233)；在本文研究时间范围外的4月10日—4月12日，该地再次发生由聚餐导致共计10人被陆续确诊的事件(http:∥news.sina.com.cn/c/2020-04-13/doc-iircuyvh7520426.shtml)，说明本文方法能有效提前发现此类空间异常并为相关部门决策提供有效支持。此外，如图10(a)、(b)中长沙市累计和新增确诊病例数在2月4日前在流空间中显著异常偏低，然而随疫情发展该地新增确诊病例数持续上升并呈现异常偏高的态势，说明本文方法探测结果可以为疫情的未来潜在爆发风险预警提供先验信息。

(5) 2月4日后，累计确诊病例数的流空间异常区域在显著减少，再次印证了传播中心人群流出是关键时间点前疫情在我国快速蔓延的关键因素。研究期间内，南阳、温州、台州和宁波等地累计和新增确诊病例数持续异常偏高，其中南阳市是连接武汉市与豫西的重要交通枢纽以及武汉外来务工人员的重要来源地。据统计，武汉市台州籍从商人员达6万余人[29]，经商群体日常接触人员繁杂且春节期间大量台州籍经商人员返乡；截至2月4日，宁波市以及与其流空间邻近的佛山市累计确诊病例数分别为120人和49人，统计两市官方渠道发布的新冠肺炎疫情信息[30-31]可知，两地武汉输入型病例数量接近(分别为25人和29人)，而与输入型病例有流行病学关联的则分别为50人和1人，且百度迁徙信息显示1月10日—2月3日期间宁波市城内人群出行强度持续显著高于佛山市(图11)，为此推测与输入型病例的频繁接触导致南阳、台州、宁波等地疫情初期发展态势异常严峻。值得注意的是，温州市累计确诊病例数始终显著高于武汉人口迁出强度与之十分相近的城市(如西安、南京)，造成这种异常现象的主要原因在于温商占据武汉返温人群主体，且在疫情初始爆发地——华南海鲜市场活动频繁(http:∥www.henan.gov.cn/2020/03-07/1301598.html)。

(6) 通过百度迁徙信息可以发现河南省是武汉人群迁出比例最大地区，但由于采取了新增130所定点医院、全面排查武汉及其他各大城市外来人员等一系列强硬防控措施，使省内扩散得到有效遏制。例如图10(c)中2月4日郑州市累计确诊病例数显著低于武汉人群迁出强度与其相近的深圳、上海等城市。

图10 部分城市疫情流空间异常信息时间序列Fig.10 Time series of epidemic anomaly information of partial cities obtained in flow space

图11 宁波市和佛山市城市内部人群出行强度Fig.11 The intensity of intra-city human mobility in Ningbo and Foshan

2.3 讨 论

通过以上试验分析可以发现，传统基于欧氏空间邻近约束的异常探测策略仅能发现由传播中心人群大量流出导致的空间异常区域，这种具有强解释性的浅层异常难以提供深层防控决策信息，且直接采用疫情专题属性作为空间单元特征值仅能发现与专题属性平均值差异较大的全局异常区域。与此相比，本文方法在流空间中构建区域间邻近关系，并构造疫情属性值空间局部变化梯度变量作为异常度指标，可以准确发现由于区位条件、职业组成、居民疫情防控意识等传播中心人群流动之外的因素引发的空间异常，真正为疫情分区分级的精准防控提供风险评估决策支持。

然而，本文方法仍然存在以下问题和局限性：①该方法主要考虑传播中心人群流出强度对疫情跨区域传播的影响，难以适应多中心传播或非传播中心区域间人群流动驱动下的疫情扩散情景；②单来源信息(如本文采用的百度迁徙数据)由于不可避免地存在采样偏差，导致基于此度量的人群跨区域流动强度有偏，另外，百度迁徙数据不区分交通出行方式，从而无法有效反映不同出行场景中疫情传播环境差异性对疫情扩散过程造成的差异化影响；③除传播中心人群流出强度因素外，缺乏考虑区域内其他疫情相关因素(如区域内部人群活动强度、确诊病例空间动态分布等)，难以对空间异常进行更深层的挖掘与解释；④虽然本文方法探测结果能够从空间域视角支持疫情的早期预警，但仍然难以在时间域确定疫情在潜在高风险区域的具体爆发时间或时段。

3 结 论

本文关注流行病在空间区域的发展态势，针对现有欧氏空间邻近约束的异常探测方法难以有效识别人群跨区域流动之外的其他因素导致的潜在空间异常这一问题，从流空间邻近关系约束的视角，采用一种流空间距离加权的疫情专题属性二阶差异作为空间单元疫情特征变量，提出改进的Moran’s I指数模型，用于深层探测疫情专题属性在流空间具有显著局部偏离的空间异常区域。试验结果表明，本文方法可以准确探测疫情不同发展阶段由各类因素导致的空间异常以及蕴含潜在疫情暴发风险的异常区域，能够科学指导各地政府相关部门在疫情不同发展阶段展开针对性的分区分级防控，同时能够在城市尺度上指导疫情潜在爆发风险预警。

下一步研究工作将进一步融合不同平台(如铁路、长途客车、航空等)记录的任意区域间人群迁徙信息，构建多元出行场景下人群流动蕴含的区域流空间邻近关系，全面获取更精细的疫情分布空间异常区域探测结果。另外，需要关联社交媒体签到、用户移动轨迹等个体粒度时空大数据，从舆论扩散的空间分异性、城市内部确诊病例活动规律与人群聚集性等多个视角进一步剖析疫情扩散的时空分异特征，并将此作为一种先验信息嵌入复杂网络或深度学习模型实现疫情演化过程的动态预测。