基于貝葉斯公式的無參考信號Lamb波時間反轉(zhuǎn)損傷識別方法

1、基于貝葉斯公式的無參考信號Lamb波時間反轉(zhuǎn)損傷識別方法 繆傲（武漢大學(xué)土木建筑工程學(xué)院， 武漢 430072）摘 要 現(xiàn)有的主動Lamb波損傷監(jiān)測大多采用基于參考信號的差信號方法獲取損傷散射信號，然后用確定性的方法進行損傷定位，因此在適用性方面受到很大影響。本文引入Lamb波的時間反轉(zhuǎn)過程，提出了一種無參考主動Lamb波損傷識別方法，且考慮損傷識別過程中的不確定性因素，通過貝葉斯公式，推導(dǎo)出損傷位置、波速等未知參數(shù)的聯(lián)合后驗概率分布后，利用蒙特卡洛馬爾科夫鏈（MCMC）方法對未知參數(shù)進行采樣估計，所得到的馬爾可夫鏈的極限分布即為未知參數(shù)的后驗分布。最后，通過在鋁板上的數(shù)值仿真結(jié)果表明，該方法

2、能夠較為準(zhǔn)確的識別出損傷的位置、大致范圍等特征。 關(guān) 鍵 詞 損傷監(jiān)測；時間反轉(zhuǎn)；無參考；貝葉斯；MCMC1 前言V1V2V12V11V21V21V22Lamb波作為一種超聲導(dǎo)波，在結(jié)構(gòu)中衰減慢，能夠傳播較遠的距離，并且對結(jié)構(gòu)中的小損傷敏感，所以在結(jié)構(gòu)健康監(jiān)測領(lǐng)域中的應(yīng)用越來越廣泛。由于損傷會引起結(jié)構(gòu)中傳播的Lamb波的散射，故損傷散射信號直接與損傷位置有關(guān)?，F(xiàn)有的主動Lamb波損傷監(jiān)測方法大多基于參考信號，即以健康狀態(tài)響應(yīng)信號為基準(zhǔn)，通過差信號的方法獲取損傷散射信號，但真實結(jié)構(gòu)和外界環(huán)境的變化對其影響很大，因此在適用性方面受到了很大的影響，且傳統(tǒng)的基于損傷散射信號傳播時間的損傷定位方法中，例

3、如，脈沖回波法、橢圓定位法、四點圓弧定位法等，把損傷散射信號傳播時間、一定頻厚積下的波速都作為確定值，但實際的量測誤差和損傷定位中不確定性是不可避免的。例如，在損傷識別前，損傷的尺寸、程度是我們所不知道的，損傷對波速變化產(chǎn)生的影響也是不確定的，這將反過來影響通過理論計算所得到得損傷散射信號傳播時間，再者，用做激勵和傳感的壓電片的尺寸也沒有考慮，這都給損傷的識別帶來誤差。除此之外，由于Lamb波的頻散特性，損傷散射信號波包在傳播過程中會發(fā)生畸變，從而也會對損傷散射信號傳播時間產(chǎn)生影響。在這種情況下，概率的方法在損傷識別過程中能夠考慮各種不確定性因素，因此，比用確定性的方法更合適。針對上述的諸多問

4、題，本文首先根據(jù)板結(jié)構(gòu)中主動Lamb波時間反轉(zhuǎn)理論，分析了存在單損傷板結(jié)構(gòu)中Lamb波時間反轉(zhuǎn)聚焦過程，提出了Lamb波無基準(zhǔn)損傷散射信號傳播時間的提取。然后從概率的角度，考慮損傷識別過程中模型本身和量測所存在的不確定性和誤差，把利用時間反轉(zhuǎn)所得到的傳感器列陣中損傷散射信號傳播時間D作為損傷位置、波速等未知參數(shù)的樣本信息，再結(jié)合其先驗信息，采用貝葉斯公式，得出未知參數(shù)的聯(lián)合后驗概率分布。最后通過MCMC方法，得到單個未知參數(shù)的后驗分布。2基于貝葉斯方法的時間反轉(zhuǎn)損傷定位2.1單損傷板結(jié)構(gòu)中Lamb波的時反聚焦時間反轉(zhuǎn)法是法國科學(xué)家Fink最先由光學(xué)應(yīng)用引入到Lamb波領(lǐng)域中，并開展了大量的理論

