何雷宇 嚴(yán) 星 白超英*③

(①長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院地球物理系，陜西西安 710054； ②新疆財(cái)經(jīng)大學(xué)網(wǎng)絡(luò)與實(shí)驗(yàn)教學(xué)中心，新疆烏魯木齊 830011； ③長(zhǎng)安大學(xué)計(jì)算地球物理研究所，陜西西安 710054)

1 引言

層析成像始于醫(yī)學(xué)中的CT技術(shù)，Aki等[1]首次將其應(yīng)用于地震學(xué)研究中。經(jīng)過(guò)四十多年的發(fā)展和完善，地震層析成像已成為研究地下地質(zhì)結(jié)構(gòu)行之有效的主要技術(shù)途徑之一，并且涌現(xiàn)出多種成像的方法技術(shù)[2]。

傳統(tǒng)射線類走時(shí)成像方法一般采用規(guī)則網(wǎng)格單元(2D中的矩形網(wǎng)格單元、3D中的立方體網(wǎng)格單元)模型參數(shù)化方式，這種簡(jiǎn)單模型參數(shù)化可以很容易地通過(guò)網(wǎng)格節(jié)點(diǎn)的編號(hào)確定該節(jié)點(diǎn)的空間坐標(biāo)，從而節(jié)約計(jì)算機(jī)內(nèi)存空間。然而這種簡(jiǎn)單的模型參數(shù)化形式也存在諸多問(wèn)題。例如：若要對(duì)復(fù)雜區(qū)域進(jìn)行較為精確的描述，則必須進(jìn)行網(wǎng)格單元的精細(xì)剖分，從而增加計(jì)算時(shí)間。更為重要的是規(guī)則網(wǎng)格單元參數(shù)化對(duì)起伏地表及地下不規(guī)則波阻抗界面無(wú)法進(jìn)行準(zhǔn)確的擬合，從而造成正演中不必要的計(jì)算誤差，進(jìn)而導(dǎo)致反演中出現(xiàn)假象。

另一方面，從反演角度來(lái)說(shuō)，規(guī)則網(wǎng)格的精細(xì)剖分，意味著更多的未知模型參數(shù)參與反演，從理論上就需要更多的射線進(jìn)行覆蓋，故而求解的方程組的維數(shù)也將大大增加，數(shù)據(jù)對(duì)模型的約束也就越差。此外，更多的未知參數(shù)意味著參數(shù)矢量中欠定的和位于零空間的分量數(shù)目更多，性態(tài)大大變壞； 層析反演難度的增加，通常需加入正則化約束以改變方程組的病態(tài)，而這又增加了存儲(chǔ)量和計(jì)算量[3]。

上述問(wèn)題的有效解決方案之一是采用不規(guī)則網(wǎng)格單元(包括三角網(wǎng)格或四面體網(wǎng)格)進(jìn)行模型參數(shù)化。例如可根據(jù)模型中不同區(qū)域物性的復(fù)雜程度選擇不規(guī)則網(wǎng)格單元的大小，同時(shí)不規(guī)則網(wǎng)格單元剖分對(duì)起伏地表和地下波阻抗界面(含不規(guī)則異常體)的描述更加精確，從而確保了正演所需的計(jì)算精度。反演中由于未知模型參數(shù)相對(duì)較少，則所需射線也相對(duì)較少，故而待求解的方程組維數(shù)也大大減小，從總體上減少了反演的存儲(chǔ)量[3]。

三維情況下使用四面體網(wǎng)格單元的另一個(gè)優(yōu)勢(shì)在于：實(shí)際中地震臺(tái)站的分布是不均勻的，如北京圈附近的臺(tái)站密度與我國(guó)西部地區(qū)的臺(tái)站密度、陸地的臺(tái)站密度與海洋區(qū)域的臺(tái)站密度是不同的，加之天然地震的發(fā)生本身就具有不均勻性，兩者結(jié)合造成了地震數(shù)據(jù)覆蓋密度的極度差異性，所以有必要對(duì)數(shù)據(jù)密度大的區(qū)域進(jìn)行更細(xì)致的剖分以增加數(shù)據(jù)的利用效率。

