許金鑫，李慶武，管志強(qiáng)，王肖霖

（1 南京船舶雷達(dá)研究所，南京211106）

（2 河海大學(xué)物聯(lián)網(wǎng)工程學(xué)院，江蘇常州213002）

0 引言

作為研究武器內(nèi)部結(jié)構(gòu)的重要手段，高能閃光X 射線照相技術(shù)可以獲得目標(biāo)內(nèi)爆和演化過程的清晰圖像。高能閃光X 射線照相中目標(biāo)的密度重建屬于典型的高維非線性反演問題，密度重建過程會遇到一些難題，包括數(shù)據(jù)維度高、投影圖像中存在噪聲和模糊等，影響圖像重建的質(zhì)量［1-3］。

由于需要有效地提供重建結(jié)果的不確定度以更好地量化武器的內(nèi)部信息，基于統(tǒng)計(jì)的馬爾科夫鏈蒙特卡羅（Markov Chain Monte Carlo，MCMC）算法逐漸得到了研究，MCMC［4-5］作為一類重要的隨機(jī)采樣算法，能對后驗(yàn)分布提供全面的描述。目前高維數(shù)據(jù)下MCMC 重建算法的研究主要集中在線性模型上，非線性重建方面的研究較少。由于忽略了系統(tǒng)模糊對重建結(jié)果的影響，線性重建結(jié)果在目標(biāo)界面區(qū)域易受模糊的影響。因此需在更符合實(shí)際成像過程的非線性模型上構(gòu)建貝葉斯重建模型，對后驗(yàn)分布進(jìn)行有效采樣并提供相應(yīng)的不確定度分析。高能閃光X 射線投影數(shù)據(jù)維度較高，相比于線性模型非線性成像模型涉及系統(tǒng)模糊的計(jì)算且無法像線性模型一樣推導(dǎo)出目標(biāo)參數(shù)的均值和方差形式從而直接采樣，直接對非線性成像模型進(jìn)行MCMC 求解存在計(jì)算量大以及推理復(fù)雜等問題。

文獻(xiàn)［6］研究了一種先隨機(jī)后優(yōu)化（Random Then Optimizion，RTO）的非線性重建算法，通過求解振蕩的優(yōu)化問題對近似的后驗(yàn)分布進(jìn)行采樣。進(jìn)一步，文獻(xiàn)［7］在RTO 算法的基礎(chǔ)上構(gòu)造非線性分層貝葉斯模型，避免了正則項(xiàng)參數(shù)的人為選取。文獻(xiàn)［8］將RTO 算法擴(kuò)展到非高斯先驗(yàn)，通過變量變換建立非高斯先驗(yàn)與高斯先驗(yàn)?zāi)Ｐ烷g的聯(lián)系。針對RTO 算法在處理高維數(shù)據(jù)時(shí)存在計(jì)算量大的問題，文獻(xiàn)［9］引入一種新的子空間加速策略使得RTO 算法的計(jì)算復(fù)雜度與目標(biāo)參數(shù)的維度成線性比例關(guān)系。在此基礎(chǔ)上，文獻(xiàn)［10］將RTO 和偽邊際MCMC 結(jié)合以提高對模型參數(shù)魯棒的采樣性能。此外，隨機(jī)最大后驗(yàn)（Random Maximum Posterior，rMAP）算法［11］避免了隨機(jī)最大似然（Randomized Maximum Likelihood，RML）［12-13］算法需要評估正向模型二階導(dǎo)數(shù)帶來的昂貴計(jì)算。

RTO 算法及其改進(jìn)方案更多地為了提高樣本統(tǒng)計(jì)效率或計(jì)算效率，在降低樣本統(tǒng)計(jì)誤差方面的研究較少。多級技術(shù)［14-16］和多保真度技術(shù)［17-18］通過分析優(yōu)化問題從而實(shí)現(xiàn)給定替代模型集的最佳分配利用。文獻(xiàn)［17］研究了一種模型管理策略，在誤差和成本方面平衡了高保真和替代模型的模型評估數(shù)量，從而加速估計(jì)具有昂貴計(jì)算成本的高保真度模型的樣本統(tǒng)計(jì)量。該算法中優(yōu)化的目標(biāo)函數(shù)建立在多個(gè)替代模型的樣本統(tǒng)計(jì)差之上，最終得到的樣本融合統(tǒng)計(jì)量相比于單一的高保真度模型估計(jì)的統(tǒng)計(jì)量在精度上的提升并不明顯。

