劉曉佳 李 劍 孫澤鵬 馬明星 魏曉曼

（中北大學(xué)信息探測(cè)與處理山西省重點(diǎn)實(shí)驗(yàn)室 太原 030051）

1 引言

爆炸原理是物體在某極短的時(shí)間范圍內(nèi)，釋放射出來大量熱能量，產(chǎn)生高溫，同時(shí)破壞性極強(qiáng)，難以通過直接觀察的方式研究。隨著戰(zhàn)斗部研制向著智能化，一體化發(fā)展，對(duì)于爆炸場(chǎng)的壓力毀傷能力等參量研究具有重要意義。但是由于傳感器節(jié)點(diǎn)布設(shè)難度大等問題，有限的數(shù)據(jù)無法直接、全面地判斷彈藥的毀傷能力，因此沖擊波場(chǎng)重建對(duì)于彈藥研發(fā)過程和定性驗(yàn)證是至關(guān)重要的。

目前，沖擊波超壓場(chǎng)重建的方法有層析成像法，插值重建法等，白苗苗等［1］針對(duì)爆炸過程中獲得的數(shù)據(jù)較少，對(duì)反演速度和精度進(jìn)行分析，加入先驗(yàn)值進(jìn)行EM 反演；郭亞麗等［2］針對(duì)沖擊波峰值超壓經(jīng)驗(yàn)公式的求解的局限性，利用廣義逆方法進(jìn)行系數(shù)矩陣求解，對(duì)網(wǎng)格單元加重權(quán)，對(duì)沖擊波速度場(chǎng)參數(shù)進(jìn)行反演，再進(jìn)而根據(jù)沖擊波峰值的超壓系數(shù)與沖擊波速度分布的函數(shù)關(guān)系可得到超壓分布，其重建的結(jié)果遠(yuǎn)遠(yuǎn)優(yōu)于經(jīng)驗(yàn)公式；楊志等［3］采對(duì)沖擊波壓力場(chǎng)進(jìn)行插值，對(duì)比多種方法，表明B樣條插值算法效果最好；趙化彬等［4］提出了一種基于非均勻的有理B 樣條（NURBS）“蛛網(wǎng)”插值的場(chǎng)重建方法；張曉光［5］將多元多項(xiàng)式理論應(yīng)用到?jīng)_擊波峰值插值理論計(jì)算中，擬合出波后衰減波形進(jìn)行彈藥沖擊波重建；姚悅［6］利用Biharmonic 樣條曲面插值算法進(jìn)行超壓峰值重建；郭晶［7］提出了一種利用攝像機(jī)對(duì)沖擊波場(chǎng)進(jìn)行拍照記錄得到?jīng)_擊波波陣面的傳播速度對(duì)當(dāng)量進(jìn)行估算，分析爆炸過程。目前科研人員集中研究二維沖擊波場(chǎng)重建方法。由于爆炸沖擊波是各向異性，常規(guī)的二維傳播模型無法有效表征三維空間的傳播特性。因此，三維沖擊波超壓場(chǎng)重建是大威力彈藥毀傷威力評(píng)估時(shí)亟待解決的瓶頸問題。

針對(duì)上述問題，本文在二維走時(shí)層析成像模型的基礎(chǔ)上，結(jié)合射線追蹤理論，建立了三維沖擊波超壓場(chǎng)重建模型，由于三維超壓場(chǎng)重建問題是大型稀疏矩陣求解問題，且傳統(tǒng)重建算法對(duì)邊緣噪聲的抑制效果較差，在層析成像的基礎(chǔ)上結(jié)合了壓縮感知理論中的字典學(xué)習(xí)（Dictionary Learning）聯(lián)合全變分（Total Variation）正則化，提出了TV-DL 迭代修正算法，為三維沖擊波超壓場(chǎng)重建及實(shí)現(xiàn)其全方位可視化提供了一種新方法，提高了沖擊波超壓場(chǎng)的重建精度。

2 三維沖擊波超壓場(chǎng)重建原理

炸藥在空氣中快速爆炸，瞬時(shí)所產(chǎn)生爆炸的爆炸高壓一般可達(dá)1010Pa，空氣爆炸的瞬時(shí)初始?jí)毫Υ蠹s為105Pa，爆炸產(chǎn)物極速膨脹，對(duì)空氣進(jìn)行壓縮，形成爆炸沖擊波，由于沖擊波傳播速度較快（瞬時(shí)10-6s）可以認(rèn)為是沿直射線方向傳播的，走時(shí)可以通過速度和路徑長(zhǎng)度之間的函數(shù)來表示，如下式：

