區域尺度地震體波和面波走時聯合成像：進展與展望

方洪健，劉 影，姚華建，張海江

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

2 南方海洋科學與工程廣東省實驗室（珠海），珠海 519080

3 中國科學技術大學地球和空間科學學院 地震與地球內部物理實驗室，合肥 230026

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

5 中國科學院比較行星學卓越創新中心，合肥 230026

0 引言

通過觀測到的天然地震的地震波數據重建地下的三維波速結構，又稱為地震層析成像，是地震學研究的最核心的問題之一.基于獲得的波速模型可進一步研究區域/全球構造演化、地震災害評估、斷層形態和地震發震構造等.地震層析成像的基礎是地震學的正反演理論.目前，由于計算機計算能力的大幅提升，研究者已經可以相對準確地模擬地震波在復雜介質中的傳播過程，并且在全球尺度可模擬到較高的頻率（如,Bishop et al.,2021;Mohr et al.,2017;Sun et al.,2021;Zhang et al.,2012）.地震學的正問題是相對穩定的.換句話說，在給定地下介質模型以及地震震源模型的情況下，天然地震學研究常用的一些地震震相可以被準確地且唯一地模擬出來.所以，目前對正演問題的研究大多針對于如何提高模擬高頻地震波的效率、如何考慮更加切合實際的介質（如模擬多次散射及衰減）以及如何更好地處理邊界等.然而，通過觀測到的地震數據確定地下的介質模型（地震學反問題）卻存在諸多問題，其中最為突出的是反演結果的非唯一性.著名地球物理學家John Claerbout 曾經說過：“Inverse theory is the fine art of dividing by zero (inverting a singular matrix)”（Claerbout,2014）.換言之，地球物理學家如果想要獲得可靠的結果，那么就必須處理好除零，或者處理好反演系統零空間的問題.關于如何減小或消除區域尺度地震層析成像中的零空間，即如何改善區域尺度層析成像波速模型，多數研究主要針對以下幾個方面.

首先，可以通過在正演過程中引入更好的物理近似以減小反演的零空間.由于效率較高，目前絕大部分的區域及全球尺度的地震層析成像還是基于射線追蹤的走時層析成像.射線追蹤常用的方法，如偽彎曲法（Um and Thurber,1987）、基于圖像的方法（Moser,1991）等，在復雜模型中通常會由于方法的局限性獲得不準確的走時及射線信息.考慮到基于射線追蹤獲得走時及射線路徑的局限性，基于快速行進法（Fast Marching Method）直接求解程函方程獲得走時場的方法被提出并被應用到地震成像中（de Kool et al.,2006;Lan and Zhang,2011;Rawlinson and Sambridge,2004;White et al.,2020）.該方法求解過程相對于最初的基于有限差分的方法（Vidale,1990）更為穩定.射線路徑可通過追蹤走時場的梯度獲得.該種方法在復雜的介質中也可以獲得較為準確的走時場，并且可以擴展到多種震相以及有地形起伏的各向異性介質（如,de Kool et al.,2006;Lan and Zhang,2011）.盡管基于求解程函方程的走時成像方法較基于傳統射線追蹤的方法有一定的改進，其更為準確的射線路徑以及計算的走時可在一定程度上減小反演的零空間，但其仍然是基于地震波波動方程的高頻近似，即認為觀測到的走時信息只對射線經過區域的地震波波速敏感.這導致基于線性化方法求解反問題時如果模型網格劃分較密，那么獲得的結果將會由于反演系統的欠定性而存在非唯一性，即反演系統具有較大的零空間.考慮到基于射線的成像方法的局限性，Dahlen 等（2000）和Hung 等（2000）發展了有限頻層析成像方法.有限頻成像將觀測數據從單一走時擴展到不同頻率下的走時，并且考慮不同頻率的走時對地下結構的敏感度隨著頻率的變化而變化，即低頻走時的敏感度相對高頻的走時的敏感度會擴展到較大范圍.值得注意的是，基于有限頻層析成像方法對傳統射線經過區域無敏感度，其敏感核形狀與香蕉-甜甜圈（banana-donut）相似.基于有限頻層析成像的敏感核不再是一條無限細的射線，即在同樣網格化的情況下，其反演系統的靈敏度矩陣相比于傳統的走時成像方法更加稠密，具有相對較小的零空間，所以在相同數據覆蓋下，有限頻層析成像理論上會獲得更穩定的結果.此外，有限頻成像方法還考慮了地震波的波前愈合現象.Maceira 等（2015）通過兩種成像方法的對比發現有限頻層析成像對低速異常的成像效果優于傳統走時成像方法.基于有限頻成像方法，Montelli 等（2004）首次觀測到全球地幔柱的分布及形態，盡管當時該方法較傳統方法的優勢存在一定的爭議（van der Hilst and de Hoop,2005;de Hoop and van der Hilst,2005）.最初的有限頻成像方法在求解敏感度矩陣時采用了動態射線追蹤的方法，其本質上仍屬于波動方程的高頻近似.Tromp 等（2005）提出了基于波動方程的伴隨層析成像方法.該方法通過直接求解波動方程，所引入的近似相對較少，可以較為準確地計算觀測數據對地下結構的敏感度.由于需要強大的計算資源支撐，該種方法目前的應用主要在構建區域及全球尺度模型，并且其觀測數據主要為相對低頻的波形數據（如,Bozda? et al.,2016;Fichtner et al.,2009;Lei et al.,2020;Tape et al.,2009;Zhu et al.,2012）.

