GNSS 坐标时间序列在地表形变监测中的应用

王龙超

(中国铁路设计集团有限公司，天津 300251)

随着GNSS 观测站多年的持续运行，各个观测网络已经积累了大量的观测数据，GNSS 坐标时间序列已成为地壳形变监测、全球板块运动动力学研究、区域构造运动、断层滑动和地震形变监测等地学研究的主要方法之一[1]。

近年来，中国大陆西部区域的地质构造活动和地表形变呈现异常活跃的态势，每年都会发生各种强烈的地质灾害活动。 因此，有必要对该区域地表形变进行研究，深入揭示中国西部地区块体运动的趋势特征和成因[2]。

已有学者对区域地表形变进行了相关研究，张风霜等利用中国陆态网基准站多年累积的观测数据，对中国大陆GNSS 基线时间序列进行分析，认为中国大陆的GNSS 基线总体上呈现NE 向缩短，NW 向伸长的运动趋势； GPS 基线时间序列可以捕捉强震孕震形变信息，同时说明了大范围区域构造活动的剧烈变化对地震的孕育发生有一定的推动作用[3]。

曹里等提出利用多维度应变参数分析中国川滇地区的应变特征，基于去除欧亚板块背景场的精确速度场，利用Delaunay 三角形构网的方法求得川滇地区的各应变参数。 验证了川滇地区地震频繁发生在最大剪切应变率出现高值和面膨胀率梯度较大区域[4]。

李延兴等提出了板内块体的整体旋转与均匀应变模型，并基于中国大陆地区速度场，计算了大陆8 个块体的应变参数。 结果表明，块体应变参数参考框架无关，计算的应变参数是唯一的[5]。

选用中国西部陆态网80 个GNSS 连续观测站2010 年6 月28 日至2018 年8 月14 日的观测资料，获得各个测站的GNSS 坐标时间序列结果，基于上述结果计算研究区域的速度场、应变速率场和基线时间序列，进而研究分析监测区域的地表形变[6-7]。

1 GNSS 数据处理和坐标时间序列分析

1.1 GNSS 数据收集

选用中国西部陆态网80 个GNSS 连续观测站2010 年6 月28 日至2018 年8 月14 日的观测资料，并对数据进行精密处理，从而获得中国大陆西部各地块的构造形变特征。 图1 为中国西部陆态网GNSS 连续观测站点的分布示意，其中，红色和黑色线条分别表示区域界线、地块界线。

1.2 GNSS 数据处理

基线处理软件采用应用较为广泛的GAMIT 10.7 高精度数据处理软件，网平差采用MIT 和SIO 共同研发的GLOBK 软件，该软件采用平差方法是最小二乘估计和卡尔曼滤波法，将全球IGS-H 文件合并进行联合平差，获得这些站点在ITRF2014 框架下的坐标。 所选观测站2010 年6 月28 日至2018 年8 月14 日单天解的NRMS 值的统计结果见图2。 可以看出GNSS 观测网的质量良好，基线解算质量较高。

图1 GNSS 连续观测站站点分布

图2 单天解NRMS 值统计结果

以单周解为例，在解算的历元范围内根据统计，N方向的平均精度为3.5 mm，有60%的测站精度优于2.3 mm，95%的测站精度优于5.2 mm；E 方向平均精度为4.7 mm，有60%的测站精度优于3.0 mm，95%的测站精度优于7.4 mm；U 方向平均精度为8.7 mm，有60%的测站精度优于9.4 mm，95%的测站精度优于13.7 mm。

1.3 GNSS 坐标时间序列分析

GNSS 坐标时间序列不仅携带了代表构造运动的线性信号，也包含代表周期运动、各种因素引起的突变项、噪声信号等非线性信号。 为更加准确提取和分析GNSS 坐标时间序列中的各种信息，对GNSS 观测站的单个坐标分量的坐标时间序列进行建模，可采用如式(1)所描述的非线性模型进行分析。

