基于動態時間規整算法的縱波與轉換波時間域匹配

張一琛 ，陳雙全 *，靳松 ，李向陽 ,

1 中國石油大學(北京)油氣資源與探測國家重點實驗室，北京 102249 2 中國石油大學(北京)CNPC物探重點實驗室，北京 102249 3英國聯邦地質調查局(BGS)，愛丁堡 EH14 4AP

0 引言

多波地震勘探技術能夠采集到更加豐富的地震波場信息，并能利用多波地震數據的走時、振幅、以及縱、橫波之間的時差、速度比和其他各向異性參數，對油氣藏的巖石物性和流體性質進行更為精確的成像描述，從而最大限度地消除單一縱波儲層預測的多解性，提高非構造油氣藏探測精度，有助于研究穿透氣云的成像和預測裂縫的發育程度[1-2]。目前，多波地震勘探主要是利用轉換波數據與縱波數據進行聯合反演與解釋，開展油藏的儲層描述與油氣檢測。其中多波地震勘探的一個重要環節是縱波與轉換波數據在不同時間域的匹配處理，這是縱波與轉換波資料進行聯合反演與解釋的關鍵一步。縱波與轉換波數據匹配效果的好壞和精度的高低將直接影響后續多波聯合反演與解釋的效果。

早期的縱波與轉換波數據的匹配主要是依靠人工層位標定法。解釋人員結合測井資料進行地震剖面的層位標定，通過對比相似的反射層位來匹配同相軸。然而該方法具有極大的局限性，僅能適用于地層反射特征明顯的區域，一旦地層構造復雜，反射波特征不明顯時，匹配結果就會出現較大的隨意性和誤差。

之后，為了盡量避免在數據匹配過程中的人為因素的干擾和誤差，研究人員提出了利用縱波與轉換波數據的運動學與動力學特征進行基于數據處理的匹配方法。該類方法的主要思路是先給出初始縱橫波速度比模型，然后結合相關性來校正模型或者利用最小二乘等迭代方法更新速度比模型，之后再憑借精確的速度比模型將轉換波匹配到縱波時間域。Gaiser(1996)利用縱橫波速度比(γ0)和最大相關法[3]進行多波資料匹配，Chan(1998)對常數γ0值多次試驗修正后實現時間對數域縱橫波匹配[4]，該匹配方法只適應于特定目的層。Van Dok和Gaiser等(2001)利用最大相似性原理[5]，掃描縱波與轉換波資料的γ0譜，拾取平均γ0值來匹配圖像。Garotta等(2002)通過相關求取縱橫波速度比[6]，在時間域匹配縱波與轉換波。Fomel等(2003)利用最小二乘法優化全局最優解來實現縱橫波自動匹配[7]。Van Dok和Kristiansen(2003)通過手動拾取具有明顯地質特征的層位實現時間域的大致匹配[8]。Nahm等(2003)通過對準主要探區斷層瞬時相位時間切片來進行數據的匹配[9]。Nickel and Sonneland(2004)采用最小二乘法經過多次快速迭代實現自動同相軸匹配[10]。Xu和Stewart(2005)利用最大相關理論[11]，在大套地層對準的原則上實現縱波和轉換波的時間匹配。Fomel等(2005)利用最小二乘法，振幅和頻率平衡匹配法，剩余γ0掃描法，和時間匹配微調的方法[12]來實現縱橫波的匹配。Fomel和Jin(2009)介紹了一種基于兩種數據局部相似性的方法進行時移地震數據匹配[13]。Wang等(2009)對比縱波和轉換波的同層反射，輔以相似性最大原理求取γ0來實現數據匹配[14]。Bansal 等(2009)在縱波和轉換波時間精確匹配的基礎上，利用中、遠偏移距地震道信息，拓展轉換波由于動校正拉伸或固有因素導致衰減的高頻信息[15]，提高分辨率，從而實現二者頻帶的匹配。陳雙全等(2010)利用縱橫波速度比掃描方法校正數據偏差，并根據構造成像的標志性層位確定大致時間及速度范圍選定計算時窗，進行速度比掃描得到最佳速度比算子，從而完成匹配[16]。Zhou等(2010)通過控制層位求取γ0值匹配蘇里格地區的縱波和轉換波剖面，并用全波測井數據進行匹配微調[17]。Chen等(2013)參考VTI介質各向異性特性，實現了轉換波疊前雙參數掃描結合縱橫波速度比模型[18]的疊前匹配方法。然而，上述提出的方法大致都會有兩大問題：(1)相關法中的相關時窗不能適應于變化劇烈的時移量，且計算效率、精度較低。(2)建立初始速度比模型本身就是一個難題，模型若偏差太大，迭代結果也不會理想。

