混合域高分辨率双曲Radon变换及其在多次波压制中的应用

巩向博,韩立国,王升超

(吉林大学地球探测科学与技术学院,吉林长春130021)



混合域高分辨率双曲Radon变换及其在多次波压制中的应用

巩向博,韩立国,王升超

(吉林大学地球探测科学与技术学院,吉林长春130021)

提出了一种快速且分辨率较高的双曲Radon变换方法。首先将地震道集数据沿时间轴进行平方变换,时间平方与截距时间平方呈线性关系,使得正、反Radon变换算子在频率域解耦;其次最优化反演求解目标函数时,在时间域内稀疏约束模型,而Radon变换算子在频率域实现,即混合域方式提高Radon域内分辨率;最后在反演计算过程中引入快速收缩迭代软阈值算法(FISTA)对稀疏反演优化问题进行求解,加快了反演收敛速度。模型和实际海上数据试算结果证明了混合域双曲Radon变换方法的计算效率和分辨率都较高。

双曲Radon变换;高分辨率;多次波压制;稀疏约束;混合域

Radon变换利用各类地震波在运动学上的差异,将其从时间空间域变换到时间慢度域加以区分,已广泛应用于地震资料处理中,诸如:波场分离[1-4]、数据重建[5-7]、多次波压制[8-11]、地震相关噪声压制[11-12]、平面波分解[13-14]、偏移成像[15]、面波提取[16]等。根据积分路径不同,Radon变换可以分为线性Radon变换、抛物Radon变换、多项式Radon变换[17]、双曲Radon变换、椭圆Radon变换、各向异性Radon变换[10]等。另外,根据时间(t)与截距时间(τ)的对应关系,Radon变换分为时不变Radon变换,如线性Radon变换、抛物Radon变换、多项式Radon变换等;或是时变Radon变换,如双曲Radon变换、椭圆Radon变换等。时不变与时变Radon变换的计算区别在于能否用傅里叶内核的Radon变换算子求解反问题。时不变Radon变换可在频率域求解,降低其反演计算时的算子矩阵维数,并有显式的矩阵运算公式;而时变Radon变换只能在时间域求解,常因变换算子矩阵过大限制了其应用。

沿着时间方向对地震数据进行傅里叶变换,使得Radon算子在频率域解耦,应用最小平方反演方法重获Radon变换的结果[18]。若要提高Radon域内的分辨率,可采用迭代更新权系数矩阵的方法,将模型矩阵方差作为稀疏权加入到最小平方反演之中,即迭代重加权方法(IRLS),这也成为工业界广泛使用的方法[19]。熊登等[8]和TRAD等[20]使用了l1-l2范数的联合反演方法,以l2范数约束变换误差,l1范数约束模型稀疏度作为正则性条件进行反演,并提出了混合域Radon变换算法,即变换算子仍为频率域算子,但稀疏权矩阵在时间域内更新,这样更好地提高了Radon域内的稀疏度与分辨率。LU[21]和ZHANG等[22]也使用了混合域Radon变换,并应用了迭代收缩的l1-l2范数联合反演方法,相比于IRLS方法,有着更好的反演收敛速度。ABBAD等[23]基于奇异值分解技术提出了频率域快速算法,即λ-f域Radon变换方法。近年来,多位学者研究了λ-f域Radon变换方法的波场分离、数据重建及多次波的压制[4,24-25]。

双曲Radon变换等由于其时变性,变换算子只能是时间域内的叠加算子。一般而言,时间域Radon变换计算量非常庞大,尤其是地震采样点数较多的时候,以地震记录道数为nx,样点数为nt,斜率个数为np为例,其一次时间域变换就需要nt·nx·np次加法,若高分辨率算法则需要多次迭代,更增加了计算的复杂性。TRAD等[26]采用了预条件的共轭梯度算法,基于迭代重加权思想求解时间域的稀疏双曲Radon变换;NG等[27]提出了一种新的时间域迭代Radon变换,首先累加最大能量的同相轴,然后反变换提取此同相轴,并在原数据中将其减去,依次迭代累加,能够得到高分辨率的时间域结果;WANG等[28]基于迭代相减方法,发展了一种贪婪反演算法,强同相轴与海上实际数据多次波压制结果都证明了这是一种稳健的方法;HU等[29]利用降阶近似法加速了时间域双曲Radon变换的计算效率。但这些方法都无法解决随着数据维数增加,双曲Radon变换计算量大幅度增加的难题。

