页岩储层中，微纳米孔隙是其孔隙结构的主体，页岩气大多赋存于微纳米孔隙中。作为衡量和评价页岩储层优劣的重要指标，页岩储层微观孔隙结构的特殊性和重要性引起了国内外的重视(Loucks et al., 2009; Javadpour, 2009; 滕吉文和刘有山, 2013)。对于页岩的直接数字成像是最为有效的微观孔隙结构研究方法，例如，透射电子显微镜(TEM)、扫描电子显微镜(SEM)、CT扫描技术、扫描探针显微镜技术等。聚焦离子束扫描电子显微镜(FIB-SEM)技术已经广泛应用于页岩的微结构成像(Chalmers et al., 2012; Loucks et al., 2009; Ambrose et al., 2010; Passey et al., 2010; Schieber, 2010; Sondergeld et al., 2010; Zhang and Klimentidis, 2011)，但是FIB-SEM技术刻画岩石的三维微结构特征是有损的且耗时。CT成像技术是一种无损的、高分辨率、高效率、制样简单的微观成像技术，能够方便快捷的获得页岩内部结构图像，在页岩微观结构分析中已经得到良好的应用(Panahi et al., 2012; Iglauer et al., 2013; Riepe et al., 2011; Mayo et al., 2015; Peng et al., 2015; 刘学锋等, 2013; 邹才能等, 2011, 2012; 崔景伟等, 2012)。通过这些新发展的先进三维成像技术，人们研究了微纳米孔隙的大小、类型、形状、连通性等特征，比较分析了北美主要页岩储层在微纳孔隙结构特征的异同，有效指导了页岩气的勘探开发。我们研究了基于同步辐射X射线的页岩的微纳米结构CT成像技术，取得了高分辨率的三维成像效果(Wang et al., 2016)。在获得页岩微纳米成像数据之后，另一个重要的科学问题是如何对这些三维数据体进行孔隙结构成像。

图像分割是指将图像分成各具特征的目标区域并提取出感兴趣的目标的技术和过程。图像分割方法有很多种，适用于不同的问题，总的来说可以分为：按利用区域内灰度值相似性的基于区域的分割方法，利用区域内的灰度值的不连续性的基于边界的分割方法，以及基于特殊理论的分割方法(何俊等, 2009)。图像分割方法已经在医疗CT图像、模式识别等领域得到了广泛的应用。如前所述，我们研究的页岩成像数据是CT成像数据，获得的是页岩内部成分对高能同步辐射X射线的衰减，是物理量，因而传统的分类方法很难获得有物理意义的分类结果。此外，传统X射线CT扫描并不总是足以区分岩石组成成分，例如当不同的成分在同一能量X射线下显示出相似的X射线吸收特性时。而基于多能量(光谱)的X射线CT扫描的数据约束建模(DCM)可以用来表征材料微观组成分布(Yang, 2012)。在该方法中，假设体素的总体积是其构成各组分的子体积的总和，并假设当X射线束通过体素时，总衰减是其组成成分的衰减的总和。后一个假设在体素大小很小时是有效的。基于这两条假设，采用多能量X射线CT扫描数据作为约束，Yang (2012)建立了非负约束的线性或非线性优化问题并采用相应的优化方法求解。优化过程则是最小化预期与X射线吸收系数之间的差异。然而这类方法要求反复求解最优化问题实现对每个页岩切片微结构的图像分割，计算效率较低。神经网络作为近年来热门数据挖掘技术，具有很强的泛化能力和数值逼近能力，在图像处理领域已经取得了瞩目的成就(Bao et al., 1998; Chen et al., 1990)。

