姚華建，羅 松，李 成，胡少乾，方洪健

1 中國科學技術(shù)大學地球和空間科學學院，合肥 230026

2 蒙城地球物理國家野外科學觀測研究站，蒙城 233500

3 安徽理工大學深部煤礦采動響應(yīng)與災害防控國家重點實驗室，淮南 232001

4 中國地質(zhì)大學（武漢）地球物理與空間信息學院，武漢 430074

5 中山大學地球科學與工程學院，珠海 519080

0 引言

面波是沿地球表面?zhèn)鞑サ囊环N地震波，與在地球內(nèi)部傳播的體波相對應(yīng)，是研究地殼上地幔速度結(jié)構(gòu)和地震震源機制的一種非常重要的地震波.根據(jù)偏振方向的差異，面波可以主要分為瑞利（Rayleigh）面波和勒夫（Love）面波.瑞利面波的質(zhì)點運動軌跡呈現(xiàn)為橢圓形，且垂向和徑向的瑞利面波相位相差90°；勒夫面波在與傳播方向垂直的水平面上運動.瑞利面波傳播的速度主要對地下不同深度介質(zhì)沿垂直方向偏振的SV 波的速度（VSV）比較敏感，對縱波速度和密度有較弱的敏感性.勒夫面波傳播的速度主要對地下不同深度沿水平方向偏振的SH 波的速度（VSH）敏感，對密度具有較弱的敏感性.不同周期面波的敏感深度不一樣，通常周期越長的面波穿透深度越深；因此相對于短周期面波，長周期面波對深部介質(zhì)具有更強的敏感性.且由于不同深度介質(zhì)的速度存在差異，因此導致不同周期的面波具有不同的傳播速度，這就是面波的頻散.通過提取面波的頻散信息，我們可以反演地殼上地幔的橫波速度結(jié)構(gòu).

時頻分析方法（例如,Dziewonski et al.,1969;Landisman et al.,1969;Levshin and Ritzwoller,2001;Yao et al.,2005）被廣泛應(yīng)用于測量事件到臺站或者雙臺之間的面波頻散曲線，包括群速度（U）頻散和相速度（C）頻散.波形反演方法也常被用于獲得事件到臺站之間的相速度頻散曲線（例如,Trampert and Woodhouse,1995;Van Heijst and Woodhouse,1999;Woodhouse and Dziewonski,1984），但需要考慮地震時間的初始相位問題.基于地震事件得到的頻散曲線一般周期在10 s 之上，長周期可達300 s，但10 s 以內(nèi)的短周期頻散曲線難以提取，主要是由于短周期面波受復雜上地殼結(jié)構(gòu)的散射及較強的衰減影響較大.因此，基于地震事件提取的面波頻散通常難以獲得上地殼或中上地殼的精細結(jié)構(gòu).本世紀初，一系列研究發(fā)現(xiàn)連續(xù)背景噪聲互相關(guān)方法可以有效地恢復臺站之間的中短周期面波信號（例如,Sabra et al.,2005;Shapiro and Campillo,2004;Shapiro et al.,2005;Yao et al.,2006），從而大大促進了中短周期面波頻散曲線的可靠提取，顯著提升了中上地殼（甚至地殼淺層）結(jié)構(gòu)成像的可靠性和分辨率.

1950 年代起，單臺法（或雙臺法）已經(jīng)被廣泛用于測量地震事件到臺站（或雙臺之間）的面波頻散曲線，并被用于反演路徑上的或不同區(qū)域平均的地殼上地幔層狀橫波速度結(jié)構(gòu)模型（例如,Brune and Dorman,1963;Knopoff et al.,1966;Press,1957）.從1980 年代起，基于面波頻散數(shù)據(jù)地震學家發(fā)展了面波層析成像方法（例如,Barmin et al.,2001;Montagner and Jobert,1988;Montagner and Nataf,1986,1988;Ritzwoller and Levshin,1998;Tanimoto,1986a,1986b;Trampert and Woodhouse,1995），即兩步法面波成像：第一步先通過面波走時成像反演獲得不同周期的二維相速度（或群速度）分布圖，第二步是基于每個格點的純路徑頻散反演格點下方一維層狀橫波速度模型，然后將所有格點的一維模型拼合起來得到研究區(qū)域的三維橫波速度模型.基于地震面波提取的中長周期面波頻散和背景噪聲面波提取的中短周期面波頻散可以有效地融合，得到從短周期到長周期的面波頻散資料，從而有效約束地殼上地幔的三維橫波速度結(jié)構(gòu)（例如,Yang et al.,2008;Yao et al.,2006,2008）.目前無論是三維各向同性橫波速度結(jié)構(gòu)的面波成像，還是方位各向異性橫波速度結(jié)構(gòu)的面波成像（例如,Montagner and Jobert,1988;Yao et al.,2010），以及徑向各向異性橫波速度結(jié)構(gòu)的面波成像（例如,Huang et al.,2010;Luo et al.,2013;Shapiro et al.,2005），主要采用兩步法模式.

面波兩步法成像中的第一步，即面波二維相速度或群速度成像，通?；诿娌ㄉ渚€理論，且絕大部分研究直接假設(shè)面波沿大圓路徑傳播，即不同周期面波均沿著從事件到臺站之間的最短大圓路徑傳播，忽略路徑上的結(jié)構(gòu)異常對面波傳播路徑的影響.當面波二維相速度分布橫向變化不大時（例如低于10%），基于大圓路徑假設(shè)與基于面波射線路徑路徑追蹤方法得到的二維相速度成像結(jié)果幾乎相同，但當相速度橫向變化較大時（例如高于20%），通常會造成基于大圓路徑假設(shè)得到的二維相速度成像結(jié)果存在一定的偏差（李想等，2015）.因為短周期面波對地殼淺層結(jié)構(gòu)非常敏感，當研究區(qū)域淺層結(jié)構(gòu)橫向變化很大時，短周期二維相（群）速度成像需要采用基于面波射線路徑追蹤（Rawlinson and Sambridge,2004）的走時成像方法（Young et al.,2011），即需要考慮到射線路徑彎曲對面波走時計算的影響.

當面波頻率較低（或波長較長），以及異常體的橫向尺度與面波波長相當時，通常需要考慮面波傳播的有限頻效應(yīng)，才能夠更為準確地反演結(jié)構(gòu)異常.一些研究基于散射理論提出面波走時擾動計算的二維有限頻敏感核（sensitivity kernel）方法（例如,Spetzler et al.,2002;Yoshizawa and Kennett,2005），相對于基于射線路徑的面波走時擾動計算方法可以更好地預測面波走時異常，目前已經(jīng)被應(yīng)用到區(qū)域和全球范圍內(nèi)的二維面波相速度成像研究中（例如,Ritzwoller et al.,2002;Yang and Forsyth,2006;Yoshizawa and Kennett,2004）.值得注意的是，二維敏感核計算時需要計算面波走時，而面波走時通常是基于大圓路徑計算得到的.但當介質(zhì)存在較強的橫向非均勻性時，除了有限頻效應(yīng)外，由于射線路徑彎曲所導致的走時異常是不能被忽略的（Yang and Hung,2005）.

當臺站分布較為密集且均勻時，通過事件（或背景噪聲虛擬源）到臺站的面波數(shù)據(jù)和臺陣分析方法可以有效構(gòu)建不同周期的二維相速度分布圖，例如采用基于面波走時的程函方程成像方法（Lin et al.,2009），基于面波走時和振幅的亥姆霍茲方程成像方法（Helmholtz tomography,Lin and Ritzwoller,2011），以及面波波形梯度成像方法（wave gradiometry,Liang and Langston,2009）.程函方程成像方法通過計算走時的梯度場可以獲得面波的傳播方向和路徑，即該方法考慮了面波傳播路徑的彎曲效應(yīng).由于亥姆霍茲方程成像和面波波形梯度成像方法還使用了面波振幅或波形信息，因此可以同時考慮面波的射線路徑彎曲效應(yīng)和有限頻效應(yīng)，在理論上具有更高的成像精度.

無論是基于面波走時的成像反演或臺陣面波成像方法，構(gòu)建了不同周期的二維相速度（或群速度）分布圖后，都需要進一步得到反演區(qū)域各網(wǎng)格點的純路徑頻散，然后通常采用單點頻散曲線反演的方法獲得格點下方的一維層狀橫波速度模型.由于每個格點的層狀橫波速度反演是獨立的，格點與相鄰格點之間的模型并沒有橫向平滑約束，可能會導致最后拼合得到的三維模型出現(xiàn)一些橫向假異常.此外，當有的周期路徑分布較稀疏時，往往難以得到該周期可靠的相速度（或群速度）分布圖，這些數(shù)據(jù)便無法用于反演三維橫波速度模型.

