基于多路径Radon变换的地震数据噪声压制和波型分离方法研究

唐欢欢，毛伟建*

1 中国科学院精密测量科学与技术创新研究院计算与勘探地球物理研究中心，武汉 430077 2 大地测量与地球动力学国家重点实验室,武汉 430077

0 引言

抗噪性能良好的Radon变换被广泛应用于图像分析和信号重构等方面，它是一种积分变换，在地震勘探数据处理领域，根据积分路径的不同，Radon变换可分为线型、双曲型以及抛物型(Thorson,1978;Hampson,1986;Yilmaz,1989).线性Radon变换(Linear Radon,L_Radon)又称为倾斜叠加，它是将时空域的同相轴振幅沿着偏移距方向以一定的斜率进行叠加运算，因此线性同相轴在线性Radon变换域的系数表现为一个强能量点，对应的纵横坐标分别表示同相轴的时间截距及速度信息；而随机噪声经过线性叠加后在变换域分布无规律且能量极弱，该特征有利于随机噪声压制(Turner,1990;巩向博等，2009；Sabbione et al.,2015)、面波压制(李卫忠等，1998；Wang,1999;张智等，2017)以及波场分离(Moon et al.,1986;王维红等，2006；李志娜等，2014；Lin and Sacchi,2020).但在油气勘探地震数据中，有效信号为反射波，其时距曲线为双曲型，经过线性Radon变换后，其变换系数分布并不集中，失去了变换的稀疏性特征且物理意义不明显.Hampson(1987)将双曲Radon变换(Hyperbolic Radon,H_Radon)引入到地震数据重建中，但双曲Radon变换在频率域不解耦，只能在时间域进行计算导致效率较低(Spitzer et al.,2001;石颖和王维红，2012).双曲型同相轴经过部分动校正后可近似为抛物型，而抛物Radon变换(Parabolic Radon,P_Radon)在频率域是解耦的，且在频率域具有共轭对称性使其计算效率较高，因此抛物Radon变换被广泛地应用于多次波压制(Foster and Mosher,1992;熊登等，2009；巩向博等，2014；刘菲菲等，2017)及插值重建(Sacchi and Ulrych,1995).Sacchi和Ulrych (1995)提出的利用贝叶斯参数反演求解抛物Radon变换系数方法提高了抛物Radon变换的分辨率，有利于提高多次波压制(Zhang and Wang,2006;Song et al.,2015;Sun et al.,2019)和数据重建效果(刘喜武等，2004；Wang et al.,2019).Johansen(1995)利用二阶正交多项式对叠后水平同相轴地震数据进行拟合，正交多项式系数可以表征地震数据振幅随偏移距变化情况，该方法缺少同相轴走时特征的相关参数，因此该方法针对的是水平同相轴数据，根据此特征，薛亚茹等(2014)将正交多项式引入到Radon变换中，得到具有保幅特征的2D高阶Radon变换；唐欢欢和毛伟健(2014)、唐欢欢等(2018，2020)将其推广到3D，提出基于3D高阶Radon变换的地震数据保幅重建方法.在改善积分路径方面，牛滨华(2001)提出的多项式Radon变换，将积分路径从线型或抛物型变为2阶多项式曲线，提高了浅层双曲型同相轴数据在Radon域的稀疏程度，该变换的优势在于不增加计算量的情况下可以对复杂形态同相轴数据进行稀疏表示，但不能同时对线性、抛物型及复杂形态同相轴数据进行稀疏表示；Trad等(2001)提出了将双曲Radon变换算子和线性Radon变换算子相结合的混合Radon变换，该方法缺少这两种算子的交叉项无法对复杂同相轴进行精确拟合，且在时间域进行计算效率较低；鲁娥和李庆春(2013)将其应用到地震数据去噪中来，但该变换与积分路径相关的变换参数只有线性项，导致在变换域出现混叠现象.

以上国内外学者分别从积分路径的形态、计算效率、分辨率以及保幅性方面对Radon变换进行改进，但现有的Radon变换应用于地震数据处理中仍然存在一个问题，即当线性同相轴、双曲型同相轴以及复杂形态同相轴同时存在时，不同形态同相轴在Radon域系数不能同时达到最佳稀疏状态，只有同相轴形态和积分路径一致的数据在Radon域系数才能集中在一点，而其他的系数则是发散状态，在空间上有一定的延展性，这将降低随机噪声压制效果以及波型分离精度.

