0 引言

随着能源需求的增加，页岩油气作为一种典型的非常规油气资源得到了越来越多的关注.我国在四川盆地古生界海相页岩和陆相页岩的页岩气勘探开发中均取得了重大突破(崔景伟等，2012).页岩气主要赋存于页岩基质的微孔、夹层、裂缝中，以吸附式、活塞式和置换式聚集(滕吉文和刘有山，2013).不同于传统能源“圈闭成藏”的特征，页岩油气致密储集层具有纳米级孔喉网络系统，兼具生、储、盖等多项功能(邹才能等，2012)，且其中的纳米级孔喉储层中的流体多为非达西渗流，传统“渗透率”等参数已不足以描述其渗透能力.因此页岩微纳孔隙的连通性以及定量识别，对于孔隙分类、孔隙演化规律的认识甚至赋集油气有效性的判别均有重要意义.孔隙特征研究一直受到广泛关注，油气储层孔隙表征的常规方法有高分辨率的场发射扫描电镜(Sondergeld et al., 2010)、原子力显微镜(Javadpour，2009)、氩离子抛光(王羽等, 2015, 2017; Wang et al., 2018)、聚焦离子束系统(Zhang and Klimentidis, 2011)以及不同于光学辐射的压汞分析(田华等，2012)等流体渗入方法.

X射线计算机断层成像技术(X-ray CT)作为一种重要的医学诊疗手段和新型的无损检测技术，由于其独特的成像优势，在医学和工业界应用很广.高能同步辐射光的引入极大地促进了基于CT的，纳米级尺度下三维孔隙结构的成像研究.在我国非常规油气储层研究的历程中，油气纳米孔的首次发现，也是利用了纳米CT的重构技术(邹才能等，2011).如何提高页岩纳米CT成像的实验精度，进而明确其适用范围是页岩纳米孔隙表征问题研究中的热点和难点之一.已有人利用微纳米CT技术对页岩样本，尤其是重庆漆辽龙马溪组页岩样本进行成像，积累了大量经验，并与其他方法进行了比对(王羽等, 2016；Wang et al., 2017a；Wang et al., 2018)，其中大部分的CT重构部分采用的是基于软件的滤波反投影算法；王彦飞(Wang, 2017)等对于页岩的微米CT问题提出了稀疏正则化方法并利用梯度下降方法，使得图像分辨率更高且噪声更少，有效的提高了孔隙的分割精度.

CT重构的基本物理原理可用Beer定律表示：物体结构信息以物体对光的吸收系数的形式表征，光在物体中传播一小段距离Δx后，可用如下公式表示强度的损失(Natter, 2001)：ΔI/I=f(x)Δx，其中ΔI是光强的衰减量，I代表光强.如果用一束光沿着一条线通过物体后，强度可以表示为I₁/I₀= ，其中I₁是光线穿过物体后的光强.我们需要从上述积分中获取表征物体结构信息的吸收函数f.在实际求解过程中，上述公式在数学上可以理解为Radon变换，因此利用CT方法反演结构信息可以理解为一个求解Radon逆变换的过程.在纳米尺度CT成像方面，可以采用除了基于吸收还有包括基于相位的成像方法(唐巍和王彦飞, 2017)，两者各有优缺点可以互为补充，但这两种方法都会用到求解Radon逆变换的步骤，因此准确高效的做好数值Radon逆变换对于解决CT纳米孔隙成像问题有着重要意义.

传统的CT重建方法是以目标为中心，扫描不同角度下的投射数据，利用重建算法还原出目标的断层图像.1985年，Tuy和Smith给出了精确重建的充分必要条件(Tuy 1983；Smith, 1985).可以证明，由完全投影数据可以重建唯一的图像，但当投影数据不满足Tuy-Smith完备性条件时，现有的经典滤波反投影重建算法无法精确重建断层图像(杨富强等，2014).此外，我们还可以使用基于统计模型的方法，如Bayes方法(Herman and Lent, 1976)，最大熵方法(Skilling and Bryan, 1984)等；迭代重建的方法，如ART算法(Gordon et al., 1970)，SIRT算法(Hudson and Larkin, 1994)等，由于迭代重建方法在有噪声数据和有限角度数据情况下的优秀表现，在越来越多领域得到了人们的重视.

