基于彈性變分模態分解的癲癇腦電信號分類方法*

景鵬 張學軍2)? 孫知信

1) (南京郵電大學電子與光學工程學院,微電子學院,南京 210023)

2) (南京郵電大學,射頻集成與微組裝技術國家地方聯合工程實驗室,南京 210023)

3) (南京郵電大學,江蘇省郵政大數據技術與應用工程研究中心,南京 210003)

4) (南京郵電大學,國家郵政局郵政行業技術研發中心(物聯網技術),南京 210003))

癲癇腦電信號分類對于癲癇診治具有重要意義.為了實現病灶性與非病灶性癲癇腦電信號的分類,本文利用彈性網回歸重構變分模態分解算法,提出彈性變分模態分解算法并將其應用到所提癲癇腦電信號分類方法中.該方法先將原信號分割成多個子信號,并對各子信號進行彈性變分模態分解,然后從分解后的不同變分模態函數中提取精細復合多尺度散布熵作為特征,最后利用支持向量機進行分類.針對癲癇腦電的公共數據集,最終的實驗結果表明,準確率、靈敏度和特異度三個性能指標分別達到92.54%,93.22%和91.86%.

1 引 言

癲癇是一種由于人腦神經系統紊亂造成腦功能障礙的疾病,據世界衛生組織2019年統計,全世界有超5000萬人患有癲癇[1],對具有抗藥性的癲癇患者需要進行侵入性手術治療,切除大腦中的部分致癇區域,因此手術前準確識別大腦中的致癇區域至關重要.腦電圖 (electroencephalograph,EEG)記錄了人腦皮層發出的輕微電流信號,癲癇發作時人腦中樞神經系統會出現異常使得同步神經元突然放電,相應的腦電圖就會呈現出異常波動,從人腦致癇區域獲得的腦電信號稱為病灶性腦電信號,從非致癇區域獲得的腦電信號稱為非病灶性腦電信號[2],因此對病灶性腦電信號和非病灶性腦電信號進行分類,有助致癇區域的識別.

目前國內外對于癲癇腦電信號的分類技術主要由特征提取與特征分類組成,其本質可歸結為模式識別問題[3].特征提取先要選擇合適的特征,如高階矩[4]、概率密度函數[5]、自回歸模型系數[6]、各類熵[7-10]等均已成功應用到癲癇腦電信號的分類中.文獻[11]提出了一種新的復雜性度量特征,即精細復合多尺度散布熵(refined composite multiscale dispersion entropy,RCMDE),其作為特征在處理非線性不平穩生物信號時效果較好.此外,特征提取通常在時域、頻域或時頻域內進行,其中又以時頻域分析法應用較多,如短時傅里葉變換[12]、小波變換[13]、經驗模態分解等,經驗模態分解作為一種自適應時頻處理方法應用較廣,但存在端點效應、模態混疊等缺陷[14].由文獻[14]提出的變分模態分解 (variational mode decomposition,VMD)是一種新的基于完全非遞歸思想的時頻域信號分析方法,其主要思想是將信號分解成若干個圍繞在各個中心頻率附近的窄帶變分模態分量且中心頻率不斷變化,與經驗模態分解通過不斷循環篩分獲取本征模態函數不同,VMD通過尋找約束變分模型的最優解自適應獲取變分模態函數(variational mode function,VMF).VMD 在構建初始是針對盲信號,使用Tikhonov正則化思想去構建一個約束變分方程[14],但本質是通過對最小二乘回歸加殘差平方和二次懲罰項,從而達到收縮估計系數的目的,這也可以理解為使用了線性規劃中的嶺回歸(ridge regression)來構建約束方程,雖然嶺回歸與最小二乘回歸相比降低了估計值的方差,提高了估計精度,但是它的回歸結果中包含所有的預測變量,沒有進行變量選擇,因此會影響模型的準確性[15].而套索回歸(LASSO regression)使用一次懲罰項來對變量進行收縮,相較于嶺回歸收縮程度要小,能選出更精確的模型,但是如果預測變量具有群組效應或者預測變量相關性較強時,則會導致估計不穩定[15].Zou和Hastie[16]提出的彈性網回歸(elastic net regression)方法綜合了嶺回歸和套索回歸的思想,兼有套索回歸和嶺回歸的優點.

