基于小波分析方法的脑电诱发电位单导少次提取

王永轩 邱天爽 刘 蓉

（大连理工大学 电子信息与电气工程学部，大连 116024）

引言

脑电诱发电位（evoked potential，EP）包含了有关神经系统传导通路的信息，因此有着重要的临床应用价值。由头皮表面测得的EP信号叠加了自发脑电（electroencephalogram，EEG）等背景噪声［1］，信噪频谱重叠且信噪比非常低，提取难度较大。因为EP信号对外界刺激具有锁时性，所以叠加平均法（ensemble average，EA）和加权叠加平均法被广泛应用在临床EP提取中［1］。该方法简单且易于实现，但需要几百次甚至上千次刺激才能得到较为理想的EP信号。为了实现快速提取、不丢失细节以及不受EP信号的时变影响，少次提取成为一个重要的课题。此外，诱发电位是由特定位置的大脑皮层细胞产生的反应，传导到头皮不同位置时相当于经过了不同的系统，因此诱发电位的单导提取更显得尤为重要。

在过去的几十年里EP信号少次提取的研究取得了长足的进步［2-4］，包括PSR算法［5］、自适应滤波法［6］、高阶累积量法［7］、神经网络法［8-9］、独立分量分析［10-13］、小波分析［14-17］等等。PSR算法通过FFT将时域信号变换到频域再进行处理，最终设计出最佳滤波器对观测信号进行滤波。自适应滤波法使用自适应算法来改变滤波器的参数和结构，这样可以跟踪观测信号的时变特征。高阶累积量法是利用高阶累积量对高斯噪声及部分有色噪声不敏感的特性，设计出一个有限冲激响应匹配滤波器，从而达到增强有用信号同时抑制噪声的作用。神经网络是非线性处理过程，能够适应诱发电位的非线性特性。独立分量分析是伴随着解决盲源分离问题发展起来的多导信号处理方法，它是以信号之间的独立性为前提，可以分离叠加在一起的独立信号。小波分析是一种时频分析方法，可以同时反映出信号的全局特征和局部特征，这样可以根据诱发电位的性质对小波系数进行选取。针对诱发电位与自发脑电的特性，本研究提出了使用预白化与小波分析相结合的方法来实现对脑电诱发电位的单导少次提取。

1 自发脑电建模与预白化

在头皮表面获得的观测信号除了诱发电位和自发脑电以外，还包括工频、心电、肌电、眼电等干扰以及各种噪声。由于目前对脑电信号的去噪方法已经比较完善，所以可以假设非脑电的信号已被除去，观测信号只是诱发电位与自发脑电的叠加。因此观测信号的模型可表示为

式中，x（n）、s（n）和v（n）分别表示观测到的带噪信号、纯净EP信号和自发脑电，n为采样点序号，N为单次EP信号长度。在给定的采样频率下，自发脑电表现为有色噪声，用AR模型表示为

式中，P为滤波器的阶数，w（n）为激励白噪声。滤波器的系统传递函数为

在实际处理中为了得到自发脑电的AR模型，需要一段在刺激前记录的观测信号，即纯净的自发脑电，令其为v0（m），（m=0，1，…，M-1），这样可以由v0（m）得到自发脑电AR模型的参数估计。自发脑电在采样频率为fs=1kHz时典型的AR模型为［18］

为了应用小波分析方法去噪，首先要对观测信号进行滤波处理，使其中的自发脑电白化为白噪声。将观测信号经过一个系统传递函数为G1（z）=1/G0（z）=A（z）的滤波器，则其中的自发脑电将被白化为白噪声，对应的差分方程为

经过了白化滤波器后，诱发电位成为期望信号u（n）

同时观测信号成为y（n），并有

由于s（n）与v（n）是相互独立的，所以u（n）与w（n）也是独立的。同时由于v（n）是零均值的平稳有色噪声，所以w（n）是零均值的平稳白噪声。经过白化后就可以通过小波分析方法从y（n）中估计出u（n），最后由u（n）恢复s（n）。

2 小波分析与去噪

2.1 小波分析原理

小波分析是一种时频分析方法，该方法需首先确定一个母小波，再通过母小波的伸缩平移族得到一组基。作为一种分解分离方法，小波分析包括分解、处理和综合三个过程。分解就是将原始信号向这组基上进行投影，得到不同尺度和平移参量下的系数，即小波系数。然后根据需要对小波系数进行加工处理。将处理后的小波系数进行反变换即为综合过程。小波分析的另一重要概念是多分辨率分析，可使用Mallat快速算法逐级计算。一个确定的小波分析相当于一个双通道滤波器组，包括低通滤波器H0（z）和高通滤波器H1（z），在不同的分辨率级下滤波器组是相同的。分解过程中的第j+1级的尺度系数cj+1（k）和小波系数dj+1（k），可由第j级的尺度系数cj（k）按式（8）得到，逐级递推计算可得到所有的小波系数。

