吳端坡, 王紫萌, 蔣鐵甲, 董 芳, 馮 維, 劉兆霆
(1.杭州電子科技大學通信工程學院,杭州 310018;2.浙江大學醫學院附屬兒童醫院,杭州 310052;3.浙大城市學院信息與電氣工程學院,杭州 310015;4.浙江環瑪信息科技有限公司,杭州 310052)
隨著經濟的發展和社會的進步,人們對智慧醫療、智能醫療的需求日益增加,醫學單一學科發展已無法滿足人們的需求,隨著人工智能技術、機器人技術和納米載體技術等在醫學方面的廣泛應用,學科交叉與發展得到了加強,醫工結合應運而生[1]。醫工結合以醫學為基礎,以工學為手段,學科之間相互交叉與滲透,極大地推動了醫療水平發展[2]。
癲癇是大腦神經元突發性異常放電導致短暫的大腦功能障礙的一種慢性疾病,僅中國癲癇患者就有900百萬人之多[3-4]。棘波是癲癇患者腦電信號(Electroencephalogram,EEG)中的典型特征波,其波形尖銳,突出于背景波,具有高幅和瞬變的特性[5-6]。Nishida等[7]首次用數學形態學的方法檢測出EEG中的棘波,但沒有考慮單一開-閉或者閉-開運算引發的統計偏倚現象。隨機森林(Random Forest,RF)是貝爾實驗室的Tin Kam Ho在1995年提出的一種機器學習算法,包含了多個決策樹,其輸出結果由這些決策樹投票得到[8]。相對其他集成學習算法,RF具有極高的準確率,能夠處理大量高維度數據,抗過擬合能力強,常用于EEG信號檢測與分類研究[9]。
Matlab是矩陣運算與信號檢測領域常用的數據分析平臺,被科研技術人員、高校師生廣泛使用[10-11],《Matlab與仿真》是作者所在學院一直開設的專業必修課程,旨在培養學生的創新意識和規范編程的能力。本文在上述研究的基礎上,根據癲癇棘波的波形特征,對腦電信號進行形態學濾波,并結合RF分類模型,通過Matlab仿真實驗實現癲癇棘波自動檢測,提高學生的動手能力和獨立思考的能力。該實驗對教學硬件要求較低,能夠在高校大規模推廣。
如圖1所示,本文通過國際10-20電極分布系統采集真實病例的EEG信號[12]。該數據庫中EEG信號的采樣頻率為1 kHz,所有數據都已經由專業醫生進行了棘波放電標記。在EEG信號采集過程中,由于設備本身、肌肉活動和眨眼等因素容易產生噪聲和偽跡,在本研究中采用5階巴特沃斯帶通濾波器對EEG進行預處理,濾波范圍為1.6~70 Hz。

圖1 國際10-20電極分布系統
本文采用了5例患者的EEG信號進行試驗,分別命名為S1~S5。每個患者EEG信號的相關信息見表1。每個EEG數據中包含16個通道(Fp1、Fp2、F3、F4、F7、F8、T3、T4、T5、T6、C3、C4、P3、P4、O1和O2),本文分別使用其中能量最大的一個雙極導聯和平均導聯EEG進行實驗。圖2顯示了患者S1的部分癲癇棘波波形。

表1 訓練/測試數據集的相關信息

圖2 癲癇棘波信號
形態學濾波器是一種非線性濾波器,它利用信號的幾何特征將信號與預先設定的結構元素進行匹配,提取與結構元素相似的信號特征[13-14]。形態學濾波的基本操作包括腐蝕、膨脹、形態開運算和形態閉運算:

式中: f(n)為輸入的時間序列;g(n)為結構元素;f(n)、g(n)的序列長度分別為N和M,并且N>>M;?、⊕、。、·分別為腐蝕、膨脹、形態學開操作和形態學閉操作。
實驗以Matlab為仿真平臺,采用形態學濾波與RF模型的癲癇棘波檢測算法,主要包括候選棘波檢測、消除偽棘波、棘波形態特征提取和模型分類4部分,具體流程如圖3所示。

