脉搏信号是脉管搏动的轨迹，综合了心脏射血活动和脉搏沿血管树传播途中携带的各种信息。因而该信号的采集对于提前预知患者潜在的心血管病灶和针对心血管疾病采取有效治疗具有指导作用^[1-3]。然而采集到的脉搏信号经常会因为采集动作及仪器本身误差而出现不同程度的噪声，影响了信号的特征提取分析。因此，如何有效地减小噪声带来的影响，是当前脉诊客观化研究的一大热点。

针对这个问题，目前主要以小波降噪为理论基础开展研究，并在此基础上提出了多种改进方案，例如经验模态分解（empirical mode decomposition，EMD）、集合经验模态分解（ensemble empirical mode decomposition，EEMD）、互补集合经验模态分解（complementary ensemble empirical mode decomposition，CEEMD）等。但EMD方法容易产生模态混叠的问题，EEMD由于加入了白噪声可能对源信号产生一定的影响^[4]。张培玲等^[5]提出改进小波阈值结合完全自适应噪声集合经验模态分解（CEEMDAN）的方法有效地解决了噪声传递的影响，但依然无法完全克服模态混叠的问题。2014年，Dragomiretskiy等^[6]提出了变分模态分解（variational mode decomposition, VMD），该算法采用自适应非递归模型，很好地解决了经验模态分解所共有的模态混叠和端点效应问题。但是在VMD算法中，需要提前设置相关参数，不同的参数设置会带来不一样的结果^[7-8]。因此，如何获得两个参数的最优解并应用到脉搏信号中，是当前的一大问题。崔善政等^[9]提出了一种基于VMD去除心电信号基线漂移的方法，将VMD算法成功应用于心电脉搏信号的降噪研究中，并有效避免了经验模态分解带来的影响，但该方法通过设定频率阈值进行低频降噪，存在容易滤除脉搏波中有效成分的问题。韦海成等^[10]通过VMD在不同K值下的中心频率分布选取最佳K值，但未考虑到惩罚因子的影响，且基线用EMD去除会存在一定的模态混叠现象。

为解决上述问题，本文将鲸鱼优化算法应用于VMD参数选取中，对VMD中的K值与惩罚因子进行双参数优化。相比于其他传统的优化算法，鲸鱼优化算法需要设置的参数少，收敛速度快，具有较好的全局寻优能力，被广泛用在故障特征、水文分析、电磁信号等多个领域，且取得了不错的效果。故本文提出一种基于鲸鱼算法优化VMD参数选择过程，以排列熵作为选取噪声模态分量的评判标准，然后对含噪分量采用非局部均值滤波去噪方法，实现对脉搏信号噪声的有效滤除，最后将所有模态分量进行重构，获得降噪后的脉搏信号。

1 理论原理 1.1 变分模态分解

VMD是一种自适应非递归的模态分解模型，能够将原始信号分解成一组具有确定中心频率和一定带宽的模态分量。其借助了HHT（Hilbert-Huang transform）变换、频率混合等方法，通过迭代搜寻变分模型，相较于EMD能够有效地减小模态混叠影响，具有更强的鲁棒性^[11]。

模态分量变分约束模型的表达式如下：

