随着GPS观测技术的发展和定位精度的不断提高, GPS技术在测绘等领域中的应用越来越广泛^{[1, 2]}。但GPS时间序列数据处理工作非常繁琐、耗时，且很多重复性工作，而实际数据处理过程中缺少相应的功能较完备的GPS时间序列分析软件。本文在分析GPS时间序列数据处理及分析基本流程的基础上，开发了基于MATLAB的GPS时间序列数据处理与分析软件。

1 软件的设计

高精度地球动力学研究往往以大尺度、长周期、等间距采样的台站GPS时间序列为研究对象，其台站GPS时间序列数据处理与分析是一个复杂的过程，主要包括GPS单日解坐标序列解算、构造运动分析、负载效应改正及分析、GPS时间序列时空滤波、坐标序列噪声模型等5个部分，其数据处理与分析基本流程如图 1所示。

1.1 软件总体结构设计

软件主要包含GPS软件安装及配置、数据下载及归档整理、数据预处理、空间分析、坐标时间转换、时频分析、图表绘制、空间滤波分析等7大模块，各个模块包含若干下拉菜单，实现不同的功能，各模块之间既相互独立，又存在联系，即某一模块的输出结果为另一模块的输入数据，通过模块之间的集成与融合，形成一个完备的、可靠的、易操作的GPS时间序列分析软件。

1.2 软件功能模块设计

1) 软件安装配置及数据下载模块。软件安装及配置模块是本软件的基础功能之一，该模块主要介绍国际上主流高精度GPS后处理软件GAMIT、BERNESE、GIPSY及QOCA的安装及配置方法，并给出了一些常见问题解决方案^[3]。数据下载及归档整理模块提供了GPS时间序列涉及的原始数据下载接口，给出了GPS相关数据的链接入口。此外，考虑到GPS数据处理及时间序列要求数据文件严格按照一定格式进行存储，模块支持数据文件批量整理。

2) 数据预处理模块。数据预处理模块主要是对GPS数据(原始观测数据、单日解序列)进行必要的处理，主要包括基于TEQC的GPS数据质量分析(格式转换、数据编辑、质量检查、单点定位、绘制天空图等)。另一个功能就是对单日解序列进行插值拟合，对丢失的数据进行插值^[4]。

3) 空间分析模块。空间分析模块是软件的核心模块之一^[5]。该模块通过空间分析模块，为共模误差的分离提供一些依据，模块功能主要包含站点坐标时间序列的相关项分析(自相关、互相关)、站点之间距离计算、大尺度GPS网尺度划分、阶跃探测及改正、负载效应分析等，通过对多个因子进行分析，为后续的共模误差滤波分析提供依据。

4) 空间滤波模块。目前常用的滤波方法有基于主分量分析(principal component analysis, PCA)的空间滤波方法、基于频率域的独立成分分析方法和卡尔曼滤波等降维分析方法，通过空间滤波模块可以采用这3种方法进行独立分析^{[6, 7]}。

5) 时频分析模块。时频分析模块提供了功率谱分析(FFT (fast Fourier transform)和Lomb算法)和极大似然估计两大子菜单。在功率谱分析中，傅里叶变换(Fourier transform)是常用的一个方法，是信号在时间域和频率域转换的桥梁^{[8, 9]}。

6) 绘图模块。绘图模块具备了强大的时间序列相关图形绘制，主要功能包括坐标时间序列绘制(N、E、U)、站点分布图绘制、PCA空间响应绘制、阶跃及其影响结果绘制、梯度色阶图绘制(如B、L、WRMS绘制)、站点相关性绘制等。

7) 坐标时间转换模块。坐标时间转换模块的主要功能是对不同坐标时间系统进行转换，满足实际需求。

2 软件的实现 2.1 空间分析模块功能的实现

空间分析模块主要包含站点坐标时间序列的相关项分析、站点之间距离计算、大尺度GPS网尺度划分、阶跃探测及改正、负载效应分析。

相关性分析是指对两个或多个具备相关性的变量元素进行分析，从而衡量两个变量因素的相关密切程度。相关性系数的公式为^[10]：

