陳 飛, 于海英, 高麗娜,3

(1．中國科學(xué)技術(shù)大學(xué) 地球與空間科學(xué)學(xué)院地震與地球內(nèi)部物理實(shí)驗(yàn)室，合肥 230026；2．上海市地震局地震監(jiān)測(cè)中心，上海 201203；3．吉林大學(xué) 地球探測(cè)科學(xué)與技術(shù)學(xué)院，長(zhǎng)春 130026)

0 引言

大約45 Ma前，印度-歐亞大陸的碰撞造成青藏高原的隆起和抬升[31-32]，印度-歐亞大陸之間的連續(xù)匯聚導(dǎo)致碰撞帶地殼發(fā)生了不小于1 500 km的縮短。如何調(diào)節(jié)或吸收如此大規(guī)模的匯聚量，一直是地球科學(xué)界爭(zhēng)論的熱點(diǎn)。而四川盆地，西邊與青藏高原東部邊界龍門山斷裂帶相接，東邊毗鄰揚(yáng)子克拉通西緣，因此四川地區(qū)既是檢驗(yàn)不同補(bǔ)償模式和理解陸內(nèi)形變機(jī)制最為理想的研究“窗口”，同時(shí)也是探索高原側(cè)向增生、研究盆山相互作用，探討2008年汶川大地震和2013年蘆山地震發(fā)生的深部構(gòu)造背景的關(guān)鍵地區(qū)。

前人在青藏高原的東南邊緣已經(jīng)進(jìn)行了廣泛的地球物理成像研究，包括大地電磁成像[15,49]，地震體波層析成像[3-4,7-9,27-28,39]，面波層析成像[6,44-46]，接收函數(shù)分析[43]，以及接收函數(shù)和面波的聯(lián)合反演[35]。然而，在大多數(shù)情況下，不同研究得到的模型存在一定的不一致性，特別是在龍門山斷裂帶區(qū)域[1-2,4,9,27-28,40-42,49-50]，因此對(duì)四川地區(qū)的結(jié)構(gòu)進(jìn)行成像有進(jìn)一步研究的必要。

由于地球物理反演固有的多解性，為了減少解的非唯一性，一條重要的途徑就是增加約束條件，即利用多種地球物理數(shù)據(jù)進(jìn)行聯(lián)合反演。重力-地震數(shù)據(jù)聯(lián)合反演成像思想早在1980年就被Lines提出[20]，隨后Lees等[21]提出了P波走時(shí)和重力聯(lián)合成像方法，Obrebski等[26]進(jìn)行了遠(yuǎn)震體波和面波數(shù)據(jù)的聯(lián)合反演，獲得了中國大陸巖石圈的剪切波速度模型。Zhang等[48]開發(fā)了一種聯(lián)合反演方案，結(jié)合局部或區(qū)域地震體波到時(shí)和面波頻散數(shù)據(jù)，以更好地提高地下的剪切波速度結(jié)構(gòu)。通過進(jìn)一步考慮面波頻散數(shù)據(jù)對(duì)Vp和Vs的靈敏度，F(xiàn)ang等[16]提出了一種新的聯(lián)合反演算法，可以使用體波到時(shí)和面波頻散數(shù)據(jù)來同時(shí)改善Vp和Vs模型。通常面波層析成像不能很好地分辨地下的尖銳界面，相比之下，接收函數(shù)對(duì)速度在深度上的尖銳變化很敏感。因此，為了更好地估計(jì)地下介質(zhì)在垂向上的速度變化，接收函數(shù)和面波頻散的聯(lián)合反演已經(jīng)逐漸得到較廣泛地應(yīng)用,在中國區(qū)域，已有不少研究組做過區(qū)域尺度的面波和接收函數(shù)的聯(lián)合反演研究[11,43]。

噪聲面波成像目前已經(jīng)得到廣泛地應(yīng)用，但是成像的橫向分辨率取決于地震臺(tái)站的空間分布。相比較而言，由于重力異常隨著距離的增加而振幅顯著減小，所以重力數(shù)據(jù)對(duì)觀測(cè)點(diǎn)附近的介質(zhì)有很好地約束，也即重力通常在淺部具有較高的橫向分辨率。但是重力反演與任何其他勢(shì)場(chǎng)方法一樣，重力異常的解釋有著顯著的不確定性。事實(shí)上，淺部弱的密度異常和深部強(qiáng)的密度異?？梢栽诘乇懋a(chǎn)生相同的重力信號(hào)。為了解決上面這些問題，得到可靠的地下速度結(jié)構(gòu)，基于速度和密度之間的經(jīng)驗(yàn)關(guān)系,Maceira等[24]提出了將面波頻散和重力異常觀測(cè)值結(jié)合到一個(gè)系統(tǒng)進(jìn)行聯(lián)合反演的算法，以獲得剪切波速度-密度自洽的高分辨三維橫波速度模型。他們使用GRACE衛(wèi)星[38]的重力數(shù)據(jù)以及高分辨率的面波數(shù)據(jù)[23]，進(jìn)行聯(lián)合反演得到了西部塔里木盆地的3D橫波速度模型，該模型可以同時(shí)比較好地?cái)M合該地區(qū)的面波頻散和重力異常兩種不同的數(shù)據(jù)。

