CN102937721A - 利用初至波走时的有限频层析成像方法 - Google Patents

Abstract

本发明提供了一种初至波有限频地震走时层析成像的方法，包括：(a)通过进行野外勘探采集，得到原始采集的数据，并在原始炮道集上拾取初至波走时；(b)根据观测系统的高程与坐标信息建立初始速度模型；(c)在初始速度模型上进行动力学射线追踪，计算旅行时场与格林函数；(d)根据旅行时场与格林函数计算层析矩阵，并进行求解得到模型更新量；(e)对初始速度模型进行更新并判断是否满足反演精度。

Description

利用初至波走时的有限频层析成像方法

技术领域

本发明涉及地球物理勘探技术，更具体地讲，涉及一种能够提高表层速度反演精度的初至波有限频地震走时层析成像方法。

背景技术

表层速度结构是影响我国山地油气勘探的一个主要问题。剧烈的地表起伏会产生很多复杂的次生波动效应，如面波、转换波、多次波等。表层介质的强烈横向非均匀和小尺度异常体使波前断裂，同时产生大量散射波将有效反射信号湮没。如果没有合适的手段处理这些次生干扰会对地震数据处理的整个流程产生严重负面影响。这些处理手段包括表层静校正、基准面延拓、去噪等，但任何处理手段，甚至自起伏地表的偏移程序与速度分析都依赖于一个精确的表层速度模型。该表层速度模型同时可以通过照明分析辅助进行地震观测系统设计与优化。

现有的表层速度反演方法通常包括以下几种：

1、近地表调查手段，包括小折射、微测井等；

2、大炮初至折射层析；

3、基于射线理论的初至波走时层析；

4、基于射线理论的初至波井间层析；

5、基于射线理论的浅层反射波走时层析；

6、基于射线理论的浅层反射波井间层析；

7、全波形反演。

以上反演方法中，第1种反演深度浅，且只能得到有限区域的表层速度。第2种是目前广泛使用的方法之一，但该方法仅在有连续的高速折射层存在的情况下才能使用，且得到的速度是平均速度。第3、4种也是常用的表层速度反演方法，但它们基于高频射线理论，只能够得到介质的低波数成份信息，反演精度不高。第5、6种使用浅层反射信息，只在表层有连续的反射界面，且反射走时可以在地震记录中被准确拾取到的情况下才能适用。第7中反演精度与分辨率最高，但其对初始模型要求的强烈依赖性、对噪音的高敏感度、对传播算子的高要求限制了其在实际中的应用范围与效果。

发明内容

因此，本发明的目的在于克服在现有技术中的上述和其他缺点，提供一种获取具有高精度的表层速度反演方法。

为了实现上述目的，提供了一种初至波有限频地震走时层析成像的方法，包括：(a)通过进行野外勘探采集，得到原始采集的数据，并在原始炮道集上拾取初至波走时；(b)根据观测系统的高程与坐标信息建立初始速度模型；(c)在初始速度模型上进行动力学射线追踪，计算旅行时场与格林函数；(d)根据旅行时场与格林函数计算层析矩阵，并进行求解得到模型更新量；(e)对初始速度模型进行更新并判断是否满足反演精度。

附图说明

图1是示出根据本发明实施例的初至波有限频地震走时层析成像的方法的流程图。

具体实施方式

图1是示出根据本发明实施例的初至波有限频地震走时层析成像的方法的流程图。

如图1所示，在步骤S10，通过进行野外勘探采集，得到原始采集的数据，并在原始炮道集上拾取初至波走时。

初至波是激发源信号经过地层之后最先返回地面的信号，读出其在原始炮道集上最先返回地面信号的时间值就是初至波走时，采用能量比法或图像边缘检测方法可以自动拾取初至波走时。单道能量边界检测公式，即下面的等式1：

S_i＝|(B/A)×(B-A)|……等式1

其中，

，即，A为当前点以前同道的n个点能量之和，B为当前点以后同道的n个点能量之和。S_i为第i个样点求取的边界特征值。S_i越大代表i点越具有边界特征。对于单道记录而言，求取的边界特征值中最大值对应的点就是该道能量的边界。

