快速检索        
  武汉大学学报·信息科学版  2018, Vol. 43 Issue (10): 1451-1457

文章信息

马俊, 姜卫平, 邓连生, 周伯烨
MA Jun, JIANG Weiping, DENG Liansheng, ZHOU Boye
GPS坐标时间序列噪声估计及相关性分析
Estimation Method and Correlation Analysis for Noise in GPS Coordinate Time Series
武汉大学学报·信息科学版, 2018, 43(10): 1451-1457
Geomatics and Information Science of Wuhan University, 2018, 43(10): 1451-1457
http://dx.doi.org/10.13203/j.whugis20160543

文章历史

收稿日期: 2017-07-11
GPS坐标时间序列噪声估计及相关性分析
马俊1 , 姜卫平2 , 邓连生3 , 周伯烨2     
1. 武汉大学测绘学院, 湖北 武汉, 430079;
2. 武汉大学卫星导航定位技术研究中心, 湖北 武汉, 430079;
3. 湖北理工学院光谷北斗国际学院, 湖北 黄石, 435003
摘要:比较分析了极大似然估计法、最小二乘方差分量估计法以及最小范数二次无偏估计法对GPS站坐标时间序列噪声的估计效果,确定最小范数二次无偏估计法为最优的噪声方差估计方法。在此基础上对中国及其周边区域20个IGS站坐标时间序列中各方向噪声方差进行一元线性回归分析。结果表明,中国及其周边区域IGS站不同方向的噪声间具有中等强度以上的相关性,其中N方向闪烁噪声与其他方向闪烁噪声的相关性要强于白噪声。水平方向的噪声振幅能够解释垂直方向噪声振幅变化的40%~60%,而N方向噪声振幅能够解释E方向噪声振幅变化的60%~80%,获得的线性回归方程具有使用价值。
关键词GPS站坐标时间序列     最小范数二次无偏估计     噪声估计     回归分析    
Estimation Method and Correlation Analysis for Noise in GPS Coordinate Time Series
MA Jun1, JIANG Weiping2, DENG Liansheng3, ZHOU Boye2     
1. School of Geodesy and Geomatics, Wuhan University, Wuhan 430079, China;
2. GNSS Research Center, Wuhan University, Wuhan 430079, China;
3. Optics Valley BeiDou International School, Hubei Polytechnic University, Huangshi 435003, China
First author: MA Jun, PhD candidate, specializes in the noise analysis of GPS coordinates time series. E-mail: yangzhiqou.student@sina.com
Corresponding author: DENG Liansheng, PhD, lecturer. E-mail: dlseng_2000@whu.edu.cn
Foundation support: The National Science Foundation for Distinguished Young Scholars, No. 41525014; the National Natural Science Foundation of China, No. 41374033
Abstract: This paper firstly compares and analyzes the effect of maximum like lihood estimation (MLE) and least square variance estimation (LS-VCE) and minimum norm quadratic unbiased estimation (MINQUE) in the variance component estimation of GPS coordinate time series noise. After determining the minimum norm of two unbiased estimation method for the optimal estimation of noise variance, the correlations between the noises in each direction are analyzed by means of unitary linear regression, and the linear regression equations are determined.Experimental results show that, compared with the least squares variance component estimation method and maximum likelihood estimation method, MINQUE has a better estimation effect for the noise variances of GPS coordinate time series.In addition, the influence of noise on the estimation of motion parameters of the station can be reduced by using long span time series.Moreover, the same kind of noise in GPS coordinate time series of each direction has more significant correlation and the correlation between flicker noise in north and other directional flicker noises are better than that of white noises. The 40%-60% of the variance of the vertical noise can be explained by the noise variance in the horizontal direction, and north noise variance can explain the 60%-80% of change of the noise variance in east.The obtained linear regression equation has practical value.
Key words: GPS coordinate time series     least square variance estimation     noise estimation     regression analysis    

