石油物探  2020, Vol. 59 Issue (6): 961-969  DOI: 10.3969/j.issn.1000-1441.2020.06.014
0
文章快速检索     高级检索

引用本文 

霍志周, 刘喜武, 莫莉, 等. 高精度结构曲率提取方法在潜江凹陷构造解释中的应用[J]. 石油物探, 2020, 59(6): 961-969. DOI: 10.3969/j.issn.1000-1441.2020.06.014.
HUO Zhizhou, LIU Xiwu, MO Li, et al. Application of high-precision structural curvature extracting on the structural interpretation of Jianghan-Qianjiang Depression[J]. Geophysical Prospecting for Petroleum, 2020, 59(6): 961-969. DOI: 10.3969/j.issn.1000-1441.2020.06.014.

基金项目

国家科技重大专项(2017ZX05049002-004)和国家自然科学基金(41504092, 41774135)共同资助

作者简介

霍志周(1979—), 男, 博士, 高级工程师, 主要从事页岩油气地球物理预测方法研究工作。Email:kboby246@163.com

文章历史

收稿日期:2019-04-29
改回日期:2020-08-05
高精度结构曲率提取方法在潜江凹陷构造解释中的应用
霍志周1 , 刘喜武1 , 莫莉2 , 钟庆良2 , 张颖燕2     
1. 中国石油化工股份有限公司石油勘探开发研究院, 北京 100083;
2. 中国石油化工股份有限公司江汉油田分公司物探研究院, 湖北 武汉 430034
摘要:广泛用于描述地质体几何形态变化的结构曲率属性可通过求取地层视倾角的偏导数获得, 因此地层视倾角估计的精度直接影响结构曲率估计的精度。江汉潜江凹陷工区目的层断层发育、高陡构造较多、地震资料信噪比低, 常规结构曲率提取方法无法有效提取目的层的结构曲率并用于断层及高陡构造等复杂构造带的解释。为此, 采用了一种高精度结构曲率提取方法来提取研究区实际资料的结构曲率。该方法首先利用三维数据的瞬时相位构建梯度结构张量, 并在此基础上估计地层视倾角; 然后采用多窗估计技术来提高地层视倾角的估计精度; 最后以高精度地层视倾角为基础进行地层结构曲率估计(包含最大正曲率和最大负曲率)。将高精度结构曲率提取方法及商业软件曲率模块用于江汉潜江凹陷的实际资料处理, 结果表明:高精度结构曲率提取方法所提取的地层结构曲率在抗噪性以及断层结构等的刻画精度方面明显优于商业软件曲率模块所提取的地层结构曲率, 同时对高陡倾角构造及大断层等构造复杂带的刻画精度更高, 更有利于构造解释。该方法可应用于其它类似地区。
关键词复地震道分析    梯度结构张量    特征分解    多窗分析    视倾角    结构曲率    潜江凹陷    
Application of high-precision structural curvature extracting on the structural interpretation of Jianghan-Qianjiang Depression
HUO Zhizhou1, LIU Xiwu1, MO Li2, ZHONG Qingliang2, ZHANG Yingyan2     
1. Sinopec Petroleum Exploration and Production Research Institute, Beijing 100083, China;
2. Geophysical Research Institute of Jianghan Oilfield, SINOPEC, Wuhan 430034, China
Foundation item: This research is financially supported by the National Science and Technology Major Project(Grant No.2017ZX05049002-004)and the National Natural Science Foundation of China (Grant No.41504092, 41774135)
Abstract: Structural curvature is widely used to characterize geometric changes of the geological target.Usually, the curvature can be obtained by computing the partial derivatives of the apparent seismic dips.Therefore, the precision of apparent dip estimation affects the accuracy of structural curvature extraction.The Jianghan-Qianjiang depression exhibits well-developed faults and steep structures.Additionally, the signal-to-noise ratio of one 3D post-stack dataset in the Jianghan-Qianjiang depression is very low.However, the commonly-used structural curvature extracting method cannot effectively extract the curvature of target horizon, and therefore another method based on high-precision dip estimation was proposed to extract some structural curvature attributes in the Jianghan-Qianjiang depression.First, instantaneous phase was used as the fundamental dataset to reduce the influence of the lateral variation in amplitude.Then, this phase was used to construct a gradient structure tensor (GST) matrix with three non-negative eigenvalues.After that, the seismic volumetric dip was calculated from the dominant eigenvector, and a similarity measure was constructed based on these three eigenvalues.Based on the similarity measure, the dip estimation was improved using multi-window technology.Finally, the structural curvatures were extracted by computing the partial derivatives of the seismic volumetric dip.This method was applied to one 3D post-stack dataset in the Jianghan-Qianjiang depression.The results showed that the extracted structural curvatures had better anti-noise performance and fault characterization than other commercial software.Moreover, the extracted structural curvatures characterized some complex structures such as certain large faults and steep structures, unlike other methods.The results show that the proposed curvature extracting method provides a strong basis in the structural interpretation of the Jianghan-Qianjiang depression and can be extended to other regions with similar structures.
Keywords: complex trace analysis    gradient structure tensor    eigendecomposition    multi-window    apparent dip    structural curvature    Qianjiang Depression    

