圖像熵各向異性擴散保邊濾波方法及在斷層識別中的應(yīng)用

李 軍 張軍華* 劉 楊 楊 勇 杜玉山

(①中國石油大學(華東)地球科學與技術(shù)學院，山東青島 266580； ②中國石化勝利油田分公司勘探開發(fā)研究院，山東東營 257015)

0 引言

斷層的精細識別與描述對于斷塊型油氣藏的勘探開發(fā)具有重要意義。斷塊型油氣藏斷層復(fù)雜，人工解釋斷層的工作量與難度都較大，需借助相干、邊緣檢測、曲率等[1]地震屬性技術(shù)或者不同手段結(jié)合的綜合解釋方法[2]進行斷層的追蹤與識別。斷層解釋精度與地震資料信噪比有關(guān)，一般的地震數(shù)據(jù)受到地質(zhì)條件、采集條件和資料處理等多個因素的影響，信噪比通常較低，特別是中深層資料。為了提高地震資料的信噪比，同時保護斷層等地質(zhì)體邊緣信息，需要對地震數(shù)據(jù)做具有保邊作用的去噪處理。

使用雙邊濾波、Radon變換、F-X濾波等[3-5]傳統(tǒng)濾波方法去噪，并不能有效保護斷層等地質(zhì)體的邊界及方向信息。近年來，考慮地質(zhì)體構(gòu)造和裂縫發(fā)育方向的保邊濾波技術(shù)得到了很大的發(fā)展。趙明章等[6]利用構(gòu)造導(dǎo)向濾波技術(shù)提高地震資料信噪比，落實復(fù)雜斷塊圈閉；尹川等[7]利用基于傾角控制的構(gòu)造導(dǎo)向濾波技術(shù)在斷層識別中取得了很好的效果；徐紅霞等[8]基于構(gòu)造導(dǎo)向濾波技術(shù)將多屬性分析技術(shù)應(yīng)用在斷溶體預(yù)測中，提高了預(yù)測精度和實際鉆井吻合率；問雪等[9]利用構(gòu)造導(dǎo)向平滑方法進行了斷層解釋工作，提高了解釋效率及精度；Fehmers等[10]首次將各向異性擴散濾波技術(shù)應(yīng)用于地震數(shù)據(jù)去噪，在保護斷層等有用的構(gòu)造信息的同時，提高了數(shù)據(jù)的信噪比；Weickert[11]和Lavialle等[12]對一致性增強濾波進行了改進，能更好地兼顧去噪與保護斷層等構(gòu)造邊界信息；楊培杰等[13]提出了一種方向性邊界保護濾波方法，在提高信噪比同時，較好地保持了斷層信息；孫夕平等[14]討論了如何選取各向異性擴散濾波的參數(shù)，提出了具有旋轉(zhuǎn)不變特性的改進方法；王旭松等[15]結(jié)合傾角信息對各向異性擴散濾波構(gòu)造張量矩陣進行了簡化，提高了效率；張爾華等[16]將各向異性擴散濾波方法拓展到三維，在三維實際數(shù)據(jù)去噪中取得了較好的效果；嚴哲等[17]在各向異性擴散濾波算法中將相干值作為一個參數(shù)引入具體算法中，提高了對橫向不連續(xù)點的保護性能；楊千里等[18]和姚振岸等[19-20]應(yīng)用各向異性擴散濾波進行去噪，斷層形態(tài)更清晰、準確。

本文在前人研究基礎(chǔ)之上，結(jié)合圖像處理中熵的概念及特征[21-24]，提出了一種基于熵的各向異性擴散保邊濾波方法。在地震數(shù)據(jù)連續(xù)性較好的區(qū)域，熵值較大；而在包含斷點等特殊地質(zhì)構(gòu)造區(qū)域，其熵值較小。根據(jù)熵的特征，其大小只與數(shù)據(jù)的總體結(jié)構(gòu)有關(guān)，與具體的振幅值無關(guān)。在計算數(shù)據(jù)的各向異性擴散幅度時，本文方法能在剔除噪聲異常點影響情況下，定量控制保護特殊地質(zhì)邊界的二階導(dǎo)數(shù)所占比例，從而更準確地保持斷層等邊界信息。

1 方法原理

1.1 各向異性擴散濾波

各向異性擴散濾波的擴散過程等價于物理學熱傳導(dǎo)過程，類似于熱傳導(dǎo)方程，可以將擴散過程統(tǒng)一表示為[24-25]

(1)

式中：U表示地震數(shù)據(jù)的振幅信息；c為擴散系數(shù)；t為擴散時間。

(1)當系數(shù)c為常數(shù)時，擴散方程是線性各向同性方程，等同于對數(shù)據(jù)做高斯濾波，能夠提高信噪比，但不能保護好斷層等地質(zhì)體的邊界及方向信息；

