CN103876733B - 用于心肺系统相位同步分析的系统和方法 - Google Patents

Abstract

本发明提供了一种用于心肺系统的相位同步分析的系统和方法，用于评估心率与呼吸同步的程度，即心肺耦合强度；该系统包括心电信号采集模块、心电信号预处理模块、心电信号R波峰检测模块、呼吸信号采集模块、呼吸信号预处理模块、呼吸信号波峰检测模块和相位同步分析模块；其中，心电信号采集模块采集心电信号；心电信号R波峰检测模块采用基于小波变换的R波检测算法检测R波峰；呼吸信号采集模块采集呼吸信号；呼吸信号波峰检测模块采用阈值法检测呼吸波峰点；相位同步分析模块基于相位同步图计算心肺同步时间，作为心肺耦合强度的量化指标。

Description

用于心肺系统相位同步分析的系统和方法

技术领域

本发明涉及一种用于心率与呼吸相位同步分析的系统和方法。

背景技术

心血管系统和呼吸系统受到神经和体液的调控，以各自的节律工作，维持了人体正常的生命活动。但它们之间并不是相互独立没有影响的，而是密切配合，共同完成生理功能。

在中枢的调节和控制下，呼吸活动和心脏活动一样有自己的运动节律，而且这种节律的深度和频率会随着机体内外的环境而改变。同时，节律性呼吸不仅受到中枢的调控，还受到来自呼吸器官本身以及血液循环等其他器官感受器传入冲动的反射性调节，例如化学感受性呼吸反射、肺牵张反射和呼吸肌本体感受性反射。

呼吸过程由外呼吸（包括肺通气和肺换气）、内呼吸（即组织细胞和毛细血管血液之间的气体交换过程）、气体在血液中的运输三个环节组成，呼吸运动通过这三个环节使人体获取新陈代谢所需要的O₂，向外界排出代谢产物CO₂。其中，后两个环节均在血液循环系统中完成，因此，呼吸运动是在神经和体液的调控下以及心血管系统的协调配合下完成的。同时，呼吸运动又会通过神经和体液的调节来影响心血管系统的活动。呼吸性窦性心率不齐是呼吸影响心血管系统活动的很好的例子。

心血管系统受到神经调节、体液调节和自身调节，这些调节是通过各种反射来实现的，调节的因素很多，调节的过程也很复杂。以下以窦性心率不齐为例，描述呼吸过程中，各种影响因素调节心血管系统活动的过程。

人体在做吸气运动时，胸廓扩张，这一活动使支气管和细支气管的平滑肌层中的肺牵张感受器兴奋，发放冲动并将冲动传入到呼吸中枢。中枢将传入的信息进行整合，抑制迷走神经的活动，增强交感神经活性，从而使心率加快，血压升高。

血压升高后，位于颈动脉窦和主动脉弓血管外膜下的动脉压力感受器兴奋，发放冲动并将冲动传入到中枢，中枢再次将传入信息整合后，增强迷走神经的活性，抑制交感神经的活性，使心率减慢，血压降低；与此同时，血压升高引起心房和心室的血容量增多，使位于心房、心室和肺循环大血管壁的心肺感受器兴奋，它们发放的冲动经中枢整合后，进一步降低交感神经活性，增强迷走神经活性，使心率减慢，血压下降。

心率减慢引起血液中的PO₂（氧分压）降低，PCO₂升高，H⁺浓度升高，颈动脉体和主动脉体化学感受器感受到这些变化，发放冲动并传入中枢，使呼吸加快，心率加快。

上述这些感受器的活动与神经中枢共同作用，调节呼吸和心血管系统的活动，使心率与呼吸同步。研究心血管系统和呼吸系统的同步在本领域具有重要的意义。

发明内容

本发明的目的在于提供一种用于心肺系统相位同步分析的系统和方法。

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

一种用于心肺系统的相位同步分析的系统，所述系统包括心电信号采集模块、心电信号预处理模块、心电信号R波峰检测模块、呼吸信号采集模块、呼吸信号预处理模块、呼吸信号波峰检测模块和相位同步分析模块；

心电信号采集模块采集心电信号；

心电信号预处理模块包括基线漂移处理模块和数字滤波模块，数字滤波模块进行低通滤波处理、高通滤波处理及/或50/60HZ的陷波处理；

心电信号R波峰检测模块采用基于小波变换的R波检测算法检测R波峰，基于小波变换的R波检测算法采用了支持紧支集且具有一阶消失矩的二次B样条小波，二次B样条小波的等效滤波器系数如下：

h₁=0.3750,h₂=0.1250,h₃=0.0000

g₁=0.5798,g₂=0.0869,g₃=0.0061

h_k=h_1-k,g_k=-g_1-k,若k>3,h_k=g_k=0,

呼吸信号采集模块采集呼吸信号；

呼吸信号预处理模块对呼吸信号进行滤波；所述呼吸信号预处理模块使用Chebyshev低通滤波器滤除高频噪声和干扰，设置截止频率fc=0.04Hz，通带纹波0.01dB；所述呼吸信号预处理模块在低通滤波的基础上，再对呼吸信号做自适应滤波；使用Widrow-HoffLMS算法来实现自适应滤波器，取N=71，μ=-0.0001，其中N和μ为Widrow-HoffLMS算法中自适应滤波器的参数；

呼吸信号波峰检测模块采用阈值法检测呼吸信号的波峰点；

相位同步分析模块基于相位同步图计算心肺同步时间，作为心肺耦合强度的量化指标。

R波检测过程如下：

1）特征尺度的选择：

选用21-24四个尺度，采用二次B样条小波，把心电信号进行四层分解，得到小波变换成分和光滑信号

2）R波模极大值列的确定：

