黃國嬌， 白超英， 錢衛(wèi)

1 河海大學(xué)地球科學(xué)與工程學(xué)院地質(zhì)科學(xué)與工程系， 南京 211100 2 長安大學(xué)地質(zhì)工程與測繪學(xué)院地球物理系， 西安 710054

?

球坐標(biāo)系下多震相走時三參數(shù)同時反演成像

黃國嬌1， 白超英2， 錢衛(wèi)1

1 河海大學(xué)地球科學(xué)與工程學(xué)院地質(zhì)科學(xué)與工程系， 南京 211100 2 長安大學(xué)地質(zhì)工程與測繪學(xué)院地球物理系， 西安 710054

球坐標(biāo)系下多震相走時三參數(shù)(速度、震源位置和反射界面)同時反演需要解決兩個關(guān)鍵問題：(1)球坐標(biāo)系下3D速度模型中多次透射、反射(折射)及轉(zhuǎn)換波精確、快速的射線追蹤；(2)同時反演時三種不同參數(shù)間的強(qiáng)耦合問題.為此，我們將直角坐標(biāo)系下分區(qū)多步不規(guī)則最短路徑算法推廣至球坐標(biāo)系中，進(jìn)行區(qū)域或者全球尺度的多震相射線追蹤.然后將其與適合多參數(shù)同時反演的子空間算法相結(jié)合，形成一種球坐標(biāo)系下聯(lián)合多震相走時三參數(shù)同時反演的方法技術(shù).與雙參數(shù)(速度和反射界面或速度和震源位置)同時反演的數(shù)值模擬對比分析顯示：三參數(shù)與雙參數(shù)的同時反演結(jié)果大體接近，并且它們對到時數(shù)據(jù)中可容許的隨機(jī)噪聲不太敏感.結(jié)果說明本文中的同時反演成像為一種提高成像分辨率，同時反演速度、震源位置和反射界面的有效方法.

球坐標(biāo)系; 分區(qū)多步不規(guī)則最短路徑算法; 多震相走時; 三參數(shù)同時反演; 子空間法

1 引言

自1976年Aki和Lee將醫(yī)學(xué)層析成像技術(shù)應(yīng)用于地震層析成像研究以來，地震層析成像技術(shù)經(jīng)過近40年的不斷發(fā)展和完善，已成為研究地球內(nèi)部結(jié)構(gòu)、以及介質(zhì)各向異性行之有效的方法技術(shù)之一.其成像尺度涵蓋范圍很廣，小到巖石組分，大至全球構(gòu)造.對于區(qū)域或全球?qū)游龀上穸裕厍蚯什豢珊雎?，地球表面不能簡單視作平面來處?Snoke and Lahr, 2001), 而必須采用球坐標(biāo)系才能保證所需的計算精度.同時，地球內(nèi)部的結(jié)構(gòu)是復(fù)雜變化的，需用三維球坐標(biāo)系下的速度模型來描述.在速度異常(-2%)的南太平洋地區(qū)，Zhao和 Lei (2004)分別采用3D速度模型和1D速度模型計算走時和射線路徑，結(jié)果顯示PcP震相的射線路徑改變了77 km，其他震相的射線路徑偏差也達(dá)到了幾十千米，走時也存在數(shù)秒的差異.因此，如果所研究區(qū)域達(dá)到一定的尺度時(例如區(qū)域、特別是全球走時成像)，不考慮地球曲率的影響和介質(zhì)的橫向非均勻性，計算的射線路徑和走時就會包含較大的誤差.

目前應(yīng)用較廣的走時層析成像，高頻假設(shè)等因素的制約導(dǎo)致其成像的空間分辨率較低.因此，在現(xiàn)有地震和地震臺網(wǎng)分布情況下，采用多震相走時成像是提高分辨率非常有效的途徑.多震相走時成像較為典型的一個例子是Zhao等(2005)利用美國Landers地震余震進(jìn)行多震相(S，SmS，sSmS)走時同時反演S波速度和Moho面，其中僅采用相距200 km的兩個觀測臺站和其間發(fā)生的180個余震，其結(jié)果表明使用多震相走時成像其空間分辨率可提高至臺站間距的1/8～1/6(即25～35 km).

在多震相走時成像中，由于反射或反射轉(zhuǎn)換波震相的走時由震源參數(shù)、3D速度分布和反射界面共同決定，且這三類參數(shù)反演前往往是未知的，并存在強(qiáng)耦合關(guān)系(如：震源位置較小的變化將引起射線路徑、以及射線入射到界面上反射點位置的較大變化，反過來說，反射點位置較小的變化也將同樣引起射線路徑的變化等)，因此要求反演算法應(yīng)具有同時反演三參數(shù)(速度、震源參數(shù)、反射界面)的能力.

