梅文娟,劉 震*,朱靜怡,杜 立
(1. 電子科技大學自動化工程學院 成都 611731;2. 電子科技大學生命科學與技術學院 成都 611731)
自2020 年1 月以來,新型冠狀病毒疫情迅速蔓延,各地確診病例增加。疫情受到各地政府高度關注,相繼啟動重大突發(fā)衛(wèi)生事件I 級響應并采取不同程度的隔離措施。2020 年1 月30 日,國際衛(wèi)生組織將此次疫情列為“國際公眾衛(wèi)生緊急事件”。在加強醫(yī)療援助及實施隔離措施的同時,全球科研人員針對新冠病毒從流行病學[1]、病理學[2]和藥理學[3]3 個領域進行了分析,為抗擊疫情提供了有力的學術支持。其中,從流行病學角度分析疫情的發(fā)展規(guī)律并預測疫情發(fā)展趨勢對制定和實施合理的干預防控措施具有重要意義。
為評估疫情發(fā)展,現(xiàn)有的方法可分為統(tǒng)計學方法與動力學方法兩種。統(tǒng)計學模型在信息不全的情況下通過某一小樣本的情況對總體進行預測。文獻[4]根據(jù)武漢國際機場每日人流量、該樣本下的確診人數(shù)以及機場往日的流量數(shù)據(jù)對武漢潛在感染人數(shù)進行估計。然而,由于小樣本與總體樣本具有傳播特性差異,該方法的估計結果可能存在較大的偏差。同時,統(tǒng)計學模型無法反映疫情傳播的變化趨勢??紤]人口遷徙因素的影響,文獻[5]根據(jù)武漢遷徙數(shù)據(jù),分析了全國50 個城市感染新型冠狀病毒的確診比率,并利用Bootstrap 方法對確診人數(shù)進行了穩(wěn)健性估計。實驗分析證明了病毒的二代傳播在不同區(qū)域存在一定的差異性,社會活動對疫情發(fā)展存在一定的影響。
與統(tǒng)計學模型不同,動力學模型[6]基于病毒寄生宿主后各個狀態(tài)間的因果關系,利用對該病毒已知的信息與各個狀態(tài)下人群的歷史數(shù)據(jù)構建用于描述傳播機理的微分方程,從而得出對有效再生數(shù)、確診人數(shù)等指標的預測評估值。根據(jù)對疫情發(fā)展中不同群體的劃分情況不同,用于描述傳染病傳播的經(jīng)典數(shù)學模型有SIR 模型和SEIR 模型兩種形式。
SIR 模型[7]假設易感人群(susceptible)、感染人群(infectious)和康復人群(recovered)3 部分人群在病毒發(fā)展規(guī)律下以一定概率向其他狀態(tài)轉移形成“易感態(tài)?感染態(tài)?康復態(tài)”的動力學模型,可實現(xiàn)對病毒的傳染趨勢在一定精度內(nèi)的評估和預測。文獻[8]利用機器學習方法估計SIR 模型,其估計結果基本與實際數(shù)據(jù)吻合。
進一步,SEIR 模型[9]在SIR 模型的基礎上考慮到潛伏期導致感染過程存在遲滯性。因此,SEIR 模型對動力學系統(tǒng)進行細化,串聯(lián)了潛伏期人群(exposed)?;谠撃P?,文獻[10]對病毒基本再生數(shù)進行了初步預測。預測過程證明在SEIR模型下基本再生數(shù)的時間敏感性。此外,文獻[11]基于該疾病的臨床特征、傳染性特征以及政府有關的干預措施在SEIR 模型中增加無癥狀感染人群(pre-symptomatic)和住院人群(hospitalized)兩種狀態(tài),達到細化傳播過程描述的效果,借助蒙特卡洛模擬,認為在當時控制下的基本再生數(shù)高達6.47。文獻[12]就潛伏期長度對疫情的影響做了深入分析,得到在不同的潛伏期假設下,預測趨勢的增長速度,拐點和峰值都存在差別,潛伏期的不同會影響疫情趨勢變化。而從實際的效果來看,潛伏期受到人類宿主生理系統(tǒng)的影響存在個體差異,潛伏期存在的影響呈現(xiàn)出一定的不確定性。
盡管SIR 模型與SEIR 模型對基本再生數(shù)和疫情早期趨勢有較好的估計,然而基于動力學的預測模型在實際應用中仍存在局限性[13]。首先,倉室模型無法對開放式流動環(huán)境下的病毒傳播做出準確估計;第二,該類模型對相關參數(shù)的評估缺乏外界環(huán)境應力影響的引入;第三,對于疾病傳播能力及治愈概率的常數(shù)假設與實際狀況不符。因此,該類模型無法對疫情趨勢做長期準確的分析。
考慮到人為及環(huán)境影響因素,一些學者基于環(huán)境容納量的概念,引入針對疾病傳播的Logistic 增長模型[14]和指數(shù)增長模型[15]。該類模型反映了醫(yī)療和防控對于疫情發(fā)展存在一定時滯性影響[16],并且結合此類因素給出更符合實際的估計。
綜上,本文結合動力學模型對傳播機理描述的準確性和Logistic 增長模型對外界因素估計的有效性,將基于機理描述的微分方程與基于數(shù)據(jù)驅動的極限學習機[17]融合,構建用于實時預測的極限IR 模型。在該模型中,對動力學方程進行改進,突破倉室模型下僅考慮密閉環(huán)境的局限性。同時將傳播最終狀態(tài)進一步分化為治愈態(tài)和死亡態(tài),并預測從感染態(tài)到這兩種狀態(tài)概率的時變規(guī)律,用于進行更加精確的預測。最后,借助極限學習機,生成對確診人數(shù)動態(tài)變化的預測模型,用于疫情趨勢的實時預測和分析。實驗證明,極限IR 算法可實現(xiàn)準確的實時預測。
由實際的防控情況,新型冠狀肺炎的防控過程如圖1 所示。不考慮二次感染,每個感染新型冠狀病毒的個體需要經(jīng)歷易感態(tài)、潛伏態(tài)、感染態(tài)和移除態(tài)4 種狀態(tài)才能排除其對倉室模型中各個人群的變化影響。從實際的防控過程和病毒的病理特征中可以看出,潛伏態(tài)具有以下特點:1) 在實際的防控過程中,潛伏態(tài)群體因為無明顯發(fā)病癥狀難以進行準確的統(tǒng)計;2) 從病理學角度,病毒在感染宿主后已經(jīng)具備傳染力[18],且潛伏期的個體差異性較大(中位值為4 天,最長可達14 天[19]),在估計上的不確定度較大;3) 盡管潛伏期傳染能力遠弱于發(fā)病后的傳染能力,然而由于早期潛伏期人群基數(shù)很大,小概率的傳染仍然會對新增患者群體造成較大的影響。綜上,在病毒傳染機制、環(huán)境中人群個體機制等信息尚不明確的情況下,對于潛伏態(tài)的判定及潛伏態(tài)與感染態(tài)的轉換關系存在分析上的不確定性,不利于定量方法的研究。因此,本文采用SIR 模型作為基礎,簡化了從“易感態(tài)?潛伏態(tài)?感染態(tài)”這一步驟,直接考慮易感態(tài)與感染態(tài)之間的傳遞關系。
利用SIR 模型可以有效地描述病毒的傳播過程。根據(jù)SIR 模型的特性,病毒傳播過程的動力學方程如下:
式中,S 表示為易感態(tài)人群規(guī)模;I 表示為感染態(tài)人群規(guī)模;R 為恢復態(tài)人群規(guī)模;β 指感染人群將易感人群轉換為新增病例的概率;r 表示為現(xiàn)有感染人群轉換為死亡人群或治愈人群的概率。3 種狀 態(tài)的變化特征如圖1 所示。
圖1 新冠肺炎的防控過程
如圖2 所示,鑒于新冠病毒暫時未發(fā)現(xiàn)有明顯變異,康復人群不會發(fā)生第二次感染,因此在SIR模型中認為康復人群不會再次轉變?yōu)橐赘腥巳?,即?0。
圖2 SIR 模型狀態(tài)轉移圖
由于實際的防疫工作中,對于疫情的通報要求以及對于疫情的預測需求均具有一定的時間間隔,因此對動力學方程進行離散化,得到:
式中, ?S(t)、 ?I(t)和 ?R(t)分別為S、I、R 的數(shù)據(jù)差值。根據(jù)實際情況,對于數(shù)據(jù)更新的時間間隔需求為1 天,因此t 取正整數(shù)表示更新的天數(shù)。定義2020 年1 月23 日為第一個時間點,此時t=1。
由于病毒的傳播能力和對于疾病的治療能力受到外界應力的影響,相關指標會隨著疫情發(fā)展變化,本文假設傳播能力β 與概率r 均為關于時間的函數(shù)。由式(5)和式(6),時變函數(shù)β(t)和r(t)的表達式如下:
根據(jù)官方數(shù)據(jù),疫情發(fā)生早期14 天時間間隔下β(t)和r(t)的變化趨勢如圖3 與圖4 所示。在疫情發(fā)生早期,由于疾病的傳播不受人為干預,傳播能力產(chǎn)生波動且由于死亡率偏高引起r(t)較高。隨著防疫力度加大,疾病的傳播能力逐步受到限制。同時,隨著醫(yī)學對疾病防治的效果逐漸顯現(xiàn),死亡率降低,治愈率升高,r(t)呈現(xiàn)出下凹的變化趨勢。
圖3 病毒傳播能力變化趨勢
圖4 概率r(t)變化趨勢
綜上,疫情的變化不僅與病毒的固有特性相關,同時受到外界因素對疫情變化存在一定的時滯性影響。為了更好地分析疫情的變化趨勢,需要構建更為合適的預測模型。在SIR 的基礎上,將恢復態(tài)分化為死亡態(tài)(R1)和痊愈態(tài)(R2)兩組狀態(tài),并假設I 態(tài)分別以r1(t)和r2(t)的概率向R1態(tài)和R2態(tài)轉移。本文借助極限學習機構建關于感染態(tài)的自循環(huán)函數(shù)并對由I 狀態(tài)到R1狀態(tài)及R2狀態(tài)的轉移概率進行預測,提出可用于疫情趨勢實時預測的極限IR 模型,其中各狀態(tài)的傳遞關系如圖5 所示。
圖5 極限IR 模型總體框架
考慮到外界因素的時滯性影響和一定程度的不確定性,采用3 組極限學習機分別對隨時間變化凈新增感染人數(shù)、死亡率和治愈率進行實時預測評估。由圖5 給出的傳遞關系,各個狀態(tài)的動力學方程如下所示:
式中, r1(t)和 r2(t)分別為t 時刻下的治愈率和死亡率;F(I(t?1))為根據(jù)極限學習機構建的自相關函數(shù);N(t)為受到突發(fā)事件影響造成的感染人數(shù)變化,具有不可預測性,因此在構建預測模型的過程中不予考慮。綜上,預測模型的構建轉換為對r1(t)、 r2(t)和F(I(t?1))的預測模型構建問題。
根據(jù)上一節(jié)的分析,本文方法通過極限學習機對疾病的治愈率和死亡率進行預測,其模型結構如圖6 所示。
圖6 極限學習機預測模型框圖
該預測模型為單隱藏層網(wǎng)絡結構。為得到時刻t 較準確的預測估計,利用歷史有效數(shù)據(jù)對極限學習機進行訓練。考慮到數(shù)據(jù)的可靠性,假設有效數(shù)據(jù)起始時間為t0,則訓練數(shù)據(jù)的輸入為:
為確保模型的準確性,采用單步預測的方式,其對應的期望輸出為:
考慮r1(t)受到外界因素的時滯性影響,采用Sigmoid 函數(shù)作為隱藏層的核函數(shù),其表達式如下所示:
同時,考慮到多因素影響,假設概率的變化趨勢為多組logistic 回歸的加權和。由極限學習機的性質,其回歸方程如下:
根據(jù)上式,得到對于歷史數(shù)據(jù)的模型輸出為:
由最小均方準則,計算的到模型的輸出層權重Λ={λ1,λ2,···,λN}如下所示:
在實際的疫情趨勢預測中,疫情的相關數(shù)據(jù)以一天為間隔進行更新。故t+1 輸出層可由上一時刻數(shù)據(jù)更新得出:
式 中, K(t)=(HT(t ?1)H(t ?1))?1; h(r1(t ?1))為t時刻的隱藏層輸入;r ?1(t)為預測模型的輸出。
根據(jù)預測模型得到治愈率的預測值,進而根據(jù)t?1 時刻的感染人數(shù)給出此時的治愈人數(shù)預測值:
同理,對死亡率時間序列建立預測模型,進而得到死亡人數(shù)的預測值:
忽略突發(fā)性因素對模型的影響,本文對凈感染增長人數(shù)定義如下:
由圖5,可以看出 ?I?(t)是僅與I(t)有關的時間函數(shù)。從流行病傳播的角度,? I?(t)的變化趨勢與病毒固有的傳播特性與人為干預有關,因此本文提出的模型構建出第三組極限學習機用于預測? I?(t)的變化趨勢,預測模型對應的方程如下:
與治愈率和死亡率的變化趨勢不同,病毒的傳播與人類采取的防疫措施強度高度相關,因此對于傳播趨勢的分析具有一定的時效性,一部分歷史數(shù)據(jù)并不能反映當前的傳播能力。因此,在預測模型訓練過程中,利用時間窗對訓練數(shù)據(jù)的規(guī)模限制,確保預測模型的準確性。
利用本文提出的極限IR 預測模型,采用自2020 年1 月23 以來國家衛(wèi)生健康委員會公布的全國累積確診人數(shù)、累積死亡人數(shù)和累積治愈人數(shù)對模型的實際預測效果進行驗證。同時,通過與基于蒙特卡洛方法的SEIR 模型(SEIR-MC)和基于蒙特卡洛馬爾科夫方法的SIR 模型(SIR-MCMC)的預測結果進行比較,驗證算法的實時預測效果。
利用極限學習機,對治愈人群的變化趨勢和死亡人群的變化趨勢的預測結果如圖7 所示。由圖7a和圖7c 可以看出治愈率總體隨時間提升,而死亡率隨時間下降,受到實際的臨床經(jīng)驗進展影響變化趨勢發(fā)生波動。通過預測曲線,極限IR 模型能夠有效地對兩組概率的變化趨勢進行有效地預測,產(chǎn)生的預測時間曲線與實際的概率變化大致相同。因此,圖7b 和圖7d 中模型提供了治愈人數(shù)和死亡人數(shù)準確的估計,預測趨勢與實際的變化趨勢基本一致。
圖7 無風險人群預測效果比較
為說明極限IR 預測模型對現(xiàn)有確診人數(shù)的預測效果,圖8 展示了通過極限IR 模型進行的單步和多步預測結果。
圖8 極限IR 模型單步及多步預測結果
盡管在前10 天時預測結果發(fā)生了較大波動,單步預報值總體上能夠很好地描述疫情發(fā)展趨勢,且在第10 天后預測結果基本不受早期疫情傳染能力大幅波動的影響。尤其在疫情發(fā)展第10 天至第20 天,模型對現(xiàn)有確診人數(shù)的預測值與實際的確診人數(shù)基本吻合。此外,利用該模型得到的3 步預測值與5 步預測值也較好地預測了感染人數(shù)的變化趨勢。
為進一步說明極限IR 算法的效果,表1 比較了SEIR-MC 模型、SIR-MCMC 模型以及極限IR模型在2020 年2 月7 日至2 月16 日間的預測效果。由圖1 所示,SEIR-MC 模型在2020 年2 月7 日至2 月9 日的預測誤差較小,而在2020 年2 月10 日后疾病防控活動對疫情發(fā)展影響變強,利用該方法的預測誤差逐漸加大。另一方面,采用蒙特卡洛馬爾科夫方法在一定程度上可以實現(xiàn)對參數(shù)的動態(tài)評估,然而伴隨2020 年2 月12 日臨床診斷結果加入醫(yī)學診斷中的舉措,該方法對于參數(shù)的評估產(chǎn)生誤差。與SEIR-MC 和SIR-MCMC 不同,極限IR 模型利用神經(jīng)網(wǎng)絡實現(xiàn)動力學模型中時變參數(shù)的精確估計,因此從10 天數(shù)據(jù)的總體效果看,該模型能實現(xiàn)精度更高的實時預測,其誤差可以控制在10%以內(nèi)。
表1 2020 年2 月7 日至2 月16 日各模型的預測效果比較
本文針對冠狀病毒已知的傳播規(guī)律,提出用于疫情趨勢實時預測的極限IR 模型。通過建立時變假設下的傳染趨勢,死亡率和治愈率預測模型,實現(xiàn)動力學模型的及時更新,提供確診人數(shù)、新增死亡人數(shù)和新增治愈人數(shù)等指標預測值。通過實驗驗證,本文提出方法能夠達到較高的實時預測精度,為疫情趨勢分析提供有效的數(shù)據(jù)分析支持。
盡管模型在實時預測上存在很好的效果,但針對傳統(tǒng)的動力學模型做出了如下簡化:
1)考慮到潛伏期估計上的不確定性,在模型不將潛伏狀態(tài)作為獨立分析的狀態(tài)進行考慮。
2)對傳染病模型中的傳染態(tài)界定做出了調(diào)整,理由如下:第一,從疾病發(fā)生機理角度,將SIR 傳播模型中的感染態(tài)替換為傳染態(tài),可以避免因流感產(chǎn)生相似癥狀對模型中人群估計產(chǎn)生的干擾;第二,從防控過程角度,因確診人數(shù)在疫情防控工程中具有觀測精確的特征,通過確診人數(shù)可以為傳染態(tài)基數(shù)估計提供有力的數(shù)據(jù)依據(jù)。
3)機器學習方法的引入在增加了短期預測精度的同時增加了模型的復雜性并使擬合參數(shù)的價值降低。本文提出的方法通過引入極限學習機對IR 模型中的相關時變參數(shù)做定量趨勢預測。從模型的復雜性來看,盡管模型引入了相較于SIR 模型,SEIR 模型更加復雜的機器學習融合模型,但極限學習機作為單隱層少節(jié)點的神經(jīng)網(wǎng)絡模型,其空間復雜性能夠滿足實時預測的應用需求,且因訓練無迭代過程,模型預測的快速性得到保證。從擬合參數(shù)價值角度,機器學習融合的算法在治愈率、死亡率和病毒傳播能力的定性分析上能提供準確的動態(tài)變化信息,存在一定的分析價值,然而對于基本再生數(shù)這一類病毒機理性參數(shù),受到機器學習對對象機理簡化的影響,其分析需要進行進一步的研究。
綜上,盡管方法在傳染病機理傳遞的細節(jié)量化上不是至善之策,但對于疫情趨勢的及時把控、傳染態(tài)群體的統(tǒng)計特征變化及疫情發(fā)展的實施預測是現(xiàn)有較好的可取之舉。