GPS原始数据处理时虽然改正了固体潮、海潮和极潮等引起的地面点形变，但没有对大气、非潮汐海洋和陆地水引起负载进行改正，导致GPS站坐标时间序列中包含了明显的季节性变化(Tregoning and van Dam，2005; Van Dam et al.，2007).GPS成果归算时，需要利用Helmert转换(一般为七参数)将“无基准约束”或“松弛约束”的解转换到指定的全球框架下(Heflin et al.，1992; Zumberge et al.，1997; 陈俊勇，2008)，例如国际地球参考框架 ITRF(或IGS框架)(Altamimi et al.，2011).Helmert 转换是一种刚性转换，两个网应该是几何相似性的.但目前为止，由IERS提供的ITRF是线性框架(即由参考历元的位置和线性速度来实现)，并未顾及地表负载引起的季节性变化.因此“无基准约束”或“松弛约束”的解与线性框架不是几何相似的，若仍然利用Helmert转换，尚未模型化的负载形变会影响框架转换参数的估值，而且这种影响的大小与地面网的分布有关(Collilieux et al.，2009; Tregoning et al.，2009;Collilieux et al.，2012).Tregoning等研究了大气负载对框架转换参数和站坐标的影响，但他们采用的是地球物理模型计算得到的地表形变(时间跨度较短)，而且只考虑了大气负载的影响，没 有考虑海洋和陆地水的影响(Tregoning et al.，2009). Collilieux等分析了SLR站的不均匀分布对平移参数的影响(Collilieux et al.，2012).与上述文献不同，本文首先从理论上解释尚未模型化的季节性地表负载及GPS不均匀分布的测站网对框架转换参数(包括平移、尺度和旋转参数)及站坐标产生影响的机制，然后利用GRACE重力场系数仿真的GPS地表形变来分析负载形变的影响量级，最后结合实际IGS测站分布给出框架转换时参数选取和测站选取的若干建议.

GPS原始观测数据的日常处理没有对非潮汐大气和海洋(即由大气和海洋的质量重新分布引起的地表形变)、陆地水引起的地表负载进行改正的主要原因是：(1)目前海底压强模型采用反变大气压模型(Inverted Barometer)或者非反变大气压模型，对长期变化(尤其是季节性变化)是有效的，但没有顾及海洋对大气的补偿在大于两星期的频段上的变化(Van Dam and Wahr，1987)；(2)地表质量负载改正模型不够精确，不同模型参数和输入数据得到的大气负载差异的RMS可达毫米级(Van Dam et al.，2002)；(3)不能保证地表质量负载不受固体潮、海潮等改正模型的不完善的影响；(4)非潮汐大气和海洋、陆地水的地表质量负载模型的总质量变化不闭合(Clarke et al.，2005).

受轨道动力学约束，若采用“无基准约束”，同时计算GPS卫星轨道和测站坐标，理论上得到是基于CM(Center of Mass of the Earth System)原点的站坐标.ITRF框架(或IGS框架)的原点是地球形状中心CF(Center of Surface Figure).不考虑外力作用时，CM原点是空间固定的.地表负载使CF相对CM发生相对运动，即为“地心运动”(Blewitt and Clarke，2003; Dong et al.，2003).利用Helmert参数转换可以将基于CM的瞬时站坐标转换到CF框架下，网平移法将转换参数的平移量视为地心运动.本文定义CF相对于CM的运动为地心运动，平移量的相反数即为地心运动.利用Helmert七参数转换将基于CM的坐标转换到已知CF框架的公式为

若观测方程的权阵为 P，则

利用公式推导从理论上解释“负载效应”和“网效应”的产生机制.假定两个测站网分别基于CM原点和CF原点，不存在旋转和尺度上的差异.t₀时刻无地表质量负载时，CM和CF框架下的测站i坐标分别为 X_CMⁱ(t₀)和X_CFⁱ(t₀).t时刻CM和CF框架下的负载形变分别为 ε_CMⁱ(t)和ε_CFⁱ(t)，那么t时刻的站坐标 X_CMⁱ(t)和X_CFⁱ(t)分别为

