基于EEMD的脉率变异性时频域研究

0 引言

心率变异性(Heart Rate Variability，HRV)是指逐次心搏间期之间几十毫秒的微小变异[1]，主要是由于心脏窦房结自律活动通过交感和迷走神经活动等因素的调节作用而产生。通过对HRV的生理和病理意义的研究，表明HRV信号中包含着大量的有关心血管调节的重要信息，目前HRV已经成为自主神经系统的活动情况的重要检测手段，且可为各种心血管疾病的预防和诊断提供了很好的参数依据[2-3]。近年来，许多研究已经证明，在一定程度上脉率变异性(Pulse Rate Variability，PRV)和HRV一样，能够准确地反映出自主神经系统的活动情况[4-5]。

传统的HRV(PRV)分析方法主要有：时域分析、频域分析、非线性分析[6-7]。其中功率谱密度(Power Spectral Density，PSD)估计是一种常用的HRV(PRV)频域分析方法，它是将一定时间内记录的相邻RR间期(RR Interval，RRI)序列经过快速傅里叶变换(Fast Fourier Transform，FFT)或自回归(Auto-regressive，AR)，得到以序列的频率为横坐标，功率谱密度为纵坐标的HRV(PRV)谱。通过研究PSD中功率随频率的变化关系，可以了解自主神经对心率的调节作用，能够反映交感神经和迷走神经的活动水平，对病人的监护和治疗具有一定的指导意义[8]。典型的HRV(PRV)的PSD曲线可以分三个频段，其中极低频（VLF）频率范围在（0.003～0.04） Hz，其功率主要反映人体心率变化受温度及肾血管紧张素系统调节作用；低频段（LF）频率范围在（0.04～0.15） Hz，其功率反映交感和迷走神经的双重调节，交感神经占主导地位；高频段（HF）频率范围在（0.15～0.4）Hz，其功率反映迷走神经的活动，与人体的呼吸调节有关。总频率功率（TP）是自主神经系统对心血管系统影响的评判标准；LF/HF是反映交感和迷走神经活动均衡性的指标[9-11]。FFT频谱估计虽然算法简单，但分析HRV(PRV)时需对信号进行周期延拓，短程HRV(PRV)分析分辨率低，且会出现能量泄露现象。AR频谱估计法所需HRV(PRV)信号较短，短程HRV(PRV)分析分辨率较高，但确定阶数较为困难。

近年来，随着信号处理技术的发展，时频分析方法（如，短时傅里叶变换、小波变换等）由于其具有能够反映信号的时频特性的优点，因此被广泛应用到各个领域。但此类信号处理方法一般都要求被处理的信号是平稳的，而实际上很多人体生理信号包括HRV(PRV)信号都是非线性、非平稳的信号，如果把此类信号当做平稳信号来处理，必有一定的局限性。而集合经验模态分解(Ensemble Empirical Mode Decomposition，EEMD)能够根据信号自身的时间尺度特征将信号自适应地分解成高频到低频不同时间尺度分量，是一种非常适合非平稳、非线性信号的时频域分析方法[12]。由此本文提出一种先对HRV(PRV)信号进行EEMD多层分解，计算出各个本征模态函数(Intrinsic Mode Function，IMF)的能量，再根据分解的IMF分量的频谱特征，对IMF分量的能量进行重构，得到能够表示HRV(PRV)的高频段、低频段、极低频段能量的指标。

1 方法

1.1 脉搏波信号的采集系统

本文的脉搏波信号采集系统是以光电容积脉搏波描记法(Photoplethysmograph，PPG)[13-14]为基本原理，主要由反射式脉搏波光电传感器、信号调理模块、模数转换模块、微控制器、信号隔离模块、串口转USB口模块等部分构成，如图1所示。通过反射式脉搏波光电传感器采集的脉搏波信号，先由信号调理模块进行滤波和放大，经模数转换得到数字信号送微控制器处理后，再通过USB接口送上位机做进一步的滤波处理、数据存储和波形显示等。本脉搏波信号采集系统各个模块的直流供电均直接取自于上位机的USB接口的+5 V电源，省去了电源模块。同时为了保证使用者的安全，在采集系统中加入了信号隔离模块。

图1 脉搏波信号采集系统框图Fig.1 Block diagram of pulse wave signal acquisition system

图2 脉搏波峰值序列Fig.2 Pulse wave peak sequence

