基于貝葉斯優(yōu)化SEIR模型的新型冠狀病毒預(yù)測與分析

王 慶,陳俊熹,熊東平,何嘯峰

(1.南華大學(xué) 計算機學(xué)院，湖南 衡陽 421001;2.南華大學(xué) 附屬南華醫(yī)院，湖南 衡陽 421000)

新型冠狀病毒肺炎(Corona Virus Disease 2019，COVID-19)，簡稱“新冠肺炎”，是指2019年新型冠狀病毒感染導(dǎo)致的肺炎[1-2]。國內(nèi)外均受到了新型冠狀病毒帶來的巨大干擾，而該病毒與SARS[3]、中東呼吸綜合征(MERS)[4]、HIV[5]等傳染性疾病有所不同。SARS這類病毒主要表現(xiàn)為發(fā)病2～3天之后出現(xiàn)明顯發(fā)熱，發(fā)熱程度比較高，病人呼吸困難，病情相對比較重，致死率相對較高。而新型冠狀病毒肺炎的潛伏期一般在14天之內(nèi),病毒傳染能力很強，傳播能力比SARS、埃博拉病毒、H7N9強,傳染性較高導(dǎo)致防疫非常困難。通過國家衛(wèi)生健康委員會[6]公布的數(shù)據(jù)，截至3月21日24時，全國現(xiàn)有確診人數(shù)5 549人，現(xiàn)有疑似人數(shù)118人，累計治愈人數(shù)72 244人，尚在醫(yī)學(xué)觀察10 071人。由于國家采取了有效的管控策略，我國的疫情形勢正在逐漸好轉(zhuǎn)。目前國外疫情則令人關(guān)注，疫情形勢相對較重，各個國家正在積極地應(yīng)對，同時我國也派出了相關(guān)專家援助共同對抗疫情。

1 相關(guān)工作

從疫情開始至今，已有相關(guān)研究團隊對新型冠狀病毒展開建模并進行分析，范如國[7]等對新冠病毒展開建模，對中國全境、湖北省尤其武漢市數(shù)據(jù)進行研究，分析了不同潛伏期天數(shù)對感染人數(shù)以及疫情拐點到來的時間影響。耿輝[8]等研究了國家采取的防護措施對病毒傳播帶來的影響，對比了停工停學(xué)不限制出行和限制出行對易感者、潛伏者、感染者和移除者的變化趨勢。吉兆華[9]等用室模型擬合了全國及廣東、河南、江蘇、湖北、浙江5個省份1月份報告發(fā)病情況，根據(jù)模型預(yù)測病例數(shù)對實際病例數(shù)的解釋度和分析群體防控度。張琳[10]則采用一般增長模型分3個階段擬合了確診人數(shù)，發(fā)現(xiàn)確診人數(shù)在經(jīng)歷了初期(1月15日～1月27日)的無障礙指數(shù)增長，中期(1月27日～2月6日)的次指數(shù)增長后，已在2月6日進入了次線性增長階段。周濤[11]等在對武漢新型冠狀病毒感染肺炎基本再生數(shù)進行了預(yù)測，估計COVID-19的基本再生數(shù)在2.8～3.3之間。隨著COVID-19疫情的發(fā)展與防控策略的不斷加強，新型冠狀病毒肺炎傳播機理還需進一步研究。因此，研究針對如何科學(xué)高效地管控疫情，分析了國內(nèi)外早期的疫情數(shù)據(jù)，對已有的數(shù)據(jù)進行數(shù)學(xué)建模，進而預(yù)測未來疫情的發(fā)展趨勢，力圖為疫情防控提供理論支持和實踐指導(dǎo)。

2 基于貝葉斯優(yōu)化SEIR模型的新冠肺炎疫情分析

2.1 數(shù)據(jù)來源

