優化變分模態分解的超聲多普勒測流信號誤差模型研究

趙軍華，戴聰聰，李 叢，馮 陽，鄧 權，張清波

（深圳市宏電技術股份有限公司，感知研發中心，廣東 深圳 518000）

0 引 言

近年來，隨著科技的飛速發展，農田灌溉、水流監測、家用水表、工廠用水監測等不同領域都廣泛使用了超聲多普勒流量計，其測量精度和誤差也直接影響產品的工作效率[1]。在實際應用環境中，由于液體密度的不均勻性，發射的超聲波產生折射和衍射現象，導致回波信號中會夾雜著各類不同的噪聲干擾，尤其是受白天陽光影響直接導致信號完全被噪聲淹沒[2]，因此測量液體流量的超聲回波信號大多是低信噪比、非線性、非平穩信號[3,4]。如何提取回波信號中的有效信號成為實現超聲流量計高精度、高可靠性的關鍵。目前常用的去噪方法主要有小波閾值去噪[5,6]、經驗模態分解(EMD)[7,8]，從多尺度將信號分解成不同頻率的分量、具有良好時頻局部化特性的小波閾值去噪由于小波基的選取問題，不具備自適應分解能力，雖然EMD 通過添加白噪聲解決了小波閾值去噪自適應問題，能自適應分解多個本征模態函數分量(IMF)，但仍存在模態混合、缺乏嚴格的理論公式、算法效率低、抗噪性能差等諸多問題。因此，EMD 并不能精確地分離出有用信號[9]。

2014 年，Xie[10]提出了一種新的時頻分析法—變分模態分解(VMD)，該方法以強大的數學基礎來解析信號瞬時頻率，具有實際的物理意義，通過迭代搜索模型最優解來實現信號精準分離的同時利用自身的Wiener 濾波特性可以獲得不錯的去噪效果，克服了EMD 的模態混疊問題，所以VMD 也被廣泛使用在非平穩信號處理中[11-13]。然而，VMD 算法受其懲罰因子和分解層數參數影響，設置不當極易導致過分解問題，還需進一步優化[14,15]。海鷗算法(SOA)具有較強的全局搜索能力，它可以適用于參數尋優[16]，相對于其他優化算法如麻雀、蝙蝠等具有較強的魯棒性，但海鷗算法在更新種群位置時極易被束縛，導致在后續迭代過程中陷入局部最優，造成結果的多樣性[17,18]。

因此利用海鷗算法結合VMD 對超聲多普勒信號進行降噪處理，引入柯西變異算子(CV)，設置隨機擾動解決海鷗算法迭代過程中陷入局部最優的問題，優化VMD 參數以包絡譜值為自適應度函數，將云相似度計算作為有效分量篩選標準。并針對VMD 出現的中低頻振蕩現象引入奇異頻譜分析(SSA)，進一步提高回波信號信噪比。

1 算法原理

1.1 VMD原理

VMD 算法可將任何信號從低頻到高頻進行有效分離。在這個算法中，各模態分量及其對應中心頻率的確定是一個約束問題，可引入二次懲罰因子和拉格朗日乘數，使問題不受約束，在頻域內不斷更新模態分量、中心頻率和拉格朗日乘數，直到獲得變分約束問題的最優解為止。

1.2 CVSOA優化VMD參數

海鷗算法因具有超強的全局搜索能力和易于操作實現的特點，使得它可適用于參數尋優。其基本原理是計算“海鷗的最佳運動位置”作為問題解，算法主要有兩種行為：

（1）遷徙行為。海鷗個體更新探索新位置，個體之間避免碰撞且與最優個體保持密切聯系。

式中：Cs表示海鷗種群所有個體相互分散的位置；Ps表示單個海鷗所處位置；x表示當前迭代次數；A表示在一定空間內海鷗的遷徙行為；fc表示常數可調整海鷗位置。

更新完海鷗位置后，海鷗個體保證不會碰撞且會根據最優海鷗個體方向移動。

式中：Ms表示海鷗個體所處位置Ps向記錄的當前最優值位置Pbs移動路徑；B的取值是隨機的，它可以適當地平衡搜索和擴展行為；rd為一個隨機數，范圍為[0,1]。

同時，海鷗個體根據收斂方向向新位置進行移動。

式中：Ds表示海鷗個體新位置。