不規(guī)則網(wǎng)格單元在地震數(shù)值模擬中有較為廣泛的應(yīng)用。正演方面，如在波場(chǎng)模擬中使用的不規(guī)則網(wǎng)格[4，5]、變網(wǎng)格[6]和貼體網(wǎng)格[7]等。而在地震射線追蹤研究方面，李波濤等[8]將三角剖分與波前重建法相結(jié)合，實(shí)現(xiàn)了初至波射線路徑和走時(shí)的快速計(jì)算； Yu等[9]引入單元密度控制函數(shù)，實(shí)現(xiàn)了三角網(wǎng)格下初至波射線追蹤的全局算法； Lan等[10]借助于流體力學(xué)中的貼體網(wǎng)格和坐標(biāo)變換，引入Lax-Friedrichs快速掃描法求解與地形有關(guān)的程函方程以獲取地震波初至走時(shí)場(chǎng)。

上述幾位學(xué)者的研究主要是解決了2D模型中不規(guī)則網(wǎng)格單元模型參數(shù)化下初至波走時(shí)的計(jì)算與射線追蹤問(wèn)題。Vinje等[11]在開(kāi)放模型的波前構(gòu)造法中將三角網(wǎng)格應(yīng)用于反射界面的描述和波前面的描述，并實(shí)現(xiàn)了反射波射線追蹤。趙瑞等[12]將梯形網(wǎng)格與矩形網(wǎng)格結(jié)合，實(shí)現(xiàn)了多震相射線追蹤； Bai等[13]提出一種基于三角(2D)或四面體(3D)單元剖分下的最短路徑射線追蹤算法(Triangular Shortest-path Method，TSPM)，并實(shí)現(xiàn)了多震相地震射線路徑的追蹤與走時(shí)計(jì)算； 李曉玲等[14]在趙瑞等[12]研究的基礎(chǔ)上，將混合網(wǎng)格下的射線追蹤算法推廣至各向異性TI介質(zhì)中。但混合網(wǎng)格僅可模擬起伏地表，無(wú)法對(duì)模型內(nèi)部進(jìn)行差異化剖分，優(yōu)勢(shì)有限； 隨后，李興旺等[15]實(shí)現(xiàn)了四面體網(wǎng)格參數(shù)化下TI介質(zhì)中的多震相射線追蹤。

基于不規(guī)則網(wǎng)格單元剖分下的層析成像研究最早可追溯到20世紀(jì)80年代。Tarantola等[16]意識(shí)到在地震層析成像中使用固定的均勻網(wǎng)格單元的一些局限性并提出了block-less模型參數(shù)化方法； Sambridge等[17]率先將Delaunay四面體和Voronoi多面體應(yīng)用于地震層析成像領(lǐng)域； Gudmundsson等[18]利用Delaunay四面體和Voronoi多面體進(jìn)行模型參數(shù)化，使用初至P波層析成像改進(jìn)了區(qū)域上地幔地球參考模型； Curtis等[19]通過(guò)對(duì)給定的一組點(diǎn)采用Delaunay三角化來(lái)描述模型的速度結(jié)構(gòu)，進(jìn)而研究了2D情況下的井間反演問(wèn)題； B?hm等[20]用Delaunay三角形和Voronoi多面體在3D反射層析成像中發(fā)展了一種自適應(yīng)反演方法； Sambridge等[21]使用自適應(yīng)Delaunay四面體網(wǎng)格單元進(jìn)行了全球P波層析成像研究。

