CN110384490B - 基于相空间的bcg信号心率提取方法 - Google Patents

Abstract

本发明公开了一种基于相空间的BCG信号心率提取方法，利用相空间将BCG时间序列信号从一维数据转换成多维数据，对相空间上的数据，通过SVD进行降维提取主成分投影到二维平面，通过聚类在二维平面上的点来过滤掉噪声，通过对数据Z中剩余的点的波峰查找，计算得到BCG信号X的瞬时心率Z_rr，本发明具有能从更高的维度表达数据，数据信息量更高，处理数据更精准有效，抗形变能力好，弥补了波形法的缺陷，提高了BCG在心率检测上的实用性的特点。

Description

基于相空间的BCG信号心率提取方法

技术领域

本发明涉及心率提取技术领域，尤其是涉及一种抗干扰能力强，抗形变能力强，心率波峰更准确的基于相空间的BCG信号心率提取方法。

背景技术

ECG信号心率检测，心电图(Flectrocardiogram，ECG)：指心脏在每个心动周期中，由起搏点、心房、心室相继兴奋，伴随着生物电的变化，通过心电描记器从体表引出多种形式的电位变化的图形。简单来说心电图就是利用心电图机从体表记录心脏每个心动周期所产生的电活动变化图形的技术。一般通过电极片贴在人身体的各个部位(四肢和胸部)收集ECG信号。

PPG信号心率检测，光电容积脉搏波描记法(PPG)，目前市场上能看到采用这种技术的穿戴设备就有：AppleWatch、三星GalaxyGear S2、Moto 360、Microsoft Band等。

就拿AppleWatch来说，测量心率时底部的表盘会发出绿色的灯光，并且测量的时候手腕最好保持不动否侧会影响测量结果。

BCG信号心率检测，BCG信号是基于PVDF(压电薄膜)传感器的信号采集系统，可以采集到心脏收缩过程中的心冲击图信号(原理是心脏的收缩会将压力传导到压电薄膜传感器上，传感器将压力信号转换为电信号，这就是所谓的心冲击图信号)，我们把这种信号叫做BCG信号。

基于压电薄膜传感器的BCG信号采集设备，只需放置于人的胸部附近位置即可收集到BCG信号，适用场景如晚上睡觉时候，人平躺于床上，传感器放置于床单下胸部附近位置即可。

BCG信号现有的心率提取算法(称为BCG信号波形学心率提取方法，最常见的方案)主要结合了波形学，信号频谱分析，形态学这几个方式的处理手段。先做带通滤波(二阶Butterworth滤波器0.5Hz～3.6Hz)，过滤掉非心率频段的信号数据。通过查找极大值点，找出一些波峰，通过观察调试测试出正常心率的波峰幅值范围，剔除掉幅值过低(环境噪声可能引起过低的小波峰)和过高(比如体动可能引起过高的波峰)的波峰.

现有的技术的缺点如下：

ECG方案，只有去医院检测，需要在身上贴电极片不方便，不能长期监测等。

PPG方案，需要随身佩戴，不方便，对光照敏感，准确度低，鲁棒性不好。

BCG方案(波形法)，对于现实里的BCG信号效果并不是很好，主要原因就是基于波形的分析过分依赖波形，导致其抗干扰能力和抗形变能力较差，当波形模式发生改变时准确率会发生显著的下降，只能在非常规范的场景中使用，而现实中的BCG信号包含体动，呼吸，和很多环境噪声信号，导致了这种波形法方案相对于神经网络的BCG信号心率提取方法准确度低，鲁棒性也不好。

发明内容

本发明的发明目的是为了克服现有技术中的基于波形法的BCG信号心率提取方法上准确度低，鲁棒性差的不足，提供了一种抗干扰能力强，抗形变能力强，心率波峰更准确的基于相空间的BCG信号心率提取方法。

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

一种基于相空间的BCG信号心率提取方法，包括如下步骤：

(1-1)对BCG信号X进行滤波处理：

X＝{x₁，x₂，x₃，…，x_n}，i1＝1，2，3，…，n；

(1-2)对滤波处理后的BCG信号X进行分段；

将BCG信号X分成k段，k＝n/1000，设定分段数据的序号为i2，每一段数据表示为S_i2，1≤i2≤k，则S_i2＝{x_1000i2+1，x_1000i2+2，x_1000i2+3，…，x_1000i2+1000}；

