基于高維連續(xù)小波變換的地震資料不連續(xù)性檢測方法研究

王曉凱， 高靜懷， 陳文超， 楊長春， 朱振宇

1 西安交通大學(xué)計(jì)算地球物理科學(xué)系， 西安 710049 2 中國科學(xué)院地質(zhì)與地球物理研究所， 北京 100029 3 西安交通大學(xué)電子與信息工程學(xué)院， 西安 710049 4 中海油研究總院， 北京 100027

?

基于高維連續(xù)小波變換的地震資料不連續(xù)性檢測方法研究

王曉凱1,2,3， 高靜懷1,3*， 陳文超3， 楊長春2， 朱振宇4

1 西安交通大學(xué)計(jì)算地球物理科學(xué)系， 西安 710049 2 中國科學(xué)院地質(zhì)與地球物理研究所， 北京 100029 3 西安交通大學(xué)電子與信息工程學(xué)院， 西安 710049 4 中海油研究總院， 北京 100027

斷層等不連續(xù)性結(jié)構(gòu)的檢測在地震勘探中具有重要的意義.本文簡單介紹了高維連續(xù)小波變換的基本理論，包括：高維連續(xù)小波變換的定義、高維小波變換系數(shù)的切片以及快速實(shí)現(xiàn)方法等，利用典型的合成信號(hào)說明了高維連續(xù)小波變換具有更好的方向選擇性(相比于常用的高維張量積小波變換)，將高維連續(xù)小波變換引入地震資料不連續(xù)性檢測，并提出了利用小尺度高維連續(xù)小波變換系數(shù)檢測地震資料不連續(xù)性的方法，同時(shí)給出了方法的實(shí)現(xiàn)流程.合成信號(hào)以及實(shí)際三維地震數(shù)據(jù)處理效果驗(yàn)證了本文所提方法的有效性.

高維連續(xù)小波變換； 快速傅里葉變換； 不連續(xù)性

1 引言

隨著油氣勘探從構(gòu)造油氣藏向地層-巖性等隱蔽油氣藏的轉(zhuǎn)移，對(duì)斷層、裂縫以及地質(zhì)體邊緣等不連續(xù)性結(jié)構(gòu)的刻畫有著更加重要的意義.由于目前三維地震勘探的大規(guī)模應(yīng)用以及高密度地震數(shù)據(jù)采集逐步進(jìn)入工業(yè)生產(chǎn)，三維疊后數(shù)據(jù)量巨大，利用人工提取斷層等不連續(xù)結(jié)構(gòu)不僅工作量巨大，而且受主觀因素的影響造成結(jié)果的不確定性.

相干體技術(shù)是近些年來發(fā)展起來的一項(xiàng)強(qiáng)大的地震屬性解釋技術(shù)，可大大縮短三維地震數(shù)據(jù)的解釋周期(苑書金, 2007)、確定斷層的空間分布以及目的層的巖性變化、判斷數(shù)據(jù)采集以及處理的質(zhì)量(Chopra and Marfurt, 2007).Bahorich等人于1995年提出了基于歸一化道互相關(guān)的相干體算法(第一代相干體算法)，雖然該項(xiàng)技術(shù)在計(jì)算標(biāo)準(zhǔn)道相干性時(shí)只用了三道地震數(shù)據(jù)，但實(shí)際應(yīng)用效果立刻在業(yè)界引起了巨大的反響(Bahorich and Farmer, 1995; Bahorich et al., 1995).陸文凱等人將第一代相干體算法進(jìn)行發(fā)展，提出基于超道技術(shù)和高階統(tǒng)計(jì)量的相干體算法(Lu et al., 2005).Marfurt等人于1998年提出基于相似度量的相干體算法(第二代相干體算法)，該算法利用分析窗內(nèi)平均道的能量與所有道能量的比值來度量相干性(Marfurt et al., 1998).該方法提高了相干體的信噪比，但增加了計(jì)算量，同時(shí)降低了空間分辨率.Gersztenkorn與Marfurt于1999年將第一代相干體算法推廣到任意J道(Gersztenkorn and Marfurt, 1996, 1999; Marfurt et al., 1999)，提出基于特征結(jié)構(gòu)的相干體算法(第三代相干體算法).第三代相干體算法運(yùn)算量實(shí)際應(yīng)用效果好，抗噪性能強(qiáng)，分辨率高，但是其運(yùn)算復(fù)雜程度高，計(jì)算耗時(shí)，限制了其在工業(yè)界的大規(guī)模應(yīng)用.多位學(xué)者采用各種加速方法來提高第三代相干體算法的計(jì)算速度，取得了明顯的加速效果(孫夕平和杜世通, 2000; 葉增爐和何建軍, 2006; Wang et al., 2012).Marfurt與Kirlin將第三代相干體算法與雷達(dá)信號(hào)處理中的空間譜估計(jì)技術(shù)相結(jié)合，提出了基于MUSIC的相干體算法(Marfurt and Kirlin, 2000; 宋建國等, 2005).Cohen和Coifman于2002年提出局部結(jié)構(gòu)熵算法(Cohen and Coifman, 2002)，從而來度量地震數(shù)據(jù)體中的不連續(xù)性結(jié)構(gòu)，大大提高了計(jì)算效率.Randen等人提出基于梯度結(jié)構(gòu)張量技術(shù)來估計(jì)地震數(shù)據(jù)的紋理屬性(Randen et al., 2000; Chopra and Marfurt, 2007; 張軍華等, 2007).以此為基礎(chǔ)，多位學(xué)者提出基于梯度結(jié)構(gòu)張量的相干算法(Randen et al., 2000; Bakker, 2003).單一的相干屬性有時(shí)不能滿足解釋人員的需要，因此也有諸多許多學(xué)者提出將相干體技術(shù)與瞬時(shí)屬性聯(lián)合應(yīng)用提高地質(zhì)目標(biāo)體的分辨率(王西文等, 2002a, 2002b; 張軍華等, 2004; Zhou et al., 2007; 苑書金, 2007; 周艷輝和高靜懷, 2007).

