基于MUSIC及波束成形算法的非線性Lamb波微裂紋成像與評估

張 凡，徐 旭，李義豐

(南京工業大學 計算機科學與技術學院，南京 211899)

板狀結構在生產建設中被廣泛應用于機械、化工、建筑、航空航天等領域[1-3]。在長期承受載荷的應力作用下，易產生疲勞損傷，繼而發展成為有損于結構強度的宏觀損傷，引發結構破壞，生產事故。因此，結構體的早期微裂紋檢測是安全生產的重要保障手段[4-5]。

Lamb波是一種在薄板中傳播的導波，具有傳播距離長和對內部缺陷敏感的特點，被廣泛地應用于板結構的無損檢測領域中。傳統的Lamb波檢測技術多基于線性理論，根據缺陷引起回波的聲速、衰減、相位、反射系數等特性的變化來實現損傷的檢測，但對于早期的亞波長微裂紋不敏感[6-10]。結構中的微裂紋具有接觸非線性效應，會在透射和反射波中產生高次諧波，這為微裂紋 定位和評估提供了新的可能性。

傳統的基于Lamb波的損傷定位成像算法多采用飛行時間(time of flight,TOF)或基于損傷指數的概率損傷成像算法。Sohn等[11]通過使用16個PZT貼片傳感器陣列并基于時間反轉方法實現了對復合板中的缺陷檢測。Huang等使用概率損傷算法高效地實現了板中損傷定位成像。Caterina等[12]對“T”型碳纖維增強層合(carbonfiber reinforced polymer，CFRP)結構的常規損傷進行了定位成像研究。此類算法都是基于稀疏陣列構建的，具有構建簡單和計算量小的優點。TOF類算法需要人為地區分激勵信號、損傷的回波信號以及邊界反射信號波包，并通過測量損傷回波信號的飛行時間，計算損傷和激勵源的位置，實現損傷定位成像。但是，微裂紋的回波幅值較小，易淹沒在系統噪聲和其他回波信號之中，難以準確地分別損傷回波信號，故而該類方法無法實現微小損傷的精確定位成像[13-16]。而概率損傷成像方法僅能對布置傳感器時規劃的檢測路徑上的損傷實現精確成像，而對規劃路徑外的損傷則不具備高精度的成像定位能力[17]。

近年來，基于非線性Lamb波的檢測研究越來越引起研究者的關注。Jiao等通過對微裂紋混頻非線性技術的研究，使用非線性Lamb波信號對板中微裂紋情況進行了檢測；祝偉光等[18]利用非線性參數，對微裂紋生長方向的識別進行了研究；Hong等和Wang等[19]則分別對非線性Lamb波使用概率損傷算法和橢圓定位算法進行了微裂紋定位研究，但其算法定位準確性較差，定位精度不高。Lamb波在傳播經過微裂紋時，會與損傷相互作用產生新的散射波。因此，通過布置傳感器陣列，依托波達方向估計(direction of arrival,DOA)算法來實現損傷的定位成像。波達方向估計是一類使用陣列信號確定外界發射信源和陣列之間相對的方向信息的陣列信號處理方法，被廣泛的應用于現代雷達、通信、汽車、聲納中。此類算法中最廣泛的應用是利用噪聲子空間正交性的多重信號分類(multiple signal classification,MUSIC)算法[20-24]，具有高分辨力、強抗噪聲的特點。使用MUSIC算法進行陣列信號處理可以有效的解決非線性Lamb波的低幅值和易淹沒在系統的非線性噪聲中的問題，實現對微裂紋的高精度定位成像。

另外，微裂紋的損傷評估是結構健康檢測的重要內容。鑒于微裂紋非線性效應產生的二次諧波容易被各種噪聲的干擾所淹沒，本文在使用多重信號分類算法獲得損傷的準確定位后，進一步采用波束成形算法，實現對信號的空間濾波，進而提高損傷方向的信號強度，并根據接收波的非線性參數，完成對板中微裂紋情況的定量評估[25-26]。

1 理論基礎

