CN101732050A - 一种基于光电容积波的呼吸率监测方法 - Google Patents

Abstract

本发明公开了一种基于光电容积波的呼吸率监测方法，其特征在于，包括下述步骤：首先采集被测部位皮肤表层微血管床血液容积变化的光电容积波数字信号，该信号是未经滤波处理、包含呼吸信息的、采样频率在500Hz及以上的原始光电容积波信号；接着对采集得到的光电容积波采用数学形态学方法进行预处理；然后所得的预处理信号采用特征点定位方法检测容积脉搏波每拍的峰值点和谷值点；计算每拍的时间间期变量和每拍峰值点与谷值点强度的差值变量；最后对得到的时间间期变量和幅度变量进行呼吸节拍检测，计算呼吸率。

Description

一种基于光电容积波的呼吸率监测方法

技术领域

本发明涉及一种数字信号处理方法，具体涉及一种基于高采样率的光电容积波进行呼吸率监测的方法。

背景技术

光电容积波描记法(PhotoPlethysmograph，PPG)是应用光电转换技术在皮肤表面测量微血管床血液容积变化的无创检测方法。基于此方法的光电容积波已被广泛应用于血氧饱和度和心率的监测。由于人体不同系统间的相互作用和影响，呼吸运动同样也会引起外周循环的变化，因此能够利用光电容积波信号实现呼吸监测。和现有的呼吸监测方法，如胸阻抗及口鼻气流等法相比，基于光电容积波技术的呼吸监测更加简便和高效，它仅用单一的光电传感器就能同时完成呼吸率，心率以及血氧饱和度等多项指标的监测，在睡眠呼吸障碍等症的早期检测与监护应用方面有着明显的技术优势。

光电容积波信号包含着频率对应于脉率的交流成分和调制于其上的低频基线成分，其中基线中包含由呼吸引起的强度变量已得到共识。以瑞典linkoeping大学相关研究小组为首的各国研究人员自1992年始对从光电容积波中提取呼吸信息进行了大量研究工作，从光电容积波的基线中获得呼吸引起的强度变量并进一步检测呼吸率的算法不断被提出。中国专利ZL03808327.2则就通过对光电容积波的基线信号分析获得呼吸率等生理参数的方法申请了专利保护。上述研究成果的共性都是通过对光电容积波中基线信号的分析获取呼吸信息。从技术实现上则是通过滤波或时频分析技术从光电容积波中提取与呼吸频率范围相对应的信号成分，而后对该成分做进一步的分析处理来达到呼吸率检测的目的。

然而，由于人体系统的复杂性，光电容积波基线的成分并非单一。以成年人为例，不仅包含有一般频率小于0.4Hz的呼吸成分；还包含有与自主神经系统对微血管控制活动相关、频率在0.1Hz左右的Meyer波；以及与温度调节相关、频率在0.04Hz左右的极低频成分。相对而言，极低频的温度漂移较易消除，然而和呼吸频率极为接近甚至相互混叠的Meyer波成分对呼吸波的影响则不容忽视，尤其当两种波形频率混叠的时候，采用滤波方法难以将处于相同频带的两种不同波形分别提取出来。中国专利ZL03808327.2中尽管指出了呼吸波和Meyer各自对应的相位变化特点，在明确获得两种不同波形的条件下，可通过相位关系来判断其为呼吸波或Meyer波，但是该发明却并没有对如何将频率混叠的两个波形分隔开提出明确的技术解决方案。

除Meyer波的干扰外，运动伪差对光电容积波信号的影响同样不容忽视。在光电容积波的采集过程中，轻微的肢体运动甚至极微小的手指等局部运动都会在光电容积波的强度上有所反映。当从基线中提取呼吸引起的强度变量时，运动伪差直接对光电容积波的强度的改变会极大降低呼吸率检测的正确性甚至导致无法检测呼吸率。

在早期的心电信号研究中已证实呼吸对心血管系统的影响包括呼吸对心率的频率调制，即众所周知的呼吸性窦性心律不齐，以及呼吸对每搏心输出量的幅度调制。在光电容积波中同样也存在受呼吸调制的脉率变量以及每拍脉搏的幅度变量(见图2)。各种不同的呼吸调制方式中，脉率和幅度变量比强度变量有着更好的抗干扰能力。具体而言，之前提到Meyer波是微血管收缩和舒张运动在血液容积变化上的反映，也即是说Meyer波主要在强度和振幅变量中得以体现，但是对脉率间期而言，此方面的影响甚微。因此当呼吸频率居于0.05-0.15Hz时，脉率间期的变化能更准确地反映呼吸运动。而当呼吸频率高于此区域时，由于幅度变量是峰值点与谷值点的强度差值，其抗干扰能力比强度变量更优。因此针对不同的呼吸频率选择脉率间期或幅度变量能获得比强度变量更好的呼吸检测效果。

