CN111025387A - 一种页岩储层的叠前地震多参数反演方法 - Google Patents

一种页岩储层的叠前地震多参数反演方法 Download PDF

Info

Publication number
CN111025387A
CN111025387A CN201911315267.9A CN201911315267A CN111025387A CN 111025387 A CN111025387 A CN 111025387A CN 201911315267 A CN201911315267 A CN 201911315267A CN 111025387 A CN111025387 A CN 111025387A
Authority
CN
China
Prior art keywords
matrix
data
model
inversion
reflection coefficient
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
CN201911315267.9A
Other languages
English (en)
Other versions
CN111025387B (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.)
Hohai University HHU
Original Assignee
Hohai University HHU
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 Hohai University HHU filed Critical Hohai University HHU
Priority to CN201911315267.9A priority Critical patent/CN111025387B/zh
Publication of CN111025387A publication Critical patent/CN111025387A/zh
Application granted granted Critical
Publication of CN111025387B publication Critical patent/CN111025387B/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/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

本发明公开了一种适用于页岩等具有强水平层理结构储层的叠前地震反演方法。不同于以往基于等效VTI介质反射系数公式的反演方法，本发明选用了递归矩阵算法作为叠前反演的正演驱动。递归矩阵方法将波动方程精确解与递归数据模式相结合，可有效模拟页岩等具有强层状结构内部复杂的波场效应。基于递归矩阵的叠前反演方法基于各向同性假设，可减少目标参数，提高反演稳定性，可以用于估算页岩等层状结构储层的弹性参数，消除页岩内部复杂波场对下覆地层的影响，恢复下覆地层属性特征。模型参数的合成记录测试证实了本发明方法的有效性。

Description

