淺層地下爆炸層析成像信息優(yōu)化獲取技術(shù)*

白苗苗,王黎明,高 蕊,范學(xué)軍,王 可

(中北大學(xué)信息與通信工程學(xué)院,太原 030051)

?

淺層地下爆炸層析成像信息優(yōu)化獲取技術(shù)*

白苗苗,王黎明*,高 蕊,范學(xué)軍,王 可

(中北大學(xué)信息與通信工程學(xué)院,太原 030051)

基于地震層析成像理論提出了地下爆炸沖擊波超壓時空場重建技術(shù)。從模型離散化和傳感器布局兩方面進(jìn)行成像信息優(yōu)化,提高地下爆炸場層析成像信息的有效性。模型離散化方面,根據(jù)沖擊波的的衰減規(guī)律,提出了不規(guī)則網(wǎng)格劃分方法;傳感器布局方面,根據(jù)數(shù)學(xué)上欠定問題理論與地震波傳輸規(guī)律,提出了指導(dǎo)傳感器布局的指標(biāo)及優(yōu)化算法,實現(xiàn)傳感器合理布局。仿真表明,層析成像信息優(yōu)化獲取技術(shù)可以最大限度地提高數(shù)據(jù)利用率,使得反演過程更穩(wěn)定,反演誤差更小。

傳感器;優(yōu)化布局;地震層析成像;超壓場重建;地下爆炸

武器在地下爆炸,他們的破壞作用主要是爆炸產(chǎn)生的沖擊波和應(yīng)力波來摧毀目標(biāo)。準(zhǔn)確地評價這些高效能武器系統(tǒng)的殺傷力是武器系統(tǒng)在研制和靶場驗收等過程中急需解決的問題,因此需要準(zhǔn)確可靠的測試手段為武器系統(tǒng)的研究提供依據(jù)。

目前爆炸測試的主要方法是在不同爆心距布放傳感器獲取爆炸信息,或者進(jìn)行數(shù)值模擬[1-2],這種局部的方法只能測得空間有限位置的爆炸參數(shù)值,對于不規(guī)則彈藥造成的不均勻毀傷,局部點測試方法不能全面獲取爆炸信息,從而不能準(zhǔn)確評估爆炸威力。本文基于地震層析成像的理論,提出地下爆炸場層析成像,即在有效測試區(qū)域內(nèi),以有限測試點代替無限測試點,實現(xiàn)地下爆炸的超壓場時空分布,為武器系統(tǒng)毀傷效能評估提供依據(jù)[3]。

地下爆炸場層析成像歸結(jié)為單一激勵源稀疏射線層析成像問題,其射線呈稀疏分布,層析成像方程是欠定的、病態(tài)的。此類層析成像問題,信息的有效利用是非常關(guān)鍵的。本文主要從模型離散化和傳感器優(yōu)化布局[4]兩方面研究地下爆炸層析成像信息優(yōu)化技術(shù),最大限度地提高信息的有效性,降低數(shù)據(jù)的冗余性,從數(shù)學(xué)根源上改善層析成像問題。

1 層析成像信息優(yōu)化獲取依據(jù)及評價方法

1.1 理論依據(jù)

本文采用層析成像的方法重建地下爆炸超壓場,即在中心炸點周圍布放M個傳感器,爆炸后沖擊波從中心炸點到傳感器形成M條傳播路徑射線;根據(jù)沖擊波的傳輸規(guī)律,將所要重建的區(qū)域劃分為N個不規(guī)則的離散網(wǎng)格單元,通過對各個傳感器陣元采集獲得的沖擊波信號進(jìn)行分析及特征提取,得到?jīng)_擊波到達(dá)時間(走時),采用走時層析[5-6]方法重建測試區(qū)域內(nèi)各網(wǎng)格單元沖擊波速度。最后根據(jù)速度與沖擊波峰值超壓的關(guān)系,將每個網(wǎng)格的速度值轉(zhuǎn)換成超壓值,重建出整個二維面上的平面沖擊波超壓值。

地下爆炸場層析成像歸結(jié)為求解矩陣方程:

DS=T

(1)

式中:T=(t1,t2…tM)′為各條射線走時,是一個M維列向量,為試驗測試值;S=(s1,s2…sN)′為待求離散單元慢度值,是一個N維的未知的列向量;D為M×N階稀疏矩陣。

