GPS具有高精度、连续自动监测以及实时动态的特点，已经成为一种非常有效的高精度动态变形监测技术^{[1, 2]}。高精度变形监测基线较短 (小于3 km)，通常采用双差技术消除轨道误差、电离层延迟、对流层延迟、钟差等误差项。然而，由于多路径误差在测站间的相关性不强，不能通过求差的方式予以消除。因此，在短基线的动态变形监测中，得到的动态三维坐标序列的误差可认为主要为多路径效应以及随机噪声。

大量研究表明，多路径效应与天线的位置、周围环境以及卫星的几何结构密切相关^[1-4]。恒星日滤波利用GPS卫星在一个恒星日内 (23 h 56 m 04 s)，卫星几何结构与前一个恒星日相同的特点，对动态变形监测中的坐标残差序列进行滤波处理，削弱多路径误差影响，达到提取形变信号的目的。

随着对恒星日滤波算法研究的深入^[5-12]，国内外学者发现GPS卫星轨道周期并非严格遵循设计的恒星日时间，并且各颗卫星实际运行周期也各不相同。常规恒星日滤波算法主要致力于采用卫星星历对卫星重复周期的精确求解，忽略了坐标残差时间序列实际的重复性，滤波达不到最佳效果。本文从坐标残差时间序列本身出发，采用滑动窗口分段求取N、E、U 3个方向上的序列重复时间，通过相关系数法计算各颗卫星的实际时延，进而对残差序列进行滤波处理，有效地消除了动态变形监测中多路径效应的影响。

1 坐标残差时间序列相关性分析

