考虑地震子波拉伸效应的变顶点稀疏双曲线Radon变换方法

摘要

本发明公开了一种考虑地震子波拉伸效应的变顶点稀疏双曲线Radon变换方法。该方法利用地震同相轴的倾角及其导数来计算双曲线同相轴的顶点位置和速度，并利用聚类算法来确定最终的双曲线同相轴个数，降低所估计的同相轴参数个数；利用所计算的同相轴顶点位置和速度，采用开窗分割地震道的方法来构建变顶点双曲线Radon变换算子，消除Radon反变换产生的地震子波拉伸效应；然后在优化问题中引入Radon域模型和数据拟合项的L1范数最小化约束，来提高Radon域模型的稀疏性和对非高斯拟合误差的适应性；另外，本发明方法利用乘数交替方向方法求解L1‑L1范数最小化优化问题，在整个迭代过程中只需一次矩阵求逆，计算复杂度较低。

说明书

技术领域

本发明属于地震勘探技术中的地震信号处理领域，具体涉及一种考虑地震子波拉伸效应的变顶点稀疏双曲线Radon变换方法。

背景技术

Radon变换广泛应用于地震信号处理中，如多次波压制、地震数据规则化和波场分离等，其主要思想是将地震数据从时间-偏移距域映射到Radon域。Radon变换按照积分轨迹区分为三大类：线性Radon变换、抛物线Radon变换和双曲线Radon变换。

其中，双曲线Radon变换可直接应用于未做动校正的共中心点道集或共炮点道集，并利用双曲线轨迹来构造变换算子，能够更好地表示地震数据。

由于双曲线Radon变换的时变特性，计算双曲线Radon变换需要对大型稀疏双曲线Radon变换算子进行矩阵求逆，导致时间域双曲线Radon变换过程中的巨大计算量。另外，传统的双曲线Radon变换将顶点位置限制在零偏移距处，利用不同的速度参数来表征地震道集中的双曲线轨迹。当双曲线同相轴的顶点位于非零偏移距位置，如道集中存在微曲多次波、绕射多次波，传统的双曲线Radon变换进行多次波压制或数据规则化时，其效果会变差。

双曲线Radon变换的应用效果依赖于变换算子是否能有效表征道集中存在的双曲线同相轴。对于双曲线同相轴的顶点沿偏移距移动的情况，变顶点双曲线Radon变换方法通过将Radon域模型的维度空间从2D(顶点时间和速度)扩展为3D(顶点时间、速度和顶点偏移距)，能够更有效地表示顶点位于非零偏移距的双曲线同相轴。由于Radon域模型的空间维度增加一维，导致在Radon域模型中需要求解的未知参数个数变多，增加了Radon正反变换的计算量。另外，变顶点双曲线Radon变换方法对应的Radon变换方程组为欠定方程组，需要在相应的优化问题中引入Radon域模型的稀疏性约束来计算Radon域模型。传统的变顶点稀疏双曲线Radon变换方法采用迭代重加权最小二乘算法求解L1范数最小化优化问题，在每一步迭代中需要计算一次矩阵求逆，计算复杂度较高。因此，传统变顶点稀疏双曲线Radon变换的计算量较大，这限制了该变换方法的实用性。另外，传统的双曲线Radon变换在构建双曲线Radon变换算子进行正反变换时，并没有考虑地震子波的存在，在利用Radon反变换进行地震数据重构过程中，会出现地震子波拉伸效应，造成地震信号失真。

发明内容

针对现有技术中存在的上述技术问题，本发明提出了一种考虑地震子波拉伸效应的变顶点稀疏双曲线Radon变换方法，其采用如下技术方案：

考虑地震子波拉伸效应的变顶点稀疏双曲线Radon变换方法，包括如下步骤：

a、设置变量初始值：需要设置初始值的变量包括阻尼因子β、乘数交替方向方法的阈值χ和φ、乘数交替方向方法的最大迭代次数N；

b、从所有地震道集中选取某一个地震道集，利用复地震道分析方法估计地震同相轴的倾角p_r，并计算倾角导数p_rr；

c、利用下式计算双曲线同相轴的顶点位置(τ,x₀)和速度v：

x₀＝x_r-tp_rv²，其中，τ表示顶点位置中的时间位置，x₀表示顶点位置中的空间位置，t表示时间，x_r表示检波点的空间位置；

