面向心电信号去噪的正交小波构造方法

朱俊江 张远辉

(中国计量大学机电学院，杭州 310018)

面向心电信号去噪的正交小波构造方法

朱俊江*张远辉

(中国计量大学机电学院，杭州 310018)

为更好地去除心电信号的肌电干扰，设计出专门用于心电信号处理的正交小波(心电小波)。首先，为满足正交性要求，对小波滤波器系数进行参数化设计；其次，用参数化的小波滤波器系数逼近QRS模板并构造出目标函数，在此基础上以消失距作为约束条件得到优化模型；最后，采用拉格朗日乘子法对模型进行求解。从光滑性、与心电信号的相似性、重构误差以及频响特性等方面，对构造出的小波性能进行评估，结果表明，相对于db和sym系列的小波，心电小波可以获得更好的综合性能。当将该算法用于MIT-BIH数据库中的心电信号去噪时，能够抑制加性噪声的干扰，使信噪比提高10.32dB。根据应用对象设计正交小波的方法，为其他专用小波设计提供新思路。

心电信号；优化；小波；去噪

引言

临床上，心电信号是心血管疾病判断和治疗重要的辅助工具，准确处理心电信号意义重大。然而，心电信号是一种弱电信号，幅值通常为毫伏级，在采集过程中由仪器本身、人体活动等因素带来的噪声对心电信号的影响不可忽略。心电信号的干扰主要有以下几种：第一种是工频干扰，频率为50 Hz左右，其干扰幅度可达心电信号峰值的50%；第二种是基线漂移，会导致心电波形整体上下波动甚至扭曲，能量主要集中在0.1 Hz左右；第三种是肌电干扰，是肌电收缩所产生的噪声，其特点在于频率分布范围广，频谱特性接近于白噪声，因此成为心电信号去噪中的一个难点。

针对肌电干扰的处理，有学者提出采用数字滤波器[1]和自适应滤波器[2]的方法对数字化的ECG信号进行去噪，其缺点在于忽略了噪声和有用信息在频率上是重叠的，在滤除噪声的同时也会将有用信息滤除。主分量分析法[3-4]基于噪声之间是不相关的思想，通过剔去信号当中相关性较小的成分，完成信号去噪，其不足在于去噪前必须将信号进行对齐或者采用多通道的信号作为处理对象，否则将影响去噪效果，甚至去除掉有用信息。经验模态分解(EMD)[5]将信号分解为若干个按频率高低排列的本证模态函数(intrinsic mode function)，从而根据信号特点自适应性的对噪声进行剔除，其不足是当心电信号中存在异常干扰时会造成模态混叠，影响去噪效果。小波变化可以提供信号在时间和尺度上的分析，具有较强的时频局部化能力，因此成为一种最为常用的心电信号去噪方法[6-7]。基于小波去噪方法的处理，关键是母小波的选择和阈值方法的采用。从目前的研究来看，db小波、coif小波、sym小波以及bior小波都能对心电信号进行分解和重构。最常用的一种小波为db4，因为它和心电信号在形状上最为接近。然而，心脏是一个复杂系统，导致心电信号具有很强的非线性和高度的个性化特点。因此，设计出心电信号专用小波意义重大。

在小波构造方面，原始的Daubechies家族小波构造思路是在满足正交性的前提下具备最大的消失矩，因此其对称性很差。Coiflet小波系统通过减小消失矩获得更多的自由度来满足其他的一些要求，在文献[8]中证明采用因子化方法可以得到一个更加对称的小波系统。近年来，采用优化法构建更加适合具体应用的小波系统的思路被提出。在文献[9]中提出了一种弱化的对称滤波器(partial symmetry)，在对称系数中添加了非对称系数，这些非对称系数相比对称系数值非常小，用来满足小波的正交性要求。文献[10]则构造出一种完全对称的小波，并在此条件下对正交性进行优化。这两种优化方法的思路是：在完全对称或几乎对称的前提下去满足正交性。由于正交多辨析有着严格的数学定义，因此难以在该定义下建立具备特殊性能小波的数学模型。本研究采用一种结构化参数来满足正交性要求，并将QRS模板作为逼近目标，用直接优化法构建出适用于心电信号处理的心电小波。

