新冠肺炎傳播預測優化:基于少量早期數據的多階段SEIRr模型

賽 斌, 宋 兵, 譚索怡, 歐朝敏, 周 濤, 張 偉, 呂 欣

(1.國防科技大學 系統工程學院,湖南 長沙 410073; 2.電子科技大學 大數據研究中心,四川 成都 611731; 3.四川大學華西醫院 華西生物醫學大數據中心,四川 成都 610041)

1 引言

2019年12月8日，湖北省武漢市發現不明原因引起的肺炎，后被證實該疾病由新型冠狀病毒(SARSCoV-2)引起，2020年2月12日世界衛生組織(WHO)將由這一類病毒引起的疾病命名為COVID-19。與SARS不同，COVID-19潛伏期具有一定的傳染性，并且存在大量無癥狀感染者，這給疫情防控工作帶來了巨大挑戰。截至2022年6月10日，全球COVID-19感染人數已超5.117億例，累計死亡人數超632.8萬，眾多國家的社會和經濟因此遭受了巨大損害，人類社會也因此面臨著前所未有的挑戰。當前我國疫情處于平穩態勢，但仍存在著反彈風險，堅持“外防輸入，內防反彈”的總策略，減少疫情對社會各方面的影響，研究如何根據疫情形勢制定防疫措施具有重要意義。在面對突發傳染病時，掌握疫情發展態勢，合理有效的防疫政策能夠控制疫情發展，將影響降到最低，也能夠最大程度地保障社會各方面的正常運轉。由此可見，準確合理地預測傳染病發展趨勢，提前規劃應對措施或制定防控策略，才能最大化地實現“動態清零”策略。

自疫情暴發以來，基于傳染病傳播動力學模型分析與預測COVID-19的傳播特征及發展趨勢成為了研究人員的主流研究方法。在疫情發生初期，Wu等通過使用中國疾病預防控制中心發布的2019年12月31日到2020年1月28日期間醫學統計數據，結合官方航空指南的每月航班預訂數量與騰訊位置大數據得到的人口流動數據構建SEIR模型，估計得到COVID-19的基本再生數為2.68(95%Crl 2.47～2.86)，疫情倍增時間為6.4天，推斷COVID-19在中國內地主要城市的感染人數將呈現指數增長并預測武漢市將在2020年1月25日會有75815位新冠肺炎感染患者。Zhou等基于2020年1月25日前的已有的醫學統計數據及國際同行估計的感染人數，參考SARS流行病的參數構建SEIR模型估計COVID-19基本再生數，得出COVID-19基本再生數在2.9～3.2之間的結論。Shen等基于新冠疫情的傳播機制，考慮了防控措施，建立傳染病傳播動力學模型，預測出COVID-19暴發規模將在2020年2月5日到達峰值。Tang等基于似然估計和模型仿真的估計結果表明，武漢市的控制繁殖數可能高達6.47(95%CI 5.71～7.23)，該文章也預測了在最嚴格的管控措施下，中國內地將于2020年1月23日起的兩周內達到疫情峰值最小。王聰等利用文獻統計的最新流行病學數據估計疫情早期基本再生數，結果顯示武漢市在2020年1月21日前后的傳染病參數

R

為2.55(95%CI 2.55～2.65)。

可見，基于病理學的傳染病傳播動力學模型的COVID-19基本參數信息分析與疫情發展態勢預測研究對新冠肺炎流行因素的掌握以及防控策略和措施的制定具有重要作用，并且對于早期疫情的控制具有指導性作用。隨著互聯網及大數據等技術的迅猛發展，也有學者將SEIR模型結合人口遷徙網絡或機器學習等方法研究傳染病的傳播機制，例如，Yang等基于對潛伏期患者具有感染性的考慮優化了SEIR模型，并結合LSTM深度學習方法進行建模，預測得到中國內地的疫情規模將于2020年2月28日到達頂峰，然后在4月底逐漸平緩。Jia等利用手機定位大數據分析人口流動時空特征，然后基于疫情暴發前武漢輸入到全國各地的人口流動數據，構建了“時空風險源模型”，該模型最多可以提前兩周預測出新冠疫情暴發的時間、強度和地理分布。譚索怡等對北京市COVID-19的密切接觸者進行數據分析時發現，新增密切接觸者的增長率與5～6天后的新增確診病例增長率變化趨勢一致，由此推斷中國內地新增確診人數在2020年2月13日左右開始出現下行。但是，由于在疫情暴發初期醫學檢測能力不足,以及對新型冠狀病毒肺炎傳播機制、傳播能力認知的不足，導致了研究人員在進行COVID-19疫情動力學建模研究時缺乏足夠的醫學統計數據以及對新冠肺炎傳播機制的錯誤認識，因此難以準確估計疫情發展態勢。

