CN100533065C - 基于多天体路标的星际巡航自主导航方法 - Google Patents

Abstract

基于多天体路标的星际巡航自主导航方法，它涉及一种深空探测星际巡航中的自主轨道确定方法。本发明的技术方案如下：规划导航用小行星列表，处理导航小行星图像，基于多天体路标进行自主轨道确定，利用得到的多颗小行星的形心确定飞行器轨道，为了最小化舍入误差以及保证算法的数值稳定性，利用基于UD协方差分解的递推加权最小二乘算法来确定探测器的轨道。本发明的方法在没有测量误差和小行星星历误差的情况下的探测器轨道确定误差很小，基于多颗小行星图像的自主光学导航算法可以完全精确地确定探测器的轨道。在仿真假定条件下的自主轨道确定结果准确，位置误差接近100km，速度误差在0.3m/s范围内，可以满足探测器巡航段对轨道精度的要求。

Description

基于多天体路标的星际巡航自主导航方法

技术领域

本发明涉及一种深空探测星际巡航中的自主轨道确定方法。

背景技术

深空探测自主导航技术是航天科技发展的一项重要技术，其根本目的在于实现深空探测轨道确定的自主性，减少操作复杂性及降低任务花费。特别是在星际巡航阶段，自主导航技术尤为重要。自主导航按其获得信息的方法可分为三大类：1)惯性导航：连续地将测量的加速度或速度对时间积分得到位置信息；2)天文导航/GPS：直接由测量信息获得位置信息；3)光学导航/雷达图形匹配导航：处理测量信息获得导航信息。在现有的光学导航方法中，与星际巡航导航相关联的有以下几种：1)利用星敏感器进行自主导航和星光模拟半事物仿真；2)利用星敏感器星光折射或地平仪测量星光角距进行定轨；3)基于信息融合即利用直接敏感地平和星光折射敏感地平相结合，对卫星进行自主导航；4)基于小天体的图像信息的自主导航或日地月信息的自主导航方法。现有的光学导航方法随着时间的推进，取得了一定的进展，但还存在导航过程复杂、精度低、不易操作的缺点。

发明内容

本发明的目的是为解决现有的光学导航方法还存在导航过程复杂、精度低、不易操作的问题，提供一种基于多天体路标的星际巡航自主导航方法。本发明具有导航过程简单、易实现、精度高的优点。本发明的技术方案由以下步骤实现：第一、规划导航用小行星列表：

导航相机参数和导航小行星选择标准

与飞行器距离为0～1e6km，与飞行器的相对速度为0～7km/s，

绝对星等为0～12，太阳角为±40度，

规划小行星列表依据的是从美国JPL实验室得到的全部小行星的资料列表，包含有小行星的星历，星等的物理特性。根据小行星星历，采用二体轨道计算方法得到小行星当前时刻的位置和速度。小行星与飞行器的距离和相对速度由飞行器的标称轨道和小行星的位置速度得到。从而选取满足约束要求的星体。导航初始时刻，飞行器对日定向，且体坐标系与轨道坐标系重合。从而采用双矢量定姿法，期望姿态由飞行器位置和小行星位置得到。根据初始姿态和期望姿态，以及转动速度计算出姿态转换时间，选定时间最短的多颗小行星当作导航星，形成导航小行星列表。第二、处理导航小行星图像：其目的是为了得到小行星图像的形心。形心确定方法分两步：

首先确定“亮心”。

$p_{\mathrm{cb}}=\frac{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{ip}}_{\mathrm{ij}}}{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{p}_{\mathrm{ij}}},$ $l_{\mathrm{cb}}=\frac{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{jp}}_{\mathrm{ij}}}{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{p}_{\mathrm{ij}}}$

式中，p_ij是对应像素的亮度值(i是像元，j是像线)，p_cb和l_cb是计算的导航小行星图像亮心，m×n＝1024×1024是相机的分辨率。

然后根据亮心确定“形心”。

p₀＝p_cb-γR_ccosφ，l₀＝l_cb-γR_csinφ

式中， $\mathrm{\φ}=a\tan \left(\frac{A_{\mathrm{cy}}}{A_{\mathrm{cx}}}\right),$ A_c是相机坐标系太阳的方向，由指向太阳的惯性视线矢量A_I经过坐标转换得到；A_c＝T_CIA_I，其中T_CI为惯性坐标系到相机坐标系的转换矩阵。R_c是小行星半径对应的像素， $R_c=\frac{{\mathrm{RK}}_{x}f}{\mathrm{\ρ}},$ R是小行星的半径，ρ是探测器与小行星的距离，f是相机的焦距，K_x是从长度单位到像素的转换；γ是偏差因子，取值在0和1之间，表示目标半径误差ΔR引起的形心偏差，其具体表达式为

