姚 藝,錢亞冠,陳曉霞,徐慶華,方科彬

(浙江科技學(xué)院 理學(xué)院,杭州 310023)

新型冠狀病毒肺炎(以下簡稱新冠肺炎)發(fā)生以來,對疫情未來趨勢的預(yù)測一直是學(xué)術(shù)界關(guān)注的熱點(diǎn)。針對此次疫情不少研究者對其展開探討。Zhao等[1]基于2020年1月10日至24日的公開數(shù)據(jù)進(jìn)行指數(shù)增長趨勢曲線擬合,由此判斷新型冠狀病毒早期傳播能力接近或略高于SARS。倪順江[2]把復(fù)雜網(wǎng)絡(luò)理論和流行病學(xué)研究相結(jié)合,建立了傳染病動力學(xué)模型,試驗(yàn)發(fā)現(xiàn)傳染病感染人數(shù)初期時(shí)以指數(shù)形式增長。這些研究都認(rèn)為新冠肺炎初期傳播具有指數(shù)增長的特點(diǎn),對判斷新冠肺炎初期傳播的規(guī)模和特點(diǎn)具有重要意義。但傳統(tǒng)的曲線擬合完全是基于數(shù)據(jù)趨勢來進(jìn)行預(yù)測,無法考慮傳染病的傳播速度、傳播模式及各種防控措施的實(shí)施等動態(tài)信息,預(yù)測效果在初期過后并不可靠。張琳[3]利用一般增長模型,分3個(gè)階段擬合新冠肺炎確診人數(shù)(初期,無障礙指數(shù)增長階段;中期,次指數(shù)增長階段;后期,次線性增長階段)。而次指數(shù)增長方式的內(nèi)在機(jī)制,主要有空間異質(zhì)性、群聚感染及感染參數(shù)的時(shí)間異質(zhì)性[4]等幾方面解釋。在SIR(susceptible-infective-removal)模型中,考慮空間異質(zhì)性[5]、傳染率的時(shí)變性[6]及異質(zhì)混合模型[7],都能夠呈現(xiàn)出次指數(shù)的增長方式。Yang等[8]運(yùn)用長短期記憶(long-short-term-memory,LSTM)模型預(yù)測中國疫情將在2月底達(dá)到高峰,并通過機(jī)器學(xué)習(xí)算法展示若取消湖北省的交通封閉措施,將導(dǎo)致湖北省在3月中旬出現(xiàn)第二次高峰。

2020年3月17日,意大利新冠肺炎死亡人數(shù)超過武漢[9],成為疫情嚴(yán)重國家。這引起了研究者們的關(guān)注。陳曉平[10]通過聚類分析方法研究各國新冠肺炎疫情發(fā)展的嚴(yán)重程度,將國外新冠肺炎疫情國家分為主要疫區(qū)、重點(diǎn)疫區(qū)、分散疫區(qū)3類,其中意大利屬于主要疫區(qū)。劉波[11]認(rèn)為,意大利新冠肺炎疫情是全球化背景下的集體問題,這場疫情既考驗(yàn)歐盟處理公共突發(fā)事件的能力也給歐洲乃至世界產(chǎn)生了深刻而復(fù)雜的影響。

現(xiàn)有的新冠肺炎疫情預(yù)測方法包括曲線擬合(curve fitting)、傳染病動力學(xué)模型(epidemic dynamics model)及人工智能算法(artificial-intelligence,AI)三大類[12]。疫情預(yù)測的主流算法有改進(jìn)的傳染病動力學(xué)模型和LSTM模型[13]。為了更好地預(yù)測意大利新冠肺炎各階段的傳播情況,本研究選用logistic(邏輯回歸)模型、SIR模型和LSTM模型進(jìn)行預(yù)測。通過對試驗(yàn)結(jié)果的分析和對比,提取各模型中有價(jià)值的預(yù)測信息,取長補(bǔ)短,為此次疫情提供更為全面細(xì)致的預(yù)測結(jié)果。

1 預(yù)測模型

1.1 logistic模型

logistic[14]模型是Verhulst-Pear在修正密度方程時(shí)提出的,常用于描述種群、傳染病增長及商品銷售量預(yù)測等方面。logistic函數(shù)如下:

(1)

式(1)中:P(t)為t時(shí)刻的人口數(shù);P0為初始容量;r為人口增長率;K為環(huán)境最大容納量。

使用logistic模型預(yù)測的主要步驟如下:

1) 利用最小二乘法構(gòu)造損失函數(shù):

(2)

式(2):Pt為t時(shí)刻實(shí)際確診人數(shù);P(t)為t時(shí)刻模型預(yù)測的確診人數(shù)。

