心电信号平滑分解阈值去噪方法

张淼，魏国

(哈尔滨工业大学 电气工程及自动化学院，黑龙江 哈尔滨 150001)

心电图(electrocardiogram, ECG)是人类心脏生理电信号的记录，直接反映心脏的健康状况。对心律失常、室性早搏、心肌缺血和心肌梗死等心脏疾病的诊断具有重要意义。然而，实测心电图总是混有各种噪声：肌电噪声、电极接触噪声、运动伪迹、基线漂移以及电力线干扰等。这些噪声会影响心电图的真实形态，甚至淹没那些重要的心电图细节，导致医生误诊。因此，寻找一个好的心电图降噪方法一直是科学研究者们的奋斗目标。

现有的心电去噪方法主要包括基于傅里叶分析的经典数字滤波、自适应滤波[1-3]、神经网络及支持向量机学习算法[4]、小波去噪算法[5-6]、经验模态分解算法[7-11]等。由于心电信号具有非线性非平稳特性，傅里叶分析不能很好地解决频谱随时间变化的问题。神经网络学习方法则需要大量的实例信号用于训练，而处于不同时间和空间的受试者，其心电信号有很大差异，训练结果的普适性较差。小波分析和经验模态分解为非线性非平稳信号的分析和心电信号去噪提供了有力工具，并成为了目前心电去噪领域的热点。然而，小波函数的选择对小波分析去噪效果有很大影响，目前还没有针对心电信号的小波函数[5]。相对而言，经验模态分解(empirical mode decomposition, EMD)的自适应性，更有利于心电信号的去噪。但是，EMD在分解心电信号时存在模态混叠现象[7, 12]，影响信号分析的成功率。总体经验模态分解(ensemble empirical mode decomposition, EEMD)虽然解决了EMD的模态混叠问题，但是总体数量的增加使得重复分解的计算量成倍增长[13]。而辅助噪声的引入，又产生了噪声残留和信号重构误差的问题。

本文提出平滑分解(smoothing decomposition, SD)方法，用于心电信号去噪，意在避免模态混叠和基本函数的选择，简化分解过程以提高计算效率。SD基于平滑滤波原理，通过改变平滑点数，将含有噪声的心电信号分解为从高频到低频的平滑分解分量(smoothing decomposition component, sdc)。进而建立分界阈值，将sdc分成3组：含高频噪声sdc组、有用信号sdc组和低频噪声sdc组。提出时变阈值(time-varying threshold, TVth)和区间阈值(interval threshold, Ith)概念，对含高频噪声sdc组进行阈值去噪。最后将低频噪声sdc组中的信号作为基线漂移和运动伪迹剔除，用去噪后的含高频噪声sdc组和有用信号sdc组进行信号重构获得去噪心电信号。通过对仿真心电信号和MIT-BIH心律不齐数据库的真实心电信号进行去噪实验，验证了本文方法的有效性，去噪效果良好。

1 平滑分解阈值去噪方法

1.1 平滑分解

根据奈奎斯特定理，采样频率必须大于等于有用信号最高频率的2倍。假设对心电信号的采样满足奈奎斯特采样定理，则实测信号采样频率的1/2为有用心电信号的最高频率。通过三点平滑滤波，可以将频率大于1/2采样频率的信号滤除，将滤除的信号定义为第1阶平滑分解分量(sdc1)。而剩余信号中的最高频率则为前一阶提取频率的下限。继续按剩余信号最高频率的1/2进行新的平滑滤波提取，平滑滤波点数为前一级平滑滤波点数的2倍减1。直到平滑点数大于1/2数据长度且小于数据长度时分解结束。平滑分解方法描述如下。

设实际检测到的心电信号为s(n)，n=1，2，…，N为样本序号，N为信号长度。记第k阶平滑分解分量为sdck(n)，剩余信号记为rk(n)，k=1，2，…，M。M=int(lb2(N-1))为最大分解阶数。

1)当k=1时

(1)

r1(n)=s(n)-sdc1(n)

(2)

2)当k≥2时

(3)

rk(n)=rk-1(n)-sdck(n)

(4)

在以上各式中，当i<1或i>N时，s(i)和rk-1(i)取相应的数据端点值。

3) 信号重构

(5)

1.2 sdc分组

借鉴小波分解和经验模态分解去噪方法的经验，需将sdc分成3组：含高频噪声sdc组、有用信号sdc组和低频噪声sdc组。对于含高频噪声sdc组，主要包含随机噪声、肌电噪声和有用信号的高频成分，采取阈值去噪法进行去噪[9-11, 14]。对于有用信号sdc组，噪声很小可以忽略。而低频噪声sdc组的主要成分为基线漂移和运动伪迹，在心电信号重构时被剔除。

