基于結(jié)構(gòu)導(dǎo)向的梯度屬性邊緣檢測技術(shù)

宋建國，孫永壯，任登振

1 中國石油大學(xué)（華東）地球科學(xué)與技術(shù)學(xué)院，青島 266580

2 中海石油（中國）有限公司上海分公司研究院，上海 200030

3 中海油田服務(wù)股份有限公司，天津 300451

1 引 言

地震資料的不連續(xù)信息在斷層、河道、鹽丘等地質(zhì)體的解釋過程中具有重要地位，因此，實現(xiàn)對地震數(shù)據(jù)不連續(xù)性信息的檢測是一項重要任務(wù).目前，相干體、曲率體技術(shù)［1－2］等很多方法都能夠?qū)崿F(xiàn)地震資料不連續(xù)性檢測.

邊緣檢測是起源于圖像處理領(lǐng)域的一門重要技術(shù)，主要用于檢測圖像中灰度值突變區(qū)域.引申到地震勘探領(lǐng)域，地震資料中的不連續(xù)性區(qū)域反映在圖像中即為邊緣特征［3］.因此將邊緣檢測技術(shù)用于地震資料不連續(xù)信息的檢測非常有意義.

邊緣檢測技術(shù)在地震資料處理中的應(yīng)用已經(jīng)有了十幾年的歷史.眾多學(xué)者對邊緣檢測技術(shù)的實際應(yīng)用做了大量研究，提出了多種基于不同算法的、適用于不同尺度地質(zhì)體的邊緣檢測技術(shù)，如基于高階偽希爾伯特變換算法的邊緣檢測技術(shù)、小尺度特征地震邊緣檢測技術(shù)［4－5］等.國內(nèi)邊緣檢測技術(shù)方面的文章多數(shù)應(yīng)用于二維地震圖像，包括薄地層識別、剖面同相軸自動識別及檢測、初至波同相軸提取、裂縫檢測以及砂體邊界識別等［3，6－11］，對于其在三維地震資料中的應(yīng)用未作出深入研究；國外方面，一些學(xué)者提出將邊緣檢測技術(shù)應(yīng)用到三維地震資料處理中的方法，涉及到鹽丘邊界識別、斷層邊界識別等［12－14］.

梯度極大值邊緣檢測具有原理簡單、計算量小的優(yōu)點(diǎn)，被廣泛應(yīng)用于圖像處理領(lǐng)域，代表算子有Sobel算子、Facet模型曲面擬合算子等.上述常規(guī)算子如果直接用于三維地震數(shù)據(jù)地質(zhì)體邊界識別，并不能取得良好的效果，其原因有兩點(diǎn)：1）常規(guī)算子將地震資料的不連續(xù)性信息不加篩選的全部檢測出來，這導(dǎo)致地震資料中連續(xù)性很好的地層邊界信息對小斷層、微裂縫邊界的檢測有著很強(qiáng)的干擾.2）常規(guī)邊緣檢測算子檢測梯度值為絕對梯度值，這使得小地質(zhì)體邊界的弱反射信息淹沒在強(qiáng)背景地層邊界反射信息中.

Y.Luo（1996）提出了差分算子思想，該思想對沿層邊緣檢測的實現(xiàn)有著重要的指導(dǎo)意義.基于該思想，Aqrawi（2011）利用互相關(guān)方法求取地層傾角，然后將傾角約束Sobel邊緣檢測算子應(yīng)用于斷層和鹽丘邊界檢測中，取得很好效果.本文借鑒Aqrawi（2011）思想，利用梯度結(jié)構(gòu)張量方法［15］估計地層傾角和方位角，給出結(jié)構(gòu)導(dǎo)向梯度屬性算子計算公式，通過與結(jié)構(gòu)導(dǎo)向方差屬性檢測結(jié)果對比分析，驗證了該方法的可行性.

2 方法原理

2.1 常規(guī)Facet模型曲面擬合算子

梯度極大值邊緣檢測技術(shù)通過利用差分模板對原始地震圖像進(jìn)行檢測來獲得梯度模值圖像.由于邊緣存在的區(qū)域通常表現(xiàn)為局部梯度極大值，可以根據(jù)模值圖像中模值的大小判斷邊緣發(fā)育情況.常規(guī)5×5Sobel差分模板如圖1所示.Sobel模板數(shù)值的確定往往憑借經(jīng)驗，因此不易進(jìn)行擴(kuò)展.

