基于改進(jìn)MCMC的疫情人群管控模型

鄭辰彥,王露露，胡益?zhèn)?/p>

(華東交通大學(xué) 經(jīng)濟(jì)管理學(xué)院，南昌 330013)

近年來(lái)重大突發(fā)疫情在全球范圍內(nèi)頻發(fā)，如2003年的“非典”(SARS)、2009年的美國(guó)大流感(H1N1)、2015年的中東呼吸綜合征(MERS)、2018—2019年的埃博拉(Ebola)和當(dāng)前全球蔓延的新冠肺炎(NCP)，這些傳染性極強(qiáng)、致死率極高的疫情頻發(fā)說(shuō)明人類與病毒的斗爭(zhēng)將是長(zhǎng)久的、持續(xù)的。2019年發(fā)生的新冠肺炎疫情因傳播速度快、傳播規(guī)律復(fù)雜、潛伏期長(zhǎng)等特征迅速席卷全球，這必然給人民生命財(cái)產(chǎn)安全帶來(lái)巨大的損失，各國(guó)政府管理部門面臨艱難的挑戰(zhàn)。世界各國(guó)由于社會(huì)文化、政府體制、政府能力各不相同，導(dǎo)致各國(guó)政府處理疫情的措施也不一樣，取得的效果也不一樣。中國(guó)在疫情暴發(fā)后，很快地控制了疫情的發(fā)展，從當(dāng)前的效果來(lái)看，應(yīng)該說(shuō)是世界上疫情控制得最好的國(guó)家。美國(guó)和印度是采用市場(chǎng)管理的辦法，現(xiàn)在看來(lái)是處理得最差的國(guó)家。韓國(guó)等結(jié)合其國(guó)情，采取了人群管制的措施，對(duì)感染者人群進(jìn)行了隔離治療，也取得了較好的成績(jī)。此次新冠肺炎，是人類首次遇到的傳染性強(qiáng)、致死率高的傳染病。與其他傳染病不同，新冠肺炎的傳播規(guī)律復(fù)雜，除了存在易感者、潛伏者、感染者和康復(fù)者，甚至還存在大量的疑似患者，這些疑似者中就有一定的感染者，導(dǎo)致部分疑似者存在著傳染性。另外，潛伏者也存在傳染性。對(duì)于快速傳播的高致病性傳染病來(lái)說(shuō)，信息獲取的滯后所造成的后果是不可掌握的。因此，在不同的疫情管控策略下，找到符合新冠肺炎傳播規(guī)律的傳染病模型來(lái)預(yù)測(cè)疫情的發(fā)展趨勢(shì)，進(jìn)而能夠及時(shí)發(fā)現(xiàn)、追蹤和預(yù)測(cè)傳染病的傳播規(guī)律，為實(shí)行有效的公共衛(wèi)生干預(yù)措施提供數(shù)據(jù)支撐。這在沒(méi)有治療的特效藥與接種疫苗的情況下，具有十分重要的理論意義與現(xiàn)實(shí)指導(dǎo)意義。與本文相關(guān)的研究主要集中于傳染病模型的建立、傳染病模型的參數(shù)估計(jì)和傳染病臨界閾值的推導(dǎo)。

在傳染病模型方面，Chen等研究了具有感染年齡和飽和發(fā)病率的SIR流行模型[1]。Cai等研究了具有比率依賴的傳染病模型的全局動(dòng)力學(xué)和相應(yīng)的SIRS隨機(jī)微分方程[2]。Wu等在2019-nCOV的疫情背景下建立了SEIR傳播動(dòng)力學(xué)模型[3]。每種疫情的發(fā)生發(fā)展總有一些特殊的情景，便有學(xué)者將經(jīng)典SEIR模型根據(jù)實(shí)際情況進(jìn)行改進(jìn)，Liu等提出了一種具有年齡依賴性潛伏期和復(fù)發(fā)性疾病的SEIR流行模型[4]。Alonso-Quesada等研究了一類SEIR流行病模型的離散化與控制問(wèn)題[5]。趙英英等考慮標(biāo)準(zhǔn)發(fā)生率和信息干預(yù)因素，建立了SIRS傳染病模型[6]。王改霞等構(gòu)建了 SIQRS型傳播模型來(lái)分析無(wú)病平衡點(diǎn)的局部穩(wěn)定性以及接種最優(yōu)策略[7]。梅珊等考慮個(gè)人行為、交通信息和地理信息因素建立了空氣傳播動(dòng)力學(xué)模型[8]。李勇建等建立了Petri網(wǎng)傳染病傳播模型[9]。

在傳染病模型的參數(shù)估計(jì)方面，George等以模型擬合值和實(shí)際傳染病數(shù)據(jù)之間的隨機(jī)誤差平方和最小化原則，估計(jì)出傳染病模型的相關(guān)參數(shù)[10]。Aranda等建立了寨卡病毒感染者進(jìn)化的數(shù)學(xué)模型，使用哥倫比亞2016年寨卡流行的真實(shí)數(shù)據(jù)對(duì)傳播動(dòng)力學(xué)模型進(jìn)行參數(shù)估計(jì)[11]。Li等提出了一種基于遺傳算法的參數(shù)估計(jì)方法來(lái)估計(jì)埃博拉病毒模型的參數(shù)[12]。Korolev通過(guò)幾種非線性SUR估計(jì)出流行病學(xué)模型的參數(shù)[13]。Kim等建立了SEIQR傳染病模型預(yù)測(cè)韓國(guó)疫情發(fā)展，通過(guò)最小二乘法估計(jì)模型的參數(shù)[14]。李倩等基于新型冠狀病毒疫情背景，建立了SEIR型疫情傳播動(dòng)力學(xué)模型，使用了最小二乘法對(duì)模型參數(shù)進(jìn)行了估計(jì)[15]。崔景安等通過(guò)實(shí)時(shí)疫情數(shù)據(jù)對(duì)模型參數(shù)進(jìn)行估計(jì)[16]。郭尊光等建立了SEIR反應(yīng)擴(kuò)散傳染病模型，并用最小二乘優(yōu)化方法對(duì)模型中的傳染率參數(shù)進(jìn)行估計(jì)[17]。