1.2.1 含高频噪声sdc组分界

文献[15]指出，高斯白噪声的固有模态函数(intrinsic mode function, imf)仍然服从高斯分布。imf的能量随阶数k的增加近似按1/2k倍的规律减少，而sdc的特性与imf相似，这一点在后面的实验结果分析中可以得到证实。通常当信号的均值趋于零时，方差与平均能量具有相同的变化规律。因此，用各阶sdc的总体方差与相应阶次中包含的噪声估计方差进行依次比较，当两者的比值达到极大值时，表明sdc中噪声能量占比最小。将噪声能量占比最小的sdc阶数定义为含高频噪声sdc组的分界阈值(high-frequency boundary threshold, HBth)，阶数小于等于HBth的所有sdc为含高频噪声sdc组。实际应用中，当sdc的阶数较高时，对噪声方差的估计会有较大的偏差。因此，必须对参与噪声能量占比计算的sdc做一定的限制。具体的HBth计算方法如下：

(6)

(7)

(8)

(9)

首先计算σdiff(k)的第1个极小值，对应的sdc阶数记为Kdiff，然后在[2,Kdiff]区间内计算σratio(k)的第1个极大值，对应的sdc阶数记为Kratio。含高频噪声sdc组的分界阈值HBth由下式确定:

(10)

1.2.2 低频噪声sdc组分界

根据实测心电信号的R峰个数和位置，计算出平均心率周期。将各阶sdc的平均周期与心电平均周期比较，小于等于心电平均周期的sdc中，最大的sdc阶数定义为低频噪声sdc组的分界阈值(low-frequency boundary threshold, LBth)。阶数大于LBth的所有sdc确定为低频噪声sdc组。

设心电信号的平均周期为PECG，sdc的平均周期为Psdc，分别由(11)和(12)计算。

(11)

(12)

式中：NR为R峰个数；tR(1)和tR(NR)分别为第1个和最后一个R峰对应的时间；Nk为sdck的极大值点或极小值点个数；tPeakk(1)和tPeakk(Nk)分别为sdck的第1个极值点和最后一个极值点对应的时间。低频噪声sdc组的分界阈值LBth为：

(13)

1.3 阈值去噪方法

1.3.1 时变阈值去噪

基于通用阈值去噪方法的思想[16-17]，提出时变阈值(time-varying threshold，TVth)的概念，用于适应肌电噪声统计特性的时变性。由于通用阈值去噪会使噪声大的地方产生噪声残留，噪声小的地方会造成心电信号的损失。因此，需要建立自适应变化的时变阈值进行去噪。时变阈值以心电平均周期为时间窗口，在窗口内按通用阈值法计算阈值，当时间窗连续通过信号长度时，得到时变阈值函数曲线，具体方法如下：

时间窗内信号的噪声方差估计为：

(14)

式中PECG由(11)式确定。

计算时变阈值曲线：

(15)

用时变阈值进行去噪：

(16)

式中：下标k代表第k阶分解分量；系数C由经验确定，本文取为1.5。

图1为通用阈值和时变阈值的比较。其中水平虚线为通用阈值线，加粗实线为时变阈值线，信号为含有肌电噪声的心电信号。可以看出，时变阈值随肌电噪声的大小而自适应变化。

图1 通用阈值和时变阈值的比较Fig.1 Comparison of universal and time-varying threshold

1.3.2 区间阈值去噪

(17)

式中参数2C+1为平滑点数，由经验选取，本文取C=3。

区间阈值计算如下：

(18)

区间阈值法剔除残留噪声：

sdcId,k(n)=sdcTVd,k(n)·Ith(n)

(19)

k=1,2，…，HBth-1

1.4 去噪心电信号的重构

由时变阈值和区间阈值去噪后的含高频噪声sdc组和有用信号sdc组相叠加，重构去噪心电信号:

(20)

1.5 SD阈值去噪方法的计算流程

本文提出的SD阈值去噪方法的去噪流程如图2所示。

图2 本文方法去噪流程Fig.2 Flow chart of the proposed method for denoising

2 去噪结果评估方法

2.1 信噪比提升

信噪比提升定义为去噪后信号的信噪比与去噪前信号的信噪比之差，表明去噪方法的去噪能力和去噪效果。设SNRimp、SNRd、SNRn分别为信噪比提升、去噪后信噪比和去噪前信噪比:

(21)

(22)

