鞠文博, 袁 林, 李 庆, 杨 波

(1.成都信息工程大学电子工程学院/物理场生物效应及仪器四川省高校重点实验室,四川 成都 610225;2.成都信息工程大学计算机学院,四川 成都 610225)

0 引言

胎心监护是正确评估胎儿宫内状态的重要检测手段,对及时诊断胎心失常、降低围产儿死亡率有极其重要的意义[1]。20 世纪中叶,临床使用胎儿心电图(fetal electrocardiogram,FECG)观察胎儿心脏信息。测量FECG 的最准确方法是将头皮电极直接置于孕妇体内进行胎儿心电图监测,得到的信号较为纯净,但这仅适用在分娩中监测,不能反复使用[2-4]。并且侵入式的监测会导致未确诊的胎儿缺氧,因此,无害的、非侵入式的监护方法显得尤为重要。腹部心电图技术直接将母体腹部电信号作为测量信号,与传统的胎心监护手段相比,这种方法操作简单,可以在孕后期使用,孕妇主观上更容易接收,能更好保护孕妇及胎儿。

非侵入式测量操作更加简便,但会受到其他生理信号和噪声的干扰,如:母体心电(maternal electrocardiogram,MECG)、基线漂移、肌电干扰和50Hz 工频干扰及其谐波干扰等。非侵入式方法检测得到母体心电信号幅值一般为0.5 ～2 mV,胎儿心电信号幅值约为10 ～50 μV,MECG 幅值远大于FECG 且频谱上也存在重叠,传统的滤波方法难以去除[5-6],因此在临床研究中去除母体心电干扰是重要也是必要的一步。为解决无创提取胎儿心电这一难题,人们相继提出了自适应滤波器；线性分解,如小波分解[7]、奇异值分解[8]、盲或半盲源分离[9]；非线性分解如子空间分离[10]等方法。其中盲源分离(blind source separation,BSS)引申出独立成分分析的快速不动点(FastICA)和特征矩阵的联合近似对角化(JADE)两种算法,JADE 在准确性方面优于FastICA,但在计算时间方面却不尽人意；而FastICA 可以通过许多微调步骤得到与JADE 相当的结果且受输入信号质量的影响较低[11],更加符合临床实验场景。两种方法的共性均是首先从AECG 信号中移除MECG 和增强FECG,再从剩余分量中提取FECG[12],而肌电干扰噪声会降低算法性能。

本文提出一种基于主成分分析(principal component analysis,PCA)和相关性检验多通道提取胎心率算法。首先通过样本熵对多通道原始信号进行质量评估,去除质量较差的通道信号；然后利用PCA 从评估后的混合信号中提取主成分得到最优AECG；其次,Pan-Tompkins 初步定位R 峰位置,并通过极性检测消除因波形倒置而出现的误检与漏检,对母体R 峰位置进行校正,利用相干平均法对母体R 峰位置序列进行重构,构造长度为两个R-R 间期的梯形窗,使之能包括一个完整的QRS 波群,对质量评估后满足条件的每路AECG 进行加窗处理,重采样至长度与重构后的母体QRS 波群相同,并在规定阈值范围内进行序列微调,以实现与MECG 最相关；最后,将AECG 与MECG相减得到的残差即为FECG,检测胎儿QRS 波群进而计算得到胎儿心率,通过分析相关性能指标结果,本方法适用于多通道腹部信号胎心电提取。

1 方法

1.1 数据集

数据来自于PhysioNet 提供的开源数据库Challenge 2013 Traning Set A。包括75 组胎儿心电图记录,每条记录包含4 路无创孕妇腹部电信号,均以每秒1000 样本点采样,并对每个胎儿的QRS 波群进行参考标注。

1.2 方法描述

算法流程如图1 所示。该方法由4 个主要步骤组成:(1)腹部心电信号预处理；(2)母体QRS 波群定位；(3)母体ECG 信号抵消；(4)胎儿QRS 波群检测。其中步骤(1)数据预处理阶段包括样滤波与信号质量评估；步骤(2)和步骤(3)可以拆分为主成分分析提取母体ECG、母体R 峰检测与校正以及母体ECG 信号重构；相关性检验校正胎儿QRS 波群在步骤(4)中进行。

图1 算法流程图

2 预处理

2.1 去噪

