2. 上海市深空探测技术重点实验室,上海 200240
2. Shanghai Key Laboratory of Deep Space Technology,Shanghai,201109
相比于近地卫星,深空探测器飞行距离远、器地时延大、未知和不确定性因素多,存在天体遮挡与跟踪盲区,对导航的连续性、自主性等提出了更高的要求。为了确保未来深空探测重大工程任务的顺利实施,提高深空任务的成功率,降低工程技术风险,深空自主导航是必须且亟待突破的关键技术之一[1-2]。
当前深空探测自主导航技术大多通过天文测角或测距信息来实现导航状态的实时估计,存在不能直接获得速度测量信息、受目标天体观测条件约束等问题,在导航的连续性与实时性方面受限[3-6]。然而,天文光学信息中蕴含的光谱特征及其频移量包含了探测器的速度信息,基于此,本文提出了利用天文恒星光谱的直接测速导航方法,将其与测角导航相结合,形成组合导航系统,实现深空探测连续自主、实时高精度的导航[7]。
在深空测角测速组合导航系统多源信息融合过程中,要求各敏感器数据必须是统一时间基准,但在实际系统中,由于各敏感器时间基准误差、各敏感器采样周期不统一、数据传输延迟等都会造成时间不同步,而时间误差对位置和速度测量信息会带来很大的影响[8],因而,深空组合导航系统在多元信息融合前需要进行时间配准,以形成时间上统一的观测点,从而提高导航精度。本文分析了时间误差在深空测角测速组合导航系统位置估计和速度估计中的作用机理,研究了基于内插外推方法的时间配准方法,实现了测角敏感器与测速敏感器量测信息的同步。
1 测角测速组合导航系统原理 1.1 测角导航原理在巡航段末期,深空探测器与火星的距离越来越近,火星的视星等将随之逐步减小。星上测角导航敏感器通过测量获得火星、太阳的光学图像,进而提取火星、太阳相对于探测器的视线方向矢量,然后结合导航滤波算法实现自主导航。测角导航原理如图 1 所示。
根据图 1 中各视线矢量的几何关系,探测器在日心惯性坐标系中的位置可表示为
$\mathit{\boldsymbol{r}} = - \left\| {{\mathit{\boldsymbol{r}}_{\rm sm}}} \right\|\sqrt {\frac{{1 - {{\left( {{\mathit{\boldsymbol{l}}_{{\rm{sm}}}} \cdot {\mathit{\boldsymbol{l}}_{{\rm{pm}}}}} \right)}^2}}}{{1 - {{\left( {{\mathit{\boldsymbol{l}}_{{\rm{ps}}}} \cdot {\mathit{\boldsymbol{l}}_{{\rm{pm}}}}} \right)}^2}}}} \cdot {\mathit{\boldsymbol{l}}_{{\rm{ps}}}}$ | (1) |
![]() |
图 1 测角导航原理示意图 Fig. 1 Schematic diagram of angle measuring navigation |
其中:lps为太阳相对探测器的视线方向矢量;lpm为火星相对探测器的视线方向矢量,可通过导航敏感器测量获得;lsm为由火星相对太阳的方向矢量,可通过星历解算获得。lps和lpm观测量模型为
$\begin{array}{*{20}{c}}{\widetilde {{\mathit{\boldsymbol{l}}_{{\rm{ps}}}}} = {\mathit{\boldsymbol{l}}_{{\rm{ps}}}} + {\mathit{\boldsymbol{V}}_{{\rm{ps}}}}}\\{\widetilde {{\mathit{\boldsymbol{l}}_{{\rm{pm}}}}} = {\mathit{\boldsymbol{l}}_{{\rm{pm}}}} + {\mathit{\boldsymbol{V}}_{{\rm{pm}}}}}\end{array}$ | (2) |
其中:Vps、Vpm为视线矢量观测噪声。
取r的观测值为Z,则观测方程有
$\begin{aligned}& \mathit{\boldsymbol{Z}} = - \left\| {{\mathit{\boldsymbol{r}}_{{\rm{sm}}}}} \right\|\sqrt {\frac{{1 - {{\left( {{\mathit{\boldsymbol{l}}_{{\rm{sm}}}} \cdot \widetilde {{\mathit{\boldsymbol{l}}_{{\rm{pm}}}}}} \right)}^2}}}{{1 - {{\left( {\widetilde {{\mathit{\boldsymbol{l}}_{{\rm{ps}}}}} \cdot \widetilde {{\mathit{\boldsymbol{l}}_{{\mathop{\rm pm}\nolimits} }}}} \right)}^2}}}} \cdot \widetilde {{\mathit{\boldsymbol{l}}_{{\rm{ps}}}}} =\\& - \left\| {{\mathit{\boldsymbol{r}}_{{\rm{sm}}}}} \right\|\sqrt {\frac{{1 - {{\left[ {{\mathit{\boldsymbol{l}}_{{\rm{sm}}}} \cdot \left( {{\mathit{\boldsymbol{l}}_{{\rm{pm}}}} + {\mathit{\boldsymbol{V}}_{{\rm{pm}}}}} \right)} \right]}^2}}}{{1 - {{\left[ {\left( {{\mathit{\boldsymbol{l}}_{{\rm{ps}}}} + {\mathit{\boldsymbol{V}}_{{\rm{ps}}}}} \right) \cdot \left( {{\mathit{\boldsymbol{l}}_{{\rm{pm}}}} + {\mathit{\boldsymbol{V}}_{{\rm{pm}}}}} \right)} \right]}^2}}}} \cdot \left( {{\mathit{\boldsymbol{l}}_{{\rm{ps}}}} + {\mathit{\boldsymbol{V}}_{{\rm{ps}}}}} \right)\end{aligned} $ | (3) |
对式(3)在lps和lpm处泰勒展开得
$\mathit{\boldsymbol{Z}} = \mathit{\boldsymbol{r}} + \left( {\frac{{\partial \mathit{\boldsymbol{r}}}}{{\partial {\mathit{\boldsymbol{l}}_{{\rm{ps}}}}}}{\mathit{\boldsymbol{V}}_{{\rm{ps}}}} + \frac{{\partial \mathit{\boldsymbol{r}}}}{{\partial {\mathit{\boldsymbol{l}}_{{\rm{pm}}}}}}{\mathit{\boldsymbol{V}}_{{\rm{pm}}}}} \right) + o\left( \cdot \right)$ | (4) |
忽略观测量中的高阶小量,由参考天体视线方向矢量信息得到探测器位置矢量的观测方程可表示为
$\mathit{\boldsymbol{Z}} = h\left( \mathit{\boldsymbol{X}} \right) + \mathit{\boldsymbol{V}}\left( t \right)$ | (5) |
深空探测器在接收恒星光谱时,若深空探测器相对于恒星的位置是静止的,那么在探测器上接收的恒星光谱就是恒定的,若深空探测器相对于恒星的位置是变化的(视向方向接近或远离),那么所接收的光谱和相对静止时相比就会有波长的漂移,波长的漂移表现在光谱上就是谱线的移动[5-6]。根据多普勒原理,波长的漂移量与深空探测器相对于恒星位置静止时的波长之比等于视向速度与光速之比。
假设探测器、恒星在日心黄道惯性坐标系下的运动速度为v、vstar,探测器相对恒星视线方向为lstar,那么单颗恒星相对探测器的运动速度大小为
${\left( {{\mathit{\boldsymbol{v}}_{{\rm{star}}}} - \mathit{\boldsymbol{v}}} \right)^{\rm{T}}}{\mathit{\boldsymbol{l}}_{{\rm{star}}}} = {v_{{\rm{spe}}}}$ | (6) |
其中:vstar、lstar分别恒星在惯性参考坐标系下的速度矢量、视线方向矢量,可由星表获得。
![]() |
图 2 测速导航原理示意图 Fig. 2 Schematic diagram of velocity measuring navigation |
如图 2所示,建立以探测器相对太阳及其他两颗恒星的视向速度为量测的观测方案,经整理得到
$\left[ {\begin{array}{*{20}{c}}{\mathit{\boldsymbol{l}}_{{\rm{star1}}}^{\rm{T}}}\\[8pt]{\mathit{\boldsymbol{l}}_{{\rm{star2}}}^{\rm{T}}}\\[8pt]{\mathit{\boldsymbol{l}}_{{\rm{star3}}}^{\rm{T}}}\end{array}} \right]\mathit{\boldsymbol{v}} = \left[ {\begin{array}{*{20}{c}}{\mathit{\boldsymbol{l}}_{{\rm{star1}}}^{\rm{T}}{\mathit{\boldsymbol{v}}_{{\rm{star1}}}}}\\[8pt]{\mathit{\boldsymbol{l}}_{{\rm{star2}}}^{\rm{T}}{\mathit{\boldsymbol{v}}_{{\rm{star2}}}}}\\[8pt]{\mathit{\boldsymbol{l}}_{{\rm{star3}}}^{\rm{T}}{\mathit{\boldsymbol{v}}_{{\rm{star3}}}}}\end{array}} \right] - \left[ {\begin{array}{*{20}{c}}{{v_{{\rm{spe1}}}}}\\[8pt]{{v_{{\rm{spe2}}}}}\\[8pt]{{v_{{\rm{spe3}}}}}\end{array}} \right]$ | (7) |
分析式(7)可知,为了能够完全反演得到探测器在惯性参考坐标系下的速度矢量,v的系数矩阵[lstar1 lstar2 lstar3]必须可逆,即探测器相对3颗恒星的视线方向必须非共面。
综上,通过测量3颗及以上恒星的视向速度即可反演得到探测器在惯性参考坐标系下的速度矢量。
2 组合导航系统建模及时间配准算法 2.1 状态模型根据轨道动力学,建立深空探测器自主导航系统的状态模型。在火星探测巡航段,摄动因素主要考虑太阳、地球、火星等天体摄动,太阳光压摄动。设探测器在J2000日心黄道惯性坐标系下的位置矢量为
$\begin{array}{*{20}{c}}{\mathit{\boldsymbol{\dot r}} = \mathit{\boldsymbol{v}}}\\[8pt]{\mathit{\boldsymbol{\ddot r}} \!=\! - \displaystyle\frac{{{u_{\rm{s}}}}}{{{r^3}}}\mathit{\boldsymbol{r}} \!+\! \sum\limits_{i = 1}^N {{\mu _i}\left[ { - \frac{{{\mathit{\boldsymbol{r}}_i}}}{{r_i^3}} + \displaystyle\frac{{{\mathit{\boldsymbol{r}}_{si}}}}{{r_{si}^3}}} \right] \!+\! \eta {P_{{\rm{SR}}}}A{U^2}{C_{\rm{R}}}\left( {\displaystyle\frac{{{A_{\rm{R}}}}}{m}} \right)\displaystyle\frac{\mathit{\boldsymbol{r}}}{{{r^3}}}} \!+\! {\mathit{\boldsymbol{a}}_T}}\end{array}$ | (8) |
其中:第一项中μS为太阳引力常数,表示以太阳为引力中心的引力加速度项,第二项中ri是第i个行星在日心黄道惯性系下的位置矢量,μi是对应的行星引力常数,N = 1,2,3…表示行星的编号,rsi为行星相对探测器的位置矢量;第三项中η是阴影因子,PSR是距离太阳1 AU处的光压,CR为探测器的表面反射系数,AR是垂直于太阳光线方向探测器横截面积,m是探测器的质量是太阳光压摄动加速度项;最后一项表示推力加速度项。
综上,可以将探测器轨道动力学方程表述为如下一般形式
$\mathit{\boldsymbol{\dot X}}\left( t \right) = \mathit{\boldsymbol{f}}\left( {\mathit{\boldsymbol{X}}\left( t \right),t} \right) + \mathit{\boldsymbol{W}}\left( t \right)$ | (9) |
其中:W(t)表示状态模型噪声,一般作高斯白噪声处理。
2.2 量测模型 2.2.1 测角量测方程在火星巡航段末期及捕获段,可以利用观测太阳视线矢量和火星视线矢量作为导航的量测,即
${\mathit{\boldsymbol{Z}}_l} = \mathit{\boldsymbol{r}} + \left( {\frac{{\partial \mathit{\boldsymbol{r}}}}{{\partial {\mathit{\boldsymbol{l}}_{{\rm{ps}}}}}}{\mathit{\boldsymbol{V}}_{{\rm{ps}}}} + \frac{{\partial \mathit{\boldsymbol{r}}}}{{\partial {\mathit{\boldsymbol{l}}_{{\rm{pm}}}}}}{\mathit{\boldsymbol{V}}_{{\rm{pm}}}}} \right)$ | (10) |
其中:r是在参考坐标系下探测器的位置矢量;lps、lpm分别是太阳及火星视线方向矢量观测值;Vps、Vpm是量测噪声,一般作高斯白噪声处理。
双视线矢量观测可以间接得到探测器的位置和速度信息,是一种简单的测角导航方案,导航精度非常依赖敏感器的测量精度。
2.2.2 测速量测方程以探测器相对太阳及两颗恒星的视向速度大小为量测量,则量测方程可表示为
${\mathit{\boldsymbol{Z}}_v} = \left[ {\begin{array}{*{20}{c}}{{v_{{\rm{spe1}}}}}\\[8pt]{{v_{{\rm{spe2}}}}}\\[8pt]{{v_{{\rm{spe3}}}}}\end{array}} \right] + \mathit{\boldsymbol{V}} = \left[ {\begin{array}{*{20}{c}}{\displaystyle\frac{{\mathit{\boldsymbol{v}} \cdot \mathit{\boldsymbol{r}}}}{r}}\\[8pt]{\left( {{\mathit{\boldsymbol{v}}_{{\rm{star1}}}} - \mathit{\boldsymbol{v}}} \right) \cdot {\mathit{\boldsymbol{l}}_{{\rm{star1}}}}}\\[8pt]{\left( {{\mathit{\boldsymbol{v}}_{{\rm{star2}}}} - \mathit{\boldsymbol{v}}} \right) \cdot {\mathit{\boldsymbol{l}}_{{\rm{star2}}}}}\end{array}} \right] + \mathit{\boldsymbol{V}}$ | (11) |
其中:vspe1、vspe2、vspe3是由光谱仪解算得到的探测器相对恒星的视向速度值;r、v是探测器在参考惯性坐标系下的位置和速度;vStar、lStar是恒星在参考惯性坐标系下的速度和视线方向矢量,具体数据可由星表获得;V是量测噪声,一般作高斯白噪声处理。
2.3 时间配准方法时间配准指的是将各敏感器不同步的量测信息同步到统一基准时标下,并将不同步的信息配准到同一融合时刻。目前常用的时间配准算法有泰勒展开修正法、最小二乘法、内插外推法等[9-11]。
泰勒展开修正法通过对各敏感器在同一时间片内的量测数据进行泰勒展开,将高精度观测时间上的数据反演到低精度时间点上,其算法为:先取定时间片T;再将各敏感器观测数据按测量精度进行增量排序;最后将各高精度观测数据分别向最低精度时间点泰勒展开,从而形成一系列等间隔的目标观测数据以进行融合处理。泰勒展开修正法需要一阶导数,算法相对复杂。
最小二乘配准法采用最小二乘规则将第二类敏感器的n次测量值融合成一个虚拟的测量值作为时刻k第二类敏感器的测量值,然后同第一类敏感器的测量值进行融合,从而得到时刻k两敏感器测得的目标状态的融合值。需假定两类敏感器的采样周期之比n为整数。最小二乘法要求敏感器之间必须有1对n的严格对应关系,起始采样点也必须相同,而且采样精度是系统最低的敏感器采样率。
内插外推法采用在同一时间片内对各种敏感器采集的目标观测数据进行内插、外推,将高精度观测时间上的数据推算到低精度时间点上,其算法为:先取定时间片T;再将各敏感器观测数据按测量精度进行增量排序;最后将各高精度观测数据分别向最低精度时间点内插、外推,从而形成一系列等间隔的目标观测数据以进行融合处理。
鉴于测速敏感器可直接获得探测器速度信息,进而进行位置推算,考虑采用内插外推法进行深空探测测角测速组合导航系统时间配准。内插外推法以拉格朗日三点差值法为代表,假设tmi – 1、tmi、tmi + 1时刻测量数据为βi – 1、βi、βi + 1,则运用拉格朗日三点插值法计算出ti时刻的测量值为
$\begin{aligned}{{\bar \beta }_i} & = \displaystyle\frac{{\left( {{t_i} - {t_{mi}}} \right)\left( {{t_i} - {t_{mi + 1}}} \right)}}{{\left( {{t_{mi - 1}} - {t_{mi}}} \right)\left( {{t_{mi - 1}} - {t_{mi + 1}}} \right)}}{\beta _{i - 1}}+\\[5pt] & \displaystyle\frac{{\left( {{t_i} - {t_{mi - 1}}} \right)\left( {{t_i} - {t_{mi + 1}}} \right)}}{{\left( {{t_{mi}} - {t_{mi - 1}}} \right)\left( {{t_{mi}} - {t_{mi + 1}}} \right)}}{\beta _i}+\\[5pt]& \displaystyle\frac{{\left( {{t_i} - {t_{mi - 1}}} \right)\left( {{t_i} - {t_{mi}}} \right)}}{{\left( {{t_{mi + 1}} - {t_{mi - 1}}} \right)\left( {{t_{mi + 1}} - {t_{mi}}} \right)}}{\beta _{i + 1}}\end{aligned} $ | (12) |
深空天文测角测速组合导航系统结构如图 3所示。
![]() |
图 3 深空天文测角测速组合导航系统联邦滤波器框图 Fig. 3 Block diagram of federated filter of integrated navigation system based on celestial angle and velocity mesurement federated filter |
在Matlab环境下对组合导航系统进行仿真。系统模型状态方程中考虑太阳、火星及地球的引力摄动以及太阳光压摄动影响,仿真中的行星星历采用DE421;恒星数据由依巴谷星表提供,其中两颗恒星在星表中的编号为45 348和172 167。仿真起始时间为2021-1-15 05:46:07 UTC、终止时间为2021-1-25 18:40:00 UTC,真实轨道数据由STK产生,初始轨道误差Δr = [1 000 1 000 1 000] km、Δv = [100 100 100] m/s。假设测角敏感器采样率为60 s,测角精度为2 ′′,测速敏感器采样率为600 s,测速精度为1 m/s,测角敏感器的量测数据领先于测速敏感器。
图 4 给出了无时间配准情况下深空天文测角测速组合导航系统误差曲线。由图 4可知,位置误差呈现出一种无特征无规律的状态,整体发散。图 5 给出了采用内插外推时间配准算法下深空天文测角测速组合导航系统误差曲线。由图 5可知,使用时间配准算法后,组合导航系统收敛,内插外推时间配准算法有效抑制了时间误差,提高了组合导航系统的可靠性与稳定性。
![]() |
图 4 深空测角测速组合导航系统误差(无时间配准情况) Fig. 4 Error of integrated navigation system based on celestial angle and velocity mesurement(time free registration) |
![]() |
图 5 深空测角测速组合导航系统误差(有时间配准情况) Fig. 5 Error of integrated navigation system based on celestial angle and velocity mesurement(time registration) |
本文基于内插外推时间配准算法,将测角导航敏感器与测速导航敏感器量测数据统一到同一时间基准上,有效抑制了深空天文测角测速组合导航系统时间误差的影响,提高了深空测角测速组合导航系统的导航精度。仿真结果表明,该算法简单可行,在滤波精度、数据平稳性方面表现良好。
[1] |
张伟, 张恒. 天文导航在航天工程应用中的若干问题及进展[J]. 深空探测学报, 2016,3(3),204-213. Zhang W,Zhang H.Research on problems of celestial navigation in space engineering[J].Journal of Deep Space Exploration,2016,3(3),204-213. http://d.wanfangdata.com.cn/Periodical/sktcxb201603002 |
[2] |
陈晓, 尤伟, 黄庆龙. 火星探测巡航段天文自主导航方法研究[J]. 深空探测学报, 2016 ,3(3),214-218. Chen X,You W,Huang Q L.Research on celestial navigation for Mars missions during the interplanetary cruising[J].Journal of Deep Space Exploration,2016,3(3):214-218. http://d.wanfangdata.com.cn/Periodical/sktcxb201603003 |
[3] | Lightsey E G, Mogensen A, Burkhart P D. Real-time navigation for Mars missions using the Mars network[J]. Journal of Spacecraft and Rockets, 2008, 45(3): 519-533. doi: 10.2514/1.30974 |
[4] | Harlander J, Reynolds R J, Roesler F L. Spatial heterodyne spectroscopy for the exploration of diffuse interstellar emission lines at far-ultraviolet wavelengths[J][J]. Astrophysical Journal, 2013, 396(2): 730-740. |
[5] | Englert C R, Babcock D D, Harlander J M. Doppler asymmetric spatial heterodyne spectroscopy(DASH): concept and experimental demonstration[J]. Applied Optics, 2007, 46(29): 7297-7307. doi: 10.1364/AO.46.007297 |
[6] | 房建成, 宁晓琳. 深空探测器自主天文导航方法[M]. 西安: 西北工业大学出版社, 2010. |
[7] | 张伟, 陈晓, 尤伟. 光谱红移自主导航新方法[J]. 上海航天, 2013, 30(2): 32-33. Zhang W, Chen X, You W. New autonomous navigation method based on redshift[J]. Aerospace Shanghai, 2013, 30(2): 32-33. |
[8] |
袁赣南, 袁克非,等. 基于分段重叠的时间配准方法[J]. 敏感器与微系统, 2012, 31(8):48-56. Yuan G N,Yuan K F,et al.Time registration method based on segmentation overlap[J].Transducer and Microsystem Technologies, 2012, 31(8):48-56. http://d.wanfangdata.com.cn/Periodical/cgqjs201208015 |
[9] |
敬忠良, 肖刚, 李振华. 图像融合——理论与应用[M]. 北京: 高等教育出版社, 2010. Jing Z L, Xiao G, Li Z H.Image fusion: theory and applications[M].Beijing: Higher Education Press, 2010. |
[10] | Haykin S.Adaptive Filter Theory, Fourth Edition[M].[S. l.]: Pearson Education, 2006: 366-390. |
[11] | Yim J R, Crassidis J L, Junkins J L. Autonomous orbit navigation of interplanetary spacecraft[C]//AIAA/AAS Astrodynamics Specialist Conference. Reston, VA: AIAA, 2000: 53-61. |