d、在利用步骤c得到双曲线同相轴的顶点位置和速度后，利用聚类算法来检测双曲线同相轴的个数，确定最终的同相轴顶点位置和速度；

e、利用步骤d得到的同相轴顶点位置和速度，确定双曲线同相轴的轨迹，然后，利用开窗分割地震道的方法沿着垂直该轨迹的方向选取数据窗口，该数据窗口的中心与双曲线同相轴的轨迹重合，并且在垂直双曲线同相轴轨迹的方向，数据窗口的长度选为地震子波长度，最后，采用该数据窗口所对应的多条双曲线轨迹来构建变顶点双曲线Radon变换算子L；

f、利用变顶点双曲线Radon变换算子L得到如下的数学模型：d＝Lm；

其中，d表示时间域的原始地震数据，m表示时间域的Radon域模型；

对Radon域模型施加L1范数最小化约束，同时，对拟合误差施加L1范数最小化约束，对应的L1-L1范数最小化优化问题为：

g、计算逆矩阵并计算矩阵利用乘数交替方向方法来求解上述的L1-L1范数最小化优化问题，得到Radon域模型的估计结果；

h、判断所有地震道集是否处理完毕；如果否，返回执行上述步骤b至步骤g；如果全部处理完毕，则输出处理结果。

所述步骤g中利用乘数交替方向方法求解L1-L1范数最小化优化问题的计算过程如下：

g1、在每一步迭代中，Radon域模型的计算结果为：

其中，和为中间变量，并且

g2、判断迭代次数n是否达到最大迭代次数N；如果没达到，继续执行步骤g3；如果达到，输出Radon域模型的估计结果m⁽ⁿ⁺¹⁾；

g3、利用距离算子计算中间变量和

$z_1^{(n+1)}={\mathrm{prox}}_{\chi }^1[{\mathrm{Lm}}^{(n+1)}-d+b_1^{\left(n\right)}],z_2^{(n+1)}={\mathrm{prox}}_{\phi }^1[m^{(n+1)}+b_2^{\left(n\right)}],$

其中，定义距离算子为：定义距离算子

v＝{v_i,j}，C＝max(|v_i,j|)，

g4、对中间变量和计算如下：

并返回到步骤g1。

本发明具有如下优点：

本发明方法利用地震同相轴的倾角及其导数来计算双曲线同相轴的顶点位置和速度，并利用聚类算法来确定最终的双曲线同相轴个数，降低所估计的同相轴参数个数；利用所计算的同相轴顶点位置和速度，采用开窗分割地震道的方法来构建变顶点双曲线Radon变换算子，消除Radon反变换产生的地震子波拉伸效应；然后在优化问题中引入模型空间(Radon域模型)和数据拟合项(Radon反变换残差)的L1范数最小化约束，来提高Radon域模型的稀疏性和对非高斯拟合误差的适应性；另外，本发明方法利用乘数交替方向方法(Alternating Direction Method of Multipliers)求解优化问题，乘数交替方向方法利用距离算子求解L1范数最小化问题，该方法在整个迭代过程中只需一次矩阵求逆，计算复杂度较低。

附图说明

图1为本发明中考虑地震子波拉伸效应的变顶点稀疏双曲线Radon变换方法的流程图；

图2为本发明中缺失近偏移距道的3个满足双曲线时距曲线关系的同相轴示意图；

图3为利用本发明方法估计的倾角示意图；

图4为利用本发明方法所估计倾角的导数示意图；

图5为本发明中聚类方法检测到的双曲线同相轴的个数示意图；

图6为本发明中Radon正变换后的结果示意图；

图7为本发明中Radon反变换后的结果示意图。

具体实施方式

下面结合附图以及具体实施方式对本发明作进一步详细说明：

考虑地震子波拉伸效应的变顶点稀疏双曲线Radon变换方法，包括如下步骤：

a、设置变量初始值：需要设置初始值的变量包括阻尼因子β、乘数交替方向方法的阈值χ和φ、乘数交替方向方法的最大迭代次数N；

b、从所有地震道集中选取某一个地震道集，利用复地震道分析方法估计地震同相轴的倾角p_r，并计算倾角导数p_rr；