在步骤S20，根据观测系统的高程与坐标信息建立初始速度模型。

初始速度模型对所考察区域进行均匀矩形网格剖分，并认为均匀矩形网格中的每个单元内的速度值是常数，这样，速度分布可由n维模型空间中的向量s^T＝(s₁，…，s_n)确定，其中分量下标表示单元编号，每个单元内的速度值可以由地层的厚度除以地震波在介质中传播的旅行时表示，一般是几百到几千米/秒的数量级范围。

在步骤S30，在初始速度模型上进行动力学射线追踪，计算旅行时场与格林函数。

在初始速度模型上沿着已知的一条中心射线路径，计算层析核函数，利用带有振幅层析核函数的等式2进行动力学射线追踪。

Δτ＝∫_VK_T(r)Δs(r)dr……等式2

其中，Δτ为有限频地震波传播的走时延迟，K_T(r)为走时层析核函数，V为中心射线对应的有限频率，Δs(r)为有限频内r处的慢度扰动。

当介质背景场比较平缓时可以采用等式4或5计算有限频(例如二维三维介质)层析核函数。否则可以采用式3计算有限频层析核函数。

$K_T(r,\mathrm{\ω})=\frac{2\mathrm{\ω}}{v_0\left(r\right)}\·\mathrm{Im}\left[\frac{G_0(g,r)u_0(r,s)}{u_0(g,s)}\right]$ ……等式3，

其中，ω为圆频率，Im表示取复数的虚部，G₀(g，r)为无扰动速度场v₀(r)中r点在g处的格林函数，u₀(r，s)为无扰动速度场v₀(r)中s点在r处的波场，u₀(g，s)为无扰动速度场v₀(r)中s点在g处的波场。

$K_T^{2D}(r,\mathrm{\ω})=\sqrt{\frac{t(g,s)\mathrm{\ω}}{2\mathrm{\π}{v}^{2}t(r,s)t(r,g)}}\sin (\mathrm{\ω\Δt}+\frac{\mathrm{\π}}{4})$ ……等式4

$K_T^{3D}(r,\mathrm{\ω})=\frac{t(g,s)\mathrm{\ω}}{2\mathrm{\π}{v}^{2}t(r,s)t(r,g)}\sin \left(\mathrm{\ω\Δt}\right)$ ……等式5

在二维、三维情况下，有限频的范围分别可以采用等式6和等式7进行计算。

$\mathrm{\Δt}=t(r,s)+t(r,g)-t(g,s)\≤\frac{3}{8}T$ ……等式6，

$\mathrm{\Δt}=t(r,s)+t(r,g)-t(g,s)\≤\frac{4}{8}T$ ……等式7，

其中，t(r，s)、t(r，g)分别表示激发点s与检波点g到空间位置r的绕射路径走时，t(g，s)表示激发点s到检波点g的中心射线路径走时。T为地震波主频周期。

格林函数在地震工程学中常用来计算震源机制。格林函数就是点源产生的场。在采用等式3时，格林函数的求取可以根据动力学射线追踪的结果，应用高斯束公式进行计算。由于其计算属于本领域的现有技术，因此不在此赘述。

根据初始速度模型网格大小与网格内的速度正演得到地震波在经过每个网格的走时，一系列的走时在整个研究空间内组成旅行时场。

在步骤S40，根据旅行时场与格林函数计算层析矩阵，并进行求解得到模型更新量。

其中，根据旅行时场与格林函数计算的层析矩阵如等式8所示：

$\underset{r}{\mathrm{\Σ}}\frac{2\mathrm{\ω}}{v_0\left(r\right)}\·\mathrm{Im}\left[G_0(g,r)u_0(r,s){\stackrel{\‾}{u}}_0(g,s)\right]\·\mathrm{\Δs}\left(r\right)\·\mathrm{ds}={\left|\right|u_0(g,s)\left|\right|}^2\mathrm{\Δ\τ}$ ……等式8

