融合ARIMA和RBFNN 的血糖預(yù)測

王延年，雍永強，賈曉燦，李全忠，孫長青＋

（1.鄭州大學(xué) 信息工程學(xué)院，河南 鄭州450001；2.鄭州大學(xué) 公共衛(wèi)生學(xué)院，河南 鄭州450001；3.河南省人民醫(yī)院 內(nèi)分泌科，河南 鄭州450003）

0 引 言

CGMS （continous glucose monitoring system）借助 埋在臍周皮下的探頭檢測葡萄糖濃度，能連續(xù)獲得患者血糖的變化及波動。國內(nèi)外學(xué)者基于CGMS所提供的數(shù)據(jù)，提出了多種預(yù)測方法，常見的有自回歸 （auto－regressive，AR）模型、自回歸移動平均 （auto－regressive moving average，ARMA）模型［3］、支持向量機［4］、極限學(xué)習(xí)機［5］、神經(jīng)網(wǎng)絡(luò)［6］等。這些方法都只利用CGMS所提供的數(shù)據(jù)，沒有考慮外部因素 （飲食、藥物注射、運動等）對人體血糖的影響。

S.Pappada等［7］提出的神經(jīng)網(wǎng)絡(luò)方法和Eleni I.Georga等［8］提出的支持向量回歸方法不僅采用CGMS所提供的數(shù)據(jù)，而且結(jié)合了人體的生理模型。結(jié)合生理模型的預(yù)測方法考慮比較全面，但需要大量病理學(xué)和生理學(xué)的知識，所提算法和預(yù)測模型復(fù)雜龐大，預(yù)測有一定的延遲。

本文基于CGMS 所提供的血糖數(shù)據(jù)，提出一種融合ARIMA 和RBF神經(jīng)網(wǎng)絡(luò)模型預(yù)測血糖的綜合預(yù)測方法。針對外部事件對血糖的影響，提出一種發(fā)現(xiàn)和處理算法，在血糖變化受外部因素影響時能自動調(diào)整未來一段時間內(nèi)的血糖預(yù)測值，保障預(yù)測精度。

1 方 法

采用ARIMA 模型對未來血糖值進行預(yù)測，并計算出誤差項，用RBFNN 算法對輸入值、趨勢值、預(yù)測殘差、預(yù)測殘差的趨勢值等進行學(xué)習(xí)、網(wǎng)絡(luò)訓(xùn)練與擬合，最后將ARIMA 計算出的預(yù)測值與RBF 算法得到的修正值組合得到準(zhǔn)確的預(yù)測。同時針對外部因素對血糖變化的影響，提出一種發(fā)現(xiàn)和處理算法。

1.1 組合預(yù)測模型

融合ARIMA 和RBFNN 的組合預(yù)測模型的具體結(jié)構(gòu)如圖1所示。

圖1 融合ARIMA 和RBFNN 的組合預(yù)測模型

基于組合預(yù)測模型的血糖預(yù)測步驟為：

（1）用ARIMA 模型對從CGMS 采集到的血糖序列y（t）進行建模預(yù)測。PH 時間段后的預(yù)測結(jié)果為（t＋PH｜t）。

（2）利用 （1）得到結(jié)果得到ARIMA 產(chǎn)生的預(yù)測殘差：e（t）＝y(tǒng)（t）－（t｜t－PH）。

（3）將e（T），（1－z－Tm）e（t），y（t），（1－z－Tm）y（t）作為RBFNN 的輸入，得到PH 時間段后的預(yù)測結(jié)果（t＋PH ｜t）。

1.2 數(shù)據(jù)預(yù)處理

從CGMS采集的血糖數(shù)據(jù)由于存在噪聲干擾需要進行平滑濾波處理?？柭鼮V波具有計算量小、實時性高、效果好等優(yōu)點，而且采用Kalman濾波對CGMS采集的數(shù)據(jù)進行平滑處理能有效地減少信號中的噪聲，同時減少預(yù)測值和實際值之間的時延［9］。本文采用Kalman濾波對CGMS采集的數(shù)據(jù)進行平滑處理。

1.3 ARIMA模型

ARIMA 的一般表達式為

