基于形態(tài)學(xué)濾波器與RF模型的棘波檢測仿真實驗

吳端坡， 王紫萌， 蔣鐵甲， 董 芳， 馮 維， 劉兆霆

（1.杭州電子科技大學(xué)通信工程學(xué)院，杭州 310018；2.浙江大學(xué)醫(yī)學(xué)院附屬兒童醫(yī)院，杭州 310052；3.浙大城市學(xué)院信息與電氣工程學(xué)院，杭州 310015；4.浙江環(huán)瑪信息科技有限公司，杭州 310052）

0 引 言

隨著經(jīng)濟的發(fā)展和社會的進(jìn)步，人們對智慧醫(yī)療、智能醫(yī)療的需求日益增加，醫(yī)學(xué)單一學(xué)科發(fā)展已無法滿足人們的需求，隨著人工智能技術(shù)、機器人技術(shù)和納米載體技術(shù)等在醫(yī)學(xué)方面的廣泛應(yīng)用，學(xué)科交叉與發(fā)展得到了加強，醫(yī)工結(jié)合應(yīng)運而生［1］。醫(yī)工結(jié)合以醫(yī)學(xué)為基礎(chǔ)，以工學(xué)為手段，學(xué)科之間相互交叉與滲透，極大地推動了醫(yī)療水平發(fā)展［2］。

癲癇是大腦神經(jīng)元突發(fā)性異常放電導(dǎo)致短暫的大腦功能障礙的一種慢性疾病，僅中國癲癇患者就有900百萬人之多［3-4］。棘波是癲癇患者腦電信號（Electroencephalogram，EEG）中的典型特征波，其波形尖銳，突出于背景波，具有高幅和瞬變的特性［5-6］。Nishida等［7］首次用數(shù)學(xué)形態(tài)學(xué)的方法檢測出EEG中的棘波，但沒有考慮單一開-閉或者閉-開運算引發(fā)的統(tǒng)計偏倚現(xiàn)象。隨機森林（Random Forest，RF）是貝爾實驗室的Tin Kam Ho在1995年提出的一種機器學(xué)習(xí)算法，包含了多個決策樹，其輸出結(jié)果由這些決策樹投票得到［8］。相對其他集成學(xué)習(xí)算法，RF具有極高的準(zhǔn)確率，能夠處理大量高維度數(shù)據(jù)，抗過擬合能力強，常用于EEG信號檢測與分類研究［9］。

Matlab是矩陣運算與信號檢測領(lǐng)域常用的數(shù)據(jù)分析平臺，被科研技術(shù)人員、高校師生廣泛使用［10-11］，《Matlab與仿真》是作者所在學(xué)院一直開設(shè)的專業(yè)必修課程，旨在培養(yǎng)學(xué)生的創(chuàng)新意識和規(guī)范編程的能力。本文在上述研究的基礎(chǔ)上，根據(jù)癲癇棘波的波形特征，對腦電信號進(jìn)行形態(tài)學(xué)濾波，并結(jié)合RF分類模型，通過Matlab仿真實驗實現(xiàn)癲癇棘波自動檢測，提高學(xué)生的動手能力和獨立思考的能力。該實驗對教學(xué)硬件要求較低，能夠在高校大規(guī)模推廣。

1 相關(guān)理論

1.1 數(shù)據(jù)描述與數(shù)據(jù)預(yù)處理

如圖1所示，本文通過國際10-20電極分布系統(tǒng)采集真實病例的EEG信號［12］。該數(shù)據(jù)庫中EEG信號的采樣頻率為1 kHz，所有數(shù)據(jù)都已經(jīng)由專業(yè)醫(yī)生進(jìn)行了棘波放電標(biāo)記。在EEG信號采集過程中，由于設(shè)備本身、肌肉活動和眨眼等因素容易產(chǎn)生噪聲和偽跡，在本研究中采用5階巴特沃斯帶通濾波器對EEG進(jìn)行預(yù)處理，濾波范圍為1.6～70 Hz。

圖1 國際10-20電極分布系統(tǒng)

