基于廣義逆算法的沖擊波超壓場重建方法*

郭亞麗，韓 焱，王黎明

(中北大學(xué)電子測試技術(shù)國家重點實驗室，山西 太原 030051)

爆炸技術(shù)越來越多地應(yīng)用于國防及民用的各個領(lǐng)域[1]，沖擊波超壓是爆炸產(chǎn)生殺傷和破壞作用的主要因素，一直是空氣中爆炸規(guī)律研究的重要課題[2]。目前，雖然已總結(jié)了一些計算沖擊波超壓的經(jīng)驗公式[3]，但在實驗研究時，由于測試環(huán)境、測試條件以及裝藥形狀不完全相同，直接利用經(jīng)驗公式進(jìn)行計算是不準(zhǔn)確的。針對目前局部布點測試不能全面了解沖擊波傳播過程的不足以及超壓經(jīng)驗公式的局限性，本文中，擬采用網(wǎng)絡(luò)化測試技術(shù)獲取沖擊波信號，以計算機(jī)層析成像技術(shù)為基礎(chǔ)，利用加權(quán)廣義逆反演算法對爆炸沖擊波速度場進(jìn)行反演，并根據(jù)沖擊波速度與峰值超壓的關(guān)系，得到峰值超壓場的分布。

1 沖擊波超壓場重建原理

炸藥在空氣中爆炸，形成空氣沖擊波。將測試區(qū)域劃分為規(guī)則的網(wǎng)格單元，假設(shè)自由空氣中的理想沖擊波不沿網(wǎng)格單元的邊界傳播，在每個網(wǎng)格單元內(nèi)近似沿某條直線傳播，每個探測器 對應(yīng)一條射線。沖擊波到達(dá)探測器 陣元的時間為走時，沖擊波在傳輸過程的走時是速度和幾何路徑的函數(shù)。根據(jù)走時層析成像原理，有[4]

DS=T

(1)

式中：T=(t1,t2,…,tm)T為各條射線走時的m維列向量,S=(s1,s2,…,sn)T為待求離散單元慢度值，為n維未知的列向量；D為m×n階稀疏矩陣，其元素為dij，即第i條射線穿過第j個網(wǎng)格的射線長度。

2 觀測系統(tǒng)優(yōu)化布局

2.1 測試區(qū)域網(wǎng)格劃分

將測試區(qū)域劃分為若干規(guī)則的網(wǎng)格單元，每個網(wǎng)格單元內(nèi)波速視為均勻。網(wǎng)格劃分越密，層析成像的分辨率越高，但方程的未知數(shù)越多，解的不確定性也越高。網(wǎng)格劃分應(yīng)結(jié)合測試區(qū)域大小、先驗信息、實驗所用探測器 數(shù)目合理劃分。

2.2 探測器優(yōu)化布局及判斷指標(biāo)

探測器布設(shè)時，應(yīng)使射線覆蓋面廣、分布均勻；減少射線路徑矩陣中零元素的個數(shù)、降低其條件數(shù)。采用射線密度、射線正交性及矩陣D的條件數(shù)作為探測器優(yōu)化布局的判斷指標(biāo)。

射線密度代表通過各網(wǎng)格像元的射線數(shù)目，而射線正交性是通過各網(wǎng)格像元的射線之間夾角的最大正弦值來度量，射線密度小和正交性差的區(qū)域，圖像誤差大；反之，則結(jié)果比較可靠。在射線總數(shù)一定時，根據(jù)探測目標(biāo)形狀，模型分布特點等合理分布射線[5]。

3 基于加權(quán)廣義逆的沖擊波場反演算法

從矩陣的觀點求方程組(1), 也就是求系數(shù)矩陣D的逆矩陣的問題。但是對于爆炸層析成像問題，系數(shù)矩陣D一般是奇異矩陣, 不存在通常意義下的逆矩陣, 因此需要采用廣義逆理論進(jìn)行求解。對于沖擊波層析成像，除采用穩(wěn)定性好的廣義逆反演算法，還應(yīng)對反演矩陣進(jìn)行加權(quán)處理形成加權(quán)廣義逆反演，來增加反演過程中的信息量。

在廣義逆反演時，既可對觀測數(shù)據(jù)加權(quán)，也可以對模型參數(shù)加權(quán)。設(shè)A為m×n階矩陣，P、Q分別為m×m、n×n的正定矩陣[6]，若有n×m階矩陣X滿足

AXA=A，XAX=X， (PAX)T=PAX， (QXA)T=QXA

(2)