記作ARIMA（p，d，q），其中，xt是血糖值組成的時間序列，i 是自回歸 （AutoregRessive，AR）項的參數(shù)，θj是滑動平均 （MovingAverage，MA）項的參數(shù)，ω為零均值及方差為δ2的高斯白噪聲。用于對數(shù)據(jù)進行差分處理，可表示為

ARIMA 處理過程如下：

（1）由于血糖序列是一種非平穩(wěn)時間序列，首先對它進行差分處理，即

那么血糖序列的ARIMA （p，d，q）模型為

其中，｛αt｝為白噪聲序列。

（2）模型識別。p，d，q 的選擇是ARIMA 建模的關(guān)鍵。首先對血糖數(shù)據(jù)進行差分處理，每次差分比較其自相關(guān)（auto correlation function，ACF）和偏自相關(guān) （partial auto correlation function，PACF）組成的條形圖，選擇最佳的d。p，q的確定采用AIC （Akaiki information criterion）準(zhǔn)則和BIC （Bayesian information criterion）準(zhǔn)則［10］，計算AIC和BIC的方程如下［11］

式中：k——參數(shù)模型的個數(shù)，n——樣本數(shù)。

（3）估計參數(shù)并檢驗?zāi)Ｐ?。模型中的所有參?shù)的估計都采用極大似然估計法，然后進行檢驗。方法是檢驗誤差序列是否具有隨機性；如果檢驗通不過就繼續(xù)進行識別。

（4）利用選中的最佳參數(shù)組成預(yù)測模型。ARIMA 模型采用已經(jīng)采集到的在xt，xt－1，…xt－n時刻的血糖值，預(yù)測未來xt＋m時刻的血糖值。方程為xnn＋m＝E［xn＋m｜xn，xn－1，…，x1］。詳細(xì)的過程可參考文獻 ［12］。

1.4 RBF神經(jīng)網(wǎng)絡(luò)

對于線性時間序列，ARIMA 的預(yù)測簡單精準(zhǔn)。從CGMS收集到的血糖數(shù)據(jù)是線性的，但影響血糖變化的因素復(fù)雜多樣。只使用線性預(yù)測模型會帶來時延，影響預(yù)測的準(zhǔn)確性。

RBF神經(jīng)網(wǎng)絡(luò)在處理非線性問題時能力很強。不僅學(xué)習(xí)速度快，且避免局部最小問題，只有少量權(quán)值需要調(diào)整。

RBF神經(jīng)網(wǎng)絡(luò)學(xué)習(xí)過程：

（1）利用高斯函數(shù)激勵函數(shù)，輸入和輸出的非線性映射如下

其中，xp＝（xp1，xp2，…xpn）T為第p個輸入樣本，ci為隱含層高斯函數(shù)的中心，σ為高斯函數(shù)的方差，為歐式范數(shù)。

（2）隱含層和輸出層之間的線性映射，函數(shù)為

式中：h——隱含層的節(jié)點總數(shù)，ωij——隱含層到輸出層的線性加權(quán)值，yi——網(wǎng)絡(luò)中第j 個輸出節(jié)點的實際輸出，n——輸出節(jié)點數(shù)。d為期望輸出值，隱含層的σ可表示為

（3）參照樣本數(shù)據(jù)，校正輸出層和隱含層之間的參數(shù)。本文中RBF的輸入包括4個部分：

（1）ARIMA 產(chǎn)生的預(yù)測殘差

（2）ARIMA 預(yù)測殘差的趨勢（1－z－Tm）e（t）。

（3）目前的血糖值y（t），即從CGMS讀取的數(shù)據(jù)。

（4）在時間步Tm目前的血糖趨勢值（1－z－Tm）y（t）。

2 預(yù)測結(jié)果評價與分析

2.1 數(shù)據(jù)來源

本文所使用的血糖數(shù)據(jù)來自河南省人民醫(yī)院，是對病人連續(xù)監(jiān)測得到的，采用動態(tài)血糖監(jiān)測系統(tǒng) （CGMS）采集。CGMS每10秒中接收一次電流信號，每5分鐘將獲得的平均值轉(zhuǎn)換成葡萄糖值儲存起來，24小時共儲存288個葡萄糖值，共持續(xù)監(jiān)測72個小時。任意選取其中的24 個小時，其中前258個數(shù)據(jù)作為訓(xùn)練集，對血糖數(shù)據(jù)進行建模，采用最后30個樣本作為測試集。24小時所采取的某病人的血糖值如圖2所示。