自從法國女工程師Morlet在處理地震資料時(shí)首次提出小波變換的思想以后(Antoine et al., 2004)，一維小波變換在地震信號(hào)處理及解釋領(lǐng)域取得了相當(dāng)多成功的應(yīng)用(Chopra and Marfurt, 2007).1987年，Murenzi在與Grossmann以及Daubichies等人討論以后，在其博士論文中首次提出兩維連續(xù)小波變換并對(duì)其進(jìn)行了詳細(xì)的討論(Murenzi, 1990).隨后，高維連續(xù)小波變換在輪廓檢測(張長江等, 2006)、字符識(shí)別(Antoine et al., 1995a; Murenzi et al., 1999)、目標(biāo)辨識(shí)(Antoine et al., 1995b)、醫(yī)學(xué)圖像處理(Choi et al., 2001)、模式識(shí)別(Antoine et al., 1999)、天體圖像處理及目標(biāo)識(shí)別(Antoine et al., 2002)、探地雷達(dá)信號(hào)處理(陳文超等, 2000, 2003)等各個(gè)領(lǐng)域均取得了成功應(yīng)用.在高維連續(xù)小波變換中，對(duì)母小波的操作有平移、伸縮和旋轉(zhuǎn)，這三種操作由平移因子、尺度以及角度所控制.常用的二維Morlet小波具有較好的方向選擇性，因此被廣泛采用.由于高維連續(xù)小波變換的結(jié)果為四維的，數(shù)據(jù)量巨大且不易顯示，Murenzi初步探討了二維連續(xù)小波變換系數(shù)的快速實(shí)現(xiàn)算法(Murenzi et al., 1997).

二維連續(xù)小波變換在地震信號(hào)處理中也取得了許多應(yīng)用，例如：面波去除(劉財(cái)?shù)? 2003; 李媛媛等, 2004; 張華等, 2007)、噪聲衰減(高振山等, 1998; 楊立強(qiáng)等, 2005)、提高分辨率(吳愛弟和牟永光, 1997)等.Mallat等人提出利用一維小波變換系數(shù)模的極大值來檢測信號(hào)的奇異性(Mallat and Hwang, 1992).以Mallat等人提出的方法為基礎(chǔ)，Li Chun-Feng等人提出利用基于一維小波變換的Holder指數(shù)奇異性屬性刻畫地層邊界等(Li and Liner, 2008).Ouillon與Sornette等人首次將二維連續(xù)小波變換用于分析斷層等奇異性結(jié)構(gòu)，提出利用各向異性小波變換結(jié)果的最優(yōu)系數(shù)來檢測各個(gè)尺度斷層的走向(Ouillon et al., 1995, 1996).他們采用二維各向異性Mexcian Hat小波作為母小波，在給定的分辨率下保留模最大的小波系數(shù)，這樣被保留下的小尺度系數(shù)能夠反映較小斷層的結(jié)構(gòu)和走向，大尺度系數(shù)能夠反映更大區(qū)域內(nèi)斷層的結(jié)構(gòu)和走向.隨后，Darrozes與Gaillot等人將上述方法進(jìn)行改進(jìn)，利用歸一化最優(yōu)各向異性小波系數(shù)來檢測多尺度的斷層結(jié)構(gòu)，并采用玫瑰方向圖來判斷各個(gè)分辨率下的目標(biāo)區(qū)域的斷層走向(Darrozes et al., 1997; Gaillot et al., 1997, 1999).但上述分析方法采用的均為Mexican小波，方向選擇性較差.Bouchereau首次利用二維Morlet小波作為母小波，將兩維連續(xù)小波變換用于勘探地震資料的斷層檢測.她采用二維Morlet小波這種具有方向選擇性的小波作為母小波，并與自適應(yīng)濾波技術(shù)相結(jié)合，對(duì)斷層面檢測進(jìn)行了研究初探(Bouchereau, 1997).

本文直接以三維連續(xù)小波變換為基礎(chǔ)檢測地震資料的不連續(xù)性，安排如下：首先介紹高維連續(xù)小波變換的基本理論及快速實(shí)現(xiàn)方法，然后給出基于高維連續(xù)小波變換的地震資料不連續(xù)性檢測方法及實(shí)現(xiàn)流程，接著將本文所提出的方法用于合成信號(hào)及實(shí)際資料處理，最后給出結(jié)論及展望.

2 高維連續(xù)小波變換的基本理論及快速實(shí)現(xiàn)

2.1 二維連續(xù)小波變換

假設(shè)ψ(x)是二維母小波，對(duì)母小波的基本操作有三種：平移、伸縮以及旋轉(zhuǎn).經(jīng)過這三種操作的小波為

(1)

其中，x=(x,y)T為二維變量，平移由平移因子b(b=(bx,by)T)控制，伸縮由尺度因子a控制，旋轉(zhuǎn)由旋轉(zhuǎn)因子(即角度)θ控制，r-θ()表示對(duì)變量逆時(shí)針旋轉(zhuǎn)角度θ.對(duì)二維母小波的三種操作示意如圖1所示.

設(shè)f(x)是二維信號(hào)，則f(x)的二維連續(xù)小波變換(Two-DimensionalContinuousWaveletTransform, 2DCWT)定義為

CWT_2D(b,a,θ)=〈f,ψb,a,θ〉

(2)

(3)

(4)

(5)

(6)

通常情況下，取參數(shù)σ=(0,σy)T后，式(5)和式(6)簡化為

(7)

選取ε=5，采用參數(shù)σ=(0,5.33)T時(shí)二維Morlet小波的實(shí)部、虛部、包絡(luò)及二維傅里葉振幅譜如圖2所示.

當(dāng)所采用的母小波為不具有方向選擇性時(shí)，例如二維Mexican Hat小波，旋轉(zhuǎn)操作不起作用，這時(shí)二維信號(hào)f(x)的2DCWT退化為

CWT_2D(b,a)=〈f,ψb,a〉

