GNSS 时间序列分析与降噪软件的实现

熊常亮，贺小星，孙喜文，马下平

（1. 中国科学院精密测量科学与技术创新研究院，武汉 430077；2. 中国科学院大学，北京 100049;3. 华东交通大学 土木建筑学院，南昌 330013；4. 陕西铁路工程职业技术学院，陕西 渭南 714000；5. 西安科技大学 测绘科学与技术学院，西安 710054）

0 引言

随着全球卫星导航系统（global navigation satellite system, GNSS）逐步完善，全世界范围内广泛建立了GNSS 连续跟踪观测站，GNSS 观测站对国际地球参考框架的建立与维持发挥了至关重要的作用[1-3]。累计的长周期GNSS 时间序列，为地壳形变监测、全球板块运动等地学研究提供了基础数据，而剔除GNSS 时间序列中夹杂的噪声，已成为当前研究的热点[4-7]。针对GNSS 时间序列分析过程繁琐的问题，文献[8]开发了格格马特拉布（GGMatlab）软件，可对GNSS 时间序列进行站速度、阶跃等进行分析与图形绘制，但无法对时间序列进行降噪分析。文献[9]开发的吉特萨（GITSA），可对时间序列进行基于多种算法的频谱分析与回归分析，但其对矩阵实验室（matrix laboratory, MATLAB）外部工具箱具有依赖性。文献[10]基于派森（Python）语言，开发了特萨纳雷泽尔（TSAnalyzer），但不具备基准站网共模误差提取的功能。文献[11]所开发的软件，虽能够快速实现GNSS 时间序列降噪，但在数据处理算法与数据分析功能上较为单一。

考虑到GNSS 时间序列非线性分离及残余噪声影响，结合当前GNSS 时间序列分析软件存在的不足，本文在介绍GNSS 时间序列噪声处理理论基础上，基于MATLAB 平台开发了一款功能完善、界面友好、操作简单的GNSS 时间序列分析与降噪软件。该软件拥有完整的GNSS 时间序列降噪与共模误差分离功能，兼顾数据预处理、统计分析等扩展功能，且无需外部工具箱，为GNSS 时间序列数据处理与分析提供了一个有力工具。

1 GNSS 时间序列降噪分析

1.1 EMD 降噪数学模型

经验模态分解（empirical mode decomposition,EMD）方法是一种自适应的信号处理方法，可用于GNSS 时间序列噪声剔除[12-14]。它将原始信号序列分解为一系列从高频至低频的本征模态函数（intrinsic mode function, IMF）fIMF。EMD 表达式为

式中：x(t) 为原始时间序列；t为历元数；k为IMF分量的序号；m为分解的总IMF 个数；r(t) 为趋势项。

分解完成后，需判定高频噪声fIMF分量与低频信号fIMF分量的分界值K。常用的判断方法为相关系数法则，该方法计算每个fIMF分量与原始时间序列的相关系数，得到对应的相关系数矩阵。矩阵中的第一个极小值对应的fIMF即为分界fIMF，分界fIMF及分界fIMF之前的分量视为高频噪声分量，在信号重构时不再考虑，相关系数计算公式以及信号重构表达式为

式中：xˆ(t) 为降噪后时间序列；kρ为相关系数；N为总历元数。

获得重构序列后，还可依据相关系数式（2）、均方根误差（root mean square error, RMSE）SRMSE及信噪比（signal noise ratio , SNR）σSNR对重构序列进行精度评估，均方根误差sRMSE的计算公式为

信噪比的计算公式为

若使用真值与重构序列计算上述指标，则当相关系数越接近于 1，重构序列与原始序列拟合越好，降噪效果越好；均方根误差越小，重构序列与原始序列偏差越小，降噪效果越好；信噪比越大，噪声占信号的比重越小，降噪效果越好。

1.2 空间滤波数学模型

连续GNSS 时间序列中存在一种空间相关的误差，称之为共模误差（common mode error,CME）。CME 被普遍认为是GNSS 时间时间序列的主要来误差源之一，CME 的存在，极大地影响GNSS 解的精度与可靠性，尤其是对测站速度估计的影响较大。文献[15]于1997 年提出了一种依据全球定位系统（global positioning system, GPS）单日解残差时间序列进行CME 计算的方法，称之为堆栈滤波法。

