基于各向异性全变分约束的弹性波直接包络反演方法

摘要

本发明涉及一种基于各向异性全变分约束的弹性波直接包络反演方法，是通过各向异性全变分约束与弹性波直接包络反演获得弹性强散射介质的高精度纵、横波速度结构。首先，对弹性波场进行波场模式分解分别获得纵、横波场，计算包络场分别获得正传纵、横波包络场；然后，根据梯度表达式，可得到弹性波直接包络反演的纵、横波速度梯度，并可更新纵、横波速度模型；之后，对当前纵、横波速度模型施加各向异性全变分约束可得到约束后的纵、横波速度更新模型；最后，将各向异性全变分约束弹性波直接包络反演结果作为初始模型，进行各向异性全变分约束弹性波全波形反演可得到弹性强散射介质的高精度纵、横波速度模型。

权利要求

1.一种基于各向异性全变分约束的弹性波直接包络反演方法，其特征在于，是在地震数据缺失低频信息且没有模型先验信息的情况下，通过对弹性波直接包络反演过程施加各向异性全变分约束以获得弹性强散射介质的高质量大尺度纵、横波速度结构；将获得的高质量大尺度纵、横波速度结构作为各向异性全变分约束弹性波全波形反演的初始模型，获得弹性强散射介质的高精度纵、横波速度结构；

包括以下步骤：

a、安装MATLAB软件平台；

b、对数据进行静校正、去噪预处理，得到高质量的弹性波多分量观测地震数据；

c、对地震数据进行子波估计，提取每一炮数据的震源子波；

d、通过背景速度分析，得到大致的背景纵、横波速度范围，生成背景速度模型，其中不含有强散射体的任何先验信息，分别作为反演的初始纵波速度模型vp0和初始横波速度模型vs0；

e、在初始纵、横波速度模型上计算弹性波多分量模拟地震数据，对模拟地震数据取包络得到模拟包络数据 上角标i表示i方向的分量，对于二维情况，指x水平和z垂直方向；对观测地震数据取包络得到观测包络数据 通过公式(1)计算弹性波直接包络反演的目标函数σEDEI：式中，求和符号下角标sr表示对所有震源和检波点进行积分，t表示时间，T表示总的记录时间长度；

f、在初始模型上计算模拟地震波场，对模拟地震波场进行波场模式分解得到正传纵波波场和正传横波波场，分别取包络得到纵波正传包络场 和横波正传包络场g、计算模拟包络数据与观测包络数据的差，得到伴随源，伴随源反传得到伴随包络场；

对伴随包络场进行波场模式分解得到纵波伴随包络场 和横波伴随包络场h、纵波正传包络场与纵波伴随包络场进行零延迟互相关得到弹性波直接包络反演的纵波速度梯度，如公式(2)所示：式中，vp表示纵波速度，ρ表示密度，弹性波直接包络反演的横波速度梯度可由公式(3)计算：式中，vs表示横波速度，·表示点乘，μ表示剪切模量；

i、选择合适的步长，采用最速下降法对纵、横波速度模型进行更新，假设当前迭代次数为m，当前迭代更新得到的纵波速度和横波速度模型分别用 和 表示，对二者施加各向异性全变分约束的过程等价于求解公式(4)和(5)所示的最优化问题：其中，J1和J2分别表示对纵、横波速度模型施加各向异性全变分约束的目标函数，α1和α2分别表示纵、横波速度的更新步长， 和 分别为第m次迭代的纵、横波速度的更新量，λ1和λ2分别为对纵、横波速度施加各向异性全变分约束的权系数，||·||表示二范数，||·||ATV表示各向异性全变分范数，对速度模型v计算各向异性全变分范数的具体表达式如公式(6)所示：其中，nz和nx分别表示模型垂向和横向的网格点数，求解公式(4)和(5)所示优化问题，得到各向异性全变分约束后的纵、横波速度模型作为当前迭代更新后的速度模型；

j、在更新后的模型上，进行迭代停止条件判断；如果不满足停止条件，将更新后的纵、横波速度模型作为初始模型，返回第e步骤继续迭代计算；如果满足停止条件，输出结果为各向异性全变分约束后的强散射介质大尺度纵、横波速度结构vpt和vst；

