黎康毅，陳學(xué)華*，，吳昊杰，呂丙南，趙晨斐

（1.成都理工大學(xué)油氣藏地質(zhì)及開發(fā)工程國家重點(diǎn)實(shí)驗(yàn)室，四川成都 610059；2.成都理工大學(xué)地球探測與信息技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室，四川成都 610059）

0 引言

斷層和裂縫的連通網(wǎng)絡(luò)影響地下地層的形態(tài)及地震資料處理、解釋結(jié)果，要獲得精確的地質(zhì)模型，特別是在逆掩斷層及相關(guān)褶皺帶，需要精確檢測、解釋斷層與裂縫。

存在復(fù)雜斷層地區(qū)的三維地震數(shù)據(jù)質(zhì)量受采集、處理等環(huán)節(jié)的多種因素影響。在地下地質(zhì)結(jié)構(gòu)較復(fù)雜地區(qū)的地震資料往往信噪比較低，需要壓制隨機(jī)噪聲、提高地震數(shù)據(jù)的質(zhì)量，以有效實(shí)現(xiàn)可視化和屬性提取。通過應(yīng)用不同的濾波器、數(shù)據(jù)轉(zhuǎn)換和屬性組合技術(shù)，可以增強(qiáng)斷層與裂縫的地震響應(yīng)特征。Marr[1]、Witkin[2]、Koenderink[3]闡述了多尺度圖像分析的重要性，并在數(shù)字圖像處理中引入尺度空間理論；伍鵬等[4]將二維高斯迭代平滑濾波用于曲率屬性分析，但是沒有考慮裂縫方向，不能突顯不同方向的異常信息；在前人研究[5-7]的基礎(chǔ)上，Cho 等[8]、Park 等[9]應(yīng)用自適應(yīng)Gabor 濾波強(qiáng)化手指靜脈圖像并取得了不錯(cuò)的效果；Chen 等[10-12]提出了一種多尺度曲率計(jì)算方法及地震定向?yàn)V波新算法，為結(jié)合裂縫方向選取最佳濾波器提供了理論基礎(chǔ)；周元茂等[13-14]結(jié)合地震體曲率與各向異性高斯濾波方法，形成了組合型方向體曲率分析方法，可以突顯不同方向的異常信息，但是缺少自適應(yīng)性，需要人為提取裂縫方向；徐赫等[15]將旋轉(zhuǎn)窗加入希爾伯特變換，以提取地震資料中不同方向的有效信息；Jabar等[16]引入一種綜合屬性和圖像的分析技術(shù)增強(qiáng)地震斷層圖像。

本文結(jié)合前人方法，提出地震不連續(xù)信息的自適應(yīng)方向增強(qiáng)方法，根據(jù)裂縫方向選取高斯濾波器方向?qū)Φ卣饠?shù)據(jù)二次迭代濾波，可有效壓制干擾、突出地震數(shù)據(jù)的有效信息，從而增強(qiáng)特定方向的不連續(xù)信息并提取高精度地質(zhì)異常信息。

1 方法原理

本文提出的地震不連續(xù)信息的自適應(yīng)方向增強(qiáng)方法，根據(jù)裂縫走向不斷地更換各向異性高斯濾波器的方向，使數(shù)據(jù)所受的異常干擾小，同時(shí)保留更多有效信息。具體操作分為五個(gè)步驟（圖1）。

圖1 不連續(xù)信息增強(qiáng)流程

（1）在0°～180°范圍的所有給定方向的高斯窗掃描輸入地震數(shù)據(jù)，在每個(gè)位置分別選取窗內(nèi)振幅之和的最大值方向的處理結(jié)果；

（2）選定時(shí)窗掃描處理數(shù)據(jù)，將處理數(shù)據(jù)映射為L級(jí)灰度構(gòu)成地震紋理基元體；

（3）求取時(shí)窗內(nèi)數(shù)據(jù)特定方向的共生矩陣紋理參數(shù)，并判斷處理數(shù)據(jù)的特定方向；

（4）將判斷的裂縫方向代入各向異性高斯濾波方向矩陣進(jìn)行二次迭代；

（5）對(duì)處理數(shù)據(jù)采用各向異性體曲率屬性分析方法得到增強(qiáng)的最終圖像。

1.1 灰度共生矩陣

1973年，Haralick 等[17]提出了灰度共生矩陣模型(GLCM)，隨后被用于揭示圖像的紋理特征[18-22]。

其方法流程為，首先確定時(shí)窗寬度，掃描地震數(shù)據(jù)并進(jìn)行灰度處理而構(gòu)建灰度共生矩陣，設(shè)最大灰度級(jí)為L，假設(shè)水平方向(Crossline 方向)的測線數(shù)目為xmax，則地震數(shù)據(jù)的水平空間域?yàn)镹x={1,2,…,xmax}；假設(shè)垂直方向(Inline 方向)的采樣點(diǎn)數(shù)為ymax，則地震數(shù)據(jù)的垂直空間域?yàn)镹y={1,2,…,ymax}。設(shè)數(shù)據(jù)點(diǎn)對(duì)的距離為l，數(shù)據(jù)點(diǎn)對(duì)間的方向角為θ，則在(m,n)處灰度級(jí)g(m,n)=i的像素點(diǎn)與相距為l、方向?yàn)棣鹊?p,q)處灰度級(jí)g(p,q)=j的像素點(diǎn)，在圖像中出現(xiàn)的概率為P(i,j,l,θ)(圖2)，生成的灰度共生矩陣為

