金昌昆,蔣亞洲,張建中

(1.海底科學(xué)與探測(cè)技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室,山東青島266100;2.中國(guó)海洋大學(xué)海洋地球科學(xué)學(xué)院,山東青島266100)

基于曲波變換的地震反射波同相軸斜率拾取方法

金昌昆1,2,蔣亞洲1,2,張建中1,2

(1.海底科學(xué)與探測(cè)技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室,山東青島266100;2.中國(guó)海洋大學(xué)海洋地球科學(xué)學(xué)院,山東青島266100)

局部相關(guān)反射波同相軸的斜率反映了射線參數(shù)、慢度水平分量等地震波運(yùn)動(dòng)學(xué)信息,被用于地下速度建模和偏移成像等地震數(shù)據(jù)處理中。拾取反射波同相軸斜率是應(yīng)用斜率數(shù)據(jù)的必然環(huán)節(jié)。提出了一種基于曲波變換(Curvelet transform)的地震反射波同相軸斜率拾取方法。該方法根據(jù)曲波變換所具有的各向異性和局部方向性的性質(zhì),利用地震信號(hào)曲波變換后的曲波系數(shù),確定地震反射同相軸位置和延伸方向,再將延伸方向換算成斜率值。理論數(shù)據(jù)測(cè)試和實(shí)際地震資料試應(yīng)用結(jié)果表明,與常用的基于傾斜疊加以及直線—雙曲線疊加的反射波同相軸斜率拾取方法相比,該方法具有更高的精度。

曲波變換;反射波同相軸;斜率拾取;傾斜疊加;直線—雙曲線疊加

地震反射波同相軸斜率是指在地震記錄道集(如共炮點(diǎn)道集、共接收點(diǎn)道集)上,來(lái)自同一反射界面的局部相關(guān)同相軸在某道上的切線的斜率,也就是該道上的反射走時(shí)沿測(cè)線方向的梯度,它反映了射線參數(shù)和慢度水平分量。同相軸斜率包含了反射波幾何信息,因此許多學(xué)者研究了利用斜率數(shù)據(jù)進(jìn)行地震資料的處理和成像。Hua等[1-2]利用炮點(diǎn)的斜率獲取射線角,提出了計(jì)算效率很高的二維和三維偏移方法;Hertwech等[3]把斜率參數(shù)用于共反射面元疊加中;Fomel[4]和Cooke等[5]應(yīng)用斜率進(jìn)行獨(dú)立于速度模型的時(shí)間偏移成像;杜婧等[6]研究了基于斜率的VSP波場(chǎng)分離;Billette等[7]把炮點(diǎn)、接收點(diǎn)的斜率和雙程走時(shí)作為觀測(cè)數(shù)據(jù),提出基于斜率的立體層析成像方法;倪瑤等[8]對(duì)立體層析成像方法進(jìn)行了應(yīng)用測(cè)試;Guillaume等[9]研究了共偏移距矢量道集的斜率層析成像方法。

在這些方法及其應(yīng)用中,斜率的拾取是不可或缺的環(huán)節(jié),而且斜率數(shù)據(jù)的精度必然影響著它們的應(yīng)用效果。目前,斜率拾取主要采用傾斜疊加技術(shù)(slant stack)[10-11]。然而,在同相軸的線性程度不高時(shí),利用傾斜疊加技術(shù)拾取斜率的效果往往不佳;特別是在近偏移距,局部同相軸變化的線性程度往往很低,用傾斜疊加方法拾取斜率的誤差較大。為此,黃忠來(lái)等[12]提出了一種基于直線—雙曲線疊加的斜率拾取方法。在共偏移距域用線性疊加方法拾取中心點(diǎn)上的斜率值,在共中心點(diǎn)域用雙曲線疊加方法拾取半偏移距位置上的斜率值,然后把中心點(diǎn)和半偏移距上的斜率值換算成炮點(diǎn)和接收點(diǎn)上的斜率值。模型數(shù)據(jù)和實(shí)際資料的測(cè)試結(jié)果證明此方法比傾斜疊加方法具有更高的精度。

