凡友榮 楊 濤 孔華鋒
1(公安部第三研究所 上海 201204) 2(武漢商學院 湖北 武漢 430056)
根據(jù)世界衛(wèi)生組織公布的數(shù)據(jù),截至2020年5月27日,全球新冠肺炎確診病例已達548萬,超過10萬確診人數(shù)的國家包括美國、巴西、俄羅斯、英國、意大利、法國、德國、土耳其、印度、伊朗、秘魯。其中美國、巴西目前的每日新增確診數(shù)超過1萬人;意大利、法國等國家的疫情趨緩。由于印度存在人口密集、醫(yī)療條件差等問題,世界衛(wèi)生組織表示,人類對抗新冠疫情能否取得決定性勝利,未來很大程度上取決于印度控制該病毒的能力。針對全球各國不同的疫情狀況,分別進行現(xiàn)狀分析及防控措施的影響研究至關(guān)重要。
SIR模型、SEIR(Susceptible-Exposed-Infected-Recovered)模型是傳染病動力學研究的經(jīng)典模型,目前已有多項基于SIR、SEIR模型的新冠疫情研究。尹楠[1]基于SIR模型進行了有限區(qū)域內(nèi)新冠病毒傳播的模擬分析。梅文娟等[2]將極限學習機與SIR模型相結(jié)合構(gòu)建了極限IR實時預測模型。盛華雄等[3]將SIR模型與Logistic模型相結(jié)合,對不同防控措施進行了模擬和對比。喻孜等[4]使用基于時變參數(shù)的 SIR 模型預測得到了疫情拐點和最大確診數(shù)。Zhou等[5]基于SEIR模型計算得到了COVID-19的再生數(shù)。林俊鋒[6]在SEIR模型的基礎上增加了“隱形傳播者”,對隱形傳播者的數(shù)量進行了預測。范如國等[7]基于SEIR模型對3種不同潛伏期下的疫情進行了拐點預測。
本文使用的SIR-F模型在SIR模型的基礎上增加了未確診的感染者(具有傳染性),這一點與SEIR模型相同;本文模型還將“移除者”(Removed)細分為“治愈者”(Recovered)和“死亡者”(Fatal),從而能夠進一步分析治愈人數(shù)、死亡人數(shù)的變化趨勢。此外,考慮到在數(shù)據(jù)擬合的實踐中存在過擬合的現(xiàn)象,并且實際的疫情發(fā)展過程分為多個階段,每個階段的趨勢特征各不相同,因此本文先進行階段劃分,再將疫情數(shù)據(jù)分階段地進行擬合分析,從而提高分析結(jié)果的準確性。
本文的技術(shù)路線如圖1所示。首先計算各國新冠肺炎確診數(shù)量的增長因子,根據(jù)該因子將所有國家分為3個增長類型:“停滯”式、“爆發(fā)”式、“波動”式;然后在“停滯”增長態(tài)勢的國家中選擇典型例子進行案例分析,由于意大利是全球重點疫情國家之一,且其防控工作已取得顯著成效,因此選擇意大利為例,對其疫情發(fā)展過程和防控措施進行深入分析;最后針對“爆發(fā)”式增長的國家,從中識別目前爆發(fā)態(tài)勢最嚴峻的國家,并針對其進行疫情的深入分析,提出防控措施建議。
在案例分析中,首先使用Python調(diào)用Prophet的“突變點自動監(jiān)測”功能對疫情數(shù)據(jù)進行階段劃分,可用于疫情數(shù)據(jù)突變點的自動識別?;赟IR-F模型進行分階段的數(shù)據(jù)擬合、參數(shù)估計,對比不同階段的參數(shù),從而發(fā)現(xiàn)各階段疫情的變化情況,然后結(jié)合防控措施的實施時間,研究防控措施對疫情發(fā)展產(chǎn)生的影響。對于“爆發(fā)”式增長的國家,在模型擬合的基礎上進行疫情預測,并且基于不同的感染系數(shù)進行預測比對,從而模擬不同強度防控措施下的疫情發(fā)展趨勢。
本文通過增長因子來評估 “新冠肺炎”確診人數(shù)增長的態(tài)勢,增長因子的計算公式如下:
(1)
式中:GF為增長因子;ΔCn為第n天確診人數(shù)增加的數(shù)量。計算每個國家近7天的GF,進行增長情況的分類:

(2)
GC為增長態(tài)勢所屬類別,GC=1代表“爆發(fā)”的增長態(tài)勢,GC=-1代表“停滯”的增長態(tài)勢,GC=0代表“波動”的增長態(tài)勢。
由于COVID-19疫情比較特殊,有一部分感染者在死亡前仍未正式確診,將這類情況增加到經(jīng)典的SIR模型中,構(gòu)建SIR-F模型。SIR-F模型中從未感染者轉(zhuǎn)變?yōu)楦腥菊?、治愈者、死亡者的過程如圖 2 所示。

圖2 SIR-F模型中感染者的轉(zhuǎn)化過程
圖2中,S代表未感染者,S*代表未確診的感染者,I代表感染者,R代表治愈者,F(xiàn)代表死亡者,α1為S*的死亡系數(shù),α2為I的死亡系數(shù),β為感染系數(shù),γ為I的治愈系數(shù)。累計確診數(shù)量C=I+R+F。
根據(jù)上述模型設計,構(gòu)建其動力學方程如下:
(3)
(4)

(5)
(6)
式中:N為總的人口數(shù)量,N=S+I+R+F;T為疫情開始的總時長。為求解方便,將上述函數(shù)方程進行無量綱處理,設(S,I,R,F)=N×(x,y,z,w),(T,α1,α2,β,γ)=(τt,θ,τ-1k,τ-1ρ,τ-1σ),得到以下微分方程組:
(7)
在Python(版本3.7.6)環(huán)境下,使用超參數(shù)自動優(yōu)化框架Optuna對上述微分方程進行求解,從而得到各個參數(shù)(θ,k,ρ,σ)的估計值。
各國政府通過疫情防控措施改變民眾的行為,從而控制疫情的發(fā)展。本文引入感染系數(shù)β受民眾行為影響的函數(shù),如式(8)所示。設感染系數(shù)β的影響因子包括:每分鐘內(nèi)未感染者與感染者在外出過程中接觸的次數(shù)c,每周內(nèi)未感染者外出的天數(shù)g,唾液中病毒存在的可能性v,有效佩戴口罩的人數(shù)比例m,洗手降低的感染概率we,其他影響因素δ。
(8)
式(8)與SIR-F模型相結(jié)合,能夠模擬居民的外出、佩戴口罩、洗手等行為對疫情發(fā)展的影響。
本文基于的全球疫情數(shù)據(jù)來源于約翰·霍普金斯大學的系統(tǒng)科學與工程中心(https://systems.jhu.edu/ ),包含了全球217個國家在110天內(nèi)(2020年1月22日至2020年5月11日)的累計確診人數(shù)、現(xiàn)存確診人數(shù)、死亡人數(shù)、治愈人數(shù),總數(shù)據(jù)量為23 804行。各國防控措施數(shù)據(jù)來源于牛津大學的流行病預測研究組織(http://epidemicforecasting.org/),包含了96個國家在2019年12月18日至2020年4月15日期間采取的管控措施,總數(shù)據(jù)量為1 633行。
將全球疫情數(shù)據(jù)導入式(1)和式(2),計算得到96個國家的疫情增長態(tài)勢所屬類別GC,結(jié)果示例如表1所示。

表1 全球各國疫情增長態(tài)勢的分類結(jié)果
處于“爆發(fā)”增長態(tài)勢的包括印度、俄國、愛沙尼亞、巴西、肯尼亞等53個國家。將這類國家的疫情數(shù)據(jù)聚合,增長趨勢如圖3所示。

圖3 “爆發(fā)”式增長國家的疫情數(shù)據(jù)變化趨勢
處于“停滯”增長態(tài)勢的包括意大利、日本、利比亞、安哥拉、尼加拉瓜、越南等28個國家。將這類國家的確診數(shù)據(jù)聚合,增長趨勢如圖4所示。

圖4 “停滯”式增長國家的疫情數(shù)據(jù)變化趨勢
處于“波動”增長態(tài)勢的包括加拿大、以色列、美國、英國等138個國家。將這類國家的確診數(shù)據(jù)聚合,增長趨勢如圖5所示。

圖5 “波動”式增長國家的疫情數(shù)據(jù)變化趨勢
觀察圖3至圖5可知,“爆發(fā)”式呈現(xiàn)迅速上升的趨勢;“停滯”式的增長態(tài)勢已基本得到了控制;“波動”式的未來發(fā)展方向尚不明確?!氨l(fā)”式和“波動”式增長的國家是全球疫情控制的重點,因此有必要研究“停滯”式增長的國家所采取的防控措施,為其他國家的疫情防控提供有效建議。
1) 意大利疫情趨勢分析。意大利目前的疫情態(tài)勢如圖6所示。根據(jù)式(1)計算得到意大利確診人數(shù)的增長因子變化態(tài)勢如圖7所示,其增長因子GF<1的持續(xù)天數(shù)為18天。