5、和實驗研究工作。Lamb波的時反是指將感器所接收到的各個模式的信號在時間域上反轉(zhuǎn)后，在傳感器上加載，即所接收到的信號先到后發(fā)，后到先發(fā)，從而實現(xiàn)了各個模式信號在原始激勵處的自適應(yīng)聚焦。已有的對板結(jié)構(gòu)中傳播的Lamb波傳感信號時反特性研究表明：不管是在單模式還是雙模式（S0，A0模式）下，Lamb波的頻散都得到了補償，同一模式經(jīng)過前向和反向傳播后都在主波峰處下聚焦，而在向前、反向的傳播時，不同的群速度將會產(chǎn)生相應(yīng)的的時間延遲，在主波峰旁辨形成旁瓣。為了降低信號處理的復(fù)雜度，選擇合適中心頻率的窄帶信號，采用雙面激勵方式，以激發(fā)出單模式的Lamb波信號。如圖a所示，對于單模式的Lamb波信號，由A傳

6、播到B存在兩條路徑，PZTB所接收到的傳感信號為直達波信號和損傷散射信號。如圖b所示，將B處的傳感信號經(jīng)時域反轉(zhuǎn)后再次加載在PZTB上，在PZTA處得到的時反重構(gòu)信號可以表示為： 式（7）中，為B點接收到的直達波信號再次以直接波的方式到達A處，、的命名規(guī)則一樣。根據(jù)式（5）可得： 將式（8）帶入式（7）可得： 式（9）中基于壓電片的激勵模型，傳遞函數(shù)、可以簡化為由幅值和速度的方程表為： 式中：為激勵信號直線傳播至傳感器的路徑，為激勵信號經(jīng)損傷散射后再傳播至傳感器的路徑，為幅值頻散方程，為散射系數(shù)，為模式的波數(shù)?？紤]到此確定的模型中，傳播距離和窄帶激勵信號的中心頻率為確定值，故式（11）中的幅值

7、項也為確定值，為方便描述，。將式（11）帶入式（9），并進行傅立葉反變換，可得時域內(nèi)的重構(gòu)信號為： 為了獲得這個方程的近似解，通過泰勒級數(shù)在中心頻率處展開得： 帶入角頻率、波數(shù)、相速度和群速度的關(guān)系、，上式變?yōu)?其中， 帶入公式（14）、（15），公式（12）可表示為由式（17）可以看出：經(jīng)過時間反轉(zhuǎn)處理后，重構(gòu)信號主要由三部分組成，第一部分為聚焦信號的主瓣，與初始激勵信號的時反信號形狀相同，實現(xiàn)了初始激勵信號的重構(gòu)，Lamb波傳播過程中的頻散得到了補償；第二、三部分與主瓣部分有的時間延遲，對稱的出現(xiàn)在主波峰旁辨形成旁瓣。同時，重構(gòu)信號主瓣與旁瓣波峰對應(yīng)的時間差正好為損傷散射信號傳播時間于直達

8、波傳播時間的時間差。在主動監(jiān)測的傳感器列陣中，直達波傳感路徑已知，故根據(jù)時反響應(yīng)信號主、旁瓣之間的時間差就可以確定損傷散射信號的傳播時間，因而可以實現(xiàn)無參考信號的損傷定位和識別。貝葉斯分析方法貝葉斯分析方法是基于假設(shè)的先驗概率、以及把量測所得到的數(shù)據(jù)作為樣本信息而得出的。由于考慮了先驗信息和樣本信息，故未知參數(shù)的不確定性減少。隨著工程結(jié)構(gòu)的不確定性等問題中概率方法的發(fā)展，貝葉斯方法在結(jié)構(gòu)健康監(jiān)測中起著越來越重要的作用。在損傷的評價、損傷預(yù)測、傳感器優(yōu)化布置、監(jiān)測系統(tǒng)優(yōu)化設(shè)計等工程上都得到了廣泛的研究和應(yīng)用。與一般確定性的方法給出未知參數(shù)的精確解相比，貝葉斯方法的基本理念是把未知參數(shù)向量的分布看