基本再生數(shù)是傳染病領(lǐng)域的一個(gè)重要閾值參數(shù)，用來(lái)衡量傳染病的傳播能力和發(fā)展趨勢(shì)。一般在傳染病發(fā)展初期，基本再生數(shù)是一個(gè)常數(shù)，隨著管控措施的實(shí)施，基本再生數(shù)動(dòng)態(tài)變化，這時(shí)候稱之為有效再生數(shù)。在傳染病的臨界閾值方面，Lim等針對(duì)特定腸道病毒類型，依據(jù)疫情暴發(fā)的最初生長(zhǎng)階段累積報(bào)告的病例數(shù)，來(lái)估計(jì)基本再生數(shù)的具體數(shù)值[18]。Kucharski等通過(guò)統(tǒng)計(jì)隨機(jī)模擬的方法預(yù)測(cè)病例數(shù)，根據(jù)病例數(shù)落入?yún)^(qū)間的比例來(lái)估計(jì)基本再生數(shù)[19]。Zhao等首先估計(jì)代際間隔，然后通過(guò)泊松先驗(yàn)的對(duì)數(shù)似然估計(jì)來(lái)估計(jì)指數(shù)增長(zhǎng)率，最后通過(guò)推導(dǎo)公式得到疫情的基本再生數(shù)[20]。Jung等使用馬爾可夫鏈蒙特卡羅方法擬合指數(shù)增長(zhǎng)模型，然后通過(guò)擬合的參數(shù)計(jì)算基本再生數(shù)[21]。Sanche等通過(guò)估計(jì)指數(shù)增長(zhǎng)率，結(jié)合潛伏期和時(shí)間間隔來(lái)估計(jì)基本再生數(shù)[22]。Tang等使用疫情數(shù)據(jù)擬合代際時(shí)間間隔的Gamma分布，并將數(shù)據(jù)代入到更新方程來(lái)估計(jì)基本再生數(shù)[23]。武文韜等通過(guò)有效接觸率和感染期，對(duì)廣東省新冠肺炎疫情的基本再生數(shù)進(jìn)行了估計(jì)[24]。李盈科等通過(guò)群體的生育率以及矩量生成函數(shù)推導(dǎo)出基本再生數(shù)[25]。謝家榮等構(gòu)建了滾動(dòng)SEIR傳播動(dòng)力學(xué)模型，通過(guò)新增感染人數(shù)和累計(jì)感染人數(shù)表達(dá)式推導(dǎo)出傳染病的基本再生數(shù)[26]。

梳理上述文獻(xiàn)發(fā)現(xiàn)，現(xiàn)有主流文獻(xiàn)中的傳染病傳播動(dòng)力傳播模型中均沒(méi)有考慮疑似者的存在，還假設(shè)潛伏者沒(méi)有傳染性?，F(xiàn)實(shí)中由于技術(shù)的短缺，存在大量無(wú)法確診的疑似者，這些疑似者中有一部分人具有傳染性。同時(shí)，潛伏者也有傳染性。本文同時(shí)考慮了這些具體的因素而建立新的傳染病傳播動(dòng)力學(xué)模型。因本模型考慮的問(wèn)題更符合現(xiàn)實(shí)，使用本模型與傳統(tǒng)經(jīng)典模型相比，在同樣的情況下對(duì)疫情發(fā)展的預(yù)測(cè)會(huì)更準(zhǔn)確一些。其次，在現(xiàn)有經(jīng)典文獻(xiàn)中，主要從潛伏期、感染期以及指數(shù)增長(zhǎng)率這些因素和定義角度來(lái)計(jì)算傳染病的基本再生數(shù)，這些方法適用于比較簡(jiǎn)單的傳染病模型?，F(xiàn)有文獻(xiàn)幾乎沒(méi)有學(xué)者從無(wú)病平衡點(diǎn)的局部穩(wěn)定性的角度，采用基本再生矩陣推導(dǎo)更復(fù)雜情況下的傳染病基本再生數(shù)的表達(dá)式。因基本再生矩陣更適合復(fù)雜情況下的傳染病傳播模型，得到的基本再生數(shù)表達(dá)式也更能描述傳染病的傳播能力。最后，現(xiàn)有文獻(xiàn)主要是從最小二乘法、智能算法等方法來(lái)估計(jì)傳染病模型的參數(shù)，使用馬爾可夫鏈蒙特卡羅方法對(duì)傳染病模型參數(shù)進(jìn)行估計(jì)的學(xué)者較少。馬爾可夫鏈蒙特卡羅可以通過(guò)隨機(jī)變量的先驗(yàn)分布，結(jié)合數(shù)據(jù)，得到最接近真實(shí)數(shù)據(jù)的參數(shù)后驗(yàn)分布，但馬爾可夫鏈蒙特卡羅算法在執(zhí)行過(guò)程中，無(wú)法監(jiān)測(cè)其收斂性程度，因此，有必要設(shè)計(jì)一個(gè)統(tǒng)計(jì)指標(biāo)對(duì)馬氏鏈的收斂性進(jìn)行判別。故本文在這些方面做了一些探索性的工作。本文具體研究思路是，在人群管制模式下，以潛伏者和部分疑似者均具有傳染性為前提假設(shè)，構(gòu)建了傳染病傳播動(dòng)力學(xué)模型，并使用基本再生矩陣?yán)碚撏茖?dǎo)出疫情的有效再生數(shù)。同時(shí)設(shè)計(jì)改進(jìn)的馬爾可夫鏈蒙特卡羅算法估計(jì)模型的參數(shù)，并進(jìn)行馬氏鏈?zhǔn)諗吭\斷，然后進(jìn)行了仿真模擬，最后分析初始易感者數(shù)量和有效接觸率對(duì)疫情趨勢(shì)和有效再生數(shù)的影響。綜上，本文與上述經(jīng)典文獻(xiàn)相比，有以下3點(diǎn)創(chuàng)新：①根據(jù)新冠病毒的傳播特性，將潛伏者和部分疑似者均具有傳染性這兩個(gè)特征納入傳染病模型中，將SEIR模型拓展為SEDIQR模型，彌補(bǔ)了以前經(jīng)典文獻(xiàn)中模型的不足；②考慮新因素的前提下，使用基本再生矩陣?yán)碚撏茖?dǎo)疫情的有效再生數(shù)表達(dá)式，這種方法更能揭示疫情發(fā)展的內(nèi)在規(guī)律和影響疫情發(fā)展的關(guān)鍵因素；③將改進(jìn)的馬爾可夫鏈蒙特卡羅方法和Geweke方法相結(jié)合，這樣更能得到符合真實(shí)數(shù)據(jù)的參數(shù)后驗(yàn)分布，進(jìn)而對(duì)參數(shù)進(jìn)行估算，彌補(bǔ)了馬爾可夫鏈蒙特卡羅方法收斂性無(wú)法判別的缺點(diǎn)。

