孫明志，譚 凱，魯小飛，張彩紅，李 琦，劉志軍

(中國(guó)地震局地震大地測(cè)量重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430071)

0 引言

天山山系綿延于中亞腹地，全長(zhǎng)達(dá)2 500 km，南北寬250～400 km，是世界上最大、最活躍的陸內(nèi)造山帶(鄧起東等，2000；Vinnik，2002)。新生代以來(lái)，在印度與歐亞板塊的碰撞及兩側(cè)盆地楔入背景下，天山山體復(fù)活隆升，兩側(cè)盆地凹陷，發(fā)生了強(qiáng)烈的地殼變形，發(fā)育了一系列斷裂帶，控制了一系列強(qiáng)震(Molnar，Tapponnier，1975；張培震等，1996)。研究天山地區(qū)的地殼形變與應(yīng)力分布特征，對(duì)理解該地區(qū)構(gòu)造變形的運(yùn)動(dòng)學(xué)機(jī)制及地震成因，分析未來(lái)潛在地震危險(xiǎn)性有重要的科學(xué)意義。

大量學(xué)者利用GPS等大地測(cè)量資料，采用插值方法、彈性位錯(cuò)模型等對(duì)天山地區(qū)的構(gòu)造運(yùn)動(dòng)特征進(jìn)行了研究，得到了一系列十分有價(jià)值的研究成果(王琪等，2000；牛之俊等，2007；王曉強(qiáng)等，2007；楊少敏等，2008；王偉等，2014)。隨著計(jì)算機(jī)技術(shù)和數(shù)值計(jì)算理論的發(fā)展，數(shù)值模擬逐漸成為了研究地學(xué)問題的重要手段。與傳統(tǒng)研究方法相比，數(shù)值模擬體現(xiàn)出一定的計(jì)算優(yōu)勢(shì)。有限元數(shù)值模擬以較為客觀的模型架構(gòu)和各種地質(zhì)、地震及大地測(cè)量資料為基礎(chǔ)，可以兼顧整個(gè)區(qū)域橫向及縱向上物質(zhì)的非均勻性，由此可以得到更為準(zhǔn)確、細(xì)致的構(gòu)造運(yùn)動(dòng)特征(劉峽，2007)。

近年來(lái)不少學(xué)者對(duì)天山地區(qū)的構(gòu)造變形特征進(jìn)行了數(shù)值模擬研究，如Qiu等(2009)利用有限元模擬西南天山與帕米爾—西昆侖之間的匯聚趨勢(shì)，得出了匯聚邊界向北、向東南延伸的結(jié)論，但其模型沒有考慮斷裂的影響；雷顯權(quán)等(2011)采用非連續(xù)接觸分析方法模擬天山地區(qū)的地殼運(yùn)動(dòng)變形，發(fā)現(xiàn)地殼運(yùn)動(dòng)速度由西向東、由南向北逐漸減小，然而其建立的天山中段二維彈性薄片模型，不能考慮復(fù)雜介質(zhì)分布情況，在空間上具有一定的局限性；王子韜等(2020)通過有限元計(jì)算了天山區(qū)域主要斷裂帶彈性應(yīng)變能密度增加速率及應(yīng)力積累速率，發(fā)現(xiàn)天山的斷層和地震活動(dòng)性主要受控于近南北向的主壓應(yīng)力，其建立的二維平面彈性模型將斷層處理為連續(xù)變形的軟弱帶，不能很好地體現(xiàn)斷裂帶的滑動(dòng)特征。

物理模型與地球?qū)嶋H情況的符合程度和精細(xì)化程度可能會(huì)幫助我們提高對(duì)地殼運(yùn)動(dòng)和地震活動(dòng)的了解。因此，本文納入天山地區(qū)主要活動(dòng)斷裂，根據(jù)天山地區(qū)的地塊構(gòu)造、活動(dòng)斷裂帶分布和流變結(jié)構(gòu)等資料，建立了1種彈性模型和2種黏彈性模型，共3種不同的三維塊體運(yùn)動(dòng)斷層位錯(cuò)模型。在GPS地殼運(yùn)動(dòng)觀測(cè)結(jié)果的約束下，附加斷層滑動(dòng)速率和閉鎖深度約束，使用有限元程序Pylith(Aagaard，2017)模擬天山地區(qū)的地殼運(yùn)動(dòng)變形，討論巖石圈縱向分層和黏彈性效應(yīng)對(duì)模擬地表速度場(chǎng)結(jié)果的影響，計(jì)算獲得研究區(qū)地殼速度場(chǎng)及應(yīng)力場(chǎng)，并分析其量值大小、分布形態(tài)及特征。

