王清振,姜秀娣,翁 斌,丁繼才,陳劍軍

(中海油研究總院,北京100028)

高抗噪性三維體曲率分析技術(shù)及其在高陡地層發(fā)育區(qū)的應(yīng)用

王清振,姜秀娣,翁 斌,丁繼才,陳劍軍

(中海油研究總院,北京100028)

曲率屬性可以較好地描述巖層的彎曲、褶皺、斷層和裂縫等信息,但是由于曲率屬性提取過程中需要計算地震數(shù)據(jù)的二階微分,因此對噪聲比較敏感,影響了結(jié)果的穩(wěn)定性和可信度。在模型試算的基礎(chǔ)上,提出了一種新的高抗噪性的體曲率計算方法。首先利用小波變換方法提取三維地震資料的相位余弦體;然后基于復地震道理論計算三維體視傾角屬性,對視傾角進行中值濾波以壓制微分計算中生成的噪聲;最后利用兩個方向的視傾角屬性計算得到各類高信噪比的體曲率屬性。模型實驗和實際資料處理結(jié)果表明,該方法計算的曲率屬性對噪聲不敏感,適用于低信噪比資料的處理,特別是在陡傾角地層發(fā)育的地區(qū),計算結(jié)果比C3等相干屬性假象更少,能夠真實地反映地下斷層及巖性邊界等地質(zhì)信息。

斷層檢測;曲率;相位余弦;相干;小波變換

曲率屬性分析是近年來興起的一種利用地層的彎曲程度進行不連續(xù)性檢測的方法,與相干算法相比,該方法對復雜河道、裂縫,小斷層及褶皺的刻畫能力更為優(yōu)越[1-2]。LISLE[3]對高斯曲率與裂縫的相關(guān)性進行了研究,驗證了曲率和裂縫之間的對應(yīng)關(guān)系;ROBERTS[4]對地震曲率屬性進行了系統(tǒng)的分類,并闡述了多種曲率屬性如平均曲率、高斯曲率、極大曲率、極小曲率、最大正曲率和最小負曲率的定義以及各自的特點,給出了利用解釋的地質(zhì)層位計算二維曲率的方法。但由于人工解釋的地質(zhì)層位存在偏差,計算的曲率準確性受到很大影響,容易引起構(gòu)造假象。為此AL-DOSSARY等[5]提出了三維體曲率的計算方法,它首先利用原始地震資料的振幅信息計算地層傾角,然后通過相似性掃描的方法來計算曲面的一階導數(shù),將曲率算法從二維推廣到三維。

國內(nèi)學者近年來也對曲率屬性進行了大量研究,CHEN等[6]利用小波變換對曲率屬性進行了多尺度分解,并將得到的多個尺度的曲率屬性進行融合顯示,提高了大斷層和小斷層的綜合顯示能力;GAO[7]提出了一種用于小斷層識別的梯度曲率屬性,取得了較好效果;王世星[8]和印興耀等[9]給出了一種基于傾角掃描的體曲率提取方法;QI等[10]給出了斷層控制的喀斯特地區(qū)的曲率屬性特征。楊國權(quán)等[11]、李福強等[12]分別用5×5,7×7網(wǎng)格來擬合二次曲面以期得到更高精度的曲率屬性,并將其權(quán)系數(shù)表示成微分算子的形式來提高其運算效率。

當前,三維體曲率屬性已成為繼相干技術(shù)之后的又一大熱點應(yīng)用技術(shù)[13],尤其在河道及裂縫型油氣儲層預測中效果顯著[14]。但是由于曲率計算過程中要對地震數(shù)據(jù)求兩次微分運算,而微分運算對噪聲較為敏感,因而嚴重影響了曲率屬性的實際應(yīng)用效果。本文提出了一種新的體曲率計算方法,首先利用連續(xù)小波變換方法提取三維地震資料的相位余弦體,然后基于復數(shù)道分析理論,采用CLAERBOUT[15]提出的差分算法計算瞬時頻率和x,y方向的瞬時波數(shù),進一步計算三維地震數(shù)據(jù)體的視傾角,最后應(yīng)用中值濾波后的視傾角數(shù)據(jù)計算各類曲率屬性。模型實驗和實際資料處理結(jié)果表明:①該方法對噪聲不敏感,可以應(yīng)用于低信噪比資料;②該算法計算得到的曲率屬性與C3等相干屬性相比假象更少,能夠真實地反映地下斷層及巖性邊界等地質(zhì)信息,特別適用于陡傾角地層發(fā)育的地區(qū)。

