基于X射线CT图像的泥质致密砂岩纵横波速度计算

P⁃ and S⁃Wave Velocity Calculation Using X⁃Ray CT Images for Shaly Tight Sandstone

摘要

基于三维X射线CT图像重构的数字岩心模型，是目前计算岩石纵横波速度的重要方法，但由于CT图像难以准确区分致密砂岩中泥质和岩石碎屑颗粒，利用X射线CT图像研究泥质致密砂岩纵横波速度难度较大.本研究通过建立二维数字岩心模型，开展有限元模拟，讨论泥质含量、分布形式以及微孔隙的发育程度对岩石弹性参数的影响，基于影响岩石弹性参数的主控因素，探索利用三维数字岩心模拟泥质砂岩纵横波速度的可行性.研究表明：分散泥质以及骨架结构泥质对岩石体积模量影响较小，颗粒‒颗粒接触面结构泥质对体积模量影响较大，三种泥质分布形式对剪切模量影响相当；泥质含量对岩石体积模量的影响小于泥质分布形式，而泥质含量对剪切模量的影响大于泥质的分布形式；与分散泥质相关的微孔隙对岩石弹性参数影响较小，而岩石弹性参数对与骨架结构泥质和颗粒‒颗粒接触面泥质相关的微孔隙敏感，且微孔隙增加对剪切模量的减小大于体积模量.根据上述模拟结果，针对泥质砂岩建立了基于三维分水岭算法的颗粒‒颗粒接触面泥质模型，模拟结果与实测值吻合较好.

Abstract

It is an important method to simulate P- and S-wave velocities using digital core obtained from X-ray CT images. However, since it is impossible to differentiate clay and grains, and abundant micro-pores exist in shaly tight sandstones, simulating P- and S-wave velocities using digital core is a challenge. In this study, 2D digital core models are constructed and simulated using finite element method to understand the effect of clay content and distribution, amount of micro-pores on the rock elastic properties. The results will be used to assist the construction of a 3D model that can be used to simulated P- and S-wave velocities of shaly sandstones. The results show that dispersed clay and framework clay have minor effects on the bulk modulus, while interstitial clay shows large effect on bulk modulus. Effects of the different clay distributions on shear modulus are similar. Clay distribution has larger effect on bulk modulus than clay content, whereas clay content has larger effect on shear modulus than clay distribution. Micro-pores related to dispersed clay have minor effect on rock elastic properties, however micro-pores related to framework clay and grain-grain contact clay are sensitive to rock elastic properties. In addition, micro-pores have larger effect on shear modulus than bulk modulus. Based on the above results, a 3D digital core model using 3D watershed method on the X-ray images has been built and the results show good match with the measured velocities.

Graphical abstract

关键词

数字岩心 / 泥质致密砂岩 / 有限元 / 分水岭算法 / 纵横波速度 / 地球物理 / 油气地质.

Key words

digital core / shaly tight sandstone / finite element method / watershed method / P⁃ and S⁃wave velocities / geophysics / petroleum geology

0 引言

岩石的纵横波速度与岩石的成分、孔隙度、流体类型及含量以及裂缝有密切的关系，因此纵横波速度是进行测井、地震储层评价、流体识别的重要参数（Wyllie et al.，1958；Nur and Simmons，1969；Han et al.，1986；Liu et al.， 2022； 乔汉青等，2023；程绩伟等，2024）.岩石的纵横波速度目前主要通过实验室岩心测量获取，但随着数字岩心技术的发展，在数字岩心建模的基础上，利用数值模拟算法计算岩石弹性参数，进而获取纵横波速度逐渐成为重要的研究手段.Arns et al.（2002）首次利用micro⁃CT对枫丹白露砂岩进行了三维扫描，利用有限元方法计算了体积模量与剪切模量，数值计算结果与Han et al.（1986）的实验结果一致.刘学锋（2010）、姜黎明（2012）等也对枫丹白露砂岩进行了弹性参数模拟，与实测值吻合较好.Saxena and Mavko（2016）利用岩石的二维图片模拟了岩石的弹性参数，并与三维数字岩心模拟结果相比较，发现二者有较好的一致性，认为利用二维图片模拟岩石弹性参数具有可行性.上述模拟主要针对纯净砂岩或碳酸盐岩，Shulakova et al.（2013）尝试利用数字岩心技术模拟含泥质砂岩的弹性参数，其研究利用扫描电镜图像统计了泥质含量并首先计算了基质岩石的弹性参数，然后基于CT扫描图像，应用有限元方法模拟岩石弹性参数.该方法考虑了泥质含量的影响，但未考虑其分布特征以及与泥质相关微孔隙发育的影响.总体而言，目前开展的岩石弹性特征的数字岩心模拟主要针对纯净砂岩，而对含泥质砂岩讨论相对较少.

