目前，尽管岩石物理实验可以描述多孔介质中的微观机理，但随着研究的不断深入，岩石物理实验方法已不能完全满足一些特殊储层的研究需要，如天然气水合物储层岩心研究、裂缝发育的碳酸盐岩储层的岩心获取问题、深部钻探岩心物性分析等(姚军和赵秀才，2010).随着计算机应用技术的发展，数值模拟方法逐渐成为研究岩石物理性质的重要辅助手段，而数字岩心技术是数值模拟的前提，其建模的好坏直接影响岩石物理数值计算的结果.

数字岩心建模方法主要有两种，第一种是物理实验方法构建数字岩心，即采用物理实验方法对岩心进行扫描，用实际岩心扫描数据直接建立三维数字岩心(Lymberopoulos and Payatakes，1992).第二种是数值模拟方法如高斯场法(Joshi，1974; Quiblier，1984)、模拟退火法(Hazlett，1997)、过程法(Bakke and Oren，1997)、顺序指示法(Keehm，2003)、多点地质统计法(Okabe and Blunt，2005)、马尔可夫链-蒙特卡洛法(Wu et al.，2004)等.其中，模拟退火算法是20世纪80年代初期发展起来的一种求解大规模组合优化问题的随机性方法，与高斯模拟法相比，模拟退火法的优势在于它能够将更多的岩石信息考虑进来，使所建模型与真实岩石在孔隙空间结构上更接近(姚军等，2007).过程法建模结果在孔隙连通性上优于模拟退火法(Bakke and Oren，2003)，但过程法建模过程较为复杂，更适合于成岩过程简单的岩石(聂昕，2014).多点地质统计法所建的数字岩心孔隙连通性好，但计算速度较慢且只适合于各向同性介质(张挺等，2009).可见，每种建模方法都有自身的优点与局限，基于优缺互补的原则，近期陆续发展了模拟退火法与其他建模方法相结合的混合数字岩心建模法，如将高斯模拟法的输出作为模拟退火法的输入，使得建模时间大大减少(Hidajat et al.，2002)；针对过程法的缺陷，刘学锋等将过程法和模拟退火法相结合，提出了三维数字岩心重构的混合法(刘学锋，2010)等.

模拟退火法作为较早建立数字岩心的一类方法，它能包含反映岩石空间结构特征的较多信息，与其他方法相比有不可替代的优势，但目前应用中的问题主要有建模效率低下、建模结果孔隙连通性差、存在过多的孤立点等(赵秀才等，2007).本文主要针对传统模拟退火法运行到后期出现的问题，设计了基于模拟退火算法的新系统生成方案，作为后续补充优化算法，并取得很好效果.对于其他随机建模方法，如多点地质统计法、马尔可夫链蒙特卡罗法等，均可用此后续补充优化方法来做进一步的细致优化，有较好的实用性和移植性.

假设孔隙度为φ的岩石中只考虑孔隙与岩石骨架两相系统，并用1来表示孔隙相，0代表岩石骨架.那么产生一个0~1之间的随机数r，若r在(0，φ)之间，那么此点Z(r)便为孔隙点，否则为骨架点.

对模拟三维岩心x，y，z三个方向依次扫描赋值便可得到完全随机的初始岩心模型.

多相系统中第j相的两点概率函数定义如下：

两点概率函数可认为是在系统中任取两点，这两点同时分布于同一相中的概率.对于各向同性的介质而言，s₂^j(r₁，r₂)只取决于r₁，r₂两点的距离r=|r₁-r₂|，因此(2)式可以表示为

对只考虑孔隙与骨架的两相系统而言，建模过程中通常以孔隙相为研究对象，故以S(r)来简化表示.对于各项异性系统，为了尽可能多地获取各个方向的介质信息，可选取不同的方向来计算两点概率函数值，但计算量也会增加.对于各项同性介质，为了节省运算时间只需沿一个主轴方向计算即可.

线性路径函数L⁽ⁱ⁾(r₁，r₂)是描述岩石微观孔隙结构的重要函数，其定义如下：

由以上定义可以看出，线性路径函数体现了系统的连通性信息，在各向同性介质中，线性路径函数的值取决于r₁，r₂这两点之间的距离r，故表达式可以简化为L⁽ⁱ⁾(r)，对于孔隙相为φ的系统有：L^(j)(0)=φ.

模拟退火算法(SA)是模拟物理退火过程的一种最优化算法(Metropolis et al.，1953).据统计热力学研究表明，在某一特定温度T下，物体原子能量的分布概率P满足如下关系：

当系统温度经过退火而缓慢降低后，将达到能量最低的平衡状态.SA算法就是模拟系统在退火过程中原子能量的概率分布来进行优化计算的.定义在第(K+1)次搜索时原子能量概率分布的接受准则(Metropolis准则)为

相应状态的接受准则为

满足式(8)的状态是第(k+1)次状态.这样搜索既能向性能指标优化的方向迭代，又有一定概率接受性能指标劣化的状态.开始时温度较高，接受劣态的概率大，随着迭代次数增加，系统得到改善，逐渐达到全局最优解.