面對上述問題，學者們轉換思路進行了更深入的探索。Yuan等(2008)采用模擬退火算法求得縱波和轉換波的最大相似性[19]，實現二者的時間匹配。之后，Deng等(2011)開發出快速模擬退火算法[20]的相關應用。Hale(2013)實現了利用動態時間規整方法進行多維縱波與轉換波數據的匹配[21]，之后Hale和Compton(2013)改進該方法，利用誤差累積進行時移量的平滑處理，完善了匹配效果[22]。Chen等(2014)利用疊前縱波反演得到的偽橫波和轉換波掃描估算速度比[23]實現匹配。Hyoungsu(2015)采用正則化全變差和分段線性基函數來構造時間規整最優函數[24]，該方法能在厚層計算精確的速度比和避免時差的劇烈變化。Zhang(2017)憑借P波AVO反演出具有較高相關性的屬性，然后利用DTW小范圍匹配[25]。Wang(2017)提出范數擬合線性動態規整函數，再在稀疏性約束下匹配剖面[26]。在研究過程中發現，非線性最優化求解的思路在縱橫波匹配中，相較于其他方法，有著更高的精度和效率，匹配結果更加理想。

本文正是在前人研究之上，為了提高匹配的精度和效率，利用具有非線性最優化特性的動態時間規整算法形成了基于數據驅動的縱波與轉換波數據時間域匹配方法和實際數據處理技術流程。從經典的動態時間規整算法入手，改進完善了算法流程以適用于縱波與轉換波數據的時間域匹配。在實際數據處理中，首先通過設定初始縱橫波速度比，對縱波與轉換波數據預匹配，使縱波與轉換波的主要層位足夠接近，滿足應用動態時間規整算法的條件。然后，對轉換波數據采用動態時間規整算法求取時移量，完成從轉換波時間域匹配到縱波時間域的過程。最后，利用文中提出的縱波與轉換波時間域匹配方法與數據處理流程，對一實際多波地震數據進行了縱波與轉換波數據的匹配處理。

1 縱波與轉換波時間域匹配原理

如圖1所示，對于地下界面同一反射點的縱波與轉換波數據，由于波傳播的路徑和傳播的時間不同，觀測到的同一反射點的地震波信息在縱波數據與轉換波數據中，不僅時間不同，地震波的波形特征也會產生差異。

S表示激發的地震波的炮點，P為來自同一反射點C的縱波數據記錄點，V為轉換波記錄點，O點為SP的中點位置，x表示觀測的偏移距，h為反射界面的深度，VP與VS分別為地層的縱波與橫波速度，α 為縱波入射角度與反射角度，β 為反射的轉換波角度。

對于偏移或疊加處理后的疊加剖面，縱波與轉換波數據中來自同一反射界面的反射波時間關系可以表示成：

圖1 縱波及轉換波數據記錄示意圖Fig. 1 P-wave and converted-wave (C-wave) seismic data geometry chart

縱波的雙程旅行時

和轉換波的雙程旅行時

設地層中的縱波與轉換波速度比為γ=VP/VS，則有

因此，縱波與轉換波數據在疊后剖面中進行的時間域匹配，就是根據公式(3)將轉換波數據的旅行時間域轉換到縱波旅行時間域。由于地震數據中地震了波的存在，使得這一匹配過程成為非線性反演問題，需要采用非線性方法進行處理。

2 動態時間規整算法

動態時間規整(Dynamic time warping，簡稱DTW)算法是在一定約束下，為兩個信號序列尋找對齊的最優化路徑的方法。該方法由Sakoe和Chiba(1978)提出[27]，目前此算法在語音識別領域應用廣泛[28]，其中圖2是最常見的時間序列信號匹配應用的例子。

而在地震數據處理中，Anderson和Gaby(1983)在地球物理的時間序列應用中基于動態時間規整算法提出了動態波形匹配算法[29]，之后Liner和Clapp(2004)又提出了動態處理方法來解決地震記錄各道之間的匹配問題[30]。