曲波變換是Candès等[13]于1999年提出的一種多尺度變換方法,它也是一種多分辨、帶通、具有方向性的函數(shù)分析方法。2002年,Candès等[14]又提出了第二代曲波變換的框架系,2005年提出了兩種快速離散曲波變換方法[15],即非均勻快速傅里葉變換(Unequally-spaced FFT,USFFT)和基于Wrapping算法。相對(duì)于第一代曲波變換,第二代曲波變換具有概念更簡(jiǎn)潔、運(yùn)算速度更快、冗余性更小等優(yōu)點(diǎn)。曲波變換具有的局部性、多尺度性和多方向性,可以用來(lái)較好地表征地震數(shù)據(jù)的局部特征[15],從而越來(lái)越多地被應(yīng)用于地震資料處理。Herrmann等[16]提出了全新的基于曲波變換的地震數(shù)據(jù)處理的概念;Neelamani等[17]利用曲波變換壓制地震相干噪聲和隨機(jī)噪聲;張恒磊等[18]提出在曲波域采用非線性閾值方法對(duì)疊前地震資料隨機(jī)噪聲進(jìn)行壓制;Mansour等[19]利用曲波變換實(shí)現(xiàn)了海上隨機(jī)采樣及數(shù)據(jù)重構(gòu);劉國(guó)昌等[20]利用曲波變換進(jìn)行地震數(shù)據(jù)的缺失道插值;Hubb等[21]把曲波變換應(yīng)用于地震成像中;張廣智等[22]利用曲波變換的多尺度性識(shí)別裂縫發(fā)育帶。

地震反射同相軸斜率拾取需要自動(dòng)追蹤來(lái)自同一反射界面的局部相關(guān)同相軸,并確定該同相軸的延伸方向。根據(jù)曲波變換具有的局部性和方向性的優(yōu)勢(shì),本文提出一種基于第二代曲波變換的地震同相軸斜率拾取方法,通過(guò)理論合成數(shù)據(jù)和實(shí)際地震資料的測(cè)試應(yīng)用,并與基于傾斜疊加的方法和基于直線—雙曲線疊加的方法對(duì)比,驗(yàn)證本文方法的有效性。

1 方法原理

函數(shù)f(t1,t2)∈L2的離散曲波系數(shù)[15]表示為:

(1)

圖1 離散曲波空間頻率域區(qū)域分塊示意

(2)

把每一層尺度ji上的曲波系數(shù)cjilk投影到尺度jp的網(wǎng)格gp中,尺度ji上網(wǎng)格點(diǎn)坐標(biāo)投影到尺度jp的網(wǎng)格gp的子集Ai,p(k)中[23],有:

(3)

同理,把尺度ji上的方向參數(shù)l投影到尺度jp中的方向參數(shù)Lp的子集Di,p(l)中,有:

(4)

式中:Li是尺度ji上總的方向個(gè)數(shù);Lp是尺度jp上總的方向個(gè)數(shù)。

在選取的jp層中的每個(gè)方向l和坐標(biāo)k=(k1,k2)上,計(jì)算曲波系數(shù)的模值和大小:

(5)

(5)式中計(jì)算Mlk時(shí),方向參數(shù)l只需從1變化到Lp/2即可,這是因?yàn)閘和l+Lp/2相差180°,實(shí)際表示同一個(gè)方向。一般情況下求出Mlk后,就可以得到在每一個(gè)網(wǎng)格點(diǎn)k=(k1,k2)上的主方向l0(k)以及主方向?qū)?yīng)的最大曲波系數(shù)值Mlk,即:

(6)

根據(jù)Mlk來(lái)判斷同相軸振幅極大值位置,如果Mlk0是每一道地震記錄曲波系數(shù)的模值Mlk的局部極大值,且系數(shù)Mlk0大于一定的閾值,坐標(biāo)k0便是同相軸振幅極大值位置,此時(shí)此點(diǎn)同相軸斜率值為:

(7)

式中:θl0是方向參數(shù)l0對(duì)應(yīng)的角度值;Δt是地震記錄時(shí)間采樣率;Δx是地震記錄的道間距。

算法的實(shí)現(xiàn)流程如下。

1) 讀入地震記錄,選出將要求取的同相軸斜率的地震記錄。

2) 進(jìn)行曲波變換。采用Wrapping算法對(duì)數(shù)據(jù)進(jìn)行變換,得到曲波系數(shù)。