中國(guó)國(guó)內(nèi)對(duì)不規(guī)則網(wǎng)格層析成像的研究較少。成谷等[3]對(duì)比了三角網(wǎng)格層析成像與矩形網(wǎng)格層析成像的優(yōu)缺點(diǎn)； 于師建等[22]采用三角網(wǎng)格模型參數(shù)化進(jìn)行走時(shí)層析成像研究。他們所做的研究主要是不規(guī)則網(wǎng)格參數(shù)化下初至波走時(shí)層析成像問(wèn)題，但并未引進(jìn)反射波。鑒于此，本文研究了2D三角網(wǎng)格模型參數(shù)化的多震相聯(lián)合反演及速度與反射界面同時(shí)反演問(wèn)題[23]。為了更好地解決實(shí)際地震工作的需要，還將該算法進(jìn)一步推廣至3D四面體網(wǎng)格。在此過(guò)程中，要解決的主要問(wèn)題是如何進(jìn)行模型參數(shù)化、用參數(shù)化后的網(wǎng)格進(jìn)行射線追蹤及如何求解Jacob矩陣元素并進(jìn)行合理分配。顯然，3D情況下算法實(shí)現(xiàn)的難度總體上大于2D情況。

文中基于Multistage TSPM正演算法[13]，結(jié)合共軛梯度求解帶約束的阻尼最小二乘解反演算法，同時(shí)得到了四面體網(wǎng)格單元下走時(shí)關(guān)于速度、走時(shí)關(guān)于反射界面的偏導(dǎo)計(jì)算公式； 討論并實(shí)現(xiàn)了四面體網(wǎng)格單元模型參數(shù)化下3D復(fù)雜模型中地震多震相走時(shí)同時(shí)反演成像技術(shù)。本文主要研究?jī)?nèi)容包括： ①反射界面已知時(shí)地震多震相走時(shí)聯(lián)合成像技術(shù)； ②多震相走時(shí)聯(lián)合同時(shí)反演速度模型和反射界面的技術(shù)。數(shù)值模擬實(shí)驗(yàn)結(jié)果表明，四面體網(wǎng)格單元模型參數(shù)化下的多波走時(shí)成像及同時(shí)反演成像結(jié)果較為可信，可作為起伏復(fù)雜介質(zhì)中速度場(chǎng)和反射界面重建的一種有效方法。

2 多震相射線追蹤算法

2.1 模型參數(shù)化

本文采用的網(wǎng)格由開(kāi)源的Tetgen生成[24]。利用Multistage TSPM進(jìn)行射線追蹤之前，首先要對(duì)模型進(jìn)行四面體網(wǎng)格單元參數(shù)化(圖1)。采用四面體單元進(jìn)行模型參數(shù)化時(shí)，定義四面體單元的四個(gè)角點(diǎn)為主節(jié)點(diǎn)。為了增加射線出射角的覆蓋率(即保證正演的計(jì)算精度)，在四面體單元的四條邊上等間距插入次級(jí)節(jié)點(diǎn)，四條邊插值結(jié)束后再對(duì)四個(gè)面進(jìn)行插值。四面體網(wǎng)格單元內(nèi)部沒(méi)有節(jié)點(diǎn)，但炮點(diǎn)和檢波點(diǎn)可置于其中(圖1)。模型中速度的采樣是在主節(jié)點(diǎn)上進(jìn)行，而次級(jí)節(jié)點(diǎn)上的速度則是通過(guò)速度插值得到(詳見(jiàn)下文)。

圖1 速度插值示意圖

2.2 最小走時(shí)的計(jì)算和速度插值

從炮點(diǎn)開(kāi)始計(jì)算，節(jié)點(diǎn)走時(shí)的計(jì)算公式為(兩點(diǎn)間距離除以平均速度)

(1)

