一種動態(tài)匹配特征子波拾取地震同相軸的技術

彭仁艷，徐振旺，劉文峰，余 鋒，董旭光

1.中國石化石油工程地球物理有限公司科技研發(fā)中心，江蘇 南京 210000

2.中國石油遼河油田分公司勘探開發(fā)研究院，遼寧 盤錦 124010

3.中國石化石油工程地球物理有限公司華東分公司，江蘇 南京 210000

4.中國石化石油工程地球物理有限公司國際業(yè)務發(fā)展中心，江蘇 南京 210000

引言

地震同相軸自動識別與拾取在速度分析、資料解釋等方面具有重要的應用價值，但對于大規(guī)模三維地震數(shù)據(jù)，人工拾取工作量巨大，效率低下，無法滿足工業(yè)生產的需要[1-8]。在地震剖面上自動識別與拾取同相軸方面，前人已進行了大量研究，基本思想為基于數(shù)據(jù)的空間分布特征，利用局部相似性尋找線性相關的同相軸[9-14]。根據(jù)計算局部相似性算法的不同對拾取方法進行分類，目前，工業(yè)界中使用最廣泛的是基于傾角信息的拾取算法，即根據(jù)地震數(shù)據(jù)的空間局部相似性，估計整個數(shù)據(jù)空間的傾角場，然后，根據(jù)參考道或種子點，沿著估計的傾角方向依次追蹤相鄰道，最后，將拾取的點進行模式分類，構造出相應的反射層位。其中，描述局部相關性，計算傾角的方法是整個自動拾取的關鍵步驟，其計算方法也有很多種[15]，如利用瞬時波數(shù)方向估計反射層的傾角[16]、利用平面波分解或預測計算局部空間的傾角[17]、利用結構張量的方法估計傾角[18]等。平面波分解假設同相軸是局部線性的，計算精度不高；結構張量對數(shù)據(jù)比較敏感，易受噪音的干擾。另外一種比較常用的自動拾取方法是基于多項式擬合的拾取算法[19]，此類方法大多用在成像道集上的剩余時間（深度）差拾取，通過構造多項式系數(shù)擬合剩余曲率曲線，來完成拾取的工作，這種方法基于多項式假設，不能適用于任意地震剖面拾取情況。

動態(tài)波形匹配是廣泛用于聲波信號匹配、字典匹配等領域的自動匹配算法，它提供了良好的求取兩列信號最佳匹配的方法。動態(tài)波形匹配沒有考慮地震信號的特點，本文結合地震信號的子波特征，對地震信號進行矢量化特征表達，在此基礎上結合動態(tài)波形匹配算法，提出一種同相軸自動識別與拾取方法。介紹了動態(tài)波形匹配的基本原理，在此基礎上，將子波特征矢量化的數(shù)據(jù)利用動態(tài)波形匹配算法進行層位自動追蹤，實現(xiàn)基于子波特征的動態(tài)波形匹配同相軸自動拾取。給出了理論數(shù)據(jù)測試實例和實際資料測試實例，證明了方法的正確性和實用性。

1 動態(tài)波形匹配法基本原理

動態(tài)波形匹配是互相關運算的推廣，是波形特征提取的重要算法，早期在聲波匹配、字典對比等方面應用較為普遍。近年來，地球物理學家將該算法推廣到地震數(shù)據(jù)處理中，并提出動態(tài)時間規(guī)整方法，在去除地震數(shù)據(jù)的靜校正量、校平地震道集及匹配縱橫波成像深度差等方面均得到成功應用[20-25]。動態(tài)波形匹配通過計算波形匹配最短路徑，把波形匹配最短路徑的大小作為相鄰地震道波形之間差異大小程度的度量。

動態(tài)波形匹配的目的是找到一個時移量的序列，使得兩個序列f和g 相匹配后的誤差和達到最小。令f為參考道，g 與f相匹配，定義eili表示對時間序列g 的第i個采樣點進行時移li后與時間序列f第i個采樣點之間的距離，其表達式為

動態(tài)時間規(guī)整DTW 算法通過動態(tài)規(guī)劃法進行求解，分為累積和回溯兩步計算。

累積過程中，利用校正誤差eili遞歸地計算距離累積誤差dili，假設時間應變的約束條件為|ui?ui?1|≤1，則累積過程為

在累積過程的最后一步，即i=N?1 時，計算出所有時差l對應的距離函數(shù)后，找到最小的距離，就可以得到最終的最小總距離函數(shù)

通過往前遞歸地回溯，最后可以找到每一個采樣點對應的時差，即最短路徑。

2 基于子波特征表達的動態(tài)波形匹配地震剖面自動拾取

