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

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

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

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

1 模型

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖1 Jan.24th 時(shí)間段的擬合效果圖

圖2 Jan.28th 時(shí)間段的擬合效果圖

圖3 Jan.31th 時(shí)間段的擬合效果圖

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

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

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

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

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

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

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

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

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

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

4 結(jié) 束 語

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

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

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