针对该问题，本文提出了多路径Radon变换(Multi-path Radon,M_Radon)，这里我们把积分路径仅为线性、双曲型或抛物型的Radon变换称为单路径Radon变换，多路径Radon变换是指对每个同相轴都有线型和抛物型的积分路径，且具有两种不同类型的变换参数来与积分路径相对应，然后利用最小二乘稀疏反演方法对不同形态同相轴自适应匹配最佳积分路径，得到在Radon域同时聚焦的高分辨率系数.该变换有利于提高随机噪声的压制、面波压制以及波型分离效果.

1 方法原理

在地震数据处理中，经典抛物Radon正反变换如公式(1a)、(1b)所示(Hampson,1986;Beylkin,1987)，是广义Radon变换的一种特殊形式，公式中d为地震数据，x为检波点坐标，q为抛物型积分路径的曲率，即速度平方的倒数,m为数据在Radon域系数；将积分公式中x2变为x，如公式(2a)、(2b)所示，其中p为斜率，即速度的倒数，公式为线性Radon正反变换，也称为倾斜叠加.公式(1a)、(1b)、(2a)、(2b)分别为

(1a)

(1b)

(2a)

(2b)

通过线性Radon或抛物Radon变换，理论上可以将t-x域的线性同相轴或抛物型同相轴变换为Radon域的一个点，如图1所示，图1a为一合成的线性同相轴，64个接收道，道间距为20 m，4 ms时间采样率，512个时间采样点，合成该数据的速度为v=4000 m·s-1，经过线性Radon变换后其Radon域系数如图1b所示；图1c为一抛物型同相轴，速度为v=1800 m·s-1，经过抛物Radon变换后其Radon域系数如图1d所示.从图中可以看出，t-x域的地震数据经过合适类型的Radon变换后，其Radon域系数都各自集中在一个点附近，从图1b、d上很容易找到这两点对应的变换参数分别为p=0.00025 s·m-1，q=3.08×10-7s2·m-2；再利用公式p=1/v和q=1/v2，可以反向计算出t-x域地震数据线性同相轴和抛物型同相轴的速度分别为v=4000 m·s-1和v=1800 m·s-1，与合成模型数据的速度一致，即Radon域系数分布位置具有明确的物理意义.因此，可以利用Radon变换的稀疏性和明确的物理意义特征进行随机噪声压制、波型分离、多次波压制以及数据重建.

图1 单轴数据Radon变换效果(a)t-x域线性同相轴；(b)图(a)数据线性Radon变换结果；(c)t-x域抛物型同相轴；(d)图(c)数据抛物Radon变换结果.Fig.1 Radon transform result of single-axis data(a)Linear in-phase axis in t-x domain；(b)Linear Radon transform of data in (a)；(c)Parabolic in-phase axis in t-x domain；(d)Radon transform of data in (c).

但是，在实际地震数据中，直达波、折射波和面波的时距曲线为线性同相轴，反射波的时距曲线为双曲型(经部分NMO后可近似为抛物型)，在使用线性Radon或抛物Radon变换时，不能使所有波型数据在变换域同时达到最佳稀疏状态，如图2合成数据所示，为t-x域地震数据，包含线型和抛物型同相轴；其线性Radon变换域系数如图3a所示，抛物Radon变换域系数如图3d所示；可以看出，在使用线性Radon变换时，仅有线性同相轴的Radon域系数是聚焦、稀疏的，其他抛物型同相轴Radon域系数并不稀疏；同样地，使用抛物Radon变换时，稀疏性仅对抛物型同相轴有效；这是因为当地震数据同相轴形态和积分路径形态不一致时，需要一系列的变换参数和积分路径组成的基函数来对地震数据拟合，从而在Radon域剖面上一系列连续的p值或q值位置都有系数分布，导致系数不是稀疏的，且不能从Radon域剖面得到所有同相轴的速度信息.图3b、e分别是从Radon域系数反变换回时空域的数据；图3c、f分别为经过不同变换后的数据与原始数据的差剖面，可以看出抛物Radon变换对线性同相轴有一定的正反变换误差，特别是在同相轴的起始位置；用公式(3)来计算这两种变换的相对均方误差：

图2 线性和抛物型同相轴同时存在的地震数据Fig.2 Seismic data with linear and parabolic events

图3 图2中地震数据Radon变换效果(a)线性Radon正变换系数；(b)线性Radon反变换数据；(c)线性Radon变换误差剖面；(d)抛物Radon正变换系数；(e)抛物Radon反变换数据；(f)抛物Radon变换误差剖面.Fig.3 Radon transform of seismic data in Fig.2(a)Coefficients of linear Radon forward transform；(b)Data of inverse linear Radon transform；(c)Error of linear Radon transform；(d)Coefficients of parabolic Radon inverse transform；(e)Data of inverse parabolic Radon transform；(f)Error of parabolic Radon inverse transform.

