CN106483559A - 一种地下速度模型的构建方法 - Google Patents

Abstract

本发明公开了一种地下速度模型的构建方法，属于油气勘探开发技术领域，避免了叠前地震同相轴的人工/半自动拾取的过程，加快了地下速度模型的构建过程。该地下速度模型的构建方法包括：获取共反射面元叠加算子；根据所述共反射面元叠加算子，获取地震波的运动学波场属性参数；基于建立的正则化层析方程，结合运动学波场属性参数和所建立的初始地下速度模型，构建地下速度模型。

Description

一种地下速度模型的构建方法

技术领域

本发明涉及油气勘探开发技术领域，具体地说，涉及一种地下速度模型的构建方法。

背景技术

估计地下介质的速度(尤其是其低波数部分，也称宏观速度)分布是地震勘探的核心问题之一，对地下构造偏移成像、预测油气储层参数都依赖于一个准确的速度模型。地震层析成像技术具有较高的反演精度，是估计地下速度模型的重要方法，其中应用最广泛的是反射层析成像。

反射层析具有效率高、容易实现等优点，然而其需要在叠前地震道集上拾取沿整个剖面连续分布的同相轴。这一操作十分耗时(特别是在三维情况下)，且对于低信噪比叠前地震数据而言，几乎是不可能完成的任务。

另外一类在工业中得到广泛应用的宏观速度估计技术是偏移速度分析(MVA)。该类技术将层析反演的原理与偏移成像生成的共成像点道集(CommonImage Gather，简称CIG)结合起来，降低了速度估计技术对叠前地震数据信噪比的要求，具有较高的实用价值。然而该类技术需要多次实施叠前深度偏移并人机交互拾取共成像点道集，这使得该类技术在三维情况下十分低效，一定程度上降低了其实用价值。

发明内容

本发明的目的在于提供一种地下速度模型的构建方法，避免了叠前地震同相轴的人工/半自动拾取的过程，加快了地下速度模型的构建过程。

本发明提供了一种地下速度模型的构建方法，包括：

获取共反射面元叠加算子；

根据所述共反射面元叠加算子，获取地震波的运动学波场属性参数；

基于建立的正则化层析方程，结合运动学波场属性参数和所建立的初始地下速度模型，构建地下速度模型。

可选的，基于建立的正则化层析方程，结合运动学波场属性参数和所建立的初始地下速度模型，构建地下速度模型包括：

步骤a、基于建立的正则化层析方程，结合运动学波场属性参数，计算模型更新量；

步骤b、将当前的地下速度模型叠加模型更新量，获得更新的地下速度模型；

步骤c、正演计算，获得更新的地下速度模型的运动学波场属性参数；

步骤d、判断更新的地下速度模型的运动学波场属性参数与地震波观测运动学波场属性参数之间的数据残差是否小于或等于预设残差，若是，执行步骤e，若否，执行步骤f；

步骤e、输出该更新的地下速度模型为所构建的地下速度模型；

步骤f、根据更新的地下速度模型，更新正则化层析方程，将更新的地下速度模型作为当前的地下速度模型，返回执行步骤a。

可选的，根据更新的地下速度模型，更新正则化层析方程包括：

获取数据与条件矩阵和更新的地下速度模型对应的弗雷歇(Fréchet)导数矩阵；

获取与地下速度模型空间维度相同的单位矩阵、对B样条基函数系数的更新量在横向和纵向上的一阶差分算子；

获取所设定的单位矩阵的乘数，获取两个一阶差分算子各自的乘数；

根据数据残差、弗雷歇导数矩阵、单位矩阵、对B样条基函数系数的更新量在横向和纵向上的一阶差分算子和所获取的乘数，更新正则化层析方程。

可选的，根据所述共反射面元叠加算子，获取地震波的运动学波场属性参数包括：