在堆栈滤波法中，假设基准站网中有S个站，在剔除GNSS 时间序列的均值与趋势项后，得到其残差时间序列，计算CME 的公式为

式中：εk,t为第k个基准站第t个历元的残差时间序列；εSF(t)为堆栈滤波算得的第t个历元的CME 值。

原始残差时间序列减去CME 值，即为滤波后的残差时间序列。通常以滤波前后时间序列的均方根（root mean square, RMS）SRMS值评估时间序列滤波后的稳定性和精确度。SRMS的计算公式为

2 软件模块及框架设计

2.1 基于EMD 的GNSS 时间序列降噪模块设计

软件中首先实现了基于EMD 的GNSS 时间序列降噪模型。此外，针对EMD 方法中仍存在着模态混叠、分界fIMF值判断不准等问题，影响GNSS 时间序列的精确降噪，实现了三种改进的EMD 方法：

1）削弱混叠法。该方法每次仅分离出最高频的fIMF，通过重复利用EMD 方法不断从剩余的高频fIMF分量提取低频分量，对于模态混叠问题可起到一定的抑制作用[16]。

2）周期密度法。该方法基于平均周期与能量密度的乘积为常数的特点，推导出一种能够自适应确定分界fIMF值的算法，对于分界fIMF值的判断可取得精确的结果[17]。

3）复合指标法。该方法综合考虑时间序列的平滑度与SRMSE计算出一张复合指标，用于分界fIMF值的判断，结果具有可靠性[18]。

此外， 软件中还整合了经验模态分解（ensemble empirical mode decomposition, EEMD）方法，并将EEMD 作为核心算法之一[19]。通过上述五种基于EMD 的降噪方法搭建了软件的降噪分析模块。

2.2 GNSS 时间序列共模误差分离模块设计

基于对堆栈滤波法理论的分析易得，该方法所计算的CME 值是GNSS 基准站网内所有站点在该历元的残差时间序列的平均值。因此，当基准站网的规模过大时，空间分布显然不平均，此时利用堆栈滤波法计算CME 值就失去了意义。

基于此，文献[20]将单日解误差引入堆栈滤波法，提出了加权堆栈滤波法；文献[21]将站点间相关性进行量化计算出相关系数再与堆栈滤波法相结合，提出了相关加权叠加滤波法。文献[22]采用主成分分析法（principal component analysis，PCA）与卡尔胡宁-勒夫（Karhunen-Loeve）变换相结合的方法进行CME 提取，有效抑制了部分测站的地区效应。

综上所述，软件中的空间滤波分析模块将上述四种滤波算法全部实现，使得该模块拥有解算各种尺度的GNSS 基准站网的能力。

2.3 软件框架设计

根据前述的GNSS 模块，结合MATLAB 图形用户界面进行软件开发。软件设计包括数据预处理、空间滤波、降噪分析、扩展功能等模块，软件结构见图1。

图1 软件总体结构

为了提高软件的可操作性，各模块单独设计界面，并将相关源程序、输出数据等保存到对应文件夹，既方便软件后期的维护与扩展同时兼顾软件使用时的清晰简洁。通过批处理功能提高数据处理效率。此外，计算过程中的参数、结果与图形都对应有独立的函数，进一步提高了软件的兼容性，便于更新与升级。

2.3.1 数据预处理模块

数据预处理模块主要是对GNSS 时间序列数据进行必要的处理，主要包括GNSS 时间序列的数据质量分析（格式编辑、时间序列可视化、SRMS值计算、站间相关性分析等）。除此之外还包含粗差剔除功能，可对GNSS 时间序列进行粗差探测与剔除。具体功能如下：

1）格式编辑模块。支持将中国大陆构造环境监测网络及斯克里普斯轨道和常驻阵列中心（Scripps orbit and permanent array center, SOPAC）所发布的两种*.neu 文件批量转换为软件中支持的*.mom 文件，也可将起止时间点设为参数、批量截取时间序列为给定年限的文件。