检测4个尺度上的模极大值列，以减小噪声以及高P波和高T波对R波检测的影响；

通过以下步骤确定R波对应的模极大值点：

（1）从尺度2⁴的小波变换结果中找出大于阈值ε₄的模极大值点，得到这些点的位置集合

（2）从尺度23的小波变换结果中，在邻域内，邻域选为左右各10个点，找出大于阈值ε₃且与处模极大值点同符号的模极大值点，将其位置定为若附近有几个模极大值点，则选最大的一个；若这个点（指所述最大的一个）的幅值小于1.2倍其它几个模极大值点的幅值，则选最靠近的点；若邻域内没有与处模极大值点同符号的模极大值点，则令为0；得到所有这些点位置的集合 $\left\{n_k^3{|}_{k=1...N}\right\};$

（3）使用步骤（2）的方法，找到尺度2²、2¹上的模极大值点位置的集合 $\left\{n_k^2{|}_{k=1...N}\right\}$ 和 $\left\{n_k^1{|}_{k=1...N}\right\};$

其中，不同的尺度采用的不同的幅度阈值{ε₁,ε₂,ε₃,ε₄}是根据最新检测到的小波变换模极大值来刷新的，刷新公式为：

若 $\left|W_{2^j}f\left(n_k^j\right)\right|\≥2A_j^m,$

则 $A_j^{m+1}=A_j^m,$

否则 $A_j^{m+1}=0.875A_j^m+0.125\left|W_{2^j}f\left(n_k^j\right)\right|$

而j=1,2,3,4

其中，代表检测到的小波变换模极大值；

3）奇异点奇异度的计算：

令假设a为Lipschitz指数的上限；令a≈log₂a_j+1(n_k)-log₂a_j(n_k)；通过四层小波分解，可以得到a₁、a₂和a₃；在R波峰点处必定有a₁>0，且通常情况下a₂>0，而且即使a₂<0时，a₁+a₂必定会大于0；对于大多数R波，通常有a₃<0而且a₁+a₂+a₃≤0，而对于高频噪声和剧烈干扰，a₁≤0，a₂≤0，a₃≤0，且a₁+a₂+a₃≤0；因此，从a₁+a₂+a₃的值不能分辨R波、高频噪声和干扰，而a₁+a₂有很好的分辨效果；所以，本算法在检测R波时选用了a₁、a₂，并令若a'>0，则相应的模极大值点是R波峰值点所对应的；若a'突然减小变为负值，则相应的模极大值点一定是噪声或干扰所对应的，应删除相应的模极大值列；

4）去除孤立的和多余的模极大值列：

（1）删除孤立的模极大值列；

R波在每个特征尺度上都对应于一对模极大值列，即正极大值-负极小值对；这两个模极大值点的间距在尺度2¹上比R波的宽度小；设为尺度2¹上的一个正极大值点，为同一尺度上的负极小值点，若和间距大于给定的阈值120ms，则称为孤立模极大值点，应将其删除；

（2）删除多余的模极大值列；

多余的模极大值列可以使用下面的准则来删除：

选择QRS波的能量主要集中的尺度2³上的模极大值来判别；设两个负极小值点分别为Min1和Min2，其幅值分别为A1和A2，而它们与相应的正极大值点的距离分别为L1和L2；

准则1：若则Min2是多余模极小值点；

准则2：若则Min1是多余模极小值点；

准则3：否则，若Min1，Min2在该正极大值的同侧，那么离该正极大值远的是多余模极大值点；若Min1，Min2在该正极大值点的两侧，那么该正极大值点后的那一点为多余点；

5）R波峰点检测：

检测出正极大值-负极大值对的过零点的位置得到R波峰点的位置。

在检测到一个R波之后，把其后200ms内的极值都忽略。

先对前30秒内所检测的RR间期进行平均，得到最近一段时间的平均心动周期T，若本次检测的RR间期大于1.5T，则在此间期内在尺度2³上用0.5ε₃检测模极大值；若此区间内的一对正极大值-负极小值对点之间的间隔小于140ms，则认为有漏检，检测它们之间的零交叉点，并用3点的时移修正，得到重检的R波。

呼吸信号采集模块采用呼吸绑带换能器测量呼吸运动，该呼吸绑带换能器测量呼吸时胸腔或腹腔运动；呼吸信号采集模块包括前置集成模块，用于对所述呼吸信号进行放大、模/数转换等处理。

心电信号采集模块采用一次性氯化银电极作为传感器，包括三个心电导联，心电信号采集模块包括前置集成模块，用于对所述心电信号进行阻抗匹配、放大、模/数转换等处理。

使用基于庞加莱截面的方法获得相角，即以呼吸信号的峰值点作为呼吸周期的分界点，通过在两个峰值点之间作线性插值获得呼吸信号的相位。

在每次出现R尖峰的时刻观察呼吸信号的相位，即可获得相位同步图，即心肺同步图；检测第i个和第i+m个观测点处循环相对相位的差值是否在一个比较小的阈值范围内；如果连续k个点，其中k≥2m，满足下式条件，就认为这段时间两个系统是相位同步的：

$\&Exists;k>1|{\mathrm{\&Psi;}}_{i+m}-{\mathrm{\&Psi;}}_m|<\mathrm{\&epsiv;},i\&Element;\{l,...,l+k-1|0\&le;l\&le;N_r-k+0\}.---\left(6\right)$

式中，N_r是所有观测点的数目，Ψ_i+m,Ψ_m是循环相对相位；

在检测出来的代表相位同步的观测点位置上做标记，根据标记的观测点算出总的相位同步时间。

一种用于心肺系统的相位同步分析的方法，所述方法包括下列步骤：

采集心电信号；

对心电信号进行基线漂移处理、低通滤波处理、高通滤波处理及/或50/60HZ的陷波处理；

