青藏高原及周邊2001年以來三次特大地震引起的形變場(chǎng)分布特征

張?jiān)?鄭曉靜， 王蘭民

1 蘭州大學(xué)土木工程與力學(xué)學(xué)院， 蘭州　730000 2 中國(guó)地震局蘭州地震研究所， 蘭州　730000

?

青藏高原及周邊2001年以來三次特大地震引起的形變場(chǎng)分布特征

張?jiān)?,2， 鄭曉靜1*， 王蘭民2

1 蘭州大學(xué)土木工程與力學(xué)學(xué)院， 蘭州730000 2 中國(guó)地震局蘭州地震研究所， 蘭州730000

基于地殼黏彈模型在GPS觀測(cè)資料和地震位錯(cuò)數(shù)據(jù)為約束條件下，應(yīng)用有限元數(shù)值模擬方法對(duì)2001年昆侖山8.1級(jí)、2008年汶川8.0級(jí)和2015年尼泊爾8.1級(jí)特大地震引起的地殼形變分布特征進(jìn)行了模擬計(jì)算，獲得了地震位移場(chǎng)和形變場(chǎng).這3次大地震分別發(fā)生在青藏高原的北部、東部邊界和南部邊界，盡管震級(jí)大小基本相當(dāng)，但發(fā)震區(qū)域和發(fā)震斷層性質(zhì)都各不相同，其各自產(chǎn)生的地殼形變場(chǎng)分布特征存在明顯的差異，主要表現(xiàn)為形變場(chǎng)區(qū)域大小、幅值大小等的差異，以及在不同地殼深度也存在明顯的差異，這些差異主要取決于地震斷層性質(zhì)和地下介質(zhì)結(jié)構(gòu)的不同.通過模擬計(jì)算，可以進(jìn)一步了解大地震產(chǎn)生的應(yīng)力加卸載區(qū)分布特征，對(duì)預(yù)測(cè)未來地震發(fā)生區(qū)域范圍提供重要參考依據(jù).

地震； 應(yīng)力加卸載區(qū)； 形變場(chǎng)特征

1　引言

由于青藏高原的地形地貌、深部結(jié)構(gòu)和地質(zhì)構(gòu)造極具特殊性，其動(dòng)力學(xué)和運(yùn)動(dòng)學(xué)過程相當(dāng)復(fù)雜，它的演化過程廣泛引起地學(xué)家的關(guān)注.自有限元方法發(fā)展到地學(xué)領(lǐng)域以來，在青藏地區(qū)開展了不少有限元分析和研究，獲得了一些較有價(jià)值的研究成果(Tapponnier and Molnar,1976; 傅容珊等,2000;陸師闊和蔡永恩,2004;鄧宗策等,1990;戴黎明等,2011).有限元方法已成為研究和探索該特殊地區(qū)的動(dòng)力學(xué)過程的重要方法之一.

在應(yīng)用有限元法對(duì)青藏高原進(jìn)行的構(gòu)造應(yīng)力場(chǎng)模擬方面，地學(xué)家們認(rèn)為印度板塊以一定的速度持續(xù)向北擠入，這是青藏高原隆升的主要?jiǎng)恿?傅容珊等的研究認(rèn)為擠壓過程主要是受巖石的力學(xué)特性、邊界條件以及剝蝕作用等因素的影響，并認(rèn)為高原隆升過程是一個(gè)非均勻的演化過程(傅容珊，2000；鄭洪偉等,2006).陳開平等的研究認(rèn)為印度大陸和歐亞大陸碰撞是川滇菱形塊體向東南運(yùn)動(dòng)的動(dòng)力源(陳開平和馬瑾,1995).楊立強(qiáng)等(2006)對(duì)青藏高原殼幔形變模擬分析認(rèn)為，青藏高原3D殼幔形變主要特征表現(xiàn)為主碰撞帶在緯向上存在距離上的差異，在經(jīng)向上地殼物質(zhì)存在“逃逸”現(xiàn)象.張東寧等(1995)用黏彈性冪律蠕變本構(gòu)關(guān)系探究了高原內(nèi)部張性構(gòu)造應(yīng)力分布狀態(tài)的可能機(jī)制.在此基礎(chǔ)上，以現(xiàn)代構(gòu)造活動(dòng)為主，建立了青藏地區(qū)構(gòu)造形變場(chǎng)有限元模型(張東寧等,2007).王輝等(2006)利用二維有限元模型，以GPS資料為約束條件，探討了青藏地區(qū)活動(dòng)塊體在應(yīng)變積累過程中的變形與運(yùn)動(dòng)狀態(tài).

