卫星跟踪重力场探测系统误差分析与总体设计20070129

1、卫星跟踪重力场探测系统的误差分析与总体指标设计章传银常晓涛王正涛李娜（中国测绘科学研究院，100039，北京）摘要：卫星重力场探测系统的误差分析是卫星重力数据处理、重力场恢复算法设计、 卫星重力场探测性能分析和重力场探测系统总体设计的重要依据。本文从卫星跟踪引 力位双差重力场恢复数学模型出发，通过综合低低跟踪卫星重力场探测系统的主要观 测量误差、系统结构和重力场恢复算法，导出了系统从主要观测量到重力场位系数的 误差传播算法，从而构建卫星重力场探测系统的误差分析方法。在此基础上，提出了 低低跟踪卫星重力场探测系统的性能评价算法和总体指标设计算法。关键词：卫星跟踪卫星，地球重力场，误差分析，性能评

2、价，总体指标1引言地球重力场及其时变信息是地球科学、海洋学、测绘学、空间科学的基本量。利 用地球卫星探测地球重力场的主要方式有三种：卫星跟踪技术、卫星重力梯度测量和 海洋卫星测高技术，其中卫星跟踪卫星和卫星重力梯度测量技术被认为是21世纪初最 有价值和应用前景的重力场探测技术。德国、美国的CHAMP与GRACE卫星跟踪卫 星任务的实施，开辟了探测地球重力场及其变化的新纪元。卫星重力场探测系统的数学模型可一般描述为：已知地球重力场元g与重力场位 系数X之间关系的系统模型g = g &,r）（ 1.1）式中r为重力卫星的卫星位置向量。通过构建卫星重力场探测系统的观测量与地球重 力场元之间关系的观测

3、方程L = L（g（X , r ）, r, t）（1.2）从而采用合适的参数辨识方法估计重力场位系数X。在卫星跟踪卫星地球重力场探测系统中，一种重力场恢复的空域算法是：将（1.1） 式直接代入（1.2）式，按最小二乘法估计重力场位系数X。显然，重力场元g的类型决定了系统模型和观测方程的结构。为简化观测方程的 结构，有时将重力卫星观测量按某种函数形式进行组合，称组合后的统计量为模型观 测量，相应的观测方程为模型观测方程。不难看出，若要客观全面地分析卫星重力场探测系统的地球重力场探测能力，必 须同时顾及重力卫星观测量的误差、观测方程的结构和系统模型的结构。弄清卫星重 力场探测系统从主要观测量到重力

4、场位系数的误差传播特性，对有效地设计卫星重力 场探测系统、合理地构建卫星重力数据处理和重力场恢复算法、客观地评价重力场探 测能力均具有重要的指导意义。2卫星跟踪重力场恢复的数学模型基于能量守恒原理，在地球坐标系中，重力卫星的动能可表示为:(2.1 )=V + V (t) + w 2 (x2 + y2) + j f r dt + E一01_ r2 .2式中，V为稳态地球引力位，v （t）为时变引力位；f为非保守力向量（由星载加速度 计测得），w为地球自转角速率，r （x, y, z）表示重力卫星在地球坐标系中的位置向量， E为积分常量。2.1模型观测方程的构造构建如下的双星视线坐标系：x轴为卫星

5、b指向卫星a的视线方向，x轴在卫星 a、卫星b和地心构成的平面内，垂直于x轴，并指向地心方向（近似平行于向径方向）， x与x、x构成右手直角坐标系（x轴近似平行于轨道面的法线方向）。设视线坐标3123系在地球坐标系中的方向数为分别为l、l、l。123将重力卫星a、b的动能进行星间一次差分后，再对前后历元一次差分，得L = A V k，k+1 = (k+1 rk+1 )p k+1 (k rk 君 k A V k，k+1 (t) A j f dr -w 2 (xAx + yAy )k，k+1 + 6 r 1 Ar )+G r X Ar ),k+1(2.2)）,k+1式中（f -dr）,k+1为重力

6、卫星的非保守力在一个采样间隔时间内的沿轨道积分。式（2.2）即为以引力位双差aVk,k+1作为模型观测量，由卫星跟踪视线速率按能 量积分法恢复地球重力场的模型观测方程。2.2重力场位系数的最小二乘解在重力场位系数求解过程中，通常引入一个参考重力场，利用参考重力场位系数 按下式计算（2.2）中的a Vk e，作为模型观测量L的近似值l 0乙(C cos m 入 + S sin m X) P (sin 巾) m = 0（2.3）V r rn = 2式中r（r,X,。）为地球外空间任意点，u = GM为地球引力常数的位系数；P （sin巾）为完全规格化的Legendre缔合函数，P （sin巾）为其

