引用本文
贾春梅, 姜国庆, 刘志成, 许璐. 频域稀疏双曲Radon变换去噪方法[J]. 物探与化探, 2016,40(3): 527-533.
JIA Chun-Mei, JIANG Guo-Qing, LIU Zhi-Cheng, XU Lu. Denoising method based on sparse hyperbolic Radon transform in the frequency domain[J]. Geophysical and Geochemical Exploration, 2016,40(3): 527-533.
Doi:10.11720/wtyht.2016.3.14  
Permissions
Copyright©2016, 《物探与化探》编辑部
《物探与化探》编辑部 所有
频域稀疏双曲Radon变换去噪方法
贾春梅1, 姜国庆2, 刘志成1, 许璐1
1. 中国石油化工股份有限公司 石油物探技术研究院,江苏 南京 211103
2. 江苏省地质调查研究院 基础地质研究所,江苏 南京 210018

作者简介: 贾春梅(1986-),女,物探工程师,硕士研究生,2011年毕业于中国石油大学(北京),主要从事地震资料去噪方法研究工作。E-mail:jiacm86@163.com

收稿日期: 2015-05-13
基金: 中石化科技部项目(P14127); 徐州城市地质调查、徐州地区岩溶塌陷调查项目(12120114022001)
摘要

Radon变换是地震数据处理中广泛应用的一种技术方法,但实际应用中存在分辨率及运算效率低的问题。笔者在相移双曲线的基础上,推导了时不变双曲积分路径,提出了基于贝叶斯原理的稀疏双曲Radon变换,给出离散变换的参数选取准则,并且为保证不同反射层的反射波具有最佳收敛,在变换中引入了动态时窗。频域稀疏双曲Radon变换不仅具有较高的分辨率,而且可以在频域实现,计算效率大幅提高。模型数据测试和实际资料应用结果表明,频域稀疏双曲Radon变换可以有效去除规则干扰波、压制随机噪声,提高数据信噪比,且具有较好的保真度和保幅性。

关键词: Radon变换; 稀疏Radon变换; 随机噪声压制; 时不变双曲线
中图分类号:P631.4 文献标志码:A 文章编号:1000-8918(2016)03-0527-07 doi: 10.11720/wtyht.2016.3.14
Denoising method based on sparse hyperbolic Radon transform in the frequency domain
JIA Chun-Mei1, JIANG Guo-Qing2, LIU Zhi-Cheng1, XU Lu1
1.Research Institute of Geophysical Prospecting,SINOPEC,Nanjing 211103,China
2.Basic Geological Research Institute,Geological Survey of Jiangsu Province,Nanjing 210018,China
Abstract

Radon transform is a kind of technical method widely used in seismic data processing,but there are problems of low resolution and low computational efficiency in practical application.On the basis of phase shift hyperbola,the authors deduced time-invariant hyperbolic integral path,proposed sparse hyperbolic Radon transform based on Bayes principle,provided parameter selection principle of discrete transform,and introduced dynamic time windows to ensure the condition that reflection waves of different layers have the best convergence.Sparse hyperbolic Radon transform in the frequency domain has higher resolution and it can be realized in the frequency domain,so the computational efficiency can be drastically improved.Model data test and actual application result show that sparse hyperbolic Radon transform in the frequency domain can effectively remove regular interference,suppress random noise and improve the signal-to-noise ratio;in addition,it also has better fidelity and amplitude preservation.

Keyword: Radon transform; sparse Radon transform; random noise suppressing; time invariant hyperbolic

