1 引言

进入21世纪以来，欧美发达国家先后实施了CHAMP(CHAllenging Minisatellite Payload)、GRACE(Gravity Recovery and Climate Experiment)和GOCE(Gravity field and steady-state Ocean Circulation Explorer)卫星重力测量任务.在三个新一代卫星重力计划中，卫星跟踪卫星技术得到广泛使用，卫星跟踪卫星观测数据成为解算重力场模型的主要数据源之一，包含高-低和低-低模式的卫星跟踪卫星数据.应用的基本步骤则是首先利用卫星重力观测数据确定全球重力场模型，然后以该模型为基础，开展相关的科学研究和工程应用.因此，如何利用卫星跟踪卫星观测数据解算高精度高分辨率的卫星重力场模型成为卫星重力测量任务中的关键环节.

在卫星重力场建模领域，有能量法(Wolff，1969; Jekeli，1999)、动力法(Reigber，1989)(指一般意义上的两步法)、加速度法(Reubelt et al.，2003)和短弧法(Mayer et al.，2005)等.这些方法的共同点是需要已知重力卫星的精密轨道，然后由卫星精密轨道，通过轨道摄动原理解算重力场模型.与此相对，同解法是“一步法”，直接利用原始跟踪数据同时确定重力卫星的精密轨道以及卫星重力场模型.同解法使观测数据和观测模型之间达到了更高程度的“同化”，理论上更加严密，可实现更优的参数估计.国际上的著名地学研究机构GFZ(GeoForschungsZentrum)，就是通过直接处理原始的跟踪数据，同时解算轨道和重力场模型，该方法较一般的“两步法”，观测模型更加严密，解算的精度也更高(Zhu et al.，2004).

同解法实现复杂，其具体技术路线并未见诸公开刊物，在这个特定领域，国内一直处于跟踪状态.目前，在卫星精密定轨技术方面，国内已经开展了卓有成效的工作，在运动学定轨、简化的动力学定轨、自适应定轨等方面，取得了较为全面的成果.国内展开的定轨方法研究，覆盖面广、技术路线全面，包括了非差、单差等多种技术模式，径向精度总体上已经接近3~5cm的水平.这些工作使得我国在低轨卫星精密定轨的数据处理领域，达到了国际先进水平(郑作亚，2004；彭冬菊和吴斌，2007；赵齐乐等，2008；李建成等，2009；秦显平，2009).在卫星重力场建模方面，国内多个单位采用“两步法”研究了基于卫星跟踪卫星技术的重力场建模(周旭华，2005；周旭华等，2006；肖云，2006；邹贤才，2007；张兴福，2007).在后续的模型研制上也取得了长足的进步，并且在部分方法上做了一定的改进(陈秋杰等，2013)，建立了静态和时变重力场模型(沈云中等，2015；王长青等，2015).但是，在同时完成卫星精密定轨与重力场解算的同解法上，目前国内还没有完全突破.而国际上，GRACE数据产品的官方机构CSR与GFZ，采用的都是同解法提供的产品.可以说，同解法是卫星重力场模型研制中的一项关键技术，国内目前尚未完全掌握，当前的工作主要是在仿真研究的基础上(邹贤才等，2010a)，对关键问题分别进行技术攻关.作者通过多年工作，从仿真研究开始，到动力法中的关键技术，比如并行计算、星载数据的精细处理等等，解决了一系列技术问题，不断完善方法的同时，从无到有建立了软件平台，完成了同解法平台的联调联试.

本文主要从观测模型的特征分析、数据处理的主要方案等方面展开，通过实测数据处理剖析同解法的特点.目的是希望能推动同解法在海量实测卫星观测数据处理中的应用，提高我国在卫星重力场模型研制领域的技术水平，进一步缩小与国际一流水平的差距，并促进我国对国际上新一代重力卫星数据资源的有效利用.上述研究对我国发展自主卫星重力测量系统，以及卫星重力模型的研制理论和实际应用都有积极的促进意义.

2 理论模型与方法

在分别论述各关键技术之前，先对本文的术语做简单约定，以免混淆.后文所述的动力法，相对同解法(SSM，Simultaneous Solution Method)而言，可称“两步法”.它以卫星精密轨道为观测值，利用动力学方法确定重力场模型和其他有关参数(周旭华，2005；周旭华等，2006；肖云，2006;邹贤才等，2010a).几何法，则是利用星载GPS(Global Positioning System)观测数据，按照动态定位的原理所做的精密定轨，也称为运动学定轨.另外，本文依然采用动力法中的分弧段思想，即将观测数据分解为多个弧段，每个弧段将轨道初值以及相关模型参数设置为未知参数.通过多弧段法方程综合获得最终的参数解算结果，其基本数学模型参见邹贤才等(2010b).