1 技術(shù)原理

1.1 基于小波變換相位余弦屬性提取

實地震道x(t)的小波變換定義為:

(1)

(2)

(3)

(4)

(5)

式中:Re[c(t)]和Im[c(t)]分別表示c(t)的實部和虛部。利用相位屬性可以很方便地計算出相位余弦屬性cos[φ(t)]。

因為直接利用上述公式計算的瞬時屬性與Hilbert方法計算的瞬時屬性是等價的,所以對于含噪信號,在利用小波變換得到WΨ(t,a)之后,本文首先在小波變換域確定有效信號的能量分布空間,然后再利用公式(3)到公式(5)計算瞬時屬性,從而提高了瞬時屬性的信噪比,為后續(xù)視傾角的提取提供了高品質(zhì)基礎(chǔ)資料。

1.2 基于復數(shù)道分析理論計算視傾角

由地下三維空間中任一反射點r(t,x,y)處得到的相位余弦數(shù)據(jù)可以看作是時間標量u(t,x,y),其梯度grad(u)表示該點處反射面上沿著方向矢量所在的法截面截取曲線的一階導數(shù),即反射面沿著不同方向的變化率,可以用下式表示[19-20]:

(6)

式中:px,qy,rt表示反射點r處地震記錄u沿x,y和t方向上的視傾角分量(見圖1)。圖1中,φ為方位角,θ為真傾角。

曲率是描述曲線上任意一點彎曲程度的一種屬性,在數(shù)學上可定義為曲線上某點處角度與弧長變化率之比,也可表示為該點的二階微分[21]:

(7)

式中:α是角度;s是弧長。代入視傾角px,qy,得到沿x方向和y方向的曲率分量:

(8a)

圖1 三維空間中視傾角與傾角和方位角的關(guān)系[5]

(8b)

由此可知,要計算體曲率,需要先求取視傾角數(shù)據(jù)體。有多種求取視傾角的方法,可以通過傾角掃描或梯度結(jié)構(gòu)張量方法計算[22-24],也可以采用復數(shù)道分析理論方法[25-26],利用瞬時頻率和x,y方向的瞬時波數(shù)來計算三維地震數(shù)據(jù)體的視傾角。

其中,瞬時頻率ω,x和y方向的瞬時波數(shù)kx和ky分別為:

(9)

式中:φ為瞬時相位;φx和φy為沿著x和y方向的空變相位;u為輸入的相位余弦數(shù)據(jù);uH為其時間方向的Hilbert變換。

從(9)式可看出,計算瞬時頻率、瞬時波數(shù)時需要對初始地震數(shù)據(jù)及其Hilbert變換求二階偏導數(shù),這類運算對噪聲非常敏感,因此,本文按下式計算瞬時頻率:

(10)

式中:z(t)為初始信號的解析信號;Δt為信號的采樣間隔。

為了提高計算精度,對(10)式表示的方法進行改進,利用中心差分格式,并將虛部運算化為實部,得到瞬時頻率的實現(xiàn)公式為:

u(t)-[u(t+Δt)-u(t-Δt)h(t)}/[u2(t)+

h2(t)]

(11)

式中:u(t)為輸入相位余弦數(shù)據(jù);h(t)為其相應(yīng)的Hilbert變換。

同理可得沿x和y方向的瞬時波數(shù)kx,ky的實現(xiàn)公式:

(12)

視傾角(px,qy)可利瞬時波數(shù)kx,ky和瞬時頻率ω進行計算:

(13)

由于計算過程中有微分運算,因此求得的視傾角數(shù)據(jù)常常存在較大的噪聲和突變點。為解決這個問題,本文應(yīng)用中值濾波的方法對視傾角數(shù)據(jù)進行預處理,提高了傾角屬性的信噪比,從而能夠增強后續(xù)曲率屬性提取的穩(wěn)定性和平面效果。

1.3 曲率屬性提取

在三維空間中,任一點的曲率值可以利用其相鄰道和樣點的視傾角值擬合的空間曲面方程來計算,基于最小二乘逼近原理,N次曲面擬合方程如下[4]:

(14)

當N=2時,將方程系數(shù)置換,得到二次曲面方程:

z(x,y)=ax2+by2+cxy+dx+ey+f

(15)

方程的系數(shù)為:

(16)

因此,可以利用(16)式中的5個系數(shù)計算三維地震數(shù)據(jù)中每一點的曲率。在過三維空間任一點的曲面上有很多種曲率屬性,其中比較常用的曲率屬性有平均曲率、高斯曲率、最大正曲率和最小負曲率等,其數(shù)學含義及地質(zhì)含義見表1。