根据已确定的脉搏波峰值位置，将相邻两个波峰位置对应的采样点作差，得出相邻波峰之间的采样点数N，再用采样间隔乘以N，即t=T×N，t就是所需的PPI，以PPI为纵坐标得到的序列图就是PRV信号，图3是一段5 min脉搏波提取的PRV信号序列图。很显然PRV信号在时域中是不等间隔的，如采用传统的FFT谱估计法或AR谱估计法分析PRV，还需事先对PRV信号序列进行三次样条插值和重采样操作，这种操作会对分析结果造成一定的影响。由于本文采用EEMD方法适用于不等间隔序列的分解，因此就省

图3 PPI序列图Fig.3 Diagram of PPI sequen ce

1.2 PRV信号的提取

本文采用基于可变阈值的方法提取脉搏波信号的峰值。由于人体在静息状态下的脉搏波信号的峰值波动的范围一般不会超过峰谷值高度的0.4倍，本文先将采集的5 min长度的脉搏波信号的最大峰谷值的0.4倍作为初始阈值。只要脉搏波信号中的某些点同时满足下面两个条件：① 该点的幅值加上阈值大于最大的峰值；② 该点在其邻域内（本文脉搏波信号的采样频率为500 Hz，邻域长度定义为200个数据点）是最大值，即可认为这些点是脉搏波信号的峰值。但由于个体的差异和各种干扰的影响，有可能出现极少数脉搏波信号的峰值波动会大于初始阈值，从而造成峰值漏检的情况。本文根据计算出的脉搏波的峰峰值间期(Peak-Peak Interval , PPI)大小的情况来改变初始阈值，当PPI中有出现大于平均脉动周期的3/2时，说明出现峰值漏检情况，程序将以0.1倍初始阈值为步长自动增加阈值，重新提取峰值，直到没有漏检为止。采用本文方法可以做到无漏检、误检的情况。检测的脉搏波峰值如图2所示。去了这一步。

1.3 PRV信号的EEMD分解及能量重构

经验模态分解(Empirical Mode Decomposition，EMD)作为一种很好的时频分析工具，很适用对非线性、非平稳信号的分析，但也存在一些缺点，例如信号经过分解后会产生模态混叠现象。为了消除EMD分解带来的模态混叠问题，Huang提出了EEMD方法，当要分解的信号处在白噪声频谱均匀分布的时频空间上时，不同时间尺度信号会自动分布到合适的参考尺度上，利用零均值噪声的特性，经过多次平均运算后，噪声会自动相互抵消，集成多次分解的均值即可作为最终分解结果。EEMD的实质是对时域的信号序列进行平稳化处理，其结果是信号中不同尺度的波动被逐级分解出来，产生多层尺度不同的数据序列即本征模态函数IMF，然后可进一步对分解后的各个IMF分量进行频谱分析。

EEMD 分解的步骤如下：

(1) 向待分解的信号x(t)加入均值为零、标准差为常数的白噪声n(t)，得到信号X(t)，即：

对一名受试者在坐立静息状态下采集5 min脉搏波信号，利用形态学滤波器和SG平滑滤波器分别去除脉搏波信号中的基线漂移和高频干扰信号。采用上述PRV信号提取方法从滤波后的的脉搏波信号提取PRV信号，再对PRV信号进行EEMD 7层分解，各个IMF分量及其对应的频谱图，如图4所示，从IMF频谱图可以看出PRV信号中的HF（0.15 Hz～0.4 Hz）分量主要分布在IMF1中，LF(0.04 Hz～0.15 Hz)分量主要集中在IMF2和IMF3中，VLF(0.003 Hz～0.04 Hz)分量主要分布在IMF4～IMF7中，Cr表示PRV信号的趋势。

由数字信号序列能量定义可知，对各个IMF分量平方并在时域内（0～300 s）做积分运算，即可得到各个IMF分量的能量，再根据各个IMF分量的能量分布情况进行各个频段的能量合成，可以分别计算出PRV信号中的HF（0.15 Hz～0.4 Hz）能量EHF、LF（0.04 Hz～0.15 Hz）能量ELF、VLF（0.003 Hz～0.04 Hz）能量EVLF等，如公式(7)～(13)所示。

图4 IMF分量及其频谱Fig.4 IMF components and its spectrum

2 结果

随机选取20名健康的受试者在坐立静息状态下采集5 min脉搏波信号，分别采用EEMD时频方法和AR功率谱估计方法测出各个频域参数，并对两种方法测得的参数通过SPSS 19.0统计软件进行对应相关系数计算，结果如表1所示。

可见，基于EEMD的脉率变异性时频域分析方法和基于AR功率谱估计方法测得的对应各项指标相关系数均大于0.96，具有高度相关性。

表1 两种方法的对比Tab.1 Comparison of the two methods

3 结论

