董培育　胡才博　石耀霖*

1)中國地震局地震研究所， 地震大地測量重點實驗室， 武漢　430071 2)中國科學院計算地球動力學重點實驗室， 北京　100049

?

青藏高原及周邊區(qū)域地表長期變形數(shù)值模擬

董培育1，2)胡才博2)石耀霖2)*

1)中國地震局地震研究所， 地震大地測量重點實驗室， 武漢430071 2)中國科學院計算地球動力學重點實驗室， 北京100049

印度板塊向歐亞板塊俯沖擠壓， 不僅令青藏高原上地殼在擠壓作用下發(fā)生彈性變形和運動， 且青藏高原高溫高壓下的下地殼會發(fā)生柔性流動， 并對脆性的上地殼有拖曳作用， 這2種作用一起形成現(xiàn)今的高原運動變形場。這一動力學過程已得到GPS觀測資料的證實。因此在二維平面問題中僅用上地殼在邊界作用下的彈性變形解釋是不夠的， 還要考慮柔性下地殼流動對上地殼的拖曳作用。但是拖曳力作用的大小和方向不易確定， 故文中建立了二維平面彈性有限元模型， 利用加載等效體力來模擬下地殼流動對上地殼產(chǎn)生的拖曳力。以高原內(nèi)部的GPS觀測資料為約束， 利用試錯法反演出模型中關鍵點的力， 其他位置上的力則用關鍵點上的力進行雙線性插值計算。以此來反演計算出模型區(qū)域內(nèi)的柔性下地殼的差異性流動對脆性上地殼產(chǎn)生的拖曳力(節(jié)點力的形式， 單位： N)的大小和范圍， 在86°～100°E， 26°～32°N地區(qū)主要以SE向為主， 最大達到108N； 西部局部(31°～36°N， 76°～80°E)地區(qū)有較弱的W向拖曳力， 最大為107N。文中為深入研究青藏高原及周邊區(qū)域的長期地表變形動力學機制提供了1個新的思路。

有限元數(shù)值模擬青藏高原拖曳力雙線性插值

0　引言

印度板塊長期持續(xù)向N擠壓歐亞板塊， 形成了青藏高原， 使其成為目前地球上最年輕、 形成過程最復雜的高原， 令其一直處于地球科學研究和討論的中心和熱點(許志琴等， 2011)。對于其隆升和運動變形機制， Tapponnier等(1982)認為高原被大型走滑斷裂帶切割， 在印度板塊推擠下， 各塊體沿著斷裂帶被擠出。England等(1982， 1983)利用薄板流變模型數(shù)值模擬印度板塊與歐亞板塊的碰撞形變， 且在一定區(qū)域中定量描述了高原地殼厚度和應變率的時空分布。傅容珊等(1999； 2000a， b)將大陸巖石層視為冪指數(shù)律控制的薄層， 上覆在黏滯性較低的軟流層之上蠕變流動； 數(shù)值模擬結果表明， 青藏高原隆升的主要動力來源于印度板塊對歐亞板塊的碰撞擠壓。許志琴等(1996)提出青藏高原腹地隆升的主要原因與大陸內(nèi)部深部熱驅(qū)動——地幔底辟有關。目前得到眾多學者普遍認可的主要是， 上地殼為脆性物質(zhì)， 下地殼為柔性物質(zhì)， 且溫度高， 密度低， 導電性好， 地震波速低等， 因此其黏滯系數(shù)較低， 非常柔軟。 在印度板塊擠壓下， 下地殼容易流動， 拖曳上地殼運動， 導致高原的隆升和變形(吳功建等， 1991； 石耀霖等， 1992； 曹建玲等， 2009)。