在應(yīng)用有限元法對(duì)青藏高原進(jìn)行的地震機(jī)制與地震預(yù)測(cè)模擬方面，曹雪峰等(2003)研究顯示南北地震帶上的地震存在呼應(yīng)特征，認(rèn)為這可能是區(qū)域應(yīng)力場(chǎng)平衡的結(jié)果，是地塊邊界力(外載荷)作用的整體表現(xiàn).陳祖安等(2008)模擬分析了1997年瑪尼7.5級(jí)地震對(duì)青藏地區(qū)地塊系統(tǒng)穩(wěn)定性的影響.獲得了該地震引起的地塊邊界帶上庫侖應(yīng)力的變化，其結(jié)果表明除了發(fā)震斷層兩端庫侖破裂應(yīng)力增大和更加集中外，2001年昆侖山8.1級(jí)地震發(fā)震斷層段(東昆侖斷裂中段)上的庫侖破裂應(yīng)力增加了約2 MPa，2008年改則6.9級(jí)地震發(fā)震斷層段(喀拉昆侖斷裂帶)上增加了約0.7 MPa，其應(yīng)力的加載使得這2次地震的發(fā)震斷層(接近破裂強(qiáng)度)處于失穩(wěn)狀態(tài)，加速了地震的發(fā)生.陳連旺等(2008)研究認(rèn)為，在一次強(qiáng)震發(fā)生之后，發(fā)震斷層本身是處于快速卸載的過程，而在其他潛在地震孕育活動(dòng)斷裂帶上庫侖破裂應(yīng)力處于加載的過程，可能加載作用起到主要作用，其相互作用主要表現(xiàn)為應(yīng)力的加載，一次強(qiáng)震的發(fā)生可能對(duì)下一個(gè)強(qiáng)震的孕育過程起到加速的作用，引起強(qiáng)震或多次地震的發(fā)生，直到區(qū)域應(yīng)變能達(dá)到低值水平狀態(tài)，區(qū)域地震將進(jìn)入一個(gè)新的平靜時(shí)段.

以上兩方面的大尺度研究，其有限元模型是平均化意義上的三維模型，一般把地塊當(dāng)成均勻?qū)嶓w來處理，也就是說同一地塊是用一組力學(xué)參數(shù)來描述的，但實(shí)際上地塊內(nèi)部間的非均勻性是明顯存在的，使用更為精細(xì)的模型進(jìn)行模擬計(jì)算，其結(jié)果無疑更加接近實(shí)際.目前能夠反映整個(gè)模型區(qū)域三維結(jié)構(gòu)的模型是全球地殼CRUST 1.0模型，該模型主要把地殼分為上、中、下三層，水平分辨尺度為1°×1°.CRUST 1.0全球地殼模型是本研究模型建立的重要數(shù)據(jù)源.本文以2001年11月14日昆侖山口8.1級(jí)、2008年5月12日汶川8.0級(jí)和2015年4月25日尼泊爾8.1級(jí)3次特大地震為例，應(yīng)用有限元數(shù)值模擬方法，研究震級(jí)相近的大地震在不同地區(qū)和不同構(gòu)造條件下產(chǎn)生的物理場(chǎng)特征及其差異性.

2　建立有限元黏彈模型

2.1模型區(qū)域與實(shí)體模型的建立

模型區(qū)域包括青藏地塊區(qū)的6個(gè)地塊、塔里木地塊南部、阿拉善地塊東南部、鄂爾多斯地塊西部、華南地塊西部、滇緬地塊區(qū)的滇南地塊和滇西地塊，共包括12個(gè)地塊和15條邊界帶(張培震等,2003; 張國(guó)民等,2005)，其中Ⅰ級(jí)地塊區(qū)邊界有8條，次級(jí)地塊邊界有7條.區(qū)域東西長(zhǎng)約2800 km,南北寬約1800 km.根據(jù)全球地殼模型CRUST 1.0數(shù)據(jù)，可以獲得模型區(qū)上、中、下地殼1°×1°網(wǎng)格的P波速度、S波速度、介質(zhì)密度和幾何參數(shù)等數(shù)據(jù)，利用這些參數(shù)可計(jì)算網(wǎng)格中的楊氏模量和泊松比.上地殼為完全彈性模型.中、下地殼為黏彈模型，其黏滯系數(shù)分別設(shè)置為1019和1020Pa·s，各層黏性是均勻的.地殼上、中、下三層模型由1521個(gè)實(shí)體拼接而成，完成了實(shí)體模型的建立.

2.2模型網(wǎng)格劃分

ANSYS系統(tǒng)中可提供使用的單元類型有上百種，我們所分析的模型是不規(guī)則的，希望或要求模型在施加外力后，其相應(yīng)的變形量小.為此可以優(yōu)先選用solid187單元，其單元為具有6個(gè)中間節(jié)點(diǎn)的高階三維固體結(jié)構(gòu)的單元，具有二階位移插值精度，可以更好地表示不規(guī)則模型，支持大變形和大應(yīng)變等模型.為了描述剪切松弛核函數(shù)G(t)和體積松弛核函數(shù)K(t)，我們采用的是Prony級(jí)數(shù)形式.在網(wǎng)格劃分時(shí)，采用網(wǎng)格單元平均長(zhǎng)度為20 km的四面體單元進(jìn)行模型非結(jié)構(gòu)化網(wǎng)格自動(dòng)劃分，共劃分為488707個(gè)單元，節(jié)點(diǎn)數(shù)為740841個(gè).

2.3模型斷層參數(shù)設(shè)定

