CN107229066B - 基于地面地震构造约束的vsp数据全波形反演建模方法 - Google Patents

Abstract

一种基于地面地震构造约束的VSP数据全波形反演建模方法，包括以下步骤：步骤S1：基于地面地震构造约束进行VSP数据全波形反演，获得中间速度模型M₁；步骤S2：针对中间速度模型M₁，进行没有地面地震构造约束的VSP数据全波形反演，得到最终速度模型M。采用了基于波动理论的全波形反演算子，相对于传统的VSP数据建模方式，建模精度更高，能够更加清晰的刻画模型的构造细节；引入地面地震构造约束滤波算子，不仅克服了传统全波形反演的数据不完备引起的建模噪音问题，以及低频缺失造成的反演陷入局部极值问题，而且在一定程度上解决了VSP数据角度覆盖有限的问题，提高了建模的角度范围。

Description

基于地面地震构造约束的VSP数据全波形反演建模方法

技术领域

本公开涉及石油地震勘探速度建模技术，尤其涉及一种高精度的基于地面地震构造约束的VSP(Vertical Seismic Profilling，垂直地震剖面)数据全波形反演建模方法。

背景技术

随着油气勘探的日益精细化、复杂化，高精度的速度建模在地球物理处理中扮演着越来越重要的角色。由于地震波速度不仅决定偏移成像的质量，还与岩石性质关系密切，能够反映岩石类别和富含流体情况，最终还可以影响地震解释结果的可信性，因而地震波速度是一个非常重要的地层物性参数。

近年来基于地震数据振幅、相位、波形匹配的全波形反演方法由于是基于有限频近似，未采用高频近似假设，因而理论上具有高精度刻画地下构造的能力，被业界认为是速度建模精度最高的一种方法。目前该方法在常规地面地震中已经得到了较好的应用，特别是海上地震资料的全波形反演已有大量的成功应用案例。但是采用VSP数据进行全波形反演以得到高精度速度模型目前还鲜有研究案例。

目前业内常规利用VSP数据对地下介质构造进行建模的方式有很多。在地震数据解释中，业内专家学者陈信平(1992)、周熙襄等(1992)采用零偏VSP的初至旅行时来计算不同地层间的层速度。当地下构造复杂时，采用这种方法求解的层速度误差较大。另外一种利用VSP数据进行速度建模的方法是旅行时层析。国外学者Stewart(1984)、Salo等(1989)、Lee(1990)详实系统地研究了如何应用VSP数据进行旅行时层析以获取地层层速度，此后国内学者刘清林(1996)采用多偏移距VSP数据完成了井间层析成像，朱健等(1999)利用VSP数据上行波的旅行时同时对层速度和界面形态参数进行了反演。黄光南等(2012)发展了变阻尼约束层析成像方法，并把该方法应用在了VSP数据速度建模当中，在一定程度上改善了射线不均匀覆盖带来的影响。由于VSP数据的角度覆盖有限，且射线层析类速度建模方案往往基于高频近似假设，因而采用VSP地震数据层析建模往往只能获取有限角度范围内的低波数背景速度，而无法精细刻画异常速度体的准确位置和构造形态。虽然全波形反演可以有效刻画模型的构造细节，但是当初始模型不准或数据不完备时，梯度会产生噪音，造成反演结果陷入局部极小值问题。

发明内容

本公开的目的是通过在VSP数据全波形反演的梯度中引入地面地震构造约束，从而有效地解决常规地面地震数据全波形反演面临的局部极小值问题，并在一定程度上突破了VSP数据角度覆盖局限问题。

本公开提供一种基于地面地震构造约束的VSP数据全波形反演建模方法，包括以下步骤：

步骤S1：基于地面地震构造约束进行VSP数据全波形反演，获得中间速度模型M₁；

其中所述步骤S1包括以下子步骤S11至S16：

子步骤S11：采用初始速度模型M₀以及地面地震数据进行叠前深度偏移，获得地下地质构造模式f(x)，其中x表示空间向量；

子步骤S12：构建所述地下地质构造模式f(x)的结构特征张量算子；

子步骤S13：基于时间域数据匹配的全波形反演构建VSP数据残差，基于所述VSP数据残差计算所述速度模型M₀的梯度g_v(x)，其中所述VSP数据残差是模拟的VSP数据与观测的VSP数据之差；