其次，可通過采用更加符合實際的先驗信息以及更合理的網格化方案.地震成像中最常用的先驗信息主要有兩種: 第一種為假設獲得的模型是平滑的；第二種為假設獲得的模型不能較初始模型有太大差異.在線性化反演中，這兩種信息通?？梢酝ㄟ^設計合理的矩陣并整合到靈敏度矩陣，從而減小原始靈敏度矩陣的零空間.除了這兩種先驗信息，其他的先驗信息，如假設模型梯度是稀疏的，或者模型在小波域是稀疏的等，也被應用到了地震層析成像中并取得了良好的效果（如,Chiao and Liang,2003;Fang and Zhang,2014;Loris et al.,2007;Simons et al.,2011）.關于模型的網格化，在基于梯度下降的伴隨成像中通常使用求解正演方程的網格.由于基于梯度下降的方法收斂性較慢，通常地震層析成像中可利用目標函數極小值處其導數為零這一特征，直接求解線性方程組以加速收斂.然而，在基于求解線性方程組獲得模型變化量的線性化反演中，如何網格化模型有著不可忽略的影響.如果網格較為粗糙，那么反演系統的零空間較小，反演比較穩定，但其代價是模型的分辨率較低.反之，如果網格較為精細，則需要較強的先驗信息的約束才能獲得較為穩定的結果.所以，在通過求解線性方程的反演中，模型所采用的反演網格通常與正則化項是耦合的.實際反演中反演網格的選取目前沒有統一的標準，通常的做法是嘗試多種網格化以及規則化方案以檢驗反演模型的可靠性.此外，反演過程中采用多套網格或者采用非規則網格也在區域及全球成像中取得了一定的改善（Fang et al.,2019;Sambridge and Gu?mundsson,1998;Tong et al.,2019;White et al.,2021;Zhang and Thurber,2005;Zhou,2003）.

線性化反演中通過引入先驗信息的方法可以在一定程度上減小反演系統的零空間，使得反演穩定.但其缺點是所獲得的模型對先驗信息的依賴性較大，而且獲得的模型通常是一個單一的模型，其誤差難以估測.鑒于此，Sambridge（1999）提出直接對模型空間采樣的非線性反演方法.該種方法不依賴于梯度信息，故不受初始模型選擇的影響.更重要的是，該種方法還可以獲得模型的不確定度，這對后續模型的使用至關重要.這種思想被應用到了面波成像（Bodin and Sambridge,2009）、面波接收函數數據聯合反演（Bodin et al.,2012）、面波數據直接反演（Zhang et al.,2018）等.由于這種基于對模型空間進行采樣的方法需要進行大量的正演計算，且當模型空間的維度過高（比如大于1 000）時無法保證在有限的時間對模型空間進行足夠的采樣，所以該種方法目前只在少數區域有所應用.所幸的是，近期有研究顯示通過引入梯度信息的蒙特卡羅以及機器學習的方法可以進一步減小計算量（Earp et al.,2020;Zhang and Curtis,2019,2021;Zhao et al.,2021）.考慮到模型不確定度在解釋過程中的重要性，隨著計算機計算能力的繼續提升，基于該方法的層析成像將會逐漸取代傳統的線性化反演.值得注意的是，該種方法沒有通過減小零空間來增加獲得模型可靠性，而是通過搜索可以擬合數據及符合先驗信息的多個模型，最后通過模型的統計學參數，如均值及方差來獲得穩定的模型及其不確定度.

最后，由于觀測到的地震記錄上不同震相通常對地下結構具有不同的靈敏度，所以在觀測系統無法改變的情況下，通過整合更多的地震觀測信息是減小或者消除反演系統零空間最直接的方法.如上所述，通過擬合觀測到的全部波形信息，即全波形反演，來構建地下結構模型是地震層析成像的終極目標.但是，通過近30 年的研究顯示，真正意義的全波形反演目前還很難達到.其主要原因有以下幾個方面：首先，地震波場模擬比較耗時；其次，波形信息（反演所采用的目標函數）通常與波速結構存在非常強的非線性關系，導致優化過程中如果初始模型選擇不合理容易陷入局部極值；波形反演過程中目標函數的選擇、地震震源時間函數處理等都會影響到獲得結果.此外，由于波形的振幅信息通常會受到很多因素的影響，如地震震源、衰減結構、地震波傳播時受結構的影響會出現聚焦以及散焦、場地效應以及地震儀的耦合等，以至于絕大多數波形成像還無法使用波形的絕對振幅信息.所以，目前基于波形的層析成像主要針對單一震相，如初至波、面波等.或者在同時擬合多種震相時只擬合頻率較低的波形.相比之下，聯合反演同樣遵循利用更多地震記錄上的信息來構建地下波速結構.其基本思想是將地震記錄上的信息拆分成不同震相，如P 波、S 波、面波等.然后對于不同的震相可使用效率較高的正演方法，如求解程函方程.之后將所有震相的數據整合到同一反演系統中，從而達到利用更多波形信息以減小反演系統零空間的目的.

