CN1153983C

CN1153983C - 地震信号处理和勘察的方法

Info

Publication number: CN1153983C
Application number: CNB961911727A
Authority: CN
Inventors: 库尔特・J・马富特; 库尔特·J·马富特; 鳌た铝; R·林恩·柯林; ・L・法默; 史蒂文·L·法默; ・S・巴霍利希; 迈克尔·S·巴霍利希
Original assignee: Core Lab Global N V
Current assignee: Core Lab Global N V; BP Corp North America Inc
Priority date: 1995-10-06
Filing date: 1996-09-30
Publication date: 2004-06-16
Anticipated expiration: 2016-09-30
Also published as: AR003802A1; AU7382996A; NO326418B1; CN1576885A; WO1997013166A1; US5930730A; EP0796442B1; EP0796442A1; MX9704213A; NO971801L; AU709621B2; CN1166207A; NO971801D0

Abstract

用于勘察烃类化合物的一种方法，一种图形表示和实现的物件。在本发明的一个实施例中，该方法包括步骤：采集3D地震数据；将数据分成相对小的三维单元列；确定每个单元内地震轨迹的相似性，沉降和沉降方位角；将每个单元内的沉降、沉降方位角和相似性显示成二维图形的形式。在一个实施例中，相似性是时间，单元内地震轨迹数目，单元内表现沉降和表现沉降方位角的一个函数。单元的相似性通过产生单元内的轨迹的多个相似量度并选择其中最大的量度来确定。此外，对应于单元内最大相似量度的表观沉降和表观沉降方位角，认为是该轨迹的真实沉降和真沉降方位角的估计。由色彩度，章度和亮度表示的一个彩色图被用于描述每个单元的相似性，真实沉降方位角和真实沉降。真实沉降位角在色彩标度上图形显示，真实沉降在章度标度上图形显示，相似性的最大量度在该彩色图的亮度标度上图形显示。

Description

地震信号处理和勘察的方法

本发明属于一般的地震勘察领域，特别是，识别三维结构和地层构造的方法和装置。

本申请是1995年10月6日提交的申请号为005,032的临时申请和属于Bahorich和Farmer的申请号为353,934、提交日为1994年12月12日的美国专利申请的部分继续。

在地震勘察中，地震数据是沿着在陆地上是地音探测器列或在水中是水下听音器光束导线构成的线(见图1中线10和11)采集。地音探测器和水下听音器起着接收传入地中并从地下岩石交界面反射回表面的能量的传感器的作用。能量在陆地上通常是由Vibroseis^车辆提供，该车通过以给定的间隔和频率在地表振动大地来传送脉冲。在水下，通常是用气枪源。返回地表的能量的细微变化常常反映出地层结构和贮油结构中流体内容的变化。

在完成三维(3D)地震勘察中，原理是类似的，但是，线和列布置得更密以提供更详细的地表下状况的覆盖。对于这种高密度的覆盖，有大量的数字数据需要在对其做出最后解释之前进行记录、存储和处理。处理需要大量的计算机资源和复杂的软件以增强从地表下收到的信号和对遮盖了信号的噪声进行静噪。

在处理完数据后，地质工作人员整理和解释这些3D数据立方体(见图2)形式的3D地震信息，该数据立方体有效地代表着地表下构造的一个显示。利用这一数据立方体，信息可以以各种形式显示出来。在选择的深度上可制作出水平时间切片图(见图3)。采用一个计算机工作站，解释人员还可以对场进行切片以研究在不同地震区域的贮油构造问题。采用地震或适当的数据还可制作在任何方向的垂直切片或横断面。反射层的地震采集可以被描绘出来，因此可产生时间水平图。时间水平图可转换成深度以提供一定水平上的一个真实的比例构造解释。

以往采集和处理地震数据的目的是图像化解释构造和地层的地震反射。但是，由于横断面图的地层构造的信息量有限，在传统的地震显示上常常难以探测到地层的变化。虽然对时间切片和横断面的处理提供了观察到更大的缺陷部分的一个机会，但是在没有缺陷反射记录的3D体积上识别缺陷表面还是困难的。

相关性是地震轨迹相似和不相似的一个量度。两个地震轨迹的相关性越增加，地震轨迹就越相似。在0到1的标度上指定一个相关性量度，“0”表示最缺少相似，而“1”表示完全相似(即两个一致，或时间上平移的轨迹)。对于多于两个的轨迹可按同样方法定义。

一种计算相关性的方法在属于Bahorich和Famer(转让给Amoco公司)的美国专利一，申请号353,934，申请日为1994年12月12日。与可让缺陷、凹槽、坍塌和其它沉积构造3D视化的阴影消除法不同，Bahorich和Famer设计的相关性在地震数据本身上运算。当在声阻中有足够的变化时，Bahorich和Famer发展的3D地震相关性立方体在描给地震缺陷中是极有效的。在地层高亮度细微变化中也是相当有效的(例如，弯曲分布的凹槽，点棒，峡谷，坍塌和潮水排放形式的3D图象)。

虽然Bahorich和Farmer发明的方法曾是非常成功的，但它有一些局限。Bahorich发明的一个非相关假设是地震信号的零平均假设。当相关窗超过地震波列的长度时，这一假设是近似真实的。对于含有10Hz能量成分的地震数据，这需要一个相当长的100ms的窗，这样的窗会混淆与较深和较浅的时间区有关的地层。缩短时间窗(例如32ms)产生高的垂直分辨力，但却由于地震波列增加了人工耗费。不幸的是，一个更严格的非零平均运行窗的横向相关处理的计算消耗是更高量级的。此外，如果数据被相关噪声污染，仅采用两个轨迹进行表观沉降的估计其噪声干扰相当大。

因此，需要一种能克服前面方法的缺点的方法和装置。特别是，希望能够改善分辨力和计算速度。此外，在存在相关噪声情况下改善沉降的估计也是我们非常想要的。

本发明公布了一种用于确定地表下构造、缺陷和地形轮廓的方法和完成方法的装置。在本发明的一个实施例中，该方法包括步骤：采集给定地层体积上的3D地震数据；将体积划分成一个相对小的三维单元列，在每个单元内至少有5个横向分开的通常是垂直的地震轨迹；确定每个单元内轨迹相对于两个给定方向的相似性；将每个单元的相似性显示成二维图形的形式。在一个实施例中，相似性是时间、单元内地震轨迹数目、单元内表观沉降和表观沉降方位角的一个函数；单元的相似性通过产生单元内轨迹的多个相似性的量度并选择其中最大的量度产生。此外，对应于单元内最大相似量度的表观沉降和表观沉降方位角认为是轨迹的真实沉降和真实沉降方位角的估计。最后，一种由色彩度、章度、亮度表示的彩色图被用于描绘每个单元的相似性、真实沉降方位角和真实沉降；特别是，真实沉降方位角图形表示在色彩标度上，真实沉降图形表示在章度标度上，最大的相似量度表示在彩色图的亮度标度上。

