CN103462642A - 基于时频分析的多普勒胎心音的瞬时心率检测方法及装置 - Google Patents

Abstract

本发明涉及医用监护技术领域，公开了一种基于时频分析的多普勒胎心音的瞬时心率检测方法及装置，检测方法包括：通过带通滤波器对多普勒胎心音进行预处理；对预处理后的超声多普勒胎心音进行时频分析，为了计算简单快捷，采用STFT求取超声多普勒胎心音的时频图；互相关法求胎儿心音的瞬时心率法：互相关算法中首先根据多普勒胎心音的时频图选取特征频段，再根据心音信号的先验知识选取特征模板，接着求特征频段和特征模板的互相关函数，绘制互相关曲线，最后通过计算互相关曲线的波峰值的间隔计算超声多普勒的胎儿心音信号的瞬时心率。本发明对于临床采集到的超声多普勒胎心音信号的瞬时心率的计算，方法简单，运算速度快，准确度高。

Description

基于时频分析的多普勒胎心音的瞬时心率检测方法及装置

技术领域

本发明涉及医用监护技术领域，更具体地，涉及一种基于时频分析的多普勒胎心音的瞬时心率检测方法及装置。

背景技术

胎心监护是通过监测胎心率的变化来评价宫内安危的一种常用的胎儿监护方法。在围产阶段对母体内的胎儿进行监护，可以在了解胎儿发育状况的同时，大大地减少因缺氧、缺血而出现的窘迫，降低胎儿出生缺陷，甚至胎儿死亡的情况的发生。出生缺陷目前已成为影响我国人口质量的一个严重问题。因此，密切关注胎儿心率的变化情况对提高生育质量有着重要的意义。早至19世纪初就有产科医师通过胎心听诊估计胎儿宫内状况，随着超声多普勒技术的发展，产时电子胎心监护（electronicfetal monitoring，EFM）已成为目前最常用的胎儿监护方法。超声多普勒测量方法又是目前最常用的胎儿心率测量方法。

但是，由于从超声波换能器上检测的超声多普勒声音混杂有幅度大，分布广的噪声干扰，以及胎儿身体在母体内运动时，声音信号强弱变化明显。在时域和频域上，这些干扰信号混叠在一起对胎儿心音信号的瞬时心率计算造成很大的影响。因此，研究如何准确、有效地测量母体内胎儿心音的瞬时心率具有重要的理论价值和临床意义。

对胎心的监护和胎儿瞬时心率的研究在很早就已经开始了，并且采用的处理方法有很多种，总结起来，可分为以下几种：

（1）基于匹配滤波法计算胎儿心率：该方法的基本思想是用以往多次得到的母体心电作为模板，抵消腹部采样信号中的母体心电成分，将胎儿心电分离出来。由于腹部信号与模板相减要求有较高的精度，则需在模板计算、相位、幅值修正上采取多种措施来保证减去母体心电的准确性，这是基于心电模式的方法。

（2）采用自相关的方法计算胎儿心率；总所周知，相关方法是从加性噪声中提取已知波形，特别是对于确定性的周期信号是一种好方法。用自相关方法提取胎儿心率信号，效果不够好，主要是因为胎儿心率是重复性信号，但不是确定性的周期信号，并且胎儿心音信号的波形也不是固定不变的，周期与波形均带有随机性变化，因而造成自相关函数的波形不易检测；这是基于心音模式的方法。

正常心脏有四个心音：第一(S1)、第二(S2)、第三(S3)、第四(S4)心音，但多数情况下只能听到第一和第二心音。第一心音的出现标志着心室收缩期的开始；第二心音的出现标志心室缩张期的开始。从出现第一心音到出现第二心音的期间定义为心室收缩期；从出现第二心音到下一心动周期中出现第一心音的期间定义为心室舒张期。在一个心动周期中，心音的主要成分为：第一心音、收缩期、第二心音、舒张期，通过它们可以完整的表述出心音的时间特性。对于一般人来说，正常情况下心脏的收缩期短于舒张期。胎儿的平均心率为120-160次/分，心动周期约为0.5(秒)，其中收缩期大约占0.2(秒)，舒张期大约占0.3(秒)。也就是说，一般人的心音信号中，S1音到S2音在时域上的间隔大约为0.2(秒)。