采用基于小波变换的R波检测算法检测R波峰，基于小波变换的R波检测算法采用了支持紧支集且具有一阶消失矩的二次B样条小波，二次B样条小波的等效滤波器系数如下：

h₁=0.3750,h₂=0.1250,h₃=0.0000

g₁=0.5798,g₂=0.0869,g₃=0.0061

h_k=h_1-k,g_k=-g_1-k,若k>3,h_k=g_k=0,

采集呼吸信号；

对呼吸信号进行滤波；使用Chebyshev低通滤波器滤除高频噪声和干扰，设置截止频率fc=0.04Hz，通带纹波0.01dB；在低通滤波之后，再对呼吸信号做自适应滤波；使用Widrow-HoffLMS算法来实现自适应滤波器，取N=71，μ=-0.0001，其中N和μ为Widrow-HoffLMS算法中自适应滤波器的参数；

采用阈值法检测呼吸信号的波峰点；

基于相位同步图计算心肺同步时间，作为心肺耦合强度的量化指标。

R波检测过程如下：

1）特征尺度的选择：

选用2¹-2⁴四个尺度。采用二次B样条小波，把心电信号进行四层分解，得到小波变换成分和光滑信号

2）R波模极大值列的确定：

检测4个尺度上的模极大值列，以减小噪声以及高P波和高T波对R波检测的影响；

通过以下步骤确定R波对应的模极大值点：

（1）从尺度24的小波变换结果中找出大于阈值ε₄的模极大值点，得到这些点的位置集合

（2）从尺度2³的小波变换结果中，在邻域内，邻域选为左右各10个点，找出大于阈值ε₃且与处模极大值点同符号的模极大值点，将其位置定为若附近有几个模极大值点，则选最大的一个；若这个点的幅值小于1.2倍其它几个模极大值点的幅值，则选最靠近的点；若邻域内没有与处模极大值点同符号的模极大值点，则令为0；得到所有这些点位置的集合

（3）使用步骤（2）的方法，找到尺度22、21上的模极大值点位置的集合 $\left\{n_k^2{|}_{k=1...N}\right\}$ 和 $\left\{n_k^1{|}_{k=1...N}\right\};$

其中，不同的尺度采用的不同的幅度阈值{ε₁,ε₂,ε₃,ε₄}是根据最新检测到的小波变换模极大值来刷新的，刷新公式为：

若 $\left|W_{2^j}f\left(n_k^j\right)\right|\&GreaterEqual;2A_j^m,$

则 $A_j^{m+1}=A_j^m,$

否则 $A_j^{m+1}=0.875A_j^m+0.125\left|W_{2^j}f\left(n_k^j\right)\right|$

而j=1,2,3,4

其中，代表检测到的小波变换模极大值；

3）奇异点奇异度的计算：

令假设a为Lipschitz指数的上限；令a≈log₂a_j+1(n_k)-log₂a_j(n_k)；通过四层小波分解，可以得到a₁、a₂和a₃；在R波峰点处必定有a₁>0，且通常情况下a₂>0，而且即使a₂<0时，a₁+a₂必定会大于0；对于大多数R波，通常有a₃<0而且a₁+a₂+a₃≤0，而对于高频噪声和剧烈干扰，a₁≤0，a₂≤0，a₃≤0，且a₁+a₂+a₃≤0；因此，从a₁+a₂+a₃的值不能分辨R波、高频噪声和干扰，而a₁+a₂有很好的分辨效果；所以，本算法在检测R波时选用了a₁、a₂，并令若a'>0，则相应的模极大值点是R波峰值点所对应的；若a'突然减小变为负值，则相应的模极大值点一定是噪声或干扰所对应的，应删除相应的模极大值列。

4）去除孤立的和多余的模极大值列：

（1）删除孤立的模极大值列；

R波在每个特征尺度上都对应于一对模极大值列，即正极大值-负极小值对；这两个模极大值点的间距在尺度2¹上比R波的宽度小；设为尺度2¹上的一个正极大值点，为同一尺度上的负极小值点，若间距大于给定的阈值120ms，则称为孤立模极大值点，应将其删除；

（2）删除多余的模极大值列；

多余的模极大值列可以使用下面的准则来删除：

选择QRS波的能量主要集中的尺度2³上的模极大值来判别；设两个负极小值点分别为Min1和Min2，其幅值分别为A1和A2，而它们与相应的正极大值点的距离分别为L1和L2；

准则1：若则Min2是多余模极小值点；

准则2：若则Min1是多余模极小值点；

准则3：否则，若Min1，Min2在该正极大值的同侧，那么离该正极大值远的是多余模极大值点；若Min1，Min2在该正极大值点的两侧，那么该正极大值点后的那一点为多余点；

5）R波峰点检测：

检测出正极大值-负极大值对的过零点的位置得到R波峰点的位置。

在检测到一个R波之后，把其后200ms内的极值都忽略。

先对前30秒内所检测的RR间期进行平均，得到最近一段时间的平均心动周期T，若本次检测的RR间期大于1.5T，则在此间期内在尺度2³上用0.5ε₃检测模极大值；若此区间内的一对正极大值-负极小值对点之间的间隔小于140ms，则认为有漏检，检测它们之间的零交叉点，并用3点的时移修正，得到重检的R波。

采用呼吸绑带换能器测量呼吸运动，该呼吸绑带换能器测量呼吸时胸腔或腹腔运动；采集呼吸信号还包括对所述呼吸信号进行放大、模/数转换等处理。

对呼吸信号进行滤波是使用Chebyshev低通滤波器滤除高频噪声和干扰，设置截止频率fc=0.04Hz，通带纹波0.01dB。