利用同解法处理星星跟踪数据的实质是如何利用星载GPS观测数据和低低星间观测值，同时确定卫星的轨道初值、地球重力场模型位系数、加速度计校准参数、模糊度、卫星接收机钟差、以及其他相关参数.在分别组成了相位观测数据以及星间距离变率的法方程后，可通过方差分量估计合理定权(崔希璋等，2009)，获得最终的参数解算结果.下面分别阐述星载 GPS 载波相位以及星间距离变率的观测方程.

2.1 星载GPS载波相位观测

基本的载波相位观测方程可简写为(Xu，2003)

(1)

式中，t_r和t_s分别为信号的接收和发射时刻，δt_r和δt_s分别为接收机钟和星钟的钟差.φ(t_r)为载波相位观测值，λ为波长，N为对应的整周模糊度.ρ(t_s，t_r)为发射天线相位中心到接收机天线相位中心间的几何距离；c为真空中光速.δ_iono为电离层延迟改正；δ_other为其他各项改正，包括地球自转、相位解缠(Wu et al.，1993)等，相关改正公式可参考文献.ε_LL(t_r)为载波相位观测的测量噪声.注意，本文处理的是星载GPS数据，因此无需考虑对流层改正.自然地，在实际的数据处理分析中，依然采用消电离层组合观测值.

根据星载GPS载波相位的观测方程及其非差组合观测值，以及选定的同解法参数，可以获得星载GPS载波相位观测值对应的设计矩阵.同时注意到，求偏导的基本原则是观测泛函对所有选定的参数P进行，因此设计矩阵有下列一般形式：

(2)

式中，F为观测泛函，与具体的观测值有关.但选定某个类型的观测值后，它的表达式一定是明确的，因此本文以F统一表达.ΔT、X₀、N、Acc、EGM分别为重力卫星的钟差、轨道初值(包括位置和速度)、模糊度、加速度计校准参数以及重力场模型位系数参数组成的列向量.为了建立任意历元观测值与轨道初值的关系，借助了观测历元的轨道位置和速度X_LEO.下标1，2，…，n表示历元编号.

对于上述载波观测数据的设计矩阵，它是由两部分组成的，一部分是ΔT和N.这两类参数加上卫星的位置参数，构成了几何法精密定轨的参数(李建成等，2009).另一部分是X₀、Acc和EGM.它们属于动力法的参数.因此，实现同解法的一个可能方案就是，通过几何法精密定轨与动力学方法两项技术，分别攻关后集成.这也正是本文研究采取的技术路线.

这其中有几个关键特征需要强调.在同解法的待估参数中，接收机钟差(后文简称为钟差)为每个历元设置一个.因此，在一个弧段中，钟差参数的个数等于有效观测历元数.对某个历元而言，对应某个载波观测值有且只有一个钟差参数，因此式(2) 中标示为“1”的偏导数，除对应历元的值为常数外，其余值都为零.这种特殊结构导致相位观测值对应的法方程系数阵具有一种稀疏结构，参见图 1.其中，①的对角部分表示与钟差有关的对应法方程，它是一个对角阵.按照本文的设置，弧段长度最大值设置为24 h，根据GRACE星载载波数据采样间隔10 s计算，该对角阵的维数可以达到8640.②是与模糊度参数有关的部分.对模糊度参数而言，它具有和钟差参数类似的特点.对某个特定历元的相位观测值，一定依附于某颗卫星，自然除了这颗卫星的模糊度参数外，观测值对其他模糊度参数的偏导数一定为零.因此也形成了对角阵.③是与轨道初值、重力场模型位系数、加速度计校准参数等力学参数有关的部分，一般情况下是满阵.对于②，一天的弧长包含的模糊度参数总数一般不超过500.因此，本文设计算法时利用①的稀疏性，而将②和③统一当作一般的矩阵进行处理，这样在消参数与回代求解(邹贤才，2007;邹贤才等，2010b)时，可以兼顾运行效率与实现的难易度.上述几类参数的非对角部分当作一般的矩阵处理.法方程的系数阵是对称阵，图 1展示了其下三角部分.

2.2 星间距离变率观测值的观测模型

记GRACE任务的两颗卫星分别为A和B，且它们在空间直角坐标系下的状态向量分别为X_A=和X_B=，则在空间直角坐标系下，星间距离ρ_ab与星间距离变率可以表示为：

