CN113138416B - 一种敏感核函数优化的全波形反演速度建模方法 - Google Patents

Abstract

本发明公开了一种敏感核函数优化的全波形反演速度建模方法，涉及勘探地球物理领域。反演方法在迭代早期增大正向散射角信息，用于恢复模型低波数部分，随着迭代次数的增大逐渐增强反向散射角信息，用于构建模型中的中高波数部分。该反演流程既能够自然减弱反演的非线性，避免陷入局部极值，又可以降低多尺度反演的计算成本，具有较好的应用前景。

Description

一种敏感核函数优化的全波形反演速度建模方法

技术领域

本发明涉及勘探地球物理领域，具体涉及一种敏感核函数优化的全波形反演速度建模方法。

背景技术

速度建模是地震数据处理中至关重要的一步，建模的精度直接决定了后续偏移成像处理的质量。常规速度建模方法包括偏移速度分析和射线层析反演，由于这些方法使用地震波旅行时信息，所估计的速度模型主要是低波数部分，而中、高波数部分缺失，因此很难获得全波数段高分辨率速度模型。为了解决这一难题，全波形反演应运而生，其使用地震全波场信息(旅行时、振幅、多次波、散射波等)，利用合成数据拟合观测数据来估计地下岩石属性参数，属于一种非线性参数反演方法。该方法能够获得低、中、高全波数速度信息，反演结果可对地下复杂构造形态进行高精度地刻画和描述。

全波形反演属于一种非线性反演方法，可使用全局或局部优化算法进行求解。由于全局优化算法计算成本太高，通常采用局部梯度优化算法进行求解。但在实际应用中，受不规则数据采集、有限偏移距、环境噪声干扰以及低频信号缺失等因素影响，全波形反演的非线性变得极为复杂，常规梯度算法很容易发生周期跳跃问题，进而陷入局部极小值，获得错误的反演结果。多尺度反演策略是解决局部极值的一个常用方法，即从低频到高频依次进行反演，逐步恢复低波数、中波数和高波数速度信息。成功实现多尺度全波形反演的一个基本要求是野外观测数据含有低频信号(通常要求最低频率为2～3Hz)，但常规检波器的最低频率往往大于等于5Hz，因此对于常规检波器采集的地震数据，直接应用全波形反演仍具有较大挑战性，急需研发一种适用于常规数据采集的稳定全波形反演方法。

发明内容

本发明提出了一种敏感核函数优化的全波形反演速度建模方法，在迭代早期增大正向散射角信息，用于恢复模型低波数部分，随着迭代次数的增大逐渐增强反向散射角信息，用于构建模型中的中高波数部分。

本发明具体采用如下技术方案：

一种敏感核函数优化的全波形反演速度建模方法，包括以下步骤：

(1)获取输入数据，输入数据包括初始纵波速度模型v₀(x)、初始密度模型ρ₀(x)、震源函数子波f(t)、野外观测数据d_obs(x_r,t)、预设最大迭代次数N_iter和数据冗余误差ε，其中x_r为检波点位置；

(2)使用纵波速度模型v_i(x)和初始纵波密度模型ρ₀(x)，计算纵波阻抗模型z_i(x)＝ρ₀v_i(x)，其中i＝[0,1,2...N_iter]为迭代次数，当i＝0，v_i(x)使用输入的初始速度模型v₀(x)；

(3)使用震源子波f(t)和阻抗模型，通过求解声波波动方程，计算正向延拓的压强波场p(x,t)和合成数据d_syn(x_r,t)＝p(x_r,t)δ(x-x_r)，声波波动方程为式(1)所示：

x为地下模型空间的位置坐标，x_s为震源位置，

为空间梯度算子，

为散度算子，δ(x)为狄拉克δ函数，

为二阶时间导数；

(4)使用合成数据d_syn(x_r,t)和野外观测数据d_obs(x_r,t)，计算数据残差，计算式为式(2)所示：

r(x_r,t)＝d_syn(x_r,t)-d_obs(x_r,t) (2)

(5)使用数据残差r(x_r,t)、纵波速度模型v_i(x)和纵波阻抗模型z_i(x)，通过求解共轭声波波动方程计算反向延拓共轭波场，共轭声波波动方程的表达式为(3)所示：

其中，

为共轭压强波场，T为地震记录长度；

(6)使用正向延拓的压强波场p(x,t)和反向延拓的共轭压强波场

计算速度和阻抗敏感核函数，表达式如(4)所示：

其中，K_v(x)为速度敏感核函数，K_z(x)为阻抗敏感核函数；

(7)使用速度和阻抗敏感核函数计算优化的模型梯度，计算式为(5)所示：

g＝λK_v(x)+K_z(x) (5)

其中λ为一个可调节的标量参数，用于控制速度和阻抗核函数对梯度的贡献大小，在λ设置为式(6)所示

其中i＝[0,1,2...N_iter]为迭代次数；