综合过程中第j级的尺度系数cj（k），可由第j+1级的尺度系数cj+1（k）和小波系数dj+1（k）按式（9）得到，直到恢复原始信号。

以上分解和综合中使用的滤波器组是相同的，因此这是正交小波分析。在正交小波分析中，从原始信号空间到小波系数空间的投影矩阵是正交阵。将序列y（n）、u（n）和w（n）表示为向y、u和w，再分别将它们在各个分辨率级上的所有小波系数组合在一起，用向量Y、U和W表示，则有Y=TTy、U=TTu和W=TTw，其中T为正交阵。以Harr正交小波分析为例说明，在Harr正交小波分析中有h0（n）=｛1/，1/｝，h1（n）=｛1/，-1/｝。设原始信号为y=［y（0），y（1），...y（7）］T，经过Mallat快速算法得到的第1、2、3级的小波系数和第3级的尺度系数组成的向量为

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

则T为

所有的正交小波分析都存在这样的正交阵T对应这种正交变换关系。对于白噪声w，有E（wwT）=σ2I，其中σ为白噪声功率的均方根，所以E（WWT）=E（TTwwTT）=σ2I。这说明白噪声经正交小波变换后得到的小波系数在各个分辨率级上仍表现为白噪声，即白噪声有限的能量将均匀分散在所有的小波系数上。而对于期望信号u，由于并未被白化而有一定的规则性，以及小波分析的聚焦作用，小波系数非零值将集中在少数系数上，大多数系数接近于零。这样通过对小波系数进行处理后再进行反变换，则可以在很大程度上降低噪声能量，提取出期望信号。

2.2 小波去噪

为实现小波去噪，应设计一个非线性函数对Y进行处理，使非线性函数的输出尽可能接近U，通用的方法主要是阈值法［8］，但阈值的选取对结果影响较大，难以确定合适的取值。本研究采用加权结合阈值的方法，目标是使重构信号与期望信号u的均方误差达到最小。令加权系数向量为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（n），然后由式（13）、（14）得到加权系数b（n），相乘可得U（n）的估计值（n）=b（n）·Y（n），最后由式（9）进行反变换得到u（n）的估计值（n）。在自发脑电建模与预白化部分中，已经提到需要一段在刺激前记录的观测信号，即不包含诱发电位成分的纯净的自发脑电，由这段观测信号被白化后得到的白噪声，就可以估计出上述白噪声功率均方根σ。

3 恢复EP信号

利用小波分析方法从y（n）中估计出u（n）后需要由u（n）恢复出s（n），可以采用两种方法，一是频域法，二是矩阵法，理论和仿真实验证明两种方法的结果几乎完全一样，只是在起点处有微小的差别。频域法是首先将式（6）中a（n）=｛1，a1，a2，…，aP｝补零延拓成N点，然后对a（n）和（n）进行FFT变换得到A（k）和（k），最后对（k）=（k）/A（k）进行IFFT得到（n），则（n）就是对s（n）的估计结果。矩阵法是将式（6）由卷积形式改写为矩阵形式u=As，其中系数矩阵A为

显然矩阵A是下三角矩阵，所以一定可逆，令=A-1，则（n）就是对s（n）的估计结果。在实际计算时不用求逆，对于系数为下三角矩阵的矩阵方程可以用迭代方法求解。

4 仿真实验与分析

仿真实验使用的EP信号为经过叠加平均法得到猫的体感诱发电位，有效频带为0～20Hz，峰值频率为10Hz。自发脑电为由式（4）的AR模型产生的有色噪声，有效频带为0～30Hz，峰值频率为13Hz。实验中采样频率为1kHz，采样时间为512ms，有色噪声符合自发脑电的特性，并且在频带上完全覆盖了EP信号，初始信噪比（SNRi）由0dB逐渐降低到-10dB。算法使用的小波为db3正交小波。

图1和图2分别为SNRi=0dB和SNRi=-10dB条件下小波方法和20次叠加平均法的估计结果。图3给出了在0dB条件下有色噪声和EP信号的原始信号、预白化后信号和小波系数。从图3中可以明显看出有色噪声被白化为白噪声，其小波系数也表现为白噪声特性。而EP信号经白化滤波器后得到的期望信号仍有一定规则性，其小波系数非零值将集中在少数系数上，大多数系数接近于零。图4和图5分别为随着SNRi的降低，不同方法估计结果的信噪比（SNRo）和估计结果与纯净EP信号相关系数（r）的变化情况（数据为50次实验结果的平均值）。

