振動信號包絡線的稀疏重構最優化算法研究與應用

于巖君， 葉慶衛， 陸志華， 周 宇

(寧波大學 信息科學與工程學院，浙江 寧波 315211)

HHT(Hilbert-Huang Transform)技術是由Huang等[1-2]提出的處理非線性、非平穩信號的新方法，這是繼傳統傅里葉時頻分析和小波分析方法后的具有重大突破性的應用數學方法。HHT需要經過經驗模態分解(Empirical Mode Decomposition，EMD)和Hilbert變換兩大步驟實現，其中EMD分解是其核心步驟。如今，EMD算法已在機械故障診斷、信號處理與分析、經濟金融工程等領域[3-5]得到了廣泛的應用，而包絡線的提取好壞直接影響EMD分解的結果，因此越來越多的學者對此問題做了大量的研究。鐘佑明[6]在經典三次樣條插值理論下提出Akima插值算法和分段冪函數法，通過三者之間的對比實驗，從擬合曲線的“光滑性”和“柔性”方面做了改進，但是仍無法消除端點飛翼和“過沖”、“欠沖”問題。廣義檢波濾波是除三次樣條插值法外的提取包絡的另一經典方法，此方法容易產生解調混頻效應，張帆等[7]提出采用合適的采樣頻率或提高采樣頻率可使廣義檢波解調分析中的混頻效應忽略不計，但是這種方法對于處理復雜信號效果不佳。張郁山等[8]提出了基于自回歸模型的延拓結構,采用“邊篩分邊延拓”方法將延拓得到的極值點和信號本身的極值點聯合起來構成包絡，最后篩分出IMF分量，但是自回歸模型需為預測的數據點個數設置一個較大的數值,影響篩選過程的精度,通用性比較差。

文獻[9]提出了基于稀疏復原的包絡線提取算法。稀疏復原本質是求解一個范數最小模型問題，它可以利用較少的觀測值精確地恢復出原始信號，某些信號的小波變換、離散余弦變換等都具有一定的稀疏性[10]，所以稀疏重構具有普遍適用性。雖然文獻[9]的稀疏復原方法在一定程度上抑制了端點效應，但是抗噪性能有待提高，而且在尋優改變DCT基頻帶寬度的變化因子m時，一方面需要人為挑選，降低了包絡線提取精度；另一方面尋找到的m只具有局部最優的性質，不具有全局性，所以尋優結果存在一定的片面性。

針對以上問題，本文結合粒子群算法(Particle Swarm Optimization，PSO)的獨特優勢[11]進一步對基于稀疏復原提取信號包絡線的算法做了改進。在文獻[9]的基礎上引入混合變異粒子群優化算法對變化因子m進行全局尋優，從而自適應地提取出最優的包絡線；并且將包絡線稀疏重構優化模型中的凹二次規劃問題轉換為凸二次規劃問題，利用內點法求得該模型中的稀疏解，其中，采用混合變異粒子群有效地避免了變化因子m陷入局部最優。實驗仿真與實際提取的斜拉索振動信號測試表明，基于PSO的信號包絡線稀疏重構最優化算法同樣可以克服端點飛翼問題，同時PSO的引入，可以自適應的找出最適合包絡線平穩特性的DCT基，克服了人為主觀確定的缺點，提高了包絡線的提取精度與抗噪性能。

1 稀疏復原擬合包絡線基本原理

包絡線的提取貫穿于整個EMD分解過程，后續也會影響到瞬時頻率的計算結果。但是至今在包絡線提取算法上還沒有系統的理論支撐，目前最常用的經驗算法有三次樣條插值、分段冪函數和廣義檢波濾波等。文獻[9]已經對三次樣條插值和包絡檢波原理做了詳細介紹，在此不再贅述。

文獻[9]基于壓縮感知理論，提出了利用稀疏重構的思想提取包絡線，下面以上包絡線為例，詳細說明此方法的基本原理。

