起伏地表地震波旅行時(shí)混合網(wǎng)格線性插值射線追蹤計(jì)算方法

王 琦 朱 盼 葉 佩 李 勤 李慶春*

(①長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院，陜西西安 710054； ②中國(guó)國(guó)土資源航空物探遙感中心，北京 100083； ③西安科技大學(xué)地質(zhì)與環(huán)境學(xué)院，陜西西安 710054)

1 引言

旅行時(shí)線性插值(Linear traveltime interpolation 簡(jiǎn)稱LTI)是Asakawa等[1]提出的基于線性假設(shè)網(wǎng)格單元擴(kuò)展射線追蹤算法。由于該方法計(jì)算速度快、原理簡(jiǎn)單，得到了很大發(fā)展。趙改善等[2]結(jié)合界面二次源法將該方法推廣于追蹤反射波旅行時(shí)；聶建新等[3]將旅行時(shí)二次插值與線性插值方法聯(lián)合，降低了累積誤差；張賽民等[4]用拋物線插值取代線性插值，改善了因線性引起的誤差；張東等[5]為了求得最小走時(shí)，在向前處理過(guò)程中采用了多方向循環(huán)的計(jì)算方法。

近年來(lái)多次波越來(lái)越引起人們關(guān)注，對(duì)多次波的正確模擬和認(rèn)識(shí)是壓制多次波干擾的前提。目前對(duì)多次波的射線追蹤以基于網(wǎng)格單元擴(kuò)展的方法為主。Rawlinson等[6]利用快速步進(jìn)法結(jié)合分區(qū)多步計(jì)算技術(shù)實(shí)現(xiàn)了層狀介質(zhì)中多次透射波和多次反射波的計(jì)算，De Kool等[7]將此方法擴(kuò)展到三維球面坐標(biāo)系中；唐小平等[8]實(shí)現(xiàn)了分區(qū)多步改進(jìn)的最短路徑算法，可在二維和三維層狀介質(zhì)中追蹤透射、反射(或折射)和轉(zhuǎn)換的多次波；趙瑞等[9]提出了分區(qū)多步不規(guī)則最短路徑算法，也可實(shí)現(xiàn)多次波的追蹤計(jì)算，并與分區(qū)多步FMM算法進(jìn)行對(duì)比，認(rèn)為此方法比分區(qū)多步FMM算法具有高精度和高效率的優(yōu)勢(shì)。本文為了說(shuō)明分區(qū)多步LTI算法的優(yōu)勢(shì)，與文獻(xiàn)[9]中的分區(qū)多步FMM算法和分區(qū)多步ISPM算法的計(jì)算精度和效率亦做了比對(duì)。

關(guān)于起伏地表地震波傳播問(wèn)題，近年來(lái)國(guó)內(nèi)外學(xué)者開展了相關(guān)研究。岳玉波等[10]用插值函數(shù)對(duì)起伏界面進(jìn)行擬合，利用兩點(diǎn)射線追蹤方法實(shí)現(xiàn)了均勻速度模型中初值射線路徑計(jì)算，但這種方法計(jì)算效率較低； 孫章慶等[11]在矩形網(wǎng)格中針對(duì)起伏地表改進(jìn)了快速步進(jìn)法(FMM)，實(shí)現(xiàn)了起伏地表?xiàng)l件下的旅行時(shí)計(jì)算，此方法的網(wǎng)格仍然是單一的矩形網(wǎng)格，計(jì)算過(guò)程中如果網(wǎng)格過(guò)于稀疏會(huì)使界面呈階梯狀，引起較大誤差，加密網(wǎng)格又會(huì)使計(jì)算效率大打折扣且仍不能很好地逼近起伏地表或復(fù)雜構(gòu)造形態(tài)； 孫章慶等[12]用FMM法中窄帶技術(shù)代替?zhèn)鹘y(tǒng)的LTI法的波前擴(kuò)展方式實(shí)現(xiàn)了起伏地表射線追蹤，該方法采用三角網(wǎng)與矩形網(wǎng)相結(jié)合的方式對(duì)模型進(jìn)行剖分，使起伏地表下地震波走時(shí)的計(jì)算精度得到了提高，但并未涉及到多次波和轉(zhuǎn)換波走時(shí)的計(jì)算； Lelièvre 等[13]利用基于不規(guī)則網(wǎng)格的FMM法計(jì)算三維起伏地表模型的初至波旅行時(shí)，但該方法無(wú)論是在網(wǎng)格剖分還是在旅行時(shí)計(jì)算階段，計(jì)算量都非常大，而且此基于不規(guī)則網(wǎng)格的方法僅可以處理剖分三角形為銳角的情況。Sun等[14]提出在起伏地表附近采用間距不均一的網(wǎng)格剖分模型，用FMM法實(shí)現(xiàn)了起伏地表走時(shí)計(jì)算，但該方法需要對(duì)起伏地表的走時(shí)進(jìn)行插值處理，這樣勢(shì)必會(huì)引起誤差； Lan等[15]通過(guò)求解曲線坐標(biāo)系下的程函方程，實(shí)現(xiàn)了起伏地表下各向異性介質(zhì)的射線追蹤，該方法為精確求解，需要對(duì)求解區(qū)域進(jìn)行精細(xì)剖分，會(huì)嚴(yán)重影響計(jì)算效率； 趙后越等[16]利用規(guī)則網(wǎng)格和非規(guī)則網(wǎng)格相結(jié)合的方式進(jìn)行網(wǎng)格剖分，結(jié)合最短路徑法實(shí)現(xiàn)了各向異性起伏地表?xiàng)l件下初至波射線追蹤，此方法為了提高計(jì)算精度會(huì)在網(wǎng)格邊界加入新的節(jié)點(diǎn)，但同時(shí)降低了運(yùn)算效率。此外國(guó)內(nèi)外許多學(xué)者也對(duì)基于網(wǎng)格單元擴(kuò)展的射線追蹤算法做了有益的研究[17-26]。本文從網(wǎng)格剖分方式和局部走時(shí)計(jì)算兩個(gè)方面出發(fā)開展研究，提出一種能更好地逼近真實(shí)界面和更加精確高效地計(jì)算起伏地表下地震波旅行時(shí)的方法，文中討論了該方法的穩(wěn)定性，用理論模型試算說(shuō)明該方法具有良好的適應(yīng)性。

