基于鯨魚(yú)優(yōu)化算法的VMD-NLM脈搏信號(hào)濾波算法研究

楊海馬， 陳嘉慈， 徐笑寒， 李福鳳， 宋智超， 金 焱

（1.上海理工大學(xué)光電信息與計(jì)算機(jī)工程學(xué)院，上海 200093；2.上海中醫(yī)藥大學(xué)基礎(chǔ)醫(yī)學(xué)院，上海 201203）

脈搏信號(hào)是脈管搏動(dòng)的軌跡，綜合了心臟射血活動(dòng)和脈搏沿血管樹(shù)傳播途中攜帶的各種信息。因而該信號(hào)的采集對(duì)于提前預(yù)知患者潛在的心血管病灶和針對(duì)心血管疾病采取有效治療具有指導(dǎo)作用[1-3]。然而采集到的脈搏信號(hào)經(jīng)常會(huì)因?yàn)椴杉瘎?dòng)作及儀器本身誤差而出現(xiàn)不同程度的噪聲，影響了信號(hào)的特征提取分析。因此，如何有效地減小噪聲帶來(lái)的影響，是當(dāng)前脈診客觀化研究的一大熱點(diǎn)。

針對(duì)這個(gè)問(wèn)題，目前主要以小波降噪為理論基礎(chǔ)開(kāi)展研究，并在此基礎(chǔ)上提出了多種改進(jìn)方案，例如經(jīng)驗(yàn)?zāi)B(tài)分解（empirical mode decomposition，EMD）、集合經(jīng)驗(yàn)?zāi)B(tài)分解（ensemble empirical mode decomposition，EEMD）、互補(bǔ)集合經(jīng)驗(yàn)?zāi)B(tài)分解（complementary ensemble empirical mode decomposition，CEEMD）等。但EMD方法容易產(chǎn)生模態(tài)混疊的問(wèn)題，EEMD由于加入了白噪聲可能對(duì)源信號(hào)產(chǎn)生一定的影響[4]。張培玲等[5]提出改進(jìn)小波閾值結(jié)合完全自適應(yīng)噪聲集合經(jīng)驗(yàn)?zāi)B(tài)分解（CEEMDAN）的方法有效地解決了噪聲傳遞的影響，但依然無(wú)法完全克服模態(tài)混疊的問(wèn)題。2014年，Dragomiretskiy等[6]提出了變分模態(tài)分解（variational mode decomposition,VMD），該算法采用自適應(yīng)非遞歸模型，很好地解決了經(jīng)驗(yàn)?zāi)B(tài)分解所共有的模態(tài)混疊和端點(diǎn)效應(yīng)問(wèn)題。但是在VMD算法中，需要提前設(shè)置相關(guān)參數(shù)，不同的參數(shù)設(shè)置會(huì)帶來(lái)不一樣的結(jié)果[7-8]。因此，如何獲得兩個(gè)參數(shù)的最優(yōu)解并應(yīng)用到脈搏信號(hào)中，是當(dāng)前的一大問(wèn)題。崔善政等[9]提出了一種基于VMD去除心電信號(hào)基線漂移的方法，將VMD算法成功應(yīng)用于心電脈搏信號(hào)的降噪研究中，并有效避免了經(jīng)驗(yàn)?zāi)B(tài)分解帶來(lái)的影響，但該方法通過(guò)設(shè)定頻率閾值進(jìn)行低頻降噪，存在容易濾除脈搏波中有效成分的問(wèn)題。韋海成等[10]通過(guò)VMD在不同K值下的中心頻率分布選取最佳K值，但未考慮到懲罰因子的影響，且基線用EMD去除會(huì)存在一定的模態(tài)混疊現(xiàn)象。

