球坐標系下多震相走時三參數同時反演成像

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

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

?

球坐標系下多震相走時三參數同時反演成像

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

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

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

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

1 引言

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

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

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

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

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

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

2 多震相射線追蹤算法

2.1 參數化

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

2.2 速度插值

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

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

(1)

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

(2)

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

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

2.4 計算精度和效率分析

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

圖2 球坐標系下多震相射線路徑圖淺灰色線：初至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.圖中括號內數字為當前模型參數化下節(jié)點總數.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中可以看出，當單元大小為1°×1°×50 km，加入11個次級節(jié)點(地表次級節(jié)點間隔約為8 km)時，計算這三種震相幾乎要花費近8個小時.但是，通常加入5～7個次級節(jié)點(次級節(jié)點間距為17～22.5 km)就能滿足精度需求，此時CPU時間就減少至1.0～1.5個小時.值得注意的是，反演成像中每次迭代基本上有95%的時間花費在計算這三種震相的走時和射線路徑上.因此，我們可大致估計采用上述三種震相同時反演三參數的時間.

3 反演算法

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

3.1 子空間算法(Subspace)

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

(3)

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

(4)

(5a)

(5b)

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

(6)

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

?

第三步：采用奇異值分解法(SVD)將第一步和第二步中構建的基向量{aj}標準正交化.

3.2 Frechet偏導數計算

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

(7)

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

(8)

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

在速度、界面和震源位置同時反演時，Frechet偏導矩陣包含了三部分：一是走時關于速度變化的導數，二是走時關于反射點深度變化的導數，三是走時關于震源參數的導數.

(9)

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

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

(10)

(11)

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

(12)

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

(13)

而走時關于震源深度變化的偏導數的計算可參見圖5所示.假設震源的深度擾動值為Δrh，則走時關于震源深度變化的偏導數可表示為：

(14)

圖4 走時關于反射點深度偏導數計算示意圖Fig.4 Diagrammatic showing the calculation of traveltime derivative with respect to reflector depth change

圖5 走時關于震源位置變化偏導數計算示意圖Fig.5 Diagrammatic showing the calculation of traveltime derivative with respect to source location change

同理可得走時關于震源位置其余兩個方向(Lat,Long)變化的偏導數表達式：

(15)

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

4 數值模擬分析

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

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

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

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

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

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

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

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

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

5 結論

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

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.

附中文參考文獻

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

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

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

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

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

趙瑞, 白超英. 2010. 復雜層狀模型中多次波快速追蹤: 一種基于非規(guī)則網格的最短路徑算法. 地震學報, 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.

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

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

10.6038/cjg20151016

P315

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

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