基于粒子群優(yōu)化支持向量機(jī)的癲癇發(fā)作預(yù)測*

劉金森 黃煒嘉 李效龍

（江蘇科技大學(xué)電子信息學(xué)院 鎮(zhèn)江 212100）

1 引言

癲癇是一種危害極大的神經(jīng)系統(tǒng)疾病，具有突發(fā)性與反復(fù)性，嚴(yán)重影響了病人的日常生活。越來越多的研究表明，癲癇發(fā)作具有可預(yù)測性［1］。如果能夠在癲癇發(fā)作之前發(fā)出警報，提醒病人家屬或者醫(yī)生對病人采取相應(yīng)措施來抑制癲癇發(fā)作，就能夠減緩癲癇病人的痛苦，從而降低癲癇發(fā)作給患者帶來的危害。目前，癲癇發(fā)作自動預(yù)測已經(jīng)成為腦電疾病研究領(lǐng)域的熱點(diǎn)問題。

癲癇發(fā)作分為發(fā)作前期、發(fā)作期、發(fā)作后期和發(fā)作間期四個階段［2］。癲癇預(yù)測的主要目標(biāo)是對癲癇發(fā)作前期與癲癇發(fā)作間期的腦電信號進(jìn)行區(qū)分，識別出發(fā)作前期信號并發(fā)出預(yù)警，進(jìn)而達(dá)到癲癇預(yù)測目的。

對于發(fā)作前期與發(fā)作間期的二分類問題主要分為特征提取和模式識別兩個部分。在癲癇病人的腦電信號特征提取方面，學(xué)者們做了大量的研究工作，主要分為以下四類：1）時域方法：均值、標(biāo)準(zhǔn)差［3］、Hjorth 復(fù)雜度［4］等；2）頻域方法：功率譜［5］；3）時頻分析方法：短時傅里葉變換［6］、小波變換［7］、小波包變換［8］、經(jīng)驗(yàn)?zāi)B(tài)分解［9］；4）非線性動力學(xué)方法：最大Lyapunov 指數(shù)［10］、熵［11］。時域分析方法主要是分析信號的幾何特性，如對癲癇棘波、尖波等特征波形進(jìn)行分析；頻域分析主要是通過傅里葉變換、Weich 方法、AR 模型參數(shù)方法獲得信號的功率譜，將腦電信號各頻段的功率、相對功率、絕對功率等作為特征進(jìn)行研究。然而，時域和頻域方法都只針對信號單一角度進(jìn)行度量，且通常需要假設(shè)腦電信號為平穩(wěn)信號，這與腦電信號的非平穩(wěn)性相矛盾，小波包變換具有良好的正交性和局部特性，可以對信號進(jìn)行多尺度、多分辨分析，非常適合于腦電信號等非平穩(wěn)信號的分析研究。

隨著機(jī)器學(xué)習(xí)的發(fā)展，涌現(xiàn)出越來越多的分類方法，應(yīng)用在腦電信號分析領(lǐng)域的有K-近鄰、支持向量機(jī)、線性判別式以及人工神經(jīng)網(wǎng)絡(luò)等［12］。其中，支持向量機(jī)是應(yīng)用最多的分類方法，它不涉及人工神經(jīng)網(wǎng)絡(luò)中的網(wǎng)絡(luò)結(jié)構(gòu)選擇，并且避免了“維數(shù)災(zāi)難”［13］。目前癲癇預(yù)測方面涉及到支持向量機(jī)作為分類方法，參數(shù)設(shè)置多采用人工經(jīng)驗(yàn)選取和網(wǎng)格搜索方法［14］，憑經(jīng)驗(yàn)選取具有一定的隨機(jī)性，網(wǎng)格搜索在較大范圍內(nèi)進(jìn)行搜索，存在耗時較長等問題。因此，本研究采用小波包變換和支持向量機(jī)結(jié)合來構(gòu)建癲癇發(fā)作預(yù)測模型，同時為提高分類準(zhǔn)確率、減少計算時間，采用粒子群算法對支持向量機(jī)的參數(shù)進(jìn)行優(yōu)化。

2 方法和原理

2.1 小波包系數(shù)能量

小波包變換是軸承故障特征提取中常用的信號分析方法，近年來也越來越多地應(yīng)用于腦電信號的處理當(dāng)中［15～16］。癲癇發(fā)作前期與發(fā)作間期的腦電信號不同頻段具有不同的變化趨勢，應(yīng)用小波包變換的自適應(yīng)性可以很好地提取出信號不同頻段特征。

