復雜快中子源編碼成像高效模擬方法研究

張鳳娜,張天奎,王東明,鄭向陽,胡華四

(1.西安交通大學能源與動力工程學院,710049,西安;2.中國工程物理研究院激光聚變研究中心,621900,四川綿陽;3.環境保護部核與輻射安全中心,100082,北京)

?

復雜快中子源編碼成像高效模擬方法研究

張鳳娜1,張天奎2,王東明1,鄭向陽3,胡華四1

(1.西安交通大學能源與動力工程學院,710049,西安;2.中國工程物理研究院激光聚變研究中心,621900,四川綿陽;3.環境保護部核與輻射安全中心,100082,北京)

為了解決快中子編碼成像中復雜源圖像定義與模擬效率的問題,研究了基于深度優先合并方法的復雜源圖像定義方法。采用柵元合并技術來減小源定義所需柵元數目,同時保證蒙特卡洛模擬中源抽樣的效率,從而提高復雜源圖像編碼成像模擬效率。模擬源為二值E字母時,源合并后模擬所得編碼像的計算結果最小誤差較源合并前沒有變化且滿足統計要求,計算時間則減少且為源合并前模擬時間的1/5;對16、64和256灰度階E字母源進行了編碼像的模擬計算,模擬結果的最小誤差小于1%,符合重建研究的需要;采用3種重建算法對“西安交通大學?；铡睆碗s二值源的定義和模擬進行了源區重建,進而驗證了基于深度優先合并的源定義方法的正確性。該方法可望為增進聚變源區所歷復雜過程診斷的適應性提供一種切實可行的技術途徑。

快中子編碼成像;深度優先合并;復雜灰度源;蒙特卡洛模擬;圖像重建

核聚變產物主要包括中子和α粒子、質子等帶電粒子及光子。聚變熱核反應區的形狀反映了驅動的對稱性、流體力學的不穩定性及輻射的燒蝕等許多物理因素,是判斷點火成功與否的關鍵物理參量[1]。與X射線成像相比,聚變反應產生的粒子可以更直接地反映聚變等離子區域的形狀和熱核反應的燃燒對稱性。產生于聚變高密度等離子體區域的快中子平均能量在14.1 MeV,其具有很高的穿透能力,可以從聚變壓縮區域穿出,因此采用以快中子為對象的編碼成像技術,通過對源區進行編碼成像,由編碼圖像經重建可獲得源圖像,這種重建源圖像可以清晰地反映聚變壓縮區域的尺寸、形狀和均一性等特征[2-5]。實驗研究表明,聚變的壓縮區域形狀和分布具有非對稱、非規則的特點,反映在源圖像(表征中子發射分布)上則呈現為多灰度階的復雜形狀分布[6]。然而,編碼過程模擬中多灰度階復雜源的建模和模擬過程效率較低,因此復雜源圖像研究需要發展高效的建模技術。

建立復雜形狀灰度源圖像的編碼模擬方法,對中子成像系統參數設計與重建算法研究具有重要意義[7]。重建算法建立的目的在于對聚變等離子區域進行診斷,以反映真實的聚變物理過程。等離子區域的形狀和中子強度分布是聚變物理過程中十分重要的信息,這要求重建算法的重建結果必須能夠很好地反映這些信息,其中中子強度分布反映在重建源圖像上就是像素點的灰度值。因此,成像模擬的難點就在于需要實現復雜形狀且灰度值很高的源圖像模擬。MCNP是美國Los Alamos國家實驗室應用理論物理部的Monte Carlo小組研制的蒙特卡洛程序。目前,MCNP中的復雜形狀面源是采用很薄的體源實現的,圖像像素點灰度值表示體源的厚度,灰度越高,厚度越大,相同面積(每個像素點是大小相同的正方形)下中子抽樣也就更多,這符合高像素灰度、高中子發射強度的對應關系。在MCNP輸入卡片中,如果要實現復雜形狀二值源圖像或灰度階數很高的灰度源圖像的模擬,則源柵元的定義就需要很多個基本單元體(一般采用六面體)進行疊加,這一方面導致柵元定義語句很長,超過MCNP柵元定義語句1 000個字(word)的長度限制,另一方面柵元定義的分散性使得模擬計算的效率很低。

為了解決復雜源圖像定義與模擬效率的問題,本文研究了基于深度優先合并方法的復雜源圖像定義技術。通過柵元合并技術,減小源定義所需柵元數,同時保證蒙特卡洛模擬中源抽樣的效率,從而提高復雜源圖像編碼成像模擬的效率。

1 中子編碼成像系統模型

聚變快中子編碼成像模擬是采用MCNP[8-9]模擬實現的。MCNP主要應用于復雜三維幾何結構中粒子輸運的計算,可模擬光子、中子、中子-光子耦合及光子-電子耦合等的輸運問題。中子編碼成像系統模型見圖1,模擬所用源為氘氚(DT)聚變中子,峰值能量為14.1 MeV,源區視場直徑為10 mm。成像系統的放大倍數為5,編碼孔采用60 mm長的鎢孔[10-12],圖像探測器采用塑料閃爍光纖BCF-10組成的陣列[13],長度為100 mm,排布為199×199,單根光纖的直徑為500 μm。在模擬中,整個成像系統置于空氣環境中。圖像探測器在閃爍光纖陣列的光纖芯部,用其記錄中子在閃爍光纖陣列中的能量沉積。經過上述模型模擬獲得的編碼圖像考慮了源區快中子經過編碼孔的編碼效應,也考慮了中子與閃爍光纖的作用過程。

圖1 中子編碼成像系統模型

中子編碼成像中普遍采用2種MCNP編碼成像模型:基于重復結構的能量沉積模型[14-16];基于FIR(flux image radiograph)的點通量探測模型[16-17]。采用能量沉積計數的優點是模擬過程比較全面,但是效率很低;采用FIR計數,由于不需要閃爍光纖陣列作為計數器,因此計算效率高,但無法計算到達像面的中子在閃爍光纖陣列中能量沉積的過程。在本編碼成像系統中,雖然中子與閃爍光纖陣列材料作用過程中的串擾等因素[18]也會對編碼圖像造成一定的影響,但是總體而言孔的編碼是主要的編碼過程,所以FIR計數器獲得的編碼圖像可以很好地反映真實的編碼圖像,尤其是可以計算不同編碼孔材料的散射效應和透射效應[16]。本文旨在研究復雜源圖像的模擬方法,因此為保證計算效率,模擬中采用FIR點通量探測模擬模型。

2 復雜源圖像定義方法

首先通過對MCNP的修改與重編譯,使柵元定義長度由1 000個字長增加為100 000個字長,這樣對應源定義的柵元可以包含的六面體數便由64個增加為7 692個,能滿足像素數為100×100源面的定義需要。在此基礎上,為提高計算效率,參考醫學領域的體元合并算法[19],開展了源定義的柵元內六面體的合并研究。

2.1 源柵元六面體合并的基本定義

源柵元定義中六面體合并的算法,是將灰度值作為特征值,合并相同特征值(特征值大于0)的六面體,減少源定義柵元的六面體數目,從而提升MCNP模擬的效率。

為了實現源定義的六面體與源圖像的像素對應關系,給出三對坐標(系)的數學定義。

定義1 六面體所屬的空間坐標系與像素所屬的平面坐標系。

空間坐標系見圖2a,其中y軸方向的尺寸與源圖像像素點的灰度值成正比,而六面體合并的特征值為源圖像的灰度值。另外,面源定義由體源代替,所以在y方向并無六面體的多層分布。因此,空間上六面體的合并可簡化為面源圖像上以灰度值為特征值的像素點合并。像素平面坐標系見圖2b,其中y軸垂直紙面向里,為了方便調用像素點,定義i-j坐標系,取(1,1)在圖像左上角,并設源圖像的像素點數為N×N。

(a)空間坐標系 (b)平面坐標系 圖2 六面體和像素所屬坐標系

定義2 六面體坐標與像素點坐標。

結合圖2,六面體中心在空間的坐標簡稱為六面體中心坐標(x,y,z),像素點在二維圖像中的位置為像素點的坐標(i,j)。由于像素點的灰度值g(i,j)與六面體y方向的尺寸成正比,所以兩組坐標的轉換關系見式(1)和式(2)。

