CN111580163B - 一种基于非单调搜索技术的全波形反演方法及*** - Google Patents

一种基于非单调搜索技术的全波形反演方法及*** Download PDF

Info

Publication number
CN111580163B
CN111580163B CN202010466886.4A CN202010466886A CN111580163B CN 111580163 B CN111580163 B CN 111580163B CN 202010466886 A CN202010466886 A CN 202010466886A CN 111580163 B CN111580163 B CN 111580163B
Authority
CN
China
Prior art keywords
iteration
equal
objective function
monotonic
seismic
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
CN202010466886.4A
Other languages
English (en)
Other versions
CN111580163A (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.)
Institute of Geology and Geophysics of CAS
Original Assignee
Institute of Geology and Geophysics of CAS
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 Institute of Geology and Geophysics of CAS filed Critical Institute of Geology and Geophysics of CAS
Priority to CN202010466886.4A priority Critical patent/CN111580163B/zh
Publication of CN111580163A publication Critical patent/CN111580163A/zh
Application granted granted Critical
Publication of CN111580163B publication Critical patent/CN111580163B/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/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/16—Receiving elements for seismic signals; Arrangements or adaptations of receiving elements
    • G01V1/18—Receiving elements, e.g. seismometer, geophone or torque detectors, for localised single point measurements
    • G01V1/181—Geophones
    • 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

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)
  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)

Abstract

本发明的目的是提供一种基于非单调搜索技术的全波形反演方法及***，方法包括：首先按照设定频组范围将多尺度反演频率划分N个频组；其次采用非单调搜索技术结合初始步长，根据所述地震观测数据和初始速度模型确定第j频组对应的下次迭代的速度模型；然后判断j是否大于N；如果j大于N，则将第j频组对应的下次迭代的速度模型作为最优速度模型输出；最后基于所述最优速度模型进行地震成像。本发明将非单调搜索技术与初始步长相结合确定步长，该方法既能够在保证反演计算效率的前提下，获得高精度的反演结果；又能够在初始速度模型较差的时候，一定程度上保证反演的全局收敛性，获得较好的反演结果。

Description