c、利用下式计算双曲线同相轴的顶点位置(τ,x₀)和速度v：

x₀＝x_r-tp_rv²，其中，τ表示顶点位置中的时间位置，x₀表示顶点位置中的空间位置，t表示时间，x_r表示检波点的空间位置；

d、在利用步骤c得到双曲线同相轴的顶点位置(τ,x₀)和速度v后，利用聚类算法来检测双曲线同相轴的个数，确定最终的同相轴顶点位置和速度；

e、利用步骤d得到的同相轴顶点位置和速度，确定双曲线同相轴的轨迹，然后，利用开窗分割地震道的方法沿着垂直该轨迹的方向选取数据窗口，该数据窗口的中心与双曲线同相轴的轨迹重合，并且在垂直双曲线同相轴轨迹的方向，数据窗口的长度选为地震子波长度，最后，采用该数据窗口所对应的多条双曲线轨迹来构建变顶点双曲线Radon变换算子L，来消除Radon反变换重构地震数据时出现的地震子波拉伸效应；

f、利用变顶点双曲线Radon变换算子L得到如下的数学模型：d＝Lm；

其中，d表示时间域的原始地震数据，m表示时间域的Radon域模型；

为提高Radon域模型的稀疏性，对Radon域模型施加L1范数最小化约束。同时，为有效提高本发明方法对非高斯拟合误差的适应性，对拟合误差施加L1范数最小化约束。对应的L1-L1范数最小化优化问题为：其中，λ表示正则化因子，用来均衡对Radon域模型的L1范数最小化约束和对数据拟合项的L1范数最小化约束；

g、为了求解上式中给定的优化问题，首先计算逆矩阵其中，阻尼因子β用来保证矩阵求逆的稳定性，并计算矩阵然后利用乘数交替方向方法来求解上述的L1-L1范数最小化优化问题，得到Radon域模型的估计结果；

该步骤g中利用乘数交替方向方法求解L1-L1范数最小化优化问题的计算过程如下：

g1、在每一步迭代中，Radon域模型的计算结果为：

其中，和为中间变量，并且

g2、判断迭代次数n是否达到最大迭代次数N；如果没达到，继续执行步骤g3；如果达到，输出Radon域模型的估计结果m⁽ⁿ⁺¹⁾；

g3、利用距离算子计算中间变量和

$z_1^{(n+1)}={\mathrm{prox}}_{\chi }^1[{\mathrm{Lm}}^{(n+1)}-d+b_1^{\left(n\right)}],z_2^{(n+1)}={\mathrm{prox}}_{\phi }^1[m^{(n+1)}+b_2^{\left(n\right)}],$

其中，定义距离算子为：定义距离算子

v＝{v_i,j}，C＝max(|v_i,j|)，

g4、对中间变量和计算如下：

并返回到步骤g1；

h、判断所有地震道集是否处理完毕；如果否，返回执行上述步骤b至步骤g；如果全部处理完毕，则输出处理结果。

在仿真实验中，利用模型数据验证本发明方法的有效性：

图2为3个满足双曲线时距曲线关系的同相轴，并且缺失近偏移距，横坐标Trace Number表示道号，单位为米(m)，纵坐标Time Sample Number表示时间采样点个数，采样间距为2毫秒(ms)，道间距为4.5米(m)。在该例子中，利用本发明方法来对缺失的近偏移距进行外推，来验证本发明中Radon变换的有效性。图3为利用复地震道分析方法估计的局部倾角，图4为对应的倾角导数。利用倾角及其倾角导数进行同相轴顶点位置和速度估计，并利用聚类算法检测同相轴的个数，图5为所检测到的3个同相轴顶点位置。图6为利用所构建的Radon变换算子进行Radon正变换所得到的结果。图7为利用图6中的Radon变换结果进行反变换后得到的同相轴，可以看到近偏移距得到有效的重建。注意到本发明实例中的双曲线同相轴的顶点并非完全位于零偏移距处，本发明方法能有效重建顶点不在零偏移距处的同相轴。

当然，以上说明仅仅为本发明的较佳实施例，本发明并不限于列举上述实施例，应当说明的是，任何熟悉本领域的技术人员在本说明书的教导下，所做出的所有等同替代、明显变形形式，均落在本说明书的实质范围之内，理应受到本发明的保护。