1.1 非線性Lamb波基本理論

Lamb波是一類在無牽引、均勻和各向同性板中傳播的平面應變彈性波。其傳播受Navier方程約束。

(1)

式中：λ0和μ0為材料的拉梅常數；ρ為密度；bi為每單位質量的力；u為位移；i，j為坐標索引。

針對材料的非線性，其應力應變關系可以表示為

σ=Eε(1+βε+…)

(2)

式中：E為彈性模量；β為二階彈性系數；σ和ε分別為應力、應變。忽略二次以上的高階項，代入一維波動方程，可得

(3)

式中：u為位移；x為波傳播的距離；t為時間；c為波速。展開得

(4)

對于較小的振動幅值時，設其解得形式為

u(x,t)=u(0)+βu(1)

(5)

式中，u(0)，u(1)分別為通解和二階攝動解。一般情況下，攝動解與傳播距離成正比，故有

u(1)=xf(τ)

(6)

式中：f(τ)為任意一二階可微且連續的函數；τ=t-x/c。

本文使用激勵信號為正弦信號，故可表示為

u(0)=A1cos(f1t-k1x)

(7)

將式(6)、式(7)代入式(5)求解得

(8)

因此，式(5)可寫作

(9)

設二次諧波分量為A2，則有

(10)

不同于材料本身的非線性，微裂紋引起的聲學非線性而產生高次諧波的物理機制是接觸聲學非線性。Lamb波在傳播經過微裂紋時，會牽引其界面而產生半波整流效應。這種半波整流效應源于接觸部位剛度的不對稱特性。在壓縮過程中，微裂紋受Lamb波的作用而閉合，拉伸過程中，則會使微裂紋開放，使得在壓縮過程中的剛度高于拉伸過程的。這意味著Lamb波的壓縮過程可以傳播過界面，而拉伸過程會受到界面的阻擋，而無法通過界面。該效應導致了在透射波和散射波中高次諧波的產生，會在激勵頻率的整數倍頻率處產生分量。通過比較Lamb波信號的基波和二次諧波振幅，可以定性和定量的實現對微裂紋的評估。

圖1 微裂紋的接觸聲非線性Fig.1 Contact acoustic nonlinearity of micro cracks

圖1(a)和圖1(b)中：虛線表示微裂紋的初始位置；實線表示在Lamb波的壓縮和拉伸作用下，微裂紋左側界面的實際位置。圖1(c)為Lamb波通過微裂紋后的半波整流效應示意圖。盡管實際界面間的相互作用過程非常復雜，但其反射波非線性參數的形式皆可定義為二次諧波分量與基頻分量平方的比值。故定義相對非線性參數β′，用以衡量微裂紋情況。

(11)

1.2 MUSIC算法基本理論

MUSIC算法從陣列接收信號的協方差矩陣分離出信號子空間和噪聲子空間,而后利用信號方向向量與噪聲子空間的正交性構建空間掃描譜,實現信號的波達方向估計。常用于信號接收的線狀陣列建模如圖2所示。

設線狀陣列由N個間距為d的子陣元組成，遠場有M個不相干的信源，其與陣列首陣元的法線方向所成夾角θ為入射角。陣列接收到的信號矩陣X(t)可以通過式(12)描述。

X(t)=A(θ)S(t)+N(t)

(12)

式中，A(θ)，S(t)及N(t)的表達式分別為

(13)

S(t)=[S1(t),S2(t),…,SM(t)]T

(14)

N(t)=[n1(t),n2(t),…,nN(t)]

(15)

式中：A(Θ)為線狀陣列的方向矩陣；S(t)為M個從遠場入射的信源所發射信號矩陣，并分別與陣列首個陣元的法線方向成夾角Θ=[θ1,θ2,…,θM]T；N(t)為噪聲矩陣。

圖2 線狀陣列示意圖Fig.2 Schematic diagram of linear array

接收信號的互相關矩陣為Rx=E[XXH]，其中H表示矩陣的共軛轉置。假定噪聲為互不相關的零均值加性的高斯白噪聲。將式(12)代入后，可得

