利用不連續(xù)性的各向異性擴(kuò)散濾波方法識(shí)別斷層

王 靜 張軍華 馮德永 李紅梅

(①中國(guó)石油大學(xué)(華東)地球科學(xué)與技術(shù)學(xué)院，山東青島 266580； ②中國(guó)石化勝利油田股份公司物探研究院，山東東營(yíng) 257022)

0 引言

通常斷塊型油藏的地質(zhì)特征較復(fù)雜，其地震響應(yīng)經(jīng)過(guò)處理后仍存在噪聲，特別是斷層附近的地震反射數(shù)據(jù)的噪聲能量更強(qiáng)。當(dāng)前的絕大多數(shù)斷層檢測(cè)方法主要利用斷層處的地震反射信息急速變化識(shí)別與描述斷層，對(duì)噪聲較敏感，因此需要預(yù)先壓制地震噪聲，同時(shí)保護(hù)斷層等的邊緣信息。

采用傳統(tǒng)的濾波方法(中值濾波、F-X濾波、高斯濾波等)[1-2]去除地震噪聲，不能有效保護(hù)斷層等不連續(xù)的邊界信息，因而人們廣泛研究了保邊濾波技術(shù)。目前，各向異性擴(kuò)散濾波是保邊去噪處理中應(yīng)用較廣泛的一種濾波技術(shù)，其擴(kuò)散方程的形式為偏微分方程，具有較高的計(jì)算效率和穩(wěn)定性。Perona等[3]在各向同性擴(kuò)散方程的基礎(chǔ)上提出了各向異性擴(kuò)散方程，使擴(kuò)散濾波在不同方向上具有不同的擴(kuò)散速度，被稱(chēng)為Perona-Malik(簡(jiǎn)稱(chēng)P-M)擴(kuò)散模型； Weickert[4]系統(tǒng)研究了各向異性擴(kuò)散理論，并結(jié)合相干體技術(shù)提出了相干增強(qiáng)擴(kuò)散算法，模型試算及實(shí)際數(shù)據(jù)應(yīng)用都取得了較好的效果； Fehmers等[5]首次將各向異性擴(kuò)散濾波技術(shù)用于地震勘探領(lǐng)域，在保護(hù)斷層等構(gòu)造信息的同時(shí)，壓制了噪聲； Hale[6-7]提出了基于結(jié)構(gòu)導(dǎo)向的各向異性擴(kuò)散濾波，可以有效提高資料信噪比；問(wèn)雪等[8]將構(gòu)造導(dǎo)向平滑方法應(yīng)用于斷層解釋?zhuān)岣吡藬鄬咏忉尵龋?Lavialle等[9]改進(jìn)了三維擴(kuò)散張量特征值的設(shè)計(jì)方法，提出了保護(hù)斷層信息的擴(kuò)散濾波方法，使斷層信息更突出； Chopra等[10]利用各向異性擴(kuò)散濾波對(duì)地震數(shù)據(jù)預(yù)處理，取得了很好的效果； 楊培杰等[11]提出了一種方向性邊界信息保護(hù)的濾波方法，在保護(hù)斷層信息的同時(shí)，提高了資料信噪比； 孫夕平等[12]提出了一種具有最優(yōu)旋轉(zhuǎn)不變性的各向異性擴(kuò)散濾波方法，可有效提高地震資料信噪比，突出地層接觸關(guān)系，增強(qiáng)地震數(shù)據(jù)對(duì)層序體內(nèi)部結(jié)構(gòu)的成像能力；張爾華等[13]將各向異性擴(kuò)散濾波方法拓展到三維；嚴(yán)哲等[14]在各向異性擴(kuò)散濾波算法中將相干值作為一個(gè)參數(shù)引入具體算法，提高了對(duì)橫向不連續(xù)點(diǎn)的保護(hù)性能； 楊千里等[15]、姚振岸等[16-17]應(yīng)用各向異性擴(kuò)散濾波對(duì)地震數(shù)據(jù)進(jìn)行保邊去噪，可以清晰、準(zhǔn)確地識(shí)別斷層； 李福強(qiáng)等[18]提出了基于斷層算子的各向異性擴(kuò)散濾波方法，可有效去除噪聲，自適應(yīng)增強(qiáng)斷點(diǎn)處的能量，提高了斷層識(shí)別精度； 李軍等[19]提出了圖像熵各向異性擴(kuò)散濾波方法，可定量給出構(gòu)建結(jié)構(gòu)張量矩陣需要添加的二階導(dǎo)數(shù)占比，更合理地判斷擴(kuò)散方向。