為了避免面波兩步法成像存在的問題，一些學者提出面波一步法成像，即基于測量得到的所有路徑面波頻散走時數(shù)據(jù)直接反演三維橫波速度模型的面波直接成像方法（Boschi and Ekstr?m,2002;Fang et al.,2015;Feng and An,2010）.Boschi 和Ekstr?m（2002）以及 Feng 和An（2010）提出的方法并不考慮面波路徑的彎曲效應(yīng)，在反演過程中面波始終沿大圓路徑傳播，因此對橫向非均勻性很強的速度結(jié)構(gòu)反演會存在一定的誤差，例如近地表或地殼起伏變化很大的地區(qū).而Fang 等（2015）提出的面波直接成像方法（direct surface wave tomography,DSurfTomo）考慮了不同周期面波路徑的彎曲效應(yīng)，反演是一個逐步迭代的過程，每次迭代模型更新后會重新計算不同周期的面波射線路徑及面波走時的分辨核矩陣，從而使得反演結(jié)果更為可靠和準確.基于Fang 等（2015）發(fā)展的各向同性面波直接成像方法，Liu 等（2019）進一步發(fā)展了基于瑞利面波頻散走時直接反演三維橫波速度方位各向異性模型的成像方法（DAzimSurfTomo），Hu 等（2020）進一步提出了基于瑞利面波和勒夫面波頻散走時數(shù)據(jù)直接反演三維橫波速度徑向各向異性模型的成像方法（DRadiSurfTomo）.這些新的面波直接成像方法可以避免構(gòu)建二維相速度或群速度分布圖，且同時對不同周期的面波路徑分別進行射線追蹤，避免了常規(guī)面波走時成像的大圓路徑傳播假設(shè)，更好地考慮了結(jié)構(gòu)異常對射線路徑及最終模型反演的影響，可以獲得更為可靠的三維各向同性及各向異性的橫波速度結(jié)構(gòu)模型.

本綜述論文將首先闡述面波兩步法成像的基本原理，然后詳細介紹面波直接成像方法，之后再通過幾個具體實例來展示面波直接成像方法在不同尺度結(jié)構(gòu)成像中的應(yīng)用，最后對面波直接成像方法的若干問題和發(fā)展方向進行討論與展望.

1 面波兩步法成像基本原理

在傳統(tǒng)的面波成像中我們通常假設(shè)面波沿著大圓路徑傳播，然而在強橫向非均勻介質(zhì)中，面波路徑會偏離震源和接收臺站之間的大圓路徑.在高頻近似的情況下，角頻率為 ω的面波走時可以表示為：

其中S(l,ω)是震源和接收臺站AB 之間的實際射線路徑的局部慢度.公式（1）可以表達為離散形式：

其中Sp為沿AB 的射線路徑段 ΔlAB上的相位慢度，p是路徑段的編號（如圖1a 所示）.公式（2）可以進一步表示為：

圖1 面波成像中離散化的三維網(wǎng)格點模型及面波射線路徑示意圖.（a）水平向的二維慢度網(wǎng)格點；（b）垂直向的VS 網(wǎng)格點（黑點）和插值得到的層狀模型（黑色階梯狀實線）.圖（a）中的黑色實線表示某個周期面波在AB 間的傳播路徑，路徑上的p 點的慢度由周圍4 個水平網(wǎng)格點（1、2、3、4）的雙線性值插值獲得.圖（b）中垂直格點模型通過擾動（如紅色圓點和虛線所示）來計算頻散相對于參數(shù)模型的深度敏感核（修改自Fang et al.,2015）Fig.1 Discretized 3D grid model of surface wave tomography and illustration of surface wave ray path.(a) 2D slowness grids in the horizontal direction;(b) VS grids (black dots) in the vertical direction and the interpolated layered model (black staircase lines).In (a) the black solid line represents the propagation path between two stations A and B for the surface wave at some period.The phase slowness at any point p along the path is determined from the values at four surrounding horizontal grid points (1,2,3,4) using a bilinear interpolation method.In (b) the vertical grid model is perturbed (as shown by the red dots and red dashed lines) to compute the depth sensitivity kernel of dispersion data to model parameters (modified from Fang et al.,2015)

基于測量得到的面波相速度或群速度曲線，我們可以獲得相對應(yīng)的頻散走時tAB(ω)，然后采用基于連續(xù)模型的非線性反演方法（Montagner and Nataf,1986;Torantola and Valette,1982）來反演各個周期的二維相速度或群速度分布圖.反演的目標函數(shù)定義為：

其中d0是測量得到的同一頻率、不同路徑的面波走時向量，d是通過d=g(m)公式[例如（2）式]與給定的研究區(qū)域二維相速度（或慢度）模型m計算得到的理論走時向量，m0是初始相速度（或慢度）模型向量；Cd0是描述數(shù)據(jù)不確定度的數(shù)據(jù)協(xié)方差矩陣，通常是一個對角矩陣，對角線元素為（第i條路徑走時測量誤差）的平方；Cm0是模型的先驗協(xié)方差矩陣，通常取為高斯型的協(xié)方差矩陣，使得反演得到的二維相速度或群速度模型在空間上達到一定橫向尺度的平滑.這里模型m可以僅為各向同性相速度（或群速度）參數(shù)，也可以包含方位各向異性參數(shù)（例如,Barmin et al.,2001;Montagner and Nataf,1986）.

除了采用基于連續(xù)模型的非線性反演方法來求解（4）式外，Backus-Gilbert 方法（Yanovskaya and Ditmar,1990）、基函數(shù)展開（Trampert and Woodhouse,1995）等反演方法也常用于反演二維相速度（或群速度）分布圖，這里就不一一贅述.通過不同周期反演得到的二維相速度或群速度分布圖，我們可以獲得反演區(qū)域每個格點的純路徑頻散曲線，然后基于純路徑頻散曲線進一步反演格點下方的層狀速度結(jié)構(gòu)模型，最后將所有二維面狀格點下方的一維模型拼合成為三維速度結(jié)構(gòu)模型.從單點頻散反演一維層狀速度模型通常采用非線性問題的線性迭代反演方法（例如,Herrmann,2013）或者全局優(yōu)化反演方法（例如,Shapiro and Ritzwoller,2002;Yao et al.,2008）.由于瑞利面波頻散主要對垂直方向偏振的橫波速度VSV敏感，對密度ρ和縱波速度VP的敏感度不高，通常在反演時采用經(jīng)驗關(guān)系式將 ρ和VP與VSV耦合在一起，從而僅反演層狀VSV結(jié)構(gòu)（例如,Yao et al.,2008）.

勒夫面波主要對水平方向極化的橫波速度VSH敏感，對 ρ的敏感性較小，因此可以用于反演層狀VSH速度結(jié)構(gòu)模型.在同時具有瑞利面波和勒夫面波頻散數(shù)據(jù)時，我們可以同時或分別反演獲得VSV和VSH速度結(jié)構(gòu)（例如,Shapiro and Campillo,2004;Huang et al.,2010），再進一步獲得介質(zhì)的徑向各向異性參數(shù)，例如γ=VSH/VSV，或者ξ=2(VSHVSV)/(VSH+VSV).

對于面波方位各向異性的兩步法成像，通常采用基于連續(xù)模型的非線性反演方法（Montagner and Nataf,1986;Torantola and Valette,1982）同時反演各個周期的各向同性相速度和方位各向異性參數(shù)，之后再反演每個格點上的純路徑方位各向異性頻散曲線，獲得格點下方的橫波速度及其方位各向異性參數(shù)（Montagner and Nataf,1986）.Smith 和Dahlen（1973）通過研究發(fā)現(xiàn)弱彈性各向異性介質(zhì)可以引起瑞利波和勒夫波相速度的方向依賴性，并在不考慮 4Ψ項（Ψ 為面波傳播方位角）的影響下，位置在M處、角頻率為 ω的瑞利面波相速度c(ω,M,Ψ)可表示為如下的近似表達式：

其中c0(ω)為參考相速度，α0(ω)是各向同性相速度擾動，α1(ω)和α2(ω)是方位各向異性參數(shù)，A=分別是該頻率下瑞利面波方位各向異性的強度大小和快波方向.

根據(jù)Montagner 和Nataf（1986），M處的瑞利波相速度擾動δc(ω,M,Ψ)可以近似表示為：

