劉元浩,曹婍,沈華偉,黃俊杰,程學旗
(中國科學院計算技術研究所數據智能系統研究中心,北京 100190)
自2019年底以來的幾個月內,新型冠狀病毒肺炎COVID-19 在全世界范圍內廣泛流行,截至2020 年4月7 日,全球COVID-19 累計確診人數已達1,279,722例,并仍在持續快速增長。疫情的持續蔓延對人們的生命安全造成巨大威脅,也對國家醫療建設、物資調配、隔離管控等方面帶來挑戰。在此背景下,采用數學方法對疫情傳播進行建模并對確診病例數的增長進行及時準確地預測對于疫情防控具有重要意義。一方面,對疫情傳播進行準確預測,對于醫療衛生資源的分配、防控重點的調整等具有重要的參考價值。另一方面,在防疫工作進行過程中的重要時間節點前后,對確診病例數增長趨勢的變化進行對比,能夠有效地對防控措施有效性進行合理評估。
對于疫情傳播預測,最常用的研究框架是傳染病模型。傳統的常微分方程傳染病模型假設人群總數恒定且人群均勻混合[1],通過對人群中處于各個狀態的人數及各狀態間的相互轉換速率進行建模,推算疫情發展走勢。常見的傳染病模型根據人群劃分的不同及人群轉換的不同,包括SI[2]、SIS[3]、SIR[4]、SEIR[5]等。傳染病模型從傳染病傳播動力學的角度進行考慮,能對疫情短期內發展趨勢進行較好的模擬,但其總人數恒定且人群均勻混合的理想化假設使其應用場景受到局限。
疫情感染人數預測的另一方法是時間序列預測。疫情傳播情況會隨著時間推移而不斷演變,疫情感染人群數可以形式化為一種時間序列,并采用時間序列分析與建模的方式進行預測。線性時間序列分析模型包括自回歸(Auto-regression)模型和移動平均(Moving Average)模型[6],以及基于兩者組合而成的自回歸移動平均(Auto-regression Moving Average)模型[7]。向量自回歸(Vector Autoregressive)模型等非線性時間序列模型以及基于深度神經網絡的RNN[8],LSTM[9],TCN[10]等模型也在時間序列分析問題上有優秀的表現。采用時間序列分析模型進行疫情預測,能夠通過簡單的模型建模疫情發展的時間序列當前值與序列歷史信息間的關系,對疫情走勢做出預測。但由于缺乏對疫情的傳染性、爆發性、衰減性等特性的認識與建模,對疫情確診人數的預測仍有一定的局限。此外疫情前期可用數據有限,也給時間序列模型的學習造成了很大困難。
本文采用自增強泊松過程(RPP)模型[11]對疫情確診人數變化趨勢進行預測,該模型將病毒感染人群的動態過程建模為不均勻泊松過程,通過對病毒傳染性、級聯傳染的自增強效應和病毒傳播的時效性等三個因子進行建模,對疫情傳播過程中的關鍵因子進行刻畫,以解決上述模型中出現的問題,并使用本次COVID-19 疫情傳播數據進行實驗,證明模型的有效性。
自COVID-19 疫情發生以來,世界各地學者紛紛嘗試對疫情的發展趨勢展開研究和分析。其中以SEIR 模型為代表的微分方程傳染病模型占據了疫情趨勢預測工作的主要部分。SEIR 模型將人群劃分為易感者(Susceptible)、潛伏期感染者(Exposed)、感染者(Infected)、治愈者(Recovered)四個群體,以微分方程描述四個狀態間的轉換關系。2020年1月31日,香港大學學者Joseph T Wu 應用SEIR 模型,利用武漢早期病例數,推測疫情會在4 月達到高峰[12]。肖燕妮教授團隊同樣基于SEIR 模型,考慮跟蹤隔離等管控措施,對疫情的走勢和管控舉措的有效性進行了分析[13]。2 月28 日,鐘南山院士團隊考慮地區間人口流動對SEIR 模型進行改進,通過對實施管控措施時間的調整,論證了控制措施對于減少最終COVID-19 流行病的規模是必不可少的[14]。西安交通大學[15]、北京郵電大學[16]等國內研究機構也通過傳染病動力學建模對COVID-19疫情走勢做出了預測。上述基于傳染病模型對疫情預測分析的工作被證明可以較準確地反映小范圍空間在短期內的疫情走勢,但由于這類模型對初始參數敏感,且基于人群均勻接觸的理想假設,難以應對不同地區不同時間帶來的復雜疫情發展趨勢變化。
基于時間序列分析的疫情預測分析在流感等傳染病預測領域多有應用[17][18][19]。線性時間序列模型如Pinto[20]模型假設未來的序列值為歷史序列值的線性組合,從而通過歷史確診人數對未來的確診人數進行預測。然而線性時間序列模型在應用中面臨諸多局限,基于深度學習技術的循環神經網絡(Recurrent Neural Network,RNN)[8]模型解決了這一問題,但也因對長序列進行學習時會出現梯度爆炸或梯度消失現象,從而無法對長序列建模。時序卷積網絡(Temporal Convolutional Networks,TCN)[10]模型通過因果空洞卷積的設計,提取序列局部特征的同時增大感受野,實現了對長時間序列的有效處理。上述提到的時間序列模型,能對時間序列進行建模并預測。然而這些時間序列模型對疫情傳播的傳染性、爆發性、衰減性的關鍵性質缺乏認識,且要求有足夠的訓練數據用來學習模型參數,這使得其在疫情預測應用領域存在一定局限。
我們使用疫情傳播的動態過程刻畫人群中個體被病毒感染并發病這一事件的發生過程。對于某傳染病d,我們將其在時間段[0,T]內的疾病感染人群動態變化過程表示為個體染病事件發生的時間序列:
其中nd表示T時刻內被疾病感染的人群總人數,表示第i個染病事件發生的時間。不失一般性,令0 ≤≤T。
為了建模疫情傳播的動態過程中個體染病事件發生的速率,我們考察疾病傳播過程中的三大現象:(1)病毒傳染性,即病毒自身的傳染性對最終的感染人數起決定作用;(2)級聯傳播所帶來的自增強效應,即病毒當前的感染人數越多越容易進行新的傳播感染;(3)病毒傳播的時效性,即隨著時間推移,病毒感染人群繼續感染他人的可能性會下降。綜合考慮這三個現象,我們采用自增強泊松過程(Reinforced Poisson Process,RPP)[11]來建模疾病感染人群的動態過程。具體而言,對于某個傳染病d,其感染人群的動態過程建模為一個速率為的泊松過程。其中,λd是病毒自身的傳染性,松弛函數fd(t;θd)刻畫病毒傳播的速率隨時間演變過程。θd是松弛函數的參數,id(t) 表示病毒d在時刻t已經感染的人群數量。我們假定所有的病毒在開始感染前,都有一定初始感染人數m。因此,在第i- 1 次真實感染事件發生后到第i次真實感染事件發生前的時間段內,我們有id(t) =m+i- 1(1 ≤i≤nd)。相應地,在第nd次真實感染事件發生后到時刻T之前,我們有id(t) =m+nd。
對于疫情預測,我們采用對數正態松弛函數
作為刻畫病毒傳播時效性的松弛函數。此時松弛函數的參數θd被替換為對數正態函數的均值μd和方差σd。
整個疾病感染人群的動態過程可以表示為如圖1所示的產生式概率圖模型。
圖1 疾病感染人群動態過程的產生式概率圖模型[11]
兩次連續感染事件之間的時間間隔長度服從不均勻泊松過程。因此,設第i- 1 次真實感染事件的發生時刻為,那么第i次真實感染事件在時刻發生的概率滿足:
在第nd次真實感染事件發生時刻和觀測時刻T之間沒有感染事件發生的概率為:
那么,在時間間隔[0,T]內觀測到病毒d的染病人群動態過程的似然為
其中,Fd(t; μd,σd )是松弛函數fd(t;μd,σd)的累積分布函數。
我們通過最大似然估計,學習病毒d的參數λd,μd和σd。令似然函數導數為零,可直接求得參數λd的最大似然估計值
對于μd和σd,我們使用梯度下降法最大化似然函數,梯度
其中,φ是標準正態分布的概率密度函數,
根據泊松過程的速率函數和對應的微分方程求解,我們得到病毒感染人群的預測函數:
4.1.1 SEIR型流行病模型
傳染病學模型采用肖燕妮教授團隊的工作[13],該模型在傳統SEIR模型對人群的“易感者-暴露者-感染者-治愈者”劃分的基礎上,結合COVID-19 的實際情況與諸如檢疫,隔離和治療等干預措施,將人群分為易感者(S),暴露者(E),潛伏傳染者(未表現出癥狀但有傳染性)(A),具有癥狀的傳染者(I),住院患者(H)和康復者(R),并進一步劃分出被隔離的易感者(Sq)和被隔離的暴露者(Eq)。不同人群間的狀態轉換方程如下:
通過對模型設定合適的參數和初始值來推算疫情累計確診人數C=I+H+R。
4.1.2 時間序列模型
Pinto模型:采用該模型作為線性時間序列模型的代表。該模型劃定待預測時刻前的一段時間T 作為觀測窗口,將采樣窗口劃分為大量的采樣間隔,采用每個采樣間隔內的新增確診人數作為模型的輸入,通過簡單的多元線性組合給出模型的預測值[20]。
TCN模型:非線性時間序列模型采用時序卷積網絡(TCN)模型。該模型通過卷積神經網絡自動提取疫情發展歷史序列中的重要特征,并通過因果空洞卷積提升增大了感受野,從而可以觀測更久的歷史序列[10]。
本文實驗采用中國1 月20 日至3 月15 日共計56天的COVID-19 每日確診人數[21][22]作為實驗數據。數據范圍基本涵蓋了全國自疫情開始流行至爆發到基本得到控制的全過程。同時考慮到3 月16 日后國內新增確診病例來源以境外輸入為主,因此將其排除,最大程度上避免了境外輸入病例對實驗結果的影響。
考慮到疫情傳播具有地區性,不同地區疫情出現時間存在先后差異,疫情發展速度也可能不相同。我國的疫情傳播呈現出明顯的“武漢市-湖北省-全國其他地區”地區劃分:一方面體現在地區間的隔離上,自1 月23 日起武漢開始全面封城,而湖北省也率先實行了較為嚴格的出入管制措施,最大程度上減少了感染病例的流入和流出;另一方面體現在疫情傳播的時間先后和傳播的規模上,國內疫情最先發現于湖北省武漢市,隨后蔓延至湖北省和全國其他地區,我國及時采取措施將疫情大規模傳播范圍盡可能地控制在了小范圍內,截至3 月15 日,全國近84%的確診病例發現于湖北省,而其中又有近74%的病例位于武漢市。
因此本文將全國確診人數數據劃分為“全國”、“全國(除湖北)”、“湖北(除武漢)”、“武漢”四個地區層次。
圖2是四個地區的累計確診人數隨時間變化的曲線(為保證模型預測結果體現疫情發展的真正趨勢,我們排除掉2月12日新增確診人數中的臨床診斷病例數[22])。不難發現各個地區劃分下,疫情整體發展趨勢基本一致,但存在總量和增長速度等方面的明顯差異。
圖2 各地區確診人數隨時間變化曲線
除了地區劃分,疫情趨勢在不同時間階段的表現也有差異。如圖中矩形框所標識,累計確診人數的變化在時間上較為明顯地呈現出三個階段:(1)前期——加速增長階段,在疫情流行初期,累計確診人數增速持續上升,圖線呈下凸經過矩形框的右下部分;(2)中期——增速穩定階段,隨著疫情發展與防控措施的實行,每日新增確診人數基本維持不變,圖線基本沿矩形對角線呈直線;(3)后期——增速放緩階段,后期疫情得以控制,確診人數增速迅速放緩,圖線從扁平矩形框的左上部分經過。為量化表示三個階段的特點,我們對地區a的疫情發展階段u計算平均增長系數,
其中ca(i)為地區a第i天的新增確診人數,Tu為階段u的天數。在計算時,我們將累計確診人數曲線進行了平滑處理以避免每日新增確診病例數波動的影響。四個地區的各階段劃分與平均增長系數見表1。
表1 各地區疫情發展趨勢的階段劃分
我們分別考察模型在各地區不同時間階段的預測表現,作為衡量模型在不同環境下預測能力的依據。
由于確診病例數以1天為單位時間統計,因此RPP模型的最小時間單位為1天,當天所有新增病例計為同時發生。由于疫情發展的情況會隨時間變化,為保證模型較好地反映近期疫情的走勢,我們沒有使用預測時間之前的所有數據,而是在4 ≤T≤15范圍內通過搜索確定觀測窗口T大小。初始感染人數m= 20。
SEIR模型參數值設定采用文獻[13]中的取值,該文參數由武漢市早期疫情數據模擬獲得。我們使用原文方法求取了不同地區劃分下的模型參數,取值見表2。模型的初始值獲取自國家衛健委的報道數據[22],未明確報道的狀態初值由預測時間前一段時間疫情相關數據通過最大似然估計得出。
表2 SEIR模型參數取值
時間序列模型的觀測窗口大小同樣通過搜索確定,從而選取合適的觀測歷史長度同時保證一定的訓練集體量。采樣間隔設置為1天。
我們計算預測結果的MAPE(Mean Absolute Percentage Error,平均絕對百分比誤差)以衡量模型的預測能力。MAPE的計算公式為
其中n 為預測時間段的天數,Ct為第t天的累計確診人數,為其預測值。
我們使用不同模型在全國、全國(除湖北)、湖北(除武漢)、武漢四個數據集上進行實驗,對比不同地區間的差異對模型預測效果的影響。考慮到訓練數據量和預測時段可能對各模型的預測效果產生不同的影響,因此我們分別在疫情前半段與后半段進行實驗,使用1月31日之前和2月10日之前的數據訓練模型,分別對隨后一周(即2 月1 日至2 月7 日和2 月11日至2 月17 日)的累計確診人數進行預測,誤差結果如表3、表4。
表3 各模型2月1日至2月7日預測結果MAPE
表4 各模型2月11日至2月17日預測結果MAPE
RPP 模型、Pinto 模型與TCN 模型對不同地區的疫情預測效果均比較穩定。相較于Pinto 模型與TCN模型僅對歷史確診人數序列進行分析,RPP 模型對疫情傳播的關鍵因子進行了建模,其預測結果明顯優于其他模型。
SEIR 模型的預測效果在不同地區差異較大。這是由于SEIR 模型假設人群均勻混合,在全國各地采取封城措施相互隔離的情況下,絕大多數的感染人群的活動實際上被限制在了湖北省和武漢市內,這與人群均勻混合的假設高度不符,從源頭上限制了SEIR模型的表現。
根據表1,我們從時間上將疫情的發展過程劃分為前期、中期和后期三個階段。這一部分我們從時間劃分的角度,考察模型“階段內預測”和“跨階段預測”的效果。
5.2.1 階段內預測
疫情發展的每一個階段都有其特定的發展趨勢和規律,對這些規律的把握能力是模型完成精準預測的基本要求。這一部分實驗分別使用前期,中期,后期的前半段數據作為訓練數據,預測同時期后半段的累計確診人數。以武漢地區為例,各模型預測結果如圖3、圖4、圖5。
圖3 武漢市前期累計確診人數預測結果
圖4 武漢市中期累計確診人數預測結果
圖5 武漢市后期累計確診人數預測結果
由圖可以看出,在疫情前期數據量較少、數據規律性較弱時,Pinto 模型與TCN 模型難以通過少量數據掌握疫情發展的整體趨勢,因此預測效果較差。而RPP 模型與SEIR 模型通過對疫情發展固有性質的建模,可以較好地模擬確診人數加速增長的趨勢。
中后期RPP、Pinto 與TCN 模型對確診人數增速保持穩定至減緩的趨勢能夠進行較好地模擬。SEIR 模型由于其模型假設所有人都會暴露在被傳染的風險下,因此在人口總數很大時,最終累計確診人數也會變得過高。因此我們在中后期不再對其進行對比。
我們在各地區數據集上進行了相同的階段內預測實驗,綜合平均誤差結果如表5。
表5 各模型階段內實驗累計確診人數預測MAPE
可以看出RPP 模型在各個時間階段內都能準確地對確診人數進行預測,對階段內疫情發展趨勢能夠進行較好地把握。
5.2.2 跨階段預測
由于疫情發展的不同階段趨勢各不相同,從而意味著應該采取不同的防疫措施。也就是說,模型對于趨勢轉換的準確預測能力十分重要。因此這一部分實驗分別使用前期和中期的確診人數訓練模型,預測其下一個時期的疫情趨勢。
同樣以武漢地區為例,各模型預測結果如圖6、圖7。
圖6 武漢市前期-中期累計確診人數預測結果
圖7 武漢市中期-后期累計確診人數預測結果
由圖6、圖7 可知,相較于Pinto 和TCN 模型,RPP模型在中期更能把握確診人數增速保持穩定隨后趨于下降的趨勢。而后期RPP 和Pinto 模型都能較好地模擬增速迅速下降的趨勢,而TCN模型的預測結果則傾向于保持增速持續增長。
我們在各地區數據集上進行相同的跨階段預測實驗,計算預測結果的平均增長系數(p ?_u^a ) ?,再與表1中的實際值進行對比,計算MAPE,結果如表6。
表6 各模型跨階段實驗平均增長系數的MAPE
可知RPP 模型在判斷階段變化時表現明顯優于Pinto 與TCN 模型,在各時間階段都能很好地預測疫情發展趨勢的階段性變化。
我們將本文方法投入實際應用,自1 月29 日起先后對中國、日本、韓國、意大利、美國等九個國家共12個地區的疫情確診人數進行預測,累計確診人數平均誤差率小于0.5%。預測結果發布于中科天璣智疫通線上平臺(https://ncov.ictbda.com/#/,效果如圖8)
圖8 在線系統疫情預測效果
本文應用基于自增強泊松過程(RPP)的模型來預測COVID-19的疫情確診病例數。我們的實驗結果表明,RPP模型在預測疫情確診人數的任務中明顯優于傳統的傳染病模型和時間序列分析模型。在空間上,RPP模型克服了SEIR模型基于人群均勻混合的局限,在各尺度的地理區域都有穩定且準確的預測結果。在時間上,一方面,RPP模型解決了SEIR模型在人口總數很大時累計確診數持續增長的問題;另一方面,RPP模型通過建模疫情發展過程中的關鍵因素,擺脫了時間序列分析模型僅對歷史數據建模的局限性,從而對疫情發展各個階段的疫情走勢能夠進行更精確的預測,并且能準確把握疫情發展的重要階段性變化,其結果在實際應用更具有參考價值。
本文的方法也存在進一步優化的空間,本文假設感染速率與當前感染人數成正比,并使用松弛函數從整體上描述部分感染者被隔離或被治愈等情況造成的感染者總體影響力下降。未來將考慮使用Hawkes過程進行建模,細化不同狀態感染者對疾病感染速率的影響。