基于VMD-NLMS的运动状态血氧监测算法

张 林,余成波,谭 拥,米有玉

(重庆理工大学电气与电子工程学院,重庆 401320)

1 引言

血氧饱和度(SpO2)是一项表征人体心肺功能及循环系统健康程度的重要生理参数[1]。近年来推出的可穿戴的脉搏血氧监测设备,可以通过采集PPG信号实时的监测用户的血氧饱和度,有效地避免一些突发的心血管疾病。但可穿戴设备在为用户带来便捷体验的同时,也引入了更大的运动干扰,而常见的血氧监测算法的鲁棒性较差,并不能对抗较强的运动噪声,使测量结果往往存在较大误差。因此,研究一种鲁棒性更强的血氧监测算法有很大的意义和价值。

目前去除PPG信号中的运动伪影的方法主要有自适应滤波和经验模态分解及其改进算法等。有陈真诚等人[2-4]利用EMD处理PPG信号时,由于其理论存在缺陷,分解过程易出现模态混叠、端点效应等问题[5],去噪效果并不理想。因此,有Kang等人[6-7]采用VMD处理PPG信号,克服了模态混叠和边界效应的影响,但重构时仅使用加速度信号作为参考选择IMF,使得处理后的PPG信号仍有噪声残留。此外,自适应滤波方法也常用于PPG信号处理。Han等人[8]将同步测量的加速度信号作为自适应滤波器的参考信号来滤除运动伪影,对比EMD降低了时间复杂度,但其去噪效果过于依赖参考信号的质量,因此常与其它方法结合使用[9-11],但结合算法仍使用加速度信号作为噪声参考,不仅增加了硬件成本,且构造的运动噪声参考信号较为片面,去噪效果虽有提高,但并不是十分理想。

为了有效地去除PPG信号中的运动伪影,本文提出了一种基于变分模态分解和归一化变步长的自适应滤波算法来消除运动干扰。首先通过VMD算法将预处理后的PPG信号分解为若干个本征模态分量,然后通过计算每个模态的多尺度排列熵值,选择合适的模态重构噪声参考信号,最后使用归一化的最小均方误差自适应算法滤除PPG信号中的运动伪影。

2 相关算法原理

2.1 VMD算法

变分模态分解[12]是Konstantin等人在2014年提出的一种自适应的信号处理方法,通过迭代搜寻变分模态的最优解,不断更新各模态函数及中心频率,得到若干具有一定带宽的模态函数。变分模态分解是EMD算法的改进,其优点是更好的避免了模态混叠和端点效应等问题。

VMD方法将信号分解成预先定义的K个本证模态分量IMF,这里的IMF即 被重新定义为一个AM-FM信号

uk(t)=Ak(t)cos(φk(t))

(1)

其中Ak(t)为瞬时幅值,Ak(t)≥0;φk(t)为相位;φk(t)对t求一阶微分,得到uk(t)的瞬时频率ωk(t)=φk(t)′。

假设每个IMF在信号中心频率附近具有有限的带宽,由此构造了如下的约束变分问题

(2)

为了求解式(2),引入二次惩罚项α和Lagrange乘数λ(t),把问题转换为了无约束的变分问题

L({uk},{ωk},λ)

(3)

利用交替方向乘子(Alternating Direction Method of Multipliers,ADMM)对式(3)求解得到的K个IMF的频域结果进行反傅里叶变换,得到其时域信号。

2.2 NLMS算法原理

归一化最小均方误差(NLMS)算法采用变步长的方法来缩短自适应收敛过程,解决了最小均方误差(Least Mean Square,LMS)算法在收敛速度和稳态误差对步长要求相互矛盾的问题[13],是一种改进的LMS算法。NLMS算法具体原理如下

(4)

e(n)=d(n)-y(n)=d(n)-wT(n)x(n)

(5)

(6)