(9)

2DCWT的系數(shù)共有四個(gè)參數(shù)，也就是變換的結(jié)果是四維的，那么系數(shù)無法直接顯示.為了解決這些問題，固定2DCWT系數(shù)中的兩個(gè)參數(shù)，這樣就只有兩個(gè)參數(shù)變化，可以得到2DCWT系數(shù)的二維切片，便于顯示.二維切片共有6種，最為常用且意義比較明確的有空間位置切片以及尺度-角度切片：

(a)空間位置切片：固定尺度和旋轉(zhuǎn)角度，變化空間位置.利用空間位置切片可以分析某類特征在整個(gè)二維平面的分布情況.空間位置切片在諸如圖像處理、湍流分析、斷層檢測等各個(gè)方面取得了成功的應(yīng)用.

圖1 對(duì)母小波的各種操作(a) 空間域母小波的包絡(luò)； (b) 二維傅里葉域母小波的模值； (c) 空間域經(jīng)過平移的小波包絡(luò)； (d) 二維傅里葉域經(jīng)過平移母小波的模值； (e) 空間域經(jīng)過伸縮的小波包絡(luò)； (f) 二維傅里葉域經(jīng)過伸縮母小波的模值； (g) 空間域經(jīng)過旋轉(zhuǎn)的小波包絡(luò)； (h) 二維傅里葉域經(jīng)過旋轉(zhuǎn)母小波的模值.Fig.1 Sketches showing operations on mother wavelets(a) Envelop of mother wavelet (in space domain)； (b) Modulus of mother wavelet (in 2D Fourier transform domain)； (c) Envelop of shifted wavelet (in space domain)； (d) (b) Modulus of shifted wavelet (in 2D Fourier transform domain)； (e) Envelop of dilated wavelet (in space domain)； (f) Modulus of dilated wavelet (in 2D Fourier transform domain)； (e) Envelop of rotated wavelet (in space domain)； (f) Modulus of rotated wavelet (in 2D Fourier transform domain).

(b)尺度角度切片：固定空間位置，變化尺度和旋轉(zhuǎn)角度.利用尺度-角度切片，可以觀察固定空間位置的所有尺度和角度下的2DCWT系數(shù)，在分析小波函數(shù)的尺度特性以及角度分辨率方面有著重要的作用.采用如下一個(gè)合成信號(hào)(帶有衰減的二維復(fù)合平面波)，公式為

(10)

其中cn表示第n個(gè)平面波的振幅，kn表示第n個(gè)平面波的傳播方向及頻率，In表示衰減方向及衰減因子.該二維信號(hào)包絡(luò)如圖3a所示，原點(diǎn)處的2DCWT系數(shù)尺度角度切片如圖3b所示，其中vx=cos(θ)/a，vy=sin(θ)/a，可以在此切片上清晰的看出該二維信號(hào)在原點(diǎn)附近有六個(gè)平面波分量(對(duì)應(yīng)于六個(gè)峰值).

由2DCWT的定義(如式(2)所示)，計(jì)算一個(gè)2DCWT系數(shù)需要一個(gè)二重積分，對(duì)于離散情況而言，需要一個(gè)二維累加求和.上面已經(jīng)提到，2DCWT的結(jié)果是四維的，因此，直接利用式(2)計(jì)算復(fù)雜度很高.在波數(shù)域式(2)可簡寫為

CWT(f;b,a,θ)=

(11)

CWT_2D(b,a,θ)=

(12)

其中2DFT和2DIFT2分別表示二維傅里葉變換和逆二維傅里葉變換.由于二維傅里葉變換有快速算法的支持，因此，利用上式可以快速計(jì)算2DCWT系數(shù)的空間位置切片.2DCWT系數(shù)空間位置切片的快速實(shí)現(xiàn)流程圖見圖4所示.由于有快速傅里葉變換的支持，計(jì)算一個(gè)2DCWT系數(shù)空間位置切片的運(yùn)算復(fù)雜度由O(N4)降低為O(N2logN).

2.2 三維連續(xù)小波變換

目前工業(yè)界內(nèi)使用的成果地震數(shù)據(jù)多為三維的，因此三維變換更適合處理三維地震資料，所以有必要將二維連續(xù)小波變換推廣到三維.一維連續(xù)小波變換推廣到二維連續(xù)小波變換時(shí)，平移因子由一維變成二維向量，并且增加了一個(gè)旋轉(zhuǎn)因子.將二維連續(xù)小波變換推廣到三維時(shí)，對(duì)母小波的操作仍是三種：平移、伸縮及旋轉(zhuǎn)，但是平移因子為三維向量，旋轉(zhuǎn)操作里包含兩個(gè)變量，即沿傾角方向的旋轉(zhuǎn)以及沿方位角的旋轉(zhuǎn).設(shè)f(x)為三維信號(hào)，ψ(x)為三維母小波，則的三維連續(xù)小波變換(Three-Dimensional Continuous Wavelet Transform, 3DCWT)定義為

圖2 ε=5，σ=(0,5.33)T時(shí)二維Morlet小波(a) 實(shí)部； (b) 虛部； (c) 包絡(luò)； (b) 振幅譜.Fig.2 2D Morlet with ε=5 and σ=(0,5.33)T(a) Real-part； (b) Imaginary-part； (c) Envelop； (d) Amplitude spectrum.

圖3 二維合成信號(hào)及其尺度角度切片(a) 含有六個(gè)平面波的合成信號(hào)(振幅)； (b) 原點(diǎn)處的尺度角度切片.Fig.3 2D synthetic signals and scale-angle slices of 2DCWT coefficients(a) Synthetic signal which contains six plane waves；(b) One of 2DCWT coefficients scale-angle slices at origin.

圖4 2DCWT系數(shù)空間位置切片的快速實(shí)現(xiàn)流程Fig.4 Flow chart of fast realization of 2DCWT coefficients (position slice)

CWT_3D(b,a,θ,φ)=〈f,ψb,a,(θ,φ)〉

