0 概述

脉搏信号是人体重要的生理信号, 对其进行分析可以掌握人体各器官的综合信息^[1]。利用现代电子技术采集人体的脉搏信号, 并采用信号分析方法从中提取诸如幅度、频率、形状变化等有效信息是实现脉搏可靠诊断的关键^[2]。

针对脉搏信号的特征, 常用的分析方法主要有时域分析法、频域分析法、时频特征分析法和非线性分析法等^[3]。时域分析法主要提取脉搏波的主波、潮波、重搏波等有一定生理意义的点, 如文献[4]运用时域分析方法对冠心病患者的脉搏信号进行分析。但时域分析法提取的特征参数较复杂, 且信号中幅值较小的病理信息无法被提取出来。频域分析法在一定程度上弥补了时域分析中的缺陷, 如文献[5]利用小波变换提取脉搏波在多分辨率层次上的小波熵作为频域特征。频域分析法虽然能较好地反映脉搏信号在整个频率段上的特性, 但是信号的局部特性却丢失了。时频分析法能够保留脉搏信号的局部特性, 得到信号在某一时刻的瞬时幅度和频率特征^[6], 常见的方法有短时傅里叶变换^[7]和小波分析^[8]。随着非线性分析方法的发展, 信息熵、近似熵等方法也被应用到脉搏信号的分析中, 如文献[9]用近似熵能量比来表征30 min脉搏信号波形的变异性, 从而区分冠状动脉硬化患者和健康人。上述方法虽然在脉搏信号的分析中得到了较好的应用, 但对于脉搏信号的非线性、非平稳特性仍需要进行深入的研究。

希尔伯特-黄变换(Hilbert-Huang Transform, HHT)^[10]是一种新型自适应处理非线性、非平稳信号的时频分析方法。本文对传统的HHT脉搏信号分析方法进行改进, 结合时变滤波(Time Varying Filtering, TVF)的经验模态分解(Empirical Mode Decomposition, EMD)^[11]法和相对系数法自适应提取有效的脉搏信息。

1 改进的HHT分析方法

HHT主要包括两部分:一部分是基于瞬时频率的EMD方法, 其能够将任一复杂的信号自适应地分解成一系列本征模态函数(Intrinsic Mode Function, IMF); 另一部分是希尔伯特变换(Hilbert Transform, HT), 通过对各IMF分量运用HT可以得到各分量的瞬时幅值和频率信息^[12]。传统的HHT由于模态混叠问题会导致部分IMF分量及其频谱分布失去实际的物理意义^[13], 因此, 本文对HHT进行改进。

1.1 经验模态分解

EMD以信号在时频上的局部特性为基础, 对输入信号进行“筛选”, 其“筛选”的目的是不断地将低频分量从信号中分离, 最终得到关于时间轴对称的IMF。

对于输入信号x(t), 将其经过EMD分解后得到的一组IMF和一个残余分量分别用imf_m(t)和r(t)来表示:

$ x\left( t \right) = \sum\limits_{i = 1}^m {im{f_i}\left( t \right)} + r\left( t \right) $

(1)

1.2 TVF-EMD

由于群延迟和通带波动, 传统的线性滤波器并不适合处理非线性信号。B样条是用于信号分析的重要工具, 等间隔B样条拟合具有低通滤波器性质, 其截止频率由节点间距决定; 非等间隔B样条拟合具有时变滤波器性质, 其局部截止频率由局部节点间隔决定^[14]。

TVF-EMD利用B样条最小二乘法构造时变滤波器, 根据瞬时幅度和瞬时频率自适应地设计局部截止频率。通过连续应用时变滤波器, 最终得到包含较低频率的局部窄带信号, 该信号具有与IMF相似的特性。

一个包含多分量的信号可以用2个信号的组合来表示, 即:

$ z\left( t \right) = A\left( t \right){{\rm{e}}^{{\rm{j}}\varphi \left( t \right)}} = {a_1}\left( t \right){{\rm{e}}^{{\rm{j}}{\varphi _1}\left( t \right)}} + {a_2}\left( t \right){{\rm{e}}^{{\rm{j}}{\varphi _2}\left( t \right)}} $

(2)

为了得到时变滤波器的局部截止频率, 需要找出φ₁(t)和φ₂(t)的平均频率。

对于含有2个窄带分量的信号, 有:

$ {A^2}\left( t \right) = a_1^2\left( t \right) + a_2^2\left( t \right) + 2{a_1}\left( t \right){a_2}\left( t \right)\cos \left[ {{\varphi _1}\left( t \right) - {\varphi _2}\left( t \right)} \right] $

(3)

