VTI介質(zhì)起伏界面混合網(wǎng)格旅行時(shí)線性插值計(jì)算方法

王 琦 李慶春 王芷琪

(①長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院，陜西西安 710054; ②中國(guó)航空油料有限責(zé)任公司西北公司，陜西西安 710082)

1 引言

沉積盆地中頁(yè)巖的各向異性與具有垂直對(duì)稱軸的橫向各向同性介質(zhì)(VTI介質(zhì))相似[1]，而在復(fù)雜的構(gòu)造環(huán)境中，如山麓或推覆構(gòu)造區(qū)域，還需要考慮具有傾斜對(duì)稱軸的VTI介質(zhì)(TTI介質(zhì))[2]。在中國(guó)存在著大量的復(fù)雜構(gòu)造發(fā)育地區(qū)，由于近地表地形的強(qiáng)烈起伏，加之地下構(gòu)造復(fù)雜，使得計(jì)算靜校正量、客觀識(shí)別地震波以及實(shí)現(xiàn)地下目標(biāo)的精確地震成像，都成為具有挑戰(zhàn)性的難題。

VTI介質(zhì)是一種最為常見(jiàn)的各向異性介質(zhì)。Thomsen[3]根據(jù)彈性系數(shù)給出了一套各向異性參數(shù)，并分析了具有弱各向異性的VTI介質(zhì)中地震波的傳播特征，通過(guò)對(duì)相速度計(jì)算公式的近似處理，導(dǎo)出了目前仍在廣泛使用的正常時(shí)差公式；在Thomsen研究工作的基礎(chǔ)上，Sena[4]在弱各向異性線性近似下，分別推導(dǎo)了qP波、qSV波和qSH波群速度隨群角變化的表達(dá)式。

在VTI介質(zhì)波場(chǎng)正演模擬中，射線追蹤方法具有顯示直觀、計(jì)算效率高以及對(duì)模型有較高適應(yīng)能力的優(yōu)點(diǎn)。目前眾多學(xué)者已經(jīng)將各向同性介質(zhì)中的射線追蹤算法陸續(xù)推廣到了各向異性介質(zhì)中。如Faria等[5]將旅行時(shí)非線性插值算法推廣到各向異性介質(zhì)中，計(jì)算了初至qP波旅行時(shí)，但在計(jì)算每個(gè)節(jié)點(diǎn)旅行時(shí)的時(shí)候，在每個(gè)網(wǎng)格中都要分別計(jì)算射線方向與介質(zhì)對(duì)稱軸的夾角，因此當(dāng)模型區(qū)域較大時(shí)，計(jì)算效率往往比較低；鄧懷群等[6]對(duì)旅行時(shí)非線性插值算法進(jìn)行了改進(jìn)，可用于計(jì)算VTI介質(zhì)中直達(dá)波、反射波和透射波的旅行時(shí)，但在模型剖分過(guò)程中網(wǎng)格仍為單一的矩形網(wǎng)格，在計(jì)算過(guò)程中如果網(wǎng)格過(guò)大則降低了計(jì)算精度，若網(wǎng)格太小則降低了計(jì)算效率；Alkhalifah[7]利用有限差分求解程函方程實(shí)現(xiàn)了各向異性介質(zhì)中初至波旅行時(shí)的計(jì)算，但是難以實(shí)現(xiàn)，因此，在一般情況下，VTI介質(zhì)的地震波旅行時(shí)的計(jì)算依靠在弱各向異性假設(shè)條件下的擾動(dòng)理論對(duì)各向同性參考模型的旅行時(shí)校正的方法來(lái)實(shí)現(xiàn)；Zhou等[8]介紹了群速度在任意各向異性介質(zhì)中的計(jì)算方法，并與最短路徑射線追蹤方法結(jié)合，實(shí)現(xiàn)了二維和三維任意各向異性介質(zhì)中初至波和一次反射或反射轉(zhuǎn)換波的射線追蹤，該方法計(jì)算原理簡(jiǎn)單，易于實(shí)現(xiàn)，而且可以適應(yīng)較復(fù)雜的模型，但在實(shí)現(xiàn)過(guò)程中，為了提高計(jì)算精度，在矩形網(wǎng)格邊界上或長(zhǎng)方體各個(gè)面上增加了次生節(jié)點(diǎn)，犧牲了計(jì)算效率；趙愛(ài)華等[9]對(duì)旅行時(shí)最小樹(shù)算法和地震波群速度的射線角近似表示進(jìn)行了改進(jìn)，計(jì)算了初至qP波、qSV波和qSH波的旅行時(shí)，理論上具有較高的計(jì)算效率和精度，但在較復(fù)雜的地質(zhì)模型中，往往由于網(wǎng)格剖分的原因?qū)е掠?jì)算精度或效率達(dá)不到預(yù)期效果；馬德堂等[10]利用Thomsen給出的VTI介質(zhì)中相速度、群角、相角以及群速度之間的精確函數(shù)關(guān)系結(jié)合旅行時(shí)非線性插值法，實(shí)現(xiàn)了VTI介質(zhì)初至qP波的旅行時(shí)計(jì)算，由于在各向異性介質(zhì)中地震波群速度和群角被表示成了相角的復(fù)雜關(guān)系，為了精確計(jì)算給定群角的群速度，要進(jìn)行反復(fù)搜索計(jì)算，因此計(jì)算量較大、效率較低；李建國(guó)等[11]利用傳統(tǒng)的試射法實(shí)現(xiàn)了VTI介質(zhì)中反射qP波旅行時(shí)的計(jì)算，試射法在計(jì)算過(guò)程中利用了漸近線的原理，通過(guò)一次次的調(diào)整射線參數(shù)得到最終的計(jì)算結(jié)果，但在計(jì)算復(fù)雜介質(zhì)時(shí)存在陰影區(qū)；馬德堂等[12]利用由群角反插值相角來(lái)實(shí)線局部旅行時(shí)(網(wǎng)格中某個(gè)節(jié)點(diǎn)或插值點(diǎn)的旅行時(shí))的計(jì)算，并結(jié)合旅行時(shí)非線性插值法實(shí)現(xiàn)了TTI介質(zhì)中初至qP波的計(jì)算，雖然有原理簡(jiǎn)單、易于實(shí)現(xiàn)的優(yōu)點(diǎn)，但同樣存在計(jì)算效率與計(jì)算精度相互制約的缺點(diǎn)；白海軍等[13]將波前構(gòu)建法與用相速度和群速度重構(gòu)的射線追蹤算法相結(jié)合實(shí)現(xiàn)了TTI介質(zhì)的初至波射線追蹤，但算法實(shí)現(xiàn)較為復(fù)雜，且沒(méi)有考慮網(wǎng)格剖分帶來(lái)的誤差，因此仍有一些改進(jìn)空間；李曉玲等[14]利用在混合網(wǎng)格邊界加入次級(jí)節(jié)點(diǎn)的方法求取qP波、qSV波和qSH波的群速度，結(jié)合最短路徑算法實(shí)現(xiàn)了起伏層狀VTI介質(zhì)的多次波射線追蹤，但在求取群速度時(shí)，如果次級(jí)節(jié)點(diǎn)較多，計(jì)算效率會(huì)有所下降。此外，還有學(xué)者對(duì)基于網(wǎng)格單元擴(kuò)展的射線追蹤算法做出了有益的研究[15-21]。