小波包變換滿足如下的二尺度方程：

其中：hk和gk為一組正交濾波器，滿足gk=(-1)nh1-k。

Mallet小波包快速分解算法描述如下：

信號的時域能量為

由帕斯瓦爾定理［16］可知：信號在時域能量上的能量與小波包系數(shù)能量相等，因此，信號各頻段的能量可表示為

其中：j代表小波包分解層數(shù)，k表示信號頻率成分，N 代表信號采樣點(diǎn)數(shù)，Cj,k代表小波包分解系數(shù)。

為了消除量級對數(shù)據(jù)的影響，本文采用相對能量作為特征，具體表示如下：

其中：Ej,k代表信號各頻段能量值，代表信號各頻段總能量值。

2.2 支持向量機(jī)

支持向量機(jī)主要思想如下：首先通過一個非線性映射Φ(x)將y類腦電信號特征x映射到高維空間，然后找到一個超平面使得所有的點(diǎn)到超平面的距離大于一定值［17］。算法具體步驟如下。

1）超平面的構(gòu)造歸結(jié)為凸二次規(guī)劃問題：

其中：w為超平面法向量，C>0 為懲罰參數(shù)，ξi≥0為松弛變量，b為超平面的偏置向量。

2）利用拉格朗日函數(shù)求目標(biāo)函數(shù)極值，設(shè)拉格朗日函數(shù)為

其中：0 ≤αi≤C為Lagrange系數(shù)。

3）對式（7）求解最優(yōu)化問題，得到最優(yōu)解：

2.3 粒子群算法優(yōu)化支持向量機(jī)

粒子群算法是基于對鳥類捕食行為進(jìn)行建模衍生出來的一種智能優(yōu)化算法。該算法把群體中的個體視為空間中的點(diǎn)，具有一定的初始速度，在此過程中依據(jù)自身和同伴的飛行經(jīng)驗(yàn)對自身的飛行速度、前進(jìn)方向進(jìn)行動態(tài)調(diào)整［18～19］。

粒子的速度和位置更新公式如下：

式中：xi=(x1,x2,…,xn)表示粒子在N 維空間的位置，飛行速度表示為矢量vi=(v1,v2,…,vn)表示粒子在N 為空間位置，pbest表示粒子i的個體最優(yōu)位置，Gbest表示粒子群的全局最優(yōu)位置；c1，c2為常數(shù)，表示學(xué)習(xí)因子；γ1,γ2∈( )0,1 ；w表示慣性因子。

粒子群算法優(yōu)化支持向量機(jī)的流程圖如圖1所示。

圖1 PSO優(yōu)化SVM算法流程圖

其算法具體流程如下：

1）初始化參數(shù)：種群規(guī)模N、最大迭代次數(shù)i、慣性因子w、參數(shù)c、σ的變化范圍；

2）根據(jù)種群規(guī)模隨機(jī)產(chǎn)生一組粒子(c,σ)；

3）適應(yīng)度函數(shù)為十折交叉驗(yàn)證計算出的平均準(zhǔn)確率，根據(jù)適應(yīng)度函數(shù)分別計算每個粒子的適應(yīng)度函數(shù)值；

4）每個粒子的適應(yīng)度函數(shù)值即為當(dāng)前個體最優(yōu)值，比較所有粒子的個體最優(yōu)值，值最大的作為全局最優(yōu)值，對應(yīng)粒子位置(c,σ)作為全局最優(yōu)位置；

5）根據(jù)式（12）和式（13）更新粒子速度和位置；

6）根據(jù)更新后的粒子位置重新計算粒子適應(yīng)度函數(shù)值；

7）比較更新后的適應(yīng)度值與個體極值、群體極值，若優(yōu)于更新前，則用更新后的值進(jìn)行替換；

8）判斷迭代次數(shù)是否小于Tmax，若達(dá)到，則輸出當(dāng)前全局最優(yōu)解(cbest,σbest)，否則返回步驟5）。

本文參數(shù)設(shè)置如下：種群規(guī)模N=20，最大迭代次數(shù)Tmax=200，懲罰因子c最小值cmin=0.1，最大值cmax=100，核函數(shù)σ的最小值σmin=0.01，最大值σmax=1000 ，慣性因子w=1，c1=1.5，c2=1.7。