公式（6）中A、C、F、L和N用來描述具有垂直對稱軸的橫向各向同性介質(zhì)（VTI）的彈性參數(shù)，其 中F=η(A-2L)，VPH和VPV分別是沿水平和垂直方向傳播的P 波速度，VSH和VSV分別是沿水平和垂直方向極化的S 波速度.另外的六個參數(shù)Bs,c、Gs,c和Hs,c分別給出了彈性參數(shù)A、L、F的 2ψ項方位變化.敏感核?c/?pi（pi分別是A、L或F）可由一維參考模型計算獲得（Montagner and Nataf,1986）.

Montagner 和Nataf（1986）發(fā)現(xiàn)?c/?L項在公式（6）中貢獻最大，?c/?A項在地殼中相對較大，而在上地幔中可忽略不計（圖2）.?c/?F項較小，因此在實際反演中常忽略?c/?F項.（6）式的反演可以基于線性反演方法（Montagner and Nataf,1986），首先基于某一參考模型計算?c/?pi，然后再反演（6）式中相關(guān)參數(shù)，主要是L和Gs,c；或者基于全局優(yōu)化的反演方法，例如Yao（2015）提出基于鄰域算法的方位各向異性橫波速度模型反演方法.

圖2 一維速度模型（a）及其相應(yīng)的d c/dL（實線）和d c/dA（虛線）在不同周期的深度敏感核（b）（修改自Liu et al.,2019）Fig.2 Depth sensitivity kernels (b) of dc/dL (solid line) and dc/dA (dashed line) at different periods from a 1D velocity model (a)(modified from Liu et al.,2019)

最后，垂向方向偏振的方位各向異性橫波波速可以表示為：

由于Gc,s通常遠遠小于L，即，因此，（7）式可以近似表示為：

其中ΛSV和?F分別是VSV的方位各向異性振幅大小和快軸的方位角：

2 面波直接成像方法

在第1 節(jié)中我們闡述了面波兩步法成像的基本原理和常用方法，本節(jié)我們將系統(tǒng)闡述面波直接成像方法（即面波一步法成像）的原理，主要基于我們近些年發(fā)展的考慮不同周期面波射線路徑追蹤的面波直接成像方法體系，包括直接反演三維各向同性橫波速度（VS）模型的成像方法（DSurfTomo,Fang et al.,2015）、直接反演三維橫波速度（VSV）及其方位各向異性模型的成像方法（DAzim-SurfTomo,Liu et al.,2019），以及直接反演三維橫波速度結(jié)構(gòu)及徑向各向異性模型的成像方法（DRadiSurfTomo,Hu et al.,2020）.三種面波直接成像方法的簡介如圖3 所示.

圖3 面波直接成像方法（DSurfTomo、DAzimSurfTomo和DRadiSurfTomo）的簡介圖Fig.3 The sketch diagrams of the direct surface wave tomography methods (DSurfTomo,DAzimSurfTomo,and DRadi-SurfTomo)

2.1 三維各向同性橫波速度模型反演的面波直接成像方法

我們將地下三維結(jié)構(gòu)離散為三維網(wǎng)格點，包括水平向的網(wǎng)格點（圖1a）和垂向的網(wǎng)格點（圖1b），垂向的層狀速度模型可以通過垂向網(wǎng)格點插值來獲得.源A 和臺站B 之間的面波走時（頻率為 ω）可以（3）式來表示，這里考慮由于二維面波相速度的橫向非均勻性導致的面波傳播路徑的彎曲效應(yīng)，即偏離傳統(tǒng)假設(shè)的大圓路徑（圖1a）.

對于一維層狀結(jié)構(gòu)模型（圖1b），相速度的擾動可以表示為：

其中 θk表示一維層狀參考模型，αk(z)、βk(z) 和ρk(z)分別是縱波速度、橫波速度和密度.根據(jù)Brocher（2005）提出的經(jīng)驗公式，我們可以將縱波速度和密度用橫波速度的來表示.基于（3）式，觀測面波到時與理論面波到時的殘差可以表示為：

將（12）式帶入（11）式，同時利用Brocher（2005）經(jīng)驗公式，我們得到：

我們構(gòu)建的反演問題的目標函數(shù)為：

其中，右邊第一項為數(shù)據(jù)殘差項（基于L2 范數(shù)），第二項為模型正則化項（基于L2 范數(shù)），L是模型平滑算子，這里我們通常選用一階Tikhonov 模型正則化約束，λ是用來平衡數(shù)據(jù)擬合和模型平滑的權(quán)重參數(shù)，通常采用L-curve 方法來確定（Aster et al.,2013）.公式（14）所示的反演問題可以通過LSQR 算法（Paige and Saunders,1982）求解.基于（14）式，我們可以反演求解得到三維網(wǎng)格點的橫波速度模型的擾動量，之后得到新的三維橫波速度模型，縱波速度和密度再通過與橫波速度的經(jīng)驗公式（Brocher,2005）來計算獲得.之后，我們會基于更新的三維模型重新計算不同周期的相速度分布圖，然后采用快速行進方法（Rawlinson and Sambridge,2004）來計算不同周期面波的射線路徑，之后再更新G矩陣以及反演得到下一次模型的更新量，如此迭代直到反演結(jié)果收斂.

除了常規(guī)的基于Tikhonov 模型約束的反演方法外，F(xiàn)ang 等（2015）還提出了一種基于小波基函數(shù)稀疏反演的策略來求解（14）式，可以實現(xiàn)多尺度的模型分辨.Fang 等（2015）的反演方法基于預先設(shè)定的一套三維網(wǎng)格點，Luo 等（2021）在此基礎(chǔ)上進一步提出基于多套交錯網(wǎng)格的面波直接成像方法，通過不同網(wǎng)格反演得到的模型可以獲得最終平均模型以及模型的誤差估計，測試結(jié)果也顯示基于多套交錯網(wǎng)格的反演方法有助于提升模型的分辨率.

在傳統(tǒng)的面波兩步法成像中需要構(gòu)建不同周期的二維相速度和群速度分布圖，但在有些情況下由于某些周期的路徑較少，難以獲得可靠的二維相（群）速度分布圖，因此這些周期的數(shù)據(jù)無法用于兩步法成像中.但在直接反演方法中，所有周期的頻散走時測量數(shù)據(jù)均可以用于成像中，提升了數(shù)據(jù)的利用率和路徑覆蓋密度，從而有助于提升成像的分辨率.且由于每個周期每條路徑的面波射線路徑都分別計算獲得，可以更好地考慮復雜介質(zhì)（尤其是地殼淺層）對不同周期面波走時測量的影響，有助于提升反演結(jié)果的可靠性.

雖然面波頻散對縱波速度和密度的敏感性要明顯低于對橫波速度的敏感性，但直接忽略縱波速度和密度的影響也會影響橫波速度的反演結(jié)果，尤其是瑞利面波頻散對淺層地殼縱波速度敏感性也比較大.因此，面波直接成像方法在頻散走時敏感核矩陣的構(gòu)建時考慮到了縱波速度和密度對面波頻散的影響，通過經(jīng)驗公式將縱波速度和密度與橫波速度耦合在一起，從而將縱波速度和密度對頻散數(shù)據(jù)的敏感核轉(zhuǎn)變?yōu)閷M波速度的敏感核[公式（13）]，進一步提升了橫波速度的敏感性，有助于更好地反演橫波速度結(jié)構(gòu).

基于背景噪聲方法和地震事件獲得的面波頻散資料，該面波直接成像方法目前已經(jīng)被廣泛應(yīng)用于不同尺度的三維橫波速度結(jié)構(gòu)成像研究中，為認識不同區(qū)域的地殼上地幔結(jié)構(gòu)以及地殼淺層精細結(jié)構(gòu)提供了穩(wěn)定可靠的成像方法，我們將在第3 節(jié)中做簡要介紹.

2.2 三維方位各向異性橫波速度模型反演的面波直接成像方法

三維方位各向異性橫波速度模型的直接成像方法（DAzimSurfTomo,Liu et al.,2019）是在三維各向同性橫波速度直接成像方法（Fang et al.,2015）的基礎(chǔ)上進一步發(fā)展得到的.我們在有K個水平格點（k=1,2,...,K）的規(guī)則網(wǎng)格的二維球面上對研究區(qū)域進行參數(shù)化.我們以每個網(wǎng)格點k上的一維橫向各向同性分層模型 Θk給出的各向同性相速度作為參考值，這時頻率為ω的二維相速度圖可以近似表示為：

