王　鳳，白原齊

(中北大學(xué)，山西 太原 030051)

?

基于OSEM算法的沖擊波超壓場重建技術(shù)

王鳳，白原齊

(中北大學(xué)，山西 太原 030051)

摘要：利用走時(shí)層析成像方法重建水中爆炸沖擊波超壓場，結(jié)合沖擊波傳輸特性來反演速度場完成沖擊波超壓場的重建。試驗(yàn)中爆炸點(diǎn)和傳感器數(shù)目不足給重建過程帶來困難，于是文中提出了一種基于先驗(yàn)信息的OSEM反演算法。首先建立數(shù)學(xué)模型，再運(yùn)用MATLAB進(jìn)行仿真實(shí)驗(yàn)，最后對仿真結(jié)果進(jìn)行分析。實(shí)驗(yàn)表明，本文方案相比傳統(tǒng)重建方法迭代次數(shù)更少、重建精度更高。說明該算法能有效解算這種欠定重建問題。在實(shí)際爆炸中實(shí)用性較高，可指導(dǎo)野外施工。

關(guān)鍵詞：走時(shí)層析成像；水中爆炸；超壓場重建

水下爆炸是毀傷艦船、潛艇以及大型魚雷等水中航行器的一條重要途徑，對艦船目標(biāo)在水下爆炸作用下毀傷效應(yīng)的研究歷來受到各海洋強(qiáng)國的重視。對水中爆炸沖擊波超壓場重建的研究在評估水中爆炸的威力、控制爆炸的危害等爆炸安全技術(shù)方面意義重大[1，2]。

目前我們一般用水下壓力傳感器進(jìn)行試驗(yàn)，測量出相應(yīng)的沖擊波超壓，由于該類傳感器防水性好、穩(wěn)定性好、對系統(tǒng)響應(yīng)頻率高等，而能滿足這種要求的壓力傳感器價(jià)格昂貴，為實(shí)際試驗(yàn)增加了成本費(fèi)用。試驗(yàn)所用的傳感器價(jià)格昂貴，傳感器數(shù)目越少，試驗(yàn)成本越低，傳感器數(shù)目越多，反演精度越高，同時(shí)試驗(yàn)成本也越高。

本文采用走時(shí)層析的方法重建水中爆炸超壓場，此類工程問題必然導(dǎo)致重建數(shù)據(jù)的不完全。運(yùn)用合適的迭代重建算法，使得可用最少數(shù)目的傳感器重建最高精度的沖擊波超壓場。

針對這種不完全數(shù)據(jù)重建問題，在算法研究方面，被提出了一些實(shí)用的迭代重建的算法及它們的改進(jìn)算法，如SIRT算法、SART算法[4，5]等，這兩種算法能在一定程度上提高重建精度。但其重建速度很慢。由于在實(shí)際爆炸工程應(yīng)用中，傳感器布局及測量誤差導(dǎo)致數(shù)據(jù)不完整，使得反演算法的不穩(wěn)定性和誤差大大增加。本文提出了一種基于先驗(yàn)信息的OSEM反演算法，來解決這種不完全數(shù)據(jù)走時(shí)的欠定問題。

1沖擊波超壓場重建方案

水中爆炸沖擊波峰值超壓場重建通過三個(gè)步驟，第一步對各個(gè)傳感器陣元采集獲得的沖擊波信號進(jìn)行分析及特征提取，得到?jīng)_擊波到達(dá)時(shí)間，在這里定義為走時(shí)，第二步采用走時(shí)層析成像[6-8]方法重建待測區(qū)域內(nèi)的沖擊波速度，第三步根據(jù)速度與沖擊波峰值超壓的關(guān)系，最后將網(wǎng)格的速度值運(yùn)用相應(yīng)公式轉(zhuǎn)換成超壓值，來重建出整個(gè)二維平面上的沖擊波超壓值。