2 矩形網(wǎng)格LTI方法的基本原理

圖1 矩形網(wǎng)格局部旅行時(shí)計(jì)算示意圖 (a)射線從矩形網(wǎng)格的橫邊穿過(guò)； (b)射線從矩形網(wǎng)格的縱邊穿過(guò)

(1)

其中tb≥ta。根據(jù)幾何關(guān)系，C點(diǎn)的旅行時(shí)可以表示為

(2)

(3)

可解出

(4)

進(jìn)而可以求得D點(diǎn)的坐標(biāo)

(5)

則C點(diǎn)的旅行時(shí)可表示為

(6)

在矩形網(wǎng)格中所有邊界不是垂直就是水平的，當(dāng)點(diǎn)A、B位于橫邊時(shí)，za=zb，如圖1a所示，對(duì)應(yīng)的公式可簡(jiǎn)化為

(7)

(8)

(9)

當(dāng)點(diǎn)A、B位于縱邊時(shí)xa=xb，如圖1b所示，對(duì)應(yīng)的公式可簡(jiǎn)化為

(10)

(11)

(12)

僅考慮以上兩種情況，即可得到矩形網(wǎng)格LTI方法關(guān)于D點(diǎn)的坐標(biāo)與tc的計(jì)算公式。

LTI算法的實(shí)現(xiàn)過(guò)程可分為向前和向后兩個(gè)步驟進(jìn)行。

第一步(向前處理)：從炮點(diǎn)開始，據(jù)式(6)計(jì)算炮點(diǎn)到網(wǎng)格邊界上各個(gè)節(jié)點(diǎn)的最小旅行時(shí);

第二步(向后處理)：從檢波點(diǎn)開始，按式(5)計(jì)算所有射線路徑與各個(gè)網(wǎng)格邊界的交點(diǎn)。

以上每一步的實(shí)現(xiàn)過(guò)程均在文獻(xiàn)[1]中有詳細(xì)敘述，這里不再贅述。

3 混合網(wǎng)格LTI射線追蹤方法

混合網(wǎng)格LTI方法的實(shí)現(xiàn)過(guò)程可以分為三個(gè)步驟：①模型剖分：②計(jì)算局部旅行時(shí)；③計(jì)算目標(biāo)類型地震波旅行時(shí)。

3.1 模型剖分

為了更接近實(shí)際地表或速度界面，同時(shí)又能兼顧計(jì)算效率，本文采用矩形網(wǎng)格與不規(guī)則四邊形網(wǎng)格相結(jié)合的方法對(duì)模型進(jìn)行剖分，如圖2所示。首先對(duì)整個(gè)模型區(qū)域用矩形網(wǎng)格剖分，對(duì)地表和速度界面函數(shù)分別計(jì)算與網(wǎng)格縱邊的交點(diǎn)，得到離散點(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è)模型的剖分。

上述混合網(wǎng)格在對(duì)模型剖分時(shí)會(huì)遇到三種情況：①地表或速度界面3個(gè)相鄰的離散點(diǎn)都在原矩形網(wǎng)格的同一層內(nèi)；②地表或速度界面3個(gè)相鄰離散點(diǎn)中左端或右端的離散點(diǎn)與其余兩個(gè)不在原矩形網(wǎng)格的同一層內(nèi)；③地表或速度界面3個(gè)相鄰離散點(diǎn)中兩端的離散點(diǎn)與中間的離散點(diǎn)不在原矩形網(wǎng)格的同一層內(nèi)。以上所述的三種情況分別對(duì)應(yīng)于圖2中編號(hào)為①、②、③的陰影部分，為了便于觀察，將陰影部分按順序放大，如圖3所示。情況①時(shí)直接連接三個(gè)相鄰的離散點(diǎn)即可構(gòu)成不規(guī)則四邊形網(wǎng)格(圖3中的①)；情況②時(shí)首先連接在原矩形網(wǎng)格同一層的兩個(gè)離散點(diǎn)，然后連接地表或界面與縱邊界無(wú)交點(diǎn)的原矩形網(wǎng)格的左下(或右下或左上或右上)角節(jié)點(diǎn)與中間離散點(diǎn)，構(gòu)成非規(guī)則四邊形網(wǎng)格(圖3中的②)；情況③時(shí)分別連接地表或界面與縱邊界無(wú)交點(diǎn)的原矩形網(wǎng)格的左下(或右下或左上或右上)角節(jié)點(diǎn)、右下(或左下或左上或右上)角節(jié)點(diǎn)與中間離散點(diǎn)，構(gòu)成不規(guī)則四邊形網(wǎng)格(圖3中的③)。圖3中各個(gè)顏色不同的區(qū)域組成了利用上述方法建立的混合網(wǎng)格。