(3)

本文提出的多路径Radon变换能解决该问题，其正反变换公式为

(4a)

(4b)

式中的积分路径既包含线性项又包含抛物型项，线性Radon变换参数p和抛物Radon变换参数q同时存在；因此，对于含有线型的直达波、折射、面波及双曲型反射波(可近似为抛物型)的地震数据，多路径Radon变换可以看作线性Radon基函数和抛物型Radon基函数的组合.选择p=0基函数，公式(4)可以对抛物型同相轴进行稀疏表达；选择q=0的基函数，公式(4)可以对线性同相轴进行稀疏表达；而p和q都不为0的基函数，表示公式(4)的积分路径为多项式，可以对时距曲线不是标准的线型或抛物型同相轴进行拟合.

和传统Radon变换一样，多路径Radon变换也将在频率域通过反演方法来求解Radon域系数m，对公式(4a)、(4b)作Fourier正变换，即得到频率域多路径Radon变换公式为

(5a)

(5b)

根据指数函数的性质，公式(5)可以进一步分解为

(6a)

(6b)

(7)

(8)

(9)

(10)

可以通过加权迭代方法(Sacchi,1997)进行近似求解.

(11)

(12)

(13)

2 模型数据测试

2.1 多路径Radon变换的稀疏表示能力分析

用多路径Radon变换对图2中的数据作正变换，线性积分路径参数p的最小值p0=0 s·m-1，采样间隔dp=2×10-5s·m-1，采样个数np=12；抛物型积分路径参数q的最小值q0=0 s2·m-2，采样间隔dq=2×10-8s2·m-2，采样个数nq=21，3次迭代求解多路径Radon域系数；系数分布如图4所示，图中有12列，分别表示积分路径参数取值为(p0,[q0,q1,…,q20])，(p1,[q0,q1,…,q20])，…,(p11,[q0,q1,…,q20])时的多路径Radon域系数，纵坐标表示时间，横坐标表示抛物型积分路径参数.从图4中可以看出，与传统Radon变换不同(图3a、d)，多路径Radon域系数可以同时聚焦，且独立分布在两个不同区域上，即图4中红色方框标识出的第1列及第11列区域.第一个区域是线性积分路径参数p=0×dp区，此时p=0，表示积分路径为抛物型，为q剖面，该区域的Radon系数有4个，分别对应了4个抛物型积分参数，能够与图2中时空域四个抛物型同相轴的起始时间、速度一一对应；第二个区域是线性积分路径参数p=10×dp区，该区域的Radon系数有1个，对应的抛物型积分路径参数q=0×dq，表示积分路径为线性，为p剖面，该位置的Radon系数与图2中时空域线性同相轴起始时间、速度一一对应，且速度可由线性积分路径参数p=10×dp计算得到；因此，通过多路径Radon不仅可以使时空域不同形态的同相轴在Radon域系数同时达到稀疏状态，而且能够很容易地识别这些系数对应的时空域同相轴形态、时间截距以及速度.这些优势有助于提高利用Radon变换方法进行随机噪声压制、面波压制以及波型分离等技术的处理效果.将多路径Radon域系数通过反变换，我们即可得到时空域地震数据，如图5a所示，变换误差剖面如图5b所示，利用公式(3)计算得到相对均方误差为1.36×10-4，误差极小，这是该变换能够用于地震数据处理的基础.

图4 多路径Radon正变换系数Fig.4 Coefficients in multi-path Radon forward transform

图5 多路径Radon反变换结果(a)反变换结果；(b)误差剖面.Fig.5 Result of multi-path Radon inverse transform(a)Result of inverse transform；(b)Error.

图6 三种变换的稀疏表示能力Fig.6 Sparse representation capability of three Radon transforms

2.2 多路径Radon变换在模型数据随机噪声压制及波型分离中的应用

