X射线计算机断层成像技术 (X-RayCT) 作为一种重要的医学诊疗手段和新型的无损检测技术 (李保磊等，2008)，由于其独特的成像优势，在医学和工业界应用很广.其基本物理原理可用Beer定律表示：物体的三维结构信息以物体吸收系数的形式存在，光在物体中传播一小段距离Δx后，可用如下公式表示强度的损失：

其中ΔI是光强的衰减量，I是光强.如果需要表示一束光沿着线L通过物体后强度的损失，该式可以表示为

I₁为光线穿过物体后的光强，I₀为入射光强.我们需要从上述积分中获取表征物体结构信息的吸收函数f.在实际应用中，上述公式也可以看成将一个二维空间的函数映射到它的线积分上的过程，这在数学上可以表示为Radon变换，因此利用CT方法反演结构信息可以理解为一个求解Radon逆变换的过程.现有的CT重建算法可以分为解析重建算法和迭代重建算法两类 (Natterer, 2001)，在实际应用中，通常采用滤波反投影 (FBP)、代数迭代 (ART) 方法.滤波反投影算法是一种显式的算法，但往往会出现环状伪影的现象; 代数迭代算法在结果上可以更为准确，但计算时间往往要更长.传统CT的数据通常为扇状投影数据，但对于利用同步辐射光进行的岩石CT投影重建的实际情况而言，投影数据应为单色的平行光扫描.由于同步辐射光具有高能量、高准直、精确的能量调节等优点，传统的扇状数据的投影算法并不适用，而可以利用平行光投影的处理方法进行数据处理.

传统的CT成像方法在处理小分辨率下的岩石样本时会面临诸如成像结果信噪比低，弱吸收物体难成像的困难，当物体与检测器拉开一段距离时，探测器采集得到的投影图像是样品吸收像和相位像的混合像.实验可以观测到，光强分布随着检测器与物体距离的增大而发生变化 (Bronnikov, 2006)，这个现象证明不能简单地利用Radon变换来模拟CT成像的全部过程.因此就一定要引入物体的相位信息进行结构成像.在硬X射线 (10~100 keV) 下，样品的折射率吸收系数约为相位移动系数的1000~2000倍左右 (Gureyev et al., 2013)，其中投影图像中相位衬度像会占主导优势，具体表现为边缘增强效应 (陈荣昌等，2009)，这会极大地影响最终的成像，因此我们需要引入相位校正.

利用物体相位信息的X相位衬度成像在越来越多领域受到了关注，这种方法已经在生物软组织成像 (Bravin et al., 2013)、材料科学 (Mayo et al., 2012)等领域取得了很大进展.相位衬度方法是将物体结构分布看成是三维的复折射率，n (r)=1-p (r)+iβ (r)，这里p代表相位因子，β代表吸收因子，r是物体的三维坐标.β与线性吸收系数μ有关，μ=4π/λ·β.根据不同光学装置, 相位衬度有不同的实现方法，如：同轴相位衬度成像 (in-line phase contrast)(Lee, 2015)，基于光栅干涉的相位衬度成像 (grating interferometer)(Nesterets and Wilkins, 2008)，基于检偏器的相位衬度成像 (analyzer-based) (Bravin, 2003)，以及基于编码孔径的相位衬度成像 (coded aperture)(Olivo et al., 2012)等.其中同轴相位衬度成像的光学装置相对简单，对光学装置的摆放位置不是非常敏感，且光学能量没有流失，作为代价这种方法的数据处理可能较为复杂 (Lee, 2015).

