CN111025388B - 一种多波联合的叠前波形反演方法 - Google Patents
一种多波联合的叠前波形反演方法 Download PDFInfo
- Publication number
- CN111025388B CN111025388B CN201911315597.8A CN201911315597A CN111025388B CN 111025388 B CN111025388 B CN 111025388B CN 201911315597 A CN201911315597 A CN 201911315597A CN 111025388 B CN111025388 B CN 111025388B
- Authority
- CN
- China
- Prior art keywords
- wave
- matrix
- inversion
- model
- data
- 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.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 58
- 238000004364 calculation method Methods 0.000 claims abstract description 12
- 230000004044 response Effects 0.000 claims abstract description 7
- 238000012545 processing Methods 0.000 claims abstract description 3
- 239000011159 matrix material Substances 0.000 claims description 84
- 238000012546 transfer Methods 0.000 claims description 30
- 230000005540 biological transmission Effects 0.000 claims description 14
- 239000010410 layer Substances 0.000 claims description 12
- 238000004422 calculation algorithm Methods 0.000 claims description 11
- 230000015572 biosynthetic process Effects 0.000 claims description 10
- 238000000354 decomposition reaction Methods 0.000 claims description 9
- 238000003786 synthesis reaction Methods 0.000 claims description 8
- 239000000126 substance Substances 0.000 claims description 6
- 238000004088 simulation Methods 0.000 claims description 5
- 238000012937 correction Methods 0.000 claims description 4
- 239000011229 interlayer Substances 0.000 claims description 4
- 238000011160 research Methods 0.000 claims description 4
- ZJPGOXWRFNKIQL-JYJNAYRXSA-N Phe-Pro-Pro Chemical compound C([C@H](N)C(=O)N1[C@@H](CCC1)C(=O)N1[C@@H](CCC1)C(O)=O)C1=CC=CC=C1 ZJPGOXWRFNKIQL-JYJNAYRXSA-N 0.000 claims description 3
- 238000009795 derivation Methods 0.000 claims description 3
- 238000001914 filtration Methods 0.000 claims description 3
- 230000003044 adaptive effect Effects 0.000 claims description 2
- 238000009792 diffusion process Methods 0.000 claims description 2
- 230000003068 static effect Effects 0.000 claims description 2
- 230000007704 transition Effects 0.000 claims description 2
- 238000012360 testing method Methods 0.000 abstract description 5
- 230000035945 sensitivity Effects 0.000 abstract 1
- 230000002194 synthesizing effect Effects 0.000 description 3
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000005457 optimization Methods 0.000 description 2
- 238000002310 reflectometry Methods 0.000 description 2
- 238000005070 sampling Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 1
- 230000002068 genetic effect Effects 0.000 description 1
- 230000008569 process Effects 0.000 description 1
- 238000002922 simulated annealing Methods 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
- 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
- G01V1/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
-
- 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
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
本发明公开了一种多波联合的叠前波形反演方法。叠前地震反演方法主要包括基于射线追踪类正演算子和基于波动方程类正演算子两大类。前者即为应用最为广泛的叠前AVO反演,后者则被称为叠前波形反演。不同于叠前AVO反演,叠前波形反演方法可以更好应用输入数据中的复杂波场响应,减少输入数据的处理难度,提高地震记录对目标反演参数的敏感度。目前,叠前波形反演多采用非线性反演策略,计算效率过低。本发明提出了一种线性策略的叠前反演方法,以降低计算时间。同时为提高反演精度和稳定性,本发明采用了PP和PS多波数据联合反演策略,该方法可有效提高三参数的反演精度,尤其是针对横波速度和密度估算。测井模型测试验证了方法的有效性。
Description
技术领域
本发明涉及地震勘探技术领域,属于地震资料多参数反演范畴,尤其涉及一种多波联合的叠前波形反演方法。
背景技术
目前,叠前地震反演方法主要包括两类,一是基于射线追踪类正演算子的叠前AVO反演方法,二是基于波动方程类正演算子的叠前波形反演[1]。第一类叠前AVO反演优势在于单一界面简单假设,数学形式相对简单,具有较高的计算效率,因此在实际地震勘探中应用最为广泛。不同于叠前AVO反演,叠前波形反演采用了波动方程类模拟算法作为正演算子,能更好利用实际数据中的全波场响应(如多次波和转换波等)开展地层介质参数反演估算。
波动方程类正演包括有反射率[2-5]和传递矩阵[6,7]等方法。反射率法目前被用于叠前波形反演的正演算子,并取得了理想的估算结果[8,10-11]。绝大多数现存的叠前波形反演方法采用非线性的反演方法,如选用模拟退火[9]、遗传算法等优化策略[10-13]。模型测试、2D和3D实际数据应用均表明了基于非线性优化算法的叠前波形反演在介质属性估算方面的优势。非线性反演方法不依赖于初始模型,无需求取偏导数实现简单,同时有助于获取全局最优解,避免落入局部极值。然而,该类算法需要多次重复正演过程,波动方程类正演计算复杂,因此采用非线性策略的叠前波形反演计算量过大,计算成本过高。因此,为了解决这一问题,需开展基于线性反演策略的叠前波形反演研究,降低计算成本,更易于向实际数据应用中推广。
针对各向同性介质的叠前地震反演,通常反演纵波速度、横波速度和密度三参数。输入数据不足会导致多参数反演的不稳定。为了提高反演的稳定性和精度,多波数据联合反演能有效提高参数估算精度[11,13-15]。密度是开展岩性识别和油气勘探的重要参数。Luo等指出,单纯的PP数据反演无法准确地获得横波速度和密度结果[16]。即使在较大入射角情况,由于密度对PP反射振幅不够敏感,密度估算的不精确性仍然是一个问题。转换波(PS数据)包含丰富的横波速度和密度信息。因此,PP和PS联合反演可以提供更好的反演结果。因此,为了提高叠前波形反演的横波速度和密度的估算精度,需要开展多波联合的叠前波形反演方法研究。
参考文献
[1]Simmons,J.L.,and M.M.Backus.1994,AVO modeling and the locallyconverted shear wave.Geophysics,59,no.8,1237-1248.
[2]Fryer,G.J.,and L.N.Frazer.1984,Seismic waves in stratifiedanisotropic media.Geophysical Journal of the Royal Astronomical Society,78,no.3,691-710.
[3]Fuchs,K.,and G.Müller.1971,Computation of synthetic seismogramswith the reflectivity method and comparison with observations.GeophysicalJournal International,23,no.4,417-433.
[4]Kennett,B.L.N.1983,Seismic Wave Propagation In Stratified Media:Cambridge University Press.
[5]Booth,D.C.,and S.Crampin.1983,The anisotropic reflectivitytechnique:theory.Geophysical Journal International,72,no.3,755-766.
[6]Carcione,J.M.2001a,Amplitude variations with offset of pressure-seal reflections.Geophysics,66,no.1,283-293.
[7]Carcione,J.M.2001b,Wave fields in real media:Wave propagation inanisotropic,anelastic,porous and electromagnetic media.Vol.38:Elsevier.
[8]Sen,M.K.,and P.L.Stoffa.1991,Nonlinear seismic waveform inversionin one dimension using simulated annealing.Geophysics,56,no.10,1624-1638.
[9]Pafeng,J.,S.Mallick,and H.Sharma.2016,Prestack waveform inversionof three-dimensional seismic data—An example from the Rock Springs Uplift,Wyoming,USA.Geophysics,82,no.1,B1-B12.
[10]Padhi,A.,and S.Mallick.2013,Accurate estimation of density fromthe inversion of multicomponent prestack seismic waveform data using anondominated sorting genetic algorithm.Geophysics,32,no.1,94-98.
[11]Padhi,A.,and S.Mallick.2014,Multicomponent pre-stack seismicwaveform inversion in transversely isotropic media using a non-dominatedsorting genetic algorithm.Geophysical Journal International,196(3):1600-1618.
[12]Gisolf,A.,P.Haffinger,C.Hanitzsch,P.Doulgeris,and P.Veeken.2014,Non-linear full wavefield inversion applied to carboniferous reservoirs inthe Wingate gas field(SNS,Offshore UK).76th EAGE Conference and Exhibition2014.
[13]Li,T.,and S.Mallick.2015,Multicomponent,multi-azimuth pre-stackseismic waveform inversion for azimuthally anisotropic media using a paralleland computationally efficient non-dominated sorting geneticalgorithm.Geophysical Journal International,200,no.2,1136-1154.
[14]Lu,J.,Y.Wang,J.Chen,and Y.An.2018,Joint anisotropic amplitudevariation with offset inversion of PP and PS seismic data.Geophysics,83,no.2,N31-N50..
[15]Lu,J.,Z.Yang,Y.Wang,and Y.Shi.2015a,Joint PP and PS AVA seismicinversion using exact Zoeppritz equations.Geophysics,80,no.5,R239-R250.
[16]Luo,C.,X.Y.Li,and G.Huang.2018b,Information on S-wave velocityand density obtained from multi-wave data.SEG Technical Program ExpandedAbstracts.
发明内容
发明目的:为克服以往叠前波形非线性反演方法的计算效率过低问题,本发明提出一种多波联合的叠前波形反演方法,选用了线性反演策略,利用全波场响应多波记录的叠前波形线性反演方法,更多利用多波信息和地震记录的转换波、多次波信息,提高弹性参数反演精度。
技术方案:为实现本发明的目的,本发明所采用的技术方案是:一种多波联合的叠前波形反演方法,包括以下步骤:
步骤一:获取研究工区的叠前地震数据,即包含全波场响应的PP和PS道集数据;正演方法的合成记录应准确模拟实际全波场数据,两者一致性好有利于获得真实的反演结果;广义传递矩阵算法可以精确模拟层状模型内部复杂的波场响应,因此叠前PP和PS道集中应包括PP反射波、PS反射波、层间多次波、层间转换波;优选的,叠前道集进行球面扩散补偿、静校正、动校正和表层多次波去除处理;
步骤二:从叠前道集中分别提取PP波和PS波的地震子波,建立子波矩阵;
步骤三:利用搜集到的测井数据,建立目标弹性参数深度域初始模型;
步骤四:基于深度域初始模型,利用广义传递矩阵正演算子获取PP波和PS波合成记录;
步骤五:将非线性正演算子线性化,建立目标函数,计算PP波和PS波合成记录与实际叠前地震数据的残差;
步骤六:利用初始模型参数计算广义传递矩阵正演算子中PP波和PS波合成记录对模型参数的偏导数;
步骤七:采用自适应估算方法获取最优的正则化权重参数:
步骤八:计算目标函数的一阶导数矩阵,并通过迭代方法获取伪海森矩阵;
步骤九:计算模型的更新方向,更新模型参数;
步骤十:对更新后的模型进行判别,若模型误差未降到预设范围,进行下一次反演迭代,重复步骤四到九;否则迭代停止,输出参数的反演结果。
进一步的,所述步骤二,从叠前道集中分别提取PP波和PS波的地震子波,建立子波矩阵,具体包括:
从叠前PP波角道集、PS波角道集中提取不同角度情况的地震子波序列,分别建立PP记录、PS记录的子波矩阵WPP(θ)、WPS(θ);表达式如下:
其中,θ表示角度,N表示角道集中角度数量;
建立PP和PS联合记录的子波矩阵,用于制作PP和PS联合合成记录,如下:
进一步的,所述步骤三,利用搜集到的测井数据,建立目标弹性参数深度域初始模型;方法如下:
不同于基于Zoeppritz方程的叠前AVO/AVA反演,基于广义传递矩阵的叠前波形反演需构建深度域初始模型。建立深度域初始模型具体步骤包括:
无需对测井数据进行时间-深度匹配,直接对深度域井口数据进行低通滤波,去除测井数据的中高频信息,仅保留低频趋势;获取时间域层位数据解释结果,对时间域结果进行时间-深度匹配,获得深度域层位数据;从井口低频趋势出发,对其沿层位数据进行插值处理,获得深度域2D/3D初始模型。
进一步的,所述步骤四,基于深度域初始模型,利用广义传递矩阵正演算子获取PP波和PS波合成记录;
初始模型参数包括纵波速度α、横波速度β、密度ρ,模型向量m表达式如下:
m=[α1,…αM,β1,…βM,ρ1,…ρM] (4)
其中,每种参数共有M个;
利用初始模型参数,选用广义传递矩阵方法获取频率域的PP和PS反射系数:
式中,R为频率慢度域反射透射系数向量,A1和A2为层状模型顶层半空间和底层半空间介质的属性矩阵,n表示除去顶层半空间和底层半空间介质,中间部分包含有n层介质,Tj为中间第j和第j+1两相邻地层之间传递矩阵,b为一维列向量,式(5)具体形式如下:
其中,RPP、RPS、TPP和TPS分别代表PP波反射系数、PS波反射系数、PP波透射系数和PS波透射系数;
A1和A2具体形式如下:
式中,上标(*)1和(*)2分别表示顶层和底层的参数变量;下标(*)P和(*)S则表示纵波和横波的参数变量;i为虚数单位,ω为角频率,h为地层厚度;
Tj表达式如下:
式中,hj表示第j层介质的厚度;c11、c13和c55均为弹性刚度参数,是纵波速度、横波速度和密度的函数,px和pz分别为水平和垂直慢度;
利用式(5)获得频率域反射系数矩阵RPP和RPS,利用傅里叶变换将频率域反射系数转为时间域反射系数序列:
基于褶积理论,通过反射系数序列与地震子波矩阵(见式(1)和式(2))可获取合成地震记录:
上标(*)syn表示为合成记录数据,dPP表示PP波记录数据,dPS表示PS波记录数据。
进一步的,所述步骤五,将非线性正演算子线性化,建立目标函数,计算PP波和PS波合成记录与实际叠前地震数据的残差;
广义传递矩阵为非线性正演算法,其PP和PS模拟记录为模型参数的函数:
dPP=FPP(m),dPS=FPS(m) (18)
通过泰勒级数展开并去掉高阶及误差项,则将非线性正演驱动线性化,得到如下表达式:
上标(*)syn和(*)obs表示为合成记录数据和实际叠前地震道集数据;
建立反演目标函数,具体形式如下:
其中,C△m为模型参数m的协方差矩阵,u为模型参数m的期望,α为PP和PS数据的权重参数,λ为正则化参数,PP和PS记录的misfit项,表示如下:
其中,misfitPP、misfitPS分别表示PP和PS记录的misfit项。
进一步的,所述步骤六,计算广义传递矩阵算法PP和PS记录对模型参数的偏导数,利用初始模型参数计算其初始偏导数结果;
计算广义传递矩阵的偏导数,即非线性正演算子对模型参数的一阶导数,具体形式如下:
反射系数对模型参数的偏导数,如下:
进一步的,所述步骤七,利用L曲线自适应估算策略获取最优的正则化权重参数:
自适应正则化权重参数获取包括:
其中,U、Σ和V为偏导数矩阵的奇异值分解结果,u、δ、v为矩阵参数;N为总的数值个数;给定不同正则化权重参数计算残差,绘制L曲线图,选取拐点位置的权重参数作为最优权重;残差方程具体表达如下:
其中,△mλ用奇异值分解结果表示的形式为:
进一步的,所述步骤八,计算目标函数的一阶导数矩阵,并通过迭代方法获取伪海森矩阵,具体包括:
目标函数的具体形式见式(21),目标函数对模型参数的一阶导数为:
每一次反演迭代中都需要重新求取伪海森矩阵H。本发明采用迭代计算方法更新获得伪海森矩阵,设定最大迭代次数,当达到最大迭代次数,则停止迭代,输出伪海森矩阵。首次迭代计算时,伪海森矩阵具体形式为:
Vi=I-ρi(Ji-Ji-1)·(△mi-△mi-1)Τ (43)
其中,J为目标函数的一阶导数,△m为模型参数扰动量。
进一步的,所述步骤九,计算模型的更新方向,更新模型参数,具体包括:
确定模型参数的更新方向,即梯度gk的负方向,第k次反演迭代的梯度表达形式为:
gk=[H(mk)]-1·J(mk) (44)
需按照更新方向计算得到新的模型参数,模型参数的迭代更新公式为:
mk+1=mk-α·gk (45)
其中,α为更新步长。
有益效果:与现有技术相比,本发明的技术方案具有以下有益的技术效果:
相比工业中最常用的叠前AVO反演,本发明方法可以更好利用全波场响应,是一种适用于多层地质模型的反演方法,采用PP和PS多波地震数据的联合反演技术提高弹性参数的估算精度和反演稳定性。
附图说明
图1是本发明反演算法的流程示意图;
图2是测井模型测试数据展示图;
图3是广义传递矩阵模拟得到的PP波和PS波合成记录;
图4是基于广义传递矩阵的PP数据叠前波形反演结果;
图5是基于广义传递矩阵的PS数据叠前波形反演结果;
图6是基于广义传递矩阵的PP、PS联合数据叠前波形反演结果。
具体实施方式
为便于技术人员对本技术的了解,以下将结合附图和实施例对本发明作进一步的描述,并不能用来限制本发明的保护范围。这里以一个测井模型测试为例进行说明:
图1是本发明所述一种多波联合的叠前波形反演方法的实施流程,具体包括以下步骤:
一、广义传递矩阵的PP和PS记录模拟
根据褶积理论,地震记录可由地震子波和反射系数褶积获得,等价于子波矩阵与反射系数序列的乘积。PP和PS合成记录dPP和dPS可由下式得到:
扩展形式为:
式中,d(θi)为入射角θi对应的合成数据,对于PP和PS记录具有不同时间长度,设PP记录采样点数为m1,而PS记录采样点数m2:
从入射角θi的实际地震道集中抽取子波,对于PP记录和PS记录提取的地震子波具有不同长度,设定PP记录的子波长度为l1,PS记录提取子波长度为l1。构建子波矩阵w(θi)如下:
r(θi)为入射角θi时计算得到的反射系数序列:
式中,rPP(tj,θi)和rPS(tj,θi)为广义传递矩阵获得的时间域反射系数,其为频率域反射系数序列通过傅里叶变换得到:
式中,RPP和RPS即为频率域反射系数:
其中,R为频率慢度域反射透射系数向量,b为一维列向量,矩阵形式如下:
RPP、RPS、TPP和TPS分别代表PP波反射系数、PS波反射系数、PP波透射系数和PS波透射系数;n表示除去顶层半空间和底层半空间介质,中间部分包含有n层介质;
A1和A2为层状模型顶层半空间和底层半空间介质的属性矩阵,具体形式如下:
其中,i为虚数单位,ω为角频率,h为地层厚度;
Tj为中间第j和第j+1两相邻地层之间传递矩阵:
式中,c11、c13和c55均为弹性刚度参数,是纵波速度、横波速度和密度的函数,px和pz分别为水平和垂直慢度;
反演的目标参数为纵波速度α、横波速度β和密度ρ,因此建立的初始模型参数向量m如下:
m=[α1,…αM,β1,…βM,ρ1,…ρM] (21)
其中,每种参数共有M个;
刚度参数与纵横波速度和密度之间的数学关系如下:
c11=α2ρ,c13=α2ρ-2β2ρ,c55=β2ρ (22)
二、广义传递矩阵的PP和PS正演算子偏导数
计算广义传递矩阵的偏导数,即非线性正演算子对模型参数的一阶导数,需将频率域偏导数通过傅里叶变换转变为时间域偏导数:
频率域的偏导数包括PP和PS反射系数对纵横波速度、密度三参数的偏导数:
反射系数对模型参数的偏导数,如下:
三、基于广义传递矩阵的PP和PS联合叠前波形反演
广义传递矩阵为非线性正演算法,其PP和PS模拟记录为模型参数的函数:
dPP=FPP(m),dPS=FPS(m) (39)
通过泰勒级数展开并去掉高阶及误差项,则将非线性正演驱动线性化,得到如下表达式:
建立反演目标函数,具体形式如下:
其中,C△m为目标参数m的协方差矩阵,u为模型参数m的期望。α为PP和PS数据的权重参数,λ为正则化参数。针对于PP和PS记录的misfit项,如下:
其中,U、Σ和V为偏导数矩阵的奇异值分解结果,u、δ、v为矩阵参数;N为总的数值个数;给定不同正则化权重参数计算残差,绘制L曲线图,选取拐点位置的权重参数作为最优权重。残差方程具体表达如下:
其中,△mλ用奇异值分解结果表示的形式为:
获得最优的正则化权重参数、合成记录misfit项和正演算子偏导数后,可计算获得目标函数的一阶导数矩阵:
每一次反演迭代中都需要重新求取伪海森矩阵H。本发明采用迭代计算方法更新获得伪海森矩阵,设定最大迭代次数,当达到最大迭代次数,则停止迭代,输出伪海森矩阵。本实施例设定最大迭代次数为5~10次。
首次迭代计算时,伪海森矩阵具体形式为:
Vi=I-ρi(Ji-Ji-1)·(△mi-△mi-1)Τ (53)
其中,J为目标函数的一阶导数,△m为模型参数扰动量。
基于上述结果可计算模型的更新方向,即梯度gk的负方向,第k次反演迭代的梯度表达形式为:
gk=[H(mk)]-1·J(mk) (54)
需按照更新方向计算得到新的模型参数,模型参数的迭代更新公式为:
mk+1=mk-α·gk (55)
更新模型后重新获得PP和PS合成记录,计算PP和PS记录的数据残差。若合成记录残差满足精度要求,则停止反演迭代输出模型参数。若合成记录残差不满足精度,则需要继续参与迭代,计算更新方法,更近模型参数。直到误差满足要求,输出模型结果。
图2为用于反演的测井模型数据,灰色曲线为原始测井数据,黑色曲线为通过低通滤波处理后的测井模型。图3(a)和3(b)为广义传递矩阵算法模拟的PP和PS波合成记录,同时也为后续反演测试的输入数据。图4、图5和图6均为基于广义传递矩阵的叠前波形反演结果,分别为采用PP输入数据、PS输入数据和PP-PS联合数据的反演结果;黑色实线为真实测井模型,黑色点线为输入初始模型,灰色虚线为反演结果。对比可看出PP数据可以获得较好的纵波速度结果,但横波速度和密度结果差些;PS数据则可获得较好的横波速度估算结果,而纵波速度和密度反演结果不理想;PP和PS联合反演则可明显提高三参数反演结果。说明多波联合的叠前波形反演可以获得较理想的三参数反演结果。
Claims (10)
1.一种多波联合的叠前波形反演方法,其特征在于:具体实现步骤如下:
步骤一:获取研究工区的叠前地震数据,即包含全波场响应的叠前道集数据;所述叠前道集包括PP反射波、PS反射波、层间多次波、层间转换波;
步骤二:从叠前道集中分别提取PP波和PS波的地震子波,建立子波矩阵;
步骤三:利用搜集到的测井数据,建立目标弹性参数深度域初始模型;
步骤四:基于深度域初始模型,利用广义传递矩阵正演算子获取PP波和PS波合成记录;
步骤五:将非线性正演算子线性化,建立目标函数,计算PP波和PS波合成记录与实际叠前地震数据的残差;
步骤六:利用初始模型参数计算广义传递矩阵正演算子中PP波和PS波合成记录对模型参数的偏导数;
步骤七:采用自适应估算方法获取最优的正则化权重参数:
步骤八:计算目标函数的一阶导数矩阵,并通过迭代方法获取伪海森矩阵;
步骤九:计算模型的更新方向,更新模型参数;
步骤十:对更新后的模型进行判别,若模型误差未降到预设范围,进行下一次反演迭代,重复步骤四到九;否则迭代停止,输出参数的反演结果。
2.根据权利要求1所述的一种多波联合的叠前波形反演方法,其特征在于:所述步骤一,对叠前道集进行球面扩散补偿、静校正、动校正和表层多次波去除处理。
4.根据权利要求1所述的一种多波联合的叠前波形反演方法,其特征在于:所述步骤三,利用搜集到的测井数据,建立目标弹性参数深度域初始模型;具体如下:
无需对测井数据进行时间-深度匹配,直接对深度域井口数据进行低通滤波,去除测井数据的中高频信息,仅保留低频趋势;获取时间域层位数据解释结果,对时间域结果进行时间-深度匹配,获得深度域层位数据;从井口低频趋势出发,对其沿层位数据进行插值处理,获得深度域2D/3D初始模型。
5.根据权利要求3所述的一种多波联合的叠前波形反演方法,其特征在于:所述步骤四,基于深度域初始模型,利用广义传递矩阵正演算子获取PP波和PS波合成记录;
初始模型参数包括纵波速度α、横波速度β、密度ρ,模型向量m表达式如下:
m=[α1,…αM,β1,…βM,ρ1,…ρM] (4)
其中,每种参数共有M个;
利用初始模型参数,选用广义传递矩阵方法获取频率域的PP和PS反射系数:
式中,R为频率慢度域反射透射系数向量,A1和A2为层状模型顶层半空间和底层半空间介质的属性矩阵,n表示除去顶层半空间和底层半空间介质,中间部分包含有n层介质,Tj为中间第j和第j+1两相邻地层之间传递矩阵,b为一维列向量,式(5)具体形式如下:
其中,RPP、RPS、TPP和TPS分别代表PP波反射系数、PS波反射系数、PP波透射系数和PS波透射系数;
A1和A2表达式如下:
式中,上标(*)1和(*)2分别表示顶层和底层的参数变量;下标(*)P和(*)S则表示纵波和横波的参数变量;i为虚数单位,ω为角频率,h为地层厚度;
Tj表达式如下:
其中,hj表示第j层介质的厚度;
式中,c11、c13和c55均为弹性刚度参数,是纵波速度、横波速度和密度的函数,px和pz分别为水平和垂直慢度;
利用式(5)获得频率域反射系数矩阵RPP和RPS,利用傅里叶变换将频率域反射系数转为时间域反射系数序列:
基于褶积理论,通过反射系数序列与地震子波矩阵获取合成地震记录:
上标(*)syn表示为合成记录数据,dPP表示PP波记录数据,dPS表示PS波记录数据。
6.根据权利要求5所述的一种多波联合的叠前波形反演方法,其特征在于:所述步骤五,将非线性正演算子线性化,建立目标函数,计算PP波和PS波合成记录与实际叠前地震数据的残差;具体如下:
广义传递矩阵为非线性正演算法,其PP和PS模拟记录为模型参数的函数:
dPP=FPP(m),dPS=FPS(m) (18)
通过泰勒级数展开并去掉高阶及误差项,则将非线性正演驱动线性化,得到如下表达式:
上标(*)syn和(*)obs表示为合成记录数据和实际叠前地震道集数据;
建立反演目标函数,具体形式如下:
其中,C△m为模型参数m的协方差矩阵,u为模型参数m的期望,α为PP和PS数据的权重参数,λ为正则化参数,PP和PS记录的misfit项,表示如下:
其中,misfitPP、misfitPS分别表示PP和PS记录的misfit项。
9.根据权利要求8所述的一种多波联合的叠前波形反演方法,其特征在于:所述步骤八,计算目标函数的一阶导数矩阵,并通过迭代方法获取伪海森矩阵,具体包括:
反演的目标函数如式(21)所示,目标函数对模型参数的一阶导数为:
采用迭代计算方法更新获得伪海森矩阵H,设定最大迭代次数,当达到最大迭代次数,则停止迭代,输出伪海森矩阵;首次迭代计算时,伪海森矩阵具体形式为:
Hi+1(m)-1=Vi ΤHi(m)-1Vi+ρi(△mi-△mi-1)·(△mi-△mi-1)Τ (41)
Vi=I-ρi(Ji-Ji-1)·(△mi-△mi-1)Τ (43)
其中,J为目标函数的一阶导数,△m为模型参数扰动量;I为单位矩阵。
10.根据权利要求9所述的一种多波联合的叠前波形反演方法,其特征在于:所述步骤九,计算模型的更新方向,更新模型参数,具体包括:
确定模型参数的更新方向,即梯度gk的负方向,第k次反演迭代的梯度表达形式为:
gk=[H(mk)]-1·J(mk) (44)
按照更新方向计算得到新的模型参数,模型参数的迭代更新公式为:
mk+1=mk-α·gk (45)
其中,α为更新步长。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201911315597.8A CN111025388B (zh) | 2019-12-19 | 2019-12-19 | 一种多波联合的叠前波形反演方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201911315597.8A CN111025388B (zh) | 2019-12-19 | 2019-12-19 | 一种多波联合的叠前波形反演方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN111025388A CN111025388A (zh) | 2020-04-17 |
CN111025388B true CN111025388B (zh) | 2020-10-30 |
Family
ID=70210528
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201911315597.8A Active CN111025388B (zh) | 2019-12-19 | 2019-12-19 | 一种多波联合的叠前波形反演方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN111025388B (zh) |
Families Citing this family (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111818126B (zh) * | 2020-06-08 | 2022-03-11 | 武汉大学 | 基于rfid多特征融合感知模型的物联网环境参数自适应反演方法 |
CN112765802B (zh) * | 2021-01-13 | 2022-11-29 | 陕西师范大学 | 基于高阶水波模型演化水波波形的方法 |
CN117388921A (zh) * | 2023-11-01 | 2024-01-12 | 中国矿业大学(北京) | 弹性参数的叠前反演方法、装置和电子设备 |
Family Cites Families (8)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
AU2006235820B2 (en) * | 2005-11-04 | 2008-10-23 | Westerngeco Seismic Holdings Limited | 3D pre-stack full waveform inversion |
CN102359924B (zh) * | 2011-09-19 | 2013-03-06 | 中国地质大学(北京) | 一种基于多波地震数据的煤岩强度的检测方法 |
CN102928870B (zh) * | 2012-09-21 | 2015-12-02 | 中国石油天然气股份有限公司勘探开发研究院廊坊分院 | 基于正则化的非线性地震叠前弹性参数反演方法 |
CN105093281B (zh) * | 2014-05-16 | 2018-06-29 | 中国石油化工股份有限公司 | 一种反演框架下的地震多波建模方法 |
CA2957818C (en) * | 2014-08-19 | 2023-02-21 | Cgg Services Sa | Joint inversion of compressional and shear seismic data in native time domains |
CN105954804B (zh) * | 2016-07-15 | 2017-12-01 | 中国石油大学(北京) | 页岩气储层脆性地震预测方法及装置 |
CN107656307B (zh) * | 2016-07-26 | 2019-02-19 | 中国石油化工股份有限公司 | 全波形反演计算方法及系统 |
CN109521468B (zh) * | 2018-10-24 | 2021-02-02 | 西南石油大学 | 一种基于Kalman滤波的PP-PS联合反演系统 |
-
2019
- 2019-12-19 CN CN201911315597.8A patent/CN111025388B/zh active Active
Also Published As
Publication number | Publication date |
---|---|
CN111025388A (zh) | 2020-04-17 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN111025388B (zh) | 一种多波联合的叠前波形反演方法 | |
CN106405651B (zh) | 一种基于测井匹配的全波形反演初始速度模型构建方法 | |
US7082368B2 (en) | Seismic event correlation and Vp-Vs estimation | |
CN111025387B (zh) | 一种页岩储层的叠前地震多参数反演方法 | |
CN113031068B (zh) | 一种基于反射系数精确式的基追踪叠前地震反演方法 | |
CN108646293B (zh) | 基于黏声拟微分方程的黏声起伏地表正演模拟系统及方法 | |
CN104122585A (zh) | 基于弹性波场矢量分解与低秩分解的地震正演模拟方法 | |
CN109188519B (zh) | 一种极坐标下的弹性波纵横波速度反演系统及方法 | |
CN110780351B (zh) | 纵波和转换波叠前联合反演方法及系统 | |
CN110187382B (zh) | 一种回折波和反射波波动方程旅行时反演方法 | |
CN110579795B (zh) | 基于被动源地震波形及其逆时成像的联合速度反演方法 | |
CN111045076A (zh) | 多模式瑞雷波频散曲线并行联合反演方法 | |
CN104237937A (zh) | 叠前地震反演方法及其系统 | |
CN105093278A (zh) | 基于激发主能量优化算法的全波形反演梯度算子的提取方法 | |
CN110531410A (zh) | 一种基于直达波场的最小二乘逆时偏移梯度预条件方法 | |
CN111722283A (zh) | 一种地层速度模型建立方法 | |
CN111175821B (zh) | 一种vti介质的各向异性参数分步反演方法 | |
CN111308558B (zh) | 页岩气水平井纵波时差校正方法 | |
CN111273349B (zh) | 一种用于海底浅部沉积层的横波速度提取方法及处理终端 | |
CN116088048A (zh) | 包含不确定性分析的各向异性介质多参数全波形反演方法 | |
CN113031067B (zh) | 一种基于Rytov-WKBJ近似的叠前地震反演方法 | |
CN111308550B (zh) | 一种页岩vti储层的各向异性参数多波联合反演方法 | |
CN111239805B (zh) | 基于反射率法的块约束时移地震差异反演方法及系统 | |
CN117388944A (zh) | 地质模型约束的多物性参数反演方法 | |
CN109613615B (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 |