對于區(qū)域或者全球走時層析成像，為保證計算精度和提高成像分辨率，需要進(jìn)行三維球坐標(biāo)系下多震相走時三參數(shù)同時反演成像.由于全球走時成像中需要成千上萬的走時數(shù)據(jù)(大多采用不僅僅是單一震相的走時)，且缺少高效和高精度的三維正演射線追蹤計算方法及相應(yīng)軟件，故目前采用較多的還是基于1D球坐標(biāo)系下的兩點射線追蹤算法，此算法反演中假設(shè)射線路徑不變，然后采用擾動的方法更新速度模型(Van der Hilst et al., 1997; Su et al., 1994; Grand et al., 1997).近些年來人們已意識到采用三維球坐標(biāo)系的重要性，相繼出現(xiàn)了幾種不同的射線追蹤方法，例如：Bijwaard和Spakman(1999)研發(fā)的快速運動學(xué)射線追蹤算法，在網(wǎng)格算法中引入了彎曲法，能夠追蹤3D球坐標(biāo)系下的(臨界)折射、多次反射、以及一些繞射(Pdiff)波，且計算的各主要震相絕對走時誤差約為0.1 s，但實質(zhì)上該算法依然是兩點間的射線追蹤.Zhao和Lei (2004)將球坐標(biāo)系下的偽彎曲射線追蹤算法(Zhao et al., 1992)進(jìn)行了推廣，可進(jìn)行區(qū)域或全球多震相走時反演成像.Tian等(2007)開發(fā)了一套可用于全球有限頻成像的動力學(xué)射線追蹤程序，所計算的主要震相走時絕對誤差約為0.1 s，但其算法依然是基于兩點間的射線追蹤問題(局域解問題)，計算速度仍然存在問題(每個震源和臺站都需要重新進(jìn)行射線追蹤).

目前為止，走時三參數(shù)(速度、反射界面和震源參數(shù))同時反演研究并不多見，Sambridge (1990)雖然是首位討論走時三參數(shù)同時反演的研究學(xué)者，反演采用子空間反演方法(Kennett et al., 1988)，但其正演采用了傳統(tǒng)的邊界值兩點射線追蹤算法(Sambridge and Kennett, 1990)；Preston (2003)在他的博士論文中也進(jìn)行了走時三參數(shù)同時反演的研究工作，其研究成果已在Science上發(fā)表(Preston et al.，2003)，其中利用正則化加權(quán)處理方式對三種參數(shù)進(jìn)行解耦，容易引入人為因素.在直角坐標(biāo)系下，黃國嬌和白超英(2013)采用分區(qū)多步最短路徑算法(Multiple ISPM)進(jìn)行多震相的射線追蹤，反演采用子空間算法，實現(xiàn)了3D多震相走時三參數(shù)同時反演成像.

鑒于此，針對區(qū)域或全球?qū)游龀上瘢疚膶⑸鲜龇謪^(qū)多步最短路徑算法推廣至球坐標(biāo)系中，進(jìn)行多震相走時和射線路徑的追蹤計算.然后將其與地震多參數(shù)反演算法(子空間反演算法)相結(jié)合，提出了一種3D球坐標(biāo)系下可利用地震多震相走時資料進(jìn)行三參數(shù)同時反演的方法技術(shù).為驗證該三參數(shù)同時反演算法的有效性和正確性，我們選用一個區(qū)域模型進(jìn)行數(shù)值模擬實驗，在合成理論走時數(shù)據(jù)中加入隨機(jī)噪聲并對其進(jìn)行三參數(shù)同時反演，將反演成像結(jié)果與雙參數(shù)同時反演進(jìn)行對比分析.結(jié)果表明，三參數(shù)和雙參數(shù)同時反演方法技術(shù)具有相近的反演能力，都能較好地反演出震源位置、速度分布、界面形態(tài)和位置.

2 多震相射線追蹤算法

2.1 參數(shù)化

利用分區(qū)多步不規(guī)則最短路徑算法進(jìn)行射線追蹤之前，需要將模型參數(shù)化.三維球坐標(biāo)系下模型參數(shù)化如圖1所示，首先根據(jù)速度間斷面的起伏形態(tài)將模型分為四個區(qū)(圖1a)，然后對每個區(qū)域用規(guī)則單元或者不規(guī)則單元將其參數(shù)化.地心點所在單元形成如圖1c所示的五面體，圖1d為兩極極點所在單元，而不包含地心、極點且遠(yuǎn)離起伏地表和界面的區(qū)域則可采用如圖1b所示的規(guī)則六面體單元參數(shù)化.單元的角點為主節(jié)點，主節(jié)點間等距離的插入次級節(jié)點(這點很重要，因為越往地心，相同扇形的單元弧長越短)，次級節(jié)點只分布在網(wǎng)格單元的邊上，單元內(nèi)無節(jié)點，但震源和地震臺站可以位于其中.

2.2 速度插值

速度采樣在主節(jié)點上進(jìn)行，次級節(jié)點的速度可由主節(jié)點上的速度值通過三次線性插值函數(shù)(3D)得到.當(dāng)節(jié)點(震源或者地震臺站)位于單元內(nèi)部時，需分規(guī)則單元和非規(guī)則單元兩種情況對待：

(1) 當(dāng)節(jié)點位于規(guī)則單元中，如圖1b所示，節(jié)點速度則由所在單元八個主節(jié)點上的速度值通過三線性拉格朗日插值函數(shù)給出，插值函數(shù)為：

(1)

(2) 若節(jié)點位于不規(guī)則單元內(nèi)，如圖1c—1d所示，則需先將單元剖分為幾個四面體(四面體的四個角點必須是主節(jié)點，編號1—4)，然后再由節(jié)點所在四面體的四個主節(jié)點上定義的四面體體積插值函數(shù)給出(Huang et al., 2013)：

(2)

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