$\mathrm{\γ}=\frac{3\mathrm{\π\ΔR}}{16}\left[\frac{\sin \mathrm{\α}(1+\cos \mathrm{\α})}{(\mathrm{\π}-\mathrm{\α})\cos \mathrm{\α}+\sin \mathrm{\α}}\right]$

α是太阳相位角。

最后进行形心像素的计算，

惯性视线矢量

，利用转换矩阵T_CI旋转到照相机坐标系，

${\stackrel{\→}{V}}_C=\left[\begin{array}{c}{V}_{C1}\\ V_{C2}\\ V_{C3}\end{array}\right]={\mathrm{TC}}_I{\stackrel{\→}{V}}_I$

视线矢量

确定，转换到二维的照相机焦平面中，

$\left[\begin{array}{c}x\\ y\end{array}\right]=\frac{f}{V_{C3}}\left[\begin{array}{c}{V}_{C1}\\ V_{C2}\end{array}\right]$

其中，

f 照相机焦距mm，

V_C1，V_C2，V_C3  视线矢量在相机系内分量，

x，y  视线矢量焦平面内投影mm，

接着，由光的畸变引起的x和y偏差Δx，Δy，

$\left[\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\end{array}\right]=Q{\left[\begin{array}{cccccc}{v}_1& v_2& v_3& v_4& v_5& v_6\end{array}\right]}^T$

其中，

$Q=\left[\begin{array}{cccccc}-\mathrm{yr}& {\mathrm{xr}}^2& {-\mathrm{yr}}^3& {\mathrm{xr}}^4& \mathrm{xy}& x^2\\ \mathrm{xr}& {\mathrm{yr}}^2& {\mathrm{xr}}^3& {\mathrm{yr}}^4& y^2& \mathrm{xy}\end{array}\right]$

r＝x²+y²，v′s是光畸变系数。修正的图像位置x′，y′表示为

$\left[\begin{array}{c}x\′\\ y\′\end{array}\right]=\left[\begin{array}{c}x+\mathrm{\Δx}\\ y+\mathrm{\Δy}\end{array}\right]$

最后，直角坐标系转化到像元和像线，

$\left[\begin{array}{c}p\\ l\end{array}\right]=\left[\begin{array}{ccc}{K}_x& K_{\mathrm{xy}}& K_{\mathrm{xxy}}\\ K_{\mathrm{yx}}& K_y& K_{\mathrm{yxy}}\end{array}\right]\left[\begin{array}{c}x\′\\ y\′\\ x\′y\′\end{array}\right]+\left[\begin{array}{c}{p}_0\\ l_0\end{array}\right]$

K阵是从毫米到像元/像线的转换矩阵，p₀，l₀是CCD的中心像元和像线。

第三、基于多颗小行星进行自主轨道确定：利用上述得到的多颗小行星的形心确定飞行器轨道。

1)定义参考轨道参数：

$X^{*}\left(t\right)={\left[\begin{array}{cccccc}x& y& z& \stackrel{\·}{x}& \stackrel{\·}{y}& \stackrel{\·}{z}\end{array}\right]}^T,$ 其中，r＝[x y z]^T和 $\stackrel{\‾}{v}={\left[\begin{array}{ccc}\stackrel{\·}{x}& \stackrel{\·}{y}& \stackrel{\·}{z}\end{array}\right]}^T$ 分别为探测器在日心黄道坐标系的位置和速度矢量；

更新的轨道参数：

X′(t)＝X^*(t)+ΔX(t)，其中，ΔX(t)为估计轨道修正量。由于参考轨道参数和真实轨道相差不大，那么在一段时间内轨道修正量是线性的，即在一时间点的轨道参数修正量ΔX(t₀)利用状态转换矩阵线性地映射到任何其它时间点t上，即ΔX(t)＝Φ(t)ΔX(t₀)

其中，给定时间t的状态转移矩阵Φ满足

$\stackrel{\·}{\mathrm{\Φ}}=\frac{\∂\stackrel{\·}{X}\left(t\right)}{\∂X\left(t\right)}\frac{\∂X\left(t\right)}{\∂X\left(t_0\right)}=\mathrm{A\Φ}\left(t\right)$

