基于变分模态分解的脉率变异性信号提取方法

杨文龙，李 鑫，丑丽娟,2，杨海萍，丑永新*

(1.常熟理工学院 电气与自动化工程学院， 江苏 苏州 215500；2.东北石油大学 计算机与信息技术学院，黑龙江 大庆 163711)

0 引言

人体心脏的收缩与舒张驱动血液在心血管系统周而复始地循环，引起人体皮肤浅表处(如手指、腕部、耳垂等[1])脉管的搏动，这种搏动称为脉搏。脉搏信号中含有人体心血管系统丰富的生理和病理信息，诊脉是中医进行疾病诊断的重要依据之一。从脉搏信号中提取的脉率[2]，是评价人体健康状态的重要指标之一。人体心脏搏动规律随着人体状态的变化而不断调整，因此，脉率也在不断发生变化。脉率之间的这种变化称为脉率变异性(Pulse Rate Variability,PRV)，PRV是评价人体心血管系统和神经系统的重要依据。研究表明，从PRV中提取参数可用于心源性猝死、冠心病、心力衰竭、高血压及呼吸暂停综合症等疾病的诊断[3-6]，具有十分重要的临床应用价值。

目前，已有多种PRV信号提取方法，主要可分为3种：① 基于脉搏波时域征的提取方法，这类方法通过设置幅度、时间或者差分阈值来提取脉搏波波峰或者波谷等特征[7]，进而计算脉搏信号的主波间期。这类方法的速度快，准确性高，但是抗噪性能差。② 基于脉搏信号频谱的PRV估计法[8]，从脉搏信号的频谱成分中估计PRV，抗噪性很强，但是误差太大。③ 从脉搏信号的分解分量中提取。这类方法先通过傅里叶级数[9]、小波分解[10]、EMD分解[11]、形态学滤波[12-13]等方法将脉搏信号分解成不同频段的分量，再从低频分量中提取PRV信号。这种方法的准确性和抗噪性很高，但是分量的准确提取选择决定着PRV信号的提取精度，傅里叶级数的窗宽、小波分解的频段、形态学滤波法结构元素的长度等因素选择主观，使其在实际应用中受到限制。EMD算法可实现脉搏波的自适应分解，然而其存在模态混叠。

近年来，有学者结合维纳滤波、Hilbert变换和频率混合变分问题的求解提出了变分模态分解(Variational Mode Decomposition，VMD)算法，在震动[14]、气象[15]、电力[16]、生物医学等信号处理方面表现出优越性能。因此，本文以VMD算法为核心，提出一种从脉搏信号中提取PRV信号的新算法。首先，通过VMD算法对脉搏信号分解，从低频分量中获取脉搏信号主波定位的参考信号；然后，以此为依据，对脉搏信号的波谷进行定位，计算PRV信号。设计实验，通过实测脉搏信号对所提出算法的准确性进行评估。

1 算法原理

1.1 VMD分解算法

VMD分解算法的目的是将原始信号分解为多个本征模态函数(Intrinsic Mode Function，IMF)的和，每个IMF分量为特定带宽的调幅-调频信号。分解后第k个IMF分量的定义如下：

imfk(t)=ak(t)cos(φk(t))，

(1)

式中，ak(t)为第k个IMF分量的包络，φk(t)为其相位，则相位对应的角频率为ωk(t)=φ′k(t)。

对于给定的脉搏信号ppg(t)，可以分解成K个IMF分量之和，每个IMF分量对应的约束变分模型为：

(2)

式中，{imfk}={imf1,imf2,…,imfk}为VMD分解得到所有IMF的集合，{ωk}={ω1,ω2,…,ωk}为各IMF的中心角频率集合。

对于式(2)所示问题的最优解问题，引入拉格朗日函数后，可转化为非约束变分问题，即

(3)

(4)

式中，ppg为时域表达式，PPG为其在频域内的表达式。式(4)经过Parseval傅里叶等距变换后，可以得到：

(5)

式中，Λn+1(ω)的更新方式为：

(6)

同理可得：

(7)

在此基础上，VMD分解分为以下几步：

② 根据式(5)～式(7)更新imfk、ωk和λ；

③ 重复前两步，直到满足以下约束条件：

(8)

1.2 基于VMD分解的PRV信号提取原理

