摘要:X射线光栅微分相位衬度成像技术可以观察到常规吸收衬度成像难以分辨的弱吸收物质的精细结构信息, 因而在医学、材料学等研究领域具有巨大的应用前景. 但传统的X射线光栅微分相位衬度成像技术由于采用分析光栅作为空间滤波器, 需要采用相位步进法扫描分析光栅来获得样品的多张投影图像才能够分离出样品的吸收、折射和散射信息, 因此存在样品曝光时间长、辐射剂量高以及X射线光通量利用率低等问题, 限制了其在各个学科领域的应用研究. 为克服上述问题, 本文提出一种基于免分析光栅相位衬度成像系统的一次曝光样品信息提取算法. 该算法只需要利用一块相位光栅, 进而采用高分辨探测器进行样品投影数据的一次采集即可提取样品的吸收、折射和散射信息. 理论和模拟研究结果表明: 与传统相位步进法相比, 该算法具有样品信息提取精度高, 且不受光栅的自成像周期需为探测器像素尺寸的整数倍条件的限制. 此外, 该算法还能够有效地减少对生物样品的辐射损伤, 因此在生物医学成像等研究领域中具有广泛的应用前景.
关键词: X射线相衬成像/
相位光栅/
一次曝光/
信息分离

English Abstract

全文HTML

--> --> -->

X射线自被伦琴发现以来就被广泛应用于医学影像诊断以及材料检测等领域. 但是传统的X射线成像技术主要是基于吸收衬度成像机制, 对由碳、氢、氧等元素构成的弱吸收物质的成像质量不佳. 为了克服上述困难, 发展基于X射线相位衬度成像理论的新型X射线成像技术成为近年来的研究热点. 其中, X射线光栅微分相位衬度成像技术由于对X射线光源的相干性要求较低, 可以与实验室光源相结合实现弱吸收物质的大视场、高分辨相位衬度成像, 因此被认为是最有可能在临床医学中获得应用的一种X射线相位衬度成像方法, 近年来获得了广泛的发展^[1-3]. 但是传统X射线光栅微分相位衬度成像技术需要采用分析光栅作为空间滤波器使用. 一方面, 分析光栅会阻挡X射线, 降低了整个成像系统X射线的利用效率. 另一方面, 分析光栅需要进行步进扫描来获取样品的多张投影图像才能提取样品的吸收、折射和散射信息, 因此存在样品曝光时间长、辐射剂量高等问题. 上述这些问题的存在都限制了X射线光栅微分相位衬度成像技术的应用范围, 特别是生物医学成像等对辐射剂量敏感的研究领域. 因此如何解决上述X射线光栅微分相位衬度成像技术发展过程中存在的X射线的利用效率低和样品辐射剂量高的问题, 从而实现快速、低剂量的X射线光栅微分相位衬度成像技术, 成为近年来X射线相位衬度成像领域发展的重要方向.
实现上述快速、低剂量X射线光栅微分相位衬度成像技术的发展需求可以通过两种途径来实现. 第一种途径是改造光栅的结构, 避免分析光栅的步进扫描过程从而减少样品的辐射损伤; 第二种途径是发展基于免分析光栅X射线相位衬度成像系统的一次曝光相位衬度成像技术. 在对第一种途径的探索中, Endrizzi等^[4]、Fu等^[5]、Wei等^[6]以及Ge等^[7]分别通过改造光栅的结构实现了样品一次曝光成像. 但是这种方法仅仅是避免了分析光栅的步进扫描过程, 减少了样品的辐射损伤. 其仍需要使用分析光栅, 因此无法解决分析光栅对X射线的吸收导致的整个成像系统X射线利用率低的问题. 此外经过改造的光栅结构还降低了成像装置在纵向上的分辨率. 而基于免分析光栅X射线相位衬度成像系统的一次曝光相位衬度成像技术则完全避免了分析光栅的使用. 其不需要分析光栅进行步进扫描来进行样品的数据采集, 只需要进行样品的一次投影数据采集就能够提取样品的吸收、折射和散射信息. 该方法可以有效地提高整个成像系统X射线利用率的同时解决传统X射线光栅微分相位衬度成像技术中样品曝光时间长、辐射剂量高的问题. 因此通过第二条途径实现上述快速、低剂量的X射线光栅微分相位衬度成像技术的发展需求更为有效. Balles^[8]、Marathe^[9]以及Wen^[10]、Bennett^[11]等验证了第二种途径的可行性, 此后基于免分析光栅X射线相位衬度成像系统的一次曝光成像理论的发展受到了广泛的关注. 但是目前基于该成像理论的样品信息提取算法还不完善. 例如Balles等^[8]和Marathe等^[9]使用的相位步进信息提取算法(phase-stepping algorithm, PS算法)虽然能够得到准确的样品吸收、折射和散射信息, 但PS算法需要满足光栅自成像期等于探测器像素尺寸的整数倍的限制条件. 然而这一限制条件在实际实验过程中是难以严格满足的, 其光栅自成像周期和探测器像素尺寸之间的差异会严重影响样品信息提取的精度, 进而影响X射线相位衬度成像中定量化分析的准确性和有效性; 而Wen等^[10]和Bennett等^[11]所使用的样品信息提取算法是对单张投影图直接进行空间谐波分析, 与PS算法相比较, 该算法会显著降低样品图像信息的空间分辨率.
针对基于免分析光栅X射线相位衬度成像系统中一次曝光成像理论在样品信息提取算法上存在的上述问题, 本文开展了基于免分析光栅X射线相位衬度成像系统的样品信息提取算法的研究, 提出一种基于超定方程理论的最小二乘法的样品信息提取算法. 该算法的适用性并不受光栅周期和探测器单个像素尺寸的匹配性限制. 此外, 该算法在免分析光栅X射线相位衬度成像系统中可以获得与传统PS算法具有相同空间分辨率的样品的吸收、折射和散射图像信息.

