基于CATS软件的GPS时间序列分析

叶锦强

（1.惠州市国土资源局 大亚湾经济技术开发区分局，广东 惠州 516081）

GPS已被广泛运用到地壳形变的研究上，而将GPS资料制成时间序列后，可以分析GPS的连续观测资料得到测站点与时间及空间的关系，进而推估出各种构造与地壳形变之间的关联性。

GPS的观测数据中存在噪声。数据噪声可以简单地分为2种：与时间无关的白噪声和与时间相关的有色噪声。白噪声是一种随机的信号，在功率频谱上呈现水平的功率谱密度。白噪声在光功率频谱内位于可见光波段，当眼睛接受光的三原色分量相等时为白光。白噪声之外的噪声统称为有色噪声。

在多数的GPS时间序列的研究中，仅将观测资料的误差视为与观测时间无关的白噪声，因为白噪声的数值模型和计算比较容易。但大量研究表明，如果GPS时间序列分析中不考虑与时间相关的有色噪声，将会低估地壳形变速度的误差[1-3]。尤其是在一些地壳活动速度比较低的地区测得的GPS位移量很小，若低估误差可能导致地壳变形的大小及方向上的错误，这时候考虑有色噪声的影响就格外重要。本文中将有色噪声分为闪烁噪声和随机游走噪声。前者可能来自GPS系统中未被模型化的大气效应、卫星轨道及星历误差、多路径效应、卫星海拔高度及方位角所造成的效应[3-5]；后者来自于地表潜移或降雨使仪器支柱产生布朗运动所造成的效应[6]。

本文利用时间序列分析软件CATS[7]，基于两种不同的方法来进行资料的噪声分析，分别是时间域上的最大似然估计法和频率域上的频谱分析法。

1 GPS资料处理与噪声分析方法

对于GPS的时间序列的分析，通常考虑长周期运动、同震变形以及震后变形。GPS单站、单分量坐标序列一般用下列周期性模型来表示[8]：

式中，ti以年为单位；系数a、b分别代表地壳位置和线性变换率；系数c、d、e和f描述了周年和半周年运动振幅；系数gj代表地震造成的同震偏移（offset）；H（t）为阶跃函数（heaviside step function）；vi为残差值，代表观测值与预测值间的差异，也就是本文需要分析的噪声信号。

由于研究目的是与地壳形变无关的噪声信息，所以针对GPS时间序列都必须按照式（1）事先扣除已知的信号，包括变形速率、各种原因造成的偏移量、季节性的周期变化（在此为年周期和半年周期），研究区站点分布见图1。

图1 研究地区站点分布图

图2 SUIY测站时间序列处理事例

图2a中从上至下分别是南北分量、东西分量、垂直分量的原始时间序列（红色线）及其拟合线（绿色线），图2b中从上至下分别是南北分量、东西分量、垂直分量的残差时间序列，也就是本文所要研究的噪声时间序列。

1.1 频谱分析

GPS原始时间序列扣除拟合值后的噪声时间序列，可以用频谱定性分析噪声的类型。频谱分析是一种在频率域上分析信号的方法。噪声时间序列在频谱域上的功率谱可以使用Power Law的形式表示[6]：

式中，α为频谱指数；P0为常数。α越大，表示噪声序列的时间相关性越高。

将式（2）取对数进行曲线拟合，运用最小二乘法可以得到P0及α值。α通常是介于-1～3之间的任意实数，其中整数α代表一些特殊的噪声类型：α=0时，是标准的白噪声；α=1时，是标准的闪烁噪声；α=2时，则是标准的随机游走噪声。另外，除了标准的白噪声以外，其余的部分统称为有色噪声。

1.2 最大似然估计法

最大似然法是一种非线性的最小二乘法，目的在于找到与时间序列最相近的模型参数。通过调整协方差矩阵使得似然函数取得最大值，即可得到与该时间序列最相近的噪声模型，这样就可以通过协方差矩阵估算出时间序列中的噪声振幅大小。与频谱分析相比，最大似然估计法可以定量地计算噪声大小。

最大似然估计法的方程式如下[1]：

或

式中，ln是自然对数；lik（x，C）是概略值；n是资料点数；det是行列式；C是资料的协方差矩阵。确定协方差矩阵的方法可以参考文献 [1]、[2]、[9]、[10]。

2 时间序列分析软件CATS