(2)當系數(shù)c使用圖像遞減函數(shù)[26]等改進系數(shù)來代替，擴散方程就是非線性各向同性方程，這種情況下去噪易造成階梯效應(yīng)和針孔效應(yīng)，且對強噪?yún)^(qū)域無效；

(3)使用擴散張量矩陣D替代擴散系數(shù)，則是真正意義上的各向異性擴散，各向異性信息隱藏在張量D中，在一些方向允許擴散，在其他方向如斷層等地質(zhì)體邊界區(qū)域減少或者不允許擴散，既提高了信噪比，又保護了邊界信息。

因此，三維各向異性擴散濾波方程最終可以表示為

(2)

1.2 熵的引入

熵的概念由Clausius首先給出，用來評價能量在空間中的分布均勻程度，熵值越大說明能量分布越均勻[24]。當整個空間的能量完全均勻時，熵值達到最大?！皥D像熵”是對圖像特征的一種統(tǒng)計表示，能夠評價圖像的平均信息量。相比于方差、均值等統(tǒng)計量，熵只與圖像的總體結(jié)構(gòu)有關(guān)，與具體的像素灰度值大小無關(guān)，可以避免噪聲變化對圖像梯度值的影響，從而可以更好地反映圖像的變化。圖像熵定義為

(3)

式中pi表示圖像中灰度值為i的像素所占比例。

類似于圖像紋理，地震數(shù)據(jù)具有同樣的紋理結(jié)構(gòu)。引入熵的概念，在不受噪聲干擾的情況下評價區(qū)域內(nèi)振幅的隱含信息：當熵值較大時，說明該區(qū)域為常規(guī)地層；當熵值較小時，說明該區(qū)域振幅變化大，應(yīng)包含著斷層等特殊地質(zhì)體的結(jié)構(gòu)信息。計算地震數(shù)據(jù)熵的步驟如下。

(1)將地震數(shù)據(jù)灰度值化。

(2)以目標點為中心，選擇3×3×3大小的窗口構(gòu)建一個灰度值數(shù)據(jù)子體，而后利用熵的概念，計算該子體的熵值作為該目標點的熵值，然后依次對每個點逐次掃描計算其熵值

(4)

式中：(i，j，k)表示目標點的坐標信息；pi，j，k，t表示灰度數(shù)據(jù)子體中灰度值為t的像素所占比例；ηs表示以目標點為中心構(gòu)建的窗口內(nèi)灰度值級數(shù)大小。

(3)對計算得到的所有數(shù)據(jù)點的熵值大小進行無量綱[0，1]區(qū)間的歸一化處理，得到歸一化后的熵值數(shù)據(jù)H1。

(4)計算地震數(shù)據(jù)的熵信息閾值

(5)

式中：Nx為Inline方向的道數(shù)；Ny為Crossline方向的道數(shù)；Nt為時間方向的采樣點數(shù)。

一般在斷層等特殊地質(zhì)體邊界處，其振幅二階導(dǎo)數(shù)具有局部極大值，其他平坦區(qū)域，二階導(dǎo)數(shù)值很小，與熵值呈負相關(guān)的關(guān)系。因此，在計算具有區(qū)域方向信息的結(jié)構(gòu)張量矩陣并添加二階導(dǎo)數(shù)輔助保護斷層等特殊地質(zhì)體邊界信息時，可將改造后的熵值作為二階導(dǎo)數(shù)添加量的比例系數(shù)

(6)

當計算點的熵值大于閾值時，認為該區(qū)域近似平層，不需要添加二階導(dǎo)數(shù)信息。為了后續(xù)各向異性擴散濾波迭代的穩(wěn)定性，閾值設(shè)置不宜較大。

1.3 擴散張量的構(gòu)建

各向異性擴散濾波的關(guān)鍵就是構(gòu)建擴散張量D。這就需要計算局部方向信息，然后依此作為特征向量構(gòu)建擴散張量D。一般使用具有對稱半正定的結(jié)構(gòu)張量計算方向信息，即

Sa=Ka*(Uσ?Uσ)

(7)

式中：Uσ=Kσ*U，其中Kσ為高斯核函數(shù)，“*”為褶積運算符，σ為噪聲尺度；“?”表示矩陣轉(zhuǎn)置與矩陣相乘；Kα為高斯核函數(shù)，其中α為整合尺度，一般應(yīng)大于噪聲尺度；Uσ為數(shù)據(jù)的梯度信息，即

(8)

添加二階導(dǎo)數(shù)信息后，結(jié)構(gòu)張量重新寫為