式中：i是波前點(diǎn)中當(dāng)前的次級(jí)震源；j為待計(jì)算任意節(jié)點(diǎn)；Nj為待計(jì)算節(jié)點(diǎn)的集合；ti為炮點(diǎn)到達(dá)i節(jié)點(diǎn)射線的最小走時(shí)；D(xi,xj)是節(jié)點(diǎn)i與j之間的距離；v(xi)、v(xj)分別是i和j節(jié)點(diǎn)上的速度。當(dāng)i節(jié)點(diǎn)不是主節(jié)點(diǎn)時(shí)，可通過(guò)下式求取該節(jié)點(diǎn)速度

(2)

式中：vj表示所在單元各主節(jié)點(diǎn)的速度，在三維模型中k=4；uj為點(diǎn)在四面體中的體積坐標(biāo)，例如圖1中P點(diǎn)在四面體T1T2T3T4內(nèi)的體積坐標(biāo)uj可表示為

(3)

式中V(PT2T3T4)、V(T1T2T3T4)等表示各四面體的體積。顯然，四點(diǎn)共面時(shí)體積為零。

2.3 分區(qū)多步計(jì)算原理

完成模型參數(shù)化和速度插值之后，為了減少計(jì)算量，避免對(duì)射線未經(jīng)過(guò)區(qū)域的計(jì)算，可根據(jù)反射界面的起伏形態(tài)將模型分成若干(層狀)區(qū)域。從震源所在區(qū)域開(kāi)始計(jì)算，根據(jù)給定的震相種類對(duì)特定射線進(jìn)行追蹤并計(jì)算節(jié)點(diǎn)的最小走時(shí)。以圖2模型為例，該模型共有上下兩個(gè)界面，將整個(gè)模型分為3個(gè)區(qū)域，假設(shè)震源和臺(tái)站均位于地表，現(xiàn)在需追蹤第一個(gè)反射界面的反射轉(zhuǎn)換波P1S1(數(shù)字1表示計(jì)算區(qū)域編號(hào)，上標(biāo)代表上行波，下標(biāo)代表下行波)。計(jì)算步驟為： ①選擇第一區(qū)為計(jì)算區(qū)域，調(diào)用P波速度模型從震源開(kāi)始計(jì)算直至一區(qū)所有網(wǎng)格全部計(jì)算完畢； ②保留一界面上的走時(shí)和射線路徑，同時(shí)選擇一區(qū)為計(jì)算區(qū)域，從界面一上搜索出走時(shí)最小點(diǎn)并將該點(diǎn)作為次級(jí)震源點(diǎn)，調(diào)用S波的速度模型，計(jì)算該區(qū)域的走時(shí)和射線路徑。按照上述原理可實(shí)現(xiàn)多震相地震射線的追蹤計(jì)算。

圖2 模型示意圖

3 反演算法

3.1 反演目標(biāo)函數(shù)

在速度與反射界面同時(shí)反演中，走時(shí)對(duì)界面深度的偏導(dǎo)數(shù)遠(yuǎn)小于走時(shí)對(duì)速度的偏導(dǎo)數(shù)，導(dǎo)致同時(shí)反演中界面基本不會(huì)更新。因此本文采用不同參數(shù)歸一化反演方法。多震相走時(shí)同時(shí)反演問(wèn)題可歸結(jié)為帶約束的阻尼最小二乘最優(yōu)化問(wèn)題，其目標(biāo)函數(shù)為

(4)

(μCm+ATCdA)Zmm=ATCdd

(5)

式中Cm、Cd為模型和數(shù)據(jù)空間的協(xié)方差矩陣。該式是對(duì)非線性問(wèn)題局部線性化的基本反演公式，其解具有局域解的特征。為了得到具有實(shí)際物理意義的解，可采用Zm、Cm及Cd的先驗(yàn)信息。因此，選擇不同組合的(Cm，Cd)可得到不同形式的反演解。若μ=0,Cd=I則式(5)變?yōu)榻?jīng)典的最小二乘問(wèn)題；若Cd=Cm=I，Zm=W，A=GW，這里W是表征射線寬度的帶寬矩陣，則式(5)變?yōu)镸eyerholtz等[26]提出的卷積壓制法。本文主要采用Tarantola等[27]提出的廣義帶約束的阻尼最小二乘反演解，即Cm、Cd分別為模型和數(shù)據(jù)空間的協(xié)方差矩陣的逆矩陣，其解為