几何属性和振幅类属性是地震解释中最为重要的两类地震属性[1]。几何属性主要用于增强和显示地震层位的几何形态, 包括地层倾角/方位角、连续性和结构曲率等。倾角/方位角属性体包含重要的地震地层学和地理学信息, 不但可直接用于构造解释, 还可为后续处理及解释提供基础数据, 例如导向滤波[2-3]、相干计算[4]、结构曲率计算[5]等。结构曲率可通过求取视倾角一阶导数获取, 被广泛应用于描述地质体的几何形态变化, 该属性对地层的弯曲、褶皱及断层结构等反应敏感。

目前各种地震属性提取方法已从二维拓展到三维。研究人员发现以地层层位一阶导数为基础的几何属性(包含地层倾角和方位角属性等)可以有效识别相干体等方法观察不到的细小地层结构。由于层位包含丰富的几何信息, 仅利用层位一阶导数还无法充分利用层位信息, 因此可利用以层位二阶导数为基础的曲率属性(结构曲率)来进一步刻画地层结构[6-9]。ROBERTS[10]于2001年详细介绍了结构曲率的基本理论, 同时提出第一代结构曲率分析方法, 即层面曲率属性(surface curvature attribute)的计算流程, 给出了多种类型结构曲率的详细计算公式, 并研究了不同类型结构曲率属性的相关性及其在实际资料中的应用效果。2006年, AL-DOSSARY等[5]利用二阶偏导数与一阶偏导数的关系, 直接利用地震数据体所包含的地层方位信息(倾角/方位角)给出了第二代结构曲率分析方法, 即体曲率属性, 并以此为基础采用分数阶导数在频率域实现了结构曲率属性的多尺度分析。2008年, KLEIN等[11]通过寻找时窗内各道之间最大互相关值的方法来计算三维地震数据体中任意点的结构曲率属性, 开拓了结构曲率属性提取的思路。2009年, GANGULY等[12]利用结构曲率属性进行井位部署和评价, 结果表明钻井成功率要明显高于利用其它地震属性确定的井位。陈学华等[13]于2010年采用小波包对地震数据进行分解, 进而在分解后的地震数据上实现了时间方向上的多尺度分析, 随后通过多尺度体结构曲率属性融合显示, 获得更为丰富的构造细节。CHOPRA等[14-15]创造性地将结构曲率属性与相干属性融合显示, 不但可用于构造识别和解释, 还可用于地震资料预处理以及监控地震资料处理的质量。