隨著GPS等空間大地測量技術日漸成熟， 全球有越來越多的GPS觀測臺站投入使用。中國在青藏高原上的GPS觀測臺站觀測到了大量的數(shù)據(jù)， 為研究青藏高原的運動和變形提供了重要的數(shù)據(jù)支持。 如圖2 所示為GPS觀測值， 這些數(shù)據(jù)揭示了青藏高原水平運動的主要特征。 其SN方向上， 主要向N運動， 且從南部到北部， 運動速率逐漸衰減， 顯示出SN向在縮短(Ganetal.， 2007)， 是由地殼增厚形變吸收掉了絕大部分， 另外小部分由高原內(nèi)部的一系列大型走滑斷裂帶吸收了(任金衛(wèi)等， 2003； 沈正康等， 2003； 許志琴等， 2011)。 EW方向上， 主要顯示為EW向的張開， 即東部向E運動， 西部向W運動， 符合物質(zhì)被擠壓后向兩側(cè)移動的力學特征； 但是W向的運動相對微弱， 主要為E向的運動， 且在東南緣， 伴隨有S向的運動， 顯示其為順時針方向運動。這是因為在青藏高原東部有堅硬的四川盆地阻擋， 而其東南部云南到緬甸地區(qū)恰好有1個軟弱的利于巖石流動的通道， 所以物質(zhì)往S流動， 青藏高原整體上形成了順時針旋轉(zhuǎn)的形態(tài)(朱守彪等， 2004； 曹建玲等， 2009)。 在青藏高原西北部為帕米爾高原， 在深部有部分物質(zhì)俯沖至帕米爾高原下， 西側(cè)物質(zhì)向W擠出， 但是擠出量相對微弱(雷建設等， 2002； 張家聲等， 2005； 許志琴等， 2011)。

由于在青藏高原上由GPS觀測數(shù)據(jù)揭示的長期地表變形表現(xiàn)出如此復雜多變的特點， 僅僅用上地殼在邊界作用下的彈性變形解釋是不夠的， 還要考慮柔軟的下地殼流動拖曳上地殼運動的附加作用(曹建玲等， 2009)。那么拖曳力有多大呢？方向如何？是如何分布的？石耀霖等(2000)首先提出了在二維模擬中引入等效體力部分反映某些三維作用的方法。朱守彪等(2005)利用遺傳有限單元算法反演出青藏高原受到的邊界作用力、 地形擴展力以及下地殼對上地殼的拖曳力。曹建玲等(2009)指出下地殼對上地殼的拖曳力依賴于下地殼的黏滯系數(shù)， 但是黏滯系數(shù)的確切給定尚有困難。本文將建立二維平面彈性有限元模型， 利用GPS觀測資料作為約束條件， 選取關鍵點， 利用試錯法給出關鍵點上的力， 同時雙線性插值計算其他位置上的力， 以此為基礎正演模擬得到的長期地表變形與GPS觀測資料吻合較好。

1　模型及邊界條件

1.1區(qū)域模型

本文參考曹建玲等(2009)的研究， 主要計算青藏高原下地殼對上地殼的拖曳作用， 因此選取較簡單模型， 區(qū)域范圍為75°～110°E， 25°～41°N。據(jù)前人研究(Tapponnier， 2001; 張培震等， 2003； 王輝等， 2006； Ganetal.， 2007； 聞學澤等， 2011； 徐錫偉等， 2014)， 青藏高原及周邊區(qū)域活動地塊主要分為拉薩、 羌塘、 巴顏喀拉、 柴達木、 祁連、 川滇、 滇南、 塔里木等， 各個地塊之間主要的斷裂帶有喜馬拉雅俯沖帶、 阿爾金斷裂帶、 祁連-海原斷裂帶、 昆侖山斷裂帶、 甘孜-瑪尼-玉樹斷裂帶、 龍門山斷裂帶、 鮮水河斷裂帶、 嘉黎斷裂帶等。由于一些大的斷裂帶周圍有很多小斷裂帶， 我們將模型簡化， 小斷裂帶整合為大的斷裂帶， 除了喜馬拉雅俯沖帶為20km寬外， 其他斷裂帶均為10km寬。由于本文為二維平面模型， 初步模型采用笛卡爾坐標系， 將研究區(qū)的經(jīng)緯度坐標系投影轉(zhuǎn)換到直角坐標系進行計算。模型及剖分網(wǎng)格如圖1 所示， 采用三角形網(wǎng)格， 在斷裂帶地區(qū)加密， 共有157,649個節(jié)點， 310,439個單元。

圖1　模型及網(wǎng)格Fig. 1　Model and grid.

表1　介質(zhì)參數(shù)表

Table1　Medium parameters

楊氏模量/Pa泊松比塔里木10×10100.25柴達木9×10100.25巴顏喀拉9×10100.25羌塘9×10100.25拉薩9×10100.25雅魯藏布江地塊2×10100.22喜馬拉雅俯沖帶南部(印度板塊)1.5×10100.22阿爾金斷裂帶8×10100.25祁連-海原斷裂帶6×10100.24昆侖山斷裂帶5.5×10100.23甘孜-瑪尼-玉樹斷裂帶5×10100.22龍門山斷裂帶5.0×10100.22龍門山北部斷裂帶5.5×10100.23鮮水河斷裂帶5.0×10100.22理塘斷裂帶5.5×10100.23嘉黎斷裂帶(西)7.5×10100.25嘉黎斷裂帶(東)7.0×10100.24喀喇昆侖斷裂帶(西)6.5×10100.24喀喇昆侖斷裂帶(東)8.0×10100.25喜馬拉雅俯沖帶(西)8×1090.22喜馬拉雅俯沖帶(東)1.0×10100.22