斷裂帶的描述有2種形式，一種是用曲線線條數(shù)表示，即一條斷裂是用多條線段連接而成的線條；另一種是用一系列的四邊形區(qū)域組成的一條不等寬度區(qū)帶來表示，這是一種接近斷裂帶實(shí)際寬度的描述.第一種形式的斷裂可根據(jù)張國(guó)民(2005)的對(duì)地塊邊界帶的劃分寬度進(jìn)行給定.對(duì)位于這些斷裂帶或地塊邊界帶內(nèi)所有單元的材料參數(shù)進(jìn)行弱化(陳連旺等,2011).材料參數(shù)弱化規(guī)則為：上地殼彈性層的斷層參數(shù)：楊氏模量=原楊氏模量/3,泊松比=原泊松比+0.02；中、下地殼黏彈層的斷層參數(shù)：楊氏模量=原楊氏模量/3,泊松比=原泊松比+0.02，黏性系數(shù)=原黏性系數(shù)/2.

2.4模型邊界條件

根據(jù)甘衛(wèi)軍(私人交流)提供的GPS數(shù)據(jù)，對(duì)模型區(qū)進(jìn)行數(shù)據(jù)挑選，對(duì)稀疏的區(qū)域進(jìn)行必要的數(shù)據(jù)插值，在喜馬拉雅以南地區(qū)GPS測(cè)點(diǎn)非常少，僅有幾個(gè)點(diǎn)，應(yīng)用這幾個(gè)測(cè)點(diǎn)的GPS數(shù)據(jù)，在平面上增補(bǔ)了幾個(gè)控制點(diǎn).在此基礎(chǔ)上，對(duì)GPS數(shù)據(jù)進(jìn)行二維平面插值，獲得要求密度點(diǎn)位的年位移，可對(duì)模型側(cè)邊界進(jìn)行位移邊界約束，在側(cè)邊界豎直方向位移條件是等同的.模型表層為自由面，約定下地殼底界面垂直于地表方向不產(chǎn)生位移，其分向位移為零.在模型構(gòu)造結(jié)處進(jìn)行水平位移約束，即水平位移分量為零.在模型的南邊界向北推擠的最大年位移量可接近50 mm，模型側(cè)邊界不同區(qū)面位移邊界約束大小和方向都存在一定差異.

2.5地震斷層位錯(cuò)加載

根據(jù)昆侖山口8.1級(jí)(震中90.90°E,36.40°N)、汶川8.0級(jí)(震中103.40°E,30.95°N)和尼泊爾8.1級(jí)(震中84.70°E,28.20°N)3次大地震(圖1a)的實(shí)際野外考察和同震資料，昆侖山口8.1級(jí)地震斷層為左旋走滑斷層，傾角陡(80°)，極震區(qū)平均水平位移滑動(dòng)量為6 m，斷裂長(zhǎng)度為350 km(陳文彬等,2001)；汶川8.0級(jí)地震斷層為右旋逆走滑(徐錫偉等,2008)，傾角較陡(60°)，全長(zhǎng)260 km，南端三江—映秀段，水平位移為0～3 m，垂直位移0～3.5 m.映秀—平驛段，平均水平位移為3～3.5 m，平均垂直位移3.5～3.7 m.平驛—青川，水平位移為3.5～0 m，垂直位移3.7～0 m.尼泊爾8.1級(jí)地震斷層為逆斷層，傾角為11°，全長(zhǎng)約120 km，斷層上盤向南的最大位移量為3.5 m(張貝等,2015).由于這些地震都發(fā)生在地塊的邊界帶上，邊界帶的寬度不小于60 km，滿足斷裂在深部也位于相應(yīng)的邊界帶內(nèi)的情況.地震位錯(cuò)斷層加載寬度60 km(斷層兩邊各30 km)，深度為25 km(假定震源深度為15 km)，對(duì)于走滑位錯(cuò)量，兩盤各加載一半，方向相反.對(duì)于逆沖或正斷位錯(cuò)量，只對(duì)上盤進(jìn)行加載.震源深度至上地殼頂部為均勻加載，其他方向?yàn)榫€性衰減加載(最小加載位錯(cuò)量為0.1 m).

3　計(jì)算結(jié)果分析

由于模式為黏彈性模型，在計(jì)算時(shí)采用的是時(shí)間步計(jì)算，參考時(shí)間為1900年或其他時(shí)間，參考時(shí)間之前以1萬年為起算年，從參考時(shí)間起按黏彈性模型進(jìn)行計(jì)算，每個(gè)地震需要計(jì)算兩步，地震發(fā)生之前為第一步，地震位錯(cuò)加載后為第二步，按地震發(fā)生時(shí)間依次進(jìn)行時(shí)間步計(jì)算，并對(duì)目標(biāo)地震位錯(cuò)加載計(jì)算結(jié)果與其地震加載之前的計(jì)算結(jié)果進(jìn)行差值處理，可以獲得目標(biāo)地震產(chǎn)生的位移場(chǎng)和應(yīng)力應(yīng)變場(chǎng)，下面重點(diǎn)對(duì)上地殼中的各種物理場(chǎng)進(jìn)行分析.

3.1位移分布特征