2

-->

2.1.基于免分析光栅X射线相位衬度成像系统的一次曝光相位衬度成像原理

基于免分析光栅X射线相位衬度成像系统装置如图1所示, 主要包括满足一定相干条件的X射线源、相位光栅${{\rm{G}}_1}$($\pi $相位光栅或$\pi /2$相位光栅)、以及探测器像素尺寸为$b \times b$的高分辨探测器. 本文中为便于描述, 在免分析光栅X射线相位衬度成像系统装置中相位光栅${{\rm{G}}_1}$选取为$\pi $相位光栅. 如图1(a)所示, 由于Talbot效应相位光栅${{\rm{G}}_1}$在Talbot距离D处会形成${{\rm{G}}_1}$光栅的自成像条纹${{{\rm{G}}'_1}}$, 其中自成像条纹${{\rm{G}}'_1}$的周期p为$\pi $相位光栅${{\rm{G}}_1}$周期${p_1}$的一半^[12]. 高分辨探测器则位于${{\rm{G}}_1}$的自成像位置处, 可以直接对自成像条纹进行图像采集.
图 1 基于免分析光栅X射线相位衬度成像系统的装置成像示意图
Figure1. Schematic imaging diagram of the phase contrast imaging system without analyser grating.

在实验过程中, 如图2所示当高分辨率探测器的像素尺寸b等于$p/N$时(N为正整数), 可使用高分辨探测器分别采集一张样品的投影图像和背景的投影图像. 此后将得到的样品投影图像沿${x_{\rm{g}}}$方向每N列抽取一列像素组成一张新投影图像, 则可获得N张新投影图像, 分别定义为${I_1}$, ${I_2}$, ···, ${I_N}$. 当N足够大时, 上述抽取方法获得的N张新投影图像中相同像素位置处的光强可以组成图2所示的光强变化曲线. 而当$N = 4$时可以获得${I_1}$ (黑色), ${I_2}$ (红色), ${I_3}$ (绿色), ${I_4}$ (黄色)共四张新投影图像.
图 2 由采集到的投影图像组成的N张新投影图像中同一像素点的光强随X射线的偏折角$\psi = {{{x_{\rm{g}}}} / D}$的变化曲线(角度曲线)
Figure2. Plot of the new images’ intensity oscillation (shifting curve) of single pixel as a function of the deviation angle $\psi = {{{x_{\rm{g}}}} / D}$.

这些新投影图像等效于传统的含分析光栅的X射线光栅微分相位衬度成像装置中, 使用占空比为$1/N$的分析光栅沿${x_{\rm{g}}}$方向做步进扫描运动, 并以${p / N}$为步进间距进行N次采集得到的投影图像. 因此传统的含分析光栅的X射线光栅微分相位衬度成像技术和图1所示的免分析光栅的相位衬度成像技术的原理是类似的, 都是对自成像条纹进行N点采样. 根据上述分析, 则图2中有样品时的新投影图像的光强可以表示为^[13,14]:

${I_\eta }\left( {\psi ;x,y} \right) = f\left( {\psi ;x,y} \right) \otimes {S'_\eta }\left( \psi \right),$

其中$ \otimes $表示卷积; $f\left( {\psi ;x, y} \right)$表示样品对X射线的角度调制函数, 包含吸收、折射和散射效应; ${S'_\eta }\left( \psi \right)$是无样品时的N张新投影图像的光强所形成的背景角度曲线, $ \eta =1, \;2, \;\cdots , N$表示一组新投影图像的序号; ${S'_\eta }\left( \psi \right)$与传统X射线光栅微分相位衬度成像系统中通过吸收光栅做步进扫描所获得的系统位移曲线类似; ${I_\eta }\left( {\psi ;x, y} \right)$是有样品时的N张新投影图像中相同像素点$\left( {x, y} \right)$上的光强变化曲线. 此处${S'_\eta }\left( \psi \right)$与${I_\eta }\left( {\psi ;x, y} \right)$均可以采用余弦函数来拟合^[14-16]. ${S'_\eta }\left( \psi \right)$可表示为

${S'_\eta }(\psi) \approx {\bar {\rm{S}}}\Big[ {1 - {V_0}\cos \Big( {\frac{{2\pi D}}{p}( {\psi - {\psi _0}} )} \Big)} \Big], $

其中${\psi _0}$表示背景角度曲线的初始角度; $\bar S = ( {S_{\max }} + {S_{\min }}) / 2$表示背景角度曲线的平均值; ${S_{\max }}$和${S_{\min }}$分别表示背景角度曲线的最大值和最小值; ${V_0} = {{\left( {{S_{\max }} - {S_{\min }}} \right)} / {\left( {{S_{\max }} + {S_{\min }}} \right)}}$表示背景角度曲线的可见度. 而样品对X射线的角度调制函数$f\left( {\psi ;x, y} \right)$可表示为^[17]

$f\left( {\psi ;x,y} \right) =\! \frac{{\exp \left( { \!-\! M\left( {x,y} \right)} \right)}}{{\sqrt {2\pi } \sigma \left( {x,y} \right)}}\exp \!\bigg(\! - \frac{{{{( {\psi \!- \!\theta ( {x,y} )} )}^2}}}{{2{\sigma ^2}({x,y})}}\! \bigg), $

其中M代表样品对X射线的吸收; $\theta $代表X射线经过样品后的折射角; ${\sigma ^2}$表示X射线经过样品后散射角度高斯分布的方差. 根据(1)式—(3)式, 则免分析光栅X射线相位衬度成像系统中探测器上的样品光强方程可表示为

$\begin{split}\;& {{ {I}}_\eta } = I\left( {\eta \frac{{2\pi }}{p};x,y} \right) \\ =\;& \bar S\exp \left( { - M\left( {x,y} \right)} \right)\bigg\{ 1 - {V_0}{{\rm{e}}^{ - 2{\pi ^2}{D^2}{\sigma ^2}(x,y)/{p^2}}}\\ &\times\cos \Big(\frac{{2\pi D}}{p}(\theta (x,y) \!+\! {\psi _0}(x,y)) - \eta \frac{{2\pi }}{p} \cdot b \Big) \bigg\}. \end{split} $

在实际实验过程中, 相位光栅自成像周期p是难以满足为探测器像素尺寸b的N倍的要求的, 因此为了保证样品信息提取的精度, 需要在角度曲线余弦拟合中考虑探测器的像素尺寸b和光栅自成像周期p之间的差值, 这里令探测器像素尺寸b与光栅自成像的周期的$1/N$的差值所对应的角度曲线上的相位差为$\alpha $, 即

$\alpha = {{2\pi \left( {b - {p / N}} \right)}}/{p}.$

2 -->

2.2.信息提取算法

将(5)式代入成像方程式(4)可以得到N张样品光强图的表达为

$\begin{split}{{\rm{I}}_\eta }\left( {x,y} \right) =\;& a\left\{ {1 + V\cos \left( {\varphi - \left( {\frac{{2\pi }}{N} + \alpha } \right) \cdot \eta } \right)} \right\},\\\eta =\;& 1,2, \cdots,N.\\[-10pt]\end{split}$

令

$ \begin{split}&a=\overline{S}{\rm{exp}}\left[-M\left(x,y\right)\right]\;, \\ &\varphi =\frac{2\pi D}{p}\left[\theta \left(x,y\right)+{\psi }_{0}\left(x,y\right)\right], \\&V=-{V}_{0}{\rm{e}}^{-2{\pi }^{2}{D}^{2}{\sigma }^{2}\left(x,y\right)/{p}^{2}}, \end{split}$

以及

$ \begin{split}&{x}_{1}=a, ~~ {x}_{2}=aV{\rm{cos}}\varphi , ~~ {x}_{3}=aV{\rm{sin}}\varphi.\end{split}$

联立(6)式中的N个方程, 则(6)式又可以写为

$\begin{split}&\left[ \!\!{\begin{array}{*{20}{c}} 1&{\cos \left( {\left( {\dfrac{{2\pi }}{N} \!+\! \alpha } \right) \cdot 1} \right)}&{ \!-\! \sin \left( {\left( {\dfrac{{2\pi }}{N}\! +\! \alpha } \right) \cdot 1} \right)} \\ 1&{\cos \left( {\left( {\dfrac{{2\pi }}{N} \!+\! \alpha } \right) \cdot 2} \right)}&{ \!-\! \sin \left( {\left( {\dfrac{{2\pi }}{N}\! +\! \alpha } \right) \cdot 2} \right)} \\ \vdots & \vdots & \vdots \\ 1&{\cos \left( {\left( {\dfrac{{2\pi }}{N} \!+\! \alpha } \right) \cdot \eta } \right)}&{ \!-\! \sin \left( {\left( {\dfrac{{2\pi }}{N} \!+\! \alpha } \right) \cdot \eta } \right)} \end{array}} \!\!\right]\\&\times\left[ {\begin{array}{*{20}{c}} {{x_1}} \\ {{x_2}} \\ {{x_3}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{I_1}} \\ {{I_2}} \\ \vdots \\ {{I_\eta }} \end{array}} \right], \\[-35pt]\end{split}$

进而简写为

${{A}}\hat {{x}} = \hat {{I}}, $

其中,

$\hat {{x}} = \left[ {\begin{array}{*{20}{c}} {{{\hat x}_1}} \\ {{{\hat x}_2}} \\ {{{\hat x}_3}} \end{array}} \right], $

$\hat {{I}} = \left[ {\begin{array}{*{20}{c}} {{I_1}} \\ {{I_2}} \\ \vdots \\ {{I_\eta }} \end{array}} \right].$

(10)式中的A是由采样点数N以及探测器的像素尺寸b以及相位光栅自成像的周期p决定的常数矩阵. 由于$N > 3$, 故简化后的样品投影图像的表达式(10)式是一个超定线性方程组, 为了解(10)式, 可以将其转换为一个最小二乘优化问题:

${{x}} = \mathop {\arg \min }\limits_{\hat {{x}}} {\big( {{{A}}\hat {{x}} - \hat {{I}}} \big)^2},$

其解为

${{x}} = {\rm{pinv}} \left( {{A}} \right) \cdot \hat {{I}} = \left[ {\begin{array}{*{20}{c}} {{x_1}} \\ {{x_2}} \\ {{x_3}} \end{array}} \right],$

(14)式中的${\rm{pinv}} \left( {{A}} \right)$表示矩阵A的伪逆. 将该解中的${x_1}$, ${x_2}$和${x_3}$代入(8)式可以得到:

$\begin{split}& a = {x_1}, ~~ \varphi = \arctan \left( {\frac{{{x_3}}}{{{x_2}}}} \right), ~~ V = \frac{{\sqrt {x_2^2 + x_3^2} }}{{{x_1}}}{\;}{.} \end{split} $

同理, 当没有放置样品时, 可以将探测器上获得的背景投影图像分为N张光强图像:

${\hat {{I}}_{\rm{bg}}} = \left[ {\begin{array}{*{20}{c}} {{I_{{\rm{bg}}1}}} \\ {{I_{{\rm{bg}}2}}} \\ \vdots \\ {{I_{{\rm{bg}}N}}} \end{array}} \right], $

其解为

${{{x}}_{\rm{bg}}} = {\rm{pinv}} \left( {{A}} \right) \cdot {\hat {{I}}_{\rm{bg}}} = \left[ {\begin{array}{*{20}{c}} {{x_{{\rm{bg}}1}}} \\ {{x_{{\rm{bg}}2}}} \\ {{x_{{\rm{bg}}3}}} \end{array}} \right].$

由于无样品时吸收${M_{\rm{bg}}}$, 折射${\theta _{\rm{bg}}}$和散射${\sigma _{\rm{bg}}}$都为0, 故可得

$\begin{split} &{a_{\rm{bg}}} = {x_{{\rm{bg}}1}} = \bar S , \\& {\varphi _{\rm{bg}}} = \arctan \left( {\frac{{{x_{{\rm{bg}}3}}}}{{{x_{{\rm{bg}}2}}}}} \right) = \frac{{2\pi D}}{p}{\psi _0}\left( {x,y} \right), \\& {V_{\rm{bg}}} = \frac{{\sqrt {x_{{\rm{bg}}2}^2 + x_{{\rm{bg}}3}^2} }}{{{x_{{\rm{bg}}1}}}} = - {V_0}{\;}{.} \end{split} $

因此, 结合背景扣除^[18,19], 并将(7)式分别代入(15)式和(18)式得扣除背景后样品的吸收、折射和散射信息提取公式分别为

$\begin{split} & M = - \ln \left( {\frac{{{x_1}}}{{{x_{1{\rm{bg}}}}}}} \right), \\ &\theta = \frac{p}{{2\pi D}}\left[ {\arctan \left( {\frac{{{x_3}}}{{{x_2}}}} \right) - \arctan \left( {\frac{{{x_{3{\rm{bg}}}}}}{{{x_{2{\rm{bg}}}}}}} \right)} \right], \\& {\sigma ^2} = - \dfrac{{{p^2}}}{{2{\pi ^2}{D^2}}} \left[ \ln \left(\dfrac{ \dfrac{{\sqrt {{x_2}^2 + {x_{\rm{3}}}^2} }}{{{x_1}}} }{ \dfrac{{\sqrt {{x_{2{\rm{bg}}}}^2 + {x_{{\rm{3}}{\rm{bg}}}}^2} }}{{{x_{1{\rm{bg}}}}}} } \right) \right]. \end{split} $

其中$\arctan $是$\tan$的反三角函数.
下文中为描述方便, 将上述样品信息提取算法称为基于超定方程理论的最小二乘法的样品信息提取算法(least-squares algorithm, LS算法). 而传统的光栅微分相位衬度成像中N步等间距平移分析光栅的相位步进法对应的信息提取公式为^[15,19,20]

$\begin{split}& M = - \ln \left( {{{{a_0}}}/{{{a_0}_{{\rm{bg}}}}}} \right), ~~ \theta = \frac{p}{{2{\rm{\pi }}D}}\left( {\varphi - {\varphi _{{\rm{bg}}}}} \right), \\& {\sigma ^2} = - \frac{{{p^2}}}{{2{\pi ^2}{D^2}}}\left[ {\ln \left( {{V / {{V_{\rm{bg}}}}}} \right)} \right]{\;}{.} \end{split} $

其中

$\begin{split} &{a_0} = \frac{{\rm{1}}}{N}\mathop \sum \limits_{\eta = 1}^N {I_\eta }, ~~ \varphi = {\rm{arctan}}\left( {\frac{{\displaystyle \sum \limits_{\eta = 1}^N {I_\eta }{s_\eta }}}{{ \displaystyle\sum \limits_{\eta = 1}^N {I_\eta }{c_\eta }}}} \right), \\& V = \frac{2}{{{a_0}N}}\sqrt {{{\left( {\displaystyle \sum \limits_{\eta = 1}^N {I_\eta }{s_\eta }} \right)}^2} + {{\left( {\mathop \sum \limits_{\eta = 1}^N {I_\eta }{c_\eta }} \right)}^2}} , \\ &{s_\eta } = \sin \left( {2\pi \frac{\eta }{N}} \right),{\;}{c_\eta } = \cos \left( {2{\rm{\pi }}\frac{\eta }{N}} \right){\;}{.} \end{split} $

2

-->

3.1.系统参数设计

模拟实验装置如图1所示, X射线的能量为$E = 12{\;}{\rm{keV}}$. 选用的$ \rm{\pi}$相位光栅周期${p_1}$为4.80 μm, 探测器放置在一阶Talbot距离D为27.88 mm处. 实验中设计了一个圆柱体和三张厚度不同的纸作为实验样品, 如图3所示. 其中圆柱体的材质为PMMA, 其高为$r = 0.61\;{\rm{ mm}}$, 对X射线具有吸收和折射效应, 其复折射率虚部$\beta $为$1.79 \times {10^{ - 9}}$、实部$\delta $为${{1}}{{.85}} \times {\rm{1}}{{\rm{0}}^{{{ - 6}}}}$; 而纸张对X射线只有散射效应, 其散射宽度与纸的层数c有关, c层纸对X射线的散射角度方差为$c{\sigma ^2}$, 实验中${\sigma ^2}$取${{1}}{{.67}} \times {\rm{1}}{{\rm{0}}^{{{ - 11}}}}\;{\rm{ ra}}{{\rm{d}}^2}$. 模拟实验过程为: X射线穿过圆柱体时, 吸收效应会引起X射线光强的衰减, 同时折射效应会使得X射线偏离原来的入射方向, 折射角为$\theta $. 随后穿过圆柱体的X射线会照射在纸上, 此时X射线会以折射角为中心发生小角度的前向高斯散射, 此后小角度前向散射光会被成像探测器上的像素单元所接收. 当X射线穿过不同厚度的纸张时, 小角度前向高斯散射光的散射宽度不同(散射宽度为$4\sqrt {\left( {2\ln 2} \right) \cdot c} \cdot \sigma $), 因此会造成探测器上单个像素单元接收到的X射线光子数发生变化, 即X射线光强发生变化. 上述样品仿真模型的具体描述可参见参考文献[17,21].
图 3 模拟所用的样品示意图. 样品由直径为0.461 mm PMMA圆柱状和PMMA后不同厚度(0—3层)的纸组成, 其中PMMA圆柱状对X光具有吸收和折射效应, 而纸张对X光只有散射效应
Figure3. Schematic diagram of the simulation sample. The simulation sample has a 0.46 mm diameter PMMA cylinder that combined refraction and absorption effects. The PMMA cylinder overlie paper layers (0?3 layers) that exhibit the scattering effects.

在上述模拟过程中, 假设样品厚度远小于成像距离, 则可使用基于菲涅尔衍射理论来计算X射线经过样品及${{\rm{G}}_1}$光栅后在探测器上的光强分布图像^[19,22,23]. 此外, 为了比较探测器的像素尺寸与光栅自成像周期的匹配差对样品信息提取精度的影响, 可以分别使用两个不同像素尺寸的X射线高分辨成像探测器来对模拟样品进行数据采集. 这里假设两个探测器的总视场均为$r \times r$, 且探测器的像素尺寸分别为${b_1} \!=\! 0.60$ μm ($\alpha = 0$)和${b_2} \!=\! 0.66\;$ μm ($\alpha \ne 0$). 模拟得到的样品投影图像如图4(a),(c)所示. 其中图4(a)对应的探测器像素尺寸为0.60 μm, 图4(c)对应的探测器像素尺寸为0.66 μm, 从图4(a)右下角的放大图图4(b)可以看到, 当探测器的像素尺寸与光栅自成像的周期相匹配时, 探测器采集到的样品投影图像也为周期性分布, 且投影图像的周期与光栅自成像周期相等. 而从图4(c)右下角的放大图图4(d)可以观察到, 当探测器的像素尺寸与光栅自成像周期不匹配时, 匹配差的存在使投影图中出现了莫尔条纹.
图 4 探测器上的样品投影图　(a) 探测器像素尺寸为0.60 μm (α = 0)时模拟样品的投影图像; (b) 图(a)的局部(红框内)放大图; (c) 探测器像素尺寸为0.66 μm (α ≠ 0)时模拟样品的投影图像; (d) 图(c)的局部放大图
Figure4. The projective images of the simulation sample: (a) The projective image of the simulation sample with pixel size of 0.60 μm (α = 0); (b) local enlarged drawing of Fig. (a); (c) the projective image of the simulation sample with pixel size of 0.66 μm (α ≠ 0); (d) local enlarged drawing of Fig. (c).

使用(19)式描述的基于超定方程理论的样品信息提取算法将图3中所示的不同探测器采集到的样品投影图像以$N \cdot b$为抽取周期进行列像素抽取, 则可以分别重排为$N = 4$张样品投影图像. 如图5所示, 其中探测器的像素尺寸为0.60 μm时采集的样品投影图像所组合出的4张新的投影图像如图5(a)—(d)所示, 分别为样品处于角度曲线的谷位、上坡位、峰位、下坡位时的投影图像. 从图5(b)和图5(d)中可以观察到, 两幅图像在圆柱体的左右两边缘位置分别出现亮暗的现象, 其代表了折射角方向的不同, 该现象和传统含分析光栅的X射线光栅微分相位衬度成像装置中分析光栅做4次步进扫描后在位移曲线上坡、下坡位得到的样品投影图像是吻合的^[14]. 此外图5(a)和图5(c)分别为样品位于角度曲线的谷位和峰位时的投影图像, 此时由于角度曲线是对称分布的, 因此图5(a)和图5(c)样品边界处的光强分布也是对称的.
图 5 由原始的样品投影图组成的4张新图像　(a)—(d) 探测器像素尺寸为0.60 μm (α = 0)时所得到的4张新投影图; (e)—(h) 探测器像素尺寸为0.66 μm ($\alpha \ne {\rm{0}}$)时所得到的4张新投影图
Figure5. Four new images extracted from one original projective image: (a)?d) Contain four new images with the pixel size of 0.60 μm (α = 0); (e)?(h) contain four new images with the pixel size of 0.66 μm (α ≠ 0).

当采用探测器的像素尺寸为0.66 μm时, 采集到的样品投影数据经过抽取和重排后, 组合出新的4张样品投影图像如图5(e)—(h)所示, 其4张投影图像在${x_{\rm{g}}}$方向上均存在周期性的明暗条纹, 这是因为当探测器的像素尺寸与光栅自成像周期为非整数倍时, 按周期抽取后的同一投影图像中的多列像素不再分布在角度曲线的同一位置, 从而出现了莫尔条纹.
2 -->

3.2.样品的信息提取

利用PS以及LS样品信息提取算法, 对图5中的两组数据分别进行样品的吸收、折射和散射信息的提取. 当采用的探测器像素尺寸为0.60 μm时, 样品信息提取结果如图6所示. 其中, 图6(a)—(c)为使用PS算法提取得到的样品的吸收、折射和散射信息, 图6(d)—(f)为使用LS算法提取得到的样品的吸收、折射和散射信息. 可以看到图6(a)中由于圆柱体样品是厚度不均匀的, 中心区域的光强吸收较大, 使得其在中心区域样品的吸收信息更强. 而在折射信息图6(b)中, 圆柱边界区域的折射角大且两侧折射角方向相反, 使得折射图像两侧一亮一暗. 图6(c)所示为样品的散射信息, 由于纸的厚度从上到下依次变薄, 其散射系数也不断减小, 因此图6(c)中可以看到从上到下, 图像依次由亮变暗.
图 6 探测器像素尺寸为0.60 μm时模拟样品的信息提取结果　(a)—(c) PS算法时提取的吸收、折射和散射信息; (d)—(f) LS算法时提取的吸收、折射和散射信息; (g)—(i) 在虚线位置处(PS算法(绿色)和LS算法(紫色))的提取的吸收、折射和散射信息与理论值(蓝色)的对比图
Figure6. Sample information extracted with pixel size of 0.60 μm: (a)?(c) Depict the absorption, refraction and scattering information of the simulated sample obtained by PS algorithm; (d)?(f) depict the absorption, refraction and scattering information obtained by LS algorithm; (g)?(i) are profiles of the absorption, refraction and scattering images extracted by LS (purple) and PS (green) at the dotted lines in (a)?(f), as well as the theoretical values (blue).

为了更好的分析图6中的样品信息提取的准确性, 沿图6中虚线所示方向, 分别获取了样品的吸收信息、折射信息和散射信息的强度分布曲线. 如图6(g)—(i)所示, 当探测器的像素尺寸与光栅自成像的周期相匹配时, PS算法与LS算法所得的吸收、折射和散射基本一致, 且和相应的理论值能够较好地吻合.
当模拟中所用的探测器的像素尺寸为0.66 μm时, 样品信息提取结果如图7所示. 其中, 图7(a)—(c)使用PS算法提取得到的样品的吸收、折射和散射信息, 图7(d)—(f)为使用LS算法提取得到的样品的吸收、折射和散射信息. 图7(a)—(f)与图6(a)—(f)基本一致, 不过当使用PS算法进行样品信息提取时, PS算法提取的样品吸收、折射和散射图像中出现了与光栅方向平行的条状伪影, 如图7(a)—(c)所示. 同样的, 图7(a)—(f)中虚线所示的样品的吸收信息、折射信息和散射信息的强度分布曲线的比较结果如图7(g)—(i)所示, 可以发现当探测器的像素尺寸与光栅自成像的周期不匹配时, PS算法提取的吸收、折射和散射像中出现了与光栅方向平行的条状伪影, 所得到的吸收、折射和散射信息与理论值差异较大, 而LS算法所得的吸收、折射与散射信息则与理论值吻合较好. 另外, 由图7(h)和图7(i)图中可以看到, 折射信息和散射信息在圆柱样品边界区域的提取值与理论值之间存在比较明显的偏差.
图 7 探测器像素尺寸为0.66 μm时模拟用样品的信息提取结果　(a)—(c) PS算法时提取的吸收、折射和散射信息; (d)—(f) LS算法时提取的吸收、折射和散射信息; (g)—(i) 在虚线位置处(PS算法(绿色)和LS算法(紫色))的提取的吸收、折射和散射信息与理论值(蓝色)的对比图
Figure7. Sample information extracted with pixel size of 0.66 μm: (a)?(c) Depict the absorption, refraction and scattering information of the simulated sample obtained by PS algorithm; (d)?(f) depict the absorption, refraction and scattering information obtained by LS algorithm; (g)?(i) are profiles of the absorption, refraction and scattering images extracted by LS (purple) and PS (green) at the dotted lines in Fig. (a)?(f), as well as the theoretical values (blue).

为了进一步对LS算法和PS算法两种算法进行定量比较, 这里使用平均绝对误差(mean absolute error, MAE)来评估样品信息提取值e与理论值${e_0}$的误差情况, 其平均绝对误差可表示为

$\begin{split} {}& {\rm{MAE}} =\frac{{\displaystyle\sum\limits_{y = 1}^n {\displaystyle\sum\limits_{x = 1}^m \!{\left| {e\left( {x,y} \right) \!-\! {e_0}\left( {x,y} \right)} \right|} } }}{{mn}}; \\& x = 1,2,\cdots,m; ~~ y = 1,2,\cdots,n, \end{split}$

其中$m \times n$为样品信息值的总像素数. 根据(22)式计算出圆柱样品区域内的MAE误差的结果如表1所列. 在探测器像素尺寸与光栅自成像的周期匹配时, LS算法样品信息提取值的MAE误差与PS算法样品信息提取值的MAE误差基本相等, 且误差均较小, 说明两种算法提取的样品信息值和理论值均能较好地吻合. 当探测器的像素尺寸与光栅自成像的周期不匹配时, 探测器的像素尺寸为0.66 μm时, 两种算法样品信息提取值的MAE误差均有所增加, 但是LS算法样品信息提取值的MAE误差比PS算法样品信息提取值的MAE误差小将近1个量级, 说明在周期不匹配的情况下LS算法比PS算法提取得到的样品信息更加准确.

探测器的 像素尺寸		M/${10^{ - 4}}$			$\theta $/${10^{ - 9}}\;{\rm{ rad}}$			${\sigma ^2}/{10^{ - 13}}\;{\rm{ ra}}{{\rm{d}}^{\rm{2}}}$
		PS	LS		PS	LS		PS	LS
0.60 μm		2.10	2.10		4.53	4.53		1.36	1.36
0.66 μm		43.00	6.08		115.92	13.80		35.80	5.90

表1PS和LS算法的理论值和提取值的平均绝对误差
Table1.Mean absolute error of theoretical and extracted information by PS and LS algorithm.

此外, 图6(a),(c),(d),(f)和图7(a),(c),(d),(f)中样品边界处均出现了两条纵向的明暗相间的条纹状伪影, 这主要是由两个方面的原因造成的: 1) 本文采用了菲涅耳传播理论对X射线穿过样品的物理过程进行模拟, 随着X射线传播距离的增加, 在Talbot距离D处样品边界处出现了相位的二阶导数信息, 即衍射效应, 而LS算法或PS算法无法提取出样品的相位二阶导数信息, 因此引起了样品边界处的明暗相间的条纹状伪影; 2) 在样品信息提取算法中, 需要对角度曲线$S'\left( \psi \right)$进行余弦拟合, 而实际角度曲线在上述位置存在较强高次谐波, 并非完全可以用单一余弦函数来表述, 故拟合过程存在误差, 这同样会使得样品边界部分信息提取出现上述明暗相间的条纹状伪影^[17,24].
为了进一步研究探测器的像素尺寸与光栅自成像周期不匹配时对信息提取的准确性的影响. 在图1所示成像装置中换用多个不同像素尺寸(0.54, 0.57, 0.60, 0.63, 0.66 μm)的探测器进行样品的投影数据的获取, 此后利用PS算法以及LS算法分别提取样品的吸收、折射和散射信息, 并按照图6和图7中虚线位置处提取样品的吸收信息、折射信息和散射信息的强度分布曲线, 其结果如图8所示.
图 8 不同像素尺寸的探测器提取样品的吸收、折射和散射信息　(a)—(c) 利用PS算法提取的样品吸收、折射和散射信息的强度曲线; (d)—(f) 利用LS算法提取的样品吸收、折射和散射信息的强度曲线
Figure8. The absorption, refraction, and scattering information of the simulated sample with different pixel sizes: (a)?(c) Depict the profiles of extracted absorption, refraction and scattering images with different pixel sizes by PS algorithm; (d)?(f) depict the profiles of extracted absorption, refraction and scattering images with different pixel sizes by LS algorithm.