(1-3)找出分段数据S_i2的所有波峰x_l，对x_l进行平滑处理；

(1-4)构建相空间A；

BCG信号X＝{x₁，x₂，x₃，…，x_n}的相空间重构，使用延时坐标法，

相空间重构技术的延时坐标法有两个关键的参数：嵌入的维数d和延迟时间t。

BCG信号的本质是压力传感器受到的压力，将BCG时间序列信号(一维数据)结合相空间转换为多维数据表示，根据力的原理，提取心率，相当于在相空间里提取出心脏跳动产生的力，剔除掉噪声产生的力，在相空间里识别心脏跳动力的作用。

(1-5)相空间A进行数据分段

将相空间A，以60为步进梯度，长度为60构建分段数据集Y，Y＝{Y(1)，Y(2)，…，Y(n_y)}；

Y(i5)＝y(60i5-59)，y(60i5-58)，…，y(60i5)；

其中，

的含义为对

取整；i5＝1，2，3，…，n_y；

(1-6)用奇异值分解SVD法提取主要成分信号

对于Y(i5)做svd分解，提取主要成分信号Y_svd(i5)；

(1-7)对主要成分信号进行数据连接；

对于Y_svd(i5)数据，将每一段数据的首部和上一段数据的尾部连接，得到数据Z，数据Z为n_y×2的矩阵；

将数据Z显示到二维平面上，二维平面上显示出n_y个点；

按照k-近邻进行谱聚类得到3个中心点A，B，C；

将A部分以外的数据点丢弃，数据Z变成n_y1×2的矩阵，0＜n_y1＜n_y；

(1-8)查找波峰，得到BCG信号X的瞬时心率Z_rr。

本发明利用相空间将BCG时间序列信号从一维数据转换成多维数据，对相空间上的数据，通过SVD进行降维提取主成分投影到二维平面，通过聚类在二维平面上的点来过滤掉噪声(短期力和长期力等)的处理。

通过对数据Z中剩余的点的波峰查找，计算得到BCG信号X的瞬时心率Z_rr。

与ECG信号和PPG信号的采集相比，BCG非接触式的测量心率的方法有不用携带，便捷，无创，可以长期监测等优点。

现实环境的复杂性导致了BCG收集的信号包含了体动呼吸和各种环境中的噪声，处理难度特别高，使用BCG信号波形法心率检测效果并不好。

而本发明使用相空间的方法将一维数据转换成相空间中的多维数据，能从更高的维度表达数据，数据信息量更高，处理数据更精准有效，从而可以更准确的剔除噪声，识别波峰点，在准确性，鲁棒性抗干扰，抗形变能力方面都很好。这样既弥补了波形法的缺陷。使得BCG在心率检测上实用性大大增强。

作为优选，步骤(1-1)包括如下步骤：用二阶Butterworth滤波器0.5Hz～3.6Hz对BCG信号进行滤波，过滤掉BCG信号的非心率频段的信号数据，剔除噪声。

作为优选，步骤(1-3)包括如下步骤：

找出波峰x_l后面的第一个波谷x_j，计算j-l，如果j-l＜50，则x_l为小波峰，需要平滑掉；

平滑步骤如下：用0.5，0.1，0.5分别乘以x_l-1，x_l，x_l+1：

x_l-1＝0.5x_l-1，

x_l＝0.1x_l，

x_l+1＝0.5x_l+1。

作为优选，步骤(1-4)包括如下步骤：

(1-4-1)利用自相关系数法计算相空间的延迟时间t；

计算BCG信号X的自相关函数R(t)：

当

时，则t为相空间的延迟时间t；

(1-4-2)确定嵌入维度d＝3；

(1-4-3)构建相空间；

y(i4)＝(x_i4，x_i4+t，x_i4+2t，x_i4+3t，…，x_i4+(d-1)t)

1≤i4≤n-(d-1)t

构建X的相空间A：

将(x₁，x_1+t，x_1+t-x₁)、(x₂，x_2+t，x_2+t-x₂)、(x₃，x_3+t，x_3+t-x₃)，直到(x_n-t，x_n，x_n-x_n-t)分别构成三维向量，得到(n-t)×3的相空间A。