圖2和圖3中藍(lán)色實(shí)線代表單一矩形網(wǎng)格剖分時(shí)對(duì)地表或界面的逼近，紅色實(shí)線代表在混合網(wǎng)格剖分時(shí)逼近的地表或界面?？梢钥闯龊笳邔?duì)地表或界面的擬合程度更高，剖分方法也更加合理。

圖2 混合網(wǎng)中模型剖分與局部旅行時(shí)計(jì)算示意圖

3.2 混合網(wǎng)格LTI局部旅行時(shí)計(jì)算公式的建立

通過(guò)前兩節(jié)的分析得知，適用于矩形網(wǎng)格的局部旅行時(shí)計(jì)算公式并不適用于混合網(wǎng)格，因此首先要建立針對(duì)不規(guī)則四邊形的局部旅行時(shí)計(jì)算公式。

圖4 不規(guī)則四邊形網(wǎng)格中局部旅行時(shí)計(jì)算示意圖

設(shè)A、B點(diǎn)所在直線的方程為Ax+Bz+C=0，分別將A、B點(diǎn)的坐標(biāo)代入即可得

(za-zb)x-(xa-xb)z+za(xa-xb)-

xa(za-zb)=0

(13)

上式中tb>ta， Δt=tb-ta，求得td后，由幾何關(guān)系，C點(diǎn)的旅行時(shí)可以表示為

(14)

(15)

(16)

(17)

將式(16)、式(17)代入式(15)得

(18)

再將式(16)、式(17)、式(18)代入式(14)便得到了適用于不規(guī)則四邊形網(wǎng)格局部走時(shí)的計(jì)算公式

(19)

D點(diǎn)的坐標(biāo)可由A、B點(diǎn)的坐標(biāo)線性插值得到

(20)

以上即為不規(guī)則四邊形網(wǎng)格下求解tc和(xd,zd)的基本公式。

3.3 穩(wěn)定性分析

在3.2的推導(dǎo)過(guò)程中，從tc關(guān)于r的一階偏導(dǎo)數(shù)為零，即

(21)

可以得出： Δt的符號(hào)必須與d2-r的符號(hào)一致，否則在實(shí)數(shù)域內(nèi)無(wú)極值點(diǎn)存在。此外，D點(diǎn)必須介于A、B點(diǎn)之間或與B點(diǎn)(A點(diǎn))重合，否則計(jì)算得到的結(jié)果無(wú)實(shí)際物理意義，或與前提假設(shè)不符。記點(diǎn)D與點(diǎn)E之間的距離為l，則

(22)

理論上l的取值存在7種情況：

(1)xd∈(xa,xe)，且zd∈(za,ze)，則l∈(0,d2)；

(2)xd∈(xe,xb)，且zd∈(ze,zb)，則l∈(-d1+d2,0)；

(3)xd=xe，且zd=ze，則l=0；

(4)xd=xa，且zd=za，則l=d2；

(5)xd=xb，且zd=zb，則l=-(d1-d2)；

(23)

或

(24)

(25)

當(dāng)l∈(-∞,-d1+d2)時(shí)，求得

有

(26)

故后兩種情況無(wú)實(shí)際物理意義，在這里不予考慮。在式(19)、式(20)的基礎(chǔ)上針對(duì)以上前5種情況單獨(dú)考慮，可得到每種情況的表達(dá)式

(27)

(28)

(3)xd=xe，且zd=ze時(shí)，r=d2，此時(shí)Δt=0，代入式(19)、式(20)，得

tc=ta+sd3或tc=tb+sd3

(29)

(30)

(31)

(32)

(33)

(34)

以上敘述建立了不規(guī)則四邊形網(wǎng)格局部旅行時(shí)的計(jì)算公式，并且對(duì)其穩(wěn)定性做了分析。針對(duì)每種情況分別利用相應(yīng)的公式就可以穩(wěn)定且高效地計(jì)算出tc和(xd,zd)。

3.4 目標(biāo)波型的旅行時(shí)計(jì)算

本文提出的方法與分區(qū)多步計(jì)算技術(shù)[6]結(jié)合，就可以追蹤任意多次反射(或反射轉(zhuǎn)換)或透射(或透射轉(zhuǎn)換)波。旅行時(shí)的計(jì)算分為向前處理和向后處理兩個(gè)過(guò)程，在向前處理過(guò)程中計(jì)算每個(gè)網(wǎng)格節(jié)點(diǎn)的最小旅行時(shí)，向后處理過(guò)程是對(duì)射線路徑的追蹤以及計(jì)算整條射線路徑旅行時(shí)之和。在計(jì)算反射波或透射波的旅行時(shí)時(shí)需要引入分區(qū)多步的思想，即將整個(gè)模型區(qū)域按所需追蹤的波的類型和速度特征分成幾個(gè)獨(dú)立的計(jì)算區(qū)域，然后按照所分計(jì)算區(qū)域，在每個(gè)區(qū)域內(nèi)獨(dú)立進(jìn)行計(jì)算。要運(yùn)用分區(qū)多步計(jì)算技術(shù)，首先要建立模型并網(wǎng)格化(這里采用混合網(wǎng)格進(jìn)行剖分)，速度界面也要離散成不連續(xù)的點(diǎn)。

3.4.1 模型分區(qū)

