工程科学学报 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之后的子信号中,以 wk49 [2om-4zm+4]为范围,搜索极大值(由陡脉冲特 =i2(t)+ (t)+ ik(1) k=3 性,陡脉冲干扰前支时域上表现为时长短的特点, (wk-1k 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|wk49} {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 ·
658 工程科学学报,第42卷,第5期 式(16) AR Pnoise +y(Pnoise-PORS) (16) 式中,R为当前R波的检测阈值,Pnoise和PoRs分别 是最近n个非R峰值点幅值的均值和最近n个 Time/s R峰值点幅值的均值.y为比例常数,当n取8时, (b) y取0.475效果最佳算法流程图如图7.图中, m-1和rm分别表示第m-1个和第m个已经被定位 0.5 的R波所在位置,ITmean表示连续n个RR间期长 度;tRp为不应期,此期间无R波,时长为200ms2o; max[.]表示求括号中的最大值 Time/s 4 实验和数据分析 (c) 实验部分采用MIT-BH数据库数据与实际 采集的心电信号共同验证本文中提出的算法 图8(a)所示,是天津市机器人与智能重点实验室 Time/s 研制的不可逆电脉冲消融设备样机,可以释放峰 图6R波检测的算法验证.(a)原信号:(b)重组的QRS:(c)R波增强 峰值范围1000~3000V的陡脉冲.为采集受陡脉 Fig.6 Verification of the R-wave detection algorithm:(a)original 冲相关噪声干扰的心电信号,本文采用TI公司 signal;(b)restructured QRS wave;(c)enhanced R wave ADS1298芯片,自行设计制作了心电信号采集设 处理后序列极大值的均值决定.随后,每检测到 备,如图8(b)所示.本文的实验部分,陡脉冲的峰 一个R波即对阈值进行更新.阈值更新算法如 峰值为1000V. Start Detection of peak points by threshold rs-rar1>I8p? No Yes Fa Ya-l No Take the maximum point in <1.5rTmon? -+TTmcan2:I&TTmeanp Yes Max[x v() No Whether it is R wave Yes Update threshold Insert new R wave No Whether the data transfer ends es End 图7R波检测的算法流程图 Fig.7 Flowchart of the R-wave detection algorithm
处理后序列极大值的均值决定. 随后,每检测到 一个 R 波即对阈值进行更新. 阈值更新算法如 式(16). λR = Pnoise +γ(Pnoise − PQRS) (16) λR Pnoise PQRS rm−1 rm rrmean tRP 式中, 为当前 R 波的检测阈值, 和 分别 是最近 n 个非 R 峰值点幅值的均值和最近 n 个 R 峰值点幅值的均值. γ 为比例常数,当 n 取 8 时, γ 取 0.475 效果最佳[19] . 算法流程图如图 7. 图中, 和 分别表示第 m−1 个和第 m 个已经被定位 的 R 波所在位置, 表示连续 n 个 RR 间期长 度; 为不应期,此期间无 R 波,时长为 200 ms[20] ; max[.] 表示求括号中的最大值. 4 实验和数据分析 实验部分采用 MIT−BIH 数据库数据与实际 采集的心电信号共同验证本文中提出的算法. 图 8(a)所示,是天津市机器人与智能重点实验室 研制的不可逆电脉冲消融设备样机,可以释放峰 峰值范围 1000~3000 V 的陡脉冲. 为采集受陡脉 冲相关噪声干扰的心电信号,本文采用 TI 公司 ADS1298 芯片,自行设计制作了心电信号采集设 备,如图 8(b)所示. 本文的实验部分,陡脉冲的峰 峰值为 1000 V. 0 1 2 0 2 4 6 8 (a) Time/s Amplitude/mV 0 1 2 0 2 4 6 8 (c) Time/s Amplitude/mV 0 0 2 −0.5 0.5 4 6 8 (b) Time/s Amplitude/mV 图 6 R 波检测的算法验证. (a) 原信号;(b) 重组的 QRS;(c) R 波增强 Fig.6 Verification of the R-wave detection algorithm: (a) original signal; (b) restructured QRS wave; (c) enhanced R wave Start End Detection of peak points by threshold rm−rm−1>tRP? rm−rm−1 <1.5rrmean? Whether it is R wave Whether the data transfer ends Yes Yes Yes Yes No No No No Take the maximum point in [rk−1+rrmean/2, rk−rrmean/2] Max[xmvw(rm), xmvw(rm−1)] Update threshold Insert new R wave 图 7 R 波检测的算法流程图 Fig.7 Flowchart of the R-wave detection algorithm · 658 · 工程科学学报,第 42 卷,第 5 期
姚晰童等:陡脉冲干扰下的心电信号滤波及QRS提取 659 a) b 满足{K∈[5,9]}.为进一步确定参数,设计以下实 验:选取MIT-BIH数据库中第100号数据,取 7200个采样点,在随机噪声级别均为5dB的条件 下,添加一定数量的模拟高压陡脉冲干扰.在此, 以信噪比(Signal to noise ratio,SNR)作为评价指 标,如式(17).式中,Xclean表示未添加噪声的心电 信号,Xnoise表示添加了噪声后的心电信号,Nl表示 全部参加测试的采样点数 图8实验设备展示.(a)不可逆电消融样机:(b)心电信号采集板 SNR=10log Fig.8 Experimental equipment:(a)prototype of irreversible electrical iean(/∑(ise(-Xcleam()} pulse ablation surgery.(b)ECG signal sampling circuit (17) 4.1ECG滤波算法实验 实验结果如图10所示.实验结果表明,当 4.1.1预处理对实时性的提升 K=6且a=2500时,滤波效果较为理想 实验平台基本参数:主频2.9GHz,i5处理器,运 20 *一K=5 行内存8G,Windows10操作系统.本次实验选用 6—K=6 500Hz采样频率下16.384s的包含陡脉冲干扰的 --K=7 -K=8 ECG信号为实验对象 十一K=9 10 图9对比了心电信号是否预处理对w2的影响 NS 根据VMD相关理论,分解层数K应当尽可能的少 以节约计算时间;同时表示基线信息因而w约等 于0,w2的值意味着分解的有效性.当w2<10Hz时 5001000150020002500300035004000 有效分离高压陡脉冲干扰.若不进行预处理,分解 层数K≥14时可将陡脉冲干扰与QRS波群分离.当 图10对分解层数和惩罚因子进行验证 K=14时,耗时30.1939s.如对ECG进行预处理, Fig.10 Values of K and quadratic penalty were verified through 当K≥5时即满足要求.当K=5时,仅需2.511s.显 experiments 然,对带噪心电信号预处理可以有效提高实时性 4.1.3陡脉冲干扰抑制 且并较少内存耗费 图11为对心电信号中陡脉冲干扰抑制效果 其中,图11(a)为原信号图,图11(b)和图11(c)分 60 Original signal 别是选择保留和不保留陡脉冲前支的心电信号 50 eAfter pre-processing 图.可以看到,运用本算法,对陡脉冲冲激有着明 0 显的消除效果,且最大程度上保持了心电信号应 30 有波形.图12对比了本文提出算法与小波包阈值 去燥图]、带阻滤波算法(采用FR带阻滤波器滤除 立20 0~10Hz信号)在消除陡脉冲干扰时的效果.其中 10 图l2(a)为本文提出算法对陡脉冲干扰的消除效 99999 果,图12(b)为小波包阈值去燥算法对陡脉冲干扰 123456789101112131415 的消除效果,虽然这种算法相对较好的保留了心 图9信号预处理效果对比 电信号的特征,但是对陡脉冲干扰消除并不彻底; Fig.9 Effect of signal preprocessing interference 图12(c)为采用带阻滤波器算法对陡脉冲干扰的 消除效果,此算法则完全不能保留心电信号特征, 4.1.2惩罚因子a和分解层数K的确定 且抑噪效果差 惩罚因子a和分解层数K是VMD算法两个 为进一步验证算法效果,选择MIT-BH数据 重要的参数,决定了滤波器通频带宽度1.由图9, 库中10段含噪声较少的心电数据,每段数据长度 对心电信号进行预处理的条件下,分解层数K应 为16.384s:由式(5)的数学模型,分别向心电信号
4.1 ECG 滤波算法实验 4.1.1 预处理对实时性的提升 实验平台基本参数:主频 2.9 GHz,i5 处理器,运 行内存 8 G,Windows10 操作系统. 本次实验选用 500 Hz 采样频率下 16.384 s 的包含陡脉冲干扰的 ECG 信号为实验对象. w2 u1 w1 w2 w2 图 9 对比了心电信号是否预处理对 的影响. 根据 VMD 相关理论,分解层数 K 应当尽可能的少 以节约计算时间;同时 表示基线信息因而 约等 于 0, 的值意味着分解的有效性. 当 <10 Hz 时 有效分离高压陡脉冲干扰. 若不进行预处理,分解 层数 K≥14 时可将陡脉冲干扰与 QRS 波群分离. 当 K=14 时,耗时 30.1939 s. 如对 ECG 进行预处理, 当 K≥5 时即满足要求. 当 K=5 时,仅需 2.511 s. 显 然,对带噪心电信号预处理可以有效提高实时性 且并较少内存耗费. 4.1.2 惩罚因子 α 和分解层数 K 的确定 惩罚因子 α 和分解层数 K 是 VMD 算法两个 重要的参数,决定了滤波器通频带宽度[10] . 由图 9, 对心电信号进行预处理的条件下,分解层数 K 应 xclean xnoise Nal 满足{K|K∈[5,9]}. 为进一步确定参数,设计以下实 验 : 选 取 MIT−BIH 数 据 库 中 第 100 号 数 据 , 取 7200 个采样点,在随机噪声级别均为 5 dB 的条件 下,添加一定数量的模拟高压陡脉冲干扰. 在此, 以信噪比( Signal to noise ratio, SNR)作为评价指 标,如式(17). 式中, 表示未添加噪声的心电 信号, 表示添加了噪声后的心电信号, 表示 全部参加测试的采样点数. SNR = 10log ∑ Nal i=1 x 2 clean(i)/ ∑ Nal i=1 (xnoise(i)− xclean(i))2 (17) 实验结果如图 10 所示. 实验结果表明 ,当 K=6 且 α=2500 时,滤波效果较为理想. 4.1.3 陡脉冲干扰抑制 图 11 为对心电信号中陡脉冲干扰抑制效果. 其中,图 11(a)为原信号图,图 11(b)和图 11(c)分 别是选择保留和不保留陡脉冲前支的心电信号 图. 可以看到,运用本算法,对陡脉冲冲激有着明 显的消除效果,且最大程度上保持了心电信号应 有波形. 图 12 对比了本文提出算法与小波包阈值 去燥[8]、带阻滤波算法(采用 FIR 带阻滤波器滤除 0~10 Hz 信号)在消除陡脉冲干扰时的效果. 其中 图 12(a)为本文提出算法对陡脉冲干扰的消除效 果,图 12(b)为小波包阈值去燥算法对陡脉冲干扰 的消除效果,虽然这种算法相对较好的保留了心 电信号的特征,但是对陡脉冲干扰消除并不彻底; 图 12(c)为采用带阻滤波器算法对陡脉冲干扰的 消除效果,此算法则完全不能保留心电信号特征, 且抑噪效果差. 为进一步验证算法效果,选择 MIT−BIH 数据 库中 10 段含噪声较少的心电数据,每段数据长度 为 16.384 s;由式(5)的数学模型,分别向心电信号 (a) (b) 图 8 实验设备展示. (a) 不可逆电消融样机;(b) 心电信号采集板 Fig.8 Experimental equipment: (a) prototype of irreversible electrical pulse ablation surgery; (b) ECG signal sampling circuit 1 2 3 4 5 6 7 8 9 0 10 20 30 40 50 60 K 10 11 12 13 14 15 Original signal After pre-processing Frequency/Hz 图 9 信号预处理效果对比 Fig.9 Effect of signal preprocessing interference 0 5 500 1000 1500 2000 2500 SNR/dB 3000 4000 3500 10 15 20 α K=5 K=6 K=7 K=8 K=9 图 10 对分解层数和惩罚因子进行验证 Fig.10 Values of K and quadratic penalty were verified through experiments 姚晰童等: 陡脉冲干扰下的心电信号滤波及 QRS 提取 · 659 ·
660 工程科学学报,第42卷,第5期 (a) (a) 5 Auapnujdurv 420 -2 -4 0 4 0 2 Time/s Time/s (b) (b) 4 2 一d -4 2 2 Time/s Time/s (c) (c) 420 -4 -4 2 3 0 2 3 4 Time/s Time/s 图11施放陡脉冲期间相关噪声抑制.(a)陡脉冲干扰消除前:(b) 图12陡脉冲干扰抑制效果对比.(a)本文提出算法效果;(b)小波 消除陡脉冲干扰后支:(c)完全消除陡脉冲干扰 算法效果:(c)带阻滤波器算法效果 Fig.11 Result of eliminating steep pulse interference:(a)ECG Fig.12 Comparison of the results of steep pulse noise filtering:(a) disturbed by steep pulses,(b)forehead is retained,(c)fore branches are algorithm in this article;(b)wavelet algorithm;(c)bandpass filter eliminated 以认为文献[I8]中a,b的取值是无误的 中加入不同数量的模拟高压陡脉冲干扰.分别采 再分别将信噪比为1~10dB的随机噪声和 用本文提出的算法、小波包阈值法、带阻滤波器 相同数量的模拟陡脉冲加入到等长的心电数 对噪声进行处理.实验结果如表1所示 据中,以此数据对算法进行对比验证.表2为将 表中,第一列表示加入不同数量模拟陡脉冲 本文提出的算法与wavelet去噪法、EMD去燥 干扰后的原始信号信噪比,其后三列分别表示采 法的效果对比.信噪比越大说明信号包含随机 用不同滤波算法处理后的滤波效果,以信噪比表 噪声的强度越小,也就意味着对噪声的消除效果 示.可以从以上实验结果中得出结论,本文设计算 越好 法可有效去除心电信号中的高压陡脉冲干扰,且 可以得出结论,在消除消融设备工作带来的 相比于其他算法具有明显的优势 随机噪声方面,本算法可以明显的提高信号的信 4.1.4随机噪声消除效果 噪比;相较于另外两种算法,提升效果更佳明显 阈值公式(10)由文献[18提出,式中参数a,b分 4.2QRS波群检测实验和数据分析 别取0.083和1.39.再此对参数正确性进行验证 对于R波检测准确度,采用两个指标进行评价: 采用一段来自MIT-BIH数据库第103条心电数 (I)误检率(Error detection,.ED):表达式如式 据,加入一定数量的模拟陡脉冲干扰以及5dB的 (18)所示; 随机噪声模拟真实环境.实验结果如图13所示: (2)敏感度(Sensitivity,.SE):表达式如式(19) 当a和b改变时,信噪比不同程度降低.由此,可 所示. 表1陡脉冲干扰消除效果对比 Table 1 Results obtained after a variety of filtering processes SNR of the original signal/dB Bandpass filter(SNR)/dB Wavelet packet filter(SNR)/dB Algorithm of this paper(SNR)/dB 1.481 -5369 0.199 8.022 -6.732 -11.456 -8.710 3.070 -12.110 -18.928 -13.913 -0.410
中加入不同数量的模拟高压陡脉冲干扰. 分别采 用本文提出的算法、小波包阈值法、带阻滤波器 对噪声进行处理. 实验结果如表 1 所示. 表中,第一列表示加入不同数量模拟陡脉冲 干扰后的原始信号信噪比,其后三列分别表示采 用不同滤波算法处理后的滤波效果,以信噪比表 示. 可以从以上实验结果中得出结论,本文设计算 法可有效去除心电信号中的高压陡脉冲干扰,且 相比于其他算法具有明显的优势. 4.1.4 随机噪声消除效果 阈值公式(10)由文献 [18] 提出,式中参数 a,b 分 别取 0.083 和 1.39. 再此对参数正确性进行验证. 采用一段来自 MIT−BIH 数据库第 103 条心电数 据,加入一定数量的模拟陡脉冲干扰以及 5 dB 的 随机噪声模拟真实环境. 实验结果如图 13 所示: 当 a 和 b 改变时,信噪比不同程度降低. 由此,可 以认为文献 [18] 中 a,b 的取值是无误的. 再分别将信噪比为 1~10 dB 的随机噪声和 相同数量的模拟陡脉冲加入到等长的心电数 据中,以此数据对算法进行对比验证. 表 2 为将 本文提出的算法与 wavelet 去噪法[8]、 EMD 去燥 法[9] 的效果对比. 信噪比越大说明信号包含随机 噪声的强度越小,也就意味着对噪声的消除效果 越好. 可以得出结论,在消除消融设备工作带来的 随机噪声方面,本算法可以明显的提高信号的信 噪比;相较于另外两种算法,提升效果更佳明显. 4.2 QRS 波群检测实验和数据分析 对于 R 波检测准确度,采用两个指标进行评价: (1)误检率(Error detection, ED):表达式如式 (18)所示; (2)敏感度(Sensitivity, SE):表达式如式(19) 所示. 表 1 陡脉冲干扰消除效果对比 Table 1 Results obtained after a variety of filtering processes SNR of the original signal/dB Bandpass filter(SNR)/dB Wavelet packet filter(SNR)/dB Algorithm of this paper(SNR)/dB 1.481 −5.369 0.199 8.022 −6.732 −11.456 −8.710 3.070 −12.110 −18.928 −13.913 −0.410 0 0 1 −5 5 2 3 4 (a) Time/s Amplitude/mV 0 0 1 −4 4 2 −2 2 3 4 (b) Time/s Amplitude/mV 0 0 2 1 −4 −2 4 2 3 4 (c) Time/s Amplitude/mV 图 11 施放陡脉冲期间相关噪声抑制. (a) 陡脉冲干扰消除前;(b) 消除陡脉冲干扰后支;(c) 完全消除陡脉冲干扰 Fig.11 Result of eliminating steep pulse interference: (a) ECG disturbed by steep pulses; (b) forehead is retained; (c) fore branches are eliminated 0 2 0 −2 1 −4 4 2 3 4 (a) Time/s Amplitude/mV 0 0 1 −4 4 2 −2 2 3 4 (b) Time/s Amplitude/mV 0 0 2 1 −4 −2 4 2 3 4 (c) Time/s Amplitude/mV 图 12 陡脉冲干扰抑制效果对比. (a) 本文提出算法效果;(b) 小波 算法效果;(c) 带阻滤波器算法效果 Fig.12 Comparison of the results of steep pulse noise filtering: (a) algorithm in this article; (b) wavelet algorithm; (c) bandpass filter · 660 · 工程科学学报,第 42 卷,第 5 期
姚晰童等:陡脉冲干扰下的心电信号滤波及QRS提取 661 (a) TP 15 SE=TP+FN 100% (19) 10 式中,TP表示正确检测,FP表示漏检的情况,FN 表示误检的情况,Niotal表示R波总数量. 0 0.01 006 02 03 为验证本算法的实际表现,采用8段来自不同 患者的心电信号数据.每组数据集的时间长度均 (b) 为30min.考虑相关陡脉冲疗法的治疗位置可能 15 位于躯干,因此,采用肢体导联作为心电信号来 源,避免干扰正常治疗过程.经过人工标记,共标 记16218个R峰值点.此期间,共记录释放陡脉冲 0.01 0.5 11.391.5 522次.将本文提出算法与差分阈值法0,带通滤 b 波器法的对比结果,实验数据如表3所示.经过 图13a,b最优值验证.(a)b恒定a改变:(b)a恒定b改变 对比发现,相比于本文提出的算法,其他两种算法 Fig.13 Verification of the optimal values of a and b:(a)when b is 的误检率远高于本算法,这与未对陡脉冲干扰进 constant and a is changed;(b)when a is constant and b is changed 行抑制由直接关系.可以通过数据得出结论,本文 FP+FN ED= .100% (18) 提出的算法,在应对陡脉冲干扰下的心电信号时, Ntotal 有着更好的准确度和敏感度 表2去噪效果对比 Table 2 Comparison of the results of the filtering processes SNR between ECG and random noise/dB Algorithm Random noise SNR=1 dB Random noise SNR=3 dB Random noise SNR=8 dB Random noise SNR=10 dB Wavelet 7.121 8.116 11.659 15.227 EMD网 6.129 8.900 12.387 14.981 This paper 8.039 10.122 14.649 20.836 表3与常用算法准确度对比 效果较为理想.其中,对陡脉冲相关噪声干扰下的 Table 3 Comparison of the accuracy of the algorithms ECG,提取QRS波群的准确度达到了99%以上. Algorithm Ntotal TP FP FN ED/%SE/% 参考文献 Difference threshold algorithm 14157765129612.71091.61 Bandpass filter algorithm 16218151383777036.65995.56 [1]Prausnitz M R,Bose V G,Langer R,et al.Electroporation of mammalian skin:a mechanism to enhance transdermal drug Algorithm in this article 157413601172.94299.23 delivery.Proc Natl Acad Sci USA,1993,90(22):10504 [2]Lu D H.Zhang J X,Dai Y.Analysis and compensation of stray 5结论 parameters in high-voltage pulse generator circuits.Chin/Eng, 本文以陡脉冲的应用为背景,对这种特殊环 2016.38(5):734 (吕东浒,张建勋,代煜.陡脉冲发生器电路中杂散参数的分析 境下采集的心电信号进行了研究.创新点在于:针 和补偿.工程科学学报,2016,38(5):734) 对陡脉冲干扰特性,本文提出了一种基于VMD的 [3] Gehl J.Electroporation:theory and methods,perspectives for drug 心电信号陡脉冲干扰抑制算法;基于VMD特性, delivery,gene therapy and research.Acta Physiol Scand,2003. 提出了一种阈值算法,有效的消除了电子设备带 177(4):437 来的随机噪声;并根据分割后子信号的频域特性, [4] Jin S L,Wang Y,Yin Y G.Effect of high intensity pulsed electric 提出了一种基于VMD的QRS波群提取算法.最 field pretreatment on elasticity of starch-myofibril protein mixed 终,以在天津市机器人与智能重点实验室研制的 gelatin.JJilin Univ Eng Technol Ed,2014,44(2):573 (金声琅,王莹,殷涌光.高压脉冲电场对淀粉和肌原纤维蛋白 不可逆电消融设备干扰下采集的心电信号,与来自 混合凝胶强度的影响.吉林大学学报:工学版,2014,44(2): MT-BIH数据库中的心电数据共同对算法性能进 573) 行了验证.结果证明,算法可以实现预期目标,且 [5]Feng X Q,Wang Y H,Xu F X.Advances on application of high
ED = FP+FN Ntotal · 100% (18) SE = TP TP+FN · 100% (19) Ntotal 式中,TP 表示正确检测,FP 表示漏检的情况,FN 表示误检的情况, 表示 R 波总数量. 为验证本算法的实际表现,采用 8 段来自不同 患者的心电信号数据. 每组数据集的时间长度均 为 30 min. 考虑相关陡脉冲疗法的治疗位置可能 位于躯干,因此,采用肢体导联作为心电信号来 源,避免干扰正常治疗过程. 经过人工标记,共标 记 16218 个 R 峰值点. 此期间,共记录释放陡脉冲 522 次. 将本文提出算法与差分阈值法[20] ,带通滤 波器法[21] 的对比结果,实验数据如表 3 所示. 经过 对比发现,相比于本文提出的算法,其他两种算法 的误检率远高于本算法,这与未对陡脉冲干扰进 行抑制由直接关系. 可以通过数据得出结论,本文 提出的算法,在应对陡脉冲干扰下的心电信号时, 有着更好的准确度和敏感度. 5 结论 本文以陡脉冲的应用为背景,对这种特殊环 境下采集的心电信号进行了研究. 创新点在于:针 对陡脉冲干扰特性,本文提出了一种基于 VMD 的 心电信号陡脉冲干扰抑制算法;基于 VMD 特性, 提出了一种阈值算法,有效的消除了电子设备带 来的随机噪声;并根据分割后子信号的频域特性, 提出了一种基于 VMD 的 QRS 波群提取算法. 最 终,以在天津市机器人与智能重点实验室研制的 不可逆电消融设备干扰下采集的心电信号,与来自 MIT−BIH 数据库中的心电数据共同对算法性能进 行了验证. 结果证明,算法可以实现预期目标,且 效果较为理想. 其中,对陡脉冲相关噪声干扰下的 ECG,提取 QRS 波群的准确度达到了 99% 以上. 参 考 文 献 Prausnitz M R, Bose V G, Langer R, et al. Electroporation of mammalian skin: a mechanism to enhance transdermal drug delivery. Proc Natl Acad Sci USA, 1993, 90(22): 10504 [1] Lü D H, Zhang J X, Dai Y. Analysis and compensation of stray parameters in high-voltage pulse generator circuits. Chin J Eng, 2016, 38(5): 734 (吕东澔, 张建勋, 代煜. 陡脉冲发生器电路中杂散参数的分析 和补偿. 工程科学学报, 2016, 38(5):734) [2] Gehl J. Electroporation: theory and methods, perspectives for drug delivery, gene therapy and research. Acta Physiol Scand, 2003, 177(4): 437 [3] Jin S L, Wang Y, Yin Y G. Effect of high intensity pulsed electric field pretreatment on elasticity of starch-myofibril protein mixed gelatin. J Jilin Univ Eng Technol Ed, 2014, 44(2): 573 (金声琅, 王莹, 殷涌光. 高压脉冲电场对淀粉和肌原纤维蛋白 混合凝胶强度的影响. 吉林大学学报: 工学版, 2014, 44(2): 573) [4] [5] Feng X Q, Wang Y H, Xu F X. Advances on application of high 表 2 去噪效果对比 Table 2 Comparison of the results of the filtering processes Algorithm SNR between ECG and random noise /dB Random noise SNR=1 dB Random noise SNR=3 dB Random noise SNR=8 dB Random noise SNR=10 dB Wavelet[8] 7.121 8.116 11.659 15.227 EMD[9] 6.129 8.900 12.387 14.981 This paper 8.039 10.122 14.649 20.836 表 3 与常用算法准确度对比 Table 3 Comparison of the accuracy of the algorithms Algorithm Ntotal TP FP FN ED/% SE/% Difference threshold algorithm[20] 16218 14157 765 1296 12.710 91.61 Bandpass filter algorithm[21] 15138 377 703 6.659 95.56 Algorithm in this article 15741 360 117 2.942 99.23 0.01 0.06 0.083 10 5 0 15 0.1 0.2 0.3 (a) a SNR/dB 0.01 0.5 1.39 10 5 0 15 1 1.5 2 (b) b SNR/dB 图 13 a,b 最优值验证. (a) b 恒定 a 改变;(b) a 恒定 b 改变 Fig.13 Verification of the optimal values of a and b: (a) when b is constant and a is changed; (b) when a is constant and b is changed 姚晰童等: 陡脉冲干扰下的心电信号滤波及 QRS 提取 · 661 ·
662 工程科学学报,第42卷,第5期 voltage pulsed electric field in food quality and safety.J Food Sci 动轴承故障诊断.工程科学学报,2016,38(9):1327) Biotechnol,2013,32(4):337 [14]Xue Y J,Cao J X,Wang D X,et al.Application of the variational- (冯叙桥,王月华,徐方旭.高压脉冲电场技术在食品质量与安 mode decomposition for seismic time-frequency analysis.IEEEJ 全中的应用进展.食品与生物技术学报,2013,32(4):337) Selected Topics Appl Earth Observ Remote Sens,2016,9(8):3821 [6] Maor E,Ivorra A,Leor J,et al.The effect of irreversible [15]Liu YY,Yang G L,Li M,et al.Variational mode decomposition electroporation on blood vessels.Technol Cancer Res Treat,2007, denoising combined the detrended fluctuation analysis.Signa/ 6(4):307 Proces3,2016,125:349 [7]Zhang R Y,Fang X,Liu Y Q,et al.Design of a real-time ECG [16]Cao X W,Deng Q K.Frequency analysis on the ECG waveform. filter for resource constraint computer 2012 International Chin J Med Phys,2001,18(1):46 Conference on Information and Automation (ICIA).Shenyang, (曹细武,邓亲恺.心电图各波的颜率分析,中国医学物理学杂 2012:846 志,2001,18(1):46) [8]W,Chen Y B,et al.ECG denoising method based on [17]Lu Y L,Wang R,Hu M L,et al.Development of wearable ECG adaptive wavelet threshold selection.Comput Eng Appl,2018. and SCG monitoring system.Beijing Biomed Eng,2018,37(4): 54(15):29 381 (王磊,孙玮,陈奕博,等.基于自适应小波值的心电信号降噪 (鲁亚磊,王瑞,胡敏露,等.穿戴式心电及心振信号监测系统的 方法.计算机工程与应用,2018.54(15):29) 研制.北京生物医学工程,2018,37(4):381) [9]Tan X,Chen XX,Hu X Y,et al.EMD-based electrocardiogram [18]Kopsinis Y,McLaughlin S.Development of EMD-based denoising delineation for a wearable low-power ECG monitoring device.Can methods inspired by wavelet thresholding.IEEE Trans signal J Electr Comput Eng,2014,37(4):212 [10]Dragomiretskiy K,Zosso D.Variational mode decomposition. Process,2009,57(4):1351 IEEE Trans Signal Process,2014,62(3):531 [19]Hamilton P S,Tompkins W J.Quantitative investigation of QRS [11]Chen S.Guo M Z,Zhang J X.Design and realization of high- detection rules using the MIT/BIH arrhythmia database.IEEE voltage steep pulse generator.Automnstr,2019,34(2):5 Trans Biomed Eng,1986,33(12):1157 (陈穗,郭蒙召,张建勋.高压陡脉冲发生器的设计与实现.自动 [20]Liu Z D,Chen J,Tang M F.Pulse transit time detection based on 化与仪表,2019,34(2):5) waveform time domain feature and dynamic difference threshold.J [12]Wang Y X,Markert R.Filter bank property of variational mode Biomed Eng,2017,34(3):329 decomposition and its applications.Signal Process,2016,120:509 (刘增丁,陈摸,汤敏芳.基于波形时域特征和动态差分倒值的 [13]Zhang D,Feng Z P.Fault diagnosis of rolling bearings based on 脉搏波传导时间检测算法.生物医学工程学杂志,2017,34(3): variational mode decomposition and calculus enhanced energy 329) operator.Chin J Eng,2016,38(9):1327 [21]Pan J P,Tompkins W J.A real-time QRS detection algorithm. (张东,冯志鹏.基于变分模式分解和微积分增强能量算子的滚 IEEE Trans Biomed Eng,1985,32(3):230
voltage pulsed electric field in food quality and safety. J Food Sci Biotechnol, 2013, 32(4): 337 (冯叙桥, 王月华, 徐方旭. 高压脉冲电场技术在食品质量与安 全中的应用进展. 食品与生物技术学报, 2013, 32(4):337) Maor E, Ivorra A, Leor J, et al. The effect of irreversible electroporation on blood vessels. Technol Cancer Res Treat, 2007, 6(4): 307 [6] Zhang R Y, Fang X, Liu Y Q, et al. Design of a real-time ECG filter for resource constraint computer // 2012 International Conference on Information and Automation (ICIA). Shenyang, 2012: 846 [7] Wang L, Sun W, Chen Y B, et al. ECG denoising method based on adaptive wavelet threshold selection. Comput Eng Appl, 2018, 54(15): 29 (王磊, 孙玮, 陈奕博, 等. 基于自适应小波阈值的心电信号降噪 方法. 计算机工程与应用, 2018, 54(15):29) [8] Tan X, Chen X X, Hu X Y, et al. EMD-based electrocardiogram delineation for a wearable low-power ECG monitoring device. Can J Electr Comput Eng, 2014, 37(4): 212 [9] Dragomiretskiy K, Zosso D. Variational mode decomposition. IEEE Trans Signal Process, 2014, 62(3): 531 [10] Chen S, Guo M Z, Zhang J X. Design and realization of highvoltage steep pulse generator. Autom Instrum, 2019, 34(2): 5 (陈穗, 郭蒙召, 张建勋. 高压陡脉冲发生器的设计与实现. 自动 化与仪表, 2019, 34(2):5) [11] Wang Y X, Markert R. Filter bank property of variational mode decomposition and its applications. Signal Process, 2016, 120: 509 [12] Zhang D, Feng Z P. Fault diagnosis of rolling bearings based on variational mode decomposition and calculus enhanced energy operator. Chin J Eng, 2016, 38(9): 1327 (张东, 冯志鹏. 基于变分模式分解和微积分增强能量算子的滚 [13] 动轴承故障诊断. 工程科学学报, 2016, 38(9):1327) Xue Y J, Cao J X, Wang D X, et al. Application of the variationalmode decomposition for seismic time-frequency analysis. IEEE J Selected Topics Appl Earth Observ Remote Sens, 2016, 9(8): 3821 [14] Liu Y Y, Yang G L, Li M, et al. Variational mode decomposition denoising combined the detrended fluctuation analysis. Signal Process, 2016, 125: 349 [15] Cao X W, Deng Q K. Frequency analysis on the ECG waveform. Chin J Med Phys, 2001, 18(1): 46 (曹细武, 邓亲恺. 心电图各波的频率分析. 中国医学物理学杂 志, 2001, 18(1):46) [16] Lu Y L, Wang R, Hu M L, et al. Development of wearable ECG and SCG monitoring system. Beijing Biomed Eng, 2018, 37(4): 381 (鲁亚磊, 王瑞, 胡敏露, 等. 穿戴式心电及心振信号监测系统的 研制. 北京生物医学工程, 2018, 37(4):381) [17] Kopsinis Y, McLaughlin S. Development of EMD-based denoising methods inspired by wavelet thresholding. IEEE Trans signal Process, 2009, 57(4): 1351 [18] Hamilton P S, Tompkins W J. Quantitative investigation of QRS detection rules using the MIT/BIH arrhythmia database. IEEE Trans Biomed Eng, 1986, 33(12): 1157 [19] Liu Z D, Chen J, Tang M F. Pulse transit time detection based on waveform time domain feature and dynamic difference threshold. J Biomed Eng, 2017, 34(3): 329 (刘增丁, 陈骥, 汤敏芳. 基于波形时域特征和动态差分阈值的 脉搏波传导时间检测算法. 生物医学工程学杂志, 2017, 34(3): 329) [20] Pan J P, Tompkins W J. A real-time QRS detection algorithm. IEEE Trans Biomed Eng, 1985, 32(3): 230 [21] · 662 · 工程科学学报,第 42 卷,第 5 期