香港CORS站坐标时间序列分析研究

2019-05-07陈晨魏冠军高志钰寇瑞雄

全球定位系统 2019年2期

陈晨,魏冠军,高志钰,寇瑞雄

(1.兰州交通大学测绘与地理信息学院,甘肃兰州 730070;2.地理国情监测技术应用国家地方联合工程研究中心，甘肃兰州 730070；3.甘肃省地理国情监测工程实验室,甘肃兰州 730070)

0 引言

全球定位系统(GPS)基准站坐标时间序列既可以反映出测站的线性变化,又可以反映出测站存在的非线性变化[1]．建立合适的区域连续运行参考站(CORS)噪声模型对研究CORS站的稳定性及区域地表形变规律具有重要的意义．Zhang等[2]首次将噪声分析引入到GPS坐标时间序列分析中来,研究结果表明坐标时间序列中具有“白噪声(WN)+闪烁噪声(FN)”的特征．Mao等[3]采用谱分析和极大似然估计方法分析了23个全球IGS站的噪声特性,得出最佳噪声模型为WN加FN．Langbein[4]分析了加利福尼亚南部及内华达南部地区236个连续GPS运行站的噪声模型,认为最佳噪声模型以FN或者随机游走噪声RW为主．蒋志浩等[5]对我国国家CORS网进行了数据处理,认为国家CORS网基准站坐标时间序列噪声表现为WN、FN及RW．李昭等[6]对中国区域11个IGS基准站坐标时间序列进行地表质量负荷改正后,得到各坐标分量上的最优噪声模型主要以WN加FN为主,同时也存在幂律噪声(PL)、带通噪声(BP)及一阶高斯马尔科夫噪声(FOGM)．李斐等[7]对南极半岛地区8个IGS基准站坐标时间序列进行空间滤波后认为各站点坐标时间序列中不仅存在WN,还存在有色噪声(CN),在有些站点的E方向可能还存在WN．袁林果等[8]对香港12个GPS连续运行参考站坐标时间序列分析后,认为最优噪声模型为可变自噪声(VW)加FN．

由此可见,不同地区噪声特性存在着差异,针对此项问题,本文以香港卫星参考站网(SatRef)为研究对象,来分析测站坐标时间序列的周期特性与噪声特性．在袁林果研究的基础上,选择更多的噪声模型来分析确定最优噪声模型．考虑到香港地区区域较小(1106.34平方千米),最长基线HKNP-HKWS仅为49.89 km,参考站之间具有很强的空间相关性．因此,本文决定选取石碑山(HKOH)、昂坪(HKNP)、锦田(HKKT)、蓝地(HKLT)和南丫岛(HKLM)这五个数据质量较好的站点来进行分析．利用功率谱分析的方法来获取各个参考站点的周期项．对单天解坐标分量时间序列建立拟合参数模型,采用6种噪声类型,将其组合成14种噪声模型,运用极大似然法估计参考站坐标时间序列在各个噪声模型下的噪声量级、周期性振幅、测站速度及不确定度等信息．选用Langbein提出的保守估计准则判断不同噪声模型的优劣,确定测站各分量上的最优噪声模型．比较分析在仅考虑WN和在顾及CN的影响下,参考站的线性速度以及周期振幅和相位等参数估计值的差异性,建立合适的区域CORS网速度场模型．

1 区域CORS网数据预处理

香港卫星参考站网(SatRef)截止到2018年10月1日共有18个参考站,其中T430站和大埔滘(HKTP)站的观测时间较短,所以本次数据处理不采用这两个点．选取数据长度为2014-05-01到2017-08-31共40个月的参考站网观测数据,同时引入ITRF2014参考框架下BIFS、CUSV、LHAZ、PIMO、SHAO、TWTF、TCMS和PBRI共8个IGS站的同步观测数据,起到一个约束的作用．基线解算时海潮改正模型采用FES2004模型,对流层模型采用Saastamoinen改正模型,地潮改正模型采用IERS03．检验基线解算的成果,其最短、最长基线重复率分别为0.4 mm和1.9 mm,最短基线的相对重复率为2.24×10-8量级,最长基线相对重复率为5.87×10-10量级．基线解得标准均方根误差(NRMS)值大多处于0.15～0.18之间,最大不超过0.2,表明基线解算精度很高,符合网平差要求．GAMIT基线解算完成后,利用GLOBK软件进行平差,得到各参考站N、E、U三个方向上的平均NRMS值分别为1.9 mm、1.8 mm和1.3 mm,均优于3 mm,达到了较高的解算精度．

经过GAMIT/GLOBK软件处理后得到各站的原始坐标时间序列,由于篇幅限制,此处只列出 HKOH站三个坐标分量上的原始坐标时间序列图,并对其进行了简单的线性拟合,如图1～3所示．

图1 HKOH站N方向原始坐标时间序列图

图2 HKOH站E方向原始坐标时间序列图

图3 HKOH站U方向原始坐标时间序列图

从图1～3可以得出以下结论:

1)HKOH站在N、E两个水平方向主要表现为线性变化,且有整体向东南方向运动的趋势,N方向的运行速度约为-12 mm/a,E方向上的运行速度约为32 mm/a．

2)U方向上主要呈现为周期性变化,周期性振幅具有夏季大,冬季小的季节性特征,振幅约为10 mm左右．

3)该站在U方向上要比在N、E方向上的离散程度大,这是因为高程方向上的精度要低于水平方向精度,即在水平方向上该点较稳定,在高程方向上波动较大．

2 参考站坐标时间序列功率谱分析

功率谱分析法是分析时间序列的常用方法之一．功率谱密度表示了各个频率上的能量强度.观察离散数据系列的周期性,如果某个信号具有周期特性,则其周期运动所对应的功率在全部功率中占有比重较大,在功率谱图上表现为峰值[9-10]．

地球物理现象的功率谱可以用下式表示:

P(f)=P0f-α,

(1)

其对数形式为

InP(f)=InP0-αInf,

(2)

式中:P(f)为功率谱密度;f为频率;P0为系数;α为谱指数,P0和α都为待定参数．

许多研究表明GPS坐标时间序列中存在着周期信号,表现为年周期与半年周期特性[11-13]．本文选择了SatRef中HKOH、HKKT、HKNP、HKLM和HKLT五个具有代表性的站点来探测香港地区的卫星参考站的坐标时间序列是否也具有周期特性．主要方法是利用快速傅里叶变换(FFT)来获取各个站点的功率谱密度P(f),通过分析它们的功率谱图,能够直观地揭示各参考站点N、E、U三个坐标分量上的周期特性．需要注意的是,在进行FFT的时候,原始坐标时间序列需要满足均匀采样和零均值的特性．在不满足上述条件的情况下,也可以利用周期图法进行谱分析,弊端是计算时间较长[14]．

去除掉原始坐标时间序列线性趋势,经过FFT得到各参考站的功率谱图,如图4～6所示．

图4 各基准站N方向功率谱图

图5 各基准站E方向功率谱图

图6 各基准站U方向功率谱图

图4～6中横坐标表示年周期数,纵坐标表示单位周期内的能量．从图中可以直观地看出各个参考站点在N、E、U三个坐标分量方向上均表现出一定的周期特性,且以年周期项最为明显．各个分量上还存在较为明显的半年周期项．HKNP和HKLT两个站的U方向上还存在1.5年周期项．经过统计发现,所有参考站点在U分量方向上的功率谱能量高于N、E分量方向2～5倍,这说明高程方向上的周期性变化要比水平方向变化更明显．经过统计发现,除HKOH站以外,其余各站N分量方向上年周期所对应的功率谱能量普遍大于E方向,说明在N方向的年周期变化要比E方向明显．

3 参考站坐标时间序列噪声模型建立与分析

3.1 极大似然估计法

GPS坐标时间序列噪声分析可以用极大似然估计(MLE)的方法来完成．MLE可以同时估计噪声类型、周期性振幅、测站速率及不确定度.

MLE可以估计参考站坐标时间序列残差中包含的WN与相关噪声的振幅．主要原理是对单天解坐标分量时间序列建立如下的参数模型[6,15]:

y(ti)=a+bti+csin(2πti)+dcos(2πti)+

esin(4πti)+fcos(4πti)+

(3)

NIn(2π)],

(4)

式中:N为时间序列长度;C为协方差阵;det表示矩阵C的行列式;lik为似然值．

3.2 噪声类型

大量的研究成果表明,GPS坐标时间序列中同时存在WN和CN．本文主要分析WN与CN的组合噪声模型．在实验中选取的6种噪声分别为:WN、VW、FN、RW、PL和一阶高斯马尔科夫噪声(GM)．将6种噪声分成两组组合成14种噪声模型,第一组噪声模型分别为:WN、WN+FN、WN+RW、WN+FN+RW、WN+GM、WN+RW+GM和WN+PL．第二组噪声模型分别为:VW、VW+FN、VW+RW、VW+FN+RW、VW+GM、VW+RW+GM和VW+PL．

3.3 最优噪声模型评价准则

根据极大似然估计法的原理可知,不同的噪声模型将得到不同的MLE值,一般认为,极大似然值越大,结果越可靠[3-4]．但是噪声模型中包含的未知参数越多,其MLE值也越大．因此,为了确保结果的可靠性,不能只以MLE值的大小来判定最优噪声模型．对于不同的噪声模型,蒙特·卡洛实验表明:95%的显著水平下,当两种噪声模型的MLE之差大于3.0时,两种模型具有可区分性[15]．

3.4 香港卫星参考站坐标时间序列噪声分析

3.4.1 WN模型的选取

针对SatRef中选取的五个参考站点,运用极大似然估计法得到其在14种噪声模型下的MLE值[16]．对五个站点N、E、U三个分量上在WN模型和VW模型下的平均MLE值进行比较,结果如表1所示．