3) 提取位置和方向參數(shù)。選取合適的曲波系數(shù),求取同相軸坐標(biāo)位置和對(duì)應(yīng)位置方向參數(shù)。

4) 由(7)式計(jì)算該點(diǎn)上的斜率值。

2 模型測(cè)試與實(shí)際資料試應(yīng)用

2.1 兩層水平均勻介質(zhì)模型理論數(shù)據(jù)

首先選用能夠直接求得炮點(diǎn)和接收點(diǎn)的理論斜率值的模型數(shù)據(jù)進(jìn)行測(cè)試分析,以便對(duì)用本文方法自動(dòng)拾取的同相軸斜率值與理論斜率值進(jìn)行比較。圖2所示是一個(gè)簡(jiǎn)單的兩層水平均勻介質(zhì)模型。炮點(diǎn)S坐標(biāo)為xs,接收點(diǎn)R的坐標(biāo)為xr,上層介質(zhì)速度為v1,下層介質(zhì)速度為v2,則一次反射波時(shí)距方程可表示為[24]:

(8)

這樣,接收點(diǎn)pr的斜率值為:

(9)

圖2所示模型中炮點(diǎn)坐標(biāo)S為坐標(biāo)原點(diǎn),接收點(diǎn)與炮點(diǎn)距離為x。界面深度h=200m,取v1=1000m/s,v2=2000m/s,根據(jù)公式(9)計(jì)算接收點(diǎn)上的理論斜率值pr。

圖2 兩層水平均勻介質(zhì)模型

同時(shí)基于褶積模型合成該模型的單炮地震記錄數(shù)據(jù)。其中,地震子波用Ricker子波,中心頻率30Hz,均勻布設(shè)256個(gè)接收點(diǎn),接收點(diǎn)間距10m,中心點(diǎn)放炮,每道采樣點(diǎn)數(shù)為1024,采樣周期為2ms。合成的共炮點(diǎn)地震記錄如圖3所示。對(duì)該合成記錄用傾斜疊加方法(疊加窗空間長(zhǎng)度為9道)直接拾取接收點(diǎn)的斜率值pr,如圖4中的虛線所示。再應(yīng)用本文方法確定各道同相軸的位置及斜率值pr,其中,取Lp=256,即在尺度jp中的方向數(shù)為256個(gè),拾取結(jié)果如圖4中的點(diǎn)線所示。圖5 是這兩種方法拾取的斜率值相對(duì)計(jì)算理論值的誤差曲線。從圖中可以看出,由于近偏移距處的斜率值較小以及時(shí)距曲線曲率較大,兩種方法在近偏移距處的誤差都明顯大于遠(yuǎn)偏移距處的誤差;但即使在近偏移距處,本文方法拾取斜率值的精度也明顯高于傾斜疊加方法。

圖3 兩層水平介質(zhì)模型的合成地震記錄

圖4 計(jì)算與拾取的兩層模型各接收點(diǎn)的斜率值

圖5 圖4中兩種拾取值相對(duì)于理論值的誤差

一般地,曲波變換的方向數(shù)越多,斜率拾取的精度越高。圖6所示為方向數(shù)分別取64和512時(shí)拾取結(jié)果在偏移距250～1000m的相對(duì)誤差?？梢娖凭嘈∮?00m時(shí),兩種結(jié)果的相對(duì)誤差都比較大;在300～550m,方向數(shù)為64時(shí)拾取斜率誤差(虛線)要大一些;方向數(shù)為512時(shí)拾取結(jié)果的誤差(實(shí)線)在偏移距400m左右基本穩(wěn)定,而方向數(shù)為64的結(jié)果誤差在偏移距大于500m時(shí)才趨于穩(wěn)定;當(dāng)偏移距大于550m時(shí),走時(shí)斜率較大,兩種方法拾取結(jié)果的誤差均很小。由此可以看出,當(dāng)方向數(shù)較少時(shí),所得結(jié)果的誤差較大。

圖6 取不同方向數(shù)時(shí)拾取結(jié)果的相對(duì)誤差

2.2 Marmousi模型聲波方程模擬數(shù)據(jù)