為解決上述問(wèn)題，本文將鯨魚(yú)優(yōu)化算法應(yīng)用于VMD參數(shù)選取中，對(duì)VMD中的K值與懲罰因子進(jìn)行雙參數(shù)優(yōu)化。相比于其他傳統(tǒng)的優(yōu)化算法，鯨魚(yú)優(yōu)化算法需要設(shè)置的參數(shù)少，收斂速度快，具有較好的全局尋優(yōu)能力，被廣泛用在故障特征、水文分析、電磁信號(hào)等多個(gè)領(lǐng)域，且取得了不錯(cuò)的效果。故本文提出一種基于鯨魚(yú)算法優(yōu)化VMD參數(shù)選擇過(guò)程，以排列熵作為選取噪聲模態(tài)分量的評(píng)判標(biāo)準(zhǔn)，然后對(duì)含噪分量采用非局部均值濾波去噪方法，實(shí)現(xiàn)對(duì)脈搏信號(hào)噪聲的有效濾除，最后將所有模態(tài)分量進(jìn)行重構(gòu)，獲得降噪后的脈搏信號(hào)。

1 理論原理

1.1 變分模態(tài)分解

VMD是一種自適應(yīng)非遞歸的模態(tài)分解模型，能夠?qū)⒃夹盘?hào)分解成一組具有確定中心頻率和一定帶寬的模態(tài)分量。其借助了HHT（Hilbert-Huang transform）變換、頻率混合等方法，通過(guò)迭代搜尋變分模型，相較于EMD能夠有效地減小模態(tài)混疊影響，具有更強(qiáng)的魯棒性[11]。

模態(tài)分量變分約束模型的表達(dá)式如下：

式中：k為模態(tài)分解層數(shù)；uk(t)為第k個(gè)模態(tài)分量；ωk(t)為 第k個(gè)模態(tài)分量的中心頻率；f(t)為原始信號(hào)。為求解式(1)需要引入拉格朗日算子，可將變分模型求解轉(zhuǎn)換為無(wú)約束變分模型，通過(guò)k值的計(jì)算得到各個(gè)模態(tài)分量。該算子表達(dá)式如下：

式中：α為懲罰因子；λ(t)為拉格朗日乘法算子。為得到最優(yōu)解，通過(guò)交替方向乘子法ADMM進(jìn)行迭代更新，具體步驟如下：

b.根據(jù)式（3）、式（4）、式（5）更新{uk}，{ωk}和λ；

c.不斷重復(fù)步驟b進(jìn)行迭代，直到函數(shù)滿足以下條件：

式中，ε為判斷精度。

1.2 鯨魚(yú)優(yōu)化算法

鯨魚(yú)優(yōu)化算法（whale optimization algorithm，WOA）是Mirjalili等[12]于2016年提出的，該算法從鯨魚(yú)捕食獵物的行為中獲得啟發(fā)，通過(guò)模擬座頭鯨一邊收縮包圍圈，一邊螺旋向上進(jìn)行狩獵的行為，以實(shí)現(xiàn)優(yōu)化搜索的目的。鯨魚(yú)群在每一次尋找獵物時(shí)，都會(huì)不斷地更新包圍圈，再不斷地螺旋向上，整個(gè)過(guò)程就像蜿蜒的山路，直到山頂便停止搜索。該行為可以概括為如下的數(shù)學(xué)模型：

式中：D為當(dāng)前鯨魚(yú)與獵物之間位置的參數(shù)；X(t)表 示迭代t次后當(dāng)前鯨魚(yú)的位置；X*(t)表示當(dāng)前獵物的所在位置；A·D表示鯨魚(yú)每次游動(dòng)的距離，也就是步長(zhǎng)大小。C 和A為系數(shù)向量，計(jì)算公式為

式中：a在迭代過(guò)程中會(huì)從2到0線性衰減；r 表示[0,1]之間的隨機(jī)數(shù)。

鯨魚(yú)在捕食獵物時(shí)在螺旋上升的同時(shí)，會(huì)產(chǎn)生泡泡網(wǎng)進(jìn)行攻擊，因而鯨魚(yú)算法假設(shè)該兩種方式發(fā)生的比例為1∶1，可將該方法描述為

式中：b定義了螺旋形狀的常量；l為[-1,1]的隨機(jī)數(shù)；p為[0,1]的隨機(jī)數(shù)，表示螺旋更新概率。

在狩獵階段，除了氣泡攻擊外，還需要通過(guò)隨機(jī)的方式去搜尋獵物，觸發(fā)條件為p＜0.5且|A|＞1。此時(shí)個(gè)體根據(jù)不同的鯨魚(yú)位置不斷進(jìn)行信息更新，從而找到全局最優(yōu)解，該過(guò)程可以描述為

