白云岩储层是油气勘探中重要的储层类型之一, 该类储层孔隙类型复杂多变、矿物成分多样, 其岩石物理特征受到众多因素的影响^[1]。四川盆地中部高石梯地区(川中高石梯地区)台缘带的勘探实践表明, 灯四段储层的发育与储集性能都受到白云岩中硅质矿物的影响, 优质储层仅发育在硅质含量较低的白云岩地层中^[2], 研究区硅质矿物的含量从0变化至50%以上^[3], 其成因复杂, 矿物结构特征也各不相同, 硅质矿物及其结构特征对地层的弹性性质造成影响^[4], 并且含硅质白云岩与储层均具有低阻抗的特征, 如何区分储层与硅质层是目前储层地震预测的重点和难点问题。岩石物理模型作为联系储层特征与地震响应的桥梁, 可以分析与地震特征有关的岩石物理性质及其与地震响应之间的关系, 是解决上述问题的关键手段之一^[5]。因此, 为了明确储层参数与弹性参数之间的关系, 分析硅质矿物对地震响应特征的影响, 也为模型正演及储层预测提供依据, 需要建立能够描述硅质矿物及其结构特征的岩石物理模型^[6]。

目前最常用的碳酸盐岩岩石物理模型是XU等^[7]对Xu-White模型的推广^[8], 张秉铭等^[9]也对Xu-Payne模型进行了改进, 但该类方法的重点在于孔隙结构的描述——使用平均模型计算岩石矿物的弹性模量。LI等^[10]根据孔隙特征将岩石划分成小网格来计算弹性模量, 使用了平均模型来描述岩石矿物; 曹晓初等^[11]综合使用了自相容近似(SCA)模型和微分等效介质(DEM)模型, 但未说明构建岩石基质时所使用的几何因子含义。而当硅质矿物及其结构特征发生变化时, 上述模型针对基质矿物所使用的平均模型无法对这种特征进行有效揭示, 计算结果会出现一定偏差。PRIDE等^[12]引入固结参数建立了干岩模量与基质模量之间的关系, 但仅考虑了成岩作用的因素而忽视了孔隙结构对纵、横波速度的影响。CIZ等^[13]和KUMAR等^[14]均提出了固体替换方程以描述岩石中固体包含物, 该方程目前多用于含稠油储层以替代Gassmann方程^[15]。总的来说, 现有的碳酸盐岩岩石物理模型很难描述硅质矿物及其结构特征的变化, 无法准确反映研究区储层特征。

本文在描述孔隙结构特征的岩石物理模型基础上, 利用固体替换等对岩石基质模量求取的方法进行改进, 提出了一种新的白云岩岩石物理模型, 并且将其应用到川中高石梯地区灯四段台缘带的白云岩储层, 通过横波速度预测验证模型的有效性, 以期为后续储层预测工作奠定基础。

1 灯四段白云岩特征

川中高石梯地区在灯影组沉积时期西侧为德阳-安岳裂陷, 东侧是台缘-台内的碳酸盐岩台地沉积体系, 其中发育的丘滩亚相是优质储集相带^[16]。在灯四段沉积后, 桐湾运动的抬升作用带来了长期的风化剥蚀及淋滤, 在灯影组顶部形成了一套广泛发育的岩溶风化壳, 再加上长时间的埋藏, 使得灯四段岩性及储集条件相对复杂, 因而影响储层岩石物理特征的因素众多。

灯四段台缘带的主要储集岩类型为藻白云岩, 包含的主要矿物类型有白云石、硅质矿物(主要为石英)等, 黏土平均含量小于1%。该地区灯四段具有普遍含气的特征, 其富集的优劣与储集条件密切相关, 而储集条件则与其中分布的硅质矿物呈现此消彼长的关系。研究区硅质含量与孔隙度交会(图 1a)显示硅质含量与孔隙度呈现负相关关系, 硅质含量较高的样点孔隙度相对较低、储层也不发育; 纵波速度与孔隙度交会(图 1b)显示孔隙度较高的储层具有低速特征, 但硅质含量较高且孔隙度较低的硅质层中(图中虚线框区域)也具有低速特征, 即无法在低速区分开储层与硅质层, 这是该地区目前地球物理储层预测存在多解性的主要原因。

