利用不連續性的各向異性擴散濾波方法識別斷層

王 靜 張軍華 馮德永 李紅梅

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

0 引言

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

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

本文在前人的研究基礎上，分析了結構張量特征值與三維地震圖像的局部結構特征之間的關系，利用斷層置信度參數重新構建擴散張量特征值，進而重構擴散張量，實現基于不連續性的三維各向異性擴散濾波方法，可有效提高地震資料信噪比，保護斷層等不連續的邊界信息。

1 方法原理

1.1 各向異性擴散濾波方程

各向異性擴散濾波的擴散過程等價于熱傳導過程，表示為[20]

(1)

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

(2)

式中U0(x,y,z)為原始三維地震數據，(x，y，z)為某點的空間坐標。

1.2 結構張量的構建

在構建D的過程中，一般采用具有對稱半正定的結構張量Sρ計算局部方向信息，其形式為

(3)

(4)

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

(5)

式中α∈[0,1]為穩定系數，用于控制高斯濾波后地震數據對空間的二階導數的貢獻率。

對Sρ進行特征分解，得

(6)

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

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

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

(2)λ1>λ2≈λ3≈0，在ζ1方向梯度變化較大，另外兩個特征方向的梯度變化很小且幾乎相等，指示平面線性地質體，對應模型為近似面狀線性結構；

(3)λ1>λ2>λ3≈0，其中一個方向梯度變化近似為零，指示在斷層等不連續的邊界處出現不連續的間斷信息，對應模型為線狀結構模型；

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

van Kempen等[22]根據結構張量特征值與地震圖像結構的關系，用線狀置信度量Mline和面狀置信度量Mplane表征圖像結構[18]，即

(7)

Mline∈[0,1]、Mplane∈[0,1]。Mline描述某點鄰域與線狀結構的相似程度，當Mline=1時，指示線狀結構；Mplane描述某點鄰域與面狀結構的相似程度，當Mplane=1時，指示面狀結構。由Mline、Mplane得到檢測斷層等不連續結構邊界信息的斷層置信度量Mfault[21,23]

(8)

表1為置信度參數與地層結構以及特征值之間的關系。

表1 置信度參數與地層結構以及特征值之間的關系

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

1.3 利用Mfault設計擴散張量特征值

構建擴散張量D時，為了保證擴散方向沿著構造方向，D的特征向量應與Sρ的一致，因此

(9)

式中γ1、γ2和γ3為D的三個特征值，對應的特征向量也為ζ1、ζ2和ζ3。

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

為了克服常規擴散濾波方法不能保護數據邊緣信息的缺點，在D的特征值設計中引入Mline、Mplane以及Mfault，即

(10)

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

綜上所述，利用Mfault設計D特征值的方法可在增強擴散濾波作用的同時，較好地保護斷層等特殊地質結構的邊界信息。

1.4 方法實現

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

(11)

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

本文方法的實現流程為：

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

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

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

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

(5)構建D后，根據式(11)計算第k+1次迭代結果；

(6)如果k

2 理論模型測試

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

圖1 Qdome模型立體顯示

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

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

2.1 濾波參數的選擇

2.1.1σ的影響

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

圖3 不同噪聲尺度的濾波結果

2.1.2ρ的影響

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

圖4 不同整合尺度的濾波結果

2.1.3 Δt的影響

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

圖5 不同迭代步長的濾波結果

2.1.4K的影響

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

圖6 不同最大迭代次數的濾波結果

2.2 模型試算

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

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

3 實際地震資料應用

選取斷裂眾多、結構復雜、勘探難度較大的勝利油田X斷塊地震資料作為研究數據，應用文中方法進行濾波處理并計算相干屬性。圖9為本文方法濾波前、后的Inline3570剖面。由圖可見，濾波后地震剖面的同相軸更連續、斷點更清晰(圖9b黃色和藍色箭頭處)，斷層更易識別，說明本文方法在去噪的同時能有效保護斷層等不連續性邊界的信息。圖10為濾波前、后的相干時間切片。由圖可見，利用濾波后的相干時間切片(圖10b)可準確地描述斷層形態及展布，如在濾波前的相干時間切片中一些不易識別的小斷層(圖10a中黃色和紅色箭頭處)，在濾波后的相干時間切片上得到準確、清晰地識別，這也證明了該濾波方法在提高資料信噪比和保持斷層等邊緣信息方面的有效性和和準確性。

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

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

4 結論

針對三維地震資料保邊濾波處理的需要，本文研究了三維各向異性擴散濾波方法，并采用結構張量分析了三維地質結構的局部特征。基于斷層置信度參數重新設計了擴散張量的特征值，可以控制地震數據沿不同方向的擴散強度，在各向異性擴散濾波處理中，能夠更好地保護斷層等邊界信息。理論模型測試和實際地震資料應用證明了該方法在提高資料信噪比的同時，可有效保護斷層等的邊界信息，并在地層連續區域增強同相軸的連續性，為斷層解釋提供良好的基礎數據。