CN102670188A - 一种用于单导胎儿心率检测的鲁棒自适应估计方法 - Google Patents

Abstract

本发明提出了一种用于单导胎儿心率检测的鲁棒自适应估计方法，其包括如下步骤：采集孕妇腹壁心电信号，并进行用于消除噪声和干扰的预处理；对孕妇腹壁心电信号中的母体心电信号进行建模、滤波；在滤波过程中进行自适应估计和鲁棒性分析；根据母体心电信号得到胎儿心电信号；得到胎儿心率。本发明在对母体心电信号进行滤波的过程中进行自适应估计和鲁棒性分析，提高了胎心率计算的可靠性和准确性。

Description

一种用于单导胎儿心率检测的鲁棒自适应估计方法

技术领域

本发明涉及胎儿心电检波方法，特别涉及一种用于单导胎儿心率检测的鲁棒自适应估计方法。

背景技术

心电监护是一项常用的医疗监护方式，在医院一般都配有心电监护设备，目前，心电监护设备正向便携式、家用化的方向发展，特别是家用胎儿心电监护设备是近年来的发展趋势。胎儿心电监护为胎儿的健康状况提供有用的参考信息，而胎儿心率是胎儿心电监护的基础项目之一，因此，如何准确获得胎儿心率对胎儿监护具有十分重要的价值。

作为一种无创监护手段，原始心电数据是通过在母体腹壁表面放置电极采集得到。这样采集到的腹壁心电信号(AECG)除了包含常见的各种肌电噪声，基线漂移，工频干扰、运动伪迹等，还含有母体心电信号(MECG)。胎儿心电信号(FECG)幅值一般在微伏级，而母体心电信号幅值较大，在毫伏级。所以，与MECG和各种噪声相比，胎儿心电信号显得十分微弱。如何在这样的大噪声背景下提取出胎儿心电信息，进行便携式监护是国内外学者研究的热点。

目前，采用较多的提取出胎儿心电信息的方法是母体心电抵消法，即首先预处理，然后提取出MECG，从腹壁混合信号中减去MECG，从而得到胎儿心电，再经过胎儿心电的R波检测，估计出胎儿心率。使用这种母体心电抵消法估计胎心率，最关键的一点是准确有效地估计MECG，只有准确的提取MECG，才能得到准确的胎心率。从腹壁混合信号中估计MECG，则胎儿心电和其他各种噪声都是干扰信号，问题就转换为成人的心电信号提取。目前，关于心电信号的提取，主要的方法有小波变换法、神经网络法、独立成分分析法(ICA)以及自适应滤波法等各种数字滤波。对这些方法研究分析，各自都有其优缺点。其中，小波变换法是一种时间-尺度分析方法，具有时频局域化分析能力和多分辨率分析的特点。利用小波变换多尺度、多分辨率的特性将噪声和ECG信号分解在不同尺度从而进行ECG信号和噪声的分离。虽然目前小波降噪法已经取得了一定的成果和应用，拟人脑机理和功能的计算机和人工智能处理技术，它具有非线性、高度并行性、鲁棒性和自主适应性等特点。但它在ECG信号中主要用于QRS检测中的自适应非线性预测，通过前面的采样点来预测当前的采样点。神经网络法的不足在于学习训练过程需要花费较长时间，而且要求训练样本具有一定的代表性，这在实际ECG降噪应用中很难实现。自适应滤波算法是针对信号和噪声频率变化而发展的自适应调节的滤波技术。自适应滤波一般采用最小均方误差准则作为代价函数，能够对参考输入信号进行自适应的加权处理，通过不断修正权向量，使滤波器输出逼近原始输入信号。虽然自适应滤波能够对ECG信号进行自适应处理，但这种自适应滤波方法要求额外提供参考信号，不适用于便携式的单通道的胎心率估计。

自从2003年P.E.McSharry提出了心电信号的动态模型之后，基于模型的ECG消噪方法成为近年来的研究热点。近几年，不少学者研究了利用动态ECG模型进行贝叶斯滤波提取ECG的方法，该方法在ECG降噪和FECG提取方面取得了一定的成果，但也还存在不少的问题和挑战，例如滤波不稳定，收敛性不好，有时甚至发散，而且滤波系统鲁棒性能不够高，容易受动态模型和噪声的影响。

发明内容