當源圖像單邊的像素數N為奇數時有

(1)

式中:Δs為六面體在x、z方向上劃分的最小單元尺寸(μm);k為像素點的灰度值與六面體在y方向尺寸的轉換比例尺寸(μm)。

當N為偶數時有

(2)

為了與模擬尺度相對應,式(1)和式(2)的常數取值為Δs=100 μm,k=0.01 μm。取k值很小的目的在于,即使像素點的灰度值g(i,j)很大,也能保證y較小,可以滿足體源代替面源模擬的要求。

定義3 像素合并區域坐標與對應的六面體合并區域坐標。

將像素合并區域(長方形形狀)對角線上的2個頂點的像素坐標稱作該合并區域的坐標(i1,j1,i2,j2),其中i1j2,在該區域中像素點的灰度值相同,均為g(i1,j1)。采用式(3)、式(4)的坐標變換,在MCNP中定義的六面體合并區域的坐標為(x1,y1,z1,x2,y2,z2)。

當源圖像單邊的像素數N為奇數時有

(3)

當N為偶數時有

(4)

比較式(1)和式(3)及式(2)和式(4),式(1)與式(3)是定義了六面體中心一點,而式(2)與式(4)則需要定義合并后整個六面體的邊界。

為了實現對源圖像的特征值(灰度值)合并,予以下述2個合并性質進行說明。

性質1 對于像素合并區域(i1,j1,i2,j2),若i1=1(i2=N),則該像素合并區域在i負(正)方向上不可合并,對應于六面體合并區域(x1,y1,z1,x2,y2,z2)在x負(正)方向上不可合并;若j1=N(j2=1),則該像素合并區域在j正(負)方向上不可合并,對應于六面體合并區域(x1,y1,z1,x2,y2,z2)在z負(正)方向上不可合并。

性質2 設要合并區域的特征值為g(1),則對于像素合并區域(i1,j1,i2,j2),i1>1,i21。如果存在j2≤j≤j1,i=i2+1(i=i1-1),使得g(i,j)≠g(1),則像素合并區域(i1,j1,i2,j2)在i正(負)方向上不可合并,對應于六面體合并區域(x1,y1,z1,x2,y2,z2)在x正(負)方向上不可合并;如果存在i1≤i≤i2,j=j1+1(j=j2-1),使得g(i,j)≠g(1),則像素合并區域(i1,j1,i2,j2)在j正(負)方向上不可合并,對應于六面體合并區域(x1,y1,z1,x2,y2,z2)在z負(正)方向上不可合并。

2.2 基于深度優先合并的源定義方法

在以上定義和性質的基礎上,深度優先的源定義方法是在像素區域合并基礎上展開的,并轉換成合并的六面體區域參數,從而生成MCNP的輸入卡片。深度優先像素區域合并是每次選擇可合并像素點行(或列)最大的方向為合并方向,流程見圖3。深度優先合并算法的具體步驟如下:

(1)從像素合并區域(1,1)開始找到第一個特征值(灰度值)為g值(大于0)的像素點,以該像素點為初始的像素合并區域;

(2)根據性質1、性質2判斷當前像素合并區域4個方向的可合并性,置相應的標志(真或假);

(3)判斷每個可合并性為真方向上的可合并層數,合并最大層數的方向,特別當i,j方向最大可合并層數相同時,優先合并i方向,上轉步驟(2);

(4)當4個方向可合并性均為假時,本次合并結束。由式(3)與式(4)將像素合并區域(i1,j1,i2,j2)變換為六面體合并區域(x1,y1,z1,x2,y2,z2),并在MCNP輸入文件中寫入成MCNP六面體定義語句;

(5)按順序搜索下一個g值像素點,成功轉步驟(2),否則算法結束。

上述算法完成了整個源圖像上特征值為g值的像素點合并,其他每個特征值依照上述算法合并一輪,就可以完成整個源圖像上所有大于零灰度值的像素點合并定義,從而實現任意形狀灰度源圖像定義。為了減少計算時間,模擬模型選用基于FIR點通量模型。

圖3 深度優先合并算法流程圖