圖5a所示為多次反射波在二維層狀介質(zhì)模型中的分區(qū)情況，圖中展示了一條在界面折返三次的多次反射波。根據(jù)分區(qū)多步的思想，該多次反射波旅行時(shí)的計(jì)算過(guò)程可以在四個(gè)獨(dú)立的計(jì)算區(qū)域(圖中1、2、3、4代表4個(gè)分區(qū))內(nèi)進(jìn)行。具體實(shí)現(xiàn)過(guò)程如下。該多次反射波由四段組成： 第一段是來(lái)自地表到第二個(gè)反射界面的的下行波，在計(jì)算時(shí)將第一、二層作為第一個(gè)分區(qū)來(lái)計(jì)算第一段下行波的波前； 第二段是自第二個(gè)反射界面上的反射點(diǎn)到第一個(gè)反射界面的上行波，在計(jì)算時(shí)將第二層作為第二個(gè)分區(qū)來(lái)計(jì)算第二段上行波的波前；第三段是自第一個(gè)反射界面反射點(diǎn)到第二個(gè)反射界面的下行波，在計(jì)算時(shí)將第二層作為第三個(gè)分區(qū)來(lái)計(jì)算第三段下行波的波前；第四段是自第二個(gè)反射界面反射點(diǎn)到地表的上行波，在計(jì)算時(shí)將第一、二層作為第四個(gè)分區(qū)來(lái)計(jì)算第四段上行波的波前。

圖5b所示為多次透射轉(zhuǎn)換波在二維層狀介質(zhì)模型中的分區(qū)情況，圖中展示了一條在界面兩次透射、一次反射、三次轉(zhuǎn)換(射線每透過(guò)一次界面或每反射一次就發(fā)生一次轉(zhuǎn)換)的轉(zhuǎn)換波射線。按照分區(qū)多步思想，該多次透射轉(zhuǎn)換波旅行時(shí)的計(jì)算過(guò)程可以在四個(gè)獨(dú)立的計(jì)算區(qū)域(圖中1、2、3、4代表4個(gè)分區(qū))內(nèi)進(jìn)行。假設(shè)該射線初始入射為P波，則第一段是自地面到第一個(gè)透射界面的下行P波，在計(jì)算時(shí)將第一層作為第一個(gè)分區(qū)來(lái)計(jì)算第一段下行P波的波前，此時(shí)計(jì)算所用速度為第一層的P波速度；第二段是自第一層透射點(diǎn)到第二層反射界面的下行S波，在計(jì)算時(shí)將第二層作為第二個(gè)分區(qū)來(lái)計(jì)算第二段下行S波的波前，此時(shí)計(jì)算所用速度為第二層的S波速度；第三段是自第二層反射點(diǎn)到第一層透射界面的上行P波，在計(jì)算時(shí)將第二層作為第三個(gè)分區(qū)來(lái)計(jì)算第三段上行P波的波前，此時(shí)計(jì)算所用速度為第二層的P波速度；第四段是自第一層透射點(diǎn)到地表的上行S波，在計(jì)算時(shí)將第一層作為第四個(gè)分區(qū)來(lái)計(jì)算第四段上行S波的波前，此時(shí)計(jì)算所用速度為第一層的S波速度。這樣就可以追蹤P-P波、P-S波、S-S波、S-P波。

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

3.4.2 多步計(jì)算

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

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

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

(1)計(jì)算由震源到第一個(gè)分區(qū)下界面的下行波的波前，并找出界面離散點(diǎn)上走時(shí)最小的點(diǎn)，如果是P波入射，則計(jì)算時(shí)用P波速度，如果是S波入射，則計(jì)算時(shí)用S波速度，如圖6a所示。

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

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

(4)由以上三個(gè)步驟計(jì)算波前，再結(jié)合LTI的向后處理就可以追蹤各類型波的射線路徑。

4 數(shù)值算例

為了驗(yàn)證本文提出的計(jì)算方法的計(jì)算精度和效率，首先將分區(qū)多步LTI算法分別與分區(qū)多步FMM算法和分區(qū)多步ISPM算法的計(jì)算精度及效率進(jìn)行比對(duì)，再選取三個(gè)速度模型，分別為起伏地表一維線性增加速度模型、帶有垂直斷層的層狀介質(zhì)模型和加上起伏地表與速度界面的的等價(jià)Marmousi模型，進(jìn)行試算。對(duì)第一個(gè)模型分別用分區(qū)多步矩形網(wǎng)格和混合網(wǎng)格LTI法計(jì)算多次反射波的旅行時(shí)，對(duì)比分析算法的計(jì)算精度。對(duì)第二個(gè)模型分別追蹤層間多次反射波和多次透射轉(zhuǎn)換波，對(duì)第三個(gè)模型追蹤其層間多次反射波，合成各個(gè)目標(biāo)波形的地震記錄。

4.1 方法對(duì)比及精度分析