另一方面, 实际测井数据表明硅质矿物与储层的岩石物理特征具有明显相关性。图 2是台缘带典型井(GS1井)测井数据中纵、横波速度与孔隙度、硅质含量的交会图。由图 2可见: 随着孔隙度的增加, 纵、横波速度整体呈减小趋势; 速度的分散情况随着硅质矿物含量变化呈现一定的规律性(图 2a和图 2b), 在孔隙度相近时, 岩石的纵波速度随着硅质含量的增多而略微减少, 而横波速度则随着硅质含量的增多而略微增加。孔隙度相近且硅质含量相等时(图 2c和图 2d), 纵、横波速度数据点较分散, 其中, 纵波速度相差可达0.5 km/s, 横波速度相差可达0.2 km/s; 在孔隙度较低时(小于0.5%), 相似硅质含量样点的纵波速度仍然有超出0.5 km/s的变化范围, 横波速度也有一定的分散(图 2c和图 2d中红色虚线框)。因此, 确定硅质矿物含量及分布规律对评价灯影组储集条件、研究储层地球物理响应特征具有极其重要的意义。

研究表明, 硅质矿物广泛分布于研究区的地层中, 局部地区可见硅质条带^[17]。研究区硅质矿物主要有两种成因: ①次生充填或交代的硅质矿物; ②沉积形成的硅质岩^[18]。两种不同成因的硅质矿物具有不同的结构特征: 次生的硅质矿物一般晶粒粗大, 主要是充填或者交代溶洞或裂缝(图 3a); 沉积形成的硅质选择性交代菌藻类及其它碳酸盐岩颗粒, 主要表现为细晶或微晶、放射状玉髓以及隐晶质硅质^[19](图 3b和图 3c), 宏观结构上单层厚度较薄, 一般为条纹、条带结构(图 3d)。不同成因的硅质矿物其结构特征不同, 岩石物理特征也存在差异, 在岩石物理建模时需要对这些矿物的结构特征加以区分。

2 针对硅质矿物结构特征的岩石物理模型 2.1 岩石基质的建模方法

当地震波波长远大于岩石的非均匀尺度(颗粒尺度)时, 岩石就可以等效为统计意义上的均匀介质^[2], 在地震岩石物理建模时通常使用等效介质模型的概念来描述和表征岩石特性^[20]。对于岩石这种混合物的等效模量进行估算时一般依赖于3个条件^[21]: ①各种矿物的弹性模量; ②各种矿物的体积分数; ③各个成分如何相互组合在一起的几何细节(结构特征)。若只考虑条件①和②, 则可以通过界限平均模型(Voigt-Reuss-Hill(VRH)模型、Hashin-Shtrikman界限平均)来求解; 若考虑组合细节, 则会假设一种具有特殊形状的包含物, 这类包含物模型主要有KT模型、DEM模型、SCA模型等^[22]; 在目前常用的主要用于岩石基质的3种包含物模型中加入空孔隙^[5], 硅质矿物被视为岩石基质(各种矿物的混合物), 通常不会考虑其矿物结构特征。要对混合矿物组合在一起的细节(矿物结构)进行描述, 这3种模型理论上也适用, 即将硅质视为一种包含物, 利用等效纵横比参数描述硅质矿物, 从而求取能够描述硅质矿物结构特征的岩石基质弹性模量。

图 4展示了利用这些模型向白云石中直接加入石英时的理论计算结果。其中, 在使用任意高宽比的椭球状包含物时(公式(1)), 发现随着孔隙纵横比参数的变化(0.01~0.99), 利用3种模型计算的结果变化非常小, 并且均位于VRH平均值附近。