1984年，Haralick提出了基于多項式曲面擬合的Facet模型邊緣檢測算子［16］.Facet檢測算子的基本思想可以描述為如下過程：

圖1 Sobel算5×5窗口大小兩個方向卷積模板Fig.1 5×5convolution template of Sobel operator for two direction

I（r，c）表示局部鄰域S內(nèi)數(shù)據(jù)，利用式（1）所示二元三次多項式對I（r，c）進(jìn)行最小二乘擬合求取擬合系數(shù)，利用求得的擬合系數(shù)構(gòu)建擬合函數(shù)f（r，c），最終對擬合函數(shù)f（r，c）求取梯度模值.這樣，原來直接對圖像鄰域數(shù)據(jù)進(jìn)行的模板操作轉(zhuǎn)化為對擬合函數(shù)的處理.該方法能夠在抑制噪聲的同時提高邊緣定位精度，而且算子更易擴(kuò)展.擬合函數(shù)為）

式中ki為擬合系數(shù).為了方便擬合系數(shù)的計算，通常將式（1）用正交多項式組合來代替，即：

其中g(shù)i（r，c）代表第i個正交多項式.利用最小二乘法可以得到式（3）：

可見每一個Ki都可以單獨(dú)由I（r，c）的線性組合表達(dá)，其中Gi為Ki對應(yīng)的卷積模板.對比式（1）和式（2）可知，ki可以由Ki線性表示.因此，為了求取f（r，c）的梯度幅值，只需要求得k2和k3.k2和k3對應(yīng)的卷積模板Gr和Gc則可以由Gi卷積模板線性表示.通過上述轉(zhuǎn)換和推導(dǎo)，曲面擬合方法求取梯度最終轉(zhuǎn)化為利用兩個相互正交的模板Gr和Gc對原始圖像局部鄰域進(jìn)行模板卷積的過程，5×5窗口大小模板系數(shù)如圖2所示：

利用背斜模型驗證Facet模型算子檢測效果.圖3所示為背斜模型，圖中加入了5db高斯白噪，圖3a為抽取的一條縱測線剖面，圖3b為時間切片.為了對比分析Facet檢測結(jié)果，同時利用Sobel算子進(jìn)行檢測，處理中Facet模板和Sobel模板均為3×3.

Sobel算子檢測結(jié)果和Facet算子檢測結(jié)果分別如圖4和圖5所示.對比圖4和圖5可以發(fā)現(xiàn)，F(xiàn)acet模型算子檢測結(jié)果不但具有定位精度高的優(yōu)點(diǎn)，而且相比其他算子抑噪性能更好，在噪聲壓制和定位精度上都有很大提高.

2.2 常規(guī)算子改進(jìn)

常規(guī)邊緣檢測算子輸出絕對梯度模值，導(dǎo)致地震資料中小地質(zhì)體邊界弱反射信息在地震剖面上沒有良好的響應(yīng)特征.基于此，通過引入能量歸一化因子對常規(guī)算子進(jìn)行改進(jìn)，以r方向為例，絕對梯度值計算公式由式（4）

改進(jìn)為式（5）：

圖2 Facet模型算子5×5窗口大小兩個方向卷積模板Fig.2 5×5convolution template of Facet operator for two direction

利用常規(guī)Facet算子和能量歸一化后的Facet模型算子對某地震剖面進(jìn)行處理，處理前后如圖6所示.從原始剖面上可以看到，上覆巨厚石膏夾鹽巖層經(jīng)過后期構(gòu)造運(yùn)動，地層褶皺變形，在該地層下方存在的若干礁體，由于石膏夾鹽巖層的波阻抗差太大，地震波能量難以穿透巨厚強(qiáng)反射界面，導(dǎo)致巖層下方礁體反射信息微弱.常規(guī)Facet模型算子對下層的弱信息無法突出顯示，能量歸一化Facet算子檢測結(jié)果中，強(qiáng)的上覆地層信息和弱的礁體反射信息都得到了很好的顯示.

引入能量歸一化使得弱反射信號得到了增強(qiáng)顯示，但是強(qiáng)的背景地層信息仍然存在，本文采用沿層思想，利用結(jié)構(gòu)張量矩陣估計局部地層傾角和方位角信息，然后沿著局部地層發(fā)育方向提取分析窗口，再利用能量歸一化梯度算子進(jìn)行處理，以實現(xiàn)邊界不連續(xù)性信息的識別和刻畫，在下一節(jié)對該方法展開詳細(xì)論述.

2.3 結(jié)構(gòu)導(dǎo)向梯度屬性的計算方法

為了研究相干濾波，1999年 Weickert提出了各向異性非線性擴(kuò)散濾波技術(shù)，處理過程中通過利用結(jié)構(gòu)張量的特征向量控制擴(kuò)散方向，同時利用結(jié)構(gòu)張量的特征值控制特征方向上的擴(kuò)散量，實現(xiàn)對原始地震圖像的平滑處理，增強(qiáng)具有一致性的同相軸及其所反映的重要地質(zhì)結(jié)構(gòu)，同時實現(xiàn)地震圖像邊緣信息的保護(hù).處理中在地震資料局部結(jié)構(gòu)特征分析的基礎(chǔ)上進(jìn)行擴(kuò)散參數(shù)的調(diào)整，從而實現(xiàn)邊界信息完整性的保持和噪音最大限度的消除［17］.梯度結(jié)構(gòu)張量算法所構(gòu)建的梯度結(jié)構(gòu)張量矩陣包含了豐富的局部地層信息，因此本文利用該算法估計局部地層結(jié)構(gòu)傾角和方位角信息.令f（x，y，t）代表某三維地震數(shù)據(jù)體，如疊后偏移處理后的數(shù)據(jù)體或瞬時相位等地震屬性體等，梯度結(jié)構(gòu)張量矩陣定義為

其中Gρ表示尺度為ρ的三維高斯函數(shù)，＊表示褶積，fσ＝Gσ＊f，Δfσ為梯度.

梯度結(jié)構(gòu)張量的求取過程可以用如下流程描述：

（1）利用尺度為σ的三維高斯函數(shù)對f進(jìn)行濾波，得到平滑后的數(shù)據(jù)體fσ，σ可以根據(jù)所要研究的地質(zhì)體大小和數(shù)據(jù)體噪聲水平綜合確定，σ值越小，分析得到的地質(zhì)體信息越精細(xì)，但噪聲影響越大；σ值越大，分析得到的地質(zhì)體信息越粗糙，但噪聲壓制效果越好；

（2）對fσ求取梯度矢量和梯度張量；

（3）利用尺度為ρ的高斯函數(shù)對梯度張量進(jìn)行高斯濾波，進(jìn)一步抑制噪聲，使得到的局部地層結(jié)構(gòu)信息更準(zhǔn)確.

從求取過程可以看出，最終局部地層結(jié)果的精細(xì)程度由尺度σ和ρ共同決定，一般選取較小的σ值（如0.5）和較大的ρ（如1.5）.

不難證明實對稱矩陣Tρ為半正定矩陣，它的特征值為不小于零的實數(shù)，并且不為零的特征值對應(yīng)的特征向量相互正交.利用矩陣分解，Tρ可以變?yōu)?/p>

其中λ1≥λ2≥λ3≥0為Tρ的三個特征值，其各自對應(yīng)的特征向量為V1、V2和V3.這些特征值和特征向量具有明顯的物理意義，三個特征向量可以構(gòu)成一個局部三維坐標(biāo)系，V1是局部地層結(jié)構(gòu)變化最快的方向，可以近似代表局部地層的法線方向（圖7），變化率為λ1；V2和V3所構(gòu)成的平面垂直于V1，該平面上局部地層結(jié)構(gòu)變化最快的方向為V2，變化率為λ2，變化最慢的方向為V3，變化率為λ3.

通過計算每個目標(biāo)點(diǎn)所在層位的局部法線方向V1，可以求取x（橫測線）方向和y（縱測線）方向的延遲時間τx（x，y，t）和τy（x，y，t），根據(jù)延遲時間可以提取局部地層分析窗口，在該分析窗口進(jìn)行能量歸一化算子處理，就可以消除背景地層邊界信息的干擾.

圖7 局部地層結(jié)構(gòu)的法線方向Fig.7 The normal direction of local formation structure

為了提高梯度屬性的解釋精度，抑制地震數(shù)據(jù)體中的噪聲，文中采用多窗口劃分思想進(jìn)行局部地層分析窗口的提取.1976年Kuwahara首次提出多窗口濾波思想，其后，眾多學(xué)者通過研究，給出了不同的多窗口劃分方案.以5×5窗口為例，圖8給出了Y.Luo等的窗口劃分方案和Rugumira Georgia等［18］的窗口劃分方案.本文選用Y.Luo等提出的窗口劃分方法，即圖8a所示，將5×5大小的窗口劃分為9個小多邊形，包括中央一個矩形，四個頂點(diǎn)處的四個六邊形和四條邊上的四個五邊形.

