CN107783190A - 一种最小二乘逆时偏移梯度更新方法 - Google Patents

一种最小二乘逆时偏移梯度更新方法 Download PDF

Info

Publication number
CN107783190A
CN107783190A CN201710969635.6A CN201710969635A CN107783190A CN 107783190 A CN107783190 A CN 107783190A CN 201710969635 A CN201710969635 A CN 201710969635A CN 107783190 A CN107783190 A CN 107783190A
Authority
CN
China
Prior art keywords
mrow
msub
gradient
iteration
wave field
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
CN201710969635.6A
Other languages
English (en)
Other versions
CN107783190B (zh
Inventor
方修政
吴迪
钮凤林
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China University of Petroleum Beijing
Original Assignee
China University of Petroleum Beijing
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by China University of Petroleum Beijing filed Critical China University of Petroleum Beijing
Priority to CN201710969635.6A priority Critical patent/CN107783190B/zh
Publication of CN107783190A publication Critical patent/CN107783190A/zh
Application granted granted Critical
Publication of CN107783190B publication Critical patent/CN107783190B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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. analysis, for interpretation, for correction
    • 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
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/60—Analysis
    • G01V2210/62—Physical property of subsurface
    • G01V2210/624—Reservoir parameters
    • 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/64—Geostructures, e.g. in 3D data cubes
    • 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/65—Source localisation, e.g. faults, hypocenters or reservoirs

Abstract

本发明涉及一种最小二乘逆时偏移梯度更新方法，其步骤：根据初始速度模型建立二范数意义下波场残差目标泛函；正演模拟声波方程，存储震源边界波场和最后时刻的全部区域波场；逆时重建震源波场，反传残差波场获得伴随状态变量，利用逆散射条件更新单炮目标泛函相对于模型参数的梯度，判断是否满足多炮循环条件，如果是，重复执行并累加单炮目标泛函梯度，如果否，利用目标泛函梯度计算共轭梯度方向；正演模拟震源波场和共轭梯度波场数据，判断是否满足多炮循环条件，如果是，重复并累加每炮共轭梯度波场数据，如果否，计算迭代步长；更新偏移剖面，迭代更新；判断是否满足迭代条件，如果否，终止迭代并输出迭代剖面，如果是则返回正演模拟存储边界波场。

Description

一种最小二乘逆时偏移梯度更新方法
技术领域
本发明涉及一种石油地震勘探数据高效率、低存储、高精度成像领域，特别是关于一种最小二乘逆时偏移梯度更新方法。
背景技术
随着勘探目标越来越复杂，对于成像结果的精确性、分辨率和保幅性也提出了更高的要求。常规的一些偏移方法，例如逆时偏移、Kirchhoff偏移等，其精确性通常会受到地震数据噪音、不规则采样、有限偏移距等影响。
最小二乘偏移(Least Square Migration，LSM)作为一个反演问题不断的迭代拟合正演模拟的地震记录和实际观测的地震记录，隐含考虑了Hessian矩阵对振幅的影响，可以减少常规偏移过程中产生的假象，具有更高的成像精度以及保幅性等优势，可以获得更加精确的反演目标。Tarantola(1984)理论推导了二范数意义下波场残差最小二乘目标函数，通过目标函数的下降，不断的迭代反演成像剖面，证明了常规的偏移结果只是最小二乘偏移的第一步迭代。Nemeth et al.(1999)对不完整数据进行了最小二乘Kirchhoff偏移并取得了很好的效果，进一步展示了最小二乘偏移的优势。目前，最小二乘逆时偏移(LeastSquare Reverse Time Migration，LSRTM)的研究也取得了很大进展。但是，当速度体存在强反射界面或者地震记录含有大量折射波时，基于伴随状态法求取的目标函数梯度会产生很强的低频噪音，在迭代过程中会干扰成像精度、降低收敛速度。
发明内容
针对伴随状态法求取目标函数梯度，收敛速度和计算效率受低频噪音影响的问题，本发明的目的是提供一种最小二乘逆时偏移梯度更新方法，该方法可以加快目标泛函梯度收敛速度，提高偏移剖面的更新效率。
为实现上述目的，本发明采取以下技术方案：一种最小二乘逆时偏移梯度更新方法，其特征在于包括以下步骤：1)根据真实速度模型得到初始速度模型，建立二范数意义下波场残差目标泛函；2)采用基于CUDA计算平台，CPU/GPU协同并行加速方法正演模拟声波方程，存储震源边界波场和最后时刻的全部区域波场；3)利用存储的边界波场和最后时刻的全部区域波场逆时重建震源波场，利用伴随状态法逆时反传残差波场获得伴随状态变量，利用逆散射条件更新单炮目标泛函相对于模型参数的梯度，判断是否满足多炮循环条件，如果是，重复执行步骤3)并累加单炮目标泛函梯度，如果否，利用目标泛函梯度计算共轭梯度方向；4)基于线性Born近似正演模拟震源波场和共轭梯度波场数据，判断是否满足多炮循环条件，如果是，重复步骤4)并累加每炮共轭梯度波场数据，如果否，计算迭代步长；5)更新偏移剖面和波场残差，进入下一步迭代更新；6)判断是否满足迭代条件，如果否，终止迭代并输出迭代剖面，如果是，返回步骤2)迭代更新。
进一步，所述步骤1)中，波场残差目标泛函建立过程如下：采用高斯函数对真实速度模型进行平滑处理，得到初始平滑速度模型；波场残差为反偏移方程数值模拟的地震记录和野外观测的地震记录振幅差值，二范数意义下波场残差目标泛函如下：
其中，dsyn为反偏移的理论数据，dobs为观测的去除直达波的地震记录，ng为炮点的数目，ns为检波点的数目，xr为检波点，xs为震源，w角频率。
进一步，所述步骤2)中，采用基于CUDA计算平台，CPU/GPU协同并行加速技术，把数值模拟计算区域分割成多个16×16网格节点的子区域GPU显存，每个子区域开存储空间为28×28网格节点的共享存储器模拟整个区域波场，每个时刻把子区域数据拷贝到共享存储器中，同时拷贝相邻的子区域6个网格节点数据到该共享存储器相应的位置，把该时刻共享存储器模拟出来的波场拷贝到GPU显存，同时保存该时刻的计算区域四周边界向内部计算区域连续6个网格节点波场到计算机内存中，每一个时刻重复着上述的步骤直达最大时刻，保存最大时刻的全部计算区域波场值到GPU显存中。
进一步，所述步骤3)中，具体处理过程如下：3.1)采用保存的边界波场以及最大时刻全部波场逆时重建震源波场，把每一个时刻的CPU端的边界波场拷贝到GPU显存相应的子区域边界波场处，再用步骤3)相同的方式把子区域数值拷贝到共享存储器上，利用波动方程时间可逆性，重建最大时刻到零时刻每一个历史时刻的内部计算区域的震源波场；3.2)利用伴随状态方法，反传残差波场获得每一个时刻伴随状态变量，利用逆散射成像条件求取单炮目标泛函相对于模型参数的梯度；3.3)根据目标泛函相对于模型参数的梯度，获得第k次迭代共轭梯度dk。
进一步，所述步骤3.2)中，第k次迭代梯度(逆散射成像条件)更新表达式如下：
其中，λk为第k次迭代伴随状态变量，v0(x)为背景参考速度模型，ω是圆周频率，S(ω)为震源函数，G0(xr；x,ω)为散射点x到检波点xr的格林函数，G0(x；xs,ω)为从震源xs到散射点x的格林函数，▽为关于空间步长的梯度算子，ng,ns分别代表炮点和检波点的数目。
进一步，所述步骤3.3)中，第k次迭代共轭梯度dk的计算公式为：
其中，k为迭代次数；βk为共轭梯度系数，▽J(mk)为第k次迭代目标泛函梯度。
进一步，所述步骤4)中，迭代步长具体计算过程如下：基于声波方程线性Born近似模拟计算反射地震数据，第k次迭代单炮反射波数据Lmk如下：
Lmk＝ω2∫G0(xr；x,ω)mk(x)G0(x；xs,ω)s(ω)dx
同理计算单炮共轭梯度波场数据Ldk：
Ldk＝ω2∫G0(xr；x,ω)dk(x)G0(x；xs,ω)s(ω)dx
式中，xr为检波点位置、xs为震源位置，mk为第k次迭代反射率模型，dk第k次迭代共轭梯度方向，s(ω)为震源子波，G0(xr；x,ω)为散射点x到检波点xr的格林函数，G0(x；xs,ω)为从震源xs到散射点x的格林函数，L为线性Born正演算子如下：
L＝ω2∫G0(x；x',ω)G0(x'；xs,ω)s(ω)dx'
累加迭代步长ak表达式分子和分母，累加所有炮共轭梯度波场数据和残差数据的乘积和累加所有炮共轭梯度波场数据和共轭梯度波场数据乘积，得到第k次迭代步长ak。
进一步，所述第k次迭代步长ak为：
进一步，所述步骤5)中，根据获得的第k次迭代步长ak，利用第k次迭代共轭梯度方向dk和当前迭代反射率模型mk，计算更新偏移剖面，偏移剖面的更新表达式为：
mk+1＝mk+akdk。
进一步，所述步骤6)中，给定最大迭代次数终止条件N，当迭代次数k小于迭代终止条件N，继续迭代，当迭代次数k大于等于迭代终止条件N，输出最终的迭代剖面成像结果。
本发明由于采取以上技术方案，其具有以下优点：与常规互相关成像条件最小二乘逆时偏移相比，本发明高效的最小二乘逆时偏移梯度更新方法，考虑伴随状态法目标函数梯度会受到低频噪音干扰收敛速度和计算效率，引入线性Born逆散射成像条件应用于最小二乘逆时偏移方法中，可以有效的压制低频噪音对最小二乘逆时偏移成像质量的干扰，加快收敛速度，提高剖面的更新效率。同时CPU/GPU协同并行加速技术以及共享存储器的引入可以提高计算效率。
附图说明
图1是本发明的整体流程示意图；
图2是本发明的Salt2d反演剖面结果示意图，其中，图a是Salt2d速度模型，图b是第一次迭代结果(RTM),图c是本发明20次迭代结果(LSRTM)。
具体实施方式
下面结合附图和实施例对本发明进行详细的描述。
如图1所示，本发明提供一种最小二乘逆时偏移梯度更新方法，其包括以下步骤：
1)根据真实速度模型得到较为准确的初始速度模型，建立二范数意义下波场残差目标泛函；
建立过程如下：
采用高斯函数对真实速度模型进行平滑处理，得到较为准确的初始平滑速度模型。波场残差为反偏移方程数值模拟的地震记录和野外观测的地震记录振幅差值，二范数意义下波场残差目标泛函如下：
其中，dsyn为反偏移的理论数据，dobs为观测的去除直达波的地震记录，ng为炮点的数目，ns为检波点的数目，xr为检波点，xs为震源，w为角频率。
2)正演存储边界波场：采用基于CUDA计算平台，CPU/GPU协同并行加速方法正演模拟声波方程，存储震源边界波场和最后时刻的全部区域波场；
具体如下：
采用空间12阶精度，时间2阶精度，完全匹配层(PML)吸收边界条件数值模拟声波方程，则频域常密度声波方程为：
其中，v(x)为速度模型，ω是圆周频率，P(x,ω)为频率压力波场，S(ω)为震源函数，▽为关于空间步长的梯度算子，δ(x-xs)为狄拉克δ函数。
采用基于CUDA计算平台，CPU/GPU协同并行加速技术，把数值模拟计算区域分割成多个16×16网格节点的子区域(block)GPU显存，每个子区域开存储空间为28×28网格节点的共享存储器模拟整个区域波场，每个时刻把子区域(block)数据拷贝到共享存储器(share-memory)中，同时拷贝相邻的子区域(block)6个网格节点数据到该共享存储器相应的位置，保证该共享存储器中间内部区域16×16网格节点数据可以利用有限差分数值模拟出来，把该时刻共享存储器模拟出来的波场拷贝到GPU显存，同时保存该时刻的计算区域四周边界向内部计算区域连续6个网格节点波场到计算机内存中，每一个时刻重复着上述的步骤直达最大时刻，保存最大时刻的全部计算区域波场值到GPU显存中。
3)利用存储的边界波场和最后时刻的全部区域波场逆时重建震源波场，利用伴随状态法逆时反传残差波场获得伴随状态变量，利用逆散射条件更新单炮目标泛函相对于模型参数的梯度，判断炮数Is(Is为循环的炮数)是否满足多炮循环终止条件NS(NS为循环的总炮数)，如果满足循环条件，重复执行步骤3)并累加单炮目标泛函梯度；如果不满足循环条件，则炮循环终止；利用目标泛函梯度计算共轭梯度方向；
具体如下：
3.1)采用步骤2)保存的边界波场以及最大时刻全部波场逆时重建震源波场，把每一个时刻的CPU端的边界波场拷贝到GPU显存相应的子区域(block)边界波场处，再用步骤2)相同的方式把子区域数值拷贝到共享存储器上，利用波动方程时间可逆性，重建最大时刻到零时刻每一个历史时刻的内部计算区域的震源波场。
3.2)利用伴随状态方法，反传残差波场获得每一个时刻伴随状态变量，利用逆散射成像条件求取单炮目标泛函相对于模型参数的梯度，第k次迭代梯度更新表达式如下：
其中，λk为第k次迭代伴随状态变量，v0(x)为背景参考速度模型，ω是圆周频率，S(ω)为震源函数，G0(xr；x,ω)为散射点x到检波点xr的格林函数，G0(x；xs,ω)为从震源xs到散射点x的格林函数，▽为关于空间步长的梯度算子，ng,ns分别代表炮点和检波点的数目。
3.3)根据目标泛函相对于模型参数的梯度，利用以下公式获得第k次迭代共轭梯度dk：
其中，k为迭代次数；βk为共轭梯度系数，▽J(mk)为第k次迭代目标泛函梯度。
4)基于线性Born近似正演模拟震源波场和共轭梯度波场数据，判断炮数Is(Is为循环的炮数)是否满足多炮循环终止条件NS(NS为循环的总炮数)，如果满足循环条件，重复步骤4)并累加每炮共轭梯度波场数据；如果不满足循环条件，则终止炮循环，进行迭代步长的计算；
迭代步长具体计算过程如下：
基于声波方程线性Born近似模拟计算反射地震数据(也叫反偏移方程)，第k次迭代单炮反射波数据Lmk如下：
Lmk＝ω2∫G0(xr；x,ω)mk(x)G0(x；xs,ω)s(ω)dx
同理计算单炮共轭梯度波场数据Ldk表达式如下：
Ldk＝ω2∫G0(xr；x,ω)dk(x)G0(x；xs,ω)s(ω)dx
式中：xr为检波点位置、xs为震源位置，mk为第k次迭代反射率模型，dk第k次迭代共轭梯度方向，s(ω)为震源子波，G0(xr；x,ω)为散射点x到检波点xr的格林函数，G0(x；xs,ω)为从震源xs到散射点x的格林函数，L为线性Born正演算子(反偏移算子)如下：
L＝ω2∫G0(x；x',ω)G0(x'；xs,ω)s(ω)dx'
累加迭代步长ak表达式分子和分母，累加所有炮共轭梯度波场数据和残差数据的乘积(步长表达式分子)和累加所有炮共轭梯度波场数据和共轭梯度波场数据乘积(步长表达式分母)，第k次迭代步长ak表达式如下：
5)更新偏移剖面，进入步骤6)迭代更新。
具体如下：
根据步骤4)获得的第k次迭代步长ak，利用第k次迭代共轭梯度方向dk和当前迭代反射率模型mk，计算更新偏移剖面。偏移剖面的更新表达式为：
mk+1＝mk+akdk。
6)判断是否满足迭代终止条件，如果是，终止迭代并输出迭代剖面，如果否，返回步骤2)迭代更新；
具体如下：
给定最大迭代次数N终止条件，当迭代次数k小于迭代终止条件N，继续迭代，当迭代次数k大于等于迭代次数N，输出最终的迭代剖面成像结果。
如图2所示，为本发明对Salt2d速度模型采用最小二乘逆时偏移得到的偏移剖面，图a为Salt2d速度模型。图b为逆散射条件逆时偏移剖面，可以发现低频被很好地压制。图c为逆散射条件，最大迭代次数N＝20的最小二乘逆时偏移剖面。
综上所述，目前广泛使用的伴随状态法求取目标泛函梯度会产生低频噪音，干扰迭代收敛速度和计算效率。本发明利用线性Born逆散射成像条件迭代更新最小二乘逆时偏移梯度，可以有效的压制低频噪音对最小二乘逆时偏移成像质量的干扰，加快收敛速度，提高剖面的更新效率。同时CPU/GPU协同并行加速技术以及共享存储器的引入可以提高计算效率。
上述各实施例仅用于说明本发明，各个步骤都是可以有所变化的，在本发明技术方案的基础上，凡根据本发明原理对个别步骤进行的改进和等同变换，均不应排除在本发明的保护范围之外。