SNRimp=SNRd-SNRn

(23)

式中：s(n)为不含噪声的源信号；sn(n)为含噪声的测量信号;sd(n)为去噪后的信号。

2.2 相似度指标

CR定义为去噪后信号与无噪声源信号之间的相关系数，用于表达去噪后信号与源信号的相似程度，从波形形态方面评价去噪效果:

(24)

2.3 幅值度指标

相关系数CR的大小只能反映去噪信号和源信号之间在波形上的相似度。而心电信号的P波、T波、QRS波的大小和位置(间期)则是心脏病理判断的重要指标。为此，用去噪信号与源信号的平均能量之比，作为去噪后心电信号幅值变化的评价指标，用ER表示。当ER等于 1时，表明去噪信号的幅值与源信号一致；ER大于1时，说明去噪信号中还有残留噪声；ER小于1时，说明去噪后信号出现了有用信息的损失。能量比ER的计算式为：

(25)

3 实验验证与实验结果分析

3.1 实验数据来源

1)仿真心电信号。6个成人仿真心电信号Syn00～Syn05，由文献[18]和[19]提供的心电信号工具包(OSET)中的FECGSYN仿真生成。采样频率为360 Hz，信号长度为10 s。如图3所示。需要指出的是，本文中为了后续计算信噪比等参数，对数据均进行了统一化处理，因此纵坐标没有实际的物理含义。通过图展示的是信号形状，因而不关心振幅的单位。本文后续的图均为此情况。

2)真实心电信号。真实心电信号来源于MIT-BIH心律不齐数据库(MITDB)[20]。采样频率360 Hz，截取的时间长度为10 s。表1给出本文使用的数据记录编号、通道名和开始时间位置。

3)实验中的噪声。本文实验涉及了4种噪声，分别为基线漂移(baseline wandering, BW)、肌电噪声(electromyogram, EMG)、白噪声(white noise, WN)和混合噪声(hybrid noise, HN)。其中BW和EMG来源于MIT-BIH的NSTDB数据库[21]，采样频率360 Hz，截取长度10 s。白噪声WN由Matlab 2016的随机函数生成。而混合噪声HN则由以上3种噪声混合生成。其中BW经过高频噪声滤除，而EMG则经过基线漂移修正，如图4所示。

图3 仿真心电信号Fig.3 Simulated ECG signals

表1 从MIT-BIH心律不齐数据库选取的真实心电信号

3.2 实验结果及分析

3.2.1 平滑分解方法的验证

随机采集了30个单位方差的白噪声信号，采样频率360 Hz，采样时间10 s。分别用本文方法进行平滑分解。针对每个分解分量sdc计算峰值点个数、平均频率以及sdc能量与原始信号能量的百分比。平均结果列于表2。图5给出了某个白噪声信号的分解结果。结果显示sdc的峰值点数、平均频率和能量百分比都与文献[18]中的EMD分解结果一致。因此平滑分解(SD)方法具有与EMD相似的性能。

图4 本文实验中的4种噪声Fig.4 Four categories of noises used in experiments.

表2 30个随机白噪声平滑分解结果的平均值

3.2.2 时变阈值和区间阈值去噪效果

仿真信号Syn00叠加EMG噪声，改变信噪比从-5～30 dB。用本文的SD方法进行分解后，分别用通用阈值法、时变阈值法、区间阈值法和时变+区间阈值法进行去噪实验。结果列于表3～5。表中数据显示，时变阈值法去噪效果好于通用阈值法，时变阈值+区间阈值法最好。

表3 Syn00+EMG的去噪结果(信噪比提升)Table 3 Denoising results of Syn00+EMG(SNRimp)

图5 白噪声的平滑分解分量(前5 s，1 800采样点)Fig.5 Smooth decomposition component of white noise (first 5 s, 1 800 sampling points)

表4 Syn00+EMG的去噪结果(相关系数)Table 4 Denoising results of Syn00+EMG(CR)

表5 Syn00+EMG的去噪结果(能量比)Table 5 Denoising results of Syn00+EMG(ER)

图6为Syn00+EMG信号在信噪比5 dB时，应用不同阈值法得到的去噪信号波形。结果显示，通用阈值法和时变阈值法，在EMG变化较大的地方，存在噪声残留。区间阈值法则在Q和S波附近有噪声残留。这是因为区间阈值不能去除波群范围内的噪声所致。而时变+区间阈值法综合了两者优势，能够有效弥补其不足，因此去噪效果最好，既没有残留噪声，又能很好地保留信号波形的原有特征，波形失真度最小，适合于肌电噪声的消除。