(6)

該解的約束條件是

式(4)可用迭代的共軛梯度法進(jìn)行求解[28，29]，其關(guān)鍵是如何求解四面體單元中具有偏導(dǎo)性質(zhì)的Jacobi矩陣中的元素。

3.2 Jacobi矩陣的求解

根據(jù)射線理論可知，三維情況下沿某一射線路徑Rj的走時(shí)tj可用下式表示

(7)

式中：vc(x,y,z)為四面體單元內(nèi)速度函數(shù)； ds為射線的微分元； 模型參數(shù)化后，射線路徑Rj上的走時(shí)可寫成求和形式

(8)

式中：n為該射線所穿過(guò)的四面體單元總數(shù)；Rj,k表示穿過(guò)第k個(gè)四面體單元的射線長(zhǎng)度；vc,k(x,y,z)表示該四面體單元的速度分布。

速度與界面同時(shí)反演時(shí)的Jacobi矩陣包括兩部分：走時(shí)對(duì)速度變化的偏導(dǎo)數(shù)和走時(shí)對(duì)反射點(diǎn)深度變化的偏導(dǎo)數(shù)

(9)

式中：vk是第k個(gè)四面體單元的速度分布；tj是第j條射線的走時(shí)；zk是第k個(gè)反射點(diǎn)的深度值。

射線穿過(guò)某個(gè)四面體單元時(shí)，其走時(shí)對(duì)速度的一階偏導(dǎo)可用通過(guò)該四面體單元射線兩端點(diǎn)的平均值來(lái)代替[30]

(10)

其中

式中：ki和kj是穿過(guò)第k個(gè)四面體單元內(nèi)射線的兩個(gè)端點(diǎn)；vb(k)是第k個(gè)四面體單元內(nèi)四個(gè)主節(jié)點(diǎn)的速度值。

(11)

圖3 穿過(guò)某一個(gè)四面體網(wǎng)格的射線示意圖

式(10)中，界面p點(diǎn)處走時(shí)關(guān)于界面深度的導(dǎo)數(shù)為[20]

(12)

式中：γ1為入射向量與界面法向量的夾角；γ2為出射向量與界面法向量的夾角；v1,v2分別為入射向量側(cè)和出射向量側(cè)的速度。實(shí)際計(jì)算時(shí)，只需記錄入射點(diǎn)，反射點(diǎn)和出射點(diǎn)便可得到入射向量和出射向量，再求取在反射界面上反射點(diǎn)處的法向量即可得到走時(shí)對(duì)界面深度的偏導(dǎo)。

本文采用式(11)計(jì)算走時(shí)對(duì)速度的偏導(dǎo)，用式(12)計(jì)算走時(shí)對(duì)界面深度的偏導(dǎo)。在速度與界面的同時(shí)反演中，為了避免因反射界面更新引起某些區(qū)域的過(guò)度(或欠)更新，筆者對(duì)兩種不同類型的Jacobi矩陣元素進(jìn)行了歸一化處理。

4 數(shù)值模擬分析

4.1 計(jì)算效率分析

已知反演的耗時(shí)中正演的占比超過(guò)90%，所以同時(shí)反演的計(jì)算效率主要取決于選用的射線追蹤方法。為此對(duì)比了立方體模型參數(shù)化和四面體模型參數(shù)化下單炮正演的CPU計(jì)算時(shí)間，所選立方體模型的尺寸為80km×80km×48km。將主節(jié)點(diǎn)間距固定為8.0km，將次級(jí)節(jié)點(diǎn)間距從4.0km到0.25km逐步減半縮小，并記錄不同模型參數(shù)化下的計(jì)算時(shí)間(所用計(jì)算機(jī)：Dell，Intel(R) Core(TM)2 Duo CPU E7200，主頻2.52GHz)。

