工程科学学报 Chinese Journal of Engineering 陡脉冲干扰下的心电信号滤波及QRS提取 姚晰童代煜张建勋葛锦涛陈通杨灏 ECG filtering and QRS extraction under steep pulse interference YAO Xi-tong.DAI Yu.ZHANG Jian-xun.GE Jin-tao,CHEN Tong.YANG Hao 引用本文： 姚晰童，代煜，张建勋，葛锦涛，陈通，杨灏.陡脉冲干扰下的心电信号滤波及QRS提取U.工程科学学报，2020,42(5：654- 662.doi:10.13374/i.issn2095-9389.2019.06.20.004 YAO Xi-tong,DAI Yu,ZHANG Jian-xun,GE Jin-tao,CHEN Tong,YANG Hao.ECG filtering and QRS extraction under steep pulse interference[J].Chinese Journal of Engineering,2020,42(5):654-662.doi:10.13374/j.issn2095-9389.2019.06.20.004 在线阅读View online::htps:/ldoi.org10.13374.issn2095-9389.2019.06.20.004 您可能感兴趣的其他文章 Articles you may be interested in 隧道地质预报探地雷达信号干扰消除方法 Research on the interference elimination method of GPR signal for tunnel geological prediction 工程科学学报.2020.42(3：390 https:doi.org10.13374.issn2095-9389.2019.04.12.001 基于光电容积脉搏波的有限空间生理疲劳测量 Confined space physiological fatigue measurement based on photoplethysmography pulse wave signal 工程科学学报.2018.40(10：1215htps:ldoi.org10.13374.issn2095-9389.2018.10.008 基于高阶同步压缩变换的行星齿轮箱声音信号共振频带特征提取 Acoustic Signal Analysis in Resonance Region for Planetary Gearbox Fault Diagnosis Based on High-Order Synchrosqueezing Transform 工程科学学报.优先发表htps:/1doi.org/10.13374j.issn2095-9389.2019.07.18.002 降压式18脉冲自耦变压器优化设计 Optimal design of a new step-down 18-pulse autotransformer 工程科学学报.2017,393：456 https::/1oi.org/10.13374.issn2095-9389.2017.03.019 基于GPR反射波信号多维分析的隧道病害智能辨识 An intelligent identification method to detect tunnel defects based on the multidimensional analysis of GPR reflections 工程科学学报.2018,40(3：293 https:loi.org10.13374j.issn2095-9389.2018.03.005 基于管道流体信号的自振射流特性检测方法 Detection method of the self-resonating waterjet characteristic based on the flow signal in a pipeline 工程科学学报.2019,41(3：377 https:1doi.org/10.13374.issn2095-9389.2019.03.011
陡脉冲干扰下的心电信号滤波及QRS提取 姚晰童 代煜 张建勋 葛锦涛 陈通 杨灏 ECG filtering and QRS extraction under steep pulse interference YAO Xi-tong, DAI Yu, ZHANG Jian-xun, GE Jin-tao, CHEN Tong, YANG Hao 引用本文: 姚晰童, 代煜, 张建勋, 葛锦涛, 陈通, 杨灏. 陡脉冲干扰下的心电信号滤波及QRS提取[J]. 工程科学学报, 2020, 42(5): 654- 662. doi: 10.13374/j.issn2095-9389.2019.06.20.004 YAO Xi-tong, DAI Yu, ZHANG Jian-xun, GE Jin-tao, CHEN Tong, YANG Hao. ECG filtering and QRS extraction under steep pulse interference[J]. Chinese Journal of Engineering, 2020, 42(5): 654-662. doi: 10.13374/j.issn2095-9389.2019.06.20.004 在线阅读 View online: https://doi.org/10.13374/j.issn2095-9389.2019.06.20.004 您可能感兴趣的其他文章 Articles you may be interested in 隧道地质预报探地雷达信号干扰消除方法 Research on the interference elimination method of GPR signal for tunnel geological prediction 工程科学学报. 2020, 42(3): 390 https://doi.org/10.13374/j.issn2095-9389.2019.04.12.001 基于光电容积脉搏波的有限空间生理疲劳测量 Confined space physiological fatigue measurement based on photoplethysmography pulse wave signal 工程科学学报. 2018, 40(10): 1215 https://doi.org/10.13374/j.issn2095-9389.2018.10.008 基于高阶同步压缩变换的行星齿轮箱声音信号共振频带特征提取 Acoustic Signal Analysis in Resonance Region for Planetary Gearbox Fault Diagnosis Based on High-Order Synchrosqueezing Transform 工程科学学报.优先发表 https://doi.org/10.13374/j.issn2095-9389.2019.07.18.002 降压式18脉冲自耦变压器优化设计 Optimal design of a new step-down 18-pulse autotransformer 工程科学学报. 2017, 39(3): 456 https://doi.org/10.13374/j.issn2095-9389.2017.03.019 基于GPR反射波信号多维分析的隧道病害智能辨识 An intelligent identification method to detect tunnel defects based on the multidimensional analysis of GPR reflections 工程科学学报. 2018, 40(3): 293 https://doi.org/10.13374/j.issn2095-9389.2018.03.005 基于管道流体信号的自振射流特性检测方法 Detection method of the self-resonating waterjet characteristic based on the flow signal in a pipeline 工程科学学报. 2019, 41(3): 377 https://doi.org/10.13374/j.issn2095-9389.2019.03.011
工程科学学报.第42卷.第5期：654-662.2020年5月 Chinese Journal of Engineering,Vol.42,No.5:654-662,May 2020 https://doi.org/10.13374/j.issn2095-9389.2019.06.20.004;http://cje.ustb.edu.cn 陡脉冲干扰下的心电信号滤波及QRS提取 姚晰童)，代煜区，张建勋，葛锦涛，陈通，杨灏2) 1)南开大学机器人与信息自动化研究所，天津3003102)南开大学电子信息与光学工程学院，天津300310 通信作者，E-mail:daiyu@nankai.edu.cn 摘要为消除陡脉冲带来的干扰，分析了陡脉冲干扰的特点，建立了陡脉冲噪声数学模型，提出了基于变分模态分解 (Variational mode decomposition,VMD)的心电信号滤波算法，提取叠加在心电信号中陡脉冲干扰分量、识别陡脉冲干扰分量 并剔除陡脉冲干扰分量：为减少VMD分解层数、提高实时性并减少内存消耗，提出了心电信号预处理算法；针对医疗环境中 的随机噪声伴随陡脉冲出现的情况，分析了VMD后子信号中随机噪声的特点，提出了基于VMD子信号能量估计的阈值去 噪算法；利用变分模态分解的带通滤波器组特性，提出了基于变分模态分解子信号重组的Q$波群检测算法，配合滤波算法 以提高心电信号特征检测精度.以添加了高斯白噪声和模拟陡脉冲干扰的MT-BH数据库心电信号和医疗环境中采集的心 电信号为实验对象，分别实现对滤波算法和QS波群检测算法的定量对比分析. 关键词心电信号：陡脉冲干扰：变分模态分解：MT-BH数据库：QRS波群 分类号TN911.72 ECG filtering and QRS extraction under steep pulse interference YAO Xi-tong,DAI Yu,ZHANG Jian-xun,GE Jin-tao,CHEN Tong,YANG Hao2) 1)Institute of Robotics Automatic Information System,Nankai University,Tianjin 300310,China 2)College of Electronic Information and Optical Engineering.Nankai University,Tianjin 300310,China Corresponding author,E-mail:daiyu@nankai.edu.cn ABSTRACT Applying a steep pulse voltage of appropriate amplitude to a cell membrane can induce transient and reversible breakdown of the membrane,which has broad application prospects in biomedicine and clinical fields.However,the noise generated by the steep pulse seriously interferes with a patient's electrocardiogram(ECG)signal resulting in decrease in the accuracy of the ECG feature point detection algorithm.Thus,doctors are unable to understand the state of the patient during treatment,thus limiting complete benefits of the therapy.To eliminate the interference caused by the steep pulse,we analyzed the characteristics of steep pulse interference and established the mathematical model of steep pulse noise.Moreover,we proposed an ECG signal filtering algorithm based on variational mode decomposition(VMD)to extract the steep pulse interference component superimposed on the ECG signal. The proposed algorithm could identify and eliminate the steep pulse interference component.We also designed an ECG signal preprocessing algorithm to reduce the decomposition layer of the VMD algorithm,which improved the real-time performance and reduced the memory consumption.To identify the random noise in the medical environment accompanied by the occurrence of steep pulses,we analyzed the characteristics of random noise in the sub-signal after VMD.Further,we proposed a threshold denoising algorithm based on VMD for sub-signal energy estimation.On the basis of the characteristics of a band-pass filter bank with VMD,we proposed a QRS complex detection algorithm based on VMD sub-signal recombination.Combined with the filtering algorithm,the proposed algorithm was able to improve the accuracy of ECG signal detection.By conducting experiments on ECG signals from the MIT-BIH database with Gaussian white noise and simulated steep pulse interference and those collected in the medical environment,we 收稿日期：2019-06-20 基金项目：国家重点研发计划资助项目(2017YFC0110402:天津市自然科学基金项目资助项目(18 JCYBJC18800)
陡脉冲干扰下的心电信号滤波及 QRS 提取 姚晰童1)，代 煜1) 苣，张建勋1)，葛锦涛1)，陈 通1)，杨 灏2) 1) 南开大学机器人与信息自动化研究所，天津 300310 2) 南开大学电子信息与光学工程学院，天津 300310 苣通信作者，E-mail：daiyu@nankai.edu.cn 摘 要 为消除陡脉冲带来的干扰，分析了陡脉冲干扰的特点，建立了陡脉冲噪声数学模型，提出了基于变分模态分解 （Variational mode decomposition, VMD）的心电信号滤波算法，提取叠加在心电信号中陡脉冲干扰分量、识别陡脉冲干扰分量 并剔除陡脉冲干扰分量；为减少 VMD 分解层数、提高实时性并减少内存消耗，提出了心电信号预处理算法；针对医疗环境中 的随机噪声伴随陡脉冲出现的情况，分析了 VMD 后子信号中随机噪声的特点，提出了基于 VMD 子信号能量估计的阈值去 噪算法；利用变分模态分解的带通滤波器组特性，提出了基于变分模态分解子信号重组的 QRS 波群检测算法，配合滤波算法 以提高心电信号特征检测精度. 以添加了高斯白噪声和模拟陡脉冲干扰的 MIT−BIH 数据库心电信号和医疗环境中采集的心 电信号为实验对象，分别实现对滤波算法和 QRS 波群检测算法的定量对比分析. 关键词 心电信号；陡脉冲干扰；变分模态分解；MIT−BIH 数据库；QRS 波群 分类号 TN911.72 ECG filtering and QRS extraction under steep pulse interference YAO Xi-tong1) ，DAI Yu1) 苣 ，ZHANG Jian-xun1) ，GE Jin-tao1) ，CHEN Tong1) ，YANG Hao2) 1) Institute of Robotics & Automatic Information System, Nankai University, Tianjin 300310, China 2) College of Electronic Information and Optical Engineering, Nankai University, Tianjin 300310, China 苣 Corresponding author, E-mail: daiyu@nankai.edu.cn ABSTRACT Applying a steep pulse voltage of appropriate amplitude to a cell membrane can induce transient and reversible breakdown of the membrane, which has broad application prospects in biomedicine and clinical fields. However, the noise generated by the steep pulse seriously interferes with a patient ’s electrocardiogram (ECG) signal resulting in decrease in the accuracy of the ECG feature point detection algorithm. Thus, doctors are unable to understand the state of the patient during treatment, thus limiting complete benefits of the therapy. To eliminate the interference caused by the steep pulse, we analyzed the characteristics of steep pulse interference and established the mathematical model of steep pulse noise. Moreover, we proposed an ECG signal filtering algorithm based on variational mode decomposition (VMD) to extract the steep pulse interference component superimposed on the ECG signal. The proposed algorithm could identify and eliminate the steep pulse interference component. We also designed an ECG signal preprocessing algorithm to reduce the decomposition layer of the VMD algorithm, which improved the real-time performance and reduced the memory consumption. To identify the random noise in the medical environment accompanied by the occurrence of steep pulses, we analyzed the characteristics of random noise in the sub-signal after VMD. Further, we proposed a threshold denoising algorithm based on VMD for sub-signal energy estimation. On the basis of the characteristics of a band-pass filter bank with VMD, we proposed a QRS complex detection algorithm based on VMD sub-signal recombination. Combined with the filtering algorithm, the proposed algorithm was able to improve the accuracy of ECG signal detection. By conducting experiments on ECG signals from the MIT–BIH database with Gaussian white noise and simulated steep pulse interference and those collected in the medical environment, we 收稿日期: 2019−06−20 基金项目: 国家重点研发计划资助项目 (2017YFC0110402)；天津市自然科学基金项目资助项目 (18JCYBJC18800) 工程科学学报，第 42 卷，第 5 期：654−662，2020 年 5 月 Chinese Journal of Engineering, Vol. 42, No. 5: 654−662, May 2020 https://doi.org/10.13374/j.issn2095-9389.2019.06.20.004; http://cje.ustb.edu.cn
姚晰童等：陡脉冲干扰下的心电信号滤波及QRS提取 655· compared and analyzed the filtering algorithm and QRS complex detection algorithm. KEY WORDS ECG signal;steep pulse interference;variational mode decomposition;MIT-BIH database;QRS complex 陡脉冲是指短时高压电脉冲.通过向细胞施 充电过程：后支以前支的终点为起点，频率低，时 加陡脉冲可以使得细胞出现瞬时的电穿孔.基于 间长，对心电信号的波形影响较大，其数学模型如 这种特性，陡脉冲广泛的应用于靶向药物传递川、 式(2)： 癌症治疗等多个领域，其可行性已经在人体实验 Urall(t)=se- (2) 中得到了证实)除此之外，陡脉冲在材料领域、 式中，Ua为陡脉冲冲激下降沿的幅值，s为陡脉 食品工业问也有着广泛的应用 冲冲激的峰值电压，R、C分别表示人体等效电阻 心电是心脏活动的重要表现，是医生对患者 和等效电容，o为陡脉冲冲激后支起点时间.此 生理状态进行判断的重要佐证.然而，心电信号 外，在医疗环境中，随机噪声伴随陡脉冲干扰心电 (Electrocardiogram,ECG)受到陡脉冲的干扰，波形 信号.随机噪声的能量均匀的分布在整个频带上， 可能发生变异，这导致病患的生理状况无法被掌 在此以nw表示.经过以上分析可知，受陡脉冲相 握，影响治疗效果.以不可逆电消融疗法阿为例， 关噪声干扰的心电信号的数学模型如式(3)： 心电信号的作用除监测患者生理状态外，还包括 Xreal(t)=xideal(f)+npluse(t)+nw(t) (3) 同步陡脉冲的施放，如果不将陡脉冲的施放时间 与心电周期同步，手术范围将被严格的限制在远 其中，xrea与Xideal分别表示实际心电信号和理想心 离心脏的20cm的范围外 电信号，npluse表示陡脉冲干扰 目前，对心电信号的滤波算法的研究，主要针 (a) 对高频噪声、基线漂移、工频干扰，对于陡脉冲干 扰的滤波算法研究较为缺乏.由于陡脉冲产生的 相关噪声与心电信号主成分频带重合，经典的数 字滤波器门虽然简单易实现，但是效果较差：小波 算法图对心电信号滤波是基于信号与噪声频谱分 离特性，运用于此场景效果不理想，且小波基的选 (b) 择影响最终结果；基于EMD的滤波算法9改进了 1 小波算法的缺陷，但是对分解信号的进一步处理 0 '*rtr* 需要依赖经验，且对采样频率和噪声敏感.因此， 0 12 3 456 本文针对陡脉冲干扰的特性，基于VMD1o算法提 Time/s 出了解决方案，实现了对陡脉冲干扰下心电信号 图1受陡脉冲干扰的ECG.(a)受到陡脉冲干扰的心电信号：(b)未 的滤波并满足了心电周期同步需求 受陡脉冲干扰的心电信号 Fig.I ECG disturbed by steep pulses:(a)ECG disturbed by steep 1陡脉冲干扰下心电信号的分析 pulses;(b)ECG not disturbed by steep pulses 陡脉冲具有如下特征：脉冲幅值要求达到 2基于VMD的陡脉冲干扰抑制 kVcm,脉冲上升沿时间为ns级别，脉冲宽度为 s级别.可用表达式(1)表示： VMD算法最早由Dragomiretskiy等于20l4年提 ∫S,h<t≤h+4 出，可以将任意信号x()分解成离散数量K的子信 6 pluse(t)）= (1) 0, others u(Band-limited intrinsic mode functions,BLIMF). 式中，6puse()表示高压陡脉冲的时域曲线，S为陡 每一个子信号的能量都集中在相应的中心频率wk 脉冲幅值，4为陡脉冲脉宽，h为陡脉冲施放时间 周围，具有特殊的稀疏特性其基本表达式如式 点.陡脉冲对ECG的干扰如图1所示.时域上可 (4).这是一种后验的、非递归算法，并具有带通滤 将陡脉冲干扰分为两个部分，即：前支和后支.其 波器属性.目前，这种算法被广泛的应用在轴承 中，前支频率高，时间短，位于心电信号中的不应 故障检测，地震信号分析，经济走势预测等 期，相当于陡脉冲对心电采集设备等效RC电路的 领域
compared and analyzed the filtering algorithm and QRS complex detection algorithm. KEY WORDS ECG signal；steep pulse interference；variational mode decomposition；MIT−BIH database；QRS complex 陡脉冲是指短时高压电脉冲. 通过向细胞施 加陡脉冲可以使得细胞出现瞬时的电穿孔. 基于 这种特性，陡脉冲广泛的应用于靶向药物传递[1]、 癌症治疗[2] 等多个领域，其可行性已经在人体实验 中得到了证实[3] . 除此之外，陡脉冲在材料领域[4]、 食品工业[5] 也有着广泛的应用. 心电是心脏活动的重要表现，是医生对患者 生理状态进行判断的重要佐证. 然而，心电信号 （Electrocardiogram, ECG）受到陡脉冲的干扰，波形 可能发生变异，这导致病患的生理状况无法被掌 握，影响治疗效果. 以不可逆电消融疗法[6] 为例， 心电信号的作用除监测患者生理状态外，还包括 同步陡脉冲的施放，如果不将陡脉冲的施放时间 与心电周期同步，手术范围将被严格的限制在远 离心脏的 20 cm 的范围外. 目前，对心电信号的滤波算法的研究，主要针 对高频噪声、基线漂移、工频干扰，对于陡脉冲干 扰的滤波算法研究较为缺乏. 由于陡脉冲产生的 相关噪声与心电信号主成分频带重合，经典的数 字滤波器[7] 虽然简单易实现，但是效果较差；小波 算法[8] 对心电信号滤波是基于信号与噪声频谱分 离特性，运用于此场景效果不理想，且小波基的选 择影响最终结果；基于 EMD 的滤波算法[9] 改进了 小波算法的缺陷，但是对分解信号的进一步处理 需要依赖经验，且对采样频率和噪声敏感. 因此， 本文针对陡脉冲干扰的特性，基于 VMD[10] 算法提 出了解决方案，实现了对陡脉冲干扰下心电信号 的滤波并满足了心电周期同步需求. 1 陡脉冲干扰下心电信号的分析 陡脉冲具有如下特征 ：脉冲幅值要求达 到 kV·cm−1，脉冲上升沿时间为 ns 级别，脉冲宽度为 μs 级别[11] . 可用表达式（1）表示： δpluse(t) = { S, h < t ⩽ h+∆ 0, others （1） 式中， δpluse (t) 表示高压陡脉冲的时域曲线，S 为陡 脉冲幅值，Δ 为陡脉冲脉宽，h 为陡脉冲施放时间 点. 陡脉冲对 ECG 的干扰如图 1 所示. 时域上可 将陡脉冲干扰分为两个部分，即：前支和后支. 其 中，前支频率高，时间短，位于心电信号中的不应 期，相当于陡脉冲对心电采集设备等效 RC 电路的 充电过程；后支以前支的终点为起点，频率低，时 间长，对心电信号的波形影响较大，其数学模型如 式（2）： Ufall(t) = se − t−t0 RC （2） 式中，Ufall 为陡脉冲冲激下降沿的幅值，s 为陡脉 冲冲激的峰值电压，R、C 分别表示人体等效电阻 和等效电容，t0 为陡脉冲冲激后支起点时间. 此 外，在医疗环境中，随机噪声伴随陡脉冲干扰心电 信号. 随机噪声的能量均匀的分布在整个频带上， 在此以 nw 表示. 经过以上分析可知，受陡脉冲相 关噪声干扰的心电信号的数学模型如式（3）： xreal(t) = xideal(t)+npluse(t)+nw(t) （3） 其中， xreal 与xideal 分别表示实际心电信号和理想心 电信号，npluse 表示陡脉冲干扰. 2 基于 VMD 的陡脉冲干扰抑制 uk wk VMD 算法最早由Dragomiretskiy 等于2014 年提 出[10] ，可以将任意信号 x(t) 分解成离散数量 K 的子信 号 （Band-limited intrinsic mode functions, BLIMF）. 每一个子信号的能量都集中在相应的中心频率 周围，具有特殊的稀疏特性[10] . 其基本表达式如式 （4）. 这是一种后验的、非递归算法，并具有带通滤 波器属性[12] . 目前，这种算法被广泛的应用在轴承 故障检测[13] ，地震信号分析[14] ，经济走势预测[15] 等 领域. Time/s Amplitude/mV 0 1 −1 0 1 2 2 3 4 5 6 (a) Amplitude/mV −1 0 1 2 Time/s 0 1 2 3 4 5 6 (b) 图 1 受陡脉冲干扰的 ECG. （a）受到陡脉冲干扰的心电信号；（b）未 受陡脉冲干扰的心电信号 Fig.1 ECG disturbed by steep pulses: (a) ECG disturbed by steep pulses; (b) ECG not disturbed by steep pulses 姚晰童等： 陡脉冲干扰下的心电信号滤波及 QRS 提取 · 655 ·
656 工程科学学报，第42卷，第5期 (a) min{∑1la,[(60+)*k(0]e-jw》 2 (uk).(wel πt (4) ∑鉴14)=0 10 1214 (b) 式中，K表示分割层数，表示第k个子信号，wk表 示，的中心频率，x表示原信号，6()为狄拉克函数， a,代表梯度运算，*代表卷积运算，j为虚数符号 1214 (c) 2.1心电信号的预处理 根据香农采样定律，信号的频带宽度是采样 8101214 频率的0.5倍.以500Hz的采样频率为例，实际信 号的宽度为0~250Hz.而心电信号的主要成分集 0.2r -4 中于0.05~49HzI6因此，不对50~250Hz的信 -0.2 0 2468101214 号进行分析，以提高实时性.在此，采用了一个 (e) 0.2 4阶R滤波器对这部分信号进行衰减，截止频率 4+44件 为50Hz. -0.2 0 2468101214 2.2陡脉冲干扰的抑制 0.2 对式(2)进行傅里叶变换，可得其幅频曲线 -g4++H+ A(0,如式(5) 0246 8101214 5e是 (g) 0.2 A(f)= (5) 1、2 0 V(RC)+2xf)2 -0.2L 0 6 8 101214 Time/s 式中，∫表示频率.由于心电采样设备通常要求大 图3VMD6阶分解子信号.(a)原信号：(b):(c)2:(d):(e): 输入电阻7，因此1/RC为常数且较小，在此为方 (f)l5:(g)46 便分析将其忽略.幅频特性如图2. Fig.3 Use of VMD to divide the signal into six layers:(a)original signal;(b):(c)2:(d);(e)4,(日s,(g)6 1.5 设心电采集系统共模抑制比为-100dB,陡脉冲 的峰值大于1000V,理论上产生噪声的幅值约大于 0.5 10mV.实际测得幅值由于损耗而略小.以肢体导 联为例，P波幅值不超过0.25mV,T波幅值不超过 0 5101520253035404550 0.5mV16,远低于同频段陡脉冲干扰幅值.基于以上 Frequency/Hz 分析，本文设计了一种自适应阈值，对陡脉冲干扰在2 图2衰减模型单边幅颜图 中的陡脉冲千扰分量进行消除.阈值入s知如式(6). Fig.2 Diagram of a single-sided amplitude-frequency attenuation model Asp mean(sort(lu2D)[1 10]) (6) 可以看到，陡脉冲干扰的后支能量主要分布 上式表示对子信号2中绝对幅值最大的10个采样 在10Hz前.VMD算法基于对混合信号的主成分 点sort2lD[1:10取均值.sort(0表示对输入信号序 估计，结合其带通滤波器属性，包含陡脉冲干扰后 列的绝对值降序排列，mean(~)表示均值运算 支的子信号对应中心频率w应小于10Hz.由于此 u2(乙)，ku2(rl<sp i2()= (7) 部分干扰对心电波形干扰较大，是抑制的重点， 0. lu2(r川≥dsp 图3展示了当=6时，各子信号时域图.其中，子 式(7)为阈值使用方法.式中，为采样点区间，包 信号u的中心频率w接近于0，包含直流分量、基 含相邻两个零交叉点或两个近似零交叉点（所谓 线漂移、陡脉冲干扰的低频干扰分量，将抛弃， 近似零交叉点，即距离实际零交叉点最近的采样 可以有效地消除基线漂移和直流分量.2的中心 点)之间的全部采样点.r代表了之之间的峰值点 频率满足w2<10Hz,包含陡脉冲干扰分量、P波、 其中，将峰值大于阈值的的起点记作zom.处理前 T波分量6.完全抛弃2不利于恢复心电信号. 后效果如图4所示，其结果用2表示
min {uk },{wk } { ∑K k=1 ||∂t[(δ(t)+ j πt ) ∗ uk(t)]e−jwk t ||2 2 } ∑K k=1 uk(t) = x(t) （4） uk wk uk δ(t) ∂t 式中，K 表示分割层数， 表示第 k 个子信号， 表 示 的中心频率，x 表示原信号， 为狄拉克函数， 代表梯度运算，*代表卷积运算，j 为虚数符号. 2.1 心电信号的预处理 根据香农采样定律，信号的频带宽度是采样 频率的 0.5 倍. 以 500 Hz 的采样频率为例，实际信 号的宽度为 0～250 Hz. 而心电信号的主要成分集 中于 0.05～49 Hz[16] . 因此，不对 50～250 Hz 的信 号进行分析，以提高实时性. 在此，采用了一个 4 阶 IIR 滤波器对这部分信号进行衰减，截止频率 为 50 Hz. 2.2 陡脉冲干扰的抑制 对式（ 2）进行傅里叶变换，可得其幅频曲线 A(f)，如式（5）. A(f) = se − t0 RC √ ( 1 RC ) 2 +(2π f) 2 （5） 式中，f 表示频率. 由于心电采样设备通常要求大 输入电阻[17] ，因此 1/RC 为常数且较小，在此为方 便分析将其忽略. 幅频特性如图 2. wk u1 w1 u1 u2 w2 u2 可以看到，陡脉冲干扰的后支能量主要分布 在 10 Hz 前. VMD 算法基于对混合信号的主成分 估计，结合其带通滤波器属性，包含陡脉冲干扰后 支的子信号对应中心频率 应小于 10 Hz. 由于此 部分干扰对心电波形干扰较大，是抑制的重点. 图 3 展示了当 K=6 时，各子信号时域图. 其中，子 信号 的中心频率 接近于 0，包含直流分量、基 线漂移、陡脉冲干扰的低频干扰分量. 将 抛弃， 可以有效地消除基线漂移和直流分量. 的中心 频率满足 <10 Hz，包含陡脉冲干扰分量、P 波、 T 波分量[16] . 完全抛弃 不利于恢复心电信号. u2 λsp 设心电采集系统共模抑制比为−100 dB，陡脉冲 的峰值大于 1000 V，理论上产生噪声的幅值约大于 10 mV. 实际测得幅值由于损耗而略小. 以肢体导 联为例，P 波幅值不超过 0.25 mV，T 波幅值不超过 0.5 mV[16] ，远低于同频段陡脉冲干扰幅值. 基于以上 分析，本文设计了一种自适应阈值，对陡脉冲干扰在 中的陡脉冲干扰分量进行消除. 阈值 如式（6）. λsp = mean(sort(|u2|)[1 : 10]) （6） u2 sort(|u2|)[1 : 10] || 上式表示对子信号 中绝对幅值最大的 10 个采样 点 取均值. sort( ) 表示对输入信号序 列的绝对值降序排列，mean(·) 表示均值运算. uˆ2( −→z )= u2( −→z ), |u2(r)| < λsp 0, |u2(r)| ⩾ λsp （7） −→z −→z −→z zon uˆ2 式（7）为阈值使用方法. 式中， 为采样点区间，包 含相邻两个零交叉点或两个近似零交叉点（所谓 近似零交叉点，即距离实际零交叉点最近的采样 点）之间的全部采样点. r 代表了 之间的峰值点. 其中，将峰值大于阈值的 的起点记作 . 处理前 后效果如图 4 所示，其结果用 表示. 5 10 0 0.5 1 1.5 15 20 25 30 Frequency/Hz Noise amplitude/dB 35 40 45 50 图 2 衰减模型单边幅频图 Fig.2 Diagram of a single-sided amplitude–frequency attenuation model −2 0 2 (a) 0 2 4 6 8 10 12 14 −1 0 1 (b) 0 2 4 6 8 10 12 14 −1 0 1 (c) 0 2 4 6 8 10 12 14 −0.2 0 0.2 (d) 0 2 4 6 8 10 12 14 −0.2 0 0.2 (e) 0 2 4 6 8 10 12 14 −0.2 0 0.2 (f) 0 2 4 6 8 10 12 14 −0.2 0 0.2 (g) 0 2 4 6 8 10 12 14 Amplitude/mV Time/s 图 3 VMD 6 阶分解子信号. （a）原信号；（b）u1；（c）u2；（d）u3；（e）u4； （f）u5；（g）u6 Fig.3 Use of VMD to divide the signal into six layers: (a) original signal; (b) u1 ; (c) u2 ; (d) u3 ; (e) u4 ; (f) u5 ; (g) u6 · 656 · 工程科学学报，第 42 卷，第 5 期
姚晰童等：陡脉冲干扰下的心电信号滤波及QRS提取 657 (a) (a 0.5 0.2 年作+ -0.5 6 0 Time/s Time/s (b) (b) 0.5 0.2 事水 02 0 2 6 0 4 Time/s Time/s 图4？中干扰分量消除效果.(a)阈值处理前：(b)阈值处理后 图5消除随机噪声分量之后的子信号.(a)4:(b)5 Fig.4 Interference component elimination effect in(a)before using Fig.5 Sub-signal after eliminating random noise components:(a): thresholds;(b)after using thresholds (b)us 对于陡脉冲干扰前支，由于其在时域上与后 xfilter(t)=xreal(t)-fipluse(t)-fw(t) 支存在紧密关联，在2之后的子信号中，以 wk<12,w+1≥12 wk≤49.w+1>49 [2om-4zm+4]为范围，搜索极大值（由陡脉冲特 =i2(t)+ (t)+ ik(1) k=3 性，陡脉冲干扰前支时域上表现为时长短的特点， (wk-1<12,多12} (11) 因此，4取500Hz采样频率下的5个采样点，即 式中，iw是对心电信号中随机噪声的估计，fipluse是 10ms).将极大值所属脉冲置零以完全消除陡脉 算法对陡脉冲冲激的估计.，表示消除随机噪声 冲冲激 后的k 2.3消除随机噪声 随机噪声由电子器件产生，在频域分布均匀， 3基于VMD算法的QRS提取 根据VMD算法的带通滤波器组特性2，中的随 在心电信号中，QRS波群的能量集中在8~ 机噪声分量的幅值远远小于原信号随机噪声幅 16Hz之间%如式(12)所示，则可以提取QRS波 值.受基于wavelet的阈值滤波算法图启发，本文 群主要成分 提出了一种基于VMD的阈值法去随机噪声算法 由于P、T波幅值小，采用阈值消除相应子信号中 oRs0=∑40 (12) 8≤w联I61 的随机噪声分量可能得不偿失.因此需，对应的 在重组信号xQRs的基础上，对QRS波群进行 满足式(8)： 检测，可避免P、T波等信号对检测过程的干扰 Won≤wk≤wof (8) 为了突出R波的位置，对重组信号进行预处理.预 式中，wom和wo分别对应P、T波频带的上限与心 处理算法如式(13)、(14)、（15)所示，依次对重组 电信号主要频带的上限，分别取12Hz和49Hz6 信号序列表示差分、平方、积分运算.xdf、xg、mvw 对满足条件的子信号施加阈值，如式(9) 分别是差分、平方、积分运算的输出.式中，T表 [sgn((iT)u(iT)引-)lk(iT引>k ing(iT)= 示采样周期，W表示移动滑窗宽度，在此取80ms lk(iT引≤ 效果较为理想，这与R波的实际宽度有关咧效果 (9) 如图6(c)所示 其中，i表示采样点序号，T为采样周期，sgn()代表 xdifr(t)=xoRs(t)-xoRs(t-iT) (13) 阶跃函数，代表，的阈值.如下式(10)所示 xsg(t)=xdifr(t) (14) A=rabi (sum(ug(iT))-In N) (10) i=W-1 式中，N代表，中采样点的数量，a和b均为常数 Xmvw(t)= xsg(t-iT)/W (15) τ为比例因子.施加阈值后效果如图5. 最终，抑制噪声后的心电信号用xe表示.如 处理之后的心电信号，采用阈值法对R波峰 式(11)所示： 值点进行提取.初始阈值o由最初10个1s内预
u2 [zon −∆t , zon +∆t ∆t 对于陡脉冲干扰前支，由于其在时域上与后 支 存 在 紧 密 关 联 ， 在 之 后 的 子 信 号 中 ， 以 ] 为范围，搜索极大值（由陡脉冲特 性，陡脉冲干扰前支时域上表现为时长短的特点， 因此， 取 500 Hz 采样频率下的 5 个采样点，即 10 ms）. 将极大值所属脉冲置零以完全消除陡脉 冲冲激. 2.3 消除随机噪声 uk uk wk 随机噪声由电子器件产生，在频域分布均匀， 根据 VMD 算法的带通滤波器组特性[12] ， 中的随 机噪声分量的幅值远远小于原信号随机噪声幅 值. 受基于 wavelet 的阈值滤波算法[8] 启发，本文 提出了一种基于 VMD 的阈值法去随机噪声算法. 由于 P、T 波幅值小，采用阈值消除相应子信号中 的随机噪声分量可能得不偿失. 因此需 对应的 满足式（8）： won ⩽ wk ⩽ woff （8） won woff uk 式中， 和 分别对应 P、T 波频带的上限与心 电信号主要频带的上限，分别取 12 Hz 和 49 Hz[16] . 对满足条件的子信号 施加阈值，如式（9）. uˆk(iT) = { [sgn(uk(iT))](|uk(iT)|−λk) |uk(iT)| > λk 0 |uk(iT)| ⩽ λk （9） λk uk 其中，i 表示采样点序号，T 为采样周期，sgn(·) 代表 阶跃函数， 代表 的阈值. 如下式（10）所示. λk = τab k 2 (sum(uk(iT))·lnN) 1 2 （10） 式中，N 代表uk 中采样点的数量，a 和 b 均为常数[18] . τ 为比例因子. 施加阈值后效果如图 5. 最终，抑制噪声后的心电信号用xfilter 表示. 如 式（11）所示： xfilter(t) = xreal(t)−nˆpluse(t)−nˆw(t) = uˆ2(t)+ {k|wk<12 ∑,wk+1⩾12} k=3 uk(t)+ {k|wk⩽49 ∑,wk+1>49} {k|wk−1<12,wk⩾12} uˆk(t) （11） nˆw nˆpluse uˆk uk 式中， 是对心电信号中随机噪声的估计， 是 算法对陡脉冲冲激的估计. 表示消除随机噪声 后的 . 3 基于 VMD 算法的 QRS 提取 在心电信号中，QRS 波群的能量集中在 8～ 16 Hz 之间[16] . 如式（12）所示，则可以提取 QRS 波 群主要成分. xQRS(t) = ∑ {k|8⩽wk⩽16} uk(t) （12） xQRS xdiff xsq xmvw 在重组信号 的基础上，对 QRS 波群进行 检测，可避免 P、T 波等信号对检测过程的干扰. 为了突出 R 波的位置，对重组信号进行预处理. 预 处理算法如式（13）、（14）、（15）所示，依次对重组 信号序列表示差分、平方、积分运算. 、 、 分别是差分、平方、积分运算的输出. 式中，T 表 示采样周期，W 表示移动滑窗宽度，在此取 80 ms 效果较为理想，这与 R 波的实际宽度有关[19] . 效果 如图 6（c）所示. xdiff(t) = xQRS(t)− xQRS(t−iT) （13） xsq(t) = xdiff 2 (t) （14） xmvw(t) = i=∑ W−1 i=0 xsq(t−iT)/W （15） λ0 处理之后的心电信号，采用阈值法对 R 波峰 值点进行提取. 初始阈值 由最初 10 个 1 s 内预 0 0 2 −0.5 0.5 4 6 (a) Time/s Amplitude/mV 0 −0.5 0.5 (b) Amplitude/mV Time/s 0 2 4 6 图 4 u2 中干扰分量消除效果. （a） 阈值处理前；（b） 阈值处理后 Fig.4 Interference component elimination effect in u2 : (a) before using thresholds; (b) after using thresholds 0 0 2 −0.2 0.2 4 6 (a) Time/s Amplitude/mV 0 −0.2 0.2 (b) Amplitude/mV Time/s 0 2 4 6 图 5 消除随机噪声分量之后的子信号. （a） u4；（b） u5 Fig.5 Sub-signal after eliminating random noise components: (a) u4 ; (b) u5 姚晰童等： 陡脉冲干扰下的心电信号滤波及 QRS 提取 · 657 ·