劉振球 嚴(yán) 瓊 左佳鷺 方綺雯 張鐵軍△

醫(yī)學(xué)研究中的時(shí)間序列是指將某種疾病或者某個(gè)現(xiàn)象的某一指標(biāo)，在每個(gè)時(shí)間點(diǎn)上的數(shù)據(jù)，按照時(shí)間先后順序排列而成的一維數(shù)據(jù)。由于疾病的發(fā)生發(fā)展多與時(shí)間相關(guān)，尤其是傳染病，因此時(shí)間序列往往蘊(yùn)含著一定的規(guī)律。在以往的流行病學(xué)研究中，研究者們習(xí)慣于使用滑動(dòng)平均模型(ARIMA)或者季節(jié)性差分自回歸滑動(dòng)平均模型(SARIMA)描述傳染病的時(shí)間分布特征，并對(duì)未來的發(fā)病情況進(jìn)行預(yù)測(cè)[1]。然而，這些經(jīng)典的回歸分析均暗含一個(gè)假設(shè)，即數(shù)據(jù)是平穩(wěn)的。如果數(shù)據(jù)非平穩(wěn)，極有可能導(dǎo)致“虛假回歸”，基于此做出來的預(yù)測(cè)模型往往與真實(shí)數(shù)據(jù)相差較大。

傳染病的發(fā)病率受到多種因素的影響，比如人口數(shù)、社會(huì)經(jīng)濟(jì)狀況、季節(jié)因素、其他傳染病的發(fā)病情況等[2]，因此由傳染病的發(fā)病率數(shù)據(jù)構(gòu)成的時(shí)間序列，往往也是不平穩(wěn)的。所以，需要尋找一種更好的方法來避免或者減弱各種因素對(duì)序列平穩(wěn)性的影響?；诖?，研究者通過大量的實(shí)驗(yàn)，證實(shí)經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition,EMD)具有類似小波變換中的二進(jìn)濾波器特性，通過分解、重組，實(shí)現(xiàn)對(duì)數(shù)據(jù)的去噪[3]。與此同時(shí)，隨著計(jì)算機(jī)技術(shù)的不斷發(fā)展，機(jī)器學(xué)習(xí)方法被廣泛運(yùn)用在醫(yī)學(xué)領(lǐng)域。以往研究證實(shí)，神經(jīng)網(wǎng)絡(luò)、隨機(jī)森林、支持向量機(jī)等一系列算法模型對(duì)于時(shí)間序列的預(yù)測(cè)，其效果要優(yōu)于傳統(tǒng)的統(tǒng)計(jì)模型[4-5]。本文擬將EMD與BP神經(jīng)網(wǎng)絡(luò)模型結(jié)合，實(shí)現(xiàn)對(duì)傳染病時(shí)間序列數(shù)據(jù)的預(yù)測(cè)并評(píng)估其效果。

原理和方法

1.經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)

EMD是由美國科學(xué)院院士黃鍔博士提出的一種信號(hào)分解方法[6]，廣泛運(yùn)用于自適應(yīng)性時(shí)間序列分析模型。它尤其適合于分析非線性、非平穩(wěn)的信號(hào)序列，具有很高的信噪比。該方法的核心是經(jīng)驗(yàn)?zāi)Ｊ椒纸猓磳⒃緩?fù)雜的信號(hào)數(shù)據(jù)分解為有限個(gè)本征模函數(shù)(intrinsic mode function,IMF)，所分解出來的各個(gè)IMF分量包含了原信號(hào)的不同時(shí)間尺度的局部特征信號(hào)。基于迭代的思想，EMD在分解過程中先找出信號(hào)的全部極值點(diǎn)，利用三次樣條插值求出上下包絡(luò)曲線，并定義信號(hào)的局部包絡(luò)均值為慢振蕩分量。通過不斷減去慢振蕩分量來“篩選”出快振蕩分量。自適應(yīng)地得到一組IMF信號(hào)。一維信號(hào)EMD的具體流程如下:

(1)上包絡(luò)線：找到原始數(shù)據(jù)X(t)的所有極大值點(diǎn)，采用三次樣條插值函數(shù)擬合形成原始數(shù)據(jù)的上包絡(luò)線e(t)max。

(2)下包絡(luò)線：找到原始數(shù)據(jù)X(t)的所有極小值點(diǎn)，同樣采用三次樣條插值函數(shù)擬合形成原始數(shù)據(jù)的下包絡(luò)線e(t)min。

(3)求均值：假設(shè)上下包絡(luò)線的平均值記為m(t)，則m(t)=(e(t)max+e(t)min)/2