k、以vpt和vst作为初始模型，进行各向异性全变分约束弹性波全波形反演，得到最终反演结果，即强散射介质的高精度纵、横波速度结构。

说明书

基于各向异性全变分约束的弹性波直接包络反演方法

技术领域

[0001] 本发明涉及一种在弹性波直接包络反演中施加各向异性全变分约束以提升强散射介质纵、横波速度结构反演效果的方法。

背景技术

[0002] 强散射介质通常能够形成油气资源的良好封存介质，因此，强散射介质的高精度参数建模方法已经在勘探地球物理领域引起高度重视。全波形反演方法是目前勘探地震领域参数建模精度最高的方法，其能充分利用地震波场的运动学和动力学信息。近年来，人们希望利用全波形反演方法对实际勘探中遇到的强散射介质(如盐丘等)进行高精度参数建模，以提高强散射体及其屏蔽区的成像质量。但是，常规全波形反演方法基于波恩近似这一弱散射近似，且实际地震数据往往缺失有效低频信息，这使得直接利用常规全波形反演方法进行强散射介质的参数建模存在一定难度。近年来，国内外学者从不同角度对全波形反演的方法理论框架进行了改进，以适应强散射介质的参数反演，这些研究绝大多数均是基于声波介质，主要包括Laplace‑Fourier域波形反演方法、标准设置方法、全变分约束方法、深度学习方法、以及直接包络反演方法等。关于弹性强散射介质的多参数建模研究仍处于起步探索阶段，现有方法对于弹性强散射体内部以及下方速度建模效果均不十分理想。总之，当前缺乏一种能够在低频数据缺失且无先验信息情况下对强散射介质进行高精度纵、横波速度建模的方法。

[0003] 直接包络反演方法能够在地震数据缺失低频且无先验信息情况下，对强散射介质进行大尺度参数建模，且其具有相对较高的计算效率。直接包络反演方法与全波形反演方法的最大区别是其定义了直接包络敏感核函数(区别于常规包络反演的波形敏感核函数)，能够将低频包络数据扰动直接映射为强散射介质的大尺度参数扰动。但是，当前的弹性波直接包络反演方法在对弹性强散射介质进行纵、横波速度建模时，难以获得理想的强散射体内部速度信息，从而影响强散射体下边界及其下方的速度建模效果。

发明内容

[0004] 本发明的目的就是针对上述现有技术的不足，提供一种新的强散射介质弹性波多参数反演方法，进而解决油气资源地震勘探中强散射储层纵、横波速度高精度建模问题，为复杂强散射储层成像提供高精度的纵、横波速度模型。

[0005] 本发明的思想是：充分利用各向异性全变分约束能够锐化模型边界和增强层内均匀性的优势，将其引入弹性波直接包络反演过程中。在弹性波直接包络反演每次迭代中，对迭代更新后的纵、横波速度模型分别施加各向异性全变分约束，使得强散射体内部趋于均匀且边界更为明显。最终结合各向异性全变分约束弹性波全波形反演技术获得强散射介质的高精度纵、横波速度模型，克服现有技术的不足。

[0006] 本发明的目的是通过以下技术方案实现的：

[0007] 首先，准备好预处理后的弹性波多分量观测地震数据；利用观测地震数据估计震源子波，给出纵、横波速度的初始模型(无需先验信息)；在初始速度模型上进行正演模拟，得到模拟地震数据与正传地震波场；对正传地震波场进行波场模式分解，分别得到纵、横波正传波场，求取相应的包络场，分别得到纵、横波正传包络场；对观测包络数据和模拟包络数据作差得到包络伴随源；在初始纵、横波速度模型上进行伴随包络场计算与模式分解，得到纵、横波伴随包络场；利用正传包络场与伴随包络场分别计算纵、横波速度梯度；求取步长，并进行纵、横波速度模型的更新；对更新后的纵、横波速度模型分别施加各向异性全变分约束，获得约束后的纵、横波速度模型；进行迭代更新，直到满足停止条件，得到强散射介质的大尺度纵、横波速度结构；以强散射介质的大尺度纵、横波速度模型为初始模型，进行各向异性全变分约束弹性波全波形反演，得到强散射介质的高精度纵、横波速度建模结果。