由于t₀时刻CM和CF框架是重合的，X_CMⁱ(t₀)= X_CFⁱ(t₀).t时刻受地表质量重新分布的影响，CF原点相对于CM原点发生相对运动，即地心运动为 G(t)，那么任意时刻

下面将详细推导负载效应和网效应对七参数的影响.为推导方便，做一些符号简化.

令

取单位权时，式(3)可以简化为

首先分析负载效应和网效应对平移参数的影响.由于GPS卫星轨道受复杂的非保守力影响，利用网平移法得到的平移参数与“真实地心”存在差异，更准确地说利用网平移法得到的地心运动为“似地心运动”.

理想情况下，利用式(6)可得到准确的平移参数，即 T =- G .但实际数据处理时，基于CM原点的瞬时坐标 X_CMⁱ包含了地表质量负载，而CF框架下的长期坐标 X_CFⁱ并不包含地表质量负载，此时

利用式(9)求得的平移量受CF框架下的地表质量负载的影响，即负载效应的影响.负载效应会导 致平移参数向负载形变大的方向移动.若测站数量 少(对地表形变采样不足)或分布不均匀，会进一步加剧这种影响，即为网效应.

对于尺度参数，若CF框架下的站坐标已包含地表质量负载的影响时

若测站沿各轴分布是均匀的，即

最终可得

这样利用Helmert转换求得的尺度参数为0，与假设一致.注意此时测站必须近似均匀分布，才能保证 .也就是说，即使没有负载效应的影响，不均匀分布的地面网也会导致尺度参数估值不为0.因此从理论上讲，不均匀测站分布对尺度的影响比平移参数还要显著.进一步，若顾及CF框架下的坐标并不包含地表负载形变，那么

即使测站近似均匀分布，受负载效应的影响，仍然有 U ^T l ≠0.同理对于旋转参数，不考虑CF框架下的地表质量负载的影响时，

若测站沿各轴均匀分布

这样利用Helmert转换求得旋转参数为0，与假设一致.而且根据式(15)，由于负号的存在，不同轴的网效应有可能相互抵消，因此不均匀的地面网对旋转参数的影响比对尺度参数的影响要小得多.与平移参数类似，旋转参数实际上还受CF框架下的负载效应的影响，同样由于负号的存在，这种影响仍然比尺度参数要小.需要注意的是，由式(7)可知，由于(A ^T A)^-1的影响，七参数不仅受负载效应及网效应的影响，七参数之间也会相互影响.

本文利用GRACE重力场球谐系数仿真计算由地表负载引起的GPS测站的形变，并基于仿真数据研究负载效应和网效应对框架转换及测站坐标的影响.采用CSR提供的RL05版本的重力场球谐系数作为仿真输入数据，截止至60阶.时间跨度为2002年8月—2014年4月.由于GSM系数不包含大气和海洋的高频非潮汐变化，首先利用GAC系数进行恢复.考虑到GRACE数据处理采用的是CM框架，一阶球谐系数为0.需要利用Swenson提供的一阶项进行恢复(Swenson et al.，2008)，得到CF框架下的1～60阶的重力场球谐系数.GRACE RL04版本的ΔC₂₀包含明显的由K₂(3.73a)和S₂(0.44a)潮引起的混叠信号(Chen et al.，2009).虽然RL05版本的ΔC₂₀质量有所提高，但仍比SLR差.因此仍然需要用CSR提供的ΔC₂₀替换RL05的系数.另外由于GRACE系数的高阶噪声较大，采用500km的高斯滤波进行去噪.最后将重力球谐系数转换为表面密度球谐系数，两者间的转换关系参见相关文献(Wahr et al.，1998).