2.4 輸出結(jié)果后處理

為了減少錯誤警報，還需要對SVM 輸出結(jié)果進(jìn)行判決分析。本文提出一種基于閾值的判斷方法，具體過程如下：

1）統(tǒng)計每隔1min 的時間窗，支持向量機(jī)的輸出為發(fā)作前期的標(biāo)簽數(shù)目：

其中：numpreictal[t]表示第t 個時間窗內(nèi)的發(fā)作前期樣本個數(shù)；labeli表示第i個樣本的SVM輸出標(biāo)簽。

2）使用numpreictal[t]的值可以判斷在過去的這1min 內(nèi)信號處于發(fā)作前期還是發(fā)作間期，具體如下：

其中：P值為一個閾值，P∈[1,6]，6表示以10s為窗口大小，1min 內(nèi)的窗口數(shù)。當(dāng)numpreictal[t]≥p時，state[t]=1，即認(rèn)為在過去的1min內(nèi)，信號處于發(fā)作前期狀態(tài)；否則處于發(fā)作間期。

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

3.1 實(shí)驗(yàn)數(shù)據(jù)與處理

本文采用CHB-MIT 癲癇開源數(shù)據(jù)庫（https：//www.physionet.org/content/chbmit/1.0.0/），該數(shù)據(jù)庫記錄了23 位癲癇病人的數(shù)據(jù)，所有信號以16 比特分辨率進(jìn)行采樣，采樣頻率為256Hz。

根據(jù)本次實(shí)驗(yàn)數(shù)據(jù)集的特點(diǎn)以及相應(yīng)實(shí)驗(yàn)驗(yàn)證，本文將每次癲癇發(fā)作前10min 定為發(fā)作前期。同時為了實(shí)現(xiàn)更快地對數(shù)據(jù)進(jìn)行處理，減少計算時間，本文采用滑動窗口技術(shù)對信號進(jìn)行分段處理。發(fā)作前期、發(fā)作間期樣本均采用10s 無重疊窗口，窗口步長為10s。

3.2 特征提取結(jié)果與分析

已有實(shí)驗(yàn)表明Daubechies 小波家族對腦電數(shù)據(jù)具有很好的處理效果，經(jīng)試驗(yàn)采用db5 小波對信號分析，五層分解后共得到25=32 個頻率成分，各節(jié)點(diǎn)及其頻率成分如表1所示。

表1 小波包變換分解結(jié)果

腦電的分析范圍多為0 ～32Hz，所以對分解后得到前8個頻段進(jìn)行分析，其相應(yīng)頻譜如圖2所示。

圖2 第五層前8個節(jié)點(diǎn)頻譜圖

對經(jīng)過小波包變換分解后各頻段分解系數(shù)求取相對能量值來反映病人在發(fā)作前期、發(fā)作間期的信號變化情況，以第一位病人的第一次發(fā)作為例，圖3、圖4分別為該病人第一次發(fā)作前期一個10s窗口的小波包系數(shù)相對能量情況。

圖3 發(fā)作前期小波包系數(shù)相對能量

圖4 發(fā)作間期小波包系數(shù)相對能量

通過對比兩圖，可以發(fā)現(xiàn)發(fā)作前期與發(fā)作間期的小波包系數(shù)能量占比有較為明顯的差異，癲癇發(fā)作前期各頻段能量由大到小依次為S5_0、S5_1、S5_6、S5_7、S5_5、S5_4、S5_3、S5_2，癲癇發(fā)作間期各頻段能量由大到小依次為S5_0、S5_2、S5_1、S5_3、S5_4、S5_5、S5_6、S5_7。結(jié)果表明，小波包系數(shù)能量比值能夠很好地反映發(fā)作前期與發(fā)作間期腦電信號不同頻率的變化情況，由此證實(shí)了小波包系數(shù)相對能量作為特征區(qū)分兩類信號的有效性。

3.3 識別結(jié)果與分析

本文選取5 位病人前兩次發(fā)作前期10min、發(fā)作間期10min 的腦電數(shù)據(jù)，對樣本數(shù)據(jù)進(jìn)行特征提取。將其構(gòu)成240×8 維的特征向量作為訓(xùn)練數(shù)據(jù)。選擇第三次癲癇發(fā)作提取出來的特征向量作為測試數(shù)據(jù)進(jìn)行分析。得到粒子群優(yōu)化支持向量機(jī)的適應(yīng)度曲線如圖5所示。

