CN113156500A - 数据驱动的快速构造约束叠前地震多道反演方法 - Google Patents

数据驱动的快速构造约束叠前地震多道反演方法 Download PDF

Info

Publication number
CN113156500A
CN113156500A CN202110339971.9A CN202110339971A CN113156500A CN 113156500 A CN113156500 A CN 113156500A CN 202110339971 A CN202110339971 A CN 202110339971A CN 113156500 A CN113156500 A CN 113156500A
Authority
CN
China
Prior art keywords
seismic
inversion
data
operator
formula
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
CN202110339971.9A
Other languages
English (en)
Other versions
CN113156500B (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 East China
Original Assignee
China University of Petroleum East China
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 East China filed Critical China University of Petroleum East China
Priority to CN202110339971.9A priority Critical patent/CN113156500B/zh
Publication of CN113156500A publication Critical patent/CN113156500A/zh
Application granted granted Critical
Publication of CN113156500B publication Critical patent/CN113156500B/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/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
    • 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/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
    • G01V1/364—Seismic filtering
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/30—Noise handling
    • G01V2210/32—Noise reduction
    • G01V2210/324—Filtering
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/50—Corrections or adjustments related to wave propagation
    • G01V2210/51—Migration
    • G01V2210/512—Pre-stack
    • 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

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明公开了一种数据驱动的快速构造约束叠前地震多道反演方法。该方法包括：步骤1，估算地层倾角，并在旋转坐标系下分别构建平行于地层倾角和垂直于地层倾角的构造算子；步骤2，基于地震数据的互相关系数描述地下介质的地震反射特征，并构建数据局部优化算子；步骤3，构建数据驱动的快速构造约束叠前地震多道反演方法的目标泛函；步骤4，构建数据驱动的快速构造约束叠前地震反演的优化算法，并进一步验证反演方法的可行性。本发明从地震反射系数近似方程出发，考虑地下介质局部连续性和地震资料本身信噪比对反演结果的影响，能实现地下介质各弹性参数的稳定反演。

Description

数据驱动的快速构造约束叠前地震多道反演方法
技术领域
本发明针对油气地震勘探领域，涉及复杂地质构造区域和地震资料质量低的叠前地震多道反演，特别是高倾角构造区的弹性参数确定性预测。
背景技术
叠前地震反演作为油气地球物理勘探领域获取地下介质弹性、物性、岩性和含流体性质的关键技术，是储层定量预测和含油气性检测的重要途径。基于地震振幅随偏移距或角度变化(Amplitude variation with offset/angle)信息，叠前地震反演通过地下介质弹性、物性参数与地层反射系数的定量关系，实现储层岩性圈闭刻画和含油气性展布研究。然而，叠前地震反演具有强不确定性和不稳定性，反演参数的空间连续性有待改善。基于此，将多道地震处理技术应用于叠前地震反演的研究越来越广泛。多道地震处理技术考虑地震数据相邻地震道的相关性，能够有效压制随机噪声，提高地震解释的稳定性。多道地震处理技术在地震勘探领域已经被广泛研究，例如：噪声压制、地震资料处理和地震资料插值等。
近年来，多道地震处理技术由于其稳定性强而发展迅速，尤其是多道地震反演技术。一维横向无约束多道地震反演(laterally unconstrained inversion(1D-LUI))是一种经典的多道反演方法，与单道地震反演相比，能有效压制噪声。实际上，1D-LUI仅求解反演参数的全局最优解，未考虑相邻地震道的内在联系，反演参数的连续性有待提高。与之相比，一维横向约束多道反演技术(laterally constrained inversion(1D-LCI))具有更好的横向连续性，该技术假设地下介质为水平层状分布，通常用一个全变差算子最为横向约束。1D-LCI在地下结构简单、成层分布时，通常能取得好的应用效果，然而当地下结构复杂时，该方法容易出现地层边界过于平滑或边界模糊的现象，降低反演参数的分辨率，尤其在高倾角区域。基于此，基于地下构造或地震反射特征的横向约束方法相继被提出。目前，多道地震反演的研究趋向于发展更合理、更符合实际情况的横向约束或空间约束。然而，横向约束算子引入目标泛函通常采取降维策略，即将地震数据逐道排列成一列向量，这不可避免的会产生大型矩阵，严重影响了计算机存储和运算效率，制约其广泛实际应用。为克服这一问题，许多相应的快速算法相继发展，如Block Coordinate Dencent(BCD)快速算法，该算法的思想是在每次迭代中逐道估计局部最优解直至收敛，这使得实现大规模的地震资料多道反演相对容易。然而，BCD算法在每次迭代中循环估计各个地震道，解决了计算机存储问题，但计算效率提高并不明显，并且逐道估计局部最优解来逼近全局最优解，不能保证反演参数为全局的最优解，且地震反演的精度仍需要进一步得到改善。
发明内容
本发明要解决的就是复杂构造区域叠前地震反演缺乏稳定性和连续性的问题。
为实现上述目的，解决现有方法存在的关键问题，本发明采用下述技术方案：一种数据驱动的快速构造约束叠前地震多道反演方法，其包括如下步骤：
步骤1，估算地层倾角，并在旋转坐标系下分别构建平行于地层倾角和垂直于地层倾角的构造算子；
步骤2，基于地震数据的互相关系数描述地下介质的地震反射特征，并构建数据局部优化算子；
步骤3，构建数据驱动的快速构造约束叠前地震多道反演方法的目标泛函；
步骤4，构建数据驱动的快速构造约束叠前地震反演的优化算法。
优选的，所述步骤1包括：
将地层倾角定义为水平方向与地震数据变化梯度最小方向的夹角，则估算地层倾角的表达式为：
φ(x,z)＝tan-1(rx(x,z)/rz(x,z)) (1)，
式(1)中，φ为地层倾角，x、z分别为横向、纵向位置，rx(x,z)和rz(x,z)表示地震数据沿x、z方向的一阶偏导数；
利用地层倾角构建旋转算子的表达式为：
Figure BDA0002999164650000031
Figure BDA0002999164650000032
式(2)和(3)中，
Figure BDA0002999164650000033
和
Figure BDA0002999164650000034
为旋转坐标系下的旋转算子，
Figure BDA0002999164650000035
为第i采样点估计的地层倾角；
设Rx和Rz分别为沿x和z方向的全变差算子，将Rx和Rz旋转，得到平行和垂直于地层倾角方向的构造算子的表达式为：
Figure BDA0002999164650000036
Figure BDA0002999164650000037
式(4)和(5)中，Rparl和Rperp分别为平行和垂直于地层倾角方向的构造算子。
优选的，所述步骤2包括：
根据互相关的定义，计算互相关系数的公式C表示为：
Figure BDA0002999164650000038
式(6)中，s(i,j)表示第j道的第i个采样点的地震记录，j'为第j道的相邻地震道，w为用于相关系数计算的时窗，k表示相邻地震道内用于相关性分析的采样点与分析点的上下漂移时间，k限定在[-33]内；
构建数据局部优化算子H的表达式为：
Figure BDA0002999164650000041
式(7)中，Hi,j和Ci,j分别为数据局部优化算子H和互相关系数C在第j道的第i个采样点值，C0为互相关系数的阈值。
优选的，所述步骤3包括：
叠前地震多道正演模型表示为：
Figure BDA0002999164650000042
式(8)中，S和m为二维地震记录和待反演参数，Si为第i道地震数据，i＝1,...,n，每道地震数据包含的角度个数为h，G为正演算子，mi(i＝1,2,…,n)表示由纵波阻抗、横波阻抗和密度的自然对数组成的第i道待反演参数，Si，G和mi具体表达式为：
Figure BDA0002999164650000043
Figure BDA0002999164650000044
式(9)、(10)中，c1＝1+tan2(θ)，c2＝-8γ2tan2(θ)，
Figure BDA0002999164650000045
θi(i＝1,2,…,h)为不同角度的入射角，矩阵W表示由地震子波时移得到的子波核矩阵，W(θi)为不同角度地震子波核矩阵(i＝1,2,…,h)，IP、IS和ID分别为纵波阻抗、横波阻抗和密度的自然对数，L为一阶差分算子；
地层倾角为零的特殊情况时，地层水平成层分布，则有：
Rparl＝Rx,Rperp＝Rz (11)，
Rparl和Rperp分别为横向和纵向全变差算子，标准构造约束多道反演的目标泛函表示为：
Figure BDA0002999164650000051
式(12)中，
Figure BDA0002999164650000052
分别为式(8)中矩阵S、m展开的对角矩阵，
Figure BDA0002999164650000053
为式(10)中矩阵G的对角矩阵，λ＞0、α＞0分别为平行、垂直地层倾角的正则化参数；
基于全变差正则化的多道地震反演方法是基于地层水平成层分布，所述基于全变差正则化的多道地震反演方法对反演参数在纵向和横向上分别采用全变差算子进行约束，所述基于全变差正则化的多道地震反演方法的目标泛函为：
Figure BDA0002999164650000054
式(13)中，
Figure BDA0002999164650000055
为佛罗贝尼乌斯范数，
在式(13)的基础上，根据构造算子的定义，采用阿达玛乘积算子将构造算子引入，构建快速构造约束叠前地震多道反演方法的目标泛函表达式为：
Figure BDA0002999164650000056
式(14)中，
Figure BDA0002999164650000057
为阿达玛乘积算子，表示两个相同规模的矩阵对应元素相乘；Qcos和Qsin为与地震数据相同规模的旋转算子，其形式为：
Figure BDA0002999164650000058
Figure BDA0002999164650000059
式(15)、(16)中，
Figure BDA00029991646500000510
为第i道第j个采样点的地层倾角；
通过阿达玛乘积算子将数据局部优化算子H引入式(14)，得到数据驱动的快速构造约束叠前地震多道反演方法的目标泛函：
Figure BDA0002999164650000061
式(17)中，算子H控制了每个采样点地震数据对反演的贡献，本发明进一步根据相邻地震道反演，通过
Figure BDA0002999164650000062
项恢复地震道每个采样点处的反演参数，β＞0为地震反演横向约束的正则化算子。
优选的，所述步骤4包括：
将式(14)改写成一个约束极小化问题，表示为：
Figure BDA0002999164650000063
式(18)中，X1和X2为中间变量，表示反演参数沿纵向和横向的一阶偏导；
根据交替方向乘子算法，式(18)的拉格朗日形式为：
Figure BDA0002999164650000064
式(19)中，Z1和Z2分别为与X1和X2对应的对偶变量，采用交替方向乘子算法对式(19)进行求解，其中每个变量的更新迭代规则为：
Figure BDA0002999164650000065
Figure BDA0002999164650000066
Figure BDA0002999164650000067
Figure BDA0002999164650000068
Figure BDA0002999164650000071
根据式(20)至式(24)交替估计m、X1、X2、Z1和Z2这五个参数；
将式(17)改写成一个约束极小化问题，表示为：
Figure BDA0002999164650000072
式(25)中，
Figure BDA0002999164650000073
为中间变量，根据交替方向乘子算法将式(25)进一步写成拉格朗日的形式：
Figure BDA0002999164650000074
式(26)中，Z3为与
Figure BDA0002999164650000075
对应的对偶变量，采用交替方向乘子算法对式(26)进行求解，其中每个变量的更新迭代规则为：
Figure BDA0002999164650000076
Figure BDA0002999164650000077
Figure BDA0002999164650000078
Figure BDA0002999164650000079
Figure BDA00029991646500000710
Figure BDA00029991646500000711
Figure BDA00029991646500000712
根据式(27)至式(33)，交替估计m、X1、X2、
Figure BDA00029991646500000713
Z1、Z2和Z3这七个参数，直至目标泛函收敛。
进一步验证反演方法的可行性，可行性验证包括两个部分：
第一步：数值实验测试方法可行性，设置理论弹性模型，基于常规数值模型建立合成地震记录，首先在无噪声干扰的情况下，分别采用快速构造约束叠前地震反演方法和数据驱动的快速构造约束叠前地震反演方法进行反演，测试反演结果与真实模型的差异，然后在合成地震记录上添加噪声进行反演，测试反演方法的稳定性；
第二步：实际地震资料处理测试本发明方法的有效性，选取地质构造复杂的二维地震剖面，同样分别采用快速构造约束叠前地震反演方法和数据驱动的构造约束叠前地震反演方法进行反演，通过与实际测井数据进行对比，能够提高反演参数的稳定性和可靠性。
本发明方法从地震反射系数近似方程出发，充分考虑地下介质局部连续性和地震资料本身信噪比对反演结果的影响，提高了地下介质各弹性参数反演的稳定性、连续性和可靠性，能实现地下介质各弹性参数稳定反演。
附图说明
附图是用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例一起用于解释本发明，并不构成对本发明的限制。为了更清楚地说明本发明的具体实施过程或技术方案，附图如下：
图1本发明公开的一种数据驱动的快速构造约束叠前地震多道反演方法的流程示意图。
图2估算地层倾角和旋转坐标系下构建构造算子示意图。
图3a marmousi2纵波阻抗真实模型示意图。
图3b marmousi2横波阻抗真实模型示意图。
图4a受噪音干扰的合成小角度地震数据示意图。
图4b根据受噪音干扰的地震记录估算的地层倾角示意图。
图4c无噪音干扰的合成小角度地震记录示意图。
图4d根据无噪音干扰的地震记录估算的地层倾角示意图。
图5a基于常规构造约束反演的纵波阻抗示意图。
图5b基于常规构造约束反演的横波阻抗示意图。
图5c基于快速构造约束反演的纵波阻抗示意图。
图5d基于快速构造约束反演的横波阻抗示意图。
图5e基于数据驱动的快速构造约束反演的纵波阻抗示意图。
图5f基于数据驱动的快速构造约束反演的横波阻抗示意图。
图6a常规构造约束反演方法相对误差随迭代次数的变化示意图。
图6b常规构造约束反演、快速构造约束反演和数据驱动的快速构造约束反演方法相对误差对比示意图。
图7a实际资料基于常规构造约束反演的纵波阻抗示意图。
图7b实际资料基于常规构造约束反演的横波阻抗示意图。
图7c实际资料基于快速构造约束反演的纵波阻抗示意图。
图7d实际资料基于快速构造约束反演的横波阻抗示意图。
图7e实际资料基于数据驱动的快速构造约束反演的纵波阻抗示意图。
图7f实际资料基于数据驱动的快速构造约束反演的横波阻抗示意图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
本发明充分考虑地下介质局部连续性和地震资料本身信噪比对反演结果的影响，研发了一种数据驱动的快速构造约束叠前地震多道反演技术(详见图1流程示意图)，提高了地下介质各弹性参数反演的稳定性、连续性和可靠性。
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
本发明一种数据驱动的快速构造约束叠前地震多道反演方法，包括如下步骤：
步骤1，估算地层倾角，并在旋转坐标系下分别构建平行于地层倾角和垂直于地层倾角的构造算子；
步骤2，基于地震数据的互相关系数描述地下介质的地震反射特征，并构建数据局部优化算子；
步骤3，构建数据驱动的快速构造约束叠前地震多道反演方法的目标泛函；
步骤4，构建数据驱动的的快速构造约束叠前地震反演的优化算法。
步骤1中，将地层倾角定义为水平方向与地震数据变化梯度最小方向的夹角，则估算地层倾角的表达式为：
φ(x,z)＝tan-1(rx(x,z)/rz(x,z)) (1)，
式(1)中，φ为地层倾角，x、z分别为横向、纵向位置，rx(x,z)和rz(x,z)表示地震数据沿x、z方向的一阶偏导数。这种计算地层倾角的方法对信噪比高的地震数据通常能够准确的估计地层倾角，而当地震数据信噪比低(信噪比小于4)时，对地层构造描述精度有待改善，该种方法反正切函数会放大噪声的影响，对地震数据质量有一定的依赖性。通常采用平滑的方式，降低噪声影响，但这是以牺牲一些细节构造信息为代价。
本发明对地震多道反演添加平行和垂直地层倾角方向的构造约束，以保留反演参数的构造特征，这通常通过构造算子实现。估算地层倾角就是为了构建更合理的构造算子，在估算地层倾角之后，利用旋转坐标系(图2)可以构建构造算子。图中Xr、Zr坐标轴为旋转之后的坐标轴，旋转角度为地层倾角φ。
构造算子的构建可以从一阶全变差算子Rx和Rz(沿x和z方向的全变差算子)出发，利用地层倾角可以构建如下的旋转算子：
Figure BDA0002999164650000101
Figure BDA0002999164650000111
式(2)和(3)中，
Figure BDA0002999164650000112
和
Figure BDA0002999164650000113
为旋转坐标系下的旋转算子，
Figure BDA0002999164650000114
为第i采样点估计的地层倾角。将Rx和Rz旋转，可以得到平行和垂直于地层倾角方向的构造算子，其具体形式为：
Figure BDA0002999164650000115
Figure BDA0002999164650000116
式(4)和(5)中，Rparl和Rperp分别为平行和垂直于地层倾角方向的构造算子。构造算子的精度直接关系到反演参数的稳定性和准确性。在估算地层倾角之前进行随机噪声压制、连续性处理等，可以有效地提高倾角估计的准确性。
步骤2中，基于地震数据互相关评价地震数据每个采样点的可靠性，利用相关系数控制每个采样点地震数据对反演参数的贡献，以降低噪音干扰严重、极度不连续或死道等采样点处地震数据对反演的贡献，提高地震反演的可靠性和连续性。根据互相关的定义，计算互相关系数的公式C表示为：
Figure BDA0002999164650000117
式(6)中，s(i,j)表示第j道的第i个采样点的地震记录，j′为第j道的相邻地震道，w为用于相关系数计算的时窗，k表示相邻地震道内用于相关性分析的采样点与分析点的上下漂移时间，根据地震数据在局部一般不会发生剧烈抖动，将k限定在[-3 3]内。本发明利用互相关系数控制每个采样点地震数据对反演的贡献，这一过程通过数据局部优化算子H实现，算子H的表达式为：
Figure BDA0002999164650000118
式(7)中，Hi,j和Ci,j分别为数据局部优化算子H和互相关系数C在第j道的第i个采样点值，C0为互相关系数的阈值。将互相关系数小于该值的采样点视为不可靠的点，控制该点地震数据对反演参数不做贡献，根据相邻地震道估计这些采样点的反演参数。
步骤3中，常规多道地震反演方法易产生大型矩阵，不利于计算机运算求解。本发明为避免大型矩阵的产生，提高反演效率，发展了一种引入构造算子的方法。
以反演纵波阻抗、横波阻抗和密度参数为例，叠前地震多道正演模型通常可以表示为：
Figure BDA0002999164650000121
式(8)中，S和m为二维地震记录和待反演参数，Si为第i道地震数据i＝1,...,n，每道地震数据包含的角度个数为h，G为正演算子，mi(i＝1,2,…,n)表示由纵波阻抗、横波阻抗和密度的第i道自然对数组成的待反演参数。Si，G和mi具体表达式为：
Figure BDA0002999164650000122
Figure BDA0002999164650000123
式(9)、(10)中，c1＝1+tan2(θ)，c2＝-8γ2tan2(θ)，
Figure BDA0002999164650000124
θi(i＝1,2,…,h)为不同角度的入射角，矩阵W表示由地震子波时移得到的子波核矩阵，W(θi)为不同角度地震子波核矩阵(i＝1,2,…,h)，IP、IS和ID分别为纵、横波阻抗和密度的自然对数，L为一阶差分算子。
首先考虑地层倾角为零
Figure BDA0002999164650000125
的特殊情况，构造算子为：
Rparl＝Rx,Rperp＝Rz (11)
即地层水平成层分布，构造算子分别为横向和纵向全变差算子。此时，标准构造约束多道反演的目标泛函可以表示为：
Figure BDA0002999164650000131
式(12)中，
Figure BDA0002999164650000132
分别为式(8)中矩阵S、m展开的对角矩阵，
Figure BDA0002999164650000133
为式(10)中矩阵G的对角矩阵，λ＞0、α＞0分别为平行、垂直地层倾角的正则化参数。
然而，基于全变差正则化(Total variation regularization,简称TV正则化)的多道地震反演方法是基于地层水平成层分布，该方法对反演参数在纵向和横向上分别采用全变差算子进行约束，其反演的目标泛函为：
Figure BDA0002999164650000134
式(13)中，
Figure BDA0002999164650000135
为佛罗贝尼乌斯范数(Frobenius范数，简称F-范数)。可见，式(12)和式(13)是等价的，不同的是，式(13)具有极高的反演效率。本发明在式(13)的基础上，根据构造算子的定义，采用阿达玛乘积算子(即Hadamard乘积算子)将构造算子引入目标泛函，其表达式为：
Figure BDA0002999164650000136
式(14)中，“
Figure BDA0002999164650000137
”为Hadamard乘积算子，表示两个相同规模的矩阵对应元素相乘。Qcos和Qsin为与地震数据相同规模的旋转算子，式(2)和式(3)中的
Figure BDA0002999164650000138
和
Figure BDA0002999164650000139
实际为矩阵Qcos和Qsin展开的对角阵，
Figure BDA00029991646500001310
和
Figure BDA00029991646500001311
矩阵规模远大于Qcos和Qsin，严重影响反演效率，Qcos和Qsin的形式为：
Figure BDA00029991646500001312
Figure BDA00029991646500001313
式(15)、(16)中，
Figure BDA0002999164650000141
为第i道第j个采样点的地层倾角。
式(14)为本发明发展的快速构造约束叠前地震多道反演方法的目标泛函，该方法利用二维多道正演模型的运算优势，在全变差正则化的基础上采用Hadamard乘积算子引入构造算子，避免了大型矩阵的产生。
通过Hadamard乘积算子将数据局部优化算子H引入目标泛函(14)，可得到：
Figure BDA0002999164650000142
式(17)中，算子H控制了每个采样点地震数据对反演的贡献，本发明进一步根据相邻地震道反演，通过
Figure BDA0002999164650000143
项恢复地震道每个采样点处的反演参数，β＞0为地震反演横向约束的正则化算子。式(17)即为数据驱动的快速构造约束叠前地震多道反演方法的目标泛函，该方法能有效降低地震数据质量对反演结果的影响。
步骤4中，直接求解式(17)具有一定困难，首先研究式(14)的求解策略。从形式上看，快速构造约束多道地震反演的目标泛函(14)与常规多道地震反演的目标泛函是不同的，它包括Frobenius范数和Hadamard乘积算子两种特殊运算符，这增加了目标泛函的求解难度。本发明从两种运算符的基本定义出发，基于交替方向乘子算法，即ADMM(Thealternating direction method of multipliers,简称ADMM)算法，详细推导了其求解过程。
直接对式(14)进行求解是非常困难的，本发明首先将其改写成一个约束极小化问题，可以表示为：
Figure BDA0002999164650000144
式(18)中，X1和X2为中间变量，表示反演参数沿纵向和横向的一阶偏导。根据交替方向乘子算法，式(18)的拉格朗日形式为：
Figure BDA0002999164650000151
式(19)中，Z1和Z2分别为与X1和X2对应的对偶变量，交替方向乘子算法在对偶交替估计过程中使反问题的解逐步逼近最优解，具有更高的计算效率和精度。采用交替方向乘子算法对式(19)进行求解，其中每个变量的更新迭代规则为：
Figure BDA0002999164650000152
Figure BDA0002999164650000153
Figure BDA0002999164650000154
Figure BDA0002999164650000155
Figure BDA0002999164650000156
根据式(20)～(24)交替估计m、X1、X2、Z1和Z2这五个参数，通常就能够收敛。该算法不涉及大型矩阵，且在每次迭代中对所有地震道同时估计，避免了块坐标下降算法逐道估计局部极小值的不足，提高了反演效率和精度。
快速构造约束叠前地震多道反演算法的伪代码如下所示：
Figure BDA0002999164650000157
Figure BDA0002999164650000161
数据驱动的快速构造约束叠前地震多道反演的目标泛函同样包含Frobenius范数和Hadamard乘积算子两种特殊运算符，难以直接求解。与式(14)的求解过程一样，本发明基于交替方向乘子算法详细推导了求解过程。
首先将式(17)改写成一个约束极小化问题，可以表示为：
Figure BDA0002999164650000162
式(25)中，
Figure BDA0002999164650000163
为中间变量，根据交替方向乘子算法将式(25)进一步写成拉格朗日的形式：
Figure BDA0002999164650000164
式(26)中，Z3为与
Figure BDA0002999164650000165
对应的对偶变量，采用交替方向乘子算法对式(26)进行求解，其中每个变量的更新迭代规则为：
Figure BDA0002999164650000166
Figure BDA0002999164650000167
Figure BDA0002999164650000171
Figure BDA0002999164650000172
Figure BDA0002999164650000173
Figure BDA0002999164650000174
Figure BDA0002999164650000175
该算法在快速构造约束多道反演方法的基础上降低了地震资料质量的影响。根据上述公式交替估计m、X1、X2、
Figure BDA0002999164650000176
Z1、Z2和Z3这七个参数，直至目标泛函收敛。
数据驱动的快速构造约束叠前地震多道反演算法的伪代码如下所示：
Figure BDA0002999164650000177
Figure BDA0002999164650000181
进一步验证反演方法的可行性，可行性验证包括两个步骤：
第一步：数值实验测试方法可行性，设置理论弹性模型，基于常规数值模型建立合成地震记录，首先在无噪声干扰的情况下，分别采用快速构造约束叠前地震反演方法和数据驱动的快速构造约束叠前地震反演方法进行反演，测试反演结果与真实模型的差异，然后在合成地震记录上添加噪声进行反演，测试反演方法的稳定性；
第二步：实际地震资料处理测试本发明方法的有效性，选取地质构造复杂的二维地震剖面，同样分别采用快速构造约束叠前地震反演方法和数据驱动的构造约束叠前地震反演方法进行反演，通过与实际测井数据进行对比，能够提高反演参数的稳定性和可靠性。
为进一步说明本发明的可行性和有效性，下面列举两个实施例：
实施例1：数值模型测试，详见图3、图4、图5和图6。
本发明选取部分Marmousi-2纵、横波阻抗(图3a和图3b)和密度模型对该方法进行测试。用于测试的真实模型共包含851道，每道606个采样点，采样间隔为2ms。利用30Hz雷克子波和真实阻抗模型合成地震记录，在合成地震记录中随机去除50％地震道的部分地震数据，并添加10％的高斯噪声以测试反演方法的稳定性。由于地震数据具有带限性质，反演结果缺乏低频成分，本文将真实阻抗模型进行0-5Hz低通滤波得到低频阻抗模型，用来补充反演结果的低频成分。
本文基于叠后地震记录(图4c)估计地层倾角和构建构造算子。本发明提出的数据驱动的快速构造约束叠前地震多道反演方法需利用互相关来评价地震数据的合理性，首先对合成的角度道集进行相关性分析，并计算互相关系数。为方便成图，现有的常规构造约束反演方法可以简称为SCI(Structurally Constrained Inversion)，本发明将提出的快速构造约束反演简称为FSCI(Fast Structurally Constrained Inversion)，其目标泛函表示为式(17)，本发明进一步提出了数据驱动的快速构造约束叠前地震多道反演简称为DFSCI(Data-driven Fast Structurally Constrained Inversion)，其目标泛函表示为式(25)。将这三种方法的反演连续性、可靠性和效率分别进行了对比分析。
图5为三种方法纵、横波阻抗的反演结果。从整体上看，与真实阻抗模型相比，三种方法均与真实模型相近，表明了本发明反演测试工作是可靠的。另外，常规构造约束反演方法(SCI)与快速构造约束反演方法(FSCI)反演结果基本一样，构造特征明显，地层边界清晰，尤其在高倾角区域，这说明快速构造约束反演能够保留常规构造约束反演方法的优势。与前两种方法相比，数据驱动的快速构造约束反演方法(DFSCI)具有更高的连续性和信噪比(图5e和图5f)，表明本发明提出的数据驱动的快速构造约束叠前地震多道反演方法能够保证反演参数的稳定性、可靠性和连续性。
图6对三种方法的反演效率进行了对比，从图6a中可看出，常规构造约束反演的收敛效率比较低，需要数十次甚至上百次迭代才能收敛，且从图6b可以看出，该方法迭代一次需要更长的时间，而本发明提出的两种方法收敛速度极快，收敛误差小，也说明了这两种方法具有高的计算精度和效率。
实施例2：实际资料处理，详见图7。
本发明在常规构造约束反演方法(SCI)的基础上发展了快速构造约束反演方法(FSCI)，为改善对地震资料质量的影响又进一步提出了数据驱动的快速构造约束反演方法(DFSCI)。并且分别利用三种方法进行了叠前地震反演测试。图7a和7b为常规构造约束反演方法反演的纵横波阻抗，可以看出，该方法能够保证地层边界清晰，构造特征明显。图7c和7d为本文提出的快速构造约束反演方法反演的纵横波阻抗，与常规方法相比，该方法的构造特征和地层边界更为清晰，这是因为该方法为全局优化算法，具有更高的计算精度。图7e和7f在保证地层边界和构造特征的同时，进一步提高了反演参数的横向连续性和信噪比，这表明本文提出的两种方法具有更好的实际应用效果。

Claims (5)

1.一种数据驱动的快速构造约束叠前地震多道反演方法，其特征在于，其包括如下步骤：
步骤1，估算地层倾角，并在旋转坐标系下分别构建平行于地层倾角和垂直于地层倾角的构造算子；
步骤2，基于地震数据的互相关系数描述地下介质的地震反射特征，并构建数据局部优化算子；
步骤3，构建数据驱动的快速构造约束叠前地震多道反演方法的目标泛函；
步骤4，构建数据驱动的快速构造约束叠前地震反演的优化算法。
2.如权力要求1所述数据驱动的快速构造约束叠前地震多道反演方法，其特征在于：
所述步骤1包括：
将地层倾角定义为水平方向与地震数据变化梯度最小方向的夹角，则估算地层倾角的表达式为：
φ(x,z)＝tan-1(rx(x,z)/rz(x,z)) (1)，
式(1)中，φ为地层倾角，x、z分别为横向、纵向位置，rx(x,z)和rz(x,z)表示地震数据沿x、z方向的一阶偏导数；
利用地层倾角构建旋转算子的表达式为：
Figure FDA0002999164640000011
Figure FDA0002999164640000012
式(2)和(3)中，
Figure FDA0002999164640000013
和
Figure FDA0002999164640000014
为旋转坐标系下的旋转算子，
Figure FDA0002999164640000015
为第i采样点估计的地层倾角；
设Rx和Rz分别为沿x和z方向的全变差算子，将Rx和Rz旋转，得到平行和垂直于地层倾角方向的构造算子的表达式为：
Figure FDA0002999164640000021
Figure FDA0002999164640000022
式(4)和(5)中，Rparl和Rperp分别为平行和垂直于地层倾角方向的构造算子。
3.如权力要求2所述数据驱动的快速构造约束叠前地震多道反演方法，其特征在于：
所述步骤2包括：
根据互相关的定义，计算互相关系数的公式C表示为：
Figure FDA0002999164640000023
式(6)中，s(i,j)表示第j道的第i个采样点的地震记录，j'为第j道的相邻地震道，w为用于相关系数计算的时窗，k表示相邻地震道内用于相关性分析的采样点与分析点的上下漂移时间，k限定在[-3 3]内；
构建数据局部优化算子H的表达式为：
Figure FDA0002999164640000024
式(7)中，Hi,j和Ci,j分别为数据局部优化算子H和互相关系数C在第j道的第i个采样点值，C0为互相关系数的阈值。
4.如权力要求3所述数据驱动的快速构造约束叠前地震多道反演方法，其特征在于：
所述步骤3包括：
叠前地震多道正演模型表示为：
Figure FDA0002999164640000031
式(8)中，S和m为二维地震记录和待反演参数，Si(i＝1,…,n)为第i道地震数据，每道地震数据包含的角度个数为h，G为正演算子，mi(i＝1,2,…,n)表示由纵波阻抗、横波阻抗和密度的自然对数组成的第i道待反演参数，Si，G和mi具体表达式为：
Figure FDA0002999164640000032
Figure FDA0002999164640000033
式(9)、(10)中，c1＝1+tan2(θ)，c2＝-8γ2tan2(θ)，
Figure FDA0002999164640000034
θi(i＝1,2,…,h)为不同角度的入射角，矩阵W表示由地震子波时移得到的子波核矩阵，W(θi)(i＝1,2,…,h)为不同角度地震子波核矩阵，IP、IS和ID分别为纵波阻抗、横波阻抗和密度的自然对数，L为一阶差分算子；
地层倾角为零的特殊情况时，地层水平成层分布，则有：
Rparl＝Rx,Rperp＝Rz (11)，
Rparl和Rperp分别为横向和纵向全变差算子，标准构造约束多道反演的目标泛函表示为：
Figure FDA0002999164640000035
式(12)中，
Figure FDA0002999164640000036
分别为式(8)中矩阵S、m展开的对角矩阵，
Figure FDA0002999164640000037
为式(10)中矩阵G的对角矩阵，λ＞0、α＞0分别为平行、垂直地层倾角的正则化参数；
基于全变差正则化的多道地震反演方法是基于地层水平成层分布，所述基于全变差正则化的多道地震反演方法对反演参数在纵向和横向上分别采用全变差算子进行约束，所述基于全变差正则化的多道地震反演方法的目标泛函为：
Figure FDA0002999164640000041
式(13)中，
Figure FDA0002999164640000042
为佛罗贝尼乌斯范数，
在式(13)的基础上，根据构造算子的定义，采用阿达玛乘积算子将构造算子引入，构建快速构造约束叠前地震多道反演方法的目标泛函表达式为：
Figure FDA0002999164640000043
式(14)中，
Figure FDA0002999164640000044
为阿达玛乘积算子，表示两个相同规模的矩阵对应元素相乘；Qcos和Qsin为与地震数据相同规模的旋转算子，其形式为：
Figure FDA0002999164640000045
Figure FDA0002999164640000046
式(15)、(16)中，
Figure FDA0002999164640000047
为第i道第j个采样点的地层倾角；通过阿达玛乘积算子将数据局部优化算子H引入式(14)，得到数据驱动的快速构造约束叠前地震多道反演方法的目标泛函：
Figure FDA0002999164640000048
式(17)中，算子H控制了每个采样点地震数据对反演的贡献，进一步根据相邻地震道反演，通过
Figure FDA0002999164640000049
项恢复地震道每个采样点处的反演参数，β＞0为地震反演横向约束的正则化算子。
5.如权力要求4所述数据驱动的快速构造约束叠前地震多道反演方法，其特征在于：
所述步骤4包括：
将式(14)改写成一个约束极小化问题，表示为：
Figure FDA0002999164640000051
式(18)中，X1和X2为中间变量，表示反演参数沿纵向和横向的一阶偏导；
根据交替方向乘子算法，式(18)的拉格朗日形式为：
Figure FDA0002999164640000052
式(19)中，Z1和Z2分别为与X1和X2对应的对偶变量，采用交替方向乘子算法对式(19)进行求解，其中每个变量的更新迭代规则为：
Figure FDA0002999164640000053
Figure FDA0002999164640000054
Figure FDA0002999164640000055
Figure FDA0002999164640000056
Figure FDA0002999164640000057
根据式(20)至式(24)交替估计m、X1、X2、Z1和Z2这五个参数；
将式(17)改写成一个约束极小化问题，表示为：
Figure FDA0002999164640000061
式(25)中，
Figure FDA0002999164640000062
为中间变量，根据交替方向乘子算法将式(25)进一步写成拉格朗日的形式：
Figure FDA0002999164640000063
式(26)中，Z3为与
Figure FDA0002999164640000064
对应的对偶变量，采用交替方向乘子算法对式(26)进行求解，其中每个变量的更新迭代规则为：
Figure FDA0002999164640000065
Figure FDA0002999164640000066
Figure FDA0002999164640000067
Figure FDA0002999164640000068
Figure FDA0002999164640000069
Figure FDA00029991646400000610
Figure FDA00029991646400000611
根据式(27)至式(33)，交替估计m、X1、X2、
Figure FDA00029991646400000612
Z1、Z2和Z3这七个参数，直至目标泛函收敛；
进一步验证反演方法的可行性，可行性验证包括两个部分：
第一步：数值实验测试方法可行性，设置理论弹性模型，基于常规数值模型建立合成地震记录，首先在无噪声干扰的情况下，分别采用快速构造约束叠前地震反演方法和数据驱动的快速构造约束叠前地震反演方法进行反演，测试反演结果与真实模型的差异，然后在合成地震记录上添加噪声进行反演，测试反演方法的稳定性；
第二步：实际地震资料处理测试本发明方法的有效性，选取地质构造复杂的二维地震剖面，同样分别采用快速构造约束叠前地震反演方法和数据驱动的构造约束叠前地震反演方法进行反演，通过与实际测井数据进行对比，能够提高反演参数的稳定性和可靠性。
CN202110339971.9A 2021-03-30 2021-03-30 数据驱动的快速构造约束叠前地震多道反演方法 Active CN113156500B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110339971.9A CN113156500B (zh) 2021-03-30 2021-03-30 数据驱动的快速构造约束叠前地震多道反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110339971.9A CN113156500B (zh) 2021-03-30 2021-03-30 数据驱动的快速构造约束叠前地震多道反演方法

Publications (2)

Publication Number Publication Date
CN113156500A true CN113156500A (zh) 2021-07-23
CN113156500B CN113156500B (zh) 2022-08-23

Family

ID=76885666

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110339971.9A Active CN113156500B (zh) 2021-03-30 2021-03-30 数据驱动的快速构造约束叠前地震多道反演方法

Country Status (1)

Country Link
CN (1) CN113156500B (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113917538A (zh) * 2021-09-18 2022-01-11 中国石油大学(华东) 一种地震波形约束的叠前地震反演方法、设备及存储介质
CN114966878A (zh) * 2022-04-12 2022-08-30 成都理工大学 一种基于混合范数和互相关约束的三维重力场反演方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103293551A (zh) * 2013-05-24 2013-09-11 中国石油天然气集团公司 一种基于模型约束的阻抗反演方法及系统
US20160116637A1 (en) * 2013-06-04 2016-04-28 Total S.A. Method of constraining seismic inversion
CN111694056A (zh) * 2020-06-03 2020-09-22 西安交通大学 一种压制地震资料异常噪声的方法、存储介质及设备
CN112162318A (zh) * 2020-09-29 2021-01-01 地球脉动(无锡)科技有限公司 基于倾角约束的多道反褶积处理方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103293551A (zh) * 2013-05-24 2013-09-11 中国石油天然气集团公司 一种基于模型约束的阻抗反演方法及系统
US20160116637A1 (en) * 2013-06-04 2016-04-28 Total S.A. Method of constraining seismic inversion
CN111694056A (zh) * 2020-06-03 2020-09-22 西安交通大学 一种压制地震资料异常噪声的方法、存储介质及设备
CN112162318A (zh) * 2020-09-29 2021-01-01 地球脉动(无锡)科技有限公司 基于倾角约束的多道反褶积处理方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
HAITHAM HAMID ET AL.: "Prestack structurally constrained impedance inversion", 《GEOPHYSICS》 *
印兴耀等: "地震数据互相关驱动的多道反演方法", 《地球物理学报》 *
程三: "基于构造约束的地震信号多道反演方法研究", 《中国优秀博硕士学位论文全文数据库（硕士) 基础科学辑》 *

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113917538A (zh) * 2021-09-18 2022-01-11 中国石油大学(华东) 一种地震波形约束的叠前地震反演方法、设备及存储介质
CN114966878A (zh) * 2022-04-12 2022-08-30 成都理工大学 一种基于混合范数和互相关约束的三维重力场反演方法

Also Published As

Publication number Publication date
CN113156500B (zh) 2022-08-23

Similar Documents

Publication Publication Date Title
US10295683B2 (en) Amplitude inversion on partitioned depth image gathers using point spread functions
Wu et al. Horizon volumes with interpreted constraints
US8363509B2 (en) Method for building velocity models for pre-stack depth migration via the simultaneous joint inversion of seismic, gravity and magnetotelluric data
US20190353814A1 (en) Direct Migration of Simultaneous-Source Survey Data
US9869783B2 (en) Structure tensor constrained tomographic velocity analysis
CA2683618C (en) Inverse-vector method for smoothing dips and azimuths
US11294087B2 (en) Directional Q compensation with sparsity constraints and preconditioning
CN113156500B (zh) 数据驱动的快速构造约束叠前地震多道反演方法
US10310117B2 (en) Efficient seismic attribute gather generation with data synthesis and expectation method
CN111522063B (zh) 基于分频联合反演的叠前高分辨率流体因子反演方法
Yang et al. Data-driven fast prestack structurally constrained inversion
Sun et al. Intelligent ava inversion using a convolution neural network trained with pseudo-well datasets
CN114861515A (zh) 层速度数据体的计算方法、装置、设备及介质
CN110208858B (zh) 基于叠前反演的“甜点”概率直接估算方法及系统
US20230140656A1 (en) Method and system for determining seismic processing parameters using machine learning
Yang et al. Fast spatial structure constrained 3D prestack inversion based on the Sylvester equation
RU2126984C1 (ru) Способ определения глубинно-скоростных параметров среды и построения ее изображения по сейсмическим данным - система prime
CN111239805B (zh) 基于反射率法的块约束时移地震差异反演方法及系统
Li et al. Progressive multitask learning for high-resolution prediction of reservoir elastic parameters
CN111308550B (zh) 一种页岩vti储层的各向异性参数多波联合反演方法
CN114740528A (zh) 一种超微分拉普拉斯块约束的叠前多波联合反演方法
CN107678065A (zh) 提高地震分辨率的保构造井控空间反褶积方法和装置
Li et al. Seismic full-waveform inversion based on superresolution for high-precision prediction of reservoir parameters
CN112147700A (zh) 速度异常区的低频模型构建方法及系统
CN112462428B (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