由于胎儿心音信号属于非平稳信号，常规傅里叶变换方法不能刻画任一时刻的频率成分，致使无法对其进行全面的分析。时频分析是非平稳信号分析的一个有力工具，该方法能够将一维时域信号变换到二维的时频平面，提供了时间域和频率域的联合分布信息，清楚的描述了信号频率随时间变换的关系。

发明内容

本发明的主要目的在于克服现有胎儿心率检测方法的不足，提供一种基于时频分析的多普勒胎心音的瞬时心率检测方法，该检测方法将胎儿心音的时间域和频率域的分布信息联合起来，再利用了心音信号的先验信息——心音信号中S1音和S2音在观测信号时域上的间隔，从而检测胎心音的瞬时心率。

为解决上述技术问题，本发明的技术方案如下：

一种基于时频分析的多普勒胎心音的瞬时心率检测方法，包括如下步骤：

S1.信号预处理：对采集到的多普勒胎心音进行带通滤波，带通滤波器的通带频率为f_L到f_H；

S2.时频分析：对于经过步骤S1预处理的多普勒胎心音进行时频分析，得到时频图；

S3.选取特征频段和模板：从时频图中选取信号中的特征频段，特征频段为f_cL到f_cH；从时频图选取包含S1音和S2音特征的时频块，所述时频块的时间间隔为0.2(秒)＜t₀＜0.5(秒)；

S4.求互相关函数：对特征频段和模板求互相关函数，根据互相关函数的结果绘出相关曲线；

S5.用阈值检测法求互相关曲线的峰值；

S6.对峰值做差分，求瞬时心率值；根据瞬时心率值绘制出瞬时心率图。

更进一步的，所述步骤S1中f_L为50Hz，f_H为250Hz，即采用带通频带为50-250Hz的带通滤波器。

更进一步的，所述步骤S2采用短时傅里叶变换STFT进行时频分析，短时傅里叶变换STFT定义如下：

$s(w,t)=\frac{1}{2\mathrm{\π}}{\∫}_{-\∞}^{+\∞}{e}^{-\mathrm{iw}}x\left(\mathrm{\τ}\right)h(\mathrm{\τ}-t)\mathrm{d\τ}---\left(1\right)$

其中h(t)是窗函数，x(τ)是信号，τ为信号自变量，t是时间变量，w为频率自变量，即沿时间轴移动分析窗，得到二维的时频图由s(w,t)表示。

更进一步的，所述步骤S3中特征频段取200-400Hz，即f_cL为200Hz，f_cH为400Hz。

特征模版可以通过将一维胎心音信号投影到二维时频平面后，根据特征频段（fcl（200Hz）到fch（400Hz））上胎心音信号的特征，提取出特征明显，时间长度为0.2到0.5秒的时频块作为特征模板。特征模板也可以通过提取一段完整，正常的胎心音，投影到二维时频平面生成，但生成的模板时间长度应在0.2到0.5秒内。

更进一步的，所述步骤S4中采用二维互相关函数绘制互相关曲线，公式如下：

$C(i,j)=\underset{m=0}{\overset{\mathrm{Ma}-1}{\mathrm{\Σ}}}\underset{n=0}{\overset{\mathrm{Na}-1}{\mathrm{\Σ}}}A(m,n)\·\mathrm{conj}\left(B(m+i,n+j)\right)---\left(2\right)$

其中A为Ma行Na列的矩阵，B为Mb行Nb列的矩阵，conj(B)表示B的共轭复数，0≤i<Ma+Mb-1，并且0≤j<Na+Nb-1，C(i,j)表示互相关曲线。

