高阶3D Radon变换及其数据重建应用

薛亚茹, 陈健升, 钱步仁

(中国石油大学地球物理与信息工程学院,北京 102249)



高阶3D Radon变换及其数据重建应用

薛亚茹, 陈健升, 钱步仁

(中国石油大学地球物理与信息工程学院,北京 102249)

针对高分辨率Radon变换无法较好地保持地震波振幅变化信息问题,在高阶Radon变换和3D Radon变换基础上,根据地震波传播的横向连续性提出高阶3D Radon变换。该方法首先将正交多项式拟合延拓至正交多项式面拟合,模拟平面波沿不同传播方向的振幅变化特性。与3D Radon变换结合,实现地震波在不同传播路径下的正交多项式面拟合,构成高阶3D Radon变换。该方法既考虑了地震波的传播路径,又考虑了地震波传播振幅变化,因此具有保幅以及高分辨率的特点。应用于地震数据重建表明,该方法具有良好的重建效果,并具有一定的抗噪性,且能够较好地保留地震波在不同方向上的振幅变化特性。

Radon变换; 正交多项式变换; 保幅; 数据重建

引用格式:薛亚茹,陈健升,钱步仁. 高阶3D Radon变换及其数据重建应用[J].中国石油大学学报(自然科学版),2016,40(3):69-76.

XUE Yaru, CHEN Jiansheng, QIAN Buren. High order 3D Radon transform and its application in data reconstruction[J].Journal of China University of Petroleum(Edition of Natural Science), 2016,40(3):69-76.

Radon变换在地震数据重建、多次波压制等方面得到广泛应用[1-2],其基本思想是将不同速度的同相轴映射在Radon域不同速度参数上,从而实现不同传播速度波的分离。由于有限的地震采集系统,同相轴变换到Radon域后并没有映射到单独的一个速度参数上,而是具有一定的速度参数范围,降低了同相轴在Radon域的分辨率。为了提高Radon变换分辨率,诸多学者们提出基于最小熵和Bayes理论的反演方法及其各种改进方法[2-4],将不同速度信号的先验知识与反演矩阵结合,实现了高分辨率Radon变换,并能够对采集系统外数据进行较好的外推。高分辨率Radon变换提高了Radon域参数的分辨率,但是无法保持同相轴的振幅变化特性,这是因为Radon变换的实质是将地震数据分解成一系列常振幅同相轴的叠加。如果分辨率达到最高,同相轴变换到Radon域仅用一个曲率参数表征,那么反变换后该同相轴振幅是一个常数,地震波传播振幅随炮检距变化的信息就会丢失,因此Radon变换分辨率与振幅变化信息的保持性相互制约。针对Radon变换分辨率与保幅问题,薛亚茹等[5]将正交多项式变换与Radon 变换结合,提出了高阶Radon变换。该变换利用正交多项式变换描述同相轴振幅变化信息[6],利用Radon变换描述同相轴路径参数,在常规Radon变换仅描述同相轴振幅叠加特性的基础上,增加了同相轴振幅变化的斜率和曲率信息,这些信息不仅提高了Radon变换的分辨率,而且保留了地震波传播振幅变化特性。上述高阶Radon变换仅考虑了地震波在一个方向上的传播特性,随着3D采集系统的发展,3D Radon变换处理高密度、宽方位地震数据与2D方法相比具有明显优势。Donati等[7]提出了3D tau-p地震数据重建;Hugonnet[8]提出了3D高分辨率抛物Radon变换滤波器,该滤波器利用椭圆模型处理宽方位角地震数据;Zhang[9]提出3D高分辨率tau-p变换,引入新的算法提高计算效率。这些3D Radon变换仍是以常振幅同相轴作为变换基函数。唐欢欢等[10]在此基础上将高阶Radon变换与3D变换结合,提出高阶3D变换,对两个空间方向分别做多项式拟合。数据重建是地震资料处理的一个关键环节,目前常用的高维数据重建主要有基于Fourier变换的方法以及奇异谱分解等方法[11-12]。这些方法均假设同相轴是线性或拟线性,利用预测算子实现插值。由于实际同相轴往往不满足该条件,因此基于Fourier变换的重建方法通常采用分窗处理，而Radon变换考虑了地震波的传播路径特点,在这一点上优于基于Fourier变换的方法。笔者针对3维地震数据体,将沿测线的多项式拟合方法延拓到针对3维采集的多项式曲面拟合,进一步完善高阶3D Radon变换,将该方法应用于数据重建。