在低通滤波之后，再对呼吸信号做自适应滤波；使用Widrow-HoffLMS算法来实现自适应滤波器，取N=71，μ=-0.0001，其中N和μ为Widrow-HoffLMS算法中自适应滤波器的参数。

采用一次性氯化银电极作为传感器，用三个心电导联采集心电信号；在采集心电信号后对所述心电信号进行阻抗匹配、放大、模/数转换等处理。

使用基于庞加莱截面的方法获得相角，即以呼吸信号的峰值点作为呼吸周期的分界点，通过在两个峰值点之间作线性插值获得呼吸信号的相位。

在每次出现R尖峰的时刻观察呼吸信号的相位，即可获得相位同步图，即心肺同步图；检测第i个和第i+m个观测点处循环相对相位的差值是否在一个比较小的阈值范围内；如果连续k个点，其中k≥2m，满足下式条件，就认为这段时间两个系统是相位同步的：

$\&Exists;k>1|{\mathrm{\&Psi;}}_{i+m}-{\mathrm{\&Psi;}}_m|<\mathrm{\&epsiv;},i\&Element;\{l,...,l+k-1|0\&le;l\&le;N_r-k+0\}.---\left(6\right)$

式中，N_r,是所有观测点的数目，Ψ_i+m,Ψ_m是循环相对相位。

在检测出来的代表相位同步的观测点位置上做标记，根据标记的观测点算出总的相位同步时间。

本发明的有益效果体现在：

本发明提供了一种用于心肺系统的相位同步分析的系统和方法，本发明可以获得心肺同步时间，可以反映心率与呼吸同步的程度，即心肺耦合强度。

附图说明

图1为用于心肺系统的相位同步分析的系统的框图；

图2为原始ECG信号和4层小波分解结果；

图3为R波峰点检测算法流程图；

图4为R波峰点检测检测效果；

图5为零相位滤波的实现；

图6为具有线性相位的自适应滤波器；

图7为同步图的原理；

图8为心肺系统相位同步图。

具体实施方式

下面结合附图和实施例对本发明作详细说明。

图1示出了本发明用于心肺系统的相位同步分析的系统的框图。如图1所示，用于心肺系统的相位同步分析的系统包括心电信号采集模块、心电信号预处理模块、心电信号R波峰检测模块、呼吸信号采集模块、呼吸信号预处理模块、呼吸信号波峰检测模块和相位同步分析模块。

一、心电信号采集模块和心电信号预处理模块

心电信号采集模块采集心电信号。该心电信号采集模块采用一次性氯化银电极作为传感器，包括三个心电导联。心电信号采集模块包括前置集成模块，用于对所述心电信号进行阻抗匹配、放大、模/数转换等处理。

心电信号预处理模块包括基线漂移处理模块和数字滤波模块。数字滤波模块进行低通滤波处理、高通滤波处理及/或50/60HZ的陷波处理。

二、心电信号R波峰检测模块

心电信号R波峰检测模块采用基于小波变换的R波检测算法。该算法采用了支持紧支集且具有一阶消失矩的二次B样条小波。二次B样条小波的等效滤波器系数如下：

h₁=0.3750,h₂=0.1250,h₃=0.0000

g₁=0.5798,g₂=0.0869,g₃=0.0061

h_k=h_1-k,g_k=-g_1-k,若k>3,h_k=g_k=0,

具体的R波检测过程如下：

1）特征尺度的选择：

小波变换将信号分解为不同频段的成分。高频分量和噪声主要落在小尺度上，低频分量和噪声主要落在大尺度上。对于不同的人，ECG信号中QRS波的频谱稍有不同，但能量主要集中在尺度2³和2⁴上，且在2³上能量最大。高频分量较多的QRS波，2²上能量最大，而低频分量较多的QRS波，2⁴上能量最大。因此，本发明中选用2¹-2⁴四个尺度。采用二次B样条小波，把心电信号进行四层分解，得到小波变换成分和光滑信号见图2。

2）R波模极大值列的确定：

从图2可见，每个R波在每个特征尺度上都会生成一对小波变换模极大值，即正极大值-负极小值对。高频噪声会在尺度2¹和2²上生成模极大值，而低频的高P波或高T波会在尺度2³和2⁴上生成模极大值，因此，本发明检测这4个尺度上的模极大值列，以减小噪声以及高P波和高T波对R波检测的影响。

通过以下步骤确定R波对应的模极大值点：

（1）从尺度24的小波变换结果中找出大于阈值ε₄的模极大值点，得到这些点的位置集合

（2）从尺度2³的小波变换结果中，在邻域内（本发明中“邻域”选为左右各10个点），找出大于阈值ε₃且与处模极大值点同符号的模极大值点，将其位置定为若附近有几个模极大值点，则选最大的一个。若这个点的幅值小于1.2倍其它几个模极大值点的幅值，则选最靠近的点。若邻域内没有与处模极大值点同符号的模极大值点，则令为0。得到所有这些点位置的集合

（3）使用步骤（2）的方法，找到尺度2²、2¹上的模极大值点位置的集合 $\left\{n_k^2{|}_{k=1...N}\right\}$ 和 $\left\{n_k^1{|}_{k=1...N}\right\}.$

在该算法中，不同的尺度采用的不同的幅度阈值{ε₁,ε₂,ε₃,ε₄}是根据最新检测到的小波变换模极大值来刷新的，刷新公式为：

若 $\left|W_{2^j}f\left(n_k^j\right)\right|\&GreaterEqual;2A_j^m,$

则 $A_j^{m+1}=A_j^m,$

否则 $A_j^{m+1}=0.875A_j^m+0.125\left|W_{2^j}f\left(n_k^j\right)\right|$

而j=1,2,3,4