為了更好地得到四川盆地的巖石圈Vs結(jié)構(gòu)，在四川盆地區(qū)域進(jìn)行了面波頻散和衛(wèi)星重力數(shù)據(jù)的聯(lián)合反演。在反演過程中，Vs模型受到了由He等[19]提供的關(guān)于中國大陸地區(qū)高分辨率莫霍面起伏數(shù)據(jù)的約束。通過聯(lián)合反演得到了一個(gè)可靠的高分辨率速度-密度自洽的三維剪切波速度模型。該模型與已發(fā)布的四川地區(qū)的速度模型在大的特征上基本一致，但該模型進(jìn)一步提高了主要特征的分辨率，并且能夠同時(shí)比較好地?cái)M合面波和重力數(shù)據(jù)。

1 數(shù)據(jù)和方法

1.1 數(shù)據(jù)

本研究采用的面波數(shù)據(jù)是8 s～70 s周期范圍內(nèi)的瑞利波相速度，該組面波頻散數(shù)據(jù)來自Shen等[34]的最新面波層析成像研究。他們使用收集了中國多個(gè)密集臺(tái)陣(CEArray, China Array, NECESS, PASSCAL, GSN)及其周邊地區(qū)(Korean Seismic Network, F-Net, KNET)總共2 000多個(gè)地震臺(tái)站的波形數(shù)據(jù)，進(jìn)行了背景噪聲面波層析成像，并且整合了華南、東北部分地區(qū)已有的地震面波層析成像的結(jié)果[52-53]，得到了整個(gè)區(qū)域內(nèi)8 s～70 s的瑞利波相速度和群速度分布。其模型中的斷層特征與已知的地質(zhì)和構(gòu)造特征具有非常好的相關(guān)性, 如圖1所示，短周期(10 s～20 s)的相速度分布反映了上中地殼S波速度的平均變化情況，其低速異常區(qū)與沉積盆地的分布密切相關(guān)，高速區(qū)與隆起構(gòu)造相關(guān)；中等周期(30 s～40 s)的相速度分布主要反映了研究區(qū)下地殼、殼幔過渡帶S波速度的變化情況；長(zhǎng)周期(50 s～70 s)的相速度分布圖主要反映了上地幔頂部S波速度的橫向變化。

研究采用的布格重力異常數(shù)據(jù)是從Bureau Gravimétrique International (bgi.omp.obs-mip.fr)提供的全球模型EMG2008提取得到的(圖2(a))。具有更好的準(zhǔn)確性和分辨率。因?yàn)楸狙芯康某上衲繕?biāo)區(qū)域是巖石圈的Vs結(jié)構(gòu)，所以為了去除地幔深部大尺度結(jié)構(gòu)和地幔對(duì)流對(duì)重力數(shù)據(jù)的影響，我們采用了中值濾波的方法來過濾長(zhǎng)波長(zhǎng)重力異常，即對(duì)每個(gè)節(jié)點(diǎn)去除節(jié)點(diǎn)4° 范圍內(nèi)重力數(shù)據(jù)的平均值(圖2(b))。

1.2 聯(lián)合反演方程

為了能夠?qū)崿F(xiàn)聯(lián)合反演，面波反演方程組與重力反演方程組需要統(tǒng)一到一個(gè)系統(tǒng)中。面波頻散數(shù)據(jù)主要與橫波速度有關(guān)，重力數(shù)據(jù)經(jīng)過速度與密度的轉(zhuǎn)換，也與橫波速度有關(guān)。用m0表示初始S波速度模型，聯(lián)合反演系統(tǒng)方程表示為式(1)[24]。

(1)

圖1 10 s、20 s、30 s、40 s、50 s、70 s的瑞利面波相速度分布圖Fig.1 Rayleigh wave phase velocity map for 10 s, 20 s, 30 s, 40 s,50 s, 70 s, periods(a)10 s; (b)20 s; (c)30 s; (d)40 s; (e)50 s; (f)70 s

圖2 四川地區(qū)重力異常分布Fig.2 Comparison of observed gravity anomaly data and predicted gravity data(a)原始布格重力異常；(b)平滑濾波后的重力異常分布；(c)基于面波單獨(dú)反演的模型預(yù)測(cè)的重力異常分布；(d)基于聯(lián)合反演的Vs模型預(yù)測(cè)的重力異常分布

