基于各向異性Radon變換的疊前地震數(shù)據(jù)重建

王升超,韓立國,鞏向博

(吉林大學地球探測科學與技術(shù)學院,吉林長春130012)

?

基于各向異性Radon變換的疊前地震數(shù)據(jù)重建

王升超,韓立國,鞏向博

(吉林大學地球探測科學與技術(shù)學院,吉林長春130012)

基于常規(guī)雙曲Radon變換進行疊前地震數(shù)據(jù)重建時,即使分辨率很高也難以得到準確的中遠偏移距數(shù)據(jù),因為地震信號同相軸在偏移距較遠時已經(jīng)不再滿足雙曲線假設。提出基于各向異性Radon變換,即由偏移距、慢度、非橢圓率三參數(shù)控制的各向異性Radon變換進行疊前地震數(shù)據(jù)重建的方法,考慮了偏移距較遠時地震數(shù)據(jù)的非雙曲線同相軸,通過稀疏約束控制Radon變換提高了地震數(shù)據(jù)的重建精度。針對時變Radon變換耗時較長的問題,在算子求解過程中采用快速迭代軟閾值算法(Fast iterative shrinkage-thresholding algorithm,FISTA)加快了收斂速度,提高了計算效率。模擬數(shù)據(jù)和實際數(shù)據(jù)的應用結(jié)果表明,各向異性Radon變換方法在提高數(shù)據(jù)重建準確度和計算效率方面效果良好。

各向異性Radon變換;遠偏移距;數(shù)據(jù)重建

地震數(shù)據(jù)受到施工條件、采集孔徑、儀器故障等影響,難免會出現(xiàn)缺失及不規(guī)則的情況,對后續(xù)處理效果產(chǎn)生很大影響,特別是表面多次波衰減、波動方程偏移等基于多道處理的算法[1-4]。KABIR等[5]利用拋物Radon變換[6]重建近偏移距缺失數(shù)據(jù),取得了良好效果,但只能處理規(guī)則采樣數(shù)據(jù),且需要多次迭代。SACCHI等[7]提出高分辨率Radon變換方法,并將其成功應用于近偏移距地震道重建,但在提高分辨率的同時增加了計算量。王維紅等[8]提出并應用加權(quán)拋物線Radon變換法進行數(shù)據(jù)規(guī)則化和重建,在近偏移距地震數(shù)據(jù)重建中取得了良好效果。目前應用Radon變換進行數(shù)據(jù)重建的研究主要集中在提高分辨率和減小計算量兩個方面,且多針對各向同性條件下的近偏移距地震數(shù)據(jù)。海上地震數(shù)據(jù)多為遠偏移距,且深層陡傾角成像需要準確的遠偏移距數(shù)據(jù),準確重建偏移距較遠的地震數(shù)據(jù)同樣十分重要[9-10]。

相比于拋物Radon變換,雙曲Radon變換分辨率更高,進行數(shù)據(jù)重建具有更高的精度。但目前使用雙曲Radon變換重建遠偏移距地震數(shù)據(jù)面臨以下問題:①常規(guī)雙曲Radon變換公式基于水平或中等傾角地質(zhì)構(gòu)造的反射同相軸在CMP道集上呈雙曲線走時的假設,而偏移距較遠的反射波同相軸已經(jīng)不能用簡單的雙曲線走時公式進行準確描述;②雙曲Radon變換具有時變屬性,變換算子只能在時間域表示,利用共軛梯度法求解計算量巨大,這也是時間域求解Radon變換的固有問題[11-13]。本文提出各向異性Radon變換數(shù)據(jù)重建方法,將常規(guī)雙曲變換公式改為由偏移距、慢度、非橢圓率三參數(shù)控制的各向異性變換公式[14];用FISTA算法[15]取代常規(guī)共軛梯度求解方法;采用稀疏約束控制各向異性Radon變換,通過正則化參數(shù)λ平衡反演誤差與模型稀疏化程度[16]。理論模型和實際資料重建結(jié)果表明,本文方法對中遠偏移距地震數(shù)據(jù)有著很好的重建效果。