目前基于網(wǎng)格單元波前擴(kuò)展的多次波波前追蹤方法中較為成熟的主要有分區(qū)多步FMM算法和分區(qū)多步ISPM算法。為了驗(yàn)證分區(qū)多步LTI算法的計(jì)算精度和計(jì)算效率，對(duì)同一個(gè)速度模型分別用以上三種方法，在4種網(wǎng)格間距下計(jì)算反射波走時(shí)。模型尺寸為100km×40km，速度為4000m/s，地表水平，在30km處有一條水平反射界面。在模型參數(shù)化時(shí)，分區(qū)多步FMM算法中正方形網(wǎng)格尺寸分為四種尺度，分別為1000m×1000m、500m×500m、250m×250m和125m×125m；分區(qū)多步ISPM算法中網(wǎng)格尺寸為4000m×4000m，并在網(wǎng)格單元邊界上相應(yīng)加入3、7、15和31個(gè)次級(jí)節(jié)點(diǎn)；分區(qū)多步LTI算法中網(wǎng)格大小和分區(qū)多步FMM算法中設(shè)置相同，這樣保證了以上三種方法中網(wǎng)格節(jié)點(diǎn)的間距相同。炮點(diǎn)在模型的左上角坐標(biāo)為(0，0)，100個(gè)檢波器等間距(1000m)排列在地面，最小炮檢距為1000m。圖7和圖8分別給出了隨炮檢距變化的分區(qū)多步FMM算法(圖7a)，分區(qū)多步ISPM算法(圖7b)和分區(qū)多步LTI算法(圖8)相對(duì)于解析解在4種不同節(jié)點(diǎn)距下反射波走時(shí)的相對(duì)誤差。表1給出了三種方法分別所用CPU耗時(shí)。從圖7、圖8和表1中可見，除了網(wǎng)格間距為1000m×1000m的情形外分區(qū)多步LTI算法的精度均優(yōu)于分區(qū)多步FMM和分區(qū)多步ISPM算法，在計(jì)算效率方面，分區(qū)多步LTI算法明顯優(yōu)于分區(qū)多步FMM算法而低于分區(qū)多步ISPM算法。

圖7 隨著炮檢距變化的FMM算法
(a)和ISPM算法(b)數(shù)值解相對(duì)于解析解在4種不同網(wǎng)格距下反射波旅行時(shí)的相對(duì)誤差(引自文獻(xiàn)[9])

圖8 隨著炮檢距變化的分區(qū)多步LTI算法數(shù)值解相對(duì)于解析解在4種不同網(wǎng)格距下反射波旅行時(shí)的相對(duì)誤差

表1 均勻速度模型中FMM、ISPM和分區(qū)多步LTI算法 在4種不同網(wǎng)格間距下CPU耗時(shí)對(duì)比表

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

圖9a是一個(gè)大小為400m×400m、每一層的速度皆為一維線性增加的三層速度模型，其速度表達(dá)式如圖9a中所列，其中h1、h2、h3表示每一層的厚度。炮點(diǎn)置于(200m，-20m)處，在地表上布設(shè)101個(gè)檢波點(diǎn)，橫向道間距均為4m，模型中剖分網(wǎng)格尺寸(在混合網(wǎng)格中僅指橫向網(wǎng)格間距)分別為0.1m×0.1m、0.2m×0.2m、0.4m×0.4m、1.0m×1.0m、2.0m×2.0m和4.0m×4.0m， 分別利用本文提出的方法和完全矩形網(wǎng)剖分的方法計(jì)算自炮點(diǎn)經(jīng)界面三次折返到起伏地表檢波點(diǎn)的射線路徑的旅行時(shí)。

Kim[27]通過(guò)改變節(jié)點(diǎn)間距的方法來(lái)做誤差分析，此方法的前提是誤差收斂。認(rèn)為隨著網(wǎng)格間距的縮小，誤差也會(huì)隨之減小，最終收斂于最小值；反之，隨著網(wǎng)格間距增大，誤差也會(huì)隨之向上浮動(dòng)。

本文采用此方法進(jìn)行誤差分析。為了驗(yàn)證本文所采用的方法比矩形網(wǎng)LTI法在進(jìn)行多次反射波射線追蹤時(shí)的數(shù)值穩(wěn)定性以及射線路徑的可靠性更好，本文通過(guò)改變網(wǎng)格單元大小進(jìn)行誤差分析，以最小網(wǎng)格單元尺寸為0.1m×0.1m的結(jié)果作為參考值，計(jì)算出其他五種網(wǎng)格單元尺寸相對(duì)于單元尺寸為0.1m×0.1m時(shí)各檢波點(diǎn)的多次反射波旅行時(shí)的相對(duì)誤差。

從圖9b(采用分區(qū)多步混合網(wǎng)LTI算法)、圖9c(采用分區(qū)多步矩形網(wǎng)LTI算法)和表2中可以看出：兩種算法計(jì)算多次反射波旅行時(shí)所得的誤差都隨網(wǎng)格大小縮小而逐漸收斂。但相比較而言，利用分區(qū)多步混合網(wǎng)格LTI算法所計(jì)算的旅行時(shí)的誤差總體上要低于分區(qū)多步矩形網(wǎng)格LTI算法，除網(wǎng)格尺寸為4.0m×4.0m以外，其余網(wǎng)格尺寸下計(jì)算的相對(duì)誤差基本在0.3%以內(nèi)，且在相同網(wǎng)格尺寸下，分區(qū)多步混合網(wǎng)格LTI算法的計(jì)算精度高于分區(qū)多步矩形網(wǎng)格LTI算法，在保證精度相當(dāng)?shù)臈l件下，分區(qū)多步混合網(wǎng)格LTI算法的計(jì)算效率高于分區(qū)多步矩形網(wǎng)格LTI算法。據(jù)此說(shuō)明分區(qū)多步混合網(wǎng)格LTI算法是穩(wěn)定的，進(jìn)行多次反射波模擬的射線路徑是可靠的。

4.2 二維層狀介質(zhì)多次反射波模擬

