程萬港，陳萬銘，陳燁麗，夏志樂

（臺(tái)州學(xué)院 電子與信息工程學(xué)院，浙江 臨海 317000）

0 引言

新型冠狀病毒肺炎（以下簡稱新冠肺炎）是由2019新型冠狀病毒所感染導(dǎo)致的肺炎，經(jīng)常伴有發(fā)熱、干咳、乏力等癥狀。因此，它的臨床癥狀與普通肺炎以及普通感冒極其相似。但新冠肺炎與普通流感不同的是，其對(duì)人體有極大的危害性。由于新冠肺炎起病隱匿，潛伏時(shí)間是時(shí)變不確定的，一般為1～14天。早期容易與感冒以及普通肺炎混淆，極易導(dǎo)致漏診或誤診乃至死亡，而且它具有極強(qiáng)的傳播性。自疫情暴發(fā)以來，截至2020年7月24日，全球累計(jì)確診新冠肺炎病例15，641，288例，治愈926，505例，累計(jì)死亡635，598例。目前，雖然中國疫情得到了很大程度的控制，但國外感染人數(shù)仍在持續(xù)增加，嚴(yán)重影響人類的健康和世界經(jīng)濟(jì)的發(fā)展。因此對(duì)新冠病毒的傳染機(jī)理和調(diào)控措施的研究顯得尤為迫切。

目前已有眾多的學(xué)者致力于新冠肺炎疫情預(yù)測及防控的研究，建立了相應(yīng)的新冠肺炎數(shù)學(xué)模型，為新冠肺炎的調(diào)控起到了積極的作用。Qun等［1］搜集了肺炎人口統(tǒng)計(jì)學(xué)特征等信息，估計(jì)了疫情倍增時(shí)間和基本生殖數(shù)；Wan等［2］研究了武漢封城數(shù)據(jù)，基于SEIR建模方法，預(yù)測了武漢市疫情感染人數(shù)峰值；Joseph等［3］預(yù)測了新冠肺炎在國內(nèi)以及全世界范圍內(nèi)的公共衛(wèi)生風(fēng)險(xiǎn)和藥物預(yù)防干預(yù)措施；Smirnova等［4］比較了正則化算法，修改了TSVD，用有限數(shù)據(jù)進(jìn)行預(yù)測，確定最有效的穩(wěn)定策略；Ivanov等［5］研究了兩類SEIR模型，一類由常微分方程組描述，另一類以一階差分方程的離散模型描述，預(yù)測了疫情的高峰日和最大感染人數(shù)；Bentout等［6］運(yùn)用SEIR模型對(duì)阿爾及利亞進(jìn)行預(yù)測，利用最小二乘法和最小化殘差平方和的方法擬合數(shù)據(jù)，估計(jì)了流行病相關(guān)參數(shù)和基本繁殖數(shù)；蔡潔等［7］考慮了新冠肺炎疫情期間城市的管控力度，基于SEIR模型模擬出武漢市疫情發(fā)展情況；許家雪等［8］在傳統(tǒng)SEIR傳染病模型基礎(chǔ)上，結(jié)合實(shí)際情況建立疫情傳播的微分方程，并對(duì)疫情的傳播規(guī)模進(jìn)行了預(yù)測。以上文獻(xiàn)表明，SEIR模型能夠清晰地描述出病毒傳播的邏輯關(guān)系，并能對(duì)疫情的發(fā)展趨勢(shì)進(jìn)行較為準(zhǔn)確的預(yù)測。但是傳統(tǒng)的SEIR模型沒有綜合考慮疫情地區(qū)人口流動(dòng)、感染人員可能死亡、潛伏時(shí)滯可能是時(shí)變的等重要的影響因素。

綜上所述，本文對(duì)傳統(tǒng)的SEIR模型進(jìn)行了改進(jìn)，給出了SEIR非線性時(shí)變時(shí)滯模型，所提模型綜合考慮了疫情區(qū)域人員流動(dòng)、感染者可能死亡、潛伏時(shí)變時(shí)滯等不確定因素。因此，可以更精確有效地處理現(xiàn)實(shí)問題。對(duì)于疫情區(qū)域人員流動(dòng)以及感染人員死亡等影響因素，本文采用混合型高斯混合模型進(jìn)行建模。根據(jù)Expectation-Maximization（EM）算法，確定模型中未知參數(shù)的最優(yōu)取值。最后，以中國湖北省疫情數(shù)據(jù)對(duì)改進(jìn)后模型進(jìn)行仿真，且與元胞自動(dòng)機(jī)模型進(jìn)行了比較。結(jié)果表明，本文所提供的模型在運(yùn)行時(shí)間和預(yù)測準(zhǔn)確性等方面明顯優(yōu)于傳統(tǒng)SEIR模型和元胞自動(dòng)機(jī)模型［9］-［11］。

