2D/3D 起伏地表多震相地震波走時(shí)的因式分解程函方程算法

張?jiān)?，李夕海＊，白超英，牛超，王藝婷，曾小?/p>

（1.火箭軍工程大學(xué)，陜西西安 710025； 2.長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院地球物理系，陜西西安 710054； 3.長(zhǎng)安大學(xué)計(jì)算地球物理研究所，陜西西安 710054）

0 引言

基于幾何地震學(xué)的地震波走時(shí)計(jì)算技術(shù)具有算法便捷、計(jì)算效率高等優(yōu)點(diǎn)而被廣泛應(yīng)用于地震層析成像［1-4］、地震定位［5--7］等領(lǐng)域。地震波走時(shí)計(jì)算方法主要可分為三種：一是基于運(yùn)動(dòng)學(xué)方程的射線追蹤類(lèi)方法［8-9］，但當(dāng)遇到復(fù)雜介質(zhì)時(shí)，這類(lèi)方法容易陷入局域解，存在陰影區(qū)問(wèn)題； 二是基于費(fèi)馬原理的最短路徑算法（Shorest Path Method，SPM）［10-11］，這類(lèi)方法通過(guò)采用大量的網(wǎng)格線段近似逼近地震波射線，為保證結(jié)果精度，需用比其他方法更龐大數(shù)量的網(wǎng)絡(luò)節(jié)點(diǎn)，這也導(dǎo)致了更高的計(jì)算空間要求； 三是基于有限差分求解程函方程類(lèi)方法，它有效避免了運(yùn)動(dòng)學(xué)方程類(lèi)方法的不足，且在計(jì)算精度與計(jì)算效率方面具有顯著優(yōu)勢(shì)，自提出以來(lái)即發(fā)展迅速，作為正演算法多次解決地震波走時(shí)反演問(wèn)題［12-14］，近年有不少學(xué)者將其推廣至各向異性介質(zhì)［15-20］。

快速行進(jìn)法（Fast Marching Method，F(xiàn)MM）［21］與快速掃描法（Fast Sweeping Method，F(xiàn)SM）［22］均是通過(guò)有限差分求解程函方程，從而得到全局地震波走時(shí)的計(jì)算方法。FMM 采用窄帶技術(shù)近似波前擴(kuò)展，即隱含了這樣一個(gè)假設(shè)：群速度給出的能量傳播方向與相速度給出的波前傳播方向相同。然而，該假設(shè)僅對(duì)各向同性介質(zhì)是成立的，對(duì)于各向異性介質(zhì)則不成立。雖然近年來(lái)也有不少學(xué)者通過(guò)加入其他約束性條件將其推廣至各向異性介質(zhì)［23-25］，但過(guò)多的約束條件也使得算法變得復(fù)雜且低效。因此，對(duì)于各向異性介質(zhì)而言，人們更愿意采用FSM。FSM 由Zhao［22］率先提出，通過(guò)沿交替方向掃描計(jì)算域求解程函方程，在有限次數(shù)的迭代中收斂得到數(shù)值解，比FMM 更靈活且穩(wěn)健性更高。有關(guān)FMM 與FSM 的詳細(xì)比較，可參考文獻(xiàn)［26］。盡管FSM 在各向異性介質(zhì)及穩(wěn)健性上明顯優(yōu)于FMM，但研究表明［27-28］，對(duì)于復(fù)雜模型或網(wǎng)格節(jié)點(diǎn)較多的模型而言，F(xiàn)SM 計(jì)算效率低于FMM。主要原因是對(duì)于較復(fù)雜模型FSM 需更多的迭代次數(shù)才能收斂，對(duì)于具有大量網(wǎng)格節(jié)點(diǎn)的復(fù)雜模型，此問(wèn)題更嚴(yán)重。如Li等［29］為保證反演效率，選擇FMM 作為正演算法以解決反演算法的效率問(wèn)題；Benaichouche 等［30］則選用FSM 作為正演地震波走時(shí)計(jì)算方法，然而在反演計(jì)算雅可比矩陣時(shí)，選用FMM求解原始程函方程，這樣雖然提高了效率，但在一定程度上犧牲了計(jì)算精度。

