基于SSA-PSO-GSA-SVR的全球新冠肺炎疫情分析及預(yù)測(cè)

成云云， 白艷萍， 續(xù) 婷， 譚秀輝， 程 蓉

(中北大學(xué) 理學(xué)院， 山西 太原 030051)

0 引言

時(shí)至今日，全球報(bào)告的新增病例數(shù)仍持續(xù)上升，截止到目前為止，全球累計(jì)確診病例已超過(guò)4.4億，僅在過(guò)去一周(2月22日28日)，全球新增近700多萬(wàn)例新病例，死亡人數(shù)新增近4萬(wàn)人[1]，相當(dāng)于每分鐘就有近700人被感染新冠病毒，嚴(yán)重危害著人們的身體與身心健康.因此，對(duì)疫情傳播及發(fā)展趨勢(shì)進(jìn)行全面、準(zhǔn)確、科學(xué)的分析和預(yù)測(cè)，可以為新冠肺炎疫情監(jiān)測(cè)預(yù)警系統(tǒng)的建設(shè)奠定理論基礎(chǔ)，有助于各國(guó)家更好的采取預(yù)防控制措施.

疫情爆發(fā)初期，國(guó)內(nèi)外就有諸多學(xué)者對(duì)新冠肺炎疫情進(jìn)行了分析與預(yù)測(cè)[2-5]，有關(guān)疫情傳播及預(yù)測(cè)的研究方法主要由傳染病動(dòng)力學(xué)模型和時(shí)間序列機(jī)器學(xué)習(xí)模型為主要框架.在動(dòng)力學(xué)模型中，SEIR模型作為經(jīng)典的傳染病模型被廣泛使用[6-8].但使用傳染病動(dòng)力學(xué)模型需要收集和分析大量相關(guān)的參數(shù)，如R0、就診率、隔離率等.對(duì)于系統(tǒng)中人口數(shù)量較多的情況，由于疫情數(shù)據(jù)的有限性與不完整性，使得動(dòng)力學(xué)模型受到一定的限制.而在機(jī)器學(xué)習(xí)模型中，作為人工智能的核心研究領(lǐng)域，它可以讓計(jì)算機(jī)具備從海量原始數(shù)據(jù)中提取特征因素、主動(dòng)尋求系統(tǒng)內(nèi)部規(guī)律的能力，目前已被廣泛應(yīng)用于傳染病預(yù)測(cè)領(lǐng)域[9-11].因此，將其應(yīng)用到新冠肺炎預(yù)警系統(tǒng)中可以突破傳統(tǒng)傳染病模型的局限性，解決其建模步驟繁雜、考慮因素單一以及提取空間特征受限等問(wèn)題，從而達(dá)到實(shí)時(shí)、動(dòng)態(tài)預(yù)警的目的[12-14].

然而，新冠肺炎傳播受多種因素交互影響，具有復(fù)雜的非線(xiàn)性特征，且數(shù)據(jù)樣本點(diǎn)較少.針對(duì)以上特性，結(jié)合SVR在解決小樣本問(wèn)題上的優(yōu)勢(shì)將其作為基礎(chǔ)預(yù)測(cè)模型，而單一的SVR模型尋優(yōu)過(guò)程耗時(shí)較長(zhǎng)，更新速度較慢，對(duì)模型的預(yù)測(cè)精度產(chǎn)生了一定的影響.因此，本文采用PSO-GSA聯(lián)合算法優(yōu)化SVR的懲罰參數(shù)C和核函數(shù)參數(shù)g，并利用SSA對(duì)原始序列進(jìn)行去噪得到聯(lián)合實(shí)時(shí)動(dòng)態(tài)預(yù)測(cè)模型(SSA-PSO-GSA-SVR).

1 基本原理

1.1 奇異譜分析

