基于尾波干涉法反演結(jié)構(gòu)微小損傷

趙元明，郭懷攀

(中國飛行試驗研究院，陜西 西安 710000)

飛行試驗中每次飛行前后機務(wù)人員必須目視檢查飛機的進氣道、 掛架、 艙蓋等關(guān)鍵部位是否有損傷，但肉眼檢查并不能發(fā)現(xiàn)某些極小的損傷和暗紋.尾波干涉法為檢測這類損傷提供了可能性.多重散射波是振動信號或聲音信號在散射介質(zhì)中傳播時，由于介質(zhì)的非均勻性產(chǎn)生的散射現(xiàn)象.用脈沖信號對被測介質(zhì)進行激勵時，檢測到的響應(yīng)信號是脈沖直達波并且后面拖著長長的“尾巴”，這個“尾巴”稱為尾波.與直達波相比，由于多重散射波在介質(zhì)中來回多次傳播，對介質(zhì)的微小變化進行了重復(fù)的采樣，可以把微小變化進行放大，因此對檢測介質(zhì)的微小變化具有更好的靈敏度，這就是尾波干涉原理.

Pacheco & Snieder在2005年研究得出多重散射波對于散射介質(zhì)的微小變化的敏感性并非在空間中均勻分布，并根據(jù)振源和監(jiān)測點的位置、 散射介質(zhì)物理特性、 尾波的觀察時間等變量推導(dǎo)出敏感核的解析式[1].本文對尾波干涉理論及其敏感核構(gòu)造進行簡要介紹，并基于MATLAB建立數(shù)值仿真平臺，成功地反演出介質(zhì)的微小變化的大小、 形狀及位置分布，驗證了算法及程序的有效性，為尾波干涉法在工程實踐中的推廣及應(yīng)用奠定了基礎(chǔ).

1 尾波干涉

Snieder在總結(jié)前人工作基礎(chǔ)上提出尾波干涉(Coda Wave Interferometry，CWI)原理[2]，尾波是振動信號或聲音信號在散射介質(zhì)中經(jīng)過多次散射的結(jié)果，比直達波對介質(zhì)有更多的采樣，所以對介質(zhì)的微小變化更為敏感，可以識別直達波所不能識別的微小變化.

(1)

式中：δt(t)為介質(zhì)速度擾動前后同一監(jiān)測點接收到波形的走時差，s1和s2分別為振(聲)源和監(jiān)測點的位置，K(s1,s2,x0,t)為敏感核，其表達式為

p(x0,s2,t-τ)dτ,

(2)

式中：p為兩點之間的波場強度，是時間和位置的函數(shù).在均勻的二維散射介質(zhì)中，強度p的解析解有兩種，一種是擴散方程的解

(3)

式中：l為平均自由程，r為振源到接收點的距離，c為波速，Θ為Heaviside函數(shù).另一種解析解是散射方程的解

(4)

式中：D為擴散系數(shù).把兩種解析解帶入敏感核進行對比，證實擴散方程的解更能準(zhǔn)確地描述敏感核的分布[3], 因此在本文中，采用擴散方程的解構(gòu)造敏感核.圖 1 所示為當(dāng)t=1 s和t=3 s時, 平均自由程均為500 m，臺站間距均為3 km時敏感核的分布.從圖 1 中可以看出敏感核呈馬鞍形分布，分別在振(聲)源和檢測點處達到峰值，隨著與這兩點距離的變大，敏感核相應(yīng)的降低.對比t=1 s和t=3 s時刻的敏感核分布可以看出，t越大，敏感核分布越廣，說明從振(聲)源發(fā)出的信號到達接收點之前經(jīng)過的空間范圍越廣.

圖 1 尾波觀察時間為1 s和3 s時的敏感核分布Fig.1 Distribution of sensitive nucleus when the observation time of tail wave are 1 s and 3 s

m=m0+CmGT(GCmGT+Cd)-1(d-Gm0),

(5)

式中：Cd為測得的數(shù)據(jù)的對角協(xié)方差矩陣；GT為G的轉(zhuǎn)置，G=Δvk(r′,r,s,t)；m0為初始模型，為零矩陣；Cm為平滑矩陣[4],

式中：λ為相關(guān)長度，λ0為網(wǎng)格間距，δm為先驗標(biāo)準(zhǔn)差.

2 數(shù)值仿真平臺