用向量x(t)=[x1(t),x2(t),…,xN(t)]T表示所要求解的上包絡線信號，其中N表示采樣點的個數，用yd(t)=[xd1,xd2,…,xdk]表示從原始信號x(t)中提取的所有極大值點所構成的信號，其中K表示極大值點的個數。由于x(t)不是嚴格的K-稀疏，所以要將其投影在構造的基矩陣上進行稀疏表示，即：

x(t)=Ψs

(1)

s為稀疏基Ψ下的K稀疏信號，Ψ為N×N的余弦變換矩陣。稀疏基的構建是整個稀疏復原過程中的關鍵技術之一，根據包絡線的低頻特性和離散余弦變換的“能量集中”特性，經過離散余弦變換后信號的能量將集中在低頻部分，所以選擇了結構簡單的離散余弦變換基(DCT基)。為了提高包絡線的重構精度與效率，文獻[9]定義變化因子m，主要用于改變稀疏基的頻帶寬度，減小重構信號的頻率之差。其表達式如下

k=1,2,…,Nn=0,1,…,N-1

(2)

由于包絡線信號x(t)都經過其極值點，所以將極大值點所組成的信號yd(t)作為信號x(t)壓縮采樣后得到的M維觀測信號，其過程可描述為

yd(t)=Φ·x(t)

(3)

矩陣ΦM×N稱為觀測矩陣(又稱測量矩陣)，其中M在數值上與極大值點的個數K是相等的。由于所求稀疏解s稀疏度越高恢復出的x(t)效果越好，所以由式(1)、(3)可得包絡信號x(t)重構模型為

.t.yd(t)=Φ·Ψ(m)·s

(4)

基于信號包絡線的低頻特性，應該使恢復出的包絡線足夠光滑，從而根據頻域特征提取更精確的信息，基于此設計了一階微分矩陣D(N-1)，定義信號的平穩度函數E(t)來體現信號包絡的變化程度。其中：

(xN-1-xN)2

(5)

因此提取的最佳信號包絡線是滿足E(t)最小的稀疏重構獲得的曲線，綜上所述，基于稀疏復原的信號包絡線提取的數學模型為

(6)

2 基于粒子群的包絡線稀疏復原算法

2.1 稀疏復原擬合包絡線的改進模型

針對以上問題，本文對此做了兩點改進。首先將式(6)中的非凸模型轉化為凸二次規劃問題，從而利用內點法求得稀疏解；其次為了得到全局最優m值，本文將DCT基中的m放在分子上，范圍取[0,1]之間，利用混合變異粒子群算法對m自主尋優，避免了m陷入局部最優的問題。同時在尋優過程中，通過改變粒子群算法中的適應度函數，直接提取出誤差最小的包絡線信號。

針對式(6)中的l0范數最小化問題，匹配追蹤類算法是目前比較成熟的求解算法，包括MP、OMP、ROMP等。Donoho在文獻[12]提出滿足RIP條件下可將l0模型轉化為最小化l1模型，針對l1模型一般用凸松弛算法求解，包括BP、內點法等。因此式(6)中優化模型可改為

(7)

式(7)中的λ稱為正則化參數，此模型所求解一般是稀疏的(其證明詳見文獻[14])，由此符合稀疏復原的要求。由于上述l1-正則化模型是不可微的，所以受文獻[13]的啟示，通過引入線性不等式約束將稀疏重構模型轉化為一個凸二次規劃問題，這樣便可以用內點法求解此凸二次問題，重構模型修正為

(8)

若采用內點法求解上述模型中的稀疏解s，首先要構造懲罰函數，由此將不等式約束的優化問題轉換為無約束優化問題，其中懲罰函數構造形式如下

(9)

然后利用懲罰函數構造增廣目標函數，即無約束優化問題的目標函數

(10)

對于求解式(10)無約束優化問題本文選擇了具有超線性收斂速度且結構簡單的共軛梯度法，在每次迭代過程中根據凸二次函數f(s,u)的Hesse陣G產生的共軛方向以及由精確線搜索方法所確定的搜索步長來求解正定二次目標函數極小點，此極小點即為模型中的稀疏解s。