1.2材料參數(shù)

參考前人的研究成果(Huangetal.， 2003； 王輝等， 2006； 黃建平等， 2008； 陳連旺等， 2011)， 并且考慮到各斷裂帶內(nèi)部較軟弱， 給出相對于地塊較低的楊氏模量E。 在最近幾十年來較活躍的， 發(fā)生過地震的斷裂帶相對弱一些， 而較穩(wěn)定的斷裂帶， 稍硬一些。由于印度板塊擠壓歐亞板塊， 青藏高原南緣抬升， 因此藏南地區(qū)軟一些， 包括喜馬拉雅俯沖帶南部的印度板塊也較軟， 而北部有堅硬的塔里木盆地， 所以藏北地區(qū)硬一些。綜合考慮各方面因素， 給出各個地塊和斷裂帶的楊氏模量和泊松比， 如表1 所示。

1.3邊界條件

模型采用青藏高原及周邊的726個觀測臺站上觀測得到的相對于穩(wěn)定的歐亞板塊的GPS水平觀測數(shù)據(jù)(Ganetal.， 2007)(特別指出， 本文給定的GPS速度場及模擬計算結果均為相對于穩(wěn)定的歐亞參考框架)， 如圖2 中的黑線所示； 并將其均勻插值到模型邊界上， 得到有限元模型的邊界條件， 如圖2 中的紅線所示。

圖2　GPS臺站觀測數(shù)據(jù)以及插值邊界條件Fig. 2　GPS observation data and boundary condition via interpolation.

1.4初步結果

在只有GPS邊界條件的情況下， 不考慮下地殼對上地殼的拖曳力作用， 利用以上有限元模型計算研究區(qū)域的長期速度場分布， 將初步計算結果與GPS觀測值之間進行對比， 如圖3 所示， 主要差異為：

(1)在青藏高原東部， 觀測值明顯大于計算值。觀測值顯示地表變形有強烈的SE向旋轉(zhuǎn)， 而實際計算結果只能顯示非常微弱的旋轉(zhuǎn)， 甚至沒有旋轉(zhuǎn)。王輝等(2006)的青藏高原二維數(shù)值模擬結果顯示， 只有觀測的GPS速率邊界條件， 需要在東南邊施加大于觀測值的S向位移， 才能夠跟觀測值吻合。曹建玲等(2009)利用三維流變模型模擬青藏高原GPS位移， 指出在東南緣加載應力邊界條件， 結果與觀測值之間的誤差較小。朱守彪等(2004)對川滇地區(qū)的計算結果表明， 加載拖曳力后， 能夠與GPS觀測值大體吻合。

(2)在西部地區(qū)， 觀測值的相對計算結果有微弱的W向偏轉(zhuǎn)。朱守彪等(2005， 2006)的研究結果指出青藏高原的西部地區(qū)有W向的拖曳力。黃建平等(2008)計算出的地幔對流產(chǎn)生的拖曳力場顯示， 在青藏高原西部地區(qū)存在有W向的拖曳力。

圖3　僅加載邊界條件計算結果與觀測值對比Fig. 3　Comparison between observation and simulation results with only boundary condition.

觀測結果和模擬計算結果對比表明， 僅僅采用GPS邊界條件得到的計算結果不能夠與實際觀測值吻合， 因為沒有考慮柔性的下地殼對上地殼的拖曳力作用。Copley等(2007)模擬了青藏高原地表的GPS速度， 指出主要是下地殼流動， 拖曳上地殼運動， 導致該區(qū)域的構造變形。曹建玲等(2009)的研究結果指出GPS繞喜馬拉雅東構造結的順時針旋轉(zhuǎn)速度大于GPS觀測到的地表速度， 該速度之差每年達到數(shù)毫米； 在柔性下地殼區(qū)域， 下地殼對上地殼有拖曳力作用， 大小一般在0.01～0.1MPa之間。以上結果表明， 在模擬青藏高原及其鄰區(qū)長期變形時， 下地殼對上地殼的拖曳力是不可以忽略的。

2　拖曳力計算

