基于经验模态分解和形态滤波的血压测量研究

甘永进，莫 沛，陈时东，郑金存，*，蒋曲博

(1. 玉林师范学院 电子与通信工程学院，广西 玉林 537000；2. 广西科技大学鹿山学院 电气与计算机工程系，广西 柳州 545000；3. 桂林电子科技大学 电子工程与自动化学院，广西 桂林 541004)

0 引言

目前，示波型电子血压计是市场的主流[1]。然而因心率异常导致振荡波不规则以及血管形变恢复需要一定时间，示波型电子血压计无法实现连续测量每搏血压的变化[2]，而连续血压测量在家庭保健和临床诊断中意义重大。通过光电容积脉搏波(Photoplethysmography,PPG)完成对每个心动周期血压值的检测，为实现血压的无创连续测量提供一定的参考价值。但是PPG信号容易受到采集系统带来的白噪声、工频干扰[3]以及由肌电干扰等因素引起的基线漂移等噪声的干扰，对后续血压测量的精度带来影响。

为实现由PPG信号达到连续精准地测量血压的目的，本文通过TI系列AFE44X0芯片结合反射式传感器DCM03设计光电容积脉搏波采集系统，设计经验模态分解(Empirical Mode Decomposition,EMD)结合形态学滤波的方法对原始PPG信号进行自适应滤波以及基线去除工作，对滤波后的PPG信号进行二次微分得到加速脉搏波，采用Adaboost算法对加速脉搏波进行分类，对分类后的加速脉搏波进行PPG信号传导时间(Pulse Wave Transit Time， PWTT)的提取，最后建立不同的血压测量模型。

1 光电容积脉搏波去噪算法

1.1 EMD去噪

经验模态分解由数据时间尺度特征来进行信号分解以得到一系列的本征模函数(Intrinsic Mode Function，IMF)，其表征信号的不同振荡模式，包含不同时间尺度[4-5]信号局部特征。其中，各个IMF中零点和极点个数相等或相差1，且局部极大值和极小值确定的信号包络线的均值为0[6]，使得定义瞬时频率的同时可避免由信号不对称引起的瞬时频率波动[7]。

本文采用光电反射式传感器DCM03结合血氧模拟前端AFE44X0系列芯片进行反射式光电容积脉搏波的采集。反射式传感器DCM03的光发射器发出光线照射到指尖后发生漫反射，DCM03的光接收器接收指尖血液的光感应信号，实现光电信号转换，再由AFE44X0进行I-V转换、放大滤波及A/D转换等处理得到PPG数字信号并输出至MCU[8]。将PPG数字信号进行EMD分解， 噪声分布在各个IMF分量上，故将所有的IMF分量进行软阈值去噪，可有效地对原始PPG信号的噪声进行抑制，其步骤如下:

①对原始PPG信号进行EMD分解，得到各个IMF分量，如图1所示；

②对各IMF分量采用Stein无偏似然估计准则进行阈值估计；

③使用wthresh进行去噪，得到去噪后各IMF分量如图2所示；

④将去噪后的各IMF分量进行重构得到去噪后的PPG信号。

图1 容积脉搏波的EMD分解图
Fig.1 EMD decomposition diagram of PPG

图2 去噪后的IMF分量
Fig.2 IMF component after denoising

1.2 基于形态学的基线漂移滤除

设采样得到的原始一维脉搏信号为f(n)，定义域为F={0,1,…,N-1}，结构元素为g(n)，定义域为G={0,1,…,M-1}，且M

(1)

(2)

据此，由腐蚀与膨胀运算组成开、闭运算为

(f∘g)(n)=(f⊖g⊕g)(n),

(3)

(f·g)(n)=(f⊕g⊖g)(n)。

(4)

一般，闭运算用于填充细小空洞，实现平滑或抑制信号波谷噪声；开运算可用于断开窄小的连接，消除微小的尖刺，滤除信号峰值噪声，平滑信号边界轮廓。常采用形态开、闭的级联形式对信号进行处理。传统形态的开-闭和闭-开运算以不同顺序级联开闭运算，定义为

FOC(f,g)=f∘g·g,

(5)

FCO(f,g)=f·g∘g。

(6)

