CN105445789A - 基于多次反射折射波约束的三维菲涅尔体旅行时层析成像方法 - Google Patents

Abstract

一种基于多次反射折射波来约束的三维地震菲涅尔体旅行时层析成像方法，包括：通过原始地震记录拾取初至波时间；建立离散模型；计算离散介质网格节点上的波前传播时间；向震源方向逐单元追踪射线路径；计算菲涅尔体；求解层析反演方程；将速度模型转换得到最终的层析速度模型；利用多次反射折射波旅行时公式来计算某控制点风化层厚度；计算多个控制点的风化层厚度，实现精确控制层析模型风化层厚度要求。其主要优势和特点包括：技术效果的可靠性。多次反射折射波可以精确求取风化层厚度，解决了层析反演在近地表多解性问题，使反演速度模型更加精确。操作简单易实现。该方法流程及计算简单，运算速度快，适合用于三维地震资料处理。

Description

基于多次反射折射波约束的三维菲涅尔体旅行时层析成像方法

技术领域

本发明属于油气勘探地震资料处理领域，是一种地震资料表层速度建模的有效方法。

现有技术

目前，近地表速度建模方法有多种，室外施工方法主要有小折射、微测井等，室内主要有基于勘探炮集记录初至的层析反演等，但这些表层建模方法在复杂表层下都有相应的局限性，如小折射施工要求地形平坦，有一个稳定的连续折射界面等，这在山地和黄土塬等复杂地区都是无法满足的；微测井尽管比较准确，但成本太高。

层析成像是目前精度较高的复杂表层速度建模方法，特别是菲涅尔体旅行时层析成像方法，其菲涅尔体路径相当于物理路径，更合乎地震波实际传播规律，提高了反演稳定性和分辨率，但该方法也是在线性近似的前提下反演得到的，其反演结果必然存在多解性或错误结果。

菲涅尔体层析成像主要包括初至拾取、初始模型建立和反演等步骤，成像的目的是获取有关的近地表信息，比如风化层厚度等，从而计算出准确的静校正量。从理论上讲，反演问题是一个非线性问题，菲涅尔体层析反演也不例外，但在实际计算中，各种反演方法都假设其为一线性模型，这样可以简化算法，减少工作量，但同时给反演带来多解性或错误结果。

多次反射折射波是复杂表层下一类主要的地震波动现象。特别是对于南方碳酸盐岩山区和西部山前带，由于风化层很薄，且风化层常与下覆的高速层地形走势一致，很容易导致多次反射折射波的产生。可以充分利用多次反射折射波旅行时公式准确计算风化层厚度，来约束菲涅体层析反演结果，以提高其准确性。

为解决层析反演的多解性问题，这里引入多次反射折射波，多次反射折射波是复杂表层下一类主要的地震波动现象，可以通过其旅行时公式来精确计算风化层厚度，约束菲涅尔体层析成像结果，该方法目前未见相应研究文献。

发明内容

本发明目的是针对现有技术存在的问题，提供一种精确计算风化层厚度，约束菲涅尔体层析成像结果的基于多次反射折射波来约束的三维地震菲涅尔体旅行时层析成像方法。

本发明具体技术方案：

(1)通过原始地震记录拾取初至波时间。

(2)建立离散模型，利用规则网格把表层模型离散成若干立方体单元，用网格节点处的速度值表示整个模型速度的变化。

(3)正向计算初至波旅行时，根据程函方程，使用基于水平集的GMM(GroupMarchingMethod)波前扩展算法，从震源开始，逐步计算离散介质网格节点上的波前传播时间。

(4)反向确定射线路径，利用正向步骤计算出的网格节点上的波前传播时间，从接收点开始，向震源方向逐单元追踪射线路径。

(5)计算菲涅尔体，先求出菲涅尔体半径，再以射线为中心轴线，用菲涅尔体半径确定菲涅尔体空间分布。下式是计算菲涅尔体半径的近似式：