上述大多數(shù)算法僅針對(duì)水平地表?xiàng)l件下各向異性介質(zhì)中的初至波、一次反射波或反射轉(zhuǎn)換波的旅行時(shí)進(jìn)行了計(jì)算，但對(duì)于復(fù)雜地質(zhì)模型會(huì)導(dǎo)致陰影區(qū)域且計(jì)算效率低。而對(duì)于起伏地表?xiàng)l件下基于網(wǎng)格單元擴(kuò)展的射線追蹤算法，多數(shù)仍以單一的網(wǎng)格形狀進(jìn)行模型剖分，若網(wǎng)格間距過(guò)大則計(jì)算精度降低，若網(wǎng)格間距過(guò)小則計(jì)算效率降低，混合網(wǎng)格是有效解決這一問(wèn)題的途徑，如李曉玲等[14]提出的分區(qū)多步不規(guī)則網(wǎng)格最短路徑(ISPM)算法，計(jì)算了各向異性介質(zhì)混合網(wǎng)格起伏地表?xiàng)l件下的多次波的旅行時(shí)。

LTI算法是Asakawa等[22]提出的基于線性假設(shè)的網(wǎng)格單元擴(kuò)展射線追蹤算法。由于該方法計(jì)算速度快、精度高、原理簡(jiǎn)單，是傳統(tǒng)的有限差分解程函方程方法的一種高級(jí)形式，且計(jì)算精度高于傳統(tǒng)的有限差分解程函方程法。LTI方法近年來(lái)得到了很大發(fā)展，不少學(xué)者對(duì)該方法進(jìn)行了大量研究和改進(jìn)。Li等[23]對(duì)LTI算法的計(jì)算公式進(jìn)行了改進(jìn)，使該方法對(duì)平面波假設(shè)的依賴性降低，提高了計(jì)算精度；趙改善等[24]結(jié)合界面二次源法將該方法推廣，可以用于追蹤反射波旅行時(shí)，彌補(bǔ)了僅能計(jì)算初至波旅行時(shí)的缺陷；Cardaerlli等[25]將該方法用于橢圓各向異性介質(zhì)中地震波旅行時(shí)的計(jì)算；聶建新等[26]將旅行時(shí)二次插值與線性插值方法聯(lián)合，降低了累積誤差；Kumar等[27]將該算法進(jìn)一步推廣到TTI介質(zhì)中；張賽民等[28]用拋物線插值取代線性插值，減小了因線性插值引起的誤差；張東等[29，30]為了求得網(wǎng)格加點(diǎn)的最小旅行時(shí)，在向前處理過(guò)程中采用了多方向循環(huán)的計(jì)算方法，彌補(bǔ)了計(jì)算過(guò)程中射線方向沒(méi)有考慮完全導(dǎo)致計(jì)算節(jié)點(diǎn)旅行時(shí)并不一定是最小值的缺陷，并將該方法擴(kuò)展到三維介質(zhì)。

綜上所述，近年來(lái)各向異性介質(zhì)起伏界面或地表地震波傳播問(wèn)題倍受關(guān)注，是現(xiàn)階段復(fù)雜環(huán)境地震勘探必須面對(duì)的重要問(wèn)題。研究混合網(wǎng)格剖分，是提高基于網(wǎng)格單元擴(kuò)展的射線追蹤算法在起伏界面或地表?xiàng)l件下各向異性介質(zhì)各種地震波旅行時(shí)計(jì)算精度和效率的有效途徑。本文利用基于網(wǎng)格單元擴(kuò)展的LTI射線追蹤算法和Sena[4]推導(dǎo)的適用于弱各向異性介質(zhì)的群速度計(jì)算公式，在向前和向后處理過(guò)程中結(jié)合分區(qū)多步計(jì)算技術(shù)，實(shí)現(xiàn)了一種可以計(jì)算起伏地表(界面)地震初至波以及多種類型后續(xù)波的VTI介質(zhì)射線追蹤算法。在利用LTI算法計(jì)算時(shí)，充分考慮了射線的傳播方向，并在計(jì)算過(guò)程中采用全方位循環(huán)的方法計(jì)算各個(gè)網(wǎng)格節(jié)點(diǎn)的旅行時(shí)，以保證每一個(gè)節(jié)點(diǎn)的旅行時(shí)都是最小值。