在本发明的另一个实施例中，公布的完成方法的装置包括一个计算机可读的具有让计算机完成地震探测处理的程序的媒体。在一个实施例中，计算机采集给定地震体积上的三维地震数据并且该媒体指令计算机完成下述工作：将该体积分成一个相对小的三维单元列，每个单元由其内部的至少5个横向分开的通常是垂直的地震轨迹表征；在每个单元内确定轨迹相对于两个给定方向的相似性；将每个单元的相似性存储起来用于显示成二维图形的形式。在一个实施例中，媒体程序将相似性定义成时间、单元内地震轨迹数目、单元内轨迹表观沉降和表观沉降方位角的一个函数；单元的相似性通过产生轨迹的多个相似量度和选择最大的量度来确定。此外，对应于单元内最大相似量度的表观沉降和表观沉降方位角被认为是该轨迹的真实沉降和真实沉降方位角的估计。该计算机包括产生由色彩度、章度和亮度表示的彩色显示的装置；该媒体具有将真实沉降方位角在色彩标度上图形表示、将真实沉降在章度标度上图形显示、将最大相似量度在亮度标度上图形显示的程序。

本发明的方法特别适合于解释3D地震体积内的缺陷平面和探测细微的3D地层构造。这是因为被缺陷线切断的地震轨迹通常具有与缺陷两侧轨迹不同的地震轨迹。沿着时间切片测量多通道相关性或轨迹的相似性揭示了沿这些缺陷线的低相关性的轮廓。这样的测量可揭示传统地震断面上不能马上明了的关键性的地表下细节。通过沿一系列的时间切片计算轨迹的相似性，这些缺陷轮廓还可识别出缺陷平面或表面。

本发明的方法提出了一种在噪声环境中通常比估计地震相关性的三轨迹横向相关方法更可靠的多轨迹相似性方法。此外，本发明提出的这一相似方法提供了：

对于高质量的数据，比三轨迹地震相关性的横向相关量度更高的垂直分辨力；

图形表示相关事件的3D实体角(沉降/方位角)的能力；

将复杂“轨迹”特性的概推广到复杂“反射层”特性的概念的能力；

通过将这些增强的复杂轨迹特性与相关性和实体角相结合，提供了适合于地质统计学分析方法的量化3D地震地层数据特性的基础。此外，相对于所选区域的沉降图的地震相关性可进行下列分析：

在详细选择之前结构和地层构架的分析

包括浅区、深区和主要感兴趣区附近的整个数据体积的结构和地层构造的分析。

不能被所选择的波峰和波谷表示的细微构造的分析。

所选构造或层系边界顶部和底部的内部中的构造分析。

与相关性相结合，相关地震反射事件的实体角沉降的数据立方体可让人迅速看到结构以及地震数据和被解释的层系边界之间的地层关系(如上复层和平移断层)。

从下面本发明的详细说明，其中的实施例、权利要求以及附图就可以明显看到本发明的众多的其它先进性和优点。

图1表示的是本发明从地表下采集3D地震数据用于处理的地音探测器的布置示意图。

图2是从采用图1布置采集的数据中获得的信息的图形表示。

图3是按以前的方法处理3D地震数据的一个水平时间切片(t＝1200ms)的图形表示。

图4A到4H显示了可能用于地震相关性，沉降和沉降方位角的运行窗分析的各种分析窗(计算星迹)；

图5是本发明采用以分析点为中心的椭圆分析窗进行处理的图形表示。

图6A和6B是当走向和沉降平行于采集轴的勘察分析以及当说明缺陷垂直于主反射面走向和沉降(p₀，q₀)时，所用的一个矩形沉降/方位棋盘格的例子。

图7A到7C是3个实体角沉降/方位角空间的图形表示。

图8A到8D描述了3D地震特性(φ，c，d)在3D色彩空间(H，L，S)上的图形表示。

图9显示了通过图8A的色彩半球的四个相关性值的四个表面。

图10到10C表示了图3地震数据的普通垂直切片。

图11A到11C表示了本发明方法应用于对应于图10A到10C的数据，所获得的地震特性、沉降、沉降方位角和相关性。

图12A和12B是穿过产生图11A和11B的沉降方位角立方体的时间切片(t＝1200ms和t＝1600ms)

图13A和13B是相关性的灰度显示。

图14A到14C表示了对应于图10A到10C的数据的相关性切片。

图15A和15B表示了本发明的一种相似算法和沉降/方位角算法的应用结果。

图16A和16B是完成本发明一个实施例的步骤的流程图。

虽然本发明容许有多种不同的实例，但图中所表示以及这里将要详细描述的只是本发明的一个特殊实施例。但是应当认识到，这里揭示的内容是从例举本发明原理的角度考虑的，并不想使本发明仅限于所例举说明的具体实施例。

本方法的第一步(见图16A)是获得一组分布于地层三维体积中的地震信号轨迹形式的地震数据。获得这些数据并将其缩减成数字形式以用于3D地震数据处理的方法对于那些熟悉此专业的人来说是皆知的。

相似处理

下一步是产生一个“相关立方体”。这一步骤是通过对3D地震数据作一多轨迹相似算法来完成。这一算法可有多种形式。无论什么形式，其作用是比较3D地震体积内相邻区域的地震数据的相似性。这个值(或特性)可作为地质构造中信号不连续性以及跨越断层和风化不一致区的信号不连续性的可靠估计。

我们定义了一个分析网格(或计算星迹)，该网格可以是具有“J”个轨迹的椭圆的或矩形的形式，并以一个给定的输出轨迹为中心(见图4A到4H)。

在该图中， X”表示分析窗的中心，而“○”表示在相似计算中所用的增加的轨迹。具有相同轨迹间隔(Δx＝Δy)用于分析数据的最小尺度的环形和矩形窗表示在图4A和4D中。横线/走向(Y)方向的轨迹间隔为顺线/沉降(X)方向的两倍(Δy＝2Δx)的分析数据用的最小环形和矩形窗表示在图4B和4E中。这样的不等距离间隔，通常用于勘察走向方向地质缓慢变化的情况。较大的分析窗用于反射沉降和方位角分辨力较高的情况，或者用于增加不完善数据区的信噪比，大的分析窗表示在图4C和4F中。