子步骤S14：基于各向异性扩散方程构建沿所述地下地质构造模式f(x)的结构特征张量算子平滑的地面地震构造约束滤波算子，对所述梯度g_v(x)沿所述地面地震构造约束方向进行平滑滤波，得到基于地面地震构造约束的速度模型梯度G_v(x)；

子步骤S15：通过P-L-BFGS算法更新所述速度模型M₀，获得更新的速度模型M₀’；

子步骤S16：判断是否完成给定次数的迭代，若是则将更新的速度模型作为中间速度模型M₁，并继续到步骤S2，否则返回子步骤S11，以更新的速度模型M₀’作为子步骤S11中的初始速度模型；

步骤S2：针对所述中间速度模型M₁，进行没有地面地震构造约束的VSP数据全波形反演，得到最终速度模型M。

优选地，采用克希霍夫偏移算子、单程波偏移算子或者逆时偏移算子进行叠前深度偏移。

优选地，所述结构特征张量算子包括所述地下地质构造模式f(x)在每个点的一阶对称半正定结构张量矩阵的特征向量和特征值。

优选地，采用声波正演模拟方程模拟VSP数据和正向传播波场，其中所述声波正演模拟方程表示为公式(3)：

其中，

p(x,t)＝[v_x(x,t),v_y(x,t),v_z(x,t),p(x,t)]^T代表正向传播波场，

s(x,t)＝[0,0,0,s(xs,t)]^T代表震源矢量，x_s为震源位置坐标，

b代表地下介质密度的倒数，

κ＝v²/b代表体积模量，

v代表地下介质的速度参数。

优选地，根据所述VSP数据残差求解逆时传播的伴随波场，以及基于所述伴随波场根据公式(4)求解所述速度模型M₀的梯度g_v(x)：

其中，NS表示震源个数，T代表波场传播的终止时刻，Q(x)是所述逆时传播的伴随波场的分量。

优选地，所述地面地震构造约束滤波算子表述为公式(5)：

其中，D(x)表示扩散张量，且扩散张量D(x)与所述结构特征张量算子的特征值一致，α代表光滑因子。

优选地，采用非精确线性搜索法计算对所述速度模型M₀进行更新的步长。

优选地，采用P-L-BFGS算法求解海塞矩阵的逆矩阵，并基于所述逆矩阵更新所述速度模型M₀。

优选地，所述步骤2包括以下子步骤S21至S23：

子步骤S21：基于所述时间域数据匹配的全波形反演算子构建VSP数据残差，基于所述VSP数据残差计算所述中间速度模型M₁的梯度；

子步骤S22：通过P-L-BFGS算法更新所述速度模型M₁，获得更新的速度模型M₁’；

子步骤S23：判断是否完成给定次数的迭代，若是则将所述更新的速度模型M₁’作为最终速度模型M并输出；否则返回子步骤S21，以更新的速度模型M₁’作为子步骤S21中的中间速度模型。

与现有技术相比，本发明的有益效果是：采用了基于波动理论的全波形反演算子，相对于传统的VSP数据建模方式，建模精度更高，能够更加清晰的刻画模型的构造细节；引入地面地震构造约束滤波算子，不仅克服了传统全波形反演的数据不完备引起的建模噪音问题，以及低频缺失造成的反演陷入局部极值问题，而且在一定程度上解决了VSP数据角度覆盖有限的问题，提高了建模的角度范围；地面地震构造约束的VSP数据全波形反演与不加地面地震构造约束的VSP数据全波形反演串联策略的实施，可以进一步提高VSP数据建模的精度。

附图说明

通过结合附图对本公开示例性实施例进行更详细的描述，本公开的上述以及其它目的、特征和优势将变得更加明显，其中，在本公开示例性实施例中，相同的参考标号通常代表相同部件。

图1示出了根据示例性实施例的基于地面地震构造约束的VSP数据全波形反演建模方法的流程图；

图2示出了根据示例性实施例的基于地面地震构造约束的VSP数据全波形反演建模方法的步骤1的详细流程图；

图3示出了根据示例性实施例的基于地面地震构造约束的VSP数据全波形反演建模方法的步骤2的详细流程图；

图4示出了根据示例性实施例的复杂模型的逆时偏移结果；

图5示出了根据示例性实施例的复杂构造模型偏移结果及其结构特征张量算子的方向性示意图；

图6示出了根据现有技术的复杂构造的真实速度模型；

图7示出了根据现有技术的常梯度初始速度模型；