1 三維有限元模型的建立

1.1 幾何模型構(gòu)建

考慮到天山地區(qū)的GPS測(cè)點(diǎn)分布及地震活動(dòng)性(圖1)，將研究區(qū)選取為(75°～94°E，38°～47°N)，包含了中國(guó)境內(nèi)天山及境外天山的一部分，區(qū)域面積約為1 565 124 km。大陸板塊內(nèi)部以塊體運(yùn)動(dòng)為特征(鄧起東等，2002)。本文參考張培震等(2003)對(duì)西域活動(dòng)地塊區(qū)的劃分，將研究區(qū)劃分為塔里木、天山和準(zhǔn)噶爾3個(gè)活動(dòng)塊體。由一系列斷裂耦合成的南天山山前斷裂帶和北天山山前斷裂帶為研究區(qū)內(nèi)的重要活動(dòng)構(gòu)造帶，作為劃分塊體的邊界帶。

圖1 天山地區(qū)GPS測(cè)點(diǎn)分布與M≥5.0地震活動(dòng)特征

模型厚度取為80 km，將縱向上分為上地殼、中地殼、下地殼和巖石圈上地幔4層介質(zhì)，以構(gòu)建均勻分層模型和非均勻分層模型。均勻分層模型4層介質(zhì)的厚度均為20 km，非均勻分層模型每層厚度依據(jù)全球地殼模型Crust 1.0(Laske，2012)，由每個(gè)塊體內(nèi)的不同單元的各層厚度取均值確定。由于納入模型的斷裂的傾角均較大(馮先岳，1986)，故將其作一定簡(jiǎn)化，固定為70°。另外，研究區(qū)的面積并不太大，可以將研究區(qū)域的表面作為平面處理。將關(guān)鍵點(diǎn)經(jīng)緯度坐標(biāo)投影轉(zhuǎn)換成笛卡爾坐標(biāo)后使用Coreform Trelis軟件建立幾何模型并進(jìn)行網(wǎng)格劃分，網(wǎng)格劃分時(shí)斷層面格網(wǎng)邊長(zhǎng)設(shè)為5 km，剩余部分按照距斷層面的距離以1.07的漸變梯度進(jìn)行劃分。均勻分層模型包含985 052個(gè)四面體單元，177 174個(gè)節(jié)點(diǎn)；非均勻分層模型包含1 121 486個(gè)四面體單元，201 561個(gè)節(jié)點(diǎn)。圖2展示了三維非均勻分層模型的幾何圖形及網(wǎng)格劃分。

圖2 三維非均勻分層模型的幾何圖形及格網(wǎng)劃分

1.2 模型約束條件與物性參數(shù)

GPS觀測(cè)得到的天山地區(qū)現(xiàn)今地殼運(yùn)動(dòng)結(jié)果，可為有限元模型提供可靠的運(yùn)動(dòng)學(xué)邊界條件約束。Wang和Shen(2020)收集和處理了多個(gè)來(lái)源的數(shù)據(jù)，并整合以往研究中的速度解，得到了統(tǒng)一參考框架下的中國(guó)大陸及周邊區(qū)域GPS觀測(cè)結(jié)果，這些數(shù)據(jù)精度較高、密度較大，還給出了中國(guó)境外天山西部的GPS速度結(jié)果。因此，本文采用這些可靠且豐富的數(shù)據(jù)作為模型的邊界條件。首先，選取靠近邊界且不處于斷層處的GPS測(cè)點(diǎn)，將其GPS速度值賦給離站點(diǎn)最近的模型地表側(cè)邊邊界節(jié)點(diǎn)。然后，對(duì)已獲得速度值的節(jié)點(diǎn)進(jìn)行插值，并規(guī)定邊界速度不隨深度改變，進(jìn)而獲取剩余所有側(cè)面節(jié)點(diǎn)的速度值(圖3)。對(duì)于模型底面，規(guī)定其在垂向上的位移為零，水平向不作約束。最終，有限元模型上表面為自由表面，側(cè)面在水平向上由GPS觀測(cè)結(jié)果約束而在垂向上自由，下表面在垂向上固定而在水平向上自由。

