GNSS测站速度估计模型重建

徐克科,伍吉仓

(1.河南理工大学 测绘与国土信息工程学院，河南 焦作 454000;2.同济大学 测绘与地理信息学院，上海 200092)

GNSS测站速度估计模型重建

徐克科1,伍吉仓2

(1.河南理工大学 测绘与国土信息工程学院，河南 焦作 454000;2.同济大学 测绘与地理信息学院，上海 200092)

获取准确的GNSS测站速度对于研究全球板块运动、地壳形变、地震活动及其地球动力学过程至关重要。为此，以GNSS基线解算后的基线向量作为观测值，重建了最小二乘综合速度解算模型和卡尔曼滤波速度估计模型。模型中，考虑了测站坐标、运动速度、年、半年周期项一同作为参数，在基线解网平差的同时，一并求解获取速度估值。同时，利用坐标时序分析的方法，顾及白噪声和幂律噪声的影响，对单日解坐标时间序列重建了时序速度拟合估计模型，以获取长期趋势项作为速度值。基于3种模型，对川滇地区中国地壳运动观测网络2010—2014年21个GPS基准站的观测数据进行了速度估计与对比分析。结果表明：3种模型所估计测站速度非常接近，差异基本处于0～1 mm/a范围之内；速度估值中误差均在亚毫米水平。由此得出，3种速度估计模型具有本质的一致性，均可正确估计测站运动速度，能够满足高精度地壳形变研究的需要。

GNSS；速度；基线向量；卡尔曼滤波；时序分析

要正确理解全球板块运动、地壳形变构造运动过程及其动力学机制，首先必须先获取高精度、可靠的测站运动速度[1-6]。如果只观测两期，通过GNSS高精度处理软件经赫尔模特坐标变换可得到这两期的ITRF框架下的坐标，两期的坐标相减除以两期间隔时间便可得到测站实际运动速度。然而，地壳形变的研究通常需要经过长时间连续站连续观测或区域站多期的持续观测。随着GNSS观测资料的持续积累，通过GNSS长期的观测数据解算的单日解坐标时间序列可以估计出更高精度的测站运动速度。目前，GNSS测站速度估计的方法主要可归结于两大类。一是将测站坐标和速度作为参数，建立关于基线向量的观测方程，在对基线解网平差的同时一起求解[7-8]。二是对网平差后单日解测站坐标时间序列进行时序分析，拟合出线性趋势项作为速度[9-12]。有关研究表明，当GNSS坐标时间序列跨度大于4.5 a时，拟合中是否有季节项的加入对速度估计误差影响可以忽略不计；低于2.5 a，对速度估计影响较大，不可接受；当时间跨度是整数年加上半年时，季节项对速度场的影响最小[13-14]。GAMIT软件的输出结果文件H文件提供了单日解基线向量参数估计结果及其协方差阵[9，15]，这为利用合并单天解估计速度提供了有利条件。对此，以基线向量作为观测值，引入测站坐标、运动速度、年、半年周期项一同作为参数，重建了基线向量最小二乘综合解算模型、基线向量卡尔曼滤波解算模型；顾及白噪声和幂律噪声的影响构建了坐标时序分析速度拟合估计模型。利用川滇地区中国地壳运动观测网络2010—2014年GPS基准站观测数据，分别基于这3种模型进行了GNSS测站运动速度的估计和对比分析。

1 基线向量最小二乘综合解算模型

(1)

(2)

(3)

上述的观测方程中，待求参数为测站i和j在参考时刻m的坐标、测站i，j的运动速度和年、半年周期项系数，观测值为基线三维向量，权阵为GNSS基线解算的协方差阵的逆阵。这只是基于其中一条基线所列出的误差方程。同理，可列出所有时段整个GNSS网所有基线的观测方程并联立求解。因基线解算结果单天解是松弛解，联合解的站坐标组成的网没有固定的参考基准，因此需要确定参考框架。采用若干稳定的IGS站的站坐标和速度进行拟稳平差或重心基准平差。

2 基线向量卡尔曼滤波估计模型

