基于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)。