王永軒 邱天爽 劉 蓉

（大連理工大學(xué) 電子信息與電氣工程學(xué)部，大連 116024）

引言

腦電誘發(fā)電位（evoked potential，EP）包含了有關(guān)神經(jīng)系統(tǒng)傳導(dǎo)通路的信息，因此有著重要的臨床應(yīng)用價值。由頭皮表面測得的EP信號疊加了自發(fā)腦電（electroencephalogram，EEG）等背景噪聲［1］，信噪頻譜重疊且信噪比非常低，提取難度較大。因為EP信號對外界刺激具有鎖時性，所以疊加平均法（ensemble average，EA）和加權(quán)疊加平均法被廣泛應(yīng)用在臨床EP提取中［1］。該方法簡單且易于實現(xiàn)，但需要幾百次甚至上千次刺激才能得到較為理想的EP信號。為了實現(xiàn)快速提取、不丟失細(xì)節(jié)以及不受EP信號的時變影響，少次提取成為一個重要的課題。此外，誘發(fā)電位是由特定位置的大腦皮層細(xì)胞產(chǎn)生的反應(yīng)，傳導(dǎo)到頭皮不同位置時相當(dāng)于經(jīng)過了不同的系統(tǒng)，因此誘發(fā)電位的單導(dǎo)提取更顯得尤為重要。

在過去的幾十年里EP信號少次提取的研究取得了長足的進(jìn)步［2-4］，包括PSR算法［5］、自適應(yīng)濾波法［6］、高階累積量法［7］、神經(jīng)網(wǎng)絡(luò)法［8-9］、獨立分量分析［10-13］、小波分析［14-17］等等。PSR算法通過FFT將時域信號變換到頻域再進(jìn)行處理，最終設(shè)計出最佳濾波器對觀測信號進(jìn)行濾波。自適應(yīng)濾波法使用自適應(yīng)算法來改變?yōu)V波器的參數(shù)和結(jié)構(gòu)，這樣可以跟蹤觀測信號的時變特征。高階累積量法是利用高階累積量對高斯噪聲及部分有色噪聲不敏感的特性，設(shè)計出一個有限沖激響應(yīng)匹配濾波器，從而達(dá)到增強有用信號同時抑制噪聲的作用。神經(jīng)網(wǎng)絡(luò)是非線性處理過程，能夠適應(yīng)誘發(fā)電位的非線性特性。獨立分量分析是伴隨著解決盲源分離問題發(fā)展起來的多導(dǎo)信號處理方法，它是以信號之間的獨立性為前提，可以分離疊加在一起的獨立信號。小波分析是一種時頻分析方法，可以同時反映出信號的全局特征和局部特征，這樣可以根據(jù)誘發(fā)電位的性質(zhì)對小波系數(shù)進(jìn)行選取。針對誘發(fā)電位與自發(fā)腦電的特性，本研究提出了使用預(yù)白化與小波分析相結(jié)合的方法來實現(xiàn)對腦電誘發(fā)電位的單導(dǎo)少次提取。

1 自發(fā)腦電建模與預(yù)白化

在頭皮表面獲得的觀測信號除了誘發(fā)電位和自發(fā)腦電以外，還包括工頻、心電、肌電、眼電等干擾以及各種噪聲。由于目前對腦電信號的去噪方法已經(jīng)比較完善，所以可以假設(shè)非腦電的信號已被除去，觀測信號只是誘發(fā)電位與自發(fā)腦電的疊加。因此觀測信號的模型可表示為

式中，x（n）、s（n）和v（n）分別表示觀測到的帶噪信號、純凈EP信號和自發(fā)腦電，n為采樣點序號，N為單次EP信號長度。在給定的采樣頻率下，自發(fā)腦電表現(xiàn)為有色噪聲，用AR模型表示為

式中，P為濾波器的階數(shù)，w（n）為激勵白噪聲。濾波器的系統(tǒng)傳遞函數(shù)為

在實際處理中為了得到自發(fā)腦電的AR模型，需要一段在刺激前記錄的觀測信號，即純凈的自發(fā)腦電，令其為v0（m），（m=0，1，…，M-1），這樣可以由v0（m）得到自發(fā)腦電AR模型的參數(shù)估計。自發(fā)腦電在采樣頻率為fs=1kHz時典型的AR模型為［18］