根据Farrell的负荷理论(Farrell et al.，1972)，可由负荷Love数和表面密度球谐系数计算测站的 地表负载形变.一阶负荷Love数与框架原点有关，二阶及以上负荷Love数与框架原点无关.CF框架下的一阶负荷Love数为h₁^CF=-0.268，l₁^CF=-0.134；CM 框架下的一阶负荷Love数为h₁^CM=-1.290，l₁^CM=-0.887(Blewitt，2003).二阶及以上负荷Love数采用Farrell提供的数值(Farrell，1972).由此可分别计算得到CM和CF框架的下地表形变.同时还可得到CF原点和CM原点间的相对平移量(即地心运动)，作为真值供比较分析.为了分析网效应对框架转换的影响，本文采用了IGS跟踪站和全球均匀分布的格网两组测站网，并计算了两组测站网基于CM和CF原点的形变时间序列.由于IGS框架是GNSS精密用户使用的参考框架，IGS精密产品都统一于IGS框架，精密用户的成果也都归算于IGS框架或者至少提供与IGS框架的转换参数(Rebischung et al.，2012).因此第一组测站网选取IGb08框架的所有232站(用“sta232”表示)和91个核心站(用“sta91”表示)，如图 1所示.232站是91站的全球加密，特别是欧洲和北美地区测站非常密集；第二组由10°间隔的全球均匀分布的格网组成(用“grid”表示)，格网点的高程统一设置为0.

图 1

Fig. 1

图 1 IGS框架点全球分布图(91个核心站用红色标记) Fig. 1 Site distribution of IGS reference frame station (91 core stations are shown in red)

选用IGb08为Helmert转换时的已知CF框架.我们知道IGb08是线性参考框架，第一组测站网在已知框架下任意历元的坐标可由IGb08框架提供的测站在参考历元的位置和线性速度得到.认为10°格网点的线性速度为0，可仿真得到10°格网点在已知CF框架下的坐标.为后文叙述方便，将这两组测站对应的已知CF线性框架统称为IGb08.然后分别利用两组测站网在CF框架下的负载形变，可以得到包含负载形变的站坐标时间序列，作为瞬时CF框架.

本文共分三种策略计算转换参数，分别为：(1)以包含地表负载形变的、基于CF原点的站坐标时间序列为基准(即瞬时CF框架)，将基于CM原点的 站坐标时间序列转换到瞬时CF框架下，标记为“cm_cf”；(2)以IGb08为基准，将基于CM原点的站坐标时间序列转换到长期CF框架下，标记为“cm_igb08”；(3)以IGb08为基准，将基于CF原点的瞬时站坐标转换到长期CF框架下，标记为“cf_igb08.

由(3)式可计算平移参数，取其相反数即为地心运动.图 2给出了“cm_igb08”和“cm_cf”两种情况下计算得到的地心运动(用GX、GY和GZ表示)，并与真值进行比较.结合表 1给出的地心运动的计算值与真值的RMS可以看出，无论采用均匀分布的格网，还是实际IGS站，“cm_cf”恢复的地心运动都比“cm_igb08”要好.因为采用“cm_cf”策略时，采用的是基于CF原点的瞬时框架(包含了负载形变)作为基准.因此理论上讲“cm_cf”恢复的地心运动不受负载效应的影响.采用10°格网，GX、GY和GZ计算值与真值的RMS最小，分别为0.49 mm，0.35 mm和0.46 mm.IGb08框架并不包含地表负载形变，由(9)式可知，“cm_igb08”计算的 地心运动会受地表负载的影响.因此“cm_igb08_grid”和“cm_cf_grid”间的差异在一定程度上代表了负载效应对地心运动的影响：对Y方向和Z方向的影响较大，RMS大约为0.1 mm.

表 1(Table 1)

