金昌昆,蔣亞洲,張建中

(1.海底科學與探測技術教育部重點實驗室,山東青島266100;2.中國海洋大學海洋地球科學學院,山東青島266100)

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

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

(1.海底科學與探測技術教育部重點實驗室,山東青島266100;2.中國海洋大學海洋地球科學學院,山東青島266100)

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

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

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

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

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

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

1 方法原理

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

(1)

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

(2)

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

(3)

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

(4)

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

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

(5)

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

(6)

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

(7)

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

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

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

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

3) 提取位置和方向參數(shù)。選取合適的曲波系數(shù),求取同相軸坐標位置和對應位置方向參數(shù)。

4) 由(7)式計算該點上的斜率值。

2 模型測試與實際資料試應用

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

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

(8)

這樣,接收點pr的斜率值為:

(9)

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

圖2 兩層水平均勻介質模型

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

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

圖4 計算與拾取的兩層模型各接收點的斜率值

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

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

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

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

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

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

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

圖8 圖7箭頭所示同相軸對應的斜率值

2.3 實際地震資料試應用

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

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

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

圖11 圖10箭頭所示同相軸對應的斜率值

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

3 結論與認識

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

[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] 杜婧,王尚旭,劉國昌,等.基于局部斜率屬性的VSP波場分離研究[J].地球物理學報,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] 倪瑤,楊鍇,陳寶書.立體層析反演方法理論分析與應用測試[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] 黃忠來,張建中,蔣亞洲.基于直線和雙曲線疊加的斜率自動拾取方法[J].地球物理學進展,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] 張恒磊,劉天佑,張云翠.基于高階相關的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] 劉國昌,陳小宏,郭志峰.基于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變換的多尺度性識別裂縫發(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—),男,博士在讀,研究方向為地震資料處理和層析反演。

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

國家自然科學基金項目(41074077、41230318)資助。

P631

A

1000-1441(2015)04-0414-06

10.3969/j.issn.1000-1441.2015.04.007