基于上述討論，本文結(jié)合實(shí)際估計(jì)的模糊核，提出一種基于隨機(jī)擾動優(yōu)化和多模型融合的高能閃光X射線圖像非線性重建算法。首先，在非線性正向模型下構(gòu)建分層貝葉斯模型，采用弱信息先驗(yàn)定義與噪聲和先驗(yàn)精度參數(shù)有關(guān)的超參數(shù)以獲得相對于共軛先驗(yàn)更合理的重建結(jié)果。然后，基于隨機(jī)擾動的RTO 算法研究一種適用于高維數(shù)據(jù)的求解模型，結(jié)合雅可比矩陣投影約束該優(yōu)化模型的求解，并設(shè)計(jì)相應(yīng)的提議分布在M-H 框架中進(jìn)行樣本的接受判別計(jì)算，以降低樣本統(tǒng)計(jì)偏差。最后，研究一種新穎的多模型融合方案，有效提高非線性模型樣本估計(jì)的效率和目標(biāo)高密度、邊界區(qū)域重建的準(zhǔn)確度。

1 加速隨機(jī)擾動的非線性分層MCMC 重建

1.1 非線性重建模型及其最小二乘解

1.2 基于弱信息超先驗(yàn)的分層貝葉斯模型

1.3 加速隨機(jī)擾動優(yōu)化求解

式中，|Jx|涉及協(xié)方差矩陣的求逆以及一階導(dǎo)數(shù)間相乘，高維度模型下的計(jì)算仍復(fù)雜。由于正向投影過程是對線吸收系數(shù)矩陣的每一行依次掃描產(chǎn)生投影數(shù)據(jù)，線性模型下相鄰行之間相互獨(dú)立。受模糊核的影響，非線性模型下投影圖像中任意一行數(shù)據(jù)受線吸收系數(shù)矩陣中當(dāng)前行數(shù)據(jù)以及以模糊核半徑覆蓋范圍內(nèi)的上、下行數(shù)據(jù)的影響，此時(shí)?f(x)T?f(x)為對稱矩陣且數(shù)值集中在對角線附近區(qū)域，如圖1所示。較亮區(qū)域數(shù)值非零且以塊狀依次分布在對角線上，稱為“有效矩陣塊”。假設(shè)線吸收系數(shù)矩陣的行、列數(shù)均為M，模糊核大小k×k，則圖1（a）中任意一個(gè)有效矩陣塊的大小為(M/2)×(M/2)（紅圈區(qū)域），非線性模型中任意一行線吸收系數(shù)則對應(yīng)((2k-1)M/2) ×((2k-1)M/2)大小的有效矩陣（紅色正方形區(qū)域，紅圈區(qū)域?yàn)橹行挠行Ь仃噳K）。因此，可通過任意一行線吸收系數(shù)對應(yīng)的有效矩陣來替代完整的一階導(dǎo)數(shù)乘積近似計(jì)算|Jx|，精度矩陣的行、列數(shù)也由原先的(M/2)×M降至(M/2)×(2k-1)，有效減少高維數(shù)據(jù)下接受率的計(jì)算量。

圖1 不同正向模型的一階導(dǎo)數(shù)乘積Fig.1 First derivative products of different forward models

2 最小方差約束的多模型蒙特卡羅估計(jì)

2.1 最小方差約束的融合模型

2.2 融合模型最優(yōu)系數(shù)管理

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

為驗(yàn)證本文非線性重建算法的有效性，在常密度、變密度法國實(shí)驗(yàn)客體（French Test Object，F(xiàn)TO）透射率圖像［22］上進(jìn)行實(shí)驗(yàn)。以相對誤差（Relative Error，RE）［4］、峰值信噪比（Peak Signal-to-Noise Ratio，PSNR）和結(jié)構(gòu)相似度（Structural Similarity Index Measurement，SSIM）［23］作為客觀評價(jià)指標(biāo)來分析算法的重建精度。樣本統(tǒng)計(jì)效率通過積分自相關(guān)時(shí)間（Integral Autocorrelation Time，IACT）［6］來評估，IACT 通過馬爾科夫鏈的自相關(guān)函數(shù)（Autocorrelation Function，ACF）［6］計(jì)算得到，表示獲得一個(gè)獨(dú)立樣本所需的MCMC 步驟數(shù)量。