1 傳染病動(dòng)力學(xué)模型

1.1 模型描述

根據(jù)人群管制下疫情的現(xiàn)實(shí)傳播規(guī)律，將人群分為易感者、潛伏者、疑似者、隔離者、感染者、康復(fù)者6個(gè)類別。假設(shè)潛伏者在傳染期內(nèi)有傳染性，并考慮疑似者的情況，認(rèn)為部分疑似者也有傳染性。假設(shè)潛伏者和部分疑似者的傳染性與感染者的一樣強(qiáng)，首先感染者、疑似者和潛伏者以相同傳染性的易感者，將易感者轉(zhuǎn)化為成潛伏者，然后潛伏者以一定的比例成為未確診的疑似者和確診的感染者。由于采取了防控措施對(duì)疑似者和感染者進(jìn)行部分隔離，因此，疑似者和感染者中分別有一部分人群進(jìn)入到隔離者類別中，且隔離者沒(méi)有傳染性。隨著對(duì)疑似者和感染者的隔離治療，最終隔離者將會(huì)以一定的比例成為康復(fù)者和死亡者。人群管制下傳染病傳播動(dòng)力學(xué)模型示意圖如圖1所示。

S、E、D、I、Q、R分別表示易感者、潛伏者、疑似者、感染者、隔離者和康復(fù)者人數(shù)圖1 疫情的傳播特征

1.2 模型建立

設(shè)N為某區(qū)域內(nèi)的總?cè)藬?shù)，S(t)、E(t)、D(t)、I(t)、Q(t)、R(t)分別為t時(shí)刻(天)的易感者、潛伏者、疑似者、感染者、隔離者和康復(fù)者人數(shù)，β(t)表示隨時(shí)間變化的疫情有效接觸率，它是關(guān)于時(shí)間t的分段函數(shù)，α為疑似者中未確診感染者所占的比例，?為潛伏者轉(zhuǎn)為疑似者的比例，σ為潛伏者轉(zhuǎn)為感染者的比例，θ表示疑似者轉(zhuǎn)為隔離者的轉(zhuǎn)移率，η表示感染者轉(zhuǎn)為隔離者的轉(zhuǎn)移率，γ表示隔離者的退出率，δ表示疫情的平均死亡率。不考慮出生率、死亡率和人口遷移對(duì)疫情傳播的影響。那么人群管制下傳染病傳播動(dòng)力學(xué)模型：

(1)

(2)

(3)

(4)

(5)

(6)

首先定義疫情的有效接觸率函數(shù)，假設(shè)t=τ1之前，疫情的有效接觸率為一個(gè)常數(shù)，即β(t)=β0，在t=τ1之后到t=τ2之前，假設(shè)有效接觸率快速增加，為β(t)=β1，在t=τ2之后，政府意識(shí)到疫情的嚴(yán)重性，如果不加以控制，疫情就會(huì)蔓延，因此，對(duì)部分疑似者和感染者進(jìn)行隔離，此時(shí)人們之間的接觸減少，有效接觸率開(kāi)始下降。假設(shè)有效接觸率以速率κ進(jìn)行指數(shù)式減少，那么整個(gè)疫期內(nèi)的有效接觸率可以用如下的分段函數(shù)來(lái)表示，即

(7)

1.3 傳染病的有效再生數(shù)

在人群管制模式下，基本再生數(shù)是動(dòng)態(tài)變化的，這時(shí)可以用有效再生數(shù)來(lái)對(duì)疫情的傳播潛力進(jìn)行衡量。

結(jié)合人群管制下的傳染病傳播動(dòng)力學(xué)模型以及動(dòng)態(tài)有效接觸率函數(shù)，可以通過(guò)基本再生矩陣的方法求得人群管制條件下的有效再生數(shù)。

設(shè)x=(x1,x2,…,xn)t為每個(gè)人群類別在t時(shí)刻的個(gè)體數(shù)量，且有xi≥0，i=1,2,…,n，無(wú)病平衡點(diǎn)為Xs={x≥0|xi=0,i=1,2,…,m}，m表示有感染個(gè)體的人群類別數(shù)量。Fi(x)為人群類別i中新感染個(gè)體的比率，Vi+(x)表示以各種方式轉(zhuǎn)換到人群類別i的轉(zhuǎn)換率，Vi-(x)表示從人群類別i轉(zhuǎn)移到其他類別的轉(zhuǎn)換率，Vi(x)=Vi-(x)-Vi+(x)為轉(zhuǎn)移比率。那么傳染病傳播動(dòng)力學(xué)模型可以表示為

(8)

如果X0為模型式(8)的無(wú)病平衡解，那么對(duì)于導(dǎo)數(shù)dF(X0)和dV(X0)，有[27]

有效再生數(shù)可以通過(guò)傳染病的無(wú)病平衡點(diǎn)推導(dǎo)得出，即通過(guò)求解下一代矩陣的最大特征值就可以得到有效再生數(shù)的表達(dá)式。

R(t)=ρ(FV-1)

(9)

