在GPS短基线测量中，通过卫星、历元之间的差分以及双频组合等原理，可以消除GPS测量的大部分误差，如卫星轨道误差、卫星钟差、对流层延迟、电离层延迟等^{[1, 2]}。这些误差可以根据测站、卫星之间误差的相关性进行消除，然而多路径误差与测站的周围环境、接收机本身的结构、性能相关，不同测站之间的多路径误差没有相关性，不能采取上述方法进行消除。现阶段主要通过接收机信号处理、天线结构设计、观测数据后期处理^[3-5]3种方法来减弱多路径误差。在GPS变形监测中，建筑物周围环境基本保持不变，多路径效应具有很强的周期性，可以通过滤波的方法去除观测噪声，提取多路径效应的周期变化。奇异谱分析 (singular spectrum analysis, SSA) 是从Karhumen-Loeve理论上发展起来的，是从时间序列相空间动力结构出发，以经验正交的方式展开时间序列子空间，不基于正弦波假定前提，不需要先验信息，能够准确识别时间序列中的周期信号^{[6, 7]}。

在使用GPS进行高层建筑变形监测时，根据奇异谱分析理论提取多路径效应误差，可以提高变形监测的精度，获取建筑物真实的变形信息，提高GPS技术在变形监测中的适用性。

1 奇异谱分析原理

奇异谱分析是将序列长度为N的一维时间序列按照给定的嵌套维数及一定的时间滞后量排列成二维的时滞矩阵，在时滞矩阵中对时间序列相空间进行主成分分析，获取矩阵的特征值和特征向量^{[8, 9]}。特征向量通常具有正交特性，对应着原有时间序列中的不同波动信号。将给定的时间序列x按照嵌套维数M、时滞为1排列成时滞矩阵，即

$\mathit{\boldsymbol{X}}{\rm{ = }}\left[ {\begin{array}{*{20}{c}} {{x_1}}& \cdots &{{x_{N - M + 1}}}\\ \vdots &{}& \vdots \\ {{x_M}}& \cdots &{{x_N}} \end{array}} \right] \equiv \left[ {\begin{array}{*{20}{c}} {{x_{10}}}& \cdots &{{x_{1,N - M}}}\\ \vdots &{}& \vdots \\ {{x_{M0}}}& \cdots &{{x_{M,N - M}}} \end{array}} \right]$

(1)

其中，N为序列长度；M为嵌入维数。大量实践证明，M一般取N/3左右，分析效果较为理想^[10]。矩阵X的第i个状态量为：

${\mathit{\boldsymbol{X}}_i}{\rm{ = }}\left[ {\begin{array}{*{20}{c}} {{x_{i + 1}}}\\ {{x_{i + 2}}}\\ \vdots \\ {{x_{i + M}}} \end{array}} \right] \equiv \left[ {\begin{array}{*{20}{c}} {{x_{1i}}}\\ {{x_{2i}}}\\ \vdots \\ {{x_{Mi}}} \end{array}} \right],i = 0,1,2, \cdots ,N - M$

(2)

共N-M+1个状态，矩阵中的元素与原来的时间序列元素的对应关系为x_ji=x_j+1。式 (1) 的变量间协方差是原序列x_i不同滞后元素的自协方差。构造滞后自协方差阵T_x：

${\mathit{\boldsymbol{T}}_x} = \left[ {\begin{array}{*{20}{c}} {C\left( 0 \right)}&{C\left( 1 \right)}&{C\left( 2 \right)}& \cdots &{C\left( {M - 1} \right)}\\ {C\left( 1 \right)}&{C\left( 0 \right)}&{C\left( 1 \right)}& \cdots & \vdots \\ {C\left( 2 \right)}&{C\left( 1 \right)}&{C\left( 0 \right)}& \cdots & \vdots \\ \vdots & \vdots & \vdots &{C\left( 0 \right)}&{C\left( 1 \right)}\\ {C\left( {M - 1} \right)}&{C\left( {M - 2} \right)}& \cdots &{C\left( 1 \right)}&{C\left( 0 \right)} \end{array}} \right]$

(3)

T_x是Toeplitz矩阵，同时也是实对称矩阵。主对角线元素是时间序列x的方差 (或时滞为0的自协方差)。C(j) 为时间序列x滞后为j的自协方差，0≤j≤M－1。C(j) 用Yule-Walke估计法得到^[11]：

$C\left( j \right) = \frac{1}{{N - j}}\sum\limits_{i = 1}^{N - j} {{x_i}{x_{i + j}}} ,j = 0,1,2, \cdots ,M - 1$

(4)

然后根据公式：

${\mathit{\boldsymbol{T}}_x}{\mathit{\boldsymbol{E}}^k} = {\lambda _k}{\mathit{\boldsymbol{E}}^k},k = 1,2,3, \cdots ,M$

(5)

求得T_x的特征值和特征向量。T_x的特征向量E^k就是M个分量构成的一个时间序列，它反映时间序列x中的时间演变型。定义状态向量在第M个特征向量上的投影为：

$a_i^k = {\mathit{\boldsymbol{X}}_i}{\mathit{\boldsymbol{E}}^k} = \sum\limits_{j = 1}^M {{x_{i + j}}E_j^k} ,0 \le i \le N - M$

(6)

可以由其中一部分特征向量和时间系数来重建x的成分：

