天山地區(qū)現(xiàn)今地殼形變及構造應力的三維數(shù)值模擬*

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

(中國地震局地震大地測量重點實驗室，湖北 武漢 430071)

0 引言

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

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

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

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

1 三維有限元模型的建立

1.1 幾何模型構建

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

圖1 天山地區(qū)GPS測點分布與M≥5.0地震活動特征

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

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

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

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

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

考慮到斷層活動在地殼運動變形過程中起到至關重要的作用(鄭勇等，2007)，除邊界約束外，還需要對模型中的斷層進行約束。Pylith軟件通過自動創(chuàng)建“內(nèi)聚單元”實現(xiàn)斷層兩側頂點的相對運動，從而控制斷層滑動，由此可以創(chuàng)建包含斷層的不連續(xù)體模型。參考牛之俊等(2007)和劉代芹等(2016)給出的南、北天山斷裂帶的閉鎖深度及滑動速率結果，經(jīng)過網(wǎng)格搜索法實驗，計算地表GPS站速度模擬值與觀測值間的加權均方根誤差，最終確定南天山斷裂帶的閉鎖深度為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流變包含與時間無關的彈性行為和與時間有關的黏性行為。Pylith軟件由介質(zhì)密度、橫波速度和縱波速度定義彈性，由黏滯系數(shù)定義黏性。 Crust 1.0為分辨率1°×1°的全球地殼模型，將全球劃分為64 800個單元，垂向上分為9層，每層均給出P波和S波速度以及密度參數(shù)。石耀霖和曹建玲(2008)對實驗室流變實驗結果應用于估算巖石圈等效黏滯系數(shù)中的多種影響因素進行了討論，并以溫度和應變速率的研究成果為基礎，對中國大陸地殼和上地幔等效黏滯系數(shù)做出了較為詳細的估計。本文參考這些結果，確定研究區(qū)介質(zhì)的物性參數(shù)。考慮到中地殼及下地殼的黏滯系數(shù)隨深度增加而減小，在確定介質(zhì)分界處及底面深度處的黏滯系數(shù)后，規(guī)定其余深度黏滯系數(shù)由此進行線性插值得到，以盡量實現(xiàn)與真實情況的吻合。表1展示了非均勻分層模型的物性參數(shù)。

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

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

2 模擬速度場結果及其特征

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

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

表2 3種地表GPS站模型模擬速度與觀測速度對比

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

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

3 模擬應力率場結果及其特征

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

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

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

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

4 討論

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

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

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

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

5 結論

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

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

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

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