融合粒子群与双边长短窗能量差的心电P波检测

2022-05-18徐礼胜苏宇剑谭浚宜方喜冬

东北大学学报（自然科学版） 2022年4期

徐礼胜, 苏宇剑, 谭浚宜, 方喜冬

(1. 东北大学医学与生物信息工程学院, 辽宁沈阳 110169; 2. 沈阳东软智能医疗科技研究院有限公司, 辽宁沈阳 110167)

心电图分析可用于诊断心血管疾病,大多数诊断所需数据都可以从心电图波形的间隔、幅度和形状中获得[1],所以准确且可靠提取心电图特征十分重要.心电信号是在皮肤表面获得的由心脏极化和去极化产生的电信号,这种电信号具有周期性,每个心电周期按照前后顺序分为P波、QRS波群、T波和U波.其中:P波代表心房去极化,QRS波群代表心室去极化,T波代表心室去极化后再极化[1].虽然心电图中P波幅值小、时长短,但P波蕴含了丰富的病理信息.P波是判断室性早搏和正常心拍的主要判据之一,是心房功能能否正常的主要判断依据.P波起点到QRS波群的起点间的时长是判断窦性心动过速、预激综合征等的重要指标[2],因此,心电P波的准确提取对生理健康判断具有重要意义.

P波幅值较低,在不同导联中的形态也有所不同,心电P波按照波峰的形态可分为正向、倒置与双向三种[3].目前较好的P波提取方法有自适应滤波法、高斯拟合法、神经网络法等,但这些检测方法都以正向P波为主[4-6].Mohamed等[4]对心电波形首先采用带通滤波,然后使用两个移动的平均滤波器检测P波和T波,该方法的缺点是滤波器需要精准参数,而且参数无法自适应[4].Rao等[5]使用基于动态规划的参数混合高斯拟合法,先预设一种与P波类似的高斯模型,并与心电信号对比,与其模型最相近的部分为心电P波.这种方法只能检测出正向P波,对于双向P波则无法检测.Lenis等[6]运用了无相位小波变换对心电P波进行提取分析,并取得良好结果,但未考虑到P波有正向、负向和双向三种形态.目前检测P波的方法都是从波形时域形态或频域方面检测,但P波波形有三种不同形态,很难使用一种方法同时检测所有形态的P波.因为心电低频噪声与P波在频域重合,所以频域法需要心电信号有非常好的质量.本文提出了一种融合粒子群和双边长短窗能量差检测多形态心电P波的算法,此算法通过能量变换检测出所有形态的心电P波,并具有一定的抗噪能力.

1 方法与步骤

本文算法主要通过双边长短窗能量差法提取心电P波起点与终点的特征,并使用粒子群来优化其参数,最终确定P波的起点与终点.算法具体步骤如图1所示.

图1 算法具体步骤

首先对信号进行预处理,包括去除基线漂移和高频噪声;再根据R波波峰将心电信号分离为单个心电周期,并找出预选区域.根据Q(t)找到P波的起点与终点,在计算过程中使用粒子群算法得到最优解.在图1中,ag,Lg分别代表粒子群算法的初始解,其中ag为短窗长,Lg为长窗长比短窗长的倍数.a和L分别为每次通过粒子群算法迭代所得的短窗长、长窗长比短窗长的倍数.Xs和Xe分别为通过双边长短窗平均能量差法计算所得心电P波的起点与终点.图1中的目标函数是使用当前迭代所得参数提取P波的起点与终点,将其结果与金标准对比,计算标准差,其结果越小则参数优越性越高.

1.1 双边长短窗平均能量差法

在双边长短窗平均能量差法的基础上,对长短时窗能量差法加以改进,基础的长短时窗能量差法广泛用于微地震波起点的检测[7-8].微地震波是指在山体施工、土木建设、能源开发或汽油采集等过程中引起地下应力场变化,导致岩体破裂,发生微小震动而产生的地震波,其幅值低、频率低、含有大量高频噪声,与人体心电类似[9].本研究将基础的长短时窗平均能量差法加以改进,用于检测心电P波的起点与终点.