奇異譜分析(Singular Spectrum Analysis,SSA)可以提取出代表原時(shí)間序列不同成分的信號(hào)，如長(zhǎng)期趨勢(shì)信號(hào)、周期信號(hào)、噪聲信號(hào)等[15,16]，是一種研究非線(xiàn)性時(shí)間序列的有效方法.假設(shè)新冠肺炎疫情數(shù)據(jù)時(shí)間序列為XT=(x1,x2,…,xT)，其中數(shù)據(jù)序列x1,x2,…,xT為等間隔采樣，長(zhǎng)度T為110，具體計(jì)算步驟如下：

(1)嵌入

將一維的新冠序列數(shù)據(jù)映射到多維軌跡矩陣上，如式(1)所示：

(1)

式(1)中：d=10為嵌入維度，以此得到一個(gè)Xi=[xi,xi+1,…,xi+d-1]的信號(hào)，該步驟的結(jié)果是得到一個(gè)軌跡矩陣Xi，且Xi為Hankel矩陣，表明每一條沿副對(duì)角線(xiàn)的元素都相等.

(2)奇異值分解(SVD)

計(jì)算XXT并求得其d個(gè)特征值λ1,λ2,…,λd，對(duì)應(yīng)的特征向量為U1,U2,…,Ud，且令d′=max(i,λi>0)，由此得到的主成分子空間包含了新冠數(shù)據(jù)序列的不同的特征，其中每個(gè)特征值都代表對(duì)應(yīng)的SVD分量對(duì)原始序列的貢獻(xiàn)率，則

(2)

(3)分組

將式(2)劃分成m個(gè)不同的組以表示不同的新冠肺炎疫情數(shù)據(jù)序列的節(jié)律、噪聲等，并將每組內(nèi)所包含的矩陣相加，則

X=XI1+XI2+…+XIm

(3)

XI=Xi1+Xi2+…+Xip

(4)

(4)重構(gòu)

運(yùn)用對(duì)角平均將矩陣還原為長(zhǎng)度為N的時(shí)間序列，設(shè)Y為d×K的矩陣，其元素為yij.令d*=min{d,K},K*=max{d,K}.將矩陣Y變換為Y={y1,y2,…,yN}，式(5)為對(duì)角平均公式：

(5)

1.2 萬(wàn)有引力搜索算法

萬(wàn)有引力搜索算法(GSA)[17,18]是一種模擬粒子之間萬(wàn)有引力現(xiàn)象的元啟發(fā)式智能優(yōu)化算法，其原理為：把搜索空間中運(yùn)動(dòng)的粒子看成是優(yōu)化問(wèn)題的潛在解，當(dāng)群體之間存在慣性質(zhì)量最大的粒子時(shí)，其他慣性質(zhì)量較小的粒子不斷趨向靠近慣性質(zhì)量最大的那個(gè)粒子，從而使算法收斂至最佳解.如圖1所示，粒子Mi所受到各類(lèi)粒子的合力F1的方向偏向于慣性質(zhì)量最大的粒子M4.

圖1 萬(wàn)有引力作用圖

(6)

式(6)中：λ為[0,1]中的一個(gè)隨機(jī)數(shù)，β為[0,1]中的一個(gè)均勻隨機(jī)分布.根據(jù)牛頓定律，得到粒子慣性質(zhì)量的更新公式：

(7)

(8)

式(8)中：fiti(t)為粒子i在t次迭代時(shí)的適應(yīng)度函數(shù)值，best(t)和worst(t)分別為最佳和最劣的適應(yīng)度函數(shù)值，定義如下：

(9)

以SVR的均方誤差MSE作為適應(yīng)度函數(shù)，首先初始化每個(gè)粒子的位置和速度，并設(shè)置算法的最大迭代次數(shù)和參數(shù)；然后計(jì)算粒子的適應(yīng)度值，根據(jù)文中公式(9)找到best(t)與worst(t)來(lái)更新Mi(t)；其次計(jì)算粒子所受到的合力、加速度，再由公式(6)分別更新粒子的速度和位置；以此迭代循環(huán)，直到新冠肺炎疫情預(yù)測(cè)值的均方誤差得到最小的值，終止迭代.