本文研究基于神经网络的页岩微纳米孔隙微结构的CT图像分割方法。我们的研究思路是：第一步，获得不同能量的同步辐射X射线扫描页岩样本得到的CT成像数据(Wang et al., 2016)作为BP神经网络训练数据的输入数据；第二步，建立基于不同能量同步辐射X射线扫描CT图像的非线性数据约束模型(DCM)，求解上述模型获得页岩内部孔隙的图像分割结果，作为BP神经网络训练的标签数据；第三步，建立多层神经网络和误差逆传播(BP)正则化模型；第四步，研究求解上述神经网络模型极小化的最优化算法，得到训练好的神经网络模型；第五步，给出基于BP神经网络的页岩微纳米孔隙分析结果。

1 基本理论 1.1 人工神经网络模型

人工神经网络(Artificial Neural Networks)指的是由具有适应性的简单单元组成的广泛并行互连的网络，它的组织能够模拟生物神经系统对真实世界物体所作出的交互反应(Kohonen, 1988)。人工神经网络中最基本的组成是神经元。在生物神经网络中，每个神经元与其他神经元相连，当它被激活时，就会向相连的神经元发送化学物质，从而改变这些神经元内的电位；如果某神经元的电位超过了一个“阈值”，那么它就会被激活，向其他神经元发送化学物质(周志华, 2016)。

M-P神经元模型是最常用的神经元模型(如图 1所示) (McCulloch and Pitts, 1943)。MP神经元接收到来自其他神经元传递过来的输入信号x_i(i=1, 2, …n)，这些输入信号经过加权求和之后与神经元的阀值进行比较，最后通过“激活函数”处理，产生本神经元的输出，常用的激活函数如图 2所示。

将上述的神经元按照一定的层次结构连接起来就构成了神经网络。其中，最简单的人工神经网络模型—感知机模型由两层神经元组成：输入层和输出层。输入层接收外界输入信号，输出层是M-P神经元接收来自输入层的信号。给定训练样本(x, y)，可以通过训练获得权重w_i(i=1, 2, …, n)以及阈值θ。将阈值θ视作输入固定为-1的连接权重w_n+1，则可以将权重和阈值统一起来学习，权重可以通过以下方式迭代更新：

(1)

(2)

其中，η∈(0, 1)称为学习率，感知机的输出为，权重值将迭代更新直至=y(周志华，2016)。

值得指出的是，感知机只拥有一层功能神经元，限制了感知机的学习能力。感知机可以用于解决线性可分问题，但是对于非线性可分问题，感知机的学习过程不能收敛到正确解。这就需要考虑多层神经网络。

常用的多层前馈神经网络中，输入层与输出层之间的神经元称为隐含层(Hidden layer)。常见的多层前馈神经网络是如图 3所示的结构，每层神经元与下层神经元全连接，。其中输入层神经元用于接收外界输入；隐层与输出层负责处理功能。神经网络的训练过程，就是根据训练数据来调整神经元之间的权重以及每个功能神经元的阈值，来达到学习的目的。本文采用三层前馈神经网络模型，通过下文所述的BP误差反向传播算法进行训练。

1.2 误差反向传播(BP)算法

训练多层前馈神经网络的算法中，误差反向传播(BP)算法是迄今最成功的神经网络学习算法之一，且在实际任务中广泛应用。常用的“BP网络”通常是指用BP算法训练的多层前馈神经网络。

对训练样本(x_k, y_k)∈T，T为训练集，假定神经网络的输出为_k=(₁^k, ₂^k, …_l^k)，则网络在(x_k, y_k)上的均方误差为

(3)

BP算法训练BP网络的步骤如下：

算法1. (BP算法)

1) 权值初始化

网络训练开始时网络权值为未知数，本文采用较小的随机数作为各层权值的初始值。

2) 网络正向传播

计算各层神经元的输出值

(4)

其中，f是第j层激活函数，为图 2所示的阶跃函数或者sigmoid函数；w_hj表示权值，θ_j表示阈值，y_j表示第j层神经元的输出。通常将将权重和阈值统一起来学习，变量统一用w来表示。

3) 网络反向传播，更新权值

BP算法可以基于梯度下降算法，以目标的负梯度方向对参数进行调整。对公式(3)的误差E_k，给定学习率η，按如下公式迭代更新：

