快速3d自由表面多次波预测方法
【技术领域】
[0001] 本发明属于地震勘探技术中地震信号处理领域，具体涉及一种基于混合频率-时 间域抛物线稀疏反演的快速3D自由表面多次波预测方法。
【背景技术】
[0002] SRME(Surface Related Multiple Elimination)广泛用来压制海洋地震勘探中 的自由表面多次波。2D SRME并不能有效地适应地下介质在横测线方向存在倾角或拖缆 米集存在羽角漂移的情况（Dragoset，B.，E.Verschuur，I.Moore and R. Bisley, 2010, A perspective on 3D surface-related multiple elimination:Geophysics, 75, 75A245-7 5A261.)。考虑到地震波场的3D传播特性，3D SRME比2D SRME能更好地预测自由表面多次 波的旅行时。然而，直接应用3D SRME需要地震波场在横测线和纵测线方向有足够大的采样 密度。实际3D地震采集系统在炮线之间或拖缆线之间有比较大的间距，这会限制3D SRME 方法的直接和广泛应用。
[0003] 目前，已经提出多个方法来解决3D SRME方法面临的数据稀疏采样问题。这些 方法包括在应用3D SRME之前进行数据重建的方法、即时内插方法和横测线多次波贡献 道集（Multiple Contribution Gather/MCG)重建方法等。基于双曲线稀疏反演的横测 线 MCG 重建方法（van Dedem, E. J.，2002, 3D surface-related multiple prediction:PL D.thesis, Delft University of Technology.)假设横测线MCG中的同相轴满足双曲 线时距曲线规律，每一个同相轴的顶点位置表征多次波的真实旅行时。双曲线稀疏反 演方法通过对每一个曲率参数，沿着同相轴顶点的时间和空间方向施加稀疏约束，来估 计同相轴顶点位置。后来，抛物线稀疏反演方法（Hokstad，K，. and R. So 11 i e，2006, 3D surface-related multiple elimination using parabolic sparse inversion:Geophys ics，71，V145-V152.)被提出来，并假设横测线MCG中的同相轴满足抛物线时距曲线规律。 抛物线稀疏反演方法在频率域进行Radon变换，具有比双曲线稀疏反演方法更高的计算效 率。频率域抛物线稀疏反演方法将数据空间和模型空间的关系表达为：
[0010] 其中，d表示数据空间（横测线MCG)，m表示模型空间，L表示Radon变换算子，Nq 表示总的曲率参数个数，Ny表示总的顶点个数，Nd表示横测线MCG中总的地震道数。。。表 示曲率参数，表示顶点位置，y k表示第k道的位置到相应同相轴顶点位置的偏移距。
[0011] 频率域抛物线稀疏反演方法的优化问题为：
[0013] 其中，λ表示正则化因子。频率域抛物线稀疏反演方法并不能在时间方向施加稀 疏约束，相比于时间域双曲线稀疏反演方法，会降低多次波预测的精度。
[0014] 另外，频率域抛物线稀疏反演方法采用迭代重加权最小二乘算法来估计横测线 MCG的同相轴顶点位置。迭代重加权最小二乘算法在预测每一道中的多次波时，需要在每一 次迭代对于每一个频率值，进行一次矩阵-矩阵相乘和矩阵求逆运算，计算量较大。
【发明内容】
[0015] 本发明的目的在于提出一种快速3D自由表面多次波预测方法，该方法利用基于 快速迭代收缩阈值算法的混合频率-时间域抛物线稀疏反演方法来估计横测线MCG的同相 轴顶点位置，并对同相轴顶点位置进行振幅和相位校正后，沿横向相加得到3D自由表面多 次波的预测道，在所有横测线MCG处理完毕后得到所有道的3D自由表面多次波预测结果。
[0016] 为了实现上述目的，本发明采用如下技术方案：
[0017] 快速3D自由表面多次波预测方法，包括如下步骤：
[0018] a对3D地震数据缺失的近偏移距进行外推，并对缺失的地震炮集进行内插，然后 构建相应的横测线MCG ;对横测线MCG进行同相轴顶点估计的数学模型表示为：
[0020] 其中，F表示沿时间方向的傅立叶正变换，F1表示相应的逆傅立叶变换，祖表示时 间域的同相轴顶点位置道集，3表示横测线MCG，L表示Radon变换算子；
[0021] b设置变量初始值；需要设置初始值的变量包括收缩阈值Sa、阻尼因子β、最大迭 代次数Ν、迭代步长l t、曲率参数个数Nq和大小、顶点参数个数N ￥和大小；
[0022] c根据曲率参数个数Nq和大小、顶点参数个数N y和大小构建Radon变换算子L，并 计算矩罔
[0023] 其中，Lh表示Radon变换算子L的共辄转置矩阵，I表示单位矩阵；
[0024] d利用步骤c得到的矩阵E对横测线MCG进行同相轴顶点估计，并得到3D自由表 面多次波的预测道Mp;
[0025] e判断横测线MCG是否全部处理完毕；如果否，返回步骤d ;如果全部处理完毕，输 出3D自由表面多次波的预测结果。
[0026] 进一步，步骤d中3D自由表面多次波的预测道Mp的具体计算过程如下：
[0027] dl设置迭代数η = 0,利用步骤c得到的矩阵E计算同相轴顶点位置的初始值
[0028] d2令η = η+l,对上一步估计的顶点位置进行收缩阈值操作：
[0030] 其中，Ta为收缩阈值算子，定义为：
[0033] 其中，表示第η步估计的顶点位置的收缩阈值结果，涵W表示第η步估 计的顶点位置，Nt表示时间采样点个数，叫.;.?表示向量租.的下标为（i, j, k)的元素，
[0035] d3对收缩阈值结果进行更新：
[0037] 其中，表示更新后的收缩阈值结果，序歹〖
[0038] d4求取更新后的同相轴顶点位置：
[0040] d5判断迭代次数η是否达到最大迭代次数N ;如没达到，返回步骤d2 ;如果达到， 利用同相轴顶点位置的估计结果，按下面公式得到3D自由表面多次波的预测道Mp:
[0042] 其中，F(qu，τ )表示振幅和相位校正因子，/7j(uV)表示同相轴顶点位置的估计 结果，qu表示曲率参数，yv表示顶点位置，τ表示旅行时，ω表示频率，i表示虚数单位。 [0043] 本发明具有如下优点：
[0044] 本发明方法采用混合频率-时间域抛物线稀疏反演方法来估计横测线MCG的同相 轴顶点位置，实现3D自由表面多次波预测。该方法对于每一个曲率参数，能沿着同相轴顶 点的时间和空间方向施加稀疏约束。在已知道位置相同的横测线MCG中，基于快速迭代收 缩阈值算法的混合频率-时间域抛物线稀疏反演方法对每一个频率值，只需计算一次矩阵 求逆。相比于传统的基于频率域抛物线稀疏反演的3D自由表面多次波预测方法，本发明方 法在降低计算量的同时，能进一步提高多次波预测的精度。
【附图说明】
[0045] 图1为本发明中快速3D自由表面多次波预测方法的流程示意图；
[0046] 图2a为原始数据的共炮点道集图；
[0047] 图2b为真实一次波的共炮点道集图；
[0048] 图3a为横测线MCG示意图；
[0049] 图3b为本发明中快速3D自由表面多次波预测方法估计的同相轴顶点示意图，其 中，由图3a顶部的黑色方框标明的地震道作为输入数据；
[0050] 图3c为传统频率域抛物线稀疏反演方法估计的同相轴顶点示意图，其中，由图3a 顶部的黑色方框标明的地震道作为输入数据；
[0051] 图4a为直接利用3D SRME方法得到的预测多次波示意图；
[0052] 图4b为利用本发明快速3D自由表面多次波预测方法得到的预测多次波示意图；
[0053] 图4c为传统的频率域抛物线稀疏反演方法得到的预测多次波示意图；
[0054] 图4d为2D SRME方法得到的预测多次波示意图；
[0055] 图5a为图2a中2050毫秒到2396毫秒的放大显示结果示意图；
[0056] 图5b为图4a中2050毫秒到2396毫秒的放大显示结果示意图；
[0057] 图5c为图4b中2050毫秒到2396毫秒的放大显示结果示意图；
[0058] 图5d为图4c中2050毫秒到2396毫秒的放大显示结果示意图；
[0059] 图5e为图4d中2050毫秒到2396毫秒的放大显示结果示意图；
[0060] 图6a为利用图4a中的预测多次波进行自适应相减的结果示意图；
[0061] 图6b为利用图4b中的预测多次波进行自适应相减的结果示意图；
[0062] 图6c为利用图4c中的预测多次波进行自适应相减的结果示意图；
[0063] 图6d为利用图4d中的预测多次波进行自适应相减的结果示意图。
【具体实施方式】
[0064] 本发明的基本思想为：将横测线MCG的同相轴顶点估计问题表征为混合频率-时 间域抛物线稀疏反演问题，构建相应的数学模型，在频率域进行Radon正反变换，并构建对 每一个曲率参数，沿同相轴顶点的时间和空间方向施加稀疏约束的优化问题，同时引入快 速迭代收缩阈值算法降低优化问题求解的计算复杂度。
[0065] 具体地，