Kaczmarz方法是一种代数迭代算法，最早由波兰数学家Kaczmarz提出(Kaczmarz, 1937).Gordon等将Kaczmarz方法应用到医学成像领域(Gordon et al., 1970).最早的CT扫描装置利用的就是代数迭代重建技术来重建图像，但由于当时的计算能力有限，成像速度较慢，后来逐渐被滤波反投影方法取代.传统Kaczmarz方法是一种求解线性相容方程组的方法，它按照矩阵行的顺序，依次把初始值投影到由矩阵每一个行向量和所对应观测值共同决定的超平面上, 这种方法可以被认为是一种凸集投影方法的特例.很明显，Kaczmarz的收敛性质依赖于行顺序，如果行顺序不好的时候，Kaczmarz的收敛结果非常慢，而且很难对于该方法做出收敛结果的分析.Herman等(Herman and Meyer, 1993)已在数值实验上观测到当随机的对矩阵的行进行选择的时候，可以得到收敛性质的显著提高.Strohmer和Vershynin在2009年提出了一种解决满秩矩阵相容问题的随机Kaczmarz方法(Strohmer and Vershynin, 2009)，并且给出了这种方法在期望上满足线性收敛性质的证明.Needell给出了在不相容问题上，这种随机Kaczmarz方法可以收敛到最小二乘解的一个邻域内的证明(Needell, 2010)，且证明该球的半径与不相容问题的噪声大小有关.Eldar和Needell (Eldar and Needell, 2011)依据Johnson-Lindenstrauss准则给出了一个更优秀的行选择办法，使得收敛率可以进一步提高，但是这种方法需要大量额外的计算量，而且并没有解决随机Kaczmarz方法在不相容问题上不收敛到最小二乘解的问题.Zouzias和Freris在Popa方法的启发下提出了一个加强RK (Rondomized Extended Kaczmarz)方法(Zouzias and Freris, 2013)，主要思想为将观测数据投影到由观测矩阵列向量张成的空间中，这使得不相容问题可以收敛到最小二乘解，但迭代需要较长时间，计算量很大.本文所需要求解的矩阵元素代表像素对于投影射线的贡献，因此是一个非常稀疏的大型矩阵，不适合求逆等运算.同样作为代数迭代方法，SIRT在每次迭代过程中使用全部的矩阵数据，这会导致收敛速度较慢的情形.而且对于CT的迭代成像而言，算子的条件数很大，使得问题成为一个不适定问题.正则化方法(王彦飞, 2007)是解决不适定问题的一种有效解决方法，它通过最小化如下的目标函数得到问题的解：，其中α是正则参数，用来平衡实际问题中近似解的精确性和稳定性.矩阵A是高维度的，且非常稀疏，这使我们只能通过迭代法对上述最小化问题进行求解，一般的正则化方法都通过求解如下的正规化方程(A^TA+αI_n)x=A^Tb进行求解，其中I_n为一个单位算子，但是当矩阵算子A的条件数很大时，会导致解的不稳定性.在本文中我们把上述求解系统转换为增广正则(Ivanov and Zhdanov, 2013)的表达形式，并利用Kaczmarz投影法对于增广表达的线性求解系统进行迭代求解，又发现这种方法对于目标解的初始值有较高的要求，因此使用FBP方法计算初值.

我们需要通过CT方法无损的获取页岩样本内部大量的微米-纳米孔隙，这些孔隙往往是页岩气的重要储集空间和渗流通道，选用好的方法对岩石内部结构进行精准成像，对于页岩气勘探层位置选取，以及后续的勘探开发有着非常重要的意义.页岩的纳米CT重构在不加先验信息的时候是一个欠定问题，在迭代的过程中我们对页岩样本的吸收系数在整体区域上进行离散化处理，迭代每次算出的结果(Wang et al., 2017b).但是实际实验过程中，每次迭代在样本周围没有物质的地方也会算出结果，这时需要加一些约束项使得周边的吸收系数为0，这样会进一步提高由该算法所能计算出的效果，在迭代的具体过程中引进了如前所述的关于随机行选择的策略使得迭代的速度得到了明显的提高，基于这些改进，本文提出了一种带约束的块状Kaczmarz正则算法.

本文分为以下几个部分，第1节主要介绍基本的Kaczmarz方法，第2节介绍本文根据问题的不适定性提出的处理方法，第3节通过两个数值模型：经典Shepp-Logan模型和页岩模型验证算法的收敛性质，第4节对实际数据进行处理，最后一部分给出讨论和总结以及对未来工作的思考.

