CN101843489A - 一种呼吸信号处理方法 - Google Patents

Abstract

本发明公开了一种呼吸信号处理方法，包括以下步骤：A、通过特定的呼吸电路获取呼吸生理信号，AD采集得到呼吸数字信号，并对该信号进行数据预处理；B、对预处理后的呼吸信号进行频谱转换得到其频域分布；C、根据呼吸信号的生理参数特征进行是否窒息判断；D、引入心脏活动生理参数，分析呼吸信号的频谱分布，进行是否受到心脏活动干扰的判断；E、对呼吸波形频域能量包络分析，寻找正确的谱峰，将该谱峰对应的频率点转换为呼吸率，并综合历史呼吸率的值得到当前的呼吸率值。本发明的处理方法可最大限度的提高呼吸率计算结果准确性和稳定性大幅提高。

Description

一种呼吸信号处理方法

技术领域

本发明涉及一种呼吸信号处理方法，尤其涉及一种利用频谱法对呼吸时域波形进行转换，对波形频谱能量分布进行分析以达到去除干扰，提高呼吸检测的稳定性和准确性的呼吸信号处理方法。

背景技术

现有呼吸测量装置常使用一种基于阻抗法的测量方法来获得呼吸波信号。人体呼吸运动时，胸壁肌肉交变张驰，胸廓交替变形，胸腔的电阻抗也随之交替变化，通过检测电阻抗的微小变化，就可以得到呼吸的变化。实际的呼吸测量过程一般是借助体表心电信号检测中贴在体表特定位置的电极将高频载波信号施加到人体胸腔，通过它可以把胸腔的呼吸变化引起的微小阻抗变化调制到高频载波信号上，然后经过一系列电路对这个载波信号进行放大、检波、解调等，就可以得到模拟的呼吸信号，再经过A/D转换得到数字的呼吸信号，最后运用呼吸算法计算出呼吸检测参数呼吸率及窒息报警信息。

平静呼吸时，新生儿的呼吸率是30～70BPM(Beats Per Minute)，成人的是12～30BPM，但是如果考虑异常情况，一般要求呼吸检测范围是8～120BPM，个别可高达150BPM。于是可以得到该呼吸率检测范围所对应的呼吸波的频率为0.125～2.5Hz。

目前市场上呼吸检测的方法主要采用波形法，该方法通过一段时间内波形的平均值(即基线值)，来判定当前呼吸波处于上升或下降趋势，用极值的方法求得波形的波峰、波谷。根据一定的阈值条件来判定有效的波峰或波谷，再根据有效波峰或波谷的周期计算波形周期，从而得到呼吸率；并且根据一段时间内波形的幅度平均值大小来进行呼吸窒息判断。

虽然波形法的计算过程具有比较直观、运算量小的优点，但在实际临床过程中发现：当病人躁动导致波形紊乱时，经常有效波形周期找不准，导致计算呼吸率错误；阻抗法检测获得呼吸波形的原理，由于必须共用心电电极片，会受到心电活动的影响(即心动干扰Cardiovascular Artifact)，获取的呼吸波形或多或少会受到心电波形的干扰，尤其当病人熟睡时安静呼吸、心电电极片的位置不佳或短暂的呼吸阻塞时，呼吸波形在时域可能夹杂着很强的心动干扰甚至完全被心电干扰淹没，此时波形法无法正确区分心动干扰波形和呼吸波形；当波形出现基线漂移时，由于上述方法计算的基线值无法很快更新，会导致波形漏检致使呼吸率值偏低。总之，利用波形法进行呼吸检测的抗扰能力不强，无法准确地判定心动干扰。

发明内容

本发明的目的是为了克服波形法这种现有技术的不足，提供一种准确性高、稳定性好的呼吸信号处理方法。

