-
随着GNSS（global navigation satellite system）观测数据的不断丰富和数据处理技术的不断发展，分析测站非线性运动规律可以获得各种地球物理现象和季节性影响，如果对这部分残差规律进行建模拟合从而消除其影响，那么可以进一步提高GPS（global positioning system）站坐标的精度[1]。尽管目前国际地球参考框架（international terrestrial reference frame，ITRF）中基准站的历元坐标和速度场已达到了毫米级，但几乎所有IGS（international GNSS service）站坐标时间序列（尤其是垂直方向）都呈现显著的非线性运动趋势，且振幅可达1~2 cm[2-5]。因此，地面点相对于参考框架的运动无法用线性速度进行完整描述，仅基于线性速率的模型来维持参考框架会有一定的局限性，而且其精度只能达到厘米级[6-7]。
要精确表达站点坐标，通常采用两种方式描述非线性运动：（1）不考虑引起测站坐标非线性变化的各种物理机制，只根据坐标时间序列本身的运动趋势建模；（2）从产生形变的物理机理入手，分析各因素影响，从而对测站坐标进行改正。国内外对IGS站坐标时间序列的分析已有不少相关研究，也取得较大进展。如文献[8]对中国IGS站坐标时间序列进行频谱分析，发现高程方向周年特性表现明显；文献[9-12]证实了大气负荷、水文负载、非潮汐海洋负载等负载变化是引起GPS测站垂向位移的主要因素。由于引起测站坐标时间序列非线性变化的各种机制比较复杂，目前，国际上还未建立包含多种机制影响的非线性变化理论改正模型，所以本文倾向于对站点本身的运动信息进行建模。
奇异谱分析（singular spectrum analysis，SSA）是时间序列分析中一种常用的方法，在气候学[13-14]、测量学[15]和海洋学[16]等领域应用广泛，它的优越性在于无需任何先验信息和假设条件，能稳定识别和强化周期信号，把最可预报的分量聚集到若干个时间序列中，以此选择若干有意义的分量重建序列，降低噪声影响。因此，特别适合分析有周期震荡的时间序列数据[17]。根据贡献率较大的前几项来对时间序列进行重构，能有效地将半年及半年以上的周期项提取出来，而短周期项如季节周期项、月周期项容易被当作噪声忽略掉。因此，需要用一定方法削弱这种影响。
为更好地从噪声信号中分离出时间序列的周期性、长期性变化信号等有用信息，提高建模精度，本文提出一种小波多尺度分解和SSA相结合的方法提取坐标时间序列中有用信息的本文思路，通过对全球11个测站20年（1999―2018年）的GPS垂向坐标时间序列进行实验分析，验证了利用本文方法能够有效地提取时间序列中的有用信息，提高GPS坐标时间序列的精度。
-
多分辨率分析是小波分析的重要概念，又称为多尺度分析。它能将信号在不同分辨率上进行分解，得到低频部分和高频部分。其中，低频部分能够反映信号的概貌；而高频部分又叫细节部分，能够刻画信号的细节信息。小波多分辨率分析主要是通过空间分解，将信号分解在子空间列中，方便进行信号分析[18-20]，其分解示意图如图 1所示。
图 1 小波多分辨率分解示意图
Figure 1. Schematic Diagram of Wavelet Multi‐resolution Decomposition
图 1中S表示原始时间序列，A1、A2、A3表示小波分解后各层时间序列的低频部分，D1、D2、D3则表示小波分解后各层时间序列的高频部分。分解后的各部分与原始信号具有以下关系：S=A3+D3+D2+D1，以此类推，若要进行进一步的分解，可以把低频部分A3分解成低频部分A4和高频部分D4。因此，可以看出小波变换是几乎无损的信号分析方法。
-
对给定的一维时间序列x=（x1，x2，x3…xn），标准的SSA算法主要由如下4个步骤构成：
1）构造时滞矩阵。选择适当的窗口长度$ L $，$ 1<L<N/2 $，一般情况下，$ L $的选取不宜超过整个数据长度N的1/3，如果根据事先经验大致确定数据的周期特征，则$ L $一般取周期的公倍数。根据向量x构建轨迹矩阵$ \mathit{\boldsymbol{X}} $为：
$$ \mathit{\boldsymbol{X}}=({x}_{ij}{)}_{i, j=1}^{L, K}=\left[\begin{array}{cccc}{x}_{1}& {x}_{2}& \cdots & {x}_{K}\\ {x}_{2}& {x}_{3}& \cdots & {x}_{K+1}\\ ⋮& ⋮& \ddots & ⋮\\ {x}_{L}& {x}_{L+1}& \cdots & {x}_{N}\end{array}\right] $$ (1) 式中，K=N$ - $L+1；$ \mathit{\boldsymbol{X}} $为$ M\times K $阶时滞矩阵，且$ \mathit{\boldsymbol{X}} $中副对角线上的元素相等，即对于$ \mathit{\boldsymbol{X}} $中的元素$ {x}_{i, j} $，有$ {x}_{i, j}={x}_{i-1, j+1} $，因此，$ \mathit{\boldsymbol{X}} $是一个Hankel矩阵。
2）奇异值分解（singular value decomposition，SVD）。对时滞矩阵$ \mathit{\boldsymbol{X}} $进行奇异值分解：
$$ \mathit{\boldsymbol{X}}=\mathit{\boldsymbol{U}}\cdot {\mathit{\pmb{\Lambda }}}^{\frac{1}{2}}\cdot {\mathit{\boldsymbol{V}}}^{\mathrm{T}} $$ (2) 式中，$ {\mathit{\boldsymbol{U}}}_{1}, {\mathit{\boldsymbol{U}}}_{2}\cdots {\mathit{\boldsymbol{U}}}_{d} $分别为对应的特征向量；$ \mathit{\boldsymbol{\Lambda }}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}[{\lambda }_{1}, {\lambda }_{2}\cdots {\lambda }_{d}] $为矩阵$ \mathit{\boldsymbol{X}}{\mathit{\boldsymbol{X}}}^{\mathrm{T}} $和矩阵$ {\mathit{\boldsymbol{X}}}^{\mathrm{T}}\mathit{\boldsymbol{X}} $的前$ d $（$ d=\mathrm{m}\mathrm{i}\mathrm{n}\left\{L, K\right\} $）个最大非零特征值构成的对角矩阵；$ {\mathit{\boldsymbol{\Lambda }}}^{\frac{1}{2}} $表示时滞矩阵$ \mathit{\boldsymbol{X}} $的奇异值。由此推知：
$$ \left\{\begin{array}{l}{\mathit{\boldsymbol{X}}}_{i}=\sqrt{{\lambda }_{i}}{\mathit{\boldsymbol{U}}}_{i}{\mathit{\boldsymbol{V}}}_{i}^{\mathrm{T}}\\ {\mathit{\boldsymbol{V}}}_{i}={\mathit{\boldsymbol{X}}}^{\mathrm{T}}{\mathit{\boldsymbol{U}}}_{i}/\sqrt{{\lambda }_{i}}\end{array}\right. $$ (3) 式中，$ \sqrt{{\lambda }_{i}} $为时滞矩阵$ \mathit{\boldsymbol{X}} $的奇异值，$ \sqrt{{\lambda }_{1}}\ge \sqrt{{\lambda }_{2}}\ge \cdots \ge \sqrt{{\lambda }_{d}}\ge 0 $为奇异谱。则时滞矩阵$ \mathit{\boldsymbol{X}} $的SVD可以表示为：
$$ \mathit{\boldsymbol{X}}={\mathit{\boldsymbol{X}}}_{1}+{\mathit{\boldsymbol{X}}}_{2}+\cdots +{\mathit{\boldsymbol{X}}}_{d} $$ (4) 3）分组。将矩阵$ {\mathit{\boldsymbol{X}}}_{i} $的下标$ \left\{\mathrm{1, 2}\cdots d\right\} $分成$ M $个互不相交的集合$ {I}_{1}, {I}_{2}\cdots {I}_{M} $，设$ I=\left[{i}_{1}, {i}_{2}\cdots {i}_{p}\right] $，那么与集合$ I $相关的矩阵$ {\mathit{\boldsymbol{X}}}_{I} $就可以表示为$ {\mathit{\boldsymbol{X}}}_{I}={\mathit{\boldsymbol{X}}}_{i1}+{\mathit{\boldsymbol{X}}}_{i2}+\cdots +{\mathit{\boldsymbol{X}}}_{ip} $，轨迹矩阵就可以表示成：
$$ \mathit{\boldsymbol{X}}={\mathit{\boldsymbol{X}}}_{I1}+{\mathit{\boldsymbol{X}}}_{I2}+\cdots +{\mathit{\boldsymbol{X}}}_{IM} $$ (5) 一般说来，在SVD中，每个矩阵$ {\mathit{\boldsymbol{X}}}_{I} $对轨迹矩阵$ \mathit{\boldsymbol{X}} $的贡献度与它的特征值有关，其关系可以用下面的式子表示：
$$ {\eta }_{I}=\sum\limits _{i\in I}{\lambda }_{i}/\sum\limits _{i=1}^{L}{\lambda }_{i} $$ (6) 4）对角平均化。对角平均化的目的就是把第3）步分解得到的矩阵$ {\mathit{\boldsymbol{X}}}_{IM} $重新转换成为长度为$ N $的新时间序列，称之为重建成分（reconstruction component，RC），所有RC之和等于原序列。假设$ \mathit{\boldsymbol{z}}=({z}_{1}, {z}_{2}\cdots {z}_{N}) $为$ \mathit{\boldsymbol{z}} $经过对角平均化得到的时间序列，则对角平均化的公式可表示为：
$$ {z}_{n}=\left\{\begin{array}{l}\frac{1}{n}\sum\limits _{j=1}^{n}{z}_{j, n-j+1}, 1\le n\le L\\ \frac{1}{L}\sum\limits _{j=1}^{L}{z}_{j, n-j+1}, L<n<K\\ \frac{1}{N-n+1}\sum\limits _{j=n-K+1}^{L}{z}_{j, n-j+1}\begin{array}{c}\begin{array}{c}, \end{array}\end{array}K\le n\le N\end{array}\right. $$ (7) 所有重建成分RC叠加后的和与原始序列$ x $相同，即：
$$ \mathit{\boldsymbol{x}}={z}_{1}+{z}_{2}+\cdots +{z}_{N}=\sum\limits _{n=1}^{N}{z}_{n} $$ (8) 截取前K个贡献大的成分近似表示原序列，则有[21]：
$$ \widehat{\mathit{\boldsymbol{x}}}={z}_{1}+{z}_{2}+\cdots +{z}_{K}=\sum\limits _{k=1}^{K}{z}_{k}, k=\mathrm{1, 2}\cdots N $$ (9) 上述步骤中，有两个参数需自行选择，分别是窗口长度$ L $和重构阶次$ K $。
-
首先，对原始坐标时间序列S进行小波分解和重构，得到低频部分和多层的高频部分。具体过程如下，对$ S=\left\{{S}_{1}, {S}_{2}\cdots {S}_{n}\right\} $进行小波多尺度分解和重构，得到下式：
$$ S={A}_{B}+{D}_{1}+{D}_{2}+{D}_{3}+\cdots +{D}_{N} $$ (10) 式中，$ {A}_{B} $为低频部分，$ {A}_{B}=\left\{{a}_{B1}, {a}_{B2}\cdots {a}_{Bn}\right\} $为S中低频信号的重构结果；$ {D}_{N} $为多层的高频部分，$ {D}_{1}=\left\{{d}_{11}, {d}_{12}\cdots {d}_{1n}\right\}, {D}_{2}=\left\{{d}_{21}, {d}_{22}\cdots {d}_{2n}\right\}\cdots {D}_{N}=\left\{{d}_{N1}, {d}_{N2}\cdots {d}_{Nn}\right\} $分别为原始坐标时间序列S中第1层到第N层的高频信号的重构结果。
那么，$ {t}_{i} $时刻的原始坐标时间序列S可以表示为：
$$ {S}_{i}={a}_{Bi}+{d}_{1i}+{d}_{2i}+\cdots +{d}_{Ni} $$ (11) 式中，$ {a}_{Bi} $为低频信号在$ {t}_{i} $时刻的值；$ {d}_{1i}, {d}_{2i}\cdots {d}_{Ni} $分别为各层高频信号在$ {t}_{i} $时刻的值。
利用SSA方法对$ {A}_{B}=\left\{{a}_{B1}, {a}_{B2}\cdots {a}_{Bn}\right\} $和$ {D}_{i}=\left\{{d}_{i1}, {d}_{i2}\cdots {d}_{in}\right\}, $i=1，2，3…N分别进行奇异谱分解和重建，得到低频信号在$ {t}_{i} $时刻的拟合值为$ {\widehat{a}}_{Bi} $，各层高频信号在$ {t}_{i} $时刻的拟合值分别为$ {\widehat{d}}_{1i}, {\widehat{d}}_{2i}\cdots {\widehat{d}}_{Ni} $。则原始坐标时间序列S在$ {t}_{i} $时刻的拟合序列为：
$$ {\widehat{S}}_{i}={\widehat{a}}_{Bi}+{\widehat{d}}_{1i}+{\widehat{d}}_{2i}+\cdots +{\widehat{d}}_{Ni} $$ (12) -
实验采用的数据来自IGS网站提供的全球测站坐标时间序列文件，大部分测站时间序列包含1999―2018年近20 a的数据，采样间隔为一周。§1.2使用SSA算法对坐标时间序列数据进行处理时，有窗口长度L和重构阶次K两个参数需要自行选择。
由于测站坐标时间序列中存在的已知周期有周年周期和半年周期，且本文所用数据的时间分辨率为一周，因此，在利用SSA方法提取坐标时间序列的信息时可以选择已知周期的最小公倍数52为窗口长度[22]。
参数K主要根据奇异值的贡献率来确定，若太小，则表示后面部分信号被当作噪声剔除掉；若太大，则会使得部分噪声被当作信号而提取出来。因此，确定窗口长度$ L $与重构阶次K是SSA的关键。
鉴于篇幅，本文以BJFS站垂直方向时间序列为例，根据SSA方法对其进行分析，并将前14阶特征值的贡献率进行了统计，如图 2和表 1所示。
图 2 BJFS站垂向时间序列及其特征值贡献率
Figure 2. Vertical Time Series and Its Eigenvalue Contribution Ratio of BJFS Station
表 1 前14阶特征值贡献率统计
Table 1. Statistics of the First 14 Order Eigenvalue Contribution Ratio
阶次 重构成分 特征值贡献率/% 阶次 重构成分 特征值贡献率/% 1 RRC1 32.30 8 RRC8 1.20 2 RRC2 30.45 9 RRC9 1.13 3 RRC3 4.71 10 RRC10 1.06 4 RRC4 4.68 11 RRC11 0.93 5 RRC5 4.19 12 RRC12 0.83 6 RRC6 1.55 13 RRC13 0.77 7 RRC7 1.36 14 RRC14 0.77 由图 2（b）和表 1可以看出，有多对特征值近似相等，一般来说，近似相等的特征值表示周期与振幅相同的重构分量。其中RRC1和RRC2的特征值贡献率较大，分别为32.30%和30.45%；RRC3、RRC4和RRC5贡献率次之，分别为4.71%、4.68%和4.19%；从第6阶开始，各重构序列的特征值贡献率较低，且阶次越高，贡献率越小。
本文选择将特征值贡献率较大的前6阶特征值的重构成分进行快速傅里叶变换(fast Fourier transform, FFT)提取周期，结果如图 3和图 4所示。
图 3 前6阶特征值的重构成分
Figure 3. Reconstruction Components of the First 6 Order Eigenvalues
图 4 前6阶重构成分的频率谱
Figure 4. Frequency Spectrums of the First 6 Order Reconstruction Components
根据图 4可知，重构序列RRC1和RRC2的周期相同，将其合并可以得到周期为1 a、振幅为4.76 mm的周期项；RRC5的周期分别为0.5 a（与RRC3的周期相同）和9 a（与RRC4的周期相同），将其合并得到周期为0.5 a、振幅为0.95 mm的周期项以及周期为9 a、振幅为0.88 mm的周期项；而RRC6对应的周期为0.3 a（即季节周期项）且振幅影响较小。上述结果与直接利用FFT提取的结果相吻合。
综上可知，在利用SSA对BJFS站垂向时间序列进行分析时，原序列的趋势项和周期项主要分布在前5个重构序列中。因此，本文认为前5阶RC为主要信息成分，可以选择这5个特征值对时间序列进行重构，拟合结果如图 5所示。
图 5 SSA的拟合结果及残差
Figure 5. Fitting Results and Residuals of Singular Spectrum Analysis
从图 5（a）可以看出，重构后的时间序列大致表现出了测站垂直方向时间序列的走势，但是拟合效果不太理想，图 5（b）中残差振幅在3 mm左右，而且还存在部分周期规律。
在对残差序列进行频谱分析时发现，时间序列中还存在一些半年以下的周期项，如图 6所示。
图 6 残差序列的频率谱
Figure 6. Frequency Spectrum of the Residual Series
因此，基于SSA的非线性运动建模方法能有效地将半年及半年以上的周期项提取出来，而短周期项如季节周期项、月周期项则容易被忽略。
-
在对原始坐标时间序列进行小波分解的过程中，还存在着小波基函数和分解层数如何确定的问题。对几种常用的小波基函数及其特性进行了统计，结果如表 2所示。
表 2 常用小波基函数及其特性
Table 2. Common Wavelet Basis Functions and Their Characteristics
小波基函数 支撑长度 消失矩阶数 对称性 特点 Haar小波 1 1 对称 在时域上是不连续的，频率的局部性差 dbN小波 2N N 近似对称 光滑性随着N的增加而增加，性能较好 symN小波 2N N 近似对称 能够减少相位失真，更适合于图像处理领域 Meyer小波 有限长度 对称 对信号进行分析会产生失真现象 bior小波 重构：2Nr+1
分解：2Nd+1Nr－1 不对称 具有线性相位性，在信号与图像的重构中应用比较普遍 由表 2可知，dbN小波具有良好的时间序列分析特性，在处理坐标时间序列数据方面优势明显。因此，本文一般选择dbN小波进行小波分解。
理论上，随着分解层数的增加，测站坐标时间序列被划分得越来越细，各层信号的平稳性和平滑性也越来越好。由于分解过程中存在计算上的误差，以及分解层数越多，会带来计算时间的增加，相应的处理效率也会随之下降，因此，应选择合适的分解层数。
以BJFS站垂直方向时间序列为例，在不同分解层数（2~6层）时，本文方法的均方根误差（root mean square error，RMSE）和平均绝对误差（mean absolute error，MAE）如表 3所示。根据表 3可知，当分解层数超过3层时，新本文方法的拟合精度随着分解层数的增加基本上没有太大的变化；当分解层数为3~4层时，本文方法的之间差异不大。考虑到测站坐标时间序列数据量很大，在保证处理效率的基础上，本文进行小波分解时，分解层数一般选择3层。
表 3 不同分解层数下拟合精度对比
Table 3. Fitting Accuracy Comparison of Different Decomposition Layers
拟合精度 分解层数 2 3 4 5 6 RMSE/mm 2.15 1.88 1.80 1.80 1.80 MAE/mm 1.68 1.49 1.43 1.44 1.44 根据表 2和表 3确定的小波基函数以及分解层数，实验利用具有正交性和高度紧支撑性的db4小波将BJFS站垂向时间序列分解3层得到各层低频信息和高频信息，如图 7所示。
图 7 小波分解后的各层低频部分和高频部分
Figure 7. Low-Frequency Part and High-Frequency Part of Each Layer After Wavelet Decomposition
a3+d1+d2+d3之和即为原序列，a3为小波分解后的低频部分，称为近似系数，反映出了时间序列的概貌，表现出较强的规律性（主要包含了趋势项和周期项）；d1、d2、d3分别为小波分解后的各层高频部分，称为细节系数，表现出的随机项成分较多。按照§2.1中的特征值选取方法，对小波分解后的低频部分a3和高频部分d1、d2、d3分别选择贡献率较大的前5阶特征值进行序列重构，结果如图 8~图 11所示。
图 8 近似系数a3的特征值贡献率及拟合效果
Figure 8. Eigenvalue Contribution Ratio and Fitting Effect of Approximation Coefficient a3
图 9 细节系数d3的特征值贡献率及拟合效果
Figure 9. Eigenvalue Contribution Ratio and Fitting Effect of Detail Coefficient d3
图 10 细节系数d2的特征值贡献率及拟合效果
Figure 10. Eigenvalue Contribution Ratio and Fitting Effect of Detail Coefficient d2
图 11 细节系数d1的特征值贡献率及拟合效果
Figure 11. Eigenvalue Contribution Ratio and Fitting Effect of Detail Coefficient d1
从图 8~图 11中可以看出，经过小波3层分解后的各层时间序列相对平稳了许多。从近似系数a3和细节系数d3中可以看出，重构序列和原始序列在大部分时间段基本上实现重合，而在细节系数d1和d2中表现出来的随机项成分较多，周期项较少。因此，各特征值贡献率差异不大，相应的拟合效果要小于其他层。将各层时间序列的拟合结果进行叠加得到最终的拟合结果，并将其与原始测站坐标时间序列进行比较，效果如图 12所示。
图 12 小波分解和SSA的拟合结果及残差
Figure 12. Fitting Results and Residuals of Wavelet Decomposition and Singular Spectrum Analysis
从图 12（a）中可以看出，经过小波多尺度分解和SSA相结合的方法进行拟合后的时间序列不仅能够反映出测站坐标时间序列的走势，而且两者之间的差异明显减小。从剩余残差信息图 12（b）中可以看出，本文方法拟合后的残差振幅在2 mm左右，拟合精度有一定程度的提高。
在对残差序列进行频谱分析时发现，一些短周期项如季节周期项、月周期项的影响明显减小，如图 13所示。
图 13 残差序列的频率谱
Figure 13. Frequency Spectrum of the Residual Series
与单纯的SSA方法相比，本文方法能更有效地对各层时间序列的周期项进行提取，减小部分短周期项如季节周期项、月周期项容易被忽略的影响。
-
本文在全球范围内低、中、高纬度各选取了3~4个测站，利用§2.1和§2.2中的方法进行分析，通过计算剩余残差序列的MAE与RMSE，对两种方法在测站垂向时间序列的拟合效果进行了评估，结果如表 4所示。
表 4 两种方法的垂直方向拟合精度对比
Table 4. Comparison of Fitting Accuracy Between Two Models in the Vertical Direction
测站名称 纬度 SSA/mm 小波SSA/mm RMSE MAE RMSE MAE ADIS 9.2°N 2.59 1.98 1.90 1.44 TUVA 8.3°S 2.94 2.17 2.23 1.71 NAUR 0.3°S 2.82 2.09 2.09 1.61 IISC 13.1°N 2.73 2.00 1.98 1.47 WUHN 30.5°N 2.90 2.31 2.11 1.66 STR1 35.2°S 2.50 1.81 1.83 1.38 BJFS 39.6°N 2.48 1.93 1.88 1.49 MAC1 54.3°S 2.12 1.65 1.65 1.27 HOFN 64.2°N 2.13 1.67 1.66 1.30 KELY 66.6°N 2.64 2.05 1.82 1.43 MAW1 67.4°S 2.19 1.68 1.46 1.14 剩余残差序列的RMSE和MAE分别为：
$$ \left\{\begin{array}{l}{R}_{\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}}=\sqrt{\frac{\sum\limits _{i=1}^{n}({Y}_{i}-{\widehat{Y}}_{i}{)}^{2}}{n}}\\ {M}_{\mathrm{M}\mathrm{A}\mathrm{E}}=\frac{1}{n}\sum\limits _{i=1}^{n}|{Y}_{i}-{\widehat{Y}}_{i}|\end{array}\right. $$ (13) 式中，$ {\widehat{Y}}_{i} $是真实值$ {Y}_{i} $的拟合值；$ n $是样本容量或者数据个数。
根据表 4和图 5、图 12可以看出，在测站垂直方向上，利用小波多尺度分解和SSA相结合的方法对坐标时间序列进行拟合时，其MAE与RMSE皆明显小于仅基于SSA方法的MAE与RMSE。
经过计算，与SSA方法的拟合结果相比，本文方法的RMSE和MAE分别降低了约26.5%和25.5%。因此，可以得出与上述分析相同的结论，即在相同特征值的情况下，小波分解和SSA的拟合效果要优于仅基于SSA的拟合效果。
由于经过小波分解后，坐标时间序列与原始序列相比在频率成分上相对单一、平滑，并且原始序列中的趋势项和周期项主要集中在低频部分中，而随机项则大部分存在于高频部分中。
基于SSA将贡献率较大的前几项特征值进行时间序列重建的方法，会使一部分短周期项被当作噪声剔除掉。而本文方法可以将时间序列进行分层，在每一层上根据所选的特征值重建序列，在一定程度上减少了部分短周期项被当作噪声剔除的概率，从而进一步提高了建模精度。可以看出，采用小波分解和SSA相结合的方法进行非线性运动建模有很好的自适应性，并且建模精度较高。
-
首先，本文提出一种小波多尺度分解和SSA相结合的非线性运动建模方法，通过小波分解将GPS坐标时间序列分解为各层高频和低频部分；然后，利用SSA对各层时间序列进行分析，截取贡献率较大的前几个特征值重构序列，并将各层时间序列的拟合结果进行叠加得到最终的拟合结果；最后，对本文方法的拟合结果进行了分析与评估，得出如下结论。
1）与单纯的SSA方法相比，小波多尺度分解和SSA相结合的非线性运动建模方法能有效地融合两种方法的优点。本文方法将原始序列分解成频率成分相对单一、平滑的各层高频部分和低频部分，在每一层上根据所选的特征值重建序列，在一定程度上减少了部分周期项（如季节周期项、月周期项）被当作噪声剔除的概率。结果表明，本文方法的建模精度提高了约26%。
2）小波分解虽然具有较大的优势，即可以将测站坐标时间序列分解到不同的尺度上，在不同尺度上分别进行SSA建模，但是要注意选择合适的分解层数。考虑到数据处理效率和精度要求，分解层数一般定为3层即可。
3）小波分解和SSA相结合的方法对坐标时间序列建模具有很好的自适应性，由于测站所处地理位置不同，坐标时间序列所表现出来的规律也不尽相同，各测站前K阶特征值贡献率的大小就存在差异。因此，不宜选取相同的特征值个数对所有的测站进行处理。如何合理地选择特征值个数是今后需要进一步探讨的方向。
Application of Wavelet Decomposition and Singular Spectrum Analysis to GNSS Station Coordinate Time Series
-
摘要: 为了有效地提取GNSS（global navigation satellite system）站坐标时间序列中的有用信息，提高坐标时间序列的建模精度，提出一种小波多尺度分解与奇异谱分析相结合的非线性运动建模方法，并利用全球11个测站20年（1999―2018年）的GPS（global positioning system）垂向坐标时间序列对所提方法进行了验证。首先，通过小波分解将坐标时间序列分解到不同尺度上; 然后，对分解后的各层高频部分和低频部分进行奇异谱分析; 最后，通过叠加合成得到原始坐标时间序列的拟合值，并对所提方法的拟合效果进行评估。结果表明，与单纯的奇异谱分析方法相比，所提方法能够更加准确地从含噪声的有限尺度时间序列中提取趋势和周期等有用信息，降低了部分周期项如季节周期项、月周期项被当作噪声剔除的概率，并且建模精度有26%的提高。Abstract:
Objectives In order to effectively extract useful information from time series of the global navigation satellite system(GNSS) sites, and improve the modeling accuracy of coordinate time series, this paper proposes a nonlinear motion modeling method combining wavelet decomposition and singular spectrum analysis. Experiments were carried out using global positioning system(GPS) vertical coordinate time series of 11 stations around the world from 1999 to 2018. Methods Firstly, the coordinate time series is decomposed into different scales by wavelet decomposition.And then the singular spectrum analysis(SSA) is performed on the high-frequency part and the low-frequency part of each layer. Finally, combine the fitting value that is the fitting coordinate time series, and the fitting effect of the new method was evaluated. Results The results show that, compared with the simple singular spectrum analysis method, the new method can extract useful information such as trend and period more accurately from the limited scale of the time series with noise. And reduce the partial period items in the singular spectrum analysis method to some extent, for example, the seasonal period item and the monthly period item are regarded as the probability of noise rejection, and the modeling accuracy is improved by 26%. Conclusions The purpose is to propose a nonlinear motion modeling method based on wavelet decomposition and singular spectrum analysis, so as to improve the modeling accuracy of coordinate time series. -
Key words:
- wavelet decomposition /
- singular spectrum analysis(SSA) /
- time series /
- nonlinear /
- contribution ratio
-
图 2 BJFS站垂向时间序列及其特征值贡献率
Figure 2. Vertical Time Series and Its Eigenvalue Contribution Ratio of BJFS Station
图 4 前6阶重构成分的频率谱
Figure 4. Frequency Spectrums of the First 6 Order Reconstruction Components
图 7 小波分解后的各层低频部分和高频部分
Figure 7. Low-Frequency Part and High-Frequency Part of Each Layer After Wavelet Decomposition
图 8 近似系数a3的特征值贡献率及拟合效果
Figure 8. Eigenvalue Contribution Ratio and Fitting Effect of Approximation Coefficient a3
图 9 细节系数d3的特征值贡献率及拟合效果
Figure 9. Eigenvalue Contribution Ratio and Fitting Effect of Detail Coefficient d3
图 10 细节系数d2的特征值贡献率及拟合效果
Figure 10. Eigenvalue Contribution Ratio and Fitting Effect of Detail Coefficient d2
图 11 细节系数d1的特征值贡献率及拟合效果
Figure 11. Eigenvalue Contribution Ratio and Fitting Effect of Detail Coefficient d1
图 12 小波分解和SSA的拟合结果及残差
Figure 12. Fitting Results and Residuals of Wavelet Decomposition and Singular Spectrum Analysis
表 1 前14阶特征值贡献率统计
Table 1. Statistics of the First 14 Order Eigenvalue Contribution Ratio
阶次 重构成分 特征值贡献率/% 阶次 重构成分 特征值贡献率/% 1 RRC1 32.30 8 RRC8 1.20 2 RRC2 30.45 9 RRC9 1.13 3 RRC3 4.71 10 RRC10 1.06 4 RRC4 4.68 11 RRC11 0.93 5 RRC5 4.19 12 RRC12 0.83 6 RRC6 1.55 13 RRC13 0.77 7 RRC7 1.36 14 RRC14 0.77 
下载: 导出CSV
表 2 常用小波基函数及其特性
Table 2. Common Wavelet Basis Functions and Their Characteristics
小波基函数 支撑长度 消失矩阶数 对称性 特点 Haar小波 1 1 对称 在时域上是不连续的，频率的局部性差 dbN小波 2N N 近似对称 光滑性随着N的增加而增加，性能较好 symN小波 2N N 近似对称 能够减少相位失真，更适合于图像处理领域 Meyer小波 有限长度 对称 对信号进行分析会产生失真现象 bior小波 重构：2Nr+1
分解：2Nd+1Nr－1 不对称 具有线性相位性，在信号与图像的重构中应用比较普遍
下载: 导出CSV
表 3 不同分解层数下拟合精度对比
Table 3. Fitting Accuracy Comparison of Different Decomposition Layers
拟合精度 分解层数 2 3 4 5 6 RMSE/mm 2.15 1.88 1.80 1.80 1.80 MAE/mm 1.68 1.49 1.43 1.44 1.44
下载: 导出CSV
表 4 两种方法的垂直方向拟合精度对比
Table 4. Comparison of Fitting Accuracy Between Two Models in the Vertical Direction
测站名称 纬度 SSA/mm 小波SSA/mm RMSE MAE RMSE MAE ADIS 9.2°N 2.59 1.98 1.90 1.44 TUVA 8.3°S 2.94 2.17 2.23 1.71 NAUR 0.3°S 2.82 2.09 2.09 1.61 IISC 13.1°N 2.73 2.00 1.98 1.47 WUHN 30.5°N 2.90 2.31 2.11 1.66 STR1 35.2°S 2.50 1.81 1.83 1.38 BJFS 39.6°N 2.48 1.93 1.88 1.49 MAC1 54.3°S 2.12 1.65 1.65 1.27 HOFN 64.2°N 2.13 1.67 1.66 1.30 KELY 66.6°N 2.64 2.05 1.82 1.43 MAW1 67.4°S 2.19 1.68 1.46 1.14
下载: 导出CSV
-
[1] 马俊, 姜卫平, 邓连生, 等. GPS坐标时间序列噪声估计及相关性分析[J]. 武汉大学学报·信息科学版, 2018, 43 (10): 1 451-1 457 doi: 10.13203/j.whugis20160543 Ma Jun, Jiang Weiping, Deng Liansheng, et al. Estimation Method and Correlation Analysis for Noise in GPS Coordinate Time Series[J]. Geomatics and Information Science of Wuhan University, 2018, 43 (10): 1 451-1 457 doi: 10.13203/j.whugis20160543 [2] 马俊, 姜卫平, 周晓慧, 等. 联合小波和方差分量估计方法分析中国IGS测站时间序列变化特征[J]. 武汉大学学报·信息科学版, 2018, 43 (4): 629-636 doi: 10.13203/j.whugis20150731 Ma Jun, Jiang Weiping, Zhou Xiaohui, et al. Coordinate Time Series of IGS Station in China Based on the Combination of Wavelet Spectral and Variance Component Estimation[J]. Geomatics and Information Science of Wuhan University, 2018, 43 (4): 629-636 doi: 10.13203/j.whugis20150731 [3] Tregoning P, Watson C. Atmospheric Effects and Spurious Signals in GPS Analyses[J]. Journal of Geophysical Research Solid Earth, 2009, 114(B9), DOI: 10.1029/2009JB006344 [4] 朱文耀, 符养, 李彦. GPS高程导出的全球高程振荡运动及季节变化[J]. 中国科学: 地球科学, 2003, 33(5): 470-481 Zhu Wenyao, Fu Yang, Li Yan. Global Elevation Oscillation Movement and Seasonal Change Derived from GPS Elevation[J]. Science in China: Earth Science, 2003, 33(5): 470-481 [5] Amiri-Simkooei A R, Tiberius C C J M, Teunissen P J G. Assessment of Noise in GPS Coordinate Time Series: Methodology and Results[J]. Journal of Geophysical Research Solid Earth, 2007, 112(B7): 1-19 [6] 姜卫平, 马一方, 邓连生, 等. 毫米级地球参考框架的建立方法与展望[J]. 测绘地理信息, 2016, 41(4): 1-6 Jiang Weiping, Ma Yifang, Deng Liansheng, et al. Establishment of mm-Level Terrestrial Reference Frame and Its Prospect[J]. Journal of Geomatics, 2016, 41(4): 1-6 [7] Bogusz J, Klos A. On the Significance of Periodic Signals in Noise Analysis of GPS Station Coordinates Time Series[J]. GPS Solutions, 2016, 20(4): 655-664 doi: 10.1007/s10291-015-0478-9 [8] 张鹏, 蒋志浩, 秘金钟, 等. 我国GPS跟踪站数据处理与时间序列特征分析[J]. 武汉大学学报·信息科学版, 2007, 32(3): 251-254 http://ch.whu.edu.cn/article/id/1848 Zhang Peng, Jiang Zhihao, Bei Jinzhong, et al. Data Processing and Time Series Analysis for GPS Fiducial Stations in China[J]. Geomatics and Information Science of Wuhan University, 2007, 32(3): 251-254 http://ch.whu.edu.cn/article/id/1848 [9] van Dam T M, Wahr J, Milly P C D, et al. Crustal Displacement due to Continent Water Loading[J]. Geophysical Research Letters, 2001, 28(4): 651‑654 doi: 10.1029/2000GL012120 [10] 王敏, 沈正康, 董大南. 非构造形变对GPS连续站位置时间序列的影响和修正[J]. 地球物理学报, 2005, 48(5): 1 045-1 052 Wang Min, Shen Zhengkang, Dong Danan. Effects of Non-tectonic Crustal Deformation on Continuous GPS Position Time Series and Correction to Them[J]. Chinese Journal of Geophysics, 2005, 48(5): 1 045-1 052 [11] 张诗玉, 钟敏, 唐诗华, 等. 我国GPS基准站地壳垂直形变的大气负荷效应[J]. 武汉大学学报·信息科学版, 2006, 31(12): 1 090-1 093 http://ch.whu.edu.cn/article/id/2624 Zhang Shiyu, Zhong Min, Tang Shihua, et al. Atmospheric Load Effect of Crustal Vertical Deformation at GPS Reference Stations in China[J]. Geomatics and Information Science of Wuhan University, 2006, 31(12): 1 090-1 093 http://ch.whu.edu.cn/article/id/2624 [12] Jiang Weiping, Li Zhao, van Dam T, et al. Comparative Analysis of Different Environment Loading Methods and Their Impacts on the GPS Height Time Series[J]. Journal of Geodesy, 2013, 87(7): 687-703 doi: 10.1007/s00190-013-0642-3 [13] Wyatt M, Kravtsov S, Tsonis A. Atlantic Multidecadal Oscillation and Northern Hemisphere's Climate Variability[J]. Climate Dynamics, 2012, 38(5): 929-949 [14] Chen Q, Dam T V, Sneeuw N, et al. Singular Spectrum Analysis for Modeling Seasonal Signals from GPS Time Series[J]. Journal of Geodynamics, 2013, 72(12): 25-35 [15] Zabalza J, Ren J, Zheng W, et al. Singular Spectrum Analysis for Effective Feature Extraction in Hyperspectral Imaging[J]. IEEE Geoscience and Remote Sensing Letters, 2014, 11(11): 1 886-1 890 doi: 10.1109/LGRS.2014.2312754 [16] Kondrashov D, Berloff P. Stochastic Modeling of Decadal Variability in Ocean Gyres[J]. Geophysical Research Letters, 2015, 42(5): 1 543-1 553 doi: 10.1002/2014GL062871 [17] 郭金运, 高文宗, 于红娟, 等. 基于奇异谱分析的静态相对重力观测重力固体潮提取[J]. 地球物理学报, 2018, 61(10) : 3 889-3 902 Guo Jinyun, Gao Wenzong, Yu Hongjuan, et al. Gravity Tides Extracted from Relative Gravimetric Data with Singular Spectrum Analysis[J]. Chinese Journal of Geophysics, 2018, 61(10): 3 889-3 902 [18] 文鸿雁. 基于小波理论的变形分析模型研究[D]. 武汉: 武汉大学, 2004 Wen Hongyan. Research on Deformation Analysis Model Based on Wavelet Transform Theory[D]. Wuhan: Wuhan University, 2004 [19] Ghaderpour E, Ince E S, Pagiatakis S. Least-Squares Cross-Wavelet Analysis and Its Applications in Geophysical Time Series[J]. Journal of Geodesy, 2018, 92(10): 1 223-1 236 doi: 10.1007/s00190-018-1156-9 [20] 范朋飞. 高精度GPS站点坐标时间序列分析与应用[D]. 西安: 长安大学, 2013 Fan Pengfei. Analysis and Application of High-Precision Coordinate Time Series of GPS Site[D].Xi'an: Chang'an University, 2013 [21] 王解先, 连丽珍, 沈云中. 奇异谱分析在GPS站坐标监测序列分析中的应用[J]. 同济大学学报(自然科学版), 2013, 41(2): 282-288 Wang Jiexian, Lian Lizhen, Shen Yunzhong. Application of Singular Spectral Analysis to GPS Station Coordinate Monitoring Series[J]. Journal of Tongji University(Natural Science), 2013, 41(2): 282-288 [22] 周茂盛, 郭金运, 沈毅, 等. 基于多通道奇异谱分析的GNSS坐标时间序列共模误差的提取[J]. 地球物理学报, 2018, 61(11): 4 383-4 395 Zhou Maosheng, Guo Jinyun, Shen Yi, et al. Extraction of Common Mode Errors of GNSS Coordinate Time Series Based on Multi-channel Singular Spectrum Analysis[J]. Chinese Journal of Geophysics, 2018, 61(11): 4 383- 4 395 -
点击查看大图
计量
- 文章访问数: 446
- HTML全文浏览量: 83
- PDF下载量: 91
- 被引次数: 0
