利用GPS资料求解中国大陆应变场的新方法*

段维波 吴 云 陈慧杰

(中国地震局地震研究所，武汉 430071)

利用GPS资料求解中国大陆应变场的新方法*

段维波 吴 云 陈慧杰

(中国地震局地震研究所，武汉 430071)

探讨利用边长变化方法解算区域应变场，并用该方法对1999—2001年中国大陆及周边区域网GPS测站边长变化量进行分析，给出了该期间的中国大陆应变率场变化。分析结果表明：对大范围、密集的GPS测站的边长变化量进行连续变形分析，既能揭示中国大陆地壳应变的总体分布，又能显示局部地区的构造获得。总体上，中国大陆构造的水平变形强度在南北地震带产生突变，呈现出西强东弱、南强北弱的态势。而昆仑山地块中部、鲜水河断裂带和云南中部地区，其应变速率最大。天山东部具有拉张环境。

面应变;应变率;主应变;区域网;GPS

1 引言

运用GPS观测数据研究地壳形变已越来越受到人们的关注［1，2］。通常利用区域大地测量网包括GPS网的数据，可以根据两期测量得到的坐标变化即位移来计算地壳应变率场［3］。目前有多位学者对中国大陆的应变率场分布进行了研究［4-11］。这些方法可以概括为局部方法和整体方法。局部算法是先划分网格(或块体)并按网格(或块体)计算应变率;整体算法是将速度场建立为点位(经纬度或高斯坐标)的某种函数，然后对该方法求偏导获取应变率场分布。虽然利用上述算法得到的应变率场可以从总体上反映整个研究区域的变形状态，但不同的方法即使利用相同数据也很难得到完全相同的结果，所以有必要探索新的、有效的应变场计算方法。

鉴于此，本文对三角形方法进行扩展，即在区域划分时不按三角形进行划分而是按长方形(或者正方形)进行划分(称之为长方形法)，利用边长的变化量计算应变场。

2 数据来源及计算结果

处理中的数据来源于中国大陆陆态观测网络数据中心以及IGS数据中心提供的79个GPS测站的观测数据(图1)，使用GAMIT/GLOBK软件(10.35版)在ITRF2000坐标框架下采用IGS精密星历进行GPS基线数据处理。基线处理控制文件取自从SIO网站下载的sestb1.和autcln.cmd文件。

图1 GPS测站分布Fig.1 Distribution of GPS sites in China

数据处理控制参数选择如下：

采样间隔：30 s;

最大历元数：2 880;

截至高度角：10°;

基线处理：松弛解;

对流层误差模型：Saastamoinen模型和缺省气象参数;

天顶延迟参数个数：25;

电离层延迟：消除电离层的LC观测量;

卫星轨道参数：9参数(6个开普勒根数和3个太阳辐压参数);

数据清理：AUTCLN。

数据处理中使用ITRF2000框架下的核心站的坐标和坐标运动速率为坐标系统控制参数，给予强约束，即东西向和南北向坐标约束为3 mm，高程向为10 mm。对其余的GPS观测站的坐标和坐标运动速率给予松约束，即东西向和南北向坐标约束为30 mm，高程向为100 mm。

研究中两个测段的基线正规解在Q文件中的NRMS值均小于0.3［12］。基线结果如表1所示(由于基线条数较多，篇幅限制，只列出其中部分基线长度)。

表1 基线长度(单位：km)Tab.1 Length of each baseline(unit：km)

3 根据边长变化量求解应变

如图2所示，在直角坐标系中，S为变形前的边长，对应的基线向量为(ΔX，ΔY，ΔZ)T，其中ΔX、ΔY、ΔZ为基线的坐标分量;S'为变形后的边长;dS、dS'分别为沿边长 S和S'的微元，由弹性力学理论［13］：

式中EKL(X，t)是拉格朗日应变张量。对公式两边沿边长S积分，假定应变张量不变，则EKL(X，t)为常量，得到：

顾及

且EKL为对称张量，式(2)可以展开成

进一步假设边长的变化量远远小于边长，即：

并且引入边长应变ε=ΔS/S，式(4)可以改写为：

这里，α、β、γ是边长S变形前的方向余弦，即：

如果仅考虑二维平面应变，式(6)可写成：

如果我们引入边长的方位角φ来代替α、β，则式(8)可写成：

图2 边长变化示意图Fig.2 Change of a side length

3 根据四边形法模型计算应变场

如图3所示，有19个测站，形成的基线条数为171条，由于每条基线在每个方框里面所分配的边长长度及其变化量无法得知，我们采用如下方法获得其值，即假设E018～JB15基线穿过了m个方框(此时m的取值与方框的大小、测站E018和测站JB15的位置有关)每个方框的长度变化为ΔSi，其中(1≤i≤m)，那么边长的总变化量为：