1 数据来源

实验数据来自于MIT-BIH心率不齐数据库(MIT-BIH arrhythmia database)，包含了48条记录，每条记录持续时长30 min。这些记录都包含两个导联的数据，其中一个为修改后的II 型导联，另外一条数据是V1、V2、V4 或 V5型导联。采样频率为360 Hz，给出的数据用一个通过频率为0.1～100 Hz的带通滤波器进行滤波，量化分辨率为200个点/mV。将所有记录中修改后的II 型导联心电信号随机分为了两个数据集，分别记为数据集A、数据集B。其中，数据集A用于构造心电信号模板，数据集B用于衡量小波去噪效果。

2 方法

2.1 心电小波构造方法

采用优化法对心电小波进行设计，步骤如图1所示。首先，为了构造出正交的小波，本研究采用一种叫晶格结构化[11-13]的方法对小波滤波器系数进行参数化；其次，以小波滤波器系数作为参数，通过逼近QRS模板构造出目标函数，通过限制小波的消失距构造出约束条件；最后，通过求解优化模型，得到心电小波。

图1 心电小波设计步骤Fig.1 The block diagram for ECG wavelets design

2.2.1 小波滤波器系数参数

假设H(z)和G(z)分别为小波正交滤波器组中的低通滤波器和高通滤波器，则

(1)

其中

(2)

归一化常量为

如此，在给定参数α=(α1,α2，…，αL)时，即可得到一个长度为N=2L的滤波器组系数。

2.2.2 优化模型

2.2.2.1 目标函数

正常人的心跳包括PQRST波。其中，在心房除极过程中，心电向量从窦房结指向房室结，除极由右心房至左心房，在心电图上反映为P波，QRS波群反映了左右心室的快速去极化过程，而T波代表心室快速复极化过程。标准的心电信号包括12种导联形式，导联形式不同，心电信号的形状也不尽相同，图2模拟了一个修正后的Ⅱ导联心跳周期。其中，QRS波群振幅最大，是衡量心跳频率最常用的特征，也是心电信号中最显著的特征。通过构造与QRS波形较为相似的小波函数，就可以得到更好的去噪效果。为此，针对数据集A，首先采用截止频率为1Hz的低通滤波器对脉搏信号进行滤波，以将信号中的基线漂移去除，采用极大值法提取到每个脉搏波的极大值点，并以此点的前0.2s为起始点、后0.3s作为终止点，对心电信号进行了截取。将每个心跳进行叠加后取平均，最后形成QRS模板。

图2 一个典型的心电周期Fig.2 A typical ECG cycle

在QRS模板上，均匀选取N个点作为小波滤波器的逼近系数，最终可以建立起优化目标函数，即

(3)

式中，fQRS(n)表示在QRS模板上选取的第n+1个点。

2.2.2.2 约束条件

为保证小波函数和尺度函数的光滑性，小波还需要具备一定的消失距。其定义为，若小波具备P阶消失距，则存在

(4)

上式对于l=0,1，…，P-1均成立。

采用晶格结构化参数构造小波时，对于系数长度为N=2L的小波，其自由度还剩L个，当考虑消失距时，每增加一个消失距，则消耗一个自由度。通过减小消失距获得更多的自由度，从而使小波满足其他的一些要求，因此构造出的优化模型如下：

(5)

2.2.2.3 求解方法

优化模型中的目标函数是二次非线性函数，而约束条件为一次非线性等式。对于这样的数学模型，可以采用拉格朗日乘子法对该数学模型进行求解，列出拉格朗日方程如下：

(6)

其中，优化目标E表示滤波器系数和逼近对象之间的误差，可以进一步表示为

(7)

cn为固定值，令

(8)

αi的导数可以表示为

(9)

式中,Hαi=∂H/∂αi，Mτ，αi=∂Mτ/∂αi，Kαi=∂K/∂αi。

