CN110888159B - 基于角度分解与波场分离的弹性波全波形反演方法 - Google Patents

基于角度分解与波场分离的弹性波全波形反演方法 Download PDF

Info

Publication number
CN110888159B
CN110888159B CN201911122015.4A CN201911122015A CN110888159B CN 110888159 B CN110888159 B CN 110888159B CN 201911122015 A CN201911122015 A CN 201911122015A CN 110888159 B CN110888159 B CN 110888159B
Authority
CN
China
Prior art keywords
wave field
velocity
wave
longitudinal
omega
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.)
Expired - Fee Related
Application number
CN201911122015.4A
Other languages
English (en)
Other versions
CN110888159A (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.)
Xian University of Technology
Original Assignee
Xian University of 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 Xian University of Technology filed Critical Xian University of Technology
Priority to CN201911122015.4A priority Critical patent/CN110888159B/zh
Publication of CN110888159A publication Critical patent/CN110888159A/zh
Application granted granted Critical
Publication of CN110888159B publication Critical patent/CN110888159B/zh
Expired - Fee Related 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
    • G01V1/282—Application of seismic models, synthetic seismograms
    • 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
    • G01V1/30—Analysis
    • G01V1/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity

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)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明公开了基于角度分解与波场分离的弹性波全波形反演方法，具体为：1、采集得到多分量观测数据；2、计算频率域地震数据；3、构建矩形网格地质模型；4、给定目标函数；5、给定用于反演的频率分量；6、对于第一个频率分量，对地震数据进行分解；7、给定用于反演的散射角范围；8、对目标函数进行优化；9、对地震数据进行分解得到无角度信息的纵波波场和横波波场；10、分别利用纵波波场和横波波场对目标函数进行优化，得到纵波速度和横波速度；11、对于其他频率分量，重复步骤9和步骤10，直到获得最终反演结果。本发明方法能够有效的避免弹性波全波形反演中的局部极值问题和多参数串扰问题。

Description