研究所使用的疫情數(shù)據(jù)來自于國家衛(wèi)生健康委員會、湖北省衛(wèi)健委[12]、湖南省衛(wèi)健委[13]、安徽省衛(wèi)健委[14]、陜西省衛(wèi)健委[15]、甘肅省衛(wèi)健委[16]以及根據(jù)github上Blankerl提供的數(shù)據(jù)接口(https://github.com/BlankerL/DXY-COVID-19-Data)所獲得的疫情數(shù)據(jù)。依據(jù)疫情發(fā)展情況的不同，國內(nèi)分析1月20日～2月16日的疫情數(shù)據(jù)，而國外研究2月22日～3月22日的疫情數(shù)據(jù)。

2.2 模型參數(shù)分析

通過模型辨識和參數(shù)估計，對SEIR模型中關(guān)鍵參數(shù)進行了分析，參數(shù)k代表一個感染者平均接觸人數(shù)，k值的變化能夠體現(xiàn)國家采取的防護措施帶給病毒傳播的影響。比如延長春節(jié)假期、企業(yè)延遲復(fù)工、學(xué)校延期開學(xué)、限制出行、居家隔離等，這些措施最大化地扼制了病毒的傳播。而影響我國疫情發(fā)展的重要影響因素就是春運，春運前后全國及各省感染者平均接觸人數(shù)隨時間的變化如圖1所示。由圖1可以看出，1月24日之前完成的春運接觸效果開始在1月25日之后顯現(xiàn)，多地確實出現(xiàn)了春運帶來的影響，k值有所上升。總體來說，全國的k值相對比較平穩(wěn)，k約為10.23人，可以作為平均接觸人數(shù)的參考。但是在1月24日之后，隨著疫情的發(fā)展，國家采取了“封城”舉措，全國多地采取一級響應(yīng)，數(shù)學(xué)建模的參數(shù)需要動態(tài)調(diào)整，采取分層抽樣來計算出新的k值。根據(jù)病毒傳染程度分成5個等級，一級為最嚴(yán)重，五級為輕微，然后進行加權(quán)平均，結(jié)果計算可得k_new約為8.70人。另外，還計算了國家采取的種種防控手段對平均接觸人數(shù)帶來的影響因子(k_changed_ratio)，利用公式k_changed-ratio=k/k_new，計算得出k_changed_ratio約為1.16。k_new的值相比春運之前的k值有所減小，說明國家采取的管控措施對病毒的抑制起到了一定的影響，并在之后效果會更加顯著。疫情嚴(yán)重程度的關(guān)鍵是防控因子的研究，為疾病預(yù)防控制部門提供重要的決策依據(jù)。

2.3 貝葉斯優(yōu)化方法

(1)貝葉斯優(yōu)化過程分析。針對研究需要解決預(yù)測疫情感染人數(shù)以及拐點到來的問題，提出了貝葉斯優(yōu)化的方法，該方法能夠有效快速地尋找到所要解決問題的最優(yōu)參數(shù)，省去了根據(jù)人工經(jīng)驗調(diào)參所花費的時間。貝葉斯優(yōu)化主要面向的問題場景是：

X*=argmaxf(x)(x∈S),

(1)

式中,S是x變量的候選集，即參數(shù)x可能取值的集合。目標(biāo)從集合S中選取一個x，使得f(x)的值最大或最小。這里f(x)具體公式形態(tài)可能無法得知，即黑盒函數(shù)。但是可以選擇一個x，通過實驗或者觀察得出f(x)的值。貝葉斯優(yōu)化的流程圖如圖2所示。

圖2 貝葉斯優(yōu)化流程圖

貝葉斯優(yōu)化[17-18]有兩個核心過程，先驗函數(shù)(Prior Function,PF)與采集函數(shù)(Acquisition Function,AC)，采集函數(shù)也叫效能函數(shù)(Utility Funtcion)。在貝葉斯決策理論的框架下，許多采集函數(shù)可以解釋為在點x評估f相關(guān)的預(yù)期損失，然后通常選擇具有最低預(yù)期損失的點。PF主要利用高斯過程回歸，AC主要包括EI、PI、UCB這幾種方法，接下來簡要介紹一下AC的幾種常用方法。

PI(Probability of Improvement):假設(shè)f′=minf這個f′表示目前已知的f的最小值。接著定義采集函數(shù)如下：

(2)

把u(x)理解成一個獎勵函數(shù)，如果f(x)不大于f′就有獎勵，反之沒有。改進的采集函數(shù)作為變量x的期望函數(shù)如下：

(3)

最后aPI(x)的最大值即可求出基于高斯分布滿足要求的變量x。

EI(Expected Improvement):上述的PI函數(shù)的缺點是有可能找到的是局部最優(yōu)點，而不是全局最優(yōu)點。而EI函數(shù)則可以找出全局最優(yōu)點。f′的定義和上述一樣，f′=minf。但是采集函數(shù)如下：

u(x)=max(0,f′-f(x)),

(4)

最終關(guān)于變量x的采集函數(shù)如式(5)所示：

(f′-u(x))Φ(f′;u(x),K(x,x))+K(x,x))N(f′;u(x),K(x,x)),

(5)

通過計算使得aEI值最大的點即為最優(yōu)點。式(5)中有兩個組成部分,要使得式(5)值最大則需要同時優(yōu)化左右兩個部分：左邊需要盡可能地減少μ(x)；右邊需要盡可能地增大方差(或協(xié)方差)K(x,x)。在探索(Exploration)和利用(Exploitation)此類問題上是一個典型理論。