圖3 天山地區(qū)有限元模型邊界條件

考慮到斷層活動(dòng)在地殼運(yùn)動(dòng)變形過程中起到至關(guān)重要的作用(鄭勇等，2007)，除邊界約束外，還需要對(duì)模型中的斷層進(jìn)行約束。Pylith軟件通過自動(dòng)創(chuàng)建“內(nèi)聚單元”實(shí)現(xiàn)斷層兩側(cè)頂點(diǎn)的相對(duì)運(yùn)動(dòng)，從而控制斷層滑動(dòng)，由此可以創(chuàng)建包含斷層的不連續(xù)體模型。參考牛之俊等(2007)和劉代芹等(2016)給出的南、北天山斷裂帶的閉鎖深度及滑動(dòng)速率結(jié)果，經(jīng)過網(wǎng)格搜索法實(shí)驗(yàn)，計(jì)算地表GPS站速度模擬值與觀測(cè)值間的加權(quán)均方根誤差，最終確定南天山斷裂帶的閉鎖深度為17 km，左旋走滑速率約束為2.8 mm/a，85°E以西區(qū)域擠壓速率約束為8.8 mm/a，85°E以東區(qū)域擠壓速率約束為3.5 mm/a；北天山斷裂帶的閉鎖深度確定為16 km，左旋走滑速率約束為2.1 mm/a，擠壓速率約束為3.3 mm/a。

Maxwell流變包含與時(shí)間無(wú)關(guān)的彈性行為和與時(shí)間有關(guān)的黏性行為。Pylith軟件由介質(zhì)密度、橫波速度和縱波速度定義彈性，由黏滯系數(shù)定義黏性。 Crust 1.0為分辨率1°×1°的全球地殼模型，將全球劃分為64 800個(gè)單元，垂向上分為9層，每層均給出P波和S波速度以及密度參數(shù)。石耀霖和曹建玲(2008)對(duì)實(shí)驗(yàn)室流變實(shí)驗(yàn)結(jié)果應(yīng)用于估算巖石圈等效黏滯系數(shù)中的多種影響因素進(jìn)行了討論，并以溫度和應(yīng)變速率的研究成果為基礎(chǔ)，對(duì)中國(guó)大陸地殼和上地幔等效黏滯系數(shù)做出了較為詳細(xì)的估計(jì)。本文參考這些結(jié)果，確定研究區(qū)介質(zhì)的物性參數(shù)?？紤]到中地殼及下地殼的黏滯系數(shù)隨深度增加而減小，在確定介質(zhì)分界處及底面深度處的黏滯系數(shù)后，規(guī)定其余深度黏滯系數(shù)由此進(jìn)行線性插值得到，以盡量實(shí)現(xiàn)與真實(shí)情況的吻合。表1展示了非均勻分層模型的物性參數(shù)。

表1 非均勻分層模型的物性參數(shù)

至此，通過建立的幾何模型、給定的約束條件和物性參數(shù)，經(jīng)有限元數(shù)值計(jì)算后，即可得到模擬結(jié)果。模型通過指定運(yùn)動(dòng)學(xué)上的邊界條件和斷層滑移速率模擬應(yīng)力絕對(duì)狀態(tài)下的擾動(dòng)，因此忽略了重力的影響。

2 模擬速度場(chǎng)結(jié)果及其特征

為研究巖石圈流變物質(zhì)黏彈性效應(yīng)對(duì)地表形變的影響，本文給均勻分層幾何模型分別賦予彈性參數(shù)和黏彈性參數(shù)，建立均勻分層彈性模型(M1)和均勻分層黏彈性模型(M2)。兩種模型唯一的區(qū)別為彈性模型將整個(gè)域視為彈性，黏彈性模型將上地殼視為彈性，中、下地殼和上地幔視為黏彈性。黏彈性模型與彈性模型對(duì)應(yīng)的彈性參數(shù)相同。另外，為非均勻分層幾何模型賦予黏彈性參數(shù)，建立非均勻分層黏彈性模型(M3)，以討論縱向介質(zhì)分層對(duì)地表形變的影響。最后，根據(jù)表1設(shè)置的黏滯系數(shù)，估計(jì)出黏彈性介質(zhì)的最大弛豫時(shí)間。3種模型的數(shù)值計(jì)算采用500 a為1個(gè)時(shí)間步，經(jīng)過100個(gè)時(shí)間步的計(jì)算，最終得到5萬(wàn)年后的地殼形變計(jì)算結(jié)果。