其中，t為沖擊波到達(dá)探測(cè)點(diǎn)的時(shí)間，即為走時(shí)；L為沖擊波到達(dá)測(cè)點(diǎn)的傳播路徑；v 為沖擊波波速；S 為慢度；r為對(duì)應(yīng)射線。

假設(shè)將重建范圍劃分為M×N×K=J個(gè)網(wǎng)格，對(duì)式（1）進(jìn)行離散化處理，如圖1所示。

圖1 三維離散化示意圖

圖2 求解三維投影矩陣流程圖

上圖為三維空間中第i 條射線路徑上射線穿過各網(wǎng)格的投影值之和，其函數(shù)表達(dá)式如下式：

其中，I 為射線總數(shù)，即測(cè)點(diǎn)個(gè)數(shù)；J 為離散網(wǎng)格總數(shù)；ti為第i 條射線的走時(shí)，即信號(hào)到達(dá)第i 個(gè)探測(cè)點(diǎn)前的時(shí)間；aij為第i 條投影射線穿過第j 個(gè)投影網(wǎng)格的投影射線長(zhǎng)度，即系數(shù)投影矩陣。

上式可總結(jié)為矩陣形式：

其中，A 為投影矩陣向量，其大小為I×(M×N×K)的二維矩陣，S 為慢度向量矩陣，其大小為(M×N×K)×1 的一維矩陣，T 為走時(shí)向量矩陣，其大小為I×1的一維矩陣。

在建立三維沖擊波超壓場(chǎng)模型的基礎(chǔ)上，矩陣A 的求解是構(gòu)建三維模型的關(guān)鍵，投影矩陣A 實(shí)際上就是I 條射線在三維空間中所有網(wǎng)格的長(zhǎng)度信息，本文基于射線追蹤最短路徑的原理計(jì)算矩陣A，具體實(shí)現(xiàn)流程如下。

在求射線與網(wǎng)格交點(diǎn)的時(shí)候，本文采用先得到射線與網(wǎng)格交點(diǎn)一軸的坐標(biāo)點(diǎn)，根據(jù)比例關(guān)系求得其他兩軸坐標(biāo)點(diǎn)，計(jì)算步驟如式（4）。

其中，steplen 為步長(zhǎng)，lx為區(qū)域長(zhǎng)度，Bx為傳感器x坐標(biāo)，By、Bz同理為傳感器的y、z坐標(biāo)。

本文投影矩陣A 為二維矩陣，索引函數(shù)就是將求得的離散網(wǎng)格內(nèi)的數(shù)據(jù)還原到對(duì)應(yīng)三維空間位置處。

3 三維重建模型迭代修正算法

三維沖擊波超壓場(chǎng)重建過程中，重建區(qū)域大，但是爆炸過程破壞性極大，導(dǎo)致可布設(shè)的測(cè)點(diǎn)少，三維重建問題是一個(gè)典型的病態(tài)稀疏矩陣的求解問題。傳統(tǒng)的迭代重建算法重建誤差較大，針對(duì)這個(gè)問題，利用壓縮感知在稀疏矩陣約束方面的優(yōu)勢(shì)［10～12］，本文在建立三維沖擊波超壓場(chǎng)模型（3）基礎(chǔ)上，利用字典學(xué)習(xí)對(duì)圖像的稀疏表征優(yōu)勢(shì)［14～15］，全變分（Total Variation，TV）最小化將重建問題轉(zhuǎn)變?yōu)榍蠼庾钚≈档膯栴}，TV 最小化理論對(duì)圖像邊緣紋理的保真特性［13］，實(shí)現(xiàn)對(duì)層析模型的正則約束。將TV 最小化與字典學(xué)習(xí)相結(jié)合，重新構(gòu)建了重建模型：

其中A為投影矩陣，S為重建矩陣，T為投影數(shù)據(jù)。第一項(xiàng)為重建數(shù)據(jù)的保真項(xiàng)，第二項(xiàng)為正則項(xiàng)，第三項(xiàng)為字典學(xué)習(xí)項(xiàng)。μ為保真項(xiàng)系數(shù)，λ和β為正則項(xiàng)系數(shù)。

