摘要:在基于光栅的X射线相衬信号的恢复方法中, 主要有相移法和傅里叶变换法两种方法. 相移法具有精度高、噪声小的优点, 但由于至少需要三幅图像才能恢复出相衬信号, 样品所受的辐射剂量大. 而傅里叶变换法只需一幅图像即可恢复出物体的相衬信号, 具有快速、实时的优点, 但恢复出的信号精度低, 易受伪影影响. 因此, 本文利用两幅图像傅里叶变换法恢复X射线相衬信号, 该方法能够有效地抑制相衬信号中由于频谱混叠所产生的伪影. 另外, 通过增加载波条纹的频率, 能够拉大频域中的频谱间隔, 从而进一步抑制伪影的产生.
关键词: X射线相衬成像/
傅里叶变换/
暗场成像

English Abstract

全文HTML

--> --> -->

X射线相衬成像是研究X射线通过物体后微弧度量级折射角的技术, 其在工业、材料和医学等领域具有潜在的应用价值^[1-8]. 在基于光栅的X射线相衬信号的恢复方法中, 傅里叶变换法是一种快速、实时的方法^[9]. 该方法在2008年由Wen等^[10,11]提出, 但恢复出来的图像较为模糊. 之后韩国的Lim等^[12-14]和德国的Seifert等^[15,16]在图像质量方面取得不错的进展. 但经过载波调制后的物体信息, 在进行傅里叶变换后不同级次的频谱会不可避免地发生混叠. 因此, 如果对频谱进行滤波的窗口宽度太小, 则高频信息会损失较多, 从而引起恢复出的相衬像变得模糊. 如果滤波的窗口宽度太大, 频谱容易发生混叠, 从而在恢复出的图像中引起严重的伪影. 对此, 本文利用两幅相差${\text{π}}$相位的图像进行傅里叶变换恢复X射线相衬信号^[17,18], 该方法能够在尽量保留物体高频信息的同时, 减少相衬信号中伪影. 和其他利用两幅图像恢复相衬信号的算法^[19-21]相比, 本文提出的算法对光栅在相移曲线中的位置没有特殊的要求, 从而降低了实验的操作难度.

2

-->

2.1.频谱混叠产生伪影的机理

利用傅里叶变换法恢复相衬信号的非相干X射线成像系统示意图如图1所示. 它由一个焦斑5 μm的点源S、周期96 μm的吸收光栅G1、物体O和像素尺寸为74.8 μm的平板探测器D所组成.
图 1 非相干X射线成像系统示意图
Figure1. Schematic diagram of incoherent X-ray imaging.

在一般情形下, 点源S照明吸收光栅G1在探测器平面上产生的载波条纹强度可以表示成:

$\begin{split}I\left( x \right) = \;&a\left( x \right) + b\left( x \right)\cos \left[ {2\pi {f_0}x + \phi \left( x \right)} \right]\\ =\;& a\left( x \right) + c\left( x \right)\exp \left( {{\rm{i}}2{\rm{\pi}} {f_0}x} \right) \\ &+ {c^*}\left( x \right)\exp \left( { - {\rm{i}}2{\rm{\pi}} {f_0}x} \right), \end{split}$

其中$a\left( x \right)$, $b\left( x \right)$和${f_0}$分别表示条纹的均值、振幅和频率; $\varphi \left( x \right)$表示样品所引起的相移, 其中$c\left( x \right)$为

$c\left( x \right) = \frac{1}{2}b\left( x \right)\exp \left[ {{\rm{i}}\phi \left( x \right)} \right].$

对(1)式做傅里叶变换, 可得一般情形下条纹的频谱:

$G\left( f \right) = A\left( f \right) + C\left( {f - {f_0}} \right) + {C^*}\left( { - f - {f_0}} \right),$

其中$G\left( f \right)$, $A\left( f \right)$和$C\left( f \right)$分别是$I\left( x \right)$, $a\left( x \right)$和$c\left( x \right)$的傅里叶变换. 在条纹的频谱中, $A\left( f \right)$, $C\left( {f - {f_0}} \right)$和${C^*}\left( { - f - {f_0}} \right)$分别代表零级、正一级和负一级频谱, 如图2所示.
图 2 一般情形下载波条纹的频谱图
Figure2. The Fourier spectrum of carrier fringe patterns in general case.