本文利用彈性網回歸對變分模態分解算法進行改進,構建一種新的時頻域信號分析方法,稱為彈性變分模態分解 (elastic variational mode decomposition,EVMD),并通過在 EVMD 域中提取精細復合多尺度散布熵作為癲癇腦電信號的特征,然后利用支持向量機 (support vector machine,SVM)實現病灶性腦電信號和非病灶性腦電信號的分類.此外,將本文所提方法應用到實際癲癇腦電數據集中進行仿真分析,以證明該方法的有效性.

2 方 法

應用本文所提方法處理癲癇腦電信號的步驟如圖1所示.先利用彈性變分模態分解算法對分割后的癲癇腦電子信號進行分解,從分解后的各VMF中提取RCMDE作為特征,最后利用支持向量機對特征進行分類,進而完成對病灶性腦電信號與非病灶性腦電信號的區分.本節將對所提方法中使用到的主要算法進行敘述.

2.1 彈性變分模態分解

圖1 彈性變分模態分解處理癲癇腦電信號的流程圖Fig.1.Architecture of processing epileptic EEG by EVMD.

本文提出的彈性變分模態分解是一種完全非遞歸式的時頻域信號分析與處理方法,與原始變分模態分解算法不同的是,EVMD利用彈性網回歸代替VMD中構建約束變分方程時使用的嶺回歸,即EVMD算法構建的約束變分方程中既存在L2正則化項,也存在L1正則化項.本節將對EVMD算法的構造與求解進行詳細敘述.

首先,通過以下步驟建立一個約束變分模型：1)通過希爾伯特變換求出每個模態函數的解析信號,進而可以獲得信號的單邊頻譜； 2)通過一個混合指數項將每個模態函數調制到對應中心頻率的基頻帶上； 3)通過彈性網回歸方法施加L1正則化與L2正則化來估計帶寬,即求信號梯度二范數的平方與梯度一范數的和.建立的約束變分模型為

其中,x為原腦電信號,{uk}={u1,u2,···,uK}為所有模態的集合,{ωk}={ω1,ω2,···,ωK}為對應中心頻率的集合,為沖激函數,?t為對t求偏導.

接著對約束變分模型求最優解.先將拉格朗日乘法算子以及損失項加入到約束變分模型中,從而將約束變分模型轉化為一個非約束變分模型,得到的增廣拉格朗日函數為

其中λ為拉格朗日乘法算子,α為二次懲罰因子,β為一次懲罰因子,為損失項.

然后利用交替方向乘子法來求解(2)式的最優解,迭代公式為：

迭代停止條件為：

其中,η為拉格朗日乘子更新參數,ε為收斂容限.

下面介紹模態uk與中心頻率ωk的詳細更新求解步驟,即迭代式(3)式與(4)式的具體推導過程.

A模態uk更新

步驟1將(3)式改寫為如下等值公式：

步驟2將(7)式進行整體傅里葉變換：

步驟 3將 (8)式中ω+ωk替換為ω,即將帶有懲罰因子的項中的ω用ω-ωk代替,得：

步驟4因為實信號具有埃米爾特對稱性,因此將(9)式改為非負頻域上的半空間積分：

步驟5將(10)式看為一個二次優化問題進行求解,解得：

B中心頻率ωk更新

步驟1因為中心頻率ωk不出現在重構保真項中,因此可將(4)式重新寫為

步驟2與求解模態uk的操作類似,將(12)式進行傅里葉變換到頻域,并最終變換成非負頻域上的半空間積分：