一种页岩储层的叠前地震多参数反演方法
技术领域
本发明涉及地震勘探技术领域，属于地震资料多参数反演范畴，尤其涉及一种页岩储层的叠前地震多参数反演方法。
背景技术
随着全球范围内对石油和天然气需求的不断增长，油气勘探方向逐步由常规油气储层转向非常规油气储层。非常规储层式相对比常规均质的孔隙型砂岩而言，可以是非常规岩性、非常规的储集空间和非常规的电性特征等。最主要的非常规储层包括有页岩储层、碳酸盐岩储层等。目前，页岩勘探取得了初步的成效，但是大部分地震勘探技术是以往针对常规孔隙型砂岩储层提出的，适用于页岩的技术尚处于发展初期，这也限制了页岩油气储层的有效开采。因此，开展符合页岩储层规律适用于页岩储层的有效的地震预测方法是目前非常规油气研究的重点与难点。
不同于常规厚层砂岩储层，受特殊沉积和应力环境的影响，页岩内部矿物颗粒多呈定向排列，水平层理发育明显，表现出很强的横向各向同性特点。对于储层信息地震预测，叠前AVO/AVA反演是应用较为广泛的一种技术手段，主要是利用地震振幅随偏移距/入射角度变化，即地震AVO/AVA信息去挖掘地下介质的弹性信息。对于横向各向同性属性很强的页岩储层，学者们通常将其按照等效VTI介质来进行研究。针对等效VTI介质，很多学者针对其反射透射系数矩阵进行了研究[1,2],并给出了多种精确表达式[3,4]。为便于实际应用，TI反射系数的近似式[5,6]得到了广泛的研究。在众多近似式中，Rüger近似表达式应用最为广泛。基于该类等效VTI介质的叠前AVO反演方法也逐渐发展，如基于Rüger近似式的各向异性五参数同时、分步反演等[7-9]。这类反演方法虽然可以较好描述页岩等强层理结构储层，但是需要同时反演五个目标参数，且五参数敏感差异很大，易引起反演的不稳定性，很难同时获得准确的五参数结果。
除上述等效VTI理论来模拟页岩储层之外，还可采用递归矩阵正演算法来模拟具有强水平层理的介质。递归矩阵方法采用了波动方程的解析解，利用了递归数学模式精确模拟层状地质模型内部的波场响应[10-13]。有效的还原页岩储层中复杂的多次波、转换波以及强烈的透射衰减等现象。该类算法基于各向同性假设，因此基于此的叠前地震反演仅需要反演纵、横波速度和密度三参数，减少目标参数提高反演稳定性。对于不关注页岩各向异性参数的情况，采用基于递归矩阵的叠前地震反演可以有效估算页岩储层的等效弹性参数结果。同时可消除页岩储层内部复杂响应对下覆底层的影响，较好恢复下覆底层的反演结果。
[1]Thomsen L.Weak elastic anisotropy[J].Geophysics,1986,51(10):1954-1966.
[2]Daley P F,Hron F.Reflection and transmission coefficients forseismic waves in ellipsoidally anisotropic media[J].Geophysics,1979,44(1):27-38.
[3]Graebner M.Plane-wave reflection and transmission coefficients fora transversely isotropic solid[J].Geophysics,1992,57(11):1512-1519.
[4]寻浩，董敏煜，牟永光.横向各向同性介质中的AVO[J].石油地球物理勘探，1997，32(1):45-56.
[5]Thomsen L A.Weak anisotropic reflections,Offset-DependentReflectivity:Theory and Practice of AVO Analysis[M].Library.seg.org,1993:103–114.
[6]Rüger A.P-wave reflection coefficients for transversely isotropicmodels with vertical and horizontal axis of symmetry[J].Geophysics,1997,62(3):713-722.
[7]侯栋甲,刘洋,任志明,等.基于贝叶斯理论的VTI介质多波叠前联合反演[J].石油物探,2014,53(3):294-303.
[8]刘炜，王彦春，谢玮.基于快速非支配排序遗传算法的VTI介质多分量叠前联合反演[J].地球物理学报，62(4)：1453-1470.
[9]Lu J,Wang Y,Chen J,et al.Joint anisotropic amplitude variationwith offset inversion of PP and PS seismic data[J].Geophysics,2018,83(2):N31-N50.
[10]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.
[11]Kennett,B.L.N.1983,Seismic Wave Propagation In Stratified Media:Cambridge University Press.
[12]Booth,D.C.,and S.Crampin.1983,The anisotropic reflectivitytechnique:theory.Geophysical Journal International,72,no.3,755-766.
[13]Carcione,J.M.2001a,Amplitude variations with offset of pressure-seal reflections.Geophysics,66,no.1,283-293.
发明内容
发明目的：针对目前层状地质模型的反演问题，本发明提出一种页岩储层的叠前地震多参数反演方法。该方法克服以往叠前AVA/AVO反演基于单一反射界面的简单假设，适用于页岩等具有强水平层理结构的地质模型，同步且准确获取弹性三参数反演结果，提高地震反演的精度。
技术方案：为实现本发明的目的，本发明所采用的技术方案是：一种页岩储层的叠前地震多参数反演方法，包括以下步骤：
步骤一：获取研究工区的叠前地震数据，此数据是包含全波效应的道集数据；
步骤二：设定需要开展反演的深度范围，利用搜集到的测井数据，建立设定深度范围内的弹性参数深度域初始模型；
步骤三：基于初始模型，利用递归矩阵计算时间角度域反射系数序列；
步骤四：基于叠前地震道集和已取得的反射系数序列，计算初始模型对应的正演模拟记录道集，并计算叠前道集和模拟记录数据残差项；
步骤五：利用初始模型参数计算递归矩阵正演的偏导数；
步骤六：基于目标函数和正演算子偏导数计算目标函数的雅可比矩阵和伪海森矩阵；
步骤七：计算模型的更新方向并更新模型参数；
步骤八：重复步骤三到七，进行下一次反演迭代，直到模型误差降到预设范围，迭代停止，输出参数的反演结果。
进一步的，所述步骤一，获取研究工区的叠前地震数据，此数据是全波场响应数据；输入地震数据需是包含全波场响应的叠前角道集记录，记录中包含有一次反射波、层间多次波、层间转换波。本发明选用递归矩阵方法作为反演的正演算子。递归矩阵是一种结合了波动方程精确解析解和递归算法的正演方法，可精确模拟水平层理结构的全波场记录。基于递归矩阵的反演适用于页岩储层的弹性参数提取。因此，输入地震道集需是全波场道集，提高反演精度。
进一步的，所述步骤二，设定需要开展反演的深度范围，利用搜集到的测井数据，建立设定深度范围内的弹性参数深度域初始模型；具体如下：
不同于射线追踪类的叠前反演方法，基于递归矩阵的叠前反演方法需构建深度域初始模型，反演结果也为深度域参数模型。建立参数深度域初始模型包括：
对深度域测井数据进行低通滤波，获得井口处低频初始模型；
通过时间深度校正后，将时间域层位数据转换为深度域数据；
基于井口初始模型和层位数据，插值处理后获得深度域初始模型。
进一步的，所述插值处理方法为克里金插值处理方法。
进一步的，所述步骤三，基于初始模型，利用递归矩阵计算时间角度域反射系数序列；具体如下：
基于初始模型参数，利用递归矩阵方法获取频率慢度域的反射系数序列，具体为：
C·R＝b (1)
R＝[RPP RPS TPP TPS]T (2)
Figure BDA0002325662850000031
Figure BDA0002325662850000041
式中R＝R(s1，f)为频率慢度域反射透射系数向量，s1表示水平慢度，f表示频率；RPP、RPS、TPP和TPS分别代表PP波和PS波的反射系数、PP波和PS波的透射系数；ηP、ξP、χP和
Figure BDA00023256628500000413
均为待求变量；Pj表示中间任意第j和第j+1两层之间的传递属性矩阵；L1和L2分别代表最顶层和最底层的介质属性矩阵，n表示n层介质：
Figure BDA0002325662850000042
Figure BDA0002325662850000043
下标(*)P和(*)S表示纵波和横波相关的参数，下标(*)1和(*)2表示第一层和最后一层的介质参数；Pj的表达式如下：
Figure BDA0002325662850000044
其中，hj表示第j层的厚度，
Figure BDA0002325662850000045
和
Figure BDA0002325662850000046
可由下方程计算得出：
Figure BDA0002325662850000047
Figure BDA0002325662850000048
Figure BDA0002325662850000049
式中，η、ξ、x和
Figure BDA00023256628500000410
的表达式如下：
Figure BDA00023256628500000411
Figure BDA00023256628500000412
x＝β2ρ(ξs1+ηs3) (13)
Figure BDA0002325662850000051
s1和s3分别表示水平慢度和垂直慢度；垂直慢度s3是水平慢度s1的函数：
Figure BDA0002325662850000052
Figure BDA0002325662850000053
上标P和S分别表示纵波和横波的相应参数；其中，A1、A2和A3分别满足如下：
Figure BDA0002325662850000054
其中，α、β和ρ分别表示纵波速度、横波速度和密度；
基于慢度与入射角θ的数学关系，将频率慢度域反射系数转为频率角度域，具体为：
Figure BDA0002325662850000055
其中，α为入射纵波的速度；通过傅里叶变换，将频率角度域反射系数转换为时间角度域反射系数，具体为：
Figure BDA0002325662850000056
其中r(θ，t)为时间角度域反射系数，R(θ，f)为频率角度域反射系数。
进一步的，所述步骤四，获得初始模型的正演地震道集包括：
利用叠前地震道集数据，提取不同入射角度的地震子波数据；
对反射系数序列和地震子波进行褶积处理得到初始模型的正演地震道集，具体为：
F(m)＝D＝W·r (20)
其中，F表示递归矩阵正演函数，其是模型参数向量m的函数；D为正演合成数据，W为子波矩阵，r为步骤三中得到时间角度域反射系数系列。
式(20)的矩阵形式如下：
Figure BDA0002325662850000057
式中，N表示正演合成数据共包含N个入射角度；d(θi)为入射角θi对应的合成数据；
d(θi)＝[d(t1，θi) d(t2，θi) … d(tM，θi)]T (22)
M表示每个角度数据对应的时间采样点数；从入射角θi的实际地震道集中抽取子波，构建子波矩阵w(θi)：
Figure BDA0002325662850000061
每个子波序列总长度为L，w(lj，θi)表示入射角为θi长度为lj时的子波振幅能量值；r(θi)为入射角θi时计算得到的反射系数序列，通过式(19)计算得到：
r(θi)＝[r(t1，θi) r(t2，θi) … r(tM，θi)]T (24)
进一步的，所述步骤五，利用初始模型参数计算递归矩阵正演的偏导数；
计算递归矩阵对模型的偏导数，其中正演算子偏导数是递归矩阵非线性正演算子对模型参数的一阶导数，具体如下：
Figure BDA0002325662850000062
Figure BDA0002325662850000063
对式(1)求取模型参数m*＝α，β，ρ的偏导数，得到：
Figure BDA0002325662850000064
式中，
Figure BDA0002325662850000065
针对不同深度的反射界面，其表达形式有所不同，第j个反射界面的
Figure BDA0002325662850000066
具体形式为：
Figure BDA0002325662850000067
其中，
Figure BDA0002325662850000068
为求导的中间过程，如下：
Figure BDA0002325662850000071
对式(8)求取模型参数的偏导数，得到
Figure BDA0002325662850000072
Figure BDA0002325662850000073
通过式(11)至(14)，可以得到偏导数
Figure BDA0002325662850000074
和
Figure BDA0002325662850000075
Figure BDA0002325662850000076
Figure BDA0002325662850000077
Figure BDA0002325662850000078
Figure BDA0002325662850000079
对式(9)求取偏导数
Figure BDA0002325662850000081
如下所示：
Figure BDA0002325662850000082
Figure BDA0002325662850000083
h为地层厚度；由式(15)和(16)求得
Figure BDA0002325662850000084
和
Figure BDA0002325662850000085
Figure BDA0002325662850000086
Figure BDA0002325662850000087
Figure BDA0002325662850000088
Figure BDA0002325662850000089
Figure BDA00023256628500000810
上述偏导数均为
Figure BDA00023256628500000811
和
Figure BDA00023256628500000812
的函数。
进一步的，所述步骤六，基于目标函数和正演算子偏导数，计算目标函数的雅可比(Jacobian)矩阵和伪海森(Pseudo Hessian)矩阵；
目标函数的具体形式为：
Figure BDA00023256628500000813
其中，d为实际地震道集数据，F(m)为初始模型的正演合成数据，Km为通过模型数据获取的协方差矩阵，u为通过模型参数的期望，λ为正则化权重参数。
目标函数的雅可比Jacabian矩阵，即为目标函数对模型参数的一阶导数，形式为：
Figure BDA00023256628500000814
其中，
Figure BDA00023256628500000815
为步骤五中获取的递归矩阵正演的偏导数，符号
Figure BDA00023256628500000816
表示函数J(m)对模型参数的一阶偏导数。
选用高斯牛顿的优化算法，获得伪海森矩阵的形式为：
Figure BDA0002325662850000091
符号
Figure BDA0002325662850000092
表示函数J(m)对模型参数的二阶偏导数。
进一步的，所述步骤七，计算模型的更新方向并更新模型参数；
模型参数向量的更新方向为计算获取的梯度的负方向；其中，梯度gk的计算公式为：
Figure BDA0002325662850000093
其中，
Figure BDA0002325662850000094
为目标函数的伪海森(Pseudo Hessian)矩阵，
Figure BDA0002325662850000095
为目标函数的雅可比(Jacabian)矩阵；
按照更新方向计算得到新的模型参数，模型参数的迭代更新公式为：
mk+1＝mk-κ·gk
其中，κ为更新步长。
有益效果：与现有技术相比，本发明的技术方案具有以下有益的技术效果：
以往适用于层状地质模型的叠前波形反演方法，采用了非线性反演策略，计算效率低，计算成本过高；本发明引入递归矩阵作为反演的正演算子，采用了线性反演策略，提高层状地质模型的反演精度，同时提高计算效率。
附图说明
图1是本发明算法的实施流程图；
图2是测试模型的参数值；
图3(a)是使用模型参数的Zoeppritz方程记录模拟结果；
图3(b)是使用模型参数的递归矩阵算法记录模拟结果；
图4是常规基于Zoeppritz方程的叠前反演方法的三参数结果与理论模型对比；
图5是基于递归矩阵的叠前反演方法的三参数结果与理论模型对比。
具体实施方式
为便于技术人员对本技术的了解，以下将结合附图和实施例对本发明的技术方案作进一步的说明，并不能用来限制本发明的保护范围。
本实施例以一个测井模型测试为例进行说明：
图1是本发明一种适用于页岩等具有强水平层理结构模型的叠前多参数反演方法的实施流程，测试模型的参数值如图2所示；具体包括以下步骤：
一、弹性各向同性介质的递归矩阵正演
通过测井数据可获得初始模型参数向量m，如下：
m＝[α1 … αn β1 … βn ρ1 … ρn]T (1)
其中，模型共有n层介质，αi、βi和ρi代表第i层的纵波速度、横波速度和密度参数。基于初始模型参数，利用递归矩阵方法获取频率慢度域的反射系数序列，首先利用入射角和模型参数计算垂直慢度s3：
Figure BDA0002325662850000101
Figure BDA0002325662850000102
上标P和S分别表示纵波和横波的相应参数；其中，A1、A2和A3如下：
Figure BDA0002325662850000103
其中，s1是水平慢度；α、β和ρ分别表示纵波速度、横波速度和密度；
利用模型弹性参数计算中间变量η、ξ、x和
Figure BDA0002325662850000104
Figure BDA0002325662850000105
Figure BDA0002325662850000106
x＝β2ρ(ξs1+ηs3) (7)
Figure BDA0002325662850000107
基于上述结果计算层状模型的顶、底两层介质的属性矩阵：
Figure BDA0002325662850000108
Figure BDA0002325662850000109
h为地层厚度；计算中间多层部分中第j层和第j+1层介质之间的传递矩阵Pj：
Figure BDA0002325662850000111
其中，hj表示第j层的厚度，
Figure BDA0002325662850000112
和
Figure BDA0002325662850000113
可由下方程计算得出：
Figure BDA0002325662850000114
Figure BDA0002325662850000115
其中，
Figure BDA0002325662850000116
基于上述结果可求取水平层状模型的反射系数R＝[RPP RPS TPP TPS]T：
R＝C-1·b (11)
其中，
Figure BDA0002325662850000117
Figure BDA0002325662850000118
式(11)中R＝R(s1，f)为频率慢度域反射透射系数向量；Pj表示中间任意第j和第j+1两层之间的传递属性矩阵，n表示n层介质。基于慢度与入射角的数学关系s1＝sinθ/α，将频率慢度域反射系数转为频率角度域，具体形式为：
Figure BDA0002325662850000119
其中，α为入射纵波的速度；通过傅里叶变换，将频率角度域反射系数转换为时间角度域反射系数，具体形式为：
Figure BDA00023256628500001110
其中r(θ，t)为时间角度域反射系数，R(θ，f)为频率角度域反射系数。
从输入叠前道集数据中估算不同入射角的子波信息，建立子波增广矩阵：
Figure BDA00023256628500001111
w(θi)为入射角θi对应的子波矩阵：
Figure BDA0002325662850000121
每个子波序列总长度为L，w(lj，θi)表示入射角为θi长度为lj时的子波振幅能量值；
计算得到正演合成数据D：
D＝W·r (22)
其中，
D＝[d(θ1) d(θ2) … d(θN)]T (23)
r＝[r(θ1) r(θ2) … r(θN)]T (24)
r为时间角度域反射系数向量。式(22)的矩阵形式如下：
Figure BDA0002325662850000122
式中，N表示正演合成数据共包含N个入射角度，d(θi)为入射角θi对应的合成数据，如下：
d(θi)＝[d(t1，θi) d(t2，θi) … d(tM，θi)]T (26)
M表示每个角度数据对应的时间采样点数；r(θi)为入射角θi时计算得到的反射系数序列，通过式(19)计算得到：
r(θi)＝[r(t1，θi) r(t2，θi) … r(tM，θi)]T (27)
二、基于递归矩阵的叠前反演方法
反演的目标函数的具体形式为：
Figure BDA0002325662850000123
其中，d为实际地震道集数据，Km为通过模型数据获取的协方差矩阵，u为通过模型参数的期望，λ为正则化权重参数。D为初始模型的正演合成数据，这里为递归矩阵方法合成得到F(m)。
梯度类反演方案需按照更新方向计算得到新的模型参数，模型参数的迭代更新公式为：
mk+1＝mk-κ·gk (29)
κ为更新步长。这里采用高斯牛顿方法，针对第k次迭代的更新方法为：
Figure BDA0002325662850000131
其中，
Figure BDA0002325662850000132
为目标函数的雅可比(Jacabian)矩阵，即为目标函数对模型参数的一阶导数，形式为：
Figure BDA0002325662850000133
其中，
Figure BDA0002325662850000134
为正演算子对模型的偏导数，即正演算子F(m)对模型参数的一阶导数，符号
Figure BDA0002325662850000135
表示函数J(m)对模型参数的一阶偏导数。
Figure BDA0002325662850000136
为目标函数的海森Hessian矩阵，即目标函数对模型参数的二阶导数。通常二阶导数很难求得，通常选用其他方法替代Hessian矩阵，即伪海森Pseudo Hessian矩阵。这里选用如下计算方法求取伪海森矩阵：
Figure BDA0002325662850000137
符号
Figure BDA0002325662850000138
表示函数J(m)对模型参数的二阶偏导数。
利用更新得到的新模型mk+1得到新的正演合成数据Dk+1。若合成结果达到误差要求，则停止迭代，输出模型结果。若合成数据未达到误差精度要求，则需要重新计算数据残差：
Δd＝d-F(m) (33)
重新计算偏导数
Figure BDA0002325662850000139
雅可比矩阵
Figure BDA00023256628500001310
和伪海森矩阵
Figure BDA00023256628500001311
从而确定模型更新方向gk，更新模型数据。直到合成数据误差到达要求，输出模型结果。
三、递归矩阵对模型的偏导数
计算递归矩阵对模型的偏导数，其中偏导数是递归矩阵非线性正演算子对模型参数的一阶导数，具体形式如下：
Figure BDA0002325662850000141
即需先计算得到频率角度域的偏导数
Figure BDA0002325662850000142
对其进行傅里叶变换转变时间角度域结果与子波矩阵相乘即可求得正演算子对模型的偏导数
Figure BDA0002325662850000143
其中，
Figure BDA0002325662850000144
对模型参数m*＝α，β，ρ的偏导数具体表现形式如下：
Figure BDA0002325662850000145
Figure BDA0002325662850000146
针对不同反射界面有所不同表达形式，第j个反射界面的具体形式为：
Figure BDA0002325662850000147
其中，
Figure BDA0002325662850000148
Figure BDA0002325662850000149
可有如下方程计算求得：
Figure BDA00023256628500001410
Figure BDA00023256628500001411
是偏导数
Figure BDA00023256628500001412
和
Figure BDA00023256628500001413
的函数，其具体表达形式如下：
Figure BDA0002325662850000151
Figure BDA0002325662850000152
Figure BDA0002325662850000153
Figure BDA0002325662850000154
式(38)中的偏导数
Figure BDA0002325662850000155
有如下表达形式：
Figure BDA0002325662850000156
Figure BDA0002325662850000157
偏导数
Figure BDA0002325662850000158
是导数
Figure BDA0002325662850000159
和
Figure BDA00023256628500001510
的函数。导数
Figure BDA00023256628500001511
Figure BDA00023256628500001512
的计算式已经给出，垂直慢度对模型参数的导数
Figure BDA00023256628500001513
可由下面的方程得到：
Figure BDA0002325662850000161
Figure BDA0002325662850000162
其中，
Figure BDA0002325662850000163
Figure BDA0002325662850000164
Figure BDA0002325662850000165
上述偏导数均为
Figure BDA0002325662850000166
和
Figure BDA0002325662850000167
的函数：
Figure BDA0002325662850000168
Figure BDA0002325662850000169
Figure BDA00023256628500001610
Figure BDA00023256628500001611
图3为精确Zoeppritz方程和递归矩阵算法两种正演方法得到的合成记录，图3(a)是使用模型参数的Zoeppritz方程记录模拟结果；图3(b)是使用模型参数的递归矩阵算法记录模拟结果；可以看出相对比单界面假设的Zoeppritz方程，递归矩阵可以模拟层间反射，箭头处指出了层间效应。图4为基于Zoeppritz的叠前反演结果，输入数据为图3(b)结果，虚线为初始模型；可以看出层间效应会干扰该类反演方法，无法很好反演顶部多层结构以及下覆地层结构。图5为基于递归矩阵的叠前反演方法，输入数据为图3(b)结果，虚线为初始模型；该方法可以有效处理多层结构带来的内部效应，较好表征多层结构同时恢复下覆底层形态。