然而现有的研究方法却极大制约了采用脉率和幅度变量进行呼吸检测的进展。前述光电容积波主要包含有基线分量和交流分量，已有的算法通常将两种成分隔离开来，将每种成分独立用于不同的研究目的。一方面通过滤波手段滤除高频的工频干扰和低频的基线信号，滤除与呼吸相关的部分，仅仅利用滤波所得的光电容积脉搏波来完成血氧饱和度和心率检测。另一方面在研究呼吸检测时，认为呼吸信息仅仅包含于通过滤波提取的0.5Hz以下的基线成分中，而将呼吸研究局限于基线成分。因此造成光电容积波中呼吸引起的强度变化成为提取呼吸信息的唯一考虑。

发明内容

本发明针对现有技术中采用光电容积波中强度变量监测呼吸率的局限性，提出了一种使用数学形态学和特征点定位等算法，从高采样率光电容积波提取呼吸引发的脉率间期变量和每拍幅度变量，并由此进行呼吸率监测的方法。该测量方法使用简便无创，价格低廉，非常符合家庭及社区日常健康监护的需要。

为达到以上目的，本发明是采用如下技术方案予以实现的：

一种基于光电容积波的呼吸率监测方法，其特征在于，包括下述步骤：

步骤1，采集被测部位皮肤表层微血管床血液容积变化的光电容积波数字信号，该信号是未经滤波处理、包含呼吸信息的、采样频率在500Hz及以上的原始光电容积波信号；

步骤2：对步骤1采集得到的光电容积波进行预处理，包括采用数学形态学方法分别获取原始光电容积波信号的上下包络，对所得的上下包络求均值，并从原始光电容积波信号中减去该均值；

步骤3：对步骤2所得的预处理信号，采用特征点定位方法检测容积脉搏波每拍的峰值点和谷值点；

步骤4：根据步骤3得到的峰值点和谷值点，计算每拍的时间间期变量和每拍峰值点与谷值点强度的差值变量，即幅度变量；

步骤5：对步骤4得到的时间间期变量和幅度变量进行呼吸节拍检测，计算呼吸率。

上述方案中，步骤2所述对光电容积波进行预处理的具体方法是：设f(n)和B(m)分别为待处理的函数以及形态学算子的结构元素，则形态学算子和运算分别定义为：

侵蚀算子：

$\left(\mathrm{f\ΘB}\right)\left(n\right)=\underset{m=0,...,M-1}{\min }\{f(n-\frac{M-1}{2}+m)-B\left(m\right)\},n=\{\frac{M-1}{2},...,N-\frac{M+1}{2}\}---\left(1\right)$

膨胀算子：

$(f\⊕B)\left(n\right)=\underset{m=0,...,M-1}{\max }\{f(n-\frac{M-1}{2}+m)+B\left(m\right)\},n=\{\frac{M-1}{2},...,N-\frac{M+1}{2}\}---\left(2\right)$

开运算：

闭运算： $f\·B=f\⊕\mathrm{B\ΘB}---\left(4\right)$

采用上述算子和运算符，设输入的原始光电容积波信号为S_Ori(n)，第一个结构元素为B₁(m₁)，第二个结构元素为B₂(m₂)，则光电容积波上包络S_up(n)的求取方法为：

S_up(n)＝S_Ori(n)·B₁(m₁)οB₂(m₂)，                (5)

下包络S_down(n)的求取方法为：

S_down(n)＝S_Ori(n)οB₁(m₁)·B2(m₂)                (6)

上下包络的均值S_Avg(n)表达式为：

S_Avg(n)＝(S_up(n)+S_down(n))÷2                    (7)

从原始光电容积波信号中减去该均值得预处理的光电容积波S_PPG(n)表达式为：

S_PPG(n)＝S_Ori(n)-S_Avg(n)                    (8)

步骤3中所述特征点定位方法是：先进行峰值点检测：为光电容积波中每一采样点设置一个标志位，当前采样点较前点小或保持不变时标志位置为-1，当前采样点较前点大时标志位置为1；找到当前采样点标志位为1而其后一采样点标志位为-1的点，判断其在局部区域是否为最大值，若是则该点判定为峰值点；然后进行谷值点检测：相邻峰值点之间最小点所对应的即为谷值点。