2) 將損失函數(shù)分別對模型參數(shù)K、P0、r求導(dǎo)建立方程組,即

求出與真實(shí)確診病例數(shù)據(jù)最接近的logistic函數(shù):

由上述logistic函數(shù)定義可知此次疫情累計(jì)確診病例人數(shù)最大會達(dá)到268 999人;模型的初始容量P0為22 060人,人口增長率r為0.01。

3) 利用P(t)函數(shù)預(yù)測在測試集上的確診病例數(shù),用平均絕對百分比誤差[15](mean absolute percentage error,MAPE)指標(biāo)計(jì)算模型精確度。

(3)

式(3)中:Yt為預(yù)測值;yt為真實(shí)值;n為預(yù)測點(diǎn)的個(gè)數(shù)。

MAPE值大小與預(yù)測精度呈反比關(guān)系,MAPE值越小,預(yù)測越準(zhǔn)確。一般情況下,如果MAPE值小于20%,則模型擬合效果好。

由于意大利于2020年5月4日實(shí)施解封,導(dǎo)致實(shí)際累計(jì)病例數(shù)在后期仍持續(xù)向上增長,而在logistic模型中,人口數(shù)在后期達(dá)到飽和狀態(tài)后會維持水平無法向上增長。因此,本研究把數(shù)據(jù)分成2020年2月21日至5月3日和5月4日至8月18日分別進(jìn)行擬合,并在兩段輸入數(shù)據(jù)里劃分訓(xùn)練集和測試集(比例為7∶3)。

1.2 SIR模型

SIR模型是一種倉室模型[16],它將發(fā)病地區(qū)的所有人分為易感者(Susceptible,S)、感染者(Infective,I)及移出者(Recovered,R)3類,并通過某一群體轉(zhuǎn)移至另一群體的傳播學(xué)機(jī)制來建立常微分方程組,以揭示疫情傳播的規(guī)律。

模型設(shè)置了如下假定:1)冠狀病毒通過飛沫直接傳播;2)總?cè)丝跀?shù)是一個(gè)常數(shù),不發(fā)生變化,任何時(shí)刻的3類人群總數(shù)不變;3)易感者為尚未感染的健康人群;4)感染者是已被感染、具有傳染性的人群;5)感染者一段時(shí)間之后可以康復(fù),且成為移出者后對于該疾病具有免疫能力。

SIR模型的常微分方程組形式[17]如下:

(4)

式(4)中:β為易感者被傳染的平均概率;γ為感染者痊愈的平均速率;t為冠狀病毒傳播天數(shù);N為人口總數(shù)。

(5)

存在。式(4)的第2個(gè)方程意味著I(t)的變化依賴S(t)的大小。如果

(6)

那么當(dāng)t→∞時(shí)可以得到I0>I(t)→0,此時(shí)疾病將消亡而不會暴發(fā)。根據(jù)基本再生數(shù)理論,定義:

(7)

基本再生數(shù)R0[18]是刻畫傳染病是否會暴發(fā)的閾值。當(dāng)R0<1時(shí),疾病不會暴發(fā),且隨時(shí)間的推移逐漸消亡;當(dāng)R0>1時(shí),疾病會暴發(fā)一段時(shí)期,在達(dá)到I(t)最高值后隨時(shí)間的推移逐漸消亡。同時(shí),我們定義相應(yīng)的當(dāng)日再生數(shù)為

(8)

式(8)中:S(t)為當(dāng)日易感者的人群數(shù)量;βt為當(dāng)日的平均傳染率;γt為當(dāng)日的平均恢復(fù)系數(shù)。

使用SIR模型預(yù)測的主要步驟如下:

1) 將輸入數(shù)據(jù)進(jìn)行規(guī)范化處理,輸入數(shù)據(jù)的區(qū)間為2020年2月21日至8月18日。

2) 利用最小二乘法構(gòu)造損失函數(shù):

(9)

式(9)中:I(t)為真實(shí)感染者病例數(shù);R(t)為真實(shí)抵抗者病例數(shù);It為模型感染者病例數(shù);Rt為模型抵抗者病例數(shù)。為了加強(qiáng)模型的擬合效果,損失函數(shù)中對確診病例數(shù)的誤差項(xiàng)進(jìn)行了加權(quán),取權(quán)重a=9。

3) 將損失函數(shù)分別對模型參數(shù)β、γ建立方程組,即

以求出一個(gè)既滿足SIR常微分方程關(guān)系又與真實(shí)確診病例數(shù)據(jù)最接近的SIR曲線,最后求得的SIR模型常微分關(guān)系如下:

此時(shí),易感者的平均被傳染率βt為0.000 001 29;感染者痊愈的平均速率γt為0.031 408 08。

4) 利用求得的SIR模型常微分關(guān)系預(yù)測確診病例數(shù)據(jù),用MAPE指標(biāo)驗(yàn)證模型的準(zhǔn)確性。

1.3 LSTM模型

LSTM模型[19]能夠?qū)﹂L時(shí)間跨度的時(shí)間序列保持良好的記憶,該模型在時(shí)間序列的預(yù)測問題上有著突出的表現(xiàn)。LSTM模型原理是充分利用過去若干天數(shù)的數(shù)據(jù)信息預(yù)測未來若干天數(shù)的數(shù)據(jù)。使用LSTM模型進(jìn)行確診病例數(shù)預(yù)測的主要步驟如下:

1) 讀取疫情數(shù)據(jù),將輸入數(shù)據(jù)進(jìn)行規(guī)范化處理并劃分為訓(xùn)練集和測試集,劃分時(shí)序窗口,設(shè)置窗口長度L=3,學(xué)習(xí)率α=0.005,迭代偏置為0.000 01,最大迭代次數(shù)T=2 000;

2) 將輸入數(shù)據(jù)經(jīng)過輸入層的所有節(jié)點(diǎn)傳遞到隱藏層,隱藏層結(jié)合Sigmiod激活函數(shù)將變量映射到[0,1]之間以決定信息是否向下層傳遞;

3) 運(yùn)用tanh函數(shù)計(jì)算當(dāng)前狀態(tài)xt、上一刻LSTM的隱藏層狀態(tài)ht-1和記憶單元狀態(tài)ct-1等參數(shù)值經(jīng)過遺忘門疊加處理形成的新記憶單元狀態(tài)值ct;

4) 利用Sigmoid函數(shù)獲得初始輸出,并利用tanh函數(shù)對數(shù)據(jù)縮放得到預(yù)測輸出ht;

6) 權(quán)重和偏置更新,基于損失函數(shù)E使用Adam梯度優(yōu)化器更新W和b,當(dāng)更新前后差值的絕對值小于ε時(shí)結(jié)束循環(huán);

7) 輸入測試數(shù)據(jù)集,驗(yàn)證模型準(zhǔn)確性;

8) 運(yùn)用訓(xùn)練后的LSTM模型對未來確診人數(shù)進(jìn)行預(yù)測,即用歷史確診人數(shù)xt,xt+1,…,xt+L-1預(yù)測下一時(shí)刻的確診人數(shù)pt+L。

2 試驗(yàn)結(jié)果分析

2.1 試驗(yàn)數(shù)據(jù)

