大孔徑厚針孔成像數(shù)值模擬研究

姚志明，段寶軍，馬繼明，宋 巖，嚴(yán)維鵬，宋顧周，韓長材

(西北核技術(shù)研究所 強(qiáng)脈沖輻射環(huán)境模擬與效應(yīng)國家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710024)

高能X/γ射線或中子具有較強(qiáng)的穿透能力，沒有合適的聚焦材料，無法由透鏡直接聚焦成像。1969年起，美國洛斯阿拉莫斯實(shí)驗(yàn)室開始發(fā)展的針孔成像技術(shù)[1]是射線源空間分布圖像測量的常用手段。針孔成像技術(shù)屬于直接成像，像與物一一對應(yīng)，所見即所得。為使針孔成像過程具有較高的空間分辨率，通常將針孔直徑設(shè)計(jì)得較小，典型值為0.5 mm。小孔徑限制了針孔對射線的接收效率，無法給出低強(qiáng)度射線源的信息。

相對于針孔成像，編碼孔成像技術(shù)擁有更大的射線接收角。編碼成像是每個物點(diǎn)在接收面上形成1個編碼孔投影像，不同物點(diǎn)產(chǎn)生的像相互錯開并疊加，物與像不一一對應(yīng)，需對編碼像進(jìn)行解碼復(fù)原。目前，適用于高能射線編碼的半影孔成像技術(shù)在慣性約束聚變的內(nèi)爆靶丸診斷中得到大量應(yīng)用[2-4]。半影孔成像將源圖像信息編碼到半影區(qū)域中，孔徑大于源尺寸，適用于百μm級尺寸的內(nèi)爆靶丸診斷。

對于cm級尺寸的射線源，例如FTO(French test object)模型[5]，其內(nèi)部高原子序數(shù)材料在高溫流體狀態(tài)下發(fā)生裂變，膨脹過程直徑可達(dá)10 cm，此時，半影孔孔徑需大于10 cm，編碼圖像接收屏尺寸需達(dá)30 cm以上。受測試系統(tǒng)空間布局限制，半影孔成像技術(shù)很難應(yīng)用于大尺寸射線源的診斷。若射線源劑量率低于傳統(tǒng)針孔成像系統(tǒng)靈敏度響應(yīng)下限，就無法給出射線源信息。本文提出介于傳統(tǒng)針孔成像和半影孔成像之間的大孔徑厚針孔成像方法，基于Geant4蒙特卡羅模擬軟件和Matlab圖像處理軟件建立大孔徑厚針孔編碼成像和圖像復(fù)原過程的數(shù)值模擬方法，并優(yōu)化復(fù)原算法的輸入?yún)?shù)，通過與傳統(tǒng)針孔成像結(jié)果的比較，評估大孔徑厚針孔成像方法的可行性。

1 大孔徑厚針孔成像原理

與半影孔成像相似，大孔徑厚針孔成像分為兩個步驟，如圖1所示。第1步是編碼成像過程。射線源經(jīng)過大孔徑厚針孔成像，由退化圖像記錄裝置記錄為數(shù)字化圖像。針孔孔徑可達(dá)mm～cm量級，與傳統(tǒng)針孔成像中使用的百μm量級的厚針孔相比[6]，射線接收效率可提高2～4個量級。理想針孔成像的空間分辨由式(1)描述，M為成像系統(tǒng)放大倍數(shù)，隨著針孔孔徑d的增大，空間分辨Δr(物面上可被分辨的兩個點(diǎn)源之間的最近距離)線性增加，射線源經(jīng)大孔徑厚針孔后得到的是模糊的像。

圖1 大孔徑厚針孔成像原理示意圖Fig.1 Schematic of large thick aperture imaging

(1)

第2步是圖像復(fù)原。數(shù)字化圖像輸入計(jì)算機(jī)后采用算法程序運(yùn)算解碼，得到空間分辨改善后的清晰圖像。圖像復(fù)原過程具有病態(tài)性，無法得到與真實(shí)圖像相同的結(jié)果，但可根據(jù)針孔點(diǎn)擴(kuò)散函數(shù)(PSF, point spread function)等先驗(yàn)知識，獲取盡可能逼近真實(shí)圖像的近似解。

2 圖像復(fù)原算法

大孔徑厚針孔成像系統(tǒng)物面上不同位置點(diǎn)源的PSF不同，屬于空間移變系統(tǒng)[4]。當(dāng)滿足近軸條件時，PSF可認(rèn)為是空間移不變的，能采用逆濾波、Wiener濾波和Lucy-Richardson等空間移不變圖像復(fù)原算法復(fù)原。退化圖像i(x,y)、源圖像o(x,y)、針孔PSFh(x,y)和附加噪聲n(x,y)之間的關(guān)系可由式(2)描述：