步骤5中对呼吸节拍的检测，采用特征点定位法：为每个呼吸节拍点设置一个标志位，当前节拍点值比前点小或保持不变时标志位置为-1，而当前节拍点值较前点大时标志位置为1；若当前节拍点值的标志位为1，而其后一点为-1，且该点为局部最大值，则设定该点为节拍起始点，找到每个呼吸节拍的起始点，则每个呼吸节拍即可得，呼吸率即为每个呼吸节拍时间间期的倒数。

本发明的发明点在于：兼顾光电容积波的两种主要成分，着重考虑如何从光电容积波中正确提取脉率间期及幅度信息。针对传统滤波技术对光电容积脉搏波产生的相、频失真问题，借鉴图像处理中广泛使用的数学形态学方法，从光电容积波信号的图形特征出发，完成了极低频基线的消除，此方法利用简便的集合运算来获取所需的信号成分，不仅降低了算法实现的计算复杂度，同时也很好地解决了相、频失真的问题。

本发明的优点在于：

使用光电容积波中的频率和幅度变量进行呼吸检测，此两种变量与传统采用的强度变量相比，抗干扰能力更强，提高了呼吸监测的正确性。

使用数学形态学方法针对时域波形的图形特点对信号进行非线性处理，去除不相关的信号成分，能够较好地抑制滤波技术引起的相、频失真，保证所提取的脉率间期变量的正确度。

采用数学形态学信号处理方法，无需采用卷积、积分等复杂的运算操作，而只需加减、集合的交并等简单运算，极大降低了计算复杂度，节省了存储空间，便于单片机实现。

本发明简洁实用，无需人工预处理过程。实施时所需资源少、软硬件开销低。既易于集成于多参数检测的脉搏血氧仪中，也可作为独立的呼吸检测装置，还可以作为中间件处理程序，便于无线检测及网络化应用。其简便的检测方式尤其适合于睡眠呼吸监测及特种行业的生理监测。

附图说明

图1为本发明方法中脉率间期及每拍幅度变量的示意图。图中圆圈所标注的点为光电容积脉搏波的峰值点，方块所标注的点为谷值点，两个相邻的峰值点间的时间间隔t2-t1即为该拍的脉率间期，而一个脉搏周期中峰值点与谷值点的强度差值I2-I1即为该拍的幅度。依次获得的每个脉搏波的时间间期和幅度序列即为光电容积波的脉率间期和每拍幅度变量。

图2为本发明方法的步骤流程示意图。

图3为图2中步骤2-5的具体流程图。

图4为本发明对某检测对象处于坐位平静状态下，左手食指采集得到的30秒钟原始光电容积波进行预处理的相关波形图；图中1)为任意选择的30秒指端原始光电容积波数字输入信号；2)和3)分别为检测得到的该光电容积波信号的下包络和上包络；4)为该光电容积波信号上下包络的均值；5)为减去上下包络均值的光电容积波信号，也即经过预处理的光电容积波信号；

图5为采用本发明对某检测对象坐位平静状态下左手食指采集而得的90s光电容积波进行信号处理的结果。图中六个波形图依次为：1)预处理光电容积波峰值点、谷值点检测结果图；2)提取的脉率间期波形图；3)每拍幅度变量波形图；4)作为对比的参考呼吸信号；5)从脉率间期变量中检测的呼吸率；6)从每拍幅度变量检测的呼吸率。

具体实施方式

以下结合附图及实施例对本发明作进一步的详细描述。

如图2所示，一种基于光电容积波的呼吸率监测方法，包括下述步骤：

步骤1，采集原始数字光电容积波信号，本发明中为通过红光或近红外光投射或反射，反映指端、腕部或前额等部位皮肤表层微血管床血液容积变化的光电容积波信号，具体而言既包括对模拟光电信号进行A/D转换后所得，也包括采用光-频率转换式所得的数字光电容积波信号。要求信号为未经过滤波等处理的采样频率在500Hz及以上的原始光电容积波信号，包含着本发明需检测的呼吸信息。

步骤2，对采集所得原始数字光电容积波信号进行上、下包络检测、对所得的上下包络求均值、并进行从光电容积波中减去该均值的预处理，如图3的流程图所示。

对原始光电容积波信号进行上下包络的检测均须用到基于数学形态学的组合运算。数学形态学中有两种最基本的算子(侵蚀和膨胀算子)和两种派生的运算符(开运算和闭运算)，设f(n)和B(m)分别为待处理的函数以及形态学算子的结构元素，则形态学算子和运算分别定义为：

侵蚀算子：

$\left(\mathrm{f\ΘB}\right)\left(n\right)=\underset{m=0,...,M-1}{\min }\{f(n-\frac{M-1}{2}+m)-B\left(m\right)\},n=\{\frac{M-1}{2},...,N-\frac{M+1}{2}\}---\left(1\right)$