2 VTI介質(zhì)起伏地表混合網(wǎng)格LTI射線追蹤方法

VTI介質(zhì)混合網(wǎng)格LTI射線追蹤方法的實(shí)現(xiàn)過(guò)程可以分為四個(gè)步驟：模型剖分、計(jì)算局部旅行時(shí)、計(jì)算全局旅行時(shí)和射線路徑、計(jì)算后續(xù)波旅行時(shí)。

2.1 模型剖分

為了更好地逼近實(shí)際地表或速度界面，同時(shí)又能兼顧計(jì)算效率，本文采用矩形網(wǎng)格與不規(guī)則四邊形網(wǎng)格相結(jié)合的方法對(duì)模型進(jìn)行剖分，如圖1所示。首先對(duì)整個(gè)模型區(qū)域用矩形網(wǎng)格剖分，分別計(jì)算地表和速度界面與網(wǎng)格縱邊的交點(diǎn)，得到離散地表點(diǎn)和速度界面點(diǎn)；然后將相鄰的界面離散點(diǎn)或網(wǎng)格節(jié)點(diǎn)依次連接，就構(gòu)成一條由多條折線段組成的近似于實(shí)際地表或速度界面的折線，把這些折線段與原矩形網(wǎng)格邊界和節(jié)點(diǎn)組合，構(gòu)成了不規(guī)則四邊形網(wǎng)格；最后對(duì)除地表以上的各個(gè)網(wǎng)格節(jié)點(diǎn)賦予相應(yīng)速度值，就完成了對(duì)整個(gè)模型的剖分。

圖1 混合網(wǎng)格模型剖分示意圖

上述混合網(wǎng)格在對(duì)模型剖分時(shí)會(huì)遇到三種情況，如圖1中的①、②、③所示的陰影部分。為了便于觀察，將陰影部分按順序放大，如圖2所示。

(1)當(dāng)?shù)乇砘蛩俣冉缑?個(gè)相鄰的離散點(diǎn)都在原矩形網(wǎng)格的同一層內(nèi)，直接連接三個(gè)相鄰的離散點(diǎn)即可構(gòu)成不規(guī)則四邊形網(wǎng)格(圖2中①)。

(2)當(dāng)左端或右端的離散點(diǎn)和另兩個(gè)點(diǎn)不在原矩形網(wǎng)格的同一層內(nèi)，首先連接在原矩形網(wǎng)格同一層的兩個(gè)離散點(diǎn)，然后連接與中間離散點(diǎn)相鄰的可與另一個(gè)離散點(diǎn)(或原矩形網(wǎng)格節(jié)點(diǎn))構(gòu)成四邊形的網(wǎng)格節(jié)點(diǎn)；依次在其他層中構(gòu)造四邊形網(wǎng)格，這樣就構(gòu)成了不規(guī)則四邊形網(wǎng)格(圖2中②)。

(3)當(dāng)3個(gè)相鄰離散點(diǎn)中兩端的離散點(diǎn)和中間的離散點(diǎn)不在原矩形網(wǎng)格的同一層內(nèi)，分別連接與中間離散點(diǎn)相鄰的同一層網(wǎng)格內(nèi)可與其余兩個(gè)離散點(diǎn)(或原矩形網(wǎng)格節(jié)點(diǎn))構(gòu)成四邊形的網(wǎng)格節(jié)點(diǎn)，依次在其他層中構(gòu)造四邊形網(wǎng)格單元，這樣就構(gòu)成了不規(guī)則四邊形網(wǎng)格(圖2中③)。

圖2中各個(gè)顏色不同的區(qū)域組成了利用上述方法建立的混合網(wǎng)格。最后將Thomsen參數(shù)δ、ε和γ以及縱橫波的背景速度和密度賦予各個(gè)網(wǎng)格節(jié)點(diǎn)上，這樣就完成了VTI介質(zhì)混合網(wǎng)格模型的建立。

由圖1和圖2可以看出，混合網(wǎng)格剖分逼近的地表或界面(紅色實(shí)線)比單一矩形網(wǎng)格剖分逼近(藍(lán)色實(shí)線)的擬合程度更高，剖分方法也更合理。

圖2 混合網(wǎng)格剖分局部放大顯示

2.2 局部旅行時(shí)計(jì)算公式

地震波運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)特征在各向異性介質(zhì)中有著較大的變化，其中速度的變化最為突出。群速度和相速度不再相同，群角和相角發(fā)生了分離，如圖3 所示。其中，相速度指波前的傳播速度，傳播方向?yàn)橄辔蛔兓羁斓姆较?，且始終垂直于波前面。而相速度矢量與介質(zhì)對(duì)稱軸的夾角稱為相角，用θ表示，VTI介質(zhì)的對(duì)稱軸與z軸平行。群速度指地震波能量的傳播方向，與介質(zhì)對(duì)稱軸的夾角稱為群角，用φ表示。在各向異性介質(zhì)中準(zhǔn)確求解群速度與相速度非常困難，經(jīng)過(guò)眾多學(xué)者的深入研究，得到了可用于數(shù)值計(jì)算的群速度的近似計(jì)算公式[3，4，31-33]。Sena[4]在弱各向異性線性近似下，分別推導(dǎo)出了qP波、qSV波和qSH波的群速度vP(φ)、vSV(φ)和vSH(φ)隨群角φ變化的表達(dá)式