利用micro⁃CT图像建立泥质砂岩数字岩心的难点主要体现在两个方面：其一是泥质难以识别.micro⁃CT技术识别矿物主要通过矿物之间密度的差异，而致密砂岩中石英颗粒与泥质的密度相当，因此难以区分颗粒和泥质，而泥质的含量（Han et al.，1986）、分布形式（Minear，1982；Dvorkin and Nur，1996；Sams and Andrea，2001；Xu et al.， 2020）等因素对岩石弹性特征影响较大.其二是泥质的存在导致大量微孔隙发育（Swanson，1985；曾文冲，1991）.微孔隙既包括CT分辨率以下的纳米尺度孔隙，也包括CT能识别的微米尺度（1~30 μm）孔隙.前者由于无法识别，前人在建模过程中没有考虑（Arns et al.，2002；Shulakova et al.， 2013），而后者的发育则会导致建模过程中网格复杂，增加数值模拟的时间.如何建立合理的含泥质致密砂岩数字岩心模型，开展数值模拟工作是研究的难点.

本研究将在对泥质的含量、分布形式以及微孔隙发育特征分析的基础上，分析影响致密砂岩弹性特征的主控因素，进而有针对性地开展致密砂岩数字岩心建模及模拟，最后分析利用数字岩心进行泥质砂岩声波速度模拟的可行性.

1 不同泥质分布形式数字岩心建模及模拟方法

1.1　建模方法

由于CT扫描得到的图像能够较好地显示颗粒的大小、接触关系以及孔隙结构特征，因此泥质的含量、分布形式以及微孔隙的发育特征是影响建模精度最重要的因素.

首先考虑泥质的含量以及分布形式的建模方法.Minear（1982）将泥质的分布形式划分为3类：层状泥质、分散泥质以及结构泥质，其中结构泥质又划分为骨架结构泥质和颗粒‒颗粒接触面结构泥质，见图1，其中灰色为岩石颗粒，黑色为泥质，白色为孔隙.层状泥质与其他泥质分布类型实际上属于不同的尺度（王克文，2007），本研究暂不讨论.对于胶结物分布形式的模拟，目前主要有两种方式.一种是以（椭）球状或（椭）圆形颗粒为基础的胶结物模拟方法（姜黎明，2012；Jin et al.，2012），但该方法主要用于模拟规则的颗粒形状，能够分析高孔高渗储层的特征.为了研究致密砂岩储层的特征，本研究将在CT扫描图像构建数字岩心的基础上，建立胶结物的分布模型，建模的方法依据Torskaya（2013）的建模方法.

建模过程中为了减少大规模的运算，采用二维图像进行模拟研究.同时，为了模拟泥质含量变化的影响，选择孔隙度较高的样品（图2）.图2a为原始图像，图2b为二值化之后的图像，其中灰色代表颗粒，黑色代表孔隙，样品面孔率为16.65%.

Torskaya（2013）模拟了三种胶结物分布形式，包括颗粒包壳模型、孔隙分布模型和喉道分布模型.模拟胶结物的分布以及体积受以下几个参数控制：（1）到孔隙中轴线最近的距离，d_s ：定义为像素点到孔隙中轴线最近的距离；（2）孔隙‒颗粒距离，d_e ：定义为孔隙像素到最近的颗粒像素的距离；（3）孔隙直径，d_m ：定义为包含该像素点的孔隙最大内切球（圆）直径.

颗粒包壳模型、孔隙分布模型和喉道分布模型三种不同分布模式的实现条件为：（1）孔隙分布模型：像素点到孔隙中轴线最近的距离大于某一个值，d_s >ε.（2）颗粒包壳模型：像素到颗粒边缘的最近距离小于某一个固定值，d_e <ε.（3）喉道分布模型：包含像素点的最大内切球的直径与像素到颗粒边缘的最近距离的乘积小于某一固定值，d_e ×d_m <ε.其中ε反映胶结物的体积（厚度）.