圖1 球坐標(biāo)系下模型參數(shù)化示意圖(圖a, 為清晰起見, 次級節(jié)點未標(biāo)出)和所用參數(shù)化單元(圖b—d), 圖b—d中大球表示主節(jié)點, 小球表示所加入的次級節(jié)點Fig.1 Model parametrization in 3D spherical coordinates (diagram a, for clear representation, the secondary nodes are not depicted) and three kinds of cells used to divide the spherical model (diagram b, c and d, in the figure the larger circles are primary nodes and the smaller circles are the secondary nodes, respectively)

分區(qū)多步計算首先是模型分區(qū)，然后按照給定的震相，自震源所在區(qū)域開始，沿著目標(biāo)射線通過的區(qū)域分段分區(qū)進(jìn)行射線追蹤計算.舉例來說，如圖1a中所示，根據(jù)反射界面將模型分為了四個區(qū)，速度不連續(xù)(間斷)面連接相鄰區(qū)域.若計算圖1a中界面一(深度660 km)的透射或反射震相時，首先自震源所在的單元開始波前擴(kuò)展計算，然后按照一定的步長進(jìn)行第一區(qū)的波前擴(kuò)展計算.等第一區(qū)所有節(jié)點上的走時計算完成之后，波前停留在第一區(qū)的速度離散界面上，并保存了震源到界面離散點上的最小走時及相應(yīng)的射線路徑.此時如果計算透射波P1P2(上、下標(biāo)分別代表上行波和下行波，1和2表示區(qū)號)，則調(diào)用第二區(qū)的P波速度(如果追蹤P1S2震相則調(diào)用第二區(qū)的S波速度)，自界面一上走時最小的點開始進(jìn)行第二區(qū)向下的波前擴(kuò)展掃描計算；若是計算反射波P1P1(或反射轉(zhuǎn)換波P1S1)，同樣自界面一上走時最小的點開始，并調(diào)用第一區(qū)相應(yīng)的速度模型向上進(jìn)行波前擴(kuò)展計算.同理，依據(jù)多次波是由簡單的入射、透射、反射及轉(zhuǎn)換波按照一定的規(guī)律組合而成的原理，即可實現(xiàn)任意多震相地震射線的追蹤計算.圖2是3D球坐標(biāo)系下采用上述分區(qū)多步不規(guī)則最短路徑算法追蹤計算的多震相射線路徑圖(所用速度模型和反射界面位置詳見第4節(jié)反演實例)，更復(fù)雜的多次波均可采用該算法進(jìn)行追蹤計算，但是需要在輸入文件中給出多次波的傳播途徑，即射線穿過哪些區(qū)域，與哪些界面相交，遇到界面是反射、折射還是轉(zhuǎn)換.

2.4 計算精度和效率分析

Huang等(2013)將該算法計算的49種不同的全球主要震相走時與AK135全球走時表(Kennet et al., 1995)進(jìn)行對比分析，討論該算法的計算精度問題.結(jié)果表明絕大多數(shù)震相的最大絕對走時誤差小于0.1 s.在全球或者區(qū)域走時層析成像中，還要求正演算法具有較高的計算效率.因此，我們分析了圖2所示模型中三種不同震相(P，P410P和P660P)在不同尺度單元參數(shù)化和次級節(jié)點情況下的計算效率，結(jié)果如圖3所示.其中震源和臺站的個數(shù)及位置分布均與第4節(jié)的反演實例中相同，單元大小分別為4°×4°×200 km(圖3a)、2°×2°×100 km(圖3b)和1°×1°×50 km(圖3c)，計算機(jī)配置為Intel i5-2300CPU處理器、2.8 GHz的主頻和4 G內(nèi)存.

圖2 球坐標(biāo)系下多震相射線路徑圖淺灰色線：初至P波；灰線：深度410 km處界面的反射波P410P；黑線：深度660 km處界面的反射波P660P.Fig.2 Multiple ray paths in a spherical coordinateThe directed arrivals (light grey lines), the reflection from the 410 km subsurface (gray lines) and the reflections from 660 km subsurface (black lines) are shown.

圖3 計算多震相(P，P410P和P660P)的CPU時間圖(a) 4°×4°×200 km；(b) 2°×2°×100 km；(c) 1°×1°×50 km.圖中括號內(nèi)數(shù)字為當(dāng)前模型參數(shù)化下節(jié)點總數(shù).Fig.3 Results of computational efficiency tests on tracing the multiple phases ( P, P410P and P660P)(a) Cell scale of 4°×4°×200 km; (b) Cell scale of 2°×2°×100 km; (c) Cell scale of 1°×1°×50 km. In the figures, the numbers in the parentheses indicate the total number of nodes for the corresponding parameterized model.

從圖3c中可以看出，當(dāng)單元大小為1°×1°×50 km，加入11個次級節(jié)點(地表次級節(jié)點間隔約為8 km)時，計算這三種震相幾乎要花費近8個小時.但是，通常加入5～7個次級節(jié)點(次級節(jié)點間距為17～22.5 km)就能滿足精度需求，此時CPU時間就減少至1.0～1.5個小時.值得注意的是，反演成像中每次迭代基本上有95%的時間花費在計算這三種震相的走時和射線路徑上.因此，我們可大致估計采用上述三種震相同時反演三參數(shù)的時間.

3 反演算法