數(shù)值仿真平臺主要由以下幾部分組成： 激勵信號、 速度場、 信號傳播方式、 監(jiān)測點的分布等，以此平臺為基礎(chǔ)，輸入激勵信號獲得響應(yīng)信號，計算獲得速度擾動的空間分布.數(shù)值仿真平臺結(jié)構(gòu)圖如圖 2 所示.

圖 2 數(shù)值仿真平臺結(jié)構(gòu)圖Fig.2 Structure diagram of numerical simulation platform

2.1 激勵源和擾動區(qū)域的設(shè)定

激勵源為一個脈沖信號，如圖 3 所示，信號持續(xù)時長為6π/100，激勵信號函數(shù)為

(6)

如圖 4 所示，在速度場的正中間施加一個速度+0.5%、 2 km×2 km的矩形速度擾動區(qū)域，根據(jù)監(jiān)測點的布設(shè)位置，在振(聲)源處發(fā)出單位脈沖激勵，11個監(jiān)測點連續(xù)記錄波形.

圖 3 激勵源波形Fig.3 Excitation source waveform

圖 4 在10 km×10 km的監(jiān)測區(qū)域正中央施加一個矩形2 km×2 km的速度擾動區(qū)域Fig.4 Add a rectangle with 2 km×2 km velocity perturbation region to the central of 10 km×10 km monitoring area

2.2 速度場和監(jiān)測點的布設(shè)

實際的傳播介質(zhì)是非均勻性的，因此在速度場中，用速度的微小波動模擬介質(zhì)的非均勻性.為了取得理想的仿真結(jié)果，假設(shè)速度場是無限大的，速度場的邊界應(yīng)該設(shè)置為吸收邊界，吸收邊界采用Clayton Engquist majda吸收邊界條件進行仿真，反射率大約為2.5%，滿足本研究的要求[5].

首先，在10 km×10 km的二維速度場中，邊界為吸收邊界，運用傳播速度不同模擬介質(zhì)的非均勻性，速度場的速度模型為200×200的矩陣，速度矩陣的均值為5 km/s, 方差為200，最小值為4 958.7 km/s，最大值為5 083.5 km/s，如圖 5 所示.

監(jiān)測點的布設(shè)原則是確保均勻地覆蓋監(jiān)測區(qū)域，通過嘗試不同的布設(shè)方案，發(fā)現(xiàn)監(jiān)測點布設(shè)位置對反演結(jié)果有一定的影響，這不在本文討論范圍，因此在本文中監(jiān)測點的位置固定不變，以消除監(jiān)測點位置不同對結(jié)果的影響，監(jiān)測測點的分布如圖 6 所示，1為振(聲)源，2～12為監(jiān)測點.

圖 5 二維速度場模型Fig.5 Two-dimensional velocity field model

圖 6 傳感器布設(shè)分布Fig.6 Sensor layout distribution

2.3 信號傳播方式

運用波動方程進行波的傳播計算.

(7)

式中：f(x,z)為振源，v為波的傳播速度，p為位移.

把等式兩邊的偏導(dǎo)數(shù)用泰勒展開式展開，整理后波動方程轉(zhuǎn)化為差分格式，從泰勒展開式的數(shù)學(xué)理論得知，展開式的階數(shù)越高越逼近于真實值，在權(quán)衡運算效率和精確度后，選擇了時間二階、 空間八階的有限差分格式.

pk+1(i,j)=2pk(i,j)-pk-1(i,j)+

s(t)*δ(i-i0)*δ(i-j0),

(8)

式中：uk(i,j)為在第k次(對應(yīng)的時間)迭代時在(i,j)處的位移；v為振動的傳播速度，是位置(i,j)的函數(shù)； Δh為網(wǎng)格間距，監(jiān)測區(qū)域內(nèi)網(wǎng)格為均勻分布，因此Δh為固定常數(shù)； Δt為采樣時間間隔，為采樣頻率f的倒數(shù).

3 仿真結(jié)果分析

圖 7(a) 為在2號監(jiān)測點記錄到的波形，最大波峰處為直達波，波峰后拖著長長的尾巴即為散射造成的尾波.由于尾波干涉法觀察的是介質(zhì)的微小變化，因此同一個監(jiān)測點處擾動前后的波形基本一致，只考慮相位差，忽略幅度的變化[6].

