組合型方位體曲率分析方法

周元茂 陳學華 羅 鑫王開華 呂丙南 李 泂

(①成都理工大學油氣藏地質及開發工程國家重點實驗室，四川成都 610059； ②成都理工大學地球探測與信息技術教育部重點實驗室，四川成都 610059； ③中國水利水電第七工程局有限公司試驗檢測研究院，四川成都 610059)

·綜合研究·

組合型方位體曲率分析方法

周元茂①②③陳學華*①②羅 鑫①②王開華③呂丙南①②李 泂①②

(①成都理工大學油氣藏地質及開發工程國家重點實驗室，四川成都 610059； ②成都理工大學地球探測與信息技術教育部重點實驗室，四川成都 610059； ③中國水利水電第七工程局有限公司試驗檢測研究院，四川成都 610059)

考慮到儲層和地質構造的復雜多樣性，為了突出儲層結構及地質構造特征及其方向性，如儲層的空間分布和內部的精細結構、斷層和裂縫系統、褶皺等，在經典的高斯濾波方法基礎上，推導了各向異性高斯迭代平滑濾波方法及其實現算法，同時建立各向異性地震體曲率分析方法，將兩者聯合應用，形成了地震資料的組合型方向體曲率分析方法。實際地震資料的處理結果表明，本文方法可在壓制地震背景干擾的同時突出具體方向的地質異常信息。

去噪 各向異性 高斯濾波 體曲率屬性

1 引言

數字濾波一直是地震勘探領域重要的研究課題之一，有效的地震資料濾波方法一直是人們的研究重點[1-4]。傳統的高斯濾波器是各向同性的，濾波時對圖像平滑區和邊緣過渡區進行了均衡處理，這樣就會導致過濾噪聲后圖像邊緣變得模糊甚至邊緣消失[5-9]。為了解決這一問題，研究者們做了大量工作，Van Ginkel等[10]提出了去卷積的方法提高高斯濾波的角分辨率，但是該方法需要用到復雜的傅里葉去卷積算法；Geusebroek等[11]提出了一種快速各向異性高斯濾波方法，通過在兩個不同方向上選用不同的高斯尺度，使濾波器在去噪的同時可以較好地保留圖像邊緣等重要信息，并利用高斯函數的可分解性將濾波器沿長軸和短軸方向分解為兩個一維濾波器與圖像進行卷積，使計算簡化，提高了效率。

地震屬性是對地震資料的幾何學、運動學、動力學及統計學特征的一種測量和描述，它反映不同的地質信息，是刻畫、描述地層結構、巖性及物性等地質特征的方法[12，13]，因此在油藏識別和儲層預測中起著重要作用。隨著20世紀80年代后期三維地震勘探的發展及普及，地震屬性分析方法逐漸從二維面空間拓展到三維體空間。 Roberts[14]詳述了曲率屬性的基本理論，提出了第一代曲率分析方法——層面曲率屬性的計算和工作流程，表明曲率屬性對斷層和裂縫走向等構造的幾何特征的提取十分有效，為曲率屬性在地震資料構造解釋中的推廣和應用奠定了基礎； Al-Dossary等[15]利用地震數據體所包含的空間方位信息，得到了第二代曲率分析方法——體曲率屬性，并在此基礎之上采用分數傅里葉變換實現了體曲率屬性的多尺度分析； Chen等[16]采用小波包分解實現了時間或深度方向上的多尺度分析，并通過多尺度體曲率屬性的融合顯示獲得了更加豐富的構造細節； Chopr等[17-19]則將曲率屬性與相干屬性融合顯示以進行構造上的識別和解釋，并將地震幾何屬性應用于地震資料預處理中，用以監控地震數據處理的質量。

上述方法在對地震數據去噪時，由于不能對地質構造信息進行選擇性識別，在地震數據去噪以及地震屬性的提取時，在識別某一具體方向(如特定方向)的地質異常信息過程中會受到其他方向細微構造(如裂縫、小斷層等)的干擾。本文的目的是在壓制地震數據背景噪聲的同時，突出特定方向的地質構造特征，提高地震資料信噪比，以便獲得更精確的地質構造信息，達到高分辨率構造精細識別效果。

2 方法原理

2.1 各向異性高斯濾波器

傳統高斯濾波的基本原理是將高斯核函數與原始信號進行卷積[20，21]。以二維高斯濾波為例，其濾波器表達式為

(1)