2001年11月14日昆侖山口8.1級(jí)地震位于巴顏喀拉地塊與柴達(dá)木地塊的邊界東昆侖帶的中段，野外考察極震區(qū)平均水平位移滑動(dòng)量為6 m，實(shí)際計(jì)算上地殼頂界面破裂面極震區(qū)的平均位移為5.9 m，水平位移量大于0.4 m的分布區(qū)域位于地震破裂帶的南北兩側(cè)，存在一定的斜對(duì)稱性，在地塊內(nèi)部存在明顯的非均勻性.1 cm以上位移量的區(qū)域基本上遍及青藏高原中部和北部地區(qū)，以及塔里木地塊的東南邊緣(圖1b)，地塊邊界為強(qiáng)衰減帶.位移方向基本以發(fā)震斷層為界，北部向西，南部向東，與同震GPS觀測(cè)一致(萬永革等,2004)，但在較遠(yuǎn)區(qū)域位移方向存在一定的變化.2008年5月12日汶川8.0級(jí)地震位于巴顏喀拉地塊與華南地塊的邊界岷江—龍門山帶上，野外考察極震區(qū)平均水平位移為3～3.5 m，平均垂直位移為3.5～3.7 m.實(shí)際計(jì)算上地殼頂界面破裂面極震區(qū)的平均位移量為3.3 m，水平位移量大于0.3 m的分布區(qū)域位于地震斷層兩側(cè)，有一定的對(duì)稱性.1 cm以上位移量的區(qū)域也位于斷層帶兩側(cè)(圖1c)，西部區(qū)域范圍大，巴顏喀拉地塊北邊界為明顯的衰減帶.東部區(qū)域同等位移量分布范圍小.位移方向?yàn)槲鞑繀^(qū)域向東、北部區(qū)域?yàn)闁|南、東部區(qū)域方向?yàn)镹WW，與GPS同震位移觀測(cè)結(jié)果基本一致(徐韶光等,2010)，南部區(qū)域以安寧河—小江帶為界位移方向反向.與昆侖山口8.1級(jí)地震類似，地塊邊界為強(qiáng)衰減帶，位移大范圍分布呈現(xiàn)出非均勻性和不對(duì)稱性.2015年4月25日尼泊爾8.1級(jí)地震位于喜馬拉雅帶的中段，發(fā)震斷層為逆沖斷層，最大位錯(cuò)量為3.5 m，實(shí)際計(jì)算接近3.6 m，位移量大于0.3 m的分布區(qū)域主要位于地震破裂帶以北的NNE區(qū)域的拉薩地塊中部地域，呈口字型分布，每經(jīng)過一條邊界帶，其位移衰減是明顯的.位移方向表現(xiàn)為極震區(qū)以南為NNE，北部地區(qū)向南或指向震中(圖1d).這次地震主要造成青藏高原向南運(yùn)動(dòng).

3.2應(yīng)變場(chǎng)分布特征

昆侖山口8.1級(jí)地震引起的加卸載區(qū)各2個(gè)，卸載區(qū)域位于破裂帶的南北兩側(cè)，等效應(yīng)變和最大剪應(yīng)變?yōu)樨?fù)值，等效應(yīng)變減少0.01×10-5～0.12×10-5，最大剪應(yīng)變減少0.01×10-5～0.07×10-5.加載區(qū)域位于破裂帶兩端，等效應(yīng)變和最大剪應(yīng)變?yōu)檎?，等效?yīng)變?cè)黾?.01×10-5～0.12×10-5，最大剪應(yīng)變?cè)黾?.01×10-5～0.07×10-5(圖2a—2b).這次地震之后的5年時(shí)間里，青藏高原70%的5級(jí)以上地震和約90%的6級(jí)以上地震發(fā)生在加載區(qū)或加卸載區(qū)交界部位.汶川8.0級(jí)地震引起的加卸載區(qū)各有2個(gè)，卸載區(qū)分別位于四川盆地、巴顏喀拉地塊東部地區(qū)和川滇地塊最北部地區(qū)，等效應(yīng)變減少0.01×10-5～0.07×10-5，最大剪應(yīng)變減少0.01×10-5～0.045×10-5.加載區(qū)分別是地震破裂帶的兩端、川滇地塊中部地區(qū)、以及北部和西北部地區(qū)，即岷江—龍門山帶中北段(編號(hào)10)、安寧河—小江帶北段(編號(hào)11)，等效應(yīng)變?cè)黾?.01×10-5～0.07×10-5，最大剪應(yīng)變?cè)黾?.01×10-5～0.045×10-5(圖2c—2d).2010年的玉樹7.1級(jí)地震、2013年的蘆山7.0級(jí)地震和岷縣漳縣6.6級(jí)地震以及強(qiáng)余震都發(fā)生在加載區(qū).尼泊爾8.1級(jí)地震引起的加卸載區(qū)各1個(gè)，極震區(qū)北部地區(qū)大范圍為卸載區(qū)，等效應(yīng)變減少0.01×10-5～0.12×10-5，最大剪應(yīng)變減少0.01×10-5～0.07×10-5.加載區(qū)主要分布于喜馬拉雅帶內(nèi)，等效應(yīng)變?cè)黾?.01×10-5～0.09×10-5，最大剪應(yīng)變?cè)黾?.01×10-5～0.045×10-5(圖2e—2f)，7級(jí)以上強(qiáng)余震都發(fā)生在加載區(qū).