接下來，以兩個合成地震記錄f[ i]，g[ i]為例，簡單介紹DTW算法匹配時間域的地震信號的過程。假設兩個地震記錄在時間域內存在一個時差u[ i ]，f[ i ]和g[ i ]的關系近似為

圖2 兩個時間序列信號的對齊Fig. 2 The alignment of two time series signals

在實際地震數據匹配問題中，f[ i ]和g[ i]是指數據中的采樣點數值。DTW算法用于匹配縱波與轉換波數據，需要完成以下三步處理得到時間域的偏移量。

2.1 對齊誤差

縱波與轉換波數據存在時移量，對齊誤差則可以定義為：

式中，i是地震道上的某一個采樣點數。l是記錄在u[ i]中的時移量。這樣就能建立一個橫坐標是轉換波數據的采樣點數，縱坐標是時移量的對齊誤差二維矩陣。

2.2 累積

累積就是對各個坐標點上的對齊誤差依次迭代累積求和，得到累積誤差d，公式如下：

由此得到各個坐標點上的累積誤差。之后在累積誤差的基礎上回溯最優化路徑，從而找到轉換波數據地震道上每個采樣點對應的時移量序列u[ i ]。

2.3 回溯

在完成累積處理后，需要遍歷i=N?1上的所有l，

從而找到最小的d[ N?1,l]，即

最小的d[ N ?1,l ]中的l便是時移序列u[ i]的最后一個數值u[ N?1]。然后迭代尋找前一個時移量，直到最后找到第一個時移量u[0]，公式如下：

需要注意的是，為了防止相鄰采樣點之間的時移量變化太大，需要添加縱向約束條件：

以此保證各采樣點的時移量縱向平滑變化。在此約束下，由后一個時移量迭代求前一個時移量時，就只需要在d[ i? 1,l?1]，d[ i? 1,l]，d[ i? 1,l+ 1]三者中比較找最小d，從而得到相應的時移量。最后，通過回溯最佳路徑就找到了每一個采樣點的時移量，依次記錄在序列u[ i]中。

3 實際數據處理流程

本文基于動態時間規整算法，形成了縱波與轉換波數據的時間域匹配方法。針對實際數據處理過程，建立了如圖3所示的數據處理流程：

(1)按照層位控制縱波與轉換波數據匹配的思路，通過設定初始縱橫波速度比，根據公式(3)將轉換波數據轉換到縱波時間域，對轉換波數據進行預匹配處理。

(2)tP表示縱波地震剖面的雙程旅行時，tC表示相應的轉換波地震剖面的雙程旅行時，tP和tC之間的關系可以記為：

圖3 實際數據處理流程圖Fig. 3 Registration processing fl ow chart for real fi eld seismic data

u(tP)表示時移序列。利用預匹配處理后的轉換波數據與縱波數據，采用動態時間規整算法求取縱波與轉換波數據之間的時移量序列u(tP)。

(3)利用計算得到的時移量序列u(tP)，對轉換波數據進行更新匹配，之后，可以計算得到縱橫波速度比模型。

(4)計算縱波數據與轉換波數據之間的相關譜，驗證最終的時間域匹配結果。

上述的(2)(3)步驟可以重復處理1～2次。

4 實際數據應用研究

在實際數據應用處理中，選用了一條二維陸上多波地震勘探數據中的縱波與轉換波數據，進行該方法與處理流程的應用研究。該縱波與轉換波地震剖面是經過常規數據處理后的疊后地震剖面，如圖4所示，左圖(圖4a)為原始的縱波地震剖面，右圖(圖4b)為原始的轉換波地震剖面?？v波數據的每道采樣點數為1000，時間采樣間隔為4 ms；轉換波數據的每道采樣點數為2000，采樣間隔為2 ms，則原始的縱波與轉換波數據的時間長度均為4 s。在匹配過程中，要將整個4 s時間長度的轉換波資料匹配到縱波時間域。根據前面的實際數據匹配處理流程，該實際數據應用研究的過程主要包括以下三個方面。

4.1 轉換波數據預匹配處理