利用延遲時間，可以提取局部沿層數(shù)據(jù).延遲時間的求取可以采用以下兩種方法：

（1）利用V1直接求取τx（x，y，t）和τy（x，y，t）；

（2）根據(jù)幾何關(guān)系，計算得到傾角數(shù)據(jù)體θ（x，y，t）和方位角數(shù)據(jù)體φ（x，y，t），將θ（x，y，t）和φ（x，y，t）轉(zhuǎn)化為沿x方 向和y方向的延遲 時 間τx（x，y，t）和τy（x，y，t），轉(zhuǎn)化公式為

在進(jìn)行社區(qū)醫(yī)院的工作中，對患者進(jìn)行自我管理模式健康教育是很有必要進(jìn)行的，社區(qū)醫(yī)院應(yīng)該充分認(rèn)識到這一點(diǎn)，做好對患者的安全教育工作，使其配合醫(yī)院的治療，從而增強(qiáng)治療療效，幫助患者恢復(fù)身體健康。

本文采用第二種計算方法，該方法能夠避免較大傾角情況下異常大值延遲時間的產(chǎn)生.

為了驗證Facet檢測算子檢測效果，同時利用方差算子、Sobel算子進(jìn)行處理.本文給出結(jié)構(gòu)導(dǎo)向方差屬性，結(jié)構(gòu)導(dǎo)向Sobel梯度屬性和結(jié)構(gòu)導(dǎo)向Facet模型梯度屬性三種屬性計算公式.

結(jié)構(gòu)導(dǎo)向方差屬性的求取采用以下方法：

（1）分別求取每個局部窗口均值，根據(jù)均值求取方差；

（2）沿法線方向?qū)Ω骶植糠治龃翱谇蟮玫姆讲罱Y(jié)果作和；

（3）按能量對結(jié)果進(jìn)行歸一化處理.

為了增加輸出結(jié)果的穩(wěn)定性，處理中采用垂向2K＋1個樣點(diǎn)的疊加求均值作為輸出，2K＋1一般取一個子波延續(xù)時間內(nèi)的樣點(diǎn)數(shù)，具體計算公式為

圖8 二維多窗口保邊濾波窗口劃分示意圖［18］（a）Y.Luo等人窗口劃分思想，在5×5鄰域內(nèi)包含四個五邊形，四個六邊形和一個矩形；（b）Rugumira Georgia等窗口劃分思想，在5×5鄰域內(nèi)包含四個六邊形和一個八邊形.Fig.8 Window dividing diagram of two－dimensional multi－window edge－preserving filter（a）Window dividing method of Y.Luo，there are four pentagon，four hexagon and one rectangle in 5×5neighborhood；（b）Window dividing method of Rugumira Georgia，there are four hexagon and one octagon in 5×5neighborhood.

具有邊緣保持效果的插值方法有多種，如基于邊緣保持濾波器的多項式擬合插值算法［19］，三次樣條插值算法等.其中三次樣條插值算法在邊緣提取過程中抑制噪聲和邊緣定位精度的綜合性能相對較為優(yōu)秀，因此本文采用三次樣條插值算法獲取近似值.D和的計算公式為

其中Δt為采樣率，（x，y，t）和（x，y，t）分別為與局部窗口中心點(diǎn)所對應(yīng)的縱測線和橫測線相同層位的延遲時間，標(biāo)識含義如下式所示：

在同樣窗口大小和窗口數(shù)目的情況下，結(jié)構(gòu)導(dǎo)向Sobel梯度屬性計算公式為

其中My和Mx代表Sobel算子兩個正交方向模板.結(jié)構(gòu)導(dǎo)向Facet模型梯度屬性計算公式為

其中Wy和Wx代表Facet模型算子兩個正交方向模板.

3 實際資料處理

為了驗證Facet算子檢測結(jié)果，選取實際地震數(shù)據(jù)進(jìn)行處理.選取來自荷蘭境內(nèi)北海F3區(qū)塊三維地震數(shù)據(jù)體的一部分，該區(qū)域內(nèi)發(fā)育大量東西走向斷層，抽取原始資料第1560ms時間切片，如圖9所示.由于原始數(shù)據(jù)體噪聲比較嚴(yán)重，為了抑制噪聲，首先進(jìn)行三維邊緣保持濾波，對濾波后數(shù)據(jù)利用結(jié)構(gòu)張量矩陣求取傾角和方位角數(shù)據(jù)體.