(3)

和

(4)

式(3) 和式(4) 即所谓的观测模型.根据观测模型，可以获得星间距离变率对卫星状态的偏导数.由式(4) 可得星间距离变率对两颗卫星状态的偏导数为(周旭华，2005;邹贤才，2007)：

(5)

(6)

因此，可写出星间距离变率对应的设计矩阵B_krr为：

(7)

式中双星状态对力学参数的偏导数，与公式(2) 中的对应部分计算方法相同.因此，星间距离变率观测数据的法方程可以在动力法工作的基础上直接完成.由此，联合解算这两类观测数据的法方程，可以获得同解法中的各项参数.

2.3 模型设置以及参数配置

本节简要介绍本文采用的基本参数、基本的力模型配置和参数估计方案.参见表 1.

对于星载加速度计数据的处理，研究团队已经在相关文章中探讨过该问题，因而表 1仅列举了参数配置方案.有兴趣的读者可以参考邹贤才等(2015) ，后文将做进一步的数值检验.

最后简要讨论一下并行计算.已有研究表明，卫星重力数据处理中，分弧段处理以及法方程综合技术非常适合在共享内存的平台上实现，并且无须设计并行算法，采用批处理即可实现.在法方程综合阶段，计算热点之一是大型线性方程组求解.注意到，同解法的参数解算中，法方程均为对称阵，且正常条件下都是正定阵，自然可以直接采用已有成果(邹贤才等，2010b).本文的所有计算工作都在武汉大学测绘学院的计算集群上实施(含7个Intel节点，共192核)，搭载CentOS LINUX操作系统，任务管理通过SLURM实现.

3 数值分析与讨论

前已述及，同解法的最大特点是能以更统一的模型对参数进行整体求解，理论模型更加严密.因而同解法的解算结果中包含了多种信息，包括精化的重力场模型、加速度计校准参数、精化的轨道初值以及与GPS观测有关的接收机钟差、模糊度等.其中有些参数的解算结果目前没有条件进行检核，比如模糊度和钟差的解算结果.有些参数则可以设计合理方案进行检核.这正是本节的主要工作，通过科学合理的检核方案，对同解法相关产品进行评价.

3.1 轨道质量评估

与星载GPS几何法定轨不同(李建成等，2009)，同解法并不是直接给出卫星的精密轨道.它先是做力模型精化，然后通过轨道积分给出轨道.这一点与简化的动力学方法(Yunck et al.，1990)类似.同时考虑到JPL(Jet Prolusion Laboratory)的简化动力学轨道(GRACE官方发布的精密轨道)质量很高，为国内外许多研究群体采用，并且可以提取到时变信号(冉将军等，2014)，而且它与同解法一样，能提供动力学意义上的卫星速度，所以本文将同解法给出的轨道(一个月的数据)与JPL的轨道进行了比较.给出轨道径向(R: Radial)、切向(T: Tangent)和法向(N: Normal)差异的统计结果.参见表 2.从该表可以看出，两套轨道在 RTN 三方向的均方根不超过3 cm.其中，轨道切向的精度在2.7 cm左右，另外两个方向优于2 cm.上述结论对两颗卫星均成立.相比已有的几何法定轨结果(李建成等，2009)，同解法的轨道精度更高，不仅如此，在顾及两套轨道误差的情况下，轨道差异的极值也只有11 cm.进一步，两套轨道的平均偏差不超过1 cm，轨道的一致性好.这些结果从侧面表明本文在力模型配置上的合理性，说明动力学模型在定轨中起到了良好的约束作用.

如果选择一个单独的弧段(弧长23小时59分)比较，不仅可以很清晰地看出与 JPL 轨道的差异，还可以发现差值的分布特点.如图 2所示.从绝对值上讲，三个方向(RTN分别用红、绿、蓝表示)都在几厘米左右，可以认为两者的轨道精度相当，但是轨道差异中存在明显的周期性现象，且与GRACE卫星的轨道周期一致，这表明轨道残差依然受到了类似重力场因素的影响.合理的解释是简化动力学方法中采用的重力场没有精化所致——JPL的精密轨道是采用简化动力学方法建立的，不是纯粹意义上的动力学定轨，没有精化重力场模型这一过程.因此，在实现“最优”轨道的同时，它在力模型上的缺点通过轨道积分，特别是较长时段的轨道积分放大，比较明显地显示出来了.