1 方法原理

1.1 各向異性Radon變換原理

任意水平層狀介質(zhì)模型反射波走時方程可由級數(shù)展開得到與層厚、層速度及偏移距2階、4階、6階等有關(guān)的復雜函數(shù)。當模型為各向同性且偏移距較小時(最大偏移距小于反射界面深度),可忽略4階及更高階項,用簡單的雙曲線公式來描述縱波反射同相軸。常規(guī)雙曲Radon變換定義為:

(1)

式中:τ是截距時間;p是慢度(速度倒數(shù));x是炮檢距;d(t,x)為道集數(shù)據(jù)。

但長偏移距縱波反射同相軸時距曲線不滿足雙曲線規(guī)律,需用非雙曲公式進行精確描述[17]。TSVANKIN[18]給出了水平單層各向異性介質(zhì)長偏移距動校時差公式,ALKHALIFAH等[19]給出了該公式用動校速度和非橢圓率參數(shù)表示的形式。權(quán)衡計算精度與效率,本文采用ALKHALIFAH等[19]給出的非雙曲時差公式定義了各向異性Radon變換,其正變換公式如下。

(2)

與常規(guī)Radon變換相比,各向異性Radon變換在包含高階項的同時,增加了非橢圓率參數(shù)η,該參數(shù)與每個采樣時刻一一對應。要想獲得高精度的各向異性Radon變換域數(shù)據(jù),需要準確的η值,因此如何求取實際數(shù)據(jù)的η值顯得尤為重要。理論模型中的η值是已知的,對于實際地震數(shù)據(jù),本文采用雙參數(shù)速度分析的方法來求取準確的η值[20-21]。雖然各向異性Radon變換公式更加復雜,但是離散求解時矩陣的維數(shù)、大小未變,計算量與常規(guī)雙曲Radon變換相當。

為避免數(shù)據(jù)重建時產(chǎn)生假頻[22],得到高分辨率的采樣結(jié)果,需要合理選擇參數(shù),下面討論參數(shù)的選取準則。參數(shù)的采樣率為:

(3)

已知p=1/v,則(3)式可化為:

(4)

令Δv=vk-1-vk,得:

(5)

當vk和vk-1比較接近時,vk≈vk-1,公式(5)變?yōu)?

(6)

設原始數(shù)據(jù)中最大有效反射速度為vmax,要使(6)式對所有的有效波場均成立,則式中vk應替換為vmax。同時,在CSP或CMP道集上做變換時,Δv選取的范圍通常為50～100m/s,因而得到臨界值關(guān)系如下:

(7)

確定了Δp的值后,根據(jù)p=1/v給出p的初始值p0和最大值pmax,就可以在不增加數(shù)據(jù)重建計算量的同時保證重建的準確性[15]。

1.2 稀疏約束反演

為了得到更準確的重建結(jié)果,本文在常規(guī)最小二乘反演的基礎上,提出了基于稀疏約束的反演方法,采用正則化參數(shù)λ來平衡模型的稀疏化與反演誤差,從而控制重建精度。高分辨率Radon變換通常先定義反變換,再通過求解最小二乘約束下的目標函數(shù)求取正變換。Radon變換反變換公式如下。

d=Lm

(8)

式中:d和m分別表示離散原始數(shù)據(jù)與Radon域數(shù)據(jù)的矩陣形式;L和LT分別為Radon變換算子和伴隨變換算子。因Radon變換的非正交性,建立最小平方意義下的線性化反演誤差目標函數(shù),其范數(shù)條件的解為:

(9)

在時間域內(nèi)進行稀疏約束可以很好地提高Radon域數(shù)據(jù)的分辨率。將Radon域數(shù)據(jù)作為稀疏條件進行約束反演,對反演誤差取l2范數(shù),對稀疏約束模型取l1范數(shù),即通過線性反演問題的l1-l2混合范數(shù)求解,建立如下的目標函數(shù):

(10)