9、成隨機變量聯(lián)合分布，然后對每個參數(shù)做出在區(qū)間范圍內(nèi)的估計。具體為：將時間反轉(zhuǎn)所得的Lamb波損傷散射信號傳播時間D作為樣本信息，綜合未知參數(shù)的先驗分布，再根據(jù)貝葉斯公式，得出后驗信息，然后根據(jù)后驗信息去推斷未知參數(shù)。對于薄板結(jié)構(gòu)，考慮具有條激勵傳感路徑的壓電片列陣，對于第（i=1,2）條傳感路徑的Lamb損傷散射信號傳播時間可以示為lamb波從驅(qū)動器傳播到損傷再到傳感器所經(jīng)歷的時間。如圖是一個簡單的模型用來計算第條傳感路徑的Lamb損傷散射信號傳播時間。在理論上，第i 條傳感路徑的Lamb損傷散射信號傳播時間可以表達為： 式（18）中和分別是第條分別激勵傳感路徑中損傷、驅(qū)動器、傳感器位置的中心

10、坐標(biāo)，是在一定的激勵頻率時lamb波損傷散射信號的傳播波速。在式（18）中，壓電片的尺寸、損傷的范圍都沒有考慮，材料屬性的不確定性也會使實際波速偏離理論值，故在公式1中沒有采用特定的頻厚積下所對應(yīng)的Lamb波波速，而是把波速也作為一個未知參數(shù)用來確定損傷的中心位置。因此在該損傷定位方法中，未知參數(shù)向量可以表示為。在下文中，為了方便，用來表示中的未知參數(shù)。在未知參數(shù)被確定后，損傷的定位則被轉(zhuǎn)化為從概率的角度來解決確定未知參數(shù)分布的問題。假定建模誤差和量測的不確定性分別用變量和表示，那么第條激勵傳感路徑的損傷散射信號的傳播時間從概率的角度可以表示為：是在未知參數(shù)向量下，利用公式（18）所計算得到的

11、第條傳感路徑所對應(yīng)的損傷散射信號的傳播時間。為了簡單和方便，通常將和假定為平均值為0，標(biāo)準(zhǔn)差分別為和的相互獨立的正態(tài)分布。由式（19）得： 變量服從正態(tài)分布，參數(shù)和的似然函數(shù)可以表示為： 式（21）中，方差，似然函數(shù)為在公式1和給定的未知參數(shù)向量下，對于所測得的損傷散射信號傳播時間 D的概率分布。 為了方便，將公式3中的求和部分記為 由貝葉斯公式，根據(jù)所測得的TOF值，結(jié)合未知參數(shù)的先驗概率密度函數(shù)和使然函數(shù)，未知參數(shù)的后驗概率密度函數(shù)可以表示為 是 的聯(lián)合后驗分布函數(shù)， 是的聯(lián)合先驗分布函數(shù)。是一個常數(shù)，在后驗分布中僅起到一個正則化因子的作用。對于每個參數(shù)的邊緣概率密度函數(shù) 表示參數(shù)向量中對

12、除外的其它參數(shù)和做一個多維積分根據(jù)所測得的損傷散射信號傳播時間，公式7給出了每個變量的邊緣分布函數(shù)之后，根據(jù)后驗分布函數(shù)，可以得到每個未知參數(shù)的后驗分布，那么損傷的位置和波速度也就能夠確定。但是在實際應(yīng)用中，由于涉及到多個參數(shù)的高維積分，不便于用分析的方法計算，使得貝葉斯方法的應(yīng)用大為受阻。運用MCMC方法，我們可以從復(fù)雜的多維積分分布中采樣，從而能夠?qū)ω惾~斯估計方法進行很好的模擬采樣由于在貝葉斯方法中為一個常數(shù)，在計算后驗分布中僅起到一個正則化因子的作用，如果將去掉，可將貝葉斯公式改寫為可將貝葉斯公式改寫為 其中符號表示左右兩邊相差一個常數(shù)因子式的右邊雖然不是正常的概率密度函數(shù)，但它確是后驗