本发明旨在至少解决现有技术中存在的技术问题，特别创新地提出了一种用于胎儿心率检测的鲁棒自适应估计方法。

为了实现本发明的上述目的，本发明提供了一种用于胎儿心率检测的鲁棒自适应估计方法，其包括如下步骤：

S1：采集孕妇腹壁心电信号，并进行用于消除噪声和干扰的预处理；

S2：对孕妇腹壁心电信号中的母体心电信号进行建模、滤波；

S3：在滤波过程中进行自适应估计和鲁棒性分析；

S4：根据母体心电信号得到胎儿心电信号；

S5：得到胎儿心率。

本发明在对母体心电信号进行滤波的过程中进行自适应估计和鲁棒性分析，提高了胎心率计算的可靠性和准确性。

本发明的附加方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。

附图说明

本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解，其中：

图1是本发明用于胎儿心率检测的鲁棒自适应估计方法的流程图；

图2是本发明采用的ECG模型的三维图；

图3是图2中所示三维图转换到相位空间的图形。

具体实施方式

下面详细描述本发明的实施例，所述实施例的示例在附图中示出，其中自始至终相同或类似的标号表示相同或类似的元件或具有相同或类似功能的元件。下面通过参考附图描述的实施例是示例性的，仅用于解释本发明，而不能理解为对本发明的限制。

在本发明的描述中，除非另有规定和限定，需要说明的是，术语“安装”、“相连”、“连接”应做广义理解，例如，可以是机械连接或电连接，也可以是两个元件内部的连通，可以是直接相连，也可以通过中间媒介间接相连，对于本领域的普通技术人员而言，可以根据具体情况理解上述术语的具体含义。

本发明公开了一种用于单导胎儿心率检测的鲁棒自适应估计方法，如图1所示，其包括如下步骤：

S1：采集孕妇腹壁心电信号，并进行用于消除噪声和干扰的预处理；

S2：对孕妇腹壁心电信号中的母体心电信号进行建模、滤波；

S3：在滤波过程中进行自适应估计和鲁棒性分析；

S4：根据母体心电信号得到胎儿心电信号；

S5：得到胎儿心率。

本发明的目的是准确有效的估计MECG，从而获得胎心率，由于母体腹壁心电信号还包含其他噪声和干扰，因此在估计MECG之前，需要对其进行预处理。预处理包括去除高频噪声、基线漂移和工频干扰。因为母体心电频率主要在0.05Hz-100Hz之间，所以去除高频噪声可以采用一个低通滤波器实现。而50Hz或60Hz的工频干扰和基线漂移采用整系数滤波器去除，整系数滤波器采用递归算法，处理速度高。

本发明对母体心电信号建模的步骤为：首先采用非线性动态模型，其基本表达式由以下三个常微分方程组成：