其中，ds为面元或体元的大小。

可以采用LSQR方法或者SIRT方法对等式8进行求解得到模型更新量Δs。

在步骤S50，对初始速度模型进行更新并判断是否满足反演精度。

利用Δs对模型进行更新，在上一轮计算结果上加上Δs作为这一轮的反演结果。设定一个门槛值作为判断是否达到反演精度的标准，当Δs小于门槛值时，说明达到反演精度，则计算结束。否则，返回步骤S30。

如上所述，根据本发明的初至波有限频地震走时层析成像方法继承了传统射线层析高效、稳定的优点，同时考虑了地震波传播的有限频效应，因此可以反演得到具有更高精度的反演结果。与传统射线层析相比，此处的高精度可以体现在两个地方。一是反演速度不再偏高，更加接近实际；二是反演结果的分辨率更高，可以得到介质的更高波数成份。传统射线层析中由于射线路径总体上是沿着高速介质传播，因此射线层析基本只能反演出高速异常体，模型更新量施加之后平滑会将高速扰动进一步弥散到低速区域，进而导致整体的反演速度偏高。而有限频层析中，由于有限频不再沿着高速区优势采样，有限频同样会将低速异常体包含其中，因此反演的结果不再偏高。射线层析基于高频假设的射线理论，导致其只有在异常体尺度大于有限频宽度时才有效，也就是说射线层析只能反演较大的异常体。有限频层析由于考虑了地震波传播的有限频特征，可以反演更小的异常体，因此有限频层析反演结果具有更高的分辨率。

根据本发明的初至波有限频地震走时层析成像方法可用于地震勘探中的非均匀各向同性速度结构反演。由于利用的是地震数据的初至波走时信息，因此反演区域不限于地表，具体反演区域受观测系统决定，如：观测系统在地表则可以反演表层速度结构，观测系统在井间则可以反演井间速度结构。反演的介质属性也不限于P波速度，具体反演属性由输入的初至数据信息决定，如输入的是P波初至走时信息则反演的是P波速度，输入的是转换波初至走时信息则反演的是转换波速度。

本发明的初至波有限频地震走时层析成像方法可被用于地震勘探中的观测系统设计与地震数据表层校正中，尤其是山地、沙漠、丘陵、黄土源、砾滩等具有复杂表层的探区。反演的表层速度应用于观测系统设计可以提高采集效率，提高表层照明均匀度。反演的表层速度也可以应用于静校正与基准面延拓等表层校正处理流程中，进而保障后续速度分析与叠加剖面的数据质量。反演的表层速度也可以应用于抗噪、起伏地表偏移成像与速度分析中。

尽管已经本发明的实施例具体显示和描述了本发明，但是本领域的技术人员应该理解，在不脱离由权利要求限定的本发明的精神和范围的情况下，可以对其进行形式和细节上的各种改变。

Claims (10)

1.一种初至波有限频地震走时层析成像的方法，包括：

(a)通过进行野外勘探采集，得到原始采集的数据，并在原始炮道集上拾取初至波走时；

(b)根据观测系统的高程与坐标信息建立初始速度模型；

(c)在初始速度模型上进行动力学射线追踪，计算旅行时场与格林函数；

(d)根据旅行时场与格林函数计算层析矩阵，并进行求解得到模型更新量；

(e)对初始速度模型进行更新并判断是否满足反演精度。

2.如权利要求1所述的方法，其中，在步骤(a)中，采用能量比法或图像边缘检测方法来拾取初至波走时。

3.如权利要求2所述的方法，其中，当采用能量比法被表示单道能量边界检测公式，如下面的等式1：

S_i＝|(B/A)×(B-A)|……等式1

其中，

，A为当前点以前同道的n个点能量之和，B为当前点以后同道的n个点能量之和，S_i为第i个样点求取的边界特征值，S_i越大代表i点越具有边界特征，对于单道记录而言，求取的边界特征值中最大值对应的点就是该道能量的边界。