圖5 粒子群優(yōu)化算法適應(yīng)度曲線

從圖5 中可以看出經(jīng)過200 次迭代處理，尋找得到的最佳c 和g 組合為（1.2399，7.5833），訓(xùn)練樣本交叉驗(yàn)證準(zhǔn)確率為98.3051%。

為了驗(yàn)證粒子群算法優(yōu)化的有效性，本文設(shè)計了支持向量機(jī)隨機(jī)參數(shù)選擇和網(wǎng)格搜索兩種方法進(jìn)行對比實(shí)驗(yàn)。其中網(wǎng)格搜索范圍設(shè)置為c,σ∈( 2-8,28)，步長大小設(shè)置為0.8。圖6 為網(wǎng)格搜索優(yōu)化支持向量機(jī)參數(shù)選擇等高線圖。

圖6 網(wǎng)格搜索等高線圖

從圖6 中可以看出，該方法尋找到的最優(yōu)參數(shù)組合為（11.3137，2），對于測試集樣本的準(zhǔn)確率為90.8333%。三種方法的準(zhǔn)確率、所需時間以及尋找到的最佳參數(shù)組合具體如表2所示。

表2 三種方法性能比較

從表2 中數(shù)據(jù)可以看出，隨機(jī)參數(shù)選方法雖然耗時較短，但分類效果最差，分類準(zhǔn)確率只有63.3%；網(wǎng)格搜索方法的分類準(zhǔn)確率與粒子群算法相近，但是耗時卻約為粒子群優(yōu)化算法的兩倍，這對于計算速度要求較高的癲癇發(fā)作預(yù)測系統(tǒng)來說，顯然不太適合。而粒子群優(yōu)化算法在處理速度和識別準(zhǔn)確率方面均達(dá)到較好效果，所以驗(yàn)證了本文方法的有效性。

3.4 癲癇發(fā)作預(yù)測結(jié)果與分析

圖7 為對第1 位病人的第三次發(fā)作預(yù)測結(jié)果，時間選取為距離癲癇發(fā)作前的10min。其中圖（a）、圖（b）、圖（c）、圖（d）分別為發(fā)作前期原始信號時域圖、支持向量機(jī)輸出結(jié)果、閾值決策結(jié)果、發(fā)作前期處理后時域圖。橫軸為時間，間隔為10s，即為本文采用的滑動窗口大小。

從圖7（b）中可以看出，支持向量機(jī)在第10s 輸出為0。如果不采用閾值決策方法對輸出結(jié)果進(jìn)行處理，該段信號就會被判定為發(fā)作間期，從而不會發(fā)出預(yù)警，這增加了癲癇預(yù)測系統(tǒng)的誤報率。若采用本文的后處理方法，如圖（c）所示對每隔1min的6 個輸出結(jié)果進(jìn)行閾值決策，閾值P 取為4。0～1min 支持向量機(jī)輸出3 個1，小于閾值4，經(jīng)決策后判決為0；2min～3min 支持向量機(jī)輸出4 個1，等于閾值4，經(jīng)決策后判決為1，即系統(tǒng)認(rèn)定為發(fā)作前期，在3min這一時刻發(fā)出預(yù)警。

表3 為5 名癲癇病人的實(shí)驗(yàn)預(yù)測結(jié)果，平均準(zhǔn)確率為83.67%，平均敏感性為75.34%，平均特異性為92%，平均誤報率為0，平均預(yù)測時間為6.8min。

表3 五位癲癇病人預(yù)測結(jié)果

4 結(jié)語

癲癇發(fā)作前期與發(fā)作間期的腦電信號不同頻段能量具有明顯不同，采用小波包變換方法可以很好地提取出這一特征。同時經(jīng)過粒子群算法優(yōu)化后的支持向量機(jī)與隨機(jī)參數(shù)方法、網(wǎng)格尋優(yōu)方法相比，識別準(zhǔn)確率更高、運(yùn)行時間較短。此外，經(jīng)實(shí)驗(yàn)驗(yàn)證了本文提出的閾值后處理方法在減少系統(tǒng)誤報率方面的有效性。實(shí)驗(yàn)結(jié)果表明，本文方法準(zhǔn)確率高，并且無誤報率，預(yù)測時間也在合理范圍內(nèi)，為進(jìn)一步研究實(shí)時性癲癇發(fā)作系統(tǒng)提供了參考價值。