閆志勛，劉怡明2，林 童，孫曉云，王明明

(1. 石家莊鐵道大學(xué) 電氣與電子工程學(xué)院，河北 石家莊 050043；2. 蘇州大學(xué) 機(jī)電工程學(xué)院，江蘇 蘇州 215137)

作為一種新興的超聲導(dǎo)波檢測(cè)手段，磁致伸縮導(dǎo)波技術(shù)在錨桿和錨固質(zhì)量無(wú)損檢測(cè)領(lǐng)域中有檢測(cè)快速、可靠性高、檢測(cè)距離長(zhǎng)等優(yōu)點(diǎn)，在各工程結(jié)構(gòu)的健康檢測(cè)和監(jiān)測(cè)中有良好的應(yīng)用前景[1-4]。由于電磁超聲信號(hào)幅值較微弱，常被噪聲所淹沒，使得采集到的信號(hào)失真[5]，因此，在磁致伸縮錨桿檢測(cè)中，對(duì)信號(hào)進(jìn)行降噪十分必要。

在磁致伸縮錨桿檢測(cè)中，檢測(cè)信號(hào)因噪聲的疊加而被視為非穩(wěn)態(tài)信號(hào)。針對(duì)此類信號(hào)，所采用的降噪方法主要有經(jīng)驗(yàn)?zāi)B(tài)分解(EMD)、小波閾值法等，但是小波閾值法的小波基和閾值選取較困難，EMD方法存在模態(tài)混疊與端點(diǎn)效應(yīng)等問題，并且抗噪性較差。文獻(xiàn)[6]中提出變分模態(tài)分解(VMD)方法，該方法相對(duì)于遞歸式模態(tài)分解方法有完善的數(shù)學(xué)理論基礎(chǔ)，能更準(zhǔn)確地分解信號(hào)，運(yùn)算效率高，抗噪性好，并且不易發(fā)生模態(tài)混疊和端點(diǎn)效應(yīng)。該方法需要預(yù)先設(shè)置懲罰因子和模態(tài)個(gè)數(shù)，分解結(jié)果在較大程度上受這2個(gè)參數(shù)的影響，然而磁致伸縮錨桿檢測(cè)信號(hào)是復(fù)雜的，模態(tài)個(gè)數(shù)和懲罰因子難以確定。

為了實(shí)現(xiàn)信號(hào)的最優(yōu)分解，針對(duì)模態(tài)個(gè)數(shù)和懲罰因子難以確定的問題，本文中利用布谷鳥搜索(CS)算法對(duì)VMD的最優(yōu)參數(shù)組合進(jìn)行搜尋[7]，從而確定VMD的最優(yōu)參數(shù)，實(shí)現(xiàn)信號(hào)的最優(yōu)分解。單獨(dú)使用VMD算法對(duì)信號(hào)進(jìn)行處理可能導(dǎo)致不能完全剝離反射波信號(hào)中的回波信息。獨(dú)立分量分析(ICA)是一種基于樣本高階統(tǒng)計(jì)信息的新的信號(hào)處理方法[8]，目前的研究大都是基于多觀測(cè)通道的，而通常采集到的錨桿磁致伸縮信號(hào)均為一維數(shù)據(jù)，不能滿足ICA對(duì)于多維數(shù)據(jù)的要求。

針對(duì)以上問題，本文中提出CS算法優(yōu)化VMD與ICA結(jié)合的方法，以達(dá)到更好的降噪效果。

1 CS算法優(yōu)化VMD參數(shù)

1.1 CS算法與VMD基本理論

1.1.1 CS算法

CS算法是Yang等[7]基于布谷鳥的巢寄生特性及萊維飛行而提出的，首先假設(shè)以下3條理想化規(guī)則：

1)每只布谷鳥一次只產(chǎn)一個(gè)卵，孵化時(shí)的鳥巢是隨機(jī)選擇；

2)最好的巢保留到下一代；