步驟3求解(13)式二次優化問題,可得：

至此,完整的EVMD算法已經介紹完畢,其具體流程可以參考算法1所示.

Algorithm 1：EVMD

2.2 精細復合多尺度散布熵

精細復合多尺度散布熵是文獻[11]在散布熵的基礎上對信號進行多尺度量化得到的,避免了散布熵因單一尺度上處理信號會出現復雜性特征提取不完全等問題.文獻[11]詳細介紹了RCMDE的原理,這里僅簡述其對腦電信號處理時的計算步驟：

步驟1假設腦電信號x經彈性變分模態分解后得到的信號uk是長度為Ψ的時間序列,則uk的第h個粗粒度近似信號為

其中,γ={1,2,···,N},N=Ψ/τ,τ為尺度因子.

步驟2使用正態累積分布函數將粗粒度近似信號a映射到b={b1,b2,···,bN}上：

其中,μ表示均值,σ表示粗粒度近似a的標準差.b的范圍從0到1.

步驟3為了實現信號的多尺度化,將b以線性變換的形式映射到{1,2,··,c}中,記為z,即

其中R(*)是取整運算,c代表類別個數.

步驟4若嵌入維數標記為m和時間延遲標記為d,則時間序列定義為

其中i={1,2,···,N-(m-1)d}.

步驟5設每一個時間序列對應一個散布模 式Θv0v1···vm-1,其 中v= 1,2,···,c.若=v0,對應的散布模式為Θv0v1···vm-1.

步驟6計算每個散布模式Θv0v1···vm-1的概率

步驟7根據香農熵的定義,信號uk的精細復合多尺度散布熵為

3 實驗結果與分析

3.1 實驗數據集

本文使用的癲癇腦電信號數據集是來自瑞士伯爾尼大學神經科提供的一個公共數據集[17].該數據集數據采自5位患有顳葉病灶癲癇的病人,共包括3750對病灶性數據段和3750對非病灶性數據段,每對數據段均包含兩列數據,分別采自同一個區域相鄰的兩個通道.數據段的采樣頻率為512 Hz,時間為 20 s[18].本文選用該數據集中全部腦電數據段作為實驗數據,即選用3750對病灶性腦電數據段與3750對非病灶性腦電數據段,每對數據段均包含兩列數據.

3.2 實驗

3.2.1 EEG 數據段分割

先將原始腦電數據經過0.5—100 Hz的帶通濾波器進行預處理,從而去除部分偽跡與噪聲.然后將每對時長為20 s的數據段分割為時長均為1 s的子數據段,且相鄰兩個子數據段之間的時間重疊率為50%,因此每個數據段共分割為39個子數據段.之所以需要對原始腦電數據進行分割,第一是因為EVMD算法與VMD一樣在處理長時間非平穩非線性的腦電信號時模態的譜帶會隨時間劇烈變化影響分解效果,而將信號分割為短時間的子信號,可以一定程度規避這種局限性.第二是因為可以減小腦電數據中隨機波動點的影響,比如假設某數據段中由于噪聲等原因出現數個非預期的隨機波動點,若對整個數據段進行分析顯然這些點會影響最終分析結果,而對若干子數據段進行分析,這些點只會影響部分子數據段的分析結果,而對最終所求結果的影響將降低.

3.2.2 EVMD 分解

因為測試時間和被試的不同,各頻帶段的腦電信號的中心頻率會隨著神經活動而發生輕微的變化,簡單地固定中心頻率的帶通濾波無法消除這種變化的影響,而EVMD旨在將信號分解為圍繞在中心頻率附近的變分模態分量,即變分模態函數VMF,且其中心頻率是不斷迭代變化的,故可以捕捉到這種變化.

