新型冠狀病毒肺炎早期時空傳播特征分析*

王聰 嚴潔 王旭 李敏

1)(四川警察學院計算機科學與技術系,瀘州 646000)

2)(四川警察學院,四川警事科學研究院,成都 610200)

3)(四川警察學院道路交通管理系,瀘州 646000)

4)(四川師范大學影視與傳媒學院,成都 610068)

5)(四川師范大學計算機科學學院,成都 610068)

通過最新公布的流行病學數據估計了易感者-感染者模型參數,結合百度遷徙數據和公開新聞報道,刻畫了疫情前期武漢市人口流動特征,并代入提出的支持人口流動特征的時域差分方程模型進行動力學模擬,得到一些推論: 1)未受干預時傳染率在一般環境下以95%的置信度位于區間[0.2068,0.2073],擬合優度達到0.999;對應地,基本傳染數R0位于區間[2.5510,2.6555];極限環境個案推演的傳染率極值為0.2862,相應的R0極值為3.1465;2)百度遷徙規模指數與鐵路發送旅客人數的Pearson相關系數達到0.9108,有理由作為人口流動的有效估計;3)提出的模型可有效推演疫情蔓延至外省乃至全國的日期,其中41.38%的預測誤差≤ 1 d,79.31%的預測誤差 ≤ 3 d,96.55%預測誤差 ≤ 5 d,總體平均誤差約為 2.14 d.

1 引 言

2020年伊始,新型冠狀病毒肺炎(COVID-19)疫情猝然爆發.臨近春節,疫情在春運因素的疊加下從九省通衢的武漢三鎮向祖國大地迅速蔓延.目前國內疫情發展已相對趨于平穩,但國際視角下發展趨勢卻日益惡化: 截至2020年3月9日,疫情已出現在全球102個國家和地區,其中意大利、伊朗和韓國等國家確診人數均已達到數千人.如不采取緊急措施,COVID-19的全球流行趨勢將日趨嚴重.考慮到疫情早期傳播過程客觀上更能反映其傳染特征,為給國內可能產生的局部復發及國際疫情防控協作提供有意義的決策參考,利用當前相對豐富的病例調查數據,挖掘疫情在爆發初期的傳染和域間傳播規律仍有必要.

COVID-19的基本數據估計是研究的基礎.來自《新英格蘭醫學雜志》的研究[1]對COVID-19早期的425名感染者的統計顯示: 疾病平均潛伏期約為5.2 d,平均傳染周期約為7.5 d,基本傳染數(basic reproduction number,R0)約為2.2.Zhou等[2]利用公開新聞報道數據代入倉室模型得到了另一個可能的結果,他們初步估計,COVID-19的 R0值約為2.8—3.3,傳播能力接近SARS.通過對“鉆石公主”號郵輪個案的分析,文獻[3]認為在95%的置信度下,R0的值約為2.28(2.06—2.52).然而該文獻對每日感染人數的泊松分布假設并未得到其他文獻印證,因乘客和工作人員陸續離船和核酸檢測周期帶來的時變因素也不能完全忽略.改進的倉室模型同樣被應用于防控效果評價[4]和拐點預測[5].但傳統的倉室模型并不能夠模擬疫情在地域間的傳播,《科學》雜志的一項研究[6]利用復雜網絡對人類遷徙流動進行了分析,指出疾病的傳播與城市間的“等效距離”,而非地理距離密切相關,而等效距離的核心構成因子是城市間的交通流量.因此,考慮人口遷徙的疫情傳播模型得到了廣泛關注,《柳葉刀》的一項研究[7]利用了航空網絡人口流量數據和騰訊移動數據,通過一種考慮人口遷移的SEIR模型[6,8]首次對本次疫情的時空傳播特征進行建模,得出了 R0的估計值為2.68,并預測疫情在全國的指數級爆發將會在1—2周后出現.然而該研究的主要問題在于,鐵路而非航空才是國內春運期間的核心運力;假定武漢至其他城市的流量相等,與真實的域間交通流特征同樣有著顯著差異.Yang等[9]基于百度遷徙[10]數據和中國疾病預防控制中心(Chinese Center for Disease Control and Prevention,CCDC)每日公布的上報病例數[11],參照SARS的流行病學調查,利用深度學習的方法預測了COVID-19在國內和幾個重點省份的發病曲線,認為疫情在4月底前難以完全平息.該研究的結果緊密依賴于百度遷徙規模指數,而該指數構造方法目前并未公開,與真實人口遷徙的關系并不完全明確.由于流行病學模型的預測結果通常顯著依賴于初始系數的取值,以上研究在參量估值上的差異可能使得使用不同來源的研究結果大相徑庭[12].更重要的問題在于,來自CCDC的針對截至2月11日中國內地所有報告病例的調查[13]證實,前期各項研究采用的上報日期數據與真實發病日期之間在時間和波形上均存在顯著差異.作為真實發病日期的一種近似,確診和上報日期存在的這種誤差對前期各項研究的可靠性也提出了挑戰.考慮到有限的醫療資源、社會承受力和行政成本[14],各級防控措施的介入時點和力度是另一個亟待研究的問題.管控介入時點的一個重要參考是疫情首次到達本地時間[15],Chen等[16]建立了全球流行病的城際傳播模型,利用航空流量數據對疾病到達航空網中任一城市的時間進行了估計.該研究將城市間的交通流量以及疾病傳染率定義為常數.而在真實環境下,傳染率等關鍵參數通常是時變[17,18]的,且對于本次疫情而言,春運因素使得交通流量的常數假定也難以成立.通過以上分析可知,目前的研究集中于疫情基本參數的估計和疫情傳播防控過程的建模和預測,且受限于早期相對匱乏的數據,進一步的有效性交叉印證和誤差校正也有其必要性.