式中： Xrand為鯨魚(yú)的隨機(jī)位置向量。

1.3 非局部均值去噪

非局部均值去噪（non-local means，NLM）通過(guò)選中待處理的像素點(diǎn)，通過(guò)設(shè)置搜尋區(qū)域，對(duì)該區(qū)域內(nèi)具有相似特征的像素點(diǎn)進(jìn)行加權(quán)平均，在盡可能減小對(duì)圖像本身?yè)p失的同時(shí)，進(jìn)行降噪濾波[13]。通過(guò)對(duì)確定的含有噪聲的模態(tài)分量（intrinsic mode function，IMF）應(yīng)用NLM算法進(jìn)行去噪，以達(dá)到消除脈搏波形噪聲的效果。其具體原理，可描述為

式中：Z(i)為歸一化參數(shù)。權(quán)重需滿足兩個(gè)條件

在搜尋像素點(diǎn)中：若兩個(gè)樣本的相似度較高，則該權(quán)值就越大；相似度低，權(quán)值就小。權(quán)值計(jì)算為

2 改進(jìn)VMD降噪算法

2.1 排列熵

為了濾除脈搏信號(hào)中存在的噪聲，更好地保留脈搏噪聲的有效信息，需要對(duì)經(jīng)由VMD分解得到的各模態(tài)分量建立評(píng)判標(biāo)準(zhǔn)。本文采用排列熵（permutation,entropy，PE）作為區(qū)分噪聲分量的標(biāo)準(zhǔn)。排列熵能夠很好地展現(xiàn)時(shí)間序列的復(fù)雜性，同時(shí)又兼具計(jì)算簡(jiǎn)單、抗噪能力強(qiáng)等優(yōu)點(diǎn)[14]。假設(shè)一時(shí)間序列X(i)={x(1),x(2),···,x(n)}，對(duì)該序列中任一元素進(jìn)行空間重構(gòu)可得到重構(gòu)矩陣如下：

式中：m為嵌入維數(shù)；τ為延遲時(shí)間。對(duì)矩陣X進(jìn)行重排后可得到

式中：j1， j2，···，jm表示各元素所在列索引。可通過(guò)排列組合定理獲得所有情況的概率，即為各模態(tài)分量的PE值，如下：

式中，pi為式(15)所列的矩陣中每行序列出現(xiàn)的概率。歸一化后，可得

綜上所述，Hp越大，反映該分量隨機(jī)性越大；相反該分量隨機(jī)性就越小。

2.2 基于WOA-VMD-NLM的聯(lián)合去噪方法

本文通過(guò)WOA對(duì)VMD的參數(shù)進(jìn)行優(yōu)化，其具體步驟如下：

a.初始化WOA參數(shù)，設(shè)定初始種群大小，模態(tài)分解數(shù)K 以及懲罰因子α等。

b.初始化種群，在指定的范圍內(nèi)隨機(jī)生成初始化群體。

c.進(jìn)行VMD分解，并計(jì)算種群適應(yīng)度。本文將樣本熵作為目標(biāo)參數(shù)。

d.開(kāi)始迭代，尋找最優(yōu)解。根據(jù)p 與A的不同取值，根據(jù)式（10）、（11）可計(jì)算出更新后的個(gè)體位置。計(jì)算更新后的適應(yīng)度值。

e.結(jié)束迭代過(guò)程。如果當(dāng)前迭代次數(shù)t并未達(dá)到初始化設(shè)定的迭代次數(shù)時(shí)，會(huì)繼續(xù)尋找當(dāng)前最優(yōu)解。否則迭代結(jié)束，獲取此時(shí)的最佳參數(shù)K 和α。

f.確定模態(tài)分量的性質(zhì)。通過(guò)計(jì)算各模態(tài)分量的排列熵，通過(guò)排列熵的分布規(guī)律確定噪聲分量。

g.對(duì)確定的噪聲分量進(jìn)行非局部均值去噪。

h.將信號(hào)主導(dǎo)的模態(tài)分量與去噪后的分量進(jìn)行重構(gòu)，獲得濾波后的脈搏信號(hào)。

2.3 降噪效果評(píng)價(jià)指標(biāo)