為了應(yīng)用小波分析方法去噪，首先要對觀測信號進(jìn)行濾波處理，使其中的自發(fā)腦電白化為白噪聲。將觀測信號經(jīng)過一個系統(tǒng)傳遞函數(shù)為G1（z）=1/G0（z）=A（z）的濾波器，則其中的自發(fā)腦電將被白化為白噪聲，對應(yīng)的差分方程為

經(jīng)過了白化濾波器后，誘發(fā)電位成為期望信號u（n）

同時觀測信號成為y（n），并有

由于s（n）與v（n）是相互獨立的，所以u（n）與w（n）也是獨立的。同時由于v（n）是零均值的平穩(wěn)有色噪聲，所以w（n）是零均值的平穩(wěn)白噪聲。經(jīng)過白化后就可以通過小波分析方法從y（n）中估計出u（n），最后由u（n）恢復(fù)s（n）。

2 小波分析與去噪

2.1 小波分析原理

小波分析是一種時頻分析方法，該方法需首先確定一個母小波，再通過母小波的伸縮平移族得到一組基。作為一種分解分離方法，小波分析包括分解、處理和綜合三個過程。分解就是將原始信號向這組基上進(jìn)行投影，得到不同尺度和平移參量下的系數(shù)，即小波系數(shù)。然后根據(jù)需要對小波系數(shù)進(jìn)行加工處理。將處理后的小波系數(shù)進(jìn)行反變換即為綜合過程。小波分析的另一重要概念是多分辨率分析，可使用Mallat快速算法逐級計算。一個確定的小波分析相當(dāng)于一個雙通道濾波器組，包括低通濾波器H0（z）和高通濾波器H1（z），在不同的分辨率級下濾波器組是相同的。分解過程中的第j+1級的尺度系數(shù)cj+1（k）和小波系數(shù)dj+1（k），可由第j級的尺度系數(shù)cj（k）按式（8）得到，逐級遞推計算可得到所有的小波系數(shù)。

綜合過程中第j級的尺度系數(shù)cj（k），可由第j+1級的尺度系數(shù)cj+1（k）和小波系數(shù)dj+1（k）按式（9）得到，直到恢復(fù)原始信號。

以上分解和綜合中使用的濾波器組是相同的，因此這是正交小波分析。在正交小波分析中，從原始信號空間到小波系數(shù)空間的投影矩陣是正交陣。將序列y（n）、u（n）和w（n）表示為向y、u和w，再分別將它們在各個分辨率級上的所有小波系數(shù)組合在一起，用向量Y、U和W表示，則有Y=TTy、U=TTu和W=TTw，其中T為正交陣。以Harr正交小波分析為例說明，在Harr正交小波分析中有h0（n）=｛1/，1/｝，h1（n）=｛1/，-1/｝。設(shè)原始信號為y=［y（0），y（1），...y（7）］T，經(jīng)過Mallat快速算法得到的第1、2、3級的小波系數(shù)和第3級的尺度系數(shù)組成的向量為

Y=［d1（0），d1（1），d1（2），d1（3），d2（0），d2（1），d3（0），c3（0）］T

則T為

所有的正交小波分析都存在這樣的正交陣T對應(yīng)這種正交變換關(guān)系。對于白噪聲w，有E（wwT）=σ2I，其中σ為白噪聲功率的均方根，所以E（WWT）=E（TTwwTT）=σ2I。這說明白噪聲經(jīng)正交小波變換后得到的小波系數(shù)在各個分辨率級上仍表現(xiàn)為白噪聲，即白噪聲有限的能量將均勻分散在所有的小波系數(shù)上。而對于期望信號u，由于并未被白化而有一定的規(guī)則性，以及小波分析的聚焦作用，小波系數(shù)非零值將集中在少數(shù)系數(shù)上，大多數(shù)系數(shù)接近于零。這樣通過對小波系數(shù)進(jìn)行處理后再進(jìn)行反變換，則可以在很大程度上降低噪聲能量，提取出期望信號。

2.2 小波去噪