對(duì)Marmousi模型用聲波方程模擬(道間距25m,采樣周期4ms)得到的地震數(shù)據(jù)[7],首先計(jì)算其理論斜率值。取出某個(gè)反射同相軸的走時(shí),在一個(gè)包含9道的窗口內(nèi)用二次曲線擬合反射走時(shí),對(duì)擬合的二次曲線求導(dǎo)獲取共炮域各接收點(diǎn)的走時(shí)沿測(cè)線方向的導(dǎo)數(shù),這樣求得的導(dǎo)數(shù)值即可看作理論斜率值,作為其它方法拾取結(jié)果比較的參考基準(zhǔn)。

然后采用3種不同方法拾取模擬數(shù)據(jù)炮點(diǎn)和接收點(diǎn)的斜率值。第1種是用常規(guī)的傾斜疊加方法(疊加窗空間長(zhǎng)度為9道,時(shí)間窗口為40采樣點(diǎn))在共炮域拾取接收點(diǎn)的斜率;第2種是文獻(xiàn)[12]提出的基于直線—雙曲線疊加的斜率拾取方法;第3種是用本文的曲波變換方法(方向參數(shù)Lp=256)進(jìn)行拾取。這里給出第32炮記錄的一個(gè)反射同相軸的拾取結(jié)果。

圖7箭頭指示拾取的反射同相軸;圖8是分別采用傾斜疊加方法、直線—雙曲線疊加方法和本文方法拾取的斜率值與理論斜率值的對(duì)比;圖9是3種方法分別拾取的斜率值相對(duì)理論斜率值的相對(duì)誤差?？梢钥闯?本文方法可應(yīng)用于復(fù)雜模型數(shù)據(jù),并且比傾斜疊加方法及直線—雙曲線疊加方法具有更高的拾取精度。

圖8 圖7箭頭所示同相軸對(duì)應(yīng)的斜率值

2.3 實(shí)際地震資料試應(yīng)用

對(duì)經(jīng)過(guò)去噪和靜校正等常規(guī)處理后的某復(fù)雜地形條件下的實(shí)際地震資料(道間距15m,采樣周期4ms),采用上述3種方法拾取同相軸斜率。其中,通過(guò)對(duì)同相軸進(jìn)行二次曲線(9點(diǎn))擬合再求導(dǎo)數(shù),并把求得的導(dǎo)數(shù)值作為與其它方法結(jié)果進(jìn)行比較時(shí)的理論斜率值;常規(guī)的傾斜疊加方法取疊加窗道數(shù)為9道;本文方法拾取時(shí)的方向參數(shù)Lp取256。這里給出第82炮記錄中的兩個(gè)反射同相軸(見圖10中箭頭所示)的斜率拾取結(jié)果。圖10 是實(shí)際地震資料第82炮的記錄;圖11是3種方法拾取結(jié)果與理論值的對(duì)比;圖12是3種方法拾取結(jié)果的相對(duì)誤差。由圖10，圖11和圖12可見,由于實(shí)際地震資料中含有噪聲干擾,傾斜疊加方法的拾取結(jié)果(圖11中綠線)明顯受到了噪聲的影響,斜率值曲線具有相對(duì)劇烈的跳躍變化;直線—雙曲線疊加斜率拾取方法比較穩(wěn)定,相對(duì)誤差平均為2.1%(圖12中藍(lán)實(shí)線);本文方法拾取的斜率值曲線(圖11中藍(lán)線)無(wú)明顯的跳躍現(xiàn)象,與理論值更符合,平均誤差為1.8%(圖12中黑實(shí)線)。

實(shí)際資料試應(yīng)用的結(jié)果表明,本文基于曲波變換的斜率拾取方法比基于傾斜疊加及直線—雙曲線疊加的斜率拾取方法具有更高的拾取精度以及較強(qiáng)的抗噪能力。

圖10 實(shí)際地震資料第82炮的地震記錄(箭頭標(biāo)注的是3種方法拾取的同相軸)

圖11 圖10箭頭所示同相軸對(duì)應(yīng)的斜率值

圖12 圖11中3種方法拾取值相對(duì)于理論值的誤差

3 結(jié)論與認(rèn)識(shí)