2.2 粒子群尋優變化因子m

本文利用粒子群算法對上述模型中的Ψ(m)進行尋優，將式(8)中的目標函數作為粒子群算法中的適應度函數，同時PSO的引入可以靈活改變稀疏復原方法，本文選擇了內點法求得稀疏解s， 最后將尋找到的最優稀疏基Ψ(m)和與之對應的稀疏解s代入x(t)=Ψ(m)·s中即可求得誤差最小的包絡信號。

用混合變異粒子群算法尋優變化因子m，能使搜索算法很快地收斂到全局最優值。由于只對m尋優，則搜索空間大小為1×1，且m的維數為一維信號。設種群粒子數為N,第i(i=1,2,…,N)個粒子的位置為xi，飛行速度為vi，第i個粒子目前為止搜索到的最優位置記為pbesti(簡記pbi)，整個種群中目前搜索到的最優位置記為gbesti(簡記gbi)，每個粒子根據以下公式更新自己的飛行速度和位置

vi(t+1)=w·vi(t)+c1·r1·(pbi(t)-xi(t))+

c2·r2·(gbi(t)-xi(t))

(11)

xi(t+1)=xi(t)+vi(t+1)

(12)

其中t為當前進化代數，w為慣性因子(一般取0～1)，c1,c2為學習因子(本文取c1=1,c2=2)，r1,r2為區間[0,1]之間的均勻隨機數。具體尋優步驟如下：

步驟1輸入所需要的參數，除了以上所提到的w,c1,c2,r1,r2,N還要輸入最大迭代次數maxD；

步驟2種群初始化，不同的粒子代表不同的變化因子m，且以隨機的方式求出每個粒子的初始位置與初始速度；

步驟3構造適應度函數，每一個m值對應一個DCT基，調用內點法(其求解過程在2.3節中詳細說明)子程序求出與DCT基對應的稀疏解s，由式(1)即可求出與之對應的包絡線x(t)，進而由式(5)得到不同的m所對應的平穩度值，用以下公式作為判斷每個粒子(即m)好壞的適應度函數

(13)

步驟4更新粒子速度，根據上面式子(13)中的適應度值找出局部最優值pbesti與全局最優值gbesti，然后根據式(11)、(12)更新每個粒子m的速度與位置；

步驟5按照文獻[14]進行粒子變異，得到變異后的粒子群；

步驟6檢驗終止條件，當達到最大迭代次數或者全局最優位置滿足最小界限時，則終止迭代；否則轉到步驟3。

2.3 基于粒子群的稀疏復原算法總體流程

利用混合變異粒子群算法在每次迭代過程中都需要由當前m值所構造的DCT基和觀測值yd(t)來重構包絡信號x(t)，所以稀疏復原方法的選擇對整體的算法流程至關重要。本文選擇了內點法，下面詳細介紹怎樣利用內點法解決式(10)中的模型。

步驟1參數設置，輸入初始點s0,u0，懲罰因子τ0>0，衰減系數ρ∈(0,1)，終止誤差0≤ε≤1，迭代次數為k=1；

步驟2根據式(9)構造懲罰函數P(s,u)，以sk-1為初始點利用共軛梯度法求解式(10)中的無約束子問題f(s,u)，求得極小值點sk；

步驟3判斷終止條件，若τk·p(s,u)≤ε，停算，輸出s*≈sk作為稀疏解；

步驟4令τk+1=ρ·τk，k=k+1，返回步驟1。

綜上所述，基于粒子群的包絡線稀疏重構最優化算法總體流程如圖1。

圖1 粒子群算法的總體流程圖

粒子群算法在每次迭代過程中通過調用上述內點法子程序對求得的m尋優，直到迭代終止時找到最優的變化因子m*,進而得到最適合包絡線變化特性的DCT(m*)基，最后利用內點法復原出原始包絡信號。內點法的引入提高了稀疏復原的速度，可用于復雜信號包絡的提取，同時利用粒子群算法解決了局部尋優變化因子m的缺陷，較大程度上提高了包絡線提取精度。