使用50 Hz 陷波滤波器消除工频干扰；基线漂移是由孕妇呼吸或运动引起的,使用巴特沃斯低通滤波器去除5 Hz以下的基线漂移干扰以及母体心电信号中的P、T波[13],实验证明去除P 波和T 波并不影响R 峰的检测；肌电属于高频干扰,使用小波去噪去除信号中的高频噪声[14]。腹部心电信号预处理前后对比如图2 所示。

图2 腹部心电预处理前后对比图

2.2 信号质量评估

信噪比(SNR)通常是作为评估信号质量好坏的指标之一,狭义上定义为输出信号的功率与同时输入的噪声功率之比。但检测信号为胎儿心电信号、母体心电信号与其他噪声的混合,且FECG 与MECG 频谱上存在重叠,因此使用信噪比无法准确判断各个通道信号质量的好坏。

近似熵(approximate entropy,ApEn)与样本熵(sample entropy,SampEn)是通过度量信号中产生新模式的概率大小来衡量时间序列复杂性[15],新模式产生的概率越大,序列的复杂性越高。熵值结果越低,序列的自我相似性就越低,信号就越简单；相反,熵值结果越高,信号就越复杂。

定义一个由N个数据组成的等样时间序列{x(n)}=x(1),x(2),…,x(n), 算法相关参数m,r,其中,m为整数,表示重构的维度；r为实数,表示“相似度”的度量值。样本熵的计算步骤如下:

(1)重构m维向量Xm(1),…Xm(N-m+1),其中Xm(i)=[x(i),x(i+1),…x(i+m-1)]T,1≤i≤N-m+1

(2)统计满足Xm(i)与Xm(j)之间的距离d小于等于r的j(1≤j≤N-m,j≠i)的数目,并记作Bi,则对于1≤i≤N-m定义:

其中,距离d由Xm(i)与Xm(j)对应元素中最大的差值的绝对值决定,即:

(3)求(r)对于所有i值的平均值,记为Bm(r),即

(4)增加维度到(m+1)维,重复步骤(2)将满足条件的个数记为Ai,则Ai(r)定义为

(5)定义Am(r)为

由于在实际计算应用过程中,N不可能为∞,因此当N取有限值时,样本熵估计为

式中,嵌入维数m一般取1 或2；相似度r的选择在很大程度上取决于实际应用场景,通常选择r=0.1·std ～0.25·std,其中std 表示原时间序列的标准差。在本次实验中取m=2,r=0.25·std。

设定一个恒定的样本熵阈值为1.5,每个通道的样本熵值与该阈值相比,将值小于1.5认定为质量良好的信号,用于后续分析；舍弃熵值大于1.5的通道信号,如果认定的质量良好的信号数量少于2 个,则保留熵值最小的2 个通道信号以供进一步分析[16]。a19 样本中4 路信号的熵值如表1 所示,4 路AECG 分别在10 s内的样本熵值如图3 所示。

表1 4 路信号样本熵值

图3 4 路AECG 分别在10 s 内的样本熵值

3 MECG 提取与重构

3.1 主成分分析

经过信号质量评估筛选后的AECG 全部用于PCA,自动保留了3 个成分,其各分量的方差贡献率分别为51.02%、47.15%和1.61%,这意味着几乎保留了原始信号的所有信息,且第一分量作为主分量包含信息最大。母体QRS 波群与胎儿QRS 波群在不同分量中的表现如图4 所示(红色框内为母体QRS 波群,黄色框内为胎儿QRS 波群)。在不失真的情况下,选择第一分量作为最优参考信号。

图4 不同分量下的MQRS 与FQRS

图5 校正前后R 峰位置对比

3.2 R 峰位置定位与校正

母体R 峰的检测方法主要用Pan-Tompkins 算法[17],相比于Q 波、P 波,R 波所在位置斜率最大。该方法首先对心电信号进行差分求导,然后进行平方,锐化信号。定义差分平方后的心电序列为X、长度为N、采样率fs,阈值thre 的计算方法为将心电按采样频率分为等长序列X(1),X(2),…,X(k),计算每个序列的均值记为m(1),m(2),…,m(k),取最小均值的0.3倍作为阈值,即:

thre=min{m(1),m(2),…,m(k)}×0.3,k=fs

