基于低频约束的高分辨率Radon变换多次波压制方法研究

2019-07-11刘仕友马继涛孙万元应明雄

物探化探计算技术 2019年3期

刘仕友，马继涛，孙万元，应明雄

(1. 中海石油(中国)有限公司湛江分公司，湛江 524057；2. 中国石油大学(北京) 地球物理学院物探系，北京 102249)

0 引言

多次波是海洋地震数据中常见的噪音，有效去除多次波是地震资料处理中的一个重要研究内容。Radon变换是工业界多次波压制最为常用的有效手段之一，同时也被广泛用于地震数据重建，波场分离，速度分析等。

奥地利数学家Radon在1917年提出Radon变换，美国斯坦福大学地球物理小组在20世纪70年代对Radon变换进行了重点研究，并取得了显著进展，为其在地球物理领域的应用做出了重要贡献。Hampson[1]将线性Radon变换改进为抛物线Radon变换，并将其应用于多次波的压制中；抛物线Radon变换利用一次波和多次波动校正速度的差异进行多次波压制。在动校正后，地震数据的一次波基本被校平，多次波同相轴则由双曲线被部分动校为抛物线，通过抛物线Radon变换可以把时空域中的曲线映射为Radon域中的点，多次波和一次波对应曲线曲率具有较大的差异，从而在Radon域可以实现一次波和多次波的分离，达到去除多次波的目的。Beylkin[2]对基于最小二乘理论的Radon变换方法进行了讨论，由于Radon变换算法自身和地震数据采集空间和时间有限等原因，Radon域的数据往往会出现能量发散的情况，这严重影响了该方法多次波压制的效果，因此发展高分辨率Radon变换算法是十分必要的。Sacchi等[3]提出高精度的Radon变换，奠定了高精度Radon变换的理论基础。迭代高分辨率算法在处理稀疏采样数据时分辨率降低，其原因在于有限的空间采样会限制Radon变换的分辨率，产生空间假频，同时该算法进行了多次迭代运算，计算量较大。Herrmann[4]提出一种不需要迭代高分辨率算法，用数据的低频部分来约束数据的抛物线分解，低频结果在运算中不会产生假频，能够在采样稀疏的情况下仍取得比较好的高分辨率效果。

国内的专家学者对高分辨率Radon变换也展开过较为深入详细的研究。吴律[5]系统阐述了Tau-p的原理及应用情况；牛滨华[6]等提出了基于多项式理论的Radon变换；许多学者对Radon变换及其提高分辨率的方法做了系统深入地研究[7-13]。

笔者对抛物线Radon变换、基于迭代的高分辨率的Radon变换，及基于低频约束的高分辨率Radon变换方法进行系统阐述，并对比了这些方法在地震数据多次波压制中的效果。模拟数据和实际数据的对比可以看出，相比于普通最小二乘Radon变换和迭代高分辨率算法，基于低频约束的高分辨率Radon变换算法分辨率更高，多次波压制更加彻底。

1 方法原理

1.1 Radon变换原理

二维抛物线Radon变换是沿着抛物线对地震数据进行叠加得到Radon域中的数据，二维Radon变换可表示为：

(1)

其中:m为Radon域数据；d为时空域地震数据；t、τ为时间(一个是Radon域的时间，一个是时空域的时间)；q为描述曲线弯曲程度的曲线参数；x为偏移距。对式(1)进行离散化，可得式(2)。

(2)

其中：Nq为Radon域数据的曲率个数。对两边做傅立叶变换，转换到频率域中进行计算，

(3)

将式(3)写为算子形式，则抛物线Radon变换可以表示为：

D=LM

(4)

其中，算子L可以写为：

(5)

因此，为得到Radon域的数据M，可以建立如式(6)目标函数进行求解。

(6)

将式(6)最小化，可以得到M的最小二乘解，即：

(7)

其中,μ是阻尼因子，取值一般在0.1～1之间。由式(7)可以看出，在进行Radon域M的计算时，对所有频率μ均为相同的定值，其目的是为了提高矩阵运算的稳定性，但该因子同时也对Radon域的数据进行了平滑，降低了该域数据的分辨率。

1.2 基于迭代的高分辨率Radon变换

因为常规Radon变换在Radon域数据能量发散，影响多次波压制的效果，因此Sacchi等发展了基于迭代的高分辨率Radon变换，该方法建立了如下的目标函数：

