正確區分癲癇發作期（seizure）與非發作期（non-seizure）對癲癇治療有著重要意義。本研究以顱內腦電信號（iEEG）作為研究對象，提出了一種基于雙密度雙樹復小波變換（DD-DT CWT）的癲癇發作期自動檢測算法。實驗數據來自美國國家衛生研究所（NINDS）設立在Kaggle上的15 719個競賽數據，處理后的數據庫由55 023段發作期數據和501 990段非發作期數據組成，每段數據長1 s，包含174個采樣點。首先對信號進行重采樣；然后利用DD-DT CWT處理腦電信號，并從中提取小波熵、方差、能量和均值共四類特征；最后使用最小二乘支持向量機（LS-SVM）學習分類，并通過比較不同小波分解層數下的實驗結果選取合適的分解層數。實驗結果表明：所提取的四類特征在發作期與非發作期存在差異，八位患者中，采用3層分解時分類的平均準確率較高，達到91.98%，靈敏度為90.15%，特異性為93.81%。本文工作表明，我們提出的算法在癲癇患者腦電信號的二分類中有優良的性能，能夠自動高效地檢測出癲癇發作期。
引用本文： 康同舟, 左潤東, 鐘嵐烽, 陳文靜, 張恒, 劉紅秀, 賴大坤. 基于雙密度雙樹復小波變換的癲癇發作期自動檢測算法. 生物醫學工程學雜志, 2021, 38(6): 1035-1042. doi: 10.7507/1001-5515.202105075 復制
引言
癲癇（epilepsy）是一種常見的神經系統疾病，目前全球范圍內共有五千多萬癲癇患者。癲癇發作時神經元異常放電引起大腦電活動的同步性增加，導致患者短暫的功能障礙[1]。癲癇不僅會造成患者生理障礙，同時會影響患者的心理健康，而且患者過早死亡率是普通人群的三倍。使用頭皮電極采集腦電圖（electroencephalogram，EEG）并結合磁共振成像、計算機斷層掃描等影像學檢查掃描大腦結構是癲癇術前診斷的常規方式。近年來，隨著新型腦電檢測技術的發展，記錄腦電信號的方式越來越精準，尤其是顱內腦電圖（intracranial EEG，iEEG）的應用，可以直接測量患者腦部異常放電情況[2]。然而，腦電圖數據量巨大，傳統的人工診斷方式會造成醫療資源的浪費，導致大量癲癇患者得不到有效的治療。通過設計性能優越的檢測算法，對患者的癲癇發作期進行自動檢測，從而協助醫生高效率高質量地對癲癇患者進行診斷，對癲癇的治療有著重要的臨床意義。
基于腦電信號的癲癇發作期自動檢測算法基本步驟是通過對原始腦電信號進行處理，提取相應特征并結合機器學習方法進行分類。小波變換因其具有良好的時頻局部化特性、多分辨率特性、去相關性等特點，在生物醫學信號處理中得到廣泛應用[3]。2000年Kingsbury[4]提出了雙樹復小波變換（dual tree complex wavelet transform，DT CWT）的概念，該小波變換具有近似位移不變性和良好方向選擇性；隨后，Selesnick等[5]提出了雙密度小波變換（double density CWT，DD CWT）的理論，并于2004年[6]提出了雙密度雙樹復小波變換（dual density dual tree CWT，DD-DT CWT）的概念。2014年，Yang等[7]使用DT CWT識別EEG信號中的運動片段，取得了91.07%的準確率。2015年，Roy等[8]使用DD CWT去除多通道EEG信號中的噪聲，結果表明在去除EEG偽影中，DD CWT方法優于獨立成分分析。上述結果證明DT CWT以及DD CWT可以用于腦電信號的識別以及噪聲處理。DD-DT CWT綜合了二者的優點，相比于DT CWT具有更好的重構性、抗混疊性和平移不變性，目前在圖像分割[9]、圖像增強[10]、信號去噪[11]等圖像處理領域和人臉識別領域[12]已有廣泛應用。
1999年，Suykens等[13]在支持向量機（support vector machine，SVM）基礎上改進得到最小二乘支持向量機（least squares SVM，LS-SVM），此方法是基于SVM理論的單隱層多層感知器分類訓練方法。Sharma等[14]使用經驗模態分解法（empirical mode decomposition method，EMD）對腦電信號進行處理，利用LS-SVM作為分類器，取得了84%的準確率和靈敏度。小波熵可以表征信號在不同分解尺度上的能量分布[15]特征。2006年Natwong等[16]從高分辨率心電圖中提取小波熵對心肌梗死患者進行檢測。Deivasigamani等[17]首先使用DT CWT對腦電信號進行分解，選取了多個特征值進行分析，并利用自適應網絡模糊推理系統分類，實驗結果準確率最高可達99%。Xiang等[18]使用模糊熵對癲癇發作和正常腦電圖信號進行檢測和分類，提取不同狀態下的腦電信號特征，并通過SVM進一步訓練和分類，對非發作期和發作期腦電信號分類的準確率為88.55%，靈敏度為90.36%，特異性為87.63%。Bhattacharyya等[19]在2017年提出了一種基于可調Q小波變換（tunable-Q wavelet transform，TQWT）的多變量子帶熵度量方法，并將特征輸入到隨機森林和LS-SVM分類，選擇最優分類，得到最高84.67%的準確率。
綜上所示，我們可以得知國內外對于癲癇自動檢測算法的研究已趨于成熟。但仍存在以下問題：① 腦電信號十分微弱，一般在微伏數量級，易受偽差干擾，尋找穩定合適的信號處理方法至關重要。② 準確的特征可以更好地反映癲癇發作期與非發作期EEG信號的區別，許多研究采用單一特征值，降低了信號的可區分性，如何選擇特征以及如何判斷特征在發作期與非發作期之間是否存在差異是類似研究中的難點。本文針對上述問題，設計并實現了基于機器學習的自動檢測算法，完成了基于患者腦部信號的癲癇發作期與非發作期的二分類。
1 方法
圖1為本文提出的基于DD-DT CWT的癲癇發作期自動檢測算法實施框圖。本方法包括如下四個步驟：① 腦電信號的預處理，將腦電數據處理為同一采樣頻率以及相同的時間片段，并將數據轉化為.mat格式方便后續分析；② 對信號進行DD-DT CWT；③ 從經過DD-DT CWT后的信號中提取小波熵、均值、方差以及能量四類特征；④ 將提取的四類特征送入LS-SVM中進行分類并進行算法性能評估。
圖1
基于DD-DT CWT的癲癇發作期自動檢測算法實施框圖
Figure1.
Principle block diagram of automatic seizure detection algorithm based on dual density dual tree complex wavelet transform
1.1 數據來源與預處理
本文數據來源于公開數據庫UPenn and Mayo Clinic’s Seizure Detection Challenge[20]。該數據庫為美國國家衛生研究所（National Institues of Health，NINDS）設立在 Kaggle 上的競賽數據。原始數據來自8名已接受癲癇手術評估的顳葉和顳葉外癲癇患者，數據集分為兩部分：訓練集（training-set）和測試集（test-set）。測量時使用EEG記錄來自沿海馬前后軸植入的深度電極和不同位置的硬膜下電極網格，癲癇腦電數據分為發作期和非發作期，前者被標注為seizure，后者被標注為non-seizure。發作期以及非發作期數據如表1所示。
表1
發作期與非發作期數據統計
Table1.
Statistics of seizure and non-seizure data
所有8名患者中癲癇發作期數據1 390個，癲癇非發作期數據14 329個，為了方便后續處理，將各個數據按固定時間為一段進行截取，從而得到新的數據集。本研究選擇了常用的數據長度和采樣頻率，綜合考慮后采取了將所有數據采樣頻率降至174 Hz的方案，即每秒數據包含174個數據點，同時按照每段數據1 s進行劃分。發作期和非發作期數據段數均足夠后續的訓練與測試，處理后8名患者中癲癇發作期數據為55 023段，癲癇非發作期數據為501 990段。
1.2 腦電一維DD-DT CWT
從癲癇發作間期到發作期，腦電圖信號的頻率成分通常會發生變化。為了分析這些變化，對腦電信號的分類提供有用的信息，我們對腦電信號采用了DD-DT CWT。DD-DT CWT從離散小波變換的基礎上改進而來，區別于后者，它具有一個尺度函數和兩個小波函數，兩個小波函數之間的偏移量為半個單位。DD-DT CWT通過兩個高通濾波器和一個低通濾波提取腦電信號的高頻成分和低頻成分。本文中，我們定義 
為尺度函數，用于分解原始腦電信號的低頻成分，同時對應的低通濾波器為 
。定義小波函數 
和小波函數 
分解原始信號的高頻成分，對應的高通濾波器分別為 
和 
。由于相鄰小波頻帶間隔小，降低了平移敏感性，因此該小波變換有著更好的近似平移不變性。
表征腦電信號低頻特征對應的尺度方程：
 |