3.2.3 与小波去噪和EMD去噪结果的比较

选择仿真信号Syn00和Syn03分别叠加5 dB的肌电噪声和白噪声进行小波阈值法和EMD阈值法实验。实验结果与本文方法进行比较，验证本文方法的有效性和效果。小波阈值去噪法选择了OriginPro 8 SR4的wtdenoise方法，以避免编程对去噪效果的影响。阈值选择为sqtwolog固定阈值。针对不同的DB小波进行实验，实验结果的局部波形如图7所示，从中不难看出不同DB小波去噪结果的不同。

图6 4种阈值去噪方法对Syn00+EMG的去噪结果(5 dB)Fig.6 Denoising results of Syn00+EMG by four threshold denoising methods (5 dB)

图7 不同小波函数对小波去噪结果的影响Fig.7 Wavelet denoising resuls by different wavelet functions

EMD阈值去噪分2种情况：原始EMD通用阈值去噪和EEMD通用阈值去噪。表6为EMD、EEMD和本文的SD方法的分解时间和循环次数比较。表中信号采样频率360 Hz，信号长度3 600点。EEMD的总体数量取为50，循环次数可由EMD循环次数乘以50估算。表6结果显示，本文SD方法的分解时间远远小于EMD和EEMD。

表6 含噪心电信号分解时间比较Table 6 Comparison of decomposition time for noisy ECG

图8为DB4小波阈值去噪、EMD通用阈值去噪、EEMD通用阈值去噪和本文SD通用阈值去噪及SD时变+区间阈值去噪的比较。图中左侧为Syn00叠加肌电噪声的结果，右侧为Syn00叠加白噪声的结果。可以看出，本文方法去噪结果好于其他分解方法，特别是时变+区间阈值方法效果更好。从去噪波形上看，EMD、EEMD和WT更适合于白噪声的去噪，而本文方法对肌电噪声去噪也适合。

3.2.4 仿真心电信号实验结果及分析

通过以上实验，验证了本文方法具有很好的去噪效果，特别是对肌电噪声的滤除具有明显优势。本节将通过仿真心电信号叠加不同信噪比的各种噪声进行去噪实验，进一步验证本文方法的有效性、适应性和去噪效果。

表7～10分别给出了基线漂移(BW)、肌电噪声(EMG)、白噪声(WN)和混合噪声(HN)与仿真信号Syn00～Syn05叠加后，用本文方法进行去噪，各信号去噪结果的平均值。表中数据表明：在信噪比从0 dB～20 dB范围内，基线漂移滤除效果最好，相关系数大于0.99，能量比在1～1.02，表明没有信号损失，噪声残留较小。肌电噪声去噪效果相对最差，但是信噪比提升大于2.44 dB，相关系数大于0.88，能量比大于0.94小于1，这说明去噪时有一定的信号损失，但损失较小，不影响结果的合理性，仍具有很好的去噪效果。白噪声去噪结果的信噪比提升和相关系数均高于肌电噪声的情况。但是，能量比小于1，并且小于肌电噪声去噪的能量比。这说明白噪声去噪时，产生了较大的幅值损失，但是相关系数大于0.954，具有优良的波形相似性。混合噪声的去噪效果介于肌电噪声和白噪声之间，这是因为多个随机变量的叠加，统计特性趋于高斯分布的原因。

表7 仿真心电信号基线漂移(BW)滤除结果Table 7 Results of eliminating BW for simulated ECG

表8 仿真心电信号肌电噪声(EMG)去噪结果Table 8 Results of eliminating EMG for simulated ECG

表9 仿真心电信号白噪声(WN)去噪结果Table 9 Results of eliminating WN for simulated ECG

图8 本文方法SD与WT、EMD、EEMD去噪波形比较(5 dB)Fig.8 Comparison of denoising waveforms among the proposed SD, WT, EMD and EEMD (5 dB)

表10 仿真心电信号混合噪声(HN)去噪结果Table 10 Results of eliminating HN for simulated ECG

图9给出仿真信号Syn00～Syn05叠加5 dB肌电噪声后，去噪信号与未去噪信号的波形比较。图10为Syn04分别与BW、HN、EMG和WN 4种噪声叠加的去噪波形。由图9和图10可以看出，去噪后波形与无噪信号(Syn00～Syn05)波形一致，P波、T波和QRS波群完整清晰。

3.2.5 真实心电信号实验结果及分析