如图8所示, 当探测器的像素尺寸与光栅自成像的周期0.60 μm之间的偏差越大, 即在$\alpha $的绝对值越大时, PS算法提取的吸收、折射和散射像中得到的与光栅方向平行的条状伪影就越强. 而LS算法的提取结果受到的探测器的像素尺寸的影响相对较小. 考虑到两种算法都是假设摇摆曲线为余弦曲线, 并进行N点采样后进行拟合. 因此, 其信息分离准确性很大程度上取决于其余弦假设的可靠性. 由于PS算法要求拟合所用余弦函数周期必须等于实际角度曲线的周期, 在免分析光栅相位衬度成像系统中PS算法拟合角度曲线所用的余弦函数的周期是N个探测器像素宽度之和$N \cdot b$, 当$\alpha = 0$时, 拟合角度曲线所用的余弦函数的周期$N \cdot b$与实际角度曲线周期p相等, 故PS算法与LS算法都能得到精确的信息提取结果, 但$\alpha \ne 0$时, 实际角度曲线的周期p不再等于N个探测器像素宽度之和$N \cdot b$, 即拟合所用余弦函数周期与角度曲线周期不匹配, 所以会影响拟合的可靠性, 并且随着$\alpha $绝对值的增大, 余弦函数周期与角度曲线周期的差异越大, 提取得到的样品信息的失真更加严重. 而LS算法中余弦函数的周期可以是任意的, 不需要是一组探测器的像素尺寸之和$N \cdot b$, 因此随着像素大小的变化, 能够较为稳定的分离得到吸收、折射和散射像.
此外, 在上述免分析光栅一次曝光相位衬度成像系统中, LS算法与PS算法对单张投影图像的处理方法相同, 均是将探测器采集到的单张原始投影图按像素抽取为N张新投影图像, 故LS算法与PS算法所得的样品吸收、折射和散射图像信息的空间分辨率是相同的, 均为探测器采集到的单张投影图空间分辨率的${1 / N}$倍.

基于免分析光栅X射线相位衬度成像系统的一次曝光成像理论由于不需要使用分析光栅, 从而避免了分析光栅对X射线的衰减, 提升了整个成像系统中X射线的利用效率, 同时一次曝光也减小了对样品的辐射损伤, 因此基于免分析光栅X射线相位衬度成像系统的一次曝光成像理论的发展是近年来光栅成像技术发展的方向, 但是在该成像理论中传统的相位步进算法(PS算法)对光栅周期与探测器像素尺寸的匹配提出了严格的限制条件, 而空间谐波分析算法会严重降低样品吸收、折射和散射图像信息的空间分辨. 这些问题的存在限制了一次曝光成像理论的应用和发展. 针对上述问题本文提出了适用于上述成像系统的基于超定方程理论的最小二乘法的样品信息提取算法(LS算法). 该算法可以获得与传统PS算法具有相同空间分辨率的样品吸收、折射和散射信息, 此外在光栅周期与探测器像素尺寸存在匹配偏差时依然可以得到可靠的信息提取结果. 因此, LS算法的使用有助于降低上述成像系统对实验装置精度的要求, 增加成像装置的实用性, 从而使得快速、低剂量光栅相衬成像技术在临床医学、材料等领域更具实用前景.