1.3 支持向量回歸

支持向量回歸(SVR)[19,20]是基于支持向量機(jī)(SVM)的一個(gè)回歸預(yù)測(cè)算法，在解決小樣本、非線(xiàn)性模式識(shí)別方面具有獨(dú)特的優(yōu)點(diǎn).其核心思想為：引入核函數(shù)將樣本空間中非線(xiàn)性回歸映射為高維空間中的線(xiàn)性回歸問(wèn)題,即尋找一個(gè)最優(yōu)分類(lèi)超平面,使得所有訓(xùn)練樣本離該最優(yōu)超平面的誤差最小.給定新冠肺炎疫情數(shù)據(jù)作為訓(xùn)練集D={(x1,y1),…,(xT,YT)}，xi∈Rd，yi∈R，i=1,2,…,T,并使用原始疫情數(shù)據(jù)訓(xùn)練SVR，得到一個(gè)最優(yōu)回歸函數(shù)：

f(x)=ωTx+b

(10)

式(10)中：w為權(quán)值向量，b為偏置向量，根據(jù)輸入的影響因素xi，由yi=f(xi)，得到新冠疫情預(yù)測(cè)數(shù)據(jù)yi，通過(guò)引入松弛因子，SVR的目標(biāo)優(yōu)化函數(shù)可表示為：

(11)

(12)

通過(guò)交叉驗(yàn)證選擇SVR的最佳參數(shù)C和g，在實(shí)驗(yàn)中將C和g在2的指數(shù)范圍網(wǎng)格內(nèi)進(jìn)行搜索，先在大范圍內(nèi)粗略尋找最佳參數(shù)，縮小C和g的取值范圍，然后在此基礎(chǔ)上進(jìn)行精細(xì)的參數(shù)選擇.

2 本文提出的SSA-PSO-GSA-SVR預(yù)測(cè)算法

綜上理論基礎(chǔ)和分析，為提高全球新冠肺炎疫情預(yù)測(cè)模型的預(yù)測(cè)精度，本文提出了一種SSA-PSO-GSA-SVR聯(lián)合預(yù)測(cè)算法，圖2為該預(yù)測(cè)模型整體框架，對(duì)原始序列進(jìn)行SSA去噪得到重構(gòu)序列，然后構(gòu)建PSO-GSA聯(lián)合算法尋找SVR最優(yōu)的懲罰參數(shù)C和核函數(shù)參數(shù)g，達(dá)到優(yōu)化目的，提出的SSA-PSO-GSA-SVR預(yù)測(cè)算法主要分為兩大塊，其算法流程如下：

圖2 SSA-PSO-GSA-SVR預(yù)測(cè)算法流程圖

(1)SSA降噪

運(yùn)用SSA算法將一維新冠疫情時(shí)間序列數(shù)據(jù)轉(zhuǎn)換為多維度序列進(jìn)行主成分分析，將有效的信號(hào)成分進(jìn)行分組處理，并根據(jù)累計(jì)貢獻(xiàn)率選擇若干分量進(jìn)行重構(gòu)，然后在重構(gòu)數(shù)據(jù)集上進(jìn)行預(yù)測(cè)訓(xùn)練.通過(guò)SSA分解重構(gòu)得到的主分量較為規(guī)則，具有更好的可預(yù)測(cè)性.

(2)PSO-GSA優(yōu)化SVR模型

根據(jù)收集到的新冠肺炎疫情數(shù)據(jù)，使用GSA-PSO-SVR預(yù)測(cè)模型進(jìn)行訓(xùn)練時(shí)，選取前94組為訓(xùn)練數(shù)據(jù)，后16組數(shù)據(jù)為測(cè)試數(shù)據(jù).通過(guò)訓(xùn)練前94組數(shù)據(jù)，得到訓(xùn)練后的預(yù)測(cè)模型，用以預(yù)測(cè)全球新冠疫情發(fā)展趨勢(shì)，通過(guò)Matlab軟件，編寫(xiě)GSA-PSO-SVR預(yù)測(cè)算法程序，以實(shí)現(xiàn)以下尋優(yōu)步驟：