Rx=E[(A(θ)S(t)+N(t))(A(θ)S(t)+N(t))H]=
AE[SSH]AH+E[NNH]=ARSAH+RN

(16)

式中：RS=E[SSH]為信號相關矩陣；RN=σ2I為噪聲相關矩陣；σ2為噪聲信號的方差；I為M×M的單位矩陣。現假定協方差矩陣Rx的特征值為{λ1,λ2,…,λN}，所對應的特征向量依次為{q1,q2,…,qN}。在多信源不相干且信源個數Mλ2>…>λM>…λN，其中,M個較大的特征值對應于信號空間，余下的N-M個較小的特征值對應于噪聲空間。信號特征值和噪聲特征值對應的特征向量分別為信號特征向量和噪聲特征向量。

根據特征值和特征向量定義

Rxqi=λiqi

(17)

對于M

Rxqk=λkqk=σ2qk

(18)

此時，將式(16)左右同乘qk可得

Rxqk=(ARSAH+σ2I)qk=ARSAHqk+σ2qk

(19)

(20)

AHqk=0

(21)

式(21)表明，噪聲特征向量和陣列方向選擇向量相正交。利用此正交特性，通過求解陣列信號的互相關矩陣對應的噪聲特征向量EN，對信源到達方位角進行搜索，生成空間譜，實現波達方向估計。空間譜函數定義為

(22)

本文使用以上的MUSIC算法對非線性Lamb波回波信號進行處理，以實現對微裂紋的高精度定位成像。

在實際檢測中，聲學非線性效應很微弱，產生的二次諧波的幅值會遠遠小于基波的幅值，極易被淹沒在噪聲中。因此，為實現對微裂紋的定性和定量分析，需要開發針對微弱的二次諧波的信號的增強算法和抗干擾檢測技術。本文通過提取非線性Lamb波的高次諧波信號及對非線性參數的分析，完成對微裂紋的定位成像和損傷程度評估。

1.3 波束成形算法

線狀陣列接收到信號后，通過將多個陣元接收到的數據進行延時、加權、求和等一系列處理，使得處理后的信號對指定的方向具有較高的增益而對其他方向的信號的增益較低，實現空間濾波。傳統波束形成通常按以相移波束形成或時延波束形成的方法，選擇合適的加權信號，使得指定方向的信號被加強而其他方向接收到的信號被減弱，這有利于從陣列信號中完成對幅值較低，易被淹沒在噪聲中的微裂紋信號的提取。

本文采用基于時延加權求和方法的波束成形算法(conventional beamforming,CBF)實現對微裂紋信號的空間濾波，通過式(23)實現對陣列信號的波束成形。

(23)

式中：n為陣列的序號；N為陣列的陣元數；θ0為所形成的波束的指向角；f和c分別為Lamb波的頻率和波速。

通過該波束成形算法的使用，可降低陣列信號中其他方向雜波所造成的干擾，確保非線性參數評估方法的準確性。

2 有限元仿真模型

本文使用有限元軟件ABAQUS 建立了三維鋁板的仿真模型，通過顯式動力學分析法對非線性Lamb波在鋁板中的傳播進行分析。鋁板的物理參數如表1所示，其幾何參數為600 mm ×600 mm ×4 mm，如圖3所示。通過嵌入Cohesive界面單元創建埋藏微裂紋模型，其上下界面設置為硬接觸，且為法向無位移的面-面相互作用。以板的幾何中心點作為原點，長寬方向分別為坐標系的x軸和y軸方向，厚度方向為z軸方向。以坐標(-50,50) 的A1點，和坐標(-50,-50)的B1點為首陣元，5 mm為陣元間距，分別設立由7個陣元構成的A、B兩個傳感器陣列。于坐標(-50,0)處施加對稱載荷，在板中激發出單一S0模態的Lamb波。激勵載荷信號為中心頻率150 kHz的經5周期漢寧窗調制的正弦波信號，時域波形和頻域頻譜如圖4所示。在本文中基頻信號的群速度約為5 440 m/s，波長為36.3 mm,與微裂紋相互作用產生的300 kHz非線性Lamb波二次諧波信號的波速約為5 210 m/s，波長為17.4 mm。

