应用引力梯度测量的火星中低轨道航天器自主导航

陈培 韩锦飞 赖玉敏 孙秀聪 谭龙玉

(1 北京航空航天大学宇航学院，北京 100191) (2 上海航天控制研究所，上海 201109) (3 上海空间智能控制技术重点实验室，上海 201109)

航天器的自主导航在航天任务中有着举足轻重的地位，对地球附近(地球中低轨道高度)和近地行星附近(近地行星中低轨道高度)的航天器导航都是如此：地球附近的航天器任务中，目前大多使用卫星定位(如GPS、北斗等)，其优点是定位精度高，且可以实现实时定位导航，但其也存在着电磁波信号易受到干扰和欺骗的问题[1]，尤其是军用航天器，为保证安全性与稳定性，其对自主导航的需求更为迫切；另外，近年来对于近地行星如火星等的探测任务逐渐增多，而目前GPS等卫星导航在其他行星附近处于完全不可用的状态，近地行星附近的探测器也需要一种更加有效且自主的导航方式。对于近地行星附近的航天器导航，目前存在的导航方式有限。常见的航天器自主导航方式有磁场导航、星光折射导航、X脉冲星导航、惯性导航等[2]。磁场导航利用地球附近每个点与地磁矢量存在的一一对应关系进行匹配导航[3]。其优点是利用了地球本身的磁场信息，可以实现实时测量解算。但地磁导航也存在着地磁场模型精度不够高的缺点，而且由于地球本身的磁场信号比较弱，测量过程中很容易受到周围磁场干扰。而且在一些行星，如金星、火星，据目前的资料显示，其上是没有磁场的，因此无法利用磁场对其附近地航天器进行导航。星光折射属于天文导航，容易受到天体可见性，以及大气模型的影响，而大气模型随着季节、维度、太阳活动而不断变化，因此利用星光折射法对地球附近的航天器进行导航精度并不高[4]；另外，由于近地行星周围没有像地球一样厚厚的大气层，因此，星光折射无法应用。X射线脉冲星导航能够提供一种独立的绝对时空基准，解决星座长时间自主运行的问题，为航天器自主导航提供了一种新的思路和可行途径。而其缺点也很明显，系统非常复杂，且只能实现相对惯性系的导航，相对于行星本体系的导航精度没有保证[5]；惯性导航的优点是在短时间内的导航精度很高，但由于众所周知的误差累积问题，长时间导航时会产生很大的导航误差，因此无法胜任长时间的自主导航任务[6]。以上导航方式都没办法对实现近地行星附近的导航，因此需要寻找一种能够完全自主且能够长时间在行星附近进行导航的方法。

近年来，世界航天强国也各自制定了相应的航天计划：美国的“空间探索计划”、俄罗斯的未来航天发展规划、欧空局制定的“曙光计划”与“欧洲阿波罗计划”等超大规模卫星探测计划等[7]。与此同时，我国在深空探测领域也取得了显著成就，对于月球以及火星的探索都展开了研究工作。但是，目前我国的航天探测任务主要以地面测控为主，有测控精度低、通讯时延长等缺点[8]。其对于地面设备的依赖使得其自主生存能力大大减弱。而且地面遥测遥控的方式无法满足交会飞越、下降着陆等任务阶段的实时性要求[9]，在此情况下，航天器的自主导航显得尤为重要。本文在上述研究基础上，针对目前近地行星附近的自主导航问题，提出了近地行星附近的引力梯度导航方法，并以火星为例，初步探索该方法应用于火星附近自主导航的可行性。

1 引力梯度匹配定位研究

随着大地测量技术的飞速发展，近地行星的引力场模型已经具有很高的精度，而且重力梯度仪也有很高的测量精度，欧空局在2009年已经将高精度的静电重力梯度仪配置于GOCE卫星上进行引力梯度测量。另外，引力梯度导航技术日趋成熟，并且已经成功应用于水下潜器和低空飞行器的惯性辅助导航系统中，并成功提高了惯性导航的精度[10-11]。但是在航天器导航上的应用却很少有人研究。笔者已经将全张量重力梯度应用到绕地飞行卫星导航系统。对于中心引力模型或只考虑J2项影响的模型，利用特征值分解可以直接获得估计位置，同时采用概率统计的方式解决了解的符号模糊问题[12]；对于考虑高阶项的EGM2008重力场模型，利用卡尔曼滤波可以获得位置的较为精确估计值[13]。以上定位方法均利用GOCE卫星实测数据进行了可行性验证，取得了成功。