圖3 癲癇棘波檢測算法具體流程圖
為提高癲癇棘波自動檢測算法的準確性,采用形態學濾波的方法進行候選棘波檢測,具體過程如下。
對單通道EEG信號進行候選檢測的第1步就是選擇一個合適的結構元素,選擇能夠最大程度反映棘波信號幾何特征的三角形結構元素。三角形結構元素如圖4所示,其中:A為結構元素的中心高度;2L為結構元素的寬度。
為提高棘波提取的準確率,進一步的抑制背景信號,使用了一種改良的形態學濾波方法,設置2個結構元素,根據數據集中所有數據的棘波持續時間來確定結構元素的寬度,經統計,由醫生標記的癲癇棘波平均持續時間為60 ms,癲癇棘波最大波幅為結構元素高度A的范圍為[130,150]。對EEG信號進行候選所用的2個結構元素表示如下:

圖4 三角形結構元素結構

式中:A1和A2分別為130和150;2L的值為60 ms×1 kHz=60,因此本實驗中取L為30。
為同時去除信號中的正、負脈沖,消除誤差偏倚,構建一個基于形態開-閉(Open-closing,OC)和形態閉-開(Close-opening,CO)級聯的濾波器[15-16]:

式中, f (k)為患者EEG信號。
由于開操作具有擴展性而閉操作具有反擴展性,因此單獨進行OC操作或者CO操作會使結果出現統計偏倚現象。針對以上現象,使用平均加權的方法消除統計偏倚。EEG信號f(k)經過OC操作和CO操作的平均加權后,得到的輸出信號即為背景信號:

將原始EEG信號與背景信號p(k)相減得到候選棘波信號:

形態學濾波后得到的結果中具有很多偽棘波,對后面基于RF的棘波檢測結果有很大的影響,因此使用閾值判斷法消除候選棘波中的偽棘波。候選棘波EEG信號的閾值:

式中:μ為每一通道的候選棘波信號的平均值;σ為每一通道的候選棘波信號的方差。
消除偽棘波后,提取每個棘波10個形態學特征,以便進行模型訓練和分類。這10個特征分為4類,包括持續時間、振幅、斜率和面積。持續時間特征包括:棘波左半波的持續時間Dl、棘波右半波的持續時間Dr和棘波持續時間Ds;振幅特征包括:棘波左半波的峰值Al、棘波右半波的峰值Ar和棘波峰值As;斜率特征包括:棘波左半波斜率Sl、棘波右半波斜率Sr和棘波的銳度Ss;面積特征為Areas,圖5所示為特征提取的參考圖。其中:Sl=Dl/Al;Sr=-Dr/Ar;Ss=Sl-Sr;

圖5 特征提取參考圖
實驗中,選擇RF分類器進行棘波檢測。RF中每棵樹使用的數據集是隨機采樣的,對數據集的適應能力很強,可以減少數據不平衡的影響;對所有樹的分類結果進行加權平均,得到RF分類器的最終預測結果,避免分類模型的過擬合。由Breiman提出的RF分類器包含多個決策樹,其輸出類別由所有樹給出的最高票數決定。通過Bootstrap技術從原始訓練樣本集N中隨機選擇k個樣本來生成新的訓練樣本集,從k個單獨的決策樹分類器中生成RF。RF分類器的本質是對決策樹算法的改進。候選棘波的形態特征向量被輸入到RF分類器中,所有可能的棘波被分為癲癇棘波和非棘波。基于RF的棘波檢測流程如圖6所示。實驗中,將每個患者EEG信號的50%用于訓練模型,剩余50%用于模型測試。

圖6 基于RF的棘波檢測
為評估本實驗所用方法的有效性和優越性,通過計算混淆矩陣來體現該算法得到的檢測結果與專家標記結果之間的差異,混淆矩陣記為

式中:TN為檢測正確的非棘波數;FN為漏檢棘波數;TP為檢測正確的棘波數;FP為誤檢棘波數。
準確性、敏感性和特異性是評估癲癇棘波檢測算法性能的3個重要指標,分別記為AC、SE和SP,可由Cm計算得到,其中:AC為檢測正確的樣本數量與所有樣本數量的比值。SE為檢測正確的棘波數目與實際所有棘波數目之比,是衡量本文所提出方法對癲癇棘波的識別能力;SP為檢測正確的非棘波與實際所有非棘波數目之比,是衡量分類器檢驗負樣本的能力。

為保證仿真實驗所驗證的算法在實際醫療實踐中具有較低的FP,采用FPM來更好地反映所提出的方法的性能。FPM每分鐘的誤檢棘波數目:

式中:L為EEG信號記錄的總持續時間。
本實驗是通過Matlab2016a軟件實現的。根據不同的目的,使用表1所示的5份EEG數據進行實驗。
圖7所示為S1EEG數據經過形態學濾波候選檢測的結果,圖中t為EEG信號持續時間,Am為EEG信號的幅值。圖7(a)是患者的真實EEG信號;圖7(b)是經過形態學濾波后得到的背景信號;圖7(c)是提取到的候選棘波信號。可見本實驗使用的形態學濾波候選棘波檢測方法得到的結果有較高的準確性。
為評價該仿真實驗所用方法的通用性,使用AC、SP、SE和FPM對本文所提出的棘波檢測方法的性能進行了評價。表2為每個患者基于單通道雙極導聯EEG信號的棘波檢測結果。表中AC的值最小為90.7,最大為96.2。對于S1-S5的EEG數據,SE從71.8(S3)到97.2(S5),SP從93.4(S2)到97.3(S1),FP最低每分鐘出現0.99個。從表3中可以看出,S3的FP遠遠大于其他患者,通過視覺分析,S3中未檢測到的棘波的幅值與背景信號的幅值差距較小,因此易產生漏檢。

圖7 形態學濾波各階段結果

表2 基于雙極導聯EEG的癲癇棘波檢測結果

表3 基于平均導聯EEG信號的檢測結果
為驗證本實驗所用方法的可行性,表3給出了使用對應的平均導聯EEG信號的棘波檢測結果。比較表2、3的結果,發現使用雙極導聯EEG的棘波檢測模型的性能優于平均導聯EEG。經過分析,主要有以下2個原因:眼電圖(Electrooculogram,EOG)的波形與棘波相似,雙極導聯EEG信號可以抵消平均導聯上存在的EOG偽跡;癲癇患者的棘波放電通常出現在2個電極之間的區域,對平均導聯EEG做差可以增強棘波幅值,使形態學濾波結果更加準確。
表4是本實驗所用方法與其他方法性能指標的比較,所有方法均使用表1給出的5個數據進行驗證。文獻[17]中采用基于模板匹配的癲癇棘波檢測方法,其SE低于本實驗所用的方法。文獻[18]中提出了一種基于人工神經網絡模型的癲癇樣放電檢測方法,其SE僅有73.0,FPM比本實驗所用方法高4.04。文獻[19]中采用張量分解提取棘波特征,其性能略低于本實驗所用的方法。由表4可見,文獻[19]中的SE最高,為99.3,但SP僅為82.8。

表4 本方法與其他方法的比較
文獻[20]中的數據集包括健康人的EEG信號,本實驗所用的數據集中僅包括癲癇患者的EEG信號,因此文獻[20]中的SE會略高于其他方法。文獻[17]中具有最優的SP值,但其SE只有69.3。文獻[19]中的SE和SP差距較小,其性能指標低于本實驗所用的方法。臨床上,醫生對SE的要求高于SP。因此,該方法在保證較高SE的同時,也大大提高了SP。
本文設計的實驗是通過Matlab仿真平臺對表1給出的5個患者的EEG數據進行癲癇棘波檢測的研究,對每組數據的各項性能指標進行了統計和分析,將抽象復雜的數據處理結果通過直觀的圖像表示出來,并將該實驗應用于信號處理教學,達到了以下教學目的。
(1)了解癲癇發作特征波,掌握形態學濾波的基本原理和應用,深入理解RF分類模型的優勢。仿真結果直觀地體現了癲癇棘波的檢測結果,使抽象的算法原理變得直觀具體,有利于學生對Matlab及數字信號處理的學習興趣,對于數字信號處理及Matlab仿真的教學具有較大的參考價值。
(2)采用Matlab作為仿真環境,要求學生熟練使用Matlab函數庫解決實際問題,掌握基于Matlab的癲癇棘波智能檢測實驗的流程。
(3)通過對基于不同癲癇棘波智能檢測方法的仿真實驗結果進行比較,掌握不同分類器的訓練過程和異同點。
通過Matlab仿真實驗可以提高學生的邏輯思維能力和創新能力,為學生參與科研工作奠定理論基礎。在以后的仿真實驗中還可以深入分析發作間期與癲癇棘波之間的關系,以及癲癇樣波形的自動分類算法,為癲癇發作檢測、癲癇發作類型分類技術提供理論與實驗支持。