表征腦電信號高頻特征對應的小波方程（其中i =1, 2）：
 |
經過變換的原始腦電信號分為實部和虛部兩部分，分別由樹 A和樹 B表示。樹A對應第一個雙密度小波變換，使用了兩個高通濾波器 
和 
來表示腦電信號高頻系數，一個低通濾波器 
來表示腦電信號低頻系數。樹B對應第二個雙密度小波變換，使用兩個高通濾波器 
和 
來表示腦電信號高頻系數，一個低通濾波器 
來表示腦電信號低頻系數。記兩個樹對應的尺度函數分別為 
和 
，四個小波函數分別為 
、
、
、
。
尺度函數和小波函數共同構成了如下的復尺度函數和復小波函數，并協同分解原始腦電信號：
 |
 |
 |
以對原始腦電信號3層分解尺度為例，可以得到如圖2的分解示意圖。X(n)表示原始信號，沿樹A和樹B方向分別進行實部和虛部分解。在分解1層時，我們得到2條低頻子帶 
、
，以及4條高頻子帶 
、
、
、
；當進行第2層分解時，我們將現有的2條低頻子帶繼續分解，產生了新的6條子帶，同時保留了1層分解時產生的4條高頻子帶；在3層分解時，對第2層分解所得低頻子帶進行分解，產生新的6條子帶，保留1層、2層所得高頻子帶。最終獲得了12條腦電高頻子帶和2條腦電低頻子帶分別表示信號中的高頻成分和低頻成分。
圖2
3層分解尺度下雙密度雙樹復小波分解示意圖
Figure2.
Schematic diagram of dual density dual tree complex wavelet decomposition under three-level decomposition scale
1.3 腦電一維DD-DT CWT特征提取
基于傳統信息熵理論，本文提出如下基于原始腦電信號一維DD-DT CWT的小波熵函數。首先對原始腦電信號進行DD-DT CWT變換，第m個分解層數、第n個子帶 
對應的小波能量 
如式(6)所示：
 |