1　混合域高分辨率双曲Radon变换

双曲Radon变换可以定义为:

(1)

式中:τ为截距时间;p为慢度(速度倒数);x为炮检距;d(t,x)为道集数据;m(τ,p)为Radon变换域内数据。

若令t′=t2,则将原数据d(t,x)变换至d′(t′,x)域,同样我们令τ′=τ2,将m(τ,p)进行坐标变换至m′(τ′,p),即原双曲Radon变换公式可变为:

(2)

这样沿着双曲线进行积分累加变为沿着抛物线进行积分累加,在新的d′(t′,x)和m′(τ′,p)坐标域内,t′与τ′呈线性变化关系,因此公式(2)可以在频率域求解。

我们定义S为时间轴平方算子;S-1为时间轴开方算子,则道集数据d(t,x)与Radon域数据m(τ,p)可分别有如下公式:

(3a)

(3b)

(3c)

(3d)

对(2)式沿时间方向进行傅里叶变换,并应用傅里叶变换的时移性质,频率域双曲Radon变换就可以写为如下的矩阵形式:

(4)

式中:m和d分别为离散数据与Radon域数据的矩阵形式;F为傅里叶正变换算子。定义LT为Radon变换频率域伴随算子:

(5)

因Radon变换的非正交性,我们建立最小平方意义下线性化反演误差的目标函数,即双曲Radon变换的频率域l2范数条件的解为:

(6)

式中:F-1为傅里叶反变换算子;L为Radon变换频率域正变换算子。

时间域内稀疏约束可以很好地提高Radon域内分辨率,我们把拉伸的时间轴道集数据作为稀疏条件进行约束反演,仍采用傅里叶内核的频率域Radon正、反变换算子,对反演误差取l2范数,对稀疏约束模型取l1范数,即求解线性反问题的l1-l2混合范数,建立如下的目标函数:

(7)

式中:λ为正则化参数,是平衡反演误差与模型稀疏化的折中参数,其值越大,Radon域内结果越稀疏。FISTA为求解线性优化问题提供了方便且快捷的解决方案,采用FISTA求解公式(7),可以提高反演的收敛速度和计算效率[33]。其中,初值为最小平方反演解,迭代更新公式如下:

(8a)

(8b)

t0=1

(8c)

当k≥0时,

(8d)

(8e)

(8f)

式中:k为当前迭代次数;soft为取软阈值算子;α为利普希茨常数,α≥max[eig(LTL)],其中max(A)为取向量A的最大值,eig(B)为取矩阵B的特征向量。经过数次迭代就可以获得混合域高分辨率的Radon变换结果。需要注意的是,此时得到的高分辨率Radon域结果,其时间轴仍为原时间轴的平方,只需将得到的结果沿时间轴开方,就得到具有实际物理意义的高分辨率双曲Radon变换结果:

(9)

类似推导,频率域双曲Radon反变换公式为:

(10)

2　多次波压制流程

相对于多次波来说,一次波有更大的正常时差差异,利用这些正常时差差异或者走时曲线曲率的不同,我们可以使用Radon变换来压制多次波,其步骤如下。

1) 将原道集数据进行时间轴平方变换。

4) 达到最大迭代次数后,将计算结果进行时间轴开方变换,得到双曲Radon域内结果。

5) 在Radon域内根据时差差异设置切除线,切掉多次波的能量。

6) 将切除后的Radon域数据定义为mp,将其进行时间轴平方变换。

7) 进行双曲Radon反变换计算。

8) 进行时间轴开方运算,输出压制多次波的道集结果。

图1　高分辨率双曲Radon变换压制多次波流程

3　效果分析

3.1模型算例

我们采用双曲同相轴来验证本文方法的正确性(图2)。设计了如图2a所示的8个截距时间不同、曲率也不同的双曲同相轴,道集数据的采样率为4ms,500个采样点,由50道组成,偏移距为0～2000m,道间距40m;图2b为Radon正变换的理想结果;沿时间轴进行平方变换,道集中小于1s的同相轴被压缩,大于1s的同相轴被拉伸,且原双曲同相轴变为抛物同相轴,时间轴位置也相应发生变化(图2c)。