在線度為1.5 mm的E字母源(見圖4a)中,有信息的像素點(即屬于發射中子區域的)數為89,所以在進行合并前需要定義89個六面體并組合為一個柵元作為源區域,而合并后只需要4個六面體,如圖4b所示的合并區域,這樣源定義所需的六面體數大大下降。同時,合并前模擬結果(見圖4c)與合并后模擬結果(見圖4d)一致。

(a)源圖像 (b)合并示意

(c)合并前模擬結果 (d)合并后模擬結果 圖4 線度為1.5 mm時E字母源圖像、合并示意與模擬結果

E字母源模擬的參數和結果見表1。采用深度優先合并后,計算結果最小誤差沒有變化且滿足統計要求,計算時間則減少為不到合并前直接模擬時間的1/5,這說明深度優先合并方法能夠大大降低模擬計算時間。

表1 線度為1.5 mm的E字母源圖像合并前后模擬參數和結果

3 復雜二值源的模擬圖像與重建結果

以“西安交通大學校徽”為源,通過深度優先合并定義,模擬得到其編碼圖像,見圖5。校徽源模擬的相關參數和結果見表2。從計算結果的最小誤差可以看到,模擬結果的漲落較小,也可從圖5b圖像很平滑得到驗證。此外,從計算時間看,盡管深度優先合并的效果使得源柵元定義所需的六面體數降低到原始的1/12,但與圖4所示簡單二值源相比,模擬復雜形狀二值源仍需大量的計算時間。

(a)源圖像 (b)編碼圖像圖5 線度為10 mm的“西安交通大學?；铡蹦M結果

表2 線度為10 mm的“西安交通大學?；铡痹磮D像模擬參數及結果

(5)

從圖6的重建圖像與表3重建結果的相關系數可以看出,RL算法與遺傳算法的重建結果均與源圖像吻合較好,從而驗證了基于深度優先合并的源定義方法的正確性。因為只有在模擬正確的基礎上,才能對校徽這樣的復雜形狀完成重建。

(a)源圖像 (b)RL算法

(c)維納濾波 (d)遺傳算法圖6 線度10 mm的“西安交通大學校徽”源圖像及其重建圖像

重建方法RL算法維納濾波遺傳算法rCC083580836206774

4 灰度源的模擬圖像

(a)16階灰度 (b)16階灰度模擬結果

(c)64階灰度 (d)64階灰度模擬結果

(e)256階灰度 (f)64階灰度模擬結果圖7 3種灰度在線度為7.15 mm時E字母源圖像及其模擬結果

不同灰度階的E字母源圖像與對應模擬獲得的編碼圖像見圖7,其中E字母源的灰度值按中心單點的高斯函數分布設計了16、64和256的3種灰度階數。從圖7中顯示,隨著源圖像灰度階數從16增加為256,編碼圖像中心強度降低,圖像灰度過渡平滑,不同區域間差異不顯著。這說明源圖像灰度階由低變高,對應源區上中子發射強度變化梯度由大到小,進而使得編碼結果的強度梯度也由大到小變化,因此高灰度階模擬結果與平滑過渡的實驗結果更為相符。當灰度階數從64增大到256時,模擬結果變化并不明顯,這表明E字母源圖像按照64階灰度定義能滿足編碼圖像平滑過渡的要求。相關的模擬參數和結果見表4。

從表4可得,隨著源圖像灰度階數的不斷增多,雖然合并前源柵元的定義所需的六面體數目變化不大,但是深度優先合并后的六面體數在3種灰度下差異很大。這主要是在較高灰度階數下,像素點灰度值細微的變化已經導致這些像素點很難合并到一個六面體定義中。隨著源柵元中定義的六面體數的增多,更高灰度階的計算時間不斷增加,不過計算結果的最小誤差始終控制在1%以下,符合重建算法的重建對象需求。

表4 3種灰度在線度為7.15 mm時E字母源圖像的模擬參數和結果

5 結論和展望