拉东变换由奥地利数学家Radon于1917年提出, 该变换在医学、天文学及物理学中得到了特别广泛的应用。20世纪70年代, 美国斯坦福大学Claerbout、Chapman等率先将Radon变换引入到地震勘探中[1]; Sacchi和Ulrych通过研究提出了利用Radon变换的稀疏解来提高Radon变换域分辨率的思想[2]; Cary总结了离散Radon变换所存在的截断效应和空间假频问题[3]; Sacchi根据稀疏解的思想给出了利用稀疏解来提高Radon变换域分辨率的算法[4]; Nowak提出了加权的Radon变换压制多次波方法, 可以很好地保持变换后的振幅, 便于AVO分析[5]; Yu在Wavelet-Radon域进行了地震数据重建和去除空间假频的研究[6]。Juefu Wang提出应用贪婪算法代替共轭梯度法来提高Radon变换的分辨率和收敛速率, 取得了较好的效果[7]。国内学者对Radon变换也进行了许多研究工作, 王维红基于部分动校正(NMO) 后反射同相轴在CMP 道集上的抛物线走时近似, 给出了加权抛物Radon 变换叠前地震数据重建方法(WPRT)[8]; 熊登发展了混合域高分辨率抛物Radon变换[9]; 石颖结合基于波动方程预测和双曲Radon变换联合压制表面多次波[10]; 鲁娥采用混合Radon变换分离面波等线性噪声, 采用自适应滤波技术在Radon域识别并剔除多次波[11]; 张振波采用高分辨率抛物线Radon变换实现对地震数据中多次波及空间假频的快速、有效压制[12]。

Radon变换是沿着特定路径(如直线、抛物线、双曲线)对介质的某个特性线进行积分, 根据积分路径的不同, Radon变换可分为线性Radon变换、抛物Radon变换和双曲Radon变换, 其中双曲Radon变换的积分路径最接近于有效波反射同相轴, 具有最高的精度。传统双曲Radon变换具有较高的分辨率, 但是其积分路径不具有时不变性, 无法转化为频率— 空间域来求解, 因此, 直接进行双曲Radon变换成本非常高, 以致难以实现。为了解决这一问题, 引入时不变双曲积分路径进行频域稀疏Radon变换, 该方法不仅克服了传统双曲Radon变换计算效率低下的问题, 而且提高了分辨率, 更加有利于随机噪声的压制。

1 理论实现
1.1 频域稀疏双曲Radon变换

在水平层状介质模型中, 常规的动校正Dix方程[13]

t=t02+x2/vrms2(1)

是在小偏移距时对于水平层状介质模型的近似得出的。在这个方程中, t是波从震源到接收点的传播时间; t0是表面到反射体的垂直反射双程旅行时; x是炮点到检波点的距离; vrms为水平层状介质的均方根速度。

将Dix公式按照高斯超几何级数序列展开简化得到相移双曲公式

t=τs+τ02+x2/v2, (2)

其中:τ 0=t0/S; τ s=τ 0(S-1); v2=S vrms2; S=μ 4/ μ22。相移双曲线就是常规动校正双曲线沿时间轴进行了τ s的时移, 相移Radon变换的积分路径是时不变的, 其求解方程对频率是解耦的, 可提高双曲Radon正变换的计算效率。另外, 对于均匀采样的地震数据, 时不变双曲Radon变换法所形成的频率— 空间域矩阵具有Toeplitz结构, 应用Levinson递推算法可大大提高矩阵求解计算的效率。

图1所示为Dix常规双曲线与相移双曲线分别与反射同相轴进行对比, 由图中可以看出Dix的常规动校正曲线只是在近偏移距比较接近真实曲线, 而相移双曲线比常规动校正双曲线更好地接近于实际反射同相轴。

图1 Dix动校正公式(左)与相移动校正公式(右)对比

为了使时不变双曲线积分路径内的参数物理意义更加清楚, 同时便于赋值和计算, 对公式(2)进行了进一步推导, 得出了新的时不变双曲积分路径:

t=τ+h(x2+zref2-zref), (3)

式中:τ 为截距时间, h为射线参数, zref为参考深度。这就是文中所应用的时不变双曲线的积分路径, 由于t和τ 具有时不变的关系, 因此可以变换到频率域进行计算, 大大减少了工作量。

频域Radon变换用矩阵形式可以写为:

M=LHD, 　D=LM, (4)

式中:M是双曲Radon正变换剖面经Fourier变换后频域值; D是时空域地震数据经Fourier变换后频域值; LH是矩阵L的共轭转置, 矩阵LH和L的元素分别可以表示为:

LHm, n=eiωhmθn, Lm, n=e-iωhnθm。(5)

所谓稀疏Radon变换, 就是构造的Radon变换域为分割的稀疏脉冲, 或称为稀疏解。对于双曲线Radon变换, 笔者采用Cauchy类准则。使目标函数 最小。式中:Mk是Radon变换结果M的元素; μ 和b是相应的分布参数。对目标函数关于M求导, 得:

LHLM-LHD+QM=0。(7)

方程(7)需要迭代求解, 经推导得到解的具体形式为

Mk=LHL+2μb+Mk-12-1LHD。(8)

式中:k为迭代次数。

理论上时间域的一个双曲同相轴经过传统的双曲Radon变换映射到Radon域会聚焦成一个点, 然而在实际情形中并不是按照这种方式。图2a显示了时间域的双曲同相轴通过传统的Radon变换得到的Radon域数据, 图中出现了水平和倾斜能量发散。水平方向的假频是由于近偏移距能量的贡献, 而倾斜方向的假频则是远偏移距数据的截断效应。图2b为时间域双曲同相轴经过频率域稀疏Radon变换计算结果, 通过对比可以看出稀疏Radon变换具有更好的收敛性。

图2 常规Radon变换(左)和稀疏Radon变换(右)结果对比

1.2 参数的选取

对于离散变换都有参数选取的问题, 如果选取不当容易造成假频, 并且得不到稀疏收敛的数据影响处理效果。对于双曲Radon变换, 射线参数h的采样间隔为[10]

Δh≤1/(θrangefmax)。(9)

一般情况下, zref选为最大偏移距xmax, 在这种假设下θ range=0.4xmax, h的采样间隔可以表示为

Δh≤5/(2xmaxfmax)。(10)

如果zref=xmax, 射线参数h的取值范围为

hrange≤2/(Δxfmax)。(11)

1.3 动态时窗的引入

通过上面的推导可以看出, 在频率域稀疏双曲Radon变换中涉及到参考深度, 也就是zref的选取。下面通过一个模型算例对zref的选取原则进行讨论。如图3a为速度模型, 图3b为根据速度模型合成的地震记录, 图3c为zref取900时变换结果, 图3d为zref取3 000时的变换结果。通过对比可以看出, 当参考深度zref取为900时, 经过稀疏双曲Radon正变换后层状模型第一层收敛较好, 二、三层收敛较差; 而当参考深度zref取为3 000时, 变换后层状模型二、三层收敛相对较好, 第一层变换后则明显发散。

图3 参考深度选取模型试算

通过模型试算可知, 在稀疏双曲Radon变换中由浅层到深层采用同一参考深度无法取得较好的收敛效果, 对于不同反射层的反射波, 如果想得到最好的收敛效果必须选取适当的参考深度, 即由浅层到深层的反射必须选取不同的参考深度。因此, 在变换中必须加入动态时窗, 对于同一时窗内的数据采用相同的参考深度, 然后将所有时窗内变换后的数据叠加就可以得到变换域数据。

在本文的变换中, 先对原始数据进行了动校正, 因此采用了矩形时窗进行数据的选取, 引入动态时窗的频率域稀疏双曲Radon变换的方法流程(图4)。

图4 频率域稀疏双曲Radon变换的方法流程

2 模型数据测试
2.1 频域稀疏Radon变换波场分离去除多次波

多次波在地震勘探中普遍存在, 对于海上勘探的影响尤其严重。由于多次波的同相轴为双曲线, 应用文中所研究的频率域稀疏双曲Radon变换能够使其能量在变换域充分收敛, 便于将其切除。图5所示为薄互层模型。图6a为根据图5所示的速度模型合成的含有层间一阶多次波和层间二阶多次波的地震记录; 对此地震记录作稀疏双曲Radon变换可以得到变换域数据, 如图6b所示, 可以看出, 经过变换, 一次波和多次波分别收敛, 图中椭圆内就是多次波收敛的能量; 在变换域将其切除然后进行反变换就可以得到去除了多次波的一次反射波(图6c); 图6d为通过变换切除的多次波。通过模型试算可以看出, 利用频域稀疏双曲Radon变换能够有效地完成有效波和多次波的波场分离, 并且具有较高的分辨率和保真度。