利用不同數據對模型具有不同靈敏度的聯合反演的方法在研究大地震斷層面滑移分布規律中已有較多應用（如,Cirella et al.,2018;Delouis et al.,2002;Inbal et al.,2017;Koketsu,2016;Konca et al.,2007;Yokota et al.,2011;Yue et al.,2014）.此外，不同種類的數據，如地震數據和大地電磁數據、地震數據和重力數據等，也可以通過結構約束，如交叉梯度，或者將不同物性之間的一些經驗關系整合到同一反演系統中來共同約束地下結構（如,Bao et al.,2021;Bennington et al.,2015;Lin and Zhdanov,2019;Maceira and Ammon,2009;Popovici et al.,2017;Roecker et al.,2004;Syracuse et al.,2016;Zhdanov et al.,2012）.對于地震成像方面已有很多綜述文章，感興趣的讀者可以參考近期的綜述文章，如Liu 和Gu（2012）、Rawlinson 等（2010）、Thurber和Ritsema（2015）等.本文主要針對區域尺度不同震相聯合反演獲得地殼上地幔三維波速結構的方法，介紹其基本原理以及目前的應用情況，著重介紹其在構建川滇公共速度模型中的應用，最后探討一些可進一步改善的方案.

1 地震數據聯合反演

1.1 區域尺度體波走時層析成像

區域尺度觀測到的地震體波（包括直達P、S 波，Pn 波和Sn 波）到時數據可以用于構建區域三維波速模型.在使用單種數據，如P 波數據進行成像時，基于線性化的體波走時層析成像方法通?？梢院喕扇缦路匠蹋?/p>

當然，S 波也可以整合到以上方程來同時獲得三維的VP和VS波速結構.其線性化的方程如下：

1.2 區域尺度面波層析成像

基于面波頻散數據的區域尺度面波成像通常分為兩步.首先將測量的不同周期的面波頻散數據進行成像獲得二維的相速度或者群速度面，然后再通過空間點上反演獲得的不同周期相速度或者群速度進行一維反演獲得地下的VS結構.最后通過整合不同空間點的一維VS結構便可獲得區域的三維VS結構.總的來說，該種成像方法將面波反演分解成兩個計算量相對較小的子過程.由于計算量較小，可以對兩步法獲得的模型進行詳細的誤差分析來評估模型的不確定度.但是兩步法反演也有著一定的不足.比如，如果某些周期頻散數據較少，盡管該周期對地下不同深度結構有著不同程度的約束，但由于不能獲得穩定的相速度或者群速度面，所以該種數據通常會被舍棄.而面波頻散數據的數量在不同周期差異很大是很普遍的.此外，通常在進行兩步法成像時并不考慮面波的非大圓路徑傳播的影響.由于短周期面波主要對地殼結構比較敏感，而地殼結構相比于上地幔結構通常具有較大的橫向差異.所以，對于利用噪聲數據互相關獲得的短周期面波數據在地殼結構復雜區域基于面波大圓路徑傳播的假設會產生較大誤差（李想等,2015;Xia et al.,2018）.考慮到兩步法面波成像的不足，Fang 等（2015）發展了基于射線追蹤的面波頻散數據直接反演的方法.其主要思想是將不同周期面波的深度敏感核和水平方向的敏感核整合到一起以建立從三維VS結構到面波頻散數據的方程.直接反演的線性化方程如下：

1.3 體波走時與面波相速度圖聯合反演

區域尺度地震圖上最顯著的震相為直達P 波、直達S 波（圖1）.由于震級較小或者較深地震面波通常不發育，所以在區域尺度上很少有直接利用區域地震面波進行成像的.近20 年來，通過對噪聲數據互相關恢復臺站之間格林函數的方法已被廣泛應用到不同尺度地殼及上地幔的地震成像中（如,Lin et al.,2007;Shapiro et al.,2005;Yang et al.,2007;Yao et al.,2006）.受噪聲源性質的影響，恢復的格林函數中面波的成分占主導（見圖1）.在地震走時層析成像中，利用直達P 波和S 波進行成像時在射線分布稀疏以及射線交叉較少區域（通常在淺部）其縱向分辨率較低.利用面波數據進行成像時在淺部具有較高的分辨率，但其縱向分辨率會隨著頻率的減小而降低.所以，聯合地震體波及面波數據進行成像可以利用兩種數據對結構敏感度具有互補的特性，降低反演系統的零空間.

圖1 以美國南加州為例的區域尺度觀測的地震波形及噪聲互相關記錄.（a）南加州地區寬頻帶地震臺站分布（藍色三角形）.（b）2022 年2 月6 日一個震級為3.1 級地震[位置見圖（a）紅色五角星]在CI.JVM 臺站垂直分量[位置見圖（a）標注的藍色三角形]不同頻段的波形記錄.其中圖（b）上波形為1.0 Hz 的高通濾波，其次為帶通濾波，頻率范圍為0.05～0.5 Hz.圖（b）下圖顯示離地震位置較近的臺站CI.LMH 及CI.JVM 臺站間互相關記錄.互相關利用了2022年2 月到3 月一個月的連續記錄.由圖可知區域尺度上地震P 波和S 波較為明顯，但面波不發育.聯合反演中的面波數據可由臺站間通過對連續記錄進行互相關獲得.全波形成像受計算資源限制，通常只擬合低頻的地震波形，而區域尺度低頻波形常常具有較低的信噪比，故需對用來反演的波形進行嚴格的質量控制Fig.1 The seismic waveforms and ambient noise cross-correlation function of two stations in Southern California.(a) The blue triangles show the distribution of broadband seismic stations in Southern California.(b) The vertical components of a magnitude 3.1 earthquake (marked as red star in the left panel) at station CI.JVM (marked as blue filled triangle in the left panel) at different frequencies.The right top waveform is high pass filtered at 1 Hz and the right middle waveform is band pass filtered at 0.05～0.5 Hz.The right bottom shows the ambient noise cross-correlation function between station CI.LMH and CI.JVM using one-month continuous records (from February to March 2022).P-wave and S-wave could be easily recognized in the earthquake recording,but surface wave is difficult to identify.For regional-scale joint inversion,the surface wave data can be obtained by ambient noise cross-correlation using continuous records.For full waveform tomography,only low frequencies are used due to the large computational cost.As the signal to noise ratio of low frequency waveform at region scale is relatively low,strict quality control should be applied on the waveforms when doing full waveform tomography