(5)

如果达到预先给定的停止条件, 则输出由权值和阈值决定的BP神经网络；否则，更新权值和阈值，返回步骤2)。

1.3 Tikhonov正则化

BP神经网络的目标是要最小化训练集T上的累积误差。(Hornik et al., 1989)证明，只需两个包含足够多神经元的隐层，多层前馈网络就能以任意精度逼近任意复杂度的连续函数。

但是，BP神经网络在最小化累计误差的过程中经常遭遇过拟合，即训练误差很小，但测试误差却可能很大。可以通过以下两种策略缓解BP算法的过拟合问题(周志华，2016)。第一种手段是将数据分成训练集和验证集，训练集用来计算梯度、更新权重和阈值，验证集用来估计测试误差，若训练集误差降低但验证集误差升高，则停止训练，同时返回具有最小验证集误差的连接权和阈值，在该策略中，迭代步数起着正则化参数的作用(Wang, 2007)。第二种处理手段是对在目标函数中对连接权和阈值施加先验约束，也就是我们通常所说的正则化(Tikhonov and Arsenin, 1977; 肖庭延等, 2003; Wang, 2007)。对于一个非线性算子方程F(x)+ε=d(ε为误差项)，设ψ(x)为模拟数据与观测数据的拟合，则一般的Tikhonov正则化形式如下：

(6)

其中，第二项是约束项，P为规范化算子(通常是正定或半正定的)，α>0为正则化参数，用于数据拟合项与正则化项进行折中。通常正则化参数可以通过先验选取或解后验偏差方程获得(王彦飞, 2007)。由于F(x)的非线性特性，通常需要做线性化处理后进行迭代计算。

1.4 Levenberg-Marquardt算法

如公式(3)所示，目标函数E_k(w)表示神经网络在训练集上的误差，是关于权重w和阈值θ的函数。此处，将权重和阈值统一起来学习，变量统一用w来表示。因而神经网络的训练过程可看作一个参数寻优过程，即在参数空间中，寻找一组最优参数w^*和θ^*使得E_k(w)最小。

基于梯度的搜索方法是使用最为广泛的参数寻优方法(袁亚湘, 1993)。在此类方法中，我们从某些初始解出发，参数沿着与误差梯度相反的方向移动，使误差函数减小，直到取得极小值，它的主要计算代价是计算目标函数的梯度。但是，这种基于梯度下降方法的线性收敛速度很慢。本文采用Levenberg-Marquardt算法(L-M)(肖庭延等, 2003)对结果与目标数据的均方误差进行最小化，得到网络中的权值参数。业已证明，L-M算法是一种正则化方法，可以用来稳定地求解参数寻优这类特殊的反问题(Wang and Yuan, 2005; 王彦飞, 2007)。下面对L-M算法作简要阐述。

设误差目标函数的极小化形式为：

(7)

其中，e(w)的每个组份可以写成e_j(w)=||y_j-_j||(j=1, 2, …, l)，e_j(w)为误差；y_j表示期望的网络输出向量；_j为实际的网络输出向量；l为输出的维度；w为网络权值和阈值所组成的向量。根据前面所述的Tikhonov正则化思想，我们建立如下的正则化模型

(8)

其中，Δw表示权值增量，参数R>0为广义球半径。利用Lagrange乘子μ求解(8)的优化算法称为L-M算法。本文，本文取P为恒等算子，Lagrange乘子μ也可以称作学习率。公式(8)是极小化问题(7)的子问题形式，并且是一种正则化算法(Wang and Yuan, 2005; 王彦飞, 2007)。

设w^k表示第k次迭代的权值和阈值所组成的向量，新的权值和阈值所组成的向量w^k+1为w^k+1=w^k+Δw。求解极小化问题(8)，得到权值增量Δw的计算公式如下(肖庭延等，2003)：

(9)

其中，I为单位矩阵；μ为用户定义的学习率；J(w)表示Jacobian矩阵，写成下面的形式

(10)