本文所建数字岩心模型是仅考虑骨架和孔隙的两相系统，因流体的储渗空间为孔隙，故以孔隙空间为研究对象计算各特征函数值.建模目标函数设置为重构系统与参考系统的自相关函数与线性路径函数差值的平方和，某时刻系统的目标函数可定义为

数字岩心建模步骤如下：(1)对实际岩心CT图像或电镜扫描图像进行二值化处理，统计建模基本资料φ、S (r)、L(r);(2)随机生成孔隙度为φ的三维数字岩心(0代表岩石骨架，1代表孔隙)，并计算初始模型目标函数E，进入优化进程;(3)在孔隙和岩石骨架空间任意选取一孔隙点和骨架点，两者互换位置生成新系统，并计算目标函数值E′;(4)新系统由Metropolis准则来判定是否取代旧系统;(5)根据内循环终止条件和外循环终止条件来决定是否终止循环;(6)最后输出数字岩心模拟结果，结束进程.

上述传统模拟退火法建立数字岩心，开始时对完全随机的模型进行优化，温度最高，由Metropolis准则可以看出，劣态系统的接受概率大，这将有效避免系统陷入局部最优，导致循环终止.随着温度的降低，迭代次数的增加，系统孔隙结构的形态逐步形成，但分布较为分散，连通性差，孤立的孔隙点或骨架点较多.如果再按传统的随机交换骨架点和孔隙点作为新系统的生成方案，就会导致原有已形成的孔隙结构遭到破坏，孔隙结构陷入形成-破坏-再形成的恶性循环，无形中增加了不必要的计算量.而孔隙中的骨架点和骨架中孤立的孔隙点都是不可能大量存在的，新系统生成方案选取的对象点主要以孤立孔隙点和骨架点为目标.故本文设计了两个阶段的优化方案，第一阶段用传统的SA方法对初始的随机模型做优化.第二阶段以择多算子选取的对象点与孔隙骨架的边界点交换生成新系统，作为补充优化阶段.

新系统生成方案的对象点选取标准是以D3Q19模型网格结构(Michalewict et al.，1996)为基础，借助二值模式下的19邻域来定义择多算子(刘学锋等，2013；赵秀才等，2008)，见图 1，网络结构的中心体素0，若在其19邻域中体素的状态值与中心体素相反的个数大于等于n，则将该中心体素状态判定为新系统生成方案选取的对象点.n的取值根据具体情况而定，若n取值过高，会导致符合条件的点越来越难选取，从而使程序在寻找孤立点上花费过多时间.若n的取值过低，则相当于传统的系统生成方案，导致优化后模型与实际岩心孔隙形态相差较大.

图 1

Fig. 1

图 1 D3Q19模型网格结构 Fig. 1 The D3Q19 lattice structure model

孔隙与骨架边界点的判定：若中心体素点的状态为孔隙，且周围18个邻域点的状态均为孔隙，则中心体素点为完全孔隙点，空间位置用(i，j，k)表示.如图 2所示，用0~5表示6个不同的方向，0代表x轴正方向，3代表x轴负方向，同理1，4和2，5分别代表y轴和z轴的正，负方向.若O点为完全孔隙点，则生成0~5之间的一个随机正整数t，例如t=5，那么计算按照z轴的负方向进行，即x、y坐标值不变，令k=k-1，若此时点(i，j，k-1)的状态仍为孔隙，则继续判定(i，j，k-2)点的状态，当找到(i，j，k-n)点的状态为骨架时循环结束，此时可以判定(i，j，k-n)点即为孔隙与骨架的边界点.若计算超出所建模型的范围且没找到满足条件的边界点，则重新随机生成一个0~5之间的整数，即选择其他方向进行计算.

图 2

Fig. 2

图 2 随机方向选择示意 Fig. 2 The choice of random direction diagram

选取新系统生成方案对象点与孔隙和骨架的边界点交换，产生新系统并计算目标函数E₁，判断是否满足Metropolis准则，若不满足则重新产生新系统，若满足则更新系统.判断是否满足终止条件，若满足则输出优化结果，循环结束.补充优化算法流程图见图 3.

图 3

Fig. 3

图 3 SA补充优化算法流程 Fig. 3 SA added optimization algorithm flow chart

图 4为某砂岩岩心CT扫描二值化图像，黑色代表岩石骨架，白色代表孔隙空间，大小为300×300像素点.提取实际岩心建模资料：孔隙度为24.93%.根据式(3)和式(5)计算数字岩心建模函数，并得到建模函数表达式，如图 5所示自相关函数与线性路径函数均随着r₁，r₂两点之间距离的增大而减小，当r=0，L(0)=S(0)=0.25=φ时，符合自相关函数和线性路径函数的性质.

图 4

Fig. 4

图 4 原始岩心二值化图 Fig. 4 The two value map of original core

图 5

Fig. 5

图 5 原始岩心的统计函数 (a) 原始岩心自相关函数; (b) 原始岩心线性路径函数. Fig. 5 Statistical functions of original core (a) Auto-correlation function of original core; (b) Linear path function of original core.

