发明内容:
本发明的目的,是为了克服现有多导生理记录仪存在价格高昂、操作复杂的缺点,解决传统分析技术如傅立叶分析、自回归分析方法完全失效的问题,提供一种利用ECG信号检测睡眠呼吸暂停的测量方法。
本发明的目的可以通过采取如下措施达到:
利用ECG信号检测睡眠呼吸暂停的测量方法,其特征在于:
利用便携式睡眠呼吸暂停检测记录分析仪对患者进行整夜的监测,记录心电、口鼻气流、血氧饱和度信号,然后利用含有睡眠呼吸暂停的一段心电信号进行如下步骤处理,
1)利用斜率阀值法进行QRS波检测,计算R-R间期时间序列,并利用最小二乘法直线拟合对R-R间期时间序列进行局部检波得到心率变异信号HRV;
2)把HRV随心跳次数的变化作为一种随机信号,选择连续的500~800个(最佳为600个)数据点作为滑动窗口,对窗口内的HRV信号进行经验模态分解EMD,将HRV信号分解为一组内在模态函数IMF;
3)对每个IMF进行Hilbert变换得到HH谱,既HRV信号的幅度和频率的时间分布;
4)根据IMF的HH谱计算各个IMF分量的瞬时频率的平均值、瞬时幅度的标准差、各IMF分量的能量与总能量的比值作为特征值,根据特征值的变化检测出睡眠呼吸暂停的时间和次数。
本发明的目的还可以通过采取如下措施达到:
本发明的一种实施方式是:进行前述第1)步聚时,为了排除QRS波检测过程中因R波漏检而引起的不正常的R-R间期,采用移动平均窗口滤波器滤除这些脉冲噪声,其方法是,取连续的41个R-R间期数据,计算除第21个数据点以外的其它40个数据的平均值,然后用第21个数据点与该局部平均值进行比较,若该数据点大于局部平均值的20%或低于局部平均值的20%,则将其移除;然后改变窗口位置,重复上述计算对所有数据点进行滤波处理。
本发明的一种实施方式是:在前述第1)步中利用最小二乘法直线拟合进行局部检波得到心率变异信号HRV的方法是,取连续的80个数据点组成滑动窗口,进行最小二乘法直线拟合,计算公式如下式:
Si=yi-(axi+b)
其中Si为检波后的HRV信号,yi是心率信号,xi是时间数据,a、b的计算公式如下:
本发明的一种实施方式是:进行前述第2)步操作时,选择连续的600个数据点作为滑动窗口,对窗口内的数据进行经验模态分解EMD,将HRV信号分解为一组内在模态函数IMF;然后对每个IMF进行Hilbert变换得到HH谱,既HRV信号的幅度和频率的时间分布。
本发明的一种实施方式是:
1)把瞬时频率的平均值在0.02Hz—0.055Hz之间的IMF分量作为敏感频段,将相应的瞬时幅度标准差和能量比作为特征值,若瞬时幅度标准差在0.02至0.6之间,能量比大于22.5%则可判为睡眠呼吸暂停;
2)步长为100个数据点,逐步改变窗口的位置,通过敏感频段内的特征值的变化实现呼吸暂停的检测。
本发明具有如下有益效果:
本发明是根据心率信号在呼吸暂停出现时呈现明显的周期性波动,信号某些频带内能量的空间分布与正常呼吸时相比会发生相应变化,分别计算出各频段的功率值,以及各频段功率占总功率的百分比作为特征值,使本不明显的信号特征在某些频带内以显著的能量变化的形式表现出来,采用普通的心电监护技术实现睡眠呼吸暂停的检测,与现行昂贵的睡眠监护仪相比,具有设计简单,使检测对患者睡眠的影响达到最小,对睡眠监护的家庭化以及睡眠呼吸障碍的早期诊断将具有极其重要的意义。
具体实施方式
本实施例利用便携式睡眠呼吸暂停检测记录分析仪对患者进行整夜的监测,记录心电、口鼻气流、血氧饱和度信号,然后利用含有睡眠呼吸暂停的一段心电信号进行上述步骤的处理,图1是经QRS波识别得到的R-R间期时间序列,图2是选择窗口宽度为80检波得到的HRV数据。
图3和图4是同一个体不同时段各600个数据点的HRV信号,通过经验模式分解而得到六阶IMF,图中x(t)表示的是原始的HRV信号,imf1—imf6表示六阶IMF,ref表示的是原始信号减去各阶IMF后的残差。从图3、4可以看出,第一阶IMF包含信号的局部高频成份,随着IMF阶数增高,IMF频率成份逐步降低。和傅立叶方法的全局频率成份分解相比,这种分解方式由局部频率成份决定,可以很好地揭示HRV信号中频率的时变特性。由图3和图4的IMF,通过(14)式可分别计算出两段HRV的HH边际谱如图5、6所示。从HH边际谱图可以看出,在低频段,发生睡眠呼吸暂停时所对应的总幅度值明显高于正常呼吸时的总幅度值,也就是说伴随呼吸暂停而产生的心率的波动,改变了HRV能量的分布。但由于图5、6表征的是整个时间跨度内信号在每个频率点上能量累积的分布情况,为了进一步确定能量值变化较大的频段,需要对各IMF分量进行功率谱分析如图7、图8所示。
利用Hilbert-Huang变换对心率变异信号进行分析和特征值提取,利用特征值的变化实现睡眠呼吸暂停检测。
本发明的理论基础如下:
1、经验模态分解
经验模态分解首先假设所采集的信号数据是由许多基本的内在模态函数(IMF)叠加而成,然后将信号分解成若干个本征模态函数。而每一个内在模态函数代表了一个简单的振动模态,它们或线性或非线性,但必须满足下面两个条件:①在整个信号长度上,极值点和过零点的数目必须相等或者至多相差一个;②在任意时刻,由极大值点定义的上包络线和由极小值点定义的下包络线的平均值为零,也就是说信号的上下包络线对称于时间轴。
在实际信号的处理过程中,完全满足第二个条件是不现实的,所以只要二者的平均值小于一个预先确定的小量即可。根据定义,可以采用如下方法分解函数:
设时间序列为X(t),则:
(1)找出X(t)的所有极大值点和极小值点,将其用三次样条函数分别拟合为原数据序列的上、下包络线。上、下包络线的均值为平均包络线m1,将原序列减去m1便可得到一个去掉低频的新序列h1,即
X(t)-m1=h1 (1)
一般h1不一定是一个平稳序列,为此需要对它重复上述过程。如果h1的平均包络线为m11,则去除该包络线所代表的低频成分后的序列为h11,即
h1-m11=h11 (2)。
重复以上过程,经过k次,直到满足以上两个限定条件,结束分解,使得h1k成为第一个IMF项,即
h1(k-1)-m1k=h1k (3)
令
C1=h1k (4)
C1为从该数据分离出的第一个IMF分量,它包括信号的最小时间尺度即最短周期的模态,令原始信号与C1的差值为剩余信号r1
即X(t)-C1=r1 (5)
由于r1包含了较长周期的组份,可将其视为原始数据并重复以上的过程获得C2,其新的差值为:
r1-C2=r2 (6)
......
rn-1-Cn=rn (7)
当rn为单调序列或者相对原始信号幅度极小可忽略不计时,认为完成提取信号内在模态的过程。
综合以上方程(5)(6)和(7),最后获得:
从方程(8)可知,原始信号可表示为n个内在模态函数Cj和一剩余信号rn之和,rn或是常数量或是一组单调数据。所有的IMF分量经过逆向叠加,最后可以还原回原始数据X(t)。
希尔伯特(Hilbert)变换
希尔伯特变换与其它变换不同,它把信号从时间域仍然变换到时间域。设X(t)为一时间序列,Y(t)是它的希尔伯特变换,即
将X(t)和Y(t)形成一复数,可以得到X(t)对应的解析信号Z(t)
将信号的瞬时幅值a(t)、瞬时相位θ(t)表示如下
其瞬间频率定义式为
对于真实信号,利用HT变换求出其共轭正交分量,然后对实信号进行解析表示,可以求出该信号的三个瞬时特征参数,即瞬时幅度、瞬时相位和瞬时频率,从而实现真正意义上瞬时参数的提取。但是,由瞬时频率的物理意义可知,并不是任意的信号都能用瞬时频率来讨论。当信号满足只含一种振动模态,没有复杂叠加波的情况才可行。也就是说对于含有复杂波系的信号序列,直接HT变换在一定程度上失去了方法上的有效性。为此,必须对含有复杂波系的信号序列进行经验模态分解(EMD),得到一系列频率成分从高到低的内在模态函数(IMFs)分量,再对每一个分量进行HT变换得到HH谱,进而得到边际谱。
2、HH谱及边际谱
对式(12)的每一个IMF作HT变换后累加得
这里省略了残差函数rn,Re表示取实部。表达(13)称为HH谱.
由式(11)(13),信号幅度a(t)与瞬时频率ω(t)都是时间的函数,因此可把幅度显示在频率-时间平面上,即构成了HH幅度谱,HH谱精确地描述了信号的幅值在整段上随频率和时间变化的规律。由于能量可用振幅的平方来描述,因此H(ω,t)也在一定程度上反映了信号能量在频率(或时间)各种尺度上的分布规律。HH谱H(ω,t)确定以后,就可以利用下式对时间积分得到HH边际谱(marginal Hilbert spectrum):
HH边际谱提供了每一个频率值所对应的总幅度值,在统计意义上表征了整个时间跨度内信号在每个频率点上能量累积的分布情况。
HHT方法边界效应的解决办法
EMD通过多次的计算包络线来逐个分解IMF。在每一次的计算过程中,要根据信号的上、下包络来计算信号的局部平均值;上、下包络是由信号的局部极大值和极小值通过三次样条插值(Cubic Spline)算法给出,它的构造影响着EMD的全过程,决定着EMD分解的结果,是HHT变换的关键问题。由于信号两端不可能同时处于极大值和极小值,因而边界处的均值需要估计,在对包络进行样条插值时,要对信号或其极值向外进行延拓,以确保包络线抵达端点,否则要么会因在端点处产生较大摆动严重影响数据的完整性,要么会因在端点处谱能量的扩散,信号由两端向数据的中心污染破坏整段数据,最终使得EMD分解失效。本文选用镜像边界延拓方法,就是根据端点和距端点最近的两个极值之间的大小关系决定镜像延拓的对称中心,当端点的值在离端点最近的两个极值之间时,以离端点最近的极值点为对称点进行镜像延拓。反之,以端点为对称点进行镜像延拓。
本发明的实现步骤如下:
(1)先对采集到的心电图(ECG)信号进行滤波,以滤除基线漂移,然后利用斜率阀值法进行QRS波检测,计算R-R间期时间序列。
(2)为了排除QRS波检测过程中因R波漏检而引起的不正常的R-R间期,我们采用移动平均窗口滤波器滤除这些脉冲噪声。其方法是取连续的41个R-R间期数据,计算除第21个数据点以外的其它40个数据的平均值,然后用第21个数据点与该局部平均值进行比较,若该数据点大于局部平均值的20%或低于局部平均值的20%,则将其移除。然后改变窗口位置,重复上述计算对所有数据点进行滤波处理。
(3)利用最小二乘法直线拟合进行局部检波得到心率变异(HRV)信号。其方法是取连续的80个数据点组成滑动窗口,进行最小二乘法直线拟合,计算公式如(16)式:
Si=yi-(axi+b) (16)
其中Si为检波后的HRV信号,yi是心率信号,xi是时间数据,a、b的计算公式如(17)(18)所示:
(4)把HRV随心跳次数的变化作为一种随机信号,选择连续的600个数据点作为滑动窗口,对窗口内的数据进行经验模态分解(EMD),将HRV信号分解为一组内在模态函数(IMF)。然后对每个IMF进行Hilbert变换得到HH谱,既HRV信号的幅度和频率的时间分布。
(5)根据IMF的HH谱计算各个IMF分量的瞬时频率的平均值、瞬时幅度的标准差、各IMF分量的能量与总能量的比值。
(6)把瞬时频率的平均值在0.02Hz—0.055Hz之间的IMF分量作为敏感频段,其相应的瞬时幅度标准差和能量比作为特征值,步长为100个数据点,逐步改变窗口的位置,如果窗口内有呼吸暂停出现时,敏感频段内的特征值会发生较大变化,通过特征值的变化实现呼吸暂停的检测。
表1:各IMF分量瞬时频率平均值
| IMF1 | IMF2 | IMF3 | IMF4 | IMF5 | IMF6 |
正常呼吸(Hz) | 0.219091 | 0.080622 | 0.039806 | 0.017478 | 0.0084146 | 0.0040864 |
呼吸暂停(Hz) | 0.2139 | 0.079383 | 0.041595 | 0.028977 | 0.062419 | 0.0050818 |
各IMF分量的功率谱反映了信号能量在频率各种尺度上的分布规律,但要准确分析各个不同对象间的差异,还必须从定量的角度来度量。因此,有必要在求得瞬时幅度分布规律的基础上来提取反映睡眠呼吸暂停特征的指标。现定义分解后各层能量为
其中N为每个IMF的数据点数。
分解后各层总能量为
其中Ej为j层的能量,n为分解的层数。定义各层的能量比为:
标准差反映的是总体的发散程度,表征了瞬时幅度a(t)随时间的随机波动或者振动的强度。表2为计算得到的睡眠呼吸暂停和正常呼吸时各IMF的能量比和瞬时幅度的标准差,表中数据基本反映了睡眠呼吸暂停引起心率的波动情况。
表2 睡眠呼吸暂停和正常呼吸时各IMF的能量比和瞬时幅度的标准差
| IMF1 | IMF2 | IMF3 | IMF4 | IMF5 | IMF6 |
呼吸暂停能量比Mj(%) | 7.19 | 20.9 | 30.62 | 19.98 | 15.34 | 5.98 |
呼吸暂停瞬时幅度标准差 | 0.01131 | 0.01985 | 0.02684 | 0.01063 | 0.00605 | 0.004619 |
正常呼吸能量比Mj(%) | 5.17 | 20.97 | 13.01 | 20.98 | 27.74 | 12.14 |
正常呼吸瞬时幅度标准差 | 0.00639 | 0.01714 | 0.01158 | 0.01498 | 0.0118 | 0.002012 |
通过多次实验对比发现,IMF3的平均瞬时频率一般在0.02Hz至0.05Hz之间,能量比和瞬时幅度的标准差在睡眠呼吸暂停出现时产生较大的变化,说明这些特征值能够灵敏地反映睡眠呼吸的变化过程。通过改变滑动窗口的位置,可以迅速准确地进行睡眠呼吸暂停的检测和定位。
需要说明的是,希尔伯特-黄变换(Hilbert-Huang Transform,简称HHT)是继小波变换之后,由美国科学家N.E.Huang等人于1998年提出的又一种主要用于非平稳信号分析的新方法,这种方法不仅可以适用于线性过程的分析,而且适用于非线性和非平稳时间序列的分析。HHT的核心思想是将信号通过经验模态分解(empirical mode decomposition,EMD),分解成数个固有模态函数(Intrinsic Mode Function,IMF)然后利用Hilbert变换构造解析信号,得出信号的瞬时频率和振幅,进而得到Hilbert谱[4]。由于EMD方法是依据信号本身的时域信息进行的时域分解,得到的IMF通常个数是有限和平稳的,而且是具有实际意义的窄带信号,基于这些IFM分量进行的Hilbert变换其结果反映了真实的物理信息,因此,其Hilbert谱也能够准确反映出信号能量、频率在空间或时间尺度上的分布。它不受傅立叶分析要求信号必须是线性、平稳和周期性的局限,又具有小波分析的全部优点,且在分辨率上消除了小波分析的模糊和不清晰,具有更准确的谱结构。因此这种基于EMD的Hilbert频谱分析方法,在非线性和非平稳过程的分析中具有很高的应用价值。