$R\left(l\right)\≈\sqrt{\frac{l(L-l)}{\mathrm{fLs}\left(l\right)}}$ ①

式中，L是激发点与接收点之间的射线长度；f为地震波主频；l为从激发点到待计算菲涅尔体半径处的射线长度，s(l)为该处的平均慢度，R(l)为该处的菲涅尔体半径。

(6)求解层析反演方程，利用菲涅尔体空间及其能流密度的分布，建立初至走时残差与网格节点慢度修正量之间满足的线性方程。下式是每对激发点和接受点可获得的反演方程：

$\underset{j=1}{\overset{J}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\ω}}_k\frac{v_k}{V}{l}_{\mathrm{ij}}\mathrm{\Δ}{s}_k=\mathrm{\Δ}{t}_i$ ②

其中，Δt_i是第i条射线对应的走时残差；Δs_k是第k个单元的慢度增量，可以用该单元节点上的慢度增量表示，从而把该式整理成节点慢度增量与走时残差的关系式；ω_k是反映各单元能流密度大小的权系数，且ω_k可通过对一条菲涅耳体内各处的影响因子归一化得到。

(7)将第(2)步的速度模型换成第(6)步的速度模型，重复第(2)-(6)步多次(优选8-12次)，即可得到最终的层析速度模型。

(8)创新性的利用多次反射折射波旅行时公式来计算某控制点风化层厚度。多次反射折射波旅行时公式表示为：

$h_1=\frac{v_1\mathrm{\Δt}}{2\cos \mathrm{\θ}}=\frac{v_1\mathrm{\Δt}}{2\sqrt{1-{\left(\frac{v_1}{v_1}\right)}^2}}=\frac{\mathrm{\Δt}}{2\sqrt{\frac{1}{{v_1}^2}-\frac{1}{{v_2}^2}}}$ ③

其中：时间间隔Δt和表层结构参数v₁,v₂，可以通过微测井或多次反射折射波炮集记录来求取v₁,v₂。

上述方案中还包括：(9)根据层析模型需要，重复第(8)步，计算多个控制点的风化层厚度，实现精确控制层析模型风化层厚度要求。

发明的效果

采用基于多次反射折射波约束的三维地震菲涅尔体旅行时层析成像方法能够比较准确的获取近地表的速度模型，其主要优势和特点包括以下几点：

第一、技术效果的可靠性。多次反射折射波可以精确求取风化层厚度，解决了层析反演在近地表多解性问题，使反演速度模型更加精确。

第二、操作简单易实现。该方法流程及计算简单，运算速度快，适合用于三维地震资料处理。

附图说明

图1一种基于多次反射折射波约束的三维菲涅尔体旅行时层析成像方法实施例的流程图。

图2为菲涅尔体层析反演速度模型。

图3为基于多次反射折射波约束的菲涅尔体层析反演速度模型。

图4为纵线121线菲涅尔层析静校正后叠加结果。

图5为纵线121线基于多次反射折射波约束的菲涅尔体层析静校正后叠加结果。

具体实施方式

下面通过优化实施例对本发明技术方案做进一步说明。

(1)通过原始地震记录拾取初至波时间。

(2)建立离散模型，利用规则网格把表层模型离散成若干立方体单元，用网格节点处的速度值表示整个模型速度的变化。

(3)正向计算初至波旅行时，根据程函方程，使用基于水平集的GMM(GroupMarchingMethod)波前扩展算法，从震源开始，逐步计算离散介质网格节点上的波前传播时间。

(4)反向确定射线路径，利用正向步骤计算出的网格节点上的波前传播时间，从接收点开始，向震源方向逐单元追踪射线路径。

(5)计算菲涅尔体，先求出菲涅尔体半径，再以射线为中心轴线，用菲涅尔体半径确定菲涅尔体空间分布。下式是计算菲涅尔体半径的近似式：

$R\left(l\right)\≈\sqrt{\frac{l(L-l)}{\mathrm{fLs}\left(l\right)}}$ ①

