動(dòng)態(tài)時(shí)間規(guī)整的斷層增強(qiáng)方法

董 林 宋維琪 胡建林 曾 超 趙寶銀 高文中

(①中國(guó)石油大學(xué)(華東)，山東青島 266580；②中國(guó)石油冀東油田勘探開發(fā)研究院物探所，河北唐山063000)

0 引言

在陸相斷陷盆地中，斷層往往控制圈閉的形成，因此準(zhǔn)確解釋斷層具有重要意義。地震剖面中的同相軸錯(cuò)斷、扭曲是明顯的斷層響應(yīng)，但是當(dāng)斷層斷距較小或地震資料信噪比較低時(shí)，斷層解釋的合理性較低，影響圈閉可靠性評(píng)價(jià)。因此，為提高斷層解釋精度，人們不斷研究有效識(shí)別斷層的方法[1-3]。

早期人們主要利用各種地震屬性體(如相干[4-6]、方差[7-8]、曲率[9-10]等屬性)識(shí)別斷層，雖然可在一定程度上反映斷層的空間特征，但對(duì)噪聲很敏感，對(duì)復(fù)雜斷裂的檢測(cè)效果不佳。因此有必要對(duì)屬性體進(jìn)行斷層增強(qiáng)處理，以更精確地識(shí)別斷層。

隨著人們不斷研究斷層識(shí)別方法，出現(xiàn)了許多針對(duì)斷層增強(qiáng)的處理方法。Barnes[11]利用組合成分濾波器有效去除屬性信息中的非斷層成分。Donias等[12]沿著斷層方向增強(qiáng)斷層屬性信息，根據(jù)斷層的平面特征有效增強(qiáng)斷層特征。Lavialle等[13]利用基于偏微分方程的非線性濾波技術(shù)，在自動(dòng)提取斷層之前去除非斷層信息。Machado[14]設(shè)計(jì)了高斯拉普拉斯算子，沿?cái)鄬臃较蛟鰪?qiáng)斷層屬性的同時(shí)抑制非斷層信息。Cohen等[15]、Wu等[16]通過掃描所有斷層走向和傾角組合平滑斷層屬性，有效增強(qiáng)斷層特征。上述斷層增強(qiáng)算法大多首先根據(jù)屬性體計(jì)算斷層方向，然后增強(qiáng)斷層方向的斷層信息，壓制非斷層方向的其他信息(如殘留地層響應(yīng)、噪聲等)，可得到理想的計(jì)算結(jié)果，但計(jì)算精度依賴于斷層方向的計(jì)算精度。

動(dòng)態(tài)時(shí)間規(guī)整(Dynamic Time Warping，DTW)算法最早由Sakoe等[17]提出并用于語(yǔ)音識(shí)別，是一種基于動(dòng)態(tài)規(guī)劃的模板匹配算法，根據(jù)兩個(gè)信號(hào)的相似性局部拉伸或壓縮信號(hào)，使兩個(gè)信號(hào)最佳匹配。在地震勘探領(lǐng)域，人們將DTW算法用于圖像對(duì)齊、斷層斷距計(jì)算等方面。Anderson等[18]將DTW算法用于地球物理問題，并稱這種算法為“動(dòng)態(tài)波形匹配”。Hale[19]利用DTW算法匹配縱波與轉(zhuǎn)換波地震數(shù)據(jù)、估計(jì)斷層斷距，取得了較好的效果。Wu等[20]利用DTW算法自動(dòng)解釋斷層，在低信噪比地震數(shù)據(jù)中準(zhǔn)確識(shí)別了斷層。吳天麒等[21]提出了基于互相關(guān)的動(dòng)態(tài)時(shí)間規(guī)整算法(CDTW)，可將角道集同相軸拉平且對(duì)波形畸變具有良好的魯棒性。