将道集数据按照稀疏约束的混合域高分辨率算法进行Radon正变换,得到图2d所示结果。再沿截距时间进行开方变换,得到高分辨率双曲Ra-don域结果(图2e)。图2f是时间域高分辨率双曲Radon变换结果,因频率域非解耦,TRAD等采用预条件共轭梯度法求解混合范数反演,即在时间域内将稀疏m作为预条件正则化因子,加入到共轭梯度反演过程中[28]。其本质是时间域叠加计算,计算量受模型空间影响较大,当道集离散样点数较大时,计算效率会很低。在Intel CPU主频1.9GHz电脑运行,本文方法整个过程需要计算的时间为2.3s,而TRAD方法的计算时间为8.6s。相对于传统时间域计算结果(图2f),本文方法的结果(图2e)分辨率更高,且其计算效率也更高。

图2　模型数据测试结果a 模型数据; b 理想Radon变换结果; c 沿时间轴进行平方变换结果; d 傅里叶内核的稀疏约束混合域Radon变换结果; e 高分辨率双曲Radon变换结果; f 时间域稀疏约束Radon变换结果

3.2海上实际数据算例

图3a为海上某长偏移距CMP道集数据,偏移距道数为76道,采样率2ms,最大偏移距8km,

时间样点数为4500,最大采集时间9s。由于海水较深,海底一次反射在2800ms处,而多次波出现在5s以下。采用本文方法进行高分辨率双曲Radon变换。首先按照一次波速度对CMP道集进行动校正,一次波被拉平,而多次波由于欠校正仍有剩余时差,通过动校正减小正常时差,以降低Radon变换压制多次波出现假频的概率。图3b是CMP道集一次波速度动校正的结果。按照本文方法进行高分辨率的双曲Radon变换,其Radon域内结果如图3c所示,图3c中5s以下的多次波已经与深层反射分离,一次波视速度较大,其分布在p值较小的范围内,而多次波视速度小,分布在p值较大的范围内。图3d为采用时间域高分辨率Radon变换方法得到的Radon域数据,对比图3c

图3　实际数据测试结果a CMP道集数据; b 一次波速度动校正之后的CMP道集; c 傅里叶内核的稀疏约束混合域双曲Radon变换结果及切除线设置; d 时间域双曲Radon变换结果; e 切除多次波能量后反Radon变换的动校正CMP道集; f 多次波衰减后的CMP道集; h 衰减掉的多次波数据

和图3d中2～5s的能量可见,图3c收敛效果更好。设置图3c中所示的切除线,将多次波的能量进行切除,反Radon变换之后就得到压制多次波后的动校正道集数据(图3e)。再按照一次波的速度进行反动校正,就恢复了原CMP道集中的一次波数据(图3f),衰减掉的多次波数据如图3g所示。

此算例在Intel CPU主频1.9GHz电脑运行,整个过程耗时8.9s。而采用TRAD等的方法,即时间域高分辨率Radon变换[28],计算过程需要195.8s。图4a和图4b分别为CMP道集压制多次波前、后的速度谱。从速度谱中也可以明显地看出,时间在5～9s,叠加速度在1.5～2.0km/s的多次波能量已经得到压制,在图4b出现了叠加速度超过3km/s的一次波能量团。

图4　CMP道集压制多次波前(a)、后(b)速度谱

4　讨论

沿时间轴进行平方变换的道集数据,其双曲走时曲线变为抛物走时曲线,地震波的运动学规律发生了变化。对于小于1s的同相轴其剩余时差会减小,而对于大于1s的同相轴其剩余时差会增大,导致易出现假频现象。即对于深层反射的多次波增加了变换的难度。这种情况下我们可以选取在一次波和多次波速度之间的动校正速度,经动校正后,一次波过校正而多次波校正不足,其剩余时差的绝对值降低。

时间轴平方变换与开方变换在离散情况下存在误差,应用余弦变换可以抵消深层反射在拉伸时间轴上的畸变。我们分别采用理论模型和实际数据验证了时间轴平方与开方变换的效果,定义振幅误差为:

(11)

其中,sum为求和运算符。采用公式(11)计算得到理论模型和实际数据的振幅误差分别为:e1=0.0283,e2=0.0226。我们发现振幅的误差主要表现在高频部分,低通滤波后的数据再进行S和S-1变换,也能降低振幅误差e。也可以采取小波变换插值的S和S-1变换降低误差,但会增加计算量。另外,因为Radon变换本身就是非正交变换,其正、反变换带来的误差也会累加到最终计算结果中。

5　结论

我们提出了一种混合域高分辨率双曲Radon变换方法,该方法分别对时间坐标轴与截距时间坐标轴进行平方变换,在变换后的拉伸时间坐标轴条件下,Radon变换算子在频率域解耦。因此,我们可以采用傅里叶内核的频率域Radon正、反变换算子进行数据变换,并在时间域内稀疏约束反演变换域空间,从而提高了Radon域内数据的分辨率。采用FISTA提高线性反演的收敛速度。模型数据和实际数据试算结果证明了本文方法既保持了双曲Radon变换的高分辨率,又有较高的计算效率。本文方法也同样适用于其它时变Radon变换方法。

[1]李远钦.一种非线性Radon变换及非零偏移距VSP波场分离[J].石油物探,1994,33(3):33-39

LIYQ.Anon-linearRadontransformandnon-zerooffsetVSPwavefieldseperation[J].GeophysicalProspectingforPetroleum,1994,33(3):33-39

[2]曾有良,乐友喜,单启铜,等.基于高分辨率Radon变换的VSP波场分离方法[J].石油物探,2007,46(2):115-119

ZENGYL,YUEYX,SHANQT,etal.VSPwavefieldseparationbasedonhigh-resolutionRadontrasformation[J].GeophysicalProspectingforPetroleum,2007,46(2):115-119

[3]冯晅,张先武,刘财,等.带有多道相关的抛物线Radon变换法分离P-P、P-SV波[J].地球物理学报,2011,54(2):304-309

FENGX,ZHANGXW,LIUC,etal.SeparatingP-PandP-SVwavebyparabolicRadontransformwithmultiplecoherence[J].ChineseJournalofGeophysics,2011,54(2):304-309

[4]李志娜,李振春,王鹏.基于改进线性Radon变换的井孔地震上下行波场分离方法[J].地球物理学报,2014,57(7):2269-2277

LIZN,LIZC,WANGP.WavefieldseparationbyamodifiedlinearRadontransforminboreholeseismic[J].ChineseJournalofGeophysics,2014,57(7):2269-2277

[5]张红梅,刘洪.基于稀疏离散τ-p变换的非均匀地震道重建[J].石油物探,2006,45(2):141-145

ZHANGHM,LIUH.Sparsenessdiscreteτ-ptransforminirregularseismictracereconstruction[J].GeophysicalProspectingforPetroleum,2006,45(2):141-145

[6]王维红,裴江云,张剑锋.加权抛物Radon变换叠前地震数据重建[J].地球物理学报,2007,50(3):851-859

WANGWH,PEIJY,ZHANGJF.PrestackseismicdatareconstructionusingweightedparabolicRadontransform[J].ChineseJournalofGeophysics,2007,50(3):851-859

[7]李晶晶,孙成禹,谢俊法,等.相对保幅的抛物线Radon变换法地震道重建[J].石油物探,2014,53(2):181-187

LIJJ,SUNCY,XIEJF,etal.Seismictracereconstructionbyrelativeamplitude-preservedparabolicRadontransform[J].GeophysicalProspectingforPetroleum,2014,53(2):181-187

[8]熊登,赵伟,张剑锋.混合域高分辨率抛物Radon变换及在衰减多次波中的应用[J].地球物理学报,2009,52(4):1068-1077

XIONGD,ZHAOW,ZHANGJF.Hybrid-domainhigh-resolutionparabolicRadontransformanditsapplicationtodemultiple[J].ChineseJournalofGeophysics,2009,52(4):1068-1077

[9]石颖,王维红.基于波动方程预测和双曲Radon变换联合压制表面多次波[J].地球物理学报,2012,55(9):3115-3125