作为优选，步骤(1-6)包括如下步骤：

设定u，sig，v＝svd(Y(i5).T)；

其中，Y(i5).T表示Y(i5)的转置，转置前矩阵Y(i5)为60×3的矩阵，转置后矩阵Y(i).T为3×60的矩阵；

U为3×3的矩阵；

Sig为包括3个元素的列向量；

V为60×60矩阵；

ut_trans＝((Y(i5)×u)×sig).T；

ut_trans是3×60矩阵，ut_trans丢弃第3行数据后，得到2×60矩阵；

Y_svd(i5)＝ut_trans.T；

ut_trans.T表示ut_trans的转置，Y_svd(i5)为Y(i5)对应的主要成分信号投影到二维平面处理后的结果；

作为优选，步骤(1-8)包括如下步骤：

将数据Z按列拆分，得到2个长度为n_y1的1维数据Z₁和Z₂；对Z₁和Z₂分别进行归一化处理，分别查找出Z₁和Z₂上的所有波峰，相邻波峰点之间的距离即为RR间期；

如分量Z₁第一个波峰点x_p序列号是p，后面相邻的波峰点x_q序列号是q；

则距离为

Z₁上的RR间隔集合表示为：

Z_{1_rr}＝{dis₁₁，dis₁₂，dis₁₃，…，diS_1n}

Z₂上的RR间隔集合表示为：

Z_{2_rr}＝{dis₂₁，dis₂₂，dis₂₃，…，diS_2n}

计算RR间隔在Z₁和Z₂上的均值：得到BCG信号X的瞬时心率Z_rr：

因此，本发明具有如下有益效果：处理数据更精准有效，更准确的剔除了噪声，识别的波峰点在准确性、鲁棒性、抗干扰和抗形变能力方面都很好，既弥补了波形法的缺陷，又使BCG在心率检测上实用性大大增强。

附图说明

图1是本发明的一种流程图；

图2是本发明的未滤波的原始BCG数据一种示意图；

图3是本发明的0.5Hz～3.6Hz带通滤波以后的BCG数据一种示意图；

图4是本发明的一种查找波峰的示意图；

图5是本发明的一种k-近邻谱聚类的一种示意图。

具体实施方式

下面结合附图和具体实施方式对本发明做进一步的描述。

如图1所示的实施例是一种基于相空间的BCG信号心率提取方法，包括如下步骤：

(1-1)对BCG信号X进行滤波处理：

X＝{x₁，x₂，x₃，…，x_n}，i1＝1，2，3，…，n；

用二阶Butterworth滤波器0.5Hz～3.6Hz对BCG信号进行滤波，过滤掉BCG信号的非心率频段的信号数据，剔除噪声；

图2是的未滤波的原始BCG数据的示意图，图3是0.5Hz～3.6Hz带通滤波以后的BCG数据的示意图；

(1-2)对滤波处理后的BCG信号X进行分段；

将BCG信号X分成k段，k＝n/1000，设定分段数据的序号为i2，每一段数据表示为S_i2，1≤i2≤k，则S_i2＝{x_1000i2+1，x_1000i2+2，x_1000i2+3，…，x_1000i2+1000}；

将BCG信号X按照每1000个数据分为一段，分成k＝n/1000段，是因为采样频率为125，即1秒钟采样125个数值，那么1000个数值就是8秒的数据；

S₁＝{x₁，x₂，x₃，…，x₁₀₀₀}

S₂＝{x₁₀₀₁，x₁₀₀₂，x₁₀₀₃，…，x₂₀₀₀}

S_k＝{x_k001，x_k002，x_k003，…，x_k000}；

(1-3)找出分段数据S_i2的所有波峰x_l，对x_l进行平滑处理；

找出波峰x_l后面的第一个波谷x_j，计算j-l，如果j-l＜50，则x_l为小波峰，需要平滑掉；

找波峰波谷的办法是一个数x_l减前一个数x_l-1，差值(x_l-x_l-1)是正数；数x_l+1减前一个数x_l，差值(x_l+1-x_l)为负数说明找到了一个波峰，波峰点数值是x_l。

当前面的差值是负数，后面的差值为正数说明找到了一个波谷。只要后一个差值和前一个差值正负号不同，就是找到了一个波峰或者波谷。