i(x,y)=h(x,y)*o(x,y)+n(x,y)

(2)

兩邊進(jìn)行Fourier變換，得到頻譜形式下的表達(dá)式(式(3))，其中，H(u,v)為調(diào)制傳遞函數(shù)，I(u,v)、O(u,v)、N(u,v)分別為退化圖像、源圖像和附加噪聲的頻譜。

I(u,v)=H(u,v)·O(u,v)+N(u,v)

(3)

2.1 理想逆濾波

(4)

2.2 Wiener濾波

(5)

其中，Sn(u,v)、So(u,v)分別為噪聲n(x,y)和真實(shí)圖像o(x,y)的功率譜，兩者的比值為NSR(noise to signal power ratio)。

2.3 Lucy-Richardson算法

Lucy-Richardson算法[8]是1種迭代方法，簡稱L-R算法。利用貝葉斯公式結(jié)合極大似然估計(jì)可得到迭代方程，如式(6)所示。其中，k為迭代次數(shù)。迭代過程中使用的h(x,y)是在近軸條件下中心軸線位置的PSF。

(6)

2.4 圖像復(fù)原效果評估

為定量比較不同復(fù)原算法的復(fù)原效果優(yōu)劣和尋找迭代復(fù)原算法的最優(yōu)迭代次數(shù)，通常引入均方根誤差σ[9]，表達(dá)式為：

(7)

其中，M和N為圖像沿x軸和y軸的像素?cái)?shù)。σ越小，圖像復(fù)原結(jié)果在均方根誤差意義上越接近源圖像，復(fù)原效果越好。

3 蒙特卡羅模擬與圖像復(fù)原

3.1 大孔徑厚針孔成像過程模擬

蒙特卡羅方法能逼真地描述射線與物質(zhì)相互作用的物理過程，Geant4軟件包含了p、e-、n和γ射線等多種粒子的物理過程。本文采用Geant4模擬大孔徑厚針孔的編碼成像過程。射線源為1.25 MeV的γ射線，形狀尺寸為6.4 cm×6.4 cm的字母A，射線源距離針孔中心1 m，退化圖像記錄屏距離針孔中心1 m。計(jì)算程序中源粒子發(fā)射過程采用偏倚抽樣技巧，只在稍大于針孔接收立體角的范圍內(nèi)發(fā)射射線，可提高計(jì)算效率，縮短計(jì)算時間。圖像記錄過程對射線的能量和種類進(jìn)行判斷，排除次級γ射線和其他射線干擾，只記錄直穿γ射線圖像。圖2a為源圖像，圖2b～d為退化圖像。為便于比較，模擬了孔徑為0.5 mm的傳統(tǒng)小孔徑厚針孔成像，如圖2e、f所示。

a——源圖像；b——5 mm針孔退化圖像；c——10 mm針孔退化圖像；d——15 mm針孔退化圖像；e——0.5 mm小孔徑針孔成像；f——0.5 mm小孔徑針孔管道因子校正圖像圖2 源圖像與針孔成像模擬結(jié)果Fig.2 Source image and simulation image of aperture

0.5 mm小孔徑針孔成像與源圖像形狀基本一致，由于厚針孔管道因子效應(yīng)，圖像中心區(qū)域成像效率較邊緣的高，管道因子校正是厚針孔成像的常用方法[6]，管道因子fp隨像素與圖像中心距離r的變化曲線如圖3所示，圖2f為校正圖像；5 mm針孔退化圖像基本保留了字母A的輪廓形狀信息，但與0.5 mm孔徑針孔成像相比，中心三角形輪廓已變得模糊；10 mm針孔退化圖像中幾乎無法識別源形狀為字母A；15 mm針孔退化圖像源形狀信息全部損失，僅能觀察到模糊的三角形亮斑。隨著針孔孔徑的增大，成像系統(tǒng)空間分辨率變差，源圖像退化嚴(yán)重，結(jié)合復(fù)原算法恢復(fù)源圖像信息是必要的。

圖3 管道因子曲線Fig.3 Curve of fp vs r

3.2 圖像復(fù)原

Matlab軟件中包含豐富的圖像處理函數(shù)，deconvwnr函數(shù)可實(shí)現(xiàn)逆濾波和Wiener濾波，輸入?yún)?shù)包括退化圖像、PSF和NSR。deconvlucy函數(shù)可實(shí)現(xiàn)L-R復(fù)原算法，輸入?yún)?shù)包括退化圖像、PSF和k。不同孔徑針孔軸線的PSF由蒙特卡羅模擬計(jì)算得到，如圖4所示。