由于地下介質的縱橫波速度比一般介于1.5到3.5之間，因此我們可以先設定縱橫波速度比γ為2，利用公式(3)建立起縱波和轉換波時間域的射線傳播雙程旅行時的對應關系，將轉換波轉換到縱波時間域。由于相同地下介質中橫波傳播速度要小于縱波速度，轉換波剖面的地質層位僅對應縱波剖面的一部分，因而選取了縱波時間3 s以內的數據進行匹配。另外，該數據的縱波與轉換波的采樣點數、采樣間隔等參數存在一定的差異，因此在預處理中采用地震道內插值和重采樣相結合的方法，使得二者的采樣點數、采樣間隔等參數統一。經過預處理，縱波和轉換波剖面的采樣點數和采樣間隔分別相同，采樣點數為2001，采樣間隔為1 ms。用γ=2進行時間域預轉換后的結果如圖5所示，圖5a為預處理后的縱波地震剖面，圖5b為預處理后的轉換波地震剖面。

圖4 原始縱波剖面(a)和原始轉換波剖面(b)Fig. 4 The original P-wave seismic section (a) and C-wave seismic section (b)

在圖5中可以發現，轉換波剖面經過最初設置的速度比γ=2的線性壓縮后，圖5a和圖5b中從上到下的主要層位依次大致對應，在時間域上有較小的時移量偏差。所以，從預處理線性壓縮后的結果中能夠初步辨識出二者大致對應的地質層位。在預處理中，預先設定速度比并進行后續的線性壓縮，其目的是為了大體上劃定范圍，將縱波與轉換波剖面初步對應起來，以便之后應用動態時間規整方法時，避免大范圍地掃描時移量而導致時移變化過大的情況發生。

4.2 應用DTW求取時移量進行時間域匹配

得到了預處理之后的縱波和轉換波剖面后，就可以應用動態時間規整算法計算轉換波剖面各道上的時移序列u[ i ]了。然后再利用時移序列對轉換波剖面進行時間規整，實現轉換波剖面精準匹配到縱波剖面時間域的結果。圖6a展示的是將動態時間規整算法應用到預處理后的轉換波剖面計算得到的各地震道的時移序列，在約束條件式(9)的限制下，縱向上每一道的時移量變化平緩。圖6b是預處理資料中線性壓縮掉的和DTW算法計算的時移量之和，即是原始數據的總時移量。

應用時移序列校正預處理后的轉換波數據，得到規整后的轉換波地震剖面(圖7b)。圖7對比了匹配后的轉換波地震剖面和縱波地震剖面，圖7a是縱波地震剖面，圖7b是經過時移量規整到縱波時間域后的轉換波地震剖面。對比圖7a和圖7b，可以看出規整后的轉換波地震剖面完全刻畫了PP波成像的主要構造。和縱波地震剖面層位整體一一對齊，二者主要的地質層位在時間域位置幾乎完全相當，匹配結果一致性很好，精度較高，具有很好的可比性，完成了PS波傳播時間到PP波傳播時間的轉化。

4.3 匹配結果分析

為了評估利用本文提出的方法和數據處理流程處理后的匹配效果，應用互相關法驗證動態時間規整算法在轉換波剖面與縱波剖面匹配中的實際效果。圖8

是本文方法匹配前后的地震道數據的互相關時延譜。圖8a是沒有應用動態時間規整算法的縱波和轉換波的互相關時延譜，圖8b是動態時間規整算法校正后的轉換波數據與縱波數據的互相關時延譜。

圖5 預處理后的縱波剖面(a)和預處理后的轉換波剖面(b)Fig. 5 The pre-processing P-wave seismic section (a) and the C-wave seismic section (b)

圖6 由動態時間規整算法計算的預處理后資料的時移量(a)和原始數據的總時移量(b)Fig. 6 The time shifts (a) computed by using DTW method for pre-processing sections and the sum of the time shifts (b) for raw data

圖7 (a)預處理后的縱波地震剖面和(b)用DTW算法計算的時移量匹配后的轉換波地震剖面對比Fig. 7 Comparison of the (a) pre-processing P-wave seismic section and (b)time matching processed C-wave seismic section by using DTW method

圖8 動態時間規整匹配前后的地震道數據的互相關時延譜比較.(a)應用動態時間規整匹配前的縱波與轉換波數據的互相關時延譜；(b)應用動態時間規整匹配后的縱波與轉換波數據的互相關時延譜Fig. 8 Comparison of correlation time shift spectrum computed by using the P- and C-wave seismic traces, (a) before using DTW method processing, (b) after using DTW method processing

