基于長(zhǎng)短時(shí)記憶神經(jīng)網(wǎng)絡(luò)的手足口病發(fā)病趨勢(shì)預(yù)測(cè)

馬停停，冀天嬌，楊冠羽*，陳 陽(yáng)，許文波，劉宏圖

（1.東南大學(xué)計(jì)算機(jī)科學(xué)與工程學(xué)院，南京 210096；2.中國(guó)疾病預(yù)防控制中心病毒病預(yù)防控制所衛(wèi)生部醫(yī)學(xué)病毒學(xué)和病毒病重點(diǎn)實(shí)驗(yàn)室，北京102206）（*通信作者電子郵箱yang.list@seu.edu.cn）

0 引言

手足口病（Hand-Foot-Mouth Disease，HFMD）是由多種腸道病毒引起的一種常見(jiàn)傳染病［1］，2008 年5 月，衛(wèi)生部將手足口病納入丙類(lèi)傳染病管理，開(kāi)始網(wǎng)絡(luò)直報(bào)。手足口病具有散播快、流行性強(qiáng)的特點(diǎn)，并可在短時(shí)間內(nèi)造成一定規(guī)模的流行，成人以隱性感染為主，常表現(xiàn)為無(wú)癥狀病毒攜帶者，多發(fā)于5 歲以下嬰幼兒［2］。其癥狀主要表現(xiàn)為發(fā)熱和手、足、口等部位的皮疹、皰疹，伴或不伴口腔潰瘍，病情嚴(yán)重者可引起心肌炎、肺水腫、無(wú)菌性腦膜腦炎等致命性并發(fā)癥，重癥病例死亡率較高，危害嚴(yán)重［3］。研究手足口病的發(fā)病特點(diǎn)以及傳播規(guī)律，并建立相應(yīng)的時(shí)間序列模型來(lái)預(yù)測(cè)發(fā)病趨勢(shì)是很有必要的，能夠?yàn)橹贫A(yù)防措施提供科學(xué)依據(jù)。

近年來(lái)越來(lái)越多的學(xué)者從事手足口病病原體分析［4］，其與氣象因素的相關(guān)性研究以及流行趨勢(shì)預(yù)測(cè)等方面的工作。Nguyen 等［5］研究了環(huán)境因素與手足口病發(fā)病的關(guān)系，一些研究表明，手足口病的季節(jié)性變化與降雨量、相對(duì)濕度、溫度和氣壓有關(guān)［6-8］。雖然許多研究報(bào)告了平均溫度或濕度與手足口病之間有顯著相關(guān)性，但這相關(guān)性沒(méi)有普適性，在其他地方?jīng)]有顯著性［9-10］。不同地區(qū)現(xiàn)有調(diào)查結(jié)果的不一致可歸因于不同地理區(qū)域的天氣條件，人口統(tǒng)計(jì)特征和衛(wèi)生條件的多樣性［11］。在預(yù)測(cè)方面，通常使用差分自回歸滑動(dòng)平均（AutoRegressive Integrated Moving Average，ARIMA）模型［12］或季節(jié)性差分自回歸滑動(dòng)平均（Seasonal AutoRegressive Integrated Moving Average，SARIMA）模型［13］，并對(duì)未來(lái)的發(fā)病情況進(jìn)行預(yù)測(cè)。Pons-salort 等使用傳播動(dòng)力學(xué)模型（Susceptible Infectious Recovered，SIR）結(jié)合基于粒子濾波的最大似然推理框架對(duì)日本20 種常見(jiàn)的病毒血清型進(jìn)行建模并進(jìn)行2 年的預(yù)測(cè)，結(jié)果較好，發(fā)現(xiàn)人類(lèi)腸道病毒引起的疾病發(fā)病率與血清型特異性免疫有關(guān)［13］。但這些方法有一定局限性，ARIMA 模型不能對(duì)非線性關(guān)系進(jìn)行建模，SIR 模型不能加入氣象因素，無(wú)法充分利用多維度輸入數(shù)據(jù)中信息。本文使用長(zhǎng)短時(shí)記憶（Long Short-Term Memory，LSTM）網(wǎng)絡(luò)模型能夠有效解決這兩個(gè)問(wèn)題。

1 數(shù)據(jù)處理

1.1 數(shù)據(jù)來(lái)源