其中，代表检测到的小波变换模极大值。这个方法可以保证，检测的QRS波幅值突然增大时不会影响以后的幅值判别。

3）奇异点奇异度的计算：

令假设a为Lipschitz指数的上限。令a≈log₂a_j+1(n_k)-log₂a_j(n_k)。通过四层小波分解，可以得到a₁、a₂和a₃。在R波峰点处必定有a₁>0，且通常情况下a₂>0，而且即使a₂<0时，a₁+a₂必定会大于0。对于大多数R波，通常有a₃<0而且a₁+a₂+a₃≤0，而对于高频噪声和干扰剧烈，a₁≤0，a₂≤0，a₃≤0，且a₁+a₂+a₃≤0。因此，从a₁+a₂+a₃的值不能分辨R波、高频噪声和干扰，而a₁+a₂有很好的分辨效果。所以，本算法在检测R波时选用了a₁、a₂，并令若a'>0，则相应的模极大值点是R波峰值点所对应的。若a'突然减小变为负值，则相应的模极大值点一定是噪声或干扰所对应的，应删除相应的模极大值列。

4）去除孤立的和多余的模极大值列：

运动伪迹和肌电噪声的频带通常与QRS波的频带重叠。因此，应从前面所得到的模极大值列集合中，剔除由伪迹或肌电噪声引入的模极大值列。

（1）删除孤立的模极大值列。

R波在每个特征尺度上都对应于一对模极大值列，即正极大值-负极小值对。这两个模极大值点的间距在尺度2¹上比R波的宽度小。设为尺度2¹上的一个正极大值点，为同一尺度上的负极小值点，若和间距大于给定的阈值（依据经验我们设间距阈值为120ms），则称ni1为孤立模极大值点，应将其删除。

（2）删除多余的模极大值列。

一个R波在每个尺度上只产生一对模极大值点。但一些带噪声的R波，在一对正极大值-负极小值对的邻域（阈值为120ms）内，会产生两对或更多的模极大值列，而其中仅有一对是有用的。多余的模极大值列可以使用下面的准则来删除：

本发明选择QRS波的能量主要集中的尺度2³上的模极大值来判别。设两个负极小值点分别为Min1和Min2，其幅值分别为A1和A2，而它们与相应的正极大值点的距离分别为L1和L2。

准则1：若则Min2是多余模极小值点；

准则2：若则Min1是多余模极小值点；

准则3：否则，若Min1，Min2在该正极大值的同侧，那么离该正极大值远的是多余模极大值点；若Min1，Min2在该正极大值点的两侧，那么该正极大值点后的那一点为多余点。

5）R波峰点检测：

由前面算法我们已经剔除了噪声和干扰以及P波、T波对应的正极大值-负极小值对，并且删除了孤立的和多余的模极大值列，获得了尺度21上ECG信号R波峰点对应的的正极大值-负极小值对。检测出这些正极大值-负极大值对的过零点的位置，就得到了R波峰点的位置。

在我们的算法波中，为了提高检测率，我们还运用了以下两条策略：

1）不应期。

一般人的心率小于300次/分钟，在一次心跳过后的一段时间内不会产生另一次心跳，也就是会有一段不应期。因此，在检测到一个R波之后，把其后200ms内的极值都忽略，可以避免因噪声干扰而造成的误检。

2）反向搜索。

在心率失常或其他情况下，R波幅度和频率会突然变小，导致模极大值点的幅度达不到阈值，导致漏检。在我们的算法中，先对前30秒内所检测的RR间期进行平均，得到最近一段时间的平均心动周期T，若本次检测的RR间期大于1.5T，则在此间期内在尺度2³上用0.5ε₃检测模极大值。若此区间内的一对正极大值-负极小值对点之间的间隔小于140ms，则认为有漏检，检测它们之间的零交叉点，并用3点的时移修正，得到重检的R波。采用这种方法可以减少绝大多数情况下的漏检。

综上，R波峰点检测算法流程如图3，检测效果如图4所示。

三、呼吸信号采集模块、呼吸信号预处理模块和呼吸信号波峰检测模块

呼吸信号采集模块采集呼吸信号。呼吸信号采集模块采用呼吸绑带换能器测量呼吸运动，该呼吸绑带换能器测量呼吸时胸腔腹腔运动。呼吸信号采集模块包括前置集成模块，用于对所述呼吸信号进行放大、模/数转换等处理。

呼吸信号预处理模块对呼吸信号进行滤波。首先使用Chebyshev低通滤波器滤除高频噪声和干扰，设计截止频率fc=0.04Hz，通带纹波0.01dB。由于Chebyshev低通滤波器具有非线性相位，这会造成后面呼吸波峰检测出错，所以这里采用一定措施，实现“零相位”滤波。让信号通过图5所示的结构：

图中最后输出y(n)的傅里叶变换

$\begin{array}{c}Y\left(e^{\mathrm{jw}}\right)=W^{*}\left(e^{\mathrm{jw}}\right)={\left[V\left(e^{\mathrm{jw}}\right)H\left(e^{\mathrm{jw}}\right)\right]}^{*}={\left[U^{*}\left(e^{\mathrm{jw}}\right)H\left(e^{\mathrm{jw}}\right)\right]}^{*}\\ =U\left(e^{\mathrm{jw}}\right)H^{*}\left(e^{\mathrm{jw}}\right)=X\left(e^{\mathrm{jw}}\right)H^{*}\left(e^{\mathrm{jw}}\right)=X\left(e^{\mathrm{jw}}\right){\left|H\left(e^{\mathrm{jw}}\right)\right|}^2\end{array}---\left(1\right)$

可见输出和输入有着相同的相位。

经过低通滤波后，滤除了呼吸信号的高频噪声，信号变得更光滑。但是，对于大骑波的抑制效果并不是很好。为了保证了后面呼吸波峰检测的准确性，在低通滤波的基础上，再对呼吸信号做自适应滤波。