采用雙線性插值，我們計算路徑上第p段的慢度：

其中 νpk是雙線性插值函數(shù)系數(shù).第k個格點上的相速度可以由一維橫向各向同性分層模型Θk正演計算獲得：

其中g(shù)(Θk,ω,ψ)是面波頻散的正演計算函數(shù)，可以使用傳播矩陣方法對具有半無限空間的分層介質(zhì)進行計算（Haskell,1953;Herrmann,2013），或采用球形介質(zhì)的自由震蕩方法（例如,Dahlen and Tromp,1999;Woodhouse,1988;Yang et al.,2010）.

根據(jù)Liu 等（2019），頻率為 ω時第i個瑞利面波走時測量可以表示為：

其中 μik是沿著第i個走時數(shù)據(jù)對應(yīng)的射線路徑上的雙線性插值系數(shù).在這里，我們假設(shè)方位各向異性相速度變化的幅度通常僅為百分之幾，一般遠小于各向同性相速度變化.因此，通過快速行進方法（Rawlinson and Sambridge,2004）計算射線路徑時僅基于二維各向同性相速度(ω)來計算頻率為ω時的射線路徑，忽略弱方位各向異性對介質(zhì)射線路徑的影響.

如果我們對方程（18）的面波走時ti進行微擾，并以各向同性相速度(ω)作為參考值，觀測值和參考各向同性模型預測值之間的走時殘差可近似為：

基于Liu 等（2019）的推導，我們最終可以得到：

等式（20）可以用d=Gm的矩陣形式來表達，其中d是所有頻率、所有路徑的瑞利面波走時的殘差向量，G是數(shù)據(jù)敏感核矩陣，模型參數(shù)向量m由下式給出：

考慮到模型參數(shù)的平滑約束，我們構(gòu)建的反演方程如下：

其中Giso和miso對應(yīng)于公式（20）中第一部分關(guān)于對VSV各向同性擾動的敏感核矩陣和橫波速度擾動向量，GAA和mAA對應(yīng)于公式（20）中第二部分關(guān)于方位各向異性參數(shù)擾動的敏感核矩陣和方位各向異性參數(shù)擾動向量，Liso和LAA分別是各向同性和方位各向異性模型參數(shù)的粗糙度矩陣，λ1和 λ2是用于平衡數(shù)據(jù)擬合和模型正則化的權(quán)重因子.Liso和LAA通常由Tikhonov 正則化的有限差分算子表示，可以將其選擇為一階或二階空間導數(shù)算子（Aster et al.,2013）.在每次迭代反演中，當各向同性的VSV（即 β）模型迭代更新后，會重新計算二維各向同性相速度分布圖以及新的射線路徑，再更新各向同性和各向異性參數(shù)的敏感核矩陣.通過最后反演得到的模型參數(shù)m，我們可以由（9-10）式計算得到橫波速度VSV的方位各向異性強度和快波方向.

考慮到各向同性橫波速度和方位各向異性參數(shù)在反演中可能存在的折衷性（trade-off），我們通常采用兩階段方式進行反演（Liu et al.,2019）：首先基于（13-14）式直接反演三維各向同性的橫波速度模型（Fang et al.,2015），然后基于得到的三維各向同性模型作為初始模型再進一步根據(jù)（20-22）式同時反演三維各向同性及方位各向異性的橫波速度模型.這種兩階段反演方式可以大大減輕各向同性橫波速度初始模型的不準確性對方位各向異性反演的影響，一些實際應(yīng)用例子證實該方法是一種有效和穩(wěn)定的反演方法.

2.3 三維徑向各向異性橫波速度模型反演的面波直接成像方法

三維徑向各向異性結(jié)構(gòu)的面波直接成像方法（DRadiSurfTomo,Hu et al.,2020）是在各向同性面波直接成像方法（DSurfTomo,Fang et al.,2015）的基礎(chǔ)上進一步發(fā)展的.該方法可以用來同時反演瑞利波和勒夫波的頻散走時數(shù)據(jù)（TR,TL）來獲取平均剪切波速度（VS）和徑向各向異性參數(shù)ξ=2(VSHVSV)/(VSH+VSV)（Hu et al.,2020）.首先，使用垂直極化剪切波速VSV和無量綱化的徑向各向異性參數(shù) γ（定義為γ=VSH/VSV）對模型進行參數(shù)化，可以唯一地將通常使用的徑向各向異性參數(shù)表達為ξ=2(γ-1)/(γ+1)（Huang et al.,2010;Xie et al.,2013）.

當反演參數(shù)為VSV和VSH時，根據(jù)面波直接反演方法（Fang et al.,2015）的（13）式，基于瑞利波和勒夫波的頻散走時數(shù)據(jù)反演三維VSV和VSH的問題可以表示為：

其中 ΔTR和 ΔTL分別是所有路徑和所有周期的瑞利波和勒夫波走時的擾動向量;GSV和GSH分別是瑞利波和勒夫波走時相對于橫波速度參數(shù)VSV和VSH的敏感核矩陣， ΔVSV和 ΔVSH分別是VSV和VSH模型參數(shù)的更新向量，在實際反演中我們會分別對 ΔVSV和ΔVSH施加空間平滑約束，通過逐步迭代來獲得最終的VSV和VSH模型結(jié)果，并計算得到徑向各向異性參數(shù)ξ.通常勒夫波信號的信噪比相對較低，因此勒夫波頻散數(shù)據(jù)要比瑞利波數(shù)據(jù)少很多，且勒夫波數(shù)據(jù)的分辨深度要相對更淺一些（圖4）.由于瑞利波和勒夫波射線路徑及敏感核的不同，可能會導致VSH反演結(jié)果在局部區(qū)域不穩(wěn)定或具有較低的分辨率，因此會導致相應(yīng)的徑向各向異性結(jié)果出現(xiàn)不穩(wěn)定性.

圖4 不同周期瑞利波和勒夫波相速度深度敏感核的對比.（a）計算深度敏感核所采用的橫波速度模型；（b）5 s、20 s、40 s 周期的瑞利面波相速度的深度敏感核；（c）5 s、20 s、40 s 周期的勒夫面波相速度的深度敏感核Fig.4 The comparison of depth sensitivity kernels of Rayleigh wave and Love wave at different periods.(a) The shear wave velocity model used to compute depth sensitivity kernels;(b) Depth sensitivity kernels of Rayleigh wave phase velocities at periods of 5 s,20 s,and 40 s;(c) Depth sensitivity kernels of Love wave phase velocities at periods of 5 s,20 s,and 40 s

當反演參數(shù)轉(zhuǎn)換為VSV和γ 時，基于γ=VSH/VSV，我們可以將ΔVSH表示成：

這里γ 和 Δγ分別為向量，算符·表示兩個向量元素之間對應(yīng)相乘，即(A·B)i=AiBi.

將（24）式帶入（23）式并進行合并，最后線性化反演問題可表達為：

與三維方位各向異性橫波速度結(jié)構(gòu)反演類似，我們需要對三維模型參數(shù) ΔVSV和 Δγ分別施加模型平滑正則化項，然后通過迭代求解方法來獲得ΔVSV和 Δγ，反演收斂時最終得到VSV和 γ，然后再計算三維平均剪切速度VS=VSV(1+γ)/2和徑向各向異性ξ.因為該方法可以直接對無量綱化的徑向各向異性參數(shù) γ施加模型正則化平滑，所以可以降低瑞利波和勒夫波路徑分布以及深度敏感核差異對反演結(jié)果的影響，使得徑向各向異性參數(shù)的反演結(jié)果更加平滑、穩(wěn)定、可靠.基于理論模型的分辨率測試結(jié)果表明，三維徑向各向異性模型的直接成像方法可以較單獨反演三維VSV和VSH模型，然后再計算得到徑向各向異性的方法具有更好的空間分辨率，尤其是在相對深部區(qū)域（Hu et al.,2020）.

3 應(yīng)用實例

在第2 節(jié)中我們介紹了基于面波頻散走時直接反演三維各向同性橫波速度、方位各向異性橫波速度以及徑向各向異性橫波速度模型的成像方法，在本節(jié)中我們將展示通過不同區(qū)域、不同頻帶的面波頻散走時資料開展不同尺度和不同深度范圍的三維橫波速度結(jié)構(gòu)直接成像的應(yīng)用實例.

3.1 三維各向同性橫波速度結(jié)構(gòu)直接成像的應(yīng)用