1 Kaczmarz算法 1.1 经典Kaczmarz算法

Kaczmarz方法旨在解决如下形式的线性方程组Ax=b的求解.其中，.本文中用a_i^T或a_i表示A的行向量，用a_(j)或a_(j)^T表示A的列向量，用b_i表示b中的第i个元素，即

经典的Kaczmarz算法采用下面的循环迭代的方式

其中，ω是松弛因子，一般取在(0, 2)之间.

令F是从到 的映射，并按照如下方式定义复合函数：

则Kaczmarz算法可以进一步定义为：选择任意的元素，按照如下循环得到迭代序列：

上述算法有很明显的几何解释，如图 1所示：f_i是把 中的一个点，投影到如下定义的超平面〈x,a_i〉=b_i.

对于CT数据而言，x是未知图像的列向量表示，b是投影数据的列向量表示，A是m×n的投影系数矩阵，每个元素a_{i, j}代表第j个像素对第i条投影射线的贡献.计算a_{i, j}的模型有好几种，常用的是以下三种：

1.面积权值：假定射线具有一定宽度，a_{i, j}设为射线与像素相交的面积；

2.距离权值：假设射线没有宽度，以像素中心的坐标与射线的距离定义权值.

3.长度权值：假设射线没有宽度，a_{i, j}设为射线与像素的交叉长度.

因此完整的Kaczmarz算法可以写成：

算法1   (Kaczmarz算法)

Step 1：输入投影数据b，观测矩阵A, 其中A矩阵的大小为m×n；

Step 2：设置初始化最大的迭代次数为K_max，令k=0，x_k=0；

Step 3：取出一个指标i=mod(k, m)，进行投影

Step 4：当k≤K_max时，令k=k+1，返回Step3;否则，结束循环，输出当前结果.

算法1有至少以下两个缺点：第一，由于是依照行顺序来选择，如果行顺序很糟糕，可能会收敛速度非常慢；第二，很难对收敛率进行理论分析.

1.2 分块Kaczmarz迭代算法

上述方法在每次迭代过程中只用到了矩阵一行的信息，在行顺序不好的时候可能会收敛的较慢，在运算上也更耗时，因此有人提出了块状的Kaczmarz投影算法(Elfving, 1980)，每一次选取矩阵A中的一些行，在实现过程中可以将矩阵A和投影值b分成几个独立的区块，如

在每次迭代过程中使用其中一个区块的数据，设置内外两层循环，外层循环指外部大循环，内层循环每次利用多行矩阵和投影数据，其中A_l为A中独立区块的一部分，则迭代过程可以表示为

其中，A_l^†是A_l的广义逆，K为外层循环的次数，l=1, 2, …, m/p，diag(·)表示对角矩阵.这种方法的好处是可以节省运算的时间，可以较快的收敛到目标解，当p等于1的时候与上述经典Kaczmarz方法相同.针对本文所要求解的问题，矩阵A是非常稀疏的，且在图像分辨率达到512×512的时候，矩阵的列数较多，在实际使用过程中用A_l^T[diag(A_lA_l^T)]^-1(b^l-A_lx_{K, l})/p替代了广义逆.

2 不适定问题的求解 2.1 问题的不适定性

需要指出，上述分析和算法仅适用于相容问题，对于不相容问题：

其中 是噪声，投影算法并不一定能收敛到最小二乘解，这一点在几何上可以如此理解：当存在噪声的时候，x_k-x与a_i不再正交，因为a_i(x_k-x)=b_i-a_ix=r_i≠0.相当于Kaczmarz迭代的每一步投影到了错误的平面.本问题的矩阵A是非常稀疏的，如果对矩阵A进行奇异值分解有：

其中，U,V为正交矩阵，Σ=diag(σ₁, …, σ_{min(m, n)})，不相容问题的最小二乘解可以表示成

对于非常小的特征值σ_i，求解时，数据b中的噪声将被放大，会出现结果非常不稳定的情况.实际的纳米数据512×512分辨率，投影仅有180个角度，因此离散化后得到的方程组是一个欠定的方程组.同时，CT的重构问题，矩阵A表示的是光线在物体内部传播过程中发生变化的物理模型，其中每个元素a_{i, j}代表第j个像素对第i条投影射线的贡献，因此对于矩阵 每一行最多有 个非零元素，这种稀疏性对于求解的稳定性是很大的考验.可见，无论从唯一性还是稳定性上均不满足适定性条件.

