王 琦 李慶春 王芷琪

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

1 引言

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

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

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

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

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

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

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

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

2.1 模型剖分

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

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

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

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

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

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

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

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

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

2.2 局部旅行時計算公式

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

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

(1)

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

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

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

(2)

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

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

(3)

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

(4)

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

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

(5)

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

(6)

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

(7)

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

(8)

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

(9)

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

(10)

C點的局部旅行時tPC、tSVC和tSHC為

(11a)

(11b)

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

2.3 全局旅行時和射線路徑的計算

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

向前處理過程包括以下步驟。

步驟1：計算炮點所在網(wǎng)格的各個節(jié)點的旅行時(圖5a)。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

向后處理過程包括以下步驟。

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

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

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

步驟4：依次將震源與插值點相連得到整條初至波射線路徑(圖7d)，將每兩個插值點(包括接收點)間的旅行時相加得到初至波的旅行時。

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

2.4 后續(xù)波的旅行時計算

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

2.4.1 模型分區(qū)

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

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

2.4.2 多步計算

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

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

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

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

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

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

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

3 數(shù)值算例

首先用具有解析解的起伏地表均勻VTI介質(zhì)模型檢驗本文算法的相對誤差；然后將本文算法分別與波動方程有限差分算法、分區(qū)多步快速步進(jìn)(FMM)算法和分區(qū)多步ISPM算法進(jìn)行精度和效率對比；最后對帶有斷層的起伏層狀模型進(jìn)行正演模擬，驗證算法對復(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ò)展過程以及三種地震波旅行時數(shù)值解與解析解的相對誤差分布。

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

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

3.2 兩層水平層狀模型

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

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

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

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

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

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

表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為在各個分區(qū)內(nèi)的各向異性介質(zhì)和各向同性介質(zhì)的多次反射、轉(zhuǎn)換波波前擴(kuò)展過程，可以看出在速度均勻的同一層內(nèi)波前呈大致平行的弧線，當(dāng)波前穿過速度界面時，就會隨速度發(fā)生變化；同時可以看出VTI介質(zhì)與各向同性介質(zhì)中的波前在傳播過程中總是在速度最小的方向相切，說明在各向異性介質(zhì)中波前擴(kuò)展方式正確。圖15為圖14a模型模擬的單炮合成記錄，包括VTI介質(zhì)qP波初至、界面Ⅰ反射qP波、界面Ⅱ反射qSH波、qSV波、VTI介質(zhì)多次反射、轉(zhuǎn)換波，各向同性介質(zhì)P波初至、界面Ⅰ反射P波、多次反射轉(zhuǎn)換波。由圖15可以看出，地震波在各向異性介質(zhì)和各向異性介質(zhì)中傳播的差異，且各向異性程度越強差異越大；在VTI介質(zhì)中傳播的qSV波和qSH波也有較大差異，僅在震源正下方的位置會重疊，這是因為在VTI介質(zhì)的對稱軸方向上qSV和qSH波的速度相同，偏離這個方向后qSH波和qSV波的各向異性效應(yīng)不同。

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

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

圖15 單炮合成記錄

4 結(jié)論

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