從圖 7(a) 中可以看到擾動前后的波形幾乎完全重合，這是因為介質(zhì)發(fā)生的是微小變化，所以對波形的幅度及相位的影響非常小.當(dāng)把波形放大觀察時，如圖 7(b) 中1～1.4 s的對比波形，仍然發(fā)現(xiàn)擾動前后波形高度重合； 而3.1～3.5 s時的波形，發(fā)現(xiàn)有很明顯的走時差，這是因為尾波把介質(zhì)的微小變化進行了放大.

圖 7 擾動前后的波形及部分波形放大圖Fig.7 Waveform before and after disturbance recorded by monitoring point and partial waveform enlargement

3.1 走時差和敏感核的計算

計算走時差常用的方法有兩種： 延展法和互相關(guān)法[7].延展法主要針對實驗數(shù)據(jù)信噪比較低的情況，通過對波形進行拉伸(壓縮)后與擾動前采集到的波形進行互相關(guān)計算獲得走時差，有較好的穩(wěn)定性； 互相關(guān)法是比較常用的一種方法，相對于延展法計算量較小，對于本文具有較高噪比的仿真數(shù)據(jù)，采用互相關(guān)法[8].選擇尾波時間段為3～7.7 s的數(shù)據(jù)計算走時差，每段數(shù)據(jù)長度為0.3 s, 重疊0.1 s, 每個檢測點可以得到15段尾波數(shù)據(jù)，用互相關(guān)法計算得到的走時差(以振源和2號監(jiān)測點為例)如圖 8 所示[9，10].

圖 8 走時差Fig.8 Time lag

以1號點位置為激勵源，其它11個監(jiān)測點記錄響應(yīng)信號，在擾動前后均可記錄到11個振動信號，對監(jiān)測點擾動前后波形進行互相關(guān)計算可得到總計11×15=165個走時差.

敏感核表示從振源發(fā)出的信號在某一時刻經(jīng)過空間某一點到達接收點的概率分布.使用散射方程的解計算敏感核.

k(r′,r,s,t)=

(9)

把10 km×10 km的觀測區(qū)域劃分為20×20的網(wǎng)格，網(wǎng)格間隔為500 m.得到振源和2號監(jiān)測點，時間為3.25 s時刻的敏感核如圖 9 分布，從圖中可以看出，敏感核為馬鞍形，在震源和監(jiān)測點處敏感核數(shù)值較大.

圖 9 t=3.25 s時振源和2號監(jiān)測點的敏感核分布Fig.9 The distribution of sensitive nucleus of vibration source and No.2 monitoring point in the tail wave at t=3.25 s

3.2 速度擾動分布

圖 4 所示在監(jiān)測區(qū)域的正中央施加矩形 +0.5% 的速度擾動區(qū)域，圖 10 為反演結(jié)果，黑色方框標(biāo)注出施加的擾動區(qū)域，顏色條表示反演得到速度擾動得百分比，從圖中看出反演結(jié)果與施加的擾動區(qū)域吻合度較好，不僅成功地定位了速度擾動區(qū)域的位置，并且反演出擾動的幅度.

這種介質(zhì)內(nèi)的振動(聲音)傳播速度的微小變化可能由于各種原因?qū)е?，例如溫度?結(jié)構(gòu)裂紋、 密度變化等，需根據(jù)具體介質(zhì)進行分析.同時也進行了其它區(qū)域的反演計算，圖 11 是在監(jiān)測區(qū)域的左上方施加+0.5%速度擾動、 右下方施加-0.5% 速度擾動的反演結(jié)果.

圖 10 正中央施加+0.5%的速度擾動的反演結(jié)果Fig.10 Inversion results of velocity perturbation with +0.5% in the center

圖 11 施加+0.5%速度擾動和-0.5%速度擾動的反演結(jié)果Fig.11 Inversion results of +0.5% velocity perturbation and -0.5% velocity perturbation

4 結(jié)束語

本文對尾波干涉理論進行了簡要介紹，并基于MATLAB建立數(shù)值仿真平臺，證明了該算法能夠準(zhǔn)確監(jiān)測到結(jié)構(gòu)微小的損傷及定位損傷的具體位置，更深層次的研究以及工程應(yīng)用需進行繼續(xù)探究，尾波干涉法是一種較新的探傷技術(shù)，尤其是對介質(zhì)的微小變化具有較高的敏感性，其在飛機試飛、 橋梁工程、 石油開采等方面有廣泛的應(yīng)用.