基于迭代阈值收缩的高分辨率Radon变换方法效果对比

2021-04-17马继涛廖震齐娇迟麟

物探与化探 2021年2期

马继涛，廖震，齐娇，迟麟

(中国石油大学(北京) 地球物理学院，北京 102249)

0 引言

Radon 变换是数学领域的一种变换方法，在地震数据处理的速度分析、波场分离和多次波压制等方面均有应用。但受地震数据采集空间有限及方法自身原因限制，其常规方法的变换域存在平滑效应，且分辨率较低，严重影响方法的效果，如何提高变换域分辨率一直是地球物理学家们研究的重点[1]。

Radon变换由Claerbout[2]为首的斯坦福大学地球物理专家于20世纪80年代将其引入地球物理领域。Thorson和Claerbout[3]使用双曲Radon变换进行速度分析和反演，该方法是在时间—空间域进行求解，需要求取大型矩阵的逆，计算量大；Hampson[4]对抛物线Radon变换的反变换进行了重新定义，并应用最小二乘方法求得地震数据的抛物线Radon变换，该方法被首次应用到多次波压制中，但由于地震数据采集的孔径效应[5]，变换域的分辨率较低。Beylkin[6]对离散Radon变换进行了系统阐述，提出在频域的最小二乘(LS)计算方法，避免了大型矩阵的求逆过程。

为提高Radon变换的分辨率，Scales等[7]提出用迭代重加权最小二乘算法[8](IRLS)来解决Thorson和Claerbout构建的稀疏逆问题RT模型，Sacchi和Ulrych[9]使用概率密度函数，基于贝叶斯原理提出了一种在频域的迭代性稀疏Radon算法FSRT[10](frequency sparse radon transform)，该方法通过增强模型的平滑度来避免放大观测中产生的随机误差，具有较强的抗噪能力，被广泛的用于地震信号的处理；但该方法在处理小时窗中时差十分接近的同相轴，和具有显著AVO效应的数据时会生成假象；Abbad等[11-12]改进了Sacchi和Ulrych提出的FSRT方法，通过引入新的变量消除了变换算子对频率的依赖性，加快了迭代速度，提高了计算效率；时间域稀疏Radon变换TSRT(time sparse radon transform)是一种比FSRT分辨率更高的变换方法，此法沿着曲率和截距时间方向上施加稀疏约束[13]，能获得更具稀疏性的Radon模型；混合时频域稀疏Radon变换TSRTMD[14](time sparse radon transform mixed domain)在频率和时间混合域中进行稀疏Radon变换，一般的TSRTMD中稀疏Radon变换通过对最小二乘算法得到的时间域Radon模型进行收敛运算，正向和反向Radon变换均在频域中实现，Lu[15]在此法基础上引入迭代对计算进行加速，通过对最小二乘算法得到的Radon模型进行迭代收敛处理，进一步使Radon模型更加接近原始值，以达到提高Radon变换分辨率的效果。经模拟数据和实际数据验证，表明此方法能够在很大程度上增强Radon变换的分辨率。该类迭代收缩算法涉及到矩阵的向量乘法运算，能够在每次迭代后进行收缩步长处理。

迭代收敛系列算法通过在时频域中进行稀疏Radon变换[16]来分离噪声信号。概括地说，此类迭代方法中，每次迭代过程中，先对地震信号的Radon域模型乘以变换算子，然后用频域的信号与其结果作差，之后乘以变换算子的伴随矩阵，最后对所得结果的标量进行收缩处理步骤，实现迭代过程中的Radon模型向原始值靠近，收敛Radon域能量团，以达到提高Radon变换分辨率的效果。迭代算法区别于一般的时域和频域Radon变换，由于其是在频域进行正向和反向Radon变换，时域上进行相关的迭代收缩步骤，所以它具有时域RT和频域RT两者的优点，比单一域的RT也更具稀疏性[17]，而且不同于单纯的LS算法处理，迭代算法在时域上对LS得到的Radon模型进行迭代，在每次迭代后能得到更为收敛的Radon模型，这对提高Radon域的数据精度方面非常有效。

迭代收敛算法系列中有各种迭代收缩方法，它们之间各有差异。例如迭代收缩阈值算法[18](iterative shrinkage thresholding,IST)，IST基础上的快速迭代收缩阈值算法(fast iterative shrinkage thresholding,FIST)，稀疏加速时不变迭代收缩阈值算法(sparse radon transform iterative shrinkage,SRTIS)等。对IST、FIST、SRTIS这3种迭代方法的系统阐述，了解该3种方法的收缩作用原理以及实现方式；利用模拟数据验证3种算法，并对3种算法的收敛速度等进行对比分析；利用实际地震资料验证3种迭代收缩算法的收敛效率和处理效果。