根據(jù)沖擊波的傳輸特性，將重建區(qū)域劃分成規(guī)則網(wǎng)格。本文中假設(shè)一條射線對應(yīng)一個(gè)傳感器并假設(shè)沖擊波不是沿網(wǎng)格單元的邊界傳播。定義沖擊波從激勵源到達(dá)探測器的時(shí)間為走時(shí)，綜上可得沖擊波在傳輸?shù)倪^程中，其走時(shí)是速度和幾何路徑的函數(shù)，如式(1)所示：

(1)

式中，r為射線；v為速度；s為慢度；L為積分路徑。 將上式離散化，根據(jù)測試區(qū)域網(wǎng)格劃分情況，對于第i條射線有：

(2)

式中，ti為第i條射線的走時(shí)，即沖擊波波陣面到達(dá)傳感器的時(shí)間；dij為第i條射線穿過第j個(gè)網(wǎng)格的射線長度；sj為第j個(gè)網(wǎng)格中的慢度；M為射線數(shù)；N為網(wǎng)格數(shù)。

將上式(2)寫成矩陣形式為：

DS=T.

(3)

其中，P=(p1，p2，…，pM)′為各條射線走時(shí)的M維列向量；X=(x1，x2，…，xN)′為待求離散單元慢度值，為N維未知的列向量；D為M×N階稀疏矩陣，其元素為dij。

矩陣方程組(3)中，因?yàn)槊總€(gè)網(wǎng)格單元只有一少部分射線穿過，所以得到的方程式中的矩陣D是一個(gè)大型的稀疏矩陣，該矩陣中非零元素比例很大，它的不完全投影使得問題欠定，同時(shí)測量過程存在一定的誤差，這些誤差決定了不相容性。方程組(3)采用OSEM算法求解。能有效改善該問題。

2重建算法及重建誤差

2.1聯(lián)合迭代算法SIRT

聯(lián)合迭代重建算法SIRT，是比較常用的迭代重建算法。其公式為：

(4)

迭代步驟分為：對向量x賦初值x=x(0)。

修正x=x(k)(k=0,1,2,…)：

1) 投影：

(5)

即計(jì)算射線1到射線I的線積分估計(jì)值。

2) 反投影：計(jì)算第j個(gè)體素的修正值，按體素進(jìn)行逐個(gè)修正

(6)

每迭代一輪，方程組更新一次。

3) 重復(fù)2)，xk收斂時(shí)停止。

2.2聯(lián)合代數(shù)重建算法SART

(7)

其中，λ是0<λ<2的一個(gè)松弛因子。

具體步驟為：

1) 對圖像向量x賦初值x=x(0)；

2) 選擇投影方向。對x=x(k)(k=0,1,2,…)進(jìn)行修正：

第一步，投影：計(jì)算該投影方向角度下的所有射線的線積分估計(jì)值，

(8)

第二步，反投影：計(jì)算J各提體素的修正值，按體素進(jìn)行逐個(gè)修正

(9)

循環(huán)到完成一輪迭代。

2.3重建誤差

如下為第i個(gè)網(wǎng)格單元的速度相對誤差：

(10)

其中，v′j為真實(shí)值，vj為重建的速度值。

定義網(wǎng)格平均相對誤差為：

(11)

速度最大相對誤差

(12)

定義第i條射線走時(shí)相對誤差為：

(13)

相應(yīng)地走時(shí)平均相對誤差：

(14)

從而得到走時(shí)最大相對誤差為：

(15)

3OSEM重建算法

為了加快收斂速度，減少運(yùn)算時(shí)間，提高重建精度，人們提出了很多快速算法，有序子集最大期望值法是很有應(yīng)用前景的一種快速迭代重建算法，它是在最大似然期望法(Maximum Likelihood Expext at ion-max-imization，ML-EM)的基礎(chǔ)上發(fā)展起來的。傳統(tǒng)的ML-EM方法計(jì)算式為：

(16)

設(shè)定重建圖像大小為n×n，360個(gè)投影角度，每個(gè)角度下n條投影射線。將投影數(shù)據(jù)按投影角度劃分成T個(gè)有序子集。 OSEM算法的迭代公式如下：Sl表示第l個(gè)子集，l=1,2,…,T。