當原始頻散數據無法獲得時，體波走時成像可通過整合公開的不同周期面波的相/群速度進行聯合反演，以共同約束地下結構（Zhang et al.,2014）.其求解方程如下：

1.4 體波走時與面波頻散數據聯合反演

考慮到傳統聯合體波走時數據和面波相速度面數據反演方法的局限性，如面波數據與體波數據類型不一致而導致很難選擇兩種數據的權重以及短周期面波復雜傳播路徑的問題，Fang 等（2016）發展了區域體波和通過噪聲數據互相關獲得的面波數據的聯合反演，該方法可綜合利用體波和面波數據對地下不同區域結構的敏感性，獲得一個可以同時擬合不同數據的統一的、自洽的模型.這種方法本質上與West 等（2004）提出的聯合反演方法相似，即基于面波直接反演的方法同時反演體波和面波數據.具體來說，可通過將面波的深度敏感核整合到求解相/群速度面時的二維敏感核，以建立三維波速結構與面波頻散數據之間的關系，這樣面波數據的反演思路就可等同于體波數據的反演，進而可將面波反演與體波反演整合到同一反演系統以達到聯合反演的目的.該種反演思想被廣泛用于區域尺度地殼及上地幔的成像（如,Feng and An,2010;Nunn et al.,2014;Obrebski et al.,2011），所用的體波及面波數據皆來自天然地震.值得指出的是，新的聯合反演方法利用基于射線追蹤的短周期面波直接反演的方法（Fang et al.,2015），所以該方法考慮了短周期面波在復雜結構中并非沿著大圓路徑傳播的問題.由于面波直接反演方程與體波反演方程的相似性，可直接將二者整合到同一反演系統.新的聯合反演需求解的線性方程為：

合成數據測試結果證明了兩種數據對地下結構具有互補的靈敏度并且獲得的模型優于單獨反演的結果（圖2）.其在美國南加州地區的應用所獲得的模型能夠很好地刻畫斷層帶結構及盆地結構等.而且，獲得的模型較南加州地區參考模型能更好地擬合地震觀測記錄，證明了聯合反演方法具有較大的應用潛力.該種方法目前已被用于構建川滇區域公共速度模型（Liu et al.,2021）、阿拉斯加地區模型（Nayak et al.,2020）、北美板塊三維波速模型（Golos et al.,2018）等.Liu 等（2021）在構建川滇區域參考模型時還考慮了地形起伏對獲得模型的影響，這為后續帶有強烈地形起伏的利用體波和短周期面波聯合成像更廣泛的應用奠定了基礎.

圖2 體波和面波數據聯合反演合成數據測試.（a）輸入模型；（b）由聯合反演獲得的模型；（c）體波數據單獨反演的模型；（d）面波數據單獨反演的模型.左列顯示VS，右列顯示VP.由圖可知聯合反演相比于體波數據單獨反演在淺部有較大改善，較面波單獨反演對VP 以及深部模型改善較大，并且異常的整體幅度與輸入模型更為接近，可更好地約束VP/VS 模型（修改自Fang et al.,2016）Fig.2 Synthetic test of joint inversion of body-and surface-wave data.(a) Input model;(b) Recovered model from joint inversion;(c) Separate inversion using body wave only;(d) Separate inversion using surface wave only are shown from the top to bottom.The left and right column show the VS and VP model,respectively.The resolution at shallow depths is improved in the joint inversion compared to separate inversioin using the body wave only.Besides,the joint inversion improves VP model significantly compared to separate inversion surface wave only.The better resolved anomaly amplitude in the joint inversion leads to more realiable VP/VS model (modified from Fang et al.,2016)

1.5 川滇地區的應用

本節主要以構建川滇地區公共模型為例討論體波面波聯合反演方法的應用.首先概括川滇地區公共速度模型構建的必要性以及聯合反演模型的特征.之后詳細介紹聯合反演中所用的體波和面波數據、模型構建、分辨率分析以及模型驗證.

