基于尾波干涉法反演結構微小損傷

趙元明，郭懷攀

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

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

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

1 尾波干涉

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

(1)

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

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

(2)

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

(3)

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

(4)

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

圖 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為測得的數據的對角協方差矩陣；GT為G的轉置，G=Δvk(r′,r,s,t)；m0為初始模型，為零矩陣；Cm為平滑矩陣[4],

式中：λ為相關長度，λ0為網格間距，δm為先驗標準差.

2 數值仿真平臺

數值仿真平臺主要由以下幾部分組成： 激勵信號、 速度場、 信號傳播方式、 監測點的分布等，以此平臺為基礎，輸入激勵信號獲得響應信號，計算獲得速度擾動的空間分布.數值仿真平臺結構圖如圖 2 所示.

圖 2 數值仿真平臺結構圖Fig.2 Structure diagram of numerical simulation platform

2.1 激勵源和擾動區域的設定

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

(6)

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

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

圖 4 在10 km×10 km的監測區域正中央施加一個矩形2 km×2 km的速度擾動區域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 速度場和監測點的布設

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

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

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

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

圖 6 傳感器布設分布Fig.6 Sensor layout distribution

2.3 信號傳播方式

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

(7)

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

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

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

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

(8)

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

3 仿真結果分析

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

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

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

3.1 走時差和敏感核的計算

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

圖 8 走時差Fig.8 Time lag

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

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

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

(9)

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

圖 9 t=3.25 s時振源和2號監測點的敏感核分布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 所示在監測區域的正中央施加矩形 +0.5% 的速度擾動區域，圖 10 為反演結果，黑色方框標注出施加的擾動區域，顏色條表示反演得到速度擾動得百分比，從圖中看出反演結果與施加的擾動區域吻合度較好，不僅成功地定位了速度擾動區域的位置，并且反演出擾動的幅度.

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

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

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

4 結束語

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