为解决上述的技术问题，本发明的构思为：常规的呼吸检测方法都是对呼吸的时域波形进行分析，本发明从另外一种角度-频域对呼吸信号进行分析。波形的频域分布往往有许多时域观察不到的特征，尤其像呼吸波形具有周期性，其频域会在相应的频率处出现能量强度很大的谱峰。同样在时域无法辨识的呼吸波形和干扰波形，在频域也会在不同的频率点出现谱峰。从而把时域对呼吸波形和干扰波形的分离转换为频域的不同谱峰的分离，后者更易实现且准确。

为了达到上述目的，本发明采用以下技术方案：

A、通过特定的呼吸电路获取呼吸生理信号，AD采集得到呼吸数字信号，并对该信号进行带通滤波的数据预处理；

B、对预处理后的呼吸信号进行频谱转换得到其频域分布；

C、根据呼吸信号的生理参数特征进行是否窒息判断。

D、引入心脏活动生理参数，分析呼吸信号的频谱分布，进行是否受到心脏活动干扰(即心动干扰)的判断；

E、对呼吸波形频域能量包络进行分析，寻找正确的谱峰，将该谱峰对应的频率点转换为呼吸率，并综合历史呼吸率的值得到当前的呼吸率值。

上述方案中，步骤A所述的带通滤波器包括双向IIR滤波器或无乘法椭圆IIR滤波器。

上述方案中，步骤B所述的将呼吸信号进行频域转换的方法包括傅立叶变换、小波变换或希尔伯特变换。

上述方案中，步骤B所采用的常用频域转换方法之一傅立叶变换包括快速傅立叶变换(FFT)或线性调频Z变换(CZT)。

上述方案中，获得呼吸信号的频域分布后，如步骤C首先设置呼吸频谱能量阈值以及时域的呼吸波形幅度限值条件，呼吸信号在一段时间内频谱能量最大值小于设定的频域能量阈值线且时域幅度小于时域阈值线则判定呼吸信号窒息。该方法综合时域和频域的条件来判定呼吸波形是否窒息，以提高呼吸窒息检测的准确性。

上述方案中，步骤D引入的心脏活动生理参数包括心率值或脉率值。

对非窒息波形进行是否受心动干扰判断，方法为，求得呼吸频谱能量的最大值和次大值，若其对应的呼吸率的值在当前的心率值/或脉率值的设定范围内，则判定呼吸波形受到不同程度的心动干扰，计算呼吸率值时剔除该频率点，从而提高呼吸检测抗心动干扰的能力。

剔除心动干扰后，实施步骤E，该步骤还包括：

a.求呼吸频谱能量的极大值点确定若干个谱峰。

b.对步骤a求得的谱峰进行分析，根据各个谱峰的能量大小比例、谱峰对应的呼吸率值和历史呼吸率值的大小比例来确定一个正确的谱峰，并将其对应的频率点转换为呼吸率。

c.把步骤b计算所得的呼吸率值和历史呼吸率值做加权平均，获得当前的呼吸率的值。

采用上述的技术方案，可以比较直观有效地去除心动干扰的影响，提高窒息判断的可靠性，频域更容易分离干扰信号和正常的呼吸信号，再加入利用历史呼吸率对谱峰的判断，这些均可提高呼吸检测的抗干扰性，提高呼吸检测的稳定性和准确性。

附图说明

图1为本发明的流程图；

图2是正常呼吸的呼吸数据；

图3是本发明处理正常呼吸数据得到的呼吸率的结果图；

图4是呼吸中含有较大的心动干扰的呼吸数据；

图5是本发明处理含有较大的心动干扰的呼吸数据得到的呼吸率的结果图；

图6是病人躁动时获得的呼吸数据；

图7是本发明处理病人躁动时的呼吸数据得到的呼吸率的结果图；

图8是呼吸中含有较大干扰的呼吸数据；

图9是本发明处理含有较大干扰呼吸数据得到的呼吸率的结果图。

具体实施方式

下面结合附图和具体实施例对本发明进行进一步详细的说明。