由于地层结构曲率可由地层视倾角的偏导数获得, 因此地层倾角估计的结果会影响地层结构曲率估计。目前估计地层倾角的方法很多, MARFURT等[16]于1998年利用基于相似度对三维数据进行倾角扫描以获得地层倾角; BARNES[17-18]利用三维复地震道分析来计算地层倾角; BAKKER等[19]和LUO等[20]利用梯度结构张量和加权梯度结构张量直接估计倾角。作为一种地层倾角的直接估计方法, 梯度结构张量可不经倾角扫描过程来获得高精度地层倾角估计, 因此被广泛用于地层倾角估计及其它地震属性提取。WU等[21]利用有方向性的梯度结构张量来改善倾角估计精度; WANG等[22]在联合复地震道分析、梯度结构张量及多窗分析的基础上给出了一种稳定的地层倾角估计方法; 王震等[23]引入梯度结构张量方法用于刻画断溶体的轮廓, 并对其单一特征值及组合特征值进行断溶体轮廓刻画效果的对比分析; 张尔华等[24]将梯度结构张量用于地层倾角估计, 以此倾角估计结果为约束提出了非线性各向异性扩散滤波器并成功用于三维叠后数据的噪声压制; 李勇等[25]提出了利用改进短时傅里叶变换来获得瞬时相位, 在WANG等[22]的工作基础上给出基于梯度结构张量的多尺度曲率属性估算方法; 王清振等[26]及彭达等[27]利用梯度结构张量构造相干度量, 并将其用于断层及盐丘检测。此外, 其它倾角估计类方法也被广泛用于地层倾角估计及曲率属性提取, 李海山等[28]将平面波分解引入结构曲率属性提取, 即从纵测线和联络测线方向利用平面波解构法估算每一个剖面的倾角, 并在此基础上计算结构曲率属性。

江汉盆地古近纪-新近纪时期为陆相盐湖盆地, 潜江组沉积时期沉积中心位于潜江凹陷, 目的层陡倾角构造及断层等复杂构造发育, 因此该地区的精细断层提取对储层刻画具有重要意义。地层结构曲率属性对于精细描述断层等构造具有明显的优势, 因此在对该工区叠后地震资料进行解释前提取地层结构曲率属性有助于提高后续人工解释的可靠性及准确度。采用常规商业软件结构曲率提取模块对该工区叠后资料进行处理, 其结果不但受噪声影响较大, 同时所提取出的结构曲率受陡倾角构造及大断层等复杂构造的屏蔽, 不能有效地刻画陡倾角及大断层等复杂构造附近区域的结构。为此, 本文以WANG等[22]提出的高精度地层倾角估计方法为基础, 针对潜江凹陷三维叠后地震资料, 进行结构曲率(包含最大正曲率和最大负曲率)估计, 并将其用于刻画陡倾角及大断层等复杂构造区域的断层结构等。

1 方法技术

计算层位数据的一阶导数可获得地层的视倾角属性, 而地层结构曲率属性作为层位的二阶导数类属性, 可通过进一步计算层位一阶导数类属性(视倾角等)的导数来获得。利用地层倾角计算地层结构曲率属性时, 地层倾角的精度及稳定性直接影响地层结构曲率属性提取的精度及稳定性。梯度结构张量类方法避免了极为耗时的倾角扫描过程, 可通过矩阵特征分解直接估计出地层倾角。但由于常规梯度结构张量类方法大量进行地震数据的偏导数等运算, 不但振幅的横向变化对地层倾角有显著的影响, 而且会放大噪声的影响。此外, 由于构建梯度协方差矩阵时有分析窗的存在, 断层等细小结构会被模糊。采用相位数据作为基础数据可降低振幅横向变化对倾角估计的影响, 且采用多窗分析技术在改善断层模糊现象的同时也可提高计算稳定性。因此, 在相位数据基础上联合梯度结构张量分析方法及多窗分析技术来估计地层倾角, 可综合各种倾角估计方法的优势, 获得相比于单项技术更为稳定且精度较高的倾角估计结果。

1.1 地层倾角估计

首先基于梯度结构张量和多窗分析实现地层倾角估计。假设三维地震资料为s(x, y, t)(其中xy表示空间两个平面坐标, t表示时间坐标), 那么每道信号的复地震道c(x, y, t)虚部H[s(x, y, t)](H[s(x, y, t)]在后文简写为sh(x, y, t))可以由Hilbert变换得到:

$ \begin{array}{*{20}{l}} {c(x,y,t) = s(x,y,t) + {\rm{i}}H[s(x,y,t)]}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = s(x,y,t) + {\rm{i}}{s_h}(x,y,t)} \end{array} $ (1)