（2）攻擊行為。海鷗進攻行為呈螺旋形下降。其行為軌跡x軸、y軸、z軸分量分別如下式所示。

式中：Ps(x)為存儲歷史最佳海鷗位置，位置更新時，算子易被Pbs束縛，易被固定在某個區域，從而影響了結果的多樣性，即陷入局部最優。

通常情況下，設置隨機可解決這個問題。Cauchy 分布的隨機序列兩翼分布較寬，零點波峰高于Guassian 分布，比Guassian分布的抗干擾能力更強。

則可利用柯西變異來改進SOA算法。計算公式為：

式中：Cauchy(0,1)為分布函數；rand(0,1)是一個均勻分布隨機數；p為隨機變異概率。

在迭代后期，進入收斂狀態所有海鷗個體攻擊到同一位置不再更新。因此，結合柯西變異策略，對其更新的位置設置擾動，逃離早熟狀態。實現步驟如下：

（1）陷入局部最優后，判斷隨機數是否小于變異概率p。

（2）對最佳海鷗個體位置Pbs進行柯西變異，并保留Pbs。

（3）判斷海鷗種群是否繼續陷入局部最優，若是，則循環步驟（1）～（3），否則找到全局最優。

利用下式對Pbs進行變異策略。

式中：rand 為一個大于0且小于1的均勻分布隨機數；Pbs為當前海鷗種群最佳個體位置。

1.3 CVSOA-VMD適應度函數

根據CVSOA 算法理念，需要依據適應度函數來評判參數組合是否為最優結果。文獻[19]提出了包絡熵的概念，包絡熵可以更好地反映信號的稀疏特性和不確定程度，數值越大，信號的不確定性也就越大。本文選用包絡熵作為CVSOA 的適應度函數，x(j)(j= 1,2,…,N)信號計算公式如下：

式中：N為信號采樣點數；a(j)為信號x(j)經Hilbert解調后得到的包絡信號。

1.4 融合VMD與SSA的降噪模型

云相似度法是基于云模型理論提出的一種衡量兩向量的相似方法。將CVSOA 用于VMD 參數優化，以包絡譜熵值作為CVSOA 的適應度函數，將信號按最優參數進行VMD 分解，得到一系列IMFs，計算其云相似度值，將閾值以下判定為雜波分量并去除，最后利用SSA 去除殘余噪聲。具體步驟流程如圖1所示。

圖1 融合VMD與SSA的降噪算法流程圖Fig.1 Flow chart of noise reduction algorithm integrating VMD and SSA

（1）通過逆向云發生器[20]得到原始信號y(t)的云向量=(Ex,En,He)。

（2）通過CVSOA-VMD 分解原始回波信號y(t)，得到IMFk(k= 1,2…,n)。

（4）求取原始信號與各云向量的余弦夾角作為云相似度值，即：

（5）采用統計分析法確定閾值大小，將低于閾值的IMF分量去除，得到有效IMF分量。

（6）對有效分量進行重構，進行SSA 二次諧波處理后得到最終去噪信號。

2 算例分析

為了驗證該方法的可行性和有效性，分別采用以下3種方案處理模擬信號和實測走車數據。方案一：通過EMD 分解信號得到IMF 組分，利用云相似度對各組分進行相關性分析，重構有效組分，進行SSA 二次濾波處理；方案二：利用SOA對VMD 參數進行優化，按最優參數分解得到IMF 成分，并將干擾成分篩選掉重構有效成分，進行SSA 二次濾波處理；方案三：利用CV 對SOA 進行優化，根據最優參數分解得到的IMF分量，計算云相似值篩選出噪聲分量，再對重構信號進行SSA二次濾波，得到最終的降噪信號。便于表述和分析，方案一、二、三分別記為EMD-SSA、SOA-VMD-SSA 和CVSOAVMD-SSA。

2.1 仿真信號分析

由于實測數據中含有一定的誤差，因此方法的有效性和可行性首先通過模擬信號進行驗證。根據超聲波時間序列的特點，根據公式(11)生成模擬信號[21]，時間為2 ms，頻率為500 kHz，加入高斯白噪聲，其信噪比為8，其含噪波形如圖2所示。

圖2 時域波形圖Fig.2 Time domain waveform

式中：A、a、f、τ分別為信號的幅值、衰減系數、特征頻率和峰值時刻。