考虑到重建系统与原始参考系统的自相关函数和线性路径函数在自变量r较小时差异较大，并随着r值的增大而减小，故设定α_i，β_i的值在r<r_max/2范围是r_max/2<r<r_max范围内的2倍来突出差异，以加快优化速度.

以二维数字岩心建模为例.如图 6b所示，以(i，j)为中心点，若其周围8个邻域点与中心点状态值相反的个数大于n，则中心点即为生成新系统所选取的对象点，由图 6a可以看出棕色的9点邻域模板，若中心点的状态为孔隙，其余的邻域点中有1个点处于孔隙相，有7个点处于骨架相，且处于骨架相的点数远大于孔隙相的点数，则可判定中心点状态为孤立点.本文n的取值分别依次为5~8，n的每次取值对应一次系统的遍历扫描，直到满足循环终止条件.系统是否更新取决于系统能量是否降低.对数字岩心模型遍历扫描的次数用循环终止条来控制，当系统能量降低到足够小时可视为模型达到最优(E< 10^-6).这样第二阶段的遍历扫描运算取代第一阶段的随机搜索，程序在第一阶段进行传统SA随机搜索优化计算时效率低下，而第二阶段的遍历扫描运算则可加快优化进程，提高建模效果.

图 6

Fig. 6

图 6 对象点选取示意图 (a) 孔隙骨架分布； (b) 九点邻域模板. Fig. 6 Target selection diagram (a) The distribution of pore and matrix; (b) The nine point field template.

二维数字岩心建模大小为300×300×300像素，将SA算法中随机生成的数字岩心称为初始岩心，生成的初始岩心孔隙度为24.72%，如图 7a所示孔隙点几乎均匀分布，无法反映岩心孔隙空间的任何特征.初始岩心经传统SA 算法处理一段时间后输出的数字岩心称第一阶段建模结果，传统SA 算法在建模参考函数的约束下对初始岩心进行有效的优化，使系统能量达到最低，由图 7b可以看出模型孔隙分布较为紊乱、孔隙空间结构较为独立、连通性差、含有较多的孤立孔隙点.这是由两方面原因造成的：第一，模拟退火建模使用的建模函数有限，只能反映有限的孔隙空间结构特征; 第二，理论上只要运行足够多的时间，建模岩心就会和真实岩心足够的接近，但本文为了提高运行效率，只将传统的SA算法运行一段时间的结果输出.经补充优化算法输出的数字岩心称第二阶段建模结果，如图 7c所示，经补充优化算法处理后的结果与实际岩心相比虽不能完全吻合，但从孔隙形态，孔隙连通性上都有很大改善.三维数字岩心建模大小为80×80×80像素，建模孔隙度为24.93%.由图 8所示，最终数字岩心模型较第一阶段结果，空隙分布较为集中，连通性较好，孤立点明显减少，建模效果明显加强.

图 7

Fig. 7

图 7 二维数字岩心建模结果对比 (a) 完全随机的模型; (b) 第一阶段建模结果; (c) 第二阶段建模结果. Fig. 7 Comparison chart of 2D digital core modeling results (a) Completely random model; (b) Modeling results from the first stage; (c) Modeling results from the second stage.

图 8

Fig. 8

图 8 三维数字岩心建模结果对比(蓝色为孔隙，红色为骨架) (a) 第一阶段数字岩心模型; (b) 第二阶段数字岩心模型. Fig. 8 Comparison chart of 3D digital core modeling results (The blue color represents the pore,the red represents the skeleton) (a) Digital core model from the first stage; (b) Digital core model from the second stage.

变差函数是与图像结构相关，评价图像性质的重要函数，定义如下：

为了评价数字岩心建模效果，分别对实际岩心、第一阶段建模结果、第二阶段建模结果的变差函数对比分析.如图 9所示，在X，Y方向上第一阶段建模结果与原始岩心的变差函数有一定的偏差.而经第二阶段处理后模型显著改善，其变差函数与真实岩心的几乎完全吻合，说明最终的模型更能反映真实岩心的孔隙空间结构特征.

图 9

Fig. 9

图 9 数字岩心建模变差函数对比 (a) X方向变差函数图; (b) Y方向变差函数图. Fig. 9 Comparison chart of variation function for digital core modeling (a) The variation function for X direction; (b) The variation function for Y direction.

(1)本文定义了孔隙骨架边界点的选取标准，并以择多算子选取的对象点与边界点交换作为新系统的生成方案，设计了用遍历扫描搜索代替传统随机搜索的后续补充优化方案，取得良好的效果.

(2)基于择多算子选取孤立点参数设置十分重要，n值依次为5~8，对n为不同的取值，依次进行优化计算.与传统的SA建立数字岩心方法相比优化了建模效率，从建模结果上看孔隙空间形态和孔隙连通性上都有较大改善.

(3)对于其他的随机法建立数字岩心方法，如多点统计法、过程法等，可用本文的设计的基于模拟退火的后续优化方案来进行更细致的优化，具有很好的实用性.