本文通過甲類傳染病防控措施尚未啟動時病例的回顧數據,首次利用病例的發病日期而非確診上報日期,給出了未受干預條件下COVID-19傳染能力估計,并利用一宗極端環境下的個案討論了疫情的極限傳染能力.在此基礎上,首次給出了疫情蔓延至全國各省區的關鍵時點估計: 利用百度遷徙的公開數據和新聞報道數據推斷湖北至各省區的真實人口遷徙流量,并與其他來源的公開報道數據進行了交叉校驗;利用估計的參量值和人口遷徙流量構造了時域差分方程傳播模型,用于推斷其他地區首例病人到達時間,并將預測結果與收集得到的全國各省區首達病例時間做對比,驗證了模型的有效性.本文的結論對復工復產復學過程中疫情局部復發的防控工作可能有所助益,也為疫情在域外的傳播防控提供了參考.

2 疫情的初始傳染情況估計

2.1 基本數據

來自CCDC最新文獻的統計[13],截至2月11日,確診病例共計44672人.從中可整理12月31日前病例發病日期統計如表1.

2020年1月1日以后的精確數據文獻[13]并未給出,本文采用兩種方法進行估計:

1)利用GetData Graph Digitizer軟件抓取文獻[13]公布的每日發病人數統計直方圖,將其轉化為近似數值;

2)利用文獻[19]提出的蒙特卡羅方法,從該文獻擬合的韋布爾分布中采樣得到發病日期與上報日期的時間差,并結合CCDC網站[11]爬取得到的各地上報數據,反推出每日發病人數的估計值.

為校驗方法的可靠性,統計兩種估計方法得到的每日累積發病人數,依據文獻[13]明確給出的1月10日,20日和31日數據推算的累積發病人數和對應日期上報CCDC的確診人數如表2所列.

由于CCDC公布的確診病例統計自1月3日開始,表中12月31日的上報人數以3日人數代替.表2有兩點值得注意.第一,文獻[13]數據與上報數據,特別是疫情前期的數據存在顯著差距.其中一個重要原因是,文獻[13]以發病日期統計病例,而CCDC網站以各地上報日期統計病例.這一差別極有可能增大前期部分研究的誤差.考慮到真實發病日期更具參考意義,在此以文獻[13]數據作為參數估計的基準數據.第二,表2的多數數據為某一區間范圍.這是因為文獻[13]百分比僅保留至小數點后一位,由此帶來的舍入誤差在早期病例數量研究中不容忽視,因此按照四舍五入的原則,反推出單項數據的上下限.從表2可以看出,無論是蒙特卡羅法亦或軟件抓取法,所得數據與文獻[13]給出的當日真實累積發病人數均較為接近,其參考價值可以互相印證.比較而言,軟件抓取法的數值均落入依據文獻[13]反推的區間內;且從2月11日的總計數據來看,軟件抓取法獲得的累積估計發病人數與截至當日的實際累積發病人數誤差僅約0.04%;受2月11日將臨床診斷納入確診依據這一重大變化的影響,12日新增確診病例高達15153人,導致蒙特卡羅法在2月11日的估計誤差達到了54.82%,因此本文采用軟件抓取法的數據為基準.