由于DTW方法能穩(wěn)定、快速地從兩個(gè)序列的誤差矩陣中求出誤差最小路徑，因此本文利用DTW方法計(jì)算斷層圖像的最優(yōu)斷層線，并對(duì)最優(yōu)斷層線平滑處理，得到的斷層圖像較原始屬性圖像的斷層特征更清晰、連續(xù)。模型和實(shí)際資料處理結(jié)果表明，所提方法能抑制斷層圖像的非斷層信息，斷層增強(qiáng)結(jié)果更清晰、連續(xù)，具有較好的魯棒性和適用性。

1 方法原理

1.1 DTW方法原理

對(duì)于給定的兩個(gè)離散序列(實(shí)際上不一定與時(shí)間有關(guān))，利用DTW方法能計(jì)算它們的相似程度，并且給出一個(gè)可最大限度降低兩個(gè)序列距離的點(diǎn)到點(diǎn)的匹配方案。設(shè)f(t)為參考序列，g(t)為待匹配序列，兩個(gè)序列之間的對(duì)齊誤差e(i，n)為序列f(i)與g(i+n)之間的歐氏距離(n為序列g(shù)(t)相對(duì)序列f(t)的時(shí)間滯后)，f(t)與g(t)匹配后的累積對(duì)齊誤差為

e(i，n)≡[f(i)-g(i+n)]2

(1)

DTW求解的目標(biāo)是找到一個(gè)時(shí)間序列

T(t)=[T(0)，T(1)，…，T(N-1)]

(2)

使f(t)與g(t)最相似，即

f(i)=g[i+T(i)]

(3)

且滿足

(4)

其中l(wèi)(t)為f(t)與g(t)對(duì)齊時(shí)(t=0，1，2，…，N-1)對(duì)應(yīng)對(duì)齊誤差矩陣中從點(diǎn)(0，0)到點(diǎn)(N-1，N-1)的任意一條路徑(N為時(shí)間序列的長(zhǎng)度)，且D[l(t)]滿足

(5)

并服從約束條件

|T(i)-T(i-1)|≤1

(6)

當(dāng)f(t)與g(t)對(duì)齊時(shí)，式(6)確保i+T(i)隨著采樣點(diǎn)的增加不會(huì)過快增加或減少。

1.2 最優(yōu)斷層線計(jì)算方法

斷層屬性體局部極大值(或極小值)處往往指示斷層，斷層線的脊線由斷層圖像中連續(xù)屬性極值點(diǎn)組成。若將斷層屬性圖像看作DTW的對(duì)齊誤差矩陣，將DTW求解對(duì)齊誤差最小值路徑相應(yīng)地變?yōu)榍蠼庾畲笾德窂?，則求取屬性圖像中的最優(yōu)斷層線變?yōu)榍笕D像中最大值路徑問題

(7)

式中：j(i)為路徑上的縱坐標(biāo)i(對(duì)應(yīng)垂直方向)對(duì)應(yīng)的橫坐標(biāo)(對(duì)應(yīng)水平方向)；C[i，j(i)]為路徑上對(duì)應(yīng)的屬性值，其約束條件為

|j(i+1)-j(i)|≤ε0<ε≤1

(8)

式中ε為計(jì)算斷層線時(shí)的斜率約束，以保證在計(jì)算最優(yōu)路徑過程中相鄰點(diǎn)之間不會(huì)偏離橫坐標(biāo)方向，從而產(chǎn)生過大的橫向位移。因此，在計(jì)算最優(yōu)斷層線時(shí)，為了確保斜率約束的作用，斷層必須是近于垂直的。文中以只有一條近垂直斷層的斷層屬性圖像(圖1b)為例，說明利用DTW計(jì)算最優(yōu)斷層線的方法，具體步驟如下。

(1)選取斷層上的一個(gè)點(diǎn)作為控制點(diǎn)，將斷層圖像中控制點(diǎn)左、右兩側(cè)錐形區(qū)域的數(shù)據(jù)置0(圖1c)，確保計(jì)算出的最優(yōu)斷層線從控制點(diǎn)通過。

(2)本文采用Hale[19]提出的三階段算法(非線性平滑、正向積累和反向追蹤)求解斷層圖像中的最大值路徑(最優(yōu)斷層線)，求解過程如下。