式中σ為高斯函數的基本方差，它決定了高斯核寬度，與平滑程度有直接關系，σ越大，高斯窗口越寬，平滑程度越好。在二維空間中，高斯函數是旋轉對稱的(圖1a)，對一個圖像的各個方向上的平滑效果

是均衡的(即各向同性)。對于式(1)，當在x、y方向選取不同的尺度時，高斯函數將不再對稱，如圖1b所示，其在x、y平面上的投影為橢圓(圖2a)。橢圓型高斯函數為

(2)

將該橢圓沿其軸線旋轉一個角度θ后，如圖2b所示，可得濾波算子

(3)

其中θ為u軸與x軸之間的夾角。將式(3)中Gθ變換到頻率域，可得

Gθ(ωx,ωy;σu,σv,θ)

(4)

圖1 二維各向同性和二維各向異性高斯濾波器

圖2 各向異性高斯濾波器橢圓及其坐標系

為了使該二維濾波器能分解為沿x方向和沿y方向的一維濾波器，將上式改寫為

Gθ(ωx,ωy;σu,σv,θ)

(5)

為了使濾波方向沿T(y=xtanφ)方向濾波而不是沿y方向濾波，如圖2c所示，在T方向上將式(5)改寫為

(6)

利用傅里葉逆變換將上式變換到時間域

(7)

(8)

其中σx為高斯函數沿x方向的標準差，即

(9)

Gθ(x,T;σu,σv,θ)=Gx(x;σx)*Gφ(T;σφ)

(10)

式中

(11)

通過將二維高斯濾波分解為兩個方向上的一維濾波以實現各向異性高斯濾波，提高了運算速度，節省了計算卷積的時間。同時采用加權平滑濾波的方法可有效解決濾波過程中出現邊緣模糊的問題。

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

在幾何地震學中，三維地震反射體在空間上的任一反射點r(z，x，y)可以認為是時間標量u(t，x，y)，那么梯度grad(u)反映的是反射面沿著不同方向的變化率，即反射面沿著方向矢量所在的法截面截取的一階導數，其結果為該反射點的視傾角向量。

(12)

式中px、py、pz分別為沿x、y和z方向的視傾角分量，如圖3所示。其中：ρ為真傾角；ζ為方位角；a為沿著反射面的單位矢量傾角；ψ為走向。

圖3 三維空間中傾角、方位角與視傾角的關系[22]

一個三維地震數據體可以先轉化為傾角數據體，再根據傾角數據體計算其中任意點的曲率。根據復數道分析中的瞬時頻率和瞬時波數計算三維地震數據體的視傾角，瞬時頻率的計算需要對初始地震數據及其希爾伯特變換進行兩次微分運算[23，24]，而微分運算對噪聲較為敏感，為了解決這一問題，Claerbout[25]提出估算瞬時頻率的方法

(13)

式中：φ為瞬時相位；s(t)為初始解析信號； Δt為信號的采樣間隔。

本文采用中心差分格式計算，將虛部運算化到實部，得到瞬時頻率計算公式

(14)

式中：g(t)為t方向初始地震信號；h(t)為其相應的希爾伯特變換。與式(13)相比，式(14)具有更高的精度和更小的計算量，更適用于較大型的三維地震數據。

同理，可得到x方向和y方向的瞬時波數kx、ky的計算公式為

(15)

式中：g(x)、g(y)分別為x方向和y方向初始地震信號；h(x)、h(y)為其相應的希爾伯特變換。

可由瞬時頻率ω和瞬時波數kx、ky計算出x方向及y方向的視傾角分量px和qy

(16)

在此基礎上，構造新的方位視傾角分量

(17)

式中β為根據需要突出的地質構造所在的方向而選取的角度。通過重新計算方位視傾角分量(pβ,qβ)得到的體曲率屬性可以反映某一特定方向上的構造信息。

對于三維空間中某點的曲率值，通過與其相鄰道和樣點的視傾角值擬合出的空間曲面方程計算得到。根據最小二乘逼近原理，得到二次曲面方程為

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

(18)

對上式方程兩邊求微分，并將pβ、qβ帶入其中，得到了趨勢面方程系數

(19)

綜上可見，三維數據體中過任意點的空間曲面可由三維時窗下各個數據點的視傾角計算出來。過空間中某點的曲面可以產生任意多個曲率。最大正曲率Kpos和最小負曲率Kneg對線性構造敏感，因此在描述斷層、裂縫等構造上最為有效[26，27]，是用于實際地震資料構造解釋的主要屬性。其表達式為