本文采用了5例患者的EEG信號進(jìn)行試驗，分別命名為S1～S5。每個患者EEG信號的相關(guān)信息見表1。每個EEG數(shù)據(jù)中包含16個通道（Fp1、Fp2、F3、F4、F7、F8、T3、T4、T5、T6、C3、C4、P3、P4、O1和O2），本文分別使用其中能量最大的一個雙極導(dǎo)聯(lián)和平均導(dǎo)聯(lián)EEG進(jìn)行實驗。圖2顯示了患者S1的部分癲癇棘波波形。

表1 訓(xùn)練／測試數(shù)據(jù)集的相關(guān)信息

圖2 癲癇棘波信號

1.2 形態(tài)學(xué)濾波器

形態(tài)學(xué)濾波器是一種非線性濾波器，它利用信號的幾何特征將信號與預(yù)先設(shè)定的結(jié)構(gòu)元素進(jìn)行匹配，提取與結(jié)構(gòu)元素相似的信號特征［13-14］。形態(tài)學(xué)濾波的基本操作包括腐蝕、膨脹、形態(tài)開運算和形態(tài)閉運算：

式中： f（n）為輸入的時間序列；g（n）為結(jié)構(gòu)元素；f（n）、g（n）的序列長度分別為N和M，并且N＞＞M；?、⊕、。、·分別為腐蝕、膨脹、形態(tài)學(xué)開操作和形態(tài)學(xué)閉操作。

2 癲癇棘波檢測仿真實驗設(shè)計

實驗以Matlab為仿真平臺，采用形態(tài)學(xué)濾波與RF模型的癲癇棘波檢測算法，主要包括候選棘波檢測、消除偽棘波、棘波形態(tài)特征提取和模型分類4部分，具體流程如圖3所示。

圖3 癲癇棘波檢測算法具體流程圖

2.1 候選棘波檢測

為提高癲癇棘波自動檢測算法的準(zhǔn)確性，采用形態(tài)學(xué)濾波的方法進(jìn)行候選棘波檢測，具體過程如下。

對單通道EEG信號進(jìn)行候選檢測的第1步就是選擇一個合適的結(jié)構(gòu)元素，選擇能夠最大程度反映棘波信號幾何特征的三角形結(jié)構(gòu)元素。三角形結(jié)構(gòu)元素如圖4所示，其中：A為結(jié)構(gòu)元素的中心高度；2L為結(jié)構(gòu)元素的寬度。

為提高棘波提取的準(zhǔn)確率，進(jìn)一步的抑制背景信號，使用了一種改良的形態(tài)學(xué)濾波方法，設(shè)置2個結(jié)構(gòu)元素，根據(jù)數(shù)據(jù)集中所有數(shù)據(jù)的棘波持續(xù)時間來確定結(jié)構(gòu)元素的寬度，經(jīng)統(tǒng)計，由醫(yī)生標(biāo)記的癲癇棘波平均持續(xù)時間為60 ms，癲癇棘波最大波幅為結(jié)構(gòu)元素高度A的范圍為［130，150］。對EEG信號進(jìn)行候選所用的2個結(jié)構(gòu)元素表示如下：

圖4 三角形結(jié)構(gòu)元素結(jié)構(gòu)

式中：A1和A2分別為130和150；2L的值為60 ms×1 kHz＝60，因此本實驗中取L為30。

為同時去除信號中的正、負(fù)脈沖，消除誤差偏倚，構(gòu)建一個基于形態(tài)開-閉（Open-closing，OC）和形態(tài)閉-開（Close-opening，CO）級聯(lián)的濾波器［15-16］：

式中， f （k）為患者EEG信號。

由于開操作具有擴展性而閉操作具有反擴展性，因此單獨進(jìn)行OC操作或者CO操作會使結(jié)果出現(xiàn)統(tǒng)計偏倚現(xiàn)象。針對以上現(xiàn)象，使用平均加權(quán)的方法消除統(tǒng)計偏倚。EEG信號f（k）經(jīng)過OC操作和CO操作的平均加權(quán)后，得到的輸出信號即為背景信號：

將原始EEG信號與背景信號p（k）相減得到候選棘波信號：

2.2 消除偽棘波

