基于最大似然屬性和拉普拉斯金字塔的斷面波增強方法建立

張冰， 徐嘉亮* ， 王維紅， 石穎， 王鵬

1 東北石油大學地球科學學院， 大慶 163318 2 油氣資源與勘探技術教育部重點實驗室(長江大學)， 武漢 430100

0 引言

斷層及斷裂構造對于改善儲層物性、聚集成藏等具有重要作用，因此在有利勘探區開發過程中，斷層的精細刻畫與有效識別占有重要地位(徐嘉亮等，2018a；余博文，2021).在復雜斷裂地區的勘探開發，由于受多期構造運動影響，斷層較為發育，斷層類型較多，在解釋過程中對斷層的識別和刻畫的精細程度要求較高，就導致斷層精細解釋和斷層平面組合變得更加困難(馬玉歌等，2020；李海晨等，2020).隨著地球物理及計算機相關技術的發展，斷層識別技術取得較大進步，出現不同的斷層識別方法.現在經過國內外學者的大量研究，已經形成較多斷層識別方法，如沿層相干屬性技術、相干體技術、方差體技術、曲率體屬性、優勢頻帶相位分析技術、螞蟻追蹤算法、邊緣檢測技術、面包切片技術等(Marfurt et al.，1998；Gersztenkorn and Marfurt，1999).相干體技術是現在使用最廣泛的斷層識別技術，第一代相干體技術(C1)主要利用互相關；第二代相干體技術(C2)是基于地震道相似性；第三代相干體技術(C3)是提取協方差矩陣的特征向量(Bahorich and Farmer，1995)；在此基礎上又衍生出一些其他相干體技術.螞蟻追蹤算法是根據螞蟻覓食行為而提出的斷層自動識別技術(陳桂和劉洋，2021).邊緣檢測技術來源于圖像處理，從最初應用到地震勘探領域到現在已經發展了較多的邊緣檢測技術，并在實際應用中取得較好的成果.Hale(2013)在研究斷層掃描、提取斷面和估算斷距時提出最大似然屬性，并在識別斷層時效果較好.此外，在成像過程中斷面波會在斷層處成像，因此增強斷面波能量可以在一定程度上對斷層進行精細識別與刻畫(徐嘉亮等，2015).

目前對于復雜構造斷層及小尺度斷層，利用識別度較高的斷面波進行斷層及斷面解釋是比較常用的方法(徐嘉亮等，2018b).如何提高斷面波識別度，降低繞射波等干擾成為提高斷面波識別度的關鍵(劉保金等，2012).本文利用斷層識別度更高的大似然屬性方法對斷層進行了有效識別，并利用拉普拉斯金字塔增強方法對識別的斷面波進行了增強，進而提高了斷面波的識別度.通過在X地區進行本文方法應用，結果表明在地震剖面上斷面波能量有所提高，得到更加清楚、連續性更好的斷裂，斷層組合關系更加清晰，取得較好的應用效果.

1 方法原理與關鍵步驟

1.1 最大似然屬性原理

最大似然法是用統計方法建立一組區分函數，在多類別分析中計算出各個類別樣本的歸屬概率，樣本屬于哪類由歸屬概率的大小決定(楊午陽，2017).

基于最大似然屬性識別斷層的原理源于對地震圖像的相似性分析.Hale提出面向斷面波的相似性理論，即應用矩形求取視窗中傾角和傾向的相似度.相似度屬性由相似性通過構造導向平滑得到，相似性計算的基礎是對地質體上沿反射界面不同方向變化率的估計，因此應先估計沿反射界面不同方向的反射斜率，然后進行最大似然屬性計算.

1.1.1 估計反射斜率

計算反射斜率方法有很多，其中梯度結構張量(GST)是一種有效的方法在精確表征和分析地震數據中反射結構的特征上.其計算原理如下：

假設原始地震數據為S(x,y,t)，其復數域的解析信號可表示為

S(x,y,t)=s(x,y,t)+jsH(x,y,t)

=A(x,y,t)ejPhs(x,y,t)，

(1)

(2)

其梯度結構張量為：

(3)

其中w(x,y,t)為權重系數，表示對梯度向量的各個分量進行光滑平均，以達到對地震數據平滑濾波處理的目的.其表達式為：

(4)

式中：σ為噪聲尺度參數，通常取值在0.1～3.0之間，x、y、t分別為地震數據中沿線方向、道方向和時間方向的變量.

經過平滑濾波后地震數據的梯度結構張量GST根據矩陣特征分解，可表示為：