图3给出了三种胶结物分布模型，A1~A8为孔隙分布模型，B1~B8为颗粒包壳模型，C1~C8为喉道分布模型.模型中灰色为岩石颗粒，黑色为孔隙，白色为胶结物.模型的面孔率在1.23%到16.65%之间，从A1~A8，B1~B8，C1~C8孔隙度依次减小.从图中可以看出，对于孔隙分布模型而言，当泥质含量较低时（A1~A4），主要为分散泥质；而当泥质含量较高时（A5~A8），主要为骨架结构泥质.对于喉道分布模型而言，主要为颗粒‒颗粒接触面结构泥质，但在致密砂岩中，由于颗粒与颗粒间为线接触，因此只能代表部分颗粒‒颗粒接触面结构泥质.颗粒包壳模型既包括分散泥质，也包括颗粒‒颗粒接触面结构泥质，其中B1主要为分散泥质，随着泥质包壳的增厚，分散泥质的比例增加.

在致密砂岩中，由于颗粒‒颗粒之间主要为线接触，因此颗粒‒颗粒接触面结构泥质实际上难以直接模拟，只能模拟喉道内的一部分.为了模拟上述泥质分布，我们采用了三维分水岭算法.分水岭（Watershed）是地形学中的一个经典概念（黎蔚和高璐，2013），在图像处理领域中，灰度图像可以被看作地形表面，图像中每个像素的灰度值代表该点的高度，每一个局部极小值（Local minimum）及其影响区域称为集水盆地（Catchments basin），一般将颗粒‒颗粒接触面位置定义为分水岭.分水岭算法通常被用于图像的分割，进而统计对象的数量、形状等参数.Lin et al.（2002）将分水岭算法用于煤岩颗粒的分割，Sheppard et al.（2004）将其应用于分割岩石颗粒或孔隙，进而分析岩石颗粒（孔隙）的大小及形状.上述研究主要应用于颗粒的分割，而如果直接提取颗粒的分水岭，则可以直接用于模拟颗粒‒颗粒接触面结构泥质.

利用分水岭算法构建颗粒接触面结构泥质的方法如下：（1）根据图2b中的二值化图像，计算颗粒和孔隙之间的距离图（图4a）.图中颜色越白代表颗粒与孔隙的距离越大，可以将距离最大的位置认为是山峰的位置；（2）利用H⁃maxima方法（Liu et al.， 2017）对不同的局部高点进行标记（图4b），两个颗粒之间如果进行了标记，则可以认为是两个颗粒，如果没有标记，则认为是一个颗粒.局部高点标记后对局部高点进行分割，不同的颜色代表不同的高点（图4c）；（3）利用基于高点标记的分水岭算法进行岩石颗粒的分割，图4d中蓝色的分割线即为分水岭的位置.（4）将获取的分水岭图像与相应的二值化图像并集即可得到颗粒‒颗粒接触面泥质（图4e），将图4e与二值化图像相加即可得到包含颗粒‒颗粒接触面泥质的模型（图4f）.此时，分水岭位置的像素控制在1个像素，可以通过改变分水岭处像素的多少，调整颗粒和颗粒接触面泥质的含量.

对比图4f和图2a可以看出，除了部分位置存在颗粒的过度分割，分水岭算法总体能较好地分割颗粒，模拟颗粒‒颗粒接触面泥质，颗粒接触面泥质的含量为4.24%.

对于致密砂岩而言，微孔隙的发育特征实际上与泥质有密切的关系，为了模拟微孔隙发育特征的影响，我们通过改变泥质的弹性参数来模拟微孔隙的影响.

1.2　模型参数选取

模型仅考虑石英和泥质两种介质，为了模拟微孔隙发育特征的影响，采用了四种不同的泥质弹性参数.泥质弹性参数的选取参考Sams and Andrea（2001）给出的垂直于层理方向黏土矿物纵横波速度与孔隙度的关系，弹性参数分别为孔隙度0%，2.5%，5%和10%时的取值，具体取值见表1.同时，给出了Xu and White（1995）以及Mavko et al.（2009）给出的泥质弹性参数，可以看出Xu and White（1995）的参数更接近泥质孔隙度为0%时的参数，Mavko et al.（2009）的泥质参数更接近泥质孔隙度为5%时的弹性参数.

1.3　模拟方法

声波特征的模拟采用有限元软件Comsol Multiphysics，利用固定力学模块可以模拟岩石的线弹性特征，根据弹性参数即可计算岩石的声波速度.图5给出了计算纵波模量（Shulakova et al.，2013）和剪切模量的边界条件.在计算纵波模量时，设置图像底部固定，两侧在x方向上位移为零，y方向可自由移动.添加顶部的正应力（形变），则纵波模量M=σ_y/ε_y .其中，σ_y 为正应力，ε_y 为正应变.在模拟剪切应变时，固定底部，顶部添加x方向应力，剪切模量为μ=σ_x/γ.其中σ_x 为剪应力，γ为剪应变.体积模量表示为：