7、导数。a和s为完全规格化将位系薮按某种方式排列成一列向量，令X为参考重力场位系数的改正数，则有 如下关于X的误差方程和法方程（2.4）V = L - L o = BX , （B t PB X = B t PV位系数改正数X的最小二乘解为（2.5）X = （B t PB ）1B t PV式中P为模型观测量的权阵，当仅利用一种跟踪模式的数据估计重力场位系数时，可 取单位阵。模型观测量的单位权中误差a 0可由下式计算(2.6)V t PVb = 0 V I - M式中I为模型观测量个数，M为位系数个数m = K 2 + 2 K - 3，K为最大位系数阶数。法方程系数阵N和位系数改正数的协因素阵Q分别

8、为N = B t PB，Q = t pb）1（ 2.7）则位系数改正数的中误差用下式表示b . = b Q,i = 1,2,. ,M（ 2.8）3低低跟踪卫星重力场的观测误差模型对于低低跟踪重力卫星，卫星轨道近似为圆轨道，则有（3.1）式并忽略加速度计在12（径 则有式中v为卫星沿轨方向的速率。将（3.1）式代入（2.8）向）、1 （轨道面法向）方向上的误差对A ,s 1的影响：(3.2)A V k, k +1 =v Ap k, k + 1 -A V k, k +1 （t） - A（ fds ）pVG . 1）+G .r .Ar）k,k+1 TOC o 1-5 h z r 233 J式中f为1

9、方向上的加速度计观测量，ds为重力卫星在一个历元间隔时间内的弧长。 13.1引力位双差的误差模型下面不考虑时变引力位模型误差，并忽略由于卫星轨道误差导致的视线坐标系方 向数的计算误差，分别导出视线速率误差b .、1方向上的加速度计测量误差b 、卫 星轨道误差和卫星速度误差对引力位双差A二,k+1的影响公式。f的加速度f求偏导差b f可表示为A V(3.3)令卫星a和卫星b在1方向上的加速度计观测误差相等，将（3.2）式对1方向上 则由于卫星加速度计测量误差导致的引力位双差a Vk, k+1的影响误b f = % 2 ds bA Vf式中为采样间隔，这里的采样间隔是指重力场恢复算法中采用的时间间

10、隔，而不是 主要载荷的采样间隔，下同。卫星速度误差对引力位双差AVk,k+1的影响误差可表示为b 勺牝 0， b vi2A VAV，b v3AVy)c 13 y 2 c %)(3.4)式中b 2、b %分别为轨道速度在地心和轨道面法向上的误差。vv(3.5)对于300600km高度的圆轨道低轨卫星，存在如下近似关系v化A v 5% 7 ,v牝A vIiiii2233 13 + ro 4式（3.10）描述了低低跟踪重力场探测系统的主要测量误差对引力位双差 “,5的总影响。右边第一项表示视线速率误差的影响，第二项表示加速度计测量误差的影响， 第三、四项分别表示卫星轨道向径方向和轨道面法向速度误差的

11、影响，第五项表示卫 星轨道在XY平面位置误差的影响。3.2引力位双差的误差影响特性由（3.10 ）可以看出，决定模型观测量-引力位双差误差的大小有两个因素：一 是视线速率、加速度计、卫星速度和卫星位置等测量误差，下称主要观测量误差；另 一个是卫星高度、视线距离、采样间隔等卫星重力场探测系统的结构参数，下称系统 的结构参数。所有主要观测量误差对引力位双差的影响都与视线距离p、卫星的地心 距,和采样间隔有关。通过比较不同观测量误差对模型观测量影响的权比，可以反映不同观测量对重力 场探测的贡献程度大小。由（3.3 ）和（3.9 ）式可知，当加速度计测量误差和视线速率 误差对模型观测量的影响相同时，下

12、式成立bb /（3.11）这说明，对于任意轨道高度、视线距离的低低跟踪卫星重力场探测系统，加速度 计的设计精度比视线速率设计精度高出两个数量级是最为合理的，超出这个范围并不 能提高重力场的探测能力。表1为按式（3.3） （3.7） （3.8 ） （3.9 ）计算出的不同观测量误差之间的权比关系。 例如，表1中第1行数据表示，当卫星轨道高度为450km、视线距离为150km、采样 间隔为10s时，1.4x1&6m/s的视线速率误差、2.3x10-8m/s2的加速度计测量误差、 3x10-5m/s的卫星速度误差和2.2m的卫星位置误差对双差引力位的影响是相同的。表1低低跟踪卫星主要观测量误差对引力

