基于時變參數(shù)-SIR 模型的COVID-19疫情評估和預(yù)測

喻 孜，張貴清，劉慶珍，呂忠全

(1. 南京林業(yè)大學理學院 南京 210037；2. 天津科技大學理學院 天津 濱海新區(qū) 300457；3. 中國醫(yī)學科學院血液病醫(yī)院 天津 和平區(qū) 300020)

截至2020 年1 月31 日晚12 點，全國確診新型冠狀病毒肺炎(COVID-19)病患接近1 萬。政府迅速采取了一系列措施來對抗疫情。從2020 年1 月23 日開始，地方政府陸續(xù)開始“封城”并通過媒體實時發(fā)布確診人數(shù)、專家通過網(wǎng)絡(luò)及媒體傳播防治疫情的注意事項、全社會自覺在家自我隔離、戴口罩出門、企業(yè)延期復(fù)工及人員聚集場所關(guān)停等。疫情的發(fā)展已經(jīng)嚴重影響到人們的日常工作與生活，并逐漸開始產(chǎn)生恐慌情緒。全社會都在關(guān)注兩個重要問題：一是確診人數(shù)將如何變化，什么時候會下降？二是迄今為止，政府所采用的防控措施，產(chǎn)生了多大的積極作用？回答這兩個問題，需要對疫情進行科學的預(yù)測和評估。

目前已經(jīng)有大量工作對COVID-19 疫情進行了研究。文獻[1]根據(jù)最早425 例確診病例數(shù)據(jù)，描述了病例特征，并估計了關(guān)鍵流行病學延遲時間分布情況，在病毒早期呈指數(shù)增長的初期，估計了傳染病倍增時間和基本再生數(shù)。文獻[2]發(fā)現(xiàn)了第二代病例，指出病毒會存在“人傳人”，對金銀潭醫(yī)院最開始收治患者時候的信息進行歸納研究，指出非華南海鮮市場暴露的病例的存在。文獻[3]基于SEIR 模型，根據(jù)2020 年1 月25 日之前的發(fā)展趨勢，估算得到了COVID-19 的再生數(shù)，并對疫情發(fā)展進行了預(yù)測。還有很多其他研究[4-9]，從不同角度對疫情的特點進行了評估。這些研究都對疫情的防治提供了重要的參考。在病毒傳播研究中[10-13]，由于SIR 動力學模型能給出較為清晰的邏輯關(guān)系和準確的趨勢預(yù)測，因而被廣泛采用。相比于之前的病毒傳播事件，COVID-19 疫情比較特殊：1) 相對于總的人口基數(shù)，目前的治愈人數(shù)和死亡人數(shù)都較低，因此出院人員(removed)比率可以忽略不記；2) 從2020 年1 月11 日，官方開始報道確診人數(shù)以來，存在人口流動，因此傳統(tǒng)SIR 模型認為總體環(huán)境是封閉(S+I+R=常數(shù))并不符合現(xiàn)實；3) 由于染病確診后將被隔離，所以確診病人并不形成新的感染；4) 此次病毒有較長的潛伏期，會帶來明顯的延遲效應(yīng)；5) 此次疫情突發(fā)，診斷力量在初期存在明顯不足。綜上所述，為了更為準確地對疫情進行預(yù)測，需要對SIR 模型進行修正。本文根據(jù)實際情況，對SIR 模型進行了適當修正，聚焦國內(nèi)確診人數(shù)的變化，對趨勢進行預(yù)測和判斷，并根據(jù)結(jié)果對現(xiàn)有政策進行評估。

1 模型

本文根據(jù)實際疫情對SIR 模型進行了一些修正，描述如下：

1) 由于不存在封閉情況，考慮開放體系。

2) 新型冠狀病毒的治愈人數(shù)和死亡人數(shù)相對較小，因此只考慮Susceptible(易感)和Infected(感染)兩類人群。

3) 目前數(shù)據(jù)以天為單位發(fā)布，因此不考慮連續(xù)變化情況，只考慮離散的方程。

4) 當確診后，病人立即隔離，因此不會作為新的感染源。因此感染源應(yīng)該是第n+1 天和第n 天診斷的病人差。假設(shè)平均染病者(I)會使k 個人成為易感者(k 為易感再生數(shù))，則每日新增的感染者為k(In+1?In)。

5) 假設(shè)每天有 μ1的概率(當日感染率)易感者會成為感染者，則當日新增病人為 μ1Sn。 其中 kμ1，描述了當日病人再生數(shù)，后文簡稱當日再生數(shù)。

6) 病情會有約10 天的潛伏期，在初期出現(xiàn)癥狀后，很多暴露者才會去尋求診斷，同時主要發(fā)病地武漢的初期診斷能力有限，這都會帶來一個延遲效應(yīng)。取τ=4，代表平均延遲天數(shù)。在平均4 天內(nèi)的潛伏者會以 μ2的概率(潛伏感染率)被確診為感染者。則每日潛伏感染人數(shù)為 μ2(Sn?Sn?τ)。 kμ2描述了由于潛伏原因?qū)е碌牟∪嗽偕鷶?shù)，后文簡稱潛伏再生數(shù)。