其中：Ds、Dg分別表示面波頻散數(shù)據(jù)對(duì)橫波速度Vs的偏導(dǎo)矩陣和重力異常對(duì)橫波速度Vs的偏導(dǎo)矩陣；ωs、ωg分別表示平衡面波數(shù)據(jù)和重力數(shù)據(jù)對(duì)模型約束的權(quán)重；rs、rg分別表示面波數(shù)據(jù)和重力數(shù)據(jù)的殘差；p(T)表示不同面波周期的權(quán)重，本研究加大了長(zhǎng)周期面波頻散數(shù)據(jù)的權(quán)重；Δ、S分別表示深度方向上和水平方向上的平滑正則化矩陣；W表示先驗(yàn)?zāi)Ｐ蚼a的約束矩陣?？梢杂肔SQR算法[29]求解上面的聯(lián)合反演系統(tǒng)方程組。同時(shí)，面波和重力的聯(lián)合反演算法可用框架表示(圖3)。

圖3 面波和重力聯(lián)合反演算法框架Fig.3 The algorithm framework of joint inversion

用DISPER80[36]求解聯(lián)合反演的面波頻散正演部分的計(jì)算。理論上，面波頻散是地下介質(zhì)橫波速度、縱波速度和密度的函數(shù)[10]。但是，面波頻散對(duì)縱波速度和密度的敏感性要遠(yuǎn)小于對(duì)橫波速度的敏感性[10,37],因此在反演中通常只考慮橫波速度對(duì)頻散曲線的影響。

1.3 重力異常的理論計(jì)算和近似

由牛頓萬有引力定律可知，一個(gè)區(qū)域V內(nèi)的密度異常在空間任何一點(diǎn)(x,y,z)處引起的重力異?？?/p>

圖4 單獨(dú)面波反演和聯(lián)合反演每個(gè)格點(diǎn)面波頻散數(shù)據(jù)殘差分布比較Fig.4 Comparison of surface wave dispersion data RMS residual from separate inversion and joint inversion(a) 單獨(dú)反演；(b)聯(lián)合反演

以表示為：

(2)

其中：G是牛頓引力常數(shù);r是空間點(diǎn)與密度異常區(qū)域V中心的空間幾何距離，一般密度異常是空間的復(fù)雜連續(xù)函數(shù)，式(2)體積分很難解析求解，通常采用數(shù)值求解。將研究區(qū)域離散成了許多個(gè)棱柱。對(duì)于其中任意大小的棱柱，其在空間一點(diǎn)產(chǎn)生的重力異?？梢杂肞louff公式表示[30]。

yiln(Rijk+xi)]

(3)

在觀測(cè)點(diǎn)的重力異?？梢园阉欣庵呢暙I(xiàn)加起來(圖4)：

(4)

將研究區(qū)域進(jìn)行精細(xì)的網(wǎng)格劃分，模型的網(wǎng)格數(shù)目巨大，這樣導(dǎo)致在計(jì)算地面上某一點(diǎn)重力異常時(shí)，其計(jì)算量非常大，計(jì)算速度較慢。為了提高重力異常正演計(jì)算的效率，在程序中采取了近似方法。因?yàn)橹亓Ξ惓ｋS著距離成二次方衰減，所以在實(shí)際重力數(shù)據(jù)計(jì)算過程中只考慮其周圍幾圈網(wǎng)格的貢獻(xiàn)，這將達(dá)到很高的精度[24]。

1.4 地震波速度與密度的關(guān)系

地震波速度與密度之間的經(jīng)驗(yàn)關(guān)系。常用的一個(gè)經(jīng)驗(yàn)關(guān)系為式(5)[12]。

α=3.125ρ-2.40

(5)

其中：α是P波速度;ρ是介質(zhì)密度。Birch關(guān)系主要適合于地殼和上地幔中密度較大、埋藏較深的巖石。另外一個(gè)非線性關(guān)系式Nafe-Drake[23]關(guān)系常用于淺層沉積物：

α=6.86-7.55ρ+2.64ρ2

(6)

Maceira and Ammon[24]利用一個(gè)分段函數(shù)式可以把這兩個(gè)經(jīng)驗(yàn)關(guān)系式結(jié)合起來:

ρ(α)=cos2φ ρ2(α)+sin2φρ1(α)

(7)

其中ρ2(α)、ρ2(α)分別是(5)式和(6)式，其中

(8)

由于我們希望得到的是橫波速度模型，這里采用Brocher[13]橫縱波速度經(jīng)驗(yàn)關(guān)系來進(jìn)行轉(zhuǎn)換：

(9)

2 四川地區(qū)的聯(lián)合成像應(yīng)用

將研究區(qū)域在水平方向上進(jìn)行了網(wǎng)格離散化。網(wǎng)格大小在經(jīng)度和緯度方向上均是0.5°，在深度方向上，對(duì)模型進(jìn)行了不均勻分層：0 km～50 km按每2.5 km劃分為一層, 50 km～100 km以下是每 5 km一層，100 km以下是每10 km一層，計(jì)算深度從地表到180 km深度為止。為了反演系統(tǒng)的穩(wěn)定性，對(duì)最終Vs模型在深度和水平方向上加了平滑約束，同時(shí)在深度方向上在Moho面附近放松了平滑約束，允許速度自由跳變。莫霍面起伏是從利用中國大陸寬頻帶永久地震臺(tái)的接收函數(shù)分析研究中獲得的[19]。采用前人在該區(qū)域的一維平均速度模型作為初始模型，先進(jìn)行多次單獨(dú)面波反演選取反演參數(shù)得到三維速度模型，再用該模型作為聯(lián)合反演的初始速度模型進(jìn)行聯(lián)合反演，得到最終的模型。對(duì)于單獨(dú)的面波反演，面波數(shù)據(jù)的RMS殘差從0.176 km/s下降到0.028 km/s(圖4(a)),但是從單獨(dú)面波反演得到的Vs模型擬合重力數(shù)據(jù)很差，重力數(shù)據(jù)的RMS殘差是44.4 mgal(圖3(c))。對(duì)于加上重力數(shù)據(jù)進(jìn)行聯(lián)合反演得到的模型，面波數(shù)據(jù)的最終RMS殘差變?yōu)?.039 km/s，比單獨(dú)面波的殘差稍大一些，但仍在研究誤差允許的范圍內(nèi)(<0.05 km/s)(圖4(b))，而聯(lián)合反演得到的模型擬合重力數(shù)據(jù)的殘差只有6.2 mgal，可見聯(lián)合反演能夠同時(shí)比較好地?cái)M合面波數(shù)據(jù)和重力數(shù)據(jù)(圖3(d))。另外我們注意到，實(shí)際上重力數(shù)據(jù)的測(cè)量誤差通常比面波的誤差要少很多，單獨(dú)面波反演有可能過度地?cái)M合了面波數(shù)據(jù)中的噪音，而聯(lián)合反演模型犧牲了一部分對(duì)誤差較大的面波數(shù)據(jù)的擬合，更好地?cái)M合測(cè)量誤差更小的重力數(shù)據(jù)，從這個(gè)角度說明聯(lián)合反演的模型更可靠。

驗(yàn)證聯(lián)合反演的可靠性，需進(jìn)行模型的恢復(fù)度測(cè)試，輸入模型為聯(lián)合反演的模型，并且用相似度來定量描述模型的恢復(fù)度情況。每個(gè)網(wǎng)格點(diǎn)相似度的計(jì)算公式為：

(10)

其中：DVr是反演恢復(fù)出來的速度異常；DVt是真實(shí)的速度異常。

從公式(10)可知，如果模型百分之百被恢復(fù)出來，則S值為“1”，如果恢復(fù)出來的模型與真實(shí)模型完全相反，則S值為“0”，即S越接近“1”,模型的恢復(fù)度測(cè)試越好，越接近“0”，模型的恢復(fù)度越差。圖5為5 km、10 km、20 km、40 km深度處的模型恢復(fù)度測(cè)試情況，圖5左邊為輸入模型，圖5中間為面波單獨(dú)反演的恢復(fù)度測(cè)試得到的模型，圖5右邊為聯(lián)合反演的恢復(fù)度測(cè)試得到的模型結(jié)果，圖6是根據(jù)圖5的模型用公式(10)計(jì)算出來的模型相似度情況。從圖5和圖6可以看出，面波單獨(dú)反演和聯(lián)合反演的模型分辨率情況都還可以，細(xì)致比較兩者的分辨率情況可以發(fā)現(xiàn)，在淺部5 km、10 km、20 km處聯(lián)合反演的模型分辨率相比面波單獨(dú)反演在有些區(qū)域有著明顯提高(比如20 km的分辨率情況)，可以看出在四川盆地里面，聯(lián)合反演的恢復(fù)情況比單獨(dú)反演的恢復(fù)情況更好一些，這是由于重力數(shù)據(jù)為反演提供了額外的數(shù)據(jù)，而深部?jī)烧咧g分辨率差異不大，這是因?yàn)橹亓Φ挠绊戨S著距離衰減很快，所以在深部對(duì)模型的分辨率改善不大是合理的。另外由圖6可以知道，不管是面波單獨(dú)反演還是聯(lián)合反演，兩者在斷層附近的分辨率都比較差，這是由于兩者均采用了平滑反演策略，在反演中加入了平滑正則化項(xiàng)約束導(dǎo)致的。

3 結(jié)果和討論

對(duì)四川地區(qū)進(jìn)行了反演成像，得到了四川地區(qū)精細(xì)的三維S波速度模型，探討聯(lián)合反演過程中重力異常數(shù)據(jù)的引入對(duì)最終三維S波速度模型的影響，詳細(xì)對(duì)比了面波單獨(dú)反演和聯(lián)合反演兩者不同深度的速度結(jié)果，如圖7所示。為了詳細(xì)研究青藏高原東南邊緣的速度結(jié)構(gòu)沿深度的變化情況，對(duì)聯(lián)合反演模型沿經(jīng)度、緯度方向及斜向截取了4條剖面(圖8)，這些剖面跨過了多個(gè)不同的構(gòu)造單元區(qū)域，可以觀察到構(gòu)造單元在深部的速度異常變化，速度剖面圖中上面的黑線代表地形起伏，下方的黑線代表Moho面深度[19]。

對(duì)比面波單獨(dú)反演以及面波和重力聯(lián)合反演的結(jié)果，可以看出兩個(gè)模型在大的特征上都給出了相似特征(圖7)，這是合理的，說明模型的穩(wěn)定性。但是速度的絕對(duì)大小有一些差異，在一些細(xì)節(jié)特征上，聯(lián)合反演給出了更好的結(jié)果。如圖7中的17.5 km速度切片所示，聯(lián)合反演在四川盆地的內(nèi)部斷層兩邊給出比較明顯的速度差異，這更符合實(shí)際的地質(zhì)特征，這在一定程度上說明了聯(lián)合反演成像的優(yōu)點(diǎn)和可靠性。3-D橫波速度模型的主要特征與圖1中觀察到的主要特征基本一致，但是提供了更精細(xì)的在深度上的速度結(jié)構(gòu)分布信息。在淺層地殼中，速度分布與地表構(gòu)造單元吻合較好，龍門山斷層兩邊速度差異比較大，東邊低速異常區(qū)與四川盆地厚度大于10 km的以古生代和中生代為主的沉積層[14]有很好的對(duì)應(yīng)性，西邊高速區(qū)與厚度達(dá)10 km的中上三疊復(fù)理石[14]對(duì)應(yīng)(圖7中5 km和10 km速度切片圖，圖8)。在中地殼中，中云南塊體表現(xiàn)為高速(圖7中10 km和17.5 km速度切片圖)，驗(yàn)證了前人的研究成果[3,7,11,22]，說明了可能與地幔物質(zhì)上涌形成的峨眉山洪流玄武巖的內(nèi)核有關(guān)[17]。在下地殼，速度分布模式與淺部的正好相反，龍門山斷層的西邊是低速，東邊是高速(圖7中40 km速度切片圖，圖8)。面波和重力聯(lián)合反演的速度模型在青藏高原中下地殼存在明顯的低速區(qū)(圖8)，這與前人在該地區(qū)的很多研究一致[18,31,35,46-47]。下地殼龍門山東邊的高速(圖7中40 km速度切片圖)反映了四川盆地堅(jiān)硬的基巖，說明四川盆地是穩(wěn)定的塊體，這種堅(jiān)硬的塊體在青藏高原的東邊形成一個(gè)“墻”，對(duì)青藏高原地下物質(zhì)向東遷移造成一定的阻擋作用，從而在此積累應(yīng)力，導(dǎo)致了汶川地震和蘆山地震的發(fā)生。

圖5 單獨(dú)反演和聯(lián)合反演在不同深度模型的恢復(fù)度測(cè)試結(jié)果Fig.5 Comparison of resolution restoration test of separate inversion and joint inversion

圖6 單獨(dú)反演和聯(lián)合反演在不同深度恢復(fù)度測(cè)試的相似度結(jié)果Fig.6 Comparison of resolution semblance of separate inversion and joint inversion

圖7 不同深度的單獨(dú)反演和聯(lián)合反演的Vs切片對(duì)比Fig.7 The Vs model comparison of surface-wave only and joint inversion at different depth

圖8 沿著AA’，BB’，CC’和DD’的Vs模型剖面圖Fig.8 Four Vs profiles from the Vs model by joint inversion

4 結(jié)論

本研究收集了四川地區(qū)的瑞利波相速度頻散數(shù)據(jù)和重力資料，并對(duì)該地區(qū)進(jìn)行了地震面波-布格重力聯(lián)合成像反演，獲得了該地區(qū)的速度-密度自洽的巖石圈三維S波速度結(jié)構(gòu)，聯(lián)合反演得到的模型能夠比較好地?cái)M合面波頻散和重力異常數(shù)據(jù)。比較兩種模型可以看出，面波反演和聯(lián)合反演的模型在大的特征上比較相似，但在一些細(xì)節(jié)上面波和重力聯(lián)合反演給出了更好地約束。模型中的速度分布差異的強(qiáng)烈性反映了該地區(qū)強(qiáng)的結(jié)構(gòu)不均勻性。

