基于激发主能量优化算法的全波形反演梯度算子的提取方法
【技术领域】
[0001] 本发明涉及石油地震勘探速度建模技术领域，尤其涉及一种基于激发主能量优化 算法的全波形反演梯度算子的提取方法。
【背景技术】
[0002] 地震波的传播速度是地震勘探中的一个基本参数，在地震资料处理和地震解释的 诸多环节中扮演着重要的角色。影响地震波传播速度的因素有很多，诸如地下介质的岩性、 密度、孔隙度、埋藏深度等等，这些因素最终体现在地震波的走时、波形、振幅和相位等属性 上变化。只有充分了解速度及影响因素的相互关系，才有可能更为准确地估计出地下的结 构和构造。我国大多数勘探区域越来越复杂，呈现出地表构造和地下构造双复杂特征。传 统的速度分析方法已经很难达到所需的精度要求，也会直接影响成像质量。常用的速度分 析方法有：叠加速度分析、偏移速度分析、层析速度反演和全波形反演。这四种速度分析方 法的准确性依次升高，复杂度也依次升高。
[0003] 全波形反演技术是目前精度很高的速度分析和反演工具，其利用叠前地震记录与 正演模拟炮记录之间的寻求最佳匹配来获取最优反演速度场。在反演过程中不仅考虑了地 震波传播的旅行时等运动学信息，而且加入了振幅、相位、波形等动力学信息，能够适应强 横向变速介质和各向异性介质的速度反演。
[0004] 全波形反演的发展大致可以分为三个阶段，第一个阶段是从20世纪80年代到90 年代，这一时期是全波形反演理论形成和初步验证阶段，Lailly、Tarantola等人建立了基 于波动方程的波形反演的理论构架，之后一批学者分别利用简单模型对反射和透射波形反 演效果进行了理论验证。Gauth i er等人率先对理论模型合成的多偏移距数据进行了声波和 弹性波全波形反演测试，反演结果既展示了全波形反演方法的应用潜力，也暴露出其根本 缺陷--对初始模型的严重依赖性。
[0005] 第二个阶段是从上世纪90年代初到21世纪初，为全波形稳定发展阶段，全波形反 演方法形成了多个分支和变种，其中以Pratt等提出的频率域波形反演方法最具代表性， Pratt提出频率域声波和弹性波波形反演理论，并利用实际井间地震数据对方法进行验证。 该方法利用稀疏矩阵直接解法求解频率域的离散化声波或弹性波波动方程，只要频率与物 性参数不改变，通过一次稀疏矩阵LU分解和多次高斯消元法中的前向、后向回代运算，即 可完成所有炮点的单频波场响应计算；应用研究方面以井间地震波形反演和海洋天然气水 合物的叠前波形反演最具影响力。此外为解决全波形反演因初始模型不准而陷入局部极小 问题，Zhou等人先后给出了波动方程走时和波形联合反演方法，以充分利用走时和波形两 类数据对不同尺度模型参数变化的敏感度的优势，使井间数据全波形反演更为稳定。
[0006] 自21世纪初至今，全波形反演达到了蓬勃发展的阶段，全波形反演逐步从二维反 演走向三维反演，从声波走向弹性波反演，从单参数反演到多参数反演，从各向同性到各向 异性反演。此外还有Laplace/Laplace-Fourier域全波形反演、各向异性全波形反演、先验 信息约束的全波形反演和基于超级炮和相位编码的全波形反演。此外随着计算机水平的进 步，全波形反演的计算效率进一步得到了提升，其实现平台逐步从单机单核走向了多核多 节点并行平台、GPU并行平台。在一定程度上，全波形反演的计算效率完全取决于大数据的 I/O交换。
【发明内容】
[0007] 为了提高时间域全波形反演的计算效率，本发明提供了一种基于激发主能量优化 算法的全波形反演梯度算子的提取方法，其利用有限差分算法计算震源波场和波场误差反 传波场，再根据激发主能量优化算法计算模型更新量的梯度。该方法可以在保障计算精度 的同时，减少大量I/O操作，大大提高全波形反演的计算效率，相比于常规梯度算法，在保 持精度的前提下，提高了近3倍的效率。
[0008] 本发明解决其技术问题所采用的技术方案是：
[0009] -种基于激发主能量优化算法的全波形反演梯度算子的提取方法，包括：
[0010] 建立全波形反演参数模型；
[0011] 利用常规的叠前单炮数据作为观测数据，建立正反演观测系统；
[0012] 将震源置于地表炮点处，利用优化有限差分算法计算该震源产生的正传波场；
[0013] 确定震源波场的正向传播过程中地下各点震源波场最大值及其对应的时刻，并保 存地下各点在激发主能量时窗内的波场；
[0014] 利用优化有限差分模拟算法计算波场残差逆时传播波场，判断地下各点不同时刻 的波场残差逆时传播波场是否位于该点激发主能量时窗之内，通过基于激发主能量优化算 法得到全波形反演梯度算子。
[0015] 前述的基于激发主能量优化算法的全波形反演梯度算子的提取方法，具体包括如 下步骤：
[0016] (1)建立全波形反演参数模型；
[0017] (2)采用常规的叠加前单炮数据作为观测数据，建立正反演观测系统；
[0018] (3)将震源置于地表炮点处，利用下式表示的优化有限差分算法计算该震源产生 的正传波场：
[0019] t-[_ m=.i
_
[0020] 其中U为震源产生的正传波场，Cni为优化有限差分系数
f为介质速 度，Λ X和Λ z为空间采样间隔，Λ t为时间采样间隔，i，k为X和z空间坐标，η为时间坐 标；
[0021] (4)确定震源波场的正向传播过程中地下各点最大值及其对应的时刻Ts ; CV
[0022] (5)定义激发主能量时窗长度为正传波场中一个子波的时间长度w，保存地下各 点在激发主能量时窗内的波场，即以地下各点的#最大值对应时刻Ts为时窗中心的整 Cf 个主能量时窗之内的波场值；
[0023] (6)利用步骤⑶中所述优化有限差分正演模拟算法得到叠前单炮合成数据RU， 并其与观测地震数据Cl cibs求差得到波场残差RU-Clcibs ;
[0024] (7)利用下式表示的优化有限差分模拟算法计算波场残差逆时传播波场：
[0025]
[0026] 其中RU-Clcibs为波场残差，I；为反传算子，Cni为优化有限差分系数，U b为波场残差逆 时传播波场，V为介质速度，Δ X和Δ z为空间采样间隔，Δ t为时间采样间隔，i, k为X和 z空间坐标，η为时间坐标；
[0027] (8)判断地下任意一点的不同时刻波场残差的逆时传播波场是否位于该点激发主 能量时窗之内，如果位于，则按下式提取基于激发主能量优化算法的用于全波形反演的梯 度算子：
[0028]
[0029] 其中RU-Clcibs为波场残差，I；为反传算子，为正演波场对时间的二阶偏导，V为 (；? 空间速度分布，Ts为激发最大振幅对应时间步，w为激发主能量时窗长度，g(v)为全波形反 演的梯度。
[0030] 上述方案更进一步包括：
[0031] 步骤（1)中全波形反演的模型参数求取采用迭代方法求解，迭代公式为
[0032] mk+1 = mk+ Δ m (I)
[0033] 其中mk+1，mk分别表示第k次迭代模型参数和第k+1次迭代模型参数，Λ m是参数 模型的更新量，其表达式为：
[0034] Am = -Qg (2)
[0035] 其中g是参数模型更新量的梯度算子；α是计算步长。所以利用先验信息建立全 波形反演参数的初始模型；
[0036] 步骤（3)中将震源置于地表炮点处，采用优化有限差分算法计算震源激发的正传 波场；
[0037] TF传姑M Π M加下的姑动方IS :<