基于上述考慮，在每日模型迭代中，各部分更新關(guān)系如下：

為了求解式(1)，需要知道(k, μ1, μ2)3 個模型參數(shù)[14-18]，以及I 和S 的初值。I(0)通過官方數(shù)據(jù)給出(數(shù)據(jù)來源于官方網(wǎng)站：http://www.nhc.gov.cn/)，剩下3 個參數(shù)以及S(0)則沒有信息。本文先給定3 個參數(shù)和S(0)的初值，求解方程，將得出的解與真實數(shù)據(jù)進行比較，通過梯度下降法尋找方差的最小值，從而確定模型參數(shù)。由于梯度下降法對初值極為敏感，因此，本文還根據(jù)經(jīng)驗在一定初值范圍內(nèi)進行多次搜索，保證得到的是最優(yōu)解。

模型參數(shù)(k, μ1, μ2)是根據(jù)一個階段的數(shù)據(jù)求出，因而能夠反應(yīng)病毒在該階段的傳播趨勢。通過截取不同階段的數(shù)據(jù)來求解不同階段的(k, μ1,μ2)，就能分析病毒演化趨勢是否在外力作用下發(fā)生改變。最后再通過模型參數(shù)(k, μ1, μ2)的變化規(guī)律，對未來進行預(yù)測。

官方從2020 年1 月11 日開始公布數(shù)據(jù)，前4 天均無確診人數(shù)變化。本文選擇2020 年1 月15 日作為數(shù)據(jù)研究的起始點，然后變化數(shù)據(jù)的結(jié)束點來截取不同時間段的數(shù)據(jù)進行求解。為了描述方便，后文對時間段使用結(jié)束日期的英文簡寫來標識。例如1 月15 日至1 月24 日這個時間段，簡稱為Jan.24th 時間段，1 月15 日至1 月28 日這個時間段，簡稱為Jan.28th 時間段。Jan.24th 時間段是文章研究的第一個時間段，此后每增加1 天構(gòu)成一個新的時間段。這是基于以下兩點考慮：1) 從Jan.24th 時間段開始，已經(jīng)積攢了足夠多的數(shù)據(jù)(10 天的數(shù)據(jù))；2) 政府從1 月11 日開始公布確診人數(shù)并采取措施，趨勢的形成需要一個時間，并且病毒的最長潛伏期約為14 天。需要指出，在理想狀況下，S(0)應(yīng)該是一個確定值，然而，在開放體系下，通過現(xiàn)有的數(shù)據(jù)無法評估S(0)的真實值。因此只能通過式(1)的解與數(shù)據(jù)擬合來反向求解S(0)。這會造成不同階段的S(0)略有不同。然而，由于本文采用了大范圍梯度下降搜索，因此每一組參數(shù)都是在與現(xiàn)有數(shù)據(jù)擬合最優(yōu)情況下得到的，因此文章求解得到的(k, μ1, μ2)能體現(xiàn)數(shù)據(jù)背后隱藏的病毒演化趨勢。

2 模型參數(shù)的影響

2.1 易感再生數(shù)、當日感染率、潛伏感染率的影響

圖1～圖3 給出了不同時間段的擬合效果，其中曲線是模型計算結(jié)果，圓點代表官方統(tǒng)計數(shù)據(jù)。從圖中發(fā)現(xiàn)模型的計算結(jié)果和實際公布的數(shù)據(jù)吻合較好，說明模型是合理的。同時也發(fā)現(xiàn)易感再生數(shù)、當日感染率、潛伏感染率參數(shù)也確實發(fā)生了變化，并且是在不斷減小，這和政府的防控措施及民眾的居家防護有直接關(guān)系。圖4 給出了3 個時間段確診人數(shù)隨時間變化的趨勢圖。從圖中可以清楚發(fā)現(xiàn)，疾病的傳播趨勢確實是在減緩的，如果按照Jan.24th 時間段的疫情發(fā)展趨勢，那么在2020 年1 月31 日的確診人數(shù)將會達到近40 000，同樣，如果按照Jan.28th 時間段的疫情發(fā)展趨勢，那么1 月31 日的確診人數(shù)將會達到15 000 左右，明顯高于實際疫情的發(fā)展情況。所以圖4 表明政府的防控措施和民眾的居家防護對于疫情的控制確實效果明顯(確診人數(shù)按照趨勢變化減少了將近1/2)。

圖1 Jan.24th 時間段的擬合效果圖

圖2 Jan.28th 時間段的擬合效果圖

圖3 Jan.31th 時間段的擬合效果圖

圖4 不同時間段確診人數(shù)的走勢圖

2.2 易感再生數(shù)、當日再生數(shù)、潛伏再生數(shù)的變化趨勢