对零偏地震数据上的一样点，设定不同的运动学属性参数对，所述运动学属性参数对包括零偏射线由地下反射点传播至地表的角度与局部波前曲率半径；

利用所述共反射面元叠加算子，得到各运动学属性参数对在该样点的叠加能量谱；

获取叠加能量谱值最大的运动学属性参数对；

根据所获取的运动学属性参数对，基于射线理论，获取适合于后续建立正则化层析方程的地震波运动学波场属性参数。

可选的，共反射面元叠加算子涉及地表中点和半偏移距下地震波传播的走时关系式，该表达式不依赖于地下宏观速度模型。

可选的，运动学波场属性参数包括正入射点波到达地表的传播时间、波前曲率、横向坐标以及慢度矢量的水平分量。

本发明带来了以下有益效果：本发明实施例提供了一种地下速度模型的构建方法，该构件方法中利用共反射面元叠加算子来提取、观测地震波场的运动学波场属性参数，并用于层析速度反演，避免了叠前地震同相轴的人工/半自动拾取，理论上较反射层析更加适应于低信噪比叠前地震数据。

本发明的其它特征和优点将在随后的说明书中阐述，并且，部分地从说明书中变得显而易见，或者通过实施本发明而了解。本发明的目的和其他优点可通过在说明书、权利要求书以及附图中所特别指出的结构来实现和获得。

附图说明

为了更清楚地说明本发明实施例中的技术方案，下面将对实施例描述中所需要的附图做简单的介绍：

图1是本发明实施例中的构建速度模型的流程示意图；

图2是本发明实施例中的由s₀到s₁的透射过程的示意图；

图3是本发明实施例中的层析迭代的过程示意图；

图4是本发明实施例中的真实理论模型的示意图；

图5是本发明实施例中的初始速度模型的示意图；

图6是本发明实施例中的层析迭代1次后的速度模型的示意图；

图7是本发明实施例中的层析迭代10次后的速度模型的示意图；

图8是本发明实施例中的层析迭代20次后的速度模型的示意图；

图9是本发明实施例中的误差泛函下降曲线的示意图；

图10是本发明实施例中的部分炮道集的示意图；

图11是本发明实施例中的Kirchhoff叠前深度偏移结果的示意图。

具体实施方式

以下将结合附图及实施例来详细说明本发明的实施方式，借此对本发明如何应用技术手段来解决技术问题，并达成技术效果的实现过程能充分理解并据以实施。需要说明的是，只要不构成冲突，本发明中的各个实施例以及各实施例中的各个特征可以相互结合，所形成的技术方案均在本发明的保护范围之内。

本发明提供了一种地下速度模型的构建方法，如图1所示，该方法包括：

步骤S101、获取共反射面元叠加算子。

二维情况下，对地震数据中的某一零偏样点(ξ_m,T)，其中T表示正入射点波的传播时间，ξ_m是该零偏样点的地表中点坐标。共反射面元叠加算子是该零偏样点邻域范围内的中点坐标ξ与半偏移距h的函数，可由运动学属性参数对(α,R_NIP)表示如下：

其中，v₀为近地表速度，R_N是正入射点波的波前曲率半径。对于任一零偏样点(ξ_m,T)，其运动学波场属性d可通过对叠前数据进行多参数相干分析得到。式(1)反映了中点-半偏移距坐标系下地震波场的走时关系，可以看出共反射面元叠加方法只需要得知近地表速度，是一种完全数据驱动的方法。值得注意的是运动学属性参数对有明确的物理意义，其中α,R_NIP分别对应于零偏射线由地下反射点传播至地表的角度与局部波前曲率半径，对应于法向入射点(Normal IncidentPoint，简称NIP)波层析的数据分量，这里法向入射点波即为零偏射线。

步骤S102、根据所述共反射面元叠加算子，获取地震波的运动学波场属性参数。

