应用主频带约束高精度 Radon 变换的多次波压制方法

摘要： Radon 变换是地震资料处理中常用的一种算法，通过特定路径求和实现地震数据插值、多次波压制和波场分离。然而，由于采集数据和算法限制，变换域的分辨率较低，提高Radon 域分辨率一直是研究的热点。最常用的提高分辨率方法是多次迭代重加权算法，通过迭代更新权重，将加权值聚焦到地震数据曲率上，但很难将加权值聚焦到真实曲率位置。文中提出了一种优化提高变换域分辨率的新方法，即在地震数据的主频带内求取一个加权矩阵，并将加权值聚焦至真实地震曲率位置。首先，计算地震数据的主频，并以主频为中心频率，向低频和高频方向各取一个范围作为约束频带； 对主频带内的地震数据采用低频约束策略，从低频到高频迭代计算加权矩阵； 从主频带的最后一个加权矩阵得到最终的加权矩阵，并应用于所有其他频率的计算。主频带数据的信噪比高、振幅强，因此，该方法更加稳定，得到的加权矩阵可以显著增强变换域的分辨率。此外，与其他迭代方法相比，它避免了每个频率加权矩阵的迭代过程，计算效率更高。合成数据和实际数据的测试结果证明了该方法的在多次波压制方面的有效性和优势。

关键词： Radon 变换，多次波压制，高分辨率，主频带，约束，效率

中图分类号：P631 文献标识码：A DOI：10. 13810/j. cnki. issn. 1000‐7210. 2024. 06. 011

0 引言

Radon 变换是一种广泛应用于医学CT 图像重建和地震数据处理的数学方法，在地震同相轴识别、多次波压制、地震数据重建、速度分析和波场分离等方面有着广泛的应用。在地震资料处理领域，Radon变换可分为线性变换、抛物线变换和双曲线变换。

20 世纪 70 年代，美国斯坦福大学地球物理小组对Radon 变换进行深入研究，并取得了重大进展，为其在地球物理领域的应用做出了重要贡献。

Hampson[1]将线性Radon 变换改进为抛物线Radon变换，并应用于多次波压制。经过正常时差校正后，一次波被校正为平直的同相轴，而多次波则由双曲线校正为抛物线，可用抛物线方程描述。由于变换域内多次波和一次波曲率不同，可以在变换域通过切除实现二者的分离。

Beylkin[2]讨论了基于最小二乘理论的Radon 变换。然而，由于地震数据采集的空间和时间有限，变换后的Radon 域的聚焦点经常出现发散假象，极大地影响了多次波压制的效果。因此，研发高分辨率 Radon变换可以改善同相轴的可分离性、消除发散假象。

提高 Radon 域的分辨率的途径可以大致分为如下三类。第一类算法是在时间域通过大尺度求逆运算提高分辨率。时域方法具有一定的优势，比如在数据和模型上都能处理时变问题，在消除混叠假象方面效果显著，通常能获得分辨率更高的结果。Thorson 等[3]给出了一种采用随机反演的高分辨率时域Radon 变换算法，同时增强了时间和曲率方向的稀疏性。Stoffa 等[4]提出了一种基于平面波分解、利用相似加权导出的窗滤波器，可以有效消除混叠假象，提高分辨率。Yilmaz 等[5]采用时域离散平面波分解技术计算每个点的斜率，迭代估计地震数据中的倾角范围，并应用到变换中，有效消除了混叠假象。

然而，时域高分辨率方法计算效率较低，很快被Hampson[1]的频率域快速、时不变抛物 Radon 变换所取代。但Hampson 的频域最小二乘方法由于时空域地震数据有限，在变换域存在截断效应，分辨率较低。为此，Sacchi 等[6]提出了一种基于贝叶斯框架的优化反演算法、利用最大熵函数计算先验概率、采用以柯西范数为目标函数约束项的非线性优化方法，并假设残差项和模型参数项分别服从高斯分布和柯西分布，利用迭代重加权最小二乘 （IRLS）法得到稀疏解。IRLS 法是频率域上常用的高分辨率Radon 变换算法之一，属于高分辨率算法的第二类。Sacchi 等[6]的方法在处理稀疏采样数据时会出现空间混叠现象，从而大大降低分辨率，且计算成本比传统的最小二乘法高。Cary[7]强调了时域算法在时间和曲率两个方向都具有稀疏性的优势，同时建议结合频域算法改善Radon 变换的保幅性。Herrmann等[8]提出了一种利用低频数据约束高频抛物分解的非迭代高分辨率方法，避免了空间混叠，同时实现了在空间域稀疏情况下的高分辨率结果。然而，该方法需先计算低频结果，并以此作为权重约束高频计算，计算量仍然较大。Chen 等[9]提出了一种利用峰值频率求取加权矩阵的非迭代方法。然而，该方法的结果受所选频率的影响较大，并且在变换域的分辨率不够高。所有频域方法在频域计算的模型权重在时间域都耦合在一起，在不同时间对所有同相轴应用相同的权重。因此，频率域高分辨率算法的分辨率只在 Radon 参数，即曲率方向有所提高。此外，频域高分辨率Radon 变换算法对反演过程正则化所需的超参数比较敏感，必须经过大量的测试才能确定。如果每个频率的超参数均是默认的，即使有一个频率的参数选择不当，也会影响整个Radon 域的结果，这使得该方法的实现同样具有挑战性。