K=M-4μ/3,（1）

纵横波速度分别为：

V_{\mathrm{p}}=\sqrt[]{\frac{M}{\rho }},（2）

V_{\mathrm{s}}=\sqrt[]{\frac{\mu }{\rho }},（3）

其中V_p 和V_s 分别为纵横波速度，km/s，ρ为岩石密度，g/cm³.

2 模拟结果分析

研究首先模拟了不同胶结物分布模式的岩石弹性参数变化规律.其中除了孔隙分布模型、包壳模型和喉道分布模型，还包括包含颗粒‒颗粒接触面泥质的喉道分布模型.图6给出了不同胶结物分布形式的体积模量和剪切模量与孔隙度的关系，泥质弹性参数选表1中泥质5的取值，图中结点对应于图3中A1~A8，B1~B8和C1~C8的图像有限元模拟结果.结果显示，不同泥质分布模型，岩石的弹性参数差异较大.从图6中可以看出，喉道分布模型以及含颗粒‒颗粒接触面泥质的喉道分布模型体积模量均大于孔隙分布模型和颗粒包壳模型.孔隙分布模型与颗粒包壳模型在相同孔隙度条件下体积模量相当，后者略高于前者.

2.1　泥质分布形式对岩石物理性质的影响

从泥质的分布形式来看，以分散泥质为主的A1~A4和B1比与其孔隙度相当的以颗粒‒颗粒接触面泥质为主的C1~C4低10 GPa左右，说明颗粒‒颗粒接触面泥质能明显增加岩石的体积模量，分散泥质增加对岩石体积模量影响较小.以骨架结构泥质为主的A5到A7的体积模量也小于泥质含量相当的以颗粒‒颗粒接触面泥质的C5到C7的体积模量，但其随孔隙度减小的增加速度大于A1~A4，说明骨架结构泥质对岩石体积模量的影响也较小，但强于分散泥质.对比喉道分布模型和含颗粒‒颗粒接触面泥质的喉道分布模型，后者略低于前者，且曲线几乎平行，说明颗粒‒颗粒接触面泥质与喉道分布泥质对体积模量的影响相当.同时，如果仅分析A1到A4以及C1到C4体积模量的变化，二者体积模量变化差值为4.6 GPa和5.4 GPa.并且从C5到C8体积模量随孔隙度减小而增加的速度减小，说明泥质含量对体积模量的影响小于泥质的分布形式.当泥质含量很高时，包括A8，B6~B8以及C8，孔隙多被泥质充填，不同泥质分布形式体积模量趋于一致.

对比孔隙分布模型、颗粒包壳模型以及喉道分布模型，三种类型剪切模量相当，说明泥质的分布形式对剪切模量几乎没有影响，剪切模量主要受泥质含量控制.同时，含颗粒‒颗粒接触面泥质的喉道分布模型体积模量低于其他三种模式.其中，喉道分布模型和含颗粒‒颗粒接触面泥质的喉道分布模型的泥质含量差别为4.24%，剪切模量差值为 2.38 GPa，体积模量差值为1.13 GPa，说明颗粒‒颗粒接触面泥质对剪切模量的影响大于体积模量.

2.2　泥质弹性参数差异对岩石物理性质的影响

除了模拟泥质分布形式的影响，还模拟了泥质的弹性参数的差异对岩石弹性参数的影响，能够间接分析微孔隙对岩石弹性参数的影响.图7分别为孔隙分布模型、颗粒包壳模型、喉道分布模型以及含颗粒‒颗粒接触面泥质的喉道分布模型的弹性参数与孔隙度的关系，且考虑不同的泥质弹性参数.以分散泥质为主的A1~A4和B1，不同的泥质弹性参数条件下体积模量均相同，说明与分散泥质相关的微孔隙对岩石的弹性参数影响小.与骨架结构泥质和颗粒‒颗粒接触面泥质相关的微孔隙发育则会略微减小体积模量，而对剪切模量影响较大.

通过模拟可以得到以下规律：（1）泥质的分布形式对岩石体积模量影响较大，而对剪切模量影响较小，其中分散泥质和骨架泥质对岩石体积模量影响较小，喉道分布泥质和颗粒‒颗粒接触面泥质则会明显增加体积模量.（2）泥质含量对岩石体积模量的影响小于泥质分布形式，而泥质含量对剪切模量的影响大于泥质的分布形式.（3）与分散泥质相关的微孔隙对岩石弹性参数没有影响，而骨架结构泥质和颗粒‒颗粒接触面泥质则对微孔隙的发育敏感，且微孔隙增加对剪切模量的减小大于体积模量.

