0 引言

岩石物理是联系地下岩石弹性性质与储层物理性质之间关系的桥梁，定量表征了岩石体积模量、剪切模量等弹性参数与孔隙度、泥质含量、含水饱和度等物性参数之间的关系，为地震的正演模拟、反演解释提供合理的岩石物理模型，从而减少地震解释的多解性(Yin et al., 2015, 2017；Zong et al., 2015, 2017).大量理论岩石物理模型在解决各种特征(岩性、孔隙度、流体类型、裂缝等)储层岩石的问题上已经得到广泛应用(Mavko et al., 1998)，然而，由于储层参数的数目通常高于地震参数数目，储层参数的定量预测通常成为欠定问题，具有非唯一解(Johansen et al., 2013).同时，理论岩石物理模型表征的物性参数与弹性参数的关系是强非线性的，直接推导其逆运算解析表达式很困难(刘倩等，2016).因此，稳定、灵活的储层参数预测方法成为了储层定量描述和监测的关键.

为解决上述反演问题，Johansen等(2004)提出了一种岩石物理理论约束下物性参数(孔隙度、矿物含量、饱和度等)重采样的储层参数定量预测方法，该方法提供了一个从地震参数直接定量预测岩性、储层特征的反演思路，并且对于岩石物理模型的类型、输入弹性参数的选择以及之间的相关性都具有很强的灵活性，将其称之为岩石物理逆建模(Inverse Rock Physics Modeling Method, IRPM)；Moyano等(2011)应用传统交会图法和逆建模法，分别定性和定量地对比了两种砂泥岩模型对数据的解释能力；Johansen等(2013)将该方法应用于实际工区，根据预测结果评估了不同岩石物理建模方法的有效性，结果表明相比Hashin-Shtrikman-Walpole模型，选择DEM模型计算骨架模量可以提供更为稳定的弹性参数组合，预测结果更为精确可靠；Moyano等(2015)又利用了井数据的空间相关性，约束了储层参数预测范围，减少预测的不确定性；Bredesen等(2015)基于地震反演数据，利用逆建模方法识别具有相似AVO响应的富有机质页岩和含气砂岩；刘倩(刘倩等，2016; Liu et al., 2017)提出了能够进行四参数预测的扩展逆建模方法，并在目标函数中加入模型约束项，提高预测稳定性.

但是，岩石的弹性性质除了与骨架、流体的宏观弹性性质有关，还受到骨架、孔隙微观结构及孔隙尺度上各结构的相互作用等微观因素影响(Zhang and Toksöz, 2012).前人对于逆建模方法的研究均基于理论岩石物理模型构建相应的弹性参数组合，建模过程是将岩石的结构特征、孔隙度演变过程进行理想化假设，与地下真实岩石微观结构并不一致，所建立的弹性参数与物性参数间的关系与实际情况不符，因此导致反演结果存在误差.随着计算机技术与成像技术的发展，数字岩石物理技术作为传统岩石物理方法的补充手段，利用X射线CT扫描对地下岩石进行数字化成像，构建反映实际微观结构的数字岩心，并通过有限元(Arns et al., 2002; knackstedt et al., 2003; Grechka et al., 2006; Sain, 2010; 刘学峰，2010；张晋言和孙建孟，2012)、有限差分(Saenger et al., 2000, 2004, 2011；Saenger and Shapiro, 2002)等方法开展岩心尺度数值模拟，定量研究数字岩心的岩石微观结构与宏观弹性参数之间关系，弥补了传统岩石物理方法的诸多不足(刘向君等，2014).

基于前人的研究成果，本文将逆建模思想应用于数字岩心的微观领域，发展了一种数字岩心逆建模理论下的储层参数定量预测方法.首先利用数字岩心技术对岩石微观孔隙结构精确表征的优势，通过有限元数值模拟，建立了数字岩心储层参数与弹性参数的直接对应关系；将有限元模拟结果利用插值算法建立弹性参数的三维数据集，实现了弹性模量的有限元数值解的快速构建；最后通过逆建模思想，在弹性参数三维数据集中提取与输入弹性参数值相对应的模型参数组合，实现储层参数的有效预测.数字岩心逆建模方法基于对岩石实际微观结构的数字化精确建模，避免了理论岩石物理建模过程中理想化的等效过程，使物性参数的预测更为精确.最后模型测试结果表明，数字岩心逆建模方法的准确性和稳定性良好，该方法能够获得合理的数字岩心储层参数.