1 迭代收缩算法原理

设m为Radon域模型，d为时空域地震数据，x为偏移距，q为描述抛物线曲线弯曲程度的曲线参数，t、τ为时间，对地震道集d(x,t)，可将其视为由m沿特定曲率的抛物线叠加求和[19]而得，该变换的离散形式如下：

(1)

式(1)可写成算子形式D=LM，利用最小二乘方法求解M，可得：

M=(LLH+μI)-1LD，

(2)

式中，M为Radon域模型，L为变换算子，μ为阻尼因子，D为时空域数据。

由于地震数据采集空间有限，求和过程只能在有限的时空域进行，最小二乘法求解得到的Radon域数据存在剪刀状发散的假象，一次波和多次波无法聚焦，使得能量团之间存在相互重叠，从而会导致在Radon域进行切除时对有效信号造成损伤，导致地震数据失真。因此提高Radon域数据的分辨率一直是该方法研究的热点。常规提高分辨率的算法是在频率域进行的，本文在此讨论对比几种基于迭代收缩阈值的时间域提高分辨率的算法。

1.1 迭代收缩阈值算法原理

常规的迭代收缩阈值算法(IST)原理如下：

mk+1=

Sλ{mk+αF-1[(LTL)-1LT(F[d]-LF[mk])]}(3)

(4)

其中，mk+1是第k+1次迭代后的Radon模型，mk是第k次迭代后的Radon模型，α是步长，Sλ为阈值算子，F和F-1为傅立叶正、反变换算子，λ为阈值运算所对应的参数。|m|为Radon模型的绝对值。

从式(3)、(4)可以看出，与最小二乘算法相比(式2)，IST在时间域多了一个迭代的步骤，且每一步迭代仅包含1次简单的阈值收缩操作，因此能够大大提高计算效率。在每次阈值处理过程中，增强能量团中心，减弱边缘的发散能量，可以使得模型本身能量团更为集中，由此消除平滑效应带来的分辨率降低的影响。

1.2 快速迭代收缩阈值算法原理

迭代收缩阈值算法尽管可以提高分辨率，但其收敛速度较为缓慢。因此，在常规迭代收缩阈值算法的基础上，Beck、Teboulle等[20]提出了一种改进的快速迭代收缩阈值算法FIST。此改进方法大大加快了迭代的收敛速度，其公式与IST相比增加了一迭代计算流程，如式(5)：

mk+1=

(6)

式中，初始值m0=0，t0=1，并且tk由下式计算得到：

(7)

从以上迭代步骤可以看出，FIST在原有算法的基础上增加了一重迭代计算，增加的算式旨在对Radon域数据模型迭代收敛前先进行残差校正，因此能使相同迭代次数下FIST加快减小与原始数据模型的差距，增大收敛速度。

1.3 稀疏加速时不变迭代收缩阈值算法

FIST虽然在收敛速度上提升了不少，但由于此方法只添加了一重迭代步骤，并没有从本质上改变收敛效果。基于TSRTMD，Lu[21]给出了一种改进的SRTIS算法，此法具备更高的迭代收缩效率，不仅能进一步加快收敛速度，而且也能极大地改善收敛效果，使得变换域的分辨率更高。

迭代收缩算法通常应用全局阈值，在Radon变换中，Radon域模型通常沿着时间维度衰减，因此为遵循模型衰减，将模型收缩用二维方式来展示：

mk=Tα{mk-1+2tF-1[(LTL)-1LT(F[d]-LF[mk-1])]}

(8)

其中，t是步长，k为当前迭代次数，Tα为一收缩算子，计算公式如下：

(10)

(11)

当均值滤波器大小与模型m相同时，二维收缩将达到最好，阈值将降低到全局阈值。

SRTIS是一个混合的时频域Radon变换，该法与其他迭代算法类似，均在频域中完成对地震数据的正向和反向Radon变换，在时域中对LS求得的Radon模型进行迭代收缩步骤，其不同点在于式(9)的收缩算法中，SRTIS先对初始Radon模型m0的绝对值作二维均值滤波处理，再将得到的结果归一化后作为收缩算子Tα的参数值，这种处理方式不同于上述的IST和FIST中预先设定的固定阈值参数，它利用每次迭代后Radon模型的更新来进行动态阈值，能够保证每次迭代后的收缩处理更加贴合原始数据，因此其收敛效率和计算效果大大优于IST。

2 方法应用效果

2.1 模拟数据

本文利用一个简单的模拟地震数据对比3种方法的处理效果，该数据有80道，道间距为20 m，每道采样点为200个，采样率为4 ms。图1a给出了时空域的模拟数据，图中可以看出，数据中5条曲线振幅较强，振幅值为1，两条曲线振幅较弱，振幅值为0.1。图1b为对应的最小二乘Radon变换方法的变换结果。

