CN113031067B

CN113031067B - 一种基于Rytov-WKBJ近似的叠前地震反演方法

Info

Publication number: CN113031067B
Application number: CN202110207893.7A
Authority: CN
Inventors: 黄广谭
Original assignee: Zhejiang University ZJU
Current assignee: Zhejiang University ZJU
Priority date: 2021-02-24
Filing date: 2021-02-24
Publication date: 2022-05-27
Anticipated expiration: 2041-02-24
Also published as: CN113031067A

Abstract

本发明公开了一种基于Rytov‑WKBJ近似的叠前地震反演方法。包括以下步骤：步骤1：获取储层的共中心点叠前地震数据、叠前角道集数据、叠后地震数据和测井数据；步骤2：获得深度域地震数据；步骤3：利用地震层位数据和储层的测井数据构建深度域低频模型；步骤四：构建目标函数；步骤5：利用高斯‑牛顿法求解目标函数获得模型更新量，并更新模型扰动量；步骤6：若叠前地震角道集数据和合成地震数据两者的误差水平达到预设残差阈值，输出反演结果，反之进入步骤7；步骤7：更新迭代次数，返回步骤4进入下一次迭代，直至迭代终止。本发明应用于叠前地震反演中，可有效且快速的模拟强非均质性、复杂地质构造储层的复杂地震波场。

Description

一种基于Rytov-WKBJ近似的叠前地震反演方法

技术领域

本发明属于地震勘探技术领域，具体涉及一种基于Rytov-WKBJ近似的叠前地震反演方法。

背景技术

随着全球范围内对石油和天然气需求的不断增长，油气勘探方向逐步由传统的简单构造油气储层转向复杂地质构造油气储层^[1]。这些复杂地质构造油气储层往往具有较强的非均匀性及横向不连续性。目前常规地震反演算法仍主要针对地质构造相对简单的储层开展研究，因此如果将该类算法用于预测地质构造相对复杂的储层时，则会引起构造和岩性预测的偏差。对于构造较为复杂的储层而言，由于可能存在断层、裂缝、以及特殊地质体等原因，对储层的准确定位和描述提出了更好的要求^[2-3]。因此，针对目前复杂地质构造油气储层反演预测的技术瓶颈，需要开展针对性理论研究，以提高该类储层地震反演的预测精度，为后期的精准开发提供可靠的数据支撑。

不同于构造较为简单、地层平缓的地层，复杂构造储层常发育有断层、裂缝以及特殊地质体。而正是这些非均质体的存在，地震波场表现出较强的散射特征^[4-5]。在传统的地震储层弹性参数预测中，通常基于Zoeppritz及其近似式的AVO/AVA反演是目前应用最为广泛的技术方法。然而这种基于半空间无限介质、单一界面假设条件的方法，并不能用于模拟非均匀性较强的复杂地质构造储层的地震波场。而基于一维波动方程的反射率法反演，也无法克服这一缺陷。因此，基于散射积分方程求解的正演模拟方法被提出应用于地震反演算法当中。原则上讲，Green函数只能用于描述均一介质条件下点源所产生的波场，而对于复杂背景条件下Green函数是没有解析波场的^[6]。WKBJ近似的提出使得复杂构造条件下，求解散射积分方程成为可能，进而提出了Lippmann-Schiwinger方程^[7]。然而求解Lippmann-Schiwinger方程，需要消耗大量的计算量。因此Born近似的引入使得该方法得以在工业生产中得到应用。然而Born近似是通过对波场振幅的泰勒展开并省去高阶项得到的线性化公式，其地震波场模拟精度不高^[8-10]。因此基于Born-WKBJ近似的地震波正演模拟，很难用于准确模拟强扰动所引起的波场^[11-13]。不准确的波场模拟也会导致最终反演结果出现偏差。

模拟地震波在非均匀性强的复杂地质构造介质传播是影响地震反演结果质量的重要因素。在Lippman-Schwinger方程的基础上，引入更为精确的近似替代一阶Born近似是提高地震波模拟精度的关键因素。