在上地殼，速度分布與地表的地質(zhì)特征基本吻合，四川盆地低速反映了厚的沉積層，青藏高原東邊高速區(qū)與厚度達(dá)10 km的中上三疊復(fù)理石有關(guān)。青藏高原東部中下地殼存在明顯的低速區(qū)，并以四川盆地為界，這反映了青藏高原東南邊緣存在中下地殼流。本研究得到的新模型是有一定的地質(zhì)意義，可以為進(jìn)一步研究該區(qū)域提供重要的參考資料。

參考文獻(xiàn)：

[1] 王衛(wèi)民，趙連鋒，李娟，等．四川汶川 8．0 級(jí)地震震源過程[J]．地球物理學(xué)報(bào)，2008，51(5)： 1403-1410.

WANG W M, ZHAO L F, LI J,et al．Rupture process of the Ms 8.0 Wenchuan earthquake of Sichuan,China[J]. Chinese J.Geophys, 2008,51(5):1403-1410. (In Chinese)

[2] 張竹琪, 張培震, 王慶良. 龍門山高傾角逆斷層結(jié)構(gòu)與孕震機(jī)制[J].地球物理學(xué)，2010， 53 (9): 2068-2082.

ZHANG Z Q, ZHANG P Z , WANG Q L.The structure and seismogenic mechanism of Longmenshan high dip-angle reverse fault[J]. Chinese J.Geophys,2010, 53 (9): 2068-2082. (In Chinese)

[3] 王椿鏞，吳建平，樓海,等．川西藏東地區(qū)的地殼 P 波速度結(jié)構(gòu)[J]．中國科學(xué)(D 輯)， 2003,33(增刊):181-189.

WANG C Y, WU J P, LOU H, et al. P-wave crustal structure in western Sichuan and eastern Tibetan region[J]. Science in China, Series D, 2003,33 (supp.): 181-189. (In Chinese)

[4] 吳建平，黃媛，張?zhí)熘?等.汶川8.0級(jí)地震余震分布及周邊區(qū)域P波三維速度結(jié)構(gòu)研究[J]. 地球物理學(xué)報(bào),2009,51 (2)：320-328.

WU J P, HUANG Y, ZHANG T Z, et al. Aftershock distribution of the Ms8.0 Wenchuan earthquake and 3-D P-wave velocity structure in and around source region[J]. Chinese Journal of Geophysics, 2009,52(1):102-111. (In Chinese)

[5] 房立華，吳建平，王未來，等.四川蘆山Ms7.0級(jí)地震及其余震序列重定位[J].科學(xué)通報(bào),2013(20)：1901-1909.

FANG L H, WU J P, WANG W L, et al. Relocation of the mainshock and aftershock sequences of Ms7.0 Sichuan Lushan earthquake[J]. Chin Sci Bull,2013(20):1901-1909. (In Chinese)

[6] 李昱，姚華建，劉啟元,等．川西地區(qū)臺(tái)陣環(huán)境噪聲瑞利波相速度層析成像[J]．地球物理學(xué)報(bào)，2010, 53(4)：842-852.

LI YU,YAO H J,LIU Q Y,et al.Phase velocity array tomography of Rayleigh waves in western Sichuan from ambient seismic noise[J]. Chinese J.Geophys, 2010, 53(4)：842-852. (In Chinese)

[7] 吳建平，明躍紅，王椿鏞. 云南數(shù)字地震臺(tái)站下方的S波速度結(jié)構(gòu)研究[J]. 地球物理學(xué)報(bào)，2001,44(2)：228-230.

WU J P, MING Y H, WANG C Y. S wave velocity structure beneath digital seismic stations of Yunnan Province inferred from teleseimic receiver function modeling[J]. Chinese J. of Geophy, 2001,44(2): 223-232. (In Chinese)

[8] 郭飚，劉啟元，陳九輝，等．川西龍門山及鄰區(qū)地殼上地幔遠(yuǎn)震P波層析成像[J]．地球物理學(xué)報(bào)，2009,52(2)：346-355.

GUO B, LIU Q Y, CHEN J H, et al. Teleseismic P-wave tomography of the crust and upper mantle in Longmenshan area, west Sichuan[J]. Chinese J.Geophys, 2009,52(2)：346-355. (In Chinese)

[9] 胥頤，黃潤(rùn)秋，李志偉，等. 龍門山構(gòu)造帶及汶川震源區(qū)的S波速度結(jié)構(gòu)[J]. 地球物理學(xué)報(bào)，2009, 52 (2)： 329-338.