式中， $A=\frac{\∂\stackrel{\·}{X}\left(t\right)}{\∂X\left(t\right)},$ Φ(t₀)＝I_6×6。通过这些方程的积分来获得Φ的时间函数。

***的状态方程写为：

$\stackrel{\·}{X}=\left[\begin{array}{c}\stackrel{\‾}{v}\\ -\frac{{\mathrm{\μ}}_s}{r^3}\stackrel{\‾}{r}+\underset{i=1}{\overset{n_p}{\mathrm{\Σ}}}{\mathrm{\μ}}_i[\frac{{\stackrel{\‾}{r}}_{\mathrm{ri}}}{r_{\mathrm{ri}}^3}-\frac{{\stackrel{\‾}{r}}_{p1}}{r_{p1}^3}]-\frac{\mathrm{AG}}{{\mathrm{mr}}^3}\stackrel{\‾}{r}+\frac{k}{m}\stackrel{\‾}{T}+\stackrel{\‾}{a}\end{array}\right]$

2)确定观测矩阵：

对于某一观测历元，观测矩阵为

$H=\left[\begin{array}{cccccc}\∂p/{\∂x}_1& \∂p/{\∂x}_2& \∂p/{\∂x}_3& 0& 0& 0\\ \∂l/{\∂x}_1& \∂l/{\∂x}_2& \∂l/{\∂x}_3& 0& 0& 0\end{array}\right]$

式中，与速度分量有关的偏微分都是零，这是由于p，l仅与拍照瞬时探测器相对小行星的位置有关，而与探测器的速度无关。

3)统一观测历元：

由于得到的多颗小行星像元、像素是不同时刻的，所以需要利用状态转移矩阵把各个观测历元的观测矩阵变换到同一观测历元，有

$\tilde{H}=\mathrm{H\Φ}$

其中就是在某时刻的观测偏微分矩阵。

4)实现轨道确定：

为了最小化舍入误差以及保证算法的数值稳定性，利用基于UD协方差分解的递推加权最小二乘算法来确定探测器的轨道。其中的测量残差Y是图像处理得到的与计算预测得到的小行星中心对应的像素差，这样就利用多颗小行星的观测数据(小行星中心的像素)，通过多步递推来修正某观测历元的轨道参数。首先给出一个协方差矩阵P₀，观测加权矩阵W，矢量Y为观测质心与标称轨迹预报计算质心的差值。初始点状态批量滤波方程求解矢量

和正式协方差矩阵P：

$\hat{q}={[P_0^{-1}+{\tilde{H}}^{T}W\tilde{H}]}^{-1}{\tilde{H}}^T\mathrm{WY}$

$P={[P_0^{-1}+{\tilde{H}}^{T}W\tilde{H}]}^{-1}$

其中，

$W=\left[\begin{array}{cc}1/{\mathrm{\σ}}_0^2& 0\\ 0& 1/{\mathrm{\σ}}_0^2\end{array}\right]$

${\mathrm{\σ}}_0=\frac{{\tan }^{-1}(R/\mathrm{\ρ})}{13\×{10}^{-6}}$

R小行星假定半径ρ小行星范围。

本发明的有益效果：

为验证本发明的效果，以探测Ivar1627小行星的一段巡航段轨道为例，进行数学仿真，仿真参数：

1)探测器的标称轨道利用数值积分得到。在J2000.0日心黄道惯性坐标系，

探测器的初始位置[1.13984×10¹¹-1.14516×10¹¹-6.73821×10⁶]m，

探测器的初始速度为[2.71473×10⁴1.883397×10⁴-2.59150]m/s；

2)测量精度：姿态误差方差为10^-12rad²，像素误差为：像元0.1，像线0.1；

3)测量频率：测量数据以450秒采样间隔输出，利用选择和规划的小行星图像进行轨道确定；

4)初始误差：探测器位置在各方向误差的方差为1×10¹⁴m²，各方向速度误差方差为10⁴m²/s²；

5)星历误差：小行星的星历在各个方向误差的方差为1×10¹⁰m²。

参考轨迹状态参数初值为

X_参考＝X_标称+[1×10⁷ 1×10⁷ 1×10⁷ -10 10 10]′

协方差阵

$P_0=\left[\begin{array}{cccccc}{10}^{14}& 0& 0& 0& 0& 0\\ 0& {10}^{14}& 0& 0& 0& 0\\ 0& 0& {10}^{14}& 0& 0& 0\\ 0& 0& 0& 100& 0& 0\\ 0& 0& 0& 0& 100& 0\\ 0& 0& 0& 0& 0& 100\end{array}\right]$

