CN111999764B - 基于时频域目标函数的盐下构造最小二乘逆时偏移方法 - Google Patents

Abstract

本发明公开了一种基于时频域目标函数的盐下构造最小二乘逆时偏移方法，利用地震信号的时频域局部化特征和弱地震信号的时频域振幅相位信息，提高盐下构造的成像精度。首先，对地震信号进行时频变换，构建时频域最小二乘偏移目标函数；其次，引入振幅调节因子，增强时频域目标函数中相位信息的权重，减弱地震数据振幅信息对深部成像结果影响，提高盐下构造对应的弱地震信号的成像精度；最后，推导时频域振幅相位联合目标函数对应的梯度算子，并利用优化算法对扰动模型进行更新迭代。通过盐丘模型的数值测试，验证了本专利方法可以获得盐下构造的高精度成像结果。

Description

基于时频域目标函数的盐下构造最小二乘逆时偏移方法

技术领域

本发明涉及一种地下构造逆时偏移方法，具体说是一种基于时频域目标函数的盐下构造最小二乘逆时偏移方法。

背景技术

随着油气勘探技术的不断发展，对地下构造探测精度需求逐渐提高，如今正从构造勘探阶段逐步转向岩性勘探阶段。偏移方法在地下构造成像中占有重要的地位，多年来为了不断地提高地下构造的成像精度，从基于射线理论的Kirchhoff偏移和波束偏移方法，逐渐发展为基于波动理论的偏移方法(Baysal等，1983)。在复杂构造成像过程中，逆时偏移方法具有成像精度高的优点，能够适用于横向速度变化剧烈的区域，获取地下模型高精度反射系数(Xu等，2011)。但是逆时偏移算子是正演算子的共轭转置，而不是它的逆(Claerbout 1992)。因此在观测范围有限的情况下，基于互相关成像条件的逆时偏移方法会产生低频噪音，降低分辨率，无法获得准确的反射系数。

针对逆时偏移成像过程中存在的问题，发展了最小二乘逆时偏移方法，并在反演的理论框架下，利用反偏移数据来不断地与观测数据相匹配，并通过优化算法更新迭代，最终获取地下高精度反射系数信息。关于最小二乘偏移方法，首先由Nemeth等(1999)在最小二乘目标函数的框架下实现了Kirchhoff偏移。随后，Dai等(2011，2012)利用多源地震数据进行最小二乘逆时偏移测试，很大程度上提高了计算效率。任浩然等(2013)将Hessian算子应用于最小二乘逆时偏移中，得到相对振幅保真的成像结果，改善了地震偏移成像的精度。Tan等(2014)在原始最小二乘逆时偏移方法的基础上，提出了波场分离成像条件，并通过更新震源波场实现陡倾角断层的高精度成像。为了克服初始模型依赖问题，刘玉金等(2015)提出扩展成像条件下的最小二乘逆时偏移方法，测试结果表明该方法可以获得更准确的振幅属性信息。此外，为了减弱地震数据振幅信息对成像结果的影响，在互相关目标函数的基础上提出了很多改进的偏移方法(Zhang等，2015；Liu等2016；李庆洋等，2016；Liu等，2017；Yi等，2019)。陈生昌和周华敏(2018)在反射波动方程的基础上，考虑了入射波传播算子与反射波传播算子对偏移成像作用的不同。方修政等(2018)指出基于常规互相关成像条件的最小二乘逆时偏移方法梯度含有很强的低频噪音，为此提出了基于逆散射成像条件的最小二乘逆时偏移方法方法。巩向博等(2019)利用多源地震数据，在稀疏约束的理论框架下，针对小尺度地质构成像问题造展开研究。Liu等(2020)在单步最小二乘逆时偏移方法的理论框架下，加入了Gabor反褶积滤波，增强模拟数据与观测数据的可匹配程度，进而提高盐下构造的成像精度。Li等(2020)通过引入权重因子衰减地震数据中较强的反射波信息，进而实现增强地下的散射信号，有效地提高了深部构造的成像质量。经过多年的研究，最小二乘逆时偏移方法技术得到了很快地发展，逐渐成为当前的研究热点问题。但是地震波场在地下传播是一个极其复杂的过程，受到粘弹介质的影响，致使在数据的匹配过程中的振幅信息对最小二乘逆时偏移结果影响较为严重。

发明内容