FMM 是通過(guò)有限差分求解程函方程，采用窄帶近似波前而擴(kuò)展起來(lái)的一種方法，Rawlinson 等［31］將FMM 與分區(qū)多步計(jì)算技術(shù)（Multistage Computational Technique）相結(jié)合，實(shí)現(xiàn)了2D 復(fù)雜地質(zhì)結(jié)構(gòu)下的多震相地震波走時(shí)計(jì)算（Multistage FMM）。隨后，de Kool 等［32］考慮地球曲率影響，又將該算法推廣至球坐標(biāo)系下進(jìn)行多震相地震波走時(shí)計(jì)算。孫章慶等［33］基于FMM，提出不等距迎風(fēng)差分格式（Uneven Grid Upwind Finite-Difference Formula）計(jì)算3D 起伏地表?xiàng)l件的直達(dá)波走時(shí)，隨后又引入群行進(jìn)法（Group Marching Method）［34］以提高FMM 計(jì)算效率［35］。然而，震源點(diǎn)奇異性問(wèn)題一直伴隨著FMM 的發(fā)展。為解決此問(wèn)題，Rawlinson 等［31］采用源點(diǎn)網(wǎng)格加密技術(shù)，分析表明，經(jīng)過(guò)對(duì)震源附近網(wǎng)格細(xì)化后的FMM 可明顯降低源點(diǎn)附近的計(jì)算誤差，但相對(duì)于整個(gè)計(jì)算區(qū)域，源點(diǎn)附近的計(jì)算誤差仍然較大，因此震源區(qū)網(wǎng)格細(xì)化并不能根本解決震源點(diǎn)的奇異性問(wèn)題。此外，震源區(qū)網(wǎng)格細(xì)化的同時(shí)也意味著額外增加了計(jì)算量。Alkhalifah 等［36］則在球坐標(biāo)系下求解程函方程，但這種算法復(fù)雜，實(shí)際應(yīng)用中難以推廣。Fomel等［37］提出將地震波走時(shí)因式分解，用FSM 求解因式分解程函方程以解決震源奇異性問(wèn)題，數(shù)值模擬表明，求解因式分解程函方程數(shù)值解的準(zhǔn)確性高于直接計(jì)算原始程函方程，且震源附近改善效果更明顯。周小樂(lè)等［38］推導(dǎo)建立了曲線坐標(biāo)系因式分解程函方程，結(jié)合FSM，在解決震源奇異性問(wèn)題的同時(shí)，可得到起伏地表?xiàng)l件的地震波走時(shí)。

鑒于FMM 對(duì)復(fù)雜模型計(jì)算效率的優(yōu)勢(shì)及求解因式分解程函方程可解決震源奇異性問(wèn)題，本文提出了一種起伏地表?xiàng)l件下基于FMM 求解因式分解程函方程的數(shù)值算法，同樣可從根源上解決震源奇異性問(wèn)題。與傳統(tǒng)FMM 相同，本文所提算法可在保證O（NlgN，N為網(wǎng)絡(luò)節(jié)點(diǎn)數(shù)） 計(jì)算成本內(nèi)使用一階或二階有限差分公式求解因式分解程函方程。同時(shí)，考慮到加入后續(xù)震相約束可顯著提高反演分辨率，本文又結(jié)合分區(qū)多步計(jì)算技術(shù)實(shí)現(xiàn)了多震相地震波走時(shí)的計(jì)算。

文中首先闡釋如何建立因式分解程函方程； 然后介紹利用FMM 求解因式分解程函方程； 最后通過(guò)2D/3D 均勻模型、速度線性增加模型和起伏界面模型驗(yàn)證了所提算法的有效性與魯棒性。為表述方便，將本文所提算法稱(chēng)為基于因式分解程函方程FMM，簡(jiǎn)稱(chēng)FMM-F。

1 因式分解程函方程的導(dǎo)出

限于篇幅，在推導(dǎo)因式分解程函方程過(guò)程中，僅展示2D 直角坐標(biāo)系情形。對(duì)于3D，只需進(jìn)行類(lèi)似的擴(kuò)展即可。首先引入2D 直角坐標(biāo)系程函方程

式中T（x，z）和S（x，z）均為空間位置坐標(biāo)（x，z）的函數(shù)，分別表示震源點(diǎn)O到空間點(diǎn)（x，z）的地震波走時(shí)與該點(diǎn)的波慢度。將T（x，z）分解為一個(gè)關(guān)于點(diǎn)（x，z）到震源點(diǎn)O的距離函數(shù)T0與一個(gè)地震波走時(shí)擾動(dòng)量T1的乘積形式，即

這樣，根據(jù)鏈?zhǔn)角髮?dǎo)法則，標(biāo)準(zhǔn)程函方程式（1）便可寫(xiě)成

式中T0可通過(guò)解析求得，T1則需通過(guò)數(shù)值求解得到，再通過(guò)式（2）便可精確求得地震波走時(shí)，從而規(guī)避震源奇異性問(wèn)題。本文選擇T0為震源(x0，z0)與待求點(diǎn)（x，z）之間的距離函數(shù)，即

則對(duì)于震源點(diǎn)(x0，z0)有