[0008] 本发明所述的基于各向异性全变分约束的弹性波直接包络反演方法是通过MATLAB平台实现的；

[0009] 本发明所述的基于各向异性全变分约束的弹性波直接包络反演方法，包括以下步骤：

[0010] a、安装MATLAB软件平台；

[0011] b、对数据进行静校正、去噪预处理，得到高质量的弹性波多分量观测地震数据；

[0012] c、对地震数据进行子波估计，提取每一炮数据的震源子波；

[0013] d、通过背景速度分析，得到大致的背景纵、横波速度范围，生成背景速度模型，其中不含有强散射体的任何先验信息，分别作为反演的初始纵波速度模型vp0和初始横波速度模型vs0；

[0014] e、在初始纵、横波速度模型上计算弹性波多分量模拟地震数据，对模拟数据取包络得到模拟包络数据 上角标i表示i方向的分量，对于二维情况，指x(水平)和z(垂直)方向；对观测地震数据取包络得到观测包络数据 通过公式(1)计算弹性波直接包络反演的目标函数σEDEI：

[0015]

[0016] 式中，求和符号下角标sr表示对所有震源和检波点进行积分，t表示时间，T表示总的记录时间长度；

[0017] f、在初始模型上计算模拟地震波场，对模拟地震波场进行波场模式分解得到正传纵波波场和正传横波波场，分别取包络得到纵波正传包络场 和横波正传包络场[0018] g、计算模拟包络数据与观测包络数据的差，得到伴随源，伴随源反传得到伴随包络场；对伴随包络场进行波场模式分解得到纵波伴随包络场 和横波伴随包络场[0019] h、纵波正传包络场与纵波伴随包络场进行零延迟互相关得到弹性波直接包络反演的纵波速度梯度，如公式(2)所示：

[0020]

[0021] 式中，vp表示纵波速度，ρ表示密度。弹性波直接包络反演的横波速度梯度可由公式(3)计算：

[0022]

[0023] 式中，vs表示横波速度，·表示点乘，μ表示剪切模量；

[0024] i、选择合适的步长，采用最速下降法对纵、横波速度模型进行更新。假设当前迭代次数为m，当前迭代更新得到的纵波速度和横波速度模型分别用 和 表示，对二者施加各向异性全变分约束的过程等价于求解公式(4)和(5)所示的最优化问题：

[0025]

[0026]

[0027] 其中，J1和J2分别表示对纵、横波速度模型施加各向异性全变分约束的目标函数，α1和α2分别表示纵、横波速度的更新步长， 和 分别为第m次迭代的纵、横波速度的更新量，λ1和λ2分别为对纵、横波速度施加各向异性全变分约束的权系数，||·||表示二范数，||·||ATV表示各向异性全变分范数。对速度模型v计算各向异性全变分范数的具体表达式如公式(6)所示：

[0028]

[0029] 其中，nz和nx分别表示模型垂向和横向的网格点数。求解公式(4)和(5)所示优化问题，得到各向异性全变分约束后的纵、横波速度模型作为当前迭代更新后的速度模型；

[0030] j、在更新后的模型上，进行迭代停止条件判断；如果不满足停止条件，将更新后的纵、横波速度模型作为初始模型，返回第e步骤继续迭代计算；如果满足停止条件，输出结果为各向异性全变分约束后的强散射介质大尺度纵、横波速度结构vpt和vst；

[0031] k、以vpt和vst作为初始模型，进行各向异性全变分约束弹性波全波形反演，得到最终反演结果，即强散射介质的高精度纵、横波速度结构。

[0032] 与现有技术相比，本发明的有益效果在于：本发明将各向异性全变分约束引入弹性波直接包络反演过程中，能够促使强散射体内部速度均匀、边界信息突出，从而提升强散射体边界刻画以及下方屏蔽区的速度建模效果，进而获得弹性强散射介质的高精度纵、横波速度反演结果。