XU Y, HUANG R Q, LI Z W, et al.S-wave velocity structure of the Longmen Shan and Wenchuan earthquake area[J]. Chinese J.Geophys, 2009, 52 (2)： 329-338. (In Chinese)

[10] AKI, K., P. G. RICHARDS. Quantitative Seismology: Theory and Methods[M].San Francisco：W. H. Freeman, 1980.

[11] BAO X, SONG X, XU M, et al. Crust and upper mantle structure of the North China Craton and the NE Tibetan Plateau and its tectonic implications[J]. Earth and Planetary Science Letters, 2013,369:129-137.

[12] BIRCH, F. The velocity of compressional waves in rocks to 10 kilobars, part 2[J]. J. Geophys. Res., 1961,66:2199-2224.

[13] BROCHER, T. M. Empirical relations between elastic wavespeeds and density in the Earth’s crust, Bull[J]. Seismol. Soc. Am., 2005,95:2081-2092.

[14] BURCHFIEL BC, ROYDEN LH, VAN DER HILST RD, et al. A geological and geophysical context for the Wenchuan earthquake of 12 May 2008, Sichuan, People's Republic of China[J]. GSA Today,2008, 18: 4-11.

[15] BAI D, UNSWORTH MJ, MEJU MA, et al. Crustal deformation of the eastern Tibetan plateau revealed by magnetotelluric imaging[J]. Nat Geosci,2010, 3(5):358-362.

[16] FANG HONGJIAN, HAIJIANG ZHANG, HUAJIAN YAO, et al. A new three-dimensional joint inversion algorithm of body-wave and surface-wave data and its application to the Southern California Plate Boundary Region[J]. Journal of Geophysical Research-Solid-Earth,2016,121:101-114.

[17] HE B, XU Y, CHUNG S, et al.Sedimentary evidence for a rapid, kilometer-scale crustal doming prior to the eruption of the Emeishan flood basalts[J]. Earth Planet Sci Lett,2003, 213(3):391-405.

[18] HU S, HE L, WANG J. Heat flow in the continental area of China: a new data set[J]. Earth Planet Sci Lett ,2000,179(2):407-419.

[19] HE R X, SHANG C, YU H, et al. Detailed moho depth mapping of continental China by receiver function analysis[J]. Geophysical Journal International,2014,199, 1910-1918.

[20] LINES, L. R., A. K. SCHULTZ, AND S. TREITEL. Cooperative inversion of geophysical data[J].Geophysics, 1988,53:8-20.

[21] LEES, J. M., AND J. C. VANDECAR. Seismic tomography constrained by Bouguer gravity anomalies: Applications in western Washington, Pure Appl[J]. Geophys., 1991,135:31-52.

[22] LIU Q, VAN DER HILST RD, LI Y, et al. Eastward expansion of the Tibetan Plateau by crustal flow and strain partitioning across faults[J]. Nat Geosci ,2014,7(5):361-365.

[23] MACEIRA, M., S. R. TAYLOR, C. J. AMMON, et al. High-resolution Rayleigh wave slowness tomography of central Asia[J].J. Geophys. Res., 2005,110,B06304:2106-2122.

[24] MACEIRA, M., AND C. J. AMMON. Joint inversion of surface wave velocity and gravity observations and its application to central Asian basins shear velocity structure[J].J. Geophys. Res.,2009,114,B02314:2156-2202.

[25] NAFE, J. E., C. L. DRAKE. Physical properties of marine sediments[J].In The Sea, 1963(3)：794-815.

[26] OBREBSKI, MATHIAS, RICHARD M. ALLEN, et al. Shear wave tomography of China using joint inversion of body and surface wave constraints[J]. Journal of Geophysical Research: Solid Earth,2012, 117:0148-0227.

[27] PEI S, SU J, ZHANG H, et al. Three-dimensional seismic velocity structure across the 2008 Wenchuan Ms 8.0 earthquake, Sichuan, China[J]. Tectonophysics,2010, 491(1):211-217.

[28] PEI S, ZHANG H, SU J, et al. Ductile gap between the Wenchuan and Lushan earthquakes revealed from the two-dimensional Pg seismic tomography[J]. Scientific reports, 2014(4):6489:1-6.

[29] PAIGE CC, SAUNDERS MA. LSQR: sparse linear equations and least squares problems[J]. ACM Trans math Soft,1982, 8(2):195-209.

[30] PLOUFF, D.Gravity and magnetic fields of polygonal prisms and application to magnetic terrain corrections[J].Geophysics, 1976,41:727-741.

[31] ROYDEN LH, BURCHFIEL BC, KING RW,et al. Surface deformation and lower crustal flow in eastern Tibet[J]. Science,1997, 276(5313):788-790.