图3 四边形法模型展示Fig.3 Show og the model of quadrilateral method

即某条基线对应的线应变为：

令：

将上式改写为矩阵形式：AX=ε

s11，s12，…，sCn1为基线1，2，…，穿过第一个方框时分配的长度，s12，s22，…，sCn2为基线1，2，…，穿过第二个方框时分配的长度，依次类推。

每条基线穿过哪几个方框，在穿过的方框中所分配的长度，ACn×M可根据程序计算出来，故ACn×M已知。

其中，Cn=n(n-1)/2，M=k×l×3。n为测站个数，k为纬度划分个数，l为经度划分个数。M则为区域被划分后的方框总个数。

对于多条基线有：

利用最小二乘法可求解超定方程组(13)，可求解出每个方框的(E11，E22，E12)

式中ε1为最大主应变，ε2为最小主应变。

根据求出的ε1，2、rmax、εarea画出的等值线图如图4、5、6所示。

图4显示的是172个3°×3°网格内平均最大主应变率。该应变率图表明，构造应力相对集中于青藏及其周缘的断裂带、天山和川滇地区。整个形变区内，以喜马拉雅构造带的挤压应变为重要特征，拉张的成分很小。天山地区的应变基本类似于喜马拉雅，也具有明显的挤压特征。其他地壳的应变一般同时含有挤压和拉张成分。从主应变率的大小分析，塔里木盆地，阿拉善以及华南地球的平均主压应变率很小，即相邻点位之间不发生相对位移，其内部结构完整。青藏高原内部的拉萨、羌塘、祁连山地块、阿尔金断裂带和天山，其内部相邻点位之间发生了相对运动，内部也积累了一定的构造应力。就其范围而言，形变、应变规模不大。地块内部变形的方向和幅度之差远远小于地块间的差异。例如具有地块边界性质的鲜水河断裂带、昆仑山断裂和喜马拉雅，相邻点位之间相对位移比较显著，构造应力相对集中，内部构造活动强烈，也是强震最为频繁的地带。

整个区域的面膨胀率分布见图5，面膨胀率的正值表示地壳伸展，图中用实线表示;负值表示地壳压缩，图中用虚线表示。此图反映的局部特征十分明显，地壳水平拉伸和压缩一目了然。面膨胀的高压缩区仍在喜马拉雅一带，昆仑山一带水平压缩高值区和水平伸展高值区交替出现。川滇地区则主要呈现的是地壳拉伸状态，华南地区则有水平压缩低值区与拉伸低值区交替出现。整个天山地区以压缩为主，但是天山东部有小部分地区呈现水平拉伸态势，这与该地区的挤压构造特征有所出入。由于GPS更多地反映现今变形，因此，GPS得到的应变场与长期构造应变不一致是可以由多种原因的，就GPS观测而言，更多的观测，更精确的结果，也许能为进一步澄清两者的不同提供帮助。

中国大陆地区的最大剪应变率(图6)以南北地震带为界，呈西强东弱态势。藏南，羌塘地块中西部部分地区以及新疆西北地区属最大见应变率高值区。其中藏南地区最大剪应变率量级为24×10-8/ a;天山东部、安宁河-小江断裂带以及西藏中部部分地区的最大剪应变率相对也较强，量级为14× 10-8/a。塔里木盆地地区的剪应变率表现的很小，最大值约为7.2×10-8/a。其中华南地块的剪应变率最小，其最大剪应变率小于4.8×10-8/a。

图4 研究区域最大、最小主应变Fig.4 Maxium and minium main strains of China mainland

图5 研究区域面膨胀率结果Fig.5 Surface expansion ratio of China mainland

图6 研究区域剪应变率结果Fig.6 Shear strain rate of China mainland

4 结论

1)利用中国大陆及周边的GPS观测资料构建的中国大陆形变场，各类应变场。在数千千米的尺度上得到的应变场真实，直观，丰富，与其他资料间接得到的同类结果相比，具有更为坚实的资料基础，并少受主观假设的制约。比较精细地展示了中国大陆的构造变形与应变率总特点，即西强东弱，南强北弱态势。与长期以来形成的认识大体一致，并用大地测量观测的角度进一步加深，细化这种基本认识，突出了现今地壳观测对中国大陆应变场变化研究的作用。当然，从另外一方面来说，GPS资料解算的结果不一定与地质时间尺度的平均运动完全吻合。两者的不同既有可能是地质资料的经度、代表性方面的问题，也有可能是GPS观测的时间太短，精度不够，更有可能是地壳变动的周期变化所致。