基于角度分解与波场分离的弹性波全波形反演方法
技术领域
本发明属于地球物理勘探技术领域，具体涉及一种基于角度分解与波场分离的弹性波全波形反演方法。
背景技术
石油和天然气是事关国家经济发展及国家安全的重要战略资源，随着石油工业的发展和勘探开发的不断深入，从构造勘探阶段逐渐走向岩性勘探阶段，勘探难度逐渐加大。为了顺应这种要求，全波形反演迅速发展起来，并成为当今地球物理界的研究热点。全波形反演是进行地层参数估计的有效方法，不但能够为地震成像提供高分辨率高精度的速度模型，而且具有同时反演各种不同参数的能力。然而全波形反演的推广应用也受到一些问题的影响。首先，全波形反演对初始模型的依赖性极强。作为一个高度非线性问题，全波形反演通常采用迭代优化算法进行求解，而迭代的顺利进行需要一个足够准确的初始模型，即大尺度背景介质参数。大尺度的背景介质与地震数据中的低频分量相对应，然而低频分量在地震信号的采集过程中往往是缺失的，这就导致全波形反演容易陷入局部极值。近年来，学者们发展了不同的方法来克服局部极值问题，如拉普拉斯-傅里叶域全波形反演方法，包络反演方法，波场积分方法等通过某种非线性变换制造出低频信息进而用于反演。层析成像全波形反演方法通过在扩展域对问题进行求解从而得到大尺度背景介质。自适应全波形反演方法将波场匹配问题转化成一个滤波器系数的优化问题来对反演问题进行求解。上述方法均在一定程度上缓解了全波形反演中的局部极值问题。
目前，全波形反演在声波领域得到了较为广泛的发展。声波全波形反演在实现上相对来说较为简单，然而声波近似只考虑了波场的纵波分量而忽略了波场的横波分量。与之相比，采用弹性波动方程能够更加准确得刻画波场的传播特征，有利于更加精确得了解地层介质信息，因此研究弹性波全波形反演具有重要意义。然而在弹性介质情况下，存在多种波场的相互耦合，导致多参数反演的相互串扰问题，为弹性全波形反演带来困难。学者们研究了不同的方法来解决弹性波反演中的参数串扰问题，如交替反演，分级反演，利用Hessian矩阵进行解耦等。
对于弹性波全波形反演，在缺失低频的情况下，如何同时有效解决局部极值问题和多参数串扰问题是提高弹性波全波形反演的关键。
发明内容
本发明的目的是提供一种基于角度分解和波场分离的弹性波全波形反演方法，该方法通过采用慢度域的局部分解变换，能够同时实现波场在角度上的分解和不同波场模式的分离，从而达到同时克服反演的局部极值问题和多参数串扰问题的效果。
本发明所采用的技术方案是，基于角度分解与波场分离的弹性波全波形反演方法，具体按照以下步骤实施：
步骤1、设置炮点和检波点，采集得到多分量观测数据
Figure BDA0002275708580000021
其中t表示时间变量；xr,xs分别表示检波点和炮点的位置；
Figure BDA0002275708580000022
为包含x分量和z分量的多分量波场向量：
Figure BDA0002275708580000023
步骤2、根据傅里叶变换计算频率域地震数据u(ω,xr；xs)，其中ω表示频率变量；u(ω,xr；xs)仍然包含两个分量：u(ω,xr；xs)＝[ux(ω,xr；xs),uz(ω,xr；xs)]；
步骤3、构建矩形网格地质模型，设定模型大小为：横向Nx个网格点，纵向Nz个网格点，并设定正演模拟的空间采样间隔、时间采样间隔dt、最大采样时间Nt，所述空间采样间隔包括横向采样间隔dx和纵向采样间隔dz；
步骤4、给定纵波速度模型vP(x)和横波速度模型vS(x)的初始纵波速度模型vP0(x)和初始横波速度模型vS0(x)，给定目标函数J(vP(x),vS(x))；
步骤5、给定用于反演的频率分量ω1,ω2,...ωN；
步骤6、对于第一个频率分量ω1即最低频率分量，对地震数据进行分解，得到频率角度域纵波波场uP(θ,x,ω)和横波波场uS(θ,x,ω)，其中θ表示角度分量；
步骤7、对于第一个频率分量ω1，给定用于反演的散射角范围；
步骤8、在所给定散射角范围内，针对纵波速度vP(x)和横波速度vS(x)，分别利用纵波波场uP(θ,x,ω)和横波波场uS(θ,x,ω)通过迭代优化算法对目标函数J(vP(x),vS(x))进行优化，得到大尺度纵波速度vP1(x)和大尺度横波速度vS1(x)；
步骤9、对于第二个频率分量ω2，对地震数据进行分解，并对分解后数据进行角度叠加，得到无角度信息的纵波波场uP(x,ω)和横波波场uS(x,ω)；
步骤10、针对纵波速度vP(x)和横波速度vS(x)，分别利用纵波波场uP(x,ω)和横波波场uS(x,ω)通过迭代优化算法对目标函数J(vP(x),vS(x))进行优化，得到纵波速度vP2(x)和横波速度vS2(x)；
步骤11、对于其他频率分量，重复步骤9和步骤10，直到获得纵波速度vPN(x)和横波速度vSN(x)，即最终反演结果。
本发明的特点还在于，
步骤2中，频率域地震数据u(ω,xr；xs)根据如下傅里叶变换公式计算得到：
Figure BDA0002275708580000041
式(1)中，T为波场最大观测时间。
步骤4中，设定目标函数J(vP(x),vS(x))为观测波场与计算波场之间残差的二范数，具有如下形式：
Figure BDA0002275708580000042
式(2)中，uobs(ω)和ucal(ω)分别表示频率域观测波场和频率域计算波场，其中频率域计算波场通过如下方式得到：
首先根据时间域弹性波动方程得到时间域计算波场
Figure BDA0002275708580000043
再对
Figure BDA0002275708580000044
进行步骤2中所述傅里叶变换得到ucal(ω)；上标t表示转置，*表示取共轭。
步骤4中，所述时间域计算波场
Figure BDA0002275708580000045
满足如下弹性波动方程：
Figure BDA0002275708580000046
式(3)中，δ(x-xs)fx(t)和δ(x-xs)fz(t)表示位于xs处的波场x分量和z分量震源函数；ρ表示地层密度，设置为常数；λ和μ表示拉梅系数，与纵波速度vP和横波速度vS之间存在如下关系：
Figure BDA0002275708580000051
步骤6中，对最低频率数据进行分解，得到频率角度域纵波波场uP(θ,x,ω)和横波波场uS(θ,x,ω)的具体公式如下：
Figure BDA0002275708580000052
Figure BDA0002275708580000053
式(5)中，uP(n,x,ω)和uS(n,x,ω)为分解得到的慢度域纵波波场和横波波场；
Figure BDA0002275708580000054
表示慢度向量，其中
Figure BDA0002275708580000055
为指向波场传播方向的单位向量，n为慢度向量n的绝对值；W(x′-x)表示以坐标点x为中心的采样窗；“×”和“·”分别表示叉乘运算和点乘运算；
式(6)中，nP和nS表示纵波慢度的绝对值和横波慢度的绝对值，即
Figure BDA0002275708580000056
其中
Figure BDA0002275708580000057
和
Figure BDA0002275708580000058
表示采样窗内的纵波平均速度和横波平均速度。结合式(5)和式(6)可以得到角度域纵波波场uP(θ,x,ω)和横波波场uS(θ,x,ω)。
步骤6具体实施方式如下：
根据式(5)和式(6)分别对震源侧波场
Figure BDA0002275708580000059
和接收点侧波场
Figure BDA00022757085800000510
进行分解，其中两个上标箭头为区分震源侧波场和接收点侧波场，得到震源侧角度域纵波波场
Figure BDA00022757085800000511
和横波波场
Figure BDA00022757085800000512
以及接收点侧角度域纵波波场
Figure BDA00022757085800000513
和横波波场
Figure BDA00022757085800000514
其中θi和θg分别表示震源侧波场传播方向和接收点侧波场传播方向，波场散射角通过如下公式计算得到：
Θ＝|θg-θi| (7)
式(7)中Θ表示波场散射角。
步骤7中，选取散射角范围依次为0°～30°，0°～60°和0°～180°。
步骤8中，分别利用纵波波场uP(θ,x,ω)和横波波场uS(θ,x,ω)通过迭代优化算法对目标函数J(vP(x),vS(x))进行优化的具体实施方式如下：
在弹性波传播过程中，纵波和横波经过散射会产生转换波，整个波场中存在多种模式波，即PP波，PS波，SP波，SS波；为了避免在反演过程中vP和vS的相互串扰问题，反演vP时只使用PP波，而反演vS时只使用SS、SP和PS波；为了避免局部极值问题，从小散射角到大散射角依次进行反演；因此在考虑角度分解与波场分离的情况下，vP和vS的梯度公式如下：
Figure BDA0002275708580000061
Figure BDA0002275708580000062
Figure BDA0002275708580000063
Figure BDA0002275708580000071
式(8)～(11)中，
Figure BDA0002275708580000072
和
Figure BDA0002275708580000073
分别对应于纵波速度vP(x)和横波速度vS(x)的角度滤波器，控制散射角在选定的范围之内；“Re”表示取实部；计算出梯度之后，利用如下准则对纵波速度vP(x)和横波速度vS(x)进行迭代更新：
Figure BDA0002275708580000074
式(12)中，vP(x)n,vP(x)n+1分别为第n步和第n+1步迭代的纵波速度值，vS(x)n,vS(x)n+1分别为第n步和第n+1步迭代的横波速度值；αn为迭代步长，为一个正常数；
通过上述方法计算速度梯度并进行迭代，在迭代收敛之后得到大尺度纵波速度vP1(x)和大尺度横波速度vS1(x)。
步骤9中，对分解后数据进行角度叠加，得到无角度信息的纵波波场和横波波场的公式如下：
Figure BDA0002275708580000075
Figure BDA0002275708580000076
式(13)和式(14)中，
Figure BDA0002275708580000077
和
Figure BDA0002275708580000078
表示震源侧无角度信息的纵波波场和横波波场；
Figure BDA0002275708580000079
和
Figure BDA00022757085800000710
表示接收点侧无角度信息的纵波波场和横波波场；
Figure BDA00022757085800000711
和
Figure BDA00022757085800000712
分别表示对震源侧波场传播角度求和以及对接收点侧波场传播方向求和。
步骤10中，分别利用纵波波场uP(x,ω)和横波波场uS(x,ω)通过迭代优化算法对目标函数J(vP(x),vS(x))进行优化的具体实施方式如下：在仅考虑波场分离的情况下，vP和vS的梯度公式如下：
Figure BDA0002275708580000081
Figure BDA0002275708580000082
Figure BDA0002275708580000083
Figure BDA0002275708580000084
计算出梯度之后，利用与式(12)相同的准则对纵波速度vP(x)和横波速度vS(x)进行迭代更新；通过上述方法计算速度梯度并进行迭代，在迭代收敛之后得到第二个频率分量下的纵波速度vP2(x)和大尺度横波速度vS2(x)。
本发明一种基于角度分解与波场分离的弹性波全波形反演方法的有益效果是：首先，从线性初始模型出发，对于地震数据中的最低频率分量进行分解变换，得到不同角度信息下的纵波波场与横波波场；其次，基于分解后的波场得到包含不同散射角信息以及不同波场模式的速度梯度，并得到大尺度纵波速度与大尺度横波速度；再次，对于从第二个频率分量开始的其他频率分量，利用不含有角度信息的纵波波场和横波波场得到不含有角度信息的不同波场模式的速度梯度，并对速度进行迭代更新；最后，遍历所有频率分量直到得到最终反演的纵波速度和横波速度。
与常用的反演方法相比，该方法能够有效的避免弹性波全波形反演中的局部极值问题和多参数串扰问题。该方法具有上述优势的原因在于：首先，因为在反演过程中考虑到散射角信息，利用散射角与不同尺度的介质扰动之间的关系，避免了反演中的局部极值问题；其次，因为对波场进行了分离，针对纵波速度和横波速度的散射特性对这两者进行解耦，所以能够避免这两者之间的相互串扰。
附图说明
图1是本发明方法的流程图；
图2是目标区真实纵波速度模型；
图3是目标区真实横波速度模型；
图4是线性初始纵波速度模型；
图5是线性初始横波速度模型；
图6是使用传统方法反演得到的纵波速度模型；
图7是使用传统方法反演得到的横波速度模型；
图8是使用本发明方法反演得到的大尺度纵波速度模型；
图9是使用本发明方法反演得到的大尺度横波速度模型；
图10是使用本发明方法反演得到的最终纵波速度模型；
图11是使用本发明方法反演得到的最终横波速度模型。
具体实施方式
下面结合附图和具体实施方式对本发明进行详细说明。
本发明基于角度分解与波场分离的弹性波全波形反演方法，如图1所示，具体步骤分别如下：
步骤1、设置炮点和检波点，对原始地震数据进行采集得到多分量观测数据
Figure BDA0002275708580000101
并对采集到的地震数据进行常规预处理，包括动静校正、噪声消除等处理，其中t表示时间变量；xr,xs分别表示检波点和炮点的位置；
Figure BDA0002275708580000102
为包含x分量和z分量的多分量波场向量：
Figure BDA0002275708580000103
步骤2、根据傅里叶变换计算频率域地震数据u(ω,xr；xs)，其中ω表示频率变量；u(ω,xr；xs)仍然包含两个分量：u(ω,xr；xs)＝[ux(ω,xr；xs),uz(ω,xr；xs)]；频率域地震数据根u(ω,xr；xs)据如下傅里叶变换公式计算得到：
Figure BDA0002275708580000104
式(1)中，T为波场最大观测时间。
步骤3、构建矩形网格地质模型，设定模型大小为：横向Nx个网格点，纵向Nz个网格点，并设定正演模拟的空间采样间隔(包括横向采样间隔dx和纵向采样间隔dz)、时间采样间隔dt、最大采样时间Nt。实际中可以根据地震数据有效频带范围、地震记录时间长度、地震数据观测系统确定相关参数；选定网格参数的标准是使得基于该网格进行有限差分正演模拟时，可以不仅满足稳定性条件而且有效压制数值频散。
步骤4、给定纵波速度模型vP(x)和横波速度模型vS(x)的初始纵波速度模型vP0(x)和初始横波速度模型vS0(x)，给定目标函数J(vP(x),vS(x))；设定目标函数J(vP(x),vS(x))为观测波场与计算波场之间残差的二范数，具有如下形式：
Figure BDA0002275708580000105
式(2)中，uobs(ω)和ucal(ω)分别表示频率域观测波场和频率域计算波场，其中频率域计算波场通过如下方式得到：
首先根据时间域弹性波动方程得到时间域计算波场
Figure BDA0002275708580000111
再对
Figure BDA0002275708580000112
进行步骤2中所述傅里叶变换得到ucal(ω)；上标t表示转置，*表示取共轭。
步骤4中，时间域计算波场
Figure BDA0002275708580000113
满足如下弹性波动方程：
Figure BDA0002275708580000114
式(3)中，δ(x-xs)fx(t)和δ(x-xs)fz(t)表示位于xs处的波场x分量和z分量震源函数；ρ表示地层密度，设置为常数；λ和μ表示拉梅系数，与纵波速度vP和横波速度vS之间存在如下关系：
Figure BDA0002275708580000115
步骤5、给定用于反演的频率分量ω1,ω2,...ωN，N取7，ω1,ω2,...ω7分别对应设置为5Hz,8Hz,12Hz,16Hz,20Hz,25Hz,30Hz；
步骤6、对于第一个频率分量ω1(即最低频率分量)，对地震数据进行分解，得到频率角度域纵波波场uP(θ,x,ω)和横波波场uS(θ,x,ω)，其中θ表示角度分量；
对最低频率分量5Hz数据进行分解，得到频率角度域纵波波场uP(θ,x,ω)和横波波场uS(θ,x,ω)的具体公式如下：
Figure BDA0002275708580000121
Figure BDA0002275708580000122
式(5)中，uP(n,x,ω)和uS(n,x,ω)为分解得到的慢度域纵波波场和横波波场。
Figure BDA0002275708580000123
表示慢度向量，其中
Figure BDA0002275708580000124
为指向波场传播方向的单位向量，n为慢度向量n的绝对值。W(x′-x)表示以坐标点x为中心的采样窗。“×”和“·”分别表示叉乘运算和点乘运算。
式(6)中，nP和nS表示纵波慢度的绝对值和横波慢度的绝对值，即
Figure BDA0002275708580000125
其中
Figure BDA0002275708580000126
和
Figure BDA0002275708580000127
表示采样窗内的纵波平均速度和横波平均速度。结合式(5)和式(6)可以得到角度域纵波波场uP(θ,x,ω)和横波波场uS(θ,x,ω)。
步骤6具体实施方式如下：
根据式(5)和式(6)分别对震源侧波场
Figure BDA0002275708580000128
和接收点侧波场
Figure BDA0002275708580000129
进行分解(其中两个上标箭头为区分震源侧波场和接收点侧波场)，得到震源侧角度域纵波波场
Figure BDA00022757085800001210
和横波波场
Figure BDA00022757085800001211
以及接收点侧角度域纵波波场
Figure BDA00022757085800001212
和横波波场
Figure BDA00022757085800001213
其中θi和θg分别表示震源侧波场传播方向和接收点侧波场传播方向。波场散射角可以通过如下公式计算得到：
Θ＝|θg-θi| (7)
式(7)中Θ表示波场散射角。
步骤7、对于第一个频率分量ω1，给定用于反演的散射角范围；选取散射角范围依次为0°～30°，0°～60°和0°～180°；
步骤8、在所给定散射角范围内，针对纵波速度vP(x)和横波速度vS(x)，分别利用纵波波场uP(θ,x,ω)和横波波场uS(θ,x,ω)通过迭代优化算法对目标函数J(vP(x),vS(x))进行优化，得到大尺度纵波速度vP1(x)和大尺度横波速度vS1(x)，具体实施方式如下：
在弹性波传播过程中，纵波和横波经过散射会产生转换波，整个波场中存在多种模式波，即PP波(纵波-纵波)，PS波(纵波-横波)，SP波(横波-纵波)，SS波(横波-横波)。不同模式的波对不同介质参数敏感程度不同，即经过不同参数散射会产生不同模式的波。纵波速度vP主要产生PP波，而横波速度vS产生SS波、SP波、PS波以及PP波，其中SS波占主要地位。为了避免在反演过程中vP和vS的相互串扰问题，反演vP时只使用PP波，而反演vS时只使用SS、SP和PS波。为了避免局部极值问题，从小散射角到大散射角依次进行反演。因此在考虑角度分解与波场分离的情况下，vP和vS的梯度公式如下：
Figure BDA0002275708580000131
Figure BDA0002275708580000132
Figure BDA0002275708580000141
Figure BDA0002275708580000142
式(8)～(11)中，
Figure BDA0002275708580000143
和
Figure BDA0002275708580000144
分别对应于纵波速度vP(x)和横波速度vS(x)的角度滤波器，控制散射角在选定的范围之内；“Re”表示取实部。计算出梯度之后，利用如下准则对纵波速度vP(x)和横波速度vS(x)进行迭代更新：
Figure BDA0002275708580000145
式(12)中，vP(x)n,vP(x)n+1分别为第n步和第n+1步迭代的纵波速度值，vS(x)n,vS(x)n+1分别为第n步和第n+1步迭代的横波速度值；αn为迭代步长，为一个正常数；
通过上述方法计算速度梯度并进行迭代，在迭代收敛之后可以得到大尺度纵波速度vP1(x)和大尺度横波速度vS1(x)；
步骤9、对于第二个频率分量ω2，对地震数据进行分解，并对分解后数据进行角度叠加，得到无角度信息的纵波波场uP(x,ω)和横波波场uS(x,ω)；
由于已经得到大尺度纵波速度vP1(x)和大尺度横波速度vS1(x)，有效避免了局部极值问题，后续各频率分量的反演中不需要再分角度反演，只需考虑不同模式波场的分离。对分解后数据进行角度叠加，得到无角度信息的纵波波场和横波波场的公式如下：
Figure BDA0002275708580000151
Figure BDA0002275708580000152
式(13)和式(14)中，
Figure BDA0002275708580000153
和
Figure BDA0002275708580000154
表示震源侧无角度信息的纵波波场和横波波场；
Figure BDA0002275708580000155
和
Figure BDA0002275708580000156
表示接收点侧无角度信息的纵波波场和横波波场；
Figure BDA0002275708580000157
和
Figure BDA0002275708580000158
分别表示对震源侧波场传播角度求和以及对接收点侧波场传播方向求和；
步骤10、针对纵波速度vP(x)和横波速度vS(x)，分别利用纵波波场uP(x,ω)和横波波场uS(x,ω)通过迭代优化算法对目标函数J(vP(x),vS(x))进行优化，得到纵波速度vP2(x)和横波速度vS2(x)，具体实施方式如下：
在仅考虑波场分离的情况下，vP和vS的梯度公式如下：
Figure BDA0002275708580000159
Figure BDA00022757085800001510
Figure BDA00022757085800001511
Figure BDA0002275708580000161
计算出梯度之后，利用与式(12)相同的准则对纵波速度vP(x)和横波速度vS(x)进行迭代更新；
通过上述方法计算速度梯度并进行迭代，在迭代收敛之后可以得到第二个频率分量下的纵波速度vP2(x)和大尺度横波速度vS2(x)。
步骤11、对于其他频率分量，重复步骤9和步骤10，直到获得纵波速度vPN(x)和横波速度vSN(x)，即最终反演结果。
模型算例
本发明的具体实施过程被应用于某弹性地质模型。该模型横向共有9200米，纵向共有3000米。目标区的纵波速度模型和横波速度模型分别如图2和3所示。设置炮点个数为23，检波点个数为230。炮点和检波点均等间隔分布于目标区顶端。将目标区在横向和纵向上进行等间隔网格化分，横向划分为2300个网格，纵向划分为750个网格，横向采样间隔和纵向采样间隔均为4m。记录波场数据的时间采样间隔为1ms。
纵波速度模型和横波速度模型的初始模型设定为线性初始模型，分别如图4和5所示，其中纵波速度模型从浅到深逐渐从1500m/s线性增长到4500m/s，横波速度模型从浅到深逐渐从400m/s线性增长到2500m/s。
采用传统反演方法得到的纵波速度模型和横波速度模型如图6和图7所示。可以看到传统反演方法得到的速度模型受到局部极值问题和多参数串扰问题的影响。
经过本发明所提出的方法的步骤1到步骤8进行反演得到的目标区大尺度纵波速度模型和大尺度横波速度模型如图8和图9所示。步骤5中，N取7，ω1,ω2,...ω7分别对应设置为5Hz,8Hz,12Hz,16Hz,20Hz,25Hz,30Hz；步骤7中，对于第一个频率分量ω1，给定用于反演的散射角范围；选取散射角范围依次为0°～30°，0°～60°和0°～180°；
大尺度的速度模型即长波长速度分量，在反演过程中，长波长分量的准确获取是整个参数反演的关键，因为它保证了迭代能够收敛到正确的全局最优解以及后续精细结构的进一步反演。
基于上述所获得的大尺度速度模型，使用本发明中所提出的方法反演得到的精细纵波速度模型和横波速度模型如图10和图11所示。可以看出所得纵波速度和横波速度都与真实模型接近，说明了本发明中所提出方法的有效性。