GST(x,y,t)=

(5)

式中：λ1、λ2、λ3為梯度結構張量GST的特征值，且λ1>λ2>λ3>0，α1、α2、α3為梯度結構張量GST的特征值λ1、λ2、λ3對應的特征向量.即x、y方向的視傾角θx、θy可表示為：

(6)

(7)

式(6)和(7)的α1x(x,y,t)、α1y(x,y,t)、α1t(x,y,t)分別為式(5)中得到的最大特征值λ1對應的特征向量α1(x,y,t)分別在x、y和t方向上的三個分量.

1.1.2 相似性屬性的計算

根據估計的反射斜率，由式(8)計算相似性：

(8)

式中：u為矩形方向窗口中的地震數據；K為矩形方向窗口中地震數據的道數；下標k表示時窗中第k道信號；xk、yk分別為第k道信號在x、y方向上的距離；θx、θy分別為x、y方向上的視傾角；θxxk+θyyk為第k道信號在時間軸的時移量.

然而，應用此方法計算的相似性，在分子和分母較小的地方會有很大變化.針對這種不穩定性需要進行相應的平滑處理，平滑處理后即得到相似屬性.假設時窗大小為(2ω+1)ms，Z=ω/Δt，Δt表示采樣率，此時式(8)可變為：

(9)

式(9)即為進行構造導向平滑處理后得到的計算穩定性增強的相似屬性.

1.1.3 最大似然屬性的計算

通過矩形方向窗口對斷層傾角和傾向進行掃描得到的相似系數計算似然屬性，如圖1所示(以2D地震數據為例).假設斷層傾角θ掃描范圍為[θmaxmin]，掃描間隔為Δθ，斷層傾向φ掃描范圍為[φmaxmin]，掃描間隔為Δφ.在求取最大似然屬性過程中，對于某一個采樣點，分別對每一個掃描的傾角和傾向計算相似性，直到沿著相似度為最小值，即得到最大似然屬性，公式為：

Maximun likelihood=max[1-semblance8].

(10)

由公式(10)可知，最大似然屬性的數值范圍在0～1之間.這樣通過最大似然屬性對斷層及斷裂的有效識別分析，就可以為解釋斷層提供更好的支持.

1.2 建立拉普拉斯金字塔增強原理

拉普拉斯金字塔增強是空域濾波能量增強中銳化濾波增強的一種，其主要是在圖像的拉普拉斯金字塔分解過程中加入映射函數，得到由高斯金字塔演化而來的中間過程的拉普拉斯金字塔，最后得到需要的細節進行增強.

1.2.1 拉普拉斯金字塔分解

地震疊后數據的拉普拉斯金字塔構建是從高斯金字塔演變而來，因此首先應對地震信號進行高斯金字塔分解.假設原地震信號為F，并以原信號F作為高斯金字塔最底層地震信號G0；再對其進行低通濾波平滑和向下采樣，常采用水平方向和垂向均為1/2，即得到高斯金字塔第一層；不斷重復上述過程，就獲得縮小尺寸的地震信號；層級越高，地震信號越小，從而構成高斯金字塔，其表達式為：

圖1 最大似然屬性計算示意圖(a) 斷面位置及采樣點； (b) 掃描斷面計算相似性； (c) 最小相似性對應位置.Fig.1 Schematic diagram of maximum likelihood attribute calculation(a) Section position and sampling point; (b) Similarity of scanning section calculation; (c) Minimum similarity corresponding position.

(11)

式中：N為層數；Rl為高斯金字塔第l層的行數；Cl為高斯金字塔第l層的列數；ω(m,n)是5×5的二維高斯內核，其表達式為：

(12)

ω(m,n)實際是高斯低通濾波器,經過n次重復之后得到Gn，此時Gn只包含級少數的像素點.

1≤l≤N,0≤i≤Rl,0≤j≤Cl，

(13)

式中：N為層數；ω(m,n)是5×5的二維高斯內核.拉普拉斯金字塔第l層地震信號Ll表達式為：

(14)

通過對地震信號的拉普拉斯金字塔分解，可以把原地震信號分解到不同的空間頻帶上，根據需要對相應的細節進行增強，如圖2所示.通過拉普拉斯金字塔也可對增強后的地震信號進行重建，其表達式為：

圖2 構造過程(a) 高斯金字塔； (b) 對高斯金字塔進行變換； (c) 拉普拉斯金字塔.Fig.2 Construction process(a) Gaussian pyramid; (b) Transform of Gaussian pyramid; (c) Laplace pyramid.

(15)

當l=0時，得到重建后地震信號，假如在變換過程中沒有做相關參數修改，得到的重建信號與原始地震信號完全相同.

1.2.2 地震數據的重映射

在對地震剖面進行初始變化時，首先假設一個能量參數ε，如果地震剖面斷面部分的能量變化小于ε，就認為是細節部分，如果大于ε，就認為是邊緣部分.在重映射算法中，當計算輸出金字塔系數時，某點強度(強度是指單通道像素值的大小)可表示為g0=Gl0(x0,y0).在構建拉普拉斯金字塔過程中地震斷面波能量在ε附近的被認為是細節處理，其表達式為：

φd(i)=g0+sign(i-g0)εfd(|i-g0|/ε)，

(16)

式中：g0為當前要處理的像素點；sign(i-g0)為亮度差異(i-g0)的符號值；ε可通過限定細節部分像素點的個數≤nd來確定當前層中ε的值；fd(x)=ex是映射函數，用于控制增強細節部分；|i-g0|為像素點i與數據周圍像素點的差異值；下標d表示細節detail的首字母.

對重映射處理后的地震剖面構造其拉普拉斯金字塔，然后再對增強后的地震剖面進行重建，即可得到斷面波能量增強的地震剖面，有利于后續的斷層解釋與有利圈閉的識別.

1.3 關鍵步驟

在實際應用過程中，通過最大似然屬性進行斷裂識別和拉普拉斯金字塔斷面增強包括以下步驟：

1.3.1 斷裂特征分析與成像加強

通過對工區內地震數據的分析，明確目的層內斷裂構造在地震上的反射特征.了解工區的構造發育史，和斷裂形成的時期，以及斷裂特征類型，以便于后續相關參數的合理選擇.

最大似然屬性斷裂識別和拉普拉斯金字塔增強斷面的理論基礎仍然是在地震剖面上反射波同相軸的不連續性，噪聲、較低的分辨率和地層巖性的突然變化也可能導致反射波同相軸不連續，因此在地震數據處理過程中要注重疊前相關噪聲的有效壓制和提高地震數據的分辨率，在Kirchhoff 疊前時間偏移過程中，要注重偏移孔徑和反假頻因子的合理選擇，在偏移成像后通過隨機噪聲壓制和疊后提高分辨率處理技術，來進一步提高斷面波成像效果.

1.3.2 最大似然屬性斷裂識別

最大似然屬性原理在上文中已經敘述，從地震數據體中計算最大似然屬性時，有兩個關鍵參數：矩形窗口的水平計算因子、垂直計算因子.水平計算因子和垂直計算因子越大，計算周期越長，計算出最大似然屬性連續性越好.

1.3.3 拉普拉斯金字塔增強斷面

根據最大似然屬性對斷層的識別，通過拉普拉斯金字塔進行斷面增強，首先對地震信號進行高斯金字塔搭建，再通過從每一層減去高斯金字塔上一層信號且進行內核卷積的地震信號，就得到拉普拉斯金字塔，在搭建拉普拉斯金字塔過程中加入一個映射函數來增強斷面部分，并用拉普拉斯金字塔對地震信號進行重建，即得到斷面增強后的地震數據.

2 模型試算

本文以圖3所示的帶有兩條斷層的四層速度模型進行方法測試.該模型在水平和垂直方向上分別有500個網格點，網格間距為2 m，炮間距為20 m，每炮共有400個檢波器接收記錄.第一炮與最后一炮的位置分別在200 m和800 m的位置.該觀測系統能夠保證左側斷層滿覆蓋接收，右側斷裂不滿覆蓋接收，進而驗證在滿覆蓋與不滿覆蓋情況下本文方法的實用性.

圖3 帶有兩條斷層的正演速度模型Fig.3 Forward velocity model with two faults

由于正演模型斷層傾角較大，斷面波能量較低.圖4a為Kirchhoff疊前時間偏移原始剖面，圖4b為利用本文方法增強斷層后的疊前時間偏移剖面，兩條斷層能量均有所增強，連續性變好且識別度變高，更有助于利用斷面波進行斷層解釋.

3 實際應用對比

本文將最大似然屬性和拉普拉斯金字塔增強在X地區進行應用，并取得較好效果.X地區主控及分支斷層在剖面上組合成花狀構造、階梯狀構造和反Y字型構造，同時還伴有一些小的走滑斷裂.在原疊加剖面上這些構造位置關系不明確，造成斷層解釋和平面組合困難.此外，由于河道交錯疊置復雜，河道邊界的異常響應也干擾斷層識別.為落實斷層解釋，采用最大似然屬性識別斷層和拉普拉斯金字塔增強斷面，并取得一定效果.圖5a為增強前測線A的偏移剖面，斷面處連續性較差且斷層搭接關系不明確，且存在假斷層響應影響斷層解釋精度，圖5b是增強后測線A的疊加剖面，斷面波能量有所增強，斷層連續性較好且搭接關系明確.圖6a和圖6b為圖5斷層增強前后剖面的局部方法，由局部放大比較中能夠明顯看出利用本文提出的方法能夠增強斷面波的能量，進而提高斷層的識別度.

圖5 增強前后測線A的疊加剖面對比(a) 原始疊加剖面； (b) 增強后疊加剖面.Fig.5 Comparison of stacked profiles A before and after enhancement(a) Original section； (b) Enhanced section.

圖6 增強前后測線A的疊加剖面對比(a) 原始疊加剖面； (b) 增強后疊加剖面.Fig.6 Comparison of stacked profile A before and after enhancement(a) Original section； (b) Enhanced section.

圖7和圖8是對測線B疊前時間偏移剖面的增強及其布局方法的比較，可看到增強前后斷面波能量明顯提高，斷層連續性較好且搭接關系更明確.

圖7 增強前后測線B的疊加剖面對比(a) 原始疊加剖面； (b) 增強后疊加剖面.Fig.7 Comparison of stacked profile B before and after enhancement(a) Original section； (b) Enhanced section.

圖8 增強前后測線B的疊加剖面對比(a) 原始疊加剖面； (b) 增強后疊加剖面.Fig.8 Comparison of stacked profile B before and after enhancement(a) Original section； (b) Enhanced section.

圖9a為該區T2(館陶組)層段常規相干體水平切片，可以發現X地區主斷裂為近東西向展布，但該相干體橫向分辨能力較低，斷層邊界模糊影響其真實性.圖9b表示斷面波照明度增強后最大似然水平切片，最大似然體橫向分辨能力較高，次級斷裂與主斷裂的接觸關系較清楚，斷層細節更加清楚、連續，且刻畫斷層邊界更加清晰，斷層組合關系明確，有利于斷層識別，更加符合斷裂的實際地質特征.

圖9 X地區地震資料平面分析(a) 常規相干體水平切片； (b) 斷面波照明度增強后相干體水平切片.Fig.9 X area seismic data plane analysis(a) Horizontal slice of conventional coherent body； (b) Horizontal slice of coherent body after enhancing fault-surface wave illumination.

4 討論及結論

通過分析最大似然屬性識別斷層和基于拉普拉斯金字塔算法增強斷面識別度的原理，以及在X地區的應用，與直接應用常規地震資料對斷層精細解釋和識別相比，應用此方法，可以顯著提高斷面波在疊加剖面上的振幅能量，有利于后續分析斷層位置以及斷層展布特征、斷層平面組合關系，較好地實現斷層精細解釋和斷裂系統的有效識別.然而在最大似然屬性識別斷裂時，要注意矩形窗口的水平計算因子、垂直計算因子的合理選擇，保證運算時間的合理性，且計算的最大似然屬性具有較好的連續性，同時，在研究中發現最大似然屬性在某些高振幅的地質體邊界處仍具有一定異常響應，影響斷層識別，需要進一步改善提高.此外在基于拉普拉斯金字塔算法增強斷面時，由于地震圖像中存在與斷面部分強度變化相似的異常響應，在拉普拉斯金字塔算法增強斷面時異常相應也會得到增強，對于該算法的不足之處，還需要針對原地震剖面中斷層的特點進行改進，以獲得更好的斷面增強效果.

本文分析了最大似然屬性識別斷層和基于拉普拉斯金字塔算法增強斷面的原理，并在X地區進行應用，即通過最大似然屬性對工區數據體進行掃描，并計算數據體采樣點之間的相似性，以獲得工區內最可能發育斷層的位置及概率，然后應用拉普拉斯金字塔增強算法進行斷面增強，即在構造拉普拉斯金字塔的過程中加入一個映射函數，使斷面得到增強.在X地區的實際應用表明，最大似然屬性在疊加剖面上識別斷層位置比較準確，再通過拉普拉斯金字塔算法增強斷面，可得到斷面波能量較高的疊加剖面，在剖面上斷層平面組合關系更明確，提高了斷層解釋精度，為后續有利圈閉識別和井位優選提供了依據.