對(duì)于任意一點(diǎn)，有

則對(duì)于震源點(diǎn)(x0，z0)，由因式分解程函方程式（3），可得

上述T0(x0，z0)與T1(x0，z0)即為本算法的初始化條件。

2 數(shù)值求解T1

傳統(tǒng)FMM 采用有限差分近似走時(shí)梯度，利用窄帶模擬波前擴(kuò)展，采用堆排序算法尋找波前中走時(shí)最小點(diǎn)作為次級(jí)震源。本文算法同樣基于以上技術(shù)，這些技術(shù)已在許多FMM 算法相關(guān)文獻(xiàn)［21，31-32］中有詳細(xì)闡述，此處不再贅述。本文算法的關(guān)鍵在于如何利用FMM 求解因式分解程函方程，下面分別闡述兩種方法求解T1。

2.1 規(guī)則網(wǎng)格求解T1

首先，對(duì)因式分解程函方程進(jìn)行離散，則有（為表述方便，用點(diǎn)（x，z）的編號(hào)（i，j）表示其位置）

式中：max（A，B，C）表示求取A、B、C中最大值；和分別表示空間點(diǎn)（i，j）處向l負(fù)方向和l正方向的n階差分算子。數(shù)值計(jì)算時(shí)分別采用一階或二階差分進(jìn)行計(jì)算。嚴(yán)格地講，二階差分屬于混合差分，因?yàn)楫?dāng)二階差分條件不滿足時(shí)，即采用一階差分。以下以x方向?yàn)槔榻B一階和二階差分求解因式分解程函方程過(guò)程。

規(guī)則網(wǎng)格（圖1）中，假定此時(shí)（i，j）點(diǎn)的走時(shí)為待求點(diǎn)走時(shí)，引入符號(hào)Tknow表示T值已知、Tunknow表示T值未知。對(duì)于一階差分，可得到式（9）； 根據(jù)鏈?zhǔn)角髮?dǎo)法，則有式（10）

圖1 規(guī)則網(wǎng)格走時(shí)擾動(dòng)因子T1求解示意圖

式中：T1(i-1，j)和T1(i+1，j)分別為(i-1，j)點(diǎn)和(i+1，j)點(diǎn)的走時(shí)擾動(dòng)因子，可結(jié)合T(i-1，j)和T(i+1，j)由式（2）求出；T0(i，j)可由式（4）解析求出；h為網(wǎng)格剖分間距；。這樣，式（10）中就只有一個(gè)待求未知量T1(i，j)。

對(duì)于二階差分，實(shí)際為混合階數(shù)差分（當(dāng)不滿足二階差分條件時(shí)，即運(yùn)用一階差分）。與一階差分相似，仍需先判斷待求點(diǎn)(i，j)鄰近點(diǎn)狀態(tài)。同樣，以x方向?yàn)槔?，具體如下

式中各變量物理意義與一階差分（式（9））類(lèi)似。

z方向差分算子的推導(dǎo)，與x方向類(lèi)似。至此，再結(jié)合窄帶技術(shù)［21］，便可實(shí)現(xiàn)直達(dá)波走時(shí)計(jì)算。

2.2 不規(guī)則網(wǎng)格求解T1

為提高走時(shí)反演成像的分辨率，采用多震相走時(shí)聯(lián)合反演是一個(gè)不錯(cuò)的選擇，這就要求正演算法具有計(jì)算多震相地震波走時(shí)的能力。分區(qū)多步計(jì)算技術(shù)自提出以來(lái)，就與多種網(wǎng)格類(lèi)地震波走時(shí)計(jì)算方法相結(jié)合，用于模擬多震相地震波傳播。本算法同樣作為網(wǎng)格單元類(lèi)算法，也可結(jié)合分區(qū)多步計(jì)算技術(shù)實(shí)現(xiàn)多震相地震波走時(shí)的計(jì)算。

然而實(shí)際地表（或地下界面）多呈不規(guī)則性，若用規(guī)則網(wǎng)格剖分，則地下界面（或地表）附近會(huì)產(chǎn)生不規(guī)則網(wǎng)格（圖2b 中M1M2界面）。這樣，上述規(guī)則網(wǎng)格的迎風(fēng)差分公式顯然已不適用，因此需要針對(duì)界面附近不規(guī)則網(wǎng)格，推導(dǎo)建立新的局部走時(shí)計(jì)算公式。圖2不規(guī)則網(wǎng)格中的點(diǎn)可劃分為兩類(lèi)：一類(lèi)是界面附近點(diǎn)（C點(diǎn)）； 另一類(lèi)是界面點(diǎn)（M點(diǎn)）。下面分別針對(duì)這兩類(lèi)點(diǎn)建立局部走時(shí)計(jì)算公式。

