基于Wood阻抗法的悬浮模式天然气水合物饱和度预测方法*

2019-10-24鲍祥生

中国海上油气 2019年5期

鲍祥生侯斌

(1.广东石油化工学院广东茂名 525000；2.油气藏地质及开发工程国家重点实验室(成都理工大学)四川成都 610059；3.江苏油田江苏扬州 225009)

天然气水合物是由水和天然气在低温高压条件下形成的一种笼形化合物，是一种类冰状的物质[1],该物质没有固定的化学分子式，是一种非化学计量的混合物[2]，主要分布在陆上的永久冻土带和海域的海底沉积物中[3]。水合物燃烧的产物是二氧化碳和水，因此属于一种高效洁净的非常规能源，是未来化石能源的重要接替能源[4]。水合物大部分分布于海域，全世界范围内海域水合物分布占到了98%，研究表明水合物主要分布在受到温压条件控制的稳定带内[5-11]。但稳定带只是决定水合物存在的空间范围，并不对水合物富集程度起决定性影响。水合物在稳定带某个地区富集还受到气源条件、气体运移条件、储集条件等因素制约。

目前预测水合物富集有许多地震方法[1,5-6,8-10,12-13]，如BSR法、振幅空白带法、属性预测法、波阻抗预测法、水合物饱和度预测法等。依据量化特征可将这些方法分为3类：①定性预测方法，如BSR法、振幅空白带法，该类方法具有使用简单直观、易于理解等优点，但定量化程度不够，不利于判断不同位置富集程度差异；②半定量预测方法，如属性预测法，该类方法具有一定的量化特征，部分方法已能初步反映出富集程度差异性，但富集程度差异也是相对的；③定量预测方法，如波阻抗预测法、水合物饱和度预测法，该类方法能够绝对反映水合物的富集程度差异。在这3类方法中，定量预测方法是实际应用中预测水合物富集最受欢迎的一类方法，其中水合物饱和度预测法与水合物资源量计算最相关[14-15]，所以该方法对一个地区的水合物商业化进程具有重要意义，目前针对该方法的研究已成为水合物定量预测的热点课题。由文献[14-19]可知，Timur 时间平均方程法(简称Timur法)、Wood 修正方程法(简称Wood法)、Gassmann法是目前被用来预测水合物饱和度最经典的3种方法，但有各自的适用条件[20]。

水合物在海域沉积物中以3种主要赋存方式[19-21]存在，即悬浮模式、颗粒接触模式及胶结模式。依据胡高伟等[22]研究成果可知，我国神狐海域富集水合物[23]是以悬浮模式赋存的。叠后地震反演可直接反演出波阻抗，不同饱和度水合物会导致波阻抗发生变化，尽管目前Wood法被认为可以较好地预测悬浮模式的水合物饱和度[20]，但并没有给出不同饱和度水合物与波阻抗之间的关系，因此很难直接将其应用到实际的地震资料中来进行水合物饱和度预测。针对这一问题，本文在对Wood法进行推导分析简化的基础上，形成了一种用于悬浮模式水合物饱和度预测的新方法，即Wood阻抗法。

1 基于Wood阻抗法的水合物饱和度计算原理

1.1 Wood法

求取水合物饱和度的Wood法由2个方程构成[17,19]，即

(1)

ρb=(1-Sh)φρw+φShρh+(1-φ)ρm

(2)

式(1)(2)中:Vb、Vpw、Vph、Vpm分别表示含水合物沉积物、海水、纯水合物、基质的纵波速度；φ表示孔隙度；Sh表示水合物在孔隙空间中所占的比例;ρb、ρw、ρh、ρm分别表示含水合物沉积物、海水、纯水合物、基质的密度。

1.2 基质密度和纵波速度计算

一般来说,构成海底沉积物的基质由多种物质组成，其密度表示为

(3)

式(3)中:ρm为基质的密度；fi为第i种矿物组分在基质中的体积百分比；ρi为组成基质中第i种物质的密度;n表示组成基质所包含物质的数量。

基质的纵波速度可表示为[24]

(4)

式(4)中：Km、Gm分别为基质的体积模量和剪切模量。

基质的体积模量和剪切模量可依据 Hill平均公式来求取[25],即

(5)

(6)

式(5)中：Ki为第i种矿物组分的体积模量；Km为基质的体积模量；Gi为第i种矿物组分的剪切模量；Gm为基质的剪切模量。

1.3 纯水合物及海水的纵横波速度计算

依据文献[24]中有关介质弹性参数之间的关系，对于纯水合物的纵横波速度可直接依据弹性参数并利用式(7)、(8)来进行计算。

(7)

(8)

式(7)、(8)中：E为杨氏模量；ρ为物质的密度；σ为泊松比；Vph为纵波速度；Vs为横波速度。