代数迭代算法有可能会出现半收敛的现象，即在有噪声的情况下，计算误差随着迭代次数的增加而先减少后逐渐增加，即便迭代次数足够多，也不一定能收敛到可靠的结果(王彦飞，2007).事实上，无论对于Kaczmarz投影方法，还是块状投影方法，均可观察到半收敛(Elfving et al., 2014)现象，因此需要正则化手段处理.

2.2 增广正则矩阵

Tikhonov提出可以通过求解如下的约束最小二乘问题得到问题的逼近解(王彦飞，2007)：

(1)

在实际计算过程中，极小化问题(1)的求解可以表示成以下增广正则矩阵-向量方程组(Augmented regularized normal equations)的形式

(2)

定义分别为对应(2)的增广正则矩阵、新的未知向量和右端项，方程组(2)可以简写为

(3)

其中，I_m∈R^m×m、I_n∈R^n×n分别是单位矩阵.可以证明，对于任意的α＞0，矩阵 是非奇异的，而且解满足 =(y_*^T,x_*^T)^T，其中，x_*=(A^TA+αI_n)^-1A^Tb，y_*=ω^-1r_*，r_*=b-Ax_*.

方程组(3)左边矩阵的谱条件数是

其中σ_max和σ_min分别是矩阵A的最大最小奇异值.对于增广矩阵(2)而言，谱范数会小一些(Ivanov and Zhdanov, 2013)

上述增广矩阵可以写成两行的形式，我们需要同时满足这两个线性代数方程组：

(4)

(5)

对于方程组(5)式可采用Kaczmarz的方式进行迭代求解，为表达方便，令A^ω= A^T ωI_n，令A_i^ω表示A^ω的第i行；在每一轮的迭代过程中的第i步可以选取(5)式中矩阵A^ω的第[i]行进行投影，其中，[i]∈(1, …, n)，即：

(6)

其中，(e₍1), …,e₍n))=I_n.现在我们定义

(7)

对于迭代格式(6)，可以根据，经过适当变换，继续分解为两步：

(8)

(9)

其中

(10)

对于增广矩阵，在满足一定初始条件的时，只需对方程组(5)进行Kaczmarz迭代即可始终保持方程组(4)的正确性.可以证明(Ivanov and Zhdanov, 2013)，当初始值满足

(11)

时，对于任意初值x₁及迭代的任意步均有

(12)

其中，r_i+1和x_i+1是根据迭代公式(8)和(9)计算出来的.

上述定理同时保证了在选择好的初值后仅使用方程组(5)构造迭代格式，则问题就可以收敛为正则解.

2.3 带约束的块状随机Kaczmarz的正则化方法

对于纳米CT成像问题，所能得到的观测角度只有180组，我们希望得到512×512的分辨率，这时的问题是一个欠定问题，理论的广义解是最小模数解，但很可能与精确解相差较大，除了在迭代过程中可以通过一些正则技巧改善解的稳定性，根据样本吸收值对解进行一定约束也可以改善问题的不适定性.以实际问题为例，样本吸收值的大致范围可以对问题提供一个有效约束.同时我们注意到虽然解的空间是n维，n等于一个方向像素个数的平方，但其实这时的目标区域并不全是样本，可以利用先验信息得到样本的基本轮廓，对于轮廓外部的像素点强制附零即

同时，在迭代过程中一个物理上的约束是：解的物理意义是样本对于光的吸收，不可能出现小于0的情况即

可以把这些约束总结成 中的约束集C，可以把任意的x投影到这个约束集里：

约束集可以在每一轮或每几轮迭代过程中约束求出的解，可以更好的聚焦目标区域.此外，在实际计算过程中可以观察到上述方法对于迭代结果的初值要求比较高，在零初值的情况下，会有局部有一些“亮斑”式伪影.因此使用较快速的FBP算法计算出的结果作为该方法的初值.

此外，在计算过程中Kaczmarz可能存在的问题是收敛速度会受到行顺序的影响，因此在行顺序选择上可以采用一些随机的策略.以的概率选择第i行进行迭代，其中‖· ‖_F为Frobenius范数，定义为

数学上已有证明，令x为上述问题的最小二乘解，上述算法可以在期望意义下收敛，对于随机的初始值在第k步的平均误差有(Needell, 2010)：

其中κ(A)=‖A‖_F‖A^-1‖₂为归一化条件数.行选择的一个最优决策方法为(Needell, 2010)

