基于動(dòng)態(tài)時(shí)間規(guī)整算法的縱波與轉(zhuǎn)換波時(shí)間域匹配

張一琛 ，陳雙全 *，靳松 ，李向陽(yáng) ,

1 中國(guó)石油大學(xué)(北京)油氣資源與探測(cè)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 102249 2 中國(guó)石油大學(xué)(北京)CNPC物探重點(diǎn)實(shí)驗(yàn)室，北京 102249 3英國(guó)聯(lián)邦地質(zhì)調(diào)查局(BGS)，愛(ài)丁堡 EH14 4AP

0 引言

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

早期的縱波與轉(zhuǎn)換波數(shù)據(jù)的匹配主要是依靠人工層位標(biāo)定法。解釋人員結(jié)合測(cè)井資料進(jìn)行地震剖面的層位標(biāo)定，通過(guò)對(duì)比相似的反射層位來(lái)匹配同相軸。然而該方法具有極大的局限性，僅能適用于地層反射特征明顯的區(qū)域，一旦地層構(gòu)造復(fù)雜，反射波特征不明顯時(shí)，匹配結(jié)果就會(huì)出現(xiàn)較大的隨意性和誤差。

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

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

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

1 縱波與轉(zhuǎn)換波時(shí)間域匹配原理

如圖1所示，對(duì)于地下界面同一反射點(diǎn)的縱波與轉(zhuǎn)換波數(shù)據(jù)，由于波傳播的路徑和傳播的時(shí)間不同，觀測(cè)到的同一反射點(diǎn)的地震波信息在縱波數(shù)據(jù)與轉(zhuǎn)換波數(shù)據(jù)中，不僅時(shí)間不同，地震波的波形特征也會(huì)產(chǎn)生差異。

S表示激發(fā)的地震波的炮點(diǎn)，P為來(lái)自同一反射點(diǎn)C的縱波數(shù)據(jù)記錄點(diǎn)，V為轉(zhuǎn)換波記錄點(diǎn)，O點(diǎn)為SP的中點(diǎn)位置，x表示觀測(cè)的偏移距，h為反射界面的深度，VP與VS分別為地層的縱波與橫波速度，α 為縱波入射角度與反射角度，β 為反射的轉(zhuǎn)換波角度。

對(duì)于偏移或疊加處理后的疊加剖面，縱波與轉(zhuǎn)換波數(shù)據(jù)中來(lái)自同一反射界面的反射波時(shí)間關(guān)系可以表示成：

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

縱波的雙程旅行時(shí)

和轉(zhuǎn)換波的雙程旅行時(shí)

設(shè)地層中的縱波與轉(zhuǎn)換波速度比為γ=VP/VS，則有

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

2 動(dòng)態(tài)時(shí)間規(guī)整算法

動(dòng)態(tài)時(shí)間規(guī)整(Dynamic time warping，簡(jiǎn)稱(chēng)DTW)算法是在一定約束下，為兩個(gè)信號(hào)序列尋找對(duì)齊的最優(yōu)化路徑的方法。該方法由Sakoe和Chiba(1978)提出[27]，目前此算法在語(yǔ)音識(shí)別領(lǐng)域應(yīng)用廣泛[28]，其中圖2是最常見(jiàn)的時(shí)間序列信號(hào)匹配應(yīng)用的例子。

而在地震數(shù)據(jù)處理中，Anderson和Gaby(1983)在地球物理的時(shí)間序列應(yīng)用中基于動(dòng)態(tài)時(shí)間規(guī)整算法提出了動(dòng)態(tài)波形匹配算法[29]，之后Liner和Clapp(2004)又提出了動(dòng)態(tài)處理方法來(lái)解決地震記錄各道之間的匹配問(wèn)題[30]。

接下來(lái)，以?xún)蓚€(gè)合成地震記錄f[ i]，g[ i]為例，簡(jiǎn)單介紹DTW算法匹配時(shí)間域的地震信號(hào)的過(guò)程。假設(shè)兩個(gè)地震記錄在時(shí)間域內(nèi)存在一個(gè)時(shí)差u[ i ]，f[ i ]和g[ i ]的關(guān)系近似為

圖2 兩個(gè)時(shí)間序列信號(hào)的對(duì)齊Fig. 2 The alignment of two time series signals

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

2.1 對(duì)齊誤差

縱波與轉(zhuǎn)換波數(shù)據(jù)存在時(shí)移量，對(duì)齊誤差則可以定義為：

式中，i是地震道上的某一個(gè)采樣點(diǎn)數(shù)。l是記錄在u[ i]中的時(shí)移量。這樣就能建立一個(gè)橫坐標(biāo)是轉(zhuǎn)換波數(shù)據(jù)的采樣點(diǎn)數(shù)，縱坐標(biāo)是時(shí)移量的對(duì)齊誤差二維矩陣。