[32] ROYDEN LH, BURCHFIEL BC, VAN DER HILST RD.The geological evolution of the Tibetan Plateau[J]. Science,2008, 321(5892):1054-1058.

[33] SHAPIRO, N.M., CAMPILLO, M., STEHLY, et al. High-resolution surface-wave tomography from ambient seismic noise[J]. Science,2005(307):1615-1618.

[34] SHEN, W., RITZWOLLER, M.H., KANG, D., et al. A seismic reference model for the crust and uppermost mantle beneath China from surface wave dispersion[J]. Geophysical Journal International, 2018,206(2):954-979.

[35] SUN X, BAO X, XU M, et al.Crustal structure beneath SE Tibet from joint analysis of receiver functions and Rayleigh wave dispersion[J]. Geophys Res Lett,2014,41(5):1479-1484.

[36] Saito, M., DISPER80: A subroutine package for the cal- culation of seismic normal mode solutions, in D. J. Doornbos (ed.), Seismological Algorithms: Computational Methods and Computer Programs, Academic Press[J].New York, 1988:293-319.

[37] TAKEUCHI, H., M. SAITO, Seismic surface waves, in Methods of Computational Physics, edited by B. A. Bolt[J].Academic, New York,1972:217-295.

[38] TAPLEY, B..GGM02—An improved Earth gravity field model from GRACE[J]. J. Geod., 2005,79:467-478.

[39] WANG M, HUBBARD J, PLESCH A, et al.Three‐dimensional seismic velocity structure in the Sichuan basin, China[J]. J Geophys Res Solid Earth,2016,121:1007-1022.

[40] WANG Z, FUKAO Y, PEI S. Structural control of rupturing of the Mw7. 9 2008 Wenchuan Earthquake, China[J]. Earth Planet Sci Lett ,2009,279(1):131-138.

[41] WANG Z, HUANG R, PEI S.Crustal deformation along the Longmen-Shan fault zone and its implications for seismogenesis[J]. Tectonophysics，2014，610:128-137.

[42] WANG Z, SU J, LIU C, et al.New insights into the generation of the 2013 Lushan Earthquake (Ms 7.0), China[J]. J Geophys Res Solid Earth ,2015,120: 3507-3526.

[43] XU L, RONDENAY S, VAN DER HILST RD. Structure of the crust beneath the southeastern Tibetan Plateau from teleseismic receiver functions[J]. Phys Earth Planet Inter ,2007,165(3):176-193.

[44] YAO H, BEGHEIN C, VAN DER HILST RD. Surface wave array tomography in SE Tibet from ambient seismic noise and two-station analysis-II. Crustal and upper-mantle structure[J].Geophys J Int,2008，173(1):205-219.

[45] YAO H, VAN DER HILST R D, MONTAGNER J P.Heterogeneity and anisotropy of the lithosphere of SE Tibet from surface wave array tomography[J]. J GeophysRes,2010, 115:1016-1125.

[46] YAO H.Lithospheric structure and deformation in SE Tibet revealed by ambient noise and earthquake surface wave tomography: Recent advances and perspectives[J]. Earthq Sci,2012(25):371-383.

[47] YANG Y, RITZWOLLER MH, ZHENG Y, et al. A synoptic view of the distribution and connectivity of the mid-crustal low velocity zone beneath Tibet[J]. J Geophys Res,2012, 117,B04303:2156-2202.

[48] ZHANG H, MACEIRA M, ROUX P, et al.Joint inversion of body-wave arrival times and surface-wave dispersion for three-dimensional seismic structure around SAFOD[J]. Pure Appl Geophys,2014,171(11):3013-3022.

[49] ZHAO G, UNSWORTH MJ, ZHAN Y,et al. Crustal structure and rheology of the Longmenshan and Wenchuan Mw 7.9 earthquake epicentral area from magnetotelluric data[J]. Geology ,2012,40(12):1139-1142.

[50] ZHAO Z, FANG J, ZHEN S .Crustal velocity and relocated events in Longmen Shan fault zone[J]. Acta Seismol Sin, 1997(19):615-622.

[51] ZHENG, Y.F., XIAO, W.J., ZHAO, G. Introduction to tectonics of China[J]. Gondwana Res,2013(23):1189-1206.

[52] ZHENG, Y., W. SHEN, L. ZHOU, Y. YANG, et al. Crust and uppermost mantle beneath the North China Craton, northeastern China, and the Sea of Japan from ambient noise tomography[J].J. Geophys. Res.,2011,116, B12312:2156-2202.

[53] ZHOU, L., J. XIE, W. SHEN, et al. The structure of the crust and uppermost mantle beneath South China from ambient noise and earthquake tomography[J].Geophysical Journal International, 2012,189(3):1565-1583.