为了进一步验证该思路，本文采用高精度的星间距离变率观测值对轨道进行检核.具体过程如下.简化动力学轨道与同解法的动力学轨道都能给出卫星的位置和速度，因此可以利用轨道导出星间距离变率观测值.通过导出量与实测值比较，顾及实测值达到10^-6m·s^-1左右的精度，可以对轨道和星间距离变率之间的一致性进行检验.结果参见图 3.该图中，振幅较大(红色)曲线是由JPL精密轨道导出的星间距离变率与观测值之间的差异(JPLD)，振幅较小(绿色)曲线是由同解法精密轨道导出的星间距离变率与观测值之间的差异(SSMD).很明显，同解法给出的精密轨道与星间距离变率观测值之间的一致性要好得多.数值上，JPLD分布在-7.3×10^-6/6.5×10^-6之间，RMS(root mean squares)为2.5×10^-6，SSMD分布在-1.3×10^-6/1.4×10^-6之间，RMS为2.7×10^-7.相较而言，SSMD不仅更加稳定，而且与星间距离变率更吻合，RMS要高一个量级.对于其他弧段，有类似的结果.因此，本文认为同解法给出的轨道、速度与精密星间距离变率观测值之间的一致性更好，这也从侧面验证了同解法在精密定轨的同时精化重力场模型的重要优势.

上述轨道质量检核是利用星间距离变率对同解法和简化动力学轨道位置和速度的综合检验.卫星激光测距(SLR，Satellite Laser Ranging)是一种轨道质量检核的保留手段，因此本文也实施了SLR检验.做地面跟踪站速度场、潮汐、偏心改正，卫星质心、相对论、大气折射以及时距偏差等改正(Hulley and Pavlis，2007)，得到星地距离检核卫星轨道，情况如表 3和表 4.以观测数据最为丰富的跟踪站7090(澳大利亚堪培拉站)为例，双星的径向精度在2 cm以内，进一步验证了表 2中结果的可靠性与轨道质量.

3.2 加速度计校准结果检核

加速度计校准是卫星重力数据处理中的一项基本工作，但目前没有准确可靠的校准参数作为外部检核标准.即使GRACE团队发布了技术文档(Bettadpur，2009)，也只能作参考.根据GRACE官方手册，卫星星固坐标系x轴实际上是两颗卫星之间的连线.因此这两颗卫星各自的星固系x轴重合，但是方向相反.沿轨的非保守力源主要是大气阻力.这两颗卫星的空间距离在200 km左右，因此大气阻力之间的相关性很强.在准备本文工作时，国外也正好有相关的成果发表(Kim and Tapley，2015).为此，下面展示这种相关性，从侧面验证同解法的非保守力模型精化效果.

利用同解法给出的加速度计校准参数对原始的加速度计数据进行校准，并表示在卫星的星固系中，这样就得到了x方向的非保守力.两颗卫星的空间相关性使得它们各自受到的非保守力也存在很强的相关性，同时注意到两颗卫星星固系x轴方向相反，因此可以预见它们的x方向的非保守力大小接近，符号相反.这一特性，为实测数据的处理结果所证实.参见图 4.其中，红、绿分别表示两颗卫星的非保守力，蓝色为其差别.

从该图可以看出，两颗卫星沿轨方向的非保守力振幅接近，但是相位相反.通过求和之后的分布(蓝色曲线)分析，可以更加清楚地看出这个特点.结合它的基础工作(邹贤才等，2015)，可以认为同解法的加速度计校准内符合也非常好.这里也同时给出GRACE 双星加速度计的比例系数供参考，其估计策略请参见表 1.GRACE-A 卫星的三轴比例系数及其精度分别为(0.9701，0.9255，0.9374) 、(7×10^-5，3×10^-3，4×10^-4)，GRACE-B卫星的比例系数为(0.9549，0.9418，0.9309) ，精度指标相同.

3.3 动力学模型配置检验

与3.2节不同，本节在更宏观的角度研究力模型配置的合理性，并主要通过重力场模型精度检核展开.目前国外发布的月重力场模型并不是都有精度信息，很难简单的直接应用误差阶方差以及GPS水准的检核方法，因此也做了同解法重力场模型(SSM_EGM)与两步法结果的对比.

基本考虑如下.同解法的观测模型涵盖了传统的动力学方法(“两步法”)的所有参数，作者设计的软件也支持运行模式上的切换.如果将它的输入数据修改为轨道，就可以用两步法反演重力场模型及相关参数.首先，在分析中测试利用同解法给出的纯动力法轨道，按照两步法反演重力场.结果是反演的重力场模型与同解法确定的重力场模型一致，其精度指标接近计算机的截断误差(在10^-14~10^-15之间).对于同一套软件，力模型配置可以完全保持不变，参见表 1(注意，相关动力学参数仍设置为待解参数进行估计，并不是简单的固定为同解法的结果).因此它只能看作是软件测试中的正反算.但是这个现象说明，如果基于JPL精密轨道的两步法反演结果不好，则说明双方在力模型配置上的差异比较明显.