(8)使用基于敏感核函数优化的模型梯度g和抛物线拟合方法计算模型更新步长α，更新模型为(7)所示：

v_i+1(x)＝v_i(x)+αg(x) (7)；

(9)重复步骤(2)～(8)进行迭代，直到达到最大迭代次数i＝N_iter，或者数据残差逼近预设的冗余误差，数据残差计算公式为式(8)所示：

本发明具有如下有益效果：

本申请使用敏感核函数优化的模型梯度进行更新，主要解决全波形反演建模的周期跳跃问题，避免反演算法陷入局部极值，进而收敛至准确解。

相对常规全波形反演方法，本申请记载的方法能够在早期迭代过程中使用速度敏感核占优的低波数信息，在后续迭代过程中逐渐增加阻抗敏感核信息，恢复中高波数部分，最终获得全波数段模型信息。

本申请记载的建模方式能够降低全波形反演的非线性，进而避免周期跳跃问题。相对常规多尺度全波形反演算法，本申请记载的方法不需要在不同频段进行反演，总的迭代次数相对较少，具有较高的计算效率。

附图说明

图1为Marmousi模型真实速度场；

图2为全波形反演需要的Marmousi模型初始速度场；

图3为Marmousi模型在第一次迭代中的常规梯度；

图4为Marmousi模型在第一次迭代中的由速度敏感核计算的梯度；

图5为Marmousi模型在第一次迭代中的由阻抗敏感核计算的梯度；

图6为使用常规全波形反演获得Marmousi速度模型；

图7为使用常规多尺度全波形反演获得Marmousi速度模型；

图8为使用敏感核函数优化的全波形反演速度建模方法的基于优化敏感核函数的全波形反演获得Marmousi速度模型；

图9为三种全波形反演方法数据收敛曲线对比结果示意图。

具体实施方式

下面结合附图和具体实施例对本发明的具体实施方式做进一步说明：

常规全波形反演使用纵波速度作为反演参数，局部寻优算法中计算的速度梯度包含所有散射角信息，笼统地使用这些信息进行速度模型更新，容易使得反演过程陷入局部极值采用拉格朗日乘数法推导了速度和阻抗敏感核函数，通过分析发现速度敏感核函数主要为正向散射角信息，主要反映了地下参数模型的低波数部分，而阻抗敏感核函数主要为反向散射角信息，主要反映了参数模型的中高波数部分。

基于上述分析，提出了敏感核函数优化的全波形反演速度建模方法，在迭代早期增大正向散射角信息，用于恢复模型低波数部分，随着迭代次数的增大逐渐增强反向散射角信息，用于构建模型中的中高波数部分，该反演流程既能够自然减弱反演的非线性，避免陷入局部极值，敏感核函数优化的全波形反演速度建模方法包括以下步骤：

(1)获取输入数据，输入数据包括初始纵波速度模型v₀(x)、初始密度模型ρ₀(x)、震源函数子波f(t)、野外观测数据d_obs(x_r,t)、预设最大迭代次数N_iter和数据冗余误差ε，其中x_r为检波点位置；

(2)使用纵波速度模型v_i(x)和初始纵波密度模型ρ₀(x)，计算纵波阻抗模型z_i(x)＝ρ₀v_i(x)，其中i＝[0,1,2...N_iter]为迭代次数，当i＝0，v_i(x)使用输入的初始速度模型v₀(x)；

(3)使用震源子波f(t)和阻抗模型，通过求解声波波动方程，计算正向延拓的压强波场p(x,t)和合成数据d_syn(x_r,t)＝p(x_r,t)δ(x-x_r)，声波波动方程为式(1)所示：

x为地下模型空间的位置坐标，x_s为震源位置，

为空间梯度算子，

为散度算子，δ(x)为狄拉克δ函数，

为二阶时间导数；

(4)使用合成数据d_syn(x_r,t)和野外观测数据d_obs(x_r,t)，计算数据残差，计算式为式(2)所示：

r(x_r,t)＝d_syn(x_r,t)-d_obs(x_r,t) (2)

(5)使用数据残差r(x_r,t)、纵波速度模型v_i(x)和纵波阻抗模型z_i(x)，通过求解共轭声波波动方程计算反向延拓共轭波场，共轭声波波动方程的表达式为(3)所示：

其中，

为共轭压强波场，T为地震记录长度；

(6)使用正向延拓的压强波场p(x,t)和反向延拓的共轭压强波场

计算速度和阻抗敏感核函数，表达式如(4)所示：

其中，K_v(x)为速度敏感核函数，K_z(x)为阻抗敏感核函数；

(7)使用速度和阻抗敏感核函数计算优化的模型梯度，计算式为(5)所示：

g＝λK_v(x)+K_z(x) (5)

其中λ为一个可调节的标量参数，用于控制速度和阻抗核函数对梯度的贡献大小，在λ设置为式(6)所示

其中i＝[0,1,2...N_iter]为迭代次数；

