廖思佩　侯　強(qiáng)　杜永超

1　中國(guó)地質(zhì)大學(xué)(武漢)機(jī)械與電子信息學(xué)院，武漢市魯磨路388號(hào)，430074

基于GPS形變資料的川滇地區(qū)應(yīng)力場(chǎng)數(shù)值模擬研究

廖思佩1侯強(qiáng)1杜永超1

1中國(guó)地質(zhì)大學(xué)(武漢)機(jī)械與電子信息學(xué)院，武漢市魯磨路388號(hào)，430074

摘要：基于近10 a的GPS坐標(biāo)時(shí)間序列獲得川滇地區(qū)速度場(chǎng)，結(jié)合地質(zhì)資料將該區(qū)劃分為7個(gè)塊體，并使用ABAQUS建立三維有限元模型，提取速度場(chǎng)邊界值作為模型約束條件，模擬得到應(yīng)力場(chǎng)。位移場(chǎng)數(shù)值模擬結(jié)果顯示，模型較好地反映了川滇地區(qū)西北部青藏高原推擠、東側(cè)揚(yáng)子塊體阻擋的格局，同時(shí)位移方向在菱形塊體及拉薩塊體內(nèi)部繞東構(gòu)造結(jié)順時(shí)針旋轉(zhuǎn)，與實(shí)際相符。應(yīng)力場(chǎng)模擬結(jié)果表明，龍門山斷裂帶、鮮水河斷裂帶和安寧河斷裂帶交界區(qū)域南北兩側(cè)、小金河-麗江斷裂帶西段與瀾滄江斷裂帶交界處南北側(cè)、南汀河斷裂帶西段南北側(cè)、紅河斷裂帶中段東西側(cè)均存在明顯的應(yīng)力分布不均現(xiàn)象，地震危險(xiǎn)性較大。

關(guān)鍵詞：GPS資料；速度場(chǎng)；應(yīng)力場(chǎng)；有限單元建模；數(shù)值模擬

研究地殼應(yīng)力狀態(tài)對(duì)了解地震觸發(fā)機(jī)制、評(píng)估地震危險(xiǎn)性具有重要意義[1]。前人基于GPS形變資料，采用有限元方法對(duì)不同區(qū)域構(gòu)造應(yīng)力場(chǎng)進(jìn)行了數(shù)值模擬研究[2-4]。本文在此基礎(chǔ)上，采用近期GPS觀測(cè)數(shù)據(jù)和間接平差法估算測(cè)站速度，并插值得到速度場(chǎng)；結(jié)合地質(zhì)資料，使用ABAQUS建立三維有限元模型，利用速度邊界值作為位移邊界條件約束模型，模擬得到川滇地區(qū)構(gòu)造應(yīng)力場(chǎng)，并對(duì)該地區(qū)未來地震危險(xiǎn)性作出相應(yīng)分析。

1GPS數(shù)據(jù)處理

1.1數(shù)據(jù)獲取

GPS數(shù)據(jù)來自中國(guó)地震局GNSS服務(wù)平臺(tái)提供的基礎(chǔ)數(shù)據(jù)產(chǎn)品(http:∥www.cgps.ac.cn)。中國(guó)大陸構(gòu)造環(huán)境監(jiān)測(cè)網(wǎng)絡(luò)共建設(shè)有260個(gè)基準(zhǔn)站和2 000多個(gè)區(qū)域站，其中81個(gè)基準(zhǔn)站和472個(gè)區(qū)域站位于研究區(qū)(96°～108°E，21°～35°N)內(nèi)。服務(wù)平臺(tái)僅提供基準(zhǔn)站的坐標(biāo)時(shí)間序列，可下載的81個(gè)測(cè)站中有4個(gè)來源于中國(guó)地殼運(yùn)動(dòng)觀測(cè)網(wǎng)絡(luò)。對(duì)昆明基準(zhǔn)站的坐標(biāo)時(shí)間序列進(jìn)行一元線性回歸分析，即可得到臺(tái)站所在位置的年平均速度。

1.2一元線性回歸分析