2.2 累積

累積就是對(duì)各個(gè)坐標(biāo)點(diǎn)上的對(duì)齊誤差依次迭代累積求和，得到累積誤差d，公式如下：

由此得到各個(gè)坐標(biāo)點(diǎn)上的累積誤差。之后在累積誤差的基礎(chǔ)上回溯最優(yōu)化路徑，從而找到轉(zhuǎn)換波數(shù)據(jù)地震道上每個(gè)采樣點(diǎn)對(duì)應(yīng)的時(shí)移量序列u[ i ]。

2.3 回溯

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

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

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

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

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

3 實(shí)際數(shù)據(jù)處理流程

本文基于動(dòng)態(tài)時(shí)間規(guī)整算法，形成了縱波與轉(zhuǎn)換波數(shù)據(jù)的時(shí)間域匹配方法。針對(duì)實(shí)際數(shù)據(jù)處理過(guò)程，建立了如圖3所示的數(shù)據(jù)處理流程：

(1)按照層位控制縱波與轉(zhuǎn)換波數(shù)據(jù)匹配的思路，通過(guò)設(shè)定初始縱橫波速度比，根據(jù)公式(3)將轉(zhuǎn)換波數(shù)據(jù)轉(zhuǎn)換到縱波時(shí)間域，對(duì)轉(zhuǎn)換波數(shù)據(jù)進(jìn)行預(yù)匹配處理。

(2)tP表示縱波地震剖面的雙程旅行時(shí)，tC表示相應(yīng)的轉(zhuǎn)換波地震剖面的雙程旅行時(shí)，tP和tC之間的關(guān)系可以記為：

圖3 實(shí)際數(shù)據(jù)處理流程圖Fig. 3 Registration processing fl ow chart for real fi eld seismic data

u(tP)表示時(shí)移序列。利用預(yù)匹配處理后的轉(zhuǎn)換波數(shù)據(jù)與縱波數(shù)據(jù)，采用動(dòng)態(tài)時(shí)間規(guī)整算法求取縱波與轉(zhuǎn)換波數(shù)據(jù)之間的時(shí)移量序列u(tP)。

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

(4)計(jì)算縱波數(shù)據(jù)與轉(zhuǎn)換波數(shù)據(jù)之間的相關(guān)譜，驗(yàn)證最終的時(shí)間域匹配結(jié)果。

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

4 實(shí)際數(shù)據(jù)應(yīng)用研究

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

4.1 轉(zhuǎn)換波數(shù)據(jù)預(yù)匹配處理

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

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

在圖5中可以發(fā)現(xiàn)，轉(zhuǎn)換波剖面經(jīng)過(guò)最初設(shè)置的速度比γ=2的線(xiàn)性壓縮后，圖5a和圖5b中從上到下的主要層位依次大致對(duì)應(yīng)，在時(shí)間域上有較小的時(shí)移量偏差。所以，從預(yù)處理線(xiàn)性壓縮后的結(jié)果中能夠初步辨識(shí)出二者大致對(duì)應(yīng)的地質(zhì)層位。在預(yù)處理中，預(yù)先設(shè)定速度比并進(jìn)行后續(xù)的線(xiàn)性壓縮，其目的是為了大體上劃定范圍，將縱波與轉(zhuǎn)換波剖面初步對(duì)應(yīng)起來(lái)，以便之后應(yīng)用動(dòng)態(tài)時(shí)間規(guī)整方法時(shí)，避免大范圍地掃描時(shí)移量而導(dǎo)致時(shí)移變化過(guò)大的情況發(fā)生。

4.2 應(yīng)用DTW求取時(shí)移量進(jìn)行時(shí)間域匹配

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

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

4.3 匹配結(jié)果分析

為了評(píng)估利用本文提出的方法和數(shù)據(jù)處理流程處理后的匹配效果，應(yīng)用互相關(guān)法驗(yàn)證動(dòng)態(tài)時(shí)間規(guī)整算法在轉(zhuǎn)換波剖面與縱波剖面匹配中的實(shí)際效果。圖8

是本文方法匹配前后的地震道數(shù)據(jù)的互相關(guān)時(shí)延譜。圖8a是沒(méi)有應(yīng)用動(dòng)態(tài)時(shí)間規(guī)整算法的縱波和轉(zhuǎn)換波的互相關(guān)時(shí)延譜，圖8b是動(dòng)態(tài)時(shí)間規(guī)整算法校正后的轉(zhuǎn)換波數(shù)據(jù)與縱波數(shù)據(jù)的互相關(guān)時(shí)延譜。

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