图8和图9分别显示根据现有技术，采用VSP数据以常梯度初始速度模型进行传统全波形反演所得到的梯度结果和最终速度模型结果；

图10示出了根据示例性实施例的采用地面地震构造约束滤波算子滤波后的速度模型梯度；

图11示出了根据示例性实施例的基于地面地震构造约束的VSP数据全波形反演速度建模结果；以及

图12示出了根据示例性实施例的基于地面地震构造约束的VSP数据全波形反演联合无构造约束的VSP数据全波形反演的速度建模结果。

具体实施方式

下面将参照附图更详细地描述本公开的优选实施例。虽然附图中显示了本公开的优选实施例，然而应该理解，可以以各种形式实现本公开而不应被这里阐述的实施例所限制。相反，提供这些实施例是为了使本公开更加透彻和完整，并且能够将本公开的范围完整地传达给本领域的技术人员。

图1-3示出了根据示例性实施例的基于地面地震构造约束的VSP数据全波形反演建模方法的流程图。如图1-3所示，根据示例性实施例的基于地面地震构造约束的VSP数据全波形反演建模方法包括以下步骤：

步骤S1：基于地面地震构造约束进行VSP数据全波形反演，获得中间速度模型M₁

步骤S1包括以下子步骤S11至S16：

子步骤S11：采用初始速度模型M₀以及地面地震数据进行叠前深度偏移，获得地下地质构造模式f(x)

采用初始速度模型M₀以及地面地震数据进行叠前深度域偏移，可以得到地下地质构造模式f(x)，其中x表示空间向量。在进行第一次运算时，初始速度模型M₀可以选择本领域现有的速度模型，例如常梯度速度模型，在经过后续的迭代更新后，M₀是经过更新的速度模型。在进行叠前深度偏移时采用的成像算子可以是克希霍夫(Kirchhoff)偏移算子，也可以是单程波偏移算子或逆时偏移算子。图4显示根据示例性实施例的复杂模型的逆时偏移结果，经过逆时偏移可以清楚获取地下地质构造模式。

子步骤S12：构建地下地质构造模式f(x)的结构特征张量算子

基于图像结构特征理论，在图像学领域中，可以把子步骤S12中通过偏移成像得到的地下地质构造模式f(x)看作一幅图像，该图像具有方向性纹理结构，通过求解纹理结构可以得到结构特征张量算子，其包括地下地质构造模式f(x)的在每个点的一阶对称半正定结构张量矩阵的特征向量和特征值，用于表征图像的局部结构信息，即图像特征的主要变化方向及变化速率。

对于二维偏移成像数据体，它在每个点的一阶对称半正定结构张量矩阵表述为公式(1)：

根据矩阵的特征分解定理把结构张量矩阵进行特征值分解，得到公式(2)：

T＝λ_uuu^T+λ_vvv^T (2)

其中，λ_u≥λ_v≥0表示特征值，u，v分别为两个特征值对应的特征向量。特征值及特征向量刻画了局部图像的特征。当两个特征值都为0时，图像区域为常数；当两个特征值相等并都大于零时，表示图像为各向同性；当两个特征值不相等并均大于零时，表示图像为各向异性，且局部图像区域内存在主方向，较大特征值对应的特征向量表示图像梯度变化最快的方向，其与图像局部的线性特征方向相垂直，相对而言较小特征值对应的特征向量代表了图像局部的线性方向。图5显示了复杂构造模型偏移成像结果及其结构特征张量算子的方向性示意图。

子步骤S13：基于时间域数据匹配的全波形反演构建VSP数据残差，基于VSP数据残差计算速度模型M₀的梯度g_v(x)

采用声波正演模拟方程模拟VSP数据和正向传播波场，然后计算VSP观测数据和模拟的VSP数据在时间空间域中的数据残差，即VSP数据残差。其中声波正演模拟方程表述为公式(3)：

其中：

p(x,t)＝[v_x(x,t),v_y(x,t),v_z(x,t),p(x,t)]^T代表正向传播波场，

s(x,t)＝[0,0,0,s(xs,t)]^T代表震源矢量，x_s为震源位置坐标，

b为地下介质密度的倒数，

κ＝v²/b为体积模量，

v是地下介质的速度参数。

采用声波正演模拟方程模拟VSP数据后，可以计算VSP观测数据和模拟的VSP数据在时间空间域中的数据残差，即VSP数据残差。然后，根据VSP数据残差求解逆时传播的伴随波场，以及基于所述伴随波场求解速度模型M₀的梯度g_v(x)，速度模型M₀的梯度g_v(x)的表达式为公式(4)：