對(duì)于某一測(cè)站的GPS時(shí)間序列，設(shè)含n天單日坐標(biāo)解，ti(i=1,2,…,n)為第i天的歷元，yi(i=1,2,…,n)為ti歷元的位移觀測(cè)值，采用直線y=at+b對(duì)該位移時(shí)間序列進(jìn)行擬合，其中a(為擬合直線的斜率，即測(cè)站的年平均速度)和b是回歸參數(shù)。

根據(jù)間接平差原理，對(duì)時(shí)間序列構(gòu)造誤差方程：

(1)

(2)

(3)

式中，NA=ATPA，U=ATPl。求解該方程得：

(4)

1.3歐拉矢量法去除板塊整體運(yùn)動(dòng)

利用上述方法估算得到ITRF2008框架下各GPS測(cè)站的年平均速度值，但該值包含大陸地殼運(yùn)動(dòng)的整體速度，應(yīng)予以消除。國(guó)際上較為通用的NNR_Nuvel1A模型[5]主要采用國(guó)外數(shù)據(jù)，中國(guó)大陸內(nèi)部變形存在許多局部特征，據(jù)此楊元喜等[6]提出適于中國(guó)大陸速度場(chǎng)的自適應(yīng)擬合推估法，并給出具體的板塊運(yùn)動(dòng)歐拉矢量，如表1所示。

表1　中國(guó)大陸整體歐拉矢量

板塊運(yùn)動(dòng)歐拉矢量與測(cè)站站心水平運(yùn)動(dòng)分量存在如下關(guān)系[6]：

(5)

式中，vn、ve為測(cè)站站心水平運(yùn)動(dòng)的N方向和E方向速度分量，λ和φ為測(cè)站的經(jīng)緯度，R表示地球半徑。

將表1中楊元喜等計(jì)算得到的歐拉矢量和81個(gè)測(cè)站的經(jīng)緯度代入式(5)，得到每個(gè)測(cè)站的整體速度，并用回歸分析中計(jì)算的年平均速度值減去該速度，得到各測(cè)站的相對(duì)速度。

1.4樣條插值得到速度場(chǎng)

薄板樣條插值以最小曲率面來擬合控制點(diǎn)，具有連續(xù)、光滑的數(shù)學(xué)特性，且計(jì)算簡(jiǎn)便[7]。使用Matlab提供的曲面擬合工具箱中的薄板樣條函數(shù)對(duì)81個(gè)測(cè)站的N方向和E方向速度進(jìn)行插值，結(jié)果見圖1(單位:mm/a)。

圖1　川滇地區(qū)N方向和E方向速度場(chǎng)等值線Fig.1　Velocity field of direction N and E in Sichuan-Yunnan region

由圖1(a)看出，小金河-麗江斷裂帶北側(cè)附近區(qū)域北向速度呈現(xiàn)高負(fù)值異常，峰值達(dá)-16 mm/a；整個(gè)川滇菱形地塊及其西側(cè)部分區(qū)域的南向速度幾乎都在10 mm/a以上，明顯大于周圍的南向速度，呈現(xiàn)出向南逃逸的趨勢(shì)；菱形地塊西北區(qū)域速度方向朝北，其原因是青藏高原和印度板塊的擠壓作用。圖1(b)中西北部東向速度非常大，主要是受青藏高原北東向推擠的影響，由于穩(wěn)定的揚(yáng)子塊體對(duì)其阻擋，東向速度在岷江斷裂、龍門山斷裂附近急速減小，以致于在西南部呈現(xiàn)反向運(yùn)動(dòng)；龍門山斷裂帶中段存在明顯的高速異常；南汀河斷裂帶中段以及紅河斷裂帶南段則存在一定的低速異常。

將N方向和E方向速度進(jìn)行綜合，得到總的速度大小及方向,如圖2(圖2(b)中圓形代表川滇地區(qū)6級(jí)以上地震，星形代表7級(jí)以上地震,速度等值線單位:mm/a)所示。研究區(qū)域西北部具有較大的東向和北東向速度，至揚(yáng)子塊體急速減小，并在川滇菱形塊體內(nèi)部逐漸向南偏轉(zhuǎn)。橫向觀察研究區(qū)域，呈現(xiàn)西大東小的格局；縱觀菱形塊體及其西側(cè)區(qū)域，位移方向自北向南逐漸偏轉(zhuǎn)，顯示了以東構(gòu)造結(jié)(26°N，96°E)為軸旋轉(zhuǎn)的特性。地震多數(shù)發(fā)生于速度變化強(qiáng)烈、梯度明顯的區(qū)域。