本文在前人的研究基礎(chǔ)上，分析了結(jié)構(gòu)張量特征值與三維地震圖像的局部結(jié)構(gòu)特征之間的關(guān)系，利用斷層置信度參數(shù)重新構(gòu)建擴(kuò)散張量特征值，進(jìn)而重構(gòu)擴(kuò)散張量，實(shí)現(xiàn)基于不連續(xù)性的三維各向異性擴(kuò)散濾波方法，可有效提高地震資料信噪比，保護(hù)斷層等不連續(xù)的邊界信息。

1 方法原理

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

各向異性擴(kuò)散濾波的擴(kuò)散過(guò)程等價(jià)于熱傳導(dǎo)過(guò)程，表示為[20]

(1)

式中：U為輸入的地震數(shù)據(jù)振幅信息；t為擴(kuò)散時(shí)間；g為擴(kuò)散系數(shù)。g的取值有三種情況： ①當(dāng)g取常數(shù)時(shí)，式(1)為線(xiàn)性各向同性方程，對(duì)數(shù)據(jù)進(jìn)行濾波后，能夠提高信噪比，但不能保護(hù)斷層等地質(zhì)體的邊界及方向信息； ②當(dāng)g為非負(fù)單調(diào)遞減函數(shù)[3]時(shí)，式(1)為非線(xiàn)性各向同性方程，濾波后容易產(chǎn)生階梯效應(yīng)和針孔效應(yīng)，而在強(qiáng)噪聲區(qū)域無(wú)去噪效果； ③使用擴(kuò)散張量D作為擴(kuò)散系數(shù)[4]，可以根據(jù)各向異性信息對(duì)數(shù)據(jù)進(jìn)行自適應(yīng)濾波，其三維各向異性擴(kuò)散濾波方程為

(2)

式中U0(x,y,z)為原始三維地震數(shù)據(jù)，(x，y，z)為某點(diǎn)的空間坐標(biāo)。

1.2 結(jié)構(gòu)張量的構(gòu)建

在構(gòu)建D的過(guò)程中，一般采用具有對(duì)稱(chēng)半正定的結(jié)構(gòu)張量Sρ計(jì)算局部方向信息，其形式為

(3)

(4)

式中(x0,y0,z0)為參考點(diǎn)空間坐標(biāo)。在斷層位置處經(jīng)過(guò)高斯濾波后地震數(shù)據(jù)對(duì)空間的二階導(dǎo)數(shù)一般具有局部極大值，將該信息添加到Sρ的計(jì)算公式中，可將式(3)改寫(xiě)為

(5)

式中α∈[0,1]為穩(wěn)定系數(shù)，用于控制高斯濾波后地震數(shù)據(jù)對(duì)空間的二階導(dǎo)數(shù)的貢獻(xiàn)率。

對(duì)Sρ進(jìn)行特征分解，得

(6)

式中：λ1、λ2和λ3為Sρ的三個(gè)特征值，且λ1>λ2>λ3；ζ1、ζ2和ζ3分別為λ1、λ2和λ3對(duì)應(yīng)的特征向量，ζ1指示梯度變化最大的方向，ζ3指示梯度變化最小的方向；ε為橫向不連續(xù)因子，在地層連續(xù)處ε趨近于1，在斷層位置處ε趨近于0。

Bakker[21]指出，對(duì)于三維地震數(shù)據(jù)來(lái)說(shuō)，當(dāng)結(jié)構(gòu)張量的特征值λi(i=1，2，3)為0時(shí)，表示圖像沿該特征向量的方向移不變。當(dāng)沿一個(gè)方向移不變時(shí)，為線(xiàn)狀線(xiàn)性結(jié)構(gòu)；當(dāng)沿兩個(gè)方向移不變時(shí)，為面狀線(xiàn)性結(jié)構(gòu)。根據(jù)特征值的大小及其對(duì)應(yīng)的不同的特征向量可以將地震圖像結(jié)構(gòu)表述為以下四類(lèi)：