圖8a中較大的相關系數分散在各個時移量的位置上，說明動態時間規整算法前轉換波剖面和縱波剖面沒有達到匹配的效果，平均相關系數為-0.06，相關性極差。圖8b中，在整個時間段上，最深的紅色集中在零時延處，即零時延處幾乎都達到較大的相關系數，整體相關系數大多分布在0.68左右，表明利用動態時間規整算法匹配后的轉換波地震剖面與縱波剖面具有很好的一致性，匹配效果良好。其他時延位置較大的相關性的產生是由于地震子波的影響所致。比較動態規整時間域匹配前后的互相時延譜，可以定量確定，經過規整后的轉換波地震剖面與縱波地震剖面的一致性得到極大的提升。

5 討論與結論

本文利用非線性反演的DTW方法，直接從縱波與轉換波數據之間的時移量入手，進行反演求取整個數據的時移量序列，完成了縱波與轉換波數據的時間域匹配。經實際多波地震數據的應用，結果表明，文中提出的縱波與轉換波數據時間域匹配方法和技術流程，可以很好地解決縱波與轉換波數據匹配過程中的非線性問題，而且，匹配過程中不需要進行層位的控制或約束，就可以達到很好的匹配效果。眾所周知，PP波與PS波的匹配是一個系統的復雜工程，不僅包括二者的運動學特征匹配，還包括動力學特征匹配與綜合解釋等。本文所展示的是如何在時間域完成匹配，解決二者的運動學非線性優化問題，還未考慮PP波與PS波的動力學特征匹配過程，這是本文下一步研究的方向。

本文取得如下幾點認識：

(1)經典的動態時間規整算法需要做小小的改動引入時移量，便能夠利用非線性全局最優化的思想找到縱波和轉換波匹配的最佳路徑。

(2)通過設定初始的縱橫波速度比對縱波和轉換波資料進行預處理，可以在整體上有效地找到二者對應的標志性層位，為之后應用動態時間規整算法更精細地依靠波形相似性計算時移量奠定基礎。該方法的逐步精細化的研究思路具有很好的推廣性。

(3)在計算轉換波每一道的時移序列時，設定的時移量變化約束很好地限制了縱向時移量變化范圍，使時移量在縱向上平滑變化。

(4)動態時間規整算法能夠在縱橫波匹配的應用中取得良好的效果。該方法跳出了以往利用人工層位拾取和利用縱橫波速度比結合相關法或迭代法匹配的思路，避免了人為誤差，相關窗口不適應時移量劇烈變化，計算效率較低和初始速度比模型難建立等一系列問題。轉而利用非線性最優化的DTW算法計算的各個采樣點的時移量直接校正轉換波地震剖面，直觀、便捷地實現了縱波與轉換波數據匹配的目的。同時可以利用計算的時移序列推導出地層的縱橫波速度比，為之后的地層速度建模提供了參考。

[1] LI X Y. Fracture detection using azimuthal variation of P-wave move out from orthogonal seismic survey lines[J]. Geophysics, 1999,64(4): 1193-1201.

[2] LI X Y, ZHANG Y G. Seismic reservoir characterization: How can multicomponent data help?[J]. Journal of Geophysics and Engineering, 2011, 8(2): 123-141.

[3] GAISER J E. Multicomponent VP/VScorrelation analysis[J]. Geophysics, 1996, 61(4): 1137-1149.

[4] CHAN W. Analyzing converted-wave seismic data: Statics, interpolation, imaging, and PP correlation[M]. Calgary: University of Calgary, 1998.

[5] VAN D R, GAISER J. Stratigraphic description of the Morrow formation using mode-converted shear waves: Interpretation tools and techniques for three land surveys[J]. The Leading Edge, 2001, 20(9): 1042-1047.

[6] GAROTTA R, GRANGER P Y, DARIU H. Combined interpretation of PP and PS data provides direct access to elastic rock properties[J]. The Leading Edge, 2002, 21(6): 532-535.

[7] FOMEL S, BACKUS M M. Multicomponent seismic data registration by least squares[C]. SEG Technical Program Expanded Abstracts,Dallas, 2003: 781-784.

[8] VAN D R, KRISTIANSEN P. Event registration and VP/VScorrelation analysis in 4C processing[C]. SEG Technical Program Expanded Abstracts, Dallas, 2003: 785-788.