3 實驗仿真

為了驗證稀疏復原方法提取包絡線的有效性，本文使用的仿真信號為

x(t)=1.8e-3tcos(2πf1t+1.3)+

0.6e-3tcos(2πf2t+0.7)+

0.1e-2tcos(2πf3t+0.7)

(14)

仿真信號中的模態頻率分別取f1=0.8 Hz、f2=3 Hz、f3=10 Hz，時間t的取值范圍為[0,3]s，信號采樣點數為500個。為了突出本文改進算法的優越性，結合粒子群尋優算法對式(14)中的信號采用內點法復原包絡信號，與三次樣條插值算法和分段冪函數算法對比，其結果如圖2，同時與文獻[9]的結果(見圖3)相對比。

圖2 基于粒子群的內點法復原方法提取信號上包絡

由實驗結果可知，基于粒子群的內點法復原信號上包絡線明顯地抑制了端點效應問題，通過與分段冪函數法結果對比，明顯的解決了端點處過沖與欠沖問題，而且自適應地找到改變DCT基頻帶寬度的變化因子m，從而自適應地尋找到一條符合包絡線緩慢變化特性的最佳曲線。與文獻[9]的結果相對比，本文提取包絡更加平滑，其中dd=0.02，而文獻[9]中dd=0.04，而且文獻[9]中通過m=m+1人為篩選取得m=18時得到最優包絡，而本文利用混合變異粒子群算法自適應地在[0,1]之間尋找全局最優m值(m=4.707 426×10-2)取其倒數大概是m=23，不僅減小了人為確定所帶來的誤差，解決了局部尋優的片面性，而且提高了提取精度。為了更好地觀察粒子群尋優過程，圖4給出了每次迭代過程中適應度函數所對應的值。

圖3 基于稀疏復原方法(OMP)提取信號上包絡線

圖4 內點法復原包絡時粒子群算法尋優m過程

從圖4中可以看出，在0～10次迭代過程中，尋優速度較快，在10～25次時，開始逐漸靠近最優值，25～40次時，尋優值基本上保持不變。所以，利用粒子群對變化因子m尋優不僅可以精確地找到m值以提高包絡線的提取精度，而且可以快速地找到變化因子m以提高包絡線的提取速度。同時，從統計出的粒子群尋優程序運行的平均時間(見表1)來分析本文算法的運算復雜度。本數據是通過大量的實驗仿真運行統計得到的，其運行環境如下：MATLAB 2014b，Intel Core i3處理器，3.3 GHz的CPU處理器，4 GB內存。

從表格中可以看到，當粒子數目一定時，隨著迭代次數的增加，程序運行的平均時間逐漸減少，從而說明達到一定的迭代次數時可以找到最優的粒子，當超過這個迭代次數時，粒子繼續保持最優值。通過大量實驗，我們得到這個迭代次數大概為30，所以無需過多迭代次數即可得到全局最優值；還可以看到，隨著粒子數目的增多，程序運行的平均時間逐漸增加，但是通過仿真數據觀察到，當粒子數目達到一定數量時，提取的包絡誤差保持不變，從而說明并不是粒子數量越多越好，通過大量實驗得到最佳粒子數目為20，所以本文中的仿真取粒子數N=20。

表1 粒子群尋優算法的時間統計

為了檢驗改進算法對信號處理的抗噪性能，仍然采用式(14)仿真信號x(t)，然后給x(t)添加不同信噪比的高斯白噪聲n(t),得到加噪信號x1(t)，然后對x1(t)用本文提出的基于粒子群的稀疏優化算法提取包絡線。圖5中給出的是信噪比取15 dB時的提取結果,從圖中可以看出，本文方法可以自適應地提取出最佳包絡線，消除了端點效應問題，而且除端點外提取的結果幾乎和三次樣條插值、分段冪函數算法相吻合。與文獻[9]中取同樣信噪比的結果(見圖6)相比，本文提出的方法不僅具有一定的抗噪性，而且提高了包絡線的精度，明顯減小了提取誤差(dd=0.13)，同時對比圖6中的m=2,本文對m進行了全局尋優。