不确定度和观测加权矩阵

$\begin{array}{cc}{\mathrm{\σ}}_0=0.001& W=\left[\begin{array}{cc}1/{0.001}^2& 0\\ 0& 1/{0.001}^2\end{array}\right]\end{array}$

m＝500Kg；k＝1；A＝50(m²)；G＝4.65e-6(N/m²)

在没有测量误差和小行星星历误差的情况下，观测到小行星个数与轨道估计误差之间的关系如表1所示。可以看出，随着观测到小行星个数的增加，轨道估计精度逐步提高；对于能观测到7个小行星时，4次轨道确定后，位置精度能达到4.3m，速度精度0.0013m/s，非常接近轨道真值；对于能观测到8个小行星时，2次轨道确定后就非常接近轨道真值；但对于能观测到6个以下小行星时，轨道确定收敛速度慢，且不能精度确定轨道。这些表明，至少需要观测到7个小行星才能进行轨道确定。

表1 观测小行星个数与轨道估计误差的关系

在没有测量误差和小行星星历误差的情况下的探测器轨道确定误差如附图2所示，可以看出，基于多颗小行星图像的自主光学导航算法可以完全精确地确定探测器的轨道。在仿真假定条件下的自主轨道确定结果如附图3所示，可以看出，位置误差接近100km，速度误差在0.3m/s范围内，可以满足探测器巡航段对轨道精度的要求。

附图说明

图1是相机拍照的几何约束坐标图，图2是无测量误差和星历误差情况下的轨道确定误差曲线图，图3是利用多颗小行星图像的轨道确定误差曲线图。

具体实施方式

具体实施方式一：(参见图1～图3)本实施方式的步骤如下：

第一、规划导航用小行星列表：

1)导航星选择标准，

导航相机参数和导航小行星选择标准

 

选择标准	范围
		与飞行器距离	0～1000000km
与飞行器的相对速度	0～7km/s
		绝对星等	0～12
太阳角	±±40度

2)规划小行星列表第一步：从美国JPL实验室可以得到小行星列表，包含有小行星的星历，星等等物理特性。采用二体轨道计算方法得到小行星当前的位置和速度。与飞行器的距离和相对速度可以由标称轨道和小行星的星历得到。

3)规划小行星列表第二步：飞行器对日定向，且体坐标系与轨道坐标系重合(见附图1)。采用双矢量定姿法，期望姿态由飞行器位置和小行星位置得到。初始姿态和期望姿态计算出转换时间，选定时间最短的12颗；

第二、处理导航小行星图像：

其目的是为了得到小行星图像形心。形心确定方法分两步：

第一步，确定“亮心”。

$p_{\mathrm{cb}}=\frac{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{ip}}_{\mathrm{ij}}}{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{p}_{\mathrm{ij}}},$ $l_{\mathrm{cb}}=\frac{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{jp}}_{\mathrm{ij}}}{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{p}_{\mathrm{ij}}}$

式中，p_ij是对应像素的亮度值(i是像元，j是像线)，p_cb和l_cb是计算的导航小行星图像亮心，m×n＝1024×1024是相机的分辨率。

第二步，确定“形心”。

p₀＝p_cb-γR_ccosφ，l₀＝l_cb-γR_csinφ

式中， $\mathrm{\φ}=a\tan \left(\frac{A_{\mathrm{cy}}}{A_{\mathrm{cx}}}\right),$ A_c是相机坐标系太阳的方向，可由指向太阳的惯性视线矢量A_I经过坐标转换得到；A_c＝T_CIA_I，其中T_CI为惯性坐标系到相机坐标系的转换矩阵。R_c是小行星半径对应的像素， $R_c=\frac{{\mathrm{RK}}_{x}f}{\mathrm{\ρ}},$ R是小行星的半径，ρ是探测器与小行星的距离，f是相机的焦距，K_x是从长度单位到像素的转换；γ是偏差因子，取值在0和1之间，表示目标半径误差ΔR引起的形心偏差，其具体表达式为

$\mathrm{\γ}=\frac{3\mathrm{\π\ΔR}}{16}\left[\frac{\sin \mathrm{\α}(1+\cos \mathrm{\α})}{(\mathrm{\π}-\mathrm{\α})\cos \mathrm{\α}+\sin \mathrm{\α}}\right]$

α是太阳相位角。

第三步，像素的计算，

如果得到了一个惯性视线矢量

利用转换矩阵T_CI旋转到照相机坐标系， ${\stackrel{\→}{V}}_C$