印度板塊與歐亞板塊碰撞導致了青藏高原的隆升以及高原內部和川滇地區大地震的頻繁發生，高原演化模型至今存在諸多爭議.一個較為準確的高分辨率三維巖石圈波速模型是進一步認識該區域的構造演化及進行地震災害評估分析的基礎.因此，川滇地區公共殼幔模型的構建是目前一大研究熱點.現已有研究利用不同的數據資料（如體波到時、面波頻散、面波振幅比、遠震P 波接收函數等），通過不同的地震層析成像方法（如地震體波走時層析成像、面波層析成像、接收函數、聯合反演等），獲得了川滇地區地殼及上地幔的結構（波速模型、Moho 面起伏模型等）.目前，大部分模型是基于單一數據反演所構建的，因此通常只能擬合單一類型的數據，各個模型之間的一致性還有待驗證.此外，受數據和方法的限制，川滇地區已有三維速度模型的橫向分辨率最高可達50～100 km，垂向分辨率約為10～20 km.如前文所述，相對于單一數據反演，地震體波及面波聯合反演可以在擬合兩種數據的同時，基于兩種數據對于介質速度結構具有互補的敏感性以獲得更為可靠的、分辨率更高的三維VP及VS模型.

Liu 等（2021）從中國地震臺網中心收集了2008 年10 月至2016 年9 月期間川滇地區230 個固定臺站的體波到時數據，共計約20 萬地震事件、94 萬P 波到時、92 萬S 波到時.具體數據質量控制和篩選過程如下：（1）保留初至到時，去除明顯不符合時距曲線的數據；（2）僅保留被10 個以上臺站同時記錄到的地震事件；（3）在此基礎之上，考慮到川滇地區地震事件分布極不均勻，將整個區域劃分為多個體積相等的長方體，每個長方體僅隨機選擇一個地震事件進行保留.最終保留了約33 000個地震事件、39 萬P 波初至到時、37 萬S 波初至到時，并構建了約138 萬地震對雙差走時數據.對于面波數據，Liu 等（2021）收集了約8 100 條瑞利面波相速度頻散數據，頻帶范圍為5～50 s.該數據包括從124 個固定臺站提取的6 910 條瑞利面波相速度頻散（Yang et al.,2019），從35 個臨時臺陣提取的544 條瑞利面波相速度頻散（Yao et al.,2010），以及從16 個臨時臺站和云南地區49 個固定臺站提取的640 條瑞利面波相速度頻散（Qiao et al.,2018）.

由于目前聯合反演方法還是基于線性化反演思路，因此初始模型的選擇會在一定程度上影響反演的結果，過于偏離真實模型的初始模型會導致反演陷入局部極值，獲得誤差較大的反演結果.本研究以Yang 等（2019）通過聯合反演獲得的川滇地區三維VS模型為基礎，通過地殼VP與VS經驗公式Brocher（2005）轉換獲得三維VP模型，并應用高斯平滑，構建了一個包含先驗信息且較為平滑的三維VP及VS模型.在體波及面波聯合反演程序中，Fang 等（2016）將模型設置成了相對于平均海平面的水平層狀等間隔網格.本研究在反演過程中采用了多尺度網格反演策略，首先將水平方向的網格設置為0.5°×0.5°進行反演，然后以此結果作為初始模型，再進行水平網格為0.25°×0.25°的反演.考慮到川滇地區的實際地形，將深度方向的網格點范圍設置為-5～100 km，網格間隔隨深度增加而增加（50 km 以上為5 km 間隔，50～70 km 為10 km 間隔，70～100 km 深度為30 km 間隔）.由于川滇地區高程起伏較大，在進行體波走時和面波頻散正演及敏感度矩陣計算時都考慮了實際的高程信息.其中，對于體波走時，在正演和敏感度矩陣計算時均考慮臺站高程.對于面波頻散走時，在正演頻散和計算深度敏感核時，考慮該地理位置的真實高程進行計算.為避免聯合反演被單一數據所主導，在反演過程中，每次迭代前通過估算體波和面波數據的誤差，來確定兩種數據的權重.最終，經過11 次迭代，體波走時的均方根殘差從0.93 s 下降至0.43 s，面波頻散走時的均方根殘差從4.16 s 下降至2.06 s.

基于川滇地區豐富的地震體波及背景噪聲面波數據，利用體波及面波聯合反演方法構建了川滇地區公共波速模型1.0 版本.該版本包括了地殼及上地幔頂部（深度70 km 以上）的三維VP及VS模型（圖3、4）.在地殼淺部（5 km 左右），地殼結構比較復雜.四川盆地顯示相對較低的VS，沿著龍門山斷裂西北顯示條帶狀的高速異常.四川盆地內部的速度異常有著橫向不均勻性，川西、川東之間存在較強的速度結構差異.在20 km 左右，以龍門山斷層為界，VP及VS模型在青藏高原和四川盆地分別為大范圍低速和高速異常.低VP及VS異常位于龍門山—麗江—小金河斷裂帶以西以及小江斷裂帶附近.值得注意的是，這兩個低速異常并不相通，而是被位于安寧河斷裂帶、則木河斷裂帶下方的高速異常所隔開.該高速異常的位置與峨眉山大火成巖省內帶相對應.因此，Liu 等（2021）認為位于西側的大范圍中下地殼低速異?？赡芊从沉饲嗖馗咴镔|的東南向擠出，而位于東側的小范圍中下地殼低速異常很可能并非來源于青藏高原.張智奇等（2020）結合上地幔三維VS分布，認為小江斷裂帶附近的小范圍低速異?？赡苁怯捎陂L英質地殼在較高溫下發生塑性變形甚至部分熔融而形成的.在下地殼至上地幔頂部深度（40 km 以下），VP及VS的高低速異常分布與Moho 面深度相對應.從地震重定位結果來看，大部分地震發生于高速異常內部或高低速異常邊界，3 級以上地震基本位于大型斷裂帶上.

