[go: up one dir, main page]
More Web Proxy on the site http://driver.im/

CN115480294A - 基于线性伯恩近似约束的非线性最小二乘逆时偏移方法 - Google Patents

基于线性伯恩近似约束的非线性最小二乘逆时偏移方法 Download PDF

Info

Publication number
CN115480294A
CN115480294A CN202110602166.0A CN202110602166A CN115480294A CN 115480294 A CN115480294 A CN 115480294A CN 202110602166 A CN202110602166 A CN 202110602166A CN 115480294 A CN115480294 A CN 115480294A
Authority
CN
China
Prior art keywords
order
wave field
coordinate position
wavefield
scattered
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.)
Granted
Application number
CN202110602166.0A
Other languages
English (en)
Other versions
CN115480294B (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 Petroleum and Chemical Corp
Geophysical Research Institute of Sinopec Shengli Oilfield Co
Original Assignee
China Petroleum and Chemical Corp
Geophysical Research Institute of Sinopec Shengli Oilfield Co
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 Petroleum and Chemical Corp, Geophysical Research Institute of Sinopec Shengli Oilfield Co filed Critical China Petroleum and Chemical Corp
Priority to CN202110602166.0A priority Critical patent/CN115480294B/zh
Publication of CN115480294A publication Critical patent/CN115480294A/zh
Application granted granted Critical
Publication of CN115480294B publication Critical patent/CN115480294B/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
    • G01V1/30—Analysis
    • G01V1/301—Analysis for determining seismic cross-sections or geostructures
    • 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/303—Analysis for determining velocity profiles or travel times
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/62—Physical property of subsurface

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)
  • Radar Systems Or Details Thereof (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明提供一种基于线性伯恩近似约束的非线性最小二乘逆时偏移方法，包括：步骤1，在给定速度扰动的情况下，通过观测数据和一阶散射波波动方程重构一阶散射波场；步骤2，推导基于该一阶散射波场计算速度扰动的更新梯度；步骤3，对基于线性伯恩近似约束的非线性LSRTM的反演问题凸性进行讨论；步骤4，根据一阶散射波场表达式计算的速度扰动更新梯度，得到一个高精度和高保幅的地震剖面。该基于线性伯恩近似约束的非线性最小二乘逆时偏移方法保证了一阶散射波场重构的准确性，从而使根据一阶散射波场表达式计算的速度扰动更新梯度也更加准确，使我们可以得到一个高精度和高保幅的地震剖面。

Description

基于线性伯恩近似约束的非线性最小二乘逆时偏移方法
技术领域
本发明涉及石油地球物理勘探技术领域，特别是涉及到一种基于线性伯恩近似约束的非线性最小二乘逆时偏移方法。
背景技术
最小二乘逆时偏移通过最小化模拟数据和观测数据之间的差可以给出高分辨率和高保幅的成像结果，而且相比于普通的逆时偏移方法不易产生低频噪音和假象。LSRTM理论是以Born近似下的一阶散射波假设为基础的，其对应的散射波表达式只能描述弱散射势和小散射体下的波场。但地震波实际传播过程是一个非线性过程，包含多种波现象，而且地下介质也十分复杂，弱散射势和小散射体条件一般难以满足。当利用全波场信息进行反演时，强散射势和大散射体产生的高阶散射波及直达波、回转波等多种波会与一阶散射波产生假象，影响成像质量。
在申请号：CN202010756625.6的中国专利申请中，涉及到一种优化的最小二乘逆时偏移成像方法，具体按照以下步骤实施：步骤1、采用波动方程对地震数据进行偏移与反偏移；步骤2、采用共轭梯度法构建最小二乘框架，并引入优化因子，对共轭梯度法的收敛速度进行提升；步骤3、采用GPU的CUDA计算架构进一步提升优化的LSRTM成像的计算效率。
在申请号：CN201811051319.1的中国专利申请中，涉及到一种基于L-BFGS算法的最小二乘偏移成像方法及装置，方法包括：根据目标勘探区的声波速度模型和声波速度模型的共炮点记录进行逆时偏移成像，得到初始成像结果；根据初始成像结果计算目标勘探区的合成地震记录和原始观测记录的记录残差；根据记录残差计算目标泛函关于初始成像结果的梯度；根据初始成像结果和目标泛函关于初始成像结果的梯度，通过L-BFGS算法得到第二成像结果；判断记录残差和L-BFGS算法的计算次数是否满足预设终止条件。
在申请号：CN201610855634.4的中国专利申请中，涉及到一种最小二乘偏移成像优化方法及系统。该方法包括：基于时窗内地震信号能量信息测度、时窗内地震信号对比度测度及时窗内结构相似度测度，获得地震数据结构相关系数；基于所述地震数据结构相关系数，建立由所述地震数据结构相关系数加权的最小二乘偏移成像模型；以及对由所述地震数据结构相关系数加权的最小二乘偏移成像模型求解，获得地震反射系数成像。
以上现有技术均与本发明有较大区别，未能解决我们想要解决的技术问题，为此我们发明了一种新的基于线性伯恩近似约束的非线性最小二乘逆时偏移方法。
发明内容
本发明的目的是提供一种基于线性伯恩近似约束的非线性最小二乘逆时偏移方法，用于得到一个高精度和高保幅的地震剖面的基于线性伯恩近似约束的非线性最小二乘逆时偏移方法。
本发明的目的可通过如下技术措施来实现：基于线性伯恩近似约束的非线性最小二乘逆时偏移方法，该基于线性伯恩近似约束的非线性最小二乘逆时偏移方法包括：
步骤1，在给定速度扰动的情况下，通过观测数据和一阶散射波波动方程重构一阶散射波场；
步骤2，推导基于该一阶散射波场计算速度扰动的更新梯度；
步骤3，对基于线性伯恩近似约束的非线性最小二乘逆时偏移LSRTM的反演问题凸性进行讨论；
步骤4，根据一阶散射波场表达式计算的速度扰动更新梯度，得到一个高精度和高保幅的地震剖面。
本发明的目的还可通过如下技术措施来实现：
在步骤1中，在地下介质的传播的标量场对应的齐次亥姆霍兹Helmholtz方程为：
Figure BDA0003092997860000031
其中，
Figure BDA0003092997860000032
为坐标位置，ω为角频率，Δ为哈密顿算子，u为波场值，Γ为正演算子，v为速度值。
引入背景速度场
Figure BDA0003092997860000033
得到：
Figure BDA0003092997860000034
其中，
Figure BDA0003092997860000035
为坐标位置，ω为角频率，Δ为哈密顿算子，u为波场值，v为速度值，v0为背景速度。
定义散射势为：
Figure BDA0003092997860000036
其中，
Figure BDA0003092997860000037
为坐标位置，ω为角频率，v为速度值，V为散射势，v0为背景速度，
Figure BDA0003092997860000038
即为LSRTM感兴趣的速度扰动。
根据散射势，上式化简为：
Figure BDA0003092997860000039
其中，
Figure BDA00030929978600000310
为坐标位置，ω为角频率，Δ为哈密顿算子，u为波场值，v0为背景速度场，V为散射势。
其特征在于，在步骤1中，假设背景介质中传播的入射波场满足齐次Helmholtz方程：
Figure BDA00030929978600000311
其中，
Figure BDA00030929978600000312
为坐标位置，ω为角频率，Δ为哈密顿算子，u为波场值,v0为背景速度场。
则式(5)化简为一个非齐次的Helmholtz方程，描述了在背景介质中传播的入射波场遇到散射体之后产生的散射波场：
Figure BDA0003092997860000041
其中，
Figure BDA0003092997860000042
为坐标位置，ω为角频率，Δ为哈密顿算子，u为波场值,v0为背景速度场，us为散射波场，V为散射势。
定义格林Green函数：
Figure BDA0003092997860000043
其中，
Figure BDA0003092997860000044
为坐标位置，ω为角频率，Δ为哈密顿算子，G为格林函数，δ为脉冲函数，
Figure BDA0003092997860000045
为移动坐标位置。
利用Green函数得到散射波场的解析解：
Figure BDA0003092997860000046
其中,us为散射波场，V为散射势，u为波场值，G为格林函数，
Figure BDA0003092997860000047
为坐标位置，
Figure BDA0003092997860000048
为移动坐标位置，ω为角频。
由于上式是一个体积分的表达式，所以散射波场由散射势及散射体大小控制。
在步骤1中，在一阶伯恩Born近似的假设下，即散射势较弱，散射体较小，使得产生的散射波场远小于背景波场：
Figure BDA0003092997860000049
得到一阶Born近似下的散射波场：
Figure BDA00030929978600000410
其中，
Figure BDA00030929978600000411
为一阶散射波场，V为散射势，u为波场值，G为格林函数，
Figure BDA0003092997860000051
为移动坐标位置，ω为角频率，
Figure BDA0003092997860000052
为检波点坐标位置，Ω为积分体。
对应的全波场为：
Figure BDA0003092997860000053
其中，u为波场值，u0为背景波场值，V为散射势，G为格林函数，ω为角频率，
Figure BDA0003092997860000054
为检波点坐标位置，Ω为积分体,
Figure BDA0003092997860000055
为移动坐标位置。
将公式(10)带入到公式(8)中，得到二阶散射波表达式；重复上述过程即可得到多阶散射波场表达式，用求和形式简单表达N阶散射波场：
Figure BDA0003092997860000056
其中，us为散射波场，
Figure BDA0003092997860000057
为第i阶散射波场，N为总阶数，
Figure BDA0003092997860000058
为坐标位置，ω为角频率。
在步骤2中，当采用一阶Born近似下的散射波方程作为正演算子时，L-2范数下的LSRTM目标泛函表示为：
Figure BDA0003092997860000059
其中，J为目标泛函，
Figure BDA00030929978600000510
为一阶模拟散射波场，
Figure BDA00030929978600000511
为一阶散射波场，dn为第n阶背景波场。
其对应的速度扰动更新的梯度表示为：
Figure BDA00030929978600000512
其中，
Figure BDA00030929978600000513
为速度扰动，J为目标泛函，
Figure BDA00030929978600000514
为坐标位置，ω为角频率，
Figure BDA00030929978600000515
为一阶模拟散射波场，
Figure BDA00030929978600000516
为一阶模拟波场，dn为第n阶背景波场，u0为背景波场值，V为散射势，G为格林函数，
Figure BDA0003092997860000061
为检波点坐标位置。
分析可知，上述梯度包含两部分：模拟背景波场与模拟的一阶散射波场与观测数据中的一阶散射波场之间的残差所构成的梯度：
Figure BDA0003092997860000062
其中，
Figure BDA0003092997860000063
为一阶梯度，
Figure BDA0003092997860000064
为坐标位置，ω为角频率，u0为背景波场值，G为格林函数，
Figure BDA0003092997860000065
为检波点坐标位置，
Figure BDA0003092997860000066
为移动坐标位置，
Figure BDA0003092997860000067
为一阶模拟散射波场，
Figure BDA00030929978600000612
为一阶模拟波场。
模拟背景波场与观测数据中的背景波场和高阶Born近似下的观测数据所构成的梯度：
Figure BDA0003092997860000068
其中，gn为第n阶梯度，ω为角频率，u0为背景波场值，
Figure BDA0003092997860000069
为坐标位置，
Figure BDA00030929978600000610
为检波点坐标位置，dn为第n阶背景波场，G为格林函数，
Figure BDA00030929978600000611
为移动坐标位置。
梯度第一项描述了正演算子所描述的一阶散射波对应的速度扰动，而第二项在一阶散射波假设下表现为成像噪音，即在一阶Born近似下，观测数据的其他分量由于无法被准确描述，而与背景波场产生了对速度扰动梯度的干扰。
在步骤2中，为了压制观测数据中背景波场与高阶散射波场产生的噪音，将线性伯恩近似方程作为约束项加入到非线性LSRTM的目标函数中，加强非线性LSRTM中一阶散射波的约束，以获得高精度，高分辨率的剖面。
在步骤2中，改进后的LSRTM目标函数表示为：
Figure BDA0003092997860000071
其中，λ为约束因子，根据观测数据中的背景波场和高阶散射波场的强度进行调整，J为目标泛函，
Figure BDA0003092997860000072
为坐标位置，ω为角频率，
Figure BDA0003092997860000073
为一阶模拟散射波场，
Figure BDA0003092997860000074
为一阶模拟波场，dn为第n阶背景波场，u0为背景波场值，V为散射势，G为格林函数,
Figure BDA0003092997860000075
为检波点坐标位置,
Figure BDA0003092997860000076
为移动坐标位置。
添加约束后的LSRTM分解为两个反演问题：首先针对一阶散射波场进行反演，得到相对准确的一阶散射波场；然后在此基础上计算速度扰动的更新梯度；令公式(16)对于一阶散射波场的导数为0，可以得到一阶散射波场的重构表达式：
Figure BDA0003092997860000077
其中，λ为约束因子，根据观测数据中的背景波场和高阶散射波场的强度进行调整，J为目标泛函，
Figure BDA0003092997860000078
为坐标位置，ω为角频率，
Figure BDA0003092997860000079
为一阶模拟散射波场，
Figure BDA00030929978600000710
为一阶模拟波场，dn为第n阶背景波场，u0为背景波场值，V为散射势，G为格林函数,
Figure BDA00030929978600000711
为检波点坐标位置,
Figure BDA00030929978600000712
为移动坐标位置。
Figure BDA00030929978600000713
其中，λ为约束因子，
Figure BDA00030929978600000714
为坐标位置，ω为角频率，us 1为一阶模拟散射波场，
Figure BDA0003092997860000081
为一阶模拟波场，dn为第n阶背景波场，u0为背景波场值，V为散射势，G为格林函数,
Figure BDA0003092997860000089
为移动坐标位置。
分析式(18)可看出，重构的一阶散射波场主要由两部分组成：观测数据构成的波场和模拟的一阶散射波场，两者的比重由约束因子控制：当约束因子较大时，一阶散射波对重构散射波场的影响变大，从而减弱了背景波场和高级散射波对一阶重构散射波的影响。
在步骤2中，在准确重构散射波场之后，得到目标函数对速度扰动的梯度：
Figure BDA0003092997860000082
其中，λ为约束因子，J为目标泛函，
Figure BDA0003092997860000083
为坐标位置，ω为角频率，u0为背景波场值，V为散射势，
Figure BDA0003092997860000084
为移动坐标位置，G为格林函数,
Figure BDA0003092997860000085
为速度扰动，
Figure BDA0003092997860000086
为一阶模拟散射波场，u为波场值，u0为背景波场值。
进而对速度扰动进行更新：
Figure BDA0003092997860000087
其中，α为更新步长,J为目标泛函,
Figure BDA0003092997860000088
为速度扰动。。
在步骤3中，针对一般的反演问题，其目标函数对与反演参数的二阶导数即Hessian矩阵的性质决定了反演问题的凸性：当Hessian矩阵时，Hessian矩阵为半正定的，即目标函数是全凸的，通过局部寻优法寻找到目标函数的全局极值；反之，即目标函数是局部凸的，反演问题受局部极值的影响，难以收敛到全局极值。
在步骤3中，基于线性伯恩近似约束的非线性LSRTM对一阶散射波场的Hessian矩阵为：
Figure BDA0003092997860000091
其中，λ为约束因子，J为目标泛函，
Figure BDA0003092997860000092
为坐标位置，
Figure BDA0003092997860000093
为坐标位置，ω为角频率，u0为背景波场值，V为散射势，G为格林函数,
Figure BDA0003092997860000094
为移动坐标位置,
Figure BDA0003092997860000095
为一阶模拟散射波场，
Figure BDA0003092997860000096
为一阶模拟波场，dn为第n阶背景波场，
Figure BDA0003092997860000097
为检波点坐标位置。
由一阶散射波场的表达式可看出，波场受模型参数、散射势及观测系统影响；不同观测系统下的波场是不相关的，则：
Figure BDA0003092997860000098
其中，H为二阶梯度，r为检波点位置，λ为约束因子，J为目标泛函，
Figure BDA0003092997860000099
为坐标位置，
Figure BDA00030929978600000910
为坐标位置，ω为角频率，
Figure BDA00030929978600000911
为一阶模拟散射波场。
即基于线性伯恩近似约束的非线性LSRTM针对一阶散射波场表现出强凸性，用局部寻优的方法得到较为准确的一阶散射波场；此外，Born近似下的一阶散射波场与速度扰动之间是线性的，因此基于重构的一阶散射波场得到准确的速度扰动更新梯度。
本发明中的基于线性伯恩近似约束的非线性最小二乘逆时偏移方法，提出将线性伯恩近似作为约束项加入到非线性LSRTM的目标函数中，在反演过程中通过增强线性伯恩近似的约束来压制由弱散射势和小散射体假设引起的假象。基于线性伯恩近似约束的反演方法将非线性LSRTM问题分解为两个反演问题：首先在一阶散射波场空间寻找最优解，即在给定速度扰动的情况下，通过观测数据和一阶散射波波动方程尽可能重构一个准确的一阶散射波场。然后基于该一阶散射波场计算速度扰动的更新梯度。新的目标函数对一阶散射波场的Hessian矩阵表现出强凸性，保证了一阶散射波场重构的准确性，从而使根据一阶散射波场表达式计算的速度扰动更新梯度也更加准确，使我们可以得到一个高精度和高保幅的地震剖面。
本发明的有益效果是：本发明提供了基于线性伯恩近似约束的非线性最小二乘逆时偏移方法，通过观测数据和一阶散射波波动方程尽可能重构一个准确的一阶散射波场。然后基于该一阶散射波场计算速度扰动的更新梯度。新的目标函数对一阶散射波场的Hessian矩阵表现出强凸性，保证了一阶散射波场重构的准确性，从而使根据一阶散射波场表达式计算的速度扰动更新梯度也更加准确，使我们可以得到一个高精度和高保幅的地震剖面。
附图说明
图1为本发明的具体实施例中Sigsbee2A模型的示意图；
图2为本发明的具体实施例1中反演用炮记录的示意图；
图3为本发明的具体实施例1中Sigsbee2A偏移速度模型的示意图；
图4为本发明的具体实施例1中迭代10次的LSRTM结果的示意图；
图5为本发明的具体实施例1中针对图4盐上部分的放大图；
图6为本发明的具体实施例1中针对图4盐下部分的放大图；
图7为本发明的具体实施例2中反演用炮记录的示意图；
图8为本发明的具体实施例2中Sigsbee2A偏移速度模型的示意图；
图9为本发明的具体实施例2中迭代10次的LSRTM结果的示意图；
图10为本发明的具体实施例2中针对图4盐上部分的放大图；
图11为本发明的具体实施例2中针对图4盐下部分的放大图；
图12为本发明的具体实施例3中反演用炮记录的示意图；
图13为本发明的具体实施例3中Sigsbee2A偏移速度模型的示意图；
图14为本发明的具体实施例3中迭代10次的LSRTM结果的示意图；
图15为本发明的具体实施例3中针对图4盐上部分的放大图；
图16为本发明的具体实施例3中针对图4盐下部分的放大图；
图17为本发明的基于线性伯恩近似约束的非线性最小二乘逆时偏移方法的一具体实施例的流程图。
具体实施方式
应该指出，以下详细说明都是示例性的，旨在对本发明提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本发明所属技术领域的普通技术人员通常理解的相同含义。
需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本发明的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作和/或它们的组合。
如图17所示，图17为本发明的基于线性伯恩近似约束的非线性最小二乘逆时偏移方法的流程图。该基于线性伯恩近似约束的非线性最小二乘逆时偏移方法包括了以下步骤：
步骤101，在给定速度扰动的情况下，通过观测数据和一阶散射波波动方程尽可能重构一个准确的一阶散射波场；
在地下介质的传播的标量场对应的齐次Helmholtz方程为：
Figure BDA0003092997860000111
其中，
Figure BDA0003092997860000112
为坐标位置，ω为角频率，Δ为哈密顿算子，u为波场值，Γ为正演算子，v为速度值。
引入背景速度场
Figure BDA0003092997860000113
我们可以得到：
Figure BDA0003092997860000114
其中，
Figure BDA0003092997860000115
为坐标位置，ω为角频率，Δ为哈密顿算子，u为波场值，v为速度值，v0为背景速度。
定义散射势为：
Figure BDA0003092997860000121
其中，
Figure BDA0003092997860000122
为坐标位置，ω为角频率，v为速度值，V为散射势，v0为背景速度，
Figure BDA0003092997860000123
即为LSRTM感兴趣的速度扰动。根据散射势，上式可以化简为：
Figure BDA0003092997860000124
其中，
Figure BDA0003092997860000125
为坐标位置，ω为角频率，Δ为哈密顿算子，u为波场值，v0为背景速度场，V为散射势。
假设背景介质中传播的入射波场满足齐次Helmholtz方程：
Figure BDA0003092997860000126
其中，
Figure BDA0003092997860000127
为坐标位置，ω为角频率，Δ为哈密顿算子，u为波场值,v0为背景速度场。
则式(5)化简为一个非齐次的Helmholtz方程，描述了在背景介质中传播的入射波场遇到散射体之后产生的散射波场：
Figure BDA0003092997860000128
其中，
Figure BDA0003092997860000129
为坐标位置，ω为角频率，Δ为哈密顿算子，u为波场值,v0为背景速度场，us为散射波场，V为散射势。
定义Green函数：
Figure BDA0003092997860000131
其中，
Figure BDA0003092997860000132
为坐标位置，ω为角频率，Δ为哈密顿算子，G为格林函数，δ为脉冲函数，
Figure BDA0003092997860000133
为移动坐标位置。
利用Green函数可以得到散射波场的解析解：
Figure BDA0003092997860000134
其中,us为散射波场，V为散射势，u为波场值，G为格林函数，
Figure BDA0003092997860000135
为坐标位置，
Figure BDA0003092997860000136
为移动坐标位置，ω为角频率。
由于上式是一个体积分的表达式，所以散射波场由散射势及散射体大小控制。在一阶Born近似的假设下，即散射势较弱，散射体较小，使得产生的散射波场远小于背景波场：
Figure BDA0003092997860000137
我们可以得到一阶Born近似下的散射波场：
Figure BDA0003092997860000138
其中，
Figure BDA0003092997860000139
为一阶散射波场，V为散射势，u为波场值，G为格林函数，
Figure BDA00030929978600001310
为移动坐标位置，ω为角频率，
Figure BDA00030929978600001311
为检波点坐标位置，Ω为积分体。
对应的全波场为：
Figure BDA00030929978600001312
其中，u为波场值，u0为背景波场值，V为散射势，G为格林函数，ω为角频率，
Figure BDA00030929978600001313
为检波点坐标位置，Ω为积分体,
Figure BDA00030929978600001314
为移动坐标位置。
将公式(10)带入到公式(8)中，可以得到二阶散射波表达式。重复上述过程即可得到多阶散射波场表达式，用求和形式简单表达N阶散射波场：
Figure BDA0003092997860000141
其中，us为散射波场，
Figure BDA0003092997860000142
为第i阶散射波场，N为总阶数，
Figure BDA0003092997860000143
为坐标位置，ω为角频率。
步骤102，推导基于该一阶散射波场计算速度扰动的更新梯度；
推导基于该一阶散射波场计算速度扰动的更新梯度，包括：
当采用一阶Born近似下的散射波方程作为正演算子时，L-2范数下的LSRTM目标泛函可以表示为：
Figure BDA0003092997860000144
其中，J为目标泛函，
Figure BDA0003092997860000145
为一阶模拟散射波场，
Figure BDA0003092997860000146
为一阶散射波场，dn为第n阶背景波场。
其对应的速度扰动更新的梯度可以表示为：
Figure BDA0003092997860000147
其中，
Figure BDA0003092997860000148
为速度扰动，J为目标泛函，
Figure BDA0003092997860000149
为坐标位置，ω为角频率，
Figure BDA00030929978600001410
为一阶模拟散射波场，
Figure BDA00030929978600001411
为一阶模拟波场，dn为第n阶背景波场，u0为背景波场值，V为散射势，G为格林函数，
Figure BDA00030929978600001412
为检波点坐标位置。
分析可知，上述梯度包含两部分：模拟背景波场与模拟的一阶散射波场与观测数据中的一阶散射波场之间的残差所构成的梯度：
Figure BDA00030929978600001413
其中，
Figure BDA00030929978600001414
为一阶梯度，
Figure BDA00030929978600001415
为坐标位置，ω为角频率，u0为背景波场值，G为格林函数，
Figure BDA00030929978600001416
为检波点坐标位置，
Figure BDA00030929978600001417
为移动坐标位置，
Figure BDA00030929978600001418
为一阶模拟散射波场，
Figure BDA0003092997860000151
为一阶模拟波场。
模拟背景波场与观测数据中的背景波场和高阶Born近似下的观测数据所构成的梯度：
Figure BDA0003092997860000152
其中，gn为第n阶梯度，ω为角频率，u0为背景波场值，
Figure BDA0003092997860000153
为坐标位置，
Figure BDA0003092997860000154
为检波点坐标位置，dn为第n阶背景波场，G为格林函数，
Figure BDA0003092997860000155
为移动坐标位置。
显然梯度第一项描述了正演算子所描述的一阶散射波对应的速度扰动，而第二项在一阶散射波假设下表现为成像噪音，即在一阶Born近似下，观测数据的其他分量由于无法被准确描述，而与背景波场产生了对速度扰动梯度的干扰。为了压制观测数据中背景波场与高阶散射波场产生的噪音，将线性伯恩近似方程作为约束项加入到非线性LSRTM的目标函数中，加强非线性LSRTM中一阶散射波的约束，以获得高精度，高分辨率的剖面。改进后的LSRTM目标函数可以表示为：
Figure BDA0003092997860000156
其中，λ为约束因子，根据观测数据中的背景波场和高阶散射波场的强度进行调整，J为目标泛函，
Figure BDA0003092997860000157
为坐标位置，ω为角频率，
Figure BDA0003092997860000158
为一阶模拟散射波场，
Figure BDA0003092997860000159
为一阶模拟波场，dn为第n阶背景波场，u0为背景波场值，V为散射势，G为格林函数,
Figure BDA00030929978600001510
为检波点坐标位置,
Figure BDA00030929978600001511
为移动坐标位置。
添加约束后的LSRTM可以分解为两个反演问题：首先针对一阶散射波场进行反演，得到相对准确的一阶散射波场；然后在此基础上计算速度扰动的更新梯度。令公式(16)对于一阶散射波场的导数为0，可以得到一阶散射波场的重构表达式：
Figure BDA0003092997860000161
其中，λ为约束因子，根据观测数据中的背景波场和高阶散射波场的强度进行调整，J为目标泛函，
Figure BDA0003092997860000162
为坐标位置，ω为角频率，
Figure BDA0003092997860000163
为一阶模拟散射波场，
Figure BDA0003092997860000164
为一阶模拟波场，dn为第n阶背景波场，u0为背景波场值，V为散射势，G为格林函数,
Figure BDA0003092997860000165
为检波点坐标位置,
Figure BDA0003092997860000166
为移动坐标位置。
Figure BDA0003092997860000167
其中，λ为约束因子，
Figure BDA0003092997860000168
为坐标位置，ω为角频率，
Figure BDA0003092997860000169
为一阶模拟散射波场，
Figure BDA00030929978600001610
为一阶模拟波场，dn为第n阶背景波场，u0为背景波场值，V为散射势，G为格林函数,
Figure BDA00030929978600001611
为移动坐标位置。
分析式(18)可以看出，重构的一阶散射波场主要由两部分组成：观测数据构成的波场和模拟的一阶散射波场，两者的比重由约束因子控制：当约束因子较大时，一阶散射波对重构散射波场的影响变大，从而减弱了背景波场和高级散射波对一阶重构散射波的影响。
在准确重构散射波场之后，可以得到目标函数对速度扰动的梯度：
Figure BDA00030929978600001612
其中，λ为约束因子，J为目标泛函，
Figure BDA00030929978600001613
为坐标位置，ω为角频率，u0为背景波场值，V为散射势，
Figure BDA0003092997860000171
为移动坐标位置，G为格林函数,
Figure BDA0003092997860000172
为速度扰动，
Figure BDA0003092997860000173
为一阶模拟散射波场，u为波场值，u0为背景波场值。
进而可以对速度扰动进行更新：
Figure BDA0003092997860000174
其中，α为更新步长,J为目标泛函,
Figure BDA0003092997860000175
为速度扰动。。
步骤103，对基于线性伯恩近似约束的非线性LSRTM的反演问题凸性进行讨论；针对一般的反演问题，其目标函数对与反演参数的二阶导数(Hessian矩阵)的性质一般决定了反演问题的凸性：当Hessian矩阵时，Hessian矩阵为半正定的，即目标函数是全凸的，可以通过局部寻优法寻找到目标函数的全局极值；反之，即目标函数是局部凸的，反演问题受局部极值的影响，难以收敛到全局极值。基于线性伯恩近似约束的非线性LSRTM对一阶散射波场的Hessian矩阵为：
Figure BDA0003092997860000176
其中，λ为约束因子，J为目标泛函，
Figure BDA0003092997860000177
为坐标位置，
Figure BDA0003092997860000178
为坐标位置，ω为角频率，u0为背景波场值，V为散射势，G为格林函数,
Figure BDA0003092997860000179
为移动坐标位置,
Figure BDA00030929978600001710
为一阶模拟散射波场，
Figure BDA00030929978600001711
为一阶模拟波场，dn为第n阶背景波场，
Figure BDA00030929978600001712
为检波点坐标位置。
由一阶散射波场的表达式可以看出，波场受模型参数、散射势及观测系统影响。不同观测系统下的波场是不相关的，则：
Figure BDA0003092997860000181
其中，λ为约束因子，J为目标泛函，
Figure BDA0003092997860000182
为坐标位置，
Figure BDA0003092997860000183
为坐标位置，ω为角频率，
Figure BDA0003092997860000184
为一阶模拟散射波场。
即基于线性伯恩近似约束的非线性LSRTM针对一阶散射波场表现出强凸性，可以用局部寻优的方法得到较为准确的一阶散射波场。此外，Born近似下的一阶散射波场与速度扰动之间是线性的，因此我们可以基于重构的一阶散射波场得到准确的速度扰动更新梯度。
步骤104，根据一阶散射波场表达式计算的速度扰动更新梯度，得到一个高精度和高保幅的地震剖面。
实施例1：
在应用本发明的具体实施例1中，本发明基于一次散射波约束的最小二乘逆时偏移，首先用于Sigsbee2A模型(图1)测试其适应性。Sigsbee2A模型包含一个高速的盐丘，因此无法满足弱散射势和小散射体的假设，可以很好的说明本发明方法的有效性。模型大小为401*1601，网格间距纵向为21m，横向22m。同样采用时间域声波波动方程合成地震记录，采样时间7s，采样间隔1ms。表面放置了80炮，炮点间隔440m，表面全孔径接收。图2展示了用于反演的一个炮记录，该合成记录包含直达波，回转波等多种波现象，由于盐丘的存在以及深层模型的复杂性，其炮记录还包括层间多次波，高阶散射波等现象。用于反演的背景速度如图3所示，模型除大型盐丘没有明显的散射体。在进行反演测试时，为了体现约束因子对反演的重要性，分别用普通的LSRTM，约束因子为1000的LSRTM和约束因子为100000的LSRTM进行反演，其最终反演结果如图4所示，图4(a)为常规LSRTM结果；(b)为基于一阶散射波约束的LSRTM结果(约束因子为1000)；(c)为基于一阶散射波约束的LSRTM结果(约束因子为100000)。同样由于直达波及回转波能量过强，深层的反演结果并不明显。
为了更加明显的对比两种方法的成像结果。首先截取盐上一部分进行比较：由于直达波和回转波的影响，普通LSRTM给出的反演结果中有大量的干扰噪音，浅层画弧现象严重。引入一次散射波的约束之后，干扰噪音明显减弱，反演结果的分辨率也进一步提高。随着进一步加强一阶散射波的约束，由直达波和回转波等背景波场产生的噪音已基本消除，反演结果浅层的层位更加明确。图5为针对图4盐上部分的放大图。
然后选取盐下一部分方法进行精确比较：由于观测数据的复杂性，普通的LSRTM结果中包含大量的噪音，而且并没有反演出散射点，而增加了约束的LSRTM的反演结果保真度更好，分辨率也有所提高，但同样没有反演出散射点；通过进一步加强一阶散射波的约束，LSRTM最终给出了较好的反演结果：盐下的层位和散射点位置准确，分辨率较高。图6为针对图4盐下部分的放大图。
最小二乘逆时偏移通过最小化模拟数据和观测数据之间的差可以给出高分辨率和高保幅的成像结果，而且相比于普通的逆时偏移方法不易产生低频噪音和假象。LSRTM理论是以Born近似下的一阶散射波假设为基础的，其对应的散射波表达式只能描述弱散射势和小散射体下的波场。但地震波实际传播过程是一个非线性过程，包含多种波现象，而且地下介质也十分复杂，弱散射势和小散射体条件一般难以满足。当利用全波场信息进行反演时，强散射势和大散射体产生的高阶散射波及直达波、回转波等多种波会与一阶散射波产生假象，影响成像质量。
因此，考虑到上述技术存在的缺陷，本发明提出了将线性伯恩近似的表达式作为约束项加入到非线性LSRTM的目标函数中，在反演过程中通过增强一阶散射波的约束来压制由弱散射势和小散射体假设引起的假象。
实施例2：
在应用本发明的具体实施例2中，本发明基于一次散射波约束的最小二乘逆时偏移，首先用于Sigsbee2A模型(图1)测试其适应性。Sigsbee2A模型包含一个高速的盐丘，因此无法满足弱散射势和小散射体的假设，可以很好的说明本发明方法的有效性。模型大小为201*801，网格间距纵向为42m，横向44m。同样采用时间域声波波动方程合成地震记录，采样时间7s，采样间隔1ms。表面放置了80炮，炮点间隔440m，表面全孔径接收。图7展示了用于反演的一个炮记录，该合成记录包含直达波，回转波等多种波现象，由于盐丘的存在以及深层模型的复杂性，其炮记录还包括层间多次波，高阶散射波等现象。用于反演的背景速度如图8所示，模型除大型盐丘没有明显的散射体。在进行反演测试时，为了体现约束因子对反演的重要性，分别用普通的LSRTM，约束因子为1000的LSRTM和约束因子为100000的LSRTM进行反演，其最终反演结果如图9所示，图9(a)为常规LSRTM结果；(b)为基于一阶散射波约束的LSRTM结果(约束因子为1000)；(c)为基于一阶散射波约束的LSRTM结果(约束因子为100000)。同样由于直达波及回转波能量过强，深层的反演结果并不明显。
为了更加明显的对比两种方法的成像结果。首先截取盐上一部分进行比较：由于直达波和回转波的影响，普通LSRTM给出的反演结果中有大量的干扰噪音，浅层画弧现象严重。引入一次散射波的约束之后，干扰噪音明显减弱，反演结果的分辨率也进一步提高。随着进一步加强一阶散射波的约束，由直达波和回转波等背景波场产生的噪音已基本消除，反演结果浅层的层位更加明确。图12为针对图9盐上部分的放大图。
然后选取盐下一部分方法进行精确比较：由于观测数据的复杂性，普通的LSRTM结果中包含大量的噪音，而且并没有反演出散射点，而增加了约束的LSRTM的反演结果保真度更好，分辨率也有所提高，但同样没有反演出散射点；通过进一步加强一阶散射波的约束，LSRTM最终给出了较好的反演结果：盐下的层位和散射点位置准确，分辨率较高。图11为针对图9盐下部分的放大图。
实施例3：
在应用本发明的具体实施例3中，本发明基于一次散射波约束的最小二乘逆时偏移，首先用于Sigsbee2A模型(图1)测试其适应性。Sigsbee2A模型包含一个高速的盐丘，因此无法满足弱散射势和小散射体的假设，可以很好的说明本发明方法的有效性。模型大小为101*401，网格间距纵向为84m，横向88m。同样采用时间域声波波动方程合成地震记录，采样时间7s，采样间隔1ms。表面放置了80炮，炮点间隔440m，表面全孔径接收。图8展示了用于反演的一个炮记录，该合成记录包含直达波，回转波等多种波现象，由于盐丘的存在以及深层模型的复杂性，其炮记录还包括层间多次波，高阶散射波等现象。用于反演的背景速度如图14所示，模型除大型盐丘没有明显的散射体。在进行反演测试时，为了体现约束因子对反演的重要性，分别用普通的LSRTM，约束因子为1000的LSRTM和约束因子为100000的LSRTM进行反演，其最终反演结果如图14所示，图14(a)为常规LSRTM结果；(b)为基于一阶散射波约束的LSRTM结果(约束因子为1000)；(c)为基于一阶散射波约束的LSRTM结果(约束因子为100000)。同样由于直达波及回转波能量过强，深层的反演结果并不明显。
为了更加明显的对比两种方法的成像结果。首先截取盐上一部分进行比较：由于直达波和回转波的影响，普通LSRTM给出的反演结果中有大量的干扰噪音，浅层画弧现象严重。引入一次散射波的约束之后，干扰噪音明显减弱，反演结果的分辨率也进一步提高。随着进一步加强一阶散射波的约束，由直达波和回转波等背景波场产生的噪音已基本消除，反演结果浅层的层位更加明确。图15为针对图14盐上部分的放大图。
然后选取盐下一部分方法进行精确比较：由于观测数据的复杂性，普通的LSRTM结果中包含大量的噪音，而且并没有反演出散射点，而增加了约束的LSRTM的反演结果保真度更好，分辨率也有所提高，但同样没有反演出散射点；通过进一步加强一阶散射波的约束，LSRTM最终给出了较好的反演结果：盐下的层位和散射点位置准确，分辨率较高。图16为针对图14盐下部分的放大图。
最后应说明的是：以上所述仅为本发明的优选实施例而已，并不用于限制本发明，尽管参照前述实施例对本发明进行了详细的说明，对于本领域技术人员来说，其依然可以对前述实施例记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
除说明书所述的技术特征外，均为本专业技术人员的已知技术。

Claims (10)

1.基于线性伯恩近似约束的非线性最小二乘逆时偏移方法，其特征在于，该基于线性伯恩近似约束的非线性最小二乘逆时偏移方法包括：
步骤1，在给定速度扰动的情况下，通过观测数据和一阶散射波波动方程重构一阶散射波场；
步骤2，推导基于该一阶散射波场计算速度扰动的更新梯度；
步骤3，对基于线性伯恩近似约束的非线性LSRTM的反演问题凸性进行讨论；
步骤4，根据一阶散射波场表达式计算的速度扰动更新梯度，得到一个高精度和高保幅的地震剖面。
2.根据权利要求1所述的基于线性伯恩近似约束的非线性最小二乘逆时偏移方法，其特征在于，在步骤1中，在地下介质的传播的标量场对应的齐次Helmholtz方程为：
Figure FDA0003092997850000011
其中，
Figure FDA0003092997850000012
为坐标位置，ω为角频率，Δ为哈密顿算子，u为波场值，Γ为正演算子，v为速度值；
引入背景速度场
Figure FDA0003092997850000013
得到：
Figure FDA0003092997850000014
其中，
Figure FDA0003092997850000015
为坐标位置，ω为角频率，Δ为哈密顿算子，u为波场值，v为速度值，v0为背景速度；
定义散射势为：
Figure FDA0003092997850000016
其中，
Figure FDA0003092997850000017
为坐标位置，ω为角频率，v为速度值，V为散射势，v0为背景速度，
Figure FDA0003092997850000021
即为LSRTM感兴趣的速度扰动；
根据散射势，上式化简为：
Figure FDA0003092997850000022
其中，
Figure FDA0003092997850000023
为坐标位置，ω为角频率，Δ为哈密顿算子，u为波场值，v0为背景速度场，V为散射势。
3.根据权利要求2所述的基于线性伯恩近似约束的非线性最小二乘逆时偏移方法，其特征在于，在步骤1中，假设背景介质中传播的入射波场满足齐次Helmholtz方程：
Figure FDA0003092997850000024
其中，
Figure FDA0003092997850000025
为坐标位置，ω为角频率，Δ为哈密顿算子，u为波场值,v0 为背景速度场；
则式(5)化简为一个非齐次的Helmholtz方程，描述了在背景介质中传播的入射波场遇到散射体之后产生的散射波场：
Figure FDA0003092997850000026
其中，
Figure FDA0003092997850000027
为坐标位置，ω为角频率，Δ为哈密顿算子，u为波场值,v0 为背景速度场，us为散射波场，V为散射势；
定义Green函数：
Figure FDA0003092997850000028
其中，
Figure FDA0003092997850000029
为坐标位置，ω为角频率，Δ为哈密顿算子，G为格林函数，δ为脉冲函数，
Figure FDA0003092997850000031
为移动坐标位置；
利用Green函数得到散射波场的解析解：
Figure FDA0003092997850000032
其中,us为散射波场，V为散射势，u为波场值，G为格林函数，
Figure FDA0003092997850000033
为坐标位置，
Figure FDA0003092997850000034
为移动坐标位置，ω为角频；
由于上式是一个体积分的表达式，所以散射波场由散射势及散射体大小控制。
4.根据权利要求3所述的基于线性伯恩近似约束的非线性最小二乘逆时偏移方法，其特征在于，在步骤1中，在一阶Born近似的假设下，即散射势较弱，散射体较小，使得产生的散射波场远小于背景波场：
Figure FDA0003092997850000035
得到一阶Born近似下的散射波场：
Figure FDA0003092997850000036
其中，
Figure FDA0003092997850000037
为一阶散射波场，V为散射势，u为波场值，G为格林函数，
Figure FDA0003092997850000038
为移动坐标位置，ω为角频率，
Figure FDA0003092997850000039
为检波点坐标位置，Ω为积分体；
对应的全波场为：
Figure FDA00030929978500000310
其中，u为波场值，u0为背景波场值，V为散射势，G为格林函数，ω为角频率，
Figure FDA00030929978500000311
为检波点坐标位置，Ω为积分体，
Figure FDA00030929978500000312
为移动坐标位置；
将公式(10)带入到公式(8)中，得到二阶散射波表达式；重复上述过程即可得到多阶散射波场表达式，用求和形式简单表达N阶散射波场：
Figure FDA00030929978500000313
其中，us为散射波场，
Figure FDA00030929978500000314
为第i阶散射波场，N为总阶数，
Figure FDA00030929978500000315
为坐标位置，ω为角频率。
5.根据权利要求1所述的基于线性伯恩近似约束的非线性最小二乘逆时偏移方法，其特征在于，在步骤2中，当采用一阶Born近似下的散射波方程作为正演算子时，L-2范数下的LSRTM目标泛函表示为：
Figure FDA0003092997850000041
其中，J为目标泛函，
Figure FDA0003092997850000042
为一阶模拟散射波场，
Figure FDA0003092997850000043
为一阶散射波场，dn为第n阶背景波场；
其对应的速度扰动更新的梯度表示为：
Figure FDA0003092997850000044
其中，
Figure FDA0003092997850000045
为速度扰动，J为目标泛函，
Figure FDA0003092997850000046
为坐标位置，ω为角频率，
Figure FDA0003092997850000047
为一阶模拟散射波场，
Figure FDA0003092997850000048
为一阶模拟波场，dn为第n阶背景波场，u0为背景波场值，V为散射势，G为格林函数，
Figure FDA0003092997850000049
为检波点坐标位置；
分析可知，上述梯度包含两部分：模拟背景波场与模拟的一阶散射波场与观测数据中的一阶散射波场之间的残差所构成的梯度：
Figure FDA00030929978500000410
其中，
Figure FDA00030929978500000411
为一阶梯度，
Figure FDA00030929978500000412
为坐标位置，ω为角频率，u0为背景波场值，G为格林函数，
Figure FDA00030929978500000413
为检波点坐标位置，
Figure FDA00030929978500000414
为移动坐标位置，
Figure FDA00030929978500000415
为一阶模拟散射波场，
Figure FDA00030929978500000416
为一阶模拟波场；
模拟背景波场与观测数据中的背景波场和高阶Born近似下的观测数据所构成的梯度：
Figure FDA00030929978500000417
其中，gn为第n阶梯度，ω为角频率，u0为背景波场值，
Figure FDA0003092997850000051
为坐标位置，
Figure FDA0003092997850000052
为检波点坐标位置，dn为第n阶背景波场，G为格林函数，
Figure FDA0003092997850000053
为移动坐标位置；
梯度第一项描述了正演算子所描述的一阶散射波对应的速度扰动，而第二项在一阶散射波假设下表现为成像噪音，即在一阶Born近似下，观测数据的其他分量由于无法被准确描述，而与背景波场产生了对速度扰动梯度的干扰。
6.根据权利要求5所述的基于线性伯恩近似约束的非线性最小二乘逆时偏移方法，其特征在于，在步骤2中，为了压制观测数据中背景波场与高阶散射波场产生的噪音，将线性伯恩近似方程作为约束项加入到非线性LSRTM的目标函数中，加强非线性LSRTM中一阶散射波的约束，以获得高精度，高分辨率的剖面。
7.根据权利要求6所述的基于线性伯恩近似约束的非线性最小二乘逆时偏移方法，其特征在于，在步骤2中，改进后的LSRTM目标函数表示为：
Figure FDA0003092997850000054
其中，λ为约束因子，根据观测数据中的背景波场和高阶散射波场的强度进行调整，J为目标泛函，
Figure FDA0003092997850000055
为坐标位置，ω为角频率，
Figure FDA0003092997850000056
为一阶模拟散射波场，
Figure FDA0003092997850000057
为一阶模拟波场，dn为第n阶背景波场，u0为背景波场值，V为散射势，G为格林函数,
Figure FDA0003092997850000058
为检波点坐标位置，
Figure FDA0003092997850000059
为移动坐标位置；
添加约束后的LSRTM分解为两个反演问题：首先针对一阶散射波场进行反演，得到相对准确的一阶散射波场；然后在此基础上计算速度扰动的更新梯度；令公式(16)对于一阶散射波场的导数为0，可以得到一阶散射波场的重构表达式：
Figure FDA0003092997850000061
其中，λ为约束因子，根据观测数据中的背景波场和高阶散射波场的强度进行调整，J为目标泛函，
Figure FDA0003092997850000062
为坐标位置，ω为角频率，
Figure FDA0003092997850000063
为一阶模拟散射波场，
Figure FDA0003092997850000064
为一阶模拟波场，dn为第n阶背景波场，u0为背景波场值，V为散射势，G为格林函数,
Figure FDA0003092997850000065
为检波点坐标位置,
Figure FDA0003092997850000066
为移动坐标位置；
Figure FDA0003092997850000067
其中，λ为约束因子，
Figure FDA0003092997850000068
为坐标位置，ω为角频率，
Figure FDA0003092997850000069
为一阶模拟散射波场，
Figure FDA00030929978500000610
为一阶模拟波场，dn为第n阶背景波场，u0为背景波场值，V为散射势，G为格林函数,
Figure FDA00030929978500000611
为移动坐标位置；
分析式(18)可看出，重构的一阶散射波场主要由两部分组成：观测数据构成的波场和模拟的一阶散射波场，两者的比重由约束因子控制：当约束因子较大时，一阶散射波对重构散射波场的影响变大，从而减弱了背景波场和高级散射波对一阶重构散射波的影响。
8.根据权利要求7所述的基于线性伯恩近似约束的非线性最小二乘逆时偏移方法，其特征在于，在步骤2中，在准确重构散射波场之后，得到目标函数对速度扰动的梯度：
Figure FDA00030929978500000612
其中，λ为约束因子，J为目标泛函，
Figure FDA0003092997850000071
为坐标位置，ω为角频率，u0为背景波场值，V为散射势，
Figure FDA0003092997850000072
为移动坐标位置，G为格林函数，
Figure FDA0003092997850000073
为速度扰动，
Figure FDA0003092997850000074
为一阶模拟散射波场，u为波场值，u0为背景波场值；
进而对速度扰动进行更新：
Figure FDA0003092997850000075
其中，α为更新步长，J为目标泛函，
Figure FDA0003092997850000076
为速度扰动。
9.根据权利要求1所述的基于线性伯恩近似约束的非线性最小二乘逆时偏移方法，其特征在于，在步骤3中，针对一般的反演问题，其目标函数对与反演参数的二阶导数即Hessian矩阵的性质决定了反演问题的凸性：当Hessian矩阵时，Hessian矩阵为半正定的，即目标函数是全凸的，通过局部寻优法寻找到目标函数的全局极值；反之，即目标函数是局部凸的，反演问题受局部极值的影响，难以收敛到全局极值。
10.根据权利要求9所述的基于线性伯恩近似约束的非线性最小二乘逆时偏移方法，其特征在于，在步骤3中，基于线性伯恩近似约束的非线性LSRTM对一阶散射波场的Hessian矩阵为：
Figure FDA0003092997850000077
其中，λ为约束因子，J为目标泛函，
Figure FDA0003092997850000078
为坐标位置，
Figure FDA0003092997850000079
为坐标位置，ω为角频率，u0为背景波场值，V为散射势，G为格林函数，
Figure FDA00030929978500000710
为移动坐标位置，
Figure FDA00030929978500000711
为一阶模拟散射波场，
Figure FDA00030929978500000712
为一阶模拟波场，dn为第n阶背景波场，
Figure FDA00030929978500000713
为检波点坐标位置；
由一阶散射波场的表达式可看出，波场受模型参数、散射势及观测系统影响；不同观测系统下的波场是不相关的，则：
Figure FDA0003092997850000081
其中，λ为约束因子，J为目标泛函，
Figure FDA0003092997850000082
为坐标位置，
Figure FDA0003092997850000083
为坐标位置，ω为角频率，
Figure FDA0003092997850000084
为一阶模拟散射波场；
即基于线性伯恩近似约束的非线性LSRTM针对一阶散射波场表现出强凸性，用局部寻优的方法得到较为准确的一阶散射波场；此外，Born近似下的一阶散射波场与速度扰动之间是线性的，因此基于重构的一阶散射波场得到准确的速度扰动更新梯度。
CN202110602166.0A 2021-05-31 2021-05-31 基于线性伯恩近似约束的非线性最小二乘逆时偏移方法 Active CN115480294B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110602166.0A CN115480294B (zh) 2021-05-31 2021-05-31 基于线性伯恩近似约束的非线性最小二乘逆时偏移方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110602166.0A CN115480294B (zh) 2021-05-31 2021-05-31 基于线性伯恩近似约束的非线性最小二乘逆时偏移方法

Publications (2)

Publication Number Publication Date
CN115480294A true CN115480294A (zh) 2022-12-16
CN115480294B CN115480294B (zh) 2024-12-10

Family

ID=84419514

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110602166.0A Active CN115480294B (zh) 2021-05-31 2021-05-31 基于线性伯恩近似约束的非线性最小二乘逆时偏移方法

Country Status (1)

Country Link
CN (1) CN115480294B (zh)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20170115418A1 (en) * 2014-06-30 2017-04-27 Cgg Services Sas Seismic data processing using matching filter based cost function optimization
CN107783190A (zh) * 2017-10-18 2018-03-09 中国石油大学(北京) 一种最小二乘逆时偏移梯度更新方法
KR20210042720A (ko) * 2019-10-10 2021-04-20 인하대학교 산학협력단 속도 모델을 위한 탄성파 파형역산 장치 및 방법

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20170115418A1 (en) * 2014-06-30 2017-04-27 Cgg Services Sas Seismic data processing using matching filter based cost function optimization
CN107783190A (zh) * 2017-10-18 2018-03-09 中国石油大学(北京) 一种最小二乘逆时偏移梯度更新方法
KR20210042720A (ko) * 2019-10-10 2021-04-20 인하대학교 산학협력단 속도 모델을 위한 탄성파 파형역산 장치 및 방법

Non-Patent Citations (9)

* Cited by examiner, † Cited by third party
Title
KE CHEN 等: "Time-domain elastic Gauss–Newton full-waveform inversion: a matrix-free approach", GEOPHYSICAL JOURNAL INTERNATIONAL, vol. 223, no. 2, 21 July 2020 (2020-07-21), pages 1007 - 1039 *
任志明: "声波和弹性波波动方程有限差分正反演方法研究", 中国博士学位论文全文数据库 基础科学辑, no. 2, 15 February 2018 (2018-02-15), pages 011 - 80 *
刘玉金 等: "扩展成像条件下的最小二乘逆时偏移", 地球物理学报, vol. 58, no. 10, 31 October 2015 (2015-10-31), pages 3771 - 3782 *
刘玉金: "复杂介质条件下基于反演的地震成像方法研究", 中国博士学位论文全文数据库 基础科学辑, no. 1, 15 January 2018 (2018-01-15), pages 011 - 129 *
孙小东 等: "基于共轭梯度法和互相关的最小二乘逆时偏移及应用(英文)", 应用地球物理(英文版), vol. 14, no. 003, 30 September 2017 (2017-09-30), pages 381 - 386 *
徐凯: "基于 GPU/CPU 协同加速的最小二乘逆时偏移", 中国优秀硕士学位论文全文数据库 基础科学辑, no. 4, 15 April 2018 (2018-04-15), pages 011 - 563 *
方修政 等: "基于逆散射成像条件的最小二乘逆时偏移", 地球物理学报, vol. 61, no. 09, 30 September 2018 (2018-09-30), pages 3770 - 3782 *
李武群 等: "二阶Born近似及其在逆散射GRT反演中的应用", 地球物理学进展, vol. 32, no. 02, 30 April 2017 (2017-04-30), pages 645 - 656 *
王华忠 等: "最小二乘叠前深度偏移成像理论与方法", 石油物探, vol. 56, no. 02, 31 March 2017 (2017-03-31), pages 159 - 170 *

Also Published As

Publication number Publication date
CN115480294B (zh) 2024-12-10

Similar Documents

Publication Publication Date Title
CN108873066B (zh) 弹性介质波动方程反射波旅行时反演方法
CN110187382B (zh) 一种回折波和反射波波动方程旅行时反演方法
CA2972033C (en) Multistage full wavefield inversion process that generates a multiple free data set
CN107894618B (zh) 一种基于模型平滑算法的全波形反演梯度预处理方法
CN109557582B (zh) 一种二维多分量地震资料偏移成像方法及系统
US11635540B2 (en) Methods and devices performing adaptive quadratic Wasserstein full-waveform inversion
CN111290016A (zh) 一种基于地质模型约束的全波形速度建模反演方法
Feng et al. New dynamic stochastic source encoding combined with a minmax-concave total variation regularization strategy for full waveform inversion
CN104730572B (zh) 一种基于l0半范数的绕射波成像方法及装置
Zhang et al. Simultaneous interval-Q estimation and attenuation compensation based on fusion deep neural network
CN114966831A (zh) 一种基于速度-衰减解耦的粘声全波形反演方法
CN115480294A (zh) 基于线性伯恩近似约束的非线性最小二乘逆时偏移方法
Shi et al. Prestack correlative elastic least-squares reverse time migration based on wavefield decomposition
Holliger et al. Sensitivity of the lateral correlation function in deep seismic reflection data
Hadden et al. Anisotropic full-waveform inversion of crosshole seismic data: A vertical symmetry axis field data application
Shustak et al. Q-factor inversion using reconstructed source spectral consistency
CN111025393B (zh) 针对含薄煤层地层的储层预测方法、装置、设备及介质
Ernst et al. Reduction of near-surface scattering effects in seismic data
Wang et al. Rock Fracture Monitoring Based on High‐Precision Microseismic Event Location Using 3D Multiscale Waveform Inversion
Huang et al. A fast least-squares reverse time migration method using cycle-consistent generative adversarial network
Kotsi Time-lapse seismic imaging and uncertainty quantification
Li et al. A characterization method for cavity Karst reservoir using local full-waveform inversion in frequency domain
Hu et al. Waveform energy focusing tomography with passive seismic sources
Chen Accurate background velocity model building method based on iterative deep learning in sparse transform domain
Liu Direct Waveform Inversion Using Explicit Time-Space Causality Principle

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