試驗(yàn)數(shù)據(jù)來源于意大利民事保護(hù)部門的疫情通報(bào)(http://www.salute.gov.it/portale/home.html),時(shí)間范圍為2020年2月21日至8月18日,包括現(xiàn)存確診人數(shù)、累計(jì)確診人數(shù)、死亡人數(shù)和治愈人數(shù)4個(gè)指標(biāo),數(shù)據(jù)質(zhì)量良好,無缺失值。其中,logistic模型使用的數(shù)據(jù)為2020年2月21日至8月18日的累計(jì)確診人數(shù);SIR模型使用的數(shù)據(jù)為2020年2月21日至8月18日的現(xiàn)存確診人數(shù)和治愈人數(shù);LSTM模型使用的數(shù)據(jù)為2020年2月21日至8月18日的累計(jì)確診人數(shù)和現(xiàn)存確診人數(shù)。

2.2 模型精度比較

為了避免試驗(yàn)的偶然性,我們進(jìn)行了重復(fù)試驗(yàn),取均值作為該模型的MAPE指標(biāo)值,訓(xùn)練結(jié)果見表1。

表1 MAPE值比較Table 1 Comparison of MAPE values

從試驗(yàn)結(jié)果中可以看到3種模型的MAPE值均遠(yuǎn)小于20%,即3種模型預(yù)測結(jié)果均有效。其中,LSTM模型的MAPE值最小,為0.367 5%。3種模型得到的MAPE指標(biāo)值均非常小,其原因可能是分母實(shí)際確診病例的數(shù)值過大。模型SIR的MAPE指標(biāo)值與其他兩種模型相差20倍之多。導(dǎo)致這一差異的原因可能是SIR模型無法將突發(fā)的“解封”措施考慮在內(nèi),使得預(yù)測效果在后期并不理想。而Logistic模型則通過分段預(yù)測的方法大大提高了其預(yù)測精度。

通過試驗(yàn)發(fā)現(xiàn),Logistic、SIR和LSTM模型在預(yù)測中均各有優(yōu)勢與不足。因此,研究者需要綜合考慮,提取模型各自最有價(jià)值的預(yù)測信息以得到更為全面的預(yù)測結(jié)果,3種模型的優(yōu)缺點(diǎn)對比見表2。

表2 3種模型的優(yōu)缺點(diǎn)對比Table 2 Comparison of advantages and disadvantages of three models

2.3 緩增期的模型預(yù)測

由logistic模型預(yù)測得到的確診病例數(shù)結(jié)果如圖1所示。紅色的虛線為模型感染病例數(shù)的預(yù)測值,模型預(yù)測區(qū)間為2020年2月21日至12月15日,藍(lán)色的實(shí)線為截至2020年8月18日官方公布的實(shí)際數(shù)據(jù)。

圖1 logistic模型累計(jì)確診病例趨勢預(yù)測Fig.1 Trend prediction of cumulative confirmed cases by Logistic model

根據(jù)上文試驗(yàn)分析模型預(yù)測在2020年9月29日后疫情進(jìn)入緩增期,隨后確診病例數(shù)繼續(xù)向上增長,最終的累計(jì)確診人數(shù)K會達(dá)到268 999人。從圖1中可以看出,模型的預(yù)測曲線與實(shí)際數(shù)據(jù)的總體增長趨勢一致,在疫情后期兩者的偏差隨著時(shí)間節(jié)點(diǎn)的推移而逐漸變大。

2.4 當(dāng)日再生數(shù)的模型預(yù)測

由SIR模型預(yù)測得到的確診病例數(shù)結(jié)果如圖2所示。在SIR模型中用I類人群擬合實(shí)際確診病例數(shù)。I(t)、R(t)為真實(shí)的確診人數(shù)和移出者人數(shù);St、It、Rt分別為模型預(yù)測的易感者人數(shù)、感染者人數(shù)和移出者人數(shù)。從圖2中可以看出,模型預(yù)測在2020年4月21日出現(xiàn)拐點(diǎn),現(xiàn)存感染人數(shù)的峰值可達(dá)到約124 554人(略高于真實(shí)數(shù)據(jù))。但在2020年7月10日后,模型預(yù)測曲線開始偏離真實(shí)數(shù)據(jù),可見,在疫情后期SIR模型的預(yù)測效果不佳。

圖2 SIR模型現(xiàn)存確診病例趨勢預(yù)測Fig.2 Trend prediction of existing confirmedcases by SIR model

以2020年2月11日至3月8日為第一個(gè)時(shí)間段,每增加一天構(gòu)成一個(gè)新的時(shí)間段,用SIR模型求解其當(dāng)日再生數(shù)Rt,得到了當(dāng)日感染再生數(shù)隨時(shí)間推移的演化趨勢圖,預(yù)測區(qū)間為2020年3月8日至7月10日,如圖3所示。

圖3 當(dāng)日再生數(shù)隨時(shí)間的演化情況Fig.3 Change of intraday regenerationnumber over time

從圖3中可以發(fā)現(xiàn),在2020年3月8日意大利全國“封城”之后當(dāng)日基本再生數(shù)上下不穩(wěn)定地波動,總體向上增長。在2020年4月4日時(shí)當(dāng)日再生數(shù)達(dá)到峰值,之后開始快速下降,隔離措施產(chǎn)生了作用,且當(dāng)日再生數(shù)在2020年5月1日至7月10日一直小于1,疫情得到有效控制。

2.5 隔離措施對疫情傳播的影響

從2020年3月8日起,意大利正式簽署“封城令”,開始實(shí)施歐洲史上規(guī)模最大的封城措施。為驗(yàn)證物理隔離對阻斷新冠肺炎傳播是否有效果。本研究通過提高傳染率來模擬意大利不采取封城措施的情形,并與SIR模型預(yù)測到的感染者人數(shù)放在一起做對比分析,結(jié)果如圖4所示。

圖4 防控隔離措施對疫情的影響Fig.4 Impact of physical isolation measureson epidemic

圖4中,I(t)為官方提供的確診人數(shù);It為模型預(yù)測的感染者人數(shù)。SIR模型參數(shù)β為易感者被感染的平均概率。采取隔離措施相當(dāng)于降低傳染率β,兩者高度相關(guān)。將模型中的參數(shù)β提高到原來的2倍來模擬沒有采取“封城”措施時(shí)感染者病例數(shù)隨時(shí)間變化的曲線It(2β)。

通過分析我們發(fā)現(xiàn),目前的防控措施對疫情起到了很好的抑制作用。如果不采取相應(yīng)隔離措施,現(xiàn)存感染人數(shù)會比目前高出30 000人左右,也就是說,如果沒有封鎖措施,將導(dǎo)致感染人數(shù)峰值更高且更快達(dá)到峰值,造成更大的生命財(cái)產(chǎn)損失以及社會恐慌。

2.6 集中收治對感染人數(shù)峰值的影響

在新冠肺炎傳播的高峰期,醫(yī)療人員和物資的供應(yīng)是一個(gè)關(guān)鍵問題。醫(yī)院的收治效率與SIR模型中的恢復(fù)系數(shù)γ呈正相關(guān)。本研究通過降低恢復(fù)系數(shù)來模擬醫(yī)院收治效率不高的情況,結(jié)果如圖5所示。I(t)為官方提供的確診病例數(shù);It為模型預(yù)測的感染者病例數(shù);It(γ/2)為模型中的恢復(fù)系數(shù)γ減少到原來一半時(shí)的感染者病例數(shù)。

圖5 醫(yī)院收治效率對疫情的影響Fig.5 Impact of hospital admission efficiencyon epidemic

從圖5中可以看出,集中收治在疫情高峰期可以避免約3萬人被感染,在疫情后期現(xiàn)存病例數(shù)也明顯低于未采取集中收治的狀況。綜上,集中收治對感染者病例數(shù)的迅速回落發(fā)揮了關(guān)鍵作用。

2.7 意大利疫情后期的增長趨勢

由LSTM模型預(yù)測得到的確診病例數(shù)結(jié)果如圖6～7所示,模型預(yù)測的區(qū)間為2020年2月11日至8月24日。紅色的虛線為模型在訓(xùn)練集上對累計(jì)和現(xiàn)存感染人數(shù)的預(yù)測值,綠色的虛線為模型在測試集上對累計(jì)和現(xiàn)存感染人數(shù)的預(yù)測值,藍(lán)色的實(shí)線為截至2020年8月18日官方公布的數(shù)據(jù)。

圖6 LSTM模型累計(jì)確診病例趨勢預(yù)測Fig.6 Trend prediction of cumulative confirmed cases by LSTM model

從圖7中可以看出,在2020年4月21日現(xiàn)存確診病例數(shù)達(dá)到峰值后開始迅速下降,隔離措施初見成效。但疫情好轉(zhuǎn)的同時(shí)也給意大利經(jīng)濟(jì)帶來重創(chuàng),為此意大利在2020年5月4日被迫實(shí)行解封措施。從圖7可以看出,確診病例數(shù)在達(dá)到低點(diǎn)后又開始繼續(xù)上升且增長速率越來越快,進(jìn)入二次增長階段。如果意大利政府在這種情況下放任不管、不采取任何防控措施,勢必會導(dǎo)致疫情二次爆發(fā)。

圖7 LSTM模型現(xiàn)存確診病例趨勢預(yù)測Fig.7 Trend prediction of existing confirmed cases by LSTM model

3 結(jié) 論

本研究根據(jù)意大利“封城”和“解封”兩個(gè)重要時(shí)間節(jié)點(diǎn),將疫情傳播劃分為初期(指數(shù)上升階段)、中期(次指數(shù)增長階段)和后期(二次增長階段)3階段,并利用logistic回歸模型、SIR模型和LSTM模型,對意大利疫情中后期傳播過程中感染者的發(fā)展趨勢進(jìn)行建模和分析。通過試驗(yàn)驗(yàn)證,擬合結(jié)果與實(shí)際數(shù)據(jù)高度吻合,從而可為疫情趨勢分析提供有效的數(shù)據(jù)分析支持。由此得出如下結(jié)論:

1) 意大利的緩增期出現(xiàn)于2020年9月29日左右,預(yù)計(jì)最終累計(jì)確診病例數(shù)將達(dá)到26.89萬人。當(dāng)日再生數(shù)Rt在不斷減小,這與政府的防控措施及民眾的居家防護(hù)有直接關(guān)系。

2) 根據(jù)模型預(yù)測結(jié)果,在疫情后期意大利確診人數(shù)會進(jìn)入二次增長階段,建議意大利推遲全部企業(yè)復(fù)工復(fù)產(chǎn)及學(xué)校開學(xué)的時(shí)間。目前,意大利疫情控制還不穩(wěn)定,甚至海外一些國家不僅僅是面對新冠疫情,例如“西尼羅病毒”和“沙門氏菌病毒”席卷美國,新加坡受到“登革熱”疫情的困擾。為了防止這些傳染病流入中國,我們還需加強(qiáng)出入境管理。