三參數(shù)的同時反演是一項很困難的工作，因為三者具有不同的物理量綱，并且具有強(qiáng)耦合關(guān)系.若在同一個方程中求解這三種不同的模型參數(shù)，解會偏向系數(shù)大的一方，很可能導(dǎo)致反演的失敗.因此，如何解耦成為三參數(shù)同時反演中急需解決的問題.直角坐標(biāo)系中，我們將共軛梯度法求解帶約束的阻尼最小二乘算法(DMNCL-CG)和子空間法(Subspace)進(jìn)行了對比，結(jié)果表明Subspace更適合于多參數(shù)的同時反演.因此，本文采用Subspace算法對球坐標(biāo)系中三參數(shù)進(jìn)行同時反演.

3.1 子空間算法(Subspace)

多震相走時同時反演可歸結(jié)為求下列目標(biāo)函數(shù)的極小值問題(Rawlinson and Sambridge, 2003)：

(3)

子空間法是在模型的若干個子空間同時沿著多個方向去尋找目標(biāo)函數(shù)S(m)的最小值.其算法的核心是將多種參數(shù)的模型擾動量用p維子空間的一組基向量{aj}，j=1,2,…,p來表示，即

(4)

(5a)

(5b)

其中Frechet偏導(dǎo)數(shù)矩陣G=?g/?m.將表達(dá)式(5a)(5b)代入公式(4)中就可得到反演每次迭代時模型的更新量：

(6)

第二步：若p≥1，則采用下列公式構(gòu)建其余的基向量：

?

第三步：采用奇異值分解法(SVD)將第一步和第二步中構(gòu)建的基向量{aj}標(biāo)準(zhǔn)正交化.

3.2 Frechet偏導(dǎo)數(shù)計算

三維球坐標(biāo)系中，地球內(nèi)任一點P可表示為(r,θ,φ)，實際應(yīng)用中通常用徑向深度、緯度和經(jīng)度來表示，記為(z,Lat,Long)，并且有r=R-z(R=6371.5 km為地球半徑)；θ=90°-Lat和φ=Long.根據(jù)射線理論，由第i個震源經(jīng)反射點到達(dá)第j個檢波器的走時可表述為下列積分：

(7)

其中ds為沿路徑Rij上的線積分元，而Vc(r,θ,φ)則為沿射線路徑上的速度分布函數(shù).當(dāng)模型網(wǎng)格單元化后，沿該射線路徑Rij上的走時可寫成求和形式：

(8)

其中：n為射線穿過的單元總數(shù)，Rij,k、Vc,k(r,θ,φ)分別為射線穿過第k個單元的長度和速度分布.

在速度、界面和震源位置同時反演時，F(xiàn)rechet偏導(dǎo)矩陣包含了三部分：一是走時關(guān)于速度變化的導(dǎo)數(shù)，二是走時關(guān)于反射點深度變化的導(dǎo)數(shù)，三是走時關(guān)于震源參數(shù)的導(dǎo)數(shù).

(9)

其中vk和rk分別是第k個單元內(nèi)的速度分布和第k個反射點的深度，Δvk和Δrk則分別是相應(yīng)速度和反射點深度的擾動量，Δxk(k=1,2,3)是第k個震源分別在r，θ和φ方向的擾動量.

(9)式中第一項，即有關(guān)走時關(guān)于速度的一階偏導(dǎo)數(shù)計算如下(為簡便起見，這里略去震源編號i，則tj是第j條射線Rj的走時)：

(10)

(11)

其中rn是離散界面節(jié)點的深度坐標(biāo)，rint是反射點的深度坐標(biāo)(有關(guān)wj,wj+1，wn和wr，見圖4所示)，球坐標(biāo)系中wr=(1,0,0)，?rint/?rn根據(jù)界面節(jié)點上的深度插值函數(shù)計算.vj和vj+1分別是界面兩側(cè)的速度值.當(dāng)僅考慮反射時，wj+1·wn=-wj·wn和vj=vj+1，則(16)式簡化為：

(12)

(9)式的第三項，即走時關(guān)于震源參數(shù)(位置和發(fā)震時刻)變化的偏導(dǎo)數(shù)可由以下步驟計算，首先走時(t=ta-to)關(guān)于發(fā)震時刻的導(dǎo)數(shù)為：

(13)

而走時關(guān)于震源深度變化的偏導(dǎo)數(shù)的計算可參見圖5所示.假設(shè)震源的深度擾動值為Δrh，則走時關(guān)于震源深度變化的偏導(dǎo)數(shù)可表示為：

(14)

圖4 走時關(guān)于反射點深度偏導(dǎo)數(shù)計算示意圖Fig.4 Diagrammatic showing the calculation of traveltime derivative with respect to reflector depth change

圖5 走時關(guān)于震源位置變化偏導(dǎo)數(shù)計算示意圖Fig.5 Diagrammatic showing the calculation of traveltime derivative with respect to source location change

同理可得走時關(guān)于震源位置其余兩個方向(Lat,Long)變化的偏導(dǎo)數(shù)表達(dá)式：

(15)

其中αh如圖5所示，代表射線與r軸的夾角；βh和γh則是射線在θφ平面投影后分別與θ和φ軸的夾角.

4 數(shù)值模擬分析

