λ-f域三维抛物Radon变换多次波压制方法

井洪亮 张少华 方云峰* 马光凯 张 茜 杨雪霖

(①东方地球物理公司物探技术研究中心，河北涿州 072751； ②东方地球物理公司，河北涿州 072750； ③东方地球物理公司研究院处理中心，河北涿州 072751)

0 引言

自20世纪70年代起，Radon变换引入地球物理勘探，经过几十年的发展，已逐步用于地震资料处理的诸多方面，如地震道重建、VSP上、下波场分离及多次波压制等[1-6]。尤其在多次波压制方面，Radon变换因其应用简便、操作灵活，受到业界的广泛关注。Hampson[7]首次利用τ-p变换压制多次波，但由于变换域能量聚焦不足，存在严重的能量“拖尾”现象，多次波压制效果并不理想； 随后，Kostov[8]、Sacchi[9]提出了高分辨率Radon变换，有效地解决了“剪刀状”能量发散问题，提高了多次波压制精度； Wood等[10]推出一种混合域高分辨率Radon变换，充分利用时间域的高分辨率和频率域可解耦的特点，进一步提高了Radon域内解的稀疏度和分辨率，但计算成本偏高； Abbad等[11]基于奇异值分解技术提出了频率域快速求解算法，即λ-f域Radon变换，通过引入变量λ消除变换算子与频率的关系，显著提高了计算效率。

中国学者也深入研究了Radon变换多次波压制方法。刘喜武等[12]利用稀疏约束共轭梯度算法实现高分辨率Radon变换，并成功用于多次波压制； 王维红等[13]提出了一种加权抛物Radon变换，利用变化的权系数提高抗假频能力； 熊登等[14]利用混合域Radon变换在时间域进行稀疏约束，提高了Radon域内能量团的分辨率； 石颖等[15-18]基于二维λ-f域Radon变换优化和改进地震道重建、多次波压制和波场分离等方法。近年来，在Radon变换的保幅处理方面也取得一定突破。薛亚茹等[19-22]基于正交多项式变换和常规Radon变换的不同特征，将二者结合形成高阶Radon变换，通过构建多阶Radon子区间，使每个子区间保留不同的振幅变化信息，从而具有AVO特性； 唐欢欢等[23]将高阶二维Radon变换成功拓展到三维，并取得很好的保幅效果。王亮亮等[24]实现了快速Radon变换算法； 张全等[25]提出了CPU-GPU异构并行算法，加速比达到约30倍，提高了Radon变换计算效率。

在实际资料处理中，二维抛物Radon变换压制多次波的方法已经相当成熟，但是没有考虑地震子波的空间传播特性，忽略了波场随空间曲率的变化信息，在处理宽方位三维地震数据时，往往多次波压制效果不佳，存在一定的多次波残余。为此，本文基于二维抛物Radon变换，借鉴Abbad等[11]的方法，深入研究了常规三维抛物Radon变换，并首次提出了λ-f域三维抛物Radon变换多次波压制方法。结果表明：常规三维抛物Radon变换多次波压制方法的计算效率过低，不利于大规模实际应用；λ-f域三维抛物Radon变换多次波压制方法重新定义了变换算子的构建方式，消除了算子对地震数据频率的依赖，使整个变换过程仅需计算一次变换算子及其广义逆算子，大大减少了矩阵求逆运算的次数，明显提高了计算效率。同时，根据一次波和多次波映射到λ-f域的空间分布特征，设计了三维椎体滤波器，使一次波能量投影到保留区，多次波能量投影到切除区，避免了人工定义切除区域的过程，更便于实际应用。因此，针对海量三维地震数据处理，本文所提方法可以极大地节省人工成本。

1 方法原理

1.1 常规三维抛物Radon变换原理

二维抛物Radon变换是沿着抛物线路径对地震数据积分、求和的过程。类似地，三维抛物Radon变换的物理意义是沿着抛物面对地震数据求和处理，其正、反变换的离散形式分别为

(1)

(2)

式中：d(t,x,y)为三维地震数据体，t为时间坐标，x、y分别为横、纵向炮检距；m(qx,qy,τ)为变换域的模型空间，qx和qy分别为横、纵向离散曲率参数，τ为变换域零炮检距截距时间。

对式(1)、式(2)两边进行傅里叶变换，分别得到频率域的正、反变换

(3)

(4)

通常情况下，应用Radon变换压制多次波，多是对正变换的模型空间进行处理，再反变换回时空域，而Radon域能量团的收敛性直接影响多次波去除效果。因此，本文着重推导了三维抛物Radon正变换的实现过程。根据式(4)可知，在频率域内，不同频率切片的三维地震数据具有解耦性，对某个频率，其矩阵简化形式为