Claims (9)

1.基于角度分解与波场分离的弹性波全波形反演方法，其特征在于，具体按照以下步骤实施：
步骤1、设置炮点和检波点，采集得到多分量观测数据
Figure FDA0003021427830000011
其中t表示时间变量；xr,xs分别表示检波点和炮点的位置；
Figure FDA0003021427830000012
为包含x分量和z分量的多分量波场向量：
Figure FDA0003021427830000013
步骤2、根据傅里叶变换计算频率域地震数据u(ω,xr；xs)，其中ω表示频率变量；u(ω,xr；xs)仍然包含两个分量：u(ω,xr；xs)＝[ux(ω,xr；xs),uz(ω,xr；xs)]；
步骤3、构建矩形网格地质模型，设定模型大小为：横向Nx个网格点，纵向Nz个网格点，并设定正演模拟的空间采样间隔、时间采样间隔dt、最大采样时间Nt，所述空间采样间隔包括横向采样间隔dx和纵向采样间隔dz；
步骤4、给定纵波速度模型vP(x)和横波速度模型vS(x)的初始纵波速度模型vP0(x)和初始横波速度模型vS0(x)，给定目标函数J(vP(x),vS(x))；
步骤5、给定用于反演的频率分量ω1,ω2,...ωN；
步骤6、对于第一个频率分量ω1即最低频率分量，对地震数据进行分解，得到频率角度域纵波波场uP(θ,x,ω)和横波波场uS(θ,x,ω)，其中θ表示角度分量；
步骤7、对于第一个频率分量ω1，给定用于反演的散射角范围；
步骤8、在所给定散射角范围内，针对纵波速度vP(x)和横波速度vS(x)，分别利用纵波波场uP(θ,x,ω)和横波波场uS(θ,x,ω)通过迭代优化算法对目标函数J(vP(x),vS(x))进行优化，得到大尺度纵波速度vP1(x)和大尺度横波速度vS1(x)；
步骤9、对于第二个频率分量ω2，对地震数据进行分解，并对分解后数据进行角度叠加，得到无角度信息的纵波波场uP(x,ω)和横波波场uS(x,ω)；
步骤10、针对纵波速度vP(x)和横波速度vS(x)，分别利用纵波波场uP(x,ω)和横波波场uS(x,ω)通过迭代优化算法对目标函数J(vP(x),vS(x))进行优化，得到纵波速度vP2(x)和横波速度vS2(x)；
步骤11、对于其他频率分量，重复步骤9和步骤10，直到获得纵波速度vPN(x)和横波速度vSN(x)，即最终反演结果。
2.根据权利要求1所述的基于角度分解与波场分离的弹性波全波形反演方法，其特征在于，步骤2中，频率域地震数据u(ω,xr；xs)根据如下傅里叶变换公式计算得到：
Figure FDA0003021427830000021
式(1)中，T为波场最大观测时间。
3.根据权利要求2所述的基于角度分解与波场分离的弹性波全波形反演方法，其特征在于，步骤4中，设定目标函数J(vP(x),vS(x))为观测波场与计算波场之间残差的二范数，具有如下形式：
Figure FDA0003021427830000022
式(2)中，uobs(ω)和ucal(ω)分别表示频率域观测波场和频率域计算波场，其中频率域计算波场通过如下方式得到：
首先根据时间域弹性波动方程得到时间域计算波场
Figure FDA0003021427830000023
再对
Figure FDA0003021427830000024
进行步骤2中所述傅里叶变换得到ucal(ω)；上标t表示转置，*表示取共轭。
4.根据权利要求3所述的基于角度分解与波场分离的弹性波全波形反演方法，其特征在于，步骤4中，所述时间域计算波场
Figure FDA0003021427830000031
满足如下弹性波动方程：
Figure FDA0003021427830000032
式(3)中，δ(x-xs)fx(t)和δ(x-xs)fz(t)表示位于xs处的波场x分量和z分量震源函数；ρ表示地层密度，设置为常数；λ和μ表示拉梅系数，与纵波速度vP和横波速度vS之间存在如下关系：
Figure FDA0003021427830000033
5.根据权利要求4所述的基于角度分解与波场分离的弹性波全波形反演方法，其特征在于，步骤6中，对最低频率数据进行分解，得到频率角度域纵波波场uP(θ,x,ω)和横波波场uS(θ,x,ω)的具体公式如下：
Figure FDA0003021427830000034
Figure FDA0003021427830000035
式(5)中，uP(n,x,ω)和uS(n,x,ω)为分解得到的慢度域纵波波场和横波波场；
Figure FDA0003021427830000036
表示慢度向量，其中
Figure FDA0003021427830000037
为指向波场传播方向的单位向量，n为慢度向量n的绝对值；W(x′-x)表示以坐标点x为中心的采样窗；“×”和“·”分别表示叉乘运算和点乘运算；
式(6)中，nP和nS表示纵波慢度的绝对值和横波慢度的绝对值，即
Figure FDA0003021427830000041
其中
Figure FDA0003021427830000042
和
Figure FDA0003021427830000043
表示采样窗内的纵波平均速度和横波平均速度；结合式(5)和式(6)可以得到角度域纵波波场uP(θ,x,ω)和横波波场uS(θ,x,ω)；
步骤6具体实施方式如下：
根据式(5)和式(6)分别对震源侧波场
Figure FDA0003021427830000044
和接收点侧波场
Figure FDA0003021427830000045
进行分解，其中两个上标箭头为区分震源侧波场和接收点侧波场，得到震源侧角度域纵波波场
Figure FDA0003021427830000046
和横波波场
Figure FDA0003021427830000047
以及接收点侧角度域纵波波场
Figure FDA0003021427830000048
和横波波场
Figure FDA0003021427830000049
其中θi和θg分别表示震源侧波场传播方向和接收点侧波场传播方向，波场散射角通过如下公式计算得到：
Θ＝|θg-θi| (7)
式(7)中Θ表示波场散射角。
6.根据权利要求5所述的基于角度分解与波场分离的弹性波全波形反演方法，其特征在于，步骤7中，选取散射角范围依次为0°～30°，0°～60°和0°～180°。
7.根据权利要求5所述的基于角度分解与波场分离的弹性波全波形反演方法，其特征在于，步骤8中，分别利用纵波波场uP(θ,x,ω)和横波波场uS(θ,x,ω)通过迭代优化算法对目标函数J(vP(x),vS(x))进行优化的具体实施方式如下：
在弹性波传播过程中，纵波和横波经过散射会产生转换波，整个波场中存在多种模式波，即PP波，PS波，SP波，SS波；为了避免在反演过程中vP和vS的相互串扰问题，反演vP时只使用PP波，而反演vS时只使用SS、SP和PS波；为了避免局部极值问题，从小散射角到大散射角依次进行反演；因此在考虑角度分解与波场分离的情况下，vP和vS的梯度公式如下：
Figure FDA0003021427830000051
Figure FDA0003021427830000052
Figure FDA0003021427830000053
Figure FDA0003021427830000054
式(8)～(11)中，
Figure FDA0003021427830000055
和
Figure FDA0003021427830000056
分别对应于纵波速度vP(x)和横波速度vS(x)的角度滤波器，控制散射角在选定的范围之内；“Re”表示取实部；计算出梯度之后，利用如下准则对纵波速度vP(x)和横波速度vS(x)进行迭代更新：
Figure FDA0003021427830000061
式(12)中，vP(x)n,vP(x)n+1分别为第n步和第n+1步迭代的纵波速度值，vS(x)n,vS(x)n+1分别为第n步和第n+1步迭代的横波速度值；αn为迭代步长，为一个正常数；
通过上述方法计算速度梯度并进行迭代，在迭代收敛之后得到大尺度纵波速度vP1(x)和大尺度横波速度vS1(x)。
8.根据权利要求7所述的基于角度分解与波场分离的弹性波全波形反演方法，其特征在于，步骤9中，对分解后数据进行角度叠加，得到无角度信息的纵波波场和横波波场的公式如下：
Figure FDA0003021427830000062
Figure FDA0003021427830000063
式(13)和式(14)中，
Figure FDA0003021427830000064
和
Figure FDA0003021427830000065
表示震源侧无角度信息的纵波波场和横波波场；
Figure FDA0003021427830000066
和
Figure FDA0003021427830000067
表示接收点侧无角度信息的纵波波场和横波波场；
Figure FDA0003021427830000068
和
Figure FDA0003021427830000069
分别表示对震源侧波场传播角度求和以及对接收点侧波场传播角度求和。
9.根据权利要求8所述的基于角度分解与波场分离的弹性波全波形反演方法，其特征在于，步骤10中，分别利用纵波波场uP(x,ω)和横波波场uS(x,ω)通过迭代优化算法对目标函数J(vP(x),vS(x))进行优化的具体实施方式如下：在仅考虑波场分离的情况下，vP和vS的梯度公式如下：
Figure FDA00030214278300000610
Figure FDA0003021427830000071
Figure FDA0003021427830000072
Figure FDA0003021427830000073
计算出梯度之后，利用与式(12)相同的准则对纵波速度vP(x)和横波速度vS(x)进行迭代更新；通过上述方法计算速度梯度并进行迭代，在迭代收敛之后得到第二个频率分量下的纵波速度vP2(x)和大尺度横波速度vS2(x)。
CN201911122015.4A 2019-11-15 2019-11-15 基于角度分解与波场分离的弹性波全波形反演方法 Expired - Fee Related CN110888159B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911122015.4A CN110888159B (zh) 2019-11-15 2019-11-15 基于角度分解与波场分离的弹性波全波形反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911122015.4A CN110888159B (zh) 2019-11-15 2019-11-15 基于角度分解与波场分离的弹性波全波形反演方法

