CN106033124A - 一种基于随机最优化的多震源粘声最小二乘逆时偏移方法 - Google Patents

Abstract

本发明提供一种基于随机最优化的多震源粘声最小二乘逆时偏移方法，该基于随机最优化的多震源粘声最小二乘逆时偏移方法包括：1）读取野外观测记录及预设参数合成超道集；2）采用当前反射系数模型，通过多震源激发正演模拟超道集，计算数据残差；3）根据数据残差计算更新梯度；4）通过随机最优化思想修改梯度并计算更新步长；5）由梯度及更新步长更新反射系数模型。本发明提供的方法提出了将随机最优化思想推广到相位编码的粘声最小二乘逆时偏移中，通过加权平均之前的梯度，减小了梯度的随机波动，得到了较好的效果。

Description

一种基于随机最优化的多震源粘声最小二乘逆时偏移方法

技术领域

本发明涉及油气物探工程领域，特别是涉及到一种基于随机最优化的多震源粘声最小二乘逆时偏移方法。

背景技术

实际地下介质的粘滞性是普遍存在的，地震波在黏弹性介质中的传播主要表现为速度频散与振幅衰减。不考虑黏滞性的叠前成像算法不仅会使成像位置发生偏离，而且还会引起成像振幅的欠估计，严重影响甚至误导随后的地震数据处理、解释等工作。随着油气勘探开发的深入，勘探精度要求逐渐提高，地震波成像也逐步从构造成像向岩性成像发展。然而，目前常规的偏移成像方法还无法满足岩性油气藏勘探开发的需求，究其原因是由于常规偏移算子是正演算子的共轭转置，而不是其逆算子。而基于反演思想的最小二乘偏移相对于常规偏移来说，具有更高的成像分辨率、振幅保真性和均衡性以及压制偏移噪音等优势，越来越受到学者的重视，然而计算量过于庞大限制了其进一步推广应用。

由于最小二乘逆时偏移(LSRTM)的计算量与炮数成线性关系，因而通过相位编码技术将多个炮集组合成一个超道集，可有效减小计算量。然而研究发现相位编码算法的目标泛函是真实目标泛函的随机无偏估计，其梯度也是如此。由于相位编码LSRTM的梯度是随机的，因而其步长也应该是随机的。然而，目前相位编码LSRTM求解时仍然采用与传统LSRTM算法相同的确定性最优化解法，例如最速下降法、共轭梯度法等，忽略了梯度和步长的随机性。为此我们发明了一种新的基于随机最优化的多震源粘声最小二乘逆时偏移方法，解决了以上技术问题。

发明内容

本发明的目的是提供一种考虑实际地下粘滞性的高效率和高精度的基于随机最优化的多震源粘声最小二乘逆时偏移方法。

本发明的目的可通过如下技术措施来实现：基于随机最优化的多震源粘声最小二乘逆时偏移方法，该基于随机最优化的多震源粘声最小二乘逆时偏移方法包括：1)读取野外观测记录及预设参数合成超道集；2)采用当前反射系数模型，通过多震源激发正演模拟超道集，计算数据残差；3)根据数据残差计算更新梯度；4)通过随机最优化思想修改梯度并计算更新步长；5)由梯度及更新步长更新反射系数模型。

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

在步骤1，输入初始反射系数模型、偏移速度场、观测数据、品质因子、迭代终止的阈值及偏移参数，初始反射系数模型的值为0，即第1次迭代与常规粘声逆时偏移等价。

在步骤2，计算数据残差时，基于标准线性固体模型的扰动波场的控制方程为：

$\frac{1}{v_0^2}\frac{{\∂}^2{p}_s}{\∂t^2}=\frac{\∂}{\∂t}\[\mathrm{\ρ}\left(1-\underset{i=1}{\overset{I}{\Σ}}\left(1-\frac{{\mathrm{\τ}}_{\mathrm{\ϵ}i}}{{\mathrm{\τ}}_{\mathrm{\σ}i}}\right)e^{-t/{\mathrm{\τ}}_{\mathrm{\σ}i}}\right)H\left(t\right)\]*\▿\·(\frac{1}{\mathrm{\ρ}}\▿p_s)-m\left(x\right)\frac{{\∂}^2{p}_0}{\∂t^2}---\left(1\right)$

其中，p_s为扰动波场，v₀为背景速度，ρ为密度，I为标准线性固体的个数，τ_εi,τ_σi为松弛时间；H(t)为单位阶跃函数，▽为梯度算子，▽·为散度算子，*为时间上的卷积算子，m(x)为模型参数，即反射系数模型，p₀为背景波场，即在背景介质中传播的波，其控制方程为：