${\stackrel{\→}{V}}_C=\left[\begin{array}{c}{V}_{C1}\\ V_{C2}\\ V_{C3}\end{array}\right]={\mathrm{TC}}_I{\stackrel{\→}{V}}_I$

视线矢量

一旦确定，需要转换到二维的照相机焦平面中，

$\left[\begin{array}{c}x\\ y\end{array}\right]=\frac{f}{V_{C3}}\left[\begin{array}{c}{V}_{C1}\\ V_{C2}\end{array}\right]$

其中，

f 照相机焦距mm，

V_C1，V_C2，V_C3 视线矢量在相机系内分量，

x，y 视线矢量焦平面内投影mm，

接着，由光的畸变引起的x和y偏差Δx，Δy，

$\left[\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\end{array}\right]=Q{\left[\begin{array}{cccccc}{v}_1& v_2& v_3& v_4& v_5& v_6\end{array}\right]}^T$

其中，

$Q=\left[\begin{array}{cccccc}-\mathrm{yr}& {\mathrm{xr}}^2& {-\mathrm{yr}}^3& {\mathrm{xr}}^4& \mathrm{xy}& x^2\\ \mathrm{xr}& {\mathrm{yr}}^2& {\mathrm{xr}}^3& {\mathrm{yr}}^4& y^2& \mathrm{xy}\end{array}\right]$

r＝x²+y²，v′s是光畸变系数。修正的图像位置x′，y′表示为：

$\left[\begin{array}{c}x\′\\ y\′\end{array}\right]=\left[\begin{array}{c}x+\mathrm{\Δx}\\ y+\mathrm{\Δy}\end{array}\right]$

最后，直角坐标系转化到像元和像线，

$\left[\begin{array}{c}p\\ l\end{array}\right]=\left[\begin{array}{ccc}{K}_x& K_{\mathrm{xy}}& K_{\mathrm{xxy}}\\ K_{\mathrm{yx}}& K_y& K_{\mathrm{yxy}}\end{array}\right]\left[\begin{array}{c}x\′\\ y\′\\ x\′y\′\end{array}\right]+\left[\begin{array}{c}{p}_0\\ l_0\end{array}\right]$

K阵是从毫米到像元/像线的转换矩阵，p₀，l₀是CCD的中心像元和像线。

第三、基于多天体路标(多颗小行星)进行自主轨道确定：利用上述得到的多颗小行星的形心及背景恒星确定的姿态，确定轨道。

1)定义参考轨道参数：

$X^{*}\left(t\right)={\left[\begin{array}{cccccc}x& y& z& \stackrel{\·}{x}& \stackrel{\·}{y}& \stackrel{\·}{z}\end{array}\right]}^T,$ 其中，r＝[x y z]^T和 $\stackrel{\‾}{v}={\left[\begin{array}{ccc}\stackrel{\·}{x}& \stackrel{\·}{y}& \stackrel{\·}{z}\end{array}\right]}^T$ 分别为探测器在日心黄道坐标系的位置和速度矢量；

更新的轨道参数：

X′(t)＝X^*(t)+ΔX(t)，其中，ΔX(t)为估计轨道修正量。如果参考轨道参数和真实轨道相差不大，那么在一段时间内轨道修正量是线性的，即在一时间点的轨道参数修正量ΔX(t₀)可以利用状态转换矩阵线性地映射到任何其它时间点t上，即ΔX(t)＝Φ(t)ΔX(t₀)

其中，给定时间t的状态转移矩阵Φ满足

$\stackrel{\·}{\mathrm{\Φ}}=\frac{\∂\stackrel{\·}{X}\left(t\right)}{\∂X\left(t\right)}\frac{\∂X\left(t\right)}{\∂X\left(t_0\right)}=\mathrm{A\Φ}\left(t\right)$

式中， $A=\frac{\∂\stackrel{\·}{X}\left(t\right)}{\∂X\left(t\right)},$ Φ(t₀)＝I_6×6。***的状态方程可写为：

$\stackrel{\·}{X}=\left[\begin{array}{c}\stackrel{\‾}{v}\\ -\frac{{\mathrm{\μ}}_s}{r^3}\stackrel{\‾}{r}+\underset{i=1}{\overset{n_p}{\mathrm{\Σ}}}{\mathrm{\μ}}_i[\frac{{\stackrel{\‾}{r}}_{\mathrm{ri}}}{r_{\mathrm{ri}}^3}-\frac{{\stackrel{\‾}{r}}_{p1}}{r_{p1}^3}]-\frac{\mathrm{AG}}{{\mathrm{mr}}^3}\stackrel{\‾}{r}+\frac{k}{m}\stackrel{\‾}{T}+\stackrel{\‾}{a}\end{array}\right]$

