CN1170170C

CN1170170C - 用于地震分析的频谱解译

Info

Publication number: CN1170170C
Application number: CNB971920958A
Authority: CN
Inventors: 格雷戈里A・帕蒂卡; 格雷戈里A·帕蒂卡; M・格里德利; 詹姆斯M·格里德利
Original assignee: BP Corp North America Inc
Current assignee: BP Corp North America Inc
Priority date: 1996-12-06
Filing date: 1997-11-12
Publication date: 2004-10-06
Anticipated expiration: 2017-11-12
Also published as: WO1998025161A3; US5870691A; CN1210591A; EG21105A; NO328897B1; EP0882245A2; RU2187828C2; WO1998025161A2; NO983523D0; NO983523L; AR011288A1

Abstract

本发明涉及地震数据的处理方法，它提供对细微的薄层调谐效应震波和其他种类水平岩层间断情况的改进的量化和显化指示。来自薄层的反射具有指示层位厚度的频率范围上的特征显示：反射在其频谱的振幅上具有周期性的曲线下凹排列，并被与薄层厚度成反比例的空间距离相间隔。该特征显示通过三维量值用于薄层反射道及对厚度和水平延伸进行评价。本发明通过一种从频谱范围白噪化在频谱中来展示突出地质信息的新颖方法被增强。

Description

用于地震分析的频谱解译

技术领域

本发明一般涉及一种对细微地震引起的薄层调谐效应量化和显化的方法。在此公开的方法是通过将地震反射信号以不连续的傅里叶(或其他规范)方法，将取决于待测的薄层厚度的长度进行变换，分解为频率分量来实现的。通过所述不连续的转换进行分解之后，所包含的因素被条理化并以展现及增强表现薄层反射情况的特征频率范围的形式显示，使得其他方法无法观测到的地下层位厚度的变化显示出来。本发明使地震解译者分析和绘制作为空间位置，传播时间和频率的函数的地下地质及地层特征达到一种前所未有的程度。

背景技术

由于大多数规范的地球物理探测是一门相对年轻的科学，该领域一些最早的工作可追溯到1920年，以及有名的CMP大致日期仅在1950年。然而，从其创始年起，它已成为石油工业寻找地下石油沉积卓著的方式。尽管地球物理勘探通常被重力、磁力和地震三大领域所围绕，现今，却是地震法排除其他方式占据绝对支配地位。事实上，简单统计区域中地震队的数量已成为可接受的估量整个石油工业状况的方式。

一种地震测量，体现为试图将地球表面下的地层通过将声能送入地下，并接收由地下岩层返回的回声进行描绘。下送声能的声源，可以是爆炸，或是地上的地震震动器，以及水下环境的气锤。地震测量中，能量移动穿过待研究的层地质结构上面的地层。每当能源引爆，由广泛分布于地面上的许多位置进行记录。众多的引爆/记录组合便结合产生出连续的可延伸数英里的地下剖面。在二维(2-D)地震测量中，记录位置通常沿直线设置，有别于此，三维(3-D)测量中，记录位置呈栅格状分布于整个面上。简单概念上，二维地震线可以想象为给出直接位于记录位置之下的地层横截面图像。三维测量产生一个数据立方体或容量，即至少在概念上，产生一个位于地层测绘面积之下的三维图像。显然可以很容易从三维数值中提取二维测线。

地震测量包括极大量单独的地震记录或地震道。在典型的二维测量中，通常有数万个道，而在三维测量中，单独的道数量将达到数百万个。地震道是声能从地下不均匀反射的数字记录。部分反射不时发生由于地下物质的声学阻抗产生的变化。数字取样通常达到0.004秒(4毫秒)的间隔，尽管亦通用2毫秒和1毫秒的间隔。因此，地震道上的各取样点与传播时间相结合，能量反射情形下具有两路传播时间。再有，地震测量中各道的表面位置被仔细地记录，以及其本身通常作为道的一部分(作为部分的道引导信息)。这使得地震信息包含于与最后相关于特定地下层位的道之中，由此提供了一种地震数据在图上定位和显示轮廓的方式，并从其中提取特征亦即测绘。送入地下的信号称为地震波形或地震子波。不同地震波形的产生取决于震源是气锤、炸药还是振动器。“特征波源”或“脉冲源”的概念通常用于描述已记录的特殊地震波形的地震特征。

地面产生的震源从出原点立即向外向下传播，之后遇到并穿过下面的岩体部分。两类不同的岩体之间的界面上具有位势，可发生震波的反射，在界面上反射回的地震能量大小取决于两部分岩体之间声阻抗的对比，确定这个对比值的常用测量方法是反射系数。反射系数可看成反射波与入射波的振幅之比。按岩体的性质，

其中，岩体的声阻定义为岩石密度(ρ₁和ρ₂分别为上、下岩体的密度)乘以波在岩体中的速度，V₁和V₂为上、下岩体中对应的速度(严格说，上式只在波从岩石界面垂直入射岩石界面时才精确。但在工程上，如果入射波到达界面时，距离垂直方向偏转不超过20°，通常也算满足垂直的条件。)

在地表记录到的反射能量可用下层反射函数表达为震波的褶积，该函数即所谓“褶积模型”。简而言之，褶积模型试图将地表记录到的震动信号描述为下行的原波的褶积，这里要用到反射函数，该函数代表了地下不同岩层之间界面处的反射系数。用公式表示为：

x(t)＝w(t)*e(t)+n(t)

其中，x(t)是记录到的地震波曲线，w(t)是震源波，e(t)是地层的反射函数，n(t)是随机环境噪声，“*”表示褶积运算。此外，褶积模型还要求如下条件：(1)震源波在地下运行时保持不变(即保持恒定)(2)地表记录到的地震道可以表示为各个界面均在地下的源波的若干单个褶积的算术和(叠加法则。即波的反射和传播是线性系统。)尽管人们普遍认为褶积模型不能完整地描述波动的机械特性，但这个模型在多数场合下是足够精确的，这样该模型就非常实用。在地球物探学会Ozdogan Yilmaz的《地震数据处理》第2.2章中，讨论到褶积模型的一些性质，现将这部分内容引入本文，此供参考。

对于石油公司寻找钻井位置来说，经过正确采集、处理的地震数据，能够为勘探人员提供大量的信息。例如，一份震动剖面图为勘探人员提供了地下岩层结构的总体情况，显示出地下岩层与碳氢化合物的蕴藏情况相关的重要特征，如断层，褶皱，背斜，不整合，地下盐丘及矿脉等许多特征。在对震动数据进行计算机处理过程中，可由程序自动对波速进行估算，临近地表的非均匀速度也能检测显示出。有时，可直接应用震动数据直接估算岩石的孔隙率，水饱和性和碳氢化合物组分。一些不太直观的震波特征如相位、波峰振幅、波峰至波谷的比值等多种参数，根据经验，都能与已知的碳氢化合物蕴藏量有关，而且这种关系也适用于为了新的勘测目的而采集到的震动数据。简单地说，无需钻一眼井，震动数据即可提供地下结构和地层方面一些最翔实的信息。

然而，在以下这个重要场合下，震动数据也有其局限性：一些非常“薄”的岩石部分常常很难探测清楚。特别是，尽管在岩石层相当“厚”时，震动反射数据能够提供描绘地下结构的近似“地质断面图”，但当岩层较“薄”时，由该结果得出的震动图象非常不清晰。这种现象被本领域技术人员称为震动清晰度问题。

震动清晰度在本文中用来指一条震动道中的垂直清晰度，并且放宽设定为是对于两个地下振动反射体之间的最小差别而言的，这两个反射体应能从震动记录上分辨出是两个不同的界面，而不是单个合成的反射面。可以解释为：在震动视图上，地下的物体可以很方便地看作是两部分的合成。一束可辨的反射，来自物体的顶面；第二束可辨的反射，很可能来自该物体的另一端，即底部。在理想情况下，在记录到的震波曲线上，物体的上、下两面都显示出为可区分的独立反射体，各自可分别进行“时间选取”(time picked)，(即进行标记)。在两个标记时刻期间所采集到的数据中，包含着介于其中的岩石部分的信息。另一方面，如果震动物体不是足够厚，则从物体上下两面的反射就会重叠，这样就导致两个反射波之间的干涉，使下面的图形产生混乱。这种图形混乱就是本领域技术人员所称的“薄层”问题的现象之一。