式中，L是激发点与接收点之间的射线长度；f为地震波主频；l为从激发点到待计算菲涅尔体半径处的射线长度，s(l)为该处的平均慢度，R(l)为该处的菲涅尔体半径。

(6)求解层析反演方程，利用菲涅尔体空间及其能流密度的分布，建立初至走时残差与网格节点慢度修正量之间满足的线性方程。下式是每对激发点和接受点可获得的反演方程：

$\underset{j=1}{\overset{J}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\ω}}_k\frac{v_k}{V}{l}_{\mathrm{ij}}\mathrm{\Δ}{s}_k=\mathrm{\Δ}{t}_i$ ②

其中，Δt_i是第i条射线对应的走时残差；Δs_k是第k个单元的慢度增量，可以用该单元节点上的慢度增量表示，从而把该式可以整理成节点慢度增量与走时残差的关系式；ω_k是反映各单元能流密度大小的权系数，且ω_k可通过对一条菲涅耳体内各处的影响因子归一化得到。

(7)将第(2)步的速度模型换成第(6)步的速度模型，重复第(2)-(6)步10次左右，即可得到最终的层析速度模型。

(8)创新性的利用多次反射折射波旅行时公式来计算某控制点风化层厚度。多次反射折射波旅行时公式表示为：

$h_1=\frac{v_1\mathrm{\Δt}}{2\cos \mathrm{\θ}}=\frac{v_1\mathrm{\Δt}}{2\sqrt{1-{\left(\frac{v_1}{v_1}\right)}^2}}=\frac{\mathrm{\Δt}}{2\sqrt{\frac{1}{{v_1}^2}-\frac{1}{{v_2}^2}}}$ ③

其中：时间间隔Δt和表层结构参数v₁,v₂，可以通过微测井或多次反射折射波炮集记录来求取v₁,v₂

(9)根据层析模型需要，重复第(8)步，计算多个控制点的风化层厚度，实现精确控制层析模型风化层厚度要求。

下面结合附图和应用实施例更进一步说明。

本次研究对新疆某工区为靶区进行了应用，该工区为山前带，多次反射折射波发育，用该方法反演了近地表速度模型，以验证本方法的效果，具体流程见图1。

首先，依据步骤1，从实际资料求取初至时间。

其次，依据步骤2，根据炮点和接收点的分布确定速度模型的范围，在接收线方向10280m，在炮线方向3480m，在垂直方向1020m。用40×40×40的网格离散化该模型。地震波主频取为30Hz，初至模型取为速度为2000m/s的匀速模型。

然后，依据步骤3，采用程函方程正向计算初至波前走时。

接着，依据步骤4，利用正向步骤计算出的网格节点上的波前传播时间，从接收点开始直至到震源位置为止，获得相应的射线路径。

接着，依据步骤5，计算菲涅尔体。在正演计算得到射线路径后，沿射线路径逐点确定射线的菲涅尔体半径，到射线的垂直距离小于该半径的点组成菲涅尔体。

接着，依据步骤6，求解层析反演方程，得到速度模型。至此，第一次循环结束。

接着，依据步骤7，将步骤2的速度模型换成步骤6的速度模型，频率提高为40Hz，循环步骤2至步骤6，如此循环10次，每次循环将频率提高10Hz，就可以得到较为准确的速度模型。图2是第10次的反演出来的速度模型。从模型上可以看出，在桩号3928位置处风化层中有一个3米左右低速薄层，这给出一个疑问，是否准确？

接着，依据步骤8，在该桩号3928位置处，利用多次反射折射波旅行时公式来计算风化层厚度。表1是桩号3928处的各种方法速度模型参数表。因为微测井数据是实际测出来的，是相对真实可靠的，而多次反射折射波计算的风化层厚度和微测井数据相近，所以也是准确的。可以判定层析反演的速度模型可能由于计算误差或错误出现了一个高速薄层。

最后步骤9，依次计算多个控制点的风化层厚度，实现精确控制层析模型风化层厚度要求，得到精确速度模型如图3。