2）粗差剔除模块。包含3Sigma、5Sigma、3IQR以及MAD 四种粗差剔除算法[23]，在完成粗差剔除后，软件将自动绘制图形、标定粗差点位置，并给出粗差剔除率报告。

3）参数计算模块。可计算SRMS值与相关系数，SRMS值反映时间序列的稳定性的重要指标，见式（7）。相关系数则反映不同站点时间序列间相关性的指标，软件中设计有相应模块来计算该指标，其计算表达式为

计算完成后，相关系数矩阵将以文本形式进行保存。

2.3.2 降噪分析模块

降噪分析模块作为软件的核心模块之一，内置的算法为上述的EMD、EEMD 以及最新的三种基于EMD 的改进算法。该模块可实现完整的GNSS时间序列降噪流程，包括数据导入与仿真、基于EMD 的时间序列降噪以及精度评估等。其中数据仿真方面支持周期性信号、时变季节性信号以及多种噪声模型。同时可绘制一系列与降噪分析相关的图形，如EMD 分解后的相关系数值、降噪前后的时间序列对比图、EMD/EEMD 分解后的每个fIMF等。降噪分析模块也将输出计算过程中的重要参数供算法研究。除此之外，拥有基于快速傅里叶变换的频谱分析功能，可计算GNSS 时间序列的功率谱密度，评估频域中信号的主导周期。

2.3.3 空间滤波分析模块

空间滤波分析模块为GNSS 基准站网时间序列空间滤波而设计，拥有完善的滤波分析功能。核心算法为堆栈滤波法、加权堆栈滤波法、相关加权叠加滤波法以及PCA 方法，因此可适用于各种规模的GNSS 基准站网解算。该模块拥有数据导入导出、站点SRMS值计算、CME 提取及图形绘制等功能。可快速准确地实现数据可视化，输出如滤波前后时间序列、CME、SRMS等一系列图形与表格。

2.3.4 扩展功能模块

软件的扩展功能模块用于增强GNSS 时间序列的表达与具体应用，分为站点搜索、阶跃分析以及统计分析三部分功能：

1）站点检索。基于用户自建GNSS 站点库，根据中心站经纬度或点名，搜寻附近站点并输出距离。再者，GNSS 与卫星激光测距等空间观测技术的并置站时间序列被广泛运用[24-25]。基于此，站点搜索模块中设计有GNSS 站与验潮站并置站搜寻功能，可快速输出站点数据库中的并置站。

2）阶跃分析。可对仿真时间序列添加阶跃，也可以依据数据中心发布的参数，对含有阶跃的实测时间序列进行修正，提高GNSS 时间序列的精确度。

3）统计分析。运用箱型图与小提琴图的算法，对GNSS 时间序列进行可视化表达，直观地反映数据的分布情况。

3 软件实现与算例分析

3.1 软件界面

在MATLAB 的图形用户界面面板，分别给每个功能模块建立独立界面，再通过接口与主界面相连，各模块界面如图2 所示。

3.2 算例分析

3.2.1 时间序列降噪分析

选用一组仿真时间序列，在本软件中进行时间序列的降噪分析。降噪方法选用前述的五种方法，以综合比对不同方法的降噪差异。因仿真数据真值已知，可基于三个传统评价指标定量反映各方法的降噪效果以及降噪分析模块的稳定性与可靠性。仿真数据表达式为

式中：1f、f2、f3为三个周期分量；δ为6 dB 的白噪声。采样点设置为2 000 个，各个分量及噪声的图形见图3。

图3 仿真数据各分量及噪声

采用软件中的统计分析功能对数据进行表达。图4 为基于箱型图与小提琴图的统计结果，图中直观地表达了时间序列的典型代表值、异常值及数据整体分布。

图4 仿真信号的统计分析结果

除统计分析方法外，频谱分析也是GNSS 时间序列常用的分析方法之一。软件的频谱分析功能可利用快速傅里叶变换计算仿真信号的功率谱密度，如图5 所示。

图5 仿真信号的功率谱密度