本发明所要解决的技术问题在于，克服现有技术存在的缺陷，利用地震信号的时频域局部化特征，提出了一种基于时频域目标函数的盐下构造最小二乘逆时偏移方法，充分利用弱地震信号的时频域相位信息，提高地构造的成像精度。

本发明基于时频域目标函数的盐下构造最小二乘逆时偏移方法，其步骤是：

步骤1.使用MATLAB软件及其地震数据处理软件Crews工具包对地震数据进行预处理，并将预处理后的地震数据作为时频域最小二乘逆时偏移的数据输入。

步骤2.构建初始速度模型，作为波动方程正演模拟的模型输入。

步骤3.定义观测系统，提取震源子波；利用初始速度模型进行波动方程正演模拟，存储正传波场，并与观测数据的反传波场互相关，获得逆时偏移成像结果；在速度模型上正演模拟对应的声波方程为：

其中x,z为坐标轴，s表示慢度场，u地震波场，f为震源子波，t为时间；综合计算效率和模拟精度，优选用有限差分法进行波动方程正演模拟；

步骤4.计算扰动介质产生的扰动地震波场，扰动介质表示为：

将慢度场的平方s²表示为背景慢度的平方

与扰动慢度平方Δs²的叠加。对应的总地震波场同样分解为背景地震波场u₀和扰动地震波场u_s的叠加，即：

u＝u₀+u_s

其中背景地震波场u₀与总地震波场u都满足波动方程：

当扰动地震波场u_s＜＜u₀时，可以用背景地震波场代替总地震波场u₀≈u₀+u_s，即：

则要想求解扰动场u_s，必须先求解背景场u₀，然后将扰动项作为新的震源继续传播。

步骤5.利用Gabor变换获得观测数据和模拟数据的时频域振幅相位信息。其中观测数据和模拟数据Gabor变换为：

其中h(τ-t)表示高斯窗函数；t表示时间；τ为中间变量；ω为角频率，u(t)和d(t)分别表示时间域模拟数据和观测数据；

和

表示时频域模拟数据和观测数据，F_h[·]表示对地震数据作用的Gabor变换算子；

步骤6.构建基于时频域振幅相位的最小二乘逆时偏移目标函数：

其中i表示虚数；ns为震源数目；nr为检波器数目；模型扰动

和

表示模拟数据和观测数据在时频域的相位信息；

和

表示观测数据和模拟数据在时频域的振幅信息；ε∈[0,1]表示权重因子，用来控制目标函数中波形振幅和波形相位的比例。当ε＝0时，目标函数为纯相位目标函数，完全忽略了振幅的影响；当ε＝1时，目标函数将退化为常规最小二乘目标函数。为了获得目标函数对模型扰动的更新量，需要求目标函数关于模型参数的偏导数：

其中

*为复共轭；Re[·]取数据的实部；

是数据的残差。其中时频域模拟数据的绝对值对速度参数的偏导数可以表示为：

则目标函数关于模型参数的偏导数可以表示为：

由于

存在如下变换关系：

目标函数的梯度算子可以表示为：

改变公式对应的积分顺序，最终目标函数的梯度算子为：

步骤7.利用Gabor反变换将时频域的梯度表达式转换到时间域中，其中为了对公式做进一步的简化，地震观测数据和模拟数据Gabor逆变换可以表示为：

其中

表示Gabor逆变换算子。因此将中对应的部分，目标函数的梯度可以进一步在时空域简化为：

步骤8.定义伴随震源，并将其反传至模型空间。对应的伴随震源为：

步骤9.利用零延迟互相关方法获得最小二乘逆时偏移的成像梯度。根据波动方程的表达式，Born正演模拟可以用矩阵的形式来表示：

A_su_s＝A₀u₀

其中

要想获得扰动地震波场相对于模型扰动的偏导数，只需要对上式两边同时计算其关于模型扰动的偏导数：

由于A_s和u₀与模型扰动无关，则上式变为：

同样将梯度写成矩阵的形式为：

步骤10.利用L-BFGS优化算法计算模型扰动的更新方向。迭代公式如下：

m_k+1＝m_k-α_kH_kg_k

其中m_k为第k步模型扰动的更新量，α_k为步长，H_k为近似海森矩阵的逆，

为模型扰动的更新梯度。在L-BFGS优化更新中，只需要保存少数的向量对，用于更新Hessian矩阵，其更新公式如下：

H_k+1＝V_k ^TH_kV_k+ρ_ks_ks_k ^T

