伍华杰 甘佺 吕永强 汪启宏 伍一凡 宋婧 凤麟核团队

1（中国科学院合肥物质科学研究院 合肥 230031）

2（中国科学技术大学 合肥 230026）

3（中科超睿技术有限公司 青岛 266041）

4（中子科学国际研究院 青岛 266199）

5（山东省中子科学技术重点实验室 青岛 266199）

中子照相技术作为一种对物体内部结构进行无损检测的手段，因其具有穿透能力强［1］、轻元素灵敏、成分识别精准［2］及抗干扰强等独特的优势，在航天航空［3-4］、锂电池［5-6］、混凝土建筑和非开箱安检［7］等方面检测中都具有不可替代的重要价值。随着中子照相技术的发展，包括基于反应堆、大型加速器、小型加速器和微型中子源等的中子照相系统陆续投入研发和市场，因此中子照相无损检测具有良好的发展前景。虽然中子照相系统逐渐向微小型化发展，但由于基本的成像结构未变，中子照相图像依然存在几何不锐度，亟须开展关于几何不锐度修正的方法研究。

中子照相图像几何不锐度的修正，属于中子照相图像的复原问题，且一直是一个难点。金炜等［8］提出一种基于小波降噪的Richardson-Lucy（wavelet-RL）算法，假设中子照相图像中的少量噪声为高斯白噪声，采用BivaShrink双变量去噪模型，刻画小波变换后图像相同位置处不同尺度子带系数间的相关性，显著提升图像的去噪性能，从而提高图像几何不锐度的修正效果。Souza 等［9-10］提出摇摆曲线（Rocking Curve，RC）的概念，RC 半峰全宽（Full Width at Half Maximum，FWHM）表示几何不锐度，并采用Richardson-Lucy（RL）算法进行复原，图像的前景边缘部分变得更加锐利，但是未考虑γ白斑噪声（Gamma Spots Noise，GSN）的去除。乔双等［11］将调制核引入到正则化RL方法中，提出一种用于中子照相图像去模糊的方法SK-TV-RL，能够抑制中子照相图像弱噪声的放大，并能复原图像中的细节信息。Yazid 等［12］对RTP（PUSPATI TRIGA Mark II）反应堆得到的中子照相图像，采用中值滤波与RL结合的方法进行复原，然而直接采用中值滤波去除GSN会损坏图像前景的全局细节与边缘信息，降低图像的整体复原效果。

目前，有关几何不锐度修正的研究大多集中于含有弱噪声的中子照相图像，对于含有强噪声的中子照相图像的研究较少，而在采用RL方法修正几何不锐度的过程中，中子照相图像中的GSN难以去除并会不断放大，严重影响图像的质量。因此，本文提出一种基于改进RL 的中子照相图像几何不锐度修正算法，并与以往有关中子照相图像的几何不锐度修正算法作对比，验证采用改进的RL算法进行几何不锐度修正的效果。

本文采取的方法主要分为两个过程：第一步是根据中子照相图像的几何不锐度，对用于表示几何不锐度的点扩散函数（Point Spread Function，PSF）建立数学模型；第二步是结合第一步得出的PSF，采用改进的RL 方法对图像中的前景边缘细节进行复原，以达到更好的几何不锐度修正效果，提升图像的分辨率。

1 中子照相图像复原方法

实际上，中子照相系统只能获取到降质后的中子照相图像，可以大致将中子照相图像的降质过程描述为一个离散线性系统［13］：

式中：Idegraded表示降质中子照相图像；I表示原始清晰中子照相图像；k表示PSF；n表示加性噪声；*为卷积。中子照相图像的复原需要同时考虑PSF和加性噪声，既要准确估计PSF，也应尽可能地去除加性噪声的干扰。

本文提出的方法中，首先根据中子源出口到像探测器的距离、中子源出口直径和被检品中心到像探测器的距离，精准地估计出中子照相图像的PSF，然后根据中子照相图像中GSN的特点，提出一种改进的RL（improved-RL）算法。算法总体流程如图1所示。

图1 算法总体流程图Fig.1 Block diagram of the improved-RL algorithm

1.1 中子照相图像的PSF估计