UCB(Upper Confidence Bound):UCB可以簡單地理解為上置信邊界，UCB通常采用最大化f而不是最小化f來描述。但是在最小化的情況下，采集函數(shù)將采取以下形式：

aUCB(x)=u(x)-βσ(x),

(6)

(2)貝葉斯優(yōu)化結(jié)果分析。為了提高SEIR模型計算效率和分析精確度，需要優(yōu)化α和β這些主要參數(shù)。常用的調(diào)參方式有網(wǎng)格搜索(Grid Search)和隨機搜索(Random Search)。Grid Search是全空間掃描，但較慢。而Random Search雖然快，但可能錯失空間上的一些重要的點，精度不夠。因此，研究選擇貝葉斯優(yōu)化方法來調(diào)整參數(shù)。該算法速度較快、效果較好。貝葉斯優(yōu)化利用平滑性而無需計算梯度，可處理大量變量并行優(yōu)化。研究分別分析我國和意大利的優(yōu)化參數(shù)α和β，將二者迭代優(yōu)化次數(shù)均設(shè)置為1 000次，經(jīng)過1 000次的迭代優(yōu)化之后可以得出我國和意大利α和β的參數(shù)結(jié)果隨著時間步的變化分別如圖3和圖4所示。由圖3可以看出，我國疫情隨著參數(shù)α和β的迭代次數(shù)不斷增加，α在0.26～0.30越來越密集，β在0.050～0.100越來越密集，等到1 000次迭代優(yōu)化之后，可以得出α和β的最優(yōu)解分別為0.28和0.08。由圖4可以看出，意大利疫情隨著參數(shù)α和β的迭代次數(shù)不斷增加，α在0.16～0.20越來越密集，而β在0.050～0.075越來越密集，等到1 000次迭代優(yōu)化之后，可以得出α和β的最優(yōu)解分別為0.18和0.06。

圖3 中國α、β隨著時間步t的變化圖

圖4 意大利α、β隨著時間步t的變化圖

2.4 基于貝葉斯優(yōu)化的SEIR模型

在傳染病SEIR模型中，種群(Population)內(nèi)的N個個體的狀態(tài)可分為如下幾類：(1)易染狀態(tài)S(Susceptible)，即健康狀態(tài)，可被感染的個體；(2)潛伏狀態(tài)E(Exposed)，處于傳染病潛伏期的個體；(3)感染狀態(tài)I(Infected)，處于感染狀態(tài)的個體還能夠感染健康狀態(tài)的個體；(4)移除狀態(tài)R(Removed,Refractory or Recovered)，為被隔離或因病愈而具有免疫力的人。傳統(tǒng)的SEIR模型只是考慮了有潛伏者(Exposed)這種人群，但是并沒有考慮到潛伏者也具有傳染性。為了適應(yīng)COVID-19疫情發(fā)展與不斷加強的防控策略，研究中的SEIR模型考慮潛伏者的傳染性，根據(jù)春運防護前后病毒傳播機制不同的實際問題，模型分為春運前后兩個階段，利用貝葉斯優(yōu)化方法更新參數(shù)，改進的SEIR模型如圖5所示，模型的微分方程如式(7)所示。

圖5 改進的SEIR模型示意圖

(7)

3 結(jié)果與分析

新冠病毒治愈率、死亡率、重癥率變化如圖6所示。由圖6可知，我國新冠肺炎的整體趨勢隨著時間的推移、治療經(jīng)驗和手段的不斷豐富、全國醫(yī)院緊急馳援武漢，治愈人數(shù)比例有了明顯的上升，重癥和死亡人數(shù)比例在不斷下降，感染人數(shù)不再上升，綜合表明病毒疫情正在被逐漸控制。將所有省份按確診人數(shù)進行新型冠狀病毒感染情況等級的劃分，一級為最嚴(yán)重，五級為最輕微，采用分層抽樣的方法抽取了湖北、湖南、安徽、陜西和甘肅5個省份來進行分析。研究將針對各個省份的確診人數(shù)、治愈率、確診人數(shù)增長率等方面進行分析。4個省在1月20日～2月16日確診人數(shù)變化圖如圖7所示。由圖7可知，湖南和安徽的確診人數(shù)相較其他兩省增長比較明顯，原因和每個省份的地理位置密切相關(guān)，前面兩省與湖北省接壤，地理位置比較靠近，所以病毒傳播的比較快速，而后兩省離得相對較遠(yuǎn)，所以傳播速度慢，同時各個省份積極應(yīng)對此次疫情，從圖7中還可以看出，確診人數(shù)慢慢變得平穩(wěn)起來。

圖6 我國新冠病毒治愈率、死亡率、重癥率變化圖圖7 4個省在1月20日～2月16日期間確診人數(shù)變化圖