隨著COVID-19相關特性研究的不斷深入，有研究表明COVID-19的潛伏期患者也具有感染性，并且存在無癥狀感染者，而大量已有的建模預測研究忽略了這些特征。所以，考慮到各種實際情況時，一般都需要在傳統SEIR模型上進行必要的修改以貼合新冠肺炎實際傳播情況，例如，北大統計學團隊的Yan等細分人群種類建立vSIADR模型估計有效再生數并對各國疫情進行波次劃分，研究了25個國家的控制措施對COVID-19傳播的效果；López和Rodo等改進SEIR模型建立了SEIQRDC模型解釋了潛伏期的感染傳播情況。魏永越等開發SEIR(+CAQ)動態模型較為準確地預測了中國內地主要疫區的累計確診病例數。然而，由于沒有全面考慮到疫情傳播特征與機制以及對疾病的認識的動態變化造成的防控、診斷及治療的優化等問題，以上研究方法存在或預測準確率低，或沒有全面考慮COVID-19傳播特征，或用于擬合數據多、預測周期短等問題。因此，在面對重大新發突發傳染病時，構建更為合理的傳播動力學模型，采用多階段建模方法預測疫情發展趨勢并分析傳染病傳播機制勢在必行，這將有利于相關部門制定更為科學合理的防控措施,各地醫療衛生機構提前做好應急準備工作，有效提升疫情防控效率。

為解決上述問題，本文在綜合考慮COVID-19的傳播特性的基礎上，構建了多階段SEIRr模型預測COVID-19傳播趨勢，并基于國內外不同階段真實醫學統計數據進行實證分析，驗證了模型的有效性、可靠性及普適性。其主要創新點如下：(1)充分考慮了COVID-19傳播特性，即潛伏期患者傳染性與自治愈等特點，增強了傳染病傳播動力學模型的穩健性。(2)創新性地應用了梯度下降算法進行參數估計，提升了傳播動力學模型參數估計的有效性。(3)增加了對模型輸入初始醫學統計數據初始值的估計，降低了模型對前期數據的質量要求。本文提出的方法和模型具有很好的可解釋性，可以為相關部門制定科學合理的疫情防控措施提供一定的理論支撐。

2 模型與仿真

2.1 SEIRr模型建模

經典的SEIR模型將人群分為了

S

(易感者，未感染但存在感染風險人群)、

E

(潛伏者，感染但不不具有傳染性的未發病狀人群)、

I

(感染者,已發病的具有傳染性人群)以及

R

(康復者，指病愈或因病死亡人群)4類。然而，據現有研究表明：COVID-19的潛伏期人群具有傳染性，且由于疫情初期檢測力度不足導致部分無癥狀患者和自治愈人群被選擇性忽略了；另外，現實情況下并不是每個人與病毒攜帶者接觸后就一定會被感染。因此，

E

類人群中存在一定比例的無癥狀感染者、自治愈人群以及無效接觸者，他們將直接從潛伏期人群中移除出去。

(1)

根據此次COVID-19疫情的特征，本文做出如下假設:

(1)假設考察地區的人口總量

N

不變，即不考慮人群生死及遷移等種群動力因素。

(2)假設考察期間，病毒不發生變異，疫情防控政策維持相對穩定，即保證傳染率不發生變化。

(3)假設COVID-19的傳染率為

β

，則在

t

時刻單位時間內被感染者感染的人數為

β

·

I

(

t

)

S

(

t

)

/N

(

t

)。(4)考慮潛伏期人群具有傳染能力，其傳染能力為感染者的

k

倍(

k

<1)。(5)假設

S

與病毒攜帶者接觸后可以轉變為無癥狀、無效接觸、自治愈和感染但未被確診四類人群中的一種，其中前三類人群將以速率

μ

從

E

中移除成為

r

類，感染未確診人群將以速率

α

轉為感染者

I

，最終被移除成為

R

類。則在

t

