一种应力水平各向异性介质方位地震反演方法及系统

1.本发明涉及地震储层识别领域，尤其涉及一种应力水平各向异性介质方位地震反演方法及系统。


背景技术：

2.从地表地震数据中检测地震各向异性有助于表征地下裂缝发育程度，指导地下裂缝位置预测。然而，初始各向异性弹性介质受到应力作用后，刚度张量的对称性将发生变化，变化后的介质并不完全等同于无应力的各向异性介质。因此，在研究裂缝和应力引起的各向异性对地震波速度、振幅或岩石其他物理性质的影响时，两者的解耦变得格外重要。同时，地应力场的估算也是地下油气勘探开发的关键。目前所采用的地应力估算方法通常是根据地面反射地震数据，利用地震速度和地应力之间的经验近似方程得到的，这种方法往往局限于经验模型的准确性或适用性。对于经验模型不准确或不适用的情况下，以往的地应力估算方法准确性较差。


技术实现要素：

3.针对现有技术存在的不足，本发明的目的是提供一种应力水平各向异性介质方位地震反演方法及系统，利用非线性声弹性理论构建了水平单轴或双轴应力hti介质的有效刚度张量，并对裂缝和应力诱导的各向异性参数进行了解耦，用avaz地震反演方法对未知模型参数进行估计，利用不同方位角的地震反射振幅数据进行反演，能够提高地下裂缝位置预测的准确性。
4.为了实现上述目的，本发明是通过如下的技术方案来实现：
5.第一方面，本发明的实施例提供了一种应力水平各向异性介质方位地震反演方法，包括：
6.建立线性滑移边界条件下的hti介质模型，并对hti介质模型施加平行于对称轴的单轴小应力；
7.对水平单轴应力hti介质中裂缝和应力诱导的刚度张量进行解耦；
8.建立解耦裂缝和应力诱导各向异性参数的pp波avaz近似方程；
9.利用不同方位之间的反射系数振幅差异进行pp波avaz反演，形成应力水平各向异性介质方位地震反演方案。
10.作为进一步的实现方式，利用非线性声弹性理论表示水平单轴应力作用下的有效弹性刚度矩阵；基于小应变假设，将应力各向异性介质中的有效弹性刚度张量进行简化。
11.作为进一步的实现方式，基于弱各向异性背景假设和水平单轴小应力诱导的各向异性假设(忽略了计算中出现的各向异性与应力相关的近似高阶无穷小项)调整水平单轴应力作用下的有效弹性刚度张量的表达式。
12.作为进一步的实现方式，利用渐近射线理论和固定相位法将线性化的pp波avaz反射系数的变化用散射函数表示：
[0013][0014]
其中，δ表示无应力或应力岩石弹性性质的扰动，δρ表示密度的扰动，η
ij
与入射p波和散射p波的慢度矢量和偏振矢量相关；表示水平单轴应力hti介质中的刚度张量扰动，θ表示入射角，表示方位角。
[0015]
作为进一步的实现方式，根据弱无应力背景弹性性质假设、弱无应力hti介质各向异性假设以及hti介质施加水平单轴小应力假设近似表示刚度张量扰动；并根据刚度张量扰动得到解耦的裂缝和应力各向异性参数的pp波avaz近似方程。
[0016]
作为进一步的实现方式，利用卷积模型生成合成方位地震数据：
[0017]
δd＝gm,
[0018]
其中，δd表示不同方位角之间的差，m表示未知的模型参数；
[0019][0020]
w表示小波向量；
[0021]
根据贝叶斯反演框架，合成方位地震数据可以结合先验信息和似然函数进行求解，即未知模型的后验概率分布函数与未知模型参数的先验信息和似然函数的乘积成正比。
[0022]
作为进一步的实现方式，目标函数φ(m)表示为：
[0023][0024]
其中，表示未知模型参数的正则化系数，m
i0
表示未知模型参数的低频初始模型，表示未知模型参数的方差，表示使用数据的方差，δd表示不同方位角之间的差，m表示未知的模型参数。
[0025]
第二方面，本发明的实施例还提供了一种应力水平各向异性介质方位地震反演系统，包括：
[0026]
介质模型构建模块，被配置为：建立线性滑移边界条件下的hti介质模型，并对hti介质模型施加平行于对称轴的单轴应力；
[0027]
解耦模块，被配置为：对水平单轴应力hti介质中裂缝和应力诱导的刚度张量进行解耦；
[0028]
近似方程构建模块，被配置为：建立解耦裂缝和应力诱导各向异性参数的pp波avaz近似方程；
[0029]
反演方案生成模块，被配置为：利用不同方位之间的反射系数振幅差异进行pp波avaz反演，形成应力水平各向异性介质方位地震反演方案。
[0030]
第三方面，本发明的实施例还提供了一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现所述的一种应力水平各向异性介质方位地震反演方法中的步骤。
[0031]
第四方面，本发明的实施例还提供了一种计算机设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行程序时实现所述的一种应
力水平各向异性介质方位地震反演方法中的步骤。
[0032]
本发明的有益效果如下：
[0033]
(1)本发明利用非线性声弹理论来描述应力hti介质中的p波传播特征，研究了应力介质中近水平单轴应力作用下裂缝hti介质界面的弱各向异性近似pp波方位反射系数；利用方位地震振幅差数据反演了解耦的裂缝和应力诱导的各向异性参数，消除了各向同性背景弹性性质的影响。
[0034]
(2)本发明的pp波avaz近似反演方法能够在中等噪声条件下对未知模型参数进行合理估计；该反演方法在反演解耦裂缝和应力诱导的各向异性参数方面的可靠性较高，在地下宽方位地震资料中对于各向异性应力场的估计将具有更广泛的应用。
附图说明
[0035]
构成本发明的一部分的说明书附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。
[0036]
图1是本发明根据一个或多个实施方式的流程图；
[0037]
图2(a)和图2(b)是本发明根据一个或多个实施方式的合成方位地震数据集；
[0038]
图3是本发明根据一个或多个实施方式的基于无噪声合成数据的pp波avaz反演解耦裂缝和应力诱导各向异性参数反演结果；
[0039]
图4是本发明根据一个或多个实施方式的基于噪声合成数据的pp波avaz反演解耦裂缝和应力诱导各向异性参数反演结果；
[0040]
图5(a)是本发明根据一个或多个实施方式的方位角为20
°
、入射角为10
°
的部分叠后地震资料；
[0041]
图5(b)是本发明根据一个或多个实施方式的方位角为20
°
、入射角为20
°
的部分叠后地震资料；
[0042]
图5(c)是本发明根据一个或多个实施方式的方位角为20
°
、入射角为30
°
的部分叠后地震资料；
[0043]
图6(a)是本发明根据一个或多个实施方式的方位角为65
°
、入射角为10
°
的部分叠后地震资料；
[0044]
图6(b)是本发明根据一个或多个实施方式的方位角为65
°
、入射角为20
°
的部分叠后地震资料；
[0045]
图6(c)是本发明根据一个或多个实施方式的方位角为65
°
、入射角为30
°
的部分叠后地震资料；
[0046]
图7(a)是本发明根据一个或多个实施方式的方位角为110
°
、入射角为10
°
的部分叠后地震资料；
[0047]
图7(b)是本发明根据一个或多个实施方式的方位角为110
°
、入射角为20
°
的部分叠后地震资料；
[0048]
图7(c)是本发明根据一个或多个实施方式的方位角为110
°
、入射角为30
°
的部分叠后地震资料；
[0049]
图8(a)是本发明根据一个或多个实施方式的方位角为155
°
、入射角为10
°
的部分叠后地震资料；
[0050]
图8(b)是本发明根据一个或多个实施方式的方位角为155
°
、入射角为20
°
的部分叠后地震资料；
[0051]
图8(c)是本发明根据一个或多个实施方式的方位角为155
°
、入射角为30
°
的部分叠后地震资料；
[0052]
图9(a)是本发明根据一个或多个实施方式的方位角为20
°
到110
°
、入射角为10
°
的方位振幅差地震数据；
[0053]
图9(b)是本发明根据一个或多个实施方式的方位角为20
°
到110
°
、入射角为20
°
的方位振幅差地震数据；
[0054]
图9(c)是本发明根据一个或多个实施方式的方位角为20
°
到110
°
、入射角为30
°
的方位振幅差地震数据；
[0055]
图10(a)是本发明根据一个或多个实施方式的方位角为65
°
到155
°
、入射角为10
°
的方位振幅差地震数据；
[0056]
图10(b)是本发明根据一个或多个实施方式的方位角为65
°
到155
°
、入射角为20
°
的方位振幅差地震数据；
[0057]
图10(c)是本发明根据一个或多个实施方式的方位角为65
°
到155
°
、入射角为30
°
的方位振幅差地震数据；
[0058]
图11(a)是本发明根据一个或多个实施方式的基于pp波avaz反演的裂缝各向异性参数反演法向裂缝弱度结果；
[0059]
图11(b)是本发明根据一个或多个实施方式的基于pp波avaz反演的裂缝各向异性参数反演切向裂缝弱度结果；
[0060]
图12(a)是本发明根据一个或多个实施方式的基于pp波avaz反演的p波应力诱导各向异性参数结果；
[0061]
图12(b)是本发明根据一个或多个实施方式的基于pp波avaz反演的s波应力诱导各向异性参数结果；
[0062]
图13是本发明根据一个或多个实施方式的基于pp波avaz反演的近井地震数据解耦裂缝和应力诱导各向异性参数反演结果。
具体实施方式
[0063]
实施例一：
[0064]
本实施例提供了一种应力水平各向异性介质方位地震反演方法，如图1所示，包括：
[0065]
建立线性滑移边界条件下的hti介质模型，并对hti介质模型施加平行于对称轴的单轴小应力(小于10mpa，约为密西西比时代伯利亚砂岩单轴断裂强度的四分之一)；
[0066]
对水平单轴应力hti介质中裂缝和应力诱导的刚度张量进行解耦；
[0067]
建立解耦裂缝和应力诱导各向异性参数的pp波avaz近似方程；
[0068]
利用不同方位之间的反射系数振幅差异进行pp波avaz反演，形成应力水平各向异性介质方位地震反演方案。
[0069]
具体的，水平单轴应力作用下的有效弹性刚度矩阵可以用非线性声弹性理论来表示，其中包含三阶弹性常数来描述弹性特性与地应力之间的关系，基于小应变假设，应力各
向异性介质中的有效弹性刚度张量可简化为：
[0070]cijkl
≈δ
ik
t
jl
+a
ijkl
+a
ijklmnemn
,(1)
[0071]
其中，
[0072][0073][0074]
其中，δ
ip
是克罗内克符号，t
jl
表示预先存在的主应力；a
ijkl
表示无应力状态下的四阶张量，a
ijklmn
表示无应力状态下六阶张量，用于描述线性和非线性形变；e
mn
表示从无应力状态到形变状态的应变张量；ρ表示密度，αi表示在应力岩石中沿坐标轴传播的p波速度，β
ij
表示在应力岩石中沿坐标轴传播的s波速度。
[0075]
在应力hti介质中的有效弹性刚度张量c
ijkl
是不对称的(c
ijkl
≠c
jikl
，c
ijkl
≠c
ijlk
，当i＝k且j＝l时，c
ijkl
≠c
jilk
)，这与固有各向异性和裂缝诱导的hti介质中不同。
[0076]
由于t
jl
的数值远比a
ijkl
和a
ijklmn
小得多，即|a
ijklmn
|》》|a
ijkl
|》》|t
jl
|。忽略式(1)中的t
jl
，基于弱各向异性背景假设和水平单轴小应力诱导的各向异性假设水平单轴应力作用下的有效弹性刚度张量可进一步表示为：
[0077][0078]
其中，
[0079][0080][0081][0082][0083]c44
≈μ,(9)
[0084]c66
≈μ(1-δ
t
)+μγs.(10)
[0085]
其中，λ表示无应力岩石的第一拉梅常数，μ表示无应力岩石的第二拉梅常数，m＝λ+2μ表示无应力岩石的p波模量；δn表示无应力hti介质的无量纲法向裂缝弱度，δ
t
表示无应力hti介质的无量纲切向裂缝弱度；γ
p
表示与二阶刚度张量a
ij
(即voigt符号表示的四阶张量a
ijk
)和三阶刚度张量a
ijkl
(即voigt符号表示的六阶张量a
ijklmn
)相关的p波应力诱导各向异性参数，γs表示与二阶刚度张量a
ij
(即voigt符号表示的四阶张量a
ijk
)和三阶刚度张量a
ijkl
(即voigt符号表示的六阶张量a
ijklmn
)相关的s波应力诱导各向异性参数：
[0086][0087][0088]
式中，k
p
＝2a
155
/a
33
和ks＝a
456
/a
55
表示与三阶刚度张量和二阶刚度张量之比有关的两个常数，分别用于控制应力诱导的p波和s波各向异性参数；t
11
表示施加在岩石上的水平单轴应力。
[0089]
公式(5)-(10)是在水平单轴应力hti介质中导出的解耦裂缝和应力诱导的刚度张量。此处忽略了公式(4)中由t
jl
表示的与有效弹性刚度张量不对称性相关的近似项。因此，应力四阶刚度张量可以用voigt符号表示为二阶矩阵
[0090]
利用渐近射线理论和固定相位法，将线性化的pp波的avaz反射系数的变化用散射函数来表示：
[0091][0092]
其中，θ表示入射角，表示方位角，ρ表示密度，散射函数可以表示为：
[0093][0094]
式中，符号δ表示无应力或应力岩石弹性性质的扰动，δρ表示密度的扰动，η
ij
与入射p波和散射p波的慢度矢量和偏振矢量相关；为水平单轴应力hti介质中的刚度张量扰动，根据弱无应力背景弹性性质假设、弱无应力hti介质各向异性假设以及hti介质施加水平单轴小应力假设，刚度张量扰动可以近似为：
[0095][0096]
其中，
[0097][0098][0099]
[0100][0101]
δc
44
≈δμ,(20)
[0102]
δc
66
≈δμ-μδδ
t
+μδγs.(21)
[0103]
将公式(16)-(21)代入公式(13)，可以得到解耦的裂缝和应力各向异性参数的pp波avaz近似方程：
[0104][0105]
其中，
[0106][0107][0108][0109][0110][0111][0112][0113]
式中的式中的和分别表示p波和s波速度的平均值以及反射界面上的密度均值；δα，δβ和δρ分别表示纵波速度、横波速度以及纵波密度的扰动值，δδn和δδ
t
分别表示法向裂缝弱度和切向裂缝弱度的扰动值；δγ
p
和δγs分别表示引入的p波和s波应力诱导各向异性参数的扰动值；为方位角，是观测方位角与裂缝法线之差(即)。
[0114]
利用不同方位间的振幅差进行pp波avaz反演，可以消除各向同性背景对pp波方位振幅数据的影响，直接获得解耦的裂缝和应力各向异性参数。当有m个界面，n个入射角，l个方位角时，不同方位角之间的pp波反射系数可以写成矩阵形式：
[0115][0116]
其中，
[0117][0118][0119][0120][0121][0122][0123][0124][0125][0126][0127][0128][0129][0130][0131]
式中，t1和tm分别是第一个和最后一个反射界面，符号t表示矩阵的转置。
[0132]
之后，利用卷积模型以如下矩阵形式生成合成方位地震数据：
[0133]
δd＝gm,
ꢀꢀꢀ
(35)
[0134]
其中，
[0135]
δd＝[wδr
pp
]
mnl
×1,
[0136][0137][0138]
其中，δd表示不同方位角之间的差，w表示小波向量，m表示未知的模型参数。
[0139]
根据贝叶斯反演框架，公式(35)的目标函数可以结合先验信息和似然函数进行求解，即未知模型的后验概率分布函数(pdf)p(m|δd)与未知模型参数的先验pdfp(m)和似然函数p(δd|m)的乘积成正比：
[0140]
p(m|δd)
∝
p(m)p(δd|m),
ꢀꢀꢀ
(37)
[0141]
式中的未知模型参数的先验pdfp(m)通过柯西分布建模给出：
[0142][0143]
似然函数p(δd|m)可表示为带有l2范数约束的高斯分布：
[0144][0145]
其中，和分别表示未知模型参数和使用数据的方差。
[0146]
因此，可以推导出未知模型参数的后验pdfp(m|δd)为：
[0147][0148]
求出的最终目标函数φ(m)可以表示为：
[0149][0150]
本实施例还加入了低频初始模型，进一步将目标函数写为：
[0151][0152]
式中为未知模型参数的正则化系数，m
i0
为未知模型参数的低频初始模型。
[0153]
使用迭代重加权最小二乘(irls)算法求解公式(42)来获得未知模型参数m的合理估计。
[0154]
本实施例利用非线性声弹理论来描述应力hti介质中的p波传播特征，研究了应力介质中近水平单轴应力作用下裂缝hti介质界面的弱各向异性近似pp波方位反射系数；利用方位地震振幅差数据反演了解耦的裂缝和应力诱导的各向异性参数，消除了各向同性背景弹性性质的影响。
[0155]
实施例二：
[0156]
本实施例使用合成的方位角地震数据集来验证实施例一反演方法，将35hz的ricker子波与pp波方位角反射系数卷积生成合成数据，图2(a)为无噪声合成方位地震数据集；然后在无噪声数据中加入高斯噪声产生噪声数据，图2(b)为信噪比(snr)为2的合成方位地震数据集，4个方位角分别为20
°
，65
°
，110
°
，155
°
。然后使用方位差数据来估计解耦的裂缝和应力诱导的各向异性参数，并在20
°
～110
°
和65
°
～155
°
之间分别选择一个振幅差数据，得到最大的方位差振幅。
[0157]
图3和图4分别为使用pp波avaz反演方法的无噪声和有噪声合成数据的裂缝和应力各向异性参数反演结果。在无噪声数据的情况下，所有反演的解耦裂缝和应力诱导各向异性参数与真实模型参数吻合良好。但是，使用有噪声数据反演的解耦裂缝和应力各向异性参数的精度低于无噪声情况，最终反演的未知模型参数仍能满足实际应用的需要。因此，实施例一的反演方法可以为估计中等噪声情况下的解耦裂缝和应力诱导各向异性参数提供一种可行的方法。
[0158]
除此之外，本实施例还使用从中国西南地区获得的真实数据集来验证所提出的反演方法。目标层为含页岩气储层，包含大量发育的高角度裂缝。同时经历多次构造活动，水平构造应力是地震各向异性的来源之一。使用简化的水平单轴应力hti模型来处理方位角pp波地震数据，对其进行振幅保存、去噪、速度分析、偏移等处理，然后进行排序生成方位角叠加地震数据。图5(a)-图5(c)、图6(a)-图6(c)、图7(a)-图7(c)、图8(a)-图8(c)为平均方位角为20
°
、65
°
、110
°
和155
°
的部分叠加地震数据，以及三个入射角分别为10
°
、20
°
和30
°
的叠加数据。分别处理后的方位角地震数据具有良好的信噪比，可很好地用于解耦裂缝和应力诱导各向异性参数的方位角地震反演。图9(a)、图9(b)、图9(c)分别为方位角20
°
和110
°
，入射角10
°
、20
°
和30
°
的方位振幅差地震数据。图10(a)、图10(b)、图10(c)分别为方位角65
°
和155
°
，入射角10
°
、20
°
和30
°
的方位振幅差地震数据。
[0159]
之后使用生成的方位振幅差数据来反演解耦的裂缝和应力诱导的各向异性参数。图11(a)、图11(b)分别为基于pp波avaz反演的裂缝诱导法向弱度参数和切向弱度参数结果。可以看出，位于a井的含气目标储层红色方块与裂缝弱度高值几乎完全吻合。图12(a)、图12(b)分别为基于pp波avaz反演的p波和s波应力各向异性参数结果。相比之下位于a井的含气目标储层应力诱发各向异性参数值较低，有利于形成全方位的拉张压裂和裂缝网络，促进页岩气的储采。
[0160]
为了进一步验证所提反演方法的可靠性，本实施例利用图13所示的近井地震数据，给出了基于pp波avaz反演的裂缝和应力各向异性参数反演结果。结果表明，该反演方法与真实测井数据具有较好的一致性，可以为具有hii对称性的页岩气储层提供合理可靠的解耦裂缝和应力诱发各向异性参数估计。
[0161]
本实施例的pp波avaz近似反演方法能够在中等噪声条件下对未知模型参数进行合理估计；通过在页岩气储层采集的野外数据集，进一步证明了反演方法在反演解耦裂缝和应力诱导的各向异性参数方面的可靠性，在地下宽方位地震资料中对于各向异性应力场的估计将具有更广泛的应用。
[0162]
实施例三：
[0163]
本实施例提供了一种应力水平各向异性介质方位地震反演系统，包括：
[0164]
介质模型构建模块，被配置为：建立线性滑移边界条件下的hti介质模型，并对hti介质模型施加平行于对称轴的单轴应力；
[0165]
解耦模块，被配置为：对水平单轴应力hti介质中裂缝和应力诱导的刚度张量进行解耦；
[0166]
近似方程构建模块，被配置为：建立解耦裂缝和应力诱导各向异性参数的pp波avaz近似方程；
[0167]
反演方案生成模块，被配置为：利用不同方位之间的反射系数振幅差异进行pp波avaz反演，形成应力水平各向异性介质方位地震反演方案。
[0168]
实施例四：
[0169]
本实施例提供了一种计算机可读存储介质，其上存储有计算机程序，该程序被处理器执行时实现实施例一所述的一种应力水平各向异性介质方位地震反演方法中的步骤。
[0170]
实施例五：
[0171]
本实施例提供了一种计算机设备，包括存储器、处理器及存储在存储器上并可在
处理器上运行的计算机程序，所述处理器执行程序时实现实施例一所述的一种应力水平各向异性介质方位地震反演方法中的步骤。
[0172]
以上所述仅为本技术的优选实施例而已，并不用于限制本技术，对于本领域的技术人员来说，本技术可以有各种更改和变化。凡在本技术的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本技术的保护范围之内。