3 三维数字岩心预测声波速度

3.1　三维数字岩心模型构建

数值模拟结果显示，泥质的分布位置对岩石弹性参数有较大影响，然而，CT图像无法区分泥质和砂岩颗粒，因此利用CT图像模拟含泥砂岩的声波特征存在较大的难度.Shulakova et al. （2013）提出首先利用自洽模型计算颗粒和泥质构成的固体介质的弹性参数，其中泥质含量通过SEM图像确定，该方法没有考虑泥质分布形式的影响.

根据数值模拟结果实际上可以将泥质的分布形式分为两类，一类为分散泥质和骨架结构泥质，二者对岩石弹性特征影响较小；另一类为喉道分布泥质和颗粒‒颗粒接触面泥质，对岩石弹性影响较大.对于第一类，可以直接将其作为岩石颗粒处理，第二类可以利用分水岭算法确定颗粒‒颗粒接触面泥质.为了验证该思路的可行性，从鄂尔多斯盆地定北地区太原组选取了3块致密砂岩样品（D20，D105，D6）开展了micro⁃CT扫描，同时对样品制取了铸体薄片，利用AP⁃608覆压孔渗仪开展了孔隙度和渗透率测试.

三个样品的孔隙度分别为8.36%，7.73%和5.54%，空气渗透率分别为2.439 mD，0.058 mD和0.041 mD.micro⁃CT扫描采用通用公司的Nanotom纳米CT，分辨率分别为1.518 μm，0.9 μm以及0.875 μm.从铸体薄片可以看出（图8），D20样品仅含少量分散泥质，泥质含量为0.8%，D105和D6样品泥质含量较高，分别为19%和18%，孔隙和喉道几乎被泥质完全充填，D105样品存在钙质胶结物的溶蚀.从与之对应的CT切片显示，D20样品为纯净砂岩，D105和D6样品颗粒间仍发育大量微孔隙，特别是泥质并未完全充填粒间孔，泥质主要呈分散泥质形式，而不是骨架结构泥质.

利用CT图像建立数字岩心的过程依次为图像滤波、图像二值化、表征单元体确定.本研究图像滤波采用各向异性扩散滤波（Soille，2003），图像二值化采用直方图阈值法（姜黎明，2012），表征单元体选取采用从小到大取立方体，逐渐扩大，并计算不同边长立方体的孔隙度，作出孔隙度随数字岩心图像边长变化规律曲线，同时改变立方体的中心，重复采样，最终所有孔隙度随数字岩心图像边长变化规律曲线趋于稳定的最小岩心边长即为表征单元体尺寸（姜黎明，2012）.对三个样品的分析发现（图8），表征单元体为 1 100×1 100×1 100，利用数字岩心计算的三个样品的孔隙度分别为8.9%、5.2%和4.3%，可以看出纯净砂岩数字岩心孔隙度略大于实测孔隙度，可能受非均质性的影响.含泥质砂岩数字岩心孔隙度小于实测孔隙度，样品D105由于钙质胶结溶蚀孔无法识别，与实测孔隙度的差值（2.53%）大于样品D6与实测孔隙度的差值（1.24%）.

分水岭算法主要用于确定颗粒‒颗粒接触面泥质的位置，但无法确定颗粒‒颗粒接触面泥质的含量.从图4中可以看出，分水岭分割界限一般为1个像素，少量位置为多个像素，而泥质含量对岩石的弹性性质有较大影响.从铸体薄片图8B1和图8C1可以看出，颗粒‒颗粒间泥质厚度变化较大，图中标出的两处分别为12 μm和21 μm，部分颗粒‒颗粒直接接触.从CT扫描样品图8C2中则可以看出，颗粒‒颗粒接触面位置厚度相当，标出位置为13 μm.由于CT得到的颗粒‒颗粒接触面泥质厚度较为一致，且D105样品和D6样品泥质含量相当，那么可以在数字岩心计算中首先尝试利用相同的颗粒‒颗粒接触面泥质厚度.

图9为样品D105图像处理的过程，图像处理分为4步：（1）对原始图像进行滤波、重采样，重采样之后立方体像素边长为3.6 μm，然后二值化处理，删除孤立的小孔隙，并做平滑处理（图9a，9b）；（2）利用基于标记的分水岭算法得到颗粒‒颗粒接触面泥质，分水岭位置一般为1个像素，少量位置存在多个像素（图9c）；（3）对分水岭位置的像素进行膨胀运算（刘学锋，2010），即可得到宽度为三个像素的颗粒‒颗粒接触面泥质（图9d）.