图1用通常方法介绍了在褶积模型原理下，薄层问题是如何产生的。首先考虑图1a中描述的“厚”层反射。图的左侧部分代表从地表产生的源波，源波沿路径P₁向下传播，穿过地层径直到达以“A”标出的岩石部界面(注意，图上的波动路径实际上是垂直的，为叙述方便起见描绘成有一倾斜角度。这与工程上通常的实际情况一致。)在图1a中，下行的震波到达界面“A”后，一部分能量沿路径P₂向回反射，在地表被记录为反射项R₁。注意，R₁波与源波比较，相位相反，因此在界面“A”即指负的反射系数。此处的相位反转仅作为举例说明，本领域技术人员反射系数的正负都是可能的。

下行能量的其余部分(在一部分能量从界面“A”反射之后)继续穿过厚层，到达位于厚石层底部的界面“B”。到达“B”界面后，部分波动能量继续沿路径P₅向地层深处传播，而其余能量沿路径P₄反射回到地面，记录为反射项R₂。注意，从界面“B”处的反射在时间上晚于从界面“A”处的反射。两次反射之间的精确的时间间隔与两个界面之间介质层的厚度和波在其中的穿行速度有关，越是厚层和/或慢速，上下两个反射的时间间隔就越长。岩层的时标厚度即为震波横穿其间所需的时间。

地面上合成的厚层反射，即实际记录到的内容，是在考虑到两次反射间的时间间隔的条件下，两个反射的算术和(叠加)。因两个返回的波在时间上不重迭，故所记录到的震动记录上清楚地将两次反射显示为两个独立地层的图形。(注意，图中并未准确描述出两次反射之间的时间间隔。本领域技术人员知道，准确的时间间隔应是岩层的时标厚度的二倍。)

再看图1b，描述的是薄层反射。同样，地表发出的源波沿路径P₆，到达标记为“C”的岩石部界面(同上，图中的波动路径实际上是垂直的)。如图1b所示，当下行震波到达界面“C”后，一部分能量沿路径P7向回反射至地表，记录作反射R3，其余部分的下行能量继续穿过薄层达到界面“D”。到达D界面后，部分波动能量沿路径P10继续向地层深处传播，而其余能量沿路径P9反射回到地面，记录为反射项R4。

同样要注意到，从界面“D”的反射在时间上晚于从界面C处的反射，但在薄层的情况下，两次反射之间的时间间隔较小，因为波动从界面“D”反射前所经过的距离较短。实际上，两交反射的时间间隔短到使两个返回(上行)波重叠。由于薄层反射的合成同样是两个反射波的算术和，使得实际记录到的信号无法清楚地表达上面反射或下面反射中的任何一个，信号的解释也相应地很难。这个不确定的合成反射是薄层问题的典型例子。

显然，对于石油公司的勘探人员而言，地下勘测目标的厚度是个具有重要经济意义的因素，因为假设其它条件相同，则岩石部分越厚，其蕴藏的碳氢化合物的体积也就越大。由于精确地确定岩层厚度有如此的重要性，人们竭力用种种方法来完善薄层问题也就不足为奇了。

应用最广泛的第一种方法是缩短震波的波长，长波得到的结果往往差于短波。在数据处理阶段，运用熟知的信号处理技术会使记录到的地震波形显著缩短。本领域技术人员通过实例明白，传统预测性的重叠合法可用以白化波谱，减小有效波长。同样，常规的波处理技术，包括源信号的重叠合法及其它多种方式，均可替代使用以期得到一个相同的最终结果：较紧密的波形。尽管以上任一种处理都可使震动断面图的特征产生显著变化，且能明显缩短其波长，但通常还必须采取进一步的措施。

但无论震波有多密集，总会有一些具经济价值的岩层，因为太薄而无法圆满解决上述问题。即使是最好的信号处理专家最后也只能是推迟这种不可避免的问题发生。于是，便用到了其它更宽领域的方法，这些方法目标更集中在对合成反射特征的分析上。这些方法通常是基于观察的结果：即使只有一个合成反射，且岩层厚度不能直接观察到，仍可从记录到的震动数据中找到信息，这些震动数据可以间接用来估计板质部分的实际厚度。

图4a举例介绍了一种常见的“狭缩”震动模型，其中所研究的地层(此处其厚度用穿行时间测定，而不是用长度)厚度逐渐减小(即“狭缩”)直至在最左端消失。图4b是用这种模型计算合成的地震道族，描述了无噪声条件下震波的褶积及岩层表面的界面。注意到图4b最右侧边缘，第一条震道上记录的合成信号使反射体可以明显地分为两部分，上部发出负反射，底部发出正反射。在图4b中左移，上下部分各自的反射开始合成并成一个单独的合成反射，最后当中间的厚度减小到零时反射消失。但需注意，即使在退化为单个反射之后，此合成反射的性质仍在继续变化。因此，即使几乎没有直观特征表明该反射来自两个界面，但在岩层厚度减少时，合成反射所呈现出的变化仍说明，在适当条件下，这些反射中所包含的信息可提供与薄层厚度有关的一些内容。

Widess在1973年进行的先驱工作(见Widess的《薄层厚度》地球物理，卷38，第1176-1180页)已经创立出一种薄层分析的简易方法，其中依赖于薄层合成反射波的波峰-波谷振幅的校准曲线得到完善，连同波峰-波谷间的时间间隔，就可以估算薄层的近似厚度。(参见Neidell和Poggiagliomi的《地层模型及说明一地球物理学原理与技术》，“震动地层学在碳氢化合物勘探中的应用”，AAPG Memoir 26.1977)。在薄层问题的情况下，建立起“调谐”振幅是进行校准曲线处理的必要步骤，即在岩层厚度达到使上下两部分的反射发生最大正干涉时，便出现了调谐振幅。至少在理论上，谐振厚度仅依赖于主波的波长λ，当上下两面的反射系数同号时为λ/2，异号时为λ/4。

由于校准曲线型方法的适用性，现已被成功用于多种勘测设备。然而这些由振幅和时间决定的校准曲线方法，要在很大程度上依赖于准确的震动处理建立正确的波动相位，并要控制相应的逐线震波幅度。震动处理方面的技术人员知道，要得到一个能真正保持各处都有相应振幅的震动断面图，其难度很大。进一步说，上述揭示的由校准曲线确定的方法，在检测大范围三维勘查中薄层反射信号时，并不能很好地适用：该法最适合用于单个震波线上独立的反射体的场合。很难由一条震波线得出校准曲线，而要找到一条在整个三维结构震动数据中都有效的校准曲线要更难了。

迄今为止，在震动数据处理和说明领域，需要一种方法从以传统方法得到的震动数据中提取有用的薄层信息已是众所周知的了，但这些震动数据正面临着上述问题。进而，这种方法最好还能为后续的震动地层分析、结构分析提供有用的内容。因此，如本发明人所见，现在人们也需认识到，要找到一种能提出并解决上述难题的震动数据处理方法，目前存在着对这种方法的迫切需求，并且已经存在一段时间了。

但在开始叙述本发明之前，需注意和记住，下面的叙述连同附图，不应被视为只将该发明局限于所展示、描述的例子(或优选实例)。这是因为本领域的技术人员也有能力在本文所附权利要求范围内，提出同样形式的发明。最后，尽管可以根据褶积模型的不同方面对本发明作出解释，但下述方法并不依赖于记录到的地震道的某个特殊模型，而且在偏差较大情况下与标准的褶积模型同样有效。

发明内容