表1 各類曲率屬性數(shù)學含義及地質(zhì)含義

2 數(shù)值實驗

為了驗證上述技術(shù)對河道和斷層的檢測能力以及對資料信噪比的適應(yīng)性,本文對包含河道和推覆構(gòu)造的推覆體模型進行了實驗。縱測線和橫測線都是401道,縱向185ms,采樣間隔1ms。圖2a是第320條橫測線(Xline320)速度剖面,可以看出推覆體模型含有高角度地層、多條斷層和河道等地質(zhì)信息,可以用來檢測算法對傾斜地層、斷層和小河道的檢測能力。圖2b是應(yīng)用主頻為30Hz的子波得到的正演模型地震剖面,可以看到,由于主頻較低,河道信息在地震記錄上不是很明顯。用常規(guī)C3相干算法和本文提出的算法分別進行了不連續(xù)性檢測,圖2c是本文方法計算的最小負曲率結(jié)果,圖2d是C3相干算法檢測結(jié)果,可以看出,C3相干算法受地層結(jié)構(gòu)影響較為強烈,大傾角地層存在較強的層間不連續(xù)性信息,影響了平面上斷層和河道的檢測。通過對比可以看出,圖2d中a點和b點處的不連續(xù)性是由大傾角地層引起的,不是真實的斷層和巖性邊界,c點和e點是真實的斷層信息。圖2c中消除了傾斜地層的影響,只有c點和e點處的斷層被檢測出來,如實地反映了真實的斷層和河道等地質(zhì)體邊界信息。其中d點處雖然曲率較大也被檢測出來,說明曲率屬性既能檢測斷層,也能檢測褶皺等地質(zhì)信息。

圖2 推覆體模型及其不連續(xù)性檢測效果a 三維模型Xline320線的速度剖面; b 對應(yīng)的地震剖面; c 本文方法計算的最小負曲率100ms處時間切片; d C3相干算法計算結(jié)果100ms處時間切片

為檢驗本文曲率分析技術(shù)的抗噪性,對正演地震記錄加入部分高斯白噪聲,并基于小波變換求取了相位余弦屬性,然后對原始數(shù)據(jù)(圖3a)、含噪數(shù)據(jù)(圖3b)和相位余弦數(shù)據(jù)(圖3c)分別用本文算法和C3相干算法進行不連續(xù)性檢測。圖3d,圖3e和圖3f是C3相干算法計算結(jié)果,圖3g,圖3h和圖3i是本文算法計算結(jié)果。對比發(fā)現(xiàn),C3相干算法對噪聲要敏感得多,無論是含噪數(shù)據(jù)還是相位余弦數(shù)據(jù),檢測結(jié)果中河道邊界都非常模糊,特別是左下方河道與斷層交錯的位置(紅色方框內(nèi)),河道邊界不可辨識。而本文方法基于相位余弦數(shù)據(jù)計算的結(jié)果信噪比較高,在本文如此大的噪聲背景下仍然能夠較好地刻畫出河道的邊界信息。

圖3 基于推覆體模型的本文曲率分析技術(shù)抗噪性實驗a 原始數(shù)據(jù); b 含噪數(shù)據(jù); c 相位余弦數(shù)據(jù); d 原始數(shù)據(jù)C3相干算法計算結(jié)果; e 含噪數(shù)據(jù)C3相干算法計算結(jié)果; f 相位余弦數(shù)據(jù)C3相干算法計算結(jié)果; g 原始數(shù)據(jù)本文算法計算結(jié)果; h 含噪數(shù)據(jù)本文算法計算結(jié)果; i 相位余弦數(shù)據(jù)本文算法計算結(jié)果

3 實際應(yīng)用

利用南海某油田斷層和高陡地層比較發(fā)育的實際地震資料進行試驗處理,分別提取C3相干屬性和曲率屬性。圖4是三維地震數(shù)據(jù)的一條任意測線剖面,左側(cè)地層較平,中間發(fā)育高陡地層,右側(cè)發(fā)育大量斷層。圖5是4100ms處C3相干時間切片和地震任意測線的三維顯示結(jié)果,可以看出,中部高陡地層發(fā)育區(qū)(箭頭指示區(qū))用C3相干算法識別為斷層,所以該區(qū)域顯示出一團強烈的不相干性,無法判斷斷層位置,右側(cè)斷層發(fā)育區(qū)用C3相干算法識別的結(jié)果與實際剖面不太吻合,因為該區(qū)部分地層傾斜角度也比較大。圖6是曲率K1屬性的切片和地震剖面共同顯示的結(jié)果,可以看出,中部陡傾角地層發(fā)育區(qū)(箭頭指示區(qū))曲率切片非常干凈,沒有顯示出不連續(xù)性,和地層實際情況吻合較好,右側(cè)斷層發(fā)育區(qū)切片和剖面在斷點處較一致。說明曲率屬性在陡傾角地層發(fā)育的地區(qū)有著非常好的應(yīng)用效果，而用相干屬性方法則無法達到這個效果。