1　高阶2D Radon变换基本原理

1.1正交多项式变换

通常地震数据d(t,x)在某时刻振幅沿偏移距变化可以用正交多项式拟合[6],即

(1)

式中,pj(x)为j阶代数多项式,{pj(x),j=0,1,…，N}为由N+1个偏移距坐标x确定的单位正交多项式集;c(t,j)为t时刻振幅在j阶正交多项式的分解系数,称之为正交多项式系数谱。利用最小二乘法求得系数c(t,j):

(2)

在正交多项式系数谱中,

(3)

它表示水平方向各地震道振幅的叠加;c(t,1)是地震道数据与一阶正交多项式的内积运算,表示了各地震道数据振幅变化的斜率信息,同理c(t,2)表示各地震道振幅变化的曲率信息。对于经过动校正后的地震同相轴,前几阶正交多项式系数包含了振幅的主要能量信息,高阶系数表征了振幅变化的细节和噪声,因此滤除高阶系数可以较好地压制噪声。Johansen等[6]提出该方法并实现数据去噪和AVO特征剖面分析。

1.2高阶Radon变换

Radon变换根据积分路径不同,可以分为线性Radon变换、抛物Radon变换和双曲Radon变换。以抛物Radon变换为例,给出高阶抛物Radon变换的数学模型及物理解释,同理可推导出高阶线性和双曲Radon变换。

Radon变换对的时域表达式为

(4)

式中,d(t,x)为地震数据;m(τ,q)为地震数据在Radon域的参数;q为积分路径的曲率参数。式(4)对应的频率域矩阵形式为

(5)

其中

(6)

式中,L为Radon算子;LH为L共轭转置矩阵。L算子中每一列表示不同偏移距在某一曲率路径上的地震道相移,其幅度均是1,因此L算子在重构数据时各地震道振幅随炮检距是不变的,这也是高分辨率Radon变换无法保留振幅随炮检距变化信息的主要原因。

对比式(2)和式(4),可以看到抛物Radon变换沿不同曲率的抛物路径对地震数据叠加,获得了数据振幅叠加特性和路径曲率参数;正交多项式变换沿水平方向对地震道数据进行正交多项式变换,获得了地震道振幅的叠加、斜率和曲率等多项式特性,而没有同相轴的路径方向参数。将Radon变换的方向特性和正交多项式变换的振幅特性结合,得到同时表征路径和振幅特性的变换,即高阶Radon变换

(7)

式中,mj(τ,q)为在时刻τ沿曲率q的抛物路径上地震数据振幅的第j阶正交多项式分解系数,即实现地震数据任意时刻沿不同曲率路径的正交多项式变换;m0(τ,q)为不同时刻沿不同曲率路径对地震道的振幅叠加,与常规Radon变换仅有一个比例因子之差。m1(τ,q)是不同曲率地震道振幅变化的斜率,称之为斜率剖面;同理 m2(τ,q)描述地震道的振幅变化曲率特性,称之为曲率剖面,依次类推可以获得其他更高阶的多项式特性剖面[13]。式(7)是对地震道数据振幅的分解过程,根据{pj(x),j=0,1,…，N}的单位正交特性,易得地震道的重构表示为

(8)

上述表达式清楚地表明高阶Radon变换将地震道数据表示成一系列不同曲率路径的多项式振幅的叠加,与实际地震道振幅变化更接近。以二阶Radon变换为例,式(8)的矩阵形式为

(9)

式中,L0、L1、L2分别称为零阶、一阶和二阶Radon算子,其中L0与式(6)定义相同,L1、L2是将L0算子的振幅分别用一阶和二阶多项式加权得到,从而获得振幅随炮检距变化的变换基函数。