目前对于火星引力场的研究，如文献[14]研究了目前两种最新的火星引力场模型JGMRO-120D和GMM-3，引力异常和星历积分两个方面的分析结果表明，对于10 m量级的定轨精度，这两个模型均可满足导航需求。由于JGMRO-120D模型可以直接从NASA官网下载，模型文件中含有火星引力场模型的时变修正公式，且此修正公式相对于GMM-3模型更为简单，使用也更加方便，因此本文采用JGMRO-120D模型，并在使用的过程中，考虑了火星引力场的时变特性，加入了时变修正项。通过仿真算例，实现利用引力梯度对航天器进行自主导航，为近地行星附近的自主导航方式的研究开辟一条新的道路。下面首先对引力梯度的基本概念加以阐释。引力梯度是引力加速度对位置向量的导数，其单位是E(Eötvös, 1E=10-9s-2)，也是引力势对位置向量的二阶导数。它是一个二维矩阵，称为引力梯度张量(gravity gradient tensor，GGT)矩阵。

(1)

式中：g是引力加速度向量；r是位置向量；(gx,gy,gz)和(x,y,z)分别是g和r在行星固连坐标系下对应方向上的分量。

火星周围空间各处引力场强度分布不同，它是由火星本身决定的。和火星相固连的每一点都有唯一的GGT，即位置和引力梯度张量是一一对应的。重力梯度仪并不能直接测量在火星固连坐标系下的引力梯度Vij，而只能测量在重力梯度仪本体坐标系(Gradiometer Reference Frame ，GRF)下的引力梯度。两个坐标系的定义分别如图1所示：O-XYZ表示火星固连坐标系，O为火星质心，OZ由地心指向协议北极点，OX由火星质心指向火星历元平赤道面与地球J2000历元平赤道面的交点，OY由右手空间直角坐标系可得。而GRF则与重力梯度仪安装角度有关，这里可认为其定义见图1，O′表示重力梯度仪中心，O′Z′沿着OO＇方向并远离火星质心，O′X′表示沿着轨道切线方向并与速度方向相同，而O′Y′则垂直于轨道平面[14]。以VM和VG分别表示火星固连坐标系和重力梯度仪本体坐标系下的GGT，则VM可以用下面的转换方式转换到VG。

(2)

图1 火星固连坐标系(MCF)和重力梯度仪坐标系(GRF)Fig.1 Mars-centered Mars-fixed reference frame and gravity gradiometer reference frame

由于GGT是对称的，所以重力梯度仪只需输出VG中的6个元素，因此在构造观测方程时，引入下面两个列向量

(3)

则观测方程可以写成

(4)

(5)

2 引力梯度匹配定位方法

与地球引力场模型相比，火星引力场模型的确定更为复杂，其中涉及到火星引力场季节性变化。火星非球形引力场模型与地球不同，其非球形引力位中的球谐项系数基本都要比地球相应值大一个量级，尤其是赤道椭率项，其大小接近于火星动力学扁率项 。由于火星引力场的这些特性，使得基于火星引力场的定位更加困难，因此，本文在建模时加入了对火星重力场模型时变量的修正。下面介绍一种将特征值分解和最小二乘迭代结合起来进行位置估计的方法。由于引力梯度是位置向量的非线性函数，故所要研究的问题是非线性问题。这里用最小二乘迭代的方法对位置进行求解。在非线性最小二乘问题中，初值的选取是比较关键的。在本方法中，初值用特征值分解法给出。在中心引力场模型中，将引力加速度对位置向量求一阶导数，得到中心引力场下GGT在行星固连坐标系下的表达式

(6)

式中：μ表示火星引力常数；r表示重力梯度仪的火心距，r=|r|；对求得的GGT进行特征值分解得到如下

(7)

Γ中对角线元素即为GGT的特征值，Φ中的3个列向量分别对应其特征向量。从Γ中可以看出，特征值只取决于r，无论其在坐标系中的任何方位，只要距离球心的距离相等，GGT分解出的特征值就相等；对应的特征向量只取决于φ和λ，二者分别代表了行星上的经度和纬度,见图1，无论其距离球心的距离如何变化，只要它所对应的经纬度没有发生变化，特征向量就不变。事实上，由球坐标系与笛卡尔坐标系之间的关系可知，Φ中的第三个特征向量对应着它们特征值中r的方向，是位置向量的单位向量，表示指向。即

r=rΦ3

(8)

式中：Φ3表示特征向量矩阵Φ的第3列。

用ξ表示特征值中的最大值，η为对应的特征向量。由式(7)可知

(9)

求得位置向量为

(10)

这里需要注意的是，特征向量存在方向上的不确定性，也就是说通过特征值分解法我们解算出的应该是两个关于球心对称的位置向量。这样在给定初值的时候就会出现两个位置初值，从而造成结果的歧义性。要消除这一问题，如果知道大致的初始位置，则可以采用比较法排除差距较大的值。

用特征值分解法给出位置初值后，采用非线性最小二乘法迭代求解位置，进一步提高位置解的精度。假设Q为重力梯度仪测量噪声方差矩阵。用非线性最小二乘迭代方法求解的迭代方程