圖3 VTI介質(zhì)混合網(wǎng)格中由規(guī)則邊界計(jì)算局部旅行時(shí)示意圖(a)射線從網(wǎng)格橫邊穿過(guò)； (b)射線從網(wǎng)格縱邊穿過(guò)

(1)

式中：δ、ε和γ為T(mén)homsen參數(shù)；vP0、vS0分別為各向異性介質(zhì)對(duì)稱軸方向上的相速度；ρ為介質(zhì)的密度。

本文提出的混合網(wǎng)格由矩形和不規(guī)則四邊形組成，因此在計(jì)算局部旅行時(shí)要對(duì)網(wǎng)格邊界(主要針對(duì)旅行時(shí)已知的邊界)加以區(qū)分。如果網(wǎng)格的AB邊界與原矩形網(wǎng)格邊界重合，如圖3所示，則稱此邊界為規(guī)則邊界； 反之，如圖4所示，則稱此邊界為不規(guī)則邊界。混合網(wǎng)格LTI算法仍然滿足矩形網(wǎng)格LTI算法的前提假設(shè)。

如果AB邊為規(guī)則的水平邊界，如圖3a所示，射線穿過(guò)AB邊界的D點(diǎn)到達(dá)C點(diǎn)，則C點(diǎn)的旅行時(shí)tC可表示為

(2)

式中：s代表C點(diǎn)附近的平均慢度，在VTI介質(zhì)中不再是一個(gè)常數(shù)而是群角φ的函數(shù)； Δt=tB-tA為B點(diǎn)與A點(diǎn)地旅行時(shí)差；l、d1、d2、d3的含義如圖3所示。由幾何關(guān)系可將式(1)表示成關(guān)于l的函數(shù)

圖4 由不規(guī)則邊界計(jì)算局部旅行時(shí)示意圖

(3)

式中：sP0和sS0分別為P波和S波的垂直相慢度；l∈(0，d2)。將式(3)代入式(2)，可得

(4)

式中tPC、tSVC和tSHC分別為穿過(guò)D點(diǎn)的射線以sP、sSV和sSH為慢度到達(dá)C點(diǎn)的局部旅行時(shí)。

根據(jù)費(fèi)馬原理，tPC、tSVC和tSHC在D點(diǎn)應(yīng)滿足?tPC/?l=0、?tSVC/?l=0和?tSHC/?l=0，可得

(5)

為了求解式(5)，構(gòu)造三個(gè)關(guān)于l的非線性函數(shù)fHP(l)、fHSV(l)和fHSH(l)

(6)

式中Δt與l同號(hào)，且l∈(0，d2)，利用二分法求出以上三個(gè)函數(shù)等于零時(shí)的正實(shí)數(shù)根，即可求得局部旅行時(shí)tPC、tSVC和tSHC。解得l后即可求出D的坐標(biāo)

(7)

同理,當(dāng)AB邊為規(guī)則的垂直邊界時(shí)(圖3b)，通過(guò)與上述過(guò)程類似的推導(dǎo)過(guò)程，構(gòu)造三個(gè)關(guān)于l的非線性函數(shù)fVP(l)、fVSV(l)和fVSH(l)

(8)

同樣可用二分法求出以上三個(gè)函數(shù)等于零時(shí)的正實(shí)數(shù)根，即可求得局部旅行時(shí)tPC、tSVC和tSHC。D點(diǎn)的坐標(biāo)可表示為

(9)

當(dāng)AB邊為不規(guī)則邊界時(shí)(圖4)，將AB邊界離散成若干個(gè)離散點(diǎn)，將每一個(gè)離散點(diǎn)都考慮為D點(diǎn)，且D點(diǎn)的局部旅行時(shí)通過(guò)A、B點(diǎn)的旅行時(shí)線性插值得到，則C點(diǎn)附近的慢度sP、sSV和sSH可表示為

(10)

C點(diǎn)的局部旅行時(shí)tPC、tSVC和tSHC為

(11a)

(11b)

以式為在混合網(wǎng)格VTI介質(zhì)中計(jì)算局部旅行時(shí)tPC、tSVC和tSHC和插值點(diǎn)坐標(biāo)(xD，zD)的計(jì)算公式。

2.3 全局旅行時(shí)和射線路徑的計(jì)算

LTI算法的實(shí)現(xiàn)分為向前和向后處理兩個(gè)步驟，向前處理過(guò)程得到全局旅行時(shí)，向后處理過(guò)程得到射線路徑。

向前處理過(guò)程包括以下步驟。

步驟1：計(jì)算炮點(diǎn)所在網(wǎng)格的各個(gè)節(jié)點(diǎn)的旅行時(shí)(圖5a)。

步驟2：計(jì)算炮點(diǎn)網(wǎng)格所在列的每個(gè)網(wǎng)格上的各個(gè)節(jié)點(diǎn)的局部旅行時(shí)(圖5b)。

步驟3：從左向右逐步計(jì)算炮點(diǎn)網(wǎng)格右側(cè)的每一列上所有網(wǎng)格節(jié)點(diǎn)的旅行時(shí)，具體如下：