采用信噪比SNR、均方根誤差RMSE兩個(gè)指標(biāo)來(lái)評(píng)判該降噪算法的效果。信噪比反映了輸出信號(hào)與噪聲信號(hào)的功率比，可以用來(lái)表明該算法的降噪性能[15]。均方根誤差可用于對(duì)該信號(hào)在降噪前后差異性的一個(gè)評(píng)價(jià)。3個(gè)指標(biāo)的公式如下：

式中：Ps為 有效信號(hào)；Pn為噪聲信號(hào)；x(t)為原始信號(hào)； x(t)’為降噪后信號(hào)；N為采樣個(gè)數(shù)。因而當(dāng)信噪比越大、均方根誤差越小時(shí)，表明信號(hào)的降噪性能越為優(yōu)異。

3 實(shí)驗(yàn)分析

3.1 仿真實(shí)驗(yàn)

為了驗(yàn)證該算法在脈搏波上的可行性，以干凈的脈搏波信號(hào) x(t)作為有效信號(hào)，為了模擬信號(hào)的肌電噪聲，以40 dB信噪比的高斯噪聲λ進(jìn)行替代，采用頻率為0.1 Hz的正弦波信號(hào)用于模擬脈搏信號(hào)的基線漂移，如式（21）所示。圖1為干凈的信號(hào)及其仿真信號(hào)。

圖1 干凈脈搏波及其含噪仿真信號(hào)Fig.1 Clean pulse wave and its noisy simulation signal

現(xiàn)將圖1中的脈搏仿真信號(hào)進(jìn)行基于鯨魚(yú)優(yōu)化算法選取VMD的參數(shù)。每一次的迭代過(guò)程，都是在不斷改變選取的參數(shù)值，以求達(dá)到最優(yōu)解。該樣本的參數(shù)優(yōu)化過(guò)程如圖2和圖3所示。

可以看出，當(dāng)?shù)螖?shù)達(dá)到3次時(shí)，鯨魚(yú)優(yōu)化過(guò)程趨于穩(wěn)定，最終獲取的最優(yōu)參數(shù)組合為K =10，α=5 356。最終的分解效果如圖4所示。

圖2 仿真信號(hào)模態(tài)個(gè)數(shù)優(yōu)化Fig.2 Optimization of the number of IMFs

圖3 懲罰因子優(yōu)化過(guò)程Fig.3 Optimization process of penalty factor

由圖4可知，原始的脈搏信號(hào)經(jīng)由VMD分解后產(chǎn)生10個(gè)模態(tài)分量，其中部分模態(tài)屬于高頻噪聲信號(hào)，部分屬于低頻基漂噪聲，而對(duì)于不同的脈搏信號(hào)，VMD分解的模態(tài)分量數(shù)是不同的，因而需要一個(gè)指標(biāo)來(lái)衡量分解出的模態(tài)屬性。本文采用排列熵對(duì)分解的模態(tài)進(jìn)行量化觀察，來(lái)判別噪聲分量和有效信號(hào)的模態(tài)。在排列熵計(jì)算中，需要設(shè)定嵌入維數(shù)m和時(shí)間延遲τ，根據(jù)經(jīng)驗(yàn)，本文取m=3，τ=10[16]。如表1所示為IMF分量的PE值分布，經(jīng)過(guò)對(duì)多段脈搏信號(hào)的計(jì)算，最終選定當(dāng)PE值高于1.1時(shí)為高頻噪聲信號(hào)。此外根據(jù)頻譜圖可知，經(jīng)由VMD分解后的第一個(gè)模態(tài)的峰值頻率為0.122 Hz，屬于加入的基線漂移信號(hào)，可認(rèn)為是脈搏信號(hào)的基線漂移。

在提取完噪聲信號(hào)后，對(duì)高頻噪聲模態(tài)分量進(jìn)行NLM非局部均值濾波，然后將基線漂移分量進(jìn)行剔除，最后將處理后的模態(tài)分量與其他模態(tài)分量進(jìn)行重構(gòu)，獲得降噪后的信號(hào)。

為直觀體現(xiàn)該方法的降噪效果，本文將小波閾值降噪、EEMD-小波閾值和CEEMDAN-小波閾值降噪算法進(jìn)行了比較，如圖5所示。