采用VMD 對含噪仿真信號進行分析，選取包絡譜熵值作為CVSOA 的適應度函數，利用CVSOA 對懲罰因子α和模態分解個數K值進行尋優?？紤]到分解層數K設置過大而引起的過分解現象，設定K值運算最大值為9。圖3 是海鷗種群進化迭代變化的對比曲線圖，圖3 中為兩種方法的全局最優包絡譜值。

圖3 迭代次數與包絡譜熵值關系圖Fig.3 Relation between iteration times and envelope spectrum entropy

由圖3 可以看出：CVSOA 算法優化VMD 在第5 次迭代收斂，包絡譜熵值為1.942，最佳參數組合為[4 365,6]，SOA算法優化VMD 在第6 次迭代收斂，包絡譜熵值為1.913，最佳參數組合為[8 327,6]，分別利用最優參數對仿真信號進行VMD 分解，如圖4所示。

圖4 EMD、VMD模擬信號分量圖Fig.4 EMD and VMD analog signal component diagrams

從圖4 中的對比分析可以看出：VMD 分解模態從低頻到高頻分布，低頻振蕩也主要集中在前幾階模態中。EMD 分解模態從高頻到低頻排列，信號能量主要集中在高階模態中，多個分量出現嚴重的模態混疊現象，分解出多個原本屬于同一分量部分信息的低頻分量，這是由于EMD 終止條件的不合理，導致了過分解；SOA-VMD 分解的IMF2、IMF3 分量產生頻率重疊部分。相比之下，CVSOA-VMD 分解的各分量頻率單一，模態混疊現象的抑制效果更好，優于EMD 和SOAVMD 方法。以求準確篩選信號和噪聲分量，分別計算各分量的云相似度值，圖5 為3 種分解算法各分量云相似度值，圖5中3種方法都取IMF1～IMF6分量。

圖5 3種分解算法各分量云相似度值Fig.5 Cloud similarity values of each component of three decomposition Algorithms

由圖5分析可知：階數分量靠前的云相似度值更大，云相似度值越大的分量與原始信號波形相似度也越大，可以看出CVSOA-VMD 分解的IMF3分量與無噪原始信號最接近，SOAVMD 分解的IMF2、IMF3 分量與無噪原始信號最接近。3 種分解算法云相似度值都呈現先增大后減小的趨勢，將IMF 分量小于閾值0.4 視為噪聲成分。對有效成分進行重構，采用SSA進行二次諧波降噪處理，上述3種方案降噪處理后的結果如圖6所示。

圖6 3種算法降噪信號圖Fig.6 Three algorithms for denoising signals

對比圖6可以看出：CVSOA-VMD-SSA 有效抑制了諧波干擾，幅值及趨勢變化得到了完整地保留，相比EMD-SSA、SOA-VMD-SSA 更平滑，與無噪原始信號存在差異最小，變化規律趨勢也更接近。加入不同信噪比的白噪聲，以均方根誤差和信噪比等兩個指標進行評價，定量分析和說明上述3種方案對模擬信號的降噪效果，計算結果如圖7所示。

圖7 信噪比及均方根誤差Fig.7 The signal-to-noise ratio and root mean square error

由圖7可知：隨著白噪聲信噪比的增加，去噪信號的信噪比逐漸增大，均方根誤差逐漸減小，CVSOA-VMD-SSA 較SOA-VMD-SSA、EMD-SSA 均有改善，信噪比可達30.78 dB，均方根誤差最低達到0.01，說明與原信號的差異最小，去噪效果更好，更能反映原信號的特征信息。

3 走車實驗分析

3.1 水槽走車實驗

水槽走車試驗是以水槽中的水保持靜止、走車帶動流速儀在水槽中以設定的勻速運動為原則設計的，實驗裝置如圖8 所示。走車速度可通過PLC 上位機操作設定在全范圍內的任意恒定速度值，走車裝置最大額定速度可達3 m/s，水槽長度為30 m，寬度為1 m，深度為1.2 m。此次走車實驗速度設置為0.3 m/s，走車距離為30 m，流量計緊扣固定在走車裝置上，將流量計的外接線通過485 轉232 串口連接到電腦上，串口助手進行數據接收。首先對串口接收的數據進行去趨勢分析，圖9 為走車實測的超聲回波時域波形圖。