一种基于非单调搜索技术的全波形反演方法及***
技术领域
本发明涉及地震勘探技术领域，特别是涉及一种基于非单调搜索技术的全波形反演方法及***。
背景技术
全波形反演(Full waveform inversion，FWI)综合利用地震波场中的全波信息，通过拟合实际观测波场和预测波场来定量修正初始模型，以重建地下介质的速度和密度等参数模型，可以为地震勘探和地球内部结构探索提供高精度的纵横波速度和密度等丰富的弹性参数信息。
全波形反演总体上是一个非线性数据拟合过程，这个过程以迭代的方式重复进行，直到观测数据与合成地震数据构成的目标函数足够小为止。首先需要给出一个地下地层参数的初始估计值，求取目标函数及其梯度，并结合反演优化算法构造搜索方向，计算步长，通过不断更新模型参数使目标函数逐次减小，最终获得满意的结果。对于非线性的地震波全波形反演问题，迭代步长必须采用沿着搜索方向的线性搜索方法求取。线性搜索方法主要包括非精确线性搜索和精确线性搜索两大类方法。精确线性搜索方法是从目标函数出发进行求解，获得最优的精确步长，可以说不同的目标函数对应不同的精确步长计算方法。抛物线类方法和直接求解法是广泛使用的精确线性搜索方法，但两者都是基于二范数型目标函数进行求解获得。针对Huber范数和Hybrid范数也给出了混合范数相适应的精确步长计算方法。虽然精确线性搜索方法在反演精度和收敛速度上有一定优势，但其需要针对不同的目标函数计算不同的精确步长，通用性较差，如互相关目标函数、包络型目标函数等目前具有较好发展前景的目标函数，尚未有研究给出其精确步长计算方法。
相对于精确线性搜索方法，非精确线性搜索方法却对目标函数具有广泛的适用性。非精确线性搜索方法主要利用判断条件和初始步长来确保每次迭代过程中目标函数在搜索方向上有足够的下降量，其总体上可以使目标函数下降快，在每一步并不要求目标函数达到精确最小，仅需要合适的下降量即可。目前，在地震全波形反演中传统的非精确线性搜索方法通常使用单调搜索技术要求目标函数逐次下降，如经典的Wolfe判断条件。虽然这类方法也能使反演收敛，但是迫使目标函数逐次下降，会破坏反演优化方法的全局收敛性质，并且当反演陷入局部极小值附近，将会出现小步长和之字形现象，降低反演的计算效率和精度，目标函数逐次下降不利于保护反演的全局收敛性。
发明内容
基于此，本发明的目的是提供一种基于非单调搜索技术的全波形反演方法及***，将非单调搜索技术与初始步长相结合确定步长，以提高全波形反演的计算效率和精度。
为实现上述目的，本发明提供了一种基于非单调搜索技术的全波形反演方法，所述方法包括：
步骤S1：利用检波器采集地震观测数据；
步骤S2：按照设定频组范围将多尺度反演频率划分N个频组；其中，N为大于等于1的正整数；
步骤S3：采用非单调搜索技术结合初始步长，根据所述地震观测数据和初始速度模型确定第j频组对应的下次迭代的速度模型；
步骤S4：判断j是否大于N；如果j大于N，则将第j频组对应的下次迭代的速度模型作为最优速度模型输出；如果j小于或等于N，则令j＝j+1，将第j频组对应的下次迭代的速度模型作为初始速度模型，并返回步骤S3；
步骤S5：基于所述最优速度模型进行地震成像。
可选的，所述采用非单调搜索技术结合初始步长，根据所述地震观测数据和初始速度模型确定第j频组对应的下次迭代的速度模型，具体包括：
步骤S31：采用有限差分法，根据初始速度模型求解声波方程，计算第j频组第k次迭代的模拟波场数据；
步骤S32：根据所述地震观测数据和所述模拟波场数据构建目标函数值；
步骤S33：利用伴随状态法，根据所述模拟波场数据计算梯度；
步骤S34：利用有限内存拟牛顿法，根据所述梯度构造搜索方向；
步骤S35：采用非单调搜索技术结合初始步长，根据所述目标函数值、所述搜索方向和所述梯度确定地震全波形反演的步长；
步骤S36：根据所述搜索方向和所述步长确定第j频组对应的下次迭代的速度模型；
步骤S37：判断迭代次数是否大于或等于迭代次数阈值；如果迭代次数k小于迭代次数阈值，则令k＝k+1，并返回“步骤S31”；如果迭代次数k大于或等于迭代次数阈值，则执行“步骤S4”；
和/或，计算连续目标函数相对差值，并判断所述连续目标函数相对差值是否小于第一设定系数；如果所述连续目标函数相对差值大于或等于设定系数，则令k＝k+1，返回“步骤S31”；如果所述连续目标函数相对差值小于设定系数，则执行“步骤S4”。
可选的，所述采用非单调搜索技术结合初始步长，根据所述目标函数值、所述搜索方向和所述梯度确定地震全波形反演的步长，具体公式为：
Figure BDA0002512960350000031
其中，mk+1为第k+1次迭代的速度模型，αk为第k次迭代的步长，dk为第k次迭代的搜索方向，c1为常数，Ck+1＝(ηkQkCk+f(mk+1))/Qk+1,Qk+1＝ηkQk+1，ηk为非线性程度，取值为0或1，Qk为所选用的连续目标函数的数量，Ck为Qk次目标函数的平均值，f(mk+1)为k+1次迭代的目标函数值，gk表示第k次迭代的梯度，
Figure BDA0002512960350000032
表示第k次迭代的初始步长。
可选的，所述利用伴随状态法，根据所述模拟波场数据计算梯度，具体公式为：
Figure BDA0002512960350000033
其中，gk表示第k次迭代的梯度，mk为第k次迭代的速度模型，T表示观测时间长度，λk表示第k次迭代的反传残差波场，pk(xr,t|xs)表示t时刻由位于xs处的震源在xr处计算的模拟波场数据。
可选的，所述根据所述地震观测数据和模拟波场数据构建目标函数值，具体公式为：
Figure BDA0002512960350000041
其中，f(mk)表示第k次迭代的目标函数值，mk为第k次迭代的速度模型，xs表示第s个震源位置，xr表示第r个检波点位置，t表示时间，pk(xr,t|xs)表示t时刻由位于xs处的震源在xr处计算的模拟波场数据，
Figure BDA0002512960350000042
表示t时刻由位于xs处的震源在xr处的观测波场数据。
本发明还提供一种基于非单调搜索技术的全波形反演***，所述***包括：
采集模块，用于利用检波器采集地震观测数据；
分组模块，用于按照设定频组范围将多尺度反演频率划分N个频组；其中，N为大于等于1的正整数；
速度模型确定模块，用于采用非单调搜索技术结合初始步长，根据所述地震观测数据和初始速度模型确定第j频组对应的下次迭代的速度模型；
判断模块，用于判断j是否大于N；如果j大于N，则将第j频组对应的下次迭代的速度模型作为最优速度模型输出；如果j小于或等于N，则令j＝j+1，将第j频组对应的下次迭代的速度模型作为初始速度模型，并返回步骤S3；
地震成像模块，用于基于所述最优速度模型进行地震成像。
可选的，所述速度模型确定模块，具体包括：
模拟波场数据确定单元，用于采用有限差分法，根据初始速度模型求解声波方程，计算第j频组第k次迭代的模拟波场数据；
目标函数值构建单元，用于根据所述地震观测数据和所述模拟波场数据构建目标函数值；
梯度计算单元，用于利用伴随状态法，根据所述模拟波场数据计算梯度；
搜索方向构造单元，用于利用有限内存拟牛顿法，根据所述梯度构造搜索方向；
步长确定单元，用于采用非单调搜索技术结合初始步长，根据所述目标函数值、所述搜索方向和所述梯度确定地震全波形反演的步长；
速度模型确定单元，用于根据所述搜索方向和所述步长确定第j频组对应的下次迭代的速度模型；
判断单元，用于判断迭代次数是否大于或等于迭代次数阈值；如果迭代次数k小于迭代次数阈值，则令k＝k+1，并返回“模拟波场数据确定单元”；如果迭代次数k大于或等于迭代次数阈值，则执行“判断模块”；
和/或，计算连续目标函数相对差值，并判断所述连续目标函数相对差值是否小于第一设定系数；如果所述连续目标函数相对差值大于或等于设定系数，则令k＝k+1，返回“模拟波场数据确定单元”；如果所述连续目标函数相对差值小于设定系数，则执行“判断模块”。
可选的，所述采用非单调搜索技术结合初始步长，根据所述目标函数值、所述搜索方向和所述梯度确定地震全波形反演的步长，具体公式为：
Figure BDA0002512960350000051
其中，mk+1为第k+1次迭代的速度模型，αk为第k次迭代的步长，dk为第k次迭代的搜索方向，c1为常数，Ck+1＝(ηkQkCk+f(mk+1))/Qk+1,Qk+1＝ηkQk+1，ηk为非线性程度，取值为0或1，Qk为所选用的连续目标函数的数量，Ck为Qk次目标函数的平均值，f(mk+1)为k+1次迭代的目标函数值，gk表示第k次迭代的梯度，
Figure BDA0002512960350000052
表示第k次迭代的初始步长。
可选的，所述利用伴随状态法，根据所述模拟波场数据计算梯度，具体公式为：
Figure BDA0002512960350000053
其中，gk表示第k次迭代的梯度，mk为第k次迭代的速度模型，T表示观测时间长度，λk表示第k次迭代的反传残差波场，pk(xr,t|xs)表示t时刻由位于xs处的震源在xr处计算的模拟波场数据。
可选的，所述根据所述地震观测数据和模拟波场数据构建目标函数值，具体公式为：
Figure BDA0002512960350000054
其中，f(mk)表示第k次迭代的目标函数值，mk为第k次迭代的速度模型，xs表示第s个震源位置，xr表示第r个检波点位置，t表示时间，pk(xr,t|xs)表示t时刻由位于xs处的震源在xr处计算的模拟波场数据，
Figure BDA0002512960350000061
表示t时刻由位于xs处的震源在xr处的观测波场数据。
根据本发明提供的具体实施例，本发明公开了以下技术效果：
本发明的目的是提供一种基于非单调搜索技术的全波形反演方法及***，方法包括：首先按照设定频组范围将多尺度反演频率划分N个频组；其次采用非单调搜索技术结合初始步长，根据所述地震观测数据和初始速度模型确定第j频组对应的下次迭代的速度模型；然后判断j是否大于N；如果j大于N，则将第j频组对应的下次迭代的速度模型作为最优速度模型输出；最后基于所述最优速度模型进行地震成像。本发明将非单调搜索技术与初始步长相结合确定步长，该方法既能够在保证反演计算效率的前提下，获得高精度的反演结果；又能够在初始速度模型较差的时候，一定程度上保证反演的全局收敛性，获得较好的反演结果。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。
图1为本发明实施例基于非单调搜索技术的全波形反演方法流程图；
图2为本发明实施例基于非单调搜索技术的全波形反演***结构图；
图3为本发明实施例Marmousi模型示意图；
图4为本发明实施例光滑模型示意图；
图5为本发明实施例第一反演结果示意图；
图6为本发明实施例Overthrust模型示意图；
图7为本发明实施例均匀模型示意图；
图8为本发明实施例第二反演结果示意图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明的目的是提供一种基于非单调搜索技术的全波形反演方法及***，将非单调搜索技术与初始步长相结合确定步长，以提高全波形反演的计算效率和精度。
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
图1为本发明实施例基于非单调搜索技术的全波形反演方法流程图，如图1所示，本发明公开一种基于非单调搜索技术的全波形反演方法，所述方法包括：
步骤S1：利用检波器采集地震观测数据。
步骤S2：按照设定频组范围将多尺度反演频率划分N个频组；其中，N为大于等于1的正整数。
步骤S3：采用非单调搜索技术结合初始步长，根据所述地震观测数据和初始速度模型确定第j频组对应的下次迭代的速度模型。
步骤S4：判断j是否大于N；如果j大于N，则将第j频组对应的下次迭代的速度模型作为最优速度模型输出；如果j小于或等于N，则令j＝j+1，将第j频组对应的下次迭代的速度模型作为初始速度模型，并返回步骤S3。
步骤S5：基于所述最优速度模型进行地震成像。
在步骤S2之前还包括：对所述地震观测数据进行去噪和振幅能量补偿处理。
下面对各个步骤进行详细论述：
步骤S3：采用非单调搜索技术结合初始步长，根据所述地震观测数据和初始速度模型确定第j频组对应的下次迭代的速度模型，具体包括：
步骤S31：采用有限差分法，根据初始速度模型求解声波方程，计算第j频组第k次迭代的模拟波场数据；所述模拟波场数据为正传波场。
步骤S32：根据所述地震观测数据和所述模拟波场数据构建目标函数值，具体公式为：
Figure BDA0002512960350000081
其中，f(mk)表示第k次迭代的目标函数值，mk为第k次迭代的速度模型，xs表示第s个震源位置，xr表示第r个检波点位置，t表示时间，pk(xr,t|xs)表示t时刻由位于xs处的震源在xr处计算的模拟波场数据，
Figure BDA0002512960350000082
表示t时刻由位于xs处的震源在xr处的观测波场数据。
步骤S33：利用伴随状态法，根据所述模拟波场数据计算梯度，具体公式为：
Figure BDA0002512960350000083
其中，gk表示第k次迭代的梯度，mk为第k次迭代的速度模型，T表示观测时间长度，λk表示第k次迭代的反传残差波场，pk(xr,t|xs)表示t时刻由位于xs处的震源在xr处计算的模拟波场数据。
步骤S34：利用有限内存拟牛顿法，根据所述梯度构造搜索方向，具体公式为：
dk＝-Hk·gk (3)；
其中，dk表示第k次迭代的搜索方向，gk表示第k次迭代的梯度，Hk表示第k次迭代的海森矩阵的逆矩阵。
步骤S35：采用非单调搜索技术结合初始步长，根据所述目标函数值、所述搜索方向和所述梯度确定地震全波形反演的步长，具体公式为：
Figure BDA0002512960350000084
其中，mk+1为第k+1次迭代的速度模型，αk为第k次迭代的步长，dk为第k次迭代的搜索方向，c1为常数，Ck+1＝(ηkQkCk+f(mk+1))/Qk+1,Qk+1＝ηkQk+1，ηk为非线性程度，取值为0或1，Qk为所选用的连续目标函数的数量，Ck为Qk次目标函数的平均值，f(mk+1)为k+1次迭代的目标函数值，gk表示第k次迭代的梯度，
Figure BDA0002512960350000091
表示第k次迭代的初始步长。
当ηk＝0时(4)式的判断条件退化为传统的单调判断条件，ηk＝1时Ck表示近Qk次迭代情况下目标函数的平均值。当然，除了可以选取连续迭代情况下一定数量的目标函数平均值进行非单调条件判断，还可以尝试这些目标函数的最大值作为Ck。非单调搜索技术不要求目标函数单调下降，允许在反演过程中出现目标函数局部上升的现象，不仅充分利用了迭代过程中产生的相关目标函数信息，也在一定程度上保护了反演的全局收敛性质。
步骤S36：根据所述搜索方向和所述步长确定第j频组对应的下次迭代的速度模型，具体公式为：
mk+1＝mk+αkdk (5)；
其中，mk+1为第k+1次迭代的速度模型，mk为第k次迭代的速度模型，αk为第k次迭代的步长，dk为第k次迭代的搜索方向。
步骤S37：判断迭代次数是否大于或等于迭代次数阈值；如果迭代次数k小于迭代次数阈值，则令k＝k+1，并返回“步骤S31”；如果迭代次数k大于或等于迭代次数阈值，则执行“步骤S4”。
或，步骤S37：计算连续目标函数相对差值，并判断所述连续目标函数相对差值是否小于第一设定系数；如果所述连续目标函数相对差值大于或等于设定系数，则令k＝k+1，返回“步骤S31”；如果所述连续目标函数相对差值小于设定系数，则执行“步骤S4”。
或，步骤S37：判断迭代次数是否大于或等于迭代次数阈值；如果迭代次数k小于迭代次数阈值，则令k＝k+1，并返回“步骤S31”；如果迭代次数k大于或等于迭代次数阈值，则执行“步骤S38”。
步骤S38：计算连续目标函数相对差值，并判断所述连续目标函数相对差值是否小于第一设定系数；如果所述连续目标函数相对差值大于或等于设定系数，则令k＝k+1，返回“步骤S31”；如果所述连续目标函数相对差值小于设定系数，则执行“步骤S4”。
步骤S4：判断j是否大于N；如果j大于N，则将下次迭代的速度模型作为最优速度模型输出；如果j小于或等于N，则令j＝j+1，将下次迭代的速度模型作为初始速度模型，并返回步骤S3。
步骤S5：基于所述最优速度模型进行地震成像。
如图2所示，本发明还提供一种基于非单调搜索技术的全波形反演***，所述***包括：
采集模块1，用于利用检波器采集地震观测数据。
分组模块2，用于按照设定频组范围将多尺度反演频率划分N个频组；其中，N为大于等于1的正整数。
速度模型确定模块3，用于采用非单调搜索技术结合初始步长，根据所述地震观测数据和初始速度模型确定第j频组对应的下次迭代的速度模型。
判断模块4，用于判断j是否大于N；如果j大于N，则将第j频组对应的下次迭代的速度模型作为最优速度模型输出；如果j小于或等于N，则令j＝j+1，将第j频组对应的下次迭代的速度模型作为初始速度模型，并返回步骤S3。
地震成像模块5，用于基于所述最优速度模型进行地震成像。
作为一种实施方式，本发明所述速度模型确定模块3，具体包括：
模拟波场数据确定单元，用于采用有限差分法，根据初始速度模型求解声波方程，计算第j频组第k次迭代的模拟波场数据。
目标函数值构建单元，用于根据所述地震观测数据和所述模拟波场数据构建目标函数值。
梯度计算单元，用于利用伴随状态法，根据所述模拟波场数据计算梯度。
搜索方向构造单元，用于利用有限内存拟牛顿法，根据所述梯度构造搜索方向。
步长确定单元，用于采用非单调搜索技术结合初始步长，根据所述目标函数值、所述搜索方向和所述梯度确定地震全波形反演的步长。
速度模型确定单元，用于根据所述搜索方向和所述步长确定第j频组对应的下次迭代的速度模型。
第一判断单元，用于判断迭代次数是否大于或等于迭代次数阈值；如果迭代次数k小于迭代次数阈值，则令k＝k+1，并返回“模拟波场数据确定单元”；如果迭代次数k大于或等于迭代次数阈值，则执行“判断模块”。
或，第二判断单元，用于计算连续目标函数相对差值，并判断所述连续目标函数相对差值是否小于第一设定系数；如果所述连续目标函数相对差值大于或等于设定系数，则令k＝k+1，返回“模拟波场数据确定单元”；如果所述连续目标函数相对差值小于设定系数，则执行“判断模块”。
或，第三判断单元，用于判断迭代次数是否大于或等于迭代次数阈值；如果迭代次数k小于迭代次数阈值，则令k＝k+1，并返回“模拟波场数据确定单元”；如果迭代次数k大于或等于迭代次数阈值，则执行“第四判断单元”；
第四判断单元，用于计算连续目标函数相对差值，并判断所述连续目标函数相对差值是否小于第一设定系数；如果所述连续目标函数相对差值大于或等于设定系数，则令k＝k+1，返回“模拟波场数据确定单元”；如果所述连续目标函数相对差值小于设定系数，则执行“判断模块”。
作为一种实施方式，本发明所述采用非单调搜索技术结合初始步长，根据所述目标函数值、所述搜索方向和所述梯度确定地震全波形反演的步长，具体公式为：
Figure BDA0002512960350000111
其中，mk+1为第k+1次迭代的速度模型，αk为第k次迭代的步长，dk为第k次迭代的搜索方向，c1为常数，Ck+1＝(ηkQkCk+f(mk+1))/Qk+1,Qk+1＝ηkQk+1，ηk为非线性程度，取值为0或1，Qk为所选用的连续目标函数的数量，Ck为Qk次目标函数的平均值，f(mk+1)为k+1次迭代的目标函数值，gk表示第k次迭代的梯度，
Figure BDA0002512960350000112
表示第k次迭代的初始步长。
作为一种实施方式，本发明所述利用伴随状态法，根据所述模拟波场数据计算梯度，具体公式为：
Figure BDA0002512960350000113
其中，gk表示第k次迭代的梯度，mk为第k次迭代的速度模型，T表示观测时间长度，λk表示第k次迭代的反传残差波场，pk(xr,t|xs)表示t时刻由位于xs处的震源在xr处计算的模拟波场数据。
作为一种实施方式，本发明所述根据所述地震观测数据和模拟波场数据构建目标函数值，具体公式为：
Figure BDA0002512960350000121
其中，f(mk)表示第k次迭代的目标函数值，mk为第k次迭代的速度模型，xs表示第s个震源位置，xr表示第r个检波点位置，t表示时间，pk(xr,t|xs)表示t时刻由位于xs处的震源在xr处计算的模拟波场数据，
Figure BDA0002512960350000122
表示t时刻由位于xs处的震源在xr处的观测波场数据。
本发明利用非单调搜索技术求解最优化问题不要求目标函数单调下降，允许在反演过程中出现目标函数局部上升的现象，充分利用了迭代过程中产生的相关目标函数信息。非单调搜索技术配合反演优化算法求解最优化问题，可以减少反演迭代过程中的目标函数和梯度的计算次数，提高计算效率。同时，在一定程度上可以保护反演的全局收敛性质。
图3为本发明实施例Marmousi模型示意图，图4为本发明实施例光滑模型示意图，本发明将Marmousi模型作为真实参数模型，将光滑模型作为初始速度模型进行反演，具体反演结果如图5所示，从图5中可以看出Marmousi模型的整体结构和精细构造都得到很好的反演。本试验共利用95次迭代以及105次额外目标函数计算，即可获得高精度的反演结果。该实例证明了本方法在保证反演计算效率的情况下，可以获得高精度的反演结果。
图6为本发明实施例Overthrust模型示意图，图7为本发明实施例均匀模型示意图，本发明将Overthrust模型作为真实参数模型，将均匀模型作为初始速度模型进行反演，具体反演结果如图8所示，从图8中可以看出非单调搜索技术依旧可以恢复重建Overthrust模型的中浅部的主要构造。这个例子从一定程度上说明了当初始模型较差的时候，非单调搜索技术可以保护全波形反演的全局收敛性。
本说明书中各个实施例采用递进的方式描述，每个实施例重点说明的都是与其他实施例的不同之处，各个实施例之间相同相似部分互相参见即可。
本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处。综上所述，本说明书内容不应理解为对本发明的限制。