13、位双差影响的权比关系轨道高 度（km）视线距 离（km）采样间 隔（s）视线速率误 差（10-7m/s）加速度计测量误 差（10-8m/s2）速度误差 （10-5m）最大位置误 差（m）4501501013.92.313.02.19304.60.770.822201013.92.312.01.93304.60.770.812501501013.72.282.82.24304.60.760.852201013.72.281.92.97304.60.760.83由表1不难看出：（1）通过分析引力位双差的误差特性，可以确定重力卫星主要观测量之间的误差匹配关系，作为设计主要载荷技术指标的基本依据。（2

14、）不同结构参数的卫星跟踪重力场探测系统，主要观测量误差对引力位双差 影响的权比关系几乎相同。（3）当低低跟踪视线速率误差在10-6m/s时，加速度计测量精度高于10-8m/s2, 卫星速率精度高于10-5m/s,或者卫星轨道位置精度高于80cm，均不能改善重力场的 探测能力。同理，当加速度计的测量精度在10-8m/s时，若使低低跟踪视线速率精度高 于10-6m/s,卫星速率精度高于10-5m/s,或者卫星轨道位置精度高于80cm，均不能改 善重力场的探测能力。（4）由于高精度的低低跟踪视线速率的抵消作用，卫星在沿轨方向上的速度误 差影响可以忽略；卫星轨道位置在地球自转轴方向的误差影响为零，XY

15、面上的位置 误差影响相对于其他观测量的误差影响小得多。由（2.8）式可知，位系数的误差等于模型观测量的单位权中误差和位系数协因素 平方根的乘积，因此，引力位双差误差的大小不能直接反映重力场的探测能力，表1 中主要观测误差的列变化不能用于说明重力场的探测能力。例如，我们不能从表1中 得出“采样间隔越小重力场探测能力越大”这样的结论。4低低跟踪卫星重力场探测系统的结构特性由（2.8）式可知，为评价重力场的探测能力，除了需要求得模型观测量单位权中 误差，即引力位双差的误差外，必须计算位系数的协因素阵，两者综合在一起才能评 价重力场的探测能力。4.1位系数的协因素阵及其计算法方程系数阵N或位系数的协因

16、素阵Q描述了模型观测量与位系数之间的关系， 因此，分析法方程系数阵N或协因素阵Q结构，可以揭示重力场探测系统对地球重力 场位系数的灵敏度，位系数估计解的稳定性，以及模型观测量对位系数的贡献。这是 卫星重力场探测系统设计和重力场探测性能分析的重要技术手段。在低低跟踪卫星重力场恢复数学模型中，位系数的协因素矩阵完全由系统的结构 参数和模型观测量的类型决定，与模型观测量本身的大小及其误差的大小均无关。因 此，可以事先根据系统结构参数和模型观测量的类型来计算位系数的协因素阵。o- lu7 p r 设两颗同名重力卫星的轨道为圆轨道，卫星轨道倾角i为89，卫星质心的地心 距为两颗重力卫星间隔为p，采样间隔

17、为c。则重力卫星轨道角速度o$、轨道周 期八沿轨方向线速度7、采样间隔时间内重力卫星轨迹的弧长次、以及同一历元卫 星a与卫星b的角距差可用下式表示：（4.1）p、 = 2 arctg 2 r )卫星升交点赤径的变化率为Q = -7.2722 x 10 -4cos i（4.2）卫星轨道地心纬度和经度按下式计算。=arcsin(sin i sin 0 )(4.3)一兀(一兀 V 0一兀 /2)X = Xo + arctg (cos itg 0 ) 一 ( 一。)t + 0(一兀 / 2 0兀 / 2) 兀(兀 /2 0 n )式中0为t时刻卫星与升交点的角距，X0为升交点的经度。对于任意一种卫星跟

18、踪卫星重力场探测系统的结构参数（轨道高度、低低跟踪距 离、采样间隔、卫星运行时间等），可事先按（4.1）和（4.3）式计算两颗同名重力卫 星a、b的轨道，再按（2.7）式计算与之一一对应的法方程矩阵或协因素矩阵。取卫星轨道倾角i = 89，轨道高度为450km，视线距离为220km，采样间隔30 秒，卫星运行时间为30天，位系数协因素阵的对角线分量的平方根如图1所示。3.0E-007 2.9E-007 28E-007 2.7E-007 2.6E-007 2.5E-007 2.4E-007 2.3E-007 2 2E-OO7 2.1E-007 20E-007 1.9E-007 1 8E-007