由式(1)～式(6)可知，存在感染者的人群類別有4個(gè)，分別為E、D、I和Q，因此m=4，那么F、V均為四階矩陣，對(duì)于疫情的無(wú)病平衡點(diǎn)有E=D=I=Q=0,設(shè)人群管制下疫情的無(wú)病平衡點(diǎn)為x0，則x0=(S0,0,0,0,0)，為了推導(dǎo)的方便，這里假設(shè)模型中的有效接觸率為β，那么有

同樣計(jì)算有

那么FV-1的特征值為

(10)

因此FV-1的譜半徑為

(11)

將有效接觸率β表示成式(7)所示的依時(shí)間變化的形式，以反映不斷變化的公共衛(wèi)生干預(yù)措施，那么人群管制下的疫情有效再生數(shù)可以表示為

(12)

2 傳染病模型參數(shù)估計(jì)

假設(shè)人群管制下潛伏者的傳染性與感染者的傳染性相同，潛伏者轉(zhuǎn)為感染者的比例σ、感染者轉(zhuǎn)為隔離的轉(zhuǎn)移率η和疫情的平均死亡率δ是已知的(見(jiàn)文獻(xiàn)[14])，在上述參數(shù)已知的情況下，合理地估計(jì)剩余參數(shù)β0、β1、k、α、?、θ和γ是對(duì)傳染病進(jìn)行預(yù)測(cè)分析的前提基礎(chǔ)。

2.1 Geweke收斂診斷

將馬爾可夫鏈分成若干段，在馬氏鏈前一部分和后一部分漸進(jìn)獨(dú)立的假設(shè)條件下，Geweke構(gòu)造的檢驗(yàn)統(tǒng)計(jì)量近似服從正態(tài)分布，比較序列前一部分和后一部分均值的差異，然后構(gòu)造Geweke檢驗(yàn)統(tǒng)計(jì)量，如果檢驗(yàn)統(tǒng)計(jì)量的值介于-2到2之間，認(rèn)為馬氏鏈?zhǔn)鞘諗康摹eweke檢驗(yàn)統(tǒng)計(jì)量為

(13)

式中：E(xs)和V(xs)表示序列開(kāi)始部分的均值和方差；E(xe)和V(xe)表示序列結(jié)尾部分的均值和方差。通過(guò)Python的pymc3包的arviz.geweke函數(shù)得到收斂診斷結(jié)果，如圖2所示。

圖2 傳染病模型參數(shù)馬氏鏈的Geweke得分

由圖2可知，對(duì)于人群管制下傳染病傳播動(dòng)力學(xué)模型中的被估參數(shù)，其馬氏鏈的20個(gè)序列段的Geweke得分均在-2到2之間波動(dòng)，且每個(gè)序列的Geweke得分絕對(duì)值差異也比較小，說(shuō)明序列的開(kāi)始部分和結(jié)尾部分參數(shù)的均值差異較小，也即各參數(shù)的馬氏鏈達(dá)到平穩(wěn)狀態(tài)，馬爾可夫鏈蒙特卡羅算法達(dá)到收斂。

2.2 傳染病模型參數(shù)估計(jì)結(jié)果

馬爾可夫鏈蒙特卡羅(MCMC)是一種將馬爾可夫鏈和蒙特卡羅方法相結(jié)合的一種無(wú)監(jiān)督機(jī)器學(xué)習(xí)算法，是以馬爾可夫鏈為概率模型的蒙特卡羅方法，馬爾可夫鏈蒙特卡羅方法構(gòu)造一條馬爾可夫鏈，使其平穩(wěn)分布為目標(biāo)分布，進(jìn)而得到平穩(wěn)分布的樣本，然后依據(jù)樣本進(jìn)行蒙特卡羅模擬，從而對(duì)隨機(jī)變量的估計(jì)量進(jìn)行近似數(shù)值計(jì)算。本小節(jié)采用Metropolis-Hastings采樣(簡(jiǎn)稱M-H采樣)。M-H采樣是在MCMC采樣的基礎(chǔ)上進(jìn)行了改進(jìn)，從而解決了MCMC采樣可能存在接受概率較低的問(wèn)題。在對(duì)接受概率進(jìn)行改進(jìn)之后，就可以采用M-H算法對(duì)所需樣本進(jìn)行采樣，采樣算法如下：

1)首先構(gòu)造一條馬爾可夫鏈，使平穩(wěn)分布為目標(biāo)分布p(x)，設(shè)置收斂步數(shù)為n1，迭代步數(shù)為n2。

2)從簡(jiǎn)單的概率分布函數(shù)中獲得初始樣本值x0。

3)從t=0到t=n2-1進(jìn)行迭代：

第1步，從條件概率分布(建議分布)Q(x|xt)中獲得樣本值x*。

第2步，從均勻分布U(0,1)獲得一個(gè)隨機(jī)樣本值u。

第4步，如果不滿足第3步，就拒絕xt+1=x*，則xt+1=xt。

因此樣本集合(xn1+1,xn1+2,…,xn2)即為所需樣本集。

參數(shù)β0、β1、k、α、?、θ和γ的先驗(yàn)分布采用正態(tài)分布，建議分布選用多元正態(tài)分布，通過(guò)Python軟件執(zhí)行M-H算法過(guò)程，算法迭代了90 000次，經(jīng)過(guò)了75 000次的燃燒期。程序在Python的jupyter notebook上運(yùn)行了60 min之后得到結(jié)果。各參數(shù)的迭代軌跡如圖3所示，各參數(shù)的結(jié)果見(jiàn)表1。

圖3 人群管制下傳染病模型參數(shù)的馬氏鏈軌跡