面波直接成像方法（Fang et al.,2015）已經(jīng)被廣泛應(yīng)用于不同尺度三維各向同性橫波速度結(jié)構(gòu)的反演，包括幾百至上千千米的區(qū)域尺度的地殼結(jié)構(gòu)（例如,臺灣海峽,Chen et al.,2016b;喜馬拉雅造山帶,Singer et al.,2017;福建—臺灣海峽地區(qū),Zhang Y Y et al.,2018;長江中下游成礦帶,Luo et al.,2019;郯廬斷裂帶中南段,Bem et al.,2020;日本中部,Nimiya et al.,2020;青藏高原,Huang et al.,2020），數(shù)十千米至數(shù)千米的斷裂帶尺度的上地殼結(jié)構(gòu)（例如,2016 年臺灣MW6.4 級梅隆地震震源區(qū),Kuo-Chen et al.,2017;郯廬斷裂帶肥東段,Gu et al.,2019;龍門山斷裂帶,Zhang Y T et al.,2020;郯廬斷裂帶廬江段,Li C et al.,2020;程海斷裂帶,Yang H F et al.,2020;郯廬斷裂帶巢湖段,Luo et al.,2021），數(shù)千米到數(shù)十千米的城市和盆地尺度的地殼淺層結(jié)構(gòu)（例如,臺北盆地,Fang et al.,2015;臺灣宜蘭盆地,Chen et al.,2016a;合肥市區(qū),Li et al.,2016;李玲利等,2020;洛陽盆地,Zhou et al.,2018），以及數(shù)千米到數(shù)十千米的資源勘探尺度近地表或地殼淺層結(jié)構(gòu)（例如,呼圖壁儲氣庫區(qū)域,王娟娟等,2018;重慶某頁巖氣開發(fā)水壓致裂區(qū),Liu et al.,2018;新疆Cu-Ni 礦集區(qū),Du et al.,2020;遼寧白云金礦區(qū),Li X et al.,2020）.除了采用背景噪聲面波數(shù)據(jù)外，中長周期的地震面波和中短周期的背景噪聲面波也聯(lián)合在一起，用于直接反演地殼上地幔三維橫波速度結(jié)構(gòu)，例如福建—臺灣海峽地區(qū)（Zhang Y Y et al.,2020）和青藏高原東南緣（Zhang Z Q et al.,2020）.

由于地殼淺部介質(zhì)的橫向變化非常大，例如盆地地區(qū)的近地表橫波速度僅為數(shù)百米每秒，但基巖出露的地區(qū)近地表橫波速度可超過3 km/s，因此在結(jié)構(gòu)復雜區(qū)域的地殼淺部結(jié)構(gòu)成像時必須考慮不同周期射線路徑的彎曲，否則僅采用大圓路徑的面波成像會出現(xiàn)較大的偏差.例如，在臺北盆地的地殼淺層噪聲面波成像研究中，由于地殼淺層速度結(jié)構(gòu)差異很大，導致射線路徑存在較大的彎曲（Fang et al.,2015），且相同臺站對、不同周期的面波射線路徑會存在明顯差異，反映出地下介質(zhì)隨深度也存在不同的變化模式.合肥市區(qū)的近地表噪聲面波成像研究中（Li et al.,2016），如果假設(shè)面波沿大圓路徑傳播，那么成像結(jié)果并不能揭示穿過市區(qū)的蜀山斷層的位置（圖5b），但如果考慮不同周期面波射線路徑的彎曲，則可以清晰揭示蜀山斷層的位置（圖5a）.因此，在面波直接成像方法中考慮不同周期面波射線路徑的彎曲效應(yīng)非常重要，這也是保障結(jié)構(gòu)復雜地區(qū)面波各向同性速度結(jié)構(gòu)可靠成像的重要基礎(chǔ)，同時也是保障后續(xù)各向異性成像結(jié)果可靠的重要條件.

圖5 考慮面波射線路徑彎曲（a）和大圓路徑假設(shè)下（b）的300 m 深度處的橫波速度反演結(jié)果.（a,b）圖中三角形表示本研究中使用的地震臺站位置，灰線表示合肥市區(qū)內(nèi)的主要道路.黑色實線表示蜀山斷層（SSF）位置.紅色直線PP'表示橫跨斷層的一條剖面的位置，其速度剖面如（c）所示.（c）中黑色三角形指示根據(jù)速度差異推測的蜀山斷層（SSF）的位置（修改自 Li et al.,2016）Fig.5 Shear wave velocity inversion results at the depth of 300 m considering the surface wave ray path bending effect (a) and with the great-circle propagation hypothesis (b).In (a,b) the triangles show the stations used in the study,gray lines show the main roads in the Hefei city,the black line shows the location of the Shushan Fault (SSF),and the red line PP' represents the profile location across the fault,with its velocity profile shown as (c).In (c) the black triangle indicates the location of the Shushan fault (SSF) (modified from Li et al.,2016)

3.2 三維方位各向異性橫波速度結(jié)構(gòu)直接成像的應(yīng)用

基于區(qū)域固定和流動寬頻帶臺站的中短周期的背景噪聲面波頻散數(shù)據(jù)，三維方位各向異性橫波速度結(jié)構(gòu)的面波直接反演方法已經(jīng)成功應(yīng)用于數(shù)百千米尺度的地殼和上地幔頂部成像研究中，包括云南地區(qū)（Liu et al.,2019）、四川盆地及其周緣地區(qū)（Zhang Z Q et al.,2022）、渤海灣地區(qū)（Li et al.,2022）、郯廬斷裂帶中南段中部地區(qū)（Bem et al.,2022）.這些研究揭示了地殼和上地幔頂部隨深度和區(qū)域變化的方位各向異性特征，通過與遠震和近震橫波分裂方法得到的方位各向異性結(jié)果對比，更好地揭示了隨深度變化的變形模式.例如在云南地區(qū)方位各向異性的成像結(jié)果（Liu et al.,2019）如圖6 所示，地殼方位各向異性的快波方向主要呈現(xiàn)為南北向，主要受南北向的斷裂帶分布和地殼物質(zhì)運移的控制；但上地幔頂部的快波方向主要呈近東西向，與遠震橫波分裂測量的快波方向類似，顯示出上地幔物質(zhì)變形與地殼物質(zhì)變形的模式顯著不同，可能主要受上地幔軟流圈物質(zhì)的東西向流動控制，這與印度板塊在緬甸下方的東向俯沖及西向后撤可能密切相關(guān).基于區(qū)域固定和流動寬頻帶臺站的中短周期背景噪聲面波頻散和中長周期地震面波頻散的聯(lián)合成像，Luo 等（2022）獲得了郯廬斷裂帶中南段及鄰近區(qū)域的地殼上地幔三維方位各向異性的橫波速度結(jié)構(gòu)，對認識郯廬斷裂帶的起源、演化和區(qū)域變形機理起到了重要的促進作用.

圖6 云南地區(qū)三維方位各向異性橫波速度結(jié)構(gòu)模型.（a-c）圖中的紅色短棒顯示方位異性的大小和快波方向；圖（c）中的藍色短線顯示遠震橫波分裂測量結(jié)果（常利軍等，2015），黑色箭頭指示板塊絕對運動方向（Argus et al.,2011）（修改自Liu et al.,2019）Fig.6 3D azimuthally anisotropic shear wave velocity model in the Yunnan area.In (a-c) the red bars give the amplitude and fast direction of azimuthal anisotropy.The blue bars show the teleseismic shear wave splitting measurements (Chang et al.,2015)and the black arrow indicates the absolute plate motion (Argus et al.,2011) (modified from Liu et al.,2019)

通常固定或流動寬頻帶臺站的間距比較大，一般在數(shù)十千米，因此難以提取短周期面波頻散數(shù)據(jù)，所以難以獲得地殼淺部或近地表的方位各向異性.而地殼淺部的方位各向異性對于認識區(qū)域最大主壓應(yīng)力方向、構(gòu)造（例如斷層）走向分布等具有重要的作用.隨著短周期密集臺陣的廣泛布設(shè)，臺站之間的間距通常在數(shù)千米，甚至達數(shù)百米，因此可以有效提取從幾赫茲到幾秒周期頻帶的短周期頻散曲線，從而可以用于反演地殼淺層的三維方位各向異性，對認識淺部的構(gòu)造和變形模式具有重要的意義.基于在郯廬斷裂帶安徽巢湖段和山東濰坊段布設(shè)的短周期密集臺陣，臺站間距約為2～5 km，可以有效提取到0.5 s 周期的短周期面波頻散曲線（靳佳琪等,2022;Luo and Yao,2021），最終反演獲得了郯廬斷裂帶不同區(qū)段的上地殼精細三維方位各向異性橫波速度模型.