PSO算法在求解優(yōu)化問(wèn)題時(shí)易陷入局部最優(yōu)，GSA算法具有局部?jī)?yōu)化能力差、收斂過(guò)早等缺點(diǎn).為解決兩種算法在單獨(dú)求解優(yōu)化問(wèn)題時(shí)的局限性，構(gòu)建了一種更為高效的算法-PSO-GSA聯(lián)合算法，其速度和位置具體更新公式為：

Vi(t+1)=ωVi(t)+c1r1aci(t)+

c2r2(Pbest-Xi(t))

(13)

Xi(t+1)=Xi(t)+Vi(t+1)

(14)

由PSO-GSA算法原理：在迭代初期，粒子的慣性權(quán)重ω比較大，具有較好的全局搜索能力，可以較快地確定最優(yōu)解的大致區(qū)域；隨著迭代次數(shù)的累加，ω的值越來(lái)越小，加強(qiáng)了粒子的局部開(kāi)發(fā)能力，從而粒子群在可行解區(qū)域內(nèi)通過(guò)微調(diào)即可找到全局最優(yōu)解.

3 新冠肺炎傳播仿真模型的數(shù)據(jù)驗(yàn)證

3.1 數(shù)據(jù)樣本

本文采用的實(shí)驗(yàn)數(shù)據(jù)為世界衛(wèi)生組織(World Health Organization,WHO)官網(wǎng)公布的全球新冠肺炎疫情數(shù)據(jù)，按照世衛(wèi)組織工作區(qū)域?qū)⑷騽澐譃榱鶄€(gè)地區(qū)：美洲區(qū)域、歐洲區(qū)域、東地中海區(qū)域、東南亞區(qū)域、非洲區(qū)域、西太平洋區(qū)域，如圖3所示.以國(guó)際周為時(shí)間單位(周日為一周的開(kāi)始)，2020年1月26日的疫情數(shù)據(jù)為分析原點(diǎn)，以此類(lèi)推至2022年2月22日至2月28日為第110周，實(shí)驗(yàn)采用前85%數(shù)據(jù)進(jìn)行模型訓(xùn)練，15%數(shù)據(jù)用于模型的測(cè)試.根據(jù)新冠肺炎每周新增確診以及死亡病例用來(lái)度量疫情的發(fā)展趨勢(shì)，定義增量病例為每周累計(jì)確診及死亡人數(shù)相對(duì)前一時(shí)間單位的變化量.

圖3 全球世衛(wèi)組織六個(gè)區(qū)域總部分布圖

3.2 評(píng)價(jià)指標(biāo)

為對(duì)比評(píng)估各個(gè)模型的預(yù)測(cè)效果，選取平均絕對(duì)百分比誤差(MAPE)作為評(píng)價(jià)指標(biāo)，計(jì)算公式如下：

(15)

3.3 預(yù)測(cè)結(jié)果對(duì)比分析

首先對(duì)全球六個(gè)區(qū)域原始新增確診及死亡數(shù)據(jù)序列進(jìn)行奇異譜分解，由于監(jiān)測(cè)數(shù)據(jù)共110組，為避免過(guò)度分解，設(shè)置窗口長(zhǎng)度d為10.目前重構(gòu)信號(hào)分量的選取方法有很多，本文以貢獻(xiàn)率為依據(jù)選取重構(gòu)分量的個(gè)數(shù)，其中λi代表新增確診病例序列貢獻(xiàn)率，βi代表新增死亡病例序列貢獻(xiàn)率，如表1所示.

表1 六個(gè)地區(qū)原始序列組別貢獻(xiàn)率