圖1　三次特大地震位置及其在上地殼產(chǎn)生的同震水平位移分布圖(位移單位：m)(a) 震中分布，紅色實(shí)心圓標(biāo)出了3個(gè)地震的位置；(b)、(c)和(d)分別是2001年昆侖山口8.1級(jí)、2008年汶川8.0級(jí)和2015年尼泊爾8.1級(jí)的同震位移； 1.喜馬拉雅帶,2.喀喇昆侖—嘉黎帶,3.瑪尼—玉樹帶,4.東昆侖帶,5.西秦嶺—德令哈帶,6.西昆侖帶,7.阿爾金帶,8.海原—祁連帶,9.賀蘭帶,10.岷江—龍門山帶,11.安寧河—小江帶,12.鮮水河帶,13.三江帶,14.瀾滄江帶,15.紅河帶.Fig.1　Three earthquakes location and co-seismic horizontal displacement patterns caused by three earthquakes in the upper crust (displacement unit: m)(a) The epicenter distribution； (b), (c) and (d) co-seismic displacement of 2001 Kunlun Mountain MS8.1, 2008 Wenchuan MS8.0 and 2015 Nepal MS8.1; 1. Himalayan zone, 2. Karakorum-Jiali belt, 3. Mani-Yushu belt, 4. East Kunlun Mountains, 5. West Qinling-Delhi belt, 6. West Kunlun belt, 7. Altyn belt, 8. Haiyuan-Qilian belt, 9. Helan belt, 10. Minjiang-Longmen Shan belt, 11. Anning River-Little River belt, 12. Xianshuihe belt, 13. Sanjiang belt, 14. Lancang River belt, 15. Red River belt.

圖2　三次特大地震在上地殼產(chǎn)生的同震等效應(yīng)變和最大剪應(yīng)變分布圖(×10-5)(a)—(b)分別是2001年昆侖山口8.1級(jí)地震的同震等效應(yīng)變和最大剪應(yīng)變； (c)—(d)分別是2008年汶川8.0級(jí)地震的同震等效應(yīng)變和最大剪應(yīng)變； (e)—(f)分別是2015年尼泊爾8.1級(jí)地震的同震等效應(yīng)變和最大剪應(yīng)變.Fig.2　Distribution of co-seismic equivalent strain and maximum shear strain caused by three earthquakes in the upper crust (×10-5)(a)—(b) Co-seismic equivalent strain and the maximum shear strain of 2001 Kunlun Mountain MS8.1 earthquake, respectively； (c)—(d) Co-seismic equivalent strain and the maximum shear strain of 2008 Wenchuan MS8.0 earthquake, respectively； (e)—(f) Co-seismic equivalent strain and the maximum shear strain of 2015, Nepal MS8.1 earthquake, respectively.

3.3地震破裂帶應(yīng)力應(yīng)變場(chǎng)分布特征

一次大地震的發(fā)生，其發(fā)震斷層兩盤會(huì)產(chǎn)生相對(duì)運(yùn)動(dòng)，運(yùn)動(dòng)的位移稱為地震位錯(cuò).我們可以通過野外考察的地震位錯(cuò)數(shù)據(jù)來模擬計(jì)算深部地震時(shí)的應(yīng)力場(chǎng)狀態(tài)，以及地震產(chǎn)生的應(yīng)變場(chǎng)分布，幫助我們了解需要多大的應(yīng)力才能產(chǎn)生實(shí)際大小的位錯(cuò).也就是說發(fā)生相等大小的地震需要釋放多大的應(yīng)力和應(yīng)變，用地震時(shí)的應(yīng)力和應(yīng)變分別與地震前的應(yīng)力和應(yīng)變相減，在斷層處的值為正值；如果用地震前的值減去地震后的值則結(jié)果為負(fù)值，這個(gè)量值我們可以理解為是地震釋放的相關(guān)量值.這3次大地震所產(chǎn)生的等效應(yīng)力和等效應(yīng)變、最大剪應(yīng)力和最大剪應(yīng)變?nèi)鐖D3和圖4所示，3次地震的震源體大小存在區(qū)別，尤其在不同深度上，應(yīng)力和應(yīng)變分布差異明顯.汶川地震的應(yīng)力應(yīng)變分布狀態(tài)在不同深度上差異不大，這一點(diǎn)可以說明震級(jí)相當(dāng)?shù)倪@3次大地震在水平方向產(chǎn)生位移場(chǎng)、應(yīng)力和應(yīng)變場(chǎng)存在較大的變化，汶川8.0級(jí)地震所產(chǎn)生的應(yīng)力應(yīng)變場(chǎng)分布于整個(gè)地殼，而昆侖山口大地震和尼泊爾大地震所產(chǎn)生的應(yīng)力主要分布于中、上地殼，應(yīng)變場(chǎng)在下地殼明顯減弱.這可能與地震斷層性質(zhì)和所處的空間介質(zhì)結(jié)構(gòu)有關(guān).一次特大地震發(fā)生，釋放的應(yīng)力可達(dá)10 MPa以上，產(chǎn)生的最大剪應(yīng)力為6 MPa以上，產(chǎn)生的等效應(yīng)變接近3.0×10-4，最大剪應(yīng)變可大于1.0×10-4(圖3和圖4).

圖3　三次大地震破裂帶不同深度層的等效應(yīng)力和等效應(yīng)變分布圖Fig.3　Distribution of equivalent stress and equivalent strain at different depths of three major earthquake rupture belts

圖4　三次大地震破裂帶不同深度層的最大剪應(yīng)力和最大剪應(yīng)變分布圖Fig.4　Distribution of the maximum shear stress and maximum shear strain at different depths of three major earthquake rupture belts

4　結(jié)論與討論