為實現(xiàn)小波去噪，應(yīng)設(shè)計一個非線性函數(shù)對Y進(jìn)行處理，使非線性函數(shù)的輸出盡可能接近U，通用的方法主要是閾值法［8］，但閾值的選取對結(jié)果影響較大，難以確定合適的取值。本研究采用加權(quán)結(jié)合閾值的方法，目標(biāo)是使重構(gòu)信號與期望信號u的均方誤差達(dá)到最小。令加權(quán)系數(shù)向量為b=［b（0），b（1），..，b（N-1）］T，則與u的均方誤差為

令?F（b）/?b=0，可得

因為E（Y2（n））=E（（U（n）+W（n））2）=U2（n）+σ2，所以式（12）可近似為

再定義閾值為σ，則有

在實際計算時正交陣T是無需知道的，通過式（8）即可由y（n）得到Y(jié)（n），然后由式（13）、（14）得到加權(quán)系數(shù)b（n），相乘可得U（n）的估計值（n）=b（n）·Y（n），最后由式（9）進(jìn)行反變換得到u（n）的估計值（n）。在自發(fā)腦電建模與預(yù)白化部分中，已經(jīng)提到需要一段在刺激前記錄的觀測信號，即不包含誘發(fā)電位成分的純凈的自發(fā)腦電，由這段觀測信號被白化后得到的白噪聲，就可以估計出上述白噪聲功率均方根σ。

3 恢復(fù)EP信號

利用小波分析方法從y（n）中估計出u（n）后需要由u（n）恢復(fù)出s（n），可以采用兩種方法，一是頻域法，二是矩陣法，理論和仿真實驗證明兩種方法的結(jié)果幾乎完全一樣，只是在起點處有微小的差別。頻域法是首先將式（6）中a（n）=｛1，a1，a2，…，aP｝補零延拓成N點，然后對a（n）和（n）進(jìn)行FFT變換得到A（k）和（k），最后對（k）=（k）/A（k）進(jìn)行IFFT得到（n），則（n）就是對s（n）的估計結(jié)果。矩陣法是將式（6）由卷積形式改寫為矩陣形式u=As，其中系數(shù)矩陣A為

顯然矩陣A是下三角矩陣，所以一定可逆，令=A-1，則（n）就是對s（n）的估計結(jié)果。在實際計算時不用求逆，對于系數(shù)為下三角矩陣的矩陣方程可以用迭代方法求解。

4 仿真實驗與分析

仿真實驗使用的EP信號為經(jīng)過疊加平均法得到貓的體感誘發(fā)電位，有效頻帶為0～20Hz，峰值頻率為10Hz。自發(fā)腦電為由式（4）的AR模型產(chǎn)生的有色噪聲，有效頻帶為0～30Hz，峰值頻率為13Hz。實驗中采樣頻率為1kHz，采樣時間為512ms，有色噪聲符合自發(fā)腦電的特性，并且在頻帶上完全覆蓋了EP信號，初始信噪比（SNRi）由0dB逐漸降低到-10dB。算法使用的小波為db3正交小波。

圖1和圖2分別為SNRi=0dB和SNRi=-10dB條件下小波方法和20次疊加平均法的估計結(jié)果。圖3給出了在0dB條件下有色噪聲和EP信號的原始信號、預(yù)白化后信號和小波系數(shù)。從圖3中可以明顯看出有色噪聲被白化為白噪聲，其小波系數(shù)也表現(xiàn)為白噪聲特性。而EP信號經(jīng)白化濾波器后得到的期望信號仍有一定規(guī)則性，其小波系數(shù)非零值將集中在少數(shù)系數(shù)上，大多數(shù)系數(shù)接近于零。圖4和圖5分別為隨著SNRi的降低，不同方法估計結(jié)果的信噪比（SNRo）和估計結(jié)果與純凈EP信號相關(guān)系數(shù)（r）的變化情況（數(shù)據(jù)為50次實驗結(jié)果的平均值）。

小波分析的優(yōu)勢在于它是時頻分析，這樣可以將期望信號變換到少數(shù)系數(shù)上，實現(xiàn)了能量集中。在本算法中預(yù)白化是必不可少的，只有在確保噪聲是白噪聲的前提下小波分析才能發(fā)揮它的優(yōu)勢，因為白噪聲的小波系數(shù)分散在所有值上，這樣才能突出期望信號的少數(shù)系數(shù)。同時由實驗結(jié)果也可以看出估計的效果受制于信噪比，這是因為雖然白噪聲小波系數(shù)均勻分散在所有值上，但當(dāng)噪聲足夠強時，仍可以淹沒期望信號的少數(shù)較大系數(shù)，從而使期望信號系數(shù)估計的誤差增大。通過整個估計過程可知，對觀測信號小波系數(shù)的非線性處理是至關(guān)重要的環(huán)節(jié)，如果已知更多的先驗信息就可以更靈活地處理小波系數(shù)，從而獲得更好的效果。同時也應(yīng)注意到正交小波的選取也很重要，如果選取的小波在某一分辨率級上與期望信號非常接近，則期望信號的小波系數(shù)會更加集中，這樣重構(gòu)信號與期望信號的均方誤差會更小。