目前高阶Radon变换在数据重建和多次波压制等领域得到一定的应用[5,14],但是在实际处理时发现,上述Radon变换仅能处理2维数据,对3维数据只能分剖面处理,因此有必要研究高阶3D Radon变换。

2　高阶3D Radon变换原理

2.1多项式曲面拟合

考虑到实际3维地质构造,为了更好地模拟波场振幅传播特点,本文中将常规多项式拟合延拓至多项式面拟合,即分别沿纵、横测线方向用正交多项式拟合地震波的振幅变化。设沿纵、横测线方向地震道振幅分别用正交多项式{pi(x),i=0,1,…,N}和{pj(y),j=0,1,…,M}拟合,则多项式{pi(x),i=0,1,…,N}和{pj(y),j=0,1,…,M}张量形成多项式曲面{pij(x,y),i=0,1,…,N,j=0,1,…,M},每一个曲面pij(x,y)由相互正交的多项式张量形成,因此{pij(x,y),i=0,1,…,N,j=0,1,…,M}两两相互正交,形成正交多项式特征曲面集合。设3维地震数据d(t,x,y),不同时刻的时间切片用正交多项式特征曲面拟合得到

(10)

式中,cij(t)为t时刻振幅随偏移距沿x方向的i阶和沿y方向的j阶正交多项式分解系数,这些系数构成了正交多项式谱。式(10)的矩阵形式表示为

D=PxCPy.

(11)

式中,Px和Py分别表示正交多项式集{pi(x),i=0,1,…,N}和{pj(y),j=0,1,…,M}构成的多项式矩阵,其张量系数cij(t)形成了NM个正交多项式系数谱,N和M分别表示两个方向的拟合阶次。

2.23D Radon变换

3D Radon变换同时沿纵、横测线方向对地震数据进行叠加,其变换对表示为[8]

(12)

式中,d(t,x,y)为3D地震数据;qx、qy为纵、横测线的路径曲率参数。式(12)的频率域矩阵表示为

(13)

式中,Lx、Ly为路径积分函数。由于Lx、Ly不是正交矩阵,因此通常采用最小二乘解

(14)

上述最小二乘解是在与原数据均方误差最小准则下得到的,受有限采集空间的影响,最小二乘解的Radon参数聚焦性较弱,其高分辨率解可以通过添加稀疏约束，利用加权迭代或迭代收缩算法求得[2,9]。

2.3高阶3D Radon变换

类似高阶2D Radon变换式(9),将式(13)中3D Radon变换算子Lx、Ly与式(11)中正交多项式集Px、Py结合,分别构成两个测线方向的高阶Radon算子,用矩阵表示为

(15)

式中,Lx0、Lx1、Lx2为经过各阶正交多项式{pi(x),i=0,1,…,N}加权的x方向Radon算子,同理得到y方向高阶Radon算子;mij为数据D在x方向的i阶、y方向的j阶的Radon参数。高阶3DRadon变换将3维地震数据沿某一曲率的切片数据用各阶正交多项式面拟合,涵盖了振幅变化信息,具有较好的保幅性质,但同时高阶3DRadon变换的反演参数规模是纵、横测线多项式拟合阶次的乘积,比2DRadon反演参数大大增加,因此提高计算效率是高阶3DRadon变换的研究重点。

2.4高阶高分辨Radon变换

高阶3D Radon变换实质由正交多项式变换和3D Radon变换构成,其中正交多项式变换是一个完备正交变换,3D Radon变换是一个欠定问题,这两个变换的结合构成超完备的变换空间,属于信号过完备稀疏表征的范畴。该类问题是近几年信号处理领域的研究热点,过完备表征是指信号变换空间的维数远远大于信号的维数,从过完备空间中寻找最佳匹配信号的基函数构成子空间,从而获得信号的稀疏表征。高阶3D Radon变换就是从不同曲率的正交多项式面空间中寻找与实际数据更加匹配的基函数,从而获得高分辨结果,因此如何寻找到最佳匹配的基函数是高分辨率研究的主要难点。目前匹配追踪和基追踪是解决该类问题常用的两种方法[15-16]。本文中采用匹配追踪方法实现高分辨率高阶3D Radon变换。