其中N代表當前子帶下小波系數個數，
為小波系數。
J為分解層數，則總能量為：
 |
相對小波能量定義為：
 |
根據Shannon熵的計算公式和定義，此處小波熵為：
 |
除了小波熵，本文結合腦電信號的一維DD-DT CWT與傳統特征提取方法，提取的其他三種特征分別為均值（Mean）、方差（SD）和能量（Energy）
均值（Mean）計算公式如式（10）所示：
 |
其中，h代表低頻子帶系數，g代表高頻子帶系數，a代表低頻子帶系數個數，b代表高頻子帶系數個數
方差（SD）計算公式如式（11）所示：
 |
其中SD低和SD高分別代表低頻子帶方差和高頻子帶方差，n代表方差總數。
能量（Energy）計算公式如式（12）所示：
 |
腦電信號經過3層分解尺度下的DD-DT CWT后，可以得到12條高頻子帶和2條低頻子帶，四種類型特征值對應的數目分別為：均值4個，方差4個，能量4個，小波熵4個，共計16個特征值。
1.4 LS-SVM分類器
LS-SVM是基于SVM理論上的單隱層多層感知器分類訓練方法。相比于SVM使用一個超平面將兩組數據分開，LS-SVM回歸采用一個超平面對已知數據進行擬合，該分類器的核心思想是將數據從輸入空間映射到特征空間中，然后在這個特征空間里面構造最優的決策函數，使得最終模型的輸出誤差達到最小化。其最大優勢是運行速度遠高于SVM，但其解不具有稀疏性。
對于訓練集，其中輸入的特征值共包含N維數據，在本文中即為16維特征值，標記分類結果，取值為1（陽性）或 ?1（陰性），分別代表癲癇發作期與非發作期。LS-SVM通過如下公式求解最優：
 |
w為系數向量，確定了一個在維度為d的空間中的分類超平面，b為一個常量，γ為懲罰系數，
為分類誤差，求解以上最小化函數可以得到：
 |
其中 
表示核函數，計算 
和 
在特征空間中的內積。核函數的主要作用是將非線性可分問題轉換為線性可分問題。本文使用的核函數為式（15）所示的徑向基函數核：
 |
2 結果與討論
2.1 評估標準
本文根據文獻[21-24]將腦電信號劃分為1 s片段，最終得到了癲癇發作期數據55 023段以及非發作期數據501 990段，該數據集符合臨床樣本比例。本文采用十折交叉驗證法對算法進行性能評估，將訓練集隨機分成十等分，九份樣本作為訓練集，一份作為測試集，以體現算法的整體性能。本文采用二分類機器學習中腦部疾病檢測常用的評價指標：靈敏度（sensitivity，Sen）、特異性（specificity，Sp）、準確率（accuracy，Acc），分別如式（16）~（18）所示。
 |
 |
 |