对于杨氏模量和泊松比，可利用物质的体积模量和剪切模量来进行求取,即

(9)

(10)

式(9)、(10)中：K、G分别为物质的体积模量和剪切模量。

对于海水的纵波速度，一般取1 480 m/s;由于海水为液体，横波在其中无法传播，所以海水中的横波速度为0。

1.4 Wood阻抗法推导

将式(1)两边同时乘以1/ρb，则有

(11)

令含水合物沉积物的纵波阻抗Zb=ρbVb,海水的纵波阻抗Zw=ρwVpw,纯水合物的纵波阻抗Zh=ρhVph；基质的纵波阻抗Zm=ρmVpm,则式(11)可表示为

(12)

由文献[11]、[26]可知,海水的密度一般可取1.03 g/cm3,纯水合物的密度一般可取0.9 g/cm3。由式(2)可知

ρb/ρw=(1-Sh)φ+φShρh/ρw+(1-φ)ρm/ρw

≈(1-Sh)φ+0.87φSh+0.97(1-φ)ρm

=0.97(1-φ)ρm+φ-0.13φSh

(13)

ρb/ρh=(1-Sh)φρw/ρh+φSh+(1-φ)ρm/ρh

≈1.15(1-Sh)φ+φSh+1.11(1-φ)ρm

=1.11(1-φ)ρm+1.15φ+0.15φSh

(14)

ρb/ρm=(1-Sh)φρw/ρm+φShρh/ρm+(1-φ)

≈(1-φ)+1.03φ/ρm-0.13φSh/ρm

(15)

令Cbw=ρb/ρw,Cbh=ρb/ρh,Cbm=ρb/ρm,则有Cbh>Cbw>1>Cbm。由于ρb一般大于1.5 g/cm3，而基质密度比ρb大，但最大一般不超过3 g/cm3，所以最小的系数Cbm都大于0.5;而φSh一般在0.1左右，式(13)～(15)中φSh所在项往往在0.01左右，因此相对于大于0.5的数值而言，该数值可忽略不计。鉴于φSh所在项可以忽略，则Cbw、Cbh、Cbm可近似表示为

Cbw≈0.97(1-φ)ρm+φ

(16)

Cbh≈1.11(1-φ)ρm+1.15φ

(17)

Cbm≈(1-φ)+1.03φ/ρm

(18)

由公式(16)～(18)可知，Cbw、Cbh、Cbm可被认为是与孔隙度和基质密度有关的系数。

将系数Cbw、Cbh、Cbm代入式(12)，可得

(19)

式(19)即为本文提出的Wood阻抗法求取水合物饱和度公式。

由此可见，在已知一个地区含水合物沉积物的纵波阻抗、海水的纵波阻抗、纯水合物的纵波阻抗、基质的纵波阻抗、孔隙度以及基质密度的情况下，则可求出水合物饱和度。其中，含水合物沉积物的纵波阻抗可以通过纵波阻抗反演求得，海水的纵波阻抗、纯水合物的纵波阻抗则依据弹性参数测量可以获得，基质的纵波阻抗、孔隙度以及基质密度一般依据一个地区钻井资料可以获得。

2 Wood阻抗法实验误差分析

为验证Wood阻抗法的可靠性，这里以一组数据为例来进行误差分析。在应用方法前先做一些基本的数据假设：

1)假设海域含水合物沉积物的基质由粉砂和黏土2种物质组成，它们在基质中的比例分别是75%、25%；

2)假设悬浮模式水合物饱和度为30%；

3)假设水合物沉积物的孔隙度为40%。

表1列出了一些研究所涉及物质的弹性参数[11,26]。由表1和式(3)可算出基质的密度约为2.63 g/cm3；由表1和式(5)、(6)可算出基质的体积模量和剪切模量分别约为33.94、19.32 GPa；由基质的密度、体积模量和剪切模量，利用式(4)可算出基质的纵波速度约为4762.34 m/s；由表1和式(9)、(10)可算出水合物的杨氏模量和泊松比分别约为6.3 GPa、0.31；由表1和式(7)以及计算的杨氏模量和泊松比可算出水合物的纵波速度为3 126.94 m/s；由式(16)～(18)可算出Cbw、Cbh、Cbm分别约为1.93、2.21、0.76；由式(2)可算出含水合物沉积物的密度约为1.97 g/cm3；由式(1)可算出含水合物沉积物的纵波速度约为1 855.96 m/s；结合已算出的数据可知含水合物沉积物的纵波阻抗Zb、流体的纵波阻抗Zw、纯水合物的纵波阻抗Zp、基质的纵波阻抗Zm约为3 651.61、1 524.4、2 814.25、12 536.88(m·s-1)(g·cm-3)。由式(19)可算出预测的水合物饱和度是28.5%，与实际已知值30%非常接近，说明本文提出的Wood阻抗法具有一定的适用性。

