CN109116413A - 成像域立体层析速度反演方法 - Google Patents

Abstract

本发明提供一种成像域立体层析速度反演方法，包括：步骤1，对输入的初始速度场做叠前深度偏移，得到深度偏移剖面，并抽取角度域共成像点道集；步骤2，在角度域共成像点道集中，拾取与之对应位置处的剩余曲率与深度残差，并转换为走时残差；步骤3，计算立体层析数据空间；步骤4，在当前速度模型中进行射线追踪，得到模型中走时信息和反演核函数；步骤5，利用得到的立体层析核函数以及真实的数据空间与模型中数据空间的残差，建立反演方程组，并计算得到模型的更新量，完成此次迭代。该成像域立体层析速度反演方法可以获得更为准确的立体层析数据空间，提高了反演的精度与稳定性，在今后立体层析反演的实践中具有更好的利用价值。

Description

成像域立体层析速度反演方法

技术领域

本发明涉及油气勘探地震资料处理领域，特别是涉及到一种成像域立体层析速度反演方法。

背景技术

地震速度的反演是地震数据处理重要环节。立体层析反演是针对传统反射层析在低信噪比数据中同相轴拾取困难这一问题提出的。该方法重新定义了层析反演的数据空间和模型空间，即利用地震波走时，炮、检点位置以及炮、检点处射线出射方向同时约束速度模型来增强反演的稳定性和精度。在立体层析速度反演中，数据空间拾取的准确性是影响反演效果的最重要因素之一。通常立体层析速度反演所需的数据空间是基于数据域在共炮点与共检波点道集上对局部同相轴做倾斜叠加得到的。这是一种稳定的但是效率较低的数据空间拾取方法，对于二维数据具有可接受的工作量，但对于三维数据则很难接受。并且，利用叠前道集进行倾斜叠加，在交互拾取时对于一次波是有效的，但对于复杂构造中存在的绕射波无法分辨。在立体层析中，获取稳定准确的数据空间是保证反演准确性的前提。因此，针对数据空间拾取不准确与工作量大的问题。为此我们发明了一种新的成像域立体层析速度反演方法，解决了以上技术问题。

发明内容

本发明的目的是提供一种能有效的在成像域来反演速度场中低波数成分的成像域立体层析速度反演方法。

本发明的目的可通过如下技术措施来实现：成像域立体层析速度反演方法，该成像域立体层析速度反演方法包括：步骤1，对输入的初始速度场做叠前深度偏移，得到深度偏移剖面，并抽取角度域共成像点道集；步骤2，在角度域共成像点道集中，拾取与之对应位置处的剩余曲率与深度残差，并转换为走时残差；步骤3，计算立体层析数据空间；步骤4，在当前速度模型中进行射线追踪，得到模型中走时信息和反演核函数；步骤5，利用得到的立体层析核函数以及真实的数据空间与模型中数据空间的残差，建立反演方程组，并计算得到模型的更新量，完成此次迭代。

本发明的目的还可通过如下技术措施来实现：

在步骤2中，在叠前深度偏移剖面中拾取感兴趣的初始反射点位置(x₀,z₀)以及构造倾角ξ，并把初始反射点位置对应到角度域共成像点道集中；然后自动拾取剩余深度Δz以及剩余曲率tanφ。

在步骤2中，获取炮、检点处的射线参数p_r、p_s：

其中θ为入射角，θ_s、θ_r分别为射线在炮点、检波点处的出射角度，v为反射点处速度；

由走时残差与深度残差的转化示意图，求取走时残差为：

Δt＝2sΔz cosξcosθ

其中s表示慢度。

在步骤3中，由炮、检点处的射线参数p_r、p_s，得到炮点和检波点处射线出射时的水平慢度矢量p'_sx、p'_rx为：

p'_sx＝p_s·sinθ′_s

p'_rx＝p_r·sinθ′_r

