一種針對EEG 偽跡的多窗卷積壓制算法

孫望強，劉亮，郭慶功

（四川大學電子信息學院，成都610065）

0 引言

腦電信號的各種波形傳達了具有臨床價值的信息，對腦電信號的分析在神經科學、認知科學等領域都得到了廣泛的應用[1-2]。測量腦電信號的過程中，腦電信號會受到電噪聲等外部干擾和眼電、心電等生理偽跡信號的干擾[3-5]。眼電等偽跡干擾幅值會比腦電信號大10 倍之多，很大程度上會“淹沒”神經信號[6-7]，因此腦電信號的信號SAR 很低，會使得對腦電信號的分析變得困難。

Nolan H.等人[8]基于統計特征對神經信號和偽跡干擾成分進行區分并有效刪除偽跡成分，Petr Nejedly、Radüntz T 等人[9-10]將機器學習算法應用于腦電特征中實現偽跡的自動識別消除，Soomro MH[11]基于經驗模式分解結合ICA 的算法用于眼電偽跡消除，并對比了偽跡前后腦電信號的SAR 和相關系數。但這些研究都未在少通道（如8 導）采集系統下進行分析，且都是通過粗刪除的方式對判斷為偽跡的成分進行處理。

在實際應用中，考慮腦電采集設備的成本問題，會選用8 導腦電信號采集設備。多通道采集系統測量的數據量較多，有效地區分腦電偽跡和神經信號之后刪除偽成分可能不會導致神經信號泄露的問題，但在少通道采集系統（8 道）下，即使很好地判斷出“偽跡”成分，這些成分中勢必會含有一部分神經信號成分，如果通過粗刪除的方式去除，會導致神經信號泄露的問題[12]。

為了解決這個問題，本文在Delpozo-Ba 等人[13]研究基礎上，結合信號時、空特征，在局部化的時間片段上進行偽跡特征判別，對分類出的含有偽跡干擾較多的整導成分不進行粗刪除處理，而是利用不同的窗函數優化組合并卷積壓制偽跡特征，結果證明該算法可以消除眼電等偽跡且降低神經信號泄露的影響，提高處理之后腦電信號的SAR。

1 偽跡消除算法原理

1.1 ICA在偽跡消除算法中的應用

腦電信號分析中，最廣泛使用的模型是假設被測的腦電信號是多個腦源與偽跡產生的電信號的線性組合，并瞬間傳播到頭皮區域[14-16]。可以用數學方程式（1）表示：

X=As+V （1）

其中X=[x(1 ),x( 2 ),…,x( N )]是一個n 行N 列的腦電信號數據矩陣，行n 表示測量的腦電信號的通道數（采集設備導數），列N 表示采樣的數據點數；A 是一個n×n 的未知混合矩陣，V 代表外部噪聲。ICA 的核心目標是估計解混矩陣W，得到觀察信號的源成分，模型方程式如公式（2）所示：

S=WX （2）

其中S 是一個n×N 的未知源成分矩陣，由神經信號成分與偽跡成分組成。為了確?；镜腎CA 模型能被估計，要做一些假設和約束：第一，源成分S 被假設是統計獨立的；第二，源成分具有非高斯的分布；第三，假定混合矩陣A 是方陣，且可逆；第四，假定所有的混合變量X 與獨立成分S 都是零均值；第五，觀察信號X 必須可以白化。基于這些約束，利用ICA 解混觀察信號X 得到源成分，并通過偽跡特定的統計特征、時間特征、空間特征等可將源成分S 分類成偽跡成分和神經信號成分。

1.2 實現代碼多窗卷積組合優化的局部偽跡壓制算法原理

步驟一，通過ADJUST[17]算法分離觀察信號，并根據偽跡時、空特征分類得到干凈成分。步驟二，將觀察信號通過ICA 處理之后，在局部化時間片段上根據突發性電壓定位偽跡特征，并利用窗函數，結合步驟一得到的干凈成分插值壓制局部偽跡。一方面在更精細的時間片段上處理，降低了偽跡的影響，另一方面保留了大量的有用的神經信號。流程如圖1 所示。

眨眼等偽跡可以定位到相應的每一個時刻，表現為測量信號波形的突發性，可以用局部時刻突發性的電壓來定位這些特征時刻。為了量化特征，首先通過ICA 分離觀察信號得到的成分數據并對其求樣本Z 分數，Z 分數是一種對原始數據的均值和方差進行標準化的常用統計數據分析方法。其定義如公式（3）所示：

這里X 表示數據，E(X)表示其期望值，σ(x)表示其標準差，其中X 表示ICA 分離之后的成分的每一個電壓值和一階導數值。這樣就可以利用Z 分數來量化特征，表示為A[c,k]，標識第k 時刻成分c 的特征。

圖1 算法流程圖

由于一些小的脈沖干擾噪聲可能會造成觀察信號異常的數據點，可以用窗函數卷積運算來穩定所提取到的特征，更穩定的特征可以表示為B[c,k]，其定義如公式（4）：

這里A[c,k]表示由Z 分數量化的特征，*表示卷積運算，M 代表用于卷積的窗函數，L 表示整個的卷積范圍，公式（4）的分母是一個歸一化因子，它有效地在每一時刻將積分轉換為加權平均，這樣做可以抵消掉卷積的邊緣效應。