椭圆和矩形分析窗以分析点为中心，由一个主轴和短轴定义，主轴的方位表示于图4A和4H中。采集轴(x、y)从北-东轴(x′、y′)旋转φ0角度。这样的对称窗对裂纹探测非常有用。

如果我们把(x、y)轴的原点设在含J个地震轨迹u_j(t，x_j，y_j)的分析窗的中心，那么我们定义相似参数σ(τ，p，q)为：

σ (τ, p, q) = \frac{{Σ_{j = 1}^{J} u [τ - ({px}_{j} + {qy}_{j}), x_{j}, y_{j}]}^{2}}{J Σ_{j = 1}^{J} {u [τ - ({px}_{j} + {qy}_{j}), x_{j}, y_{j}]}^{2}} - - - (1)

其中三变量(τ，p，q)定义了时间τ时的一个局部二维事件，p和q为x和y方向以ms/m量度的表观沉降。由于p＝d sin φ，q＝d cos φ，这里d是真实沉降，φ是沉降方位角，因此得到：

u_j(τ，p，q，x，y)＝u_j[τ-d(x sin φ+y cos φ)，x，y].

熟悉本专业的人会认出，在方程(1)的分母J起着正则化因子的作用。分子代表平均能量，分母的总和项代表轨迹的总能量。事实上，方程(1)表示了相关能量与非相关能量的比

目的是完成顺线和横线方向表观沉降(P，q)的同步2维探测(见图5)。但是方程(1)给出的相似估计对于小的但相关的地震事件值是不稳定的，例如如果我们沿着一个平面相关波列的零交叉线求和就可能出现这种情况。为避免这一点，我们评估的时间T时的相关性c(τ，p，q)和表观沉降(p，q)是一个时间窗或者半长度K＝w/Δt样本2wms高度上垂直分析窗上的平均相似。

c (τ, p, q) = \frac{Σ_{k = - k}^{- k} {Σ_{j = 1}^{J} u [τ + kΔt - ({px}_{j} - {qy}_{j}), x_{j}, y_{j}]}^{2}}{J Σ_{k = - k}^{- k} Σ_{j = 1}^{J} {u [τ - kΔt - ({px}_{j} - {qy}_{j}), x_{j}, y_{j}]}^{2}} - - - (2)

总的来说，我们不知道但希望估计出与一个假定的3D反射事件局部沉降和方位角有关联的(p，q)值。

在本发明处理方法的一个实施例中，我们通过对所有可能的表现沉降进行强力探测来估计(p，q)(见图6A和6B)。我们假设程序可以从通常的地震数据(例如垂直数据切片)显示中估计出最大的真实沉降，d_max(以ms/m量度)，因此沉降限制为

\sqrt{p^{2} + q^{2}} \leq + d_{\max}

如果x_max和y_max是矩形分析窗的半宽和半长，并且f_max是地震数据中最高的短时的频率成分，那么每周期两点采样数据的Nyquist准则限制了表观沉降增量Δp和Δq如：

x_maxΔp≤1/(2f_max)，and y_max Δq≤1/(2f_max).

应当注意到Nyquist准则对于地震数据的线性运算有效，而方程(2)是非线性的。在实践中，我们发现Δp和Δq限制在Nyquist采样准则所要求的一半以获得一个相关倾斜事件的精确相似是必需的。

因此，我们对地震反射层的表观沉降的估计的研究减小到n_p*n_q

离散表观沉降对(p_l，q_m)上的相似c(P_l，q_m)的计算，其中：

n_p＝(2d_max/Δp)+1，

n_q＝(2d_max/Δq)+1

表观沉降对(p_l，q_m)在下式情况下可认为是反射层表观沉降的一个估计：

c(p，q)≥c(p_l，q_m) (3)

对于所有的-n_p≤l≤+n_p，-n_q≤m≤+n_q

估计的表观沉降

与估计的真实沉降和沉降方位

间的关系由下面简单的几何关系式表示：

\hat{p} = \hat{d} \sin \hat{φ}

\hat{q} = \hat{d} \sin \hat{φ}

其中以ms/m量度，角

从正x¹(或北)轴顺时针量度。当顺线采集方向x与N-S(x′)轴不一致时需做一φ₀角的简单座标旋转(见图4G)。

实体角离散和显示

优化角度离散化的重要性有两个原因：减少计算费用，限制能够用商业解释工作站软件显示的颜色的数目(例如目前的Landmaxk的“Seisworks”的64和Geoquest的“IESZ”系统的32)。

图7A表示的是在69角矩形网格中采用等增量Δp和Δq对表观沉降的离散。图7B表示的是在97角辐射状网格中采用等增量ΔdΔφ进行的离散。明显，我们不希望对10个不同方位角进行d＝0ms/m的沉降采样。图7C的“中国棋盘格”式网格更接近代表一个相等的因而更经济的具有最少点数(即61角)的(d，φ)表面采样。图7A和7C的每一网格代表一实角ΔΩ的近似的相等片。对于图7C所示角离散，环形分析半径为a，所选择的增量沉降Δd为

aΔd \leq \frac{1}{2 f_{\max}} . . . . . . . . (4)

显示

虽然独立地用图表示出相似性，沉降和方位角是可能的，但明显后两个特性是相互偶合的。此外，我们对这些估计的信心正比于相关性/相似性。其他人(见Bucher等的并转让给Amoco公司的美国专利4,970,699，“采色图形显示地质数据的方法”)已表明彩色HLS(色彩，亮度，章度)模型在显示多成份地震特性中是相互有效的(又见Foley，J.D和Van Dam，A.1981，交互作用图形基础，Addison-wesley，Reading.MA)。

参考图8A到8D，在这个方案中，我们直接将方位角φ图形表示在色彩轴H上：

H＝φ

其中H(通常称作“色轮”)和φ都在-180和+180度之间变化(见图8B)。蓝色对应于北，橙红色对应于东，黄对应于南，墨绿对应于西。对应于雾沉降的方位角是任意指定的一个0度(北)值，因此显示为蓝色。

下一个，我们将平均的相似/相关性c，图形表示在亮度轴L上(见图8C)：

L＝αc

其中0≤L≤100

0≤c≤1.0

α是一个小于100的比例常数，这是因为靠近L＝0(黑)和L＝100(自)的色彩和章度的变化难于区别。白色，或L＝100，对应于高相似或c＝1，而黑色，或L＝0对应于低相似，c＝0。中等的相似对应于中等的灰色调(例如如银色，灰色和碳灰)。亮度表示了照明量，代表了范围从黑到白的灰度。

最后，将沉降d图形表示在章度轴S上：