V_k＝I-ρ_ky_ks_k ^T

s_k＝m_k+1-m_k，y_k＝g_k+1-g_k

其中H_k+1是根据向量对{s_k,y_k}和H_k计算得到；H_kg_k的乘积可以通过梯度g_k与向量对{s_k,y_k}之间一系列向量的内积与向量的和来获得的。其中的近似Hessian矩阵的逆矩阵H_k需满足以下更新公式：

步骤11.判断是否满足终止条件，若满足则输出基于时频域目标函数的最小二乘逆时偏移结果。若不满足终止条件，将当前的成像结果继续作为下一次循环的初始扰动模型，直到满足终止条件。

本发明方法，充分利用弱地震信号的时频域振幅相位信息，提高地下构造的成像精度。首先，对地震信号进行时频变换，构建时频域最小二乘偏移目标函数；其次，引入振幅调节因子，增强时频域目标函数中相位信息的权重，减弱地震数据振幅信息对反演结果影响，提高盐下构造对应的弱地震信号的成像精度；最后，推导时频域振幅相位联合目标函数对应的梯度算子，并利用优化算法对反射系数模型进行更新迭代。通过盐丘模型的数值测试，验证了本专利方法可以获得盐下构造的高精度成像结果。

附图说明

图1是本发明基于时频域目标函数的盐下构造最小二乘逆时偏移方法流程框图。

图2是盐丘模型，(a)真实盐丘速度模型；(b)初始盐丘速度模型；(c)真实模型扰动。

图3是偏移结果，(a)逆时偏移结果；(b)最小二乘逆时偏移结果；(c)基于时频域目标函数的最小二乘逆时偏移结果。

具体实例

由于盐丘速度模型中存在强速度扰动构造，对地震波的传播具有很强的屏蔽作用，因此很难获得盐丘下部构造的高精度成像结果。为了解决这个问题，本专利使用基于时频域目标函数的最小二乘逆时偏移成像方法，充分利用弱地震信号的时频域相位信息，提高盐下构造的成像精度。真实的盐丘速度模型如图2a所示，初始速度模型如图2b所示，真实的模型扰动如图2c所示。在地表均匀分布50个震源，每个震源对应着600个检波器。地震数据记录时间为5s，时间间隔为2ms。震源采用主频为8Hz的雷克子波，同时利用L-BFGS优化算法来计算下降方向，并对模型扰动进行更新迭代。

本专利使用模型数据进行测试，首先利用真实速度模型(图2a)正演获得观测数据，并将其视为实际生产过程使用的地震数据。根据前述技术方案流程，首先对地震数据进行预处理，构建初始速度模型(图2b)，定义观测系统，以及提取地震数据对应的震源子波信息。为了获得与观测数据相匹配的地震数据，利用Born正演模拟，并在Gabor时频变换的框架下获得时频域模拟数据和观测数据。在优化理论下，构建给予时频域振幅相位的最小二乘逆时偏移目标函数，并求取目标函数相对于模型扰动的偏导数，进而获得模型扰动的更新量。最后利用L-BFGS优化算法在初始模型的基础上不断更新模型扰动，获得高精度的成像结果(图3)。

模型参数如下：

表1：基于时频域目标函数的最小二乘逆时偏移方法测试参数

将图3a与图3b对比可以看出，最小二乘逆时偏移结果相对于常规的逆时偏移结果在深部构造成像上有了明显的改善，尤其是盐丘下部的弱散射构造的成像结果明显优于常规逆时偏移方法。同时最小二乘逆时偏移结果在浅层区域有效地减弱了震源采集脚印的影响，成像结果在能量方面有了很好地均衡，更接近真实模型扰动。其中逆时偏移成像结果在振幅上与真实模型扰动相差较远，这主要是因为地震波在传播过程中能量迅速衰减，深部反射信号较弱，成像效果较差。图3c是本专利提出的基于时频域目标函数的最小二乘逆时偏移成像结果。将图3c与图3b和3a相对比可以看出，基于时频域目标函数的最小二乘逆时偏移成像结果在深部区域成像方面有了明显地改善。由于时频域振幅相位信息具有与地下速度模型更好的线性对应关系，极大程度上减弱了最小二乘逆时偏移对初始速度模型的依赖。此外，在时频域目标函数中引入权重因子，不仅增强了反演过程的稳定性，而且缓解扰动波场匹配过程中出现的周波跳跃现象。数值试验结果表明，基于时频域目标函数的最小二乘逆时偏移方法在深部构造高精度成像上具有一定的优势。