圖10～12為濾波后數(shù)據(jù)體結(jié)構(gòu)導(dǎo)向方差屬性、Sobel梯度屬性和Facet模型梯度屬性1560ms時間切片.圖13為局部放大對比圖.三種屬性的計算采取相同參數(shù)，窗口大小5×5，垂向時窗包含15個樣.

由圖10方差屬性結(jié)果可以看出，原始資料中的大斷層信息基本得到展現(xiàn)，獲取的屬性在噪聲抑制方面效果較好.對比圖11Sobel梯度屬性，Sobel屬性斷層細(xì)節(jié)展現(xiàn)更好，一些弱的斷層得到較好顯示.將以上兩種屬性與本文提出的Facet模型梯度屬性對比分析，F(xiàn)acet梯度屬性具有更好的噪聲抑制效果，細(xì)節(jié)信息最為豐富.

圖9 1560ms處的時間切片F(xiàn)ig.9 Time slice of 1560ms

圖10 結(jié)構(gòu)導(dǎo)向的方差屬性在1560ms處的時間切片F(xiàn)ig.10 Structure directed variance attribute slice of 1560ms

觀察圖13局部放大對比圖，F(xiàn)acet模型梯度屬性對原始斷層展示最為清晰，尤其對于反射信息較弱的小斷裂發(fā)育區(qū)，刻畫信息最為清晰，箭頭所指區(qū)域給出良好展示.

圖13 原始切片和三種屬性不連續(xù)性檢測結(jié)果的局部放大圖Fig.13 Enlargement of original slice and three kinds of un－continuity detection

4 結(jié) 論

為了實現(xiàn)儲層小斷裂、微裂縫等地震不連續(xù)性信息的邊緣檢測，本文對常規(guī)邊緣檢測算子進(jìn)行了改進(jìn)，一是引入能量歸一化因子，突顯了地震資料中小地質(zhì)體邊界弱反射信息，二是利用局部地層結(jié)構(gòu)傾角和方位角信息約束，使邊緣檢測算子在增強(qiáng)弱信息的同時消除背景地層信息的干擾.利用結(jié)構(gòu)導(dǎo)向方差屬性算子，結(jié)構(gòu)導(dǎo)向Sobel梯度屬性算子和結(jié)構(gòu)導(dǎo)向Facet模型梯度屬性算子對實際數(shù)據(jù)進(jìn)行處理，證明結(jié)構(gòu)導(dǎo)向Facet模型梯度屬性邊緣檢測算子檢測結(jié)果最精細(xì)，包含邊界信息最豐富，對小斷層、微裂縫的識別最準(zhǔn)確，可以作為一種地震資料精細(xì)解釋工具.

（

）

［1］ Bahorich M，F(xiàn)armer S.3－D seismic discontinuity for faults and stratigraphic features：The coherence cube.TheLeading Edge，1995，14（10）：1053－1058.

［2］ Wang Y S，Yang P J.Orientation steering fault preserving filtering.SEGTechnicalProgramExpandedAbstracts，2010，29（1）：1550－1554.

［3］ 孫夕平，杜世通.邊緣檢測技術(shù)在河道和儲層小斷裂成像中的應(yīng)用.石油物探，2003，42（4）：469－472.Sun X P，Du S T.Using edge detection technique to image channels and minor faults.GeophysicalProspectingfor Petroleum（in Chinese），2003，42（4）：469－472，doi：1000－1441（2003）04－0469－04.

［4］ 陳學(xué)華，賀振華，黃德濟(jì).地震資料的高階偽希爾伯特變換邊緣檢測.地球物理學(xué)進(jìn)展，2008，23（4）：1106－1110.Chen X H，He Z H，Huang D J.Seismic data edge detection based on higher－order pseudo Hilbert transform.Progressin Geophysics（in Chinese），2008，23（4）：1106－1110.Doi：1004－2903（2008）04－1106－05.

［5］ 孫夕平，周超.小尺度邊緣特征地震檢測技術(shù)研究.石油地球物理勘探，2011，46（1）：121－125.Sun X P，Zhou C.Small－scale edge characteristic seismic detection technique.OilGeophysicalProspecting（in Chinese），2011，46（1）：121－125，doi：1000－7210（2011）01－0121－125.

［6］ 唐向宏，賀振華.圖像邊緣檢測方法在薄地層識別中的應(yīng)用.石油地球物理勘探，2002，37（2）：163－165.Tang X H，He Z H.Study for application of image edges detection to thin－layer recognition.OilGeophysicalProspecting（in Chinese），2002，37（2）：163－165.