一般利用相位信息对物体进行结构反演时存在两种最基本的假设：即弱吸收 (μ≈0)(Bronnikov, 1999)和吸收相位耦合(μ∝β) (Wu and Liu, 2005; Paganin et al., 2002; Beltran et al., 2010)，对于研究物体为岩石的情况，物体对光的吸收显然是不可避免的，因此应采取第二种假设.相位信息对于岩石成像的作用是一种“校正”，通过选择合适的相位移动吸收比，可以较好地提高信噪比，同时减少相衬成像的边缘增强效应.通常的校正过程采用的是频域的算法 (Bronnikov, 1999; Wu and Liu, 2005; Paganin et al., 2002; Beltran et al., 2010).在国内刘慧强等(2012) 利用同步辐射光，基于修正的Bronnikov算法对有机复合样品的成像进行研究，叶琳琳等(2013) 利用相位吸收二重性Paganin算法对植物器官的显微层析实验数据实现相位重建和优化处理.在频域上求解TIE方程是因为，直接利用TIE进行相位恢复其结果将受到低频噪声的严重干扰 (Paganin et al., 2004)，频域求解的过程相当于利用正则化的方法对这些频率的噪声进行压制.Sixou (2015) 对此问题在频域求解进行过正则化的理论证明和分析.

本文通过构建算子方程，提出了一种在空间域解决相位恢复的算法.该算法可以把问题表示为u_e=f+e，其中是含二阶差分项的算子，u_e是观测数据，f是待恢复的只含吸收信息的投影数据，e是观测噪声.该模型在数学上可以转换成一个变分问题进行表达，并通过正规化方程求出问题的解.本文在数值上探讨相位恢复问题的不适定性并提出了一种基于空间域的正则化相位恢复的计算方法.

本文第2节介绍TIE方程的基本理论并导出我们的空间域正则化算法，第3节是数值实验以及实验的结果分析，最后给出讨论和总结以及对未来工作的思考.

X射线同轴相位衬度层析方法具有装置简单、成像直观、分辨率高等优点，其原理如图 1所示.首先约定系统笛卡儿坐标为 (x₁, x₂, x₃)，检测器的坐标为 (x, y)，其中y与x₃方向相同，几何构造与坐标约定如图 1所示.

图 1

Fig. 1

图 1 待测物体坐标系和观测平面 Fig. 1 Coordinate system for the measured object and the plane of observation

同时约定入射光波场U_i为平面波 (考虑到同步辐射光的高准直性，这个假设是可以实现的)，经过物体后的波场U_θ可以表示为

(1)

其中，T_θ (x, y) 为物体对光的透射函数，它包含了对光强度的吸收和相位的延迟，可以表示为

其中，吸收和相位延迟满足Beer定律，有如下积分公式：

约定μ是物体的线性吸收系数，g (x₁, x₂, x₃)=real[n(x₁, x₂, x₃)]-1，其中n是物质的复折射率，δ(·) 为狄拉克函数，上述两个式子描述了二维的Radon变换.在距离物体为d的平面可以将光的强度描述为 (Bronnikov, 2002):

(2)

式中^{* *}代表二维卷积，h_d代表按如下定义的菲涅尔传播子：

如果考虑近场条件λd＜＜D² (D为物体尺度，λ为光波长)，该式可退化为

(3)

其中▽²为拉普拉斯算子，为方便计算，可以将其表示为▽²U_θ=(a_θ+ib_θ).其中

如果设I_θ⁰为在0平面的光强分布，可以对 (3) 式做如下简化：

上述约等号成立的条件是近场条件，即λd＜＜D².若进一步假设物体的μ_θ变化不大，可以得到在近场时候相位与光强的关系式：

(4)

由 (4) 式可知，投影图像由吸收信息I_θ⁰和相位信息φ_θ的拉普拉斯变换组成，因此在图像上表现为边缘增强效应.如果对样品做单一性假设，令T (r) 为样品的投影厚度，可以得到表征投影厚度的TIE方程 (Paganin et al., 2002)：

(5)

(5) 式等号右侧为观测值，左侧表征相位的影响.通常的做法是对 (5) 式等号左右两边做傅里叶变换，可以得到 (Paganin et al., 2002)：

其中，k_⊥指的是平面频率域坐标，“⊥”表示与光线传播轴相垂直的平面.这相当于在频域做一次滤波，反演得到投影厚度T (r) 进而用滤波反投影的方法 (FBP) 求出模型参数.本文没有进行傅里叶变换，利用传统的频率域方法求解TIE方程.我们提出直接在空间域将 (5) 式等号左侧的运算写成算子运算的形式，先求出e^{-μT (r)}，再用传统的滤波反投影方法求出样品的投影厚度.因而，上述公式可以简写为

