CN105342594B - 用于家庭监护的单通道母体腹壁胎儿心率鲁棒估计方法 - Google Patents
用于家庭监护的单通道母体腹壁胎儿心率鲁棒估计方法 Download PDFInfo
- Publication number
- CN105342594B CN105342594B CN201510883842.0A CN201510883842A CN105342594B CN 105342594 B CN105342594 B CN 105342594B CN 201510883842 A CN201510883842 A CN 201510883842A CN 105342594 B CN105342594 B CN 105342594B
- Authority
- CN
- China
- Prior art keywords
- msub
- mrow
- msubsup
- fetal
- mtd
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Expired - Fee Related
Links
- 238000000034 method Methods 0.000 title claims abstract description 57
- 210000002458 fetal heart Anatomy 0.000 title claims abstract description 52
- 238000012544 monitoring process Methods 0.000 title claims abstract description 31
- 230000008774 maternal effect Effects 0.000 title claims abstract description 18
- 230000003187 abdominal effect Effects 0.000 title claims abstract description 15
- 230000001605 fetal effect Effects 0.000 claims abstract description 126
- 210000002784 stomach Anatomy 0.000 claims abstract description 56
- 210000003754 fetus Anatomy 0.000 claims description 44
- 238000004422 calculation algorithm Methods 0.000 claims description 20
- 238000000605 extraction Methods 0.000 claims description 13
- 230000009897 systematic effect Effects 0.000 claims description 13
- 238000002565 electrocardiography Methods 0.000 claims description 9
- 238000000354 decomposition reaction Methods 0.000 claims description 8
- 238000001914 filtration Methods 0.000 claims description 8
- 239000011159 matrix material Substances 0.000 claims description 8
- 238000000718 qrs complex Methods 0.000 claims description 8
- 238000001514 detection method Methods 0.000 claims description 6
- IAPHXJRHXBQDQJ-ODLOZXJASA-N jacobine Natural products O=C1[C@@]2([C@H](C)O2)C[C@H](C)[C@](O)(C)C(=O)OCC=2[C@H]3N(CC=2)CC[C@H]3O1 IAPHXJRHXBQDQJ-ODLOZXJASA-N 0.000 claims description 3
- 238000005070 sampling Methods 0.000 claims description 3
- 230000005611 electricity Effects 0.000 claims description 2
- 241000208340 Araliaceae Species 0.000 claims 1
- 235000005035 Panax pseudoginseng ssp. pseudoginseng Nutrition 0.000 claims 1
- 235000003140 Panax quinquefolius Nutrition 0.000 claims 1
- 235000008434 ginseng Nutrition 0.000 claims 1
- 210000004761 scalp Anatomy 0.000 description 30
- 230000008859 change Effects 0.000 description 10
- 238000010586 diagram Methods 0.000 description 10
- 238000012935 Averaging Methods 0.000 description 6
- 238000000926 separation method Methods 0.000 description 5
- 210000001015 abdomen Anatomy 0.000 description 4
- 238000005562 fading Methods 0.000 description 4
- 238000012545 processing Methods 0.000 description 4
- 230000001360 synchronised effect Effects 0.000 description 4
- 230000003044 adaptive effect Effects 0.000 description 3
- 230000036541 health Effects 0.000 description 3
- 230000008569 process Effects 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 201000004569 Blindness Diseases 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 2
- 238000004070 electrodeposition Methods 0.000 description 2
- 238000007689 inspection Methods 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- 238000005312 nonlinear dynamic Methods 0.000 description 2
- 230000009984 peri-natal effect Effects 0.000 description 2
- 238000011158 quantitative evaluation Methods 0.000 description 2
- 230000011218 segmentation Effects 0.000 description 2
- 230000001568 sexual effect Effects 0.000 description 2
- 210000004291 uterus Anatomy 0.000 description 2
- 208000036818 High risk pregnancy Diseases 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 231100001261 hazardous Toxicity 0.000 description 1
- 210000003128 head Anatomy 0.000 description 1
- 210000001161 mammalian embryo Anatomy 0.000 description 1
- 238000010197 meta-analysis Methods 0.000 description 1
- 210000003205 muscle Anatomy 0.000 description 1
- 230000035935 pregnancy Effects 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 210000001519 tissue Anatomy 0.000 description 1
- 238000002604 ultrasonography Methods 0.000 description 1
Landscapes
- Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)
Abstract
本发明提供了一种用于家庭监护的单通道母体腹壁胎儿心率鲁棒估计方法，将母体心电作为胎儿心电的一部分进行整体建模，从波形形态角度建立腹壁胎儿心电信号状态空间模型，运用强跟踪滤波方法实现胎儿心电R波与母体心电波形的递推估计，并对估计的胎儿心电R波波形进行峰值定位，从而获取胎儿瞬时心率，本方法适用于在计算、存储能力受限的本地监护终端上直接提取胎儿心率，且操作简单，孕妇完全可以独立完成监护过程，提升了用户体验。
Description
技术领域
本发明涉及胎儿监护领域，具体涉及一种用于家庭监护的单通道母体腹壁胎儿心率鲁棒估计方法。
背景技术
围产期，孕妇需要定期到医院进行例行性胎儿健康检查。对于高危妊娠孕妇，到医院检查次数则更为频繁。这类孕妇大多出行不便，而且作为人口大国，我国产科门诊病人众多。因此，研制胎儿心率家庭监护系统，让孕妇可以足不出户实现胎儿的动态监护，以缓解产科门诊医师的压力，降低高危患者出行风险，是当前母婴监护领域研究的重大课题。
目前，基于多普勒超声的胎儿心率(fetal heart rate,FHR)监护是临床上胎儿健康检查的基本技术手段。然而，基于超声探测的胎心率监护方式对胎位十分敏感，必须将探头实时对准胎儿心脏瓣膜，通常需要专业医师进行操作。并且这种方式属于有源探测，对胎儿健康存在潜在的危害，一次连续监护时间通常不超过半小时。上述因素，限制了多普勒超声胎心率监护在远程家庭监护中的应用。
从孕妇腹壁采集胎儿心脏电信号(fetal electro-cardiogram,FECG)实现胎儿心率监护，由于完全无创，操作简单，可在整个围产期随时进行，是最适用于远程家用的胎儿无创监护方法。但是，胎儿心电通过经过多种肌层组织传输至体表，致使腹壁采集的FECG信噪比极低，其频谱与母体心电信号(Maternal ECG,MECG)、电极接触噪声、肌电噪声等干扰信号大部重叠，分离提取胎儿心电图算法异常复杂。这导致目前的胎儿远程监护还需要将腹壁胎儿心电图传输到监护中心进行分析与处理，再将提取的胎儿心率回传至客户端。整个监护过程，需要实时接入通信网络，系统实现复杂，资源消耗大，成本高。并且孕妇完全暴露在电磁环境中，用户体验差，不易被用户所接受。
从母体腹壁心电中分离提取出微弱的胎儿心电成分是电生理信号处理领域的经典课题，经历了从最初的自适应噪声抵消到以奇异值分解(SVD)即主元分析(PCA)、独立元分析(ICA)为代表的“全盲”分离，再到最近提出的以约束ICA、周期元分析(PICA)为代表的“半盲”分离的研究过程。由于腹壁FECG本身幅度微弱(比母体干扰低10至30倍)且逐搏间动态变化(包括心率、波形形态)，一旦胎儿心电被母体QRS波群完全覆盖，现有分离方法通常表现出很差的鲁棒性，常常导致胎儿心电的严重衰减，甚至完全消失。同时，现有的提取方法需要进行采集通道的冗余配置，导致操作复杂，孕妇需在家人协作下才能完成胎儿监护过程。
发明内容
本申请通过提供一种用于家庭监护的单通道母体腹壁胎儿心率鲁棒估计方法，以解决现有技术中远程家用胎儿无创监护系统实现复杂，资源消耗大，鲁棒性差以及操作复杂而导致的孕妇无法独立完成监护的技术问题。
为解决上述技术问题，本申请采用以下技术方案予以实现：
一种用于家庭监护的单通道母体腹壁胎儿心率鲁棒估计方法，包括如下步骤：
S1：信号预处理：对腹壁胎儿的心电信号进行预处理，以去除心电信号中的基线漂移和工频干扰；
S2：构建单通道腹壁胎儿心电信号的状态空间模型，通过系统参数初始化得到模型所需相关参数的初始值；
S3：强跟踪波形估计：根据胎儿心电信号的状态空间模型及提取的模型参数，运用强跟踪滤波方法，以实现胎儿心电R波与母体心电各个波形的递推估计；
S4：胎儿R波定位：采用递归Kalman滤波方法，实现R波动态定位，从而获取胎儿瞬时心率；
其中，步骤S2中构建单通道腹壁胎儿心电信号的状态空间模型的方法如下：
首先，将母体心电作为胎儿心电的一部分，构建一种六状态的胎儿心电信号动态模型：
式中，为k时刻胎儿心电的瞬时相位，为k时刻胎儿心电的瞬时心率，为k时刻母体心电的瞬时相位，为k时刻母体心电的瞬时心率，δ为采样时间间隔，k+1时刻胎儿心电R波用一个高斯函数进行建模，表示为Rk+1，分别为胎儿心电R波对应高斯核函数的幅度、宽度和中心位置，k+1时刻母体心电P波和T波分别用两个高斯核函数进行建模，表示为Pk+1和Tk+1，k+1时刻母体心电QRS波群采用三个高斯核函数进行建模，表示为Ck+1，αi,k和bi,k，i∈{P1,P2,Q,R,S,T1,T2}分别为母体心电各个波形对应高斯核函数的幅度和宽度，分别为母体心电各个波形对应高斯核函数的中心位置，为高斯噪声；
由所述胎儿心电信号动态模型得到系统的观测模型为：
式中，为胎儿心电相位量测噪声，为母体心电相位量测噪声，Sk为胎儿心电信号测量样本，υS,k为量测噪声；
令Xk表示状态矢量，Yk表示观测向量，Uk表示系统噪声，Vk表示量测噪声，函数f和g分别描述系统方程和观测方程，则腹壁胎儿心电信号的状态空间模型为：
本发明以胎儿心率估计为目标，只需跟踪胎儿R波的位置，因此，仅就胎儿心电R波建模，以单一的高斯核函数描述，然而对于腹壁胎儿心电图，受胎位、胎龄、子宫肌电及电极配置等诸多因素的影响，母体心电波形逐拍间通常存在较大的变化，尤其是P波、T波。为了更好的跟踪母体心电波形的动态变化，本发明将其分为P波、QRS复波及T波三个部分分别进行建模，并扩展到状态向量中，以此约束相应的波形参数，从而构建一种六状态的胎儿心电信号动态模型。
进一步地，系统参数初始化是进行波形估计的前提，即确定状态矢量的处置与方差以及系统噪声与观测噪声的均值与方差，步骤S2中系统参数初始化包括母体参数初始化和胎儿参数初始化，其中，母体参数初始化具体包括：对预处理之后的心电信号进行R波探测，获得母体瞬时心率的均值与方差，同时将心电信号从时域转换为相位域，从而获得母体心电的瞬时相位运用单通道奇异值分解方法获取母体平均心拍及平均心拍的高斯核参数；
胎儿参数初始化具体包括：从预处理之后的心电信号中去除提取的母体心电成分，获得胎儿心电初步估计，再进行R波探测，获得胎儿心率的均值与方差，同时将心电信号从时域转换为相位域，从而获得胎儿心电的瞬时相位运用单通道奇异值分解方法获取胎儿平均心拍及平均心拍的高斯核参数。
由于单通道电极配置方案需要随胎位情况而变化，而不同电极位置获取的胎儿心电参数差异很大。通常情况下，每次进行单通道胎心率监护，都需要进行系统参数的初始化。但由于参数初始化时间短，通常只需要大约10个心拍，所以只在开始阶段存在10秒左右的时延，电极配置稳定后，即可实现在线估计，算法实时性不受影响。因此，作为一种优选的技术方案，步骤S22中系统参数初始化过程中选取10秒的心电数据。
从腹壁胎儿心电信号的状态空间模型可见，系统动力学模型具有明显的非线性。当上述系统模型具有足够的精度，扩展Kalman滤波(Extended Kalman Filtering，EKF)可以给出比较准确的状态估计值。因此，在成人心电波形估计、微小特征提取及分段中，EKF估计方法获得了深入的研究与应用。然而，在腹壁心电处理中，由于胎儿宫内旋转等引起腹壁胎儿心电模型结构的摄动，导致胎儿心电信号动态模型与所描述的非线性腹壁心电系统不能完全匹配。加之系统外部干扰噪声子宫收缩引起的肌电噪声、母体随机性动作引起的电极接触噪声等非平稳。这就要求滤波系统具备很强的自适应跟踪能力，而扩展Kalman滤波框架之下的心电信号动态估计方法鲁棒性差。一方面，EKF对观测模型不确定性带来的粗差干扰鲁棒性很差；另一方面，当系统达到稳态时,其增益趋于极小值,此时将丧失对突变状态的跟踪能力，难以处理系统模型自身的摄动问题；同时扩展Kalman滤波算法对初始参数的依赖性也较大。
针对非线性动态估计中系统模型的不确定性问题，本发明采用一种强跟踪滤波方法(Strong Tracking Filter,STF)，该方法通过引进时变的渐消因子,动态调节增益矩阵,使得滤波器对模型不确定性具有较好的鲁棒性而且计算复杂度与EKF相当，从而提高对心电波形的动态跟踪能力。
进一步地，步骤S3中强跟踪滤波方法具体为：
对于胎儿心电信号的状态空间模型的非线性状态方程，扩展Kalman滤波递归估计算法为：
式中，为状态矢量的先验估计值，为状态矢量的后验估计值，Kk为可调增益阵，为先验估计的协方差矩阵，为后验估计的协方差矩阵，Rk为观测噪声方差，Qk为状态噪声方差，rk为更新后的信息，即新息，Ak、Bk、Ck、Dk为胎儿心电信号状态空间模型的Jacobin矩阵，
在扩展Kalman滤波递归估计算法上，强迫不同时刻的新息矢量rk正交，即同时满足以下两个方程：
求得满足上式的可调增益阵Kk，即可得到强跟踪滤波方法。
与现有技术相比，本申请提供的技术方案，具有的技术效果或优点是：实现了单通道胎儿心率的鲁棒在线估计，适用于在计算、存储能力受限的本地监护终端上直接提取胎儿心率，同时保证服务质量；另一方面，单通道采集的电极配置操作简单，孕妇完全可以独立完成监护过程，提升了用户体验。
附图说明
图1为本发明的胎儿心率估计流程图；
图2为原始腹壁胎儿心电图与头皮胎儿心电图对比图；
图3为母体心电图与头皮胎儿心电图EKF波形估计对比图；
图4为母体心电图与头皮胎儿心电图强跟踪波形估计对比图；
图5为基于强跟踪滤波算法的头皮胎儿心电R波定位图与腹壁胎儿心电R波定位图对比图；
图6为基于强跟踪滤波算法的头皮胎儿心电瞬时心率图与腹壁胎儿心电瞬时心率图对比图；
图7为基于强跟踪滤波算法的腹壁胎儿瞬时心率图与头皮胎儿瞬时心率图；
图8为基于强跟踪滤波算法的腹壁胎儿瞬时心率与头皮胎儿瞬时心率对比图；
图9为基于强跟踪滤波算法的腹壁胎儿瞬时心率与头皮胎儿瞬时心率相关性图；
图10为基于强跟踪滤波算法的腹壁胎儿心率与头皮胎儿心率分段平均对比图；
图11为基于强跟踪滤波算法的腹壁胎儿心率与头皮胎儿心率分段平均误差对比图。
具体实施方式
本申请实施例通过提供一种用于家庭监护的单通道母体腹壁胎儿心率鲁棒估计方法，以解决现有技术中远程家用胎儿无创监护系统实现复杂，资源消耗大，鲁棒性差以及操作复杂而导致的孕妇无法独立完成监护的技术问题。
为了更好的理解上述技术方案，下面将结合说明书附图以及具体的实施方式，对上述技术方案进行详细的说明。
实施例
一种用于家庭监护的单通道母体腹壁胎儿心率鲁棒估计方法，如图1所示，包括如下步骤：
S1：信号预处理：对腹壁胎儿的心电信号进行预处理，以去除心电信号中的基线漂移和工频干扰；
S2：构建单通道腹壁胎儿心电信号的状态空间模型，通过系统参数初始化得到模型所需相关参数的初始值；
S3：强跟踪波形估计：根据胎儿心电信号的状态空间模型及提取的模型参数，运用强跟踪滤波方法，以实现胎儿心电R波与母体心电各个波形的递推估计；
S4：胎儿R波定位：采用递归Kalman滤波方法，实现R波动态定位，从而获取胎儿瞬时心率；
其中，步骤S2中构建单通道腹壁胎儿心电信号的状态空间模型的方法如下：
本发明以胎儿心率估计为目标，只需跟踪胎儿R波的位置，因此，仅就胎儿心电R波建模，以单一的高斯核函数描述，然而对于腹壁胎儿心电图，受胎位、胎龄、子宫肌电及电极配置等诸多因素的影响，母体心电波形逐拍间通常存在较大的变化，尤其是P波、T波。为了更好的跟踪母体心电波形的动态变化，本发明将其分为P波、QRS复波及T波三个部分分别进行建模，并扩展到状态向量中，以此约束相应的波形参数，从而构建一种六状态的胎儿心电信号动态模型：
式中，为k时刻胎儿心电的瞬时相位，为k时刻胎儿心电的瞬时心率，为k时刻母体心电的瞬时相位，为k时刻母体心电的瞬时心率，δ为采样时间间隔，k+1时刻胎儿心电R波用一个高斯函数进行建模，表示为Rk+1，分别为胎儿心电R波对应高斯核函数的幅度、宽度和中心位置，k+1时刻母体心电P波和T波分别用两个高斯核函数进行建模，表示为Pk+1和Tk+1，k+1时刻母体心电QRS波群采用三个高斯核函数进行建模，表示为Ck+1，αi,k和bi,k，i∈{P1,P2,Q,R,S,T1,T2}分别为母体心电各个波形对应高斯核函数的幅度和宽度，分别为母体心电各个波形对应高斯核函数的中心位置，为高斯噪声；
由所述胎儿心电信号动态模型得到系统的观测模型为：
式中，为胎儿心电相位量测噪声，为母体心电相位量测噪声，Sk为胎儿心电信号测量样本，υS,k为量测噪声；
可见，系统的观测模型是线性模型，而动态模型为非线性模型，令Xk表示状态矢量，Yk表示观测向量，Uk表示系统噪声，Vk表示量测噪声，函数f和g分别描述系统方程和观测方程，则腹壁胎儿心电信号的状态空间模型为：
进一步地，系统参数初始化是进行波形估计的前提，即确定状态矢量的处置与方差以及系统噪声与观测噪声的均值与方差，考虑到参数初始化算法的简洁性，且初始参数精度要求不高，无需准确分离出所有心拍，因此，本实施例中采用单通道奇异值分解SVD方法，提取母体心电与胎儿心电，获取其平均模型，采用静态高斯函数建模方法，从而获得高斯核参数估计值，步骤S2中系统参数初始化包括母体参数初始化和胎儿参数初始化，其中，母体参数初始化具体包括：对预处理之后的心电信号进行R波探测，获得母体瞬时心率的均值与方差，同时将心电信号从时域转换为相位域，从而获得母体心电的瞬时相位运用单通道奇异值分解方法获取母体平均心拍及平均心拍的高斯核参数；
胎儿参数初始化具体包括：从预处理之后的心电信号中去除提取的母体心电成分，获得胎儿心电初步估计，再进行R波探测，获得胎儿心率的均值与方差，同时将心电信号从时域转换为相位域，从而获得胎儿心电的瞬时相位运用单通道奇异值分解方法获取胎儿平均心拍及平均心拍的高斯核参数。
由于单通道电极配置方案需要随胎位情况而变化，而不同电极位置获取的胎儿心电参数差异很大。通常情况下，每次进行单通道胎心率监护，都需要进行系统参数的初始化。但由于参数初始化时间短，通常只需要大约10个心拍，所以只在开始阶段存在10秒左右的时延，电极配置稳定后，即可实现在线估计，算法实时性不受影响。因此，作为一种优选的技术方案，步骤S22中系统参数初始化过程中选取10秒的心电数据。
从腹壁胎儿心电信号的状态空间模型可见，系统动力学模型具有明显的非线性。当上述系统模型具有足够的精度，扩展Kalman滤波(Extended Kalman Filtering，EKF)可以给出比较准确的状态估计值。因此，在成人心电波形估计、微小特征提取及分段中，EKF估计方法获得了深入的研究与应用。然而，在腹壁心电处理中，由于胎儿宫内旋转等引起腹壁胎儿心电模型结构的摄动，导致胎儿心电信号动态模型与所描述的非线性腹壁心电系统不能完全匹配。加之系统外部干扰噪声子宫收缩引起的肌电噪声、母体随机性动作引起的电极接触噪声等非平稳。这就要求滤波系统具备很强的自适应跟踪能力，而扩展Kalman滤波框架之下的心电信号动态估计方法鲁棒性差。一方面，EKF对观测模型不确定性带来的粗差干扰鲁棒性很差；另一方面，当系统达到稳态时,其增益趋于极小值,此时将丧失对突变状态的跟踪能力，难以处理系统模型自身的摄动问题；同时扩展Kalman滤波算法对初始参数的依赖性也较大。
针对非线性动态估计中系统模型的不确定性问题，本发明采用一种强跟踪滤波方法(Strong Tracking Filter,STF)，该方法通过引进时变的渐消因子,动态调节增益矩阵,使得滤波器对模型不确定性具有较好的鲁棒性而且计算复杂度与EKF相当，从而提高对心电波形的动态跟踪能力。
进一步地，步骤S3中强跟踪滤波方法具体为：
对于胎儿心电信号的状态空间模型的非线性状态方程，扩展Kalman滤波递归估计算法为：
式中，为状态矢量的先验估计值，为状态矢量的后验估计值，Kk为可调增益阵，为先验估计的协方差矩阵，为后验估计的协方差矩阵，Rk为观测噪声方差，Qk为状态噪声方差，rk为更新后的信息，即新息，Ak、Bk、Ck、Dk为胎儿心电信号状态空间模型的Jacobin矩阵，
在扩展Kalman滤波递归估计算法上，强迫不同时刻的新息矢量rk正交，即同时满足以下两个方程：
求得满足上式的可调增益阵Kk，即可得到强跟踪滤波方法。
然而，对于非线性系统，利用正交原理,让残差序列时时正交很难精确满足。为减少计算量，通常采用近似算法。本实施例中提出一种带多重时变渐消因子的次优强跟踪滤波算法,采用时变的渐消因子对过去的数据进行渐消，通过实时调整状态预报误差的协方差阵，来减弱老数据对当前滤波值的影响。
预测误差协方差阵的计算过程如下：
Λ(k+1)＝diag{λ1(k+1),…,λn(k+1)}
式中，遗忘因子0.95≤ρ≤0.995。
为了进一步验证本发明的显著效果，本实施例进行了进一步的实验与分析。
目前，可供腹壁胎儿心电信号处理研究的开源数据库有三个：比利时鲁汶大学的Daily数据库、西班牙巴伦西亚大学的非入侵式胎儿心电数据库和波兰西里西亚大学的腹壁与直接胎儿心电数据库。前两个数据库包含若干条腹壁胎儿心电数据采集记录，是最为常用的实测数据，用于验证胎儿心电分离提取算法的效能。但这两个数据库都没有给出标准的胎儿心电波形或者胎儿瞬时心率。这导致基于实测数据的算法性能分析还处于肉眼观察的水平，算法性能的量化研究必须采用人工合成数据来实现。
最近，波兰西里西亚大学提供的数据库在Physionet上进行公开，该数据库首次给出了从胎儿头皮同步记录的胎儿心电信号，数据采集于5位孕妇，没有胸部电极，四个腹壁电极以肚脐为中心进行配置，头皮电极进行同步采集，采样率1KHz，分辨率16位，孕期38-41周，记录时长5分钟。头皮胎儿心电直接从胎儿体表采集，能够获得清晰完整的胎儿心电波形，如图2所示，通过定位R波可获得准确的胎儿瞬时心率，这为验证腹壁胎儿心率估计算法提供了一个很好的参考。为此，本实施例以该数据库中头皮胎儿心电数据为参考，利用腹壁心电数据记录就本文提出的胎儿心率估计算法进行量化评估。
本实施例选取第四号记录中胎儿心电最弱的第一通道进行算法测试。图2(a)、(b)分别给出了原始的腹壁胎儿心电信号和同步采集的头皮胎儿心电波形。可以看到，原始腹壁心电存在严重的肌电噪声与基线漂移，与同步采集的头皮胎儿心电信号相对应，可以看到腹壁采集的胎儿心电十分微弱，特别是第3至第6个胎儿心拍几乎完全被肌电噪声与母体心电所淹没，如图2(a)中圆圈所示。
首先从母体心电与胎儿心电R波的波形估计比较EKF递归估计算法与强跟踪滤波算法的性能。图3、图4分别显示了腹壁胎儿心电EKF波形估计与强跟踪波形估计的对比图。图3(a)中显示母体心电EKF波形估计结果图。从图3样本点6000附近可以看到，由于母体心电本身QRS波群逐搏变化以及肌电噪声的影响，导致母体心电与胎儿心电的波形跟踪出现较大误差，EKF估计的胎儿R波存在多个相似的脉冲波形，对R波定位存在很大干扰，严重影响瞬时胎心率估计(图中圆圈所示)。另一方面，强跟踪估计对母体心电波形的变化以及干扰噪声有很强的鲁棒性，可以准确估计出胎儿心电与母体心电波形，如图4圆圈中所示。同时，即使当胎儿R波完全被母体心电T波或QRS复波所覆盖，本发明提出的强跟踪滤波方法可以有效分离母体心电与胎儿心电R波，如图4方框所示。
进一步，本实施例从胎儿心率入手量化评估EKF递归估计算法与强跟踪滤波算法的性能。图5给出了强跟踪滤波算法下头皮胎儿心电R波定位与腹壁胎儿心电R波定位的对比图，可见二者都可以准确定位R波峰值。图6为强跟踪滤波算法下头皮胎儿瞬时心率与腹壁胎儿瞬时心率的对比图,从图6中可以看到，单通道获取的腹壁胎儿瞬时心率围绕头皮胎儿瞬时心率上下波动，二者最大差值不超过4拍/分(beat per minute,BPM)。它们的平均心率十分接近，仅相差0.0785BPM。
在该数据库中随机选取1分钟的单通道腹壁胎儿心电数据与头皮心电数据，进行强跟踪胎儿心率估计。去掉前后端的抖动部分，图7(a)为1分钟内连续100个心拍的腹壁胎儿瞬时心率估计情况，图7(b)为1分钟内连续100个心拍的头皮胎儿瞬时心率估计情况，图8为腹壁胎儿瞬时心率与头皮胎儿瞬时心率对比图。从图7(a)和图7(b)可以看到，总体上，估计的腹壁胎儿心率很好地跟踪了头皮胎儿心率，但存在一定的波动，特别是在第40心拍与第60心拍之间波动较大，误差都在5BPM以内。图9为腹壁胎心率与头皮胎心率的相关性测试，可见腹壁胎心率与头皮胎心率有很好地相关性。
由于腹壁胎儿心电瞬时心率存在一定的波动，为此将每10秒内的胎儿心率进行平均，对腹壁胎儿分段平均的胎心率与头皮胎儿心率再进行对比分析。选取5分钟单通道腹壁胎儿心电与头皮胎儿心电数据进行实验。图10给出了5分钟的实验结果，从图10可以看到分段平均之后，腹壁胎儿心率更为平滑，与头皮胎儿心率更为接近。图11给出了头皮胎儿心率与腹壁胎儿心率的误差图，图11中给出的误差为标准差。可见，两者的分段平均误差都在2BPM以内。
综合上述，从胎儿R波波形跟踪和胎儿心率估计两方面实验结果来看，本发明提出的胎儿心电鲁棒处理方法可以有效地从单通道母体腹壁心电中可靠地恢复胎儿R波波形，实现胎儿R波的准确定位，获得的胎儿瞬时心率及分段平均心率与同步采集的头皮胎儿心率十分接近，绝大部分误差低于2BPM，最大误差不超过5BPM，适用于家庭胎儿心率动态监护。
本申请的上述实施例中，通过提供一种用于家庭监护的单通道母体腹壁胎儿心率鲁棒估计方法，将母体心电作为胎儿心电的一部分进行整体建模，从波形形态角度建立腹壁胎儿心电信号状态空间模型，运用强跟踪滤波方法实现胎儿心电R波与母体心电波形的递推估计，并对估计的胎儿心电R波波形进行峰值定位，从而获取胎儿瞬时心率。
应当指出的是，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的普通技术人员在本发明的实质范围内所做出的变化、改性、添加或替换，也应属于本发明的保护范围。
Claims (4)
1.一种用于家庭监护的单通道母体腹壁胎儿心率鲁棒估计方法，其特征在于，包括以下步骤：
S1：信号预处理：对腹壁胎儿的心电信号进行预处理，以去除心电信号中的基线漂移和工频干扰；
S2：构建单通道腹壁胎儿心电信号的状态空间模型，通过系统参数初始化得到模型所需相关参数的初始值；
S3：强跟踪波形估计：根据胎儿心电信号的状态空间模型及提取的模型参数，运用强跟踪滤波方法，以实现胎儿心电R波与母体心电各个波形的递推估计；
S4：胎儿R波定位：采用递归Kalman滤波方法，实现R波动态定位，从而获取胎儿瞬时心率；
其中，步骤S2中构建单通道腹壁胎儿心电信号的状态空间模型的方法如下：
首先，将母体心电作为胎儿心电的一部分，构建一种六状态的胎儿心电信号动态模型：
式中，为k时刻胎儿心电的瞬时相位，为k时刻胎儿心电的瞬时心率，为k时刻母体心电的瞬时相位，为k时刻母体心电的瞬时心率，δ为采样时间间隔，k+1时刻胎儿心电R波用一个高斯函数进行建模，表示为Rk+1， 分别为胎儿心电R波对应高斯核函数的幅度、宽度和中心位置，k+1时刻母体心电P波和T波分别用两个高斯核函数进行建模，表示为Pk+1和Tk+1，k+1时刻母体心电QRS波群采用三个高斯核函数进行建模，表示为Ck+1，αi,k和bi,k，i∈{P1,P2,Q,R,S,T1,T2}分别为母体心电各个波形对应高斯核函数的幅度和宽度，分别为母体心电各个波形对应高斯核函数的中心位置，为高斯噪声；
由所述胎儿心电信号动态模型得到系统的观测模型为：
式中，为胎儿心电相位量测噪声，为母体心电相位量测噪声，Sk为胎儿心电信号测量样本，υS,k为量测噪声；
令Xk表示状态矢量，Yk表示观测向量，Uk表示系统噪声，Vk表示量测噪声，函数f和g分别描述系统方程和观测方程，则腹壁胎儿心电信号的状态空间模型为：
2.根据权利要求1所述的用于家庭监护的单通道母体腹壁胎儿心率鲁棒估计方法，其特征在于，步骤S2中系统参数初始化包括母体参数初始化和胎儿参数初始化，其中，母体参数初始化具体包括：对预处理之后的心电信号进行R波探测，获得母体瞬时心率的均值与方差，同时将心电信号从时域转换为相位域，从而获得母体心电的瞬时相位运用单通道奇异值分解方法获取母体平均心拍及平均心拍的高斯核参数；
胎儿参数初始化具体包括：从预处理之后的心电信号中去除提取的母体心电成分，获得胎儿心电初步估计，再进行R波探测，获得胎儿心率的均值与方差，同时将心电信号从时域转换为相位域，从而获得胎儿心电的瞬时相位运用单通道奇异值分解方法获取胎儿平均心拍及平均心拍的高斯核参数。
3.根据权利要求1所述的用于家庭监护的单通道母体腹壁胎儿心率鲁棒估计方法，其特征在于，步骤S3中强跟踪滤波方法具体为：
对于胎儿心电信号的状态空间模型的非线性状态方程，扩展Kalman滤波递归估计算法为：
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<msubsup>
<mover>
<mi>X</mi>
<mo>^</mo>
</mover>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
<mo>-</mo>
</msubsup>
<mo>=</mo>
<mi>f</mi>
<mrow>
<mo>(</mo>
<msubsup>
<mover>
<mi>X</mi>
<mo>^</mo>
</mover>
<mi>k</mi>
<mo>+</mo>
</msubsup>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>K</mi>
<mi>k</mi>
</msub>
<mo>=</mo>
<msubsup>
<mi>P</mi>
<mi>k</mi>
<mo>-</mo>
</msubsup>
<msubsup>
<mi>C</mi>
<mi>k</mi>
<mi>T</mi>
</msubsup>
<msup>
<mrow>
<mo>&lsqb;</mo>
<msub>
<mi>C</mi>
<mi>k</mi>
</msub>
<msubsup>
<mi>P</mi>
<mi>k</mi>
<mo>-</mo>
</msubsup>
<msubsup>
<mi>C</mi>
<mi>k</mi>
<mi>T</mi>
</msubsup>
<mo>+</mo>
<msub>
<mi>D</mi>
<mi>k</mi>
</msub>
<msub>
<mi>R</mi>
<mi>k</mi>
</msub>
<msubsup>
<mi>D</mi>
<mi>k</mi>
<mi>T</mi>
</msubsup>
<mo>&rsqb;</mo>
</mrow>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msubsup>
<mi>P</mi>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
<mo>-</mo>
</msubsup>
<mo>=</mo>
<msub>
<mi>A</mi>
<mi>k</mi>
</msub>
<msubsup>
<mi>P</mi>
<mi>k</mi>
<mo>+</mo>
</msubsup>
<msubsup>
<mi>A</mi>
<mi>k</mi>
<mi>T</mi>
</msubsup>
<mo>+</mo>
<msub>
<mi>B</mi>
<mi>k</mi>
</msub>
<msub>
<mi>Q</mi>
<mi>k</mi>
</msub>
<msubsup>
<mi>B</mi>
<mi>k</mi>
<mi>T</mi>
</msubsup>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>r</mi>
<mi>k</mi>
</msub>
<mo>=</mo>
<msub>
<mi>Y</mi>
<mi>k</mi>
</msub>
<mo>-</mo>
<mi>g</mi>
<mrow>
<mo>(</mo>
<msubsup>
<mover>
<mi>X</mi>
<mo>^</mo>
</mover>
<mi>k</mi>
<mo>-</mo>
</msubsup>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msubsup>
<mover>
<mi>X</mi>
<mo>^</mo>
</mover>
<mi>k</mi>
<mo>+</mo>
</msubsup>
<mo>=</mo>
<msubsup>
<mover>
<mi>X</mi>
<mo>^</mo>
</mover>
<mi>k</mi>
<mo>-</mo>
</msubsup>
<mo>+</mo>
<msub>
<mi>K</mi>
<mi>k</mi>
</msub>
<msub>
<mi>r</mi>
<mi>k</mi>
</msub>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msubsup>
<mi>P</mi>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
</mrow>
<mo>+</mo>
</msubsup>
<mo>=</mo>
<msubsup>
<mi>P</mi>
<mi>k</mi>
<mo>-</mo>
</msubsup>
<mo>-</mo>
<msub>
<mi>K</mi>
<mi>k</mi>
</msub>
<msub>
<mi>C</mi>
<mi>k</mi>
</msub>
<msubsup>
<mi>P</mi>
<mi>k</mi>
<mo>-</mo>
</msubsup>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
式中，为状态矢量的先验估计值，为状态矢量的后验估计值，Kk为可调增益阵，为先验估计的协方差矩阵，为后验估计的协方差矩阵，Rk为观测噪声方差，Qk为状态噪声方差，rk为更新后的信息，即新息，Ak、Bk、Ck、Dk为胎儿心电信号状态空间模型的Jacobin矩阵，
在扩展Kalman滤波递归估计算法上，强迫不同时刻的新息矢量rk正交，即同时满足以下两个方程：
<mfenced open = "" close = "">
<mtable>
<mtr>
<mtd>
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<mi>E</mi>
<mo>&lsqb;</mo>
<msub>
<mi>X</mi>
<mi>k</mi>
</msub>
<mo>-</mo>
<msub>
<mover>
<mi>X</mi>
<mo>^</mo>
</mover>
<mi>k</mi>
</msub>
<mo>&rsqb;</mo>
<msup>
<mrow>
<mo>&lsqb;</mo>
<msub>
<mi>X</mi>
<mi>k</mi>
</msub>
<mo>-</mo>
<msub>
<mover>
<mi>X</mi>
<mo>^</mo>
</mover>
<mi>k</mi>
</msub>
<mo>&rsqb;</mo>
</mrow>
<mi>T</mi>
</msup>
<mo>=</mo>
<mi>min</mi>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mi>E</mi>
<mo>&lsqb;</mo>
<msub>
<mi>r</mi>
<mrow>
<mi>k</mi>
<mo>+</mo>
<mi>j</mi>
</mrow>
</msub>
<msup>
<msub>
<mi>r</mi>
<mi>k</mi>
</msub>
<mi>T</mi>
</msup>
<mo>&rsqb;</mo>
<mo>=</mo>
<mn>0</mn>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
</mtd>
<mtd>
<mrow>
<mi>k</mi>
<mo>=</mo>
<mn>0</mn>
<mo>,</mo>
<mn>1</mn>
<mo>,</mo>
<mn>2</mn>
<mo>,</mo>
<mn>...</mn>
</mrow>
</mtd>
<mtd>
<mrow>
<mi>j</mi>
<mo>=</mo>
<mn>0</mn>
<mo>,</mo>
<mn>1</mn>
<mo>,</mo>
<mn>2</mn>
<mo>,</mo>
<mn>...</mn>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
求得满足上式的可调增益阵Kk，即可得到强跟踪滤波方法。
4.根据权利要求1所述的用于家庭监护的单通道母体腹壁胎儿心率鲁棒估计方法，其特征在于，步骤S2中系统参数初始化过程中选取10秒的心电数据。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510883842.0A CN105342594B (zh) | 2015-12-04 | 2015-12-04 | 用于家庭监护的单通道母体腹壁胎儿心率鲁棒估计方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510883842.0A CN105342594B (zh) | 2015-12-04 | 2015-12-04 | 用于家庭监护的单通道母体腹壁胎儿心率鲁棒估计方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN105342594A CN105342594A (zh) | 2016-02-24 |
| CN105342594B true CN105342594B (zh) | 2018-01-30 |
Family
ID=55318775
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201510883842.0A Expired - Fee Related CN105342594B (zh) | 2015-12-04 | 2015-12-04 | 用于家庭监护的单通道母体腹壁胎儿心率鲁棒估计方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN105342594B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107684423A (zh) * | 2016-08-05 | 2018-02-13 | 深圳先进技术研究院 | 一种胎儿心电分离方法及装置 |
| CN108836315A (zh) * | 2018-07-23 | 2018-11-20 | 智联时空(北京)科技有限公司 | 用于胎儿心电监测的智能腰带 |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101790346A (zh) * | 2007-07-24 | 2010-07-28 | 皇家飞利浦电子股份有限公司 | 监测胎儿心率的方法 |
| TW201114406A (en) * | 2009-10-21 | 2011-05-01 | Ind Tech Res Inst | Apparatus and method for maternal-fetal surveillance |
| CN102670188A (zh) * | 2012-04-25 | 2012-09-19 | 重庆大学 | 一种用于单导胎儿心率检测的鲁棒自适应估计方法 |
| CN104382589A (zh) * | 2014-12-09 | 2015-03-04 | 南京大学 | 基于部分按段重采样的胎儿心电图分离提取方法 |
| KR101504487B1 (ko) * | 2014-05-23 | 2015-03-23 | 광주과학기술원 | 태아의 심박동수를 실시간으로 측정하기 위한 시스템 |
| CN104462867A (zh) * | 2014-12-11 | 2015-03-25 | 中国人民解放军重庆通信学院 | T波动态建模与多通道融合估计方法 |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20130018273A1 (en) * | 2011-07-15 | 2013-01-17 | Tufts University | Systems and methods for analysis of fetal heart rate data |
-
2015
- 2015-12-04 CN CN201510883842.0A patent/CN105342594B/zh not_active Expired - Fee Related
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101790346A (zh) * | 2007-07-24 | 2010-07-28 | 皇家飞利浦电子股份有限公司 | 监测胎儿心率的方法 |
| TW201114406A (en) * | 2009-10-21 | 2011-05-01 | Ind Tech Res Inst | Apparatus and method for maternal-fetal surveillance |
| CN102670188A (zh) * | 2012-04-25 | 2012-09-19 | 重庆大学 | 一种用于单导胎儿心率检测的鲁棒自适应估计方法 |
| KR101504487B1 (ko) * | 2014-05-23 | 2015-03-23 | 광주과학기술원 | 태아의 심박동수를 실시간으로 측정하기 위한 시스템 |
| CN104382589A (zh) * | 2014-12-09 | 2015-03-04 | 南京大学 | 基于部分按段重采样的胎儿心电图分离提取方法 |
| CN104462867A (zh) * | 2014-12-11 | 2015-03-25 | 中国人民解放军重庆通信学院 | T波动态建模与多通道融合估计方法 |
Non-Patent Citations (1)
| Title |
|---|
| 徐雯，等.基于扩展卡尔曼和多项式网络的单导.《浙江省信号处理学会2015年年会会议论文集-信号处理与大数据》.2015,127-133. * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN105342594A (zh) | 2016-02-24 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US11476000B2 (en) | Methods and systems using mathematical analysis and machine learning to diagnose disease | |
| Andreotti et al. | Robust fetal ECG extraction and detection from abdominal leads | |
| Behar et al. | Combining and benchmarking methods of foetal ECG extraction without maternal or scalp electrode data | |
| US11234613B2 (en) | Respiration estimation method and apparatus | |
| Liu et al. | A multi-step method with signal quality assessment and fine-tuning procedure to locate maternal and fetal QRS complexes from abdominal ECG recordings | |
| US9451901B2 (en) | Method and device for evaluation of myocardial damages based on the current density variations | |
| Tavakolian et al. | Estimation of hemodynamic parameters from seismocardiogram | |
| Vaneghi et al. | A comparative approach to ECG feature extraction methods | |
| Jamshidian-Tehrani et al. | Fetal ECG extraction from time-varying and low-rank noninvasive maternal abdominal recordings | |
| Su et al. | Recovery of the fetal electrocardiogram for morphological analysis from two trans-abdominal channels via optimal shrinkage | |
| CN105342594B (zh) | 用于家庭监护的单通道母体腹壁胎儿心率鲁棒估计方法 | |
| Keenan et al. | Detection of fetal arrhythmias in non-invasive fetal ECG recordings using data-driven entropy profiling | |
| CN116172539A (zh) | 基于机器学习的生命体征检测方法、系统、设备及介质 | |
| Wang et al. | Contactless radar heart rate variability monitoring via deep spatio-temporal modeling | |
| Sæderup et al. | Comparison of cardiotocography and fetal heart rate estimators based on non-invasive fetal ECG | |
| CN103637796B (zh) | 基于广义特征值最大化的胎儿心电信号自适应盲提取方法 | |
| Li et al. | Robust adaptive fetal heart rate estimation for single-channel abdominal ECG recording | |
| Marzbanrad et al. | Classification of Doppler ultrasound signal quality for the application of fetal valve motion identification | |
| CN103876731A (zh) | 一种胎儿心电信号提取装置及方法 | |
| Katebi et al. | Unsupervised hidden semi-Markov model for automatic beat onset detection in 1D Doppler ultrasound | |
| Kang et al. | Non-Contact Realtime Vital Signs Monitoring System Based on Millimeter Wave FMCW Radar | |
| Wang et al. | Fetal electrocardiogram extraction algorithm in noise: using BSE | |
| Zhidong | Noninvasive diagnosis of coronary artery disease based on instantaneous frequency of diastolic murmurs and SVM | |
| Zamrath et al. | Robust and computationally efficient approach for Heart Rate monitoring using photoplethysmographic signals during intensive physical excercise | |
| Wan et al. | A method for extracting fECG based on ICA algorithm |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| C10 | Entry into substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20180130 Termination date: 20201204 |