中子半影成像的软件系统.doc

中子半影成像的软件系统赵宗清1,2丁永坤1郝轶聃1 袁永腾1 李志超1(1.中国工程物理研究院激光聚变研究中心，四川绵阳，621900；2.清华大学工程物理系，北京 100084)摘要：中子半影成像的软件系统包括成像过程模拟和图像重建算法研究。目前，我们已分别发展了Monte-Carlo方法和解析模拟方法，对成像的物理过程进行了细致的模拟，两种方法得到了非常一致的比对结果。我们发展了维纳滤波法和Richardson-Lucy方法用于图像重建，并分别采用两种重建算法对模拟成像的结果进行了解码。通过数值实验发现，Richardson-Lucy方法得到的结果与源符合得更好，较维纳滤波法更具可行性。据此详细考察了散射中子、编码孔加工误差及整个成像系统的瞄准精度对编码图像的影响，为成像系统硬件部分的设计提供了重要的参考意见。关键词：蒙特卡罗方法；解析模型；维纳滤波；Richardson-Lucy方法引言：中子成像技术，可以获得内爆中子产生区域图像，从而验证内爆动力学计算的正确与否，了解靶丸压缩与产额之间的关系，获得各种黑腔、靶丸设计的性能参数，达到优化靶及黑腔设计等目的。美国劳伦斯里弗莫尔实验室（LLNL）的D.Ress，R.A.Lerche等人，在实验上采用针孔成像技术，首次获得了内爆中子图像，证实了中子成像技术的可行性，该结果发表在1988年的Science上 1。此后，美国，法国和日本的科学家2-8在对中子成像的模拟，解谱和实验等方面进行了深入的研究。其中，2002年，L.Disdier 3,4等人在Omega装置上利用中子半影成像技术实现了20m的分辨率。而在即将建成的美国NIF装置和法国LMJ装置上要实现510m的分辨率。在神光II激光装置上，中子产额较低，不能满足中子成像的要求 10。但随着我国ICF事业的发展，在神光III原型装置上DT中子产额有望达到1010/发以上，这为中子成像技术的发展提供了有利的条件。因此，为了提高我国ICF核诊断的水平，我们将大力发展中子半影成像技术。中子半影成像系统的建设，包括硬件平台建设和软件平台建设两部分。硬件平台的建设包括编码孔的设计加工，图像记录系统研制和瞄准技术研究等。目前我们在编码孔加工技术和瞄准技术方面还有一定的技术瓶颈，需要进行突破。同时，为了考察中子成像技术的可行性，并且为硬件设计提供定量化的参考意见，软件平台的建设必须先于硬件平台建设。图1 中子半影成像硬件系统2成像过程的数值模拟2.1 蒙特卡罗模拟11Monte-Carlo方法是一种随机抽样方法。利用已知的中子反应截面数据, 模拟各种微观物理过程, 通过概率抽样对源粒子的行为进行跟踪,决定每次碰撞后次级粒子的运动方向和速率。我们采用MCNP作为计算程序。MCNP是一个大型多功能的蒙特卡罗计算程序, 可处理复杂场所三维几何结构的中子光子耦合输运问题。MCNP具有较强的通用性, 在源描述、空间物质的几何分布上具有很大的灵活性,可处理任意三维几何结构问题, 适用面宽。对中子的输运问题,MCNP详细处理了各种微观物理过程,适用于中子成像的要求。利用MCNP模拟中子成像，能使我们更好地理解中子在成像系统中的输运过程，同时能够得到图像重建所需要的点扩散函数，为图像重建算法的研究提供条件。2.2 解析模拟12通过物理分析和几何建模，我们得到成像传递函数，由两部分的乘积组成，即几何衰减因子(1)及吸收衰减因子(2)其中为中子束在锥孔材料中的吸收路程，为材料的中子吸收系数。详细内容可参看文献12。1.3 两种模拟方法的比较图2 蒙特卡罗方法和解析方法得到的点扩散函数图2是蒙特卡罗模拟和解析模拟分别得到的点扩散函数。从图2中可以看出，蒙特卡罗模拟得到的图像考虑了散射等物理因素的影响，存在着统计误差，而解析法得到的图像非常光滑。Monte-Carlo方法的优势在于通过对单中子事例的微观物理过程的计算来实现宏观模拟，因此其结果在包括噪声分布等的各方面均非常接近真实物理过程。而解析方法的优势在于，忽略噪声、对中子束的几何光学近似及采用各种宏观物理现象的唯象理论，使得计算效率极高。因此，我们的发展方向定位为，采用解析方法对中子成像系统进行定性或半定量估计，并在此基础上采用Monte-Carlo方法做进一步的细致研究，并且采用并行机提高Monte-Carlo的计算效率。3 图像重建3.1 成像模型中子半影成像，主要有两个步骤5：第一步：由中子源发射中子，经过半影孔形成编码像，编码像由亮区、半影区和本底组成。源的亮度分布o(x,y)通过半影孔点扩散函数h(x,y)之后，像的亮度分布i(x,y)为：（3） 其中，*代表二维卷积，n (x,y)为干扰噪声。 第二步：由编码像经过图像反演得到源的亮度分布，即中子通量分布。图像反演的基本思想是：先分别对i(x,y)，h(x,y)和n(x,y)做傅立叶变换得到其频域分布I(kx,ky)，H(kx,ky)和N(kx,ky)，相应地： （4）通过一定的图像反演算法，就可以得到o(kx,ky)，再进行逆傅立叶变换便可以得到源的中子强度的二维空间分布。3.2 维纳滤波维纳滤波就是寻找传输函数为Hw(kx,ky)的滤波器，使得恢复出的目标估计值Ow(x,y)与真实的O(x,y)的均方差达到最小，即： （5）经过推导，可得 （6）其中，为编码孔成像系统的传递函数的复共轭，和分别是靶图像和噪声的功率谱密度。3.3 Richardson-Lucy方法14,15如果成像方程服从泊松分布： （7）O是真实的源分布，是PSF（从物体真实位置j发出的光散射到观察像素i上的份额）。I是无噪声污染的退化图像。观测计数和期望计数的联合概率L为： （8）经过一系列的推导，可以得到R-L迭代公式为： （9）从方程（9）可以明显看出，如果Richardson-Lucy迭代收敛16，意味着随着迭代进行，修正因子将趋于单位值。因此，对于泊松分布的数据，它必然收敛于最大似然解。中子成像是把通过CCD记录中子通过编码孔传输以后所成的编码像，图像分布近似服从于泊松分布，因此，Richardson-Lucy方法是可行的。3.4 解谱结果比较图3 解谱结果比较 由图3可见，Richardson-Lucy非线性迭代法对噪声有很好的抑制作用，得到了比较理想的结果，与源符合得较好。而维纳滤波得到的结果（图3d），重建源图像有一定的展宽，周围出现了一些调制。与同一数据的维纳滤波结果相比，Richardson-Lucy方法具有明显的优越性17。RL方法作为一种非线性迭代法，迭代次数过多会导致高频放大而使结果不可信；但由于数学上的困难，目前尚无法给出一个迭代终止的有效判据，尤其是ICF中子源无法定标，使得我们需要进一步寻找其它更加合适的解码方法，如遗传算法等。4 系统参数优化通过数值模拟，我们详细考察了散射中子，锥孔加工精度和瞄准精度对成像质量的影响。下面给出一些定量结果，更为具体的内容可参考文献11-13。1) 加工精度:双锥的加工椭圆度好于0.99。2) 瞄准精度:瞄准时，视场中心误差控制在50微米之内；编码孔旋转角度不小于0.1mrad;物距误差控制在300微米之内。3) 编码孔优化设计:为了有效避免编码孔弹性散射的影响，编码孔外径为1cm左右。同时，为了增加信号的对比度，编码孔材料选择为至少5cm长的钨合金。5 结论我们发展了两种模拟方法和两种图像重建算法，建立了中子半影成像的软件系统。利用软件平台，我们对系统参数对成像质量的影响进行了定量的考察，为硬件设计提供了参考意见，有利于中子成像技术的发展。参考文献1 D.Ress, R.A.Lerche, R.J.Ellis et al,Neutron imaging of laser fusion targets, J. Science, 1988, v241,p956,.2 K.A. Nugent, coded imaging of thermonuclear neutrons, J.Review of Scientific Instruments, 1988,v59,p1658,3 L.Disdier, A.Rouyer, D.C.Wilson et al, high-resolution neutron imaging of laser imploded DT targets, J. Nuclear instruments and Methods in Physics Research A , 2002,v489, p496,4 L.Disdier, A.Rouyer, A.Fedotoff et al, neutron imaging of ICF target plasmas J.Review of Scientific Instruments, 2003,v 74, p1832, 5 G. P. Grim, G. L. Morgan, M. D. Wilke et al, Progress on neutron pinhole imaging for inertial confinement fusion experiments, J.Review of Scientific Instruments, 2004,v75, p3573, 6 L. Disdier, A. Rouyer, I. Lantujoul et al, Inertial confinement fusion neutron images J.Physics Of Plasmas, 2006, v13, p056317, 7 O.Delaga, R.A.Lerche, T.C.Sangster et al, SiRinc: A code for assessing and optimizing the neutron imaging diagnostic capabilities in inertial confinement fusion experiments J. Review of Scientific Instruments, 2001,v72, p869, 8 R.A.Lerche, N.Izumi, R.K.Fisher et al, neutron images recorded with high-resolution bubble detectors J. Review of Scientific Instruments, 2003,v74, p1709, 9 Y.-W, Chen, R.Kodama and Z.Nakao, penumbral imaging of laser-imploded targets and its blind reconstruction, IEEE transactions on nuclear science, 1998,v45,p992,10 Zhao Zongqing, Chen JiaBin, Ding Yongkun et al, Exploration of the yield-ratio method at Shenguang II laser facility J. Rev. Sci. Instrum., 2006, V77, P073502,11 赵宗清，丁永坤，刘东剑等，中子半影成像的数值模拟，J.强激光与粒子束，2006，18(7)，120312 郝轶聃，缪文勇，赵宗清等，中子半影成像中椭圆误差的解析计算，J.强激光与粒子束，2007，19(3)，50713 赵宗清，丁永坤，郝轶聃等，瞄准误差对中子成像的影响，中国科学E辑(投稿)14 William Hadley Richardson, Bayesian-based Iterative method of image restoration J. Journal of optical society of America, 1972,V62,P55,15 L.B.Lucy, An iterative technique for rectification of observed distributions J.The Astronomical journal, 1974 ,V79,P745,16 Shepp, L. A., & Vardi, Maximum Likelihood Reconstruction in Emission TomographyJ. IEEE Trans Med Imag, 1982, MI-1:11312217 Zhao Zongqing, Ding Yongkun, Dong Jianjun et al, Richardson-Lucy method for decoding X-ray ring code image J. Plasma Phys. Control. Fusion, 2007,V49 ,P11451150,Software system of neutron penumbral imagingZhao Zong-qing1,2 Ding Yong-kun1 Hao Yi-dan1 Yuan Yong-teng1 Li Zhi-chao11 Laser Fusion Research Center, China Academy of Engineering Physics, Mianyang, 621900, China2 Department of Engineering Physics, Tsinghua University, Beijing, 100084, ChinaAbstract:Software system of neutron penumbral imaging includes simulation of imaging process and research of image reconstruction algorithms. Monte-Carlo method and analytic method are used to simulate t