一般地，中子照相图像的PSF 估计分为三个过程：第一步是计算中子照相图像的几何不锐度；第二步是为PSF 选择合适的分布模型；第三步是根据几何不锐度和分布模型估计出PSF。

通常情况下，在中子照相系统中，通过中子源出口的中子束穿过被检品某点后，在像探测器处会产生一圆形光斑，如图2 所示。其中，E表示被检品边缘上的点，C表示被检品的中心点。该圆形光斑被称为几何不锐度，其直径大小一般为［14］：

图2 中子照相系统结构示意图Fig.2 Structural diagram of the neutron imaging system

式中：L表示中子源出口到像探测器距离；D表示中子源出口直径；d为被检品中心到像探测器的距离。

与X 射线衍射场相似，中子源各点处发射的中子束强度分布形状为RC［15］。同理，所有穿过被检品边缘某一点后的中子束强度分布形状也为RC。将RC 的FWHM 表示为中子照相图像的几何不锐度G［8］。如图3所示，其中，W表示RC的FWHM。用于表示PSF 的函数模型包括高斯分布［16］、柯西分布［12，17］等。本文采用二维柯西分布模型来表示PSF，表达式为：

图3 摇摆曲线示意图Fig.3 Sketch of the rocking curve

式中：Cauchy(x，y)表示二维柯西分布函数；ξ表示为分散度参数。几何不锐度G与参数ξ的关系为：

G确定后，根据式（4）可得参数ξ的值，代入式（3）中可得表示中子照相图像的PSF。

1.2 改进的RL算法

RL 算法是Richardson 和Lucy 提出的一种非盲反卷积图像复原方法，算法表达式如式（5）所示：

式中：It表示第t次迭代的结果；k表示PSF；kT表示k在水平方向的翻转。

由于照相得到的中子照相图像中存在灰度值较高的GSN，若直接采用RL算法进行复原，得到的复原图像会出现严重的振铃效应。因此，建立如下数学模型：

式中：γ(·)为被GSN 污染的区域，γ͂(·)为未被GSN 污染的区域，在该数学模型下需要检测出所有的γ(·)。由于GSN具有孤立性，可采用拉普拉斯算子检测出γ(·)，并配合中值滤波去除γ(·)内的GSN，具体方法如式（7）所示：

式中：M(·)表示中值滤波；U(x，y)表示点(x，y)的邻域；∇2Idegraded(x，y)表示Idegraded二阶梯度之和；T为指定阈值。

上述过程只会对含有GSN 的区域进行去噪处理，即使在多次使用的情况下，仍然不会损坏图像中的高频信息。因此，将以上步骤引入到RL 算法中，以对图像复原过程中的振铃效应进行抑制。改进的RL算法的表达式为：

改进的RL 算法在中子照相图像反卷积前对GSN进行去除，抑制因GSN产生的振铃效应。

2 实验结果和分析

本文中子照相系统的中子源类型为热中子（能量区间为0.005～0.5 eV），成像端采用LiF 转换屏和电荷耦合元件（Charge Coupled Device，CCD）结合的数字成像方式。其中，该CCD 相机为Atik Infinity Camera，分辨率为1 392×1 040。为降低CCD相机的暗电流，将CCD相机放置于密封的真空箱内，温度降至-10 ℃。另外，在真空箱外包裹一层铅屏蔽层以防止中子束或次级粒子的辐射。本文选取线对样品进行实验［18-19］，线对样品实物如图4 所示。线对样品拥有多组粗细不同的线状区域，有利于检测中子照相系统的分辨率。

图4 线对样品实物图Fig.4 Photograph of the line-pair sample

首先，测量中子照相系统的参数值，得到中子源出口到像探测器距离L与中子源出口直径D的比值，即L/D=80，被检品中心到像探测器的距离d=2 cm，转换屏大小为20 cm×20 cm。拍照方式设置为单次曝光，曝光300 s后记录并保存图像。由于保存后的图像尺寸较大，且线对样品的特征信息主要集中在图像的中心区域，需要进行裁剪，裁剪后的图像尺寸为512 pixel×512 pixel；然后根据式（2）计算出中子照相图像的几何不锐度，再根据式（3）、（4）计算出对应的PSF；最后对线对中子照相图像分别采用RL 算法、小波降噪RL算法、调制核正则化RL算法、文献［12］方法和本文改进的RL算法进行图像复原，并绘制复原后图像的行灰度曲线，实验结果如图5 所示。从实验结果可看出，相对于RL 算法、小波降噪RL 算法和调制核正则化RL 算法，经本文改进的RL算法处理后，复原后图像噪声的波动幅度明显减小；相对于文献［12］方法，复原后图像行灰度曲线的峰谷值更大、波峰更窄，几何不锐度修正程度有较大的提升。