其中，NS表示震源个数，T代表波场传播的终止时刻，Q(x)是逆时传播的伴随波场的分量。

子步骤S14：基于各向异性扩散方程构建沿地下地质构造模式f(x)的结构特征张量算子平滑的地面地震构造约束滤波算子，并对子步骤S13中得到的速度模型M₀的梯度g_v(x)沿地面地震构造约束方向进行平滑滤波，得到基于地面地震构造约束的速度模型梯度G_v(x)

实际地震资料往往缺少道集或存在较强的噪音干扰，而且其低频成分信息往往不可靠，因而采用这种数据进行传统全波形反演时会使反演陷入局部极值，反演结果中出现较强的噪音，而这种噪音是每次迭代的梯度中的噪音引起的。

图8和图9分别显示了采用VSP数据以常梯度初始速度模型(图7)进行传统全波形反演所得到的梯度结果和最终速度模型结果，从图中可以看出梯度图像中存在较强的噪音干扰，而反演结果陷入了局部极值，与真实速度模型(图6)差距较远。

为了提升VSP数据反演结果的质量，改善有限的角度覆盖范围，需要对速度模型M₀的梯度g_v(x)进行预处理。在此预处理方式采用沿地面地震构造约束方向对梯度g_v(x)进行平滑滤波，地震构造约束方向即子步骤12中求解出的较小特征值对应的特征向量代表的局部线性方向。

地面地震构造约束滤波算子采用近似的各向异性扩散方程来构建，如公式(5)所示：

其中，D(x)表示扩散张量，为保证扩散方向沿着地面地震构造方向进行，扩散张量D(x)应与子步骤S12求解出的结构特征张量算子的特征值一致，G_v(x)表示经过滤波后的速度模型梯度，α表示光滑因子，α越大，原始图像沿构造方向扩散程度越大，滤波后的图像越平滑，当α＝0时代表扩散程度为0，此时并未对原始梯度做任何预处理操作。

经过双线性变换之后可得到近似各向异性方程的紧致矩阵算子表达式：

ΗG_v(x)＝g_v(x)

其中Η是与扩散张量有关的稀疏矩阵。

采用共轭梯度法求解上述线性方程组，可得到沿地面地震构造约束方向经过滤波的速度模型梯度G_v(x)。图10显示了采用地面地震构造约束滤波算子滤波后的速度模型梯度，其不仅消除了原始梯度的噪音，而且保留了地面地震的构造形态。

子步骤S15：通过P-L-BFGS算法更新速度模型M₀，获得更新的速度模型M₀’

在对梯度进行滤波之后，采用非精确线性搜索法求解对速度模型进行更新的步长，采用最优化理论的P-L-BFGS优化迭代算法近似求解海塞矩阵的逆矩阵，并基于该逆矩阵更新速度模型M₀，获得更新的速度模型M₀’。

子步骤S16：判断是否完成给定次数的迭代，若是则继续到步骤S2；若未完成给定次数的迭代，则返回子步骤S11，以更新的速度模型作为子步骤S11中的初始的速度模型。

图11为基于地面地震构造约束的VSP数据全波形反演速度建模结果。

步骤2：针对步骤1中获得的中间速度模型M₁，进行没有地面地震构造约束的VSP数据全波形反演，得到最终速度模型M

步骤2包括以下子步骤：

子步骤S21：基于时间域数据匹配的全波形反演算子构建VSP数据残差，基于VSP数据残差计算中间速度模型M₁的梯度；

子步骤S22：基于最优化理论，通过P-L-BFGS算法更新速度模型M₁，获得更新的速度模型M₁’；

子步骤S23：判断是否完成给定次数的迭代，将更新的速度模型M₁’作为最终速度模型M并输出；若未完成给定次数的迭代，则返回子步骤S21，以更新的速度模型M₁’作为子步骤S21中的中间速度模型。

子步骤S21至S23的细节与前述子步骤S13、S15和S16相同，只是基于中间速度模型M₁来执行，在此不再赘述。

图12为联合反演得到的最终速度建模结果，其速度模型的细节刻画更为清晰、更加准确。

上述技术方案只是本发明的一种实施例，对于本领域内的技术人员而言，在本发明公开的原理的基础上，很容易做出各种类型的改进或变形，而不仅限于本发明上述具体实施例的描述，因此前面的描述只是优选的，而并不具有限制性的意义。