(1)從上向下由該列左邊界上的旅行時(shí)已知的節(jié)點(diǎn)計(jì)算右邊界上的旅行時(shí)未知節(jié)點(diǎn)的旅行時(shí)，此時(shí)假設(shè)到達(dá)該列右邊界節(jié)點(diǎn)的射線僅來(lái)自左邊界(圖5c)，黑色虛線表示射線；

(2)從下向上由網(wǎng)格下邊界上的節(jié)點(diǎn)旅行時(shí)計(jì)算上邊界上的節(jié)點(diǎn)旅行時(shí)，若計(jì)算得到的旅行時(shí)小于之前計(jì)算所得的旅行時(shí)，則更新該節(jié)點(diǎn)的旅行時(shí)，使其取得最小值(圖5d)；

(3)從上向下由網(wǎng)格上邊界上的節(jié)點(diǎn)旅行時(shí)計(jì)算網(wǎng)格下邊界上的節(jié)點(diǎn)旅行時(shí)，若計(jì)算得到的旅行時(shí)小于之前計(jì)算所得，則更新該節(jié)點(diǎn)的旅行時(shí)，使其取得最小值(圖5e)；

(4)循環(huán)執(zhí)行(1)、(2)和(3)，直到整個(gè)模型最右邊一列，就可以得到炮點(diǎn)網(wǎng)格所在列右側(cè)每個(gè)節(jié)點(diǎn)的旅行時(shí)。

步驟4：炮點(diǎn)所在列左側(cè)網(wǎng)格節(jié)點(diǎn)的旅行時(shí)的計(jì)算方法與步驟3類似，計(jì)算完成后就得到了整個(gè)模型區(qū)域中每個(gè)網(wǎng)格節(jié)點(diǎn)的旅行時(shí)(圖5f)。

雖然通過(guò)上述步驟計(jì)算得到了整個(gè)模型所有網(wǎng)格節(jié)點(diǎn)的旅行時(shí)，但是該計(jì)算過(guò)程沒(méi)有考慮到在復(fù)雜介質(zhì)中會(huì)遇到射線逆向傳播的現(xiàn)象，仍需按行掃描重新計(jì)算每個(gè)節(jié)點(diǎn)的旅行時(shí)。

步驟5：計(jì)算炮點(diǎn)網(wǎng)格所在行的每個(gè)網(wǎng)格上的各個(gè)節(jié)點(diǎn)的局部旅行時(shí)(圖6a)，計(jì)算結(jié)果與按列掃描的計(jì)算結(jié)果相比較，取最小值。

步驟6：炮點(diǎn)網(wǎng)格所在行的所有網(wǎng)格節(jié)點(diǎn)的旅行時(shí)已經(jīng)得到更新，從上向下逐步計(jì)算炮點(diǎn)網(wǎng)格下側(cè)每一行上所有網(wǎng)格節(jié)點(diǎn)的最小旅行時(shí)，具體如下：

(1)從右向左由該行上邊界上的節(jié)點(diǎn)的旅行時(shí)計(jì)算下邊界上的節(jié)點(diǎn)的旅行時(shí)，此時(shí)假設(shè)到達(dá)該行下邊界節(jié)點(diǎn)的射線僅來(lái)自上邊界(圖6b)；

(2)從左向右由網(wǎng)格左邊界上的節(jié)點(diǎn)旅行時(shí)計(jì)算右邊界上的節(jié)點(diǎn)旅行時(shí)，若計(jì)算得到的旅行時(shí)小于之前的旅行時(shí)，則更新該節(jié)點(diǎn)的旅行時(shí)，使其取得最小值(圖6c)；

圖5 向前處理中旅行時(shí)按列掃描計(jì)算網(wǎng)格節(jié)點(diǎn)旅行時(shí)(a)計(jì)算炮點(diǎn)網(wǎng)格節(jié)點(diǎn)旅行時(shí)； (b)計(jì)算炮點(diǎn)網(wǎng)格所在列節(jié)點(diǎn)的旅行時(shí)； (c)僅考慮射線來(lái)自網(wǎng)格左邊界點(diǎn)的情況；d)僅考慮射線來(lái)自網(wǎng)格下邊界的情況； (e)僅考慮射線來(lái)自網(wǎng)格上邊界的情況； (f)網(wǎng)格所有節(jié)點(diǎn)計(jì)算完成

圖6 向前處理過(guò)程按行掃描計(jì)算網(wǎng)格節(jié)點(diǎn)最小旅行時(shí)(a)計(jì)算炮點(diǎn)網(wǎng)格所在行節(jié)點(diǎn)的最小旅行時(shí)；(b)僅考慮射線來(lái)自網(wǎng)格上邊界；(c)僅考慮射線來(lái)自網(wǎng)格左邊界；(d)僅考慮射線來(lái)自網(wǎng)格右邊界

(3)從右向左由網(wǎng)格右邊界上的節(jié)點(diǎn)旅行時(shí)計(jì)算網(wǎng)格左邊界上的節(jié)點(diǎn)旅行時(shí)，若計(jì)算得到的旅行時(shí)小于之前的旅行時(shí)，則更新該節(jié)點(diǎn)的旅行時(shí)，使其取得最小值(圖6d)；

(4)循環(huán)執(zhí)行(1)、(2)和(3),直到整個(gè)模型最下邊一行，得到炮點(diǎn)網(wǎng)格所在行下側(cè)每個(gè)節(jié)點(diǎn)的旅行時(shí)。

步驟7：計(jì)算炮點(diǎn)所在行上側(cè)網(wǎng)格節(jié)點(diǎn)的旅行時(shí)，計(jì)算方法與步驟6類似，計(jì)算完成后得到了整個(gè)模型區(qū)域中每個(gè)網(wǎng)格節(jié)點(diǎn)的最小旅行時(shí)。

