参考站长期稳定性计算和变化规律研究是卫星大地测量中的重要研究内容，是监测板块运动、建立坐标参考框架的重要内容。

对于参考站长期稳定性的研究，国内外学者主要利用参考站长时间位置序列进行分析。文献[1]由两个相距50 m的GPS参考站1998~2001年3 a的时间序列，解算出其相对位置变化在0.3 mm以内。文献[2]则分析了3 a的某GPS站点日坐标时间序列，得出该点沉降了3 cm或者说是1 cm/a。文献[3]利用两个GPS参考站点的一年的时间序列，通过一定范围内参考站、CORS、IGS 跟踪站三者线性速度的一致性，确定参考站是否存在局部平面变形或沉降,量化参考站的长期稳定性。

上述研究主要侧重于位置时间序列中的线性变化规律，然而最新研究表明，参考站位置时间序列中不仅存在线性变化，而且存在显著的上下波动^{[4, 5]}。这种非线性变化的起因一般认为是由技术系统误差、非构造形变以及随机因素(噪声^[5])三部分组成。相比线性变化，对非构造形变等造成的非线性变化的研究并不多见。本文利用江西省连续运行参考网(Jiangxi continuously operating reference stations，JXCORS)参考站3 a的位置时间序列，对参考站长期位置变化规律进行了研究，着重对其非构造形变的规律和起因进行了剖析，力图寻找适合于省级系统长期位置变化分析方法。

JXCORS由覆盖江西省的62个永久性GNSS连续运行参考站组成，绝大部分站点为楼顶站(52个，楼层小于4层)，少数为土层站。建站参考了《IGS基准站建立规范》、《中国地壳形变监测网络基准站建立规范》等，接收机为天宝R3，天线为Zephyr Geodetic II。

利用GAMIT/GLOBK(10.4)求取JXCORS参考站位置时间序列。首先，利用GAMIT软件，引入BJFS、SHAO、WUHN、KUNM、TNML 5个IGS站点作为约束，进行基线计算；然后，用 GLOBK基线解算获取的松弛解h文件进行平差，从而得到各站点的天解位置时间序列；通过坐标转换，得到各个站点的N、E、U三个分量的时间序列(N、E、U为站心坐标系坐标分量)。图 1为DUCH的位置时间序列。

图 1 DUCH站的位置时间序列 Fig. 1 Position Time Series of DUCH

图选项

从图 1可以看出，参考站位置既有线性变化，也有非线性变化，且具有相似的趋势，说明有相似的影响因素。

参考站的线性变化可以分解为水平方向和垂直方向。水平方向上主要表现为大尺度下的构造运动，比较一定范围内参考站及其他跟踪站之间线性速度的一致性，可以确定参考站点是否存在局部平面变形^[3]。垂直方向上主要通过速率大小反映地面沉降或抬升。

赣江断裂带横贯江西全省，是地表上一个醒目的低谷带(赣江)，由3条断裂组成。水平和垂直速度场可能与该断裂带有一定关系，在文中的各图中各以黑线示出。

本文主要通过ITRF框架和欧亚框架下的水平速度场来分析参考站水平方向上的线性变化，在欧亚框架下与附近的陆态工程站点的水平速度场进行比较，分析参考站的局部平面变形。

按照GLOBK平差策略，当一个GPS测站的数据至少有2个历元且间隔不少于0.9 a时，可以估计速度。GAMIT/GLOBK解算结果描述的是测站的绝对位移和速度。为了获取JXCORS网的速度场，本文采用获取坐标单日解的中国及其周边的5个IGS站BJFS、SHAO、WUHN、KUNM、TNML进行基准定义，对JXCORS网基线解算的松弛结果h文件进行平差，获取ITRF2005框架下的速度场，得到的水平速度场见图 2。

图 2 参考站在ITRF2005框架下的水平速度场 Fig. 2 Horizontal Velocity Field Under ITRF2005

图选项

计算可知，62个站点E方向平均速度为33.667 mm/a，N方向的平均速度为-11.687 mm/a，水平方向的平均速度为35.643 mm/a，优势方向SEE109.13°。

采用GLOBK提供的速度场参考框架转换工具CVFRAME进行转换，表 1给出了进行框架转换的欧亚板块欧拉矢量参数ω_x、ω_y、ω_z。欧亚板块的欧拉矢量参数是由模型NNR-NUVEL-1A所确定的，该模型可以直接用作背景场，进而可以方便地计算区域变形场。将ITRF2005框架下JXCORS参考站的水平速度转换到以欧亚板块为运动背景场的框架下，得到的水平速度场如图 3所示。

