技术领域
[0001] 本
发明涉及
地震勘探技术领域，属于地震资料多参数反演范畴，尤其涉及一种多波联合的叠前波形反演方法。
背景技术
[0002] 目前，叠前地震反演方法主要包括两类，一是基于射线追踪类正演算子的叠前AVO反演方法，二是基于
波动方程类正演算子的叠前波形反演[1]。第一类叠前AVO反演优势在于单一界面简单假设，数学形式相对简单，具有较高的计算效率，因此在实际
地震勘探中应用最为广泛。不同于叠前AVO反演，叠前波形反演采用了波动方程类模拟
算法作为正演算子，能更好利用实际数据中的全波场响应(如多次波和转换波等)开展
地层介质参数反演估算。
[0003] 波动方程类正演包括有反射率[2-5]和传递矩阵[6,7]等方法。反射率法目前被用于叠前波形反演的正演算子，并取得了理想的估算结果[8，10-11]。绝大多数现存的叠前波形反演方法采用非线性的反演方法，如选用模拟
退火[9]、
遗传算法等优化策略[10-13]。模型测试、2D和3D实际数据应用均表明了基于非线性
优化算法的叠前波形反演在介质属性估算方面的优势。非线性反演方法不依赖于初始模型，无需求取偏导数实现简单，同时有助于获取全局最优解，避免落入局部极值。然而，该类算法需要多次重复正演过程，波动方程类正演计算复杂，因此采用非线性策略的叠前波形反演计算量过大，计算成本过高。因此，为了解决这一问题，需开展基于线性反演策略的叠前波形反演研究，降低计算成本，更易于向实际数据应用中推广。
[0004] 针对各向同性介质的叠前地震反演，通常反演纵
波速度、横波速度和
密度三参数。输入数据不足会导致多参数反演的不稳定。为了提高反演的
稳定性和
精度，多
波数据联合反演能有效提高参数估算精度[11，13-15]。密度是开展岩性识别和油气勘探的重要参数。
Luo等指出，单纯的PP数据反演无法准确地获得横波速度和密度结果[16]。即使在较大入射
角情况，由于密度对PP反射振幅不够敏感，密度估算的不精确性仍然是一个问题。转换波(PS数据)包含丰富的横波速度和密度信息。因此，PP和PS联合反演可以提供更好的反演结果。因此，为了提高叠前波形反演的横波速度和密度的估算精度，需要开展多波联合的叠前波形反演方法研究。
[0005] 参考文献
[0006] [1]Simmons,J.L.,and M.M.Backus.1994,AVO modeling and the locally converted shear wave.Geophysics,59,no.8,1237-1248.
[0007] [2]Fryer,G.J.,and L.N.Frazer.1984,Seismic waves in stratified anisotropic media.Geophysical Journal of the Royal Astronomical Society,78,no.3,691-710.
[0008] [3]Fuchs,K.,and G.Müller.1971,Computation of synthetic seismograms with the reflectivity method and comparison with observations.Geophysical Journal International,23,no.4,417-433.
[0009] [4]Kennett,B.L.N.1983,Seismic Wave Propagation In Stratified Media:Cambridge University Press.
[0010] [5]Booth,D.C.,and S.Crampin.1983,The anisotropic reflectivity technique:theory.Geophysical Journal International,72,no.3,755-766.[0011] [6]Carcione,J.M.2001a,Amplitude variations with offset of pressure-seal reflections.Geophysics,66,no.1,283-293.
[0012] [7]Carcione,J.M.2001b,Wave fields in real media:Wave propagation in anisotropic,anelastic,porous and electromagnetic media.Vol.38:Elsevier.[0013] [8]Sen,M.K.,and P.L.Stoffa.1991,Nonlinear seismic waveform inversion in one dimension using simulated annealing.Geophysics,56,no.10,1624-1638.[0014] [9]Pafeng,J.,S.Mallick,and H.Sharma.2016,Prestack waveform inversion of three-dimensional seismic data—An example from the Rock Springs Uplift,Wyoming,USA.Geophysics,82,no.1,B1-B12.
[0015] [10]Padhi,A.,and S.Mallick.2013,Accurate estimation of density from the inversion of multicomponent prestack seismic waveform data using a nondominated sorting genetic algorithm.Geophysics,32,no.1,94-98.
[0016] [11]Padhi,A.,and S.Mallick.2014,Multicomponent pre-stack seismic waveform inversion in transversely isotropic media using a non-dominated sorting genetic algorithm.Geophysical Journal International,196(3):1600-1618.[0017] [12]Gisolf,A.,P.Haffinger,C.Hanitzsch,P.Doulgeris,and P.Veeken.2014,Non-linear full wavefield inversion applied to carboniferous reservoirs in the Wingate gas field(SNS,Offshore UK).76th EAGE Conference and Exhibition 2014.
[0018] [13]Li,T.,and S.Mallick.2015,Multicomponent,multi-azimuth pre-stack seismic waveform inversion for azimuthally anisotropic media using a parallel and computationally efficient non-dominated sorting genetic algorithm.Geophysical Journal International,200,no.2,1136-1154.
[0019] [14]Lu,J.,Y.Wang,J.Chen,and Y.An.2018,Joint anisotropic amplitude variation with offset inversion of PP and PS seismic data.Geophysics,83,no.2,N31-N50..
[0020] [15]Lu,J.,Z.Yang,Y.Wang,and Y.Shi.2015a,Joint PP and PS AVA seismic inversion using exact Zoeppritz equations.Geophysics,80,no.5,R239-R250.[0021] [16]Luo,C.,X.Y.Li,and G.Huang.2018b,Information on S-wave velocity and density obtained from multi-wave data.SEG Technical Program Expanded Abstracts.
发明内容
[0022] 发明目的：为克服以往叠前波形非线性反演方法的计算效率过低问题，本发明提出一种多波联合的叠前波形反演方法，选用了线性反演策略，利用全波场响应多波记录的叠前波形线性反演方法，更多利用多波信息和地震记录的转换波、多次波信息，提高弹性参数反演精度。
[0023] 技术方案：为实现本发明的目的，本发明所采用的技术方案是：一种多波联合的叠前波形反演方法，包括以下步骤：
[0024] 步骤一：获取研究工区的叠前地震数据，即包含全波场响应的PP和PS道集数据；正演方法的合成记录应准确模拟实际全波场数据，两者一致性好有利于获得真实的反演结果；广义传递矩阵算法可以精确模拟层状模型内部复杂的波场响应，因此叠前PP和PS道集中应包括PP反射波、PS反射波、层间多次波、层间转换波；优选的，叠前道集进行球面扩散补偿、静校正、动校正和表层多次波去除处理；
[0025] 步骤二：从叠前道集中分别提取PP波和PS波的地震子波，建立子波矩阵；
[0026] 步骤三：利用搜集到的
测井数据，建立目标弹性参数深度域初始模型；
[0027] 步骤四：基于深度域初始模型，利用广义传递矩阵正演算子获取PP波和PS波合成记录；
[0028] 步骤五：将非线性正演算子线性化，建立目标函数，计算PP波和PS波合成记录与实际叠前地震数据的残差；
[0029] 步骤六：利用初始模型参数计算广义传递矩阵正演算子中PP波和PS波合成记录对模型参数的偏导数；
[0030] 步骤七：采用自适应估算方法获取最优的正则化权重参数：
[0031] 步骤八：计算目标函数的一阶导数矩阵，并通过
迭代方法获取伪海森矩阵；
[0032] 步骤九：计算模型的更新方向，更新模型参数；
[0033] 步骤十：对更新后的模型进行判别，若模型误差未降到预设范围，进行下一次反演迭代，重复步骤四到九；否则迭代停止，输出参数的反演结果。
[0034] 进一步的，所述步骤二，从叠前道集中分别提取PP波和PS波的地震子波，建立子波矩阵，具体包括：
[0035] 从叠前PP波角道集、PS波角道集中提取不同角度情况的地震子波序列，分别建立PP记录、PS记录的子波矩阵WPP(θ)、WPS(θ)；表达式如下：
[0036]
[0037]
[0038] 其中，θ表示角度，N表示角道集中角度数量；
[0039] 建立PP和PS联合记录的子波矩阵，用于制作PP和PS联合合成记录，如下：
[0040]
[0041] 进一步的，所述步骤三，利用搜集到的测井数据，建立目标弹性参数深度域初始模型；方法如下：
[0042] 不同于基于Zoeppritz方程的叠前AVO/AVA反演，基于广义传递矩阵的叠前波形反演需构建深度域初始模型。建立深度域初始模型具体步骤包括：
[0043] 无需对测井数据进行时间-深度匹配，直接对深度域井口数据进行低通滤波，去除测井数据的中高频信息，仅保留低频趋势；获取时间域层位数据解释结果，对时间域结果进行时间-深度匹配，获得深度域层位数据；从井口低频趋势出发，对其沿层位数据进行插值处理，获得深度域2D/3D初始模型。
[0044] 进一步的，所述步骤四，基于深度域初始模型，利用广义传递矩阵正演算子获取PP波和PS波合成记录；
[0045] 初始模型参数包括纵波速度α、横波速度β、密度ρ，模型向量m表达式如下：
[0046] m＝[α1,…αM,β1,…βM,ρ1,…ρM] (4)
[0047] 其中，每种参数共有M个；
[0048] 利用初始模型参数，选用广义传递矩阵方法获取
频率域的PP和PS反射系数：
[0049]
[0050] 式中，R为频率慢度域反射透射系数向量，A1和A2为层状模型顶层半空间和底层半空间介质的属性矩阵，n表示除去顶层半空间和底层半空间介质，中间部分包含有n层介质，Tj为中间第j和第j+1两相邻地层之间传递矩阵，b为一维列向量，式(5)具体形式如下：
[0051]
[0052] 其中，RPP、RPS、TPP和TPS分别代表PP波反射系数、PS波反射系数、PP波透射系数和PS波透射系数；
[0053] A1和A2具体形式如下：
[0054]
[0055]
[0056] 式中，上标(*)1和(*)2分别表示顶层和底层的参数变量；下标(*)P和(*)S则表示纵波和横波的参数变量；i为虚数单位，ω为角频率，h为地层厚度；
[0057] Tj表达式如下：
[0058]
[0059]
[0060]
[0061]
[0062]
[0063]
[0064] 式中，hj表示第j层介质的厚度；c11、c13和c55均为弹性
刚度参数，是纵波速度、横波速度和密度的函数，px和pz分别为
水平和垂直慢度；
[0065] 利用式(5)获得频率域反射系数矩阵RPP和RPS，利用傅里叶变换将频率域反射系数转为时间域反射系数序列：
[0066]
[0067]
[0068] 基于褶积理论，通过反射系数序列与地震子波矩阵(见式(1)和式(2))可获取合成地震记录：
[0069]
[0070] 上标(*)syn表示为合成记录数据，dPP表示PP波记录数据，dPS表示PS波记录数据。
[0071] 进一步的，所述步骤五，将非线性正演算子线性化，建立目标函数，计算PP波和PS波合成记录与实际叠前地震数据的残差；
[0072] 广义传递矩阵为非线性正演算法，其PP和PS模拟记录为模型参数的函数：
[0073] dPP＝FPP(m)，dPS＝FPS(m) (18)
[0074] 通过泰勒级数展开并去掉高阶及误差项，则将非线性正演驱动线性化，得到如下表达式：
[0075]
[0076] 其中， 为广义传递矩阵合成记录对模型参数的偏导数，△m为模型扰动量；△dPP和△dPS为PP和PS记录的扰动量，即合成数据与实际记录的残差：
[0077]
[0078] 上标(*)syn和(*)obs表示为合成记录数据和实际叠前
地震道集数据；
[0079] 建立反演目标函数，具体形式如下：
[0080]
[0081] 其中，C△m为模型参数m的协方差矩阵，u为模型参数m的期望，α为PP和PS数据的权重参数，λ为正则化参数，PP和PS记录的misfit项，表示如下：
[0082]
[0083]
[0084] 其中，misfitPP、misfitPS分别表示PP和PS记录的misfit项。
[0085] 进一步的，所述步骤六，计算广义传递矩阵算法PP和PS记录对模型参数的偏导数，利用初始模型参数计算其初始偏导数结果；
[0086] 计算广义传递矩阵的偏导数，即非线性正演算子对模型参数的一阶导数，具体形式如下：
[0087]
[0088]
[0089]
[0090] 反射系数对模型参数的偏导数，如下：
[0091]
[0092] 其中， 为求导的中间过程偏导数矩阵：
[0093]
[0094]
[0095]
[0096] 式(29)和(30)中的偏导数 为：
[0097]
[0098]
[0099]
[0100]
[0101] 式(30)至(34)中的 为慢度参数对模型的偏导数；此外，为刚度参数对模型参数的导数； 为密度偏导数。
[0102] 进一步的，所述步骤七，利用L曲线自适应估算策略获取最优的正则化权重参数：
[0103] 自适应正则化权重参数获取包括：
[0104] 对PP和PS合成记录的正演算子偏导数，即 和 进行奇异值分解：
[0105]
[0106]
[0107] 其中，U、Σ和V为偏导数矩阵的奇异值分解结果，u、δ、v为矩阵参数；N为总的数值个数；给定不同正则化权重参数计算残差，绘制L曲线图，选取拐点
位置的权重参数作为最优权重；残差方程具体表达如下：
[0108]
[0109] 其中，△mλ用奇异值分解结果表示的形式为：
[0110]
[0111] 进一步的，所述步骤八，计算目标函数的一阶导数矩阵，并通过迭代方法获取伪海森矩阵，具体包括：
[0112] 目标函数的具体形式见式(21)，目标函数对模型参数的一阶导数为：
[0113]
[0114] 其中，符号 表示对模型参数求取一阶导数； 为步骤六得到的PP和PS正演算子对模型参数的偏导数；
[0115] 每一次反演迭代中都需要重新求取伪海森矩阵H。本发明采用迭代计算方法更新获得伪海森矩阵，设定最大迭代次数，当达到最大迭代次数，则停止迭代，输出伪海森矩阵。首次迭代计算时，伪海森矩阵具体形式为：
[0116]
[0117] 式中，符号 表示对模型参数求取二阶导数；第i次迭代的伪海森矩阵形式如下：
[0118]
[0119]
[0120] Vi＝I-ρi(Ji-Ji-1)·(△mi-△mi-1)Τ (43)[0121] 其中，J为目标函数的一阶导数，△m为模型参数扰动量。
[0122] 进一步的，所述步骤九，计算模型的更新方向，更新模型参数，具体包括：
[0123] 确定模型参数的更新方向，即梯度gk的负方向，第k次反演迭代的梯度表达形式为：
[0124] gk＝[H(mk)]-1·J(mk) (44)
[0125] 需按照更新方向计算得到新的模型参数，模型参数的迭代更新公式为：
[0126] mk+1＝mk-α·gk (45)
[0127] 其中，α为更新步长。
[0128] 有益效果：与
现有技术相比，本发明的技术方案具有以下有益的技术效果：
[0129] 相比工业中最常用的叠前AVO反演，本发明方法可以更好利用全波场响应，是一种适用于多层地质模型的反演方法，采用PP和PS多波地震数据的联合反演技术提高弹性参数的估算精度和反演稳定性。
附图说明
[0130] 图1是本发明反演算法的流程示意图；
[0131] 图2是测井模型测试数据展示图；
[0132] 图3是广义传递矩阵模拟得到的PP波和PS波合成记录；
[0133] 图4是基于广义传递矩阵的PP数据叠前波形反演结果；
[0134] 图5是基于广义传递矩阵的PS数据叠前波形反演结果；
[0135] 图6是基于广义传递矩阵的PP、PS联合数据叠前波形反演结果。
具体实施方式
[0136] 为便于技术人员对本技术的了解，以下将结合附图和
实施例对本发明作进一步的描述，并不能用来限制本发明的保护范围。这里以一个测井模型测试为例进行说明：
[0137] 图1是本发明所述一种多波联合的叠前波形反演方法的实施流程，具体包括以下步骤：
[0138] 一、广义传递矩阵的PP和PS记录模拟
[0139] 根据褶积理论，地震记录可由地震子波和反射系数褶积获得，等价于子波矩阵与反射系数序列的乘积。PP和PS合成记录dPP和dPS可由下式得到：
[0140]
[0141] 扩展形式为：
[0142]
[0143] 式中，d(θi)为入射角θi对应的合成数据，对于PP和PS记录具有不同时间长度，设PP记录
采样点数为m1，而PS记录采样点数m2：
[0144]
[0145]
[0146] 从入射角θi的实际地震道集中
抽取子波，对于PP记录和PS记录提取的地震子波具有不同长度，设定PP记录的子
波长度为l1，PS记录提取子波长度为l1。构建子波矩阵w(θi)如下：
[0147]
[0148]
[0149] r(θi)为入射角θi时计算得到的反射系数序列：
[0150]
[0151]
[0152] 式中，rPP(tj,θi)和rPS(tj,θi)为广义传递矩阵获得的时间域反射系数，其为频率域反射系数序列通过傅里叶变换得到：
[0153]
[0154]
[0155] 式中，RPP和RPS即为频率域反射系数：
[0156]
[0157] 其中，R为频率慢度域反射透射系数向量，b为一维列向量，矩阵形式如下：
[0158]
[0159] RPP、RPS、TPP和TPS分别代表PP波反射系数、PS波反射系数、PP波透射系数和PS波透射系数；n表示除去顶层半空间和底层半空间介质，中间部分包含有n层介质；
[0160] A1和A2为层状模型顶层半空间和底层半空间介质的属性矩阵，具体形式如下：
[0161]
[0162]
[0163] 其中，i为虚数单位，ω为角频率，h为地层厚度；
[0164] Tj为中间第j和第j+1两相邻地层之间传递矩阵：
[0165]
[0166]
[0167]
[0168]
[0169]
[0170] 式中，hj表示第j层介质的厚度；而σ和δ均为φ和 的函数，形式如下：
[0171]
[0172] 式中，c11、c13和c55均为弹性刚度参数，是纵波速度、横波速度和密度的函数，px和pz分别为水平和垂直慢度；
[0173] 反演的目标参数为纵波速度α、横波速度β和密度ρ，因此建立的初始模型参数向量m如下：
[0174] m＝[α1,…αM,β1,…βM,ρ1,…ρM] (21)[0175] 其中，每种参数共有M个；
[0176] 刚度参数与纵横波速度和密度之间的数学关系如下：
[0177] c11＝α2ρ，c13＝α2ρ-2β2ρ，c55＝β2ρ (22)[0178] 二、广义传递矩阵的PP和PS正演算子偏导数
[0179] 计算广义传递矩阵的偏导数，即非线性正演算子对模型参数的一阶导数，需将频率域偏导数通过傅里叶变换转变为时间域偏导数：
[0180]
[0181]
[0182] 频率域的偏导数包括PP和PS反射系数对纵横波速度、密度三参数的偏导数：
[0183]
[0184] 反射系数对模型参数的偏导数，如下：
[0185]
[0186] 其中， 为求导的中间过程偏导数矩阵：
[0187]
[0188]
[0189]
[0190] 式(28)和(29)中的偏导数 为：
[0191]
[0192]
[0193]
[0194]
[0195] 式(30)至(33)中的 如下：
[0196]
[0197]
[0198] 此外，刚度参数的偏导数 和密度偏导数 为：
[0199]
[0200]
[0201]
[0202] 三、基于广义传递矩阵的PP和PS联合叠前波形反演
[0203] 广义传递矩阵为非线性正演算法，其PP和PS模拟记录为模型参数的函数：
[0204] dPP＝FPP(m)，dPS＝FPS(m) (39)
[0205] 通过泰勒级数展开并去掉高阶及误差项，则将非线性正演驱动线性化，得到如下表达式：
[0206]
[0207] 其中， 为广义传递矩阵合成记录对模型参数的偏导数，△m为模型扰动量；△dPP和△dPS为PP和PS记录的扰动量，即合成数据与实际记录的残差：
[0208]
[0209] 建立反演目标函数，具体形式如下：
[0210]
[0211] 其中，C△m为目标参数m的协方差矩阵，u为模型参数m的期望。α为PP和PS数据的权重参数，λ为正则化参数。针对于PP和PS记录的misfit项，如下：
[0212]
[0213]
[0214] 自适应正则化权重参数获取，需对PP和PS合成记录的正演算子偏导数，即和 进行奇异值分解：
[0215]
[0216]
[0217] 其中，U、Σ和V为偏导数矩阵的奇异值分解结果，u、δ、v为矩阵参数；N为总的数值个数；给定不同正则化权重参数计算残差，绘制L曲线图，选取拐点位置的权重参数作为最优权重。残差方程具体表达如下：
[0218]
[0219] 其中，△mλ用奇异值分解结果表示的形式为：
[0220]
[0221] 获得最优的正则化权重参数、合成记录misfit项和正演算子偏导数后，可计算获得目标函数的一阶导数矩阵：
[0222]
[0223] 其中，符号 表示对模型参数求取一阶导数； 为PP和PS正演算子对模型参数的偏导数。
[0224] 每一次反演迭代中都需要重新求取伪海森矩阵H。本发明采用迭代计算方法更新获得伪海森矩阵，设定最大迭代次数，当达到最大迭代次数，则停止迭代，输出伪海森矩阵。本实施例设定最大迭代次数为5～10次。
[0225] 首次迭代计算时，伪海森矩阵具体形式为：
[0226]
[0227] 式中，符号 表示对模型参数求取二阶导数；第i次迭代的伪海森矩阵形式如下：
[0228]
[0229]
[0230] Vi＝I-ρi(Ji-Ji-1)·(△mi-△mi-1)Τ (53)[0231] 其中，J为目标函数的一阶导数，△m为模型参数扰动量。
[0232] 基于上述结果可计算模型的更新方向，即梯度gk的负方向，第k次反演迭代的梯度表达形式为：
[0233] gk＝[H(mk)]-1·J(mk) (54)
[0234] 需按照更新方向计算得到新的模型参数，模型参数的迭代更新公式为：
[0235] mk+1＝mk-α·gk (55)
[0236] 更新模型后重新获得PP和PS合成记录，计算PP和PS记录的数据残差。若合成记录残差满足精度要求，则停止反演迭代输出模型参数。若合成记录残差不满足精度，则需要继续参与迭代，计算更新方法，更近模型参数。直到误差满足要求，输出模型结果。
[0237] 图2为用于反演的测井模型数据，灰色曲线为原始测井数据，黑色曲线为通过低通滤波处理后的测井模型。图3(a)和3(b)为广义传递矩阵算法模拟的PP和PS波合成记录，同时也为后续反演测试的输入数据。图4、图5和图6均为基于广义传递矩阵的叠前波形反演结果，分别为采用PP输入数据、PS输入数据和PP-PS联合数据的反演结果；黑色实线为真实测井模型，黑色点线为输入初始模型，灰色虚线为反演结果。对比可看出PP数据可以获得较好的纵波速度结果，但横波速度和密度结果差些；PS数据则可获得较好的横波速度估算结果，而纵波速度和密度反演结果不理想；PP和PS联合反演则可明显提高三参数反演结果。说明多波联合的叠前波形反演可以获得较理想的三参数反演结果。