為檢驗該三參數(shù)同時反演算法的有效性及正確性，我們選擇一個尺度為20°×20°×700km的區(qū)域模型，如圖6所示.圖6a是經(jīng)度-深度剖面的起伏速度圖，速度沿著緯度方向不變化.模型中引入兩個起伏反射界面模擬地幔中兩個速度不連續(xù)面(如圖6b所示)，深度分別位于410km和660km處，起伏幅度均為±20km.模型參數(shù)化選用1°×1°×50km的單元，次級節(jié)點間距為10km(相當(dāng)于在單元面的不同方向加了9個次級節(jié)點).30個震源隨機(jī)分布在深度20km至100km之間(圖7a是震源在經(jīng)緯度面上的投影；圖7b則是震源在經(jīng)度-深度剖面上的投影)，441個檢波器均勻(間距為1°×1°)分布于地表(模型頂面).

本文的同時反演中采用了3種不同的震相：直達(dá)P波、反射界面410km和660km的一次反射波(P410P、P660P).為檢驗反演算法抗噪聲能力，三種震相走時數(shù)據(jù)中加入了隨機(jī)噪聲(直達(dá)P波、P410P和P660P分別加入±0.2s、±0.3s和±0.3s的隨機(jī)誤差).選用背景層狀速度作為初始速度模型，兩個初始界面均為水平(深度分別為410km和660km).為了模擬實際地震定位中存在誤差這一現(xiàn)實，初始震源位置的選擇則是分別對實際震源位置在三個方向(經(jīng)度、緯度、深度)進(jìn)行擾動后得到，其相應(yīng)擾動最大幅度分別為±2°、±2°和±10km.

圖6 (a)某一緯度時, 經(jīng)度-深度剖面速度圖; (b)兩個起伏反射界面圖Fig.6 (a) Velocity model showing one cross section and (b) two undulating subsurface interfaces

圖7 震源在經(jīng)緯度剖面(a)和經(jīng)度-深度剖面(b)投影Fig.7 The sources projected on the horizontal plane (a) and the longitude-depth cross section (b)

為了進(jìn)行系統(tǒng)的對比研究，模擬實驗中進(jìn)行了三種不同的反演，其一是速度和界面的同時反演(簡記為：V+R反演)，其二是速度和震源位置的同時反演(簡記為：V+S反演)，其三是速度、界面和震源位置的同時反演(簡記為：V+R+S反演).不失一般性，三種不同的反演選用相同的反演參數(shù)，即阻尼因子和子空間維數(shù)分別為0.001和6(對2～8之間逐一嘗試，選擇較合適的子空間維數(shù)為6).迭代均為15次，相應(yīng)的走時殘差分別由初始時的1.834s，9.15s，9.269s下降至0.087s(V+R反演), 0.25s(V+S反演), 0.289s(V+R+S反演).

圖8給出了三種不同反演算法在Lat=10°垂直剖面上速度場的反演結(jié)果.從中可以看出，三種不同的反演算法重建的速度場分布與真實速度(圖8a)具有相似的異常形態(tài)，但是異常起伏的幅度似乎未完全恢復(fù)，且模型下部的更新相對較差，其主要原因是模型下部只有P660P震相的射線穿過，射線覆蓋率和交叉概率有限所致.對比圖8b—8d，會發(fā)現(xiàn)相對于雙參數(shù)的反演結(jié)果(圖8b和8c)，三參數(shù)同時反演的結(jié)果(圖8d)要略差些，但是差異不明顯，仍然能夠反演出異常的形態(tài)和結(jié)構(gòu).

為更加直觀顯示反演更新的速度場與真實速度場之間的差異，圖9給出了與圖8對應(yīng)的垂直剖面上反演的速度場與真實速度場之間的百分比差異.從圖9中可以看出，雖然三種反演方法都能重建出速度異常的形態(tài)和結(jié)構(gòu)，但是速度異常幅度并未完全恢復(fù).并且，可以明顯地看出雙參數(shù)同時反演結(jié)果要略好于三參數(shù)同時反演結(jié)果，特別是模型下部速度異常的重建結(jié)果.

圖10為兩種反演算法對反射界面的更新結(jié)果.可以看出，兩種反演算法都能較好地更新反射界面的形態(tài)和位置，與真實反射界面(圖6b)很相似，但僅從圖10中無法判斷哪種算法的更新結(jié)果更好.因此，為進(jìn)一步評價反射界面的更新結(jié)果，我們給出了兩種反演算法更新的反射界面與真實反射界面間的平均差異(表1).由于三參數(shù)(V+R+S)同時反演的難度遠(yuǎn)大于雙參數(shù)(V+R)同時反演的難度，故從表1中可知，雙參數(shù)(V+R)同時更新反射界面的結(jié)果要略好于三參數(shù)(V+R+S)同時更新的結(jié)果，但通過圖10的對比可得出三參數(shù)也能較好地更新反射界面幾何形態(tài)和位置.

圖8 三種反演算法在Lat=10°垂直剖面的成像結(jié)果(a) 真實模型; (b) V+R; (c) V+S; (d) V+R+S.Fig.8 The reconstructed velocity fields in Lat=10 ° cross section with three kinds of inversion methods(a) Real velocity model; (b) V+R; (c) V+S; (d) V+R+S.