Publications (2)

Publication Number Publication Date
CN110888159A CN110888159A (zh) 2020-03-17
CN110888159B true CN110888159B (zh) 2021-08-06

Family

ID=69747684

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911122015.4A Expired - Fee Related CN110888159B (zh) 2019-11-15 2019-11-15 基于角度分解与波场分离的弹性波全波形反演方法

Country Status (1)

Country Link
CN (1) CN110888159B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113917522B (zh) * 2020-07-10 2024-03-19 中国石油化工股份有限公司 用于指导采集观测系统设计的地震正演方法
CN114721044B (zh) * 2022-04-21 2023-03-10 湖南工商大学 一种多频率接收函数和振幅比联合反演地壳结构的方法及系统
CN116540297B (zh) * 2023-05-06 2024-03-08 中国科学院地质与地球物理研究所 一种弹性波地震数据全波形反演方法、系统及设备

Family Cites Families (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20140372043A1 (en) * 2013-06-17 2014-12-18 Wenyi Hu Full Waveform Inversion Using Perfectly Reflectionless Subgridding
US11175421B2 (en) * 2014-01-10 2021-11-16 Cgg Services Sas Device and method for mitigating cycle-skipping in full waveform inversion
CN105891888B (zh) * 2016-03-28 2017-03-08 吉林大学 多域分频并行多尺度全波形反演方法
EP3685193B1 (en) * 2017-09-21 2022-10-26 Chevron U.S.A. Inc. System and method for improved full waveform inversion
US11353613B2 (en) * 2017-11-17 2022-06-07 Cgg Services Sas Seismic exploration using image-based reflection full waveform inversion to update low wavenumber velocity model
CN108181657B (zh) * 2017-12-28 2019-02-12 中国石油大学(北京) 全波形反演梯度计算中分离偏移和层析成像模式的方法
CN110441816B (zh) * 2019-09-20 2020-06-02 中国科学院测量与地球物理研究所 不依赖子波的无串扰多震源全波形反演方法及装置

Also Published As

Publication number Publication date
CN110888159A (zh) 2020-03-17

Similar Documents

Publication Publication Date Title
CN108873066B (zh) 弹性介质波动方程反射波旅行时反演方法
CN110888159B (zh) 基于角度分解与波场分离的弹性波全波形反演方法
CN110007340B (zh) 基于角度域直接包络反演的盐丘速度密度估计方法
CN108646293B (zh) 基于黏声拟微分方程的黏声起伏地表正演模拟系统及方法
CN109188519B (zh) 一种极坐标下的弹性波纵横波速度反演系统及方法
CN111239819B (zh) 一种基于地震道属性分析的带极性直接包络反演方法
CN110187382B (zh) 一种回折波和反射波波动方程旅行时反演方法
MX2011003850A (es) Estimado de señal de dominio de imagen a interferencia.
CN112327358A (zh) 一种粘滞性介质中声波地震数据正演模拟方法
CN109946742B (zh) 一种TTI介质中纯qP波地震数据模拟方法
Huang et al. A multi-block finite difference method for seismic wave equation in auxiliary coordinate system with irregular fluid–solid interface
CN110658558A (zh) 吸收衰减介质叠前深度逆时偏移成像方法及系统
CN111665556A (zh) 地层声波传播速度模型构建方法
Cao et al. 2.5-D modeling of seismic wave propagation: Boundary condition, stability criterion, and efficiency
CN111505714B (zh) 基于岩石物理约束的弹性波直接包络反演方法
CN111257930B (zh) 一种黏弹各向异性双相介质区域变网格求解算子
Raknes et al. Challenges and solutions for performing 3D time-domain elastic full-waveform inversion
CN114624766B (zh) 基于行波分离的弹性波最小二乘逆时偏移梯度求取方法
CN108680957B (zh) 基于加权的局部互相关时频域相位反演方法
CN113866823A (zh) 一种粘声各向异性介质中的正演成像方法
CN115373022A (zh) 一种基于振幅相位校正的弹性波场Helmholtz分解方法
CN116068621A (zh) 一种基于刚度矩阵分解的各向异性介质正演方法及系统
CN111665546B (zh) 用于可燃冰探测的声学参数获取方法
MX2011003852A (es) Atributos de procesamiento de imagen con inversion de tiempo.
Ouyang et al. Scattering pattern analysis and generalized Radon transform inversion for acoustic VTI media

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
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20210806