本发明人找到一种新方法，利用离散的傅里叶变换，通过传统的二维和三维地震数据，模拟描绘薄层范围和其它岩层的横向间断范围。特别是，本发明的动机来源于观察到，在指示地层厚度的领域内，来自薄层的反射具有一种特征图象：匀质的薄层会使合成反射的振幅谱线中产生一列周期性的凹槽，相邻凹槽之间的距离反比于薄层的时标厚度。进而，如果能正确列出傅里叶变换系数，解读器就可以利用这个特征图象，在三维空间内追踪薄层反射，估测薄层的厚度，并能达到迄今还不可能的精度。更一般地，这里所介绍的方法还可用于检查、确定所测位置岩石群的垂直及水平间断性。同时，为强调谱线中所提供的地质信息而出现的一种新的频域白化方法，也进一步提高了本发明的实用性。最后，本发明亦可直接用于揭示与所关心的地下结构、地层方面特征相关的一些震动特性，即向勘测人员提供数值，使他们可以制图并用以判断地下碳氢化合物及其它矿物的储量。

本发明的技术方案是这样的：一种通过定位周期性频域凹槽来显示薄层的方法，其在处理包含分布在预定的地层体内的三维地震道集的地震数据中实现，该预定的地层体包括有助于产生、迁移、累积或呈现碳氢化合物的薄层，该三维地震道集由至少两个不同的地震线组成，每个地震道与所述地震线中的至少一个对应，该地震道包括至少由时间、位置及振幅表征的数字样本，该方法包括以下步骤：(a)选择一部分所述地层体和其中包含的地震道和数字样本，以此在所述地层体内确定一个三维抽样层，所述选择的地震道在至少两个不同的地震线上与至少两个不同的预定的地震线对应；(b)在所述三维抽样层内利用离散的正交变换对所选择的地震道的至少一部分和所述选择的数字样本进行变换，并从如此变换的地震道产生变换系数；(c)将所述变换系数组织成为调谐立方体，用于地震勘探，所述调谐立方体包括所述变换的三维阵列；(d)使用计算机制作和交互显示所述调谐立方体的变换系数的至少一部分，并标识一个或多个地下地质特征，该地下地质特征有助于在所述预定的地层体内产生、迁移、累积或呈现碳氢化合物。

基于大致背景技术，本发明选择采用较简短的离散傅里叶变换，来确定地震道的频率组成。如本领域技术人员所熟知，对时间级数的傅里叶变换计算中，即使全部为实数值，结果也将出现复数傅里叶变换系数组，这些复数形如“A+Bi”，其中“i”为“虚”数或负一的平方根。进而，已知表达式“A+Bi”可等价写为：

A+Bi＝re^iθ

其中，

r = | A + Bi | = \sqrt{A^{2} + B^{2}},

及

θ = ta n^{- 1} (\frac{A}{B})

量θ称为复数A+Bi的相角(或仅称为“相”)，量“r”为模，表达式|A+Bi|为复数模的计算式；也称为绝对值。对每个变换系数计算其模，即可得到频谱。频谱中每个系数的数值大小正比于该频率在原始数据中的强度。最后，在对一些特殊的时间级数进行完傅里叶变换，最终的复系数级数就可以落入频率域，而未经变换的数据则被认为落入时间域。

再回到本发明的讨论，本发明依赖于大致观察得到的结果，即整条震道所作的傅里叶变换得到的频谱，跟源波的频谱很相似，越短的频谱窗就能反映出更多下面的地质信息。这是由于长分析窗中包含了大量地质方面的变量，这些变量在长期内会呈现出一种“白”(或者说随机的，无关的)反射函数，具有“平坦的”振幅波谱。这样，从整条地震道上计算得到的频谱，就要在很大程度上依赖源波的频率内容(实例参见地球物理勘探学会Ozdogan Yilmaz的《震动数据处理》(1987)第2.2.1节，其阐述的内容引入本文，此供参考)。另一方面，在分析窗足够短使地层反射函数非白时，所得傅里叶频谱既受波动的影响，又受当地地质条件的影响。在这样小的窗口下，地质因素就充当了滤波器，衰减源波的谱线，增加非稳定的短窗谱线。

图2大致描述了上述构思。其中绘出了一个典型的地震道和由其运算得出的几个频谱，图的上部是整条地震道的傅里叶变换频谱。谱线的形状基本上就是典型的激波谱线。而下部绘制的由计算得出的短窗谱线，则是非稳定的，能反映出可在极短的时间内发生显著变化的地下地质特征。

图3介绍了本发明这种观察方法的重要性，图中大致绘出了两种有代表性的谱线。上部的频谱图代表典型的宽频带源波，而下部的频谱图则以更常见的方式表达了合成薄层反射的频域图象。在后一种情况下，薄层的地质状况被用来作为频域滤波器，并在反射波的频率内容上作出了标记。如图中所大致介绍的，本发明人发现，匀质薄层在加入凹槽或窄消频带后，对反射波的振幅波谱产生了影响，由此显出特殊的形状。在匀质层中，速度和密度处处保持恒定。这里所述凹槽之间的距离等于薄层时标厚度的倒数。时标厚度是指波沿一个方向穿过岩层所需的时间长短(岩层的厚度除以速度)。这样，就可以利用振幅波谱中被衰减的频率来确定某个薄层的反射，并测量其厚度。

在图4中，将前一段中得到的结果推广，用于分析简化的二维地质模型，该模型中薄层的频域表达已经得出。图4a表达了典型的“狭缩”反射函数(地质模型)，图4c用灰度标记表达了由这种模型计算出的傅里叶变换频谱振幅。这个标记由对模型进行五十等分得到时间级数而产生。每部分中只有两个非零值：一个对应于上层的反射系数，另一个对应于下层的反射系数。非零值位于其它各处均为零的震道上，每个非零值所在位置分别用时间值反映了上下两部分反射体的位置。然后即可对时间级数进行标准的离散型傅里叶变换，再对每个系数求出复数振幅。

图4c上较浅的部分对应于振幅谱线的较大值，较深的部分对应于较小值。因此，振幅谱线中的“凹槽”在图上就表示为较深色的数值。从最严格的意义上讲，此图表示了地质上的傅里叶变换；特别，在本例中还表示了在波上所带的特征信号。此图中涉及到本发明的最重要之处在于，当模型厚度减小时，凹槽之间的距离增加。进而，对于厚度一定的模型，凹槽为周期性的，其周期等于该处时标厚度的倒数。故，如果能在震动勘测图上确定这个信号，即周期性频域凹槽的位置，这就是薄层存在的有力证据。

按照本发明的一个目的，已经提出了一种解译包括薄层情况在内的震动数据的方法，其中数据被分解为傅里叶变换的二维或三维数列，以此提供了一种改进的描绘所述薄层范围的图形的方法。本实例利用一个单独的傅里叶变换窗口，并对每个地震道横断所测区域的部分使用该窗口。图5中利用三维地震数据对本实例作了大致介绍，但本领域技术人员会想到，如果将此方法应用于二维地震道上，会有利于提高其中所含薄层反射的可视性。

第一步先将空间上相关的地震道集中在一起，仅为介绍起见，这些震道也许只是一条或几条发射记录，一个恒定的偏移选排，一个CMP选排，一幅VSP勘测图，一条二维地震线，一条从三维地震勘测图中提取出的二维叠加地震线，或更合理地，是一个三维勘测图中的三维块。此外，本发明还可用于数据变换后的二维线或三维体地震勘测图，即图中的“偏移量”或空间轴(三维坐标中的“X”轴或“Y”轴)，已定向，代替纵轴或“时间”轴。更一般地，任何三维数字量均可由这里揭示的方法进行处理。为了明了，下文将所述的纵轴称为“时间”轴，当然，本领域技术人员会懂得，数字样本不应以时间单位来划分。无论怎样选取，对于在地下地质特征方面隐含着下层空间相关性的地震道组，应用这里提出的发明都是最有效的。同样只为说明起见，虽然可想而知可以利用任何同类的空间相关地震道组，但下面的讨论将只依据叠加三维勘测图中包含的震道进行表述。

如图5所示，下一步从特定的三维体范围内选取一个抽样层，作为实例的该抽样层可以象图5中那样，是由两个选定反射体作为界面的波状区域。这样，反射体在分析之前即可按需要压平或数据化(即按时间将单震道上下移动达到平整)，并有可能复原重组。习惯上，可以确定一个特定的时间间隔(比如从2200毫秒至2400毫秒)，这样形成一个“立方体”，或更准确地说，在三维体内的一个地震数据“箱”；一个子体。此外，可以靠确定亦线“线上”或“相交线”的界限来限定抽样层的水平延伸范围。其他确定抽样层的方法当然也是可以的，这些方法本发明人均已考虑过。

