中国境内COVID-19的空间格局演变及聚集性分析

陈进钊，梅志雄，吕佳慧

(华南师范大学 地理科学学院，广东 广州 510631)

0 引言

新型冠状病毒肺炎(Corona Virus Disease 2019，COVID-19)是一种潜伏期长、具极强传染性的疾病。2020年1月，随着中国春运人口大规模流动，COVID-19在短时间内从湖北省武汉市蔓延至中国多地，引起中国乃至全球的高度关注。相关领域学者在COVID-19基因组表征、溯源分析、流行特征、风险评估、检测手段、病患治疗、临床表现、药物研发、致病机理、传播途径等方面进行了研究[1]，为了解、检测、防控COVID-19及对病患的救治提供了较好指导。

传染病的聚集性特征是疾病空间流行规律研究的重要组成部分[2]，对其进行研究有助于更好地把握传染病发展及其科学防控和评估已采取的防控措施等具有重要意义[3]。一些学者研究了SARS[4]、手足口病[5]、甲型H1N1[6]、流感[7]等传染病的时空分布及聚集特征，挖掘出各传染病的高发聚集时间和聚集区域。这些成果为本文提供了有益借鉴。目前，也有少数学者从时空分布视角对COVID-19进行了研究，如曹培明等[8]、王宣焯等[9]、王刚等[10]、刘逸等[11]、刘郑倩等[12]利用GIS制图和一般统计方法，描述性分析了COVID-19在省、市内部的时空特征及其变化趋势。这些研究为相应区域COVID-19防控提供了一定的参考，但缺乏从空间聚集演化和时空聚集性方面进行探讨，未能揭示COVID-19的空间集聚特征和时空聚集热点区。苏理云等[13]虽对中国各省COVID-19累计确诊人数进行了空间自相关分析，探讨了其空间集聚特征，但未进行时空聚集性热点探测，也未揭示更小尺度如地市尺度上疫情聚集性特征。刘勇等[14]、简子菡等[15]分别从区县尺度和镇域尺度对河南省的时空扩散特征、空间格局演化特征进行分析，但其研究时间段皆较短，且未能进行全国范围的研究。还有学者从医药学视角发表了COVID-19代表性研究论文[16]，但缺乏地理学视角下的分析。

从地理学视角，挖掘中国境内COVID-19的空间聚集演化和时空聚集规律对有效实施联防联控有重要意义。为此，本文对2020年1月11日至7月31日中国地级市的COVID-19疫情数据进行空间自相关分析和时空聚集性探测，挖掘COVID-19疫情在中国各地市的空间聚集演化规律及时空爆发热点情况，为全面了解该疫情时空分布特点、制定行之有效的防控政策提供一定的科学参考。

1 数据与研究方法

1.1 数据来源及处理

本文以中国地级市为基本空间单元(台湾省、香港和澳门的数据不全而不包含在研究范围内)，把北京、上海、天津、重庆几个直辖市分别作为一个单元，最后形成362个研究单元。COVID-19数据来源于国家、各省和自治区的卫生与健康委员会官网，由于各级卫健委从2019年12月31日开始公布疫情数据，但2020年1月1日～1月10日未公布数据，因此本次研究时段选择从2020年1月11日到7月31日。所需人口数据来源于中华人民共和国第六次人口普查数据，行政区划图取自中科院资源环境数据云平台(http：//www.resdc.cn/)。

为了方便表达时间跨度达半年的疫情空间聚集演化和时空聚集特征的变化，从1月11日开始，以7天为时间间隔选取时间段，选取至2月29日；从3月份开始以月份为单位选取时间段，共选取得到12个时段。本文分别对中国地市每一时段内COVID-19新增患病率进行空间自相关分析以探测其空间集聚规律。通过时空扫描统计量方法进行时空集聚性分析时，则使用每日COVID-19新增病例数组成的序列，由于时空扫描统计量方法不允许单日确诊人数小于0，故将研究数据中由于核减产生的新增确诊病例小于0的疫情数据设置为0，以便适用于时空扫描统计量方法分析。

1.2 研究方法

1.2.1 空间自相关方法

空间自相关度量指标有多种，本文采用常用的Moran’sI指数，通过GeoDa软件依据Queen法(适合面状对象)构建空间权重矩阵，进而计算全局和局部Moran’sI指数。

全局Moran’sI反映区域单元观测值在整个研究区的空间自相关特征，其计算公式为[17]：

(1)

全局Moran’sI虽可度量空间中相似属性的聚集程度，但不能指出聚集区的具体位置。局部空间自相关可揭示空间异质，识别聚集区的空间位置和范围[18]。局部Moran’sI计算公式为[19]：

(2)

式中：各指标的含义同式(1)。

局部空间自相关计算结果可划分为高—高、低—低、高—低、低—高4种类型。高—高(低—低)型表示某单元与其相邻单元间存在正空间自相关，为高(低)值在局部空间上的聚集模式，Ii为正；高—低(低—高)型表示该单元与其相邻单元存在负空间自相关，高(低)值单元被低(高)值单元包围，Ii为负[20]。

1.2.2 时空扫描统计量方法

时空聚集性分析可从时间和空间两个维度对疫情的聚集进行精确定位，并定量计算出疫情的相对危险程度。本文通过SaTScan软件的回顾性时空扫描方法进行时空聚集性探测，可挖掘COVID-19疫情的时空聚集变化趋势及其高发区域和时段。

时空扫描统计量方法是基于一个移动的圆柱形活动窗口，通过改变窗口在空间上的圆心和半径以及在时间上的扫描时段，根据扫描圆柱体内部事件数和根据某一理论分布而得到假设事件数构造对数似然比(Log Likelihood Ratio，LLR)。LLR可用于评价扫描圆柱体内部时空的发生异常情况[21]。基于离散型Poisson模型假设，其LLR为[22]：

(3)

式中：c为窗口内发病数；C为研究范围总发病数；E[c]为基于无效假设由协变量校正过的窗口内预期发病数。I为扫描指示值，进行高发病率聚集扫描时，若窗口内实际发病数高于预期发病数，则I=1，反之I=0。

选取LLR值最大的窗口为一级高发病聚集窗口，确定该窗口所包括的单元为一类时空聚集区。其他具有统计意义的窗口，根据LLR大小对它们进行排序，确定这些窗口对应的单元为二类时空聚集区。最后计算探测到的时空聚集区的相对危险程度(Risk Ratio，RR)[23]。

(4)

式中：各指标的含义同式(3)。通过RR值可评估各聚集区的相对风险大小。

本文采用基于Poisson概率模型的时空扫描统计量的方法，时间单位为1天，时间扫描下限为3天、上限为研究时长的50%，空间扫描上限为10%总人口所对应的区域。时间趋势调整方法采用模型自动计算的对数线性趋势法，该方法可根据数据中观察到在时间上发病数的整体变化趋势来对假设理论分布进行调整。

2 结果分析

2.1 COVID-19的空间聚集格局演变

2.1.1 总体空间聚集格局的演变

利用式(1)计算各时段中国地市COVID-19新增患病率的全局Moran’sI值和相关指标(表1)。由表1知，全局Moran’sI值在1月11日～1月18日时段不显著且值较小，但在1月19日开始至2月29日的各时段均显著为正，表明1月19日～2月29日各时段中国地市COVID-19患病率总体上具正空间集聚特征。但随时间的推移，其聚集程度呈先增强后减弱趋势。1月26日～2月1日时段全局 Moran’sI值达到最大值，说明该时段中国地市COVID-19新增患病率的总体空间差异最小，COVID-19疫情从武汉市向中国其余地区快速蔓延使得多地市COVID-19患病率都经历了较快增长。2月2日后全局 Moran’sI值逐渐下降，说明中国地市COVID-19每一时段内新增患病率的总体空间差异又逐步扩大，区域极化效应逐渐显现，这表明2月2日后中国统一实施严格防控措施，如限制交通和人员流动与聚集等逐渐显现出效果，只有少数疫情核心地市如湖北省部分地市仍然形式较严峻外，中国其它地市疫情总体上得到逐步控制。从3月份开始，全局Moran’sI值不显著且值较小。这表明疫情得到有效的控制，因此从地理视角看，并未呈现显著的聚集分布。

表1 中国地市COVID-19患病率的全局Moran’s ITab.1 Global Moran’s I of the prevalence of COVID-19 in China

2.1.2 局部空间聚集格局的演变

上文显示中国地市COVID-19各时段内患病率的空间聚集性显著时段为1月19日～2月29日的6个时间段，故进一步对该6个时间段的局部聚集特征进行分析，利用式(2)计算各地市COVID-19新增病例患病率的局域Moran’sI值，并通过GIS软件得到局部空间聚集演变图(图1)。

从局部空间聚集时空分布看，高—高聚集区集中分布在武汉市及其周围地市，其原因是中国境内COVID-19疫情始于武汉市，在1月23日武汉封城前武汉周围地市与武汉市人员流动规模大而受到武汉市疫情扩散影响较大，造成武汉周围一些地市较多患者被确诊从而形成以武汉为中心的高—高聚集区。高—高聚集区范围于1月26日～2月1日时段达到最大、2月2日后逐步收缩。低—低聚集区1月19日～2月8日主要分布在中国西部和东北、西南、华北的部分地市。 2月9日～2月22日低—低聚集区散布在中国西部以及东北和华南等部分地市，对应的地市数量减少。2月23日～2月29日，中国未出现明显的低—低聚集区。

低—高区在1月19日～1月25日时段主要分布在湖北省内及与湖北省相邻的部分地市，且与高—高区相互交错分布，在1月26日～2月1日时段与高—高区交错分布态势消失且有整体外移趋向，这表明此时段疫情在原交错地带迅速蔓延，原低—高区向高—高区转变，而非交错地带受影响相对较小而依旧为低—高区；低—高区范围2月2日后范

围呈现收缩态势，随后趋于稳定。高—低区1月19日～1月25日零星分布在中国多地，如四川省甘孜藏族自治州、云南省丽江市、甘肃省兰州市、陕西省铜川市、山东省日照市和内蒙古锡林郭勒盟，1月26日后，高—低区仅在个别时段零星出现，如2月2日～2月8日的贵州省毕节市，2月9日～2月29日的四川省甘孜藏族自治州和2月16日～2月22日的山东省济宁市。高—低区和低—高区表示高低(低高)患病率的不同值在空间上聚集的效应，是空间异质区。高—低型地市具强烈的空间极化特征，低—高型地市受到周边高值地市的极化影响，如低—高区基本都出现在高—高聚集区周围，表明这些地市在外防疫情高发地输入病例的工作效果显著。另外，研究期内高—低区较长时间出现在四川甘孜藏族自治州，其累计患病率达7.14人/10万人，明显高于其周边相邻地市，该州的道孚县占该州93.6%的确诊病例，据相关报道这与疫情初期该县未能有效落实各项防控措施有关。位于高原上而非位于东部发达地区或湖北省周边地区的甘孜藏族自治州出现这种情况，也警示面对如此高传染性的疾病，任何地区防控措施不许有半点懈怠，即便偏远的地区也会存在疫情高发的风险。

2.2 COVID-19时空聚集性分析

采用SaTScan软件对各地市每日COVID-19新增病例数组成的序列进行时空聚集性分析，探测到具有统计学意义的时空聚集区共18个，LLR最大的时空聚集区为一类时空聚集区(1个，用A标识)，其余为二类时空聚集区(17个)。按聚类时间的起始时间对二类聚集区进行先后排序，使用B1，B2，…，B17进行标识。3月16日国务院联防联控机制新闻发布会上指出 “外防输入” 已经成为疫情防控的重中之重。为了更好展示境外输入时期COVID-19的时空聚集情况。以3月16日为时间点将研究时段分割为两个时间段，分别作图展示(图2)。

2.2.1 COVID-19时空聚集区描述

如图2(a)和表2所示。一类时空聚集区(A区)包括湖北省的武汉、鄂州和孝感3市，其聚类时间段为1月27日至3月1日，根据式(3)和式(4)得到其相对危险程度(RR)为376.78、对数似然比(LLR)高达234 872.51。3月16日前的二类时空聚集区中，除B2、B3和B6聚集区离COVID-19爆发地武汉市相对较近外，监测到的其他二类时空聚集区分散在其他省域的少数地市、为本次疫情在中国范围的散布型高发时空热点区。其中，B2区聚类时段为1月28至2月13日，位于湖北省东部、安徽大部分地市以及浙江西部、福建西北部、江苏西部、江西北部的部分地市，其RR为2.88、LLR为1 886.00；B3区聚类时段为1月29日至2月14日，位于湖北中西部、重庆市以及陕西省东南部、贵州省东北部、湖南省西北部和四川省东部的部分地市，其RR为4.39、LLR为4 491.00；B6区聚类时段为2月1日至2月6日，位于河南省大部分地市和山西省东南部、安徽省西北部、湖北省北部的部分地市，其RR为1.83、LLR为182.30。其余的散布型高发时空热点区中，B1区聚类时段为1月27日至2月7日，位于浙江省温州市，其RR为3.72；B4区聚类时段为1月29日至2月5日，位于广东省深圳市，其RR为2.57；B5区聚类时段为2月1日至2月12日，位于安徽省蚌埠市，其RR为3.99；B7区聚类时段为2月4日至2月12日，位于海南省三亚市和保亭黎族苗族自治县，其RR为5.55；B8区聚类时段为2月9日至2月19日，位于四川省甘孜藏族自治州，其RR为6.10；B9区聚类时段为2月20日至2月22日，位于山东省济宁市，其RR为15.48。

3月16日后的二类时空聚集区在空间上呈现散布式分布，聚类时段总体上相比于3月16日前的二类时空聚集区跨度长，呈现复杂多变的特点。如图2(b)显示，B10区聚类时间段为3月21日至4月12日，位于上海市，其RR为2.70；B11区聚类时间段为5月31日至7月26日，位于成都市，其RR为3.82；B12区聚类时间段为6月5日至7月31日，位于广州市，其RR为10.11；B13区聚类时间段为6月12日至6月29日，位于北京市，其RR为64.25；B14区聚类时间段为6月16日至7月4日，位于兰州市，其RR为10.11；B15区聚类时间段为7月5日至7月26日，位于内蒙古中部，陕西省北部，宁夏回族自治区北部地区，其RR为8.56；B16区聚类时间段为7月17日至7月31日，位于乌鲁木齐市，其RR为2 337.98；B17区聚类时间段为7月23日至7月31日，位于大连市，其RR为265.87。

2.2.2 COVID-19时空聚集性特征及原因分析

对于3月16日前各时空聚类区，根据各RR值大小可知：一类时空聚集区的相对危险程度远大于其余时空聚集区，疫情极为严峻。在3个邻近一类时空聚集区的二类时空聚集区(B2、B3、B6)中，位于武汉以西的B3区的疫情最为严峻，而位于武汉以北的B6区的疫情严重程度最低，位于武汉以东的B2区介于两者之间；其余散布型高发时空热点区中，B9区(山东省济宁市)发生监狱聚集性疫情爆发事件，疫情严重程度明显比其余二类聚集区高，B1区(温州市)、B5区(蚌埠市)、B7区(三亚市和保亭黎族苗族自治县)和B8区(甘孜藏族自治州)的疫情严重程度相对也较高。而B4区(深圳市)的疫情严重程度在散布型高发时空热点区中相对较低。

结合国家及各级卫健委公布的每一病例的具体来源，B10(上海市)，B11(成都市)，B12(广州市)，B14(兰州市)和B15时空聚集区的病例来源皆为境外输入，为境外输入型时空聚集区。境外输入型时空聚集区的病例来源大多比较固定且单一，B10聚集区中来自于欧美国家的输入型病例占91.98%，B11聚集区中来自埃及的输入型病例占84.62%，B14聚集区中来自于沙特阿拉伯的输入型病例占85.71%，B12聚集区中来自于南亚、东南亚的病例达75%。这与国际航班的开通情况有关，且时空聚集区RR值和聚类时段与航班来源国的疫情有较强的相关性。B15聚集区的形成原因主要是呼和浩特白塔国际机场在北京市爆发本土疫情后，接受了原本飞往北京的航班。B13(北京市)，B16(乌鲁木齐市)和B17(大连市)皆为本土病例，为扩散型时空聚集区。扩散型时空聚集区的RR值均高于境外输入型时空聚集区。其中最严重的B16(乌鲁木齐市)的疫情爆发和一起聚集性活动关联，且疫情发展较快；B13(北京市)和B17(大连市)的疫情爆发于当地的菜市场，与受污染的进口海鲜有关。总体而言，3月份后中国进入严防境外输入时期，从地理视角看在中国范围内疫情趋于分散性，局部性且形势多变。相比于境外人员，对境外进口的海鲜进行全面排查难度较大，容易出现关联的本土病例，并存在短时期内出现爆发性疫情的可能性。

表2 全国地市COVID-19疫情时空聚集性分析结果Tab.2 Spatio-temporal cluster results of COVID-19 epidemic at prefecture level in China

3 结论与讨论

3.1 结论

(1)1月19日开始至2月29日的各时段中国地市COVID-19患病率的全局空间集聚在均显著为正，中国地市COVID-19患病率总体呈正空间聚集特征，且聚集程度呈先升后降趋势。其余时间段中国地市COVID-19患病率的全局空间集聚性不显著。

(2)中国地市COVID-19患病率的局部高—高聚集区主要分布在武汉市及其周围地市，1月19日～1月25日呈以武汉为中心向外扩张，高—高聚集区范围于1月26日～2月1日时段达到最大、2月2日后逐步收缩。低—低聚集区1月19日～2月8日主要分布在中国西部和东北、西南、华北的部分地市。2月9日～2月22日低—低聚集区散布在中国西部以及东北和华南等部分地市，对应的地市数量减少。2月23日～2月29日中国未出现明显的低—低聚集区。低—高区主要分布于高—高区周围，1月19日～1月25日与高—高区相互交错分布，1月26日～2月8日该交错分布态势消失并呈整体外移趋势，随后保持稳定。高—低区在1月19日～1月25日零星分布在中国多地，随后仅在个别时段零星分布。

(3)中国地市COVID-19疫情一类时空聚集区位于湖北省武汉市、鄂州市和孝感市，聚类时间段为1月27日～3月1日。该时空聚集区的相对危险程度远大于其余时空聚集区，且对应的聚类时间较长，为本次疫情最严峻的区域。二类时空聚集区中在3月16日前的共有9个，在3月16日后的共有8个。3月16日前的时空聚集区中有3个(B2、B3、B6)位于一类时空聚集区周围，其中B2和B3的相对危险程度较高，且持续时间长，为疫情初期防控救助的重点区域。其余6个二类时空聚集区聚类时间较短且包含的地市较少，皆为散布型高发时空热点区。3月16日后的8个时空聚集区中，3个(B13，B16和B17)为扩散型时空聚集区，其余皆为境外输入型时空聚集区。扩散型时空聚集区的相对危险程度更高，且与本地聚集活动或境外受污染海鲜的输入有关；而境外输入型时空聚集区的相对危险程度较低，且与境外病例来源国的疫情严重程度有较强的相关性。输入病例来源较为单一，与航班开通情况相关，时空上整体呈现多变态势。

3.2 讨论

本文研究结果可为了解本次COVID-19疫情中国范围内的区域聚集性发展变化，从而制定行之有效的联防联控政策提供一定的参考。但是，本文研究数据仅采集中国境内的病例数据，未能深刻揭示中国疫情与世界疫情发展的关系。在外防输入时期，分析中国和世界的疫情时空关联性，寻找其中的时空关联性对外防输入工作有重要意义。因此，后续需要用关联思维从地理视角审视疫情能得到更深刻的研究结果。另外，时空扫描统计量方法应用于COVID-19分析时，采用不同参数产生时空聚集区的空间范围和聚类时长也有一定差异，采用何种参数最有利于揭示疫情的时空热点也是需要进一步研究的问题。本文采用一种基于模型自动计算的对数线性趋势法对病例数在时间维度上对理论假设病例数进行调整，可将时空扫描统计量方法拓展至较长时间段的，病例总数变化范围较大的疫情聚集性探测，且得到良好的效果。但随着时间的推移，RR值呈现过度偏高的现象。从方法上，如何对理论假设病例数进行调整也是一个需要进一步研究的问题。