圖10為在(200m，-20m)處單炮激發(fā)、11道接收的多次反射波波前在每個(gè)計(jì)算區(qū)域內(nèi)的擴(kuò)展過(guò)程和多次波的射線路徑。從圖10b～圖10e波前在各個(gè)計(jì)算區(qū)域中的擴(kuò)展過(guò)程和方式可以看出，在速度均勻的同一區(qū)域內(nèi)波前大致呈相互平行的弧線。當(dāng)波前穿過(guò)速度分界面時(shí)，就會(huì)隨速度的改變發(fā)生變化，這是由速度分界面兩邊的速度變化引起的，這種現(xiàn)象符合地震波在層狀介質(zhì)中的傳播規(guī)律。圖10g為根據(jù)圖10a模型模擬計(jì)算的單炮合成地震記錄，在合成記錄中分別顯示了直達(dá)波、來(lái)自Ⅰ與Ⅱ波阻抗界面的一次反射P波以及多次反射P波。從圖10g中可以看出， ②與②*相比完全不一致，如A、B兩處相差較大，這是由于檢波點(diǎn)分別位于水平地表和起伏地表引起的，即如果將反射波靜校正到虛擬的水平地表(浮動(dòng)基準(zhǔn)面)，地震記錄就會(huì)出現(xiàn)上述問(wèn)題，不能客觀地反映界面的實(shí)際位置。

表2 圖9a所示模型由分區(qū)多步混合網(wǎng)格LTI和分區(qū)多步矩形網(wǎng)LTI計(jì)算所得的 多次反射波旅行時(shí)的誤差、相對(duì)誤差和CPU在5種不同網(wǎng)格尺寸下對(duì)比表

圖10 層間多次波波前在各分區(qū)中傳播過(guò)程及射線路徑(波前間隔為0.01s)

(a)模型參數(shù)和分區(qū)編號(hào)； (b)下行P波波前在第一個(gè)分區(qū)內(nèi)傳播； (c)上行P波波前在第二個(gè)分區(qū)內(nèi)傳播； (d)下行P波波前在第三個(gè)分區(qū)內(nèi)傳播； (e)上行P波波前在第四個(gè)分區(qū)內(nèi)傳播； (f)層間多次波射線路徑； (g)單炮合成記錄

4.3 二維層狀介質(zhì)多次透射、轉(zhuǎn)換、一次反射波模擬

圖11為在(200m，-20m)處單炮激發(fā)、11道接收的多次透射轉(zhuǎn)換波波前在各個(gè)分區(qū)中的傳播過(guò)程以及多次透射轉(zhuǎn)換波的射線路徑。從圖11b～圖11e中波前在各個(gè)分區(qū)中的傳播過(guò)程可以看出，當(dāng)波前以橫波速度傳播時(shí)，波前的等時(shí)線較密(如圖11c、圖11e)，當(dāng)波前以縱波速度傳播時(shí)波前等時(shí)線相對(duì)較為稀疏(如圖11b、圖11d)，這是因?yàn)闄M波傳播速度相對(duì)于縱波傳播速度較慢，符合地震波在介質(zhì)中的傳播規(guī)律。從圖11f中可以看出：在各個(gè)速度均勻的分層中，射線路徑是直線，但是通過(guò)速度分界面時(shí)射線就會(huì)發(fā)生偏折。若入射的是S波，透射或反射為P波時(shí)，相應(yīng)的出射角或反射角會(huì)增大；若入射的是P波，透射或反射為S波，相應(yīng)的出射角或反射角會(huì)減小，這符合斯奈爾定律?？梢哉f(shuō)明在計(jì)算時(shí)分區(qū)徹底，計(jì)算方法正確，方法靈活性較強(qiáng)，可以追蹤到層狀介質(zhì)中各種透射、轉(zhuǎn)換、反射波的射線路徑和旅行時(shí)。圖11g為在圖11a中模擬的單炮合成記錄。記錄中分別顯示了各道直達(dá)波、分別來(lái)自Ⅰ與Ⅱ波阻抗界面的P-P波、P-S波、圖10f中所示的多次反射P波以及圖11f中所示的多次透射轉(zhuǎn)換波。

4.4 Marmousi模型起伏檢波點(diǎn)反射波及多次波模擬

圖12為在人為劃定劇烈起伏地表和兩條反射界面(圖12a中加粗的黑色實(shí)線)后的Marmousi模型中，在(2500m，-500m)處單炮激發(fā)、11道接收的多次反射P波波前在各個(gè)分區(qū)中的傳播及射線路徑。由圖12b～圖12e中可以看出，波前在擴(kuò)展過(guò)程中在速度較高的區(qū)域等時(shí)線被拉伸，在速度較低的區(qū)域等時(shí)線被壓縮，且波前等時(shí)線連續(xù)，滿足地震波的傳播規(guī)律，說(shuō)明基于混合網(wǎng)格剖分的LTI方法對(duì)構(gòu)造復(fù)雜的速度模型有很強(qiáng)的適應(yīng)能力，適用于復(fù)雜速度模型的多次波模擬。圖12g給出了在圖12a的模型中模擬的單炮合成記錄。在合成記錄中可以看到初至波、分別來(lái)自于第Ⅰ與第Ⅱ波阻抗界面的一次反射P波和圖12f所示的多次反射P波。

圖11 多次透射轉(zhuǎn)換波波前在各分區(qū)中的傳播過(guò)程及射線路徑(波前間隔為0.01s) (a)模型參數(shù)和分區(qū)編號(hào)； (b)下行P波波前在第一個(gè)分區(qū)內(nèi)傳播； (c)下行S波波前在第二個(gè)分區(qū)內(nèi)傳播；(d)上行P波 波前在第三個(gè)分區(qū)內(nèi)傳播； (e)上行S波波前在第四個(gè)分區(qū)內(nèi)傳播； (f)多次透射轉(zhuǎn)換波射線路徑； (g)單炮合成記錄