向后處理過(guò)程包括以下步驟。

步驟1：跟向前處理過(guò)程類似，這里將接收點(diǎn)作為2.2節(jié)中的C點(diǎn)，將接收點(diǎn)附近的每一條網(wǎng)格邊界作為2.2節(jié)中的AB邊，即可計(jì)算出每一條網(wǎng)格邊界上的插值點(diǎn)，即D點(diǎn)。利用t=tmin+ds計(jì)算出接收點(diǎn)到其附近所有網(wǎng)格邊界插值點(diǎn)的旅行時(shí)(圖7a)，其中tmin為插值點(diǎn)的旅行時(shí)，d為接收點(diǎn)到插值點(diǎn)的距離，s為插值點(diǎn)的慢度；

步驟2：根據(jù)步驟1得到了接收點(diǎn)附近所有網(wǎng)格邊界上的插值點(diǎn)坐標(biāo)，那么到達(dá)接收點(diǎn)的射線一定經(jīng)過(guò)了這些插值點(diǎn)，在這些點(diǎn)中滿足t=tmin+ds取最小值的點(diǎn)就是要求的插值點(diǎn)；

步驟3：將步驟2得到的插值點(diǎn)坐標(biāo)作為新的接收點(diǎn)，重復(fù)步驟1、步驟2直到插值點(diǎn)出現(xiàn)在震源所在網(wǎng)格的邊界或節(jié)點(diǎn)上為止(圖7b和圖7c)；

步驟4：依次將震源與插值點(diǎn)相連得到整條初至波射線路徑(圖7d)，將每?jī)蓚€(gè)插值點(diǎn)(包括接收點(diǎn))間的旅行時(shí)相加得到初至波的旅行時(shí)。

圖7 向后處理過(guò)程計(jì)算初至波射線路徑(a)計(jì)算接收點(diǎn)到其所在網(wǎng)格節(jié)點(diǎn)的最小旅行時(shí)，求得射線與網(wǎng)格的交點(diǎn)； (b)和(c)將確定的交點(diǎn)作為新的接收點(diǎn)，求下一個(gè)交點(diǎn)直到震源所在網(wǎng)格； (d)連接震源與計(jì)算得到的各個(gè)交點(diǎn)，得到射線路徑

2.4 后續(xù)波的旅行時(shí)計(jì)算

將VTI介質(zhì)起伏地表初至波混合網(wǎng)格LTI算法與分區(qū)多步計(jì)算技術(shù)[12]結(jié)合，可以實(shí)現(xiàn)任意多次反射(或反射轉(zhuǎn)換)或透射(或透射轉(zhuǎn)換)波的追蹤。在計(jì)算時(shí)需要引入分區(qū)多步的思想(模型分區(qū)和多步計(jì)算)，即將整個(gè)模型區(qū)域按所需追蹤的波的類型和速度特征分成相應(yīng)幾個(gè)獨(dú)立的計(jì)算區(qū)域，在每個(gè)區(qū)域內(nèi)獨(dú)立進(jìn)行計(jì)算。分區(qū)多步計(jì)算首先要建立模型并網(wǎng)格化，速度界面也要離散成不連續(xù)的點(diǎn)。

2.4.1 模型分區(qū)

圖8所示為多次反射、透射轉(zhuǎn)換波在二維層狀介質(zhì)模型中的分區(qū)情況，圖中顯示了一條在界面折返三次、轉(zhuǎn)換三次的多次反射轉(zhuǎn)換波。根據(jù)分區(qū)多步的思想，該多次反射轉(zhuǎn)換波旅行時(shí)的計(jì)算過(guò)程可以在四個(gè)獨(dú)立的計(jì)算區(qū)域(每個(gè)標(biāo)號(hào)代表一個(gè)分區(qū))內(nèi)進(jìn)行。具體實(shí)現(xiàn)過(guò)程為：該多次反射轉(zhuǎn)換波由四段組成，第一段是來(lái)自地表到第二個(gè)反射界面的下行qP波，將第一、第二層作為第一個(gè)分區(qū)計(jì)算第一段下行qP波的波前，采用第一、二層的qP波速度；第二段是自第二個(gè)反射界面上的反射點(diǎn)到第一個(gè)反射界面的上行qSV波，將第二層作為第二個(gè)分區(qū)來(lái)計(jì)算第二段上行qSV波的波前，采用第二層的qSV波速度；第三段是自第一個(gè)反射界面反射點(diǎn)到第二個(gè)反射界面的下行qP波，將第二層作為第三個(gè)分區(qū)計(jì)算第三段下行qP波的波前，采用第二層的qP波速度；第四段是自第二個(gè)反射界面反射點(diǎn)到地表的上行qSH波，將第一、第二層作為第四個(gè)分區(qū)計(jì)算第四段上行qSH波的波前，采用第一、二層的qSH波速度。

圖8 二維層狀介質(zhì)模型多次反射、透射轉(zhuǎn)換波分區(qū)示意圖

2.4.2 多步計(jì)算

圖9所示為反射(或反射轉(zhuǎn)換)或透射(或透射轉(zhuǎn)換)波的波前擴(kuò)展示意圖。具體實(shí)現(xiàn)過(guò)程為：

(1)計(jì)算由震源到第一個(gè)分區(qū)下界面的下行波的波前，并找出界面離散點(diǎn)上旅行時(shí)最小的點(diǎn)(圖9a)。