公式（4）得到穩定特征之后，需要創建一個向量來標記第n 個時刻的成分c 的特征是不是屬于偽跡。對每個時刻的某個成分的每一個特征設置一個閾值，設置閾值的目的是為了分類含有局部偽跡特征的成分和神經信號成分。為了更加方便與簡潔，采用二分類閾值為p。其決策邏輯函數如公式（5）所示：

式（5）中自變量x 就是需要判斷的特征，也就是由式（5）得到的穩定特征，這里需要傳入的是max([B[c ,k-d],…,B[c,k+d]])，其中d 表示每個時間點附近需要的一個緩沖區域，使得這一區域內的特征被上面的邏輯函數f(x)分類為噪聲。緩沖區域取0.09s 是驗證為比較合適的，在這個0.09s 的區域會產生連續的噪聲，可以與相鄰的時刻很好的分離開。對每個成分每個時刻的特征B[c,k]用邏輯判斷之后，用F[c,k]表示，F[c,k] 的值用1 和0 分別表示這個特征是否是偽跡。

經過式（5）處理之后，需要對那些被標注為1 的特征時間段進行處理，不直接刪除這些時間片段，而是通過步驟一獲得的干凈成分按比例插值處理，也就是按照邏輯函數F[c,k]來混合兩種成分。其中控制混合兩種成分的信號為Y[c,k]，定義如公式（6）所示：

其中M1 是長度為N 的窗函數，用于對檢測區域的邊緣進行調整，從而平滑混合中的過渡，最后得到干凈的成分R[c,k]。

Y[ c,k] 是由式（7）得到的控制信號，G[ c,k] 是基于時、空特征得到的相對干凈的成分，D[ c,k] 是由觀察信號盲源分離得到的成分，最后混合G[ c,k] 和D[c ,k] 得到保留更多有用信號的干凈成分R[c ,k] ，這樣得到的結果不僅減少了神經信號的泄露，而且成功地壓制了偽跡。

2 實驗結果分析

本文腦電信號采集系統是ANT-NEURO EEG/ERP系統，給被試人員分別采用8 電極、32 電極采集被試腦電信號。其中A/D 采樣頻率設置為1000Hz，頭皮阻抗小于5KΩ，受試者在頭腦清醒的情況下接受“怪球”實驗。對采集的腦電數據利用EEGLAB、MNE-Python 聯合進行算法編程與分析，8 電極、32 電極圖分別如圖2、圖3 所示。

圖2 8道電極圖

圖3 32 道電極圖

其中FPz、Fz、Cz、Pz、Oz 等分別表示腦的前額區域、額區、中央區、頂區和枕區，8 電極采集系統采用“平均”參考電極，32 道采集系統采用M1、M2 為參考電極，在分析時以刺激呈現200ms 為基線，每一個epcoh分析時程是-200ms 到800ms。

首先將采集信號用FIR 濾波器進行1hz-30Hz 的濾波，隨后進行基線校準和分段處理，可以得到圖4（8電極），圖5（32 電極）的觀察信號波形。

圖4 經過分段處理后的8道觀察信號

圖5 經過分段處理后的32道觀察信號

圖4、圖5 中可以很明顯地看到眼電等偽跡干擾，需要壓制這些偽跡造成的干擾。首先是要獲得觀察信號的源成分，因此我們采用ICA 算法對圖4、圖5 中的觀察信號進行盲源分離，分別得到8 導、32 導時間波形成分圖6、圖7 所示。

圖6 8電極源成分時間波形圖

圖7 32電極源成分時間波形圖

利用ADJUST[17]算法，基于時空特征，可以得到相對干凈的成分，但會造成神經信號泄露。當得到相對干凈的成分之后，結合最優的多窗卷積壓制算法，將相對干凈的成分按比例局部插值到偽跡成分以達到壓制偽跡的目的，這樣做不僅能得到更加干凈的神經信號成分，而且能保留更多的神經信號，實現偽跡噪聲降低和神經信號泄漏之間的均衡。利用多窗卷積壓制算法得到的腦源成分（8、32）如圖8、9 所示，對比圖6、圖7發現大部分偽跡成分已經被成功壓制。

圖8 處理之后8電極干凈源成分

圖9 處理之后32電極干凈源成分

最后利用混合矩陣可恢復觀察信號，如圖10、圖11 所示。

圖10 恢復的8電極觀察信號

圖11 恢復的32電極觀察信號

對比圖4、圖5 得到的觀察信號，可以看到圖10、圖11 中眨眼等偽跡干擾已經成功被壓制。由以上分析可得，利用本文算法處理之后，所采集的腦電信號的偽跡干擾確實成功被壓制。為了進一步量化算法去噪效果，采用SAR（噪聲偽跡比）衡量，其定義如公式（8）所示。

其中mean 代表求數據的均值，X 代表原始未校正的EEG 信號。SAR 體現了消除偽跡的效果，8 電極、32電極SNR 結果分別如表1、表2 所示。

表1 8 電極采集系統SAR 對比表

表2 32 電極采集系統SAR 對比表

3 結語

在本文中，基于腦電信號時、空特征和窗函數卷積組合優化的方式來對局部偽跡進行壓制。在8、32 電極采集系統下，未處理腦電信號SAR 分別是-19.765dB，-19.616dB，通過本文提出的方法處理后腦電信號的SAR 分別是2.832dB，2.743dB，且在8 電極采集系統下效果較FastICA、EMD-ICA 處理后的SAR 有明顯提高。實驗結果表明該算法不僅能夠壓制眨眼偽跡干擾，還能較好地保留神經信號的成分，實現了偽跡噪聲降低和神經信號泄露之間的均衡。