图5 线对样品的中子照相图像、复原结果与对应的行灰度曲线图(a) 原始中子照相图像，(b) RL算法，(c) 小波降噪RL，(d) 调制核正则化RL，(e) 文献[12]方法，(f) 本文方法，(g) 原始中子照相图像对应行灰度曲线图，(h) RL算法对应行灰度曲线图，(i) 小波降噪RL对应行灰度曲线图，(j) 调制核正则化RL对应行灰度曲线图，(k) 文献[12]方法对应行灰度曲线图，(l) 本文方法对应行灰度曲线图Fig.5 Neutron photographic image of line-pair sample, restored results and corresponding line grey profile(a) Original neutron photographic image, (b) RL algorithm, (c) Wavelet-RL, (d) SK-TV-RL, (e) Method in Ref.[12], (f) Our method,(g) Corresponding line grey profile of original neutron photographic image, (h) Corresponding line grey profile of RL algorithm, (i)Corresponding line grey profile of wavelet-RL, (j) Corresponding line grey profile of SK-TV-RL, (k) Corresponding line grey profile of method in Ref.[12], (l) Corresponding line grey profile of our method

为客观衡量几何不锐度修正效果，通过无参考图像质量评价指标平均梯度（Average Gradient，AG）和空间频率（Spatial Frequency，SF）对中子照相图像复原质量进行分析，如表1 所示。其中，AG 可以反映图像中前景信息的纹理变化特征，SF能够反映图像的整体活跃程度，均可衡量图像的清晰程度。

表1 文献[12]方法和改进的RL算法复原图像后的平均梯度和空间频率Table 1 AGs and SFs of restored images based on the algorithm presented in Ref.[12] and improved-RL algorithm

AG的定义如下：

式中：M和N分别表示图像的宽度和高度；D为梯度算子。

SF的表达式如下：

一般来说，AG和SF值越大，表明图像的前景纹理层次越多、信息结构越复杂。图5中可看出，原始中子照相图像、经RL算法、小波降噪RL算法和调制核正则化RL 算法复原后的中子照相图像仍然存在大量的GSN，而当图像中存在较多的GSN 时，指标会异常增大，影响AG 和SF 对图像前景信息复原效果的判定，因此，只计算文献［12］方法和本文方法复原后图像的AG和SF指标，如表1所示。

从表1可以看出，与文献［12］方法相比，经过本文提出的算法复原后，线对中子照相图像AG 和SF值分别提高了60.23%和29.90%，图像锐度有明显的提升。

因此，通过结合图5 和表1 的结果可看出，本文中改进的RL 算法相比于RL 算法、小波降噪RL 算法和调制核正则化RL 算法，能抑制GSN 在中子照相图像复原过程中产生的振铃效应，更有利于图像前景细节信息的复原；相比于文献［12］方法，有更好的几何不锐度修正效果。

3 结语

本文提出了一种基于改进RL 的中子照相图像几何不锐度修正算法对几何不锐度进行修正。通过准确计算表示几何不锐度的PSF，结合改进的RL算法，修正了因几何不锐度引起的图像降质，同时避免了GSN对图像复原过程的影响。实验结果表明，本文提出的方法可以获得高质量的复原图像，能够进行高质量的中子照相无损检测，对提升中子照相系统的检测效率与精度意义重大。

作者贡献声明伍华杰负责图像数据的获取、算法的设计、程序的编写、文章的编写与修订；甘佺负责实验装置的搭建、文章的编写与修订；吕永强提供部分算法思路、文章的修订；汪启宏负责提供部分算法思路；伍一凡负责验证算法的正确性；宋婧负责提供实验装置的指导。