S＝100d/d_max

章度(S＝0)和色彩选择是任意的；我们可以将给定相似性的具有值(1+＝0，S＝100)的这一特性显示成白色，青蓝，纯蓝，深蓝和黑色。章度表示色彩缺少白光的冲淡。一个全章度的色彩没有白色加入，增加白色“冲淡”了色彩而没有改变其色彩(见图8D)。

图9表示的是通过图8A所示(φ，c，d)的3维(H，L，S)色彩半球的四个常数表面，对应于c＝1.00，c＝0.75，c＝0.5和c＝0.00。

附录1非常详细地描述了配色法。

HLS色彩模型的优点是：方位角是循环的并且图形的循环色彩轮平滑；对应于d＝0的方位无意义；对于浅沉降所有方位角光滑覆盖到灰色；弱区域沉降和方位角估计不准确的地方和低相似性(跨缺陷区)由深色彩表示。

数学方法的实现

例如，Landmark和Geo Quest解释工作站(见图16B)可通过加入工作为地震体积的处理过的数据用来观察和解释缺陷和地层构造。可视化软件(例如：Landmark的Seiscube软件)可用来对地震体积迅速的切片以帮助理解复杂的缺陷关系。

计算程序

编制了一FORTRAN77程序来完成计算和提供前述显示所需的信息。另外的详细情况在附录2中给出。每个轨迹U_MN通过其顺线和横线指标M和N进入。用户对输入数据中每个点/轨迹指定一矩形或椭圆空间分析窗或单元(见图4G)。分析窗的主轴和短轴，a和b规定为a＝aplength和b＝apwidth。主轴φ_a的方位规定为φ_a＝apazim矩形分析窗(图4H)由命令行指定的-R表示。与分析窗中心相关(并且对应于落在窗内的轨迹)的2J指标做成简单的表列，其中m_(j)，n_(j)分别表示x和y方向轨迹指标(相对于分析轨迹U_MN)。程序完成顺线和横线方向表观沉降(p，q)的同步2维搜索，其中(p²+q²)^1/2＜+smax。增量dp和dq的选择使得数据采样点在分析窗边界上每周期＝1/(fref)为四个点。为了解释，将每个表观沉降对(p，q)用真实(时间或深度)沉降d和沉降方位角φ的球座标系来表示是方便的。对试算的沉降和方位角(见图5)，将分析窗的数据插值到分段时间，τ-px-qy，其实质是“平滑”数据。分析点这一试算沉降的相似性定义成分析窗中这些平滑轨迹的相似。

对时间域的数据，我们按下式平滑分析点(M，N)的第j个轨迹：

u_f(τ，p，q，x，y)＝u[τ-(px+qy)]＝u(τ-d(xsinφ+ycosφ)].

其中x和y是从分析窗中心量起的距离。这可表示为：