算法2. (基于L-M算法的BP神经网络计算过程)

输入：训练集数据D={(x_k, y_k)}_k=1^l，给出训练误差允许值ε，常数μ₀和β(0 < β < 1)；

迭代计算：

1) 采用较小的随机数初始化权值和阈值向量，令k=0, μ=μ₀；

2) 计算网格输出及误差指标函数E(w^k)；若E(w^k) < ε，转到7)；

3) 计算Jacobian矩阵J(w^k)；

4) 计算Δw；

5) 更新权值w^k+1=w^k+Δw；

6) 计算误差指标函数E(w^k+1)；

7) 若E(w^k+1) < E(w^k)，则令k:=k+1，μ:=μβ转到2)，否则令k:=k+1，μ:=μ/β，转到2)；

8) 输出网络权值，停止迭代。

算法2的第4步要求正确地求解公式(9)，即保证公式(9)中求逆矩阵的正定性。因此算法2需要输入合适的初始参数值μ₀，这通常通过重复倍增μ₀的值得到。

2 训练BP神经网络 2.1 生成训练数据

建立训练数据是实现BP神经网络的重要步骤，这依赖于研究的问题。我们建立的训练数据的输入数据是基于多能量(光谱)的X射线CT扫描的页岩重构后的成像数据(Wang et al., 2017; 王羽等2015)，标签数据是页岩内部孔隙结构图像。标签数据可以通过以下描述的方法计算获得。

在多能量X射线CT扫描数据约束下，假设体素的总体积是其构成材料的子体积的总和，并假设当X射线束通过体素时，总衰减是其组成成分的衰减的总和，则可以建立基于数据约束的DCM数学模型来表征材料微观组成分布(Yang, 2012)。

基于重构后的页岩的多能量X射线CT扫描数据，定义一个立方体网格，体素位于每个网格点上。在第n个体素(n=1, 2, …, N，其中N是立体网格数据上体素的总数目)，构建这个体素目标函数如下：

(11)

这个方程相当于最小化预期的线性吸收系数与测量的吸收系数之间的差异。其中ξ是正则化参数，Ω_n是解的先验约束(肖庭延等, 2003)。对于每个X射线束能量(谱)l，预期的线性吸收系数与CT重构后的线性吸收系数的差异可以表达成如下公式，

(12)

其中，v_n^(m)表示第m个组分的体积分数；μ_n^l表示由第l个CT数据集得到的在这个体素内的总线性吸收系数；μ^{(m, l)}表示第m个组分在第l个CT数据集对应的X射线束下的单一物质吸收系数。

我们考虑下述的极小化问题

(13)

通过调节体积分数变量v_n^(m)(m=0, 1, …, M)，使目标函数T_n达到最小。

先验约束条件可以表示成如下形式(Yang, 2012)，

(14)

其中，S^m表示物质成分的自能量(化学特性)，I_k^{(m₁, m₂)}代表距离为k的相邻体素成分m₁和m₂之间的界面自由能，当k=0时，表示体素内部。特别的，当k=1时，表示最临近的体素；当k=2时，表示次临近的体素，以此类推。最大相互作用范围为K。相邻距离为k的相邻体素的数目为N^(k)，并且相邻距离为k的第j个邻接向量表示为n_j^(k)。这些参数由用户定义，可以基于先验信息和其他测量手段获得(Zschornack, 2007)。

通过非线性规划算法求解约束优化问题(13)(袁亚湘, 1993)，得到训练数据v_n。逐体素对页岩内部孔隙含量进行计算，得到结果作为神经网络训练时的输出数据。

2.2 创建并训练神经网络

