任意旋轉(zhuǎn)加窗希爾伯特變換的地震資料體邊緣檢測(cè)方法

徐 赫 陳學(xué)華* 呂丙南 黎康毅 徐 斌

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

0 引言

圖像分析中的邊緣檢測(cè)技術(shù)被廣泛用于地震數(shù)據(jù)裂縫檢測(cè)。目前地震數(shù)據(jù)邊緣檢測(cè)方法主要有基于圖像空間域微分算子(Prewitt算子、 Robert算子、Canny算子、Sobel算子、高斯偏導(dǎo)濾波器(LOG)等)的邊緣檢測(cè)法及以變換域小波分析為基礎(chǔ)的多尺度邊緣檢測(cè)法[1]。實(shí)際勘探經(jīng)驗(yàn)表明，油氣藏中普遍分布孔隙和裂縫，因此可利用地震資料檢測(cè)邊緣特征。Wu等[2-4]、丁燕等[5]利用深度學(xué)習(xí)算法提取三維地震圖像中的斷層及其產(chǎn)狀信息。

希爾伯特變換廣泛應(yīng)用于檢測(cè)地震資料不連續(xù)邊緣信息。黨志敏等[6]、熊曉軍等[7]利用廣義希爾伯特變換檢測(cè)含噪地震資料的邊緣信息。陳學(xué)華等[9-10]將高階偽希爾伯特變換用于地震資料邊緣檢測(cè)。謝靜等[11]提出了時(shí)間域加窗希爾伯特變換邊緣檢測(cè)方法?；谌S地震資料在深度方向的多尺度特征，李斌等[12]提出了多尺度加窗希爾伯特變換的地震資料體邊緣檢測(cè)方法。在前人的研究[13-17]基礎(chǔ)上，Lyu等[18-19]實(shí)現(xiàn)了基于二維希爾伯特變換的三維地震資料體邊緣檢測(cè)方法。利用上述方法檢測(cè)裂縫時(shí)，往往沒(méi)有考慮地震資料不連續(xù)信息在方向上的多樣化與復(fù)雜性。為此，周元茂等[20]聯(lián)合應(yīng)用各向異性高斯迭代平滑濾波方法與地震體曲率分析方法，形成了組合型方向體曲率分析方法。

基于希爾伯特變換的邊緣檢測(cè)算法用于實(shí)際地震資料時(shí)，提取的往往是裂縫綜合信息，信息繁雜且會(huì)弱化部分有效信息，難以分辨某一特定方向的數(shù)據(jù)特點(diǎn)，不利于分析裂縫走向。因此，考慮到斷層、裂縫等不連續(xù)信息的復(fù)雜性，為了分析地震資料中各個(gè)方向的不連續(xù)特征，本文構(gòu)建了一種任意角度旋轉(zhuǎn)的加窗希爾伯特變換(arbitrarily rotated windowed Hilbert transform，ARWHT)邊緣檢測(cè)算子，將其用于三維地震資料體邊緣檢測(cè)，不僅能提取地震資料在任意方向的不連續(xù)特征，還可突出局部異常信息，減少噪聲的影響，從而獲得更精確的儲(chǔ)層裂縫走向及分布信息。

1 方法原理

1.1 ARWHT邊緣檢測(cè)算子

為了充分提取、分析三維地震資料中任意方向的不連續(xù)信息并突出局部化特征，本文構(gòu)建了ARWHT邊緣檢測(cè)算子。

若將希爾伯特算子所有點(diǎn)表示為一個(gè)幾何信息矩陣Q， 將其與方向檢測(cè)模板T相乘﹐就得到旋轉(zhuǎn)變換后的幾何信息矩陣Q′，即

Q′=QT

(1)

其中

設(shè)加窗后的希爾伯特算子序列的樣點(diǎn)坐標(biāo)為P[n，h′(n)]，將其逆時(shí)針旋轉(zhuǎn)角度θ后，新算子序列的樣點(diǎn)坐標(biāo)為P′[nR，h′R(n)]，則有