為了更清楚地展示從初始階段(Jan.24th 時間段)起的趨勢變化情況，本文從第一個階段(2020 年1 月15 日?1 月24 日)開始，每增加一天構(gòu)成一個新的時間段進行求解，得到了易感再生數(shù)、當日感染再生數(shù)和潛伏感染再生數(shù)隨時間的演化趨勢。結(jié)果如圖5～圖7 所示，圖中的橫坐標代表了時間段結(jié)束日期，圖中x 的值取時間段結(jié)束日期與1 月23 日間隔的天數(shù)(例如1 月24 日對應(yīng)x=1，1 月25 日對應(yīng)x=2)。從圖中可以看到3 個參數(shù)都有下降趨勢。同時，趨勢每天都在改變。說明政府的防控措施和民眾的居家隔離確實產(chǎn)生了強力效果。需要注意，2020 年1 月27 日，多地醫(yī)療工作者馳援武漢，使得檢測能力得到大幅提升，所以導(dǎo)致1 月27 和28 日的確診人數(shù)大幅上升，出現(xiàn)了明顯的和前期不一樣的易感再生數(shù)、當日再生數(shù)和潛伏再生數(shù)。這些參數(shù)的波動，也從側(cè)面證明了本文模型能夠敏感地體現(xiàn)外力對趨勢的改變。另一方面，也說明參數(shù)固定不變的動力學模型無法真實體現(xiàn)此次 COVID-19 疫情的變化情況。

圖5 易感再生數(shù)隨著時間的演化情況(去除了27 日的波動較大點)

圖6 當日再生數(shù)隨著時間的演化情況

3 基于最優(yōu)參數(shù)的疫情預(yù)測

為了結(jié)合現(xiàn)有趨勢對未來進行預(yù)測，本文去除波動數(shù)據(jù)的影響對參數(shù)進行了擬合，從而得到了參數(shù)的變化趨勢方程如下：

式中，x 代表與1 月23 日間隔的天數(shù)，從1 月24 日開始取1。參數(shù)按照式(2)的規(guī)律變化，代入式(1)中進行求解，預(yù)測1 月31 日以后的數(shù)據(jù)，如圖8 所示。從圖中可以看出26 天后，也就是2020 年2 月9 日左右應(yīng)該會出現(xiàn)確診總數(shù)的峰值，同時確診人數(shù)也達到最大的54 000 左右。隨后確診人數(shù)將會下降。

圖7 潛伏再生數(shù)隨著時間的演化情況(去除27，28 日上下兩個波動點后，發(fā)現(xiàn)直線擬合效果最佳，因此用直線擬合)

圖8 圖中給出了在模型最優(yōu)參數(shù)條件下的模型估計值

4 結(jié) 束 語

本文基于修正SIR 模型研究了COVID-19 疫情發(fā)展到2020 年2 月1 日(政府隔天公布數(shù)據(jù))為止的趨勢和政策的效果。模型的參數(shù)根據(jù)實時數(shù)據(jù)擬合得到，因而修正過后的時變參數(shù)-SIR 模型結(jié)果能夠反應(yīng)不同階段疫情的趨勢變化。通過趨勢變化能夠解讀政府的防控措施產(chǎn)生的效果。結(jié)果表明，在2020 年1 月27 日和28 日，由于全國醫(yī)療系統(tǒng)馳援武漢使得確診能力得到了增強，一些積壓和潛伏病人得到收治，因此模型參數(shù)出現(xiàn)擾動。其他時間段，參數(shù)變化趨勢明顯。2020 年1 月24 日之后，政府的的防控措施有效地降低了病毒蔓延趨勢，使得染病人數(shù)的變化趨勢逐日下降。易感再生數(shù)、當日感染率和潛伏感染率都大幅度降低?；诋斍暗内厔輰σ咔榘l(fā)展進行了預(yù)測。結(jié)果表明，在2 月9 日左右，疫情會達到高峰，隨后確診人數(shù)將會出現(xiàn)下降。

從動力學方程的解和實際數(shù)據(jù)的比較來看，本文提出的時變參數(shù)-SIR 模型能夠及時反應(yīng)現(xiàn)有趨勢變化，甚至體現(xiàn)一些突發(fā)情況(2020 年1 月27 日全國醫(yī)療系統(tǒng)馳援武漢)帶來的影響。因此，本文基于現(xiàn)有數(shù)據(jù)對趨勢變化的分析結(jié)論應(yīng)該是有效的。從預(yù)測角度來說，動力學預(yù)測方案本身就會存在一定的不確定性。然而，本文并沒有采用固定參數(shù)去求解動力學方程，而是基于前期積累的趨勢，采用動態(tài)參數(shù)求解，這在一定程度上減少了動力學方程的不確定性。因而，本文預(yù)測的趨勢應(yīng)該能夠提供有益的參考價值。本文基于COVID-19 病毒的流行特點，對SIR 模型進行的理論修正，以及通過數(shù)據(jù)確定動態(tài)參數(shù)進行預(yù)測的方案，會對病毒的研究提供積極的理論參考。

本研究工作還得到天津科技大學青年教師創(chuàng)新基金(2014CXLG23)的資助，在此表示感謝。