時刻由感染者轉為康復者的人數

R

(

t

)=

γ

·

I

(

t

);由潛伏者轉為移除者的人數

r

(

t

)=

μ

·

E

(

t

)。

基于以上假設，本文構建了考慮疫情防控措施影響的SEIR模型——SEIRr模型，圖1為SEIRr模型的示意圖。模型的微分方程表達式如公式(1)所示，微分方程中的各個參數的取值參考及其具體含義參照表1。

圖1 SEIRr模型示意圖

2.2 參數設計與估計方法

利用傳染病傳播模型進行仿真預測的必要前提是得到有效的參數值，在本模型中需要求解的參數包括

β

、

k

、

μ

、

α

、

γ

。由于

α

和

γ

分別是潛伏期患者發病速率(平均潛伏時長

T

的倒數)與平均康復速率(平均治愈時長

T

的倒數)，二者受醫療資源與診療水平限制，本文假設短時間內二者為常量。據Li等，Guan等估計，

T

與

T

分別約為5.4天與13天，故

α

=1

/

5.4，

γ

=1

/

13。因此，待估參數還有

β

、

k

、

μ

。此外，為進一步提高SEIRr模型的預測效果，本文在估計上述未知參數的同時通過梯度下降算法自動調整微分方程輸入初始值

S

(0)、

E

(0)、

I

(0)、

R

(0)，從而有效克服因疫情暴發初期對新發傳染病認知不足造成的數據統計偏差和數據不足等問題。

為有效、準確估計模型參數，本文采用批量梯度下降算法進行參數估計，各參數值在梯度下降算法擬合優化過程中損失函數最小時取得。梯度下降算法的損失函數如下

(2)

其中代表實際值與預測值間的均方誤差為L2正則項(可以有效避免過擬合問題)，

θ

代表

β

,

k

,

μ

,

α

,

γ

等待擬合參數。

3 實證分析

3.1 數據來源

本文所采用的用于模型擬合與預測的真實醫學統計數據包含國內和國外兩部分。其中國內數據獲取于中國國家衛生健康委員會官方網站及各地方衛生健康委員會官方網站的每日疫情通報數據，國外數據主要來源于美國霍普金斯統計的各國每日累計確診人數、累計死亡人數以及累計治愈人數等流行病學數據。由于本文研究模型的假設條件之一為傳播環境穩定，故武漢數據取2020年1月23日宣布“封城”以后進行研究，其他疫區的醫學統計數據均取自于2020年1月24日全國各地實施公共衛生安全一級響應措施以后。同理，本文所研究的國外的醫學統計數據也均取自于傳染病傳播環境相對穩定之后的數據。用于模型擬合的數據為疫情前期10天左右的每日現有確診人數，其余數據用于檢驗模型擬合效果。

3.2 武漢疫情預測分析

首先，考慮到疫情初期病毒檢測水平的不足而造成確診病例數堆積，2020年2月12日中國衛生健康委員會首次將臨床診斷病例納入統計導致武漢市當天累計確診人數激增，這使得武漢市在該日前后兩天統計數據出現明顯斷層。為消除數據問題造成的預測偏差，本文對2020年2月12日新增加的臨床診斷病例進行平均平滑處理，將其平均分配到1月23日至2月12日期間，圖2a為數據處理前后的對比圖。然后，我們利用梯度下降算法迭代運算5000次對平滑處理后得到的每日確診人數進行擬合，估計相應的傳播動力學參數。SEIRr模型參數初始值設定以及通過批量梯度下降算法擬合得到的結果如表1所示。

表1 武漢市SEIRr模型參數設置詳情

詳情參數含義初始值參數設置依據擬合結果MCMC抽樣N艙室中的接觸人數2000000文獻[28]2001272-S(0)初始易感染人群1997455S=N-E-I-R1998364-E(0)初始潛伏期人群1500估計值1502-I(0)初始感染人群991湖北省衛健委數據再平滑處理1349-R(0)初始移除類人群(包括治愈與死亡)54湖北省衛健委數據57-β易感者與感染者的接觸速率4.5按經驗隨機取值4.3081.023kk·β表示易感者與潛伏期人群的接觸速率1按經驗隨機取值0.2020.015μ潛伏期人群中的無效接觸、無癥狀感染、自治愈人群轉化為r的速率2.9按經驗隨機取值3.2360.079α潛伏期人群被檢測為感染者的檢測速率1/5.4文獻[22,23]1/5.41/5.4γ感染者被治愈速率1/13文獻[23]1/131/13