膨胀算子：

$(f\⊕B)\left(n\right)=\underset{m=0,...,M-1}{\max }\{f(n-\frac{M-1}{2}+m)+B\left(m\right)\},n=\{\frac{M-1}{2},...,N-\frac{M+1}{2}\}---\left(2\right)$

开运算：

闭运算： $f\·B=f\⊕\mathrm{B\ΘB}---\left(4\right)$

采用上述算子和运算符，设输入的原始光电容积波信号为S_Ori(n)，第一个结构元素为B₁(m₁)，第二个结构元素为B₂(m₂)

则光电容积波上包络S_up(n)的求取方法为：

S_up(n)＝S_Ori(n)·B₁(m₁)οB₂(m₂)，                (5)

下包络S_down(n)的求取方法为：

S_down(n)＝S_Ori(n)οB₁(m₁)·B2(m₂)                (6)

本发明中B₁(m₁)和B₂(m₂)为长度不同，取值均为零的离散序列。设光电容积波采样率为F(Hz)，则B₁的长度约为1×F，而B₂的长度取为1.5B₁。该长度随脉率的改变略有调整。

上下包络的均值S_Avg(n)表达式为：

S_Avg(n)＝(S_up(n)+S_down(n))÷2                (7)

经过预处理的光电容积波S_PPG(n)表达式为：

S_PPG(n)＝S_Ori(n)-S_Avg(n)                (8)图4中1)

为30秒的采样率为500Hz，从左手食指采集得到的原始光电容积波信号；2)为对原始信号依次采用数学形态学开、闭运算获得的光电容积波信号下包络波形；3)为依次采用数学形态学闭、开运算获得的原始光电容积波信号上包络波形；4)为上下包络的均值；5)为原始光电容积波减去上下包络均值后所得的预处理波形。

步骤3，获取处理后光电容积波每拍的峰值点和谷值点。如图1，点(t1，I1)和(t2，I2)分别为一个脉搏周期中的峰值点和谷值点，每拍时间间期即为(t2-t1)ms，每拍幅度为(I1-I2).步骤2中所得的光电容积波已消除极低频的温度漂移。步骤3中对每拍峰值点和谷值点的检测采用特征点定位法。先进行峰值点检测：为光电容积波中每一采样点设置一个标志位，当前采样点较前点小或保持不变时标志位置为-1，而当前采样点较前点大时置为1；找到当前点标志位为1而其后一点标志位为-1的点，判断其在局部区域是否为最大值，若是则该点判定为峰值点。谷值点检测：相邻峰值点之间最小点所对应的即为谷值点。图1中，红色圆圈标注即为峰值点，红色方框标注即为谷值点。如图5中1)所示为对时长90s，呼吸频率控制在10次/分钟的预处理光电容积波峰值点和谷值点的检测结果，上排红色圆圈为检测到的该序列所有峰值点，下排黑色圆圈为检测得到的该序列所有谷值点。

步骤4中每拍的时间间期为所测得的每拍峰值点和谷值点所对应的时间坐标之差值，而每拍峰值点与谷值点的强度差值即为幅度。时间差值与幅度变量的采样频率与脉率相当。依次计算光电容积波每拍脉搏波的时间间期和幅度值，形成新的时间间期序列和幅度序列，即为从光电容积波中提取的由呼吸引起的脉率间期变量和幅度变量。如图5中2)所示为检测得到的脉率间期变量；3)为检测得到的每拍幅度变量。

步骤5中对呼吸节拍的检测，仍采用特征点定位法。每个呼吸节拍的起始点可假设为：在起始点之前，后一点较前一点的变化趋势都是逐步变小或者保持不变，而在起始点则是比前一点增大。因此为每个点设置一个标志位，当某点值比前点小或保持不变时，该点的标志位置为-1，而较前点大时标志位置为1；若当前点的标志位为1，其后点为-1，且该点为局部最大值，则设定该点为节拍起始点，找到每个节拍的起始点，则每个呼吸节拍即可得，呼吸率即为每个呼吸节拍时间间期的倒数。本发明中采用在0.15Hz以下取每拍时间间期变量检测得到的呼吸率，而在0.15Hz以上则取幅度变量检测得到的呼吸率。图5中为基于脉率间期变量检测得到的呼吸率；6)为基于每拍幅度变量检测得到的呼吸率。图中原始光电容积波是在每分10次的呼吸频率下采集得到，即此时的呼吸频率高于0.15Hz，采用6)所示基于幅度变量的呼吸率检测结果。本算法检测的呼吸率也可计为当前时刻往前一分钟内的呼吸节拍数。也可将监测得到的呼吸波形以及呼吸率以合适的方式显示在显示屏上或保存于存储介质。