这步挑选，提取与抽样层有关的数据称为对数据“建立子集”(见图5)。挑选抽样层的一个指导原则是希望是保持该层(在时间上)尽量短。由于考虑到将要用长窗傅里叶变换谱来集中子波，用短窗傅里叶变换谱来容纳地质方面的信息，这样的选层方法与前面所采取的总的基本原则是相符的。还有一个“隐藏的”扩展窗口。通常可以自动、随意地用于傅里叶变换窗口：使窗口长度延长到2的幂。这样对窗口的延长是为了计算方便，窗口的长度为2的幂，便于进行简称FFT的快速傅里叶变换算法。但是，本发明人明确建议不要使用这种在工程上普遍使用的方法，而倾向于使用更常规的离散型傅里叶变换算法，虽然在计算上效率不高，但可以使分析窗的长度保持尽可能最小的数值。对于当今计算机强大的运算能力，没有理由不在抽样层内只进行数据的变换。

图5中，用于本实例“运算”步骤包括至少一步运算过程：对抽样层进行离散型傅里叶变换的计算。所得系数，抽样层的谱分解，都被作为输出的谱分解块(“调谐立方体”)的一部分存储起来，准备后面的观测。注意到每条地震道的输出调谐立方体中，将有一震道(即傅里叶变换系数集)经处理后作为输入量的一部分。还要注意到，这里优先采用的输出排列中，在数据体的水平面薄片内部包括有与单个普通傅里叶频率相应的系数。

还可以按需要，使“运算”步骤中包括其它运算过程，以便改进输出量及后继分析的质量。首先，在作变换计算之前，可以将权函数加于抽样层内地震数据上。权函数的作用是在傅里叶分析的窗口内修剪数据，减小由“棚车”状分析窗所导致的频域畸变。本领域技术人员熟悉在变换前使用加权函数。这里揭示的发明所优选的加权函数取高斯函数的形式，该函数应用于此种场合从许多方面说都是最好的。也就是说，其它多种加权函数也都是可行的。

另外，由于勘探人员往往最关心振幅谱线，在将转换系数存入备用存储区后，振幅谱线可由转换系数计算得出。相谱或其它一些导出性质也都同样可由转换系数计算得出，然后保存。实际上本发明人也已用到了这些计算。

最后，作为运算步骤的一部分，对输出值体的每个平面(即频率)都可单独进行频率定标。如图10所示，本发明人发现，如果对输出值体中每个频率薄片都分别定标，使其具有相同的均值，然后进行检测，这样做会有许多好处。这种定标法仅是许多种可用方法之一，但发明人倾向于此法，因其虽损失了一般的子波信息，却改进了所存频谱中的地质内容。

谱线经计算、保存后，即可备用于薄层的地球物理勘探。注意，当数据随后罗列出后，各条谱线应按与其它谱线有相同的空间相关性进行组织和检测，正如由之计算导出它们的震道那样，这一点非常重要。即，变换前数据中呈现的空间相关性，在处理转换系数时必须保持。这里倾向使用的检测转换系数的方法是，首先，当然在假设输入数据是从原始的三维体中得来的条件下，将这些数据转换形成三维“体”(调谐立方体)。但要注意，竖轴(“Z”)不再象变换前那样表示“时间”，而是按规定，表示频率单位，如存于其中的傅里叶变换系数那样。

图5中上一步所介绍的调谐立方体，现在可用检测传统的三维震动数据体的任何方式来检测。前面谈到，本发明人已发现，在系数体检测连续的水平薄片，是确定和显化薄层效应的一种较好方法。注意，在调谐立方体处理中，水平薄片代表着与单个傅里叶频率对应的全部系数，因此在整个面内频率保持恒定。进而，发明人倾向于沿数据体作一系列动态的水平检测，协助对体内所含数据进行分析。在抽样层是单条地震测线的局部而不是测体的局部的情况下，最终的显示内容即使并不是数据“立方体”，这里也仍被称为调谐立方体，呈现为傅里叶变换谱按其初始空间相关性所显示的一族空间相关的地震道。

使连续的水平薄片在谱线体中移动，是检测、分析变换系数的一个较好方法，这种动态化最好能在高速工作站的计算机显示器上实现。本领域技术人员熟知，这种以在数据体中交互式追踪形式的动态化，是检测大量数据时既快又省力的方法。数据体的检测可在水平方向、垂直方向或倾斜方向的薄片上进行，每种都能得到相同的检测。而更重要的是，在本发明所述范围内，对连续的水平薄片依次进行快速检测，提供了一种分析大量数据并从中确定薄层反射的有特征的方法，其细节容后讨论。注意，对于这里揭示的方法，在对薄片进行动态化和检测时，最好将这些薄片按频率顺序(严格地升序或降序)进行排列。

按照本发明的第二个目的，提供了一种地震数据处理系统，可以改善薄层情况下震动记录的显示状况，其中用一系列重叠的短窗口傅里叶变换，将数据分解成为一系列傅里叶变换二维线或三维体，由此提供薄层部分的改进描图法。图6介绍了本实例应用于三维地震数据的大致情况。但本领域技术人员会意识到，同样的方法也可以用于二维地震线，改进该处所含薄层反射的可视性。如图6和上文揭示的，本实例的第一步包括解译器描绘抽样震层的即时边界。如前所述，这样的描图会得到单条地震线的震动数据立方体或矩形片。

本实例中，对每条震道并不是应用一个单窗口傅里叶变换，而是用到了一系列重叠的短窗口傅里叶变换。根据应用场合不同，窗口长度和重叠的个数也各不相同，但同样，窗口长度不必等于2的幂，而应进行选择，以便能对下面的地质状况作最佳描绘。注意，在进行变换之前，权重应可选择性地加到每个短窗内的数据上，并且同样，高斯权重是较佳的选择。

如图6指出的，经过每个短窗傅里叶变换计算，由此得到的系数分别存放于单个的调谐立方体中，这些调谐立方体应分别保持相伴于产生它们的短窗。注意，在本例中，有多少重叠的分析窗口，就会有多少调谐立方体产生。如果需要，则应对每个调谐立方体中的每个频率面分别进行定标。

每个由活动窗口产生的短窗调谐立方体，应象对前面第一个实例所建议的，可以被单独精确检验。同样，应通过水平薄片或恒频图，能对每个立方体进行检测，这样提供了一种利用频率使地质变化显化的方法。进而，由于有了在震道上不同时间点处算出的一个调谐立方体族，实际上也就产生了一个增大地下深度的调谐立方体族。

最后，按照本发明的第三个目的，提供了一种地震数据处理系统，可以改善薄层情况下震动记录的显示状况，其中用一个短窗傅里叶变换，将数据分解成为一系列傅里叶变换二维线或三维体，然后重组为单频率调谐立方体，由此提供了对薄层部分的改进描述。

如图7所大致介绍的，本实例中的第一步恰与前两个实例中的相同：首先解译数据，然后建立子集。之后，从抽样层中的地震数据来计算一系列重叠的短窗傅里叶变换；在计算变换之前，还可任选地对每个窗口内应用加权或减弱。在前面实例中，每个短窗变换得到的系数都进行累加。然而在本例中，不是将所得的傅里叶变换系数当作调谐立方体，而是数据重组为单频能量立方体，这样可以在水平面或垂直面内查找薄层影响的证据。

更特殊地，本发明人所考虑的重组，从概念上包含着从所有调谐立方体中提取各个对应于特定频率的水平薄片。其后，将这些单个的同频薄片“叠放”在一起，最上边的一片含有从最上边活动窗口计算出的系数，第二片含有从顶部下方第一个活动窗口计算得出的系数，等等。注意，重组之后，系数体组织成为“X-Y”和时间的单元。这是由于纵轴表示产生特定系数的活动窗口的“时间”顺序。

为了利用按上述步骤所构成单频调谐立方体所带的信息，地震解译员会选择一种频率及与其相应的震动体(如可能选择对应于10Hz和/或11Hz的系数值，等等)。每个“恒频立方体可从平面视图或水平视图，或其它任意方式进行检测，这样就提供了一种方法，可对特定频率显化水平延伸的地质变化。