设计如图6的自适应滤波器，使用Widrow-HoffLMS算法来实现自适应滤波器，获得的结果y(n)比输入x(n)延迟(N-1)/2个单位。取N=71，μ=-0.0001，N和μ为Widrow-HoffLMS算法中自适应滤波器的参数。将经过切比雪夫滤波器滤波的呼吸信号，经过图6给出的自适应滤波器滤波后，再通过消除时延程序运算，得到平滑的呼吸波形。

呼吸信号波峰检测模块采用阈值法检测呼吸波峰点。

四、相位同步分析模块

1.相位同步定义

相位同步定义为：

在上式（2）中，n、m是整数，Φ₁、Φ₂分别是两个振荡器的相位，定义在整个实轴上。是相对相位，const是正的足够小的常量。

如果在一段时间内，上式（2）成立，则认为Φ₁和Φ₂在这段时间T内发生了n:m锁相，即Φ₁和Φ₂之间出现了n:m的相位同步，此时|nΦ₁-mΦ₂|在某个平均值附近小范围地波动。

有时，在受干扰系统中，相对相位会以2π的整数倍移动。对此，可以这样理解，两个振荡器的相位差在一个比较稳定的值周围小范围波动，但是受到噪声干扰后，从一个物理上比较稳定的状态跳变到了另一个比较稳定的状态（i是整数），并在这个新的稳定状态周围小范围波动，直到发生下一次移动。对此，定义循环相对相位(cyclicrelativephase)：

2.相位同步图

由相位同步的定义可知，计算相对相位并得出循环相对相位Ψ_n,m，以时间为横轴画出循环相对相位Ψ_n,m，当两个系统的相位以n:m比率同步时，会有一段水平的稳定期。出现水平稳定期的地方，两系统就处于相位同步状态，并可以计算同步时间。但是这个方法存在两个问题：需要通过尝试才能获得整数n和m，而没有通用的方法；当两个系统之间存在多种同步模式时，不能在一个图中体现出来，也不能发现同步状态之间的跳跃。

相位同步图可以克服这些缺陷。这个方法采用一个“相位频闪观测仪”，在第一个振荡器的循环相位达到特定值θ的时间t_k去观察第二个振荡器的相位φ₂(t_k)，即当φ₁(t_k)mod2π=θ时，相对t_k画出ψ(t_k)，得到同步图，其中：

$\mathrm{\&psi;}\left(t_k\right)=\frac{1}{2\mathrm{\&pi;}}\left[{\mathrm{\&phi;}}_2\left(t_k\right)\mathrm{mod}2\mathrm{\&pi;}\right]---\left(4\right)$

在第一个振荡器相位到达θ时，且没有噪声干扰的情况下，如果两个系统出现n:1同步，同步图上就会出现n条水平线。如果有噪声，水平线会被扰乱，但是当两系统同步时依然能观察到水平带。

对于n:m比率的同步，将第二个振荡器相邻的m个振荡合并为一个周期，这些新周期的每一个周期内的相位从0增加到2πm，见式（5），这样就可以分析n:m比率的同步：

${\mathrm{\&psi;}}_m\left(t_k\right)=\frac{1}{2\mathrm{\&pi;}}\left[{\mathrm{\&phi;}}_2\left(t_k\right)\mathrm{mod}2\mathrm{\&pi;m}\right]---\left(5\right)$

图7给出了同步图的原理。

在心肺耦合研究中，ECG信号的R尖峰具有奇异性，可以以R尖峰作为特征点，将ECG信号简化为点过程。然后，在每次出现R尖峰的时刻观察呼吸信号的相位，即可获得心肺同步图（cardiorespiratorysynchrogram,CRS）。

3.同步图的量化

对同步图的量化的主要目的是检测同步图中的水平线结构，以此量化同步强度。我们假设相位同步图包含n条水平线，两系统以m:n同步模式同步。那么，在相位同步图出现水平线的时间段，第二个振荡器的每个观测周期中（m个周期），都会有n个观测点，且对应的第i个观测点处的循环相对相位都近似相等。这样才会出现n条水平线。因此，检测第i个和第i+m个观测点处循环相对相位的差值是否在一个比较小的阈值范围内。如果连续k个点，一般情况下取k≥2m，满足下面式（6）的条件，就认为这段时间两个系统是相位同步的：

$\&Exists;k>1|{\mathrm{\&Psi;}}_{i+m}-{\mathrm{\&Psi;}}_m|<\mathrm{\&epsiv;},i\&Element;\{l,...,l+k-1|0\&le;l\&le;N_r-k+0\}.---\left(6\right)$

式中，N_r是所有观测点的数目。

在检测出来的代表相位同步的观测点位置上做标记，就可以计算出总的相位同步时间。相位同步时间是反映两个弱耦合系统耦合强度的指标。

4.基于庞加莱截面的方法获得相角

基于庞加莱截面的方法是，选择一个与混沌吸引子的所有轨线横截相交庞加莱截面，那么，连续两个庞加莱截面之间的每一个混沌吸引子可以定义相位。混沌系统的相位是时间的线性函数，可以用线性插值的方法计算相位，若混沌吸引子连续与庞加莱截面相交，则相位增加2π，下式（7）是线性迭代的公式：

$\mathrm{\&phi;}\left(t\right)=2\mathrm{\&pi;}\frac{t-t_n}{t_{n+1}-t_n}+2\mathrm{\&pi;n},t_n\&le;t\&le;t_{n+1}---\left(7\right)$

其中，t_n是第n次与庞加莱截面相交的时间。

5.基于同步图的心肺耦合分析