圖1 SNR=0dB條件下的估計結(jié)果。（a）有色噪聲；（b）純凈EP信號；（c）單次帶噪信號；（d）小波方法估計結(jié)果；（e）20次帶噪信號；（f）20次疊加平均結(jié)果Fig.1 The estimated result at SNRi=0dB.（a）colored noise；（b）clean EP signal；（c）one-sweep noisy signal；（d）estimated result with wavelet method；（e）twenty-sweep noisy signal；（f）estimated result with twenty-EA method

PSR算法不具備時頻分析與聚焦信號能量的功能，所以在系數(shù)處理和信號重建中只能消除有限的噪聲，在信號與噪聲頻譜重疊的情況下算法性能會退化。自適應(yīng)濾波法通常是以疊加平均法得到的波形做為參考信號，其性能取決于參考信號和EP的相關(guān)程度。獨立分量分析是少次提取EP信號的新方法之一，但只能處理多導(dǎo)腦電信號，同時還要做很多假設(shè)，包括信號源必須是獨立的，信號源數(shù)目要不多于觀測信號路數(shù)，觀測信號是信號源的線性混合，混和矩陣是列滿秩的等等。目前用于EP信號提取的小波分析方法通常需要已知EP信號的部分先驗信息，并且需要人為選擇小波系數(shù)，這在應(yīng)用中是不盡可行的。本算法雖不需要EP信號的先驗信息，但需要假設(shè)觀測信號只是誘發(fā)電位與自發(fā)腦電的疊加，并沒有考慮其他高斯/非高斯、周期/非周期性的各種不同噪聲。仿真實驗表明本算法效果與20次疊加平均效果相當(dāng)。

圖2 SNR=-10dB條件下的估計結(jié)果。（a）有色噪聲；（b）純凈EP信號；（c）單次帶噪信號；（d）小波方法估計結(jié)果；（e）20次帶噪信號；（f）20次疊加平均結(jié)果Fig.2 The estimated result at SNRi=-10dB.（a）colored noise；（b）clean EP signal；（c）single sweep noisy signal；（d）estimated result with wavelet method；（e）twenty-sweep noisy signal；（f）estimated result with twenty-EA method

圖3 SNR=0dB條件下有色噪聲和EP信號的原始信號、白化后信號和小波系數(shù)。（a）有色噪聲；（b）EP信號；（c）白化后得到的白噪聲；（d）期望信號；（e）白噪聲小波系數(shù)；（f）期望信號小波系數(shù)Fig.3 The original signals，whitened signals and wavelet coefficients of the colored noise and EP signal at SNRi=0dB.（a）colored noise；（b）clean EP signal；（c）whitened noise；（d）expected signal；（e）the coefficients of whitened noise；（f）the coefficients of expected signal

圖4 隨著SNRi的降低不同方法估計結(jié)果的信噪比（SNRo）的變化情況Fig.4 The relationship between SNRowith SNRi by different method

圖5 隨著SNRi的降低不同方法估計結(jié)果的相關(guān)系數(shù)（r）的變化情況Fig.5 The relationship between correlation coefficient r with SNRiby different method

5 結(jié)論