$ \left\{\begin{array}{l} P=\frac{1}{3} T_{i i j j} \\ Q=\frac{1}{5}\left(T_{i j i j}-\frac{1}{3} T_{i i j j}\right) \end{array}\right. $

(1)

式中: P和Q为待求解的椭球状形状因子; T_iijj和T_ijij分别为与孔隙纵横比、基质和包含物模量有关的几何因子以及与孔隙纵横比有关的几何因子, 反映第i种孔隙对岩石基质的影响, 具体计算方式可参考文献[21]。

与此同时, 对于研究区的硅质矿物而言, 其结构形态也更接近于硬币状缝隙的假设。所以, 在描述矿物结构特征时, 需要使用硬币状缝隙的包含物形状因子:

$ \left\{\begin{array}{l} P^{(\mathrm{m}, \mathrm{i})}=\frac{K_{\mathrm{m}}+\frac{4}{3} G_{\mathrm{i}}}{K_{\mathrm{i}}+\frac{4}{3} G_{\mathrm{i}}+{\rm{ \mathsf{ π} }} \alpha \beta_{\mathrm{m}}} \\ Q^{(\mathrm{m}, \mathrm{i})}=\frac{1}{5}\left[1+\frac{8 K_{\mathrm{m}}}{4 G_{\mathrm{i}}+{\rm{ \mathsf{ π} }} \alpha\left(G_{\mathrm{m}}+2 \beta_{\mathrm{m}}\right)}+\right. \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \left.2K \frac{K_{\mathrm{i}}+\frac{2}{3}\left(G_{\mathrm{i}}+G_{\mathrm{m}}\right)}{K_{\mathrm{i}}+\frac{4}{3} G_{\mathrm{i}}+{\rm{ \mathsf{ π} }} \alpha \beta_{\mathrm{m}}}\right] \end{array}\right. $

(2)

式中: α为硬币状包含物的等效纵横比, 取值区间为(0, 1);K_m和G_m为基质矿物的弹性模量(K, G分别代表体积模量和剪切模量, 下同); K_i和G_i为包含物的弹性模量; β=G(3K+G)/(3K+4G), 其下标代表对应的弹性模量, 例如β_m代表式中K_m和G_m。

由图 4可以发现, 随着石英含量的变化, 大部分的模型结果都位于Hill平均值之上; 对于体积模量K, 3种模型的结果随α的增大均会很快超出Voigt上限值, 并且KT模型和DEM模型的结果位于Hill平均值和Voigt上限之间的α范围更窄; 对于剪切模量G, 由于白云石和石英的剪切模量相近(本文使用岩石物理弹性参数见表 1), VRH模型中Voigt上限与Reuss下限范围很窄, 3种包含物模型结果也都随着α的增大很快超出Voigt上限值。

VRH模型中Voigt上限与Reuss下限描述了矿物组合时的2种极限状态, 多种矿物混合时的模量值一般应在该界限值范围之内。因此, 直接利用这3种常用包含物模型向白云石中加入硅质矿物时限制较多, 硅质矿物的等效纵横比只能在较小的范围内变化, 存在适用性问题。所以要描述硅质矿物的结构特征, 可以参照处理孔隙结构的方法: 将硅质矿物视为一种包含物, 先使用一种模型向其中加入“空孔隙”(硅质矿物所占体积百分数), 然后利用固体替换方程将硅质矿物加入其中。

2.2 利用固体替换求取岩石基质弹性模量

含硅质白云岩主要组成矿物有白云石、硅质矿物(主要为石英)以及极少的黏土等。因此在使用固体替换方程时, 岩石骨架矿物只有白云石一种, 硅质矿物(体积分数为V_si)和黏土(体积分数为V_sh)都作为需要考虑矿物结构特征的固体包含物加入其中。对于空固体孔隙的加入, 由于固体包含物的体积分数可从0.1%变化至30%以上, 而KT模型计算高孔隙时、SCA模型计算较小孔隙纵横比时会产生不可靠的极值^[5], 因此在固体替换时使用了适用性更广的经典DEM模型, 同时为了避免DEM模型计算中不同矿物加入顺序导致的误差, 可以将硅质矿物和含量极少的黏土视为具有相同矿物结构特征的混合物(等效纵横比参数一致)。对于含硅质白云岩, 将硅质矿物和黏土视为包含物, 进行固体替换的具体步骤如下。

1) 利用DEM模型(公式(3))向纯白云石加入硅质矿物和黏土两种矿物所占体积分数(V_si+V_sh)的空孔隙, 得到固体岩石骨架的有效弹性模量值。该步骤中使用的孔隙形状因子(P, Q)为硬币状裂隙(公式(2))。

$ \left\{\begin{array}{l} \left(1-\varphi_{\mathrm{m}}\right) \frac{\mathrm{d}\left(K_{*}\right)}{\mathrm{d}\left(\varphi_{\mathrm{m}}\right)}=\left(K_{\mathrm{i}}-K_{*}\right) P^{(*, \mathrm{i})} \\ \left(1-\varphi_{\mathrm{m}}\right) \frac{\mathrm{d}\left(G_{*}\right)}{\mathrm{d}\left(\varphi_{\mathrm{m}}\right)}=\left(G_{\mathrm{i}}-G_{*}\right) Q^{(*, \mathrm{i})} \end{array}\right. $

(3)

式中: K_*和G_*分别为待求解的固体岩石骨架的有效弹性模量值; K_i和G_i分别为包含物弹性模量, 此处为空孔隙, 该参数设为0;φ_m为固体包含物体积分数, 即V_si+V_sh, 由ECS元素俘获测井数据给出; P^{(*, i)}和Q^{(*, i)}分别为硬币状缝隙的包含物形状因子, 具体计算方式见公式(2), 角标*表示将公式(2)中K_m, G_m替换成K_*, G_*, 公式(2)中α即变为描述硅质矿物结构特征的等效纵横比参数α_si。由于P和Q系数中包含待求解参数, 该式一般是耦合的, 可使用四阶Runge-Kutta法来求取其近似解, 初始解为固体岩石骨架矿物的弹性模量值, 即白云岩体积模量值。

2) 利用VRH平均计算硅质矿物和黏土混合弹性模量值, 得到固体包含物的弹性模量值。

3) 利用固体替换方程计算出固体饱和岩石的弹性模量值, 该值就是前文所述的“新的描述硅质矿物特征的岩石基质弹性模量”。步骤中的固体替换方程为各向异性Gassmann方程推广到弹性固体填充孔隙空间条件下的结果^[13]。该方程包含有与孔隙中固体相关的新定义参数^[15], 实际应用中, 对方程进行了简化, 最终固体替换方程可以写成:

$ \left\{\begin{array}{l} K_{\mathrm{m}_{-} \mathrm{s}}^{-1}=K_{\mathrm{m}_{-} \mathrm{d}}^{-1}-\frac{K_{\mathrm{m}_{-} \mathrm{d}}^{-1}-K_{0_{-} \mathrm{m}}^{-1}}{\varphi_{\mathrm{m}}\left(K_{\mathrm{f}_{-} \mathrm{m}}^{-1}-K_{0_{-} \mathrm{m}}^{-1}\right)+\left(K_{\mathrm{m}_{-} \mathrm{d}}^{-1}-K_{0_{-} \mathrm{m}}^{-1}\right)} \\ G_{\mathrm{m}_{-} \mathrm{s}}^{-1}=G_{\mathrm{m}_{-} \mathrm{d}}^{-1}-\frac{G_{\mathrm{m}_{-} \mathrm{d}}^{-1}-G_{0_{-} \mathrm{m}}^{-1}}{\varphi_{\mathrm{m}}\left(G_{\mathrm{f}_{-} \mathrm{m}}^{-1}-G_{0_{-} \mathrm{m}}^{-1}\right)+\left(G_{\mathrm{m}\_{\mathrm{d}}}^{-1}-G_{0\_{\mathrm{m}}}^{-1}\right)} \end{array}\right. $

(4)

式中: K_{m_s}, G_{m_s}为待求解的固体饱和岩石的弹性模量值; K_{0_m}, G_{0_m}为固体岩石骨架矿物—白云岩的弹性模量值; K_{m_d}, G_{m_d}为固体岩石骨架的有效弹性模量值, 分别等于公式(3)的中K_*和G_*; K_{f_m}, G_{f_m}为固体包含物的弹性模量值。

图 5显示了固体替换的理论计算结果, 变化的参数为石英含量(V_si)及包含物等效纵横比参数。在不含黏土的(V_sh=0)模拟中(图 5a), K_{m_s}值不会小于Reuss下限, 并且仅在α值较大时超出Voigt上限值(随着V_si增大, α的上限值约从0.95逐渐减小为0.70, 当α≈-0.154 9V_si²+0.391 5V_si+0.271 7时, K_{m_s}=K_VRH平均); G_{m_s}与VHR平均计算结果基本相同(图 5b)。在含5%黏土(V_sh=5%)的情况下, K_{m_s}的变化特征与不含黏土时特征基本一致(图 5c), G_{m_s}仅在V_si < V_sh时会在等效纵横比值较大时略微高于Voigt上限值(图 5d); 同时, G_{m_s}的模拟结果出现先增大后又逐渐减小的规律(图 5d), 经实验分析发现, 该特征是由于泥质含量大于硅质含量而导致。泥质矿物具有更低的剪切模量值, 当其含量对包含物造成的影响较大时, 会对计算结果造成更大的影响。对于研究区而言, 硅质含量普遍远高于泥质含量, 且泥质含量大部分小于1%, 因此该方法比直接使用包含物模型计算得到的结果能够更好地描述硅质矿物的结构特征。

2.3 岩石物理建模流程

根据含硅质白云岩的特征, 研究中将岩石等效为岩石基质、岩石骨架、混合流体和饱和岩石4个部分, 与其它模型相比, 本研究主要针对岩石基质模量的求取进行了改进。新提出的针对硅质白云岩矿物结构特征的建模方法是利用多个岩石物理模型组合来构建的(图 6), 具体流程如下。

1) 求取岩石基质模量(K_{m_s}, G_{m_s})。使用固体替换方法向“纯白云石”中加入“硅质矿物和黏土”。该步骤中的“纯白云石”实际上指的是不具有结构特征的骨架矿物(例如白云石、石灰石), “硅质矿物和黏土”是指需要描述结构特征的部分矿物(也就是作为包含物的部分, 例如硅质、黏土等)。对于具有更多矿物的情况, 该步骤中应当依据矿物是否具有结构特征, 分别利用VRH平均计算骨架矿物以及包含物矿物的模量。