多路径Radon变换的强稀疏表示能力使其在随机噪声压制和波型分离问题中具有天然优势，先以模型数据为例进行展示.图7a为干净的原始数据，包含待分离的线性和抛物型两个同相轴，同相轴在远偏移距形态接近且有重叠.地震数据共有64道，道间距为20 m，每道有300个时间采样点，采样间隔为4 ms.将采取常用的傅里叶变换(FK)、传统的线性Radon、抛物Radon以及本文提出的多路径Radon变换方法对其进行分离.图7b、c分别为该数据在线性Radon变换域和抛物Radon变换域的系数分布图，可以看出在线性Radon变换中仅有线性同相轴的系数分布是稀疏、聚焦的，如图7b中白色箭头所示，而抛物性同相轴的系数分布在空间上具有延展性，使二者的系数有部分重叠，很难将其完全分离开来；类似的情况也出现在抛物Radon变换系数中，如图7c所示；图7d、e、f分别为线性同相轴的FK谱、抛物同相轴的FK谱以及二者混合一起后的FK谱，很明显这两个轴的FK谱能量差异极大，重叠在一起后的FK谱中这两个轴的能量也有部分重叠，同样也不能将其完全分离开；图7g为多路径Radon域系数，可以看出线性同相轴和抛物同相轴的系数分布在不同区域且同时稀疏、聚焦，很容易将其分离.

图7 时空域数据在不同变换域系数(a)原始数据；(b)线性Radon变换域系数；(c)抛物Radon变换域系数；(d)线性轴FK谱；(e)抛物轴FK谱；(f)原始数据FK谱；(g)多路径Radon变换域系数.Fig.7 Coefficients of temporal-spatial data in different transform domains(a)Original data；(b)Coefficients in linear Radon transform domain；(c)Coefficients in parabolic Radon transform domain；(d)FK spectrum of linear event；(e)FK spectrum of parabolic event；(f)FK spectrum of original data；(g)Coefficients of multi-path Radon transform domain.

表1 四种变换波型分离效果比较Table 1 Comparison of data separation by 4 transforms

为了测试本文方法抗噪能力，在图7a原始干净数据上添加了一定的随机噪声，数据的信噪比为-16.09 dB，如图9a所示.根据上文三种变换的稀疏表示能力曲线来确定噪声压制的阈值，首先将图7变换域中前30%的强能量系数保留，其他系数进行压制，再进行下一步的分离.我们比较了本文提出的多路径Radon变换与传统线性Radon、抛物Radon变换方法的去噪效果.图9b、c分别为数据在线性Radon、抛物Radon变换域系数，可以看出，由于随机噪声的存在，使得线性同相轴和抛物同相轴的变换域系数更加难以区分，进而影响下一步的分离；图9d为数据在多路径Radon变换域系数，可以看出，即使存在能量较强的随机噪声，由于不同类型同相轴在该变换域系数都同时稀疏、聚焦，从而仍然能够较容易将其区分并分离；线性Radon变换方法的噪声压制及波型分离结果如图8e、f所示，抛物Radon变换方法的噪声压制及波型分离结果如图9g、h所示,本文提出的多路径Radon变换方法的噪声压制及波型分离结果如图9i、j所示，从图中可以很明显看出，线性Radon变换仅对线性同相轴的噪声压制效果较好，抛物Radon变换仅对抛物同相轴的噪声压制有效，而多路径Radon变换方法对线性同相轴和抛物同相轴的分离及噪声压制同时有效，同相轴分离完全、无振幅损失且噪声压制干净.三种方法噪声压制及分离后数据的信噪比如表2所示，本文提出的多路径Radon变换方法在去噪及分离后，信噪比达到28 dB左右，与原始-16.09 dB含噪数据相比，信噪比提高了44 dB，其效果明显优于线性Radon和抛物Radon变换.

图8 波型分离效果(a)(b)傅里叶变换波型分离后的数据；(c)(d)线性Radon变换波型分离后的数据；(e)(f)抛物Radon变换波型分离后的数据；(g)(h)多路径Radon变换波型分离后的数据.Fig.8 Waveform separations by different transforms(a)and (b)Fourier transform；(c)and (d)Linear Radon transform；(e)and (f)Parabolic Radon transform；(g)and (h)Multi-path Radon transform.

图9 去噪及分离效果(a)含噪数据;(b)线性Radon变换系数;(c)抛物Radon变换系数;(d)多路径Radon变换系数;(e)(f)线性Radon变换去噪及分离结果;(g)(h)抛物Radon变换去噪及分离结果;(i)(j)多路径Radon变换去噪及分离结果.Fig.9 Denoised and separated data by multi-path Radon transform(a)Noised data；(b)Coefficients of linear Radon transform；(c)Coefficients of parabolic Radon transform；(d)Coefficients of multi-path Radon transform；(e)and (f)Denoised and separated data by linear Radon transform；(g)and (h)Denoised and separated data by parabolic Radon transform；(i)and (j)Denoised and separated data by multi-path Radon transform.

表2 含噪数据噪声压制及波型分离效果比较Table 2 Comparison of denoised and separated data by three Radon transforms