圖4 仿真信號(hào)的分解效果圖Fig.4 Decomposition effect diagram of simulation signal

表1 各模態(tài)分量排列熵分布Tab.1 Arrangement entropy distribution of each IMF

圖5（a）是本文所采用的方法得出的濾波效果，相較于圖1中的原始信號(hào)，該方法較好地保持了原有信號(hào)的信號(hào)特征，并且獲得了更為平滑的數(shù)據(jù)曲線。圖5（b）中是小波閾值降噪法，右圖中可以看出該降噪后信號(hào)存在明顯的波形失真。圖5（c）采用的是EEMD結(jié)合小波閾值的方法[16]，該方法在整體的信號(hào)平滑上同樣有著不錯(cuò)的表現(xiàn)，但其對(duì)于脈搏波的基線漂移濾除效果不佳。圖5（d）采用的是CEEMDAN結(jié)合小波閾值的去噪效果[14]，可以看出該方法對(duì)脈搏信號(hào)的基線漂移進(jìn)行了一定的降噪處理，但降噪后的脈搏波依然存在一定的毛刺信號(hào)，尤其是在脈搏周期的峰谷點(diǎn)處依然存在毛刺噪聲，這會(huì)影響對(duì)脈搏特征點(diǎn)的提取，在局部降噪表現(xiàn)上不佳。因此從圖像上看，本方法相較于其他方法具有更好的濾波效果。

為了量化說(shuō)明本算法的濾波效果，首先采用相關(guān)系數(shù)進(jìn)行評(píng)價(jià)，與標(biāo)準(zhǔn)的干凈信號(hào)進(jìn)行對(duì)比。如表2所示相關(guān)系數(shù)結(jié)果，可以從中看出，改進(jìn)的VMD算法與原始信號(hào)的相關(guān)程度更高，從而說(shuō)明本方法對(duì)于原始干凈信號(hào)具有較好的還原度，保留了較多的細(xì)節(jié)。

3.2 實(shí)例信號(hào)分析

為了進(jìn)一步驗(yàn)證本文提出算法對(duì)實(shí)測(cè)脈搏信號(hào)的降噪效果，本文采用自行研制的便攜式脈搏信號(hào)采集裝置作為本文的脈搏信號(hào)樣本來(lái)源，并挑選了10名被試者參與脈搏信號(hào)的采集。由于脈搏信號(hào)受抖動(dòng)、呼吸等影響較為明顯，會(huì)產(chǎn)生一定的基線漂移，故本實(shí)驗(yàn)選擇安靜的環(huán)境，采用自主研制的便攜式脈搏采集裝置分別在同一時(shí)間段、同一房間內(nèi)進(jìn)行實(shí)驗(yàn)。測(cè)試者在平靜狀態(tài)下進(jìn)行實(shí)驗(yàn)，戴上采脈裝置，手臂保持水平放置，盡量讓橈動(dòng)脈與心臟處于同一水平面。本次實(shí)驗(yàn)

圖5 仿真信號(hào)去噪對(duì)比圖Fig.5 Comparison diagram of simulated signal denoising

采集1 min的脈搏信號(hào)，并取其中的1000個(gè)數(shù)據(jù)點(diǎn)作為實(shí)驗(yàn)對(duì)象，其余數(shù)據(jù)可持續(xù)驗(yàn)證實(shí)驗(yàn)。如圖6所示為本次采集的某位被測(cè)者的脈搏信號(hào)。為方便計(jì)算，將脈搏信號(hào)進(jìn)行歸一化處理。

表2 相關(guān)系數(shù)比較Tab.2 Comparison of correlation coefficients

同樣對(duì)脈搏信號(hào)進(jìn)行WOA-VMD分解以及噪聲分量評(píng)判，經(jīng)由VMD分解后產(chǎn)生11個(gè)模態(tài)分量，如圖7所示為分解效果圖，表3為實(shí)例信號(hào)的排列熵分布。

圖6 采集的原始數(shù)據(jù)Fig. 6 The original collected data

圖7 實(shí)例信號(hào)分解效果圖Fig.7 Decomposition effect diagram of sample signal

表3 各模態(tài)分量排列熵分布Tab.3 Arrangement entropy distribution of each IMF