更进一步的，所述步骤S5中用阈值检测法求互相关曲线的峰值，

阈值为：threshold=param×max{R(n)}，param的取值范围为0-1，R(n)表示互相关曲线。

更进一步的，所述步骤S6求瞬时心率的计算公式为：

或

本发明的又一目的是提出一种应用基于时频分析的多普勒胎心音的瞬时心率检测方法的装置，采用该检测装置能够准确的从胎儿心音信号中获取瞬时心率，该检测装置包括：

信号预处理模块，用于对采集到的多普勒胎心音进行带通滤波，其中带通滤波器的通带频率为f_L到f_H；

时频分析模块，用于对预处理的多普勒胎心音进行时频分析，得到时频图；

特征频段和模板的选取模块，用于从时频图中选取信号中的特征频段和特征频段和模板，其中特征频段为f_cL到f_cH，时频块的时间间隔为0.2(秒)＜t₀＜0.5(秒)；

互相关模块，用于对特征频段和模板求互相关函数，并根据互相关函数的结果绘出相关曲线；

峰值提取模块，用于采用阈值检测法提取互相关曲线的峰值；

瞬时心率图绘制模块，用于对峰值做差分，求瞬时心率值；并根据瞬时心率值绘制出瞬时心率图。

更进一步的，所述信号预处理模块的带通滤波器的带通频带为50-250Hz，即f_L为50Hz，f_H为250Hz。

更进一步的，所述时频分析模块采用采用短时傅里叶变换STFT进行时频分析，短时傅里叶变换STFT定义如下：

$s(w,t)=\frac{1}{2\mathrm{\&pi;}}{\&Integral;}_{-\&infin;}^{+\&infin;}{e}^{-\mathrm{iw}}x\left(\mathrm{\&tau;}\right)h(\mathrm{\&tau;}-t)\mathrm{d\&tau;}---\left(1\right)$

其中h(t)是窗函数，x(τ)是信号，τ为信号自变量，t是时间变量，w为频率自变量，即沿时间轴移动分析窗，得到二维的时频图由s(w,t)表示。

与现有技术相比，本发明技术方案的有益效果是：本发明的瞬时心率检测方法是利用时频分析将一维非平稳的胎儿心音信号转换到能清楚的描述信号频率随时间变换的二维时频平面，再在二维时频平面上利用S1音和S2音的先验信息提取特征模板，计算特征模板与特征频段的归一化互相关曲线从而计算出瞬时心率，本发明的检测方法准确率比常规方法高。

附图说明

图1为本发明瞬时心率检测方法的流程图。

图2为临床采集的超声多普勒胎儿心音信号示意图。

图3为利用STFT对胎儿心音信号进行时频变化后的二维时频平面图。

图4为归一化互相关曲线示意图。

图5为采用本发明的检测方法检测到的胎儿心音信号的瞬时心率示意图。

具体实施方式

下面结合附图和实施例对本发明的技术方案做进一步的说明。

如图1所示，本发明基于时频分析的多普勒胎心音的瞬时心率检测方法，包括如下步骤：

S1.信号预处理，主要对采集到的多普勒胎心音进行带通滤波，带通滤波器的通带频率为f_L到f_H，其中采集到的多普勒胎心音信号如图2所示；在本实施例中带通滤波器的带通频带为50-250Hz，即f_L为50Hz，f_H为250Hz；

S2.时频分析，对于经过步骤S1中预处理的多普勒胎心音进行时频分析，本发明中采用短时傅里叶变换方法进行时频分析，得到的二维时频平面图如图3所示；短时傅里叶变换STFT是一种时频分析方法，其定义如下：

