VTI介質起伏界面混合網格旅行時線性插值計算方法

王 琦 李慶春 王芷琪

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

1 引言

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

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

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

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

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

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

2 VTI介質起伏地表混合網格LTI射線追蹤方法

VTI介質混合網格LTI射線追蹤方法的實現過程可以分為四個步驟：模型剖分、計算局部旅行時、計算全局旅行時和射線路徑、計算后續波旅行時。

2.1 模型剖分

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

圖1 混合網格模型剖分示意圖

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

(1)當地表或速度界面3個相鄰的離散點都在原矩形網格的同一層內，直接連接三個相鄰的離散點即可構成不規則四邊形網格(圖2中①)。

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

(3)當3個相鄰離散點中兩端的離散點和中間的離散點不在原矩形網格的同一層內，分別連接與中間離散點相鄰的同一層網格內可與其余兩個離散點(或原矩形網格節點)構成四邊形的網格節點，依次在其他層中構造四邊形網格單元，這樣就構成了不規則四邊形網格(圖2中③)。

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

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

圖2 混合網格剖分局部放大顯示

2.2 局部旅行時計算公式

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

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

(1)

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

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

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

(2)

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

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

(3)

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

(4)

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

根據費馬原理，tPC、tSVC和tSHC在D點應滿足?tPC/?l=0、?tSVC/?l=0和?tSHC/?l=0，可得

(5)

為了求解式(5)，構造三個關于l的非線性函數fHP(l)、fHSV(l)和fHSH(l)

(6)

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

(7)

同理,當AB邊為規則的垂直邊界時(圖3b)，通過與上述過程類似的推導過程，構造三個關于l的非線性函數fVP(l)、fVSV(l)和fVSH(l)

(8)

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

(9)

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

(10)

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

(11a)

(11b)

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

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

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

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

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

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

步驟3：從左向右逐步計算炮點網格右側的每一列上所有網格節點的旅行時，具體如下：

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

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

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

(4)循環執行(1)、(2)和(3)，直到整個模型最右邊一列，就可以得到炮點網格所在列右側每個節點的旅行時。

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

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

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

步驟6：炮點網格所在行的所有網格節點的旅行時已經得到更新，從上向下逐步計算炮點網格下側每一行上所有網格節點的最小旅行時，具體如下：

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

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

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

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

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

(4)循環執行(1)、(2)和(3),直到整個模型最下邊一行，得到炮點網格所在行下側每個節點的旅行時。

步驟7：計算炮點所在行上側網格節點的旅行時，計算方法與步驟6類似，計算完成后得到了整個模型區域中每個網格節點的最小旅行時。

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

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

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

步驟3：將步驟2得到的插值點坐標作為新的接收點，重復步驟1、步驟2直到插值點出現在震源所在網格的邊界或節點上為止(圖7b和圖7c)；

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

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

2.4 后續波的旅行時計算

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

2.4.1 模型分區

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

圖8 二維層狀介質模型多次反射、透射轉換波分區示意圖

2.4.2 多步計算

圖9所示為反射(或反射轉換)或透射(或透射轉換)波的波前擴展示意圖。具體實現過程為：

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

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

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

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

圖9 波前分區擴展示意圖

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

3 數值算例

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

3.1 均勻VTI介質模型

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

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

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

3.2 兩層水平層狀模型

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

表1 兩層水平層狀VTI介質模型參數

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

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

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

圖13 不同網格間距兩層水平層狀各向同性介質模型的分區多步LTI算法數值解的相對誤差

表2 兩層水平層狀各向同性介質模型參數

表3 不同網格間距分區多步FMM、ISPM

注：表中FMM和ISPM算法數據引自文獻[34]

3.3 二維起伏層狀VTI介質

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

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

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

圖15 單炮合成記錄

4 結論

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