3)寄生巢數(shù)量是固定的，寄生宿主發(fā)現(xiàn)外來(lái)鳥蛋的概率也是固定的，設(shè)為Pa∈[0,1]。

基于以上理想化規(guī)則，布谷鳥對(duì)鳥巢的位置和路徑更新是按照萊維飛行模式進(jìn)行的，路徑和位置的更新公式為

(1)

1.1.2 VMD算法

VMD算法是對(duì)變分問題的求解，分解步驟如下。

1)構(gòu)建目標(biāo)函數(shù)。

①進(jìn)行Hilbert變換，尋求每個(gè)模態(tài)函數(shù)uk(t)解析信號(hào)，并得到單邊頻譜

(2)

式中：δ(t)為沖擊函數(shù)；k為分解的本征模態(tài)分量(IMF)的層數(shù)。

②通過加入指數(shù)項(xiàng)，估計(jì)每個(gè)模態(tài)的中心頻率e-jωkt，其中ωk為中心頻率，把每個(gè)模態(tài)的頻譜調(diào)制到相應(yīng)的基頻帶，即

(3)

③估計(jì)每個(gè)模態(tài)信號(hào)帶寬，從而獲得由變分問題組成的目標(biāo)函數(shù)，即

min{uk},{ωk}∑Kk=1??tδ(t)+jπt[]uk(t){}e-jωkt22{},∑Kk=1uk=f(t),ì?í????(4)

2)引入拉格朗日算子λ(t)和二次懲罰因子α，得到增廣拉格朗日函數(shù)

L({uk},{ωk},λ)=

(5)

式中：〈·〉為內(nèi)積；K為當(dāng)前分層數(shù)。

(6)

式中arg min表示式子達(dá)到最小值時(shí)的變量的取值。

利用Parseval-Plancherel傅里葉等距變換后，將頻率ω替代為ω-ωk，并將得到的結(jié)果轉(zhuǎn)換為負(fù)頻率區(qū)間積分的形式，則優(yōu)化問題的解為

(7)

根據(jù)同樣的過程，求得中心頻率更新方式

(8)

(9)

綜上所述，VMD的實(shí)現(xiàn)過程如下。

2)根據(jù)式(7)—(9)，在頻域內(nèi)更新uk、ωk、λ；

1.1.3 峭度準(zhǔn)則

峭度是描述波形尖峰度的一個(gè)無(wú)量綱參數(shù)，對(duì)沖擊信號(hào)特別敏感。峭度值C的定義為

(10)

式中：μ、σ分別為信號(hào)x的均值和標(biāo)準(zhǔn)差[9]。當(dāng)峭度值C=3時(shí)，定義分布曲線具有正常峰度，即零峭度，如果信號(hào)中疊加了大量沖擊脈沖信號(hào)，則信號(hào)的峭度值會(huì)明顯增大。利用VMD算法對(duì)錨桿的磁致伸縮信號(hào)進(jìn)行分解后，如果IMF中更多的是反射波信息，波形中出現(xiàn)規(guī)律性沖擊脈沖，則分量信號(hào)峭度值較大。

1.2 CS-VMD參數(shù)優(yōu)化

現(xiàn)實(shí)中磁致伸縮信號(hào)是復(fù)雜的，在進(jìn)行VMD運(yùn)算時(shí)，核心問題是如何選取合適的K及α。如果僅改變一個(gè)參數(shù)，另一個(gè)參數(shù)設(shè)置為常數(shù)，這樣得到的只是局部最優(yōu)結(jié)果[10]。為了得到一組最佳參數(shù)，本文中利用CS算法對(duì)VMD的輸入?yún)?shù)K和α進(jìn)行尋優(yōu)，在進(jìn)行參數(shù)優(yōu)化時(shí)，選用峭度值作為適應(yīng)度函數(shù)。

在VMD運(yùn)算后，計(jì)算所有IMF的峭度值，將峭度值的最大值稱為局部極大值，用maxEIMF表示，與該值所對(duì)應(yīng)的分量稱為局部最佳分量。

VMD算法的α和K的優(yōu)化步驟如下：

