基于階段式SIR-F模型的新冠肺炎疫情評估及預測

凡友榮 楊 濤 孔華鋒

1(公安部第三研究所 上海 201204) 2(武漢商學院 湖北 武漢 430056)

0 引 言

根據世界衛生組織公布的數據，截至2020年5月27日，全球新冠肺炎確診病例已達548萬，超過10萬確診人數的國家包括美國、巴西、俄羅斯、英國、意大利、法國、德國、土耳其、印度、伊朗、秘魯。其中美國、巴西目前的每日新增確診數超過1萬人；意大利、法國等國家的疫情趨緩。由于印度存在人口密集、醫療條件差等問題，世界衛生組織表示，人類對抗新冠疫情能否取得決定性勝利，未來很大程度上取決于印度控制該病毒的能力。針對全球各國不同的疫情狀況，分別進行現狀分析及防控措施的影響研究至關重要。

SIR模型、SEIR(Susceptible-Exposed-Infected-Recovered)模型是傳染病動力學研究的經典模型，目前已有多項基于SIR、SEIR模型的新冠疫情研究。尹楠[1]基于SIR模型進行了有限區域內新冠病毒傳播的模擬分析。梅文娟等[2]將極限學習機與SIR模型相結合構建了極限IR實時預測模型。盛華雄等[3]將SIR模型與Logistic模型相結合，對不同防控措施進行了模擬和對比。喻孜等[4]使用基于時變參數的 SIR 模型預測得到了疫情拐點和最大確診數。Zhou等[5]基于SEIR模型計算得到了COVID-19的再生數。林俊鋒[6]在SEIR模型的基礎上增加了“隱形傳播者”，對隱形傳播者的數量進行了預測。范如國等[7]基于SEIR模型對3種不同潛伏期下的疫情進行了拐點預測。

本文使用的SIR-F模型在SIR模型的基礎上增加了未確診的感染者(具有傳染性)，這一點與SEIR模型相同；本文模型還將“移除者”(Removed)細分為“治愈者”(Recovered)和“死亡者”(Fatal)，從而能夠進一步分析治愈人數、死亡人數的變化趨勢。此外，考慮到在數據擬合的實踐中存在過擬合的現象，并且實際的疫情發展過程分為多個階段，每個階段的趨勢特征各不相同，因此本文先進行階段劃分，再將疫情數據分階段地進行擬合分析，從而提高分析結果的準確性。

1 技術路線

本文的技術路線如圖1所示。首先計算各國新冠肺炎確診數量的增長因子，根據該因子將所有國家分為3個增長類型：“停滯”式、“爆發”式、“波動”式；然后在“停滯”增長態勢的國家中選擇典型例子進行案例分析，由于意大利是全球重點疫情國家之一，且其防控工作已取得顯著成效，因此選擇意大利為例，對其疫情發展過程和防控措施進行深入分析；最后針對“爆發”式增長的國家，從中識別目前爆發態勢最嚴峻的國家，并針對其進行疫情的深入分析，提出防控措施建議。

在案例分析中，首先使用Python調用Prophet的“突變點自動監測”功能對疫情數據進行階段劃分，可用于疫情數據突變點的自動識別。基于SIR-F模型進行分階段的數據擬合、參數估計，對比不同階段的參數，從而發現各階段疫情的變化情況，然后結合防控措施的實施時間，研究防控措施對疫情發展產生的影響。對于“爆發”式增長的國家，在模型擬合的基礎上進行疫情預測，并且基于不同的感染系數進行預測比對，從而模擬不同強度防控措施下的疫情發展趨勢。

2 模型構建

2.1 增長態勢識別

本文通過增長因子來評估 “新冠肺炎”確診人數增長的態勢，增長因子的計算公式如下：

(1)

式中：GF為增長因子；ΔCn為第n天確診人數增加的數量。計算每個國家近7天的GF，進行增長情況的分類：

(2)