表1 鋁板材料參數Tab.1 Material parameters of aluminum plate

圖3 三維鋁板及埋藏裂紋示意圖Fig.3 The schematic of 3D aluminum plate model and buried crack damage modes

圖4 激勵信號Fig.4 Excitation signal

3 仿真結果與分析

3.1 微裂紋定位成像

于板中(87.5,7.5)處設立埋藏微裂紋(見圖3)。零厚度Cohesive界面單元的長度l為5 mm，高度h為3 mm。首先，對無損模型和含微裂紋模型的有限元仿真結果進行對比分析，驗證Lamb波的非線性效應。圖5(a)為微裂紋模型中陣列A首陣元所接收信號的時域信號，圖5(b)為首陣元所接收回波信號的頻譜。可見，微裂紋板相較于無損板，二倍頻處產生了由微裂紋而引起的非線性信號。

圖5 傳感器陣列接收的信號Fig.5 Signal received by sensor array

接收到的回波信號為同時包含基波和二次諧波分量混疊波形，為了實現損傷定位，需要濾除不能表征微裂紋的基波信號。本文設計并使用巴特沃爾茲帶通濾波器提取二次諧波信號，回波信號的波包通過濾波器前后對比時域波形圖及頻譜圖，如圖6所示。圖6表明了帶通濾波器在有效濾除了基波信號的同時，較好地保留了二次諧波分量。

濾波后的微裂紋回波信號經過MUSIC算法處理后，可獲得微裂紋位置與陣列的法線方向夾角。為實現損傷定位成像，需通過兩個不同位置的傳感器陣列，分別獲得微裂紋相對于兩個陣列的不同方向角，而后利用交叉定位的方法來實現對微裂紋的定位。由以A1點和B1點為首陣元的兩個線狀陣列A、B，分別計算得到這兩個陣列所接收信號的空間譜，結果如圖7所示，兩虛線分別表示微裂紋相對于陣列A、B的實際入射方向角。

圖6 濾波前后首陣元接收到的信號Fig.6 Original and filtered signals received by the first element of array

圖7 不同陣列回波信號處理后的空間譜Fig.7 Spatial spectrums calculated from different arrays’ echo signal

利用不同位置陣列所提取到的空間譜，分別對位于(72.5,2.5) 、(22.5,-2.5) 和(87.5,7.5)處的微裂紋進行MUSIC算法定位成像，結果分別如圖8(a)～圖8(c)所示，其中‘+’號標注為實際損傷的位置。圖8表明通過使用MUSIC算法處理非線性Lamb波信號可以獲得優秀的微裂紋定位成像。

在實際檢測中，微裂紋的回波信號幅值較低，易被淹沒在系統噪聲和其他回波干擾信號中，故難以使用傳統檢測方法實現定位。對無噪聲的仿真信號中添加高斯白噪聲以模擬真實環境中接收到的信號。信噪比為25 dB時的信號如圖9所示，可見回波信號被淹沒在噪聲中。在此情況中，傳統檢測方法不具備準確區分回波信號和雜波信號的能力。

圖8 基于MUSIC算法的非線性Lamb波微裂紋成像圖Fig.8 Micro crack imaging by using nonlinear Lamb wave based on MUSIC algorithm

圖9 25 dB信噪比時的回波信號Fig.9 Echo signal at SNR of 25 dB

本文的MUSIC算法具備較強的抗噪聲能力，在上述信噪比為25 dB的情況下，使用MUSIC算法對位于(87.5,7.5)處的微裂紋進行定位成像，結果如圖10所示。該圖表明在低信噪比情況下，基于MUSIC算法的非線性Lamb波對微裂紋仍具備較好的微裂紋定位成像能力。

圖10 信噪比為25 dB時微裂紋定位成像圖Fig.10 Imaging of micro crack location at SNR of 25 dB

3.2 微裂紋尺寸評估