2.3 多路径Radon变换在实际数据噪声压制及波型分离中的应用

现将该方法应用于塔里木某区域采集的实际数据面波压制测试中，如图10a所示，共71道，道间距为30 m,每道625个时间采样,采样间隔为4 ms，强烈的面波表现为线性，将近偏移距的有效信号湮没.根据上文模型数据的测试结果，可以发现在线性和抛物型同相轴同时存在的地震数据中，线性Radon变换的处理效果优于抛物Radon变换，因此在实际数据面波压制测试中，将仅提供线性Radon和本文提出的多路径Radon变换的处理效果对比.二者面波压制效果分别如图10b、c所示，从图中可以看出这两种方法都能够将线性面波压制，但多路径Radon变换方法能够得到分辨率更高的近偏移距反射波同相轴，如图中A白色椭圆区域所示；同时，对于深层较弱信号，如2.4 s深度B白色椭圆区域所示，本文方法恢复的同相轴更为连续；从去除的面波数据图10e、f中红色箭头处，可以进一步发现线性Radon变换方法去除了更多有效信号.

图10 面波压制效果(a)原始数据；(b)线性Radon变换面波压制结果；(c)多路径Radon变换面波压制结果；(d)留白区域；(e)(f)分别为(b)(c)去除的面波.Fig.10 Suppression of surface waves by different transforms(a)Original data；(b)Linear Radon transform；(c)Multi-path Radon transform；(d)Whitespace；(e)and (f)Surface waves removed by (b)and (c),respectively.

接着将多路径Radon变换方法应用于VSP数据的上下行波分离中，原始VSP数据如图11所示，该数据共201道，道间距为20 m，每道1905个时间采样，采样间隔为4 ms，该数据较为复杂，每个同相轴都不是标准的线性或双曲型.图12为三种Radon变换方法的上下行波分离效果：图12a、b分别为线性Radon变换方法分离的下行波和上行波；图12c、d分别为抛物Radon变换方法分离的下行波和上行波；图12e、f分别为多路径Radon变换方法分离的下行波和上行波.可以看出，多路径Radon变换方法的分离效果最佳，能够实现上下行波的有效分离；线性Radon变换方法次之，其振幅存在一定的损失，且在上下行波的起始交叉点位置有部分假频，如图12a、b中红色箭头位置所示；抛物Radon变换方法的效果最差，在近偏移距处倾角较大的同相轴没有得到有效恢复，且存在抛物型脚印，如图12c、d中红色虚框处所示；图13为近偏移距0～1000 m数据分离效果的放大图，可以更明显地看出多路径Radon变换方法上下行波分离的更为彻底且无信号损失，该实例进一步显示了本文方法的有效性.

图11 原始VSP数据Fig.11 Original VSP data

图13 VSP数据上下行波分离近偏移距放大效果(a)(b)线性Radon变换近偏移距分离效果；(c)(d)抛物Radon变换近偏移距分离效果；(e)(f)多路径 Radon变换近偏移距分离效果.Fig.13 Enlargement of separation of near-offset upward and downward waves in VSP data (a)and (b)Linear Radon transform；(c)and (d)Parabolic Radon transform；(e)and (f)Multi-path Radon transform.

3 结论

本文提出的基于最小二乘稀疏反演多路径Radon变换方法，能够使不同形态同相轴的Radon域系数自动分离，且保持系数的稀疏性.该方法解决了在地震数据中同时存在线性和抛物型同相轴时，传统Radon变换域系数不能同时准确地进行稀疏表达的问题.本文从变换的稀疏程度、计算效率以及处理效果方面做了详细的分析，与传统线性Radon、抛物Radon相比，在不增加额外计算量的情况下，多路径Radon变换能够更好地识别不同系数所对应的同相轴形态、时间截距以及速度.这些优势使Radon变换在随机噪声压制、面波压制以及波型分离中获得更佳的处理效果.

需要注意的是，在地质结构特别复杂区域，地震剖面中的同相轴不满足线性及双曲型走时特征时，需要采样间隔较小的p参数和q参数对其进行拟合，即p和q的采样个数增加，使多路径Radon变换算子规模成倍加大，从而导致该变换的计算效率降低，且变换域系数不再集中于某些点而是在一定范围内进行展布，失去稀疏特征.因此，综合考虑处理效果和计算效率，多路径Radon变换在同相轴有明显线性和双曲型形态的地震数据中更具优势.

致谢感谢评审专家提出的中肯意见和建议，使论文得以完善和丰富.