在一般情形下, 假定各级次频谱是完全分开的, 利用合适宽度和位置的滤波器分别与条纹不同级次的频谱相乘, 可以得到不同级次的频谱. 在得到正一级频谱后, 做傅里叶逆变换可得

$\frac{1}{2}b\left( x \right)\exp \left[ {{\rm{i}}\phi \left( x \right)} \right]\exp \left( {{\rm{i}}2{\rm{\pi}} {f_0}x} \right).$

对光路中没有样品的条纹图重复上述操作可得

$\frac{1}{2}{b_{\rm{r}}}\left( x \right)\exp \left( {{\rm{i}}2{\rm{\pi}} {f_0}x} \right),$

其中${b_{\rm{r}}}\left( x \right)$表示没有样品时条纹的振幅, 将上述两式相除再取其角度即可得到物体的相位:

${\phi _1}\left( x \right) = {\rm{angle}}\left( {\frac{{b\left( x \right)\exp \left[ {{\rm{i}}\phi \left( x \right)} \right]\exp \left( {{\rm{i}}2{\rm{\pi}} {f_0}x} \right)}}{{{b_{\rm{r}}}\left( x \right)\exp \left( {{\rm{i}}2{\rm{\pi}} {f_0}x} \right)}}} \right),$

其中“angle”表示取角度. 然而, 我们在实验中发现, 由于物体的频谱相当地宽, 不同级次的频谱会不可避免地发生混叠, 如图3所示.
图 3 实际情况下发生的频谱混叠
Figure3. Spectrum aliasing between different harmonic peaks in practice.

当频谱发生混叠时, 各个不同级次的频谱将带有其他级次频谱的分量. 由于零级频谱在整个频谱图中所占能量最大, 因此零级频谱对一级频谱所造成的混叠也最大. 此时, 将正一级频谱提取出来并做傅里叶逆变换后可得

${a_1}\left( x \right) + \frac{1}{2}b\left( x \right)\exp \left[ {{\rm{i}}\phi \left( x \right)} \right]\exp \left( {{\rm{i}}2{\rm{\pi}} {f_0}x} \right), $

其中${a_1}$是滤波后残留在正一频的零级频谱分量. 对无物体的载波条纹图重复上述操作后可得

${a_2}\left( x \right) + \frac{1}{2}{b_{\rm{r}}}\left( x \right)\exp \left( {{\rm{i}}2{\rm{\pi}} {f_0}x} \right),$

其中${a_2}$是滤波后残留在正一频的零级频谱分量, (7)式除以(8)式后取其角度即可得到有频谱混叠时恢复得到的相位:

$\begin{split}&{\phi _{{\rm{er}}}}( x )= \\\;& {\rm{angle}} \left[ {\frac{{{a_1}( x ) \!+\! 0.5b( x )\exp [ {{\rm{i}}\phi ( x )} ]\exp ( {{\rm{i}}2{\rm{\pi}} {f_0}x} )}}{{{a_2}( x ) + 0.5{b_{\rm{r}}}( x )\exp ( {{\rm{i}}2{\rm{\pi}} {f_0}x} )}}} \right].\end{split}$

下面, 用数值模拟来对比有无频谱混叠对恢复出来相位的影响. 各数值计算参数如下: ${a_1} = 0.08$, ${a_2} = 0.09$, $b = 0.8$, ${b_r} = 1$, $\varphi \left( x \right) = 1$.
图4(a)—(c)分别是载波条纹周期与探测器像素尺寸的比值r为3, 4和5时, 利用(9)式恢复出来的相位. 当零级频谱和一级频谱发生混叠时, 恢复出来的相位是周期性的, 其周期和比值r相同. 如果零级频谱和一级频谱没有发生混叠, 即当${a_1} = $$ 0$, ${a_2} = 0$时, 利用(9)式恢复出的相位是准确且无周期的, 如图5所示. 据此, 我们可以初步判断出, 零级频谱和一级频谱混叠是傅里叶变换法中产生伪影的一个重要原因.
图 4 频谱混叠对恢复出相位的影响　(a), (b)和(c)分别代表载波条纹周期与探测器像素尺寸比值r为3, 4和5时的情形
Figure4. The impact of spectrum aliasing on phase retrieval. (a), (b) and (c) denote the cases, in which the ratios of the carrier fringe period to size of detector pixel are 3, 4 and 5, respectively.