武漢市現有確診人數的預測如圖2b所示，圖中綠色圓圈為參數擬合用數據，其余真實數據(黑色圓圈)作為對比評價SEIRr模型預測效果，最終得到了

R

=0.975的擬合優度。從擬合得到的結果可見：①武漢市每日現有確診人數的預測結果與實際基本符合；②潛伏期人群相較于感染期人群的傳染率不高，約為感染者的20

%

左右，這與中國-世界衛生組織新型冠狀病毒肺炎(COVID-19)聯合考察報告所描述的“感染者在出現癥狀前1～2天才能檢測到核酸病毒判斷潛伏期患者在出現癥狀1～2天才具有傳染性”的研究是相吻合的；③從圖2c可以看出后期治愈速率是大于1

/

13的，從而導致疫情進度存在10天左右的滯后，這與醫療資源與診療手段的更新相關。

圖2 武漢市疫情預測

從結論③中可以看出，隨著各界對COVID-19認知水平的提高，武漢市針對COVID-19的診斷與治療效果存在顯著提升，然而在新冠肺炎暴發初期發表的論文基本都忽略了該因素對COVID-19發展態勢的影響，這也導致其預測COVID-19的發展趨勢與實際情況存在較大的偏差。本文考慮改變后期(2020年2月27日以后)參數進行仿真模擬預測疫情發展態勢以驗證我們的推測，當后期治愈速率更新為1

/

10時的預測效果最佳，其擬合優度提升到

R

=0.996,和圖2d所示。為進一步驗證SEIRr模型的預測效果和批量梯度下降算法可靠性，本文使用控制變量法進行了必要的對比試驗。首先進行了馬爾科夫蒙特卡洛方法(MCMC)和批量梯度下降算法的參數估計效果對比實驗：為保證與批量梯度下降算法參數一致，MCMC確定模型輸入初始值以及參數

α

和

γ

與梯度下降算法一致的基礎上，使用了8000次NUTS采樣，最終確定

k

、

μ

、

β

三個參數值，預測值如表1所示。結果表明，與批量梯度下降相比，MCMC僅能估計三個參數，且估計精度與抽樣次數相關，抽樣次數越多則精度越高；估計參數的預測效果如圖2(f)所示，批量梯度下降方法的結果顯然優于MCMC方法。然后是SEIR和SEIRr的模型預測效果對比實驗：在該實驗中兩模型的參數估計方法均采用批量梯度下降算法迭代運算5000次對武漢市1月23日至3月2日期間的每日確診人數進行擬合并得到模型參數，然后將參數代入相應模型進行仿真實驗，預測效果如圖2(e)所示。結果表明，SEIR模型存在疫情峰值時間預測偏早，疫情峰值偏高等問題，SEIRr模型預測效果顯然比傳統SEIR模型更好。

3.3 國內主要地區疫情預測

為了驗證模型的普適性，本文用同樣的方法對國內的廣東、河南、湖南、重慶、浙江等地區在2020年1月24日至2020年2月3日期間，以及上海在2022年3月11日至2022年6月10日期間每日現有確診人數進行擬合，從而得到相應的傳染病傳播動力學方程的參數。與2020年原始新冠病毒特性不同的是，此輪上海疫情是由變異毒株——“奧密克戎”引起，據復旦大學研究成果表明，其平均潛伏期時長和在院治療時長分別為1.2天和5.6天。因此，上海疫區的參數

α

與

γ

分別設置為1

/

1.2和1

/

5.6。通過仿真預測得到的結果如圖3中的紅色實線所示，各地區擬合優度也都在90

%

以上，與武漢市情況類似，國內其他地區的后期治愈速率同樣得到了提升，治愈速率更新后的擬合效果如圖3中的藍色實線所示(其中2022年上海疫情僅用單階段模型即可得到相對準確可靠的擬合結果)。這表明，在傳染病傳播環境穩定的條件下，SEIRr模型的預測精度和準確性具有可靠性；隨著人們對病毒認識的提高和診療水平的提升，使得COVID-19患者的平均治愈速率能夠提升并維持相對穩定。

圖3 中國主要地區疫情發展趨勢預測

3.4 國外疫情趨勢擬合與預測