第三类高分辨率 Radon 方法是时频域混合方法。该类算法计算效率高、分辨率高。时频混合域高分辨率算法在传统的频域最小二乘算法基础上，增加一个额外的时域阈值收缩步骤。这一步骤有助于缓解截断效应引起的边缘能量泄漏，使变换域能量向聚焦点中心收敛，从而提高 Radon 变换的分辨率，如Lu[10]利用混合域算法提出的稀疏Radon 变换，在时域通过简单的收缩处理提升变换域的稀疏性，且在频域实现正、反Radon 变换，兼顾了时域和频域算法的优势。Wang 等[11]在时频混合域提出了鲁棒时不变的 Radon 变换，并在时域对数据不拟合和 Radon 模型进行双L1 范数稀疏约束压制异常值，在频域进行正、反变换，利用一维交替分裂Bregman算法对Radon 模型进行阈值收缩。Li 等[12]提出了一种快速迭代收缩阈值算法，并将其应用于自适应多次波压制。Liu 等[13]将基于L1 范数约束的高分辨率Radon 变换与形态成分分析方法相结合，对地震数据中的多次波和折射波进行分离。Geng等[14‐15]利用 L1‐2 范数约束在时频混合域提出了稀疏Radon 变换，该变换比 L1 范数约束具有更强的稀疏性，并采用三维拉格朗日乘子进一步将该算法扩展到三维空间。薛亚茹等[16] 将IRLS 和阈值收缩（ISTA）法相结合，给出了一种加权迭代软阈值方法。与频域高分辨率算法类似，这种混合域高分辨率算法面临时域收缩参数设置困难的问题，收缩参数设置不合理会导致算法不收敛。张全等[17]将全局贪婪优化算法同迭代阈值收缩结合，引入Radon 反问题求解，提高了收敛速率，但需要设置迭代步长和算法的重启条件，复杂度较高。井洪亮等[18] 通过将传统Radon 变换的离散曲率参数q 更改为耦合频率的参数λ，变换到另外的λ⁃f 模型域，并设计三维体锥形滤波器对多次波进行切除，减少了计算量，但滤波器形式复杂，对于近炮检距处小时差的多次波压制效果差。Li 等[19]提出了一种结合了频率—慢度域Radon 算子的鲁棒正交匹配追踪算法，使用交替乘子法求解，用L1 ‐L1 范数稀疏约束增强了Radon 变换的群稀疏性，并应用于混合震源数据分离和三维地震数据重建。

基于前人的研究，本文给出了一种利用主频带约束的高分辨率 Radon 变换算法。通过将单频扩展到一个主频带，可得到一个稀疏约束矩阵； 与单频相比，该加权约束矩阵能更精确地描述地震数据的特征，将地震数据的能量聚焦于主要曲率参数。利用地震数据的主频带计算稀疏加权矩阵，并约束求解所有频率，可以提高计算精度和效率。在该算法中，所有频率的加权矩阵是相同的。在该主频带内，采用Hermann 等[8]的高分辨率策略由低频到高频更新加权矩阵。主频带内的数据振幅强、信噪比高，因此能对地震数据进行精准描述。该算法参数设置也非常简单，可以快速生成算法所需的加权矩阵。与Chen 等[9]的峰值频率约束的高分辨率方法相比，本方法更加稳定，分辨率更高； 与 IRLS 算法相比，本方法具有较高的计算效率和较强的抗噪性。分别使用二维合成数据和实际数据验证了本文算法的优势。