19、1.7E-OO7 1 6E-007 1.5E-007 1.4E-007 1.3E-007 1.2E-007 1 1E-007 1 OE-0079 0E-008 8 0E-O08 7.0E-008 60E-008 5.0E-008 4.0E-008 3.0E-008 2.0E-008 1.0E-008图1位系数协因素阵对角线分量的平方根4.2地球重力场元的精度评定方法用下式计算： 2 = $n C 2nCnmm = 0不难看出，给定引力位双差误差，对于任意事先设计的系统结构参数，能计算一 一对应的地球重力场元的误差，包括位系数中误差和阶方差、重力异常误差和大地水 准面误差等。由（2.8）式计算位

20、系数的中误差。、气，则位系数误差的阶方差可 nmnm+ 2 ）（ 4.4）Snm.第n阶大地水准面及其重力异常的误差为（4.5） u u (n - 1)N R yA gR 2式中R为地球平均半径，Y为椭球面上的正常重力。累计到K阶的大地水准面和地面重力异常的误差为(4.6)=雳（ N），叮广声匕） =2l = 2比较重力场误差曲线和Kula曲线，可以用于评价恢复重力场位系数的最大可接 收的阶数。例如，德国空间中心（GFZ）综合利用2004年1月到2005年12月两年的地面 GPS跟踪站、GRACE卫星星载GPS跟踪数据、GRACE卫星之间低低跟踪数据、加 速度计数据、以及LAGEOS地面激光跟

21、踪数据等，将重力场位系数、GPS卫星轨道、 GRACE卫星轨道参数同时解算，获得150阶卫星重力场模型EIGEN-GL04S1，其位 系数的阶误差曲线如图2所示，图中的阶误差是指阶方差的平方根。图2卫星重力场模型EIGEN-GL04S1的位系数阶误差曲线由图2可以看出，GRACE重力卫星任务恢复重力场的最大可接收阶数为100阶。 EIGEN-GL04S1模型的大地水准面累计误差如图3。由图3可以看出，GRACE任务取 有效阶数约为90阶比较可靠。4.3根据协因素阵确定卫星运行最短时间长度随着重力卫星在轨运行时间的增大，协因素阵越趋完善，当在轨运行时间增加到 一定长度时，协因素阵结构改善的效果越

22、来越小，由于观测数据本身存在误差，因此， 当在轨运行时间到达某一长度时，再增大运行时间，难以进一步改善协因素阵的结构。 通常将这一时间长度成为卫星运行最短时间长度，它是重力卫星任务设计过程中的一 项重要技术指标。由（2.4）可以看出，即使重力卫星的在轨运行时间超过卫星运行最短时间长度， 位系数的协因素大小依然会随重力卫星运行时间长度的增加会不断变化，但若对位系 数的协因素进行归一化处理，则当重力卫星的在轨运行时间超过卫星运行最短时间长 度，归一化的位系数协因素不再随运行时间长度增加而发生明显变化。因此，归一化 的位系数协因素阵表示了卫星重力场探测系统的系统结构。参考（4.4 ）式，称下式为重力

23、场n阶谐分量归一化的协因素Q =乂 % ；气，式中 m = Y Q + Q,）（4.7）C nm nmm = 0n - 2 m - 0由（4.7 ）不难看出，卫星运行最短时间长度与所要求的重力场最大恢复阶数K 有关。对于给定的重力场探测系统和最大恢复阶数，当卫星运行时间长度增大到一定 数值时，各阶归一化的协因素不再发生较大变化，则对应的卫星运行时间长度，就是 卫星运行最短时间长度。例如，对于“卫星轨道倾角i = 89。，轨道高度为450km，低低跟踪视线距离为 220km，采样间隔30秒”的低低跟踪重力场探测系统，要求位系数的最大阶数为90阶。 不同卫星运行时间对应的位系数阶归一化协因素变化如

24、图4所示。4.4重力场探测系统的系统特性分析低低跟踪卫星重力场探测系统重力场恢复的最大阶数K与轨道高度、低低跟踪距 离存在如下近似关系：K -兀r / p（ 4.8）例如，GRACE 卫星，r 6828 km， p - 220 km,因此，k m 95。重力场的n阶位系数阶方差（注意与误差阶方差的区别），即重力场n阶谐分量 的平均信号强度或谱密度，通常用下式表示p = C 2 + S 2 ）（ 4.9）nmnnmm 0对于不同结构参数的重力场探测系统，即使模型观测量相同，其数值大小也不具 有可比性。但通过对谐分量的谱密度进行归一化处理，可使归一化后的谐分量谱密度 p与模型观测量大小无关。因此，