圖9 三種不同反演算法的反演結(jié)果在Lat=10°垂直剖面上與真實速度的百分比誤差(a)初始速度模型; (b) V+R反演結(jié)果; (c) V+S反演結(jié)果; (d) V+R+S反演結(jié)果.Fig.9 The difference (in percentage) between the initial and true velocity model in Lat=10° cross section, between the true velocity and the inverted velocity with three different inversion methods in the same cross section(a) Initial velocity model; (b) Results for the V+R inversion; (c) Results for the V+S inversion; (d) Results for the V+R+S inversion.

圖10 兩種不同反演算法更新的反射界面(a) V+R反演結(jié)果; (b) V+R+S反演結(jié)果.Fig.10 The updated reflected interfaces by the two different inversion algorithms(a) Results for the V+R inversion; (b) Results for the V+R+S inversion.

圖11 兩種不同反演算法得到震源的位置(a, d, g)、水平位置(b, e, h)和深度(c, f, i)分別相對于真實震源位置的收斂情況(a, b, c)初始值; (d, e, f) V+S反演結(jié)果; (g, h, i) V+R+S反演結(jié)果.DZ表示在深度方向上距離真實震源的距離.Fig.11 The convergence to the real source locations (a, d, g), the real horizontal position (b, e, h) and the real focus depths (c, f, i) by two different algorithms(a, b, c) Initial; (d, e, f) Results for the V+S inversion; (g, h, i) Results for the V+R+S inversion.

表1 界面深度平均絕對誤差統(tǒng)計表Table 1 Mean absolute error statistics of the inverted depth for the two interfaces

圖11是兩種不同反演算法對震源位置的更新結(jié)果.這里我們用三種收斂情況來評價震源位置的更新結(jié)果，更新的震源位置相對于真實震源位置距離的收斂情況(圖11左)、更新震源的水平位置相對于真實震源水平位置距離的收斂情況(圖11中)和更新的震源深度相對于真實震源深度的收斂情況(圖11右).從圖11中可以看出，雙參數(shù)的反演結(jié)果(圖11d, 11e, 11f)要略好于三參數(shù)的反演結(jié)果(圖11g, 11h, 11i).無論是雙參數(shù)(V+S)的同時反演結(jié)果(圖11d—11f)，還是三參數(shù)(V+R+S)的同時反演結(jié)果(圖11g—11i)，其中對震源水平位置(圖11e和11h)的更新結(jié)果均明顯好于震源深度(圖11f和11i)的結(jié)果.可知在該模型下存在震源深度的定位誤差大于水平面內(nèi)的定位誤差的問題，其主要原因是對于區(qū)域或者全球成像而言，雖然利用了反射波走時資料，但模型尺度相對較大，走時對震源深度較小的變化不敏感，從而導(dǎo)致震源深度的更新相對較差.然而，在初始震源偏離真實震源較遠(yuǎn)(此例中最大距離約為320 km)的情況下，三參數(shù)的同時反演中，對震源的更新結(jié)果還是較為理想的.

5 結(jié)論

對于區(qū)域或全球?qū)游龀上?，為保證計算精度、提高成像分辨率，本文將直角坐標(biāo)系下的多震相射線追蹤正演(Multistage ISPM)算法推廣至球坐標(biāo)系中，并與比較流行的多參數(shù)反演算法(Subspace)相結(jié)合，提出一種3D球坐標(biāo)系中多震相走時三參數(shù)同時反演的方法技術(shù).數(shù)值模擬實例表明，三參數(shù)同時反演方法技術(shù)和雙參數(shù)反演方法具有相近的反演能力，都能較好地反演出震源位置、速度分布、界面形態(tài)和位置.另外，無論是雙參數(shù)或是三參數(shù)同時反演算法均對走時數(shù)據(jù)中的隨機(jī)誤差不太敏感，這點在實際應(yīng)用中也是至關(guān)重要的.因此，本文提出的球坐標(biāo)系下的三參數(shù)同時反演算法是一種利用多震相走時資料同時反演速度、震源位置以及反射界面行之有效且實用性較強(qiáng)的方法技術(shù).

Aki K, Lee W H K. 1976. Determination of three-dimensional velocity anomalies under a seismic array using first P arrival times from local earthquakes: 1. A homogeneous initial model.J.Geophys.Res., 81(23): 4381-4399.

Bai C-Y, Greenhalgh S. 2005. 3-D non-linear travel-time tomography: Imaging high contrast velocity anomalies.PureAppl.Geophys., 162(1): 2029-2049.

Bai C Y, Huang G J, Zhao R. 2010. 2-D/3-D irregular shortest-path ray tracing for multiple arrivals and its applications.Geophys.J.Int., 183(3): 1596-1612.

Bai C Y, Huang G J, Li Z S. 2011. Simultaneous inversion combining multiple-phase travel times within 3D complex layered media.ChineseJ.Geophys. (in Chinese), 54(1): 182-192, doi: 10.3969/j.issn.0001-5733.2011.01.019.

Bijwaard H, Spakman W. 1999. Fast kinematic ray tracing of first- and later-arriving global seismic phase.Geophys.J.Int., 139(2): 359-369.

De Kool M, Rawlinson N, Sambridge M. 2006. A practical grid-based method for tracking multiple refraction and reflection phases in three-dimensional heterogeneous media.Geophys.J.Int., 167(1): 253-270.

Grand S P, Van der Hilst R D, Widiyantoro S. 1997. Global seismic tomography: A snapshot of convection in the Earth.GSAToday, 7(4): 1-7.