表1 物质的弹性参数Table 1 Elastic parameters of matter

3 应用效果分析

3.1 基础数据准备

以神狐海域数据为参照建立二维模型。陈芳等[27]指出,神狐海域水合物样品附近沉积物的组成按照粒度大小分为砂、粉砂、黏土等3种;沉积物中最主要的颗粒类型是粉砂，其次是黏土，再次是砂，其中粉砂的平均含量大致占75%，黏土的平均含量大致占20%，砂的平均含量大致占5%。本次应用中，参照上述沉积物粒度组成分析成果，设置用于构成沉积物基质的各物质比例是：粉砂75%、黏土20%、砂5%。

图1、2分别是神狐海域钻获水合物的S 井的纵波速度和密度曲线。从图1可以看出，该井纵波速度曲线可大致分为3段:第1段为低速层，深度介于100～198 m，纵波速度一般低于1 900 m/s；第2段为高速层，深度介于 198～218 m，纵波速度明显增加，最高达到 2 300 m/s，已被确认为含水合物层[11]；第3段为中低速度层，深度介于198～245 m，纵波速度相对含水合物层出现了明显的减小，大致在1900 m/s左右。从图2可以看出，该井密度曲线也可大致分为3段：第1段为中低密度段，深度介于100～200 m，密度主要在1.98 g/cm3左右；第2段为低密度段，深度介于200～210 m，密度主要在1.90 g/cm3左右；第3段为高密度段，深度介于210～245 m，密度主要在2.03 g/cm3左右。从图1、2的对比可以看出，水合物高速层并不完全表现为低密度，210 m之上的水合物层密度主要表现为低值，但大于210 m时主要表现为高值。陈芳等[27]对S井开展了沉积物组成成分分析，指出水合物层下段高密度处的砂岩及粗砂岩含量明显高于上部低密度处的砂岩及粗砂岩含量，水合物层下段高密度处的黏土含量要低于上部低密度处的黏土含量，因此推测水合物层表现不同的密度差异性可能主要是岩性差异造成的。尽管水合物层不同位置的密度存在一定的差异性，但一般来说，密度与纵波速度相比对纵波阻抗的贡献要小得多，因此时常可在研究时忽略密度变化对阻抗的影响。

图1 S井纵波速度随深度的变化Fig.1 Variation of P-wave velocity with depth in Well S

图2 S井密度随深度的变化Fig.2 Variation of density with depth in Well S

结合图1，本文将低速段称为上覆非水合物层，高速段称为水合物层，中低速段称为下伏非水合物层，表2列出了S井不同层段密度和速度的平均取值。

表2 S井不同层段密度和速度的平均值Table 2 Average density and velocity in different sections of Well S

3.2 二维地质模型描述

图3 是依据海水弹性参数和表2数据构建的二维地质模型，该模型共包括4层，第1、2、3、4层分别是海水层、上覆非水合物层、水合物层和下伏非水合物层，其纵波速度和密度数据见表3。依据式(1)及表1、2可估算第2层和第4层的非水合物层孔隙度分别约为36%、28%。一般来讲，海底沉积物中的孔隙度往往随深度的增加逐渐减小，模型中水合物层介于上覆非水合物层和下伏非水合物层之间，所以其孔隙度一般应介于这2个非水合物层的孔隙度之间，因此本文取这2个非水合物层的孔隙度平均值32%作为水合物层(即第3层)的孔隙度。

图3 含水合物层的二维地质模型Fig.3 2D geological model of gas hydrate-bearing sediments

表3 含水合物层的二维地质模型的纵波速度和密度数据Table 3 P-wave velocity and density data for 2D gelogical model of gas hydrate-bearing sediments

模拟海上产生零偏移距地震记录，假设震源子波是雷克子波，地震子波主频为50 Hz，则图3 对应的零偏移距记录如图4所示。图5为依据图4振幅特征解释的3个层位，从上到下依次为B1、B2、B3，其中B1代表海底的反射，具有正极性反射特征且振幅最强；B2代表上覆非水合物层与水合物层之间的反射，具有正极性反射特征且振幅相对较强；B3代表水合物与下伏非水合物层的反射，具有负极性反射特征且振幅相对较弱。图6是反演的纵波阻抗剖面，B2和B3层位之间水合物层具有较高的纵波阻抗，在4 100 (m·s-1)(g·cm-3)左右。基于图6，再依据式(19)预测出水合物饱和度(图7)，预测的水合物饱和度大致介于28%～38%。而Zhang Haiqi等[23]利用氯离子法估算水合物含量富集层段水合物含量主要介于30%～40%，可见本文方法预测结果可信，说明该方法具有一定的适用性。