［7］ 李紅星，劉財，陶春輝.圖像邊緣檢測方法在地震剖面同相軸自動檢測中的應(yīng)用研究.地球物理學(xué)進(jìn)展，2007，22（5）：1607－1610.Li H X，Liu C，Tao C H .The study of application of edge measuring technique to the detection of phase axis of the seismic section.ProgressinGeophysics（in Chinese），2007，22（5）：1607－1610，doi：1004－2903（2007）05－1607－04.

［8］ 茍量，彭真明.小波多尺度邊緣檢測及其在裂縫預(yù)測中的應(yīng)用.石油地球物理勘探，2005，40（3）：309－313.Gou L，Peng Z M.Multi－scale edge detection of wavelet and application in fracture prediction.OilGeophysicalProspecting（in Chinese），2005，40（3）：309－313.

［9］ 李輝峰，鄒強(qiáng)，金文昱.基于邊緣檢測的初至波自動拾取方法.石油地球物理勘探，2006，41（2）：150－155.Li H F，Zou Q，Jin W Y.Method of automatic first breaks pick－up based on edge detection.OilGeophysicalProspecting（in Chinese），2006，41（2）：150－155.

［10］ 許輝群，桂志先.邊緣檢測技術(shù)砂體邊界識別方法研究.石油天然氣學(xué)報，2009，31（5）：75－77.Xu H Q，Gui Z X.Method to identify sandstone boundary with edge detection.JournalofOilandGasTechnology（J.JPI）（in Chinese），2009，31（5）：75－77.Doi：1000－9752（2009）05－0075－03.

［11］ 熊會軍，管業(yè)鵬，于蘊(yùn)杰等.基于圖像邊緣檢測方法提取地震剖面同相軸.地球物理學(xué)進(jìn)展，2009，24（6）：2250－2254.Xiong H J，Guan Y P，Yu Y J，et al.Extraction of cophasal axes on seismic sections based on the edge detection method.ProgressinGeophysics（in Chinese），2009，24（6）：2250－2254.doi：10.3969／j.issn.1004－2903.2009.06.045.

［12］ Luo Y，Higgs W G，Kowalik W S.Edge detection and stratigraphic analysis using 3Dseismic data.SEGTechnical ProgramExpandedAbstracts，1996，15（1）：324－327.

［13］ Aqrawi A A，Boe T H，Barros S.Detecting salt domes using a dip guided 3DSobel seismic attribute.SEGTechnical ProgramExpandedAbstracts，2011，30（1）：1014－1018.

［14］ Aqrawi A A，Boe T H.Improved fault segmentation using a dip guided and modified 3DSobel filter.SEGTechnical ProgramExpandedAbstracts，2011，30（1）：999－1003.

［15］ Weickert J.Coherence－enhancing diffusion filtering.International JournalofComputerVision，1999，31（2－3）：111－127.

［16］ Haralick R M.Digital step edges from zero crossing of second directional derivatives.IEEETransactionsonPatternAnalysis andMachineIntelligence，1984，6（1）：58－68.Doi：0162－8828／84／0100－005801.00.

［17］ 王緒松，楊長春.對地震圖像進(jìn)行保邊濾波的非線性各向異性擴(kuò)散算法.地球物理學(xué)進(jìn)展，2006，21（2）：452－457.Wang X S，Yang C C.An edge－preserving smoothing algorithm of seismic image using nonlinear anisotropic diffusion equation.ProgressinGeophysics（in Chinese），2006，21（2）：452－457，doi：1004－2903（2006）02－0452－06.

［18］ Georgia R，高靜懷，汪玲玲等.基于方向模板的保持邊緣平滑后地震圖像的不連續(xù)性檢測.地球物理學(xué)進(jìn)展，2008，23（5）：1398－1405.Georgia R，Gao J H，Wang L L，et al.Discontinuities detection on edge preserved smoothed seismic image using directional masks.ProgressinGeophysics（in Chinese），2008，23（5）：1398－1405，doi：1004－2903（2008）05－1398－08.

［19］ 陸艷洪，陸文凱，翟正軍.一種邊緣保持的地震數(shù)據(jù)插值方法.地球物理學(xué)報，2012，55（3）：991－997.Lu Y H，Lu W K，Zhai Z J.An edge－preserving seismic data interpolation method.ChineseJournalofGeophysics（in Chinese），2012，55（3）：991－997，doi：10.6038／j.issn.0001－5733.2012.03.029.