[1]Sen,M.K.,and I.G.Roy,2003,Computation of differential seismogramsand iteration adaptive regularization in prestack waveform inversion:Geophysics,68,2026-2039.

[2]Abubakar,A.,P.M.van den Berg,and S.Y.Semenov,2004,A robustiterative method for born inversion:IEEE Transactions on Geoscience andRemote Sensing,42,342-354.

[3]Abubakar,A.,T.M.Habashy,P.M.van den Berg,and D.Gisolf,2005a,Thediagonalized contrast source approach:an inversion method beyond the bornapproximation:Inverse problems,21,685.

[4]陈生昌,曹景忠,马在田.稳定的Born近似叠前深度偏移方法[J].石油地球物理勘探,2001,36(3):291-296.

[5]陈生昌,曹景忠,马在田.基于Rytov近似的叠前深度偏移方法[J].石油地球物理勘探,2001,36(006):690-697.

[6]Fu,L.,2005,Rough surface scattering:Comparison of variousapproximation theories for 2D SH waves:Bulletin of the Seismological Societyof America,95,646-663.

[7]Gisolf,A.,P.Ha_nger,and P.Doulgeris,2017,Reservoir-oriented wave-equation-based seismic amplitude variation with o_set inversion:Interpretation,5,no.3,SL43-SL56.

[8]丁鹏程,李振春,张凯,等.基于Born/Rytov近似的波动方程层析速度反演方法研究[C]//中国地球科学联合学术年会.2015.

[9]Feng,R.,S.M.Luthi,D.Gisolf,and S.Sharma,2017,Obtaining a high-resolution geological and petrophysical model from the results of reservoir-orientated elastic wave-equation based seismic inversion:PetroleumGeoscience,23,376-385.

[10]Haffinger,P.,A.Gisolf,and P.v.d.Berg,2013,Towards high resolutionquantitative subsurface models by full waveform inversion:Geophysical JournalInternational,193,788-797.

[11]Huang,G.,C.Chen,Xiaohong and Luo,and X.Li,2019,Prestack waveforminversion by using an optimized linear inversion scheme:IEEE Transactions onGeoscience and Remote Sensing,57,no.8,5716-5728..

[12]丁伟,李振春,刘玉莲.基于Born/Rytov近似的联合叠前深度偏移方法[J].石油物探,2003(01):31-36.

[13]张凯,尹正,李振春,等.基于Bom/Rytov近似的波动方程层析速度反演方法研究[J].Applied Geophysics,2013(03):314-322.

发明内容

本发明针对非均匀性强的复杂地质构造地质模型的反演问题，提出了一种适用于近地表和浅部强非均匀性、复杂地质构造储层的叠前AVA反演方法，提高反演结果的精度。本发明属于叠前地震资料弹性参数反演范畴，克服传统叠前AVA/AVO反演基于单一反射界面的简单假设条件，可以适用于近地表及浅部复杂地质构造储层，高精度、高分辨率提取弹性三参数。

目前，叠前地震反演方法大多都基于水平层状或仅具有轻微地层起伏的油藏储层开展储层弹性参数提取的工作。而对强非均质、复杂地质条件的油藏储层没有开展针对性研究。传统方法无法考虑强非均质体所引起地震波散射，进而造成反演结果存在误差。本发明不仅可以克服传统方法对非均匀性、复杂地质构造介质的缺陷，有效提高反演结果的质量。同时通过线性化近似简化地震波场模拟，大大提高反演计算效率。

本发明采用的技术方案包括以下步骤：

步骤一：获取储层的共中心点叠前地震数据、叠前角道集数据d_obs、叠后地震数据和测井数据，叠前角道集数据d_obs包含全波场数据；

