CN114690242A - 一种低噪音的最小二乘逆时偏移方法 - Google Patents

一种低噪音的最小二乘逆时偏移方法 Download PDF

Info

Publication number
CN114690242A
CN114690242A CN202210201779.8A CN202210201779A CN114690242A CN 114690242 A CN114690242 A CN 114690242A CN 202210201779 A CN202210201779 A CN 202210201779A CN 114690242 A CN114690242 A CN 114690242A
Authority
CN
China
Prior art keywords
function
gradient
reverse time
reflection coefficient
time migration
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.)
Pending
Application number
CN202210201779.8A
Other languages
English (en)
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 CN202210201779.8A priority Critical patent/CN114690242A/zh
Publication of CN114690242A publication Critical patent/CN114690242A/zh
Pending legal-status Critical Current

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
    • 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
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/15—Correlation function computation including computation of convolution operations

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Remote Sensing (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Computational Mathematics (AREA)
  • Pure & Applied Mathematics (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geophysics (AREA)
  • Acoustics & Sound (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Geology (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Algebra (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Computing Systems (AREA)
  • Image Processing (AREA)

Abstract

本发明公开了一种低噪音的最小二乘逆时偏移方法，包括：建立最小二乘逆时偏移的函数；通过对最小二乘逆时偏移的函数添加稀疏约束项来压制噪音，经过柯西约束项稀疏约束后的目标函数；确定最小二乘逆时偏移的函数经过柯西约束项稀疏约束后的目标函数的参数，并对其梯度进行优化和采用预条件非线性共轭梯度法对最小二乘逆时偏移的函数进行最小化求解。本发明公开的一种低噪音的最小二乘逆时偏移方法，给出了一种新的确定柯西项中正则化参数和超参数的方法，可得到更好的成像结果和加快反演的收敛速度，给含有柯西约束项的目标函数提供了一种确定步长的方法。

Description

一种低噪音的最小二乘逆时偏移方法
技术领域
本发明涉及勘探地球物理和反演综合技术领域，具体涉及一种低噪音的最小二乘逆时偏移方法。
背景技术
偏移成像技术是一项地震数据处理中的关键技术。它通过将地表接收到的地震反射事件恢复到地下的真实位置从而对地下构造进行成像。高分辨率的偏移成像结果可以帮助地质解释人员有效的寻找到地下油气储层。因此，为了得到高分辨率的地下成像结果，偏移技术先是从基于射线追踪的克希霍夫偏移发展到高斯束偏移，然后从基于单程波波动方程的单程波偏移到基于双程波波动方程的逆时偏移。目前，逆时偏移由于其具有对复杂构造成像的能力而成为最受欢迎的偏移技术之一。但是其只采用了正演算子的伴随算子而不是逆算子，因此它的成像结果依旧具有振幅不准确，照明不均衡等问题。为了进一步提高成像质量，最小二乘逆时偏移通过将逆时偏移和最小二乘反演相结合。它既有逆时偏移所拥有的能够对复杂构造进行成像的优势，又兼具最小二乘反演振幅更均衡、分辨率更高的优势。但是最小二乘反演和其它地球物理反演方法一样也是一个病态问题，这是由于地震数据的频带有限，以及检波器的偏移距和观测方位有限所导致。该现象会使得所成像的反射轴周围存在很多的旁瓣，以及会拟合零空间数据从而产生高频噪音。除此之外，由于最小二乘逆时偏移的梯度是正传波场和反传波场进行互相关。因此会在强反射界面上方产生大量的低频噪音以及部分错误成像，这些低频噪音会将真正的成像界面所掩盖。其中低频噪音在最小二乘逆时偏移的迭代过程中会影响迭代的收敛速度和成像结果，错误成像和旁瓣这一类高频噪音则会随着迭代不断的被加强从而影响最终的成像质量。
发明内容
本发明的目的在于提供一种低噪音的最小二乘逆时偏移方法，用以解决上述噪音在最小二乘逆时偏移的迭代过程中会影响迭代的收敛速度和成像结果的问题。
本发明提供一种低噪音的最小二乘逆时偏移方法，包括：
建立最小二乘逆时偏移的函数F(m)；
通过对最小二乘逆时偏移的函数F(m)添加稀疏约束项R(m)来压制噪音，最小二乘逆时偏移的函数F(m)经过柯西约束项稀疏约束后的目标函数J(m)的表达式为：
J(m)＝F(m)+βR(m) (式2)
式中，J(m)为最小二乘逆时偏移的函数F(m)经过柯西约束项稀疏约束后的目标函数，β表示正则化参数，R(m)表示柯西约束项；
确定最小二乘逆时偏移的函数F(m)经过柯西约束项稀疏约束后的目标函数J(m)的参数，并对其梯度进行优化和采用预条件非线性共轭梯度法对最小二乘逆时偏移的函数F(m)进行最小化求解。
具体地，所述柯西约束项R(m)的表达式为：
Figure BDA0003527704760000021
式中，N表示反射系数模型m元素的总数量，λk表示控制反射系数模型m稀疏性的超参数，k表示反射系数模型m元素的序号，mk表示反射系数模型m中的第k个元素。
进一步地，所述对最小二乘逆时偏移的函数F(m)进行最小化求解具体包括以下步骤：
步骤A：确定最小二乘逆时偏移的函数的正则化参数β和控制反射系数模型m稀疏性的超参数λk来求取目标函数J(m)的梯度，并且对最小二乘逆时偏移的函数的梯度进行优化；
步骤B：采用预条件非线性共轭梯度法对反射系数模型m进行迭代更新。
进一步地，所述步骤A包括以下步骤：
步骤A1：通过照明补偿系数确定柯西约束项R(m)中控制反射系数模型m稀疏性的超参数λk，其中，照明补偿系数的表达式为：
Figure BDA0003527704760000022
式中，M(x)为照明补偿系数，S表示震源的正传波场，t表示波场传播的当前时刻，x表示反射点位置，xs表示震源的位置，Tmax表示波场传播的最大时间；
步骤A2：通过乘法损失函数Jn(m)来确定正则化参数β，其中，乘法损失函数Jn(m)的表达式为：
Jn(m)＝(F(m)-ε)Cn(m) (式5)
式中，Jn(m)表示乘法损失函数，F(m)表示最小二乘逆时偏移的函数，ε表示松弛因子，Cn(m)表示乘法损失函数的约束项；
步骤A3：预设乘法损失函数Jn(m)中的松弛因子ε；
步骤A4：计算乘法损失函数Jn(m)沿着反射系数模型m的梯度
Figure BDA0003527704760000031
其计算表达式为
Figure BDA0003527704760000032
式中，ε为松弛因子，F(m)为最小二乘逆时偏移的函数，R(mn)为每次迭代后反射系数模型计算的柯西约束项进行常数归一化的函数值；
Figure BDA0003527704760000033
为最小二乘逆时偏移的函数F(m)对反射系数模型m的梯度；
Figure BDA0003527704760000034
为柯西约束项R(m)对反射系数模型m的梯度；
步骤A5：修正乘法损失函数Jn(m)沿着反射系数模型m的梯度
Figure BDA0003527704760000035
进一步地，所述步骤A5包括以下步骤：
通过新的近似更新梯度来消除第一项梯度
Figure BDA0003527704760000036
中的低频噪音以及部分错误成像，所述新的近似的梯度可以表示为：
Figure BDA0003527704760000037
式中，S表示震源的正传波场；t表示波场传播的当前时刻；x表示反射点位置；xr表示检波器的位置；xs表示震源的位置；G(x|xr,t)＊Δd(xr,t)表示残差记录Δd(xr,t)的反传波场；Δd(xr,t)表示残差记录；
Figure BDA0003527704760000038
表示先对S求取时间导数，然后沿空间Z方向做希尔伯特变换；Hz(G(x|xr,t)＊Δd(xr,t))表示对G(x|xr,t)＊Δd(xr,t)沿空间Z方向做希尔伯特变换；Ht(Δd(xr,t))表示对Δd(xr,t)沿时间t方向做希尔伯特变换；G(x|xr,t)＊Ht(Δd(xr,t))表示Ht(Δd(xr,t))的反传波场；Hz(G(x|xr,t)＊Ht(Δd(xr,t)))表示对G(x|xr,t)＊Ht(Δd(xr,t))沿空间Z方向做希尔伯特变换。
进一步地，所述步骤B包括以下步骤：
步骤B1：输入初始反射系数模型m0；
步骤B2：计算乘法损失函数Jn(m)沿着反射系数模型m的梯度
Figure BDA0003527704760000041
并据此通过非线性共轭梯度法求取初始乘法损失函数Jn(m0)沿着初始反射系数模型m0的梯度
Figure BDA0003527704760000042
的共轭梯度g0；
步骤B3：计算非线性共轭梯度法的更新步长α；
步骤B4：通过更新步长α和共轭梯度g0对初始反射系数模型m0进行更新，其更新表达式为：
m1＝m0+αg0 (式9)
式中，m1表示更新后的反射系数模型；m0为初始反射系数模型；α为非线性共轭梯度法的更新步长，g0为初始乘法损失函数Jn(m0)沿着初始反射系数模型m0的梯度
Figure BDA0003527704760000043
的共轭梯度；
步骤B5：重复步骤B1至步骤B3，直到达到迭代次数或者满足终止条件。
具体地，所述步骤B3包括以下步骤：
步骤B31：先给定一个较小的步长α0来对初始反射系数模型m0进行扰动得到扰动后的模型m0+α0g0，然后对扰动后的模型m0+α0g0计算获取目标模型值Jn(m0+α0g0)，根据目标模型值Jn(m0+α0g0)和未扰动的目标函数值Jn(m0)以及零点确定一条抛物线；
步骤B32：通过二分法对所述抛物线求取目标模型值Jn(m0+αg0)最小时的步长，即为非线性共轭梯度法的更新步长α。
本发明还公开了一种低噪音的最小二乘逆时偏移装置，包括
第一处理单元，用于建立最小二乘逆时偏移的函数F(m)；
第二处理单元，用于通过对最小二乘逆时偏移的函数F(m)添加稀疏约束项R(m)来压制噪音，最小二乘逆时偏移的函数F(m)经过柯西约束项稀疏约束后的目标函数J(m)的表达式为：
J(m)＝F(m)+βR(m) (式2)
式中，J(m)为最小二乘逆时偏移的函数F(m)经过柯西约束项稀疏约束后的目标函数，β表示正则化参数，R(m)表示柯西约束项；
第三处理单元，用于确定最小二乘逆时偏移的函数F(m)经过柯西约束项稀疏约束后的目标函数J(m)的参数，并对其梯度进行优化和采用预条件非线性共轭梯度法对最小二乘逆时偏移的函数F(m)进行最小化求解。
本发明还公开了一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现上述的方法的步骤。
本发明还公开了一种计算机设备，包括存储器、处理器以及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现上述的方法的步骤。
与现有技术相比，本发明的有益效果是：
本发明公开的一种低噪音的最小二乘逆时偏移方法，通过将柯西约束引入到最小二乘逆时偏移来压制由于地震数据的频带有限，以及检波器的偏移距和观测方位有限所导致旁瓣和噪音，并且给出了一种新的确定柯西项中正则化参数和超参数的方法，可以做到更有效和更精确的确定这两个参数；除此之外，还对最小二乘逆时偏移的梯度进行了修正来压制低频噪音和错误成像，从而得到更好的成像结果和加快反演的收敛速度；最后给含有柯西约束项的目标函数提供了一种新的快速确定步长的方法。
附图说明
图1是本发明实施例1提供的真实速度模型与偏移速度模型的对比图。其中，图1(a)为真实速度模型，用来求解双程波波动方程来得到观测地震记录如图2所示，图1(b)为偏移速度模型，用来作为最小二乘逆时偏移的偏移速度模型；
图2是通过对本发明实施例1提供的真实速度模型求解双程波波动方程所产生的炮点位于0.274km和13.99km处的观测地震记录；
图3是本发明实施例1提供的对反射系数模型m进行柯西稀疏性约束的超参数λk的深度与距离的关系图；
图4是本发明实施例1提供的传统最小二乘逆时偏移和本发明方法的偏移结果对比图，其中，图4(a)为传统最小二乘逆时偏移的成像结果；图4(b)是本发明所提供的低噪音最小二乘逆时偏移的成像结果。
具体实施方式
以下实施例用于说明本发明，但不用来限制本发明的范围。
实施例1
实施例1提供一种低噪音的最小二乘逆时偏移方法，包括以下步骤：
首先，建立最小二乘逆时偏移的函数F(m)，其表达式为：
Figure BDA0003527704760000061
式中，F(m)为最小二乘逆时偏移的函数；m表示反射系数模型，是一个矩阵；d和dobs分别表示预测数据和观测数据。
其中，观测数据dobs通过对真实速度模型求解双程波波动方程正演得到，预测数据d通过对偏移速度模型和反射系数模型m采用柯西霍夫近似正演得到，如图1和2所示。
然后，通过对最小二乘逆时偏移的函数F(m)添加稀疏约束项R(m)来压制由于数据的频带和检波点的偏移距有限所引起的旁瓣和拟合零空间数据所产生的噪音，最小二乘逆时偏移的函数F(m)经过柯西约束项稀疏约束后形成目标函数J(m)，其表达式为：
J(m)＝F(m)+βR(m) (式2)
式中，J(m)为最小二乘逆时偏移的函数F(m)经过柯西约束项稀疏约束后的目标函数，F(m)表示最小二乘逆时偏移的函数，β表示正则化参数，R(m)表示柯西约束项，其表达式为：
Figure BDA0003527704760000062
式中，N表示反射系数模型m元素的总数量，λk表示控制反射系数模型m稀疏性的超参数，k表示反射系数模型m元素的序号，mk表示反射系数模型m中的第k个元素。
最后，确定最小二乘逆时偏移的函数F(m)经过柯西约束项稀疏约束后的目标函数J(m)的参数，并对其梯度进行优化和采用预条件非线性共轭梯度法对最小二乘逆时偏移的函数F(m)进行最小化求解，具体包括以下步骤：
步骤A：确定最小二乘逆时偏移的函数的正则化参数β和控制反射系数模型m稀疏性的超参数λk来求取目标函数J(m)的梯度，并且对最小二乘逆时偏移的函数F(m)的梯度进行优化，具体包括以下步骤：
步骤A1：通过照明补偿系数确定柯西约束项R(m)中控制反射系数模型m稀疏性的超参数λk；
通过几次无约束的最小二乘逆时偏移得到一个偏移后的反射系数模型m，然后设置控制反射系数模型m稀疏性的超参数λk的值介于模型中噪音值和真实反射值之间。但是由于反射系数模型m中不同空间位置的噪音值和真实反射值并不一致，因此对整个反射系数的模型空间只设置一个超参数λk并不合适。因此需要借助照明补偿系数来进一步微调反射系数m不同空间位置所对应的超参数λk的取值，照明补偿越强的位置，超参数λk的值越小，反之则越弱。照明补偿系数也被称之为海森矩阵的对角阵，其表达式为
Figure BDA0003527704760000071
式中，M(x)为照明补偿系数，S表示震源的正传波场，t表示波场传播的当前时刻，x表示反射点位置，xs表示震源的位置，Tmax表示波场传播的最大时间。
步骤A2：通过乘法损失函数Jn(m)来自动化确定正则化参数β，其中，乘法损失函数Jn(m)的表达式为：
Jn(m)＝(F(m)-ε)Cn(m) (式5)
式中，Jn(m)表示乘法损失函数，F(m)表示最小二乘逆时偏移的函数，ε表示松弛因子，Cn(m)表示乘法损失函数的约束项；
其中，乘法损失函数的约束项Cn(m)的表达式为：
Figure BDA0003527704760000072
式中，mn表示第n次迭代的反射系数模型；
R(m)表示柯西约束项；
R(mn)表示每次迭代后第n次反射系数模型mn计算的柯西约束项进行常数归一化的函数值。
步骤A3：预设乘法损失函数Jn(m)中的松弛因子ε；
由于乘法损失函数Jn(m)包含有松弛因子ε，松弛因子ε可以作为常数，因此需要预设松弛因子ε，一般可以设置为无约束迭代几次后目标函数值的十分之一；
步骤A4：计算乘法损失函数Jn(m)沿着反射系数模型m的梯度
Figure BDA0003527704760000073
其计算表达式为：
Figure BDA0003527704760000074
式中，ε为松弛因子，F(m)为最小二乘逆时偏移的函数，R(mn)为每次迭代后反射系数模型计算的柯西约束项进行常数归一化的函数值；
Figure BDA0003527704760000075
为最小二乘逆时偏移的函数F(m)对反射系数模型m的梯度；
Figure BDA0003527704760000081
为柯西约束项R(m)对反射系数模型m的梯度；
其中，最小二乘逆时偏移的函数F(m)对反射系数模型m的梯度
Figure BDA0003527704760000082
和柯西约束项R(m)对反射系数模型m的梯度
Figure BDA0003527704760000083
共同组成了乘法损失函数Jn(m)沿着反射系数模型m的梯度
Figure BDA0003527704760000084
的两项变量。
步骤A5：修正乘法损失函数Jn(m)沿着反射系数模型m的梯度
Figure BDA0003527704760000085
由于松弛因子ε、每次迭代后反射系数模型计算的柯西约束项进行常数归一化的函数值R(mn)和最小二乘逆时偏移的函数F(m)可以作为常数，因此
Figure BDA0003527704760000086
可以认为是常数，故乘法损失函数Jn(m)沿着反射系数模型m的梯度
Figure BDA0003527704760000087
的变量包含有两项：
第一项是最小二乘逆时偏移的函数F(m)对反射系数模型m的梯度
Figure BDA0003527704760000088
第二项是柯西约束项R(m)对反射系数模型m的梯度
Figure BDA0003527704760000089
其中最小二乘逆时偏移的函数F(m)对反射系数模型m的梯度通常会在强反射界面上方产生大量的低频噪音以及部分错误成像。
为此，通过新的近似更新梯度来消除第一项梯度
Figure BDA00035277047600000810
中的低频噪音以及部分错误成像，这种新的近似的梯度可以表示为：
Figure BDA00035277047600000811
式中，S表示震源的正传波场；t表示波场传播的当前时刻；x表示反射点位置；xr表示检波器的位置；xs表示震源的位置；G(x|xr,t)＊Δd(xr,t)表示残差记录Δd(xr,t)的反传波场；Δd(xr,t)表示残差记录；
Figure BDA0003527704760000091
表示先对S求取时间导数，然后沿空间Z方向做希尔伯特变换；Hz(G(x|xr,t)＊Δd(xr,t))表示对G(x|xr,t)＊Δd(xr,t)沿空间Z方向做希尔伯特变换；Ht(Δd(xr,t))表示对Δd(xr,t)沿时间t方向做希尔伯特变换；G(x|xr,t)＊Ht(Δd(xr,t))表示Ht(Δd(xr,t))的反传波场；Hz(G(x|xr,t)＊Ht(Δd(xr,t)))表示对G(x|xr,t)＊Ht(Δd(xr,t))沿空间Z方向做希尔伯特变换。
步骤B：采用预条件非线性共轭梯度法对反射系数模型m进行迭代更新。
步骤B1：输入初始反射系数模型m0；
步骤B2：计算乘法损失函数Jn(m)沿着反射系数模型m的梯度
Figure BDA0003527704760000092
并据此通过非线性共轭梯度法求取初始乘法损失函数Jn(m0)沿着初始反射系数模型m0的梯度
Figure BDA0003527704760000093
的共轭梯度g0；
步骤B3：通过改进的步长搜索方法来计算非线性共轭梯度法的更新步长α，具体步骤如下：
步骤B31：先给定一个较小的步长α0来对初始反射系数模型m0进行扰动得到扰动后的模型m0+α0g0，然后对扰动后的模型m0+α0g0计算获取目标模型值Jn(m0+α0g0)，根据目标模型值Jn(m0+α0g0)和未扰动的目标函数值Jn(m0)以及零点确定一条抛物线；
步骤B32：通过二分法对所述抛物线求取目标模型值Jn(m0+αg0)最小时的步长，即为非线性共轭梯度法的更新步长α；
步骤B4：通过更新步长α和共轭梯度g0对初始反射系数模型m0进行更新，其更新表达式为：
m1＝m0+αg0 (式9)
式中，m1表示更新后的反射系数模型；m0为初始反射系数模型；α为非线性共轭梯度法的更新步长，g0为初始乘法损失函数Jn(m0)沿着初始反射系数模型m0的梯度
Figure BDA0003527704760000094
的共轭梯度。
步骤B5：重复步骤B2至步骤B3，直到达到最大迭代次数或者满足终止条件。
下面通过具体实施例详细说明去除最小二乘逆时偏移中部分成像噪音的方法。
由图1可知，真实速度模型包含有巨大的盐丘，这对包括逆时偏移在内的偏移算法都提出了巨大的挑战，主要原因之一是盐丘会阻挡了地震波的穿透，导致盐丘下部区域照明较弱，使得该区域通常难以成像。其中，图1(a)中的观测数据集一共包含有160炮，采用拖缆的方式进行激发，炮点间距为137.16m，检波器的间隔为46.72m，最小和最大偏移距分别为0m和7978m。
图2所示的是炮点位于0.274km和13.99km处的单炮记录。
参考图3，由于盐丘下部区域照明较弱，因此难以成像，为了避免柯西约束将该区域的有效信号进一步压制，需要设置更大的λk值来减少对该区域的稀疏约束；
由图4(a)可知，传统最小二乘逆时偏移的方法得到的成像结果存在大量低频噪音(图4a中黑色实心箭头所示)和包含错误成像在内的旁瓣、划弧等高频噪音(图4a中黑色空心箭头所示)；结合图4(b)可知，对比本发明所提供的低噪音最小二乘逆时偏移方法得到的成像结果，可以看出，本发明最终提供的最终偏移结果可以有效的压制图4a中所示的噪音。
虽然，上文中已经用一般性说明及具体实施例对本发明作了详尽的描述，但在本发明基础上，可以对之作一些修改或改进，这对本领域技术人员而言是显而易见的。因此，在不偏离本发明精神的基础上所做的这些修改或改进，均属于本发明要求保护的范围。

Claims (10)

1.一种低噪音的最小二乘逆时偏移方法，其特征在于，包括：
建立最小二乘逆时偏移的函数F(m)；
通过对最小二乘逆时偏移的函数F(m)添加稀疏约束项R(m)来压制噪音，最小二乘逆时偏移的函数F(m)经过柯西约束项稀疏约束后的目标函数J(m)的表达式为：
J(m)＝F(m)+βR(m) (式2)
式中，J(m)为最小二乘逆时偏移的函数F(m)经过柯西约束项稀疏约束后的目标函数，β表示正则化参数，R(m)表示柯西约束项；
确定最小二乘逆时偏移的函数F(m)经过柯西约束项稀疏约束后的目标函数J(m)的参数，并对其梯度进行优化和采用预条件非线性共轭梯度法对最小二乘逆时偏移的函数F(m)进行最小化求解。
2.根据权利要求1所述的方法，其特征在于，
所述柯西约束项R(m)的表达式为：
Figure FDA0003527704750000011
式中，N表示反射系数模型m元素的总数量，λk表示控制反射系数模型m稀疏性的超参数，k表示反射系数模型m元素的序号，mk表示反射系数模型m中的第k个元素。
3.根据权利要求1所述的方法，其特征在于，
所述对最小二乘逆时偏移的函数F(m)进行最小化求解具体包括以下步骤：
步骤A：确定最小二乘逆时偏移的函数F(m)的正则化参数β和控制反射系数模型m稀疏性的超参数λk来求取目标函数J(m)的梯度，并且对该目标函数J(m)的梯度进行优化；
步骤B：采用预条件非线性共轭梯度法对反射系数模型m进行迭代更新。
4.根据权利要求3所述的方法，其特征在于，所述步骤A包括以下步骤：
步骤A1：通过照明补偿系数确定柯西约束项R(m)中控制反射系数模型m稀疏性的超参数λk，其中，照明补偿系数的表达式为：
Figure FDA0003527704750000012
式中，M(x)为照明补偿系数，S表示震源的正传波场，t表示波场传播的当前时刻，x表示反射点位置，xs表示震源的位置，Tmax表示波场传播的最大时间；
步骤A2：通过乘法损失函数Jn(m)来确定正则化参数β，其中，乘法损失函数Jn(m)的表达式为：
Jn(m)＝(F(m)-ε)Cn(m) (式5)
式中，Jn(m)表示乘法损失函数，F(m)表示最小二乘逆时偏移的函数，ε表示松弛因子，Cn(m)表示乘法损失函数的约束项；
步骤A3：预设乘法损失函数Jn(m)中的松弛因子ε；
步骤A4：计算乘法损失函数Jn(m)沿着反射系数模型m的梯度
Figure FDA0003527704750000021
其计算表达式为：
Figure FDA0003527704750000022
式中，ε为松弛因子，F(m)为最小二乘逆时偏移的函数，R(mn)为每次迭代后反射系数模型计算的柯西约束项进行常数归一化的函数值；
Figure FDA0003527704750000023
为最小二乘逆时偏移的函数F(m)对反射系数模型m的梯度；
Figure FDA0003527704750000024
为柯西约束项R(m)对反射系数模型m的梯度；
步骤A5：修正乘法损失函数Jn(m)沿着反射系数模型m的梯度
Figure FDA0003527704750000025
5.根据权利要求4所述的方法，其特征在于，所述步骤A5包括以下步骤：
通过新的近似更新梯度来消除第一项梯度
Figure FDA0003527704750000026
中的低频噪音以及部分错误成像，所述新的近似的梯度可以表示为：
Figure FDA0003527704750000027
式中，S表示震源的正传波场；t表示波场传播的当前时刻；x表示反射点位置；xr表示检波器的位置；xs表示震源的位置；G(x|xr,t)*Δd(xr,t)表示残差记录Δd(xr,t)的反传波场；Δd(xr,t)表示残差记录；
Figure FDA0003527704750000031
表示先对S求取时间导数，然后沿空间Z方向做希尔伯特变换；Hz(G(x|xr,t)*Δd(xr,t))表示对G(x|xr,t)*Δd(xr,t)沿空间Z方向做希尔伯特变换；Ht(Δd(xr,t))表示对Δd(xr,t)沿时间t方向做希尔伯特变换；G(x|xr,t)*Ht(Δd(xr,t))表示Ht(Δd(xr,t))的反传波场；Hz(G(x|xr,t)*Ht(Δd(xr,t)))表示对G(x|xr,t)*Ht(Δd(xr,t))沿空间Z方向做希尔伯特变换。
6.根据权利要求2所述的方法，其特征在于，所述步骤B包括以下步骤：
步骤B1：输入初始反射系数模型m0；
步骤B2：计算乘法损失函数Jn(m)沿着反射系数模型m的梯度
Figure FDA0003527704750000032
并据此通过非线性共轭梯度法求取初始乘法损失函数Jn(m0)沿着初始反射系数模型m0的梯度
Figure FDA0003527704750000033
的共轭梯度g0；
步骤B3：计算非线性共轭梯度法的更新步长α；
步骤B4：通过更新步长α和共轭梯度g0对初始反射系数模型m0进行更新，其更新表达式为：
m1＝m0+αg0 (式9)
式中，m1表示更新后的反射系数模型；m0为初始反射系数模型；α为非线性共轭梯度法的更新步长，g0为初始乘法损失函数Jn(m0)沿着初始反射系数模型m0的梯度
Figure FDA0003527704750000034
的共轭梯度；
步骤B5：重复步骤B1至步骤B3，直到达到迭代次数或者满足终止条件。
7.根据权利要求6所述的方法，其特征在于，所述步骤B3包括以下步骤：
步骤B31：先给定一个较小的步长α0来对初始反射系数模型m0进行扰动得到扰动后的模型m0+α0g0，然后对扰动后的模型m0+α0g0计算获取目标模型值Jn(m0+α0g0)，根据目标模型值Jn(m0+α0g0)和未扰动的目标函数值Jn(m0)以及零点确定一条抛物线；
步骤B32：通过二分法对所述抛物线求取目标模型值Jn(m0+αg0)最小时的步长，即为非线性共轭梯度法的更新步长α。
8.一种低噪音的最小二乘逆时偏移装置，其特征在于，包括
第一处理单元，用于建立最小二乘逆时偏移的函数F(m)；
第二处理单元，用于通过对最小二乘逆时偏移的函数F(m)添加稀疏约束项R(m)来压制噪音，最小二乘逆时偏移的函数F(m)经过柯西约束项稀疏约束后的目标函数J(m)的表达式为：
J(m)＝F(m)+βR(m) (式2)
式中，J(m)为最小二乘逆时偏移的函数F(m)经过柯西约束项稀疏约束后的目标函数，β表示正则化参数，R(m)表示柯西约束项；
第三处理单元，用于确定最小二乘逆时偏移的函数F(m)经过柯西约束项稀疏约束后的目标函数J(m)的参数，并对其梯度进行优化和采用预条件非线性共轭梯度法对最小二乘逆时偏移的函数F(m)进行最小化求解。
9.一种计算机可读存储介质，其上存储有计算机程序，其特征在于，所述计算机程序被处理器执行时实现权利要求1-7中任一项所述的方法的步骤。
10.一种计算机设备，包括存储器、处理器以及存储在存储器上并可在处理器上运行的计算机程序，其特征在于，所述处理器执行所述计算机程序时实现权利要求1-7任一所述的方法的步骤。
CN202210201779.8A 2022-03-02 2022-03-02 一种低噪音的最小二乘逆时偏移方法 Pending CN114690242A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210201779.8A CN114690242A (zh) 2022-03-02 2022-03-02 一种低噪音的最小二乘逆时偏移方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210201779.8A CN114690242A (zh) 2022-03-02 2022-03-02 一种低噪音的最小二乘逆时偏移方法

Publications (1)

Publication Number Publication Date
CN114690242A true CN114690242A (zh) 2022-07-01

Family

ID=82138103

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210201779.8A Pending CN114690242A (zh) 2022-03-02 2022-03-02 一种低噪音的最小二乘逆时偏移方法

Country Status (1)

Country Link
CN (1) CN114690242A (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115951401A (zh) * 2022-07-19 2023-04-11 中山大学 成像条件驱动的最小二乘逆时偏移成像方法、设备及存储介质

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115951401A (zh) * 2022-07-19 2023-04-11 中山大学 成像条件驱动的最小二乘逆时偏移成像方法、设备及存储介质
CN115951401B (zh) * 2022-07-19 2023-09-15 中山大学 成像条件驱动的最小二乘逆时偏移成像方法、设备及存储介质

Similar Documents

Publication Publication Date Title
US5995905A (en) Source signature determination and multiple reflection reduction
AU2016331881B8 (en) Q-compensated full wavefield inversion
US9625593B2 (en) Seismic data processing
US20120275267A1 (en) Seismic Data Processing
US11294087B2 (en) Directional Q compensation with sparsity constraints and preconditioning
US10788597B2 (en) Generating a reflectivity model of subsurface structures
US20170031045A1 (en) Method and apparatus for modeling and separation of primaries and multiples using multi-order green's function
US11733413B2 (en) Method and system for super resolution least-squares reverse time migration
US20180059276A1 (en) System and method for focusing seismic images
CN114690242A (zh) 一种低噪音的最小二乘逆时偏移方法
Feng et al. New dynamic stochastic source encoding combined with a minmax-concave total variation regularization strategy for full waveform inversion
CN111665556B (zh) 地层声波传播速度模型构建方法
Mosher et al. Migration velocity analysis using common angle image gathers
CN110888158B (zh) 一种基于rtm约束的全波形反演方法
CN111665546B (zh) 用于可燃冰探测的声学参数获取方法
Gao et al. SREMI: Super-resolution electromagnetic imaging with single-channel ground-penetrating radar
Gholami et al. Multipliers waveform inversion
Yang et al. 3D surface-related multiples elimination based on improved apex-shifted Radon transform
CN111665549A (zh) 地层声波衰减因子反演方法
Shao et al. FWI-driven orthorhombic modeling with OBS data for reservoir imaging offshore Trinidad
Liu* et al. The research of seismic data processing sequence in gas cloud zone-applied in Bohai P oilfield
Zhu et al. Amplitude and phase in angle-domain common-image gathers
AlAli et al. High-resolution beam tomography on 3D land data applications
CN111665551B (zh) 用于桥梁基底探测的声学参数获取方法
Zhang et al. CASE STUDY-HIGH RESOLUTION VELOCITY MODEL BUILDING TO IMPROVE CARBONATE REEFS IMAGE IN TZ AREA

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