姚志明，段寶軍，馬繼明，宋 巖，嚴維鵬，宋顧周，韓長材

(西北核技術研究所 強脈沖輻射環(huán)境模擬與效應國家重點實驗室，陜西 西安 710024)

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

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

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

1 大孔徑厚針孔成像原理

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

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

(1)

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

2 圖像復原算法

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

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

(2)

兩邊進行Fourier變換，得到頻譜形式下的表達式(式(3))，其中，H(u,v)為調制傳遞函數，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)和真實圖像o(x,y)的功率譜，兩者的比值為NSR(noise to signal power ratio)。

2.3 Lucy-Richardson算法

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

(6)

2.4 圖像復原效果評估

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

(7)

其中，M和N為圖像沿x軸和y軸的像素數。σ越小，圖像復原結果在均方根誤差意義上越接近源圖像，復原效果越好。

3 蒙特卡羅模擬與圖像復原

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

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

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

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

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

3.2 圖像復原

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

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

1) 理想逆濾波

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

2) Wiener濾波

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

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

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

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

3) L-R算法

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

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

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

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

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

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

4 結語

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

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