1 前言

為了便于后面進(jìn)一步分析，本部分對(duì)傳統(tǒng)的SEIR傳染病模型進(jìn)行了簡單的描述，然后給出基本假設(shè)。

1.1 傳統(tǒng)SEIR模型的描述

傳統(tǒng)的SEIR傳染病模型如下：

其中，S(t)，E(t)，I(t)和R(t)分別表示易感者，潛伏者，感染者和康復(fù)者在t時(shí)刻的人數(shù)；q是常數(shù)，表示潛伏期的時(shí)間；α表示潛伏者轉(zhuǎn)為感染者的概率；γ表示感染者康復(fù)的概率；r1表示每個(gè)感染者平均每天接觸的人數(shù)；r2表示每個(gè)潛伏者平均每天接觸人數(shù)；β1表示患者的傳染概率；β2表示潛伏者傳染概率；β3表示新冠肺炎患者轉(zhuǎn)陰率;N表示疫情地區(qū)總?cè)藬?shù)。

模型（1）反映的傳染病傳播路徑的示意圖如圖1表示：

圖1 傳統(tǒng)的SEIR傳染病模型的示意圖

注1：從模型（1）和圖1可以看出，傳統(tǒng)的SEIR傳染病模型沒有考慮到疫情地區(qū)人員流動(dòng)和死亡等因素的影響。實(shí)際中，疫情地區(qū)人員的流動(dòng)以及受疫情感染人員的死亡對(duì)疫情的防控以及疫情的發(fā)展預(yù)測有著極大的影響。為了使傳統(tǒng)模型（1）能更好的解決實(shí)際問題，本文在原有模型的基礎(chǔ)上增加了人員流動(dòng)和死亡兩個(gè)重要的影響因素。

注2：模型（1）中考慮的潛伏時(shí)間僅為常數(shù)，而在實(shí)際中，潛伏時(shí)間可能因人而異，甚至是未知時(shí)變的。本文將研究將潛伏時(shí)滯為常數(shù)推廣到時(shí)變的情況。

1.2 基本假設(shè)

假設(shè)1：假設(shè)新型冠狀肺炎只存在人傳人的現(xiàn)象。

假設(shè)2：不考慮外界環(huán)境因素對(duì)病死率參數(shù)的影響。

假設(shè)3：假設(shè)感染者康復(fù)后自身攜帶抗體，不會(huì)再出現(xiàn)二次感染的情況。

假設(shè)4：假設(shè)未被發(fā)現(xiàn)的病毒攜帶者和其他人接觸的機(jī)會(huì)相等。

注3：因截至目前，只發(fā)現(xiàn)人傳人現(xiàn)象，所以假設(shè)1是合理的；因目前還未發(fā)現(xiàn)外界環(huán)境因素如空氣質(zhì)量、溫度、濕度等對(duì)疫情的傳播及感染有重要影響，所以假設(shè)2是合理的；因至今還未出現(xiàn)感染者康復(fù)后二次感染的情況，所以假設(shè)3是合理的；在未加有效控制的情況下，人們生活正常，人與人接觸機(jī)會(huì)相等，因此假設(shè)4是合理的。

2 模型建立

2.1 對(duì)于傳統(tǒng)模型的改進(jìn)

本部分內(nèi)容我們對(duì)模型（1）進(jìn)行改進(jìn)。為了更加深入地研究各種因素對(duì)新冠病毒感染性的影響，在傳統(tǒng)的SEIR模型的基礎(chǔ)上，我們?cè)黾恿顺鞘腥丝谶w入與遷出等變量。因?yàn)閷?duì)疫情采取了防控隔離措施，并且康復(fù)者具有抗體，不考慮再次被感染的情況，故本文在傳統(tǒng)的SEIR模型中增加了康復(fù)者、外界人口、死亡人員這三個(gè)群體。

改進(jìn)后的SEIR傳染病模型為：

其中，k表示死亡率；d(t)為潛伏時(shí)變時(shí)滯；疫情地區(qū)人員流動(dòng)情況用以下的高斯混合模型來描述：

其中，ai,bi,ci為可調(diào)參數(shù)，在文中后面部分將會(huì)給出確定其最優(yōu)值的方法。

模型（2）的示意圖如圖2所示：