x(n)输入信号,d(n)为参考信号,w(n)为滤波器抽头系数,e(n)为误差信号。μ为步长因子,为了保证自适应滤波器的性能,取值范围应满足0≤μ≤2;γ是为了避免xT(n)x(n)过小导致步长太大而预先设置的矫正值,一般取0≤γ≤1。

3 VMD-NLMS改进算法

本文提出的改进算法主要分为两个部分:基于VMD的运动噪声参考信号重构和基于NLMS的PPG信号的运动伪影消除。

3.1 基于VMD的参考噪声重构

用VMD算法对预处理后的PPG信号进行分解,得到K个中心频率由低到高的模态分量,而利用VMD方法进行参考噪声重构的算法思想就是从K个本征模态分量中选出包含运动伪影的模态,然后进行重构。因此,重构的首要问题就是建立相关的准则判断并选取出运动伪影对应的模态分量。本文通过计算IMF的多尺度排列熵MPE来选择重构信号的IMF。

多尺度排列熵是一种衡量时间序列复杂度的非线性动力学特征指标[14-15],是在排列熵基础上的改进和优化,通过对时间序列的粗粒化,得到不同时间尺度的熵值信息,文献[16]中描述了其原理及计算方法。而对IMF分量来说,计算得到的MPE值越小,说明该尺度下的时间序列越规则,越接近PPG信号的原始成分,因为无噪的PPG信号为一个周期信号;MPE值越大,时间序列越随机,越接近噪声成分。计算得到各IMF分量的MPE值后,选择合适的尺度因子得到判断IMF是否为噪声分量的阈值,然后根据设定的阈值选取合适的IMF分量进行重构即可得到噪声参考信号。

3.2 基于NLMS的运动伪影消除

将重构的噪声参考信号和原PPG信号共同作为NLMS自适应滤波器的输入,得到去除运动伪影后的PPG信号。整个算法流程图如图1所示。

图1 VMD-NLMS算法流程图

4 数据的获取与预处理

4.1 PPG信号采集设备的搭建

本文利用SFH7050传感器、MPU6050六轴加速度传感器、AFE4404模拟前端、主控芯片STM32F103及蓝牙传输模块等模块搭建了PPG信号的采集设备如图2所示。其中SFH7050传感器搭载了血氧检测所需的波长为940 nm的红外光LED和波长为660 nm的红光LED[17],而模拟前端AFE4404高度集成了LED驱动电路、放大电路、滤波电路、采样电路等血氧前端采集必需的功能模块。

图2 PPG信号采集设备

为了探究不同运动状态下的PPG信号,本文用设计的PPG信号采集装置采集了常见的三种运动状态:步行、跑步和骑行时腕部所产生的PPG信号,且在采集过程中,尽量保证手腕自然摆动,采样的频率设置为256 Hz。如图3为采集的行走状态下的双路PPG信号。

图3 步行状态下的双路PPG信号

由图3可知,采集到的未经处理的PPG信号中噪声较多。因此,若使用未处理的PPG信号进行血氧饱和度值的计算,结果与真实值会出现很大的误差,甚至出现错误的结果。

4.2 PPG信号预处理

由于人体每分钟的心率正常为40～200 次/分钟,对应的频率范围大约为0.7～3.5 Hz。而运动伪影也在这个频率范围附近。为了充分保留有用的信息并去掉其它一些不感兴趣的噪声如基线漂移和高频噪声等,首先对PPG信号进行预处理,采用长度为8 秒,步长为2 秒的滑动时间窗口截取PPG信号;然后使用截止频率为0.7～3.5 Hz的带通滤波器对信号进行处理。由于巴特沃斯带通滤波器在其通频带内有平稳的幅频特性,滤波器本文选择了巴特沃斯带通滤波器,设置通带频率为0.7 Hz～3.5 Hz,阻带截止频率为0.5 Hz～5.0 Hz,通带最大衰减系数为3 dB,阻带最小衰减系数为18 dB,PPG信号的预处理结果如图4所示。

图4 预处理前后的PPG信号