本文根据EEMD适用于分析非线性、非平稳信号的特点，提出一种基于EEMD的计算HRV(PRV)频域参数的方法，即先对HRV(PRV)信号进行EEMD分解成各个模态函数，再由各个模态函数频谱的特点进行能量重构，即可得到HRV(PRV)的几个频域参数。本文方法省去了传统HRV(PRV)频域分析方法所必须的插值和重采样环节。通过和AR频谱估计方法相对比，结果表明本文方法能够有效、准确地反映人体交感神经和迷走神经活动，为临床HRV(PRV)频域分析提供一种新的方法。在评价心血管疾病的预后，预测急性心血管事件和评估人体精神压力方面具有重要意义。

[1] Jia R T, Liu B. Human daily activity recognition by fusing accelerometer and multi-lead ECG data[C]//2013 IEEE International Conference on Signal Processing, Communication and Computing (ICSPCC). Kunming: 2013: 1-4.

[2] Scheff J, Corbett S, Calvabo S, et al. The relationship between autonomic function and heart variability in human endotoxemia [J]. J Crit Care, 2013, 28(6): e32.

[3] Guo Y, Palmer J L, Starasser F, et al. Heart rate variability as ameasure of autonomic dysfunction in men with advanced Cancer[J]. Eur J Cancer Care (engl), 2013, 22(5): 612-615.

[4] Lu S, Zhao H, Ju K. Can photoplethysmography variability serve as an alternative approach to abtain heart rate variability intormation[J]. J Clin Monit coput, 2008, 22(1): 23-29.

[5] 石波, 陈法胜, 陈建方, 等. 脉率变异性分析在穿戴式智能设备中的应用[J]. 中国医疗器械杂志, 2015, 39(2): 95-97.

[6] Chourasia V S, Mittra A K. Wavelet-based denoising of fetal phonocardiographic signals[J]. Int J Med Eng Inform, 2010, 2(2): 139-150.

[7] 黄晓林. 心率变异性的分析方法研究[D]. 南京: 南京大学，2009.

[8] Brown T E, Beightol L A, Koh J, et al. Important influence of respiration on human R-R interval power spectral is largely ignored[J]. Appl Physiol, 1993, 75(5): 2310-2317.

[9] 董红生. 心电波形检测与心率变异性分析方法研究[D]. 兰州: 兰州理工大学, 2012.

[10] Lucini D, Norbiato G, Clerici M, et al. Hemodynamic and autonomic adustments to real Life stress conditions in humans[J]. Hypertension, 2002, 39(1): 184-188.

[11] Fujiwara K, Miyajima M, Yamakawa T, et al. Epileptic seizure prediction based on multivariate statistical process control of heart rate variability features[J]. IEEE Trans Biomed Eng, 2016, 63(6): 1321-1332.

[12] 刘俏俏, 谭志向, 张懿, 等. 基于经验模态分解的胎心音信号除噪方法[J]. 生物医学工程学杂志, 2015, 32(4): 740-745.

[13] John A C, Damianou D. The wavelength dependence of the photoplethysmogram and its implication to pulse oximetry [C]// 14th Annual International Conference of the IEEE Engineering in Medicine and Biology Society. PARIS: 1992: 1173-1181.

[14] 罗志昌, 张松, 杨益民. 脉搏波的工程分析与临床应用[M]. 北京:科学出版社, 2006.

Time-frequency Research of Pulse Rate Variability Based on the EEMD

【 Writers 】ZHANG Genxuan1, SHI Bo1, ZHANG Sai1, TSAU Young2
1 Department of Medical Imaging, Bengbu Medical College, Bengbu, 233030 2 Dimetek Digital Medical Technologies Ltd, Shenzhen, 518067

pulse rate variability, heart rate variability, ensemble empirical mode decomposition

R318.6

10.3969/j.issn.1671-7104.2017.04.005

1671-7104(2017)04-0251-04

2016-12-07

安徽省高等学校自然科学研究一般项目(KJ2015B062by)

石波，E-mail: shibohome@qq.com

【 Abstract 】In this paper, based on the feature of ensemble empirical mode decomposition(EEMD) which is very suitable for analyzing non-stationary and nonlinear signals, the pulse rate variability (PRV) signal is decomposed into multiple intrinsic mode function (IMF) components by EEMD, to calculate the energy of each IMF component respectively, and high-frequency energy(EHF), low-frequency energy(ELF) and very low-frequency energy(EVLF) of PRV were reconstructed from the energy of the IMF components according to the characteristics of each IMF component spectrum. This method is compared with the traditional AR power spectrum estimation method through the experiment. The results show that the correlation coefficient of each frequency domain parameter corresponding to PRV obtained by the two methods is higher than 0.96, which indicates that the method can truly reflect the sympathetic nerve and vagus nerve activity.