图 3 参考站在欧亚框架下的水平速度场 Fig. 3 Horizontal Velocity Field Under Eurasia

图选项

计算得，62个站点在欧亚框架下的的E方向平均速度为8.657 mm/a，N方向的平均速度为0.573 mm/a，水平方向的平均速度为8.675 mm/a，优势方向为SEE89.56°。

为对成果进行比较，本文以陆态工程中江西省附近的18个陆态网站点的测站速度(欧亚框架)为参照基准^[6]，站点分布见图 4。在欧亚板块下，将其速度与本文计算的JXCORS参考站的速度进行比较。

图 4 JXCORS站点(三角)和文献[6]中的18个站点(圆形)分布图 Fig. 4 Distribution of JXCORS Stations(Triangle) and 18 Stations in Reference [6](Circle)

图选项

表 2给出了本文计算得到的速度与文献[6]中结果的差值，二者水平运动方向差异在1.63°，数值相差1 mm/a内，说明了局部平面变形很小，参考站在水平方向上稳定性良好。

从图 5可以看出，在垂直方向上，大部分测站的垂向速度是隆升的。计算可知，上升方向的平均速度为1.887 mm/a，下降方向的平均速度为0.906 mm/a。在垂直方向上，JXCORS 62个参考站的速度普遍很小。根据以上统计，JXCORS参考站在垂直方向上比较稳定。

图 5 参考站在ITRF2005框架下的垂向速度场 Fig. 5 Vertical Velocity Field Under ITRF2005

图选项

参考站的位置变化包含了线性变化和非线性变化，后者包含了技术系统误差和非构造形变以及噪声的影响，主要表现为周期变化^{[7, 8, 9, 10, 11, 12]}。已有研究表明，GPS位置时间序列存在的大部分周期信号被认为是真实的测站实际位移^{[13, 14]}，而这种位移部分由非构造形变造成^{[1, 10, 15]}。为了研究这种非线性变化的本质，本文首先实现形变信号与噪声的有效分离，建立最优噪声模型，进而根据非构造形变分析非线性变化的成因。

许多研究表明，GPS参考站的位置时间序列中不仅包含白噪声，还包含有色噪声^{[4, 11, 12, 16]}。利用最大似然估计计算不同噪声模型下的最大似然值(MLE)，通过最大似然值的比较来研究噪声特性。主要考虑白噪声(white noise,WN)、闪烁噪声(fliker noise,FN)、随机漫步噪声(random walk noise,RWN)、高斯马尔可夫噪声(Gauss Markov noise,GM)。蒙特卡罗模拟实验表明，当两种噪声模型的最大似然值之差大于3.0时，两种模型具有可区分性(95%的显著水平)。

图 6为按照上述方案计算的 “白噪声+闪烁噪声”与“白噪声”的最大似然值之差。

图 6 JXCORS参考站坐标分量的最大似然值之差 Fig. 6 Difference Between Maximum Likelihood Estimate (MLE) of JXCORS

图选项

在绝大部分测站上，“WH+FN”模型的MLE值大于“WH+RWN”和“WH+GM”；而“WH+FN+RWN”模型的MLE值大于“WH+FN”模型，且“WH+FN+RWN”中已经包涵WH噪声和FN噪声，可以完全代替“WH+FN”模型。所以，本文选取“WH+FN+RWN”为JXCORS参考站的最佳噪声模型。

引起地壳非构造形变的地球物理因素主要包括两大类^[5]：第一类是潮汐形变，第二类是地球表面流体圈中的大气和各态水的质量迁移引起的地表质量负荷变化，主要包括大气、非潮汐海洋、积雪和土壤水等质量负荷，也称为环境负载。第二类形变在进行JXCORS数据处理时没有进行改正，可能是导致参考站时间序列中非线性变化的部分来源，特别是垂直方向上时间序列呈现非线性变化的原因。本文主要对环境负载进行分析。