表 1 地心运动计算值与真值的RMS Table 1 Averaged RMS of the differences between recovered geocenter motion and true value without_scale/mm with_scale/mm GX GY GZ GX GY GZ cm_cf cm_cf grid 0.49 0.35 0.46 0.49 0.35 0.46 sta232 0.49 0.35 0.46 0.49 0.35 0.46 sta91 0.49 0.35 0.46 0.49 0.35 0.46 cm_igb08 cm_igb08 grid 0.51 0.43 0.55 0.51 0.43 0.56 grid_75 0.50 0.43 0.49 0.50 0.43 0.49 sta232 0.60 0.38 0.68 0.57 0.38 0.58 sta91 0.53 0.45 0.48 0.53 0.45 0.49 注： “without_scale”表示无尺度参数的六参数Helmert转换；“with_scale”表示七参数Helmert转换(后同).

表 1 地心运动计算值与真值的RMS Table 1 Averaged RMS of the differences between recovered geocenter motion and true value

图 2

Fig. 2

图 2 仿真计算得到的地心运动时间序列(不估计尺度参数)，时间间隔2012年8月—2014年4月 Fig. 2 Geocenter motion time series derived from simulated surface deformation (no scale parameter) from August 2012 to April 2014

由式(9)还可知，网效应对地心运动产生影响的根本原因是因为没有顾及CF框架下的负载形变，而负载形变的大小与测站位置有关.“cm_cf”转换时不受负载效应的影响，因而测站分布对地心运动也几乎没有影响：表 1中10°格网、232站和91站恢复的地心运动与真值的RMS几乎相同.不估计尺度参数时，“cm_igb08”情况下，测站分布对地心运动的影响比较明显，均匀分布格网的RMS总体比“sta232”和“sta91”要好.虽然“sta91”比“sta232”测 站数量少，但“sta91”X和Z方向的RMS比“sta232”要好. 这是因为“sta91”在X和Z方向上的分布更加均匀.232站分布在X、Y和Z轴正半轴的比例分别为59%，58%和69%；而91核心站在X、Y和Z轴正半轴的比例分别为51%，58%和52%.“sta232”和“sta91”Y方向的分布都比较均匀，而“sta232”的GY的RMS比“sta91”小，说明在保证测站均匀分布的前提下，尽可能多的测站有利于地心运动的恢复.即使采用“sta232”仍可恢复90%以上的地心运动的方差(Variance Explained)，采用“sta91”则至少可解释95%以上的地心运动的方差.

另外，考虑到GRACE球谐系数在高纬地区的质量不好，本文还利用纬度-75°~75°的10°格网组成的测站网计算了地心运动，结果见表 1中“grid_75”所示.“grid_75”与真值的RMS比“grid”要小，特别是Z方向.这是因为“grid_75”剔除了南北两极格网点对Helmert转换的影响，而极区仿真的地表负载形变的质量不好.

理论上尺度参数只与GM和C有关，而实际上还受一些技术相关性的误差的影响(Petit and Luzum，2010).例如GPS尺度就呈现出明显的季节性变化，图 3灰线所示的是IGS分析中心的周解与IGS框架的尺度参数(此处给出的是cod/mit/gfz三个分析中心的平均值，用“ac_igs”表示).

图 3

Fig. 3

图 3 尺度参数时间序列 Fig. 3 Scale time series

根据式(10)的推断，由于平移参数不为0，即使没有负载效应的影响，测站要近似均匀分布，尺度参数才为0.而实际上采用“cm_cf”策略时，无论采用哪种测站网，尺度参数均小于0.064 mm(约0.01ppb). 因此在这种情况下，即使测站网不严格满足均匀分布，平移参数对尺度参数的实际影响也是很小的，统计意义上可以认为尺度参数为0.

另外采用相同的测站网时，“cf_igb08”和“cm_igb08”两种策略下计算得到的尺度参数几乎相同，如图 3所示.“cf_igb08”转换时，对于式(10)，l_i=-ε_CFⁱ；对于“cm_igb08”，l_i=-(G+ε_CFⁱ).根据前面给出的结论，平移参数不会通过测站分布对尺度参数产生影响.由此不难推断相同测站网的“cf_igb08”和“cm_igb08”的尺度参数是非常接近的，不同测站网求解的尺度参数差异很大，因为影响尺度参数的地表负载与测站分布密切相关.