本研究提出了一種基于小波分析的方法用于腦電誘發(fā)電位單導(dǎo)少次提取。該算法首先將有色噪聲進(jìn)行預(yù)白化，確保噪聲的小波系數(shù)分散在所有系數(shù)上，而期望信號的小波系數(shù)非零值集中在少數(shù)系數(shù)上。這樣通過對小波系數(shù)進(jìn)行處理后再進(jìn)行反變換，最后恢復(fù)誘發(fā)電位，可以在一定程度上降低噪聲能量。實驗表明，采用加權(quán)結(jié)合閾值的小波分析方法，可以抑制自發(fā)腦電的干擾，與初始信噪聲比相比，估計結(jié)果的信噪比獲得了提高。更重要的是本算法擺脫了對EP信號先驗信息的依賴性，無需人為選擇小波系數(shù)，符合實際應(yīng)用背景。理論分析和仿真試驗均表明了該算法的應(yīng)用于EP信號單導(dǎo)少次提取問題有效性。將在下一步的研究工作中繼續(xù)完善該算法，并將其應(yīng)用于實際測試EP信號提取問題中。

［1］潘映輔.臨床誘發(fā)電位學(xué)（第二版）［M］.北京：人民衛(wèi)生出版社，2000.

［2］朱常芳，胡廣書.誘發(fā)電位快速提取算法的新進(jìn)展［J］.國外醫(yī)學(xué)生物醫(yī)學(xué)工程分冊，2000，23（4）：211-216.

［3］張金鳳，邱天爽.誘發(fā)電位波形提取方法及進(jìn)展［J］.北京生物醫(yī)學(xué)工程，2003，22（4）：3030-307.

［4］陶彩林，鄒凌，馬正華，等.誘發(fā)電位單次提取技術(shù)的研究進(jìn)展［J］.生物醫(yī)學(xué)工程學(xué)雜志，2009，26（5）：1158-1161.

［5］Kavch M.A new method for the estimation of average evoked responses［J］.IEEE Trans on System，Man and Cybernetics，1978，8（5）：414-417.

［6］Stefanos DG，Perttu OR，Mika PT，et al.Single-trial dynamical estimation of event-related potentials：A kalman filter-based approach［J］.IEEE Transactions on Biomedical Engineering，2005，52（8）：1397-1406.

［7］Gharieb RR，Cichocki A.Noise reduction in brain evoked potentialsbased on third-order correlations［J］.IEEE Transactions on Biomedical Engineering，2001，48（5）：501-512.

［8］Qiu Wei，Kenneth SM，F(xiàn)rancis HY，et al.Adaptive filtering of evoked potentials with radial basis function neutral network prefilter［J］.IEEE Transactions on Biomedical Engineering，2002，49（3）：225-231.

［9］Qiu Wei，Chang Chunqi，Liu Wenqing，et al.Real-time datareusing adaptive learning of a radial basis function network for tracking evoked potentials［J］.IEEE Transactions on Biomedical Engineering，2006，53（2）：226-227.

［10］Aapo H，Erkki O.Independent component analysis：algorithms and applications［J］.Neural Networks，2000，13（4-5）：411-430.

［11］洪波，唐慶玉，楊福生.ICA在視覺誘發(fā)電位的少次提取與波形分析中的應(yīng)用［J］.中國生物醫(yī)學(xué)工程學(xué)報，2000，19（3）：334-341.

［12］畢曉輝，邱天爽，朱勇，等.基于經(jīng)驗?zāi)Ｊ椒纸夂酮毩⒎至糠治龅膯蔚郎俅蜤P信號提?。跩］.中國生物醫(yī)學(xué)工程學(xué)報，2008，27（6）：817-821.

［13］師黎，鐘麗輝，王瑞.基于虛擬通道ICA-WT大鼠視覺誘發(fā)電位少次提?。跩］.中國生物醫(yī)學(xué)工程學(xué)報，2010，29（3）：379-383.

［14］Bartnik EA，Blinowska KJ，Durka PJ.Single evoked potential reconstruction by means of wavelet transform［J］.Biological Cybemetics，1992，67：175-181.

［15］Quiroga R，Garcia H.Single-trial event-related potentials with wavelet denoising［J］.Clin Neurophysiol，2003，114：376-390.

［16］Wang Zhisong，Alexander M，David A，et al.Single-trial evoked potential estimation using wavelets［J］.Computers in Biology and Medicine，2007，37（4）：463-473.

［17］鄒凌，陶彩林，王正洪.基于小波變換的誘發(fā)電位信號去噪研究［J］.計算機(jī)應(yīng)用與軟件，2010，27（1）：85-87.

［18］Yu Xiaohu，He Zhenya，Zhang Yisheng.Time-varying adaptive filters for evoked potential estimation［J］.IEEE Transactions on Biomedical Engineering，1994，41（11）：1062-1071.