下地殼對上地殼的拖曳力如何計算呢？采用朱守彪等(2004)文章中提出的方法， 在彈性學平衡方程中令水平剪應力分量隨深度變化為常量， 那么在二維問題中， 該部分即可作為等效體力， 代表特定深度內(nèi)下地殼對上地殼的拖曳剪切力作用。如式(1)所示：

(1)

考慮在上1節(jié)的有限元模型中加載拖曳力； 除此之外， 幾何模型、 材料模型和邊界條件均不變。由于模型是二維平面的， 拖曳力在有限元計算中以節(jié)點力或集中力的形式出現(xiàn)。我們采用試錯法調(diào)整計算關鍵點上的拖曳力。關鍵點的選取如圖4b所示， 共26個。

圖4 為區(qū)域劃分簡圖， 4個子區(qū)域組成1個大區(qū)域。如圖4a中所示以1號點為中心點， 共有4個子區(qū)域， 且1號、 2號、 3號和4號， 4個點組成的是其中1個子區(qū)域。每1個子區(qū)域(由4個點組成)作為1個雙線性插值區(qū)， 有4個角點的拖曳力值， 就可以插值求得整個子區(qū)域內(nèi)的其他任意點上的拖曳力。同樣地， 可以得到其他3個子區(qū)域上的值。這樣就得到了整個大區(qū)域上的拖曳力值。參考該示意圖， 我們將模型劃分區(qū)域， 主要觀察哪些地區(qū)的觀測值和模擬值差異大， 以此作為分區(qū)依據(jù)。

圖4　拖曳力區(qū)域劃分圖Fig. 4　Sketch of region division.a 分區(qū)簡單示意圖； b 真實分區(qū)圖， 圓圈內(nèi)的數(shù)字表示子區(qū)域號， 方框內(nèi)的數(shù)字表示關鍵點號

將每個大區(qū)域內(nèi)的每1個子區(qū)域(即4個點組成的區(qū)域)作為1個雙線性插值區(qū)， 區(qū)域內(nèi)各個點的節(jié)點力表示為

(2)

式(2)中的a、b、c、d、m、n、o、p由下面的方程組(3)解出

(3)

式(3)中的x1、 y1、 x2、 y2、 x3、 y3、 x4、 y4分別為4個點的坐標， Fx1、 Fx2、 Fx3、 Fx4和Fy1、 Fy2、 Fy3、 Fy4分別為根據(jù)觀察和參考前人工作給出的試探節(jié)點力值(x為EW向，y為SN向)， 在試驗中不斷改進。然后根據(jù)式(1)即可得到各個區(qū)域內(nèi)任意1個節(jié)點上的力。最終得到的區(qū)域內(nèi)關鍵點的拖曳力值如表2 所示。插值到模型區(qū)域內(nèi)各節(jié)點上的拖曳力如圖5 所示。

表2　各個子區(qū)域內(nèi)關鍵點的拖曳力值(子區(qū)域內(nèi)關鍵點的排列順序如圖4a所示)

Table2　Drag force in different key nodes of each subregion (the order of key nodes in each subregion is shown as in Fig. 4a)

N

圖5　加載至模型中的拖曳力Fig. 5　Drag force in this model.

如圖6 所示， 加載拖曳力后有限元模型計算得到的速度場值與GPS觀測值之間吻合較好， 只有個別地區(qū)仍存在誤差。這已經(jīng)是經(jīng)過多次試驗得到的最好的結果。 按照式(4)計算均方根誤差， 式中的n是所有加載拖曳力的節(jié)點的個數(shù)，Δi是每個點的計算值和觀測值之差，error是均方根誤差。結果如表3 所示， 加載拖曳力之后的均方根誤差明顯小于無拖曳力的情況。

(4)

圖6　邊界條件和拖曳力共同作用的計算結果與臺站觀測對比Fig. 6　Comparison of observation and simulation results with boundary condition and drag force.a 觀測臺站上的計算結果與觀測值對比； b 均勻插值分布的結果對比

表3　計算結果與觀測值均方根誤差對比

Table3　Comparison of MSR error between simulations and observations

EW向分量/mm·a-1SN向分量/mm·a-1矢量模/mm·a-1均勻插值無拖曳力5.335,12.434,45.773,4有拖曳力1.711,72.140,32.685,2觀測臺站無拖曳力4.658,63.097,95.281,9有拖曳力2.194,22.565,43.153,6

圖7　計算的應變率場Fig. 7　Strain rate field by calculation.