(4)構(gòu)建新數(shù)據(jù)：將原始數(shù)據(jù)X(t)減去m(t)，得到新數(shù)據(jù)d(t)，即：d(t)=X(t)-m(t)

(5)判斷d(t)是否滿足IMF的條件，若滿足，則令ci(t)=d(t)；否則以d(t)作為原始數(shù)據(jù)重新回到第一步，繼續(xù)進(jìn)行篩選，重復(fù)k次后，直到d(t)k滿足IMF的條件。

(6)篩選出來的第一個(gè)IMF表示為c1，可以先從原始數(shù)據(jù)中分離出來，即：r1=X(t)-c1

然后將r1作為分解的數(shù)據(jù)，重復(fù)上述過程，得到c2，即：r2=r1-c2，…，rn=rn-1-cn

直到滿足以下兩個(gè)條件之一，分解結(jié)束。

①rn或者cn小于給定的值；

②rn為單調(diào)函數(shù)，無法再從中篩選出更多的IMF。

一般情況下,以連續(xù)兩次迭代得到di(t)和di-1(t)的歸一化均方差來判斷di(t)是否為IMF分量。歸一化均方差定義為:

其中，SD取值一般為0.2～0.3。

2.基于EMD的BP神經(jīng)網(wǎng)絡(luò)

人工神經(jīng)網(wǎng)絡(luò)(artificial neural networks,ANN)作為一種機(jī)器學(xué)習(xí)方法，具有很強(qiáng)的學(xué)習(xí)能力和適應(yīng)能力，在非線性系統(tǒng)的預(yù)測(cè)研究中得到廣泛的應(yīng)用。BP神經(jīng)網(wǎng)絡(luò)是一種按誤差逆?zhèn)鞑ニ惴ㄓ?xùn)練的多層前饋型網(wǎng)絡(luò)，是目前應(yīng)用最為廣泛的神經(jīng)網(wǎng)絡(luò)模型之一[7]。它能學(xué)習(xí)和貯存大量的“輸入—輸出”模式的映射關(guān)系，而且無需事先揭示描述這種關(guān)系的數(shù)學(xué)方程。它的學(xué)習(xí)規(guī)則是使用最速下降法，通過反向傳播來不斷調(diào)整網(wǎng)絡(luò)的權(quán)值和閾值，使網(wǎng)絡(luò)的誤差平方和最小。BP神經(jīng)網(wǎng)絡(luò)模型的拓?fù)浣Y(jié)構(gòu)包括輸入層(input layer)，隱層(hide layer)以及輸出層(output layer)[8]。其中任一神經(jīng)元的輸出Hj可以表示為：

公式中，ωij表示當(dāng)前層和上一層之間的聯(lián)合權(quán)重，i為當(dāng)前層單元，j為上一層單元，n為輸入的層數(shù)，bj是閾值，f表示非線性激活函數(shù)，通常使用sigmoid方程(f(t)=1/(1+e-x))作為激活函數(shù)[9]。

對(duì)于非線性系統(tǒng)，BP神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)有著非常明顯的優(yōu)勢(shì)，但是在復(fù)雜的非線性系統(tǒng)中，非平穩(wěn)因素給預(yù)測(cè)帶來了一定的困難。由于傳染病的發(fā)病率數(shù)據(jù)通常都是非線性非平穩(wěn)的時(shí)間序列數(shù)據(jù)，因此利用神經(jīng)網(wǎng)絡(luò)對(duì)其進(jìn)行預(yù)測(cè)，效果不甚理想。為了提高預(yù)測(cè)的精度，我們使用EMD對(duì)原始數(shù)據(jù)進(jìn)行分解，以降低非平穩(wěn)性對(duì)預(yù)測(cè)精度的影響。對(duì)各個(gè)IMF進(jìn)行單獨(dú)預(yù)測(cè)之后，再將預(yù)測(cè)結(jié)果進(jìn)行合成，得到最終的預(yù)測(cè)結(jié)果。

3.統(tǒng)計(jì)方法

R語言作為一門開源的統(tǒng)計(jì)性計(jì)算機(jī)語言，有專門為經(jīng)驗(yàn)?zāi)B(tài)分解而設(shè)計(jì)的擴(kuò)展包EMD，我們可以利用該擴(kuò)展包快速進(jìn)行時(shí)間序列數(shù)據(jù)的經(jīng)驗(yàn)?zāi)B(tài)分解；除此之外，R語言中也提供了大量的用于BPANN的擴(kuò)展包，比如nnet，AMORE以及caret等。下面將以EMD包結(jié)合nnet和caret包對(duì)數(shù)據(jù)進(jìn)行分析。