則稱X為A的加權(quán)廣義逆A+，A+可以表示為A+=Q-1(PAQ-1)+P。

采用自然權(quán)矩陣，設(shè)P、Q均為對角陣，對稱正定，矩陣P的對角線元素為當(dāng)前模型的各射線的走時；矩陣Q的對角線元素為各射線在第j個網(wǎng)格上的總貢獻(xiàn)，表示為射線穿越任何網(wǎng)格單元的長度和與相應(yīng)單元速度乘積，J.G.Berryman[7]稱其為射線分布矩陣。

這種權(quán)矩陣的最大特點是物理意義明確,表征了以下基于物理原理的判斷：(1)沖擊波作為信息載體,隨傳播距離的增加,衰減幅度增大，波形信號減弱，故射線路徑較短者,其信噪比應(yīng)較高,也較接近真實路徑,故對小走時測值加重權(quán)。(2)一般射線越密集,重建圖像就越精確?；谶@樣的考慮,對眾多射線穿過的網(wǎng)格單元加重權(quán)。于是，加權(quán)廣義逆走時層析成像算法寫為S=Q-1(PDQ-1)+PT。

圖1 爆炸場速度分布模型Fig.1 The model for velocity distribution in explosion field

圖2 網(wǎng)格及探測器分布示意圖 Fig.2 Schematic distribution of grids and detectors

4 數(shù)值模擬及實驗結(jié)果

4.1 數(shù)值模擬

設(shè)置10 m×10 m自由場空氣爆炸沖擊波速度分布模型,見圖1。根據(jù)實驗條件，將測試區(qū)域劃分為10×10的網(wǎng)格單元，炸點位于測試區(qū)域中心，見圖2。根據(jù)對稱性原理，對1/4區(qū)域進(jìn)行速度場反演。選擇的探測器數(shù)為13，選取2種不同的布設(shè)方式,見圖3。圖3(a)中探測器節(jié)

圖3 探測器的2種不同的布設(shè)方式Fig.3 Two different distributions of detectors

圖4 不同的探測器布設(shè)方式下的射線密度分布Fig.4 Half-line density distributions with two different distributions of detectors

圖5 不同的探測器布設(shè)方式下的射線正交性分布Fig.5 Half-line orthogonality distribution with two different distributions of detectors

圖6 速度分布反演的初始模型Fig.6 The initial model for velocity inversion

圖7 各個網(wǎng)格單元速度相對誤差Fig.7 The relative error of velocity in each grid

點基本呈均勻分布，其坐標(biāo)分別為:S1(0,4.9)、S2(0,4.15)、S3(0,3.32)、S4(0,2.5)、S5(0,1.66)、S6(0,0.83)、S7(0,0)、S8(0.83,0)、S9(1.66,0)、S10(2.5,0)、S11(3.32,0)、S12(4.15,0)、S13(4.5,0)。圖3(b)為按照探測器布設(shè)指標(biāo)優(yōu)化布局，各節(jié)點坐標(biāo)為:S1(0,4.9)、S2(0,3.8)、S3(0,3)、S4(0,2)、S5(0,1.3)、S6(0,1)、S7(0,0.5)、S8(1.2,0)、S9(2,0)、S10(2.5,0)、S11(3,0)、S12(3.5,0)、S13(4.9,0)。

2種探測器布局方式下,測試區(qū)域射線密度及正交性分布見圖4～5。由圖4～5可以看出:射線總數(shù)一定時，探測器優(yōu)化布局方式下射線分布更合理，正交性更好；同時，矩陣D的條件數(shù)為35.61，遠(yuǎn)小于均勻布局方式下的2.74×1016?？梢?，探測器優(yōu)化布局使得方程更穩(wěn)定。采用相同初始模型(圖6)和加權(quán)廣義逆算法，對2種布局下的爆炸場速度進(jìn)行反演重建，相對于真實模型(圖1)，探測器均勻布局和優(yōu)化布局的網(wǎng)格平均相對誤差分別為8.47%和3.32%，各個網(wǎng)格ng單元速度相對誤差εv見圖7，反演結(jié)果的二維圖像見圖8。由反演誤差及圖像分布均可得到：采用探測器優(yōu)化布局得到的反演結(jié)果誤差更小，更接近真實模型。

圖8 不同的探測器布局方式下的反演結(jié)果Fig.8 Inversion results with two different distributions of detectors

圖9 實驗中探測器布局Fig.9 Detector distribution in experiment

圖10 沖擊波峰值超壓場的重建結(jié)果 Fig.10 The reconstruction results for peak overpressure of shock wave