$\frac{1}{v_0^2}\frac{{\∂}^2{p}_0}{\∂t^2}=\frac{\∂}{\∂t}\[\mathrm{\ρ}\left(1-\underset{i=1}{\overset{I}{\Σ}}\left(1-\frac{{\mathrm{\τ}}_{\mathrm{\ϵ}i}}{{\mathrm{\τ}}_{\mathrm{\σ}i}}\right)e^{-t/{\mathrm{\τ}}_{\mathrm{\σ}i}}\right)H\left(t\right)\]*\▿\·(\frac{1}{\mathrm{\ρ}}\▿p_0)+f---\left(2\right)$

其中，f为经过编码的震源项，

公式(1)可用矩阵算子形式表示为:

p_s＝Lm (3)

其中，L为粘声介质扰动波场的线性正演算子。

在步骤2，在计算扰动波场和背景波场时，松弛时间τ_σ、τ_ε的计算公式如公式(4)所示：

$\left\{\begin{array}{c}{\mathrm{\τ}}_{\mathrm{\σ}}=\frac{\sqrt{Q^2+1}-1}{wQ}\\ {\mathrm{\τ}}_{\mathrm{\ϵ}}=\frac{\sqrt{Q^2+1}+1}{wQ}\end{array}\right.---\left(4\right)$

其中，w为圆频率，Q为品质因子。

在步骤3，通过数据残差计算更新模型的梯度方向，计算公式如公式(5)所示:

g＝L^*(Lm-p_obs) (5)

其中，g为梯度，p_obs为观测记录，L^*为正演算子的共轭转置，即逆时偏移算子，L为扰动波场的线性正演算子，m为反射系数模型。

该基于随机最优化的多震源粘声最小二乘逆时偏移方法还包括，在步骤3之后，判断梯度是否满足迭代终止条件，即梯度的模小于预设的阈值，若满足，输出当前反射系数模型，流程结束；否则流程进入到步骤4。

在步骤4，将随机最优化思想推广到相位编码的最小二乘逆时偏移算法中修改梯度，随机最优化方法需要加权平均之前的梯度，因此在前几次迭代时不需要修改梯度，修改后的梯度如公式(6)所示，

$\stackrel{\‾}{g^k}=\frac{{\mathrm{\Σ}}_{i=k-j}^{i=k}{e}^{a(i-k)}{g}^k}{{\mathrm{\Σ}}_{i=k-j}^{i=k}{e}^{a(i-k)}}---\left(6\right)$

其中，为修改后的梯度，g为梯度，k为当前迭代次数；j为加权的前期迭代次数，综合考虑效果和效率，令其等于10；e为自然常数，a为衰减因子，取为0.4。

在步骤4，由修改后的梯度计算更新步长，如公式(7)所示，

${\mathrm{\α}}^k=\frac{{\left(g^k\right)}^{*}{g}^k}{{\left({\mathrm{Lg}}^k\right)}^{*}{\mathrm{Lg}}^k}---\left(7\right)$

其中，α^k为第k次迭代的更新步长，g^k为第k次迭代的梯度，L为扰动波场的线性正演算子。

在步骤5，由公式(8)更新反射系数模型，

m^k+1＝m^k-α^kP^kg^k (8)

其中，P^k为预处理算子，m^k为第k次迭代的反射系数模型，α^k为第k次迭代的更新步长，g^k为第k次迭代的梯度，采用背景波场的能量来近似Hessian矩阵的对角元素，在减少计算量的同时加速了收敛速度，更新反射系数模型后流程返回到步骤2。

本发明中的基于随机最优化的多震源粘声最小二乘逆时偏移方法，是一种有效提高计算效率、快速压制串扰噪声的，且补偿粘滞性的真振幅最小二乘逆时偏移方法。该方法将实际地下的粘滞性引入到反演成像的最小二乘逆时偏移中，不仅能很好的避免常规粘声介质成像的不稳定性，且有效补偿了由粘滞性造成的能量的吸收衰减。考虑到基于常规相位编码的多震源最小二乘逆时偏移方法压制串扰较慢，将随机最优化思想推广到多震源最小二乘逆时偏移中，有效提高了计算效率和收敛速度。

附图说明

图1为本发明的基于随机最优化的多震源粘声最小二乘逆时偏移方法的一具体实施例的流程图；

图2为本发明的一具体实施例中Marmousi模型速度场、扰动模型及品质因子模型的示意图；

图3为本发明的一具体实施例中基于随机最优化的多震源粘声LSRTM反演成像结果的示意图；

图4为本发明的一具体实施例中数据残差曲线的示意图。

具体实施方式

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

如图1所示，图1为本发明的基于随机最优化的多震源粘声最小二乘逆时偏移方法的流程图。

在步骤101，读取野外观测记录及预设参数。输入初始反射系数模型、偏移速度场、观测数据、品质因子、迭代终止的阈值及偏移参数。初始模型的值为0，即第1次迭代与常规粘声逆时偏移等价。流程进入到步骤102。

在步骤102，通过相位编码合成超道集。流程进入到步骤103。

在步骤103，采用当前反射系数模型，通过多震源激发正演模拟扰动波场，并计算数据残差，基于标准线性固体(GSLS)模型的扰动波场的控制方程为：

$\frac{1}{v_0^2}\frac{{\∂}^2{p}_s}{\∂t^2}=\frac{\∂}{\∂t}\[\mathrm{\ρ}\left(1-\underset{i=1}{\overset{I}{\Σ}}\left(1-\frac{{\mathrm{\τ}}_{\mathrm{\ϵ}i}}{{\mathrm{\τ}}_{\mathrm{\σ}i}}\right)e^{-t/{\mathrm{\τ}}_{\mathrm{\σ}i}}\right)H\left(t\right)\]*\▿\·(\frac{1}{\mathrm{\ρ}}\▿p_s)-m\left(x\right)\frac{{\∂}^2{p}_0}{\∂t^2}---\left(1\right)$

其中，p_s为扰动波场，v₀为背景速度，ρ为密度，I为标准线性固体的个数，τ_εi,τ_σi为松弛时间；H(t)为单位阶跃函数，▽为梯度算子，▽·为散度算子，*为时间上的卷积算子，m(x)为模型参数，即反射系数模型。p₀为背景波场，即在背景介质中传播的波，其控制方程为：

$\frac{1}{v_0^2}\frac{{\∂}^2{p}_0}{\∂t^2}=\frac{\∂}{\∂t}\[\mathrm{\ρ}\left(1-\underset{i=1}{\overset{I}{\Σ}}\left(1-\frac{{\mathrm{\τ}}_{\mathrm{\ϵ}i}}{{\mathrm{\τ}}_{\mathrm{\σ}i}}\right)e^{-t/{\mathrm{\τ}}_{\mathrm{\σ}i}}\right)H\left(t\right)\]*\▿\·(\frac{1}{\mathrm{\ρ}}\▿p_0)+f---\left(2\right)$

其中，f为经过编码的震源项。

公式(1)可用矩阵算子形式表示为:

p_s＝Lm (3)

在计算扰动波场和背景波场时，松弛时间的计算公式如公式(4)所示：

$\left\{\begin{array}{c}{\mathrm{\τ}}_{\mathrm{\σ}}=\frac{\sqrt{Q^2+1}-1}{wQ}\\ {\mathrm{\τ}}_{\mathrm{\ϵ}}=\frac{\sqrt{Q^2+1}+1}{wQ}\end{array}\right.---\left(4\right)$

其中，w为圆频率，Q为品质因子。

流程进入到步骤104。

在步骤104，通过数据残差计算更新模型的梯度方向，计算公式如公式(5)所示。

g＝L^*(Lm-p_obs) (5)

其中，g为梯度，p_obs为观测记录，L^*为正演算子的共轭转置，即逆时偏移算子。流程进入到步骤105。

在步骤105，判断梯度是否满足迭代终止条件，即梯度的模小于预设的阈值，若满足，输出当前反射系数模型，流程进入到步骤109；否则流程进入到步骤106。

步骤106，将随机最优化思想推广到相位编码的最小二乘逆时偏移算法中修改梯度，随机最优化方法需要加权平均之前的梯度，因此在前几次迭代时不需要修改梯度。修改后的梯度如公式(6)所示，

$\stackrel{\‾}{g^k}=\frac{{\mathrm{\Σ}}_{i=k-j}^{i=k}{e}^{a(i-k)}{g}^k}{{\mathrm{\Σ}}_{i=k-j}^{i=k}{e}^{a(i-k)}}---\left(6\right)$

其中，为修改后的梯度，k为当前迭代次数；j为加权的前期迭代次数，综合考虑效果和效率，令其等于10；e为自然常数，a为衰减因子，取为0.4。流程进入到步骤107。

步骤107，由修改后的梯度计算更新步长，如公式(7)所示，

${\mathrm{\α}}^k=\frac{{\left(g^k\right)}^{*}{g}^k}{{\left({\mathrm{Lg}}^k\right)}^{*}{\mathrm{Lg}}^k}---\left(7\right)$

其中，α^k为第k次迭代的更新步长，g^k为第k次迭代的梯度。流程进入到步骤108。

步骤108，由公式(8)更新反射系数模型，

m^k+1＝m^k-α^kP^kg^k (8)

其中，P^k为预处理算子，最优的预处理算子应为Hessian矩阵的逆，但计算量巨大，直接求取并不现实，本发明中用背景波场的能量来近似Hessian矩阵的对角元素，在减少计算量的同时加速了收敛速度。流程返回到步骤103。

步骤109，输出最终偏移成像结果。

实施例一

以Marmousi模型为例测试本发明所述方法，速度场如图2(a)所示，对应的扰动模型如图2(b)所示,品质因子模型如图2(c)所示。从图2可以看出，Marmousi模型中断块发育，速度变化剧烈，可用于检验成像算法的优劣。模型参数如下：横向网格点数为737，网格间距12.5m，纵向网格点数为375，网格间距8m。计算参数为：在地表以50m间隔均匀分布185炮，每炮都是737个检波点全接收，震源为主频20Hz的雷克子波，时间采样间隔0.6ms。

图3为成像测试结果，其中，图3(a)为采用传统极性编码得到的多震源粘声LSRTM结果，图3(b)为采用随机最优化极性编码的多震源粘声LSRTM成像结果，图3(c)采用随机最优化极性编码的多震源声波LSRTM成像结果。对比图3(a)和图3(b)可知，随机最优化算法相比传统算法成像质量有了明显改善，低频噪音及串扰噪音更弱，分辨率略高。对比图3(b)和图3(c)可以看出，由于考虑了地下介质的衰减特性，本发明算法相比随机最优化多震源声波LSRTM成像振幅更加均衡，尤其是在深部地区的照明补偿效果及分辨率更高。

相应的三种算法的数据残差曲线如图4所示(图中纵坐标为对数坐标)，其中“I”线为采用传统极性编码的多震源粘声LSRTM收敛曲线，“II”线为本发明算法，“III”线为采用随机最优化极性编码的多震源声波LSRTM收敛曲线。从图4可以清楚发现，本发明算法的误差收敛更快，在第30次迭代时的数据残差已经优于其他两种算法迭代第60次的数据残差，即采用本算法所需迭代次数更少，从而可显著节约计算量、提高计算效率。

Claims (9)

1.基于随机最优化的多震源粘声最小二乘逆时偏移方法，其特征在于，该基于随机最优化的多震源粘声最小二乘逆时偏移方法包括：

1)读取野外观测记录及预设参数合成超道集；

2)采用当前反射系数模型，通过多震源激发正演模拟超道集，计算数据残差；