u^f _{M+m(j)，N+n(j)}(τ，p，q)＝u_{M+m(j)，N+n(j)}[τ-(pn_(j)Δx+qm(j)Δy]

其中Δx和Δy是顺线和横线轨迹间隔。

对于深度域数据我们采用下式平滑第j轨迹：

u_f(ζ，p，q，x，y)＝u[ζ-(px+qy)]＝u[ζ-d(x sinφ+y cosφ)].

对所有后面的沉降和方位用下式计算相似：

σ (τ, p, q) = \frac{{(Σ_{j = 1}^{J} [u_{f} (t, p, q, x_{j}, y_{j})])}^{2}}{J Σ_{j = 1}^{J} [u_{f} (t, p, q, x_{j}, y_{j})]^{2}} - - - (5)

如同速度分析，每个沉降方位和分析点的相似通过形成一从-K到+K的部分和的运行窗时间积分来平滑，这里K＝apheight/dt。因此我们定义相关性，c(τ，p，q)为：

c (τ, p, q) = \frac{Σ_{- K}^{- K} {(Σ_{j = 1}^{- J} [u_{f} (t, p, q, x, y)])}^{2}}{J Σ_{- K}^{+ K} Σ_{j = 1}^{J} [u_{f} (t, p, q, x, y)]^{2}} - - - (6)

具有最大的(运行时间窗积分的)相关性c的沉降和方位角对Ω＝(d，φ)被认为是分析点处相关性c、沉降和方位角的估计。

实施例：

图11A到11C采用基于式(6)表示的相关性算法的相似方法和图8和9描给的色彩显示技术，显示了对应于图10A到10C的3维地震特性(φ，c，d)。输入数据的暂时采样间隔为4ms，顺线轨迹间隔Δx＝12.5m，横线轨迹间隔Δy＝25m，顺线采集方向沿南-北轴。对于图11A到11C，胶用了a＝b＝60m的环形分析窗或单元(见图4A)，以便包括计算中的总共11个轨迹。最大搜索沉降(见图7C)为d_max＝0.25ms/m，引起61个搜索角。采用的暂时积分时间是w＝16ms，或K＝4，因此平均化了9个样点相似计算。

在图10A和10B中，线AA’和BB’选作通过盐拱中心南到北和西到东的垂直切片。线CC’是南到北线的偏移，表示了垂直切片上径向缺陷的面貌。在图11A到11C中，盐拱内部由黑色代表，对应着普遍低相关性的区域。低相关性区域对应着线CC’上看到的径向缺陷。相关性，平沉降由浅灰色代表，主要位于盐拱以外的区段，特别是在线CC’上。盐拱北侧上的蓝色(在北-南线AA’上可见)对应着向北急剧沉降的沉积物(d＝d_max)。这些沉降从盐拱向外递渐变浅，因此在变缓之前，开始显示为蓝色(章度，S＝100.0)，藏蓝(S＝0.75)和钢蓝(S＝0.50)；最后显示为灰色(S＝0.0)。盐拱南侧的黄色(在AA’线可见)对应着向南急剧沉降的沉积物。盐拱东侧的橙红色(显示在东-西线BB’上)对应着向东急剧倾斜的沉积物。这些沉降从盐拱向外也逐渐变浅，一开始显示为橙红色(S＝100.00)，再就是赭色(S＝50.0)，最后为灰色，对应于平缓的沉降。最后，盐拱西侧森林绿(显示在AA’线上)对应于向西的急剧沉降的沉积物。这些沉降也从盐拱向外逐渐平缓，并和图9中图例西部表示的色彩来显示。北-南线CC’与盐拱径向不一致。因此，绘出了不同缺陷块由平面向外的旋转，其绿色块对应着向西南的沉降，青蓝色块对应着西北方向的沉降。

由于这些3维特性对输入地震体积上的每一点都计算，因此可显示成水平的特性时间切片(见图12A和12B)；这些对应于一个未处理的地震数据时间切片。盐拱内部以及径向缺陷显示成黑色，对应于数据非相关区。由于在t＝1200ms处盐拱的径向近似对称(见图12A)，位于拱侧的沉降沉积物也以一个水平的简单的方式向外辐射，使得其方位角非常接近图9左侧上的色彩图例。在t＝1600ms处，这一形式的对称性略少(见图12B)，这里向南要比向北的沉降更浅。此外，相关数据的内部块可从盐拱内部看到。

图9显示的色彩图例只允许四个相关“桶”。为了更详细检验相关性，可将相关性按单个特性绘出。这表示在图13A和13B中，其中所有184种色彩都应用于图8C显示的简单的灰度。在这一显示中，最大的相关性(c＝1.0)设为白色；最小的相关性(c＝0.0)设为黑色。虽然盐拱内部表示成一高非相关性区，但这一显示更好地表示出径向缺限形式的细微。特别的是，从盐拱辐射出的缺陷被显示出，其中我们去掉了一些分叉点。除相关特性更连续的重新分级外，感觉的部分差异是由于人类视网膜采用不同的接收器(锥形细胞相对于杆形细胞)来看彩色、黑色和白色。男性和女性之间在区分绿色和蓝色的能力上也存在着生理差异。因此，男性解释人员常常更喜欢图13A，13B和图15A表示的简单的单特性相关性显示，而不喜欢图11A到12B和图15B的多特性(φ，c，d)显示。事实上，这些显示是问侯性的：3D分量显示对于识别相邻旋转的缺陷块间相互冲突的沉降方位角的面貌有作用；单组分显示对于增强区分它们的边界或不相关缺陷的不连续性有作用。

处理方法上的考虑

仔细研究图13A和13B发现环状的非相关能量环绕着盐拱。为了研究这些现象的原因，将对应于图10A到10中的地震数据的单分量相关性立方体进行垂直切片。这些切片显示在图14A到14C。盐拱内部明显非相关。一个非相关的海底峡谷构造(参见Nissen等、“3D地震相关性技术在识别和描绘坍塌构造中的应用”，1995 SEG ExpandedAbstracts，1532-1534页)显示在盐拱的北部。如果图10A到10C所显示的地震数据被图14A到14C所显示的相关区覆盖，就可看到在图10A到10C中地震反射事件的零交叉情况下，图14A到14C的低相关性的区域之间是接近一致的。如果假定整个数据的地震噪声水平是固定的但不相关，那么这一现象很容易理解。对于表观沉降与地震反射层的强振幅的峰或谷重合的分析点(这样，信号能量的估计相对于不相关噪声是高的)，可以预想到信噪比是高的，导致高相关性的估计。但是如果分析点使得表观沉降与同样地震反射层的零交叉重合，则使信号与不相关噪声相比较低，可以预想到信噪比很低，导致相关性的低估计。

我们已发现三种增加信噪比的方法：第一种更适合于结构分析；第二种更适合于地层学分析；第三种适合于前二者。

对于急剧沉降的缺陷(从垂直向计小于45°)，可通过简单地增加方程(2)给出的垂直分析窗w的大小来增加信噪比。可以观察到两种效果。第一，对应于反射层零交叉点的结构漏失随垂直积分窗的增加而减少。第二，由于很少缺陷是完全垂直，缺陷的边界分辨力随垂直窗的增大而降低。w＝16ms的分析窗(可包含数据中峰值30Hz能量的一个完整循环)似互是一好的折衷方案。

第二个增加信噪比的方法(同等地适合于地层和结构分析)是沿一要解释的地层区域提取相关性。如果这一地层区域具有一地震数据的极值，例如峰或谷，只有那些具有相对高的信噪比的数据有选择地显示出来。显然，提取对应于零交叉点数据的相关性会大大地恶化相关性显示。这一方法的一个更经济的方案是首先沿感兴趣的区域平滑数据，然后只沿所选的区域计算地震特性。这一方法对自动(和人工！)挑选者在放大上有点更灵敏，因为选择中的循环跳过有些随机性，因此总是出现非相关性。

浅层构造(例如浅的沟槽；对应于再形成三角洲的浅冲积沟槽构造，小的阶梯中的缺陷)在要解释的地质层上面或下面不存在，因此，其所在地层上面或下面所包含的数据增加了无关联的幅值变化，从而使这些不连续性看起来更相关，因而变模糊。如果在解释区域上面或下面的时间样点含有独立的，或许是强的振幅不连续，假定地层区域含有不同地质时间产生的不同地层区域的混合构造，这些不连续将影响大的分析窗的分析。

第三种方法是将原始的地震轨迹u_j集推广成一分析轨迹V_j集，定义为：

v_j(t)≡v_j(t)+iu_j ^H(t)

其中u_j ^H(t)是u_j(t)的积分或者Hilbert变换，i表示 σ(τ，p，q)和c(τ，p，q)的计算完全与式(1)和式(2)类，这里注意V_j ²定义为：

v_j ²≡v_jv_j ^*≡(u_j+iv_j)(u_j-iv_j).

第三种方法避免了在式(1)的相似估计中在其它强反射层的“零交叉”处的数值不稳定性。

水平分析窗的作用

检查式(2)，很明显分析计算费用随分析包含的轨迹数线性增加。但是，通过比较11轨迹迹相关时间切片的相似性与3轨迹横相关时间切片的相似性，(这里每个都有一相同的W＝32ms的垂直分析窗)，可使人相信计算中增加更多的轨迹可增加信噪比。一般来说，信噪比随分析窗的增大而增大。但是，总体相关性有点减少(可看到较少的白色)，因为当增加窗尺寸时，由常数(p，q)平面事件确定的可能的曲面反射层的近似降低了。一般来说，在达到局部平面反射层近似不再起作用的一点之前，沉降/方位角的估计的信噪比是随计算中轨迹的数目增加的。

结论

本发明的3维相似技术提供了一种优秀的地震相关性的量度。通过采用任意大小的分析窗，我们能够平衡最大边界分辨力和信噪比之间的相矛盾的要求，这一点对于采用固定的三轨迹横向相关性技术是不可能的。精确的相关性量度可通过短时(垂直)积分窗获得，该窗与数据中最短周期同一量级，而零平均横向相关技术更偏向于采用比数据中最长周期更大的积分窗。因此，相似方法要比横向相关性方法产生更少的垂直向地质模糊，即使对于大的空间分析窗(见图15A和15B)。与相关估计同等重要，相似方法提供了一种估计每个反射事件的3D实体角(沉降和方位角)的直接手段。这些实体角图形可以与由构造边界定义的传统时间结构图有关，也可以无关。同Bahorich和Farmer的基本的相关性(例如横向相关)方法一样，在对用于地质沉降粗略观察的数据进行解释之前，就可获得瞬时沉降/方位角立方体的估计。在这种勘察模式中，相关性和瞬时沉降/方位角立方体可让用户早在项目解释阶段就可在重要的结构或沉积构造上挑出重要的沉降和走向线。在一种解释方法中，这些沉降和方位角可与构造和或层系边界有关，使人能够用图形表示三维结构内部的层次和超度形式。最后，估计了数据立方体中每一点的瞬时沉降和方位角，就可将常规的地震轨迹特性应用到局部平面反射层，因此大大地增加了信噪比。

从前面的描述，可以看到对那些精通此行的人明显可有各种变化，选择和修改。因此，这些描述只能当作是说明，是为了教那些精通此行的人实现本发明的方式。其它算法也可用来量度地震数据区附近的相似性，或产生“不连续立方体”。此外，相等的计算也可代替上面说明和描述的计算。例如，在沉降和方位角(d，φ)上的搜索可代替表观沉降P和Q的搜索。计算相似的逆可用来获得一种类似于照片页片的显示。本发明的一些特性也可独立于本发明的其它特性进行应用。例如，在估计了实体角(沉降和方位角)后，可获得一更光滑和更可靠的常规复杂轨迹特性的多轨迹估计(Taner，M.T.，Koehleo F.，和Sheriff，R.E；1979；“复杂的地震轨迹分析”；Geophysics，44，1041-1063)。代替计算单个轨迹的这些特性，还可计算分析窗内轨迹角堆积的特性。即可以计算：

a_i(τ，p，q)＝{[U(τ，p，q)]²+[U^H(τ，p，q)]²}^1/2，

ψ_i(τ，p，q)＝tan^-1{U^H(τ，p，q)/U(τ，p，q)}，

f_{i} = \frac{dΨ}{dτ} = \frac{U (τ, p, q) \frac{{&PartialD; U}^{H}}{&PartialD; τ} (τ, p, q) + U^{H} (τ, p, q) \frac{&PartialD; U}{&PartialD; τ (τ, p, q)}}{[U (τ, p, q)]^{2} + [U^{H} (τ, p, q)]^{2}}

b_{j} (τ, p, q) = \frac{| U (τ, p, q) \frac{&PartialD; U}{&PartialD; τ} (τ, p, q) + U^{H} (τ, p, q) \frac{{&PartialD; U}^{H}}{&PartialD; τ} (τ, p, q) |}{[U (τ, p, q)]^{2} + [U^{H} (τ, p, q)]^{2}}

其中：

U(τ，p，q)是

Σ_{j = 1}^{J} u [τ - ({px}_{j} + {qy}_{j}), x_{j}, y_{j}]}^{2}