D20样品由于泥质含量极低，直接进行第一步操作即可，而对于样品D105和D6，则需要执行上述3步操作.对于D105和D6样品，颗粒‒颗粒接触面泥质分别为10.8 μm和10.5 μm，符合CT图像中颗粒‒颗粒接触面泥质的宽度.构建的三维模型见图10，其中黄色为颗粒，绿色为泥质，D105和D6的泥质含量分别为17.1%和15.0%.

除了颗粒‒颗粒接触面泥质的宽度，其岩石力学参数也是影响岩石声学特征的重要参数.D105和D6数字岩心孔隙度与实测孔隙度的差值分别为2.53%和1.24%，可以认为该孔隙度主要为颗粒‒颗粒间泥质孔隙度，其中D105也包括部分方解石溶蚀孔.如果假设实测孔隙度与数字岩心孔隙度的差值为泥质孔隙度，颗粒‒颗粒间泥质的孔隙度则分别为14.7%和8.3%，而样品D105由于存在方解石的溶蚀孔，其实际孔隙度应低于14.7%，为了方便，取泥质孔隙度为10%时的岩石力学参数作为颗粒‒颗粒间泥质的力学参数.

3.2　岩心纵横波速度测试

岩心纵横波速度测试除了上述3块样品，还另外选取了53块不同孔隙度的样品，这些样品测试可以反映研究区纵横波速度测试结果的趋势，检验3块micro⁃CT扫描样品纵横波测量的准确性.

实验采用目前常用的脉冲穿透法，该方法一般在柱样两端放置两个声波探头，其中一个为激发探头，另一个为接收探头，岩样和声波探头之间需要添加耦合剂（选用凡士林），然后通过示波器等记录纵横波通过岩样的时间，根据岩样长度和声波穿过的时间，即可计算纵横波速度.测试的流程如下：

（1）将2.5 cm直径的岩心切成4~6 cm的柱样，样品的长度尽量长，保证较长的声波传播时间，减小实验误差.在切割过程中，保证样品端面与样品长轴垂直，以确保样品与探头之间的耦合.为了保证样品端面与样品长轴垂直，样品测3次长度，3次长度测量的误差在0.05 mm以内.

（2）对岩心进行洗油、洗盐，然后在60 ℃温度条件下烘干24 h，确保岩心的干燥.

（3）在样品的端部涂凡士林作为耦合剂，先加轴压到5 MPa，再加围压到5 MPa，然后同时加围压和轴压到35 MPa，内压为大气压.

（4）首先开展纵波测试，根据纵波接收波形，确定首波到达时间，根据样品的长度计算纵波速度.然后进行横波测试，根据横波波形计算首波达到时间，计算纵波速度.在测试过程中，需要保证纵波速度大于横波速度，且纵横波速度比一般要大于1.5，以保证测量结果的准确性.

3.3　模型验证及讨论

利用二维数字岩心构建中提到的方法对3块数字岩心进行有限元模拟，图11为D20样品三维数字岩心（a）挤压和（b）剪切条件下米塞斯应力分布.图11a挤压过程中，保持数字岩心底部固定，侧面仅允许Z轴方向变形，顶部施加向下的压应力.根据垂向施加的应力和相应的应变即可计算纵波模量.图11b剪切模拟中，保持数字岩心底部固定，侧面仅允许X轴方向变形，顶部施加X轴方向的剪切应力.根据水平的剪应力和剪应变，即可计算剪切模量.根据纵波模量和剪切模量，利用公式（2）和公式（3）分别计算纵横波速度，即可与实测岩石纵横波数据进行对比.

从图12的模拟结果可以看出，3个样品的纵横波速度与实测值均较为吻合，其中D20的数字岩心模拟纵横波速度与实测值吻合最好，其他2个样品横波模拟值与实测值吻合略差，但对比其他样品的横波测试结果，总体仍在较窄的范围内变化，说明利用数字岩心技术能够较为精确模拟岩石纵横波速度.同时，对比D20和D105，二者孔隙度相当，但由于D105泥质含量较高，横波速度小于前者，说明泥质对于声波速度的影响极为重要，考虑颗粒‒颗粒接触面泥质的模型总体有较高的精度.

