秦喜文,郭 宇,董小剛,郭佳靜,袁 迪
(長春工業(yè)大學 a.數學與統(tǒng)計學院;b.研究生院,長春 130012)
癲癇是短暫和意外的腦電紊亂,是一種因腦部損傷而嚴重影響人類健康的腦部疾病。癲癇發(fā)作的跡象不僅可在癲癇患者的大腦中找到,也可在正常人的大腦中找到[1]。腦電信號是非線性、非平穩(wěn)的時序信號,可通過頭皮或顱內電極上的傳感器檢測,這些信號是神經元膜電位非常豐富的外部表現[2]。準確評估、術前評估、癲癇預防以及緊急警報都依賴于癲癇發(fā)作的快速檢測。醫(yī)生可通過監(jiān)測腦電信號評估大腦的狀態(tài)。因此,腦電信號在癲癇發(fā)作的檢測和鑒別中起著重要作用。傳統(tǒng)的癲癇檢測方法是神經學專家對長期腦電圖記錄的目視掃描。然而,這是一項耗時的任務,并且由于大量的腦電圖數據和不同神經學專家的臨床判斷標準不同,診斷可能不準確[3]。因此,開發(fā)智能型癲癇自動檢測系統(tǒng)顯得尤為迫切,具有重要的實際意義。
1875年,英國醫(yī)生Caton[4]發(fā)現動物大腦存在電活動現象。 1934年,經Adrian等[5]證實,人類的腦電信號獲得了世界各界的認同。在此之后,眾多學者開始對腦電圖進行相關研究,并把這部分知識應用于藥理學、生物學和醫(yī)學等諸多領域。Das等[6]從經驗模態(tài)分解和離散小波變換兩方面對局灶性和非局灶性腦電圖進行了綜合分析。Li等[7]提出了一種基于經驗模態(tài)分解和支持向量機的癲癇發(fā)作期腦電特征提取和模式識別方法。還有一些學者將時頻分析與機器學習相結合識別診斷腦電信號,Zhuang等[8]給出了一種基于經驗模態(tài)分解(EMD:Empirical Mode Decomposition)的特征提取和腦電信號識別方法;Chen等[9]提出了一種利用經驗模態(tài)分解和近似熵的腦電特征提取方法,設計了一種將深度信念網絡與支持向量機相結合的分類方法,提取特征向量,識別出人類快樂、平靜、悲傷和恐懼4種主要情緒;陸苗等[10]利用改進的經驗模態(tài)分解處理腦電信號,用所得固有模態(tài)函數(IMF:Intrinsic Mode Function)求取能量熵等特征值提高分類精度;韓敏等[11]針對單一極限學習機在癲癇腦電信號研究中分類結果不穩(wěn)定、泛化能力差的缺陷,提出一種基于互信息的AdaBoost極限學習機分類算法。隨著計算機科學技術的不斷發(fā)展,腦電信號的分析方法層出不窮,越來越多樣化,陸續(xù)出現了非線性分析法[12]、時頻分析法[13-14]、高階譜分析法[15]等系列現代化分析方法[16]。筆者針對非平穩(wěn)、非線性的腦電信號,提出一種基于局部均值分解和迭代隨機森林相結合的腦電信號分類方法。利用局部均值分解,將腦電信號分解成若干個乘積函數分量和一個殘余分量,然后對分解后的所有分量進行特征提取,最后使用隨機森林和迭代隨機森林等方法進行分類。
局部均值分解(LMD:Local Mean Decomposition)方法可逐步地將調頻信號與調幅包絡信號分離。通過平滑處理原始信號,將平滑后的信號從原始信號中減去,然后對包絡函數進行估計,對結果進行幅度解調實現分離。通過使用由原始信號的連續(xù)極值之間的時間延遲加權的移動平均可獲得包絡估計和原始信號的平滑版本。如果得到的信號不具有平坦包絡,則重復該過程,直到獲得具有平坦包絡的調頻信號。然后,從調頻信號中計算瞬時頻率,將所有包絡估計函數相乘得最終包絡,將該包絡乘以頻率調制信號以形成乘積函數,從原始信號中減去該乘積函數,對得到的信號重復整個過程,以產生具有相關包絡和瞬時頻率的第2乘積函數。繼續(xù)該分解,直到剩余信號不再包含振蕩,得到一個單調函數為止。
LMD算法步驟如下[17]。
1)從原信號中記錄全部局部極值點ni,對信號中相鄰的2個局部極值點取均值
(1)
將所有局部均值mi用直線串聯起來,然后利用滑動平均法對局部均值線段進行平滑處理,得局部均值函數m11(t)。
2)利用局部均值點,求出倆相鄰極值點之間的包絡估計值
ai=|ni-ni+1|/2
(2)
將全部相鄰2個包絡估計值ai用直線連接,然后采用滑動平均方法進行平滑處理,得包絡估計函數a11(t)。
3)將局部均值函數m11(t)從原信號x(t)中剔除,得
h11(t)=x(t)-m11(t)
(3)
4)對h11(t)進行解調,用h11(t)除以包絡估計函數a11(t),得
s11(t)=h11(t)/a11(t)
(4)
計算s11(t)的包絡估計函數a12(t),若a12(t)≠1,說明s11(t)并不是一個純調頻信號,將重復上述解調過程,直到使s1n(t)成為一個純調頻信號為止,即-1≤S1n(t)≤1,且其包絡函數a1(n+1)(t)=1,有
(5)
其中
(6)
迭代終止條件為
(7)
a1n(t)≈1
(8)
5)在上述迭代過程中,將求得的所有包絡估計函數相乘,得包絡信號
(9)
6)用包絡信號a1(t)和純調頻信號s1n(t)相乘,可得原信號的乘積函數(PF:Product Function)
P1(t)=a1(t)s1n(t)
(10)
7)將一個分量P1(t)從原信號x(t)中剔除,獲得一個全新的信號u1(t),再將u1(t)作為原始數據重復上述步驟,循環(huán)k次,直到uk為一個單調函數為止
(11)
所有的乘積函數分量和uk重組,可求出原始信號
(12)
這說明LMD分解沒有造成原信號的丟失。
對局部均值分解后得到的各階乘積函數提取重要特征,其中包括均值、標準差等一系列簡單統(tǒng)計量,以及能量熵和信息熵。
1.2.1 變異系數
變異系數(CV:Coefficient of Variation)是衡量資料中各觀測值變異程度的統(tǒng)計量,為標準差與平均數的比值
C=σ/μ
(13)
其中
其中l(wèi)為數據長度,即分解后的乘積函數長度。
1.2.2 波動指數
波動指數是描述信號波動幅度大小的一個統(tǒng)計量,在實驗中可見癲癇發(fā)作間期信號的波動劇烈程度小于發(fā)作期,且相對穩(wěn)定。其定義為
(14)
1.2.3 能量熵
在不同狀態(tài)下,腦電振動信號的頻率分布會發(fā)生變化,故癲癇發(fā)作期與發(fā)作間期腦電信號的能量也會發(fā)生變化。一般情況下,發(fā)作間期腦電信號中的能量分布是均勻的、不確定的,因而大于癲癇發(fā)作期腦電信號。因此可用LMD能量熵判斷患者是否處于癲癇狀態(tài)。
乘積函數分量的能量熵通過以下步驟計算[18]。首先,計算第i個PF的能量
(15)
其中t1和t2分別為信號開始和結束時間。之后計算整個PF的能量
(16)
最后,定義整個PF的能量熵為
(17)
1.2.4 信息熵
信息熵是一個數學的抽象概念,是為了消除數據不確定性的度量。癲癇發(fā)作期信號的信息熵通常比發(fā)作間期的信息熵低。其定義為
H(x)=E[I(x)]=E[log2(1/p(xi))]=-∑p(xi)log2(p(xi)),i=1,…,n
(18)
其中x表示隨機變量,與之相對應的是所有可能輸出的集合,定義為符號集;p(x)表示輸出概率函數。變量的不確定性越大,熵也就越大。
利用LMD分解波恩癲癇研究中心數據庫中癲癇發(fā)作間期與發(fā)作期數據,計算分解后乘積函數的特征值。
迭代隨機森林(iRF:iterative Random Forest)在隨機森林的基礎上通過對選定的特征進行迭代重新賦權(iterative re-weighting),得到一個帶有特征權重的隨機森林,然后將泛化的隨機交叉樹作用于帶有特征權重的隨機森林上,識別特征的高階交互作用,保證其具有較高的預測能力。迭代隨機森林的具體工作流程如下[18]。
1)對隨機森林進行迭代重新賦權,給定迭代次數K,在數據集D上,通過迭代生成K個具有特征權重的隨機森林,記為RF(ωk),k=1,…,K。初始權重k=1時,ω=(1/p,…,1/p),并計算p個特征的重要性,即基尼純度的平均減少,記為υ(1)=(υ1,…,υp)。令ω(k)=υ(k-1),即用隨機森林特征的重要性作為其權重。
2)將泛化的隨機交叉樹作用于RF(ωK),產生一組交叉作用集S。
3)計算Bagged穩(wěn)定得分,使用“外層”(out layer)自助法,以評價重現交叉作用的穩(wěn)定性。生成自助抽樣數據集D(b),b=1,…,B,在D(b)上擬合隨機森林RF(ωK),并且使用泛化隨機交叉樹識別交互作用集S(b),給出交叉作用集S的穩(wěn)定分數公式
(19)
實驗數據取自德國波恩大學癲癇研究中心數據庫。該數據庫是目前癲癇腦電信號研究最為普遍采用的腦電數據庫,其中包含5個子集,分別標號為A、B、C、D、E。每個子集包含100個實驗樣本,每個樣本采樣時間為23.6 s、頻率為173.61 Hz,每個樣本信號包含4 097個樣本點。筆者選取數據庫中數據集D和數據集E,分為癲癇發(fā)作間期和發(fā)作期進行分析。圖1給出了發(fā)作間期與發(fā)作期的原始癲癇腦電信號。
a 癲癇發(fā)作間期 b 癲癇發(fā)作期
圖1a為一組癲癇發(fā)作間期的腦電信號樣本,可見其信號波動程度較大,但具有一定的規(guī)律性。圖1b為一組癲癇發(fā)作期的腦電信號樣本,可見該信號波幅較小,且不穩(wěn)定。由圖1可知,兩種腦電信號在波動幅度和穩(wěn)定性方面具有明顯差異。
圖2和圖3給出了LMD分解后癲癇發(fā)作期和發(fā)作間期的腦電信號PF分量。癲癇發(fā)作間期腦電信號經LMD分解后得到4個乘積函數和1個殘余分量(見圖2);癲癇發(fā)作期腦電信號經LMD分解后得到4個乘積函數和1個殘余分量(見圖3)。
圖2 發(fā)作間期EEG分解后的PF分量 圖3 發(fā)作期EEG分解后的PF分量Fig.2 PF of EEG after decomposition in interictal period Fig.3 PF of EEG after decomposition in ictal period
對癲癇腦電信號發(fā)作間期和發(fā)作期的200個樣本信號,經過局部均值分解后提取所有乘積函數的特征向量。表1為癲癇腦電信號發(fā)作間期的特征樣本,表2為癲癇腦電信號發(fā)作期的特征樣本。
表1 癲癇腦電信號發(fā)作間期特征樣本
表2 癲癇腦電信號發(fā)作期特征樣本
筆者分別選取支持向量機、隨機森林及迭代隨機森林3種機器學習方法對局部均值分解后所求特征向量進行分類,并比較分類結果(見表3)。實驗樣本共1 284個,選取900個作為訓練集,另外384個作為測試集,結果表明,LMD與迭代隨機森林相結合的算法對癲癇發(fā)作間期和發(fā)作期EEG的乘積函數分類效果最佳,正確識別率為98.177%。
表3 癲癇腦電信號分類結果比較
在臨床醫(yī)學上,癲癇腦電信號的自動檢測和分類有著重要意義。Zhang等[19]提出了一種基于LMD和SVM結合的腦電信號分類方法,證明了在診斷腦電信號方面,腦電信號自適應分解與機器學習相結合的全新算法具有可行性和有效性。
在處理非平穩(wěn)和非線性信號時,LMD相對于其他時頻分析方法具有很多優(yōu)勢。LMD采用平滑處理的方法形成局部均值函數和局部包絡函數,這既保持了EMD分解的優(yōu)勢,又可避免像EMD分解中采用3次樣條函數形成上下包絡時產生的過包絡、欠包絡現象,此外LMD端點效應的擴散程度遠小于EMD。此方法除了分解腦電信號,還可應用于其他領域,如機械故障、金融股市和環(huán)境質量等。在目前眾多分類算法中,iRF不僅具有和RF一樣的抗噪音和抗干擾能力,且經過多次迭代賦權,可充分利用每棵樹的分類能力,使整體分類效果得到提升。
通過對德國波恩大學癲癇研究中心的腦電數據進行實驗分析,將局部均值分解方法的自適應分解特性與迭代隨機森林的多次賦權相結合,使LMD-iRF算法分類結果的準確率高于LMD-SVM和LMD-RF算法。此方法不僅可用于生物醫(yī)學信號分析處理,還可應用在環(huán)境監(jiān)測、機械故障診斷和金融分析等領域,是一種具有較高精度的分類算法。
筆者利用局部均值分解處理癲癇腦電信號,并對分解的信號進行特征提取,再對所提取特征使用迭代隨機森林方法進行分類,完成對癲癇腦電信號的自動檢測。實驗結果表明,該方法可對癲癇發(fā)作期和發(fā)作間期的腦信號準確分類,分類準確度高達98.177%,為癲癇的快速準確檢測奠定了良好的基礎。