(2)將(1)中計(jì)算出的旅行時(shí)最小的界面離散點(diǎn)作為新的震源。如果需要計(jì)算反射波的波前，則由界面上新的震源向上計(jì)算上行波在該分區(qū)內(nèi)的波前(圖9b)；如果需要追蹤轉(zhuǎn)換波，則采用轉(zhuǎn)換波的速度即可計(jì)算其相應(yīng)的波前。如果需要計(jì)算透射波的波前，由界面上新的震源向下計(jì)算下行波在該分區(qū)內(nèi)的波前(圖9c)；如果需要追蹤轉(zhuǎn)換波，則采用轉(zhuǎn)換波的速度即可計(jì)算其相應(yīng)的波前。

(3)結(jié)合(1)、(2)計(jì)算波前，可以追蹤任意多次反射(或反射轉(zhuǎn)換)或透射(或透射轉(zhuǎn)換)波。

(4)由以上三步計(jì)算波前，再結(jié)合LTI的向后處理,追蹤任意多次反射轉(zhuǎn)換波的射線路徑。

圖9 波前分區(qū)擴(kuò)展示意圖

(a)波前由震源擴(kuò)展到速度界面離散點(diǎn)；(b)若反射，由界面離散點(diǎn)旅行時(shí)計(jì)算上行波波前；(c)若透射，由界面離散點(diǎn)旅行時(shí)計(jì)算下行波波前

3 數(shù)值算例

首先用具有解析解的起伏地表均勻VTI介質(zhì)模型檢驗(yàn)本文算法的相對(duì)誤差；然后將本文算法分別與波動(dòng)方程有限差分算法、分區(qū)多步快速步進(jìn)(FMM)算法和分區(qū)多步ISPM算法進(jìn)行精度和效率對(duì)比；最后對(duì)帶有斷層的起伏層狀模型進(jìn)行正演模擬，驗(yàn)證算法對(duì)復(fù)雜模型的適應(yīng)能力。

3.1 均勻VTI介質(zhì)模型

VTI均勻介質(zhì)模參數(shù)為：模型尺寸為400m×400m；網(wǎng)格尺寸為1m×1m；vP0=2700m/s；vS0=1500m/s；ε=0.15；δ=0.08；γ=0.338；ρ=2.10g/cm3。圖10為在均勻VTI介質(zhì)模型中三種地震波(qP、qSV 和qSH)的波前擴(kuò)展過(guò)程以及三種地震波旅行時(shí)數(shù)值解與解析解的相對(duì)誤差分布。

從圖10a～圖10c可以看出：三種地震波的各向異性效應(yīng)都比較明顯，且在每種地震波的傳播過(guò)程中，紅線(各向異性介質(zhì))與黑線(各向同性介質(zhì))總是在速度最慢的方向相切。根據(jù)式(1)可知，在各向異性介質(zhì)中群速度最慢的方向的速度與在各向同性介質(zhì)中是相等的；qSV與qSH的旅行時(shí)擴(kuò)展特征不同，前者在對(duì)角線方向傳播最快、在垂直和水平方向上傳播最慢，而后者在水平方向傳播最快、在垂直方向上傳播最慢。從圖10d～圖10f可以看出,由本文方法計(jì)算的三種波的旅行時(shí)與解析解的相對(duì)誤差均較小，說(shuō)明算法較精確。

圖10 均勻VTI介質(zhì)模型qP波、qSV波和qSH波旅行時(shí)等值線(單位ms)及相對(duì)誤差(a)、(d)qP波； (b)、(e)qSV波； (c)、(f)qSH波黑線為各向同性介質(zhì)的旅行時(shí)等值線，紅線為VTI介質(zhì)的旅行時(shí)等值線(單位：ms)

3.2 兩層水平層狀模型

設(shè)計(jì)兩層水平層狀VTI介質(zhì)模型，尺寸為1000m×1000m，各層參數(shù)如表1所示。源位于模型表面(500m，0)處。圖11a和圖11b分別為由本文算法和有限差分算法計(jì)算的模型表面101道接收的地震記錄，道間距為10m。圖11c為圖11a和圖11b的疊合顯示，可以看出兩種算法計(jì)算的初至波、反射波和轉(zhuǎn)換波的記錄都比較吻合，說(shuō)明本文算法的計(jì)算結(jié)果比較準(zhǔn)確。

表1 兩層水平層狀VTI介質(zhì)模型參數(shù)

目前基于網(wǎng)格單元波前擴(kuò)展的波前追蹤方法中，較為成熟的主要有分區(qū)多步FMM算法和分區(qū)多步ISPM算法，為了驗(yàn)證分區(qū)多步LTI算法的計(jì)算精度和計(jì)算效率，對(duì)文獻(xiàn)[34]的兩層水平層狀各向同性速度模型(表2)分別用以上三種方法、在4種網(wǎng)格間距下計(jì)算反射波旅行時(shí)。模型尺寸為100km×40km，地表水平，在深度30km處有一水平反射界面。在模型參數(shù)化時(shí)，分區(qū)多步FMM算法中正方形網(wǎng)格邊長(zhǎng)分別為1000m、500m、250m和125m；分區(qū)多步ISPM算法中網(wǎng)格尺寸為4km×4km，并在網(wǎng)格單元邊界上相應(yīng)加入3、7、15和31個(gè)次級(jí)節(jié)點(diǎn)；分區(qū)多步LTI算法中網(wǎng)格大小與分區(qū)多步FMM算法中設(shè)置相同。炮點(diǎn)在模型的左上角，100個(gè)檢波器等間距(1km)排列在地面，最小炮檢距為1.0km。圖12和圖13分別給出了隨炮檢距變化的分區(qū)多步FMM算法(圖12a)、分區(qū)多步ISPM算法(圖12b)和分區(qū)多步LTI算法(圖13)相對(duì)于解析解在4種不同網(wǎng)格距下反射P波旅行時(shí)的相對(duì)誤差。表3給出了三種方法的CPU耗時(shí)。由圖12和表3可以看出，除了網(wǎng)格尺寸為1000m×1000m的情形外，分區(qū)多步LTI算法的精度均優(yōu)于分區(qū)多步FMM和分區(qū)多步ISPM算法，在計(jì)算效率方面，分區(qū)多步LTI算法明顯優(yōu)于分區(qū)多步FMM算法，但低于分區(qū)多步ISPM算法。