因此，可以采用牛顿插值法对式(9)进行迭代求解。需要注意的是，优化模型是一个非凸问题，所以必须合理地设定初值。为了保证小波滤波器性能，采用文献[14]中的方法对算法的初值进行设定。

2.2 基于小波的心电信号去噪方法

3 结果

3.1 心电小波构造结果

图1提供了设计心电小波的优化算法，该算法以滤波器长度和消失距的个数为输入参数。通过改变滤波器长度和消失距，就能生成不同的滤波器系数，进而构造出不同的小波函数。

例1：首先，根据式(1)、(2)的参数化方法，构造出滤波器系数长度为2L=8的小波滤波器系数。其次，在QRS模板上均匀选取8个点作为逼近对象，并分别设定消失距为2和3，构造出优化目标和约束条件；采用拉格朗日乘子法对优化目标进行求解，最终得到目标函数值分别为0.091 3和0.368 7；将求解得到的未知数αi代入式(1)、(2)，可以得到两组高、低通滤波器系数，如表1所示。利用这两组系数生成小波函数，分别为QRS(8,2)和QRS(8,3)。

例2：首先，根据式(1)、(2)的参数化方法，构造出滤波器系数长度为2L=12的小波滤波器系数。其次，在QRS模板上均匀选取12个点作为逼近对象，并分别设定消失距为3、4和5，构造出优化目标和约束条件；采用拉格朗日乘子法对优化目标进行求解，优化的目标结果分别为0.733、2.967 2和5.263 4；将求解得到的未知数αi代入式(1)、(2)，可以得到3组高低通滤波器系数，如表1所示。利用这3组系数生成小波函数，分别为QRS(12,3)、QRS(12,4)、QRS(12,5)。

在例1和例2中得到高、低通滤波器系数如表1所示。在此基础上，可以按照文献[8]中提供的方法分别生成小波系统的尺度函数和小波函数，其中QRS(8,3)的尺度函数和小波函数如图3所示(虚线部分显示的是db4小波)。

3.2 心电信号去噪实验

采用本文第2.2节所提的去噪方法对数据集B进行去噪，采用信噪比和最小二乘误差作为衡量手段，并以db4、db6作为对比，得到心电小波族的去噪效果，如表2所示。母小波选择QRS(8,3)时，信号的去噪效果如图4所示。

表1 心电小波滤波器系数

图3 db4(虚线)和QRS(8,3)(实线)小波函数和尺度函数。(a)小波函数；(b)尺度函数Fig.3 Wavelet function and scale function of db4 (dash line) and QRS(8,3) (solid line). (a) Wavelet function; (b) scale function

输入噪声/dBQRS(8，2)QRS(8，3)QRS(12，3)db4db6SNRMSESNRMSESNRMSESNRMSESNRMSE-207265546×10-46850597×10-46147596×10-47026595×10-46070735×10-4-127852526×10-47329601×10-47421578×10-47311596×10-46244724×10-4-68869428×10-48199514×10-48194473×10-46554503×10-46211626×10-4

图4 两种小波去噪效果对比。(a)含噪信号；(b)QRS(8,2)小波的去噪效果；(c)db4的去噪效果Fig.4 Denoise result under two different wavelets.(a) Signal with noise; (b) Denoise result of QRS(8,2); (c) Denoised result of db4

4 讨论和结论

从图4中可以看出，采用Q(8,2)小波对心电信号进行去噪，得到的Q波和S波失真较db4有所减小。表3表明，不同的母小波导致了不同的去噪效果。在对心电信号中的白噪声进行消除时，采用QRS(8,2)小波可比db4小波平均提高信噪比10.32 dB，并获得较小的MSE值，表明心电小波在心电信号去噪方面的优势。

采取特殊的小波系数参数化方法，可以通过牺牲正交小波的消失距来提高心电小波与心电信号主特征相似性，促使心电小波在处理心电信号时取得更好的效果。相比多分辨率分析方法，直接优化构造小波的方法采用了优化算法，直接对小波滤波器系数进行求解，构造思路更为简单，因此具有更高的灵活性，能够满足根据特殊应用设计专用小波的需求。用直接优化法构造出心电专用小波，不仅在给定的应用中能够取得较好的效果，而且为其他专用小波设计提供了新思路。