1)非線性平滑

f(0，j)=g(0，j)

(9)

f(i，j)=g(i，j)+

(10)

式中：d=?1/ε?為1/ε的取整；i=1，2，3，…，N-1；g(i，j)為輸入屬性圖像；f(i，j)為正向積累(從上至下)計(jì)算的圖像值，橫坐標(biāo)i對(duì)應(yīng)圖像的行方向，縱坐標(biāo)j對(duì)應(yīng)圖像的列方向。

相似地，反向積累(從下至上)計(jì)算的圖像值為

圖1 最優(yōu)斷層線計(jì)算結(jié)果(a)原始地震圖像； (b)斷層屬性圖像； (c)控制點(diǎn)設(shè)置； (d)非線性平滑； (e)計(jì)算出的最優(yōu)斷層線； (f)斷層增強(qiáng)結(jié)果

b(N-1，j)=g(N-1，j)

(11)

b(i，j)=g(i，j)+

(12)

式中i=N-2，N-3，…，0。圖像的雙邊平滑為

s(i，j)=f(i，j)+b(i，j)-g(i，j)

(13)

在正向、反向計(jì)算過程中加入低斜率約束平滑，使平滑后的斷層圖像更平滑、連續(xù)，確保計(jì)算出的斷層線更準(zhǔn)確。為了確保計(jì)算出的斷層線經(jīng)過控制點(diǎn)，同樣在平滑后圖像中將控制點(diǎn)左、右錐形區(qū)域的數(shù)據(jù)置0(圖1d)。

2)正向積累

采用正向積累(式(10))計(jì)算平滑后圖像的正向積累圖像a。

3)反向追蹤

首先找出a中第N-1行中的最大值的點(diǎn)，其位于第l列。然后對(duì)a從第N-1行中的最大值的點(diǎn)開始向前依次搜索前一行的最大值，直到找到最大值路徑上的所有點(diǎn)(圖1e中的黑線)，得到最大值路徑，即

j(N-1)=l

(14)

j(i-1)=

(15)

式中i=N-1，N-2，…，1。

(3)計(jì)算出最優(yōu)斷層線后，僅保留原始斷層圖像中最優(yōu)斷層線上的屬性值，再沿?cái)鄬泳€對(duì)斷層屬性值平滑處理，即可得到最終斷層增強(qiáng)結(jié)果(圖1f)。若經(jīng)過控制點(diǎn)的最大值路徑記為L(zhǎng)(i)=[i，j(i)]，i=0，1，2，…，N-1，則該斷層線上的新屬性值為

Cnew[i，j(i)]=〈C[L(i)]〉

(16)

式中：C[L(i)]為斷層線的原始屬性值；〈·〉表示沿路徑的一維高斯平滑。

1.3 斷層增強(qiáng)處理流程

實(shí)際斷層屬性圖像中的斷層不能直接用最優(yōu)斷層線計(jì)算方法增強(qiáng)處理。本文的斷層增強(qiáng)處理方法是將復(fù)雜斷層分解為局部窗口內(nèi)的簡(jiǎn)單斷層，然后在局部窗口內(nèi)計(jì)算最優(yōu)斷層線，最后綜合局部最優(yōu)斷層線得到最終的復(fù)雜斷層增強(qiáng)結(jié)果。以原始斷層圖像(圖2a)中的斷層為例，說明斷層增強(qiáng)的具體實(shí)現(xiàn)步驟。

1.3.1 斷層圖像濾波處理

為了準(zhǔn)確選取分布在斷層上的種子點(diǎn)，首先采用組合成分濾波器[9]對(duì)斷層屬性圖像濾波處理，去除斷層圖像中的非斷層信息。

圖2 斷層圖像濾波及種子點(diǎn)選取效果(a)原始斷層圖像； (b)去除非斷層信息斷層圖像； (c)去除非斷層信息前選擇的種子點(diǎn)； (d)去除非斷層信息后選擇的種子點(diǎn)