瞬时振幅A(x, y, t)及瞬时相位P(x, y, t)可由下面的公式得到:

$ A(x,y,t) = \sqrt {{s^2}(x,y,t) + s_h^2(x,y,t)} $ (2)
$ P(x,y,t) = {\rm{arctan}}\left[ {\frac{{{s_h}(x,y,t)}}{{s(x,y,t)}}} \right] $ (3)

基于相位数据P(x, y, t)计算相位梯度V(x, y, t):

$ \begin{array}{*{20}{l}} {\mathit{\boldsymbol{V}}(x,y,t) = \left[ {\frac{{\partial P(x,y,t)}}{{\partial x}},\frac{{\partial P(x,y,t)}}{{\partial y}},} \right.}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {{\left. {\frac{{\partial P(x,y,t)}}{{\partial t}}} \right]}^{\rm{T}}}} \end{array} $ (4)

由于相位存在缠绕问题, 对相位直接计算偏导数会导致不稳定现象出现, 但可通过计算(3)式右端项的偏导数来避免相位缠绕问题, 其中瞬时相位对x的偏导数如下:

$ \begin{array}{l} \frac{{\partial P(x,y,t)}}{{\partial x}} = \frac{{\partial {\rm{arctan}}\left[ {\frac{{{s_h}(x,y,t)}}{{s(x,y,t)}}} \right]}}{{\partial x}}\\ = \frac{{s(x,y,t)\frac{{\partial {s_h}(x,y,t)}}{{\partial x}} - {s_h}(x,y,t)\frac{{\partial s(x,y,t)}}{{\partial x}}}}{{{A^2}(x,y,t)}} \end{array} $ (5)

瞬时相位对yt的偏导数可用类似方法求得。为简化符号, 记∂P/∂x=∂P(x, y, t)/∂x, ∂P/∂y=∂P(x, y, t)/∂y, ∂P/∂t=∂P(x, y, t)/∂t。考虑到振幅类属性的稳定性较高, 在振幅较大处具有较高的信噪比(对随机噪声而言), 因此引入A(x, y, t)做为加权项, 可以在分析窗内构建基于相位的梯度结构张量如下:

$ \begin{array}{l} \begin{array}{*{20}{l}} {{\bf{ST}}(x,y,t) = \sum {{A^2}} (x,y,t)V(x,y,t){\mathit{\boldsymbol{V}}^{\rm{H}}}(x,y,t)}\\ {{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} = \sum {{A^2}} (x,y,t) \cdot } \end{array}\\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \left[ {\begin{array}{*{20}{l}} {{{\left( {\frac{{\partial P}}{{\partial x}}} \right)}^2}}&{\frac{{\partial P}}{{\partial x}}\frac{{\partial P}}{{\partial y}}}&{\frac{{\partial P}}{{\partial x}}\frac{{\partial P}}{{\partial t}}}\\ {\frac{{\partial P}}{{\partial y}}\frac{{\partial P}}{{\partial x}}}&{{{\left( {\frac{{\partial P}}{{\partial y}}} \right)}^2}}&{\frac{{\partial P}}{{\partial y}}\frac{{\partial P}}{{\partial t}}}\\ {\frac{{\partial P}}{{\partial t}}\frac{{\partial P}}{{\partial x}}}&{\frac{{\partial P}}{{\partial t}}\frac{{\partial P}}{{\partial y}}}&{{{\left( {\frac{{\partial P}}{{\partial t}}} \right)}^2}} \end{array}} \right] \end{array} $ (6)

上述矩阵为对称半正定的, 因此其所有特征值均为非负的。对对称半正定矩阵进行特征分解:

$ \begin{array}{l} {\bf{ST}}(x,y,t) = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{v}}_1}}&{{\mathit{\boldsymbol{v}}_2}}&{{\mathit{\boldsymbol{v}}_3}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{u_1}}&0&0\\ 0&{{u_2}}&0\\ 0&0&{{u_3}} \end{array}} \right] \cdot \\ {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\left[ {\begin{array}{*{20}{l}} {{\mathit{\boldsymbol{v}}_1}}&{{\mathit{\boldsymbol{v}}_2}}&{{\mathit{\boldsymbol{v}}_3}} \end{array}} \right]^{\rm{T}}} \end{array} $ (7)