圖2 改進(jìn)的SEIR傳染病模型

其中β5表示遷入率，β4表示遷出率，它們的差值β4-β5可用高斯型函數(shù)（3）進(jìn)行模擬。

注4：與模型（1）相比，模型（2）加入感染人員死亡以及外界人口遷入遷出等因素，當(dāng)參數(shù)β4、β5以及k極小時(shí)，模型（2）退化為模型（1），因此模型（1）是模型（2）的一個(gè)特例。從數(shù)學(xué)角度看，f(t)是非線性函數(shù)，所以模型（1）的簡單線性疊加不可能獲得模型（2）。

注5：f(t)采用的函數(shù)為高斯型混合函數(shù)，其中部分參數(shù)的最優(yōu)解可以通過3.2給出的最大期望算法來獲得。

2.2 參數(shù)的優(yōu)化與選取

高斯混合模型（Gaussian Mixed Model）指的是多個(gè)高斯分布函數(shù)的線性組合。主要應(yīng)用于聚類分析以及概率密度估計(jì)上，具有對(duì)離散型數(shù)據(jù)擬合精度高且適用性強(qiáng)的優(yōu)點(diǎn)。適當(dāng)?shù)倪x取GMM中的可調(diào)參數(shù)，可以很好地描述疫情地區(qū)人員流動(dòng)的情況。下面給出GMM中可調(diào)參數(shù)的具體優(yōu)化方法。

GMM模型的EM優(yōu)化原理［11］：GMM（Gaussian Mixture Model）高斯混合模型由K個(gè)Gaussian分布組成，每個(gè)Gaussian稱為一個(gè)“高斯分量”，這些高斯分量的線性組合就組成了GMM的概率密度函數(shù)：

其中πk表示每個(gè)高斯分量被選中的概率（權(quán)值），μk表示均值，Σk表示方差。

GMM的對(duì)數(shù)似然函數(shù)為：

由于對(duì)數(shù)函數(shù)里面有加權(quán)，無法用直接求導(dǎo)的方法求解該方程的求最大值。為了解決這個(gè)問題，我們采取EM算法進(jìn)行優(yōu)化參數(shù)，下面給出算法的具體步驟：

步驟1：估計(jì)數(shù)據(jù)由每個(gè)高斯分量生成的概率（并不是每個(gè)高斯分量被選中的概率）;對(duì)于每個(gè)數(shù)據(jù)來說，它的第k個(gè)高斯分量生成的概率為

步驟2：估計(jì)每個(gè)高斯分量的參數(shù)：假設(shè)上一步得到的 是正確的“數(shù)據(jù)x由高斯分量k生成的概率”，亦可當(dāng)作該高斯分量在生成這個(gè)數(shù)據(jù)上所做的貢獻(xiàn)。集中考慮所有的數(shù)據(jù)點(diǎn)，實(shí)際上可以看作高斯分量生成了γ(1,k)x1,……,γ(N,k)xN這些點(diǎn)。由于每個(gè)高斯分量都是一個(gè)標(biāo)準(zhǔn)的Gaussian分布，可以容易求出最大似然所對(duì)應(yīng)的參數(shù)值如下：

步驟3：重復(fù)迭代前兩步，直到似然函數(shù)的值收斂為止。

算法檢驗(yàn)：采用以上算法得出GMM中項(xiàng)數(shù)n選為3是效果最好，即：

其中利用MATLAB編程得到的部分?jǐn)?shù)據(jù)如表1所示：

表1 各個(gè)參數(shù)的選取以及其置信區(qū)間

結(jié)果顯示錯(cuò)誤平方和為SSE=9.117，多重測定系數(shù)R-square=0.8984，表明該函數(shù)對(duì)此數(shù)據(jù)具有較好的擬合效果，擬合圖如圖3所示：

圖3 函數(shù)擬合圖

對(duì)于易感者引入了人口流動(dòng)以β3-β4表示人口凈遷入率，潛伏周期q=d(t)，對(duì)模型（2）進(jìn)行優(yōu)化得：

感染者以k的概率死亡，引入死亡者人數(shù)的變化方程：

由上可知，改進(jìn)后的SEIR傳染病模型為：

注6：如果所給數(shù)據(jù)不同，參數(shù)優(yōu)化的結(jié)果也會(huì)不同。

3 仿真研究

3.1 改進(jìn)模型數(shù)據(jù)仿真

