優化變分參數與改進小波軟閾值重構濾波算法

汪 梅，王 將，李遠成，董立紅，馬 天，李銘禹

(1.西安科技大學 計算機科學與技術學院，陜西 西安 710054；2.西安科技大學 電氣與控制工程學院，陜西 西安 710054)

0 引 言

礦工的不安全情緒可能會導致礦難的發生，識別礦工情緒能夠判斷其當前的情緒狀態，為煤礦的安全生產減少人為因素的安全威脅[1]。腦電信號以其客觀性和不易隱藏性的特點，被廣泛應用于情緒識別領域[2-3]。但腦電信號由于受到在采集過程中外界環境干擾和采集設備等因素的影響，而被引入噪聲導致信號質量下降[4]。這些噪聲有可能影響信號的主要特征，對于后續的特征分析和情緒識別有一定影響[5-7]。因此，在利用腦電信號識別礦工情緒的過程中，對于噪聲的濾除就顯得尤為必要。

經典腦電濾波方法主要包括傅里葉分解或小波分解及重構。傅里葉分解由于其構造函數為周期性的正、余弦波，導致其對非周期性或局部特征較明顯的信號處理效果較差[8]。小波分解具有良好的時頻分析能力，可以很好地分辨信號的突變部分[9-11]。近年來，經驗模態分解及其改進方法在腦電濾波算法研究中日漸增多[12-14]。相比于小波分解，經驗模態分解不需要預先設定母小波和分解層次就可以自適應地對非線性、非平穩信號進行分解處理，但此方法缺乏良好的數學理論，并且對采樣數據和噪聲都很敏感[15-18]。為此，變分模態分解(variation mode decomposition，VMD)以其完備的數學理論支持和較好的噪聲魯棒性，已廣泛應用在生物電信號濾波領域。

KAUR等將VMD算法分別與離散小波變換和小波包變換結合起來對腦電信號進行濾波，發現基于小波包變換的VMD法性能更優[19]。DORA等利用VMD提取腦電信號中的眼電分量，并結合回歸的方法獲得干凈的腦電信號[20]。XIAO等通過結合VMD和小波閾值的方法對肌電信號進行濾波，實驗濾波效果要優于單一的小波閾值法和經驗模態分解法[21]。盧莉蓉等則利用VMD和小波軟閾值方法來去除心電信號中肌電的干擾[22]。但是VMD的模態數和帶寬分別取決于分解個數K和懲罰因子α的預設值，而且K和α值的大小都會影響濾波效果。

針對上述VMD在對腦電信號濾波時，VMD的分解效果會受到參數K和α選取的影響，導致信號濾波效果差的問題，提出一種優化變分參數與改進小波軟閾值重構濾波算法。該算法結合了VMD能提供有效頻率劃分和小波閾值時頻分析能力強的優點，改進了由K和α選取不當造成的影響和傳統小波閾值存在逼近程度較差或平滑性不足問題。文中創新點包括以下4個方面：第1，給出烏燕鷗參數優化的VMD算法；第2，擴展相關系數差值比的方法來區分有效分量和含噪分量；第3，利用改進的小波軟閾值處理含噪分量得到去噪分量；第4，提出優化變分參數與改進小波軟閾值重構的濾波算法。最后，對模擬信號和情緒腦電信號進行濾波處理，檢驗所提方法的濾波性能。

1 烏燕鷗參數優化變分模態分解

1.1 變分模態分解算法

VMD是一種自適應的信號處理方法，與經驗模態分解法相比，VMD提供了更有效的頻率劃分，可以更好地避免模態混疊等問題[23-24]。非平穩多分量的信號f(t)通過VMD可以得到K個具有特定中心頻率ωk的模態分量uk。為了確定uk和ωk，需要構造一個約束變分問題，見式(1)。

(1)

式中 {uk}，{ωk}分別為模態分量集合和中心頻率集合。

引入拉格朗日乘數λ求解上式，可得到

(2)

(3)

(4)

(5)

式中n為迭代次數；α為懲罰因子；τ為保真系數；ε為收斂精度判據。

1.2 烏燕鷗優化算法

1.2.1 遷移行為

遷移過程中烏燕鷗應該滿足以下3個條件。

條件1：避免碰撞

通過引入參數SA來計算烏燕鷗個體的新位置C(z)，以避免相鄰個體之間的碰撞

C(z)=SA×P(z)

(6)