Claims (6)

1.一种基于非单调搜索技术的全波形反演方法，其特征在于，所述方法包括：
步骤S1：利用检波器采集地震观测数据；
步骤S2：按照设定频组范围将多尺度反演频率划分N个频组；其中，N为大于等于1的正整数；
步骤S3：采用非单调搜索技术结合初始步长，根据所述地震观测数据和初始速度模型确定第j频组对应的下次迭代的速度模型；
步骤S4：判断j是否大于N；如果j大于N，则将第j频组对应的下次迭代的速度模型作为最优速度模型输出；如果j小于或等于N，则令j＝j+1，将第j频组对应的下次迭代的速度模型作为初始速度模型，并返回步骤S3；
步骤S5：基于所述最优速度模型进行地震成像；
采用非单调搜索技术结合初始步长，根据所述地震观测数据和初始速度模型确定第j频组对应的下次迭代的速度模型，具体包括：
步骤S31：采用有限差分法，根据初始速度模型求解声波方程，计算第j频组第k次迭代的模拟波场数据；
步骤S32：根据所述地震观测数据和所述模拟波场数据构建目标函数值；
步骤S33：利用伴随状态法，根据所述模拟波场数据计算梯度；
步骤S34：利用有限内存拟牛顿法，根据所述梯度构造搜索方向；
步骤S35：采用非单调搜索技术结合初始步长，根据所述目标函数值、所述搜索方向和所述梯度确定地震全波形反演的步长，具体公式为：
Figure FDA0003051353540000011
其中，mk+1为第k+1次迭代的速度模型，αk为第k次迭代的步长，dk为第k次迭代的搜索方向，c1为常数，Ck+1＝(ηkQkCk+f(mk+1))/Qk+1,Qk+1＝ηkQk+1，ηk为非线性程度，取值为0或1，Qk为所选用的连续目标函数的数量，Ck为Qk次目标函数的平均值，f(mk+1)为k+1次迭代的目标函数值，gk表示第k次迭代的梯度，
Figure FDA0003051353540000021
表示第k次迭代的初始步长；
步骤S36：根据所述搜索方向和所述步长确定第j频组对应的下次迭代的速度模型；
步骤S37：判断迭代次数是否大于或等于迭代次数阈值；如果迭代次数k小于迭代次数阈值，则令k＝k+1，并返回“步骤S31”；如果迭代次数k大于或等于迭代次数阈值，则执行“步骤S4”；
或，步骤S37：计算连续目标函数相对差值，并判断所述连续目标函数相对差值是否小于第一设定系数；如果所述连续目标函数相对差值大于或等于设定系数，则令k＝k+1，返回“步骤S31”；如果所述连续目标函数相对差值小于设定系数，则执行“步骤S4”；
或，步骤S37：判断迭代次数是否大于或等于迭代次数阈值；如果迭代次数k小于迭代次数阈值，则令k＝k+1，并返回“步骤S31”；如果迭代次数k大于或等于迭代次数阈值，则执行“步骤S38”；
步骤S38：计算连续目标函数相对差值，并判断所述连续目标函数相对差值是否小于第一设定系数；如果所述连续目标函数相对差值大于或等于设定系数，则令k＝k+1，返回“步骤S31”；如果所述连续目标函数相对差值小于设定系数，则执行“步骤S4”。
2.根据权利要求1所述的基于非单调搜索技术的全波形反演方法，其特征在于，所述利用伴随状态法，根据所述模拟波场数据计算梯度，具体公式为：
Figure FDA0003051353540000022
其中，gk表示第k次迭代的梯度，mk为第k次迭代的速度模型，T表示观测时间长度，λk表示第k次迭代的反传残差波场，pk(xr,t|xs)表示t时刻由位于xs处的震源在xr处计算的模拟波场数据。
3.根据权利要求1所述的基于非单调搜索技术的全波形反演方法，其特征在于，所述根据所述地震观测数据和模拟波场数据构建目标函数值，具体公式为：
Figure FDA0003051353540000031
其中，f(mk)表示第k次迭代的目标函数值，mk为第k次迭代的速度模型，xs表示第s个震源位置，xr表示第r个检波点位置，t表示时间，pk(xr,t|xs)表示t时刻由位于xs处的震源在xr处计算的模拟波场数据，
Figure FDA0003051353540000032
表示t时刻由位于xs处的震源在xr处的观测波场数据。
4.一种基于非单调搜索技术的全波形反演***，其特征在于，所述***包括：
采集模块，用于利用检波器采集地震观测数据；
分组模块，用于按照设定频组范围将多尺度反演频率划分N个频组；其中，N为大于等于1的正整数；
速度模型确定模块，用于采用非单调搜索技术结合初始步长，根据所述地震观测数据和初始速度模型确定第j频组对应的下次迭代的速度模型；
判断模块，用于判断j是否大于N；如果j大于N，则将第j频组对应的下次迭代的速度模型作为最优速度模型输出；如果j小于或等于N，则令j＝j+1，将第j频组对应的下次迭代的速度模型作为初始速度模型，并返回步骤S3；
地震成像模块，用于基于所述最优速度模型进行地震成像；
所述速度模型确定模块，具体包括：
模拟波场数据确定单元，用于采用有限差分法，根据初始速度模型求解声波方程，计算第j频组第k次迭代的模拟波场数据；
目标函数值构建单元，用于根据所述地震观测数据和所述模拟波场数据构建目标函数值；
梯度计算单元，用于利用伴随状态法，根据所述模拟波场数据计算梯度；
搜索方向构造单元，用于利用有限内存拟牛顿法，根据所述梯度构造搜索方向；
步长确定单元，用于采用非单调搜索技术结合初始步长，根据所述目标函数值、所述搜索方向和所述梯度确定地震全波形反演的步长；
速度模型确定单元，用于根据所述搜索方向和所述步长确定第j频组对应的下次迭代的速度模型；
第一判断单元，用于判断迭代次数是否大于或等于迭代次数阈值；如果迭代次数k小于迭代次数阈值，则令k＝k+1，并返回“模拟波场数据确定单元”；如果迭代次数k大于或等于迭代次数阈值，则执行“判断模块”；
或，第二判断单元，用于计算连续目标函数相对差值，并判断所述连续目标函数相对差值是否小于第一设定系数；如果所述连续目标函数相对差值大于或等于设定系数，则令k＝k+1，返回“模拟波场数据确定单元”；如果所述连续目标函数相对差值小于设定系数，则执行“判断模块”；
或，第三判断单元，用于判断迭代次数是否大于或等于迭代次数阈值；如果迭代次数k小于迭代次数阈值，则令k＝k+1，并返回“模拟波场数据确定单元”；如果迭代次数k大于或等于迭代次数阈值，则执行“第四判断单元”；
第四判断单元，用于计算连续目标函数相对差值，并判断所述连续目标函数相对差值是否小于第一设定系数；如果所述连续目标函数相对差值大于或等于设定系数，则令k＝k+1，返回“模拟波场数据确定单元”；如果所述连续目标函数相对差值小于设定系数，则执行“判断模块”；
所述采用非单调搜索技术结合初始步长，根据所述目标函数值、所述搜索方向和所述梯度确定地震全波形反演的步长，具体公式为：
Figure FDA0003051353540000041
其中，mk+1为第k+1次迭代的速度模型，αk为第k次迭代的步长，dk为第k次迭代的搜索方向，c1为常数，Ck+1＝(ηkQkCk+f(mk+1))/Qk+1,Qk+1＝ηkQk+1，ηk为非线性程度，取值为0或1，Qk为所选用的连续目标函数的数量，Ck为Qk次目标函数的平均值，f(mk+1)为k+1次迭代的目标函数值，gk表示第k次迭代的梯度，
Figure FDA0003051353540000051
表示第k次迭代的初始步长。
5.根据权利要求4所述的基于非单调搜索技术的全波形反演***，其特征在于，所述利用伴随状态法，根据所述模拟波场数据计算梯度，具体公式为：
Figure FDA0003051353540000052
其中，gk表示第k次迭代的梯度，mk为第k次迭代的速度模型，T表示观测时间长度，λk表示第k次迭代的反传残差波场，pk(xr,t|xs)表示t时刻由位于xs处的震源在xr处计算的模拟波场数据。
6.根据权利要求4所述的基于非单调搜索技术的全波形反演***，其特征在于，所述根据所述地震观测数据和模拟波场数据构建目标函数值，具体公式为：
Figure FDA0003051353540000053
其中，f(mk)表示第k次迭代的目标函数值，mk为第k次迭代的速度模型，xs表示第s个震源位置，xr表示第r个检波点位置，t表示时间，pk(xr,t|xs)表示t时刻由位于xs处的震源在xr处计算的模拟波场数据，
Figure FDA0003051353540000054
表示t时刻由位于xs处的震源在xr处的观测波场数据。
CN202010466886.4A 2020-05-28 2020-05-28 一种基于非单调搜索技术的全波形反演方法及*** Expired - Fee Related CN111580163B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010466886.4A CN111580163B (zh) 2020-05-28 2020-05-28 一种基于非单调搜索技术的全波形反演方法及***

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010466886.4A CN111580163B (zh) 2020-05-28 2020-05-28 一种基于非单调搜索技术的全波形反演方法及***