用于呼吸信号测量的监护设备主要包括上位机和下位机，下位机主要是通过硬件电路来获取呼吸数据，上位机主要是接受来自下位机的数据，显示呼吸波形、呼吸率以及窒息报警信息。本发明主要在下位机中实现，若下位机单片机计算能力不足，本发明内容可移植至上位机中实现。该方法的计算流程如图1所示，其主要步骤包括：

A、对呼吸数字信号进行带通滤波数据预处理；

B、对预处理后的呼吸信号进行频谱转换得到其频域分布；

C、设定呼吸信号频域能量阈值线和呼吸信号时域波形幅度阈值线，呼吸信号在一段时间内频谱能量最大值小于设定的频域能量阈值线且时域波形幅度小于时域阈值线则判定呼吸信号窒息。

D、利用呼吸频谱能量的最大值和次大值对应的呼吸率值和当前的心率值(但不限于)限定条件来判定呼吸波形是否存在心动干扰。

E、对呼吸波形频域能量包络分析，寻找正确的谱峰，将该谱峰对应的频率点转换为呼吸率，把该呼吸率和历史呼吸率的值做平均加权获取当前的呼吸率的值。

本发明实例中步骤A所用的带通滤波器为无乘法IIR滤波器(但不限于)。

本发明实例中步骤B涉及的频谱转换的方法使用线性调频Z变换(Chirp Z-Transform，CZT)(但不限于)。其理论推导如下：

已知x(n)(0≤n≤N-1)是有限长序列，其Z变换为：

$X\left(z\right)=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}x\left(n\right)z^{-n}---\left(1\right)$

为了适应z可以沿z平面更一般的路径取值，所以沿z平面上的一段螺线作等分角的抽样，z的这些抽样点z_k为：

z_k＝AW^-k    k＝0，1，…，M-1    (2)

其中A和W分别为：

$\left\{\begin{array}{c}A=A_0{e}^{j{\mathrm{\θ}}_0}\\ W=W_0{e}^{-j{\mathrm{\φ}}_0}\end{array}\right.---\left(3\right)$

于是可得：

$z_k=A_0{e}^{j{\mathrm{\θ}}_0}*W_0^{-k}{e}^{j{\mathrm{\φ}}_{0}k}$ k＝0，1，…，M-1    (4)

其中：

A₀：表示抽样点的矢量半径长度。

θ₀：表示起始抽样点z₀的相角。

φ₀：表示两相邻抽样点之间的角度差。

W₀：表示螺线的伸展率。

将(4)式中的z_k带入(1)式可得：

k＝0，1，.…，M-1    (5)

为了方便计算，对CZT的计算公式经过一系列的变换，它可以改写为线性卷积的形式：

$X\left(z_k\right)=W^{\frac{k^2}{2}}\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}g\left(n\right)h(k-n)=W^{\frac{k^2}{2}}[g\left(k\right)*h\left(k\right)]$ k＝0，1，...，M-1    (6)