SHIY,WANGWH.Surface-relatedmultiplesuppressionapproachbycombiningwaveequationpredictionandhyperbolicRadontransform[J].ChineseJournalofGeophysics,2012,55(9):3115-3125

[10]巩向博,韩立国,李洪建.各向异性Radon变换及其在多次波压制中的应用[J].地球物理学报,2014,57(9):2928-2936

GONGXB,HANLG,LIHJ.AnisotropicRadontransformanditsapplicationtodemultiple[J].ChineseJournalofGeophysics,2014,57(9):2928-2936

[11]巩向博,韩立国,王恩利,等.压制噪声的高分辨率Radon变换法[J].吉林大学学报(地球科学版),2009,39(1):152-157

GONGXB,HANLG,WANGEL,etal.DenoisingviahighresolutionRadontransform[J].JournalofJilinUniversity(EarthScienceEdition),2009,39(1):152-157

[12]薛昭,董良国,单联瑜.Radon变换去噪方法的保幅性理论分析[J].石油地球物理勘探,2012,47(6):858-867

XUEZ,DONGLG,SHANLY.AmplitudepreservationtheoreticalanalysisofRadontransformsde-nosingmethod[J].OilGeophysicalProspecting,2012,47(6):858-867

[13]巩向博,吕庆田,韩立国,等.起伏地表地震波场角度域照明分析[J].吉林大学学报(地球科学版),2013,43(2):610-615

GONGXB,LVQT,HANLG,etal.Methodofangle-domainilluminationanalysisonruggedsurface[J].JournalofJilinUniversity(EarthScienceEdition),2013,43(2):610-615

[14]王雄文,王华忠.基于压缩感知的高分辨率平面波分解方法研究[J].地球物理学报,2014,57(9):2946-2960

WANGXW,WANGHZ.Aresearchofhigh-resolutionplane-wavedecompositionbasedoncompressedsensing[J].ChineseJournalofGeophysics,2014,57(9):2946-2960

[15]黄新武,吴律,宋炜.拉冬投影法三维叠前深度偏移[J].地球物理学报,2004,47(2):321-326

HUANGXW,WUL,SONGW.3Dprestackdepthmigrationwithradonprojection[J].ChineseJournalofGeophysics,2004,47(2):321-326

[16]潘冬明,胡明顺,崔若飞,等.基于拉东变换的瑞雷面波频散分析与应用[J].地球物理学报,2010,53(11):2760-2766

PANDM,HUMS,CUIRF,etal.DispersionanalysisofRayleighsurfacewavesandapplicationbasedonRadontransform[J].ChineseJournalofGeophysics,2010,53(11):2760-2766

[17]牛滨华,孙春岩.多项式Radon变换[J].地球物理学报,2001,44(2):263-271

NIUBH,SUNCY.PolynomialRadontransform[J].ChineseJournalofGeophysics,2001,44(2):263-271

[18]BEYLIKNG.DiscreteRadontransform[J].IEEETransactionsonAcoustics,Speech,andSignalProcessing,1987,35(2):162-172

[19]SACCHIMD,UlrychTJ.High-resolutionvelocitygathersandoffsetspacereconstruction[J].Geophysics,1995,60(4):1169-1177

[20]TRADDO,UlrychTJ,SacchiMD.LatestviewsofthesparseRadontransform[J].Geophysics,2003,68(1):386-399

[21]LUW.Anacceleratedsparsetime-invariantRadontransforminthemixedfrequency-timedomainbasedoniterative2Dmodelshrinkage[J].Geophysics,2013,78(4):V147-V155

[22]ZHANGY,LUW.2Dand3Dprestackseismicdataregularizationusinganacceleratedsparsetime-invariantRadontransform[J].Geophysics,2014,79(5):V165-V177

[23]ABBADB,URSINB,PORSANIMJ.Afast,modifiedparabolicRadontransform[J].Geophysics,2011,76(1):V11-V24

[24]石颖,张振,李婷婷,等.λ-f域加权抛物Radon变换地震数据重建方法研究[J].地球物理学进展,2014,29(4):1752-1757

SHIY,ZHANGZ,LITT,etal.Investigationonseismicdatareconstructionbyλ-fdomainweightedparabolicRadontransformmethod[J].ProgressinGeophysics,2014,29(4):1752-1757

