王永軒 邱天爽 劉 蓉

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

引言

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

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

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

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

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

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

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

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

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

同時(shí)觀測(cè)信號(hào)成為y（n），并有

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

2 小波分析與去噪

2.1 小波分析原理

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

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

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

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

則T為

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

2.2 小波去噪

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

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

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

再定義閾值為σ，則有

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

3 恢復(fù)EP信號(hào)

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

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

4 仿真實(shí)驗(yàn)與分析

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

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

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

圖1 SNR=0dB條件下的估計(jì)結(jié)果。（a）有色噪聲；（b）純凈EP信號(hào)；（c）單次帶噪信號(hào)；（d）小波方法估計(jì)結(jié)果；（e）20次帶噪信號(hào)；（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í)頻分析與聚焦信號(hào)能量的功能，所以在系數(shù)處理和信號(hào)重建中只能消除有限的噪聲，在信號(hào)與噪聲頻譜重疊的情況下算法性能會(huì)退化。自適應(yīng)濾波法通常是以疊加平均法得到的波形做為參考信號(hào)，其性能取決于參考信號(hào)和EP的相關(guān)程度。獨(dú)立分量分析是少次提取EP信號(hào)的新方法之一，但只能處理多導(dǎo)腦電信號(hào)，同時(shí)還要做很多假設(shè)，包括信號(hào)源必須是獨(dú)立的，信號(hào)源數(shù)目要不多于觀測(cè)信號(hào)路數(shù)，觀測(cè)信號(hào)是信號(hào)源的線性混合，混和矩陣是列滿秩的等等。目前用于EP信號(hào)提取的小波分析方法通常需要已知EP信號(hào)的部分先驗(yàn)信息，并且需要人為選擇小波系數(shù)，這在應(yīng)用中是不盡可行的。本算法雖不需要EP信號(hào)的先驗(yàn)信息，但需要假設(shè)觀測(cè)信號(hào)只是誘發(fā)電位與自發(fā)腦電的疊加，并沒(méi)有考慮其他高斯/非高斯、周期/非周期性的各種不同噪聲。仿真實(shí)驗(yàn)表明本算法效果與20次疊加平均效果相當(dāng)。

圖2 SNR=-10dB條件下的估計(jì)結(jié)果。（a）有色噪聲；（b）純凈EP信號(hào)；（c）單次帶噪信號(hào)；（d）小波方法估計(jì)結(jié)果；（e）20次帶噪信號(hào)；（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信號(hào)的原始信號(hào)、白化后信號(hào)和小波系數(shù)。（a）有色噪聲；（b）EP信號(hào)；（c）白化后得到的白噪聲；（d）期望信號(hào)；（e）白噪聲小波系數(shù)；（f）期望信號(hào)小波系數(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的降低不同方法估計(jì)結(jié)果的信噪比（SNRo）的變化情況Fig.4 The relationship between SNRowith SNRi by different method

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

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

［2］朱常芳，胡廣書(shū).誘發(fā)電位快速提取算法的新進(jìn)展［J］.國(guó)外醫(yī)學(xué)生物醫(yī)學(xué)工程分冊(cè)，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在視覺(jué)誘發(fā)電位的少次提取與波形分析中的應(yīng)用［J］.中國(guó)生物醫(yī)學(xué)工程學(xué)報(bào)，2000，19（3）：334-341.

［12］畢曉輝，邱天爽，朱勇，等.基于經(jīng)驗(yàn)?zāi)Ｊ椒纸夂酮?dú)立分量分析的單道少次EP信號(hào)提?。跩］.中國(guó)生物醫(yī)學(xué)工程學(xué)報(bào)，2008，27（6）：817-821.

［13］師黎，鐘麗輝，王瑞.基于虛擬通道ICA-WT大鼠視覺(jué)誘發(fā)電位少次提?。跩］.中國(guó)生物醫(yī)學(xué)工程學(xué)報(bào)，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ā)電位信號(hào)去噪研究［J］.計(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.