Claims (9)

1.一种页岩储层的叠前地震反演方法，其特征在于：具体实现步骤如下：
步骤一：获取研究工区的叠前地震数据，此数据是包含全波效应的道集数据；
步骤二：设定需要开展反演的深度范围，利用搜集到的测井数据，建立设定深度范围内的弹性参数深度域初始模型；
步骤三：基于初始模型，利用递归矩阵计算时间角度域反射系数序列；
步骤四：基于叠前地震道集和已取得的反射系数序列，计算初始模型对应的正演模拟记录道集，并计算叠前道集和模拟记录数据残差项；
步骤五：利用初始模型参数计算递归矩阵正演的偏导数；
步骤六：基于目标函数和正演算子偏导数计算目标函数的雅可比矩阵和伪海森矩阵；
步骤七：计算模型的更新方向并更新模型参数；
步骤八：重复步骤三到七，进行下一次反演迭代，直到模型误差降到预设范围，迭代停止，输出参数的反演结果。
2.根据权利要求1所述的方法，其特征在于：所述步骤一中，叠前地震数据需是全波场响应记录，其中包括有一次反射波、层间多次波、层间转换波。
3.根据权利要求1所述的方法，其特征在于：所述步骤二中，建立参数深度域初始模型包括：
对深度域测井数据进行低通滤波，获得井口处低频初始模型；
通过时间深度校正后，将时间域层位数据转换为深度域数据；
基于井口初始模型和层位数据，插值处理后获得深度域初始模型。
4.根据权利要求1所述的方法，其特征在于：所述步骤三中，获得时间角度域反射系数序列包括：
基于初始模型参数，利用递归矩阵方法获取频率慢度域的反射系数序列，具体为：
C·R＝b
R＝[RPP RPS TPP TPS]T，
Figure FDA0002325662840000011
Figure FDA0002325662840000012
式中R＝R(s1,f)为频率慢度域反射透射系数向量，s1表示水平慢度，f表示频率；RPP、RPS、TPP和TPS分别代表PP波和PS波的反射系数、PP波和PS波的透射系数；ηP、ξP、χP和
Figure FDA0002325662840000013
均为待求变量；Pj表示中间任意第j和第j+1两层之间的传递属性矩阵；L1和L2分别代表最顶层和最底层的介质属性矩阵，n表示n层介质；
基于慢度与入射角θ的数学关系，将频率慢度域反射系数转为频率角度域，具体为：
Figure FDA0002325662840000021
其中，α为入射纵波的速度；通过傅里叶变换，将频率角度域反射系数转换为时间角度域反射系数，具体为：
Figure FDA0002325662840000022
其中r(θ,t)为时间角度域反射系数，R(θ,f)为频率角度域反射系数。
5.根据权利要求4所述的方法，其特征在于：所述步骤四中，获得初始模型的正演地震道集包括：
利用叠前地震道集数据，提取不同入射角度的地震子波数据；
对反射系数序列和地震子波进行褶积处理得到初始模型的正演地震道集，具体为：
F(m)＝D＝W·r
其中，F表示递归矩阵正演函数，其是模型参数向量m的函数；D为正演合成数据，W为子波矩阵，r为步骤三中得到时间角度域反射系数系列。
6.根据权利要求5所述的方法，其特征在于：所述步骤五中，计算递归矩阵正演的偏导数，其中该导数是递归矩阵非线性正演算子对模型参数的一阶导数，具体如下：
Figure FDA0002325662840000023
Figure FDA0002325662840000024
不同深度的反射界面
Figure FDA0002325662840000025
表达形式不同，第j个反射界面的
Figure FDA0002325662840000026
具体形式为：
Figure FDA0002325662840000027
7.根据权利要求6所述的方法，其特征在于：所述步骤六中，目标函数为：
Figure FDA0002325662840000028
其中，d为实际地震道集数据，F(m)为初始模型的正演合成数据，Km为通过模型数据获取的协方差矩阵，u为通过模型参数的期望，λ为正则化权重参数；
目标函数的雅可比矩阵，即为目标函数对模型参数的一阶导数，形式为：
Figure FDA0002325662840000029
其中，
Figure FDA0002325662840000031
为步骤五中获取的递归矩阵正演的偏导数，符号
Figure FDA0002325662840000032
表示函数J(m)对模型参数的一阶偏导数。
选用高斯牛顿的优化算法，获得伪海森矩阵的形式为：
Figure FDA0002325662840000033
符号
Figure FDA0002325662840000034
表示函数J(m)对模型参数的二阶偏导数。
8.根据权利要求7所述的方法，其特征在于：所述步骤七中，模型参数向量的更新方向为计算获取的梯度的负方向；其中，梯度gk的计算公式为：
Figure FDA0002325662840000035
其中，
Figure FDA0002325662840000036
为目标函数的伪海森矩阵，
Figure FDA0002325662840000037
为目标函数的雅可比矩阵；
按照更新方向计算得到新的模型参数，模型参数的迭代更新公式为：
mk+1＝mk-κ·gk
其中，κ为更新步长。
9.根据权利要求3所述的方法，其特征在于：所述插值处理采用克里金插值处理方法。
CN201911315267.9A 2019-12-19 2019-12-19 一种页岩储层的叠前地震多参数反演方法 Active CN111025387B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911315267.9A CN111025387B (zh) 2019-12-19 2019-12-19 一种页岩储层的叠前地震多参数反演方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911315267.9A CN111025387B (zh) 2019-12-19 2019-12-19 一种页岩储层的叠前地震多参数反演方法

