基于貝葉斯優化SEIR模型的新型冠狀病毒預測與分析

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

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

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

1 相關工作

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

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

2.1 數據來源

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

2.2 模型參數分析

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

2.3 貝葉斯優化方法

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

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

(1)

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

圖2 貝葉斯優化流程圖

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

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

(2)

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

(3)

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

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

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

(4)

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

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

(5)

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

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

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

(6)

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

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

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

2.4 基于貝葉斯優化的SEIR模型

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

圖5 改進的SEIR模型示意圖

(7)

3 結果與分析

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

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

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

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

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

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

圖10 中國1月25日前、后確診人數及拐點擬合結果

圖11 意大利確診人數及拐點預測結果

4 總結

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