形態(tài)學(xué)濾波后得到的結(jié)果中具有很多偽棘波，對后面基于RF的棘波檢測結(jié)果有很大的影響，因此使用閾值判斷法消除候選棘波中的偽棘波。候選棘波EEG信號的閾值：

式中：μ為每一通道的候選棘波信號的平均值；σ為每一通道的候選棘波信號的方差。

2.3 棘波形態(tài)特征提取

消除偽棘波后，提取每個棘波10個形態(tài)學(xué)特征，以便進(jìn)行模型訓(xùn)練和分類。這10個特征分為4類，包括持續(xù)時間、振幅、斜率和面積。持續(xù)時間特征包括：棘波左半波的持續(xù)時間Dl、棘波右半波的持續(xù)時間Dr和棘波持續(xù)時間Ds；振幅特征包括：棘波左半波的峰值A(chǔ)l、棘波右半波的峰值A(chǔ)r和棘波峰值A(chǔ)s；斜率特征包括：棘波左半波斜率Sl、棘波右半波斜率Sr和棘波的銳度Ss；面積特征為Areas，圖5所示為特征提取的參考圖。其中：Sl＝Dl/Al；Sr＝－Dr/Ar；Ss＝Sl－Sr；

圖5 特征提取參考圖

2.4 隨機森林分類

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

圖6 基于RF的棘波檢測

3 實驗結(jié)果與分析

3.1 性能評估指標(biāo)

為評估本實驗所用方法的有效性和優(yōu)越性，通過計算混淆矩陣來體現(xiàn)該算法得到的檢測結(jié)果與專家標(biāo)記結(jié)果之間的差異，混淆矩陣記為

式中：TN為檢測正確的非棘波數(shù)；FN為漏檢棘波數(shù)；TP為檢測正確的棘波數(shù)；FP為誤檢棘波數(shù)。

準(zhǔn)確性、敏感性和特異性是評估癲癇棘波檢測算法性能的3個重要指標(biāo)，分別記為AC、SE和SP，可由Cm計算得到，其中：AC為檢測正確的樣本數(shù)量與所有樣本數(shù)量的比值。SE為檢測正確的棘波數(shù)目與實際所有棘波數(shù)目之比，是衡量本文所提出方法對癲癇棘波的識別能力；SP為檢測正確的非棘波與實際所有非棘波數(shù)目之比，是衡量分類器檢驗負(fù)樣本的能力。

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

式中：L為EEG信號記錄的總持續(xù)時間。

3.2 結(jié)果分析

本實驗是通過Matlab2016a軟件實現(xiàn)的。根據(jù)不同的目的，使用表1所示的5份EEG數(shù)據(jù)進(jìn)行實驗。

圖7所示為S1EEG數(shù)據(jù)經(jīng)過形態(tài)學(xué)濾波候選檢測的結(jié)果，圖中t為EEG信號持續(xù)時間，Am為EEG信號的幅值。圖7（a）是患者的真實EEG信號；圖7（b）是經(jīng)過形態(tài)學(xué)濾波后得到的背景信號；圖7（c）是提取到的候選棘波信號?？梢姳緦嶒炇褂玫男螒B(tài)學(xué)濾波候選棘波檢測方法得到的結(jié)果有較高的準(zhǔn)確性。

為評價該仿真實驗所用方法的通用性，使用AC、SP、SE和FPM對本文所提出的棘波檢測方法的性能進(jìn)行了評價。表2為每個患者基于單通道雙極導(dǎo)聯(lián)EEG信號的棘波檢測結(jié)果。表中AC的值最小為90.7，最大為96.2。對于S1-S5的EEG數(shù)據(jù)，SE從71.8（S3）到97.2（S5），SP從93.4（S2）到97.3（S1），F(xiàn)P最低每分鐘出現(xiàn)0.99個。從表3中可以看出，S3的FP遠(yuǎn)遠(yuǎn)大于其他患者，通過視覺分析，S3中未檢測到的棘波的幅值與背景信號的幅值差距較小，因此易產(chǎn)生漏檢。

圖7 形態(tài)學(xué)濾波各階段結(jié)果

表2 基于雙極導(dǎo)聯(lián)EEG的癲癇棘波檢測結(jié)果

表3 基于平均導(dǎo)聯(lián)EEG信號的檢測結(jié)果