圖5 信噪比15 dB時基于粒子群的內點法復原上包絡

Fig.5 The interior-point methods based on PSO for extraction the upper envelope of the signal added noise when SNR is 15 dB

當噪聲強度增強時，選取信噪比為10 dB的加噪信號再次進行仿真，由圖7、圖8可以看出，本文方法仍具有一定的抗噪性能，不僅可以解決端點效應問題，而且誤差仍可以控制在一定的范圍之內，其中dd=0.22，而文獻[9]中的dd=1.41。隨著信噪比的降低，本文提取信號包絡線的精度沒有明顯降低，相比文獻[9]有了更強的抗噪性能。

圖6 信噪比15 dB時基于稀疏復原(OMP)提取信號上包絡

Fig.6 The OMP based on sparse recovery for extraction the upper envelope of the signal added noise when SNR is 15 dB

圖7 信噪比10 dB時基于粒子群的內點法復原上包絡

圖8 信噪比10 dB時基于稀疏復原(OMP)復原上包絡

4 實際斜拉索振動信號的應用測試

在上述實驗仿真中，為驗證本文算法提取的包絡線的抗噪性能，人為添加了高斯白噪聲，但是自然界中噪聲模型很復雜，為了驗證包絡線提取在EMD算法中的應用，對實際采樣信號進行了仿真與EMD分解。本文的實測信號采集于寧波市某一斜拉索大橋，此橋總長67 m，由102根直徑0.15 m拉索構成支撐系統；采用WS-ZHT2振動設備和雙傳感器，并且將雙傳感器裝置在拉索和梁端的鉸支上[15-16]，如圖9所示。

圖9 傳感器安裝圖

為了驗證本文提出的改進后的基于稀疏復原提取包絡線方法的有效性，本實驗仿真中采樣頻率fs=100 Hz采樣512個點。圖10是采用本文改進算法對采集信號的包絡線提取結果，并且與三次樣條插值結果對比，由圖看出，基于粒子群的稀疏優化方法仍可以精確地提取采集信號的包絡線，并且數值保持在允許誤差之內，再次驗證了本文提出方法的可行性和實用性。

圖10 采用基于粒子群的稀疏優化和三次樣條插值、分段冪 函數法對實際采集信號提取上包絡的比較

Fig.10 The comparison chart on extracted the envelope of actual acquisition signal based on sparse optimization and based on cubic spline interpolation method and subsection power function method

EMD分解將非平穩、非線性信號分解為有限個本征模函數IMF，這是Hilbert Huang變換(HHT)的關鍵步驟。圖12是采用基于粒子群的稀疏優化方法提取采集信號包絡線的EMD分解結果，其中IMFi表示本征模函數的階數，隨著階數的增加，頻率依次降低。與文獻[9]中的結果(見圖11)對比可以看出基于粒子群的稀疏復原方法經EMD分解得到的IMF分量在端點處更加精確，圖形表現的更加平穩，更接近于正弦函數，這一點在低頻分量上表現的更為明顯。

圖11 基于稀疏復原算法對實測信號EMD分解結果

圖12 基于粒子群的稀疏復原對實測信號EMD分解結果

5 結 論

EMD分解的關鍵步驟是信號上下包絡線的提取，本文將包絡提取模型中的凹問題轉換為凸二次規劃問題，利用內點法求得模型中的稀疏解，并且將粒子群優化算法與稀疏復原算法相結合，利用混合變異粒子群算法，將改變DCT基頻帶寬度的變化因子m作為粒子屬性，通過多次迭代獲得全局最優m，最終自適應地提取出最優的包絡線信號。經驗證該方法可以有效解決端點效應問題，提高包絡線的提取精度與抗噪性能，粒子群自適應尋優方法的引入，減少了人為選擇所產生的提取誤差，從實際橋梁振動信號的模態分解過程中，證明了本文方法的有效性和實用性。