首先，由最后1個(gè)時(shí)間步的模型計(jì)算結(jié)果計(jì)算震間速度場(chǎng)。然后，從模擬速度場(chǎng)結(jié)果中提取模型所含地表中424個(gè)GPS站的速度值，與實(shí)際觀測(cè)的GPS數(shù)據(jù)進(jìn)行比較，結(jié)果見表2。在3種模型中，超過95%的站點(diǎn)速度模擬結(jié)果與觀測(cè)結(jié)果差異在4 mm/a以內(nèi)，速度方位模擬結(jié)果與觀測(cè)結(jié)果差異在20°以內(nèi)，且都基本服從正態(tài)分布。由此可見，3種模型均能較好地模擬出地表速度場(chǎng)，非均勻分層模型的模擬結(jié)果略優(yōu)于均勻分層的模擬結(jié)果，說明介質(zhì)厚度的改變對(duì)本文模擬地表速度場(chǎng)的影響并不大。此外，相較于純彈性模型，黏彈性模型得到了更優(yōu)的模擬結(jié)果。下面僅對(duì)擬合地表速度場(chǎng)效果最優(yōu)的非均勻分層黏彈性模型(以下簡(jiǎn)稱M3模型)的模擬結(jié)果展開分析。

表2 3種地表GPS站模型模擬速度與觀測(cè)速度對(duì)比

從使用M3模型得到的GPS站點(diǎn)震間速度場(chǎng)與實(shí)測(cè)速度場(chǎng)的對(duì)比圖(圖4a)及速度殘差圖(圖4b)可以看出，兩者符合程度較好。圖4c是使用M3模型給出的地表速度場(chǎng)。從大小來(lái)看，天山地區(qū)形變速率呈現(xiàn)由南向北、由西向東逐漸減小的態(tài)勢(shì)。這是因?yàn)樵搮^(qū)域從南向北、從西向東逐漸遠(yuǎn)離碰撞板塊邊界，板塊的推擠效應(yīng)逐步減弱。研究區(qū)運(yùn)動(dòng)速率最大位于模型西南角，即帕米爾—西昆侖弧形斷裂以南區(qū)域，達(dá)20 mm/a。在模型所含區(qū)域內(nèi)，南天山匯聚速率為12～15 mm/a，中天山匯聚速率為6～9 mm/a，北天山匯聚速率為4～5 mm/a。從方向來(lái)看，由研究區(qū)西南部的NNW向逐漸向東、向北過渡到NNE及NE向，從塔里木地區(qū)到準(zhǔn)噶爾地區(qū)速度由順時(shí)針方向轉(zhuǎn)為逆時(shí)針方向，驗(yàn)證了塔里木地塊的順時(shí)針旋轉(zhuǎn)運(yùn)動(dòng)。

圖4 天山地區(qū)使用M3模型模擬速度場(chǎng)

3 模擬應(yīng)力率場(chǎng)結(jié)果及其特征

現(xiàn)今構(gòu)造應(yīng)力場(chǎng)是現(xiàn)今地殼變形和構(gòu)造運(yùn)動(dòng)的動(dòng)力因素，是地震活動(dòng)的內(nèi)在直接原因之一(鄧起東等，2000)。模擬得出的應(yīng)力實(shí)際上是應(yīng)力的平均年增量，即應(yīng)力率，它反映了現(xiàn)今天山地區(qū)應(yīng)力場(chǎng)的基本特征。由1900—2020年天山地區(qū)的地震目錄可知，該地區(qū)大部分地震的震源深度在10～20 km。為研究天山地區(qū)地殼內(nèi)部的應(yīng)力狀態(tài)，本文提取地表以下15 km的應(yīng)力模擬結(jié)果，計(jì)算水平最大、最小主應(yīng)力率。圖5為天山地區(qū)地下15 km的水平方向主應(yīng)力率圖和最大剪應(yīng)力率云圖。

