针对微尺寸X 射线源的非相干全息层析成像*

陈纪辉 王峰 理玉龙 张兴 姚科 关赞洋 刘祥明

1) (中国工程物理研究院,激光聚变研究中心,绵阳 621900)

2) (复旦大学,核物理与离子束应用教育部重点实验室,上海 200433)

1 引言

惯性约束聚变(inertial confinement fusion,ICF)是实现可控核聚变反应的一条重要途径[1].对ICF 中高温等离子体的X 射线辐射的空间分辨研究,可以提供ICF 靶内由于做功和能量(电子热传导、辐射热传导)输运导致的流体状态的空间演化信息[2,3].目前,用于X 射线成像的设备主要有:X 射线针孔相机[4]、KB 反射镜[5,6]、Wolter 显微镜[7]以及球面弯晶[8].上述成像诊断仪器采用的都是直接成像方式,即通过光透射、反射、衍射等传输过程,在探测平面直接反映成像目标的光辐射空间分布的高分辨二维图像.然而,通过这种方式采集到的二维图像丢失了沿探测方向上的强度分布信息,其包含的光强分布图像实际上是在一定轴距范围内光强的沿轴积分叠加.综上所述,现有的诊断手段还无法通过单次测量实现对射线源的层析成像.

全息技术是实现对物体三维成像的一条有效途径.相干全息术利用相干性极强的光源(通常是激光)作为参考光,在照明物体后,将物光波前的振幅和相位以干涉条纹的形式记录为全息图.之后再通过适当的方式,从这幅包含物光“全部信息”的全息图中重建出物体的三维形貌[9].传统的相干全息术对于照明光源的相干性要求极高,而非相干光源在自然界中更为普遍,这在很大程度上限制了其应用范围.1961 年,Mertz 和Young[10]将全息技术推广至非相干光领域.他们将针孔相机中的透光小孔替换为菲涅耳波带片,使得不同位置的光透过波带片后在底片上产生放大倍数不同、中心位置不同的波带片形貌投影.这些叠加投影形成的全息图反映了原光源的空间信息,通过适当的重建方法就能还原出其在各个位置的光强分布.在非相干全息理论被提出后不久,核医学[11,12]、ICF 成像诊断[13,14]以及天文学[15]等领域都采用了这种具有高收光效率的“编码成像”方式.近年来,非相干全息技术被广泛应用于无透镜成像领域,并衍生出了菲涅耳相机的理念[16,17].2020 年,Wu 等[18]在反演菲涅耳相机底片(非相干全息图)的过程中引入了信号传输中压缩感知优化算法,获得了高分辨率、低背景噪声的二维重建图像.2023 年,Soltau 等[19]使用菲涅耳波带片作为编码孔径完成了高光通量的X 射线全场荧光成像.值得注意的是,到目前为止的非相干全息成像(编码成像)研究大都集中于二维光源图像的采集与重建,而早期利用Gabor 波带片进行三维编码成像的尝试,受到了非聚焦像的严重干扰[20,21].

本文在前人工作的基础上,介绍了一种基于非相干全息成像机制与压缩感知重建算法,针对微尺寸X 射线源的层析成像技术;并通过模拟三维光源的照明全息记录和层析重建阐明了该诊断手段的应用价值.模拟结果显示了在选用的实验参数体系下,非相干全息成像系统能对轴向尺寸为16 mm的光源进行层析成像,以4 mm 为步长的重建图像中几乎不存在来自其他层面信号的干扰.

2 非相干全息层析成像机制

2.1 非相干全息成像理论模型

非相干全息的成像过程如图1 所示,空间体光源可以看作是数个点光源的集合,位于不同空间位置的点光源独立透过波带片,并在探测器上留下不同中心位置和放大倍数的波带片投影(当像距远小于波带片焦距时,不考虑衍射带来的影响),各个投影叠加形成的非相干全息图由下式表示:

图1 非相干全息成像示意图Fig.1.Procedure of incoherent holography.

其中,h是探测器上获取的光强分布;(x′,y′,z′)表示物方区域的坐标;t'是波带片投影在探测面上的强度分布;o是待测体光源的光强分布,M是成像缩放倍数,且该光源延轴向的固有长度为z0.这里考虑将正弦波带片作为“编码孔径”,其在探测面上的投影可以写为

其中,r*是经过投影放大后的波带片最内环半径:

式中,d是波带片距离探测器之间的距离(像距),z是光源与波带片之间的距离(物距),r1是波带片孔径的最内环半径.根据(1)式的表述,对于二维光源,非相干全息图可以视为光源平面上的光强分布与该物距下产生的波带片投影的卷积;而三维光源产生的全息图则是通过所有物距平面产生的二维全息图累加形成.通常,在空间域较为复杂的卷积和退卷积运算可以转化为在频域中更为简洁的表述,即