濟(jì)南市和廣州市的手足口病電子病例數(shù)據(jù)來(lái)源于中國(guó)疾病預(yù)防控制中心（Chinese Center for Disease Control and Prevention，CDC），電子病例數(shù)據(jù)包括病人年齡、所在地區(qū)及確診日期（對(duì)患者的個(gè)人信息已進(jìn)行脫敏處理），癥狀嚴(yán)重的患者，對(duì)其致病的病原體進(jìn)行檢測(cè)，并分為三類(lèi)EV71、CVA16和其他腸道病毒。數(shù)據(jù)范圍2010 年1 月1 日至2018 年12 月31 日。本文選取濟(jì)南市和廣州市手足口病進(jìn)行研究，兩市的發(fā)病模式不同，具有代表性，我國(guó)其他地區(qū)發(fā)病情況與兩市相近。

相關(guān)的氣象數(shù)據(jù)來(lái)源于中國(guó)氣象科學(xué)數(shù)據(jù)共享服務(wù)網(wǎng)，包括濟(jì)南市和廣州市每日的平均氣溫和相對(duì)濕度。

1.2 統(tǒng)計(jì)分析

濟(jì)南市是山東省省會(huì)，地處中國(guó)華東地區(qū)，位于北緯36°40′，東經(jīng)117°00′，常住人口746.04 萬(wàn)人。據(jù)統(tǒng)計(jì)，2010—2018 年濟(jì)南市共報(bào)告HFMD 111 638 例，男女發(fā)病比為1.46，發(fā)病年齡主要集中在0～6 歲（95.45%），以散居兒童為主（61.77%）。

廣州市是廣東省省會(huì)，地處中國(guó)南部，位于北緯23°06′，東經(jīng)113°15′，常住人口1 490.44萬(wàn)人。2010—2018年廣州市共報(bào)告HFMD 542 759例，男女發(fā)病比為1.6，發(fā)病年齡主要集中在0～6歲（96.66%），以散居兒童為主（75.46%）。

對(duì)于每日發(fā)病數(shù)據(jù)，以周為單位進(jìn)行累計(jì)。對(duì)于氣象數(shù)據(jù)，以周為單位取平均值。圖1 為2010—2018 年廣州市與濟(jì)南市每周發(fā)病數(shù)統(tǒng)計(jì)。兩市手足口病呈現(xiàn)不同的發(fā)病模式，濟(jì)南市全年發(fā)病曲線呈現(xiàn)一年一個(gè)高峰期，峰值出現(xiàn)在5—7月份；而廣州市全年發(fā)病曲線呈現(xiàn)一年兩個(gè)高峰期，第一個(gè)高峰期峰值在5—6月份，第二個(gè)高峰期峰值在8—9月份。

表1 是濟(jì)南市和廣州市氣象因素的統(tǒng)計(jì)值。濟(jì)南市全年平均氣溫15℃左右，廣州市全年平均氣溫22℃左右，廣州市平均相對(duì)濕度高于濟(jì)南市。

圖1 濟(jì)南市和廣州市2010—2018年每周發(fā)病數(shù)統(tǒng)計(jì)Fig.1 Weekly incidence statistics of Jinan and Guangzhou from 2010 to 2018

表1 濟(jì)南市和廣州市2010—2018年氣象因素統(tǒng)計(jì)Tab.1 Statistics of meteorological factors in Jinan and Guangzhou from 2010 to 2018

1.3 相關(guān)性分析

氣象因素對(duì)傳染病的發(fā)生、傳播起著至關(guān)重要的作用，其不僅影響個(gè)體的自身免疫能力，還會(huì)影響病原微生物的繁殖和傳播能力。HFMD 作為一種腸道病毒導(dǎo)致的傳染病，適宜的溫度和濕度必然會(huì)導(dǎo)致腸道病毒在外環(huán)境中大量繁殖并提高病毒的存活能力。各氣象因素與手足口病發(fā)病數(shù)的秩相關(guān)系數(shù)（又稱(chēng)Spearman 相關(guān)系數(shù)）如表2 所示，各相關(guān)系數(shù)顯著水平p值均小于0.05，具有統(tǒng)計(jì)學(xué)意義。結(jié)果表明平均氣溫與發(fā)病數(shù)存在強(qiáng)相關(guān)性，平均相對(duì)濕度與發(fā)病數(shù)存在弱相關(guān)性。