$ \begin{array}{l} \varphi '\left( t \right) = \frac{1}{{{A^2}\left( t \right)}}\left\{ {{{\varphi '}_1}\left( t \right)\left\{ {a_1^2\left( t \right) + {a_1}\left( t \right){a_2}\left( t \right)\cos \left[ {{\varphi _1}\left( t \right) - } \right.} \right.} \right.\\ \;\;\;\;\;\;\;\;\;\;\;\left. {\left. {{\varphi _2}\left( t \right)} \right]} \right\} + {{\varphi '}_2}\left( t \right)\left\{ {a_2^2\left( t \right) + {a_1}\left( t \right){a_2}\left( t \right)\cos \left[ {{\varphi _1}\left( t \right) - } \right.} \right.\\ \;\;\;\;\;\;\;\;\;\;\;\left. {\left. {{\varphi _2}\left( t \right)} \right]} \right\} + \frac{1}{{{A^2}\left( t \right)}}\left\{ {{{a'}_1}\left( t \right){a_2}\left( t \right)\sin \left[ {{\varphi _1}\left( t \right) - {\varphi _2}\left( t \right)} \right] - } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\left. {{{a'}_2}\left( t \right){a_1}\left( t \right)\sin \left[ {{\varphi _1}\left( t \right) - {\varphi _2}\left( t \right)} \right]} \right\} \end{array} $

(4)

通过对A(t)的局部最小值和最大值应用插值, 可以得到2个信号分量的瞬时频率φ′₁(t)和φ′₂(t), 进而求得2个分量的平均频率:

$ {{\varphi '}_{{\rm{bis}}}}\left( t \right) = \frac{{{{\varphi '}_1}\left( t \right) + {{\varphi '}_2}\left( t \right)}}{2} $

(5)

n阶B样条空间中的任意信号可表示为:

$ g_m^n\left( t \right) = \sum\limits_{k = - \infty }^\infty {c\left( k \right){\beta ^n}\left( {\frac{t}{m} - k} \right)} $

(6)

其中, βⁿ(t)为B样条基函数, m为节点间隔, c(k)为控制系数。拟合的结果由n、m及c(k)来确定。

对于给定信号x(t), 其控制系数c(k)利用最小二乘法得到:

$ \min \varepsilon _m^2 = \min \sum\limits_{t = - \infty }^{ + \infty } {{{\left[ {x\left( t \right) - g_m^n\left( t \right)} \right]}^2}} $

(7)

信号的节点由h(t)=cos[∫φ′_bis(t)dt]计算得到, 且局部截止频率将根据节点分布逐渐收敛。B样条拟合的频率响应特性如图 1所示。从图 1(a)可以看出, 阶数n越大, B样条拟合的结果越接近理想的低通滤波器。图 1(b)的结果表明该滤波器的截止频率约为1/(2m)。因此, 时变滤波器的局部截止频率与平均频率是一致的, 将B样条最小二乘的拟合结果记为m(t)。

EMD的停止条件只是通过判断标准偏差指标, 并没有直接对IMF的2个条件进行检查。TVF-EMD将定义更严格的停止准则。对于式(3)中的2个信号分量, 其瞬时频率的加权平均值和瞬时带宽分别由式(8)和式(9)得到。

$ {\varphi _{{\rm{avg}}}}\left( t \right) = \frac{{a_1^2\left( t \right){{\varphi '}_1}\left( t \right) + a_2^2\left( t \right){{\varphi '}_2}\left( t \right)}}{{a_1^2\left( t \right) + a_2^2\left( t \right)}} $

(8)

$ B\left( t \right) = \sqrt {\frac{{a{'}_1^2\left( t \right) + a{'}_2^2\left( t \right)}}{{a_1^2\left( t \right) + a_2^2\left( t \right)}} + \frac{{a_1^2\left( t \right)a_2^2\left( t \right){{\left[ {{{\varphi {'}}_1}\left( t \right) - {{\varphi {'}}_2}\left( t \right)} \right]}^2}}}{{{{\left[ {a_1^2\left( t \right) + a_2^2\left( t \right)} \right]}^2}}}} $

(9)

一般情况下, 当信号的瞬时带宽B(t)足够小时, 该信号被认为是局部窄带信号。由于B(t)是一个绝对量, 因此本文采用一个相对标准来衡量它偏离φ_avg(t)的程度。

$ \theta \left( t \right) = \frac{{B\left( t \right)}}{{{\varphi _{{\rm{avg}}}}\left( t \right)}} $

(10)

当θ(t)变小时, 信号的瞬时带宽也跟着变小。设定一个阈值ξ, 当θ(t)≤ξ时, 即可认为该信号为局部窄带信号, 即IMF分量。对于输入信号x(t), 如果θ(t)≤ξ(本文ξ=0.1), 则x(t)被视为一个IMF分量; 否则, 令x(t)=x(t)-m(t), 并重复分解过程。

1.3 有效分量的提取

相关系数是回归分析方法的基本统计量, 它主要用来表示2个变量的线性相关性^[15]。由于TVF-EMD分解出的IMF个数会受到原始信号点数的影响, 为了避免TVF-EMD过分解的影响, 本文引入了相关系数法。以分解出的各阶IMF与原始信号之间的相关系数通过阈值标准来判断IMF是真实的信号分量或是无实际意义的虚假分量。

本文方法对所有的IMF和原始信号进行归一化处理, 从而避免剔除IMF幅值较小但又属于真实信号分量的情况。各阶IMF与原始信号的归一化相关系数为:

$ {R_i} = \frac{{\sum\limits_{n = 1}^N {\left( {{x_n} - \bar x} \right)\left( {im{f_{in}} - \overline {im{f_i}} } \right)} }}{{\sqrt {\sum\limits_{n = 1}^N {{{\left( {{x_n} - \bar x} \right)}^2}} } \sqrt {\sum\limits_{n = 1}^N {{{\left( {im{f_{in}} - \overline {im{f_i}} } \right)}^2}} } }} $

(11)

将相关系数的标准差作为判定的阈值th, 即:

$ th = {\left( {\frac{1}{{m - 1}}\sum\limits_{i = 1}^m {{{\left( {{R_i} - \bar R} \right)}^2}} } \right)^{\frac{1}{2}}} $

(12)

当R_i>th时, 对应的IMF为有效分量; 反之, 该IMF为无效分量。

1.4 改进的HHT流程

本文改进的HHT方法首先对输入信号进行TVF-EMD分解得到一组IMF, 然后利用相关系数法提取有效的IMF信号分量, 最后对提取的IMF进行HT得到信号的Hilbert谱和边际谱, 其流程如图 2所示。

2 实验与结果分析 2.1 脉搏信号的分解

本文实验使用压力脉搏传感器来采集正常人体桡动脉处的脉搏信号, 采样点数为2 000个。采集的脉搏信号如图 3所示。

将脉搏信号分别进行EMD和TVF-EMD分解, 结果如图 4和图 5所示。由于篇幅所致, 图 5只列出了主要的分解结果。从图 4和图 5可知, EMD方法得到了7阶IMF分量, 而TVF-EMD方法得到了18阶IMF分量。TVF-EMD方法对脉搏信号具有更好的分解性能, 有效解决了模态混叠问题。在TVF-EMD分解中, IMF1~IMF7对应脉搏信号中的高频噪声干扰, 这些分量的幅值小、频带宽, 没有脉搏信号的实际特征。IMF8~IMF17对应一定周期性的脉搏信号, 体现出心脏收缩和舒张的主要特征, 并且这些分量还包含脉搏信号中主波、潮波和重搏波的主要特征, 是脉搏信号特征提取的关键部分。TVF-EMD能够将脉搏信号更准确地分解成具有特定含义的IMF分量。

2.2 脉搏信号有效IMF的提取

应用相关系数法得到18阶IMF分量与原始脉搏信号之间的归一化相关系数如表 1所示, 通过式(12)计算判定阈值th为0.062 4。

图 6为各阶IMF对应的相关系数与标准差的关系图。根据相关系数大于标准差的准则, 提取出有效的IMF分量为IMF10、IMF11、IMF12、IMF13、IMF16和IMF17。

2.3 脉搏信号的Hilbert谱

将分解出的信号分量进行HT可以得到具有时频特征的Hilbert谱。脉搏信号经传统HHT方法分解后和改进HHT方法分解后得到的Hilbert谱分别如图 7和图 8所示。

对比图 7和图 8可以看到, 传统HHT方法得到的Hilbert谱出现了较为严重的失真, 谱线有严重的震荡和断裂现象, 且无法得到频率较高的分量; 而改进HHT方法分解后得到的Hilbert谱明显含有噪声产生的高频信号的谱线, 即信号中的干扰成分被更好地分离。由图 8可知, 通过设定阈值对有效分量提取后即可将信号中的干扰成分去除。在0 Hz~20 Hz的范围内有效模态分量的谱线更为完整, 没有出现明显的失真, 即对应着脉搏信号3个主要波段。因此, 通过脉搏信号的Hilbert谱可以观察到脉搏信号的本质特征, 具有临床意义。

2.4 脉搏信号的Hilbert边际谱

Hilbert边际谱能够反映出脉搏信号的幅值在整个频域上的变化情况。脉搏信号在传统HHT方法和改进的HHT方法下得到的Hilbert边际谱分别如图 9和图 10所示。

从图 9和图 10可以看到, 脉搏信号的频率主要分布在0 Hz~20 Hz的范围内, 且在0 Hz~10 Hz的频率范围内信号分量最为密集。虽然2种方法均体现了脉搏信号的本质特征, 但是, 利用改进的HHT得到的Hilbert边际谱不仅有效克服了高频噪声对信号的影响, 而且信号主成分的图谱特征更加明显。由图 10可知, 脉搏波的主波、潮波和重搏波的频率范围分别在2 Hz、3 Hz和5 Hz附近。通过分析脉搏信号Hilbert边际谱的频率域是否出现异常的频率分量, 可以对人体的健康状况做出初步诊断。

3 结束语

脉搏信号含有大量的人体生理信息, 其特征参数与人体心血管系统的健康状况密切相关。针对脉搏信号的非平稳特性以及传统HHT中EMD存在的模态混叠问题, 本文对HHT脉搏信号分析方法进行改进。实验结果表明, 改进方法可大幅提升信号分解性能, 并且有效抑制模态混叠现象。通过改进方法得到的Hilbert谱和Hilbert边际谱能够准确反映出脉搏信号的时频特征, 这些特征为脉搏信号的模式识别提供了可靠依据。