(5)

式中：Lqx=e-ifqxx2为与x有关的算子；Lqy=e-ifqyy2为与y有关的算子。采用最小二乘方法，使式(5)的目标函数最小化，即

(6)

分别求取矩阵Lqx和Lqy的逆，可以得到Radon域数据M。因此在最小二乘意义下，式(5)的解为

(7)

上标“H”表示取共轭。为了防止奇异值对矩阵求逆过程产生异常影响，引入阻尼因子μ(通常取值为0.01～1.00)，主要用于增强矩阵求逆的稳定性。最终解的表达式为

(8)

1.2 λ-f域三维抛物Radon变换原理

针对常规三维抛物Radon变换计算效率偏低的问题，本文借鉴Abbad等[11]提出的改进二维Radon变换方法，重新定义变换算子的构建方式，引入两个数学运算变量λx=qxf和λy=qyf(其物理意义为曲率与频率的乘积)，可将常规三维抛物Radon变换的f-qx-qy域转换到一个全新的λx-λy-f域。该法继承了二维λ-f域Radon变换思路，创新实现了三维算法，故简称为λ-f域三维抛物Radon变换方法。

针对任意频率分量，将变量λx和λy引入变换公式，则改进三维抛物Radon正、反变换为

(9)

(10)

类似于常规三维抛物Radon变换方法，针对某个频率切片，将式(10)的矩阵形式简化为

(11)

式中Lλx、Lλy为重构的变换算子，具体形式为

式中：m和p分别为三维地震道集的纵测线和非纵测线的道数；n和k为λ-f域的三维模型空间的维度。通过对变换算子Lλx、Lλy重构，Lλx由炮检距xm和λxn(λx的离散值)计算获得，Lλy由炮检距yp和λyk(λy的离散值)计算获得，从而消除了对频率f的依赖。

采用最小二乘求解方式得到λ-f域三维抛物Radon变换的解

(12)

以模拟的三维地震记录为例，对比、分析常规三维抛物Radon正变换及λ-f域三维抛物Radon正变换的谱分布特征及聚焦性(图1)。利用常规三维抛物Radon变换将合成数据(图1a)变换到三维Radon域，并抽取四条同相轴对应最大能量团的垂向(图1b)和横向(图1c)切片。理论上，当截距时间和介质速度横向变化不大时，同相轴映射到变换域应呈“点”状分布，但受采集孔径限制，形成空间截断效应，能量团存在“拖尾”现象，能量泄露造成收敛不足(图1b、图1c)。对图1a进行λ-f域三维抛物Radon正变换，抽取变换域内不同位置的纵向切片(图1d)，由于地震数据中每条同相轴的曲率qx、qy基本固定，因此λx和λy与频率f成正比，四条抛物同相轴映射到λ-f域呈以不同曲率为斜率的“窄直线”状分布，且当计算频率极低时，λx和λy趋于0，即所有直线近似过原点。抽取λ-f域内不同频率的横向切片(图1e)发现，四条同相轴在λ-f域能量团较收敛，且随着频率增加，更容易区分能量团的相对位置关系，便于后续的多次波分离。

图1 合成炮记录Radon正变换结果

1.3 λ-f域离散采样条件

为了避免λ-f域三维抛物Radon变换出现假频而影响算法精度，需要将所有同相轴对应的Nλ(λx和λy的个数)值准确映射到变换域。若Nλ不足使同相轴能量映射不完整，造成能量泄露； 若Nλ过大会引入一定的随机噪声，降低资料信噪比。选取合适的Nλ值是必要的，决定该值的两个重要参数是采样间隔和扫描范围。

根据常规三维抛物Radon变换的曲率选取准则，由λx=qxf、λy=qyf可知，λ-f域重采样间隔Δλx和Δλy应满足

(13)

(14)

式中：xmax、xmin和ymax、ymin分别为横向和纵向最大、最小炮检距； Δqx、Δqy分别为横、纵向曲率采样间隔。

根据式(11)，重构的Radon算子Lλx、Lλy的相位误差必须满足

(15)

(16)

因此，根据qx、qy及地震数据最大频率fmax，可以确定参数λx和λy的扫描范围

(17)

(18)

由式(13)～式(18)可见，Δλx、Δλy及λx、λy的扫描范围仅与纵、横向最大、最小炮检距以及空间采样间隔Δx、Δy有关，而与频率无关。