表2 氣象因素與手足口病秩相關(guān)系數(shù)Tab.2 Rank correlation coefficients between meteorological factors and HFMD

2 網(wǎng)絡(luò)模型

本章介紹了如何將預(yù)測(cè)問(wèn)題轉(zhuǎn)化為監(jiān)督學(xué)習(xí)問(wèn)題和LSTM模型的原理，以及模型多步預(yù)測(cè)的方式。

2.1 問(wèn)題轉(zhuǎn)化

本文使用滑動(dòng)窗口將預(yù)測(cè)問(wèn)題作為監(jiān)督學(xué)習(xí)問(wèn)題進(jìn)行研究。手足口病的周發(fā)病數(shù)按時(shí)間排列可看作時(shí)間序列，用s=[s[0]，s[1]，…，s[T]]表示，每周的氣象數(shù)據(jù)與之對(duì)應(yīng)。在建模預(yù)測(cè)中，需要將數(shù)據(jù)作為向量傳送給模型，該向量由固定時(shí)間步長(zhǎng)的數(shù)據(jù)組成，可由序列s使用滑動(dòng)窗口的方式獲得，設(shè)滑動(dòng)窗口的大小為nt。給定固定長(zhǎng)度的過(guò)去值，建立模型的目的是預(yù)測(cè)時(shí)間序列未來(lái)的n0個(gè)值。這樣，給定離散時(shí)間t的輸入矢量定義為xt=[s[t-nt+1]，s[t-nt+2]，…，s[t]]，需要做的是推斷接下來(lái)的n0周的發(fā)病數(shù)，可以將輸出定義為yt=[s[t+1]，s[t+2]，…，s[t+n0]]。

2.2 LSTM模型

LSTM 神經(jīng)網(wǎng)絡(luò)最早由Hochreiter 等提出，并由Graves 進(jìn)行改進(jìn)［15］，是基于RNN 的一種完善，解決RNN 中易出現(xiàn)的梯度消亡問(wèn)題。LSTM單元結(jié)構(gòu)［16］如圖2所示。

圖2 LSTM單元結(jié)構(gòu)Fig.2 Structure of LSTM cell

圖2中存在3個(gè)控制門(mén)，分別是輸入門(mén)、輸出門(mén)和遺忘門(mén)，3個(gè)門(mén)的輸出分別連接到1個(gè)乘法單元上，從而分別控制網(wǎng)絡(luò)的輸入、輸出以及存儲(chǔ)單元的狀態(tài)。LSTM用兩個(gè)門(mén)來(lái)控制單元狀態(tài)c的內(nèi)容，遺忘門(mén)決定了上一時(shí)刻的單元狀態(tài)ct-1有多少保留到當(dāng)前時(shí)刻ct。另一個(gè)是輸入門(mén)，它決定了當(dāng)前時(shí)刻的網(wǎng)絡(luò)的輸入xt有多少保存到單元狀態(tài)ct。LSTM 使用輸出門(mén)來(lái)控制單元狀態(tài)ct有多少輸出到LSTM 的當(dāng)前輸出值ht。LSTM單元更新公式如下。

遺忘門(mén)的公式為：

其中：Wfx和Wfh是遺忘門(mén)的權(quán)重矩陣，bf為偏置項(xiàng)。

輸出門(mén)為：

其中：Wix和Wih是輸入門(mén)的權(quán)重矩陣，bi為偏置項(xiàng)。

當(dāng)前時(shí)刻的單元狀態(tài)由上一次的單元狀態(tài)按元素乘以遺忘門(mén)，再加上當(dāng)前輸入的單元狀態(tài)按元素乘以輸入門(mén)，符號(hào)°表示按元素乘。

輸出門(mén)控制了長(zhǎng)期記憶對(duì)當(dāng)前輸出的影響：

LSTM最終的輸出，由輸出門(mén)和單元狀態(tài)共同確定：

2.3 預(yù)測(cè)方式