3.1 提議分布以及雅可比矩陣投影分析

為驗(yàn)證提議分布和雅可比矩陣投影約束共同作用下樣本估計(jì)的可靠性，本小節(jié)在常密度FTO 圖像上計(jì)算樣本估計(jì)的精度和統(tǒng)計(jì)效率，并將RTO-MH 算法［6］得到的樣本估計(jì)作為參考。在高維數(shù)據(jù)下不同的初始值x0會影響LM 算法的收斂速度以及重建結(jié)果的準(zhǔn)確度，因此通過實(shí)驗(yàn)分析初始值設(shè)置對重建性能的影響。定義本文算法在采樣過程中不考慮樣本的接受判別，即統(tǒng)計(jì)所有采樣樣本的算法為RML，對分辨率為100×100 的常密度FTO 透射率圖像進(jìn)行非線性重建。圖2 給出了三種設(shè)置下超參數(shù)λ、δ和目標(biāo)參數(shù)的ACF 及IACT 值，橫、縱坐標(biāo)分別表示滯后時(shí)間和自相關(guān)系數(shù)，其中第一行數(shù)據(jù)對應(yīng)非線性模型下的線性解，即模糊矩陣B為單位矩陣，每次采樣中目標(biāo)參數(shù)的初始值均為0.01。第二行中模糊矩陣B仍為單位矩陣，但當(dāng)前采樣下LM 算法的初始值為上一次采樣得到的樣本值。第三行數(shù)據(jù)則在非線性解下得到，B為標(biāo)準(zhǔn)模糊矩陣，目標(biāo)參數(shù)的初始值均為0.01?？梢钥闯?，在高維數(shù)據(jù)下，本章算法依舊取得了與RTO-MH 算法相近的樣本計(jì)算效率，且優(yōu)于RML 算法。RML 算法由于缺少了樣本的接受判別過程對不穩(wěn)定的樣本值也進(jìn)行了統(tǒng)計(jì)導(dǎo)致樣本統(tǒng)計(jì)效率較低。與RTO-MH 和RML 算法相比，本章算法可得到更低的IACT 值，擁有更高的統(tǒng)計(jì)效率。受目標(biāo)參數(shù)維度較高的影響，相比于噪聲精度參數(shù)，先驗(yàn)精度參數(shù)的采樣樣本具有更高的相關(guān)性。此外，從圖2 中第三行與前兩行結(jié)果對比可知，實(shí)際求解模型的準(zhǔn)確性也會影響目標(biāo)參數(shù)特別是超參數(shù)的統(tǒng)計(jì)效率。

圖2 各算法樣本估計(jì)的ACF 和IACTFig.2 ACF and IACT of sample estimation for each algorithm

圖3 為各算法樣本估計(jì)的RE 值和單次采樣時(shí)間。從圖3（a）可以看出，圖2 中第三行所代表的非線性求解比第一、二行代表的線性求解更加準(zhǔn)確，說明所構(gòu)建模型準(zhǔn)確的重要性。從圖3（b）單次樣本計(jì)算時(shí)間中可以看出，相比于線性求解，非線性求解需要更多的計(jì)算量以及內(nèi)存需求，高效的線性求解也推動了本文多模型融合算法的研究與實(shí)施。從圖3（b）中同一顏色的柱狀圖可以看出，本文算法由于簡化了隨機(jī)擾動的目標(biāo)函數(shù)并對接受率的計(jì)算進(jìn)行優(yōu)化，在計(jì)算效率上相比于RTO-MH 算法得到較明顯的提升，且與RML 算法相近。對比圖3 中灰、綠色柱狀圖可知，對于RTO-MH 和RML 算法，雖然上一次采樣的樣本估計(jì)值作為當(dāng)前采樣中LM 算法初始值的設(shè)置可以得到相比于參數(shù)初始值為0.01 更準(zhǔn)確的重建結(jié)果，但需要更多的計(jì)算時(shí)間。因此，本文LM 算法中目標(biāo)參數(shù)初始值設(shè)置為x0=0.01，適用于隨機(jī)擾動優(yōu)化模型的求解，提議分布和雅可比矩陣投影約束也適用于FTO 透射率圖像的非線性重建。