圖4 實際地震剖面

圖5 4100ms處C3相干時間切片與地震剖面

圖6 4100ms處曲率K1時間切片與地震剖面

4 結(jié)束語

本文利用小波變換計算相位余弦屬性,在此基礎(chǔ)上計算多種曲率屬性和角度屬性,進行斷層、河道等地質(zhì)體邊緣信息的檢測。該項技術(shù)能夠充分利用三維地震資料的空間信息計算地層的三維彎曲程度。與傳統(tǒng)相干算法相比,能夠有效地去除高傾角地層對斷層檢測的影響,得到合理可靠的檢測結(jié)果,特別適用于在陡傾角地層發(fā)育的地區(qū)進行不連續(xù)性檢測。

此外,由于采用了相位余弦計算、視傾角中值濾波等處理技術(shù),該算法具有良好的抗噪性,適用于低信噪比數(shù)據(jù)。當然,在進行曲率分析之前,如果能進行必要的保邊緣去噪處理,檢測效果可能會更好。

致謝:感謝中國石油大學(北京)袁三一博士提供推覆體模型數(shù)據(jù),感謝成都理工大學陳學華教授在曲率認識方面進行的有益探討。

[1] CHOPRA S,MARFURT K J.Curvature attribute applications to 3D surface seismic data[J].The Leading Edge,2007,26(4):404-414

[2] CHOPRA S,MARFURT K J.Emerging and future trends in seismic attributes[J].The Leading Edge,2008,27(3):298-318

[3] LISLE R J.Detection of zones of abnormal strains in structures using Gaussian curvature analysis[J].AAPG Bulletin,1994,78(12):1811-1819

[4] ROBERTS A.Curvature attributes and their application to 3D interpreted horizons[J].First Break,2001,19(2):85-100

[5] AL-DOSSARY S,MARFURT K J.3D volumetric multispectral estimates of reflector curvature and rotation[J].Geophysics,2006,71(5):41-51

[6] CHEN X H,YANG W,HE Z H,et al.The algorithm of 3D multi-scale volumetric curvature and its application[J].Applied Geophysics,2012,9(1):65-72

[7] GAO D L.Integrating 3D seismic curvature and curvature gradient attributes for fracture characterization:methodologies and interpretational implication[J].Geophysics,2013,78(2):O21-O31

[8] 王世星.高精度地震曲率體計算技術(shù)與應(yīng)用[J].石油地球物理勘探,2012,47(6):965-972 WANG S X.High-precision calculation of seismic volumetric curvature attributes and its application[J].Oil Geophysical Prospecting,2012,47(6):965-972

[9] 印興耀,高京華,宗兆云.基于離心窗傾角掃描的曲率屬性提取[J].地球物理學報,2014,57(10):3411-3421 YIN X Y,GAO J H,ZONG Z Y.Curvature attribute based on dip scan with eccentric window[J].Chinese Journal of Geophysics,2014,57(10):3411-3421

[10] QI J,ZHANG B,ZHOU H L,et al.Attribute expression of fault-controlled Karst-Fort Worth Basin,TX[J].Expanded Abstracts of 84thAnnual Internat SEG Mtg,2014:1522-1527

[11] 楊國權(quán),劉延利,張紅文.曲率屬性計算方法研究及效果分析[J].地球物理學進展,2015,30(5):2282-2286 YANG G Q,LIU Y L,ZHAGN H W.The calculation method of curvature attributes and its effect analysis[J].Process in Geophysics,2015,30(5):2282-2286

[12] 李福強,賀振華,文曉濤,等.改進的曲率計算方法及其效果分析[J].石油物探,2012,51(2):146-150 LI F Q,HE Z H,WEN X T,et al.An improved curvature calculation algorithm and its effect analysis[J].Geophysical Prospecting for Petroleum,2012,51(2):146-150

[13] KLEIN P,RICHARD L,JAMES H.3D curvature attributes:a new approach for seismic interpretation[J].First Break,2008,26(4):105-112

