基于MATLAB-GUI的GNSS基准站网坐标序列数据处理软件 | ![]() |
全球导航卫星系统(global navigation satellite system,GNSS)基准站坐标变化序列反映了在地球物理负荷等各种因素作用下GNSS基准站位置的变化特性,对于地壳形变、陆地水等各类负荷变化的反演分析有重要意义。随着全球和区域高精度GNSS基准站观测网络的建立,观测精度不断提高,GNSS基准站网坐标序列已经成为研究地壳形变、地表质量迁移变化及其对全球气候变化相应的不可或缺的数据。目前,解算GNSS基准站网坐标主要采用Bernese、GAMIT/GLOBK以及PANDA等软件,但处理分析基准站网的坐标序列缺乏功能完善且界面友好的软件。国际上最早由Herring[1]开发GGMATLAB的Tsview模块,实现了GNSS基准站坐标序列的可视化及偏移分析;Williams[2]开发的CATS软件增加了频谱分析以及噪声分析功能;Goudarzi等[3]开发的GITSA增加了趋势和回归分析功能。目前,国内有关GNSS坐标时间序列分析软件的研制较为缺乏,舒颖等[4]基于MATLAB平台开发了GNSS时间序列数据处理与分析软件。上述软件不但功能单一,而且只能处理单一格式的基准站数据。此外,GGMATLAB和GITSA软件安装过程复杂且严重依赖MATLAB函数包,而CATS软件操作繁琐且易受操作系统的限制。因此,本文介绍了项目组开发的一个功能更强且界面友好的GNSS基准站网坐标序列数据处理和分析软件MTSA。该软件可实现区域GNSS基准站网的坐标序列数据自动处理,软件界面良好, 操作简单,具有一定的实用价值。
1 软件设计 1.1 设计工具简介MTSA基于MATLAB开发环境下的脚本语言和图形用户界面(graphical user interface,GUI)进行开发。GUI是由窗口、光标、按键、菜单等对象构成的一个用户界面,有别于传统的MFC(C++)、WPF(C#)或者wxPython(Python)等界面框架,GUI除了支持一般的Windows基本控件,同时提供MATLAB数学库的接口以及高质量的图形输出。
1.2 软件界面设计在MATLAB控制台上输入guide命令即可调出GUI设计面板,具体步骤如下:①构建figure窗口作为主窗口并设置相关属性,后加入的子对象可以依据主窗口的句柄值添加到其中; ②根据软件功能需求向主窗口中添加菜单栏、工具栏、按钮、标签、坐标轴和状态栏等控件,构建子窗口作为响应窗口,并编写callback函数来实现相应的功能; ③使用GUI自带的布局器对界面进行布局,通过样式表等手段美化界面。
在设计GNSS基准站网坐标序列数据处理软件界面时,主要通过Tab页构建界面总体框架,软件的总体框架如图 1所示,主要包括图形用户界面控制、输入输出、数据预处理、时空滤波、时频分析和实用工具6个模块。图 1中,缩写的含义为PCA(principle component analysis)和KLE(karhunen-loeve expansion)。
![]() |
图 1 软件总体框架 Fig.1 Overall Framework of Software |
2 软件功能介绍
开发的软件的主功能界面如图 2所示,主菜单包含文件、工具、关于3项,各项包含若干子菜单,工具栏提供了一些便捷功能。其中,
![]() |
图 2 软件主界面 Fig.2 Main Interface of Software |
1) 输入输出模块。该模块用于软件内部模块与模块间的数据传送同时也是软件与外部进行数据交换的接口。软件能够导入原始的坐标序列数据,在设计时已将该模块映射到其余功能模块中并且内置了图形输出模块,用户可以有选择地输出各个模块的分析后数据。
2) 图形用户界面控制模块。该模块用于软件与用户进行交互操作,包含菜单栏、工具栏、主功能界面和状态栏。菜单栏、工具栏前面已有介绍,MTSA的主功能界面如图 2所示,软件左上为一个项列表,成功导入数据后,列表中会显示所有测站名,可以自由地选择测站进行分析。软件左下为一个文本框,用于记录每个操作的返回信息,软件的中间区域为图形显示窗口,用于输出分析结果。软件的右侧为一系列按钮、标签、单选框等控件的组合,构成了软件的主功能菜单。
3) 数据预处理模块。由于测站的迁移、天线相位中心的变化以及各种外界环境影响,导致观测序列不完整且存在异常值。MTSA采用两种方法探测异常值:①阈值法,根据用户给定的限值(如3倍中误差),超出限值之外的当作异常值剔除; ②IQR(inter quartile range)法则, 基于统计学的稳健孤立值探测法,根据统计指标构造构造相应的探测区间。软件可以单独选择某一个法则或者结合两个法则进行异常值的探测,同时提供了3种插值方法。
已有研究表明,在对区域GNSS网站点坐标序列进行分析时,不同站点坐标时间序列中包含一种与时空相关联的误差, 即所谓的共模误差[5],分析各站点形变特征时应去除。MTSA提供了Stacking、PCA和KLE 3种滤波方法提取共模误差,选择滤波方法后点击提取,数据预处理分页上会显示所提取的共模误差。
4) 时频分析模块。本模块为软件功能的核心模块,MTSA提供了诸如谱分析、小波多分辨率分析、经验模态分析等基本分析功能。
在经典谱分析中,常采用快速傅里叶变换将原本难以处理的时域信号转换成易于分析的频域信号,众多研究表明,GNSS坐标序列数据除含有白噪声外,还包含有色噪声[6-8]。可以通过对GNSS功率谱拟合得到的谱指数确定噪声类型,可通过方差分量估计定量求出噪声分量。软件中选好谱分析类型后点击分析按钮,则可得当前测站的3个方向的功率谱周期图和拟合图,并显示谱指数和噪声分量值。在多分辨率分析模块中,MTSA提供了小波多分辨率分析和经验模态分析。与谱分析不同的是,多分辨率分析能够同时在时域和频域上识别时间序列信号的特性。小波分析法通过选定基函数,使用Mallat算法将任何一个平方可积空间中的信号按照指定分辨率分解成一系列正交小波子空间的分信号,根据子信号的低频部分和高频部分完全重构。软件右下方提供了多分辨率分析的模块,只要选择好小波基函数和分解层数,点击小波按钮,则在中间多分辨率分析分页上自动绘出分解图。与小波分析不同,经验模态分解无需设定基函数,就可以根据时间序列本身的时间特性自适应地将时间序列信号分解成若干本征模函数,每个函数分量对应着不同的频率特征[9]。
5) 实用工具模块。本模块提供了GNSS坐标序列格式转换和GNSS时间格式转换两个实用工具。与统一使用RINEX格式的GNSS数据不同,不同分析中心提供的GNSS坐标序列格式有多种,MTSA目前能够支持以下3种格式:美国板块边界观测实验室提供的pos格式;美国宇航局喷气推进实验室提供的tseries格式;CATS软件中使用的neu格式。
进行格式转换时,软件能够自动识别原始数据格式,只需选择转换模式以及目标格式,软件会在指定位置生成转换文件,该模块同时支持公历日、儒略日、GPS(global positioning system)周秒的相互换算。
3 算例分析本文使用GAMIT对苏州CORS(continuously operating reference stations)网进行基线解算,借助GLOBK提供的计算时间序列的子程序glred对基线解算获得的松弛解文件进行平差,得到约束解文件,从而获得苏州CORS网各站点的坐标时间序列,选取CHAS等10个站2009-2013年的坐标序列数据作为算例。
软件成功导入数据后,依次进行异常值探测、断点检测、插值、共模误差分析等数据预处理工作。考虑到GNSS坐标序列的偏态特性,本文采用改进的IQR法则[10]进行异常值的探测,该方法弥补了标准IQR法则要求数据正态分布的不足,引入非正态数据分布偏态因子对GNSS坐标序列的偏态性进行度量,改进标准IQR法则中异常值的探测区间,循环遍历每个数据值,若在探测区间外则剔除。以CHAS站为例,图 3中的红圈为采用改进的IQR法则探测出的异常值。分别使用Stacking、PCA和KLE共3种滤波方法提取10个测站的共模误差,其中CHAS站的3个方向共模误差如图 4所示。
![]() |
图 3 异常值检测 Fig.3 Outlier Detection |
![]() |
图 4 共模误差 Fig.4 Common Mode Error |
由图 4可知,3种方法得到的结果类似。但Stacking法得到的共模误差与PCA和KLE法差别较大,主要原因是Stacking法假定所有站点具有相同的共模误差,其本质是所有站点残差序列的加权平均值,然而这种假定是不合理的,因为不同站点的共模误差分布不可能完全一致,因此该方法有极大的局限性;而PCA法舍弃了共模误差均匀分布的假设,转而让数据本身去揭示共模误差的空间分布,因此得到的共模误差更合理;KLE法和PCA法类似,只是以标准化后的协方差矩阵取代原协方差矩阵,这种变换能够有效地抑制强烈的局部噪声,但同时也会造成信号失真[11]。依次对所有测站各个方向坐标序列进行上述预处理,计算处理前后坐标序列的RMS(root mean square)值,图 5给出了10个测站预处理前后RMS。从图 5中可看出,经预处理后数据质量明显提高。对CHAS站坐标序列进行快速傅里叶变换算得功率谱作为纵轴,以周期作为横轴,即得功率谱密度图,如图 6所示,周期图的峰值对应的横坐标为时间序列可能存在的周期,从图 6中可以看出,CHAS站3个方向坐标分量存在明显的半年和一年周期。图 7为CHAS站小波多分辨率分析结果。
![]() |
图 5 预处理前后RMS Fig.5 RMS Values Before and After Preprocessing |
![]() |
图 6 周期功率谱图 Fig.6 Periodic Power Spectrum |
![]() |
图 7 小波多分辨率分析 Fig.7 Wavelet Multi-resolution Analysis |
4 结束语
随着空间大地测量观测技术的提高,目前已经积累了多年的GNSS基准站坐标时间序列数据,利用这些数据可进行大量地球物理现象研究。近年来,随着GNSS坐标时间序列分析理论算法的蓬勃发展,有必要开发一款高质量的GNSS坐标时间序列自动处理分析软件。本文介绍了一款基于MATLAB及其图形用户接口平台开发的GNSS基准站网坐标序列数据分析处理软件。该软件界面良好,操作简便,功能完善,具有一定的实用性,可用于GNSS坐标时间序列的数据预处理、噪声分析以及信号分析。相比直接使用代码进行处理的方法,图形界面软件操作更为简便、高效、直观,从而为GNSS基准站坐标序列的数据处理和信号提取提供了一个强有力的工具。
[1] |
Herring T. MATLAB Tools for Viewing GPS Velocities and Time Series[J]. GPS Solutions, 2003, 7(3): 194-199. DOI:10.1007/s10291-003-0068-0 |
[2] |
Williams S D P. CATS: GPS Coordinate Time Series Analysis Software[J]. GPS Solutions, 2008, 12(2): 147-153. DOI:10.1007/s10291-007-0086-4 |
[3] |
Goudarzi M A, Cocard M, Santerre R, et al. GPS Interactive Time Series Analysis Software[J]. GPS Solutions, 2013, 17(4): 595-603. DOI:10.1007/s10291-012-0296-2 |
[4] |
舒颖, 花向红, 贺小星, 等. GPS时间序列数据处理与分析软件的设计及实现[J]. 测绘地理信息, 2017, 42(5): 29-32. |
[5] |
周茂盛, 郭金运, 沈毅, 等. 基于多通道奇异谱分析的GNSS坐标时间序列共模误差的提取[J]. 地球物理学报, 2018, 61(11): 4 383-4 395. |
[6] |
Mao A, Harrison C G A, Dixon T H. Noise in GPS Coordinate Time Series[J]. Journal of Geophysical Research Solid Earth, 1999, 104(B2): 2 797-2 816. DOI:10.1029/1998JB900033 |
[7] |
Amiri-Simkooei A R. Noise in Multivariate GPS Position Time-Series[J]. Journal of Geodesy, 2009, 83(2): 175-187. |
[8] |
Amiri-Simkooei A R. Non-negative Least-Squares Variance Component Estimation with Application to GPS Time Series[J]. Journal of Geodesy, 2016, 90(5): 451-466. DOI:10.1007/s00190-016-0886-9 |
[9] |
张恒璟, 程鹏飞. 基于经验模式分解的CORS站高程时间序列分析[J]. 大地测量与地球动力学, 2012, 32(3): 129-134. |
[10] |
杨凯钧, 袁鹏, 秦昌威, 等. 顾及偏态的IQR法则在GPS时间序列异常值探测中的应用[J]. 大地测量与地球动力学, 2015, 35(5): 793-796. |
[11] |
胡守超, 伍吉仓, 孙亚峰. 区域GPS网三种时空滤波方法的比较[J]. 大地测量与地球动力学, 2009, 29(3): 95-99. |