匹配追踪算法是一种迭代算法,其基本原理是每次迭代时找出与原始数据最佳匹配的基函数,在该子空间进行变换;对残差数据重复上述匹配步骤,直到匹配误差最小,匹配的所有基函数构成了反演空间,其解也是稀疏的,因此在匹配追踪算法中关键是如何找到匹配基函数。

在高阶3D Radon变换中,每一个正交多项式面是相互正交的,根据Pasval定理,正交多项式系数矩阵表明了原数据在正交多项式域的能量分布,即

(16)

式(16)表明在不同曲率的路径上同相轴的能量分布情况。当曲率参数与实际数据曲率匹配时,正交多项式谱集中在几个低阶系数,这时上述能量谱出现局部极大值;当曲率参数偏离实际数据曲率时,正交多项式谱向高阶蔓延,导致能量聚焦性变弱,上述低阶能量系数谱不能完全表征信号能量,因此其值比真实值偏小。通过上述分析,在正交多项式能量谱中根据局部极值可以找到与数据同相轴相匹配的曲率参数,由这些参数构成Radon子空间,实现该参数空间的正交多项式拟合,即实现了一个匹配追踪,通过对残差数据迭代进行上述匹配,直至误差达到精度要求,最终得到高阶3D Radon变换的高分辨率解。综上所述,实现高分辨率Radon变换流程如下:

(1) 用原始数据D作为初始化残差数据Dresi,初始化Radon参数M为零。

(2) 重复以下步骤直至满足停止条件。

②利用式(16)计算Radon能量分布,根据局部极大值选取Radon域子空间Si;

最后所有的子空间Si构成了高阶稀疏3DRadon变换。

3　理论模型与实际数据处理

3.1理论模型实验分析

高阶3D Radon变换是3D Radon变换的改进与发展,该方法能够更好地保留地震数据的振幅变化特性。本文中以地震数据去噪重建为例验证本文方法的有效性。

图1(a)和3(a)是褶积模型合成3维地震记录横、纵测线剖面,子波主频为40 Hz,道间距25 m,测线间距25 m,4个同相轴振幅随偏移距发生变化,用于模拟地震波在不同方向传播时的振幅变化特性。在原始合成记录加入随机噪声,如图1(d)和3(d)所示,信噪比为0 dB。为验证本方法去噪重建效果,沿横测线方向每10条测线保留1条,测线间距达到250 m,最后仅保留了6条测线数据,如图1(d)所示。采用高阶3D Radon变换和高阶2D Radon变换对横测线剖面的重建结果如图1(b)、1(c)所示,图1(e)、1(f)为其对应的重建误差剖面。通过对比可以看到,高阶2D Radon变换对稀疏采样数据能够较好地实现数据内插,振幅变化也有体现,但是其外推数据误差较大,同时由于没有纵测线方向的约束,仍有部分噪声未被压制,重建数据信噪比为4.5 dB。高阶3D Radon变换噪声压制明显,如图1(b)所示;残差数据图1(e)中没有同相轴痕迹,体现了较好的去噪重建效果,重建数据信噪比为19 dB。

图2为横测线地震数据的F-K谱。图2(a)为原始数据的F-K谱,图2(d)是测线间距250 m时的F-K谱,可以看到大量的假频。高阶3D Radon及2D Radon重建数据的F-K谱分别如图2(b)、2(c)所示,图2(e)、2(f)分别为其与真实频谱之间的误差。通过对比可以看到,3D变换方法几乎重建了原始数据,而且噪声压制效果较好;在高阶2D Radon方法的残差谱中,仍残留了同相轴的部分信息,这与图1(f)中仍有同相轴痕迹是一致的。

图1　横测线数据重建方法比较Fig.1　Crossline data reconstruction comparison of different methods

图2　横测线数据重建F-K谱比较Fig.2　F-K spectrum comparison