表1 WN模型与VW模型MLE值对比

从表1中可以看出,简单WN模型的MLE值大于VW模型的MLE值,根据最优噪声模型评价准则可知,就WN模型而言,WN模型比VW模型更适合于SatRef．此外,利用MLE计算得到的五个参考站分量VW+FN模型的MLE值明显大于其他组合模型,这与袁林果等对SatRef的噪声模型研究结果一致[8]．但此后的研究结果认为VW模型仅能反映测站分量的质量好坏,并不能作为参考站最优噪声模型[6]．根据上述两点原因,本文决定采用WN与其组合模型来分析SatRef的最优噪声模型．

3.4.2 最优噪声组合模型的确定

将WN+FN、WN+RW、WN+FN+RW、WN+GM、WN+RW+GM和WN+PL六种噪声模型的MLE值与WN的MLE值作差,得到MLE值的差值图如图7所示.

从图7中可以看出:6种噪声组合模型与WN模型MLE值的差值都大于0,表明SatRef的坐标时间序列中不仅含有WN,还存在CN．WN+FN、WN+GM、WN+RW+GM和WN+PL四种噪声模型与WN模型的差值大于其他两种噪声模型,说明这4种噪声模型较其它模型更优越,且这四种噪声模型的MLE值相差不大．

通过3.3节的介绍可知,不能只以MLE值的大小来判定参考站点最优噪声模型,为了确保结果的可靠性,本文选用Langbein提出的保守估计准则判断不同噪声模型的优劣．首先,选取WN+FN和WN+RW两种组合模型MLE值较大者作为零假设,然后将WN+FN+RW、WN+GM和WN+PL噪声模型的MLE值与零假设作比较,若MLE差值大于2.6,则拒绝零假设,认为该模型更优,否则接受零假设,认为所选的复杂模型无效．若WN+FN+RW、WN+GM和WN+PL均优于零假设,则选择MLE较大的作为最优模型,根据此方法来寻找各个参考站点N、E、U分量上的最优噪声模型．本文将接受WN+RW+GM组合模型的阈值设为5.2,接受其余噪声组合模型的阈值设为2.6．表2示出参考站各分量的最优噪声模型．

图7 五个参考站N/E/U方向MLE值差值

参考站点N分量E分量U分量 HKOHWN+FNWN+FNWN+PL HKKTWN+FNWN+FNWN+GM HKNPWN+GMWN+GMWN+PL HKLMWN+PLWN+FNWN+FN HKLTWN+GMWN+FNWN+FN

从表2中可以看出SatRef坐标时间序列的噪声特性呈现多样化的特征,15个坐标分量上WN+FN组合模型所占比例最大,为53.3%．WN+GM和WN+PL组合模型所占比例分别为26.7%和20%．由此可见,尽管香港地区区域较小,参考站点之间相距较近,但不同的站点的坐标时间序列仍表现出不同的噪声特性．目前对产生这种差异性的原因尚不明确,推测可能是与框架的定义以及地球物理效应的空间相关性有关,有的测站可能存在较大的系统误差[5-6]．

4 顾及有色噪声对参数估计的影响

利用CATS软件进行极大似然估计时,会得出公式(3)中的各个参数．为顾及CN对参数估计的影响,在ITRF2014框架下,本文对SatRef中选取的五个参考站点采用WN+FN噪声模型估计其线性速率、年周期振幅和相位三个参数,并与WN模型进行比较[7,17-18]．

表3所示为参数估计结果,可以得出以下两个结论:

1)在WN+FN和WN两种不同的噪声模型下所估计的参数具有明显的差异．三个参数中,速率最大差异在HKLM站的N方向,为2.64 mm/a;年周期振幅最大差异在HKNP站的N方向,为0.54 mm;相位最大差异在HKLT站的U方向,为2.18 rad．

2)比较在WN+FN和WN两种噪声模型下所得到的线性速率、年周期振幅和相位三个参数估计的最大不确定度,前者约为后者的10倍、5倍和8倍．由此说明,仅考虑WN所获取的参考站线性速率精度并不能反映实际的SatRef速度场的精度,同时也在很大程度上低估了年周期振幅和相位的参数估计不确定度．

结合结论1)、2)可知,在进行参考站坐标时间序列分析时,CN对参数估计的影响是不容忽视的．

5 顾及CN的水平速度场分析

由第4章可以得出在WN+FN噪声模型下HKOH、HKKT、HKNP、HKLM和HKLT五个站点的线性速率估值及其精度,其中N方向的平均速率约为-11.46 mm/a,E方向的平均速率约为31.45 mm/a,U方向的平均速率为-0.18 mm/a,水平速度场图如图8所示．由文献[19]和[20]可知位于华南块体的SHAO、WUHNS、XIAM的平均速度分别为为35.50 mm/a、35.82 mm/a和34.54 mm/a[19-20],结合SatRef的水平速度为33.47±1.5 mm/a可知,华南块体整体稳定,运动速度与趋势基本一致．