所述步骤一中，叠前角道集数据包含全波场响应的叠前角道集数据，即叠前角道集数据中包含有一次反射波、多次反射波、潜行波、绕射波等。本发明基于Rytov-WKBJ近似求解散射积分方程作为地震波场正演模拟的正演算子，散射积分方程的求解可以有效模拟复杂地质构造储层的地震全波场，是影响非均匀性强、复杂地质构造波场模拟的重要因素。

步骤二：对共中心点叠前地震数据进行动校正，获取均方根速度信息，对时间域的叠前角道集地震数据进行时-深转换，获得深度域地震数据；

基于波动方程的地震波场正演算子，将深度域的参数模型转化为时间域的地震数据。因此在合成角道集数据时需要提供深度域的弹性参数低频模型，作为正演算子的背景参数，同时也是反演过程中的初始模型。因此需要首先对深度域地震数据进行层位拾取，为低频模型构建奠定基础。

步骤三：对油气储层的叠后地震数据进行层位拾取，获得油气储层段附近的地震层位数据，并利用地震层位数据和储层的测井数据，基于反距离加权或克里金插值方法构建深度域低频模型；

所述低频模型由储层的拉梅常数λ₀、μ₀以及密度参数ρ₀构成；其中，λ₀、μ₀和ρ₀为将低频测井数据(低频滤波后的测井曲线)中的拉梅常数λ、μ和密度ρ进行空间插值获取的结果；

所述步骤三中，地震层位信息是进行测井数据插值时重要的横向约束条件，是低频模型构建的重要组成要素。该过程需消耗一定的人工拾取成本，其准确性至关重要，将直接影响低频模型构建的准确性，最终影响反演结果质量。

步骤四：构建目标函数；

定义目标函数如下：

其中，d_obs为叠前角道集数据(实际地震道集数据)，λ为拉梅常数；P_s为合成地震数据即散射波场；将低频模型作为背景模型m₀；χ_i-1表示反演模型m_i-1与背景模型m₀的差，即模型扰动量；i为迭代次数；

所述步骤四中的模型扰动量χ_i-1为反演模型m_i-1与背景模型m₀之间的差；背景模型m₀＝{Λ_r，M_r，ρ_r}；反演模型m_i-1＝{Λ，M，ρ}；模型扰动量χ_i-1通过下述计算得到：

其中，

ρ为密度参数；

对于第1次迭代运算，反演模型m_i-1即为背景模型m₀。

所述步骤四中的合成地震数据P_s由模型扰动量χ和背景地震波场P₀计算得到，具体过程如下：

1)通过推导得到合成地震数据P_s和全波场P之间的关系如下：

其中，j为复数单位，w表示角频率，Δr表示空间任意点到震源的距离，χ_i-1为模型扰动量，P和P_s分别代表全波场和合成地震数据；

全波场P等于背景地震波场P₀与合成地震数据P_s的加和，即：

P＝P₀+P_s

其中，背景地震波场P₀为背景模型所产生的地震波场；

所述背景地震波场P₀为基于背景模型m₀所产生的波场，背景地震波场P₀的散射积分方程为：

P₀＝∫G·(r-r_s)·Sdr

其中，P₀表示背景地震波场，G表示格林函数，r表示地震波接收点相对震源激发位置的传播距离，r_s表示震源位置，S表示震源处的波场值。

2)令

将其带入公式

中，并根据P＝P₀+P_s得到全波场：

其中，I表示单位对角矩阵；

3)由于直接求解步骤2)的全波长过于占用计算成本，故通过引入Rytov近似计算合成地震数据P_s：

首先对公式

进行泰勒展开得到：

保留一阶近似，则：

P＝P₀[1+φ(r)]

其中，φ表示相位角，r表示地震波接收点相对震源激发位置的传播距离；

然后根据Rytov近似，令P＝P₀exp[φ(r)]，并将

代入公式

后得线性化近似式：

可以发现相位角函数是一个与传播距离r和模型扰动量χ_i-1的函数，这里为了表示方便我们用

表示相位角随模型扰动量的变化，φ(r)表示相位角随传播距离的变化，将