圖8 實驗裝置圖Fig.8 Experimental device

圖9 實測原始信號Fig.9 Original signal

3.2 信號降噪處理

接下來對圖9 采樣的時域信號利用CVSOA 對VMD 參數進行尋優，最優參數組合分別為[3 894,6]、[7 652,6]，按最優參數組合對信號進行VMD 分解，計算各分量的云相似度值，結果如表1所示。

表1 不同分量云相似度值Tab.1 Cloud similarity values of different components

從表1 中可以看出，IMF1～IMF5 分量的相似度值差異明顯，差值范圍在0.32～0.52 左右。為進一步觀察云相似度值的分布規律和門檻選擇的可靠性及普遍性，取超聲多普勒信號10 組，分別進行CVSOA-VMD 分解，再計算云相似度值，限于篇幅取前5 組。如表2 所示，是每組信號的最優參數組合。

表2 各樣本最優參數組合Tab.2 Optimal parameter combination of each sample

由表2 可知，第6 組信號分解模態個數為5，故曲線擬合均取前5個分量的云相似度值，如圖10所示。

圖10 云相似度值擬合曲線圖Fig.10 Cloud similarity value fitting curve

由圖10 可知：IMF3 成分雖有一定波動，但各值均較大，說明與原信號最有關聯，含原始信息量最多，能完全反映原信號的特征信息，而IMF1、IMF2、IMF4、IMF5 的云相似度值均較低，表明與原始信號沒有關聯。根據多次實驗統計分析確定云相似度閾值為0.6。因此，選擇以上包含最多原始信號特征的IMF成分進行重構，并對重構信號進行SSA二次諧波干擾抑制，如圖11所示，得到最終降噪信號。

通過對比圖11(a)、圖11(b)分析可知，直接小波閾值去噪信號更光滑，但是信號出現失真較為嚴重，EMD-SSA 去噪不徹底，還殘留少部分噪聲，而圖11(c)、圖11(d)基本去除了存在的噪聲，但SOA-VMD-SSA 去噪波峰波谷值不能很好保留，相比之下CVSOA-VMD-SSA 能夠很好的保留波峰波谷值且噪聲抑制效果更好。計算其信噪比和均方根誤差，以進一步量化CVSOA-VMD-SSA 算法的去噪效果，結果如表3所示。

由表3 分析可知，CVSOA-VMD-SSA 信噪比達21.71 dB，相比直接小波閾值、EMD-SSA、SOA-VMD-SSA 分別提高13.35、8.90、5.47 dB；均方根誤差最低達0.010 4，表明CVSOA-VMD-SSA 方法的去噪效果更優。為更進一步驗證去噪算法對流速探測精度的影響，表4給出了文中所述的幾種算法在探測流速在0.1～0.6 m/s范圍內的測速誤差。

表4 不同走車流速探測精度對比 m/sTab.4 Comparison of the detection accuracy of different vehicle flow rate

由表4分析可知：CVSOA-VMD-SSA 相比于其他算法更利于隨機噪聲的抑制，隨著流速的增加，幾種算法誤差也逐漸增大，但對比不同走車實測流速下的探測精度，CVSOAVMD-SSA 誤差最小，范圍在0.01～0.03m/s，為實際工程應用提供理論支撐。

4 結 語

本文利用CVSOA 對VMD 的參數進行優化，提出了將CVSOA-VMD 與SSA結合的方法，通過對模擬仿真信號和不同流速走車試驗時間序列的降噪分析，驗證了該方法的有效性和可靠性，并與WD、EMD-SSA、SOA-VMD 等方法對比分析，得出以下結論：

（1）模擬信號結果顯示：CVSOA-VMD-SSA 能有效地去除噪聲干擾，克服模態混疊及SOA 易陷入局部最優解問題，降噪信號的幅值及變化趨勢得到完整地保留，與無噪原始信號存在差異最小，從評價指標看，CVSOA-VMD-SSA 明顯優于其他幾種算法。

（2）走車試驗數據表明：采用多組信號統計分析可確定云相似度閾值為0.6，選用有效分量重構，并進行SSA 二次濾波降噪，CVSOA-VMD-SSA 能有效抑制噪聲，對比不同走車流速探測精度，CVSOA-VMD-SSA 誤差最小，范圍在0.01～0.03 m/s，可為實際工程應用提供理論支撐。