1 方法原理 1.1 弹性模量的有限元模拟

弹性模量有限元模拟是应用有限元法计算载荷下的数字岩心应力应变场，反推弹性参数(Garboczi and Day, 1995; Bohn and Garboczi, 2003).由于结构本身是由离散的像素组成，无需进行网格离散；单个非边界像素包含八个结点，p分量的单元位移分布近似函数u_p由八个结点p分量位移的线性插值计算：

(1)

u_rq表示第r个结点位移的q分量，N_{p, rq}为立方体单元形状矩阵，即形函数，其中，N_{1, r1}=N_{2, r2}=N_{3, r3}=N_r，N_{1, r2}=N_{1, r3}=N_{2, r1}=N_{2, r3}=N_{3, r1}=N_{3, r2}=0.

将式(1)代入像素的正、切应变表达式(2)(3)：

(2)

(3)

考虑到当p和q相等时，N_{p, rq}的值为N_r，否则为0，得到单元r结点的p分量位移u_rp对像素的正、切应变的表达形式：

(4)

(5)

其中，是一个3×8的矩阵，将8个结点三分量上的位移u_rp和像素的正、切应变联系起来，每个像素内的能量为：

(6)

其中，p，q，r，s的取值范围是1到3，应变张量含有6个独立分量，分别是：ε_xx，ε_yy，ε_zz，ε_xz，ε_yz，ε_xy，将弹性张量C_pqrs简化为C_αβ，下标α，β的取值范围是1到6.弹性势能方程(7)可转化为：

(7)

运用微分矩阵L_pq，像素应变可表达如下：

(8)

其中，S_{α, rq}(x, y, z)=L_αpN_{p, rq}(x, y, z).

将式(8)代入弹性势能方程(7)，得到：

(9)

则各像素能量E_n可用结点位移的形式表达为：

(10)

其中，D=∫d³r[S_{α, rq}]^TC_αβ[S_{β, sq}]为有限元方程的刚度矩阵.

对于边界的像素，其某些相邻像素由于周期性边界条件会位于岩心的另一端，而这些位于另一端的像素在长度上会有一个跳跃，结点之间存在位移跨度，因此引入边界后，像素的能量表达式为：

(11)

其中

(12)

δ_r为边界的校正因子，b为全局数组，C为全局常量.将系统内各像素的能量相加，组装三维数字岩心系统弹性自由能E_n.利用共轭梯度算法，使系统弹性自由能取极小值，求解系统内各结点位移分量，应保证总弹性势能梯度足够小，结点位移才接近精确值.运用结点位移使用插值法求解三维数字岩心的平均应变、应力，从而得到体积模量和剪切模量.

1.2 数字岩心三维数据集构建

数字岩心有限元数值模拟方法建立了从储层参数到弹性参数的转换关系，用矢量形式表示为：

(13)

其中，m=(m₁, m₂, m₃, …)为储层物性参数，如孔隙度、泥质含量、饱和度等；d=(d₁, d₂, d₃, …)为弹性参数，如纵横波速度、体积模量、剪切模量等；F为数字岩心弹性参数的有限元模拟.

然而，有限元模拟弹性模量时间长、效率低.本文基于Fontainebleau砂岩数字岩心，计算了不同孔隙度、泥质含量或饱和度下数字岩心的体积模量和剪切模量，如图 1所示.

大量前期测试结果发现，对于同一地区、同类型数字岩心，弹性参数的数值模拟结果为储层物性参数单调变化的函数，因此通过插值多项式拟合弹性参数随物性参数连续变化的关系.拉格朗日插值算法公式结构紧凑，便于应用于理论分析，其基本表达式如下：