Claims (9)

1.一种基于地面地震构造约束的VSP数据全波形反演建模方法，包括以下步骤：

步骤S1：基于地面地震构造约束进行VSP数据全波形反演，获得中间速度模型M₁；

其中所述步骤S1包括以下子步骤S11至S16：

子步骤S11：采用初始速度模型M₀以及地面地震数据进行叠前深度偏移，获得地下地质构造模式f(x)，其中x表示空间向量；

子步骤S12：构建所述地下地质构造模式f(x)的结构特征张量算子；

子步骤S13：基于时间域数据匹配的全波形反演构建VSP数据残差，基于所述VSP数据残差计算所述速度模型M₀的梯度g_v(x)，其中所述VSP数据残差是模拟的VSP数据与观测的VSP数据之差；

子步骤S14：基于各向异性扩散方程构建沿所述地下地质构造模式f(x)的结构特征张量算子平滑的地面地震构造约束滤波算子，对所述梯度g_v(x)沿所述地面地震构造约束方向进行平滑滤波，得到基于地面地震构造约束的速度模型梯度G_v(x)；

子步骤S15：通过P-L-BFGS算法更新所述速度模型M₀，获得更新的速度模型M₀’；

子步骤S16：判断是否完成给定次数的迭代，若是则将更新的速度模型作为中间速度模型M₁，并继续到步骤S2，否则返回子步骤S11，以更新的速度模型M₀’作为子步骤S11中的初始速度模型；

步骤S2：针对所述中间速度模型M₁，进行没有地面地震构造约束的VSP数据全波形反演，得到最终速度模型M。

2.根据权利要求1所述的基于地面地震构造约束的VSP数据全波形反演建模方法，其中采用克希霍夫偏移算子、单程波偏移算子或者逆时偏移算子进行叠前深度偏移。

3.根据权利要求1所述的基于地面地震构造约束的VSP数据全波形反演建模方法，其中所述结构特征张量算子包括所述地下地质构造模式f(x)在每个点的一阶对称半正定结构张量矩阵的特征向量和特征值。

4.根据权利要求1所述的基于地面地震构造约束的VSP数据全波形反演建模方法，其中采用声波正演模拟方程模拟VSP数据和正向传播波场，其中所述声波正演模拟方程表示为公式(3)：

其中，

p(x,t)＝[v_x(x,t),v_y(x,t),v_z(x,t),p(x,t)]^T代表正向传播波场，

s(x,t)＝[0,0,0,s(x_s,t)]^T代表震源矢量，x_s为震源位置坐标，

b代表地下介质密度的倒数，

κ＝v²/b代表体积模量，

v代表地下介质的速度参数。

5.根据权利要求4所述的基于地面地震构造约束的VSP数据全波形反演建模方法，其中根据所述VSP数据残差求解逆时传播的伴随波场，以及基于所述伴随波场根据公式(4)求解所述速度模型M₀的梯度g_v(x)：

其中，NS表示震源个数，T代表波场传播的终止时刻，Q(x)是所述逆时传播的伴随波场的分量。

6.根据权利要求3所述的基于地面地震构造约束的VSP数据全波形反演建模方法，其中所述地面地震构造约束滤波算子表述为公式(5)：

其中，D(x)表示扩散张量，且扩散张量D(x)与所述结构特征张量算子的特征值一致，α代表光滑因子。

7.根据权利要求1所述的基于地面地震构造约束的VSP数据全波形反演建模方法，其中采用非精确线性搜索法计算对所述速度模型M₀进行更新的步长。

8.根据权利要求7所述的基于地面地震构造约束的VSP数据全波形反演建模方法，其中采用P-L-BFGS算法求解海塞矩阵的逆矩阵，并基于所述逆矩阵更新所述速度模型M₀。

9.根据权利要求1所述的基于地面地震构造约束的VSP数据全波形反演建模方法，其中所述步骤S2包括以下子步骤S21至S23：

子步骤S21：基于所述时间域数据匹配的全波形反演算子构建VSP数据残差，基于所述VSP数据残差计算所述中间速度模型M₁的梯度；

子步骤S22：通过P-L-BFGS算法更新所述速度模型M₁，获得更新的速度模型M₁’；

子步骤S23：判断是否完成给定次数的迭代，若是则将所述更新的速度模型M₁’作为最终速度模型M并输出；否则返回子步骤S21，以更新的速度模型M₁’作为子步骤S21中的中间速度模型。