﻿ 基于鲸鱼优化算法的VMD-NLM脉搏信号滤波算法研究
 上海理工大学学报  2022, Vol. 44 Issue (6): 553-561 PDF

1. 上海理工大学 光电信息与计算机工程学院，上海 200093;
2. 上海中医药大学 基础医学院，上海201203

VMD-NLM pulse signal filtering algorithm based on whale optimization algorithm
YANG Haima1, CHEN Jiaci1, XU Xiaohan1, LI Fufeng2, SONG Zhichao1, JIN Yan1
1. School of Optical-Electrical and Computer Engineering, University of Shanghai for Science and Technology, Shanghai 200093, China;
2. College of Basic Medicine, Shanghai University of Traditional Chinese Medicine, Shanghai 201203, China
Abstract: To solve the problem that the collected pulse wave signal is easily disturbed by high frequency noise and baseline drift, which make subsequent pathological research and analysis difficult and the measurement accuracy of inferior quality, an improved filtering algorithm combining variational mode decomposition and non-local mean noise reduction was proposed. Aiming at the problem that the selection of variational modal parameters has different effects on the results, the whale optimization algorithm was used to adaptively select appropriate parameters, and the modal components were screened according to the permutation entropy results. The experimental results show that the signal-to-noise ratio and root mean square error of the improved signal are better than other noise reduction methods after being processed by the filtering algorithm.
Key words: pulse signal     whale optimization algorithm     variational mode decomposition     permutation entropy     non-local mean filtering

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)

 $\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)

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)

1.2 鲸鱼优化算法

 $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 {A}}=2 {\boldsymbol{{a}}} \cdot {\boldsymbol{r}}-{\boldsymbol{{a}}}$ (8)
 ${\boldsymbol C}=2{\boldsymbol{r}}$ (9)

 ${\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)

 $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)

1.3 非局部均值去噪

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

 $\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 排列熵

 $\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)

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

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

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

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

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

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

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

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

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

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

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

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

2.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)

3 实验分析 3.1 仿真实验

 图 1 干净脉搏波及其含噪仿真信号 Fig. 1 Clean pulse wave and its noisy simulation signal
 $y(t) = x(t) + 0.5\sin (0.2{\text{π}} t) + \lambda$ (21)

 图 2 仿真信号模态个数优化 Fig. 2 Optimization of the number of IMFs

 图 3 惩罚因子优化过程 Fig. 3 Optimization process of penalty factor

 图 4 仿真信号的分解效果图 Fig. 4 Decomposition effect diagram of simulation signal

 图 5 仿真信号去噪对比图 Fig. 5 Comparison diagram of simulated signal denoising

3.2 实例信号分析

 图 6 采集的原始数据 Fig. 6 The original collected data

 图 7 实例信号分解效果图 Fig. 7 Decomposition effect diagram of sample signal

 图 8 实例信号去噪效果对比 Fig. 8 Comparison of sample signal denoising

 图 9 不同方法的信噪比 Fig. 9 Signal-to-noise ratio of different methods

 图 10 不同方法的均方根误差 Fig. 10 Root mean square error of different methods
4 总　结

 [1] 郭睿, 王忆勤, 燕海霞, 等. 中医常见脉象的血液动力学参数分析[J]. 上海中医药大学学报, 2010, 24(6): 26-29. DOI:10.16306/j.1008-861x.2010.06.009 [2] PRADA E J A, GALLEGO C A B, GARCÍA J F C. On the development of an efficient, low-complexity and highly reproducible method for systolic peak detection[J]. Biomedical Signal Processing and Control, 2021, 68: 102606. DOI:10.1016/j.bspc.2021.102606 [3] 田泽懿, 张磊, 单新治, 等. 基于脉搏波传导时间的血压检测研究进展[J]. 光学仪器, 2020, 42(1): 88-94. [4] 韩庆阳, 王晓东, 李丙玉, 等. EEMD在同时消除脉搏血氧检测中脉搏波信号高频噪声和基线漂移中的应用[J]. 电子与信息学报, 2015, 37(6): 1384-1388. DOI:10.11999/JEIT141390 [5] 张培玲, 李小真, 崔帅华. 基于改进小波阈值-CEEMDAN算法的ECG信号去噪研究[J]. 计算机工程与科学, 2020, 42(11): 2067-2072. DOI:10.3969/j.issn.1007-130X.2020.11.020 [6] DRAGOMIRETSKIY K, ZOSSO D. Variational mode decomposition[J]. IEEE Transactions on Signal Processing, 2014, 62(3): 531-544. DOI:10.1109/TSP.2013.2288675 [7] KUMAR A, ZHOU Y Q, XIANG J W. Optimization of VMD using kernel-based mutual information for the extraction of weak features to detect bearing defects[J]. Measurement, 2021, 168: 108402. DOI:10.1016/j.measurement.2020.108402 [8] 顾文景, 周丽. 改进VMD算法在颤振试验信号模态参数辨识中的应用[J]. 振动工程学报, 2021, 34(2): 292-300. DOI:10.16385/j.cnki.issn.1004-4523.2021.02.009 [9] 崔善政, 郭艳珍, 梁钊, 等. 变分模态分解在去除心电图信号基线漂移中的应用[J]. 电子测量与仪器学报, 2018, 32(2): 167-171. DOI:10.13382/j.jemi.2018.02.024 [10] 韦海成, 蔡坤, 赵静. 改进VMD算法消除脉搏波基线漂移研究[J]. 电子测量与仪器学报, 2020, 34(8): 144-150. DOI:10.13382/j.jemi.B1902831 [11] 吴龙文, 聂雨亭, 张宇鹏, 等. 基于变分模态分解的自适应滤波降噪方法[J]. 电子学报, 2021, 49(8): 1457-1465. DOI:10.12263/DZXB.20190972 [12] MIRJALILI S, LEWIS A. The whale optimization algorithm[J]. Advances in Engineering Software, 2016, 95(12): 51-67. [13] 齐德明. 基于改进型的非局部均值滤波算法在医学图像处理中的研究与应用[J]. 计算机应用与软件, 2021, 38(9): 256-261,279. DOI:10.3969/j.issn.1000-386x.2021.09.040 [14] 郝东昊, 鲁华祥, 陈刚, 等. 抑制心电信号中工频干扰的滤波采样系统设计[J]. 传感器与微系统, 2018, 37(4): 94-97. DOI:10.13873/J.1000-9787(2018)04-0094-04 [15] 李志军, 张鸿鹏, 王亚楠, 等. 排列熵—CEEMD分解下的新型小波阈值去噪谐波检测方法[J]. 电机与控制学报, 2020, 24(12): 120-129. DOI:10.15938/j.emc.2020.12.015 [16] 陈真诚, 吴贤亮, 赵飞骏. EEMD结合小波阈值的光电容积脉搏波信号降噪[J]. 光学 精密工程, 2019, 27(6): 1327-1334.