13、分布的主要部分，故只要式的右邊已知，就可以從后驗分布中采樣，從而避免了計算的問題。MCMC方法蒙特卡洛模擬中，在后驗分布中取樣，當(dāng)這些樣本獨立時，利用大數(shù)定律樣本均值會收斂到期望值。但在貝葉斯分析領(lǐng)域中，后驗分布是復(fù)雜的、高維分布，很難得到獨立的樣本。目前MCMC方法已經(jīng)成為一種處理復(fù)雜統(tǒng)計問題的特別流行的工具，尤其在經(jīng)常需要復(fù)雜的高維積分運算的貝葉斯分析領(lǐng)域，通過合理的定義和實施，MCMC方法總能得到一條具有平穩(wěn)分布的馬爾科夫鏈，該馬爾可夫鏈的平穩(wěn)分布就是我們所需要的后驗分布。下面給出M-H算法1、構(gòu)造一個輔助的概率密度函數(shù)，通常被稱為“提議函數(shù)”；2、在第n次迭代時，根據(jù)從提議函數(shù)中產(chǎn)生一

14、個新狀態(tài)；3、計算接受概率 概率置我們稱之為接受提議，以概率置，我們稱之為拒絕提議。重復(fù)上述過程，當(dāng)n足夠時，知道得到的馬爾可夫鏈?zhǔn)諗康揭粋€平穩(wěn)分布。在利用馬爾可夫鏈進行抽樣時，在收斂之前的一段時間，各個狀態(tài)的邊際分布還不能認為是穩(wěn)定分布的，所以在進行估計之前，應(yīng)該把前面的這n-1次迭代值去掉，這個過程稱之為burn-in。提議函數(shù)可以是任意的分布形式，但是選擇合適的提議函數(shù)得到的馬爾可夫鏈能夠更好地收斂，通常，選取的提議函數(shù)滿足對稱型，即，這將使接受概率變得更簡單些。本文采用基于高斯分布和均勻分布的MCMC算法從后驗分布中對未知參數(shù)向量進行采樣。與傳統(tǒng)的MH算法相比，本文所提供算法結(jié)合吉布斯

15、采樣的理念：如果在其它參數(shù)值已知時，那么某個參數(shù)的條件分布可以作為后驗分布來使用。根據(jù)其它變量的當(dāng)前值，依次對分布的每個變量生成一個樣本序列。在此算法中，引入優(yōu)化參數(shù)L，調(diào)節(jié)調(diào)節(jié)接受概率，使得馬爾可夫鏈經(jīng)過適當(dāng)?shù)牡螖?shù)后收斂。結(jié)合未知參數(shù)的先驗分布、量測所得到得損傷散射信號傳播時間，運用MCMC算法，可以得到聯(lián)合概率分布中各個位置參數(shù)的后驗分布概率。參數(shù)可以直接從伽馬分布中進行采樣，對于參數(shù)的估計，具體算法如下：設(shè)定總的迭代次數(shù)NT，在收斂之前的一段時間，其分布不能認為是穩(wěn)定分布，所以在估計時，應(yīng)該把前面NB的迭代值去掉。1、初始迭當(dāng)n=0時，； 從中采樣 2、n=1時， 從U(0,1)中產(chǎn)