動態(tài)波形匹配通過定義點與點之間距離，求解點與點之間的最短路徑。該方法沒有考慮地震信號的特征。根據(jù)褶積模型，地震數(shù)據(jù)可以視為反射系數(shù)序列與子波的褶積加隨機噪音組成，即

地震剖面拾取就是識別主要的反射同相軸，即提取主要反射系數(shù)序列。地震子波及噪音的存在對自動拾取帶來較大的干擾，增加了自動拾取的難度。地震子波數(shù)學上就是一小段波動序列，由于薄層調諧效應、吸收衰減以及噪音的影響，在實際地震剖面中子波發(fā)生較大變化。為了減少子波對自動拾取的影響，對地震子波進行特征提取，將地震子波抽象為由波峰位置、波寬、振幅及特征斜率（主瓣與副瓣幅值比）等屬性組成的序列。抽象過程表示為矢量形式為

式中：p--子波波峰/波谷極值位置；q--子波波寬（子波持續(xù)時間/深度范圍）；a--子波波峰/波谷極值大??；v--子波特征斜率；iii，jjj，kkk，lll--單位向量。

式（9）為子波特征矢量。對地震剖面s(t) 進行子波提取，并用子波特征矢量表示

式中：

s(t)--長度為t的地震剖面；ccctc長度為tc的特征矢量剖面數(shù)據(jù)，即選取適當振幅（波峰/波谷）閾值，通過子波特征矢量表達將長度為t的地震剖面轉化為長度為tc的矢量數(shù)據(jù)。

式（10）將地震數(shù)據(jù)轉化為子波特征矢量數(shù)據(jù)，將一維地震道轉化為四維特征序列數(shù)據(jù)。需要指出的是，式（10）中tc通常遠小于t，雖然將一維數(shù)據(jù)擴展為四維，但是轉換后的數(shù)據(jù)規(guī)模得到大大壓縮，能夠顯著提高數(shù)據(jù)的后續(xù)處理效率。子波特征矢量表達對地震剖面進行了特征的提取，降低了剖面中子波以及噪音的影響，對地震同相軸的識別和提取具有重要作用。

為了進行地震同相軸的自動拾取，需要對特征表達后的地震數(shù)據(jù)進行自動追蹤。定義emn,klnk為特征表達矢量數(shù)據(jù)序列cmn和cm+k,n+lnk之間的距離，根據(jù)矢量距離公式，定義為

式中：

ωp，ωq，ωa，ωvp,q,a,v 相應的加權系數(shù)；m--道號；j對應道的樣點數(shù)；k--道號移動量；lnk采樣點移動量，通常，|k|=1 表示相鄰兩道進行波形匹配；|l|＜?--相鄰子波的距離，根據(jù)不同的數(shù)據(jù)特征進行選擇，一般大于同一地震道相鄰子波的最大間距。

在實際處理中，會逐地震道依次迭代處理。不失一般性，后面的公式中忽略m，并令k=1，emn,klnk簡化為en,ln。

用動態(tài)規(guī)劃法求解滿足en,ln的最短路徑。首先，遞歸地計算距離累積誤差dn,ln，誤差積累過程，表達式如下

式中：ε--誤差閾值。

當兩個子波之間的誤差小于閾值ε 時，認為該兩個子波處在同一個同相軸上，令其累積誤差最小，就可以把這兩個子波先匹配起來。當輸入數(shù)據(jù)信噪比低，子波特征變化不明顯時，可適當減小ε 值。對于每一個參考位置子波，都存在一個最短路徑找到另一道上與其波形特征最為逼近的子波，因此，這里把一定動態(tài)范圍內的匹配距離存儲起來，以便于后續(xù)的誤差回代。

下一步為誤差回代，找到一個時移量的序列，使得兩個序列相匹配后的誤差和達到最小，找到最短路徑。

式中：

L子波特征化后地震道樣點數(shù)。

至此，兩個序列中的子波動態(tài)匹配工作已經完成。對整個地震剖面，依次進行序列的兩兩匹配，把相互匹配的子波連接起來，就構成了同相軸，從而實現(xiàn)地震剖面的自動拾取。

3 數(shù)值實驗

3.1 性能測試

首先，采用理論數(shù)據(jù)測試本文方法的正確性以及抗噪性；圖1a 為理論地震剖面，構造一高斯子波，使其在橫向上沿正弦函數(shù)分布，縱向上子波間的間隔為50 個網(wǎng)格，整個剖面共有11 個平行同相軸。圖1b 是利用本文方法在圖1a 所示數(shù)據(jù)上自動拾取的結果（圖中紅色圓點）?？梢钥闯?，圖1b 中的層位均是由同一個子波構造而成，信噪比很高，非常易于拾取，當前大部分自動拾取方法都能做到。

圖1 合成地震剖面及自動拾取結果覆蓋于地震剖面上Fig.1 The synthetic seismic section and the automatic picking result covering on the seismic section