其中f_i为选择第i层进行迭代.但是如果想要实现上述算法，相当于在算法运行过程中每次采样都要被替换为进行m次投影计算，这对于实际的CT问题显然运算量过大.对于本问题而言，采用随机策略除了可以提高收敛速度外还可以有效的避免一些由于顺序选择所带来的迭代初期的较大区域性伪影，这一点在后面的数值实验中有所体现.根据随机的思想和约束正则化的思想，本文提出的带约束的Kaczmarz正则化方法如下：

算法2  (带约束的块状Kaczmarz正则算法)

Step 1：输入投影数据b，观测矩阵A, 其中A矩阵的大小为m×n；

Step 2：利用FBP方法求出初步结果x^fbp，给定正则参数ω＞0，定义最大的迭代次数为K₂，k=0，x^k=x^fbp；

Step 3：根据先验信息或通过阈值的方法得出对结果的约束集C；

Step 4：x⁰代入(7)式求出r₀，并根据(10)式求出ρ₀；

Step 5：方程组(5)中矩阵以 的概率分布对所有的行向量采样，每次取出部分指标组合成指标集，利用块状投影方法，结合公式(8)和(9)进行迭代，每一次块状迭代后，把当前结果x^k投影到约束集x^k=P_C{x^k}中；

Step 6：如果k等于K₂，输出结果；否则返回Step 4.

就随机选择的Kaczmarz算法而言，具体块大小有：

在采用分块随机Kaczmarz算法的计算过程中，算法的数值稳定性与收敛性质依赖于分块的方法.分块问题中如何选择子矩阵大小与子矩阵的条件数有关(Needell and Tropp, 2014)，具体的，对于任意子矩阵A_t有参数：

其中λ_min，λ_max分别代表矩阵的最小最大的特征值，令x为上述问题的最小二乘解, 随机选择原则下平均收敛速率可表示为其中σ_min²(A)为矩阵A最小奇异值的平方，p为每个子矩阵的行数.需要指出，现有理论仅能提供基于平均收敛率的结果.本问题中矩阵的最优分割意味着需要计算所有任意分割情况下子矩阵的特征值，计算量太大，随机取部分特定大小的子矩阵观察后可发现，随着子矩阵规模的递增，β有递增的趋势，但m/p必然有递减的趋势，经过一定的数值模拟发现，选取5~10行为子矩阵的规模为宜.

3 数值实验

Shepp-Logan模型是CT重构方法所研究的经典模型(Wang et al., 2017本条文献指代信息不明确)；页岩模型是针对实际问题提出的新的数值模型.在处理结果的呈现方法上，本节与第4节的结果均是以图像灰度表示不同像素点位置的物质对光线的吸收系数，进而在图像上表示页岩不同组分所构成的内部结构.两组模型分别进行没有噪声情况下和有噪声情况下的数值实验，为了比较算法之间的性质，成像质量用归一化均方差(NMSE, normalized mean squared error) (简称均方差)描述.误差越小，成像质量越好，NMSE定义为

主要比较本文所提出的算法2与传统的FBP方法.在实际CT实验中所采集的投影图像受噪声的影响，除去成像系统中的坏点，其噪声主要表现为高斯噪声，数值实验过程中加入一定能量的高斯噪声.

3.1 Shepp-Logan数值模型

Shepp-logan模型采用128×128的分辨率，模型如图 2所示.

当数据无噪声时，不同方法重构结果如图 3所示.

图 4给出了无噪声情况下归一化均方差随每百次迭代次数变化情况.从归一化均方差上看，算法1收敛速度较慢，从图 3b上看存在区域性的伪影，FBP算法在左下角有亮条和暗条的伪影(图 3a)，这些现象都没有在算法2的计算结果(图 3c)中观察到，算法2的结果除了残差较低外，图像显示上效果上也较好，与模型基本一致.

在加入0.8%高斯噪声的时候，可以得到的结果如图 5所示.图 6给出了有噪声情况下归一化均方差随每百次迭代次数变化情况.

在有噪声情况下图 5的三个子图分别出现了一定程度的模糊，图 5b的区域性伪影依旧存在.从归一化均方差上看，算法2是表现最好的算法，均方差为0.05，FBP算法为0.075由于噪声的存在，如果继续迭代顺序Kaczmarz算法会显露出“半收敛”的迹象.正则化方法最后收敛的结果与FBP方法近似.需要指出，这里使用的Shepp-logan模型对于页岩样本而言过于简单，很多伪影的排列以及由算法带来的残差可能与图形的形状有关，因此有必要使用更符合页岩样本的模型来验证算法的准确性.