2) 计算干岩石骨架弹性模量。对于缝洞型储层, 一般需要分析多种孔隙类型对岩石物理特征的影响, 而DEM模型在加入多种孔隙时添加顺序会影响结果, SCA模型不适于硬币状裂隙的计算, 因此这一步中使用KT模型向岩石基质中加入空孔隙, 计算干岩石骨架的弹性模量:

$ \left\{\begin{array}{l} \left(K_{\mathrm{dry}}-K_{\mathrm{m}_{-} \mathrm{s}}\right) \frac{3 K_{\mathrm{m}_{\mathrm{s}}}+4 G_{\mathrm{m}_{-} \mathrm{s}}}{3 K_{\mathrm{dry}}+4 G_{\mathrm{m}_{-} \mathrm{s}}}=\sum\limits_{i=1}^{M} w_{i}\left(K_{i}-K_{\mathrm{m}_{-} \mathrm{s}}\right) P \\ \left(G_{\mathrm{dry}}-G_{\mathrm{m}\_\mathrm{s}}\right) \frac{G_{\mathrm{m}\_{\mathrm{s}}}+\zeta_{\mathrm{m}\_{\mathrm{s}}}}{G_{\mathrm{dry}}+\zeta_{\mathrm{m}_{-} \mathrm{s}}}=\sum\limits_{i=1}^{M} w_{i}\left(G_{i}-G_{\mathrm{m}_{-} \mathrm{s}}\right) Q \end{array}\right. $

(5)

式中: K_dry, G_dry分别为待求解的干岩石骨架的弹性模量; ζ_{m_s}=(G_{m_s}/6)×(9K_{m_s}+8G_{m_s})/(K_{m_s}+2G_{m_s}); M为包含物数量; w_i为每种孔隙的体积百分比, 需迭代计算。研究表明, 岩石物理模型中更关注的是孔隙的形状特征(即孔隙纵横比)而不是其成因分类^[5], 结合缝洞型碳酸盐岩储层的孔隙结构特征, 本文模型将孔隙分为硬孔隙(孔隙纵横比为0.8, 体积分数为w_c)和软孔隙(孔隙纵横比为0.02, 体积分数为w_s, 2种孔隙体积分数之和为孔隙度)。

3) 使用Wood方程计算混合流体的体积模量:

$ K_{\mathrm{f}}=\left(\frac{S_{\mathrm{w}}}{K_{\mathrm{w}}}+\frac{S_{\mathrm{g}}}{K_{\mathrm{g}}}\right)^{-1} $

(6)

式中: K_f为混合流体体积模量; S_w和S_g分别为水、气的饱和度, S_g=1-S_w; K_w和K_g分别为水、气的体积模量。

4) 使用Gassmann方程计算饱和岩石的体积模量:

$ \left\{\begin{array}{l} K_{\mathrm{sat}}=K_{\mathrm{dry}}+\frac{\left(K_{\mathrm{m}_{-} \mathrm{s}}-K_{\mathrm{dry}}\right)}{K_{\mathrm{m}_{-} \mathrm{s}}\left(1-\varphi-\frac{K_{\mathrm{dry}}}{K_{\mathrm{m}_{-} \mathrm{s}}}+\varphi \frac{K_{\mathrm{m}_{-} \mathrm{s}}}{K_{\mathrm{f}}}\right)} \\ G_{\mathrm{sat}}=G_{\mathrm{dry}} \end{array}\right. $

(7)

式中: K_sat和G_sat分别为流体饱和岩石的弹性模量; φ为岩石的孔隙度。

5) 利用弹性参数计算饱和流体岩石的纵波速度和横波速度:

$ v_{\mathrm{Pc}}=\sqrt{\frac{3 K_{\mathrm{sat}}+4 G_{\mathrm{sat}}}{\rho_{\mathrm{sat}}}} $

(8)

$ v_{\mathrm{Sc}}=\sqrt{\frac{G_{\mathrm{sat}}}{\rho_{\mathrm{sat}}}} $

(9)

式中: v_Pc和v_Sc分别指计算得到的纵、横波速度; ρ_sat为饱和流体岩石的密度, 可以通过加权平均计算或密度测井获取。

在本文方法中, 无法直接获取的参数共有两种, 分别是第1)步中描述矿物结构特征的等效纵横比α_si和第2)步中描述岩石孔隙结构特征的两种体积分数w_c和w_s。目前, 多数方法都会依据已有的纵波速度对这些参数进行反演计算, 然后将反演得到的参数再代入模型中即可计算得到横波速度。由于目前多数学者认为孔隙结构是影响碳酸盐岩弹性模量的一种重要因素^[23-24], 即以横波速度为目标时, 横波速度主要受剪切模量影响, 而石英的剪切模量对岩石的剪切模量影响较小, 因此孔隙对横波速度的影响更大。通过理论计算可以知道随着体积分数w_c或α_si的变化呈单调增大或减少, 同时迭代计算两个参数会变得更复杂; 图 2的线性拟合结果也显示孔隙度与纵、横波速度之间的相关性略好于硅质矿物, 因此为了减少计算量, 文中使用了相对简单的优先迭代孔隙体积分数的反演方法, 在孔隙度较低或者其它变化孔隙度体积分数仍然不能得到准确解的情况下, 再迭代硅质矿物的等效纵横比参数以获取更准确的解。

图 7为本文方法的结构参数迭代反演计算流程图, 图中目标函数f为:

$ f=\left|{v}_{\mathrm{Pc}}-v_{\mathrm{P}_{-} \text {measured }}\right| $

(10)

式中: v_{P_measured}为测井测量得到的纵波速度值。

本文采用模拟退火迭代反演方法求取各参数, 具体步骤如下。