圖3 川滇公共速度模型（1.0 版本）在不同深度（相對于平均海平面）的VP 及VS 分布.其中，黑色實線為區域主要斷層分布，紫色虛線為峨眉山大火成巖省的內、中、外帶，紅色圓點為4.0 級以上地震的重定位結果（修改自Liu et al.,2021）Fig.3 The community velocity model (V1.0) of southwest China at different depths (the depth sea level is 0 km).Solid black lines represent major faults in this region.Purple dashed lines mark the inner,intermediate,and outer boundary of the Emeishan large igneous province.The relocated events with magnitude larger than 4.0 are shown in red circles (modified from Liu et al.,2021)

為了驗證模型的可靠性，本研究進行了棋盤測試，以評價在模型當前觀測系統及網格設置下的模型分辨能力.為更好地模擬真實射線路徑，本研究選擇了最終反演所獲得的三維VP及VS模型作為初始模型，在水平和垂直方向加以±5%的擾動.考慮到數據在不同深度的分辨率能力不同，設置了隨深度而變化的速度異常尺度：中上地殼的速度異常橫向和垂向尺度約為50 km 和10 km，中下地殼的速度異常橫向和垂向尺度約為75 km 和20 km.分辨率測試結果顯示，在路徑覆蓋較好的區域，測試結果顯示基本能夠恢復所設置的速度異常，模型的橫向分辨率為50～75 km，縱向分辨率為10～20 km，相比于該區域已有模型具有較高的分辨率.為定量分析模型的不確定性，本研究通過隨機抽取數據重復反演來進行自舉法分析.結果顯示，在絕大多數區域，VP及VS模型的標準差在0.05 km/s 以內，自舉法分析獲得的平均模型與反演所獲得的模型基本一致.盡管自舉法分析所重復的反演次數有限，并不能真正地統計分析模型的不確定性，但這仍在一定程度上說明了通過當前數據和方法反演獲得的模型是較為穩定可靠的.

為進一步驗證模型的可靠性，該研究將所獲得的模型與川滇地區已有的三個模型進行對比，分別是：Xin 等（2019）通過體波走時單獨反演獲得的中國大陸巖石圈參考模型1.0 版本，Shen 等（2016）通過面波頻散單獨反演獲得的中國大陸地殼及上地幔三維VS模型，Yang 等（2019）通過面波頻散、ZH 振幅比和接收函數聯合反演獲得的川滇地區地殼及上地幔頂部三維VS模型.另外，選取了三種不同類型的數據（面波單點頻散、ZH 振幅比以及氣槍P 波到時），利用上述幾種模型分別進行正演，對比不同模型對于不同數據的擬合情況.結果顯示，不同模型之間存在整體一致性，但仍存在一定的差異，尤其是在淺部較為明顯.除了使用接收函數數據的Yang 等（2019）模型，其他模型在Moho 面處都是較為平滑的.相比于Shen 等（2016）和Xin等（2019），Liu 等（2021）和Yang 等（2019）對于面波頻散數據的擬合較好.Liu 等（2021）和Xin 等（2019）對于ZH 振幅比的擬合基本處于同一水平，優于Shen 等（2016），但次于Yang 等（2019）.值得注意的是，這里用于正演的面波單點頻散數據和ZH 振幅比數據均來自于Yang 等（2019）.Liu 等（2021）和Xin 等（2019）對于氣槍P 波到時數據的擬合程度基本一致.這說明通過體波及面波聯合反演所獲得的模型（Liu et al.,2021）對未參與反演的體波或者面波數據的擬合并不比通過單一數據反演所獲得的模型差（Shen et al.,2016;Xin et al.,2019）.

此外，我們還對比了利用本文介紹的兩種不同體波面波聯合反演方法獲得的模型的差異.圖4g-4l展示了Gao 等（2017）通過聯合體波走時和面波相速度圖的反演所獲得的青藏高原東南緣（101°E～105.5°E,27°N～32.5°N）三維VP及VS模型，其使用的是Zhang 等（2014）所發展的體波走時與面波相速度圖聯合反演方法.其VP模型橫向分辨率約為50 km，VS模型橫向分辨率在中上地殼可達50 km.對比發現，兩個模型在中上地殼（30 km）以上整體具有較好一致性，如龍門山斷裂帶下方在地殼淺部的低速異常，麗江—小金河斷裂帶以西的中下地殼低速異常等.但由于所使用的數據與方法不同，兩個模型仍然具有一定的差異，尤其在下地殼深度.例如，在Gao 等（2017）的VP及VS模型中，Moho 面深度附近的速度起伏較大，且VP及VS模型差異較大，對反映下地殼流的低速異常成像效果也有待提高.