表1 2019年累積發病人數Table 1.The cumulative number of confirmed cases in 2019.

表2 2020年各時間段累積發病人數Table 2.The cumulative number of confirmed cases in each time slots in 2020.

2.2 一般情況下的傳染率與基本傳染數估計

倉室模型是目前流行病傳播研究的主流方法之一.模型通常包括易感者(susceptible,S)、感染者(infected,I)、潛伏期感染者(exposed,E)和康復者(recovered,R)四類人群,根據考察目的的不同,可分別構成SI,SIS,SIR,SEIR等模型.考慮到疫情初期尚不存在,或極少存在康復者,且潛伏期感染者同樣存在傳染的可能,簡化討論起見,將潛伏者與感染者合并計算,并以SI模型模擬疫情初期發病情況.在此假定下,潛伏期對模型的核心影響可近似表達為時域上的平移,不影響其他參數的實質性討論.令總人口為N;人群中初始感染者人數為 i0;傳染率,即每名傳染者每天使得易感者致病的概率為 β;則在時刻t,感染者的數量 it可表達為

可見這是一個Logistic函數.將首次出現發病病例的2019年12月9日作為坐標原點,并考慮到我國在1月21日將COVID-19列為乙類法定傳染病,并采取甲類傳染病的防控措施,利用2020年1月21日及以前的數據進行估計.將N定義為武漢城市長期居住人口,公開資料顯示約為1418.65萬人.注意到截至12月31日疫情雖已傳播至湖北部分市州,然而由于此時病例主要集中在武漢,且總人口基數足夠大,使用武漢亦或湖北的人口數據為基數對估計結果不產生決定性影響.將表1和表2數據代入(1)式做回歸分析,可得在不同時間點上,β的擬合值,95%置信區間與對應的擬合優度(goodness of fit,R2)如表3所列.

從表3中可以看出,即使在疫情早期數據不夠充分,且受隨機因素影響波動較大的情況下,模型的擬合優度最低亦能達到0.868;隨著日期的推移,可以獲取的數據日漸豐富,傳染率 β 的取值從0.2213逐漸收斂到0.2060—0.2070,擬合優度也隨之上升至0.999.可以認為SI模型能夠有效模擬疫情早期的爆發情況.抽取表3中 β 在幾個時間點上,包括12月31日,1月3日,1月9日和1月20日的取值對疫情做擬合,可得圖1.

表3 傳染率 β 可能的取值Table 3.Possible values of the infection rate.

圖1 不同的 β 取值下實際累積發病人數與預測發病人數的對比(a)2019年;(b)2020年Fig.1.Comparing the actual cumulative number of cases and its estimations according to different β :(a)Year 2019;(b)year 2020.

注意圖1縱軸取對數.為便于觀察,將幾個重要的時間節點數據列于表4中.

表4 重要時間節點的累積發病人數Table 4.Cumulative confirmed cases in key time nodes.

由表4可見,早期數據中存在的較大波動可能導致顯著的計算誤差.而代入1月10日和20日的數據后,無論取上限還是下限所受波動均在可接受的范圍內.因此有理由將該表的結尾,即1月20日的 β 取值和置信區間作為傳染率的有效估計.

一般認為基本傳染數(basic reproduction number)R0的一個可行的計算公式為[20]

其中 Tg為平均感染天數,即病例從具備傳染能力到康復或死亡而失去傳染能力的平均時長.由此可得到 R0在疫情早期的一般估計.文獻[1]給出的Tg可能的取值為7.5 d.《新英格蘭醫學雜志》最近的研究[21]對 1099個病例進行了調查,顯示COVID-19潛伏期中位數和四分位距(interquartile range,IQR)為 4.0 d(2.0—7.0 d);自出現癥狀至發展為肺炎的中位數和IQR為3.0 d(1.0—6.0 d).由于COVID-19在潛伏期同樣具有傳染性,且從IQR來看,這兩個數據均存在一定程度的右偏,有理由認為 Tg的取值為(4.0+3.0)d+右偏量,其后病人應已得到收治.通過不同文獻的相互印證,可以推定 Tg=7.5的取值在當前研究階段是有依據的.將1月21日估計的 β 取值和置信區間代入(2)式,估計對應的 R0為2.5525,95%置信區間為[2.5510,2.6555].