圖3 各算法樣本估計(jì)的RE 值和重建時(shí)間Fig.3 The RE of sample estimation and reconstructed time by each algorithm

3.2 多模型融合性能分析

定義模糊核1 和模糊核2，其中模糊核1 通過函數(shù)1/(i2+j2+1)產(chǎn)生，i、j分別表示水平、垂直方向上的數(shù)值；模糊核2 通過文獻(xiàn)［23］復(fù)原算法在高能閃光X 射線靜態(tài)投影圖像上估計(jì)得到。本文中模糊核內(nèi)各元素值均通過除以元素之和進(jìn)行歸一化處理。根據(jù)模糊核1、模糊核2 以及均值為0 方差為3 的高斯模糊核正向計(jì)算得到的透射率圖像，對采用標(biāo)準(zhǔn)核和估計(jì)核下的重建結(jié)果進(jìn)行實(shí)驗(yàn)，其中估計(jì)核通過文獻(xiàn)［23］算法計(jì)算得到。

多模型融合前后重建結(jié)果RE 值如圖4所示，其中S 和E 分別表示標(biāo)準(zhǔn)核和估計(jì)核，F(xiàn) 表示使用多模型融合?？梢钥闯觯诓捎脴?biāo)準(zhǔn)核的融合結(jié)果中，如圖中的黑色實(shí)線和虛線所示，低噪聲下多模型融合的效果不如噪聲較嚴(yán)重時(shí)明顯。這主要是因?yàn)樵诘驮肼暻闆r下，由于使用準(zhǔn)確的模糊核，其非線性重建結(jié)果較準(zhǔn)確，特別在無噪聲情況下重建結(jié)果接近于真實(shí)值，此時(shí)多模型融合無法有效發(fā)揮其作用。而在估計(jì)核的融合結(jié)果中，由于估計(jì)的模糊核相比于標(biāo)準(zhǔn)核存在一定誤差，導(dǎo)致在無噪聲情況下重建結(jié)果一定程度偏離真實(shí)值，此時(shí)高密度區(qū)域的非線性和線性重建結(jié)果通過加權(quán)融合可得到更準(zhǔn)確的融合結(jié)果，且隨著噪聲的增加，該融合作用越明顯。值得注意的是，從圖中紅色和黑色實(shí)線可以看出，在噪聲相對較嚴(yán)重的情況下，估計(jì)核下的融合結(jié)果比采用標(biāo)準(zhǔn)核的融合結(jié)果更加準(zhǔn)確，說明多模型融合機(jī)制以及非線性模型中模糊核估計(jì)的有效性。

圖4 多模型融合前后的RE 值Fig.4 The RE involving multi-model fusion

3.3 與現(xiàn)有重建算法對比實(shí)驗(yàn)

本小節(jié)在FTO 圖像上進(jìn)行實(shí)驗(yàn)，分別在常、變密度的FTO 透射率圖像上添加0%、0.25%、0.75%、1.25%等級的高斯噪聲，與現(xiàn)有的基于不確定度分析的重建算法LRIS_Gamma［24］、LRIS_Jeffrey［24］、SPA［25］、S_RTO［26］以及TLE_Gibbs［27］進(jìn)行比較，以驗(yàn)證非線性重建算法的有效性。LRIS_Gamma 和LRIS_Jeffrey 算法基于低秩理論計(jì)算目標(biāo)參數(shù)閉合解，SPA 為基于參數(shù)分裂和增廣的采樣算法，S_RTO 非線性算法在RTO 算法的基礎(chǔ)上近似表示雅可比矩陣的分解因子并引入低秩方法進(jìn)行因子分解，有效降低算法的計(jì)算復(fù)雜度。TLE_Gibbs 為基于兩級不確定度約束的高效MCMC 線性采樣算法。Proposed（T_F）算法為本文的非線性重建算法，其中T_F 表示采用真實(shí)模糊核，且引入多模型融合處理；I_nF 表示非線性模型下的線性解，即模糊矩陣B=I；Proposed 表示本文采用估計(jì)模糊核的非線性重建算法。