图3(b)、3(c)为2D和3D方法对某一缺失测线的重建结果剖面,图3(e)、3(f)为其对应的误差剖面,测线缺失前的模拟加噪数据剖面如图3(d)所示。可以看到,高阶3D Radon方法不仅恢复了同相轴振幅连续变化特征,同时较好地压制了噪声;高阶2D Radon方法虽然采用多项式拟合恢复了同相轴振幅,但是由于缺失了一个方向的约束导致残留了较多的噪声。

本文中也对非规则缺失数据进行了重建实验,将图4(a)中横测线原始数据随机缺失45道,如图4(d)所示,最大测线间距300 m。图4(b)、4(c)分别为高阶3D Radon及高阶2D Radon重建结果,其对应的重建误差分别如图4(e)、4(f)所示。可以看出高阶3D Radon变换对非规则缺失数据仍表现出较好的重建效果,2D方法中仍有部分随机噪声未被压制。

图4　横测线非规则缺失数据重建方法比较Fig.4　Crossline irregular data reconstruction comparison of different methods

3.2实际数据处理分析

传统数据采集时两炮之间时间间隔要足够大,以保证炮记录不互相混叠,所以通常炮点采样密度低。以某一海洋拖缆数据为例,实现高阶3D Radon变换炮间插值。图5是5个共炮点记录,炮间距25 m,排列成3维数据体,沿测线方向构成共炮点记录,沿炮集方向构成了共偏移距道集。由于两个方向同相轴的路径特点不一样,在实际处理中采用了混合Radon变换策略,即对共炮点记录采用双曲Radon变换,共偏移距剖面采用了线性Radon变换。利用第1、3、5共炮点记录插值重构第2、4共炮点记录结果如图6所示,高阶3D Radon方法重建信噪比是21.5 dB,2D Radon方法重建信噪比是14.6 dB。

图5　某海洋5炮数据记录Fig.5　Five shot gathers of marine data

图6　共炮点记录重建结果Fig.6　Reconstructed results of common shot gathers

从图7所示重构误差可以看到高阶3D Radon方法对同相轴的振幅信息保持更好,能量损失更小。

图7　共炮点记录重建误差Fig.7　Error of reconstructed results of common shot gathers

4　结　论

地震数据去噪重建是地震数据处理中的关键环节。基于实际3维数据在不同方向上的连续性,提出高阶3D Radon变换,该方法将正交多项式面拟合嵌入高阶Radon变换,不仅考虑了常规Radon变换的叠加特性,还增加了振幅的斜率和曲率信息,更加接近实际数据的真实情况。针对反演参数过多的问题,采用匹配追踪算法,降低反演空间维数,并且仅采用最小二乘方法求解Radon参数,避免了加权迭代等方法中反复修正反演矩阵的过程,算法简单易实现。模拟数据测试以及实际数据处理表明该方法通过对3维数据的两个方向添加约束,对于大距离缺失数据重建效果较好。高阶3D Radon的优势是提高分辨率的同时可实现振幅保真,该方法也可推广至多次波压制等相关应用。

[1]HAMPSON D. Inverse velocity stacking for multiple elimination[J]. Journal of Canada Social Exploration Geophysics,1986,22:44-55.

[2]SACCHI M, ULRYCH T. High-resolution velocity gathers and offset space reconstruction[J]. Geophysics,1995,60(4):1169-1177.

[3]THORSON R, CLAERBOUT J. Velocity-stack and slant-stack stochastic inversion[J]. Geophysics,1985,50(12):2727-2741.

[4]TRAD T, ULRYCH T, SACCHI M. Latest views of the sparse Radon transform[J]. Geophysics,2003,68(1):386-399.

[5]XUE Y R, MA J T, CHEN X H. High-order sparse Radon transform for AVO-preserving data reconstruction[J]. Geophysics,2014,79(2):V13-V22.

[6]JOHANSEN T A, BRULAN L, LUTRO J. Tracking the AVO by using orthogonal polynomials[J]. Geophysical Prospecting,1995,43:245-261

[7]DONATI M S, MARTIN N W, STEWART R R. 3D τ‐p filtering: proceedings of the 58th EAGE Annual International Conference and Exhibition, 1996[C/OL]. [2015-11-10]. http://earthdoc.eage.org/ publication/publicationdetailspublication=13553.