一般來(lái)說(shuō)，無(wú)論采用何種模型參數(shù)化方式，當(dāng)主節(jié)點(diǎn)的間距固定時(shí)，其正演的計(jì)算效率主要取決于次級(jí)節(jié)點(diǎn)的間距，但是隨著次級(jí)節(jié)點(diǎn)間距的減小，兩種模型參數(shù)化下的CPU時(shí)間差將大幅減小(表1，CPU時(shí)間差從次級(jí)節(jié)點(diǎn)間距為4.0km時(shí)的近200倍降到了次級(jí)節(jié)點(diǎn)間距為0.25km時(shí)的不到2倍)。從總節(jié)點(diǎn)數(shù)和CPU時(shí)間的關(guān)系圖(圖4)中也可得出類似結(jié)論。兩種模型參數(shù)化下的CPU時(shí)間都與次級(jí)節(jié)點(diǎn)間距成近似指數(shù)型的關(guān)系，但是立方體模型曲線的斜率更大。隨著總節(jié)點(diǎn)數(shù)目的增加，兩條曲線最后完全發(fā)散。

從圖4和表1可見(jiàn)，在需要高精度計(jì)算的情況(即更多的次級(jí)節(jié)點(diǎn))下，上述兩種模型參數(shù)化的計(jì)算效率處于同一水平。由于CPU時(shí)間中90%以上都被正演所消耗，因此這個(gè)結(jié)論也同樣適用于后面的同時(shí)反演部分。

圖4 CPU時(shí)間隨節(jié)點(diǎn)總數(shù)的變化關(guān)系曲線

次級(jí)節(jié)點(diǎn)間距km立方體模型節(jié)點(diǎn)總數(shù)/個(gè)四面體模型節(jié)點(diǎn)總數(shù)/個(gè)立方體模型計(jì)算時(shí)間/s四面體模型計(jì)算時(shí)間/s4.00 5135 23848 0.03 5.422.0025827864740.116.031.001156913503751.4812.030.50489339136527524.7079.970.2520123155471683390.33723.03

4.2 多震相走時(shí)聯(lián)合反演及同時(shí)反演

為了驗(yàn)證基于四面體網(wǎng)格單元下的反演算法的有效性，我們選擇了一個(gè)比較復(fù)雜的速度模型(圖2，模型尺度為80km×80km×48km)。地表起伏約為2km，共計(jì)3757個(gè)主節(jié)點(diǎn)。次級(jí)節(jié)點(diǎn)間距為1.0km，總共有613928個(gè)節(jié)點(diǎn)。共剖分出18432個(gè)四面體網(wǎng)格單元。在z=-30km附近有一個(gè)用正弦函數(shù)生成起伏界面，最大起伏為2km； 在z=-46km附近有一個(gè)上拱的弧形界面。30個(gè)震源(圖2中的灰色圓球)隨機(jī)分布在10～20km的深度范圍內(nèi)，289個(gè)檢波器(圖2中的倒三角形)均勻分布在地表(注：以下所有數(shù)值模擬實(shí)驗(yàn)中均采用此炮—檢排列)。圖5給出了單炮情況下的射線路徑分布圖。真實(shí)速度模型如圖6所示，可見(jiàn)水平面及垂直剖面上速度分布特征。下面討論初至波反演、初至波與反射波聯(lián)合反演、速度與反射界面同時(shí)反演。

4.2.1 多震相走時(shí)聯(lián)合反演成像