圖6 意大利疫情數(shù)據(jù)趨勢圖

圖7 意大利增長因子變化趨勢圖
使用Prophet 對意大利的累計確診數(shù)據(jù)進行突變點的識別,如圖8所示。

圖8 意大利累計確診數(shù)據(jù)的突變點識別
由圖8可知,存在3個突變點:“2020-03-17”、“2020-04-05”、“2020-04-25”,則意大利的疫情發(fā)展過程可以劃分為4個階段,如表2所示。每個階段的變化規(guī)律有明顯的變化,因此下文將分階段地進行擬合分析。

表2 意大利疫情發(fā)展的4個階段
2) SIR-F模型擬合及參數(shù)估計。將4個階段的真實數(shù)據(jù)分別導入SIR-F模型進行擬合,擬合效果如圖9所示。

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

表3 意大利4個階段的參數(shù)估計值
3) 感染率降低的原因分析。意大利在3月采取了多項限制居民出行的管控措施,如表4所示。

表4 意大利3月采取的限制出行措施
考慮到政府的管控措施從開始執(zhí)行到產(chǎn)生實際的疫情控制效果存在時間延遲,因此認為第1階段為實施管控措施之前,第2、第3、第4階段為實施管控措施之后。因此本文假設居民在第1階段的出行情況與平時相同,設意大利民眾平時每周出行5.6天,即g1=5.6。
根據(jù)2020年3月13日報道,由于出行限制、公共活動場所的關(guān)閉,使得居民相互之間的潛在接觸次數(shù)降低了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??梢娒癖娒恐艿某鲂刑鞌?shù)已經(jīng)下降到了0.51天。
上述分析過程從實際的疫情數(shù)據(jù)反推得到了實施管控后的每周出行天數(shù),可知民眾的每周出行天數(shù)減少了90.9%,使得感染系數(shù)降低了92.6%,證明了限制出行的管控措施對新冠病毒的擴散起到了明顯的抑制作用;同時證明了本文的分析模型能夠較靈敏地模擬現(xiàn)實中防疫管控下新冠肺炎的傳播情況。
除了上述管控措施,意大利的疫情防控和診療工作均參照了中國方案,這也是意大利疫情控制效果較好的重要原因。意大利從5月4日已進入抗疫和恢復經(jīng)濟的并行階段,新聞顯示62%的意大利人對于復工感到害怕,因此出行習慣的恢復還需要時間。若出行量逐漸增加,應從其他方面控制感染率,例如提高戴口罩、勤洗手習慣的普及率等,從而防止疫情再次反彈。
印度的增長因子GF>1的持續(xù)天數(shù)為全球最長(66天),其增長態(tài)勢極其嚴峻,因此下文將印度作為“爆發(fā)”式增長的典型例子進行深入分析。
1) 印度疫情趨勢分析。印度目前的疫情態(tài)勢如圖10、圖11所示。