2.2 ARIMA預(yù)測結(jié)果

將288個數(shù)據(jù)樣本輸入到MATLAB R2012A 中，采用ARIMA 的數(shù)據(jù)分析功能，得到血糖值自相關(guān)圖和偏自相關(guān)圖，如圖3所示。從圖3可知，血糖值序列有較高的自我相關(guān)性。需要進行差分處理，使血糖值序列變成平穩(wěn)的時間序列。取差分階數(shù)d＝1時，血糖序列基本平穩(wěn)，所以最佳差分階數(shù)為d＝1。其1階偏相關(guān)圖和自相關(guān)圖如圖4所示。

圖2 24小時采取的病人的血糖值

圖3 血糖值的ACF圖和PACF圖

圖4 1階偏相關(guān)和自相關(guān)圖

2.3 ARIMA模型參數(shù)的識別

采用從低階到高階逐步試探法來識別模型的參數(shù)，嘗試采用不同的p和q值，比較AIC和BIC取得最小值。

從表1的模型指標(biāo)對比可以看出，模型ARIMA （0，1，1）比較適合預(yù)測血糖值。采用一步預(yù)測法進行預(yù)測，預(yù)測結(jié)果如圖5所示：ARIMA 模型用于對血糖的預(yù)測，能夠好地把握血糖的變化趨勢，但精度不高。

表1 血糖序列的AIC和BIC檢測值

圖5 ARIMA 模型的預(yù)測結(jié)果

2.4 ARIMA－RBF組合預(yù)測

根據(jù)ARIMA （0，1，1）預(yù)測結(jié)果和血糖實際值得到血糖值的殘差序列，計算出ARIMA 預(yù)測誤差的趨勢和當(dāng)前血糖值的趨勢，結(jié)合當(dāng)前的血糖值對RBF進行建模和訓(xùn)練，得到預(yù)測值。組合模型的預(yù)測結(jié)果如圖6所示。

從圖6可知，ARIMA－RBFNN 組合預(yù)測模型不僅能夠準(zhǔn)確把握血糖的變化趨勢，而且能比較準(zhǔn)確地預(yù)測出短時間能的血糖值。

2.5 模型預(yù)測性能對比

為了能評價和比較預(yù)測結(jié)果，本文使用兩個性能指標(biāo)

（1）均方誤差 （mean squared error，MSE）

圖6 ARIMA－RBFNN 組合預(yù)測結(jié)果

（2）平均絕對相對誤差 （mean absolute error，MAE）

從表2的對比結(jié)果可知，ARIMA－RBFNN 的血糖預(yù)測精度要高于單一的ARIMA 模型的預(yù)測精度，預(yù)測誤差有較大的降低，是一種有效的血糖預(yù)測方法。用ARIMARBFNN 模型進行血糖值的預(yù)測，只需要利用CGMS提供的歷史數(shù)據(jù)，相對于結(jié)合各種生理模型的預(yù)測方法［5，6］，ARIMA－RBF神經(jīng)網(wǎng)絡(luò)組合預(yù)測算法簡單可靠。

表2 兩種預(yù)測方法誤差分析

2.6 針對外部因素對血糖的影響的奇異點發(fā)現(xiàn)和處理算法

通過預(yù)測結(jié)果的觀察，某些時刻預(yù)測值與真實值差距較大，經(jīng)過分析研究，差距較大的時間段大多有飲食、胰島素的注入、運動等事件的發(fā)生。如圖7 所示。預(yù)測時間段受飲食的影響，飲食后的預(yù)測精度大大降低。

定義滿足如下條件的點為奇異點

式中：threshold（k）——判定某點是否為奇異點的門檻，大小可調(diào)，當(dāng)k時刻真實血糖值y（k）與預(yù)測的血糖值（k｜k－PH）滿足上述條件時為奇異點。

通過實驗和觀察發(fā)現(xiàn)，奇異點一旦出現(xiàn)，將會持續(xù)3到5個采集點?？梢栽陬A(yù)測過程中不斷加入實際血糖值來判斷奇異點，為了提高預(yù)測的精度，提出了如下的算法：