圖2　川滇地區(qū)速度場(chǎng)及震中分布Fig.2　Velocity field and epicenter distribution in Sichuan-Yunnan region

2川滇地區(qū)有限元模型建立

2.1主要斷裂和塊體分區(qū)

圖3標(biāo)注了川滇地區(qū)的主要活動(dòng)斷裂帶[8]?；趫D3，選取分區(qū)控制點(diǎn)，要求既要能反映模型主要特征，又不可過于簡(jiǎn)略以致與實(shí)際不符;連接各控制點(diǎn)逐個(gè)建立塊體，同時(shí)垂向拉伸20 km。模型中除龍門山斷裂帶以外均設(shè)置為直立斷層，考慮龍門山斷裂帶淺陡深緩的特征(上地殼淺部?jī)A角約有60°～80°，深部則放緩至30°～40°)，本文設(shè)置龍門山斷裂傾角為70°[9]。

Ⅰ馬爾康塊體：Ⅰ1松潘-甘孜塊體；Ⅰ2秦嶺塊體Ⅱ川滇菱形塊體：Ⅱ1滇西北塊體；Ⅱ2滇中塊體Ⅲ揚(yáng)子塊體?、舻嵛髂蠅K體　Ⅴ拉薩塊體圖3　川滇地區(qū)主要斷裂帶及塊體分布Fig.3　Distribution of the major fault zones and blocks partition in Sichuan-Yunnan region

2.2介質(zhì)參數(shù)

建模時(shí)選用彈性模型，同時(shí)需對(duì)塊體進(jìn)行介質(zhì)分區(qū)，斷裂帶作為弱化帶處理，其楊氏模量為所在塊體楊氏模量的1/3。各塊體介質(zhì)參數(shù)如表2所示[1, 10-11]。

表2　塊體介質(zhì)參數(shù)

2.3網(wǎng)格劃分

ABAQUS劃分網(wǎng)格有多種策略，本文選擇全局種子尺寸為 30 km，將斷裂附近邊界種子尺寸適當(dāng)縮小。采用自由網(wǎng)格劃分技術(shù)，通過四面體單元來劃分網(wǎng)格，共建立了29 867個(gè)網(wǎng)格、56 898個(gè)節(jié)點(diǎn)，如圖4所示。

圖4　川滇地區(qū)上地殼三維有限元模型Fig.4　Three dimensional finite element model of upper crust in Sichuan-Yunnan region

2.4分析步驟及相互作用

由于時(shí)間在靜態(tài)分析過程中沒有實(shí)際物理意義，因此將時(shí)間長(zhǎng)度定為1。一般來說，可建立3個(gè)分析步驟：首先添加模型的初始應(yīng)力場(chǎng)以及一直存在的載荷或約束；然后建立接觸，使模型能夠平穩(wěn)過渡；最后開始加載邊界條件，此時(shí)將初始時(shí)間增量定為0.1，這樣就可以使載荷平滑地加載到模型中，從而避免了計(jì)算不收斂的問題。

定義相互作用時(shí)首先選定作用屬性為“接觸”，接下來添加切向行為，即沿接觸面的行為，選用罰函數(shù)摩擦模型，并將摩擦系數(shù)設(shè)置為0.6。

2.5位移邊界條件

由GPS數(shù)據(jù)處理結(jié)果可得到川滇地區(qū)速度場(chǎng)，選取模型各個(gè)邊界上的網(wǎng)格節(jié)點(diǎn)作為控制點(diǎn)，通過樣條插值得到這些點(diǎn)的速度值，并逐一添加至模型中，如圖5所示。考慮時(shí)間步長(zhǎng)，模型設(shè)置位移條件而非速度條件，邊界約束施加的是年位移量[11]。模型上表面為自由表面，底部垂直方向固定，水平方向自由。

圖5　模型邊界條件Fig.5　Boundary conditions of the model

3模擬結(jié)果分析

3.1位移場(chǎng)模擬結(jié)果