Publications (2)

Publication Number Publication Date
CN111580163A CN111580163A (zh) 2020-08-25
CN111580163B true CN111580163B (zh) 2021-08-10

Family

ID=72125676

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010466886.4A Expired - Fee Related CN111580163B (zh) 2020-05-28 2020-05-28 一种基于非单调搜索技术的全波形反演方法及***

Country Status (1)

Country Link
CN (1) CN111580163B (zh)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112084655A (zh) * 2020-09-08 2020-12-15 南京众诚土地规划设计咨询有限公司 一种基于非单调线搜索的探地雷达参数反演方法
CN113552620A (zh) * 2021-09-07 2021-10-26 中国地震局地球物理研究所 一种适用于波形反演的步长计算方法及***
CN115184986B (zh) * 2022-06-28 2023-09-29 吉林大学 不依赖震源的全局包络互相关全波形反演方法

Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107450102A (zh) * 2017-07-28 2017-12-08 西安交通大学 基于分辨率可控包络生成算子的多尺度全波形反演方法
CN108549100A (zh) * 2018-01-11 2018-09-18 吉林大学 基于非线性高次拓频的时间域多尺度全波形反演方法
CN110058307A (zh) * 2019-05-05 2019-07-26 四川省地质工程勘察院 一种基于快速拟牛顿法的全波形反演方法
CN110058302A (zh) * 2019-05-05 2019-07-26 四川省地质工程勘察院 一种基于预条件共轭梯度加速算法的全波形反演方法
WO2019186286A1 (en) * 2018-03-27 2019-10-03 King Abdullah University Of Science And Technology Robust full waveform inversion of seismic data method and device
WO2020009752A1 (en) * 2018-07-02 2020-01-09 Exxonmobil Upstream Research Company Full wavefield inversion with an image-gather-flatness constraint