(17)

3.1建立數(shù)學(xué)模型

建立仿真模型是一個(gè)自由介質(zhì)的沖擊波速度分布模型，假設(shè)根據(jù)沖擊波傳輸特性重建區(qū)域均勻?qū)ΨQ，所以取四分之一區(qū)域研究。如圖1所示，激勵源位于重建區(qū)域中心位置，探測器遍布于整個(gè)邊界并且盡可能多方位覆蓋被測區(qū)域。先把該區(qū)域劃分為5×5的網(wǎng)格單元，使它的分辨率為1 m，傳感器分布在爆炸點(diǎn)的下邊界和左邊界，本文用到了13個(gè)傳感器放在邊界，多余的傳感器放置在內(nèi)部，用來采集信息作為先驗(yàn)信息。

圖1　重建區(qū)域網(wǎng)格均勻離散圖

3.2仿真結(jié)果及分析

基于以上模型參數(shù)，下面主要從平均誤差和計(jì)算速率(迭代次數(shù))兩個(gè)方面來分析。

(a) 慢度平均誤差與迭代次數(shù)的變化關(guān)系

(b) 慢度最大誤差與迭代次數(shù)變化關(guān)系

(a) 理想假設(shè)模型

(b) SIRT算法重建結(jié)果

(c) SART算法重建結(jié)果

(d) OSEM算法重建結(jié)果

迭代次數(shù)OSEMSARTSIRT316783收斂精度(%)OSEMSARTSIRT6.3448.93913.04

通過以上分析可知：在沖擊波場重建中，創(chuàng)痛算法和本文算法的精度的變化趨勢都達(dá)到穩(wěn)定趨勢，OSEM算法的迭代次數(shù)最少，即該算法的重建速度最快；OSEM算法的收斂精度最高，即用該算法重建誤差最小。所以本文選用OSEM算法以最高的精度和最快的速度解決欠定問題。

3.3對OSEM算法進(jìn)行改進(jìn)

本文通過加入先驗(yàn)信息，即采集重建區(qū)域內(nèi)部傳感器的信息作為先驗(yàn)信息，重建爆炸場。一般先驗(yàn)信息有三種: 第一種是能得模型參數(shù)大概分布特點(diǎn)；第二種是由先驗(yàn)認(rèn)識得部分參數(shù)的一定取值范圍；第三種是能獲得少部分參數(shù)的精確值。這里采用第三種，直接把采集到的信息作為特定參數(shù)的精確值，將其加入到方程(3)。實(shí)際試驗(yàn)應(yīng)用中，我們把反演區(qū)域進(jìn)行個(gè)別點(diǎn)取樣，然后測得該位置波速作為第3種先驗(yàn)信息，相當(dāng)于在線性方程式(3)下增加約束方程組

WS=S′ .

(18)

式中，W為k×n維矩陣，k為被約束的參數(shù)個(gè)數(shù)，且wi=[000…100]，設(shè)被約束參數(shù)的位置為1，其他為0。s′為采樣點(diǎn)測出的觀測值，于是式(3)改為：

(19)

式中，D0是無樣點(diǎn)先驗(yàn)信息的距離矩陣，T0是走時(shí)矩陣。將上式改寫為DS=T再進(jìn)行計(jì)算。

3.3.1仿真結(jié)果及分析

圖4　本文方法的重建平均誤差的變化

圖5　本文方法的各網(wǎng)格絕對誤差

圖6　本文方法的OSEM算法重建結(jié)果

由圖4、圖5可知，隨著采樣點(diǎn)由3個(gè)增加到第6個(gè)，反演誤差逐漸減小。加入3個(gè)先驗(yàn)值時(shí)反演誤差接近5%，加入4個(gè)先驗(yàn)值時(shí)和再加一個(gè)是誤差變化不大，都接近3%，加入6個(gè)先驗(yàn)值時(shí)誤差減小到低于2%。