組合成分濾波器是一種組合濾波器，相當(dāng)于構(gòu)建了一個(gè)協(xié)方差矩陣

(17)

其中

(18)

(19)

對(duì)式(17)進(jìn)行特征分解得到相應(yīng)的特征值λ1、λ2、λ3，其中λ1≥λ2≥λ3，對(duì)應(yīng)的特征向量分別為e1、e2、e3。向量e1和e2能確定當(dāng)前窗口數(shù)據(jù)的最佳擬合平面(斷層平面)，e3垂直于最佳擬合平面。根據(jù)特征向量分別計(jì)算當(dāng)前點(diǎn)面狀構(gòu)造置信度Λ=(λ2-λ3)/λ1、當(dāng)前點(diǎn)所在斷層主方向及當(dāng)前點(diǎn)到擬合斷層面的歐式距離R。再根據(jù)這3個(gè)參數(shù)設(shè)計(jì)3個(gè)帶通濾波器剔除非高角度的擬合平面及Λ較小、R較大的點(diǎn)。由這3個(gè)濾波器組成1個(gè)組合濾波器可剔除斷層圖像中的非斷層信息(圖2b)。

1.3.2 種子點(diǎn)選取

為了獲取均勻分布在斷層上的種子點(diǎn)(圖2c、圖2d)，在去除非斷層信息的斷層圖像中，使用非極大值抑制篩選fast(Features from Accelerated Segment Test)特征點(diǎn)的方法獲取種子點(diǎn)。該方法先采用非極大值抑制的方法保留去除非斷層信息后斷層屬性圖像脊線上的值，將其他區(qū)域的數(shù)據(jù)置零，得到細(xì)化的屬性圖像。然后從細(xì)化的屬性中選擇種子候選點(diǎn)(保留斷層圖像屬性值中大于某個(gè)閾值的所有像素點(diǎn)作為候選點(diǎn))。再按屬性值從大到小的順序依次選擇所有滿足條件的候選點(diǎn)作為種子點(diǎn)。在選擇種子點(diǎn)時(shí)，計(jì)算當(dāng)前候選點(diǎn)和所有候選點(diǎn)之間的歐氏距離，只有候選的種子點(diǎn)與已選擇的種子點(diǎn)的最小距離大于預(yù)定義的距離閾值r時(shí)，當(dāng)前候選點(diǎn)才被選為種子點(diǎn)。

1.3.3 種子點(diǎn)斷層方向計(jì)算

由于斷層不一定是近于垂直的，因此計(jì)算前需要估算每個(gè)種子點(diǎn)的斷層方向(圖3a)，這樣可以沿著斷層方向選取數(shù)據(jù)進(jìn)行計(jì)算，確保在計(jì)算最優(yōu)路徑時(shí)局部窗口內(nèi)斷層是近于垂直的(圖3b)。結(jié)合先驗(yàn)信息與計(jì)算出的斷層方向，進(jìn)一步去掉方向與斷層方向不一致的種子點(diǎn)。

為了更好地估計(jì)多方向相交的斷層方向，采用方向掃描的方法大致估算每個(gè)種子點(diǎn)斷層的走向(傾角)。這里定義每個(gè)掃描方向的置信度為

(20)

式中：C(θ，Pi3)為斷層屬性值，θ為當(dāng)前掃描的角度，Pi3為當(dāng)前掃描方向的點(diǎn)；H為傾角掃描數(shù)量。從0°～180°每隔一定角度計(jì)算當(dāng)前角度的置信度，Q(θ)最大時(shí)對(duì)應(yīng)的θ即為當(dāng)前種子點(diǎn)所在斷層的斷層傾角(圖3a中的紅線)。由于本文研究的斷層為高角度斷層且不需要精確估計(jì)斷層方向，因此掃描間隔(10°)較大，且傾角掃描范圍為60°～120°。

圖3 斷層方向計(jì)算及沿?cái)鄬臃较蜻x取數(shù)據(jù)效果(a)斷層方向掃描； (b)沿?cái)鄬舆x取計(jì)算數(shù)據(jù) 沿紅框的局部窗口選擇數(shù)據(jù)，斷層近于垂直