從圖5a可以看出，研究區(qū)大部分區(qū)域的主壓應(yīng)力率大于主拉應(yīng)力率。天山兩側(cè)褶皺帶的主壓應(yīng)力率大于天山內(nèi)部主壓應(yīng)力率，且呈現(xiàn)與山脈走向相正交的特征，與Chen等(2005)的研究結(jié)論相似。主壓應(yīng)力率較大的地區(qū)集中分布在天山褶皺帶上，主壓應(yīng)力率最大的區(qū)域?yàn)槟咸焐脚c帕米爾高原相接區(qū)，達(dá)4.780 kPa/a，此區(qū)域也是發(fā)生過幾次強(qiáng)震的地震活躍區(qū)，其未來(lái)地震危險(xiǎn)性值得重點(diǎn)關(guān)注。天山西段(75°～84°E)平均主壓應(yīng)力率為1.896 kPa/a，天山東段(85°～94°E)平均主壓應(yīng)力率為1.478 kPa/a。另外，塔里木地塊和準(zhǔn)噶爾地塊作為較穩(wěn)定塊體，其應(yīng)力率大致均勻分布且明顯小于天山的應(yīng)力率，平均主壓應(yīng)力率分別為0.471 kPa/a、0.235 kPa/a。圖5b中的五角星代表發(fā)生于1900—2020年、震源深度介于10～20 km、≥5地震，由此可以看出，最大剪應(yīng)力率大的地區(qū)發(fā)震次數(shù)多，地震活動(dòng)性強(qiáng)，如南、北天山西段，最大剪應(yīng)力率較大，地震頻發(fā)；最大剪應(yīng)力率小的地區(qū)地震少發(fā)，如南、北天山東段，塔里木和準(zhǔn)噶爾，最大剪應(yīng)力率小，發(fā)震的次數(shù)少。

圖5 天山地區(qū)地下15 km截面模擬應(yīng)力率場(chǎng)

研究區(qū)主壓應(yīng)力方向大多近N-S向，絕大部分最大、最小主應(yīng)力接近水平。南天山從東到西主壓應(yīng)力總體方向?yàn)镹NW向，主壓應(yīng)力軸集中在NNW15°～NNE10°，與李杰等(2012)根據(jù)GPS結(jié)果計(jì)算的南天山主壓應(yīng)變軸方向相符。北天山從東到西主壓應(yīng)力方向從NNE向過渡為NNW向，主壓應(yīng)力軸集中在NNW16°～NNE26°。

4 討論

本文使用3種模型對(duì)地表速度場(chǎng)進(jìn)行擬合，結(jié)果表明，使用M3模型計(jì)算得到的天山地區(qū)震間速度場(chǎng)與GPS觀測(cè)結(jié)果有較好的一致性，說明本文建立的有限元模型是合理的，模擬結(jié)果是可靠的。此外，非均勻分層模型的模擬結(jié)果優(yōu)于均勻分層模型模擬結(jié)果，并且黏彈性模型提高了對(duì)震間地殼形變的擬合程度，這與前人(王輝等，2007；Li，2015)得出的實(shí)驗(yàn)結(jié)果相一致。王輝等(2007)建立了川滇地區(qū)均勻分層和不均勻分層模型，探討莫霍面起伏對(duì)模擬速度場(chǎng)結(jié)果的影響，結(jié)果表明不均勻分層模型的結(jié)果較均勻分層模型的結(jié)果有局部改善，但總體改進(jìn)不多。Li等(2015)使用有限元模型來(lái)研究黏彈性對(duì)震間形變的控制作用，結(jié)果證實(shí)巖石圈流變物質(zhì)的黏彈性效應(yīng)對(duì)震間地表形變場(chǎng)有很大的影響。