式中：

r探测器在日心坐标系内位置矢量；v探测器在日心坐标系内速度矢量；r_pi第i个摄动行星体在日心坐标系内的位置矢量，仿真中考虑了，水、金、地、火、木五大行星以及月球的摄动影响，r_ri探测器相对第i个摄动行星体的位置矢量，其中r_ri＝r_pi-r；μ_δ太阳引力常数GM；μ_i第i个摄动天体的引力常数；n_p摄动天体的数量；A探测器的有效面积；G日光流量常数；T推进***推力矢量；k推力比例系数，近似为1；m探测器质量；a探测器总的附加加速度；

在上式中，右边第一项是右中心体太阳引力引起的加速度；第二项是第三体引力加速度的总和；第三项是太阳光压；第四项是推进***的推进加速度；最后一项表示作用在探测器上的其他附加加速度。前两项引力加速度的作用是直接的，但是作用在探测器上的非引力因素是值得讨论的。对于太阳光压，式中对于探测器采用了简单的球形模型，但在实际中，探测器的特征面积主要是由两个太阳帆板决定的，而探测器主体只占非常小的一部分。

2)确定观测矩阵：

对于某一观测历元，观测矩阵为

$H=\left[\begin{array}{cccccc}\∂p/{\∂x}_1& \∂p/{\∂x}_2& \∂p/{\∂x}_3& 0& 0& 0\\ \∂l/{\∂x}_1& \∂l/{\∂x}_2& \∂l/{\∂x}_3& 0& 0& 0\end{array}\right]$

式中，与速度分量有关的偏微分都是零，这是由于p，l仅与拍照瞬时探测器相对小行星的位置有关，而与探测器的速度无关。

3)统一观测历元：

由于得到的多颗小行星像元、像素是不同时刻的，所以需要利用状态转移矩阵把各个观测历元的观测矩阵变换到同一观测历元，有

$\tilde{H}=\mathrm{H\Φ}$

4)实现轨道确定：

为了最小化舍入误差以及保证算法的数值稳定性，利用基于UD协方差分解的递推加权最小二乘算法来确定探测器的轨道。其中的测量残差Y是图像处理得到的与计算预测得到的小行星中心对应的像素差，这样就可以利用多颗小行星的观测数据(小行星中心的像素)，通过多步递推来修正某观测历元的轨道参数。首先给出一个协方差矩阵P₀，观测加权矩阵W，矢量Y为观测质心与标称轨迹预报计算质心的差值。初始点状态批量滤波方程求解矢量

和正式协方差矩阵P：

$\hat{q}={[P_0^{-1}+{\tilde{H}}^{T}W\tilde{H}]}^{-1}{\tilde{H}}^T\mathrm{WY}$

$P={[P_0^{-1}+{\tilde{H}}^{T}W\tilde{H}]}^{-1}$

其中，

$W=\left[\begin{array}{cc}1/{\mathrm{\σ}}_0^2& 0\\ 0& 1/{\mathrm{\σ}}_0^2\end{array}\right]$

${\mathrm{\σ}}_0=\frac{{\tan }^{-1}(R/\mathrm{\ρ})}{13\×{10}^{-6}}$

R小行星假定半径ρ小行星范围。

Claims (1)

1、一种基于多天体路标的星际巡航自主导航方法，基于多天体路标的星际巡航自主导航方法的技术方案由以下步骤实现：第一、规划导航用小行星列表；

导航相机参数和导航小行星选择标准：与飞行器距离为0～1e6km，与飞行器的相对速度为0～7km/s，绝对星等为0～12，太阳角为±40度；

规划小行星列表依据的是从美国JPL实验室得到的全部小行星的资料列表，包含有小行星的星历，星等的物理特性，根据小行星星历，采用二体轨道计算方法得到小行星当前时刻的位置和速度，小行星与飞行器的距离和相对速度由飞行器的标称轨道和小行星的位置速度得到，从而选取满足约束要求的星体，导航初始时刻，飞行器对日定向，且体坐标系与轨道坐标系重合，从而采用双矢量定姿法，期望姿态由飞行器位置和小行星位置得到，根据初始姿态和期望姿态，以及转动速度计算出姿态转换时间，选定时间最短的多颗小行星当作导航星，形成导航小行星列表；其特征在于，第二、处理导航小行星图像：形心确定方法分两步：