(1)λ1≈λ2≈λ3≈0，三個(gè)特征向量方向的變化很小或者基本沒(méi)有變化，指示平滑區(qū)域，表示地質(zhì)體呈各向同性；

(2)λ1>λ2≈λ3≈0，在ζ1方向梯度變化較大，另外兩個(gè)特征方向的梯度變化很小且?guī)缀跸嗟龋甘酒矫婢€(xiàn)性地質(zhì)體，對(duì)應(yīng)模型為近似面狀線(xiàn)性結(jié)構(gòu)；

(3)λ1>λ2>λ3≈0，其中一個(gè)方向梯度變化近似為零，指示在斷層等不連續(xù)的邊界處出現(xiàn)不連續(xù)的間斷信息，對(duì)應(yīng)模型為線(xiàn)狀結(jié)構(gòu)模型；

(4)λ1>λ2>λ3>0，三個(gè)特征向量方向的梯度變化都不同，并且均不為零，說(shuō)明區(qū)域各向異性特征明顯，表示地質(zhì)體的不整合結(jié)構(gòu)。

van Kempen等[22]根據(jù)結(jié)構(gòu)張量特征值與地震圖像結(jié)構(gòu)的關(guān)系，用線(xiàn)狀置信度量Mline和面狀置信度量Mplane表征圖像結(jié)構(gòu)[18]，即

(7)

Mline∈[0,1]、Mplane∈[0,1]。Mline描述某點(diǎn)鄰域與線(xiàn)狀結(jié)構(gòu)的相似程度，當(dāng)Mline=1時(shí)，指示線(xiàn)狀結(jié)構(gòu)；Mplane描述某點(diǎn)鄰域與面狀結(jié)構(gòu)的相似程度，當(dāng)Mplane=1時(shí)，指示面狀結(jié)構(gòu)。由Mline、Mplane得到檢測(cè)斷層等不連續(xù)結(jié)構(gòu)邊界信息的斷層置信度量Mfault[21,23]

(8)

表1為置信度參數(shù)與地層結(jié)構(gòu)以及特征值之間的關(guān)系。

表1 置信度參數(shù)與地層結(jié)構(gòu)以及特征值之間的關(guān)系

在地震剖面上，如果地震同相軸的連續(xù)性較好，則Mline趨于0、Mplane趨于1、Mfault趨于0；同理，在存在斷層等構(gòu)造的同相軸不連續(xù)區(qū)域，Mline趨于1、Mplane趨于0、Mfault趨于1。

1.3 利用Mfault設(shè)計(jì)擴(kuò)散張量特征值

構(gòu)建擴(kuò)散張量D時(shí)，為了保證擴(kuò)散方向沿著構(gòu)造方向，D的特征向量應(yīng)與Sρ的一致，因此

(9)

式中γ1、γ2和γ3為D的三個(gè)特征值，對(duì)應(yīng)的特征向量也為ζ1、ζ2和ζ3。

在三維各向異性擴(kuò)散濾波中，擴(kuò)散濾波的強(qiáng)度取決于γ1、γ2和γ3。特征值越接近0，擴(kuò)散濾波的強(qiáng)度越??；特征值趨于1，擴(kuò)散濾波強(qiáng)度最大。

為了克服常規(guī)擴(kuò)散濾波方法不能保護(hù)數(shù)據(jù)邊緣信息的缺點(diǎn)，在D的特征值設(shè)計(jì)中引入Mline、Mplane以及Mfault，即

(10)

式(10)表明： ①沿梯度變化最大的方向，即沿ζ1方向，γ1=a，該方向的擴(kuò)散系數(shù)趨于0，即擴(kuò)散濾波的強(qiáng)度較小，可以保護(hù)斷層等不連續(xù)結(jié)構(gòu)的邊界信息； ②沿ζ2和ζ3的方向，即垂直于最大梯度的方向，可以根據(jù)Mline、Mplane以及Mfault調(diào)整特征值。如果Mfault趨于0，地層為面狀結(jié)構(gòu)，γ2和γ3的值趨于1，沿ζ2和ζ3組成的平面擴(kuò)散模型濾波強(qiáng)度大； 如果Mfault趨于1，γ1≈γ2≈a，沿ζ1、ζ2方向擴(kuò)散濾波強(qiáng)度小，可以有效保護(hù)斷層等不連續(xù)的邊界信息。