(见式(1)分子)；

U^H(τ，p，q)是Hilbert变换，或U(τ，p，q)的积分分量

a_i(τ，p，q)是包迹或瞬时振幅；

Ψ_i(τ，P，q)是瞬时相位；

f_i(τ，P，q)是瞬时频率；

b_i(τ，p，q)是瞬时带宽(见Cohen，L.；1993，“瞬态分析”；Proc.IEEE intConf.Acoust，Speech Signal Processing，4，105-109)。

除这些“瞬时”特性外，还提出了其它一些特性用以表征轨迹包迹给定凸起内的信号为包迹τe峰值处该特性的信号。这些特性包括(见Bodine，J.H.；1994；“地震特性的波形分析”发表在SEG第54届国际年会，Atlanta，GA，USA)：

波列包迹：

a_r(τ，p，q)＝a_j(τ_e，p，q)

波列相位：

ψ_{r} (τ, p, q) = ψ_{j} (τ_{e}, p, q)

波列频率：

f_r(τ，p，q)＝f_i(τ_e，p，q)

波列带宽：

b_r(τ，p，q)＝b_i(τ_e，p，q)

零相位分量：

90度相位分量：

U^{90} (τ, p, q) = - Sin [ψ_{r} (τ, p, q)] U (τ, p, q) + Cos [ψ_{r} (τ, pq,)] U^{H} (τ, p, q)

以及偏斜度、上升时间和响应长度。由于混合沿实际沉降方向出现，事件的慢变化幅值，相位，频率和带宽成分应保留。此外，相关性/相似/相似性的计算可让人完成类似地震区域的“构造分析”。构造分析结合“群分析”导出区段分析。在其它应用中，这可让人制作地质相关性和外插地质表面特性。此外，相关性的确定可以用于对后堆迭和前堆迭地震反演施加一预先的限制。因此，应当认识到在不脱离所附权利要求所界定的本发明的精神实质的情况下，本发明可以有各种修正，选择，变化和改变。当然，所附权利要求中将要包括所有这样的涉及权利要求范围的修改。

附录1

多特性HLS校正

色彩是纯色的或100％章度的色彩，对应于下列1994年无毒性96色“Crayola”标准。

方向	φ(色彩)	Crayola色彩
方向	φ(色彩)	Crayola色彩	N	0	蓝色
NNE	30	铅色	N	0	蓝色
NNE	30	铅色	ENE	60	深红色
E	90	赭色	ENE	60	深红色
E	90	赭色	ESE	120	红色
SSE	150	桔红色	ESE	120	红色
SSE	150	桔红色	S	180	黄色
SSW	210	石灰绿	S	180	黄色
SSW	210	石灰绿	WSW	240	绿色
W	270	森林绿	WSW	240	绿色
W	270	森林绿	WNW	300	蓝绿色
NNW	330	蔚蓝色	WNW	300	蓝绿色
NNW	330	蔚蓝色	N	360	蓝色

部分50％章度对应于“不鲜明的”或“模糊的”色彩

方向	φ(色彩)	Crayola色彩
方向	φ(色彩)	Crayola色彩	N	0	军蓝
NE	45	暗褐色	N	0	军蓝
NE	45	暗褐色	E	90	褐红色
SE	135	深棕色	E	90	褐红色
SE	135	深棕色	S	180	金色
SW	225	橄榄色	S	180	金色
SW	225	橄榄色	W	270	茶绿色
NW	315	钢蓝色	W	270	茶绿色
NW	315	钢蓝色	N	360	军蓝

0％章度对应于无色彩

方向	φ(色彩)	Crayola色彩
方向	φ(色彩)	Crayola色彩	N	0	灰
E	90	灰	N	0	灰
E	90	灰	S	180	灰
W	270	灰	S	180	灰
W	270	灰	N	360	灰