首先确定“亮心”，

$p_{\mathrm{cb}}=\frac{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{ip}}_{\mathrm{ij}}}{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{p}_{\mathrm{ij}}},$        $l_{\mathrm{cb}}=\frac{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\mathrm{jp}}_{\mathrm{ij}}}{\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{p}_{\mathrm{ij}}}$

式中，p_ij是对应像素的亮度值，i是像元，j是像线，p_cb和l_cb是计算的导航小行星图像亮心，m×n＝1024×1024是相机的分辨率；

然后根据亮心确定“形心”，

p₀＝p_cb-γR_ccosφ，l₀＝l_cb-γR_csinφ

式中， $\mathrm{\φ}=a\tan \left(\frac{A_{\mathrm{cy}}}{A_{\mathrm{cx}}}\right),$ A_c是相机坐标系太阳的方向，由指向太阳的惯性视线矢量A_I经过坐标转换得到；A_c＝T_CIA_I，其中T_CI为惯性坐标系到相机坐标系的转换矩阵，R_c是小行星半径对应的像素， $R_c=\frac{{\mathrm{RK}}_{x}f}{\mathrm{\ρ}},$ R是小行星的半径，ρ是探测器与小行星的距离，f是相机的焦距，K_x是从长度单位到像素的转换；γ是偏差因子，取值在0和1之间，表示目标半径误差ΔR引起的形心偏差，其具体表达式为

$\mathrm{\γ}=\frac{3\mathrm{\π\ΔR}}{16}\left[\frac{\sin \mathrm{\α}(1+\cos \mathrm{\α})}{(\mathrm{\π}-\mathrm{\α})\cos \mathrm{\α}+\sin \mathrm{\α}}\right]$

α是太阳相位角，

最后进行形心像素的计算，

惯性视线矢量

利用转换矩阵T_CI旋转到照相机坐标系，

${\stackrel{\→}{V}}_C=\left[\begin{array}{c}{V}_{C1}\\ V_{C2}\\ V_{C3}\end{array}\right]=T_{\mathrm{CI}}{\stackrel{\→}{V}}_I$

视线矢量

确定，转换到二维的照相机焦平面中，

$\left[\begin{array}{c}x\\ y\end{array}\right]=\frac{f}{V_{C3}}\left[\begin{array}{c}{V}_{C1}\\ V_{C2}\end{array}\right]$

其中，

f照相机焦距mm，

V_C1，V_C2，V_C3视线矢量在相机系内分量，

x，y视线矢量焦平面内投影mm，

接着，由光的畸变引起的x和y偏差Δx，Δy，

$\left[\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\end{array}\right]=Q{\left[\begin{array}{cccccc}{v}_1& v_2& v_3& v_4& v_5& v_6\end{array}\right]}^T$

其中，

$Q=\left[\begin{array}{cccccc}-\mathrm{yr}& {\mathrm{xr}}^2& -{\mathrm{yr}}^3& {\mathrm{xr}}^4& \mathrm{xy}& x^2\\ \mathrm{xr}& {\mathrm{yr}}^2& {\mathrm{xr}}^3& {\mathrm{yr}}^4& y^2& \mathrm{xy}\end{array}\right]$

r＝x²+y²，v₁ v₂ v₃ v₄ v₅ v₆是光畸变系数；修正的图像位置x′，y′表示为

$\left[\begin{array}{c}x\′\\ y\′\end{array}\right]=\left[\begin{array}{c}x+\mathrm{\Δx}\\ y+\mathrm{\Δy}\end{array}\right]$

最后，直角坐标系转化到像元和像线，

$\left[\begin{array}{c}p\\ l\end{array}\right]=\left[\begin{array}{ccc}{K}_x& K_{\mathrm{xy}}& K_{\mathrm{xxy}}\\ K_{\mathrm{yx}}& K_y& K_{\mathrm{yxy}}\end{array}\right]\left[\begin{array}{c}x\′\\ y\′\\ x\′y\′\end{array}\right]+\left[\begin{array}{c}{p}_0\\ l_0\end{array}\right]$

K阵是从毫米到像元/像线的转换矩阵，p₀，l₀是CCD的中心像元和像线；

第三、基于多颗小行星进行自主轨道确定：利用上述得到的多颗小行星的形心确定飞行器轨道；

1)定义参考轨道参数：