在所有上述实例中，由于变换前的初始震道是空间上相关的，这个事实便为本发明提供了新的用途，认识到这一点很重要。特别是，众所周知，短窗傅里叶变换系数本身带有大量杂波，与长窗变换相比频率清晰度很差。为提高变换后数值的可靠性，本发明人所采用的一种方法是对转换前的数据值加上高斯权重。不过，本发明人所采取的另一同样重要的步骤是，在数据体内显示出象输入数据那样具有同样空间相关性的系数。由于这样显示的震道包含空间相关的信息，将它们邻接显示可以使在观察时消除杂波，得到潜在的清晰信号信息。

最后，尽管本发明中是根据离散的傅里叶变换形式进行讨论的，在实际中，有多种恰好也可以这样用的离散时间数据变换，傅里叶变换只是其中之一。一般的步骤为(1)计算短窗变换(2)将所得变换系数统一为数据体(3)从该体中查找薄层影响的迹象；这些步骤均用除傅里叶变换之外的其它多种离散数据变换来完成。如果所作变换并非傅氏变换，则应按原函数相同对系数一起进行分组，形成调谐立方体。这样，当下文用到“单频调谐立方体”时，本发明人的意指不仅仅是由传统傅里叶变换系数形成的调谐立方体，而是由具有相同原函数的系数组成的任一立方体。

本领域技术人员会明白，离散傅里叶变换只是多种离散的线性保积变换之一，这些变换满足以下性质：(1)它们为线性算子(2)可逆，(3)它们的原函数形成标准正交集。用公式表示，如果x(k)，k＝1，L，代表时间级数，X(n)是其n次变换值，n＝1，L，则时间级数的正变换对于这种同类的变换，一般可写为：

X (n) = Σ_{k = 0}^{L - 1} x (k) A (k; n)

其中A(k；n)为正向变换核或原函数族。进而，存在逆变换，可将变换后所得值变回到原数值：

X (k) = Σ_{n = 0}^{L - 1} X (n) B (k; n)

其中B(k；n)为逆变换核。正交化条件是指两个不同的原函数其内积必须为零，每个原函数的模必须为单位一。此项条件可简要概括为下式

Σ_{n = 0}^{L - 1} A (j; n) A * (k; n) = δ (j - k)

Σ_{k = 0}^{L - 1} A (k; n) A * (k; m) = δ (n - m)

其中，

A^*(k；n)表示A(k；n)的共轭复数。对于离散傅里叶变换，与对长度L进行正变换相对应的原函数，通常选为复指数函数集：

A (k; n) = {e^{- 2 πikn / L}, k = 0, L - 1}

这样就有L原函数(或此例中的基向量)，对每个“n”值均有一个原函数，其中，

n = - \frac{1}{2}, \cdot \cdot \cdot, 0, \cdot \cdot \cdot, \frac{L}{2} - 1

概括起来：每个变换系数X(n)由与特定原函数相对应的数据窗口算得；集中与特定抽样层对应的所有变换系数，构成调谐立方体；将这些系数存入备用存储区，这些系数的空间排列应与经计算得出窗口的那些震道相同。

借助另一个具体例子，本领域技术人员明白，离散傅里叶变换也可代之以离散的沃尔什变换，沃尔什系数也类似进行分组、显示、分析。按上面揭示的方式，可以在一系列重叠的活动窗口内进行沃尔什变换，由此而得的系数经过组织，存储至调谐立方体。当然，本领域技术人员用这些系数代表类似的量“序列”，而不是用这些所得变换系数代表频率。这样，“单序列”调谐立方体可由沃尔什变换系数形成，其方式与构成傅里叶调谐立方体的方式完全类似。下文中将这些“同频(或更一般地说，同原函数)调谐立方体”称为单正交原函数调谐立方体。

最后，尽管离散的傅里叶变换具有正交原函数集的特征，但在变换前对原函数所施加的非无效权函数则破坏了其正交性。按传统理论，对窗口所加权函数被视为加于原函数，而不是加于数据，以保护下面数据的整体性。但是，正交的原函数在加权后通常不能再正交。亦即，实际上无论是对数据还是对原函数加权函数，变换所得最终结果恰好相同。

为避免对离散的正交交换加权时产生的较小的理论困难，要选择一个不受此影响的正交变换/加权的二者结合。按照例子，对局部进行的余弦(或局部的正弦)变换是一种离散的正交变换，其中选用的权函数呈平滑、精细递减形式，这样就在频率清晰度有所下降的情况下，保持了原函数的正交性，进而，本身余弦/正弦变换蕴含的基本原理也为通向普遍的子波变换领域架起了理论之桥。

附图说明

图1为概述薄层问题原理的示意图；

图2示出了典型的地震道，并对比了由其计算得出的长、短窗谱线；

图3概述了如何在频域上表达对薄层发射地震子波的反应情况；

图4包括简化的震动狭缩模型，其产生的褶积反应，及这种褶积反应的频域表示；

图5为介绍本实例大致方法的示意图；

图6为介绍本发明所选实例如何用于勘探的示意图；

图7为介绍另一实例的示意图；

图8为介绍本实例的流程图；

图9为介绍移动恒频薄片过程中，薄层形状的示意图；

图10介绍为恒频薄片定标以改进变换后数据的地质内容的一般方法；

图11为介绍另一实例的示意图。

具体实施方式

本发明提供了一种利用离散傅里叶变换处理地震数据的方法，由此改进了其作为薄层探测器的用途。

按照本发明的第一种方案，提供了一种利用离散的傅里叶变换改进和检测地层效应的方法。其中对抽样层确定的窗口进行了一次傅里叶变换，由此得到的系数用新的方式显示出。如图5中大致介绍的，令x(k，j，nt)代表一个三维地震数据体，其中k＝1，K，而j＝1，J代表在给定三维体内确定特定震道的指针。仅举例，这些指针也许是沿着线上和与线相交的位置数字，尽管其它定位法也是可能的。变量“nt”用来表示每条地震道内的时间(或深度)位置，nt＝0，NTOT-1，在单条震道上所有样本点的数字。连续值x(k，j，nt)(即抽样率)之间的时间间隔用Δt表示，其中Δt按习惯以毫秒为单位。这样，三维体内每条震道均包括(NTOT)^*Δt毫秒个数据记录，第一个样本按习惯取自“零”时刻。需要声明，本领域技术人员知道，根据这里揭示的发明，一些尤其适于分析的地震数据并不是按“时间”排序的。仅举例说明，利用深度迁移程序处理过的震动数据样本是按照深度增量Δz的顺序存储到地震道中的。然而，本发明能够而且已经被完全按此种方式应用于这类数据中。这样，在下文中，使用Δt(及“时间”)具有更广泛的意义，表示连续的数字样本的间隔，而不必关心这种间隔采用的是什么计量方式。

开始时，勘探人员和震动解译员要从三维体中选出抽样层。仅举例说明，这一步可通过数字化时间选取(“拾取”)震动现象来完成，选取范围可以是数字化图表，或更普遍地，从震动工作站中选取。在选取到一例后，勘探人员就试图精确定出从每条地震道上所呈现出的同一反射体的特征(如波峰，波谷，零交叉等)其最终目标是产生包含时间和表面位置信息的计算机文件，这些信息在二维面和三维体内追踪着地震活动。如图11所示，对于给定的这类信息，计算机程序可设计为，在数据体内为任意一条震道读出选取到的数据并找出抽样层；和/或演示本发明的方法。上述程序可由磁盘、磁带、光盘或CD-ROM装入计算机中。

作为替换手段，解译人员还可以在整个体内设定所研究震动的起止时间范围，形成一个抽样“立方体”，该“立方体”在通常意义上代表初始三维勘测体的三维子体。为介绍起见，尽管本领域技术人员会意识到，后面讨论的同样技术可简单地应用于非恒时窗口，在下面的讨论中还是要假设已提取出了三维子立方体。同样为了介绍所公开的技术，需假设提取后的时标抽样层范围，要从三维子体的第一个样本到最后一个，即下文所述的样本号“N”。同样，尽管本领域技术人员会意识到在通常情况下，抽样层的延伸范围只是三维体的一部；但在后文中还是要假设抽样层出现于子体中的每条震道上。