由Geweke收斂診斷結(jié)果和各參數(shù)的馬氏鏈軌跡圖都可以看出，在經(jīng)過(guò)75 000次的燃燒期之后，各參數(shù)的馬氏鏈均達(dá)到平穩(wěn)狀態(tài)，可以使用平穩(wěn)分布的樣本估算參數(shù)的后驗(yàn)均值。由表1可知，β0、β1、k、α、?、θ和γ的后驗(yàn)均值分別為β0=0.129 1，β1=0.767 4，k=0.165 5，α=0.212 3，?=0.481 6，θ=0.508 9，γ=0.037 0。β0、β1、k、α、?、θ和γ的95%HDI分別為[0.126,0.132]、[0.766,0.768]、[0.161,0.169]、[0.211,0.213]、[0.477,0.492]、[0.5,0.519]、[0.036,0.038]。HDI稱為最高密度區(qū)間，用于估計(jì)貝葉斯統(tǒng)計(jì)的可信區(qū)間，為最小寬度的貝葉斯可信區(qū)間(BCI)，一般可以表示統(tǒng)計(jì)的不確定性。

表1 模型參數(shù)估計(jì)結(jié)果

3 實(shí)證分析

以人群控制典型韓國(guó)的數(shù)據(jù)為例。在疫情開(kāi)始第18天(2020年2月9號(hào)，τ1=18)舉行禮拜儀式，第31名感染者參與到這次儀式中，導(dǎo)致大規(guī)模感染，感染者迅速增加，到第27天(2020年2月18號(hào)，τ2=27)，第31名感染者確診，在這之后政府開(kāi)始采取防控措施應(yīng)對(duì)疫情。

選用2020年1月22日—4月30日的韓國(guó)的累計(jì)病例數(shù)、累計(jì)治愈人數(shù)、累計(jì)死亡人數(shù)數(shù)據(jù)進(jìn)行實(shí)證分析。數(shù)據(jù)來(lái)源于今日頭條公布的疫情數(shù)據(jù)，根據(jù)公布的數(shù)據(jù)畫出韓國(guó)累計(jì)病例數(shù)、累計(jì)治愈人數(shù)、累計(jì)死亡人數(shù)、現(xiàn)有確診人數(shù)隨傳播時(shí)長(zhǎng)的變化趨勢(shì)如圖4所示，新增確診人數(shù)如圖5所示。

圖4 累計(jì)病例、現(xiàn)有確診、累計(jì)治愈和累計(jì)死亡人數(shù)

圖5 新增確診人數(shù)

由圖4可知，韓國(guó)在疫情開(kāi)始第27天(2月18日)之后，累計(jì)病例數(shù)和現(xiàn)有確診人數(shù)開(kāi)始迅速增加，在疫情持續(xù)大約50天，在3月13日現(xiàn)有確診人數(shù)達(dá)到峰值，大約為7 400人左右，之后現(xiàn)有確診人數(shù)開(kāi)始下降，在4月10日之后，累計(jì)確診病例超過(guò)10 000人，之后基本保持平穩(wěn)。在疫情開(kāi)始3月1日之后開(kāi)始有治愈患者出院，隨后治愈患者增加的幅度先增加后減小，原因是醫(yī)院患者先增加后逐漸減小。死亡者的變化幅度不是很大，基本維持在一條水平直線狀態(tài)，可見(jiàn)韓國(guó)政府的防控措施取得了很大的效果。

由圖5可知，在疫情開(kāi)始一個(gè)月之前，新增確診人數(shù)整體變化不大，在一個(gè)月之后，新增人數(shù)迅速增加，在3月1日左右，新增人數(shù)達(dá)到峰值，約為810人，之后新增人數(shù)整體開(kāi)始出現(xiàn)下滑，在4月10日之后，新增人數(shù)總體上保持平穩(wěn)的狀態(tài)。

將模型計(jì)算結(jié)果與實(shí)際值進(jìn)行擬合，得到模型預(yù)測(cè)的住院患者與現(xiàn)有確診人數(shù)的擬合圖，如圖6所示，并得到模型預(yù)測(cè)的感染者人數(shù)和康復(fù)者人數(shù)變化趨勢(shì)，如圖7所示。

圖6 確診人數(shù)模型值與實(shí)際值擬合結(jié)果

圖7 模型預(yù)測(cè)的感染者和康復(fù)者

由圖6可知，韓國(guó)預(yù)計(jì)在3月下旬，現(xiàn)有確診人數(shù)達(dá)到峰值，約為7 000人左右，韓國(guó)約在3月1日左右疫情開(kāi)始出現(xiàn)拐點(diǎn)，在6月初，隨著大量的住院患者被治愈，韓國(guó)的現(xiàn)有確診人數(shù)逐漸減小，接近于0。模型得出的現(xiàn)有確診人數(shù)的預(yù)測(cè)曲線與實(shí)際情況基本一致。由此說(shuō)明，使用人群管制下的傳染病傳播動(dòng)力學(xué)模型來(lái)預(yù)測(cè)韓國(guó)疫情，具有比較好的預(yù)測(cè)效果，從而驗(yàn)證了所建立的模型的有效性。

由圖7(a)可知，韓國(guó)的感染者在2月中旬之前并沒(méi)有明顯的變化，在2月中旬之后，感染者迅速增加，在3月1日左右，感染者人數(shù)到達(dá)峰值，超過(guò)1 000人，之后開(kāi)始下降，大約在4月10日，感染者人數(shù)減少到0。由圖7(b)可知，在2月下旬之前，還未有患者被治愈，在2月下旬之后逐漸有康復(fù)者出現(xiàn)，并且康復(fù)者增加的速率先增大后減小，6月初，康復(fù)者人數(shù)超過(guò)10 000人，之后進(jìn)入平穩(wěn)階段。

將估計(jì)的參數(shù)和已知參數(shù)代入到人群管制下的有效再生數(shù)和有效接觸率表達(dá)式中，可以得到R(t)和β(t)隨時(shí)間t的變化趨勢(shì)，如圖8所示。

圖8 疫情有效再生數(shù)和有效接觸率的變化趨勢(shì)

