2012, Vol. 55
2. 武汉大学测绘学院,武汉 430079;
3. 中国测绘科学研究院,北京 100039;
4. 导航与位置服务国家测绘地理信息局重点实验室,北京 100048
2. School of Geodesy and Geomatics, Wuhan University, Wuhan 430079, China;
3. Chinese Academy of Surveying and Mapping, Beijing 100039, China;
4. Key Laboratory of Navigation & Location Based Service of NASG,Beijing 100048, China
我国于2008年7月1日启用了中国地心坐标系统2000(CGCS2000)和相应的坐标参考框架[1-3]，国家CORS(全球导航卫星系统连续运行参考站)站运动特征分析是维持CGCS2000(中国大地坐标系统2000)框架精确性、动态性和稳定性的重要基础.目前在分析CGCS2000 下我国CORS 站运动特征时，只考虑了地壳构造运动引起的线性运动速率，没有顾及由非构造形变等因素引起的CORS 站的非线性运动，特别是在CORS高度方向的周期性变化(年振幅可达1~2cm)，导致国家CORS站的速度估值和误差估计存在偏差，降低了CGCS2000 坐标框架精确性、动态性，也使得目前无法实现周期性稳定的CGCS2000坐标参考框架，影响了地壳运动监测和相关地学研究成果的可靠性.
国家CORS站的非线性运动主要受到地壳非构造形变信息、框架点公共误差信息和观测误差信息的组合影响.随着CGCS2000下国家CORS网坐标时间序列资料的不断积累，使得研究国家CORS站的非线性运动成为可能，特别是分析站点高度方向的周期性运动成为可能.本文研究基于计算获取的CGCS2000下我国国家CORS网1999年至2009年坐标时间序列，在分析其中蕴藏的地壳非构造形变信息、框架点公共误差信息和包含的噪声性质基础上，精确描述了CGCS2000下我国CORS站运动特征.研究方法首先采用国际卫星对地观测数据及相关地球物理模型，分别计算了大气负荷、海洋非潮汐、积雪和土壤水等4 项质量负荷效应造成的地壳非构造形变，并以此分析和模制了这些非构造形变对国家CORS站坐标时间序列的影响.其次采用主成分空间滤波算法(PCA)提取质量负荷改正后的国家CORS 网坐标时间序列中包含的坐标框架公共误差(CME)的时空特性.最后采用最大似然法和加权最小二乘法评定了考虑不同噪声影响的CGCS2000框架下的国家CORS网年速度估值和实际精度.
2 国家CORS站数据处理本文对国家CORS网的数据处理采用3个步骤.
第1步，利用GAMIT 软件获得国家CORS 网点坐标和卫星轨道的单日松弛解；
第2 步，利用GLOBK 软件将获得的国家CORS网单日松弛解结果和全球IGS连续运行站的松弛解合并，得到一个包含国内CORS网和全球IGS(InternationalGNSSService)站的单日松弛解结果；
第3步，利用GLOBK 软件平差处理上一步合并后的单日松弛解，并选取全球IGS站在ITRF1997框架的坐标和速度值作为约束(选取方法见参考文献[4])，采用相似变换方法求解出国家CORS网成果在CGCS2000坐标框架中的单日解坐标值.
本文作者采用上述方法对近10年CORS站数据进行处理，获得了国家CORS 站29 站的坐标时间序列，CORS站坐标时间序列水平分量的重复性优于5mm, 垂直分量的重复性优于12mm, 结果达到了较高的精度，为开展CORS站运动速度估值和精度估计奠定了较好的基础.
3 地壳非构造形变分析目前估计国家CORS站运动速度参数时，通常由CORS站日解坐标值组成CORS 站原始坐标时间序列，利用纯数学意义上的回归分析方法建立CORS站坐标变化的经验模型，由于采用谐波分析方法不能够完整提取非构造形变信息，不利于分离坐标时间序列中包含的地球质心周期性运动和GNSS技术系统偏差的影响.因此分析地球物理因素对国家CORS站非构造形变的影响，并对其影响进行定量改正是分离非构造形变影响的根本途径，是维持CGCS2000坐标框架的精确性、动态性的必要保证[5-6].
引起地壳(表)非构造形变的地球物理因素主要包括两大类[7].第一类是潮汐形变，包括固体潮、海洋潮和极潮.这类形变目前已经建立了较为精确的计算模型，在CORS 站数据处理时做了相应的改正.固体潮和极潮改正采用IERS2003标准模型，海洋潮模型采用进行了地球质心改正的FES2004 模型.第二类是地球表面流体圈中的大气和各态水的质量迁移引起的地表质量负荷变化，主要包括大气、非潮汐海洋、积雪和土壤水等质量负荷.地表质量负荷效应引起的地表形变在CORS 站数据处理中未作改正，是国家CORS站时间序列中非线性变化部分的主要来源，特别是高度方向表现出周期性变化的主要原因.
本文采用负荷格林函数方法分别计算了大气负荷、海洋非潮汐、积雪和土壤水等4项质量负荷效应造成的地壳非构造形变，并以此分析和模制了这些非构造形变对国家CORS站坐标时间序列的影响.文中大气负载的计算取自美国国家环境预报中心(NCEP)提供的采样率为6小时的格网地表气压资料(下载地址为:ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis/)；积雪和土壤湿度资料取自美国国家海洋与大气局的日平均解数据(下载地址为:ftp://ftp.cdc.noaa.gov/Datasets/ncep.reanalysis.dailyavgs)；非潮汐海洋资料取自美国国家航空航天局(NASA)下喷气动力实验室(JPL)的ECCO 数据(http://ecco.jpl.nasa.gov/external/rsync.php).
地表质量负荷主要作用于径向，主要影响CORS站的垂向位置变化.由于篇幅所限，本文只给出了具有区域代表性的北京房山(BJFS)、海口(HAIK)、哈尔滨(HRBN)、拉萨(LHAS)、武汉(WUHN)、乌鲁木齐(URUM)6个站的高度分量质量负荷改正时间序列(见图 1).
|
图 1 BJFS、HAIK、HRBN、LHAS、WUHN、URUM高度分量时间序列的负荷改正 S1表示国家CORS站髙度分量原始时间序列，S1中红线代表对S1的谐波拟合，32、33、34、35分别表示大气、土壤水、积雪和海洋非潮汐负荷变化引起的髙度分量变化，S6表示负荷改正后的坐标时间序列，S6中蓝线代表对S6的谐波拟合. Fig. 1 Loading corrected for coordinate time series in the height direction: BJFS，HAIK and HRBN S1 denote the original coordinates time series, S2» S3，S4» S5 denote the position perturbations caused by the atmospheric, soil moisture, snow and non-tidal ocean mass loading respectively.Red curves and blue curves overlaying on S1 and S6 are annual and semi-annualcorrections obtained from empirical fitting. |
我国CGCS2000参考框架是在全球ITRF1997参考框架基础上实现的.全球ITRF1997 参考框架存在非线性变化，这样就不可避免地导致CGCS2000参考框架也存在着非线性的变化，因此我国CORS站在CGCS2000 框架下的原始坐标时间序列中包含时空相关的公共误差，公共误差是国家CORS网数据处理的主要的误差源.公共误差的存在使得我国CORS 站的公共变化特征掩盖了局部变化特征，影响了CGCS2000 框架下我国CORS站速度估值和误差估计，本文作者采用主成分空间滤波算法(PCA)去除了我国CORS网公共误差，定量分析了由公共误差引起的CORS站非线性运动.
主成分分析(PrincipalComponent Analysis, PCA)是研究如何将多指标问题转化为较少的综合指标的一种重要的统计方法，它是将高维空间的问题转化到低维空间去处理，使问题变得比较简单、直观，而且这些较少的综合指标之间互不相关，能提供原有指标的绝大部分信息[8].对于一个共有n站，观测了m天的国家CORS 网时间序列，矩阵Xij=(x1，x2，…，xj)(其中i= 1，2，…，m；j= 1，2，…，n)中的每一列包含着站点某一个方向的去趋势化和去均值化后的值，每一行表示在一个给定的历元，所有站点的某个分量的值.X的协方差矩阵B定义为
|
(1) |
它是(n×n)的对称矩阵，可以分解为B=VΛVT，式中，特征向量构成的矩阵VT 是一个(n×n)的正交矩阵，Λ 矩阵是由k个非零对角元素λ 组成的特征值矩阵.在多数情况下，矩阵B为满秩矩阵，即k=n，所以矩阵B的特征值-特征向量对可以写成:(λ1，v1)，(λ2，v2)，…，(λn，vn).针对Xij的任一行向量，第k个主成分表示为
|
(2) |
式中vjk表示第k个特征向量的第j个分量.同时矩阵Xij可以由主成分表示为
|
(3) |
该方法中主成分代表共性误差的时间特性，特征向量代表相应的空间特性.第一主成分对应最大的方差(特征值)，具有整网最多的信息，即携带了原始数据最大化的方差信息，它往往反映整个网的共同变化趋势；越靠后的主成分分量，所携带的整个网的信息越少.
5 国家CORS站噪声估计目前在估计国内CORS站运动速度参数时，没有考虑有色噪声对我国CORS站速度估计的影响.估计的CORS站速度值存在偏差，速度精度会被高估几倍、甚至一个数量级.
本文首先采用地表质量负荷改正减小了国家CORS站原始坐标时间序列中的地壳非构造形变信息，其次采用PCA 空间滤波方法去除了CGCS2000框架下公共误差信息，特别是去除了国家CORS网共有的非线性变化信息.在此基础上最后采用了白噪声、闪烁噪声和随机游走噪声组合的误差模型，采用最大似然准则将函数模型的未知参数和噪声分量一起进行了估计，并采用加权最小二乘估计方法，估计了CGCS2000下我国CORS 站运动速度参数和速度精度.建立了拟合函数模型如下[9]:
|
(4) |
式中ti是以年为单位的时间序列中解算历元，a为常数项，b为线性速度，Oj为在Tj时刻CORS站的位移，H表示阶梯函数.c和d组合表示全年性周期运动，e和f组合表示半年性周期运动.最后一项表示发生在历元Tj处大小为Oj的j0 个偏移常量，vi代表残差序列.
5.1 最大似然估计在有色噪声估计方法中，我们假定公式(1)中误差序列v由三个变量的线性组合构成，
|
(5) |
式中α(t)代表同一分布的单位权方差随机变量，β(t)、ρ(t)代表时间相关的随机变量序列，v(t)是时间序列的残差序列，比例因子a代表白噪声数值，ba≠0、ca≠0 代表闪烁噪声和随机游走噪声数值，顾及有色噪声的时间序列的协方差矩阵可以表示为
|
(6) |
式中Cx是时间序列的协方差矩阵，a2、ba2和ca2分别是时间序列的白噪声分量、闪烁噪声和随机游走噪声分量参数估值，I是单位权矩阵，Jf(t)和Jrw(t)分别是闪烁噪声和随机游走噪声协方差矩阵[10].由Cx协方差矩阵构造如下最大似然函数:
|
(7) |
式中，N为时间序列的长度，C是协方差矩阵，$\hat{v}$为公式(4)采用此处定义的协方差矩阵C用加权最小二乘法求得的残差.
5.2 加权最小二乘估计方法本文作者利用最大似然估计方法得到的国家CORS站各坐标分量白噪声和有色噪声估计值，重新定义坐标时间序列协方差阵(即公式(3)中的协方差矩阵C)，对PCA 滤波和三角函数拟合后的坐标时间序列拟合残差结果重新定权，采用加权最小二乘估计方法，估计了CGCS2000下我国CORS站运动速度参数和速度精度.
6 结果与分析本文通过计算得到了31个国家CORS站质量负荷改正和空间滤波前后的运动特征信息，由于篇幅所限，表 1-表 3只列出了其中6个站结果.
|
表 1 大气、土壤水、积雪和海洋非潮汐质量负载引起的测站垂向周年运动的振幅和相位 Table 1 The annual vertical movement caused by atmospheric pressure, soil moisture and snow, non-tidal ocean mass loading |
|
|
表 2 测站垂向位置周年及半周年变化的振幅和初相 Table 2 Amplitudes and initial phases of annual and semi-annual variations of station vertical positions |
|
|
表 3 CORS站坐标重复性统计(mm) Table 3 Repeatabilities of station positions (mm) |
表 1列出了4项负荷效应对CORS站垂直位置的影响，从表 1 可以看出土壤水负荷效应对CORS站垂直位置影响最大，平均影响为6.08 mm, 乌鲁木齐站达到9.82 mm；大气负荷效应平均影响为3.98mm, 哈尔滨站达到5.26mm；积雪和土壤水的影响要小很多，平均在亚毫米级水平.
6.2 国家CORS网中的公共误差图 2为经过PCA 分析方法得到的地表质量负荷改正后国家CORS站的公共误差时间序列，公共误差序列在北分量周期振幅0.8 mm, 东分量周期振幅2.8mm, 高度分量周期振幅10.1mm.本文应用频谱周期图方法得到的公共误差序列频谱周期图见图 3.从图 3可以看出国家CORS 站的公共误差时间序列在北分量上表现出了周年和几年的周期性变化信息，但是信号强度均较弱，北分量表现出的几年的周期性信息，需要进一步研究分析.东、高度两个分量上均具有明显周年性变化信号，特别是高度方向上表现出信号强度较强的周年信号，其信号强度比水平方向大10倍左右.
|
图 2 CORS站公共误差序列 Fig. 2 Time series of CME |
|
图 3 CORS站公共误差序列周期图 Fig. 3 Periodic map of power spectra of CME |
表 2列出了6个国家CORS站负荷改正和空间滤波前后周年和半周年信号的振幅和相位估值，从表 2可以看出:原始坐标时间序列高度分量的周年信号比较明显，周年项平均振幅为8.91 mm, 最大振幅为12.52 mm；负荷改正后高度分量周年信号强度减弱，周年项平均振幅为3.89 mm, 最大振幅为6.74mm.负荷改正后高度方向振幅减少32%~88%，平均振幅减少近41%.负荷改正和空间滤波后时间序列的周年项平均振幅为1.99 mm, 最大振幅为5.39 mm, 比原始时间序列振幅平均减少83.7%，比负荷改正后高度方向振幅减少60%.说明了国家CORS 站高度方向的周年性变化主要是由质量负荷效应和CGCS2000 框架下公共误差引起，特别是框架公共误差引起的非线性变化不容忽视，通过质量负荷改正和PCA 分析方法大大消弱国家CORS站高度方向的周期性变化.
从表 3可以看出:负荷改正后北、东和高方向平均坐标重复性值为2.89mm、4.10mm 和8.11mm, 相对于原始坐标时间序列平均坐标重复性值分别减小了26.9%、20.7%和21.4%.负荷改正和空间滤波后平均坐标重复性值为2.73 mm、3.73 mm、6.64mm, 相对于原始坐标时间序列平均坐标重复性值分别减小了30.9%、27.8% 和35.6%.说明负荷改正和空间滤波后国家CORS 站坐标精确性和可靠性得到了较大提高.
6.3 国家CORS站噪声估计表 4列出了负荷改正和空间滤波后国家CORS站坐标时间序列在北、东、高三个分量上的功率谱指数估值和白噪声和有色噪声参数估值.通过表 4 统计结果不仅证明了CORS 站噪声不满足纯白噪声性质，白噪声也不是噪声的主要成分，而且证明了国家CORS站的噪声性质更接近于白噪声、闪烁噪声和随机游走噪声的叠加.
|
|
表 4 基于FN+VW+RW模型的CORS站坐标时间序列谱指数和噪声参数估值(单位:mm) Table 4 Results of spectral index and estimated parameters of noise property based on FN+VW+RW model (unit:mm) |
表 5列出了采用FN+VW+RW 模型估计的CORS站水平速度及其精度，为便于比较分析，表 5也列出了仅考虑可变白噪声的估计结果.
|
|
表 5 基于FN+VW+RW和VW模型估计的CORS站水平速度及其精度 Table 5 Estimated velocities and their uncertainties (Iff) on horizontal component based on the VW+FN+RW and VW model |
考虑可变白噪声加有色噪声模型的参数估计量方差不会随着观测量的增多逐步减小，可以更加客观地描述国家CORS站运动参数.从表 5所列结果可以看出，考虑可变白噪声、闪烁噪声和随机游走噪声的CORS 站速度值和精度与仅考虑可变白噪声模型的结果有较大不同:速度估值偏差小于5% 的站有76.8%，在5% ~10% 的站有12.5%，大于10%的站有8.9%；在速度值精度评估中误差放大1~2倍的站有37.5%，放大2~3倍的站有23.2%，放大3倍以上的站有14.3%.
7 结论(1) 本文介绍了提取CORS站坐标时间序列中蕴藏的非线性运动信息的科学方法和结果，为该领域的进一步深入研究奠定了基础.采用本方法可以分离地球物理因素引起的地表质量负荷变化和去除CGCS2000下我国CORS 网公共误差信息，减小国家CORS网整体非线性运动信息.负荷改正和空间滤波后平均坐标重复性值相对于原始坐标时间序列分别减小了30.9%、27.8%和35.6%.说明负荷改正和空间滤波后国家CORS 站坐标精确性和可靠性得到了较大提高.
(2) 在负荷改正和空间滤波后的国家CORS站坐标时间序列中包含有本地不同噪声信号影响，在水平和高度分量上均表现出不同的噪声性质，本文作者采用了功率谱分析方法证明了在CGCS2000框架下我国CORS 站坐标时间序列误差不具备纯白噪声的特性，白噪声、闪烁噪声和随机游走噪声的噪声性质是国家CORS 站坐标时间序列的基本特征.
(3) 目前在CGCS2000 坐标系下的国家CORS站速度估值是在白噪声的假定下估计的，只考虑白噪声性质时，会导致CORS 站速度估值存在偏差，也会使速度估值的方差偏小，没有反映实际精度.在利用国家CORS网技术维护和更新我国CGCS2000坐标框架和进行相关领域科学研究时，有必要采用较为严格的误差模型代替白噪声模型，对CORS站运动特征参数进行科学估计.
| [1] | 陈俊勇, 杨元喜, 王敏, 等. 2000国家大地控制网的构建和它的进步. 测绘学报 , 2007, 36(1): 1–8. Chen J Y, Yang Y X, Wang M, et al. Establishment of 2000 national geodetic control network of China and its technologically progress. Acta Geodaetica et Cartographica Sinica (in Chinese) , 2007, 36(1): 1-8. |
| [2] | Yang Y X. Chinese geodetic coordinate system 2000. Chinese Science Bulletin , 2009, 54: 2714-2721. |
| [3] | 陈俊勇. 中国现代大地基准——中国大地坐标系统2000(CGCS2000)及其框架. 测绘学报 , 2008, 37(3): 269–271. Chen J Y. Chinese modern geodetic datum: Chinese Geodetic Coordinate System 2000 (CGCS2000) and its frame. Acta Geodaetica et Cartographica Sinica (in Chinese) , 2008, 37(3): 269-271. |
| [4] | 秘金钟, 蒋志浩, 张鹏, 等. IGS跟踪站与国内跟踪站联合处理的框架点选择研究. 武汉大学学报(信息科学版) , 2007, 32(8): 704–706. Bei J Z, Jiang Z H, Zhang P, et al. On framework sites selection for unite-processing of IGS CORS and domestic CORS. Geomatics and Information Science of Wuhan University (in Chinese) , 2007, 32(8): 704-706. |
| [5] | 张鹏, 蒋志浩, 秘金钟, 等. 我国GPS基准站数据处理与时间序列特征分析. 武汉大学学报 (信息科学版) , 2007, 32(3): 251–254. Zhang P, Jiang Z H, Bei J Z, et al. Data processing and time series analysis for GPS fiducial stations in China. Geomatics and Information Science of Wuhan University (in Chinese) , 2007, 32(3): 251-254. |
| [6] | 蒋志浩, 张鹏, 秘金钟, 等. 基于CGCS2000的中国地壳水平运动速度场模型研究. 测绘学报 , 2009, 38(6): 471–476. Jiang Z H, Zhang P, Bei J Z, et al. The model of crustal horizontal movement based on CGCS2000 frame. Acta Geodaetica et Cartographica Sinica (in Chinese) , 2009, 38(6): 471-476. |
| [7] | 王敏, 沈正康, 董大南. 非构造形变对GPS 连续站位置时间序列的影响和修正. 地球物理学报 , 2005, 48(5): 1045–1052. Wang M, Shen Z K, Dong D N. Effects of non-tectonic crustal deformation on continuous GPS position time series and correction to them. Chinese J. Geophys. (in Chinese) , 2005, 48(5): 1045-1052. |
| [8] | Dong D, Fang P, Bock Y, et al. Spatiotemporal filtering using principal component analysis and Karhunen-Loeve expansion approaches for regional GPS network analysis. J. Geophys. Res. , 2006, 111(B3): B03405. DOI:10.1029/2005JB003806 |
| [9] | Williams S D P, Bock Y, Fang P, et al. Error analysis of continuous GPS position time series. J. Geophys. Res. , 2004, 109: B03412. DOI:10.1029/2003JB002741 |
| [10] | Williams S D P. The effect of coloured noise on the uncertainties of rates estimated from geodetic time series. J. Geodesy , 2003, 76(9-10): 483-494. DOI:10.1007/s00190-002-0283-4 |