圖2 模型界面及網(wǎng)格剖分示意圖（a）及其局部放大圖（b）

（1）C點(diǎn)為當(dāng)前待計(jì)算點(diǎn)

此時(shí)以C點(diǎn)為中心的z方向出現(xiàn)CM≠h，即在C點(diǎn)的x軸方向與z軸上半軸方向仍可采用常規(guī)迎風(fēng)差分格式，而z軸下半軸常規(guī)迎風(fēng)差分格式顯然已不能適用，為此建立C點(diǎn)z軸向下方向的不等距一階、二階迎風(fēng)差分公式（記|CM|=h1）

而二階差分算子的推導(dǎo)過(guò)程如下，利用泰勒函數(shù)將T1(M)與T1(D)在C點(diǎn)附近展開(kāi)

將式（15）代入式（14），可得不等距二階迎風(fēng)差分公式。

此外，由于|CM|≠BC，導(dǎo)致T(M)與T(B)的大小已無(wú)可比性，則規(guī)則網(wǎng)格的迎風(fēng)差分條件在此已不適用，因此建立如下不等距迎風(fēng)差分格式

（2）M點(diǎn)為當(dāng)前待計(jì)算點(diǎn)

此時(shí)M點(diǎn)周?chē)乃膫€(gè)網(wǎng)格均為非正交，導(dǎo)致程函方程不成立。Rawlinsnon 等［31］采用三角網(wǎng)格縫合界面，而后在三角網(wǎng)格中求解程函方程，但實(shí)現(xiàn)較為困難。本文綜合不等距迎風(fēng)差分公式、費(fèi)馬原理與惠更斯原理，求解界面走時(shí)。對(duì)于M點(diǎn)的z軸下方向，同樣有不等距一階、二階差分公式

式中h2為M、D之間的距離。對(duì)于M點(diǎn)z軸上方向，表達(dá)式類(lèi)似。

此外對(duì)于M點(diǎn)x方向，由于界面起伏，x方向常規(guī)差分公式已不再適用。本文采用惠更斯原理，即將M1、M2視為子波源，就有

式中：li為Mi與M之間的距離；MiM為Mi與M之間的波慢度。

最后，再利用費(fèi)馬原理，得到

綜上所述，由于起伏界面的存在，可將界面上的點(diǎn)分為三類(lèi)：第一類(lèi)為規(guī)則網(wǎng)格中節(jié)點(diǎn)，采用常規(guī)迎風(fēng)差分公式（式（8）～式（10））計(jì)算修正因子T1，然后根據(jù)式（2）計(jì)算該點(diǎn)地震波走時(shí)； 第二類(lèi)為界面附近節(jié)點(diǎn)（C），對(duì)于該點(diǎn)遠(yuǎn)離界面方向和x方向，走時(shí)擾動(dòng)因子的計(jì)算方法與規(guī)則網(wǎng)格相同，對(duì)于該點(diǎn)的界面方向，采用式（13）計(jì)算該方向的走時(shí)擾動(dòng)因子，最后根據(jù)費(fèi)馬原理，由式（16）計(jì)算該點(diǎn)走時(shí)；第三類(lèi)為界面上點(diǎn)（M），對(duì)于該點(diǎn)z軸方向，建立不等距迎風(fēng)差分公式（式（17））計(jì)算該點(diǎn)在z方向的走時(shí)擾動(dòng)因子，最后結(jié)合費(fèi)馬原理、惠更斯原理，由式（19）計(jì)算M點(diǎn)地震波走時(shí)。

3 多震相地震波走時(shí)計(jì)算方法

基于前面構(gòu)建的起伏地表模型下規(guī)則網(wǎng)格與不規(guī)則網(wǎng)格的迎風(fēng)差分格式，再結(jié)合分區(qū)多步技術(shù)與窄帶技術(shù)，便可實(shí)現(xiàn)多震相地震波走時(shí)計(jì)算。所謂分區(qū)多步計(jì)算技術(shù)即指按照地下起伏界面的具體情況，將模型分成獨(dú)立的層狀（或起伏層狀）區(qū)域，相鄰區(qū)域由地下界面相連接［18］（圖2a）。窄帶技術(shù)根據(jù)模型中網(wǎng)格節(jié)點(diǎn)所處狀態(tài)將所有節(jié)點(diǎn)分為三種：“完成點(diǎn)”、“窄帶點(diǎn)”及“遠(yuǎn)離點(diǎn)”。其中“完成點(diǎn)”表示走時(shí)計(jì)算已結(jié)束的網(wǎng)格點(diǎn)；“窄帶點(diǎn)”表示當(dāng)前波前上的網(wǎng)格節(jié)點(diǎn)，其走時(shí)大小還可更新；“遠(yuǎn)離點(diǎn)”表示走時(shí)計(jì)算還未實(shí)施的網(wǎng)格節(jié)點(diǎn)。