$\left\{\begin{array}{c}\\ \stackrel{\·}{x}=\mathrm{\αx}-\mathrm{\ωy}\\ \stackrel{\·}{y}=\mathrm{\αy}+\mathrm{\ωx}\\ \stackrel{\·}{z}=-\underset{i\∈\{P,Q,R,S,T\}}{\mathrm{\Σ}}{a}_i\mathrm{\Δ}{\mathrm{\θ}}_i\exp (-\frac{\mathrm{\Δ}{{\mathrm{\θ}}_i}^2}{2{b_i}^2})-(z-z_0)\end{array}\right.---\left(1\right)$

其中，

a_i是高斯函数的幅值，

b_i是每个波形分散程度，

θ_i为波形的中心位置，

$\mathrm{\α}=1-\sqrt{x^2+y^2},$

为ECG信号的相位，

θ_i为波形的中心位置，

ω是在x-y平面上沿着圆周运动的角速度，映射到ECG信号中为瞬时心率，

z轴值以时间为参变量绘图描述出实际ECG的幅值变化，

z₀通过式子Z₀(t)＝Asin(2πf₂t)耦合到呼吸频率f₂，用以表征基线漂移的情况。根据现有的心电数据，将a_i，b_i和θ_i对应赋值得到ECG模型的三维图如图2所示。

由于以上模型是三维的，没有直观地对应实际的ECG信号，不利于医生参考。为了使得ECG模型更加直观，更适用于ECG信号的处理，需要在三个常微分方程的基础上进行改进。对一个ECG信号建模，首先在时域进行R波定位，确定每个R波的位置之后，将R-R间期映射到[-π，π]之间，这样就将ECG信号从时域转换到相位空间，一个ECG周期对应相位[0，2π]。这就得到简化改进后的ECG模型：

其中，k为时间序列值，

为ECG信号的相位，

θ_i为波形的中心位置，

Z_k为ECG信号的幅度值，

b_i为每个波形分散程度，

δ为采样周期，

η为模型误差或其他噪声源引起的加性随机白噪声。

ECG信号的幅度值Zk是以

为变量的函数，即幅度值Z_k在[0，2π]之间跟随

变化描绘出PQRST波的波形。(2)式中使用

i∈{P，Q，R，S，T}代替a_i，这样模型参数和ECG波形就一一映射。式中以采样周期δ将原动态模型离散化，根据上述方法对ECG信号建模有利于简化ECG信号估计框架。从三维模型图转换到相位空间，如图3所示。

本发明采用17状态变量的滤波框架，除了把

Zk作为状态变量之外，还把模型参数α_i，b_i，θ_i也当作状态变量进行估计，这样更能接近实际的心电数据，更具有灵活性。因此，状态变量可以写为

i∈{P，Q，R，S，T}，状态方程为：

对(3)式进行线性化处理，写为X_k+1＝F·X_k+W，其中，F为系统矩阵，即对状态方程求雅克比矩阵求得，W＝[w_a w_b w_θ w_ω η]^T为过程噪声，在本实施方式中，假设为高斯白噪声。

观测方程为：

由于观测量只有φ_k和Z_k，所以观测模型也可以表示为：Y_k＝H·X_k+V，其中 $H=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]$ 为观测矩阵， $V=\left[\begin{array}{c}{v}_1\\ v_2\end{array}\right]$ 为观测噪声。

本发明中孕妇腹壁心电信号包含了多种噪声和干扰，这些噪声和干扰因人因环境而异，不适合采用预先设定的噪声统计信息。本发明的滤波框架，在对MECG信号进行估计时，肌电噪声、胎儿心电和运动伪迹等都可以视为观测噪声。因此，在滤波的同时对观测噪声方差阵进行自适应估计，以提高滤波效率。本发明对观测噪声协方差阵估计采用移动窗口估计法。观测噪声协方差矩阵由观测残差估计得到，V_k表示观测残差向量，则其协方差矩阵

的估值可以由以下方法确定： ${\widehat{\mathrm{\Σ}}}_{V_k}=\frac{1}{N}\underset{j=0}{\overset{N}{\mathrm{\Σ}}}{V}_{k-j}{V}_{k-j}^T,$

其中，N为窗口宽度，则k时刻观测向量的协方差矩阵∑_k的估值为：滤波过程中不断地根据观测残差计算观测噪声协方差阵并不断更新，以达到自适应估计的目的。

另外，从母体腹壁心电信号中估计MECG，主要的噪声包含基线漂移，工频干扰，肌电噪声还有胎儿心电信号。其中，肌电噪声并不是完全表现为高斯型的噪声，它具有一定的非高斯性，在幅值较大处对MECG信号造成一定程度的冲击，胎儿心电信号的QRS波相对于MECG信号也表现为一个个冲击信号。所以肌电噪声和胎儿心电信号对MECG的估计有较严重的干扰，因此，本发明在滤波过程中引入鲁棒抗差思想，在本实施方式中，采用抗差最小二乘法，使得滤波过程能滤除严重的冲击干扰而有效地估计MECG。

抗差最小二乘估计的关键是建立合适的权函数，在保证较强的抗差性的同时又有较高的效率。抗差最小二乘估计中有多种极值函数和权因子方案，本发明中选择的是经典的Huber法。等价权由原观测权和权因子决定，而权因子由极值函数ρ(v)规定，Huber法规定的ρ(v)函数为：