其中θ′_s、θ′_r分别为射线分别在炮点、检波点处的出射角度。

在步骤3中，通过下面的修正公式对炮、检点处射线的出射信息进行修正，得到正确的

其中ψ＝cosξ·cosθ。

在步骤3中，由走时残差与深度残差的转化示意图，得到角道集中走时残差与深度残差的转换关系式为：

Δt＝2sΔz cosξcosθ

通过修正公式以及转换关系式，得到炮、检点坐标S、R，炮、检点处正确的地表射线出射信息以及正确的双程旅行时t。

在步骤4中，利用常速度梯度法进行射线追踪，并用哈密顿Hamiltonian系统对射线路径进行描述：

其中x表示射线传播位置信息，n表示射线的方向信息，t表示射线的传播时间，v表示介质中波的传播速度，利用龙格－库塔Runge-Kutta算法求解上式，得到模型中走时信息和反演核函数。

在步骤5中，选取最小平方准则建立目标函数：

其中Δd＝d_true-g(m)，C_D表示数据空间的协方差矩阵，d表示数据空间，m表示模型空间，g(m)表示数据空间与模型空间关系的非线性算子，由核函数表达，d_true表示真实数据空间，d_mod表示正演模拟得到的数据空间。

在步骤5中，根据目标函数，采用局部优化求取极值的方法，得到与其对应的反演方程：

其中G为核函数，利用最小二乘正交分解法来求解上述矩阵方程，得到的Δm即为该次迭代模型的更新量。

该成像域立体层析速度反演方法还包括，在步骤5之后，步骤6，判断是否进行下一次迭代。

在步骤6中，当判断需要进行下一次迭代时，流程返回到步骤1，重复此过程。

在步骤6中，当判断不需要进行下一次迭代时，跳出此循环，得到最终的层析速度场。

本发明中的成像域立体层析速度反演方法，针对传统反射层析在低信噪比数据中同相轴拾取困难这一问题，首先采用在偏移过程中包含地下射线角度信息的DSR(双平方根方程)叠前深度偏移方法提取出共成像点道集，然后从中获取剩余曲率与局部半偏移距射线参数，进而获得走时残差以及相对应炮点和接收点处的射线参数，再通过特定的修正公式，计算得到真实的立体层析数据空间。在初始的成像域道集中进行剩余曲率与构造倾角的拾取，并运用特定的修正公式修正得到更为准确的立体层析数据空间，从而建立立体层析反演方程，在成像域实现了立体层析速度反演，达到较高的速度反演精度。本发明利用此方法获得了更为准确的立体层析数据空间，提高了反演的精度与稳定性。这种方法将在今后立体层析反演的实践中具有更好的利用价值。

附图说明

图1为本发明的一具体实施例中立体层析的数据空间与模型空间的示意图；

图2为本发明的一具体实施例中走时残差与深度残差的转化示意图；

图3为本发明的一具体实施例中真实速度模型的示意图；

图4为本发明的一具体实施例中初始速度模型的示意图；

图5为本发明的一具体实施例中速度场中第400CDP点处速度为真实速度0.85倍时的角度域共成像点道集形态的示意图；

图6为本发明的一具体实施例中速度场中第400CDP点处速度为真实速度1.0时的角度域共成像点道集形态的示意图；

图7为本发明的一具体实施例中速度场中第400CDP点处速度为真实速度1.075倍时的角度域共成像点道集形态的示意图；

图8为本发明的一具体实施例中初始偏移剖面中初始反射点拾取和构造倾角ξ的示意图；

图9为本发明的一具体实施例中10次迭代后得到的速度场的示意图；

图10为本发明的一具体实施例中22次迭代后得到的速度场的示意图；

图11为本发明的一具体实施例中初始速度场偏移剖面的示意图；

图12为本发明的一具体实施例中立体层析迭代22次后速度场偏移剖面的示意图；

图13为本发明的一具体实施例中初始角道集的示意图；

