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

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

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

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

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

引言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

式中：

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

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

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

式中：

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

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

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

式中：ε--誤差閾值。

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

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

式中：

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

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

3 數(shù)值實驗

3.1 性能測試

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

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

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

式中：

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

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

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

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

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

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

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

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

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

圖6 5 dB+10 個網(wǎng)格含噪+時移地震剖面及自動拾取結(jié)果覆蓋于地震剖面上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)格含噪+時移地震剖面及自動拾取結(jié)果覆蓋于地震剖面上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 個采樣點(diǎn)，圖8a、圖8b 分別為該探區(qū)Inline 測線方向和Crossline 測線方向的成像剖面以及本文算法自動拾取效果的展示（圖中紅色圓點(diǎn)所示）。從圖中可以看出，通過給定合適的參數(shù)，文中拾取算法能夠較為準(zhǔn)確地識別出地震同相軸位置，為后續(xù)基于成像道集的層析成像、層位解釋提供了一定幫助。

圖8 東部某探區(qū)實際地震剖面同相軸自動拾取結(jié)果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ù)需要注意。振幅閾值越大，子波特征表達(dá)后，保留能量越大的特征子波，因此，自動拾取結(jié)果越稀疏；而動態(tài)匹配范圍，則是控制匹配追蹤時縱向掃描的范圍影響自動拾取的精度與效率，特別是如果引入種子點(diǎn)后，拾取特定地震同相軸時，動態(tài)匹配范圍過大，可能會導(dǎo)致錯誤的拾取結(jié)果。

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

5 結(jié)論

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

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

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