1)初始化鳥巢個(gè)數(shù)、更新概率、迭代次數(shù)等參數(shù)，初始化群體，隨機(jī)產(chǎn)生J個(gè)鳥巢的初始位置yi(i=1,2,…,J)，確定適應(yīng)度函數(shù)為峭度值；

2)在不同參數(shù)條件下對(duì)信號(hào)進(jìn)行VMD運(yùn)算，計(jì)算每組參數(shù)相應(yīng)的適應(yīng)度值maxEIMF并記錄當(dāng)前的最優(yōu)解；

3)保留上一代最佳IMF的鳥巢，利用式(1)對(duì)鳥巢位置進(jìn)行更新；

4)比較新的鳥巢適應(yīng)度值與上一代鳥巢適應(yīng)度值，選取更好的鳥巢位置并保留；

5)比較隨機(jī)數(shù)r和發(fā)現(xiàn)概率Pa，若r>Pa，則更新鳥巢位置，反之保留位置不變；

6)若未滿足終止條件，則返回步驟2)；

7)輸出全局最優(yōu)位置。

算法流程如圖1所示。

圖1 布谷鳥搜索(CS)算法優(yōu)化變分模態(tài)分解(VMD)算法流程圖

2 VMD與ICA聯(lián)合降噪方法

2.1 ICA

ICA是在源信號(hào)的先驗(yàn)信息和信道傳輸參數(shù)均未知的情況下，利用隨機(jī)矩陣對(duì)信號(hào)進(jìn)行混合，對(duì)混合后的信號(hào)進(jìn)行線性變換，使得變換后的信號(hào)盡量逼近源信號(hào)。設(shè)Y(t)=(Y1(t),Y2(t),…,YM(t))T為M個(gè)觀測(cè)信號(hào)，設(shè)源信號(hào)S(t)=(S1(t),S2(t),…,SN(t))T為N個(gè)未知獨(dú)立源信號(hào)，各分量之間互相獨(dú)立，觀測(cè)到的線性混合模型為

Y(t)=AS(t)，

(11)

式中A為滿秩的M×N型混合矩陣(M≥N)。ICA算法的目的是找到分離矩陣W，使得

(12)

在混合矩陣A及源信號(hào)S都未知的情況下，估計(jì)出的各獨(dú)立分量的排列順序及幅值具有不確定性，但是不影響信號(hào)整體的判斷。為了提高計(jì)算的收斂速率，通常采用快速ICA(FastICA)算法。

2.2 排列熵算法

排列熵(permutation entropy,PE)是由Bandt等[11]提出的一種衡量序列復(fù)雜度和動(dòng)力學(xué)突變的方法。計(jì)算步驟如下。

給定時(shí)間序列{X(i),i=1,2,…,n}進(jìn)行相空間重構(gòu)，可得到重構(gòu)矩陣[12]為

j=1,2,…,Q,

(13)

式中：m為嵌入維數(shù)；T為延遲時(shí)間；Q為重構(gòu)分量個(gè)數(shù)。對(duì)X(i)中第j個(gè)向量升序排列，得

x[i+(j1-1)T]≤x[i+(j2-1)T]≤…≤

x[i+(jm-1)T]。

(14)

如果重構(gòu)向量存在相等的情況，則按原順序排列，因此，對(duì)于任意一個(gè)時(shí)間序列X(j)，符號(hào)序列為S(l)=(j1,j2,…,jm)其中l(wèi)=1,2,…,q，且q≤m!。每種符號(hào)序列出現(xiàn)概率記為P1,P2,…,Pq，此時(shí)，時(shí)間序列X(i)的q種不同符號(hào)序列的PE可定義為

(15)

當(dāng)Pj=1/m!時(shí)，Hp(m)最大，對(duì)應(yīng)值為ln(m!)。對(duì)Hp(m)進(jìn)行歸一化處理，即

0≤Hp=Hp/ln(m!)≤1。

(16)

