汪劍眉,李 鋼
(北京郵電大學(xué)經(jīng)濟(jì)管理學(xué)院 北京 海淀區(qū) 100083)
新型冠狀病毒肺炎(COVID-19,以下簡稱新冠肺炎)全球確診病例已超過百萬例,受疫情影響的國家和地區(qū)已達(dá)200 多個。疫情在中國得到有效控制,目前境內(nèi)確診病例僅有零星散發(fā)。在積極干預(yù)的背后,流行病傳播動力學(xué)研究為趨勢判斷和防控決策提供了重要參考。早期研究受限于時間跨度短、干預(yù)成效未完全顯現(xiàn)以及對病毒特性的認(rèn)知局限。目前國內(nèi)疫情態(tài)勢漸趨明朗,具備了更完整的數(shù)據(jù)觀測周期和更豐富的臨床認(rèn)知佐證,使研究具有更多優(yōu)化和拓展空間。
本文以SIR 經(jīng)典倉室模型為邏輯架構(gòu),設(shè)計了分群、分階段的非均勻感染力擴(kuò)展模型,對全國、湖北、武漢3 層次分形網(wǎng)絡(luò)進(jìn)行建模訓(xùn)練,通過降階得到更有效的子過程迭代,在參數(shù)設(shè)置和精度優(yōu)化方面有所突破。在各干預(yù)杠桿強度上下浮動20%的場景設(shè)定下模擬疫情發(fā)展進(jìn)程,評估中國防控模式的成效,并為我國和全球其他國家的后續(xù)防控工作提供有益參考。
在建模和參數(shù)研究方面,文獻(xiàn)[1]應(yīng)用SEIR倉室模型以多重擬合確定參數(shù),結(jié)合LSTM 模型對疫情趨勢進(jìn)行綜合預(yù)測,在參數(shù)估計方面更貼近臨床實踐;文獻(xiàn)[2]對倉室模型進(jìn)行了擴(kuò)展,增加無癥狀感染者A、被隔離易感者Sq等,采用M-H算法進(jìn)行參數(shù)估計建模;文獻(xiàn)[3]在這一研究的基礎(chǔ)上通過歐拉數(shù)值方法進(jìn)行參數(shù)優(yōu)化;文獻(xiàn)[4]根據(jù)吉布斯采樣確定基本再生數(shù),以國際確診病例反推武漢感染人數(shù),并針對感染力進(jìn)行差異化分場景仿真;文獻(xiàn)[5]通過蒙特卡羅方法模擬聯(lián)合分布進(jìn)行參數(shù)推斷和預(yù)測,將研究限定在2020 年1 月23 日前武漢的無干預(yù)自由傳播期,這種階段劃分有效區(qū)隔了病毒傳播早期與中后期的顯著不同。
在控制措施效果評估方面,文獻(xiàn)[6]主要圍繞人口流動性和確診感染病例數(shù)之間相關(guān)性的變化,評估交通管制、切斷本地傳播鏈等舉措成效;文獻(xiàn)[7]將有效再生數(shù)作為評價防控效果的核心指標(biāo),基于蒙特卡洛方法以患者癥狀出現(xiàn)時間和世代間隔對有效再生數(shù)進(jìn)行估計;文獻(xiàn)[8]研究了個人防護(hù)、早期診斷收治和后期診斷收治3 種控制策略組合下病毒的擴(kuò)散速度,通過構(gòu)造目標(biāo)函數(shù)并轉(zhuǎn)化為哈密爾頓函數(shù)極小化最優(yōu)控制問題,發(fā)現(xiàn)全面組合策略更有效。
上述研究為模型的優(yōu)化和突破提供啟發(fā):1)在檢疫隔離和集中收治等政策下各群體感染力不均勻,應(yīng)考慮分群。2) 干預(yù)政策出臺時間點附近疫情曲線呈現(xiàn)明顯彎折,參數(shù)具有顯著分段特征,各時段感染力不均勻,應(yīng)考慮分階段。3) 建模主體選擇方面,根據(jù)國家衛(wèi)健委疫情通報數(shù)據(jù)計算,2020 年3 月25 日之前湖北確診病例中70%以上來自武漢,全國確診病例中70%以上來自湖北,局部與整體具有分形網(wǎng)絡(luò)自相似量化統(tǒng)計特征,如結(jié)合分形架構(gòu),可在參數(shù)估計時充分利用全國、省、市3 層次數(shù)據(jù)。4) 根據(jù)國家和湖北衛(wèi)健委疫情通報數(shù)據(jù)計算,未感染者始終在國內(nèi)人口中占據(jù)主體,武漢約99.56%、湖北約99.89%、全國約99.94%,如能將感染者和未感染者的迭代過程分離,可有效提高模型精度和效率。5) 可將接觸、隔離、收治和治愈視為干預(yù)杠桿,模擬并觀測其在加強和減弱兩個方向上調(diào)節(jié)時疫情態(tài)勢變化,從而描繪可干預(yù)空間、評估干預(yù)成效并實現(xiàn)橫向可比。
基于經(jīng)典SIR 倉室模型[9]進(jìn)行擴(kuò)展,如圖1 所示。為體現(xiàn)感染力異質(zhì)性,根據(jù)人群流動自由度分群:
圖1 SIR 非均勻感染力擴(kuò)展模型的邏輯架構(gòu)
S 為易感者(susceptible),包括自由易感者Sf和隔離易感者Sq,對新冠肺炎人群普遍易感。
I 為感染者(infectious),包括自由感染者If、隔離感染者Iq和收治感染者Ih,按其是否被隔離分別包含于If和Iq中,其中If主要包含潛伏期感染者及始終無癥狀感染者[10-11]。
R 為移除者(removed),包括治愈者和病亡者。
感染路徑僅有一條,即自由感染者If傳染自由易感者Sf。人群轉(zhuǎn)化路徑包括以隔離效率系數(shù)q 監(jiān)測并隔離有密切接觸史或疑似癥狀的自由易感者Sf,如其未感染則成為隔離易感者Sq,并在隔離結(jié)束后再次成為自由易感者Sf(隔離解除率λ 為隔離期的倒數(shù));如其已感染則成為隔離感染者Iq。自由易感者Sf以接觸效率系數(shù)c 和日感染率β 被感染,轉(zhuǎn)化為自由感染者If。If和Iq繼而以收治效率系數(shù)δ 轉(zhuǎn)化為收治感染者Ih。三類感染者分別以治愈效率系數(shù)γ 及日病亡率α 成為移除者R。考慮分形架構(gòu),令i={1,2,3}分別表示武漢、湖北、全國3 個層次。模型微分方程如下:
感染力為每位感染者的日均感染人數(shù),即FOI=βcSf/N=R/DI(在干預(yù)初期為Ro/DI)。
病毒特性參數(shù)包括基礎(chǔ)再生數(shù)Ro、世代間隔DI、日感染率β、日病亡率α,為不變參數(shù),由病毒本身特性決定[12-13],不隨階段改變而改變,其中前3 項參數(shù)在武漢、湖北、全國3 個層次保持一致,日病亡率α 考慮供給保障能力有所不同而保留差異。干預(yù)杠桿參數(shù)包括接觸效率系數(shù)c、隔離效率系數(shù)q、收治效率系數(shù)δ、治愈效率系數(shù)γ,為可變參數(shù),在不同階段、不同層次具有不同的參數(shù)值,從而模擬空間分布上感染力的非均勻性。
一個完整的疫情周期需要經(jīng)歷自由感染者高峰、隔離和收治感染者高峰、移除者新增高峰,根據(jù)峰值點數(shù)據(jù)觀測,將觀察期劃分干預(yù)初期(2020年1 月20 日?2 月4 日)、探索期(2020 年2 月5 日?2 月17 日)、 掌 控 期(2020 年2 月18 日?2 月27 日)、平 穩(wěn) 期(2020 年2 月28 日 之 后) 4 個 階段,從而模擬時間分布上感染力的非均勻性。
時間步t 表示從起始日(1 月20 日)開始計算的第t 天,每日樣本矩陣X=(Sf, Sq, If, Iq, Ih, R)。第t+1 步狀態(tài)變量僅由第t 步狀態(tài)變量決定,視為有限狀態(tài)空間內(nèi)的馬爾可夫過程,x(t+1)=x(t)A,x(t)為n 維狀態(tài)向量,A 為狀態(tài)轉(zhuǎn)移矩陣。因自由易感者人口占比超過99.5%,Sf/N≈1,則有:
分塊降階,則:
由此分離為兩個子過程,即易感者迭代(Sf, Sq)A1、感染與移除者迭代(If, Iq, Ih, R) A4,對于疫情的模擬主要基于后者。兩個子過程數(shù)據(jù)量級相差百倍以上(易感者與感染和移除者的數(shù)量比值武漢超過200 倍,湖北超過800 倍,全國超過17 000倍),分離迭代有利于提高模型精度和計算效率。但子過程轉(zhuǎn)移矩陣首行元素和不為1,不再是完整的馬爾可夫過程。
數(shù)據(jù)來源為國家衛(wèi)健委、湖北衛(wèi)健委、武漢衛(wèi)健委發(fā)布數(shù)據(jù),并進(jìn)行了兩次修正。一是2020 年2 月12 日湖北確診方式增加“臨床診斷”,導(dǎo)致數(shù)據(jù)躍升,按確診和臨床病例比例對前序數(shù)據(jù)進(jìn)行修正;二是2020 年1 月24 日9 家持有資質(zhì)的醫(yī)院和武漢市疾控中心獲得檢測授權(quán),短期內(nèi)大量積壓待檢病例得到確診導(dǎo)致數(shù)據(jù)躍升,按平均增速對前序數(shù)據(jù)進(jìn)行修正。
在獲得病毒特性參數(shù)估計值的基礎(chǔ)上,干預(yù)杠桿參數(shù)估計轉(zhuǎn)化為求各階段線性定常系統(tǒng)狀態(tài)轉(zhuǎn)移矩陣的方程解析解。
通過分形理論的支持,如拓?fù)鋵W(xué)和動力學(xué)領(lǐng)域時變曲線的無標(biāo)度屬性,可以綜合利用武漢、湖北、中國3 個層次的統(tǒng)計數(shù)據(jù)進(jìn)行參數(shù)估計。結(jié)合擬合特征,應(yīng)用監(jiān)督式機器學(xué)習(xí)GPR 高斯過程回歸方法,將各時間步解向量作為樣本對回歸模型進(jìn)行訓(xùn)練,得到4 個階段的參數(shù)估計如表1 所示。
接觸效率系數(shù)在干預(yù)初期受春節(jié)影響明顯上浮,在交通管制、社區(qū)封閉等最大限度控制人口流動和聚集的應(yīng)急干預(yù)下出現(xiàn)幾何級數(shù)下降。隔離效率系數(shù)總體平穩(wěn)略增,武漢以外區(qū)域優(yōu)于武漢,湖北以外區(qū)域優(yōu)于湖北。擬合周期主體處于醫(yī)療資源非擠兌期,全國實行統(tǒng)一診療方案、防控方案,收治效率系數(shù)在探索期出現(xiàn)短期提升,源于增加“臨床診斷”和全面排查應(yīng)收盡收。治愈效率系數(shù)隨診療方案不斷成熟而顯著改善,通過提高醫(yī)療資源周轉(zhuǎn)率,避免收治和隔離環(huán)節(jié)擁塞。
表1 干預(yù)杠桿參數(shù)估計均值
基于參數(shù)估計計算模型在各時期3 個區(qū)域?qū)哟蔚母腥玖等绫? 所示。
表2 各階段非均勻感染力FOI
模型對疫情發(fā)展實際趨勢擬合良好,以此作為基礎(chǔ)場景,如圖2 所示。
圖2 基礎(chǔ)場景模型預(yù)測與實際值比較
分析疫情態(tài)勢對干預(yù)杠桿參數(shù)變化的敏感度SAF[15]。態(tài)勢變量包括感染者峰值規(guī)模NP、峰值周期TP、終值規(guī)模NF、終值周期TF。設(shè)杠桿參數(shù)為x,有:
以現(xiàn)有干預(yù)強度為基礎(chǔ)場景,增強或減弱單項杠桿的強度,使其變化幅度達(dá)到20%,發(fā)現(xiàn)敏感度在趨強和趨弱兩個方向上存在不對稱的邊際效應(yīng),弱干預(yù)場景下邊際損失遞增,強干預(yù)場景下邊際收益遞減,詳見圖3, If1/ If2/If3分別為武漢/湖北/全國自由感染者, Iq1/Iq2/Iq3為分別為武漢/湖北/全國隔離感染者, Ih1/Ih2/Ih3分別為武漢/湖北/全國收治感染者, R1/R2/R3為分別為武漢/湖北/全國移除者。
圖3 干預(yù)杠桿參數(shù)的強干預(yù)場景與弱干預(yù)場景模擬分布
弱干預(yù)場景中,如人口流動性干預(yù)降低20%,全國可達(dá)21 萬累計確診規(guī)模;如對密切接觸者和疑似者的隔離效率降低20%,全國可達(dá)15 萬累計確診規(guī)模;如集中收治效率降低20%,全國可達(dá)10 萬累計確診規(guī)模。強干預(yù)場景下,除降低接觸和提高收治會使峰終值降至基礎(chǔ)場景的37%和80%以外,加強其他兩個杠桿對峰終值影響不大。表3 為干預(yù)杠桿參數(shù)的多場景調(diào)控幅度模擬值。
表3 干預(yù)杠桿參數(shù)的多場景調(diào)控幅度模擬值
表4 干預(yù)杠桿參數(shù)的敏感度
4 類干預(yù)杠桿中,降低接觸是有效抑制疫情的核心杠桿,敏感度是其他干預(yù)杠桿的3.5 倍以上。隔離、收治、治愈屬于防御型杠桿,可維持疫情態(tài)勢平穩(wěn)發(fā)展。
表4 為干預(yù)杠桿參數(shù)的敏感度值。弱干預(yù)敏感度提示了防控的重要風(fēng)險點,降低接觸是防止疫情失控的最基礎(chǔ)手段;其次是隔離和收治,可避免峰終值過載;治愈可避免峰值周期不斷拖延。強干預(yù)敏感度表明了主動防控的提升空間,降低接觸是控制疫情規(guī)模和周期收效最為顯著的手段,及時收治可以在一定程度上抑制終值規(guī)模。
根據(jù)模型的子過程解構(gòu),最大限度控制人口流動和聚集,對有效遏制疫情發(fā)揮了關(guān)鍵作用;通過增加“臨床診斷”、全面排查、加大床位供給實現(xiàn)應(yīng)收盡收,在短期內(nèi)有力扭轉(zhuǎn)了增長勢頭。
根據(jù)敏感度分析,強弱場景邊際效應(yīng)不對稱,全程減弱干預(yù)會引發(fā)較大疫情風(fēng)險,但加強干預(yù)只會帶來較小疫情改善,證明中國方案比較接近準(zhǔn)帕累托最優(yōu)。雖然在降低接觸和集中收治方面理論上存在優(yōu)化空間,但實際受限于干預(yù)初期對病毒特性的發(fā)現(xiàn)與認(rèn)知、探索期對醫(yī)療需求的剛性兌付、掌控期對社會管理的托底投入以及平穩(wěn)期防疫與經(jīng)濟(jì)發(fā)展的再平衡。整體而言,中國模式干預(yù)杠桿組合有力有效。