本文研究的问题中BP神经网络的输入为2*1的向量，是由能量为18keV以及能量为20keV的X射线下生成的CT数据经过重构得到，在第n个体素下的数值组成的向量。标签数据为一个标量值，是利用DCM模型通过非线性规划算法计算出的页岩组份微结构图像数据在体素n下的数值(即在第n个体素下的孔隙体积百分比)。选择用三层的前馈神经网络计算。其中，输入层维度为2，隐层的神经元个数均为10，输出层神经元个数为1；输入层、隐藏层、输出层之间均为全连接；隐藏层和输出层的激活函数均选择tansig函数。如下图 4所示(图中，x表示输入数据，w表示权重值，(表示阈值，f是激活函数，y是各层的输出)。

下面给出多能谱页岩CT数据微纳米孔隙结构训练的流程。

算法3：基于神经页岩微纳米孔隙结构训练流程

输入：训练集数据D={(x_k, y_k)}_k=1^m，其中，x_k为多能CT图像第k个体素的值组成的列向量，y_k是第k个体素内孔隙含量值。

训练：给出训练误差允许值ε，常数μ₀和β(0 < β < 1)，

1) 并且用较小的随机值初始化权值和阈值，令k=0, μ=μ₀；

2) 对训练数据集进行随机选取，选取其中70%为训练数据，30%为验证数据，30%为测试数据；

3) 采用算法2对本神经网络模型进行训练，得到训练好的网络模型。

2.3 实验结果

通过上海光源同步辐射CT成像装置，在18keV和20keV能量的X射线下对页岩样本进行了同步辐射扫描，并经过图像重构算法进行重构得到页岩的CT图像(Wang YF et al., 2017; 王羽等2015)。选取其中的一片数据，通过DCM模型分析其孔隙结构，得到训练数据。通过算法3对图 4中的神经网络进行训练，得到训练好的神经网络模型。图 5、图 6、图 7所示为网络训练过程中的误差以及迭代性能分析。图 5所示为训练数据、验证数据、测试数据的误差分布，可以看出，误差均分布在0附近。图 6所示为样本均方误差随迭代次数的增加逐渐减少的示意图，从图中可以看出，采用正则化L-M方法，误差下降较快，最终稳定在10^-4量级，并且验证集与测试集与训练集的迭代误差同步递减，没有出现过拟合现象。图 7所示为训练集、验证集与测试集的数据与标签数据的拟合情况对比，从图中可以看出，绝大多数样本数据都能很好的拟合标签数据，少数样本会出现较大误差。

测试实验中仍然采用同一区域的岩石样本在18keV和20keV的X射线下的CT数据(201*201)作为输入。得到的结果如图 8所示，其中图 8a, b分别为在18keV和20keV能量的X射线下，页岩样本同步辐射扫描重构后的图像；图 8c为通过神经网络计算所得的孔隙结构图像；图 8d为标签数据。通过计算，图 8c, d的均方误差在10^-4左右，达到了精度要求。

从图 8中可以看出，神经网络可以很好的给出页岩内部孔隙结构的图像分割结果。取同一块页岩样本不同切片的图像数据进行多次测试实验证明，均方误差均在10^-4量级。并且一旦训练好网络，对同类型的页岩样本进行孔隙结构分析，对与201*201大小的图片数据，需要0.1s左右时间。如果采用DCM模型用非线性规划计算，用时是神经网络算法的数百倍，因而后者显著提高了计算效率。

3 结论

本文根据不同物质在不同能量的X射线下的吸收系数不同，对采用不同能量下的同步辐射光CT图像进行图像分割得到页岩内部孔隙分布图像。并将BP神经网络算法用于页岩多能谱CT图像的孔隙结构分析中，可以通过多能谱CT图像直接得到页岩内部孔隙分布情况，定量描绘孔隙的大小、分布、以及联通情况。相比于传统页岩CT图像分割方法(如求解DCM模型的优化算法)，计算效率得到了显著地提高。同时，神经网络模型可扩展性强，通过对不同页岩样本图像分割方法或页岩数字成像数据的学习，可以优化神经网络刻画页岩微纳米孔隙结构的精度和准确性。

致谢 感谢中国科学院先导科技专项(B类)提供页岩样品和上海光源提供同步辐射数据。感谢审稿专家的宝贵建议。