其中，TP代表腦電數據真實情況為癲癇發作期，分類器判定為癲癇發作期的數據個數；FP代表腦電數據真實情況為癲癇非發作期，分類器判定為癲癇發作期的數據個數；FN代表腦電數據真實情況為癲癇發作期，分類器判定為癲癇非發作期的數據個數；TN代表腦電數據真實情況為癲癇非發作期，分類器判定為癲癇發作期的數據個數。
2.2 癲癇發作期與非發作期特征差異性分析
本文針對腦電信號微弱、易受干擾以及個體差異大等問題，將DD-DT CWT用于一維腦電信號處理。在分解層數J為3時隨機抽取部分連續數據計算各類特征值，其對比圖如圖3所示，紅色為發作期，藍色為非發作期，橫坐標表示對應的樣本，縱坐標表示樣本對應的特征值。由圖3可見，癲癇患者的發作期與非發作期腦電信號特征提取后所得的均值、方差、能量、小波熵均具有一定的差異性。
圖3
發作期和非發作期各樣本特征值對比
Figure3.
Comparison of feature values between seizure and non-seizure periods
為了進一步對比四類特征中每個特征的分類效果，本文采用t檢驗對特征進行統計學檢驗，如表2所示。首先對八位患者的腦電數據進行DD-DT CWT，分解層數J設定為3層；接著計算各類特征對應的特征值并形成特征矩陣；最后在癲癇發作期與非發作期之間對每個特征進行t檢驗，檢驗水準為0.05。每類特征包含4維特征矩陣，對于均值、方差以及能量三類特征，其中第1維特征代表第3層低頻子帶均值，其余3維特征代表不同分解層數下的高頻子帶均值。分析表2中各特征P值，僅第1層高頻子帶均值以及高頻子帶 
小波熵得到的P值略大于0.05，其余特征值得到的P值均小于0.05，即這些特征在發作期與非發作期之間的差異均具有統計學意義。上述結果證明本文選取的特征能夠用于區分癲癇發作期與非發作期之間的腦電信號。本文結合多種特征值完成對癲癇自動檢測算法的設計，為癲癇的特征提取方法提供了新的研究思路。
表2
分解層數J為3時各特征統計學檢驗的P值
Table2.
Statistical test for each feature when decomposition level J is 3
2.3 最優分解層數的探討
由于該小波變換分解層數的可調性，不同分解尺度代表著對信號高頻細節的提取程度不同，也會進一步影響到算法的工作效率和實驗結果。為確定該算法最合適的分解層數，本文選用了相同的訓練集數據和測試集數據在不同分解尺度下進行實驗，比較實驗所得結果。由于分解層數的更換，四種特征值對應的特征值數目會發生改變，如表3所示。
表3
不同分解尺度下特征值數目
Table3.
Number of eigenvalues under different decomposition scales
隨著分解層數的增多，由于分解后所得子帶數目的增加，前三種特征值（均值、方差、能量）的個數每層增加1，而小波熵的特征值數目不隨分解層數的增加而增加。
通過調整算法中分解層數參數（J = 1, 2, 3, 4）運行算法，我們得出了分解層數J分別為1、2、3、4的平均實驗結果以及每位患者分類層數為3時的分類結果，如表4、表5所示。
表4
不同分解層數分類器平均測試結果以及總耗時
Table4.
Average test results and total time consumption of clas sifiers with different decomposition levels
表5
分解層數為3時每位患者測試結果
Table5.
Test results of each patient when the number of decom position layers is 3
從以上各層分解實驗結果可得每位患者的最優分解層數，如表6所示。整體上看，分解層數越多，算法的性能就越好。但由于患者之間的個體差異，不同患者采用該算法的效果有一定差距，為了更完善地討論如何選用最合適的分解層數，我們認為應當從兩方面入手進一步分析。第一是關注不同層數平均準確率；第二是應當重視各個患者準確率達到最大值時算法所處的分解層數。
表6
每位患者的最優分解層數
Table6.
Optimal decomposition levels per patient
基于以上兩點，共有7名患者在選用算法分解層數J為3或4時準確率最大，所以分解層數J應考慮設置為3或4。就整體的平均準確率而言，分解層數J為3或4層時的準確率相差0.43%。從實際情況出發，分解層數每增加一層，特征維度增加3，進而導致算法的運算復雜度增加。通過分析表4可知，隨著分解層數的增加，算法運行時間增加，如果需要檢測的腦電信號數據量龐大（例如很多患者有長達數小時的腦電數據），設計算法時需要優先考慮算法運行過程的參數量以及運行效率。相比于分解層數J為1、2、3層，在4層分解時，分解層數的增大和特征值數目的增加都會在不同程度上增加算法的計算量與運行時間，從而造成運算資源浪費，因此選用3層分解層數更具有實際意義，而從追求算法性能角度來說，4層分解層數下本算法的性能最為優異。
2.4 癲癇發作期自動檢測算法的性能對比
如表7所示，關于癲癇發作期自動檢測的算法主要包括文獻[14]、[19]、[25]、[26]，這些算法均使用了LS-SVM作為分類器，同時選用的特征值均含有熵。本文實驗結果所顯示的準確率、靈敏度和特異性均高于這四項研究，證明了本文提出的基于一維DD-DT CWT小波熵的合理性。值得注意的是，本方法與其他作者所使用的數據集不同，但本方法所用數據庫數據量（發作期55 023段，非發作期501 990段，每段1 s）大于其他作者的數據庫數據量（發作期3 750段，非發作期3 750段，每段20 s）。同時，大部分患者一次癲癇發作持續時間小于兩分鐘，數據段劃分過長可能會導致本算法在應用于臨床數據時檢測出的發作期存在部分非發作期片段，因此本文綜合考慮了其他研究[21-23]的時間劃分方式，最終決定選擇將數據劃分為1 s片段。另一方面，癲癇患者腦電信號中非發作期時長遠大于發作期時長[24]，數據集符合真實樣本比例。本文將該小波變換用于對一維腦電信號的處理，實驗結果的成功證明該方法對于腦電信號分析的可行性。
表7
本文方法與近期國際報道的算法性能對比表
Table7.
Comparison of the performance of this method and the algorithm reported recently
盡管本文基于DD-DT CWT實現了對腦電癲癇發作期信號的自動識別，然而本文算法的結果是在樣本比例接近臨床情況下不區分患者的結果，主要適用針對某個患者長期監測的場景。隨著后續數據采集過程中患者數量的增加，還有待探索泛化性能更佳的算法，使得算法性能提升的同時滿足不同患者的檢測需求，并減少因患者個體差異而造成的分類結果不穩定問題。
3 結論
本文提出了一種基于DD-DT CWT的顱內腦癲癇發作自動檢測算法，準確地實現了對癲癇發作期的自動識別與分類。本方法利用DD-DT CWT對經過預處理的信號進行分析，并從中提取了四類特征（小波熵、方差、能量以及均值），t檢驗結果表明在癲癇發作期與非發作期之間四類特征存在差異，采用LS-SVM對特征進行分類，可完成對癲癇發作期的準確識別。相較于已報道的癲癇發作自動檢測算法，本文方法在保證較高檢測精度的基礎上，綜合應用了多種腦電特征。本文工作可輔助臨床醫生劃定癲癇發作期與非發作期，具有潛在的臨床應用價值。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
引言
癲癇（epilepsy）是一種常見的神經系統疾病，目前全球范圍內共有五千多萬癲癇患者。癲癇發作時神經元異常放電引起大腦電活動的同步性增加，導致患者短暫的功能障礙[1]。癲癇不僅會造成患者生理障礙，同時會影響患者的心理健康，而且患者過早死亡率是普通人群的三倍。使用頭皮電極采集腦電圖（electroencephalogram，EEG）并結合磁共振成像、計算機斷層掃描等影像學檢查掃描大腦結構是癲癇術前診斷的常規方式。近年來，隨著新型腦電檢測技術的發展，記錄腦電信號的方式越來越精準，尤其是顱內腦電圖（intracranial EEG，iEEG）的應用，可以直接測量患者腦部異常放電情況[2]。然而，腦電圖數據量巨大，傳統的人工診斷方式會造成醫療資源的浪費，導致大量癲癇患者得不到有效的治療。通過設計性能優越的檢測算法，對患者的癲癇發作期進行自動檢測，從而協助醫生高效率高質量地對癲癇患者進行診斷，對癲癇的治療有著重要的臨床意義。
基于腦電信號的癲癇發作期自動檢測算法基本步驟是通過對原始腦電信號進行處理，提取相應特征并結合機器學習方法進行分類。小波變換因其具有良好的時頻局部化特性、多分辨率特性、去相關性等特點，在生物醫學信號處理中得到廣泛應用[3]。2000年Kingsbury[4]提出了雙樹復小波變換（dual tree complex wavelet transform，DT CWT）的概念，該小波變換具有近似位移不變性和良好方向選擇性；隨后，Selesnick等[5]提出了雙密度小波變換（double density CWT，DD CWT）的理論，并于2004年[6]提出了雙密度雙樹復小波變換（dual density dual tree CWT，DD-DT CWT）的概念。2014年，Yang等[7]使用DT CWT識別EEG信號中的運動片段，取得了91.07%的準確率。2015年，Roy等[8]使用DD CWT去除多通道EEG信號中的噪聲，結果表明在去除EEG偽影中，DD CWT方法優于獨立成分分析。上述結果證明DT CWT以及DD CWT可以用于腦電信號的識別以及噪聲處理。DD-DT CWT綜合了二者的優點，相比于DT CWT具有更好的重構性、抗混疊性和平移不變性，目前在圖像分割[9]、圖像增強[10]、信號去噪[11]等圖像處理領域和人臉識別領域[12]已有廣泛應用。
1999年，Suykens等[13]在支持向量機（support vector machine，SVM）基礎上改進得到最小二乘支持向量機（least squares SVM，LS-SVM），此方法是基于SVM理論的單隱層多層感知器分類訓練方法。Sharma等[14]使用經驗模態分解法（empirical mode decomposition method，EMD）對腦電信號進行處理，利用LS-SVM作為分類器，取得了84%的準確率和靈敏度。小波熵可以表征信號在不同分解尺度上的能量分布[15]特征。2006年Natwong等[16]從高分辨率心電圖中提取小波熵對心肌梗死患者進行檢測。Deivasigamani等[17]首先使用DT CWT對腦電信號進行分解，選取了多個特征值進行分析，并利用自適應網絡模糊推理系統分類，實驗結果準確率最高可達99%。Xiang等[18]使用模糊熵對癲癇發作和正常腦電圖信號進行檢測和分類，提取不同狀態下的腦電信號特征，并通過SVM進一步訓練和分類，對非發作期和發作期腦電信號分類的準確率為88.55%，靈敏度為90.36%，特異性為87.63%。Bhattacharyya等[19]在2017年提出了一種基于可調Q小波變換（tunable-Q wavelet transform，TQWT）的多變量子帶熵度量方法，并將特征輸入到隨機森林和LS-SVM分類，選擇最優分類，得到最高84.67%的準確率。
綜上所示，我們可以得知國內外對于癲癇自動檢測算法的研究已趨于成熟。但仍存在以下問題：① 腦電信號十分微弱，一般在微伏數量級，易受偽差干擾，尋找穩定合適的信號處理方法至關重要。② 準確的特征可以更好地反映癲癇發作期與非發作期EEG信號的區別，許多研究采用單一特征值，降低了信號的可區分性，如何選擇特征以及如何判斷特征在發作期與非發作期之間是否存在差異是類似研究中的難點。本文針對上述問題，設計并實現了基于機器學習的自動檢測算法，完成了基于患者腦部信號的癲癇發作期與非發作期的二分類。
1 方法
圖1為本文提出的基于DD-DT CWT的癲癇發作期自動檢測算法實施框圖。本方法包括如下四個步驟：① 腦電信號的預處理，將腦電數據處理為同一采樣頻率以及相同的時間片段，并將數據轉化為.mat格式方便后續分析；② 對信號進行DD-DT CWT；③ 從經過DD-DT CWT后的信號中提取小波熵、均值、方差以及能量四類特征；④ 將提取的四類特征送入LS-SVM中進行分類并進行算法性能評估。
圖1
基于DD-DT CWT的癲癇發作期自動檢測算法實施框圖
Figure1.
Principle block diagram of automatic seizure detection algorithm based on dual density dual tree complex wavelet transform
1.1 數據來源與預處理
本文數據來源于公開數據庫UPenn and Mayo Clinic’s Seizure Detection Challenge[20]。該數據庫為美國國家衛生研究所（National Institues of Health，NINDS）設立在 Kaggle 上的競賽數據。原始數據來自8名已接受癲癇手術評估的顳葉和顳葉外癲癇患者，數據集分為兩部分：訓練集（training-set）和測試集（test-set）。測量時使用EEG記錄來自沿海馬前后軸植入的深度電極和不同位置的硬膜下電極網格，癲癇腦電數據分為發作期和非發作期，前者被標注為seizure，后者被標注為non-seizure。發作期以及非發作期數據如表1所示。
表1
發作期與非發作期數據統計
Table1.
Statistics of seizure and non-seizure data
所有8名患者中癲癇發作期數據1 390個，癲癇非發作期數據14 329個，為了方便后續處理，將各個數據按固定時間為一段進行截取，從而得到新的數據集。本研究選擇了常用的數據長度和采樣頻率，綜合考慮后采取了將所有數據采樣頻率降至174 Hz的方案，即每秒數據包含174個數據點，同時按照每段數據1 s進行劃分。發作期和非發作期數據段數均足夠后續的訓練與測試，處理后8名患者中癲癇發作期數據為55 023段，癲癇非發作期數據為501 990段。
1.2 腦電一維DD-DT CWT
從癲癇發作間期到發作期，腦電圖信號的頻率成分通常會發生變化。為了分析這些變化，對腦電信號的分類提供有用的信息，我們對腦電信號采用了DD-DT CWT。DD-DT CWT從離散小波變換的基礎上改進而來，區別于后者，它具有一個尺度函數和兩個小波函數，兩個小波函數之間的偏移量為半個單位。DD-DT CWT通過兩個高通濾波器和一個低通濾波提取腦電信號的高頻成分和低頻成分。本文中，我們定義 為尺度函數，用于分解原始腦電信號的低頻成分，同時對應的低通濾波器為 。定義小波函數 和小波函數 分解原始信號的高頻成分，對應的高通濾波器分別為 和 。由于相鄰小波頻帶間隔小，降低了平移敏感性，因此該小波變換有著更好的近似平移不變性。
表征腦電信號低頻特征對應的尺度方程：
|
表征腦電信號高頻特征對應的小波方程（其中i =1, 2）：
|
經過變換的原始腦電信號分為實部和虛部兩部分，分別由樹 A和樹 B表示。樹A對應第一個雙密度小波變換，使用了兩個高通濾波器 和 來表示腦電信號高頻系數，一個低通濾波器 來表示腦電信號低頻系數。樹B對應第二個雙密度小波變換，使用兩個高通濾波器 和 來表示腦電信號高頻系數，一個低通濾波器 來表示腦電信號低頻系數。記兩個樹對應的尺度函數分別為 和 ，四個小波函數分別為 、、、。
尺度函數和小波函數共同構成了如下的復尺度函數和復小波函數，并協同分解原始腦電信號：
|
|
|
以對原始腦電信號3層分解尺度為例，可以得到如圖2的分解示意圖。X(n)表示原始信號，沿樹A和樹B方向分別進行實部和虛部分解。在分解1層時，我們得到2條低頻子帶 、，以及4條高頻子帶 、、、；當進行第2層分解時，我們將現有的2條低頻子帶繼續分解，產生了新的6條子帶，同時保留了1層分解時產生的4條高頻子帶；在3層分解時，對第2層分解所得低頻子帶進行分解，產生新的6條子帶，保留1層、2層所得高頻子帶。最終獲得了12條腦電高頻子帶和2條腦電低頻子帶分別表示信號中的高頻成分和低頻成分。
圖2
3層分解尺度下雙密度雙樹復小波分解示意圖
Figure2.
Schematic diagram of dual density dual tree complex wavelet decomposition under three-level decomposition scale
1.3 腦電一維DD-DT CWT特征提取
基于傳統信息熵理論，本文提出如下基于原始腦電信號一維DD-DT CWT的小波熵函數。首先對原始腦電信號進行DD-DT CWT變換，第m個分解層數、第n個子帶 對應的小波能量 如式(6)所示：
|
其中N代表當前子帶下小波系數個數，為小波系數。
J為分解層數，則總能量為：
|
相對小波能量定義為：
|
根據Shannon熵的計算公式和定義，此處小波熵為：
|
除了小波熵，本文結合腦電信號的一維DD-DT CWT與傳統特征提取方法，提取的其他三種特征分別為均值（Mean）、方差（SD）和能量（Energy）
均值（Mean）計算公式如式（10）所示：
|
其中，h代表低頻子帶系數，g代表高頻子帶系數，a代表低頻子帶系數個數，b代表高頻子帶系數個數
方差（SD）計算公式如式（11）所示：
|
其中SD低和SD高分別代表低頻子帶方差和高頻子帶方差，n代表方差總數。
能量（Energy）計算公式如式（12）所示：
|
腦電信號經過3層分解尺度下的DD-DT CWT后，可以得到12條高頻子帶和2條低頻子帶，四種類型特征值對應的數目分別為：均值4個，方差4個，能量4個，小波熵4個，共計16個特征值。
1.4 LS-SVM分類器
LS-SVM是基于SVM理論上的單隱層多層感知器分類訓練方法。相比于SVM使用一個超平面將兩組數據分開，LS-SVM回歸采用一個超平面對已知數據進行擬合，該分類器的核心思想是將數據從輸入空間映射到特征空間中，然后在這個特征空間里面構造最優的決策函數，使得最終模型的輸出誤差達到最小化。其最大優勢是運行速度遠高于SVM，但其解不具有稀疏性。
對于訓練集，其中輸入的特征值共包含N維數據，在本文中即為16維特征值，標記分類結果，取值為1（陽性）或 ?1（陰性），分別代表癲癇發作期與非發作期。LS-SVM通過如下公式求解最優：
|
w為系數向量，確定了一個在維度為d的空間中的分類超平面，b為一個常量，γ為懲罰系數， 為分類誤差，求解以上最小化函數可以得到：
|
其中 表示核函數，計算 和 在特征空間中的內積。核函數的主要作用是將非線性可分問題轉換為線性可分問題。本文使用的核函數為式（15）所示的徑向基函數核：
|
2 結果與討論
2.1 評估標準
本文根據文獻[21-24]將腦電信號劃分為1 s片段，最終得到了癲癇發作期數據55 023段以及非發作期數據501 990段，該數據集符合臨床樣本比例。本文采用十折交叉驗證法對算法進行性能評估，將訓練集隨機分成十等分，九份樣本作為訓練集，一份作為測試集，以體現算法的整體性能。本文采用二分類機器學習中腦部疾病檢測常用的評價指標：靈敏度（sensitivity，Sen）、特異性（specificity，Sp）、準確率（accuracy，Acc），分別如式（16）~（18）所示。
|
|
|
其中，TP代表腦電數據真實情況為癲癇發作期，分類器判定為癲癇發作期的數據個數；FP代表腦電數據真實情況為癲癇非發作期，分類器判定為癲癇發作期的數據個數；FN代表腦電數據真實情況為癲癇發作期，分類器判定為癲癇非發作期的數據個數；TN代表腦電數據真實情況為癲癇非發作期，分類器判定為癲癇發作期的數據個數。
2.2 癲癇發作期與非發作期特征差異性分析
本文針對腦電信號微弱、易受干擾以及個體差異大等問題，將DD-DT CWT用于一維腦電信號處理。在分解層數J為3時隨機抽取部分連續數據計算各類特征值，其對比圖如圖3所示，紅色為發作期，藍色為非發作期，橫坐標表示對應的樣本，縱坐標表示樣本對應的特征值。由圖3可見，癲癇患者的發作期與非發作期腦電信號特征提取后所得的均值、方差、能量、小波熵均具有一定的差異性。
圖3
發作期和非發作期各樣本特征值對比
Figure3.
Comparison of feature values between seizure and non-seizure periods
為了進一步對比四類特征中每個特征的分類效果，本文采用t檢驗對特征進行統計學檢驗，如表2所示。首先對八位患者的腦電數據進行DD-DT CWT，分解層數J設定為3層；接著計算各類特征對應的特征值并形成特征矩陣；最后在癲癇發作期與非發作期之間對每個特征進行t檢驗，檢驗水準為0.05。每類特征包含4維特征矩陣，對于均值、方差以及能量三類特征，其中第1維特征代表第3層低頻子帶均值，其余3維特征代表不同分解層數下的高頻子帶均值。分析表2中各特征P值，僅第1層高頻子帶均值以及高頻子帶 小波熵得到的P值略大于0.05，其余特征值得到的P值均小于0.05，即這些特征在發作期與非發作期之間的差異均具有統計學意義。上述結果證明本文選取的特征能夠用于區分癲癇發作期與非發作期之間的腦電信號。本文結合多種特征值完成對癲癇自動檢測算法的設計，為癲癇的特征提取方法提供了新的研究思路。
表2
分解層數J為3時各特征統計學檢驗的P值
Table2.
Statistical test for each feature when decomposition level J is 3
2.3 最優分解層數的探討
由于該小波變換分解層數的可調性，不同分解尺度代表著對信號高頻細節的提取程度不同，也會進一步影響到算法的工作效率和實驗結果。為確定該算法最合適的分解層數，本文選用了相同的訓練集數據和測試集數據在不同分解尺度下進行實驗，比較實驗所得結果。由于分解層數的更換，四種特征值對應的特征值數目會發生改變，如表3所示。
表3
不同分解尺度下特征值數目
Table3.
Number of eigenvalues under different decomposition scales
隨著分解層數的增多，由于分解后所得子帶數目的增加，前三種特征值（均值、方差、能量）的個數每層增加1，而小波熵的特征值數目不隨分解層數的增加而增加。
通過調整算法中分解層數參數（J = 1, 2, 3, 4）運行算法，我們得出了分解層數J分別為1、2、3、4的平均實驗結果以及每位患者分類層數為3時的分類結果，如表4、表5所示。
表4
不同分解層數分類器平均測試結果以及總耗時
Table4.
Average test results and total time consumption of clas sifiers with different decomposition levels
表5
分解層數為3時每位患者測試結果
Table5.
Test results of each patient when the number of decom position layers is 3
從以上各層分解實驗結果可得每位患者的最優分解層數，如表6所示。整體上看，分解層數越多，算法的性能就越好。但由于患者之間的個體差異，不同患者采用該算法的效果有一定差距，為了更完善地討論如何選用最合適的分解層數，我們認為應當從兩方面入手進一步分析。第一是關注不同層數平均準確率；第二是應當重視各個患者準確率達到最大值時算法所處的分解層數。
表6
每位患者的最優分解層數
Table6.
Optimal decomposition levels per patient
基于以上兩點，共有7名患者在選用算法分解層數J為3或4時準確率最大，所以分解層數J應考慮設置為3或4。就整體的平均準確率而言，分解層數J為3或4層時的準確率相差0.43%。從實際情況出發，分解層數每增加一層，特征維度增加3，進而導致算法的運算復雜度增加。通過分析表4可知，隨著分解層數的增加，算法運行時間增加，如果需要檢測的腦電信號數據量龐大（例如很多患者有長達數小時的腦電數據），設計算法時需要優先考慮算法運行過程的參數量以及運行效率。相比于分解層數J為1、2、3層，在4層分解時，分解層數的增大和特征值數目的增加都會在不同程度上增加算法的計算量與運行時間，從而造成運算資源浪費，因此選用3層分解層數更具有實際意義，而從追求算法性能角度來說，4層分解層數下本算法的性能最為優異。
2.4 癲癇發作期自動檢測算法的性能對比
如表7所示，關于癲癇發作期自動檢測的算法主要包括文獻[14]、[19]、[25]、[26]，這些算法均使用了LS-SVM作為分類器，同時選用的特征值均含有熵。本文實驗結果所顯示的準確率、靈敏度和特異性均高于這四項研究，證明了本文提出的基于一維DD-DT CWT小波熵的合理性。值得注意的是，本方法與其他作者所使用的數據集不同，但本方法所用數據庫數據量（發作期55 023段，非發作期501 990段，每段1 s）大于其他作者的數據庫數據量（發作期3 750段，非發作期3 750段，每段20 s）。同時，大部分患者一次癲癇發作持續時間小于兩分鐘，數據段劃分過長可能會導致本算法在應用于臨床數據時檢測出的發作期存在部分非發作期片段，因此本文綜合考慮了其他研究[21-23]的時間劃分方式，最終決定選擇將數據劃分為1 s片段。另一方面，癲癇患者腦電信號中非發作期時長遠大于發作期時長[24]，數據集符合真實樣本比例。本文將該小波變換用于對一維腦電信號的處理，實驗結果的成功證明該方法對于腦電信號分析的可行性。
表7
本文方法與近期國際報道的算法性能對比表
Table7.
Comparison of the performance of this method and the algorithm reported recently
盡管本文基于DD-DT CWT實現了對腦電癲癇發作期信號的自動識別，然而本文算法的結果是在樣本比例接近臨床情況下不區分患者的結果，主要適用針對某個患者長期監測的場景。隨著后續數據采集過程中患者數量的增加，還有待探索泛化性能更佳的算法，使得算法性能提升的同時滿足不同患者的檢測需求，并減少因患者個體差異而造成的分類結果不穩定問題。
3 結論
本文提出了一種基于DD-DT CWT的顱內腦癲癇發作自動檢測算法，準確地實現了對癲癇發作期的自動識別與分類。本方法利用DD-DT CWT對經過預處理的信號進行分析，并從中提取了四類特征（小波熵、方差、能量以及均值），t檢驗結果表明在癲癇發作期與非發作期之間四類特征存在差異，采用LS-SVM對特征進行分類，可完成對癲癇發作期的準確識別。相較于已報道的癲癇發作自動檢測算法，本文方法在保證較高檢測精度的基礎上，綜合應用了多種腦電特征。本文工作可輔助臨床醫生劃定癲癇發作期與非發作期，具有潛在的臨床應用價值。
利益沖突聲明：本文全體作者均聲明不存在利益沖突。