接下来，对零偏地震数据上任一样点(ξ,T)，通过给定最大、最小值范围，等间隔地选取从而设定不同的运动学属性参数对(α,R_NIP)_i(i为任意正整数)。对每组运动学属性参数对，利用式(1)描述的共反射面元叠加算子进行叠加，得到各组运动学属性参数对的叠加能量谱。

对样点(ξ,T)，比较不同运动学属性参数对的叠加能量谱值，选取最大的叠加能量值对应的运动学属性参数对作为该点的运动学属性参数对(α,R_NIP)。

根据射线理论，可知有如下两个定量关系式：

其中，式(2)中的p^ξ为射线理论中慢度矢量的水平分量，代表了零偏数据的同相轴斜率，M为地震波场的局部波前曲率。

结合式(2)和(3)，对零偏地震数据上的任意一个样点可以得到一组运动学波场属性参数d＝(T,M,ξ,p^ξ)，即该运动学波场属性参数包括对NIP波到达地表的传播时间T、波前曲率M、地表中点横坐标ξ以及慢度矢量的水平分量p^ξ。本发明中，运动学波场属性参数的提取过程完全自动化，并作为运动学波场属性参数作为正则化层析方程的建立与求解的输入。

由于引入了同相轴斜率与波前曲率参数用来约束速度模型，本发明实施例中，反演得到的速度模型具有较高的精度，满足后续偏移成像的要求。

步骤S103、基于建立的正则化层析方程，结合运动学波场属性参数和所建立的初始地下速度模型，构建地下速度模型。

本发明的层析反演属于数据拟合类的速度估计方法，即找到一个如式(5)定义的速度模型m，使得正演计算得到的数据d_mod＝f(m)与式(4)定义的拾取数据d的误差最小化。f表示层析的非线性正算子，即反射点至地表方向的初值射线追踪。显然，对于一个零偏射线的反射，假如反射点的坐标、出射角以及速度模型正确时，由反射点以出射角出射射线，当到达观测面时，零偏射线刚好到以正确的出射方向、正确的走时到达正确的地表位置。

本发明实施例中的层析数据分量和模型分量如下：

d＝(T,M,ξ,p^ξ)_i,i＝1,...,n_data (4)

m＝[(x,z,θ)_i,i＝1,...,n_data；v_jk,j＝1,...,n_x,k＝1,...,n_z] (5)

式(4)中数据分别为拾取到的n_data个运动学波场属性参数。同样，每个数据点d对应于一个零偏反射过程，在模型中对应于式(5)的一组模型参数(x,z,θ)_i，即反射点横坐标、纵坐标、出射角。

以二范数衡量数据残差，本发明实施例中的构建地下速度模型的过程可归结为如下泛函的求极值问题：

其中，Δd(m)＝d-f(m)。对角阵也称为数据协方差矩阵，起加权不同数据分量的作用。

由于反问题固有的病态性与多解性，本发明实施例中，引入额外的先验模型进行约束，以保证反演结果的可靠性。改造式(7)的数据残差泛函如下：

其中为较小的正实数(这三个参数分布在0到1之间，这三个参数也需要用户通过测试自定义给出，当反演的初始模型(形式如式(5)定义)越好、数据品质越高则这几个参数给定越小，反之给定越大)，用于加权数据拟合项与规则化项之间的权重；m_r为某给定的先验模型，先验模型的形式与式(5)定义的速度模型一致，实施反演时需要用户自定义给出先验模型。给定先验模型的方法很多，本发明实施例取先验模型m_r＝0。

为阻尼项，这项的引入使得每次迭代的模型更新量较小。和分别约束速度模型在横向和纵向上的光滑程度，D₁、D₂分别是对B样条基函数系数更新量在横向和纵向上的一阶差分算子。

对一个有7个数的数列而言，其一阶差分算子定义为如下6×7的矩阵，

同理，对于有N个点的数列，一阶差分算子是形如上式的(N-1)×N矩阵。