GPS坐标时间序列广泛应用于参考框架建立、地壳形变监测等研究。由于GPS坐标时间序列中的噪声严重影响了测站运动参数(包括速度、周期振幅等)的估计精度,且不同方向噪声之间存在着一定的相关性,因此需要采用有效的方法准确估计出噪声的振幅,确定不同方向噪声之间的相关性和交互影响, 进而建立严密的三维噪声模型。这有助于准确获取测站运动参数及其不确定度, 并提高GPS坐标时间序列的分析效率。

常用的GPS坐标时间序列噪声估计方法有极大似然估计法(maximum likelihood estimation,MLE)、最小二乘方差分量估计法(least square variance estimation,LS-VCE)。MLE方法通过构造一个与噪声协方差阵以及观测残差有关的对数似然函数,并寻找使函数值最大的协方差阵,可以同时估计噪声方差、周期振幅、测站速度及不确定度[1-2]。Amiri-Simkooei利用LS-VCE获得了IGS站坐标时间序列噪声的方差以及协方差,并分析了不同测站相同方向以及同一测站不同方向噪声的相关性,发现同一测站不同方向时间序列所含噪声弱相关或者不相关[3]。然而,以上两种方法须进行多次的搜索迭代过程才能获得最佳估值,运算时间较长。此外,Niu等人利用Allan方差分析法分析了采样率为1 Hz、时间长度为24 h的GPS位置时间序列的噪声特性及其振幅,取得了较好的效果[4],但该方法噪声方差估计的精度与分组的数目以及样本长度有关,对于由单天解和周解组成的采样间隔较长的时间序列精度不佳。除此之外,还有M估计、卡尔曼滤波等方法用于时间序列中方差分量估计的抗差滤波[5-6]

不同学者对以上方法研究颇多,得到了一些有用的结论,但并未对GPS坐标时间序列中噪声方差的最优估计方法形成一致意见,且忽略了各方向噪声之间的相关性。针对以上问题,本文首先比较分析了MLE、LS-VCE以及最小范数二次无偏估计法(minimum norm quadratic unbiased estimation,MINQUE)对GPS坐标时间序列噪声方差的估计效果,进而确定最优的噪声估计方法。在此基础上, 通过一元线性回归分析了中国及其周边区域20个IGS站坐标时间序列不同方向噪声方差之间的相关性及相互影响,并获得了具有实用价值的线性回归方程。

1 GPS坐标时间序列噪声估计方法

目前已有较多的文献论述MLE和LS-VCE,本文在此不再赘述,仅对MINQUE进行介绍。如果同类观测值的方差存在着不同因素的不同影响,MINQUE可用于估计这些不同因素的方差分量,因此,对于含有白噪声和有色噪声的GPS坐标时间序列,可利用该方法同时获得这两种噪声的方差,进而通过加权最小二乘准确获得GPS测站的运动参数。GPS单站、单分量的运动模型为[3]

$ \begin{array}{*{20}{c}} {y\left( {{t_i}} \right) = {x^{\left( 1 \right)}} + {x^{\left( 2 \right)}}{t_i} + \sum\limits_{k = 2}^q {\left( {{x^{\left( {2k - 1} \right)}}\cos {\omega _k}{t_i} + } \right.} }\\ {\left. {{x^{2k}}\sin {\omega _k}{t_i}} \right) + {v_i}} \end{array} $ (1)

式中,y(ti)为GPS单站单分量的坐标时间序列;x(1)x(2)分别为初始位置和速率;q为时间序列中包含的周期信号数;x(2k-1)x(2k)是调和函数的系数,用于描述周期信号的振幅;vi为误差。式(1)的矩阵形式及随机模型为:

$ \mathit{\boldsymbol{y}} = \mathit{\boldsymbol{Ax}} + \mathit{\boldsymbol{\varepsilon }} $ (2)
$ \mathit{\boldsymbol{D}}\left( \mathit{\boldsymbol{y}} \right) = {\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{y}}} = \sigma _w^2\mathit{\boldsymbol{I}} + \mathit{\boldsymbol{\sigma }}_\mathit{\boldsymbol{k}}^2{\mathit{\boldsymbol{Q}}_\mathit{\boldsymbol{k}}} $ (3)