Patent Citations (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107450102A (zh) * 2017-07-28 2017-12-08 西安交通大学 基于分辨率可控包络生成算子的多尺度全波形反演方法
CN108549100A (zh) * 2018-01-11 2018-09-18 吉林大学 基于非线性高次拓频的时间域多尺度全波形反演方法
WO2019186286A1 (en) * 2018-03-27 2019-10-03 King Abdullah University Of Science And Technology Robust full waveform inversion of seismic data method and device
WO2020009752A1 (en) * 2018-07-02 2020-01-09 Exxonmobil Upstream Research Company Full wavefield inversion with an image-gather-flatness constraint
CN110058307A (zh) * 2019-05-05 2019-07-26 四川省地质工程勘察院 一种基于快速拟牛顿法的全波形反演方法
CN110058302A (zh) * 2019-05-05 2019-07-26 四川省地质工程勘察院 一种基于预条件共轭梯度加速算法的全波形反演方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
基于初始模型优化的稀疏高效全波形反演方法;段超然;《中国博士学位论文全文数据库 基础科学辑》;20181215(第12期);第2,4章 *

Also Published As

Publication number Publication date
CN111580163A (zh) 2020-08-25

Similar Documents

Publication Publication Date Title
CN111580163B (zh) 一种基于非单调搜索技术的全波形反演方法及***
CN106526674B (zh) 一种三维全波形反演能量加权梯度预处理方法
CN103630933B (zh) 基于非线性优化的时空域交错网格有限差分方法和装置
CN107843925B (zh) 一种基于修正相位的反射波波形反演方法
CN107783185B (zh) 一种层析静校正的处理方法及装置
WO2023087451A1 (zh) 基于观测数据自编码的多尺度无监督地震波速反演方法
SG184803A1 (en) Artifact reduction in method of iterative inversion of geophysical data
CN110618453A (zh) 一种基于改进阻尼最小二乘法的波阻抗反演方法
CN103454677B (zh) 基于粒子群与线性加法器结合的地震数据反演方法
CN105319581A (zh) 一种高效的时间域全波形反演方法
CN104237937B (zh) 叠前地震反演方法及其***
CN113486591B (zh) 一种卷积神经网络结果的重力多参量数据密度加权反演方法
CN106483559A (zh) 一种地下速度模型的构建方法
CN110879412A (zh) 地下横波速度反演方法、装置、计算设备及存储介质
CN104965222A (zh) 三维纵波阻抗全波形反演方法及装置
CN111158059B (zh) 基于三次b样条函数的重力反演方法
CN109541691B (zh) 一种地震速度反演方法
CN109239776B (zh) 一种地震波传播正演模拟方法和装置
CN104730572A (zh) 一种基于l0半范数的绕射波成像方法及装置
CN108508481B (zh) 一种纵波转换波地震数据时间匹配的方法、装置及***
CN117331119A (zh) 一种面向隧道探测的快速地震波走时计算方法
CN111273346B (zh) 去除沉积背景的方法、装置、计算机设备及可读存储介质
CN111208568B (zh) 一种时间域多尺度全波形反演方法及***
Wang et al. Memory optimization in RNN-based full waveform inversion using boundary saving wavefield reconstruction
CN115421195A (zh) 地震探测中速度场的生成方法、装置、设备及存储介质

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: 20210810

CF01 Termination of patent right due to non-payment of annual fee