表4　主要斷裂帶的走向、 主應變率方向及滑動方向

Table4　Strike and direction of the principle strain rate and slip of major faults

斷裂帶主壓應變率方向主張應變率方向二者比較滑動方向阿爾金斷裂帶NE-SWNW-SE相當左旋海原斷裂帶NE-SWNW-SE壓大于拉左旋昆侖山斷裂帶NE-SWNW-SE相當左旋甘孜玉樹斷裂帶NE-SWNW-SE相當左旋鮮水河斷裂帶NWW-SEENNE-SSW相當左旋龍門山斷裂帶NWW-SEENNE-SSW相當右旋(逆掩不明顯)嘉黎斷裂帶NNE-SSWNWW-SEE壓稍小于拉右旋喀喇昆侖-雅魯藏布江斷裂帶近SN近EW壓遠大于拉右旋喜馬拉雅斷裂帶近SN近EW壓遠大于拉逆掩

計算出的應變率場如圖7 所示， 與前人文中的應變率場基本一致， 具有可比性(沈正康等， 2003； Ganetal.， 2007)。圖7 顯示喜馬拉雅處于SN向擠壓狀態(tài)， 對應逆掩斷層； 藏南地區(qū)EW向拉應變率明顯大于SN向壓應變率， 與該區(qū)SN走向的正斷層吻合； 高原內(nèi)部的數(shù)條走滑斷層， 其右旋或者左旋滑動性質(zhì)， 計算與觀測吻合。主要斷裂帶的詳細情況如表4 所示。

3　結論與討論

模擬計算依賴于材料參數(shù)和邊界條件， 需要參考更多的資料對材料參數(shù)的差異性分布進行約束。本文依靠GPS數(shù)據(jù)插值得到邊界條件， 而研究區(qū)域的GPS臺站數(shù)不足， 尤其在青藏高原中部地區(qū)臺站較少， 插值得到的數(shù)據(jù)精度會受到很大的影響。模型西北部即阿爾金斷裂帶和喀喇昆侖斷裂帶交界的區(qū)域構造非常復雜， 在帕米爾高原還發(fā)生著雙向俯沖(許志琴等， 2011)， 還需要進一步細化研究。

另外， 地球曲率的影響不容忽視， 石耀霖等(2006)指出在利用位移計算應變時， 如果采用平面坐標， 尤其在高緯度地區(qū)， 會產(chǎn)生較大的系統(tǒng)誤差。由于本文為初步模型， 暫且考慮直角坐標系。同時， 雖然本文采用的二維平面模型能夠?qū)栴}簡化， 但是不能考慮真實的三維結構， 如下地殼黏滯系數(shù)， 不能與垂向的GPS實測值進行對比等。因此在深入計算中， 需考慮三維球坐標， 以及更加精細的物質(zhì)結構和邊界條件等。

本文采用加載等效體力的方式來模擬柔性下地殼對脆性上地殼的拖曳剪切力作用。選取關鍵點， 給出試探值， 利用試錯法， 給出最優(yōu)結果， 同時在其他位置點上利用雙線性插值計算其等效體力。最終得到區(qū)域內(nèi)的上地殼下部受到的拖曳力的大小和方向， 以此為基礎模擬得到的速度場結果能夠與GPS觀測有較好的吻合， 誤差較小， 為深入研究青藏高原及其鄰區(qū)地表變形提供了新的思路。

曹建玲， 石耀霖， 張懷， 等. 2009. 青藏高原GPS位移繞喜馬拉雅東構造順時針旋轉(zhuǎn)成因的數(shù)值模擬 [J]. 科學通報， 54(8): 1398—1410.

CAO Jian-ling， SHI Yao-lin， ZHANG Huai，etal. 2009. Numerical simulation of GPS observed clockwise rotation around the eastern Himalayan syntax in the Tibetan plateau [J]. Chinese Science Bulletin， 54(8): 1398—1410(in Chinese).

陳連旺， 詹自敏， 葉際陽， 等. 2011. 流變特性對青藏高原構造變形影響的數(shù)值模擬 [J]. 大地測量與地球動力學， 31(3): 8—14.

CHEN Lian-wang， ZHAN Zi-min， YE Ji-yang，etal. 2011. Numerically modeling the influence of rheological properties on tectonic deformation of Tibet Plateau [J]. Journal of Geodesy and Geodynamics， 31(3): 8—14(in Chinese).

傅容珊， 黃建華， 徐耀民. 2000a. 印度與歐亞板塊碰撞的數(shù)值模擬和現(xiàn)代中國大陸形變 [J]. 地震學報， 22(1): 1—7.