利用GAMIT软件的输出结果的H文件提供的单日解基线向量参数估计结果及其协方差阵，基于动态卡尔曼滤波理论，将测站坐标、测站速度、年、半年周期项作为具有方差信息的准观测值，建立了动态卡尔曼滤波漂移速度估计的状态模型。对一系列的单日解进行卡尔曼滤波合并。并通过对全球IGS 核心站的约束来进行参考框架的定义和七参数转换，从而获取了ITRF参考框架下的测站的坐标和运动速度[9，15]。因测站坐标参数是随机变化的，采用马尔科夫过程来描述这些参数是随机漫步过程。将测站运动速度、年、半年周期项系数作为恒定参数，测站坐标作为时变参数，列出关于基线向量的卡尔曼滤波模型的观测方程和状态方程。观测方程同第1节基线向量最小二乘综合解算模型中的观测方程，见式(3)。建立以测站坐标、速度、年、半年周期项系数作为状态向量的状态转移方程，见式(4)。

(4)

Dt=

(5)

从k到k+1时刻的预报公式为

(6)

(7)

从k到k+1时刻的修正公式为

(8)

(9)

(10)

式中，K为卡尔曼滤波增益矩阵。

上述卡尔曼滤波估计的时间顺序是逐渐增加的，直到最后参考历元时刻的参数最佳估计值和协方差矩阵，这是前向卡尔曼滤波。还需要再运行后向卡尔曼滤波估计进行回代。后向卡尔曼滤波估计与前向卡尔曼滤波估计公式相同，只是时间顺序相反。

3 坐标时序拟合估计模型

对于长时间连续GNSS坐标时间序列，为获取GNSS测站精确的速度场信息，可以对此采用时序分析的方法进行速度拟合。GNSS单站、单分量坐标时间序列速度拟合函数模型可表示为

y(ti)=a+bti+csin(2πti)+dcos(2πti)+

esin(4πti)+fcos(4πti)+

(11)

式中:ti为以年为单位的GPS站点单日历元，H表示Heaviside阶梯函数，参数a表示站点的起始位置，b表示线速度，c和d组合表示年周期项系数，e和f表示半年周期项系数，第七项表示发生在历元Tgj处大小为gj的ng个偏移常量，即站点由于地震或者更换接收机天线等原因产生的后续历元位置整体偏移。

在速度估计时，考虑到有色噪声对速度估计的影响，通过构造最大似然函数估计谱指数。对于幂律谱噪声，协因数阵Jpl可表示成转换矩阵与其转置的乘积[11,16]，即

(12)

转换矩阵T(a)可表示为

(13)

(14)

式中:σwh，σpl为白噪声和有色噪声分量；Rwh，Jpl为白噪声和有色噪声协方差矩阵；Rwh为单位阵。构造噪声分量和谱指数的极大似然函数为

(15)

式中，n表示时间序列的总长度，其对数形式为

(16)

采用单纯形法求解上式中的最大值，得到相应的似然值及白噪声与有色噪声分量σwh，σpl。研究结果表明，原始时间序列的最优噪声模型与PL+VW噪声模型的极大似然值相差甚微，幂律噪声的水平可以反映时间序列的整体噪声水平[17]。因为闪烁噪声和随机漫步噪声均属于幂律噪声类型，在此，采用PL+VW噪声模型进行了速度拟合。

4 实例分析

利用川滇地区中国地壳运动观测网络21个基准站(包括框架站)自2010—2014年的观测数据，经GAMIT10.4软件处理后分别得到了单日基线解和单日坐标解，部分站的结果见图1和图2。因不同软件，不同方法的解算结果之间的比较，更能客观地评价数据处理结果的外部一致性与可靠性。为此，同时利用BERNESE5.2软件也进行了解算，作了比较。由图1基线结果看，两种软件解算的DX、基线长度DL离散度均在5 mm之内，DY，DZ离散度在10 mm之内；两者基线时间序列整体吻合度较高，变化趋势较一致。由图2坐标解算结果看，两种软件网平差后单日解坐标时间序列离散度接近，基本都在10 mm之内，中误差均优于3 mm；两者坐标时序整体趋势变化一致。可见，数据处理过程中所解算基线向量和坐标结果可靠，为后续测站运动速度估计的准确性提供了保障。

图1 基线SCLH-SCXJ、SCDF-SCYX三分量与基线长度单日解时间序列

图2 测站SCLH、SCXJ、SCDF和SCYX单日解坐标时间序列及中误差