[1] HUANG N E, SHEN Z, LONG S R, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J]. Proceedings of the Royal Society A Mathematical Physical & Engineering Sciences, 1998, 454(1971):903-995.

[2] DELECHELLE E, LEMOINE J, NIANG O. Empirical mode decomposition: An analytical approach for sifting process[J]. IEEE Signal Processing Letters, 2005, 12(11):764-767.

[3] CHENG J, YU D, YU Y. The application of energy operator demodulation approach based on EMD in machinery fault diagnosis[J]. Mechanical Systems & Signal Processing, 2007, 21(2):668-677.

[4] FLEUREAU J, NUNES J C, KACHENOURA A, et al. Turning tangent empirical mode decomposition: a framework for mono-and multivariate signals[J]. IEEE Transactions on Signal Processing, 2011, 59(3):1309-1316.

[5] 楊云飛. 基于EMD分解技術的不同市場原油價格相關性分析及預測研究[D].武漢：華中科技大學,2011.

[6] 鐘佑明. 希爾伯特-黃變換局瞬信號分析理論的研究[D]. 重慶：重慶大學,2002..

[7] 張帆, 丁康. 廣義檢波解調分析的三種算法及其局限性研究[J]. 振動工程學報, 2002, 15(2):243-248.

ZHANG Fan，DING Kang.Research on the three algorithms and limitations of generalized detection-filtering demodulation analysis[J].Journal of Vibration Enginering，2002, 15(2):243-248.

[8] 張郁山,梁建文,胡聿賢. 應用自回歸模型處理EMD方法中的邊界問題[J]. 自然科學進展, 2003,13(10):48-53.

ZHANG Yushan，LIANG Jianwen，HU Yuxian. The processing of end effects in EMD method by autoregressive model.[J]. Progresses in Nature Science，2003,13(10):48-53.

[9] 徐靜妹，葉慶衛，王曉東，等. 斜拉索振動信號的包絡線稀疏復原算法[J]. 振動與沖擊，2015，34(23)：187-193.

XU Jingmei，YE Qingwei，WANG Xiaodong，et al. Envelope sparse recovery algorithm for stay cable vibration signals[J]. Journal of Vibration and Shock，2015，34(23)：187-193.

[10] CANDES E J, TAO T. Near-optimal signal recovery from random projections: universal encoding strategies?[J]. IEEE Transactions on Information Theory, 2007, 52(12):5406-5425.

[11] BOERINGER D W, WERNER D H. Particle swarm optimization versus genetic algorithms for phased array synthesis[J]. IEEE Transactions on Antennas & Propagation, 2004, 52(3):771-779.

[12] MALLAT S G, ZHANG Z. Matching pursuits with time-frequency dictionaries[J]. IEEE Transactions on Signal Processing, 1993, 41(12):3397-3415.

[13] KIM S J, KOH K, LUSTIG M,et al. An interior-point method for large-scalel1—regularized least squares[J]. IEEE Journal of Selected Topics in Signal Processing, 2007, 8(3):1519-1555.

[14] 陳君波, 葉慶衛, 周宇,等. 一種新的混合變異粒子群算法[J]. 計算機工程與應用, 2007, 43(7):59-61.

CHEN Junbo，YE Qingwei，ZHOU Yu，et al. New multi-mutation particle swarm optimization algorithm[J]. Computer Engineering and Applications, 2007, 43(7): 59-61.

[15] YE Q, FENG Z, YUAN D. Nonlinear vibration signal tracking of large offshore bridge stayed cable based on particle filter[J]. Polish Maritime Research, 2016, 22(4):70-75.

[16] YE Q, FENG Z, ZHAN Q. Vibration signal fundamental frequency estimation of bridge stayed cable based on era algorithm.[J].Journal of the Balkan Tribological Association 2016, 22(2):1335-1343.