Huang G J, Bai C Y. 2010. Simultaneous inversion with multiple traveltimes within 2-D complex layered media.ChineseJ.Geophys. (in Chinese), 53(12): 2972-2981, doi: 10.3969/j.issn.0001-5733.2010.12.021.

Huang G J, Bai C Y, Greenhalgh S. 2013. Fast and accurate global multiphase arrival tracking: the irregular shortest-path method in a 3-D spherical earth model.Geophys.J.Int., 194(3): 1878-1892.

Huang G J, Bai C Y. 2013. Simultaneous inversion of three model parameters with multiple classes of arrival times.ChineseJ.Geophys. (in Chinese), 56(12): 4215-4225, doi: 10.6038/cjg20131224.

Kennett B L N, Sambridge M S, Williamson P R. 1988. Subspace methods for large inverse problems with multiple parameter classes.Geophys.J.Int., 94(2): 237-247.

Kennett B L N, Engdahl E R, Buland R. 1995. Constraints on seismic velocities in the earth from traveltimes.Geophys.J.Int., 122(1): 108-124.

Moser T J. 1991. Shortest path calculation of seismic rays.Geophysics, 56(1): 59-67.

Preston L A. 2003. Simultaneous inversion of 3D velocity structure, hypocenter locations, and reflector geometry in Cascadia [Ph. D. thesis]. Washington: University of Washington.

Preston L A, Creager K C, Grosson R S, et al. 2003. Intraslab earthquakes: Dehydration of the Cascadia slab.Science, 302(5648): 1197-1200.

Rawlinson N, Sambridge M. 2004. Multiple reflection and transmission phases in complex layered media using a multistage fast marching method.Geophysics, 69: 1338-1350.

Rawlison N, Sambridge M S. 2003. Seismic traveltime tomography of the crust and lithosphere.AdvancesinGeophysics, 46: 81-198. Sambridge M S. 1990. Non-linear arrival time inversion: Constraining velocity anomalies by seeking smooth models in 3-D.Geophys.J.Int., 102(3): 653-677. Sambridge M S, Kennett B L N. 1990. Boundary value ray tracing in a heterogeneous medium: A simple and versatile algorithm.Geophys.J.Int., 101(1): 157-168.

Snoke J A, Lahr J C. 2001. Locating earthquakes: At what distance can the earth no longer be treated as flat?Seism.Res.Lett., 72(5): 538-541.

Su W J, Woodward R L, Dziewonski A M. 1994. Degree 12 model of shear velocity heterogeneity in the mantle.J.Geophys.Res., 99(B4): 6945-6980.

Tang X P, Bai C Y. 2009a. Multiple ray tracing in 2D layered media with the shortest path method.ProgressinGeophysics(in Chinese), 24(6): 2087-2096, doi: 10.3969/j.issn.1004-2903.2009.06.022.

Tang X P, Bai C Y. 2009b. Multiple ray tracing within 3-D layered media with the shortest path method.ChineseJ.Geophys. (in Chinese), 52(10): 2635-2643, doi: 10.3969/j.issn.0001-5733.2009.10.024.

Tian Y, Huang S H, Nolet G, et al. 2007. Dynamic ray tracing and traveltime corrections for global seismic tomography.J.Comp.Phys., 226(1): 672-687.

Van der Hilst R D, Widiyantoro S, Engdahl E R. 1997. Evidence for deep mantle circulation from global tomography.Nature, 386(6625): 578-584.

Vidale J E. 1998. Finite-difference calculation of travel times.Bull.Seism.Soc.Am., 78(6): 2062-2076.

Vidale J E. 1990. Finite-difference calculation of traveltimes in three dimensions.Geophysics, 55(5): 521-526.

Zhao D P, Hasegawa A, Horiuchi S. 1992. Tomographic imaging of P and S wave velocity structure beneath northeastern Japan.J.Geophys.Res., 97(B13): 19909-19928.

Zhao D P, Lei J S. 2004. Seismic ray path variations in a 3D global velocity model.Phys.EarthPlanet.Inter., 141(3): 153-166.

Zhao D P, Todo S, Lei J S. 2005. Local earthquake reflection tomography of the Landers aftershock area.EarthPlanet.Sci.Lett., 235(3-4): 623-631.

Zhao R, Bai C Y. 2010. Fast multiple ray tracing within complex layered media: The shortest path method based on irregular grid cells.ActaSeismologicaSinica(in Chinese), 42(4): 433-444.

附中文參考文獻(xiàn)

白超英, 黃國嬌, 李忠生. 2011. 三維復(fù)雜層狀介質(zhì)中多震相走時聯(lián)合反演成像. 地球物理學(xué)報, 54(1): 182-192, doi: 10.3969/j.issn.0001-5733.2011.01.019.

黃國嬌, 白超英. 2010. 二維復(fù)雜層狀介質(zhì)中地震多波走時聯(lián)合反演成像. 地球物理學(xué)報, 53(12): 2972-2981, doi: 10.3969/j.issn.0001-5733.2010.12.021.

黃國嬌, 白超英. 2013. 多震相走時聯(lián)合三參數(shù)同時反演成像. 地球物理學(xué)報, 56(12): 4215-4225, doi: 10.6038/cjg20131224.