常规GPS变形监测过程中，每个恒星日内测站周围的环境基本保持不变，卫星与测站的几何构型基本相同，并且数据处理时采用的处理策略也完全相同。因此，在相邻两个恒星日之间，测站的动态坐标时间序列表现出极大的相似性。如图 1为利用UNAVCO网站 (ftp://data-out.unavco.org/pub) 提供的1 Hz数据求得测站MIDA在2012年年积日为002-003、UTC为4:00~8:00两天里4 h的N方向坐标残差时间序列。

从图 1中可以看出，两天里该测站的坐标残差时间序列表现出极强的相关性，并且年积日003天的序列相对于002天有一定的负时延。由于在短时间内，测站的位置并未发生实质性的移动，那么坐标时间序列表现出的波动性主要为多路径效应以及高频的随机噪声。如果测站在监测过程中受到瞬时外部环境的激励影响产生瞬时移动，那么通过这种重复性去除多路径效应，能够很好地提取出形变信号。

本文采用相关系数法来研究两恒星日之间坐标序列的相关延迟。相关系数是变量之间相关程度的指标，考察的是两个变量之间的关联程度。也就是说，当某一个变量发生变化时，另一个变量会发生什么变化。相关系数计算公式为：

$ r = \frac{{\sum\limits_{i = 1}^n {({X_i}- \overline X )({Y_i}- \overline Y )} }}{{\sqrt {\sum\limits_{i = 1}^n {{{({X_i}- \overline X )}^2} \cdot \sum\limits_{i = 1}^n {{{({Y_i} - \overline Y )}^2}} } } }} = \frac{{\sum {XY - \frac{{\sum X \cdot \sum Y }}{n}} }}{{\sqrt {\left[{\sum {{X^2}-\frac{{{{(\sum X )}^2}}}{n}} } \right]\left[{\sum {{Y^2}-\frac{{{{(\sum Y )}^2}}}{n}} } \right]} }} $

(1)

式中, X、Y分别是X、Y的平均值。

由于动态解算的坐标时间序列受到高频噪声的影响，两序列之间求得的相关性系数不准确，因此需要对其进行消噪处理。本文采用移动平均的方法来平滑序列^[13]，剔除随机噪声。移动平均的具体公式为：

$ \begin{gathered} \left( {{y_i} = {x_{i- [n/2]}} + {x_{i - [n/2] + 1}} + \cdots } \right. \hfill \\ \left. { + {x_i} + \cdots + {x_{i + [n/2] - 1}} + {x_{i + [n/2]}}} \right)/n -\overline x \hfill \\ \end{gathered} $

(2)

式中，x_i为第i历元的值；[]表示取整；n为移动平均周期；x为x_i总的平均值。应该注意到，序列的移动平均周期应与多路径效应的影响周期相一致。本文经过大量试验分析，认为取n=10时效果最好。

利用上述移动平均方法对图 1两序列进行平滑，然后将年积日002天N、E、U 3个方向以236为中心左右两边各取30个历元作为负延迟进行平移，并采用相关系数法求两序列相关系数，获取最大相关系数时的相关延迟。

从表 1中可以看出，N、E、U 3个方向上，相邻两个恒星日的坐标残差序列均表现为极强相关性。相关系数最大时的相关延迟与设计延迟236 s左右是相符的，而且与Choi^[5]和Elósegui^[7]等给出的修正后的恒星日周期也具有较好的一致性。另外，N、E、U 3个方向上相关系数最大时的相关延迟会出现不同的情况。因此，在求取序列的相关延迟时并不能单纯地考虑卫星的平均轨道周期，对3个方向上都进行同一个相关延迟，而应结合残差序列的实际情况进行具体分析。

2 恒星日滤波优化算法

由于测站在不同时段观测到的卫星各不相同，因此通过利用全部卫星的平均轨道周期以及整个恒星日序列相关性求得的相关延迟对于每一时段都不是最佳的延迟量。本文提出对残差序列进行分段，逐时段求取最大相关系数时的延迟量。文献[10]中指出，对于一天的GPS观测数据，24 h计算得到的延迟时间均不相同，且有的时段显著低于246 s，而1 h内卫星数变化相对较为平缓，因此本文算法中数据处理时将序列以1 h为单位进行滤波处理。在求取每时段最大相关系数的延迟量时，可利用该时段观测到的GPS卫星的平均轨道重复周期求得的相关延迟为中心，两边各取30个历元作为候选相关延迟量，对序列进行平移后计算相关系数，选取得到的最大相关系数时的相关延迟即为该时段最优相关延迟。具体计算流程见图 2。

3 算例分析

为了评估分段相关性优化恒星日滤波算法的有效性，采用本文提出的方法对实测数据进行滤波分析。试验数据对应为2012年002、003两天UTC为4:00~8:00间4 h的观测值，采样率为1 Hz。高频GPS数据处理软件为GAMIT/GLOBK中的TRACK运动学分析模块。为了对比优化算法的改进效率，本文分别利用3种方案对实测数据进行恒星日滤波。

方案1：传统的恒星日滤波算法，即利用相关延迟的理论值236 s对数据进行平移、滤波；

方案2：采用Choi^[5]提出的改进恒星日滤波算法，即利用32颗卫星计算得到的平均轨道重复周期求取的相关延迟对数据进行平移、滤波；

方案3：本文提出的顾及分段相关性的恒星日滤波算法。

分别利用上述3种方案进行滤波处理后，得到滤波后的坐标残差时间序列，如图 3所示，限于篇幅，仅展示N方向上的滤波效果。

从图 3可以看出，3种恒星日滤波方法均能很好地消除部分长周期误差，坐标序列的稳定性明显提高，表现出显著的随机性。相比较前两种方案，方案3的效果更好。

表 2为3种滤波处理后的坐标残差时间序列的精度统计。从表 2可以看出，3种滤波方案在一定程度上均很好地改善了坐标残差时间序列的精度。方案3由于从坐标残差时间序列本身出发，分段求取的最优相关延迟，从而更为精细地对残差序列进行滤波，取得了更优的滤波效果。在N、E、U 3个方向上，滤波后的精度，方案3较方案1、方案2分别提高了7.7%和6.1%、5.0%和3.2%、11.2%和4.8%。

为了进一步得出优化后的恒星日滤波在不同频率域上的滤波效果，本文对N、E、U 3个方向上的滤波结果进行了谱分析，如图 4所示。

由图 4可见，优化的恒星日滤波算法在0.000 1~0.01 Hz的频率上均可有效地提高定位精度，即可降低100~10 000 s周期的噪声。

4 结束语

本文提出了一种分段相关性的恒星日滤波优化算法。该算法直接以坐标残差序列为研究对象，对其进行分段求取最优相关延迟，然后进行恒星日滤波，消除多路径效应影响。通过算例验证，该算法较常规恒星日滤波算法取得更好的效果。然而，该算法仅适用于低动态的变形监测中，对于监测对象始终处于高动态的情况不适用，有待进一步研究。