(11)

(12)

高阶次引力场模型的Hk表达式比较复杂，这里不再详细列出。对于初值，Hk可由式(12)计算，在迭代过程中，Hk由高阶次的公式得到。位置估计的协方差矩阵可以由协方差公式得到

(13)

由于估计的协方差表示估计结果在各个方向上对真值的偏离程度，可以用其来描述最终定位的精度，同样，我们也可以用其平方根来表示精度。这里我们用其平方根，即使用标准差来表示位置估计的精度。方法流程如图2所示。

图2 行星附近引力梯度导航方法流程图

3 仿真分析

基于以上原理，在MATLAB软件中，以火星为例进行仿真验证。卫星轨道高度设定为距离火星表面300 km，轨道偏心率为0.05。动力学模型为NASA发布的最新的火星引力场JGMRO-120D模型，其最高阶数为120。火星重力场模型具有时变特性，在JGMRO-120D模型文件中，给出了模型的时变修正项

dJ3=2.958 705 874 08×10-9×sin [w×(t-t0)]

(14)

式中：dJ3表示火星重力场模型的J3项修正系数;w为火星自转的平均角速度，w=191.39 ()/a；(t-t0)为距离1999年1月1日的时间长度，单位为a。在仿真时已经加入了火星重力场的时变特性。考虑到计算速度，仿真采用前20阶球谐项系数，假定重力梯度仪固连坐标系与轨道坐标系重合。设定仿真采样频率为0.1 Hz，即采样间隔为10 s，一共采集6000个历元的数据。仿真分两次分别添加两种强度不同的噪声，一种与当前重力梯度仪的测量噪声相当，每个方向都设为正态分布的高斯白噪声，其标准差为0.1 E(100 mE)；另一种是标准差为10 mE的高斯白噪声。

当测量噪声为100 mE量级时，测量噪声协方差阵为

(15)

由特征值分解法可以得到位置解，由于实际解算时采用的是高阶引力场模型，因此需要将此位置解作为初始值利用最小二乘法进行迭代求解。将估计值与仿真的真值比较，可以得到在火星固连坐标系下位置估算的3个方向误差。对于估计方差，采用单历元的测量数据进行分析，每个历元有6个测量数据，迭代得到对应的方差和标准差，如图3所示。

图3 噪声100 mE时位置估算偏差及相应的标准差Fig.3 Position estimation error and corresponding standard deviation in noise 100mE

仿真结果表明，在300 km高度处，测量噪声在100 mE量级时，x方向估计误差约为154.3 m，y方向误差约为170.3 m，z方向估计误差约为139.3 m，三维误差约为268.9 m。

随着重力梯度仪测量精度的提高，测量噪声越来越小，未来有可能将测量噪声降低到10 mE以下。这里仿真得到测量噪声为10 mE时的结果：当测量噪声设置为10 mE时，用同样的方法得到火星固连坐标系下仿真定位结果，如图4所示。

图4 噪声10 mE时位置估算偏差及相应的标准差Fig.4 Position estimation error and corresponding standard deviation in noise 10mE

仿真结果表明，在测量噪声为10 mE量级时，x方向估计误差约为15.5 m，y方向误差约为16.7 m，z方向估计误差约为14.1 m，三维估计误差约为26.8 m，不同测量噪声水平下解算结果如表1所示。

表1 不同测量噪声水平下系统定位误差

从表1中分析可知，当重力梯度仪测量噪声水平在10 mE时，由此方法可以得到30 m以内的定位结果，当重力梯度仪测量噪声水平在100 mE时，定位误差也能控制在300 m以内。由于目前重力梯度仪的测量水平越来越高，美国马里兰大学的Paik团队甚至研制出了测量水平在1 mE量级的超导重力梯度仪[15]，由此可以看到，此种方法在已知高精度引力场的行星附近进行导航是完全可行的。

4 结束语

在近地行星附近，引力梯度导航是一种实用性非常强的导航方法。本文通过理论分析了引力梯度进行近地行星导航的可行性，给出了一种基于引力梯度的特征值分解导航算法，最后以火星为例进行了数值仿真。结果表明，在当前100 mE测量噪声水平，定位误差达到268 m左右，具有较高的精度；如果未来10 mE水平的重力梯度仪能够得到应用，则可以得到更加精确的导航结果。与其它方法相比，本方法是首次利用引力梯度对火星附近的航天器进行导航验证。方法整体思路简单，易于实现，且不用依赖任何地基设备即可对行星附近航天器进行长时间、全天候的实时自主导航，大大增加了此种导航方法在行星附近航天器上的应用价值。另外，文中采取的模型仅为较低阶的火星引力场模型，利用仿真数据验证方法的可行性与可能达到的定位精度。在实际航天任务中，可采用更高阶的火星引力场模型，以进一步提高定位精度。虽然本文以火星为例，但是此种方法适用于地球、类地行星及其它所有引力场接近有能力场且比较稳定的各类较大质量星体的中低轨道卫星的自主导航，应用前景十分广阔。