根据反演的局部最优化理论，式(7)刻画的最优化问题可等效于如下线性方程组的迭代求解，通过多次迭代达到反演解收敛的目的，其中方程组(8)称为本发明实施例层析反演的正则化层析方程。

式(8)中，Δd是反演的数据残差，Δd(m)＝d-f(m)；I是一个与模型空间维度相同的单位矩阵其中N_{mod el}是本发明层析反演所用到的模型分量的维度；F是Fréchet(弗雷歇)导数矩阵，是一个N_data×N_{mod el}的矩阵，其中N_data是本发明实施例中层析反演所用到的数据分量的维度；D是数据预条件矩阵，用来加权本发明层析反演的不同量纲的数据空间分量，对每个运动学波场属性参数(T,M,ξ,p^ξ)_i，D由如下形式定义：

本发明实施例中取(σ_T,σ_M,σ_ξ,σ_p)＝(10³,10⁹,10,10⁶)。

其中f表示透射射线追踪算子(在油气地球物理勘探领域较为成熟)，m是式(5)定义的模型分量；F为在m_n处取值的Fréchet导数矩阵，其中F_ij表示Fréchet导数矩阵的第i行第j列的元素，具体如下：

本发明实施例中层析的Fréchet导数的计算基于射线扰动理论与射线传播矩阵的性质，考虑到射线扰动理论在射线中心坐标系下的简洁形式，本发明实施例在该坐标系下经过推导，得到了层析反演所需的Fréchet导数。

如图2所示，考虑如下一个由s₀到s₁的透射过程，有两边全微分得：

由式(1)即可计算出

的计算可以由射线扰动理论求得。射线扰动理论本质上是建立了观测射线场的扰动量、初始射线场以及速度的扰动量之间的线性关系，由该线性关系可以换算出(ξ,p^ξ,M)的扰动量与(x₀,z₀,θ,v)的扰动量的线性关系式。

射线扰动理论在射线中心坐标系下的形式最为简单，此处推导时先在射线中心坐标系下应用射线扰动理论，然后根据直角坐标与射线中心坐标的几何关系导出(ξ,p^ξ,M)与(x₀,z₀,θ,v)的线性关系式，即可求得

根据射线扰动理论，扰动坐标q与相应的慢度矢量分量p可由下式所示：

其中：Π(s₁,s₀)称为从s₀点传播到s₁点的射线传播矩阵，其为一2×2的矩阵，记为：

(Q₁(s₁,s₀),P₁(s₁,s₀))^T与(Q₂(s₁,s₀),P₂(s₁,s₀))^T分别为傍轴射线追踪系统初值为(1,0)^T与(0,1)^T的解，物理上分别对应于归一化的平面波震源与点震源。射线中心坐标系下的傍轴射线追踪系统为：

射线传播矩阵有如下性质：

Π(s₁,s)＝Π(s₁,s₀)Π(s₀,s)＝Π(s₁,s₀)Π^-1(s,s₀)

可知：

由A(2)、A(3)得：

其中，

考虑到射线中心坐标系与直角坐标系的几何关系有：

q₀＝cosθΔx-sinθΔz A(5)

s₀＝sinθΔx+cosθΔz

上述的几何关系A(5)代入A(4)得到：

另外，由得，

而根据射线理论，有：

S＝ΔS₁(Δv)+ΔS₂(Δq,Δp) A(10)

根据A(6)、A(7)即可求得导数联立A(5)、A(8)至A(12)可求得本发明实施例中用B样条插值表达速度模型，反演时在迭代中更新的速度参数为B样条基函数的系数，因此需要求取数据分量对B样条基函数系数的导数，这只需要将代入式A(6)、A(7)、A(11)即得到相应于B样条基函数系数的全微分表达式即可。其中，式A(12)中的v_q、v_qq、v_qqq分别代表速度模型在垂直于射线传播方向上的一阶方向导数、二阶方向导数与三阶方向导数。