图5 速度模型

图6 频域稀疏Radon变换波场分离去除多次波模型试算

2.2 频域稀疏Radon变换压制随机噪声模型试算

地震勘探中的随机噪声通常被认为是白噪, Radon变换作为典型的叠加积分变换, 本身就具有压制随机噪声的作用。经过积分运算以后, 具有双曲传播特征的波场能量完成收敛, 通常认为变换域的小能量就是随机噪声叠加后的能量, 通过设置阈值对小能量进行压制, 然后进行反变换便可完成对随机噪声的压制。文中采用了硬阈值的方法, 设wj, k为Radon变换域数据, w˙j, k为用阈值方法得到的变换数据, λ 为阈值, 硬阈值函数如图7a所示。

图7b为CDP模型数据, 其采样间隔 2 ms, 采样点数2 048, 采样道数120道, 在合成CDP数据中加入了随机噪声, 信噪比为1。对模型数据进行稀疏Radon变换并对比和分析加阈值与不加阈值的去噪效果。图7c为稀疏Radon反变换得到的时间域数据, 与图7b对比可以看出稀疏Radon变换提高了信噪比。图7d为处理中去掉的随机噪声, 从图中可以看出单纯进行稀疏Radon变换不会损失有效能量。

为了进一步提高数据的信噪比, 对Radon域数据进行硬阈值处理, 将小于阈值的能量置零, 然后进行反变换, 得到如图7e的数据, 从图中可以看出数据的信噪比得到了大幅度提高。图7f为硬阈值处理后的误差数据, 经过阈值处理后浅层和大偏移距处有能量损失。

图7 频域稀疏Radon变换压制随机噪声模型试算

3 实际资料应用

通过模型试算可知稀疏双曲Radon变换对于波场分离和压制随机噪声都有较好的效果, 下面通过实际地震资料去噪来验证方法的实用性。图8a为某地区地震资料中抽取的CMP道集, 该原始数据采样间隔1 ms, 经抽稀为采样间隔2 ms; 原始记录CMP覆盖次数104, 采样点数2 501, 抽取99道, 2 000 个采样点进行处理, 验证稀疏双曲Radon变换的去噪效果。对图8a所示的CMP道集做稀疏Radon变换, 得到Radon域数据(图8b), 由图中可以清楚地看出1 000 ms和2 000 ms附近的双曲同相轴得到了较好的收敛; 对于Radon域数据进行硬阈值处理, 以压制随机噪声叠加收敛的小能量, 得到如图8c所示的处理后Radon域数据; 对于处理后的Radon域数据进行稀疏Radon反变换得到了去除噪声后的时空域CMP数据(图8d)。通过对比可以看出, 处理后的数据信噪比得到了明显的提高, 有效地去除了随机噪声, 并且在1 000 ms和2 000 ms附近的反射同相轴连续性变好, 一些若隐若现的弱反射同相轴也得到了加强, 有效地提高了信噪比, 使得反射同相轴的追踪更加容易。图8e为硬阈值处理后的误差剖面, 从图中可以看出随机噪声被去除, 并且误差剖面中没有明显的有效信号的影子, 有效反射波的波组特征未发生变化, 充分证明了该方法具有良好的保真度和保幅性, 验证了本文所用方法在实际地震资料去噪中的可行性和有效性。

图8 频域稀疏Radon变换实际资料去噪

4 结论