利用非線性MUSIC算法，可以從回波信號中獲得微裂紋的位置信息。為了減少噪聲和其他方向反射Lamb波的干擾，本文使用波束成形算法對信號進行空間濾波，以提取獲得指定方向角上的回波信號。通過分析特定角度的回波信號，在微裂紋成像的基礎上，進一步實現對微裂紋情況的評估。

對坐標位于(72.5,2.5)處的微裂紋進行分析。首先，對上文中的微裂模型進行有限元分析，并通過波束成形算法對陣列信號進行空間濾波。圖11為對陣列所接收原始信號中的基頻信號及二倍頻信號分別進行濾波及波束成形處理后所得到的波形圖。

圖11 波束成形處理后的波形圖Fig.11 Waveform obtained after beamforming processing

本文分兩種情況，即分別單獨改變Cohesive界面單元的高度(h)、厚度(δ)兩組模型進行仿真分析。第一組仿真分析為保持Cohesive界面單元的長度l與厚度δ不變(h=5.0 mm,δ=300 nm)，選擇不同高度(h=1.5 mm,h=2.0 mm,h=2.5 mm,h=3.0 mm,h=3.5 mm)的裂紋模型；第二組仿真保持微裂紋長度l與高度h不變(l=5.0 mm,h=3.0 mm)，選擇不同厚度(δ=100 nm,δ=150 nm,δ=200 nm,δ=300 nm,δ=400 nm)的裂紋模型進行分析。

圖12為進行第一組仿真分析得到的結果，圖12(a)對不同高度的微裂紋模型所接收的信號波束成形處理后得到的頻譜圖，圖12(b)為對不同高度的微裂紋采用式(11)計算所得的非線性參數。圖12(b)表明隨著微裂紋的高度增加，回波信號的非線性參數會隨之而單調遞減。根據接觸非線性理論，當微裂紋的高度增加時，只有在微裂紋的端點附近處的界面會受到Lamb波的作用，產生持續的拍擊。微裂紋高度上的變化并不會引起拍擊面積的大幅度變化，在頻域上表現是二次諧波的幅值隨微裂紋高度不同的變化較小。而基波主要是由裂紋的反射產生的，其幅值隨著裂紋高度的增加而增加，且變化明顯。因此，表現為隨著微裂紋高度的增加，回波信號的非線性參數單調遞減。

圖12 對不同高度的微裂紋回波信號波束成形處理Fig.12 Beamforming processing of echo signal with different height h of microcrack

圖13為進行第二組仿真分析得到的結果，圖13(a)為對不同厚度的微裂紋模型所接收信號波束成形處理后得到的頻譜圖，圖13(b)為對不同厚度的微裂紋采用式(11)計算而得的非線性參數。該圖表明，隨著微裂紋的厚度增加，回波信號的非線性參數會隨之而單調遞減。根據接觸非線性理論，當微裂紋的厚度增加時，基波幅值在增加，而能產生接觸非線性的區域將減少，導致了二次諧波幅值的降低，因此，隨著微裂紋厚度的增加，回波信號的非線性參數單調遞減。

圖13 對不同厚度微裂紋回波波束成形處理Fig.13 Beamforming processing of echo signal with different deepths d of microcrack

4 結 論

本文研究了非線性Lamb波檢測中微裂紋定位成像和評估問題。利用有限元模擬分析方法，建立了包含埋藏微裂紋的薄鋁板結構模型。利用多輸入信號算法，實現了對埋藏微裂紋的定位成像，并通過向接收信號中人工添加高斯白噪聲來模擬實際檢測信號，驗證了該MUSIC算法具有非常好的抗干擾損傷定位成像能力。繼而結合波束成形算法，對陣列信號進行空間濾波處理，可有效提高微裂紋所在方向回波信號的幅值；通過頻譜分析可知該回波信號的非線性參數與微裂紋高度呈遞減關系，與微裂紋厚度呈下降關系。說明回波信號的非線性參數可以有效反映板中損傷的情況，通過分析該非線性參數可實現對微裂紋損傷情況的評估。