双边长短时窗平均能量差法的原理是通过能量变化来判断信号的边界.基础的长短时窗平均能量差法通过计算,将起点前没有信号的部分判断为低能量区,将起点后含有信号的部分判断为高能量区.当长短窗都位于起点前,t点位于起点前时由于长窗和短窗的区域都为噪声,能量值偏低,长窗与短窗的平均能量差值接近于0;继续向后滑动,一部分长窗包含信号,进入了高能量区,差值开始减少;继续滑动,全部长窗进入高能量区,短窗在低能量区达到极小值;再继续滑动,短窗进入高能量区,长窗不变,差值变大.具体计算方式如下:

1) 先选定一组长短窗,短窗包含l个采样点,长窗包含L个采样点,短窗在长窗的内部前侧,长、短窗交界处为一个采样点,此点既存在于长窗内,又存在于短窗内,长窗与短窗组合成组合窗Wf.Wf的移动方式如图2所示,图中虚线为长窗,实线为短窗,共同组成组合窗Wf.

2) 以信号第1个点为起点,计算短窗与长窗内采样点的平均能量差Qf:

(1)

式中:E(ε)为采样点的能量值,数值上等于采样点数值的平方;Qf为t点的平均能量差.向后移动时窗,步长为一个采样点,计算下一个采样点的平均能量差值,遍历信号中所有采样点,得到Qf(t).Qf(t)的变化过程分为3个阶段:①初始点没有进入到组合窗内,所以Qf(t)的变化是噪声的变化,其幅值改变很小;②初始点到达组合窗内,长窗能量值逐渐增加,Qf(t)减小;③组合窗的起点与信号初始点重合,Qf(t)成为极小值,随组合窗向后移动,Qf(t)逐渐恢复至第一阶段的状态,Qf(t)与P波起点的关系如图3所示,可知Qf(t)的极小值点与心电P波初始点的对应关系.

图2 Wf移动方式示意图

图3 Qf(t)与P波起点的关系示意图

基础的长短时窗能量差法可以检测出波形的初始点,但由于算法本身限制,基础的长短时窗能量差法无法计算P波最后长为L的部分,导致P波的终点无法识别,这是算法中的缺点.本研究将组合窗的位置改变,将短窗放置在长窗内部后侧,长度不变,创造一组新的长短组合窗Wb,新的组合窗中长短窗终点重合,具体放置位置如图4所示.此组合窗内平均能量的计算过程与基础长短时窗能量差法类似,也是将短窗内采样点的平均能量减去长窗内采样点的平均能量,计算结果为Qb，如图5所示.此组合窗只可以从信号的第L个点开始计算,同样向后滑动,每次滑动一个采样点至P波结束.计算式为

(2)

图4 Wb移动方式示意图

图5 Qb(t)与P波终点的关系示意图

在Wb的移动过程中,经历起点时和基础长短时窗Wf不一样,Wb的短窗与长窗同时接触信号的起点,在此之前新长短时窗Wb内的平均能量差Qb是噪声的变化,其值几乎不变.在Wb接触信号初始点后,Qb的数值会增加,因为短窗大部分位于高能量区内,而长窗大部分位于低能量区内,随组合窗Wb向后滑动,短窗逐渐退出高能量区域,长窗依旧保持在高能量区内,Qb开始下降.当Wb的终点与P波截止点重合时,Qb达到极小值.再向后滑动长窗逐渐退出高能量区,Qb开始增加.

Qb可以将心电P波的截止点和起始点都标记出来,P波的起点即为Qb开始变化的起点,P波的终点即为Qb的极小值点.但Qb开始变化的起点不容易测量,而且,P波前可能存在一定量噪声,会干扰Qb开始变化的起点,导致检测结果不准确.