1) 先依据实际资料情况设定一个合适的初始α_{si_0}(矿物迭代次数用下角标_t表示, 下角标_0代表初始值): 该值选取的方式可以依据储层实际情况估算, 也可以通过固定孔隙参数然后经快速检索找到一个误差最小的值作为初始值(快速检索时α_si可设为0.01, 0.05, 0.10, 0.20, 0.50, 0.75, 0.99, 然后代入上一个采样点处的孔隙百分比计算, 这种方式更能体现地层中矿物的特征)。初始w_{c_0}(由于w_s=1-w_c, 可视为同一个迭代参数, 孔隙迭代次数用下角标_n表示)可设置成50%或者其它经验值, 然后计算初始目标函数值f_{_0}, 令当前解f=f_{_0}转到步骤2)。

2) 扰动产生新解w_{c_n}(此时孔隙迭代次数n=n+1), 计算新目标函数值f_{_n+t}, 转到步骤3)。参数扰动在实际中建议进行优化设计, 例如当f_{_n+t}大于实数a时, 孔隙体积分数以0.5%进行变化, 否则以0.01%变化, a的取值应保证当变化0.5%的体积分数时, 纵波速度变化小于a值。

3) 判断f_{_n+t}是否小于当前解f, 是则转到步骤4), 否则转到步骤5)。

4) 接受w_{c_n}, α_{si_t}为w_c, α_si的解, 然后转到步骤6)。当α_si未被扰动时, 则解为α_{si_0}。

5) 按Metropolis准则接受新解, 即当rand < exp[－(f_{_n+t}－f)/t_max－t]时接受w_{c_n}, α_{si_t}分别为dw_c, α_si的解(rand为随机数), 否则维持原解。然后转到步骤6)。

6) 判断n是否小于n_max(孔隙度迭代次数上限)。是则转到步骤2), 否则转到步骤10)。

7) 扰动产生新解α_{si_t}(此时矿物迭代次数t=t+1), 计算目标函数f_{_n+t}, 然后转到步骤8)。

8) 判断t是否小于t_max(矿物迭代次数上限), 是则转到步骤3), 否则转到步骤9)。

9) 令α_{si_t}=α_si(即令当前被扰动的α_{si_t}还原为目前最优解), 令孔隙迭代次数n等于0, 即重置孔隙迭代次数, 然后转到步骤2)。

10) 判断在当前解情况下是否小于给定的极小值e(e依据精度要求设定, 建议设为10~50 m/s, 精度越高, w_c扰动间隔需要设置得更小), 是则转到步骤11), 否则转到步骤7)。为了避免始终无法找满足条件的解, 建议在此处额外设置一个循环计数器(若计数扰动w_c的次数, 上限应设为2n_max), 当达到计数上限仍然未找到满足输出解时则仍转到步骤11)。

11) 结束运算, 返回w_c, α_si的解, 计算v_Sc并输出。

3 应用效果分析

基于研究区的实际测井资料, 应用本文方法对川中高石梯地区的GS1井进行了横波速度预测。图 8给出了该井横波速度预测时所需的多矿物测井解释估算矿物含量、孔隙度和密度的测井曲线以及3种模型的横波速度预测结果(模型中使用物理参数见表 1), 3种模型分别为: ①Xu-Payne模型(孔隙类型分类方法基于KUMAR等^[14]), 该模型虽然将孔隙进行了分类, 但是分类方法依据的是纵波速度与孔隙度之间的关系模板, 对于缝洞型储层该模板适用性较差(该方法只存在孔隙纵横比0.02和0.1以及0.1和0.8两种组合); ②张秉铭^[9]等针对复杂缝洞型碳酸盐岩储层的改进模型, 该模型将孔隙简化为软、硬孔隙, 并且依据纵波速度直接迭代计算孔隙的百分比参数, 避免了图版划分的适用性问题, 但该方法不能分析硅质矿物对岩石物理特征的影响; ③本文的岩石物理模型, 本文模型在描述软、硬孔隙的基础上, 加入了一个描述矿物结构特征的硅质等效纵横比参数, 使模型的适用性和物理含义更广泛。

