卫星精密定轨与卫星重力场模型确定的理论和方法
Ø 低轨卫星精密定轨算法的改进目前在低轨卫星运动学精密定轨中,广泛采用的参数估计方法是最小二乘算法,在利用最小二乘进行运动学定轨时,需要同时处理一天的观测数据,对计算机内存要求高,且效率低,对于1Hz采样率的星载GPS数据来说,基于最小二乘的运动学定轨方法难以在普通PC机上处理。为此,本项目提出了一种改进的Kalman滤波算法,设计并实现了基于Kalman滤波的运动学定轨方法,通过对GRACE卫星的实测数据计算表明,基于Kalman滤波的运动学定轨方法,不仅对计算机要求低,计算效率高,而且可以达到与最小二乘法相当的定轨精度,可以作为今后高采样率星载GPS数据的定轨方法。
为了验证改进算法的有效性,本项目采用的数据为GRACE卫星2009年10月1日至10月7日期间1个星期的星载GPS数据,其采样率为10s。在精密定轨中,GPS卫星轨道及钟差采用了IGS发布的精密星历及30s间隔的卫星钟差产品(http://igscb.jpl.nasa.gov/components/prods.html),卫星天线相位中心采用igs05.atx发布的模型,相位缠绕改正为Wu等的模型,接收机天线相位中心偏差改正采用了GFZ发布的数据。另外,对于低轨卫星的相对论效应和相位缠绕可以改正也可以不改正,如果不改正则会被吸收到接收机钟差中,也即本课题的处理方式。
利用改进的Kalman滤波算法对GRACE卫星1个星期的数据进行了定轨计算,并将定轨结果(称为“WHU”轨道),与GFZ发布的简化动力学轨道作为参考轨道(“REF”轨道)进行了比较,其轨道差值的精度统计如图1所示。从图中可以看出,三个坐标方向的精度均小于4cm,其中X轴方向定轨精度优于3cm,Z轴方向精度在2cm到4cm之间。
图 1 运动学轨道精度统计图
可以看出,卫星轨道的整体精度较高,且误差变化很平稳,其中轨道中存在的一些异常变化,并不是Kalman滤波引起的,其主要原因是相邻历元见卫星观测图形结构发生了变化,比如跟踪到了新的卫星,对于这种情况,即使采用最小二乘,也无法避免轨道中的“突刺”现象。
Ø 基于空域最小二乘法求解纯GOCE卫星重力场模型
深入研究了利用空域最小二乘法确定GOCE卫星重力场的理论、方法和技术,包括有色噪声的处理方法、正则化方法等;基于实测数据确定了一个220阶次的纯GOCE卫星重力模型,与EIGEN-6C2模型以及美国和中国GPS水准数据比较检验的结果说明GOWHU01S与国际同类模型精度一致。本课题主要使用的数据是GOCE提供的EGG_NOM_2,SST_PKI_2,SST_PRD_2,SST_PRM_2,SST_PCV_2,GG_CCD_1i产品,数据周期为近两年,采样率1s,根据数据的分类,数据处理方法如下:
1)采用点域加速度方法(PAA: Point-wise Acceleration Approach)组成SST的法方程,最高阶次恢复到130阶次;
2)引力梯度张量对角线的三个分量分别组成法方程,最高阶次为220;
3)联合SGG和SST对应的法方程进行重力场模型的解算,模型的最高阶次选择220阶,同时采用Tikhonov正则化技术,正则化仅作用于受极空白影响较大的引力场系数(m<20);四类观测值(SGG:Vxx,Vyy,Vzz和SST)的相对权比为3:3:6:1;
选择了目前国际上最新发布的仅用卫星重力数据求解的模型GOTIM04S、GODIR04S、GOCO03S来比较。为了方便比较,将本课题采用了正则化的SGG和SST联合解算得到的卫星重力模型命名为GOWHU01S,图 2是GOTIM04S、GODIR04S、GOCO03S、GOWHU01S求解模型与EIGEN-6C2模型差异的阶RMS,从图中可以看出GODIR04S和GOCO03S在中低阶次部分(小于160)与EIGEN-6C2较为接近,主要是因为这两个模型与EIGEN-6C2一样包含了GRACE观测数据的贡献,同时EIGEN-6C2采用了直接解方法的建立SGG观测数据的法方程;GOWHU01S在大约20阶次到170阶次范围内较GOTIM04S更接近于EIGEN-6C2,但高于170阶次GOTIM04S更接近于EIGEN-6C2,原因在于本课题并没有对高于180阶次的系数采用正则化约束;从这个图上基本能看出,本课题解算的GOWHU01S模型是可靠的,与仅用GOCE观测数据解算的模型GOTIM04S差异不大。
图 2 GOTIM04S、GODIR04S、GOCO03S、GOWHU01S与EIGEN-6C2模型差异的阶RMS
为了能够更可靠的评价GOWHU01S模型的精度,下面采用外符合检验的方法对解算模型进行精度评定,同时选择了GOTIM04S、GODIR04S、EGM2008模型,将其与中国和美国的GPS水准数据进行比较,表18给出了与中国和美国GPS水准比较的统计结果,需要说明这里模型比较的时候均是采用零潮汐模型。这里比较的时候没有考虑截断误差的影响,因为GPS水准数据我们可以认为包含全频段的大地水准面信息,从表中可以看出,当最高阶次取相等数值时,GOWHU01S与GOTIM04S和GODIR04S与GPS水准之差的标准差基本相等,当模型取200阶次时,标准差的差异在毫米量级,所有卫星重力模型的精度均优于EGM2008,在中国区域几个卫星重力模型对应标准差的差异小于1mm,当取到220阶次的时候卫星重力模型的精度远高于EGM2008模型的精度。
表格 1 GOWHU01S、GOTIM04S、GODIR04S、EGM2008四个模型与中国649个GPS水准点比较的统计结果(单位:厘米)
Model |
Degree |
Mean |
Max |
Min |
RMS |
STD |
GOWHU01S |
200 |
-0.50 |
295.2 |
-323.2 |
57.14 |
56.97 |
220 |
-0.58 |
258.4 |
-237.2 |
51.26 |
50.97 |
GOTIM04S |
200 |
-0.49 |
323.3 |
-303.3 |
57.10 |
56.93 |
220 |
-0.57 |
237.6 |
-233.5 |
51.23 |
50.96 |
250 |
0.68 |
145.7 |
-189.2 |
44.19 |
43.70 |
GODIR04S |
200 |
0.47 |
321.6 |
-303.4 |
57.15 |
57.00 |
220 |
0.55 |
232.2 |
-237.1 |
51.22 |
50.97 |
260 |
0.66 |
131.3 |
-187.8 |
43.69 |
43.22 |
EGM2008 |
200 |
0.47 |
383.1 |
-288.2 |
60.29 |
60.15 |
220 |
0.58 |
291.1 |
-258.0 |
54.81 |
54.55 |
2190 |
-23.94 |
153.4 |
-173.0 |
33.91 |
24.04 |
Ø
基于GOCE和GRACE卫星重力场模型位系数的谱组合法研究了联合GOCE和GRACE重力场模型位系数的谱组合法,建立了基于位系数误差阶方差定权、误差方差定权和块对角误差协方差定权的谱组合计算公式。以GRACE重力场模型ITG-GRACE2010s和GOCE重力场模型GO_CONS_GCF_2_DIR_R1为例进行了谱组合法计算与验证,结果表明:三种谱组合方法建立的联合重力场模型精度均优于GRACE或GOCE单独求解的模型,它们的整体精度均优于由GOCE和GRACE联合求解的GOCO01S模型;基于块对角误差协方差定权的谱组合法相对最优,它考虑了同次m位系数之间的相关性。
从图 3可以看出,在上述三种谱组合法计算的重力场模型低阶部分,起主要贡献的是GRACE重力场模型;而中高阶部分,对谱组合计算模型起主要贡献的是GOCE重力场模型。同时,GRACE+GOCE-1和GRACE+GOCE-2的结果非常接近,大约在100阶至150阶,它们的精度略低于GOCO01S模型,并且GRACE+GOCE-2模型精度略优于GRACE+GOCE-1模型;而GRACE+GOCE-3在约130阶前和GOCO01S模型精度相当,在约130阶以后优于GOCO01S模型精度。可以看出,基于位系数块对角误差协方差进行谱组合求解模型的精度要高于其它两种谱组合方法的计算结果,且优于GOCO01S模型,验证了位系数谱组合方法的正确性和有效性。
图 3 GRACE和GOCE谱组合求解模型位系数阶方差
为了对谱组合计算模型的精度做进一步的分析,采用美国本土5377个GPS水准点数据进行外部检核。表格 2给出了GPS水准外部检核的统计结果,可以看出三种谱组合模型的精度不仅优于GRACE或GOCE数据单独求解的模型,并且优于GRACE和GOCE联合求解的GOCO01S模型,三种谱组合法相比GOCO01S的标准差小了大约4cm,相比单独解ITG-GRACE2010S的标准差小了12cm,相比GOCE的空域解、时域解和直接解的标准差分别小了大约6cm、4cm和0.5cm。同时,比较三种谱组合方法计算模型的精度可知,基于位系数块对角误差协方差定权的谱组合模型精度最优,这是因为它考虑了同次
m位系数之间的相关性;其次为基于位系数误差方差定权的谱组合模型,它考虑了单个位系数的误差特性;基于误差阶方差定权的谱组合模型精度最低,其原因是它将同阶位系数采用误差阶方差进行统一定权,忽略了同阶位系数之间精度的差异。
表格 2 GPS水准外部检核结果(单位:m)
地球重力场模型(阶数) |
最大值 |
最小值 |
平均值 |
STD |
ITG-GRACE2010S(180) |
3.628 |
-1.386 |
1.0296 |
0.6069 |
GO_CONS_GCF_2_SPW_R1(210) |
3.721 |
-1.116 |
1.0343 |
0.5433 |
GO_CONS_GCF_2_TIM_R1(224) |
3.640 |
-0.897 |
1.0382 |
0.5294 |
GO_CONS_GCF_2_DIR_R1(240) |
3.329 |
-0.601 |
1.0271 |
0.4883 |
GRACE+GOCE-1(240) |
3.296 |
-0.617 |
1.0214 |
0.4839 |
GRACE+GOCE-2(240) |
3.290 |
-0.615 |
1.0219 |
0.4838 |
GRACE+GOCE-3(240) |
3.302 |
-0.586 |
1.0223 |
0.4835 |
GOCO01S(224) |
3.585 |
-0.919 |
1.0394 |
0.5253 |
Ø
卫星重力场模型反演的动力积分法建立了动力积分法反演卫星重力场模型的理论模型和算法,研制了动力积分法计算软件WHUDynaSST。该软件的主要功能包括:基于IERS2003实现地固系与惯性系的转换;各种时间系统的相互转换;给出了现有的保守力计算模型,可完成卫星在任意时刻任意位置的保守力分析;给出了动力法加速度计校准模型,可实现加速度计的精密校准;利用Gauss-Jackson数值积分器实现了任意弧长、任意步长轨道积分的高效稳定算法,并能够完成任意阶次的变分方程解算;基于OpenMP实现了设计矩阵的构建、法方程的形成和解算等过程的并行算法,特别是针对最小二乘解算中的大规模矩阵求逆问题,基于MPI实现了高阶稠密对称正定矩阵的Gauss-Jordan并行求逆算法。利用GRACE卫星观测数据初步验证了动力积分法及其计算软件的正确性和有效性。
精度高、算法稳定、计算速度快的积分器在人造卫星动力学轨道确定、地球重力场模型解算等研究中有着十分重要的作用。针对卫星轨道数值积分、变分方程解算等问题,研究了Gauss-Jackson积分器的原理和计算流程,提出了移位重排方式来优化其存储方式,采用开普勒轨道、状态转移矩阵(表格3)等多种参数评估其性能,并与Runge-Kutta、Adams-Cowell等数值积分器进行了比较。计算结果表明:由于对启动点引入中值改正,Gauss-Jackson数值积分器的计算精度高、速度快,可为卫星轨道数值积分和变分方程求解等问题提供稳定、高效的算法。
表格 3 二体问题状态转移矩阵解析解与数值积分解的差值
二体问题状态转移矩阵Gauss-Jackson积分解与解析解残差(t=86400 s) |
0.24556D-10 |
0.37059D-12 |
0.11039D-11 |
0.52826D-08 |
-0.19878D-07 |
-0.39901D-07 |
-0.14196D-09 |
0.12218D-12 |
-0.90217D-12 |
-0.12972D-08 |
-0.82629D-07 |
-0.10249D-06 |
-0.17503D-09 |
-0.84832D-12 |
-0.27989D-12 |
-0.15989D-08 |
-0.10232D-06 |
-0.12554D-06 |
0.25039D-12 |
0.78540D-15 |
0.10278D-14 |
0.24787D-11 |
0.14334D-09 |
0.17403D-09 |
0.13989D-13 |
0.75463D-15 |
0.12798D-14 |
0.42775D-11 |
-0.14211D-10 |
-0.29502D-10 |
0.15293D-13 |
0.84112D-15 |
0.16535D-14 |
0.53597D-11 |
-0.18844D-10 |
-0.36522D-10 |
动力积分法是将卫星星历扰动作为地球重力场信息的泛函,通过建立轨道摄动与地球重力场位系数之间的关系,精密获取地球重力场模型的方法,具有反演精度高、获取高精度动力学轨道等优点,缺点是计算较为耗时,且对摄动力模型精度、加速度计测量和校准精度等极其敏感。针对最小二乘解算中的大规模矩阵求逆问题,基于MPI实现了高阶稠密对称正定矩阵的Gauss-Jordan并行求逆算法,减少了计算耗时(表4);并通过优化矩阵读写、存储等方式降低了单个计算节点的内存耗用量,拓展了算法的可移植性。针对加速度计精密校准问题,基于一步法,同时求解高精度的校准因子与位系数。为验证程序的可靠性,采用2009年3月1日至2009年3月31日的GRACE观测数据,获取了截断阶次为100的重力场模型,其与EIGEN-04C的大地水准面互差如图 4所示。
表格 4 不同进程数对应并行求逆的时间表
u |
进程数 |
1 |
8 |
16 |
24 |
32 |
6561 |
计算时间(s) |
2706.35 |
669.09 |
323.93 |
238.15 |
159.91 |
加速比 |
1.00 |
4.04 |
8.35 |
11.36 |
16.92 |
效率(%) |
100.00 |
50.56 |
52.21 |
47.35 |
52.89 |
14641 |
计算时间(s) |
12931.49 |
2712.32 |
1291.70 |
995.91 |
735.92 |
加速比 |
1.00 |
4.77 |
10.01 |
12.98 |
17.57 |
效率(%) |
100.00 |
59.60 |
62.57 |
54.10 |
54.91 |
图 4 GRACE卫星重力场模型的大地水准面误差