動態時間規整的斷層增強方法

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

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

0 引言

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

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

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

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

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

1 方法原理

1.1 DTW方法原理

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

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

(1)

DTW求解的目標是找到一個時間序列

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

(2)

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

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

(3)

且滿足

(4)

其中l(t)為f(t)與g(t)對齊時(t=0，1，2，…，N-1)對應對齊誤差矩陣中從點(0，0)到點(N-1，N-1)的任意一條路徑(N為時間序列的長度)，且D[l(t)]滿足

(5)

并服從約束條件

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

(6)

當f(t)與g(t)對齊時，式(6)確保i+T(i)隨著采樣點的增加不會過快增加或減少。

1.2 最優斷層線計算方法

斷層屬性體局部極大值(或極小值)處往往指示斷層，斷層線的脊線由斷層圖像中連續屬性極值點組成。若將斷層屬性圖像看作DTW的對齊誤差矩陣，將DTW求解對齊誤差最小值路徑相應地變為求解最大值路徑，則求取屬性圖像中的最優斷層線變為求取圖像中最大值路徑問題

(7)

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

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

(8)

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

(1)選取斷層上的一個點作為控制點，將斷層圖像中控制點左、右兩側錐形區域的數據置0(圖1c)，確保計算出的最優斷層線從控制點通過。

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

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)為正向積累(從上至下)計算的圖像值，橫坐標i對應圖像的行方向，縱坐標j對應圖像的列方向。

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

圖1 最優斷層線計算結果(a)原始地震圖像； (b)斷層屬性圖像； (c)控制點設置； (d)非線性平滑； (e)計算出的最優斷層線； (f)斷層增強結果

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)

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

2)正向積累

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

3)反向追蹤

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

j(N-1)=l

(14)

j(i-1)=

(15)

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

(3)計算出最優斷層線后，僅保留原始斷層圖像中最優斷層線上的屬性值，再沿斷層線對斷層屬性值平滑處理，即可得到最終斷層增強結果(圖1f)。若經過控制點的最大值路徑記為L(i)=[i，j(i)]，i=0，1，2，…，N-1，則該斷層線上的新屬性值為

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

(16)

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

1.3 斷層增強處理流程

實際斷層屬性圖像中的斷層不能直接用最優斷層線計算方法增強處理。本文的斷層增強處理方法是將復雜斷層分解為局部窗口內的簡單斷層，然后在局部窗口內計算最優斷層線，最后綜合局部最優斷層線得到最終的復雜斷層增強結果。以原始斷層圖像(圖2a)中的斷層為例，說明斷層增強的具體實現步驟。

1.3.1 斷層圖像濾波處理

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

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

組合成分濾波器是一種組合濾波器，相當于構建了一個協方差矩陣

(17)

其中

(18)

(19)

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

1.3.2 種子點選取

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

1.3.3 種子點斷層方向計算

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

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

(20)

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

圖3 斷層方向計算及沿斷層方向選取數據效果(a)斷層方向掃描； (b)沿斷層選取計算數據 沿紅框的局部窗口選擇數據，斷層近于垂直

1.3.4 種子點最優路徑計算

對于每個種子點，本文選取一個以種子點為中心的矩形窗口，在這個窗口內計算通過種子點的最優路徑。假設斷層傾角掃描的方向為u，垂直于斷層傾向為v，在選擇矩形窗口時，使矩形的長邊平行于u。在這樣的局部窗口內斷層是近于垂直的，此時在每個種子點局部窗口內計算斷層線的方法與最優斷層線計算方法一致。

1.3.5 斷層增強處理

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

(21)

2 模型數據試算

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

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

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

由未加噪模型方差屬性(圖4b)的增強結果(圖4c)、相似屬性(圖4d)的增強結果(圖4e)可見，相交斷層的交切關系清楚，在保持斷層趨勢的同時，斷層增強結果更連續、清晰。由加噪模型相似屬性(圖5b)的增強結果(圖5c)可見，完全去除了非斷層信息，相交斷層的交切關系清楚，斷層增強結果更連續、清晰。上述結果說明本文方法的抗噪性較好，具有較強的魯棒性。

此外，對比圖4e和圖5c可見，兩者存在細小的局部差異，局部斷層線相差1個像素點。這是由于噪聲影響相似屬性的局部計算結果所致，說明本文算法的增強結果依賴于屬性體的斷層檢測精度，即屬性體的斷層檢測精度越高，增強結果越接近于斷層的真實形態。

3 實際資料處理

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

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

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

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

4 結論

基于DTW的斷層增強方法不僅適用于簡單組合斷層的增強處理，而且也適用于復雜組合斷層的增強處理。模型試算和實際地震資料處理結果表明，該方法的斷層增強效果較好，具有實際應用價值。

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

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

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

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