[14] BUCK D M,ALAM A,TAYLOR J D.Fractured reservoir prediction from 3D seismic volumetric curvature and low frequency imaging[J].Expanded Abstracts of 77thAnnual Internat SEG Mtg,2007:422-426

[15] CLAERBOUT J F.Fundamentals of geophysical data processing[M].Palo Alto:McGraw-Hill Book Company,1976:1-269

[16] GAO J H,DONG X,WANG W B,et al.Instantaneous parameters extraction via wavelet transform[J].IEEE Transactions on Geoscience and Remote Sensing,1999,37(2):867-870

[17] 高靜懷,吳茜,陳文超,等.小波變換域地震資料瞬時頻率分析方法[J].石油物探,2007,46(6):534-540 GAO J H,WU Q,CHEN W C,et al.Instantaneous frequency analysis of seismic data in wavelet transform domain[J].Geophysical Prospecting for Petroleum,2007,46(6):534-540

[18] 高靜懷,汪文秉,朱光明.小波變換與信號瞬時特征分析[J].地球物理學報,1997,40(6):821-833 GAO J H,WANG W B,ZHU G M.Wavelet transform and instantaneous attributes analysis[J].Chinese Journal of Geophysics,1997,40(6):821-833

[19] 楊威.曲率屬性分析及其應(yīng)用[D].成都:成都理工大學,2011 YANG W.Analysis of curvature attribute and its application[D].Chengdu:Chengdu University of Technology,2011

[20] 楊威,賀振華,陳學華.結(jié)構(gòu)方位濾波在體曲率屬性中的應(yīng)用[J].石油物探,2011,50(1):27-32 YANG W,HE Z H,CHEN X H.Application of structure-oriented filtering to volumetric curvature calculation[J].Geophysical Prospecting for Petroleum,2011,50(1):27-32

[21] THOMAS G B,FINNEY R L.Calculus and analytical geometry[M].Upper Saddle River:Addison Wesley Publishing Company,1972:1-506

[22] MARFURT K J,KIRLIN R L,FARMER S,et al.3-D seismic attributes using a semblance-based coherency algorithm[J].Geophysics,1998,63(4):1150-1165

[23] WANG X K,CHEN W C,GAO J H,et al.Estimate dip by combining gradient structure tensor and multi-window[J].Expanded Abstracts of 85thAnnual Internat SEG Mtg,2015:1722-1727

[24] WANG X K,YANG C C,CHEN W C,et al.Estimating the dip by constructing gradient structure tensor on instantaneous phase [J].Expanded Abstracts of 84thAnnual Internat SEG Mtg,2014:1580-1584

[25] LUO Y,HIGGS W G,KOWALIK W S.Edge detection and stratigraphic analysis using 3D seismic data[J].Expanded Abstracts of 66thAnnual Internat SEG Mtg,1996:324-327

[26] BARNES A E.Theory of 2-D complex seismic trace analysis[J].Geophysics,1996,61(1):264-272

(編輯:朱文杰)

A 3D curvature attribute analysis method with excellent anti-noise property suitable for high steep formation

WANG Qingzhen,JIANG Xiudi,WENG Bin,DING Jicai,CHEN Jianjun

(CNOOCResearchInstitute,Beijing100028,China)

Curvature attributes is useful in delineating faults,fractures,flexures and folds.However,in order to get curvature attributes,it is necessary to calculate the second-order derivative of seismic data.The derivative calculation process is sensitive to noise and the results are unstable and unreliable.In this paper we present a new 3D curvature attribute analysis method with excellent anti-noise property.Firstly,wavelet transform is used to extract phase-cosine attribute from seismic data.Then apparent dips are calculated based on complex seismic trace theory and filtered using media filter.Finally,many types of curvature attributes with high S/N data are got by combining apparent dips in both directions.Applications of this method to both synthetic and real seismic data show that this method has two advantages.First,it is insensitive to random noise and adapt to low S/N data.Second,it can effectively suppress the pseudo boundary information caused by discontinuity between the tilted strata,so this method is more suitable for steep dipping strata compared with C3 method.

fault detection,curvature,phase cosine,coherence,wavelet transform

2016-03-17;改回日期:2016-10-11。

王清振(1983—),男,碩士,主要從事地震資料解釋、屬性分析、反演等技術(shù)研究及相關(guān)軟件研發(fā)工作。

國家科技重大專項(2016ZX05024-001)資助。

P631

A

1000-1441(2017)04-0559-08

10.3969/j.issn.1000-1441.2017.04.012

This research is financially supported by National Key Project of Science and Technology (2016ZX05024-001).