代入P＝P₀[1+φ(r)]中，进而通过公式P＝P₀+P_s得到合成地震数据P_s：

令F＝GjwΔrP₀，对合成地震数据P_s进行简化得到：

步骤五：利用高斯-牛顿法求解步骤四的目标函数，获得模型更新量g_i，并更新模型扰动量；

所述步骤五中模型扰动量的更新值通过下述计算得到：

1)通过高斯-牛顿法求取模型更新量：

其中，

表示目标函数对模型扰动量χ_i-1的导数；

2)通过迭代更新获得模型扰动量的更新值：

χ_i＝χ_i-1+αg_i

其中，α表示步长，χ_i表示第i次迭代的模型扰动量，g_i表示第i次迭代模型更新量。

步骤六：根据叠前地震角道集数据d_obs和合成地震数据P_s，计算两者的误差水平是否达到预设的残差阈值；若达到预设的残差阈值，将更新的模型扰动量χ_i转换为对应的反演模型m_i，并输出反演结果m_i；反之，进入步骤七；

所述步骤六中，利用L2范数定义合成地震数据与叠前地震角道集数据的误差水平：

Error＝||d_obs-P_s||₂

其中，d_obs为实际地震道集数据，Error表示误差水平。

步骤七：更新迭代次数，返回步骤四进入下一次迭代，直至达到预设残差要求，则终止迭代过程。

本发明的有益效果是：

区别于传统基于等效反射系数公式的反演方法，本发明选用了Rytov-WKBJ近似算法求解散射积分方程作为叠前反演的正演算子。克服Lippmann-Schwinger方程的非线性及巨大计算量，同时基于散射积分方程方法将波动方程可以克服传统反射系数公式无法模拟非均质所引起的散射扰动的缺陷。因此可有效且快速的模拟强非均质性、复杂地质构造储层的复杂地震波场。

本发明创新性地引入Rytov近似提高地震波模拟精度和计算效率。该方法可很好的应用于强非均质、复杂地质条件的油气藏储层弹性参数反演中，可以很好的还原地下介质真实情况，从而降低反问题求解的不适定性，增强反演稳定性，提高反演结果的准确性，可以用于估算复杂地质构造储层的弹性参数，恢复地下地层属性特征。

附图说明

图1为本发明的实施流程图。

图2为模型测试的测井曲线及相对于背景参数的扰动量。

图3(a)自上至下分别采用反射率法、1阶Born-WKBJ近似、2阶Born-WKBJ近似、1阶Rytov-WKBJ近似以及2阶Rytov-WKBJ近似的近角度地震数据模拟结果。

图3(b)自上至下分别采用反射率法、1阶Born-WKBJ近似、2阶Born-WKBJ近似、1阶Rytov-WKBJ近似以及2阶Rytov-WKBJ近似的远角度地震数据模拟结果。

图4(a)是基于1阶Born-WKBJ近似求解散射积分方程获得方程的叠前反演三参数结果与理论模型对比示意图。

图4(b)是基于1阶Rytov-WKBJ近似求解散射积分方程获得方程的叠前反演三参数结果与理论模型对比示意图。

图4(c)是基于2阶Rytov-WKBJ近似求解散射积分方程获得方程的叠前反演三参数结果与理论模型对比示意图。

具体实施方式

下面结合附图对本发明作进一步的描述，并不能用来限制本发明的保护范围。

以一组测井数据测试为例进行说明，如图1所示，本发明的具体实施流程包括以下步骤：

一、基于Rytov-WKBJ近似的合成地震角道集数据的正演模拟

对储层的共中心点叠前地震数据进行动校正，获取均方根速度信息，对时间域的叠前角道集地震数据进行时-深转换，获得深度域地震数据；对油气储层的叠后地震数据进行层位拾取，获得油气储层段附近的地震层位数据，并利用地震层位数据和储层的测井数据，基于反距离加权或克里金插值方法构建深度域低频模型。