检验思路如下.由于无法确切知道JPL在定轨中采用的力学模型配置，特别是简化动力学定轨(Yunck et al.，1990)中有可能视需要加入一些所谓的脉冲加速度信号，所以力模型上的差异使得本文无法完全恢复它在定轨中采用的重力场模型.但是，根据简化动力学方法定轨的基本原理，它需要选定某个重力场模型(显然它也是主要摄动力)，辅以其他保守力模型(包括各种潮汐以及非潮汐时变模型)和非保守力模型(可由加速度计提供，或者结合经验模型等等)，完成精密定轨.因此，通过简化动力学轨道反演重力场(记为 JPL_EGM)时，一方面，不会获得前述正反算中的内符合精度.另一方面，如果在力模型的配置上，本文的配置方案与JPL的方案不出现大的偏差，按照前述同解法轨道反演重力场的分析，也应该能很好地反演重力场模型，并能通过 GPS 水准进行独立的外部检核.本文依然采用美国的6169个GPS水准点.结果表明，SSM_EGM的精度为0.757 m，而JPL_EGM的精度为0.762 m(参与比较的截断阶次为96阶，参见表 1). CSR最新的GGM05S模型对应指标为0.755 m，顾及到 GGM05S 使用了多年的数据，上述精度表明本文解算的重力场模型质量也是可靠的，与国际机构其他发布模型的分阶次比较结果详见表 5，总体上看，随着重力场模型阶次的提高，模型大地水准面的精度不断提高；对不同的模型而言，在本文选择的截断阶上，精度相当.

最后，本文给出同解法重力场模型与国际上几种主要的重力场模型之间差异的对比，以及与 CSR 月重力场模型(同阶次模型UTCSR)之间的对比结果，参见图 5.基本的做法是假定EIGEN6C4的精度“最高”，然后各个模型分别与 EIGEN6C4作差，把差值看作位系数的误差，计算所有阶次上大地水准面的误差，从而表示重力场模型的精度.从作差后大地水准面每阶的差异来看，EGM2008与EIGEN6C4之间吻合较好，GGM05S和EIGEN6C4符合得最好，并且静态模型之间的符合程度一致优于月解模型SSMD和UTCSR.对于两个月解模型，SSMD在13阶以内，特别是J₂项上表现要好.在60阶以内，两者与静态模型的符合程度，各有优势，但是在20～35阶左右，UTCSR的差异更加稳定.60阶之后，UTCSR 模型和SSMD的表现比较接近.总体上讲，结合GPS水准检验与大地水准面误差分析，再综合前述的积分轨道的质量检验(参见图 3、表 3和表 4)，本文认为SSMD与国际上的同类型月解模型是相当的.

4 结论与展望

本文以方法研究为主，提出了同解法的一种实现方案，并结合GRACE一个月的实测数据，在高性能计算平台上实现了高质量卫星精密轨道与卫星重力场模型的同时解算.对动力法精密轨道、重力场模型以及相关附属产品做了较为全面的评价与分析.结果表明，卫星精密轨道与JPL轨道产品互检，RTN三方向的精度均优于3 cm；与SLR比较，径向精度优于2 cm.并且，本文通过星间距离变率检核了同解法的纯动力学轨道与JPL简化动力学轨道，结果表明同解法的动力学轨道与高精度距离变率之间符合更好.同时，本文还通过GRACE任务双星之间非保守力的空间相关性，以及GNSS水准对本文的力模型配置进行了检验，进一步证明了力模型设置以及相关产品的可靠性.

综上所述，本文的研究已经涉及了同解法流程中的各个环节.产品质量分析表明，本文已经突破了同解法的各项关键技术，在国内率先实现了从仿真研究到实测飞行任务数据分析的重要进展.

在同解法领域，未来仍然有许多问题可以做深入研究.比如，可以引入模糊度固定技术，探讨如何进一步提高轨道质量以及重力场模型的精度，甚至从方法上对星星跟踪观测模型严密化，在更高的层次上建立统一的观测模型，深化同解法及其产品在科学研究中的作用.上述工作涉及到理论、方法研究以及计算方案的大调整，值得做进一步的深入研究.

致谢

感谢JPL提供的GRACE任务观测数据以及武汉大学测绘学院计算中心的大力支持.