（1）根據(jù)公式判斷采集點是否是奇異點；

（2）如果是奇異點，將當(dāng)前時刻未來4個時間段血糖預(yù)測更新為

圖7 受飲食影響的ARIMA－RBFNN 預(yù)測結(jié)果

（3）在發(fā)現(xiàn)奇異點時，其5個時間段后繼續(xù)執(zhí)行判斷奇異點，如果再次發(fā)現(xiàn)，則回到 （1）；如果沒有發(fā)現(xiàn)，則不做任何處理。

采用ARIMA－RBF預(yù)測模型對奇異點進行預(yù)測后的效果對比如圖8所示，進行預(yù)測后的MSE 和MAE 為0.0069和0.0616，而不是改進前的0.0609和0.1416。

圖8 改進后與改進前的對比

3 結(jié)束語

本文基于CGMS，提出一種融合ARIMA 和RBFNN 的組合預(yù)測模型，并與ARIMA進行對比。在組合預(yù)測的基礎(chǔ)上，針對外部因素對血糖的影響，提出一種發(fā)現(xiàn)和處理算法，能夠在血糖受外部因素影響時調(diào)整未來一段時間的血糖預(yù)測值。實例結(jié)果表明，該組合預(yù)測模型結(jié)合奇異點發(fā)現(xiàn)和處理算法能夠顯著提高預(yù)測能力和預(yù)測精度，很好地反應(yīng)血糖上升或下降的趨勢，具有很好的推廣和應(yīng)用價值。

［1］Ceriello A，Novials A，Ortega E，et al.Evidence that hyperglycemia after recovery from hypoglycemia worsens endothelial function and increases oxidative stress and inflammation in healthy control subjects and subjects with type 1diabetes［J］.Diabetes，2012，61 （11）：2993－2997.

［2］Georga E I，Protopappas V C，Polyzos D，et al.A predictive model of subcutaneous glucose concentration in type 1diabetes based on random forests ［C］//Engineering in Medicine and Biology Society.Annual International Conference of the IEEE，2012：2889－2892.

［3］Eren－Oruklu M，Cinar A，Quinn L，et al.Estimation of future glucose concentrations with subject－specific recursive linear models［J］.Diabetes Technology ＆ Therapeutics，2009，11（4）：243－253.

［4］Wang Z，Lai L，Xiong D，et al.Study on predicting method for acute hypotensive episodes based on wavelet transform and support vector machine［C］//3rd International Conference on Biomedical Engineering and Informatics.IEEE， 2010：1041－1045.

［5］Mo X，Wang Y，Wu X.Hypoglycemia prediction using extreme learning machine（ELM）and regularized ELM ［C］//Control and Decision Conference.IEEE，2013：4405－4409.

［6］Allam F，Nossair Z，Gomma H，et al.Prediction of subcutaneous glucose concentration for type－1diabetic patients using a feed forward neural network ［C］//International Conference on Computer Engineering ＆Systems.IEEE，2011：129－133.

［7］Pappada S M，Cameron B D，Rosman P M，et al.Neural network－based real－time prediction of glucose in patients with insulin－dependent diabetes［J］.Diabetes Technology ＆ Therapeutics，2011，13 （2）：135－141.

［8］Georga E I，Protopappas V C，Ardigo D，et al.Multivariate prediction of subcutaneous glucose concentration in type 1diabetes patients based on support vector regression ［J］.IEEE Journal of Biomedical and Health Informatics，2013，17 （1）：71－81.

［9］Facchinetti A，Sparacino G，Cobelli C.An Online SelfTunable Method to Denoise CGM Sensor Data ［J］.IEEE Trans Bio－Med Eng，2010，57 （3）：634－641.

［10］Ye R，Suganthan P N，Srikanth N，et al.A hybrid ARIMADENFIS method for wind speed forecasting ［C］//IEEE International Conference on Fuzzy Systems.IEEE，2013：1－6.

［11］Shumway R H，Stoffer D S.Time series analysis and its applications：With R examples［M］.Springer，2010.

［12］Box G E P，Jenkins G M，Reinsel G C.Time series analysis：Forecasting and control［M］.John Wiley ＆Sons，2013.