在降噪分析模块中，使用五种降噪方法，对信号作降噪分析。其中，EEMD 的整合数设为300，噪声与信号的标准差之比设定为0.2[19]。在全部计算完成后，从精度评估报告中读取其相关系数、均方根误差及信噪比三个评价指标值，如表1 所示。

表1 降噪后信号与真实信号的各评价指标

因仿真信号真值已知，故表1 中各指标值是由降噪后数据与真实信号计算而得。相关系数越大，则降噪后序列与原始序列拟合越好；sRMSE值越小，则信号偏差越小；σSNR值越大，则降噪效果越好。基于此，对表1 分析可得，五种方法对于序列降噪都取得了优异的效果。其中，效果最好的是周期密度法与复合指标法，二者精确地判定了分界fIMF值，精度指标值最高，其次是EEMD 法、削弱混叠法及EMD 法。

除使用各评价指标进行定量反映降噪效果外，软件绘制相关图形，直观反映各方法的降噪效果，见图6。其中，该实验中周期密度法与复合指标法的结果一致，故仅列周期密度法的时间序列对比图。上述结果验证了降噪分析模块对于GNSS 时间序列降噪分析的准确性。

图6 不同方法降噪前后信号序列对比

3.2.2 基准站共模误差分离

以美国南加州的一个小区域GNSS 基准站网，时间序列跨度为2005—2011 年的数据进行基准站共模误差分离，以测试软件空间滤波模块的可靠性，站点名称及具体经纬度如表2 所示。

表2 站点经纬度

对GNSS 时间序列进行去均值、趋势项及粗差处理得到残差时间序列，之后对缺失历元进行插值补全。最后采用前文所述中的四种误差分离方法计算出其中的共模误差，用原始残差序列减去共模误差得到滤波后的残余时间序列。在软件中选择堆栈滤波、加权堆栈滤波、相关加权叠加滤波以及PCA 四种方法对原始序列残差时间序列进行CME 提取。不同方法进行共模误差提取前后的残差时间序列SRMS值如图7 至图9 所示。

图7 东（E）方向CME 提取前后的 SRMS值

图8 北（N）方向CME 提取前后的 SRMS值

图9 天（U）方向CME 提取前后的 SRMS值

在CME 提取后，8 个站点在E、N、U方向的SRMS值都有明显降低，同在一个数量级，表明时间序列的稳定性得到了提高。特别是经PCA 方法处理后的时间序列，取得了精度更高的结果，而其他方法的效果整体相当。限于篇幅，仅在图10 至图13 列出KTBW 站经四种不同滤波方法处理后的效果图，实线为原始残差时间序列，点线为滤波后残差时间序列，虚线为对应的CME 值。

图10 KTBW 站堆栈滤波前后时间序列及CME 值

图11 KTBW 站加权堆栈滤波前后时间序列及CME 值

图12 KTBW 站相关加权叠加滤波前后时间序列及CME 值

图13 KTBW 站PCA 滤波前后时间序列及CME 值

4 结束语

根据 GNSS 时间序列数据降噪处理相关原理，开发了基于MATLAB 的GNSS 时间序列分析与降噪软件。该软件考虑GNSS 时间序列分析与降噪的功能需求，充分利用MATLAB 在科学计算与数据可视化上的优势，实现了 EMD 与EEMD 算法模型，并实现了三种基于EMD 的改进算法，联合时间序列仿真、输入输出、精度评估等功能一起组成降噪分析模块。另一方面，将堆栈滤波、加权堆栈滤波、相关加权叠加滤波以及PCA 滤波法与数据批处理功能结合，组成空间滤波分析模块。最后，添加各类数据预处理及统计分析功能，共同组成功能完善、结构完整的GNSS 时间序列分析与降噪软件。根据仿真数据与实测数据测试结果表明：1）软件操作简单，具有良好的稳定性；2）所实现的三种改进算法在传统EMD 方法的基础上得出了更精确的GNSS时间序列；3）空间滤波分析模块快速实现GNSS时间序列的CME 提取，以可视化的形式给出数据处理结果，极大地提高了GNSS 时间序列分析的效率。