圖4 川滇公共速度模型（1.0 版本）在不同縱切面的VP 及VS（a-f），以及Gao 等（2017）在相應縱切面的VP 及VS（gl）.剖面位置如圖3b 所示，紅色十字為4.0 級以上地震的重定位結果.斷層位置以箭頭標識：CHF：程海斷裂；LZJF：綠汁江斷裂；XJF：小江斷裂帶；RRF：紅河斷裂帶；DLSF：大涼山斷裂；LMSF：龍門山斷裂帶；LXJF：麗江—小金河斷裂帶Fig.4 Vertical profiles of the community velocity model (V1.0) of southwest China (a-f) and the model from Gao et al.(2017) (g-l).The locations of profiles are shown in Fig.3b.The plus signs are relocated earthquakes with magnitude larger than 4.0.The faults are marked by arrows and labeled: CHF: Chenghai fault;LZJF: Lvzhijiang fault;XJF: Xiaojiang fault;RRF: Red River fault;DLSF: Daliangshan fault;LMSF: Longmenshan fault;LXJF: Lijiang-Xiaojinhe fault

綜上所述，相較于單一數據反演，通過體波走時和面波頻散走時聯合反演不僅可以同時獲得三維VP及VS模型，并且所獲得的三維VP及VS模型具有更高的分辨率和可靠性.此外，聯合反演所獲得的模型也能較好地同時擬合兩種類型的數據.通過更可靠的VP及VS模型可以獲得泊松比模型，可為定量刻畫中下地殼物質組分及熔融程度提供新的約束.

2 展 望

隨著地震記錄數據的指數級增長以及計算機計算能力的提升，地震層析成像將會是人類認識地球內部過程不可或缺的工具.而如何更充分地挖掘利用地震記錄獲得高精度的地下結構將會是層析成像方法研究的核心內容.聯合地震記錄上不同震相走時信息的反演具有原理的簡單以及容易實現的優點，且其反演相比于全波形反演受初始模型影響相對較小.此外，由于走時信息通常是從高頻的地震圖上獲得，而波形反演所利用的波形信息以及有限頻層析成像中的到時信息一般通過相對低頻的波形獲得，在區域尺度成像中目前還沒有關于二者成像精度直接的對比.在南加州區域，基于走時信息的反演在某些區域獲得的波速模型對觀測數據的擬合甚至優于基于波形反演獲得的模型（Fang et al.,2016）.特別是基于機器學習算法的發展來識別地震以及進行震相的拾?。∕ousavi et al.,2020;Zhu and Beroza,2019）可以獲得大量的走時信息，這將推進基于走時信息的聯合反演的廣泛應用.所以，全波形反演在區域尺度成像中仍然不能取代基于不同震相走時信息聯合反演.二者的結合還可能獲得更高分辨率地下結構，如通過聯合反演獲得的模型還可作為波形反演的初始模型來進一步擬合更高頻的信息，這樣便可以繼續改善所獲得的模型.類似的思想其實已經在利用遠震數據進行上地幔成像時有所體現，如Rawlinson 和Fishwick（2012）發現利用面波成像獲得的模型作為初始模型可以保留數據約束較小區域的大尺度特征.此外，目前區域尺度的利用t*數據的衰減成像以及利用橫波分裂的各向異性成像均依賴于獲得的各向同性模型，所以聯合反演獲得的模型還能為其他模型的建立提供基礎.今后的聯合反演方法的進一步發展及應用可能體現在以下幾個方面.

2.1 反演模型誤差估計

基于線性化反演的零空間可通過設置恰當的反演網格來減少或消除，如設置隨著射線覆蓋而變化的網格（Chiao and Liang,2003;White et al.,2021;Zhou,2003），但是由于觀測數據總是會有誤差，這種誤差將會傳遞到所獲得的模型.如何刻畫模型的誤差/不確定度對后續結果解釋至關重要.由于觀測數據的誤差通常很難估計，基于對模型空間采樣的非線性反演方法，如基于馬爾科夫鏈的蒙特卡羅方法可通過假設誤差分布遵循一定的物理模型，而同時獲得多個可以擬合數據的模型以及誤差分布的模型，進而刻畫獲得模型的不確定度（Galetti et al.,2017;Zhang and Curtis,2020）.另外，在地震活動性較高的區域，反演模型可通過對數據進行采樣，如某些地震在空間上相鄰，那么反演時便可通過隨機選取其中某一地震進行多次反演，進而估算由數據誤差引起的模型的不確定度（Fang et al.,2019）.

2.2 降低不同種類參數反演中的耦合

目前基于走時信息獲得各向異性成像的方法已被廣泛應用于區域各向異性結構成像（如,Eberhart-Phillips and Mark Henderson,2004;Huang et al.,2011;Wang and Zhao,2008），其公式推導主要基于弱各向異性模型（Backus,1965;Smith and Dahlen,1973）.即各向異性反演過程中假設波速隨方位的變化可以分解為各向同性的波速加上一些隨方位變化的改正量，然后同時反演各向同性波速以及控制改正量的參數，進而同時獲得各向同性波速以及各向異性結構.不難看出，各向同性參數與各向異性參數間有著強烈的耦合，即在各向異性強的地區如果只做各向同性結構成像，那么各向異性分量可能會投射到各向同性模型中，造成假象.反之，如果在各向異性較弱的區域反演考慮各向異性，由于反演參數過多，那么數據觀測誤差就有可能傳播到獲得的各向異性模型.而這些假象通常很難通過合成數據測試進行評估，所以目前還沒有比較好的方法消除其影響.