图14为本发明的一具体实施例中立体层析速度场角道集的示意图；

图15为本发明的一具体实施例中真实速度场角道集的示意图；

图16为本发明的一具体实施例中距离在3.5km处层析速度、初始速度与真实速度对比图；

图17为本发明的成像域立体层析速度反演方法的一具体实施例的流程图。

具体实施方式

为使本发明的上述和其他目的、特征和优点能更明显易懂，下文特举出较佳实施例，并配合附图所示，作详细说明如下。

如图17所示，图17为本发明的成像域立体层析速度反演方法的流程图。

步骤101，对输入的初始速度场做叠前深度偏移，得到深度偏移剖面，并抽取角度域共成像点道集(ADCIG)。

步骤102，在偏移剖面中拾取初始反射点位置以及相应位置处的构造倾角，得到炮、检点处的射线参数。在角度域共成像点道集(ADCIG)中，拾取与之对应位置处的剩余曲率与深度残差，并转换为走时残差。

首先在叠前深度偏移剖面中拾取感兴趣的初始反射点位置(x₀,z₀)以及构造倾角ξ，并把初始反射点位置对应到角度域共成像点道集中。然后自动拾取剩余深度Δz以及剩余曲率tanφ。根据图1所示几何关系可以得到炮、检点处的射线参数p_r、p_s：

其中θ为入射角，θ_s、θ_r分别为射线在炮点、检波点处的出射角度，v为反射点处速度。

由图2所示走时残差与深度残差的转化示意图，可得走时残差为：

Δt＝2sΔz cosξcosθ

其中s表示慢度。

步骤103，利用修正公式以及步骤102中拾取到的构造倾角和剩余曲率信息，计算得到准确的立体层析数据空间。

由炮、检点处的射线参数p_r、p_s，得到炮点和检波点处射线出射时的水平慢度矢量p'_sx、p'_rx为：

p'_sx＝p_s·sinθ′_s

p'_rx＝p_r·sinθ′_r

其中θ′_s、θ′_r分别为射线分别在炮点、检波点处的出射角度。

然后通过下面的修正公式对炮、检点处射线的出射信息进行修正，得到正确的

其中ψ＝cosξ·cosθ

由图2所示走时残差与深度残差的转化示意图，可得角道集中走时残差与深度残差的转换关系式为：

Δt＝2sΔz cosξcosθ

通过修正公式以及转换关系式，自然可以得到炮、检点坐标S、R，炮、检点处正确的地表射线出射信息以及正确的双程旅行时t。

步骤104，在当前速度模型中进行射线追踪，得到模型中走时信息和反演核函数。

利用常速度梯度法进行射线追踪，并用Hamiltonian系统对射线路径进行描述：

其中x表示射线传播位置信息，n表示射线的方向信息，t表示射线的传播时间，v表示介质中波的传播速度。利用Runge-Kutta算法求解上式，得到模型中走时信息和反演核函数。

步骤105，利用得到的立体层析核函数以及真实的数据空间与模型中数据空间的残差，建立反演方程组，并计算得到模型的更新量，完成此次迭代。

选取最小平方准则建立目标函数：

其中Δd＝d_true-g(m)，C_D表示数据空间的协方差矩阵，d表示数据空间，m表示模型空间，g(m)表示数据空间与模型空间关系的非线性算子，由核函数表达，d_true表示真实数据空间，d_mod表示正演模拟得到的数据空间。

根据目标函数，采用局部优化求取极值的方法，可得到与其对应的反演方程：

其中G为核函数。利用最小二乘正交分解法(LSQR)来求解上述矩阵方程，得到的Δm即为该次迭代模型的更新量。

步骤106，根据角道集的拉平程度以及速度精度的要求，判断是否进行下一次迭代，如继续迭代返回步骤101，重复此过程；否则跳出此循环，得到最终的层析速度场。