爆炸場層析成像本就是一個病態(tài)方程組的求解問題,由式(1)可知,影響方程解的唯一因素是系數(shù)矩陣D,因此要改善病態(tài)不適定反演問題的穩(wěn)定性就要改善矩陣D的結(jié)構(gòu)。矩陣D的元素dij代表第i條射線在第j個網(wǎng)格單元內(nèi)的長度,故矩陣D的結(jié)構(gòu)取決于反演區(qū)域的模型離散和測試射線也就是傳感器的布局。因此,我們需要從模型離散化和測試傳感器布局兩方面進(jìn)行優(yōu)化設(shè)計。

圖1 地下爆炸場層析成像信息優(yōu)化獲取技術(shù)結(jié)構(gòu)圖

1.2 層析成像信息優(yōu)化獲取方案及評價方法

1.2.1 評價指標(biāo)

從矩陣的觀點求解方程(1),也就是求系數(shù)矩陣D的逆矩陣問題,但是對于地下爆炸成像問題,系數(shù)矩陣D一般是奇異矩陣,不存在普通意義下的逆矩陣,因此需要采用廣義逆理論進(jìn)行求解。通過對矩陣D進(jìn)行廣義逆分解發(fā)現(xiàn),矩陣D的奇異值和秩在求逆中擔(dān)任重要因素,因此可以將矩陣D的奇異值和秩作為評價反演問題穩(wěn)定性的主要指標(biāo),故定義評價指標(biāo)為:

(2)

式中:λ1為DTD的最大特征值;λi為DTD的所有特征值,P為矩陣D的秩。E由兩部分組成,第1部分反映了D中奇異值的相對分布情況,第2部分反映了D中零空間的相對大小情況。E越大,矩陣D的數(shù)學(xué)性質(zhì)越好,反演結(jié)果也越穩(wěn)定可靠。

條件數(shù)是評價方程組解質(zhì)量的又一個重要指標(biāo),條件數(shù)越大,方程組的病態(tài)程度越大。故定義評價準(zhǔn)則,

E2=cond(D)

(3)

射線密度小和正交性差的區(qū)域,反演誤差大;反之,則結(jié)果比較可靠。因此可以將射線平均密度和平均正交性作為布局準(zhǔn)則:

(4)

式中:ρj和Oj分別表示第j個網(wǎng)格內(nèi)射線密度和射線正交性,k1和k2的取值根據(jù)具體模型的數(shù)值ρj和Oj大小決定。在進(jìn)行層析成像信息優(yōu)化獲取時,盡可能增大射線密度和射線正交性,即E3越大越好[7-8]。

1.2.2 評價方法

由以上分析可知,由矩陣特征值和秩組成的評價函數(shù)E1、矩陣條件數(shù)E2。及射線覆蓋程度均可作為判斷矩陣D及方程組穩(wěn)定性的指標(biāo),E1數(shù)值越大,反演結(jié)果越穩(wěn)定,E2。越小,方程組越接近良性,同時,在進(jìn)行射線布設(shè)時,盡可能增大E3。本文將E1、E2、E3進(jìn)行加權(quán),形成一個綜合評價函數(shù),即:

E=ω1E1+ω2/E2+ω3E3

(5)

式中:ω1、ω2、ω3的取值,根據(jù)E1、E2和E3的數(shù)量級范圍進(jìn)行選擇,E越大越好。

由于影響病態(tài)欠定方程組的解的因素比較多,從觀測系統(tǒng)優(yōu)化布局方面而言,可以將評價函數(shù)E作為評判方程組穩(wěn)定性的指標(biāo)之一,在試驗之前指導(dǎo)觀測系統(tǒng)優(yōu)化布局。本文采用評價函數(shù)E對模型離散化方式和射線分布方式進(jìn)行評價,從而優(yōu)化試驗設(shè)計。在具體實施時,首先根據(jù)模型特點劃分網(wǎng)格,其次,根據(jù)評價函數(shù)E對發(fā)射器和接收器進(jìn)行優(yōu)化布局。具體步驟如下:

①根據(jù)模型特點劃分網(wǎng)格,利用評價函數(shù)對網(wǎng)格進(jìn)行評價,得到最佳網(wǎng)格模型;

