一種基于加權(quán)多準則的火炮爆炸場中高溫火球三維溫度場的迭代重建算法

劉曉蔚 祝海江* 劉興旺 李小春 李卓明

(1.北京化工大學(xué) 信息科學(xué)與技術(shù)學(xué)院， 北京 100029； 2.中國電子科技集團公司 第四十一研究所， 青島 266555)

引 言

爆炸爆溫是指炸藥爆炸時放出的能量將爆炸產(chǎn)物加熱到的最高溫度，它是評估炸藥性能的一個重要指標。爆炸過程具有破壞力強、電磁干擾力強、溫度高、溫度動態(tài)范圍大等特點，爆炸場溫度數(shù)據(jù)的采集與測量難度高，目前尚沒有非常有效的測試技術(shù)及設(shè)備，因此，基于爆炸溫度場的非接觸測量重建方法研究具有重要的實際應(yīng)用價值。近年來，爆炸溫度場的測量與重建得到了許多研究人員的關(guān)注[1-5]。郭曉杰等[1]結(jié)合高速攝像技術(shù)和輻射測溫法實現(xiàn)了瞬態(tài)高溫場的測試。范航[4]研究了爆炸溫度場的代數(shù)重建(algebraic reconstruction techniques，ART)算法、聯(lián)合代數(shù)重建(simultaneous algebraic reconstruction technique, SART)算法與乘法代數(shù)重建(multiplicative algebraic reconstruction technology，MART)算法仿真，并給出了一種改進的多準則迭代重構(gòu)(multi criteria iterative reconstruction，MCIR)算法進行溫度場分布的仿真重建。徐維錚等[5]根據(jù)自主開發(fā)的數(shù)值計算程序開展密閉空間內(nèi)炸藥爆炸溫度場的數(shù)值計算，初步探討了內(nèi)爆炸溫度場演化過程及其分布規(guī)律。

隨著計算機技術(shù)和數(shù)字圖像技術(shù)的進步，基于圖像的溫度場三維重建方法研究得到了迅速發(fā)展。溫度場三維重建方法主要結(jié)合在醫(yī)學(xué)診斷、無損檢測等領(lǐng)域常用的計算機斷層成像技術(shù)來實現(xiàn)基于火焰圖像的溫度場三維重建。Gordon等[6]提出了經(jīng)典的代數(shù)重建技術(shù)，然而經(jīng)典ART方法由于其邊緣效應(yīng)，重建質(zhì)量不高。Gilbert[7]提出了一種優(yōu)化的ART方法以降低對噪聲的敏感度，也即同時迭代重建(simultaneous iterative reconstruction technology，SIRT)算法，該方法中的像素值在每次迭代中均使用平均像素值進行更新，從而達到更好的降噪效果，但收斂速度降低。 Andersen等[8]提出了ART的另一種變形方法，該方法通過二次優(yōu)化SIRT方法有效地提高了ART算法的收斂速度，即聯(lián)合代數(shù)重建算法。SART是一種在圖像三維重建中應(yīng)用較為廣泛的迭代重建算法，許多研究人員將其應(yīng)用于醫(yī)學(xué)計算機斷層掃描(CT)圖像的重建[9-14]、微波圖像重建[15-17]以及電力系統(tǒng)的電容層析成像[11]等方面。冀東江等[9]對SART算法進行改進，提出了在分母加上懲罰項的懲罰SART算法，其懲罰項的思想是減小畸變處的迭代步長以減弱畸變。喬全邦等[10]提出了一種基于模糊熵的自適應(yīng)聯(lián)合代數(shù)算法應(yīng)用于CT圖像重建；該方法采用模糊熵對圖像進行邊緣提取，并根據(jù)邊緣一致性原則構(gòu)造單調(diào)遞增函數(shù)，將該單調(diào)遞增函數(shù)定義為迭代步長的松弛算子，從而進行圖像重建，它解決了重建圖像邊緣模糊的問題，還可以抑制振鈴效應(yīng)。針對電容層析成像中圖像重建算法的非線性和病態(tài)問題，宋亞杰等[18]將SART算法應(yīng)用于電容層析成像，提高了圖像重建速度，并且重建圖像的邊緣信息保真度較高。Han 等[15]提出了一種在SART算法中插入非線性因子的改進方法，該方法在特定的84個投影點的情況下可以取得較好的收斂效果，但在不滿足特定條件的情況下算法不收斂。Donato等[11]通過使用雙邊正則化濾波器對CT圖像進行濾波以改進圖像的重建質(zhì)量，并在顯卡(GPU)上對乳腺樣本進行了圖像重建。Tiwari等[12]使用各向異性擴散作為正則化方法改進了SART的迭代方法，提高了CT圖像的重建質(zhì)量。Aprilliyani等[16]和Ramdani等[17]將SART算法應(yīng)用于微波系統(tǒng)中腫瘤的檢測或癌癥的無創(chuàng)層析成像系統(tǒng)的開發(fā)上，利用3 GHz的微波信號，通過特定位置的發(fā)射端和接收端來采集數(shù)據(jù)，成功重構(gòu)了圖像。Zhang等[13]改進了SART算法中投影矩陣的計算方法，使其計算速度得到顯著提高。Sabah等[14]討論了簡單反投影、濾波反投影和迭代算法的性能參數(shù)，并認為與迭代重建中的其他方法相比，SART算法改善了圖像質(zhì)量，減少了重建誤差并產(chǎn)生了條紋痕跡。