3)根据数据残差计算更新梯度；

4)通过随机最优化思想修改梯度并计算更新步长；

5)由梯度及更新步长更新反射系数模型。

2.根据权利要求1所述的基于随机最优化的多震源粘声最小二乘逆时偏移方法，其特征在于，在步骤1，输入初始反射系数模型、偏移速度场、观测数据、品质因子、迭代终止的阈值及偏移参数，初始反射系数模型的值为0，即第1次迭代与常规粘声逆时偏移等价。

3.根据权利要求1所述的基于随机最优化的多震源粘声最小二乘逆时偏移方法，其特征在于，在步骤2，计算数据残差时，基于标准线性固体模型的扰动波场的控制方程为：

$\frac{1}{v_0^2}\frac{{\∂}^2{p}_s}{\∂t^2}=\frac{\∂}{\∂t}\[\mathrm{\ρ}\left(1-\underset{i=1}{\overset{I}{\Σ}}\left(1-\frac{{\mathrm{\τ}}_{\mathrm{\ϵ}i}}{{\mathrm{\τ}}_{\mathrm{\σ}i}}\right)e^{-t/{\mathrm{\τ}}_{\mathrm{\σ}i}}\right)H\left(t\right)\]*\▿\·(\frac{1}{\mathrm{\ρ}}\▿p_s)-m\left(x\right)\frac{{\∂}^2{p}_0}{\∂t^2}---\left(1\right)$