将Qb与Qf组合在一起,先将P波的预选区域通过Wf遍历一遍得到Qf(t),再使用新的组合窗Wb遍历一次P波预选区域得到Qb(t),最后将其结果融合，得Q(t)=Qb(t)-Qf(t)，将Wf和Wb的结果融合在一起.Qf(t)在P波起点时是极小值,Qb(t)在起点之后有上升趋势,所以P波的起点在Q(t)为波峰的极大值点;Qf(t)在P波终点时因公式与算法中的缺陷为0,而P波终点在Qb(t)中表示为波谷的极小值点.依据Q(t)的计算过程,P波的起点Q(t)表示为极大值点,P波的终点Q(t)表示为极小值点,这使P波的边界提取十分简单.

采集得到的心电信号通常会有不同程度的噪声,即使通过之前的去噪手段,也无法将噪声完全去除,因此,在所有的原波形上直接判断P波的边界是不现实的.依据能量变化速度将P波起点与终点检测的困难变换为检测波峰与波谷的问题,这不仅降低了检测P波的难度,还可以同时检测出正向、负向和双向3种不同形态的P波.基于被初步处理后的心电信号的噪声能量低及信号能量高的特点,将信号分为有信号的高能量区和没信号的低能量区.低能量区在Q(t)中体现为接近基线,高能量区在Q(t)中体现为较规律的波形,这样就可以据此来判断P波的起点与终点.

1.2 粒子群参数优化原理及具体步骤

粒子群算法具有进化计算和群体智能的特点,主要过程是通过个体间的协同与竞争,实现复杂空间中最优解的搜索[10].算法具体步骤如下:首先确定粒子的初始解和初速度.初始解是将粒子均匀分布在可选择解的范围内,每个粒子的初速度都是随机的.通过目标函数计算每个粒子的适应度,目标函数使用当前解提取P波的起点与终点.将结果与金标准对比,计算标准差,结果越小则参数越优.在解空间内选取适应度最小的解,作为当前迭代的全局最优解.每个粒子的结果与上一次迭代结果进行对比,选出粒子本身的历史最优解.当找到足够好的最优解或达到了最大迭代次数时,算法完成[11].迭代式为

Vid(t+1)=∂Vid(t)+C1rand(0,1)(Pid-Xid(t))+C2rand(0,1)(Pgd-Xid(t))，

(3)

Xid(t+1)=Xid(t)+Vid(t+1) .

(4)

式中:∂为惯性因子,其值为非负;Vid(t)和Xid(t)分别为第i个粒子在第t次迭代中第d维的速度和位置;C1和C2为加速度常数,通常设置为2;rand(0,1)为[0,1]区间上的随机数;Pid为第i个粒子中第d维度的个体最优解;Pgd为第d维度全局最优解[12].粒子维度指的是有多少个参数需要优化,在本文中粒子维度为2.根据实验经验将粒子数设为36,最大迭代次数为500.

使用手动标注的金标准(gn)与本文方法估计的位置(xn)之间的误差的方差作为粒子群算法的目标函数:

(5)

2 实验及结果

2.1 实验数据

使用两个数据库,分别为Lobachevsky University Electrocardiography Database(LUE)[13]和QT Database(QT)[14],分割出共含有6 299个心动周期的心电波形,进行P波提取算法的训练和验证.LUE中心电数据为12导联,采样频率为500 Hz,在此数据库中,P波的形态包括正向、负向和双向.QT数据库中心电数据包括2个导联,P波包括正向与负向.文献[5-6]的检测方法新颖,检测结果误差小.在LUE数据库中使用了1 892个含正向P波的单周期心电波形,1 397个含负向P波的单周期心电波形,1 074个含双向P波的单周期心电波形.从QT数据库中使用了1 936 个周期的单周期心电波形,将数据的70%作为训练集,30%作为测试集,具体如表1所示.

表1 实验数据集

2.2 实验结果

使用手动标注的金标准(gn)与本文方法估计的位置(xn)之间的误差标准差S来评估精度.S是所有标记出的起点及终点的标准差,pe是数据采样周期,St为以时间为单位的误差.金标准由高级医生和两个初级医生共同标注.误差式为

(6)

2.2.1 双边长短窗平均能量差法计算结果