将低频模型m₀作为背景模型m₀，可得背景模型m₀所产生的背景地震波场P₀：

P₀＝∫G·(r-r_s)·Sdr (1)

其中，P₀表示背景模型所引起的地震波场，G表示格林函数，r表示空间任意点位置，r_s表示震源位置，S表示震源处的波场值。可推导出散射波场P_s和全波场P之间的关系：

其中，j为复数单位，w表示角频率，Δr表示空间任意点到震源的距离，χ表示反演模型与低频模型的差，即模型扰动量，P代表全波场地震数据。而全波场数据等于背景模型所产生的地震波场与散射波场的加和，即

P＝P₀+P_s (3)

因此令

可得全波场值：

其中，I表示单位对角矩阵。

但直接求解方程(4)过于占用计算成本，故引入Rytov近似对算法予以改进。首先对式(4)进行泰勒展开可得：

保留一阶近似，则P＝P₀[1+φ(r)]。

根据Rytov近似，令P＝P₀exp[φ(r)]，其中φ表示相位角，r表示地震波距离源点处的传播距离，将

代入公式(5)可得线性化近似式：

φ(r)＝GiwΔr·χ (6)

进而通过公式(3)，得到散射波场即合成地震数据

P_s＝GjwΔr·χ_i-1P₀ (7)

简化为

图2展示了Marmousi模型的第一道模型数据(包括体积模量的倒数Λ，剪切模量的倒数M，以及密度参数ρ)，及真实模型相对于背景模型的扰动χ_Λ，χ_M，χ_ρ。实线表示真实模型数据，虚线表示背景模型数据。图3则为采用几种常用的散射积分方程求解得到的地震波场值。实线为由反射率法获取的标准模拟结果，该方法可获得一维波动方程的解析解，是一维模型正、反演问题中地震数据模拟的标准。通过对比可以发现，本发明提出的方法(2阶Rytov-WKBJ近似)的准确性在模拟强扰动以及深部地区具有更好的模拟精度。

二、基于Rytov-WKBJ近似的叠前反演方法

1)定义目标函数如下：

其中，d_obs为叠前角道集数据，λ为拉梅常数；P_s为合成地震数据即散射波场；

将低频模型作为背景模型m₀；χ_i-1表示反演模型m_i-1与背景模型m₀的差，即模型扰动量；i为迭代次数；背景模型m₀＝{Λ_r，M_r，ρ_r}；反演模型m_i-1＝{Λ，M，ρ}；模型扰动量χ_i-1通过下述计算得到：

其中，

ρ为密度参数。

其中，λ、μ为拉梅常数，λ₀、μ₀和ρ₀为将低频滤波后的测井曲线(拉梅常数λ、μ和密度参数ρ)进行空间插值获取的结果；

2)利用高斯-牛顿法求取模型更新量：

其中，

表示目标函数对参数χ_i-1的导数。

模型参数的迭代更新公式为：

χ_i＝χ_i-1+αg_i (12)

其中α表示步长，χ_i表示第i次迭代的更新量，g_i表示模型更新量。

3)根据叠前地震角道集数据d_obs和合成地震数据P_s，计算两者的误差水平是否达到预设的残差阈值；若达到预设的残差阈值，根据公式(10)将将更新的模型扰动量χ_i转换为对应的反演模型m_i，并输出反演结果m_i；反之，返回步骤1)中的公式(9)。

利用L2范数定义合成地震数据与实际地震道集数据的误差水平控制迭代循环：

Error＝||d_obs-P_s||₂ (13)

其中，d_obs为实际地震道集数据，Error表示误差水平。

图4表示利用传统Born-WKBJ近似和本发明的Rytov-WKBJ近似得到的反演结果。观测数据由反射率法获取，该方法可获得一维波动方程的解析解，是一维模型正、反演问题中地震数据模拟的标准；其中，图中虚线表示理论模型，实现表示反演模型。比较可以发现，将本发明提出的Rytov-WKBJ近似求解散射积分方程的方法应用于叠前地震数据反演问题，可以获得更高精度的反演结果。