脉搏信号经过VMD分解后，可得到一系列IMF分量，不同的IMF分量中包含脉搏信号不同的频段。其中，imf1中含有脉搏信号的基波分量，基波分量的主要成分为脉搏信号的主波，即图1中的虚线信号。该信号不含脉搏信号的高频噪声、重搏波和潮波等波谷定位的干扰信息，相比直接从脉搏信号中定位波谷的难度明显降低。同时，imf1能自适应地随着脉搏信号趋势改变而变化，其两个峰值(图1中黑色圆圈)之间的最小值正好对应脉搏信号的波谷。因此，imf1分量可作为PRV信号检测的重要依据。

图1 从imf1分量中定位波谷的示意图Fig.1 Process of locating troughs fromimf1 component

依据图1所示的原理，提出图2所示的PRV信号检测流程。

图2 算法流程Fig.2 Process of proposed method

首先，对脉搏信号进行VMD分解得到imf1分量；对imf1分量进行20点平滑滤波，去除残余高频噪声的干扰；然后，计算imf1分量的极大值，存储极大值位置,以相邻极大值位置为参考，在原始脉搏信号中计算最小值位置即为脉搏信号波谷位置；对波谷位置进行一阶差分，计算PRV信号，PRV信号计算过程见式(9);最后，对PRV信号的提取结果进行评估。

(9)

式中，fs为脉搏信号采样频率，PPI(j)为由波谷位置一阶差分得到的脉搏信号间期，PRV(j)为第j个PRV信号样本点。

2 实验数据

采用自主研制的可穿戴式生理信号采集系统获取实验数据，采集系统如图3所示。系统由上位机和下位机组成，上下位机通过Zigbee组网进行无线传输。下位机分为两个模块，分别为放置在胸部的心电采集模块和放置在腕部的脉搏信号采集模块。上位机实时接收下位机数据，对其进行存储与处理。

图3 脉搏信号采集过程Fig.3 Recording process of pulse signal

招募年龄为21～26岁的健康在校大学生(男生18名，女生13名)作为实验对象。数据实验前在常熟市第一人民医院进行体检，确保实验对象的健康状态。在常熟第一人民医院伦理委员会的监督之下，实验对象佩戴实验设备，处于静坐状态，采集各路信号4 min，采样频率为1 000 Hz。

3 实验结果及分析

对实验数据按照图2流程处理，其中，VMD分解的参数配置为：分解层数为5，τ=0.001，迭代停止误差ε=0.1。对于如图4所示的脉搏信号，进行VMD分解后得到IMF分量，如图5所示，分解后得到5个IMF分量。相比于图4，imf1分量含有脉搏信号的主波成分，imf2、imf4和imf5中含有噪声成分，而imf3中含有重搏波和潮波成分。

图4 原始脉搏信号Fig.4 A raw pulse signal

(a)imf1

(b)imf2

(c)imf3

(d)imf4

(e)imf5图5 脉搏信号的VMD分解结果Fig.5 VMD decomposition results of a raw pulse signal

因此，选用imf1分量按照图2所示流程计算PRV信号，结果如图6所示。可以看出，只有第一个波谷没有定位出来，其他波谷都可以准确定位，这是由于要从imf1分量的两个极大值之间定位波谷导致的。图7为从图2脉搏信号对应完整信号中提取的PRV信号，可以看出所提方法可以准确提取PRV信号。

图6 脉搏波波谷定位结果Fig.6 Pulse wave and its trough locations

图7 PRV信号提取结果Fig.7 A PRV signal extracted from pulse signal

将所提算法用于31组实验数据的PRV信号提取，检出数、漏检数、误检数和准确率如表1所示。

表1 PRV信号提取准确率Tab.1 Accuracy of PRV signal extraction

可以看出，31组数据中，25组PRV信号的提取准确率为100%，其余各组的准确率为97%以上，其中有4组的准确率在99.5%以上。只有1组的准确率为97.35%，这是由于原始脉搏信号中运动伪迹造成的干扰段导致。总体而言，所提算法的准确率很高，这是因为通过VMD分解后，原始脉搏信号的高频噪声和重搏波、潮波等干扰特征分解到其他IMF分量中，对主波波谷定位的干扰降低，从而得到较高的准确率。

4 结束语

PRV信号准确提取十分重要，本研究提出一种基于VMD分解的脉率变异性提取算法，可准确地从原始脉搏信号中提取PRV信号。具有抗干扰性强不需要设置阈值等优点，可用于日常工程生活环境下PRV信号的准确提取。然而，VMD分解速度较慢，耗用系统资源较多，严重限制了所提算法的应用场景。因此，未来工作将对VMD分解算法进行改进，提高其运算速度，进而拓阔所提算法的应用领域。