盡管SART的改進算法對于抑制圖像畸變、提高重建質(zhì)量起到了一定的作用，但傳統(tǒng)的畸變校正方法只對邊緣區(qū)域的抑制作用明顯，在校正畸變時并不遵循某一最優(yōu)準則，并且由于爆炸場高溫數(shù)據(jù)獲取困難的特點，在投影角度較少的情況下，重建圖像依然會出現(xiàn)明顯的邊緣效應(yīng)，并且邊緣效應(yīng)會使重建圖像和原始圖像的輻射積分誤差趨近于零，導(dǎo)致迭代無法進行，使得重建圖像中間區(qū)域失真嚴重，無法達到重建要求。因此，針對爆炸場溫度瞬時性的特點，本文借鑒懲罰思想并融合均勻性準則對SART算法的迭代步長進行改進，所采用的變迭代步長方法能有效地抑制畸變和邊緣效應(yīng)，提高圖像重建質(zhì)量。

1 改進的SART算法

1.1 懲罰SART算法

SART算法是溫度場重建中一種常見的迭代重建算法[6]，該算法的基本思想是在同一截面內(nèi)對通過某一像素的不同方向的射線校正值進行累加來校正該像素值，其迭代公式為

(1)

SART算法可以使理論計算值逼近實際測量值，但是當(dāng)投影數(shù)據(jù)含有噪聲時，重建圖像中的噪聲隨迭代次數(shù)不斷累加，導(dǎo)致圖像重建效果不佳。為了解決這一問題，冀東江等[9]通過在權(quán)值中增加懲罰項來改進SART算法，其迭代公式為

(2)

其中

(3)

1.2 算法的改進

懲罰SART算法的思想是通過犧牲迭代速度，減小校正值以減弱畸變，提高重建質(zhì)量，但該方法對所有區(qū)域均使用統(tǒng)一的迭代步長，導(dǎo)致在邊緣區(qū)域畸變依然明顯。因此，本文采用變迭代步長的懲罰方式，以有效抑制邊緣區(qū)域的畸變，在SART算法滿足的最小二乘準則中加入均勻性準則，提高重建質(zhì)量。

考慮到實際測量時相鄰像素間的灰度值是接近的，即參數(shù)變化較為平緩，所以針對像素間的關(guān)系提出均勻性準則，即

(4)

迭代系數(shù)λ體現(xiàn)了校正值的權(quán)重，其值一般取0～1。λ值較大時，由于迭代較快，畸變也會被放大；λ值太小，則會使迭代速度過慢。采用變迭代步長是希望在無畸變處λ能取較大值，不影響迭代速度，而在有畸變處λ能取較小值，以減小迭代步長，抑制畸變。為實現(xiàn)這種效果，本文改進了迭代系數(shù)表達式，如式(5)所示。

(5)

其中

(6)

圖1 3×3的懲罰區(qū)域示意圖Fig.1 3×3 penalty area diagram