[0033] 有以下优点：1.本发明通过引入各向异性全变分约束，使得弹性波直接包络反演过程中更易获得高质量的强散射体速度信息。2.本发明涉及的各向异性全变分约束弹性波直接包络反演能够更好的刻画强散射体边界信息，尤其是强散射体下边界及其下方屏蔽区速度信息。3.将基于各向异性全变分约束的弹性波直接包络反演与弹性波全波形反演方法串联，能够兼顾介质大尺度宏观构造与小尺度细节构造的反演，并最终获得强散射介质的高精度纵、横波速度结构。4.本发明可以在地震数据低频信息缺失且无模型先验信息情况下进行弹性强散射介质的高精度多参数建模。

附图说明

[0034] 图1基于各向异性全变分约束的弹性波直接包络反演方法的流程图；

[0035] 图2真实速度模型图；

[0036] (a)真实纵波速度模型图，(b)真实横波速度模型图；

[0037] 图3震源子波及其频谱图；

[0038] (a)震源子波图(b)震源子波频谱图；

[0039] 图4初始速度模型图；

[0040] (a)初始纵波速度模型图，(b)初始横波速度模型图；

[0041] 图5常规弹性波全波形反演结果；

[0042] (a)常规弹性波全波形反演纵波速度结果；(b)常规弹性波全波形反演横波速度结果；

[0043] 图6无约束弹性波直接包络反演结果；

[0044] (a)无约束弹性波直接包络反演纵波速度结果；(b)无约束弹性波直接包络反演横波速度结果；(c)纵波速度最终反演结果；(d)横波速度最终反演结果；

[0045] 图7基于各向异性全变分约束的弹性波直接包络反演结果；

[0046] (a)基于各向异性全变分约束的弹性波直接包络反演纵波速度结果；(b)基于各向异性全变分约束的弹性波直接包络反演横波速度结果；(c)纵波速度最终反演结果；(d)横波速度最终反演结果。

具体实施方式

[0047] 下面结合附图和实例对本发明进一步的详细说明。

[0048] 本发明所述的基于各向异性全变分约束的弹性波直接包络反演，包括以下步骤：

[0049] a、在win7或Linux系统下安装MATLAB软件平台，要求采用MATLAB R2016a及以上版本，并且已配备并行工具包(Parallel Computing Toolbox)。

[0050] b、进行数据预处理，对数据进行静校正处理，校正起伏地表对反射同相轴的影响；对数据进行去噪处理，去除微震、低频和高频背景噪声及其他随机噪声；去除干扰波，包括声波、面波、工业电干扰、虚反射、多次反射、侧面波、底波、交混回响和鸣震等。最终高质量的弹性波多分量观测地震数据。

[0051] c、对地震数据进行子波估计，估计方法可以采用直达波估计法和自相关法等，提取出每一炮数据的震源子波。

[0052] d、通过背景速度分析，得到大致的背景纵、横波速度范围，生成背景速度模型，其中不含有强散射体的任何先验信息，分别作为反演的初始纵波速度模型vp0和初始横波速度模型vs0。

[0053] e、在初始纵、横波速度模型上计算弹性波多分量模拟地震数据，对模拟数据取包络得到模拟包络数据 上角标i表示i方向的分量，对于二维情况，指x(水平)和z(垂直)方向；对观测地震数据取包络得到观测包络数据 通过公式(7)计算弹性波直接包络反演的目标函数σEDEI：

[0054]

[0055] 式中，求和符号下角标sr表示对所有震源和检波点进行积分，t表示时间，T表示总的记录时间长度。

[0056] f、在初始模型上计算模拟地震波场，对模拟地震波场进行波场模式分解得到正传纵波波场和正传横波波场，取包络得到纵波正传包络场 和横波正传包络场 其中，波场模式分解采用公式(8)和(9)的方法：

[0057]

[0058]

[0059] 其中，u表示弹性波位移矢量，up表示纵波波场，us表示横波波场，ρ表示密度，λ和μ表示拉梅常数，▽表示梯度运算，▽·表示散度运算，▽×表示旋度运算。

[0060] g、计算模拟包络数据与观测包络数据的差，得到伴随源fs，即

[0061]

[0062] 伴随源反传得到伴随包络场；对伴随包络场进行波场模式分解得到纵波伴随包络场 和横波伴随包络场 波场模式分解方法见公式(8)和(9)。

[0063] h、纵波正传包络场与纵波伴随包络场进行零延迟互相关得到弹性波直接包络反演的纵波速度梯度，如公式(11)所示：