FU Rong-shan， HUANG Jian-hua， XU Yao-min. 2000a. Numerical simulation of the Inida-Asia collision and deformation of China Mainland [J]. Acta Seismologica Sinica， 22(1): 1—7(in Chinese).

傅容珊， 李力剛， 黃建華， 等. 1999. 青藏高原隆升過程的三階段模式 [J]. 地球物理學報， 42(5): 609— 616.

FU Rong-shan， LI Li-gang， HUANG Jian-hua，etal. 1999. Three-step model of the Qinghai-Xizang Plateau uplift [J]. Chinese Journal of Geophysics， 42(5): 609— 616(in Chinese).

傅容珊， 徐耀民， 黃建華， 等. 2000b. 青藏高原擠壓隆升過程的數(shù)值模擬 [J]. 地球物理學報， 43(3): 346—355.

FU Rong-shan， XU Yao-min， HUANG Jian-hua，etal. 2000b. Numerical simulation of the compression uplift of the Qinghai-Xizang plateau [J]. Chinese Journal of Geophysics， 43(3): 346—355(in Chinese).

黃建平， 傅容珊， 鄭勇， 等. 2008. 地幔對流拖曳力對中國大陸巖石層變形的影響 [J]. 地球物理學報， 51(4): 1048—1057.

HUANG Jian-ping， FU Rong-shan， ZHENG Yong，etal. 2008. The influence of mantle convection to the lithosphere deformation of China mainland [J]. Chinese Journal of Geophysics， 51(4): 1048—1057(in Chinese).

雷建設， 周蕙蘭， 趙大鵬. 2002. 帕米爾及鄰區(qū)地殼上地幔P波三維速度結構的研究 [J]. 地球物理學報， 45(6): 802—811.

LEI Jian-she， ZHOU Hui-lan， ZHAO Da-peng. 2002. 3-D velocity structure of P-wave in the crust and upper-mantle beneath Pamir and adjacent region [J]. Chinese Journal of Geophysics， 45(6): 802—811(in Chinese).

任金衛(wèi)， 馬宗晉. 2003. 東亞地區(qū)現(xiàn)代地殼運動特征與構造變形 [J]. 地學前緣， 10(特刊): 58— 65.

REN Jin-wei， MA Zong-jin. 2003. Crustal movement and tectonic deformation of eastern Asia [J]. Earth Science Frontiers， 10: 58— 65(in Chinese).

沈正康， 王敏， 甘衛(wèi)軍， 等. 2003. 中國大陸先進構造應變率場及其動力學成因研究 [J]. 地學前緣， 10(特刊): 93—100.

SHEN Zheng-kang， WANG Min， GAN Wei-jun，etal. 2003. Contemporary tectonic strain rate field of Chinese continent and its geodynamic implications [J]. Earth Science Frontiers， 10: 93—100(in Chinese).

石耀霖， Assumpcao M. 2000. 巴西構造應力場的遺傳有限單元法反演 [J]. 地球物理學報， 2000， 43(2): 166—174.

SHI Yao-lin， Assumpcao M. 2000. Genetic algorithm-finite element inversion of stress field of Brazil [J]. Chinese Journal of Geophysics， 43(2): 166—174(in Chinese).

石耀霖, 朱守彪. 2006. 用GPS位移資料計算應變方法的討論 [J]. 大地測量與地球動力學, 26(1) : 1—8.

SHI Yao-lin, ZHU Shou-biao. 2006. Discussion on method of calculating strain with GPS displacement data [J]. Journal of Geodesy and Geodynamics, 26(1): 1—8(in Chinese).

石耀霖， 朱元清， 沈顯杰. 1992. 青藏高原構造熱演化的主要控制因素 [J]. 地球物理學報， 35(6): 710—720.

SHI Yao-lin， ZHU Yuan-qing， SHEN Xian-jie. 1992. Tectonic processes and thermal evolution of the Qinghai-Xizang(Tibetan)Plateau [J]. Chinese Journal of Geophysics， 35(6): 710—720(in Chinese).

王輝， 張國民， 石耀霖， 等. 2006. 青藏活動地塊區(qū)運動與變形特征的數(shù)值模擬 [J]. 大地測量與地球動力學， 26(2): 15—23.

WANG Hui， ZHANG Guo-min， SHI Yao-lin，etal. 2006. Numerical simulation of movement and deformation of Qinghai-Tibet Plateau [J]. Journal of Geodesy and Geodynamics， 26(2): 15—23(in Chinese).