GC為增長態勢所屬類別，GC=1代表“爆發”的增長態勢，GC=-1代表“停滯”的增長態勢，GC=0代表“波動”的增長態勢。

2.2 SIR-F模型

由于COVID-19疫情比較特殊，有一部分感染者在死亡前仍未正式確診，將這類情況增加到經典的SIR模型中，構建SIR-F模型。SIR-F模型中從未感染者轉變為感染者、治愈者、死亡者的過程如圖 2 所示。

圖2 SIR-F模型中感染者的轉化過程

圖2中，S代表未感染者，S*代表未確診的感染者，I代表感染者，R代表治愈者，F代表死亡者，α1為S*的死亡系數，α2為I的死亡系數，β為感染系數，γ為I的治愈系數。累計確診數量C=I+R+F。

根據上述模型設計，構建其動力學方程如下：

(3)

(4)

(5)

(6)

式中：N為總的人口數量，N=S+I+R+F；T為疫情開始的總時長。為求解方便，將上述函數方程進行無量綱處理，設(S,I,R,F)=N×(x,y,z,w)，(T,α1,α2,β,γ)=(τt,θ,τ-1k,τ-1ρ,τ-1σ)，得到以下微分方程組：

(7)

在Python(版本3.7.6)環境下，使用超參數自動優化框架Optuna對上述微分方程進行求解，從而得到各個參數(θ,k,ρ,σ)的估計值。

2.3 防控措施的影響函數

各國政府通過疫情防控措施改變民眾的行為，從而控制疫情的發展。本文引入感染系數β受民眾行為影響的函數，如式(8)所示。設感染系數β的影響因子包括：每分鐘內未感染者與感染者在外出過程中接觸的次數c，每周內未感染者外出的天數g，唾液中病毒存在的可能性v，有效佩戴口罩的人數比例m，洗手降低的感染概率we，其他影響因素δ。

(8)

式(8)與SIR-F模型相結合，能夠模擬居民的外出、佩戴口罩、洗手等行為對疫情發展的影響。

3 模擬仿真

3.1 數據來源

本文基于的全球疫情數據來源于約翰·霍普金斯大學的系統科學與工程中心(https://systems.jhu.edu/ )，包含了全球217個國家在110天內(2020年1月22日至2020年5月11日)的累計確診人數、現存確診人數、死亡人數、治愈人數，總數據量為23 804行。各國防控措施數據來源于牛津大學的流行病預測研究組織(http://epidemicforecasting.org/)，包含了96個國家在2019年12月18日至2020年4月15日期間采取的管控措施，總數據量為1 633行。

3.2 全球疫情的增長態勢分析

將全球疫情數據導入式(1)和式(2)，計算得到96個國家的疫情增長態勢所屬類別GC，結果示例如表1所示。

表1 全球各國疫情增長態勢的分類結果

處于“爆發”增長態勢的包括印度、俄國、愛沙尼亞、巴西、肯尼亞等53個國家。將這類國家的疫情數據聚合，增長趨勢如圖3所示。

圖3 “爆發”式增長國家的疫情數據變化趨勢

處于“停滯”增長態勢的包括意大利、日本、利比亞、安哥拉、尼加拉瓜、越南等28個國家。將這類國家的確診數據聚合，增長趨勢如圖4所示。

圖4 “停滯”式增長國家的疫情數據變化趨勢

處于“波動”增長態勢的包括加拿大、以色列、美國、英國等138個國家。將這類國家的確診數據聚合，增長趨勢如圖5所示。

圖5 “波動”式增長國家的疫情數據變化趨勢

觀察圖3至圖5可知，“爆發”式呈現迅速上升的趨勢；“停滯”式的增長態勢已基本得到了控制；“波動”式的未來發展方向尚不明確。“爆發”式和“波動”式增長的國家是全球疫情控制的重點，因此有必要研究“停滯”式增長的國家所采取的防控措施，為其他國家的疫情防控提供有效建議。

3.3 “停滯”式增長案例分析

1) 意大利疫情趨勢分析。意大利目前的疫情態勢如圖6所示。根據式(1)計算得到意大利確診人數的增長因子變化態勢如圖7所示，其增長因子GF<1的持續天數為18天。