圖7 展示了郯廬斷裂帶巢湖段上地殼不同深度的三維方位各向異性橫波速度分布圖，模型所揭示的異常特征與主要地質(zhì)單元具有較好的一致性.各向同性結(jié)構(gòu)從西至東整體呈現(xiàn)“交錯”模式，即合肥盆地和巢湖東支主要呈現(xiàn)明顯的低速異常，這與這些區(qū)域自白堊紀以來的松散沉積直接有關(guān)，而郯廬斷裂帶、張八嶺隆起以及銀屏山地區(qū)主要為高速異常，與這些區(qū)域廣泛出露的巖漿巖、變質(zhì)巖以及早期沉積巖有關(guān).方位各向異性結(jié)構(gòu)揭示合肥盆地快波方向整體呈現(xiàn)北西—南東方向，反映主要受控于大別造山帶節(jié)理.與合肥盆地相鄰的郯廬斷裂帶在淺地殼表現(xiàn)為弱方位各向異性，隨著深度的增加，在張八嶺隆起區(qū)域呈現(xiàn)沿斷裂帶方向（即北北東）的快波異常，可能與郯廬斷裂帶大規(guī)模左旋走滑所產(chǎn)生的構(gòu)造節(jié)理有關(guān).銀屏山地區(qū)的快波方向與山體走向平行，反映主要受控于中生代華北華南板塊匯聚擠壓所產(chǎn)生的華南褶皺系.新模型所揭示的方位各向異性速度異常特征為研究郯廬斷裂巢湖段淺部地殼變形結(jié)構(gòu)及多期次構(gòu)造演化活動提供了一個重要的觀測窗口.

圖7 郯廬斷裂帶巢湖段主要地質(zhì)單元、地震臺站分布和成像結(jié)果.（a）研究區(qū)域主要地質(zhì)單元和地震臺站分布.其中黑色三角形為地震臺站位置，五角星表示主要的城市位置，深灰色實線為主要斷裂位置，包括六安斷裂（LAF）、巢湖斷裂（CLF）、照明山斷裂（ZMSF）、滁河斷裂（CHF）、嘉山—廬江斷裂（JSLJF）、池河—太湖斷裂（CHTHF）、盛橋—柏山斷裂（SQBSF）.由西至東，主要地質(zhì)單元包括合肥盆地、張八嶺隆起、郯廬斷裂帶、巢湖、銀屏山（YPM）.（b-d）深度分別為1 km、4.5 km 和8 km 深度的各向同性和方位各向異性橫波速度分布.其中小短線的長度和方位分別指示方位各向異性的強度和方位角，灰色小短線表示該方位各向異性具有較大誤差（修改自Luo and Yao,2021）Fig.7 Distribution of main geological units,seismic stations,and tomographic results in the Chao Lake segment of the Tanlu fault zone.(a) Distribution of major geological units and seismic stations within the study area.The black triangle represents the location of the seismic station.The pentagram represents the location of the main city.The dark gray solid lines indicate the main fault locations,including the LuAn fault (LAF),the Chao Lake fault (CLF),the Zhaomingshan fault (ZMSF),the Chuhe fault (CHF),the Jiashan-Lujiang fault (JSLJF),the Chihe-Taihu fault (CHTHF),and the Shenqiao-Baishan fault (SQBSF).From west to east,the main geological units are the Hefei basin,the Zhangbaling uplift,the Tanlu fault zone,the Chao lake,and the Yinping mountain (YPM).(b-d) Slices of isotropic and azimuthally anisotropic shear wave velocity at depths of 1 km,4.5 km,and 8 km,respectively.The length and azimuth of the short line indicate the strength and azimuth of the anisotropy,respectively.The gray short line indicates that the anisotropy has a large error (modified from Luo and Yao,2021)

3.3 三維徑向各向異性橫波速度結(jié)構(gòu)直接成像的應(yīng)用

基于固定臺站和布設(shè)的流動寬頻帶臺站，三維徑向各向異性橫波速度結(jié)構(gòu)的直接成像方法已經(jīng)成功地應(yīng)用于青藏高原東構(gòu)造結(jié)南迦巴瓦地區(qū)（Hu et al.,2020）、小江斷裂帶及鄰區(qū)（Yang Y et al.,2020）以及福建和臺灣海峽地區(qū)（Zhang Y Y et al.,2022），獲得了不同研究區(qū)域地殼和上地幔頂部高分辨率的三維徑向各向異性橫波速度結(jié)構(gòu)模型，對認識區(qū)域地殼構(gòu)造變形特征起到了重要的作用.

例如，在喜馬拉雅東構(gòu)造結(jié)區(qū)域，基于流動寬頻帶地震數(shù)據(jù)提取得到了5～40 s 的面波相速度頻散數(shù)據(jù)（Hu et al.,2020）.由于水平分量的噪聲互相關(guān)數(shù)據(jù)的信噪比要低于垂直分量的噪聲互相關(guān)數(shù)據(jù)，因此最終獲得的勒夫波相速度頻散走時數(shù)據(jù)僅約為瑞利波相速度頻散走時數(shù)據(jù)的一半.利用這些面波頻散數(shù)據(jù)直接反演方得到了該區(qū)域地殼的三維橫波速度以及徑向各向異性模型（圖8），檢測板測試表明該模型的橫向分辨率在1.5°左右.

結(jié)果顯示，該區(qū)域橫波速度和徑向各向異性結(jié)構(gòu)隨深度存在明顯的變化.在5 km 深度附近，南迦巴瓦峰附近顯示為低速且負各向異性（VSV＞VSH），該負各向異性結(jié)構(gòu)從南迦巴瓦峰向東沿著東嘉犁斷裂分布.相反，邊壩—洛隆斷裂顯示為高速正各向異性.在25 km 深度附近，該研究區(qū)域整體顯示為正各向異性.然而，在該深度的速度模型上，沿著東經(jīng)95°，研究區(qū)域被分為東西兩個部分，呈現(xiàn)西低東高的特征.在35 km 深度附近，該區(qū)域可以大致分為A、B、C 三部分（見圖8e 和8f）.其中，A 區(qū)域為低速正各向異性，B 區(qū)域為高速正各向異性，C 區(qū)域為低速負各向異性.結(jié)合前人在該區(qū)域的地震及其它地球物理研究結(jié)果，新模型表明：（1）淺層5 km 附近負徑向各向異性結(jié)構(gòu)與較強地震活動性，以及較強應(yīng)力積累有較好的空間對應(yīng)關(guān)系；（2）中地殼東經(jīng)95°附近可能存在強度較大的結(jié)構(gòu)，阻斷了中下地殼軟弱物質(zhì)自西向東的運移；（3）班公—怒江縫合帶以北下地殼（C 區(qū)域）的負各向異性可能與印度板塊和亞歐板塊的擠壓碰撞有關(guān).

圖8 喜馬拉雅東構(gòu)造結(jié)5 km、25 km 和35 km 深度的平均橫波速度結(jié)構(gòu)Vs（a,c,e）及其相對應(yīng)深度的徑向各向異性結(jié)構(gòu)ξ（b,d,f）.白線和灰線代表劃分的不同子區(qū)域間的邊界.（a）中灰色三角為南迦巴瓦峰的位置，粉色三角為臺站位置；（c）和（e）中的紫線代表25 km 和35 km 深度地殼中帶狀軟弱物質(zhì)的位置.（a）和（b）中英文縮寫為：BNS：班公—怒江縫合帶；ITS：印度—雅魯藏布江縫合帶；JLF：嘉犁斷裂；BB-LLF：邊壩—洛龍斷裂；CR：Comei 裂谷；MLF：米林斷裂；MTF：墨脫斷裂；NB：南迦巴瓦（修改自 Hu et al.,2020）Fig.8 Average shear wave velocity VS (a,c,e) and corresponding radial anisotropy ξ (b,d,f) around the eastern Himalayan syntaxis shown at depths of 5 km,25 km,and 35 km,respectively.White and gray dashed lines denote the boundaries for different subregions.Pink triangles in (a) denote the stations used in this study.The gray triangle in (a) denotes the location of Namche Barwa peak.Purple lines in (c) and (e) denote the proposed channelized weak zones at depths of 25 km and 35 km,respectively.Abbreviations in (a) and (b) are as follows: BNS,Bangong-Nujiang Suture;ITS,Indus-Tsangpo Suture;JLF,Jiali Fault;BB-LLF,Bianba-Luolong Fault;CR,Comei Rift;MLF,Mainling Fault;MTF,Motuo Fault;NB,Namche Barwa (modified from Hu et al.,2020)