2.3 極端環境下的個案分析

截至目前本文的分析僅涉及了一般情況.全國范圍內的復工復產目前已陸續展開,復學工作也在謹慎研究中.在特殊情況下,如環境更封閉(工廠車間)、人口密度更大(中小學與大中專院校)、人際接觸更頻繁(辦公環境與寫字樓)等,一旦爆發疫情,后果仍不堪設想.因此本文專門針對一宗封閉環境下中小人群聚集性疫情的個案進行分析.

據2月21日新聞報道[22],山東任城監獄發生聚集性疫情.監獄是典型的同時具備封閉環境、高人口密度和人際頻繁接觸等特征的研究樣本,在此將其作為極端情況進行分析.相對于文獻[3]研究的“鉆石公主”號郵輪個案,監獄環境更加封閉,人口密度更高,且任城監獄的個案核酸檢測周期較短,可以將“鉆石公主”號郵輪因檢測周期和乘客在長達10日內陸續離船導致的時變因素最大限度地排除.由于監獄在押人員數量和人員編制不屬于政府公開信息,本文根據公開報道[23]進行側面推算.該報道稱,截至2月20日該監獄確診病例207例,核酸檢測人數為2077人,因此本文假定N=2077,it=207;同時有報道稱2月初個別人出現咳嗽癥狀,因此假定 i0=1,坐標原點為2月1日,可得到t=19.考慮到缺乏數據的交叉印證,在此并未沿用文獻[3]中以泊松分布擬合每日發病人數的方法,而采用了點估計,將t=19時的確診人數 it=207直接代入(1)式進行求解,可得 β ≈ 0.2862;代入(2)式可得對應的 R0達3.1465,顯著高于一般情況下的傳染率.在此條件下,患病人數翻倍時間約為2.42 d.考慮到2077人的數量級與多數中小學或中等規模的工廠接近,這為我們敲響了警鐘,日常生產生活中的疫情防控工作仍不能放松.

2.4 簡要小結

本節估計,在95%置信區間內,未加人工干預時COVID-19的傳染率 β ∈[0.2068,0.2073],其極大似然估計為 0.2070;進而推論 R0在 [2.5510,2.6555]間波動,極大似然估計為2.5525,與學界對SARS較為樂觀的估計相近[24,25];作為一個可能的極限值,在極端利于病毒傳播的環境下,R0可能達到3.1465,與SARS較為悲觀的估計相近[26].雖較文獻[2]的估計略低,但仍可得到相似的推論,即COVID-19屬于傳染性中等偏高的疾病.

3 疫情的時空遷徙特征

回顧SARS及本次疫情的流行過程,一個共同的特點都是從一點爆發,隨著人員遷徙流動迅速蔓延至全國.顯然,如果能夠在疫情傳播的早期切斷域間傳染途徑,經濟和社會承受的壓力會小得多.因此研究疫情的時空遷徙特征,特別是爆發初期的特征十分必要.本節依賴前文的基礎參量估計和百度公開的遷徙數據,構造了簡易的時域差分方程模型,模擬了疫情早期的地域遷徙狀況,并估計了各省級行政區域的首例病人輸入時間.

3.1 遷徙流量的基礎數據估計

百度遷徙網站[10]開放了2020年至今全國所有省區和地級市的遷徙規模指數,以及每個區域的Top 100人口流量百分比.其中遷徙指數和人口流量百分比均分別公布了遷入和遷出數據.通過對該網站的爬取,將元旦至春節百度遷徙規模指數與2019年同期數據進行對比,結果如圖2所示.

考慮歷年春運工作均圍繞農歷新年安排,因此以農歷日期為基準進行統計.最新研究[27]認為雖然在關閉離漢通道前產生短暫的恐慌性遷出人流,但今年武漢春運期間的遷徙流量與往年相比沒有顯著的統計學差異,可以認為離漢人員主體是正常春節返鄉人流.然而由于遷徙指數的構成目前并未公開,有必要對該指數與遷出人口的對應關系進行驗證.

