基于迭代阈值收缩的高分辨率Radon变换方法效果对比
中国石油大学(北京) 地球物理学院,北京 102249
The comparison of effects of high-resolution Radon transform based on iterative shrinkage thresholding
College of Geophysics,China University of Petroleum(Beijing),Beijing 102249,China
责任编辑: 叶佩
收稿日期: 2019-12-17 修回日期: 2020-11-22 网络出版日期: 2021-04-20
基金资助: |
|
Received: 2019-12-17 Revised: 2020-11-22 Online: 2021-04-20
作者简介 About authors
马继涛(1983-),男,山东临清,博士,副教授,硕士生导师,2009年获中国石油大学(北京)地质资源与地质工程专业博士学位,现主要从事地震数据处理相关的教学和科研工作。
提高Radon变换分辨率的研究一直是地震数据处理方法研究的热点之一。常用的提高分辨率的算法是在频率域进行的,然而此方法在频域中计算的模型权重是偶然耦合的,会对所有同相轴施加相同的权重,导致高能量同相轴地震道生成假像。此文给出了3种时间域提高Radon变换分辨率的方法,IST(iterative shrinkage thresholding)、FIST(fast iterative shrinkage thresholding)和SRTIS(sparse radon transform iterative shrinkage),并对3种算法的收敛效率和计算效果进行了对比分析。理论和实际数据的测试表明,SRTIS方法在效率和效果的对比上均优于其他两种,并且具有较好的去除多次波能力。
关键词:
The study of the resolution improvement of Radon transform is one of the research hotspots in seismic data processing area.The commonly-used resolution improvement methods are carried out in the frequency domain.However,the weighting of the Radon model in the frequency domain is coupled,and it will impose the same weight on all the events,leading to the artifacts generated by high energy seismic events.This paper presents three resolution improvement methods of Radon transform in the time domain:Iterative Shrinkage Thresholding (IST),Fast Iterative Shrinkage Thresholding (FIST) and Sparse Radon Transform Iterative Shrinkage (SRTIS),with a comparison of their compute efficiencies and results.Synthetic data and real data test results show that SRTIS is superior to the other two methods in computation effect and efficiency,and it has a better multiple attenuation capability.
Keywords:
本文引用格式
马继涛, 廖震, 齐娇, 迟麟.
MA Ji-Tao, LIAO Zhen, QI Jiao, CHI Lin.
0 引言
Radon 变换是数学领域的一种变换方法,在地震数据处理的速度分析、波场分离和多次波压制等方面均有应用。但受地震数据采集空间有限及方法自身原因限制,其常规方法的变换域存在平滑效应,且分辨率较低,严重影响方法的效果,如何提高变换域分辨率一直是地球物理学家们研究的重点[1]。
为提高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)可写成算子形式D=LM,利用最小二乘方法求解M,可得:
式中,M为Radon域模型,L为变换算子,μ为阻尼因子,D为时空域数据。
由于地震数据采集空间有限,求和过程只能在有限的时空域进行,最小二乘法求解得到的Radon域数据存在剪刀状发散的假象,一次波和多次波无法聚焦,使得能量团之间存在相互重叠,从而会导致在Radon域进行切除时对有效信号造成损伤,导致地震数据失真。因此提高Radon域数据的分辨率一直是该方法研究的热点。常规提高分辨率的算法是在频率域进行的,本文在此讨论对比几种基于迭代收缩阈值的时间域提高分辨率的算法。
1.1 迭代收缩阈值算法原理
常规的迭代收缩阈值算法(IST)原理如下:
其中,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):
式中,Sλ为阈值算子,其计算与IST算法相同。
式中,初始值m0=0,t0=1,并且tk由下式计算得到:
从以上迭代步骤可以看出,FIST在原有算法的基础上增加了一重迭代计算,增加的算式旨在对Radon域数据模型迭代收敛前先进行残差校正,因此能使相同迭代次数下FIST加快减小与原始数据模型的差距,增大收敛速度。
1.3 稀疏加速时不变迭代收缩阈值算法
FIST虽然在收敛速度上提升了不少,但由于此方法只添加了一重迭代步骤,并没有从本质上改变收敛效果。基于TSRTMD,Lu[21]给出了一种改进的SRTIS算法,此法具备更高的迭代收缩效率,不仅能进一步加快收敛速度,而且也能极大地改善收敛效果,使得变换域的分辨率更高。
迭代收缩算法通常应用全局阈值,在Radon变换中,Radon域模型通常沿着时间维度衰减,因此为遵循模型衰减,将模型收缩用二维方式来展示:
其中,t是步长,k为当前迭代次数,Tα为一收缩算子,计算公式如下:
其中,α为阈值参数,通常取小于1的正值。K为设定的最大迭代次数,m是Radon域模型,
其中,
当均值滤波器大小与模型m相同时,二维收缩将达到最好,阈值将降低到全局阈值。
SRTIS是一个混合的时频域Radon变换,该法与其他迭代算法类似,均在频域中完成对地震数据的正向和反向Radon变换,在时域中对LS求得的Radon模型进行迭代收缩步骤,其不同点在于式(9)的收缩算法中,SRTIS先对初始Radon模型m0的绝对值作二维均值滤波处理,再将得到的结果归一化后作为收缩算子Tα的参数值,这种处理方式不同于上述的IST和FIST中预先设定的固定阈值参数,它利用每次迭代后Radon模型的更新来进行动态阈值,能够保证每次迭代后的收缩处理更加贴合原始数据,因此其收敛效率和计算效果大大优于IST。
2 方法应用效果
2.1 模拟数据
图1
图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
图2
IST不同迭代次数结果对比
a—迭代次数为40次;b—迭代次数为80次;c—迭代次数为160次
Fig.2
The iterative results with different iterations of the IST transform
a—the 40 iterations;b—the 80 iterations;c—the 160 iterations
从图2中可以看出,IST对Radon域的能量团起到了一定的收敛作用,但收敛速度缓慢。
之后对IST方法迭代得到的结果进行多次波压制,即将动校时差小于0.04 s的Radon域数据(一次波)切除掉,得到多次波压制结果,如图3所示。
图3
图3
IST方法多次波压制结果对比
a—最小二乘方法;b—迭代次数为40次;c—迭代次数为160次
Fig.3
The comparison of multiple suppression results using IST method
a—the LS transform;b—the 40 iterations;c—the 160 iterations
比较图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时能达到较好效果。
2.1.2 FIST方法处理结果
FIST方法能降低估计Radon域模型的计算成本,同时加快迭代收敛的速度。基于FIST方法,利用式(5)~(7)对最小二乘方法的Radon变换结果进行迭代收缩处理,迭代次数为10、20、40次时的结果如图4所示。从图中可见,与IST 方法相比,FIST收敛速度确实加快了。
图4
图4
FIST迭代结果
a—迭代次数为10次;b—迭代次数为20次;c—迭代次数为40次
Fig.4
The iterative result of the FIST transform
a—the 10 iterations;b—the 20 iterations;c—the 40 iterations
FIST方法还在一定程度上消除了IST方法中多次波压制残余的现象,在图5中可以看出,随着迭代次数的增加,压制结果中的多次波残余逐渐减小,但仍有部分残余存在。
图5
图5
FIST多次波压制结果对比
a—迭代次数为10次;b—迭代次数为20次;c—迭代次数为40次
Fig.5
The comparison of multiple suppression results using FIST method
a—the 10 iterations;b—the 20 iterations;c—the 40 iterations
2.1.3 SRTIS方法处理结果
图6
图6
SRTIS迭代结果
a—迭代次数为5次;b—迭代次数为15次;c—迭代次数为30次
Fig.6
The iterative result of the SRTIS transform
a—the 5 iterations;b—the 15 iterations;c—the 30 iterations
图7
图7
SRTIS多次波压制结果对比
a—迭代次数为5次;b—迭代次数为15次;c—迭代次数为30次
Fig.7
The comparison of multiple suppression results using SRTIS method
a—the 5 iterations;b—the 15 iterations;c—the 30 iterations
通过误差能量计算式R(r)=‖r-
图8
图8
SRTIS迭代结果与原始数据模型的误差能量
Fig.8
The error energy between the iterative result of SRTIS and the original data
2.1.4 3种方法的结果对比
为对比3种方法的数据处理能力,本文对3种迭代算法的收敛效率进行了比较,3种算法1次到20次迭代后的模型mk与原始模型m的误差能量结果如图9所示。图中可见SRTIS下降速率最快,FIST次之,IST最慢。说明SRTIS方法的收敛性优于其他两种,改进后的FIST优于IST。
图9
图9
3种迭代算法的收敛效率对比
Fig.9
The compasion of three iterative algorithms' contraction efficiency
为比较三者之间的收敛速度,用模拟数据以及一总长度为7 s、采样间隔为0.004 s的CDP道集计算达到10-1级别的精度标准时所耗费的时间,结果见表1所示。
表1 3种迭代算法的收敛结果对比
Table 1
算法 | 模拟数据 | 实际数据 | ||
---|---|---|---|---|
达到同精度时间/s | 所需迭代次数/次 | 达到同精度时间/s | 所需迭代次数/次 | |
IST | 1.78 | 81 | 55.77 | 87 |
FIST | 1.32 | 35 | 22.52 | 25 |
SRTIS | 1.34 | 37 | 28.26 | 31 |
从表1中可见IST方法的收敛速度非常缓慢,而FIST与SRTIS收敛速度没有很大差别,但是在保幅性能上,FIST并不能很好的达到要求,在实际资料的处理中,较易破坏同相轴,因此综合考虑,SRTIS的收敛效率较高。
此外,用一相近抛物线试验这些方法的压制干扰能力,其动较时差为0.028 s,压制动校时差大于0.02 s的抛物线,两者相差仅为0.008 s。4种方法中,最小二乘算法残余严重;IST能在一定程度上切除掉干扰抛物线,但是在尾端仍有明显的残余;FIST尾端的残余现象有得到改善;SRTIS压制效果十分明显,压制后已经很难观测到干扰的抛物线痕迹。
以上结果均说明SRTIS优于其他两种迭代收缩算法,能较大程度上提高Radon变换的分辨率,且能较为彻底地对数据中剩余时差较小的多次波进行压制,因此方法对压制剩余时差较小的层间多次波较为有利。
2.2 实际地震资料处理结果
为进一步验证迭代收缩算法的实用性,本文引入实际地震资料,基于3种迭代收缩方法与传统最小二乘算法进行对比。实际资料为墨西哥湾的某多次波干扰严重的CDP道集。
图10
图11
上述方法的处理效果中,LS方法利用最小二乘原理取得Radon域的变换结果,但最小二乘算法仅为局部收敛,因此估计的多次波与原始数据仍有较大差异,压制的时候会影响到一次波的反变换结果;IST虽然提高了分辨率,但其公式中并未将每次计算迭代后的残差量带回,因此导致所估计多次波与原数据有一定差异,且IST方法远偏移距箭头所指一次波也受到影响;FIST方法虽能进一步提高分辨率,也考虑了残差量,但因其迭代程度相对过深,导致Radon域数据失真;SRTIS因为将每次迭代后的数据模型如式(10)归一化后代入到收缩阈值步骤中,在一定程度上能够约束模型的收敛方向,使得其更靠近真解,因此估计的多次波较为真实,压制效果也与实际情况较为吻合,既提高了Radon域分辨率,也得到了相对较好的压制结果,近偏移距和远偏移距效果均不错。且从图11压制效果中可见,图11e的结果最好的保持了原数据的振幅特性,而图11c、d中一次波的横向振幅特性遭到破坏,导致一次波水平同相轴较多,不符实际。
为突出振幅能量,本文将处理后的CDP道集叠加用于对比,叠加道集效果如图12所示。
图12
图12
实际数据CDP道集叠加对比
a—原始叠加道集;b—LS Radon结果叠加道集;c—IST Radon 结果叠加道集;d—FIST Radon结果叠加道集;e—SRTIS Radon结果叠加道集
Fig.12
The comparison of real data CDP gather stack results
a—the stack of the original data;b—the stack of LS Radon result;c—the stack of IST Radon result;d—the stack of FIST Radon result;e—the stack of SRTIS Radon result
对比图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在实际地震资料的多次波压制应用中依然强于其他两种迭代方法。
3 结论
本文基于IST、FIST和SRTIS这3种迭代收缩Radon变换算法压制多次波,并给出了算法实现的计算步骤,模型数据和实际资料的处理验证了这3种方法的正确性和可行性,得出如下结论:
1)传统LS方法在用于压制多次波方面上存在局限性,而迭代收缩算法通过迭代中的阈值收缩操作来逼近原始值,在每次阈值处理过程中,增强能量团中心,削弱边缘的发散能量,使模型能量团集中,由此消除平滑效应带来的降低分辨率的影响。
2)SRTIS的压制效果均优于IST和FIST,其先对初始Radon模型m0的绝对值作二维均值滤波处理,再将得到的结果归一化后作为收缩算子Tα的参数值。不同于IST和FIST的固定阈值参数,SRTIS利用了每次迭代后Radon模型的更新来进行动态阈值,能够保证每次迭代后的收缩处理更加贴合原始数据。
3)参数的选择对于算法的实现效果有很大影响,因此在使用算法处理实际数据时需先对测试数据大小进行评估,才能达到较好的收缩阈值效益。
参考文献
Radon变换压制层间多次波技术在高石梯—磨溪地区的应用
[J].
The application of Radon transform to suppress interbed multiples in Gaoshiti—Moxi region
[J].
高分辨率Radon变换及其在地震资料处理中的应用
[D].
High resolution Radon transform and its application in seismic data processing
[D].
Velocity stack and slant stochastic inversion
[J].
Inverse velocity stacking for multiple elimination
[J].
高分辨率抛物线拉冬变换多次波压制技术
[J].
High resolution parabolic radon transform multiple wave suppression technique
[J].
Discrete Radon transform
[J].
Fast lp solution of large,sparse,linear systems:Application to seismic travel time tomography
[J].
The simplest discrete Radon transform
[C]
High-resolution velocity gathers and offset space reconstruction
[J].
Fast wavelet transforms and numerical algorithms
[J].
A fast,modified parabolic Radon transform
[J].
A wide-angle view at iterated shrinkage algorithms
[C]
Applications of time-domain high-resolution Radon demultiple
[C]
Latest views of the sparse Radon transform
[J].DOI:10.1190/1.1543224 URL [本文引用: 1]
A time-domain high-resolution Radon transform based on iterative model shrinkage
[C]
An iterative thresholding algorithm for linear inverse problems with a sparsity constraint
[J].
L1-L2 optimization in signal and image processing
[J].
An EM algorithm for wavelet based image restoration
[J].
De-multiple via a fast least squares hyperbolic Radon transform
[C]
A fast iterative shrinkage-thresholding algorithm for linear inverse problems
[J].
An accelerated sparse time-invariant Radon transform in the mixed frequency-time domain based on iterative 2D model shrinkage
[J].
/
〈 |
|
〉 |