(14)

其中:

(15)

l_i(x)称为插值基函数，将Fontainebleau砂岩数字岩心的有限元计算结果代入表达式，即可构建适用于该地区数字岩心的插值基函数.

当研究多种因素对弹性性质的影响，需输入大量数字岩心弹性参数作为插值节点进行二维或三维插值，为避免高阶拉格朗日插值表现的数值不稳定性，本文使用改进的重心拉格朗日插值法，数值稳定性好，插值近似精度高，其表达式为：

(16)

其中重心权定义为：

(17)

逆建模方法对输入参数d₀=(K₀, μ₀, v_s0, v_p0, ρ₀, f₀, …)的选择具有灵活性，本文将数字岩心的体积模量K₀，剪切模量μ₀和密度ρ₀作为逆建模输入弹性参数；为避免反演欠定性问题，本文选取需要预测的三个储层参数：孔隙度ϕ，泥质含量V_sh和含水饱和度S_w构建数字岩心三维数据集，则式(13)可表示为：

(18)

本文利用Fontainebleau纯砂岩数字岩心(只含孔隙和骨架)模型，基于其原始微观孔隙结构，分别对原孔隙部分的孔隙流体组成和原纯石英骨架部分的固体矿物组成进行改变.从骨架中选取一定量的像素表征泥质成分，构建不同泥质含量的数字岩心(Yin et al., 2017)；从孔隙中选取一定量的像素表征水，其余孔隙像素表征油，构建不同含水饱和度的数字岩心.孔隙介质中，大多数组成岩石的矿物是水润湿性的，润湿相存在于岩石孔隙的边角处，而非润湿相存在于孔隙空间的中心(Blunt, 2001)，因此本文在含单一流体数字岩心的基础上，将不同体积分数的水加入到数字岩心孔隙边角处，使数字岩心模型与岩石微观结构特征相一致.按照一定比例对孔隙空间、孔隙中的流体成分及骨架中的矿物成分的像素进行改变，从而构建三参数(ϕ, V_sh, S_w)同步变化的三维数字岩心.

通过上述建模方法，本文构建了一系列具有等间距储层物性参数的数字岩心，取值范围分别为0≤ϕ≤0.45(孔隙空间占岩心的体积分数)，0≤V_sh≤1(泥质占岩石骨架的体积分数)，0≤S_w≤1(水占孔隙空间的体积分数)，取值间距分别为10%，20%，20%.插值点处的弹性参数均通过多矿物成分数字岩心有限元法求取，其他不同储层物性参数的数字岩心弹性模量通过拉格朗日插值算法获得，从而建立不同弹性参数的数字岩心三维数据集，如图 2所示.

从图中可以看出，当孔隙度为0，弹性参数仅随泥质含量的增大而减小，与含水饱和度无关；当泥质含量与孔隙度一定(孔隙度不为0)，体积模量与密度均随含水饱和度的增大而增加，而剪切模量基本不受影响；当泥质含量与含水饱和度一定，弹性参数均随孔隙度增大而减小.数字岩心三维数据集建立了数字岩心储层物性参数与有限元数值模拟弹性参数的直接对应关系，为逆建模预测提供了数据基础；基于数字岩心数据集，该类数字岩心只需通过数据库对应搜索的方式，即可找到与其储层参数相匹配的弹性参数，相比于有限元直接计算方法，计算效率得到很大提高.

1.3 数字岩心逆建模方法

数字岩心三维数据集可视为基于数字岩心有限元模拟理论的不同弹性参数d_i(本文中指体积模量、剪切模量和密度)的采样标量场D_i(m)，它为弹性参数d_i取值范围内的数值提供了与其相对应的储层模型参数组合m_di(ϕ, V_sh, S_w).数字岩心逆建模思想就是在弹性参数d_i的标量场D_i(m)中进行重采样，提取与输入弹性参数值d_i0相对应的所有孔隙度、泥质含量和含水饱和度的模型参数组合.