考慮到鐵路運力是春運期間域間流動的主體,且歷年春運中鐵路運力的占比總體穩定,從公開的新聞報道中抽取了武漢市三大火車站(武漢站、漢口站、漢陽站)春節前的發送旅客人數.由于文獻[27]驗證了不同年份春運期間人口流量的相似性,豐富樣本起見,同樣抽取了2019年春運的報道,得到統計結果如表5所列.

限于信息采集手段,表5的數據缺失了部分日期的發送旅客人數.

圖2 武漢“封城”前夕遷徙指數與2019年的對比(a)遷入規模指數;(b)遷出規模指數Fig.2.Comparing the migration index of Wuhan before the spring festival with the same period of 2019:(a)Inner migration index;(b)outer migration index.

如表5所列,遷徙指數與三大火車站發送旅客人數之間的Pearson相關系數達到0.8408,二者呈強正相關;考慮到鐵路的極限運力和反應滯后性(春運期間的車票預訂因素),且1月22日宣布即將關閉離漢通道可能導致恐慌,將該點數據作為異常點剔除,重新計算的Pearson相關系數達到0.9108,故可以認為二者極強正相關.因此可以粗略假定區域j的遷徙指數Mj與實際遷徙人口pj滿足線性關系:

一個重要的問題是參數k的取值.針對湖北省的全部市州和直屬行政區(天門、潛江、神農架林區等),有

來自《荊州日報》和湖北省春運辦的新聞報道稱[28]: “春運首周,···截至 17日零時,全省鐵、水、公、空發送旅客1050.7萬人次”,可知

統計報道同期湖北省所有直屬區域的遷出指數,可得

于是得到比例系數k的估計約為0.100679.

進一步,武漢市春運開始至封城期間的遷出指數合計107.116992,遷入指數合計56.4331536,可估計武漢市在此期間人口凈流出為503.4815萬人.該估計亦可與武漢市政府在1月26日的新聞發布會上公開的遷離人數相互印證[29],可以認為該推論具備一定的參考價值.

3.2 疫情到達時間推演

防控決策時需要回答的重要問題是,疫情將于或已于何時到達本地.在此,提出了一個簡易的時域差分模型,利用遷徙數據,推導疫情首次到達特定行政區域的可能時間.

首先給出SI模型的差分形式:

當相對于 it-Δt和β,N取值極大,即易感者人數占絕對優勢時,有便于計算起見,(7)式可近似為