為了進一步驗證模型的普適性，本文挑選了傳染病傳播環境都是較為穩定的意大利、德國、日本、瑞士四個國外疫區疫情發展趨勢進行預測。其中意大利和德國均采取了類似于中國的封城措施，意大利在2020年3月4日宣布全國學校停課，3月8日實施“封城”措施，德國則是在3月16日由總理默克爾正式宣布實施“緊急防疫措施”。而日本則是在4月7日實行“軟封城”措施，該措施更多地是依靠居民的自覺性。瑞士實行的是“群體免疫政策”，政府沒有采取任何強制措施對COVID-19進行限制。

使用同樣的方法對國外數據進行參數估計并預測疫情發展趨勢，參數設置與預測結果如表2與圖4所示。實驗結果顯示，各國實際與預測每日現存確診人數之間的擬合優度均在0.97以上。這表明，基于前期的公布的少量醫學統計數據，利用梯度下降算法對SEIRr模型求解，也可以有效預測國外疫情發展趨勢。

表2 國外疫情SEIRr模型擬合結果

國家βkμ第一階段α第二階段α第一階段γ第二階段γ德國3.6320.6954.2131/5.41/5.41/141/14日本3.3570.03033.1081/5.41/5.41/141/14意大利3.0250.4163.0171/121/101/151/13瑞士4.7660.4123.0861/101/81/161/13

進一步，本文對比了上述四個國家擬合所得的傳播動力學模型參數。可以看出，實施“群體免疫”政策的瑞士傳染率最高；而意大利的防控措施效果較好，其傳染率是所有國家中最低的；德國和日本的傳染率較為相近；意大利與瑞士的檢測速率與治愈速率均較小，本文推斷這是由于意大利感染人數較多所致，瑞士則是因為實行“群體免疫”政策而導致的檢測水平滯后，康復周期變長；但兩參數值在第二階段是有增大的，這進一步佐證了“隨著對疾病認知的提高，醫療資源也必將得到優化”的觀點。

圖4d為意大利疫情趨勢預測結果，該結果所用的擬合數據為3月21至4月5日之間的每日現有確診人數。為說明使用近期數據能夠更好地預測疫情發展的變化趨勢，本文使用不同時期的數據對意大利疫情趨勢進行預測，其結果如圖5所示。由圖5可見，隨著參數學習所用的醫學統計病例數所處時間段的推移，SEIRr模型預測準確率不斷提升；從總體上來看，前期預測結果較實際醫學統計數據高但各時間段所推測出的疫情峰值與實際日期均高度重合。這有可能是由于意大利在疫情初期存在檢測力度不足而導致前期公布的每日現有確診人數偏低。

圖4 國外部分國家疫情發展趨勢預測

圖5 不同擬合數據下的意大利疫情預測結果

4 結論與展望

本文結合COVID-19傳播特征提出了一種考慮潛伏期患者傳染性、感染期人群自治愈能力、無癥狀感染者以及無效接觸人群的SEIRr模型。模型基于疫情早期少量醫學統計數據，使用批量梯度下降算法估計傳播動力學參數及初始感染人數從而預測疫情發展趨勢。研究得到以下的結論：(1)本文的各項實驗結果表明，通過估計模型初始值可以有效克服因早期醫學統計數據誤差以及對傳播特征的錯誤理解等原因而造成預測偏差問題。(2)相較于傳統SEIR模型，SEIRr模型假設更貼合COVID-19實際傳播情況，實現了對國內外多個地區的疫情發展趨勢的精準預測，其擬合優度提高到了0.95以上。(3)相較于基于統計學的MCMC抽樣的參數估計方法，梯度下降算法具有更好的擬合效果，能夠為其他傳播動力學參數估計提供有效借鑒。

基于本文相關研究結果，得到的傳染病防控管理啟示有：(1)從源頭抓起，堅決管制隔離病毒攜帶者和密切接觸者，防止疫情擴散。(2)降低傳播率

β

，輿論引導民眾做好個人防護措施，政府加大監管力度，切斷病毒傳播途徑。(3)提高參數

α

的值，在條件允許的情況下，適當增加核酸檢測強度，及時發現病毒攜帶者，盡早隔離。(4)提高治愈速率

γ

，加大相關科研的投入，研發疫苗和治療藥物，提升醫療資源配置，以應對大規模的疫情反彈風險。