本文建立了基于深度優先合并方法,采用FIR記錄點通量的模擬模型,實現了復雜二值源圖像和灰度源圖像的定義和模擬。在二值E字母源下,計算結果最小誤差沒有變化且滿足統計要求,計算時間減少為不到合并前直接模擬時間的1/5。對“西安交通大學校徽”進行了復雜二值源的模擬,3種重建方法的結果表明,基于深度優先合并的源定義方法是正確的,對16、64和256的3種灰度階數的E字母源圖像的模擬計算的最小誤差小于1%。這樣,基于深度優先合并方法的復雜源圖像定義方法可進行復雜源圖像(二值圖像與灰度圖像)的定義,進而能對中子及其他粒子編碼成像實現高效模擬計算。

此外,聚變反應區的燃燒過程是一個時空耦合的局部與整體均在演化的過程。該演化過程在時間上歷經點火—燃燒—最終沉寂等階段,同時在空間上對應歷經等離子體區域尺寸由小到大再到小的變化。在此期間,中子發射的空間強度分布變化十分劇烈,表現為源圖像的形狀和灰度分布變化十分復雜。可見,采用深度優先合并方法定義的復雜形狀源正好適用于對這樣復雜變化過程的編碼成像實施正過程模擬,進而實現對時空耦合源區各時段圖像的重建,所得各時段源圖像的局部與整體的邊緣可望更為清晰,此有利于正確認識聚變反應區時空耦合演化過程。

[1] PFALZNER S. 慣性約束聚變導論 [M]. 崔旭東, 譯. 北京: 原子能出版社, 2011: 1-20.

[2] VOLEGOV P L, DANLY C R, FITTINGHOFF D N, et al. Self characterization of a coded aperture array for neutron source imaging [J]. Rev Sci Instrum, 2014, 85(12): 123506.

[3] VOLEGOV P, DANLY C R, FITTINGHOFF D N, et al. Neutron source reconstruction from pinhole imaging at National Ignition Facility [J]. Rev Sci Instrum, 2014, 85(2): 023508.

[4] WILSON D C, ARAGONEZ R J, ARCHULETA T N, et al. Comparing neutron and X-ray images from NIF implosions [J]. EPJ Web of Conferences, 2013, 59: 04002.

[5] GRIM G P, ARCHULETA T N , ARAGONEZ R J, et al. Summary of the first neutron image data collected at the National Ignition Facility [J]. EPJ Web of Conferences, 2013, 59: 13017.

[6] HURRICANE O A, CALLAHAN D A, CASEY D T, et al. Fuel gain exceeding unity in an inertially confined fusion implosion [J]. Nature, 2014, 506(7488): 343-348.

[7] 張天奎. 快中子編碼成像技術中的圖像重建方法研究 [D]. 西安: 西安交通大學, 2013.

[8] BRIESMEISTER J F. MCNP: a general Monte Carlo N-particle transport code, Version 4C [R]. Los Alamos, USA: Los Alamos National Laboratory, 2000.

[9] SHULTIS J K, FAW R E. An MCNP primer [M]. Manhattan, USA: Kansas State University, 2006: 16-20.

[10]吳岳雷. 中子半影成像系統關鍵部件研究及研制 [D]. 西安: 西安交通大學, 2010.

[11]鄭向陽. 中子針孔半影成像系統設計方法與針孔試制工藝研究 [D]. 西安: 西安交通大學, 2007.

[12]李成. 慣性約束聚變中子編碼孔試制工藝研究 [D]. 西安: 西安交通大學, 2009.

[13]秦娟. 用于瞬態輻射過程診斷的閃爍光纖成像陣列樣品試制 [D]. 西安: 西安交通大學, 2008.

[14]西安交通大學核科學與技術系. MCNP-3B/PC程序使用說明書 [R]. 西安: 西安交通大學核科學與技術系, 2004.

[15]BARRERA C. Experimental component characterization, Monte Carlo-based image generation and source reconstruction for the Neutron Imaging System of the National Ignition Facility [D]. Berkeley, USA: University of California at Berkeley, 2007.

[16]X-5 Monte Carlo Team. MCNP: a general Monte Carlo N-particle transport code, version 5 [R]. Los Alamos, USA: Los Alamos National Laboratory, 2003.