實(shí)例分析

1.數(shù)據(jù)來源

本文選取全國傳染病網(wǎng)絡(luò)直報(bào)系統(tǒng)中2004年1月份到2016年12月份全國艾滋病的發(fā)病數(shù)據(jù)，具體資料來源國家衛(wèi)計(jì)委網(wǎng)站http://www.moh.gov.cn/zwgk/yqbb3/ejlist.shtml.

2.EMD包簡介

EMD包中提供了extractimf函數(shù),它從一個(gè)給定的數(shù)據(jù)中分解所有的IMF。在extractimf函數(shù)中,參數(shù)含義如下：

(1)residue：規(guī)定時(shí)間段內(nèi)觀察到的數(shù)據(jù)或信號(hào)。

(2)tt：觀察的時(shí)間。

(3)tol：停止篩選過程的容忍值。

(4)max.sift：最大分解次數(shù)。

(5)stoprule：篩選過程的停止條件,stoprule參數(shù)可以選擇“type1”或“type2”?！皌ype1”指定當(dāng)所有t時(shí)間滿足|di(t)|

(6)boundary：邊界調(diào)整。為了消除信號(hào)的邊界效應(yīng),必須進(jìn)行邊界調(diào)整。Huang通過在信號(hào)兩邊重復(fù)添加人為波段來擴(kuò)展原有信號(hào)，添加部分通過靠近邊界的極值形成的固有模式構(gòu)建[10]。Zeng和He提出兩種擴(kuò)展信號(hào)的方式,它們分別通過對(duì)稱和自反方式添加信號(hào)，被稱為偶擴(kuò)展和奇擴(kuò)展[11]?？梢酝ㄟ^指定boundary參數(shù)為“wave”、“symmetric”、“periodic”、“evenodd”或“none”來改變邊界調(diào)整的方式。

(7)sm、spar、check：分別表示包絡(luò)線的構(gòu)建方式、構(gòu)建包絡(luò)線參數(shù)及是否顯示每次循環(huán)。

最終，通過整合extrema和extractimf函數(shù)就可以實(shí)現(xiàn) EMD，EMD包也提供了emd函數(shù)直接實(shí)現(xiàn)整個(gè)過程。

3.EMD結(jié)果

將艾滋病發(fā)病例數(shù)的數(shù)據(jù)作為信號(hào)進(jìn)行EMD后，得到的結(jié)果如圖1所示。由圖可知，共分解得到4個(gè)IMF，除了四條展示于圖中的曲線外，還有四個(gè)相應(yīng)的矩陣，可供下一步進(jìn)行BPANN分析。圖1中最下方的紅色曲線，即最終EMD分解后剩余的r，以residue表示。這條曲線是將原始數(shù)據(jù)中的IMF剔除后的剩余部分，EMD將原本參差不齊的數(shù)據(jù)分解成了一條代表時(shí)間趨勢(shì)的平滑曲線。

圖1 2004-2016年艾滋病發(fā)病例數(shù)EMD分解結(jié)果

圖2中實(shí)線為2004年1月至2016年12月艾滋病發(fā)病人數(shù)的原始數(shù)據(jù)，虛線為EMD分解后的residue曲線。由圖可知，在2011年之前，我國艾滋病的發(fā)病總體處于一個(gè)緩慢增長的時(shí)期，在2011年到 2014年，經(jīng)歷了一個(gè)快速增長期，2015年之后，增長趨于平緩。

圖2 2004-2016我國艾滋病發(fā)病例數(shù)時(shí)間序列圖

4.BPANN預(yù)測(cè)結(jié)果

以2004-2015的數(shù)據(jù)作為訓(xùn)練集，2016年的數(shù)據(jù)作為測(cè)試集。在R語言caret包中對(duì)EMD分解出的四個(gè)IMF以及殘差分別進(jìn)行神經(jīng)網(wǎng)絡(luò)模型訓(xùn)練和預(yù)測(cè)，并將IMF的預(yù)測(cè)結(jié)果和殘差的預(yù)測(cè)結(jié)果進(jìn)行合成，得到原始數(shù)據(jù)的預(yù)測(cè)結(jié)果。預(yù)測(cè)結(jié)果見圖3，圖中黑色曲線為原始數(shù)據(jù)，依次是對(duì)2016年的原始數(shù)據(jù)進(jìn)行EMD后得到的4個(gè)IMF和殘差，紅色曲線BPANN模型對(duì)以上五組數(shù)據(jù)進(jìn)行預(yù)測(cè)所得到的預(yù)測(cè)值。此外，SARIMA(2,1,1)(2,1,1)12模型同樣被用來對(duì)2016年的原始數(shù)據(jù)進(jìn)行預(yù)測(cè)。進(jìn)行歸一化后的五個(gè)神經(jīng)網(wǎng)絡(luò)模型的基本參數(shù)見表1。此外，我們采用均方根誤差(RMSE)作為衡量模型擬合優(yōu)度的指標(biāo)。RMSE的值越小，則表明模型擬合優(yōu)度越高。