求解目標(biāo)函數(shù)（5）步驟如下：

步驟一：固定D和αj，更新重建矩陣S，目標(biāo)函數(shù)為

為求解式（6），引入輔助變量dx，dy，bx，by，得到目標(biāo)函數(shù)（9）：

步驟二：固定S，更新D和αj，對(duì)應(yīng)的目標(biāo)函數(shù)如下：

首先進(jìn)行字典更新，固定αj更新D，使用K-SVD 方法從重建數(shù)據(jù)S中學(xué)習(xí)字典。然后進(jìn)行稀疏表示，固定D更新αj，采用OMP 方法來更新稀疏表示系數(shù)αj。

4 實(shí)驗(yàn)仿真與結(jié)果分析

為驗(yàn)證本文的三維沖擊波超壓場(chǎng)重建模型的可行性，采用ART 算法、SART 算法和TV-DL 方法進(jìn)行三維仿真重建實(shí)驗(yàn)，比較上述三種算法各自的重建仿真效果以及實(shí)現(xiàn)沖擊波三維重建的可視化。

設(shè)置大小為24m×24m×12m 的矩形區(qū)域?yàn)橹亟▍^(qū)域，區(qū)域被均勻地劃分為1m×1m×1m 的網(wǎng)格，重建模型中添加了7.5%的隨機(jī)噪聲，炸點(diǎn)位置坐標(biāo)為（0，0，0），仿真以整個(gè)重建區(qū)域的四分之一為例，待重建區(qū)域大小為12m×12m×12m，被劃分為1728 個(gè)網(wǎng)格，布設(shè)了36 個(gè)傳感器（即36 個(gè)測(cè)點(diǎn)），圖3為傳感器布設(shè)，圖4為各射線與三維網(wǎng)格的交點(diǎn)圖。在圖4中，共有36 條不同顏色的點(diǎn)線，每個(gè)顏色的點(diǎn)線代表著一條射線（即爆炸點(diǎn)到一個(gè)測(cè)點(diǎn)的傳播路徑）與途經(jīng)網(wǎng)格的交點(diǎn)。

圖3 傳感器分布圖

圖4 射線與網(wǎng)格交點(diǎn)

仿真實(shí)驗(yàn)的TNT 當(dāng)量為5kg，重建模型選用Henrych 沖擊波超壓經(jīng)驗(yàn)公式，根據(jù)沖擊波峰值超壓與速度的關(guān)系將模型轉(zhuǎn)化為沖擊波波速，將對(duì)沖擊波超壓場(chǎng)重建問題轉(zhuǎn)化為對(duì)沖擊波速度場(chǎng)重建問題。迭代終止修正條件若只是滿足給定的迭代次數(shù)或終止修正值與初始值的誤差精度，迭代結(jié)果會(huì)出現(xiàn)局部收斂的情況，本文的迭代終止條件采用修正后的到達(dá)時(shí)間與初始到達(dá)時(shí)間的估計(jì)誤差小于0.1%，圖5為重建初始模型。

圖5 三維沖擊波重建初始模型

三種方法分別迭代300 次進(jìn)行反演重建，三維重建結(jié)果和Z=0平面的速度值分布如圖6所示。

圖6 三種算法分別迭代300次的重建結(jié)果

圖6（a）～（b）分別為ART 迭代算法的三維重建結(jié)果和Z=0平面的重建速度分布；圖6（c）～（d）分別為SART迭代算法的三維重建結(jié)果和Z=0平面的重建速度分布；圖6（e）～（f）分別為TV-DL 迭代修正算法的三維重建結(jié)果和Z=0平面的重建速度分布。

由圖6可以看出，與ART 算法和SART 算法的三維重建結(jié)果相比，在重建區(qū)域邊緣采用本文方法進(jìn)行的三維重建結(jié)果與采用傳統(tǒng)迭代算法的三維重建結(jié)果相比更為光滑，為了更直觀地觀察三種方法對(duì)于邊緣噪聲的抑制情況，本文觀察了三種重建結(jié)果的Z=0 平面的速度值分布，觀察比較圖6（b）、（d）、（f）可以明顯看出，在有效抑制邊緣噪聲問題上，TV正則化具有較大優(yōu)勢(shì)。