我們應(yīng)用有限元方法對(duì)青藏高原及周緣的三次特大地震進(jìn)行數(shù)值模擬和分析，可得到如下的結(jié)論.

(1) 這3次特大地震產(chǎn)生的位移場(chǎng)、應(yīng)力應(yīng)變場(chǎng)在空間分布上存在異同，昆侖山地震斷層為純走滑型，其位移場(chǎng)、應(yīng)力應(yīng)變場(chǎng)呈4象限分布，與左旋走滑斷層模型吻合；尼泊爾地震斷層為逆沖型，其位移場(chǎng)、應(yīng)力應(yīng)變場(chǎng)呈2象限分布，與斷層模型吻合；汶川地震斷層為逆走滑型，其位移場(chǎng)、應(yīng)力應(yīng)變場(chǎng)呈4象限分布，但分布特征較昆侖山和尼泊爾地震復(fù)雜.

(2) 一次大地震的發(fā)生，會(huì)對(duì)周圍應(yīng)力應(yīng)變場(chǎng)環(huán)境引起改變，形成明顯的加卸載區(qū)(正負(fù)值相間的區(qū)域)，應(yīng)力應(yīng)變加載區(qū)是強(qiáng)余震和未來幾年內(nèi)的強(qiáng)烈地震主要活動(dòng)地區(qū).

(3) 一次大地震的發(fā)生，在極震區(qū)可產(chǎn)生10 MPa以上的應(yīng)力降，應(yīng)變的改變量可達(dá)2.5×10-4以上.不同地震斷層性質(zhì)和不同深部結(jié)構(gòu)對(duì)地震產(chǎn)生的位移場(chǎng)、應(yīng)力應(yīng)變場(chǎng)在空間分布上存在明顯的差異.

在應(yīng)用有限元方法對(duì)地學(xué)模型進(jìn)行數(shù)值模擬時(shí)，我們應(yīng)該看到，計(jì)算結(jié)果依賴于地球內(nèi)部結(jié)構(gòu)、邊界約束條件和地震位錯(cuò)等數(shù)據(jù)的合理性和精度.隨著地球內(nèi)部的精細(xì)探測(cè)以及高密度和高精度地殼運(yùn)動(dòng)觀測(cè)，可以獲得精度更高的地殼結(jié)構(gòu)模型和地殼運(yùn)動(dòng)的時(shí)間位移模型，同時(shí)可以應(yīng)用現(xiàn)代觀測(cè)手段和測(cè)量技術(shù)得到更高精度的地震斷層位錯(cuò)數(shù)據(jù).以此為基礎(chǔ)，可對(duì)每次大地震進(jìn)行模擬計(jì)算，獲得應(yīng)力應(yīng)變加卸載分布區(qū)域，為未來地震發(fā)生研究提供區(qū)域依據(jù).

致謝中國(guó)地震局地質(zhì)研究所甘衛(wèi)軍研究員提供了寶貴的GPS資料；在有限元軟件的使用過程中得到了中國(guó)地震局蘭州地震研究所楊興悅副研究員的大力幫助；匿名審稿專家提出的十分寶貴的修改意見和建議，在此一并表示感謝！

Cao X F, Yin J Y, Yang L M. 2003. Study of simulation for the evolution of stress field on north-south seismic belt in China.NorthwesternSeismologicalJournal(in Chinese), 25(3): 221-225.

Chen K P, Ma J. 1995. Numerical analysis of tectonic deformation by continental collision between India and Eurasia.SeismologyandGeology(in Chinese), 17(3): 277-284.

Chen L W, Zhang P Z, Lu Y Z, et al. 2008. Numerical simulation of loading/unloading effect on Coulomb failure stress among strong earthquakes in Sichuan-Yunnan area.ChineseJ.Geophys. (in Chinese), 51(5): 1411-1421.

Chen L W, Zhan Z M, Ye J Y, et al. 2011. Numerically modeling the influence of rheological properties on tectonic deformation of Tibet plateau.JournalofGeodesyandGeodynamics(in Chinese), 31(3): 8-14.

Chen W B, Xu X W, Zhang Z J, et al. 2001. Preparatory field investigation on the surface deformation zone of the Qinghai-XinjiangMS8.1 earthquake on Nov. 14, 2001.NorthwesternSeismologicalJournal(in Chinese), 23(4): 313-317.

Chen Z A, Lin B H, Bai W M. 2008. 3-D numerical simulation on influence of 1997 Mani earthquake occurrence to stability of tectonic blocks system in Qingzang and Chuandian zone using DDA+FEM method.ChineseJ.Geophys. (in Chinese), 51(5): 1422-1430.

Dai L M, Li S Z, Tao C H, et al. 2011. Three-dimensional numberical modeling of activity of the Longmenshan fault zone driven by the India Plate.ProgressinGeophys. (in Chinese), 26(1): 41-51, doi: 10.3969/j.issn.1004-2903.2011.01.004.

Deng Z C, Cui J W, Feng X F, et al. 1990. The analysis of tectonic dynamics and deformation mechanics of Qinghai-Xizang plateau by elastic finite element method.BulletinoftheChineseAcademyofGeologicalSciences(in Chinese), 21: 65-67.

Fu R S, Xu Y M, Huang J H, et al. 2000. Numerical simulation of the compression uplift of the Qinghai-Xizang plateau.ChineseJ.Geophys. (in Chinese), 43(3): 346-355.