由圖8(a)可知，在疫情開(kāi)始20 d之前，人群管制下的有效再生數(shù)是一個(gè)恒定不變的常數(shù)。在20 d之后到27 d左右，有效再生數(shù)也是一個(gè)恒定不變的常數(shù)，但比其疫情開(kāi)始20 d之前的有效再生數(shù)高。在第1階段有效再生數(shù)接近于1，在第2階段有效再生數(shù)迅速增加，將近于6。此時(shí)政府應(yīng)該保持高度重視，如果不采取措施加以應(yīng)對(duì)，疫情有大范圍爆發(fā)的風(fēng)險(xiǎn)。從此之后，政府開(kāi)始意識(shí)到疫情的嚴(yán)重性，采取了疫情防控措施，有效再生數(shù)開(kāi)始下降。在疫情開(kāi)始大約37 d左右，有效再生數(shù)降到1水平以下。說(shuō)明政府的防控策略有了比較明顯的效果，疫情有逐漸緩解最后達(dá)到消亡的趨勢(shì)。

由圖8(b)在疫情開(kāi)始20 d之前，人群管制下的動(dòng)態(tài)有效接觸率是一個(gè)常數(shù)。在疫情開(kāi)始20 d到大約37 d左右，動(dòng)態(tài)有效接觸率有所上升，但也是一個(gè)常數(shù)，之后動(dòng)態(tài)有效接觸率開(kāi)始下降。

比較圖8(a)和圖8(b)以發(fā)現(xiàn)，人群管制下疫情的有效再生數(shù)和動(dòng)態(tài)有效接觸率有著密切的聯(lián)系，并且兩者之間的變化趨勢(shì)是一致的，而有效再生數(shù)能夠反映疫情的傳播能力。因此，如何控制減小動(dòng)態(tài)有效接觸率從而減小疫情的有效再生數(shù)，降低疫情的傳播能力和傳染風(fēng)險(xiǎn)，是政府管理部門和醫(yī)療專家值得思考的問(wèn)題。

4 疫情的管控力度分析

4.1 管控力度對(duì)疫情發(fā)展趨勢(shì)的影響

傳染病的傳播方式主要為接觸傳染，一般是感染者接觸易感者進(jìn)行傳染的。因此，只有阻斷人與人之間的傳播，實(shí)施一系列公共衛(wèi)生干預(yù)措施，才能真正有效地控制和遏制新冠病毒的爆發(fā)。初始易感者是指在某個(gè)疫情開(kāi)始的節(jié)點(diǎn)，在一個(gè)傳染病系統(tǒng)中所存在的易感者數(shù)量，其數(shù)量的多少對(duì)于疫情的發(fā)展具有重要的影響。接觸數(shù)是衡量接觸程度的重要指標(biāo)，能夠從側(cè)面反映傳播風(fēng)險(xiǎn)的大小，減小這兩個(gè)指標(biāo)的值對(duì)于減小傳播率、降低傳播風(fēng)險(xiǎn)進(jìn)而降低疫情的蔓延勢(shì)頭具有重要的意義。在傳播概率不變的情況下，接觸數(shù)與有效接觸率成正比。不同的初始易感者數(shù)量S0和有效接觸率β1對(duì)于確診人數(shù)和感染者發(fā)展趨勢(shì)的影響程度如圖9和圖10所示。易感者的管控力度對(duì)疫情峰值和峰值時(shí)間的影響見(jiàn)表2，有效接觸率的管控力度對(duì)疫情峰值和峰值時(shí)間的影響見(jiàn)表3。

圖9 不同系數(shù)的S0對(duì)韓國(guó)現(xiàn)有確診人數(shù)和感染者的影響

圖10 不同系數(shù)的β1對(duì)韓國(guó)現(xiàn)有確診人數(shù)和感染者的影響

表2 易感者的管控力度對(duì)疫情峰值和峰值時(shí)間的影響

表3 有效接觸率的管控力度對(duì)疫情峰值和峰值時(shí)間的影響

1)由圖9(a)和表2可知，若將初始易感者數(shù)量擴(kuò)大為原有初始易感者數(shù)量的1.1倍時(shí)，現(xiàn)有確診人數(shù)在前期增加的幅度最大，且峰值最高，超過(guò)了25 000人。在疫情開(kāi)始約51 d到達(dá)峰值，峰值時(shí)間是最晚的，且現(xiàn)有確診人數(shù)達(dá)到峰值后下降的速率也最快。若將初始易感者的數(shù)量縮小為原有初始易感者數(shù)量的0.7倍時(shí)，其現(xiàn)有確診人數(shù)的峰值是最小的，只有100多人。在疫情的第48天到達(dá)峰值，峰值時(shí)間是最早的。且疫情前期現(xiàn)有確診人數(shù)的增加幅度最慢，到達(dá)峰值后下降的速率也最慢。由此可以發(fā)現(xiàn)這樣的規(guī)律，隨著對(duì)易感者管控力度的加強(qiáng)，確診人數(shù)峰值時(shí)間是越來(lái)越早的，峰值人數(shù)是越來(lái)越小的。前期確診人數(shù)增加的幅度越來(lái)越小，到達(dá)峰值后人數(shù)下降的幅度也越來(lái)越小。從圖9的曲線形狀來(lái)看，初始易感者數(shù)量對(duì)于確診人數(shù)的影響是非常敏感的，較小的變動(dòng)對(duì)于其人數(shù)的變化是非常大的。

由圖9(b)和表2可知，若將初始易感者數(shù)量擴(kuò)大為原有初始易感者數(shù)量的1.1倍時(shí)，感染者的人數(shù)在峰值之前的增長(zhǎng)幅度最快，峰值也最高，將近2 000人。在疫情第40天到達(dá)峰值，并且在到達(dá)峰值之后的感染者下降速率也最快。若將初始易感者數(shù)量縮小為原來(lái)的0.7倍時(shí)，感染者在峰值之前和峰值之后的增長(zhǎng)和下降速率都是最慢的，且在疫情第37天到達(dá)峰值。并且隨著對(duì)易感者管控力度的加強(qiáng)，感染者峰值時(shí)間越來(lái)越早，峰值人數(shù)也越來(lái)越小。初始易感者的數(shù)量對(duì)于感染者的影響也是非常敏感的。