②根據(jù)試驗所用發(fā)射器和接收器數(shù)目,隨機(jī)給出一種布局方式,計算矩陣D,并逐一判斷矩陣D的列向量是否為零向量;同時,計算矩陣D的秩;

③如果矩陣D的所有列向量都不為零向量,并且矩陣D是滿秩的,則將該組布局方式記錄作為優(yōu)化布局的初始模型;否則,返回②;

④當(dāng)初始模型達(dá)到設(shè)定數(shù)量時,將E作為目標(biāo)函數(shù),采用優(yōu)化算法,進(jìn)行優(yōu)化布局。

2 模型離散化和傳感器優(yōu)化布局

2.1 模型離散化原理

模型離散化即是將測試區(qū)域離散化為若干網(wǎng)格單元,每個網(wǎng)格單元內(nèi)有均勻的速度分布,網(wǎng)格劃分方式影響矩陣D的結(jié)構(gòu)。網(wǎng)格劃分越密,層析成像的分辨率就越高,但方程的未知數(shù)就越多,解得不確定性也就越高。網(wǎng)格劃分應(yīng)結(jié)合爆炸沖擊波傳輸規(guī)律、測試區(qū)域大小、先驗信息、試驗所用傳感器數(shù)目合理劃分。

地下爆炸時,爆炸波在距離爆炸點不同區(qū)段上出現(xiàn)沖擊波、應(yīng)力波和彈性波,如圖2所示。其中沖擊波是一種強(qiáng)壓縮波,波前、波后介質(zhì)的狀態(tài)參數(shù)急劇變化。實質(zhì)上沖擊波是一種特殊的壓縮波,一般壓縮波壓力是連續(xù)升高,而沖擊波壓力是急劇躍升。當(dāng)發(fā)生劇烈爆炸后,沖擊波的陣面也和爆轟波陣面一樣,表示截止?fàn)顟B(tài)特征的壓力、質(zhì)點運(yùn)動速度、密度和其他參數(shù)發(fā)生突躍。波陣面之前介質(zhì)的參數(shù)和未激發(fā)之前一樣,而波陣面之后,它們發(fā)生不連續(xù)變化。沖擊波以超聲速度傳播,沖擊波的波速、壓力和能量隨著距離很快地衰減[9]。

圖2 各階段上地震爆炸波形態(tài)

圖3 沖擊波峰值超壓、速度隨距離的變化

沖擊波峰值超壓隨距離衰減波形如圖3所示,雖然各經(jīng)驗公式計算結(jié)果有差異,但可以看出,在近場區(qū),超壓衰減較快,隨比例距離的增大峰值超壓趨于平緩。若采用均勻矩形網(wǎng)格劃分方法,往往導(dǎo)致震源處超壓和速度變化較快的區(qū)域出現(xiàn)較大的誤差,故對沖擊波峰值超壓、速度進(jìn)行重建時,不宜采用單一網(wǎng)格尺度,而應(yīng)根據(jù)沖擊波超壓、速度衰減規(guī)律,結(jié)合自適應(yīng)網(wǎng)格方法建立多尺度網(wǎng)格。

沖擊波峰值超壓與距離的關(guān)系可以用下式表示:

(6)

在實際應(yīng)用中由于實際資料的限制得到的射線條數(shù)是有限的,要想獲得相對高的成像質(zhì)量網(wǎng)格不能劃分太細(xì),否則有些網(wǎng)格內(nèi)的射線覆蓋次數(shù)會過低甚至為零;同時,若網(wǎng)格劃分太粗則不能保證正演時射線追蹤的計算精度,同時,也使得近場區(qū)反演誤差增大。

本文在分析沖擊波超壓衰減特點的基礎(chǔ)上,借鑒自適應(yīng)網(wǎng)格劃分原理,提出了不規(guī)則網(wǎng)格劃分方法,即根據(jù)爆心距離劃分區(qū)域并采用不同尺度的網(wǎng)格,在超壓衰減快的近場區(qū),采用細(xì)網(wǎng)格劃分方法,網(wǎng)格尺寸較小且網(wǎng)格分布稠密;在遠(yuǎn)場區(qū)采用粗網(wǎng)格劃分方法,網(wǎng)格尺寸大且分布稀疏。為避免矩陣D的某些行向量線性相關(guān),在對稱區(qū)域也采用不同的網(wǎng)格尺度,同時,網(wǎng)格劃分應(yīng)根據(jù)試驗所用傳感器數(shù)目,兼顧層析成像的分辨率和精度要求。