[25]王维红,张振,石颖,等.λ-f域抛物Radon变换多次波压制方法[J].东北石油大学学报,2015,39(1):17-22

WANGWH,ZHANGZ,SHIY,etal.InvestigationofmultiplesuppressionbyparabolicRadontransforminλ-fdomain[J].JournalofNortheastPetroleumUniversity,2015,39(1):17-22

[26]TRADDO,ULRYCHTJ,SACCHIMD.Accurateinterpolationwithhigh-resolutiontime-variantRadontransforms[J].Geophysics,2002,67(2):644-656

[27]NGM,PERZM.HighresolutionRadontransforminthet-xdomainusing“intelligent”prioritizationoftheGauss-Seidelestimationsequence[J].ExpandedAbstractsof74thAnnualInternatSEGMtg,2004:2160-2163

[28]WANGJF,NGM,PERZM.SeismicdatainterpolationbygreedylocalRadontransform[J].Geophysics,2010,75(6):WB225-WB234

[29]HUJ,FOMELS,DEMANETL,etal.AfastbutterflyalgorithmforgeneralizedRadontransforms[J].Geophysics,2013,78(4):U41-U51

[30]FOSTERDJ,MOSHERCC.SuppressionofmultiplereflectionsusingtheRadontransform[J].Geophysics,1992,57(3):386-395

[31]BICKELSH.FocusingaspectsofthehyperbolicRadontransform[J].Geophysics,2000,65(2):652-655

[32]ZHANGY,LUW.A3Dconstrainedhybridlinear-pseudohyperbolicRadontransform[J].ExpandedAbstractsof84thAnnualInternatSEGMtg,2014:4350-4354

[33]BECKA,TEBOULLEM.Afastiterativeshrinkage-thresholdingalgorithmforlinearinverseproblems[J].SIAMJournalonImagingSciences,2009,2(1):183-202

[34]YILMAZO.Velocity-stackprocessing[J].GeophysicalProspecting,1989,37(4):357-382

[35]THORSONJR,CLAERBOUTJF.Velocity-stackandslant-stackstochasticinversion[J].Geophysics,1985,50(12):2727-2741

[36]WANGYH.MultipleattenuationcopingwiththespatialtruncationeffectintheRadontransformdomain[J].GeophysicalProspecting,2003,51(1):75-87

[37]SCHONEWILLEMA,AARONPA.Applicationsoftime-domainhigh-resolutionRadondemultiple[J].ExpandedAbstractsof77thAnnualInternatSEGMtg,2007:2565-2569

(编辑：陈杰)

High-resolution hyperbolic Radon transform and its application in multiple suppression

GONG Xiangbo,HAN Liguo,WANG Shengchao

(CollegeofGeo-ExplorationScienceandTechnology,JilinUniversity,Changchun130021,China)

We propose a fast and high-resolution hyperbolic Radon transform method.Firstly,the square transform is implemented on the seismic gather data along the time axis.Since the time square (t2) and intercept time square (τ2) has a linear relationship,the decoupling can be implemented for forward and inverse Radon transform algorithm in the frequency domain.Secondly,the sparse inverse problem is solved in the time domain and the Radon transforms are implemented in the frequency domain,so that the resolution of Radon transform is improved in the mixed frequency-time domain.Thirdly,based on the fast iterative shrinkage-thresholding algorithm (FISTA),the sparse inversion can be accelerated.Synthetic and real offshore seismic data testing shows that the new hyperbolic Radon transform is robust and effective.

hyperbolic Radon transform,high-resolution,multiple suppression,sparse constraint,mixed domain

2016-01-04;改回日期：2016-06-22。

巩向博(1981—),男,博士,副教授,主要从事勘探地震数据处理方法研究。

国家自然科学基金项目(41674119)和国家高技术研究发展计划(863计划)(2014AA06A605)联合资助。

P631

A

1000-1441(2016)05-0711-08

10.3969/j.issn.1000-1441.2016.05.010

This research is financially supported by the National Natural Science Foundation of China (Grant No.41674119) and the National High-tech R&D Program (863 Program) (Grant No.2014AA06A605).