董培育　胡才博　石耀霖*

1)中國地震局地震研究所， 地震大地測量重點(diǎn)實(shí)驗(yàn)室， 武漢　430071 2)中國科學(xué)院計(jì)算地球動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室， 北京　100049

?

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

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

1)中國地震局地震研究所， 地震大地測量重點(diǎn)實(shí)驗(yàn)室， 武漢430071 2)中國科學(xué)院計(jì)算地球動(dòng)力學(xué)重點(diǎn)實(shí)驗(yàn)室， 北京100049

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

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

0　引言

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

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

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

1　模型及邊界條件

1.1區(qū)域模型

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

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

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

Table1　Medium parameters

楊氏模量/Pa泊松比塔里木10×10100.25柴達(dá)木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)部較軟弱， 給出相對(duì)于地塊較低的楊氏模量E。 在最近幾十年來較活躍的， 發(fā)生過地震的斷裂帶相對(duì)弱一些， 而較穩(wěn)定的斷裂帶， 稍硬一些。由于印度板塊擠壓歐亞板塊， 青藏高原南緣抬升， 因此藏南地區(qū)軟一些， 包括喜馬拉雅俯沖帶南部的印度板塊也較軟， 而北部有堅(jiān)硬的塔里木盆地， 所以藏北地區(qū)硬一些。綜合考慮各方面因素， 給出各個(gè)地塊和斷裂帶的楊氏模量和泊松比， 如表1 所示。

1.3邊界條件

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

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

1.4初步結(jié)果

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

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

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

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

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

2　拖曳力計(jì)算

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

(1)

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

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

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

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

(2)

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

(3)

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

表2　各個(gè)子區(qū)域內(nèi)關(guān)鍵點(diǎn)的拖曳力值(子區(qū)域內(nèi)關(guān)鍵點(diǎn)的排列順序如圖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 所示， 加載拖曳力后有限元模型計(jì)算得到的速度場值與GPS觀測值之間吻合較好， 只有個(gè)別地區(qū)仍存在誤差。這已經(jīng)是經(jīng)過多次試驗(yàn)得到的最好的結(jié)果。 按照式(4)計(jì)算均方根誤差， 式中的n是所有加載拖曳力的節(jié)點(diǎn)的個(gè)數(shù)，Δi是每個(gè)點(diǎn)的計(jì)算值和觀測值之差，error是均方根誤差。結(jié)果如表3 所示， 加載拖曳力之后的均方根誤差明顯小于無拖曳力的情況。

(4)

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

表3　計(jì)算結(jié)果與觀測值均方根誤差對(duì)比

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觀測臺(tái)站無拖曳力4.658,63.097,95.281,9有拖曳力2.194,22.565,43.153,6

圖7　計(jì)算的應(yīng)變率場Fig. 7　Strain rate field by calculation.

表4　主要斷裂帶的走向、 主應(yīng)變率方向及滑動(dòng)方向

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

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

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

3　結(jié)論與討論

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

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

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

曹建玲， 石耀霖， 張懷， 等. 2009. 青藏高原GPS位移繞喜馬拉雅東構(gòu)造順時(shí)針旋轉(zhuǎn)成因的數(shù)值模擬 [J]. 科學(xué)通報(bào)， 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. 流變特性對(duì)青藏高原構(gòu)造變形影響的數(shù)值模擬 [J]. 大地測量與地球動(dòng)力學(xué)， 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]. 地震學(xué)報(bào)， 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]. 地球物理學(xué)報(bào)， 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]. 地球物理學(xué)報(bào)， 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. 地幔對(duì)流拖曳力對(duì)中國大陸巖石層變形的影響 [J]. 地球物理學(xué)報(bào)， 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).

雷建設(shè)， 周蕙蘭， 趙大鵬. 2002. 帕米爾及鄰區(qū)地殼上地幔P波三維速度結(jié)構(gòu)的研究 [J]. 地球物理學(xué)報(bào)， 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)代地殼運(yùn)動(dòng)特征與構(gòu)造變形 [J]. 地學(xué)前緣， 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ìn)構(gòu)造應(yīng)變率場及其動(dòng)力學(xué)成因研究 [J]. 地學(xué)前緣， 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. 巴西構(gòu)造應(yīng)力場的遺傳有限單元法反演 [J]. 地球物理學(xué)報(bào)， 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ì)算應(yīng)變方法的討論 [J]. 大地測量與地球動(dòng)力學(xué), 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. 青藏高原構(gòu)造熱演化的主要控制因素 [J]. 地球物理學(xué)報(bào)， 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. 青藏活動(dòng)地塊區(qū)運(yùn)動(dòng)與變形特征的數(shù)值模擬 [J]. 大地測量與地球動(dòng)力學(xué)， 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).

聞學(xué)澤， 杜方， 張培震. 2011. 巴顏喀拉塊體北和東邊界大地震序列的關(guān)聯(lián)性與2008年汶川地震 [J]. 地球物理學(xué)報(bào)， 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. 青藏高原亞東—格爾木地學(xué)斷面綜合地球物理調(diào)查與研究 [J]. 地球物理學(xué)報(bào)， 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. 青藏高原塊體運(yùn)動(dòng)模型與地震活動(dòng)主體地區(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. 青藏高原北部隆升的深部構(gòu)造物理作用： 以 “格爾木—唐古拉山”地質(zhì)及地球物理綜合剖面為例 [J]. 地質(zhì)學(xué)報(bào)， 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. 印度-亞洲碰撞大地構(gòu)造 [J]. 地質(zhì)學(xué)報(bào)， 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)今大陸深俯沖： 地震構(gòu)造和動(dòng)力學(xué)解釋 [J]. 巖石學(xué)報(bào)， 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. 中國大陸的強(qiáng)震活動(dòng)與活動(dòng)地塊 [J]. 中國科學(xué)(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. 用遺傳有限單元法反演川滇下地殼流動(dòng)對(duì)上地殼的拖曳作用 [J]. 地球物理學(xué)報(bào)， 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. 青藏高原地形擴(kuò)展力以及下地殼對(duì)上地殼拖曳力的遺傳有限單元法反演 [J]. 北京大學(xué)學(xué)報(bào)(自然科學(xué)版)， 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ū)構(gòu)造應(yīng)力場成因的研究 [J]. 中國科學(xué)(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改回。

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

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

P315.2

A

0253-4967(2016)02-410-13

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

doi:10.3969/j.issn.0253- 4967.2016.02.014