图 5 无频谱混叠时恢复出来的相位分布
Figure5. The phase retrieval when no spectrum aliasing.

2 -->

2.2.伪影抑制算法

前面我们指出, 零级频谱和一级频谱的混叠, 是傅里叶变换法产生伪影的重要原因. 因此, 为了减少频谱混叠, 需要尽量减少零级频谱对一级频谱的影响. 一般可以通过减小滤波窗口的宽度以抑制零频对一频的影响. 但该方法的主要问题是, 由于滤波窗口宽度的减小, 高频信息会损失较多, 从而导致恢复出的图像比较模糊. 对此, 我们采用两幅${\text{π}}$相差图像的傅里叶变换来消除零频对一频的影响, 该方法可以在尽量保留高频信息的同时, 大幅减少伪影的出现. 该方法需要采集两幅有物体的图像和两幅无物体的图像, 数据采集过程为: 先采集一幅背景图像${I_{1{\rm{b}}}}$, 然后在光路中放上物体, 采集一幅带物体的图像${I_{1{\rm{o}}}}$. 接着垂直于栅线方向移动光栅半个周期的距离, 然后采集第二幅带物体的图像${I_{2{\rm{o}}}}$. 随后把物体从光路中拿走, 采集第二幅没有物体的图像${I_{2{\rm{b}}}}$. 上述采集方式和多步相移法^[22]相比, 没有空程差, 因此可以降低误差. 该伪影抑制方法推导过程如下:
在第一幅有物体的图像${I_{1{\rm{o}}}}$中, 条纹分布可写成:

${I_{{\rm{1o}}}}\left( x \right) = a\left( x \right) + b\left( x \right)\cos \left[ {2{\rm{\pi }}{f_0}x + \phi \left( x \right)} \right].$

由于第二幅有物体的图像${I_{2{\rm{o}}}}$和第一幅有物体的图像${I_{1{\rm{o}}}}$差了半个光栅周期, 因此条纹之间的相位差为${\text{π}}$:

$\begin{split}{I_{2{\rm{o}}}}\left( x \right) =\;& a\left( x \right) + b\left( x \right)\cos \left[ {2{\rm{\pi }}{f_0}x + \phi \left( x \right) + {\rm{\pi }}} \right]\\ =\;& a\left( x \right) - b\left( x \right)\cos \left[ {2{\rm{\pi }}{f_0}x + \phi \left( x \right)} \right].\end{split}$

将(10)式减去(11)式即可消除零频对恢复信号的影响, 得到无零频的条纹分布:

${I_{\rm{o}}} = 2b\left( x \right)\cos \left[ {2{\rm{\pi }}{f_0}x + \phi \left( x \right)} \right].$

对(12)式做傅里叶变换, 然后提取正一频, 再做傅里叶逆变换, 即可得到正一频的图像:

$b\left( x \right)\exp \left[ {{\rm{i}}\phi \left( x \right)} \right]\exp \left( {{\rm{i}}2{\rm{\pi }}{f_0}x} \right).$

对两幅无物体的图像${I_{1{\rm{b}}}}$和${I_{2{\rm{b}}}}$重复上述操作可得无物体的正一频图像:

${b_{\rm{r}}}\left( x \right)\exp \left( {{\rm{i}}2{\rm{\pi }}{f_0}x} \right).$

最后用(13)式除以(14)式后取角度, 即可恢复出伪影得到抑制的相衬信号, 即:

${\phi _2} ( x ) \!=\! {\rm{angle}} \left[ {\frac{{b ( x )\exp [ {{\rm{i}}\phi ( x )} ]\exp ( {{\rm{i}}2{\rm{\pi }}{f_0}x} )}}{{{b_r} ( x )\exp ( {{\rm{i}}2{\rm{\pi }}{f_0}x} )}}} \right].$