(21)

3 模型試算與分析

為了驗證各向異性高斯濾波的效果，設計一個切片模型(圖4)。圖4a為無噪聲的原始模型，其中，虛線可表示地震資料中儲層的斷層或裂縫，圖4b為加入隨機噪聲后的模型。使用各向異性高斯濾波算法對圖4b的數據進行濾波處理，得到四個方向的濾波結果(圖5)。圖中θ表示選取的濾波方向(如箭頭所示)，其中，圖5a為θ=0°時的結果，從圖中可見，圖4b中存在的大量隨機噪聲已明顯減少，原本橫向的虛線變成連續的實線且更為粗壯，其他方向的虛線則被明顯壓制。同理，圖5b～圖5d是方向角分別為45°、90°和135°時的濾波結果，突出的信息與濾波取向相對應。說明該方法不僅能達到去噪的效果，而且能更好地反映構造中不連續性信息，證明了方法的可行性。

圖6為各向異性高斯濾波后求取的各向異性曲率(最小負曲率)屬性結果。在求取各向異性體曲率屬性時，選取與高斯濾波相對應的方位角。由圖可見，各向異性體曲率屬性效果明顯，對應方向的有效信息得到了進一步加強，證明各向異性高斯濾波與各向異性體曲率相結合的方法是可行的。

圖4 模型切片

圖5 經不同取向的各向異性高斯濾波處理后的結果

4 實例分析

為了驗證基于地震資料的組合型方向體曲率分析方法的應用效果，對LH地區的三維疊后地震數據體進行了高斯濾波及體曲率屬性提取。該地區受多期構造作用的控制，斷裂系統發育。圖7a為原始地震數據目的層振幅切片，該目的層位于珠江組上段，屬于典型的碳酸鹽巖儲層，儲層內部空隙發育，孔滲關系復雜，非均質性強； 圖7b為濾波前求取的最小負曲率沿層切片，由圖可見，隨機噪聲干擾嚴重，給后期的構造解釋、斷裂及展布規律的識別、儲層預測等帶來了困難； 圖7c為地震數據體經傳統高斯濾波后的最小負曲率沿層切片，可看出原先圖7b中的噪聲已被明顯壓制，圖中的構造信息更加清晰。

利用本文提出的算法處理該三維疊后地震數據體，得到不同濾波方向的三維地震體曲率屬性(最小負曲率)沿層切片。圖8a～圖8c是各向異性高斯迭代平滑濾波與各向異性體曲率相結合的處理結果，其各向異性高斯濾波的角度θ分別為0°、0°、90°，以及各向異性體曲率的方位β分別為0°、90°、90°。

圖6 不同取向的各向異性體曲率屬性(最小負曲率)結果

圖7 LH地區的振幅和體曲率屬性的沿層切片

圖8 組合型方向體曲率處理后的沿層切片

對比圖8a與圖8c，二者差別明顯，著重突出的構造特征方向不同。圖8a的濾波方向角為0°，從圖中可見，該切片橫向不連續性且構造異常特征得到了明顯的加強，其背景噪聲干擾相對于圖7b得到了明顯壓制。圖8c的濾波角度為90°，該切片著重突出了縱向構造特征，其縱向斷層及裂隙連續性得到明顯加強，橫向構造及背景噪聲干擾大幅減弱。圖8b為二者不同方向相結合的處理結果，在過濾大部分背景噪聲的同時突出主要地質構造特征(如較大的斷層、褶皺、落水洞等)。相對于傳統高斯濾波及體曲率分析方法，本文提出的組合型方向體曲率分析方法有其獨特的優勢，在壓制背景噪聲干擾的同時，能夠更準確地突出地震資料中特定方向上的地質異常信息，大幅減弱了背景噪聲及其他方向地質信息的干擾，體現出該算法的可行性和優越性。

5 結論

