一种非均质水合物储层的地震响应模拟方法与流程

[0001]
本发明属于天然气水合物地震勘探技术领域，具体涉及一种天然气水合物储层地震反射波形响应的模拟方法，特别涉及一种非均质水合物储层的地震响应模拟方法。


背景技术：

[0002]
天然气水合物地层的地震波形识别特征为似海底反射现象(bsr，bottom simulating reflector)，通常出现在海域环境的斜坡及海底沉积地层。但是，并非所有bsr出现的情况都对应天然气水合物储层。因此，在岩石物理模型的基础上建立天然气水合物的地质地球物理模型，应用地震反射理论与数值模拟技术进行天然气水合物地层地震反射波形响应研究具有重要意义。
[0003]
zhang等(2015)在权重方程的基础上，对bsr进行地震ava(amplitudes versus angle，即反射波振幅随入射角变化)反演预测水合物的存在。宋海斌等(2002)基于岩石物性模型，对bsr的ava特征进行分析表明，不同水合物饱和度和游离气饱和度的不同组合表现出形态相似但反射系数值各有不同的特征。petersen等(1993)对水合物地层的剪切波进行avo分析表明，水合物的存在可增加bsr之上沉积地层的剪切模量。carcione等(2001)给出了传播矩阵理论中，地震p波入射时随频率变化的反射、透射系数向量的求解过程。berge等(1999)和doverkin等(1999)的研究表明，当游离气层中气体向上运移并且在合适的温压条件下形成水合物时，游离气层顶界面处于富气相环境，水合物饱和度较高；而远离游离气层顶界面的区域则逐渐向富水相环境过渡，水合物饱和度逐渐降低。
[0004]
但是，至今还未见在岩石物理模型的基础上，建立水合饱和度随地层深度非均匀变化的地质地球物理模型，并基于传播矩阵理论实现频率域岩石物理模型与地震响应计算的无缝连接，计算天然气水合物地层地震反射波形响应的相关报导。


技术实现要素：