Sa=Ka*[Uσ?Uσ+

(9)

(10)

式中：ε代表橫向不連續(xù)因子，對于連續(xù)地層ε→1，在斷層位置ε→0，本文選取歸一化相干值作為連續(xù)因子；β1、β2、β3為D的特征值，可以寫為

(11)

式中：b為一個接近于0的正數(shù)，代表擴散強度；C通常設(shè)置為1；k為一致性參數(shù)，有

k=(β1-β2)2+(β1-β3)2+(β2-β3)2

(12)

在一致性較高區(qū)域k遠大于C，β2=β3≈1，擴散模型沿ξ2、ξ3方向擴散；在各向同性區(qū)域(k近似于0)，三個方向擴散強度都很小。

得到擴散張量D后， 就可以通過迭代法求解式(2)，最終得到擴散后的數(shù)據(jù)

Un+1=Un+Δt·div(D·Un)

(13)

式中：Un和Un+1分別為第n次和第n+1次迭代得到的數(shù)據(jù)； Δt為迭代步長。

1.4 方法的實現(xiàn)

基于熵的三維各向異性擴散濾波的具體實現(xiàn)過程為：

(1)定義迭代次數(shù)為n=0，并設(shè)置最大迭代次數(shù)為N，此時，原始數(shù)據(jù)U0(x，y，z)為迭代初始值；

(2)使用噪聲尺度為σ的高斯核函數(shù)對數(shù)據(jù)進行預(yù)處理；

(3)計算熵值，重構(gòu)結(jié)構(gòu)張量矩陣，求取矩陣的特征值及對應(yīng)的特征向量；

(4)計算擴散張量的特征值，確定橫向連續(xù)性因子ε；

(5)構(gòu)建擴散張量矩陣D，根據(jù)式(13)計算第1次迭代結(jié)果，更新迭代次數(shù)n；

(6)如果n

2 理論模型測試

Marmousi模型是一個包含復(fù)雜構(gòu)造特征的經(jīng)典地質(zhì)模型，本文截取其中包含斷層構(gòu)造的部分數(shù)據(jù)并添加了一些隨機噪聲后的模型(圖1a，信噪比為4)進行測試。

為了對比分析不同方法去噪保邊的效果，分別采用雙邊濾波、Kuwaharab濾波與本文方法對模型進行處理，結(jié)果如圖1b～圖1d所示?？梢钥闯觯弘p邊濾波方法雖在一定程度上消除了噪聲，但平滑作用較強，邊緣保持不理想；Kuwaharab濾波雖然能保持邊緣信息，但去噪效果比雙邊濾波差；而本文方法消除噪聲效果良好，同時較好保護了邊緣特征，濾波后地質(zhì)邊界清晰，同相軸連續(xù)性也得到了強化，說明了本文方法保邊去噪的有效性。

圖1 含噪Marmousi模型及不同方法去噪結(jié)果對比

3 實際資料應(yīng)用

選取斷裂眾多、結(jié)構(gòu)復(fù)雜、勘探難度較大的勝利油田M斷塊油藏實際資料驗證本文方法的應(yīng)用效果。圖2為應(yīng)用基于圖像熵的各向異性擴散保邊濾波對地震數(shù)據(jù)進行去噪處理前、 后Inline4080剖面，可見： 濾波后同相軸更連續(xù)，斷點更清晰(黑色箭頭所示)，斷層更易識別，說明本文方法在提高資料信噪比的同時也保護了斷層等地質(zhì)體的邊緣信息。

對應(yīng)用本文濾波方法前、后的地震數(shù)據(jù)體分別計算最大正曲率屬性。圖3為1.55s時間切片，圖4為沿T4的層位屬性切片。可見，無論是時間切片還是沿層切片，濾波后都提高了分辨率，斷層的形態(tài)及展布更符合地質(zhì)規(guī)律，濾波前一些不易識別的斷層(紅色箭頭所示)在濾波后其形態(tài)及展布更清晰，對后續(xù)的油藏開發(fā)具有重要意義。

圖2 本文方法濾波前(a)、后(b)Inline4080剖面對比

圖3 濾波前(a)、后(b)曲率屬性1.55s時間切片對比

圖4 濾波前(a)、后(b)沿T4層曲率屬性切片對比

4 結(jié)論

提高地震數(shù)據(jù)信噪比的同時保持斷層等地質(zhì)體的邊界信息，對于后續(xù)的地震精細解釋具有非常重要的意義。根據(jù)數(shù)據(jù)的熵信息，在構(gòu)建攜帶方向信息的結(jié)構(gòu)張量矩陣時，添加合適比例的二階導(dǎo)數(shù)信息，在各向異性擴散濾波過程中，能更好地保護邊界信息。模型測試及實際地震資料的應(yīng)用結(jié)果表明，濾波后的地震資料信噪比明顯提高，在連續(xù)地層區(qū)域，地震同相軸連續(xù)性增強，斷層等地質(zhì)體邊界得到有效保護，斷點更清晰。但由于需要計算地震數(shù)據(jù)熵信息，導(dǎo)致計算量增大，耗時較多，且對一些能量較弱區(qū)域去噪會造成細節(jié)丟失。因此，如何提高計算效率且不造成細節(jié)丟失是需要重點改進的方向。