(2)

其中

h′(n)=h(n)·w(n)

(3)

式中：h(n)為原始希爾伯特變換算子，n為采樣點(diǎn)號(hào)；w(n)為對(duì)算子作局部化處理的窗函數(shù)；a、b、c、d為設(shè)置的方向檢測(cè)模板初始形式的組成參數(shù)。將式(2)改為矩陣形式

(4)

設(shè)點(diǎn)P與坐標(biāo)原點(diǎn)O的連線OP與x軸的夾角為α，則有

(5)

(6)

(7)

為強(qiáng)調(diào)地震資料的局部化特征、突出不連續(xù)信息，本文選擇高斯窗對(duì)希爾伯特變換算子進(jìn)行局部化處理，使檢測(cè)對(duì)象的有效邊緣信息更突出、非邊緣信息的差異更明顯。高斯窗函數(shù)的表達(dá)式為

(8)

式中σ為決定窗函數(shù)時(shí)間寬度的尺度因子，即σ越大，窗函數(shù)越寬。圖1為ARWHT算子。

將ARWHT算子用于邊緣檢測(cè)，得到一個(gè)新的基于ARWHT的邊緣提取公式，邊緣提取結(jié)果為

(9)

其中

nR=ncosθ-h′(n)sinθ

(10)

h′R(n)=nsinθ+h′(n)cosθ

(11)

式中N為x(nR)的長(zhǎng)度。

1.2 基于ARWHT的地震資料體邊緣檢測(cè)算法的實(shí)現(xiàn)

由式(9)～式(11)計(jì)算任意特定方向的邊緣提取結(jié)果，進(jìn)而由

(12)

圖1 ARWHT算子(a)不旋轉(zhuǎn)； (b)逆時(shí)針旋轉(zhuǎn)45°； (b)逆時(shí)針旋轉(zhuǎn)90°； (d)逆時(shí)針旋轉(zhuǎn)135°

對(duì)于所有的算子旋轉(zhuǎn)角度θ∈[0°，180°]，取每個(gè)空間位置的D(m0，k，θ)最大值，由

(13)

得到融合了所有方向的基于ARWHT的三維地震資料體邊緣檢測(cè)結(jié)果。

2 模型試算與分析

為了驗(yàn)證本文算法的效果，設(shè)置一個(gè)正八邊形模型(圖2a)，利用ARWHT算子提取各方向的邊緣信息(圖2b、圖2c)。結(jié)果表明：逆時(shí)針旋轉(zhuǎn)45°的ARWHT算子檢測(cè)結(jié)果(圖2b)著重突出了0°、90°、135°的邊緣信息(135°邊緣信息絕對(duì)值最大)，有效過(guò)濾了45°的邊緣信息；逆時(shí)針旋轉(zhuǎn)135°的ARWHT算子檢測(cè)結(jié)果(圖2c)著重突出了0°、45°、90°的邊緣信息(45°邊緣信息絕對(duì)值最大)，有效過(guò)濾了135°邊緣信息。

圖2 正八邊形模型及其ARWHT算子邊緣檢測(cè)結(jié)果(a)正八邊形模型； (b)逆時(shí)針旋轉(zhuǎn)45°； (c)逆時(shí)針旋轉(zhuǎn)135°

為了驗(yàn)證本文方法提取復(fù)雜邊緣特征信息的有效性，對(duì)經(jīng)典的Lena圖像提取邊緣信息(圖3)。結(jié)果表明：①原始圖像(圖3a)背景中含有較多不同角度的線條，有利于驗(yàn)證算法的有效性。②本文算法弱化了與算子旋轉(zhuǎn)角度一致的邊緣信息，同時(shí)突顯了與其垂直的邊緣信息(圖3b、圖3c、圖3e、圖3f)，且邊緣信息提取效果優(yōu)于常規(guī)希爾伯特變換(圖3d)。如0°算子對(duì)帽子的陰影部分表現(xiàn)更好(圖3b紅色長(zhǎng)方形框處)，90°算子提取的帽檐上邊緣明顯更連續(xù)(圖3c黃色箭頭處)，且唇部線條也更突出。③由45°(圖3e)和135°(圖3f)算子檢測(cè)結(jié)果可見(jiàn)，對(duì)于背景中不同方向的邊緣信息，不同旋轉(zhuǎn)角度算子側(cè)重提取不同方向的邊緣信息(綠色和藍(lán)色箭頭處)。④由于眉毛角度的不同，不同旋轉(zhuǎn)角度的算子檢測(cè)結(jié)果對(duì)眉毛的表現(xiàn)效果也有差異(紅色橢圓框處)。