在圖1a 的基礎上，對其添加高斯白噪聲，所使用的信噪比公式為

式中：

r(t)--原始有效信號；n(t)--高斯白噪聲；T地震數(shù)據(jù)縱向采樣點數(shù)。

按式（14）添加信噪比，如圖2a 所示，加入隨機噪音后的信噪比為5 dB。由于噪音的影響，剖面變得更加模糊，層位也不如圖1a 明顯，但是依然有明顯地震同相軸的模式。采用本文自動拾取算法，選取振幅閾值為0.8（整個地震剖面最大振幅的0.8倍）、縱向動態(tài)匹配范圍10 個采樣點，作用于圖2a的數(shù)據(jù)上，拾取結果見圖2b，紅色圓點表示拾取到的地震同相軸。拾取的同相軸總體上正弦函數(shù)變化的趨勢不變，局部在縱向上隨機抖動。這是由于噪音的影響使得子波的形狀也變得上下抖動，波峰極值的位置發(fā)生偏移，拾取的結果才會出現(xiàn)圖2b 中所示的情況。進一步降低信噪比，如圖3a 所示為信噪比1 dB 的合成地震剖面，由于信噪比更低，剖面變得更加模糊，相同參數(shù)下自動拾取的地震同相軸的連續(xù)性自然不如圖2b，但是依然能夠看到明顯的同相軸趨勢。對比圖1b、圖2b 和圖3b 中的拾取結果，可見本文自動拾取算法具有一定的抗噪能力。

對圖1a 中高斯子波的位置（時間）和相位做隨機擾動，最大擾動范圍為8 個網(wǎng)格，擾動結果見圖4a，相鄰道的層位隨機錯動，相應的子波的相位也發(fā)生變化，這都給自動拾取帶來了一定難度。在圖4a 上采用本文方法進行自動拾取，選取振幅閾值和動態(tài)匹配范圍與圖2a 一致，拾取結果如圖4b 中紅色圓點所示，所有的層位都被拾取出來了，因為層位的錯動拾取結果不如圖1b 中那樣光滑，層位線段表現(xiàn)的曲折不平，局部出現(xiàn)較大的抖動。在此基礎上，加大最大擾動范圍為16 個網(wǎng)格，對于復雜的實際數(shù)據(jù)來說，同一層位的相鄰點錯動距離也很少達到16 個網(wǎng)格，一般情況下很少出現(xiàn)這種劇烈的陡構造。擾動結果見圖5a，相鄰道的層位隨機錯動，子波的相位也會發(fā)生更大變化，使得自動拾取難度更大。同樣參數(shù)設置下，使用本文算法拾取圖5a中的地震同相軸，拾取結果如圖5b 中紅色圓點所示，可以看出，所有的層位都被拾取出來了，因為層位的錯動更加嚴重，拾取結果不如圖4b 中那樣光滑，層位線段嚴重的曲折不平，局部位置出現(xiàn)較大的抖動，甚至出現(xiàn)了尖刺和溝谷。通過這兩組實驗說明，本文的方法對于子波位置和相位擾動具有一定的適應性，而且在實際數(shù)據(jù)中相鄰點出現(xiàn)16 個網(wǎng)格的錯動是比較罕見的，因此，即使在復雜情形下，本文拾取算法的精度仍然能夠得到滿足。

圖2 5 dB 含噪地震剖面及自動拾取結果覆蓋于地震剖面上Fig.2 The noisy seismic section(5 dB)and the automatic picking result covering on the seismic section

圖3 1 dB 含噪地震剖面及自動拾取結果覆蓋于地震剖面上Fig.3 The noisy seismic section(1 dB)and the automatic picking result covering on the seismic section

圖4 8 個網(wǎng)格時移地震剖面及自動拾取結果覆蓋于地震剖面上Fig.4 The time-lapse seismic profile(8 grids)and automatic picking results covering on the seismic section

圖5 16 個網(wǎng)格時移地震剖面及自動拾取結果覆蓋于地震剖面上Fig.5 The time-lapse seismic profile(16 grids)and automatic picking results covering on the seismic section

前文只考慮了子波位置和相位隨機擾動和隨機噪音的單一影響，將二者同時考慮進行拾取方法的穩(wěn)定性測試也是很有必要的。在地震數(shù)據(jù)處理中，位置擾動對應著子波道間時差的變化，噪音大小則決定了數(shù)據(jù)的質量（信噪比和分辨率），二者是影響拾取結果好壞的關鍵因素。