25、归一化的谐分量谱密度可用于评价不同结构参数的 重力场探测系统的性能。具体算法如下：对于模型观测量类型相同的任意两个重力场探测系统，令模型观测量是均值为 零、具有一定误差的统计量，按（2.5）式计算位系数的估值，并计算各阶归一化的谐 分量谱密度p。通过比较各阶谐分量的谱密度来分析那一种结构参数对应的系统对第 n阶位系数的探测能力更好。例如，取卫星轨道倾角i = 89。，分如下4种情况比较低低跟踪卫星重力场探测系 统的重力场探测性能：情况1：轨道高度h - 450km，视线距离p = 220km情况2：轨道高度h - 250km，视线距离p - 220km情况3：轨道高度h - 450km，视线距

26、离p = 500km情况4：轨道高度h - 250km，视线距离p - 500km按以上4种情况，分别计算直到90阶的归一化谐分量谱密度，结果如图6所示。图3四种系统的“阶误差”曲线由图3可以直接得出如下结论：5低低跟踪卫星重力场探测系统的性能评价与总体指标设计第3、4节阐述了低低跟踪卫星重力场探测系统从主要观测量误差到重力场位系 数的误差传播特性，本节以此为依据，进一步给出低低跟踪卫星重力场探测系统的性 能评价算法和总体指标设计算法。5.1低低跟踪卫星重力场探测系统的性能评价（1）卫星重力场探测性能的评价问题卫星重力场探测性能的评价问题是以重力场误差理论为基础的误差分析问题。对 于低低跟踪卫

27、星重力场探测系统，可做如下形式的一般描述：对于事先设计好的低低跟踪卫星重力场探测系统，即已知低低跟踪卫星重力场探 测系统的结构参数和主要测量指标，确定该重力场探测系统的重力场最大恢复阶数、 各阶位系数精度、相应阶次重力场元（大地水准面、地面重力异常）的精度。（2）低低跟踪卫星重力场探测性能的评价算法由主要测量指标按（3.10）式计算引力位双差的误差。由引力位双差误差和系统结构参数按4.2节计算位系数的中误差。由位系数中误差按（4.4）式计算位系数阶方差，绘制位系数阶误差曲线和Kula 曲线，确定重力场最大恢复阶数。由位系数阶误差和重力场最大恢复阶数按（4.5）式计算相应阶大地水准面和地 面重力

28、异常误差。例如，设卫星轨道倾角、89，轨道高度h = 450km，低低跟踪视线距离p = 220 km，采样间隔T =30s,卫星最短运行时间为45天；低低跟踪视线速率误差3x10-7m/s,重 力卫星加速度测量误差为5 x 10-9m/S2,卫星轨道速度误差为1 x 10-5m/s,卫星轨道位 置误差为0.3m。按（3.10 ）式求得引力位双差误差为2.76x10-3m2/s2, 按 4.3节计算位 系数的阶误差曲线和各阶大地水准面累计误差（单位：米）如图4所示。图4位系数的阶误差和各阶大地水准面累计误差曲线图由图4可知，重力场最大恢复阶数为76。虽然直到90阶的大地水准面误差还不 到0.0

29、1米，但不能认为重力场的恢复阶数大于76阶。但一般认为，取两者的最小值 为最大恢复阶数是合适的。图4隐含另一层信息，即当系统结构参数不变的情况下， 可通过适当提高主要观测量的精度，来改善重力场的探测能力。若要与其他已有重力场模型进行比较，也可通过计算已有重力场模型的位系数阶 误差，然后比较各阶位系数阶误差来考察。5.2低低跟踪卫星重力场探测系统的总体指标设计（1）重力场探测系统的总体指标设计问题重力场探测系统的总体指标设计问题是重力场误差分析的反问题。对于低低跟踪 卫星重力场探测系统，可作如下形式的一般描述：重力场探测能力的用户需求：最大重力场恢复阶数为K,相应分辨率的地面重力 异常误差为C K或大地水准面误差为Q K。确定低低跟踪卫星重力场探测系统的结构参数和主要测量指标。（2）低低跟踪卫星重力场探测系统的总体设计算法选择一组较为客观的且基本满足（4.8 ）式要求的轨道高度、低低跟踪距离和采 样间隔，轨道倾角小于并接近90 （倾角大小影响星下点轨迹，因此对确定卫星所需的 最短运行时间有影响），按4.1节计算位系数的协因