基于自適應低秩稀疏分解的超聲缺陷回波檢測方法*

周航銳孫 堅*徐紅偉繆存堅

(1.中國計量大學機電工程學院，浙江 杭州 310018；2.浙江省特種設備檢驗研究院，浙江杭州 310012)

超聲技術在無損檢測中得到了廣泛的應用，但由于電子噪聲和相干噪聲的存在(材料的晶粒散射回波和噪聲影響存在)，在實際工程使用超聲對金屬材料中的微小缺陷檢測時仍具有挑戰(zhàn)性。國內外學者提出了許多從含噪超聲信號中檢測缺陷回波的方法，如小波變換[1-2]、經(jīng)驗模態(tài)分解[3]、稀疏表示[4]等。

近年來，低秩稀疏分解(Low-Rank Sparse Decomposition，LRSD)作為一種新的矩陣分解理論被提出并廣泛用于視頻前背景分離[5]，圖像去噪[6]、語音增強[7]和故障診斷[8]等領域。

LRSD算法在充分考慮數(shù)據(jù)矩陣特征的基礎上，可以魯棒地將矩陣分解為低秩部分和稀疏部分，又稱為魯棒主成分分析(Robust Principal Component Analysis，RPCA)[9]。傳統(tǒng)的RPCA模型是非凸的，且是一個NP(non-deterministic polynomial)難問題，一般不易直接求解。在RPCA的基礎上，一系列改進的算法被眾多學者提出，Zhou等人考慮到在干擾較多的環(huán)境中RPCA模型分解得到的結果不理想，提出去分解(Go Decomposition，GoDec)算法認為觀測矩陣是由前景、背景和噪聲組成，同時使用了雙邊隨機投影(bilateral random projections，BRP)算法加速矩陣前背景分離[10]。該算法雖然在一定程度上加快了模型求解的速度，但是存在難于選擇合適參數(shù)對秩、稀疏度范圍約束的缺陷。由此Zhou等人優(yōu)化了GoDec算法，提出半軟去分解(Semi-Soft GoDec，SSGD)算法[10]，稀疏度約束值可以通過軟閾值處理的方式自適應確定，但是該算法依舊需要對背景矩陣的秩進行估計。在背景矩陣的秩估計方面，Li Li等人通過計算待分解矩陣的信息熵線性擬合估計背景矩陣的秩，但信息熵與背景矩陣的秩不完全滿足線性關系會影響秩估計結果[11]。Jain P等人提出基于奇異值投影(Singular Value Projection，SVP)的秩估計方法在背景矩陣秩較低時會出現(xiàn)迭代不連貫[12]。

本文在李安冬等人研究低秩度參數(shù)和重構精度關系的基礎上[13]，將基于誤差重建的背景矩陣秩估計方法引入到半軟去分解算法中，提出了一種自適應低秩稀疏分解的超聲缺陷回波檢測方法，對金屬材料超聲缺陷檢測信號進行數(shù)據(jù)處理，可以有效減少無用信息并提高缺陷識別率。該方法首先對含噪信號進行短時傅里葉變換(short time Fourier transform，STFT)并提取幅度譜和相位譜；然后采用誤差重建的方法估計幅度譜背景矩陣的秩，基于估計的低秩度參數(shù)執(zhí)行SSGD算法將帶噪信號幅度譜分解為稀疏、低秩和噪聲3部分并舍棄噪聲部分；最后通過時頻掩蔽分離出缺陷信號幅度譜并結合相位譜逆STFT變換得到缺陷回波信號。

1 自適應低秩矩陣分解

1.1 SSGD算法

假設受到噪聲干擾的矩陣為Y，在其滿足低秩和稀疏的優(yōu)化準則下，GoDec算法可以有效的估計出Y中的低秩部分L和稀疏部分S，如下式所示

式中:N代表噪聲部分，rank(L)表示矩陣L的秩，card(S)表示矩陣S的基數(shù)，即非零元素的個數(shù)。對于式(1)可以通過最小化分解誤差求解:

對于式(2)可以交替求解以下兩個參數(shù)的子問題:

式中:‖·‖F(xiàn)表示Frobenius范數(shù)。針對GoDec算法的參數(shù)r選擇過大會導致一部分噪聲泄露到稀疏部分S中的缺陷，Zhou等人提出了半軟去分解算法[10]:

式中:λ是軟閾值參數(shù)，‖·‖1表示l1范數(shù)。不同于GoDec算法對低秩部分L和稀疏部分S都進行硬閾值處理，SSGD對式(1)中S項采用軟閾值處理，通過軟閾值λ自適應確定約束“card(S)<＝r”中的稀疏度r的值。對于式(4)可以通過式(3)表述的交替優(yōu)化方法求解。此外，為了避免對大規(guī)模矩陣進行奇異值分解，SSGD算法引入雙邊隨機投影(Bilateral Random Projections，BRP)來求解式(3)，從而進一步提高了計算效率。

1.2 低秩矩陣構建

低秩稀疏分解算法處理的對象是含有潛在低秩和稀疏結構特征的矩陣，考慮到超聲信號的時頻譜圖能夠突出信號在不同時間下頻率能量分布特征，而缺陷檢測中的非聲學噪聲信號多以自相似的背景混響為主，在時頻譜圖上具有相似或重復的結構特點，超聲檢測中缺陷反射信號則為不平穩(wěn)的信號，其在時頻域的稀疏性是一個廣泛接受的觀點[4]。因此，本文采用短時傅里葉變換構建低秩矩陣。

通過有限長的窗口移動加權對非平穩(wěn)的原始信號進行交疊分段得到短時平穩(wěn)的數(shù)據(jù)幀，然后對預處理后的每一幀信號進行快速傅里葉變換(fast Fourier transformation，F(xiàn)FT)，其計算公式為

式中:N為信號長度，W為幀長，h(n)(n＝0，…，W-1)為歸一化窗函數(shù)，R＝t i-t i-1為相鄰兩幀交疊的點數(shù)，計算得到的?Y即為描述原始信號的時頻譜。

1.3 自動秩估計

SSGD算法在求解時頻譜的低秩和稀疏成份前，需要預估低秩成份的秩作為迭代求解的先驗信息[13]，準確的秩估計值可以避免分解出的稀疏成份受到噪聲成份的干擾。但是受采樣率和雜波的影響，背景實際的秩往往難以準確估計。當預估背景矩陣的秩明顯小于矩陣實際的秩時，分解得到的低秩成份L會有大量干擾信號泄露到稀疏成份S和噪聲成份N中，則缺陷回波信號可能會被虛假回波干擾導致微小信號漏檢；當設定的秩值逐漸遞增但依舊小于實際的秩之前，低秩成份L中包含的信號量會逐漸增加，記重建誤差為‖N‖F(xiàn)＝‖Y-(S+L)‖，其變化量為ΔN，此時重建誤差的變化量較大；當背景矩陣的秩設定接近或者大于實際的秩時，分解得到的低秩成份L包含了大量具有重復結構特征的干擾信號，此時重建誤差較之前有明顯的下降后，同時其變化量也趨于穩(wěn)定，但是隨著背景矩陣的秩設定越來越大時，稀疏部分S中的缺陷回波信號又會混入低秩成份中導致漏檢。因此根據(jù)信號重建誤差的變化趨勢自適應確定SSGD算法中的參數(shù)k，可以保證幅度譜中的缺陷回波和干擾信號有效區(qū)分。逐步遞增設定的秩值并分解幅度譜，得到重構信號和重構誤差。當重建誤差平穩(wěn)時，有‖ΔN‖F(xiàn)≤η。選擇重建誤差平穩(wěn)后對應的第一個秩值作為k值輸入到SSGD算法。算法具體流程如圖1所示。