式中,H是h的二维傅里叶变换,而T和O分别是t和o在给定物距上的二维傅里叶变换.使用传输算子A表征物空间光强分布与非相干全息图之间的映射:

其中,F-1代表二维傅里叶逆变换算符,Ti代表(4)式中T沿Z方向的第i个分量.Fi意为取三维矩阵中沿Z方向第i个分量进行二维傅里叶变换.因此,h与o之间的运算关系可以表示为

2.2 压缩感知重建模型

上述三维非相干全息成像过程可以看作是一个对数据的压缩采样过程,即使用具有(Nx×Ny)个探测单元的探测器去采集体量为(Nx×Ny×Nz)的目标信号.因此,在反解二维全息图重建三维光源光强分布的过程中,不可避免地会遇到方程个数远小于未知数的情况.根据压缩感知理论[22,23],这种求解“欠定方程组”的问题需要引入正则条件来对方程的解空间范围施加约束,继而求出既满足全息映射过程,又符合正则约束的最优解.本文三维图像信号通过最小化目标函数重建:

其中,第1 项为残差项,表征反解出的三维图像经过非相干全息投影后与实际采集到的全息图之间的差异;第2 项为正则项,对于微尺寸X 射线源的三维图像,使用它的全变分范数(total variation norm,TV norm)作为正则项来表征图像信号的稀疏性;β为正则参数,用于平衡目标函数中残差项和正则项的比重.采用两步迭代收缩阈值算法[24](TWIST)对(7)式进行求解,以获得不同物距层面上的光强分布.

3 模拟成像与重建结果

为了验证非相干全息层析成像技术的可行性,模拟了针对微尺寸X 射线源的成像过程.这里假设待测信号是分布在三维空间的I,C,F 三个独立的字母形貌光源,其轴向间距均为8 mm,水平间距为0.20 mm;每个独立字母光源占据的平面空间约为0.16 mm × 0.20 mm;在距离字母光源“C”50 mm 的位置上放置一片正弦波带片作为收光孔径,其最内环半径为0.08 mm,环带数为256 环,波带片总半径为1.28 mm;探测器位于距离波带片远离光源一面500 mm 的位置上,用于采集字母光源产生的非相干全息图.其中,探测器的探测面尺寸为16 mm × 16 mm,像素点数目为2000 × 2000,单位像素点尺寸为8 μm × 8 μm;字母光源辐射的X 射线能量设定为15 keV,上述波带片对该能量X 射线的焦距约为7.75×104mm.在本模拟成像系统的物像距和波带片参数下,其衍射效应可忽略不计,符合非相干全息的几何光学成像机制;模拟成像的三维示意图以及在上述参数下产生的非相干全息图如图2 所示,图3 是模拟成像系统二维侧视图.

图2 针对微尺寸X 射线源的成像示意图与二维全息图Fig.2.Imaging system for microscale X-ray source and the corresponding 2D hologram.

图3 模拟成像系统二维侧视图Fig.3.2D lateral view of the simulative imaging system.

本文采用早期Gabor 波带片编码成像实验用到的背传输算法[19],维纳滤波算法[21]以及前文理论部分介绍的压缩感知模型三种重建算法求解全息图.图4 展示了对本次模拟生成的非相干全息图的层析重建结果,可知上述3 种重建方法都能在确定的物距上反解出相应的字母光源,这说明非相干全息技术具有三维成像的潜力.对比这3 种算法的重建结果可以发现: 背传输算法计算出光源信号的边缘还存在着大量串联噪声,这会在实际应用中对信号识别造成影响;维纳滤波算法虽然能在一定程度上识别出待测图像信号,但由于重建出的图像对比度过低,使得目标源强分布混杂于背景噪声之中.不能有效地凸显重建有效信号的强度特征;而基于压缩感知算法的重建结果显示: 其非信号区域的串联噪声显著降低,光源物距平面上的光强分布基本上都集中在信号区域.除了存在光源信号的物距平面,我们还额外计算了距字母光源前后4 mm平面上的光源分布.从图4(b)和图4(c)的重建结果可以看出,利用背传输算法和维纳滤波算法计算得到的非信号平面上,仍然存在着大量来自其他平面上光源信号的非聚焦赝像.如前所述,在压缩感知重建的模型((7)式)中,引入TV norm 作为正则项,以此来表征重建图像在梯度域的稀疏度.由于上述非聚焦赝像无法在特定的物平面内重建成为具有明显边界的清晰图像,这些在梯度域中不稀疏的干扰项,会在求解(7)式极小值的过程中不断被消除.因此,通过压缩感知算法的重建,非信号平面内的背景噪声会随着算法迭代次数的增加而有效地被抑制.在图4(d)—(f)中,分别展示了TWIST 算法迭代100 次、500 次与2000 次的重建结果.