綜上所述，利用Mfault設(shè)計(jì)D特征值的方法可在增強(qiáng)擴(kuò)散濾波作用的同時(shí)，較好地保護(hù)斷層等特殊地質(zhì)結(jié)構(gòu)的邊界信息。

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

計(jì)算出D后，通過(guò)迭代法求解三維各向異性擴(kuò)散方程，最終得到擴(kuò)散濾波后的地震數(shù)據(jù)

(11)

式中：U(k)和U(k+1)分別為第k和k+1次迭代計(jì)算得到的數(shù)據(jù)；Δt為迭代步長(zhǎng)。在實(shí)際應(yīng)用中為了保證迭代過(guò)程的收斂以及穩(wěn)定性，Δt要足夠小，但是若Δt過(guò)小，會(huì)增加迭代次數(shù)，從而增加計(jì)算時(shí)間。通常選擇0.05≤Δt≤1即可滿(mǎn)足地震數(shù)據(jù)處理要求。

本文方法的實(shí)現(xiàn)流程為：

(1)定義初始迭代次數(shù)k為0，用原始三維數(shù)據(jù)U0(x,y,z)作為迭代初始值并設(shè)置最大迭代次數(shù)為K；

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

(3)重構(gòu)Sρ并求取特征值λ1、λ2和λ3以及對(duì)應(yīng)的特征向量ζ1、ζ2和ζ3；

(4)利用Mfault設(shè)計(jì)D的特征值，并確定橫向連續(xù)性因子ε；

(5)構(gòu)建D后，根據(jù)式(11)計(jì)算第k+1次迭代結(jié)果；

(6)如果k

2 理論模型測(cè)試

為驗(yàn)證基于Mfault的各向異性擴(kuò)散濾波方法的應(yīng)用效果，選用三維Qdome模型(圖1)[26]進(jìn)行驗(yàn)證。根據(jù)Qdome原始地震數(shù)據(jù)(圖2a)的頻帶對(duì)隨機(jī)噪聲進(jìn)行低通濾波，然后將其加入到理論模型中得到加噪模型(圖2b)。

圖1 Qdome模型立體顯示

通過(guò)前文分析可知，各向異性擴(kuò)散濾波涉及的參數(shù)較多，主要包括高斯核函數(shù)的噪聲尺度σ、整合尺度ρ、迭代步長(zhǎng)Δt以及最大迭代次數(shù)K。為了研究各個(gè)參數(shù)對(duì)濾波效果的影響，利用圖2b對(duì)參數(shù)選擇進(jìn)行試驗(yàn)。

圖2 加噪前(a)、后(b)Inline70剖面

2.1 濾波參數(shù)的選擇

2.1.1σ的影響

圖3展示了不同噪聲尺度的濾波結(jié)果。由圖可見(jiàn)： ①在其他濾波參數(shù)固定的情況下，隨σ取值增大，信噪比隨之變大(圖3a、圖3b)，但當(dāng)信噪比達(dá)到最大值(SNR=10.8947)后，隨σ取值增大，信噪比開(kāi)始下降(圖3c)； ②σ取值過(guò)大，斷層邊界模糊(圖3d箭頭處)。因此當(dāng)σ∈[1.0，2.0]時(shí)濾波效果較好。

圖3 不同噪聲尺度的濾波結(jié)果

2.1.2ρ的影響

圖4為不同整合尺度的濾波結(jié)果。由圖可見(jiàn)：在其他濾波參數(shù)固定的情況下，隨ρ取值增大，信噪比隨之變大(圖4a、圖4b)，但當(dāng)信噪比達(dá)到最大值(SNR=10.8947)后，開(kāi)始出現(xiàn)不穩(wěn)定變化(圖4c、圖4d)。因此，ρ取值不能太大。

圖4 不同整合尺度的濾波結(jié)果

2.1.3 Δt的影響

圖5為不同迭代步長(zhǎng)的濾波結(jié)果。由圖可見(jiàn)：Δt過(guò)小，導(dǎo)致計(jì)算時(shí)間增加； Δt過(guò)大，資料信噪比并未明顯提高，而且造成地震數(shù)據(jù)的能量損失，斷層邊緣處出現(xiàn)模糊(圖5d紅色箭頭處)。因此，當(dāng)Δt∈[0.05，1.00]時(shí)濾波效果較好，且計(jì)算效率較高。

