CN110687600B - 一种基于声-弹耦合方程的弹性波最小二乘逆时偏移方法 - Google Patents

一种基于声-弹耦合方程的弹性波最小二乘逆时偏移方法 Download PDF

Info

Publication number
CN110687600B
CN110687600B CN201911016294.6A CN201911016294A CN110687600B CN 110687600 B CN110687600 B CN 110687600B CN 201911016294 A CN201911016294 A CN 201911016294A CN 110687600 B CN110687600 B CN 110687600B
Authority
CN
China
Prior art keywords
wave
equation
elastic
component
elastic coupling
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
CN201911016294.6A
Other languages
English (en)
Other versions
CN110687600A (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.)
Nanjing University of Information Science and Technology
Original Assignee
Nanjing University of Information Science and Technology
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 Nanjing University of Information Science and Technology filed Critical Nanjing University of Information Science and Technology
Priority to CN201911016294.6A priority Critical patent/CN110687600B/zh
Publication of CN110687600A publication Critical patent/CN110687600A/zh
Application granted granted Critical
Publication of CN110687600B publication Critical patent/CN110687600B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30Analysis
    • G01V1/306Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/36Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/30Noise handling
    • G01V2210/32Noise reduction
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/50Corrections or adjustments related to wave propagation
    • G01V2210/51Migration
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/60Analysis
    • G01V2210/62Physical property of subsurface
    • G01V2210/622Velocity, density or impedance
    • G01V2210/6226Impedance

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

本发明提供一种基于声‑弹耦合方程的弹性波最小二乘逆时偏移方法,首先,输入初始纵横波速度和密度初始模型,利用声‑弹耦合传播方程及其伴随方程,合成炮端正传波场和检端反传波场,并根据纵、横波速度和密度的成像条件求取多参数像;然后,根据声‑弹耦合方程的一阶Born散射方程,反偏移合成反射波数据;再次,根据四分量波形残差目标函数衡量成像结果,若不满足则利用对角Hessian预条件的CGLS算法,迭代更新纵横波速度和密度成像结果,直至满足条件,返回最终结果。

Description