(a)模型參數(shù)和分區(qū)編號(hào)； (b)下行P波波前在第一個(gè)分區(qū)內(nèi)傳播； (c)上行P波波前在第二個(gè)分區(qū)內(nèi)傳播； (d)下行P波波前在第三個(gè)分區(qū)內(nèi)傳播； (e)上行P波波前在第四個(gè)分區(qū)內(nèi)傳播； (f)層間多次波射線路徑; (g)帶地形的Marmousi模型單炮合成記錄

5 結(jié)論

本文提出了用在矩形中加入不規(guī)則四邊形網(wǎng)格后形成的混合網(wǎng)格來(lái)取代單一矩形網(wǎng)格對(duì)模型進(jìn)行剖分，提高了對(duì)起伏地表和速度界面的處理能力，在矩形網(wǎng)格局部旅行時(shí)計(jì)算公式的基礎(chǔ)上，建立了適用于混合網(wǎng)格局部旅行時(shí)計(jì)算的公式，并證明了該公式的穩(wěn)定性。利用分區(qū)多步混合網(wǎng)格LTI法實(shí)現(xiàn)了相關(guān)類型地震波(初至波、一次反射波、多次反射波、多次透射轉(zhuǎn)換波)旅行時(shí)的計(jì)算和射線路徑的追蹤。

通過(guò)對(duì)帶地形的層狀模型和速度不均勻模型的試算，不但說(shuō)明了該方法的計(jì)算結(jié)果精度高于同等網(wǎng)格間距下的矩形網(wǎng)格LTI方法，而且證明了該方法對(duì)復(fù)雜速度模型具有良好的適應(yīng)能力。為起伏地表和復(fù)雜構(gòu)造情況下地震波旅行時(shí)的計(jì)算提供了一種新的方法。

[1] Asakawa E and Kawanaka T.Seismic ray tracing using linear traveltime interpolation.Geophysical Prospecting,1993,41(4):99-111．

[2] 趙改善,郝守玲.基于旅行時(shí)線性插值的地震射線追蹤算法．石油物探,1998,32(2):14-24. Zhao Gaishan,Hao Shouling.Seismic ray tracing algorithm based on the linear traveltime interpolation.GPP,1998,32(2):14-24．

[3] 聶建新,楊慧珠.地震波旅行時(shí)二次/線性聯(lián)合插值法．清華大學(xué)學(xué)報(bào),2003,43(11):1495-1498. Nie Jianxin,Yang Huizhu.Quadrarical linear travel time interpolation of seismic ray-tracing.Journal of Tsinghua University,2003,43(11):1495-1498．

[4] 張賽民,周竹生,陳靈君等.對(duì)旅行時(shí)進(jìn)行拋物型插值的地震射線追蹤方法．地球物理學(xué)進(jìn)展,2007,22(1):43-48. Zhang Saimin,Zhou Zhusheng,Chen Lingjun et al.Seismic ray tracing method of applying parabolic interpolation to travel time.Progress in Geophysics,2007,22(1):43-48．

[5] 張東,謝寶蓮,楊艷等.一種改進(jìn)的線性旅行時(shí)插值射線追蹤算法．地球物理學(xué)報(bào),2009,52(1):200-205. Zhang Dong,Xie Baolian,Yang Yan et al.A ray tra-cing method based on improved linear traveltime interpolation.Chinese Journal of Geophysics,2009,52(1):200-205．

[6] Rawlinson N,Sambridge M.Multiple reflection and transmission phases in complex layered media using multistage fast marching method.Geophysics,2004,69(5):1338-1350．

[7] De Kool M,Rawlinson N,Sambridge M.A practical grid based method for tracking multiple refraction and reflection phases in three-dimensional heterogeneous media.Geophysical Journal International,2006,167(1):253-270．

[8] 唐小平,白超英.最短路徑算法下三維層狀介質(zhì)中多次波追蹤．地球物理學(xué)報(bào),2009,52(10):2635-2643. Tang Xiaoping,Bai Chaoying.Multiple ray tracing within 3-d layered media with the shortest path me-thod.Chinese Journal of Geophysics,2009,52(10):2635-2643．

[9] 趙瑞,白超英.復(fù)雜層狀模型中多次波快速追蹤——一種基于非規(guī)則網(wǎng)格的最短路徑算法．地震學(xué)報(bào),2010,32(4):433-444. Zhao Rui,Bai Chaoying.Fast multiple ray tracing within complex layered media: the shortest path method based on irregular grid cells.Acta Seismologica Sinica,2010,32(4):433-444．

[10] 岳玉波.起伏地表下初值射線追蹤的實(shí)現(xiàn)．勘探地球物理進(jìn)展,2007,30(5):388-391. Yue Yubo.Realization of initial value ray tracing for rugged topography.Progress in Exploration Geophy-sics,2007,30(5):388-391．

[11] 孫章慶,孫建國(guó).波前快速推進(jìn)法起伏地表地震波走時(shí)計(jì)算．勘探地球物理進(jìn)展,2007,30(5):92-97. Sun Zhangqing,Sun Jianguo.Traveltime computation using fast marching method from rugged topography.Progress in Exploration Geophysics,2007,30(5):92-97．