(13)

常用的三維小波有三維Mexican-hat小波、三維Morlet小波以及三維Cauchy小波.三維Mexican-hat小波不具有方向選擇性，因此我們只介紹具有方向選擇性的三維Morlet小波.三維Morlet小波的空間域及三維傅里葉域表達(dá)式為

(14)

(15)

由于3DCWT的系數(shù)為六維，也存在兩維連續(xù)小波變換的問題.與兩維連續(xù)小波變換相同，也有兩種常用系數(shù)“切片”：空間位置“切片”(a、θ及φ固定)以及尺度方向“切片”(b固定).空間位置“切片”可以通過類似于圖4的流程利用FFT技術(shù)快速實(shí)現(xiàn)，如圖5所示.

圖5 3DCWT系數(shù)空間位置切片的快速實(shí)現(xiàn)流程Fig.5 Flow chart of fast realization of 3DCWT coefficients (position slice)

3 基于高維連續(xù)小波變換的地震資料橫向不連續(xù)性檢測方法

一維連續(xù)小波變換及其模極大值已被廣泛地用于一維信號(hào)的奇性檢測，其小尺度系數(shù)模能夠精確地定位奇性位置(Mallat, 2008).目前業(yè)內(nèi)使用的地震資料均為二維/三維，直接利用一維連續(xù)小波變換不能充分利用道與道之間的相關(guān)性來檢測地震資料的橫向不連續(xù)性，因此有必要采用高維變換來檢測地震資料的橫向不連續(xù)性.

3.1 高維連續(xù)小波變換角度分辨率分析

不具有方向選擇性高維Mexican-Hat小波不能區(qū)分地震資料的橫向不連續(xù)性和縱向不連續(xù)性.具有一定方向選擇性的二維張量積小波變換具有三個(gè)方向，但所具有的三個(gè)方向的方向選擇性較差，這樣會(huì)使傾斜方向的不連續(xù)性在水平細(xì)節(jié)分量、垂直細(xì)節(jié)分量上均有體現(xiàn).合成信號(hào)如圖6a所示，包含有垂直、水平及傾斜三個(gè)方向的不連續(xù)性.二維張量積小波變換最小尺度的水平、垂直以及傾斜分量分別如圖6b、圖6c以及圖6d所示，傾斜方向的不連續(xù)性不但在傾斜分量上(圖6d)有所體現(xiàn)，在水平分量及垂直分量上也有所體現(xiàn).因此，在實(shí)際應(yīng)用中二維張量積小波變換僅能區(qū)分水平和垂直兩個(gè)方向的不連續(xù)性，不能區(qū)分介于垂直和水平方向之間的傾斜方向不連續(xù)性，這對(duì)檢測地震資料橫向不連續(xù)性是遠(yuǎn)遠(yuǎn)不夠的.

圖6 合成測試信號(hào)的二維張量積小波分解結(jié)果(a) 合成的測試信號(hào)； (b) 二維張量積小波變換-水平分量； (c) 二維張量積小波變換-垂直分量； (d) 二維張量積小波變換-傾斜分量.Fig.6 2D wavelet decomposition of synthetic signal(a) Synthetic test signal； (b) Coefficients of 2D tensor-product wavelet transform (in horizontal direction)； (c) Coefficients of 2D tensor-product wavelet transform (in vertical direction)；(d) Coefficients of 2D tensor-product wavelet transform (in oblique direction).

圖7 合成測試信號(hào)的2DCWT結(jié)果(小尺度)(a) -45°；(b) 0°；(c) 45°度；(d) 90°.Fig.7 2DCWT coefficients of synthetic signal (small scale)(a) -45 degree； (b) 0 degree； (c) 45 degree； (d) 90 degree.

采用具有方向選擇性的二維Morlet小波作為母小波，對(duì)合成信號(hào)(圖6a)做2DCWT，小尺度下-45°、0°、45°、90°的2DCWT系數(shù)如圖7a、圖7b、圖7c以及 圖7d所示，傾斜方向的不連續(xù)性只出現(xiàn)在圖7a及圖7c中，并不影響其他方向變換的結(jié)果.因此，相比于二維張量積小波變換，2DCWT(采用具有方向選擇性的母小波，例如二維Morlet小波)具有很好的方向選擇性，更適合檢測指定方向的不連續(xù)性.由于常用的疊后地震數(shù)據(jù)為三維的，因此，本節(jié)采用具有方向選擇性的三維Morlet小波作為母小波，利用小尺度3DCWT系數(shù)的模極大值來刻畫三維地震資料的橫向不連續(xù)性.

3.2 原理介紹

一維連續(xù)小波變換的小尺度系數(shù)模極大值來可被用來定位/度量信號(hào)中的不連續(xù)性(Mallat, 2008).在小尺度下，小波具有很好的時(shí)間分辨率，因此在度量不連續(xù)性時(shí)均在小尺度下進(jìn)行.一維連續(xù)小波變換中對(duì)母小波的操作為平移和伸縮，在三維連續(xù)小波變換中，對(duì)母小波的操作增加了旋轉(zhuǎn)操作(包含傾角和方位角)，另外平移因子變?yōu)橐粋€(gè)高維矢量.因此，在利用三維連續(xù)小波變換度量地震資料不連續(xù)性時(shí)，固定在小尺度下，需要對(duì)傾角及方位角進(jìn)行掃描以得到小尺度系數(shù)模極大值來度量/定位地震資料的不連續(xù)性，公式為

(16)

由于斷層等不連續(xù)性結(jié)構(gòu)表現(xiàn)為空間方向上的不連續(xù)性，因此在角度掃描時(shí)可以限制角度搜索范圍.3.3 實(shí)現(xiàn)流程

采用小波系數(shù)模極大值來度量地震資料的不連續(xù)性時(shí)，由于各處振幅不相同，因此會(huì)造成強(qiáng)振幅區(qū)域不連續(xù)度量較強(qiáng)，而弱振幅區(qū)域不連續(xù)性度量較弱.因此我們利用瞬時(shí)相位數(shù)據(jù)體作為三維連續(xù)小波變換的輸入，來檢測三維地震資料的不連續(xù)性.將利用高維連續(xù)小波變換檢測地震資料橫向不連續(xù)性的步驟總結(jié)如下：