一种基于声-弹耦合方程的弹性波最小二乘逆时偏移方法
技术领域
本发明涉及一种面向海底四分量地震数据的成像方法,具体为一种基于声-弹耦合方程的弹性波最小二乘逆时偏移方法属于地震波成像技术领域。
背景技术
海底地震勘探逐渐成为海上油气勘探、开发和监测的一种重要地球物理手段。与常规海面拖缆观测相比,海底地震勘探具有信噪比高、检波器可精确定位、可满足重复性观测需求、易开展宽方位观测且受海况影响小等诸多优势。到目前为止,海底四分量(4C)地震数据的采集和前期处理技术,例如去噪、重定位、重定向、水陆检合并等,已经比较成熟,但是,面向4C地震数据的处理手段仍然非常匮乏,尤其是高精度的地震成像方法遇到了较大的挑战。
弹性波逆时偏移(ERTM)从观测的地震数据中提取地下构造信息,且不受地层倾角、偏移孔径的限制,是当今地震勘探最准确、最常用的成像方法。但是,ERTM通常建立在标准位移(速度)-应力弹性波方程之上,该方程不能正确模拟压力分量数据,无法实现4C数据的有效利用。通过将压力和应力张量的关系引入弹性波方程,声-弹耦合方程解决了海底固-液分层介质中地震波的精确模拟,并且可以同时提供位移和压力多种分量记录。目前,基于该方程开发的ERTM方法可以证明了可以有效压制反传过程中产生的非物理的成像假象。
然而,实际地震勘探中观测孔径是有限的,分布是不规则的,子波频带也是有限的,因此通常只能得到分辨率较低、串扰噪音严重、保幅性较差的成像结果。最小二乘偏移(LSM)方法将地震偏移转化为最小二乘优化问题,通过最佳匹配观测数据获取分辨率更高、成像假象更少且振幅保真好的成像结果。由于其优异的反演性能,被广泛应用于声波、弹性波介质成像中。但是,近几年发展的弹性波最小二乘逆时偏移(ELSRTM)方法只能用于三分量地震数据上,尚未开发面向海底四分量地震数据的弹性波最小二乘逆时偏移算法。
发明内容
针对现有海底四分量地震数据ERTM方法存在空间分辨率较低以及串扰假象严重等问题,本申请在声-弹耦合方程的框架下,提供一种面向海底四分量地震数据的最小二乘逆时偏移方法。
一种基于声-弹耦合方程的弹性波最小二乘逆时偏移方法,具体包括以下步骤:
(1)输入纵波速度、横波速度和密度的偏移模型;
(2)根据二维声-弹耦合方程计算炮端正传波场;
(3)输入多分量地震观测数据,计算多分量伴随震源,并根据声-弹耦合伴随方程计算检端反传波场;
(4)根据炮端正传波场和检端反传波场以及成像条件,计算多参数像;
(5)根据炮端正传波场、多参数像以及声-弹耦合方程的一阶Born散射方程计算反偏移波场,提取多分量合成数据;
(6)根据4C数据最小二乘偏移的目标泛函求取当前目标函数,若目标函数小于阀值(理论模型阀值可以设为1%),则输出当前多参数像作为最终成像结果,进入步骤(9);若不满足,进入步骤(7);
(7)根据最小二乘共轭梯度算法进行迭代,计算更新后的多参数像;
(8)根据对角Hessian算子,计算更新后的预条件多参数像,并返回(5);
(9)根据输出的最终多参数像,计算纵、横波阻抗的反射率。
进一步地,在步骤(2)中,二维声-弹耦合方程如下:
Figure BDA0002245813590000021
Figure BDA0002245813590000022
Figure BDA0002245813590000023
Figure BDA0002245813590000024
Figure BDA0002245813590000025
其中,正传弹性波场包括5个分量
Figure BDA0002245813590000026
依次代表:水平和垂直位移分量、压力分量以及与横波相关的正应力分量和切应力分量,fP代表正传时施加在压力分量上的子波,λ和μ为拉梅常数,ρ为密度;根据该方程计算炮端正传弹性波场,并在检波点处抽取多分量合成数据,包括多分量
Figure BDA0002245813590000031
水平分量
Figure BDA0002245813590000032
和垂直分量
Figure BDA0002245813590000033
注:本申请采用弹性波交错网格高阶差分技术进行地震波正传和后文提到的反传及反偏移模拟,并采用卷积最佳匹配层技术实现吸收边界。
进一步地,在步骤(3)中,输入多分量地震观测数据,包括压力分量
Figure BDA0002245813590000034
水平分量
Figure BDA0002245813590000035
和垂直分量
Figure BDA0002245813590000036
并利用步骤2中得到的炮端正传多分量合成数据,计算多分量伴随震源压力分量f’p,水平分量f’x和垂直分量f’z,其计算公式如下,
Figure BDA0002245813590000037
并根据声-弹耦合伴随方程,计算检端反传波场
Figure BDA0002245813590000038
Figure BDA0002245813590000039
Figure BDA00022458135900000310
Figure BDA00022458135900000311
Figure BDA00022458135900000312
其中,反传波场也包括的5个分量
Figure BDA00022458135900000313
依次代表:水平和垂直位移分量、压力分量以及与横波相关的正应力分量和切应力分量。
进一步地,在步骤(4)中,根据炮端正传波场和检端反传波场以及成像条件,计算多参数像,其中拉梅常数和密度的成像公式为:
Figure BDA00022458135900000314
Figure BDA00022458135900000315
Figure BDA00022458135900000316
根据链式法则,纵、横波速度和密度的成像公式如下,
δα=2αρ·δλ
δβ=-4βρ·δλ+2βρ·δμ。
δρ=(α2-2β2)·δλ+β2·δμ+δρ
进一步地,步骤(5)中,声-弹耦合方程的一阶Born散射方程为:
Figure BDA0002245813590000041
Figure BDA0002245813590000042
Figure BDA0002245813590000043
Figure BDA0002245813590000044
Figure BDA0002245813590000045
其中,
Figure BDA0002245813590000046
为反偏移波场的5个分量;
并根据炮端正传波场和多参数像,计算反偏移波场,在检波点处抽取合成的反射波数据压力分量
Figure BDA0002245813590000047
和位移分量
Figure BDA0002245813590000049
其中位移分量
Figure BDA00022458135900000410
包括水平分量
Figure BDA00022458135900000411
和垂直分量
Figure BDA00022458135900000412
进一步地,步骤(6)中,4C数据最小二乘偏移目标泛函E(m)定义为:
Figure BDA00022458135900000413
其中,m代表模型参数,du和dP分别表示位移分量和压力分量,上标syn和obs分别代表合成数据和观测数据,ε是权重系数,ε∈[0,1],用于调节压力数据和位移数据;ζ是尺度因子,用于平衡位移分量和压力分量之间的量纲差异,可以利用步骤2中获取的位移分量ux、uz和压力分量P进行计算,公式如下:
Figure BDA00022458135900000414
若目标函数小于阀值(理论模型阀值可以设为1%),则输出当前多参数像作为最终成像结果,进入步骤(9);若不满足,进入步骤(7);
进一步地,在步骤(7)中,最小二乘共轭梯度CGLS算法如下:
下列变量中,xk代表当前多参数像,xk+1代表更新后(即下一代)多参数像,其余变量均为数学算法的中间变量不具备明确的物理意义;
在第1次迭代前,需设定初值x0=0;r0=Hx0-mmig,p0=-r0,k=0;
然后按下列公式求解更新后的成像结果xk+1
Figure BDA0002245813590000051
xk+1=xkkpk,
rk+1=rkkHpk=rkk[LT(Lpk)],
Figure BDA0002245813590000052
pk+1=-rk+1k+1pk,
需要注意的是,在每一步迭代过程中,Hessian向量乘Hpk需要被拆解为一步反偏移Lpk,实现方式已在步骤5中阐述,偏移过程LT(Lpk)的实现方式已在步骤3和4中阐述。
进一步地,步骤(8)中,采用对角Hessian算子Hdiag对步骤(7)中更新后的多参数像xk+1进行预条件处理并替换,
xk+1=(Hdiag)-1xk+1
对角Hessian计算式为
Figure BDA0002245813590000053
其中,
Figure BDA0002245813590000061
Figure BDA0002245813590000062
Figure BDA0002245813590000063
将预条件后的多参数像xk+1返回步骤(5),并令k=k+1,进行下一次迭代;
进一步地,在步骤(9)中,纵波阻抗Ip、横波阻抗Is的反射率计算公式如下,
Figure BDA0002245813590000064
本发明基于声-弹耦合方程的弹性波最小二乘逆时偏移方法,首先,输入初始纵横波速度和密度初始模型,并利用声-弹耦合传播方程及其伴随方程,合成炮端下行数据和检端上行数据;然后,根据纵横波速度和密度的成像条件进行成像,并利用反偏移合成反射波数据;再次,根据四分量波形残差目标函数衡量成像结果,若不满足则利用对角Hessian预条件的CGLS算法,迭代更新纵横波速度和密度成像结果,直至满足条件,返回最终结果。
本发明相比现有技术具有如下有益效果:
本发明将LSRTM方法引入声-弹耦合方程框架中,发展一套新的基于声-弹耦合方程的弹性波最小二乘逆时偏移方法(AEC-LSRTM),该方法可以同时匹配海底观测得到的位移分量和压力分量,为四分量地震数据处理提供必要的方法,从而准确、快速地获取多种高分辨率的成像结果。
本发明利用对角Hessian预条件算子提高共轭梯度迭代的收敛效率,最终获取高分辨率的纵、横波速度和密度成像结果,实现弹性阻抗反射率的定量估计,为储层描述提供保幅性更好的成像剖面。
本发明建立了四分量地震数据的目标函数,并通过尺度因子和权重系数调控不同分量间的有效利用,实现了海底四分量地震数据的同时反传,因而相较于弹性波LSRTM方法并不需要额外的计算代价。
本发明通过最小二乘共轭梯度算法实现最小二乘成逆时偏移,压制目前四分量ERTM像中存在的背向散射噪音和耦合串扰噪音,获取提高成像精度,更准确地获取弹性阻抗反射率的定量估计。
另外,本发明所提出的方法,不仅可以用于处理海底观测的四分量地震数据,而且可以适用于海上多种混合勘探数据的联合成像,为海洋油气资源采集、勘探提供帮助,对推动多分量地震勘探的发展具有现实意义。
附图说明
图1为AEC-LSRTM方法流程图(也作摘要附图);
图2为Marmousi2局部纵(a)横波速度(b)和密度(c)反射率模型;
图3为海底4C地震观测数据(a)压力分量、(b)水平分量、(c)垂直分量;
图4为声-弹耦合传播方程模拟得到0.6s时刻正传波场压力分量快照水平分量(a)垂直分量(b)和压力分量(c);
图5为纵(a)、横(b)波速度和密度(c)RTM像;
图6为最终像预测的数据和观测数据之间的残差(a)压力分量、(b)水平分量、(c)垂直分量;
图7为纵(a)、横(b)波速度和密度(c)LSRTM像;
图8为纵(a)、横(b)波阻抗像;
图9为纵(a)、横(b)波阻抗反射率模型抽线,实线代表真实模型,虚线代表AEC-LSRTM反演结果;
图10为收敛曲线;
图11为层状模型实验中纵(a)、横(b)波速度和密度(c)的真实模型;
图12为AEC-LSRTM方法得到的纵(a)、横(b)波速度和密度(c)成像结果;
图13为ELSRTM方法得到的纵(a)、横(b)波速度和密度(c)成像结果;
图14为目标函数收敛曲线;
图15为Sigbee2A模型实验中仅利用海上密集拖缆压力数据和海底稀疏OBN的混合观测数据得到的纵(a)、横(b)波速度和密度(c)LSRTM像;
图16为Sigbee2A模型实验中仅利用海底稀疏OBN数据得到的纵(a)、横(b)波速度和密度(c)LSRTM像;
图17为OBN观测方式和混合观测方式预测的压力分量;
图18为OBN观测方式和混合观测方式预测的水平和垂直分量;
图19为拖缆观测方式、OBN观测方式和混合观测方式收敛曲线。
具体实施方式
下面结合附图和具体实施方式对本发明作进一步说明。
图1,展示了本发明所提出方法的流程图,下面以Marmousi2模型为例,介绍本发明的流程:
(1)输入纵、横波速度和密度的偏移模型。在该实例中,真实反射率模型如图2所示,模型网格设置为601x201,空间采样间隔为6.25m。该实验在海底均匀布设了61炮,每炮601道,采用主频为20Hz的Ricker子波,激发纵波源,时间记录1.6s,采样间隔为0.4ms,观测的4C数据如图3所示;
(2)根据二维声-弹耦合方程计算炮端正传波场,采用空间8阶,时间2阶的有限差分算法,获取下行波场
Figure BDA0002245813590000081
其快照如图4所示,并以此确定尺度因子ζ为5.68e-8
(3)输入多分量观测的反射数据,计算多分量伴随震源,并根据声-弹耦合伴随方程,计算检端反传波场
Figure BDA0002245813590000082
(4)根据炮端正传波场和检端反传波场以及成像条件,计算多参数像,如图5所示。图中可见,模型的主要界面被正确成像,但是在速度像中振幅随深度衰减严重,且浅层存在较强的背向散射,而密度像质量相对较好,另外多参数像受参数耦合效应的影响,在结构不一致区域均出现了串扰脚印;
(5)根据炮端正传波场、多参数像以及声-弹耦合方程的一阶Born散射方程计算反偏移波场,提取多分量合成数据,图6展示了第30代预测反射波数据和观测数据的残差,可见经过30次预条件CGLS算法的迭代更新,观测数据已经得到了非常好的匹配;
(6)根据4C数据最小二乘偏移的目标泛函求取当前目标函数,若目标函数小于阀值(理论模型阀值可以设为1%),则输出当前多参数像作为最终成像结果,进入步骤(9);若不满足,进入步骤(7),本实验总共迭代了30次;
(7)根据最小二乘共轭梯度算法进行迭代,计算更新后的多参数像;
(8)根据对角Hessian算子,计算更新后的预条件多参数像,并返回(5)。最终的速度和密度像如图7所示,与图5相比,利用最小二乘共轭梯度迭代后的多参数像得到了明显的改善,背向散射和耦合脚印均得到了明显的压制,同时不同深度界面的振幅更加均衡,对局部精细构造有更好的分辨能力;
(9)根据输出的最终多参数像,计算纵、横波阻抗的反射率,如图8所示。纵、横波阻抗像显示出最好的成像质量,可以更准确地模型深部的刻画,压制偏移中存在的假象。其中,在水平某位置抽取了垂向剖面,可以看出反演得到的像可以很好地匹配真实模型。图10展示了方法收敛曲线,可以看到AEC-LSRTM可以快速、稳健地达到收敛。这些证据都说明了AEC-LSRTM方法在处理4C地震数据,解决常规RTM方法中存在的振幅不均衡、偏移噪音以及耦合效应等问题中的有效性。
为了进一步说明AEC-LSRTM方法的优势,本申请在简单层状模型(如图11所示)上,将其和ELSRTM方法进行了对比。两种方法得到的纵、横波速度和密度像分别展示在图12和13,其收敛曲线展示在图14中。从图中可以看出,AEC-LSRTM方法具有更高的成像分辨率(如图中圆圈所示),更准确地刻画模型精细构造(如箭头所示),同时拥有更快的收敛效率。
另外,本申请提出的AEC-LSRTM方法不仅可以用于处理海底4C地震数据,同时适用于多种复杂的勘探任务中,这是目前常规声波、弹性波LSRTM方法所不能完成的。例如,一般海底OBN观测方式下,采样相对稀疏,AEC-LSRTM方法可以利用海上拖缆观测的密集数据,对稀疏的OBN数据进行补充,建立一种立体的多分量混合观测方式。
为了说明混合观测方式的有效性,将AEC-LSRTM方法分别应用于海底OBN观测数据和混合观测数据中(如图15和16所示)。可以看出,结合海上拖缆数据的混合观测方式得到的像可以明显改善由于海底观测数据过于稀疏导致不良的效果(如图15和16所示),同时更加准确地预测多分量的地震数据(如图17和18所示),更快速地进行收敛(如图19所示)。从以上实验可以看出,利用混合观测方式可以更有效地适应海洋地震勘探采集技术的发展和需求,有针对性地设计的海上、海底立体观测方式,从而满足不同的海洋勘探需求,对推动海底地震勘探方法和技术的发展具有现实意义。