其中 $\left\{\begin{array}{c}g\left(n\right)=x\left(n\right)A^{-n}{W}^{\frac{n^2}{2}}\\ h\left(n\right)=W^{-\frac{n^2}{2}}\end{array}\right.$ n＝0，1，...，N-1    (7)

当M＝N，

W₀＝1，

这一特殊情况时，各z_k就均匀等间隔地分布在单位圆上，即求序列的傅立叶变换。通过设置θ₀，φ₀和采样点数就可以只对一定频率范围内的信号进行频谱分析。

由以上理论分析，根据呼吸波检测的频率范围集中在0.125～2.5Hz的特点，本发明实例确定CZT变换的频率范围为0～3Hz。

为了满足呼吸检测±1BPM分辨率的要求，可以得到CZT的分辨率计算公式为：其中f₁＝0Hz，f₂＝3Hz，得到CZT的最小点数为256。本发明实例用于呼吸的频域转换的采样率f_s＝25Hz，为了保证用于一次频谱变换的数据至少包含2个呼吸周期，对于如6BPM的低频信号，2个波形周期为20s的数据，本发明实例确定CZT的点数N＝512。

对某单一频率的信号进行CZT变换后，对其进行频谱分析，在该频率处会出现一个频谱的最大峰值；对多个频率的信号叠加后所形成的信号进行CZT变换，然后对其进行频谱分析，在各个频率处的频谱能量值都会出现一个局部极大值，在所有局部极大值中，最大的局部极大值所对应的频率表示该频率的信号在时域内叠加时的信号最强。

基于上面的理论基础，设定呼吸信号频域能量阈值线和呼吸信号时域波形幅度阈值线，呼吸信号在一段时间内频谱能量最大值小于设定的频域能量阈值线且时域波形幅度小于时域阈值线则判定呼吸信号窒息。

其后对有呼吸波形的情况进行分析，步骤如下：

首先求得频谱能量最大值和此大值所对应的呼吸率值，若频谱能量最大值所对应的呼吸率的值在心率的值±r的范围内，说明呼吸中有心动干扰很大，甚至淹没了正常的呼吸波形，此时建议进行心动干扰的报警，提示可能需要重新调整电极片的位置以获取更好的呼吸波形。如果频谱能量的次大值对应的呼吸率值在心率值的±r的范围内，而频谱能量的最大值对应的呼吸率值并不在该心率值的±r范围内，说明呼吸波形受心动干扰影响的强度并没有达到严重的程度，剔除该频率点即可，不进行报警。

其次对呼吸波形频谱能量包络中的各个局部极大值查找，确定几个可能的谱峰。

对上述求得的谱峰进行分析，根据谱峰的能量大小比例、谱峰对应的呼吸率值和历史呼吸率的值的比例关系进行剔除，以确定一个正确的谱峰。原则上，真实的呼吸波形的频域分布应该在频谱能量的最大值，次大值，第三大值上。

把求得的正确谱峰所对应的呼吸率值和历史呼吸率值做加权平均，获得当前的呼吸率值，以进一步提高呼吸率计算的稳定性。

为了更形象地阐述频谱法计算呼吸率，下面分别以实际临床中获取的正常呼吸波形、含有心动干扰的呼吸波形，含有躁动干扰的呼吸波形和呼吸中含有很大的干扰为例，分别说明频谱法在处理不同类型呼吸数据时呼吸率计算的准确性。

对于图2所示的正常呼吸数据，呼吸波形规则、没有干扰，该周期波形在频域中有很集中的能量分布，如图3中所示。图3星号处的频谱能量是所有频谱能量包络的局部极大值中最大的，把它所对应的频率点转换成每分钟的呼吸率，可以得到如图2所示呼吸波形的呼吸率的值。

对于图4所示的呼吸波中含有较大的心动干扰的情况，利用频谱法计算会在其频谱上出现两个较大的峰，其中较大的峰值处所对应的频率为呼吸波的频率，而次大的峰值所对应的频率为心动的频率，分别在图5中用星号标出了呼吸率和用箭头指出了心动干扰所对应的频谱峰值。心电波形和心脏活动引起的干扰在频域比时域容易分离，该图形象地说明了频谱法对于呼吸中含有心动干扰的呼吸率的计算具有明显优势。

对于图6所示的呼吸波中含有躁动干扰的情况，虽然呼吸波形比较混乱，形状也不规则，但是频谱法从频域上非常容易地实现了呼吸率的计算，呼吸率如图7所示，而波形法对于这种呼吸波，计算的误差会比较大。

对于图8所示的呼吸波含有很大干扰的情况，呼吸波形非常紊乱，而且形状很不规则。利用频谱法进行计算时，频谱能量最大值和次大值的频谱能量大小比例相当。但能量最大值所对应的呼吸率很低，根据和历史呼吸率值的比例关系选择次大值所对应的谱峰作为当前正确的谱峰，该结果也比较符合实际，呼吸率的计算结果如图9所示。

通过上面的实例可以说明：当正常呼吸或者呼吸中所含干扰的强度比呼吸信号弱时，频谱法能准确计算出呼吸率；当呼吸中有躁动或者其他较大的干扰时，利用本发明所设定的限制条件，也能较准确地计算出呼吸率。而且本发明中还加入了用历史呼吸率的数据对当前呼吸率进行判断的条件，从而使呼吸率的计算结果更准确和稳定，使呼吸率的计算结果更为可靠。

Claims (16)

1.一种呼吸信号处理方法，包括以下步骤：

A、通过特定的呼吸电路获取呼吸生理信号，AD采集得到呼吸数字信号，并对该信号进行数据预处理；

B、对预处理后的呼吸信号进行频谱转换得到其频域分布；

C、根据呼吸信号的生理参数特征进行是否窒息判断；

D、引入心脏活动生理参数，分析呼吸信号的频谱分布，进行是否受到心脏活动干扰(即心动干扰)的判断；

E、对呼吸波形频域能量包络进行分析，寻找正确的谱峰，将该谱峰对应的频率点转换为呼吸率，并综合历史呼吸率的值得到当前的呼吸率值。

2.根据权利要求1所述的呼吸信号处理方法，其特征是所述的步骤A中数据预处理为对呼吸信号进行带通滤波。

3.根据权利要求2所述的呼吸信号处理方法，其特征是所述的带通滤波为双向零相位IIR滤波器。

4.根据权利要求2所述的呼吸信号处理方法，其特征是所述的带通滤波为无乘法椭圆IIR滤波器。

5.根据权利要求1所述的呼吸信号处理方法，其特征是所述的频谱转换的方法为傅立叶变换。

6.根据权利要求1所述的呼吸信号处理方法，其特征是所述的频谱转换的方法为小波变换。

7.根据权利要求1所述的呼吸信号处理方法，其特征是所述的频谱转换的方法为希尔伯特变换。

8.根据权利要求5所述的呼吸信号处理方法，其特征是所述的傅立叶变换为快速傅立叶变换(FFT)。

9.根据权利要求5所述的呼吸信号处理方法，其特征是所述的傅立叶变换为线性调频Z变换(CZT)。

10.根据权利要求1所述的呼吸信号处理方法，其特征是所述的窒息判断方法为设定呼吸信号频域能量阈值线和呼吸信号时域波形幅度阈值线，呼吸信号在一段时间内频谱能量最大值小于设定的频域能量阈值线且时域波形幅度小于时域波形幅度阈值线则判定呼吸信号窒息。

11.根据权利要求1所述的呼吸信号处理方法，其特征是所述的引入心脏活动参数为心率值。

12.根据权利要求1所述的呼吸信号处理方法，其特征是所述的引入心脏活动参数为脉率值。

13.根据权利要求1所述的呼吸信号处理方法，其特征是所述的心动干扰判定方法为求呼吸频谱能量的最大值和次大值，若其对应的呼吸率值在当前的心率值/或脉率值的设定范围内，则判定呼吸波形受到不同程度的心动干扰。

14.根据权利要求1所述的呼吸信号处理方法，其特征是所述的对呼吸波形频域能量包络分析的方法为对求呼吸信号频谱能量的极大值点确定若干个谱峰。

15.根据权利要求1所述的呼吸信号处理方法，其特征是所述的寻找正确的谱峰为对求得的谱峰进行分析，根据各个谱峰的能量大小比例、各个谱峰对应的呼吸率值和历史呼吸率值的大小比例来确定一个正确的谱峰，并将其对应的频率点转换为呼吸率。

16.根据权利要求1所述的呼吸信号处理方法，其特征是所述求当前呼吸率值的方法为求得的谱峰进行分析，根据各个谱峰的能量大小比例、各个谱峰对应的呼吸率值和历史呼吸率值的大小比例来确定一个正确的谱峰，并将其对应的频率点转换为呼吸率，并将获得的呼吸率值和历史呼吸率值做加权平均。

Images