專利名稱：基于小波變換的地震相干體計算方法
技術領域：
本發明涉及油氣田勘探技術領域，屬于地震資料解釋范疇，具體地說是一種基于小波變換的地震相干體計算方法。
背景技術：
地震相干是指相鄰地震道之間地震屬性，如波形、振幅、頻率、相位等相似程度的測量。
一般情況下，現在所作的相千都是基于振幅的計算，利用多道相似性將三維振幅數據體經計算轉化為相關系數數據體，在顯示上強調不相關異常，突出不連續性。它的前提在于假設地層是連續的，地震波有變化、也是漸變的，因此相鄰道、線之間是相似的。當地層連續性遭到破壞而發生變化時，如斷層、尖滅、侵入、變形等，導致地震道之間的波形特征發生變化，進而導致局部道與道之間的相關性表現邊緣相似性的突變，地層邊界、特殊巖性體的不連續性會得到低相關值的輪廓。
通過計算地震相干數據體突出那些不相干的地震數據和提取三維相關屬性體，由縱向和橫向上局部的波形相似性可以得到三維地震相關性的估計值，從而把三維反射振幅數據體轉換成三維相似系數或相關值的數據體。
地震相干的最有效和最重要的應用之一是同相軸的檢測，特別是當同相軸的幅值較小并隱藏在噪聲中時，除檢測外，該測量還對同相軸的強度賦予一個定量的值，若其估值是一個能量歸一化的度量值，則很容易轉換為信噪比，因此相干可以對數據質量給出評估，它的計算可以在相干較弱或被噪聲干擾的情況下提供出數據相似性的定量值，互相關則是它的計算基礎。目前主要有如下幾種計算方法。
1) Cl相干算法
目前的應用軟件，如Landmark, GeoFrame等算法大部分都是以經典的歸一化互相關為基礎的計算(Bahorich和Farmer, 1995， 1996),稱之為第一代算法Cl。它利用主測線和聯絡測線方向的相關系數合成主聯方向的相關系數。Cl算法具有計算量小、易于實現的優點,但是，受資料的限制大。
設相鄰兩地震道為x(n)， y(n)時窗長度為K，那么二者在縱測線/延遲的互相關函數為<formula>formula see original document page 4</formula>假設P和g分別代表地層的傾角和方位角，那么對三維三道和多道情況分別用如下公式計算其相干值-<formula>formula see original document page 5</formula>則ci算法的相干值為<formula>formula see original document page 5</formula>
2) C2相干算法
該算法和C3算法引入了協方差矩陣，使其可對任意道數進行相似分析，估計其相干性。 C2相干算法除了在噪聲環境下更穩健地測量相干之外，垂直分析時窗能被限制在只有幾個時 間樣點范圍內，能夠精確做出薄而小的地層特征圖。
假設有J道地震記錄在所選的分析時窗內，坐標(x戶為)處的值為W;。以/ =""為中 心的一對視傾角(p， 9)的2M+1個采樣點，這2M+l個采樣點對應著一個JxJ的協方差矩 陣C:
<formula>formula see original document page 5</formula>這里=、.(附"-戸,.-肌.)，是地震道沿著f = mAf -戸;-肌.的內插值。那么C2相干 算法也可用公式描述如下<formula>formula see original document page 5</formula>
式中向量<formula>formula see original document page 5</formula>力表示協方差矩陣的跡。那么，C2相干算法得到的相
干值為<formula>formula see original document page 5</formula> 3) C3相干算法
C3相干算法借助C2相干算法中引入的協方差矩陣來實現,其分辨率較Cl、 C2算法更高， 同時不需要層位的約束。假設^(乂-l，2,…，J)是協方差矩陣C的第j個特征值，其中A,是最 大值。C3相干算法的實現如下那么C3相千算法得到的相干值為.-
C3 = max C3
地震解釋過程中最大的問題仍是斷層解釋，M.Bahorich等于1995年提出地震相干數據的應用方法，從而使得斷層的自動解釋成為可能。K丄Marfurt等于1998年提出了基于相似性的相干算法計算地震屬性，該方法穩定性強，對斷層的劃分精度很高。我國一些學者引入該項技術后，在解釋斷層上取得了經驗和成果，并在此基礎上提出了改進。
M.Bahorich等提出的方法是在地震振幅數據體上作相干分析，振幅縱橫強弱關系，信噪比都會直接影響計算出的相干體的質量，尤其是斷層兩側地震信號信噪比低，在計算出的相干體上斷點模糊，噪聲干擾很大，而且能量也較弱。K丄Marfurt等提出采用地震信號和地震信號Hilbert變換來計算相干體，增強穩定性，并且取得了明顯效果。但是，Hilbert變換的質量在計算相干體時非常重要，該變換對噪聲十分敏感，如果原始地震數據中噪聲太大，很可能計算出Hilbert變換的誤差很大，以致影響K丄Marfurt等的方法應用效果。
小波變換方法在20世紀80年代興起后，迅速用于地震資料處理領域。小波變換之所以有優于Fourier變換的特點在于它可以研究信號的局部特征，而Fourier變換是著重研究信號的整體特征。小波函數可以根據信號特征進行構造，高靜懷在Morlet小波基礎上構造出了適應地震資料處理的改進Morlet小波函數，王西文在改進的Morlet小波基礎上構造出了適應高分辨地震資料處理的導數小波函數，這就為小波變換的時間-尺度域分析、分離信號和噪聲以及分頻處理帶來極大的方便。小波變換計算出的Hilbert變換有較好的抗噪性，這為準確地提取瞬時特征參數奠定了基礎。
M.Bahorich等和K丄Marfurt等提出的方法都是在全頻帶上計算相干體，而未考慮到突出一定頻帶相干體的問題，尤其是突出用以解釋小斷層的高頻段相干體信息。
發明內容
本發明要解決的技術問題在于提供一種一種基于小波變換的地震相干體計算方
法。應用本發明，能提高時頻分析精度、提高計算出的相干體質量，突出被忽略的小斷層信息，提高地震解釋精度。
本發明主要包括如下步驟1)對現有技術獲得的地震資料進行去噪音、拓寬頻譜預處理，以提高地震資料的信噪比，拓寬地震資料的有效頻帶寬度；2) 對步驟1所處理過的地震資料進行精細頻譜分析，確定地震資料有效頻譜范圍；
3) 對地震資料進行層位解釋，并對拾取的地震解釋層位進行內插、平滑等處理；
4) 以地震數據、解釋的地震層位為輸入，利用高分辨導數小波函數，對輸入地震數據的某一個地震道做小波變換，獲得小波變換分頻結果；其中導數小波函數如式(1):
&、 嚴
式(1)中m為角頻率，C為常數，t為時間；
5) 利用式(2):
cg -i 。
計算瞬時小波域的瞬時振幅或/和瞬時頻率或/和瞬時相位；
式(2)中:Cg為常數;S(b,a)為地震道的小波變換；H[S(b)]為S(b)的Hilbert變換；
6) 利用基于小波瞬時相位的相干體算法或基于小波變換的相干體算法計算瞬時相干體；其中，基于小波瞬時相位的相干體算法如式(3):
「 max A (6， a, x;, )lmax (6， a, x,., )
式(3)中max^ (6，。，/^,，乂.)禾口 max/ y (6，<3，附^,，乂)分另U表示在延遲為/禾口附時，px和A,的最大值。a為尺度因子定義在頻帶["W,w]， ,M，2...N-1，頻率隨著序號Z的增大頻率增高；a地震道位置(x,.,")和(x,+;, y,O，地震道延遲/的互相關系數；a,是在地震道位置在(A，》)和(x,、 ；;, +/)，地震道延遲w的互相關系數；是上述沿縱測線方向即地震道延遲/和沿橫測線方向即地震道延遲"！的3D相關系數；
基于小波變換的相干體算法如式(4):
<formula>formula see original document page 7</formula>
式(4)中Afa為采樣時間間隔；aEa" aw；7)利用式(5):<formula>formula see original document page 8</formula>
重構瞬時相干體，獲得高頻或低頻相干數據體；式(5)中cf,是重構系數，d,E[O， l];
8) 對每一個地震道重復上述4一7步驟，得到所有地震道的相干計算結果；
9) 利用所求取的瞬時相干數據進行制圖，提供給地震資料解釋人員，用于斷層識別等研究。本發明采用模擬地震子波的小波函數分頻計算瞬時相位，提高了時頻分析精度；根據
M.Bahorich等方法用分頻的瞬時相位計算相千體，通過重構系數，利用基于小波瞬時相位的相干體算法對一定頻帶相干體放大或減小，來突出特定頻段相干體，易于突出被忽略的小斷層信息，提高地震解釋精度；直接用小波分頻變換的實、虛部，根據基于小波變換的相干體算法，克服了直接計算Hilbert變換帶來的誤差，提高了相干體計算質量。
本發明具體實現原理如下
1)理論基礎
地震信號S (t)的小波變換可以表示為
<formula>formula see original document page 8</formula>
式中6為時間因子，且t Z)GR; a為尺度因子，且"￡/ \{0},& (6， a)為小波變
換的實部，& "， a)為小波變換的虛部，&如果小波函數滿足
"、
是小波函數'
g (,) EL1 (i ， (1L2 (7 ，g (to) ELI (R\{0}， dco/a)) HL2{R\{0}， dco/oo}Cg邦
則
式中，Cg為常數；S(Z)，")為地震道的小波變換；/^(6)]為S(6)的Hilbert變換;可分別對應不同尺度因子a，即相當不同頻率下定義瞬時振幅，瞬時相位，瞬時頻率。 以及不同尺度因子a的阻尼瞬時頻率。
小波函數選取如下式所示改進的Morlet小波，該小波函數用常數c控制高斯函數來調制 小波函數g(t)的特征，可最大程度地模擬地震子波。
.,—w
式中，c為常數；m為角頻率。
為了更好的模擬高頻率地震信號，采用如下的高分辨導數小波函數，主要用于高頻段地 震信號的分解。
& \ 嚴
式中.* m為角頻率，c為常數，t為時間。 2)基于瞬時相位的相干體算法
據此如果在地震道位置(X,， 乂)和(X,+/, 乂)，地震道延遲/的互相關系數^，在地震
道位置在(x;，和Gc,, 乂+/)，地震道延遲w的互相關系數& ，那么可將上述沿縱測線 方向(地震道延遲/)和沿橫測線方向(地震道延遲m)的3D相關系數/^定義為
、/ max / x (6,", x,.，_y,.)|max / v (6, a， ；c;, )
其中maxpx (6,a，/^,乂)禾卩max/^ (6，a，m^,，乂)分別表示在延遲為/和w時，/^禾卩a,的 最大值。"為尺度因子定義在頻帶[",,",+小z'=l，2...N-l，頻率隨著序號/的增大頻率增高。 3)基于小波變換的相干體算法
相似算法是穩定的算法，其先定義一個以分析點為中心的J道橢圓或矩形分析時窗。若取 分析點坐標(x， y)為局部中心，則定義相似系數為。(b,a,p,q)。式中給出的相似估計，但 是對一些小的相干地震同相軸是不穩定的。因此，在一個高度為2w(v^KAb)ms的垂直時窗 上計算一個平均相似系數，該平均相似系數，即相干估計值定義為
乂=1
式中，a6為采樣時間間隔；flE[a,， ",+/]。 4)重構相干體
9l )所計算的相干體僅是在 一定頻帶下的相干體，因此可以重構相干體，其定義為
<formula>formula see original document page 10</formula>式中，d,是重構系數，cf,e[O， 1]。
根據基于小波變換的相干體算法(簡稱本發明相干算法2)及其Sh (6, a)和S, (jb, a) 的定義，用計算出的相干體，其重構相干體定義為
C* a,戶，《)=Z《C[6，", p,《]
/V-l
圖1 ( a)是實施例1基于小波瞬時相位的相干體算法計算得到相干切片，
圖1 (b )是實施例1基于小波變換的相干體算法計算得到相干切片，
圖2 (a )是實施例2館陶組底界高頻相干，
圖2 ( b)是實施例2館陶組底界高頻相干與斷層疊合
圖3 ( a)是實施例2館陶組底界低頻相干，
圖3 ( b)是實施例2館陶組底界低頻相干與斷層疊合，
圖4 ( a)是實施例2原館陶組底界構造圖，
圖4 (b)是實施例2新館陶組底界構造圖，
圖5 ( a )是實施例3周清莊油田的3D地震數據中的一條剖面，
圖5 ( b )是實施例3該地震剖面采用K丄Marfurt方法計算相干剖面，
圖5 ( c )是實施例3該地震剖面采用K丄Marfurt方法利用小波計算的相干剖面，
圖5 (d )是實施例3該地震剖面采用本發明方法計算的高頻瞬時相位剖面，
圖5 (e)是實施例3該地震剖面采用本發明方法計算出低頻瞬時相位剖面，
圖5 ( f)是實施例3該地震剖面采用本發明方法重構瞬時相位剖面。
具體實施例方式
實施例l
本發明包括如下步驟1 、對利用現有技術，通過野外高分辨率地震采集設備所獲得的原始地震數據進行去噪音、 拓寬頻譜預處理，以提高地震資料的信噪比，拓寬地震資料的有效頻帶寬度；
2、 對步驟1所處理過的地震資料進行精細頻譜分析，確定地震資料有效頻譜范圍；
3、 對地震資料進行層位解釋，并對拾取的地震解釋層位進行內插、平滑等處理；
4、 以地震數據、解釋的地震層位為輸入，利用高分辨導數小波函數，對輸入地震數據的 某一個地震道做小波變換，獲得小波變換分頻結果；其中導數小波函數如式(1):
5g 廣 — ——="m _ c / )e e
&\ 嚴
式(1)中m為角頻率，c為常數,
S、由式(2):
為時間；
丄
f雄,《4)+z卿)〗
計算瞬時小波域的瞬時振幅、瞬時頻率或瞬時相位；
式(2)中Cg為常數;S(b,a)地震道的小波變換；H[S(b)]為S(b)的Hilbert變換
6、利用基于小波瞬時相位的相干體算法或基于小波變換的相干體算法計算瞬時相干體; 其中，基于小波瞬時相位的相干體算法如式(3):
= 、( max px (6， a， x, ， 乂 )Jmax (6， a， x,. ， )
式(3)中maxpx (jb,a,/,Aj/,)禾卩maxjcv (&a,nv^)/,)分別表示在延遲為/和m時，px 和Py的最大值；a為尺度因子定義在頻帶[^a,w]， /'=1,2...N-1,頻率隨著序號,'的增大頻率增 高；px是在地震道位置(jc,， 乂)和(x,+/，》 )，地震道延遲/的互相關系數；Py是在地震 道位置(義,.，y,)和(x,.， "+;)，地震道延遲m的互相關系數；Pxy是上述沿縱測線方向和沿 橫測線方向的3D相關系數；即P xy是地震道延遲I和地震道延遲m的3D相關系數；
基于小波變換的相干體算法如式(4):
<formula>formula see original document page 11</formula>
式(4)中Ab為采樣時間間隔 7、利用式(5):重構造瞬時相干體，獲得高頻或低頻相干數據體；式(5)中d,是重構系數，"e[O， l];
8、 對每一個地震道重復上述4一7步驟，得到所有地震道的相干計算結果；
9、 利用所求取的瞬時相干數據進行制圖，提供給地震資料解釋人員，用于斷層識別等研 究。
圖1 (a )為沿層采用本發明相干算法1即基于小波瞬時相位的相干體算法所獲得的高 頻相干切片；圖l (b)為沿層采用本發明相干算法2即基于小波變換的相干體算法所獲得 的高頻相干切片。由圖可以得出如下結論，該發明計算相干切片具有識別斷層或不均質體分 辨率高的特點。
實施例2 圖2 (a) _圖4 (b)為本發明在大港油田唐東一馬南地區的應用情況。其 步驟與實施例l相同。
圖2 (a )是實施例2館陶組底界高頻相千，該結果可以很好的識別微小斷裂；圖中非 常清楚地顯示出主要控制油田構造斷裂有兩組， 一組以F1為代表的北北東向斷裂； 一組是 以F2為代表的西西東向斷裂。圖中，A處是一組4-5條斷層，北北東向小斷層交匯在F2斷 裂上，而且斷裂展布清晰，分辯率高；B處是近東西向斷裂交匯在F2斷裂上，交匯點清晰； C處是北北東向斷裂交匯在F2斷裂上，斷裂點清晰；D處是一組3條雁行式排列北北東向 斷裂，分辯很清楚；E處是一條北北東向斷裂與一條北東向斷裂交匯點。
圖2 (b)實施例2館陶組底界高頻相干與構造圖的疊合，疊合結果說明該相干結果與 解釋結果一致，符合地質認識，而且左側三維數據解釋斷裂與右圖相干體切片斷裂線非常吻 合。尤其是A， B， C， D， E處斷裂交匯點清晰。
圖3 (a)實施例2館陶組底界低頻相千，該結果可以很好的識別較大斷裂，圖中Fl 和F2斷裂非常清楚。但是，A處北北東向這組小斷裂較高頻相干體分辯不清，B， C處斷裂 交匯處與高頻相干切片相比模糊不清，D處這組雁行排列斷裂間距較大，低頻相干體仍能分 辨清楚，E處斷裂交匯點模糊不清。
圖3 (b)實施例2館陶組底界低頻相干與構造圖的疊合，疊合結果說明該相干結果與 解釋結果一致，大斷裂格局非常清楚符合地質認識。
圖4 (a)和圖4 (b)分別是原沙三I底構造圖和是原沙三IV底構造圖。兩圖相比，斷裂特征相同，說明斷裂是繼承性的，這點與實際不符。圖中出現近似直角弧形斷裂，在塑 性地層中不可能出現這種斷裂，說明斷裂解釋有誤。圖中很難分析出構造受力方向和主要控 制構造斷裂體系，說明斷裂組合有問題。因此，這種構造圖無法為精細儲層預測研究提供正 確的構造信息；相反，還會引到出錯誤結論。
實施例3 圖5 (a)—圖5 (f)為本發明在大港油田周清莊應用的結果。其步驟與實 施例1相同。圖5 ( a )為實施例3周清莊油田的3D地震數據中的一條剖面，圖5 ( a ) 中所示的S31是一層生物灰巖儲層，S31至S34是一套河流相的砂泥巖互層沉積，其中砂巖 也是非常好的油氣儲層。二套油氣儲層受小斷裂控制。該油田已處于開發階段的中后期，研 究受大斷裂控制的小斷裂對油氣儲層控制作用十分重要。研究區斷裂發育，主要發育一組控 制箕狀斷陷正斷層，從地震剖面信息上很難解釋出圖中所示這組發育的正斷層。圖中所示斷 裂是經過精細構造研究后的解釋成果。
圖5 ( b )實施例3該地震剖面采用K丄Marfurt等的方法用Hilbert變換計算相干剖面， 由于地震信號信噪比較低，這使Hilbert變換誤差增大，以致在計算的相干剖面中斷層斷點非 常模糊，難以解釋出圖中所示這組發育的正斷層。
圖5 ( c )為實施例3該地震剖面采用K丄Marfort等的方法利用小波計算的相干剖面， 小波變換方法得到Hilbert變換，再計算相干剖面即基于小波變換的相干體算法。很明顯，圖 5(c)中斷層斷點較明顯。但是，該圖在用于小斷層解釋時，分辨率太低。
圖5 (d)為實施例3該地震剖面采用本發明方法計算的高頻瞬時相位剖面，首先利用 模擬地震子波的高分辨導數小波函數計算出高頻瞬時相位。很明顯，圖5 (d )中所示斷層 斷點清晰，尤其是在CDP1-CDP46;時間2100ms-2500ms之間小斷層斷點非常清楚。這為精 細構造解釋和小斷層展布的研究提供了基礎資料。
圖5 ( e)為實施例3該地震剖面采用本發明方法計算出低頻瞬時相位剖面，首先利用 本發明方法模擬地震子波的小波函數計算出低頻瞬時相位。很明顯，圖5 (e)中所示在 CDP1-CDP46;時間2100ms-2900ms之間，控制箕狀凹陷大斷層斷點非常清楚。
圖5 (f )為實施例3該地震剖面采用本發明方法重構瞬時相位剖面，取高頻重構系數 d產0.7;低頻重構系數df0.3得到重構相干剖面。圖5 (f )中所示高頻信息增強， 一些小 斷層顯示更明顯，很容易解釋出圖5 (f)中所示這組發育的正斷層。但是，在高頻信息放 大的同時，相應的高頻噪聲也變大。
總之，本發明給出基于小波變換的二種地震數據的相干體算法。
基于小波瞬時相位的相干體算法是用模擬地震子波的小波函數或高分辨導數小波函數分頻計算瞬時相位，根據分頻瞬時相位計算相干體，再進行重構。該算法可以提高斷層解釋 的分辨率，其應用效果是明顯的。
基于小波變換的相干體算法是用小波變換得到的實、虛部即Hilbert變換，直接計算相干 體。該算法可以克服直接計算Hilbert變換帶來的誤差，提高計算出的相干體質量。基于小波 瞬時相位的相干體算法可以得到高、低頻瞬時相位信息，有利于精細構造解釋，而且計算 出的相干體在小斷層分辨程度方面也較基于小波變換的相干體算法效果明顯。根據本發明的 研究實例，采用分頻重構的相千體，易于突出被忽略的小斷層信息。實際應用也證明基于小 波瞬時相位的相干體算法是十分有效的。
權利要求
1、一種基于小波變換的地震相干體計算方法，包括如下步驟1)對現有技術獲得的地震資料進行去噪音、拓寬頻譜預處理，以提高地震資料的信噪比，拓寬地震資料的有效頻帶寬度；2)對步驟1所處理過的地震資料進行精細頻譜分析，確定地震資料有效頻譜范圍；3)對地震資料進行層位解釋，并對拾取的地震解釋層位進行內插、平滑等處理；4)以地震數據、解釋的地震層位為輸入，利用高分辨導數小波函數，對輸入地震數據的某一個地震道做小波變換，獲得小波變換分頻結果；其中導數小波函數如式(1)式(1)中m為角頻率，c為常數，t為時間；5)利用式(2)計算瞬時小波域的瞬時振幅或/和瞬時頻率或/和瞬時相位；式(2)中Cg為常數；S(b，a)地震到的小波變換；H[S(b)]為S(b)的Hilbert變換；6)利用基于小波瞬時相位的相干體算法或基于小波變換的相干體算法計算瞬時相干體；其中，基于小波瞬時相位的相干體算法如式(3)式(3)中max ρx(b，a，l，xi，yi)和max ρy(b，a，m，xi，yi)分別表示在延遲為l和m時，ρx和ρy的最大值；a為尺度因子定義在頻帶[ai，ai+1]，i＝1，2...N-1，頻率隨著序號i的增大頻率增高；ρx是在地震道位置(xi，yi)和(xi+1，yi)，地震道延遲l的互相關系數；ρy是在地震道位置(xi，yi)和(xi，yi+1)，地震道延遲m的互相關系數；ρxy是上述沿縱測線方向和沿橫測線方向的3D相關系數；即ρxy是地震道延遲1和地震道延遲m的3D相關系數；基于小波變換的相干體算法如式(4)式(4)中Δb為采樣時間間隔；a∈[ai，ai+1]；7)利用式(5)重構瞬時相干體，獲得高頻或低頻相干數據體；式(5)中di是重構系數，di∈
；8)對每一個地震道重復上述4—7步驟，得到所有地震道的相干計算結果；9)利用所求取的瞬時相干數據進行制圖，提供給地震資料解釋人員，用于斷層識別等研究。
全文摘要
基于小波變換的地震相干體計算方法，步驟如下1)對地震資料進行預處理，2)對處理過的地震資料進行精細頻譜分析，確定地震資料有效頻譜范圍；3)對地震資料進行層位解釋，并進行內插、平滑等處理；4)以地震數據、解釋的地震層位為輸入，利用高分辨導數小波函數，對輸入地震數據的某一個地震道做小波變換，獲得小波變換分頻結果；5)計算瞬時小波域的瞬時振幅或/和瞬時頻率或/和瞬時相位；6)利用基于小波瞬時相位的相干體算法或基于小波變換的相干體算法計算瞬時相干體；7)重構瞬時相干體，獲得高頻或低頻相干數據體；8)對每一個地震道重復上述4-7步驟，得到所有地震道的相干計算結果；9)利用所求取的瞬時相干數據進行制圖。
文檔編號G01V1/28GK101545984SQ20091013838
公開日2009年9月30日 申請日期2009年5月5日 優先權日2009年5月5日
發明者劉軍迎, 楊午陽, 王西文 申請人:中國石油集團西北地質研究所