其中u1, u2u3ST(x, y, t)的特征值, 且满足u1u2u3, v1, v2v3为对应的特征向量。u1v1分别为主特征值及主特征向量, 则在x方向和y方向的视倾角px(x, y, t)和py(x, y, t)可由v1估计得到:

$ {{p_x}(x,y,t) = \frac{{{v_{1x}}(x,y,t)}}{{{v_{1t}}(x,y,t)}}} $ (8)
$ {{p_y}(x,y,t) = \frac{{{v_{1y}}(x,y,t)}}{{{v_{1t}}(x,y,t)}}} $ (9)

利用上述主分量分析方法可稳定地估计地层倾角, 提高地层倾角的估计精度, 为后续处理(如计算结构曲率及振幅曲率)提供可靠的倾角数据体。同时基于上述矩阵特征分解的结果可以构建相似度量:

$ C(x,y,t) = \frac{{{u_1} - {u_2}}}{{{u_1} + {u_2}}} $ (10)

当分析窗内的反射波形完全一致时, 上述相似度量为1;当分析窗内的反射波形不一致时, 上述相似度量小于1。

尽管利用(7)式和(8)式可以精确估计出地层倾角, 但当分析窗内包含有断层时, 由于空间分析窗的存在, 断层会被模糊。在这种情况下, 估计得到的倾角不能精确地反映分析窗内的真实反射倾角。上述断层模糊现象可以通过多窗分析技术[29]来改善。在当前分析点下, 除了选取一个以分析点为中心的分析窗外, 还可以选取很多包含有分析点的非中心点相邻分析窗。若选取的分析窗大小为3×3×3, 共有1个中心点分析窗及26个非中心点分析窗, 如图 1所示, 黑点表示当前分析点, 图 1a表示1个中心点分析窗, 图 1b表示3个典型的非中心点分析窗。在众多分析窗内均利用(10)式计算相似度量。由于所有的分析窗都包含有当前分析点, 则选择相似度量最大的分析窗作为最终使用的分析窗, 在该分析窗内计算地层沿x方向和y方向的视倾角, 可大大减少断层的模糊现象。

图 1 分析窗大小为3×3×3时的中心点分析窗(a)和3个典型的非中心点分析窗(b)

图 2a为一个二维合成地震剖面, 包含有振幅的横向变化和断层。图 2b为基于原始数据构建梯度结构张量得到的倾角估计, 虽然能够刻画地层倾角, 但是左侧振幅横向变化对倾角估计结果影响很大, 同时断层也被模糊; 而图 2c为本文所使用的基于瞬时相位构建梯度结构张量并结合多窗分析技术得到的倾角估计结果, 可以看出, 左侧振幅横向变化几乎没有影响地层倾角的估计, 同时断层模糊现象也得到了明显改善。

图 2 基于原始合成数据及不同方法估计的地层倾角 a二维合成地震剖面; b基于原始数据构建梯度结构张量得到的倾角估计; c基于瞬时相位构建梯度结构张量并结合多窗分析得到的倾角估计
1.2 结构曲率估计

采用前述方法获得地层倾角(x方向的视倾角及y方向的视倾角)的高精度估计, 然后利用x方向视倾角px(x, y, t)及y方向视倾角py(x, y, t)的偏导数来获得地层结构曲率。其中, 在断层及裂缝解释中常用的两个地层结构曲率(最大正曲率kpos(x, y, t)及最大负曲率kneg(x, y, t))可以通过如下公式求取:

$ \begin{array}{*{20}{l}} {{k_{{\rm{ pos }}}}(x,y,t) = [a(x,y,t) + b(x,y,t)] + }\\ {\sqrt {{{[a(x,y,t) - b(x,y,t)]}^2} + {c^2}(x,y,t)} } \end{array} $ (11)
$ \begin{array}{*{20}{l}} {{k_{{\rm{ neg }}}}(x,y,t) = [a(x,y,t) + b(x,y,t)] - }\\ {\sqrt {{{[a(x,y,t) - b(x,y,t)]}^2} + {c^2}(x,y,t)} } \end{array} $ (12)