為進(jìn)行定量分析，表1 與表2 給出了無加權(quán)中值濾波后處理情況下常、變密度FTO 透射率圖像重建結(jié)果的PSNR 和SSIM 值?？梢钥闯?，由于模型的準(zhǔn)確構(gòu)建非線性重建算法的性能優(yōu)于線性重建算法。除Proposed（T_F）算法使用真實(shí)模糊核外，本文算法（采用估計(jì)核）的非線性重建結(jié)果在所有算法中具有最佳的指標(biāo)值，在噪聲干擾較嚴(yán)重情況下，SSIM 值甚至優(yōu)于Proposed（T_F）算法。此外，隨著噪聲的增加，在1.25%噪聲水平下，本文算法整體重建性能超越Proposed（T_F）算法。

表1 常密度FTO 重建結(jié)果的PSNR 和SSIM 值Table 1 The PSNR and SSIM of reconstructed results of constant density FTO images

表2 變密度FTO 重建結(jié)果的PSNR 和SSIM 值Table 2 The PSNR and SSIM of reconstructed results of variable density FTO images

圖5、6 分別對比了1.25%噪聲水平下變、常密度FTO 透射率圖像的重建結(jié)果。從圖中可以明顯地看出，SPA、LRIS_Gamma 以及LRIS_Jeffrey 算法受模型不準(zhǔn)確的影響重建結(jié)果中仍存在較嚴(yán)重的噪聲。和其他線性算法相比，TLE_Gibbs 重建結(jié)果在常密度區(qū)域相對平滑，但在中垂線區(qū)域連貫性較差。S_RTO 算法的重建結(jié)果中噪聲較少，但由于該算法未對受不合理采樣的超參數(shù)影響而未收斂的目標(biāo)解進(jìn)行有效約束，重建均值中仍存在邊緣模糊的問題。Proposed（I_nF）算法重建結(jié)果的中垂線部分誤差較大，且與線性結(jié)果類似，均存在邊緣模糊問題。Proposed（T_F）以及本文算法（采用估計(jì)的模糊核）的重建結(jié)果與真實(shí)值最為接近，能較好地抑制噪聲，在高密度區(qū)域能夠得到較高精度的重建結(jié)果，且本文算法重建結(jié)果在圖像中垂線附近優(yōu)于Proposed（T_F）算法。

圖5 1.25%噪聲等級下變密度FTO 透射率圖像重建結(jié)果Fig.5 Reconstructed results of FTO transmission image with variable density（1.25% noise level）

圖6 1.25%噪聲等級下常密度FTO 透射率圖像重建結(jié)果Fig.6 Reconstructed results of FTO transmission image with constant density（1.25% noise level）

3.4 非密高能閃光X 射線圖像重建

圖7（a）為4 MeV 能級下的非密高能閃光X 射線靜態(tài)圖像，圖像分辨率為2 048×2 048，每像素對應(yīng)的目標(biāo)尺寸為0.096 mm。藍(lán)色虛線框內(nèi)的實(shí)驗(yàn)?zāi)繕?biāo)為倒置的錐臺，材質(zhì)為錫，錐臺的標(biāo)準(zhǔn)密度為7.3 g/cm3。紅色實(shí)線框中3# ～4#為標(biāo)記塊，材料為錫，采用LM 算法擬合透射率-面密度曲線［28］。受算法采樣效率以及計(jì)算平臺內(nèi)存的限制，線性重建實(shí)驗(yàn)在降采樣后分辨率為350×350 的圖像上開展，非線性重建實(shí)驗(yàn)由于需要消耗更多的內(nèi)存空間在分辨率為200×200 的圖像上開展。重建結(jié)果如圖8所示，可以看出，GPSR 算法重建結(jié)果仍存在中心豎直剖線區(qū)域連貫性較差的問題，LRIS_Jeffrey 和LRIS_Gamma 算法重建結(jié)果中目標(biāo)常密度區(qū)域較平滑，但在錐臺等腰區(qū)域較為模糊。SPA 與TLE_Gibbs 算法重建結(jié)果相近，重建目標(biāo)邊緣較LRIS 算法更清晰，但在錐臺下表面附近的常密度區(qū)域，其重建結(jié)果不如本文非線性算法重建結(jié)果平滑。與線性重建結(jié)果相比，非線性算法能夠較明顯地對圖像背景噪聲進(jìn)行有效抑制。此外，本文非線性算法重建結(jié)果在錐臺等腰區(qū)域具有較好的視覺效果，在統(tǒng)一的濾波算法的作用下六種算法的重建圖像均比較平滑，得益于提出的多模型融合策略，本文非線性算法的重建結(jié)果中常密度區(qū)域的平滑性最好。