图1 时空域模拟数据(a)及对应的最小二乘Radon变换结果(b)Fig.1 The simulated data in the time domain(a) and the corresponding result of the LS Radon transform(b)

2.1.1 IST方法处理结果

对最小二乘Radon变换的结果，基于IST方法，利用式(3)、式(4)对其进行迭代收缩处理，迭代次数为40、80、160次时的结果如图2所示。

从图2中可以看出，IST对Radon域的能量团起到了一定的收敛作用，但收敛速度缓慢。

之后对IST方法迭代得到的结果进行多次波压制，即将动校时差小于0.04 s的Radon域数据(一次波)切除掉，得到多次波压制结果，如图3所示。

比较图3中的多次波压制结果，IST反变换得到的Radon域模型仍与原数据存在较大差异，导致多次波压制结果存在残余。

从式(4)可以看出，λ对mk的输出结果有所影响，λ数值应比mk小两个数量级，且一定范围内，λ越大，Radon域收敛性越强，具有更高的分辨率，但会导致时间域多次波压制后的残余更严重。因此选取λ为0.001、0.001 1、0.001 3、0.001 5、0.001 7、0.001 9，迭代次数为160次时进行对比，且此数据中，mk的数量级为10-2。结果显示，将λ值取到0.001 9时能达到较好效果。

a—迭代次数为40次;b—迭代次数为80次;c—迭代次数为160次a—the 40 iterations;b—the 80 iterations;c—the 160 iterations图2 IST不同迭代次数结果对比Fig.2 The iterative results with different iterations of the IST transform

a—最小二乘方法;b—迭代次数为40次;c—迭代次数为160次a—the LS transform;b—the 40 iterations;c—the 160 iterations图3 IST方法多次波压制结果对比Fig.3 The comparison of multiple suppression results using IST method

2.1.2 FIST方法处理结果

FIST方法能降低估计Radon域模型的计算成本，同时加快迭代收敛的速度。基于FIST方法，利用式(5)～(7)对最小二乘方法的Radon变换结果进行迭代收缩处理，迭代次数为10、20、40次时的结果如图4所示。从图中可见，与IST 方法相比，FIST收敛速度确实加快了。

FIST方法还在一定程度上消除了IST方法中多次波压制残余的现象，在图5中可以看出，随着迭代次数的增加，压制结果中的多次波残余逐渐减小，但仍有部分残余存在。

2.1.3 SRTIS方法处理结果

基于SRTIS方法，利用式(8)～(11)对最小二乘的Radon变换结果进行迭代收缩处理，迭代次数为5、15、30次时的结果如图6所示，对应的多次波压制结果图7所示。

a—迭代次数为10次;b—迭代次数为20次;c—迭代次数为40次a—the 10 iterations;b—the 20 iterations;c—the 40 iterations图4 FIST迭代结果Fig.4 The iterative result of the FIST transform

a—迭代次数为10次;b—迭代次数为20次;c—迭代次数为40次a—the 10 iterations;b—the 20 iterations;c—the 40 iterations图5 FIST多次波压制结果对比Fig.5 The comparison of multiple suppression results using FIST method

a—迭代次数为5次;b—迭代次数为15次;c—迭代次数为30次a—the 5 iterations;b—the 15 iterations;c—the 30 iterations图6 SRTIS迭代结果Fig.6 The iterative result of the SRTIS transform

a—迭代次数为5次;b—迭代次数为15次;c—迭代次数为30次a—the 5 iterations;b—the 15 iterations;c—the 30 iterations图7 SRTIS多次波压制结果对比Fig.7 The comparison of multiple suppression results using SRTIS method

图4和图6对比可见，SRTIS在收敛速度上略快于FIST，在迭代次数15次时，SRTIS已经能看到清晰的数据点，能量充分聚集，平滑效应被消除，但两个弱振幅点并不能显示，在迭代次数达到30次时能看到7个数据点。从图7可见在迭代次数为30次时多次波残余基本消失，说明SRTIS的压制效果最佳。在上述SRTIS结果中，α和t均选择同一数值，α=0.04，t=0.8，此参数设置为经过验证后的最佳试验数值。

2.1.4 3种方法的结果对比

为对比3种方法的数据处理能力，本文对3种迭代算法的收敛效率进行了比较，3种算法1次到20次迭代后的模型mk与原始模型m的误差能量结果如图9所示。图中可见SRTIS下降速率最快，FIST次之，IST最慢。说明SRTIS方法的收敛性优于其他两种，改进后的FIST优于IST。

图9 3种迭代算法的收敛效率对比Fig.9 The compasion of three iterative algorithms' contraction efficiency