本发明方法中的呼吸率算法可直接加载于脉搏血氧仪的单片机程序中，使其能够检测更多生理指标；也可作为专门用于呼吸监测的可佩戴式独立装置；还可保持现有的脉搏血氧仪功能不变，而只通过有线或无线的通信方式将原始光电容积波信号发送至PDA，手机或掌上电脑，呼吸监测算法在以上中间设备的软件平台上来完成，其检测结果既可在上述终端平台上显示，也可通过网络发往上一级监护网络，有利于医疗监护信息网络化管理。

Claims (4)

1.一种基于光电容积波的呼吸率监测方法，其特征在于，包括下述步骤：

步骤1，采集被测部位皮肤表层微血管床血液容积变化的光电容积波数字信号，该信号是未经滤波处理、包含呼吸信息的、采样频率在500Hz及以上的原始光电容积波信号；

步骤2：对步骤1采集得到的光电容积波进行预处理，包括采用数学形态学方法分别获取原始光电容积波信号的上下包络，对所得的上下包络求均值，并从原始光电容积波信号中减去该均值；

步骤3：对步骤2所得的预处理信号，采用特征点定位方法检测容积脉搏波每拍的峰值点和谷值点；

步骤4：根据步骤3得到的峰值点和谷值点，计算每拍的时间间期变量和每拍峰值点与谷值点强度的差值变量，即幅度变量；

步骤5：对步骤4得到的时间间期变量和幅度变量进行呼吸节拍检测，计算呼吸率。

2.如权利要求1所述的基于光电容积波的呼吸率监测方法，其特征在于，步骤2中对光电容积波进行预处理的具体方法是：设f(n)和B(m)分别为待处理的函数以及形态学算子的结构元素，则形态学算子和运算分别定义为：

侵蚀算子：

$\left(\mathrm{f\ΘB}\right)\left(n\right)=\underset{m=0,...,M-1}{\min }\{f(n-\frac{M-1}{2}+m)-B\left(m\right)\},n\{\frac{M-1}{2},...,N-\frac{M+1}{2}\}---\left(1\right)$

膨胀算子：

$(f\⊕B)\left(n\right)=\underset{m=0,...,M-1}{\max }\{f(n-\frac{M-1}{2}+m)+B\left(m\right)\},n=\{\frac{M-1}{2},...,N-\frac{M+1}{2}\}---\left(2\right)$

开运算：

闭运算： $f\·B=f\⊕\mathrm{B\ΘB}---\left(4\right)$

采用上述算子和运算符，设输入的原始光电容积波信号为S_Ori(n)，第一个结构元素为B₁(m₁)，第二个结构元素为B₂(m₂)，则光电容积波上包络S_up(n)的求取方法为：

S_up(n)＝S_Ori(n)·B₁(m₁)оB₂(m₂)，                (5)

下包络S_down(n)的求取方法为：

S_down(n)＝S_Ori(n)оB₁(m₁)·B₂(m₂)                (6)

上下包络的均值S_Avg(n)表达式为：

S_Avg(n)＝(S_up(n)+S_down(n))÷2                    (7)

从原始光电容积波信号中减去该均值得预处理的光电容积波S_PPG(n)表达式为：

S_PPG(n)＝S_Ori(n)-S_Avg(n)                         (8)。

3.如权利要求1所述的基于光电容积波的呼吸率监测方法，其特征在于，步骤3中所述特征点定位方法是：先进行峰值点检测：为光电容积波中每一采样点设置一个标志位，当前采样点较前点小或保持不变时标志位置为-1，当前采样点较前点大时标志位置为1；找到当前采样点标志位为1而其后一采样点标志位为-1的点，判断其在局部区域是否为最大值，若是则该点判定为峰值点；然后进行谷值点检测：相邻峰值点之间最小点所对应的即为谷值点。

4.如权利要求1所述的基于光电容积波的呼吸率监测方法，其特征在于，步骤5中对呼吸节拍的检测，采用特征点定位法：为每个呼吸节拍点设置一个标志位，当前节拍点值比前点小或保持不变时标志位置为-1，而当前节拍点值较前点大时标志位置为1；若当前节拍点值的标志位为1，而其后一点为-1，且该点为局部最大值，则设定该点为节拍起始点，找到每个呼吸节拍的起始点，则每个呼吸节拍即可得，呼吸率即为每个呼吸节拍时间间期的倒数。

Images