唐小平, 白超英. 2009a. 最短路徑算法下二維層狀介質(zhì)中多次波追蹤. 地球物理學(xué)進(jìn)展, 24(6): 2087-2096, doi: 10.3969/j.issn.1004-2903.2009.06.022.

唐小平, 白超英. 2009b. 最短路徑算法下三維層狀介質(zhì)中多次波追蹤. 地球物理學(xué)報, 52(10): 2635-2643, doi: 10.3969/j.issn.0001-5733.2009.10.024.

趙瑞, 白超英. 2010. 復(fù)雜層狀模型中多次波快速追蹤: 一種基于非規(guī)則網(wǎng)格的最短路徑算法. 地震學(xué)報, 42(4): 433-444.

(本文編輯 何燕)

Simultaneous inversion of three model parameters with multiple phases of arrival times in spherical coordinates

HUANG Guo-Jiao1, BAI Chao-Ying2, QIAN Wei1

1DepartmentofGeologyScience&Engineering,SchoolofEarthSciencesandEngineering,HohaiUniversity,Nanjing211100,China2DepartmentofGeophysics,CollegeofGeologyEngineeringandGeomatics,Chang′anUniversity,Xi′an710054,China

To guarantee computational precision and improve the resolution of seismic traveltime tomography at a regional or global scale, we conduct a 3D simultaneous inversion of three different parameters (velocity, hypocenter location and reflector geometry) in spherical coordinates with multi-phase travel times. For this task, there are two key issues: (1) It needs an efficient and accurate arrival tracking algorithm for multiply transmitted, reflected (or refracted) and converted waves in a 3D variable velocity model with embedded velocity discontinuities (or subsurface interfaces), and (2) A subdimensional inversion solver is required which can easily search for different types of model parameters to balance the tradeoff between the different types of model parameter updated in the simultaneous inversion process.For these purposes, we first extend a popular grid/cell-based wavefront expanding ray tracing algorithm (the multistage irregular shortest-path ray tracing method, shorted as ISPM), which previously worked only in Cartesian coordinates at a local scale, to spherical coordinates appropriate to a regional or global scale. In 3D spherical coordinates a trapezoidal prism cell is used to divide the spherical earth model, except for the global and other irregular subsurface interfaces, where a trapezoidal cone is used. And in one previous paper we have discussed the computational accuracy of the multistage ISPM algorithm against the analytic solution of AK135 Travel Time Table for 49 kinds of global phases and concluded that the computational accuracy can be tuned to within the 0.1 s absolute time error when it works in the spherical coordinates. We then incorporate the subspace method to formulate a simultaneous inversion algorithm, in which the multiple classes of arrivals (including direct and reflected arrivals from different velocity discontinuities) can be used to simultaneously update the velocity fields, hypocenter location and reflector geometries.In order to illustrate the performance of inversion for three model class parameters, we have selected a regional model at a scale of 20°×20°×700 km. In the numerical experiment, different phases are used. For comparison, we design three different schemes. The first one is to simultaneously update both the velocity field and the source locations (referred as V+S inversion), the second is to simultaneously invert both the velocity field and the reflector positions (referred as V+R inversion), while the third scheme updates all three model parameter classes simultaneously (referred as V+R+S inversion). To account for picking uncertainty, random noise was added to the multiple sets of arrivals according to the expected picking errors. The inversion results show that the reconstructed velocity fields exhibit similar anomalous structures to the true field, but only partially being recovered. Interestingly, the recovered velocity fields are nearly the same regardless whether the two parameter class inversion or three parameter class inversion is attempted. The updated two subsurface interfaces are almost identical, no matter whether the V+R or the V+R+S simultaneous inversion algorithm is used. For the hypocenter location, the V+S inversion has a slightly better convergence to the real source locations than the V+R+S inversion, due to only two classes of the model parameters being updated in the simultaneous inversion process, but the difference is not so significant.The above results verify that the V+R+S simultaneous inversion scheme is feasible and applicable due to the almost equal capability with the constrained two model parameter class inversions (V+S or V+R). All three scheme algorithms are quite tolerant to modest errors in the travel time picks, which makes them robust enough for practical applications. The results show that the V+R+S simultaneous inversion method is an efficient and feasible scheme for updating the velocity field, reflector geometry and hypocenter location.

Spherical coordinates; Multistage irregular shortest-path method; Multi-phase travel times; Simultaneous inversion of three model parameters; Subspace method

10.6038/cjg20151016.

Huang G J, Bai C Y, Qian W. 2015. Simultaneous inversion of three model parameters with multiple phases of arrival times in spherical coordinates.ChineseJ.Geophys. (in Chinese),58(10):3627-3638,doi:10.6038/cjg20151016.

江蘇省自然科學(xué)基金項目(BK20150799)、國家自然科學(xué)基金項目(41504038)、河海大學(xué)中央高?；究蒲袠I(yè)務(wù)費項目(2014B13814)、教育部博士學(xué)科專項基金項目(20110205110010)資助.

黃國嬌，女，1987年4月生于湖北省隨州市，主要從事地震正、反演方法及應(yīng)用研究.E-mail: jiaojiao.1986520@163.com

10.6038/cjg20151016

P315

2014-12-22，2015-06-16收修定稿

≤≥? ?黃國嬌， 白超英， 錢衛(wèi). 2015. 球坐標(biāo)系下多震相走時三參數(shù)同時反演成像.地球物理學(xué)報，58(10):3627-3638,