圖1 低秩度參數(shù)估計流程圖

1.4 增強信號波形重構

對時頻譜取模即可得到幅度譜Y，在使用SSGD算法估計出Y中分別代表稀疏特征和重復特征的幅度譜S和L后，需要采用二值時頻掩蔽法[14]處理以得到最終需要的缺陷回波幅度譜。二值時頻掩蔽法假設S和L的混合幅度譜中每個時頻點只可能包含一種信源，故可以通過保留只包含稀疏成份的時頻區(qū)域而去除包含低秩成份的區(qū)域，最大限度的消除稀疏成份中的干擾信號。本文的二值時頻掩蔽計算公式如下:

式中:l和n代表幅度譜中時頻點索引。將計算出來二值時頻掩蔽值作用于估計出的低秩和稀疏部分的幅度譜，即可得到缺陷檢測信號的幅度譜:

結合原始信號STFT時記錄下的相位信息∠·，對估計出幅度譜進行逆STFT就可以得到最終的回波信號。

1.5 基于自動秩估計的SSGD矩陣分解算法

基于自適應低秩矩陣分解的超聲缺陷回波檢測步驟如下:

步驟1 根據(jù)式(5)使用漢明窗(Hamming)對觀測信號進行STFT得到時頻譜，取模得到幅度譜，記為Y。

步驟2 在各低秩度參數(shù)下計算重建誤差‖N‖F(xiàn)，記錄滿足誤差變化量小于η的第一個低秩度參數(shù)k。

步驟3 根據(jù)式(3)和式(4)在秩為k的條件下對幅度譜Y進行分解，得到噪聲部分N、低秩部分L和稀疏部分S，記去除噪聲部分N后的信號幅度譜為Y L S

步驟4 根據(jù)式(6)計算掩蔽值M，再根據(jù)式(7)對幅度譜Y LS進行二值掩蔽后可得到缺陷信號的幅度譜。

步驟5 最后對做逆STFT得到降噪后的時域回波信號。

算法流程圖如圖2所示。

圖2 方法流程圖

2 實驗

2.1 仿真缺陷回波檢測

超聲波是一種非平穩(wěn)時變信號，缺陷檢測接收的信號包括缺陷回波、材料晶粒散射形成的結構噪聲和非聲學噪聲等。在脈沖反射法中，單個缺陷回波s(t)可以表示為。

式中:β為幅值；τ為到達時間；φ為初相位；a(t)為單位峰值的包絡函數(shù)；f c為中心頻率。

考慮實際的缺陷回波包絡可能不對稱(通常為上升沿陡峭而下降沿平緩)，引入調幅回波模型(Amplitude-modulated echo model，AMEM)[15]，取

式中:α為帶寬因子；m為雙曲正切函數(shù)階次；r為不對稱因子-1

入射超聲波經(jīng)過金屬材料晶粒散射形成的結構噪聲可以建模為[16]

式中:b為大于零的常數(shù)；γ為材料衰減因子；K為為超聲探測聲束范圍內的晶粒數(shù)；σk為第k個晶粒散射回波的強度系數(shù)；ωk為第k個晶粒散射回波的頻率漂移；τk為第k個晶粒散射回波延時；ωc為超聲缺陷回波信號的中心頻率。

非聲學噪聲n(t)通常可以簡化為方差σ2的高斯白噪聲。綜上，含N個缺陷回波的檢測信號x(t)可以表示為

為簡化表述缺陷回波信號s i(t)的各個參數(shù)，記λi＝(βi，αi，r i，mi，τi，fci，φi)為s i(t)的參數(shù)向量。