將式(5)代入式(1)可得到變迭代步長SART算法公式為

(7)

綜上所述，基于變迭代步長SART算法的步驟可歸納如下。

1)將原圖像分別旋轉(zhuǎn)0°、45°、90°并向x軸投影，用得到的輻射積分圖像模擬測量值。

2)根據(jù)旋轉(zhuǎn)角度構(gòu)建相應(yīng)的投影矩陣W。

3)給定計算初值Fo。

4)確定懲罰區(qū)域長度o與懲罰系數(shù)α、β。

7)判斷式(7)計算結(jié)果是否滿足終止條件，如滿足，退出迭代；否則執(zhí)行步驟5)進行下一次迭代。

2 實驗結(jié)果與分析

2.1 實驗驗證

為檢驗本文改進的算法，采用如式(8)所示的函數(shù)模擬三維溫度場的分布情況以進行測試[4]。

(8)

所得函數(shù)圖像如圖2所示，可以看出圖形有兩個正高斯峰和1個負高斯峰，可用來模擬復(fù)雜情況下的三維溫度場分布。

圖2 模擬多峰函數(shù)圖像Fig.2 Analog image multimodal function

在仿真實驗中，采用0°、45°、90°這3個角度的輻射積分圖像進行重建，得到的輻射積分圖像如圖3所示。

圖3 輻射積分圖像Fig.3 Radiation integral images

本文采用輻射積分誤差e和峰值誤差m作為評價標準，其計算公式為

e=‖P-WF‖

(9)

m=1-max(F)

(10)

2.1.1第一組仿真實驗

在第一組仿真實驗中，將懲罰項系數(shù)設(shè)定為α=0.000 01，β=15，懲罰區(qū)域分別設(shè)定為5×5、11×11和31×31。迭代算法終止條件為輻射積分誤差e<5。

從圖4的仿真結(jié)果可以看出，變迭代步長的SART算法重建圖像的正峰和負峰都在函數(shù)所描述位置附近，方向均正確，值為0的區(qū)域畸變明顯消除。

圖4 第一組仿真實驗的三維重建圖像Fig.4 Three-dimensional reconstruction images for the first set of simulation experiments

由圖5的輻射積分誤差曲線圖可以看出懲罰區(qū)域為11×11的結(jié)果相比于懲罰區(qū)域5×5的結(jié)果收斂更快，只經(jīng)過21次迭代即可將誤差降低至5以下，且波峰和波谷的范圍更小，形狀更加突出；而當(dāng)懲罰區(qū)域擴大至31×31后，迭代加快，但邊緣效應(yīng)也更加明顯，逐漸失去抑制畸變的能力。

圖5 第一組仿真實驗的重建輻射積分誤差曲線Fig.5 Plots for reconstruction of radiation integral error for the first set of simulation experiments

由圖6的峰值誤差曲線可以看出，懲罰區(qū)域為11×11的峰值誤差在0附近，這說明重建函數(shù)最大值接近原函數(shù)。懲罰區(qū)域為5×5的峰值誤差出現(xiàn)負值，說明重建函數(shù)最大值超過了原函數(shù)最大值，這是因為懲罰區(qū)域過小，圖像對正常區(qū)域進行了過度的懲罰，抑制了正常的迭代過程，為減小輻射積分誤差而在峰值區(qū)域出現(xiàn)迭代值過大的情況。懲罰區(qū)域為31×31的峰值誤差在0.1附近，這是因為懲罰區(qū)域選擇過大導(dǎo)致峰值誤差較大，從而無明顯的懲罰作用。

圖6 第一組仿真實驗的峰值誤差曲線Fig.6 Plots for peak error for the first set of simulation experiments

由上可知，懲罰區(qū)域選擇過大，會失去懲罰作用，導(dǎo)致圖像重建效果不佳；懲罰區(qū)域選擇過小，會因為過度懲罰增加迭代次數(shù)而使峰值過大；而選擇合適的懲罰區(qū)域，則可以改善重建圖像的質(zhì)量。因此，針對不同圖像選擇合適的懲罰系數(shù)和懲罰區(qū)域可達到兼顧重建質(zhì)量和收斂速度的效果。

2.1.2第二組仿真實驗