江西是内陆省，海洋非潮汐负载未考虑，所研究的环境负载主要包括大气负载和水文负载(包括土壤湿度和积雪深度负载等)。利用QOCA采用美国国家环境预测中心(NCEP)再分析数据计算了大气负载和水文负载造成的地表位移。其中，大气负载利用的是NCEP的大气地表压力数据，时间分辨率为6 h，空间分辨率为2.5°×2.5°；积雪深度和土壤湿度负载利用的是NECP的再分析数据，时间分辨率为1 d，空间分辨率为1.875°×1.875°。地表质量负荷主要作用于径向，主要影响CORS站的垂向位置变化，本节提到的位置序列特指垂向分量。

图 7中列出了QNAN站的位置垂向分量的相对变化，包括了原始位置、大气负载、土壤湿度和积雪深度负载引起的地表位移的垂向分量。

图 7 QNAN站的位置(垂向分量)时间序列 Fig. 7 Position (Vertical Component) Time Series of QNAN

图选项

为了判断大气负载和水文负载的改正效果，本文分3种方式对原始序列进行了改正，包括大气负载改正、水文负载改正、水文负载加大气负载改正，分别计算了改正前后时间序列的RMS差。RMS差定义为：

RMS的具体计算公式见文献[5]。RMS差为正，负载改正效果良好；若为负，则表示改正效果不理想。 计算62个参考站3种方式改正前后时间序列的RMS差，分别定义为RMS-RMS_atm、 RMS-RMS_hydro 、RMS-RMS_atm-RMS_hydro，结果见表 3，篇幅原因只列出了6个测站。

由统计结果可知，经过大气负载改正后，所有测站位置序列的RMS值增大；经过水文负载改正后，76%的测站位置序列RMS值减小。这表明在江西区域大气负载计算结果不太理想，水文负载计算结果是可用的。

已有的研究表明，参考站位置序列中的周年项和半周年项周期变化^{[4, 5]}在垂向分量上最为明显。本节通过周年项和半周年项的振幅变化，研究水文负载对参考站位置变化的影响。

在“WN+FN+RWN”噪声模型下，利用CATS软件，分别计算水文负载改正前后位置序列的年周期项和半年周期项的振幅和相位，部分结果列于表 4中。

数据统计得出，29%参考站经过水文负载改正后周年振幅减小，平均周年振幅从改正前的3.49 mm增大到改正后的3.85 mm；66%参考站经过水文负载改正后半周年振幅减小，平均半周年振幅从改正前的1.97 mm减小到1.88 mm。各参考站周年项变化见图 8。

图 8 负载改正前(细)、后(粗)JXCORS参考站垂向变化的周年项(正弦函数拟合,参考点2009.0) Fig. 8 Annual Terms of Vertical Position Variations for the Fiducial Stations Before (Thin) and After(Thick) the Corrections to Effects Caused by Mass Loading(Fitting for a Yearly Sinusoidal Curve, with the Reference Time as 2009.0)

图选项

水文负载改正后，虽然周年平均振幅增大10%，但初始相位有较大的变化(1~2 mon)，所以水文负载可能是JXCORS参考站垂向分量的周年变化的成因之一，但是还有其他因素的影响，例如温度变化造成的垂向位移等，需要进一步分析；半周年项的平均振幅减小5%，部分测站减小幅度较大(40%~50%)，所以水文负载是JXCORS参考站垂向分量的半周年周期变化的成因之一，对垂向分量的半周年变化的影响与测站位置有关。

图 8中，箭头长度表示振幅,东向逆时针旋转的方位角表示初相位方向，箭头指向东、南、西、北方向依次表示极值发生在0.25、0.50、0.75、0.0 a。

本文以JXCORS为例分析了CORS参考站的长期稳定性，分别从参考站位置的线性变化和非线性变化进行研究，得出了一些结论。

1) 欧亚框架下，JXCORS参考站的水平运动方向与其附近的陆态网站点具有一致性，运动速率没有明显差异，表明JXCORS参考站不存在明显的局部平面变形；垂直方向上，抬升和沉降的平均速率较小，参考站较为稳定。

2) 为了分离形变信号和噪声，利用最大似然法获取了JXCORS参考站位置时间序列的最优噪声模型——“白噪声+闪烁噪声+随机漫步噪声”。在最优噪声模型下，根据水文负载改正前后的测站位置垂向分量的周年、半周年项的振幅和初相的变化，水文负载是JXCORS参考站位置垂向分量半周年周期变化的成因之一，可能是周年变化的成因之一，需要从更多的影响因素进一步分析。