根據2.1節對EVMD的介紹,可以看出EVMD算法在對信號進行分解之前需要分別設定： 分解模態的個數K,二次懲罰因子α,一次懲罰因子β三個參數.對于懲罰因子α,β,其值過大會引起模態重疊,較小會引入噪聲,本文根據經驗建議將α設定為 2000,β設定為20.對于分解模態的個數K,若K值過大會造成過分解,產生無用分量,同時增加計算復雜度,若K值過小會造成欠分解,使部分帶限信號分解不出來造成原信號信息的丟失.為尋求合適的K值,本文先隨機選擇一對病灶性腦電數據段,并分別進行3重、4重、5重變分模態分解,即將K值分別設為 3,4,5.圖2所示為不同K值下,各個變分模態函數分量其中心頻率ω隨迭代次數的變化曲線.可以看出,當K= 3 時,變分模態函數VMF1—VMF3的中心頻率均在40 Hz以下.當K= 4 時,變分模態函數 VMF1—VMF4 的中心頻率仍在 40 Hz以下.但當K= 5 時,可以看到VMF5的中心頻率超過了50 Hz.一方面由于腦電活動的信息大都包含在低頻帶(頻率 ＜ 40 Hz),另一方面為了防止信號的過分解或欠分解,本文選取K= 4.

圖2 不同K值下中心頻率隨迭代次數的變化曲線Fig.2.The curves of the center frequency with the number of iterations under different K values.

3.2.3 RCMDE 提取特征

根據(23)式,計算RCMDE需要分別設置：嵌入維度m,類別個數c,時間延遲d以及尺度因子τ四個參數.類別個數c若過大會導致具有較大差異的兩個量被歸為同一類,過小則導致具有較小差異的兩個量被歸為不同類,根據文獻[11]的建議,將c值設為6.對于嵌入維度m,為了保證統計可靠性,文獻 [12]中建議cm＜Ψ,其中Ψ為數據長度,本文將采樣率 512 Hz,采樣時間 20 s的數據段分割為時長為 1 s的子數據段,因此Ψ= 512,而 64= 1296 ＞Ψ,故將m值設為 3.時間延遲d為正整數,但d＞ 1 會導致模態混疊,故取 1.尺度因子τ決定信號粗粒化程度,為了選取合適的τ值,先將τ最大值設為15.

選用數據集中全部3750對病灶性腦電信號和3750對非病灶性腦電信號,對分割后的子數據段進行4重EVMD分解,從分解后得到的4個VMF(VMF1—VMF4)中分別提取15個尺度因子下的RCMDE,其特征熵值的均值加減標準差隨尺度因子變化曲線如圖3所示.圖3中不同尺度因子下的RCMDE均值用曲線相連,每個均值點上下兩個點即表示加減標準差.從圖3可以看出,從非病灶性腦電信號提取的RCMDE 特征熵值基本比病灶性腦電信號提取的RCMDE熵值大,這說明人腦非致癇區域的神經活動相對致癇區域更活躍,也說明非致癇區域提取的腦電信號相較致癇區域提取的腦電信號更隨機、更不平穩,這一結果也與文獻[18]的研究結果一致.此外,之所以會出現從不同VMF中提取的兩類信號RCMDE特征熵值的均值差異性不同的現象,是因為不同VMF是體現不同中心頻率上的腦電信號數據,而人腦非致癇區域與致癇區域在不同頻帶段上的神經活動會不同,比如經4重EVMD分解下VMF3代表的是中心頻率為20 Hz附近頻帶段的信號,此時兩類腦電信號RCMDE熵值的均值的差值與其他VMF相比較小,表明在20 Hz附近人腦非致癇區域與致癇區域的神經活動更相近,產生的腦電信號的差異性較其他頻帶段較小.但是隨著尺度因子的增大,總體非病灶性腦電信號與病灶性腦電信號提取的RCMDE特征熵值的均值的差值會有一個先增大再減小的過程,尤其是VMF2尺度因子等于15時,VMF3尺度因子大于12時,會出現病灶性腦電信號提取的RCMDE特征熵值的均值更大的反?，F象,這是因為信號被過度粗粒化從而出現失真,說明尺度因子不能過大.本文從計算量的角度綜合考慮將尺度因子τ設為7,即每個數據段提取7個尺度上的RCMDE值.