由图4可以看到,预处理前的PPG信号的上面有呼吸谐波等因素造成的表现为基线漂移的低频噪声的干扰,有工频干扰等其它一些因素造成的表现为“毛刺”状的高频噪声的干扰,还有由于采集时的运动干扰带来的表现为幅值的大幅增大或减小的运动噪声的干扰。而经过预处理后的PPG信号,从时域特征来看,高频噪声带来的信号的“毛刺”得到了很好的抑制,基线的漂移也有了一定程度的抑制,但运动噪声并未得到有效地抑制。因此,预处理后的PPG信号需要进一步的去除信号中的噪声,才能进行血氧饱和度的计算。

5 实验结果及分析

5.1 去噪结果及对比分析

为了进一步的去除PPG信号中的运动噪声,本文对预处理后的PPG信号进行VMD分解,而VMD算法的分解效果的好坏,主要取决于模态个数K和惩罚因子α的取值。K取值过大,会造成信号的过分解;而K取值过小,则会造成模态混叠,即不同中心频率的信号成分经过VMD分解后处在了同一个本征模态分量中,产生“虚假”分量。而α的取值主要影响VMD分解后的得到各模态的带宽大小模态混叠程度。如果惩罚因子α的值太小,则每个频率的带宽太大,并且将相似频率误认为同一中心频率的大带宽部分,从而将其分解为同一分量,出现模态混叠的问题。当该值太大时,部分频带宽度减小,可能会导致信号过分解,出现虚假分量。K的取值可由待处理信号的频谱中峰值的个数来确定,因此根据预处理后信号的频谱中峰的个数,本文中K设置为7。预处理后的信号经过变分模态分解后得到了IMF1-IMF7共7个本证模态分量如图5所示,各个本征模态分量的中心频率与处理前信号的峰值所对应的频率分别如表1所示,对比分解前后各模态的中心频率(cf),与分解前信号的主要频率(mf),发现频率值基本对应,可见VMD的结果可以较好反应原信号的频率特征。

表1 PPG信号的主要频率与IMF的中心频率

图5 PPG信号VMD分解得到的IMF分量

惩罚因子α的取值可通过信号分解后各模态在频域的混叠程度来判断。如图6所示为K取7,α取3300时,PPG信号经分解后各个IMF的频率-幅值图,由图可得,模态直接并没有发生混叠,所以,α可按经验取值为3300。

图6 各个IMF分量的频率-幅值图

PPG信号经过VMD分解后得到7个IMF。本文通过计算各IMF分量的MPE来确定各IMF分量所代表的频率成分。计算MPE时,嵌入维数m、时间延迟τ、尺度因子s这些参数的选择会对熵值大小造成影响。

时间延迟τ的选择对熵值的计算影响较小,本文计算了τ为1～5时,单个模态的排列熵值如图7所示。可以看到,τ的变化对熵值的大小的影响较小,因此本文中,时间延迟τ设置为1。嵌入维数m的取值一般为3～7[18]。m过大时,信号在相空间中重构时会花费更多的时间,计算时间会因此而大大增加;m过小又会使信号在相空间重构时信息太少而丢失信号的细节特征。本文统计了当时间延迟τ取1,m取3～7时计算单个排列熵所用的时间如表2。考虑到计算效率和特征损失的平衡,本文的嵌入维数设置为6。尺度因子s的值对熵值的计算有较大的影响。若s取值过小,则不能最大限度的提取信号的特征信息;若s取值过大,可能抹除信号之间的复杂度差异,还会导致熵值计算不稳定。本文计算了s从1到20时,7个IMF的熵值如图8所示,当s大于7后各个IMF的排列熵值开始重合,无法反映出IMF的表征特征;当s为4～6时,各IMF分量的排列熵值区别较为明显,其熵值如表3所示。但为了在保证熵值稳定性的同时,最大限度的提取信号的特征信息,将本文的尺度因子s设置为6。

表2 不同嵌入维数下排列熵的计算时间

表3 IMF在尺度因子不同时的多尺度排列熵

图7 单个IMF在不同时间延迟下的排列熵