为了定量表征压缩感知重建后的图像质量,使用峰值信噪比(peak signal-to-noise ratio,PSNR,单位dB)作为评价指标,其表达式定义为

其中,Omax是图像单个像素点的最大数值.MSE是均方误差,用于反映原始数据和重建数据之间的差异,其表达式可以写作:

其中,o和o*分别是源强分布矩阵和重建后的强度分布矩阵,Nx代表图像矩阵的总行数,Ny代表图像矩阵的总列数,Nz代表图像矩阵的总页数.由(9)式可以看出,MSE 值越小说明重建后的图像与源光强分布越接近.而基于MSE 计算得出的PSNR,其数值越大表明图像的失真程度越小.图5展示了TWIST 算法迭代次数与重建图像PSNR值之间的关系.可以看出,PSNR 数值在算法运行初期经历微弱的震荡后快速上升,这表明图像质量在该迭代区间(Iterations ≤ 200)内被迅速优化,并在经历大量迭代次数后逐渐趋于平稳.

图5 PSNR 与重建迭代次数间的关系Fig.5.Relationship between PSNR and reconstruction iterations.

除了对非相干全息技术应用于层析成像的可行性研究,本文还研究了不同参数的波带片对成像分辨水平的影响.与波带片聚焦成像类似,非相干全息成像的水平分辨率也与波带片的最内环半径成反比(在一定环带数目的条件下).只是在非相干全息成像中用于参考的是波带片投影的最内环半径,它不仅与波带片编码孔径相关,还受到物像距关系的调制.图6 展示了选用两种不同最内环半径波带片对前述的字母光源进行非相干全息模拟成像的重建结果对比.基于图5 所示的图像质量随迭代次数变化关系,为了使用较好且稳定的图像重建质量进行对比,两种参数的波带片重建结果均经过TWIST 算法的2000 次迭代计算得出.重建结果包含两个部分;一部分是光源在物距为54,58和62 mm 平面处的二维重建结果,另一部分是位于58 mm 平面上虚线标识区域内的一维光强分布情况.用作对比的两种波带片的总环带数均为256 环,最内环半径分别为0.08 和0.06 mm.通过对比图6(a)和图6(b)的重建结果可以发现,在选用最内环半径为0.06 mm 的波带片进行成像后,在光源信号平面内的分辨率得到了明显提高,具体表现在一维光强分布中的信号上升和下降沿所经过的距离变窄.其中,图6(a)的虚线区域由上至下共经历4 个边沿,其距离分别为80,64,64 和72 μm,而在图6(b)中相应的边沿距离为40,40,32 和40 μm.但值得注意的是,使用小内环半径的波带片后在光源信号平面的非信号区域和非信号平面内重建产生了来自其他平面光源的非聚焦像,这意味着其轴向分辨率降低,非相干全息成像的重建算法无法正确复现以4 mm 作为轴间距的层析影像.因此,非相干全息成像的轴向和水平面内分辨率相互制约,在实际应用中需要根据需求选取特定参数的波带片作为编码孔径.

图6 使用不同最内环半径波带片进行模拟成像的分辨水平对比.在距离波带片54,58 和62 mm 平面上的二维重建结果,以及在58 mm 平面内虚线标识区域的一维光强分布情况 (a) 0.08 mm;(b) 0.06 mmFig.6.Comparisons of resolution level in simulative imaging when applied FZP with different innermost radius.The 2D reconstruction result at the following objective depth: 54,58 and 62 mm,and the 1D intensity distribution of the dotted line area in 58 mm plane: (a) 0.08 mm;(b) 0.06 mm.

4 结论

本文研究了非相干全息成像机制以及三维图像重建算法.为了验证该技术的可行性,数值模拟了微尺寸体光源的非相干全息成像过程,并使用压缩感知算法对其产生的全息图进行逐层反解,得到了体光源距离探测器不同深度位置上的二维光强分布情况.结果表明: 在文中所例举的参数体系下,非相干全息技术能够对微尺寸体光源进行层析成像;相较于传统的背传输算法和维纳滤波算法,基于压缩感知算法的重建结果能有效的消除来自不同物距上光源的非聚焦赝像干扰.综上,该项技术有望用于对ICF 中黑腔等离子体分布进行层析成像,完善ICF 领域中现有的成像诊断手段.