本文提出了一種改進的各向異性高斯濾波方法。基于常規各向同性高斯濾波算法，推導了各向異性高斯迭代平滑濾波公式。將二維高斯濾波算法分解為兩個方向上的一維濾波，通過改變濾波角度來實現定向濾波的目的，同時提高了運算速度，節省了計算卷積的時間。在三維體曲率屬性的計算中加入方位參數，推導出各向異性體曲率屬性分析方法，該方法可定向選取需要突出地質構造信息。將各向異性高斯濾波方法與各向異性體曲率屬性分析方法聯合應用，形成地震資料的組合型方向體曲率屬性分析方法。該方法可對地震數據中的某一特定方向的構造信息進行兩次加強，極大地突出了該方向的地質構造特征，可以更好地反映特定方向的構造細節、不連續性，為構造解釋、斷裂及展布規律的識別、儲層預測等提供更好支持。需要指出的是，該方法在地震數據濾波及曲率屬性提取時，需要根據地質構造情況，人為地選擇濾波角度及提取曲率屬性的方位角，還無法做到自適應濾波和自適應提取曲率屬性，這點需要進一步改進。

[1] 張俊華.地震資料去噪方法——原理、算法、編程及應用.山東青島：中國石油大學出版社，2011.

[2] 楊陽，曲壽利，印興耀.沿層地震屬性數據空間濾波去噪方法.石油地球物理勘探，2007，42(2)：208-211.

Yang Yang，Qu Shouli，Yin Xingyao.Using spatial filtering for denoise of horizon-oriented seismic attri-butes data.OGP，2007，42(2)：208-211.

[3] 張華，陳小宏，楊海燕.地震信號去噪的最優小波基選取方法.石油地球物理勘探，2011，46(1)：70-75.

Zhang Hua，Chen Xiaohong，Yang Haiyan.Optimistic wavelet basis selection in seismic signal noise elimination.OGP，2011，46(1)：70-75.

[4] 羅鑫，陳學華，齊迎凱.基于匹配濾波的保幅動校拉伸校正方法.石油地球物理勘探，2016，51(4)：654-660.

Luo Xin，Chen Xuehua，Qi Yingkai.A NMO stretch correction based on an amplitude-preserved matched filtering.OGP，2016， 51(4)：654-660.

[5] 姒紹輝，胡伏原，顧亞軍.一種基于不規則區域的高斯濾波去噪算法.計算機科學，2014，41(11)：313-316.

Si Shaohui，Hu Fuyuan，Gu Yajun.Improved denoising algorithm based on non-regular area Gaussian filtering.Computer Science，2014，41(11)：313-316.

[6] 伍鵬，賀振華，陳學華.二維高斯迭代平滑濾波曲率屬性及其應用．地球物理學進展，2010，25(6)：2144-2149．

Wu Peng，He Zhenhua，Chen Xuehua et al.Curvature attribute of two-dimensional Gaussian iterated smoo-thing filtering and applications.Progress in Geophy-sics，2010，25(6)：2144-2149.

[7] 王小兵，孫久遠，湯海燕.一種高斯噪聲組合濾波方法.佳木斯大學學報(自然科學版)，2011，29(5)：696-698.

Wang Xiaobing，Sun Jiuyuan，Tang Haiyan.A mixed filter method on Gaussian noise.Journal of Jiamusi University (Natural Science Edition)，2011，29(5)：696-698.

[8] 馮芳.高斯濾波運算的幾種并行實現方法.蘭州工業學院學報，2015，22(5)：57-60.

Feng Fang.Several methods on Gaussian algorithm.Journal of Lanzhou Institute of Technology，2015，22(5)：57-60.

[9] 高萌，劉吉臻，王瑞琪.基于自適應高斯濾波的電站歷史數據穩態檢測方法.動力工程學報，2014，34(9)：708-714.

Gao Meng，Liu Jizhen，Wang Ruiqi et al.Steady-state detection of power plant history data based on adaptive Gauss filter.Journal of Chinese Society of Power Engineering，2014，34(9)：708-714.

[10] Van Ginkel M，Verbeek P W，Vliet L J.Improved orientation selectivity for orientation estimation.Proceedings of the 10th Scandinavian Conference on Image Analysis，1997，533-537.

[11] Geusebroek J M，Smeulders A W M，Weijer J V D.Fast anisotropic Gauss filtering.IEEE Transactions on Image Processing，2003，12(8)：938-943.

[12] 王世瑞，王樹平，狄幫讓.基于地震屬性特征的河道砂體預測方法.石油地球物理勘探，2009，44(3)：304-313.

Wang Shirui，Wang Shuping，Di Bangrang.Prediction of channel sand body based on seismic attributes.OGP，2009，44(3)：304-313.

[13] 李建雄，崔全章，魏小東.地震屬性在微斷層解釋中的應用.石油地球物理勘探，2011，46(6)：925-929.