利用以上模型生成超聲無損檢測信號，參數(shù)設置如下:缺陷回波個數(shù)N＝2，λ1＝(1，35，0.5，10，0.8，5，0)，λ2＝(0.8，35，0，10，3.2，5，0)構造回波信號s(t)。在信號長度范圍內材料的晶粒數(shù)量K設為2 500，衰減系數(shù)γ為10-28，σk服從瑞利分布[4]，ωk和τk服從均勻分布，構造u(t)。再添加能量大小為-3 dBW的隨機噪聲n(t)得到信噪比為-10.5 dB的含噪信號x(t)。仿真信號的采樣頻率取作f s＝500 MHz，采樣長度為2 000，仿真結果如圖3所示。

圖3 仿真信號

從圖3(b)可見，仿真信號x(t)信噪比很低，缺陷回波被淹沒。圖4為提取幅度譜后經(jīng)自適應秩估計得到的不同低秩度參數(shù)下低秩稀疏分解算法重建矩陣產(chǎn)生的誤差變化量曲線，其中誤差變化量經(jīng)歸一化處理，η根據(jù)大量仿真實驗在本文中?。?.1。圖5展示了部分低秩度參數(shù)下自適應矩陣分解得到的幅度譜低秩部分L、稀疏部分S和噪聲部分N經(jīng)逆STFT的時域波形。從圖中可以看出，低秩度參數(shù)設定低于實際值時，有大量干擾信號泄露到稀疏成份中；低秩度參數(shù)越高，分解得到的低秩矩陣包含的干擾信號越多，但稀疏矩陣中的缺陷信號幅值逐漸減弱甚至出現(xiàn)了丟失第二個幅值較小的缺陷回波的情況。所以準確估計待分解矩陣背景的秩，可以保證檢測率的同時可以有效降低泄露到缺陷回波的噪聲。對于本文的仿真信號，當k設定由2增加為3時，相應的誤差變化量首次低于η值，所以選擇低秩度參數(shù)為3進行低秩稀疏分解。對得到的低秩和稀疏部分進行二值時頻掩蔽處理，最終提取的缺陷回波信號如圖6所示，對比原始回波信號可以看到噪聲基本被消除，回波顯著，畸變程度較小。降噪后的信號幅值雖有所降低，但對于強噪聲下的缺陷檢測和定位并無影響。

圖4 重建誤差變化曲線

圖5 不同秩參數(shù)約束條件下分解得到的L、S和N對應的時域波形

圖6 去噪信號(t)和回波信號s(t)對比圖

為了評估本文所提方法在不同強度的噪聲下對缺陷檢測信號處理性能，改變仿真信號中噪聲的強度，分別得到信噪比為-12.17 dB、-8.56 dB、-5.55 dB、-3.09 dB四種不同的含噪信號x(t)，同時引入小波閾值(Wavelet thresholding denoise，WTD)，經(jīng)驗模態(tài)分解(Empirical Mode Decomposition，EMD)，EMD離散小波變換(EMD discrete wavelet transform，EMDDWT)與本文方法進行比較，處理結果如圖7所示。由時域波形圖可以看出，隨著信噪比的提升，四種方法的性能均逐步提升，但是前三種方法處理后的信號均產(chǎn)生了較大程度的畸變，尤其是在信噪比較低時，波形中仍含有較多噪聲，起振位置難以識別。其中EMD-DWT在信噪比較高時可以達到和本文所提方法相同的性能，但是其波形前后仍含有少量噪聲。本文所提方法處理后的信號整體效果最好，波形畸變程度最小，在缺陷回波之外幾乎不含有干擾雜波，信號失真更小，兩個缺陷回波的起振位置和振幅以及其他細節(jié)信息更易于識別，當信噪比處在較低水平時也依舊可以保持穩(wěn)定的性能。

圖7 不同噪聲強度下四種方法的處理結果對比

為了定量描述上述方法的性能，引入相關系數(shù)(COR)和均方根誤差(RMSE)作為評價標準。COR的計算公式為:

式中:和xorg分別是處理后的信號和原始信號，〈·〉表示內積運算。COR用來評價重構缺陷信號與原始回波信號之間的相似程度。

RMSE的計算公式為:

式中:N為用于分析的信號采樣長度。RMSE用于衡量重構缺陷信號(t)與原始回波信號s(t)之間的誤差。其中COR越大，RSME越小，去噪效果越好。去噪后信號的各項指標如表1所示。

分析表1可以得到和時域波形對比類似的結論，即相對于原含噪信號x(t)，在經(jīng)低秩矩陣分解處理后，信噪比和互相關系數(shù)的提升量最大，均方差的下降量最大，說明低秩矩陣分解可以在很大程度上去除噪聲，同時還能保持回波信號失真較小。

表1 不同噪聲強度四種方法處理結果評價指標

2.2 實際缺陷回波檢測

應用本文所提方法處理實測超聲缺陷回波信號。采用基于脈沖反射法的水浸超聲縱波對缺陷進行檢測，檢測示意圖如圖8(a)所示，收發(fā)一體式超聲探頭產(chǎn)生具有一定頻率間隔和持續(xù)時間的脈沖進入工件，在聲阻抗變化處發(fā)生不同程度的反射，其中按照入射路徑返回的部分脈沖被探頭接收。實驗中檢測對象為如圖8(b)所示的6061鋁合金，在側面依次鉆有兩個直徑2 mm的孔作為缺陷，使用其中距離上平面d＝60mm的孔F1作為代表進行檢測，材料縱波傳播聲速為c＝6 300 m/s。實驗裝置如圖9所示，超聲換能器中心頻率5 MHz，示波器采樣頻率f s為36 MHz。

圖8 缺陷檢測實驗系統(tǒng)

圖9 缺陷檢測實驗系統(tǒng)

實驗采集的超聲信號如圖10(a)所示，可以觀察到回波信號受到了較多高頻非聲學噪聲干擾，分析其頻譜(圖10(b))發(fā)現(xiàn)有大量接近或低于探頭中心頻率的頻率成分，主要由結構噪聲和缺陷回波引起[17]。

圖10 實驗采集的超聲信號及其頻譜

圖10 (a)中信號幅值首次和末次近似達到1 V對應的采樣點分別是86和120，取均值作為始波位置，記為n＝103。通過傳播聲速c和缺陷實際位置d計算得到缺陷回波的大概時間值t f＝n/f s+2d/c＝2.86μs+19.05μs＝21.91μs，底波可以通過同樣的方式計算得到，為了簡化分析，本文僅截取始波和底波間800個采樣點的數(shù)據(jù)進行分析，如圖11所示。

圖11 用于分析的信號段

最終處理結果如圖12所示，缺陷回波在圖中的采樣點對應的時間值在22.08μs附近，與實際回波時間點偏差很小。同時可以觀察到，四種方法都能較大程度抑制干擾信號，但前三種方法處理后的信號在缺陷回波外信號起伏較大，不利于波形起振位置的識別，而本文方法得到的信號無論從波形完整性還是波形平滑性上來看均優(yōu)于其他算法，去噪后信號僅含有少量的干擾信號，回波幅值和波形完整性得到較好的保留，更有利于缺陷回波的精確檢測。

圖12 不同方法處理結果對比

3 結論

針對當前基于SSGD方法分解低秩矩陣時，難以精確預估分解算法所需秩約束參數(shù)的困難，分析了SSGD模型中重建矩陣精度和秩參數(shù)的關系，對自適應確定背景矩陣秩參數(shù)進行了研究。

以STFT和低秩稀疏分解為基礎，提出了對一維超聲缺陷檢測信號執(zhí)行矩陣分解去噪的方法，改善了金屬超聲缺陷檢測中存在回波被干擾信號湮沒的狀況。結果表明將低秩稀疏分解技術應用于超聲檢測信號能獲得較理想的處理效果。

研究了本文所提方法和WTD，EMD以及EMDDWT方法在不同信噪比的仿真檢測信號下的處理效果，應用均方差和相關系數(shù)作為評價指標，量化了本文所提方法的信號處理性能，最后通過實測信號進行了驗證。