從圖中可以看到，實(shí)例信號(hào)經(jīng)由VMD分解后提取到了信號(hào)的低頻噪聲，符合仿真信號(hào)得到的結(jié)論，然后根據(jù)排列熵分布對(duì)信號(hào)分量進(jìn)行處理。同樣為了能夠反映本文提出的算法對(duì)實(shí)例信號(hào)具有同樣的效果，將去噪效果與多個(gè)算法進(jìn)行對(duì)比，如圖8所示。

圖8實(shí)例信號(hào)去噪效果對(duì)比Fig.8 Comparison of sample signal denoising

圖8 （a）是本文所采用的方法得出的濾波效果，可以看到信號(hào)在各個(gè)波峰谷點(diǎn)處較為平滑，對(duì)于基線漂移的濾除較為優(yōu)異。圖8（b）中是小波閾值降噪法，可以看出該方法對(duì)于信號(hào)降噪后存在明顯的波形失真，整體波形顯得不平滑，例如第一個(gè)波峰處與原信號(hào)不符。圖8（c）為EEMD結(jié)合小波閾值去噪，可以看到該方法對(duì)于脈搏波的基線漂移濾除效果不佳，這與仿真信號(hào)的結(jié)論一致，并且在第二個(gè)脈搏周期中的重搏波特征不明顯（圖中圈出），存在一定的有效信號(hào)損失，這會(huì)影響到后續(xù)脈搏特征提取的研究。圖6（d）采用的是CEEMDAN結(jié)合小波閾值的去噪效果。可以看出該方法對(duì)脈搏信號(hào)的基線漂移進(jìn)行了一定的降噪處理，但降噪后的脈搏波依然存在一定的毛刺信號(hào)，在局部降噪表現(xiàn)上不佳。因此從圖像上看，本方法相較于其他方法具有更好的濾波效果。

為了量化說(shuō)明本算法的濾波效果，本文采用信噪比和均方根誤差進(jìn)行說(shuō)明。其中信噪比越大，均方根誤差越小，表明信號(hào)的降噪性能越為優(yōu)異。對(duì)于采集的10組數(shù)據(jù)，利用上述的3種算法計(jì)算出被試者信號(hào)的信噪比和均方根誤差，其中由于EEMD去噪存在明顯的基線漂移，因此并未對(duì)其進(jìn)行指標(biāo)計(jì)算。如圖9和圖10所示，為每一組脈搏實(shí)例信號(hào)所計(jì)算獲得的評(píng)價(jià)指標(biāo)圖。可以看出，經(jīng)由本文提出的降噪方法得到的評(píng)價(jià)指標(biāo)在10組實(shí)驗(yàn)中均優(yōu)于其他方法，因而可認(rèn)為本文的算法具有更好的降噪效果，對(duì)于特征信息的保存更為優(yōu)異。

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

圖10 不同方法的均方根誤差Fig.10 Root mean square error of different methods

4 總結(jié)

本文提出了一種利用鯨魚(yú)優(yōu)化算法尋找VMD分解過(guò)程中所需的參數(shù)，通過(guò)與排列熵確定有效分量，結(jié)合非均值濾波方法濾除脈搏信號(hào)中存在的噪聲。該方法利用鯨魚(yú)優(yōu)化的快速、自適應(yīng)等優(yōu)點(diǎn)尋找合適的模態(tài)分解數(shù)，有效解決了VMD分解中參數(shù)設(shè)置的問(wèn)題。通過(guò)排列熵分布，能夠有效地區(qū)分出噪聲分量與有效分量。最后通過(guò)非局部均值濾波，濾除噪聲信號(hào)。實(shí)驗(yàn)證明，相較于CEEMD、傳統(tǒng)VMD等，本方法在脈搏波降噪上的性能要更為出色，對(duì)原信號(hào)具有更高的還原度，為后期結(jié)合特征提取、脈象分類等工作奠定基礎(chǔ)，對(duì)脈診客觀化具有重要意義。然而VMD的雙參數(shù)優(yōu)化往往需要多次迭代，導(dǎo)致其運(yùn)行效率低下，因此如何提高算法的運(yùn)行效率還需要進(jìn)一步深入的研究。