(6)

其中

在实际问题中，由于观测到的数据通常是带噪声的，因此 (6) 式应修改为

(7)

因此我们无法保证u_e位于算子的值域中，这就使得 (6) 式如果直接用f= ^-1_ue求解可能是不相容的，因此需要在某种度量空间下求解下述残差的极小值：

(8)

由于中含有微分算子，且微分算子的谱是无界的，用一般的差分法或有限元法所导出的微分方程的离散近似往往是病态的，而且随着离散尺寸的缩小将更加严重 (肖庭延等，2003)，由此可见相位恢复的反演问题是不适定的.求解不适定性问题的策略是用一簇相邻近的适定问题的解去逼近原问题的解，如果把原问题写为

(9)

其中是由度量空间 (F, ρ_F) 到度量空间 (U, ρ_U) 的连续算子，逆算子^-1存在、单值但未必连续，这时如果有算子R (u, α): U→F是对U中所有元素u和任意α > 0都有定义的关于u的连续算子，且f∈F有，则称算子R (u, α) 是问题 (9) 的正则算子.当这簇适定问题满足正则化条件的时候称这种解决不适定性问题的方法为正则化方法. Tikhonov通过最小化如下格式的展平泛函 (smoothing functional) 来构造正则算子 (肖庭延等，2003)：

其中Ω[f]为定义在F稠子集内的非负连续泛函，也称作稳定泛函，通常取Ω[f]=ρ_F² (f, f₀)，因此问题转化为令M^α[f, u]对f取极小，即求解一个变分问题，可以得到一组Euler方程.本文中，我们把TIE方程离散化，然后再利用上述的正则化技巧.

我们希望将算子方程离散成如下形式：

(10)

其中A是算子的离散化表达，对于中的二阶差分算子，常用的离散形式可以通过附近两点构造离散算子，但这种差分格式在精度上不高，因此考虑用以下格式进行替换：

其中，Δx=Δy为差分离散化步长.对于差分系数，可根据一维函数的泰勒展开并解相应的方程组得到：a₁=-1/12；a₂=4/3；a₃=-5/2；a₄=4/3；a₅=-1/12.

如果对于边界点的差分依旧遵守根据泰勒展开的公式推导，则会推出一个不对称的差分矩阵，而根据这个矩阵得到的算子的条件数将非常大cond (A)=O (10⁹)，如果对于投影数据边界做稳定值补值，则可以不影响结果的对拉普拉斯算子做对称化处理：根据需要补充算子边界部分使其满足五对角对称阵的形式，这样可使得矩阵的条件数降至10⁶量级.对称化后的矩阵有如下形式：

采用这种差分格式，可以构造算子A以满足较高精确度.构建了离散算子A之后问题转化为对线性方程组 (10) 的求解.

随着实际问题所要求的分辨率的提高，差分所导致的算子A条件数会越来越大，如果我们在L₂空间求解残差的极小值，便需要求解一个正规化方程组A^TAf=A^Tu_e，但正规化方程将是一个更为严重的病态问题，经统计它的谱分布不集中.我们称直接求解上述正规化方程组的方法为直接法.

对方程组 (10) 实施正则化过程，即求解下述的极小化问题：

其中α为大于0的正则参数，该问题等价于求解线性方程组：

这里α > 0称作正则参数, 并称式 (7) 中的e为误差水平.利用Tikhonov正则化方法，可以用正则解f_e^α=(αI+A^*A)^-1A^*u_e作为精确解f_T=A⁺u的近似 (王彦飞，2007).正则参数可以看成是误差水平的一个函数，并用α(e) 表示，如果正则参数具有某种先验性质，如 则=0. Groetsch (1983) 证明，对于任意的光滑性条件A⁺u∈R ((A^*A)^ν), 其中0＜ν≤1，和任意的α(e)，正则化方法所能达到收敛速度不会优于但若采用下述的迭代Tikhonov正则化方法，就可以得到更高阶的收敛速度 (Gfrerer, 1987).迭代Tikhonov正则化方法可以按如下定义：