圖10 印度疫情數(shù)據(jù)趨勢圖

圖11 印度增長因子變化趨勢圖
對印度的累計確診數(shù)據(jù)進行突變點的識別,如圖12所示,存在2個突變點(“2020-04-03”、“2020-04-28”),則印度的疫情發(fā)展過程可以劃分為3個階段,如表5所示。

圖12 印度累計確診數(shù)據(jù)的突變點識別

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

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

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

圖14 假設ρ值降至0.026的疫情預測
通過上述對比分析可知,若印度增強防控措施的力度,將感染系數(shù)在目前的基礎上降低36.59%,感染人數(shù)的最高值即可降低47.66%,到達最高感染人數(shù)的日期將推遲142天。
4) 印度疫情防控的現(xiàn)狀分析及建議。由上述模擬分析可知,印度的疫情防控將是一個漫長的過程,且感染人數(shù)十分巨大,這與印度目前的實際情況相符。雖然印度從3月25日實施了全國封鎖,采取了較嚴格的管控措施,但是由于印度存在人口基數(shù)大且密度大、核酸檢測率不高、醫(yī)療資源有限、無癥狀感染者多等問題,導致疫情防護工作很難快速取得成效。此外,從5月1日開始印度開通專列幫助農(nóng)民工返鄉(xiāng),其中不少人在抵達目的地后被檢測出新冠肺炎陽性,因此此次返鄉(xiāng)活動提高了病毒傳播的風險。以印度德里為例,在近1萬人的確診病例中,超過75%是無癥狀或輕癥感染者,且目前整體的核酸檢測率未達到0.2%,在這樣的情況下,放開公共交通會有極大的風險。
綜上,從控制疫情的角度考慮,建議印度在感染人數(shù)的峰值(2020年10月16日)到來前控制公共交通的恢復比例;出現(xiàn)社區(qū)傳染的貧民窟是疫情防控的關(guān)鍵;必須繼續(xù)執(zhí)行嚴格的限行措施才能防止疫情更加惡化。此外,防疫消毒、核酸檢測、精準隔離、醫(yī)療條件等均是阻止印度成為全球疫情“震中”地區(qū)的必要條件。
本文將SIR-F模型與疫情增長態(tài)勢的分類、防控措施影響分析方法相結(jié)合,以真實疫情數(shù)據(jù)為基礎,使用Python實現(xiàn)模擬仿真,實現(xiàn)了對疫情的趨勢分析、重點疫情地區(qū)的識別、疫情階段劃分、防控措施對疫情的影響分析、疫情預測分析。以意大利為“停滯”增長態(tài)勢的例子,以印度為“爆發(fā)”增長態(tài)勢的例子,對其疫情的發(fā)展過程進行了分階段的擬合分析,實驗證明,分階段的擬合方法能比較準確地刻畫疫情數(shù)據(jù)隨時間的變化規(guī)律,并且在邏輯上更符合客觀事實。本文分析方法同樣適用于其他國家或地區(qū),可以為新冠肺炎等類似傳染病疫情防控提供較大的參考和指導價值。
實驗結(jié)果證明了意大利的出行限制措施對疫情已經(jīng)起到了有效作用。印度的確診數(shù)量雖然低于美國等多個國家,但其爆發(fā)式增長狀況極其嚴峻,預測得到其感染人數(shù)最高值為3.84億人;若印度增強防控力度,使得感染系數(shù)降低36.59%,則其感染人數(shù)最高值將降低47.66%。最后,結(jié)合印度的最新情況,本文建議印度應重點加強貧民窟的疫情防控力度,在2020年10月16日之前繼續(xù)控制公共交通,并結(jié)合其他精準防控手段阻止其疫情繼續(xù)惡化。