如图 3所示尺度参数有明显的季节项变化.由图 4的频谱分析可以看出，“ac_igs”除了呈现出较为明显的周年变化外(约0.768 mm(0.12ppb))，半周年变化也较为明显，除此之外还存在大约120天左右的周期项；而“cm_igb08_sta232”只有明显的周年变化，没有半周年和120天的周期项.“cm_igb08_sta91”也存在明显的周年项，但振幅比“cm_igb08_sta232”要小很多，同样没有半周年和120天的周期项.已有研究表明120天的周期项可能是由太阳光压模型引起的(Arnold et al.，2014)，GRACE系数模拟的地表形变显然不包含这一误差源.虽然 “sta232”的尺度参数与“ac_igs”符合得似乎更好，但考虑到尺度参数受网效应影响较大，“ac_igs”是利用IGS框架的91个核心站得到的，而“sta232”Z方向的分布极不均匀，而且“sta91”与“grid”的频谱也更为接近，因此笔者认为“cm_igb08_sta91”能够更真实的反映未建模的地表负载对尺度参数的影响.根据图 4给出的“ac_igs”和“cm_igb08_sta91”的振幅谱，可以推断地表负载及GPS实际测站的不均匀分布可以解释30%的GPS尺度的周年变化(约0.256 mm(0.04ppb))，但不能解释GPS尺度的半周年变化.另外，“cm_igb08_grid”和“cm_igb08_sta91”频谱间的差异反映的是IGS 91个核心站的不均匀分布对尺度的影响，大约可以解释不到10%的GPS尺度的周年变化.

图 4

Fig. 4

图 4 尺度时间序列的频谱图 Fig. 4 Amplitude spectra of scale time series

如表 1所示，测站均匀或近似分布时，是否估计尺度参数对地心运动几乎没有影响，例如“grid”和“sta91”.估计尺度参数主要影响“cm_igb08_sta232”的GZ.为了解释估计尺度参数后“cm_igb08_sta232”的GX和GZ的RMS都减小，图 5给出了“cm_igb08_sta232”七参数相对于六参数得到的GX和GZ的差值(ΔGZ和ΔGX).从图中可以看出，ΔGZ和ΔGX与尺度的季节性变化是一致的.由式(9)和式(13)可知，负载效应对尺度参数和地心运动(即平移参数的相反数)的影响是相反的，“cm_igb08_sta232”的尺度参数通过式(7)中的(A ^T A)^-1影响地心运动的估值，抵消了部分负载效应对地心运动的影响，导致估计尺度参数时“cm_igb08_sta232”GZ的RMS大致提高了0.1 mm.

图 5

Fig. 5

图 5 采用232站时，七参数和六参数转换得到的地心运动的差值和尺度参数(“cm_igb08_sta232”) “ΔGX”和“ΔGZ”分别表示X方向和Z方向的差值，“scale”是七参数转换得到的尺度参数，并转换为mm(1ppb = 6.4 mm). Fig. 5 Difference of geocenter motion derived from seven-parameter transformation and six-parameter transformation as well as estimated scale parameter, when 232 stations are adopted ΔGX and ΔGZ denote the differences of geocenter motion in X and Z direction, respectively. “scale” denotes the estimated scale parameter, which is transformed to mm by factor of 6.4.

计算结果表明，无论采用哪种测站网以及是否估计尺度参数，“cm_igb08”计算得到的旋转参数都不显著，与预期相符，这里就不展开讨论了.

如表 2所示，采用均匀分布的格网时，Helmert转换之后的站坐标与CF框架下真值的残差在ENU方向上的平均值分别为0.16 mm，0.25 mm和0.25 mm(“grid”)，在忽略计算误差和模型误差的前提下，我们认为这一量级代表了负载效应的影响.如图 6—8所示和表 2所示，采用实际IGS站后，RMS都有所增加，反映的是网效应的影响.测站分布不同时，网效应的差异也十分显著.如图 8所示，232站在欧洲和北美区域的测站分布十分密集，特 别是欧洲区域.不管是否估计尺度参数，“sta232” 的U方向在欧洲区域的RMS均值可达0.8~1.2 mm，大约是“sta91”的两倍.