2

-->

3.1.两幅图像傅里叶变换法抑制伪影的效果

下面通过实验来检验上述两幅图像傅里叶变换法抑制伪影的效果. 实验装置示意图如图1所示, X射线源到探测器的距离为80 cm, 光栅离探测器的距离为45.78 cm, 此时光栅在探测器平面上的投影条纹周期与探测器像素尺寸的比值r = 3. X射线源的电压为40 kV_p, 电流为160 μA, 在每个光栅位置采集三幅图像取平均以降低噪声, 每幅图像曝光4 s, 样品为直径5 mm的有机玻璃棒(PMMA), 其离探测器的距离为53 cm. 两幅图像傅里叶变换法和单幅图像傅里叶变换法恢复出来的PMMA相衬像如图6所示. 滤波窗口都为高斯窗, 其标准差$\sigma $等于零频和一频间距的二分之一. 图6(a)是两幅图像傅里叶变换法恢复出的相衬像, 该相衬像中没有出现伪影. 而图6(b)则是单幅图像傅里叶变换法恢复出的相衬像, 该相衬像则有明显的伪影. 从图6的直观比较中可以看出, 两幅图像傅里叶变换法所恢复出的相衬像(图6(a))比单幅图像傅里叶变换法所恢复出的相衬像(图6(b))在抑制伪影方面有明显的改进.
图 6 (a) 两幅图像傅里叶变换法所恢复出的PMMA相衬像; (b)单幅图像傅里叶变换法所恢复出的PMMA相衬像
Figure6. (a) The phase-contrast image of PMMA retrieved by Fourier transform with two images; (b) the phase-contrast image of PMMA retrieved by Fourier transform with one image.

为了进一步比较两种算法对伪影的抑制效果, 需要进行定量的分析. 下面我们截取图6(a)中的白色方框区域, 对其进行按行求平均, 绘制得到图7(a)中的蓝色曲线; 在图6(b)中截取与图6(a)中白色方框所示的相同位置的区域, 同样按行取平均, 得到图7(a)中的红色曲线. 同时取图6(a)中的黑色矩形框区域, 对其进行按行求平均, 绘制得到图7(b)中的蓝色曲线; 对图6(b)取相同位置的区域, 同样按行取平均, 得到图7(b)中的红色曲线. 从图7(a)的红色曲线可以看出, 利用单幅图像傅里叶变换法恢复出的相位值有明显的周期性, 且周期恰好为3个像素宽度, 这和图4(a)的仿真结果吻合地很好. 而利用两幅图像傅里叶变换法求得的相位则无明显的周期性, 且相位值更接近理论值零. 这意味着两幅图像傅里叶变换法不仅对伪影的抑制效果更好, 而且恢复出的相位更准确. 图7(b)中红色曲线的中间部分出现了明显的震荡, 而蓝色曲线的中间部分震荡不明显. 图7(b)中红色曲线的峰峰值较蓝色曲线的峰峰值高, 说明伪影不仅仅会造成背景信号偏高, 也会使得物体的相衬信号偏高.
图 7 (a) 图6中白色方框区域按行取平均后绘制的曲线; (b) 图6中黑色方框区域按行取平均后绘制的曲线
Figure7. (a) Curves from averaging the area of white rectangle in figure 6 by row; (b) curves from averaging the area of black rectangle in Figure 6 by row.

接着计算图7(a)中背景相位的平均值和标准差, 同时计算图7(b)中物体横截面曲线的峰峰值, 并将结果列于表1中.

	背景相位 均值/rad	背景相位标 准差/rad	横截面峰 峰值/rad
单幅图像 傅里叶变换	0.3502	0.0059	0.2412
两幅图像 傅里叶变换	0.2526	0.0017	0.1112

表1两种不同傅里叶变换方法的定量比较
Table1.Quantitative comparison between two kinds of Fourier transform algorithms.