目前基于短周期密集臺陣的水平分量（切向）噪聲互相關(guān)函數(shù)提取勒夫波頻散曲線并進一步反演三維VSH結(jié)構(gòu)，以及聯(lián)合瑞利波頻散數(shù)據(jù)反演三維徑向各向異性結(jié)構(gòu)的研究還非常少.一是因為水平分量的噪聲互相關(guān)函數(shù)通常具有較低的信噪比，另外一個重要的原因是在復雜地殼淺部情況下，沿著臺站對大圓路徑方向旋轉(zhuǎn)得到的切向分量噪聲互相關(guān)函數(shù)有可能還包含了瑞利面波的成分，主要是因為復雜的淺部結(jié)構(gòu)會導致面波傳播遠偏離大圓路徑，使得簡單按照大圓路徑傳播假設(shè)旋轉(zhuǎn)得到的切向分量互相關(guān)函數(shù)（Lin et al.,2009）既包含勒夫波也包含瑞利波.一個可能的矯正方法是利用垂直分量互相關(guān)函數(shù)的瑞利波和徑向分量互相關(guān)函數(shù)的瑞利波完全相同這一性質(zhì)，先通過旋轉(zhuǎn)水平分量的互相關(guān)函數(shù)得到真正反映瑞利波真實傳播路徑的徑向互相關(guān)函數(shù)（與垂向互相關(guān)函數(shù)相同；胡亞洲，2021），然后再旋轉(zhuǎn)獲得切向分量互相關(guān)函數(shù)以及更為可靠的勒夫波，從而為后續(xù)地殼淺部的三維徑向各向異性結(jié)構(gòu)成像提供更為可靠的數(shù)據(jù)基礎(chǔ).

4 討論與展望

4.1 面波直接成像方法的頻散數(shù)據(jù)和模型參數(shù)設(shè)置問題

在背景噪聲面波頻散提取時，通常要求臺站間距離大于2～3 倍波長，以避免噪聲互相關(guān)函數(shù)零時刻附近的信號（例如遠場體波）對面波信號的干擾.例如，對間距為百千米的臺站對，通常我們可以提取到～15 s 周期的面波頻散，可以反演中地殼深度的結(jié)構(gòu).但百千米及以上間距的臺站對難以提取到很短周期（例如1～2 s 或更高頻帶）的面波頻散，主要是由于短周期或高頻面波的衰減和散射較強所導致的.隨著臺站間距的逐漸減小，我們能提取到的面波頻散的最長周期越來越小，高頻頻散數(shù)據(jù)也越來越豐富，即越來越反映地殼淺層速度結(jié)構(gòu).

通常而言，中長周期（例如10 s 周期及以上）的面波頻散提取相對比較簡單，因為擬提取的頻散曲線通常與參考頻散曲線具有較好的相似性，一般不會出現(xiàn)短周期相速度頻散提取時出現(xiàn)的周期跳躍（cycle skipping）問題.需要注意的是，隨著周期變長（例如超過30 s 周期），群速度提取的誤差會遠高于相速度提取的誤差（Bensen et al.,2008）.因此，對于尺度在數(shù)百千米甚至數(shù)千千米的區(qū)域臺陣噪聲面波成像，我們通常僅使用相速度頻散數(shù)據(jù)，而不使用群速度頻散數(shù)據(jù).

對于密集短周期臺陣，通常臺陣布設(shè)的尺度在數(shù)千米到數(shù)十千米，由于儀器頻帶寬度的限制，大部分短周期地震儀很難提取超過10 s 的面波頻散，很多間距僅為千米的臺站對僅能提取頻率為幾赫茲的高頻面波頻散.由于雙臺高頻面波相速度頻散曲線提取存在較為嚴重的周期跳躍問題，在結(jié)構(gòu)非常復雜的地殼淺部區(qū)域，面波相速度可能從數(shù)百米每秒變?yōu)閿?shù)千米每秒，因此輔助頻散提取的參考相速度頻散曲線往往與實際觀測頻散存在很大的偏差，所以雙臺短周期或高頻面波相速度頻散曲線的拾取經(jīng)常會存在很大的誤差.但群速度頻散時頻圖像在短周期相對比較簡單，通過時頻分析方法可以有效自動提取群速度頻散曲線，其誤差通常小于相速度誤差.因此，對于小尺度臺陣，當區(qū)域結(jié)構(gòu)過于復雜，相速度頻散曲線難以提取的時候，我們通常先提取群速度頻散曲線，然后基于群速度頻散數(shù)據(jù)先反演獲得區(qū)域的三維速度結(jié)構(gòu)模型，然后再基于三維模型正演計算臺站對之間的相速度頻散，以之作為參考頻散再進一步提取更為可靠的相速度頻散，以避免短周期相速度頻散提取的周期跳躍問題.在一些短周期背景噪聲成像研究中，有的僅使用群速度頻散（例如,Zhou et al.,2018），有的僅使用相速度頻散（例如,Fang et al.,2015;Gu et al.,2019;Luo and Yao,2021），有的則是相速度頻散和群速度頻散聯(lián)合一起反演三維速度模型（例如,Chen et al.,2016a;Li et al.,2016）.

為了避免相速度或群速度頻散在拾取過程中可能存在的較大誤差，例如由于周期跳躍、信噪比不高、其他信號干擾所造成的問題，對于密集布設(shè)的臺陣，我們通常采用路徑束分析（path cluster analysis）的方法（Zhang Y Y et al.,2018）.該方法要求相似路徑的臺站對具有相似的頻散曲線，通過統(tǒng)計分析的方法來剔除路徑相似但頻散出現(xiàn)異常的頻散點甚至整條頻散曲線，從而進一步保障面波直接成像方法中頻散數(shù)據(jù)的可靠性.

基于面波頻散的物理屬性、頻散資料提取的頻帶范圍以及很多不同尺度臺陣的實際噪聲成像分辨率結(jié)果，我們通常可以有效反演的橫波速度結(jié)構(gòu)的最深深度一般小于1/4～1/6 的最長臺站間距.在采用面波直接成像方法進行三維橫波速度結(jié)構(gòu)反演時，我們設(shè)置的反演最深深度通常是最長周期面波波長，但可靠的反演結(jié)果通常在1/2 波長深度以內(nèi).由于周期越短面波的深度敏感核越窄，即淺部深度分辨率越高，所以在短周期面波資料豐富時，我們通常在淺部設(shè)置更密的深度方向網(wǎng)格點；隨著深度的逐漸增加，垂直方向設(shè)置的網(wǎng)格間距也逐漸增大，這樣可以更好地提升淺部結(jié)構(gòu)的分辨率，同時又不至于使得深部分辨率不足時而導致反演系統(tǒng)不穩(wěn)定的問題.上述這些基本的認識和實際經(jīng)驗有助于我們更好地進行面波直接成像的模型參數(shù)設(shè)置，使得反演系統(tǒng)的穩(wěn)定性和分辨率達到較好的平衡.

4.2 基于面波有限頻走時和波形反演的三維速度結(jié)構(gòu)成像

目前提出的面波三維直接成像方法基于面波射線理論，所提取的面波頻散走時僅反映從源到臺站之間的射線路徑上的結(jié)構(gòu)信息.很多研究表明，面波具有顯著的有限頻效應(yīng)（Nishida,2011;Spetzler et al.,2002;Yoshizawa and Kennett,2005），面波走時不僅對面波射線路徑上的結(jié)構(gòu)敏感，還對射線路徑之外較寬的區(qū)域的結(jié)構(gòu)敏感.一些面波兩步法成像的第一步即是開展了基于有限頻理論的二維相速度分布反演（例如,Ritzwoller et al.,2002;Yoshizawa and Kennett,2004;Zhao et al.,2020），但二維敏感核的計算是基于一維均勻?qū)訝钅Ｐ偷?，因此不同周期的面波都是沿著大圓路徑傳播，這對于反演全球尺度或較大區(qū)域尺度的地殼上地幔結(jié)構(gòu)總體而言是合適的.但對于結(jié)構(gòu)復雜的地殼淺部，速度結(jié)構(gòu)不僅存在強烈的橫向變化，也存在很大的垂向變化，這會導致短周期或高頻面波路徑會強烈的彎曲.不同周期面波的敏感深度也不同，很可能會出現(xiàn)同一個源和臺站之間不同周期的面波傳播路徑也存在較大的差異.如果在考慮面波傳播走時的有限頻效應(yīng)時不考慮面波傳播的彎曲效應(yīng)，很可能導致走時異常計算存在一定的誤差.因此，一些研究在計算面波走時的二維敏感核時也同時考慮面波傳播的非大圓路徑效應(yīng)（Lin and Ritzwoller,2010;Russell and Gaherty,2021）.