式中:λ是正則化參數(shù),用于平衡反演誤差與模型的稀疏化,其值越大,Radon域內(nèi)數(shù)據(jù)越稀疏。對于規(guī)則的模擬數(shù)據(jù),λ通常取0,而對實際數(shù)據(jù),合適的λ值可以使重建結(jié)果更加準確。

1.3 各項異性Radon變換求解算法

(11a)

rk+1=rk-αkA·pk

(11b)

(11c)

(11d)

pk+1=rk+1+βkpk

(11e)

(11f)

可見共軛梯度算法的遞推步驟較多,在實際運算中收斂速度較慢。

為了提高計算效率,我們采用快速迭代軟閾值算法(FISTA)加快反演的收斂速度,迭代更新公式如下:

m0=(LTL)z-1LTd

(12a)

x0=m0

(12b)

t0=1

(12c)

當k≥0時:

(13a)

(13b)

(13c)

式中:k是當前迭代次數(shù);soft是取軟閾值算子;a是Lipschitz常數(shù),a≥max[eig(LTL)]。這樣經(jīng)過數(shù)次迭代就可以獲得混合域高分辨率的Radon變換結(jié)果。

2 模型試算及實際資料處理

2.1 理論模型

為檢驗本文各向異性Radon變換數(shù)據(jù)重建方法的效果,建立了一個遠偏移距數(shù)據(jù)模型,模擬地震記錄共81道,道間距100m,每道500個采樣點,采樣間隔4ms,各向異性參數(shù)η取固定值0.1,如圖1a所示。圖1b為該模型不規(guī)則采樣及缺失遠偏移距數(shù)據(jù)的情況,即對原始數(shù)據(jù)前50道間隔抽樣,且使第66至75道(共10道)遠偏移距數(shù)據(jù)缺失。采用高分辨率雙曲Radon變換和各向異性Radon變換分別進行數(shù)據(jù)重建,用FISTA算法迭代100次,其Radon域結(jié)果分別如圖1c和圖1d所示,重建后的結(jié)果分別如圖1e和圖1f所示?？梢钥闯?各向異性變換使Radon域內(nèi)聚焦效果更好,重建結(jié)果更接近原始數(shù)據(jù),特別是遠偏移距缺失數(shù)據(jù)的重建結(jié)果。

為進一步比較兩種方法的重建誤差,圖2a和圖2b分別給出了高分辨率雙曲Radon變換和各向異性Radon變換重建數(shù)據(jù)與原始數(shù)據(jù)之差。圖2c,圖2d,圖2e,圖2f分別給出了原始數(shù)據(jù)、不規(guī)則缺失數(shù)據(jù)、高分辨率雙曲Radon變換重建數(shù)據(jù)和各向異性Radon變換重建數(shù)據(jù)的f-k譜,對比可見,不規(guī)則缺失數(shù)據(jù)出現(xiàn)嚴重假頻,兩種Radon變換重建方法都能明顯消除這些假頻。為準確描述重建精度,定義重建誤差如下:

(14)

高分辨率雙曲Radon變換重建誤差EΔ=9.84%,各向異性Radon變換重建誤差EΔ=3.71%,各向異性Radon變換重建誤差明顯更小。

在求解各項異性Radon變換時,我們分別采用共軛梯度算法和FISTA算法進行了對比測試,運行環(huán)境為Intel(R)Core(TM)i5-2300 CPU@2.80GHz,8GB DDR3 1333MHz內(nèi)存。FISTA算法在迭代100次的情況下重建誤差為3.71%,耗時159.8s;共軛梯度法至少迭代150次才能使誤差減小到4%以下,耗時為377.5s?？梢奆ISTA算法不但計算效率高,而且收斂速度快。

圖1 模型數(shù)據(jù)重建a 原始數(shù)據(jù); b 不規(guī)則缺失地震道數(shù)據(jù); c 高分辨率雙曲Radon變換結(jié)果; d 各向異性Radon變換結(jié)果; e 高分辨率雙曲Radon變換重建結(jié)果; f 各向異性Radon變換重建結(jié)果