步驟1：記這個(gè)三維地震數(shù)據(jù)體為f(x,y,t).為簡單起見，用向量x表示三維地震數(shù)據(jù)體中某點(diǎn)的坐標(biāo)(即x=(x,y,t)T)，那么三維地震數(shù)據(jù)體可簡單記為f(x)，其中，x,y表示位置坐標(biāo)，t表示時(shí)間軸.對(duì)該三維地震數(shù)據(jù)體進(jìn)行保邊緣濾波(王偉等，2012； Zhang Bo，2016)等預(yù)處理.

步驟2：對(duì)三維地震數(shù)據(jù)體f(x)每道進(jìn)行Hilbert變換，得到復(fù)地震道F(x),并在此基礎(chǔ)上提取瞬時(shí)相位得到三維瞬時(shí)相位數(shù)據(jù)體Phase(x)，進(jìn)而獲得三維相位余弦數(shù)據(jù)體PhC為

(17)

步驟3：由于現(xiàn)有的計(jì)算機(jī)內(nèi)存有限，加之計(jì)算過程中需要大量的內(nèi)存，因此采用分塊的辦法減少對(duì)內(nèi)存的依賴性.將步驟2生成的三維相位余弦數(shù)據(jù)體進(jìn)行分為N個(gè)子塊，每個(gè)子塊記為PhC_Block(n,x)，分塊時(shí)，子塊各個(gè)方向上的點(diǎn)數(shù)盡量選擇2的整數(shù)次冪，有利于使用FFT來減少計(jì)算量.

步驟4：分別對(duì)每個(gè)三維相位余弦數(shù)據(jù)體子塊PhC_Block(n,x)進(jìn)行三維連續(xù)小波變換得到其三維連續(xù)小波變換系數(shù)CWT_3D(n,x,asmall,θ,φ)為

(18)

上述三維連續(xù)小波變換可以通過圖5所示的流程快速實(shí)現(xiàn).

步驟5：利用每個(gè)三維相位余弦數(shù)據(jù)體子塊的三維連續(xù)小波變換系數(shù)CWT_3D(n,x,asmall,θ,φ)來進(jìn)行三維地震資料的不連續(xù)性檢測.利用小尺度三維連續(xù)小波變換系數(shù)模來確定地震資料不連續(xù)性的方向及度量不連續(xù)性，公式為

(19)

步驟6：將生成的各個(gè)分塊不連續(xù)性數(shù)據(jù)體Discon(n,x)進(jìn)行拼接并輸出，得到整個(gè)地震資料的不連續(xù)性數(shù)據(jù)體，完成檢測.

4 模型及實(shí)際資料處理

首先對(duì)一個(gè)三維合成地震數(shù)據(jù)進(jìn)行測試.該三維模型含有三個(gè)水平反射界面，另外在對(duì)角線方向有一較大斷層，水平方向和垂直方向各有兩個(gè)較弱斷層，采用褶積模型產(chǎn)生三維合成數(shù)據(jù).圖8a是合成三維信號(hào)的時(shí)間切片，紅色箭頭所指分別為水平及垂直方向的兩個(gè)弱斷層，圖8c為相干體(第三代相干體算法)的時(shí)間切片(時(shí)間與圖8a對(duì)應(yīng))，圖8e為采用本方法檢測不連續(xù)性結(jié)果的時(shí)間切片(時(shí)間與圖8a對(duì)應(yīng)).圖8b是合成三維信號(hào)的一條測線，能夠看出有三個(gè)水平反射波，紅色箭頭所指為兩個(gè)較弱斷層，圖8d為相干體的一條測線(測線號(hào)與圖8b對(duì)應(yīng))，圖8f為采用本方法檢測不連續(xù)性結(jié)果的一條測線(測線號(hào)與圖8b對(duì)應(yīng)).上述兩組對(duì)比結(jié)果表明，相比較于相干體的結(jié)果，采用本文提出的方法，不但能夠刻畫明顯的不連續(xù)性，也能夠刻畫弱連續(xù)性.

圖8 合成三維地震信號(hào)橫向不連續(xù)性檢測結(jié)果(a) 合成三維信號(hào)的某時(shí)間切片； (b) 合成三維信號(hào)的某測線； (c) 相干體的對(duì)應(yīng)時(shí)間切片； (d) 相干體的對(duì)應(yīng)測線； (e) 本文方法結(jié)果的對(duì)應(yīng)時(shí)間切片； (f) 本文方法結(jié)果的對(duì)應(yīng)測線.Fig.8 Discontinuity detecting results of synthetic 3D seismic data(a) A time slice of synthetic 3D signal； (b) A line of synthetic 3D signal； (c) Detecting result of coherence algorithm (corresponding to (a))； (d) Detecting result of coherence algorithm (corresponding to (b))； (e) Detecting result of proposed method (corresponding to (a)); (f) The detecting result of proposed method (corresponding to (b)).

圖9 實(shí)際三維地震資料不連續(xù)性檢測結(jié)果時(shí)間切片1(a) 商業(yè)軟件； (b) 本文方法.Fig.9 One slice of discontinuity detecting results of real field-data(a) Commercial software. (b) Proposed method in this study.

然后我們將本文提出的基于高維連續(xù)小波變換的地震資料橫向不連續(xù)性檢測方法用于某油田的一塊三維海上實(shí)際地震資料.采用商業(yè)軟件得到的相干體結(jié)果時(shí)間切片1如圖9a所示，而采用本文方法結(jié)果的時(shí)間切片如圖9b所示.本文方法刻畫斷層更為清晰，另外弱不連續(xù)性結(jié)構(gòu)檢測能力強(qiáng)于商業(yè)軟件，如紅色方框區(qū)域內(nèi)的小河道及弱不連續(xù)性結(jié)構(gòu).在另一個(gè)時(shí)間切片2上，本文方法結(jié)果如圖10b所示，相比于上商業(yè)軟件結(jié)果(如圖10a所示)，本文方法能夠清晰的刻畫東北區(qū)域的斷層(藍(lán)色橢圓內(nèi))以及北部的細(xì)小河道(紅色橢圓內(nèi)).