圖5 不同迭代步長(zhǎng)的濾波結(jié)果

2.1.4K的影響

圖6為不同最大迭代次數(shù)的濾波結(jié)果。由圖可見(jiàn)，K越大，計(jì)算時(shí)間越長(zhǎng)，且會(huì)破壞斷層的邊界特征(圖6d紅色箭頭處)。因此，當(dāng)K∈[20，50]時(shí)，可明顯提高信噪比，且斷層邊緣特征保持較好(圖6b、圖6c)。若地震資料信噪比太低，可以增大K，以取得較好的濾波效果。

圖6 不同最大迭代次數(shù)的濾波結(jié)果

2.2 模型試算

通過(guò)參數(shù)實(shí)驗(yàn)確定了濾波參數(shù)的選擇范圍，分別應(yīng)用本文方法和各向異性擴(kuò)散濾波方法[15]對(duì)圖2b進(jìn)行去噪處理，結(jié)果表明：本文方法濾波后數(shù)據(jù)(圖7a左)的信噪比高于各向異性擴(kuò)散濾波方法(圖7b左)；本文方法濾波后地震數(shù)據(jù)基本沒(méi)有能量損失(圖7a右)，各向異性擴(kuò)散濾波方法濾波后地震數(shù)據(jù)在斷層邊緣處有輕微能量損失(圖7b右紅色箭頭處)。圖8為本文方法濾波前、后數(shù)據(jù)頻譜。由圖可見(jiàn)，濾波處理基本不影響地震數(shù)據(jù)頻帶。模型試算說(shuō)明本文方法在去噪和保護(hù)地質(zhì)體邊緣信息方面的有效性。

圖8 本文方法濾波前(a)、后(b)數(shù)據(jù)頻譜

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

選取斷裂眾多、結(jié)構(gòu)復(fù)雜、勘探難度較大的勝利油田X斷塊地震資料作為研究數(shù)據(jù)，應(yīng)用文中方法進(jìn)行濾波處理并計(jì)算相干屬性。圖9為本文方法濾波前、后的Inline3570剖面。由圖可見(jiàn)，濾波后地震剖面的同相軸更連續(xù)、斷點(diǎn)更清晰(圖9b黃色和藍(lán)色箭頭處)，斷層更易識(shí)別，說(shuō)明本文方法在去噪的同時(shí)能有效保護(hù)斷層等不連續(xù)性邊界的信息。圖10為濾波前、后的相干時(shí)間切片。由圖可見(jiàn)，利用濾波后的相干時(shí)間切片(圖10b)可準(zhǔn)確地描述斷層形態(tài)及展布，如在濾波前的相干時(shí)間切片中一些不易識(shí)別的小斷層(圖10a中黃色和紅色箭頭處)，在濾波后的相干時(shí)間切片上得到準(zhǔn)確、清晰地識(shí)別，這也證明了該濾波方法在提高資料信噪比和保持?jǐn)鄬拥冗吘壭畔⒎矫娴挠行院秃蜏?zhǔn)確性。

圖9 本文方法濾波前(a)、后(b)的Inline3570剖面

圖10 濾波前(a)、后(b)的相干時(shí)間切片(1370ms)

4 結(jié)論

針對(duì)三維地震資料保邊濾波處理的需要，本文研究了三維各向異性擴(kuò)散濾波方法，并采用結(jié)構(gòu)張量分析了三維地質(zhì)結(jié)構(gòu)的局部特征。基于斷層置信度參數(shù)重新設(shè)計(jì)了擴(kuò)散張量的特征值，可以控制地震數(shù)據(jù)沿不同方向的擴(kuò)散強(qiáng)度，在各向異性擴(kuò)散濾波處理中，能夠更好地保護(hù)斷層等邊界信息。理論模型測(cè)試和實(shí)際地震資料應(yīng)用證明了該方法在提高資料信噪比的同時(shí)，可有效保護(hù)斷層等的邊界信息，并在地層連續(xù)區(qū)域增強(qiáng)同相軸的連續(xù)性，為斷層解釋提供良好的基礎(chǔ)數(shù)據(jù)。