李成（2019）基于三維譜元法的波場模擬計算面波波場，并通過不同理論（射線理論，有限頻理論但假設(shè)面波沿大圓路徑傳播，有限頻理論但考慮面波射線路徑彎曲）計算得到的臺站處面波的走時異常，并與波形互相關(guān)方法得到的走時異常進行對比，發(fā)現(xiàn)同時考慮有限頻理論和面波射線路徑彎曲時具有最好的走時異常計算結(jié)果.基于這一分析，李成（2019）進一步提出了考慮面波有限頻效應(yīng)和射線路徑彎曲效應(yīng)的面波走時三維敏感核函數(shù)（圖9）計算方法以及面波三維直接成像方法，并初步應(yīng)用于臺北盆地地區(qū).相比于不考慮有限頻效應(yīng)的面波直接成像方法（Fang et al.,2015），新方法反演得到的地殼淺部的速度異常更強.該方法可以進一步提升面波走時三維直接成像方法在復雜地殼淺層結(jié)構(gòu)成像的精度，相比于三維波形成像具有更高的計算效率，同時也可以結(jié)合體波走時的三維敏感核，開展面波與體波走時聯(lián)合的三維成像.

圖9 臺北盆地兩個臺站間0.8 s 周期瑞利面波走時的三維敏感核函數(shù).（a）顯示了0.4 km 深度處的三維敏感核的水平剖面，黑色實線表示射線路徑.（b）顯示了三維敏感核的一條縱剖面，位置如（a）中紅色實線所示（修改自李成，2019）Fig.9 The 3D sensitivity kernel of the 0.8 s period Rayleigh wave travel-time between two stations in the Taipei Basin.(a) Shows the horizontal cross-section of the 3D sensitivity kernel at 0.4 km depth with the black curve representing the ray path.(b) Shows the vertical cross-section of the 3D sensitivity kernel with its location given as the red line in (a) (modified from Li,2019)

基于面波走時的三維直接成像方法（Fang et al.,2015;Hu et al.,2020;李成,2019;Liu et al.,2019）還并不能考慮更為復雜的波場效應(yīng)，例如多路徑效應(yīng).為了更好地利用波形信息，近來一些研究還計算了背景噪聲面波的三維敏感核（例如,Tromp et al.,2010），并直接開展基于背景噪聲面波波形數(shù)據(jù)的三維伴隨方法成像（adjoint tomography）或全波形成像，雖然計算量遠高于基于走時的面波成像方法，但可以獲得區(qū)域更高精度的各向同性三維橫波速度結(jié)構(gòu)模型（例如,Chen et al.,2014,川滇地區(qū);Gao and Shen,2014,美國Cascadia 地區(qū)）.最近，通過多分量背景噪聲面波波形的三維伴隨方法波形成像，Wang 等（2020）反演獲得了美國南加州地區(qū)三維徑向各向異性的橫波速度結(jié)構(gòu)模型.目前很多區(qū)域開展了線性密集臺陣的觀測，為了顯著降低三維全波形成像的計算量，Groos 等（2017）給出了基于線性臺陣高頻面波數(shù)據(jù)的二維全波形成像的完整流程.Zhang C 等（2018）在此基礎(chǔ)上進一步提出了線性臺陣背景噪聲面波伴隨方法的波形成像，并成功應(yīng)用于華北克拉通地區(qū)的寬頻帶線性臺陣數(shù)據(jù)，獲得了臺陣下方更為可靠的地殼和上地幔頂部二維橫波速度結(jié)構(gòu)模型.

目前國際上基于背景噪聲面波數(shù)據(jù)開展全波形成像的工作也越來越多，但針對短周期或高頻背景噪聲面波數(shù)據(jù)開展的全波形成像工作還非常少，可能與短周期背景噪聲面波波形的復雜性，以及反演需要較好的三維初始模型有一定的關(guān)系.此外，背景噪聲面波全波形成像中如何考慮噪聲源分布的影響也是目前國際上重要的研究方向.考慮到面波全波形成像通常需要較好的三維或二維初始模型，以避免反演陷入局部極小值，以及加速反演的收斂，因此可以考慮先采用基于面波走時的成像方法獲得較好的初始模型，然后再開展全波形成像.

4.3 基于頻散走時數(shù)據(jù)和其他數(shù)據(jù)聯(lián)合的三維直接成像

基于面波頻散走時數(shù)據(jù)的面波直接成像方法可以有效地獲得三維各向同性和各向異性的橫波速度結(jié)構(gòu)模型，由于面波頻散對縱波速度、密度等參數(shù)敏感性較弱，所以目前面波直接成像方法還難以反演除橫波速度之外的其他結(jié)構(gòu)參數(shù)信息.由于面波直接成像方法構(gòu)建了從面波走時到三維橫波速度結(jié)構(gòu)的正演方程，非常類似于三維體波走時成像的正演方程，因此Fang 等（2016）在三維面波直接成像（Fang et al.,2015）的基礎(chǔ)上，進一步提出了體波走時和面波頻散走時直接反演三維縱波速度和橫波速度結(jié)構(gòu)的聯(lián)合成像算法，并將該方法成功應(yīng)用于美國南加州地區(qū)以及美國大陸（Golos et al.,2018）.為了獲得更為穩(wěn)定的反演結(jié)果，該算法除了施加縱波速度及橫波速度模型參數(shù)的空間平滑之外，還增加了縱橫波速比的物理約束.為了考慮地形起伏對不同網(wǎng)格點下方的面波頻散計算的影響，Liu 等（2021）還進一步改進了模型網(wǎng)格參數(shù)化的方法，可以考慮地形起伏對面波頻散計算的影響.新方法成功應(yīng)用于我國川滇地區(qū)，構(gòu)建了川滇地區(qū)首個三維公共速度模型（Liu et al.,2021）.

體波走時和面波頻散走時直接反演三維縱波速度和橫波速度結(jié)構(gòu)后，可以通過相除的方法進一步計算得到縱橫波速比（VP/VS），但在數(shù)據(jù)覆蓋不是很好的區(qū)域得到的波速比可能存在較大的誤差.為了獲得更為可靠的三維波速比結(jié)構(gòu)，F(xiàn)ang 等（2019）進一步提出了體波走時和面波頻散走時聯(lián)合反演三維縱波速度以及三維波速比（VP/VS）結(jié)構(gòu)的新方法.新方法可以直接對三維波速比施加模型平滑，因此可以獲得更為穩(wěn)定可靠的三維波速比結(jié)構(gòu)，類似于面波頻散走時直接反演三維徑向各向異性參數(shù)的方法（Hu et al.,2020）.

不同的地震學或其他地球物理學數(shù)據(jù)對不同地下結(jié)構(gòu)參數(shù)具有不同的敏感性，如何在現(xiàn)有的三維直接成像的框架下聯(lián)合更多的數(shù)據(jù)開展精度更高、可靠性更好的三維結(jié)構(gòu)成像是未來需要進一步發(fā)展的方向.例如，聯(lián)合面波頻散走時數(shù)據(jù)和重力異常數(shù)據(jù)可以更好地直接反演三維橫波速度結(jié)構(gòu)（Du et al.,2021）.基于交替迭代的反演策略，吳萍萍等（2020）實現(xiàn)了電阻率和背景噪聲面波的三維聯(lián)合反演方法，聯(lián)合反演方法可以獲得同時滿足兩種不同類型數(shù)據(jù)的更為相似的橫波速度和電阻率結(jié)構(gòu)模型.

如何在目前面波頻散走時和體波走時（僅包括初至體波走時）聯(lián)合反演的框架下（Fang et al.,2016,2019），進一步聯(lián)合多震相體波走時數(shù)據(jù)（包括直達波、反射波等）、瑞利面波Z/H比數(shù)據(jù)、體波接收函數(shù)等數(shù)據(jù)，建立更為完善的聯(lián)合地震學反演方法，更好地反演獲得地下三維速度、波速比和界面結(jié)構(gòu)模型是今后重要的研究方向，這也為構(gòu)建更為可靠的區(qū)域公共速度模型提供重要的思路.