圖10 實(shí)際三維地震資料不連續(xù)性檢測結(jié)果時(shí)間切片2(a) 商業(yè)軟件； (b) 本文方法.Fig.10 Another slice of discontinuities detecting results of real field-data(a) Commercial software; (b) Proposed method.

5 結(jié)論

本文在介紹了高維連續(xù)小波變換基本理論的基礎(chǔ)上，提出了基于高維連續(xù)小波變換的地震資料不連續(xù)性檢測方法.該方法以瞬時(shí)相位數(shù)據(jù)體為基礎(chǔ)數(shù)據(jù)，并利用小尺度高維連續(xù)小波變換系數(shù)模值來度量不連續(xù)性.合成及實(shí)際地震資料的處理結(jié)果表明：相比較于常用商業(yè)軟件，本文所提方法的結(jié)果對(duì)細(xì)微不連續(xù)性結(jié)構(gòu)的檢測能力更強(qiáng).

Antoine J P, Demanet L, Hochedez J F, et al. 2002. Application of the 2-D wavelet transform to astrophysical images.PhysicaliaMag., 24(1): 93-116.

Antoine J P, Murenzi R, Vandergheynst P. 1999. Directional wavelets revisited: Cauchy wavelets and symmetry detection in patterns.AppliedandComputationalHarmonicAnalysis, 6(3): 314-345.

Antoine J P, Murenzi R, Vandergheynst P, et al. 2004. Two-Dimensional Wavelets and Their Relatives. Cambridge: Cambridge University Press.Antoine J P, Vandergheynst P, Bouyoucef K, et al. 1995a. Alternative representations of an image via the 2D wavelet transform: Application to character recognition.∥Proceedings of the SPIE 2488, Visual Information Processing IV. Orlando, FL, United States: SPIE.

Antoine J P, Vandergheynst P, Bouyoucef K, et al. 1995b. Target detection and recognition using two-dimensional isotropic and anisotropic wavelets.∥Proceedings of the SPIE 2485, Automatic Object Recognition V. Orlando, FL, United States: SPIE.Bahorich M, Farmer S. 1995. 3-D seismic discontinuity for faults and stratigraphic features: The coherence cube.TheLeadingEdge, 14(10): 1053-1058.

Bahorich M S, Lopez J, Haskell N L, et al. 1995. Stratigraphic and structural interpretation with 3-D coherence.∥SEG Technical Program Expanded Abstracts. SEG, 97-100.

Bakker P. 2003. Image structure analysis for seismic interpretation[Ph. D]. Dutch: Technische Universiteit Delft. Bouchereau E B. 1997. Analyse d′images par transformées en ondelettes[Ph. D]. Fran?ais: Université Joseph Fourier-Grenoble I. Chen W C, Shi Z S, Wang W B, et al. 2000. Suppressing the direct wave noise in GPR data using the 2-D continue directional wavelets.ChineseJournalofRadioScience(in Chinese), 15(3): 352-357.

Chen W C, Wang W B, Zhao R C. 2003. Ground penetrating radar (GPR) wavelet constructing and application to high resolution processing of GPR data.CoalGeology&Exploration(in Chinese), 31(5): 46-48.

Choi Y, Koo J Y, Lee N Y. 2001. Image reconstruction using the wavelet transform for positron emission tomography.IEEETransactionsonMedicalImaging, 20(11): 1188-1193.

Chopra S, Marfurt K J. 2007. Seismic Attributes for Prospect Identification and Reservoir Characterization. Tulsa: Society of Exploration Geophysicists and European Association of Geoscientists and Engineers.Cohen I, Coifman R R. 2002. Local discontinuity measures for 3-D seismic data.Geophysics, 67(6): 1933-1945.

Darrozes J, Gaillot P, De Saint-Blanquat M, et al. 1997. Software for multi-scale image analysis: The normalized optimized anisotropic wavelet coefficient method.Computers&Geosciences, 23(8): 889-895. Gaillot P, Darrozes J, Bouchez J L. 1999. Wavelet transform: a future of rock fabric analysis?JournalofStructuralGeology, 21(11): 1615-1621.

Gaillot P, Darrozes J, de Saint Blanquat M, et al. 1997. The normalised optimised anisotropic wavelet coefficient (NOAWC) method: An image processing tool for multi-scale analysis of rock fabric.GeophysicalResearchLetters, 24(14): 1819-1822.

Gao Z S, Li Q X, Liu C, et al. 1998. Application of 2D wavelet transform in coherent noise estimation.JournalofChangchunUniversityofScienceandTechnology(in Chinese), 28(4): 433-437.

Gersztenkorn A, Marfurt K J. 1996. Eigenstructure based coherence computations.∥ 1996 SEG Annual Meeting. Denver, Colorado: SEG, 328-331. Gersztenkorn A, Marfurt K J. 1999. Eigenstructure-based coherence computations as an aid to 3-D structural and stratigraphic mapping.Geophysics, 64(5): 1468-1479.

Li C F, Liner C. 2008. Wavelet-based detection of singularities in acoustic impedances from surface seismic reflection data.Geophysics, 73(1): V1-V9.

Liu C, Wang D, Yang B J, et al. 2003. Apply of 2-D wavelet transforms in surface wave eliminations of seismic exploration.ProgressinGeophysics(in Chinese), 18(4): 711-714.

Li Y Y, Li Q C, Ding L B, et al. 2004. The iterative ground-roll suppressing method based on 2D wavelet transform . The 20th Annual Meeting of Chinese Geophysical Society, Xi′an, Chinese Geophysical Society.Lu W K, Li Y D, Zhang S W, et al. 2005. Higher-order-statistics and supertrace-based coherence-estimation algorithm.Geophysics, 70(3): P13-P18. Mallat S. 2008. A Wavelet Tour of Signal Processing: The Sparse Way. 3rd ed. New York: Academic Press.