表 2(Table 2)

表 2 Helmert转换得到的站坐标与真值的RMS均值 Table 2 Averaged RMS of the differences between station coordinates derived from Helmert transformation and true value without scale/mm with scale/mm E N U E N U cm_igb08 grid 0.16 0.25 0.25 0.16 0.25 0.51 grid_80 0.16 0.26 0.24 0.16 0.26 0.50 sta232 0.29 0.44 0.42 0.25 0.34 0.87 sta91 0.31 0.34 0.26 0.32 0.37 0.61

表 2 Helmert转换得到的站坐标与真值的RMS均值 Table 2 Averaged RMS of the differences between station coordinates derived from Helmert transformation and true value

图 6

Fig. 6

图 6 采用10°格网(grid)时，Helmert转换的测站坐标与基于CF原点的站坐标真值的RMS分布 (A)：六参数Helmert转换(无尺度参数)；(B)：七参数Helmert转换. Fig. 6 RMS of differences between station coordinates derived from Helmert transformation and true value in the CF Frame, when 10°grid are used (A)Six-parameter transformation (no scale parameter);(B)Seven-parameter transformation. From the top to the bottom are the East, North and Up components, respectively.

图 7

Fig. 7

图 7 同图6，测站网为IGb08的91个均匀分布的核心站 Fig. 7 Same as Fig.6 but 91 IGb08 core reference stations are used

图 8

Fig. 8

图 8 同图6，测站网为IGb08的所有232站 Fig. 8 Same as Fig.6 but 232 IGb08 reference stations are used

对于相同的测站网，估计尺度参数时，U方向的RMS均增加.采用均匀分布的格网时，U方向的RMS增加大约0.2 mm(约0.03ppb)，与4.2节中尺度参数的量级是对应的.这一结果与预期相符，因为尺度主要影响U方向.采用实际IGS站后，U方向的RMS增大，特别是“sta232”的欧洲区域，误差均值可达1 mm.估计尺度参数后，“sta232”的E方向和N方向RMS有所减小，这是因为如4.2节所示，估计尺度参数提高了“sta232”的地心运动的精度，进而影响了站坐标的EN方向；但估计和不估计尺度参数的测站的RMS在EN方向上的空间分布仍然基本相同.为了得到准确的站坐标，应不估计尺度参数并尽量采用均匀分布的测站.例如表 2中“sta91”站，不估计尺度参数恢复的站坐标残差的RMS在ENU方向上大约为0.3 mm，比“sta232”更好，特别是U方向上.总体来看，因此相对于232站，目前IGS采用91个全球均匀分布的核心站作为框架转换的基准是合理的，特别是对于研究欧洲区域而言.

本文的仿真实验表明，基于网平移法，采用实际的IGS站至少能够恢复90%的地心运动信号.也就是说，在不考虑GPS系统误差的前提下，实际的IGS地面网已经具备了监测地心运动的能力(即使本文采用的是这种对测站分布要求较高的网平移法).对于研究地心运动而言，在选取均匀分布的测站的前提下，还应该保证尽量多的测站，以实现对地表负载形变的足够采样.框架转换时，尺度参数会吸收部分尚未模型化的地表负载，对尺度参数影响的大小还与测站分布密切相关.地表负载及GPS实际测站的不均匀分布可以解释大约30%的GPS尺度的周年变化(约0.04ppb)，但不能解释GPS尺度的半周年变化.在测站均匀分布的前提下，是否估计尺度参数对地心运动(平移参数)的恢复影响不是很大.估计尺度参数对测站的U方向影响很大，特别是测站不均匀分布时，导致测站密集区域的U方向失真.因此相对于IGb08所有框架站，目前采用91个均匀分布的核心站作为框架转换的基准是合理 的，特别是对于研究欧洲区域而言.因此框架转换 时，应尽量选取均匀分布的测站，同时不估计尺度参数.