2.2 仿真實驗

以地下爆炸沖擊波速度為模型進(jìn)行仿真。測試區(qū)域為32 m×32 m的正方形區(qū)域,假設(shè)所使用的藥包形狀規(guī)則均勻且爆炸具有對稱性,只對正方形爆炸區(qū)域的90°方向區(qū)域內(nèi)爆炸性能進(jìn)行測試,爆炸點居于測試區(qū)域原點處,傳感器分布在測試區(qū)域邊界。

采用兩種離散模型方案對測試區(qū)域進(jìn)行網(wǎng)格劃分。第1種采用均勻矩形網(wǎng)格劃分方法,網(wǎng)格數(shù)是49個,如圖4所示。第2種采用本文所提出的不規(guī)則網(wǎng)格劃分方法,如圖5所示,根據(jù)沖擊波速度隨距離大致呈指數(shù)衰減規(guī)律,在距離爆心較近的區(qū)域,速度衰減幅度大,網(wǎng)格尺寸較小且網(wǎng)格分布稠密,故劃分為近場(0～4 m)、次近場(4 m～6 m)、中場(6 m～10 m)、遠(yuǎn)場(10 m～16 m),同時避免矩陣D的某些行向量線性相關(guān),將測試區(qū)域劃分為7個區(qū)域,網(wǎng)格數(shù)是58,各個區(qū)域網(wǎng)格大小不同。

仿真結(jié)果分析如下:

①為了使得各網(wǎng)格均有射線穿過,均勻矩形網(wǎng)格所需傳感器數(shù)目最少為9,而不規(guī)則網(wǎng)格所需傳感器數(shù)目最少為5個,如圖4、圖5所示。

②相同射線分布情況下,均勻矩形網(wǎng)格對射線分布要求較高,更容易出現(xiàn)沒有射線穿過的網(wǎng)格且多數(shù)情況矩陣D不滿秩,如圖6所示,13個傳感器同樣的分布,均勻矩形網(wǎng)格不滿秩,而不規(guī)則網(wǎng)格滿秩。

③相同射線分布情況下,不規(guī)則網(wǎng)格各項指標(biāo)均優(yōu)于均勻矩形網(wǎng)格。隨機(jī)選取13個傳感器在5種不同布站方式下,比較均勻矩形網(wǎng)格和不規(guī)則網(wǎng)格所得到的矩陣的秩、評價函數(shù)E1、條件數(shù)、射線正交性和密度,仿真數(shù)據(jù)如表1所示。由表可以看出,不規(guī)則網(wǎng)格更容易滿秩且其條件數(shù)遠(yuǎn)遠(yuǎn)小于均勻矩形網(wǎng)格,同時,不規(guī)則網(wǎng)格射線正交性和射線密度也均大于均勻矩形網(wǎng)格。

圖4 均勻離散化模型

圖5 不規(guī)則離散化模型

圖6 同樣數(shù)目的傳感器在兩種模型下的分布

不同模型實驗次數(shù)RankE1E2OρE11353．0801111．080．22043．33369．7516121356．04223285．90．22713．35299．51272不規(guī)則模型31357．34051178．80．18533．23539．3326541359．3204834．00．20593．352910．3248451359．7333972．40．21793．451010．6142311258．53690．8860×10160．13902．85718．8497921160．91321．2400×10160．14112．85719．08952均勻模型31160．38192．1470×10160．14072．89809．0768941256．87291．0202×10160．14462．69398．5257951257．57730．9448×10160．15092．75518．66373

2.3 傳感器優(yōu)化布局仿真

2.3.1 優(yōu)化布局算法

對于層析成像觀測系統(tǒng),一條射線對應(yīng)著一對發(fā)射和接收器的位置,優(yōu)化射線分布,即是要優(yōu)化發(fā)射器和接收器的布局,在本文中也就是優(yōu)化爆炸點和傳感器的位置。在離散模型一定的情況下,發(fā)射器和接收器的幾何位置是影響矩陣D結(jié)構(gòu)的主要因素。在給定發(fā)射器和接收器數(shù)量的情況下,采用優(yōu)化算法,將評價函數(shù)作為目標(biāo)函數(shù),尋求發(fā)射器和接收器的最優(yōu)幾何分布。本文利用自適應(yīng)粒子群[10-12]優(yōu)化算法進(jìn)行傳感器優(yōu)化布局,算法流程如下:

①在測試范圍內(nèi),根據(jù)發(fā)射器和接收器數(shù)目,隨機(jī)給出一種布局方式,即產(chǎn)生一個d維粒子,表示為a=(xs1,…,xsm1,ys1,…,ysm1,xr1,…,xrm2,yr1,…,yrm2)xsi和ysi分別表示第i個發(fā)射器的坐標(biāo),在本文中指中心炸點坐標(biāo);xri、yri分別表示第i個接收器,也就是本文的傳感器坐標(biāo);d=2m1+2m2,m1和m2分別表示試驗所用發(fā)射器和接收器數(shù)目;

②計算矩陣D,并逐一判斷矩陣D的列向量是否為零向量;同時,計算矩陣D的秩;

③如果矩陣D的所有列向量都不為零向量,并且矩陣D是滿秩的,則將該粒子加入粒子群;否則,放棄該粒子,返回①;

④粒子群數(shù)量達(dá)到所設(shè)定數(shù)量時,根據(jù)目標(biāo)函數(shù),計算群體和個體最優(yōu)適應(yīng)度值;

⑤更新粒子位置和速度,并把速度限制在c和vmin之間,計算適應(yīng)度值;

⑥更新個體極值,對于每個粒子,將其適應(yīng)度值與所經(jīng)歷過的最好位置的適應(yīng)值進(jìn)行比較,如果當(dāng)前的適應(yīng)度值比其歷史最優(yōu)值要好,那么歷史最優(yōu)將會被當(dāng)前位置所替代;

⑦更新群體極值,對于每個粒子,將其歷史最優(yōu)適應(yīng)度值與群體最優(yōu)適應(yīng)度值進(jìn)行比較,如果該粒子的歷史最優(yōu)值比群體最優(yōu)值要好,那么群體最優(yōu)值將會被該粒子的歷史最優(yōu)值所替代;

⑧檢查終止條件,終止迭代或返回⑤,重新計算。

2.3.2 仿真結(jié)果

建立如圖5所示的不規(guī)則離散模型,隨機(jī)選取13個傳感器,根據(jù)E1、E2和E3數(shù)值大小,w1=0.8,w2=0.2,w3=15,k1=0.1,k2=0.9。采用粒子群優(yōu)化算法,其各參數(shù)為:最大迭代次數(shù)q=100;初始化群體數(shù)目M=40;學(xué)習(xí)因子c1=2,c2=2;初始/終止慣性權(quán)重wmax=0.9,wmin=0.4;粒子最大更新速度為vmax=xmax-xmin,其中xmax和xmin為所求未知數(shù)上下界;粒子速度范圍[-vmax/2,vmax/2],得到最佳評價函數(shù)E從而得到最佳射線分布。結(jié)果如表2所示。

表2 各指標(biāo)在不同射線分布下的比較

比較3種不同布局的傳感器各指標(biāo)特點,如表2所示,可以看出,采用本文提出的算法優(yōu)化布局后,各項指標(biāo)均優(yōu)于其他分布,即優(yōu)化布局后射線密度和射線正交性變大,指標(biāo)E1和指標(biāo)E2也有所改善,矩陣D的性質(zhì)改善,從而所得重建誤差減小。說明優(yōu)化布局有一定的優(yōu)越性,在實際試驗中可以利用這種方案來指導(dǎo)傳感器布局節(jié)省試驗經(jīng)費(fèi)。

圖7 傳感器均勻分布

圖8 傳感器隨機(jī)分布

圖9 傳感器優(yōu)化分布

圖10 不同布局下特征值分布

3 結(jié)論

首先,本文提出了一種地下爆炸超壓場重建的方案,基于這種方案,我們從模型離散化和傳感器優(yōu)化布局兩方面來進(jìn)行爆炸層析成像信息優(yōu)化獲取技術(shù)研究,得出如下結(jié)論:

①不規(guī)則離散模型策略可以滿足爆炸沖擊波在近場、中場、遠(yuǎn)場的不同分辨率要求,抗干擾能力強(qiáng),傳感器利用率高。