為了說明所提模型的有效性，以下以中國湖北省為例。根據(jù)相關(guān)文獻(xiàn)以及中國衛(wèi)健委發(fā)布的每日湖北省新冠肺炎新增人數(shù)、治愈人數(shù)以及湖北省人口遷移等相關(guān)數(shù)據(jù)，可選擇參數(shù)如下：β1=0.6，β2=0.125，β3=0.11，α=0.07，r1=10，r2=20，k=0.037，γ=0.0601.

部分統(tǒng)計(jì)數(shù)據(jù)如表2所示：

表2 湖北省新冠肺炎2020年1月20日—2月19日疫情統(tǒng)計(jì)數(shù)據(jù)

將以上參數(shù)代入模型（2）中，得到的結(jié)果如圖4所示，其中潛伏時(shí)間q=7，由參數(shù)r1r2反映隔離強(qiáng)度，取不同的r1r2所得到的結(jié)果如圖4和圖5所示，結(jié)果可知不同隔離強(qiáng)度下疫情得以控制的程度不同，在較高水平隔離措施下，患者及潛伏者人數(shù)遠(yuǎn)遠(yuǎn)小于在較低水平隔離措施下患者及潛伏者的人數(shù)。因此，采取封城、限制居民出去、各類隔離手段是相當(dāng)有必要的，可以很大程度上緩解疫情往壞的方向發(fā)展的趨勢(shì)。在采取高水平隔離措施的條件下，疫情暴發(fā)后15天左右新冠肺炎感染者人數(shù)達(dá)到一個(gè)峰值，約為8.5萬人，疫情暴發(fā)后四個(gè)月左右基本穩(wěn)定，得以控制；而在較低水平隔離措施條件下，感染人數(shù)的峰值將達(dá)到13.8萬人。考慮時(shí)變時(shí)滯，這里取d(t)=7+(-1)t，在較高水平隔離措施的前提下，它的峰值將達(dá)到7.5萬人，結(jié)果如圖6所示。截至目前，湖北的累計(jì)感染人數(shù)為6.8萬人，說明考慮時(shí)變時(shí)滯時(shí)，本文所用的方法同樣適用，并且更符合實(shí)際，如果能給出更為精確的時(shí)變時(shí)滯，那么預(yù)測結(jié)果將會(huì)更加準(zhǔn)確。

圖4 r1=10，r2=20時(shí)各類人群人數(shù)預(yù)測圖

圖5 r1=20，r2=30時(shí)各類人群人數(shù)預(yù)測圖

圖6 加入時(shí)變時(shí)滯后各類人群人數(shù)預(yù)測圖

3.2 與元胞自動(dòng)機(jī)模型比較

用文獻(xiàn)［9］［12］提供的元胞自動(dòng)機(jī)模型進(jìn)行仿真，參數(shù)選取與3.1所述相同，可得結(jié)果如圖7所示：

圖7 湖北省新冠肺炎疫情元胞自動(dòng)機(jī)空間模擬結(jié)果示意圖及四種人群隨時(shí)間變化結(jié)果圖

從圖6和圖7數(shù)據(jù)仿真結(jié)果與疫情實(shí)際發(fā)展數(shù)據(jù)對(duì)比來看，文中所提供的方法，在模型的精確度以及使用范圍等方面明顯優(yōu)于傳統(tǒng)的SEIR模型以及元胞自動(dòng)機(jī)模型［13-14］。改進(jìn)后的SEIR傳染病模型具有運(yùn)行時(shí)間短，參數(shù)易獲取，模型簡單易行且預(yù)測結(jié)果更精確等特點(diǎn)。當(dāng)考慮人口流動(dòng)以及隔離狀態(tài)不穩(wěn)定時(shí)，元胞自動(dòng)機(jī)模型無效或不再適用，而本文方法仍能適用。

4 結(jié)論

本文就傳統(tǒng)的SEIR模型進(jìn)行改進(jìn)，以湖北省疫情發(fā)展的參數(shù)對(duì)其進(jìn)行仿真，并將得到的結(jié)果與元胞自動(dòng)機(jī)模型仿真出來的結(jié)果進(jìn)行比較?？梢园l(fā)現(xiàn)，改進(jìn)后的SEIR模型可以更好地預(yù)測新冠肺炎疫情發(fā)展趨勢(shì)，且可通過對(duì)參數(shù)的控制提出疫情防控方面的舉措，對(duì)做好精準(zhǔn)防控具有很好的指導(dǎo)作用，可以將改進(jìn)后的新的新型冠狀肺炎SEIR傳染病模型推廣到其他仍存在疫情的國家，幫助他們更有效地做好疫情防控工作，這也是我們后續(xù)研究的方向。