圖6(a)為ABAQUS模擬的位移場(chǎng)結(jié)果，箭頭方向代表位移矢量方向，長(zhǎng)短以及灰度均代表位移大??；圖6(b)為各基準(zhǔn)站模擬結(jié)果與GPS實(shí)測(cè)值的對(duì)比，黑色箭頭為GPS實(shí)測(cè)值，灰色為模擬值。對(duì)比二者可以看出，數(shù)值模擬結(jié)果較好地反映了川滇地區(qū)西北部青藏高原推擠、東側(cè)揚(yáng)子塊體阻擋的格局，同時(shí)位移方向在菱形塊體及拉薩塊體內(nèi)部繞東構(gòu)造結(jié)順時(shí)針旋轉(zhuǎn)。

3.2應(yīng)力場(chǎng)模擬結(jié)果

川滇地區(qū)1970年以來發(fā)生的5.0級(jí)以上強(qiáng)震共613次，震源深度平均值為12.95 km，選取該深度的應(yīng)力分布圖進(jìn)行分析(圖7)。

圖7　川滇地區(qū)12.95 km深度應(yīng)力場(chǎng)模擬結(jié)果Fig.7　The simulated result of stress field at 12.95 km depth in Sichuan-Yunnan region

地震往往發(fā)生在應(yīng)力分布不均勻區(qū)域。在等效應(yīng)力云圖(圖7(a))中，龍門山斷裂帶、鮮水河斷裂帶和安寧河斷裂帶交界的“Y”型結(jié)構(gòu)南側(cè)的石棉到冕寧區(qū)域(101.5°～102.5°E，27.5°～30°N)有明顯的應(yīng)力集中現(xiàn)象，但東西兩側(cè)應(yīng)力相當(dāng)。需要注意的是，北側(cè)應(yīng)力相較于南側(cè)小了近1個(gè)量級(jí)，導(dǎo)致龍門山斷裂帶和鮮水河斷裂帶兩側(cè)應(yīng)力分布不均勻，成為地震頻發(fā)的一個(gè)因素。小金河-麗江斷裂帶西段與瀾滄江斷裂帶交界處同樣存在應(yīng)力集中現(xiàn)象，且應(yīng)力主要集中在北側(cè)，南側(cè)應(yīng)力仍然較小。這種不均勻性易誘發(fā)地震。從歷史強(qiáng)震分布來看，該區(qū)地震較少，可能是未來繼龍門山斷裂帶之后又一地震危險(xiǎn)區(qū)。南汀河斷裂帶西段南北兩側(cè)應(yīng)力梯度變化較大，紅河斷裂帶中段東西兩側(cè)同樣存在梯度明顯變化的區(qū)域，這些地方應(yīng)加強(qiáng)監(jiān)測(cè)，是地震預(yù)測(cè)的重點(diǎn)觀察區(qū)域。

觀察川滇地區(qū)最大主應(yīng)力方向(圖7(b))，同樣在“Y”型結(jié)構(gòu)下方呈現(xiàn)高壓區(qū)；拉薩塊體東側(cè)與小金河-麗江斷裂帶西段交界區(qū)域存在一定東西方向壓力，通過瀾滄江斷裂帶在北側(cè)形成近南北方向張力；滇西南塊體整體呈現(xiàn)拉張趨勢(shì)，與拉薩塊體交界處的南汀河斷裂帶西段北側(cè)承受著較大的壓力，在通過斷裂帶擠壓過程中轉(zhuǎn)化為張力；紅河斷裂帶中段存在來自滇中塊體的南向壓力，雖然數(shù)值較小，但也不容忽視。

4結(jié)語