以一组实测弹性参数(本节以K₀=18 GPa，μ₀=15 GPa，ρ₀=1.86 g·cm^-3的岩心为例进行说明)为逆建模的输入数据，分别在数字岩心体积模量、剪切模量和密度标量场中根据输入数值进行重采样，如图 3所示.

从图中可以看出，输入弹性参数K₀，μ₀和ρ₀对应数字岩心标量场中的所有模型参数组合构成了一个等值曲面，每两个单值等值曲面的交集为一条曲线，曲线上的点即为满足两个弹性参数数值的所有模型参数组合.因此，为找到与三个弹性参数均对应的模型参数组合，需将三个单值等值曲面相交，如图 4所示.

三个单值等值曲面相交可以分别得到三条交线，如图 5所示，三条等值线相交为一点(图中用蓝色圆圈标出)，该点为三组模型参数组合的空间交会，满足三个输入弹性参数K₀，μ₀和ρ₀，其坐标即为数字岩心逆建模理论下的储层物性参数.

因此，数字岩心逆建模预测储层参数需求解体积模量、剪切模量和密度对应的三组模型参数组合m_di0(ϕ, V_sh, S_w)的交集，表示为：

(19)

只需求解任意三条曲线中两条(如体积模量与密度对应模型参数组合交集，体积模量与剪切模量对应模型参数组合交集)的交点即可，建立目标函数：

(20)

其中φ_K0ρ0(m)表示体积模量K₀与密度ρ₀等值曲面的交线坐标(ϕ, V_sh, S_w)，δ为交线间的最大允许距离，即寻优范围，表征输入数据与模型参数组合之间的不确定性.δ越大，求解范围越宽.可根据实际需要设定其值的大小(本文定义为δ=0.01).求解范围内的点视为可接受点，将所有可接受点坐标求均值，即可得到三个模型参数组合的交点，即为数字岩心逆建模理论下的储层参数预测解.由此得到本节示例的最终结果为ϕ=20.60%，V_sh=49.33%，S_w=79.33%.反演的前提条件是模型，如果在给定弹性参数下没有储层参数与其对应，说明所构建储层参数与弹性参数关系不适用于该地区(Johansen et al., 2013).因此，当寻优范围过大仍不能获得最优解，则认为该弹性参数的三维数据集与实际地层情况不吻合.

地下储层岩石微观结构受多种因素控制，不同孔隙类型、矿物组分、温压条件等都会对弹性参数与物性参数之间的关系造成一定程度偏差，因此储层参数预测精度依赖于能否在上述不同地质条件下精确构建适用的弹性参数三维数据集；其次，本文提出的数字岩心逆建模预测方法利用V_p，V_s和ρ作为输入弹性参数，通常由测井资料或地震反演获得，不可避免的受到噪声、反演方法等因素影响.因此，所使用的弹性参数的准确性直接影响到后续逆建模储层物性反演的精确程度；此外，插值拟合方法中近似关系的使用，以及应用插值算法时使用的插值点个数也会对预测准确度造成一定影响.

尽管数字岩心本身即可获得孔隙度和泥质含量，但在实际工区中，储层物性特征是复杂多变的，仅对个别岩心样品进行研究无法代表该地区数字岩心物性参数的一般规律.因此需要将数字岩心有限元正演与插值拟合算法相结合，由单个数据点扩展到数据面或数据体，更加全面地建立符合储层特征的物性参数与弹性参数统计规律，从而进行数字岩心物性参数精确预测，为地震反演提供可靠的控制点处的岩石物理先验信息.

综上，数字岩心基于对岩石微观孔隙结构的精确表征，利用数值模拟方法，构建了准确的弹性参数与储层物性参数的直接对应关系，为逆建模方法提供了适用于实际工区地质条件的弹性参数三维数据集，最终实现储层物性参数的精确预测.基于数字岩心逆建模方法实现储层物性参数反演的流程图如图 6所示.

2 方法试算 2.1 数字岩心逆建模方法准确性测试

如1.2节所述，本文以孔隙度、泥质含量和含水饱和度分别为10%，20%和20%的取值间距，构建了一系列数字岩心，通过有限元算法计算其弹性模量作为插值点，利用拉格朗日插值算法求取了弹性参数随储层物性参数连续变化情况.