Publications (2)

Publication Number Publication Date
CN111025387A true CN111025387A (zh) 2020-04-17
CN111025387B CN111025387B (zh) 2020-09-22

Family

ID=70210518

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911315267.9A Active CN111025387B (zh) 2019-12-19 2019-12-19 一种页岩储层的叠前地震多参数反演方法

Country Status (1)

Country Link
CN (1) CN111025387B (zh)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112987088A (zh) * 2021-02-22 2021-06-18 成都理工大学 一种渗流介质地震横波数值模拟和成像方法
CN113031058A (zh) * 2021-02-26 2021-06-25 河海大学 一种基于反射系数精确式的页岩vti储层叠前混合反演方法
CN113156510A (zh) * 2021-04-27 2021-07-23 中南大学 一种页岩储层脆性和各向异性参数预测方法及系统
CN113589385A (zh) * 2021-08-11 2021-11-02 成都理工大学 一种基于地震散射波场分析的储层特征反演方法
CN115586573A (zh) * 2022-09-15 2023-01-10 河海大学 一种致密砂岩储层的动态约束物性参数地震反演方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4953139A (en) * 1990-02-05 1990-08-28 Mobil Oil Corporation Method for restoring and extrapolating seismic traces
CN101644783A (zh) * 2009-07-27 2010-02-10 中国石化集团胜利石油管理局 一种基于条件递推的油藏描述的方法
CN110542928A (zh) * 2018-05-28 2019-12-06 中国石油化工股份有限公司 基于vti各向异性传播矩阵的地震响应模拟方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US4953139A (en) * 1990-02-05 1990-08-28 Mobil Oil Corporation Method for restoring and extrapolating seismic traces
CN101644783A (zh) * 2009-07-27 2010-02-10 中国石化集团胜利石油管理局 一种基于条件递推的油藏描述的方法
CN110542928A (zh) * 2018-05-28 2019-12-06 中国石油化工股份有限公司 基于vti各向异性传播矩阵的地震响应模拟方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
雒聪: "《2019年油气地球物理学术年会论文集》", 30 November 2019 *