Mallat S, Hwang W L. 1992. Singularity detection and processing with wavelets.IEEETransactionsonInformationTheory, 38(2): 617-643.

Marfurt K J, Kirlin R L. 2000. 3-D broad-band estimates of reflector dip and amplitude.Geophysics, 65(1): 304-320.

Marfurt K J, Kirlin R L, Farmer S L, et al. 1998. 3-D seismic attributes using a semblance-based coherency algorithm.Geophysics, 63(4): 1150-1165. Marfurt K J, Sudhaker V, Gersztenkorn A, et al. 1999. Coherency calculations in the presence of structural dip.Geophysics, 64(1): 104-111.

Murenzi R. 1990. Ondelettes multidimensionnelles et applications à l′analyse d′images[Ph. D]. Belgique: Louvain-la-Neuve, Université Catholique de Louvain.Murenzi R, Kaplan L, Antoine J P, et al. 1997. Computational complexity of the continuous wavelet transform in two dimensions.∥Proc. Army Research Laboratory Symp. on Sensors and Electron Devices, College Park. Murenzi R, Zhai W P, Namuduri K R, et al. 1999. Scale-angle CWT features: Application in object recognition.∥Proc. SPIE 3813, Wavelet Applications in Signal and Image Processing VII. Denver, CO, USA: SPIE, 16-27.

Ouillon G, Castaing C, Sornette D. 1996. Hierarchical geometry of faulting.JournalofGeophysicalResearch-SolidEarth, 101(B3): 5477-5487.

Ouillon G, Sornette D, Castaing C. 1995. Organisation of joints and faults from 1 cm to 100 km scales revealed by optimized anisotropic wavelet coefficient method and multifractal analysis.NonlinearProcessesinGeophysics, 2: 158-177.

Randen T, Monsen E, Signer C, et al. 2000. Three-dimensional texture attributes for seismic data analysis.∥SEG Technical Program Expanded Abstracts. SEG, 668-671.

Song J G, Wang H K, Mu X. 2005. Research and implementation of cohere algorithm based on DOA.ProgressinGeophysics(in Chinese), 25(5): 1662-1669.

Sun X P, Du S T. 2000. An application of power algorithm in the seismic attribute analysis.ComputingTechniquesforGeophysicalandGeochemicalExploration(in Chinese), 22(4): 316-321.

Wang W, Gao J H, Chen W C, et al. 2012. Random seismic noise suppression via structure-adaptive median filter.ChineseJ.Geophys. (in Chinese), 55(5): 1732-1741.

Wang X K, Gao J H, Chen W C, et al. 2012. An efficient implementation of eigenstructure-based coherence algorithm using recursion strategies and the power method.JournalofAppliedGeophysics, 82: 11-18.

Wang X W, Du M J, Liu J Y, et al. 2002a. Computation of seismic coherence cube based on wavelet transform and its application.GeophysicalProspectingforPetroleum(in Chinese), 41(3): 334-338.

Wang X W, Yang K Q, Zhou L H, et al. 2002b. Methods of calculating coherence cube on the basis of wavelet transform.ChineseJ.Geophys. (in Chinese), 45(6): 847-852.

Wu A D, Mou Y G. 1997. The use of 2D wavelet transform in imporving the resolution of seismic data.GeophysicalProspectingforPetroleum(in Chinese), 36(1): 36-42.Yang L Q, Song H B, Hao T Y, et al. 2005. Method of 2-D wavelet transform in attenuating random noise.GeophysicalProspectingforPetroleum(in Chinese), 44(1): 4-6.

Ye Z L, He J J. 2006. The improvement of the coherency cube technique and its application.ComputingTechniquesforGeophysicalandGeochemicalExploration(in Chinese), 28(2): 121-123. Yuan S J. 2007. A review of seismic coherence techniques.ProgressinExplorationGeophysics(in Chinese), 30(1): 7-15.

Zhang C J, Wang X D, Zhang H R, et al. 2006. Algorithm for image edge detection based on two-dimension continuous wavelet transform.Opto-ElectronicEngineering(in Chinese), 33(4): 93-96.Zhang H, Pan D M, Zhang X Y. 2007. Application of 2-D wavelet transformation in eliminating surface wave interference.GeophysicalProspectingforPetroleum(in Chinese), 46(2): 147-150.

Zhang J H, Dong M, Zhou Z X, et al. 2007. The coherence method based on GST and its application.NaturalGasIndustry(in Chinese) , 27(S1): 381-383.

Zhang J H, Wang Y Y, Zhao Y, et al. 2004. Detection of coherent data volume by wavelet multi resolution and application.OilGeophysicalProspecting(in Chinese), 39(1): 33-36.

Zhou Y H, Gao J H. 2007. Application of local structural entropy measure on the detection of discontinuity of seismic data.CoalGeology&Exploration(in Chinese), 35(1): 71-73.

Zhou Y H, Gao J H, Chen W C. 2007. Local structural entropy based on frequency-division instantaneous phase for enhancing seismic discontinuities.∥SEG Technical Program Expanded Abstracts. SEG, 846-850.

Zhang B, Lin T F, Guo S G, et al. 2016. Noise suppression of time-migrated gathers using prestack structure-oriented filtering.Interpretation, 4(2): SG19-SG29.

附中文參考文獻(xiàn)

陳文超, 師振盛, 汪文秉等. 2000. 小波變換在去除探地雷達(dá)信號(hào)直達(dá)波的應(yīng)用. 電波科學(xué)學(xué)報(bào), 15(3): 352-357.

陳文超, 汪文秉, 趙榮椿. 2003. 探地雷達(dá)小波的構(gòu)造及在提高雷達(dá)信號(hào)分辨率中的應(yīng)用. 煤田地質(zhì)與勘探, 31(5): 46-48.