在第二組仿真實驗中，將懲罰區(qū)域設(shè)定為11×11，改變懲罰項系數(shù)進行實驗，即當(dāng)α=0.000 01時，β依次取值0.5、1、15，所得仿真結(jié)果如圖7～9所示。

圖7 第二組仿真實驗的三維重建圖像Fig.7 Three-dimensional reconstruction images for the second set of simulation experiments

圖8 第二組仿真實驗的重建輻射積分誤差曲線Fig.8 Plots for reconstruction of radiation integral error for the second set of simulation experiments

圖9 第二組仿真實驗的峰值誤差曲線Fig.9 Plots for peak error for the second set of simulation experiments

2.1.3第三組仿真實驗

在第三組仿真實驗中，將懲罰區(qū)域設(shè)定為11×11，β取值為15時，α依次取值0.000 01、0.001、0.1，所得仿真結(jié)果如圖10～12所示。

圖10 第三組仿真實驗的三維重建圖像Fig.10 Three-dimensional reconstruction images for the third set of simulation experiments

圖11 第三組仿真實驗的重建輻射積分誤差曲線Fig.11 Plots for reconstruction of radiation integral error for the third set of simulation experiments

圖12 第三組仿真實驗的峰值誤差曲線Fig.12 Plots for peak error for the third set of simulation experiments

2.2 不同重建算法的比較

本文還對比了變迭代步長SART算法、SART算法和懲罰SART算法迭代30次的重建圖像，其中變迭代步長SART算法和懲罰SART算法的懲罰區(qū)域均為11×11，懲罰項系數(shù)均為α=0.000 01，β=15，所得結(jié)果見圖13～15。

圖13 3種不同算法的三維重建圖像Fig.13 Three-dimensional reconstruction images for three different algorithms

圖14 3種不同算法的重建輻射積分誤差曲線Fig.14 Plots for reconstruction of radiation integral error for three different algorithms

圖15 3種不同算法的峰值誤差曲線Fig.15 Plots for peak error for three different algorithms

由圖13可以看出，SART算法和懲罰SART算法對圖像的還原度較高，收斂快，但存在畸變，特別是在值為0的邊緣區(qū)域畸變較為明顯。變迭代步長的SART算法不僅能還原圖像，還能明顯消除值為0的區(qū)域的畸變，使峰值區(qū)域更突出。

由圖14、15可知，在3種算法的輻射積分誤差相近的情況下，SART及懲罰SART算法的峰值誤差均在0.1以上，其值較大。變迭代步長的SART算法具有較高的收斂速度，在30次迭代后輻射積分誤差較小，峰值誤差在0.1以下，接近原函數(shù)峰值。

3 結(jié)論

針對傳統(tǒng)SART算法、懲罰SART算法在圖像重建過程中存在的誤差較大，尤其是邊緣畸變突出的問題，本文提出了一種改進的SART算法。本文所提改進算法在原有的最小二乘準則上融合了均勻性準則，即對畸變的判定與校正采用均勻性準則，對不滿足這一準則的區(qū)域進行懲罰以抑制畸變，對滿足這一準則的區(qū)域不作處理。該準則不僅適用于邊緣區(qū)域，且對圖像所有區(qū)域均適用。傳統(tǒng)的畸變校正方法只對邊緣區(qū)域的畸變抑制作用明顯，且在校正畸變時并不遵循某一準則，沒有特定的判定標準。相比之下，本文方法中對畸變懲罰系數(shù)可以人為設(shè)定，根據(jù)重建圖像的不同設(shè)定合適的判定系數(shù)，可以兼顧重建速度和質(zhì)量。實驗結(jié)果表明本文方法優(yōu)勢在于能有效抑制圖像的邊緣效應(yīng)產(chǎn)生的畸變，提高重建質(zhì)量；并且其懲罰項的系數(shù)和懲罰區(qū)域可針對不同的圖像自由選擇，選擇合適的懲罰系數(shù)和懲罰區(qū)域可以獲得更好的重建質(zhì)量和較高的收斂速度。但該改進方法是以犧牲迭代步長為代價的，因此優(yōu)化耗時較長，對于大圖像精細重建需要的時間較久，難以及時得到重建結(jié)果。