3.2 岩石数值模型

传统的FBP方法和Kaczmarz方法，伪影的严重程度与实际样品的复杂程度有关，一般内部结构越复杂或局部呈现周期性结构特征的样品，伪影程度越高.因此重构方法的有效性往往需要针对符合实际样品形态的数值模型检验，页岩数值模型在构建过程中主要参考了王羽等(2015)文中一般页岩的孔隙类型，如表 1所示.

一般地，黄铁矿在大部分光强下吸收系数都远高于石英、白云石等，在成像过程中这些高吸收的样品可能会对周围的像素点造成影响，为了更好地模拟实际数据，本文参照页岩基本成分的吸收系数和孔隙分布的实际情况构造了新的数值模型.在新的模型中，考虑到了可能的狭缝形、椭圆形、三角形、多边形以及高吸收的无规则黄铁矿晶间孔等各种情况，构建了如图 7所示的模型.

为了更好地检验新算法能否对细节部位进行精确成像，新模型采用了256×256的分辨率，内部结构的百分比如表 2所示.

对于图 7所示的模型，同样地，在无噪声情况下，比较了3种算法的结果分别是FBP方法，算法1和算法2如图 8(a, b, c)所示.无噪声情况下归一化均方差随迭代次数变化情况如图 9所示.

在没有噪声的情况下(如图 8a所示)，FBP方法由于左下角的狭缝型分布在右上角留下了一些伪影，此外由于左侧晶间孔的存在，使得图像在环绕晶间孔的上部和左部有一些条带状伪影.算法1的结果(如图 8b所示)有一些大面积的区域性伪影，实验时发现在没有噪声的情况下如果继续增加迭代次数，这种现象可以消除，但是在同等迭代次数的情况下，算法2的表现(如图 8c所示)就相对优秀.在均方误差收敛显示图中可以看到，如果选择了随机策略后的算法1在收敛速度上也有大幅度提高.如图 9所示.

以下几张图反映了图片不同像素条的具体数值情况，以图 10为例，对FBP，算法1和算法2在中间第62行的计算结果进行抽取，可绘制如图 10(b, c)，图b中黑色线为数值模型的真实值，蓝色为FBP算法的计算结果，青绿色为算法1的计算结果，红色为算法2的计算结果.图c为不同算法的结果与真实值差的绝对值对比，其中蓝色为FBP算法的计算结果，青绿色为算法1的计算结果，红色为算法2的计算结果，不难看出，红色线相对被蓝色线所包络，可以较好的展示采用算法2后解的稳定性.

在加入0.5%噪声的情况下，比较了FBP算法和算法2所得到的结果如图 11所示.有噪声情况下归一化均方差随迭代次数变化情况如图 12所示.

在归一化均方差上看，FBP方法为0.057，算法2为0.039，虽然在归一化均方差角度比较时算法2不占很大优势，但在图像上看，FBP方法处理的结果中有很多噪点.同样的，也对于其中部分计算结果进行了抽取绘制成图 13，其中黑色线为数值模型的真实值，蓝色为FBP算法的计算结果，青绿色为算法1的计算结果，红色为算法2的计算结果.右下角为不同算法的结果与真实值差的绝对值，其中蓝色为FBP算法的计算结果，青绿色为算法1的计算结果，红色为算法2的计算结果，算法1对噪声的抵抗能力最差，算法2最强.

4 实际数据

重庆漆辽剖面地理上位于重庆市石柱县六塘乡漆辽村，区域构造上属于四川盆地东缘石柱复向斜带，出露地层有上奥陶统临湘组、五峰组、观音桥组和下志留统龙马溪组.本实验所用页岩样本为龙马溪组页岩样本(本样本由中国科学院先导专项项目组提供).纳米CT吸收数据的实验试场约15 μm，每一个像素为50 nm，经过同步辐射光扫描投影，获取了180组不同角度下的投影数据，图 14为在某一角度下投影的实际数据，其中红线代表以下将以这一行的数据为例，利用180幅投影数据进行成像.