表1給出了本文所提FMM-F 算法計(jì)算直達(dá)波走時(shí)的實(shí)現(xiàn)框圖，當(dāng)要計(jì)算其他震相的地震波走時(shí)時(shí)，只需將界面上走時(shí)最小的點(diǎn)作為新的子震源，而后與計(jì)算直達(dá)波相同，計(jì)算當(dāng)前所需震相的地震波走時(shí)即可。

表1 FMM-F 算法實(shí)施過(guò)程

現(xiàn)以圖2a所示模型計(jì)算上界面反射波為例，說(shuō)明反射波走時(shí)計(jì)算原理（震源所在位置為第一分區(qū)），其他震相的走時(shí)計(jì)算也可類(lèi)推。

（1）首先設(shè)定算法初始化條件，給定T0(x0，z0)=0，T1(x0，z0)=S(x0，z0)，并根據(jù)窄帶技術(shù)將震源點(diǎn)記為“完成點(diǎn)”，即該點(diǎn)地震波走時(shí)不再更新；

（2）判斷當(dāng)前震源點(diǎn)所在網(wǎng)格是否包含界面或起伏地表：若是，則采用非規(guī)則網(wǎng)格迎風(fēng)差分格式（式（13）～式（19））計(jì)算震源周?chē)c(diǎn)走時(shí)擾動(dòng)因子T1；若否，則采用規(guī)則網(wǎng)格迎風(fēng)差分格式（式（8）～式（12））計(jì)算T1，通過(guò)式（2）和式（4）得到震源周?chē)c(diǎn)地震波走時(shí)，這些點(diǎn)構(gòu)成的面即為當(dāng)前波前所在位置，其狀態(tài)記為“窄帶點(diǎn)”；

（3）在“窄帶點(diǎn)”中搜尋走時(shí)最小點(diǎn)，將此點(diǎn)作為次級(jí)震源并將其狀態(tài)更新為“完成點(diǎn)”繼續(xù)循環(huán)第（2）步計(jì)算，直到得到震源所在分區(qū)的全局走時(shí)；

（4）讀取界面離散點(diǎn)地震波走時(shí)并將其狀態(tài)均賦為“窄帶點(diǎn)”，返回第（3）步，要注意的是由于是計(jì)算反射波，所以在第（2）步時(shí)只需考慮界面以上的網(wǎng)格節(jié)點(diǎn)即可。

若要追蹤透射波，則在第（4）步返回第（3）步時(shí)，自界面上走時(shí)最小的點(diǎn)向透射區(qū)域進(jìn)行波前擴(kuò)展即可； 若要追蹤轉(zhuǎn)換波，則自離散界面上走時(shí)最小的點(diǎn)開(kāi)始向轉(zhuǎn)換波所在區(qū)域進(jìn)行波前擴(kuò)展掃描時(shí)時(shí)，采用不同的速度模型即可。簡(jiǎn)而言之，分區(qū)多步計(jì)算的原理就是根據(jù)所要計(jì)算地震波的類(lèi)型，設(shè)計(jì)要調(diào)用的獨(dú)立區(qū)域的先后順序及對(duì)應(yīng)的速度模型，達(dá)到運(yùn)用簡(jiǎn)單組合完成復(fù)雜模型中多震相地震波走時(shí)計(jì)算之目的。

4 數(shù)值模擬及結(jié)果分析

采用幾組模型進(jìn)行數(shù)值模型計(jì)算，通過(guò)本算法與Multistage FMM 算法［31］詳盡（主要是計(jì)算精度與計(jì)算效率）對(duì)比，驗(yàn)證本文算法的優(yōu)越性。當(dāng)選用模型存在解析解時(shí)，將以解析解作為參考解； 當(dāng)選用模型不存在解析解時(shí)，將以發(fā)展較成熟的不規(guī)則最短路徑算法（ISMP）［37］作為參考解，且通過(guò)加密ISMP算法網(wǎng)格間距，確保該算法計(jì)算結(jié)果的精度。計(jì)算所用處理器配置為： Intel（R） Core（TM） i7-10750H CPU @2.60 GHz。另外，在分析算法所用CPU 時(shí)間時(shí)，均是設(shè)置程序循環(huán)執(zhí)行100 次后取均值得到每次所用CPU 時(shí)間，以避免偶然性。

4.1 模型1（均勻模型）