因形态闭运算的反扩展性和形态开运算的扩展性，式(5)和式(6)定义的两种传统的滤波器都存在统计偏移现象，即对于开-闭滤波器而言，最后的输出幅度偏小；但对于闭-开滤波器而言，最后的输出幅度偏大，在一般情况下，单独使用得到的滤波效果都不是最佳的。

欲有效地抑制采集到的脉搏信号中的不同噪声，减小最终输出的单向偏移，由两种滤波器的平均组合形式[9]，将形态开闭-闭开滤波器定义为

FOC-CO(f,g)=(f∘g·g+f·g∘g)/2。

(7)

欲将统计偏移现象更深一层地减小，根据不同尺寸的结构元素，定义广义开-闭和闭-开滤波器为

FOC(f,g1,g2)=f∘g1·g2,

(8)

FCO(f,g1,g2)=f·g1∘g2。

(9)

相比于传统的形态滤波器，广义形态滤波器能够对信号中的各种噪声进行有效的抑制[10]。为校正基线漂移[11]，本文通过尺寸不同的结构元素的开-闭和闭-开运算组合进行处理。首先，设原始脉搏血氧信号为x(n)，对受基线漂移干扰的x(n)进行广义形态闭-开运算处理；然后，将x(n)进行广义形态开-闭运算处理；之后，再把以上两个步骤的结果进行求和平均，得到基线分量；最后，将x(n)与基线作差，得到校正基线后的信号y(n)，形态学去除基线漂移[12]框图如图3所示。 由表1小波及本文算法对信号去噪后的信噪比对比知，该算法去除基线漂移效果较佳。经过EMD分解以及形态学滤波后的波形和频谱如图4、图5所示，由图4和图5可见，信号中的高频噪声以及基线漂移得到一定程度的抑制，说明了设计的算法对光电容积脉搏波消噪效果的可行性。

图3 去除基线漂移框图
Fig.3 Block diagram of baseline wander removal

表1 几种不同去基线算法的信噪比
Tab.1 SNR of different algorithms

序号方法SNR 1db2小波38.607 22db4小波38.216 33db5小波38.987 34db6小波38.527 55sym4小波38.587 76sym8小波38.497 37形态学滤波41.027 6

图4 各个阶段波形的频谱图
Fig.4 Spectrogram of each stage

图5 各个阶段波形
Fig.5 Wave form of each stage

2 血压计算方法

2.1 脉搏波传导时间计算

PPG信号经EMD和形态学滤波后，再提取PPG信号的传导时间来建立PPG信号传导时间与血压的数学关系，以此建立血压计算模型。其中，PPG信号及相应的加速脉搏波波形如6所示。

图6 PPG信号及其加速脉搏波波形
Fig.6 PPG and accelerating pulse wave

图6中，由血液微循环机理知加速脉搏波A-C段的时间tAC可以反映血液由心脏传播到指尖末端毛细血管并反射汇合的时间，而这个时间可准确的表示脉搏波传导时间。脉搏波传导时间的提取的关键是找到特征点A和C的位置。其中，特征点A是波形的极大值点，准确地识别该位置较为容易。故特征点C的位置识别成为脉搏波传导时间提取的关键。特征点C位置出现有两种情况：一种是特征点C位于B和Z之间，如图7(a)所示；另一种是特征点C位于Z和M之间，如图7(b)所示。

图7 特征点C位于不同位置的加速脉搏波
Fig.7 Acceleration pulse wave with different feature points

本文首先对经数字信号处理后的PPG信号的加速脉搏波进行分类[13]，再提取两种不同类型的PPG信号的传导时间，进而分别建立两种PPG信号(特征点C位于B和Z之间或Z和M之间)的血压计算模型。

2.1.1 加速脉搏波波形识别

为识别出如图7中的两种不同的加速脉搏波，本文利用迭代算法Adaboost进行实验。 Adaboost在训练过程中自适应调整训练样本权重，强化学习难以区分的训练样本，由加权投票机制联合弱分类器，使分类效果好的弱分类器的权重得到提高，而分类效果差的弱分类器的权重降低，形成高性能的强分类器。算法描述如下：

①输入训练集D：

D={(x1,y1),(x2,y2),…,(xn,yn)},

其中，xi∈X,yi∈Y,X为训练样本集，Y为分类类别标志，Y={-1,+1}。

②初始化权值：D1(x)=1/n。

③设置训练轮数T：fort=1,2，…，T。

