2018, Vol. 61
2. 中国科学院地质与地球物理研究所, 北京 100029
2. Institute of Geology and Geophysics, Chinese Academy of Sciences, Beijing 100029, China
岩石物理研究是沟通石油地质参数和地震响应的关键手段.由反射地震资料获取地下介质的速度、密度、泊松比、各向异性参数，进而由这些参数获得地下介质的孔隙度、裂缝发育和走向、孔隙和裂缝含流体(或气体)情况，就有望实现由反射地震资料直接识别储层，实现岩性勘探的目标.传统的岩石物理研究主要是实验测量和等效介质理论这两个途径.实验方法通过测量特定温压条件下岩心样品的岩石物理参数特征及其变化规律，应用统计方法建立岩石物理弹性参数与储层参数关系的经验模型，给出由储层参数估算弹性参数的近似算法和公式.等效介质理论基于较为严格的物理理论，获取不同矿物组分的介质，以及介质包含裂缝、孔隙和孔隙含流体(或气体)情况下的等效物理模量，以此获得裂缝、孔隙和孔隙流体对岩石物理弹性参数影响的规律.两类方法已取得了很好的结果，但在实际应用中也遇到了一些困难(Dvorkin et al., 2008).实验方法遇到的主要困难是岩心样品和尺度的限制.实际采集到的岩心样品有时不能满足对特定地质特征，如孔隙度，在特定范围内变化的要求，这就给由测量结果回归经验模型带来数据不充分的问题.地震勘探中的地震波波长是百米量级，为准确获得地震勘探频率范围的岩石物理弹性参数，理想的情况是岩心样品要能包含这一尺度下的非均匀变化特征，而这无论从实验装置和岩心采样都是很难实现的；因此就地震勘探而言，强烈依赖于其波长尺度的岩石物理弹性参数是难以通过实验直接获得的，如地震波频率范围内的衰减特征.基于等效介质理论的理论方法的主要问题是理论模型的严格假设；现行的一些等效介质理论都对介质、孔隙和裂缝分布做了均匀或均匀分布(方向随机)的假设，因此将这些理论应用于具有明显非均匀特征的储层(例如页岩储层)时将产生偏差.
作为对实验和理论方法的补充，计算岩石物理方法近年来得到重视和发展.目前比较活跃的是岩心成像与数值模拟相结合的方法(Dvorkin et al., 2011；Diaz et al., 2013).这类方法首先利用束扫描电镜(FIB-SIM)和CT成像技术，获取岩心的小尺度(微米级)精细构造，然后利用数值模拟获得岩石物理参数(渗透率、电阻率、弹性模量).尽管相对实验测量这类方法有更多的灵活性，但由于电镜尺寸和计算能力的限制，这类方法所涉及的岩心尺度是毫米量级，较实验测量的厘米量级又低一个量级.由于尺度太小，这类方法在获取弹性参数时利用的是静力学模拟方法.理论上，针对更大尺度岩心应用数值模拟方法来获取岩石物理参数是可行的，但计算能力和获取大尺度精细地质模型的困难限制了这方面的应用.已有的工作(Deeks et al., 2012；Dvorkin et al., 2012)主要是针对完全的随机模型展开，由于计算机能力的限制，也还没有达到地震波波长的尺度.
页岩储层有别于其他种类储层的一个突出特点是其成层特征导致的较地层非均匀变化更小尺度的非均匀.由于较小尺寸的岩心不能完全包括这种非均匀变化特征，使得其等效岩石物理弹性参数不能由常规的岩心测量准确得到；而由这类非均匀导致的各向异性特征也难以利用远高于地震勘探频率的高频信号来测量得到.由于基于均匀分布假设的理论方法不能准确考虑页岩储层的这类小尺度非均匀，已有的研究工作(董宁等，2014；Qin et al., 2014)只能通过假设背景介质是有规律的层状介质(即假设背景介质是VTI各向异性)，来考虑储层的这类小尺度非均匀.
为获取页岩储层的岩石物理参数，特别是各向异性参数随TOC(Total Organic Carbon)含量等地质特征的变化规律，本文发展了应用数值计算方法获取页岩储层的速度、各向异性参数的计算岩石物理系列方法.针对获得地震波波长尺度地质模型的困难，本文利用地层沉积的统计学模型(Schwarzacher，1969；Elfeki and Dekking, 2001)，发展了基于Markov链的大尺度(地震波波长尺度)页岩储层精细(毫米级)数值建模的随机建模方法；针对地震波波长尺度计算模型涉及的计算能力限制问题，本文发展了计算网格尺度的地球物理参数等效建模方法，以此大幅减少了数值模拟的计算需求.以这两项技术为基础，文中进一步给出了应用弹性波数值模拟方法求取实际地震波频带对应的岩石物理弹性参数的方法.
本文提出的系列方法较好地解决了页岩储层岩石物理参数求取面临的不同地质特征的岩心获取和较地层非均匀变化更小尺度非均匀导致的尺度效应问题，可得到随孔隙、TOC含量等特征参数变化的岩石物理弹性参数的变化规律.本文工作也丰富了计算岩石物理方法，为计算岩石物理面临的地震波波长尺度地质建模和计算能力限制问题提供了有效的解决方案.
1 大尺度页岩储层精细数值建模储层的精细数值建模是计算岩石物理方法应用的基础.为准确获得地震勘探频率范围的岩石物理弹性参数，特别是考虑储层非均匀导致的尺度效应，对百米量级的岩心样品进行计算岩石物理的数值模拟是更合理的选择；而如何获得这样较大尺度的地质模型是这一思路能否实现的关键.
页岩储层有别于其他种类储层的一个突出特点是其成层特征导致的小尺度非均匀性，这一非均匀性导致较大尺度的模型不能通过数值岩心(毫米尺度)的简单放大“复制”来得到.为此，本文发展随机过程数值建模方法，利用地层沉积的统计学模型，来获取地震波波长尺度的页岩储层精细(毫米级)数值模型；而岩心样品，则是用来获得这一随机过程的关键控制因素的依据.
地质模型的非均匀是典型的分形(Fractal)结构，可以在不同尺度下无限细分.针对页岩储层的地质特征(刘惠民等，2012；王敏等，2013)，我们在建模中将垂向最小尺度定义在与页岩储层的纹层结构同量级(毫米)，而横向非均匀的尺度定义在“10 m”量级.
本文研究基于二维模型展开.将采用二维马尔科夫链模型生成地震波波长尺度的二维精细地质模型.为此，将以一定长度的“层”和“纹层”为单元，而每个“层”和“纹层”又可以由不同的岩相构成，不同类型岩相的“层”或“纹层”单元构成了组成地质模型的基本单元.由纵向和横向转移概率矩阵联合确定不同新事件(即基本单元)的发生概率，即可由这些基本单元随机构成二维地质模型.具体算法如下：
定义二维空间位置(i，j)，其中j代表纵向变化(即沉积历史)，j增加表明深度变浅，令Zi, j表示空间位置(i，j)处的基本单元类型；Sk(k=1, n)代表定义的基本单元，pmk是纵向转移概率矩阵中(m，k)处的元素，rlk是横向转移概率矩阵中(l，k)处的元素，则Zi, j=Sk的概率为(Elfeki and Dekking, 2001)：
|
(1) |
式中N是横向的单元总数，代表模型的横向大小；l源于Zi-1, j=Sl，m源于Zi, j-1=Sm，q源于ZN, j=Sq.
在没有井约束情况，为完成二维建模，需先应用一维马尔科夫链模型生成两端，即i=1, N处的单元类型Z1, j和ZN, j.为在此后的弹性波模拟计算中在横向应用循环边界条件，我们定义Z1, j=ZN, j.Z1, j=Sk的概率为
|
(2) |
式中m源于Z1, j-1=Sm.
在应用式(1)构建二维模型时，需对底层初始化，即确定最底层的单元类型.这可再次应用一维马尔科夫链模型, Zi, 1=Sk的概率为
|
(3) |
式中q源于ZN, 1=Sq，l源于Zi-1, 1=Sl.
由式(1)—(3)计算空间位置(i，j)处Sk(k=1, n)的发生概率，将这些概率顺序排列，组成一个长度为1的数轴；利用计算机产生随机数，随机数落在数轴的那个区间中，该处就是生成该区间对应的那类基本单元.
实际工作中，对多数不能明确给出横向转移概率矩阵的情况，可用将纵向转移概率矩阵中的对角元素放大一定的倍数来近似横向转移概率矩阵，这个倍数大小恰好反映了横向均匀的程度.
上述算法将产生的是一个由不同类型岩相的“层”或“纹层”单元构成的二维模型；而对每类基本单元，我们将根据其岩相的矿物组分，由随机过程生成分辨率为0.001 m×0.001 m的非均匀构造，每个0.001 m×0.001 m的小尺度介质对应一类具体的矿物，如方解石、黏土、石英或TOC.这样就可以得到一个0.001 m×0.001 m分辨率的、由基本矿物构成的大尺度地质模型.
生成基本单元的随机过程是由其岩相的矿物组分控制，根据矿物组分的百分比组成一个长度为1的数轴；对基本单元中的随机一个离散点(按0.001 m×0.001 m离散)，利用计算机产生随机数，随机数落在数轴的那个区间中，该点就对应那个矿物；基本单元中某一点确定一类矿物后，使用的数轴也要相应调整(该类矿物的百分比减少).TOC的分布还要考虑黏土的分布和“层”或“纹层”构造的情况.
生成大尺度地质模型后，还需进一步在模型上随机添加裂缝和孔隙.孔隙的添加将依据岩心样品统计得到的不同层位、不同岩相的孔隙率，通过指定模型中各离散点对应的0.001 m×0.001 m尺寸的矿物中的孔隙度来添加孔隙；对TOC的离散点，其孔隙对应有机孔.裂缝总量的多少由统计得到的裂缝密度决定.垂直缝的方向将由主应力方向控制.
完成上述过程，就可得到包含孔隙、裂缝、有机孔的分辨率为0.001 m×0.001 m的地震波长尺度的地质模型，该模型既考虑了结构特点(层与纹层结构)，也考虑了矿物组分比例(黏土、方解石、石英和TOC)，还包括了裂缝和孔隙.
下面将根据胜利罗家的页岩储层岩心分析结果(刘惠民等，2012；王敏等，2013)，给出一个页岩储层数值建模例子.根据岩心分析，确定3种不同岩相的“层”基本单元：含泥质灰岩层、泥质灰岩层、灰质泥岩层，其尺度定义为0.1 m×10 m；确定2种不同岩相的“纹层”单元：灰质纹层、泥质纹层，其尺度定义为0.001 m×10 m；共5个基本单元.储层的岩相构成、矿物占比见表 1，“层”与“纹层”占比为4:1；裂缝主要是水平裂缝，裂缝密度是指裂缝体积之和在总体积中的占比，同时考虑了裂缝数量和大小.
|
表 1 典型页岩储层的岩相构成、矿物占比及孔隙度和裂缝密度等参数 Table 1 Rock facies composition, mineral proportion, porosity and fracture density of typical shale reservoirs |
图 1给出了生成的数值模型的纹层状部分的局部.观察图 1可知，生成的数值模拟与实际观察的页岩岩心的纹层部分很相似.若继续使用0.001 m×0.001 m分辨率表示层状部分，由于可显示的尺度太小，不能显示多个层的变化规律，因此图 2给出了用岩相表达的数值模型层状部分的局部.
|
图 1 针对表 1页岩储层生成的数值模型的纹层状部分(局部) 图中黄色样点指示方解石，深蓝色样点指示黏土，天蓝色样点指示石英，红色样点指示TOC. Fig. 1 The laminar part (local) of the numerical model generated by shale reservoir in Table 1 The yellow point indicates calcite, dark blue point indicates clay, sky blue point indicates quartz, and the red sampling point indicates TOC, respectively. |
|
图 2 针对表 1页岩储层生成的数值模型的层状部分(局部) 图中黄色部分指示含泥质灰岩，深蓝色部分指示灰质泥岩，天蓝色部分指示泥质灰岩. Fig. 2 The layered part (local) of the numerical model generated by the shale reservoir in Table 1 The yellow part indicates muddy limestone, dark blue part indicates gray mudstone, and the sky blue part indicates argillaceous limestone, respectively. |
岩石物理的等效介质理论研究表明，等效弹性参数不仅决定于介质的矿物成分、孔隙和裂缝，还受矿物组分自身的几何形状影响(Kuster and Toksoz, 1974).第1节给出了用矿物成分表达的地震波长尺度的数值模型，但没有给出各个矿物组分的几何形状和孔隙的形状；实际上，即使给出矿物组分的几何形状，若在弹性波模拟的数值离散中考虑这些形状，将需要远小于0.001 m×0.001 m尺度的计算网格.就针对地震波长尺度的模型进行数值模拟而言，这样的网格尺度需要太多的计算资源和计算时间.
为避免在生成数值模型时考虑各矿物组分的几何形状的困难，也为大幅减少数值模拟的计算规模，本文利用等效介质理论，发展了计算网格尺度的地球物理参数等效建模方法.这一方法的核心是：根据已有的计算资源，确定尽可能小的计算网格尺寸(远小于频散要求)；利用测井速度估计矿物的几何形状系数；针对每个计算网格包括的矿物成分、孔隙、裂缝、有机孔、孔隙流体(或气体)和得到的几何形状系数，生成该计算网格内的等效拉梅系数和密度.
矿物的几何形状系数估计的主要步骤是：针对数值模型的井旁部分，在测井尺度上应用等效介质方法，得到测井尺度的P波和S波速度，根据等效速度与实测速度的差异，调整几何形状系数，两者接近时的几何参数就可认为是实际的几何形状系数.
岩石物理研究已发展了多种等效介质理论方法(Mavko et al., 2009)，国内外学者已将这些等效介质理论应用于页岩储层的等效弹性参数计算(董宁等，2014；Qin et al., 2014).但本文的等效化处理仅限于各个计算网格内，因此可保留页岩储层的层状非均匀特征；而常规的基于整个地质模型的等效方法均假定各类矿物成分、孔隙、裂缝是均匀分布的，不能有效地考虑储层的非均匀性.
本文生成计算网格尺度的等效弹性参数的步骤如下：
(1) 根据TOC中有机孔的孔隙率和含流体或气体情况，利用式(4)的自洽近似(SCA)理论，求得TOC的等效弹性参数；
(2) 根据计算网格内各类矿物成分，如TOC、方解石、石英、黏土的弹性参数，将孔隙作为真空，再次利用式(4)的自洽近似(SCA)理论(Hill，1965)求得等效的干岩的压缩和剪切模量；
|
(4) |
式中KSC*和μSC*分别是等效的干岩的压缩和剪切模量；N代表共有N种不同矿物(真空作为一种)，xi是第i种矿物体积比，Ki和μi分别是第i种矿物的压缩和剪切模量，P*i和Q*i是第i种矿物分别在模量为KSC*和μSC*的背景介质中的几何系数(Berryman，1980).式(4)需采用迭代方法求解；
(3) 若孔隙中充满流体，可根据Brown-Korringa广义Gassmann理论(Brown and Korringa, 1975)实现孔隙的流体替换，得到网格内所含矿物的流体饱和的等效压缩模量和剪切模量；
(4) 根据得到的等效体积和剪切模量，以及各个计算网格中包含裂缝的数量，由Hudson理论(Hudson，1981)得到含有裂缝的等效各向异性弹性参数.
针对每个计算网格完成上述步骤，就可完成计算网格尺度的地球物理参数等效建模，这一地球物理建模过程考虑了不同矿物组分、有机孔、孔隙和裂缝等介质特征.上述等效计算采用自洽近似(SCA)理论，是因为该理论对各类矿物组分的分布没有限制.
3 弹性波数值模拟提取岩石物理弹性参数本文通过模拟平面P波和S波在页岩储层中的传播，来提取储层的P波、S波速度和各向异性参数(Zhang and Gao, 2009).等效弹性参数是频率相关的，为获得服务于地震勘探的岩石物理参数，本文的弹性波数值模拟将采用与实际地震勘探相同主频的震源.
首先生成一个100 m厚、1000 m宽的地质模型(设计较宽的模型是为了方便求取各向异性参数)，通过对各计算网格进行等效计算，得到一个100 m厚、1000 m宽的非均匀弹性介质模型；较大的模型可减少拾取地震波波峰的误差，因此通过重叠该模型，可得到一个更大的，如400 m×1000 m的模型.本文通过拾取两个时刻的波场快照上的波前，来计算储层的物理参数，因此将储层模型放置到一个均匀背景中(见图 3).在均匀介质中，既容易激发平面波，也容易准确拾取平面波的波峰.图 3中的储层模型，就是由图 1和图 2地质模型通过等效算法得到的，其计算网格的尺寸是0.02 m×0.02 m.
|
图 3 提取弹性参数的弹性介质P波速度模型 图中色标对应速度值. Fig. 3 P-wave velocity model of an elastic medium for elastic parameters extraction Color scales corresponds to velocity values. |
在模型的均匀介质部位设置平面波震源，应用弹性各向异性介质中弹性地震波模拟方法(Zhang and Verschuur, 2002)，模拟平面波传播；记录两个时刻的波场快照，分别拾取两个波场快照中平面波波峰的位置，即可求得储层的弹性参数.若设置水平的平面P波震源，在t1和t2两个时刻记录波场快照(图 4)，在两个波场快照上拾取得到平面波峰值的坐标为x1和x2, 可得到储层的垂直方向P波速度为：
|
图 4 水平平面P波入射0.1 s (a)和0.25 s (b)时刻的波场快照 Fig. 4 Snapshots of horizontal plane P wave incidence at 0.1 s (a) and 0.25 s (b), respectively |
|
(5) |
式中L是重叠储层的厚度，vP0是均匀背景的P波速度；为减少在储层界面处的反射，可令1/vP0等于储层中各单元P波速度倒数的平均值.
若设置与水平方向夹角为θ的平面P波震源，可求得储层的沿角度θ方向的P波速度：
|
(6) |
式中h2-h1是两个波前面沿垂直波前面方向的距离，如图 5所示.图 5表明，为得到大角度方向的波速，必须生成较宽的地质模型.
|
图 5 斜入射平面P波0.1 s (a)和0.29 s (b)时刻的波场快照 Fig. 5 Snapshots of the oblique plane P wave incidence at 0.1 s (a) and 0.29 s (b), respectively |
取一段模型，将该储层旋转90°放置到背景介质中，就可采用与式(5)相同的方法求得水平方向的P波速度(vP)h.
若设置水平的平面S波震源，在t1和t2两个时刻记录波场快照(图 6)，在两个波场快照上拾取得到平面波峰值的坐标为x1和x2, 可得到储层的垂直方向S波速度(vS)v.对比图 4和图 6，可观察出明显的P、S波波速差别；平面P波的波长更长.
|
图 6 水平平面S波入射0.2 s (a)和0.38 s (b)时刻的波场快照 Fig. 6 Snapshots of the horizontal plane S wave incidence at 0.2 s (a) and 0.38 s (b), respectively |
利用已求得的垂直、平行、斜入射的P波速度，再利用式(7)的VTI介质相速度公式(Tsvankin，1996)，由最小二乘方法求得储层的各向异性参数:
|
(7) |
利用本文发展的计算岩石物理系列方法，可有效解决现行岩石物理方法面临的岩心样品和岩心尺度限制的问题.通过在储层精细数值建模中人为改变某些关键参数，就可以得到特定地质特征，如TOC含量，在特定范围内变化的系列模型；而对该系列模型应用数值模拟方法提取速度和各向异性参数，就可以得到各向异性参数随TOC改变的变化规律.
以表 1所示的胜利罗家的页岩储层为模板，生成了一组TOC含量从3%变化到21%的系列储层模型，系列储层模型的矿物占比、孔隙度、裂缝密度、纹层和层的比例见表 2；地球物理参数等效建模所使用的各类矿物的岩石物理参数如表 3所示.表 2中TOC含量等矿物占比的变化，是通过改变表 1中各岩相中矿物占比实现的；如对TOC含量为15%的储层模型而言，其泥质灰岩层的矿物占比改变为方解石47%、石英14.5%、黏土19%、TOC17%.对表 2的系列储层模型而言，其用岩相表达的数值模型与图 2是一致的(但各岩相包含的矿物含量已发生变化)，而用矿物成分表达的数值模型已与图 1有明显变化，代表TOC的红色样点明显增多.
|
|
表 2 TOC含量变化的储层矿物占比等关键参数 Table 2 Proportions of the reservoir minerals with changing TOC content |
|
|
表 3 各类矿物的岩石物理参数 Table 3 Rock physical parameters of different minerals |
图 7给出了储层垂直方向P波和S波波速随TOC含量变化的曲线，从图中可看出，随TOC含量的增加，P波速度迅速降低，S波速度也降低，但降低的幅度小于P波速度；图 8给出了储层各向异性参数ε和δ随TOC含量变化的曲线，可见随TOC含量的增加，ε和δ参数也快速增加，这表明TOC含量增加导致储层的各向异性迅速增加，因此各向异性的大小是反映储层TOC含量的重要指标.需指出的是，TOC含量增加导致储层各向异性增强是由黏土(TOC赋存其中)的成层性引起的.实际上，除TOC含量外，其他因素变化(如水平裂缝增加、纹层占比增加等)也会导致储层各向异性强度变化.
|
图 7 储层垂直方向P波(实线)和S波(虚线)速度随TOC含量变化曲线 VTOC是用小数表示的TOC的含量. Fig. 7 Perpendicular P wave (solid line) and S wave (dashed line) velocity versus TOC content in reservoirs VTOC is the content of the decimal representation of TOC. |
|
图 8 储层各向异性ε(实线)和δ(虚线)参数随TOC含量变化曲线 Fig. 8 Anisotropy parameters ε (solid line) and δ (dashed line) versus TOC content in reservoirs |
本文通过发展大尺度精细地质模型数值建模、计算网格尺度的地球物理建模和地震波数值模拟提取岩石物理弹性参数等系列方法，完整地实现了由数值计算方法获得复杂储层的岩石物理弹性参数，并将这一计算岩石物理方法应用于页岩储层地球物理建模，获得了储层速度和各向异性参数随TOC含量变化的量版.相比岩心测试方法，本文方法可得到岩心测试难以求取的与尺寸效应高度相关的岩石物理参数，避免了获取不同地质特征岩心的困难.不同于基于等效介质的理论方法，本文方法可精细考虑实际储层的较地层非均匀变化更小尺度的非均匀特征，因而可得到这类小尺度非均匀导致的宏观物理参数变化.本文研究为解决计算岩石物理方法面临的地震波波长尺度地质建模和计算能力限制问题提供了有效的方案.
Berryman J G. 1980. Long-wavelength propagation in composite elastic media I. Spherical inclusions. Journal of the Acoustical Society of America, 68(6): 1809-1819. DOI:10.1121/1.385171 |
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. DOI:10.1190/1.1440551 |
Deeks J, Lumley D, Shragge J. 2012. Computational rock physics: wave propagation velocities in partially saturated rocks. //SEG Technical Program Expanded Abstracts 2012. SEG, 1-5. http://www.researchgate.net/publication/268455135_Computational_rock_physics_wave_propagation_velocities_in_partially_saturated_rocks
|
Diaz E, Walls J, Marfisi N, et al. 2013. Effective multi-scale rock characterization using digital rock physics applied to La Luna Formation, Middle Magdalena Valley Basin, Colombia. //SEG Technical Program Expanded Abstracts 2013. SEG, 2382-2386. https://www.researchgate.net/publication/269127273_Effective_multi-scale_rock_characterization_using_digital_rock_physics_applied_to_La_Luna_Formation_Middle_Magdalena_Valley_Basin_Colombia
|
Dong N, Huo Z Z, Sun Z D, et al. 2014. An investigation of a new rock physics model for shale. Chinese Journal of Geophysics (in Chinese), 57(6): 1990-1998. DOI:10.6038/cjg20140629 |
Dvorkin J, Armbruster M, Baldwin C, et al. 2008. The future of rock physics:computational methods vs. lab testing. First Break, 26(9): 63-68. |
Dvorkin J, Derzhi N, Diaz E, et al. 2011. Relevance of computational rock physics. Geophysics, 76(5): E141-E153. DOI:10.1190/geo2010-0352.1 |
Dvorkin J, Fang Q, Derzhi N. 2012. Etudes in computational rock physics:Alterations and benchmarking. Geophysics, 77(3): D45-D52. DOI:10.1190/geo2011-0236.1 |
Elfeki A, Dekking M. 2001. A Markov chain model for subsurface characterization:Theory and application. Mathematical Geology, 33(5): 569-589. DOI:10.1023/A:1011044812133 |
Hill R. 1965. A self-consistent mechanics of composite materials. Journal of the Mechanics and Physics of Solids, 13(4): 213-222. DOI:10.1016/0022-5096(65)90010-4 |
Hudson J A. 1981. Wave speeds and attenuation of elastic waves in material containing cracks. Geophysical Journal International, 64(1): 133-150. DOI:10.1111/j.1365-246X.1981.tb02662.x |
Kuster G T, Toksoz M N. 1974. Velocity and attenuation of seismic waves in two-phase media:Part Ⅰ. Theoretical formulations. Geophysics, 39(5): 587-618. DOI:10.1190/1.1440450 |
Liu H M, Zhang S P, Wang P, et al. 2012. Lithologic characteristics of Lower Es3 shale in Luojia area, Zhanghua Sag. Petroleum Geology and Recovery Efficiency (in Chinese), 19(6): 11-15. |
Mavko G, Mukerji T, Dvorkin J. 2009. The Rock Physics Handbook. 2nd ed. Cambridge: Cambridge University Press.
|
Qin X, Han D H, Zhao L X. 2014. Rock physics modeling of organic-rich shales with different maturity levels. //2014 SEG Annual Meeting. Denver, Colorado, USA: SEG, 2952-2957. https://www.researchgate.net/publication/266316943_Rock_physics_modeling_of_organic-rich_shales_with_different_maturity_levels
|
Schwarzacher W. 1969. The use of Markov chains in the study of sedimentary cycles. Journal of the International Association for Mathematical Geology, 1(1): 17-39. DOI:10.1007/BF02047069 |
Tsvankin I. 1996. P-wave signature and notation for transversely isotropic media:an overview. Geophysics, 61(2): 467-483. DOI:10.1190/1.1443974 |
Wang M, Zhu J J, Yu G H, et al. 2013. The shale lithofacies characteristics and logging analysis techniques in Luojia area. Well Logging Technology (in Chinese), 37(4): 426-431. |
Zhang J F, Verschuur D J. 2002. Elastic wave propagation in heterogeneous anisotropic media using the lumped finite-element method. Geophysics, 67(2): 625-638. DOI:10.1190/1.1468624 |
Zhang J F, Gao H W. 2009. Elastic wave modelling in 3-D fractured media:An explicit approach. Geophysical Journal International, 177(3): 1233-1241. DOI:10.1111/gji.2009.177.issue-3 |
董宁, 霍志周, 孙赞东, 等. 2014. 泥页岩岩石物理建模研究. 地球物理学报, 57(6): 1990-1998. DOI:10.6038/cjg20140629 |
刘惠民, 张守鹏, 王朴, 等. 2012. 沾化凹陷罗家地区沙三段下亚段页岩岩石学特征. 油气地质与采收率, 19(6): 11-15. |
王敏, 朱家俊, 余光华, 等. 2013. 罗家地区泥页岩岩相特征及测井分析技术. 测井技术, 37(4): 426-431. |