$X^{*}\left(t\right)={\left[\begin{array}{cccccc}x& y& z& \stackrel{.}{x}& \stackrel{.}{y}& \stackrel{.}{z}\end{array}\right]}^T,$ 其中，r＝[x y z]^T和 $\stackrel{\‾}{v}={\left[\begin{array}{ccc}\stackrel{.}{x}& \stackrel{.}{y}& \stackrel{.}{z}\end{array}\right]}^T$ 分别为探测器在日心黄道坐标系的位置和速度矢量；

更新的轨道参数：

X′(t)＝X^*(t)+ΔX(t)，其中，ΔX(t)为估计轨道修正量，由于参考轨道参数和真实轨道相差不大，那么在一段时间内轨道修正量是线性的，即在一时间点的轨道参数修正量ΔX(t₀)利用状态转换矩阵线性地映射到任何其它时间点t上，即ΔX(t)＝Φ(t)ΔX(t₀)

其中，给定时间t的状态转移矩阵Φ满足

$\stackrel{.}{\mathrm{\Φ}}=\frac{\∂\stackrel{.}{X}\left(t\right)}{\∂X\left(t\right)}\frac{\∂X\left(t\right)}{\∂X\left(t_0\right)}=\mathrm{A\Φ}\left(t\right)$

式中， $A=\frac{\∂\stackrel{.}{X}\left(t\right)}{\∂X\left(t\right)},$ Φ(t₀)＝I_6×6，通过这些方程的积分来获得Φ的时间函数，***的状态方程写为：

$\stackrel{.}{X}=\left[\begin{array}{c}\stackrel{\‾}{v}\\ -\frac{{\mathrm{\μ}}_s}{r^3}\stackrel{\‾}{r}+\underset{i=1}{\overset{n_p}{\mathrm{\Σ}}}{\mathrm{\μ}}_i[\frac{{\stackrel{\‾}{r}}_{\mathrm{ri}}}{r_{\mathrm{ri}}^3}-\frac{{\stackrel{\‾}{r}}_{\mathrm{pi}}}{r_{\mathrm{pi}}^3}]-\frac{\mathrm{AG}}{{\mathrm{mr}}^3}\stackrel{\‾}{r}+\frac{k}{m}\stackrel{\‾}{T}+\stackrel{\‾}{a}\end{array}\right]$

式中：r探测器在日心坐标系内位置矢量；v探测器在日心坐标系内速度矢量；r_pi第i个摄动行星体在日心坐标系内的位置矢量，r_ri探测器相对第i个摄动行星体的位置矢量，其中r_ri＝r_pi-r；μ_δ太阳引力常数GM；μ_i第i个摄动天体的引力常数；n_p摄动天体的数量；A探测器的有效面积；G日光流量常数；T推进***推力矢量；k推力比例系数，近似为1；m探测器质量；a探测器总的附加加速度；

2)确定观测矩阵：

对于某一观测历元，观测矩阵为

$H=\left[\begin{array}{cccccc}\∂p/\∂x_1& \∂p/\∂x_2& \∂p/\∂x_3& 0& 0& 0\\ \∂l/\∂x_1& \∂l/\∂x_2& \∂l/\∂x_3& 0& 0& 0\end{array}\right]$

式中，与速度分量有关的偏微分都是零；

3)统一观测历元：

利用状态转移矩阵把各个观测历元的观测矩阵变换到同一观测历元，有

$\tilde{H}=\mathrm{H\Φ}$

其中

就是在某时刻的观测偏微分矩阵；

4)实现轨道确定：

利用基于UD协方差分解的递推加权最小二乘算法来确定探测器的轨道，其中的测量残差Y是图像处理得到的与计算预测得到的小行星中心对应的像素差，这样就利用多颗小行星的观测数据，即小行星中心的像素，通过多步递推来修正某观测历元的轨道参数，首先给出一个协方差矩阵P₀，观测加权矩阵W，矢量Y为观测质心与标称轨迹预报计算质心的差值，初始点状态批量滤波方程求解矢量

和正式协方差矩阵P：

$\hat{q}={[P_0^{-1}+{\tilde{H}}^{T}W\tilde{H}]}^{-1}{\tilde{H}}^T\mathrm{WY}$

$P={[P_0^{-1}+{\tilde{H}}^{T}W\tilde{H}]}^{-1}$

其中，

$W=\left[\begin{array}{cc}1/{\mathrm{\σ}}_0^2& 0\\ 0& 1/{\mathrm{\σ}}_0^2\end{array}\right]$

${\mathrm{\σ}}_0=\frac{{\tan }^{-1}(R/\mathrm{\ρ})}{13\×{10}^{-6}}$

R小行星假定半径、ρ小行星范围。

Images