圖2 重建誤差對比a 高分辨率雙曲Radon變換重建誤差; b 各向異性Radon變換重建誤差; c 原始數(shù)據(jù)f-k譜; d 不規(guī)則缺失數(shù)據(jù)f-k譜; e 高分辨率雙曲Radon變換重建數(shù)據(jù)f-k譜; f 各向異性Radon變換重建數(shù)據(jù)f-k譜

2.2 實際數(shù)據(jù)處理分析

采用我國某地區(qū)實際地震資料CMP道集,驗證本文各向異性Radon變換數(shù)據(jù)重建方法的實用性。該數(shù)據(jù)共126道,道間距24m,每道626個采樣點,采樣間隔4ms,記錄長度2.5s,其中有效反射波的速度范圍在1200～1600m/s。進行各向異性Radon變換數(shù)據(jù)重建之前,首先對道集進行雙參數(shù)速度掃描,得到每個地震同相軸所對應的η參數(shù),即掃描得到的η參數(shù)與同相軸個數(shù)相同;然后將η參數(shù)按時間方向插值,得到每個離散時間采樣點對應的η值,如圖3所示。

圖3 雙譜分析所得η參數(shù)曲線

圖4a是原始實際數(shù)據(jù),圖4b是對實際數(shù)據(jù)前50道間隔抽樣并缺失第101到110道遠偏移距數(shù)據(jù)的結(jié)果,圖4c為其各向異性Radon變換的結(jié)果,圖4d 為各向異性Radon變換方法數(shù)據(jù)重建結(jié)果。可以看出實際數(shù)據(jù)重建效果較好,說明本文方法具有較強的實用性。分別對原始實際數(shù)據(jù)及其不規(guī)則缺失地震道數(shù)據(jù)和重建后數(shù)據(jù)做f-k譜分析,結(jié)果如圖5a,圖5b,圖5c所示,可見不規(guī)則缺失地震道數(shù)據(jù)出現(xiàn)了明顯假頻,利用本文方法重建數(shù)據(jù)后有效地消除了這些假頻。

圖4 實際數(shù)據(jù)處理a 原始實際數(shù)據(jù); b 不規(guī)則缺失地震道數(shù)據(jù); c 實際數(shù)據(jù)各向異性Radon變換結(jié)果; d 實際數(shù)據(jù)各向異性Radon變換重建結(jié)果

為進一步提高各向異性Radon變換重建效果,在實際數(shù)據(jù)處理中加入一定程度的稀疏約束,并取實際數(shù)據(jù)與重建數(shù)據(jù)的第101至110道進行對比。λ的取值范圍一般為0.0001～0.0100,在該區(qū)間用二分法求解λ,僅需求解幾次即可得到最佳λ值0.002。圖6a為實際數(shù)據(jù)第101至110道,圖6b為常規(guī)Radon變換的重建結(jié)果,圖6c為λ=0時的各向異性Radon變換重建結(jié)果,圖6d為λ=0.002時的各向異性Radon變換重建結(jié)果。圖6b,圖6c,圖6d重建誤差分別為11.84%,6.52%,5.04%,說明各向異性Radon變換重建效果更好,且λ=0.002時重建結(jié)果更接近實際數(shù)據(jù)。

圖5 原始數(shù)據(jù)(a),不規(guī)則缺失數(shù)據(jù)(b)及重建數(shù)據(jù)(c)f-k譜

圖6 實際數(shù)據(jù)重建位置對比a 缺失位置的實際數(shù)據(jù); b常規(guī)Radon變換重建結(jié)果; c λ=0時各向異性Radon變換重建結(jié)果; d λ=0.002時各向異性Radon變換重建結(jié)果

3 結(jié)束語

本文提出的各向異性Radon變換重建方法可以更好地恢復遠偏移距地震數(shù)據(jù),彌補了常規(guī)Radon變換方法在遠偏移距數(shù)據(jù)重建方面的不足。對于不同地區(qū)的實際數(shù)據(jù),各向異性Radon變換前需進行雙參數(shù)速度掃描,獲得相應正確的各向異性參數(shù)η,確保Radon域內(nèi)數(shù)據(jù)的高分辨率。采用稀疏約束反演的方法可以減小數(shù)據(jù)重建誤差,得到更為準確的重建結(jié)果。在時間域Radon變換的求解過程中,采用FISTA算法取代常規(guī)的共軛梯度算法可以加快收斂速度,提高運算效率。