本发明利用最小二乘QR方法(LSQR)求解矩阵方程组(8)(方程组(8)中Δm是未知量，其他均是已知量)，该方法是一种迭代的方法，可以在最小二乘意义下高效地求解大规模稀疏矩阵。

具体的，反演过程中，每轮迭代用运动学波动属性参数d来计算方程组(8)的右端项，其中f(m)是利用上一轮迭代的模型正演模拟得到的。

具体的，如图3所示，设当前的速度模型分量为m_k，则第k+1次迭代的过程为：

步骤a、基于建立的正则化层析方程，结合运动学波场属性参数，计算模型更新量Δm_k+1。

步骤b、将当前的地下速度模型m_k叠加模型更新量Δm_k+1，获得更新的地下速度模型m_k+1。

步骤c、正演更新的地下速度模型，获得更新的地下速度模型的运动学波场属性参数f(m_k+1)。

步骤d、判断更新的地下速度模型的运动学波场属性参数与地震波的运动学波场属性参数之间的数据残差是否小于或等于预设残差，若是，执行步骤e，若否，执行步骤f。

步骤e、输出该更新的地下速度模型为所构建的地下速度模型。

步骤f、根据更新的地下速度模型，更新正则化层析方程，将更新的地下速度模型作为当前的地下速度模型，返回执行步骤a。

获取数据与条件矩阵和更新的地下速度模型对应的弗雷歇导数矩阵；获取与地下速度模型空间维度相同的单位矩阵、对B样条基函数系数的更新量在横向和纵向上的一阶差分算子；获取所设定的单位矩阵的乘数，获取两个一阶差分算子各自的乘数；根据数据残差、以及上述获取的各参数，更新正则化层析方程。

显然，整个层析反演迭代过程全自动化，不需要进行人工干预，这提高了整个速度估计过程的效率，减少了人工工作量，缩短了处理周期。

基于上述构建速度模型的方法，在本发明实施例的具体实施场景中，图4为原始新疆岩性模型，模型大小为横向30公里、纵向6公里。为了验证层析对强非均质性模型的适应能力，对图4进行光滑并运用运动学与动力学射线追踪为其提供理想的数据输入。二维B样条基函数之间的间隔剖分为横向200m、纵向100m。图5为初始给定的常梯度模型，所前文所述，只需将B样条基函数系数在纵向给成相应的常梯度即可完成。由地表的数据分量在初始模型中由地表位置向下出射射线，当单程旅行时用完时停止，可以完成对反射点位置坐标和出射角的初始化。需要说明的是，该层析反演方法在理想无噪音数据中进行。

输入数据分量与初始模型分量，便可开始层析反演的模型迭代更新过程。图6为第一次迭代更新后的速度场，可以看出，表层的速度得到有效修改。随着迭代的进行，速度模型逐渐收敛，迭代10轮后的速度模型如图7所示，地表至四公里范围内的速度模型得到了较好的恢复，深部的模型由于采集孔径问题覆盖不足，未能完全恢复出来。迭代20次后的速度模型如图8所示，反演此时已基本收敛。

图9为误差泛函随着迭代次数增加的下降曲线，误差泛函取以10为底的对数，前15次误差泛函下降较快，后面的迭代中误差泛函下降十分缓慢，说明反演已收敛至极值处。本实验采用的规则化因子分别为ε_d＝0.02， 另外一种策略是在迭代过程中逐步下降规则化因子，这使得在保证反演稳定的同时更多的模型细节被恢复出来，但规则化因子也不能下降到某个指定的阈值一下，以防止反演再次出现不稳定。本实验说明本发明层析对于强非均质性模型具有较强的适应力。

图10为采用有限差分法对新疆岩性模型进行正演的炮道集(显示部分)。以反演第20次迭代的结果作为速度输入，实施Kirchhoff叠前深度偏移。偏移叠加图像如图11所示，可以看出基本的构造形态已得到较好地成像，也进一步验证了本文算法的精度。