聞學澤， 杜方， 張培震. 2011. 巴顏喀拉塊體北和東邊界大地震序列的關聯(lián)性與2008年汶川地震 [J]. 地球物理學報， 54(3): 706—716.

WEN Xue-ze， DU Fang， ZHANG Pei-zhen. 2011. Correlation of major earthquake sequences on the northern and eastern boundaries of the Bayan Har block， and its relation to the 2008 Wenchuan earthquake [J]. Chinese Journal of Geophysics， 54(3): 706—716(in Chinese).

吳功建， 高銳， 余欽范. 1991. 青藏高原亞東—格爾木地學斷面綜合地球物理調(diào)查與研究 [J]. 地球物理學報， 34(5): 552—560.

WU Gong-jian， GAO Rui， YU Qin-fan，etal. 1991. Integrated investigations of the Qinghai-Tibet Plateau along the Yadong-Golmud geoscience transect [J]. Chinese Journal of Geophysics， 34(5): 552—562(in Chinese).

徐錫偉， 程佳， 許沖， 等. 2014. 青藏高原塊體運動模型與地震活動主體地區(qū)討論： 魯?shù)楹途肮鹊卣鸬膯⑹?[J]. 地震地質(zhì)， 36(4): 1116—1134.

XU Xi-wei， CHENG Jia， XU Cong，etal. 2014. Discussion on block kinematic model and future themed areas for earthquake occurrence in the Tibetan plateau: Inspiration from the Ludian and Jinggu earthquakes [J]. Seismology and Geology， 36(4): 1116—1134(in Chinese).

許志琴， 姜枚， 楊經(jīng)綏. 1996. 青藏高原北部隆升的深部構造物理作用： 以 “格爾木—唐古拉山”地質(zhì)及地球物理綜合剖面為例 [J]. 地質(zhì)學報， 70(3): 195—206.

XU Zhi-qin， JIANG Mei， YANG Jing-sui. 1996. Tectonophysical process at depth for the uplift of the northern part of the Qinghai-Tibet Plateau: Illustrated by geological and geophysical comprehensive profile from Golmud to the Tanggula Mountains， Qinghai Province [J].Acta Geologica Sinica， 70(3): 195—206(in Chinese).

許志琴， 楊經(jīng)綏， 李海兵， 等. 2011. 印度-亞洲碰撞大地構造 [J]. 地質(zhì)學報， 85(1): 1—33.

XU Zhi-qin， YANG Jing-sui， LI Hai-bing，etal. 2011. On the tectonics of the India-Asia collision [J]. Acta Geologica Sinica， 85(1): 1—33(in Chinese).

張家聲， 單新建， 李建華， 等. 2005. 帕米爾地區(qū)現(xiàn)今大陸深俯沖： 地震構造和動力學解釋 [J]. 巖石學報， 26(4): 1215—1241.

ZHANG Jia-sheng， SHAN Xin-jian， LI Jian-hua，etal. 2005. Recent deep subducting of continental crust in Pamier: An interpretation on seismotectonics and geodynamics [J]. Acta Petrologica Sinica， 26(4): 1215—1241(in Chinese).

張培震， 鄧起東， 張國民， 等. 2003. 中國大陸的強震活動與活動地塊 [J]. 中國科學(D輯)， 33(增刊): 12—20.

ZHANG Pei-zhen， DENG Qi-dong， ZHANG Guo-min，etal. 2003. Strong earthquake activity and active block in China mainland [J]. Science in China(Ser D)， 33: 12—20(in Chinese).

朱守彪， 石耀霖. 2004. 用遺傳有限單元法反演川滇下地殼流動對上地殼的拖曳作用 [J]. 地球物理學報， 47(2): 232—239.

ZHU Shou-biao， SHI Yao-lin. 2004. Genetic algorithm-finite element inversion of drag forces exerted by the lower crust on the upper crust in the Sichuan-Yunnan area [J]. Chinese Journal of Geophysical， 47(2)： 232—239(in Chinese).

朱守彪， 石耀霖. 2005. 青藏高原地形擴展力以及下地殼對上地殼拖曳力的遺傳有限單元法反演 [J]. 北京大學學報(自然科學版)， 41(2): 225—234.

ZHU Shou-biao， SHI Yao-lin. 2005. Genetic algorithm-finite element inversion of topographic spreading forces and drag forces of lower crust to upper crust in Tibetan plateau [J]. Acta Scientiarum Naturalium Universitatis Pekinensis， 41(2): 225—234(in Chinese).