[17]陳法新, 鄭堅, 李正宏. 用MCNP程序計算針孔成像 [J]. 計算機應用與軟件, 2010, 27(7): 185-186. CHEN Faxin, ZHENG Jian, LI Zhenghong. Pinhole imaging calculated by MCNP code [J]. Computer Applications and Software, 2010, 27(7): 185-186.

[18]戰元平. 用于快中子成像的閃爍光纖陣列質子串擾抑制方法研究 [D]. 西安: 西安交通大學, 2010.

[19]施燦輝. 基于圖像的MCNP數字人體建模與仿真研究 [D]. 合肥: 合肥工業大學, 2003.

[20]RICHARDSON W H. Bayesian-based iterative method of image restoration [J]. Journal of the Optical Society of America, 1972, 62(1): 55-59.

[21]LUCY L B. Iterative technique for rectification of observed distributions [J]. The Astronomical Journal, 1974, 79(6): 745-754.

[22]ZHAO Z G, DING Y K, DONG J J, et al. Richardson-Lucy method for decoding X-ray ring code image [J]. Plasma Physics and Controlled Fusion, 2007, 49(8): 1145-1150.

[23]RESS D, LERCHE R A, ELLIS R J, et al. Neutron imaging of laser fusion-targets [J]. Science, 1988, 241(4868): 956-958.

[24]RESS D, LERCHE R A, ELLIS R J, et al. Neutron imaging of inertial confinement fusion-targets at Nova [J]. Review of Scientific Instruments, 1988, 59(8): 1694-1696.

[25]ZHANG T K, HU H S, JIA Q G, et al. Genetic algorithms applied to reconstructing coded imaging of neutrons and analysis of residual watermark [J]. Rev Sci Instrum, 2012, 83(11): 113505.

[26]金穎. 慣性約束受控熱核聚變中的中子成像技術研究 [D]. 大連: 大連理工大學, 2000: 45-63.

(編輯 苗凌)

Highly Effective Simulation Strategy for Coded Imaging of the Complex Fast Neutron Source

ZHANG Fengna1,ZHANG Tiankui2,WANG Dongming1,ZHENG Xiangyang3,HU Huasi1

(1. School of Energy and Power Engineering, Xi’an Jiaotong University, Xi’an 710049, China; 2. Laser Fusion Research Center, CAEP, Mianyang, Sichuan 621900, China; 3. Nuclear and Radiation Safety Center, MEP, Beijing 100082, China)

To solve the difficulties in complex gray-scale source definition and simulation efficiency, a strategy for defining complex gray-scale source based on the depth-first merger is proposed. A technology of cells merging is used to reduce the number of cells for the source definition and ensure the sampling efficiency of source in Monte Carlo simulation, thus the simulation efficiency of coded imaging of the complex gray-scale source can be increased. When the simulated source is the binary letter E source, the error of simulated coded image of the merged source remains same as that of the source without merging and meets the requirement of the reconstruction for statistics, while the computing time for merged source is reduced to 1/5 of the source without merging. The 16, 64 and 256 gray-scale letter E sources are defined by the definition method of complex gray-scale source based on the depth-first merger, and the corresponding simulations of coded image are also carried out. The errors of simulation results reach less than 1%, which meet the requirement of the reconstruction for statistics. The complex binary source of Xi’an Jiaotong University school badge is defined with the proposed strategy for complex gray-scale source based on the depth-first merger and the corresponding simulation of coded image is carried out, then the reconstructions of the source by three different methods are realized. The strategy for defining complex gray-scale source based on the depth-first merger is verified, so it can be expected to provide a feasible technological approach to enhance the adaptability of diagnosis of the complex capsule implosion process.

fast neutron coded imaging; depth-first merger; complex gray-scale source; Monte Carlo simulation; image reconstruction

2016-02-20。 作者簡介:張鳳娜(1986—),女,博士生;胡華四(通信作者),男,教授。 基金項目:教育部高校創新研究團隊資助項目(IRT1280);陜西省自然科學基礎研究計劃資助項目(2015JZ001);國家自然科學基金資助項目(11505166)。

時間:2016-05-10

10.7652/xjtuxb201607018

TL65+7

A

0253-987X(2016)07-118-07

網絡出版地址:http:∥www.cnki.net/kcms/detail/61.1069.T.20160510.1523.014.html