其中,

$ {a(x,y,t) = 0.5 \times \frac{{\partial {p_x}(x,y,t)}}{{\partial x}}} $ (13)
$ {b(x,y,t) = 0.5 \times \frac{{\partial {p_y}(x,y,t)}}{{\partial y}}} $ (14)
$ c(x,y,t) = 0.5 \times \left[ {\frac{{\partial {p_x}(x,y,t)}}{{\partial y}} + \frac{{\partial {p_y}(x,y,t)}}{{\partial x}}} \right] $ (15)

其它类型的结构曲率也可通过x方向视倾角px(x, y, t)及y方向视倾角py(x, y, t)的偏导数来获取[5]

2 实际应用

江汉盆地古近纪-新近纪时期为陆相盐湖盆地, 潜江组沉积时期沉积中心位于潜江凹陷, 面积约为2500km2, 陡倾角构造及断层发育。图 3为潜江凹陷潜34-10韵律顶部构造图。潜江凹陷为受北东向潜北大断裂及通海口大断裂所夹持的双断型箕状凹陷, 构造总体上表现为“一凹两斜坡”的基本格局, 地层向东、西、南抬升, 厚度逐渐变小, 中部发育的多条北东向二级断层(浩口、周矶断层等)使构造复杂化。我们选取图中红色方框内潜江凹陷王广浩区域约330km2三维叠后数据进行处理。工区内发育大量的北东向正断层, 倾向北西向, 主要发育在西部斜坡带, 同时也有少数南倾的断层。图 4a为工区的目的层位, 可见陡倾角构造(区域1)及断层(区域2)等复杂构造发育。原始数据信噪比较低, 如图 4b所示的原始数据沿目的层位切片。

图 3 潜江凹陷潜34-10韵律顶部构造
图 4 原始数据的目的层位(a)以及沿目的层位切片(b)

由于该地区地震资料信噪比低, 因此先采用常规方法对地震数据进行去噪预处理以提高信噪比, 然后利用商业软件结构曲率模块和本文的高精度结构曲率提取方法对去噪后数据体分别计算结构曲率并提取沿层切片。在利用本文方法提取地层倾角时, 考虑到对小断层刻画的需要, 空间分析窗大小选为3×3, 而时间分析窗长选为1.5倍数据主周期; 利用多窗分析提高倾角估计精度和稳定性时, 采用27个相邻分析窗(包含1个中心分析窗和26非中心分析窗)。图 5a为利用商业软件结构曲率模块计算得到最大负曲率属性沿T34层位切片。虽然在预处理环节采用了本征图像滤波方法来压制噪声以提高稳定性, 但由于常规方法在计算地层倾角及结构曲率时涉及大量求导运算, 导致计算出的最大负曲率仍受噪声影响比较严重, 如图 5a红色椭圆及红色方框区域所示, 受噪声影响区域基本观察不到任何有效结构。最大负曲率受陡倾角构造的影响, 导致陡倾角及大断层等复杂构造带的结构不清晰, 如图 5a黄色椭圆区域内结构模糊。本文方法提取出的最大负曲率属性(图 5b)的稳定性得到提升, 抗噪性改善明显, 如图 5b红色椭圆及红色方框区域所示。此外, 由于在倾角估计时采用了多窗分析技术, 因此能够明显改善陡构造及大断层等复杂构造带的刻画, 如图 5b黄色椭圆区域所示, 这为后续的陡构造及大断层等复杂构造带的结构解释以及复杂构造带附近的小断层精细解释提供了更为可靠的信息。

图 5 商业软件(a)和本文方法(b)提取的T34层最大负曲率属性沿层切片

图 6a图 6b分别为利用商业软件结构曲率模块和本文方法计算得到的最大正曲率属性沿T34层位切片。对比可以看出, 本文方法得到的最大正曲率在抗噪性及抗高陡构造屏蔽作用方面均有明显的改善。图 7图 6中方框区域的放大显示, 更加清晰地展示出本文方法对陡构造、断层等复杂构造带的刻画精度明显优于常规商业软件结构曲率模块。