式中，ti表示GNSS 观测站的单日历元；a、b 表示时间序列的初始位置和运动速度；c、d 和e、f 表示连续观测站周年项和半周年项运动系数；H(∗)表示坐标时间序列中的阶跃函数；gi表示连续观测站的坐标分量在Tp时刻的突变值；kj表示连续观测站在地震发生后的Tpost时刻的地表形变指数函数的振幅大小；τj表示GNSS 观测站在地震发生后的松弛时间常数；vi表示坐标时间序列的坐标残差。

按照上述非线性模型，将中国西部陆态网80 个GNSS 连续观测站的坐标时间序列进行建模，并对坐标时间序列进行共模误差剔除和噪声分析，分别在只考虑高斯白噪声和“白噪声+闪烁噪声”噪声类型下求得的坐标时间序列拟合参数的速度项和振幅的大小和误差[8-9]，图3 是测站SCTQ 坐标残差时间序列示意，表1 和表2 分别为速度和振幅对比结果。

图3 SCTQ 坐标残差时间序列示意

表1 速度对比结果 mm/a

续表1

表2 振幅对比结果 mm/a

从表1、表2 的对比结果可知:各坐标分量速度大小相差较小，N、E、U 方向最大相差0.21，0.42，0.29 mm/a；但对应的中误差相差较大，N、E、U 方向对应拟合不确定度差值分别为8.11，7.61，11.46 倍。 时间序列各坐标分量的振幅大小相差不大，N、E、U 方向最大相差0.15，0.18，0.25 mm/a；但对应拟合不确定度相差较大，N、E、U 方向对应拟合不确定度差值分别为3.64，3.69，3.97 倍。 对GNSS 坐标时间序列进行共模误差和噪声剔除，获得精确可靠的坐标时间序列结果。

2 GNSS 坐标时间序列在地表形变监测中的应用

在地表形变问题的研究中，通常从研究区域的速度场、基线时间序列和应变特征参数三个方面入手，分析由于地震、板块挤压碰撞引起的地表形变。 速度场能够反映研究区域块体的总体运动趋势，基线时间序列的变化特性反映了块体之间的碰撞挤压特征，应变特征参数反映了块体应变状态[10-13]，文章从上述三个方面展开研究。

2.1 水平速度场的获取

使用中国GNSS 西部陆态网80 个观测站，时间跨度8.13 年的观测数据。 在进行基线处理和网平差时，选用中国周边12 个IGS 站和稳定的陆态网连续站进行约束，最终得到陆态网连续观测站在ITRF2014 框架下的速度场，并绘制了中国西部陆态网GNSS 连续观测站速度场示意图4，蓝色箭头表示测站速度。

从图4 可以看出，中国西部地块总体上是往东南方向运动，且从南向北速率逐渐减小，从拉萨地块向北依次是羌塘地块、巴颜喀拉地块、塔里木地块、天山地块和准噶尔地块，其块体最大运动速率依次为50.39，45.91，37.91，32.11，30.83，28.39 mm/a。

图4 GNSS 连续观测站速度场示意

2.2 GNSS 基线时间序列变化分析

分析中国西部地区104 条GNSS 基线的时间序列结果的趋势性运动特征，基线分布如图5 所示，为凸显GNSS 连续观测站之间的基线运动趋势，采用最小二乘法对研究区域内的基线时间序列进行线性拟合，最后得到GNSS 连续观测站的运动速率，计算结果如表3 所示(仅列举其中20 条基线成果)。 用基线时间序列的拟合误差作为单位权中误差，用来反映拟合精度，线性拟合得到的速率可以映GNSS 基线的运动趋势特征。

图5 基线分布

表3 中国西部GNSS 基线运动速率

由表3 可知，在欧亚板块和印度洋板块碰撞和推挤力的作用下，中国大陆西部的GNSS 基线时间序列呈现出NE-SW、NW-SE 方向缩短的运动特征。 拉萨地块、羌塘地块也有NE-SW 缩短和NW-SE 方向伸长的运动趋势，缩短最快的XJRQ-XZZB 基线速率为-27.78 mm/a，伸长最快的XZRT-QHTT 基线速率为23.64 mm/a。