1.3.4 種子點(diǎn)最優(yōu)路徑計(jì)算

對(duì)于每個(gè)種子點(diǎn)，本文選取一個(gè)以種子點(diǎn)為中心的矩形窗口，在這個(gè)窗口內(nèi)計(jì)算通過種子點(diǎn)的最優(yōu)路徑。假設(shè)斷層傾角掃描的方向?yàn)閡，垂直于斷層傾向?yàn)関，在選擇矩形窗口時(shí)，使矩形的長(zhǎng)邊平行于u。在這樣的局部窗口內(nèi)斷層是近于垂直的，此時(shí)在每個(gè)種子點(diǎn)局部窗口內(nèi)計(jì)算斷層線的方法與最優(yōu)斷層線計(jì)算方法一致。

1.3.5 斷層增強(qiáng)處理

由斷層增強(qiáng)結(jié)果(圖1f)可見，平滑后路徑上的屬性分布與原始屬性分布(圖1b)一致，即原始屬性值越大，平滑后路徑上的屬性值也越大。將所有種子點(diǎn)的計(jì)算結(jié)果累加，經(jīng)歸一化得到最終的斷層增強(qiáng)結(jié)果。在完成累加計(jì)算后，經(jīng)過多個(gè)路徑的點(diǎn)可能是真正位于斷層上的點(diǎn)，往往具有較高的屬性值計(jì)算結(jié)果，而經(jīng)過路徑少的點(diǎn)一般指示噪聲區(qū)域，往往具有較低的屬性值計(jì)算結(jié)果，因此可以設(shè)置一個(gè)閾值去除非斷層上的點(diǎn)。若點(diǎn)(i，j)有m個(gè)種子點(diǎn)的最優(yōu)路徑經(jīng)過，第m個(gè)種子點(diǎn)在點(diǎn)(i，j)平滑后的新屬性值為Cnew(m)，則對(duì)所有路徑進(jìn)行平滑處理并累加后，點(diǎn)(i，j)最終的新屬性值為

(21)

2 模型數(shù)據(jù)試算

為了驗(yàn)證本文方法的可行性，采用卷積神經(jīng)網(wǎng)絡(luò)(CNN)訓(xùn)練地震斷層模型的方法[22]制作未加噪(圖4a)、加噪(圖5a)地震模型進(jìn)行驗(yàn)證，模型中存在平行和相交斷層。

圖4 本文方法對(duì)未加噪地震模型的斷層增強(qiáng)處理效果(a)地震正演模型； (b)方差屬性； (c)圖b增強(qiáng)結(jié)果； (d)相似屬性； (e)圖d增強(qiáng)結(jié)果 種子點(diǎn)的預(yù)定義歐式距離閾值r=6，傾角掃描間隔Δθ=10°，每個(gè)種子點(diǎn)的矩形窗尺度為51×21

圖5 本文方法對(duì)加噪地震模型斷層增強(qiáng)處理效果(a)地震正演模型； (b)相似屬性； (c)增強(qiáng)結(jié)果

由未加噪模型方差屬性(圖4b)的增強(qiáng)結(jié)果(圖4c)、相似屬性(圖4d)的增強(qiáng)結(jié)果(圖4e)可見，相交斷層的交切關(guān)系清楚，在保持?jǐn)鄬于厔?shì)的同時(shí)，斷層增強(qiáng)結(jié)果更連續(xù)、清晰。由加噪模型相似屬性(圖5b)的增強(qiáng)結(jié)果(圖5c)可見，完全去除了非斷層信息，相交斷層的交切關(guān)系清楚，斷層增強(qiáng)結(jié)果更連續(xù)、清晰。上述結(jié)果說明本文方法的抗噪性較好，具有較強(qiáng)的魯棒性。