根據(jù)曲波變換的各向異性和局部方向性的性質(zhì),對(duì)地震記錄進(jìn)行曲波變換,利用所得到的曲波系數(shù),可以準(zhǔn)確確定反射同相軸位置和延伸方向等參數(shù),從而提出了一種基于曲波變換的地震反射同相軸斜率拾取方法。與常用的傾斜疊加方法相比較,傾斜疊加方法拾取同相軸斜率需要選擇合適的窗口長(zhǎng)度,窗口長(zhǎng)度對(duì)拾取斜率的精度影響很大,而本文方法拾取同相軸斜率是對(duì)全部地震數(shù)據(jù)進(jìn)行曲波變換,在空間方向上無(wú)需選擇窗口。理論合成數(shù)據(jù)測(cè)試和實(shí)際地震資料試應(yīng)用的結(jié)果表明,本文提出的基于曲波變換的同相軸斜率拾取方法比常用的傾斜疊加方法和直線—雙曲線疊加斜率拾取方法具有更高的精度和更強(qiáng)的抗噪能力。此外,本文方法拾取斜率的精度與曲波變換所選取方向數(shù)的多少有關(guān),方向數(shù)越多,對(duì)同相軸方向的檢測(cè)精度越高,從而拾取斜率的精度也越高。

[1] Hua B,McMechan G A.Parsimonious 2D prestack Kirchhoff depth migration[J].Geophysics,2003,68(3):1043-1051

[2] Hua B,McMechan G A.Parsimonious 3D post-stack Kirchhoff depth migration[J].Geophysical Prospecting,2005,53(4):507-522

[3] Hertwech T,Schileicher J,Mann J.Data stacking beyond CMP[J].The Leading Edge,2007,26(7):818-827

[4] Fomel S.Velocity-independent time-domain seismic imaging using local event slopes[J].Geophysics,2007,72(3):S139-S147

[5] Cooke D,Bóna A,Hansen B.Simultaneous time imaging,velocity estimation,and multiple suppression using local event slopes[J].Geophysics,2009,74(6):WCA65-WCA73

[6] 杜婧,王尚旭,劉國(guó)昌,等.基于局部斜率屬性的VSP波場(chǎng)分離研究[J].地球物理學(xué)報(bào),2009,52(7):1867-1872 Du J,Wang S X,Liu G C,et al.VSP wavefield separation using local slope attributes[J].Chinese Journal of Geophysics,2009,52(7):1867-1872

[7] Billette F,Lambaré G.Velocity macro-model estimation from seismic reflection data by stereotomography[J].Geophysical Journal International,1998,135(2):671-690

[8] 倪瑤,楊鍇,陳寶書.立體層析反演方法理論分析與應(yīng)用測(cè)試[J].石油物探,2013,52(2):121-130 Ni Y,Yang K,Chen B S.Stereotomography inversion method:theory and application testing[J].Geophysical Prospecting for Petroleum,2013,52(2):121-130

[9] Guillaume P,Montel J P,McCarthy A,et al.Non-linear slope tomography from common offset vector volumes as applied to a high density land WAZ survey from the Sultanate of Oman[J].Expanded Abstracts of 80thAnnual Internat SEG Mtg,2010,4395-4399

[10] Lambaré G,Alerini M,Baina R,et al.Stereotomography:a semi-automatic approach for velocity macromodel estimation[J].Geophysical Prospecting,2004,52(6):671-681

[11] Lambaré G,Alerini M,Podvin P.Stereotomographic picking in practice[J].Expanded Abstracts of 74thAnnual Internat SEG Mtg,2004,2343-2346

[12] 黃忠來(lái),張建中,蔣亞洲.基于直線和雙曲線疊加的斜率自動(dòng)拾取方法[J].地球物理學(xué)進(jìn)展,2013,28(6):3007-3014 Huang Z L,Zhang J Z,Jiang Y Z.Automatic method for picking slopes of locally coherent events based on linear and hyperbolic stack[J].Progress in Geophysics,2013,28(6):3007-3014

[13] Candès E J,Donoho D L.Curvelets:a suprisingly effective nonadaptive representation for objects with edges,in curves and surface fitting:saint-malo 1999[M].Nashville:Vanderbilt University Press,2000:105-120

