景鵬 張學軍2)? 孫知信

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

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

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

4) (南京郵電大學,國家郵政局郵政行業(yè)技術研發(fā)中心(物聯(lián)網(wǎng)技術),南京 210003))

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

1 引 言

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

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

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

2 方 法

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

2.1 彈性變分模態(tài)分解

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

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

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

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

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

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

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

迭代停止條件為：

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

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

A模態(tài)uk更新

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

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

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

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

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

B中心頻率ωk更新

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

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

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

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

Algorithm 1：EVMD

2.2 精細復合多尺度散布熵

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

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

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

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

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

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

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

步驟4若嵌入維數(shù)標記為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根據(jù)香農(nóng)熵的定義,信號uk的精細復合多尺度散布熵為

3 實驗結果與分析

3.1 實驗數(shù)據(jù)集

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

3.2 實驗

3.2.1 EEG 數(shù)據(jù)段分割

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

3.2.2 EVMD 分解

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

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

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

3.2.3 RCMDE 提取特征

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

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

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

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

表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 分類結果

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

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

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

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

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

4 結 論

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