2)通过对1999—2001年区域网位移变化结果的应变场分析发现，中国大陆应变场表现了西强东弱的特性。面应变方面，以昆仑山地震震源区为界表现出西张东压的特性，符合震源机制解。最大剪应变方面，以昆仑山地震震源区附近表现为高值特性。对于东部地区，应变的值量普遍比较低。

3)此方法是对Delauney三角形法的改进，利用长方形的方法对研究区域进行等面积划分，克服了三角形法中由于三角形的形状、大小误差的影响。研究结果显示，观测点的分布以及观测点的数量对计算结果有一定的影响，分布情况越好，形成的网形越均匀，越有利对各地区的研究，观测点数量越多，划分的方框越多，越有利于得到各小区域的应变积累情况。

4)本研究在连续形变假设条件下，在结合前人研究的基础上，探索出一种新的利用GPS的位移变化量的四边形法求解中国大陆的形变场、应变场。虽然拟合结果受观测资料的精度、分布密度、分布均匀性等因素的影响，但是仍然能够较客观地反映空间分布特征和一般事实。可见随着GPS观测精度的进一步提高和观测资料的日益丰富，通过构建合适的反演模型，利用其进行地壳构造变形和地球动力学研究，以及地震预测的应用价值将进一步体现。

1 Wu Jicang，et al.Research on an intraplate movement model by inversion of GPS data in North China［J］.J Geodyn.，2001，31(5)：507-518.

2 Wu Yun，et al.The initial discussion of the crust deformation of China and its adjacent area based on GPS measurements［J］.Acta Seismologica Sinca，1999，21(5)：545-553.

3 伍吉仓，邓康伟，陈永奇.用边长变化结果计算地应变及块体划分［J］.大地测量与地球动力学，2002，(5)：35-38.

4 杨少敏，等.用双三次样条函数和GPS资料反演现今中国大陆构造形变场［J］.大地测量学与地球动力学，2002，(1)：68-75.

5 杨国华，等.由GPS观测结果推到中国大陆现今水平应变场［J］.地震学报，2002，24(4)：337-347.

6 沈正康，等.中国大陆现今构造应变率场及其动力学成因研究［J］.地学前缘，2003，10(增)：93-100.

7 张静华，等.用GPS测量结果研究华北现今构造形变场［J］.大地测量学与地球动力学，2004，(3)：40-46.

8 游新兆，等.中国大陆地壳现今运动的GPS测量结果与初步分析［J］.地壳形变与地震，2001，(3)：1-8.

9 朱守彪，等.利用广义应变花法计算地应变场［J］.防灾技术高等专科学校学报，2006，5(2)：1-5.

10 武艳强，等.利用多面函数整体求解GPS应变场的方法及应用［J］.武汉大学学报(信息科学版)，2009，34(9)：1 085-1 087.

11 武艳强，等.利用最小二乘配置在球面上整体解算GPS应变场的方法及应用［J］.地球物理学报，2009，52(7)：1 707-1 714.

12 王留朋，等.GAMIT在LINUX操作系统上的安装和使用［J］.全球定位系统，2006，(6)：43-46.

13 Eringen A C.Mechanics of continua［M］.John Wiley＆Sona，Inc.，1980.

14 杨本立，等.超定方程组最小二乘行处理法［J］.云南师范大学学报，1997，17(1)：1-4.

A NEW METHOD TO SOLVE STRAIN FIELD OF CHINA MAINLAND BY USING GPS DATA

Duan Weibo，Wu Yun and Chen Huijie
(Institude of Seismology，CEA，Wuhan 430071)

A new method for calculating strains based on side length changes is proposed.By use of this method the strains distribution of Chinese mainland in the period of 1999 to 2001 are calculated.The results show that the analysis of the continuous deformation in the large-scale and dense GPS velocity field can reveal not only the integral tectonic characters of Chinese mainland but also the tectonic deformation which have a mutation in the South-North Seismic Belt and is stronger in the west than in the east and stronger in the south than in the north.The largest strain rates are found in the areas as Kunlun block，Xianshuihe fault zone and central Yunnan，and the variation of velocity is very rapid，and extensive strain is found in the eastern part of Tianshan.

surface strain;strain rate;main strain;regional network;GPS

1671-5942(2011)Supp.-0072-05

2011-02-03

中国地震局地震研究所所长基金(IS201026031)

段维波，男，1985年生，硕士，主要从事GPS数据处理与地壳应变分析.E-mail：weibo_9999@163.com

P227;P315.72+6

A