大多數(shù)對(duì)手足口病發(fā)病趨勢(shì)預(yù)測(cè)為單點(diǎn)預(yù)測(cè)，即預(yù)測(cè)下一周的發(fā)病數(shù)，會(huì)有較好的預(yù)測(cè)精度，但這對(duì)于制定預(yù)防策略意義不是很大。通過(guò)迭代預(yù)測(cè)的方式可以得到更長(zhǎng)的預(yù)測(cè)結(jié)果。由于迭代預(yù)測(cè)將預(yù)測(cè)值當(dāng)作真實(shí)值代入模型，會(huì)造成誤差累積，隨著迭代次數(shù)增加，誤差會(huì)越來(lái)越大。在長(zhǎng)期預(yù)測(cè)與誤差累積之間做一個(gè)平衡，選擇預(yù)測(cè)后12 周的發(fā)病。LSTM的步長(zhǎng)設(shè)置會(huì)影響預(yù)測(cè)的精度，迭代預(yù)測(cè)時(shí)輸入向量中不斷補(bǔ)充預(yù)測(cè)值，導(dǎo)致真實(shí)的值數(shù)量下降，為了保證輸入向量中至少有一半的真實(shí)值，LSTM 的步長(zhǎng)設(shè)置為26 周。由于氣象因素對(duì)手足口病的影響存在一定的滯后性，且手足口病有3～7天的潛伏期，氣象因素采用前兩周的數(shù)值。將手足口病的發(fā)病序列按照2.1 節(jié)的方式轉(zhuǎn)換成輸入序列，訓(xùn)練單個(gè)模型以進(jìn)行預(yù)測(cè)下一周的發(fā)病數(shù)。在預(yù)測(cè)階段，將預(yù)測(cè)的輸出當(dāng)作正確的輸出，作為輸入的一部分遞歸反饋，重復(fù)11次此過(guò)程。

3 實(shí)驗(yàn)結(jié)果

3.1 模型實(shí)現(xiàn)

溫度和相對(duì)濕度為與氣候相關(guān)的數(shù)據(jù)類(lèi)別，與發(fā)病數(shù)據(jù)類(lèi)別不同，將其加入到模型中，需要確保不會(huì)與發(fā)病數(shù)據(jù)互相干擾，又要保證加入的氣象數(shù)據(jù)時(shí)間跨度適當(dāng)，本文選擇的方式如圖3 所示，發(fā)病數(shù)據(jù)輸入到LSTM 網(wǎng)絡(luò)中，以捕獲手足口病數(shù)據(jù)的長(zhǎng)期序列屬性，將LSTM 層的輸出結(jié)果與氣象數(shù)據(jù)進(jìn)行連接，最后通過(guò)全連接層輸出最終的預(yù)測(cè)結(jié)果。本實(shí)驗(yàn)中為氣象因素，時(shí)間跨度為兩周，每次迭代都更新。由于預(yù)測(cè)的范圍為12 周，使用氣象數(shù)據(jù)作為輸入迭代時(shí)，不能用到未來(lái)的數(shù)據(jù)，測(cè)試集2018年的氣象數(shù)據(jù)使用的是2015—2017年對(duì)應(yīng)時(shí)期的均值。

圖3 模型網(wǎng)絡(luò)結(jié)構(gòu)Fig.3 Network structure of the model

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

模型的評(píng)估使用均方誤差平方根（Rooted Mean Squared Error，RMSE）和平均絕對(duì)誤差（Mean Absolute Error，MAE）衡量。RMSE 和MAE 的數(shù)值越小，說(shuō)明模型的預(yù)測(cè)效果越好。RMSE和MAE的計(jì)算公式如下：

其中：xi為第i周實(shí)際發(fā)病數(shù)為模型預(yù)測(cè)值，n為序列的總樣本數(shù)。

3.3 實(shí)驗(yàn)與結(jié)果

在實(shí)驗(yàn)過(guò)程中，LSTM 模型用python 和tensorflow 框架實(shí)現(xiàn)，并使用GPU 進(jìn)行加速。采用自適應(yīng)學(xué)習(xí)率并將初始值設(shè)置為0.001，為了防止訓(xùn)練集過(guò)擬合，訓(xùn)練中采用L2正則化。

表3為不同模型的第12周預(yù)測(cè)結(jié)果對(duì)比，LSTM為不加氣象因素建立的模型，LSTMT 表示模型中加入溫度，LSTMR 表示模型中加入濕度，LSTMTR 表示模型中同時(shí)加入溫度和濕度。在分別加入溫度和平均濕度之后，模型的預(yù)測(cè)結(jié)果有所提升，其中加入平均溫度之后，模型預(yù)測(cè)精度提升較大。濟(jì)南市數(shù)據(jù)集上，在同時(shí)加入平均溫度和平均相對(duì)濕度時(shí)，結(jié)果最好，MAE 為74.9，RMSE 為128.3。廣州市數(shù)據(jù)集上，同時(shí)加入平均溫度和平均相對(duì)濕度時(shí)，效果不如僅加入平均溫度，但好于只加入平均相對(duì)濕度，最好結(jié)果MAE 為427.7，RMSE為604。