[9] NAHM J W, DUHON M P. Interpretation and practical applications of 4C-3D seismic data, East Cameron gas fi elds, Gulf of Mexico[J].The Leading Edge, 2003, 22(4): 300-309.

[10] NICKEL M, SONNELAND L. Automated PS to PP event registration and estimation of a high-resolution Vp-VSratio volume[C]. SEG Technical Program Expanded Abstracts, Denver, 2004: 869-872.

[11] XU C, STEWART R R. Delineating a sand reservoir using inversion of 3C-3D seismic data[C]. SEG Technical Program Expanded Abstracts, Houston, 2005: 995-998.

[12] FOMEL S, BACKUS M, FOUAD K, et al. A multistep approach to multicomponent seismic image registration with application to a West Texas carbonate reservoir study[C]. SEG Technical Program Expanded Abstracts, Houston, 2005: 1018-1021.

[13] FOMEL S, JIN L. Time-lapse image registration using the local similarity attribute[J]. Geophysics, 2009, 74(2): A7-A11.

[14] WANG J S, HUO L, TIAN C C, et al. Lithology and fluid prediction using multiwave seismic data[C]. SEG Technical Program Expanded Abstracts, Houston, 2009: 1237-1241.

[15] BANSAL R, KHARE V, JENKINSON T, et al. Correction for NMO stretch and differential attenuation in converted-wave data: A key enabling technology for prestack joint inversion of PP and PS data[J]. The Leading Edge, 2009, 28(10): 1182-1190.

[16] 陳雙全, 王尚旭, 李向陽, 等. 與 PP 波的時間匹配[J]. 石油地球物理勘探, 2010(4): 497-501. [CHEN S Q, WANG S X, LI X Y, et al. P-SV wave and PP wave time domain registration. Petroleum Geophysical Exploration, 2010(4): 497-501.]

[17] ZHOU Y J, TAO J Q, DOU Y S, et al. Quantitative interpretation and joint inversion of multicomponent seismic data: Application to the Sulige Gas Field, China[C]. SEG Technical Program Expanded Abstracts, Denver, 2010: 1646-1650.

[18] CHEN S Q, LI X M, LI X Y. Automated time-domain transform of converted waves by prestack double-parameter scanning[J]. Journal of Geophysics and Engineering, 2013, 10(4): 045010.

[19] YUAN J J, NATHAN G, CALVERT A, et al. Automated C-wave registration by simulated annealing[C]. SEG Technical Program Expanded Abstracts, Las Vegas, 2008: 1043-1047.

[20] DENG Z W, SEN M K, WANG U, et al. Prestack PP & PS wave joint stochastic inversion in the same PP time scale[C]. SEG Technical Program Expanded Abstracts, San Antonio, 2011: 1303-1307.

[21] HALE D. Dynamic warping of seismic images[J]. Geophysics, 2013, 78(2): S105-S115.

[22] HALE D, COMPTON S. Smooth dynamic warping[R]. CWP Reposts, 2013.

[23] CHEN S Q, LI X Y, TANG J. Converted-wave time domain registration using the inverted pseudo-PS-wave attribute section[J]. Journal of Geophysics and Engineering, 2014, 11(1): 015007.

[24] BAEK H, KEHO T H. Detection and characterization capabilities of time/amplitude warping[C]. SEG Technical Program Expanded Abstracts, New Orleans, 2015: 5379-5383.

[25] ZHANG T, SUN P, QIAN Z, et al. Automatic registration of PP and PS data based on seismic attributes[C]. 79th EAGE Conference and Exhibition, Paris, 2017.

[26] WANG H, SACCHI M, MA J. Linearized dynamic warping with l1-norm constraint for multi-component registration[J]. Journal of Applied Geophysics, 2017, 139: 170-176.

[27] SAKOE H, CHIBA S. Dynamic programming algorithm optimization for spoken word recognition[J]. IEEE transactions on acoustics,speech, and signal processing, 1978, 26(1): 43-49.

[28] MüLLER M. Information retrieval for music and motion[M]. Heidelberg: Springer, 2007.

[29] ANDERSON K R, GABY J E. Dynamic waveform matching[J]. Information Sciences, 1983, 31(3): 221-242.

[30] LINER C L, CLAPP R G. Nonlinear pairwise alignment of seismic traces[J]. The Leading Edge, 2004, 23(11): 1146-1150.