农汉彪, 曾巧妮
(百色学院,广西 百色 533000)
所采集的心电信号中包含有基线漂移成分。基线漂移是一种低频信号,而心电信号自身也含有十分丰富的低频成分,基漂漂移会叠加并掩盖有用的低频成分。基线漂移的存在会对后续分析、识别和诊断
在心电信号(ECG)采集的过程中,由于被试者的呼吸运动,测试电极与人体皮肤之间接触阻抗变化以及采集设备性能温度漂移等因素影响,会使得产生较大影响,为保证医学诊断的准确性,基线漂移应在心电信号预处理中予以消除。
近年来,国内外学者针对心电信号基线漂移的消除提出了很多新方法。林金朝等[1]提出了基于改进集合经验模态分解(EEMD)的消除心电信号基线漂移方法。YAO等[2]提出了基于完全集成经验模式分解自适应噪声算法(CEEMDAN)的方法,克服了传统经验模态分解(EMD)模态混叠的问题。刘春等[3]采用EEMD和经验小波变换(EWT)相结合的算法,Boda等[4]也提出了类似的利用EMD和EWT相混合方法从心电信号中抑制电力干扰(PLI)和基线漂移的方法。崔善政等[5]利用变分模态分解将心电图信号分解为一组模态分量,去除含有基线漂移成分的模态分量,重构剩余模态分量得到去除基线漂移后的心电图信号。Singhal等[6]提出了一种基于傅里叶分解法(FDM)的方法,从ECG信号中同时分离基线漂移和PLI,并获得干净的ECG数据。Romero等[7]提出了一种利用深度学习进行含有基线漂移的ECG信号滤波算法。
以上方法都趋向于将原信号进行分解,分别得到心电数据和干扰成分,而不可避免地存在频率混叠的问题,从而无法得到真正纯净的心电信号。另外,常用的算法还有中值滤波法、曲线拟合法[8]和形态学滤波[9]等。曲线拟合法通过拟合每一个心拍周期内适当数量的基点得到基线漂移干扰,基点即能够表征信号基线漂移走势的数据样本,该方法实时性、准确性好,但算法中基点的提取比较困难。针对这一问题,本文依据心电信号短时间内数据高阶统计量分布与基点之间关系提出了基于分段数据聚类提取基点从而实现基线漂移消除的方法。首先对数据进行分段并计算高阶统计量。基线漂移是缓变信号,其对短时间内分段数据的高阶统计量贡献趋于零。而当心电作用极小或为零时,数据只包含基线漂移和随机干扰,数据段的高阶统计量会趋向某些固定值(固定值和随机干扰相关)。接着对高阶统计量进行基于密度的聚类分析,可将高阶统计量趋于固定值的数据段筛选出来,筛选出数据段的中值即为基点。最后拟合基点即可得到基线漂移。
典型的心电信号如图1所示,每个周期主要由P波、QRS波群和T波组成,各个波的波形特征和之间的转换、过渡蕴含受试者的生理信息。其中,ST段是指QRS波群终点和T波起点之间所跨越的时间。S-T段期间,左右心室的肌细胞都处于兴奋期间,两者形成的综合电场在体表心电图中的贡献趋于零,导致S-T段心电信号大约处于基线的水平。另外,从上个周期的T波到下个周期的P波之间的过渡时间,也处于基线的水平。
图1 模拟心电信号波形
有上述心电信号的特点,可知连接(拟合)ST段和TP段的中值就可以得到信号的基线漂移。而由于信号的预处理优先于波段的检测与识别,即无法事先知道哪些数据处在这两个过渡段。根据信号的统计特性可以在数据采集或者处理过程中根据一个时间段内信号的统计特性对数据所在的区段进行区分。如图2所示为受到基线漂移和高斯随机噪声干扰后的心电信号,对数据进行分段处理并计算每段数据的中值、均值、方差、峭度等数学统计量。
图2 数据分段示意图
图2数据分段A-F,包含了6种可能出现的不同切分情况。其中D位于ST段,F位于TP段,假设每个区段的时间长度一致,每个区段数据的统计量如表1所示。从中值和均值上看,ABCDF 5个分段的数值都很相近而无法区分。从方差和峭度上看,D和F分段的统计值基本相同,明显区别于其他数据分段,因而可以通过分段数据的方差和峭度联合分布的密度情况对所有分段进行聚类分析,从而得到具有相同属性特征的数据分段,即可聚类出只包含基线漂移和随机噪声干扰的数据分段。
表1 分段数据统计量
在实测的心电信号中,由于受试者本身可能存在心血管疾病或其他心脏疾病,ST段和TP段并不一定是平直的,而相比其他波段,其分段数据的方差和峭度仍然存在差异,而同属于ST段和TP段的分段数据则具有相似的属性,因而同样可以通过聚类分析法方法区分识别出只包含基线漂移和随机噪声干扰的数据分段。
DBSCAN(density-based spatial clustering of applications with noise)聚类算法是一个很典型的基于密度的聚类算法,具有聚类速度快且能够有效处理噪声点和发现任意形状的空间聚类等显著优点。算法基于密度聚类的概念,要求聚类空间中的一定空间范围内所包含样本(数据点)的数目不小于某个给定阈值。算法的结果依赖于以下两个参数:
1)ε邻域:在一个样本周围邻近空间的半径;
2)minPts:邻近空间内至少包含样本个数。
若样本x1的ε邻域内至少包含minPts个样本,则x1是一个核心对象,若一个核心对象的邻域中包含了其他核心对象,则这些核心对象以及包含在它ε邻域内的所有样本构成一个类,如图3所示。即具有相同特征的对象样本之间的紧密相连的,在某类别任意样本周围一定空间范围内一定有同类别的样本存在。
图3 DBSCAN聚类算法原理示意图
依据心电信号的特性和DBSCAN聚类算法的优点,提出基于分段数据统计量聚类分析提取基准点的方法实现基线漂移消除,所提出的算法主要包含以下步骤:
1) 根据设定好的分段的大小(窗口大小)与数据重叠率(步进长度)对心电信号进行分段;设信号总数据点为N,分段数据的点数为L,步进长度为s(s≤L),则有分段数n=(N–L)/s(假设为整除的情况)。
2) 查找各个分段数据的中值与中值位置,统计分段数据的方差V,和峭度值K:
式中:L——分段数据长度;
——第j分段数据均值;
σj——第j分段数据标准差。
3) 分别将所有方差和峭度进行归一化处理以使得方差和峭度具有相同的距离观测测度。
4) 依据分段数据方差和峭度两个维度进行基于DBSCAN的聚类分析。
5) 合拼零点附近的几个分类(如果有多个分类),标记对应中值和中值位置为可用,标记其他分类和未分类的中值数据为不可用。
6) 拟合可用的中值和中值位置得到信号的基线漂移。
算法流程如图4所示。
图4 算法流程图
在数据分段过程前,需要确定分段窗口大小和步进长度两个超参数。分段数据的窗口大小直接影响着分段的数量以及提取的基线漂移的精细程度。分段数据重叠率的设置可以使得结果更加具有连续性,而如果重叠率过高则相当于移动中值滤波器。分段窗口大小,步进长度的设置与测试数据的成分结构有关,涉及数据的采样频率、主频(心率)和初始采集相位等因素。为了获取最优的分段大小和步进长度,可以采用已知数据穷举法进行实验。图5是在固定的采样频率(fs=360 Hz)且无分段数据重叠的前提下,心率、最优分段块长度之间的关系。
图5 心率与最优分块大小之间的关系
常用评价基线消除方法性能的参数有均方根误差(RMSE)、相关系数(COR)和信噪比(SNR)等。对于性能评估,RMSE值越小越好,COR和SNR越大越好。它们的计算式分别为:
式中:x——原始信号;
x′——处理后的信号;
Ps——信号的有效功率;
Pn——噪声的有效功率。
为了便于对消除方法进行分析与评估,以下将用模拟的心电信号和已知的基线漂移信号作为信号数据输入。模拟的心电信号波形如图1所示,将信号进行扩展,长度为10 000数据点,采样频率为360 Hz,信号峰值为2.0 mV,心率为60 bpm,见图6中的ECG信号。模拟的基线漂移信号由两个幅值为0.2 mV,频率分别为0.06 Hz和0.025 Hz的正弦信号相加组成,见图6中的BW信号。
图6 仿真待处理数据
为模拟实际测量时噪声干扰的情况,添加峰峰值为0.15 mV的高斯随机干扰,从而得到待处理的心电信号如图6中的ECG2P所示。
根据前述关于超参数的选择,选取了数据分块大小为160,步进长度为120。对数据进行分段后并查找各个数据段的中值得到中值数据如图7所示。
图7 仿真数据分段数据统计量
方差与峭度具有不同的示值范围,由于在后续基于密度的聚类分析中需要考虑分段数据统计量之间的距离测度,为了使得聚类分析时方差与峭度具有相同的权重,需要将所有数据段的方差和峭度进行归一化。归一化时可以将最大的方差和峭度作为1进行线性比例转化。即:
在聚类时,选用的算法是DBSCAN。在本文中,由于两个观测度已经进行了归一化处理,半径可以选用整个量值范围的1/20~1/10之间,最小集数可以选4~10之间,不同的半径和最小集数会得到不同的聚类结果。图8中是选用半径为0.05,最小集数为5时对仿真数据归一化方差与峭度聚类的结果。从中得到了3个分类,数据点相对集中没有出现独立的(未分类)数据点,由于仿真数据基线与噪声都是理想化的,所以每个分类都相对稠密。能代表基线成分的数据段应具有较小的方差和峭度,因而可以选用图中分类1和分类3所对应的数据段,舍弃分类2所对应的数据段。
图8 仿真分段数据方差-峭度聚类结果
选用分类1和分类3所对应数据段的中值进行样条曲线拟合即可得到信号中所含的基线漂移成分,结果如图9所示。
图9 仿真分段数据中值与基线漂移
为量化分析处理的效果,采用RMSE、COR和SNR三个评测参数对消除结果进行分析,并与常用基线漂移方法对同一数据进行对比分析。选用比较的方法为移动中值滤波,数学形态学滤波,小波滤波和经验模态分解(EMD)滤波。其中移动中值滤波选用窗口大小为160数据点(最优),步进长度为1数据点,中值滤波后通过0.5 Hz的低通滤波器;数学形态学滤波选用直线型结构元素,长度为55数据点(最优),数据分别经过开运算和闭运算后求和再通过0.5 Hz的低通滤波器得到平滑曲线;小波滤波采用离散Meyer小波对待处理数据进行8级(最优)分解,再重构逼近信息得到基线;经验模态分解(EMD)滤波采用基于CEEMDAN的滤波算法[2]得到。不同算法提取的基线漂移如图10所示。
图10 不同方法提取的基线漂移
图10 中本文方法、移动中值滤波、形态学滤波提取的基线漂移基本贴合目标曲线。而小波滤波和EMD滤波提取的基线漂移位于目标曲线上方,有基本固定的误差。不同方法在去除基线漂移效果评估参数详见表2。
表2 不同方法去除基线漂移效果参数
从评估参数RMSE可以看出,本文方法、移动中值滤波和形态学滤波三种方法提取的基线漂移与目标曲线相差很小,而小波滤波和EMD滤波相对较大,本文方法的均方误差最小,仅约为次好的中值滤波的31.5%。从COR参数上看,5种方法的消除基线漂移后的信号与原信号都具有非常高的相关系数。而从信噪比上看,小波滤波和EMD滤波表现比较差,本文的方法比其他最高的移动中值滤波高出约10 dB。在所有参数比较中,本文所提出方法的性能均优于其他4种方法。
为了验证基线漂移消除方法在实测信号的处理效果,选用了MIT-BIH数据中标记为117的这组数据进行分析。这组数据采样频率为360 Hz,心率约为50 bpm,数据中具有明显的基线漂移和随机干扰,数据的时域波形如图11所示。
图11 心电数据曲线MIT/BIH-117
根据数据的信息,处理的方法和仿真信号一样,分段窗口大小选用190,步进长度150。经过分段并统计得到各段的中值、方差、峭度如图12所示。
图12 分段数据统计量
实测分段数据的中值相比仿真数据没有明显基线漂移的曲线,而只是基线得基本走势。方差与峭度也相对分散。经归一化处理后,按方差和峭度两个观测维度进行DBSCAN聚类分析得到各个数据段分布如图13所示。
图13 实测数据方差-峭度聚类结果
聚类结果得到了4个分类和部分未分类的数据点。根据基线漂移的统计特征,选用分类1为可用的分段,未分类数据点和其他分类数据舍弃。采用样条曲线拟合可用的分类所对应的中值数据和中值位置即得到基线漂移成分。提取的基线漂移结果如图14所示。
图14 实测数据分段中值与基线漂移
将原信号减去所提取的基线漂移信号,即可得到处理结果,如图15所示。由于实测信号无法量化去除基线漂移效果,只能从视觉感官上进行判别和评价。从图中所示,通过所提出的方法处理后得到的数据,基本消除了基线漂移。
图15 消除基线漂移后的MIT/BIH-117数据
在实测信号分析处理过程中,由于其他干扰噪声的存在会在一定程度上影响信号基线漂移消除的效果,图16显示了不同消除方法在实测数据上处理结果的差异。本文方法能更好地贴近信号所包含的基线漂移,其他方法则会由于纯净心电信号相对基线漂移的单向性而浮在目标基线的上方。
图16 不同方法提取MIT/BIH-117基线漂移的结果
本文针对曲线拟合消除心电信号基线漂移方法中存在的提取拟合点困难的问题,提出了基于分段数据统计量聚类分析的方法来提取拟合点。文中通过一组仿真数据和一组实测数据对该方法进行了验证,从定量和定性两个方式验证方法的有效性和优越性。经数据实验验证,本文的方法较于常用方法具有明显优势,降低了均方根误差,减少信号的波形失真,提高了相关性和信噪比,从而得到纯净的心电信号,为后续分析、识别和诊断奠定了基础,保证信号数据的真实性。该方法中除了需要设置分段大小和步进大小,还需要设置聚类分析的半径和最小集合数,这些参数的设定需要了解待处理心电信号的采样频率,心率以及信号数据长度的大小,以便选择最优数值。该方法适用于原信号中存在只包含基线漂移或零点的时段,这些时段期间信号变化缓慢,被测对象未受到外界的作用力或者内部作用力相互抵消的情况,如心电信号,脉搏信号,轮轨力信号等。