在圖1a 的基礎上，對高斯子波的位置（時間）沿縱向做隨機的擾動，如圖6a 所示，最大的擾動范圍為10 個網(wǎng)格，同時加入隨機噪音，信噪比為5 dB。將本文的自動拾取算法應用在圖6a 的合成地震剖面上，參數(shù)（振幅閾值和動態(tài)匹配范圍）與前面保持一致，拾取結果如圖6b 所示，紅色圓點表示拾取到的地震同相軸。拾取的同相軸總體上保持正弦函數(shù)變化的趨勢，局部存在縱向隨機抖動，同相軸拾取結果有一定的道間時差，可見地震剖面在隨機擾動的同時加上隨機噪音的影響，使得層位曲線非常粗糙，波峰極值的位置發(fā)生偏移，拾取的結果才會出現(xiàn)圖6b 中所示的情況。在此隨機擾動的基礎上，如圖7a 所示，降低地震剖面的信噪比到1 dB。拾取結果展示在圖7b 中，紅色圓點表示拾取到的同相軸，相比圖6b 結果，發(fā)現(xiàn)拾取結果同相軸連續(xù)降低，是因為隨機噪音和位置擾動使得子波難以維持其原有的形狀和變化趨勢，子波能量分布不均衡，基于子波波形特征的自動拾取方法對子波變化比較敏感，如遇到子波變形過于劇烈的同相軸，則其拾取精度會下降。

圖6 5 dB+10 個網(wǎng)格含噪+時移地震剖面及自動拾取結果覆蓋于地震剖面上Fig.6 Noisy&time-lapse seismic section(5 dB+10 grids)and automatic picking results covering on the seismic section

圖7 1 dB+10 個網(wǎng)格含噪+時移地震剖面及自動拾取結果覆蓋于地震剖面上Fig.7 Noisy&time-lapse seismic section(1 dB+10 grids)and automatic picking results covering on the seismic section

3.2 實際資料測試

利用東部某探區(qū)實際地震成像剖面測試本文自動拾取方法的實用性。針對該探區(qū)地震剖面，選取振幅閾值為0.6、縱向動態(tài)匹配范圍10 個采樣點，圖8a、圖8b 分別為該探區(qū)Inline 測線方向和Crossline 測線方向的成像剖面以及本文算法自動拾取效果的展示（圖中紅色圓點所示）。從圖中可以看出，通過給定合適的參數(shù)，文中拾取算法能夠較為準確地識別出地震同相軸位置，為后續(xù)基于成像道集的層析成像、層位解釋提供了一定幫助。

圖8 東部某探區(qū)實際地震剖面同相軸自動拾取結果Fig.8 Automatic picking up result of in-phase axis of actual seismic section in an exploration area in the east

4 討論

通過數(shù)值實驗分析，針對低信噪比和道間時差地震數(shù)據(jù)，本文算法仍然可以實現(xiàn)地震同相軸的自動拾取工作。根據(jù)子波特征矢量的定義和算法原理，整個自動拾取過程中，振幅閾值和動態(tài)匹配范圍兩個參數(shù)需要注意。振幅閾值越大，子波特征表達后，保留能量越大的特征子波，因此，自動拾取結果越稀疏；而動態(tài)匹配范圍，則是控制匹配追蹤時縱向掃描的范圍影響自動拾取的精度與效率，特別是如果引入種子點后，拾取特定地震同相軸時，動態(tài)匹配范圍過大，可能會導致錯誤的拾取結果。

本文主要是給出了一種地震同相軸自動拾取的方法，并用于地震剖面同相軸的自動拾取中，其目的是為成像道集層析反演提供層位信息。而本文算法在地震處理解釋中的具體應用理論上應該還有一定的擴展空間。比如，引入種子點概念后，可以用于特定反射界面的追蹤、地震初至信息的獲取、角道集RMO 拾取等等。但是，難以滿足對特定斷層的解釋與識別工作。

5 結論

（1）基于褶積模型地震數(shù)據(jù)可以表達為反射系數(shù)序列和子波的褶積，地震同相軸拾取其實質是拾取主要反射系數(shù)序列的過程。

（2）子波以及噪音存在對地震同相軸的自動拾取帶來很大的干擾，影響自動拾取的精度?；谧硬ㄐ螒B(tài)特征表達，將地震數(shù)據(jù)轉化為子波特征矢量數(shù)據(jù)，減弱了子波以及噪音對自動拾取的影響。通過定義矢量距離，利用動態(tài)規(guī)劃法計算道與道之間的最小距離，從而實現(xiàn)的地震同相軸的自動拾取。同時通過特征化子波，可以對地震數(shù)據(jù)進行壓縮，有利于提升海量地震數(shù)據(jù)自動拾取的計算效率。

（3）理論和實際數(shù)據(jù)的測試結果表明，本文算法可以適應低信噪比和存在道間時差的地震數(shù)據(jù)。但是，本文未具體討論考慮種子點后算法的應用情況，這也是今后進一步研究需要探討的問題。