如图4所示，(x₁₀₃-x₁₀₂)为正数，(x₁₀₄-x₁₀₃)为正数，(x₁₀₅-x₁₀₄)为正数，(x₁₀₆-x₁₀₅)为负数；前面差值一直为正数，然后找到了第一个差值负数则x₁₀₅即为波峰点。

平滑步骤如下：用0.5，0.1，0.5分别乘以x_l--1，x_l，x_l+1：

x_l--1＝0.5x_l-1，

x_l＝0.1x_l，

x_l+1＝0.5x_l+1。

如图4所示，平滑步骤如下：

x₁₀₄＝0.5×x₁₀₄，

x₁₀₅＝0.1×x₁₀₅，

x₁₀₆＝0.5×x₁₀₆；

(1-4)构建相空间A；

(1-4-1)利用自相关系数法计算相空间的延迟时间t；

计算BCG信号X的自相关函数R(t)：

当

时，则t为相空间的延迟时间t；

(1-4-2)确定嵌入维度d＝3；

BCG信号是时间序列信号，对于每一个值，只包含位置和下一个位置数值及用于代表速度的位置和下一个位置数值之差，只能构建一个OXY的相平面，将相空间的维度d＝3。

(1-4-3)构建相空间；

y(i4)＝(x_i4，x_i4+t，x_i4+2t，x_i4+3t，…，x_i4+(d-1)t)

1≤i4≤n-(d-1)t

构建X的相空间A：

将(x₁，x_1+t，x_1+t-x₁)、(x₂，x_2+t，x_2+t-x₂)、(x₃，x_3+t，x_3+t-x₃)，直到(x_n-t，x_n，x_n-x_n-t)分别构成三维向量，得到(n-t)×3的相空间A。

(1-5)对相空间A进行数据分段

将相空间A，以60为步进梯度，长度为60构建分段数据集Y，Y＝{Y(1)，Y(2)，…，Y(n_y)}；

Y(i5)＝y(60i5-59)，y(60i5-58)，…，y(60i5)；

其中，

的含义为对

取整；i5＝1，2，3，…，n_y；

(1-6)用奇异值分解SVD法提取主要成分信号

对于Y(i5)做svd分解，提取主要成分信号Y_svd(i5)；

设定u，sig，v＝svd(Y(i5).T)；

其中，Y(i5).T表示Y(i5)的转置，转置前矩阵Y(i5)为60×3的矩阵，转置后矩阵Y(i).T为3×60的矩阵；

U为3×3的矩阵；

Sig为包括3个元素的列向量；

V为60×60矩阵；

ut_trans＝((Y(i5)×u)×sig).T；

ut_trans是3×60矩阵，ut_trans丢弃第3行数据后，得到2×60矩阵；

Y_svd(i5)＝ut_trans.T；

ut_trans.T表示ut_trans的转置；

(1-7)对主要成分信号进行数据连接；

对于Y_svd(i5)数据，将每一段数据的首部和上一段数据的尾部连接，得到数据Z，数据Z为n_y×2的矩阵；

将数据Z显示到二维平面上，二维平面上显示出n_y个点；

如图5所示，按照k-近邻进行谱聚类得到3个中心点A，B，C；经过对包含各种运动状况的BCG数据的多次试验，发现A部分为主要信号，B部分为短期力，C部分为长期力；长期力为身体活动、体动、呼吸的分量；短期力为突然的外力变化和轻微的振动；

将A部分以外的数据点丢弃，数据Z变成n_y1×2的矩阵，0＜n_y1＜n_y；

(1-8)查找波峰，得到BCG信号X的瞬时心率Z_rr。

将数据Z按列拆分，得到2个长度为n_y1的1维数据Z₁和Z₂；对Z₁和Z₂分别进行归一化处理，分别查找出Z₁和Z₂上的所有波峰，相邻波峰点之间的距离即为RR间期；

如分量Z₁第一个波峰点x_p序列号是p，后面相邻的波峰点x_q序列号是q；

则距离为

Z₁上的RR间隔集合表示为：

Z_{1_rr}＝{dis₁₁，dis₁₂，dis₁₃，…，dis_1n}

Z₂上的RR间隔集合表示为：