表1 BP神經(jīng)網(wǎng)絡(luò)模型主要參數(shù)和返回結(jié)果

*：decay ：學(xué)習(xí)速率，在0～0.1之間取值。hidden size：隱藏的神經(jīng)元個(gè)數(shù)，數(shù)值越大，模型越復(fù)雜。RMSE：均方根誤差。

圖3 BP神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)模型結(jié)果以及SARIMA預(yù)測(cè)結(jié)果與原始數(shù)據(jù)比較散點(diǎn)圖

討 論

經(jīng)驗(yàn)?zāi)B(tài)分解于1998年提出，目前已被廣泛地運(yùn)用在各種工業(yè)信號(hào)[12]和股票數(shù)據(jù)[13]的分析與預(yù)測(cè)上[14]。在醫(yī)學(xué)領(lǐng)域，EMD方法目前已經(jīng)應(yīng)用于心音、呼吸音的識(shí)別與處理上[15-16]。但是在傳染病領(lǐng)域，目前此方面的應(yīng)用尚屬空白。以往關(guān)于傳染病發(fā)病數(shù)據(jù)的時(shí)間分布研究，總是習(xí)慣于使用ARIMA模型對(duì)時(shí)間序列數(shù)據(jù)進(jìn)行處理，但是ARIMA模型需要基于平穩(wěn)數(shù)據(jù)，對(duì)于非平穩(wěn)的數(shù)據(jù)，則必須對(duì)原始的時(shí)間序列進(jìn)行差分，以獲取平穩(wěn)數(shù)據(jù)，但是差分則意味著信息的損失，而傳染病的發(fā)病數(shù)據(jù)，受到多種因素的影響，難以達(dá)到一階平穩(wěn)。因此，對(duì)于非平穩(wěn)的傳染病時(shí)間序列數(shù)據(jù)，EMD是一個(gè)更好的選擇。

本文選取了過去13年間我國艾滋病的發(fā)病例數(shù)作為研究數(shù)據(jù)，首先對(duì)艾滋病的發(fā)病趨勢(shì)進(jìn)行分解，在去除原始數(shù)據(jù)中的噪音后，所得的殘差曲線能夠準(zhǔn)確地描述艾滋病既往的發(fā)病趨勢(shì)，該結(jié)果與以往的艾滋病流行病學(xué)研究結(jié)果相吻合[17-18]。隨后運(yùn)用BP神經(jīng)網(wǎng)絡(luò)基于既往疾病數(shù)據(jù)進(jìn)行模型構(gòu)建，在此基礎(chǔ)上進(jìn)一步對(duì)2016年的疫情進(jìn)行預(yù)測(cè)，RMSE的結(jié)果表明，BP神經(jīng)網(wǎng)絡(luò)對(duì)于傳染病的時(shí)間序列數(shù)據(jù)，其預(yù)測(cè)效果明顯優(yōu)于傳統(tǒng)的SARIMA模型。

對(duì)于非平穩(wěn)的時(shí)間序列，BP神經(jīng)網(wǎng)絡(luò)預(yù)測(cè)有著明顯優(yōu)勢(shì)，但是對(duì)于影響因素眾多的復(fù)雜數(shù)據(jù)，預(yù)測(cè)結(jié)果受到多種因素的干擾，給精準(zhǔn)預(yù)測(cè)帶來了一定的困難。與EMD結(jié)合使用后，利用EMD能夠分離各個(gè)不同因素相互干擾的特性，再通過對(duì)各個(gè)分量進(jìn)行預(yù)測(cè)再合成的方法，提高了預(yù)測(cè)的精度。對(duì)于傳染病監(jiān)測(cè)獲得的時(shí)間序列數(shù)據(jù)，使用EMD結(jié)合BP神經(jīng)網(wǎng)絡(luò)的方法，既能夠清楚地識(shí)別傳染病的發(fā)病趨勢(shì)，也能精準(zhǔn)地對(duì)未來的發(fā)病情況進(jìn)行預(yù)測(cè)。但是，基于時(shí)間序列數(shù)據(jù)的預(yù)測(cè)僅僅是預(yù)測(cè)數(shù)據(jù)本身，因此獲得的R2并不高。如果其他影響因素的數(shù)據(jù)可及的話，應(yīng)用此法，可以精準(zhǔn)地預(yù)測(cè)出未來幾年內(nèi)傳染病的發(fā)病趨勢(shì)，從而準(zhǔn)確估計(jì)該傳染病帶來的疾病負(fù)擔(dān)，進(jìn)一步為衛(wèi)生政策的制定提供理論依據(jù)，將更加有利于傳染病的預(yù)防控制。