为了检验层析反演的近地表速度模型的有效性，我们利用菲涅尔体层析反演速度模型和基于多次反射折射波约束的菲涅尔体旅行时层析反演速度模型分别求取静校正量，并将静校正分别应用于实际数据，用来验证速度模型的准确性。图4是菲涅尔体层析静校正后纵线121线叠加结果，图5基于多次反射折射波约束的菲涅尔体层析静校正后纵线121线叠加结果，从图中可以看出，静校正效果改进明显，说明采用基于多次反射折射波约束三维菲涅尔体旅行时层析反演方法建立的近地表速度模型更加准确。

表1桩号3928处的各种方法速度模型参数

Claims (3)

1.基于多次反射折射波约束的三维菲涅尔体旅行时层析成像方法，其特征是包括：

(1)通过原始地震记录拾取初至波时间；

(2)建立离散模型，利用规则网格把表层模型离散成若干立方体单元，用网格节点处的速度值表示整个模型速度的变化；

(3)正向计算初至波旅行时，根据程函方程，使用基于水平集的GMM波前扩展算法，从震源开始，逐步计算离散介质网格节点上的波前传播时间；

(4)反向确定射线路径，利用正向步骤计算出的网格节点上的波前传播时间，从接收点开始，向震源方向逐单元追踪射线路径；

(5)计算菲涅尔体，先求出菲涅尔体半径，再以射线为中心轴线，用菲涅尔体半径确定菲涅尔体空间分布，下式是计算菲涅尔体半径的近似式：

$R\left(l\right)\≈\sqrt{\frac{l(L-l)}{\mathrm{fLs}\left(l\right)}}$ ①

式中，L是激发点与接收点之间的射线长度；f为地震波主频；l为从激发点到待计算菲涅尔体半径处的射线长度，s(l)为该处的平均慢度，R(l)为该处的菲涅尔体半径；

(6)求解层析反演方程，利用菲涅尔体空间及其能流密度的分布，建立初至走时残差与网格节点慢度修正量之间满足的线性方程，下式是每对激发点和接受点可获得的反演方程：

$\underset{j=1}{\overset{J}{\mathrm{\Σ}}}\underset{k=1}{\overset{K}{\mathrm{\Σ}}}{\mathrm{\ω}}_k\frac{v_k}{V}{l}_{\mathrm{ij}}\mathrm{\Δ}{s}_k=\mathrm{\Δ}{t}_i$ ②

其中，Δt_i是第i条射线对应的走时残差；Δs_k是第k个单元的慢度增量，可以用该单元节点上的慢度增量表示，从而把该式整理成节点慢度增量与走时残差的关系式；ω_k是反映各单元能流密度大小的权系数，且ω_k可通过对一条菲涅耳体内各处的影响因子归一化得到；

(7)将第(2)步的速度模型换成第(6)步的速度模型，重复第(2)-(6)步多次，即可得到最终的层析速度模型；

(8)利用多次反射折射波旅行时公式来计算某控制点风化层厚度，多次反射折射波旅行时公式表示为：

$h_1=\frac{v_1\mathrm{\Δt}}{2\cos \mathrm{\θ}}=\frac{v_1\mathrm{\Δt}}{2\sqrt{1-{\left(\frac{v_1}{v_1}\right)}^2}}=\frac{\mathrm{\Δt}}{2\sqrt{\frac{1}{{v_1}^2}-\frac{1}{{v_2}^2}}}$ ③

其中：时间间隔Δt和表层结构参数v₁,v₂，通过微测井或多次反射折射波炮集记录来求取v₁,v₂。

2.根据权利要求1所述的基于多次反射折射波约束的三维菲涅尔体旅行时层析成像方法，其特征是还包括：

(9)根据层析模型需要，重复第(8)步，计算多个控制点的风化层厚度，实现精确控制层析模型风化层厚度要求。

3.根据权利要求1或2所述的基于多次反射折射波约束的三维菲涅尔体旅行时层析成像方法，其特征是：步骤(7)重复第(2)-(6)步8-12次。