(8)

其中，R(M)是施加在Radon域数据M上的规则化因子。对式(8)最小化，可得到Radon域数据的计算公式如下：

M(ω)=(L(ω)HL(ω)+μQ(M))-1L(ω)HD(ω)

(9)

其中，施加在矩阵LHL上的规则化因子与Radon域数据M存在如下的关系：

(10)

即某个频率的加权矩阵是通过某规则，由Radon域的计算结果M给出的；如果Radon域计算结果M中的某部分能量较强，则加权矩阵Q较小，如果Radon域计算结果M中的某部分能量较弱，则加权矩阵较大，这样通过式(9)中的矩阵求逆过程，会使得Radon域能量强的点更强，能量弱的点减弱，从而使得Radon域能量聚焦，达到高分辨率的效果。该算法是一种基于贝叶斯原理的空间稀疏约束算法，通过迭代运算对Radon域内的能量进行聚焦，实际计算时，可采取如下的迭代公式进行运算：

Mk(ω)=(L(ω)HL(ω)+

λQ(Mk-1(ω)))-1L(ω)HD(ω)

(11)

其中，k为迭代次数。可以看出，某次迭代的加权矩阵Q，是由上一次迭代的运算结果M给出的。加权矩阵可以根据上一次迭代结果的Radon域数据M的大小直接赋值，例如，

(12)

其中，ε为使得求倒数稳定所施加的一个因子，可以通过式(13)进行计算。

ε=max(0.01*max|M|,1×10-6)

(13)

该算法对每个频率都需要迭代进行，运算量较大；如果想要取得较为满意的结果，一般每个频率都需要迭代数十次。如果处理的数据体很大，需要耗费大量的运算时间。

1.3 基于低频约束的高分辨率Radon 变换

由前面可知，式(9)中的加权因子Q，可以由上一次迭代的计算结果给出。Herrmann对公式(9)进行了改进，由上一个频率的计算结果计算加权因子Q[4]。该算法的一个优势是避免了每个频率为得到加权矩阵进行的迭代运算，对于每个频率，直接代入由上个频率计算结果得到的加权矩阵，就可以得到该频率的计算结果。另外，由于使用了低频的计算结果对高频的计算进行了约束，该算法可以抗假频。

非迭代的，基于低频约束的高分辨率Radon变换可表示为：

M(ωn)=(L(ωn)HL(ωn)+

μQ(M(ωn-1)))-1L(ωn)HD(ωn)

(14)

其中：Q为对角加权矩阵，由上一频率的计算结果得到；角标n表示第n个频率；式中的阻尼因子μ的取值范围变化较大，可以通过测试得到，一般的取值范围为0.01～1；Q可表示为：

Qii(ωn)=‖Mi(ωn-1)‖,i=1,2…Nq

(15)

其中，Mi为在上一频率运算得到计算结果。从式(14)、式(15)中可以看出，某个频率的加权矩阵，是由上一个频率的计算结果给出的。在没有多次波曲率先验信息的情况下，该方法更看重低频结果对整体结果的约束，低频结果在运算中不会产生假频，能够在采样稀疏的情况下仍取得比较好的高分辨率效果。

基于低频约束的高分辨率Radon变换进行多次波压制的处理步骤如下：

1)通过傅立叶变换将地震数据转换到f-x域。

2)给初始对角加权矩阵赋值为单位矩阵。

3)对低频数据进行处理，将处理结果在矩阵中保存。

4)利用上一频率的数据计算结果计算新的加权矩阵，用来进行下一个频率数据的计算，直至所有频率都计算完成。

图1 模拟数据Fig.1 Simulated data

图2 常规最小二乘算法结果Fig.2 The results of regular least-squares result(a)Radon变换结果;(b)多次波压制后的结果

5)反傅立叶变换，并对计算结果进行一次波切除。

6)进行Radon反变换，完成相减运算，得到多次波压制结果。

2 数据测试

2.1 模拟数据测试

利用一个模拟数据对三种Radon变换方法的分辨率和多次波压制效果进行了对比。首先生成一个模拟数据，该模拟数据由两个一次波和两个对应的多次波组成，一次波分别位于0.4 s、0.9 s，多次波分别位于0.8 s、1.2 s，两个多次波的动校时差量分别为330 ms、480 ms,该模拟数据如图1所示。