对于基线解，分别利用第1节基线向量最小二乘方法和第2节基线向量卡尔曼滤波方法进行了速度估计。对于单日解坐标时间序列，利用第3节时序分析的方法拟合了测站运动速度。3种方法解算的21个台站的速度结果见表1。由表1所示，3种模型所解算测站速度值非常接近，之差大都小于1 mm/a；速度估值中误差均在亚毫米水平。这表明，3种速度估计模型均可正确解算测站运动速度，能够满足高精度地壳形变研究的需要。最小二乘和卡尔曼滤波估计模型结果十分接近，两种模型具有本质的一致性。时序拟合模型结果表明，所有测站残差的有色噪声功率谱指数基本处于(-1，0)区间内，表现为闪烁噪声。不考虑时间序列中的有色噪声会导致拟合速度精度的严重高估。顾及有色噪声的影响，GPS测站拟合速度值变化不大，但是拟合速度中误差将扩大约2～8倍。同时，对于5 a的GNSS坐标时间序列，是否有季节项的加入对速度估计影响微乎其微，可以忽略。

表1 3种速度估计模型的估计结果 mm/a

5 结 论

以GAMIT解算的基线向量作为观测值，引入测站坐标、运动速度、年、半年周期项一同作为参数，重建了基线向量最小二乘综合解算模型、基线向量卡尔曼滤波解算模型。并顾及白噪声和幂律噪声的影响构建了坐标时序分析速度拟合估计模型。利用川滇地区中国地壳运动观测网络2010—2014年GPS基准站观测数据，分别基于3种模型进行了GNSS测站运动速度估计和对比分析。

1)基线向量最小二乘综合解算模型是利用三维自由网解得到的所有基线向量纳入统一综合模型作为观测值，基线向量的验后方差-协方差阵则被用来确定观测值得权阵，将参考时刻的测站坐标、速度、周年、半周年系数作为未知参数，把所有时段所有基线联合起来利用最小二乘一同求解。同时需附加一定的约束条件，如固定其中的某些点的坐标速度，或进行拟稳平差和重心基准平差等，求出最终解，这通常是法方程迭加过程。随着观测期数的增加，A阵变得很大，对计算机存储空间要求较高，数据处理变得复杂。而基线向量卡尔曼滤波估计勿需保留用过的观测值序列，按照一套递推算法，把参数估计和预报有机地结合起来，计算速度较快，节约内存。坐标时序拟合模型是对GNSS单站N，E，U三分量的坐标时间序列进行时序分析，考虑长期线性运动+年/半年周期运动+阶跃，基于谱指数计算或最小范数二次无偏估计(MINQUE)或极大似然估计，可以顾及白噪声、闪烁噪声、随机游动、幂指数噪声和带通滤波噪声等不同噪声组合，采用加权最小二乘法求取测站长时间的运动速度。同时还可以对不同噪声特性及产生原因进行深入分析。

2)3种模型速度估计结果表明，三者速度值十分接近，差异大都位于0～1 mm/a范围之内，速度值中误差均在亚毫米水平。由此验证了3种速度估计模型的一致性和可靠性，这说明3种模型所估计的测站运动速度均可以满足高精度地壳形变研究的需要。

3)时序分析结果表明，所有测站残差的有色噪声功率谱指数基本处于(-1，0)区间内。可见，5 a的数据主要表现为闪烁噪声。不考虑时间序列中的有色噪声会导致拟合速度精度的严重高估。顾及有色噪声的影响，相比只考虑白噪声，GPS测站拟合速度值变化不大，但是拟合速度中误差将扩大约2～8倍。由于随机漫步噪声的长延时相关性，而川滇地区连续站坐标时间序列的时间过短，仅5 a时间，很可能还不足以准确量化随机漫步噪声的特征。同时，对于5 a的GNSS坐标时间序列，是否有季节项的加入对速度估计影响微乎其微，可以忽略不计。

[1] 赵国强,苏小宁.基于GPS获得的中国大陆现今地壳运动速度场[J]. 地震,2014,34(1):97-103.

[2] 邹镇宇,江在森,武艳强,等. 基于GPS速度场变化结果研究汶川地震前后南北地震带地壳运动动态特征[J]. 地球物理学报,2015,58(5):1597-1609.

[3] 程鹏飞,文汉江,孙罗庆,等. 中国大陆GPS速度场的球面小波模型及多尺度特征分析[J]. 测绘学报,2015,40(10):1063-1070.

