CN113064203A - 共轭梯度归一化lsrtm方法、系统、存储介质及应用 - Google Patents
共轭梯度归一化lsrtm方法、系统、存储介质及应用 Download PDFInfo
- Publication number
- CN113064203A CN113064203A CN202110323886.3A CN202110323886A CN113064203A CN 113064203 A CN113064203 A CN 113064203A CN 202110323886 A CN202110323886 A CN 202110323886A CN 113064203 A CN113064203 A CN 113064203A
- Authority
- CN
- China
- Prior art keywords
- migration
- seismic
- conjugate gradient
- profile
- imaging
- 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
Links
- 238000000034 method Methods 0.000 title claims abstract description 72
- 238000010606 normalization Methods 0.000 title claims abstract description 38
- 238000003860 storage Methods 0.000 title claims abstract description 9
- 238000013508 migration Methods 0.000 claims abstract description 202
- 230000005012 migration Effects 0.000 claims abstract description 200
- 230000002441 reversible effect Effects 0.000 claims abstract description 143
- 238000003384 imaging method Methods 0.000 claims abstract description 120
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 11
- 238000004364 calculation method Methods 0.000 claims description 36
- 238000012545 processing Methods 0.000 claims description 16
- 238000002939 conjugate gradient method Methods 0.000 claims description 15
- 239000011159 matrix material Substances 0.000 claims description 14
- 238000001514 detection method Methods 0.000 claims description 7
- 238000004088 simulation Methods 0.000 claims description 5
- 230000036961 partial effect Effects 0.000 claims description 4
- 238000004590 computer program Methods 0.000 claims description 3
- 238000012937 correction Methods 0.000 claims description 3
- 238000012546 transfer Methods 0.000 claims description 3
- 230000017105 transposition Effects 0.000 claims 1
- 230000000694 effects Effects 0.000 abstract description 13
- 238000005286 illumination Methods 0.000 abstract description 6
- 230000006872 improvement Effects 0.000 abstract description 3
- 238000010586 diagram Methods 0.000 description 8
- 238000001914 filtration Methods 0.000 description 4
- 238000011160 research Methods 0.000 description 4
- 230000008901 benefit Effects 0.000 description 3
- 230000008859 change Effects 0.000 description 3
- 238000004519 manufacturing process Methods 0.000 description 3
- 238000003491 array Methods 0.000 description 2
- 238000009499 grossing Methods 0.000 description 2
- 238000004321 preservation Methods 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- 230000001629 suppression Effects 0.000 description 2
- 230000009471 action Effects 0.000 description 1
- 230000006978 adaptation Effects 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000009826 distribution Methods 0.000 description 1
- 230000000482 effect on migration Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 230000000670 limiting effect Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 230000002829 reductive effect Effects 0.000 description 1
- 239000004065 semiconductor Substances 0.000 description 1
- 230000003313 weakening effect Effects 0.000 description 1
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
- 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
-
- 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/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Remote Sensing (AREA)
- General Physics & Mathematics (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- Geophysics (AREA)
- Data Mining & Analysis (AREA)
- Computational Mathematics (AREA)
- Mathematical Physics (AREA)
- Theoretical Computer Science (AREA)
- Pure & Applied Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Operations Research (AREA)
- Probability & Statistics with Applications (AREA)
- Evolutionary Biology (AREA)
- Algebra (AREA)
- Bioinformatics & Computational Biology (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本发明涉及共轭梯度归一化LSRTM方法、系统、存储介质及应用,属于地球物理勘探技术领域,所述方法将常规逆时偏移结果作为最小二乘反演的输入,对其进行波恩近似正演,改进迭代算法中的梯度求取公式,提高算法收敛速度,使用震源归一化共轭度计算更新步长,不断迭代输出满足收敛条件的偏移剖面。该成像方法既能够减弱地震震源效应影响,降低偏移剖面的低频干扰,提高深部照明能量和成像精度。又通过改进使得算法收敛速度相比常规最小二乘逆时偏移算法有了大幅度提高,从而节省了大量的计算资源和计算成本。本发明为下一步的地震勘探解释工作提供了高效率、高精度的偏移成像服务。
Description
技术领域
本发明属于地球物理勘探技术领域,尤其涉及一种共轭梯度归一化LSRTM方法、系统、存储介质及应用。
背景技术
随着油气勘探的深入,勘探目标逐渐由简单构造向复杂构造演变,偏移成像方法技术也随之逐步发展,从最初的基于射线理论的Kirchhoff偏移和Beam偏移,到基于单程波和双程波的逆时偏移,再到基于反演理论的最小二乘逆时偏移。目前针对复杂构造的偏移成像方法中,基于双程波波动方程的逆时偏移方法被人们广泛的使用。相较其他常规偏移方法,它具有以下几个优点,首先就是它不受倾角的限制,可以对多种特殊波进行成像,如转换波、多次波、回转波等。其次它不受介质速度变化的影响,可以适应速度的横向变化,对于地层中的局部特殊构造也可以较好的成像。但是,逆时偏移的主要问题是:其偏移算子只是正演算子的共轭转置,并不是其逆,另外由于采集孔径有限、地下模型复杂以及地震波带宽有限等的影响,逆时偏移仅能够提供较模糊的构造信息,无法对复杂油气藏进行精细成像。
而最小二乘逆时偏移方法(Least squares inverse time migration method,简称LSRTM)是基于线性反演理论的真振幅成像方法,该方法可以克服逆时偏移技术对深部储层成像分辨率不高且振幅不均衡的问题,能够得到更精确的成像剖面。但是,常规最小二乘逆时偏移一般使用互相关共轭梯度法进行迭代,这种迭代方式受震源效应影响较大,偏移成像结果存在低频干扰,其深部成像振幅较弱,很难识别深部的有效能量信息。同时,其收敛速度较慢,会导致其计算量和存储需求都非常大,这会大大增加处理资料的成本。
发明内容
本发明要解决的技术问题在于提供一种共轭梯度归一化LSRTM方法、系统、存储介质及应用。在常规互相关共轭梯度法最小二乘逆时偏移的基础上,所述方法将最小二乘反演思路与逆时偏移相结合,首先将常规逆时偏移当作最小二乘反演的第一次结果,然后对逆时偏移成像结果进行波恩近似正演,不断迭代更新成像结果,同时在迭代的过程中,采用归一化思想对梯度算法进行优化,最终实现对复杂构造的精细成像。
本发明是通过如下技术方案来实现的:
一种共轭梯度归一化LSRTM方法,具体步骤如下:
完成地震偏移速度场数据和地震炮记录数据输入,对输入地震数据进行逆时偏移成像处理:首先基于输入的偏移速度场进行震源波场正向延拓,对地震炮记录数据进行检波波场反向延拓,得到地震波在地下传播每一时刻的正向延拓波场值和反向延拓波场值,再利用最小二乘逆时偏移互相关成像条件得到逆时偏移剖面;
进行波恩近似正演:将最小二乘反演思路与逆时偏移相结合,将常规逆时偏移结果当作最小二乘反演的第一次结果,对逆时偏移剖面进行波恩近似正演得到波恩近似正演的波场记录;
进行地震数据残差的计算:使用波恩近似正演的波场记录与实际输入的地震波场相减,得到地震数据最小二乘逆时偏移计算的残差值;
进行共轭梯度的计算:使用残差值,在常规互相关共轭梯度法最小二乘逆时偏移的基础上,采用共轭梯度归一化方法进行共轭梯度的计算;
进行更新步长的计算并完成模型的更新:使用共轭梯度值计算更新步长,并在原来偏移剖面的基础上进行模型的更新,完成归一化互相关成像条件最小二乘逆时偏移的一次迭代;
判断偏移结果是否满足输出偏移成像收敛条件:当偏移成像剖面结果满足收敛条件时,输出地震偏移剖面,如果偏移剖面不满足收敛条件时,对偏移剖面结果进行波恩近似正演并继续迭代,直到满足收敛条件时输出偏移剖面。
进一步,所述的进行波恩近似正演,使用以下公式:
d=Lm
其中m为偏移剖面或反射系数模型的矩阵形式,d为模拟数据的矩阵形式,L为波恩近似正传算子矩阵。
进一步,所述的进行共轭梯度的计算方法:
式中RS是地震记录反传波场,是震源正传波场一阶偏导使用残差值,Rres是波恩近似正演数据与采集数据残差的反传波场,为最速下降梯度,是归一化共轭梯度法修正因子,是归一化共轭梯度算子,t代表地震波传播时间,x代表地震模型横向长度,z代表地震模型纵向长度,k代表当前迭代次数,k+1代表下一次迭代次数。
进一步,所述共轭梯度归一化最小二乘法逆时偏移成像方法进行进行更新步长的计算并完成模型的更新,迭代算法如下:
式中是更新步长,m(k+1)是迭代更新后偏移剖面,m(k)是迭代更新前偏移剖面,使用共轭梯度值计算更新步长,并在原来偏移剖面的基础上进行模型的更新,是归一化共轭梯度算子,为最速下降梯度,T为矩阵转置标志。k代表当前迭代次数,k+1代表下一次迭代次数、L为波恩近似正传算子矩阵,是归一化共轭梯度算子,为最速下降梯度。
进一步,所述判断偏移结果是否满足输出偏移成像收敛条件:在地震资料计算中,设定收敛条件为迭代次数,当满足迭代次数后进行偏移成像剖面的输出。
本发明的另一目的在于提供一种计算机可读存储介质,存储有计算机程序,所述计算机程序被处理器执行时,使得所述处理器执行所述一种共轭梯度归一化LSRTM方法。
本发明的另一目的在于提供一种实施所述共轭梯度归一化LSRTM方法的系统,所述系统包括:
地震数据读入模块,用于输入地震偏移速度场数据和地震炮记录数据;
逆时偏移成像处理模块,用于对输入地震数据进行逆时偏移成像处理,首先基于输入的偏移速度场进行震源波场正向延拓,对地震炮记录数据进行检波波场反向延拓,得到地震波在地下传播每一时刻的正向延拓波场值和反向延拓波场值。再利用最小二乘互相关成像条件得到逆时偏剖面;
波恩近似正演模块,用于将最小二乘反演思路与逆时偏移相结合,将常规逆时偏移结果当作最小二乘反演的第一次结果,对逆时偏移剖面进行波恩近似正演得到波恩近似正演的波场记录;
地震数据残差计算模块,用于地震数据最小二乘逆时偏移残差值的计算,使用波恩近似正演的波场记录与实际输入的地震波场相减,得到地震数据最小二乘逆时偏移计算的残差值;
共轭梯度计算模块,用于共轭梯度的计算,使用残差值,在常规互相关共轭梯度法最小二乘逆时偏移的基础上,采用共轭梯度归一化方法进行共轭梯度的计算;
更新步长计算及模型更新模块,用于使用共轭梯度值计算更新步长,并在原来偏移剖面的基础上进行模型的更新,完成归一化互相关成像条件最小二乘逆时偏移的一次迭代;
收敛输出模块,用于判断偏移结果是否满足输出偏移成像收敛条件,当偏移成像剖面结果满足收敛条件时,输出地震偏移剖面,如果偏移剖面不满足收敛条件时,对偏移剖面结果进行波恩近似正演并继续迭代,直到满足收敛条件时输出偏移剖面。
本发明的另一目的在于提供一种地球物理勘探终端,所述地球物理勘探终端搭载实施所述共轭梯度归一化LSRTM方法的系统。
本发明与现有技术相比的有益效果:
本发明应用于模型及实际地震资料,同时综合考虑了地质、地球物理等因素的影响,相较常规最小二乘法逆时偏移成像,中深部成像精度大幅提高,同时迭代收敛速度更快,可以大幅度降低偏移的计算成本。
共轭梯度归一化最小二乘法逆时偏移成像方法从一种波场传播理论出发,充分考虑了地球物理勘探日益复杂的特点,同时结合了常规最小二乘逆时偏移成像方法分辨率高的优点。随着勘探目标的日益复杂和勘探难度的不断增加,对于复杂构造的高精度成像和保幅成像的要求也不断提高,地震勘探的数据量也在不断增加。如何提高复杂构造成像精度,如何进行更加高效的偏移成像,节省地震资料处理的时间成本,为下一步的地震资料解释提供良好的成果剖面,成为了现在的研究热点和研究难点。目前,逆时偏移被认为是适应复杂构造成像最精确的成像方法,该方法的思想是利用波动方程沿时间方向对波场进行正向延拓和反向延拓并将地下每个成像点在每一时刻的波场值都存储起来,最后选取合适的成像条件进行成像得到偏移结果来反映地下的真实构造情况。然而逆时偏移还存在低频噪音,保幅较差,对复杂构造的成像效果较差等问题。最小二乘逆时偏移是将最小二乘反演思路与逆时偏移相结合的一种方法,与常规逆时偏移相比,它的成像结果具有更高的保真度、保幅性和分辨率,对中深层储层的成像效果更好。
最小二乘偏移思想最早是由Tarantola等人提出的,之后Dai W等实现了基于相位编码的最小二乘逆时偏移方法和平面波最小二乘逆时偏移技术的研究与应用。同样的,最小二乘逆时偏移也存在着浅层低频噪音,剖面能量分布不均的问题,同时其计算量也是巨大。本发明创新性地将震源归一化互相关思想引入最小二乘逆时偏移当中,从而有效减弱了震源效应的影响,减少了浅层低频噪音的干扰;同时,优化迭代过程中的梯度公式,使得波恩近似正演的记录与观测记录可以更快的逼近,从而提高了算法的计算效率。本发明的共轭梯度归一化最小二乘法逆时偏移成像方法的实现思路是:首先输入地震数据,进行常规逆时偏移得到偏移成像剖面,接下来对其进行波恩近似正演得到波恩近似正演记录并与输入地震记录进行残差的计算,之后利用残差进行共轭梯度归一化算子的计算并更新模型,最后根据收敛条件判断是否输出偏移剖面,完成共轭梯度归一化最小二乘法逆时偏移成像。
本发明可以适用于复杂构造的偏移成像,且偏移结果不受倾角的限制,可以对多种特殊波进行成像,如转换波、多次波、回转波等;对于速度模型适应性强,可以适应速度的横向变化,对于地层中的局部特殊构造也可以很好成像。模型和实际资料试算结果表明,共轭梯度归一化最小二乘法逆时偏移成像方法可以实现对于复杂构造的精细成像,从而为后续进行高精度地震解释提供了有力保障。
本发明的方法创新性地将震源归一化互相关思想引入最小二乘逆时偏移当中,研究出一种共轭梯度归一化最小二乘法逆时偏移成像方法,该方法可以有效压制浅层低频干扰,增加深层照明,提高剖面分辨率。逆时偏移中震源归一化互相关成像条件可以有效减弱震源效应的影响,减少浅层低频噪音的干扰,增加深层照明,对偏移剖面有不错的改善作用。常规共轭梯度法最小二乘逆时偏移中的梯度在没有合适预条件算子的作用下也会有同样的问题,而归一化思想能够解决逆时偏移中震源效应问题,同样也适用于共轭梯度法的求取。共轭梯度归一化最小二乘法逆时偏移成像方法采用归一化互相关思想,将迭代算法中的梯度求取公式进行优化。模型和实际资料试算结果表明,同常规最小二乘逆时偏移相比,共轭梯度归一化最小二乘法逆时偏移成像方法可以有效降低剖面浅层低频噪声,增强地下深部照明,中、深层成像精度得到大幅提高。
本发明的方法创新性地将震源归一化互相关思想引入最小二乘逆时偏移当中,研究出一种共轭梯度归一化最小二乘法逆时偏移成像方法,其可以大幅提高迭代收敛速度,从而降低了偏移成像的计算成本。共轭梯度归一化最小二乘法逆时偏移成像方法采用归一化互相关思想,将迭代算法中的梯度求取公式进行优化,优化后的迭代算法波恩近似正演的记录相较常规最小二乘逆时偏移与观测记录更接近,在相同迭代次数下其收敛速度更快,从而节省了大量的计算时间,对于实际生产应用推广具有重要的意义。
附图说明
图1是本发明实施例提供的共轭梯度归一化最小二乘法逆时偏移成像方法流程图。
图2是本发明实施例提供的共轭梯度归一化最小二乘法逆时偏移成像方法系统的结构示意图。
图3是本发明实施例提供的共轭梯度归一化最小二乘法逆时偏移成像方法实现流程图。
图4是本发明实施例提供的Marmousi2速度模型(a)及高斯平滑后速度场(b)示意图。
图5是本发明实施例提供的常规最小二乘逆时偏移迭代十五次处理结果(a)及共轭梯度归一化最小二乘法逆时偏移成像方法迭代十五次处理结果(b)示意图。
图6是本发明实施例提供的常规最小二乘逆时偏移剖面Laplace滤波结果(a)及共轭梯度归一化最小二乘法逆时偏移成像方法处理剖面Laplace滤波结果(b)示意图。
图7是本发明实施例提供的共轭梯度归一化最小二乘法逆时偏移成像方法与常规最小二乘逆时偏移误差下降曲线对比示意图图。
图8是本发明实施例提供的实际地震资料偏移速度场示意图。
图9是本发明实施例提供的实际资料单炮地震记录示意图。
图10是本发明实施例提供的实际资料常规归一化互相关逆时偏移结果(a)及共轭梯度归一化最小二乘法逆时偏移成像方法结果(b)示意图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白,以下结合实施例和附图,对本发明进行进一步详细说明。应当理解,此处所描述的具体实施例仅仅用以解释本发明,并不用于限定本发明。
针对现有技术存在的问题,本发明提供了一种共轭梯度归一化最小二乘法逆时偏移成像方法、系统、存储介质及应用,下面结合附图对本发明作详细的描述。
如图1所示,本发明提供的共轭梯度归一化最小二乘法逆时偏移成像方法包括以下步骤:
S101:完成地震偏移速度场数据和地震炮记录数据读入,资料试算中先使用Marmousi2模型进行发明方法正确性和有效性的验证,再用某地区实际资料试算进行发明方法适应性的验证。
S102:对输入地震数据进行逆时偏移成像处理:首先基于输入的偏移速度场进行震源波场正向延拓,对地震炮记录数据进行检波波场反向延拓,得到地震波在地下传播每一时刻的正向延拓波场值和反向延拓波场值。再利用最小二乘互相关成像条件得到逆时偏移剖面。
S103:进行波恩近似正演。将最小二乘反演思路与逆时偏移相结合,将上一步逆时偏移结果当作最小二乘反演的第一次结果,对逆时偏移剖面进行波恩近似正演得到波恩近似正演的波场记录。
S104:进行地震数据残差的计算。使用波恩近似正演的波场记录与实际输入的地震波场相减,得到地震数据最小二乘逆时偏移计算的残差值。
S105:进行共轭梯度的计算。使用残差值,在常规互相关共轭梯度法最小二乘逆时偏移的基础上,采用共轭梯度归一化方法进行共轭梯度的计算。
S106:进行更新步长的计算并完成模型的更新。使用共轭梯度值计算更新步长,并在原来偏移剖面的基础上进行模型的更新,完成归一化互相关成像条件最小二乘逆时偏移的一次迭代。
S107:判断偏移结果是否满足输出偏移成像收敛条件。当偏移成像剖面结果满足收敛条件时,输出地震偏移剖面,如果偏移剖面不满足收敛条件时,对偏移剖面结果进行波恩近似正演并继续迭代,直到满足收敛条件时输出偏移剖面。在地震资料计算中,一般设定收敛条件为迭代次数,当满足迭代次数后进行偏移成像剖面的输出。
本发明提供的共轭梯度归一化最小二乘法逆时偏移成像方法,业内的普通技术人员还可以采用其他的步骤实施,图1的本发明提供的共轭梯度归一化最小二乘法逆时偏移成像方法仅仅是一个具体实施例而已。
如图2所示,本发明提供的共轭梯度归一化最小二乘法逆时偏移成像方法系统包括:
地震数据读入模块1,用于输入地震偏移速度场数据和地震炮记录数据,资料试算中先使用Marmousi2模型进行方法正确性和有效性的验证,再用某地区实际资料试算进行方法适应性的验证。
逆时偏移成像处理模块2。用于对输入地震数据进行逆时偏移成像处理,首先基于输入的偏移速度场进行震源波场正向延拓,对地震炮记录数据进行检波波场反向延拓,得到地震波在地下传播每一时刻的正向延拓波场值和反向延拓波场值。再利用最小二乘互相关成像条件得到逆时偏移剖面。
波恩近似正演模块3。用于将最小二乘反演思路与逆时偏移相结合,将常规逆时偏移结果当作最小二乘反演的第一次结果,对逆时偏移剖面进行波恩近似正演得到波恩近似正演的波场记录。
地震数据残差计算模块4。用于地震数据最小二乘逆时偏移残差值的计算,使用波恩近似正演的波场记录与实际输入的地震波场相减,得到地震数据最小二乘逆时偏移计算的残差值。
共轭梯度计算模块5。用于共轭梯度的计算,使用残差值,在常规互相关共轭梯度法最小二乘逆时偏移的基础上,采用共轭梯度归一化方法进行共轭梯度的计算。
更新步长计算及模型更新模块6。用于使用共轭梯度值计算更新步长,并在原来偏移剖面的基础上进行模型的更新,完成归一化互相关成像条件最小二乘逆时偏移的一次迭代。
判断输出模块7。用于判断是否满足输出偏移成像收敛条件。当偏移成像剖面结果满足收敛条件时,输出地震偏移剖面,如果偏移剖面不满足收敛条件时,对偏移剖面结果进行波恩近似正演并继续迭代,直到满足收敛条件时输出偏移剖面。在地震资料计算中,一般设定收敛条件为迭代次数,当满足迭代次数后进行偏移成像剖面的输出。
下面结合附图对本发明的技术方案作进一步的描述。
如图3所示,本发明提供的共轭梯度归一化最小二乘法逆时偏移成像方法包括以下步骤:
步骤100,读入地震数据。读入地震数据的偏移速度场和地震炮记录,以Marmousi2模型为例,模型水平方向长3250m,纵向深度2425m,炮间距250m,均匀分布地表,共计正演14炮;道间距5米,650道接收;接收道不移动,炮点等间隔移动。其偏移真实速度场如图4(a)所示,为了更加接近地震勘探真实情况对其进行高斯平滑,高斯平滑后的速度场如图4(b)所示。
步骤101,对输入地震偏移速度场进行逆时偏移成像处理。共轭梯度归一化最小二乘法逆时偏移成像第一次偏移的做法和常规逆时偏移一致,不同的是,最小二乘逆时偏移改变了互相关中波场的表达式。常规逆时偏移的表达式为:
m0(x,z)=LTD
式中m0代表的是逆时偏移剖面,LT代表的是近似偏移算子,D代表采集的地震数据。具体做法为,首先基于输入的偏移速度场进行震源波场正向延拓,对地震炮记录进行检波波场反向延拓,得到地震波在地下传播每一时刻的正向延拓波场值和反向延拓波场值。再利用最小二乘互相关成像条件得到逆时偏移剖面。后续将这一步生成的常规逆时偏移结果当作共轭梯度归一化最小二乘反演的第一次结果。
步骤102,利用逆时偏移结果进行波恩近似正演,一般使用以下公式:
d=Lm
其中m为偏移剖面或反射系数模型的矩阵形式,d为模拟数据的矩阵形式,L为波恩近似正传算子矩阵。
步骤103,进行地震数据残差的计算。使用波恩近似正演的波场记录与实际输入的地震波场相减,得到地震数据最小二乘逆时偏移计算的残差值。
步骤104,进行共轭梯度的计算。共轭梯度归一化最小二乘法逆时偏移共轭梯度算法为:
步骤105,进行更新步长的计算并完成模型的更新。迭代算法如下:
步骤106,对迭代更新结果进行判断,看其是否满足输出偏移成像收敛条件。当偏移成像剖面结果满足收敛条件时,输出地震偏移剖面,如果偏移剖面不满足收敛条件时,对偏移剖面结果进行波恩近似正演并继续迭代,直到满足收敛条件时输出偏移剖面。在地震资料计算中,一般设定收敛条件为迭代次数,当满足迭代次数后进行偏移成像剖面的输出。
Marmousi2模型常规最小二乘逆时偏移迭代十五次处理结果如图5(a)所示,共轭梯度归一化最小二乘法逆时偏移成像方法迭代十五次处理结果如图5(b)所示。为了更好地对比二者的效果。分别对其进行Laplace滤波处理得到图6(a)和图6(b)。可以看到无论是滤波前还是滤波后在相同迭代次数的条件下,共轭梯度归一化最小二乘法逆时偏移成像方法对于浅部低频的压制更加显著,成像效果更好。而对于构造复杂的中深部地层可以看出共轭梯度归一化最小二乘法逆时偏移成像方法对于深部的照明同样是显著强于常规最小二乘逆时偏移处理的效果。不管是在浅部的噪音压制上还是在深部剖面的照明上,共轭梯度归一化最小二乘法逆时偏移成像方法都强于常规最小二乘逆时偏移方法。
为了准确对比共轭梯度归一化最小二乘法逆时偏移成像方法的收敛速度与常规最小二乘逆时偏移的收敛速度,图7给出了是观测记录与模拟记录归一化误差下降曲线,横坐标是两种最小二乘逆时偏移中迭代次数,纵坐标是相对误差,其中相对误差=|(观测记录-模拟记录)/观测记录|,图中红色线条表示常规最小二乘法逆时偏移成像方法,蓝线表示共轭梯度归一化最小二乘法逆时偏移成像方法。从误差下降曲线可以看出,在迭代相同次数的情况下,共轭梯度归一化最小二乘法逆时偏移成像方法波恩近似正演的记录与观测记录更接近,收敛速度更快,而收敛速度更快意味着生产效率更高,对于实际生产具有重要的意义。
实际地震资料也显示出同样的应用效果。实际地震资料偏移速度场以及单炮记录如图8和图9所示,二维数据共328炮,炮点的初始位置位于最左端,炮点在检波点排列的中间且相对位置固定,检波点排列随着炮点的移动而移动。由于采用这种多次覆盖的观测系统,因此每炮的道数是不固定的。从左边开始第一炮162道接收,然后依次增加2道,一直增加到320道,然后中间有167炮满覆盖,再依次减少2道,到最后一炮160道接收。炮间距是100m,是道间距的2倍,CMP间距的4倍。通过估算取地震子波的主频为30Hz,时间采样点数为1501个,时间采样率为4ms,地震记录的总时长为6s。其归一化互相关逆时偏移结果和共轭梯度归一化最小二乘法逆时偏移成像方法结果如图10(a)和图10(b)所示。可以看到共轭梯度归一化最小二乘法逆时偏移成像方法结果。无论是在浅层低频噪音的压制方面,还是在深层构造的照明方面,处理的都很好。浅层和深层的能量较为均衡,成像较清晰,同相轴的连续性较好,信噪比和分辨率较高。
在于步骤104中,常规最小二乘法逆时偏移成像最速下降梯度g公式为:
g(k+1)=LT[Lm(k)-D]
梯度g的具体求法展开写可表示为:
式中SS是震源正传波场,Rres是波恩近似正演数据与采集数据残差的反传波场。考虑到常规共轭梯度法最小二乘逆时偏移中预条件算子难以求取,且不能很好地近似Hessian矩阵,造成震源效应对梯度的影响十分大,与逆时偏移中出现的问题相似,所以采用了归一化思想,求取震源归一化共轭梯度法来减弱震源效应的影响,提高收敛速度。改进后共轭梯度归一化最小二乘法逆时偏移最速下降梯度为:
应当注意,本发明的实施方式可以通过硬件、软件或者软件和硬件的结合来实现。硬件部分可以利用专用逻辑来实现;软件部分可以存储在存储器中,由适当的指令执行系统,例如微处理器或者专用设计硬件来执行。本领域的普通技术人员可以理解上述的设备和方法可以使用计算机可执行指令和/或包含在处理器控制代码中来实现,例如在诸如磁盘、CD或DVD-ROM的载体介质、诸如只读存储器(固件)的可编程的存储器或者诸如光学或电子信号载体的数据载体上提供了这样的代码。本发明的设备及其模块可以由诸如超大规模集成电路或门阵列、诸如逻辑芯片、晶体管等的半导体、或者诸如现场可编程门阵列、可编程逻辑设备等的可编程硬件设备的硬件电路实现,也可以用由各种类型的处理器执行的软件实现,也可以由上述硬件电路和软件的结合例如固件来实现。
以上所述,仅为本发明的具体实施方式,但本发明的保护范围并不局限于此,任何熟悉本技术领域的技术人员在本发明揭露的技术范围内,凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等,都应涵盖在本发明的保护范围之内。
Claims (8)
1.一种共轭梯度归一化LSRTM方法,其特征在于所述方法的具体步骤如下:
完成地震偏移速度场数据和地震炮记录数据输入,对输入地震数据进行逆时偏移成像处理:首先基于输入的偏移速度场进行震源波场正向延拓,对地震炮记录数据进行检波波场反向延拓,得到地震波在地下传播每一时刻的正向延拓波场值和反向延拓波场值,再利用最小二乘逆时偏移互相关成像条件得到逆时偏移剖面;
进行波恩近似正演:将最小二乘反演思路与逆时偏移相结合,将常规逆时偏移结果当作最小二乘反演的第一次结果,对逆时偏移剖面进行波恩近似正演得到波恩近似正演的波场记录;
进行地震数据残差的计算:使用波恩近似正演的波场记录与实际输入的地震波场相减,得到地震数据最小二乘逆时偏移计算的残差值;
进行共轭梯度的计算:使用残差值,在常规互相关共轭梯度法最小二乘逆时偏移的基础上,采用共轭梯度归一化方法进行共轭梯度的计算;
进行更新步长的计算并完成模型的更新:使用共轭梯度值计算更新步长,并在原来偏移剖面的基础上进行模型的更新,完成归一化互相关成像条件最小二乘逆时偏移的一次迭代;
判断偏移结果是否满足输出偏移成像收敛条件:当偏移成像剖面结果满足收敛条件时,输出地震偏移剖面,如果偏移剖面不满足收敛条件时,对偏移剖面结果进行波恩近似正演并继续迭代,直到满足收敛条件时输出偏移剖面。
2.根据权利要求1所述的一种共轭梯度归一化LSRTM方法,其特征在于所述的进行波恩近似正演,使用以下公式:
d=Lm
其中m为偏移剖面或反射系数模型的矩阵形式,d为模拟数据的矩阵形式,L为波恩近似正传算子矩阵。
5.根据权利要求1所述的一种共轭梯度归一化LSRTM方法,其特征在于所述判断偏移结果是否满足输出偏移成像收敛条件:在地震资料计算中,设定收敛条件为迭代次数,当满足迭代次数后进行偏移成像剖面的输出。
6.一种计算机可读存储介质,存储有计算机程序,其特征在于所述计算机程序被处理器执行时,所述处理器执行权利要求1所述的一种共轭梯度归一化LSRTM方法。
7.一种实施权利要求1所述共轭梯度归一化LSRTM方法的系统,其特征在于所述系统包括:
地震数据读入模块,用于输入地震偏移速度场数据和地震炮记录数据;
逆时偏移成像处理模块,用于对输入地震数据进行逆时偏移成像处理,首先基于输入的偏移速度场进行震源波场正向延拓,对地震炮记录数据进行检波波场反向延拓,得到地震波在地下传播每一时刻的正向延拓波场值和反向延拓波场值,再利用最小二乘互相关成像条件得到逆时偏剖面;
波恩近似正演模块,用于将最小二乘反演思路与逆时偏移相结合,将常规逆时偏移结果当作最小二乘反演的第一次结果,对逆时偏移剖面进行波恩近似正演得到波恩近似正演的波场记录;
地震数据残差计算模块,用于地震数据最小二乘逆时偏移残差值的计算,使用波恩近似正演的波场记录与实际输入的地震波场相减,得到地震数据最小二乘逆时偏移计算的残差值;
共轭梯度计算模块,用于共轭梯度的计算,使用残差值,在常规互相关共轭梯度法最小二乘逆时偏移的基础上,采用共轭梯度归一化方法进行共轭梯度的计算;
更新步长计算及模型更新模块,用于使用共轭梯度值计算更新步长,并在原来偏移剖面的基础上进行模型的更新,完成归一化互相关成像条件最小二乘逆时偏移的一次迭代;
收敛输出模块,用于判断偏移结果是否满足输出偏移成像收敛条件,当偏移成像剖面结果满足收敛条件时,输出地震偏移剖面,如果偏移剖面不满足收敛条件时,对偏移剖面结果进行波恩近似正演并继续迭代,直到满足收敛条件时输出偏移剖面。
8.一种地球物理勘探终端,其特征在于所述地球物理勘探终端搭载权利要求7所述的系统。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110323886.3A CN113064203B (zh) | 2021-03-26 | 2021-03-26 | 共轭梯度归一化lsrtm方法、系统、存储介质及应用 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110323886.3A CN113064203B (zh) | 2021-03-26 | 2021-03-26 | 共轭梯度归一化lsrtm方法、系统、存储介质及应用 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN113064203A true CN113064203A (zh) | 2021-07-02 |
CN113064203B CN113064203B (zh) | 2022-04-01 |
Family
ID=76563713
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202110323886.3A Active CN113064203B (zh) | 2021-03-26 | 2021-03-26 | 共轭梯度归一化lsrtm方法、系统、存储介质及应用 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN113064203B (zh) |
Cited By (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113705644A (zh) * | 2021-08-17 | 2021-11-26 | 西安交通大学 | 一种物理规律和数据双驱动的地震成像方法、系统、设备及存储介质 |
CN115166827A (zh) * | 2022-07-15 | 2022-10-11 | 中山大学 | 基于反褶积成像条件的最小二乘偏移成像方法、设备及存储介质 |
CN115951402A (zh) * | 2022-12-29 | 2023-04-11 | 中山大学 | 基于向量反射率正演算子的偏移成像方法、系统及介质 |
CN115951401A (zh) * | 2022-07-19 | 2023-04-11 | 中山大学 | 成像条件驱动的最小二乘逆时偏移成像方法、设备及存储介质 |
Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104360381A (zh) * | 2014-10-20 | 2015-02-18 | 李闯 | 一种地震资料的偏移成像处理方法 |
CN107783190A (zh) * | 2017-10-18 | 2018-03-09 | 中国石油大学(北京) | 一种最小二乘逆时偏移梯度更新方法 |
CN109407152A (zh) * | 2018-12-18 | 2019-03-01 | 吉林大学 | 基于零均值归一化互相关目标函数的时间域全波形反演方法 |
CN110058302A (zh) * | 2019-05-05 | 2019-07-26 | 四川省地质工程勘察院 | 一种基于预条件共轭梯度加速算法的全波形反演方法 |
CN110737018A (zh) * | 2019-07-09 | 2020-01-31 | 中国石油化工股份有限公司 | Vsp地震资料各向异性建模方法 |
US20210055443A1 (en) * | 2019-08-19 | 2021-02-25 | Cgg Services Sas | Prestack least-square reverse time migration on surface attribute gathers compressed using depth-independent coefficients |
-
2021
- 2021-03-26 CN CN202110323886.3A patent/CN113064203B/zh active Active
Patent Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104360381A (zh) * | 2014-10-20 | 2015-02-18 | 李闯 | 一种地震资料的偏移成像处理方法 |
CN107783190A (zh) * | 2017-10-18 | 2018-03-09 | 中国石油大学(北京) | 一种最小二乘逆时偏移梯度更新方法 |
CN109407152A (zh) * | 2018-12-18 | 2019-03-01 | 吉林大学 | 基于零均值归一化互相关目标函数的时间域全波形反演方法 |
CN110058302A (zh) * | 2019-05-05 | 2019-07-26 | 四川省地质工程勘察院 | 一种基于预条件共轭梯度加速算法的全波形反演方法 |
CN110737018A (zh) * | 2019-07-09 | 2020-01-31 | 中国石油化工股份有限公司 | Vsp地震资料各向异性建模方法 |
US20210055443A1 (en) * | 2019-08-19 | 2021-02-25 | Cgg Services Sas | Prestack least-square reverse time migration on surface attribute gathers compressed using depth-independent coefficients |
Non-Patent Citations (1)
Title |
---|
李庆洋 等: "去均值归一化互相关最小二乘逆时偏移及其应用", 《地球物理学报》 * |
Cited By (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN113705644A (zh) * | 2021-08-17 | 2021-11-26 | 西安交通大学 | 一种物理规律和数据双驱动的地震成像方法、系统、设备及存储介质 |
CN113705644B (zh) * | 2021-08-17 | 2023-09-12 | 西安交通大学 | 一种物理规律和数据双驱动的地震成像方法、系统、设备及存储介质 |
CN115166827A (zh) * | 2022-07-15 | 2022-10-11 | 中山大学 | 基于反褶积成像条件的最小二乘偏移成像方法、设备及存储介质 |
CN115166827B (zh) * | 2022-07-15 | 2023-04-28 | 中山大学 | 基于反褶积成像条件的最小二乘偏移成像方法、设备及存储介质 |
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 |
---|---|
CN113064203B (zh) | 2022-04-01 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN113064203B (zh) | 共轭梯度归一化lsrtm方法、系统、存储介质及应用 | |
RU2693495C1 (ru) | Полная инверсия волнового поля с компенсацией показателя качества | |
RU2694621C1 (ru) | Способ и устройство для обработки сейсмических данных | |
CN109669212B (zh) | 地震数据处理方法、地层品质因子估算方法与装置 | |
US10234581B2 (en) | System and method for high resolution seismic imaging | |
CN109946741B (zh) | 一种TTI介质中纯qP波最小二乘逆时偏移成像方法 | |
CN110531410B (zh) | 一种基于直达波场的最小二乘逆时偏移梯度预条件方法 | |
US11733413B2 (en) | Method and system for super resolution least-squares reverse time migration | |
CN111290019B (zh) | 一种应用于最小二乘逆时偏移的l-bfgs初始矩阵求取方法 | |
CN114460646B (zh) | 一种基于波场激发近似的反射波旅行时反演方法 | |
CN114966837A (zh) | 基于卷积型w2距离目标函数的全波形反演方法及装置 | |
CN106353798A (zh) | 多分量联合高斯束叠前逆时偏移成像方法 | |
CN113156500B (zh) | 数据驱动的快速构造约束叠前地震多道反演方法 | |
US20220283329A1 (en) | Method and system for faster seismic imaging using machine learning | |
CN110780341B (zh) | 一种各向异性地震成像方法 | |
CN116719086B (zh) | 基于点扩散函数的稀疏海底四分量数据高分辨率成像方法 | |
CN112630830B (zh) | 一种基于高斯加权的反射波全波形反演方法及系统 | |
CN112630825B (zh) | 共炮检距域Beam叠前时间偏移成像方法、系统、介质及应用 | |
CN111538083A (zh) | 基于速度梯度的崎岖海底界面的光滑处理方法 | |
Zou et al. | Log-constrained inversion based on a conjugate gradient-particle swarm optimization algorithm | |
US12000971B2 (en) | Method and system for seismic processing using virtual trace bins based on offset attributes and azimuthal attributes | |
US20230184972A1 (en) | Method and system for seismic processing using virtual trace bins based on offset attributes and azimuthal attributes | |
CN112014875B (zh) | 叠前地震反演方法及装置 | |
CN112462428B (zh) | 一种多分量地震资料偏移成像方法及系统 | |
CN116359982A (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 |