$ \left\{ \begin{gathered} \mathop {\min }\limits_{\{ {u_k}\} ,\{ {\omega _k}\} } \left\{ {\sum\limits_k {\left\| {{\partial _t}\left[ {\left( {\delta \left( t \right) + \frac{{\rm j}}{{{\text{π}} t}}} \right) {u_k}\left( t \right)} \right]{{\rm e}^{ - {\rm j}\omega t}}} \right\|} _2^2} \right\} \\ s.t.\sum\limits_k {{u_k}(t) = f(t)} \\ \end{gathered} \right. $

(1)

式中：k为模态分解层数； $ {u}_{k}\left(t\right) $ 为第k个模态分量； ${\omega}_{k}\left(t\right)$ 为第k个模态分量的中心频率； $ f\left(t\right) $ 为原始信号。为求解式(1)需要引入拉格朗日算子，可将变分模型求解转换为无约束变分模型，通过k值的计算得到各个模态分量。该算子表达式如下：

$\begin{split}& L({{{u}}_{k}\text{},{}{\omega }_{k}\text{},}\lambda )=\alpha {{\displaystyle \sum _{k}\left\Vert {\partial }_{t}\left[\left(\delta (t)+\frac{{\rm{j}}}{\text{π}t}\right) {u}_{k}(t)\right]{\rm e}^{-{\rm{j}}\omega t}\right\Vert }}_{2}^{2}+\\&\quad{\left\Vert s(t)-{\displaystyle \sum _{k}{u}_{k}(t)}\right\Vert }_{2}^{2} +\left\langle \lambda (t),f(t)-{\displaystyle \sum _{k}{u}_{k}(t)}\right\rangle\\[-12pt] \end{split} $

(2)

式中： $ \alpha $ 为惩罚因子； $ \lambda \left(t\right) $ 为拉格朗日乘法算子。为得到最优解，通过交替方向乘子法ADMM进行迭代更新，具体步骤如下：

a. 初始化 $ \left\{{u}_{k}^{1}\right\} $ ， $ \left\{{\omega }_{k}^{1}\right\} $ ， $ {\lambda }^{1} $ ，此时迭代次数为0；

b. 根据式（3）、式（4）、式（5）更新 $ \left\{{u}_{k}\right\} $ ， $ \left\{{\omega }_{k}\right\} $ 和 $ \lambda $ ；

$ \left\{ {\hat u_k^{n + 1}} \right\} = \dfrac{{\widehat f(\omega ) - \displaystyle \sum\limits_{i \ne k} {\hat u_i^{n + 1}(\omega ) + \dfrac{{{{\hat \lambda }^n}(\omega )}}{2}} }}{{1 + 2\alpha {{(\omega - \omega _k^n)}^2}}} $

(3)

$ \omega _{{k}}^{n + 1} = \frac{{ \displaystyle {\int}_0^\infty {\omega {{\left| {\hat u_k^{n + 1}(\omega )} \right|}^2}{\rm d}\omega } }}{{ \displaystyle {\int}_0^\infty {{{\left| {\hat u_k^{n + 1}(\omega )} \right|}^2}{\rm d}\omega } }} $

(4)

$ {\hat \lambda ^{{{n + 1}}}}(\omega ) = {\hat \lambda ^n} + \beta \left[ {\hat f(\omega ) - \sum\limits_k {\hat u_k^{n + 1}(\omega )} } \right] $

(5)

c. 不断重复步骤b进行迭代，直到函数满足以下条件：

$ \sum\limits_{{k}} {\frac{{\left\| {u_k^{n + 1} - u_k^n} \right\|_2^2}}{{\left\| {u_k^n} \right\|_2^2}}} < \varepsilon $

(6)

式中， $ \varepsilon $ 为判断精度。

1.2 鲸鱼优化算法

鲸鱼优化算法（whale optimization algorithm，WOA）是Mirjalili等^[12]于2016年提出的，该算法从鲸鱼捕食猎物的行为中获得启发，通过模拟座头鲸一边收缩包围圈，一边螺旋向上进行狩猎的行为，以实现优化搜索的目的。鲸鱼群在每一次寻找猎物时，都会不断地更新包围圈，再不断地螺旋向上，整个过程就像蜿蜒的山路，直到山顶便停止搜索。该行为可以概括为如下的数学模型：

$ M=\left\{\begin{array}{l}{\boldsymbol{D}}\text=\left|{\boldsymbol{C}}·{{\boldsymbol{X}}}^{*}(t)-{\boldsymbol{X}}(t)\right|\\ {\boldsymbol{X}}(t+1)={{\boldsymbol{X}}}^{*}(t)-{\boldsymbol{A}}·{\boldsymbol{D}}\end{array}\right. $

(7)

式中： $ \boldsymbol{D} $ 为当前鲸鱼与猎物之间位置的参数； $ \boldsymbol{X}\left(t\right) $ 表示迭代t次后当前鲸鱼的位置； $ {\boldsymbol{X}}^{\mathrm{*}}\left(t\right) $ 表示当前猎物的所在位置； $\boldsymbol{A}·\boldsymbol{D}$ 表示鲸鱼每次游动的距离，也就是步长大小。 $ \boldsymbol{C} $ 和 $ \boldsymbol{A} $ 为系数向量，计算公式为

$ {\boldsymbol {A}}=2 {\boldsymbol{{a}}} \cdot {\boldsymbol{r}}-{\boldsymbol{{a}}}$

(8)

$ {\boldsymbol C}=2{\boldsymbol{r}} $

(9)

式中： ${a}$ 在迭代过程中会从2到0线性衰减；r表示[0,1]之间的随机数。

鲸鱼在捕食猎物时在螺旋上升的同时，会产生泡泡网进行攻击，因而鲸鱼算法假设该两种方式发生的比例为1∶1，可将该方法描述为

$ {\boldsymbol{X}}({t}+1)=\left\{\begin{array}{l}{{\boldsymbol{X}}}^{*}(t)-{\boldsymbol{A}}·{\boldsymbol{D}}\text{}\text{}\text{}\;\;\;p < 0.5且\left|{\boldsymbol{A}}\right| < 1\\ {{\boldsymbol{D}}}^{\prime}{\rm e}^{bl}\mathrm{cos}(2{\text{π}} l)+{{\boldsymbol{X}}}^{*}(t)\text{}\text{}\;\;\;p\geqslant 0.5\text{}\end{array}\right. $

(10)

式中：b定义了螺旋形状的常量；l为[−1,1]的随机数；p为[0,1]的随机数，表示螺旋更新概率。

在狩猎阶段，除了气泡攻击外，还需要通过随机的方式去搜寻猎物，触发条件为p<0.5且|A|>1。此时个体根据不同的鲸鱼位置不断进行信息更新，从而找到全局最优解，该过程可以描述为

$ N=\left\{\begin{array}{l}{\boldsymbol{D }}=\left|{\boldsymbol{C}}·{\boldsymbol{X}}_{\rm{rand}}-{\boldsymbol{X}}(t)\right|\\ {\boldsymbol{X}}(t+1)={\boldsymbol{X}}_{\rm{rand}}-{\boldsymbol{A}}·{\boldsymbol{D}}\end{array}\right. $

(11)

式中： ${{\boldsymbol{X}}_{\rm{rand}}}$ 为鲸鱼的随机位置向量。

1.3 非局部均值去噪

非局部均值去噪（non-local means，NLM）通过选中待处理的像素点，通过设置搜寻区域，对该区域内具有相似特征的像素点进行加权平均，在尽可能减小对图像本身损失的同时，进行降噪滤波^[13]。通过对确定的含有噪声的模态分量（intrinsic mode function，IMF）应用NLM算法进行去噪，以达到消除脉搏波形噪声的效果。其具体原理，可描述为

$ \widehat X(i) = \frac{1}{{Z(i)}}\sum\limits_{j \in {\Omega _i}} {\omega (i,j)x(i)} $

(12)

式中：Z(i)为归一化参数。权重需满足两个条件

$ \left\{ \begin{gathered} 0 \leqslant \frac{{\omega (i,j)}}{{Z(i)}} \leqslant 1 \\ \sum\limits_j {\frac{{\omega (i,j)}}{{Z(i)}}} = 1 \\ \end{gathered} \right. $

(13)

在搜寻像素点中：若两个样本的相似度较高，则该权值就越大；相似度低，权值就小。权值计算为

$ \omega (i,j) = {{\rm {e}}^{\tfrac{{ - \displaystyle \sum\limits_{\delta \in \vartriangle } {{{(y(i + \delta ) - v(j + \delta ))}^2}} }}{{{h^2}}}}} $

(14)

2 改进VMD降噪算法 2.1 排列熵

为了滤除脉搏信号中存在的噪声，更好地保留脉搏噪声的有效信息，需要对经由VMD分解得到的各模态分量建立评判标准。本文采用排列熵（permutation, entropy，PE）作为区分噪声分量的标准。排列熵能够很好地展现时间序列的复杂性，同时又兼具计算简单、抗噪能力强等优点^[14]。假设一时间序列 $X(i)=\{x(1),x(2),\cdots, x(n)\}$ ，对该序列中任一元素进行空间重构可得到重构矩阵如下：

$ \boldsymbol X = \left| \begin{gathered} x(1){\text{ }}x(1 + \tau){\text{ }} \cdots {\text{ }}x(1 + (m - 1)\tau) \\ x(2){\text{ }}x(2 + \tau){\text{ }} \cdots {\text{ }}x(2 + (m - 1)\tau) \\ \vdots\qquad \vdots \qquad\qquad \qquad \vdots\;\;\;\;\;\;\;\; \\ x(K){\text{ }}x(K + \tau){\text{ }} \cdots {\text{ }}x(K + (m - 1)\tau) \\ \end{gathered} \right| $

(15)

式中：m为嵌入维数；τ为延迟时间。对矩阵X进行重排后可得到

$ x(i + ({j_1} - 1)\tau) \leqslant x(i + ({j_2} - 1)\tau) \leqslant \cdots \leqslant x(i + ({j_m} - 1)\tau) $

(16)

式中： $ {j}_{1} $ ， $ {j}_{2} $ ， $\cdots$ ， $ {j}_{m} $ 表示各元素所在列索引。可通过排列组合定理获得所有情况的概率，即为各模态分量的PE值，如下：

$ {H_p}(m) = - \sum\limits_{i = 1}^k {{p_i}\ln {p_i}} $

(17)

式中，p_i为式(15)所列的矩阵中每行序列出现的概率。归一化后，可得

$ {H_p} = \frac{{{H_p(m)}}}{{\ln (m!)}} $

(18)

综上所述， $ {H}_{p} $ 越大，反映该分量随机性越大；相反该分量随机性就越小。

2.2 基于WOA-VMD-NLM的联合去噪方法

本文通过WOA对VMD的参数进行优化，其具体步骤如下：

a.初始化WOA参数，设定初始种群大小，模态分解数 $ K $ 以及惩罚因子 $ \alpha $ 等。

b.初始化种群，在指定的范围内随机生成初始化群体。

c.进行VMD分解，并计算种群适应度。本文将样本熵作为目标参数。

d. 开始迭代，寻找最优解。根据 $ p $ 与 $\boldsymbol A$ 的不同取值，根据式（10）、（11）可计算出更新后的个体位置。计算更新后的适应度值。

e.结束迭代过程。如果当前迭代次数t并未达到初始化设定的迭代次数时，会继续寻找当前最优解。否则迭代结束，获取此时的最佳参数 $ K $ 和 $ \alpha $ 。

f.确定模态分量的性质。通过计算各模态分量的排列熵，通过排列熵的分布规律确定噪声分量。

g.对确定的噪声分量进行非局部均值去噪。

h.将信号主导的模态分量与去噪后的分量进行重构，获得滤波后的脉搏信号。

2.3 降噪效果评价指标

采用信噪比SNR、均方根误差RMSE两个指标来评判该降噪算法的效果。信噪比反映了输出信号与噪声信号的功率比，可以用来表明该算法的降噪性能^[15]。均方根误差可用于对该信号在降噪前后差异性的一个评价。3个指标的公式如下：

$ SNR = 10\lg \left(\frac{{{P_{\rm s}}}}{{{P_{\rm n}}}}\right) $

(19)

$ RMSE = \sqrt {\frac{{ \displaystyle \sum {{{[x(t) - x'(t)]}^2}} }}{N}} $

(20)

式中： ${P}_{\rm s}$ 为有效信号； ${P}_{\rm n}$ 为噪声信号；x(t)为原始信号； $ x\left(t\right)\mathrm{\text{'}} $ 为降噪后信号；N为采样个数。因而当信噪比越大、均方根误差越小时，表明信号的降噪性能越为优异。

3 实验分析 3.1 仿真实验

为了验证该算法在脉搏波上的可行性，以干净的脉搏波信号 $x(t)$ 作为有效信号，为了模拟信号的肌电噪声，以40 dB信噪比的高斯噪声 $\lambda $ 进行替代，采用频率为0.1 Hz的正弦波信号用于模拟脉搏信号的基线漂移，如式（21）所示。图1为干净的信号及其仿真信号。

$ y(t) = x(t) + 0.5\sin (0.2{\text{π}} t) + \lambda $

(21)

现将图1中的脉搏仿真信号进行基于鲸鱼优化算法选取VMD的参数。每一次的迭代过程，都是在不断改变选取的参数值，以求达到最优解。该样本的参数优化过程如图2和图3所示。

可以看出，当迭代次数达到3次时，鲸鱼优化过程趋于稳定，最终获取的最优参数组合为 $ K $ =10， $ \alpha $ =5 356。最终的分解效果如图4所示。

由图4可知，原始的脉搏信号经由VMD分解后产生10个模态分量，其中部分模态属于高频噪声信号，部分属于低频基漂噪声，而对于不同的脉搏信号，VMD分解的模态分量数是不同的，因而需要一个指标来衡量分解出的模态属性。本文采用排列熵对分解的模态进行量化观察，来判别噪声分量和有效信号的模态。在排列熵计算中，需要设定嵌入维数m和时间延迟 $ \tau $ ，根据经验，本文取m=3， $\tau $ =10^[16]。如表1所示为IMF分量的PE值分布，经过对多段脉搏信号的计算，最终选定当PE值高于1.1时为高频噪声信号。此外根据频谱图可知，经由VMD分解后的第一个模态的峰值频率为0.122 Hz，属于加入的基线漂移信号，可认为是脉搏信号的基线漂移。

在提取完噪声信号后，对高频噪声模态分量进行NLM非局部均值滤波，然后将基线漂移分量进行剔除，最后将处理后的模态分量与其他模态分量进行重构，获得降噪后的信号。

为直观体现该方法的降噪效果，本文将小波阈值降噪、EEMD–小波阈值和CEEMDAN–小波阈值降噪算法进行了比较，如图5所示。

图5（a）是本文所采用的方法得出的滤波效果，相较于图1中的原始信号，该方法较好地保持了原有信号的信号特征，并且获得了更为平滑的数据曲线。图5（b）中是小波阈值降噪法，右图中可以看出该降噪后信号存在明显的波形失真。图5（c）采用的是EEMD结合小波阈值的方法^[16]，该方法在整体的信号平滑上同样有着不错的表现，但其对于脉搏波的基线漂移滤除效果不佳。图5（d）采用的是CEEMDAN结合小波阈值的去噪效果^[14]，可以看出该方法对脉搏信号的基线漂移进行了一定的降噪处理，但降噪后的脉搏波依然存在一定的毛刺信号，尤其是在脉搏周期的峰谷点处依然存在毛刺噪声，这会影响对脉搏特征点的提取，在局部降噪表现上不佳。因此从图像上看，本方法相较于其他方法具有更好的滤波效果。

为了量化说明本算法的滤波效果，首先采用相关系数进行评价，与标准的干净信号进行对比。如表2所示相关系数结果，可以从中看出，改进的VMD算法与原始信号的相关程度更高，从而说明本方法对于原始干净信号具有较好的还原度，保留了较多的细节。

3.2 实例信号分析

为了进一步验证本文提出算法对实测脉搏信号的降噪效果，本文采用自行研制的便携式脉搏信号采集装置作为本文的脉搏信号样本来源，并挑选了10名被试者参与脉搏信号的采集。由于脉搏信号受抖动、呼吸等影响较为明显，会产生一定的基线漂移，故本实验选择安静的环境，采用自主研制的便携式脉搏采集装置分别在同一时间段、同一房间内进行实验。测试者在平静状态下进行实验，戴上采脉装置，手臂保持水平放置，尽量让桡动脉与心脏处于同一水平面。本次实验采集1 min的脉搏信号，并取其中的1000个数据点作为实验对象，其余数据可持续验证实验。如图6所示为本次采集的某位被测者的脉搏信号。为方便计算，将脉搏信号进行归一化处理。

同样对脉搏信号进行WOA-VMD分解以及噪声分量评判，经由VMD分解后产生11个模态分量，如图7所示为分解效果图，表3为实例信号的排列熵分布。

从图中可以看到，实例信号经由VMD分解后提取到了信号的低频噪声，符合仿真信号得到的结论，然后根据排列熵分布对信号分量进行处理。同样为了能够反映本文提出的算法对实例信号具有同样的效果，将去噪效果与多个算法进行对比，如图8所示。

图8（a）是本文所采用的方法得出的滤波效果，可以看到信号在各个波峰谷点处较为平滑，对于基线漂移的滤除较为优异。图8（b）中是小波阈值降噪法，可以看出该方法对于信号降噪后存在明显的波形失真，整体波形显得不平滑，例如第一个波峰处与原信号不符。图8（c）为EEMD结合小波阈值去噪，可以看到该方法对于脉搏波的基线漂移滤除效果不佳，这与仿真信号的结论一致，并且在第二个脉搏周期中的重搏波特征不明显（图中圈出），存在一定的有效信号损失，这会影响到后续脉搏特征提取的研究。图6（d）采用的是CEEMDAN结合小波阈值的去噪效果。可以看出该方法对脉搏信号的基线漂移进行了一定的降噪处理，但降噪后的脉搏波依然存在一定的毛刺信号，在局部降噪表现上不佳。因此从图像上看，本方法相较于其他方法具有更好的滤波效果。

为了量化说明本算法的滤波效果，本文采用信噪比和均方根误差进行说明。其中信噪比越大，均方根误差越小，表明信号的降噪性能越为优异。对于采集的10组数据，利用上述的3种算法计算出被试者信号的信噪比和均方根误差，其中由于EEMD去噪存在明显的基线漂移，因此并未对其进行指标计算。如图9和图10所示，为每一组脉搏实例信号所计算获得的评价指标图。可以看出，经由本文提出的降噪方法得到的评价指标在10组实验中均优于其他方法，因而可认为本文的算法具有更好的降噪效果，对于特征信息的保存更为优异。

4 总　结

本文提出了一种利用鲸鱼优化算法寻找VMD分解过程中所需的参数，通过与排列熵确定有效分量，结合非均值滤波方法滤除脉搏信号中存在的噪声。该方法利用鲸鱼优化的快速、自适应等优点寻找合适的模态分解数，有效解决了VMD分解中参数设置的问题。通过排列熵分布，能够有效地区分出噪声分量与有效分量。最后通过非局部均值滤波，滤除噪声信号。实验证明，相较于CEEMD、传统VMD等，本方法在脉搏波降噪上的性能要更为出色，对原信号具有更高的还原度，为后期结合特征提取、脉象分类等工作奠定基础，对脉诊客观化具有重要意义。然而VMD的双参数优化往往需要多次迭代，导致其运行效率低下，因此如何提高算法的运行效率还需要进一步深入的研究。