2020, Vol. 63
2. 中国科学院地球科学研究院, 北京 100029;
3. 中国石油大学(华东)深层油气重点实验室, 青岛 266580;
4. 中国科学院大学, 北京 100049
2. Institutions of Earth Science, Chinese Academy of Sciences, Beijing 100029, China;
3. Key Laboratory of Deep Oil and Gas(China University of Petroleum(East China)), Qingdao 266580, China;
4. University of Chinese Academy of Sciences, Beijing 100049, China
页岩中的TOC含量对富有机质页岩弹性性质有重要影响(如, Vernik and Nur, 1992; Carcione et al., 2011; Sayers, 2013; Zhao et al., 2016).研究TOC含量及其分布对页岩弹性模量的影响，对页岩气甜点预测具有重要的指导作用.富有机质页岩主要由矿物骨架、孔隙和有机质构成，一般而言，页岩气甜点介质的孔隙度较低，而TOC含量变化幅度较大，孔隙度与TOC含量的相关性不明显(如, Vernik and Nur, 1992; Vernik and Liu, 1997).根据岩石中TOC的赋存状态，TOC可视为岩石基质中的夹杂物，因此，有关TOC含量对岩石弹性模量的影响可视为多孔介质中的孔隙物质充填问题来研究.孔隙物质置换理论和方法广泛用于研究孔隙物质类型及其含量对岩石弹性性质的影响(如, Ciz and Shapiro, 2007; Das and Batzle, 2008; Makarynska et al., 2010; Saxena and Mavko, 2014).本文采用孔隙物质置换理论方法，研究TOC含量对岩石弹性模量的影响.
Hashin和Shtrikman(1963)首先给出了多矿物组成岩石的有效弹性模量上下限值(称为HS界限)，多孔岩石的有效弹性模量估计是近半个世纪以来的研究热点问题.对于单矿物组成的均匀各向同性饱和流体孔隙岩石，假设孔隙相互连通，孔隙压力处处平衡，可建立联系干岩弹性模量与饱和流体岩石有效弹性模量的Gassmann方程(如, Gassmann, 1951; Brown and Korringa, 1975; Zimmerman et al., 1986; Zimmerman, 1991).据此可计算多孔介质充填流体后的有效弹性模量.Brown和Korringa(1975)将Gassmann方程推广到均匀各向异性岩石，计算均匀各向异性饱和流体岩石的有效弹性模量.Ciz和Shapiro(2007)基于均匀各向同性岩石骨架，假定孔隙内部孔隙压力处处平衡，将Gassmann流体替换方程推广到了固体替换领域，得到Ciz-Shapiro(C&S)方程，并用有限元(FEM)数值模拟验证了该方程的计算精度.Makarynska等(2010)通过真实岩石的重油充填实验发现，C&S方程的计算结果与实验测量之间存在误差，这与岩石骨架非均质性或孔隙分布不均匀有关(Mavko and Saxena, 2013).Saxena和Mavko(2014)基于互易定理，通过引入孔隙应变对岩石弹性性质的影响，建立了适用于非均质岩石的固体置换方程(S&M方程).该方程通过引入描述岩石内部结构非均质性的参数，能够精确计算孔隙充填固体的岩石的有效弹性模量.
以上关于固体置换理论的研究主要针对孔隙度较大的砂岩，随着页岩气勘探技术的日益发展，有关孔隙中充填重油、沥青、TOC等固体或黏弹性材料的岩石弹性性质的研究日益增多(Vernik and Milovac, 2011; Zhao et al., 2016; Hu et al., 2018).由于页岩具有孔隙度低、孔隙形态复杂、非均质性强等特点(如, Vernik and Liu, 1997; Hornby, 1998), 固体替换理论能否适用于页岩，目前还未深入研究.本文通过数字岩心技术提取页岩孔隙和TOC形态以及分布特征，构建了不同TOC含量的岩心模型.采用FEM数值模拟(COMSOL软件)对比分析方法，研究了C&S和S&M两种固体替换方程对页岩的适用性，分析了TOC含量对页岩有效弹性模量以及页岩非均质特征的影响.本文结构如下：首先，简要介绍C&S和S&M两种替换方程，描述岩石物理实验和数字岩心基本情况.然后，通过理论模型，检验两种替换方程的计算精度.最后针对实际岩心和超声实验数据，构建不同TOC含量的岩心弹性模型，研究TOC含量变化对页岩弹性模量的影响.
1 理论方法为了研究多孔介质中的孔隙物质固体替换问题，Ciz和Shapiro(2007)基于固体骨架和孔隙物质的无穷小应变假设，并假设孔隙均匀分布，岩石内部孔隙压力处处平衡，将Gassmann流体替换方程推广到了固体替换领域，得到了如下的固体替换方程(简称C&S方程)：
|
(1) |
|
(2) |
其中Ksat和μsat分别是孔隙充填物质(固体或流体)后的岩石有效体积模量和剪切模量.sdry和μdry分别为干岩石的体积模量和剪切模量，Kif、Kgr、μif和μgr分别为孔隙物质和岩石基质颗粒的体积模量和剪切模量，ϕ为岩石的总孔隙度.通过均匀介质的FEM数值模拟验证，C&S方程在已知岩石骨架体积模量、剪切模量和孔隙度的条件下，可以精确计算孔隙充填固体后，岩石的有效体积模量和剪切模量.
当岩石骨架为非均质和孔隙分布不均时，岩石内部容易发生应力集中，孔隙压力处处平衡的假设不成立，C&S固体替换方程失效.为解决此问题，Saxena和Mavko (2014)基于互易定理，通过引入描述岩石内部结构非均质性的参数α1，α2，β1和β2, 建立了如下计算非均质岩石有效弹性模量的固体替换方程(简称S&M方程)：
|
(3) |
|
(4) |
其中，参数α1(α1≥1)、α2(α2≥0)、β1(β1≥1)和β2(β2≥ 0)与岩石内部微观结构有关，其中α1和α2分别代表孔隙物质的体积模量和剪切模量对岩石有效体积模量的影响，β1和β2分别代表孔隙物质的剪切模量以及体积模量对岩石有效剪切模量的影响.对于内部结构均匀的岩石，有α1=β1=1和α2=β2=0, S&M方程退化为C&S方程.当岩石结构均匀，孔隙物质为流体时(剪切模量为0)，两个替换方程均退化为如下的Gassmann方程：
|
(5) |
|
(6) |
其中α=1-Kdry/Kgr, 是Biot-Willis系数(Biot and Willis, 1957), M是岩石的耦合模量，可表示为
|
(7) |
其中Kf是孔隙流体的体积模量.
通过对比C&S方程与S&M方程计算结果差异，可推测S&M方程中参数α1、α2、β1和β2的变化情况，据此分析岩石非均质和孔隙充填物对岩石有效弹性模量的影响.本文假定孔隙充填不同含量的有机质，研究TOC含量变化对岩石有效弹性模量的影响.采用数值模拟的方法验证C&S方程以及S&M方程对预测TOC含量对页岩影响的有效性.C&S方程以及S&M方程是一种预测岩石有效弹性性质的方法.为突出TOC的影响，排除岩石中其他因素的影响，数值模拟验证被广泛利用.其中Ciz和Shapiro(2007)采用有限元(FEM)模拟研究了对于内部结构均匀的砂岩，C&S方程的有效性.Saxena和Mavko(2014)采用有限元模拟方法，进一步研究了砂岩中，孔隙非均匀分布时S&M方程的有效性.另外，Saenger和Shapiro(2002)通过有限差分验证了等效介质理论的有效性.Guo等(2018)通过有限元方法，验证了裂缝介质理论频散模型的正确性.采用数值模拟方法验证有关等效介质理论的相关理论是合适的，而且被广泛应用.因此，本文采用数值模拟验证C&S方程以及S&M方程的准确性是有效的.
在本研究中，HS界限，作为均匀多矿物组成的岩石弹性模量的上下界限，对于检验固体替换方程的有效性十分重要，因此，直接给出当岩石由两种矿物组成时，HS上下界限.当其结构为坚硬物质包围软弱物质时，取得其上界限.其上界限表达式为
|
(8) |
|
(9) |
当其结构为软弱物质包围坚硬物质时，取得其下界限.其下界限表达式为
|
(10) |
|
(11) |
以上界限，限定了由两种物质组成的各向同性均匀岩石弹性模量的上下界限.
2 实验测量研究TOC含量对页岩弹性性质的影响，需要分析页岩弹性模量数据及其内部结构特征.因此，我们分别进行了不同TOC含量页岩的超声波速测量实验及其数字岩心测量.
2.1 超声速度测量实验装置以及结果本研究通过超声速度测量实验测量页岩的体积模量以及剪切模量.超声速度测量实验研究中，采用的高压超声速度测量实验装置的核心部分结构如图 1所示.在该实验装置中，一个伪三轴加压装置对样品施加围压(Pc).纵、横波换能器(transducer)放置于样品两端，用于激发纵、横超声波.通过测量纵、横超声波的传播时间，结合样品长度，可获得超声波的纵、横波速度(Nur, 1969; Toksöz et al., 1976; Cheng and Toksöz, 1979; 周浩和符力耘，2018).此实验设备中纵、横波换能器主频均为1 MHz.
|
图 1 实验装置示意图 Fig. 1 The experiment setup |
利用上述实验设备，我们对一块直径为2.52 cm长度为4.95 cm的柱体页岩岩心，在稳定的温度(20 ℃室温)和不同的压力条件下，进行了干燥条件以及饱和流体条件的超声速度测量.
在实验测量过程中，岩石的围压(Pc)由5 MPa, 以10 MPa为压力间隔，逐步增大到55 MPa, 同时，分别测量页岩样品的纵、横波速度.岩石的体积模量以及剪切模量与波速之间存在如下关系(Mavko et al., 2009)：
|
(12) |
|
(13) |
其中st=dry代表干岩石的弹性模量与波速之间的关系，当st=sat代表饱和流体岩石的弹性模量与波速之间的关系.ρst为岩石密度，对于干岩石，ρdry为2.48 g·cm-3, 对于饱和流体岩石，ρsat为2.54 g·cm-3.VPst和VSst分别为样品的纵横波速度.经过方程(8)(9)的处理，超声速度测量实验结果如图 2所示.
|
图 2 实验测量弹性模量随围压变化数据 (a)体积模量；(b)剪切模量. Fig. 2 The experimental elastic moduli variation versus confining pressure (a) Bulk modulus; (b) Shear modulus. |
图 2展示了样品在干燥以及饱和水条件下的弹性模量随围压变化曲线.通过观察，我们发现，随着压力增加，页岩的弹性模量增加，与前人得到的趋势类似(如, Dewhurst and Siggins, 2006; Sarout and Guéguen, 2008).但是，通过对图 2a分析可知，在低压状态下岩石的体积模量变化相对平稳，在高压条件下，反而其增长速率较快.根据Shapiro(2003)以及David和Zimmerman(2012)的研究，对于砂岩而言，由于低压阶段裂缝大量闭合，其速度增长较快，在高压阶段，由于裂缝几乎全部闭合，因此速度增长变慢.对于本研究中的页岩样品，在高压阶段，速度增加稍大于低压阶段，说明，在加压过程中，岩石内部结构发生改变，会产生新的裂缝，导致在高压阶段，模量增加较快.因此，根据图 2中的数据，我们可知，随着压力的增加，弹性模量呈现非线性增加.这个现象说明，压力加载过程导致岩石内部结构变化，使得岩石内部非均质性发生改变(如, Fu B Y and Fu L Y, 2017, 2018).
2.2 数字岩心测量为了进一步研究页岩内部的微观结构，我们采用数字岩心技术，通过CT扫描获取页岩的数字岩心数据，提取岩石的内部结构信息.图 3所示数字岩心切片边长以及直径为1 mm，分辨率为1 μm，因此最小可分辨直径为1 μm的孔隙以及TOC.我们采用CT切片的阈值分割技术，根据不同矿物的X射线吸收率不同，对CT切片的灰度值进行阈值分割(阳树洪，2014; Zhang et al., 2016; 王志伟等，2018)，识别不同矿物.在本文中，选取像素点灰度值小于76为孔隙，灰度值在76到81之间为TOC，根据不同像素点所占面积，可预测页岩孔隙体积以及TOC体积，进而计算孔隙度以及TOC含量.通过CT扫描，得到的该样品的孔隙度为3.5%，TOC含量为6.75%.该样品矿物组成性质见表 1.
|
图 3 页岩样品的数字岩心图片 (a)正视图；(b)右视图；(c)俯视图；(d)孔隙形态图. Fig. 3 The digital cores of the shale sample (a) Front vision; (b) Right vision; (c) Top vision; (d) The shape of the pore. |
|
表 1 实验页岩样品性质 Table 1 The properties of the experimental shale sample |
我们得到的每张数字岩心切片分辨率均为1 μm.数字岩心切片如图 3所示.
图 3a—c为同一块岩心的三个方向的数字岩心切片，代表了三个方向的孔隙、矿物分布情况以及非均质特征.其中，黑色斑点为孔隙.图 3d为孔隙放大后的形态.在页岩中，TOC充填于孔隙中，并且页岩各向异性以及非均质性主要由孔隙分布决定.根据之前研究(Zhang et al., 2016; Hu et al., 2018; 王志伟等，2018), 对于本文研究的龙马溪组页岩，其TOC均以长轴状夹杂物的形态赋存于背景岩石中，并且其长轴较小(本研究中，长轴长约为10 μm)，远小于入射波波长(入射波长约为3 mm).根据Hudson(1986)、Galvin和Gurevich(2007, 2009)、Guo等(2018)的研究，对于长轴状的夹杂物对背景介质的影响，当其尺度远小于入射波波长时，均可认为夹杂物是背景介质弹性性质的一阶扰动(Fu et al., 2018)，可作为孔隙结构处理.因此，在本文中，我们将TOC作为孔隙填充物对待是合理的.我们基于数字岩心的阈值分割技术(王志伟等，2018)，同时提取孔隙分布以及TOC分布信息，计算三张切片的孔隙&TOC分布的自相关函数.
从图 4的自相关曲线可知，三张切片的微结构在较小尺度范围内差异极小(相关系数大于0.2)，随着研究尺度的增加，三张切片的孔隙和TOC分布特征差异明显增大，说明页岩中孔隙和TOC分布存在一定的各向异性.一般而言，当两个序列的相关系数小于0.3时，表征二者不相关(吴传生，2016).因此，可认为三张切片的孔隙和TOC结构和分布特征基本一致，各向异性较弱，可近似为各向同性介质.另外，根据研究，实际页岩确实存在其背景介质由于存在黏土矿物等各向异性矿物，其背景介质也存在强烈的各向异性情况(邓继新等，2004；张广智等，2015)，因此，更真实的理论模型应将背景介质视为各向异性介质，并将TOC夹杂物扰动引入其中.然而，根据Guo等(2019)的研究，在夹杂物数量密度较低的时候，夹杂物的扰动可以线性叠加于背景介质之上.因此，一种研究夹杂物扰动的方法为，首先假设背景介质为各向同性介质，研究夹杂物的影响，之后，再叠加于各向异性的背景介质之上.由于本研究关注TOC夹杂物对页岩弹性性质的影响.因此可以将夹杂物的影响独立研究，将背景介质视为各向同性介质.据此，本文采用各向同性介质的固体替换方程研究TOC含量变化对页岩弹性性质的影响.此外，由于岩心钻取原因，该实验只对样品轴向垂直于层理方向的样品进行了超声实验，在研究中，背景介质的弹性模量通过超声速度测量实验给出.
|
图 4 孔隙和TOC分布的自相关曲线 Fig. 4 Autocorrelation function of the distribution of pore and TOC |
根据前人研究(如, Todd and Simmons, 1972; Zimmerman et al., 1986; Shapiro, 2003; David and Zimmerman, 2012; Fu B Y and Fu L Y, 2017, 2018), 岩石内部孔隙度大小及其分布影响岩石的弹性性质.为了验证基于C&S和S&M固体替换方程研究孔隙度变化和孔隙非均匀分布对岩石弹性性质的影响，我们首先利用理论模型检验两个替换方程的预测精度.
3.1 孔隙均匀分布模型图 5所示为骨架结构和孔隙分布均匀的理想介质模型，边长为1 mm, 模型中均匀分布间距0.1 mm的25个圆形孔隙，形成一种孔隙的均匀分布.通过改变孔隙半径实现孔隙度变化来检验两个替换方程的预测精度.
|
图 5 岩石骨架和孔隙分布均匀的理想介质模型 Fig. 5 The theoretical model for homogenous medium with homogenous rock frame and pore distribution |
实验样品的孔隙度为3.5%，其数字岩心图片显示孔隙主要为椭圆形软孔.为确保孔隙固体替换预测的准确性，背景介质颗粒的弹性模量Kgr和μgr采用超声实验测量得到的背景弹性模量.根据David和Zimmerman(2012)的研究，当岩石受到的有效压力大于10 MPa时，岩石中孔隙纵横比小于0.01的软孔隙可认为全部闭合.根据Shapiro(2003)的研究，多孔岩石的弹性模量主要受软孔隙变形的影响，当软孔隙闭合时，超声实验测得的干岩石弹性模量可近似为岩石背景介质弹性模量.对于图 5所示的理想均匀模型，我们假定背景介质颗粒的弹性模量为Kgr=29.78 GPa, μgr=22.30 GPa.一般而言，常作为TOC的干酪根的弹性模量为Kif=2.9 GPa, μif=2.7 GPa(Mavko et al., 2009).因此，假定理想模型中孔隙完全充填TOC，利用有限元(FEM)数值模拟(COMSOL软件)，分别计算孔隙度为1%，10%，30%，50%和70%时，岩石的有效体积模量和剪切模量.具体验证的实施方案如下：
(1) 分别对岩石施加各向同性压力和纯剪切应力，采用FEM模拟计算干岩石(Kif=0 GPa, μif=0 GPa)的体积模量和剪切模量；
(2) 分别利用C&S和S&M固体替换方程，计算充填TOC后的岩石有效体积模量Ksat和剪切模量μsat;
(3) 利用FEM模拟计算岩石充填TOC后的有效体积模量和剪切模量，并与固体替换方程计算结果进行比较.
根据以上方案，两种固体替换方程和FEM数值模拟两种方法计算得到的体积模量和剪切模量随TOC含量的变化关系如图 6所示.从图中可见，三种方法计算的岩石体积模量和剪切模量曲线吻合很好，基本介于HS上下边界值之间，随着TOC含量的增加而降低.计算结果验证了理想模型岩石结构和孔隙分布均匀的假定，满足α1=β1=1和α2=β2=0.在孔隙排列方式和孔隙形态特征不变条件下，只改变孔隙大小，不会影响岩石的非均质性.前人研究(如, Hashin and Shtrikman, 1963; Shapiro, 2003; Ciz and Shapiro, 2007)表明, HS界限为已知岩石骨架物质组成后，岩石理论弹性模量的上下边界.对于均匀孔隙模型而言，弹性模量的理论计算值应分布在此范围内.图 6b所示的剪切模量当孔隙度为70%时，计算值略低于HS下界限，这是因为，HS剪切模量上下界限值是基于岩石无穷大假设计算得到的，有限尺寸岩石模型存在边界，计算结果会产生误差.根据图 6b的数据，当孔隙度为70%时，理论预测的剪切模量与HS界限的相对误差约为7.5%, 绝对误差小于0.3 GPa, 说明误差较小，可以忽略.
|
图 6 固体替换方程和FEM数值模拟计算的理想介质模型弹性模量随TOC含量变化关系 (a)体积模量；(b)剪切模量. Fig. 6 The elastic moduli variation of theoretical medium calculated by substitution equations and FEM simulation versus the fraction of TOC (a) Bulk modulus; (b) Shear modulus. |
有关C&S和S&M方程的相关研究表明，当岩石骨架非均质或孔隙分布非均匀时，C&S方程失效，S&M方程中的参数α1和β1将大于1，α2和β2将大于0.孔隙物质弹性模量对岩石有效弹性模量产生影响(Saxena and Mavko, 2014).为验证此结论，我们修改图 6均匀模型得到图 7所示的孔隙分布非均匀模型.
|
图 7 岩石骨架均匀、但孔隙分布非均匀的理想介质模型 Fig. 7 The theoretical model for heterogeneous medium with homogenous rock fame and heterogeneous pore distribution |
利用两种固体替换方程和FEM数值模拟方法计算图 7所示模型的岩石有效体积模量和剪切模量，图 8比较了孔隙度为1%，10%，30%，50%和70%时，三种方法计算的有效体积模量和剪切模量随TOC含量的变化关系.可见，S&M方程预测结果与有限元模拟结果一致，而C&S方程预测结果存在差异.由于孔隙的非均匀分布，模型均具有非均匀特征，S&M方程中的参数α1和β1均大于1.特别是剪切模量，C&S方程的预测值明显低于HS下界限，预测结果严重失真.因此，对于非均质多孔岩石，S&M固体替换方程有效弹性模量的预测精度明显优于C&S方程，其参数α1，β1，α2，β2直接反映岩石的非均质特征.
|
图 8 三种方法计算的孔隙非均匀分布理想介质模型的弹性模量随TOC含量变化关系 (a)体积模量；(b)剪切模量. Fig. 8 The elastic moduli variation of theoretical heterogeneous medium calculated by three methods versus the fraction of TOC (a) Bulk modulus; (b) Shear modulus. |
研究TOC含量对岩石弹性性质影响的理想实验方案是，选择背景弹性模量相同，但TOC含量不同的页岩组成测试样本群.由于本研究样品数量有限，为了尽量满足以上条件，可利用数字岩心修改技术，构成同一数字岩心不同TOC含量的样本群.富有机质页岩数字岩心中除了TOC所占的空间比例，骨架部分的孔隙可以按照所需比例充填TOC.
我们以图 3a的岩心切片为例，研究该数字岩心中TOC分布的尺度问题.采用指数型自相关函数构建随机介质，其自相关函数方程为
|
(14) |
|
(15) |
其中，Φ(x, y)是相关系数，(x, y)代表岩心中各点的位置坐标.图 9为该岩心切片中孔隙&TOC分布的自相关曲线及其随机介质模拟曲线，可见，当观测尺度大于10 μm时，实际和模拟的自相关曲线存在差异.然而，该差异部分的相关系数均小于0.2，即相关性很差的部分(吴传生，2016).因此，该岩心切片中孔隙&TOC分布主要以小于10 μm尺度为主，指数型自相关函数可表征其分布特征.
|
图 9 图 3a岩心切片中孔隙&TOC分布的自相关曲线及其指数型随机介质模拟结果 Fig. 9 The autocorrelation function of pore&TOC and the fitting result of exponential autocorrelation function of the rock core slice in Fig. 3a |
图 3所示数字岩心中孔隙&TOC总含量为10.25 %，通过孔隙逐步充填TOC数字岩心修补技术，使该岩心的TOC含量从0%逐步增加至10.25%, 形成同一数字岩心不同TOC含量的样本群，便于研究TOC含量变化对页岩有效弹性模量的影响.图 10为TOC含量分别为1%，5%，7%和10.25%的数字岩心切片，其中的白色斑点为TOC.根据前人研究(如, Vernik and Nur, 1992; Vernik and Liu, 1997), 页岩中孔隙度与TOC含量的相关性并不十分明显，研究TOC含量对岩石弹性模量影响时，可以忽略TOC含量的变化对孔隙度的影响.采用干酪根的体积模量和剪切模量近似表征TOC的弹性模量(Mavko et al., 2009), 令Kc=2.9 GPa, μc=2.7 GPa.对于页岩中不能闭合并且不能充填TOC的孔隙，全部作为背景介质处理，软孔隙全部闭合的页岩背景颗粒的弹性模量为Kgr=29.78 GPa, μgr=22.30 GPa.利用FEM数值模拟方法、C&S和S&M固体替换方程计算图 10岩心模型的干页岩和充填TOC后的有效弹性模量，结果如图 11所示.可见，S&M替换方程与FEM模拟结果基本一致，有效弹性模量随着TOC含量的增加而减小.C&S替换方程计算结果与FEM模拟结果存在差异，说明该岩心TOC分布存在非均质性，我们对S&M替换方程中的参数α1，α2，β1，β2进行分析，结果如图 12所示.
|
图 10 人为构建的TOC含量分别为1% (a)，5% (b)，7% (c)和10.25% (d)的数字岩心切片 Fig. 10 The artificial digital shale cores with different TOC fractions, (a) TOC fraction 1%, (b) TOC fraction 5%, (c) TOC fraction 7%, (d) TOC fraction 10.25% |
|
图 11 三种方法计算的数字岩心模型弹性模量随TOC含量变化关系 (a)体积模量；(b)剪切模量. Fig. 11 The elastic moduli variation of the digital core calculated by three methods versus the fraction of TOC (a) Bulk modulus; (b) Shear modulus. |
|
图 12 非均质参数α1，α2，β1，β2随TOC含量变化关系 (a) α1，β1；(b) α2，β2. Fig. 12 The variation of heterogeneous parameters α1, α2, β1, β2 versus content of TOC |
从图 12可知，随着TOC含量增加，参数α1, α2, β1和β2均呈现先增大后减小趋势，说明当TOC含量较小或较大时，页岩结构趋于均质，特别是参数α1和β1分别表征孔隙物质体积模量和剪切模量对岩石有效体积模量和剪切模量的影响.当TOC含量较低时，页岩整体应变导致TOC周围应变集中(如, Kubair and Bhanu-Chandar, 2008; Fu B Y and Fu L Y, 2018)，TOC的变形对页岩的有效弹性性质影响加强，参数α1和β1随着TOC含量的增加而增大.当TOC含量提高，TOC分布趋于均匀，其变形对页岩的有效弹性性质影响减弱，参数α1和β1减小.同理，我们可以观察到参数α2和β2随TOC含量的先增大后减小的过程.
4.2 压力对孔隙物质充填过程的影响一般而言，在加压超声实验中，干燥页岩和饱水页岩的弹性模量都会随着有效压力的增加而增加，但二者变化幅度存在明显差异.因此，有必要研究压力对孔隙物质充填过程的影响.如表 1所示，本文研究样品的主要矿物为石英，可采用石英为背景介质，其弹性模量为Kgr=37 GPa和μgr=44 GPa(Mavko et al., 2009).本文研究以超声实验得到的干燥页岩弹性模量作为孔隙充填过程中的干页岩弹性模量，分别用C&S和S&M替换方程计算充填后页岩的有效弹性模量，并与该样品超声实验测量值进行比较，结果如图 13所示.
|
图 13 干燥和饱水页岩弹性模量随有效压力变化：物质充填计算值与超声实验结果对比 (a)体积模量；(b)剪切模量. Fig. 13 The elastic moduli variation versus effective pressure for both dry and saturated shales: the comparison between substitution function result and experimental result (a) Bulk modulus; (b) Shear modulus. |
由图 13可知，C&S替换方程对充填后体积模量的预测值小于超声实验测量值；而对于充填后剪切模量的预测值与干燥页岩的剪切模量基本一致，大于超声实验测量值.这说明由于页岩内部结构的强非均质性，C&S方程不能预测页岩中孔隙物质充填后的弹性模量变化.S&M替换方程通过调节参数α1，α2，β1和β2，预测结果与超声实验测量值基本一致，这些参数随压力的变化如图 14所示.
|
图 14 S&M替换方程非均质参数随有效压力变化曲线 (a)参数α1，β1；(b)参数α2，β2. Fig. 14 The heterogeneous parameters of S&M substitution function variation versus effective pressure (a) α1, β1; (b) α2, β2. |
由图 14a可知，参数α1随压力的增加而增加，由于该参数表征非均质性对体积模量的影响，说明在各向同性压力作用下，孔隙形态的变化(如, Vernik and Liu, 1997; Hornby, 1998; Vernik and Milovac, 2011; Hu et al., 2018)及其周围的应力发生集中(如, Kubair and Bhanu-Chandar, 2008; Fu B Y and Fu L Y, 2018), 导致岩石内部非均质性发生较大变化.特别是充填物质受压后，其体积模量发生改变，从而影响页岩的有效体积模量.在图 14b中，参数α2为零，β1保持稳定，说明充填流体不存在剪切模量，对岩石剪切性质无影响.β2表征充填物质体积模量对岩石剪切模量的影响，其值为负，说明充填物质对岩石骨架产生了软化作用(Saxena and Mavko, 2014).根据前人研究(Gregory, 1976; Vernik and Liu, 1997; 邓继新等，2004; Kwon et al., 2004)，由于页岩样品中均含有黏土矿物，尤其是蒙脱石与水发生相互作用而膨胀, 致使页岩发生强烈的化学软化作用，削弱岩石骨架的弹性模量，使岩石整体的弹性模量降低.另外，由于加入孔隙流体之后的蒙脱石膨胀作用，将使得各向异性增加(Deng et al., 2009), 对地震波的衰减作用增加.说明页岩内部结构以及非均质特征发生强烈改变，因此β2值为负值，可以指示页岩内部非均质性的强烈变化.
5 结论和讨论本文利用数字岩心技术，构成同一数字岩心不同TOC含量的样本群，基于C&S和S&M两种固体置换理论模型，通过FEM数值模拟交叉验证，详细研究了TOC含量对页岩弹性性质的影响.骨架介质均匀的理想岩石模型测试表明，当孔隙分布均匀时，C&S和S&M两种替换方程计算的孔隙内充填TOC后的岩石弹性模量与数值模拟结果吻合；当孔隙分布非均匀时，C&S固体替换方程失效，而S&M固体替换方程预测值与数值模拟结果吻合，说明S&M固体替换方程能够较为准确地描述非均质岩石的固体充填过程，两方程的计算差异可表征岩石非均质性的存在.不同TOC含量数字岩心样本测试表明，鉴于TOC的弹性性质相对于背景介质较软，TOC含量的增加一般会导致页岩弹性模量降低.由于实际岩心孔隙及TOC分布的非均质性，C&S替换方程弹性模量预测值与FEM模拟结果存在差异，而S&M替换方程预测值与FEM模拟结果基本一致，利用该方程中的参数α1, α2, β1和β2可详细分析实际岩心孔隙及TOC分布的非均质性.研究表明C&S和S&M两种方程预测差异随着TOC含量的增加呈现先增大后减小的趋势.当TOC含量较低时，TOC作为岩石内部的夹杂物，使得岩石的非均质性增加；随着TOC含量增加到一定量后，TOC逐渐成为岩石的主要矿物，岩石矿物成分趋于单一，结构趋于均匀.压力对孔隙物质充填过程影响的研究表明，C&S方程不能预测压力条件下页岩孔隙物质充填过程的弹性模量变化，说明在各向同性压力作用下，孔隙形态的变化及其周围的应力发生集中，导致岩石内部非均质性发生较大变化.特别是充填物质受压后，其体积模量发生改变，从而影响页岩的有效体积模量.S&M替换方程通过调节参数α1, α2, β1和β2, 预测结果与超声实验测量值基本一致，特别是参数β2，表征充填物质体积模量对岩石剪切模量的影响，体现充填物质对岩石骨架产生了软化作用，削弱岩石骨架的弹性模量，使岩石整体的弹性模量降低.
致谢 感谢审稿专家对本文提出的建设性的修改意见.
Biot M A, Willis D G. 1957. The elastic coefficients of the theory of consolidation. Journal of Applied Mechanics, 24: 594-601. |
Brown R J S, Korringa J. 1975. On the dependence of the elastic properties of a porous rock on the compressibility of the pore fluid. Geophysics, 40(4): 608-616. |
Carcione J M, Helle H B, Avseth P. 2011. Source-rock seismic-velocity models:Gassmann versus Backus. Geophysics, 76(5): N37-N45. DOI:10.1190/geo2010-0258.1 |
Cheng C H, Toksöz M N. 1979. Inversion of seismic velocities for the pore aspect ratio spectrum of a rock. Journal of Geophysical Research: Solid Earth, 84(B13): 7533-7543. DOI:10.1029/JB084iB13p07533 |
Ciz R, Shapiro S A. 2007. Generalization of Gassmann equations for porous media saturated with a solid material. Geophysics, 72(6): A75-A79. DOI:10.1190/1.2772400 |
Das A, Batzle M. 2008. Modeling studies of heavy oil-in between solid and fluid properties. The Leading Edge, 27(9): 1116-1123. DOI:10.1190/1.2978973 |
David E C, Zimmerman R W. 2012. Pore structure model for elastic wave velocities in fluid-saturated sandstones. Journal of Geophysical Research: Solid Earth, 117(B7): B07210. DOI:10.1029/2012JB009195 |
Deng J X, Shi G, Liu R X, et al. 2004. Analysis of the velocity anisotropy and its affection factors in shale and mudstone. Chinese Journal of Geophysics (in Chinese), 47(5): 862-868. |
Deng J X, Wang S X, Han D H. 2009. The velocity and attenuation anisotropy of shale at ultrasonic frequency. Journal of Geophysics and Engineering, 6(3): 269-278. |
Dewhurst D N, Siggins A F. 2006. Impact of fabric, microcracks and stress field on shale anisotropy. Geophysical Journal International, 165(1): 135-148. DOI:10.1111/j.1365-246X.2006.02834.x |
Fu B Y, Fu L Y. 2017. Poro-acoustoelastic constants based on Padé approximation. The Journal of the Acoustical Society of America, 142(5): 2890-2904. DOI:10.1121/1.5009459 |
Fu B Y, Fu L Y. 2018. Poro-acoustoelasticity with compliant pores for fluid-saturated rocks. Geophysics, 83(3): WC1-WC14. |
Fu B Y, Guo J X, Fu L Y, et al. 2018. Seismic dispersion and attenuation in saturated porous rock with aligned slit cracks. Journal of Geophysical Research: Solid Earth, 123(8): 6890-6910. |
Galvin R J, Gurevich B. 2007. Scattering of a longitudinal wave by a circular crack in a fluid-saturated porous medium. International Journal of Solids and Structures, 44(22-23): 7389-7398. DOI:10.1016/j.ijsolstr.2007.04.011 |
Galvin R J, Gurevich B. 2009. Effective properties of a poroelastic medium containing a distribution of aligned cracks. Journal of Geophysical Research: Solid Earth, 114(B7): B07305. DOI:10.1029/2008JB006032 |
Gassmann F. 1951. Uber die Elastizitat poroser Medien. Vierteljahrsschrift der Naturforschenden Gesellschaft in Zurich, 96: 1-23. |
Gregory A R. 1976. Fluid saturation effects on dynamic elastic properties of sedimentary rocks. Geophysics, 41(5): 895-921. DOI:10.1190/1.1440671 |
Guo J X, Rubino J G, Barbosa N D, et al. 2018. Seismic dispersion and attenuation in saturated porous rocks with aligned fractures of finite thickness:Theory and numerical simulations-Part 1:P-wave perpendicular to the fracture plane. Geophysics, 83(1): WA49-WA62. |
Guo J X, Han T C, Fu L Y, et al. 2019. Effective elastic properties of rocks with transversely isotropic background permeated by aligned penny-shaped cracks. Journal of Geophysical Research: Solid Earth, 124(1): 400-424. DOI:10.1029/2018JB016412 |
Hashin Z, Shtrikman S. 1963. A variational approach to the theory of the elastic behaviour of multiphase materials. Journal of the Mechanics and Physics of Solids, 11(2): 127-140. DOI:10.1016/0022-5096(63)90060-7 |
Hornby B E. 1998. Experimental laboratory determination of the dynamic elastic properties of wet, drained shales. Journal of Geophysical Research: Solid Earth, 103(B12): 29945-29964. DOI:10.1029/97JB02380 |
Hu J H, Fu L Y, Wei W, et al. 2018. Stress-associated intrinsic and scattering attenuation from laboratory ultrasonic measurements on shales. Pure and Applied Geophysics, 175(3): 929-962. DOI:10.1007/s00024-017-1705-9 |
Hudson J A. 1986. A higher order approximation to the wave propagation constants for a cracked solid. Geophysical Journal International, 87(1): 265-274. DOI:10.1111/j.1365-246X.1986.tb04556.x |
Kubair D V, Bhanu-Chandar B. 2008. Stress concentration factor due to a circular hole in functionally graded panels under uniaxial tension. International Journal of Mechanical Sciences, 50(4): 732-742. DOI:10.1016/j.ijmecsci.2007.11.009 |
Kwon O, Kronenberg A K, Gangi A F, et al. 2004. Permeability of illite-bearing shale:1. Anisotropy and effects of clay content and loading. Journal of Geophysical Research: Solid Earth, 109(B10): B10205. DOI:10.1029/2004JB003052 |
Makarynska D, Gurevich B, Behura J, et al. 2010. Fluid substitution in rocks saturated with viscoelastic fluids. Geophysics, 75(2): E115-E122. |
Mavko G, Mukerji T, Dvorkin J. 2009. The Rock Physics Handbook:Tools for Seismic Analysis of Porous Media. Cambridge: Cambridge University Press.
|
Mavko G, Saxena N. 2013. Embedded-bound method for estimating the change in bulk modulus under either fluid or solid substitution. Geophysics, 78(5): L87-L99. DOI:10.1190/geo2013-0074.1 |
Nur A M. 1969. Effects of stress and fluid inclusions on wave propagation in rock[Ph. D. thesis]. Cambridge: Massachusetts Institute of Technology.
|
Saenger E H, Shapiro S A. 2002. Effective velocities in fractured media:a numerical study using the rotated staggered finite-difference grid. Geophysical Prospecting, 50(2): 183-194. DOI:10.1046/j.1365-2478.2002.00309.x |
Sarout J, Guéguen Y. 2008. Anisotropy of elastic wave velocities in deformed shales:Part 1-Experimental results. Geophysics, 73(5): D75-D89. DOI:10.1190/1.2952744 |
Saxena N, Mavko G. 2014. Exact equations for fluid and solid substitution. Geophysics, 79(3): L21-L32. DOI:10.1190/geo2013-0187.1 |
Sayers C M. 2013. The effect of kerogen on the elastic anisotropy of organic-rich shales. Geophysics, 78(2): D65-D74. |
Shapiro S A. 2003. Elastic piezosensitivity of porous and fractured rocks. Geophysics, 68(2): 482-486. |
Todd T, Simmons G. 1972. Effect of pore pressure on the velocity of compressional waves in low-porosity rocks. Journal of Geophysical Research, 77(20): 3731-3743. DOI:10.1029/JB077i020p03731 |
Toksöz M N, Cheng C H, Timur A. 1976. Velocities of seismic waves in porous rocks. Geophysics, 41(4): 621-645. DOI:10.1190/1.1440639 |
Vernik L, Nur A. 1992. Ultrasonic velocity and anisotropy of hydrocarbon source rocks. Geophysics, 57(5): 727-735. DOI:10.1190/1.1443286 |
Vernik L, Liu X Z. 1997. Velocity anisotropy in shales:A petrophysical study. Geophysics, 62(2): 521-532. DOI:10.1190/1.1444162 |
Vernik L, Milovac J. 2011. Rock physics of organic shales. The Leading Edge, 30(3): 318-323. DOI:10.1190/1.3567263 |
Wang Z W, Fu L Y, Zhang Y, et al. 2018. The ultrasonic response of numerical simulation and analysis of scattering characteristics in digital core for shale reservoir. Chinese Journal of Geophysics (in Chinese), 61(3): 1069-1082. DOI:10.6038/cjg2018K0563 |
Wu C S. 2016. Economic Mathematics-Probability Theory and Mathematical Statistics (Second Edition) (in Chinese). Beijing: Higher Education Press.
|
Yang S H. 2014. Study on the adaptive and fast algorithm of gray scale image thresholding[Ph. D. thesis] (in Chinese). Chongqing: Chongqing University.
|
Zhang G Z, Chen J J, Chen H Z, et al. 2015. Prediction for in-situ formation stress of shale based on rock physics equivalent model. Chinese Journal of Geophysics (in Chinese), 58(6): 2112-2122. DOI:10.6038/cjg20150625 |
Zhang W H, Fu L Y, Zhang Y, et al. 2016. Computation of elastic properties of 3D digital cores from the Longmaxi shale. Applied Geophysics, 13(2): 364-374. DOI:10.1007/s11770-016-0542-4 |
Zhao L X, Qin X, Han D H, et al. 2016. Rock-physics modeling for the elastic properties of organic shale at different maturity stages. Geophysics, 81(5): D527-D541. DOI:10.1190/geo2015-0713.1 |
Zhou H, Fu L Y. 2018. Scattering and intrinsic components of attenuation through the spectral ratio method in ultrasonic laboratory experiment. Chinese Journal of Geophysics (in Chinese), 61(3): 1083-1094. DOI:10.6038/cjg2018L0359 |
Zimmerman R W, Somerton W H, King M S. 1986. Compressibility of porous rocks. Journal of Geophysical Research: Solid Earth, 91(B12): 12765-12777. DOI:10.1029/JB091iB12p12765 |
Zimmerman R W. 1991. Compressibility of sandstones. Amsterdam: Elsevier.
|
邓继新, 史謌, 刘瑞珣, 等. 2004. 泥岩、页岩声速各向异性及其影响因素分析. 地球物理学报, 47(5): 862-868. DOI:10.3321/j.issn:0001-5733.2004.05.018 |
王志伟, 符力耘, 张艳, 等. 2018. 龙马溪组页岩数字岩芯超声响应数值模拟及散射特征分析. 地球物理学报, 61(3): 1069-1082. DOI:10.6038/cjg2018K0563 |
吴传生. 2016. 经济数学——概率论与数理统计(第二版). 北京: 高等教育出版社.
|
阳树洪. 2014.灰度图像阈值分割的自适应和快速算法研究[博士论文].重庆: 重庆大学.
|
张广智, 陈娇娇, 陈怀震, 等. 2015. 基于页岩岩石物理等效模型的地应力预测方法研究. 地球物理学报, 58(6): 2112-2122. DOI:10.6038/cjg20150625 |
周浩, 符力耘. 2018. 超声实验中谱比法衰减的散射与本征吸收特性. 地球物理学报, 61(3): 1083-1094. DOI:10.6038/cjg2018L0359 |