图 6 利用商业软件(a)和本文方法(b)提取的T34层最大正曲率属性沿层切片
图 7 图 6中方框区域的放大显示结果 a商业软件; b本文方法

图 8为构造图与利用本文方法得到的最大正曲率沿层切片的叠合显示, 可以看出, 采用本文方法得到的最大正曲率与传统曲率相比更加直观, 同时图 8提供了比构造图更为精细的断裂信息, 更加有利于解释断层及裂缝。

图 8 构造图与采用本文方法得到的最大正曲率沿层切片的叠合显示
3 结论

为满足江汉潜江凹陷工区油气勘探对复杂构造带的精细预测需求, 本文在已有高精度地层倾角估计的基础上, 提取地层的结构曲率属性, 并将其用于工区复杂构造带解释。实际资料处理结果表明:本文方法不但能够提高抗噪性能及改善高陡构造及大断层对周围断层的屏蔽作用, 而且大幅提升了大断层识别精度及小断层的识别能力, 效果明显优于商业软件中的常规方法。此外, 利用本文方法获得了比构造图更为精细的断裂信息, 能够清晰地刻画复杂构造带及断裂带, 为精细构造解释及裂缝预测提供更加可信的结构属性。

参考文献
[1]
TANER M T, KOEHLER F, SHERIFF R E. Complex seismic trace analysis[J]. Geophysics, 1979, 44(6): 1041-1063. DOI:10.1190/1.1440994
[2]
FEHMERS G C, HOCKER C F W. Fast structural interpretation with structure-oriented filtering[J]. Geophysics, 2003, 68(4): 1286-1293.
[3]
LAVIALLE O, POP S, GERMAIN C, et al. Seismic fault preserving diffusion[J]. Journal of Applied Geophysics, 2007, 61(2): 132-141. DOI:10.1016/j.jappgeo.2006.06.002
[4]
MARFURT K J, SUDHAKER V, GERSZTENKORN A, et al. Coherency calculations in the presence of structural dip[J]. Geophysics, 1999, 64(1): 104-111.
[5]
AL-DOSSARY S, MARFURT K J. 3D volumetric multispectral estimates of reflector curvature and rotation[J]. Geophysics, 2006, 71(5): P41-P51. DOI:10.1190/1.2242449
[6]
LISLE R J. Detection of zones of abnormal strains in structures using Gaussian curvature analysis[J]. AAPG Bulletin, 1994, 78(12): 1811-1819.
[7]
ERICSSON J B, MCKEAN H C, HOOPER R J. Facies and curvature controlled 3D fracture models in a Cretaceous carbonate reservoir, Arabian Gulf[J]. Geological Society, London, Special Publications, 1998, 147(1): 299-312. DOI:10.1144/GSL.SP.1998.147.01.20
[8]
HART B S, PEARSON R A, RAWLING G C. 3-D seismic horizon-based approaches to fracture-swarm sweet spot definition in tight-gas reservoirs[J]. The Leading Edge, 2002, 21(1): 28-35.
[9]
SIGISMONDI M E, SOLDO J C. Curvature attributes and seismic interpretation:Case studies from Argentina basins[J]. The Leading Edge, 2003, 22(11): 1122-1126. DOI:10.1190/1.1634916
[10]
ROBERTS A. Curvature attributes and their application to 3D interpreted horizons[J]. First Break, 2001, 19(2): 85-99.
[11]
KLEIN P, RICHARD L, JAMES H. 3D curvature attributes:A new approach for seismic interpretation[J]. First Break, 2008, 26(4): 105-112.
[12]
GANGULY N, DEARBORN D, MOORE M, et al. Application of seismic curvature attribute in the appraisal of the Tishrine-West field, North-East Syria[J]. CSEG Recorder, 2009, 34(6): 29-43.
[13]
CHEN X H, YANG W, HE Z H, et al. The algorithm of 3D multi-scale volumetric curvature and its application[J]. Applied Geophysics, 2012, 9(1): 65-72.
[14]
CHOPRA S, MARFURT K J. Integration of coherence and volumetric curvature images[J]. The Leading Edge, 2010, 29(9): 1092-1107. DOI:10.1190/1.3485770
[15]
CHOPRA S, MARFURT K J. Coherence and curvature attributes on preconditioned seismic data[J]. The Leading Edge, 2011, 30(4): 386-393.
[16]
MARFURT K J, KIRLIN R L, FARMER S L, et al. 3-D seismic attributes using a semblance-based coherency algorithm[J]. Geophysics, 1998, 63(4): 1150-1165.
[17]
BARNES A E. Theory of 2-D complex seismic trace analysis[J]. Geophysics, 1996, 61(1): 264-272.
[18]
BARNES A E. Weighted average seismic attributes[J]. Geophysics, 2000, 65(1): 275-285.
[19]
BAKKER P, VAN VLIET L J, VERBEEK P W. Edge-preserving orientation adaptive filtering[J]. Proceedings of IEEE-CS Conference on Computer Vision and Pattern Recognition, 1999, 535-540.
[20]
LUO Y, HIGGS W G, KOWALIK W S. Edge detection and stratigraphic analysis using 3D seismic data[J]. Expanded Abstracts of 66th Annual Internat SEG Mtg, 1996, 324-327.
[21]
WU X M, JANSON X. Directional structure tensors in estimating seismic structural and stratigraphic orientations[J]. Geophysical Journal International, 2017, 210(1): 534-548. DOI:10.1093/gji/ggx194
[22]
WANG X K, CHEN W C, ZHU Z Y. Robust seismic volumetric dip estimation combining structure tensor and multi-window technology[J]. IEEE Transaction on Geosciences and Remote Sensing, 2019, 57(1): 395-405. DOI:10.1109/TGRS.2018.2854777
[23]
王震, 文欢, 邓光校, 等. 塔河油田碳酸盐岩断溶体刻画技术研究与应用[J]. 石油物探, 2019, 58(1): 149-154.
WANG Z, WEN H, DENG G X, et al. Fault-karst characterization technology in the Tahe Oilfield, China[J]. Geophysical Prospecting for Petroleum, 2019, 58(1): 149-154.
[24]
张尔华, 王伟, 高静怀, 等. 非线性各向异性扩散滤波器用于三维地震资料噪声衰减与结构特征增强[J]. 地球物理学进展, 2010, 25(3): 866-870.
ZHANG E H, WANG W, GAO J H, et al. Non-linear anisotropic diffusion filtering for 3D seismic noise removal and structure enhancement[J]. Progress in Geophysics, 2010, 25(3): 866-870.
[25]
李勇, 张固澜, 何承杰, 等. 基于高精度时频瞬时相位谱的多尺度曲率及其应用[J]. 石油物探, 2019, 58(2): 253-261.
LI Y, ZHANG G L, HE C J, et al. Multiscale curvature via high-precision time-frequency instantaneous phase spectrum and its application[J]. Geophysical Prospecting for Petroleum, 2019, 58(2): 253-261.
[26]
王清振, 张金淼, 姜秀娣, 等. 利用梯度结构张量检测盐丘与断层[J]. 石油地球物理勘探, 2018, 53(4): 826-831.
WANG Q Z, ZHANG J M, JIANG X D, et al. Salt dome and fault detection based on the gradient-structure tensor[J]. Oil Geophysical Prospecting, 2018, 53(4): 826-831.
[27]
彭达, 肖富森, 冉崎, 等. 基于倾角导向梯度能量熵的断层检测方法[J]. 石油地球物理勘探, 2019, 54(1): 191-197.
PENG D, XIAO F S, RAN Q, et al. Fault identification based on dip-oriented gradient-energy-entropy coherence estimation[J]. Oil Geophysical Prospecting, 2019, 54(1): 191-197.
[28]
李海山, 杨午阳. 基于平面波解构的三维体曲率计算方法[J]. 石油物探, 2018, 57(6): 884-891.
LI H S, YANG W Y. Computation of 3D volumetric curvature based on plane-wave destruction[J]. Geophysical Prospecting for Petroleum, 2018, 57(6): 884-891.
[29]
MARFURT K J. Robust estimates of 3D reflector dip and azimuth[J]. Geophysics, 2016, 71(4): P29-P41.