2D 均勻模型大小為4 km×4 km，地震波傳播速度為1 km/s，震源位于（2 km，2 km）；3D 均勻模型大小為10 km×10 km×10 km，地震波傳播速度為2 km/s，震源位于（5 km，5km，5 km）。對(duì)于FMM，采用震源網(wǎng)格加密方式保證精度，加密區(qū)域范圍為[x0-1 km，x0+1 km]、[z0-1 km，z0+1 km]（2D）或[x0-1 km，x0+1 km]、[y0-1 km，y0+1 km]、[z0-1 km，z0+1 km（]3D），(x0，y0，z0）（為震源位置坐標(biāo)，加密網(wǎng)格間距為0.005 km。本文選用平均絕對(duì)誤差與平均相對(duì)誤差兩個(gè)指標(biāo)評(píng)估計(jì)算精度，平均絕對(duì)誤差計(jì)算公式為：；同時(shí)，平均相對(duì)誤差的計(jì)算公式為：。其中Texact為解析解，Tcal為通過(guò)FMM 或FMM-F 所計(jì)算的結(jié)果，N為模型中網(wǎng)格節(jié)點(diǎn)總數(shù)。

表2（2D）和表3（3D）分別給出均勻模型不同網(wǎng)格間距剖分下FMM 與FMM-F 的計(jì)算結(jié)果對(duì)比。從該表可見(jiàn)，對(duì)于均勻模型，本文所提FMM-F計(jì)算結(jié)果誤差為零（注意：誤差分析結(jié)果已消除因計(jì)算機(jī)自身小數(shù)保留位數(shù)而導(dǎo)致的計(jì)算誤差）； 同時(shí)，可看到隨著網(wǎng)格間距的減小，兩種算法的計(jì)算效率都降低了，但在相同的網(wǎng)格剖分間距下，F(xiàn)MM-F 計(jì)算效率明顯高于FMM，考慮是由于FMM 采用震源網(wǎng)格加密方式保證精度，從而使得FMM 計(jì)算效率降低。另外，當(dāng)采用二階差分時(shí)，兩種算法計(jì)算效率都降低了，但降低并不明顯。

表2 2D 均勻模型中FMM 與FMM-F 計(jì)算地震波走時(shí)誤差及CPU 時(shí)間對(duì)比

表3 3D 均勻模型中FMM 與FMM-F 計(jì)算地震波走時(shí)誤差及CPU 時(shí)間對(duì)比

4.2 模型2（速度線性增加模型）

再采用速度線性增加模型（圖3）進(jìn)行試算。2D模型（圖3a）大小為8 km×4 km，速度自地表4 km/s線性增至底界面的6 km/s，震源位于(4 km，0)； 3D模型（圖3b）大小為10 km×10 km×10 km，速度自地表2 km/s 線性增至底界面的4 km/s，震源位于（5 km，5 km，0）。該模型的地震波走時(shí)解析解由下式求出［37］

圖3 速度線性增加的2D（a）和3D（b）模型

式中：α為速度增加梯度；v0為震源處地震波傳播速度；v(x)為x點(diǎn)的速度；x0為震源位置坐標(biāo)。

表4 和表5 分別給出不同網(wǎng)格間距剖分計(jì)算2D和3D 線性模型時(shí)FMM 與FMM-F 的計(jì)算結(jié)果對(duì)比。同樣采用震源網(wǎng)格加密方式保證FMM 計(jì)算精度。水平軸加密區(qū)域范圍為[x0-1 km，x0+1 km（]2D）或[x0-1 km，x0+1 km]、[y0-1 km，y0+1 km（]3D），由于此模型震源置于地表，所以深度方向加密區(qū)域范圍為[z0，z0+1 km]，（x0，y0，z0）為震源位置坐標(biāo)，加密網(wǎng)格間距為5 m。同樣，選用平均絕對(duì)誤差與平均相對(duì)誤差兩個(gè)指標(biāo)評(píng)估計(jì)算精度。

表4 圖3a 模型中FMM 與FMM-F 計(jì)算地震波走時(shí)誤差及CPU 時(shí)間對(duì)比

表5 圖3b 模型中FMM 與FMM-F 計(jì)算地震波走時(shí)誤差及CPU 時(shí)間對(duì)比

從表4和表5可見(jiàn)，隨著網(wǎng)格間距的不斷減小，不論是FMM 還是FMM-F，其平均絕對(duì)誤差與平均相對(duì)誤差都不斷減小，計(jì)算精度不斷提高。但同時(shí)，網(wǎng)格間距的減小，導(dǎo)致網(wǎng)格節(jié)點(diǎn)數(shù)變多，計(jì)算效率也隨之下降。另外還發(fā)現(xiàn)，當(dāng)采用二階差分時(shí)，兩種算法的計(jì)算精度都有大幅度提高，且與均勻模型相同，計(jì)算效率降低并不明顯，因此建議實(shí)際應(yīng)用中采用二階差分格式。重要的是在兩種算法對(duì)比中發(fā)現(xiàn)，網(wǎng)格間距相同時(shí)，2D 或3D 情形下，無(wú)論是計(jì)算精度還是CPU 計(jì)算時(shí)間，F(xiàn)MM-F 優(yōu)于FMM。這是因?yàn)橛?jì)算精度方面，F(xiàn)MM-F 成功地避免了震源點(diǎn)的奇異性問(wèn)題； 而計(jì)算效率方面則是傳統(tǒng)FMM 為了克服震源點(diǎn)的奇異性需對(duì)震源附近網(wǎng)格細(xì)化，導(dǎo)致計(jì)算量增加。