虽然利用数字岩心技术能够较好地模拟岩石纵横波速度，但需要注意的是数字岩心技术目前主要用于模拟高有效压力条件下的岩石纵横波速度.大量的岩石纵横波速度测试结果显示，岩石的纵横波速度与压力关系密切，在低有效压力条件下，声波速度随有效压力变化极大，而在高有效压力下，声波速度随有效压力变化不大.Smith et al.（2009）对致密砂岩的实验显示，当有效压力大于28 MPa时，岩石的纵横波速度随岩石有效压力增加几乎不变，并指出裂缝型孔隙的闭合是导致纵横波速度随有效压力增加的主要因素.Arns et al.（2002）对枫丹白露砂岩开展数值模拟，与Han et al.（1986）的实验数据进行对比，其实验测试的有效压力为40 MPa.因此，数字岩心主要用于模拟高有效压力条件下的纵横波速度.

同时，该方法主要用于模拟较高泥质含量的致密砂岩.从D105和D6两个样品的分析来看，颗粒‒颗粒接触面泥质含量分别为17.1%和15.0%.对于颗粒‒颗粒接触面泥质含量较低的样品，可能会高估泥质含量，需要对方法进行改进.但该方法为利用数字岩心模拟复杂砂岩声波特征提供了建设性思路.

4 结论

（1）泥质的分布形式对岩石体积模量影响较大，而对剪切模量影响较小，其中分散泥质和骨架泥质对岩石体积模量影响较小，喉道分布泥质和颗粒‒颗粒接触面泥质则会明显增加体积模量.

（2）泥质含量对岩石体积模量的影响小于泥质分布形式，而泥质含量对剪切模量的影响大于泥质的分布形式.

（3）与分散泥质相关的微孔隙对岩石弹性参数没有影响，而骨架结构泥质和颗粒‒颗粒接触面泥质则对微孔隙的发育敏感，且微孔隙增加对剪切模量的减小大于体积模量.

（4）考虑颗粒‒颗粒接触面泥质的数字岩心能够较好模拟泥质含量较高致密砂岩的纵横波速度，为模拟复杂砂岩声波特征提出了建设性的思路.

参考文献

[1] Arns, C. H., Knackstedt, M. A., Pinczewski, W. V., et al., 2002. Computation of Linear Elastic Properties from Microtomographic Images: Methodology and Agreement between Theory and Experiment. Geophysics, 67(5): 1396. https://doi.org/10.1190/1.1512785

[2] Cheng, J.W., Zhang, F., Li, X.Y., 2024. Orthorhombic Anisotropic Rock Physics Modeling for Fractured Marine Shale Reservoir in Sichuan Basin. Earth Science, 49(1): 299-312 (in Chinese with English abstract).

[3] Dvorkin, J., Nur, A., 1996. Elasticity of High⁃Porosity Sandstones: Theory for Two North Sea Data Sets. Geophysics, 61(5): 1363-1370. https://doi.org/10.1190/1.1444059

[4] Han, D. H., Nur, A., Morgan, D., 1986. Effects of Porosity and Clay Content on Wave Velocities in Sandstones. Geophysics, 51(11): 2093-2107. https://doi.org/10.1190/1.1442062

[5] Jiang, L.M., 2012. Numerical Simulation of Acoustic and Electrical Characteristics of Natural Gas Reservoir Rocks Based on Digital Cores (Dissertation). China University of Petroleum, Dongying (in Chinese with English abstract).

[6] Jin, G. D., Patzek, T. W., Silin, D. B., 2012. Modeling the Impact of Rock Formation History on the Evolution of Absolute Permeability. Journal of Petroleum Science and Engineering, 100: 153-161. https://doi.org/10.1016/j.petrol.2012.03.005

[7] Li, W., Gao, L., 2013. Pavement Crack Detection Based on Improved Watershed Algorithm. Computer Engineering and Applications, 49(20): 263-266, 270 (in Chinese with English abstract).

[8] Lin, C. L., Miller, J. D., Luttrell, G. H., 2002. Evaluation of a CT⁃Based Coal Washability Analysis System under Simulated On⁃Line Conditions. Mining, Metallurgy & Exploration, 19(1): 9-16. https://doi.org/10.1007/BF03402894

[9] Liu, H. P., Luo, Y., Zhao, Y. C., et al., 2022. Controlling Factors and Predictions of Well⁃Connected Pore⁃Throat Volumes in Tight Oil Sandstones. Journal of Petroleum Science and Engineering, 218: 111034. https://doi.org/10.1016/j.petrol.2022.111034

[10] Liu, H. P., Rodrigues, S., Shi, F. N., et al., 2017. Coal Washability Analysis Using X⁃Ray Tomographic Images for Different Lithotypes. Fuel, 209: 162-171. https://doi.org/10.1016/j.fuel.2017.07.104