2.3 应变时间序列分析

应变场直接反映了局部应力应变的空间分布特征，从地壳内部所受的构造力和应变去分析，通过构造应变模型计算研究区域块体的应变特征参数，根据所求的应变参数特征评估地震灾害。 因此，寻找能够反映研究区域块体内部形变特征的参数和表述方法是确定区域地壳形变特征的基础。

(1) 体旋转与均匀应变模型

计算块体应变速率的方法有:多面函数法、最小二乘配置法、球谐函数法和Delaunay 三角形法。 武艳强等在研究中对比了几种方法的应用可靠性，认为各个方法都有各自应用的优点和局限性，与研究区域的空间尺度和观测数据密度有密切关系[14-17]。 大量研究表明，板块或者地块并不是纯刚体，而更接近于弹性体或粘弹性体。 试验研究和观测结果都证明岩石圈板块和板内地块不是纯刚性的而是弹塑性的，故采用文献[5]提出的块体刚性弹塑性运动应变模型，有

式中，r 为块体上矢径；Ω(kx，ky，kz)为块体的绝对欧拉矢量；λ 和h 分别表示点位的经度和纬度。

(2)应变速率场分析

结合中国西部地区的GNSS 连续观测站点，根据Delaunay 三角形构网算法，以各个块体为研究目标，构造了由68 个测站生成的39 个多边形，图6 中红色细线为中国西部地区陆态网连续观测站点Delaunay 三角形构形结果。 基于上节计算的速度场，求解研究区域的每个多边形的应变参数，绘制中国西部地区应变速率场，分析各个块体的运动趋势和产生形变的内在因素。

图6 中国西部陆态网GNSS 观测站构形结果

根据西部地区GNSS 连续观测站的水平运动速度场，利用模型计算中国西部地区各个块体的运动参数和应变参数，结果如表4 所示。 最后求得中国西部地区陆态网2010～2018 年时间跨度下的主应变速率场，如图7 所示。

图7 主应变率结果

表4 中国西部地区主应变参数结果

续表4

从表4 和图7 可以看出，天山块体和准噶尔块体区域主压应变率的方向为北西向并呈现顺时针右旋的趋势，主压应变率最大为-35.9×10-9·a-1；而主张应变率较小，最小为0.5×10-9·a-1；面膨胀率均为负值，最大为-34.4×10-9·a-1，在欧亚板块北西方向的强大推力作用下，天山块体和准噶尔块体呈现北西-南东压缩的运动趋势。 塔里木地块的主压应变方向是北西130°，该地块主压和主张应变都比较小，相对较为稳定。 祁连地块主压应变速率方向为北西方向并呈现顺时针右旋的趋势，主压应变率最大为-37.4×10-9·a-1，而主张应变率则要小于主压应变率，最小为2.6×10-9·a-1，面膨胀率都为负值，最大为-16.8×10-9·a-1，该地块整体处于压缩状态。 横断山脉核心区域应变速率场呈现弧状向东南方向发散趋势明显，且川滇地块北部地区，东西向拉伸明显，南北向挤压也比较显著，主压应变率最大为-33.7×10-9·a-1，主张应变率较最大为41.3×10-9·a-1，此处位于四川盆地和云贵高原的交界处，受到东西部两侧挤压造成该地区拉伸和收缩都比较强烈。

3 总结与展望

以中国西部陆态网GNSS 观测站积累的观测数据对地表形变进行监测作为研究主题，对GNSS 坐标时间序列进行共模误差和噪声剔除，获得精确可靠的坐标时间序列结果。 基于坐标时间序列结果总结了分析地表形变的常用方法，通过求得GNSS 连续观测站的二维平面速率场，基线时间序列和应变速率场来获取块体的运动趋势和地表形变规律。

本研究仅利用了整体旋转与均匀应变模型按照Delaunay 三角形构网规则计算了研究区域的应变参数，在今后的研究中，应对主流的最小二乘配置法、球谐函数法和多面函数法进行应参数求解进行深入研究和对比。