半参数模型在GPS季节性信号分离中的应用

卢辰龙 赵卫林 匡翠林 刘铁柱

1 郑州市市政工程勘测设计研究院，郑州市民生路1号，450052

2 中南大学地球科学与信息物理学院，长沙市麓山南路932号，410083

GPS坐标序列不仅受构造运动的影响，还受非构造运动的影响，因此GPS坐标序列中含有速率及季节性信号，尤其是在高程方向上。研究表明，GPS时间序列不仅存在异常高频周期性信号［1］，而且周期信号的振幅也是时变的［2－3］。因此，简单地使用Nikolaidis［4］提出的模型（下文称其为传统模型）去除趋势及季节性信号，可能会忽略一些有用的周期信号，引入模型误差，从而造成速率的偏差［3，5］。此外，周期信号还会影响ITRF实现的精度［6］。目前，世界上有成千上万个连续运行测站，积累了大量的GPS数据，因此在实际应用中，有必要采用新的模型对其进行处理。

半参数模型是一种既包含参数又包含非参数分量的统计回归模型，且具有较强的解释能力。Nikolaidis提出的模型分离季节性信号会造成模型残差，而半参数模型则可以弥补其系统偏差，因此，半参数模型是一种理想的分离GPS季节性信号的模型。解算半参数模型有正则化矩阵R选取和平滑因子α确定两个关键问题。由于GPS时间序列是连续的，且相邻时刻的模型误差相差不大，可采用文献［7］的方法确定正则化矩阵R。确定平滑因子α的方法已有多种，例如广义交叉证认法［8］、L 曲线法［9］等，但这些方法都存在一些弊端，例如广义交叉证认法需事先给定平滑因子α的取值范围及其步长，对于单个时间序列并无太大影响，但若处理大量数据，其效率则至关重要；L曲线法也需事先选择一些平滑因子α，且以α为自变量的加权范数应该包含L 曲线的特征点，因此也不适用于大量数据处理。基于此，本文基于噪声时间序列Hurst指数的特性，提出一种二分搜索法确定平滑因子α。通过对模拟数据及实际的GPS坐标序列进行分析，验证该方法的有效性，尤其适用于大量GPS站时间序列季节性信号的分离。

1 基于Hurst指数确定半参数模型的平滑因子

1.1 半参数方法

半参数模型为：

其中，L为n维观测向量，X为u维参数向量，A为系数矩阵，Δ为噪声，S为信号，是非随机未知量。

方程（1）的误差方程形式为：

其对应的估计准则为：

其中，R为正则化矩阵，描述了非参数的光滑性；α为平滑因子，用于调节拟合部分VTPV与光滑部分的 平 衡。

在补偿最小二乘原则下，根据拉格朗日极值法求解，可分别得到非参数分量S与参数分量X的估值：

可以看出，解算半参数模型的关键是选取合适的正则化矩阵R与平滑因子α。在实际应用中，若L是一个观测序列，且相邻时刻的模型误差相差不大时，可令正则化矩阵取：

其中，

此时，由于正则化矩阵R秩亏（rank（R）＝n－1＜n），为确保估计准则（3）有唯一解，可增加一个约束条件［7］：

1.2 基于Hurst指数确定平滑因子

Hurst指数是一种判别时间序列是否对于时间有依赖的参数。一般来说，当时间序列的Hurst指数为0.5时，表示该时间序列是布朗运动，也即白噪声；当Hurst指数＜0.5 时，则表示该时间序列是高斯随机变量过程；当Hurst指数＞0.5时，表明该时间序列是长期记忆相关的。在实际应用中，一般将Hurst指数≤0.5的时间序列数据当作噪声［10］，若半参数模型将信号完全提取出来，那么其残差序列的Hurst指数必定等于0.5。因此，可以利用Hurst指数来获取最佳的平滑因子α。Hurst指数的计算方法有多种，本文选取受数据长度及趋势或周期影响较小的聚合方差法，其基本原理如下：

1）对于给定的时间序列Xi，i＝1，2，…，N，根据长度m，将时间序列Xi分割成［N／m］＝k个子区间；

2）计算各子区间的平均值Xm（k）；

3）计算各子区间平均值的方差var（Xm）；

4）采用最小二乘拟合双对数图上（m，var（Xm））的斜率p；

5）Hurst指数为H＝（p／2）＋1。

二分搜索法的原理如图1所示。

2 模拟实验分析

2.1 模拟数据

为验证本文提出的基于Hurst指数的二分搜索法确定半参数模型平滑因子的可靠性以及半参数模型分离GPS季节性信号的有效性，通过下面公式模拟了3 000个数据：

其中，y（ti）是时变季节性信号的模拟数据；a、b、d、e是常数；ti为GPS年积日；ε（ti）为随机噪声，服从正态分布N（0，0.25）；c（ti）是振幅变化因子，其形式如下：