④将弱学习算法在权值Dt下训练，得到预测函数ht：X～{-1,+1}。

⑦训练完毕，输出函数：

为验证Adaboost算法识别不同脉搏波的可行性，本文由预测灵敏度SE、阳性率TPR以及准确率AC这3个指标来评价该算法区分两种脉搏波的性能。其中，

(10)

(11)

(12)

其中，TP和TN分别表示正确分类了第一类和第二类脉搏事件，FP和FN分别表示错误分类了第一类和第二类脉搏事件。

对250个波形信号以Adaboost算法进行识别，其中第一类(特征点C位于B和Z之间)130个，第二类(特征点C位于Z和M之间)120个，得到测量数据如表2所示。由表1知，以特征点A-Z的时间段tAZ作为分类依据，采用Adaboost算法可以准确的区分出两种不同的加速脉搏波。

表2 分类评估参数
Tab.2 Evaluation parameters

TPTNFPFN灵敏度阳性率准确率11910411160.8810.9150.892

2.1.2 特征点提取

特征点A是加速脉搏波的波峰，在波形任意位置开始，在一个周期范围内判断出幅值最大的位置即为特征点A的位置。

特征点C的位置提取分两种情况：①当C介于B与Z之间时，如B与Z间存在极大值，则极大值坐标即为C的位置；如没有极大值，则以B与Z之间斜率最小的点的坐标作为C的位置；②当C介于Z与M之间时，若B、M两点所在的直线斜率为k，以Z、M间所有点的斜率与k差值最小时的坐标作为C的位置。准确判断特征点A和C的位置后，A与C之间的采样点数与采样时间的乘积即脉搏波传导时间(PWTT)。

2.2 血压模型的建立

准确识别出两种不同的脉搏波并获取脉搏波传导时间后，由脉搏波传导时间与血压间的线性关系，进行回归分析得出回归系数及常数，建立血压测量回归方程式。经回归分析以及双弹性腔模型理论得出C点位于B和Z之间的血压方程如式(13)、(14)所示，C点位于Z与M之间血压方程如式(15)、(16)所示，其中，Ps为收缩压，Pd为舒张压。

Ps=348.852 1-1.389 2*PWTT，

(13)

Pd=20.883 6+0.402 6*Ps，

(14)

Ps=359.532 6-1.523 1*PWTT，

(15)

Pd=5.235 5+0.602 1*Ps。

(16)

3 实验分析

为验证所建立的血压估算模型的可靠性，本研究采用水银血压计与本系统对测试者同时进行检测以分析该设计的准确性。其中，以水银血压计测得数据作为真值，本设计测量数据为测值。在室温25 ℃下，对年龄(18～65)及性别各不相同的150位测试者(包括110名血压正常的测试者和40名高血压患者)进行实验，测量得到的数据及相应的误差图如图8、图9所示。部分测量数据如表3所示(表中测值和真值单位为mmHg)。

由图8、图9及表3知，本系统测得血压值与真值之间均保持较好的一致性，数据误差均保持在±5以内，优于AAMI标准。两种加速脉搏波的分类具有较高准确性，具有一定的临床意义。

表3 实验测量数据
Tab.3 Measured data

序号性别/年龄 PS测值PS真值误差PD测值PD真值误差分类结果真实类别1M/209899-16063-3112M/25102107-56263-1113F/199695+16265-3224M/29110107+36564+1115F/27116119-37473+1226M/33118116+27580-5227F/33120125-57781-4228F/35123121+27877+1219M/27106110-46265-31110F/28108109-16972-32111M/39136138-27574+11112M/45144139+59188+32213M/43139142-37681-51114F/40145142+39288+42215F/42147146+19395-222

图8 采集数据分布图
Fig.8 Distribution graph of collected data

图9 测值和真值误差图
Fig.9 Error between measure data and true data

4 结论

结合EMD及形态学方法对光电容积脉搏波进行降噪处理，对比不同滤波方法的信噪比，发现本文方法滤波效果优于小波变换。为了验证Adaboost识别不同脉搏波的有效性，本文基于预测灵敏度、阳性率及准确率3个指标进行评价，实验表明Adaboost识别不同脉搏波具有一定的可靠性。为验证血压模型的有效性，对150位测试者进行测量，实验结果证明，论文建立的血压计算模型误差在±5 mmHg以内，具有一定的准确性和可靠性，达到了较好的效果。