在本研究中，通过直接优化滤波器系数构造心电信号专用小波，主要特点在于：首先，对小波滤波器系数进行参数化，以保证构造出的小波具有正交性。其次，用参数化的小波滤波器系数逼近QRS模板，构造出了目标函数；在此基础上，以消失距作为约束条件，构造出了优化模型。最后，采用拉格朗日乘子法对模型中的非凸问题进行求解，得到了小波滤波器系数。

采用心电专用小波对心电信号去噪，可以改善去噪效果。本研究的后续工作可以从两个方面进行：一是对该专用小波在特征提取、压缩等方面继续深入研究，二是研究如何针对其他信号构造出特殊小波，以提高该信号的处理效果和效率。

[1] Singh N, Ayub S, Saini JP. Design of digital IIR filter for noise reduction in ECG signal [C]// IEEE Computer Society. Computational Intelligence and Communication Networks. Mathura:IEEE, 2013:171-176.

[2] Thakor NV, Zhu Y. Applications of adaptive filtering to ECG analysis: noise cancellation and arrhythmia detection [J]. IEEE Transactions on Biomedical Engineering, 1991,38(8):785-794.

[3] Romero I. PCA and ICA applied to noise reduction in multi-lead ECG[J]. Computing in Cardiology, 2011, 23(1):613-616

[4] Romero I. PCA-based noise reduction in ambulatory ECGs [J]. Computing in Cardiology 2010,37(1): 677-680.

[5] Blanco-Velasco M, Weng BW, Barner KE. ECG signal denoising and baseline wander correction based on the empirical mode decomposition [J]. Computers in Biology & Medicine, 2008, 38(1):1-13.

[6] Kumari R, Sadasivam V. De-noising and baseline wandering removal of electrocardiogram using double density discrete wavelet[J]. International Journal of Wavelets Multi-Resolution and Information Processing, 2011, 5(5):399-415.

[7] Qureshi SA, Masood I, Hashmi M, et al. Noise reduction of electrocardiographic signals using wavelet transforms[J]. Electronics & Electrical Engineering, 2014, 20(3):29-32.

[8] Daubechies I. Ten Lectures on Wavelets[M]. Philadelphia: Society for Industrial and Applied Mathematics, 1992:26-32.

[9] Abdelnour AF, Selesnick IW. Design of 2-band orthogonal near-symmetric CQF[C]//IEEE International Conference on Acoustics, Speech, and Signal Processing. Salt Lake City: IEEE, 2001:3693-3696.

[10] Wu Y, Zhao Y, Shen L. The Analysis and design of nearly-orthogonal symmetric wavelet filter banks[J]. Journal of Computers, 2009, 4(1):199-218.

[11] Murugesan S, Tay D. Design of almost symmetric orthogonal wavelet filter bank via direct optimization[J]. IEEE Transactions on Image Processing, 2012, 21(5):2474-2480.

[12] Murugesan S, Tay D. A new class of almost symmetric orthogonal Hilbert pair of wavelets [J]. Signal Processing, 2014, 95(2):76-87.

[13] Murugesan S, Tay D. New techniques for rationalizing orthogonal and biorthogonal wavelet filter coefficients[J]. IEEE Transactions on Circuits & Systems I, 2012, 59(3):628-637.

Construction of Orthogonal Wavelet for ECG Signal Processing

Zhu Junjiang*Zhang Yuanhui

(SchoolofMechanicalandelectricalengineering,ChinaJilianguniversity,Hangzhou310018,China)

ECG signal; optimization; wavelet; denoising

10.3969/j.issn.0258-8021. 2017. 01.014

2016-01-11，录用日期：2016-09-06

国家自然科学基金(61302191)

R318

D

0258-8021(2017) 01-0109-05

*通信作者(Corresponding author)， E-mail: zjj602@yeah.net