低亮度值对应于“黑色”；中间亮度值对应于“深”色，高亮度值对应于“淡”色。

附录2

符号对照表\semb3d[-Nfile_in][-Ofile_out][-hlsfile_hls][-tstarttstart][-tend tend][-ildmdx][-cldmdy][-aplengthaplength][-apwidthapwidth][-apheightapheight][-apazimapazim][-ilazimxazim][-clazimyazim][-dzdz][-smaxsmax][-pminpmin][-pmaxpmax][-qminqmin][-qmaxqmax][-threshthresh][-freffref][-startlinastartline][-endlineendline][-exppower][-min][-int][-R]

说明

Semb3d读入3D地震后堆迭时间或深度数据，产生相似，沉降和方位角输出数据。

命令行变量

Semb3d从命令行变量中得到全部参数。这些参数规定了输入、输出、空间分析窗以及深度离散参数。下面的命令行变量用于本发明的一个实施例。

-Nfile_in

敲-N后立刻输入输入数据组名或文件。这些输入文件应包括完整的路径名，如果文件位于一不同的路径时。例如：-N/export/data2/san_juan/time_stack告诉程序在目录/export/data2/san_juan中寻找文件“time_stack”。对于本程序，数据是以正规的二进制数据的矩形网格形式存储。轨迹数(由行头词“Num_Trc”表示)由x方向轨迹数限定。记录数(地震线由行头词“NumRec”表示)限定了y方向轨迹数。插入在无信号轨迹中的遗漏数据由无信号轨迹头的标志来标识。

-Ofile_out

在敲入-O后立刻输入output多特性数据组名或文件，特性将被一行接一行地输回到后面。没有标定的相似C范围在0.0到1.0之间。沉降值范围在0到smax之间，总为正(指向下)。时间数据单位是msec/m(msec/ft)，或深度数据为m/m(ft/ft)。方位角φ垂直于走向方向，指向最大的正沉降(指向下)的方向。方位角的值范围在0°到360°之间。定义，0°方位角输出对应于北，90°方位角输出对应于东。OMEGA＝(d，φ)的值的选择应使(当转换成8位整数时)最左边的6位对应于有效的Seiswotks色彩表。这一色彩表对应于前面描述的HLS色彩模型并由一个将扫描进OMEGA＝(d，φ)的HLS(色彩，亮度，章度)色彩图的角度图形化的程度产生。

-hlsfile_hle

-hls后输入hls表文件名，以输出一含有输出文件中每个样本色彩亮度和章度的ascii平文件。这一文件输入到一个在一工作站上正常显示所需的产生RGB(红，绿，蓝)色彩寻找表的程序中。

-tstart tstart

-tstart后输入以msec计的分析窗的beginning。

-tend tend

-tend后输入以msee计的分析窗的end。输出记录将是(tend_tstart)msec long。

-ildm dx

-ildm后敲入以m(ft)计的顺线间距(轨迹间距)

-cldm dy

-cldm后输入以m(ft)量度的模线距离(线间距)。

-dz dz

-dz后输入垂直深度样点以m(ft)表示的增量。

dz＞0的值表示数据为深度。

-aplength aplength

-aplength后输入沿所用椭圆分析窗方位方向半口径的长度(以米或英尺量度)。通过增加aplength，apwidth来增大分析窗将导致：(1)增加了角分辨力，(2)减小了空间分辨力，(3)增加了计算费用，(4)减少了整体相关性(因平面波近似无效)。

-apwidth apwidth

-apwidth后输入垂直于所用椭圆分析窗方位垂直方向的半口径宽度(以米或英尺量度)

-apheight apheight

在-apheight后输入用于相似的运行时间(深度)积分窗的以微秒(或米或英尺)量度的半长度。例子＝±2样点。增加临时积分窗apheight将导致：(1)一光滑；较少噪声的响应，(2)减少垂直分辨力，(3)对计算费用无变化。

-apazim apazim

-apazim后输入椭圆分析窗(0为北，90为东)的方位角。

-smax smax

在-smax后输入对于时间数据以msec/m(msec/ft)表示或对于深度数据以m/m(ft/ft)表示的最大沉降。当数据中没有优先的走向时推荐这样做。这一值可直接从一个显示数据段读入。smax的值的量级对于时间数据为30msec/m(10msec/ft)。在真实沉降以外增加smax的值导致同样精度结果情况下计算费用的显著增加。

-pmin pmin

在-pmin后输入对于时间数据以msec/m(msec/ft)量度，对于深度数据以m/m(ft/ft)量度的最小顺线(增加的轨迹数)沉降。当有一个起主要作用的走向平行或垂直于数据采集线时推荐如此做。这一数据可直接从一个数据显示段读入。

-pmax pmax

-pmax后输入对于时间数据以msec/m(msec/ft)量度，对于深度数据以m/m或(ft/ft)量度的最大顺线(增加的轨迹数目)沉降。当有一个起主导作用的走向平行或垂直于数据采集线时推荐这样做。这一值可直接从一显示数据段读入。

-R

输入此命令语句定义-与椭圆分析窗相对应的方向为沿方位轴的矩形窗(2*aplength长2*apwidth宽)。

-qmin qmin

在-qmin后输入对于时间数据以msec/m(msec/ft)或对于深度数据以m/m(ft/ft)量度的最小横线(增加的行数)沉降。当有一起主导作用的走向平行或垂直于数据采集线时推荐这样做。这一值可从一显示数据段中直接读入。

-qmax qmax

在-qmax后输入对于时间数据以msec/m(msee/ft)或对于深度数据以m/m(ft/ft)量度的最大横线(增加的行数)的沉降。当有一起主导作用的走向平行或垂直于数据采集线时推荐这样做。这一值可从一显九据段中直接读入。

-thresh thresh

在-thresh后输入门限值或截止相似值，在该值以下，深度和方位角认为是有效的量测；在该值以下，将显示灰阴影。一些显示软件限制显示色彩的数目。

-fref fref

在-fref后输入对于时间数据以循环/秒(Hz)或对于深度数据以循环/KM(循环/kft)表示的参照频率以用于确定要搜索的沉降的数目(例如，fref＝60Hz对于时间数据，30循环/km对于深度数据)。

-ilazim ilazim

在-ilazim后输入增加的轨迹数的方位的顺线方位角(0度为北，90°为东)。这一值如果采用，用于校正实体角输出文件。

-clazim clazim

在-clazim后输入增加的线数的方位的横线方位角(0度为北，90度为东)。这一值，如果采用，用于校正实体角输出文件。

-exp power