表5 武漢市三大火車站發送旅客人數與遷出指數Table 5.The Baidu inner migration index and the number of the travelers sent from Wuhan`s major railway stations.

其中Δt定義為一個自然日,即滿足Δ t=1 .當考慮人口遷徙因素時,假定人口遷徙概率與感染概率條件獨立,對于區域j,有

其中ij,t是區域j在第t天的感染人數;σn→j,t-1為第t—1天區域n向區域j流動的人口占區域n當日遷出總人口的比例;Mn,out,t-1為區域n在第t—1天的遷出指數;Nn為區域n的總人口數.由表3可知,由于傳染率的估計會隨著數據的豐富不斷變化,為了擬合這種數據估計上的時變特征,在此將傳染率 β改寫為βt-1,即第t—1天統計發病人數得到的傳染率估計.

根據文獻[13],截至2019年12月31日湖北省共有14個縣區報告病例.考慮到僅武漢市轄下縣區級行政區域就有13個,可以認為在疫情初期,絕大部分病例集中出現在武漢市,因此僅考慮武漢向外省遷移的情況.

疫情首次到達模擬和判決的偽Python代碼算法如下:

在此必須考慮,當疫情前期數據不夠豐富時,傳染率 β 的估計值可能不夠精確,但此時卻是實施決策的關鍵時刻,利用后期得到的傳染率代入預測在決策時點上是做不到的.因此取多個時間點的傳染率預測代入算法1進行仿真.本文取12月31日,1月3日和1月20日 β 的估計值代入計算.取1月3日的原因是,利用12月31日的數值預測,疫情將于1月4日蔓延至省外.因而3日可以作為決策前夕的關鍵時點.得到的結果如表6所列.

表6 省級區域的首例到達時間Table 6.The arrival times of each provinces.

圖3 各省區首達病例時間預測誤差Fig.3.The estimation errors of the first arrival times for each provinces.

搜索各省級區域關于首例確診病例的新聞報道,并與騰訊網“新冠肺炎確診病患活動軌跡”專題[30]提供的病例軌跡數據合并,整理了病例抵達當地的時間填入表6中.表中海南和北京進行了特別標注,原因是: 雖然海南第124號確診病例于1月2日自武漢返瓊,但發病日期為2月5日,期間間隔35 d,超過了文獻[1,21]報道的最長潛伏期,不應認定為輸入病例,故將該省首達日期認定為第157號病例抵瓊日期,即1月13日;北京市病例軌跡數據缺失,但據新聞報道[31,32],該市至遲1月8日凌晨即存在輸入病例,在此取該日期作為首達日期;西藏自治區僅1例病例且到達于1月30日,其統計意義不足,予以剔除.

各省區病例首達時間的預測誤差如圖3所示.就已搜集的信息來看,模型對疫情向省外的蔓延時點的估計較為準確.基于12月31日的數據預測,疫情將于1月4日蔓延至省外.考慮到蔓延前夕的預測對決策最為關鍵,因此以1月3日的預測數據為例分析,疫情首先于1月4日蔓延到河南省,首達省份與模型預測一致,首達時間預測誤差也僅有1 d.模型估計各省首達時間的平均預測誤差為2.14 d,其中有 41.38%的預測誤差 ≤ 1 d,79.31%的預測誤差 ≤ 3 d;除甘肅外,所有省份的首達時間預測誤差均 ≤ 5 d,反映了模型的有效性.考慮到各省最早上報病例日期為1月19日,相對在省外蔓延的確切日期滯后16 d,采用模型預測得到的理論預測值顯得更有價值.即使擇取12月31日和1月20日的模型參量估計也可達到較好效果,其中預測疫情蔓延至省外的時間誤差分別為1和2 d,首達省份預測無誤;平均預測誤差也僅為2.31和2.48 d.表6可能反映的另一個沉痛的現實是,在迅猛發展的疫情面前,強力交通管制措施采取的時間同樣顯得滯后.如果根據模型的預測,將強力交通管制措施提前到1月4日,即模型預測的即將蔓延至省外的時間施行,則可有效阻斷疫情在除河南外全部省份的蔓延.

4 結 論

本文利用最新的流行病學調查,取疫情尚未得到有力干預時的數據代入動力學模型進行擬合,給出了傳染率和基本傳染數的一般估計.為模擬復工復產復學過程中可能出現的聚集性疫情特征,利用公開報道的極端環境個案,給出了可能的傳染率和基本傳染數上限.在此基礎上,綜合流行病學調查數據、互聯網爬取數據和公開新聞報道數據,在佐證百度遷徙規模指數參考價值的前提下,估計了疫情發源地遷至各地的絕對人口數量,代入設計的時域差分模型,模擬了疫情初期在地域間的傳播過程,給出了疫情首次到達各省級區域的時間估計,并與新聞報道數據做比對,驗證了模型的有效性.

本文的創新在于: 首次以流行病學調查中公布的發病時間,而非確診和報告時間為基準進行疫情推演.相較于兩個時間節點間約10 d的錯峰,這一選擇顯然更具合理性;給出了已得到初步應用的百度遷徙規模指數的有效性驗證,側面證實了它與人口遷徙規模的粗略線性關系;利用人口遷徙數據,首次給出了疫情到達全國各省的時間預測.有限度使用新聞報道數據以改善疫情初期科研數據尚不夠豐富的現狀也是本文的一個特點.

本文還存在以下不足,將在后續工作中改進:本文依據的流行病學調查數據更多依賴于患者個人回憶,可能存在時間偏差,需要更加精確的數據進行進一步驗證;本文提出的方法對疫情爆發后強力防控措施介入條件下的推演偏差可能較為顯著,應重點思考優化方案;如何對各級力度的疫情防控措施進行量化分析與效能評價,在有效控制疫情的前提下做到松緊適度,降低對經濟社會的負面影響,也是一個值得研究的問題.