Z_{2_rr}＝{dis₂₁，dis₂₂，dis₂₃，…，dis_2n}

计算RR间隔在Z₁和Z₂上的均值：得到BCG信号X的瞬时心率Z_rr：

应理解，本实施例仅用于说明本发明而不用于限制本发明的范围。此外应理解，在阅读了本发明讲授的内容之后，本领域技术人员可以对本发明作各种改动或修改，这些等价形式同样落于本申请所附权利要求书所限定的范围。

Claims (5)

1.一种基于相空间的BCG信号心率提取方法，其特征是，包括如下步骤：

（1-1）对BCG信号X进行滤波处理：

，i1=1,2,3,…,n；

用二阶Butterworth滤波器0.5Hz_~3.6Hz对BCG信号进行滤波，过滤掉BCG信号的非心率频段的信号数据，剔除噪声；

（1-2）对滤波处理后的BCG信号X进行分段；

将BCG信号X分成k段， k = n/1000，设定分段数据的序号为i2，每一段数据表示为S_i2，1≤i2≤k，则

；

（1-3）找出分段数据S_i2的所有波峰

，对

进行平滑处理；

（1-4）构建相空间A；

（1-5）相空间A进行数据分段

将相空间A，以60为步进梯度，长度为60构建分段数据集Y ,Y={ Y(1)，Y(2)，…,Y(

)}；

Y(i5)

；

其中，

（1-6）用奇异值分解SVD法提取主要成分信号

对于Y(i5)做svd分解，提取主要成分信号Y_svd(i5)；

（1-7）对主要成分信号进行数据连接；

对于Y_svd(i5)数据，将每一段数据的首部和上一段数据的尾部连接，得到数据Z ，

的矩阵；

将数据Z显示到二维平面上，二维平面上显示出

个点；

按照k-近邻进行谱聚类得到3个中心点A、B和C，A为主要信号，B为短期力，C为长期力；将A部分以外的数据点丢弃，

的矩阵，

；

（1-8）查找波峰，得到BCG信号X的瞬时心率

2.根据权利要求1所述的基于相空间的BCG信号心率提取方法，其特征是，步骤（1-3）包括如下步骤：

找出波峰

后面的第一个波谷

，计算j-l，如果j-l＜50，则

为小波峰，需要平滑掉；

平滑步骤如下：用0.5,0.1,0.5分别乘以

，

，

：

=0.5

，

=0.1

，

=0.5

。

3.根据权利要求1所述的基于相空间的BCG信号心率提取方法，其特征是，步骤（1-4）包括如下步骤：

（1-4-1）利用自相关系数法计算相空间的延迟时间t；

计算BCG信号X的自相关函数R(t)：

；

当R(t)≤

时，则t为相空间的延迟时间t；

（1-4-2）确定嵌入维度d = 3；

（1-4-3）构建相空间；

构建

的相空间A：

将(

)、(

)、

直到 (

,

,

-

)分别构成三维向量，得到(n-t)×3的相空间A。

4.根据权利要求1所述的基于相空间的BCG信号心率提取方法，其特征是，步骤（1-6）包括如下步骤：

设定u,sig,v = svd( Y(i5).T )；

其中，Y(i5).T表示Y(i5)的转置，转置前矩阵Y(i5)为60×3的矩阵，转置后矩阵Y(i).T为3×60的矩阵；

U为3×3的矩阵；

Sig为包括3个元素的列向量；

V为60×60矩阵；

；

ut_trans = ( ( Y(i5)× u )× sig ).T；

ut_trans是3×60矩阵，ut_trans丢弃第3行数据后，得到2×60矩阵；

Y_svd(i5) = ut_trans.T ；

ut_trans.T表示ut_trans的转置。

5.根据权利要求1所述的基于相空间的BCG信号心率提取方法，其特征是，步骤（1-8）包括如下步骤：

将数据Z 按列拆分，得到2个长度为

的1维数据

和

；对

和

分别查找出

和

上的所有波峰，相邻波峰点之间的距离即为 RR 间期；

如分量

第一个波峰点

序列号是p，后面相邻的波峰点

序列号是q；

则距离为

秒；

上的RR间隔集合表示为：

上的RR间隔集合表示为：

计算RR间隔在

和

上的均值：得到BCG信号X的瞬时心率

。

Images