圖3 Lena圖像及其ARWHT算子邊緣檢測(cè)結(jié)果(a)原始圖像； (b)不旋轉(zhuǎn)； (c)逆時(shí)針旋轉(zhuǎn)90°； (d)常規(guī)方法； (e)逆時(shí)針旋轉(zhuǎn)45°； (f)逆時(shí)針旋轉(zhuǎn) 135°

3 實(shí)際地震資料處理

為了說(shuō)明本文方法在實(shí)際地震資料處理中的應(yīng)用效果，以中國(guó)KL地區(qū)的三維地震資料為例，分別使用旋轉(zhuǎn)0°、45°、90°、135°的算子提取體邊緣信息。由KL地區(qū)沿目的層振幅切片(圖4)可見(jiàn)，目的層內(nèi)河道和斷層分布極為復(fù)雜。

基于ARWHT的實(shí)際地震資料體邊緣檢測(cè)結(jié)果(圖5)表明：①不旋轉(zhuǎn)(圖5a)與逆時(shí)針旋轉(zhuǎn)90°(圖5b)算子的方向相互垂直，檢測(cè)結(jié)果差別明顯，有利于刻畫(huà)細(xì)小裂縫的方向。表現(xiàn)為：前者著重突出了縱向構(gòu)造特征，斷層及裂隙連續(xù)性明顯加強(qiáng)，顯著減弱了橫向構(gòu)造信息；后者的橫向不連續(xù)性信息以及構(gòu)造異常特征明顯加強(qiáng)，明顯壓制了縱向邊緣信息，顯示出原始地震資料中難以用肉眼分辨的大型橫向連續(xù)裂縫(紅色長(zhǎng)方形虛線框處)。另外，還可以看到一條隱蔽的橫向河道，在其他提取結(jié)果中(圖5a、圖5c、圖5d)難以分辨(藍(lán)色箭頭處)。②逆時(shí)針旋轉(zhuǎn)45°算子檢測(cè)結(jié)果(圖5c)突顯了135°的體邊緣信息(黃色箭頭處)；逆時(shí)針旋轉(zhuǎn)135°算子檢測(cè)結(jié)果(圖5d)著重表現(xiàn)了45°的不連續(xù)信息(綠色箭頭處)，還可以追蹤到一條走向大致呈45°方向的蜿蜒河道(綠色點(diǎn)狀線所示)，在其他提取結(jié)果(圖5a、圖5b、圖5c)中難以分辨。③綜合四個(gè)方向的提取結(jié)果來(lái)看，邊緣檢測(cè)結(jié)果的不同區(qū)域(黃色方形虛線框和藍(lán)色橢圓虛線框)數(shù)據(jù)的特征不同，包含的不連續(xù)信息方向特征與ARWHT算子提取邊緣的規(guī)律一致，更好地刻畫(huà)了地震異常信息的空間分布和走向特征。

在利用本文方法處理實(shí)際地震資料時(shí)，應(yīng)考慮數(shù)據(jù)的具體情況選擇算子的旋轉(zhuǎn)角度。對(duì)于已知走向信息的大斷裂或者某一區(qū)域內(nèi)分布方向大致相同的裂縫系統(tǒng)，應(yīng)選擇與其走向或分布方向垂直的算子旋轉(zhuǎn)角度提取邊緣信息，從而加強(qiáng)走向和分布方向的不連續(xù)信息特征；對(duì)于復(fù)雜多變的裂縫系統(tǒng)，其包含各個(gè)方向的信息，此時(shí)應(yīng)選擇多個(gè)旋轉(zhuǎn)角度提取邊緣信息進(jìn)行融合，以顯示更豐富的有效信息，利于地震資料精細(xì)解釋。