Lu S K, Cai Y E. 2004. Review of numerical simulation on the dynamics of Qinghai-Xizang plateau.ActaSeismologicaSinica(in Chinese), 26(5): 547-559.

Tapponnier P, Molnar P. 1976. Slip-line field theory and large-scale continental tectonics.Nature, 264(5584): 319-324.

Wan Y G, Wang M, Shen Z K, et al. 2004. Co-seismic slip distribution of the 2001 west of Kunlun mountain pass earthquake inverted by GPS and leveling data.SeismologyandGeology(in Chinese), 26(3): 393-404.

Wang H, Zhang G M, Shi Y L, et al. 2006, Numerical simulation of movement and deformation of Qinghai-Tibet plateau.JournalofGeodesyandGeodynamics(in Chinese), 26(2): 15-23.

Xu S G, Xiong Y L, Liao H, et al. 2010. Static and kinematic analysis of coseismic deformation of WenchuanMS8.0 earthquake.JournalofGeodesyandGeodynamics(in Chinese), 30(3): 27-30.

Xu X W, Wen X Z, Ye J Q, et al. 2008. TheMS8.0 Wenchuan earthquake surface ruptures and its seismogenic structure.SeismologyandGeology(in Chinese), 30(3): 597-629.

Yang L Q, Du J, Chen Y. 2006. Numerical modelling of the crust/mantle deformation in the Tibetan plateau.EarthScienceFrontiers(in Chinese), 13(5): 360-373.

Zhang B, Cheng H H, Shi Y L. 2015. Calculation of the co-seismic effect ofMS8.1 earthquake, April 25, 2015, Nepal.ChineseJ.Geophys. (in Chinese), 58(5): 1794-1803, doi: 10.6038/cjg20150529.

Zhang D N, Xu Z H. 1995. An alternative interpret ation of earthquake activity caused by normal fault of upper crust in southern Qinghai-Xizang (Tibetan) plateau.ActaSeismologicaSinica(in Chinese), 17(2): 188-195. Zhang D N, Yuan S Y, Shen Z K. 2007. Numerical simulation of the recent crust movement and the fault activities in Tibetan Plateau.ChineseJ.Geophys. (in Chinese), 50(1): 153-162.

Zhang G M, Ma H S, Wang H, et al. 2005. Boundaries between active-tectonic blocks and strong earthquakes in the China mainland.ChineseJ.Geophys. (in Chinese), 48(3): 602-610. Zhang P Z, Deng Q D, Zhang G M, et al. 2003. Active tectonic blocks and strong earthquakes in the continent of China.ScienceinChinaSeriesD:EarthSciences, 46(Supppl. 2): 13-24.

Zheng H W, Li Y D, Gao R, et al. 2006. The advance of numerical simulation in geodynamics.ProgressinGeophys. (in Chinese), 21(2): 360-369.

附中文參考文獻(xiàn)

曹雪峰, 尹京苑, 楊立明. 2003. 南北地震帶應(yīng)力場(chǎng)演變的數(shù)值模擬研究. 西北地震學(xué)報(bào), 25(3): 221-225.

陳開平, 馬瑾. 1995. 印度與歐亞大陸碰撞構(gòu)造變形數(shù)值分析. 地震地質(zhì), 17(3): 277-284.

陳連旺, 張培震, 陸遠(yuǎn)忠等. 2008. 川滇地區(qū)強(qiáng)震序列庫侖破裂應(yīng)力加卸載效應(yīng)的數(shù)值模擬. 地球物理學(xué)報(bào), 51(5): 1411-1421.

陳連旺, 詹自敏, 葉際陽等. 2011. 流變特性對(duì)青藏高原構(gòu)造變形影響的數(shù)值模擬. 大地測(cè)量與地球動(dòng)力學(xué), 31(3): 8-14.

陳文彬, 徐錫偉, 張志堅(jiān)等. 2001. 2001年11月14日青新交界MS8.1地震地表破裂帶的初步調(diào)查. 西北地震學(xué)報(bào), 23(4): 313-317.

陳祖安, 林邦慧, 白武明. 2008. 1997年瑪尼地震對(duì)青藏川滇地區(qū)構(gòu)造塊體系統(tǒng)穩(wěn)定性影響的三維DDA+FEM方法數(shù)值模擬. 地球物理學(xué)報(bào), 51(5): 1422-1430.

戴黎明, 李三忠, 陶春輝等. 2011. 印度板塊擠壓驅(qū)動(dòng)龍門山斷裂帶活動(dòng)的三維數(shù)值模型. 地球物理學(xué)進(jìn)展, 26(1): 41-51, doi: 10.3969/j.issn.1004-2903.2011.01.004.

鄧宗策, 崔軍文, 馮曉楓等. 1990. 青藏高原構(gòu)造動(dòng)力及變形機(jī)制的彈性有限元分析. 中國(guó)地質(zhì)科學(xué)院院報(bào), 21: 65-67.

傅容珊, 徐耀民, 黃建華等. 2000. 青藏高原擠壓隆升過程的數(shù)值模擬. 地球物理學(xué)報(bào), 43(3): 346-355.