$ r = \frac{{\sum {XY} - \frac{{\sum X \sum Y }}{N}}}{{\sqrt {\left( {\sum {{X^2}} - \frac{{{{\left( {\sum X } \right)}^2}}}{N}} \right)\left( {\sum {{Y^2}} - \frac{{{{\left( {\sum Y } \right)}^2}}}{N}} \right)} }} $

(1)

利用MATLAB求相关系数的Corr函数，可以快速求得某两个时间序列的相关系数，然后用Polyfit、Polyval这两个函数可以对画出的相关系数散点图进行拟合分析，得到这些散点的趋势线，在此基础上可以对其进行分析^[11]。点击主界面菜单的空间分析按钮，继而会弹出下拉菜单，单击相关性分析，会弹出选择文件对话框，选择输入文件数据(坐标时间序列)后，弹出如图 2所示的结果。为了更好地表现其相关性的趋势，模块还能对单一站点的自相关性进行分析。

区域站点分析主要包括站点之间的距离计算、大尺度GPS子网划分。站点距离计算根据站点的(B、L、H)或者(X、Y、Z)进行计算。大尺度GPS网尺度划分主要是用于GPS子网的划分，软件提供两种子网划分方法：① 基于经纬度的矩形区域划分方法，此方法是基于MATLAB GUI上实现的，读入Excel数据文件和输出Excel数据文件是通过Pushbutton控件来完成的，而显示是通过Uitable控件来实现的，数据保存在当前文件目录下。点击输入数据可读取全球的IGS站点坐标，结果显示在图 3的表格上，输入需要选取区域的经度、纬度的极值，点击结果就可将符合选取区域的IGS站的坐标提取出来，进而显示在右列的表格中，点击保存即可将提取后的区域站点保存为当前路径下的Excel表格文件。② 基于中心点及半径的圆形区域划分法，用户只需输入站点名或者经纬度及半径，即可选取所选区域内的所有IGS站。

阶跃探测是空间分析的另一个功能，GPS台站受到周围环境变化、仪器及天线变更、强震影响等，会使得GPS坐标序列产生瞬时突变，称之为阶跃。阶跃的存在使得坐标序列和后续的空间滤波存在偏差，本模块能有效地探测、分析阶跃，并对其进行改正及对结果进行绘制。点击主界面菜单的空间分析，弹出下拉菜单，单击阶跃绘制，结果如图 4所示。

2.2 时频分析功能的实现

傅里叶变换广泛应用于信号分析、光谱和声谱分析、全息技术等各个领域。MATLAB提供了一套计算快速傅里叶变换的函数，它们包括一维、二维和N维离散傅里叶变换函数fft、fft₂和fft_n。在读入站点坐标残差文件后，使用fft函数即可将时间序列转换到频率域上，进而揭示GPS信号的潜在规律^{[8, 12]}。单击主界面下的“时频分析”按钮，可进行快速傅里叶分析，并且绘制图像，经快速傅里叶变换之后的GPS信号频谱图如图 5所示。

2.3 图表绘制功能的实现

利用MATLAB强大的画图功能，使用errorbar(Time, Average, Variance)函数来实现误差棒的绘制。时间序列绘制主要用于显示数据处理结果并绘制时间序列，用户只需按照标准格式提取时间序列残差文件，其格式为“Year, X, Y, Z，σ_X, σ_Y, σ_Z”。单击“时间序列绘制”按钮，就可显示数据处理结果并绘制时间序列的误差棒，结果如图 6所示。

3 结束语

本文开发的基于MATLAB GUI的GPS时间序列分析软件，运行平台是Windows平台，通过不同模块集成，涵盖了GPS时间序列数据处理及应用分析等功能，形成了完善的GPS时间序列分析软件。运行结果表明，除了具备一般软件的界面友好的特点外，还能够方便地进行各种数据处理及图像分析，提高GPS时间序列数据处理及分析的效率，为广大GPS数据处理及时间序列分析领域人员提供了一个有力的工具。