为比较三者之间的收敛速度，用模拟数据以及一总长度为7 s、采样间隔为0.004 s的CDP道集计算达到10-1级别的精度标准时所耗费的时间，结果见表1所示。

从表1中可见IST方法的收敛速度非常缓慢，而FIST与SRTIS收敛速度没有很大差别，但是在保幅性能上，FIST并不能很好的达到要求，在实际资料的处理中，较易破坏同相轴，因此综合考虑，SRTIS的收敛效率较高。

此外，用一相近抛物线试验这些方法的压制干扰能力，其动较时差为0.028 s，压制动校时差大于0.02 s的抛物线，两者相差仅为0.008 s。4种方法中，最小二乘算法残余严重；IST能在一定程度上切除掉干扰抛物线，但是在尾端仍有明显的残余；FIST尾端的残余现象有得到改善；SRTIS压制效果十分明显，压制后已经很难观测到干扰的抛物线痕迹。

表1 3种迭代算法的收敛结果对比

以上结果均说明SRTIS优于其他两种迭代收缩算法，能较大程度上提高Radon变换的分辨率，且能较为彻底地对数据中剩余时差较小的多次波进行压制，因此方法对压制剩余时差较小的层间多次波较为有利。

2.2 实际地震资料处理结果

为进一步验证迭代收缩算法的实用性，本文引入实际地震资料，基于3种迭代收缩方法与传统最小二乘算法进行对比。实际资料为墨西哥湾的某多次波干扰严重的CDP道集。

利用最小二乘方法进行了Radon变换和多次波的压制，并基于该变换结果和上文所述的各种迭代收缩方法进行提高分辨率运算，迭代次数均为40次，得到高分辨率的Radon域结果；对所得到的Radon域数据进行一次波和多次波的分离，并将分理出的多次波变换回时空域，将其由原数据中减去，得到多次波压制后的道集数据。原始数据、LS、IST、FIST和SRTIS各方法估计的多次波模型如图10所示，多次波压制后的结果如图11所示。

图10 各方法所估计的多次波模型Fig.10 The estimated multiple model by various methods

图11 各方法多次波压制后的结果Fig.11 The multiple suppression result by various methods

上述方法的处理效果中，LS方法利用最小二乘原理取得Radon域的变换结果，但最小二乘算法仅为局部收敛，因此估计的多次波与原始数据仍有较大差异，压制的时候会影响到一次波的反变换结果；IST虽然提高了分辨率，但其公式中并未将每次计算迭代后的残差量带回，因此导致所估计多次波与原数据有一定差异，且IST方法远偏移距箭头所指一次波也受到影响；FIST方法虽能进一步提高分辨率，也考虑了残差量，但因其迭代程度相对过深，导致Radon域数据失真；SRTIS因为将每次迭代后的数据模型如式(10)归一化后代入到收缩阈值步骤中，在一定程度上能够约束模型的收敛方向，使得其更靠近真解，因此估计的多次波较为真实，压制效果也与实际情况较为吻合，既提高了Radon域分辨率，也得到了相对较好的压制结果，近偏移距和远偏移距效果均不错。且从图11压制效果中可见，图11e的结果最好的保持了原数据的振幅特性，而图11c、d中一次波的横向振幅特性遭到破坏，导致一次波水平同相轴较多，不符实际。

为突出振幅能量，本文将处理后的CDP道集叠加用于对比，叠加道集效果如图12所示。

对比图11的单道压制效果，图12叠加道集在4.3 s附近的压制效果较为明显。图12中，LS方法的叠加大致能够将3.5～3.7 s出现的干扰同相轴去除，但仍残留着多次波的干扰能量，这是该方法的理论公式所带来的误差影响，是无法消除的；IST的叠加虽然提高了分辨率，但由于在压制过程中受到多次波残余误差的影响，并不能很好地辨识出一次波的同相轴，信噪比较低，仅能在4～4.3 s模糊地看到3条同相轴，且4.7 s和5.1 s附近的同相轴聚焦并不是很好；FIST方法的叠加进一步提高了分辨率，对多次波的压制程度也更大，但因其迭代程度相对过深，造成的信噪比损失也更为严重，因此4～4.3 s处的同相轴周围干扰影响更为严重；SRTIS的叠加因为其方法在每次迭代过程中都能与原始数据模型保持一致性，所以既提高了分辨率，也并没有对信噪比造成较大破坏，平衡了分辨率和信噪比，在图中能够清楚地看到4～4.3 s的3条同相轴，且4.5 s和5.1 s附近的同相轴相对集中，于多次波压制为较好的结果。综上所述，SRTIS在实际地震资料的多次波压制应用中依然强于其他两种迭代方法。