(8)使用基于敏感核函数优化的模型梯度g和抛物线拟合方法计算模型更新步长α，更新模型为(7)所示：

v_i+1(x)＝v_i(x)+αg(x) (7)；

(9)重复步骤(2)～(8)进行迭代，直到达到最大迭代次数i＝N_iter，或者数据残差逼近预设的冗余误差，数据残差计算公式为式(8)所示：

将本申请记载的建模方法应用于国际标准Marmousi模型速度反演建模中，取得了较理想的反演结果。真实Marmousi速度模型如图1所示，输入的初始模型如图2所示，震源函数是主频为8Hz的雷克子波，观测数据集为38个炮集，均匀分布在地表，间距为250m，每个炮集包含761个检波点，间隔为12.5m。常规全波形反演第一次迭代的模型梯度如图3所示，该模型梯度正向和反向散射角信息混叠在一起，使用申请提出的速度和阻抗敏感核函数计算的梯度如图4和图5所示，可以看出速度梯度主要包含低波数信息，而阻抗敏感核函数主要包含中高波数信息，通过组合这两个梯度可以获得一个最优的模型梯度用于避免周波跳跃问题。常规全波形反演的结果如图6所示，常规多尺度全波形反演的结果如图7所示，本申请提出的全波形反演方法的结果如图8所示。可以看出，由于常规全波形反演直接使用低波数和高波数混叠的模型梯度，导致反演过程陷入了局部极值，没有正确恢复地下速度模型。常规多尺度反演虽然可以较好恢复速度模型，但是由于需要分多个频段进行波形反演，计算量较大。而本发明专利提出的新型全波形反演方法不仅能够获得准确的速度模型，而且具有和常规方法一样的计算时间，具有较好的实际应用价值。图1-8中Depth为模型深度，Distance为模型横向距离。图9中数据点用实心圆表示的实线为常规多尺度全波形反演方法(MCSV-based FWI)，数据点用三角形表示的实线为常规全波形反演方法(CSV-basedFWI)，虚线为本申请提出的基于敏感核函数优化的混合梯度全波形反演方法(HG-basedFWI)。

当然，上述说明并非是对本发明的限制，本发明也并不仅限于上述举例，本技术领域的技术人员在本发明的实质范围内所做出的变化、改型、添加或替换，也应属于本发明的保护范围。

Claims (1)

1.一种敏感核函数优化的全波形反演速度建模方法，其特征在于，包括以下步骤：

(1)获取输入数据，输入数据包括初始纵波速度模型v₀(x)、初始密度模型ρ₀(x)、震源函数子波f(t)、野外观测数据d_obs(x_r,t)、预设最大迭代次数N_iter和数据冗余误差ε，其中x_r为检波点位置；

(2)使用纵波速度模型v_i(x)和初始密度模型ρ₀(x)，计算纵波阻抗模型z_i(x)＝ρ₀v_i(x)，其中i＝[0,1,2...N_iter]为迭代次数，当i＝0，v_i(x)使用输入的初始纵波速度模型v₀(x)；

(3)使用震源函数子波f(t)和阻抗模型，通过求解声波波动方程，计算正向延拓的压强波场p(x,t)和合成数据d_syn(x_r,t)＝p(x_r,t)δ(x-x_r)，声波波动方程为式(1)所示：

x为地下模型空间的位置坐标，x_s为震源位置，

为空间梯度算子，

为散度算子，δ(x)为狄拉克δ函数，

为二阶时间导数；

(4)使用合成数据d_syn(x_r,t)和野外观测数据d_obs(x_r,t)，计算数据残差，计算式为式(2)所示：

r(x_r,t)＝d_syn(x_r,t)-d_obs(x_r,t) (2)

(5)使用数据残差r(x_r,t)、纵波速度模型v_i(x)和纵波阻抗模型z_i(x)，通过求解共轭声波波动方程计算反向延拓共轭压强波场，共轭声波波动方程的表达式为(3)所示：

其中，

为共轭压强波场，T为地震记录长度；

(6)使用正向延拓的压强波场p(x,t)和反向延拓的共轭压强波场

计算速度和阻抗敏感核函数，表达式如(4)所示：

其中，K_v(x)为速度敏感核函数，K_z(x)为阻抗敏感核函数；

(7)使用速度和阻抗敏感核函数计算优化的模型梯度，计算式为(5)所示：

g＝λK_v(x)+K_z(x) (5)

其中λ为一个可调节的标量参数，用于控制速度和阻抗敏感核函数对梯度的贡献大小，λ设置为式(6)所示

其中i＝[0,1,2...N_iter]为迭代次数；

(8)使用基于敏感核函数优化的模型梯度g和抛物线拟合方法计算模型更新步长α，更新模型为(7)所示：

v_i+1(x)＝v_i(x)+αg(x) (7)；

(9)重复步骤(2)～(8)进行迭代，直到达到最大迭代次数i＝N_iter，或者数据残差逼近预设的冗余误差，数据残差计算公式为式(8)所示：

Images