式中C(z)為不與其他烏燕鷗發生碰撞的新位置；P(z)為烏燕鷗個體當前位置；z為當前迭代次數。SA為避免碰撞的變量參數，計算公式如下

SA=Cf-(z×(Cf/Maxiterations)

(7)

z=0,1,2,…,Maxiterations

(8)

式中Cf為控制變量，使得SA的值從Cf線性遞減到0，設置Cf值為2，Maxiterations為最大迭代次數。

條件2：確定相對間距

在避免碰撞之后，當前個體會向最優個體靠攏，以獲得更好的位置

M(z)=CB×(P*(z)-P(z))

(9)

式中M(z)為當前個體位置P(z)與最優個體位置P*(z)的相對間距。CB是一個隨機變量，能夠使個體更好地進行全局搜索。CB的計算公式為

CB=0.5×Rand

(10)

式中Rand是介于[0，1]的隨機數。

條件3：向最優個體靠攏

確定無碰撞位置C(z)和相對間距M(z)后，下面開始向最優個體位置靠攏，以到達新位置D(z)。

D(z)=C(z)+M(z)

(11)

式中D(z)為烏燕鷗遷移行為后的位置。

1.2.2 攻擊行為

在烏燕鷗定位獵物后，它們會改變自身的飛行速度和角度對獵物進行攻擊，從而在空中產生螺旋狀運動軌跡。其在x′，y′和z′平面的運動行為的數學模型為

x′=Radius×sin(i)

(12)

y′=Radius×cos(i)

(13)

z′=Radius×i

(14)

Radius=u×eiv

(15)

式中Radius為螺旋運動的半徑；i為在[0，2π]范圍內的螺旋角；e為自然對數的底數；u和v為常量，設置u和v的值均為1。

烏燕鷗攻擊行為后位置的計算公式為

P(z+1)=(D(z)×(x′+y′+z′))×P*(z)

(16)

式中P(z+1)為烏燕鷗攻擊行為后的位置。

1.3 烏燕鷗參數優化變分模態分解

在VMD中，參數K和α的選取將會影響最終的分解效果。若K值過大，會造成分解過度，導致相鄰模態中心頻率的間距較近；而如果K值太小，會造成分解不足，導致分解失去了實際意義[25]。類似的，如果α選取不恰當，也會出現模態混疊的問題。因此，合理選擇K和α的值對于信號分解結果的準確性尤為重要。

由于烏燕鷗優化算法(sooty tern optimization algorithm，STOA)具有尋優能力強、精度高等特點[26]。利用STOA優化VMD以確定最佳參數組合[K，α]，并將文獻[27]所提的包絡熵作為適應度函數，將VMD參數的優化過程轉化為利用STOA尋求最小包絡熵值的過程。通過VMD將原始信號分解為K個模態分量，如果模態分量中包含的噪聲分量越多，其與原始信號相關的特征信息就越不明顯，則包絡熵越大；反之，包絡熵越小。圖1為烏燕鷗參數優化變分模態分解(STOA-VMD)算法流程。

圖1 STOA-VMD算法流程Fig.1 Flow of STOA-VMD algorithm

STOA-VMD算法的具體步驟如下。

步驟1：初始化STOA參數，并隨機形成一個數量為N的烏燕鷗種群，設置迭代次數為z，尋優維數為2，則烏燕鷗個體位置可表示為：P(z)={P1(z)，P2(z)}，其中，P1(z)=K，P2(z)=α。

步驟2：以當前烏燕鷗個體位置對應的K和α作為VMD的輸入參數對信號進行分解。

步驟3：計算不同位置的烏燕鷗個體對應的包絡熵值，更新得到當前最小包絡熵值。

步驟4：判斷是否達到迭代終止條件，如果達到，尋優停止，輸出最優的K和α，否則，令z=z+1，并更新種群位置，返回步驟2繼續進行迭代優化。

2 模態判斷和改進小波軟閾值算法

2.1 相關系數差值比的模態判斷

通過STOA-VMD將原始信號分解為K個模態分量，需要從中判斷并重構合適的分量來對信號進行濾波。相關系數可以度量原始信號與模態分量之間的相似性，相關系數值越大則相似性越強，反之則越弱。文中分別將與原始信號相關系數較大和較小的模態分量稱為有效分量和含噪分量。文獻[28]利用固定的相關系數值來區分這2類分量，難以對不同信號作出靈活的調整，具有一定的局限性。

文中提出一種相關系數差值比的判斷條件，根據相鄰模態與原始信號相關系數的差值比來確定有效分量和含噪分量的臨界點，判斷條件如下

(17)

(18)

式中Rm為原始信號f(t)與第m個分量um(t)之間的相關系數，E[·]和D[·]分別為數學期望和方差。

根據臨界點判斷條件找出有效分量和含噪分量之間的臨界點um(t)。由于模態分量按照頻率從低到高排列，而有效分量的頻帶大多分布在低頻部分。因此，對臨界點um(t)之前的有效分量保留并利用所提改進小波閾值處理其余含噪分量。

2.2 改進小波軟閾值的分解重構算法

根據小波閾值濾波的基本原理可知[29]，小波基、分解尺度的確定和閾值函數的構造都會影響其濾波效果。軟閾值和硬閾值是2種常用的閾值函數，但軟閾值法會造成重構信號與原信號之間逼近程度較差的缺點；硬閾值法的不連續則會引起重構信號平滑性不足的問題。

(19)

為了直觀地反映形狀調節因子β的作用效果，選取T=1,β=0.5,2,5,10,20,30，比較改進閾值函數與軟、硬閾值函數的特性，如圖2所示。通過調整β的值，改進閾值函數可以在軟、硬閾值之間變動，這也使得它在實際應用上更為靈活。隨著|wz|的增加，改進閾值函數越來越接近于硬閾值函數，這就彌補了軟閾值法的缺點。同時，改進閾值函數具有高階可導性，可以克服振蕩，提高信號的平滑度。

圖2 閾值函數特性比較Fig.2 Comparison of threshold function characteristics

3 優化變分與改進軟閾值重構算法

在上述理論基礎上，提出優化變分與改進軟閾值重構算法(STOA-VMD-IWS)，如圖3所示。

圖3 STOA-VMD-IWS算法流程Fig.3 Flow of STOA-VMD-IWS algorithm

首先，通過STOA-VMD算法找到最佳參數組合[K，α]，利用優化的VMD分解原始信號得到一組模態分量{u1,u2,…,uK}。然后，利用相關系數差值比的判斷條件來區分有效分量和含噪分量，并用IWS算法對含噪分量進行濾波。最后，重構有效分量和去噪分量，實現信號濾波。

4 實驗與結果分析

4.1 參數優化變分模態實驗與結果分析

實驗選取公開腦電數據集DEAP中一段預處理后的腦電信號(electroencephalogram，EEG)作為研究對象，采樣率為128 Hz，采樣個數為600。在其中分別加入強度為-10～10 dB且間隔為5 dB的高斯白噪聲，得到5種含噪EEG信號。

下面以含有10 dB高斯白噪聲的EEG信號為例，利用STOA-VMD算法尋找最佳參數組合[K，α]。設置烏燕鷗種群數為30，最大迭代次數Maxiterations為10，K的范圍為[2，10]，且K為整數，α的范圍為[500，4 000]。圖4為STOA-VMD算法的尋優過程。

圖4 STOA-VMD算法的尋優過程Fig.4 Optimization process of STOA-VMD algorithm

從圖4可以看出，當迭代數為5時，最小包絡熵值為0.969 2，此時[K，α]值為[10，3 348]。因此，以K=10，α=3 348作為VMD的輸入參數對EEG信號進行分解，結果如圖5所示。

圖5 STOA-VMD算法的分解結果 Fig.5 Decomposition results of STOA-VMD algorithm

可以看出，u1～u10的中心頻率之間相互獨立，且并未發生模態混疊的現象，證明了STOA-VMD算法的分解結果較好。為了驗證STOA參數尋優的準確性，將原始EEG信號頻譜中對應的主頻率值與u1～u10的中心頻率值進行比較，結果見表1。

表1 EEG的主頻率和u1～u10的中心頻率比較Table 1 Comparison of the main frequency of EEG and the center frequency of u1～u10

從表1可知，u1～u10的中心頻率值基本對應于EEG信號的主頻率值，能很好地反映EEG信號的頻率特性。此外，如圖5所示，u1～u10之間沒有混疊，證明STOA-VMD算法參數優化的有效性。

4.2 相關系數差值比實驗與結果分析

在得到模態分量u1～u10后，計算了u1～u10和原始EEG信號之間的相關系數，并得到相鄰模態分量相關系數差的絕對值見表2。

表2 相鄰模態分量相關系數差的絕對值Table 2 Absolute value of the correlation coefficient difference between adjacent modal components

利用公式(17)的模態判斷條件計算相關系數差值比，可以得到：|R6-R5|/|R5-R4|=0.076 5，|R5-R4|/|R4-R3|=9.136 7，則u5為有效EEG分量和含噪EEG分量的臨界點。

4.3 改進小波閾值濾波實驗與結果分析

根據相關系數差值比的模態判斷條件，得到有效EEG分量和含噪EEG分量的臨界點為u5，將臨界點前的模態分量作為有效EEG分量保留，對其余含噪EEG分量利用IWS算法處理后與有效EEG分量重構得到濾波后EEG信號，如圖6所示。

圖6 信號濾波效果Fig.6 Filtering effects of the signal

經所提STOA-VMD-IWS算法濾波后噪聲得到了明顯地抑制。為了驗證所提算法的優勢，文中將其與傳統VMD算法、STOA-VMD和小波硬閾值濾波算法(STOA-VMD-HWTF)、STOA-VMD和小波軟閾值濾波算法(STOA-VMD-SWTF)進行了對比實驗，并以信噪比(signal-to-noise ratio，SNR)和均方根誤差(root mean square error，RMSE)作為評價指標，對比結果見表3。

由表3可知，傳統VMD算法的SNR較低，RMSE較高，其濾波效果較差。在傳統VMD算法的基礎上，盡管STOA-VMD-HWTF算法和STOA-VMD-SWTF算法的濾波效果得到了一定的改善，但所提STOA-VMD-IWS算法在SNR和RMSE方面效果更好，腦電濾波效果更佳。

表3 EEG信號濾波評價指標對比Table 3 Comparison of EEG signals filtering evaluation indexes

4.4 EEG信號STOA-VMD-IWS濾波實驗分析

為了進一步驗證所提濾波算法對實測EEG信號的有效性，根據情緒二維模型理論，實驗選取煤礦井下工作環境相關的3類情緒圖片，即積極、消極和中性，在E-Prime 2.0系統上對被試進行情緒誘發，并利用Neuroscan公司的NuAmps設備完成對被試EEG信號的采集。

下面以在消極情緒圖片誘發下得到的一段EEG信號為例，圖7為STOA-VMD-IWS算法對該EEG信號的濾波效果。

圖7 4種EEG信號的濾波效果Fig.7 Filtering effects of 4 EEG signals

其中，圖7(a)的EEG信號分別來自FP1,C3,O1和T4電極，可以看出EEG信號具有大量的毛刺，說明其含有大量的高頻噪聲。圖7(b)是利用所提STOA-VMD-IWS算法對這4個電極處EEG信號進行濾波處理后的信號。經所提STOA-VMD-IWS算法處理后，高頻噪聲得到了有效抑制，也保留了EEG信號的大部分特征，可以得到其幅值變化及波動情況。

由于實驗無法預先得到“干凈”的EEG信號，導致無法利用SNR和RMSE來評價所提算法的濾波效果，故文中利用噪聲抑制比(noise suppression ratio，NSR)來評價所提算法和傳統VMD算法對4種EEG信號的濾波效果，計算結果見表4。NSR越大，說明濾波效果越好，其計算公式如下

表4 4種EEG信號的NSR計算結果Table 4 NSR calculation results of 4 EEG signals

(20)

從表4可知，2種濾波算法對C3和O1的EEG信號濾波效果要優于其他2種EEG信號，且所提算法對4種EEG信號的濾波效果均優于傳統VMD算法，證明所提算法在腦電濾波方面的有效性。

5 結 論

1)給出變分模態分解參數的烏燕鷗優化方法，即利用烏燕鷗優化算法對變分模態分量個數和懲罰因子進行優化，得到參數優化的VMD算法，解決了VMD有效分解問題。

2)擴展相關系數差值比判別方法，用于判別VMD分解后的有效分量和含噪分量，以便對含噪分量做進一步處理，解決了模態分量選取問題。

3)對小波軟閾值函數進行改進，用于對含噪分量的小波分解與重構，實現對含噪分量的去噪，解決了硬閾值不連續引起重構信號平滑性不足和軟閾值重構信號的逼近度較差問題。

4)提出優化變分參數與改進小波軟閾值重構的濾波算法，對模擬信號和情緒腦電信號進行濾波處理，信噪比最大提高了4.568 3 dB，平均提高了3.284 7 dB；均方根誤差最大降低了0.169 1，平均降低了0.069 5，解決了傳統VMD算法對腦電濾波效果較差問題，有助于提高腦電的信噪比和礦工情緒識別的準確率，對減少因礦工不安全情緒導致的礦井安全事故的發生具有重要意義。