$s(w,t)=\frac{1}{2\mathrm{\&pi;}}{\&Integral;}_{-\&infin;}^{+\&infin;}{e}^{-\mathrm{iw}}x\left(\mathrm{\&tau;}\right)h(\mathrm{\&tau;}-t)\mathrm{d\&tau;}---\left(1\right)$

其中h(t)是窗函数，x(τ)是信号。τ为信号自变量，t是时间变量，w为频率自变量。沿时间轴移动分析窗，可以得到二维的时频平面，由s(w,t)表示。

S3.选取特征频段，针对图3的时频图，选取信号中的200-400Hz的特征频段，即特征频段f_cL为200Hz，f_cH为400Hz

选取模板，针对图3的时频图，选取时频图中包含S1音和S2音特征的时频块，注意，时频块的时间间隔t₀一定要大于0.2S小于0.5S；在本实施例中，时频块的时间间隔t₀=0.4S。

S4.求互相关函数，对步骤S3中的特征频段和得到的模板求互相关函数，根据互相关函数的结果绘出互相关曲线，绘制的互相关曲线如图4所示；

其中二维互相关函数的计算公式如下：

$C(i,j)=\underset{m=0}{\overset{\mathrm{Ma}-1}{\mathrm{\&Sigma;}}}\underset{n=0}{\overset{\mathrm{Na}-1}{\mathrm{\&Sigma;}}}A(m,n)\&CenterDot;\mathrm{conj}\left(B(m+i,n+j)\right)---\left(2\right)$

其中A为Ma行Na列的矩阵，B为Mb行Nb列的矩阵，conj(B)表示B的共轭复数，0≤i<Ma+Mb-1，并且0≤j<Na+Nb-1，C(i,j)表示互相关曲线。

S5.用阈值检测法求互相关曲线的峰值；在本实施例中，阈值的求取方式为：

threshold=param×max{R(n)}，这里param取值范围为0-1，R(n)表示互相关曲线，具体实施例的也可以根据需要取0.9以上的值。

S6.对峰值做差分，求瞬时心率值；

瞬时心率的计算公式为：

或

在本实施例中采用公式（3）求取瞬时心率，根据瞬时心率值绘制出瞬时心率图，瞬时心率图如图5所示。

显然，本发明的上述实施例仅仅是为清楚地说明本发明所作的举例，而并非是对本发明的实施方式的限定。对于所属领域的普通技术人员来说，在上述说明的基础上还可以做出其它不同形式的变化或变动。这里无需也无法对所有的实施方式予以穷举。凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明权利要求的保护范围之内。

Claims (10)

1.一种基于时频分析的多普勒胎心音的瞬时心率检测方法，其特征在于，包括如下步骤：

S1.信号预处理：对采集到的多普勒胎心音进行带通滤波，带通滤波器的通带频率为f_L到f_H；

S2.时频分析：对于经过步骤S1预处理的多普勒胎心音进行时频分析，得到时频图；

S3.选取特征频段和模板：从时频图中选取信号中的特征频段，特征频段为f_cL到f_cH；从时频图选取包含S1音和S2音特征的时频块，所述时频块的时间间隔为0.2(秒)＜t₀＜0.5(秒)；

S4.求互相关函数：对特征频段和模板求互相关函数，根据互相关函数的结果绘出相关曲线；

S5.用阈值检测法求互相关曲线的峰值；

S6.对峰值做差分，求瞬时心率值；根据瞬时心率值绘制出瞬时心率图。

2.根据权利要求1所述的基于时频分析的多普勒胎心音的瞬时心率检测方法，其特征在于，所述步骤S1中f_L为50Hz，f_H为250Hz，即采用带通频带为50-250Hz的带通滤波器。

3.根据权利要求1所述的基于时频分析的多普勒胎心音的瞬时心率检测方法，其特征在于，所述步骤S2采用短时傅里叶变换STFT进行时频分析，短时傅里叶变换STFT定义如下：