为了检验数字岩心逆建模方法估算储层物性参数的准确性，本次研究选取了18个不同物性参数的数字岩心，并分为两组：插值数据点和非插值数据点，其中，插值点处的孔隙度取值分别为10%、20%、30%；泥质含量分别为20%、40%、60%；含水饱和度分别为40%、60%、80%；非插值点处的孔隙度取值分别为15%、25%、35%；泥质含量分别为30%、50%、70%；含水饱和度分别为30%、50%、70%.通过有限元计算这18个数字岩心的弹性参数，并将其作为输入数据，通过数字岩心逆建模方法估算其储层物性参数，结果分别如图 7、8所示.

从图 7、8可以发现，当数字岩心物性参数值为插值参数点时，数字岩心逆建模预测结果与数字岩心实际物性参数一致，预测准确；当数字岩心物性参数为非插值参数点时，逆建模预测结果与实际物性参数对比有微小偏差，其中，泥质含量预测基本与实际值吻合，孔隙度预测值均略小于真实值，含水饱和度预测结果不如孔隙度和泥质含量准确，但误差不超过20%.

模型试算表明，基于逆建模方法的数字岩心储层参数预测是可行的，结果较为准确，非插值数据点处的估算结果与实际值有一定偏差，可以通过增加插值数据点数目以提高数字岩心逆建模的估算准确性.

2.2 数字岩心逆建模方法稳定性测试

逆建模的输入数据通常为实际测井资料或地震反演结果，由于实际地震、测井数据均存在噪声，因此逆建模所输入的弹性参数通常是包含误差的.为了进一步检验数字岩心逆建模方法的可行性，需要对该方法进行稳定性测试.

采用2.1节所述的物性参数为插值数据点的数字岩心，通过有限元方法计算其弹性参数，然后分别向其体积模量、剪切模量和密度数据中加入1%的随机噪声，每组数据重复5次，将得到的含噪声数据进行逆建模运算，预测其孔隙度、泥质含量和含水饱和度.预测结果与实际模型参数对比结果如图 9-11所示.

从图中可以发现，当弹性参数加入随机噪声之后，基于数字岩心逆建模理论的孔隙度和泥质含量的预测结果稳定性很好，与实际数据基本一致；而含水饱和度的预测结果出现了一定的波动，说明含水饱和度对噪声的加入较为敏感，预测误差较大，相比于孔隙度和泥质含量，数字岩心逆建模方法对含水饱和度的预测可靠性较低.

3 结论

本文将理论岩石物理中的逆建模思想与数字岩石物理中的微观孔隙结构精确表征优势相结合，提出了一种数字岩心逆建模理论下的储层参数定量预测方法，实现了数字岩心储层参数的有效预测.首先通过有限元模拟和插值算法构建数字岩心弹性参数的三维标量场，实现了弹性模量的有限元数值解的快速构建；然后根据输入弹性参数，分别在数字岩心体积模量、剪切模量和密度标量场中进行重采样，提取对应的模型参数组合，三组模型参数组合的交点坐标即为数字岩心储层参数.

数字岩心基于对岩石微观孔隙结构的精确表征，利用数值模拟方法，构建了准确的弹性参数与储层物性参数的直接对应关系，为逆建模方法提供了适用于实际工区地质条件的弹性参数三维数据集，从而避免了理论岩石物理建模过程中理想化的等效过程，最终实现储层物性参数的精确预测.

准确性测试表明，数字岩心逆建模理论下的储层参数预测方法准确性较高，非插值点处含水饱和度的估算结果与实际结果有一定偏差，可以通过增加插值数据点数目提高数字岩心逆建模的估算准确性；稳定性测试表明，当输入弹性参数中加入随机噪声，含水饱和度对噪声的加入较为敏感，预测结果可靠性较低，而孔隙度和泥质含量的预测结果稳定性很好，与实际数据基本一致，从而说明基于数字岩心逆建模理论的储层参数预测方法是可行的.