图 8中纵波计算结果显示Xu-Payne模型中硅质矿物含量高的部位作为约束的纵波速度误差较大, 这是由于KUMAR等^[14]在对孔隙类型分析时未考虑硬孔隙和软孔隙的组合, 即该模型假设无法描述岩溶缝洞型储层孔隙特征。其它2种方法, 作为迭代计算约束的纵波速度与实测速度基本一致, 仅在部分测井曲线异常值处(孔隙度小于0.1%或含水饱和度为100%时)存在一定误差, 符合迭代计算的要求。但即使同样是迭代目标, 可以发现在5 100 m深度处硅质含量较高的层位, 利用本文模型迭代出的纵波速度精度也高于张秉铭模型。在横波速度预测的结果中, Xu-Payne模型在目的层顶部和底部误差较小, 中间硅质矿物发育的层位误差较大; 张秉铭模型在硅质矿物较发育的层位误差相对较小, 而其它层位误差较大; 利用本文模型的预测结果偏差整体较小, 曲线形态明显优于前2种方法并更接近于实际横波速度。

就孔隙体积分数结果而言, 利用张秉铭模型与本文模型预测结果在趋势上存在一定相似性, 这种相似反映了孔隙体积分数的准确性: 即不同类型的孔隙体积分数在理论上不变, 因此不同模型对于孔隙体积分数的变化趋势一致, 与理论假设相符。

依据硅质矿物体积分数的结果可以看出, 本文方法显示硅质矿物在5 100 m处硅质含量较高的层段处取值均偏向于0, 而在其它层段偏向于1。偏向于1说明该处属于充填原有的溶孔、溶洞或原生孔隙, 与硅质后期充填成因对应; 偏向于0说明硅质矿物成片展布与原生沉积成因对应, 或充填原有的裂缝。测井曲线显示5 100 m处硅质含量较高且连续一定厚度, 此段的硅质矿物成因应为沉积成因。因此本文方法中硅质矿物结构参数对岩石物理特征的描述是有效的。

图 9统计了不同偏差区间内的落点个数, 总体上反映出预测数据偏离实际数据的区间范围, 结果显示新模型的偏差程度在3种模型中最小。为了定量对比3种速度预测方法的差别, 分别计算了每组数据的均值、方差和Pearson相关系数(表 2), 由表 2可以看出: 本文方法得到的横波速度均值与原始数据差异最小, 相关程度最高; 同时双尾显著性检验结果表明仅有本文方法在0.05的水平下相关性显著。

除横波速度预测外, 岩石物理模型另一个更重要的用途在于对储层地震响应特征的分析。图 10展示了通过岩石物理模型补全数据后储层分类弹性模量的分布, 结果显示硅质层(硅质含量>15%)的体积模量总体最低而非储层的总体最高, 气层的剪切模量总体最低而硅质层的最高。图 10说明硅质层作为研究区的特殊岩性, 从测井岩石物理的角度具有一定的区分性, 但是其特征会对储层的辨别造成一定的影响。由于研究区储层与硅质均具有低阻抗的特征, 因此在下一步的工作中, 需要利用岩石物理模型建立起储层与硅质层的识别模板, 例如流体替代、AVO模型分析等方法找出地震属性与硅质层之间的联系, 从而在低阻抗或低速区找出真正的储层分布范围。

4 结论及建议

针对川中高石梯地区白云岩储层的特征, 综合利用固体替换方程等方法构建了一个新的岩石物理模型。新模型可在纵波速度的约束下求取孔隙的体积分数和硅质矿物的等效纵横比参数, 相较于传统模型增加了对矿物结构特征的描述方法。川中高石梯台缘带GS1井的新模型横波速度预测结果与实测速度基本吻合, 各种统计结果也表明新模型的预测结果优于Xu-Payne模型和前人提出的复杂碳酸盐岩模型。

通过本文方法可以进一步使用流体替代、AVO模型分析等手段精确分析储层、硅质层与各种地震反演属性之间的关系, 从而从低阻抗区找出真正的储层分布区域。本文给出的固体替换方法相较于普通孔隙包含物模型能够更好地描述矿物结构特征, 为复杂矿物结构特征的岩石物理模型提供了参考。在矿物更复杂的情况下和其它岩石类型中, 则需要仔细分析矿物结构及孔隙结构与岩石物理参数之间的关系。