圖11 本文算法計(jì)算的地震記錄(a)、波動(dòng)方程有限差分法計(jì)算的地震記錄(b)和二者疊合后的記錄(c)對(duì)比①為初至qP波，②為初至qSV波，③為反射qP波，④、⑤為反射轉(zhuǎn)換波，⑥為反射qSV波

圖12 不同網(wǎng)格間距兩層水平層狀各向同性介質(zhì)模型的分區(qū)多步FMM算法(a)、ISPM算法(b)反射波旅行時(shí)的相對(duì)誤差(引自文獻(xiàn)[34])

圖13 不同網(wǎng)格間距兩層水平層狀各向同性介質(zhì)模型的分區(qū)多步LTI算法數(shù)值解的相對(duì)誤差

表2 兩層水平層狀各向同性介質(zhì)模型參數(shù)

表3 不同網(wǎng)格間距分區(qū)多步FMM、ISPM

注：表中FMM和ISPM算法數(shù)據(jù)引自文獻(xiàn)[34]

3.3 二維起伏層狀VTI介質(zhì)

圖14a為二維起伏層狀VTI介質(zhì)模型的參數(shù)和分區(qū)情況。在(200m，20m)處激發(fā)，101道檢波器橫向均勻布設(shè)在起伏地表上。圖14b～圖14e為在各個(gè)分區(qū)內(nèi)的各向異性介質(zhì)和各向同性介質(zhì)的多次反射、轉(zhuǎn)換波波前擴(kuò)展過(guò)程，可以看出在速度均勻的同一層內(nèi)波前呈大致平行的弧線，當(dāng)波前穿過(guò)速度界面時(shí)，就會(huì)隨速度發(fā)生變化；同時(shí)可以看出VTI介質(zhì)與各向同性介質(zhì)中的波前在傳播過(guò)程中總是在速度最小的方向相切，說(shuō)明在各向異性介質(zhì)中波前擴(kuò)展方式正確。圖15為圖14a模型模擬的單炮合成記錄，包括VTI介質(zhì)qP波初至、界面Ⅰ反射qP波、界面Ⅱ反射qSH波、qSV波、VTI介質(zhì)多次反射、轉(zhuǎn)換波，各向同性介質(zhì)P波初至、界面Ⅰ反射P波、多次反射轉(zhuǎn)換波。由圖15可以看出，地震波在各向異性介質(zhì)和各向異性介質(zhì)中傳播的差異，且各向異性程度越強(qiáng)差異越大；在VTI介質(zhì)中傳播的qSV波和qSH波也有較大差異，僅在震源正下方的位置會(huì)重疊，這是因?yàn)樵赩TI介質(zhì)的對(duì)稱軸方向上qSV和qSH波的速度相同，偏離這個(gè)方向后qSH波和qSV波的各向異性效應(yīng)不同。

圖14 層間多次反射、透射波在各個(gè)分區(qū)中的傳播過(guò)程、射線路徑和單炮合成記錄(波前間隔為10ms)

(a)模型參數(shù)和分區(qū)編號(hào)； (b)下行P波波前在第一個(gè)分區(qū)內(nèi)傳播； (c)上行SV波波前在第二個(gè)分區(qū)內(nèi)傳播； (d)下行P波波前在第三個(gè)分區(qū)內(nèi)傳播； (e)上行SH波波前在第四個(gè)分區(qū)內(nèi)傳播； (f)層間多次波射、轉(zhuǎn)換波線路徑。紅色等時(shí)線和射線路徑為在VTI介質(zhì)中傳播，黑色等時(shí)線和射線路徑為在各向同性介質(zhì)中傳播

圖15 單炮合成記錄

4 結(jié)論

本文采用矩形網(wǎng)格和不規(guī)則四邊形網(wǎng)格組成的混合網(wǎng)格對(duì)模型進(jìn)行剖分，可以提高網(wǎng)格對(duì)起伏地表或界面的擬合程度，進(jìn)而提高計(jì)算精度。將混合網(wǎng)格中群速度與群角的關(guān)系式變換成群速度與插值點(diǎn)坐標(biāo)的關(guān)系，分別建立了針對(duì)qP、qSV和qSH波的非線性方程，并利用二分法求解，再結(jié)合分區(qū)多步計(jì)算技術(shù)，可計(jì)算VTI介質(zhì)中的初至波和多種后續(xù)波的射線路徑和旅行時(shí)。通過(guò)模型試算，與有限差分方法、分區(qū)多步FMM算法和分區(qū)多步ISPM算法進(jìn)行對(duì)比，驗(yàn)證了算法效率和精度。二維起伏層狀VTI介質(zhì)模型結(jié)果表明，本文方法適用于復(fù)雜構(gòu)造條件下VTI介質(zhì)中地震波的旅行時(shí)和射線路徑計(jì)算。