基于上述理论，以呼吸信号的峰值点作为呼吸周期的分界点，通过在两个峰值点之间作线性插值获得呼吸信号的相位。根据相位同步图原理，以快信号心电信号的每个R波峰值点为观测点观测慢信号呼吸信号的相位，获得循环相对相位。循环相对相位对2π做归一化处理后，就得到了相对相位差ψ_m(t_k)。由公式（4）和（5）可知，ψ_m(t_k)相当于ECG信号的R尖峰相对于相应呼吸周期开始点的相对距离因此可以用来替代ψ_m(t_k)：

其中R_i表示第i个R波峰值点的时间点，I_j表示第j次呼吸波峰的时间点。

根据计算出来的以时间为横轴，画出相位同步图，如图8所示。从图中可以看到在20-30s、100s和150s时间段内出现了4条水平线，说明受试者的呼吸系统和心血管系统以4:1的同步。计算同步图中出现水平线结构的时间段总和，就可以用心肺同步时间作为心肺耦合强度的量化指标。

Claims (8)

1.一种用于心肺系统的相位同步分析的系统，其特征在于：所述系统包括心电信号采集模块、心电信号预处理模块、心电信号R波峰检测模块、呼吸信号采集模块、呼吸信号预处理模块、呼吸信号波峰检测模块和相位同步分析模块；

心电信号采集模块采集心电信号；

心电信号预处理模块包括基线漂移处理模块和数字滤波模块，数字滤波模块进行低通滤波处理、高通滤波处理及/或50/60HZ的陷波处理；

心电信号R波峰检测模块采用基于小波变换的R波检测算法检测R波峰，基于小波变换的R波检测算法采用了支持紧支集且具有一阶消失矩的二次B样条小波，二次B样条小波的等效滤波器系数如下：

h₁＝0.3750,h₂＝0.1250,h₃＝0.0000

g₁＝0.5798,g₂＝0.0869,g₃＝0.0061

h_b＝h_1-b,g_b＝-g_1-b,若b>3,h_b＝g_b＝0,

呼吸信号采集模块采集呼吸信号；

呼吸信号预处理模块对呼吸信号进行滤波；所述呼吸信号预处理模块使用Chebyshev低通滤波器滤除高频噪声和干扰，设置截止频率fc＝0.04Hz，通带纹波0.01dB；所述呼吸信号预处理模块在低通滤波的基础上，再对呼吸信号做自适应滤波；使用Widrow-HoffLMS算法来实现自适应滤波器，取L＝71，μ＝-0.0001，其中L和μ为Widrow-HoffLMS算法中自适应滤波器的参数；

呼吸信号波峰检测模块采用阈值法检测呼吸信号的波峰点；

相位同步分析模块基于相位同步图计算心肺同步时间，作为心肺耦合强度的量化指标；在每次出现R波峰的时刻观察呼吸信号的相位，即可获得相位同步图，即心肺同步图；检测第i个和第i+m个观测点处循环相对相位的差值是否在一个比较小的阈值范围内；如果连续d个点，其中d≥2m，满足下式条件，就认为连续d个点对应的时间两个系统是相位同步的：

$\begin{array}{cc}\&Exists;d>1|{\mathrm{\&Psi;}}_{i+m}-{\mathrm{\&Psi;}}_m|<\mathrm{\&epsiv;},& i\&Element;\{l,...,l+d-1|0\&le;l\&le;N_r-d+1\}\end{array}.---\left(6\right)$

式中，N_r是所有观测点的数目，Ψ_i+m,Ψ_m是循环相对相位；

在检测出来的代表相位同步的观测点位置上做标记，根据标记的观测点算出总的相位同步时间，即心肺同步时间。

2.如权利要求1所述的用于心肺系统的相位同步分析的系统，其特征在于R波检测过程如下：

1)特征尺度的选择：

选用2¹-2⁴四个尺度，采用二次B样条小波，把心电信号进行四层分解，得到小波变换成分其中j＝1,2,3,4和光滑信号

2)R波模极大值列的确定：

检测4个尺度上的模极大值列，以减小噪声以及高P波和高T波对R波检测的影响；

通过以下步骤确定R波对应的模极大值点：

(1)从尺度2⁴的小波变换结果中找出大于阈值ε₄的模极大值点，得到这些点的位置集合

(2)从尺度2³的小波变换结果中，在邻域内，邻域选为左右各10个点，找出大于阈值ε₃且与处模极大值点同符号的模极大值点，将其位置定为若附近有几个模极大值点，则选最大的一个；若这个点的幅值小于1.2倍其它几个模极大值点的幅值，则选最靠近的点；若邻域内没有与处模极大值点同符号的模极大值点，则令为0；得到所有选择的点位置的集合

(3)使用步骤(2)的方法，找到尺度2²、2¹上的模极大值点位置的集合 $\left\{n_k^2{|}_{k=1...N}\right\}$ 和 $\left\{n_k^1{|}_{k=1...N}\right\};$

其中，不同的尺度采用的不同的幅度阈值{ε₁,ε₂,ε₃,ε₄}是根据最新检测到的小波变换模极大值来刷新的，刷新公式为：

若 $\left|W_{2^j}f\left(n_k^j\right)\right|\&GreaterEqual;2A_j^m,$

则 $A_j^{m+1}=A_j^m,$

否则 $A_j^{m+1}=0.875A_j^m+0.125\left|W_{2^j}f\left(n_k^j\right)\right|$

而 $\begin{array}{cc}{\mathrm{\&epsiv;}}_j=0.3A_j^{m+1}& j=1,2,3,4\end{array}$

其中，代表检测到的小波变换模极大值；

3)奇异点奇异度的计算：