图8 各IMF不同尺度下的排列熵值

本文通过验证最终的去噪效果,来确定MPE的最佳阈值。经50次试验后,用来构造噪声参考信号的IMF排列熵阈值范围分布在(0,0.2320)∪(0.3450,+∞)的次数为43次,因此将阈值设置为σmin≤0.2320和σmax≥0.3450,即当MPE≤0.2320 时,为低频噪声成分,当MPE≥0.3450 时,为高频噪声成分,此时重构的噪声参考信号最好。将重构的参考信号和预处理后的PPG信号共同输入NLMS自适应滤波器中去噪得到的干净的PPG信号。

此外,本文选择了用加速度信号作为参考选择IMF进行PPG信号重构的EMD算法(A-EMD)和VMD算法(A-VMD)对PPG信号进行处理,结果如图9所示。并计算了三种算法去噪前后的SNR和MSE如表4。

表4 三种方法处理PPG信号的SNR和MSE

图9 三种不同方法处理PPG信号的结果

由图9知,使用A-EMD方法处理PPG信号后,一些较强的运动伪影得到了滤除,但由于发生了模态混叠,导致在信号重构时丢失了部分有用的信息,在4 s～8 s出现较大的波形失真;使用A-VMD处理PPG信号后,运动伪影得到了较好的抑制,但信号丢失了一些峰值特征,且运动伪影较强时,抑制效果不够明显;使用本文方法处理PPG信号后,不仅有效的抑制了运动伪影,而且还保留了信号在一些细节处的特征。由表4可知,本文方法对比另外两种算法SNR提高了1 dB～2 dB,MSE也更小,因此去噪效果更好。图10为不同运动状态下通过本文方法对PPG信号的处理结果。

图10 VMD-NLMS处理前后的PPG信号和生成的参考信号

图10中,红线为红外光PPG信号,黑色线为红光PPG信号,可以看出,经过本文算法处理后,三种运动状态的PPG信号在保持脉搏信号时域特征的同时,较好的抑制了运动伪影。

5.2 血氧饱和度计算结果及对比分析

本文通过提取去噪后双路PPG信号的直流分量和交流分量计算血氧饱和度值。实验对象为5名身体健康的人,分别采集1分钟内的步行、跑步、骑行三种状态下左手腕部的PPG信号,然后使用三种算法处理PPG信号并计算血氧饱和度值,同时选取标准的脉搏血氧仪同步检测(同步采集时右手尽可能保持静止)的血氧饱和度值作为标准值,结果如表5-7所示。由此可计算各方法测量血氧值的平均绝对误差值(Mean Absolute Deviation,MAE)和平均绝对误差百分比(Mean Absolute Error Percentage, MAEP)如表8所示。

表5 步行状态下的血氧饱和度计算值

表6 跑步状态下的血氧饱和度计算值

表7 骑行状态下的血氧饱和度计算值

表8 血氧饱和度值测量结果的误差….(%)

由表8可知相比与A-EMD和A-VMD两种方法,本文提出的算法平均绝对误差和平均绝对误差百分比减小了2.0%～3.0%,说明使用本文算法得到的血氧值更接近于真实值。

6 结论

本文提出了一种VMD-NLMS的方法来去除PPG信号中的运动伪影。通过VMD方法构造了噪声参考信号,作为自适应滤波器的输入来去除PPG信号中的运动伪影。该方法的优点在于自适应滤波器的参考信号是利用VMD分解从PPG信号内部产生的而不是使用三轴加速度信号作为噪声参考,不仅节省了硬件成本,而且重构的噪声参考信号比用单一加速度作为参考信号更有运动噪声代表性。通过对常见的3种运动状态下PPG信号的分析和处理,计算相应的血氧饱和度,证明了本算法更强的鲁棒性和准确性;通过与基于加速度的EMD算法和基于加速度的VMD算法去噪结果的对比,证明了本算法的优越性。因此,本文提出的算法可以很好的应用于脉搏血氧值的计算。