在应用本发明的一具体实施例中，复杂断块模型是层析速度反演中的经典复杂模型，因此选择此模型来验证本文方法的可行性。复杂断块模型的真实速度场如图3所示，模型网格大小为650*550，横、纵向采样间隔分别为10米、5米。正演得到的炮记录有200炮，每炮361道接收。基于常规叠加速度分析方法得到叠加速度场并通过DIX公式进行转化得到初始层速度场，如图4所示。在基于成像域的立体层析数据空间的拾取过程中，最关键的是在角度域共成像点道集(ADCIG)中对剩余曲率(tanφ)和剩余深度(Δz)的拾取以及在共偏移距剖面中对构造倾角的拾取。

先对角道集进行分析，如图5-7所示，分别为速度场中第400CDP点处，速度为真实速度0.85倍、1.0倍和1.075倍时的角度域共成像点道集形态。由此可以看出，当速度存在误差时共成像点道集会出现弯曲的现象。并且当速度偏小时角道集向上弯曲，速度偏大时角道集向下弯曲。利用这一特征，在角度域共成像点道集(ADCIG)中即可拾取到剩余深度与剩余曲率。

如图8所示在初始偏移剖面中拾取初始反射点并用自动拾取的方法在上边拾取到构造倾角。这样通过步骤103中的修正公式以及转换关系式即可得到准确的立体层析数据空间。利用射线追踪方法进行正演得到立体层析反演核函数，并建立如步骤105所示的反演方程并求解，则可实现对初始速度场的多次迭代更新。利用本文方法，以图4作为初始模型记性速度场的迭代更新，进行10次、22次迭代后速度场结果如图9-10所示。

利用偏移剖面中偏移深度与真实深度的误差、绕射波等干扰波是否归位以及角道集的弯曲程度等(剩余曲率的大小)来判别是否进行下一次迭代。如图11-12所示分别为用初始速度场和迭代更新22次后得到速度场偏移得到的偏移剖的对比图，可以看出由于初始速度偏小，初始偏移剖面中界面深度向上移，剖面边界出现上翘现象，并且在一些陡倾构造处有较为严重的绕射等干扰波出现；而层析后得到的偏移剖面中，界面的深度与真实深度基本相同，绕射等干扰波也已经归位，可以清楚的看出地下介质的构造。同时抽取第310、320、330cdp处角道集，对比初始角道集、立体层析速度场角道集与真实速度场角道集(图13-15所示)，可以看出初始角道集中出现“笑脸”上翘现象；而利用本文提出的层析反演方法进行更新后的角道集同相轴已被拉平，并且与真实速度场偏移得到的角道集形态和同相轴深度基本一致。图16为第350cdp处初始速度、层析后速度与真实速度的对比图，可以发现层析速度与真实速度的吻合性良好。总体而言，本文提出的基于成像域的立体层析反演方法具有较高的精度，为后续的偏移成像及解释工作提供了较为准确的输入。

本发明中的成像域立体层析速度反演方法，在数据空间拾取过程中，在初始偏移剖面以及对应的共成像点道集中以每分钟50个点左右进行的拾取，本模型共拾取240个初始反射点，这种拾取效率是可以接受的；并且利用偏移剖面、角道集以及速度曲线把层析后的速度场与真实速度场进行了定性与定量的对比分析，发现得到了较好的反演结果，即可说明数据空间拾取的准确性。基于成像域的立体层析是由基于数据域的立体层析发展而来的，同时吸取了传统成像域的反射波走时层析优势并弥补了它的不足，使数据空间的拾取不需要严格沿着准确的同相轴进行，从而可以具有更高的反演精度。

Claims (12)

1.成像域立体层析速度反演方法，其特征在于，该成像域立体层析速度反演方法包括：

步骤1，对输入的初始速度场做叠前深度偏移，得到深度偏移剖面，并抽取角度域共成像点道集；

步骤2，在角度域共成像点道集中，拾取与之对应位置处的剩余曲率与深度残差，并转换为走时残差；