抽样层确定后，第二步是选择傅里叶变换窗口的长度，下文均以“L”表示。一般而言，变换窗口的长度应仅为能够反映抽样层全部所绝对必要的。习惯上按计算方便来选择傅里叶变换的长度，通常限定为2的整数幂(如32，64，128等)，这样便于使用快速简便的FFT算法，而不选用稍繁一些的混合基(mixed radix)傅里叶变换或更繁的常规离散傅里叶变换。但是在本发明中，发明人特意指出，不要将所选窗口的长度按惯例扩展为2的整数幂，而是采用更一般的离散傅里叶变换。需声明，在下面的讨论中，本领域技术人员会知道，当需要用到离散的傅里叶变换时，如条件许可便用FFT算法。否则如果所选窗口长度不是2的整数幂，就使用一般的离散傅里叶变换或混合基变量(mixedradix variant)。

在开始傅里叶变换前，应建立一个备用存储空间，用以存放算得的傅里叶系数。为了将所得变换系数存放在震道中，必须要为每条震道提供的最小辅助存储空间，其大小应能显示L个计算机字。如果震动数据值或变换结果必须按双倍(或更高)精度保存，则还需要更大的存放区。需要解释的是，对于长度L的实数时间级数进行傅里叶变换，需要L/2个复数值的存储空间，其中每一个通常应有两个计算机字的空间。实际上只有[(L/2)-1]个单值复数值，而不是L个，因为对于实数时间级数与正负频率对应的傅里叶变换系数是直接相关的：它们是一对共轭复数。此外，有两个实数值：可一起存放于一个复数值长度内，它们是：零(“dc”)赫兹系数和尼奎斯特(Nyquist)频率处的系数。最后，如果L是奇数，则单数据值的数目为[(L+1)/2]个。如果在抽样层(子立方体)中总计有J×K条地震道，则需要的备用存储区的总数，至少应该是按计算机字得出的L，J，与K的乘积。本实例中令级数A(k，j，nt)表示备用存储区容量。

如图8所示，运算的第一步是从以下子体中取得的输入震道x(j，k，nt)上提取抽样层数据值；该子体为：

y(nl)＝x(j，k，nl)，nl＝0，L-1

备选权函数为：

y(nl)＝y(nl)w(nl)，nl＝0，L-1

其中，数值y(nl)为暂存区。(注意，本实例中分析窗长度等于抽样层长度。)权函数w(t)——或也被称为数据窗——可以取许多种形式。一些较常见的数据窗口类型有Hamming，Hanning，Parzen，Bartlett和Blackman型。每种窗函数各有其优缺点。但本发明人出现，在本项应用中高斯窗口从许多角度看都是最好的。高斯权函数由下式定义：

w (nl) = σ_{3} e^{- {(nl - 1 / 2)}^{2} / σ_{2}}, nl = 0, L - 1

其中，

σ_{1} = \frac{L}{6},

σ_{2} = 2 σ_{1}^{2},

σ_{3} = \frac{1}{\sqrt{{2 πσ}_{1}}}

但通常权函数应该为实函数，且在零点无定义。

应用权函数值后，离散傅里叶变换即可按下列标准式进行计算：

X (n) = Σ_{k = 0}^{L - 1} y (k) e^{- 2 πikn / L},

n = - \frac{L}{2}, \cdot \cdot \cdot, 0, \cdot \cdot \cdot \frac{L}{2} - 1

其中X(n)表示傅里叶变换在频率f_n时的复系数，即该频率依赖于窗口长度L。通常众所周知，傅里叶变换产生的复系数，在以下傅里叶频率时提供了对谱线振幅的估值。

f (n) = \frac{n}{L (Δt / 1000)},

n = - \frac{L}{2}, \cdot \cdot \cdot, 0, \cdot \cdot \cdot \frac{L}{2} - 1

在此应注意到，地震道上的额定取样间隔Δt不应等于在勘测范围内获取数据的取样间隔。例如，实际中通常在对地震道重新取样时，选择较大的取样间隔，这样当在最高的记录频率上几乎没有有用的信息时，可以保护存储内容。另一方面，有些情况下，例如当地震道需与其它(更高的取样率)直线结合时，则可按较低的取样率重新取样。在两种情况下，数据的额定取样率可能都不能准确反映实际谱线带宽。对前式的一个简单调整就容纳了这种偶然性：

f_{n} = \frac{n}{L} F_{\max},

n = - \frac{L}{2}, \cdot \cdot \cdot, 0, \cdot \cdot \cdot \frac{L}{2} - 1

其中F_max是数据中所含的最高频率。

由于地震道是个“实的”函数(即，不是虚的)，其傅里叶变换是对称的，与正负频率对应的傅里叶系数有如下关系：

RE[X(f_n)]＝RE[X(f_-n)]

及

IM[X(f_n)]＝-IM[X(f_-n)]

其中RE[z]是提取复数z实部的函数，IM[z]是提取虚部的函数。由这个关系可知，在每个傅里叶变换窗内只得到了L/2-1个单值。因此，为了便于研究特性，下面讨论中将只考虑正频率，尽管本领域技术人员明白，只利用负频也可得到同样的结果。

下一步处理过程包括将算得的复数频率值置入备用存储级数。这些震道上充满了如下形式算出的傅里叶变换复系数：

A(j，k，L)＝X(i)，i＝0，L/2

其中，“j”和“k”相当于与原始数据震道对应的指针。实际中，级数A(j，k，i)可能从未完全在RAM(随机存储器)中保存过，但可能全部或部分地位于磁带、磁盘、光盘或其它存储设备上。另外，由于这里选用的薄层显示需要用到频率谱线，而不是复数值，故在将每个系数置入备用存储级数：

A(j，k，i)＝|X(i)|，i＝0，L/2的同时计算复数模会较为简便。然而，在许多场合下复系数是必要且有用的。于是如图8所示，复系数被优先存入了备用存储区。

上述方案要对设定的子体内的每条震道重复使用，这样备用存储级数中就装满了变换系数，供勘探人员检测。然而，在观察结果之前，要用新的方式对数据优先定标，这样，变换系数中的地质信息就变得与子波的作用相关了。图10介绍了这种包含于频域定标中的通用方法。这里揭示的定标法，能够使每个频率薄片中谱线振幅的平均值相等，由此便得到白化的子波谱线。如图8所详细介绍的，设T(j，k，i)表示将把整个调谐立方体存入其中的暂时存储级数。对于已知的频率薄片i，计算其中谱线振幅的平均值：

YAVG = \frac{1}{JK} Σ_{J = 1}^{J} Σ_{K = 1}^{K} | T (j, k, i) |

因T(j，k，i)可能取到复数值，所以谱线振幅已算出。第二步，在这个特定的频率薄片中对数字进行调整，使它们的平均值等于用户预先设定的常数，表示为变量AVG：

T (j, k, i) = \frac{AVG}{TAVG} T (j, k, i), j = 1, J, k = 1, K,

其中所设定符号用来指出，级数T(j，k，i)已经过调整。在实际中，AVG被设为特定的数值，例如100。在调谐立方体中，要对每个频率薄片(i＝0，L/2)分别重复使用这种标定法。这一过程的结果是，每个薄片具有相同的平均振幅，并完成了一种对谱线的配平。注意，这种同频率定标法只是可应用调谐立方体数据的一种定标算法，本发明人已考虑到，用其它方法也许会更好一些。举例说，不计算薄片中某项的算术平均值，而用另外的中央趋向特性或其它任意统计量(如中值，最频值，几何均值，偏离值，等等)也都是等效的。再举一例，不是将各个频率薄片中的平均值设为等于同一常数，而是设每个薄片可以等于不同的平均常数，这样就在谱线中突出了一些频率，消除了另外一些。

如果用标准的傅里叶逆变换将定标后的调谐立方体数据返回到时间域上，就得到原始输入地震道的谱平衡型式。令X(n)代表由前述处理方法得到的变换系数定标组，该组系由定标调谐立方体内位置(j，k)处取得。则输入数据的谱白化型式可由下式表示，