1) 稀疏Radon正变换能够使地震数据的同相轴按照速度差异点状收敛于不同的区域, 相对于最小二乘法消除了同相轴之间的平滑影响, 便于波场的分离。文中所用的时不变双曲积分路径更符合地下地质体反射波场的反射特点, 不但可以取得更高的分辨率, 而且可以将稀疏双曲Radon变换在频率域实现, 大大提高了计算速度, 节省了时间。模型数据测试和实际资料应用结果表明, 频域稀疏双曲Radon变换可以有效去除规则干扰波, 压制随机噪声, 提高数据信噪比, 且具有较好的保真度和保幅性。

2) 由于动校正的拉伸畸变、Radon变换近偏能量叠加及硬阈值处理的缺陷, 造成变换后浅层及大偏移距能量发散问题有待解决, 笔者正从无拉伸动校正及改用其他阈值方法进行研究改进。

The authors have declared that no competing interests exist.

参考文献
[1] Claerbout J F, Johnson A G. Extrapolation of time-dependent waveforms along their path of propagation[J]. Geophysical Journal International, 1971, 26(1): 285-293. [本文引用:1]
[2] Sacchi Mauricio D, Ulrych Tadeusz J. High-resolution velocity gathers and offset space reconstruction[J]. Geophysics, 1995, 60(4): 1169-1177. [本文引用:1]
[3] Cary P W. The simplest discrete Radon transform[C]//Expand ed Abstracts of 68th SEG Annual Internet Meeting, 1998: 1999-2002. [本文引用:1]
[4] Sacchi M D, Porsani M. Fast high-resolution parabolic Radon transforms[C]//Expand ed Abstracts of 69th SEG Annual Internet Meeting, 1999: 1477-1480. [本文引用:1]
[5] Nowak E, Imhof M. Amplitude preservation of Radon-based multiple-removal filters[J]. Geophysics, 2006, 71(5): 123-126. [本文引用:1]
[6] Yu Z, Ferguson J, McMechan G, et al. Wavelet-Radon domain dealiasing and interpolation of seismic data[J]. Geophysics, 2007, 72(2): 41-49. [本文引用:1]
[7] Wang J F. Greedy least-squares and its application in Radon transforms[C]//2009 CSPG CSEG CWLS Convention, 2009: 5-8. [本文引用:1]
[8] 王维红, 裴江云, 张剑锋. 加权抛物Radon变换叠前地震数据重建[J]. 地球物理学报, 2007, 50(3): 851-859. [本文引用:1]
[9] 熊登, 赵伟, 张剑锋. 混合域高分辨率抛物Radon变换及在衰减多次波中的应用[J]. 地球物理学报, 2009, 52(4): 1068-1077. [本文引用:1]
[10] 石颖, 王维红. 基于波动方程预测和双曲Radon变换联合压制表面多次波[J]. 地球物理学报, 2012, 55(9): 3115-3125. [本文引用:2]
[11] 鲁娥, 李庆春. 混合Radon变换地震噪声压制的应用[J]. 物探与化探, 2013, 37(4): 706-710. [本文引用:1]
[12] 张振波, 轩义华. 高分辨率抛物线拉冬变换多次波压制技术[J]. 物探与化探, 2014, 38(5): 981-988. [本文引用:1]
[13] Richard J Castle. A theory of normal moveout[J]. Geophysics, 1994, 59(6): 983-999. [本文引用:1]
[14] 石颖, 柯璇, 王维红. 一种加权抛物Radon变换地震波场重建方法[J]. 物探与化探, 2012, 36(5): 846-850. [本文引用:1]
[15] 渥·伊尔马滋. 地震资料分析——地震资料处理、反演和解释[M]. 北京: 石油工业出版社, 2006. [本文引用:1]
[16] Schonewille M A, Duijndam A J W. Parabolic Radon transform, sampling and efficiency[J]. Geophysics, 2001, 66(2): 667-678. [本文引用:1]
[17] Trad D, Ulrych T, Sacchi M D. Latest views of the sparse Radon transform[J]. Geophysics, 2003, 68: 386-399. [本文引用:1]
[18] Zhang J F, Wang W H. Multiple attenuation: an approach by incorporating multiple prediction with Radon transform[J]. Journal of Seismic Exploration, 2006, 15: 81-99. [本文引用:1]