$\mathrm{\ρ}\left(v\right)=\left\{\begin{array}{cc}{v}^2/2,& \left|v\right|\≤\mathrm{\Δ}\\ \mathrm{\Δ}\left|v\right|-\frac{1}{2}{k}^2,& \left|v\right|>\mathrm{\Δ}\end{array}\right.,$ 权因子 $\mathrm{\ω}\left(v\right)=\left\{\begin{array}{c}1,\left|v\right|\≤\mathrm{\Δ}\\ \frac{\mathrm{\Δ}}{\left|v\right|},\left|v\right|>\mathrm{\Δ}\end{array}\right.,$ 这里v为归一化的观测残差，Δ为调整系数。

本发明针对腹壁心电信号中观测噪声的冲击特性，利用鲁棒的思想来估计

这样既可以实现自适应估计的目的，又可以达到鲁棒滤波的效果，具体算法如下：

$\hat{d}\left(k\right)=\mathrm{median}\frac{|V\left(j\right)-\mathrm{median}\left(V\left(j\right)\right)|}{0.6745},$ j＝k-N+1，…，k，则 ${\widehat{\mathrm{\Σ}}}_{V_k}={\hat{d}}^2\left(k\right)$

利用观测残差的中位数估计其标准差，从而求得估计方差矩阵

本发明研究的对象是孕妇腹壁心电信号，胎儿心电信号和肌电噪声都可以视为带有冲击特性的观测噪声。因此，采用中位数抗差估计算法计算观测残差的方差阵有利于腹壁心电信号中MECG的鲁棒估计，而且由于在滤波过程对观测噪声进行了估计和修正，该方法在提高滤波过程的鲁棒性的同时也满足了自适应滤波的要求。

在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

尽管已经示出和描述了本发明的实施例，本领域的普通技术人员可以理解：在不脱离本发明的原理和宗旨的情况下可以对这些实施例进行多种变化、修改、替换和变型，本发明的范围由权利要求及其等同物限定。

Claims (6)

1.一种用于单导胎儿心率检测的鲁棒自适应估计方法，其特征在于，包括如下步骤：

S1：采集孕妇腹壁心电信号，并进行用于消除噪声和干扰的预处理；

S2：对孕妇腹壁心电信号中的母体心电信号进行建模、滤波；

S3：在滤波过程中进行自适应估计和鲁棒性分析；

S4：根据母体心电信号得到胎儿心电信号；

S5：得到胎儿心率。

2.如权利要求1所述的用于胎儿心率检测的鲁棒自适应估计方法，其特征在于，所述预处理包括去除高频噪声、基线漂移和工频干扰，所述高频噪声采用低通滤波器去除，所述工频干扰和基线漂移采用整系数滤波器去除。

3.如权利要求1所述的用于胎儿心率检测的鲁棒自适应估计方法，其特征在于，对所述母体心电信号建模的步骤为：

S21：采用非线性动态模型建立方程：

$\left\{\begin{array}{c}\\ \stackrel{\·}{x}=\mathrm{\αx}-\mathrm{\ωy}\\ \stackrel{\·}{y}=\mathrm{\αy}+\mathrm{\ωx}\\ \stackrel{\·}{z}=-\underset{i\∈\{P,Q,R,S,T\}}{\mathrm{\Σ}}{a}_i\mathrm{\Δ}{\mathrm{\θ}}_i\exp (-\frac{\mathrm{\Δ}{{\mathrm{\θ}}_i}^2}{2{b_i}^2})-(z-z_0)\end{array}\right.$

其中，

$\mathrm{\α}=1-\sqrt{x^2+y^2},$

ω是在x-y平面上沿着圆周运动的角速度，

a_i是高斯函数的幅值，

b_i是每个波形分散程度，

为ECG信号的相位，

θ_i为波形的中心位置，

z₀通过式子Z₀(t)＝Asin(2πf₂t)耦合到呼吸频率f₂，用以表征基线漂移的情况；

S22：将母体心电信号从时域转换到相位空间，

其中，k为时间序列值，

为ECG信号的相位，

Z_k为ECG信号的幅度值，

δ为采样周期，

η为模型误差或其他噪声源引起的加性随机白噪声。

4.如权利要求1所述的用于胎儿心率检测的鲁棒自适应估计方法，其特征在于，所述滤波方法为：将Z_k，模型参数α_i，b_i，θ_i当作状态变量进行估计，状态变量为：

i∈{P，Q，R，S，T}。

则状态方程为：

进行线性化处理，写为X_k+1＝F·X_k+W，

其中，F为系统矩阵，W＝[w_α w_b w_θ w_ω η]^T为过程噪声。

5.如权利要求1所述的用于单导胎儿心率检测的鲁棒自适应估计方法，其特征在于，所述自适应估计方法为移动窗口估计法。

6.如权利要求1所述的用于胎儿心率检测的鲁棒自适应估计方法，其特征在于，所述鲁棒性分析方法为抗差最小二乘法。

Images