可幸的是，目前獲得的各向異性模型與區域地質在不同程度上有著相關性，而且通過不同方法，如面波方位各向異性、接收函數方位各向異性、橫波分裂等獲得的各向異性存在可比性.所以，盡管各向異性模型誤差難以評估，研究者仍然可以通過模型研究區域構造演化等問題.雖然線性化反演時各向同性和各向異性之間的耦合難以消除，但反演過程中通?？梢圆捎靡恍┘s束以減小模型參數之間的依賴.如通?？梢韵确囱莞飨蛲阅Ｐ?，待模型收斂后再繼續反演各向異性模型（Liu et al.,2019），測試結果表明這種分步式的方法要優于同時反演所有參數.實際上，這種思路已被用于一些高精度地震層析成像中，如基于雙差走時層析成像便是首先擬合絕對到時信息，然后逐漸加入雙差走時信息，以逐步改善震源區波速結構以及地震的位置，雖然其實現方式是通過對不同種類的數據予以不同權重，然后在反演過程中隨著迭代次數的改變而改變權重的大?。╖hang and Thurber,2003）.所以，通過聯合反演獲得的各向同性模型可以為之后各向異性反演提供一個更好的初始模型，以最大限度減小模型參數之間的耦合，從而獲得更加可靠的各向異性模型.當然，各向異性參數也可以作為模型參數加入到聯合反演系統，以同時獲得各向同性模型和各向異性模型，其耦合也可通過對不同數據在不同的迭代階段設置不同的權重而得到控制.此外，各向異性和各向同性參數的耦合程度還取決于射線路徑的分布.Huang 等（2015）通過詳細的理論測試發現，在射線路徑分布單一時，各向異性結構與各向同性速度異常會產生嚴重耦合.只有當射線路徑覆蓋足夠，且方位角分布均勻時，才能獲得較為可靠的各向異性模型（Huang et al.,2015;黃周傳,2022）.

2.3 聯合反演框架

區域尺度的層析成像中，目前應用較多的聯合反演有基于面波及接收函數的聯合反演（如,Bao et al.,2015;Chen and Niu,2016;Gama et al.,2021;Guo et al.,2015;Julià et al.,2000;Obrebski et al.,2015;Shen et al.,2013;Ward et al.,2014）、基于體波及面波走時的聯合反演（Eberhart-Phillips and Fry,2017;Fang et al.,2016;Guo et al.,2018;Obrebski et al.,2011,2012;Yang et al.,2021）、遠震走時及區域地震波形聯合反演（如,Schmid et al.,2008）、遠震波形數據及噪聲互相關的聯合反演（Gao,2016;Wang et al.,2021）以及面波、接收函數及振幅比的聯合反演（Berg et al.,2018;Shen and Ritzwoller,2016;Yang et al.,2019）.其中，面波數據振幅比數據對約束近地表波速結構尤為重要，通?？梢员徽系矫娌ㄏ嗨俣然蛘呓邮蘸瘮捣囱葜幸酝瑫r約束從淺到深的波速結構（如,Chong et al.,2016;Li et al.,2016;Lin et al.,2012）.此外，近期研究表明遠震記錄中的極化方向數據也能較好地約束地表的P 波及S 波結構（Park and Ishii,2018）.這種數據已被整合到接收函數、面波振幅比數據中以共同約束淺地表結構（Xiao et al.,2021）.實際上，這些從地震記錄上提取的信息皆可整合到同一反演系統中，從而獲得一個統一的模型.具體來說，除了區域地震及互相關獲得的短周期面波數據，遠震體波以及長周期面波也可以整合到反演系統以獲得上地幔結構.如Chang 等（2010）及Liu 和Zhao（2016）等整合近震、遠震S 波走時數據以及瑞利面波相速度數據，構建了從地表到700 km 左右的高精度S 波波速模型，并論證了其相對于單獨反演的優越性.接收函數以及振幅比數據可通過單點的信息加入到反演系統，類似體波以及面波二維相/群速度圖聯合反演方式（Eberhart-Phillips and Fry,2017;Zhang et al.,2014），從而構建一個多種數據共同約束的P 波和S 波波速模型（圖5）.

圖5 地震學數據聯合反演框架Fig.5 The framework of joint inversion with different kinds of seismic data

2.4 定量化解釋

層析成像獲得的波速模型本身可以服務于其他地震學研究，如通過三維的波速模型可以更好地約束地震的震源機制、大地震的破裂過程，以及通過強地面運動模擬為地震災害評估提供參考等.在通過模型研究具體科學問題時，目前模型的解釋多數只局限于定性的描述（主要是將波速異常對應于巖性、溫度、含水量等方面），如低速異?？赡艽砹黧w、沉積物、破碎帶，高速異常通常和基底、火成巖以及穩定地塊聯系起來，當然這在一定程度上受控于反演本身的不確定性以及多解性.隨著越來越多的觀測數據的加入、更符合實際的物理模擬以及更準確的不確定度的刻畫，層析成像獲得的模型的定量化解釋將會成為一種新的趨勢，比如俯沖帶低速異?？梢院秃柯撓灯饋恚–ai et al.,2018），火山地區低速異?？梢苑赐瞥鰩r漿的熔融程度以及巖漿囊的體積（Li et al.,2016），或者在貝葉斯的框架下直接由地球物理數據定量約束地下的巖性溫度等（Afonso et al.,2013），結合其他學科觀測以及動力學等模擬，層析成像波速模型將會在定量刻畫地球內部過程和地震災害評估中發揮重要作用.