-exp后输入用于相似的非线性换算的指数。一般，大部分相似/相关性值在0.8和1.0之间。用power＝2.0换算将使这些值变换到0.64和1.0之间，用power＝4.0将使这些值变换到0.41到1.0之间，等等。这对于将数据加入解释工作站有用。

-startline startline

在-startline后输入产生的第一个输出行。

-endline endline

在-endline后输入产生的最后一个输出行。

-min

在-min后输入提取对应于搜索的最小角相似的沉降，方位角和相似的命令语句。(作为缺省，程序搜索最大的相似和相关性)。

-int

输入这一命令语句是为了比例输出使得输出数据由一范围在-128到+127之间的8字节整数表示。这对将数据加入到一个解释工作站有用。

Claims

1.一种探测地下构造、缺陷和断面的方法，包括下列步骤：

(a)采集覆盖一预定地层体积的三维地震数据；

(b)将所述体积分成一相对小的三维单元列，每个所述单元由至少5个位于其内横向分开并且一般为垂直方向的地震轨迹表征；

(c)在每个所述单元内确定所述轨迹相对于两个给定的方向的相似性；以及

(d)以用于显示的方式记录所述单元内的所述相似性。

2.如权利要求1所述的方法，其特征在于在步骤(c)中，所述给定的方向相互垂直；并且在每个单元内所述轨迹的所述相似性至少是时间、所述分析单元内的地震轨迹数目以及所述分析单元内所述轨迹的表观沉降和表观沉降方位角的一个函数。

3.如权利要求2所述的方法，其特征在于每个单元内的所述轨迹的所述相似性由每个单元内所述轨迹的相似性量度的多个计算以及选择每个单元内所述相似性的所述量度的最大值来确定；步骤(c)还包括确定相应于所述量度最大值的表观沉降和表观沉降方位角以做为所述分析单元内地震轨迹的真实沉降的估计和真实沉降方位角的估计的步骤。

4.如权利要求3所述的方法，其特征是所述相似性的多个量度中的每一个至少是所述轨迹的能量的一个函数；并且所述轨迹的能量是时间、所述分析单元内地震轨迹的数目以及所述分析单元内所述轨迹的表观沉降和表观沉降方位角的一个函数。

5.如权利要求3所述的方法，特征在于所述分析单元内的所述相似性以图形形式记录，其中图形是由色彩章度和亮度表征的彩色图；

其中所述真实沉降方位角的估计，所述真实沉降的估计和所述最大计算相似量度中之一在亮度标度，色彩标度和章度标度中之一上图形表示出；

其中所述真实沉降方位角的估计，所述真实沉降的估计和所述最大计算相似量度中之另一个在亮度标度，色彩标度和章度标度中之另一个上图形表示出；

其中所述真实沉降方位角的估计，所述真实沉降的估计和所述最大计算相似量度中之剩余的一个在亮度标度、色彩标度和章度标度中的剩余的一个上图形表示出。

6.如权利要求5所述的方法，其特征在于，步骤(d)包括将真实沉降方位角的估计在所述色彩标度上图形表示，将所述真实沉降在章度标度上图形表示，以及将所述最大的计算相似性量度在亮度轴上图形表示的步骤。

7.如权利要求1所述的方法，其特征在于所述相似性至少是下式的一个函数：

{(Σ_{j = 1}^{J} u_{j} (t, x_{j}, y_{j}))}^{2}

和

Σ_{j = 1}^{J} u_{j} {(t, x_{j}, y_{j})}^{2},

其中每个分析单元含有至少J(J≥5)个地震轨迹的部分，x和y是从分析单元中心量起的距离，并且u_j(t，x_j，y_j)是所述分析单元内地震轨迹的一个部分。

8.如权利要求7所述的方法，其特征在于所述相似性是下式的一个函数：

\frac{{(Σ_{j = 1}^{J} u_{j} (t, x_{j}, y_{j}))}^{2}}{Σ_{j = 1}^{J} u_{j} {(t, x_{j}, y_{j})}^{2}} .

9.如权利要求8所述的方法，其特征在于所述相似性是下式的函数的一个算术反函数：

\frac{{(Σ_{j = 1}^{J} u_{j} (t, x_{j}, y_{j}))}^{2}}{Σ_{j = 1}^{j} u_{j} {(t, x_{j}, y_{j})}^{2}} .

10.如权利要求7所述的方法，其特征在于所述相似性通过在从-K到+K的部分和上进行一个运行窗时间积分来确定：

\frac{Σ_{k = - k}^{+ k} {(Σ_{j = 1}^{J} u_{j} (t + kΔt, x_{j}, y_{j}))}^{2}}{Σ_{k = - k}^{+ k} Σ_{j = 1}^{J} u_{j} {(t + kΔt, x_{j}, y_{j})}^{2}}

其中K是样点时间窗的半宽。

11.如权利要求1所述的方法，其特征在于所述相似性至少是下式的一个函数：

Σ_{j = 1}^{J} u_{j} {(t, p, q, x_{j}, y_{j})}^{2}

和

Σ_{j = 1}^{J} u_{j} {(t, p, q, x_{j}, y_{j})}^{2},

其中每个分析单元含有至少J(J≥5)个地震轨迹的部分，x和y是从分析单元中心量起的距离，p和q分别是x和y方向的表观沉降，并且u_j(t，p，q，x_j，y_j)是所述分析单元内地震轨迹的一个部分。

12.如权利要求11所述的方法，其特征在于所述相似性是下式的一个函数：

\frac{{(Σ_{j = 1}^{J} u_{j} (t, p, q, x_{j}, y_{j}))}^{2}}{Σ_{j = 1}^{J} u_{j} {(t, p, q, x_{j}, y_{j})}^{2}} .

13.如权利要求11所述的方法，其特征在于所述相似性是下式的函数的一个算术反函数：

\frac{{(Σ_{j = 1}^{J} u_{j} (t, p, q, x_{j}, y_{j}))}^{2}}{Σ_{j = 1}^{J} u_{j} {(t, p, q, x_{j}, y_{j})}^{2}} .

14.如权利要求11所述的方法，其特征在于所述相似性通过在从-K到+K的部分和上进行一个运行窗时间积分来确定：

\frac{Σ_{k = - k}^{+ k} {(Σ_{j = 1}^{J} u_{j} (t, p, q, x_{j}, y_{j}))}^{2}}{Σ_{k = - k}^{+ k} Σ_{j = 1}^{J} u_{j} {(t, p, q, x_{j}, y_{j})}^{2}},

其中K是样点时间窗的半宽。

15.如权利要求11所述的方法，其特征在于p＝0和q＝0。