此外，對(duì)比圖4e和圖5c可見，兩者存在細(xì)小的局部差異，局部斷層線相差1個(gè)像素點(diǎn)。這是由于噪聲影響相似屬性的局部計(jì)算結(jié)果所致，說明本文算法的增強(qiáng)結(jié)果依賴于屬性體的斷層檢測(cè)精度，即屬性體的斷層檢測(cè)精度越高，增強(qiáng)結(jié)果越接近于斷層的真實(shí)形態(tài)。

3 實(shí)際資料處理

為了進(jìn)一步測(cè)試本文算法對(duì)實(shí)際數(shù)據(jù)的斷層增強(qiáng)效果，選取實(shí)際地震資料進(jìn)行驗(yàn)證。地震數(shù)據(jù)取自中國(guó)東部F油田，該油田處于陸相斷陷湖盆，盆地內(nèi)斷裂極其發(fā)育，對(duì)油氣的形成具有控制作用。由于斷層間距較小，小斷層發(fā)育，在相似屬性剖面上斷層響應(yīng)互相干擾，噪聲導(dǎo)致斷層特征不清晰，存在明顯的地層殘留響應(yīng)，在縱向上某些斷層呈斷續(xù)的“線段”(圖6a)。經(jīng)過斷層增強(qiáng)處理，基本去除了非斷層信息，斷層更連續(xù)、清晰(圖6b)。

圖7為F油田的斷層相似屬性水平切片及其斷層增強(qiáng)結(jié)果。由圖可見：①在相似屬性水平切片中斷層特征不明顯，噪聲導(dǎo)致部分?jǐn)鄬犹卣鞑磺逦?圖7a)；②經(jīng)過斷層增強(qiáng)處理，采用傾角約束得到的最優(yōu)斷層線是沿一定角度緩慢變化的曲線(圖7b)，不僅保持了原有斷層特征，斷層分布、延伸特征更清晰、干脆，明顯增強(qiáng)了斷層的線性特征，相應(yīng)地去掉了一些非線性(圖7a中的“圓圈狀”信息)局部特征。

圖6 F油田的二維地震剖面A的斷層相似屬性(a)及其斷層增強(qiáng)結(jié)果(b) 剖面上部存在一個(gè)花狀構(gòu)造，斷裂發(fā)育，剖面下部斷裂多為小斷裂

圖7 F油田的斷層相似屬性水平切片(a)及其斷層增強(qiáng)結(jié)果(b)

4 結(jié)論

基于DTW的斷層增強(qiáng)方法不僅適用于簡(jiǎn)單組合斷層的增強(qiáng)處理，而且也適用于復(fù)雜組合斷層的增強(qiáng)處理。模型試算和實(shí)際地震資料處理結(jié)果表明，該方法的斷層增強(qiáng)效果較好，具有實(shí)際應(yīng)用價(jià)值。

(1)根據(jù)DTW方法能穩(wěn)定、高效地從兩個(gè)序列的誤差矩陣中計(jì)算出最小誤差路徑的特點(diǎn)，將其用于計(jì)算斷層圖像的最優(yōu)斷層線，不受噪聲等非斷層信息的影響，具有較好的抗噪性和魯棒性。

(2)運(yùn)用非線性平滑、正向積累和反向追蹤等方法有效提高了最優(yōu)斷層線的求解精度。

(3)采用組合成分濾波器[11]對(duì)斷層屬性圖像濾波，能較好地去除非斷層信息。然后采用非極大值抑制篩選fast特征點(diǎn)的方法選擇種子點(diǎn)，能保證種子點(diǎn)均勻分布于斷層上，可有效減少種子點(diǎn)數(shù)量，提高算法效率。

(4)基于斷層在平面上具有線性特征、在空間上具有平面特征的假設(shè)，本文使用DTW方法計(jì)算最優(yōu)斷層線。經(jīng)過算法處理，完全去除了斷層屬性圖像中的非斷層信息，使斷層更連續(xù)、清晰，斷層分布、延伸特征更明顯，同時(shí)去除了地震屬性圖像中的一些非線性局部特征。本文算法旨在增強(qiáng)地震屬性圖像中的斷層信息，因此去除非斷層信息并不影響最終的斷層增強(qiáng)結(jié)果。