圖6 意大利疫情數據趨勢圖

圖7 意大利增長因子變化趨勢圖

使用Prophet 對意大利的累計確診數據進行突變點的識別，如圖8所示。

圖8 意大利累計確診數據的突變點識別

由圖8可知，存在3個突變點：“2020-03-17”、“2020-04-05”、“2020-04-25”，則意大利的疫情發展過程可以劃分為4個階段，如表2所示。每個階段的變化規律有明顯的變化，因此下文將分階段地進行擬合分析。

表2 意大利疫情發展的4個階段

2) SIR-F模型擬合及參數估計。將4個階段的真實數據分別導入SIR-F模型進行擬合，擬合效果如圖9所示。

圖9 意大利4個階段感染人數的擬合效果

對式(7)進行求解，計算得到各階段SIR-F模型的參數估計值，如表3所示，由該表可知參數ρ下降了，即感染系數β下降了，并且有效再生數RT明顯下降，這說明意大利的疫情控制已取得了一定成效。

表3 意大利4個階段的參數估計值

3) 感染率降低的原因分析。意大利在3月采取了多項限制居民出行的管控措施，如表4所示。

表4 意大利3月采取的限制出行措施

考慮到政府的管控措施從開始執行到產生實際的疫情控制效果存在時間延遲，因此認為第1階段為實施管控措施之前，第2、第3、第4階段為實施管控措施之后。因此本文假設居民在第1階段的出行情況與平時相同，設意大利民眾平時每周出行5.6天，即g1=5.6。

根據2020年3月13日報道，由于出行限制、公共活動場所的關閉，使得居民相互之間的潛在接觸次數降低了19%[8]，因此將采取管控措施之前的c取值為1，采取管控措施之后的c取值為0.81。

設出行限制措施只影響了式(8)中的g和c，其他值不變，得到：

(9)

將g1=5.6，c1=1，ci=0.81(i=2,3,4)，以及表3中各階段的ρ值代入式(9)，計算得到g2=2.92，g3=1.06，g4=0.51。可見民眾每周的出行天數已經下降到了0.51天。

上述分析過程從實際的疫情數據反推得到了實施管控后的每周出行天數，可知民眾的每周出行天數減少了90.9%，使得感染系數降低了92.6%，證明了限制出行的管控措施對新冠病毒的擴散起到了明顯的抑制作用；同時證明了本文的分析模型能夠較靈敏地模擬現實中防疫管控下新冠肺炎的傳播情況。

除了上述管控措施，意大利的疫情防控和診療工作均參照了中國方案，這也是意大利疫情控制效果較好的重要原因。意大利從5月4日已進入抗疫和恢復經濟的并行階段，新聞顯示62%的意大利人對于復工感到害怕，因此出行習慣的恢復還需要時間。若出行量逐漸增加，應從其他方面控制感染率，例如提高戴口罩、勤洗手習慣的普及率等，從而防止疫情再次反彈。

3.4 “爆發”式增長案例分析

印度的增長因子GF>1的持續天數為全球最長(66天)，其增長態勢極其嚴峻，因此下文將印度作為“爆發”式增長的典型例子進行深入分析。

1) 印度疫情趨勢分析。印度目前的疫情態勢如圖10、圖11所示。

圖10 印度疫情數據趨勢圖

圖11 印度增長因子變化趨勢圖

對印度的累計確診數據進行突變點的識別，如圖12所示，存在2個突變點(“2020-04-03”、“2020-04-28”)，則印度的疫情發展過程可以劃分為3個階段，如表5所示。

圖12 印度累計確診數據的突變點識別

表5 印度疫情發展的3個階段