Hp值表示時(shí)間序列的復(fù)雜度，其值越小，信號(hào)越規(guī)則，反之，隨機(jī)性越強(qiáng)。

2.3 VMD-ICA聯(lián)合降噪

本文中通過CS-VMD運(yùn)算，得到若干IMF，選取保留較多振動(dòng)信號(hào)的分量重構(gòu)成一組信號(hào)，利用FastICA算法對(duì)其與源信號(hào)組成的2路輸入通道進(jìn)行解混，得到降噪后的信號(hào)。對(duì)于由噪聲引起的毛刺，可以進(jìn)一步進(jìn)行小波閾值處理。

選取IMF重構(gòu)信號(hào)時(shí)，選用PE作為性能指標(biāo)，PE對(duì)突變的信號(hào)有有效的檢測(cè)效果[13]，因此，根據(jù)PE值進(jìn)行信號(hào)的重構(gòu)。

綜上所述，提出基于CS-VMD參數(shù)優(yōu)化與ICA聯(lián)合降噪方法，實(shí)現(xiàn)步驟如下：

1)利用CS算法優(yōu)化VMD參數(shù);

2)利用VMD方法分解回波信號(hào);

3)依據(jù)PE值對(duì)分解的各模態(tài)進(jìn)行重構(gòu)，對(duì)單通道信號(hào)進(jìn)行升維;

4)利用FastICA算法解混升維后的信號(hào);

5)對(duì)解混后的信號(hào)進(jìn)行小波閾值處理，濾除由噪聲引起的毛刺。

算法流程如圖2所示。

3 裸錨桿長(zhǎng)度實(shí)驗(yàn)分析

測(cè)試中采用的錨桿為長(zhǎng)度為3 m，直徑為22 mm。錨桿磁致伸縮信號(hào)采集示意圖如圖3所示。檢測(cè)方法如下：通過激勵(lì)線圈產(chǎn)生磁致伸縮導(dǎo)波，當(dāng)導(dǎo)波到達(dá)端部時(shí)會(huì)發(fā)生反射，對(duì)應(yīng)的反射波會(huì)有相位改變。由安放在錨桿另外一端的感應(yīng)線圈接收反射回來(lái)的信號(hào)，再通過采集模塊傳送至上位機(jī)。假設(shè)錨桿的長(zhǎng)度為L(zhǎng)0，感應(yīng)線圈接收到端面反射信號(hào)的時(shí)間為t，超聲導(dǎo)波在錨桿中的傳播波速為c，通常取值為5 000 m/s。計(jì)算錨桿長(zhǎng)度的公式為

(17)

CS—布谷鳥搜索；IMF—本征模態(tài)分量；FastICA—快速獨(dú)立分量分析。圖2 變分模態(tài)分解-獨(dú)立分量分析(VMD-ICA)聯(lián)合降噪算法流程圖

為了驗(yàn)證本文中提出方法的降噪效果，以試驗(yàn)研究中實(shí)際測(cè)得的裸錨桿磁致伸縮信號(hào)為對(duì)象進(jìn)行降噪結(jié)果對(duì)比。裸錨桿原始信號(hào)波形如圖4所示。

圖4 裸錨桿原始信號(hào)波形

由圖可知，信號(hào)中存在大量噪聲信號(hào)，端面反射信號(hào)淹沒在噪聲中，很難判斷反射時(shí)刻，這使得錨桿的計(jì)算長(zhǎng)度與實(shí)際長(zhǎng)度相差很大。

3.1 CS算法優(yōu)化VMD參數(shù)分析與驗(yàn)證

為了驗(yàn)證CS算法與傳統(tǒng)頻率觀察法相比在優(yōu)化VMD參數(shù)方面的優(yōu)勢(shì)，選取試驗(yàn)研究中的實(shí)際測(cè)得的裸錨桿磁致伸縮信號(hào)進(jìn)行以下分析：對(duì)實(shí)測(cè)信號(hào)選取不同K值進(jìn)行VMD運(yùn)算，懲罰因子α使用默認(rèn)值，即α=2 000，觀察每個(gè)模態(tài)分量的中心頻率，如表1所示。由表可知，當(dāng)K=8時(shí)，出現(xiàn)了模態(tài)分量中心頻率相近的現(xiàn)象(過分解現(xiàn)象)，因此K的取值為7。