[0064]

[0065] 式中，vp表示纵波速度，ρ表示密度。弹性波直接包络反演的横波速度梯度可由公式(12)计算：

[0066]

[0067] 式中，vs表示横波速度，·表示点乘，μ表示剪切模量。

[0068] i、选择合适的步长，采用最速下降法对纵、横波速度模型进行更新。假设当前迭代次数为m，当前迭代更新得到的纵波速度和横波速度模型分别用 和 表示，对二者施加各向异性全变分约束的过程等价于求解公式(13)和(14)所示的最优化问题：

[0069]

[0070]

[0071] 其中，J1和J2分别表示对纵、横波速度模型施加各向异性全变分约束的目标函数，α1和α2分别表示纵、横波速度的更新步长， 和 分别为第m次迭代的纵、横波速度的更新量，λ1和λ2分别为对纵、横波速度施加各向异性全变分约束的权系数，||·||表示二范数，||·||ATV表示各向异性全变分范数。对速度模型v计算各向异性全变分范数的具体表达式如公式(15)所示：

[0072]

[0073] 其中，nz和nx分别表示模型垂向和横向的网格点数。求解公式(13)和(14)所示优化问题，得到各向异性全变分约束后的纵、横波速度模型作为当前迭代更新后的速度模型。

[0074] j、在更新后的模型上，进行迭代停止条件判断；如果不满足停止条件，将更新后的纵、横波速度模型作为初始模型，返回第e步骤继续迭代计算；如果满足停止条件，输出结果为各向异性全变分约束后的强散射介质大尺度纵、横波速度结构vpt和vst。

[0075] k、以vpt和vst作为初始模型，进行各向异性全变分约束弹性波全波形反演，得到最终反演结果，即强散射介质的高精度纵、横波速度结构。

[0076] 实施例1

[0077] 本发明的整体流程如图1所示。

[0078] 假设地下真实纵、横波速度模型分别如图2a和2b所示。在真实速度模型中，背景速度较低，中部有一高速强散射盐丘体。在真实模型上进行正演模拟可以得到观测地震记录，震源子波波形和频谱分别如图3a和3b所示，为模拟实际地震采集中低频信息缺失的情况，对雷克子波进行高通滤波处理，切去了3Hz以下的低频信息，震源子波主频约为9Hz。利用观测记录进行背景速度分析，得到背景纵波速度和横波速度大致范围，建立纵、横波初始速度模型如图4a和图4b所示。初始纵、横波速度模型中均不含强散射体的先验信息。

[0079] 为对比本发明方法的反演效果，首先在初始纵、横波速度模型上进行常规弹性波全波形反演，纵、横波速度反演结果分别如图5a和5b所示。可见，由于缺失低频信息，常规弹性波全波形反演只能得到强散射体顶界面的部分信息，无法恢复强散射体内部的速度信息与强散射体的形态信息。

[0080] 然后，进行无约束的弹性波直接包络反演方法，纵、横波速度反演结果分别如图6a和6b所示。可见，虽然强散射体形态得到一定程度的恢复，但是盐丘体内部速度较真实值仍存在偏差。无约束情况下，以图6a和图6b的结果为初始模型，进行常规弹性波全波形反演，最终纵、横波速度反演结果如图6c和图6d所示。可见，最终结果盐丘内部速度均匀性不好，且盐丘体下边界刻画不充分。

[0081] 然后，利用本发明方法进行反演。以图4所示模型为初始模型，根据步骤(e)到步骤(j)，可以得到的纵、横波速度模型如图7a和图7b所示，可见，强散射体的大尺度纵、横波速度构造信息已经被成功恢复，且盐丘内速度信息较为接近真实值。以图7a和7b为初始模型，执行步骤(k)，即进行各向异性全变分约束弹性波全波形反演，最终纵、横波速度反演结果如图7c和7d所示。可见，强散射体的边界信息与内部速度被较好地反演出来。盐丘体内部、下边界和盐下区域反演效果均较无约束情况下好。总之，本发明提出方法的最终反演结果盐丘体内部速度与盐丘体边界均较为准确，其总体效果明显好于常规方法(图5和图6)。