4.如权利要求1所述的方法，其中，在步骤(b)中，初始速度模型对所考察区域进行均匀矩形网格剖分，并认为均匀矩形网格中的每个单元内的速度值是常数，速度分布由n维模型空间中的向量s^T＝(s₁，…，s_n)确定，其中分量下标表示单元编号，每个单元内的速度值由地层的厚度除以地震波在介质中传播的旅行时表示。

5.如权利要求1所述的方法，其中，在步骤(c)中，在初始速度模型上沿着已知的一条中心射线路径，计算层析核函数，利用带有振幅层析核函数的等式2进行动力学射线追踪：

Δτ＝∫_VK_T(r)Δs(r)dr……等式2

其中，Δτ为有限频地震波传播的走时延迟，K_T(r)为走时层析核函数，V为中心射线对应的有限频率，Δs(r)为有限频内r处的慢度扰动。

6.如权利要求5所述的方法，其中，当介质背景场比较平缓时采用等式4或5计算二维和三维介质的有限频的层析核函数，否则采用等式3计算有限频层析核函数：

$K_T(r,\mathrm{\ω})=\frac{2\mathrm{\ω}}{v_0\left(r\right)}\·\mathrm{Im}\left[\frac{G_0(g,r)u_0(r,s)}{u_0(g,s)}\right]$ ……等式3，

其中，ω为圆频率，Im表示取复数的虚部，G₀(g，r)为无扰动速度场v₀(r)中r点在g处的格林函数，u₀(r，s)为无扰动速度场v₀(r)中s点在r处的波场，u₀(g，s)为无扰动速度场v₀(r)中s点在g处的波场，

$K_T^{2D}(r,\mathrm{\ω})=\sqrt{\frac{t(g,s)\mathrm{\ω}}{2\mathrm{\π}{v}^{2}t(r,s)t(r,g)}}\sin (\mathrm{\ω\Δt}+\frac{\mathrm{\π}}{4})$ ……等式4

$K_T^{3D}(r,\mathrm{\ω})=\frac{t(g,s)\mathrm{\ω}}{2\mathrm{\π}{v}^{2}t(r,s)t(r,g)}\sin \left(\mathrm{\ω\Δt}\right)$ ……等式5

在二维、三维情况下，有限频的范围分别采用等式6和等式7进行计算，

$\mathrm{\Δt}=t(r,s)+t(r,g)-t(g,s)\≤\frac{3}{8}T$ ……等式6，

$\mathrm{\Δt}=t(r,s)+t(r,g)-t(g,s)\≤\frac{4}{8}T$ ……等式7，

其中，t(r，s)、t(r，g)分别表示激发点s与检波点g到空间位置r的绕射路径走时，t(g，s)表示激发点s到检波点g的中心射线路径走时，T为地震波主频周期。

7.如权利要求6所述的方法，其中，在采用等式3时，格林函数的求取根据动力学射线追踪的结果，应用高斯束公式进行计算。

8.如权利要求7所述的方法，其中，根据初始速度模型网格大小与网格内的速度正演得到地震波在经过每个网格的走时，一系列的走时在整个研究空间内组成旅行时场。

9.如权利要求8所述的方法，其中，在步骤(d)中，根据旅行时场与格林函数计算的层析矩阵如等式8所示：

$\underset{r}{\mathrm{\Σ}}\frac{2\mathrm{\ω}}{v_0\left(r\right)}\·\mathrm{Im}\left[G_0(g,r)u_0(r,s){\stackrel{\‾}{u}}_0(g,s)\right]\·\mathrm{\Δs}\left(r\right)\·\mathrm{ds}={\left|\right|u_0(g,s)\left|\right|}^2\mathrm{\Δ\τ}$ ……等式8

其中，ds为面元或体元的大小，通过对等式8进行求解得到模型更新量Δs。

10.如权利要求9所述的方法，其中，在步骤(e)中，利用Δs对模型进行更新，在上一轮计算结果上加上Δs作为这一轮的反演结果；通过设定一个门槛值作为判断是否达到反演精度的标准，当Δs小于门槛值时，说明达到反演精度，否则，返回步骤(c)。

Images