表1 不同當(dāng)前分層數(shù)K值時(shí)各模態(tài)分量的中心頻率

利用CS算法優(yōu)化VMD參數(shù)組合[K,α]。CS算法參數(shù)設(shè)置如下：更新概率Pa=0.25，鳥巢數(shù)n=20，分解層數(shù)K的搜索范圍為(2,15)，懲罰因子α的搜索范圍為(100,3 000)，最大迭代次數(shù)為50。尋優(yōu)結(jié)果為K=7，α=1 525，分解結(jié)果如圖5所示。

計(jì)算傳統(tǒng)頻率觀察法與本文中方法VMD運(yùn)算后每層IMF的PE值，如表2所示。

將傳統(tǒng)頻率觀察法與優(yōu)化后得到的參數(shù)組合應(yīng)用到本文中提出的VMD-ICA降噪方法進(jìn)行對(duì)比。選取PE值較大的IMF構(gòu)成虛擬通道，因此傳統(tǒng)頻率觀察法取IMF3—IMF6，本文中方法取IMF2—IMF6，分別重構(gòu)信號(hào)。2種方法的降噪效果對(duì)比如圖6所示，錨桿長(zhǎng)度的計(jì)算結(jié)果如表3所示。由于ICA沒有關(guān)于原始信號(hào)的任何先驗(yàn)知識(shí)，幅值具有不確定性，因此，降噪后的信號(hào)幅值可能會(huì)有改變，但對(duì)信號(hào)分析沒有影響。由表3可知，2種方法都能有效地降低噪聲影響，但是本文中提出的方法能更大限度地恢復(fù)原始信號(hào)，計(jì)算錨桿長(zhǎng)度誤差較小。

綜上所述，針對(duì)模態(tài)個(gè)數(shù)K和懲罰因子α很難確定的問題，CS算法能有效地搜尋到一組最佳輸入?yún)?shù)。

圖5 模態(tài)分量波形及頻譜

表2 各層本征模態(tài)分量(IMF)的排列熵(PE)值

(a)傳統(tǒng)頻率觀察法(b)本文中提出的方法圖6 傳統(tǒng)頻率觀察法與本文中提出方法的降噪效果

表3 傳統(tǒng)頻率觀察法與本文中提出方法的錨桿長(zhǎng)度計(jì)算結(jié)果

3.2 錨桿長(zhǎng)度實(shí)驗(yàn)分析

為了驗(yàn)證本文中提出的方法相較于其他降噪方法的優(yōu)勢(shì)，分別采用小波閾值去噪、EMD-ICA聯(lián)合降噪加小波閾值處理[13]、VMD-ICA降噪加小波閾值處理[14]及本文中提出的方法進(jìn)行裸錨桿降噪處理，結(jié)果如圖7所示。原始信號(hào)和上述4種降噪方法得到的錨桿長(zhǎng)度值和相對(duì)誤差率如表4所示。

(a)小波閾值去噪法(b)文獻(xiàn)[13]中方法(c)文獻(xiàn)[14]中方法(d)本文中提出的方法圖7 裸錨桿降噪效果對(duì)比

表4 錨桿長(zhǎng)度計(jì)算結(jié)果

文獻(xiàn)[13]中各IMF與源信號(hào)的互相關(guān)系數(shù)如表5所示，文獻(xiàn)[14]中參數(shù)組合為K=7，α=2 000，各層IMF的峭度值如表6所示。本文中提出方法的參數(shù)在3.1節(jié)中已給出，此處不再贅述。