[8]HUGONNET P, BOELLE J, MIHOUB M et al. High resolution 3D parabolic Radon filtering: proceedings of the 79th SEG Annual International Meeting,Houston,2009[C]. Houston: Society of Exploration Geophysicists, c2008.

[9]ZHANG Y Q, LU W K. 2D and 3D prestack seismic data regularization using an accelerated sparse time-invariant Radon transform[J]. Geophysics, 2014,79(1):V165-V177.

[10]唐欢欢,毛伟建. 3D高阶抛物Radon变换地震数据保幅重建[J].地球物理学报,2014,57(9):2918-2927.

TANG Huanhuan, MAO Weijian. Amplitude-preserved seismic ata reconstruction by 3D high order parabolic Radon transform[J]. Chinese Journal of Geophysics, 2014,57(9):2918-2927.

[11]黄建平,李闯,李国磊,等.基于奇异谱分析的联合去噪及规则化方法[J].地球物理学进展, 2014,29(4):1666-1671.

HUANG Jianping, LI Chuang, LI Guolei, et al. Simultaneous seisimic data denoising and regularization method based on singlular spectrum analysis [J]. Progress in Geophysics,2014,29(4):1666-1671.

[12]魏小强,雷秀丽,马庆珍.基于多道奇异谱分析的三维地震数据规则化方法[J].地球石油物理勘探,2014,49(5):846-851.

WEI Xiaoqiang, LEI Xiuli, MA Qingzhen. 3D seismic data regularization based on multi-channel singular spectrum analysis [J]. Oil Geophysical Prospection,2014,49(5):846-851.

[13]薛亚茹,唐欢欢,陈小宏.高阶高分辨率Radon变换地震数据重建方法[J].石油地球物理勘探,2014,49(1):95-100.

XUE Yaru, TANG Huanhuan, CHEN Xiaohong. Seisimic data reconstruction based on high order high resolution Radon transform[J]. Oil Geophysical Prospection,2014,49(1):95-100.

[14]薛亚茹,杨静,钱步仁.基于高阶稀疏Radon变换的预测多次波自适应相减方法[J].中国石油大学学报(自然科学版),2015,39(1):43-49.

XUE Yaru,YANG Jing,QIAN Buren. Multiples prediction and adaptive subtraction with high-order sparse Radon transform [J]. Journal of China University of Petroleum(Edition of Natural Science),2015,39(1):43-49.

[15]MALLAT S G, ZHANG Z F. Matching pursuit with time-frequency dictionaries [J]. IEEE Transaction on Signal Processing, 1993,41(12):3397-3415.

[16]CHEN S S, DONOHO D L, SAUNDERS M A. Atomic decomposition by basis pursuit[J]. SIAM Review, 2001,43:129-159.

(编辑修荣荣)

High order 3D Radon transform and its application in data reconstruction

XUE Yaru, CHEN Jiansheng, QIAN Buren

(CollegeofGeophysicsandInformationEngineeringinChinaUniversityofPetroleum，Beijing102249，China)

The high resolution Radon transform cannot preserve the event amplitude information due to the incompleteness of the transform basis functions. Based on high order 2D Radon transform and 3D Radon transform, a high order 3D Radon transform is proposed to keep the amplitude information continuously along both inline and crossline directions. The method first expands the orthogonal polynomial transform to 2 dimensions and then embeds it to 3D Radon transform to construct high order 3D Radon transform. The proposed method integrates the gradient and curvature information in the event amplitude variations into the 3D Radon transform. With these additional properties, the directional variations in the eventsamplitude can be modeled continuously in the transformation, which improves the resolution of Radon transform. The experiments of synthetic and field data show robust results in denoising and reconstruction.

Radon transform; orthogonal polynomial transform; amplitude-preservation; data reconstruction

2015-07-03

国家自然科学基金项目( 41204095)

薛亚茹(1972-) ，女，副教授，博士，研究方向为信号处理。E-mail:xueyaru@cup.edu.cn。

1673-5005(2016)03-0069-08doi:10.3969/j.issn.1673-5005.2016.03.009

TE 132

A