表3 LSTM加入不同氣象因素第12周預(yù)測(cè)結(jié)果對(duì)比Tab.3 Comparison of prediction results of the 12th week by LSTM model with different meteorological factors

對(duì)于濟(jì)南市數(shù)據(jù)集，各個(gè)模型預(yù)測(cè)結(jié)果與真實(shí)發(fā)病數(shù)曲線大致相同（如圖4），在加入氣象因素之后，可以看到模型的峰值有一定的后移，這與2018 年真實(shí)情況更為接近，同時(shí)加入溫度和濕度后效果最好，這也說(shuō)明了溫度與濕度對(duì)手足口病發(fā)病的影響；對(duì)于廣州市數(shù)據(jù)集，由于在2018 年發(fā)病人數(shù)驟降，各個(gè)模型預(yù)測(cè)峰值均高于實(shí)際值（如圖5）。相比較而言，在加入氣象因素之后，預(yù)測(cè)的峰值更接近真實(shí)值，只加入平均溫度效果最好。

圖4 4種模型在濟(jì)南市2018年的預(yù)測(cè)結(jié)果Fig.4 Prediction results of four models for Jinan in 2018

為驗(yàn)證本文模型的有效性，將本文結(jié)果與常用手足口病預(yù)測(cè)模型的結(jié)果進(jìn)行對(duì)比，如表4。表4 中SARIMA 模型和支持向量回歸（Support Vector Regression，SVR）模型的預(yù)測(cè)方式與使用LSTM模型相同。

在濟(jì)南市數(shù)據(jù)集上，使用SARIMA 模型結(jié)果略好于LSTM模型，但不如LSTM 模型加入氣象因素。SARIMA，模型依賴(lài)年度周期和前期歷史數(shù)據(jù)，當(dāng)數(shù)據(jù)集規(guī)律性較好時(shí)，如濟(jì)南市手足口病發(fā)病趨勢(shì)，結(jié)果較好。當(dāng)數(shù)據(jù)集規(guī)律性沒(méi)那么強(qiáng)時(shí)，如廣州市手足口病發(fā)病趨勢(shì)，預(yù)測(cè)結(jié)果會(huì)比較差。而LSTM模型更具有普遍性。

圖5 4種模型在廣州市2018年的預(yù)測(cè)結(jié)果對(duì)比Fig.5 Prediction results of four models for Guangzhou in 2018

表4 不同模型的第12周預(yù)測(cè)結(jié)果對(duì)比Tab.4 Comparison of prediction results of the 12th week of different models

4 結(jié)語(yǔ)

本文主要使用LSTM 模型結(jié)合溫度和濕度預(yù)測(cè)濟(jì)南市和廣州市手足口病發(fā)病趨勢(shì)，加入氣象因素能使預(yù)測(cè)結(jié)果更為準(zhǔn)確，這也說(shuō)明了氣象因素對(duì)手足口病發(fā)病有一定的影響。

本文是第一個(gè)應(yīng)用LSTM 方法并結(jié)合氣象因素對(duì)手足口病發(fā)病中期預(yù)測(cè)的，能夠?yàn)槭肿憧诓》揽靥峁﹨⒖家罁?jù)，同時(shí)為今后對(duì)手足口病發(fā)病趨勢(shì)預(yù)測(cè)的研究提供了相關(guān)理論基礎(chǔ)。但同時(shí)本研究依舊有很多待改進(jìn)的地方，手足口病由不同的腸道病毒引起，這些病毒有不同的流行周期，若對(duì)此進(jìn)行深入分析，納入模型，或許會(huì)有更好的預(yù)測(cè)效果。此外，傳染病的傳播途徑，易感人群的保護(hù)措施等通常會(huì)受到所在地區(qū)的經(jīng)濟(jì)狀況等多種社會(huì)因素的影響，這些混雜因素也尚未考慮。涉及這些方法的全面分析將有助于了解和監(jiān)測(cè)手足口病的傳播，減少感染的危險(xiǎn)，并保護(hù)兒童免受這種疾病的傷害。