②傳感器布局時,盡量使得射線路徑矩陣的秩越大越好、特征值越大越好、射線密度越大越好、射線正交性越大越好、條件數(shù)越小越好。

③重建精度很大程度上依賴于傳感器布局方式,所以,實際試驗中我們要采用優(yōu)化布局算法指導(dǎo)實際傳感器布放。節(jié)省試驗經(jīng)費(fèi)。

[1] 邱艷宇,盧紅標(biāo),蔡力艮. 爆炸沖擊波信號處理方法比較[J]. 爆破,2010,327(1):92-95.

[2] 仲倩. 爆炸源毀傷效應(yīng)測評方法研究[D]. 南京:南京理工大學(xué),2007:7.

[3] 童平. 地震層析成像方法及其應(yīng)用研究[D]. 北京:清華大學(xué),2012.

[4] 李石堅,廖備水,吳健. 面向目標(biāo)跟蹤的傳感器網(wǎng)絡(luò)設(shè)計、實現(xiàn)和布局優(yōu)化[J]. 傳感技術(shù)學(xué)報,2007,20(12):2622-2630.

[5] 劉國華,王振宇,孫堅. 彈性波層析成像及其在土木工程中的應(yīng)用[J]. 土木工程學(xué)報,2003,36(5):76-81.

[6] 成谷. 地震反射走時層析理論與應(yīng)用研究[D]. 上海:同濟(jì)大學(xué),2004:6.

[7] Zolta′n We′ber. Seismic Traveltime Tomography:A Simulated Annealing Approach[J]. Physics of the Earth and Planetary Interiors,2000,119:149-159.

[8] 白苗苗,郭亞麗,王黎明. 基于爆炸超壓場重建的傳感器優(yōu)化布局技術(shù)研究[J]. 傳感技術(shù)學(xué)報,2014,27(7):886-892.

[9] 林大超,白春華. 爆炸地震效應(yīng)[M]. 北京:地質(zhì)出版社,2007.

[10] LanLan Zhen,Ling Wang,Xiuting Wang. A Novel PSO-Inspired Probability-Based Binary Optimization Algorithm[C]//Information Science and Engineering,2008. ISISE’08. International Symposium on. IEEE,2008:248-251.

[11] 高海兵,周馳,高亮. 廣義粒子群優(yōu)化模型[J]. 計算機(jī)學(xué)報,2005,28(12):1980-1987.

[12] 劉洪波,王秀坤,譚國真. 粒子群優(yōu)化算法的收斂性分析及其混沌改進(jìn)算法[J]. 控制與決策,2006,21(6):636-645.

白苗苗(1987-),女,甘肅慶陽人,研究生在讀,研究方向為信號與信息處理、爆炸場重建、反演算法,baimiaomiao641@163.com;

王黎明(1974-),男,教授,博士生導(dǎo)師。主要從事X射線圖像處理、多維信號獲取與處理、無損檢測技術(shù)等領(lǐng)域的教學(xué)和研究工作,wlm@nuc.edu.cn。

A Technology Optimizing the Information about ShallowUnderground Explosion*

BAIMiaomiao,WANGLiming*,GAORui,FANXuejun,WANGKe

(College of Communication and Information Technolgy,North University of China,Taiyuan 030051,China)

A technology of shock wave overpressure space-time field reconstruction about underground explosion was put forward,based on the theory of seismic tomography. In order to improve the underground explosion field tomography information effectively,the paper optimized imaging information from two aspects of model discretization and sensor layout. For discrete model,the irregular grid was adopted according to the attenuation rule of shock wave. For the sensor layout,according to the mathematical underdetermined problem theory and seismic wave propagation law,put forward the index and optimization algorithm for guiding the sensor layout,realize the rational layout of sensors. Simulation results show that,tomography information optimization extraction can maximize the utilization of data,which makes the inversion process is more stable,the inversion error is small.

sensor;optimized layout;seismic tomography;overpressure space-time field;underground explosion;

項目來源:國家自然科學(xué)基金項目(61471325)

2014-11-08 修改日期:2015-03-09

C:7230

10.3969/j.issn.1004-1699.2015.06.014

TP301.6

A

1004-1699(2015)06-0858-06