聯(lián)合反演中初始速度模型選用線性增加模型，反射界面為真實(shí)界面。反演結(jié)果均以與圖6相同的形式給出(圖7為水平切片，圖8為垂直切片)。其中圖7b和圖8b為僅用初至P波的反演結(jié)果，圖7c和圖8c為初至P波和P1P反射波聯(lián)合反演的結(jié)果，圖7d和圖8d為初至P波、P1P反射波和P2P反射波聯(lián)合反演的結(jié)果，圖9 (水平切片)和圖10(垂直切片)分別為反演結(jié)果與真實(shí)模型的相對(duì)誤差。從圖7～圖10可以看出，僅用初至P波進(jìn)行成像在淺層可大體恢復(fù)速度場(chǎng)的起伏形狀，但其成像深度卻明顯不及聯(lián)合反演。在加入P1P反射波之后，由于增加了射線交錯(cuò)概率，同時(shí)又提高了射線密度，從水平切片結(jié)果(圖7c)來(lái)看，聯(lián)合反演已基本能準(zhǔn)確地反演速度場(chǎng)的異常圖像和異常幅度；從垂直切片的結(jié)果(圖8c)來(lái)看，其成像深度明顯增加；從圖9c來(lái)看，模型邊緣的速度場(chǎng)起伏的恢復(fù)效果也比初至P波的(圖9b)有明顯提升。加入P2P反射波之后，從兩個(gè)切片的結(jié)果(圖7d和圖8d)來(lái)看，速度場(chǎng)的恢復(fù)與真實(shí)模型已基本相差無(wú)幾，且從圖9d和圖10d可見(jiàn)，三種波聯(lián)合反演的精度最高。此數(shù)值模擬實(shí)驗(yàn)表明，基于四面體網(wǎng)格的多震相走時(shí)聯(lián)合反演可更為準(zhǔn)確地反演出異常體形狀和異常幅度，即效果明顯優(yōu)于單一震相的反演。

圖5 射線路徑示意圖

圖6 三維速度模型

4.2.2 多震相走時(shí)速度和界面同時(shí)反演成像

在速度與界面的同時(shí)反演中采用的初始速度模型為向下線性增加模型，實(shí)際中初始反射界面可通過(guò)一些先驗(yàn)信息獲得。這里初始反射界面我們選為水平界面(上界面為z=-30km，下界面為z=-47km)。對(duì)兩種不同的Jacobi偏導(dǎo)元素進(jìn)行了最大值歸一化處理。圖11a和圖12a為初至P波和P1P反射波聯(lián)合同時(shí)反演的速度結(jié)果，圖11b和圖12b為初至P波、P1P和P2P聯(lián)合同時(shí)反演的速度結(jié)果，圖11c、圖11d、圖12c、圖12d為對(duì)應(yīng)的相對(duì)誤差。圖13b為初至波和P1P反射波聯(lián)合同時(shí)反演的上界面結(jié)果，圖13c、圖13e分別為初至波、P1P和P2P聯(lián)合同時(shí)反演的上、下界面的反演結(jié)果。

從速度反演結(jié)果的水平切片來(lái)看，同時(shí)反演已大體恢復(fù)速度場(chǎng)的起伏形態(tài)，由于模型邊界附近相比模型中央?yún)^(qū)域數(shù)據(jù)覆蓋密度要小，而速度與界面的耦合關(guān)系又決定了界面的恢復(fù)情況必將影響速度的恢復(fù)，所以在水平切片(圖11)中可看出模型中央速度恢復(fù)較好，但邊緣區(qū)域恢復(fù)稍差，有些許失真，這種失真現(xiàn)象在雙界面的同時(shí)反演(圖11b、圖11d)中由于較多的數(shù)據(jù)覆蓋而使之得到了一些改善。

圖7 不同震相走時(shí)聯(lián)合反演成像的水平切片(z=-16km)

圖8 不同震相走時(shí)聯(lián)合反演成像的垂直剖面(y=40km)

從垂直剖面(圖12)中可看出，單界面同時(shí)反演時(shí)-30km以上的速度場(chǎng)的起伏已經(jīng)基本恢復(fù)，但在模型的邊緣有些許失真(圖12a、圖12c)。雙界面的同時(shí)反演(圖12b、圖12d)已經(jīng)基本恢復(fù)速度場(chǎng)的起伏形態(tài)，相比于單界面的情況，模型邊緣及下方速度場(chǎng)的失真情況也有了明顯改善。