Claims (9)

1.一种基于声-弹耦合方程的弹性波最小二乘逆时偏移方法,其特征在于,具体包括以下步骤:
(1)输入纵波速度、横波速度和密度的偏移模型;
(2)根据二维声-弹耦合方程计算炮端正传波场;
(3)输入多分量地震观测的反射数据,计算多分量伴随震源,并根据声-弹耦合伴随方程计算检端反传波场;
(4)根据炮端正传波场和检端反传波场以及成像条件,计算多参数像;
(5)根据炮端正传波场、多参数像以及声-弹耦合方程的一阶Born散射方程计算反偏移波场,提取多分量合成数据;
(6)根据4C数据最小二乘偏移的目标泛函求取当前目标函数,若目标函数小于阀值,则输出当前多参数像作为最终成像结果,进入步骤(9);若不满足,进入步骤(7);
其中,4C数据最小二乘偏移目标泛函E(m)定义为:
Figure FDA0002436478810000011
其中,m代表模型参数,du和dP分别表示位移分量和压力分量,上标syn和obs分别代表合成数据和观测数据,ε是权重系数,ε∈[0,1],用于调节压力数据和位移数据;ζ是尺度因子;
(7)根据最小二乘共轭梯度迭代,计算更新后的多参数像;
(8)根据对角Hessian算子,计算更新后的预条件多参数像,并返回(5);
(9)根据输出的最终多参数像,计算纵、横波阻抗的反射率。
2.根据权利要求1所述基于声-弹耦合方程的弹性波最小二乘逆时偏移方法,其特征在于,步骤(2)中,二维声-弹耦合方程如下:
Figure FDA0002436478810000021
Figure FDA0002436478810000022
Figure FDA0002436478810000023
Figure FDA0002436478810000024
Figure FDA0002436478810000025
其中,正传弹性波场包括5个分量
Figure FDA0002436478810000026
依次代表:水平和垂直位移分量、压力分量以及与横波相关的正应力分量和切应力分量,fP代表正传时施加在压力分量上的子波,λ和μ为拉梅常数,ρ为密度;根据该方程计算炮端正传弹性波场。
3.根据权利要求2所述基于声-弹耦合方程的弹性波最小二乘逆时偏移方法,其特征在于,步骤(3)中,输入多分量地震观测数据,包括压力分量
Figure FDA0002436478810000027
水平分量
Figure FDA0002436478810000028
和垂直分量
Figure FDA0002436478810000029
计算多分量伴随震源压力分量f’p,水平分量f’x和垂直分量f’z,其计算公式如下,
Figure FDA00024364788100000210
并根据声-弹耦合伴随方程,计算检端反传波场
Figure FDA00024364788100000211
Figure FDA00024364788100000212
Figure FDA00024364788100000213
Figure FDA00024364788100000214
Figure FDA00024364788100000215
其中,反传波场也包括的5个分量
Figure FDA00024364788100000216
依次代表:水平和垂直位移分量、压力分量以及与横波相关的正应力分量和切应力分量。
4.根据权利要求3所述基于声-弹耦合方程的弹性波最小二乘逆时偏移方法,其特征在于,步骤(4)中,根据炮端正传波场和检端反传波场以及成像条件,计算多参数像,其中拉梅常数和密度的成像公式为:
Figure FDA0002436478810000031
Figure FDA0002436478810000032
Figure FDA0002436478810000033
根据链式法则,纵、横波速度和密度的成像公式如下,
Figure FDA0002436478810000034
5.根据权利要求4所述基于声-弹耦合方程的弹性波最小二乘逆时偏移方法,其特征在于,步骤(5)中,声-弹耦合方程的一阶Born散射方程为:
Figure FDA0002436478810000035
Figure FDA0002436478810000036
Figure FDA0002436478810000037
Figure FDA0002436478810000038
Figure FDA0002436478810000039
其中,
Figure FDA00024364788100000310
为反偏移波场的5个分量;
并根据炮端正传波场和多参数像,计算反偏移波场,在检波点处抽取合成的反射波数据压力分量
Figure FDA00024364788100000311
和位移分量
Figure FDA00024364788100000312
6.根据权利要求1-5任一所述基于声-弹耦合方程的弹性波最小二乘逆时偏移方法,其特征在于,步骤(6)中,利用步骤(2)中获取的位移分量ux、uz和压力分量P进行计算,公式如下:
Figure FDA00024364788100000313
7.根据权利要求6所述基于声-弹耦合方程的弹性波最小二乘逆时偏移方法,其特征在于,步骤(7)中,最小二乘共轭梯度CGLS算法如下:
下列变量中,xk代表当前多参数像,xk+1代表更新后(即下一代)多参数像,其余变量均为数学算法的中间变量不具备明确的物理意义;
在第1次迭代前,需设定初值x0=0;r0=Hx0-mmig,p0=-r0,k=0;
然后按下列公式求解更新后的成像结果xk+1
Figure FDA0002436478810000041
xk+1=xkkpk
rk+1=rkkHpk=rkk[LT(Lpk)]
Figure FDA0002436478810000042
pk+1=-rk+1k+1pk
8.根据权利要求7所述基于声-弹耦合方程的弹性波最小二乘逆时偏移方法,其特征在于,步骤(8)中,采用对角Hessian算子Hdiag对步骤(7)中更新后的多参数像xk+1进行预条件处理并替换,
xk+1=(Hdiag)-1xk+1
对角Hessian计算式为
Figure FDA0002436478810000043
其中,
Figure FDA0002436478810000044
Figure FDA0002436478810000045
Figure FDA0002436478810000046
将预条件后的多参数像xk+1返回步骤(5),并令k=k+1,进行下一次迭代。
9.根据权利要求8所述基于声-弹耦合方程的弹性波最小二乘逆时偏移方法,其特征在于,在步骤(9)中,纵波阻抗Ip、横波阻抗Is的反射率计算公式如下:
Figure FDA0002436478810000051
CN201911016294.6A 2019-10-24 2019-10-24 一种基于声-弹耦合方程的弹性波最小二乘逆时偏移方法 Active CN110687600B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911016294.6A CN110687600B (zh) 2019-10-24 2019-10-24 一种基于声-弹耦合方程的弹性波最小二乘逆时偏移方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911016294.6A CN110687600B (zh) 2019-10-24 2019-10-24 一种基于声-弹耦合方程的弹性波最小二乘逆时偏移方法