16、生一個隨機數(shù)u，當(dāng)u<r時， ，調(diào)節(jié)參數(shù)，；否則 ，調(diào)節(jié)參數(shù)， 式中： k=1時，； k=2時 ， k=3時 ， ， 直接從分布中進行采樣。重復(fù)步奏直到，當(dāng)n超過Nb時，調(diào)節(jié)參數(shù)不再變化。舍棄前Nb次的迭代值，剩下的采樣 值可以作為一條收斂的馬爾可夫鏈，可以代表位 置參數(shù)的后驗分布。 數(shù)值仿真結(jié)果本文設(shè)計了相應(yīng)的驗證系統(tǒng)，通過數(shù)值仿真方法對上述理論結(jié)果進行驗證。仿真試件及壓電原件布置情況如圖所示：采用鋁板結(jié)構(gòu), 楊氏模量72Gpa,密度2700kg/m3，泊松比0.3,其尺寸為580mm* 580mm*2mm，在其相同位置的上下表面布置壓電片，構(gòu)成雙面激勵方式。以鋁板中心作為坐標(biāo)原點，選

17、用4對壓電片作為激勵/傳感器，其編號和坐標(biāo)分別為：1#（-90，-90）、2#（90，-90）、3#（90,90）、4#（-90,90）。3數(shù)值仿真結(jié)果1 2 1 1 1 200mm 1 180mm 1 200mm 1 O仿真結(jié)果1 580mm 1 580mm 1 4數(shù)值仿真結(jié)果1 In the plate view In the thickness view (a)平面視圖 (b)側(cè)視圖 鋁板與PZT布置方案的幾何模型 The geometric model of Aluminum plate and PZT layout scheme 應(yīng)用ANSYS12.0進行建模，結(jié)構(gòu)中損傷是以（20,

18、10）為圓心，直徑為8mm的圓形通孔，采用刪除有限元單元的方法模擬。劃分網(wǎng)格時，單元形狀選六面體，基體單元類型采用三維8節(jié)點實體單元solid45, 單元尺寸為1mm* 1mm*1mm，損傷周圍局部加密，壓電原件選用三維耦合場體單元solid5。為了保證壓電層與金屬面板重合部分傳遞力和位移的連續(xù)性,建模中采用黏結(jié)處理，模型的邊界全都設(shè)置為自由。在仿真實驗中，激勵信號采用正弦調(diào)制窄帶信號，表達式為：式中為激勵信號的中心頻率，取為100kHz,為階梯函數(shù)，為調(diào)制信號的波峰數(shù)，取為5，激勵信號的時頻域特性如圖所示。激勵信號 (a)時域波形（b）頻域圖Load signal(a)Time domain

19、 input signal （b） Frequency spectrum 實驗過程中，采用雙面同相激勵，垂直于板平面方向加載，此時產(chǎn)生的Lamb波以A0模式為主。逐個選取壓電原件為激勵源，其余的壓電元件作為傳感器，根據(jù)聲波互易性原理，有6個通道的信號是獨立的，表示為、 、(下標(biāo)數(shù)字表示激勵器和傳感器的編號)，然后對采集到的各監(jiān)測通道結(jié)構(gòu)響應(yīng)信號施加合適的窗函數(shù)進行截取，消除邊界反射信號，僅保留直達波信號和損傷散射信號，并對截取后的信號進行時間反轉(zhuǎn)并用作為二次激勵，得到結(jié)構(gòu)的時反聚焦信號。最后提取時各監(jiān)測通道時反重構(gòu)信號中的損傷損傷散射傳播時間用于損傷識別。取通道相應(yīng)信號進行分析Response signal of sensor波相應(yīng)信號波形截取并作時間反轉(zhuǎn)后的信號The signal after the interception and time reversal時反聚焦響應(yīng)信號Reconstructed Time reversedsingal 根據(jù)式所表述的時反聚焦信號中主、旁瓣信號對應(yīng)波峰的時間關(guān)系，提取各監(jiān)測通道下的損傷散射信號傳播時間，如表所示。對于監(jiān)測通道 、，由于損傷距離直達波路徑很近，信號散射角度大、傳播距離遠，加上Lamb波的頻散效應(yīng)導(dǎo)致的波包延拓和幅值衰減，損傷散射