高精度结构曲率提取方法在潜江凹陷构造解释中的应用

霍志周,刘喜武,莫 莉,钟庆良,张颖燕

(1.中国石油化工股份有限公司石油勘探开发研究院,北京100083;2.中国石油化工股份有限公司江汉油田分公司物探研究院,湖北武汉430034)

几何属性和振幅类属性是地震解释中最为重要的两类地震属性[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)(其中x和y表示空间两个平面坐标,t表示时间坐标),那么每道信号的复地震道c(x,y,t)虚部H[s(x,y,t)](H[s(x,y,t)]在后文简写为sh(x,y,t))可以由Hilbert变换得到:

c(x,y,t)=s(x,y,t)+iH[s(x,y,t)]

=s(x,y,t)+ish(x,y,t)

(1)

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

(2)

(3)

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

(4)

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

(5)

瞬时相位对y和t的偏导数可用类似方法求得。为简化符号,记∂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)做为加权项,可以在分析窗内构建基于相位的梯度结构张量如下:

ST(x,y,t)=∑A2(x,y,t)V(x,y,t)VH(x,y,t)

=∑A2(x,y,t)·

(6)

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

(7)

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

(8)

(9)

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

(10)

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

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

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

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

1.2 结构曲率估计

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

(11)

(12)

其中,

(13)

(14)

(15)

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

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

2 实际应用

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

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

图3 潜江凹陷潜34-10韵律顶部构造

图4 原始数据的目的层位(a)以及沿目的层位切片(b)

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

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

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

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

图7 图6中方框区域的放大显示结果a 商业软件; b 本文方法

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

3 结论

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