[1] 魏珊,陸一涵,高眉揚(yáng),等.我國戊型肝炎發(fā)病例數(shù)的時(shí)間序列分析.中國衛(wèi)生統(tǒng)計(jì),2012,29(6):808-811.

[2] 舒國通.以不同人群為基礎(chǔ)的傳染病影響因素流行病學(xué)研究.浙江大學(xué),2004.

[3] Wu ZH.A study of the characteristics of white noise using the empirical mode decomposition method.The Royal Society,2004,3(10):1597-1611.

[4] 李婉華,陳宏,郭昆,等.基于隨機(jī)森林算法的用電負(fù)荷預(yù)測(cè)研究.計(jì)算機(jī)工程與應(yīng)用,2016,52(23):236-243.

[5] 黃建國,羅航,王厚軍,等.運(yùn)用GA-BP神經(jīng)網(wǎng)絡(luò)研究時(shí)間序列的預(yù)測(cè).電子科技大學(xué)學(xué)報(bào),2009,38(5):687-692.

[6] Huang NE.A new method for nonlinear and nonstationary time series analysis:Empirical mode decomposition and Hilbert spectral analysis.//Szu HH,Vetterli M,Campbell WJ,et al.Wavelet Applications Vii.City,2000:197-209.://WOS:000088256200017.

[7] Wang SZ,Wu N,Wang L.Wind speed forecasting based on the hybird ensemble empirical mode decomposition and GA-BP neural network method.Renew Energ,2016,(94):629-636.

[8] Zhang Z.A gentle introduction to artificial neural networks.Ann Transl Med,2016,4(19):370.

[9] Huang D,Wu Z.Forecasting outpatient visits using empirical mode decomposition coupled with back-propagation artificial neural networks optimized by particle swarm optimization.PLoS One,2017,12(2):e0172539.

[10]NE Huang ZS,SR Long,ML Wu,et al.The empirical mode decomposition and Hilbert spectrum for nonlinear and nonstationary time series analysis.Proceedings of the Royal Society London A,1998,(454):903-995.

[11]Zeng MH.A simple boundary process technique for empirical mode decomposition.Proceedings of 2004 IEEE International Geoscience and Remote Sensing Symposium,2004,(6):4258-4261.

[12]趙肖宇,方一鳴,譚峰,等.EMD時(shí)頻分析拉曼光譜和近紅外光譜.光譜學(xué)與光譜分析,2016,36(2):424-429.

[13]王文波,費(fèi)浦生,羿旭明.基于EMD與神經(jīng)網(wǎng)絡(luò)的中國股票市場(chǎng)預(yù)測(cè).系統(tǒng)工程理論與實(shí)踐,2010,30(6):1027-1033.

[14]Safi K,Hutin E,Mohammed S,et al.Human static postures analysis using empirical mode decomposition.Conf Proc IEEE Eng Med Biol Soc,2016,2016:3765-3768.

[15]陳新華,成謝鋒.一種改進(jìn)型綜合去噪算法在心音信號(hào)預(yù)處理上的研究.南京郵電大學(xué)學(xué)報(bào)(自然科學(xué)版),2010,30(6):96-100,106.

[16]鐘思軍,賴克方,杜明輝,等.經(jīng)驗(yàn)?zāi)B(tài)分解在咳嗽音檢測(cè)中的應(yīng)用.北京生物醫(yī)學(xué)工程,2008,27(3):238-240,299.

[17]Zhang L,Chow EP,Jing J,et al.HIV prevalence in China:integration of surveillance data and a systematic review.Lancet Infect Dis,2013,13(11):955-963.

[18]GBD 2013 Mortality and Causes of Death Collaborators.Global,regional,and national age-sex specific all-cause and cause-specific mortality for 240 causes of death,1990-2013:a systematic analysis for the Global Burden of Disease Study 2013.Lancet,2015,385(9963):117-171.