圖6 由動(dòng)態(tài)時(shí)間規(guī)整算法計(jì)算的預(yù)處理后資料的時(shí)移量(a)和原始數(shù)據(jù)的總時(shí)移量(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)預(yù)處理后的縱波地震剖面和(b)用DTW算法計(jì)算的時(shí)移量匹配后的轉(zhuǎn)換波地震剖面對(duì)比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 動(dòng)態(tài)時(shí)間規(guī)整匹配前后的地震道數(shù)據(jù)的互相關(guān)時(shí)延譜比較.(a)應(yīng)用動(dòng)態(tài)時(shí)間規(guī)整匹配前的縱波與轉(zhuǎn)換波數(shù)據(jù)的互相關(guān)時(shí)延譜；(b)應(yīng)用動(dòng)態(tài)時(shí)間規(guī)整匹配后的縱波與轉(zhuǎn)換波數(shù)據(jù)的互相關(guān)時(shí)延譜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中較大的相關(guān)系數(shù)分散在各個(gè)時(shí)移量的位置上，說(shuō)明動(dòng)態(tài)時(shí)間規(guī)整算法前轉(zhuǎn)換波剖面和縱波剖面沒(méi)有達(dá)到匹配的效果，平均相關(guān)系數(shù)為-0.06，相關(guān)性極差。圖8b中，在整個(gè)時(shí)間段上，最深的紅色集中在零時(shí)延處，即零時(shí)延處幾乎都達(dá)到較大的相關(guān)系數(shù)，整體相關(guān)系數(shù)大多分布在0.68左右，表明利用動(dòng)態(tài)時(shí)間規(guī)整算法匹配后的轉(zhuǎn)換波地震剖面與縱波剖面具有很好的一致性，匹配效果良好。其他時(shí)延位置較大的相關(guān)性的產(chǎn)生是由于地震子波的影響所致。比較動(dòng)態(tài)規(guī)整時(shí)間域匹配前后的互相時(shí)延譜，可以定量確定，經(jīng)過(guò)規(guī)整后的轉(zhuǎn)換波地震剖面與縱波地震剖面的一致性得到極大的提升。

5 討論與結(jié)論

本文利用非線(xiàn)性反演的DTW方法，直接從縱波與轉(zhuǎn)換波數(shù)據(jù)之間的時(shí)移量入手，進(jìn)行反演求取整個(gè)數(shù)據(jù)的時(shí)移量序列，完成了縱波與轉(zhuǎn)換波數(shù)據(jù)的時(shí)間域匹配。經(jīng)實(shí)際多波地震數(shù)據(jù)的應(yīng)用，結(jié)果表明，文中提出的縱波與轉(zhuǎn)換波數(shù)據(jù)時(shí)間域匹配方法和技術(shù)流程，可以很好地解決縱波與轉(zhuǎn)換波數(shù)據(jù)匹配過(guò)程中的非線(xiàn)性問(wèn)題，而且，匹配過(guò)程中不需要進(jìn)行層位的控制或約束，就可以達(dá)到很好的匹配效果。眾所周知，PP波與PS波的匹配是一個(gè)系統(tǒng)的復(fù)雜工程，不僅包括二者的運(yùn)動(dòng)學(xué)特征匹配，還包括動(dòng)力學(xué)特征匹配與綜合解釋等。本文所展示的是如何在時(shí)間域完成匹配，解決二者的運(yùn)動(dòng)學(xué)非線(xiàn)性?xún)?yōu)化問(wèn)題，還未考慮PP波與PS波的動(dòng)力學(xué)特征匹配過(guò)程，這是本文下一步研究的方向。

本文取得如下幾點(diǎn)認(rèn)識(shí)：

(1)經(jīng)典的動(dòng)態(tài)時(shí)間規(guī)整算法需要做小小的改動(dòng)引入時(shí)移量，便能夠利用非線(xiàn)性全局最優(yōu)化的思想找到縱波和轉(zhuǎn)換波匹配的最佳路徑。

(2)通過(guò)設(shè)定初始的縱橫波速度比對(duì)縱波和轉(zhuǎn)換波資料進(jìn)行預(yù)處理，可以在整體上有效地找到二者對(duì)應(yīng)的標(biāo)志性層位，為之后應(yīng)用動(dòng)態(tài)時(shí)間規(guī)整算法更精細(xì)地依靠波形相似性計(jì)算時(shí)移量奠定基礎(chǔ)。該方法的逐步精細(xì)化的研究思路具有很好的推廣性。

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

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

[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] 陳雙全, 王尚旭, 李向陽(yáng), 等. 與 PP 波的時(shí)間匹配[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.