為進(jìn)一步分析誤差分布情況，又給出當(dāng)h=0.05 m 時(shí)兩種算法分別采用一階（圖4（2D）、圖6（3D））、二階（圖5 （2D）、圖7（3D））差分所得的絕對(duì)誤差和相對(duì)誤差分布圖。可發(fā)現(xiàn)傳統(tǒng)FMM 計(jì)算所得地震波走時(shí)在震源對(duì)角線附近區(qū)域誤差較大，當(dāng)采用一階差分格式時(shí)，其最大絕對(duì)誤差分別可達(dá)4.80 ms（2D）和20 ms （3D），最大相對(duì)誤差則分別高達(dá)7%（2D）和4% （3D）； 而FMM-F 最大絕對(duì)誤差僅分別為0.16 ms （2D）和1.5 ms （3D），而最大相對(duì)誤差也分別僅為0.016%（2D）和0.04%（3D）。當(dāng)采用二階差分格式時(shí)，兩種算法的計(jì)算精度都得到了很大提高，F(xiàn)MM 最大絕對(duì)誤差分別降至0.4 ms（2D）和1.5 ms（3D），同樣最大相對(duì)誤差也分別減小到2.7%（2D 和2%（3D）； 而FMM-F 最大絕對(duì)誤差則分別進(jìn)一步減至0.04 ms（2D）和0.1 ms（3D），最大相對(duì)誤差則分別減至0.009%（2D）和0.02%（3D）。可見(jiàn)，本文的FMM-F 不論采用一階還是二階差分格式，其計(jì)算精度高于傳統(tǒng)FMM； 同時(shí)，通過(guò)誤差分布圖可發(fā)現(xiàn)FMM-F完全解決了震源點(diǎn)的奇異性問(wèn)題。

圖4 2D 模型FMM（a、b）、FMM-F（c、d）算法一階差分計(jì)算的絕對(duì)誤差（a、c）和相對(duì)誤差（b、d）

圖5 對(duì)應(yīng)圖4 的二階差分計(jì)算結(jié)果

圖6 3D 模型中FMM（上）和FMM-F（下）算法采用一階差分格式計(jì)算誤差分布

圖7 對(duì)應(yīng)圖6 采用二階差分格式計(jì)算結(jié)果

4.3 模型3（起伏界面模型）

在上文分析了直達(dá)波的走時(shí)計(jì)算精度的基礎(chǔ)上，通過(guò)一個(gè)數(shù)值例子驗(yàn)證本文算法的多震相走時(shí)計(jì)算精度。因沒(méi)有解析解，故將Multistage ISPM 算法［39］的計(jì)算結(jié)果作為參考值。為了保證Multistage ISPM算法較高的計(jì)算精度，網(wǎng)格主節(jié)點(diǎn)間距設(shè)置為0.25 km，且在主節(jié)點(diǎn)間加入19 個(gè)次級(jí)節(jié)點(diǎn)，這樣通過(guò)Multistage ISPM 算法求取結(jié)果的精度較高，可作為參考值。所用模型如圖8a （2D）、圖9a （3D）所示，其大小分別為100 km×50 km （2D） 和100 km×100 km×50 km （3D）。

圖8 2D 起伏模型計(jì)算精度分析

圖9 3D起伏模型及不同剖面的速度分布

2D 模型地表的起伏函數(shù)設(shè)置為：z=2+2×sin (2πx/50 )，且在25～35 km 深度范圍內(nèi)含有一個(gè)傾斜斷層，斷層大約位于x=50km（圖8a 中黑線），震源置于點(diǎn)（50 km，15 km）處，并在地面以2 km 間距均勻布放51個(gè)檢波器。模型速度由地表的2 km/s線性增至斷層界面的4 km/s，而后又線性增至底界面的6 km/s（圖8a）。