Li Jianxiong，Cui Quanzhang，Wei Xiaodong.Application of seismic attributes in micro-fault interpretation.OGP，2011，46(6)：925-929.

[14] Roberts A.Curvature attributes and their application to 3D interpreted horizons.First Break，2001，19(2)：85-100.

[15] Al-Dossary S and Marfurt K J.3D volumetric multispectral estimates of reflector curvature and rotation.Geophysics，2006，71(5)：41-51.

[16] Chen X H，Yang W，He Z H.The algorithm of 3D multi-scale volumetric curvature and its application.Applied Geophysics，2012，9(1)：65-72.

[17] Chopra S and Marfurt K J.Volumetric curvature attributes add value to 3D seismic data interpretation.The Leading Edge，2007，26(7)：856-867.

[18] Chopra S，Marfurt K J.Integration of coherence and volumetric curvature images.The Leading Edge，2010，29(9)：1092-1107.

[19] Chopra S，Marfurt K J.Coherence and curvature attributes on preconditioned seismic data.The Leading Edge，2011，30(4)：386-393.

[20] 錢曉亮，郭雷，余博.基于目標尺度的自適應高斯濾波.計算機工程與應用，2010，46(12)：14-16.

Qian Xiaoliang，Guo Lei，Yu Bo.Adaptive Gaussian filter based on object scale.Computer Engineering and Applications，2010，46(12)：14-16.

[21] 章為川，程冬，朱磊.基于各向異性高斯核的多尺度角點檢測.電子測量與儀器學報，2012，26(1)：37-42.

Zhang Weichuan，Cheng Dong，Zhu Lei.Multi-scale corner detection based on anisotropic Gaussian kernels.Journal of Electronic Measurement and Instrument，2012，26(1)：37-42.

[22] Marfurt K J.Robust estimates of 3D reflector dip and azimuth.Geophysics，2006，71(4)：29-40.

[23] 李紅星，饒溯，陶春輝.廣義希爾伯特變換地震邊緣檢測方法研究.石油地球物理勘探，2015，50(3)：490-494.

Li Hongxing，Rao Su，Tao Chunhui.Seismic edge detection method based on generalized Hilbert transform.OGP，2015，50(3)：490-494.

[24] 李斌，陳學華，賀振華.基于多尺度加窗希爾伯特變換的地震資料體邊緣檢測.石油物探，2015，54(3)：345-349.

Li Bin，Chen Xuehua，He Zhenhua et al.Seismic data 3D edge detection based on multi-scale windowed Hilbert transform.GPP，2015，54(3)：354-359.

[25] Claerbout J F.Fundamentals of Geophysical Data Processing.McGraw-Hill Book Company，New York，USA，1976.

[26] 楊威，賀振華，陳學華.三維體曲率屬性在斷層識別中的應用．地球物理學進展，2011，26(1)：110-115.

Yang Wei，He Zhenhua，Chen Xuehua.Application of three-dimensional volumetric curvature attributes to fault identification.Progress in Geophysics，2011，26(1)：110-115.

[27] 楊威，賀振華，陳學華.結構方位濾波在體曲屬性中的應用.石油物探，2011，50(1)：27-32.

Yang Wei，He Zhenhua，Chen Xuehua.Application of structural azimuth filtering in volumetric curvature.GPP，2011，50(1)：27-32.

*四川省成都市成華二仙區橋東三路1號成都理工大學油氣藏地質及開發工程國家重點實驗室，610059。Email：chen_xuehua@163.com

本文于2016年10月13日收到，最終修改稿于2017年9月10日收到。

本項研究受國家自然科學基金項目(41374134、41574130)、“十三·五”國家科技重大專項(2016ZX05014-001-009)和四川省青年科技創新研究團隊項目(2016TD0023)聯合資助。

1000-7210(2017)06-1253-08

周元茂，陳學華，羅鑫，王開華，呂丙南，李泂.組合型方位體曲率分析方法.石油地球物理勘探，2017,52(6):1253-1260.

P631

A

10.13810/j.cnki.issn.1000-7210.2017.06.016

(本文編輯：宜明理)

周元茂 碩士，1991年生；2014年畢業于桂林理工大學勘查技術與工程專業，獲學士學位；2017年畢業于成都理工大學地球物理學院地球探測與信息技術專業，獲碩士學位；主要從事地震數據的去噪以及幾何屬性提取方面的研究。