2)由圖10(a)和表3可知，若將有效接觸率擴(kuò)大為初始有效接觸率的1.1倍時(shí)，現(xiàn)有確診人數(shù)在前期增加的幅度最大，且峰值最高，多于20 000人，峰值時(shí)間也是最晚的。且確診人數(shù)達(dá)到峰值后下降的速率也最快，但疫情結(jié)束的時(shí)間晚于其他接觸數(shù)下的現(xiàn)有確診人數(shù)曲線。若將有效接觸率數(shù)量縮小為原來(lái)的0.7倍時(shí)，其確診人數(shù)的峰值是最小的，約為300人左右，峰值時(shí)間是最早的。且疫情前期確診人數(shù)的增加幅度最慢，到達(dá)峰值后下降的速率也最慢。由此也可以發(fā)現(xiàn)這樣的規(guī)律，隨著對(duì)有效接觸率管控力度的加強(qiáng)，現(xiàn)有確診人數(shù)峰值時(shí)間是越來(lái)越早的，峰值人數(shù)是越來(lái)越小的，前期現(xiàn)有確診人數(shù)增加的幅度越來(lái)越小，到達(dá)峰值后人數(shù)下降的幅度也越來(lái)越小。

由圖10(b)和表3可知，若將有效接觸率擴(kuò)大為原來(lái)的1.1倍時(shí)，感染者的人數(shù)在峰值之前和峰值之后的增長(zhǎng)和下降速率都是最快的，且峰值最高，將近3 000人，在疫情開(kāi)始約51 d達(dá)到峰值。若將有效接觸率縮小為原來(lái)的0.7倍時(shí)，感染者的人數(shù)在峰值之前和峰值之后的增長(zhǎng)和下降速率都是最慢的，且峰值最小，不到100人，約在疫情開(kāi)始37 d到達(dá)峰值。因此，隨著對(duì)有效接觸率管控力度的加強(qiáng)，感染者的峰值越來(lái)越小，峰值達(dá)到時(shí)間越來(lái)越早。

3)結(jié)合圖9和圖10，表2和表3可以發(fā)現(xiàn)，當(dāng)對(duì)初始易感者和有效接觸率實(shí)施相同的管控力度時(shí)，不影響現(xiàn)有確診人數(shù)和感染者的峰值時(shí)間，但是對(duì)峰值大小卻有比較明顯的影響。當(dāng)初始易感者數(shù)量和有效接觸率呈現(xiàn)同幅度下降時(shí)，初始易感者數(shù)量對(duì)于緩解疫情發(fā)展具有更大的作用。

4.2 管控力度對(duì)有效再生數(shù)的影響

有效再生數(shù)是衡量疫情傳播能力的重要閾值參數(shù)，有效再生數(shù)大于1說(shuō)明疫情有進(jìn)一步蔓延的趨勢(shì)，有效再生數(shù)小于1說(shuō)明疫情有所緩解直到消亡。從人群管制下有效再生數(shù)的表達(dá)式可以看出，初始易感者數(shù)量和有效接觸率對(duì)有效再生數(shù)有很大的影響。不同的初始易感者數(shù)量和有效接觸率對(duì)有效再生數(shù)的影響程度如圖11和圖12所示。

圖11 不同系數(shù)的S0對(duì)有效再生數(shù)的影響

圖12 不同系數(shù)的β1對(duì)有效再生數(shù)的影響

1)由圖11可知，若將初始易感者的數(shù)量擴(kuò)大為原有初始易感者數(shù)量的1.3倍時(shí)，在疫情開(kāi)始到第18天，有效再生數(shù)是不到2的常數(shù)。在疫情開(kāi)始第18天到27天，有效再生數(shù)繼續(xù)增加，說(shuō)明在此時(shí)疫情開(kāi)始蔓延，之后由于政府采取了疫情防控措施，有效再生數(shù)開(kāi)始以指數(shù)形式下降，疫情蔓延的速度有所下降。在疫情開(kāi)始第39天之后，有效再生數(shù)處于1水平以下，說(shuō)明此時(shí)疫情開(kāi)始有所緩解并逐漸消亡。若將初始易感者的數(shù)量縮小為原有初始易感者數(shù)量的0.2倍時(shí)，在疫情開(kāi)始到18天，有效再生數(shù)小于1，此時(shí)說(shuō)明疫情并沒(méi)有流行開(kāi)來(lái)。在疫情開(kāi)始第18天到27天，有效再生數(shù)有所上升，并超過(guò)了1，之后開(kāi)始下降，在疫情開(kāi)始28天之后，有效再生數(shù)小于1。由此可以發(fā)現(xiàn)，隨著初始易感者數(shù)量占原來(lái)初始易感者數(shù)量比例的減少，有效再生數(shù)變化趨勢(shì)并未改變，但在[0,τ1]和[τ1,τ2]這兩個(gè)時(shí)間段的有效再生數(shù)隨之下降，有效再生數(shù)下降為1所需的時(shí)間也越來(lái)越短。

2)由圖12可知，若將接觸數(shù)擴(kuò)大為原有接觸數(shù)的1.3倍時(shí)，在疫情開(kāi)始到第18天，有效再生數(shù)保持在1左右。在疫情開(kāi)始第18天到第27天，有效再生數(shù)超過(guò)7，之后有效再生數(shù)開(kāi)始下降。在疫情開(kāi)始第39天之后，有效再生數(shù)下降到1以下。若將接觸數(shù)縮小為原有接觸數(shù)的0.2倍時(shí)，在疫情開(kāi)始到第18天，有效再生數(shù)接近于1，在疫情開(kāi)始第18天到第27天，有效再生數(shù)接近于1.2，此時(shí)疫情開(kāi)始蔓延。在疫情開(kāi)始第28天之后，有效再生數(shù)處于1以下。由圖形趨勢(shì)可以發(fā)現(xiàn)，隨著有效接觸率占原來(lái)有效接觸率比例的減少，在[0,τ1]和[τ1,τ2]時(shí)間段的有效再生數(shù)保持不變。但在τ2之后有效再生數(shù)開(kāi)始下降，并且有效再生數(shù)下降到1時(shí)所需的時(shí)間也越來(lái)越短。