由于常規(guī)Multistage FMM［31］采用三角網(wǎng)格離散界面，為滿足迎風(fēng)差分條件，需對(duì)鈍角三角網(wǎng)格進(jìn)行剖分，導(dǎo)致常規(guī)Multistage FMM 對(duì)界面的離散與Multistage ISPM 算法和本文算法不一致，且無(wú)可比性，故此部分不再將Multistage FMM 作為對(duì)比的算法，而僅將Multistage ISPM 算法作為參照，以驗(yàn)證本文算法的有效性。

圖8c、圖8d 分別為采用圖8a 的炮檢排列與速度結(jié)構(gòu)模型時(shí)，自震源激發(fā)后的直達(dá)波和自底界面反射上行波的走時(shí)等值線圖，可見(jiàn)Multistage ISPM 算法計(jì)算結(jié)果對(duì)應(yīng)的黑色虛線、本文算法采用一階差分計(jì)算結(jié)果對(duì)應(yīng)的紅色虛線及采用二階差分對(duì)應(yīng)的綠色虛線三種等值線的重合度較高，表明本文算法在計(jì)算直達(dá)波與反射波時(shí)精度較高。為進(jìn)一步突顯對(duì)比效果，還對(duì)等值線圖做了局部放大顯示（圖8c、圖8d右下角）及與地面檢波器觀測(cè)反射波到時(shí)的相對(duì)誤差（圖8b，計(jì)算公式為，通過(guò)局部放大圖清晰可見(jiàn)采用二階差分時(shí)的等值線與Multistage ISPM 算法的重合度更高，表明采用二階差分的計(jì)算精度更高，通過(guò)與地面檢波器的相對(duì)誤差對(duì)比，可得出同樣結(jié)論。

針對(duì)3D 模型，模擬計(jì)算了8 種震相地震波的走時(shí)，并以Multistage ISPM 算法作為標(biāo)準(zhǔn)計(jì)算了相對(duì)誤差（誤差計(jì)算公式與2D 相同，圖10）。它們分別為經(jīng)上界面反射后的純反射波P1P （圖10a）、S1S（圖10b）與反射轉(zhuǎn)換波P1S （圖10c）、S1P （圖10d），經(jīng)底界面反射后的純反射波P2P（圖10e）與反射轉(zhuǎn)換波S2P （圖10f），以及經(jīng)上界面與下界面多次透射、反射的轉(zhuǎn)換波P2S1P2S、S2S1P2P。對(duì)比圖10上圖與下圖可見(jiàn)，經(jīng)上界面反射后的地震波走時(shí)的相對(duì)誤差明顯大于經(jīng)底界面反射的地震波走時(shí)，這是緣于上界面存在一個(gè)較大孔洞，而底界面僅為一平面； 同時(shí)，8 種地震波震相的走時(shí)相對(duì)誤差均在0.4%以內(nèi)，可見(jiàn)本文算法對(duì)復(fù)雜起伏地表（含地下孔穴等劇烈起伏地質(zhì)異常體）模型中多震相地震波走時(shí)的計(jì)算精度也較高，可作為一種高效、計(jì)算精度高的多震相地震波走時(shí)場(chǎng)計(jì)算的有效方法。

圖10 3D 起伏模型地表檢波器接收的多震相走時(shí)的相對(duì)誤差分布

5 結(jié)論

由于原始程函方程的FMM 采用平面波近似波前，其波前曲率較大，導(dǎo)致震源奇異性問(wèn)題，且利用三角形（2D）或四面體（3D）縫合界面，實(shí)現(xiàn)較為困難。為此，本文發(fā)展了一種利用快速行進(jìn)法求解因式分解程函方程的多震相地震波走時(shí)計(jì)算算法。該方法采用因式分解求解程函方程解決震源奇異性問(wèn)題，以不等距差分格式精確定位地表位置，通過(guò)迎風(fēng)混合網(wǎng)格差分公式計(jì)算地震波走時(shí)，并用窄帶技術(shù)作為算法的整體波前擴(kuò)展技術(shù)。數(shù)值模擬實(shí)例表明，以解析解、Multistage ISPM 算法［31］作為參考解，與Multistage FMM 算法［18］對(duì)比分析，得到以下認(rèn)識(shí)：

（1）所提算法計(jì)算精度較高，不論是平均絕對(duì)誤差還是平均相對(duì)誤差都顯著降低；

（2）通過(guò)棄用震源網(wǎng)格細(xì)化以保證計(jì)算精度，且該算法計(jì)算效率較高；

（3）利用該算法可靈活、穩(wěn)定地計(jì)算復(fù)雜地表?xiàng)l件下的走時(shí)。