[12] 孫章慶,孫建國(guó),韓復(fù)興.復(fù)雜地表?xiàng)l件下基于線性插值和窄帶技術(shù)的地震波走時(shí)計(jì)算．地球物理學(xué)報(bào),2009,52(11):2846-2853. Sun Zhangqing,Sun Jianguo and Han Fuxing.Traveltimes computation using linear interpolation and narrow band technique under complex topographical conditions.Chinese Journal of Geophysics,2009,52(11):2846-2853．

[13] Lelièvre P G,Farquharson C G and Hurich C A.Computing first-arrival seismic traveltimes on unstructured 3-D tetrahedral grids using the fast marching method.Geophysical Journal International,2011,184(2)：885-896．

[14] Sun Jianguo,Sun Zhangqing and Han Fuxing.A finite difference scheme for solving the eikonal equation including surface to topography.Geophysics,2011,76(4):53-56．

[15] Lan Haiqing and Zhang Zhongjie.Topography-dependent eikonal equation and its solver for calculating first-arrival traveltimes with an irregular surface.Geophysical Journal International,2013,193(2):1010-1026．

[16] 趙后越,張美根.起伏地表?xiàng)l件下各向異性地震波最短路徑射線追蹤．地球物理學(xué)報(bào),2014,57(9):2910-2917. Zhao Houyue,Zhang Meigen.Tracing seismic shortest path rays in anisotropic medium with rolling surface.Chinese Journal of Geophysics,2014,57(9):2910-2917．

[17] 孫章慶.起伏地表?xiàng)l件下的地震波走時(shí)與射線路徑計(jì)算[學(xué)位論文].吉林長(zhǎng)春:吉林大學(xué),2011. Sun Zhangqing.The Seismic Traveltimes and Ray Path Computation under Undulating Earth’s Surface Condition[D].Jilin University,Changchun,Jilin,2011．

[18] 李永博,李慶春,吳瓊等.快速行進(jìn)法射線追蹤提高旅行時(shí)計(jì)算精度和效率的改進(jìn)措施．石油地球物理勘探,2016,51(3):467-473. Li Yongbo,Li Qingchun,Wu Qiong et al.Improved fast marching method for higher calculation accuracy and efficiency of traveltime.OGP,2016,51(3):467-473．

[19] 王珺,李永慶．遺傳算法和正交時(shí)頻原子相結(jié)合的地震記錄快速匹配追蹤．石油地球物理勘探,2016,51(5):881-888,893. Wang Jun,Li Yongqing.Seismic trace fast matching pursuit based on genetic algorithm and orthogonal time-frequency atom.OGP,2016,51(5):881-888,893．

[20] 秦寧,王延光,楊曉東等．非水平地表高斯束疊前深度偏移及山前帶應(yīng)用實(shí)例．石油地球物理勘探,2017,52(1):81-86. Qi Ning,Wang Yanguang,Yang Xiaodong et al.Gaussian bean prestack depth migration for undula-ting-surface area in piedmont zone.OGP,2017,52(1):81-86．

[21] 黃國(guó)嬌,白超英.二維復(fù)雜層狀介質(zhì)中地震多波旅行時(shí)聯(lián)合反演成像．地球物理學(xué)報(bào),2010,53(12): 2972-2981. Huang Guojiao,Bai Chaoying.Simultaneous inversion with multiple traveltimes within 2-d complex layered media.Chinese Journal of Geophysics,2010,53(12):2972-2981．

[22] 桑運(yùn)云,孫曉軍.起伏地表下基于拋物插值的最短路徑射線追蹤．石油物探,2014,53(2):1441-1448. Sang Yunyun,Sun Xiaojun.Shortest path raytracing based on parabolic traveltime interpolation in irregular topography.GPP,2014,53(2):1441-1448．

[23] 侯爵,張忠杰,蘭海強(qiáng)等.起伏地表下地震波傳播數(shù)值模擬方法研究進(jìn)展．地球物理學(xué)進(jìn)展,2014,29(2):488-497. Hou Jue,Zhang Zhongjie,Lan Haiqiang et al.Progress in numerical simulation of seismic wave propagation under an undulating surface.Progress in Geophysics,2014,29(2)：488-497．

[24] 孫章慶,孫建國(guó),韓復(fù)興.基于線性插值和窄帶技術(shù)的走時(shí)計(jì)算方法．石油地球物理勘探,2009,44(4):436-441. Sun Zhangqing,Sun Jianguo and Han Fuxing.Travel-time computation based on linear interpolation and narrow-band technique.OGP,2009,44(4):436-441．

[25] 葉佩,李慶春.旅行時(shí)線性插值射線追蹤提高計(jì)算精度和效率的改進(jìn)方法．吉林大學(xué)學(xué)報(bào)(地球科學(xué)版),2013,43(1):291-298. Ye Pei,Li Qingchun.Improvements of linear traveltime interpolation ray tracing for the accuracy and efficiency.Journal of Jilin University (Earth Science Edition),2013,43(1):291-298．

[26] 宋御杰.基于混合網(wǎng)格剖分的起伏地表旅行時(shí)計(jì)算方法研究[學(xué)位論文].陜西西安: 長(zhǎng)安大學(xué),2015. Song Yujie.The Travel-time Computation Method Study Under Undulating Earth’s Surface Based on Hybrid mesh[D].Chang’an University,Xi’an,Shaanxi,2015．

[27] Kim S.3-D eikonal solvers: First-arrival traveltimes.Geophysics,2002,67(4):1225-1231.