其中，p_s为扰动波场，v₀为背景速度，ρ为密度，I为标准线性固体的个数，τ_εi,τ_σi为松弛时间；H(t)为单位阶跃函数，▽为梯度算子，▽·为散度算子，*为时间上的卷积算子，m(x)为模型参数，即反射系数模型，p₀为背景波场，即在背景介质中传播的波，其控制方程为：

$\frac{1}{v_0^2}\frac{{\∂}^2{p}_0}{\∂t^2}=\frac{\∂}{\∂t}\[\mathrm{\ρ}\left(1-\underset{i=1}{\overset{I}{\Σ}}\left(1-\frac{{\mathrm{\τ}}_{\mathrm{\ϵ}i}}{{\mathrm{\τ}}_{\mathrm{\σ}i}}\right)e^{-t/{\mathrm{\τ}}_{\mathrm{\σ}i}}\right)H\left(t\right)\]*\▿\·(\frac{1}{\mathrm{\ρ}}\▿p_0)+f---\left(2\right)$

其中，f为经过编码的震源项，

公式(1)可用矩阵算子形式表示为:

p_s＝Lm (3)

其中，L为粘声介质扰动波场的线性正演算子。

4.根据权利要求1所述的基于随机最优化的多震源粘声最小二乘逆时偏移方法，其特征在于，在步骤2，在计算扰动波场和背景波场时，松弛时间τ_σ、τ_ε的计算公式如公式(4)所示：

