基于GPS形變資料的川滇地區應力場數值模擬研究

廖思佩　侯　強　杜永超

1　中國地質大學(武漢)機械與電子信息學院，武漢市魯磨路388號，430074

基于GPS形變資料的川滇地區應力場數值模擬研究

廖思佩1侯強1杜永超1

1中國地質大學(武漢)機械與電子信息學院，武漢市魯磨路388號，430074

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

關鍵詞：GPS資料；速度場；應力場；有限單元建模；數值模擬

研究地殼應力狀態對了解地震觸發機制、評估地震危險性具有重要意義[1]。前人基于GPS形變資料，采用有限元方法對不同區域構造應力場進行了數值模擬研究[2-4]。本文在此基礎上，采用近期GPS觀測數據和間接平差法估算測站速度，并插值得到速度場；結合地質資料，使用ABAQUS建立三維有限元模型，利用速度邊界值作為位移邊界條件約束模型，模擬得到川滇地區構造應力場，并對該地區未來地震危險性作出相應分析。

1GPS數據處理

1.1數據獲取

GPS數據來自中國地震局GNSS服務平臺提供的基礎數據產品(http:∥www.cgps.ac.cn)。中國大陸構造環境監測網絡共建設有260個基準站和2 000多個區域站，其中81個基準站和472個區域站位于研究區(96°～108°E，21°～35°N)內。服務平臺僅提供基準站的坐標時間序列，可下載的81個測站中有4個來源于中國地殼運動觀測網絡。對昆明基準站的坐標時間序列進行一元線性回歸分析，即可得到臺站所在位置的年平均速度。

1.2一元線性回歸分析

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

根據間接平差原理，對時間序列構造誤差方程：

(1)

(2)

(3)

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

(4)

1.3歐拉矢量法去除板塊整體運動

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

表1　中國大陸整體歐拉矢量

板塊運動歐拉矢量與測站站心水平運動分量存在如下關系[6]：

(5)

式中，vn、ve為測站站心水平運動的N方向和E方向速度分量，λ和φ為測站的經緯度，R表示地球半徑。

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

1.4樣條插值得到速度場

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

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

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

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

圖2　川滇地區速度場及震中分布Fig.2　Velocity field and epicenter distribution in Sichuan-Yunnan region

2川滇地區有限元模型建立

2.1主要斷裂和塊體分區

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

Ⅰ馬爾康塊體：Ⅰ1松潘-甘孜塊體；Ⅰ2秦嶺塊體Ⅱ川滇菱形塊體：Ⅱ1滇西北塊體；Ⅱ2滇中塊體Ⅲ揚子塊體　Ⅳ滇西南塊體　Ⅴ拉薩塊體圖3　川滇地區主要斷裂帶及塊體分布Fig.3　Distribution of the major fault zones and blocks partition in Sichuan-Yunnan region

2.2介質參數

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

表2　塊體介質參數

2.3網格劃分

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

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

2.4分析步驟及相互作用

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

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

2.5位移邊界條件

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

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

3模擬結果分析

3.1位移場模擬結果

圖6(a)為ABAQUS模擬的位移場結果，箭頭方向代表位移矢量方向，長短以及灰度均代表位移大小；圖6(b)為各基準站模擬結果與GPS實測值的對比，黑色箭頭為GPS實測值，灰色為模擬值。對比二者可以看出，數值模擬結果較好地反映了川滇地區西北部青藏高原推擠、東側揚子塊體阻擋的格局，同時位移方向在菱形塊體及拉薩塊體內部繞東構造結順時針旋轉。

3.2應力場模擬結果

川滇地區1970年以來發生的5.0級以上強震共613次，震源深度平均值為12.95 km，選取該深度的應力分布圖進行分析(圖7)。

圖7　川滇地區12.95 km深度應力場模擬結果Fig.7　The simulated result of stress field at 12.95 km depth in Sichuan-Yunnan region

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

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

4結語

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

參考文獻

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

[2]許才軍, 晁定波, 劉經南, 等. 用大地測量資料反演青藏高原構造應力場的初步嘗試[J]. 測繪學報, 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觀測結果研究華北地區現今構造應力場[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]劉峽, 傅容珊, 楊國華, 等. 用GPS資料研究華北地區形變場和構造應力場[J]. 大地測量與地球動力學, 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]楊元喜, 曾安敏, 吳富梅. 基于歐拉矢量的中國大陸地殼水平運動自適應擬合推估模型[J]. 中國科學:地球科學, 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]孫海燕, 丁咚. 薄板樣條函數及復雜曲面的數學表示[J]. 測繪工程, 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]徐錫偉, 聞學澤, 鄭榮章, 等. 川滇地區活動塊體最新構造變動樣式及其動力來源[J]. 中國科學:地球科學, 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]張竹琪, 張培震, 王慶良. 龍門山高傾角逆斷層結構與孕震機制[J]. 地球物理學報, 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]楊興悅, 陳連旺, 楊立明, 等. 巴顏喀拉塊體強震動力學過程數值模擬[J]. 地震學報, 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]馬宏生. 川滇地區強震孕育的深部動力環境研究[D]. 北京：中國地震局地球物理研究所, 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

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

DOI：10.14075/j.jgg.2016.07.019

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

中圖分類號：P315

文獻標識碼：A

項目來源：中央高校基本科研業務費(CUGL120234)。