圖2 數(shù)據(jù)點(diǎn)對(duì)之間的關(guān)系

式中L為最大灰度級(jí)。由于l、θ的選擇不同，會(huì)導(dǎo)致不同的統(tǒng)計(jì)結(jié)果，如當(dāng)θ為0°、90°時(shí)，共生矩陣元素為

選取最常用的能量E、相關(guān)性C、熵H及對(duì)比度I等4個(gè)紋理參數(shù)判斷裂縫及突出蜿蜒通道的方向，即

其中

通過時(shí)窗掃描地震數(shù)據(jù)，每掃描一次便計(jì)算一次灰度共生矩陣，并得到共生矩陣的4個(gè)紋理參數(shù)，通過

判斷時(shí)窗位置的方向。如果滿足式(6)或式(7)，則輸出時(shí)窗的數(shù)據(jù)方向；如果式(6)和式(7)所得時(shí)窗方向不同或E、C、H、I全為0，則說明時(shí)窗數(shù)據(jù)沒有明顯的特征方向，跳過并繼續(xù)掃描數(shù)據(jù)。

圖3 為一個(gè)正十六邊形模型，其中每一邊的內(nèi)角均為22.5°。使用20 個(gè)采樣點(diǎn)的時(shí)窗掃描圖3，分別計(jì)算θ為0°、22.5°、45°、67.5°、90°、112.5°、135°、157.5°的共生矩陣，通過共生矩陣得到各方向的紋理參數(shù)。將不同方向的4 個(gè)紋理參數(shù)代入式(6)、式(7)，選取式(6)、式(7)對(duì)應(yīng)的角度作為該處紋理的方向。

圖3 正十六邊形模型

圖4 為判斷紋理方向的模擬測試結(jié)果。由圖可見，準(zhǔn)確判斷了正十六邊形的各個(gè)方向，表明了方法的可行性。

圖4 判斷紋理方向的模擬測試結(jié)果

1.2 各向異性高斯濾波器

傳統(tǒng)高斯函數(shù)濾波的基本原理是將高斯核函數(shù)與原始信號(hào)卷積，以二維高斯濾波為例，濾波器表達(dá)式為

式中：x和y表示二維高斯函數(shù)的兩個(gè)方向，當(dāng)x方向和y方向呈各向同性時(shí)，高斯濾波即為各向同性濾波(圖5a)；σ為高斯函數(shù)的方差。

圖5 二維高斯濾波器

將傳統(tǒng)高斯濾波器中的x和y與方向矩陣相乘，并在在x、y方向選取不同的尺度，就得到旋轉(zhuǎn)變換后的各向異性高斯濾波器，其濾波器表達(dá)式為

式中σx和σy為高斯函數(shù)的方差，分別決定高斯核寬度和長度，與平滑程度有直接關(guān)系，σx越大，x方向高斯窗口越寬，σy越大，y方向高斯窗口越寬，平滑程度越好。當(dāng)σx≠σy時(shí)，高斯濾波即為各向異性濾波(圖5b)。

為了利用各向異性高斯濾波描述地震數(shù)據(jù)中的蜿蜒河道、裂縫和其他地質(zhì)不連續(xù)點(diǎn)，需要考慮σx、σy和θ的變化（圖6）。文中提出了一個(gè)基于各向異性高斯濾波方法的工作流程突出地震數(shù)據(jù)中的蜿蜒通道和其他地質(zhì)不連續(xù)點(diǎn)。

圖6 各向異性高斯濾波器

根據(jù)實(shí)際地震資料中的裂縫尺度，給定各向異性高斯濾波器G(x,y,σx,σy,θ)的孔徑和縱橫比。然后用一系列不同方向的濾波器處理以分析點(diǎn)為中心的地震數(shù)據(jù)

式中：Gθ(x,y,σx,σy,θk)可以視為空間分析窗口；S(τ,xi,yj)是以(xi,yj)為位置點(diǎn)、以時(shí)窗位置τ為中心的地震數(shù)據(jù)；F(τ,xi,yj,θk)為濾波后地震數(shù)據(jù)，k為橢圓濾波器的方向。

對(duì)于0o～180o范圍內(nèi)所有給定的濾波方向θk，在每個(gè)空間位置取F(τ,xi,yj,θk)的最大值，即

其中F(τ,xi,yj)突出了真實(shí)方向的主要特征，而在其他方向抑制了噪聲。

然后，定義地震方向?yàn)V波算子，得到時(shí)窗位置τ的地震數(shù)據(jù)

式中：w為時(shí)窗窗口長度的一半；M為時(shí)窗尺度。

1.3 三維各向異性體曲率分析