由圖7和表4可知，小波閾值去噪法和文獻(xiàn)[13]中方法降噪效果不明顯，波形不平滑，含噪聲較多，錨桿的計(jì)算長(zhǎng)度與實(shí)際長(zhǎng)度誤差較大；文獻(xiàn)[14]中方法波形平滑，但是誤差率最大；而本文中提出的降噪算法的降噪效果好，含噪聲較少且波形平滑，最大限度地恢復(fù)了原始信號(hào)，計(jì)算錨桿長(zhǎng)度誤差最小，僅為0.83%。

表5 各層本征模態(tài)分量(IMF)與原信號(hào)的互相關(guān)系數(shù)

表6 各層本征模態(tài)分量(IMF)的峭度

4 錨桿和錨固實(shí)驗(yàn)分析

測(cè)試中采用的錨桿長(zhǎng)度為2 m，直徑為22 mm，錨固長(zhǎng)度為1 m，錨固層材質(zhì)為混凝土，脫錨段為沙土。錨固錨桿信號(hào)采集示意圖如圖8所示。假設(shè)錨桿、錨固長(zhǎng)度分別為L(zhǎng)1、L2，超聲導(dǎo)波在錨桿、錨固中的傳播波速分別為c1、c2，感應(yīng)線圈接收到錨固底端面反射信號(hào)的時(shí)間為t1，裸錨桿波速c1一般取5 000 m/s，經(jīng)驗(yàn)證，錨固中傳播速度c2為3 000 m/s。計(jì)算錨固長(zhǎng)度的公式為

(18)

圖8 錨桿和錨固信號(hào)采集示意圖

錨桿和錨固去除電磁脈沖后的原始信號(hào)的波形如圖9所示。由圖可知，信號(hào)被噪聲所污染，缺陷回波和底端回波被噪聲淹沒。VMD參數(shù)優(yōu)化參數(shù)組合為分解層數(shù)K=9、懲罰因子α=1 424、每層IMF的PE值，如表7所示。

圖9 原始波形

表7 各層本征模態(tài)分量(IMF)的排列熵(PE)值

分別應(yīng)用4種方法得到的降噪后波形，如圖10所示。原始信號(hào)和4種降噪方法所得錨固長(zhǎng)度值和相對(duì)誤差率如表8所示。

(a)小波閾值去噪法

(b)文獻(xiàn)[13]中方法

(c)文獻(xiàn)[14]中方法

(d)本文中提出的方法圖10 采用不同方法的錨桿和錨固系統(tǒng)降噪效果

表8 錨固長(zhǎng)度計(jì)算結(jié)果

由圖10和表8可知，VMD-ICA降噪方法降噪效果不明顯，缺陷回波和底端回波不明顯，小波閾值和EMD-ICA去噪波形不平滑，缺陷回波信號(hào)和底端面回波有失真現(xiàn)象，而本文中提出的降噪算法的降噪效果好，含噪聲較少，能明顯看出缺陷回波及底端回波信號(hào)，計(jì)算錨固長(zhǎng)度相對(duì)誤差最小，僅為0.5%。

5 結(jié)論

本文中采取CS算法、VMD與ICA相結(jié)合的方式，利用CS算法對(duì)VMD的最優(yōu)參數(shù)組合進(jìn)行搜尋，借助VMD方法較好的噪聲穩(wěn)健性和自適應(yīng)分解的特性，以及ICA通過信號(hào)的不同特征，能夠分離出未知原始信號(hào)能力，更大限度地去除噪聲。通過實(shí)測(cè)信號(hào)結(jié)果分析來(lái)驗(yàn)證本文中所提方法的有效性。本文中提出的方法應(yīng)用在錨桿磁致伸縮導(dǎo)波檢測(cè)上，能更有效地保留有用的反射波信號(hào)。結(jié)果表明，相對(duì)于傳統(tǒng)小波閾值去噪法、VMD-ICA降噪以及EMD-ICA降噪方法，本文中提出方法的降噪效果有較大改善，計(jì)算錨桿長(zhǎng)度的誤差更?。辉阱^固的質(zhì)量檢測(cè)中能更有效地提取缺陷信號(hào)及回波信號(hào)。