高振山, 李勤學(xué), 劉財(cái)?shù)? 1998. 二維小波變換在消除相干噪聲中的應(yīng)用. 長春科技大學(xué)學(xué)報(bào), 28(4): 433-437.

李媛媛, 李慶春, 丁梁波等. 2004. 二維小波變換迭代法去除地震面波.∥中國地球物理學(xué)會(huì)第二十屆年會(huì). 西安: 中國地球物理學(xué)會(huì).

劉財(cái), 王典, 楊寶俊等. 2003. 二維小波變換在地震勘探面波消除中的應(yīng)用. 地球物理學(xué)進(jìn)展, 18(4): 711-714.

宋建國, 王海昆, 穆星. 2005. 基于DOA的相干算法研究與應(yīng)用. 地球物理學(xué)進(jìn)展, 25(5): 1662-1669.

孫夕平, 杜世通. 2000. 乘冪法在地震屬性分析中的應(yīng)用. 物探化探計(jì)算技術(shù), 22(4): 316-321.

王偉, 高靜懷, 陳文超等. 2012. 基于結(jié)構(gòu)自適應(yīng)中值濾波器的隨機(jī)噪聲衰減方法. 地球物理學(xué)報(bào), 55(5): 1732-1741.

王西文, 蘇明軍, 劉軍迎等. 2002a. 基于小波變換的地震相干體算法及其應(yīng)用. 石油物探, 41(3): 334-338.

王西文, 楊孔慶, 周立宏等. 2002b. 基于小波變換的地震相干體算法研究. 地球物理學(xué)報(bào), 45(6): 847-852.

吳愛弟, 牟永光. 1997. 二維小波變換在提高地震資料分辨率中的應(yīng)用. 石油物探, 36(1): 36-42.

楊立強(qiáng), 宋海斌, 郝天珧等. 2005. 基于二維小波變換的隨機(jī)噪聲壓制方法研究. 石油物探, 44(1): 4-6.

葉增爐, 何建軍. 2006. 地震相干體算法的改進(jìn)及應(yīng)用. 物探化探計(jì)算技術(shù), 28(2): 121-123.

苑書金. 2007. 地震相干體技術(shù)的研究綜述. 勘探地球物理進(jìn)展, 30(1): 7-15.

張長江, 汪曉東, 張浩然等. 2006. 圖像邊緣檢測的二維連續(xù)小波變換法. 光電工程, 33(4): 93-96.

張華, 潘冬明, 張興巖. 2007. 二維小波變換在去除面波干擾中的應(yīng)用. 石油物探, 46(2): 147-150.

張軍華, 王月英, 趙勇等. 2004. 小波多分辨率相干數(shù)據(jù)體的提取及應(yīng)用. 石油地球物理勘探, 39(1): 33-36.

張軍華, 董猛, 周振曉等. 2007. 基于GST的相干體方法研究及應(yīng)用. 天然氣工業(yè), 27(S1): 381-383.

周艷輝, 高靜懷. 2007. 局部結(jié)構(gòu)熵算法在地震數(shù)據(jù)不連續(xù)性檢測中的應(yīng)用. 煤田地質(zhì)與勘探, 35(1): 71-73.

(本文編輯 張正峰)

Detecting method of seismic discontinuities based on high dimensional continuous wavelet transform

WANG Xiao-Kai1,2,3, GAO Jing-Huai1,3*, CHEN Wen-Chao3, YANG Chang-Chun2, ZHU Zhen-Yu4

1Instituteofwaveandinformation,Xi′anJiaotongUniversity,Xi′an710049,China2InstituteofGeologyandGeophysics,TheChineseAcademyofSciences,Beijing100029,China3SchoolofElectronicandInformationEngineering,Xi′anJiaotongUniversity,Xi′an710049,China4ResearchCenterofCNOOC,Beijing100027,China

Seismic discontinuous structures, such as faults, play an important role in seismic exploration. We first introduce the basic theory of high-dimensional continuous wavelet transform, including the definition of high-dimensional continuous wavelet transform, slices of high-dimensional continuous wavelet transform coefficients and the fast realization of high-dimensional continuous wavelet transform. We use a synthetic signal example to illustrate that the high-dimensional continuous wavelet transform has good angular resolution compared with high-dimensional separable wavelet transform. Then we introduce high-dimensional continuous wavelet transform to seismic discontinuity detection for the first time and propose a seismic discontinuities detecting method based on small-scale high-dimensional continuous wavelet transform coefficients. The realization procedures of the proposed method are also presented. Finally, synthetic signal examples and field examples confirm the effeteness of the proposed method.

High dimensional CWT； Fast Fourier transform； Seismic discontinuities

10.6038/cjg20160922.

國家自然科學(xué)基金(41390454, 41504092, 41274125)、中國博士后基金(2016T90925, 2015M572566)、中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金及北京數(shù)學(xué)與信息交叉科學(xué)協(xié)同創(chuàng)新中心油氣勘探與開發(fā)平臺(tái)聯(lián)合資助.

王曉凱，男，漢族，博士，主要從事地震信號(hào)處理和地震儲(chǔ)層解釋方法研究.E-mail:xkwangxjtu@163.com

*通訊作者 高靜懷，男，教授，博士生導(dǎo)師，主要從事復(fù)雜介質(zhì)中地震波傳播及地震資料處理的理論與方法研究．E-mail:jhgao@mail.xjtu.edu.cn

10.6038/cjg20160922

P631

2014-08-14，2015-12-05收修定稿

王曉凱，高靜懷，陳文超等. 2016. 基于高維連續(xù)小波變換的地震資料不連續(xù)性檢測方法研究.地球物理學(xué)報(bào)，59(9):3394-3407,

Wang X K, Gao J H, Chen W C, et al. 2016. Detecting method of seismic discontinuities based on high dimensional continuous wavelet transform.ChineseJ.Geophys. (in Chinese),59(9):3394-3407,doi:10.6038/cjg20160922.