為了客觀闡述以上三種算法對(duì)于重建區(qū)域的重建精度，本文將采用相對(duì)誤差（RE）和均方根誤差（RMSE）作為重要評(píng)價(jià)參數(shù)，RE一般用來客觀評(píng)價(jià)三種方法的重建速度模型與初始速度模型在一個(gè)速度重建區(qū)域范圍內(nèi)的各網(wǎng)格之間的實(shí)際重建的效果，RMSE 一般用來客觀評(píng)價(jià)一個(gè)整體區(qū)域內(nèi)的實(shí)際重建結(jié)果與重建初始模型結(jié)果的偏離的程度，值越小，表明重建結(jié)果越接近實(shí)際情況。三種算法重建結(jié)果的RE曲線見圖7。

圖7 每個(gè)網(wǎng)格內(nèi)速度的相對(duì)誤差

圖7為三種迭代算法重建后的速度與初始速度在重建區(qū)域內(nèi)的1728 個(gè)網(wǎng)格中的相對(duì)誤差，可以明顯看出，在重建區(qū)域內(nèi)，相對(duì)誤差曲線呈周期性衰減，究其原因是重建網(wǎng)格的劃分是逐層排序共記12 層。由圖7可看出曲線衰減大致分為12 段，爆炸沖擊波場(chǎng)在近場(chǎng)區(qū)超壓值變化極快，與遠(yuǎn)場(chǎng)相比，重建誤差大，因此，誤差曲線呈分段衰減的趨勢(shì)。在圖7中，TV-DL算法的相對(duì)誤差基本能保持在3%左右，與ART 算法相比，降低了5%左右，與SART算法相比，降低了3%左右。進(jìn)一步用均方根誤差來作為重建后整體誤差，三種算法的RMSE 值如表1所示。

表1 重建結(jié)果的均方根誤差

已知重建區(qū)域整體的速度范圍大約在400m/s～3000m/s，由上表可知，三種算法整體的速度場(chǎng)重建的總誤差依次為45.81m/s、38.79m/s、9.85m/s，即本文方法的整體重建效果顯著優(yōu)于ART 算法和SART 算法。通過分析綜合運(yùn)用以上兩個(gè)參數(shù)，可以進(jìn)一步直接證明本文提出的方法在區(qū)域重建方面更具優(yōu)勢(shì)。從數(shù)據(jù)整體重建的效果來看，TV-DL方法所獲得的結(jié)果與真實(shí)的模型小，三維重建的誤差遠(yuǎn)比其他兩種算法小，綜合分析來看，在三維峰值超壓場(chǎng)重建模型的實(shí)際開發(fā)應(yīng)用層面上，TV-DL方法的獲得三維重建模型效果最優(yōu)。

5 結(jié)語

為得到爆炸沖擊波三維空間分布情況，并達(dá)到良好的可視化效果，首先本文深入研究了沖擊波二維走時(shí)層析成像的原理，在二維走時(shí)層析成像模型的基礎(chǔ)上，根據(jù)射線追蹤理論，構(gòu)建起了三維沖擊波峰值超壓強(qiáng)場(chǎng)的重建模型；然后，運(yùn)用了壓縮感知算法在稀疏約束方面的優(yōu)點(diǎn)，給出了一種基于TV-DL 沖擊波超壓強(qiáng)場(chǎng)迭代修正算法；開展了模擬試驗(yàn)，并對(duì)比了傳統(tǒng)重構(gòu)算法ART、SART 和TV-DL 的重構(gòu)結(jié)果。對(duì)試驗(yàn)結(jié)論分析表明，本文首次提出的三維射線走時(shí)層析成像方法已經(jīng)能夠初步實(shí)現(xiàn)三維沖擊波超壓場(chǎng)重建。TV-DL 修正代數(shù)重建精度的相對(duì)誤差基本保持在3%左右，與ART 算法相比，降低了5%左右，與SART 算法相比，降低了3%左右，本文僅對(duì)沖擊波三維重建進(jìn)行了仿真實(shí)驗(yàn)，但由于實(shí)際傳感器布設(shè)困難，后續(xù)會(huì)繼續(xù)研究在滿足重建精度的前提下減少傳感器個(gè)數(shù)，實(shí)現(xiàn)沖擊波三維重建。