Claims

1.一种基于Rytov-WKBJ近似的叠前地震反演方法，其特征在于，包括以下步骤：

所述低频模型由储层的拉梅常数λ₀、μ₀以及密度参数ρ₀构成；其中，λ₀、μ₀和ρ₀为将低频测井数据中的拉梅常数λ、μ和密度ρ进行空间插值获取的结果；

步骤四：构建目标函数；

定义目标函数如下：

J(χ_i-1)＝||d_obs-P_s||₂+λ||χ_i-1||₂

将低频模型作为背景模型m₀；χ_i-1表示反演模型m_i-1与背景模型m₀的差，即模型扰动量；i为迭代次数；

步骤七：更新迭代次数，返回步骤四进入下一次迭代，直至达到预设残差要求，则终止迭代过程；

1)合成地震数据P_s和全波场P之间的关系如下：

P_s＝jωΔr·Gχ_i-1P

其中，j为复数单位，ω表示角频率，Δr表示空间任意点到震源的距离，χ_i-1为模型扰动量，P和P_s分别代表全波场和合成地震数据；

P＝P₀+P_s

其中，背景地震波场P₀为背景模型所产生的地震波场；

2)令V(χ_i-1)＝jωΔr·χ_i-1，将其带入公式P_s＝jωΔr·Gχ_i-1P中，并根据P＝P₀+P_s得到全波场：

P＝[I-GV(χ_i-1)]^-1P₀

其中，I表示单位对角矩阵；

3)通过引入Rytov近似计算合成地震数据P_s：

首先对公式P＝[I-GV(χ_i-1)]^-1P₀进行泰勒展开得到：

保留一阶近似，则：

P＝P₀[1+φ(r)]

然后根据Rytov近似，令P＝P₀ exp[φ(r)]，并将V(χ_i-1)＝jωΔr·χ_i-1代入公式

后得线性化近似式：

φ(χ_i-1)＝GjωΔr·χ_i-1

将φ(χ_i-1)＝GjωΔr·χ_i-1代入P＝P₀[1+φ(r)]中，进而通过公式P＝P₀+P_s得到合成地震数据P_s：

P_s＝GjωΔr·χ_i-1P₀

令F＝GjωΔrP₀，对合成地震数据P_s进行简化得到：

P_s＝Fχ_i-1。

2.根据权利要求1所述的一种基于Rytov-WKBJ近似的叠前地震反演方法，其特征在于，所述步骤四中的模型扰动量χ_i-1为反演模型m_i-1与背景模型m₀之间的差；背景模型m₀＝{Λ_r，M_r，ρ_r}；反演模型m_i-1＝{Λ，M，ρ}；模型扰动量χ_i-1通过下述计算得到：

其中，

ρ为密度参数；

对于第1次迭代运算，反演模型m_i-1即为背景模型m₀。

3.根据权利要求1所述的一种基于Rytov-WKBJ近似的叠前地震反演方法，其特征在于，所述背景地震波场P₀为基于背景模型m₀所产生的波场，背景地震波场P₀的散射积分方程为：

P₀＝∫G·(r-r_s)·Sdr

4.根据权利要求1所述的一种基于Rytov-WKBJ近似的叠前地震反演方法，其特征在于，所述步骤五中模型扰动量的更新值通过下述计算得到：

1)通过高斯-牛顿法求取模型更新量：

其中，

表示目标函数对模型扰动量χ_i-1的导数；

2)通过迭代更新获得模型扰动量的更新值：

χ_i＝χ_i-1+αg_i

5.根据权利要求1所述的一种基于Rytov-WKBJ近似的叠前地震反演方法，其特征在于，所述步骤六中，利用L2范数定义合成地震数据与叠前地震角道集数据的误差水平：

Error＝||d_obs-P_s||₂

其中，d_obs为实际地震道集数据，Error表示误差水平。