本文使用中國(guó)地震局GNSS數(shù)據(jù)產(chǎn)品服務(wù)平臺(tái)提供的川滇地區(qū)81個(gè)基準(zhǔn)站的GPS坐標(biāo)時(shí)間序列數(shù)據(jù)，采用間接平差法擬合各站點(diǎn)N方向和E方向速度值，并使用歐拉矢量法去除板塊整體運(yùn)動(dòng)，最后利用樣條插值得到川滇地區(qū)速度場(chǎng)。根據(jù)塊體分區(qū)情況，建立符合川滇地區(qū)地質(zhì)特點(diǎn)的三維有限元模型，并根據(jù)速度場(chǎng)樣條插值結(jié)果設(shè)置模型邊界位移條件，將數(shù)值模擬的位移場(chǎng)與GPS實(shí)測(cè)結(jié)果對(duì)比，以此驗(yàn)證模型的合理性，最終得到應(yīng)力場(chǎng)數(shù)值模擬結(jié)果。結(jié)果顯示，川滇地區(qū)存在多處應(yīng)力集中區(qū)域，特別是龍門山斷裂帶、鮮水河斷裂帶和安寧河斷裂帶交界的“Y”型結(jié)構(gòu)南側(cè)、小金河-麗江斷裂帶西段與瀾滄江斷裂帶交界處北側(cè)。根據(jù)歷史強(qiáng)震分布，地震往往發(fā)生在速度變化明顯、應(yīng)力分布不均的地方，龍門山斷裂帶、小金河-麗江斷裂帶、南汀河斷裂帶西段以及紅河斷裂帶中段均有一定的地震危險(xiǎn)性，應(yīng)予以重視。本文的后續(xù)工作是精細(xì)化模型，并加入中下地殼和上地幔分層結(jié)構(gòu)，綜合評(píng)價(jià)川滇地區(qū)未來的地震發(fā)展趨勢(shì)。

參考文獻(xiàn)

[1]李紅. 構(gòu)造應(yīng)力場(chǎng)、活動(dòng)斷裂及區(qū)域地震活動(dòng)性的數(shù)值模擬研究[D]. 北京：中國(guó)地震局地殼應(yīng)力研究所, 2008(Li Hong. Numerical Simulation on Tectonic Stress Field,Active Fault and Regional Seismicity[D]. Beijing: Institute of Crustal Dynamics, CEA, 2008)

[2]許才軍, 晁定波, 劉經(jīng)南, 等. 用大地測(cè)量資料反演青藏高原構(gòu)造應(yīng)力場(chǎng)的初步嘗試[J]. 測(cè)繪學(xué)報(bào), 1997,26(2):3-8(Xu Caijun, Chao Dingbo, Liu Jingnan, et al. A Preliminary Attempt for Inversion of Tectonic Stress Field in Tibetan Plateau with Geodetic Data[J]. Acta Geodaetica et Cartographica Sinica, 1997,26(2):3-8)

[3]徐菊生, 袁金榮, 高士鈞, 等. 利用GPS觀測(cè)結(jié)果研究華北地區(qū)現(xiàn)今構(gòu)造應(yīng)力場(chǎng)[J]. 地殼形變與地震, 1999,19(2):83-91 (Xu Jusheng, Yuan Jinrong, Gao Shijun, et al. Research on the Precent Tectonic Stress Field in North China with GPS Observations[J]. Crustal Deformation and Earthquake, 1999,19(2):83-91)

[4]劉峽, 傅容珊, 楊國(guó)華, 等. 用GPS資料研究華北地區(qū)形變場(chǎng)和構(gòu)造應(yīng)力場(chǎng)[J]. 大地測(cè)量與地球動(dòng)力學(xué), 2006,26(3):33-39 (Liu Xia, Fu Rongshan, Yang Guohua, et al. Deformation Field and Tectonic Stress Field Constrained by GPS Observation in North China[J]. Journal of Geodesy and Geodynamics, 2006,26(3):33-39)

[5]Argus D F, G G R. No-Net-Rotation Model of Current Plate Velocities Incorporating Plate Motion Model NUVEL-1[J]. Geophysical Research Letters, 1991,18(11):2 038-2 042

[6]楊元喜, 曾安敏, 吳富梅. 基于歐拉矢量的中國(guó)大陸地殼水平運(yùn)動(dòng)自適應(yīng)擬合推估模型[J]. 中國(guó)科學(xué):地球科學(xué), 2011,41(8):1 116-1 125(Yang Yuanxi, Zeng Anmin, Wu Fumei. Horizontal Crustal Movement in China Fitted by Adaptive Collocation with Embedded Euler Vector[J]. Science China:Earth Science, 2011,41(8):1 116-1 125)

[7]孫海燕, 丁咚. 薄板樣條函數(shù)及復(fù)雜曲面的數(shù)學(xué)表示[J]. 測(cè)繪工程, 2006,15(2):7-8(Sun Haiyan, Ding Dong. Thin Plate Spline and Mathematic Description of Complicated Surface[J]. Engineering of Surveying and Mapping, 2006,15(2):7-8)