4.2 實驗及重建結(jié)果

實驗時,使用50 kg散裝TNT炸藥，炸藥密度為0.80 g/cm3，形狀接近柱形，吊離地面1.7m，進(jìn)行空中爆炸。對32 m×32 m范圍內(nèi)自由沖擊波場進(jìn)行二維重建，爆炸點居于測試區(qū)域中心，使用了30個超壓探測器 ，且探測器 與炸點在同一平面，可近似認(rèn)為是在無限空氣中的爆炸。根據(jù)有效的射線總數(shù)和探測區(qū)域的大小以及射線密度將測試區(qū)域劃分為長、寬各64×64個大小相同的網(wǎng)格，探測器 陣元分布在網(wǎng)格邊界和網(wǎng)格內(nèi)個別節(jié)點上，如圖9所示。

采用本文算法對爆炸沖擊波峰值超壓場進(jìn)行重建，同時，采用如下經(jīng)驗公式計算峰值超壓：

(6)

將某次實驗測得的沖擊波峰值超壓pm,e與重建結(jié)果pm,r及經(jīng)驗公式計算結(jié)果進(jìn)行比較，見表1。

通過多次實驗數(shù)據(jù)與重建結(jié)果的比較，得到在近場區(qū)4 m以內(nèi)，重建結(jié)果與實測值的偏差較大，接近10%，這是由于實驗時，近場區(qū)干擾較大，實驗數(shù)據(jù)本身也存在一定誤差，而且測試點數(shù)較少，導(dǎo)致重建結(jié)果與實測值的偏差增大，4 m以外的測試區(qū)域內(nèi)，采用本文算法重建得到的峰值超壓與實測值較接近，相對偏差在5%以內(nèi)，優(yōu)于經(jīng)驗公式計算結(jié)果。

實驗時，若增加超壓探測器 數(shù)目，重建精度會大大提高。峰值超壓場二維重建結(jié)果如圖10所示，從圖中可以看出，測試區(qū)域內(nèi)沖擊波峰值超壓隨距離的分布情況，距離爆炸中心5 m以內(nèi)的區(qū)域超壓衰減及其分布形狀較明顯，且其分布形式與裝藥形狀基本一致；而5 m以外的區(qū)域由于沖擊波超壓衰減幅度較大，故分布不明顯。

5 結(jié) 論

針對沖擊波超壓經(jīng)驗公式的局限性，采用網(wǎng)絡(luò)化測試技術(shù)獲取沖擊波信號，利用不完全數(shù)據(jù)重建技術(shù)得到一定區(qū)域內(nèi)爆炸沖擊波速度場及峰值超壓場分布。在一定范圍內(nèi)，重建結(jié)果與實測值較接近。然而在爆炸近場區(qū)測試時，對實驗儀器性能要求較高，所測數(shù)據(jù)本身存在一定誤差，且所測點數(shù)較少，故重建結(jié)果與實測值的偏差較大。本文中給出了一種沖擊波速度和峰值超壓的計算方法，但還需要更多的實驗數(shù)據(jù)進(jìn)行驗證。

[1] 張守中.爆炸基本原理[M].北京:國防工業(yè)出版社,1988,397-530.

[2] W E貝克.空中爆炸[M〗.江科,譯.北京:原子能出版社,1982:40-85.

[3] 李翼祺,馬素貞.爆炸力學(xué)[M].北京:科學(xué)出版社,1992,259-262.

[4] Kepler W F, Bond L J, Frangopol D M.Improved assessment of mass concrete dams using acoustic travel time tomography: II: Aapplication[J].Construction and Building Materials, 2000,14(3):147-156.

[5] 裴正林,余欽范,狄?guī)妥?井間地震層析成像分辨率研究[J].物探與化探,2006,26(3):218-224.Pei Zheng-lin, Yu Qin-fan, Di Bang-rang.Resolution in crosshole seismic tomography[J].Geophysical and Geochemical Exploration, 2006,26(3):218-224.

[6] 王振宇,劉國華,梁國錢.基于廣義逆的層析成像反演方法研究[J].浙江大學(xué)學(xué)報:工學(xué)版,2005,39(1):1-5.Wang Zhen-yu, Liu Guo-hua, Liang Guo-qian.Study on inversion methods for travel time tomography based on generalized inverse theory[J].Journal of Zhejiang University: Engineering Science, 2005,39(1):1-5.

[7] Berryman J G.Fermat’s principle and nonlinear travel time tomography[J].Physical Review Letter, 1989,62(25):2953-2956.