朱守彪， 石耀霖. 2006. 中國大陸及鄰區(qū)構造應力場成因的研究 [J]. 中國科學(D輯)， 36(12): 1077—1083.

ZHU Shou-biao， SHI Yao-lin. 2006. Causes of tectonic stress field in China and adjacent areas [J]. Science in China(Ser D)， 36(12): 1077—1083(in Chinese).

Copley A， Mckenzie D. 2007. Models of crustal flow in the Inida-Asia collision zone [J]. Geophysical Journal International， 169(2): 683— 698.

England P， Mckenzie D. 1982. A thin viscous sheet model for continental deformation [J]. Geophysical Journal International， 70(2): 295—321.

England P， Mckenzie D. 1983. Correction to： A thin viscous sheet model for continental deformation [J]. Geophysical Journal International， 73(2): 523—532.

Gan W J， Zhang P Z， Shen Z K，etal. 2007. Present-day crustal motion within the Tibetan plateau inferred from GPS measurements [J]. Journal of Geophysical Research： Solid Earth(1978—2012)， 112(B8).

Huang Z X， Su W， Peng Y J，etal. 2003. Rayleigh wave tomography of China and adjacent regions [J]. Journal of Geophysical Research: Solid Earth(1978—2012)， 108(B2).

Tapponnier P. 2001. Oblique stepwise rise and growth of the Tibet plateau [J]. Science， 294(5547): 1671—1677.

Tapponnier P， Peltzer G， Dain A Y，etal. 1982. Propagating extrusion tectonics in Asia: New insights from simple experiments with plasticine [J]. Geology， 10(2): 611— 616.

Abstract

The subduction of the Indian plate underneath Eurasian plate results not only in deformation and movement of the elastic upper crust， but also flow of the ductile lower crust in the high temperature and high pressure which drags the brittle upper crust to move at the same time. These two actions work together producing the present movement and deformation field in Tibetan plateau. The dynamics progress has been verified by GPS observation data. Therefore， in a two-dimension plain model， only the elastic deformation with the boundary action at the upper crust cannot explain the deformation well， the drag force acted on the base of upper crust by the drag of ductile flow of the lower crust also need to be considered. However， it’s hard to figure out the magnitude and direction of the drag force. Thus， we established a two-dimension plain elastic finite element model， with the equivalent-body force approach to simulate the drag force. With the internal GPS observation data of Tibetan plateau as constraint condition， we calculated inversely the drag force of key nodes in the model with trial method， and the other nodes in the model with bilinear interpolation method. Finally， we got the drag forces(nodal forces， unit: N)caused by the difference flow of ductile lower crust dragging the brittle upper crust， which are distributed mainly in the region of 86°～100°E and 26°～32°N， the direction is east and south， and the maximum reaches to 1e8N; in some areas in the western part of the study region at 31°～36°N and 76°～80°E， the direction is west， and the maximum reaches to 1e7N. All these work provides a new thought for further research on long-term dynamic mechanism of surface deformation in Tibetan plateau and its surrounding area.

NUMERICAL SIMULATION OF LONG-TERM DEFORMATION OF TIBETAN PLATEAU AND SURROUNDING AREA

DONG Pei-yu1，2)HU Cai-bo2)SHI Yao-lin2)

1)InstituteofSeismology，ChinaEarthquakeAdministration，Wuhan430071，China2)KeyLaboratoryofComputationalGeodynamics，CAS，UniversityofChineseAcademyofSciences，Beijing100049，China

finite element method， numerical simulation， Tibetan plateau， drag forces， bilinear interpolation

2015-01-19收稿， 2015-12-10改回。

中國地震局地震研究所所長基金項目(IS201526237)、 國家科技支撐項目“地震預報實用技術”(2012BAK19B035)、 深部探測數(shù)據(jù)集成與共享服務(201511028)與國家自然科學基金(41474085， 41274027)共同資助。
*

石耀霖， 男， 教授， E-mail: shiyl@ucas.ac.cn。

P315.2

A

0253-4967(2016)02-410-13

董培育， 女， 1987年生， 2015年畢業(yè)于中國科學院大學地球科學學院固體地球物理學專業(yè)， 獲博士學位， 助理研究員， 主要研究方向為地球動力學問題的數(shù)值模擬， E-mail: dongpeiyu97@163.com。

doi:10.3969/j.issn.0253- 4967.2016.02.014