Publications (2)

Publication Number Publication Date
CN110687600A CN110687600A (zh) 2020-01-14
CN110687600B true CN110687600B (zh) 2020-05-22

Family

ID=69114385

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911016294.6A Active CN110687600B (zh) 2019-10-24 2019-10-24 一种基于声-弹耦合方程的弹性波最小二乘逆时偏移方法

Country Status (1)

Country Link
CN (1) CN110687600B (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111948710A (zh) * 2020-08-05 2020-11-17 河海大学 一种利用海底多分量地震记录进行弹性矢量波成像的方法
CN111881612B (zh) * 2020-08-05 2021-03-19 武汉市政工程设计研究院有限责任公司 一种正应力和剪应力不同权重二维应力场反演方法及系统
CN112213783B (zh) * 2020-09-16 2021-11-12 河海大学 一种利用海底多分量地震记录进行各向异性地震成像方法
US11733413B2 (en) 2021-04-30 2023-08-22 Saudi Arabian Oil Company Method and system for super resolution least-squares reverse time migration

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US10234581B2 (en) * 2015-06-19 2019-03-19 Chevron U.S.A. Inc. System and method for high resolution seismic imaging
CN107589443B (zh) * 2017-08-16 2019-05-28 东北石油大学 基于弹性波最小二乘逆时偏移成像的方法及系统
WO2019182564A1 (en) * 2018-03-20 2019-09-26 Equinor Us Operations Llc. Sub-surface structure determination using seismic migration
CN110133720B (zh) * 2019-06-04 2020-02-18 南京信息工程大学 一种基于统计岩石物理模型的横波速度预测方法
CN110135112A (zh) * 2019-06-04 2019-08-16 南京信息工程大学 一种基于粒子滤波算法的页岩各向异性估计方法

Also Published As

Publication number Publication date
CN110687600A (zh) 2020-01-14

Similar Documents

Publication Publication Date Title
CN110687600B (zh) 一种基于声-弹耦合方程的弹性波最小二乘逆时偏移方法
US10802167B2 (en) Seismic acquisition method and apparatus
US6654693B2 (en) Angle dependent surface multiple attenuation for two-component marine bottom sensor data
EP2601542B1 (en) Method for separating independent simultaneous sources
MX2007002733A (es) Sistema para la atenuacion de las multiples de los fondos de agua en los datos sismicos registrados mediante sensores de presion y sensores del movimiento de las particulas.
CN110058303B (zh) 声波各向异性逆时偏移混合方法
EP2730949A2 (en) Interference noise attenuation method and apparatus
NO326400B1 (no) Fremgangsmate for dempning av overflate-multipler i multidimensjonale marinseismiske data
GB2375606A (en) Angle dependent surface multiple attenuation for two - component marine bottom sensor data
CN111175822A (zh) 改进直接包络反演与扰动分解的强散射介质反演方法
Breistøl Broadband processing of conventional 3D seismic data for near surface geohazard investigation: A North Sea case study
CN116088055A (zh) 用于海洋obn数据的上下行波场分离的方法
CN116819624A (zh) 基于波场分解的一次波和多次波分离及同时成像方法
Campman et al. A 3D Imaging Approach to Short-Wavelength Statics
Przebindowska et al. 2D acoustic full waveform tomography of marine streamer data–data preparation and choice of inversion strategies
Przebindowska et al. 2D Acoustic FWT of Marine Streamer Data Problems and Data Preprocessing-PART I

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