流空間鄰近關系約束下的流行病分布空間異常探測方法

石 巖，王 達，陳袁芳，陳炳蓉，趙冰冰，鄧 敏

中南大學地球科學與信息物理學院，湖南 長沙 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指數模型，用于深層探測疫情專題屬性在流空間具有顯著局部偏離的空間異常區域。試驗結果表明，本文方法可以準確探測疫情不同發展階段由各類因素導致的空間異常以及蘊含潛在疫情暴發風險的異常區域，能夠科學指導各地政府相關部門在疫情不同發展階段展開針對性的分區分級防控，同時能夠在城市尺度上指導疫情潛在爆發風險預警。

下一步研究工作將進一步融合不同平臺(如鐵路、長途客車、航空等)記錄的任意區域間人群遷徙信息，構建多元出行場景下人群流動蘊含的區域流空間鄰近關系，全面獲取更精細的疫情分布空間異常區域探測結果。另外，需要關聯社交媒體簽到、用戶移動軌跡等個體粒度時空大數據，從輿論擴散的空間分異性、城市內部確診病例活動規律與人群聚集性等多個視角進一步剖析疫情擴散的時空分異特征，并將此作為一種先驗信息嵌入復雜網絡或深度學習模型實現疫情演化過程的動態預測。