[14] Candès E J,Donoho D L.New tight frames of curvelets and optimal representations of objects with smooth singularities[D].USA:Department of Statistics,Stanford University,2002[15] Candès E J,Demanet L,Donoho D L,et al.Fast discrete curvelet transform[J].Multiscale Modeling and Simulation,2005,5(3):861-899

[16] Herrmann F J,Wang D L,Hennerfent G.Curvelet-based seismic data processing:a multiscale and nonlinear approach[J].Geophysics,2008,73(1):A1-A5

[17] Neelamani R,Baumstein A I,Gillard D G.Coherent and random noise attenuation using the curvelet transform[J].The Leading Edge,2008,27(2):240-248

[18] 張恒磊,劉天佑,張?jiān)拼?基于高階相關(guān)的Curvelet域和空間域的傾角掃描噪聲壓制方法[J].石油地球物理勘探,2010,45(2):208-214 Zhang H L,Liu T Y,Zhang Y C.High order correlation based dip angle scanning noise elimination method in Curvelet domain and space domain[J].Oil Geophysical Prospecting,2010,45(2):208-214

[19] Mansour H,Wason H,Herrmann F J.Randomized marine acquisition with compressive sampling matrices[J].Geophysical Prospecting,2012,60(4):648-662

[20] 劉國(guó)昌,陳小宏,郭志峰.基于Curvelet變換的缺失地震數(shù)據(jù)插值方法[J].石油地球物理勘探,2011,46(2):237-246 Liu G C,Chen X H,Guo Z F.Missing seismic data rebuilding by interpolation based on Curvelet transform[J].Oil Geophysical Prospecting,2011,46(2):237-246

[21] Huub D,Hoop M V.Leading-order seismic imaging using curvelets[J].Geophysics,2007,72(6):S231-S248

[22] 張廣智,鄭靜靜,印興耀.基于Curvelet變換的多尺度性識(shí)別裂縫發(fā)育帶[J].石油地球物理勘探,2011,46(5):757-762 Zhang G Z,Zheng J J,Yin X Y.Identification technology of fracture zone and its strike based on the Curvelet transform[J].Oil Geophysical Prospecting,2011,46(5):757-762

[23] Geb?ck T,Koumoutsakos P.Edge detection in microscopy images using curvelets[J].BMC Bioinformatics,2009,10(1):1-14

[24] 陸基孟.地震勘探原理[M].北京:石油工業(yè)出版社,1982:22-23 Lu J M.The principle of seismic prospecting[M].Beijing:Petroleum Industry Press,1982:22-23

(編輯：朱文杰)

An automatic picking method for slopes of locally coherent reflection events based on Curvelet transform

Jin Changkun1,2,Jiang Yazhou1,2,Zhang Jianzhong1,2

(1.KeyLaboratoryofSubmarineGeociencesandProspectingTechniques,MinistryofEducationofChina,Qingdao266100,China;2.CollegeofMarineGeosciences,OceanUniversityofChina,Qingdao266100,China)

The slopes of locally coherent reflection events represent ray parameters or slowness horizontal components and other kinematic information,which are used in seismic data processing,such as velocity model building,migration imaging,and etc.Picking slopes of reflection events is an essential step in seismic data processing.An automatic picking method is proposed for slopes of reflection events based on Curvelet transform.Based on the characteristics of local directionality and anisotropy for Curvelet transform,the method use the Curvelet coefficients from seismic data transformed by fast discrete Curvelet transform to obtain the location and direction of locally coherent reflection events.And the slopes are finally estimated from the above directions.The tests with both synthetic and field data show that this method has a higher precision compared with either the existing slant stack method or linear-hyperbolic stack method.

Curvelet transform,reflection events,slope picking,slant stack,linear-hyperbolic stack

2014-04-17;改回日期：2014-09-15。

金昌昆(1987—),男,博士在讀,研究方向?yàn)榈卣鹳Y料處理和層析反演。

張建中(1963—),男,教授,主要從事地球物理勘探方法研究工作。

國(guó)家自然科學(xué)基金項(xiàng)目(41074077、41230318)資助。

P631

A

1000-1441(2015)04-0414-06

10.3969/j.issn.1000-1441.2015.04.007