圖3 從各 VMF 中提取的 RCMDE 熵值的均值 (± 標準差)隨尺度因子變化曲線Fig.3.The curve of mean value(± SD) of RCMDE computed from VMF.

為了更好評估從非病灶性腦電信號提取的RCMDE特征與從病灶性腦電信號提取的RCMDE特征的區分度,這里使用學生檢驗的p值來評估其統計學上的差異性,其結果如表1所列.從表2可以看出數據經4重EVMD分解后得到的4個VMF中提取的RCMDE特征p值均小于0.05,因此在統計學上具有顯著差異,說明RCMDE可以作為區別病灶性和非病灶性癲癇腦電數據的分類特征.

表1 從各 VMF 中提取的 RCMDE 特征 p 值Table 1.The p values of RCMDE computed from VMF.

表2 EVMD 與 VMD 實驗結果對比Table 2.Comparison of experimental result between EVMD and VMD.

3.3 分類結果

選用數據集中全部3750對病灶性癲癇腦電數據段和3750對非病灶性癲癇腦電數據段,按上述方法進行處理,并將得到的每段腦電數據的RCMDE特征送入SVM進行特征分類,SVM選用線性核函數,通過網格搜索法將懲罰參數設定為0.52,并重復進行10次5折交叉驗證實驗,采用準確度、靈敏度和特異度這三個指標對最終分類結果進行度量,結果如圖4所示.由圖4可知10次5折交叉驗證實驗的平均準確度、靈敏度和特異度分別可達 92.54%,93.22%,91.86%.

為了比較本文提出的EVMD算法與原始VMD在處理病灶性癲癇腦電信號與非病灶性信號分類中的性能,將本文實驗中的EVMD換成原始VMD,其余步驟不變,重復進行實驗,對比實驗的結果如表2所列.可以看出,EVMD算法相較原始VMD算法,三個性能指標均有3個百分點以上的提高,因此按本文所提方法進行病灶性癲癇腦電信號與非病灶性癲癇腦電信號的分類實驗中,本文所提EVMD算法相較原始VMD算法性能具有一定優勢.

圖4 10 次 5 折交叉驗證實驗結果折線圖Fig.4.The line chart of the results by 5-fold cross validation for 10 times.

對相同的實際癲癇腦電信號數據集,將本文方法得到的結果與其他文獻的結果相比較,結果如表3所列,可以看出本文所提方法性能優于其他方法.

表3 本文方法與其他方法對比Table 3.Comparison between proposed method and previously published methods.

4 結 論

本文利用彈性網回歸重構了變分模態分解的約束方程,提出一種新的時頻域信號分析與處理算法,稱為彈性變分模態分解,并給出了具體推導過程.提出一種基于彈性變分模態分解,從分解后的VMF中提取RCMDE特征,并利用支持向量機對其進行分類的病灶性癲癇腦電信號與非病灶性癲癇腦電信號的分類方法.應用本文提出的方法對公共癲癇腦電數據集進行處理,得到最終準確率為92.54%,靈敏度 93.22%,特異度 91.86%,并與針對同一公共數據集的其他處理方法進行比較,證明了該方法的有效性.然而,彈性變分模態分解的參數(分解模態的個數K、懲罰因子α與β)只能通過經驗設定或者根據處理信號不斷試錯設定,無法自動獲取最優參數,因為彈性變分模態分解與原始變分模態分解結構類似,所以現在針對變分模態分解參數的優化方法,如相關性分析優化[21]、粒子群算法優化[22,23]等方法,理論上應該也可以應用到本文提出的彈性變分模態算法上,因此,對彈性變分模態分解的參數優化將是后續研究的一個可能方向.