利用重力和面波聯(lián)合反演確定四川地區(qū)三維巖石圈速度結(jié)構(gòu)

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

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

0 引言

大約45 Ma前，印度-歐亞大陸的碰撞造成青藏高原的隆起和抬升[31-32]，印度-歐亞大陸之間的連續(xù)匯聚導(dǎo)致碰撞帶地殼發(fā)生了不小于1 500 km的縮短。如何調(diào)節(jié)或吸收如此大規(guī)模的匯聚量，一直是地球科學(xué)界爭論的熱點。而四川盆地，西邊與青藏高原東部邊界龍門山斷裂帶相接，東邊毗鄰揚子克拉通西緣，因此四川地區(qū)既是檢驗不同補(bǔ)償模式和理解陸內(nèi)形變機(jī)制最為理想的研究“窗口”，同時也是探索高原側(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]，因此對四川地區(qū)的結(jié)構(gòu)進(jìn)行成像有進(jìn)一步研究的必要。

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

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

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

1.1 數(shù)據(jù)

本研究采用的面波數(shù)據(jù)是8 s～70 s周期范圍內(nèi)的瑞利波相速度，該組面波頻散數(shù)據(jù)來自Shen等[34]的最新面波層析成像研究。他們使用收集了中國多個密集臺陣(CEArray, China Array, NECESS, PASSCAL, GSN)及其周邊地區(qū)(Korean Seismic Network, F-Net, KNET)總共2 000多個地震臺站的波形數(shù)據(jù)，進(jìn)行了背景噪聲面波層析成像，并且整合了華南、東北部分地區(qū)已有的地震面波層析成像的結(jié)果[52-53]，得到了整個區(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波速度的變化情況；長周期(50 s～70 s)的相速度分布圖主要反映了上地幔頂部S波速度的橫向變化。

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

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

為了能夠?qū)崿F(xiàn)聯(lián)合反演，面波反演方程組與重力反演方程組需要統(tǒng)一到一個系統(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)基于面波單獨反演的模型預(yù)測的重力異常分布；(d)基于聯(lián)合反演的Vs模型預(yù)測的重力異常分布

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

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

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

1.3 重力異常的理論計算和近似

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

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

以表示為：

(2)

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

yiln(Rijk+xi)]

(3)

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

(4)

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

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

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

α=3.125ρ-2.40

(5)

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

α=6.86-7.55ρ+2.64ρ2

(6)

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

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

(7)

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

(8)

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

(9)

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

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

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

(10)

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

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

3 結(jié)果和討論

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

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

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

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

圖7 不同深度的單獨反演和聯(lián)合反演的Vs切片對比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ù)和重力資料，并對該地區(qū)進(jìn)行了地震面波-布格重力聯(lián)合成像反演，獲得了該地區(qū)的速度-密度自洽的巖石圈三維S波速度結(jié)構(gòu)，聯(lián)合反演得到的模型能夠比較好地擬合面波頻散和重力異常數(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]．地球物理學(xué)報，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級地震余震分布及周邊區(qū)域P波三維速度結(jié)構(gòu)研究[J]. 地球物理學(xué)報,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].科學(xué)通報,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ū)臺陣環(huán)境噪聲瑞利波相速度層析成像[J]．地球物理學(xué)報，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ù)字地震臺站下方的S波速度結(jié)構(gòu)研究[J]. 地球物理學(xué)報，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é)報，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] 胥頤，黃潤秋，李志偉，等. 龍門山構(gòu)造帶及汶川震源區(qū)的S波速度結(jié)構(gòu)[J]. 地球物理學(xué)報，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.