式中,A是设计矩阵; ε是误差;D(y)为协方差矩阵;σw2σk2分别为白噪声和有色噪声的方差;I为单位阵;Qk为有色噪声的协方差矩阵。若仅考虑周年和半周年信号,则q=3,对于采样率为1/d的时间序列,矩阵Ax的具体形式如下:

$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{A}} = }\\ {\left[ {\begin{array}{*{20}{c}} 1&{{t_1}}&{\sin \left( {{\rm{2 \mathsf{ π} }}{t_1}} \right)}&{\cos \left( {{\rm{2 \mathsf{ π} }}{t_1}} \right)}&{\sin \left( {{\rm{4 \mathsf{ π} }}{t_1}} \right)}&{\cos \left( {{\rm{4 \mathsf{ π} }}{t_1}} \right)}\\ 1&{{t_2}}&{\sin \left( {{\rm{2 \mathsf{ π} }}{t_1}} \right)}&{\cos \left( {{\rm{2 \mathsf{ π} }}{t_2}} \right)}&{\sin \left( {{\rm{4 \mathsf{ π} }}{t_2}} \right)}&{\cos \left( {{\rm{4 \mathsf{ π} }}{t_2}} \right)}\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots \\ 1&{{t_n}}&{\sin \left( {2{\rm{ \mathsf{ π} }}{t_n}} \right)}&{\cos \left( {2{\rm{ \mathsf{ π} }}{t_n}} \right)}&{\sin \left( {4{\rm{ \mathsf{ π} }}{t_n}} \right)}&{\cos \left( {4{\rm{ \mathsf{ π} }}{t_n}} \right)} \end{array}} \right]} \end{array} $
$ \mathit{\boldsymbol{x}} = \left[ {{x^{\left( 1 \right)}}{x^{\left( 2 \right)}}{x^{\left( 3 \right)}}{x^{\left( 4 \right)}}{x^{\left( 5 \right)}}{x^{\left( 6 \right)}}} \right] $

其中, cos (2πti)、sin (2πti)和cos (4πti)、sin (4πti)分别为描述周年信号和半周年信号的三角函数;x(3)x(4)x(5)x(6)分别为周年和半周年项三角函数的系数。根据式(3)和分数差分的性质,可将幂律噪声看成是某一卷积函数与一白噪声卷积后得到的,因此有色噪声的协方差阵Qk可以用下式表示:

$ {\mathit{\boldsymbol{Q}}_k} = \mathit{\boldsymbol{T}}{\mathit{\boldsymbol{C}}_w}{\mathit{\boldsymbol{T}}^{\rm{T}}} $ (4)

式中,Cw为白噪声的协方差矩阵,即Cw= IT为白噪声到有色噪声的转换矩阵,其具体表达式为:

$ \mathit{\boldsymbol{T}} = \left[ {\begin{array}{*{20}{c}} {{\psi _0}}&0&0& \cdots &0\\ {{\psi _1}}&{{\psi _0}}&0& \cdots &0\\ {{\psi _2}}&{{\psi _1}}&{{\psi _0}}& \cdots &0\\ \vdots&\vdots&\vdots&\ddots&\vdots \\ {{\psi _{n - 1}}}&{{\psi _{n - 2}}}&{{\psi _{n - 3}}}& \cdots &{{\psi _0}} \end{array}} \right] $ (5)

其中,${\psi _n} = \frac{{-K/2\left( {1-K/2} \right) \cdots \left( {n-1 - K/2} \right)}}{{n!}} = \frac{{{\psi _{n - 1}}}}{n}\left( {n - 1 - K/2} \right)$ψ0=1,K为谱指数。根据式(4)、式(5)构造出有色噪声的协方差矩阵后,即可利用MINQUE估计出噪声的方差。文献[7]直接给出了MINQUE的估计公式:

$ \mathit{\boldsymbol{\hat \sigma }} = {\mathit{\boldsymbol{S}}^{ - 1}}\mathit{\boldsymbol{W}} $ (6)
$ {\mathit{\boldsymbol{S}}_{kl}} = {\rm{tr}}\left( {\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_l}\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_k}} \right) $ (7)
$ {\mathit{\boldsymbol{w}}_k} = {\mathit{\boldsymbol{y}}^{\rm{T}}}\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{Q}}_k}\mathit{\boldsymbol{Cyl}},k = 1,2 $ (8)
$ \mathit{\boldsymbol{M}} = \sum\limits_{k = 1}^2 {{\mathit{\boldsymbol{Q}}_k}} ,{\mathit{\boldsymbol{Q}}_1} = \mathit{\boldsymbol{I}},{\mathit{\boldsymbol{Q}}_2} = {\mathit{\boldsymbol{Q}}_k}, $
$ \mathit{\boldsymbol{C}} = {{\mathit{\boldsymbol{M'}}}^{ - 1}} - {\mathit{\boldsymbol{M}}^{ - 1}}\mathit{\boldsymbol{A}}\left( {{\mathit{\boldsymbol{A}}^{\rm{T}}}{\mathit{\boldsymbol{M}}^{ - 1}}\mathit{\boldsymbol{A}}} \right){\mathit{\boldsymbol{A}}^{\rm{T}}}{\mathit{\boldsymbol{M}}^{ - 1}}, $

其中$\mathit{\boldsymbol{\hat \sigma = }}\left[{\mathit{\boldsymbol{\sigma }}_w^2\;\mathit{\boldsymbol{\sigma }} _k^2} \right]$,计算前须赋予其初值$\mathit{\boldsymbol{\hat \sigma = }}\left[{1\;\;1} \right]$。从以上公式可以看出,MINQUE无需迭代或者搜寻过程即可直接获得噪声的方差分量。方差分量估值的方差为:

$ {\mathop{\rm var}} \left( {\mathit{\boldsymbol{\hat \sigma }}} \right) = 2\sigma _0^2{\rm{diag}}\left( {{\mathit{\boldsymbol{S}}^{ - 1}}} \right) $ (9)

式中,σ02为单位权方差的估值。

2 仿真数据实验结果

为比较MLE、LS-VCE以及MINQUE对GPS站坐标时间序列中噪声的估计效果,必须将实验数据中已知噪声的标准差或方差作为比较的真值。虽然有些研究机构给出了GPS站的运动参数(包括速度以及周年与半周年项的正、余弦系数),但并未提供GPS坐标时间序列噪声的估计结果,因此实测数据并不能满足本次实验的要求。本文采用仿真数据比较分析这3种方法对时间序列噪声的估计结果。

为使仿真数据能够最大程度地接近真实的GPS站坐标时间序列,本文根据BJFS、WUHN、GUAO、KUNM、SHAO以及CHAN站的运动参数模拟出包含白噪声以及闪烁噪声的时间序列。首先,在白噪声+闪烁噪声的组合模型下,利用LS-VCE对这6个测站在2010-2015年间的坐标时间序列噪声的方差进行估计,然后将估值稍微调整作为实验结果比较的真值。由于有些测站在不同方向具有相似的噪声方差估值,因此实验仅采用各测站差异较大的噪声估值。最后,根据JPL(Jet Propulsion Laboratory)提供的测站在N、E、U方向的运动参数( https://sideshow.jpl.nasa.gov/post/series.html)与调整后的噪声方差估值,通过式(1)模拟出能够反映测站运动特征的12类仿真数据:S1NS1ES1US2ES2US3ES3US4ES4US5US6ES6U。其中,S1至S6表示分别根据以上6个测站的运动参数模拟而成,下标N、E、U表示仿真数据对应的方向。每类仿真数据包含100个时间序列,每个时间序列由初始位置、线性信号、周年信号、半周年信号、零均值白噪声以及闪烁噪声叠加构成,采样点数为2 190,采样间隔为1 d。其中闪烁噪声是由白噪声通过式(5)转换获得。此外,虽然每类仿真数据下的时间序列具有相同的运动参数,但由于每次生成的随机白噪声不同,因此这些时间序列也不相同。本文分别采用以上3种方法在相同的运算环境下估计每类仿真时间序列中噪声的方差,并求其均值。各类仿真数据具体的参数设置以及噪声估计结果的均值如表 1表 2所示。

表 1 仿真时间序列噪声参数 Table 1 Parameters of Simulated Time Series

白噪声方
差/mm2
闪烁噪声
方差/mm2
S1N 0.8 1.2
S1E 1.9 1.1
S1U 4.6 11.2
S2E 1.2 2.5
S2U 10.0 23.5
S3E 0.5 1.0
S3U 3.2 8.5
S4E 1.4 1.9
S4U 21.0 12.0
S5U 12.9 7.8
S6E 0.2 2.5
S6U 4.8 12.2
注:闪烁噪声方差单位与时间相关,通常也表示为(mm·a$^{-\frac{1}{4}}$)2
表 2 各类仿真时间序列参数估计结果均值 Table 2 Parameter Estimations for Simulated Time Series
时间
序列
MLE LS-VCE MINQUE
白噪声方差/mm2 闪烁噪声方差/mm2 白噪声方差/mm2 闪烁噪声方差/mm2 白噪声方差/mm2 闪烁噪声方差/mm2
S1N 0.81±0.02 0.95±0.13 0.80±0.04 1.12±0.38 0.81±0.07 1.10±0.39
S1E 0.91±0.02 0.89±0.14 0.90±0.05 1.07±0.39 0.90±0.07 1.05±0.39
S1U 4.65±0.04 10.02±0.10 4.60±0.28 11.10±3.00 4.60±0.07 11.10±0.39
S2E 1.01±0.02 1.76±0.15 1.00±0.06 1.98±0.58 1.00±0.07 2.01±0.38
S2U 10.15±0.07 20.52±0.47 10.04±0.60 22.83±6.31 10.04±0.07 22.76±0.39
S3E 0.51±0.01 0.90±0.10 0.50±0.03 1.00±29.00 0.50±0.07 1.00±0.39
S3U 3.25±0.04 7.47±0.27 3.22±0.20 8.22±2.17 3.22±0.07 8.18±0.39
S4E 1.01±0.02 0.80±0.14 1.00±0.05 0.99±0.39 1.00±0.07 0.99±0.39
S4U 21.21±0.08 7.89±0.51 21.01±1.03 11.65±6.04 21.01±0.07 11.68±0.39
S5U 13.04±0.06 5.52±0.50 12.91±0.64 7.76±3.88 12.91±0.07 7.81±0.39
S6E 0.20±0.02 2.41±0.09 0.20±0.02 2.49±0.38 0.20±0.07 2.50±0.39
S6U 4.83±0.05 11.38±0.32 4.77±0.29 12.48±3.26 4.77±0.07 12.50±0.39
注:闪烁噪声方差单位与时间相关,通常也表示为(mm·a$^{-\frac{1}{4}}$)2

分析表 2可知,以上3种方法都能够较为准确地估计出白噪声的方差。MLE对闪烁噪声方差的估值明显小于真值,尤其是对S4U类时间序列噪声的估值与真值相差4.11 mm2。因此MLE闪烁噪声方差估值的准确度较小,而其他两种方法闪烁噪声方差估值的准确度较大。此外,绝大部分LS-VCE噪声方差估值的中误差大于其他两种方法。其中LS-VCE闪烁噪声方差估值的中误差显著大于其他两种方法,尤其是S2US4U类时间序列的中误差分别达到6.31 mm2以及6.04 mm2。另外,MINQUE对3类仿真时间序列噪声方差估值的中误差皆相等,这是由MINQUE的性质不变性所决定的。分析式(6)-(9)可知,MINQUE方差估值的中误差与设计矩阵B、观测值数n、噪声协方差矩阵Qk有关,而与观测值无关。比较3种方法的运算时间,MLE与MINQUE的运算时间几乎相等,而LS-VCE的运算时间是MINQUE的2倍。这是由于LS-VCE需要进行多次迭代过程才能获得最佳估值,因此其运行时间远大于MINQUE。

综上所述,对白噪声而言,MLE、LS-VCE以及MINQUE方差估值的准确度相当,而对于闪烁噪声而言,MLE方差估值的准确度小于其他两种方法。LS-VCE噪声方差估值尤其是闪烁噪声方差估值的精度受噪声强度的影响较大。此外,MINQUE闪烁噪声方差估值的精度以及运算效率皆高于LS-VCE,因此,MINQUE与其他两种方法相比,具有更加可靠的噪声方差估计结果和更好的运算效率。在分析数目较多的GPS站坐标时间序列时,尤其是时间跨度较长的垂直向时间序列时,宜采用MINQUE法估计噪声的方差。

3 GPS坐标时间序列噪声回归分析

研究GPS站坐标时间序列不同方向噪声之间的相关性有助于分析测站的三维运动特征,而精确地估计噪声的方差分量是分析噪声相关性的前提。基于以上研究,本文提出利用MINQUE分析由SOPAC(Scripps Orbitand Permanent Array Center)提供的中国及其周边区域20个IGS站2004-2016年内,时间跨度为4~6 a的单天GPS站坐标时间序列(ftp://garner.ucsd.edu/archive/garner/timeseries/measures/ats/Global/),在此基础上通过回归分析研究不同方向噪声的相互关系。当前国内外学者普遍认为GPS站坐标时间序列噪声特性的最优随机模型为白噪声+闪烁噪声[8],文献[9]认为滤波前中国陆态网GPS测站的主要有色噪声为闪烁噪声。虽然有些测站存在其他最佳噪声组合模型,但考虑到该类测站较少,且受样本数目所限,本文基于白噪声和闪烁噪声的组合模型估计出噪声的方差。如无特别说明,本部分内容中的“噪声”均指噪声的方差。图 1为中国及周边区域部分IGS站的白噪声和闪烁噪声估值的柱形图。

图 1 IGS站白噪声以及闪烁噪声方差估值 Figure 1 White and Flicker Noise Variances of IGS Stations

图 1可以看出,对于同类噪声,水平方向噪声估值的差异较小,而垂直方向的噪声估值远大于水平方向,此外,绝大部分测站的闪烁噪声的估值比白噪声大。造成垂直方向噪声强度远大于水平方向的原因主要有:①垂直方向卫星分布不对称,在确定平面位置时可以通过选择观测时间段及卫星来保证卫星分布的基本对称,从而消除或削弱距离测量偏差及大气延迟误差、星历误差等误差对平面位置的影响。然而对于高程测量来说, 所有被观测的卫星均处在地平面以上,卫星分布总是不对称的, 许多不对称的误差难以消除。②对流层延迟改正后的残差的影响,由于受多种因素的影响,对流层延迟改正并不完善,致使残留下来的误差主要影响高程分量的精度。为了更好地探究GPS站坐标时间序列中各方向噪声的相互关系,本文引入一元线性回归分析法对以上测站噪声估值进行分析。一元线性回归分析主要用于两个变量之间线性关系的研究,既可以用于分析和解释变量间的关系,又可以用于预测和控制。由于根据样本数据导出的线性回归方程受到抽样误差的影响,因此本文在显著性水平0.05下对回归模型进行了线性假设的显著性检验,即检验所确定的变量之间线性关系的显著性,根据模型以及自变量X估计的因变量Y的有效性。

表 3表 4分别为水平方向以及水平方向与垂直方向噪声的一元线性回归分析结果,表中fw分别表示闪烁噪声和白噪声,下标N、E以及U分别表示北方向、东方向以及垂直向。斜率以及截距皆为回归方程的参数,其中斜率为线性回归方程自变量的系数,P值是由检验统计量的样本观察值得出的原假设可被拒绝的最小显著性水平,表示拒绝原假设H0的强度,P值越小,拒绝H0的依据越强、越充分[10],即使常数项在模型中不显著,也会在模型中保留常数项,因为去掉常数项可能会给模型带来不利影响,如方程的拟合优度将会减小。相关系数用来表示自变量和因变量的线性相关程度,其绝对值越大,说明两个变量的线性关系越明显。拟合优度表征因变量的变异中有多少百分比可由自变量来解释,拟合优度越大,自变量对因变量的解释程度越高,自变量引起的变动占总变动的百分比越高,回归方程参考价值越高。

表 3 水平方向噪声一元线性回归分析结果 Table 3 Results of Linear Regression Analysis of Horizontal Direction Noises
自变量wN,因变量wE 自变量fN,因变量fE
斜率 截距 斜率 截距
模型参数 0.67 0.27 1.13 -0.04
P 1.08×10-4 0.07 9.24×10-7 0.85
相关系数 0.76 0.86
拟合优度/% 57 75
表 4 水平方向与垂直方向噪声一元线性回归分析结果 Table 4 Results of Linear Regression Analysis of Horizontal and Vertical Noises
自变量wN,因变量wU 自变量wE,因变量wU 自变量fN,因变量fU 自变量fE,因变量fU
wN 截距 wE 截距 fN 截距 fE 截距
模型参数 6.28 2.43 8.03 1.11 9.70 0.68 6.61 4.12
P 4.8×10-3 0.23 9.47×10-4 0.58 4.44×10-5 0.81 6.21×10-4 0.14
相关系数 0.60 0.68 0.78 0.70
拟合优度/% 36 46 61 49

分析表 3可知,水平方向白噪声相关系数大于0.7,闪烁噪声相关系数大于0.8,这表明水平方向同类噪声线性关系比较明显,其中白噪声具有较为显著的线性关系,闪烁噪声具有较高的线性关系。此外,水平方向自变量系数的P值远小于显著性水平0.05,因此拒绝原假设,即认为自变量系数异于0, 回归效果是显著的,所求得的水平方向噪声的线性回归方程具有实用价值。由表 3中拟合优度可知,N方向噪声能够解释E方向噪声超过50%的变化,其中N方向闪烁噪声引起的变化占E方向闪烁噪声变化的75%。

表 4可以看出,水平方向与垂直方向噪声之间的相关性介于0.6~0.8之间,具有显著的线性关系,其中,N方向与U方向闪烁噪声的相关系数接近0.8,说明它们之间的线性关系较为显著,另外,表 4中自变量的P值皆小于显著性水平0.05,因此拒绝原假设,认为回归效果是显著的,此时所求得的水平方向与垂直方向噪声的线性回归方程具有实用价值。根据表 4中拟合优度可知,水平方向噪声能够解释垂直方向噪声超过30%的变化,其中N方向闪烁噪声引起的变化占U方向闪烁噪声变化的61%,此外,水平方向闪烁噪声对垂直方向闪烁噪声变化的解释程度高于白噪声。

图 2为不同方向噪声之间的散点图及其一元线性回归直线。由图 2可知,水平方向上的回归直线拟合度较好,拟合直线的标准差较小。原因之一是水平方向噪声强度比垂直方向小很多;其次,N、E方向闪烁噪声的回归直线拟合度较高,散点能够较为密集地分布在直线周围,说明N方向与U方向的闪烁噪声有较强的相关性。

图 2 各分量噪声方差的散点图及其拟合直线 Figure 2 Scatter Plots of the Noise Variances in Different Components and Its Linear Regression

综上所述,不同方向的GPS站坐标时间序列噪声方差之间都存在中等程度以上的线性相关以及较显著的回归效果,获得的线性回归方程具有实用价值。其中水平方向噪声之间、N方向与U方向闪烁噪声之间的线性相关程度较强。N方向噪声可较大程度地解释E方向噪声的变化。此外,N方向闪烁噪声对U方向闪烁噪声变化的解释程度仅次于其对E方向的解释程度,而其他水平方向的噪声引起的变化占U方向噪声变化的百分比接近50%。

4 结语

本文首先通过仿真数据比较分析了MLE、LS-VCE以及MINQUE对GPS站坐标时间序列噪声方差的估计效果,进而确定了最优噪声估计方法。在此基础上,利用MINQUE估计中国及其周边区域的20个IGS站坐标时间序列白噪声以及闪烁噪声的方差,并采用一元线性回归分析不同方向噪声的相互关系,建立了线性回归方程并检验其有效性,得出如下结论:

1) MINQUE对GPS站坐标时间序列中噪声方差的估计精度及效率均优于MLE和LS-VCE。

2) 中国及其周边区域的IGS站坐标时间序列不同方向噪声之间都存在中等程度以上的线性相关以及较显著的回归效果,获得的一元线性回归方程可用于分析不同方向噪声之间相互依赖的定量关系,具有一定的实用价值。其中N方向噪声可较大程度地解释E方向噪声的变化,且N方向闪烁噪声对U方向闪烁噪声变化的解释程度略逊于其对E方向的解释程度。此外,其他水平方向的噪声引起的变化占U方向噪声变化的百分比接近50%。

3) 在分析测站运动特征时,可根据时间序列中N方向噪声的方差估值及与其他方向噪声方差估值的线性回归方程获得其他方向噪声的方差,从而了解整个时段内测站三维噪声强度, 进而快速、准确地估计出测站的三维运动参数及其不确定度。这对合理解释观测值的误差来源以及基准站非线性变化的物理机制具有积极的意义。

下一步研究应选择更多的GPS测站,根据构造运动强度以及测站所在区域对GPS站坐标时间序列噪声进行分类分析,构建严密的三维噪声模型。

参考文献
[1] Williams S D P. CATS:GPS Coordinate Time Series Analysis Software[J]. GPS Solutions, 2008, 12(2): 147–153 DOI:10.1007/s10291-007-0086-4
[2] Williams S D P. The Effect of Coloured Noise on the Uncertainties of Rates Estimated from Geodetic Time Series[J]. Journal of Geodesy, 2003, 76(9): 483–494
[3] Amiri-Simkooei A R. Noise in Multivariate GPS Position Time-Series[J]. Journal of Geodesy, 2009, 83(2): 175–187 DOI:10.1007/s00190-008-0251-8
[4] Niu X, Chen Q, Zhang Q, et al. Using Allan Variance to Analyze the Error Characteristics of GNSS Positioning[J]. GPS Solutions, 2014, 18(2): 231–242 DOI:10.1007/s10291-013-0324-x
[5] Tong Haibo, Zhang Guozhu. Robust Positioning Algorithm with Modified M-estimation for Multiple Outliers[J]. Acta Geodaetica et Cartographica Sinica, 2014, 43(4): 366–371 ( 仝海波, 张国柱. 改进M估计的抗多个粗差定位解算方法[J]. 测绘学报, 2014, 43(4): 366–371. )
[6] Zhao Changsheng, Tao Benzao. Kalman Filtering of Linear System with Colored Noises[J]. Geomatics and Information Science of Wuhan University, 2008, 33(2): 180–182, 207 ( 赵长胜, 陶本藻. 有色噪声作用下的卡尔曼滤波[J]. 武汉大学学报·信息科学版, 2008, 33(2): 180–182, 207. )
[7] Cui Xizhang, Yu Zongchou, Tao Benzao, et al. Ge-neral Surveying Adjustment[M]. Wuhan: Wuhan University Press, 2009 ( 崔希璋, 於宗俦, 陶本藻, 等. 广义测量平差[M]. 武汉: 武汉大学出版社, 2009 )
[8] Mao A, Harrison C G A, Dixon T H. Noise in GPS Coordinate Time Series[J]. Journal of Geophysical Research:Solid Earth, 1999, 104(B2): 2797–2816 DOI:10.1029/1998JB900033
[9] Wang W, Zhao B, Wang Q, et al. Noise Analysis of Continuous GPS Coordinate Time Series for CMONOC[J]. Advances in Space Research, 2012, 49(5): 943–956 DOI:10.1016/j.asr.2011.11.032
[10] Sheng Zhou, Xie Shiqian, Pan Chengyi. Probability Theory and Mathematical Statistics[M]. Beijing: Higher Education Press, 2008 ( 盛骤, 谢式千, 潘承毅. 概率论与数理统计[M]. 北京: 高等教育出版社, 2008 )