利用基于最小二乘的抛物线Radon变换对该数据进行多次波压制，在进行多次波切除时，选择的Radon域切除的q值为50 ms，以减小对一次波的损害。该方法在Radon域的变换结果，以及多次波压制后的结果如图2所示。

图3 迭代方法提高分辨率结果Fig.3 The high resolution results based iterative method(a)Radon域结果；(b)多次波压制后的结果

从图2中可以明显看出，基于最小二乘理论的抛物线Radon变换在变换域具有明显的剪刀状发散现象，由于剪刀状能量发散，导致其压制结果具有明显的多次波残余；另外在近偏移距处还产生了较为明显的假象。

我们利用基于迭代的高分辨率Radon变换对该模拟数据进行了多次波压制，其结果如图3所示。

从图3中可以看出，基于迭代的高分辨率Radon变换在Radon域，没有了剪刀状发散现象，分辨率有了明显提高，但该结果在一次波下方的1.4 s处，以及多次波附近的0.9 s和1.2 s处，出现了假象。我们将迭代次数分别选取为10、20、40次，阻尼因子分别选取为0.01、0.05、0.1，该假象仍然存在。因此可以认为该假象有可能是由剪刀状发散所引起的，因为该迭代提高分辨率的方法，初始迭代值是最小二乘Radon变换的结果，该方法在进行迭代时，误将剪刀状发散的能量认为是有效能量进行进一步的聚焦，从而产生了本不存在的能量假象。该方法多次波压制过程中采用了和最小二乘Radon相同的切除参数，结果如图3(b)所示，可以看出，压制效果有了明显改善，但远偏移距的多次波仍有一些残余。

采用基于低频约束的抛物线Radon变换对该数据进行了多次波压制处理(图4)。从图4可以看出，基于低频约束的抛物线Radon变换变换域的分辨率很高，而且不存在假象；多次波压制后的结果优于其他两种方法的结果。

图4 低频约束提高分辨率结果Fig.4 The high results based on lower frequency constraint(a) Radon域结果；(b)多次波压制结果

图5 实际数据炮集Fig.5 The real shot gather(a)带有多次波的实际数据；(b)局部放大图

在计算时间上，对于此模拟数据，基于最小二乘的抛物线Radon变换和基于低频约束的高分辨率Radon变换分别耗时0.40 s和0.51 s，基于迭代的高分辨率Radon变换运算时间和所选取的迭代次数有关，本文运算中选取迭代次数与运算时间的关系如表1所示，从表1可以看出，基于迭代的高分辨率Radon变换运算时间在迭代次数是1时，运算时间和其他两种方法相比略大，但随着迭代次数增加，运算所需时间也随着大大增加。对比可以看出，基于低频约束的高分辨率Radon变换计算速度快，计算结果最优。

表1 基于迭代的高分辨率Radon变换运算时间和迭代次数表Tab.1 The iteration number and computation time of iterative based high resolution method

2.2 实际数据测试

利用墨西哥湾的一个CDP道集对三种抛物线Radon变换算法进行了测试(图5)，从图5中可以看出，一次波已经被校平，未校平的多次波在4 s以下。为更好地对结果进行对比，选取了该道集的3.2 s～4.8 s处进行了放大。

图6 实际数据最小二乘结果Fig.6 The least-squares results of real data(a)Radon域结果；(b)多次波压制结果

图7 实际数据迭代方法提高分辨率结果Fig.7 The real data high resolution results based iterative method(a)Radon域结果；(b)多次波压制结果

图8 实际数据低频约束提高分辨率结果Fig.8 The real data high results based on lower frequency constraint(a)Radon域结果；(b)多次波压制结果

我们分别利用三种抛物线Radon变换方法对CDP道集进行了处理(图6、图7、图8)。从图6～图8中可以看出，基于低频约束的抛物线Radon变换在变换域具有更高的分辨率，基于最小二乘理论的抛物线Radon变换，和其他两种方法相比，其压制结果在中远偏移距处有少量的多次波残余(箭头所指)。基于迭代的高分辨率Radon变换(迭代20次)和基于低频约束的高分辨率Radon变换具有相近的多次波压制结果，但后者计算速度要快很多，对于该CDP道集，前者用时23.68 s，而后者仅用时3.49 s。在处理大批量实际数据时，基于低频约束的高分辨率Radon变换算法具有明显的优势。