图2为模拟数据，其中图2（a）、（b）分别为固定振幅的年周期、半年周期信号与加入振幅变化因子后的年周期、半年周期信号，图2（c）为周期信号叠加得到的时变季节性信号，图2（d）为加入噪声后的时变季节性信号。

图1 二分法确定平滑因子α的流程Fig.1 The flowchart of dichotomy to determine smoothing parameterα

图2 季节性信号坐标序列模拟数据Fig.2 Simulated coordination time series with seasonal signals

2.2 基于Hurst指数的二分搜索法的可行性检验

为验证基于Hurst指数的二分搜索法确定半参数模型平滑因子的可行性，首先采用不同的平滑因子α代入半参数模型，计算残差序列的Hurst指数。图3是不同的平滑因子α对应的半参数模型残差序列的Hurst指数。

图3 不同平滑因子对应的半参数模型残差序列的Hurst指数Fig.3 The Hurst index of semi－parametric model residual series for different smoothing parameter

从图3 可以看出，平滑因子α在接近0 时，Hurst指数出现递减，这是由于平滑因子α极小，造成信号被过度提取（图4），从而导致其残差序列极度平稳，使得Hurst指数在估计时呈现下降，但该变化并不影响本文提出的二分搜索法。当平滑因子α超过某个值之后，其对应的Hurst指数呈现逐渐递增，这说明采用二分搜索法是合理的，能够确定半参数模型的平滑因子α。上述实验分析证明二分搜索法是可行的，但并未能反映出二分搜索法确定的α是较优的。为此，本文又分别计算了不同平滑因子α对应模拟数据分离的SRMS与NRMS，其定义见文献［11］。

图4 不同平滑因子对应的均方根误差Fig.4 Different smooth factors corresponding to the root mean square error

从图4看出，当平滑因子α极小时，其信号被过度提取，从而使得噪声也被分离到系统误差中，使得SRMS值增大，NRMS值噪声变小；当平滑因子α极大时，系统误差并没有被充分提取，从而使得SRMS值与NRMS值都增大。从各分图的局部放大图可以看出，当平滑因子α在100左右时，其Hurst指数为0.35，对应的SRMS值最小，约为0.076，其NRMS值为0.485，与模拟噪声仅差0.015，说明此时分离得到的信号与模拟信号最接近，一般认为该平滑因子α的分离效果最好。但SRMS值是在信号已知情况下得到的，而在实际应用中，信号一般是未知的，因此，无法通过SRMS值确定平滑因子α。二分搜索法确定的平滑因子为171.56，尽管与模拟最优平滑因子α相差较大，但从图4的局部放大图中可以看出，当平滑因子α为171.56时，其对应的SRMS值为0.084，NRMS 值为0.493，与模拟最优平滑因子α对应的SRMS 值与NRMS 值相差微乎其微，与模拟信号及噪声也十分接近。

图5 模拟数据的拟合残差Fig.5 The fitting residuals of simulated data

2.3 季节性信号提取比较分析

为比较半参数模型相对于传统模型的优越性，分别采用两种模型提取季节性信号，其结果如图5所示，其中，图5（a）表示两种模型提取的信号与模拟信号的残差序列，图5（b）表示提取的信号与加噪模拟信号的残差序列。由图5 可以看出，传统模型最小二乘拟合提取季节性信号有较大的系统偏差，尤其是在振幅变化处；由基于Hurst指数二分搜索法确定光滑因子α的半参数模型则基本完整地提取出其季节性信号，避免了系统偏差残留，且模拟季节性信号的残差波动平稳，基本为0。这充分说明，在振幅变化的季节性信号中，再使用传统模型最小二乘拟合已经无法满足其实际需要，而半参数模型则能很好地解决系统偏差问题，但其平衡因子α的选取是关键。而本文提出的基于Hurst指数的二分搜索法则弥补了该缺陷，能够极大地提高其效率，从而可以推广到GPS时间序列分析应用中。

3 GPS实测坐标序列分析

本实验采用IGS 站ZIMM 的GPS 坐标序列，该站坐标时间序列不仅跨度较长，而且季节性信号呈现明显的振幅变化。图6 为ZIMM 站U方向的坐标序列半参数模型与传统参数模型的对比。从图6（b）可以看出，半参数模型充分地分离出GPS坐标序列中的趋势信号以及季节性信号，克服了传统参数模型描述实际模型不准确的缺陷。图6（b）的残差序列整体平稳，相对于图6（c）没有显著的模型残差，这也再次说明本文提出的基于Hurst指数的二分搜索法确定半参数模型的平滑因子α是有效的。传统参数模型由于没有顾及到季节性信号的时变特性，其拟合的残差序列出现了明显的系统偏差。同时，GPS坐标时间序列不仅包含年周期、半年周期信号，可能还存在其他频率的周期信号，而传统参数模型并未建模，从而可能会造成较大的系统偏差，半参数模型的非参数部分则恰好可以弥补该部分系统偏差。因此，采用半参数模型建模相对于传统参数模型具有显著的优势。