Claims (10)

1.一种最小二乘逆时偏移梯度更新方法，其特征在于包括以下步骤：
1)根据真实速度模型得到初始速度模型，建立二范数意义下波场残差目标泛函；
2)采用基于CUDA计算平台，CPU/GPU协同并行加速方法正演模拟声波方程，存储震源边界波场和最后时刻的全部区域波场；
3)利用存储的边界波场和最后时刻的全部区域波场逆时重建震源波场，利用伴随状态法逆时反传残差波场获得伴随状态变量，利用逆散射条件更新单炮目标泛函相对于模型参数的梯度，判断是否满足多炮循环条件，如果是，重复执行步骤3)并累加单炮目标泛函梯度，如果否，利用目标泛函梯度计算共轭梯度方向；
4)基于线性Born近似正演模拟震源波场和共轭梯度波场数据，判断是否满足多炮循环条件，如果是，重复步骤4)并累加每炮共轭梯度波场数据，如果否，计算迭代步长；
5)更新偏移剖面和波场残差，进入下一步迭代更新；
6)判断是否满足迭代条件，如果否，终止迭代并输出迭代剖面，如果是，返回步骤2)迭代更新。
2.如权利要求1所述的一种最小二乘逆时偏移梯度更新方法，其特征在于：所述步骤1)中，波场残差目标泛函建立过程如下：
采用高斯函数对真实速度模型进行平滑处理，得到初始平滑速度模型；波场残差为反偏移方程数值模拟的地震记录和野外观测的地震记录振幅差值，二范数意义下波场残差目标泛函如下：
<mrow> <mi>J</mi> <mrow> <mo>(</mo> <mi>m</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>1</mn> <mn>2</mn> </mfrac> <munder> <mo>&amp;Sigma;</mo> <mrow> <mi>n</mi> <mi>g</mi> <mo>,</mo> <mi>n</mi> <mi>s</mi> </mrow> </munder> <msub> <mo>&amp;Integral;</mo> <mi>w</mi> </msub> <mi>d</mi> <mi>w</mi> <mo>|</mo> <mo>|</mo> <msub> <mi>d</mi> <mrow> <mi>s</mi> <mi>y</mi> <mi>n</mi> </mrow> </msub> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>r</mi> </msub> <mo>;</mo> <msub> <mi>x</mi> <mi>s</mi> </msub> <mo>,</mo> <mi>w</mi> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>d</mi> <mrow> <mi>o</mi> <mi>b</mi> <mi>s</mi> </mrow> </msub> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>r</mi> </msub> <mo>;</mo> <msub> <mi>x</mi> <mi>s</mi> </msub> <mo>,</mo> <mi>w</mi> <mo>)</mo> </mrow> <mo>|</mo> <msup> <mo>|</mo> <mn>2</mn> </msup> </mrow>
其中，dsyn为反偏移的理论数据，dobs为观测的去除直达波的地震记录，ng为炮点的数目，ns为检波点的数目，xr为检波点，xs为震源，w角频率。
3.如权利要求1所述的一种最小二乘逆时偏移梯度更新方法，其特征在于：所述步骤2)中，采用基于CUDA计算平台，CPU/GPU协同并行加速技术，把数值模拟计算区域分割成多个16×16网格节点的子区域GPU显存，每个子区域开存储空间为28×28网格节点的共享存储器模拟整个区域波场，每个时刻把子区域数据拷贝到共享存储器中，同时拷贝相邻的子区域6个网格节点数据到该共享存储器相应的位置，把该时刻共享存储器模拟出来的波场拷贝到GPU显存，同时保存该时刻的计算区域四周边界向内部计算区域连续6个网格节点波场到计算机内存中，每一个时刻重复着上述的步骤直达最大时刻，保存最大时刻的全部计算区域波场值到GPU显存中。
4.如权利要求1所述的一种最小二乘逆时偏移梯度更新方法，其特征在于：所述步骤3)中，具体处理过程如下：
3.1)采用保存的边界波场以及最大时刻全部波场逆时重建震源波场，把每一个时刻的CPU端的边界波场拷贝到GPU显存相应的子区域边界波场处，再用步骤3)相同的方式把子区域数值拷贝到共享存储器上，利用波动方程时间可逆性，重建最大时刻到零时刻每一个历史时刻的内部计算区域的震源波场；
3.2)利用伴随状态方法，反传残差波场获得每一个时刻伴随状态变量，利用逆散射成像条件求取单炮目标泛函相对于模型参数的梯度；
3.3)根据目标泛函相对于模型参数的梯度，获得第k次迭代共轭梯度dk。
5.如权利要求4所述的一种最小二乘逆时偏移梯度更新方法，其特征在于：所述步骤3.2)中，第k次迭代梯度(逆散射成像条件)更新表达式如下：
<mfenced open = "" close = ""> <mtable> <mtr> <mtd> <mrow> <mo>&amp;dtri;</mo> <mi>J</mi> <mrow> <mo>(</mo> <msub> <mi>m</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mo>=</mo> <munder> <mo>&amp;Sigma;</mo> <mrow> <mi>n</mi> <mi>g</mi> <mo>,</mo> <mi>n</mi> <mi>s</mi> </mrow> </munder> <msub> <mo>&amp;Integral;</mo> <mi>&amp;omega;</mi> </msub> <mo>(</mo> <mo>-</mo> <mfrac> <msup> <mi>&amp;omega;</mi> <mn>2</mn> </msup> <mrow> <msubsup> <mi>v</mi> <mn>0</mn> <mn>2</mn> </msubsup> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <msub> <mi>&amp;lambda;</mi> <mi>k</mi> </msub> <msup> <mrow> <mo>&amp;lsqb;</mo> <msub> <mi>G</mi> <mn>0</mn> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>;</mo> <msub> <mi>x</mi> <mi>s</mi> </msub> <mo>,</mo> <mi>&amp;omega;</mi> <mo>)</mo> </mrow> <mi>S</mi> <mrow> <mo>(</mo> <mi>&amp;omega;</mi> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> <mo>*</mo> </msup> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mo>+</mo> <mo>&amp;dtri;</mo> <msup> <mrow> <mo>&amp;lsqb;</mo> <msub> <mi>G</mi> <mn>0</mn> </msub> <mrow> <mo>(</mo> <mi>x</mi> <mo>;</mo> <msub> <mi>x</mi> <mi>s</mi> </msub> <mo>,</mo> <mi>&amp;omega;</mi> <mo>)</mo> </mrow> <mi>S</mi> <mrow> <mo>(</mo> <mi>&amp;omega;</mi> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow> <mo>*</mo> </msup> <mo>&amp;CenterDot;</mo> <mo>&amp;dtri;</mo> <msub> <mi>&amp;lambda;</mi> <mi>k</mi> </msub> <mi>d</mi> <mi>&amp;omega;</mi> </mrow> </mtd> </mtr> </mtable> </mfenced>
其中，λk为第k次迭代伴随状态变量，v0(x)为背景参考速度模型，ω是圆周频率，S(ω)为震源函数，G0(xr；x,ω)为散射点x到检波点xr的格林函数，G0(x；xs,ω)为从震源xs到散射点x的格林函数，▽为关于空间步长的梯度算子，ng,ns分别代表炮点和检波点的数目。
6.如权利要求4所述的一种最小二乘逆时偏移梯度更新方法，其特征在于：所述步骤3.3)中，第k次迭代共轭梯度dk的计算公式为：
<mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>d</mi> <mi>k</mi> </msub> <mo>=</mo> <mo>-</mo> <mo>&amp;dtri;</mo> <mi>J</mi> <mrow> <mo>(</mo> <msub> <mi>m</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mo>+</mo> <msub> <mi>&amp;beta;</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <msub> <mi>d</mi> <mrow> <mi>k</mi> <mo>-</mo> <mn>1</mn> </mrow> </msub> <mo>;</mo> <mrow> <mo>(</mo> <mi>k</mi> <mo>&gt;</mo> <mn>0</mn> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>d</mi> <mn>0</mn> </msub> <mo>=</mo> <mo>-</mo> <mo>&amp;dtri;</mo> <mi>J</mi> <mrow> <mo>(</mo> <msub> <mi>m</mi> <mn>0</mn> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> </mfenced>
其中，k为迭代次数；βk为共轭梯度系数，▽J(mk)为第k次迭代目标泛函梯度。
7.如权利要求1所述的一种最小二乘逆时偏移梯度更新方法，其特征在于：所述步骤4)中，迭代步长具体计算过程如下：
基于声波方程线性Born近似模拟计算反射地震数据，第k次迭代单炮反射波数据Lmk如下：
Lmk＝ω2∫G0(xr；x,ω)mk(x)G0(x；xs,ω)s(ω)dx
同理计算单炮共轭梯度波场数据Ldk：
Ldk＝ω2∫G0(xr；x,ω)dk(x)G0(x；xs,ω)s(ω)dx
式中，xr为检波点位置、xs为震源位置，mk为第k次迭代反射率模型，dk第k次迭代共轭梯度方向，s(ω)为震源子波，G0(xr；x,ω)为散射点x到检波点xr的格林函数，G0(x；xs,ω)为从震源xs到散射点x的格林函数，L为线性Born正演算子如下：
L＝ω2∫G0(x；x',ω)G0(x'；xs,ω)s(ω)dx'
累加迭代步长ak表达式分子和分母，累加所有炮共轭梯度波场数据和残差数据的乘积和累加所有炮共轭梯度波场数据和共轭梯度波场数据乘积，得到第k次迭代步长ak。
8.如权利要求7所述的一种最小二乘逆时偏移梯度更新方法，其特征在于：所述第k次迭代步长ak为：
<mrow> <msub> <mi>a</mi> <mi>k</mi> </msub> <mo>=</mo> <mfrac> <mrow> <msup> <mrow> <mo>(</mo> <msub> <mi>Ld</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mrow> <mo>(</mo> <msub> <mi>d</mi> <mrow> <mi>o</mi> <mi>b</mi> <mi>s</mi> </mrow> </msub> <mo>-</mo> <msub> <mi>Lm</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> </mrow> <mrow> <msup> <mrow> <mo>(</mo> <msub> <mi>Ld</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mi>T</mi> </msup> <mrow> <mo>(</mo> <msub> <mi>Ld</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> </mrow> </mfrac> <mo>.</mo> </mrow>
9.如权利要求1所述的一种最小二乘逆时偏移梯度更新方法，其特征在于：所述步骤5)中，根据获得的第k次迭代步长ak，利用第k次迭代共轭梯度方向dk和当前迭代反射率模型mk，计算更新偏移剖面，偏移剖面的更新表达式为：
mk+1＝mk+akdk。
10.如权利要求1所述的一种最小二乘逆时偏移梯度更新方法，其特征在于：所述步骤6)中，给定最大迭代次数终止条件N，当迭代次数k小于迭代终止条件N，继续迭代，当迭代次数k大于等于迭代终止条件N，输出最终的迭代剖面成像结果。
CN201710969635.6A 2017-10-18 2017-10-18 一种最小二乘逆时偏移梯度更新方法 Active CN107783190B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710969635.6A CN107783190B (zh) 2017-10-18 2017-10-18 一种最小二乘逆时偏移梯度更新方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710969635.6A CN107783190B (zh) 2017-10-18 2017-10-18 一种最小二乘逆时偏移梯度更新方法

Publications (2)

Publication Number Publication Date
CN107783190A true CN107783190A (zh) 2018-03-09
CN107783190B CN107783190B (zh) 2019-06-11

Family

ID=61434573

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710969635.6A Active CN107783190B (zh) 2017-10-18 2017-10-18 一种最小二乘逆时偏移梯度更新方法

Country Status (1)

Country Link
CN (1) CN107783190B (zh)

Cited By (17)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108845355A (zh) * 2018-09-26 2018-11-20 中国矿业大学(北京) 地震偏移成像方法及装置
CN109239771A (zh) * 2018-08-10 2019-01-18 杭州电子科技大学 一种基于非均匀背景介质的弹性波成像方法
CN108802821B (zh) * 2018-05-28 2019-11-08 中国石油天然气股份有限公司 一种三维起伏地表地震资料偏移成像方法、装置及系统
CN110873894A (zh) * 2018-09-04 2020-03-10 中国石油化工股份有限公司 基于高斯束反偏移的炮记录获取方法及系统
CN111190224A (zh) * 2020-01-09 2020-05-22 中国石油大学(华东) 一种基于三维地震波反向照明的动态采样全波形反演系统及方法
CN111665556A (zh) * 2019-03-07 2020-09-15 中普宝信(北京)科技有限公司 地层声波传播速度模型构建方法
CN112130199A (zh) * 2020-07-31 2020-12-25 西安工程大学 一种优化的最小二乘逆时偏移成像方法
CN112630830A (zh) * 2019-10-08 2021-04-09 中国石油化工股份有限公司 一种基于高斯加权的反射波全波形反演方法及系统
CN112630824A (zh) * 2019-10-09 2021-04-09 中国石油化工股份有限公司 一种地震成像中的离散点扩散函数生成方法及系统
CN112649862A (zh) * 2019-10-12 2021-04-13 中国石油化工股份有限公司 基于地层结构信息分离的断溶体识别方法和装置
CN112698389A (zh) * 2019-10-22 2021-04-23 中国石油化工股份有限公司 一种地震资料反演成像方法及装置
CN112965102A (zh) * 2021-02-07 2021-06-15 中国石油大学(华东) 一种基于阻抗敏感核函数的最小二乘逆时偏移成像方法
CN113064203A (zh) * 2021-03-26 2021-07-02 中国海洋大学 共轭梯度归一化lsrtm方法、系统、存储介质及应用
CN113449440A (zh) * 2021-08-30 2021-09-28 中国海洋大学 一种最小二乘逆时偏移的梯度道集相关加权预处理方法
CN115166827A (zh) * 2022-07-15 2022-10-11 中山大学 基于反褶积成像条件的最小二乘偏移成像方法、设备及存储介质
CN115951401A (zh) * 2022-07-19 2023-04-11 中山大学 成像条件驱动的最小二乘逆时偏移成像方法、设备及存储介质
CN115951402A (zh) * 2022-12-29 2023-04-11 中山大学 基于向量反射率正演算子的偏移成像方法、系统及介质

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5999486A (en) * 1998-07-23 1999-12-07 Colorado School Of Mines Method for fracture detection using multicomponent seismic data
CN105005076A (zh) * 2015-06-02 2015-10-28 中国海洋石油总公司 基于最小二乘梯度更新速度模型的地震波全波形反演方法
CN105487112A (zh) * 2014-09-18 2016-04-13 中国石油化工股份有限公司 一种构建地层反射系数的方法
CN105652321A (zh) * 2015-12-30 2016-06-08 中国石油大学(华东) 一种粘声各向异性最小二乘逆时偏移成像方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5999486A (en) * 1998-07-23 1999-12-07 Colorado School Of Mines Method for fracture detection using multicomponent seismic data
CN105487112A (zh) * 2014-09-18 2016-04-13 中国石油化工股份有限公司 一种构建地层反射系数的方法
CN105005076A (zh) * 2015-06-02 2015-10-28 中国海洋石油总公司 基于最小二乘梯度更新速度模型的地震波全波形反演方法
CN105652321A (zh) * 2015-12-30 2016-06-08 中国石油大学(华东) 一种粘声各向异性最小二乘逆时偏移成像方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
王华忠 等: "最小二乘叠前深度偏移成像理论与方法", 《石油物探》 *

Cited By (25)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108802821B (zh) * 2018-05-28 2019-11-08 中国石油天然气股份有限公司 一种三维起伏地表地震资料偏移成像方法、装置及系统
CN109239771A (zh) * 2018-08-10 2019-01-18 杭州电子科技大学 一种基于非均匀背景介质的弹性波成像方法
CN109239771B (zh) * 2018-08-10 2020-01-31 杭州电子科技大学 一种基于非均匀背景介质的弹性波成像方法
CN110873894A (zh) * 2018-09-04 2020-03-10 中国石油化工股份有限公司 基于高斯束反偏移的炮记录获取方法及系统
CN108845355A (zh) * 2018-09-26 2018-11-20 中国矿业大学(北京) 地震偏移成像方法及装置
CN111665556A (zh) * 2019-03-07 2020-09-15 中普宝信(北京)科技有限公司 地层声波传播速度模型构建方法
CN111665556B (zh) * 2019-03-07 2023-05-02 中普宝信(北京)科技有限公司 地层声波传播速度模型构建方法
CN112630830A (zh) * 2019-10-08 2021-04-09 中国石油化工股份有限公司 一种基于高斯加权的反射波全波形反演方法及系统
CN112630830B (zh) * 2019-10-08 2024-04-09 中国石油化工股份有限公司 一种基于高斯加权的反射波全波形反演方法及系统
CN112630824A (zh) * 2019-10-09 2021-04-09 中国石油化工股份有限公司 一种地震成像中的离散点扩散函数生成方法及系统
CN112630824B (zh) * 2019-10-09 2024-03-22 中国石油化工股份有限公司 一种地震成像中的离散点扩散函数生成方法及系统
CN112649862A (zh) * 2019-10-12 2021-04-13 中国石油化工股份有限公司 基于地层结构信息分离的断溶体识别方法和装置
CN112698389A (zh) * 2019-10-22 2021-04-23 中国石油化工股份有限公司 一种地震资料反演成像方法及装置
CN112698389B (zh) * 2019-10-22 2024-02-20 中国石油化工股份有限公司 一种地震资料反演成像方法及装置
CN111190224A (zh) * 2020-01-09 2020-05-22 中国石油大学(华东) 一种基于三维地震波反向照明的动态采样全波形反演系统及方法
CN112130199A (zh) * 2020-07-31 2020-12-25 西安工程大学 一种优化的最小二乘逆时偏移成像方法
CN112965102A (zh) * 2021-02-07 2021-06-15 中国石油大学(华东) 一种基于阻抗敏感核函数的最小二乘逆时偏移成像方法
CN113064203A (zh) * 2021-03-26 2021-07-02 中国海洋大学 共轭梯度归一化lsrtm方法、系统、存储介质及应用
CN113449440A (zh) * 2021-08-30 2021-09-28 中国海洋大学 一种最小二乘逆时偏移的梯度道集相关加权预处理方法
CN115166827B (zh) * 2022-07-15 2023-04-28 中山大学 基于反褶积成像条件的最小二乘偏移成像方法、设备及存储介质
CN115166827A (zh) * 2022-07-15 2022-10-11 中山大学 基于反褶积成像条件的最小二乘偏移成像方法、设备及存储介质
CN115951401A (zh) * 2022-07-19 2023-04-11 中山大学 成像条件驱动的最小二乘逆时偏移成像方法、设备及存储介质
CN115951401B (zh) * 2022-07-19 2023-09-15 中山大学 成像条件驱动的最小二乘逆时偏移成像方法、设备及存储介质
CN115951402A (zh) * 2022-12-29 2023-04-11 中山大学 基于向量反射率正演算子的偏移成像方法、系统及介质
CN115951402B (zh) * 2022-12-29 2023-07-18 中山大学 基于向量反射率正演算子的偏移成像方法、系统及介质

Also Published As

Publication number Publication date
CN107783190B (zh) 2019-06-11

Similar Documents

Publication Publication Date Title
CN107783190B (zh) 一种最小二乘逆时偏移梯度更新方法
CN106526674B (zh) 一种三维全波形反演能量加权梯度预处理方法
Higuera et al. Three-dimensional numerical wave generation with moving boundaries
CN105319581B (zh) 一种高效的时间域全波形反演方法
CN103135132B (zh) Cpu/gpu协同并行计算的混合域全波形反演方法
CN106908835B (zh) 带限格林函数滤波多尺度全波形反演方法
CN107894618B (zh) 一种基于模型平滑算法的全波形反演梯度预处理方法
CN110058302A (zh) 一种基于预条件共轭梯度加速算法的全波形反演方法
CN105137486A (zh) 各向异性介质中弹性波逆时偏移成像方法及其装置
CN104977607B (zh) 利用变步长网格声波波场模拟的时间域全波形反演方法
CN105388520B (zh) 一种地震资料叠前逆时偏移成像方法
CN110187382B (zh) 一种回折波和反射波波动方程旅行时反演方法
CN110058307A (zh) 一种基于快速拟牛顿法的全波形反演方法
CN106932820A (zh) 基于时域伪谱方法的声波方程逆时偏移成像方法
CN109946742A (zh) 一种TTI介质中纯qP波地震数据模拟方法
CN108680968B (zh) 复杂构造区地震勘探数据采集观测系统评价方法及装置
CN105182414B (zh) 一种基于波动方程正演去除直达波的方法
CN106353798A (zh) 多分量联合高斯束叠前逆时偏移成像方法
CN109239776B (zh) 一种地震波传播正演模拟方法和装置
CN106597535A (zh) 一种提高弹性波逆时偏移计算率和空间分辨率的方法
CN108445532B (zh) 一种深度域反偏移方法及装置
Jin et al. 2D multiscale non‐linear velocity inversion
Wang et al. Acoustic reverse time migration and perfectly matched layer in boundary-conforming grids by elliptic method
Wang et al. Improved hybrid iterative optimization method for seismic full waveform inversion
Bartana* et al. GPU implementation of minimal dispersion recursive operators for reverse time migration

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