从表1的定量数据中可以计算出, 单幅图像傅里叶变换法恢复的背景相位均值和标准差分别比两幅图像傅里叶变换法恢复的背景相位均值和标准差高出了38.6%和247%, 单幅图像傅里叶变换法恢复的物体横截面峰峰值比两幅图像傅里叶变换法恢复的物体横截面峰峰值高出了117%, 这就定性地证实了两幅图像傅里叶变换法恢复出来的相位要比单幅图像傅里叶变换法恢复出来的相位更加准确和稳定. 不过, 相对于单幅图像傅里叶变换法而言, 两幅图像傅里叶变换法增加了曝光时间和辐射剂量, 不利于快速成像. 但相对于相移法来说, 两幅图像傅里叶变换法的曝光时间仍然可以降低至少33%.
2 -->

3.2.不同频谱间隔对伪影产生的影响

当频谱间隔较小时, 两幅图像傅里叶变换法无法完全抑制频谱混叠产生的伪影. 因此, 选择合适的频谱间隔对于抑制伪影是非常重要的. 下面我们通过实验对比不同频谱间隔对伪影产生的影响. 实验装置示意图如图1所示, X射线源的电压为50 kV_p, 电流为160 μA. X射线源离探测器80 cm, 样品鸡翅离探测器12 cm. 采集光栅在两个不同位置的图像, 每个位置采集三幅图像取平均去噪, 每幅图像曝光4 s. 实验结果如图8所示.
图 8 载波条纹周期与探测器像素尺寸比值r = 3时, 鸡翅的频谱　(a)、相衬像(c)和暗场像(e). r = 5时, 鸡翅的频谱(b)、相衬像(d)和暗场像(f)
Figure8. Fourier spectrum (a), phase-contrast image (c) and dark-field image (e) of a chicken wing when r = 3. Fourier spectrum (b), phase-contrast image (d) and dark-field image (f) of a chicken wing when r = 5.

图8(a)、图8(c)和图8(e)分别是载波条纹周期与探测器像素尺寸比值r = 3时, 鸡翅的频谱图、相衬像和暗场像, 此时光栅离探测器45.78 cm. 图8(b)、图8(d)和图8(f)分别是r = 5时, 鸡翅的频谱图、相衬像和暗场像, 此时光栅离探测器59.47 cm. 平板探测器共有864 × 1536个像素, 图8(a)中一级频谱中心和零级频谱中心间隔288个像素, 而图8(b)中一级频谱中心和零级频谱中心间隔173个像素. 图8(a)和图8(b)使用的滤波窗口均为标准差$\sigma $ = 144个像素的高斯窗. 可以看出, 使用频谱间隔较宽的频谱图恢复出的相衬像图8(c)相对于使用频谱间隔较窄的频谱图恢复出的相衬像图8(d), 伪影明显减少. 这说明在相同的滤波窗口宽度下, 拉大频谱的间隔能有效地减少频谱混叠. 然而, 频谱的间隔并不能无限大, 其最大间隔为频谱扩展方向上像素数量的一半. 从图8(e)与图8(f)的比较中可以看出, 拉大频谱间隔对于减少暗场像的伪影同样有效, 这点可以从将(9)式中的取角度操作改为取模操作得知. 将(9)式中的取角度操作改为取模操作后, 我们将获得和暗场像相关的条纹振幅的比值, 该比值同样会受到滤波后残留的零频分量${a_1}\left( x \right)$, ${a_2}\left( x \right)$的影响.

针对单幅图像的傅里叶变换法恢复物体相位信息过程中容易出现伪影的问题, 本文分析了伪影产生的机理, 提出了抑制伪影的算法. 在利用单幅图像的傅里叶变换恢复物体相位信息过程中, 由于物体零级频谱和一级频谱的混叠, 造成一级频谱滤波后仍残留有零级频谱分量, 最终导致恢复出来的相位图像中出现了和载波条纹相同频率的伪影. 据此, 我们利用两幅相位差为${\text{π}}$的图像进行傅里叶变换去除零频, 最终有效地抑制了伪影的产生. 同时, 通过增大载波条纹的频率以拉大频谱间隔, 可以进一步地抑制频谱混叠和伪影. 本文提出的两幅图像傅里叶变换算法虽然是在非相干X射线成像系统上得到验证的, 但对于Talbot-Lau干涉仪, 理论上同样可以应用该算法进行伪影的抑制.