本文应用小波变换方法从背景噪声中提取出脑干听觉诱发电位(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次叠加后得到，进一步的研究中我们可以通过自适应的方法减少参考信号的学习时间, 从而真正地减少刺激次数。