具體來(lái)看，天山地區(qū)速度場(chǎng)從南向北、由西到東逐漸減小，運(yùn)動(dòng)速率最大達(dá)20 mm/a，位于模型西南角。在模型所含區(qū)域內(nèi)，南天山匯聚速率為12～15 mm/a，中天山匯聚速率為6～9 mm/a，北天山匯聚速率為4～5 mm/a，這與牛之俊等(2007)利用GPS數(shù)據(jù)得出的天山地殼縮短速率相近。塔里木地塊速度場(chǎng)方向順時(shí)針旋轉(zhuǎn)，準(zhǔn)噶爾地塊速度場(chǎng)方向逆時(shí)針旋轉(zhuǎn)。本文模擬得到的天山地區(qū)現(xiàn)今構(gòu)造應(yīng)力場(chǎng)以近南北向水平擠壓應(yīng)力為主，與龍海英等(2007)利用震源機(jī)制解給出的解釋一致。研究區(qū)主壓應(yīng)力率從南到北呈現(xiàn)不均勻分布特征，主要分布在天山南北兩側(cè)褶皺帶，最大值處于南天山與帕米爾高原及西昆侖交界地帶。主壓應(yīng)力率場(chǎng)成近N-S向分布。天山西段(75°～84°E)平均主壓應(yīng)力率較大。塔里木地塊和準(zhǔn)噶爾地塊作為較穩(wěn)定塊體，其應(yīng)力率大致均勻分布且明顯小于天山的應(yīng)力率。

本文模擬得到的研究區(qū)地殼形變速度場(chǎng)和應(yīng)力分布結(jié)果顯示，擠壓速率在天山內(nèi)部較為均勻，而在跨越斷裂帶時(shí)快速減小，表明斷裂帶起到了吸收形變的重要作用，最大剪應(yīng)力結(jié)果也表明在這些區(qū)域積累了大量能量，這些現(xiàn)象與該區(qū)域的地震活動(dòng)分布呈現(xiàn)出很好的一致性。因此，南、北天山山前斷裂帶的蠕滑對(duì)天山地區(qū)現(xiàn)今形變場(chǎng)和應(yīng)力場(chǎng)起到一定的調(diào)整作用，天山地區(qū)地殼的縮短變形主要發(fā)生在天山南北兩側(cè)褶皺帶，這影響了天山的地震活動(dòng)性分布。

歐亞板塊與印度板塊之間的碰撞擠壓是引起天山現(xiàn)今構(gòu)造變形的直接原因。本文的模型受代表帕米爾高原推擠和塔里木順時(shí)針旋轉(zhuǎn)的運(yùn)動(dòng)學(xué)邊界條件的約束，兩者共同給天山地區(qū)提供了西強(qiáng)東弱的推擠作用力。本文模擬得到了與前人研究結(jié)果一致的由西到東的差異性構(gòu)造運(yùn)動(dòng)特征。因此，天山地區(qū)的現(xiàn)今構(gòu)造變形特征是帕米爾高原向北推擠和塔里木塊體順時(shí)針旋轉(zhuǎn)運(yùn)動(dòng)共同作用的結(jié)果，但兩者誰(shuí)起到主導(dǎo)作用，還需進(jìn)一步探討。

5 結(jié)論

本文通過建立3種不同的三維有限元模型，以GPS觀測(cè)結(jié)果為邊界約束，對(duì)斷層附加滑移速率及閉鎖深度進(jìn)行約束，獲取了天山地區(qū)的形變場(chǎng)及應(yīng)力場(chǎng)，所得結(jié)論如下：

(1)模型縱向分層的不同對(duì)模擬地表速度場(chǎng)結(jié)果的影響并不大，而考慮黏彈性效應(yīng)可以提高對(duì)震間地殼形變的擬合程度。

(2)主壓應(yīng)力率最大的地區(qū)是南天山與帕米爾高原及西昆侖交界地帶，南、北天山西段的主壓應(yīng)力率也較大，這些地區(qū)未來(lái)的地震危險(xiǎn)性值得重點(diǎn)關(guān)注。天山地區(qū)主壓應(yīng)力場(chǎng)方向近N-S向，呈現(xiàn)出以近南北向水平擠壓應(yīng)力為主的特征。

(3)南、北天山山前斷裂帶的蠕滑對(duì)天山地區(qū)現(xiàn)今形變場(chǎng)和應(yīng)力場(chǎng)起到一定的調(diào)整作用。天山地區(qū)地殼的縮短變形主要發(fā)生在天山南北兩側(cè)褶皺帶，這影響了天山地區(qū)的地震活動(dòng)性分布。