Claims (2)

1.一种基于时频域目标函数的盐下构造最小二乘逆时偏移方法，其步骤是：

步骤1.使用MATLAB软件及其地震数据处理软件Crews工具包对地震数据进行预处理，并将预处理后的地震数据作为时频域最小二乘逆时偏移的数据输入；

步骤2.构建初始速度模型，作为波动方程正演模拟的模型输入；

步骤3.定义观测系统，提取震源子波；利用初始速度模型进行波动方程正演模拟，存储正传波场，并与观测数据的反传波场互相关，获得逆时偏移成像结果；在速度模型上正演模拟对应的声波方程为：

其中x,z为坐标轴，s表示慢度场，u为地震波场，f为震源子波，t为时间；

步骤4.计算扰动介质产生的扰动地震波场，扰动介质表示为：

将慢度场的平方s²表示为背景慢度的平方

与扰动慢度平方Δs²的叠加；对应的总地震波场同样分解为背景地震波场u₀和扰动地震波场u_s的叠加，即：

u＝u₀+u_s

其中背景地震波场u₀与总地震波场u都满足波动方程：

当扰动地震波场u_s＜＜u₀时，用背景地震波场代替总地震波场u₀≈u₀+u_s，即：

步骤5.利用Gabor变换获得观测数据和模拟数据的时频域振幅相位信息；其中观测数据和模拟数据Gabor变换为：

其中h(τ-t)表示高斯窗函数；t表示时间；τ为中间变量；ω为角频率，u(t)和d(t)分别表示时间域模拟数据和观测数据；

和

表示时频域模拟数据和观测数据，F_h[·]表示对地震数据作用的Gabor变换算子；

步骤6.构建基于时频域振幅相位的最小二乘逆时偏移目标函数：

其中i表示虚数；ns为震源数目；nr为检波器数目；模型扰动

和

表示模拟数据和观测数据在时频域的相位信息；

和

表示观测数据和模拟数据在时频域振幅信息；ε∈[0,1]表示权重因子，用来控制目标函数中振幅和相位的比例；当ε＝0时，目标函数为纯相位目标函数，完全忽略了波形振幅的影响；当ε＝1时，目标函数将退化为常规最小二乘目标函数；求目标函数关于模型参数的偏导数，获得目标函数对模型扰动的更新量：

步骤7.利用Gabor反变换将时频域的梯度表达式转换到时间域中，其中为了对公式做进一步的简化，地震观测数据和模拟数据Gabor逆变换表示为：

其中

表示Gabor逆变换算子；因此，目标函数的梯度进一步在时空域简化为：

步骤8.定义伴随震源，并将其反传至模型空间；对应的伴随震源为：

步骤9.利用零延迟互相关方法获得最小二乘逆时偏移的梯度；根据波动方程的表达式，

Born正演模拟用矩阵的形式来表示：

A_su_s＝A₀u₀

其中

计算其关于模型扰动的偏导数：

同样将梯度写成矩阵的形式为：

步骤10.利用L-BFGS优化算法计算模型扰动的更新方向；迭代公式如下：

m_k+1＝m_k-α_kH_kg_k

其中m_k为第k步模型扰动的更新量，α_k为步长，H_k为近似海森矩阵的逆，

为模型扰动的梯度；在L-BFGS优化更新中，只需要保存少数的向量对，用于更新Hessian矩阵，其更新公式如下：

H_k+1＝V_k ^TH_kV_k+ρ_ks_ks_k ^T

V_k＝I-ρ_ky_ks_k ^T

s_k＝m_k+1-m_k，y_k＝g_k+1-g_k

其中H_k+1是根据向量对{s_k,y_k}和H_k计算得到；H_kg_k的乘积通过梯度g_k与向量对{s_k,y_k}之间一系列向量的内积与向量的和来获得；其中的近似Hessian矩阵的逆矩阵H_k需满足以下更新公式：

步骤11.判断是否满足终止条件，若满足则输出基于时频域目标函数的最小二乘逆时偏移结果；若不满足终止条件，将当前的成像结果继续作为下一次循环的初始扰动模型，直到满足终止条件。

2.根据权利要求1所述的基于时频域目标函数的盐下构造最小二乘逆时偏移方法，其特征是：采用有限差分法进行波动方程正演模拟。

Images