從表1可看出，前2個(gè)或3個(gè)重構(gòu)分量的累計(jì)貢獻(xiàn)率均已經(jīng)達(dá)到了97%，剩余組貢獻(xiàn)率減小趨勢(shì)較為緩慢且均小于1，因此選取前2個(gè)或3個(gè)分量重構(gòu)監(jiān)測(cè)序列，剩余組為噪聲序列，根據(jù)圖4和圖5可得重構(gòu)后的序列趨勢(shì)與原始序列趨勢(shì)保持一致.

(a)美洲地區(qū)

(a)美洲地區(qū)

在這一部分，首先以六個(gè)區(qū)域新冠肺炎確診及死亡觀測(cè)樣本集為案例驗(yàn)證了SVR模型在疫情趨勢(shì)預(yù)測(cè)中的有效性,并對(duì)比基于粒子群優(yōu)化支持向量回歸(PSO-SVR)和萬(wàn)有引力搜索算法優(yōu)化支持向量回歸(GSA-SVR)以及基于粒子群和萬(wàn)有引力搜索算法優(yōu)化支持向量回歸(PSO-GSA-SVR)的預(yù)測(cè)性能，最后將PSO-GSA-SVR聯(lián)合預(yù)測(cè)算法應(yīng)用到基于奇異譜分析(SSA)重構(gòu)后的序列，構(gòu)建SSA-PSO-GSA-SVR預(yù)測(cè)模型.

結(jié)合經(jīng)驗(yàn)公式和多次試驗(yàn)對(duì)比，設(shè)置PSO算法和GSA算法的參數(shù)，種群大小和最大迭代次數(shù)分別設(shè)置為20和200，G0設(shè)置為10，為驗(yàn)證PSO-GSA-SVR算法的穩(wěn)定性,對(duì)模型運(yùn)行20次后取平均.選取SVR的MSE為適應(yīng)度函數(shù)，測(cè)試和驗(yàn)證PSO-GSA算法對(duì)SVR參數(shù)的尋優(yōu)性能，當(dāng)達(dá)到最大迭代次數(shù)時(shí)，輸出最小MSE值對(duì)應(yīng)的C和g(最優(yōu))，并將其賦給SVR進(jìn)行預(yù)測(cè)，最終，五種模型對(duì)各區(qū)域確診病例與死亡病例的預(yù)測(cè)結(jié)果如圖6和圖7所示.

(a)美洲區(qū)域確診病例預(yù)測(cè)結(jié)果

(a)美洲區(qū)域死亡病例預(yù)測(cè)結(jié)果

由圖6和圖7可知，基于不同區(qū)域新冠肺炎疫情數(shù)據(jù)的五種算法都能有效地預(yù)測(cè)其發(fā)展趨勢(shì)，且訓(xùn)練后基于SSA-PSO-GSA-SVR模型的擬合曲線(xiàn)與原實(shí)際新冠疫情數(shù)據(jù)趨勢(shì)擬合效果最好，其逼近程度遠(yuǎn)高于其他四種預(yù)測(cè)模型，表明SSA-PSO-GSA-SVR預(yù)測(cè)模型對(duì)疫情趨勢(shì)發(fā)展具有較好的預(yù)測(cè)性能.