從圖13中可見(jiàn)，無(wú)論是初至P波與P1P波同時(shí)反演，還是初至P波與P1P、P2P波同時(shí)反演，上界面的形態(tài)都已基本恢復(fù)。由于在雙界面的同時(shí)反演中模型邊界附近的數(shù)據(jù)密度大于單一界面同時(shí)反演情況下的數(shù)據(jù)密度，而大數(shù)據(jù)密度意味著對(duì)速度和界面更好的約束，也意味著更好的反演結(jié)果，這一點(diǎn)可從圖13b與圖13c的對(duì)比中看出，無(wú)論是幅值的恢復(fù)情況，還是界面形態(tài)的相似度，圖13c的結(jié)果都要比圖13b好。從雙界面的同時(shí)反演結(jié)果中已基本可見(jiàn)下界面近似弧形的形態(tài)特征，但其幅值的恢復(fù)不如上界面，這是因?yàn)樵谏辖缑嬷嫌兄边_(dá)波信息，可幫助更加準(zhǔn)確地反演出速度分布，進(jìn)而更好地反演上界面的形態(tài)，而處在上界面和下界面之間的區(qū)域則沒(méi)有這個(gè)優(yōu)勢(shì)，所以在雙界面的同時(shí)反演中，上界面的反演結(jié)果優(yōu)于下界面。從這個(gè)數(shù)值模擬實(shí)驗(yàn)中可得出以下結(jié)論：基于四面體網(wǎng)格的同時(shí)反演得到的速度和界面結(jié)果較為可信。

圖9 反演水平切片與真實(shí)模型的相對(duì)誤差(z=-16km)

圖10 反演垂直剖面和真實(shí)模型的相對(duì)誤差(y=40km)

圖11 不同震相走時(shí)聯(lián)合同時(shí)反演成像水平切片及相對(duì)誤差(z=-16km)

5 結(jié)果分析與討論

本文基于四面體模型參數(shù)化下的多震相射線追蹤技術(shù)，結(jié)合共軛梯度法求解帶約束的阻尼最小二乘問(wèn)題反演算法，討論并實(shí)現(xiàn)了四面體模型參數(shù)化下多震相走時(shí)聯(lián)合及同時(shí)反演技術(shù)。此外還對(duì)四面體模型參數(shù)化下的計(jì)算精度和CPU時(shí)間做了詳細(xì)討論，可知在總節(jié)點(diǎn)數(shù)目大致相同的情況下，其CPU時(shí)間遠(yuǎn)小于立方體模型下的計(jì)算時(shí)間。在同時(shí)反演中為了平衡速度和界面對(duì)走時(shí)的影響，本文對(duì)兩種不同類型的Jacobi元素進(jìn)行了最大值歸一化處理。數(shù)值模擬實(shí)驗(yàn)表明，四面體模型參數(shù)化下的多震相走時(shí)聯(lián)合成像技術(shù)具有很好的反演能力，能較準(zhǔn)確地反演速度場(chǎng)的起伏特征和分布情況；對(duì)于速度與界面的同時(shí)反演，可較好地反演速度場(chǎng)的起伏特征，并可基本恢復(fù)反射界面的幾何形態(tài)。在現(xiàn)有的地震臺(tái)站分布極不均勻的情況形下，四面體網(wǎng)格單元可對(duì)數(shù)據(jù)密度不同的區(qū)域進(jìn)行不同尺度的剖分，進(jìn)而在走時(shí)反演中更顯著地提高數(shù)據(jù)的利用效率，且四面體網(wǎng)格單元對(duì)起伏地表和異常體的擬合精度明顯優(yōu)于立方體單元網(wǎng)格，故該算法具有廣泛的應(yīng)用價(jià)值。