[11] Liu, X.F., 2010. Numerical Simulation of Elastic and Electrical Properties of Rock Based on Digital Cores (Dissertation). China University of Petroleum (East China), Qingdao (in Chinese with English abstract).

[12] Mavko, G., Mukerji, T., Dvorkin, J., 2009. The Rock Physics Handbook (2nd Ed). Cambridge University Press, Cambridge. https://doi.org/10.1017/cbo9780511626753

[13] Minear, J. W., 1982. Clay Models and Acoustic Velocities. 57th Annual Meeting, American Institute of Mining, Metallurgical and Petroleum Engineers, Dallas, https://doi.org/10.2118/11031⁃MS

[14] Nur, A., Simmons, G., 1969. The Effect of Saturation on Velocity in Low Porosity Rocks. Earth and Planetary Science Letters, 7(2): 183-193. https://doi.org/10.1016/0012⁃821X(69)90035⁃1

[15] Qiao, H.Q., Liu, C., Fang, H., et al., 2023. S⁃Wave Velocity Prediction Method of Volcanic Rock Based on Statistical Rock⁃Physics Model. Earth Science, 48(8): 2993-3006 (in Chinese with English abstract).

[16] Sams, M. S., Andrea, M., 2001. The Effect of Clay Distribution on the Elastic Properties of Sandstones. Geophysical Prospecting, 49(1): 128-150. https://doi.org/10.1046/j.1365⁃2478.2001.00230.x

[17] Saxena, N., Mavko, G., 2016. Estimating Elastic Moduli of Rocks from Thin Sections: Digital Rock Study of 3D Properties from 2D Images. Computers & Geosciences, 88: 9-21. https://doi.org/10.1016/j.cageo.2015.12.008

[18] Sheppard, A. P., Sok, R. M., Averdunk, H., 2004. Techniques for Image Enhancement and Segmentation of Tomographic Images of Porous Materials. Physica A: Statistical Mechanics and Its Applications, 339(1/2): 145-151. https://doi.org/10.1016/j.physa.2004.03.057

[19] Shulakova, V., Pervukhina, M., Müller, T. M., et al., 2013. Computational Elastic Up⁃Scaling of Sandstone on the Basis of X⁃Ray Micro⁃Tomographic Images. Geophysical Prospecting, 61(2): 287-301. https://doi.org/10.1111/j.1365⁃2478.2012.01082.x

[20] Smith, T. M., Sayers, C. M., Sondergeld, C. H., 2009. Rock Properties in Low⁃Porosity/Low⁃Permeability Sandstones. The Leading Edge, 28(1): 48-59. https://doi.org/10.1190/1.3064146

[21] Soille, P., 2003. Morphological Image Analysis-Principles and Applications. Springer, New York.

[22] Swanson, B. F., 1985. Microporosity in Reservoir Rocks-Its Measurement and Influence on Electrical Resistivity. The Log Analyst, 26(6): 42-52.

[23] Torskaya, T. S., 2013. Pore⁃Scale Analysis of Grain Shape and Sorting Effects on Fluid Transport Phenomena in Porous Media (Dissertation). The University of Texas at Austin, Austin.

[24] Wang, K.W., 2007. Reservoir Rock Electrical Properties Study Based on Percolation Network Model and Upscaling (Dissertation). China University of Petroleum (East China), Dongying (in Chinese with English abstract).

[25] Wyllie, M. R. J., Gregory, A. R., Gardner, G. H. F., 1958. An Experimental Investigation of Factors Affecting Elastic Wave Velocities in Porous Media. Geophysics, 23(3): 459-493. https://doi.org/10.1190/1.1438493

[26] Xu, M. M., Yin, X. Y., Zong, Z. Y., et al., 2020. Rock⁃Physics Model of Volcanic Rocks, an Example from Junggar Basin of China. Journal of Petroleum Science and Engineering, 195: 107003. https://doi.org/10.1016/j.petrol.2020.107003

[27] Xu, S. Y., White, R. E., 1995. A New Velocity Model for Clay⁃Sand Mixtures 1. Geophysical Prospecting, 43(1): 91-118. https://doi.org/10.1111/j.1365⁃2478.1995.tb00126.x

[28] Zeng, W.C., 1991. Oil and Gas Reservoir Well Logging Evaluation Technology. Petroleum Industry Press, Beijing (in Chinese).

基金资助

国家自然科学基金青年基金项目(41902147)

国家自然科学基金青年基金项目(41402117)