[1] 李慶忠.走向精確勘探的道路[M].北京:石油工業(yè)出版社,1994:98-112 LI Q Z.The way to obtain a batter resolution in seismic prospection[M].Beijing:Petroleum Industry Press,1994:98-112

[2] 黃新武,吳律,牛濱華.基于拋物線拉東變換的地震道重構(gòu)[J].中國礦業(yè)大學學報,2003,32(5):534-539 HUNG X W,WU L,NIU B H.Reconstruction of seismic trace by parabolic Radon transform[J].Journal of China University of Mining & Technology,2003,32(5):534-539

[3] 李遠欽.一種非線性Radon變換及非零偏移距VSP波場分離[J].石油物探,1994,33(3):33-39 LI Y Q.A non-linear Radon transform and non-zero offset VSP wave field separation[J].Geophysical Prospecting for Petroleum,1994,33(3):33-39

[4] 曾有良,樂友喜,單啟銅,等.基于高分辨率Radon變換的VSP波場分離方法[J].石油物探,2007,46(2):115-119 ZENG Y L,YUE Y X,SHAN Q T,et al.VSP wave filed separation based high resolution Radon transform[J].Geophysical Prospecting for Petroleum,2007,46(2):115-119

[5] KABIR M M N,VERSCHUUR D J.Restoration of missing offsets by parabolic Radon transform[J].Geophysical Prospecting,1995,43(3):347-368

[6] HAMPSON D.Inverse velocity stacking for multiple elimination[J].Journal of the Canadian Society of Exploration Geophysicists,1986,22(1):44-55

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

[8] 王維紅,裴江云,張劍鋒.加權(quán)拋物Radon變換疊前地震數(shù)據(jù)重建[J].地球物理學報,2007,50(3):851-859 WANG W H,PEI J Y,ZHANG J F.Prestack seismic data reconstruction using weighted parabolic Radon transform[J].Chinese Journal of Geophysics,2007,50(3):851-859

[9] 張紅梅,劉洪.基于稀疏離散τ-p變換的非均勻地震道重建[J].石油物探,2006,45(2):141-145 ZHANG H M,LIU H.Sparseness discreteτ-ptransform in irregular seismic trace reconstruction[J].Geophysical Prospecting for Petroleum,2006,45(2):141-145

[10] 李晶晶,孫成禹,謝俊法.相對保幅的拋物線Radon變換法地震道重建[J].石油物探,2014,53(2):181-187 LI J J,SUN C Y,XIE J F.Seismic trace reconstruction by relative amplitude preserved parabolic Radon transform[J].Geophysical Prospecting for Petroleum,2014,53(2):181-187

[11] 孫成禹,尚新民,石翠翠,等.影響地震數(shù)據(jù)相位特征的因素分析[J].石油物探,2011,50(5):444-454 SUN C Y,SHANG X M,SHI C C,et al.Analysis of influence factors on phase characteristics of seismic data[J].Geophysical Prospecting for Petroleum,2011,50(5):444-454

[12] 劉法啟,張關(guān)泉.借助廣義 Radon變換進行方位角校正[J].石油物探,1996,35(4):43-51 LIU F Q,ZHANG G Q.Azimuth correction with the aid of generalized Radon transform[J].Geophysical Prospecting for Petroleum,1996,35(4):43-51

[13] 張軍華,呂寧,雷凌,等.拋物線拉冬變換消除多次波的應用要素分析[J].石油地球物理勘探,2004,39(4):398-405 ZHANG J H,LV N,LEI L,et al.Analysis of application factors in multiple attenuation by parabolic Radon transform[J].Oil Geophysical Prospecting,2004,39(4):398-405

[14] 鞏向博,韓立國,李洪建.各向異性Radon 變換及其在多次波壓制中的應用[J].地球物理學報,2014,57(9):2928-2936 GONG X B,HAN L G,LI H J.Anisotropic Radon transform and its application to demultiple[J].Chinese Journal of Geophysics,2014,57(9):2928-2936

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