圖4 5、10、15 mm孔徑針孔PSFFig.4 PSF of 5, 10 and 15 mm apertures

1) 理想逆濾波

采用deconvwnr函數(shù)，當(dāng)NSR設(shè)置為0時即為理想逆濾波。復(fù)原結(jié)果如圖5所示?？煽闯?，理想逆濾波算法會引起噪聲放大，難以獲得清晰的復(fù)原圖像。

2) Wiener濾波

采用deconvwnr函數(shù)，改變輸入?yún)?shù)NSR的大小，圖像復(fù)原解與源圖像之間的σ隨之變化，如圖6所示。針孔孔徑越小，σ的極小值越小，圖像復(fù)原解越接近源圖像。針孔孔徑越大，σ取極小值時的NSR越小，即真實(shí)圖像(信號)與噪聲的功率譜之比越大。σ取極小值時的復(fù)原結(jié)果如圖7所示，相應(yīng)的σ列于表1。其中σ0.5為0.5 mm孔徑傳統(tǒng)針孔成像結(jié)果與源圖像的均方根誤差，為1.03×10-4。Wiener濾波算法可復(fù)原得到清晰的字母A的圖像，但圖像周圍存在明顯的偽影。隨著針孔孔徑的增大，偽影現(xiàn)象加重。

圖5 理想逆濾波圖像復(fù)原解Fig.5 Reconstruction result using inverse filtering method

圖6 Wiener濾波σ隨NSR的變化Fig.6 σ of Wiener filtering vs NSR

圖7 Wiener濾波圖像復(fù)原解Fig.7 Reconstruction result using Wiener filtering method

3) L-R算法

采用deconvlucy函數(shù)，改變k的大小，圖像復(fù)原解與源圖像之間的σ隨之變化，如圖8所示。針孔孔徑越小，σ的極小值越小，復(fù)原圖像越接近源圖像。針孔孔徑越大，σ為極小值所需的迭代次數(shù)越多，迭代收斂速度越慢。σ取極小值時的復(fù)原結(jié)果如圖9所示，相應(yīng)的σ列于表1。針孔孔徑為5、10、15 mm條件下，L-R算法均能復(fù)原得到清晰的字母A圖像，復(fù)原圖像中沒有明顯的偽影。

表1 圖像復(fù)原解的均方根誤差Table 1 σ of reconstruction results

圖8 L-R算法σ隨迭代次數(shù)的變化Fig.8 σ of L-R algorithm vs number of interation

圖9 L-R算法圖像復(fù)原解Fig.9 Reconstruction result using Lucy-Richardson method

3種復(fù)原算法相比，理想逆濾波算法會引起噪聲放大，難以獲得清晰的復(fù)原圖像。Wiener濾波和L-R算法均可獲得清晰的源圖像，但Wiener濾波算法偽影現(xiàn)象較嚴(yán)重。與之相比，L-R算法復(fù)原圖像中無明顯偽影，且相同孔徑條件下圖像復(fù)原解的σ極小值更小。因此，L-R算法的圖像復(fù)原效果最好。

與0.5 mm孔徑傳統(tǒng)針孔成像相比，結(jié)合L-R算法的5、10、15 mm大孔徑厚針孔成像獲取的圖像σ極小值分別為前者的1.15、1.21、1.26倍。說明在均方根誤差意義上，大孔徑厚針孔成像結(jié)果與傳統(tǒng)針孔成像接近。圖9和圖2f直觀比較可發(fā)現(xiàn)，圖9中字母A底部部分細(xì)節(jié)信息更清晰，初步驗(yàn)證了大孔徑厚針孔成像方法具有可行性。

4 結(jié)語

為解決cm級尺寸的輻射源空間分布診斷問題，本文提出了介于傳統(tǒng)針孔成像和半影孔成像之間的大孔徑厚針孔成像方法。采用蒙特卡羅模擬和圖像處理軟件實(shí)現(xiàn)了編碼成像過程和圖像復(fù)原過程的數(shù)值模擬。模擬結(jié)果表明，大孔徑厚針孔成像能獲得與傳統(tǒng)針孔成像接近的源圖像測量結(jié)果，初步驗(yàn)證了方法的可行性。

實(shí)際測量的脈沖輻射源通常具有能譜復(fù)雜、尺寸和射線定向發(fā)射強(qiáng)度連續(xù)變化的特點(diǎn)，針對具體的測量對象，優(yōu)化大孔徑厚針孔結(jié)構(gòu)，研究更高重建精度的非線性圖像復(fù)原算法是需要進(jìn)一步深入研究的問題。