x^{'} (j, k, n 1) = \frac{1}{L} * \frac{1}{w (n 1)} \cdot Σ_{k = 0}^{L - 1} X (k) e^{+ 2 πik (n 1) / L}, n 1 = 0, L - 1

其中x′(j，k，nl)代表输入数据x(j，k，nl)调整后(谱平衡)的型式。除数w(nl)在这里消除了变换前所加权函数的影响。如果正变换前未进行加权，则该项可省略。

然而，这里选取的用法，其好处并不在于对定标调谐立方体的逆变换，而在于作为勘测薄层的工具。在所有震道经过处理并放入备用存储区后，水平的(恒频)对应第i个频率的振幅片Si(j，k)即可从A(j，k，I)中提出，用以检测和/或动态化：

Si(j，k)＝|A(j，k，i)|

通过将这些薄片动态化(即连续快速地检测)薄层即可象那些在高低振幅值间连续变动的参数一样，进行重新组建。进而，对于许多类薄层，有一种特别的移动凹槽的模式，可以清晰地预示由薄层产生了参数。注意，在此处揭示的方法中，当将薄片进行动态化检测时，最好按频率对薄片排序(严格升或严格降)。

图9介绍了这种缜密移动型式的基本状态，图9a包括了透镜型地质薄层模型，图9b对所述模型进行程式化的傅里叶变换，其中只画出了凹槽。如前所述，凹槽为周期性的，其周期等于所在处模型暂时厚度的倒数。现在考虑图9a中的模型，代表了三维(圆盘状)径向对称模型的二维断面，图9b是同样径向对称的所述三维模型傅里叶变换的一维曲线组。如果按图示，将标号平面1的恒频在数据体内移动，则平面视图上将显示出与第一个凹槽对应的低振幅圆区。平面2穿过两个凹槽，呈现出两个低振幅圆区。最后，平面3包括三个低振幅圆区，对应于它所截的三个凹槽。现在如果按频率递增顺序快速连续观察这些薄片，就会有一种生出“牛眼”状图案的视觉印象，其中环从中心向外运动。这种凹槽运动图案对于薄层是细密的。

如果薄层不是圆型，可以考察一种相关图案。这样呈现出的将不是同心圆，而是一系列从厚区向薄区连续移动的凹槽。例如，将图9当作透镜型气道的截面图。从连续的频率薄片视图上看，沿其长度方向会看到一种向外运动的凹槽——即从管道中心朝向边缘的图案。

要注意到，如果薄层是非匀质的，例如，如含有速度的递增或递减，就将不再出现匀质薄层那种特征鲜明的“凹槽”图案，而显示出有一些不同的频率域。在这些情况下，确定反应特征的较好方式是建立该活动的模型，进行傅里叶变换，如前面图4所示的。利用这些信息，勘探人员可对预定的反应考察动态化的调谐立方体。

凹槽图案不仅可以高质量地指出均质薄层，而且它还可以对薄层范围进行量化测算。回到图9中，注意到凹槽由模型的最外层边缘在水平方向上限定。这样，追踪观察叠放在一起的频率薄片，并注意到凹槽运动的最外界限，就可以得到地层范围的量化估计。

以上就是在实际地震数据体内可观察到的显著视觉效果。由于典型的非薄层情况内有一些始终存在的缓慢变化的振幅谱线，薄层反映很明显，可容易地确定。注意到在本实例中，在一个窗口对整个抽样层作了计算，在抽样层内薄层的实际时标位置(即深度)并不十分重要。如果薄层位于该抽样层内，则该窗口的谱线上将显示出移动凹槽的特征图案。本领域技术人员懂得，移动某参数的时间位置不会改变其振幅谱线。进而，它只引起相位变化，如果振幅谱线已被计算、检测过，则这种相位变化就不明显。

用动态平面视图显示振幅谱线的作法，也可以代之以实例与其它任何数量的参量一使用，这些参量通过存储于调谐立方体内的复数值算出。举例说明，复变换系数的相位提供了确定薄层参数的另一种方法，更一般地，也适用于岩块的水平间断性。相位调谐立方体如下计算：

P (j, k, i) = ta n^{- 1} [\frac{IM (A (j, k, i))}{RE (A ((j, k, i))}]

其中P(j，k，i)包含初始调谐立方体中截各点上傅里叶变换的复系数的相位部分。本领域技术人员很久以来都利用相位视图帮助找出特殊反射体，相位视图体现了震动数据的连续性。但在本实例中，能指出所在处岩块水平变动的是谱线相位反应水平方向上的非连续性，这类岩层中，薄层截尾是个典型例子。从动态平面视图上看到，临近水平层边缘处的相位值显得相当“不稳定”：表现出具有异常的一阶导数。这样，薄层边缘和更通常的岩块水平间断性(如断层、裂缝、不一致性、不整合性，等等)，表现出与周围相位值形成对比的相位，因此较容易确定。这个特性既可以用于自身，确定水平边界，也可与振幅谱调谐立方体衔接衔使用，确定所在点岩块变化的存在。

最后，本发明人预见，此处揭示的调谐立方体技术可以导致对震动反射数据的新认识。这个调谐立方体(无论包含相位还是振幅数据的)可经显示和检验，与地下的岩石成分，岩石性质，地层结构及岩层地质状况建立经验关系。同样地，存储于调谐立方体中的傅里叶变换值可以被进一步利用，产生在勘探工作中有用的新震动参量。仅举例说明，可从调谐立方体数据算出的参量包括谱线振幅或相位的平均值，及其它任意数量的参量。本发明中这一点的重要性，在下文中还有详尽论述。在地震解译工作中众所周知，地震反射体性质的空间变化，通常在经验上与蓄积岩性及溶质成分有关。由于反射特征产生这些变化的精确物理原理不易理解，实用中解译员通常计算各种震动参量，然后描点绘图，寻找有预言性价值的参量。从调谐立方体计算得来的参量，表示了对反射体特性的局部分析(实际上由短窗计算得出)，因此，对解译工作的进步可能极为重要。

按照本发明中的第二个方案，提供了一种突出薄层作用的方法，这时用到的离散傅里叶变换中，从涵盖了抽样层的窗口计算得到一系列的活动短窗傅里叶变换，并用新方式显示。图6中大致介绍了这种方法，其详细过程示于图8。从概念上讲，本实例可被看作产生了一系列前面揭示的这种类型的调谐立方体，对每个用户设定的傅里叶变换的窗口位置有一个调谐立方体。

同样，X(k，j，n)表示三维地震数据体，“L”表示给定活动窗口傅里叶变换的长度。在本实例中，“L”通常要明显短于抽样层长度N。同前一样，傅里叶变换窗口长度，不是按计算简便，而是为了描绘地下特定种类薄层参数的目的来选择的。举例说，一个合理的变换长度起始点，应该是刚刚能够跨越抽样层内最“厚”的薄层的长度。注意到在例如波形不是很紧密的情况下，有必要增加其最小长度。在后一种情况下，窗口的最小长度可增加至样本中所测子波的长度。

整数变量NS表示用于连续窗口的样本增量，例如，如果NS等于1，就对抽样层内各个可能的起始点进行短窗口傅里叶变换，相继的窗口只存在一个样本之差。如果NS等于2，则相继的两个窗口共用除两个数据值外的全部相同数据值，在抽样层内每相隔一个起始点处进行变换。

本实例中傅里叶变换系数如下求得。从一条特定地震道的抽样层顶部开始，对抽样层内每个可能位置，进行一系列长度为L的活动窗傅里叶变换。如图8所示，设整数变量“M”为副变量，表示瞬时活动窗口数量。开始时设M等于1，来标明第一个活动窗口的位置。

现在，对于震动数据子体x(j，k，i)内位置(j，k)处的震道，可提取第M个活动窗口的数据，移至临时存放区，该活动窗口开始于样本数字(M-1)*NS：

y(nl)＝x(j，k(M-1)*NS+nl)，nl＝0，L-1继而进行傅里叶变换。如前面所揭示的，在进行变换之前，权函数可选择加到数据上。对于固定的M值，将前述计算加于子体内的每条震道，将在这个特定窗口位置产生调谐立方体。同样，增大M，使整个数据体再次进行此种计算，得到另一个完整的调谐立方体。这一次是对前面窗口下面开始NS样本的窗口位置计算的。

现在可将傅里叶系数放入备用存储区，等待检测。上面引入的符号要做小的调整，以保证若干个窗口可能被加到各条单独震道上。令A_M(j，k，i)代表傅里叶变换系数集的值，这些变换系数是从第“M”个算得窗口位置的抽样层内全部震道上取得的。注意到为此级数划分到的存储空间已经大大增加了。这样，存储总量依赖于为每条震道算得的活动窗口数目，即NW；并且所能存储的字总量必须至少是NW L，J和K的乘积：

存储量＝(NW)(L)(J)(K)

如前注明的，A_M(j，k，i)完全有可能从未完整保存于RAM中，而是一部分在RAM中，其余部分则存于盘中。

用上面引入的级数符号，并再假设加权数据的傅里叶变换存于X(i)，震道(i，j)第M个窗口的变换系数存放于级数位置：

A_M(j，k，i)＝X(i)，i＝0，L/2

再次，存于A_M(j，k，i)众多调谐立方体中的单个频率薄片，在用以检验模拟薄层前，用图8揭示的方法较好地进行了定标。在各种情况下，水平的频率薄片被单个定标，以使其平均值等于特定常量，因此白化了谱线。

在对抽样层内地震道进行处理后，可对每个调谐立方体单独检查，寻找薄层效应的证据。如同以前一样，通过检验与不同频率对应的一系列水平薄片，在振幅谱线中即可确定薄层效应。进而，对于与每个窗口位置对应的调谐立方体。这步可以单独进行，因此就得到一种常规指示法，显出特定薄层参量的时间和空间延伸范围。

按照本发明第三方案，提供了一种用离散傅里叶变换显化薄层效应的方法，所采用的方式在上面第二个实例中描述过；但还包括有另外的步骤，即在显示和分析前将傅里叶变换系数组成一个单频率体。图7中大致介绍了这种方法。按在上述第二个实例中所揭示的，备用存储级数A_M(j，k，i)将装满傅里叶变换系数，并较好地进行了定标。

令F(j，k，m)表示从A_M(j，k，i)中抽取的单频体，其中有L/2+1个不同值(如果L为奇数则有(L+1)/2个)，每个值代表一个长度“L”的变换产生的傅里叶频率。与第“i”个傅里叶频率对应的值从A_M(j，k，i)取得如下：

F(j，k，m)＝A_m(j，k，I)，m＝1，NW，j＝1，J，k＝1，K实际上，级数F(j，k，m)可从概念上看作是从每个活动窗口值中取得并按短窗副变量M的递增顺序将它们叠放在一起构成的。

为了识别薄层而将数据这样组织的好处是，它提供了一种确定薄层时间和空间位置的方法。解释为，如前所述，抽样层内薄层的时标位置并不影响其反应：在抽样层顶部附近和底部附近的薄层都产生相同特征的振幅谱线。从确定薄层的角度看，这一点是有好处的，但从确定碳氢化合物储量的能力方面看，这样做却并不好，通常总是地层标高较高的情况比较好。

然而，在本实例中，相同频率薄片组成的体，是以“时间”作为其纵轴的：变量M是与沿地震道下行的距离大致对应的副变量。这种组织方法也提供了另一种用途，建立薄层参数所需的大约时间。

为了介绍，假设给定的薄层活动在其频率域内有10赫兹的凹槽。则包含该地层的每个傅里叶变换短窗都显示出同样的凹槽。检查由薄片组成的10赫兹体，会有包含这个凹槽的薄片域。这样，在常数频率体内检测连续的薄片，就有可能对所研究的反射体在时间上进行局部化。更重要的是，如果已知特定凹槽发生于10赫兹频率处，则可将10赫兹调谐立方体动态化，并作为辅助手段检测确定薄层的水平延伸范围：本频率调谐立方体中所观察到的凹槽的极限，确定了薄层的极限位置。

在前面讨论中，都是按照对传统的地震数据进行处理的过程来表述的。但是，本领域技术人员明白，这里描述的发明可以为其它应用领域带来好处，用以确定除碳氢化合物外的其它地下矿藏的位置。仅举例说明，这里描述的同样手段，可用于处理和/或分析多成分地震数据、切变波数据、大地电磁数据、串井测量数据、满波形声波测井或前述各种的数字仿真模型。简而言之，这里揭示的处理方法可应用于地球物理方面的单时间级数，带来许多好处；但是最好将其应用于空间相关的时间级数族，包括薄层问题。因此在以下说明中，本领域技术人员会明白，这里用到的“地震道”一词具有广泛的含义，可一般地应用于地球物理方面的时间级数。

借助于几个优选的实例及附带的图，这项创造性的工具已作了论述，但是，除本文提及的外，本领域技术人员在不违背本发明的主旨情况下，还可能对此进行各种变化和进一步的完善。本发明主旨的范围由下面的权利要求限定。

Claims

1.一种通过定位周期性频域凹槽来显示薄层的方法，其在处理包含分布在预定的地层体内的三维地震道集的地震数据中实现，该预定的地层体包括有助于产生、迁移、累积或呈现碳氢化合物的薄层，该三维地震道集由至少两个不同的地震线组成，每个地震道与所述地震线中的至少一个对应，该地震道包括至少由时间、位置及振幅表征的数字样本，该方法包括以下步骤：

(a)选择一部分所述地层体和其中包含的地震道和数字样本，以此在所述地层体内确定一个三维抽样层，所述选择的地震道在至少两个不同的地震线上与至少两个不同的预定的地震线对应；

(b)在所述三维抽样层内利用离散的正交变换对所选择的地震道的至少一部分和所述选择的数字样本进行变换，并从如此变换的地震道产生变换系数；

(c)将所述变换系数组织成为调谐立方体，用于地震勘探，所述调谐立方体包括所述变换的三维阵列；

(d)使用计算机制作和交互显示所述调谐立方体的变换系数的至少一部分，并标识一个或多个地下地质特征，该地下地质特征有助于在所述预定的地层体内产生、迁移、累积或呈现碳氢化合物。

2.根据权利要求1的方法，其中步骤(b)所述离散的正交变换：

(i)具有多个正交原函数，和

(ii)应用于含有所述数字样本的窗口，得到与所述正交原函数相关的多个变换系数。

3.根据权利要求2的方法，其中在利用所述离散的正交变换进行变换之前，要对所述含有数字样本的窗口应用权函数。

4.根据权利要求3的方法，其中所述权函数为高斯权函数。

5.根据权利要求2的方法，其中步骤(c)包括在所述调谐立方体内对所述变换系数进行定标的步骤。

6.根据权利要求5的方法，其中在所述调谐立方体内对所述变换系数进行定标的步骤包括如下步骤：

(i)从步骤(b)的原函数中选择一个；

(ii)选择至少两个与步骤(i)中所述原函数对应的变换系数；

(iii)计算从步骤(ii)中选得的全部变换系数的复数模；

(iv)由步骤(iii)中算得的全部变换系数的复数模计算统计值；

(v)由所述统计值计算标定值；

(vi)将步骤(ii)中所述变换系数乘以所述标定值。

7.根据权利要求6的方法，其中所述统计值为全部变换系数所得复数模的算术平均值。

8.根据权利要求1的方法，其中步骤(d)含有以下步骤：

(i)选择一个正交原函数；

(ii)选择至少两个与所述正交原函数对应的变换系数；

(iii)将这样得到的各变换系数定位；

(iv)在代表所述位置处显示所述变换系数。

9.根据权利要求8的方法，其中对于多个选定的正交原函数重复进行步骤(i)至(iv)。

10.根据权利要求1的方法，其中所述步骤(b)和(c)包括下列步骤：

(i)在所述抽样层内设定窗口，所述窗口有起始样本数字并含有多个数字样本；

(ii)在所述抽样层内利用离散的正交变换对空间相关的地震道的至少一部分进行变换，所述离散的正交变换由多个正交原函数表征，并用于步骤(i)中所述数字样本的所述窗口，以产生多个与所述正交原函数对应的变换系数；

(iii)将所述多个变换系数组成调谐立方体，这里所述调谐立方体及变换系数与步骤(i)的所述窗口的所述起始样本数字对应。