$x_i^k = \left\{ {\begin{array}{*{20}{l}} {\frac{1}{M}\sum\limits_{j = 1}^M {a_{i - j}^kE_j^k} ,M \le i \le N - M + 1}\\ {\frac{1}{i}\sum\limits_{j = 1}^i {a_{i - j}^kE_j^k} ,1 \le i \le M - 1}\\ {\frac{1}{{N - i + 1}}\sum\limits_{j = i - N + M}^M {a_{i - j}^kE_j^k} ,N - M + 2 \le i \le N} \end{array}} \right.$

(7)

所有重建成分 (reconstructed components，RC)x^k之和与原序列相等。与主成分分析相似，将T_x的特征值λ_k从大到小排列，即λ₁≥λ₂≥λ₃≥…≥λ_M≥0。截取前p个较大的特征值，由其所对应的x^k之和得到的重建序列可充分反映原序列的整体特征，并且能够去除原有序列中的噪声和随机误差部分。即

$x \approx y = \sum\limits_{k = 1}^p {{x^k}} ,1 \le p \le M$

(8)

式中，p为截取特征值序列的长度。

2 仿真实验

为了验证奇异谱分析在去除噪声以及提取时间序列中波动信息的准确性和有效性，将一组谐波信号以及随机噪声叠加，利用奇异谱分析方法对仿真信号进行滤波，提取仿真信号中的谐波信号。仿真序列长度为1 000，其中三组谐波信号分别为y₁=cos (x/5)，y₂=cos (x/10)，y₃=cos (x·2π/10)；随机噪声e服从N(0, 1)，则将加噪声的仿真序列和不加噪声的仿真序列分别记为：

$\left\{ {\begin{array}{*{20}{l}} {y = {y_1} + {y_2} + {y_3} + e}\\ {w = {y_1} + {y_2} + {y_3}} \end{array}} \right.$

(9)

将加噪声的仿真序列进行滤波，N=1 000，M=350。将得到的奇异值进行大小排序，并以最大奇异值为1，其后每个奇异值与最大值相除，得到每个奇异值对应子序列的贡献率。图 1给出了仿真序列的奇异值图。

从图 1中可以看出，奇异值贡献率在前4个序列中占比较大，第5个以后，奇异值贡献率急剧下降，几乎接近于0，表明第5个以后对应的子序列在原始序列中的占比可以忽略不计，因此取前4个奇异值对应的子序列进行重建，得到重建序列，并将重建序列与原始加噪声序列做差，观察残余序列的特征。图 2给出了原始无噪声序列、加噪声序列、重建序列以及重建序列与加噪声序列做差得到的残余序列图。

从图 2中可以看出，加噪声序列的曲线很不光滑，含有大量不规则噪声。而在利用奇异谱分析进行滤波以后，重建序列与原始无噪声序列几乎相同，相差非常微小，而残余序列分布为白噪声序列，正是加入噪声的特性，因此证明了奇异谱分析在过滤不规则噪声、提取时间序列真实波动信息的有效性和准确性。

3 实例分析

在某高层建筑上布设监测点和基准点^[12]，采用不加扼流圈的GPS接收机进行测量。根据测量结果，GPS监测平面误差在X、Y方向基本一致，以Y方向的位移量为例进行奇异谱分析处理。图 3为奇异谱分析得到的奇异值贡献率图，从图中可以看出奇异值分布较好，并没有出现拖尾现象，在10以后下降到0.1以下，在20以后几乎接近于0，其奇异值对应的子序列所占比重很小，可以忽略不计。一次取前10个序列作为重建序列。

图 4给出了实测数据Y方向的奇异谱分析滤波重建的结果图，图 4(a)中折线为原始位移序列，曲线为1~10子序列组成的重建序列；图 4(b)为重建序列与原始序列的差值序列。从图 4中可以看出，原始序列含有大量白噪声以及不规则扰动，位移序列精度较差，与建筑物真实变形信息相差较大，而重建序列则较为光滑，且贴合原始序列，准确反映了建筑物的位移信息。

为统计奇异谱分析滤波后GPS坐标的真实精度，以第7~10次观测成果为例进行数据分析，将高精度测量机器人徕卡TS30观测坐标序列与重建序列进行对比，进行精度统计。表 1给出了重建序列的精度统计结果。

经计算，通过奇异谱滤波得到的重建序列误差绝对值平均值为1.02 cm；以徕卡TS30观测结果为真值，计算GPS测量成果中误差为±1.15 cm。从表 1可以看出，第8期观测中T2、T3监测点较差明显偏大，经分析，测量时受到临时构筑物遮挡，部分遮挡了GPS信号，造成误差较大。剔除第8期T2、T3两点后，重建序列误差绝对值平均值为0.91 cm，中误差为±0.84 cm，可以满足《工程测量规范》四等变形监测的精度要求^[13]。

当然，在更高精度的变形监测中，采用设置强制对中观测墩、延长观测时间、使用多频多模GPS接收机观测等方法，进一步提高观测精度，可以达到《工程测量规范》三等变形监测的精度要求^[13]。

4 结束语

本文介绍了奇异谱分析的基本原理，并利用仿真数据对其准确性和精度进行验证。结果表明，奇异谱分析能够有效去除仿真序列中的随机噪声，高精度地还原了原始序列。利用奇异谱分析对某高层建筑GPS实测Y方向位移量进行滤波，根据奇异值贡献图信息，将1~10子序列作为重建序列，并将重建序列与原始序列做差，得到残余序列。将原始序列与重建序列进行对比，可知重建序列较为光滑，更加符合建筑物的缓慢连续变形特征，奇异谱分析有效去除了观测序列中的多路径误差。与测量机器人监测结果对比进行精度验证，结果显示，重建序列精度能够满足变形监测的精度要求。