本文應用小波變換方法從背景噪聲中提取出腦干聽覺誘發電位(BAEP)信號并進行特征識別。首先我們討論了母小波和小波變換算法的選擇, 發現雙正交母小波bior5.5和穩定離散小波變換(SWT)最適合BAEP信號的小波多分辨分析。通過D6尺度小波系數的相關性分析, 發現相關性大于0.4的單次刺激記錄具有較高的信噪比, 使得僅數次刺激的疊加平均就能清晰地識別出BAEP信號的各個波。最后我們用此方法來挑選各次刺激記錄, 對每10次記錄進行疊加平均和小波濾波結合消噪, 并正確識別和計算出BAEP信號各個波的潛伏期。實驗證明本方法通過小波系數相關性分析能有效選取單次刺激的BAEP記錄, 在大大減少刺激次數的同時, 達到了更好的消噪效果。
引用本文： 田福英, 孫迎. 基于小波多分辨分析的腦干聽覺誘發電位特征提取. 生物醫學工程學雜志, 2015, 32(3): 514-519. doi: 10.7507/1001-5515.20150094 復制
引言
聽覺誘發電位(auditory evoked potential, AEP)是指聽覺系統在接受聲音刺激后，從耳蝸至各級聽覺中樞產生的相應電活動。腦干聽覺誘發電位(brainstem auditory evoked potential, BAEP)是聽覺誘發電位最早的部分，發生在聲刺激后10 ms以內，主要包含6~7個陽性波，標記為Ⅰ波~Ⅶ波，其中Ⅰ波~Ⅴ波比較穩定[1-2]。BAEP作為一個聽力評估的客觀方法已廣泛應用于臨床。由于BAEP信號不受注意力、睡意以及藥物等影響，因此非常適合嬰幼兒與某些不配合人群的聽力檢測。但是，通過電極記錄到的BAEP信號非常微弱，通常淹沒在自發腦電、肌電、眨眼以及環境、儀器本身所產生的大量背景噪聲中，信噪比非常低。提取BAEP信號最常用的方法是疊加平均，即通過施加多次刺激，然后以刺激時刻為準，計算每次刺激后記錄到的信號的總體平均，根據信噪比的不同，疊加次數通常需要達到1 000~2 000次。這種方法雖然簡單，但是完全忽略了每次刺激后的BAEP信號間的變異性[3]，而某些背景噪聲如腦電不是完全隨機的；另外, 大量的刺激也容易造成神經系統的疲勞和適應性從而影響BAEP信號。
小波變換是一種同時具有時域和頻域的良好局部化性質的時頻分析方法，它提供了信號在不同時間、不同尺度(頻率)上的信息，能有效實現時頻多分辨分析。由于小波變換對信號沒有平穩性要求，已被成功應用于心電、腦電等很多電生理信號的分析[4-5]。而近年來，如何利用小波分析工具來提取BAEP信號特征也成為一個研究熱點[6-9]。
1 數據獲取方法
1.1 BAEP信號獲取
本文所分析的數據由一臺自主研發的聽覺誘發電位采集儀記錄，與上海交通大學醫學院附屬新華醫學合作，實驗在一個隔音電屏蔽室內完成。對于臨床基于疊加平均算法的BAEP檢測儀，由于其檢測效果依賴于刺激次數，因此理論上提高刺激頻率就可以縮短檢測時間，但是當刺激頻率大于30 Hz時會引起疊加結果的反應波分化程度減弱，根據合作醫生的臨床經驗常采用22 Hz。為了滿足挖掘新的BAEP提取算法，進一步縮短檢測時間，實驗中又采用了41 Hz的刺激頻率。本次實驗的受試者為3名正常人，采用的短聲刺激頻率分別為22 Hz與41 Hz兩種，刺激聲強均為50 dB。每刺激一次記錄到一個BAEP波形，每次實驗刺激1 000次，即記錄1 000個BAEP波形。采樣頻率為100 kHz。本文分析的數據為5次實驗記錄到的5組數據。
1.2 BAEP信號的頻譜特征
應用小波變換對信號進行時頻分析時，首先需要了解信號所包含的頻率成分。圖 1所示為一例BAEP測量數據的1 000次疊加平均后的信號及其頻譜圖。
圖1
1 000次BAEP信號的疊加平均及其頻譜圖
Figure1.
Ensemble averaging of 1 000 single BAEP trails and its frequency spectrum
由于BAEP信號頻率范圍為100 Hz~3 kHz，因此我們只畫出了100 Hz~3 kHz的頻譜圖，從圖中可以看出280、450和980 Hz附近出現峰值，這正好代表著BAEP中Ⅴ波的主要頻率成分[10-11]。
2 原理與方法
2.1 小波變換基本原理
小波是一個平滑、快速衰減的振蕩波，在時域和頻域都有很好的定位。小波簇是滿足相容性條件的基本母小波ψ(t)經過位移和膨脹后的一系列小波：
| $ {{\psi }_{a, b}}\left(t \right)={{\left| a \right|}^{-1/2}}\psi \left(\frac{t-b}{a} \right) $ |
a, b∈R，a≠0分別是膨脹參數和位移參數，t表示時間。當a增大時，小波變窄。所以，每個小波都是母小波在不同時刻不同尺度上重復的唯一分析波。而對于信號S(t)∈L2(R)(實平方可積函數空間)的連續小波變換定義為函數S(t)與在不同a, b上的小波簇ψa, b的卷積：
| $ \left({{\text{W}}_{\psi }}\text{S} \right)\left(\text{a, b} \right)={{\left| a \right|}^{-1/2}}\int\limits_{-\infty }^{\infty }{\text{S}\left(\text{t} \right){{\psi }^{*}}}\left(\frac{\text{t-b}}{\text{a}} \right)d\text{t}=\left\langle \text{S}, {{\psi }_{a, b}} \right\rangle $ |
若對式(1)中小波函數的尺度參數和位移參數離散化，取aj=2-j，bj, k=2-jk，其中j, k∈Z，則得到離散的二進小波：
| $ {{\psi }_{j.k}}\left(t \right)={{2}^{j/2}}\psi \left({{2}^{j}}t-k \right)\ \ j, k\in Z $ |
那么根據式(2)，離散小波變換(discrete wavelet transform, DWT)定義為：
| $ \text{W}{{\text{T}}_{\text{S}}}\left({{2}^{j}}, k \right)=\int{s\left(t \right)\psi _{j, k}^{*}\left(t \right)}\text{d}t $ |
2.2 小波變換算法和母小波的選擇
在用小波變換方法分析信號時需要重點考慮兩個問題。一是選擇合適的小波變換算法，以減小時移變化；二是選擇合適的母小波，它關系到小波重構后波形的平滑度和相位變型問題。
對稱的母小波具有線性相位，能減小重構信號的相位扭曲程度，而母小波自身的連續性和平滑性則決定著重構信號的平滑度。Wilson等的研究證實雙正交母小波bior5.5同時具有較好的對稱性和平滑度，最適合AEP信號的多分辨分析[8]。
Wilson也指出最常用的二進DWT算法的二次采樣，將產生重構信號的移位變化，也就是說它并不是不變的，而這一點對BAEP信號各波潛伏期的研究來說是非常重要的。平穩小波變換(stationary wavelet transform, SWT)算法[12]在二進DWT算法的基礎上做了改進，它通過對濾波器插值的方法減小重構信號的移位變化性。如圖 2a，對長度為N的原始信號s，SWT首先通過s與低通濾波器Lo_D和高通濾波器Hi_D的卷積得到一尺度的近似分量cA1和細節分量cD1，與DWT的二次采樣不同，這里的cA1和cD1的長度也是N(與原始信號s長度一樣)，而不是N/2。然后用同樣的方法我們對cA1做二尺度分解得到cA2和cD2，以此類推也可以得到更高尺度的分解系數，如圖 2b，當j=0時，cA0=s。但是卷積用的濾波器是通過對上一尺度的濾波器插值后得到，如圖 2c所示，而當j=0時，F0=Lo_D，G0=Hi_D。
圖2
SWT分解算法
Figure2.
Algorithm of SWT
如圖 3所示，我們用470 Hz和280 Hz的兩個正弦波疊加后得到一個模擬的AEP信號，并采用DWT和SWT兩種算法對其進行8尺度的分解，然后對包含信號頻率成分的D7、D8兩個尺度進行重構。從圖中可以看出，用SWT算法重構后的信號幾乎與原始信號重合，平均絕對誤差(mean absolute error, MAE)為0.01，而DWT算法重構后的信號與原始信號有較大偏差，其平均絕對誤差為0.15。因此SWT比DWT更適于AEP信號的分析。
圖3
模擬BAEP信號及分別用DWT和SWT算法計算的D7、D8尺度的重構信號
Figure3.
Simulated BAEP of windowed sinusoids of 280 and 470 Hz and the D7, D8 reconstruction of the simulated BAEP by DWT and SWT respectively
2.3 相關性分析原理
相關關系是變量之間的一種非確定性關系，相關系數則表示了研究變量之間的線性相關程度。其定義為：，其中Cov(x, y)=E{[X-E(X)][Y-E(Y)]}為變量X、Y的協方差，而、分別表示變量X與Y的均方差。ρXY是一個可以用來表征X、Y之間線性關系緊密程度的量，｜ρXY｜≤1。當｜ρXY｜較大時，說明X、Y之間的線性相關程度較好；當｜ρXY｜較小時，說明X、Y線性相關的程度較差；特別當｜ρXY｜=1時，表示X、Y之間以概率1存在著線性關系。本文將計算單刺激BAEP信號與參考BAEP信號在有效小波尺度上的相關性。
3 BAEP信號提取結果與討論
3.1 BAEP信號的小波分解與重構
我們用bior5.5小波和SWT算法對1 000次疊加平均后的BAEP信號進行多分辨分析，如圖 4。
圖4
1 000次疊加平均BAEP信號的小波分解及相應尺度的信號重構
Figure4.
Wavelet coefficients for ensemble average of 1 000 BAEP trails and their corresponding reconstructed signals
由于本文中獲取BAEP信號的采樣頻率為100 kHz，那么對BAEP信號做8尺度的分解后，BAEP信號的頻率成分主要集中在D6(781~1 562 Hz)、D7(390~781 Hz)、D8(195~390 Hz)。圖 4左邊是各個尺度上的小波系數，右邊實線表示各尺度系數對應的重構信號，而虛線為1 000次疊加后的BAEP信號，在重構信號中我們可以清晰地分辨出BAEP信號的各個波，而且它們的潛伏期與原始的疊加平均BAEP信號有非常好的一致性。
3.2 單刺激BAEP信號的選取
BAEP的各個波主要出現在小波分解的D6、D7、D8和A8尺度，但是噪聲信號的出現是覆蓋整個時域和頻域的。圖 5是任意四個單刺激BAEP信號(藍色)和1 000次疊加平均BAEP信號(黑色)各尺度上的小波系數。圖中可以看出，單個刺激的BAEP信號混有強大的背景噪聲，即使在D6、D7、D8和A8尺度，其小波系數和1 000次疊加平均BAEP信號的小波系數也很難有一致性。一致性越好，則該刺激的BAEP信號混入相應頻段的噪聲越小；相反一致性越差，則該刺激的BAEP信號混入相應頻段的噪聲越大。頻段越高越容易混入噪聲。由此若以1 000次疊加平均作為一個相對干凈的BAEP參考信號，我們認為可通過單次刺激和它在D6尺度系數上的相關性，來選擇“好的”BAEP信號。于是我們計算了每個單刺激BAEP信號與1 000次疊加平均BAEP信號在D6尺度上的相關性，并按相關系數的大小進行分類。
圖5
單個刺激BAEP信號的小波系數(藍色)與1 000次疊加平均BAEP信號的小波系數(黑色)對比
Figure5.
Comparison of wavelet coefficients for single trail BAEPs (blue) and ensemble average of 1 000 BAEP trails (black)
圖 6為處于各個相關性等級的10組AEP信號的疊加平均，以及保留D6、D7、D8、A8尺度(D1~D5為高頻噪聲，系數置零)的小波重構信號。其中由于相關系數大于0.8的BAEP信號只有兩個，所以圖 6a為兩次刺激的BAEP信號的疊加平均和小波重構結果。如圖所示，相關系數大于0.4時(圖 6a、b、c、d、e)，僅用10次刺激的BAEP信號，就可以清楚地識別BAEP信號的各個波，每個波的潛伏期與1 000次疊加平均有很好的一致性。而當相關系數小于0.4時(圖 6f、g、h、i)，就很難從10次BAEP信號的處理結果中識別出各個波了。圖 7是所有D6尺度系數相關性大于0.4的記錄(229次)消噪結果與所有1 000次疊加平均的對比。可以看出，去掉D6尺度系數相關性小于0.4的記錄，雖然疊加次數大大減少，卻比1 000次的疊加結果更加清晰，容易識別。
圖6
1 000次疊加平均BAEP信號(藍)、具有近似相關性的10個BAEP信號的疊加平均信號(黑)及其小波消噪后的重構信號(紅)
Figure6.
Ensemble average of 1 000 BAEP trails (blue), average of 10 BAEP trails with similar correlation (black) and its wavelet reconstruction signal (red)
圖7
1 000次BAEP信號的疊加平均(藍)、相關性大于0.4的所有BAEP記錄(229次)的疊加平均(黑)及其小波重構信號(紅)
Figure7.
Ensemble average of 1 000 BAEP trails (blue), average of 229 BAEP trails whose correlation was larger than 0.4 (black), and its wavelet reconstruction signal (red)
3.3 BAEP提取和識別結果分析
基于以上分析結果，我們對所有實驗數據進行去噪和識別的處理。首先我們以小波分解的D6尺度系數的相關性為準則，對各次刺激記錄進行挑選，保留相關系數大于0.4的記錄，然后通過疊加平均和小波濾波方法提取出清晰的BAEP信號，最后對BAEP各個波的波峰進行準確地標識，并計算了潛伏期。表 1是分別對5次實驗中所有相關系數大于0.4的記錄的處理結果，我們計算了BAEP各波峰的潛伏期，其中最后一列為相關系數大于0.4的單刺激記錄數。為了在保證準確識別AEP各波的同時，進一步減少疊加次數，我們對每次試驗的所有D6尺度系數相關性大于0.4的記錄進行分組(每10條記錄一組)，用同樣方法對每一組數據消噪后，識別出BAEP的各個波峰并計算各波潛伏期。如表 2所示，我們計算了每次實驗中所有組BAEP各波潛伏期的均值和標準差。通過表 1與表 2可以看出，10個單刺激記錄消噪后也能正確地定位BAEP信號的各波，特別是Ⅴ波，其識別率達到100%。而標準差值不為0也正好說明每個單刺激BAEP之間是存在變異性的。
表1
BAEP信號消噪后波峰識別與潛伏期計算
Table1.
Wave identification and latency results for all BAEP trails (correlation>0.4) after de-nosing
表2
10次BAEP記錄(相關系數＞0.4)消噪后的波峰識別與潛伏期計算
Table2.
Wave identification and latency results for every 10 BAEP trails (correlation>0.4) after de-nosing
4 結論
本文介紹了一種基于小波多分辨分析的BAEP消噪方法，僅需數次刺激即可提取出清晰的BAEP信號，并正確標識臨床關心的Ⅰ~Ⅴ波，以及計算各波的潛伏期。實驗證明雙正交母小波bior5.5的選擇和SWT算法的應用非常適合BAEP信號的多分辨分析。通過單刺激記錄與“干凈BAEP”在小波信號尺度上的相關性分析，可以有效剔除“不好”的單刺激記錄，提高去噪的效率，從而減少所需刺激次數并提高檢測速率。為了挑選“好的”單刺激BAEP信號，該方法中“干凈BAEP”信號可以通過實驗初期的自學習來獲取，本文為1 000次疊加后得到，進一步的研究中我們可以通過自適應的方法減少參考信號的學習時間, 從而真正地減少刺激次數。
引言
聽覺誘發電位(auditory evoked potential, AEP)是指聽覺系統在接受聲音刺激后，從耳蝸至各級聽覺中樞產生的相應電活動。腦干聽覺誘發電位(brainstem auditory evoked potential, BAEP)是聽覺誘發電位最早的部分，發生在聲刺激后10 ms以內，主要包含6~7個陽性波，標記為Ⅰ波~Ⅶ波，其中Ⅰ波~Ⅴ波比較穩定[1-2]。BAEP作為一個聽力評估的客觀方法已廣泛應用于臨床。由于BAEP信號不受注意力、睡意以及藥物等影響，因此非常適合嬰幼兒與某些不配合人群的聽力檢測。但是，通過電極記錄到的BAEP信號非常微弱，通常淹沒在自發腦電、肌電、眨眼以及環境、儀器本身所產生的大量背景噪聲中，信噪比非常低。提取BAEP信號最常用的方法是疊加平均，即通過施加多次刺激，然后以刺激時刻為準，計算每次刺激后記錄到的信號的總體平均，根據信噪比的不同，疊加次數通常需要達到1 000~2 000次。這種方法雖然簡單，但是完全忽略了每次刺激后的BAEP信號間的變異性[3]，而某些背景噪聲如腦電不是完全隨機的；另外, 大量的刺激也容易造成神經系統的疲勞和適應性從而影響BAEP信號。
小波變換是一種同時具有時域和頻域的良好局部化性質的時頻分析方法，它提供了信號在不同時間、不同尺度(頻率)上的信息，能有效實現時頻多分辨分析。由于小波變換對信號沒有平穩性要求，已被成功應用于心電、腦電等很多電生理信號的分析[4-5]。而近年來，如何利用小波分析工具來提取BAEP信號特征也成為一個研究熱點[6-9]。
1 數據獲取方法
1.1 BAEP信號獲取
本文所分析的數據由一臺自主研發的聽覺誘發電位采集儀記錄，與上海交通大學醫學院附屬新華醫學合作，實驗在一個隔音電屏蔽室內完成。對于臨床基于疊加平均算法的BAEP檢測儀，由于其檢測效果依賴于刺激次數，因此理論上提高刺激頻率就可以縮短檢測時間，但是當刺激頻率大于30 Hz時會引起疊加結果的反應波分化程度減弱，根據合作醫生的臨床經驗常采用22 Hz。為了滿足挖掘新的BAEP提取算法，進一步縮短檢測時間，實驗中又采用了41 Hz的刺激頻率。本次實驗的受試者為3名正常人，采用的短聲刺激頻率分別為22 Hz與41 Hz兩種，刺激聲強均為50 dB。每刺激一次記錄到一個BAEP波形，每次實驗刺激1 000次，即記錄1 000個BAEP波形。采樣頻率為100 kHz。本文分析的數據為5次實驗記錄到的5組數據。
1.2 BAEP信號的頻譜特征
應用小波變換對信號進行時頻分析時，首先需要了解信號所包含的頻率成分。圖 1所示為一例BAEP測量數據的1 000次疊加平均后的信號及其頻譜圖。
圖1
1 000次BAEP信號的疊加平均及其頻譜圖
Figure1.
Ensemble averaging of 1 000 single BAEP trails and its frequency spectrum
由于BAEP信號頻率范圍為100 Hz~3 kHz，因此我們只畫出了100 Hz~3 kHz的頻譜圖，從圖中可以看出280、450和980 Hz附近出現峰值，這正好代表著BAEP中Ⅴ波的主要頻率成分[10-11]。
2 原理與方法
2.1 小波變換基本原理
小波是一個平滑、快速衰減的振蕩波，在時域和頻域都有很好的定位。小波簇是滿足相容性條件的基本母小波ψ(t)經過位移和膨脹后的一系列小波：
| $ {{\psi }_{a, b}}\left(t \right)={{\left| a \right|}^{-1/2}}\psi \left(\frac{t-b}{a} \right) $ |
a, b∈R，a≠0分別是膨脹參數和位移參數，t表示時間。當a增大時，小波變窄。所以，每個小波都是母小波在不同時刻不同尺度上重復的唯一分析波。而對于信號S(t)∈L2(R)(實平方可積函數空間)的連續小波變換定義為函數S(t)與在不同a, b上的小波簇ψa, b的卷積：
| $ \left({{\text{W}}_{\psi }}\text{S} \right)\left(\text{a, b} \right)={{\left| a \right|}^{-1/2}}\int\limits_{-\infty }^{\infty }{\text{S}\left(\text{t} \right){{\psi }^{*}}}\left(\frac{\text{t-b}}{\text{a}} \right)d\text{t}=\left\langle \text{S}, {{\psi }_{a, b}} \right\rangle $ |
若對式(1)中小波函數的尺度參數和位移參數離散化，取aj=2-j，bj, k=2-jk，其中j, k∈Z，則得到離散的二進小波：
| $ {{\psi }_{j.k}}\left(t \right)={{2}^{j/2}}\psi \left({{2}^{j}}t-k \right)\ \ j, k\in Z $ |
那么根據式(2)，離散小波變換(discrete wavelet transform, DWT)定義為：
| $ \text{W}{{\text{T}}_{\text{S}}}\left({{2}^{j}}, k \right)=\int{s\left(t \right)\psi _{j, k}^{*}\left(t \right)}\text{d}t $ |
2.2 小波變換算法和母小波的選擇
在用小波變換方法分析信號時需要重點考慮兩個問題。一是選擇合適的小波變換算法，以減小時移變化；二是選擇合適的母小波，它關系到小波重構后波形的平滑度和相位變型問題。
對稱的母小波具有線性相位，能減小重構信號的相位扭曲程度，而母小波自身的連續性和平滑性則決定著重構信號的平滑度。Wilson等的研究證實雙正交母小波bior5.5同時具有較好的對稱性和平滑度，最適合AEP信號的多分辨分析[8]。
Wilson也指出最常用的二進DWT算法的二次采樣，將產生重構信號的移位變化，也就是說它并不是不變的，而這一點對BAEP信號各波潛伏期的研究來說是非常重要的。平穩小波變換(stationary wavelet transform, SWT)算法[12]在二進DWT算法的基礎上做了改進，它通過對濾波器插值的方法減小重構信號的移位變化性。如圖 2a，對長度為N的原始信號s，SWT首先通過s與低通濾波器Lo_D和高通濾波器Hi_D的卷積得到一尺度的近似分量cA1和細節分量cD1，與DWT的二次采樣不同，這里的cA1和cD1的長度也是N(與原始信號s長度一樣)，而不是N/2。然后用同樣的方法我們對cA1做二尺度分解得到cA2和cD2，以此類推也可以得到更高尺度的分解系數，如圖 2b，當j=0時，cA0=s。但是卷積用的濾波器是通過對上一尺度的濾波器插值后得到，如圖 2c所示，而當j=0時，F0=Lo_D，G0=Hi_D。
圖2
SWT分解算法
Figure2.
Algorithm of SWT
如圖 3所示，我們用470 Hz和280 Hz的兩個正弦波疊加后得到一個模擬的AEP信號，并采用DWT和SWT兩種算法對其進行8尺度的分解，然后對包含信號頻率成分的D7、D8兩個尺度進行重構。從圖中可以看出，用SWT算法重構后的信號幾乎與原始信號重合，平均絕對誤差(mean absolute error, MAE)為0.01，而DWT算法重構后的信號與原始信號有較大偏差，其平均絕對誤差為0.15。因此SWT比DWT更適于AEP信號的分析。
圖3
模擬BAEP信號及分別用DWT和SWT算法計算的D7、D8尺度的重構信號
Figure3.
Simulated BAEP of windowed sinusoids of 280 and 470 Hz and the D7, D8 reconstruction of the simulated BAEP by DWT and SWT respectively
2.3 相關性分析原理
相關關系是變量之間的一種非確定性關系，相關系數則表示了研究變量之間的線性相關程度。其定義為：，其中Cov(x, y)=E{[X-E(X)][Y-E(Y)]}為變量X、Y的協方差，而、分別表示變量X與Y的均方差。ρXY是一個可以用來表征X、Y之間線性關系緊密程度的量，｜ρXY｜≤1。當｜ρXY｜較大時，說明X、Y之間的線性相關程度較好；當｜ρXY｜較小時，說明X、Y線性相關的程度較差；特別當｜ρXY｜=1時，表示X、Y之間以概率1存在著線性關系。本文將計算單刺激BAEP信號與參考BAEP信號在有效小波尺度上的相關性。
3 BAEP信號提取結果與討論
3.1 BAEP信號的小波分解與重構
我們用bior5.5小波和SWT算法對1 000次疊加平均后的BAEP信號進行多分辨分析，如圖 4。
圖4
1 000次疊加平均BAEP信號的小波分解及相應尺度的信號重構
Figure4.
Wavelet coefficients for ensemble average of 1 000 BAEP trails and their corresponding reconstructed signals
由于本文中獲取BAEP信號的采樣頻率為100 kHz，那么對BAEP信號做8尺度的分解后，BAEP信號的頻率成分主要集中在D6(781~1 562 Hz)、D7(390~781 Hz)、D8(195~390 Hz)。圖 4左邊是各個尺度上的小波系數，右邊實線表示各尺度系數對應的重構信號，而虛線為1 000次疊加后的BAEP信號，在重構信號中我們可以清晰地分辨出BAEP信號的各個波，而且它們的潛伏期與原始的疊加平均BAEP信號有非常好的一致性。
3.2 單刺激BAEP信號的選取
BAEP的各個波主要出現在小波分解的D6、D7、D8和A8尺度，但是噪聲信號的出現是覆蓋整個時域和頻域的。圖 5是任意四個單刺激BAEP信號(藍色)和1 000次疊加平均BAEP信號(黑色)各尺度上的小波系數。圖中可以看出，單個刺激的BAEP信號混有強大的背景噪聲，即使在D6、D7、D8和A8尺度，其小波系數和1 000次疊加平均BAEP信號的小波系數也很難有一致性。一致性越好，則該刺激的BAEP信號混入相應頻段的噪聲越小；相反一致性越差，則該刺激的BAEP信號混入相應頻段的噪聲越大。頻段越高越容易混入噪聲。由此若以1 000次疊加平均作為一個相對干凈的BAEP參考信號，我們認為可通過單次刺激和它在D6尺度系數上的相關性，來選擇“好的”BAEP信號。于是我們計算了每個單刺激BAEP信號與1 000次疊加平均BAEP信號在D6尺度上的相關性，并按相關系數的大小進行分類。
圖5
單個刺激BAEP信號的小波系數(藍色)與1 000次疊加平均BAEP信號的小波系數(黑色)對比
Figure5.
Comparison of wavelet coefficients for single trail BAEPs (blue) and ensemble average of 1 000 BAEP trails (black)
圖 6為處于各個相關性等級的10組AEP信號的疊加平均，以及保留D6、D7、D8、A8尺度(D1~D5為高頻噪聲，系數置零)的小波重構信號。其中由于相關系數大于0.8的BAEP信號只有兩個，所以圖 6a為兩次刺激的BAEP信號的疊加平均和小波重構結果。如圖所示，相關系數大于0.4時(圖 6a、b、c、d、e)，僅用10次刺激的BAEP信號，就可以清楚地識別BAEP信號的各個波，每個波的潛伏期與1 000次疊加平均有很好的一致性。而當相關系數小于0.4時(圖 6f、g、h、i)，就很難從10次BAEP信號的處理結果中識別出各個波了。圖 7是所有D6尺度系數相關性大于0.4的記錄(229次)消噪結果與所有1 000次疊加平均的對比。可以看出，去掉D6尺度系數相關性小于0.4的記錄，雖然疊加次數大大減少，卻比1 000次的疊加結果更加清晰，容易識別。
圖6
1 000次疊加平均BAEP信號(藍)、具有近似相關性的10個BAEP信號的疊加平均信號(黑)及其小波消噪后的重構信號(紅)
Figure6.
Ensemble average of 1 000 BAEP trails (blue), average of 10 BAEP trails with similar correlation (black) and its wavelet reconstruction signal (red)
圖7
1 000次BAEP信號的疊加平均(藍)、相關性大于0.4的所有BAEP記錄(229次)的疊加平均(黑)及其小波重構信號(紅)
Figure7.
Ensemble average of 1 000 BAEP trails (blue), average of 229 BAEP trails whose correlation was larger than 0.4 (black), and its wavelet reconstruction signal (red)
3.3 BAEP提取和識別結果分析
基于以上分析結果，我們對所有實驗數據進行去噪和識別的處理。首先我們以小波分解的D6尺度系數的相關性為準則，對各次刺激記錄進行挑選，保留相關系數大于0.4的記錄，然后通過疊加平均和小波濾波方法提取出清晰的BAEP信號，最后對BAEP各個波的波峰進行準確地標識，并計算了潛伏期。表 1是分別對5次實驗中所有相關系數大于0.4的記錄的處理結果，我們計算了BAEP各波峰的潛伏期，其中最后一列為相關系數大于0.4的單刺激記錄數。為了在保證準確識別AEP各波的同時，進一步減少疊加次數，我們對每次試驗的所有D6尺度系數相關性大于0.4的記錄進行分組(每10條記錄一組)，用同樣方法對每一組數據消噪后，識別出BAEP的各個波峰并計算各波潛伏期。如表 2所示，我們計算了每次實驗中所有組BAEP各波潛伏期的均值和標準差。通過表 1與表 2可以看出，10個單刺激記錄消噪后也能正確地定位BAEP信號的各波，特別是Ⅴ波，其識別率達到100%。而標準差值不為0也正好說明每個單刺激BAEP之間是存在變異性的。
表1
BAEP信號消噪后波峰識別與潛伏期計算
Table1.
Wave identification and latency results for all BAEP trails (correlation>0.4) after de-nosing
表2
10次BAEP記錄(相關系數＞0.4)消噪后的波峰識別與潛伏期計算
Table2.
Wave identification and latency results for every 10 BAEP trails (correlation>0.4) after de-nosing
4 結論
本文介紹了一種基于小波多分辨分析的BAEP消噪方法，僅需數次刺激即可提取出清晰的BAEP信號，并正確標識臨床關心的Ⅰ~Ⅴ波，以及計算各波的潛伏期。實驗證明雙正交母小波bior5.5的選擇和SWT算法的應用非常適合BAEP信號的多分辨分析。通過單刺激記錄與“干凈BAEP”在小波信號尺度上的相關性分析，可以有效剔除“不好”的單刺激記錄，提高去噪的效率，從而減少所需刺激次數并提高檢測速率。為了挑選“好的”單刺激BAEP信號，該方法中“干凈BAEP”信號可以通過實驗初期的自學習來獲取，本文為1 000次疊加后得到，進一步的研究中我們可以通過自適應的方法減少參考信號的學習時間, 從而真正地減少刺激次數。