三维抛物Radon变换采样间隔Δλx、Δλy和λx、λy的扫描范围类似于常规二维Radon变换(图2)[11]，即：在任意频率范围内，所有同相轴的曲率q的扫描范围相同(图2a)，对于低频同相轴而言，在有限的曲率范围内，频率过低可能造成采样不足，无法满足采样定理要求，出现低频同相轴能量泄露、聚焦不稳定的问题； 基于λ-f域抛物Radon变换方法，参数λ的采样间隔和扫描范围仅与炮检距、道间距有关，在一定的λ取值范围内，曲率和频率呈反比，即频率越低，曲率范围越大，低频区的覆盖区范围更宽，增加了低频段到变换域的完整映射几率(图2b)。

图2 参数变量扫描范围示意图

1.4 λ-f域三维椎体滤波压制多次波

在实际资料处理中，利用抛物Radon变换压制多次波首先要对地震数据做动校正(NMO)，然后根据一次波和多次波能量在变换域内分布位置的不同，通过手动方式定义切除区域，从而分离一次波和多次波能量团。对于海量三维地震数据而言，采用常规三维抛物Radon变换人为定义三维切除区，计算成本偏高。本文充分利用NMO后的一次波和多次波能量在λ-f域呈“窄直线”状分布特点，提出了一种椎体滤波方法，可以有效分离一次波和多次波能量，避免了人工定义切除区域的过程。实际Radon变换压制多次波的经验表明，为了满足浅层非多次波干扰区地震数据的保真性，通常在变换域内提取多次波能量，将其反变换回时空域，得到多次波模型数据，再采用自适应相减法从原始数据中去除多次波模型数据，以降低对有效波的伤害。基于此思路，在设计三维椎体滤波器时，以提取多次波能量为目标。设地震数据d(x,y,t)的λ-f域三维抛物Radon正变换为M(λx,λy,f)，定义椎体滤波器为G(λx,λy,f)，不同曲率的同相轴会映射到不同的椎体滤波区(图3)。假设一次波同相轴经NMO后均被校平，视速度v趋近于∞，由公式q≈1/(2t0v2)可知，一次波的曲率近似为0，则一次波能量从椎体顶点附近垂向分布到椎体内部的切除区域。地震数据中多次波同相轴经NMO后变为具有一定曲率(范围为qmin～qmax)的抛物面，多次波能量从顶点附近呈直线辐射至椎体保留区域，则

图3 三维数据体λ-f域椎体滤波区域示意图

(19)

(20)

2 数据实验

2.1 三维模型数据测试

为了检验λ-f域三维抛物Radon变换方法的效果，以三维正演速度模型(图4)为例进行多次波压制试算。图5展示了图4模型数据多次波压制效果。由图可见：①NMO三维模型数据的一次波同相轴均被校平，多次波同相轴仍为抛物形态(图5a)。②常规三维抛物Radon正变换抽取一阶自由表面多次波能量最强位置，当一次波和多次波的截距时间较小时，由于存在能量团“拖尾”现象，部分能量交织在一起，很难分离多次波能量，在定义切除线过程中会损伤一次波能量或残留多次波能量(图5b)。③λ-f域三维抛物Radon正变换结果的一次波能量垂直分布，而多次波能量具有一定斜率，因此根据空间位置可分离一次波和多次波(图5c)。④常规三维抛物Radon变换和λ-f域三维抛物Radon变换均能有效提取多次波模型数据，但是在近炮检距位置都残留一次波(图5d、图5e红色箭头处)。这是由于近炮检距位置下层的一次波与上层的一阶自由表面多次波的NMO时差较小，很难有效分离所致，这也是Radon类方法在近炮检距位置多次波压制效果不佳的原因。⑤常规三维抛物Radon变换多次波压制结果明显残留多次波，这是由于常规三维抛物Radon变换的能量团收敛性差，在提取多次波能量过程中，为了减少一次波能量损伤，出现部分多次波能量泄露(图5f白色箭头处)。⑥由于在λ-f域易于分离一次波和多次波能量，因此λ-f域三维抛物Radon变换较干净、彻底地压制了多次波(图5g)。

图4 三维正演速度模型

为了进一步分析常规三维抛物Radon变换和λ-f域三维抛物Radon变换的多次波压制精度，随机抽取图5a、图5f和图5g同一位置的地震道(图6)，可见，后者的多次波压制精度高于前者。

图5 图4数据多次波压制效果

图6 由图5a、图5f和图5g随机抽取的同一位置的地震道对比

表1展示了利用不同方法压制多次波的计算效率，可见在采用相同的计算平台和计算参数情况下，对同一个节点的单线程运算，常规三维抛物Radon变换压制多次波方法的计算耗时(不包括人为定义切除过程)是λ-f域三维抛物Radon变换的8倍以上。因此，针对海量地震数据去噪，后者的计算高效性更明显，具有实用价值。