3)比較圖11和圖12可知，對(duì)初始易感者數(shù)量和有效接觸率實(shí)行同樣的管控力度時(shí)，不會(huì)影響τ1時(shí)間之后有效再生數(shù)的變化，也不會(huì)影響有效再生數(shù)降為1所需的時(shí)間，但對(duì)[0,τ1]時(shí)間段的有效再生數(shù)影響卻不同。

5 結(jié)論

本文綜合考慮人群隔離、潛伏者和部分疑似者均具有傳染性的因素，構(gòu)建了人群管制模式的傳染病傳播動(dòng)力學(xué)模型，以此來(lái)預(yù)測(cè)疫情的發(fā)展趨勢(shì)和不同的管控力度對(duì)疫情發(fā)展和有效再生數(shù)的影響程度。通過(guò)有效再生數(shù)的推導(dǎo)、模型參數(shù)估計(jì)和實(shí)證分析得到以下幾點(diǎn)結(jié)論：

1)結(jié)合人群管制的特點(diǎn)，加入隔離者、考慮潛伏者和部分疑似者均具有傳染性這3個(gè)具體因素而建立的傳染病動(dòng)力學(xué)模型，能夠很好地描述人群管制下疫情的傳播規(guī)律，所建立的模型是有效且合理的。因建立的模型更符合人群管制模式下的疫情傳播，得到的預(yù)測(cè)結(jié)果相對(duì)于文獻(xiàn)[14]對(duì)韓國(guó)的疫情預(yù)測(cè)會(huì)更準(zhǔn)確一些。

2)基本再生矩陣?yán)碚撘罁?jù)無(wú)病平衡點(diǎn)的局部穩(wěn)定性，能夠很好地推導(dǎo)出人群管制下的有效再生數(shù)表達(dá)式，得出的有效再生數(shù)變化曲線也反映了人群管制下的疫情傳播能力，說(shuō)明相對(duì)于文獻(xiàn)[21-23,25]使用指數(shù)增長(zhǎng)率、感染期和潛伏期來(lái)計(jì)算基本再生數(shù)，基本再生矩陣更能夠科學(xué)合理地推導(dǎo)出傳染病的臨界閾值。

3)馬爾可夫鏈蒙特卡羅方法是將馬爾可夫鏈與蒙特卡羅相結(jié)合的一種無(wú)監(jiān)督機(jī)器學(xué)習(xí)算法，是在貝葉斯理論框架下通過(guò)計(jì)算機(jī)進(jìn)行動(dòng)態(tài)模擬的蒙特卡羅方法，主要應(yīng)用于求概率模型中未知參數(shù)的后驗(yàn)分布，從后驗(yàn)分布中采樣以構(gòu)造最接近真實(shí)數(shù)據(jù)的概率分布的一種方法。馬爾可夫鏈蒙特卡羅方法相對(duì)于最小二乘等優(yōu)化方法，能夠有效地尋找到真實(shí)的參數(shù)后驗(yàn)，其不失為一種理解性高，邏輯性強(qiáng)的參數(shù)估計(jì)方法。

4)從管控力度對(duì)疫情發(fā)展影響的分析結(jié)果可知，減小初始易感者數(shù)量和有效接觸率均能夠降低疫情現(xiàn)有確診人數(shù)和感染者的峰值，且能夠縮短到達(dá)峰值所需的時(shí)間。有效接觸率和初始易感者數(shù)量都是現(xiàn)有確診人數(shù)和感染者人數(shù)變化的重要敏感參數(shù)，兩者相輔相成。減少初始易感者的數(shù)量從另一方面也可以減少有效接觸率，減少有效接觸率也可以通過(guò)減少初始易感者的數(shù)量來(lái)實(shí)現(xiàn)，它們對(duì)疫情的發(fā)展都有重要的影響。當(dāng)初始易感者數(shù)量和有效接觸率呈現(xiàn)同幅度下降時(shí)，初始易感者數(shù)量對(duì)于緩解疫情發(fā)展具有更大的作用。從管控力度對(duì)疫情有效再生數(shù)的影響結(jié)果可知，減小初始易感者數(shù)量和有效接觸率均能夠縮短有效再生數(shù)降為1所需的時(shí)間，說(shuō)明減小初始易感者數(shù)量和有效接觸率能夠有效地緩解疫情的蔓延。從各國(guó)現(xiàn)實(shí)疫情管控效果來(lái)看，東方比西方管控得要好。原因就在于東方國(guó)家(韓國(guó)是個(gè)典型)比西方國(guó)家(美國(guó)是個(gè)典型)更重視對(duì)疫情初期對(duì)有效接觸率和初始易感者數(shù)量的控制。從某種角度來(lái)說(shuō)，美國(guó)在擁有全世界最強(qiáng)大的醫(yī)療設(shè)施條件下，對(duì)疫情管控的效果卻成為最差的國(guó)家之一，與其初期不按科學(xué)原則管控疫情有關(guān)。他們?yōu)槠渌麌?guó)家管控類似的疫情提供了反面教材，這也是世界各國(guó)應(yīng)該汲取教訓(xùn)的地方。本研究為今后在世界范圍內(nèi)控制類似疫情提供了理論和事實(shí)依據(jù)。

本文在人群管制條件下，將潛伏者和部分疑似者均具有傳染性這一特點(diǎn)考慮進(jìn)模型，比其他的學(xué)者考慮問(wèn)題更全面，更接近客觀現(xiàn)實(shí)，因此，預(yù)測(cè)的效果與現(xiàn)實(shí)結(jié)果相比較，擬合的效果更好。但沒(méi)有考慮新冠疫情中一些無(wú)癥狀感染者的情況，所以預(yù)測(cè)結(jié)果與現(xiàn)實(shí)還有一定的出入。另外，假設(shè)潛伏者和部分疑似者傳染率與感染者的傳染率相同，這一假設(shè)并沒(méi)有現(xiàn)實(shí)的數(shù)據(jù)支撐，將在后期進(jìn)一步對(duì)此展開(kāi)研究，并將新的研究完善到模型中去。