小波分析的优势在于它是时频分析，这样可以将期望信号变换到少数系数上，实现了能量集中。在本算法中预白化是必不可少的，只有在确保噪声是白噪声的前提下小波分析才能发挥它的优势，因为白噪声的小波系数分散在所有值上，这样才能突出期望信号的少数系数。同时由实验结果也可以看出估计的效果受制于信噪比，这是因为虽然白噪声小波系数均匀分散在所有值上，但当噪声足够强时，仍可以淹没期望信号的少数较大系数，从而使期望信号系数估计的误差增大。通过整个估计过程可知，对观测信号小波系数的非线性处理是至关重要的环节，如果已知更多的先验信息就可以更灵活地处理小波系数，从而获得更好的效果。同时也应注意到正交小波的选取也很重要，如果选取的小波在某一分辨率级上与期望信号非常接近，则期望信号的小波系数会更加集中，这样重构信号与期望信号的均方误差会更小。

图1 SNR=0dB条件下的估计结果。（a）有色噪声；（b）纯净EP信号；（c）单次带噪信号；（d）小波方法估计结果；（e）20次带噪信号；（f）20次叠加平均结果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算法不具备时频分析与聚焦信号能量的功能，所以在系数处理和信号重建中只能消除有限的噪声，在信号与噪声频谱重叠的情况下算法性能会退化。自适应滤波法通常是以叠加平均法得到的波形做为参考信号，其性能取决于参考信号和EP的相关程度。独立分量分析是少次提取EP信号的新方法之一，但只能处理多导脑电信号，同时还要做很多假设，包括信号源必须是独立的，信号源数目要不多于观测信号路数，观测信号是信号源的线性混合，混和矩阵是列满秩的等等。目前用于EP信号提取的小波分析方法通常需要已知EP信号的部分先验信息，并且需要人为选择小波系数，这在应用中是不尽可行的。本算法虽不需要EP信号的先验信息，但需要假设观测信号只是诱发电位与自发脑电的叠加，并没有考虑其他高斯/非高斯、周期/非周期性的各种不同噪声。仿真实验表明本算法效果与20次叠加平均效果相当。

图2 SNR=-10dB条件下的估计结果。（a）有色噪声；（b）纯净EP信号；（c）单次带噪信号；（d）小波方法估计结果；（e）20次带噪信号；（f）20次叠加平均结果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信号的原始信号、白化后信号和小波系数。（a）有色噪声；（b）EP信号；（c）白化后得到的白噪声；（d）期望信号；（e）白噪声小波系数；（f）期望信号小波系数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的降低不同方法估计结果的信噪比（SNRo）的变化情况Fig.4 The relationship between SNRowith SNRi by different method

图5 随着SNRi的降低不同方法估计结果的相关系数（r）的变化情况Fig.5 The relationship between correlation coefficient r with SNRiby different method

5 结论

本研究提出了一种基于小波分析的方法用于脑电诱发电位单导少次提取。该算法首先将有色噪声进行预白化，确保噪声的小波系数分散在所有系数上，而期望信号的小波系数非零值集中在少数系数上。这样通过对小波系数进行处理后再进行反变换，最后恢复诱发电位，可以在一定程度上降低噪声能量。实验表明，采用加权结合阈值的小波分析方法，可以抑制自发脑电的干扰，与初始信噪声比相比，估计结果的信噪比获得了提高。更重要的是本算法摆脱了对EP信号先验信息的依赖性，无需人为选择小波系数，符合实际应用背景。理论分析和仿真试验均表明了该算法的应用于EP信号单导少次提取问题有效性。将在下一步的研究工作中继续完善该算法，并将其应用于实际测试EP信号提取问题中。

［1］潘映辅.临床诱发电位学（第二版）［M］.北京：人民卫生出版社，2000.

［2］朱常芳，胡广书.诱发电位快速提取算法的新进展［J］.国外医学生物医学工程分册，2000，23（4）：211-216.

［3］张金凤，邱天爽.诱发电位波形提取方法及进展［J］.北京生物医学工程，2003，22（4）：3030-307.

［4］陶彩林，邹凌，马正华，等.诱发电位单次提取技术的研究进展［J］.生物医学工程学杂志，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，Francis 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在视觉诱发电位的少次提取与波形分析中的应用［J］.中国生物医学工程学报，2000，19（3）：334-341.

［12］毕晓辉，邱天爽，朱勇，等.基于经验模式分解和独立分量分析的单道少次EP信号提取［J］.中国生物医学工程学报，2008，27（6）：817-821.

［13］师黎，钟丽辉，王瑞.基于虚拟通道ICA-WT大鼠视觉诱发电位少次提取［J］.中国生物医学工程学报，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］邹凌，陶彩林，王正洪.基于小波变换的诱发电位信号去噪研究［J］.计算机应用与软件，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.