綜上所述，表明采樣點(diǎn)的增加伴隨著重建誤差的減小，但并表明采樣點(diǎn)數(shù)越多實(shí)驗(yàn)精度越大，因此本文的方案可以明顯地提高算法精度。但是采樣點(diǎn)數(shù)的增加意味著傳感器數(shù)目的增加，而水下傳感器成本高試驗(yàn)費(fèi)用大。實(shí)際試驗(yàn)中需要綜合考慮。圖6顯示了基于先驗(yàn)信息的OSEM算法重建結(jié)果，相比前面介紹的方法要更接近理論模型。

4結(jié)論

本文針對爆炸中不完全數(shù)據(jù)的超壓場重建問題，提出了一種基于先驗(yàn)信息的OSEM反演算法。設(shè)定初始模型后經(jīng)過MATLAB對該算法仿真，主要分析反演誤差和迭代次數(shù)兩個(gè)方面，即得到反演精度和反演速率。OSEM算法在解這種不完全數(shù)據(jù)的欠定問題中的誤差較小。對此算法進(jìn)行改進(jìn)，即加入采樣點(diǎn)作為第三種先驗(yàn)信息，發(fā)現(xiàn)加入先驗(yàn)值后的重建算法速度更快精度更高。

OSEM算法相對較穩(wěn)定，試驗(yàn)表明該方法可以獲得比較好的反演結(jié)果。但總的來說它的計(jì)算速度比較慢，并且需要加較大的存儲空間。本文的研究結(jié)果為在單一激勵源層析成像中根據(jù)要求的精度和運(yùn)算速度選擇反演算法提供了有用的參考。

參考文獻(xiàn)

[1]葉序雙，顧文彬，劉文華，等.水下爆炸研究現(xiàn)狀[J].工程爆破,1999,5(1):84-87.

[2]李磊，馮順山.水下爆炸對艦船結(jié)構(gòu)毀傷效應(yīng)的研究現(xiàn)狀與展望[J].艦船科學(xué)技術(shù),2008,30(3):26-30.

[3]Snay H G.The Scaling of Underwater Explosion Phenomena[R].AD-271468,Jun.1961.

[4]Charles L.Byrne.Applied Iterative Methods[M].A K Peters/CRC Press,2007.

[5]Tanabe K.Projection Method for Solving a Singular System of linear Equation and ITS Applications[J].Numer Math，1971,17:203-214.

[6]王家映.地球物理反演理論[M].北京:高等教育出版社,2011.

[7]GUO Cheng-hao,Liu Feng-yu.Model and Analysis of Adaptive Grid Computing[J].Computer Engineering,2008,34(8):117-119.

[8]黃靚.混凝土超聲層析成像的理論方法和實(shí)驗(yàn)研究[D].湖南:湖南大學(xué),2008.

Reconstruction Technology of Blast Wave Overpressure Field on OSEM Algorithm

Wang Feng, Bai Yuanqi

(NorthuniversityofChina,TaiyuanShanxi030051,China)

Abstract:In this paper, the technique of blast wave field reconstruction based on tomography is studied. Overpressure field is reconstructed by inverting the velocity field in the process of shock wave transmission. Since the reconstruction process is difficult due to insufficient excitation sources and detectors, an OSEM algorithm is proposed based on priori information. Appropriate models are constructed using the proposed methods, and a simulation example is put forward at last. The result reveals that compared with the traditional methods, this method has higher precision and converges faster. It also shows the validity and practicality of the developed algorithm in solving the problem of incomplete data reconstruction.

Key words:traveltime tomography; underwater explosions; reconstruction of overpressure field

中圖分類號：TP301.6

文獻(xiàn)標(biāo)識碼：A

文章編號：1674- 4578(2016)01- 0015- 04

作者簡介：王鳳(1990- )，女，湖北仙桃人，碩士研究生，主要從事信號處理、重建類算法。

收稿日期：2015-10-09