王延年，雍永強(qiáng)，賈曉燦，李全忠，孫長青＋

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

0 引 言

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

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

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

1 方 法

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

1.1 組合預(yù)測(cè)模型

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

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

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

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

（2）利用 （1）得到結(jié)果得到ARIMA 產(chǎn)生的預(yù)測(cè)殘差：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 時(shí)間段后的預(yù)測(cè)結(jié)果（t＋PH ｜t）。

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

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

1.3 ARIMA模型

ARIMA 的一般表達(dá)式為

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

ARIMA 處理過程如下：

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

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

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

（2）模型識(shí)別。p，d，q 的選擇是ARIMA 建模的關(guān)鍵。首先對(duì)血糖數(shù)據(jù)進(jìn)行差分處理，每次差分比較其自相關(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］，計(jì)算AIC和BIC的方程如下［11］

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

2 預(yù)測(cè)結(jié)果評(píng)價(jià)與分析

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

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

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

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

圖2 24小時(shí)采取的病人的血糖值

圖3 血糖值的ACF圖和PACF圖

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

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

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

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

表1 血糖序列的AIC和BIC檢測(cè)值

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

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

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

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

2.5 模型預(yù)測(cè)性能對(duì)比

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

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

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

（2）平均絕對(duì)相對(duì)誤差 （mean absolute error，MAE）

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

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

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

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

定義滿足如下條件的點(diǎn)為奇異點(diǎn)

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

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

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

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

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

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

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

圖8 改進(jìn)后與改進(jìn)前的對(duì)比

3 結(jié)束語

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

［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.