图6～图8分别为3种不同形态的心电P波样本及其对应的Q(t),心电P波的起点与Q(t)的第一个波峰对应,P波终点与Q(t)最后一个波谷对应.

图6 正向P波检测结果

图7 负向P波检测结果

图8 双向P波检测结果

2.2.2 粒子群算法迭代结果

在粒子群算法中,粒子数为36,最大迭代次数为500,每次迭代选用全局最优解作为本研究的最优参数,图9为每次迭代所得最优参数计算得到的结果与金标准的误差标准差,本文算法在100次左右取得最优解,并保持稳定.

2.2.3 不同形态心电P波检测结果

在QT数据库中使用了584个心电周期对算法进行测试,数据均为正向心电,QT数据库中的数据采样频率为250 Hz,St为7.80 ms.本文算法检测结果与金标准之间的误差直方图如图10所示,在采样频率为250 Hz时,误差最大值是16 ms,大部分在[-8, 8]ms以内.

图9 粒子群迭代结果

图10 QT数据库P波检测结果

在LUE数据库中使用了1 299个心电周期对算法进行测试,综合误差标准差为6.26.其中,569个心电周期含有正向P波,St为6.18 ms;421个心电周期含有负向P波,St为6.21 ms;309个心电周期含有双向P波,St为6.52 ms.结果如图11～图13所示.

本文算法与文献[5-6]的检测结果对比如表2所示,P波起点和终点的检测结果皆优于其他两种算法.

表2 P波的起点与终点检测结果对比

图11 正向P波检测结果

图12 负向P波检测结果

图13 双向P波检测结果

图11～图13分别为正向P波、负向P波和双向P波的误差直方图,在采样频率为500 Hz时,误差最大值为12 ms,大部分在[-8, 8]ms.将本文算法与文献[5]中基于动态规划的参数混合高斯拟合法和文献[6]中无相位平稳小波变换法进行比较,所得统计结果如表3所示.由表3可知本文方法检测结果误差较小,可以检测出双向P波的位置.

表3 P波检测标准差

2.3 不同强度噪声下本算法的精度

将LUE数据库中心电分别加入信噪比为10,20,30 dB噪声检测本算法的鲁棒性.在正向P波与负向P波中加入30,20 dB的噪声,可以辨认出P波的形状.当噪声强度为10 dB时,由于P波的幅值较小,已无法从信号整体判断出P波的准确位置.在双向P波中加入30 dB的噪声时,可以辨别P波的大致形状,当噪声强度为20,10 dB时,噪声覆盖信号,已辨别不出P波的准确位置.本文使用P波位置的检测准确率(Acc)与检测P波误差的标准差S代表本算法对噪声的鲁棒性.

P波位置的准确率为

Acc=P/Pall.

(7)

式中:P为检测出P波的数量;Pall为每种形态P波的数量.

图14～图16分别为算法在正向、负向和双向P波加入仿真白噪声后检测结果误差的变化趋势.随机白噪声用来模拟心电信号采集过程中的高频随机噪声.该随机白噪声是均值为0、方差为1的高斯随机序列.双向P波在信噪比为10 dB噪声强度下准确率降低为0,正向与负向P波在信噪比为10 dB时仍可检测出P波位置;S与信噪比具有明显的负相关性,且变化趋势相似,在信噪比提升至20 dB时,S减小缓慢.

图14 噪声对算法检测正向P波的影响

图15 噪声对算法检测负向P波的影响

2.4 不同采样频率下本算法的精度

对LUE和QT数据进行重采样,将信号的采样频率提高到1 000,2 000,3 000,4 000,5 000 Hz,并重新计算结果.重采样后并没有对P波峰值检测的准确率造成影响,起点与终点的检测结果采样点误差有所提高,但时间误差逐渐降低.本文算法的检测结果误差随采样频率的升高不断减少,并趋于平稳,在采样频率为5 000 Hz时检测结果误差降为1.1 ms.重采样后不同方法检测结果对比如表4所示,采样频率对P波检测的误差随采样频率的增加而降低.本文方法在高采样频率的情况下,性能要优于文献[5-6].