在每一个长度为i的序列范围内,将大于阈值点中的最大值作为R 波的初步参考位置,然后拟合相同长度的抛物线用于增强QRS 波,便于寻找在以R 峰为中心220 ms范围内的最大值作为QRS 波群R 峰位置[18]。

3.3 MECG 重构

基于得到的母体QRS 波群,使用相干平均法构建每个MECG 模板[19]。首先,2 个连续R 峰间期构成完整的R-R 间期,使之能够包括一个心跳过程所产生的所有事件；然后,将每个R-R 间期重采样至恒定的1K采样率,以提供与AECG 平均处理；最后,将MECG 模板重新拉伸到每个周期的实际长度,再将它们逐拍连接,获得重构的MECG 信号。

4 胎儿心率检测

4.1 相关性检验校正

在MECG 重构之后,基于相关性检验来实现预处理后的AECG 和重构的MECG 信号之间的最佳匹配,实现过程如图6 所示。首先,将重构后的MECG 信号和预处理后的AECG 信号都使用原始采样频率的10倍进行重采样,即10 KHz上采样；其次,构建一个幅值为1,时间长度为440 ms的矩形窗,使其能包括一个完整的母体QRS 波段；然后,基于相关性检验执行窗内匹配操作,通过寻找最大相关系数找到MECG 和AECG 的最佳匹配位置。

图6 MECG 匹配抵消实现过程

匹配操作首先设定一个微调阈值,步长为10 个样本点,然后将MECG 信号向左或向右移动,计算MECG和AECG 的相关系数,以达到最大相关度,从而确保从AECG 中减去重构的母体心电模板后,可以获得更清晰的胎儿心电信号。图7 显示了对于来自记录a19 的2 个QRS 波段的相关性校正过程,并给出没有通过校正的例子。

图7 2 个QRS 波段相关性校正前后MECG 与原始AECG 信号对比

4.2 FECG 检测

从AECG 信号中抵消MECG 得到的残差中,MECG的残留成分已经很少,因此将残差信号作为最终胎儿心电结果。由于胎儿心跳周期为母体的1.5 ～2 倍,相比于检测母体QRS 波群使用的220 ms,胎儿QRS 波群恒定时间范围修改为150 ms[16],再用Pan-Tompkins 算法检测校正胎儿QRS 波群,得到胎儿R 峰序列。

5 结果

提取出数据集中的fQRS 文件中的专家标注信息作为参考,与实验结果得到的胎儿R 波位置计算的F1测量值与灵敏度Se 作为评估该算法精度的两项指标,并分别列出校正前后F1和Se 的值,如表2 所示。F1测量值与Se 计算方法如下:

表2 基于数据库Challenge 2013 Set A 相关性检验校正前后的各项指标结果对比 单位:%

其中,TP 是记录标注的与实验检测相匹配的QRS 波群数量,FP 是假阳性(误检的QRS 波群)数量,FN 是假阴性(漏检的QRS 波群)数量。考虑标注信息的误差,实验结果在50 ms范围内均是认为与标注结果匹配[20]。

6 结束语

提出基于主成分分析和相关性检验的腹部胎心点提取方法。首先对4 路原始信号进行预处理,通过样本熵值评估信号质量,符合条件的信号用于后续提取主成分；然后,将第一主成分作为最优参考信号,使用斜率法检测其QRS 波群以及R 峰位置,为增加R 峰位置的准确性,对原始信号的导数进行平方处理,初步得到的R 峰位置索引,再在区间内寻找幅值最大的位置,即R 峰校正后的位置；其次,利用相干平均法重构每个MECG 模板,在以每个母体R 峰位置为中心220 ms内的矩形窗内进行相关性检验,通过寻找最大相关系数找到MECG 和AECG 的最佳匹配位置；最后,从预处理后的每路AECG 中减去各自的MECG 信号,达到去除MECG 干扰的目的,得到的残差中包含了FECG 与少量噪声混合,再用主成分分析提取较为纯净的FECG。

本文对胎儿心电的提取思想仍然是目前流行的模板相减法,优点是能够以无监督的方式提供较好的检测结果。该方法改善了使用独立成分分析算法时肌电干扰噪声对提取的信号质量影响较大的问题,并在此基础上进行了相关性检验校正,极大消除了在模板抵消时由于时域不同步而出现的误差,提高了胎儿心率检测的准确性。