[0005]
本发明的目的就在于提供一种非均质水合物储层的地震响应模拟方法，建立一种非均匀水合物饱和度的岩石物理模型，用于计算含水合物地层的速度频散与衰减，并基于传播矩阵理论开发岩石物理与地震反射波模拟一体化技术，计算非均质水合物地层的反射波形响应。
[0006]
本发明的的目的是通过以下技术方案实现的：
[0007]
一种非均质水合物储层的地震响应模拟方法，包括以下步骤：
[0008]
a、针对水合物储层、游离气地层和含水地层，获取基质骨架中各类型矿物组分的含量、弹性模量及密度，水合物的弹性模量及密度，地层水的弹性模量及密度，储层孔隙度；
[0009]
b、建立水合物地层饱和度垂向变化、地层厚度横向变化的非均匀模型，再利用bisq模型(1993)，给定输入参数，得到含天然气水合物砂岩地层的速度与衰减随频率的变化规律；
[0010]
c、利用white斑状饱和模型(1975)计算含水合物地层的下覆游离气地层的速度频
散特征；
[0011]
d、基于传播矩阵理论(carcione，2001)，计算非均质水合物地层和游离气地层的反射波形。
[0012]
进一步地，步骤b，所述非均质水合物砂岩地层的速度与衰减随频率的变化根据bisq模型(1993)计算，纵波速度、衰减系数及逆品质因子表示为：
[0013]
其中，
[0014]
其中，
[0015][0016][0017]
其中，
[0018][0019][0020]
式中，k
dry
和g
dry
分别为岩石骨架的体积模量和剪切模量，为孔隙度，ρ
s
和ρ
f
分别为固相和流体相的密度，ρ
a
为固相和流体相的耦合附加密度，r为喷流长度，k为固相体积模量，k
f
为流体相的体积模量，η为流体黏度，k为渗透率。
[0021]
进一步地，步骤c，所述游离气地层的速度频散变化根据white斑状饱和模型(1975)计算，等效平面波模量的表达式为：
[0022][0023]
其中，
[0024][0025]
[0026]
s
i
表示各层介质在特征单元中所占体积分数，也即气水饱和度；
[0027]
对于每一种不同流体饱和的层介质中，有：
[0028][0029]
是饱和介质i的等效gassmann体积模量，其中μ
dry
是干固体骨架的剪切模量；
[0030][0031]
是等效gassmann体积模量，其中k
dry
是干固体骨架的矿物颗粒的体积模量，而：
[0032][0033][0034]
k
s
是组成固体骨架的矿物颗粒的体积模量，是饱和介质i孔隙中流体相的体积模量；
[0035][0036]
r
i
表示快纵波造成的流体张力与总的法向张力之间的比值；以及，
[0037][0038]
i
i
表示慢纵波的阻抗张量，其中，κ和η
i
分别是固体骨架的渗透率和不同流体的动力学粘滞系数；
[0039][0040]
可以认为是慢纵波的复波数，其中，
[0041][0042]
为等效模量，且：
[0043][0044]
e
dry
为干固体骨架的纵波模量；
[0045]
有了介质的等效平面波模量，则介质中的纵波复速度为：
[0046][0047]
其中，ρ
e
表示特征单元的等效密度：
[0048][0049]
其中，l为特征单元的纵向总长度l＝d
1
+d
2
，ρ
i
表示介质i的密度：
[0050][0051]
其中，为介质i的孔隙度，和分别表示介质i的固体材料密度以及孔隙介质中的流体密度。
[0052]
进一步地，步骤d，所述含水合物地层和游离气地层的反射波形由传播矩阵理论(carcione，2001)计算，对于p波入射，反射、透射系数向量r＝[r
pp
,r
ps
,t
pp
,t
ps
]
t
由下式求解：
[0053][0054]
其中，矩阵a
1
、a
2
分别为与上、下层介质弹性模量有关的传播矩阵；为具有n层结构的中间层的传播矩阵，h
α
为每一层的厚度，对于单一层状介质则有n＝1；i
p
为p波入射向量，与入射介质的物性参数有关；式(16)中，上、下层介质的传播矩阵a
1
和a
2
分别为
[0055][0056][0057]
其中，i为虚数单位，ω为入射波频率，为层的总厚度；变量s
z
、β、γ、w和z有两组下标，其中下标p、s分别对应准纵波、准横波，下标1、2分别对应上、下层介质；
[0058]
传播矩阵(17)和(18)中，s
z
为垂直波慢度，通过水平波慢s
x
计算，去掉下标简写形式为
[0059][0060]
其中，
[0061]
式(19)中，s
z
表达式的两组
±
符号约定为：(+,-)：向下传播qp波(+,+)：向下传播qs波(-,-)：向上传播qp波(-,+)：向上传播qs波水平波慢度s
x
的表达式为
[0062][0063]
并有，e＝{[(c
33-c
55
)cos
2
θ-(c
33-c
55
)sin
2
θ]
2
+(c
13
+c
55
)
2
sin
2
2θ}
1/2
，其中θ为入射角，ρ为相应层位介质密度，c
11
、c
33
、c
55
和c
13
为相应层位各向异性弹性系数；
[0064]
式(16)中，中间层传播矩阵b
α
＝t(0)t-1
(h
α
)(α＝1,...,n)中矩阵t表达式为
[0065][0066]
式(16)中，p波入射向量为i
p
＝iω[β
p1
,γ
p1
,-z
p1
,-w
p1
]
t
。
[0067]
与现有技术相比，本发明的有益效果在于：
[0068]
1、本发明发现含水合物地层中，水合物的饱和度一般随地层深度变化，水合物地层的厚度也存在横向变化，由此，建立水合物地层饱和度垂向变化、地层厚度横向变化的非均匀模型更符合实际地层情况；
[0069]
2、本发明通过结合bisq理论(1993)基础上的岩石物理模型以及传播矩阵理论(carcione，2001)基础上的地震模拟方法，能够计算非均匀水合物地层的地震反射波形响应。
附图说明
[0070]
图1是水合物储层的地质地球物理模型；
[0071]
图2a-图2c是水合物饱和度随深度变化的非均匀模型以及对应的纵波速度vp频散和纵波衰减因子1/qp；
[0072]
图3a-图3b是游离气层纵波速度vp频散和纵波衰减因子1/qp；
[0073]
图4是水合物饱和度垂向变化并且储层厚度横向变化的非均匀模型；
[0074]
图5是水合物非均质模型对应的叠后地震响应。
具体实施方式
[0075]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。
[0076]
实施例1
[0077]
如图1所示，为了实现上述目的，本发明提供了一种非均质水合物储层的地震响应模拟方法，包括以下步骤：
[0078]
a、针对水合物储层、游离气地层和含水地层，获取基质骨架中各类型矿物组分的
含量、弹性模量及密度，水合物的弹性模量及密度，地层水的弹性模量及密度，储层孔隙度；
[0079]
b、建立水合物地层饱和度垂向变化、地层厚度横向变化的非均匀模型，再利用bisq模型(1993)，给定输入参数，得到含天然气水合物砂岩地层的速度与衰减随频率的变化规律；
[0080]
c、利用white斑状饱和模型(1975)计算含水合物地层的下覆游离气地层的速度频散特征；
[0081]
d、基于传播矩阵理论(carcione，2001)，计算非均质水合物地层和游离气地层的反射波形。
[0082]
步骤b，所述含天然气水合物砂岩地层的速度与衰减随频率的变化根据bisq模型(1993)计算，纵波速度、衰减系数及逆品质因子表示为：
[0083]
其中，
[0084]
其中，
[0085][0086][0087]
其中，
[0088][0089][0090]
式中，k
dry
为水合物储层岩石物理模型中干骨架的体积模量，gpa；g
dry
为水合物储层岩石物理模型中干骨架的剪切模量，gpa；k
f
为流体混合物的体积模量，gpa；k
s
为固体基质的体积模量，gpa；为水合物储层岩石物理模型的孔隙度，％；ρ
s
为固体基质的密度，kg/m
3
；ρ
f
为流体混合物的密度，kg/m
3
；ρ
a
为固体基质和流体混合物的耦合附加密度，kg/m
3
(在一具体实施方式中取值420kg/m
3
)；r为喷流长度，mm；η为流体混合物的黏度，pa
·
s；k为渗透率，m
2
；ω为角频率，hz；v
p
为纵波速度，m/s；α为衰减系数；q-1
为逆品质因子(即品质因子q的逆)。
[0091]
步骤c，所述游离气地层的速度频散变化根据white斑状饱和模型(1975)计算，等效平面波模量的表达式为：
[0092][0093]
其中，
[0094][0095][0096]
s
i
表示各层介质在特征单元中所占体积分数，也即气水饱和度。
[0097]
对于每一种不同流体饱和的层介质中，有：
[0098][0099]
是饱和介质i的等效gassmann体积模量，其中，μ
dry
是干固体骨架的剪切模量。
[0100][0101]
是等效gassmann体积模量，其中k
dry
是干固体骨架的矿物颗粒的体积模量，而：
[0102][0103][0104]
k
s
是组成固体骨架的矿物颗粒的体积模量，是饱和介质i孔隙中流体相的体积模量。
[0105][0106]
r
i
表示快纵波造成的流体张力与总的法向张力之间的比值。以及，
[0107][0108]
i
i
表示慢纵波的阻抗张量，其中，κ和η
i
分别是固体骨架的渗透率和不同流体的动力学粘滞系数。
[0109][0110]
可以认为是慢纵波的复波数，其中，
[0111][0112]
为等效模量，且：
[0113][0114]
e
dry
为干固体骨架的纵波模量。
[0115]
有了介质的等效平面波模量，则介质中的纵波复速度为：
[0116][0117]
其中ρ
e
表示特征单元的等效密度：
[0118][0119]
其中l为特征单元的纵向总长度l＝d
1
+d
2
，ρ
i
表示介质i的密度：
[0120][0121]
其中为介质i的孔隙度，和分别表示介质i的固体材料密度以及孔隙介质中的流体密度。
[0122]
步骤d，所述含水合物地层和游离气地层的反射波形由传播矩阵理论(carcione，2001)计算，对于p波入射，反射、透射系数向量r＝[r
pp
,r
ps
,t
pp
,t
ps
]
t
由下式求解：
[0123][0124]
其中，矩阵a
1
、a
2
分别为与上、下层介质弹性模量有关的传播矩阵；b
α
＝t(0)t-1
(h
α
)(α＝1,...,n)为具有n层结构的中间层的传播矩阵，h
α
为每一层的厚度，对于单一层状介质则有n＝1；i
p
为p波入射向量，与入射介质的物性参数有关。方程(16)中，上、下层介质的传播矩阵a
1
和a
2
分别为
[0125][0126][0127]
其中，i为虚数单位，ω为入射波频率，为层的总厚度。变量s
z
、β、γ、w和z有两组下标，其中下标p、s分别对应准纵波、准横波，下标1、2分别对应上、下层介质。
[0128]
传播矩阵(17)和(18)中，s
z
为垂直波慢度，通过水平波慢s
x
计算，去掉下标简写形式为
[0129]
[0130]
其中
[0131]
(19)式中s
z
表达式的两组
±
符号约定为：(+,-)：向下传播qp波(+,+)：向下传播qs波(-,-)：向上传播qp波(-,+)：向上传播qs波水平波慢度s
x
的表达式为
[0132][0133]
并有e＝{[(c
33-c
55
)cos
2
θ-(c
33-c
55
)sin
2
θ]
2
+(c
13
+c
55
)
2
sin
2
2θ}
1/2
，其中θ为入射角，ρ为相应层位介质密度，c
11
、c
33
、c
55
和c
13
为相应层位各向异性弹性系数。
[0134]
方程(16)中，中间层传播矩阵b
α
＝t(0)t-1
(h
α
)(α＝1,...,n)中矩阵t表达式为
[0135][0136]
式(16)中，p波入射向量为i
p
＝iω[β
p1
,γ
p1
,-z
p1
,-w
p1
]
t
。
[0137]
实施例2
[0138]
通过一组数据说明：
[0139]
a、针对水合物储层、游离气地层和含水地层，获取基质骨架中各类型矿物组分的含量、弹性模量及密度，水合物的弹性模量及密度，地层水的弹性模量及密度，储层孔隙度；
[0140]
b、建立水合物地层饱和度垂向变化、地层厚度横向变化的非均匀模型，再利用bisq模型(1993)，给定输入参数，得到含天然气水合物砂岩地层的速度与衰减随频率的变化规律；
[0141]
其中，含天然气水合物砂岩地层的速度与衰减随频率的变化根据bisq模型(1993)计算，纵波速度、衰减系数及逆品质因子表示为：
[0142]
其中，
[0143]
其中，
[0144]
[0145][0146]
其中，
[0147][0148][0149]
式中，k
dry
为水合物储层岩石物理模型中干骨架的体积模量，gpa；g
dry
为水合物储层岩石物理模型中干骨架的剪切模量，gpa；k
f
为流体混合物的体积模量，gpa；k
s
为固体基质的体积模量，gpa；为水合物储层岩石物理模型的孔隙度，％；ρ
s
为固体基质的密度，kg/m
3
；ρ
f
为流体混合物的密度，kg/m
3
；ρ
a
为固体基质和流体混合物的耦合附加密度，kg/m
3
(在一具体实施方式中取值420kg/m
3
)；r为喷流长度，mm；η为流体混合物的黏度，pa
·
s；k为渗透率，m
2
；ω为角频率，hz；v
p
为纵波速度，m/s；α为衰减系数；q-1
为逆品质因子(即品质因子q的逆)。
[0150]
在实施例中，bisq模型(1993)中给定参数特征喷流长度为3mm，流体混合物的粘滞系数为0.03pa*s，渗透率为100md，孔隙度为0.3，含天然气水合物岩石各种矿物体积分数、密度以及弹性模量数据参见表1，水合物、地层水密度以及弹性模量数据参见表2。
[0151]
表1
[0152]
矿物组分密度(kg/m
3
)体积模量(gpa)剪切模量(gpa)体积分数(％)石英2620364580粘土258020.96.8520
[0153]
表2
[0154]
组分密度(kg/m
3
)体积模量(gpa)剪切模量(gpa)水合物9205062.4水10402050
[0155]
水合物地层水合物饱和度随层厚度变化的非均质模型及相应的纵波速度vp频散和纵波衰减因子1/qp如见图2a-图2c所示。
[0156]
步骤c，利用white斑状饱和模型(1975)计算含水合物地层的下覆游离气地层的速度频散特征；
[0157]
其中，游离气地层的速度频散变化根据white斑状饱和模型(1975)计算，等效平面波模量的表达式为：
[0158][0159]
其中，
[0160][0161][0162]
s
i
表示各层介质在特征单元中所占体积分数，也即气水饱和度。
[0163]
对于每一种不同流体饱和的层介质中，有：
[0164][0165]
是饱和介质i的等效gassmann体积模量，其中μ
dry
是干固体骨架的剪切模量。
[0166][0167]
是等效gassmann体积模量，其中k
dry
是干固体骨架的矿物颗粒的体积模量，而：
[0168][0169][0170]
k
s
是组成固体骨架的矿物颗粒的体积模量，k
fi
是饱和介质i孔隙中流体相的体积模量。
[0171][0172]
r
i
表示快纵波造成的流体张力与总的法向张力之间的比值。以及，
[0173][0174]
i
i
表示慢纵波的阻抗张量，其中，κ和η
i
分别是固体骨架的渗透率和不同流体的动力学粘滞系数。
[0175][0176]
可以认为是慢纵波的复波数，其中，
[0177][0178]
为等效模量，且：
[0179][0180]
e
dry
为干固体骨架的纵波模量。
[0181]
有了介质的等效平面波模量，则介质中的纵波复速度为：
[0182][0183]
其中ρ
e
表示特征单元的等效密度：
[0184][0185]
其中l为特征单元的纵向总长度l＝d
1
+d
2
，ρ
i
表示介质i的密度：
[0186][0187]
其中为介质i的孔隙度，和分别表示介质i的固体材料密度以及孔隙介质中的流体密度。
[0188]
在实施例中，游离气地层孔隙中气水饱和度、密度、粘度系数以及弹性模量数据参见表3。
[0189]
表3
[0190][0191]
水合物饱和度垂向变化并且储层厚度横向变化的非均匀模型参见图4。
[0192]
步骤d，基于传播矩阵理论(carcione，2001)，计算非均质水合物地层和游离气地层的反射波形；
[0193]
其中，含水合物地层和游离气地层的反射波形由传播矩阵理论(carcione，2001)计算，对于p波入射，反射、透射系数向量r＝[r
pp
,r
ps
,t
pp
,t
ps
]
t
由下式求解：
[0194][0195]
其中，矩阵a
1
、a
2
分别为与上、下层介质弹性模量有关的传播矩阵；b
α
＝t(0)t-1
(h
α
)(α＝1,...,n)为具有n层结构的中间层的传播矩阵，h
α
为每一层的厚度，对于单一层状介质则有n＝1；i
p
为p波入射向量，与入射介质的物性参数有关。方程(16)中，上、下层介质的传播矩阵a
1
和a
2
分别为
[0196]
[0197][0198]
其中，i为虚数单位，ω为入射波频率，为层的总厚度。变量s
z
、β、γ、w和z有两组下标，其中下标p、s分别对应准纵波、准横波，下标1、2分别对应上、下层介质。
[0199]
传播矩阵(17)和(18)中，s
z
为垂直波慢度，通过水平波慢s
x
计算，去掉下标简写形式为
[0200][0201]
其中，
[0202]
(19)式中s
z
表达式的两组
±
符号约定为：(+,-)：向下传播qp波(+,+)：向下传播qs波(-,-)：向上传播qp波(-,+)：向上传播qs波水平波慢度s
x
的表达式为
[0203][0204]
并有e＝{[(c
33-c
55
)cos
2
θ-(c
33-c
55
)sin
2
θ]
2
+(c
13
+c
55
)
2
sin
2
2θ}
1/2
，其中θ为入射角，ρ为相应层位介质密度，c
11
、c
33
、c
55
和c
13
为相应层位各向异性弹性系数。
[0205]
式(16)中，中间层传播矩阵b
α
＝t(0)t-1
(h
α
)(α＝1,...,n)中矩阵t表达式为
[0206][0207]
式(16)中，p波入射向量为i
p
＝iω[β
p1
,γ
p1
,-z
p1
,-w
p1
]
t
。
[0208]
最终，得到的非均质水合物模型地震avo响应参见图5，游离气层厚度为50m。
[0209]
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。