在几何校正和角度校正之后，得到可以处理的投影数据，利用带约束的正则化Kaczmarz算法求解反演问题，得到如图 15b所示的结果，并与利用传统的FBP方法反演得到的结果(图 15a)进行比较.图像中灰度由白到黑代表物质的对光的不同吸收，黑色代表几乎没有吸收的孔隙，灰黑色代表有机质.根据切片灰度值，孔隙主要呈星点状分布于有机质中，切片图像灰度值的分布可以用来判断不同算法所获得的图像衬度优劣，我们以图 14所示XZ方向投影数据红色线所代表的第280行数据为例，从结果可以看出本文所用算法在整体上图像边界更加清晰，对比度更强，伪影更少.

对其中部分细节的计算结果进行了提取如图 16所示，其中蓝色线为FBP计算结果，红色线为算法2的计算结果.蓝色实线是FBP算法的包络线，红色虚线是算法2的包络线，从包络线可以看出FBP算法计算结果并不稳定，这部分不应该有信息，应该得到类似于固定值的结果，但其中FBP算法波动很大呈现很明显的波峰波谷.

展示孔隙细节如图 17所示，曲线右侧上面为FBP计算结果，下面为算法2计算结果，仔细观察这两种算法在计算结果上均可以发现红色圈注处的4个低吸收像素群所表示的孔隙.但从其中一条测线的数值结果，红色线有明显的4个低吸收位置点(图中以小红圈显示)分别对应着图像上的孔隙，但FBP结果(蓝色)中可以观察到多个低吸收位置点，这是FBP算法的误差导致的，如果利用数值计算孔隙时候，这种高频的低吸收位置可能会被误处理为孔隙，进而造成孔隙度计算的误差，而算法2就不存在这样的问题.

细节放大图 18可以看出，算法2得出的孔隙更清晰一些(红色圆圈标注为3个明显的低吸收孔隙，通过同样位置的细节比较可以看出，孔隙在FBP算法计算结果上看并不清楚).

5 讨论与结论

以Kaczmarz方法为代表的代数重构方法，在包括正则参数选取，松弛因子的自适应调整，矩阵的行选择策略，约束集控制等方面具有较强的可调节性.在实际数值实验的处理过程中，有人提出选择相对较小的松弛因子可以迅速找到图像中的小的细节部分，反之则会对图像较平滑的部分收敛较快，因此在迭代初期可以采用相对较小的松弛因子，这一点在数值实验过程中可以观察到；以正则参数选取为例，当原始数据噪声较大时，可以采用较大的正则参数以突出内部结构信息；以约束集控制为例，针对不同岩石样本的地质属性，设置一定的参数范围，可以使得迭代过程中的解更稳定；从较为广义的正则化理论角度，也可以根据实际问题，利用图像内禀的诸如全变差等性质调整目标函数.在本文的各例数值实验中采用的参数也不尽相同，可以呈现较好的效果.与之相比，FBP方法主要体现在滤波器，插值方法等方面的调节，我们通过不同的数值实验发现，对于很多数值模型，调整依旧很难处理方法本身带来的伪影.而Kaczmarz方法属于投影法，即把每一次的迭代值投影到真实解所存在的解空间，边界特征和点状特征都是在真实解的空间内部，因此只要算法的收敛性有所保证且迭代次数足够，总可以收敛到一个合理的解，因此对于边界特征和点状特征的保留是投影类方法本身的属性.这也很好地解释了数值实验过程中，在无噪声情况下Kaczmarz方法可以去除由FBP方法造成的伪影这一现象.正则化技巧(特别是增广正则矩阵)的引入是本文的核心，它可以有效地压制噪声，使得一些不相容、不适定的问题能求出一个相对较准确的解.从另一个侧面理解正则化，也相当于是加入了“模型是能量小的”的先验信息，与此同时还可以加一些其他先验信息如非负性等约束条件得到准确的结果.

综上所述，本文在Kaczmarz算法的基础上综合了块状迭代思想(以提高收敛速度和稳定性)、随机投影思想(提高收敛速度)，提出了一种新的基于物理约束的正则化Kaczmarz算法较稳定地获得了页岩微纳米孔隙结构成像，解决了FBP算法的伪影现象，避免了Kaczmarz方法的系统带状伪影，表现了很好的噪声抗干扰性，在有无噪声的时候都可以对图像细节进行较为精确地刻画.在接下来的工作中，我们将把算法应用到同时处理吸收数据和相位数据中，及如何将两部分的数据结合以提高成像精度，以及对于Kaczmarz算法的继续改进将做进一步的深入探讨.

致谢  感谢评审人细致的修改建议.