参考文献(References)

[1] 赵金磊. GPS定位及欺骗干扰技术[D]. 西安: 西安电子科技大学, 2014

Zhao Jinlei. GPS positioning and deception jamming technology[D]. Xi’an: Xi’an Electronic and Science University, 2014 (in Chinese)

[2] 姜竹青. 自主导航中滤波算法的研究及应用[D]. 北京：北京邮电大学，2014

Jiang Zuqing. Research and application of filtering algorithm in autonomous navigation[D]. Beijing: Beijing University of Posts and Telecommunications, 2014

[3] 寇义民. 地磁导航关键技术研究[D]. 哈尔滨：哈尔滨工业大学, 2010

Kou Yimin. Research on the key technology of geomagnetic navigation[D]. Harbin: Harbin Institute of Technology, 2010. (in Chinese)

[4] 宁晓琳, 王龙华, 白鑫贝, 等. 一种星光折射卫星自主导航系统方案设计[J]. 宇航学报, 2012, 33(11): 1601-1610

Ning Xiaolin, Wang Longhua, Bai Xinbei, et al. A scheme design of satellite autonomous navigation system based on stellar refraction[J]. Journal of Astronautics, 2012, 33(11): 1601-1610 (in Chinese)

[5] 孙景荣. X射线脉冲星导航及其增强方法研究[D]. 西安：西安电子科技大学, 2014

Sun Jingrong. X ray pulsar navigation and its enhancement method [D]. Xi’an: Xi’an Electronic and Science University, 2014 (in Chinese)

[6] 吴俊伟, 曾启明, 聂莉娟. 惯性导航系统的误差估计[J]. 中国惯性技术学报, 2002, 10(6):1-5

Wu Junwei， Tsang K M, Nie Lijuan. Estimation of the INS’s errors[J]. Journal of Chinese Inertial Technology, 2002, 10(6):1-5 (in Chinese)

[7] 韩鸿硕. 21世纪国外深空探测发展计划及进展[J]. 航天器工程, 2008, 17(3): 1-22

Han Hongsuo. 21stcentury foreign deep space exploration development plans and their progress [J]. Spacecraft Engineering, 2008, 17(3): 1-22 (in Chinese)

[8] 王大轶, 黄翔宇. 深空探测自主导航与控制技术综述[J]. 空间控制技术与应用, 2009, 35(03):6-12

Wang Dayi, Huang Xiangyu. Review of the autonomous navigation and control technology of deep space exploration[J]. Space Control and Technology and Application, 2009, 35(3): 6-12 (in Chinese)

[9] 吴伟仁, 于登云. 深空探测发展与未来关键技术[J]. 深空探测学报, 2014, 1(1): 5-17

Wu Weiren, Yu Dengyun. Development of deep space exploration and its future key technologies[J]. Journal of Deep Space Exploration, 2014, 1(1): 5-17 (in Chinese)

[10] 李兰玉. 基于全张量重力梯度的水下导航技术研究[D]. 武汉：武汉科技大学, 2014

Li Lanyu. Underwater navigation technology based on full tensor gravity gradient[D]. Wuhan:Wuhan University of Science and Technology, 2014 (in Chinese)

[11] 王文晶. 基于重力和环境特征的水下导航定位方法研究[D]. 哈尔滨工程大学, 2009

Wang Wenjing. Underwater navigation methods based on gravity and environmental features[D]. Harbin: Harbin Engineering University, 2009 (in Chinese)

[12] Chen P, Sun X, Han C. Gravity gradient tensor eigende composition for spacecraft positioning[J]. Journal of Guidance Control & Dynamics, 2015, 38: 2200-2206

[13] Sun X, Chen P, Macabiau C, et al. Autonomous orbit determination via Kalman filtering of gravity gradients[J]. IEEE Transactions on Aerospace and Electronic Systems, 2016, 52(5)

[14] 曹建峰, 刘磊, 黄勇, 等. 火星指向模型与重力场模型的发展回顾与使用[J]. 天文学进展, 2017, 35(1):127-139

Cao Jianfeng, Liu Lei, Huang Yong, et al. Review and the utilization of martian orientation model and gravity field model[J]. Progress in Astronomy, 2017, 35(1):127-139 (in Chinese)

[15] Moody M V, Paik H J, Canavan E R. Three-axis superconducting gravity gradiometer for sensitive gravity experiments[J]. Review of Scientific Instruments, 2002, 73(11): 3957-3974 (in Chinese)