2) SIR-F模型參數估計。與意大利案例的分析方法相同，將印度各階段的疫情數據分別導入SIR-F模型，參數估計結果如表6所示，可以看出參數ρ和RT存在一定程度的下降，但第3階段的值仍較高，說明目前印度的疫情擴散仍然較為迅速。

表6 印度3個階段的參數估計值

3) 印度疫情預測。以第3階段的參數值進行疫情預測，結果如圖13所示，印度感染數量將在2020年10月16日達到最高峰，最高值為3.84億人。根據紐約時報3月27日發布的一篇有關印度疫情的文章，作者Laxminarayan教授預測印度感染人數將達到3至5億，這與本文的計算結果基本符合。

圖13 目前管控措施下的印度疫情預測

若印度采取更強有力的防控措施，例如為印度貧民地區提供更多口罩、酒精等醫療物資，提高核酸檢測的覆蓋率，更精準地進行病患隔離等，將能有效降低感染系數。假設印度的ρ值從目前的0.041降低為0.026，再次進行疫情預測，結果如圖14所示。在ρ為0.026的假設下，印度感染數量的最高值降為2.01億人，最高值到來的時間推遲到了2021年3月7日。

圖14 假設ρ值降至0.026的疫情預測

通過上述對比分析可知，若印度增強防控措施的力度，將感染系數在目前的基礎上降低36.59%，感染人數的最高值即可降低47.66%，到達最高感染人數的日期將推遲142天。

4) 印度疫情防控的現狀分析及建議。由上述模擬分析可知，印度的疫情防控將是一個漫長的過程，且感染人數十分巨大，這與印度目前的實際情況相符。雖然印度從3月25日實施了全國封鎖，采取了較嚴格的管控措施，但是由于印度存在人口基數大且密度大、核酸檢測率不高、醫療資源有限、無癥狀感染者多等問題，導致疫情防護工作很難快速取得成效。此外，從5月1日開始印度開通專列幫助農民工返鄉，其中不少人在抵達目的地后被檢測出新冠肺炎陽性，因此此次返鄉活動提高了病毒傳播的風險。以印度德里為例，在近1萬人的確診病例中，超過75%是無癥狀或輕癥感染者，且目前整體的核酸檢測率未達到0.2%，在這樣的情況下，放開公共交通會有極大的風險。

綜上，從控制疫情的角度考慮，建議印度在感染人數的峰值(2020年10月16日)到來前控制公共交通的恢復比例；出現社區傳染的貧民窟是疫情防控的關鍵；必須繼續執行嚴格的限行措施才能防止疫情更加惡化。此外，防疫消毒、核酸檢測、精準隔離、醫療條件等均是阻止印度成為全球疫情“震中”地區的必要條件。

4 結 語

本文將SIR-F模型與疫情增長態勢的分類、防控措施影響分析方法相結合，以真實疫情數據為基礎，使用Python實現模擬仿真，實現了對疫情的趨勢分析、重點疫情地區的識別、疫情階段劃分、防控措施對疫情的影響分析、疫情預測分析。以意大利為“停滯”增長態勢的例子，以印度為“爆發”增長態勢的例子，對其疫情的發展過程進行了分階段的擬合分析，實驗證明，分階段的擬合方法能比較準確地刻畫疫情數據隨時間的變化規律，并且在邏輯上更符合客觀事實。本文分析方法同樣適用于其他國家或地區，可以為新冠肺炎等類似傳染病疫情防控提供較大的參考和指導價值。

實驗結果證明了意大利的出行限制措施對疫情已經起到了有效作用。印度的確診數量雖然低于美國等多個國家，但其爆發式增長狀況極其嚴峻，預測得到其感染人數最高值為3.84億人；若印度增強防控力度，使得感染系數降低36.59%，則其感染人數最高值將降低47.66%。最后，結合印度的最新情況，本文建議印度應重點加強貧民窟的疫情防控力度，在2020年10月16日之前繼續控制公共交通，并結合其他精準防控手段阻止其疫情繼續惡化。