$s(w,t)=\frac{1}{2\mathrm{\&pi;}}{\&Integral;}_{-\&infin;}^{+\&infin;}{e}^{-\mathrm{iw}}x\left(\mathrm{\&tau;}\right)h(\mathrm{\&tau;}-t)\mathrm{d\&tau;}---\left(1\right)$

其中h(t)是窗函数，x(τ)是信号，τ为信号自变量，t是时间变量，w为频率自变量，即沿时间轴移动分析窗，得到二维的时频图由s(w,t)表示。

4.根据权利要求1所述的基于时频分析的多普勒胎心音的瞬时心率检测方法，其特征在于，所述步骤S3中特征频段取200-400Hz，即f_cL为200Hz，f_cH为400Hz。

5.根据权利要求1所述的基于时频分析的多普勒胎心音的瞬时心率检测方法，其特征在于，所述步骤S4中采用二维互相关函数绘制相关曲线，公式如下：

$C(i,j)=\underset{m=0}{\overset{\mathrm{Ma}-1}{\mathrm{\&Sigma;}}}\underset{n=0}{\overset{\mathrm{Na}-1}{\mathrm{\&Sigma;}}}A(m,n)\&CenterDot;\mathrm{conj}\left(B(m+i,n+j)\right)---\left(2\right)$

其中A为Ma行Na列的矩阵，B为Mb行Nb列的矩阵，conj(B)表示B的共轭复数，0≤i<Ma+Mb-1，并且0≤j<Na+Nb-1，C(i,j)表示互相关曲线。

6.根据权利要求1所述的基于时频分析的多普勒胎心音的瞬时心率检测方法，其特征在于，所述步骤S5中用阈值检测法求互相关曲线的峰值，

阈值为：threshold=param×max{R(n)}，param取0.9，或取与0.9相近的数值，R(n)表示互相关曲线。

7.根据权利要求1所述的基于时频分析的多普勒胎心音的瞬时心率检测方法，其特征在于，所述步骤S6求瞬时心率的计算公式为：

或

8.一种应用权利要求1至7任一项所述的基于时频分析的多普勒胎心音的瞬时心率检测方法的装置，其特征在于，包括：

信号预处理模块，用于对采集到的多普勒胎心音进行带通滤波，其中带通滤波器的通带频率为f_L到f_H；

时频分析模块，用于对预处理的多普勒胎心音进行时频分析，得到时频图；

特征频段和模板的选取模块，用于从时频图中选取信号中的特征频段和特征频段和模板，其中特征频段为f_cL到f_cH，时频块的时间间隔为0.2(秒)＜t₀＜0.5(秒)；

互相关模块，用于对特征频段和模板求互相关函数，并根据互相关函数的结果绘出相关曲线；

峰值提取模块，用于采用阈值检测法提取互相关曲线的峰值；

瞬时心率图绘制模块，用于对峰值做差分，求瞬时心率值；并根据瞬时心率值绘制出瞬时心率图。

9.根据权利要求8所述的基于时频分析的多普勒胎心音的瞬时心率检测装置，其特征在于，所述信号预处理模块的带通滤波器的带通频带为50-250Hz，即f_L为50Hz，f_H为250Hz。

10.根据权利要求8或9所述的基于时频分析的多普勒胎心音的瞬时心率检测装置，其特征在于，所述时频分析模块采用采用短时傅里叶变换STFT进行时频分析，短时傅里叶变换STFT定义如下：

$s(w,t)=\frac{1}{2\mathrm{\&pi;}}{\&Integral;}_{-\&infin;}^{+\&infin;}{e}^{-\mathrm{iw}}x\left(\mathrm{\&tau;}\right)h(\mathrm{\&tau;}-t)\mathrm{d\&tau;}---\left(1\right)$

其中h(t)是窗函数，x(τ)是信号，τ为信号自变量，t是时间变量，w为频率自变量，即沿时间轴移动分析窗，得到二维的时频图由函数s(w,t)表示。

Images