CN108873066B - 弹性介质波动方程反射波旅行时反演方法 - Google Patents

弹性介质波动方程反射波旅行时反演方法 Download PDF

Info

Publication number
CN108873066B
CN108873066B CN201810666471.4A CN201810666471A CN108873066B CN 108873066 B CN108873066 B CN 108873066B CN 201810666471 A CN201810666471 A CN 201810666471A CN 108873066 B CN108873066 B CN 108873066B
Authority
CN
China
Prior art keywords
equation
wave
objective function
inversion
background
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810666471.4A
Other languages
English (en)
Other versions
CN108873066A (zh
Inventor
任志明
李振春
黄少华
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China University of Petroleum East China
Original Assignee
China University of Petroleum East China
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China University of Petroleum East China filed Critical China University of Petroleum East China
Priority to CN201810666471.4A priority Critical patent/CN108873066B/zh
Publication of CN108873066A publication Critical patent/CN108873066A/zh
Application granted granted Critical
Publication of CN108873066B publication Critical patent/CN108873066B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Investigating Or Analyzing Materials By The Use Of Ultrasonic Waves (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明公开了一种弹性介质波动方程反射波旅行时反演方法。设计新的旅行时残差目标函数；推导新目标函数下的伴随方程和背景参数梯度公式；计算背景参数的梯度；采用拟牛顿法反演算法对梯度进行处理；采用抛物线拟合法求取迭代步长；更新背景参数模型，直到满足收敛条件。本发明的有益效果是通过采用新的旅行时目标函数来提高弹性介质反射波波形反演的精度，获取可靠的纵波速度和横波速度初始模型。

Description

弹性介质波动方程反射波旅行时反演方法
技术领域
本发明属于地震波反演技术领域，涉及勘探地震和天然地震中弹性参数(纵、横波速度和密度)初始模型建立，来提高多分量弹性地震资料反演和成像的精度。
背景技术
地球物理学的基本问题就是通过研究各种地球物理场的特征来获取地球内部的结构和构造信息。地震反演通过地表或井中观测到的地震数据估算相应的地球物理参数，进而反推地下的结构形态及物质成分，可以有效地识别地质构造、预测自然灾害和勘探油气藏。随着勘探程度的不断深入，油气勘探的目标已从大型构造油气藏向中小型复杂隐蔽型油藏转变。面对日益复杂的地质目标，提高地下岩性参数的反演精度至关重要。全波形反演利用地震波的动力学和运动学信息，通过求观测数据与模拟数据误差的极小来获取地下模型参数(如速度、密度、品质因子、各向异性参数等)，是地球深部构造成像、浅层环境调查及油气勘探速度建模的有效工具。
由于存在周期跳跃(cycle-skipping)现象，全波形反演对初始模型有较强的依赖性。当初始模型与真实模型相差较远时，很难得到满意的反演结果。反射波波形反演是只采用反射波场来获取地下模型低频分量的新方法。反射波波形反演中，模型参数被拆分为背景参数和扰动参数/反射系数，反射波的振幅和相位/走时分别由反射系数和背景速度场决定，通过交替进行逆时偏移/最小二乘逆时偏移和波形反演来更新反射系数和背景速度模型。与常规全波形反演相比，反射波波形反演的敏感核/梯度更加光滑，能获得可靠的低频速度模型。反射波波形反演可以为全波形反演提供可靠的初始模型，进而提高地下模型参数的反演精度。因此发展反射波波形反演方法具有较大的理论意义和实用价值。
现有的反射波波形反演方法大多是基于声波方程设计的，只能得到纵波速度的低频模型。针对复杂介质(如弹性、各向异性)的反射波波形反演方法相对较少。与声波方程相比，弹性波方程能够较好地处理纵、横波之间的模式转换，更接近野外的实际情况。因此，弹性波全波形反演具有更好的应用前景。与此同时，弹性波全波形反演的非线性和多解性更强，对初始模型的依赖性也更强。本发明将给出一种新的弹性反射波旅行时反演方法，来同时获取纵波速度和横波速度低频模型，进而提高后续多分量地震资料全波形反演和逆时偏移的反演和成像精度。
发明内容
本发明的目的在于提供弹性介质波动方程反射波旅行时反演方法。
本发明所采用的技术方案是按照以下步骤进行：
(1)设计新的旅行时残差目标函数；
采用PP反射波旅行时残差构建目标函数；
(2)推导新目标函数下的伴随方程和背景参数梯度公式；
基于伴随方法(Adjoint method)推导伴随方程及目标函数对背景参数的梯度公式；
(3)计算背景参数的梯度；
具体包括：背景波场正向传播；扰动波场正向传播；反射波旅行时残差反向传播；扰动波场反向传播；正向和反向波场相关得到背景参数的梯度；
(4)采用拟牛顿法反演算法对梯度进行处理；
(5)采用抛物线拟合法求取迭代步长；
(6)更新背景参数模型，直到满足收敛条件。
进一步，弹性介质波动方程反射波旅行时反演步骤(1)中，目标函数为：
Figure GDA0002153109720000021
其中，τ1和τ2分别表示模拟和观测反射波水平分量和垂直分量旅行时残差。T为最大记录时间。
采用动态图像规整法(dynamic image warping)来估算局部旅行时时差。极小化函数为：
Figure GDA0002153109720000022
Figure GDA0002153109720000023
其中，Δvxobs和Δvzobs观测反射波的水平分量和垂直分量，Δvx和Δvz,模拟反射波的水平分量和垂直分量。
方程(2)关于自变量求偏导数等于0得：
Figure GDA0002153109720000031
其中，变量上方的“点”表示对时间的偏导数。
基于方程(4)和隐函数的求导法则可以得到：
Figure GDA0002153109720000032
同理可得：
Figure GDA0002153109720000033
伴随方程的虚源为目标函数对模拟波场的偏导数，则
Figure GDA0002153109720000034
Figure GDA0002153109720000035
为缓解不同波模式之间串扰对反演精度的影响，我们只采用反射纵波场(PP波)构建目标函数。本发明中通过求解纵横波分离的弹性波方程来获取PP波场。
步骤(2)中推导新目标函数下的伴随方程和背景参数梯度公式方法如下：
二维弹性波速度-应力方程：
Figure GDA0002153109720000036
Figure GDA0002153109720000037
Figure GDA0002153109720000038
Figure GDA0002153109720000039
Figure GDA00021531097200000310
其中，λ和μ为拉梅常数，ρ为密度，(vx,vz)为质点振动速度矢量，(τxx,τzz,τxz)为应力矢量。
弹性介质中，对于背景参数[λ,μ,ρ]，背景波场P＝[vx,vz,τxx,τzz,τxz]可以通过求解方程(9)得到。当存在模型扰动[Δλ,Δμ,Δρ]时，扰动/反射波场为ΔP＝[Δvx,Δvz,Δτxx,Δτzz,Δτxz]，且满足：
Figure GDA0002153109720000041
Figure GDA0002153109720000042
Figure GDA0002153109720000043
Figure GDA0002153109720000044
Figure GDA0002153109720000045
方程(10)与方程(9)相减，化简并忽略高阶微小量可得：
Figure GDA0002153109720000046
Figure GDA0002153109720000047
Figure GDA0002153109720000048
Figure GDA0002153109720000049
Figure GDA00021531097200000410
对给定的参数扰动[Δλ,Δμ,Δρ]，求解方程(9)和方程(11)可以得到反射波[Δvx,Δvz,Δτxx,Δτzz,Δτxz]，即为弹性介质中的反偏移过程。
方程(1)给出基于旅行时的目标函数，背景波场P满足方程(9)，扰动/反射波场ΔP满足方程(11)。采用拉格朗日乘子法求解该约束优化问题。目标函数变为：
Figure GDA00021531097200000411
Figure GDA0002153109720000051
其中，[ξx,ξz,ξxx,ξzz,ξxz]和[ψx,ψz,ψxx,ψzz,ψxz]为拉格朗日乘子函数
Figure GDA0002153109720000052
Figure GDA0002153109720000053
Figure GDA0002153109720000054
Figure GDA0002153109720000055
Figure GDA0002153109720000056
和
Figure GDA0002153109720000057
Figure GDA0002153109720000058
Figure GDA0002153109720000059
Figure GDA00021531097200000510
Figure GDA00021531097200000511
分步积分方程(12)并令
Figure GDA00021531097200000518
和
Figure GDA00021531097200000512
得：
Figure GDA00021531097200000513
Figure GDA00021531097200000514
Figure GDA00021531097200000515
Figure GDA00021531097200000516
Figure GDA00021531097200000517
和
Figure GDA0002153109720000061
Figure GDA0002153109720000062
Figure GDA0002153109720000063
Figure GDA0002153109720000064
Figure GDA0002153109720000065
目标函数对背景参数的梯度为：
Figure GDA0002153109720000066
Figure GDA0002153109720000067
Figure GDA0002153109720000068
进一步，步骤(3)中
a.求解方程(9)和方程(11)得到正向背景波场[vx,vz,τxx,τzz,τxz]T和扰动波场[Δvx,Δvz,Δτxx,Δτzz,Δτxz]T，初始条件为：
[vx(x,z,0),vz(x,z,0),τxx(x,z,0),τzz(x,z,0),τxz(x,z,0)]T＝0,
[Δvx(x,z,0),Δvz(x,z,0),Δτxx(x,z,0),Δτzz(x,z,0),Δτxz(x,z,0)]T＝0(18)
b.求解伴随方程(15)和(16)得到旅行时残差反向背景波场[ξx,ξz,ξxx,ξzz,ξxz]和扰动波场[ψx,ψz,ψxx,ψzz,ψxz]，终值条件为：
[ξx(x,z,T),ξz(x,z,T),ξxz(x,z,T),ξzz(x,z,T),ξxz(x,z,T)]T,
[ψx(x,z,T),ψz(x,z,T),ψxz(x,z,T),ψzz(x,z,T),ψxz(x,z,T)]T. (19)
c.通过方程(17)计算目标函数关于背景参数的梯度。
进一步，步骤(4)中处理方法如下：
采用L-BFGS法：
Figure GDA0002153109720000071
其中，Hk为海森矩阵逆的近似矩阵，直接计算Hk需要较大的计算量，这里通过几组列向量来近似Hk。
进一步，步骤(5)中求取迭代步长方法如下：
采用抛物线拟合求取迭代步长
Figure GDA0002153109720000072
其中，α1和α2为试探步长，J1和J2为相应的目标函数值，J0为当前迭代的目标函数值。
则当前迭代的最佳步长为：
Figure GDA0002153109720000073
进一步，步骤(6)中通过下式更新背景参数：
Figure GDA0002153109720000074
其中，mk和mk+1分别为当前迭代和下一次迭代的模型参数：
Figure GDA0002153109720000075
重复步骤(3)-(6)，直到满足收敛条件。初始反射系数模型往往也是不正确的。为得到满意的反演结果，在反射波波形反演中嵌套保幅性较好的最小二乘逆时偏移。通过交替进行最小二乘逆时偏移和波形反演来分别更新反射系数和背景弹性参数模型。
本发明的有益效果是通过弹性反射波旅行时反演来提高弹性参数初始模型的建模精度。发明的目的是获取精确的纵、横波速度低频模型，为后续的多分量资料弹性波全波形反演和逆时偏移提供可靠的初始模型。
附图说明
图1弹性介质波动方程反射波旅行时反演的流程图；
图2双层介质模型；
图3双层介质模型不同反演方法的敏感核；
图4Sigsbee2A模型；
图5Sigsbee2A模型不同反演方法的反演结果。
具体实施方式
下面结合具体实施方式对本发明进行详细说明。
如图1所示，为本发明实施弹性介质波动方程反射波旅行时反演的流程图，具体包括：
(1)设计新的旅行时残差目标函数；
采用PP反射波旅行时残差构建目标函数；
(2)推导新目标函数下的伴随方程和背景参数梯度公式；
基于伴随方法(Adjoint method)推导伴随方程及目标函数对背景参数的梯度公式；
(3)计算背景参数的梯度；
具体包括：背景波场正向传播；扰动波场正向传播；反射波旅行时残差反向传播；扰动波场反向传播；正向和反向波场相关得到背景参数的梯度；
(4)采用拟牛顿法反演算法对梯度进行处理；
(5)采用抛物线拟合法求取迭代步长；
(6)更新背景参数模型；
步骤(1)中设计新的旅行时残差目标函数方法如下：
与波形相比，地震波的旅行时更依赖于速度模型的长波长分量。利用旅行时残差可以一定程度提高低频模型的反演精度。本发明中采用旅行时信息进行弹性反射波反演。目标函数变为：
Figure GDA0002153109720000091
其中，τ1和τ2分别表示模拟和观测反射波水平分量和垂直分量旅行时残差。T为最大记录时间。
采用动态图像规整法(dynamic image warping)来估算局部旅行时时差。极小化函数为：
Figure GDA0002153109720000092
Figure GDA0002153109720000093
其中，Δvxobs和Δvzobs观测反射波的水平分量和垂直分量，Δvx和Δvz,模拟反射波的水平分量和垂直分量。
方程(2)关于自变量求偏导数等于0得：
Figure GDA0002153109720000094
其中，变量上方的“点”表示对时间的偏导数；
基于方程(4)和隐函数的求导法则可以得到：
Figure GDA0002153109720000095
同理可得：
Figure GDA0002153109720000096
伴随方程的虚源为目标函数对模拟波场的偏导数，则
Figure GDA0002153109720000097
Figure GDA0002153109720000098
为缓解不同波模式之间串扰对反演精度的影响，我们只采用反射纵波场(PP波)构建目标函数。本发明中通过求解纵横波分离的弹性波方程来获取PP波场。
步骤(2)推导新目标函数下的伴随方程和背景参数梯度公式方法如下：
二维弹性波速度-应力方程：
Figure GDA0002153109720000101
Figure GDA0002153109720000102
Figure GDA0002153109720000103
Figure GDA0002153109720000104
Figure GDA0002153109720000105
其中，λ和μ为拉梅常数，ρ为密度，(vx,vz)为质点振动速度矢量，(τxx,τzz,τxz)为应力矢量。
弹性介质中，对于背景参数[λ,μ,ρ]，背景波场P＝[vx,vz,τxx,τzz,τxz]可以通过求解方程(9)得到。当存在模型扰动[Δλ,Δμ,Δρ]时，扰动/反射波场为ΔP＝[Δvx,Δvz,Δτxx,Δτzz,Δτxz]，且满足：
Figure GDA0002153109720000106
Figure GDA0002153109720000107
Figure GDA0002153109720000108
Figure GDA0002153109720000109
Figure GDA00021531097200001010
方程(10)与方程(9)相减，化简并忽略高阶微小量可得：
Figure GDA00021531097200001011
Figure GDA00021531097200001012
Figure GDA0002153109720000111
Figure GDA0002153109720000112
Figure GDA0002153109720000113
对给定的参数扰动[Δλ,Δμ,Δρ]，求解方程(9)和方程(11)可以得到[Δvx,Δvz,Δτxx,Δτzz,Δτxz]，即为弹性介质中的反偏移过程。
方程(1)给出基于旅行时的目标函数，背景波场P满足方程(9)，扰动/反射波场ΔP满足方程(11)。采用拉格朗日乘子法求解该约束优化问题。目标函数变为：
Figure GDA0002153109720000114
其中，[ξx,ξz,ξxx,ξzz,ξxz]和[ψx,ψz,ψxx,ψzz,ψxz]为拉格朗日乘子函数
Figure GDA0002153109720000115
Figure GDA0002153109720000116
Figure GDA0002153109720000117
Figure GDA0002153109720000118
Figure GDA0002153109720000119
和
Figure GDA00021531097200001110
Figure GDA00021531097200001111
Figure GDA00021531097200001112
Figure GDA00021531097200001113
Figure GDA0002153109720000121
分步积分方程(12)并令
Figure GDA00021531097200001215
和
Figure GDA00021531097200001216
得：
Figure GDA0002153109720000122
Figure GDA0002153109720000123
Figure GDA0002153109720000124
Figure GDA0002153109720000125
Figure GDA0002153109720000126
和
Figure GDA0002153109720000127
Figure GDA0002153109720000128
Figure GDA0002153109720000129
Figure GDA00021531097200001210
Figure GDA00021531097200001211
目标函数对背景参数的梯度为：
Figure GDA00021531097200001212
Figure GDA00021531097200001213
Figure GDA00021531097200001214
步骤(3)中计算背景参数的梯度方法如下：
a.求解方程(9)和方程(11)得到正向背景波场[vx,vz,τxx,τzz,τxz]T和扰动波场[Δvx,Δvz,Δτxx,Δτzz,Δτxz]T，初始条件为：
[vx(x,z,0),vz(x,z,0),τxx(x,z,0),τzz(x,z,0),τxz(x,z,0)]T＝0,
[Δvx(x,z,0),Δvz(x,z,0),Δτxx(x,z,0),Δτzz(x,z,0),Δτxz(x,z,0)]T＝0(18)
b.求解伴随方程(15)和(16)得到旅行时残差反向背景波场[ξx,ξz,ξxx,ξzz,ξxz]和扰动波场[ψx,ψz,ψxx,ψzz,ψxz]，终值条件为：
[ξx(x,z,T),ξz(x,z,T),ξxz(x,z,T),ξzz(x,z,T),ξxz(x,z,T)]T,
[ψx(x,z,T),ψz(x,z,T),ψxz(x,z,T),ψzz(x,z,T),ψxz(x,z,T)]T. (19)
c.通过方程(17)计算目标函数关于背景参数的梯度。
步骤(4)中采用合适的反演算法对梯度进行处理方法如下：
对梯度进行预处理可以提高反演的精度和收敛速度。常用的方法有共轭梯度法、拟牛顿法和牛顿法等。为了兼顾反演精度和计算效率，本发明中采用L-BFGS法：
Figure GDA0002153109720000131
其中，Hk为海森矩阵逆的近似矩阵。直接计算Hk需要较大的计算量，这里通过几组列向量来近似Hk。
步骤(5)求取迭代步长方法如下：
本发明中采用抛物线拟合求取迭代步长：
Figure GDA0002153109720000132
其中，α1和α2为试探步长，J1和J2为相应的目标函数值，J0为当前迭代的目标函数值。
则当前迭代的最佳步长为：
Figure GDA0002153109720000141
步骤(6)更新背景参数模型方法如下：
基于上面的步骤，通过下式更新背景参数：
Figure GDA0002153109720000142
其中，mk和mk+1分别为当前迭代和下一次迭代的模型参数：
Figure GDA0002153109720000143
重复步骤(3)-(6)，直到满足收敛条件(如残差小于1e-5或迭代次数小于30等)。初始反射系数模型往往也是不正确的。为得到满意的反演结果，在反射波波形反演中嵌套保幅性较好的最小二乘逆时偏移。通过交替进行最小二乘逆时偏移和波形反演来分别更新反射系数和背景弹性参数模型。
本发明由于采取以上技术方案，其具有以下优点：1.可以提高纵、横波速度初始模型的建模精度。2.可以改善多分量地震资料弹性波全波形反演和逆时偏移的反演和成像精度。
下面通过几个例子来分析本发明中提出的弹性介质波动方程反射波旅行时反演方法的有效性。
首先以一个炮点-检波点对为例来说明本发明的优势。双层介质模型如图2所示，炮点和检波点位于三角形指示的位置。时间步长1ms，空间间隔10m，震源为15Hz的雷克子波，加在正应力上。图3为不同反演方法的敏感核。常规反射波波形反演的敏感核(a)λ和(b)μ；波动方程反射波旅行时反演方法的敏感核(c)λ和(d)μ；波动方程反射波旅行时反演方法的PP敏感核(e)λ和(f)μ。由图可知，新提出的反射波旅行时反演方法的敏感核比常规反射波波形反演方法的敏感核更加光滑，有利于低频速度模型的恢复。通过纵、横波分离可以有效地压制不同波模式之间串扰产生的高频噪声。
下面采用复杂的Sigsbee2A模型(如图4所示)对新提出的弹性介质波动方程反射波旅行时反演方法进行测试。时间步长1.5ms，空间间隔16m，炮点(36炮)和检波点均匀分布于地表。震源为10Hz的雷克子波，加在正应力上。初始模型为随深度线性变化的一维模型(与真实模型相差较远)。图5给出Sigsbee2A模型不同反演方法的反演结果。常规反射波波形反演结果(a)vP和(b)vS；波动方程反射波旅行时反演结果(c)vP和(d)vS；以波动方程反射波旅行时反演结果为初始模型的常规反射波波形反演结果(e)vP和(f)vS。由图可知，直接采用常规全波形反演很难得到满意的反演结果；通过本发明中提出的反射波旅行时反演方法可以获取较好的纵波速度和横波速度模型初始模型，进而采用常规全波形反演来恢复模型的高频分量。
本发明是一种新的地震波反演方法，能够获取可靠的弹性参数(纵、横波速度和密度)初始模型；能大大提高多分量地震资料弹性波全波形反演和逆时偏移的反演和成像的精度。
以上所述仅是对本发明的较佳实施方式而已，并非对本发明作任何形式上的限制，凡是依据本发明的技术实质对以上实施方式所做的任何简单修改，等同变化与修饰，均属于本发明技术方案的范围内。

Claims (1)

1.弹性介质波动方程反射波旅行时反演方法，其特征在于按照以下步骤进行：
(1)设计新的旅行时残差目标函数；
采用PP反射波旅行时残差构建目标函数；
(2)推导新目标函数下的伴随方程和背景参数梯度公式；
基于伴随方法推导伴随方程及目标函数对背景参数的梯度公式；
(3)计算背景参数的梯度；
包括：背景波场正向传播；扰动波场正向传播；反射波旅行时残差反向传播；扰动波场反向传播；正向和反向波场相关得到背景参数的梯度；
(4)采用拟牛顿法反演算法对梯度进行处理；
(5)采用抛物线拟合法求取迭代步长；
(6)更新背景参数模型，直到满足收敛条件；
步骤(1)中设计新的目标函数：
Figure FDA0002255474530000011
其中，τ1表示模拟和观测反射波水平分量旅行时残差，τ2表示模拟和观测反射波垂直分量旅行时残差，T为最大记录时间；
采用动态图像规整法来估算局部旅行时时差，极小化函数为：
Figure FDA0002255474530000012
Figure FDA0002255474530000013
其中，Δvxobs和Δvzobs分别为观测反射波的水平分量和垂直分量，Δvx和Δvz分别为模拟反射波的水平分量和垂直分量；
方程(2)关于自变量求偏导数等于0得：
Figure FDA0002255474530000014
其中，变量Δvxobs上方的·表示对时间的偏导数；
基于方程(4)和隐函数的求导法则可以得到：
Figure FDA0002255474530000021
同理可得：
Figure FDA0002255474530000022
伴随方程的虚源为目标函数对模拟波场的偏导数，则
Figure FDA0002255474530000023
Figure FDA0002255474530000024
为缓解不同波模式之间串扰对反演精度的影响，我们只采用反射纵波场PP波构建目标函数，通过求解纵横波分离的弹性波方程来获取PP波场；
步骤(2)中推导新目标函数下的伴随方程和背景参数梯度公式方法如下：
二维弹性波速度-应力方程：方程(9)；
Figure FDA0002255474530000025
Figure FDA0002255474530000026
Figure FDA0002255474530000027
Figure FDA0002255474530000028
Figure FDA0002255474530000029
其中，λ和μ为拉梅常数，ρ为密度，(vx,vz)为质点振动速度矢量，(τxx,τzz,τxz)为应力矢量；弹性介质中，对于背景参数[λ,μ,ρ]，背景波场P＝[vx,vz,τxx,τzz,τxz]可以通过求解方程(9)得到，当存在参数扰动[Δλ,Δμ,Δρ]时，扰动/反射波场为ΔP＝[Δvx,Δvz,Δτxx,Δτzz,Δτxz]，且满足：方程(10)；
Figure FDA0002255474530000031
Figure FDA0002255474530000032
Figure FDA0002255474530000033
Figure FDA0002255474530000034
Figure FDA0002255474530000035
方程(10)与方程(9)相减，化简并忽略高阶微小量可得方程(11)：
Figure FDA0002255474530000036
Figure FDA0002255474530000037
Figure FDA0002255474530000038
Figure FDA0002255474530000039
Figure FDA00022554745300000310
对给定的参数扰动[Δλ,Δμ,Δρ]，求解方程(9)和方程(11)可以得到[Δvx,Δvz,Δτxx,Δτzz,Δτxz]，即为弹性介质中的反偏移过程；
方程(1)给出基于旅行时的目标函数，背景波场P满足方程(9)，扰动/反射波场ΔP满足方程(11)，采用拉格朗日乘子法求解约束优化问题，目标函数变为：
Figure FDA00022554745300000311
其中，[ξx,ξz,ξxx,ξzz,ξxz]和[ψx,ψz,ψxx,ψzz,ψxz]为拉格朗日乘子函数
Figure FDA00022554745300000312
Figure FDA0002255474530000041
Figure FDA0002255474530000042
Figure FDA0002255474530000043
Figure FDA0002255474530000044
和
Figure FDA0002255474530000045
Figure FDA0002255474530000046
Figure FDA0002255474530000047
Figure FDA0002255474530000048
Figure FDA0002255474530000049
分步积分方程(12)并令
Figure FDA00022554745300000410
和
Figure FDA00022554745300000411
得方程(15)：
Figure FDA00022554745300000412
Figure FDA00022554745300000413
Figure FDA00022554745300000414
Figure FDA00022554745300000415
Figure FDA00022554745300000416
和方程(16)；
Figure FDA00022554745300000417
Figure FDA0002255474530000051
Figure FDA0002255474530000052
Figure FDA0002255474530000053
Figure FDA0002255474530000054
目标函数对背景参数的梯度为方程(17)：
Figure FDA0002255474530000055
Figure FDA0002255474530000056
Figure FDA0002255474530000057
步骤(3)中
a.求解方程(9)和方程(11)得到正向背景波场[vx,vz,τxx,τzz,τxz]T和扰动/反射波场[Δvx,Δvz,Δτxx,Δτzz,Δτxz]T，初始条件为：
[vx(x,z,0),vz(x,z,0),τxx(x,z,0),τzz(x,z,0),τxz(x,z,0)]T＝0,
[Δvx(x,z,0),Δvz(x,z,0),Δτxx(x,z,0),Δτzz(x,z,0),Δτxz(x,z,0)]T＝0 (18)
b.求解伴随方程(15)和(16)得到[ξx,ξz,ξxx,ξzz,ξxz]和[ψx,ψz,ψxx,ψzz,ψxz]，终值条件为：
[ξx(x,z,T),ξz(x,z,T),ξxz(x,z,T),ξzz(x,z,T),ξxz(x,z,T)]T,
[ψx(x,z,T),ψz(x,z,T),ψxz(x,z,T),ψzz(x,z,T),ψxz(x,z,T)]T (19)
c.通过方程(17)计算目标函数关于背景参数的梯度；
步骤(4)中处理方法如下：
采用L-BFGS法：
Figure FDA0002255474530000061
其中，Hk为海森矩阵逆的近似矩阵，直接计算Hk需要较大的计算量，通过几组列向量来近似Hk；
步骤(5)中求取迭代步长方法如下：
采用抛物线拟合求取迭代步长
Figure FDA0002255474530000062
其中，α1和α2为试探步长，J1和J2为相应的目标函数值，J0为当前迭代的目标函数值，计算J1和J2需要额外的四次正演运算；
则当前迭代的最佳步长为：
Figure FDA0002255474530000063
步骤(6)中通过下式更新背景参数：
Figure FDA0002255474530000064
其中，mk和mk+1分别为当前迭代和下一次迭代的背景参数：
Figure FDA0002255474530000065
为得到满意的反演结果，在反射波波形反演中嵌套保幅性较好的最小二乘逆时偏移，通过交替进行最小二乘逆时偏移和波形反演来分别更新反射系数和背景参数模型。
CN201810666471.4A 2018-06-26 2018-06-26 弹性介质波动方程反射波旅行时反演方法 Active CN108873066B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810666471.4A CN108873066B (zh) 2018-06-26 2018-06-26 弹性介质波动方程反射波旅行时反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810666471.4A CN108873066B (zh) 2018-06-26 2018-06-26 弹性介质波动方程反射波旅行时反演方法

Publications (2)

Publication Number Publication Date
CN108873066A CN108873066A (zh) 2018-11-23
CN108873066B true CN108873066B (zh) 2020-03-10

Family

ID=64295685

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810666471.4A Active CN108873066B (zh) 2018-06-26 2018-06-26 弹性介质波动方程反射波旅行时反演方法

Country Status (1)

Country Link
CN (1) CN108873066B (zh)

Families Citing this family (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110187382B (zh) * 2019-03-05 2020-10-13 中国石油大学(华东) 一种回折波和反射波波动方程旅行时反演方法
CN110261897B (zh) * 2019-04-26 2021-07-20 中国石油化工股份有限公司 基于组稀疏的叠前四参数反演方法
CN110058307B (zh) * 2019-05-05 2020-12-18 四川省地质工程勘察院集团有限公司 一种基于快速拟牛顿法的全波形反演方法
CN110285876A (zh) * 2019-07-01 2019-09-27 中国人民解放军军事科学院国防科技创新研究院 一种海洋声场全波解的获取方法
CN110764146B (zh) * 2019-10-24 2020-05-19 南京信息工程大学 一种基于声波算子的空间互相关弹性波反射波形反演方法
CN110954950B (zh) * 2019-10-31 2022-10-21 南方科技大学 地下横波速度反演方法、装置、计算设备及存储介质
CN111766628A (zh) * 2020-07-29 2020-10-13 浪潮云信息技术股份公司 一种预条件的时间域弹性介质多参数全波形反演方法
CN113126151B (zh) * 2021-03-10 2022-06-07 山东省科学院海洋仪器仪表研究所 一种基于纯波延拓方程的弹性反射波旅行时反演方法
CN113376689B (zh) * 2021-03-30 2022-04-12 中国铁路设计集团有限公司 一种考虑层间多次波的弹性反射波走时反演方法
CN113311484B (zh) * 2021-05-26 2022-02-25 中国矿业大学(北京) 利用全波形反演获取粘弹性介质弹性参数的方法及装置
CN114460646B (zh) * 2022-04-13 2022-06-28 山东省科学院海洋仪器仪表研究所 一种基于波场激发近似的反射波旅行时反演方法
CN116626751B (zh) * 2023-05-26 2024-03-19 北京中矿大地地球探测工程技术有限公司 基于多目标函数的黏弹性参数同步反演方法、装置和设备
CN116699695B (zh) * 2023-08-07 2023-11-03 北京中矿大地地球探测工程技术有限公司 一种基于衰减矫正的反演方法、装置和设备

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP3209859B1 (en) * 2014-10-24 2021-04-28 Schlumberger Technology B.V. Travel-time objective function for full waveform inversion
CN106324662B (zh) * 2015-06-15 2018-08-07 中国石油化工股份有限公司 一种针对目标层的全波形反演方法及系统
CN107843925B (zh) * 2017-09-29 2019-03-08 中国石油化工股份有限公司 一种基于修正相位的反射波波形反演方法

Also Published As

Publication number Publication date
CN108873066A (zh) 2018-11-23

Similar Documents

Publication Publication Date Title
CN108873066B (zh) 弹性介质波动方程反射波旅行时反演方法
CN108139499B (zh) Q-补偿的全波场反演
CN110187382B (zh) 一种回折波和反射波波动方程旅行时反演方法
EP1616202A2 (en) Seimsic imaging by wave migration using a krylov space expansion of the square root exponent operator
CN110007340B (zh) 基于角度域直接包络反演的盐丘速度密度估计方法
CN109188519B (zh) 一种极坐标下的弹性波纵横波速度反演系统及方法
US10310117B2 (en) Efficient seismic attribute gather generation with data synthesis and expectation method
CN110579795B (zh) 基于被动源地震波形及其逆时成像的联合速度反演方法
GB2499096A (en) Simultaneous joint estimation of P-P and P-S residual statics
AU2011312806A1 (en) Hybrid method for full waveform inversion using simultaneous and sequential source method
CN111025388A (zh) 一种多波联合的叠前波形反演方法
CN111505714B (zh) 基于岩石物理约束的弹性波直接包络反演方法
Zhong et al. Elastic reverse time migration method in vertical transversely isotropic media including surface topography
Jin et al. 2D multiscale non‐linear velocity inversion
CN111999764A (zh) 基于时频域目标函数的盐下构造最小二乘逆时偏移方法
CN114624766B (zh) 基于行波分离的弹性波最小二乘逆时偏移梯度求取方法
CN112630830A (zh) 一种基于高斯加权的反射波全波形反演方法及系统
EP4080250A1 (en) Multiple attenuation and imaging processes for recorded seismic data
CN111175822B (zh) 改进直接包络反演与扰动分解的强散射介质反演方法
Shi et al. Multiscale full-waveform inversion based on shot subsampling
Zhao et al. Double plane wave reverse time migration in the time domain
Przebindowska Acoustic full waveform inversion of marine reflection seismic data
Broggini et al. FD injection utilizing the wavefields generated by Marchenko redatuming: A target-oriented approach
Majdi et al. Application of finite difference eikonal solver for traveltime computation in forward modeling and migration
Huang Modeling and inversion of seismic data using multiple scattering, renormalization and homotopy methods

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant