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

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

1 東北石油大學(xué)地球科學(xué)學(xué)院， 大慶 163318 2 油氣資源與勘探技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室(長江大學(xué))， 武漢 430100

0 引言

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

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

1 方法原理與關(guān)鍵步驟

1.1 最大似然屬性原理

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

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

1.1.1 估計(jì)反射斜率

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

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

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

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

(1)

(2)

其梯度結(jié)構(gòu)張量為：

(3)

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

(4)

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

經(jīng)過平滑濾波后地震數(shù)據(jù)的梯度結(jié)構(gòu)張量GST根據(jù)矩陣特征分解，可表示為：

GST(x,y,t)=

(5)

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

(6)

(7)

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

1.1.2 相似性屬性的計(jì)算

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

(8)

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

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

(9)

式(9)即為進(jìn)行構(gòu)造導(dǎo)向平滑處理后得到的計(jì)算穩(wěn)定性增強(qiáng)的相似屬性.

1.1.3 最大似然屬性的計(jì)算

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

Maximun likelihood=max[1-semblance8].

(10)

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

1.2 建立拉普拉斯金字塔增強(qiáng)原理

拉普拉斯金字塔增強(qiáng)是空域?yàn)V波能量增強(qiáng)中銳化濾波增強(qiáng)的一種，其主要是在圖像的拉普拉斯金字塔分解過程中加入映射函數(shù)，得到由高斯金字塔演化而來的中間過程的拉普拉斯金字塔，最后得到需要的細(xì)節(jié)進(jìn)行增強(qiáng).

1.2.1 拉普拉斯金字塔分解

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

圖1 最大似然屬性計(jì)算示意圖(a) 斷面位置及采樣點(diǎn)； (b) 掃描斷面計(jì)算相似性； (c) 最小相似性對應(yīng)位置.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為層數(shù)；Rl為高斯金字塔第l層的行數(shù)；Cl為高斯金字塔第l層的列數(shù)；ω(m,n)是5×5的二維高斯內(nèi)核，其表達(dá)式為：

(12)

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

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

(13)

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

(14)

通過對地震信號的拉普拉斯金字塔分解，可以把原地震信號分解到不同的空間頻帶上，根據(jù)需要對相應(yīng)的細(xì)節(jié)進(jìn)行增強(qiáng)，如圖2所示.通過拉普拉斯金字塔也可對增強(qiáng)后的地震信號進(jìn)行重建，其表達(dá)式為：

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

(15)

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

1.2.2 地震數(shù)據(jù)的重映射

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

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

(16)

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

對重映射處理后的地震剖面構(gòu)造其拉普拉斯金字塔，然后再對增強(qiáng)后的地震剖面進(jìn)行重建，即可得到斷面波能量增強(qiáng)的地震剖面，有利于后續(xù)的斷層解釋與有利圈閉的識別.

1.3 關(guān)鍵步驟

在實(shí)際應(yīng)用過程中，通過最大似然屬性進(jìn)行斷裂識別和拉普拉斯金字塔斷面增強(qiáng)包括以下步驟：

1.3.1 斷裂特征分析與成像加強(qiáng)

通過對工區(qū)內(nèi)地震數(shù)據(jù)的分析，明確目的層內(nèi)斷裂構(gòu)造在地震上的反射特征.了解工區(qū)的構(gòu)造發(fā)育史，和斷裂形成的時(shí)期，以及斷裂特征類型，以便于后續(xù)相關(guān)參數(shù)的合理選擇.

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

1.3.2 最大似然屬性斷裂識別

最大似然屬性原理在上文中已經(jīng)敘述，從地震數(shù)據(jù)體中計(jì)算最大似然屬性時(shí)，有兩個(gè)關(guān)鍵參數(shù)：矩形窗口的水平計(jì)算因子、垂直計(jì)算因子.水平計(jì)算因子和垂直計(jì)算因子越大，計(jì)算周期越長，計(jì)算出最大似然屬性連續(xù)性越好.

1.3.3 拉普拉斯金字塔增強(qiáng)斷面

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

2 模型試算

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

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

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

3 實(shí)際應(yīng)用對比

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

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

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

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

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

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

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

圖9 X地區(qū)地震資料平面分析(a) 常規(guī)相干體水平切片； (b) 斷面波照明度增強(qiáng)后相干體水平切片.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 討論及結(jié)論

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

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