王 靜 張軍華 馮德永 李紅梅

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

0 引言

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

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

1 方法原理

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

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

(1)

式中：U為輸入的地震數(shù)據(jù)振幅信息；t為擴(kuò)散時間；g為擴(kuò)散系數(shù)。g的取值有三種情況： ①當(dāng)g取常數(shù)時，式(1)為線性各向同性方程，對數(shù)據(jù)進(jìn)行濾波后，能夠提高信噪比，但不能保護(hù)斷層等地質(zhì)體的邊界及方向信息； ②當(dāng)g為非負(fù)單調(diào)遞減函數(shù)[3]時，式(1)為非線性各向同性方程，濾波后容易產(chǎn)生階梯效應(yīng)和針孔效應(yīng)，而在強(qiáng)噪聲區(qū)域無去噪效果； ③使用擴(kuò)散張量D作為擴(kuò)散系數(shù)[4]，可以根據(jù)各向異性信息對數(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的過程中，一般采用具有對稱半正定的結(jié)構(gòu)張量Sρ計算局部方向信息，其形式為

(3)

(4)

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

(5)

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

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

(6)

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

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

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

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

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

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

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

(7)

Mline∈[0,1]、Mplane∈[0,1]。Mline描述某點(diǎn)鄰域與線狀結(jié)構(gòu)的相似程度，當(dāng)Mline=1時，指示線狀結(jié)構(gòu)；Mplane描述某點(diǎn)鄰域與面狀結(jié)構(gòu)的相似程度，當(dāng)Mplane=1時，指示面狀結(jié)構(gòu)。由Mline、Mplane得到檢測斷層等不連續(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è)計擴(kuò)散張量特征值

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

(9)

式中γ1、γ2和γ3為D的三個特征值，對應(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è)計中引入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è)計D特征值的方法可在增強(qiáng)擴(kuò)散濾波作用的同時，較好地保護(hù)斷層等特殊地質(zhì)結(jié)構(gòu)的邊界信息。

1.4 方法實現(xiàn)

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

(11)

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

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

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

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

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

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

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

(6)如果k

2 理論模型測試

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

圖1 Qdome模型立體顯示

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

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

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

2.1.1σ的影響

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

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

2.1.2ρ的影響

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

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

2.1.3 Δt的影響

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

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

2.1.4K的影響

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

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

2.2 模型試算

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

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

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

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

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

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

4 結(jié)論

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