[8]徐錫偉, 聞學(xué)澤, 鄭榮章, 等. 川滇地區(qū)活動(dòng)塊體最新構(gòu)造變動(dòng)樣式及其動(dòng)力來源[J]. 中國(guó)科學(xué):地球科學(xué), 2003,33(增1):151-162(Xu Xiwei, Wen Xueze, Zheng Rongzhang, et al. Pattern of Latest Tectonic Motion and Its Dynamics for Active Blocks in Sichuan-Yunnan Region[J]. Science China: Earth Science, 2003,33(Supp1):151-162)

[9]張竹琪, 張培震, 王慶良. 龍門山高傾角逆斷層結(jié)構(gòu)與孕震機(jī)制[J]. 地球物理學(xué)報(bào), 2010,53(9):2 068-2 082(Zhang Zhuqi, Zhang Peizhen, Wang Qingliang. The Structure and Seismogenic Mechanism of Longmenshan High Dip-Angle Reverse Fault[J]. Chinese J Geophys, 2010, 53(9):2 068-2 082)

[10]楊興悅, 陳連旺, 楊立明, 等. 巴顏喀拉塊體強(qiáng)震動(dòng)力學(xué)過程數(shù)值模擬[J]. 地震學(xué)報(bào), 2013,35(3):304-314(Yang Xingyue, Chen Lianwang, Yang Liming, et al. Numerical Simulation on Strong Earthquake Dynamic Process of Bayan Har Block[J]. Acta Seismologica Sinica, 2013, 35(3):304-314)

[11]馬宏生. 川滇地區(qū)強(qiáng)震孕育的深部動(dòng)力環(huán)境研究[D]. 北京：中國(guó)地震局地球物理研究所, 2007(Ma Hongsheng. Dynamics Research on Strong Shock Gestation in Sichuan-Yunnan and Its Adjacent Areas[D]. Beijing: Institute of Geophysics, CEA, 2007)

Foundation support:Fundamental Research Funds for Central Universities, No.CUGL120234.

About the first author:LIAO Sipei, postgraduate, majors in seismic risk, stress field and crust-mantle coupling, E-mail：liaosipei_h@163.com.

Numerical Simulation of the Stress Field in Sichuan-Yunnan Region with GPS Deformation Data

LIAOSipei1HOUQiang1DUYongchao1

1School of Mechanical Engineering and Electronic Information, China University of Geosciences (Wuhan),388 Lumo Road, Wuhan 430074, China

Abstract:Based on GPS coordinate time series from this decade, the velocity field of the Sichuan-Yunnan region is obtained. According to geological data, we divide this region into 7 blocks, then find the three-dimensional finite element model by ABAQUS. Through extracting the boundary value of the velocity field as model constraint conditions, we obtain numerical simulation results of displacement and stress field. The results show that our model clearly reflects the situation of pushing forces from Qinghai-Tibet plateau in the northwest and prevention from Yangzi block in the east. The displacement direction rotates clockwise on eastern Himalayan syntaxis inside rhombic block and Lhasa block. These results conform with reality. Simulation results of stress field indicate that the stress distribution is obviously uneven in some areas, such as the border region of Longmanshan, Xianshuihe and Anninghe fault zone, the border region of Xiaojinhe-Lijiang and Lancangjiang fault zone, the south and north flanks of Nantinghe fault zone, as well as the east and west flanks of Red River fault zone. Seismic risk in these areas is larger than other regions.

Key words:GPS data; velocity field; stress field; finite element modeling; numerical simulation

收稿日期：2015-07-28

第一作者簡(jiǎn)介：廖思佩，女，碩士生，主要從事地震危險(xiǎn)性分析、應(yīng)力場(chǎng)和殼幔耦合研究，E-mail：liaosipei_h@163.com。

DOI：10.14075/j.jgg.2016.07.019

文章編號(hào)：1671-5942(2016)07-0645-05

中圖分類號(hào)：P315

文獻(xiàn)標(biāo)識(shí)碼：A

項(xiàng)目來源：中央高?；究蒲袠I(yè)務(wù)費(fèi)(CUGL120234)。