陸師闊, 蔡永恩. 2004. 青藏高原動(dòng)力學(xué)數(shù)值模擬方法與研究進(jìn)展. 地震學(xué)報(bào), 26(5): 547-559.

萬永革, 王敏, 沈正康等. 2004. 利用GPS和水準(zhǔn)測(cè)量資料反演2001年昆侖山口西8.1級(jí)地震的同震活動(dòng)分布. 地震地質(zhì), 26(3): 393-404.

王輝, 張國(guó)民, 石耀霖等. 2006. 青藏活動(dòng)地塊區(qū)運(yùn)動(dòng)與變形特征的數(shù)值模擬. 大地測(cè)量與地球動(dòng)力學(xué), 26(2): 15-23.

徐韶光, 熊永良, 廖華等. 2010. 汶川地震同震形變的靜態(tài)和動(dòng)態(tài)分析. 大地測(cè)量與地球動(dòng)力學(xué), 30(3): 27-30.

徐錫偉, 聞學(xué)澤, 葉建青等. 2008. 汶川MS8.0地震地表破裂帶及其發(fā)震構(gòu)造. 地震地質(zhì), 30(3): 597-629.

楊立強(qiáng), 鄧軍, 陳赟. 2006. 青藏高原殼幔形變數(shù)值模擬研究. 地學(xué)前緣, 13(5): 360-373.

張貝, 程惠紅, 石耀霖. 2015. 2015年4月25日尼泊爾MS8.1大地震的同震效應(yīng). 地球物理學(xué)報(bào), 58(5): 1794-1803, doi: 10.6038/cjg20150529.

張東寧, 徐忠淮. 1995. 青藏高原南部上地殼正斷層地震活動(dòng)的一種可能解釋. 地震學(xué)報(bào), 17(2): 188-195.

張東寧, 袁松涌, 沈正康. 2007. 青藏高原現(xiàn)代地殼運(yùn)動(dòng)與活動(dòng)斷裂帶關(guān)系的模擬實(shí)驗(yàn). 地球物理學(xué)報(bào), 50(1): 153-162.

張國(guó)民, 馬宏生, 王輝等. 2005. 中國(guó)大陸活動(dòng)地塊邊界帶與強(qiáng)震活動(dòng). 地球物理學(xué)報(bào), 48(3): 602-610.

張培震, 鄧啟東, 張國(guó)民等. 2003. 中國(guó)大陸的強(qiáng)震活動(dòng)與活動(dòng)地塊. 中國(guó)科學(xué), 33(增刊): 12-20.

鄭洪偉, 李延棟, 高銳等. 2006. 數(shù)值模擬在地球動(dòng)力學(xué)中的研究進(jìn)展. 地球物理學(xué)進(jìn)展, 21(2): 360-369.

(本文編輯胡素芳)

The distribution characteristics of deformation field caused by three great earthquakes in the Qinghai-Tibet Plateau and its vicinity since 2001

ZHANG Yuan-Sheng1,2， ZHENG Xiao-Jing1*， WANG Lan-Min2

1DepartmentofMechanics,SchoolofCivilEngineeringandMechanics,LanzhouUniversity,Lanzhou730000,China2LanzhouInstituteofSeismology,ChinaEarthquakeAdministration,Lanzhou730000,China

Based on a viscous-elastic crustal model with constraint conditions from GPS and seismic dislocation data and adopting the finite element method, we analyzed the crustal deformation characteristics caused by 2001 Kunlun MountainMS8.1, 2008 WenchuanMS8.0 and 2015 NepalMS8.1 earthquakes, obtained the displacement and deformation field. These earthquakes have occurred in the north, east, and south margins of the Tibetan Plateau respectively. Our primary results suggest obvious differences in deformation regimes, seismogenic fault geometry and source parameters, as well as crustal deformation patters despite the substantially similar magnitudes. These differences caused by seismogenic fault nature and subsurface structures are mainly manifested by different scopes of deformed regions and different strain amplitudes. Through these simulation results, we can detect more about the distribution characteristics of stress loading and unloading area caused by large earthquakes. Moreover, it plays an important role in forecasting the earthquake possibility of seismogenic region.

Earthquake； Stress loading and unloading area； Deformation field characteristics

10.6038/cjg20161011.

國(guó)家自然基金項(xiàng)目(41574044)資助.

張?jiān)校?965年生，研究員，主要從事地球深部結(jié)構(gòu)研究.E-mail:Zhangys@gssb.gov.cn

鄭曉靜，女，1958年生，中國(guó)科學(xué)院院士，發(fā)展中國(guó)家科學(xué)院院士，從事板殼非線性固體力學(xué)、電磁固體力學(xué)和風(fēng)沙環(huán)境力學(xué)研究. E-mail: xjzheng@lzu.edu.cn

10.6038/cjg20161011

P315

2016-03-14，2016-05-25收修定稿

張?jiān)?鄭曉靜， 王蘭民. 2016. 青藏高原及周邊2001年以來三次特大地震引起的形變場(chǎng)分布特征. 地球物理學(xué)報(bào)，59(10):3637-3645,

Zhang Y S， Zheng X J， Wang L M. 2016. The distribution characteristics of deformation field caused by three great earthquakes in the Qinghai-Tibet Plateau and its vicinity since 2001.ChineseJ.Geophys. (in Chinese),59(10):3637-3645,doi:10.6038/cjg20161011.