CATS软件包含2个程序：线性部分使用最小二乘，利用式（1）计算截距、斜率、偏移以及周期信号的振幅；非线性部分则利用线性计算后的残差求得特定噪声模型的参数及振幅大小。

假设模型中只有一种噪声，则模型的协方差矩阵C=σ2J（σ表示噪声振幅大小），则式（4）可以改写为：

此时，可以由式（5）得到噪声振幅的关系式，此关系式可以令式（5）达到最大值：

当模型中不只一种噪声模型时，CATS用下面的方法来计算最大似然估计值。设2个噪声振幅分别σ1和σ2，将振幅以角度φ和一个比例常数r表示：

而协方差矩阵可以表示成：

代入式（4），就可以简单地运用角度φ计算最大似然估计法的数值。

3 最佳噪声模型分析

本文使用3个噪声模型来分析噪声时间序列：白噪声（WN）、白噪声+闪烁噪声（FL）、白噪声+随机游走噪声（RW），利用CATS软件，使用最大似然估计法在时间域上求最佳噪声模型，并使用频谱分析法从频率域上找出最佳模型。

表1中针对3种模型计算出各测站噪声三分量的MLE值。从表中可以统计出，16个GPS测站的48个分量中有90%以上的FL有较大的MLE值（相比RW和WN的MLE值）。只有极少分量（JLYJ的E分量和SUIY的N、E分量）的RW有较大MLE值，这可能是因为观测或解算中存在误差，也不排除与该站的基座有关。虽然存在上述异常情况，但整体上可以说明，白噪声+闪烁噪声是描述研究地区GPS观测网噪声时间序列的较好模型。

表1 GPS观测网中各测站的MLE值

另外需要注意的是，“陆态网络”的GPS仪器基座为水泥柱，属于深基及浅基支柱型的较稳定的仪器基座，相对来说不容易产生随机游走噪声，这与计算的结果也是一致的。

利用CATS软件计算出各测站三分量的有色噪声频谱指数，将频谱分析的模型设定为噪声加上一个Power Law的有色噪声，计算结果如表2。从表中可以看出，研究地区的GPS资料的频谱指数基本介于0～-1之间，与前人的研究结果相同[2,3]，与最大似然估计法的计算结果也相符。

表2 观测网中各测站的频谱指数

4 结 语

本文对GPS观测数据进行噪声分析。基于最大似然估计法计算结果表明，白噪声+闪烁噪声是较近似的噪声模型。少数几个噪声接近白噪声+随机游走噪声模型的测站，可能是由于仪器基座不稳定或者时间序列长度较长所产生的效应。而频谱分析所得到的结果与最大似然估计法一致。不同的噪声模型将对地壳形变速度的误差造成很大的影响，而研究更精确的模型对于反演地球自身的运动状态尤为重要。

[1]Langbein J,Johnson H.Correlated Errors in Geodetic Timeseries:Implications for Time-dependent Deformation[J].Geophys Res,1997,102(1):591-603

[2]Zhang J,Bock Y,Johnson H,et al.Southern California Permanent GPS Geodetic Array:Error Analysis of Daily Position Estimates and Site Velocities[J].Geophys Res,1997,102(8): 35-55

[3]Mao A,Harrison C G A,Dixon T H.Noise in GPS Coordinate Time Series[J].Geophys Res, 1999,104(1):2 797-2 816

[4]Williams S,Bock Y,Fang P,et al.Error Analysis of Continuous GPS Position Time Series[J].Geophys Res,2004,109(3):412-421

[5]Langbein J.Noise in Two-color Electronic Distance Meter Measurements Revisited[J].Geophys Res,2004,109(4):1-16

[6]Agnew D C.The Time-domain Behavior of Power Law Noises[J].Geophysical Research Letters,1992,19(4):333-336

[7]Williams S D P.The Effect of Coloured Noise on the Uncertainties of Rates Estimated from Geodetic Time Series[J].Geodesy, 2003,76(10):483-494

[8]Nikolaidis R.Observation of Geodetic and Deismic Deformation with the Global Positioning System[D].University of California,2002

[9]黄立人.GPS基准站坐标分量时间序列的噪声特性分析[J].大地测量与地球动力学,2006,26(2):31-34

[10]Gardner M.White and Brown Music, Fractal Curves and Oneover-fluctuations[J].Scientific American, 1978,238(4):16-32