截至目前，中國疫情形勢正在逐漸好轉(zhuǎn)，作為中國疫情發(fā)源地的武漢也已經(jīng)連續(xù)多天新增人數(shù)為0，但相同時間的國外疫情變得非常嚴(yán)重，意大利、美國、法國、德國、伊朗等國家確診人數(shù)在不斷地攀升。研究選取了疫情發(fā)展相對比較嚴(yán)重的5個國家進行分析，即意大利、美國、法國、德國和伊朗。通過分析新增人數(shù)這個指標(biāo)來反映疫情的嚴(yán)重程度。國外新冠病毒從2月22日～3月22日的新增人數(shù)趨勢變化圖如圖8所示，從圖8中可以清晰地看出，在第5天(2月27日)的時候，5個國家新增人數(shù)幾乎為0，但是從第20天(3月13日)開始，新增人數(shù)趨勢愈加明顯，5個國家中除了伊朗比較穩(wěn)定之外，其他的國家都以不同程度在增長，而新增人數(shù)最明顯的就是美國，單日新增人數(shù)能夠高達上萬人，從而可以推測出美國未來一段時間內(nèi)疫情會繼續(xù)呈爆發(fā)式增長。因此，各國政府必須采取強有力的措施來更好地應(yīng)對和管控疫情。

圖8 國外新冠病毒2月22日～3月22日新增人數(shù)變化圖 圖9 中國和意大利治愈率對比圖

據(jù)目前全球疫情數(shù)據(jù)可知，意大利現(xiàn)在是國外疫情非常嚴(yán)重的國家之一。中國和意大利治愈率對比圖如圖9所示。由圖9可以看出，中國和意大利的治愈率都在逐步地提高，但是在疫情開始的20天之內(nèi)，意大利的治愈率相比我國來說高一點。但是在第20天之后，我國的治愈率遠(yuǎn)超意大利，說明我國投入大量的醫(yī)療設(shè)備以及全國醫(yī)院馳援武漢等強有力的防控措施使疫情得到大大減緩。意大利可以借鑒我國的做法，盡早盡快控制疫情。

研究基于貝葉斯優(yōu)化SEIR模型擬合了國內(nèi)自1月20日～2月16日的數(shù)據(jù)，結(jié)果如圖10所示。通過圖10a可以看出，該模型能夠較好地擬合真實數(shù)據(jù)，預(yù)測峰值在79 107人，預(yù)測拐點在2月27日，可以較好地用來作為第二階段模型預(yù)測的依據(jù)。由圖10b可以看出，曲線擬合貼近真實數(shù)據(jù)，預(yù)測數(shù)據(jù)顯示峰值是在81 314人，疫情拐點在2月25日。春運前與春運后的與第一階段相比預(yù)測拐點稍有變化，這與治愈率γ的變化有關(guān)。隨著γ的提高，疫情的拐點提前出現(xiàn)并且峰值也將穩(wěn)定或者下降。但是總體來說，峰值穩(wěn)定在2月25日至3月上旬之間。意大利確診人數(shù)及拐點預(yù)測結(jié)果如圖11所示。由圖11可見，改進的SEIR模型對意大利確診人數(shù)的擬合結(jié)果，顯示峰值約在14萬人左右，而疫情拐點也將在4月21日左右顯現(xiàn)，這說明目前意大利的疫情形式還是比較嚴(yán)峻，需要加大防控力度，認(rèn)真對待此次疫情，才能盡早地控制住疫情。

圖10 中國1月25日前、后確診人數(shù)及拐點擬合結(jié)果

圖11 意大利確診人數(shù)及拐點預(yù)測結(jié)果

4 總結(jié)

研究基于SEIR模型來對新型冠狀病毒預(yù)測與分析，同時運用了貝葉斯優(yōu)化方法來搜索模型的最優(yōu)超參數(shù)，提高了模型預(yù)測的效果。分析了國內(nèi)外疫情的發(fā)展變化情況，分別對中國和意大利的疫情進行仿真和拐點預(yù)測，國內(nèi)的預(yù)測圖反應(yīng)了管控措施給疫情帶來的影響，使拐點提前并且在其之后感染人數(shù)下降速率加快。結(jié)果表明，模型適用性較優(yōu)，能夠較好地擬合曲線并對短期內(nèi)疫情進行預(yù)測，能夠提供一定的理論依據(jù)和指導(dǎo)意義。但是由于目前疫情還有很多未知因素，模型不可避免地會與現(xiàn)實存在一定差異，導(dǎo)致結(jié)果可能會存在一些誤差。國外形勢依舊嚴(yán)峻，應(yīng)該加大防控力度，防護措施應(yīng)做到快速并且到位，可以借鑒中國的做法，才能更好地應(yīng)對此次疫情。