在本文计算中采用了下面的方式 (肖庭延等, 2003)：

如果利用迭代Tikhonov正则化方法来求解问题 (10)，则所能达到的最佳收敛速度为 (Nashed, 1976):, 其中0＜ν≤n，此时对正则因子的选取为α(e)= c > 0, 且最优收敛速度在ν=n达到.

为了检验新算法的正确性，同时也能更好地模拟实际数据，本文的模型参数按照页岩的基本成分的吸收系数选取 (王玉丹等，2015).

模型选择以石英为背景，模型设计的情况及各成分所占百分比如表 1和图 2所示.

表 1(Table 1)

表 1 模型成分在30 keV能量下β的理论值 Table 1 Theoretical β values for materials in the model at energy E=30 keV 序号 成分 百分比 β 1 石英 75.04% 7.4×10-10 2 高岭石 3.20% 6.5×10-10 3 伊蒙间层 3.20% 7.2×10-10 4 斜长石 3.20% 9.0×10-10 5 白云石 2.92% 11.5×10-10 6 伊利石 3.02% 12.4×10-10 7 方解石 3.10% 16.0×10-10 8 绿泥石 3.04% 20.5×10-10 9 孔隙 3.27% 0

表 1 模型成分在30 keV能量下β的理论值 Table 1 Theoretical β values for materials in the model at energy E=30 keV

图 2

Fig. 2

图 2 线性吸收系数模型 Fig. 2 Model of linear absorption coefficients

在模型的构建上，根据实际纳米-CT实验经验，选取每个像素的大小为50 nm，分辨率为512×512，模型大小为25 μm，角度分辨率为1°(扫描180次) 或2°(扫描90次)，选定的光强为30 keV，在此光强下波长约为40 pm (Zschornack, 2007)，选定距离为10 cm，相位移动与吸收比例设定为1200.由于理论的算子为对应于物平面的拉普拉斯算子，为了检验算法的有效性，对模型进行了简化，上述模型我们均认为在y方向上各参数不发生变化.为了检验算法的准确性，记data为模型数据，data_res为重建数据，我们按照如下定义相对误差：

正演过程首先计算光强的衰减，模拟对模型进行如图 3所示的扫描.因为理论模型虽然在构建的时候是由圆形或椭圆形组成，但在分辨率内每个格点吸收系数恒定，因此在计算投影值时采用了逐点扫描而非传统的解析解.利用这种算法虽然在计算量上增大，但可以得到更为精确的投影值以检验算法的准确性.得到逐点扫描的数据后，还要根据公式 (2) 做卷积运算来模拟计算理论投影数据，得到如图 4的结果.

图 3

Fig. 3

图 3 观测系统示意图 Fig. 3 Schematic diagram of observational system

图 4

Fig. 4

图 4 模型正演的投影数据 Fig. 4 Projection data of forward modeling of model

如果对于未经处理的投影数据直接利用滤波反投影方法求解，残差为62.7，反演图像上只会得到如图 5所示的模糊的图像边界信息.图 6a、6b和6c分别为直解法，频域法和迭代正则化方法在无噪声情况下进行的相位恢复并经过滤波反投影之后算出的结果.其中三种方法的相对误差Err分别为0.016, 4.46, 0.016.在无噪声的数值实验中，迭代Tikhonov正则化算法的正则参数取为10^-6.

图 5

Fig. 5

图 5 未经处理直解反投影 Fig. 5 Anti-projection of direct solution without phase retrieval

图 6

Fig. 6

图 6 线性吸收系数反演结果 Fig. 6 Inversion results of linear absorption coefficients

为了更好地模拟实际问题，接下来在正演过程中加入了噪声，图 7a、7b、7c分别为直解法，频域法和迭代正则化方法在0.5%高斯白噪声情况下进行的相位恢复并经过滤波反投影之后算出的结果，此时的相对误差分别为2.25, 4.62和0.27.在带噪声的数值实验条件下，迭代Tikhonov正则化算法的正则参数取为0.0041.