[16] BECK A,TEBOULLE M.A fast iterative shrinkage-thresholding algorithm for linear inverse problems[J].SIAM Journal on Imaging Sciences,2009,2(1):183-202

[17] 尤建軍,常旭,劉伊克.VTI介質(zhì)長偏移距非雙曲動校正公式優(yōu)化[J].地球物理學報,2006,49(6):1770-1778YOU J J,CHANG X,LIU Y K.Optimization of nonhyperbolic moveout correction equation of long offset seismic data in VTI media[J].Chinese Journal of Geophysics,2006,49(6):1770-1778

[18] TSVANKIN T.Nonhyperbolic reflection moveout in anisotropic media[J].Geophysics,1994,59(8):1290-1304

[19] ALKHALIFAH T,TSVANKIN I.Velocity analysis for transversely isotropic media[J].Geophysics,1995,60(5):1550-1566

[20] ABBAD B,URSIN B,RAPPIN D.Automatic nonhyperbolic velocity analysis[J].Geophysics,2009,74(2):U1-U12

[21] 張博,韓立國,譚塵青.VTI介質(zhì)平移初至走時動校正雙參數(shù)反演[J].石油物探,2012,51(2):119-124 ZHANG B,HAN L G,TAN C Q.Shifted first-arrival travel time NMO velocity andηinversion for VTI media[J].Geophysical Prospecting for Petroleum,2012,51(2):119-124

[22] 王立歆,李強,姬小兵,等.用Radon變換法消除沙丘鳴震的應用及效果分析[J].石油物探,2002,41(1):89-91 WANG L X,LI Q,JI X B,et al.Application of Radon transform to elimination of dune ringing and its effect analysis[J].Geophysical Prospecting for Petroleum,2002,41(1):89-91

[23] 張保衛(wèi).Radon變換及其在地震數(shù)據(jù)處理中的應用[D].西安:長安大學,2007 ZHANG B W.Radon transform and its application in seismic data processing[D].Xi’an:Chang’an University,2007

[24] 熊登,趙偉,張劍鋒.混合域高分辨率拋物Radon變換及在衰減多次波中的應用[J].地球物理學報,2009,52(4):1068-1077 XIONG D,ZHAO W,ZHANG J F.Hybrid-domain high-resolution parabolic Radon transform and its application to demultiple.[J].Chinese Journal of Geophysics,2009,52(4):1068-1077

(編輯:戴春秋)

Prestack seismic data reconstruction by anisotropic Radon transform

WANG Shengchao,HAN Liguo,GONG Xiangbo

(CollegeofGeo-ExplorationScienceandTechnology,JilinUniversity,Changchun130012,China)

Although the conventional hyperbolic Radon transform owns high resolution to reconstruct seismic data,it is difficult to reconstruct the far-offset seismic data,because the event of far-offset seismic signals is no longer satisfied with the hyperbolic assumption.Therefore,we proposed a prestack seismic data reconstruction method by anisotropic Radon transform to consider the effect of far-offset factor.The formula of anisotropic Radon transform is controlled by offset,slowness and anellipticity.The travel-time curve of horizontal layer can be accurately described for far-offset by the anisotropic anellipticity parameter of the non-hyperbolic moveout formula.We also use sparsity constraint to improve the accuracy of reconstruction.During the procedure of operator solving,fast iterative shrinkage-thresholding algorithm (FISTA) is used to speed up the convergence and reduce the calculation time.The results of theoretical model and field data demonstrate that anisotropic Radon transform has good effect on the reconstruction accuracy and calculation efficiency.

anisotropic Radon transform,far-offset,seismic data reconstruction

2016-01-14;改回日期：2016-03-11。

王升超(1990—),男,碩士在讀,從事地震勘探處理技術(shù)研究。

國家自然科學基金項目(41204078)資助。

This research is financially supported by National Natural Science Foundation of China (Grant No.41204078).

P631

A

1000-1441(2016)06-0808-08

10.3969/j.issn.1000-1441.2016.06.005