為驗證本實驗所用方法的可行性，表3給出了使用對應(yīng)的平均導(dǎo)聯(lián)EEG信號的棘波檢測結(jié)果。比較表2、3的結(jié)果，發(fā)現(xiàn)使用雙極導(dǎo)聯(lián)EEG的棘波檢測模型的性能優(yōu)于平均導(dǎo)聯(lián)EEG。經(jīng)過分析，主要有以下2個原因：眼電圖（Electrooculogram，EOG）的波形與棘波相似，雙極導(dǎo)聯(lián)EEG信號可以抵消平均導(dǎo)聯(lián)上存在的EOG偽跡；癲癇患者的棘波放電通常出現(xiàn)在2個電極之間的區(qū)域，對平均導(dǎo)聯(lián)EEG做差可以增強棘波幅值，使形態(tài)學(xué)濾波結(jié)果更加準(zhǔn)確。

表4是本實驗所用方法與其他方法性能指標(biāo)的比較，所有方法均使用表1給出的5個數(shù)據(jù)進(jìn)行驗證。文獻(xiàn)［17］中采用基于模板匹配的癲癇棘波檢測方法，其SE低于本實驗所用的方法。文獻(xiàn)［18］中提出了一種基于人工神經(jīng)網(wǎng)絡(luò)模型的癲癇樣放電檢測方法，其SE僅有73.0，F(xiàn)PM比本實驗所用方法高4.04。文獻(xiàn)［19］中采用張量分解提取棘波特征，其性能略低于本實驗所用的方法。由表4可見，文獻(xiàn)［19］中的SE最高，為99.3，但SP僅為82.8。

表4 本方法與其他方法的比較

文獻(xiàn)［20］中的數(shù)據(jù)集包括健康人的EEG信號，本實驗所用的數(shù)據(jù)集中僅包括癲癇患者的EEG信號，因此文獻(xiàn)［20］中的SE會略高于其他方法。文獻(xiàn)［17］中具有最優(yōu)的SP值，但其SE只有69.3。文獻(xiàn)［19］中的SE和SP差距較小，其性能指標(biāo)低于本實驗所用的方法。臨床上，醫(yī)生對SE的要求高于SP。因此，該方法在保證較高SE的同時，也大大提高了SP。

4 結(jié) 語

本文設(shè)計的實驗是通過Matlab仿真平臺對表1給出的5個患者的EEG數(shù)據(jù)進(jìn)行癲癇棘波檢測的研究，對每組數(shù)據(jù)的各項性能指標(biāo)進(jìn)行了統(tǒng)計和分析，將抽象復(fù)雜的數(shù)據(jù)處理結(jié)果通過直觀的圖像表示出來，并將該實驗應(yīng)用于信號處理教學(xué)，達(dá)到了以下教學(xué)目的。

（1）了解癲癇發(fā)作特征波，掌握形態(tài)學(xué)濾波的基本原理和應(yīng)用，深入理解RF分類模型的優(yōu)勢。仿真結(jié)果直觀地體現(xiàn)了癲癇棘波的檢測結(jié)果，使抽象的算法原理變得直觀具體，有利于學(xué)生對Matlab及數(shù)字信號處理的學(xué)習(xí)興趣，對于數(shù)字信號處理及Matlab仿真的教學(xué)具有較大的參考價值。

（2）采用Matlab作為仿真環(huán)境，要求學(xué)生熟練使用Matlab函數(shù)庫解決實際問題，掌握基于Matlab的癲癇棘波智能檢測實驗的流程。

（3）通過對基于不同癲癇棘波智能檢測方法的仿真實驗結(jié)果進(jìn)行比較，掌握不同分類器的訓(xùn)練過程和異同點。

通過Matlab仿真實驗可以提高學(xué)生的邏輯思維能力和創(chuàng)新能力，為學(xué)生參與科研工作奠定理論基礎(chǔ)。在以后的仿真實驗中還可以深入分析發(fā)作間期與癲癇棘波之間的關(guān)系，以及癲癇樣波形的自動分類算法，為癲癇發(fā)作檢測、癲癇發(fā)作類型分類技術(shù)提供理論與實驗支持。