圖7 4 MeV 高能閃光X 射線靜態(tài)圖像Fig.7 High energy flash X-ray static image under 4 MeV

圖8 密度重建結(jié)果（4 MeV 靜態(tài)圖像）Fig.8 Results of density reconstruction（4 MeV static image）

為了更加直觀地對比不同算法的重建結(jié)果，圖9 給出了各算法在CO2、CO1 處的密度剖線對比，其中黑色實(shí)線表示標(biāo)準(zhǔn)密度值。可以看出，本文非線性重建結(jié)果的剖線在所有算法中最為平穩(wěn)且更接近標(biāo)準(zhǔn)值，即使在錐臺中垂線區(qū)域，非線性算法的重建結(jié)果相比于線性結(jié)果也更為平穩(wěn)、準(zhǔn)確，在重建精度上具有明顯優(yōu)勢。表3 列出了上述六種重建算法在不同剖線段的密度重建均值及相對誤差。對于錐臺區(qū)域，LRIS_Jeffrey 和LRIS_Gamma 算法重建精度在CO1 剖線段則高于GPSR 和SPA 算法，在CO2 剖線段與SPA 算法接近，說明低秩近似理論更適用于具有單一結(jié)構(gòu)目標(biāo)的重建。本文非線性重建算法的重建精度優(yōu)于線性重建算法，在所有算法中排名第一，其主要原因在于該非線性重建算法進(jìn)行了多物理模型的融合，在方差最小化準(zhǔn)則下合理分配加權(quán)系數(shù)，兼顧線性和非線性解的優(yōu)點(diǎn)，有效提高了非線性模型樣本估計(jì)的準(zhǔn)確性。

表3 各算法密度重建均值（g/cm3）及相對誤差（%）Table 3 Average value（g/cm3）and relative error（%）of density reconstruction of each algorithm

圖9 各算法重建的密度剖線對比Fig.9 Comparison of reconstructed density profiles of each algorithm

4 結(jié)論

本文結(jié)合估計(jì)的模糊核，提出了一種隨機(jī)擾動優(yōu)化和多模型融合的高能閃光X 射線圖像非線性重建算法。該算法從貝葉斯理論的角度考慮高能閃光X 射線圖像非線性重建問題，構(gòu)建分層貝葉斯重建模型?；陔S機(jī)擾動的RTO 算法，結(jié)合雅可比矩陣投影約束求解該非線性優(yōu)化模型，并在M-H 框架中通過對樣本概率密度函數(shù)的近似有效提高了樣本接受判別的計(jì)算效率。最后，在樣本統(tǒng)計(jì)量方差最小化準(zhǔn)則下提出了一種多模型融合方案，逐像素地對非線性和線性重建結(jié)果分配相應(yīng)的加權(quán)系數(shù)，得到融合后更準(zhǔn)確的重建結(jié)果。與現(xiàn)有重建算法的對比實(shí)驗(yàn)證明了本文非線性重建算法的有效性。

由于目標(biāo)參數(shù)維度較高，重建速度還需進(jìn)一步提升。如何引入尺度化機(jī)制減小所需求解的目標(biāo)參數(shù)的維度從而提高非線性重建的計(jì)算效率，并在替代模型較多情況下結(jié)合不確定度對各模型權(quán)重進(jìn)行約束從而獲得更準(zhǔn)確的重建結(jié)果，仍需要進(jìn)一步深入研究。