图 7

Fig. 7

图 7 噪声1%下线性吸收系数反演结果 Fig. 7 Inversion results of linear absorption coefficients with noise 1%

尽管视觉上，图 7b和7c的计算结果类似，但图 7b的值偏离真解太远.对于图 7b和7c右下角框有局部图 8a和8b.本文所进行的数值实验是在LENOVO生产的x64-based PC独立工作站上运行的，处理器为Intel (R) Xeon (R) CPU E5-2609.

图 8

Fig. 8

图 8 图 7b, 7c中方框内的细节图 Fig. 8 Details in the square of Figs. 7b and 7c

由上述反演结果可以看出，直解法的不稳定性体现在有噪声干扰的情况下，在0.01噪声尺度下几乎看不清内部的信息.在同样的条件下，处理毫米甚至微米级样本的实验中信号的信噪比不会很高，因此虽然空间域的直解法在计算速度上很快，但很难对真实数据进行成像.不过在残差的数值结果上与频域法比较，空间域的直解法保持了较大优势，而且对于模型线性吸收系数的反演结果，在数量级上也没有偏差.

频域法虽然在一定噪声尺度下依旧能够呈现信息，但该方法只是在形态、大小关系上与模型相似，只是可以刻画内部结构的边界而无法算出具有物理意义的数值，该方法在解的最终结果上与模型参数存在数量级的差异 (具体可参见图右侧色标数值).而空间域方法虽然对比度不如频率域方法，但所求值与模型真实值更为接近，可以计算出，空间域方法在残差方面要远小于频率域方法.而且如果观察原始模型孔隙位置 (即图 7b画圈位置) 可见同样的模型参数，频域法的求解结果在灰度呈现上相差很大，而且在模型参数较高的区域 (见图 8a和8b)，频域法结果在周边存在较严重的“晕”和一些线状伪影，尽管由于滤波反投影算法存在从中心向四周发散的伪影，但频域法还存在一些指向较高数值区域的线状伪影和一些其他的横竖状伪影，可见频域法的结果可能会受周边模型参数值的影响，最主要的是在残差上要大于基于空间域的直解法和Tikhonov正则化迭代解法.

本文所提出的方法，在没有噪声的时候能够和直解法、频域法一样成清晰的像；在有噪声的时候也能体现出对于噪声较好的抵御能力，同时在保证残差足够小的情况下在成像上依旧比传统的频域法好.在更换几组实验条件下 (如调整分辨率，调整样品检测器距离，调整光强等) 均可以呈现较好的图像，可见在各种情况下空间域方法表现优秀.

通过上述实验分析，本文算法较稳定地获得了页岩微纳米孔隙结构成像，对噪声表现出了较好的压制效果.理论上尚存在的问题是，在上述公式推导过程中我们假设了模型的相位移动和吸收比例是恒定的 (Paganin et al., 2002)；而在实际问题中，可能需要根据先验的物质吸收边信息合理选取多种能量，尝试不同的相移吸收比以更好地提高图像重构的信噪比 (王玉丹等，2015)，得出准确的图像.这些研究我们将在接下来针对实际数据的问题中加以解决.

相位恢复的主要作用是消除在物体吸收系数变化很大的边界可能存在的“边缘增强”效应，上述理论模型解决了二维样品的相位恢复，没有考虑物体吸收系数和相位移动因子在与切片面垂直的方向，即图 1中的y方向发生变化时的情况，因此在处理实际问题时需要交叉的数据对投影强度数据进行相位恢复，即算子A在两个方向进行差分.

对于参数数值相差不大的区域，可以通过FBP算法中滤波器的选择以增强统计图的峰值.一般汉明窗 (Hamming) 的效果要优于Ram-Lak滤波器.但即使这样，直接求解逆Radon变换，仍然存在伪影和“振铃”现象，这是FBP作为一种显式计算方法快速计算所很难解决的问题.因此为了更清晰地成像，在反投影过程以后的研究中可以考虑使用迭代优化的算法进行图像重构.

感谢审稿人耐心细致的评审意见，使本文内容得到了很大的充实.