步骤3，计算立体层析数据空间；

步骤4，在当前速度模型中进行射线追踪，得到模型中走时信息和反演核函数；

步骤5，利用得到的立体层析核函数以及真实的数据空间与模型中数据空间的残差，建立反演方程组，并计算得到模型的更新量，完成此次迭代。

2.根据权利要求1所述的成像域立体层析速度反演方法，其特征在于，在步骤2中，在叠前深度偏移剖面中拾取感兴趣的初始反射点位置(x₀,z₀)以及构造倾角ξ，并把初始反射点位置对应到角度域共成像点道集中；然后自动拾取剩余深度Δz以及剩余曲率tanφ。

3.根据权利要求2所述的成像域立体层析速度反演方法，其特征在于，在步骤2中，获取炮、检点处的射线参数p_r、p_s：

其中θ为入射角，θ_s、θ_r分别为射线在炮点、检波点处的出射角度，v为反射点处速度；

由走时残差与深度残差的转化示意图，求取走时残差为：

Δt＝2sΔzcosξcosθ

其中s表示慢度。

4.根据权利要求3所述的成像域立体层析速度反演方法，其特征在于，在步骤3中，由炮、检点处的射线参数p_r、p_s，得到炮点和检波点处射线出射时的水平慢度矢量p′_sx、p′_rx为：

p′_sx＝p_s·sinθ′_s

p′_rx＝p_r·sinθ′_r

其中θ′_s、θ′_r分别为射线分别在炮点、检波点处的出射角度。

5.根据权利要求4所述的成像域立体层析速度反演方法，其特征在于，在步骤3中，通过下面的修正公式对炮、检点处射线的出射信息进行修正，得到正确的

其中ψ＝cosξ·cosθ。

6.根据权利要求5所述的成像域立体层析速度反演方法，其特征在于，在步骤3中，由走时残差与深度残差的转化示意图，得到角道集中走时残差与深度残差的转换关系式为：

Δt＝2sΔzcosξcosθ

通过修正公式以及转换关系式，得到炮、检点坐标S、R，炮、检点处正确的地表射线出射信息以及正确的双程旅行时t。

7.根据权利要求1所述的成像域立体层析速度反演方法，其特征在于，在步骤4中，利用常速度梯度法进行射线追踪，并用哈密顿Hamiltonian系统对射线路径进行描述：

其中x表示射线传播位置信息，n表示射线的方向信息，t表示射线的传播时间，v表示介质中波的传播速度，利用龙格－库塔Runge-Kutta算法求解上式，得到模型中走时信息和反演核函数。

8.根据权利要求1所述的成像域立体层析速度反演方法，其特征在于，在步骤5中，选取最小平方准则建立目标函数：

其中Δd＝d_true-g(m)，C_D表示数据空间的协方差矩阵，d表示数据空间，m表示模型空间，g(m)表示数据空间与模型空间关系的非线性算子，由核函数表达，d_true表示真实数据空间，d_mod表示正演模拟得到的数据空间。

9.根据权利要求8所述的成像域立体层析速度反演方法，其特征在于，在步骤5中，根据目标函数，采用局部优化求取极值的方法，得到与其对应的反演方程：

其中G为核函数，利用最小二乘正交分解法来求解上述矩阵方程，得到的Δm即为该次迭代模型的更新量。

10.根据权利要求1所述的成像域立体层析速度反演方法，其特征在于，该成像域立体层析速度反演方法还包括，在步骤5之后，步骤6，判断是否进行下一次迭代。

11.根据权利要求10所述的成像域立体层析速度反演方法，其特征在于，在步骤6中，当判断需要进行下一次迭代时，流程返回到步骤1，重复此过程。

12.根据权利要求10所述的成像域立体层析速度反演方法，其特征在于，在步骤6中，当判断不需要进行下一次迭代时，跳出此循环，得到最终的层析速度场。