Cited By (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112987088A (zh) * 2021-02-22 2021-06-18 成都理工大学 一种渗流介质地震横波数值模拟和成像方法
CN113031058A (zh) * 2021-02-26 2021-06-25 河海大学 一种基于反射系数精确式的页岩vti储层叠前混合反演方法
CN113031058B (zh) * 2021-02-26 2022-06-24 河海大学 一种基于反射系数精确式的页岩vti储层叠前混合反演方法
CN113156510A (zh) * 2021-04-27 2021-07-23 中南大学 一种页岩储层脆性和各向异性参数预测方法及系统
CN113156510B (zh) * 2021-04-27 2022-07-01 中南大学 一种页岩储层脆性和各向异性参数预测方法及系统
CN113589385A (zh) * 2021-08-11 2021-11-02 成都理工大学 一种基于地震散射波场分析的储层特征反演方法
CN113589385B (zh) * 2021-08-11 2023-08-04 成都理工大学 一种基于地震散射波场分析的储层特征反演方法
CN115586573A (zh) * 2022-09-15 2023-01-10 河海大学 一种致密砂岩储层的动态约束物性参数地震反演方法

Also Published As

Publication number Publication date
CN111025387B (zh) 2020-09-22

Similar Documents

Publication Publication Date Title
CN111025387B (zh) 一种页岩储层的叠前地震多参数反演方法
US7082368B2 (en) Seismic event correlation and Vp-Vs estimation
US6901333B2 (en) Method and device for the generation and application of anisotropic elastic parameters
EP3028071B1 (en) Method and device for the generation and application of anisotropic elastic parameters in horizontal transverse isotropic (hti) media
US6839658B2 (en) Seismic processing with general non-hyperbolic travel-time corrections
US11493658B2 (en) Computer-implemented method and system employing nonlinear direct prestack seismic inversion for poisson impedance
CN113031068B (zh) 一种基于反射系数精确式的基追踪叠前地震反演方法
CN110780351B (zh) 纵波和转换波叠前联合反演方法及系统
Luo et al. Joint PP and PS pre-stack seismic inversion for stratified models based on the propagator matrix forward engine
CN101201409B (zh) 一种地震数据变相位校正方法
CN110542928A (zh) 基于vti各向异性传播矩阵的地震响应模拟方法
CN109188511A (zh) 一种砂泥岩薄互层介质多波avo联合反演方法
CN104237937A (zh) 叠前地震反演方法及其系统
CN108828660B (zh) 一种裂缝介质pp波与分裂ps波avo联合反演方法
CN111722284A (zh) 一种基于道集数据建立速度深度模型的方法
CN111175821B (zh) 一种vti介质的各向异性参数分步反演方法
CN111025388B (zh) 一种多波联合的叠前波形反演方法
CN110187389A (zh) 一种基于薄层反射理论的ava反演方法
CN111308550B (zh) 一种页岩vti储层的各向异性参数多波联合反演方法
CN113253350B (zh) 一种基于联合字典的孔隙度反演方法
CN117388944A (zh) 地质模型约束的多物性参数反演方法
CN112305595B (zh) 基于折射波分析地质体结构的方法及存储介质
Wu et al. An unsupervised inversion method for seismic brittleness parameters driven by the physical equation
CN110764146A (zh) 一种基于声波算子的空间互相关弹性波反射波形反演方法
CN113031058B (zh) 一种基于反射系数精确式的页岩vti储层叠前混合反演方法

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