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

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

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

0 引言

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

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

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

1 方法原理

1.1 各向異性擴(kuò)散濾波

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

(1)

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

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

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

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

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

(2)

1.2 熵的引入

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

(3)

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

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

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

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

(4)

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

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

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

(5)

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

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

(6)

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

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

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

Sa=Ka*(Uσ?Uσ)

(7)

式中：Uσ=Kσ*U，其中Kσ為高斯核函數(shù)，“*”為褶積運(yùn)算符，σ為噪聲尺度；“?”表示矩陣轉(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為一個(gè)接近于0的正數(shù)，代表擴(kuò)散強(qiáng)度；C通常設(shè)置為1；k為一致性參數(shù)，有

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

(12)

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

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

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

(13)

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

1.4 方法的實(shí)現(xiàn)

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

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

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

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

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

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

(6)如果n

2 理論模型測試

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

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

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

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

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

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

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

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

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

4 結(jié)論

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