[4] 魏子卿,刘光明,吴富梅. 2000中国大地坐标系:中国大陆速度场[J]. 测绘学报,2011,40(4):403-410.

[5] 张风霜,占伟.利用GNSS连续观测资料获取高精度动态速度场的研究[J].地震研究,2015,38(1):75-83.

[6] 王琪,张培震,马宗晋.中国大陆现今构造变形GPS观测数据与速度场[J].地学前缘,2002,9(2):415-429.

[7] 王解先.由GPS基线向量解算地面形变[J].同济大学学报(自然科学版),2005,33(7):967-970.

[8] 王解先.GPS精密定轨定位[M].上海:同济大学出版社,1997.

[9] HERRING T A,KING R W,MCCLUSKY S C. GLOBK Reference Manual:Global Kalman filter VLBI and GPS analysis program[R].Release 10.4.Department of Earth,Atmospheric,and Planetary Sciences,Massachusset Institute of Technology. 2010.

[10] WILLIAMS S D P.CATS:GPS coordinate time series analysis software[J].GPS solutions. 2008,12(2):147-153.

[11] LANGBEIN J.Estimating rate uncertainty with maximum likelihood:differences between power-law and flicker-random-walk models[J].Journal of Geodesy,.2012,86(9):775-783.

[12] WILLIAMS S D P.The effect of coloured noise on the uncertainties of rates estimated from geodetic time series[J].Journal of Geodesy,2003,76(9-10):483-494.

[13] BLEWITT G,LAVALLÉE D.Effect of annual signals on geodetic velocity[J].Journal of Geophysical Research:Solid Earth (1978-2012),2002,107(B7):ETG 9-1-ETG 9-11.

[14] BOS M S,BASTOS L,FERNANDES R M S.The influence of seasonal signals on the estimation of the tectonic motion in short continuous GPS time-series[J].Journal of Geodynamics, 2010,49(3):205-209.

[15] HERRING T A,KING R W,MCCLUSKY S C. GAMIT Reference Manual:GPS Analysis at MIT[R].Release 10.4.Department of Earth,Atmospheric,and Planetary Sciences,Massachusset Institute of Technology，2010.

[16] 姜卫平,周晓慧. 澳大利亚GPS坐标时间序列跨度对噪声模型建立的影响分析[J]. 中国科学(D辑：地球科学),2014,44(11):2461-2478.

[17] 吴伟伟. 华北地区GPS连续站坐标序列特征研究[D].北京：中国地震局地震预测研究所,2014.

[责任编辑：刘文霞]

Velocity estimation model reconstruction for GNSS observation station

XU Keke1,WU Jicang2

(1.School of Surveying and Land Information Engineering,Henan Polytechnic University,Jiaozuo 454000,China；2.College of Surveying and Geo-Information,Tongji University, Shanghai 200092,China)

To acquire reliable velocity of GNSS observation station is very important for the study of global plate movement,crustal deformation,seismic activity and geodynamic process. Therefore, considering the baseline vector after baselines calculation as observation value, introducing the coordinates，speed，annual and semi-annual seasonal items of GNSS stations as parameters together to be estimated, this paper reconstructs the models of least squares estimation and Kalman filter estimation, by which the velocity is determined along with GNSS network adjustment. Meanwhile, considering the effect of the white noise and power law noise, it is same to the velocity fitting estimation model by the coordinate time series using the method of timing series analysis, by which the long-term trend item is estimated and considered as velocity value. Based on three speed estimation models,the GPS data from the Crustal Movement Observation Network of China from 2010 to 2014 in Sichuan-Yunnan area is processed and analyzed.The result shows that the velocity difference between different models is less than 1mm/a and the root mean square error is at a submillimeter level.The three speed estimation models are veritied to be consistent and reliable,and the estimated velocity can satisfy the need the study of high precision crustal deformation..

GNSS; velocity;baseline vector;Kalman filter;timing series analysis

引用著录：徐克科,伍吉仓.GNSS测站速度估计模型重建[J].测绘工程，2017，26(4)：6-11.

10.19349/j.cnki.issn1006-7949.2017.04.002

2016-02-17

国家重点973基础研究发展计划资助项目(2013CB733304);国家自然科学基金资助项目(41404023)

徐克科(1979-),男，副教授，博士.

P228.4

A

1006-7949(2017)04-0006-06