在三維地震數(shù)據(jù)中，一般需要將三維地震數(shù)據(jù)體先轉(zhuǎn)化為傾角數(shù)據(jù)體，再根據(jù)傾角數(shù)據(jù)體計(jì)算任意點(diǎn)的曲率?？梢愿鶕?jù)復(fù)數(shù)道分析中的瞬時(shí)頻率和瞬時(shí)波數(shù)計(jì)算三維地震數(shù)據(jù)體的視傾角

式中：ω為瞬時(shí)角頻率；kx與ky分別為x與y方向的瞬時(shí)波數(shù)；px與qy分別為x與y方向的視傾角分量。

在此基礎(chǔ)上，構(gòu)造新的方位視傾角分量

式中β為根據(jù)地質(zhì)構(gòu)造方向而選取的角度。

根據(jù)最小二乘逼近原理，得到二次曲面方程

其中

式中f為常數(shù)。最大正曲率Kpos和最小負(fù)曲率Kneg是解釋實(shí)際構(gòu)造的主要地震屬性，其表達(dá)式為

2 實(shí)際資料處理實(shí)例

為了驗(yàn)證本文方法在LH 地區(qū)的應(yīng)用效果，分別利用一次方向各向異性高斯濾波和本文方法對(duì)三維疊后地震數(shù)據(jù)體進(jìn)行0°、90°方向的濾波。圖7 為LH地區(qū)地震振幅沿層切片。由圖可見，在多個(gè)方向存在強(qiáng)噪聲，從而難以識(shí)別有效的不連續(xù)地質(zhì)信息，也難以辨識(shí)各種尺度裂縫的空間走向。

圖7 LH 地區(qū)地震振幅沿層切片

圖8 為一次方向的各向異性高斯濾波與本文方法得到的最大正曲率屬性。由圖可見：①90°一次方向的各向異性高斯濾波(圖8a 的紅色、綠色方框處)及本文方法(圖8b 的紅色、綠色方框處)明顯增強(qiáng)了縱向不連續(xù)性信息及構(gòu)造異常特征，且圖8b 更明顯；0°一次方向的各向異性高斯濾波(圖8c 的紅色、黃色方框處)及本文方法(圖8d的紅色、黃色方框處)則著重突出了橫向構(gòu)造特征，明顯增強(qiáng)了斷層及裂隙連續(xù)性，且圖8d更明顯。②與一次方向各向異性高斯濾波(圖8a綠色箭頭處、圖8c 紅色箭頭處)相比，本文方法(圖8b綠色箭頭處、圖8d 紅色箭頭處)明顯突出了指定方向的地質(zhì)異常構(gòu)造。

圖8 一次方向的各向異性高斯濾波與本文方法得到的最大正曲率屬性

為了突出裂縫及斷層走向，使用灰度共生矩陣對(duì)一次方向的各向異性高斯濾波結(jié)果提取E、C、H、I等四個(gè)紋理參數(shù)（圖9）?？梢?，LH 地區(qū)的斷層及裂縫走向主要為90o、0o，小尺度裂縫及斷層的走向多變。圖10 為不同方法得到的實(shí)際地震數(shù)據(jù)的最大正曲率屬性。由圖可見：①原始最大正曲率屬性的隨機(jī)噪聲干擾嚴(yán)重，構(gòu)造、斷裂展布不清晰（圖10a）。②由一次方向的各向異性高斯濾波結(jié)果可大致看出斷裂、褶皺、落水洞等的輪廓，但是細(xì)節(jié)不明顯（圖10b 的紅色方框及綠色箭頭處）。③本文方法處理結(jié)果（圖10c）明顯突出了斷層及裂縫特征，并刻畫了某些細(xì)小的地質(zhì)異常構(gòu)造，如一些細(xì)小的褶皺、落水洞、隱蔽裂縫等（圖10c 的紅色方框及綠色箭頭處）。

圖9 地震數(shù)據(jù)的裂縫方向圖

圖10 不同方法得到的實(shí)際地震數(shù)據(jù)的最大正曲率屬性

3 結(jié)束語

本文提出了一種地震不連續(xù)信息的自適應(yīng)方向增強(qiáng)方法，在常規(guī)高斯濾波算法的基礎(chǔ)上，推導(dǎo)了各向異性方向迭代高斯平滑濾波公式。該方法對(duì)地震數(shù)據(jù)進(jìn)行兩次時(shí)窗掃描，自適應(yīng)地判斷裂縫等異常信息的方向，再根據(jù)判斷的方向選取匹配該方向的最優(yōu)濾波器進(jìn)行二次迭代處理，有效地壓制與裂縫等方向不同的干擾信息，極大突出了地質(zhì)構(gòu)造特征，更好地反映了特定方向的構(gòu)造細(xì)節(jié)和不連續(xù)性。同時(shí)該方法也可自行選取一個(gè)或多個(gè)特定方向，只展現(xiàn)所選方向的構(gòu)造信息。該方法可為構(gòu)造解釋、了解斷裂展布規(guī)律、儲(chǔ)層預(yù)測等提供技術(shù)支持，體現(xiàn)了算法的可行性和優(yōu)越性。