考虑地震子波拉伸效应的变顶点稀疏双曲线Radon变换方法
【专利摘要】本发明公开了一种考虑地震子波拉伸效应的变顶点稀疏双曲线Radon变换方法。该方法利用地震同相轴的倾角及其导数来计算双曲线同相轴的顶点位置和速度，并利用聚类算法来确定最终的双曲线同相轴个数，降低所估计的同相轴参数个数；利用所计算的同相轴顶点位置和速度，采用开窗分割地震道的方法来构建变顶点双曲线Radon变换算子，消除Radon反变换产生的地震子波拉伸效应；然后在优化问题中引入Radon域模型和数据拟合项的L1范数最小化约束，来提高Radon域模型的稀疏性和对非高斯拟合误差的适应性；另外，本发明方法利用乘数交替方向方法求解L1?L1范数最小化优化问题，在整个迭代过程中只需一次矩阵求逆，计算复杂度较低。
【专利说明】
考虑地震子波拉伸效应的变顶点稀疏双曲线Radon变换方法
技术领域
[0001] 本发明属于地震勘探技术中的地震信号处理领域，具体涉及一种考虑地震子波拉 伸效应的变顶点稀疏双曲线Radon变换方法。
【背景技术】
[0002] Radon变换广泛应用于地震信号处理中，如多次波压制、地震数据规则化和波场分 离等，其主要思想是将地震数据从时间-偏移距域映射到Radon域。Radon变换按照积分轨迹 区分为三大类:线性Radon变换、抛物线Radon变换和双曲线Radon变换。
[0003] 其中，双曲线Radon变换可直接应用于未做动校正的共中心点道集或共炮点道集， 并利用双曲线轨迹来构造变换算子，能够更好地表示地震数据。
[0004] 由于双曲线Radon变换的时变特性，计算双曲线Radon变换需要对大型稀疏双曲线 Radon变换算子进行矩阵求逆，导致时间域双曲线Radon变换过程中的巨大计算量。另外，传 统的双曲线Radon变换将顶点位置限制在零偏移距处，利用不同的速度参数来表征地震道 集中的双曲线轨迹。当双曲线同相轴的顶点位于非零偏移距位置，如道集中存在微曲多次 波、绕射多次波，传统的双曲线Radon变换进行多次波压制或数据规则化时，其效果会变差。
[0005] 双曲线Radon变换的应用效果依赖于变换算子是否能有效表征道集中存在的双曲 线同相轴。对于双曲线同相轴的顶点沿偏移距移动的情况，变顶点双曲线Radon变换方法通 过将Radon域模型的维度空间从2D (顶点时间和速度)扩展为3D (顶点时间、速度和顶点偏移 距），能够更有效地表示顶点位于非零偏移距的双曲线同相轴。由于Radon域模型的空间维 度增加一维，导致在Radon域模型中需要求解的未知参数个数变多，增加了Radon正反变换 的计算量。另外，变顶点双曲线Radon变换方法对应的Radon变换方程组为欠定方程组，需要 在相应的优化问题中引入Radon域模型的稀疏性约束来计算Radon域模型。传统的变顶点稀 疏双曲线Radon变换方法采用迭代重加权最小二乘算法求解L1范数最小化优化问题，在每 一步迭代中需要计算一次矩阵求逆，计算复杂度较高。因此，传统变顶点稀疏双曲线Radon 变换的计算量较大，这限制了该变换方法的实用性。另外，传统的双曲线Radon变换在构建 双曲线Radon变换算子进行正反变换时，并没有考虑地震子波的存在，在利用Radon反变换 进行地震数据重构过程中，会出现地震子波拉伸效应，造成地震信号失真。
【发明内容】
[0006] 针对现有技术中存在的上述技术问题，本发明提出了一种考虑地震子波拉伸效应 的变顶点稀疏双曲线Radon变换方法，其采用如下技术方案：
[0007] 考虑地震子波拉伸效应的变顶点稀疏双曲线Radon变换方法，包括如下步骤：
[0008] a、设置变量初始值:需要设置初始值的变量包括阻尼因子0、乘数交替方向方法的 阈值x和巾、乘数交替方向方法的最大迭代次数N;
[0009] b、从所有地震道集中选取某一个地震道集，利用复地震道分析方法估计地震同相 轴的倾角Pr，并计算倾角导数Prr;
[0010] C、利用下式计算双曲线同相轴的顶点位置(t，xq)和速度V:
，其中，T表示顶点位置中的时间位 置，XQ表示顶点位置中的空间位置，t表示时间，Xi?表示检波点的空间位置；
[0012] d、在利用步骤c得到双曲线同相轴的顶点位置和速度后，利用聚类算法来检测双 曲线同相轴的个数，确定最终的同相轴顶点位置和速度；
[0013] e、利用步骤d得到的同相轴顶点位置和速度，确定双曲线同相轴的轨迹，然后，利 用开窗分割地震道的方法沿着垂直该轨迹的方向选取数据窗口，该数据窗口的中心与双曲 线同相轴的轨迹重合，并且在垂直双曲线同相轴轨迹的方向，数据窗口的长度选为地震子 波长度，最后，采用该数据窗口所对应的多条双曲线轨迹来构建变顶点双曲线Radon变换算 子L;
[0014] f、利用变顶点双曲线Radon变换算子L得到如下的数学模型:d = Lm;
[0015]其中，d表示时间域的原始地震数据，m表示时间域的Radoni或模型；
[0016]对Radon域模型施加 L1范数最小化约束，同时，对拟合误差施加 L1范数最小化约 束，对应的L卜L1范数最小化优化问题为：argnjjn []|d -Lml + 4 ImjL];
[0017] g、计算逆矩阵丨::= (L〃L + /?丨)1 ;并计算矩阵t = El/ ;利用乘数交替方向方法来求 解上述的L1-L1范数最小化优化问题，得到Radon域模型的估计结果；
[0018] h、判断所有地震道集是否处理完毕；如果否，返回执行上述步骤b至步骤g;如果全 部处理完毕，则输出处理结果。
[0019] 所述步骤g中利用乘数交替方向方法求解L1-L1范数最小化优化问题的计算过程 如下：
[0020] gl、在每一步迭代中，Radoni或模型的计算结果为= ￡卜(广-+ -<)，
[0021] 其中，zf、b!Bi、zf和bf为中间变量，并且#=0:，bf=0，b?=0;
[0022] g2、判断迭代次数n是否达到最大迭代次数N;如果没达到，继续执行步骤g3;如果 达到，输出Radon域模型的估计结果m (n+1);
[0023] g3、利用距离算子计算中间变:
[0024] z(；+n = prox\ [Lm(n+1) -d + , zfl) = prox\ [m(n+,) + b{：]],
[0025] 其中，定义距离算子 为：/7/yav:上'],=n.mx (卜jC,0).s'_ 算子[ v];max (卜,卜 d 0,卜gn (v,.,); r 丨-1，a > 0 r , , , ,、/、 la, a> 0 , ,x
[0026] v= {vi，j}，C=max(卜i，j | )，(") = ，sgii^/)二）O' :a = Q ; 1 [-1, a<0
[0027] g4、对中间变量b(/〃U和b(_:H 15计算如下：
[0028] bf+1)=事 +Lm丨,,+1丨-d -zf,+1)，時+1)=球): + m(n+11 -4叫^
[0029] 本发明具有如下优点：
[0030] 本发明方法利用地震同相轴的倾角及其导数来计算双曲线同相轴的顶点位置和 速度，并利用聚类算法来确定最终的双曲线同相轴个数，降低所估计的同相轴参数个数;利 用所计算的同相轴顶点位置和速度，采用开窗分割地震道的方法来构建变顶点双曲线 Radon变换算子，消除Radon反变换产生的地震子波拉伸效应;然后在优化问题中引入模型 空间（Radon域模型）和数据拟合项(Radon反变换残差）的L1范数最小化约束，来提高Radon 域模型的稀疏性和对非高斯拟合误差的适应性；另外，本发明方法利用乘数交替方向方法 (Alternating Direction Method of Multipliers)求解优化问题，乘数交替方向方法利 用距离算子求解L1范数最小化问题，该方法在整个迭代过程中只需一次矩阵求逆，计算复 杂度较低。
【附图说明】
[0031]图1为本发明中考虑地震子波拉伸效应的变顶点稀疏双曲线Radon变换方法的流 程图；
[0032]图2为本发明中缺失近偏移距道的3个满足双曲线时距曲线关系的同相轴示意图；
[0033] 图3为利用本发明方法估计的倾角示意图；
[0034] 图4为利用本发明方法所估计倾角的导数示意图；
[0035] 图5为本发明中聚类方法检测到的双曲线同相轴的个数示意图；
[0036]图6为本发明中Radon正变换后的结果示意图；
[0037]图7为本发明中Radon反变换后的结果示意图。
【具体实施方式】
[0038]下面结合附图以及【具体实施方式】对本发明作进一步详细说明：
[0039]考虑地震子波拉伸效应的变顶点稀疏双曲线Radon变换方法，包括如下步骤：
[0040] a、设置变量初始值:需要设置初始值的变量包括阻尼因子0、乘数交替方向方法的 阈值x和巾、乘数交替方向方法的最大迭代次数N;
[0041 ] b、从所有地震道集中选取某一个地震道集，利用复地震道分析方法估计地震同相 轴的倾角Pr，并计算倾角导数Prr;
[0042] c、利用下式计算双曲线同相轴的顶点位置(t，XQ)和速度v:
，其中，T表示顶点位置中的时间位 置，XQ表示顶点位置中的空间位置，t表示时间，Xi?表示检波点的空间位置；
[0044] d、在利用步骤c得到双曲线同相轴的顶点位置(T，XQ)和速度v后，利用聚类算法来 检测双曲线同相轴的个数，确定最终的同相轴顶点位置和速度；
[0045] e、利用步骤d得到的同相轴顶点位置和速度，确定双曲线同相轴的轨迹，然后，利 用开窗分割地震道的方法沿着垂直该轨迹的方向选取数据窗口，该数据窗口的中心与双曲 线同相轴的轨迹重合，并且在垂直双曲线同相轴轨迹的方向，数据窗口的长度选为地震子 波长度，最后，采用该数据窗口所对应的多条双曲线轨迹来构建变顶点双曲线Radon变换算 子L，来消除Radon反变换重构地震数据时出现的地震子波拉伸效应；
[0046] f、利用变顶点双曲线Radon变换算子L得到如下的数学模型:d = Lm;
[0047] 其中，d表示时间域的原始地震数据，m表示时间域的Radon域模型；
[0048]为提高Radon域模型的稀疏性，对Radon域模型施加 L1范数最小化约束。同时，为有 效提高本发明方法对非高斯拟合误差的适应性，对拟合误差施加 L1范数最小化约束。对应 的L1-L1范数最小化优化问题为：argmjnQd-Lml +判my ;其中，A表示正则化因子，用来 均衡对Radon域模型的L1范数最小化约束和对数据拟合项的L1范数最小化约束；
[0049] g、为了求解上式中给定的优化问题，首先计算逆矩阵t = (L〃L + A丨)1 ;其中，阻尼 因子0用来保证矩阵求逆的稳定性，并计算矩阵：L = El/;然后利用乘数交替方向方法来求 解上述的L1-L1范数最小化优化问题，得到Radon域模型的估计结果；
[0050] 该步骤g中利用乘数交替方向方法求解L1-L1范数最小化优化问题的计算过程如 下：
[0051 ] gl、在每一步迭代中，Radoni或模型的计算结果为-b丨+ -bf)，
[0052] 其中，<、bf、和峻)为中间变量，并且#=0，1>丨°)=0，4°丨=0，<)=0;
[0053] g2、判断迭代次数n是否达到最大迭代次数N;如果没达到，继续执行步骤g3;如果 达到，输出Radon域模型的估计结果m (n+1);
[0054] g3、利用距离算子计算中间变量<+1>和<+1> :
[0055] zj,!+1) Lm(,,+1) - d + b[n| > = prox^ + b^ >
[0056] 其中，定义距离算子力：[v]; ; = max&. ,|-^fprox[\\\ i -max(|i； ,.|-^C,0)sgn(v.,.)； r f 1, a > 0 . . ,x a, a>0 ,、 . .
[0057] v={vij} ?C=max( | Vij | ) ? (^)：I = ，sgn(")二、a " 二 0 ; U, a = U 、 -1, a<Q
[0058] g4、对中间变量和bf0计算如下：
[0059] 丨+Lm(,I+l) -d-z(;+1 丨，<+1)=甿丨 +m!B+11 -4n+1)，并返回到步骤 gl;
[0060] h、判断所有地震道集是否处理完毕；如果否，返回执行上述步骤b至步骤g;如果全 部处理完毕，则输出处理结果。
[0061 ]在仿真实验中，利用模型数据验证本发明方法的有效性：
[0062]图2为3个满足双曲线时距曲线关系的同相轴，并且缺失近偏移距，横坐标Trace Number表示道号，单位为米(m)，纵坐标Time Sample Number表示时间采样点个数，采样间 距为2毫秒(ms)，道间距为4.5米(m)。在该例子中，利用本发明方法来对缺失的近偏移距进 行外推，来验证本发明中Radon变换的有效性。图3为利用复地震道分析方法估计的局部倾 角，图4为对应的倾角导数。利用倾角及其倾角导数进行同相轴顶点位置和速度估计，并利 用聚类算法检测同相轴的个数，图5为所检测到的3个同相轴顶点位置。图6为利用所构建的 Radon变换算子进行Radon正变换所得到的结果。图7为利用图6中的Radon变换结果进行反 变换后得到的同相轴，可以看到近偏移距得到有效的重建。注意到本发明实例中的双曲线 同相轴的顶点并非完全位于零偏移距处，本发明方法能有效重建顶点不在零偏移距处的同 相轴。
[0063]当然，以上说明仅仅为本发明的较佳实施例，本发明并不限于列举上述实施例，应 当说明的是，任何熟悉本领域的技术人员在本说明书的教导下，所做出的所有等同替代、明 显变形形式，均落在本说明书的实质范围之内，理应受到本发明的保护。
【主权项】
1. 考虑地震子波拉伸效应的变顶点稀疏双曲线Radon变换方法，其特征在于，包括如下 步骤： a、 设置变量初始值:需要设置初始值的变量包括阻尼因子β、乘数交替方向方法的阈值 X和Φ、乘数交替方向方法的最大迭代次数Ν; b、 从所有地震道集中选取某一个地震道集，利用复地震道分析方法估计地震同相轴的 倾角Pr，并计算倾角导数Prr; C、利用下式计算双曲线同相轴的顶点位置(τ，Χ〇)和速度v:其中，τ表示顶点位置中的时间位置，Χ0 表示顶点位置中的空间位置，t表示时间，xr表示检波点的空间位置； d、 在利用步骤C得到双曲线同相轴的顶点位置和速度后，利用聚类算法来检测双曲线 同相轴的个数，确定最终的同相轴顶点位置和速度； e、 利用步骤d得到的同相轴顶点位置和速度，确定双曲线同相轴的轨迹，然后，利用开 窗分割地震道的方法沿着垂直该轨迹的方向选取数据窗口，该数据窗口的中心与双曲线同 相轴的轨迹重合，并且在垂直双曲线同相轴轨迹的方向，数据窗口的长度选为地震子波长 度，最后，采用该数据窗口所对应的多条双曲线轨迹来构建变顶点双曲线Radon变换算子L; f、 利用变顶点双曲线Radon变换算子L得到如下的数学模型:d = Lm; 其中，d表示时间域的原始地震数据，m表示时间域的Radon域模型； 对Radon域模型施加 L1范数最小化约束，同时，对拟合误差施加 L1范数最小化约束，对 应的LI -L1范数最小化优化问题为：g、 计算逆矩阵Γ = (i/l1并计算矩阵；￡ ;利用乘数交替方向方法来求解上 述的L1-L1范数最小化优化问题，得到Radon域模型的估计结果； h、 判断所有地震道集是否处理完毕;如果否，返回执行上述步骤b至步骤g;如果全部处 理完毕，则输出处理结果。2. 根据权利要求1所述的考虑地震子波拉伸效应的变顶点稀疏双曲线Radon变换方法， 其特征在于， 所述步骤g中利用乘数交替方向方法求解L1-L1范数最小化优化问题的计算过程如下： gl、在每一步迭代中，Radoni或模型的计算结果为= + E卜t -bj， 其中，zf、lf)、者)和W为中间变量，并且#=0，bf):二0，#=0，bf =0; g2、判断迭代次数η是否达到最大迭代次数N;如果没达到，继续执行步骤g3;如果达到， 输出Radon域模型的估计结果m(n+1); g3、利用距离算子计算中间变量zf4和ζ?Γ11:其中，定义距离算子严为：v = 卜,| - /C', 0) sgn (v,,)，定义距离算子
【文档编号】G01V1/32GK105929447SQ201610377985
【公开日】2016年9月7日
【申请日】2016年5月31日
【发明人】李钟晓, 李振春
【申请人】中国石油大学(华东)