令假设a为Lipschitz指数的上限；令a≈log₂a_j+1(n_k)-log₂a_j(n_k)；通过四层小波分解，可以得到a₁、a₂和a₃；在R波峰点处必定有a₁>0，且通常情况下a₂>0，而且即使a₂<0时，a₁+a₂必定会大于0；对于大多数R波，通常有a₃<0而且a₁+a₂+a₃≤0，而对于高频噪声和剧烈干扰，a₁≤0，a₂≤0，a₃≤0，且a₁+a₂+a₃≤0；因此，从a₁+a₂+a₃的值不能分辨R波、高频噪声和干扰，而a₁+a₂有很好的分辨效果；所以，本算法在检测R波时选用了a₁、a₂，并令若a'>0，则相应的模极大值点是R波峰值点所对应的；若a'突然减小变为负值，则相应的模极大值点一定是噪声或干扰所对应的，应删除相应的模极大值列；

4)去除孤立的和多余的模极大值列：

(1)删除孤立的模极大值列；

R波在每个特征尺度上都对应于一对模极大值列，即正极大值-负极小值对；这两个模极大值点的间距在尺度2¹上比R波的宽度小；设为尺度2¹上的一个正极大值点，其中k≠i，为同一尺度上的负极小值点，若和其中k≠i间距大于给定的阈值120ms，则称为孤立模极大值点，应将其删除；

(2)删除多余的模极大值列；

多余的模极大值列使用下面的准则来删除：

选择QRS波的能量主要集中的尺度2³上的模极大值来判别；设两个负极小值点分别为Min1和Min2，其幅值分别为A1和A2，而它们与相应的正极大值点的距离分别为L1和L2；

准则1：若则Min2是多余模极小值点；

准则2：若则Min1是多余模极小值点；

准则3：否则，若Min1，Min2在该正极大值的同侧，那么离该正极大值远的是多余模极大值点；若Min1，Min2在该正极大值点的两侧，那么该正极大值点后的那一点为多余点；

5)R波峰点检测：

检测出正极大值-负极大值对的过零点的位置得到R波峰点的位置。

3.如权利要求2所述的用于心肺系统的相位同步分析的系统，其特征在于在检测到一个R波之后，把其后200ms内的极值都忽略。

4.如权利要求3所述的用于心肺系统的相位同步分析的系统，其特征在于先对前30秒内所检测的RR间期进行平均，得到最近一段时间的平均心动周期T，若本次检测的RR间期大于1.5T，则在此间期内在尺度2³上用0.5ε₃检测模极大值；若此区间内的一对正极大值-负极小值对点之间的间隔小于140ms，则认为有漏检，检测该一对正极大值-负极小值对点之间的零交叉点，并用3点的时移修正，得到重检的R波。

5.如权利要求1-4中任一项所述的用于心肺系统的相位同步分析的系统，其特征在于呼吸信号采集模块采用呼吸绑带换能器测量呼吸运动，该呼吸绑带换能器测量呼吸时胸腔或腹腔运动；呼吸信号采集模块包括前置集成模块，用于对所述呼吸信号进行放大、模/数转换处理。

6.如权利要求5中所述的用于心肺系统的相位同步分析的系统，其特征在于，该心电信号采集模块采用一次性氯化银电极作为传感器，包括三个心电导联，心电信号采集模块包括前置集成模块，用于对所述心电信号进行阻抗匹配、放大、模/数转换处理。

7.如权利要求6中所述的用于心肺系统的相位同步分析的系统，其特征在于，使用基于庞加莱截面的方法获得相角，即以呼吸信号的峰值点作为呼吸周期的分界点，通过在两个峰值点之间作线性插值获得呼吸信号的相位。

8.一种用于心肺系统的相位同步分析的方法，其特征在于所述方法包括下列步骤：

采集心电信号；

对心电信号进行基线漂移处理、低通滤波处理、高通滤波处理及/或50/60HZ的陷波处理；

采用基于小波变换的R波检测算法检测R波峰，基于小波变换的R波检测算法采用了支持紧支集且具有一阶消失矩的二次B样条小波，二次B样条小波的等效滤波器系数如下：

h₁＝0.3750,h₂＝0.1250,h₃＝0.0000

g₁＝0.5798,g₂＝0.0869,g₃＝0.0061

h_b＝h_1-b,g_b＝-g_1-b,若b>3,h_b＝g_b＝0,

采集呼吸信号；

对呼吸信号进行滤波；使用Chebyshev低通滤波器滤除高频噪声和干扰，设置截止频率fc＝0.04Hz，通带纹波0.01dB；在低通滤波之后，再对呼吸信号做自适应滤波；使用Widrow-HoffLMS算法来实现自适应滤波器，取L＝71，μ＝-0.0001，其中L和μ为Widrow-HoffLMS算法中自适应滤波器的参数；

采用阈值法检测呼吸信号的波峰点；

基于相位同步图计算心肺同步时间，作为心肺耦合强度的量化指标；在每次出现R波峰的时刻观察呼吸信号的相位，即可获得相位同步图，即心肺同步图；检测第i个和第i+m个观测点处循环相对相位的差值是否在一个比较小的阈值范围内；如果连续d个点，其中d≥2m，满足下式条件，就认为连续d个点对应的时间两个系统是相位同步的：

$\begin{array}{cc}\&Exists;d>1|{\mathrm{\&Psi;}}_{i+m}-{\mathrm{\&Psi;}}_m|<\mathrm{\&epsiv;},& i\&Element;\{l,...,l+d-1|0\&le;l\&le;N_r-d+1\}\end{array}.---\left(6\right)$

式中，N_r是所有观测点的数目，Ψ_i+m,Ψ_m是循环相对相位；

在检测出来的代表相位同步的观测点位置上做标记，根据标记的观测点算出总的相位同步时间，即心肺同步时间。