$\left\{\begin{array}{c}{\mathrm{\τ}}_{\mathrm{\σ}}=\frac{\sqrt{Q^2+1}-1}{wQ}\\ {\mathrm{\τ}}_{\mathrm{\ϵ}}=\frac{\sqrt{Q^2+1}+1}{wQ}\end{array}\right.---\left(4\right)$

其中，w为圆频率，Q为品质因子。

5.根据权利要求1所述的基于随机最优化的多震源粘声最小二乘逆时偏移方法，其特征在于，在步骤3，通过数据残差计算更新模型的梯度方向，计算公式如公式(5)所示:

g＝L^*(Lm-p_obs) (5)

其中，g为梯度，p_obs为观测记录，L^*为正演算子的共轭转置，即逆时偏移算子，L为扰动波场的线性正演算子，m为反射系数模型。

6.根据权利要求1所述的基于随机最优化的多震源粘声最小二乘逆时偏移方法，其特征在于，该基于随机最优化的多震源粘声最小二乘逆时偏移方法还包括，在步骤3之后，判断梯度是否满足迭代终止条件，即梯度的模小于预设的阈值，若满足，输出当前反射系数模型，流程结束；否则流程进入到步骤4。

7.根据权利要求1所述的基于随机最优化的多震源粘声最小二乘逆时偏移方法，其特征在于，在步骤4，将随机最优化思想推广到相位编码的最小二乘逆时偏移算法中修改梯度，随机最优化方法需要加权平均之前的梯度，因此在前几次迭代时不需要修改梯度，修改后的梯度如公式(6)所示，

$\stackrel{\‾}{g^k}=\frac{{\mathrm{\Σ}}_{i=k-j}^{i=k}{e}^{a(i-k)}{g}^k}{{\mathrm{\Σ}}_{i=k-j}^{i=k}{e}^{a(i-k)}}---\left(6\right)$

其中，为修改后的梯度，g为梯度，k为当前迭代次数；j为加权的前期迭代次数，综合考虑效果和效率，令其等于10；e为自然常数，a为衰减因子，取为0.4。

8.根据权利要求1所述的基于随机最优化的多震源粘声最小二乘逆时偏移方法，其特征在于，在步骤4，由修改后的梯度计算更新步长，如公式(7)所示，

${\mathrm{\α}}^k=\frac{{\left(g^k\right)}^{*}{g}^k}{{\left({\mathrm{Lg}}^k\right)}^{*}{\mathrm{Lg}}^k}---\left(7\right)$

其中，α^k为第k次迭代的更新步长，g^k为第k次迭代的梯度，L为扰动波场的线性正演算子。

9.根据权利要求1所述的基于随机最优化的多震源粘声最小二乘逆时偏移方法，其特征在于，在步骤5，由公式(8)更新反射系数模型，

m^k+1＝m^k-α^kP^kg^k (8)

其中，P^k为预处理算子，m^k为第k次迭代的反射系数模型，α^k为第k次迭代的更新步长，g^k为第k次迭代的梯度，采用背景波场的能量来近似Hessian矩阵的对角元素，在减少计算量的同时加速了收敛速度，更新反射系数模型后流程返回到步骤2。