為了進(jìn)一步分析本文方法的應(yīng)用效果及優(yōu)勢(shì)，根據(jù)式(13)融合圖5中各方向的邊緣信息提取結(jié)果(圖6c)，并與常規(guī)希爾伯特變換邊緣檢測(cè)結(jié)果(圖6a)、RMS振幅體切片(圖6b)對(duì)比、分析。結(jié)果表明：圖6c綜合了圖5各結(jié)果的優(yōu)點(diǎn)，包含了各個(gè)方向的體邊緣特征信息(圖中紅色虛線框與箭頭處)；與圖6a和圖6b相比，圖6c的有效信息更全面，不僅清楚地展示了河道和斷裂系統(tǒng)的分布特征，還刻畫(huà)了裂縫不連續(xù)地質(zhì)異常的邊界信息，有利于地震資料的精細(xì)解釋。

為了更清晰地展示河道、斷層、斷裂帶等信息，將本文方法的裂縫檢測(cè)結(jié)果(圖6c)與圖6b融合，融合結(jié)果(圖6d)不僅展示了多方向的裂縫信息，還刻畫(huà)了河道走向、裂縫邊界信息，展示了豐富的地質(zhì)信息。

綜上所述，本文方法能突顯地震資料中任意特定方向的地質(zhì)異常信息，更全面地展示了地震不連續(xù)信息的分布特征，檢測(cè)結(jié)果體現(xiàn)了算法的可行性和優(yōu)越性。

圖4 KL地區(qū)沿目的層振幅切片

圖5 基于ARWHT的實(shí)際地震資料體邊緣檢測(cè)結(jié)果(a)不旋轉(zhuǎn)； (b)逆時(shí)針旋轉(zhuǎn)90°； (c )逆時(shí)針旋轉(zhuǎn)45°； (d)逆時(shí)針旋轉(zhuǎn)135°

圖6 實(shí)際地震數(shù)據(jù)不同方法的邊緣檢測(cè)結(jié)果對(duì)比(a)常規(guī)方法； (b)RMS振幅體切片； (c)本文方法； (d)RMS+本文方法

4 結(jié)束語(yǔ)

本文在加窗希爾伯特變換理論的基礎(chǔ)上，構(gòu)建了ARWHT邊緣檢測(cè)算子，將其應(yīng)用于實(shí)際地震資料邊緣檢測(cè)，可以得到任意特定方向的不連續(xù)異常特征，實(shí)現(xiàn)了提取各方向不連續(xù)特征的三維地震資料體邊緣檢測(cè)。

地震資料的ARWHT體邊緣檢測(cè)方法在表征儲(chǔ)層和地質(zhì)結(jié)構(gòu)邊緣信息時(shí)保留了有效的方向特征。因此，可以在實(shí)際地震數(shù)據(jù)中充分挖掘并突顯任意特定方向的地質(zhì)特征。對(duì)于具有不同特征的數(shù)據(jù)，應(yīng)選擇不同旋轉(zhuǎn)角度的算子進(jìn)行邊緣檢測(cè)。實(shí)際地震數(shù)據(jù)測(cè)試表明，ARWHT體邊緣檢測(cè)方法可描述河道、斷層、斷裂帶、構(gòu)造、儲(chǔ)層裂縫分布，強(qiáng)化地震數(shù)據(jù)在任意特定方向的構(gòu)造信息，更好地突出地質(zhì)構(gòu)造特征，為解釋構(gòu)造、了解不連續(xù)地質(zhì)異常信息及其展布規(guī)律、預(yù)測(cè)儲(chǔ)層等提供依據(jù)。