综上，本发明实施例提供了一种地下速度模型的构建方法，该构件方法中利用共反射面元叠加算子来提取、观测地震波场的运动学波场属性参数，并用于层析速度反演，避免了叠前地震同相轴的人工/半自动拾取，理论上较反射层析更加适应于低信噪比叠前地震数据；同时，除反射地震走时外，将同相轴斜率与波前曲率用来约束地下速度模型，缓解了层析成像的射线多路径现象，具有更高的精度与反演稳定性。此外，本发明属于层析成像的范畴，只需要实施一次数据拾取操作，因此较偏移速度分析具有更高的计算效率。

虽然本发明所公开的实施方式如上，但所述的内容只是为了便于理解本发明而采用的实施方式，并非用以限定本发明。任何本发明所属技术领域内的技术人员，在不脱离本发明所公开的精神和范围的前提下，可以在实施的形式上及细节上作任何的修改与变化，但本发明的专利保护范围，仍须以所附的权利要求书所界定的范围为准。

Claims (6)

1.一种地下速度模型的构建方法，其特征在于，包括：

获取共反射面元叠加算子；

根据所述共反射面元叠加算子，获取地震波的运动学波场属性参数；

基于建立的正则化层析方程，结合运动学波场属性参数和所建立的初始地下速度模型，构建地下速度模型。

2.根据权利要求1所述的方法，其特征在于，基于建立的正则化层析方程，结合运动学波场属性参数和所建立的初始地下速度模型，构建地下速度模型包括：

步骤a、基于建立的正则化层析方程，结合运动学波场属性参数，计算模型更新量；

步骤b、将当前的地下速度模型叠加模型更新量，获得更新的地下速度模型；

步骤c、正演计算，获得更新的地下速度模型的运动学波场属性参数；

步骤d、判断更新的地下速度模型的运动学波场属性参数与地震波观测运动学波场属性参数之间的数据残差是否小于或等于预设残差，若是，执行步骤e，若否，执行步骤f；

步骤e、输出该更新的地下速度模型为所构建的地下速度模型；

步骤f、根据更新的地下速度模型，更新正则化层析方程，将更新的地下速度模型作为当前的地下速度模型，返回执行步骤a。

3.根据权利要求2所述的方法，其特征在于，根据更新的地下速度模型，更新正则化层析方程包括：

获取数据与条件矩阵和更新的地下速度模型对应的弗雷歇导数矩阵；

获取与地下速度模型空间维度相同的单位矩阵、对B样条基函数系数的更新量在横向和纵向上的一阶差分算子；

获取所设定的单位矩阵的乘数，获取两个一阶差分算子各自的乘数；

根据数据残差、弗雷歇导数矩阵、单位矩阵、对B样条基函数系数的更新量在横向和纵向上的一阶差分算子和所获取的乘数，更新正则化层析方程。

4.根据权利要求3所述的方法，其特征在于，根据所述共反射面元叠加算子，获取地震波的运动学波场属性参数包括：

对零偏地震数据上的一样点，设定不同的运动学属性参数对，所述运动学属性参数对包括零偏射线由地下反射点传播至地表的角度与局部波前曲率半径；

利用所述共反射面元叠加算子，得到各运动学属性参数对在该样点的叠加能量谱；

获取叠加能量谱值最大的运动学属性参数对；

根据所获取的运动学属性参数对，基于射线理论，获取适合于后续建立正则化层析方程的地震波运动学波场属性参数。

5.根据权利要求4所述的方法，其特征在于，

共反射面元叠加算子涉及地表中点和半偏移距下地震波传播的走时关系式，该表达式不依赖于地下宏观速度模型。

6.根据权利要求5所述的方法，其特征在于，

运动学波场属性参数包括正入射点波到达地表的传播时间、波前曲率、横向坐标以及慢度矢量的水平分量。