EMD-BP神經(jīng)網(wǎng)絡(luò)在傳染病發(fā)病趨勢和預(yù)測研究中的應(yīng)用*

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

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

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

原理和方法

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

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

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

(2)下包絡(luò)線：找到原始數(shù)據(jù)X(t)的所有極小值點，同樣采用三次樣條插值函數(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)篩選出來的第一個IMF表示為c1，可以先從原始數(shù)據(jù)中分離出來，即：r1=X(t)-c1

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

直到滿足以下兩個條件之一，分解結(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)作為一種機器學(xué)習(xí)方法，具有很強的學(xué)習(xí)能力和適應(yīng)能力，在非線性系統(tǒng)的預(yù)測研究中得到廣泛的應(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]。

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

3.統(tǒng)計方法

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

實例分析

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

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

2.EMD包簡介

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

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

(2)tt：觀察的時間。

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

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

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

(6)boundary：邊界調(diào)整。為了消除信號的邊界效應(yīng),必須進(jìn)行邊界調(diào)整。Huang通過在信號兩邊重復(fù)添加人為波段來擴展原有信號，添加部分通過靠近邊界的極值形成的固有模式構(gòu)建[10]。Zeng和He提出兩種擴展信號的方式,它們分別通過對稱和自反方式添加信號，被稱為偶擴展和奇擴展[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ù)就可以實現(xiàn) EMD，EMD包也提供了emd函數(shù)直接實現(xiàn)整個過程。

3.EMD結(jié)果

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

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

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

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

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

以2004-2015的數(shù)據(jù)作為訓(xùn)練集，2016年的數(shù)據(jù)作為測試集。在R語言caret包中對EMD分解出的四個IMF以及殘差分別進(jìn)行神經(jīng)網(wǎng)絡(luò)模型訓(xùn)練和預(yù)測，并將IMF的預(yù)測結(jié)果和殘差的預(yù)測結(jié)果進(jìn)行合成，得到原始數(shù)據(jù)的預(yù)測結(jié)果。預(yù)測結(jié)果見圖3，圖中黑色曲線為原始數(shù)據(jù)，依次是對2016年的原始數(shù)據(jù)進(jìn)行EMD后得到的4個IMF和殘差，紅色曲線BPANN模型對以上五組數(shù)據(jù)進(jìn)行預(yù)測所得到的預(yù)測值。此外，SARIMA(2,1,1)(2,1,1)12模型同樣被用來對2016年的原始數(shù)據(jù)進(jìn)行預(yù)測。進(jìn)行歸一化后的五個神經(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)元個數(shù)，數(shù)值越大，模型越復(fù)雜。RMSE：均方根誤差。

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

討 論

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

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

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

[1] 魏珊,陸一涵,高眉揚,等.我國戊型肝炎發(fā)病例數(shù)的時間序列分析.中國衛(wèi)生統(tǒng)計,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] 李婉華,陳宏,郭昆,等.基于隨機森林算法的用電負(fù)荷預(yù)測研究.計算機工程與應(yīng)用,2016,52(23):236-243.

[5] 黃建國,羅航,王厚軍,等.運用GA-BP神經(jīng)網(wǎng)絡(luò)研究時間序列的預(yù)測.電子科技大學(xué)學(xué)報,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時頻分析拉曼光譜和近紅外光譜.光譜學(xué)與光譜分析,2016,36(2):424-429.

[13]王文波,費浦生,羿旭明.基于EMD與神經(jīng)網(wǎng)絡(luò)的中國股票市場預(yù)測.系統(tǒng)工程理論與實踐,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)型綜合去噪算法在心音信號預(yù)處理上的研究.南京郵電大學(xué)學(xué)報(自然科學(xué)版),2010,30(6):96-100,106.

[16]鐘思軍,賴克方,杜明輝,等.經(jīng)驗?zāi)B(tài)分解在咳嗽音檢測中的應(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.