表1 利用不同方法压制多次波的计算效率

理论上，二维抛物Radon变换压制多次波方法存储变换算子的矩阵维度不大，对于每个频率，重新计算变换算子的耗时较短。因此，λ-f域二维抛物Radon变换压制多次波的效率提升并不明显。常规三维抛物Radon变换压制多次波方法存储变换算子的维度大，对于每个频率切片，需要重复进行大量的矩阵运算，尤其是求取广义逆运算，会严重制约计算效率。因此，λ-f域三维抛物Radon变换压制多次波充分体现了计算效率优势。

2.2 实际资料测试

以海洋拖揽数据为例测试多次波压制效果。图7为λ-f域三维抛物Radon变换压制海洋单炮记录多次波效果。由图可见：①NMO单炮记录中广泛存在一阶和二阶自由表面多次波，而且能量级别较大，部分有效波同相轴被掩盖(图7a)。②λ-f域三维抛物Radon正变换结果的一次波能量在λx=0、λy=0区域附近几乎垂直分布，根据多次波的最大和最小视速度范围定义椎体滤波区域，可以分离有效波和多次波(图7b)。③λ-f域三维抛物Radon变换有效恢复了能量强、曲率大的多次波(图7c)。④经多次波压制，图7a的强能量多次波均得到衰减，同时显现了部分能量较弱的一次波(图7d红色箭头处)。

图7 λ-f域三维抛物Radon变换压制海洋单炮记录多次波效果

图8为三维共炮检距剖面多次波压制效果。由图可见：原始数据存在强能量多次波，且地层具有一定倾角(图8a箭头标注处)； 常规三维抛物Radon变换有效压制了大部分较平缓的一阶、二阶自由表面多次波，但在一些构造倾角较大位置仍残留多次波模型包含两层水平层状介质，利用有限差分正演方法记录了上层的一阶和二阶自由表面多次波。模拟观测系统为规则网格，其中，纵、横向各500个网格点，网格间距为10m(图8b)； 与图8b相比，λ-f域三维抛物Radon变换干净、彻底地压制了多次波，尤其在一次波和多次波同相轴较接近的位置(t=3.3s)，且未损伤有效信号(图8c)。在相同的计算环境和计算参数情况下，常规三维抛物Radon变换处理681炮地震数据的耗时为684min，而λ-f域三维抛物Radon变换的耗时为71min，即后者的加速比提高了约9倍。

然而，尽管λ-f域三维抛物Radon变换方法的计算精度和效率高于常规三维抛物Radon变换方法，但在图8的红色方框位置仍残留多次波。为了进一步验证，抽取该位置附近的叠前CMP道集速度谱(图9)。可见，两种方法均存在明显的低速能量团，尤其在4.0s和4.7s位置附近。分析表明，在浅层(2s位置)存在绕射点，进而诱发强能量绕射多次波，其同相轴呈近似双曲线形态，不满足抛物Radon变换的假设条件，因此较难彻底压制。综上所述，利用λ-f域三维抛物Radon变换多次波压制方法去除海洋地震资料表面多次波的效果较好。

图8 三维共炮检距剖面多次波压制效果

图9 叠前CMP道集速度谱

3 结束语

λ-f域三维抛物Radon变换是对常规三维 Radon变换的改进，本文将其成功应用于三维地震资料多次波压制，理论模型和实际数据试算表明方法具有实用价值，适合于海量地震数据去噪，获得以下认识。

(1)与常规三维 Radon变换方法相比，λ-f域三维抛物Radon变换方法压制表面多次波更彻底，同相轴在λ-f域呈“窄直线”状分布，缓解了能量团“拖尾”现象，在多次波切除过程中降低了空间截断效应的影响。

(2)λ-f域三维抛物Radon变换通过对变换算子进行重构，消除了变换算子对频率的依赖，有效减少了矩阵求逆运算次数，较常规三维 Radon变换方法的计算效率提高了约8倍以上。

(3)利用λ-f域三维抛物Radon变换方法压制多次波时，需要根据多次波同相轴的视速度确定椎体滤波器的母线长度，即视速度越大，滤波器的切除线越向有效信号能量靠近。为了避免损伤有效信号，建议实际应用时不易选取过大的最大视速度。

尚需指出，λ-f域三维抛物Radon变换方法压制多次波时仍然存在Radon类方法的弊端，即当一次波和多次波速度差异较小时，在近炮检距位置无法得到较理想的多次波压制效果。