图6 ZIMM 站U 方向信号拟合及残差Fig.6 The signal fitting and residuals in U direction at ZIMM station

为进一步说明半参数模型的优势，又采用Lomb－Scargle功率谱分析法，对原始坐标序列及两种模型拟合残差序列进行分析，结果如图7所示。从图7中的竖线可以看出，半参数模型与传统参数模型都减少了年周期信号与半年周期信号的功率，但半参数模型相对于传统参数模型在低频带0.1～1.5cpy处有显著的差异，这是因为半参数模型将低频信号基本分离出来，使其残差序列减小，信噪比增大，从而使其功率谱密度减小。

图7 ZIMM 站U 方向的功率谱分析Fig.7 Power spectral density in U direction at ZIMM station

4 结 语

1）半参数模型提取周期信号的效果要优于传统参数模型。在时变周期信号的分离中，半参数模型克服了传统参数模型描述实际模型不准确的缺陷，尤其是GPS坐标序列中包含年周期信号与半年周期信号时。

2）验证了基于Hurst指数的二分搜索法确定半参数模型平滑因子的有效性与可靠性，明显提高了半参数模型数据处理的效率，为大量、自动化GPS时间序列分析提供了一种新的选择。

尽管基于Hurst指数的二分搜索法确定半参数模型平滑因子具有上述优点，但Hurst指数的计算精度受其估计方法以及时间序列跨度的影响，应予重视。

［1］田云锋.GPS坐标时间序列中的异常高频周期性噪声［J］.测绘科学，2011，36（1）：26－28（Tian Yunfeng.Anomalous High Frequency Seasonal Noises in GPS Positions Time Series［J］.Science of Surveying and Mapping，2011，36（1）：26－28）

［2］Freymueller J T.Seasonal Position Variations and Regional Reference Frame Realization［A］／／Geodetic Reference Frames［M］.Springer Berlin Heidelberg，2009

［3］卢辰龙，匡翠林，戴吾蛟，等.采用变系数回归模型提取GPS 坐标序列季节性信号［J］.大地测量与地球动力学，2014，34（5）：94－100（Lu Chenlong，Kuang Cuilin，Dai Wujiao，et al.Extracting Seasonal Signals from Continuous GPS Time Series Based on Varying－Coefficient Regression Models［J］.Journal of Geodesy and Geodynamics，2014，34（5）：94－100）

［4］Nikolaidis R.Observation of Geodetic and Seismic Deformation with the Global Positioning System［D］.San Diego：University of California，2002

［5］Blewitt G，Lavallée D.Effect of Annual Signals on Geodetic Velocity［J］.Journal of Geophysical Research，2002，107（B7）：2 145

［6］邹蓉，丁开华，杨少敏，等.利用GRACE 改正地球参考框架的季节性变化［J］.大地测量与地球动力学，2014，34（6）：50－54（Zou Rong，Ding Kaihua，Yang Shaomin，et al.Incorporating Seasonal Variations of Global and Regional ITRF Solutions Based on GRACE Measurements［J］.Journal of Geodesy and Geodynamics，2014，34（6）：50－54）

［7］孙海燕，吴云.半参数回归与模型精化［J］.武汉大学学报：信息科学版，2002，27（2）：172－174（Sun Haiyan，Wu Yun.Semiparametric Regression and Mode Refining［J］.Geomatics and Information Science of Wuhan University，2002，27（2）：172－174）

［8］胡宏昌，孙海燕.正规矩阵R及平滑因子α的选取［J］.测绘工程，2004，12（4）：5－8（Hu Hongchang，Sun Haiyan.Choice of the Regular MatrixRand Smoothing parameterα［J］.Engineering of Surveying and Mapping，2004，12（4）：5－8）

［9］王振杰，欧吉坤，曲国庆，等.用L－曲线法确定半参数模型中的平滑因子［J］.武汉大学学报：信息科学版，2004，29（7）：651－653（Wang Zhenjie，Ou Jikun，Qu Guoqing，et al.Determining the Smoothing Parameter in Semi－parametric Model Using L－Curve Method［J］.Geomatics and Information Science of Wuhan University，2004，29（7）：651－653）

［10］卢辰龙，匡翠林，易重海，等.奇异谱分析滤波法在消除GPS多路径中的应用［J］.武汉大学学报：信息科学版，2015，40（7）：924－931（Lu Chenlong，Kuang Cuilin，Yi Zhonghai，et al.Singular Spectrum Analysis Filter Method for Mitigation of GPS Multipath Error［J］.Geomatics and Information Science of Wuhan University，2015，40（7）：924－931）

［11］罗勇，匡翠林，卢辰龙，等.基于SSA 的GPS坐标序列去噪及季节信号提取［J］.大地测量与地球动力学，2015，35（3）：391－395（Luo Yong，Kuang Cuilin，Lu Chenlong，et al.GPS Coordinate Series Denoising and Seasonal Signal Extraction Based on SSA［J］.Journal of Geodesy and Geodynamics，2015，35（3）：391－395）