真实心电信号来源于MIT-BIH心律不齐数据库的临床记录，选择其中的10个记录，如表1所示。本文将其命名为MB100、MB101、MB103、MB112、MB113、MB115、MB117、MB119、MB122、MB123。由于这些临床数据不可避免地含有一定的噪声，首先对其进行高频噪声滤除和基线漂移修正，目的是为了在分析去噪效果时能够定量分析。经过去噪和修正的真实心电信号记为DMB100、DMB101、DMB103、DMB112、DMB113、DMB115、DMB117、DMB119、DMB122、DMB123。以DMB作为干净信号，估计原信号MB的信噪比，结果列于表11。图11选择给出了3个信号的噪声滤除前后的波形和滤除的噪声曲线。由表11数据和图11的波形可以证明，DMB能够反映真实的心电特性。

以DMB作为无噪真实心电信号，叠加BW、EMG、WN和HN建立含噪信号。然后应用本文方法进行平滑分解和时变+区间阈值去噪。去噪结果按10个DMB信号进行平均，结果给出在表12～15中。在信噪比0 dB～20 dB范围内，真实信号与仿真信号去噪结果一致。真实心电信号去噪结果中，基线漂移去噪效果最好，信噪比提升大于10 dB，相关系数大于0.99，表明去噪信号与源信号的相似度达到99%，能量比大于1表明去噪信号的能量大于源信号，去噪信号未丢失有用信息。能量比小于1.017表明去噪信号中残留的噪声成分不大于1.7%。肌电噪声去噪效果相对差一些，信噪比提升相对较少，相关系数大于0.89，能量比大于0.93。白噪声去噪结果在信噪比提升和相关系数上优于肌电噪声。但是，能量比大于0.87小于1，比肌电噪声的能量比低，因此有更大的幅值改变。混合噪声的去噪结果和白噪声相当。

图9 仿真心电信号的肌电噪声去噪结果(5 dB)Fig.9 Waveforms after eliminating EMG for simulated ECG (5 dB)

图10 Syn04分别叠加4种噪声时的去噪结果(5 dB)Fig.10 Waveforms after eliminating four categories of noise for Syn04 (5 dB)

图11 MB100、MB115和MB122的噪声滤除结果Fig.11 Noise filtering results of MB100, MB115 and MB122

表11 MIT-BIH真实心电信号的信噪比估计Table 11 SNR estimation for MIT-BIH real ECG

表12 真实心电信号基线漂移(BW)去噪结果Table 12 Results of eliminating BW for real ECG

表13 真实心电信号肌电噪声(EMG)去噪结果Table 13 Results of eliminating EMG for real ECG

表14 真实心电信号白噪声(WN)去噪结果Table 14 Results of eliminating WN for real ECG

表15 真实心电信号混合噪声(HN)去噪结果Table 15 Results of eliminating HN for real ECG

图12为真实心电信号叠加5 dB肌电噪声时的去噪信号波形。图13选择了真实心电信号DMB119与BW、HN、EMG和WN 4种噪声叠加时，本文方法的去噪波形。可以看出，真实心电信号去噪后的P波、T波和QRS波群还原良好。证明本文方法不仅适合于肌电噪声的去噪，而且对其他噪声的去噪效果更好。

图12 真实心电信号的肌电噪声去噪结果(5 dB)Fig.12 Waveforms after eliminating EMG for real ECG (5 dB)

图13 DMB119分别叠加4种噪声时的去噪结果(5 dB)Fig.13 Waveforms after eliminating four categories of noise for DMB119 (5 dB)

4 结论

1)通过噪声分解分析，证明了本文提出的平滑分解方法与EMD分解具有相似的滤波组特性，并具有更高的分解效率；通过对仿真心电信号和MIT-BIH真实心电信号的去噪实验表明，本文方法的sdc分组、时变阈值和区间阈值，解决了肌电噪声去噪的难题，去噪效果明显好于目前的EMD阈值去噪、EEMD阈值去噪和小波阈值去噪效果；在实验范围内，本文的平滑分解阈值去噪方法具有较高的信噪比提升，相关系数总体上大于0.9以上，去噪后的心电波形失真小，P波、T波和QRS复合波完整清晰，其位置、形态和幅值与源信号一致，达到实际临床应用的水平。

2)本文方法基于所提出的时变阈值和区间阈值方法，很好地解决了心电信号中肌电噪声问题。平滑滤波分解大大提高了EMD和EEMD的分解效率。避免了小波分解中小波函数选择的不确定性。

3)本文方法没有给出心电周期的理论确定方法，尚需观察给出。对于包含在有用信号sdc中的中低频接触噪声这一目前各种去噪方法均未涉及的难题，没有进行研究，有望于进一步研究解决。