根據(jù)表2可得，SVR預(yù)測(cè)算法的評(píng)估指標(biāo)平均絕對(duì)誤差值最大，確診病例和死亡病例的MAPE的值分別為18.979%和15.656%，預(yù)測(cè)效果較差，PSO-SVR與GSA-SVR預(yù)測(cè)模型是分別通過(guò)PSO和GSA優(yōu)化算法尋找最優(yōu)參數(shù)后進(jìn)行SVR預(yù)測(cè)，兩種算法優(yōu)化后的SVR的MAPE值與傳統(tǒng)SVR的MAPE相比，確診病例的MAPE值分別下降了20.986%和25.159%，死亡病例的MAPE值分別下降了19.711%和19.845%，這表明PSO和GSA算法對(duì)SVR參數(shù)優(yōu)化效果明顯；而與PSO-SVR與GSA-SVR預(yù)測(cè)模型相比， PSO-GSA-SVR模型的預(yù)測(cè)精度也有了一定的提高.從表2可以看出，SSA-PSO-GSA-SVR預(yù)測(cè)算法平均絕對(duì)誤差值最小，確診病例和死亡病例的MAPE的值分別為9.719%和6.139%，因此，與文獻(xiàn)[19]和文獻(xiàn)[20]提出的PSO-SVR和GSA-SVR模型的預(yù)測(cè)效果相比，本文提出的SSA-PSO-GSA-SVR模型預(yù)測(cè)精度有所提高，能夠更可靠、更準(zhǔn)確地預(yù)測(cè)新冠肺炎發(fā)展趨勢(shì).

表2 不同預(yù)測(cè)模型的評(píng)價(jià)指標(biāo)

續(xù)表2

4 結(jié)論

針對(duì)全球新冠肺炎疫情傳播問(wèn)題，本文結(jié)合群智能算法和機(jī)器學(xué)習(xí)提出了基于SSA-PSO-GSA-SVR的預(yù)測(cè)模型，用于預(yù)測(cè)各區(qū)域每周新增確診與死亡病例人數(shù)，并采用平均絕對(duì)百分比誤差評(píng)估各模型的預(yù)測(cè)效果.在仿真實(shí)驗(yàn)中，通過(guò)對(duì)比不同區(qū)域內(nèi)SVR、PSO-SVR、GSA-SVR、PSO-GSA-SVR以及SSA-PSO-GSA-SVR預(yù)測(cè)模型的MAPE值，可以看出本文提出的方法預(yù)測(cè)效果較優(yōu)，為新冠疫情的傳播及發(fā)展趨勢(shì)提供了一種有效的預(yù)測(cè)模型.同時(shí)，驗(yàn)證了支持向量回歸在解決小樣本問(wèn)題中的優(yōu)異性，采用基于奇異譜分析SSA和PSO-GSA聯(lián)合算法優(yōu)化SVR的超參數(shù)，降低了序列波動(dòng)對(duì)預(yù)測(cè)造成的影響，提高了預(yù)測(cè)模型的精度和穩(wěn)定性.

新冠肺炎傳播及發(fā)展趨勢(shì)預(yù)測(cè)是疫情監(jiān)測(cè)系統(tǒng)建設(shè)與科學(xué)采取預(yù)防控制措施的重要環(huán)節(jié).在此背景下,多種機(jī)器學(xué)習(xí)算法融合運(yùn)用以及新冠疫情信息的不確定性、間隔性、波動(dòng)性分析等都值得我們深入的思考和研究.為使本文提出的SSA-PSO-GSA-SVR成為一個(gè)更加通用的傳染病疫情傳播及趨勢(shì)預(yù)測(cè)模型，下一步的工作本課題組將會(huì)在現(xiàn)有模型的基礎(chǔ)上從以下幾方面進(jìn)行研究和改進(jìn)：

(1) 采用統(tǒng)計(jì)學(xué)方法深入分析全球各區(qū)域近期疫情嚴(yán)重程度并對(duì)整體疫情走勢(shì)進(jìn)行分類(lèi)；

(2)加入特征因素，目前的大多數(shù)新冠肺炎預(yù)測(cè)研究局限于時(shí)間序列建模，缺少分析外在因素對(duì)新冠傳播的影響，如人口流動(dòng)性、防護(hù)意識(shí)等；

(3) 擴(kuò)展數(shù)據(jù)集至各個(gè)國(guó)家，由于數(shù)據(jù)集的缺失，本文只考慮了確診病例與死亡病例，在未來(lái)工作中將嘗試提取新冠肺炎治愈病例并對(duì)其分析研究.