CN104729510A - 一种空间目标相对伴飞轨道确定方法 - Google Patents

Abstract

本发明涉及一种空间目标相对伴飞轨道确定方法，属于航天深空探测技术领域。该方法对当前时刻的轨道状态进行确定步骤如下：1)根据轨道动力学方程递推出任意时刻轨道状态预测值，并由此得到预测观测量，并得到实际观测值与预测值的差值；2)利用观测方程和1)中得到的差值解出线性方程组的解，即轨道状态估计偏差解集合；3)对2)中的解按数值递增排列，并剔除相邻两项之差大于测量误差限的解。4)根据3)得到的解集合个数，若为奇数，那么轨道状态估计偏差为3)中集合最中间的一项；若为偶数，则轨道状态估计偏差为集合最中间连续两项的平均值。5)当前时刻轨道参数估计值即为轨道状态预测值与轨道状态估计偏差之和。本发明能够降低复杂动力学模型及噪声干扰。

Description

一种空间目标相对伴飞轨道确定方法

技术领域

本发明涉及一种空间目标相对伴飞轨道确定方法，属于航天深空探测技术领域。

背景技术

随着深空探测任务的日益增多，航天器绕飞、伴飞观测目标天体已经成为未来深空科学探测的重要任务和课题。若航天器已捕获目标天体，但对其各种信息掌握较少，需采用伴飞探测，确定航天器相对目标天体的位置、速度以及姿态等参数信息。为了成功地完成科学考察任务，航天器实现伴飞探测轨道保持与控制需要很高的轨道确定精度。航天器确保较高的轨道确定精度，以便航天器在需要的位置实现准确地定点悬停或在一定区域内实现伴飞，完成后续探测任务。轨道确定精度对伴飞探测轨道控制有很大的影响，如果航天器与小天体相对信息的确定存在较大误差，会导致航天器无法进行精确地轨道控制，进而严重偏离预期轨道甚至无法继续进行伴飞探测。

由于目标天体周围引力场分布不均匀以及可能存在不稳定的自旋状态和其他未知干扰，航天器伴飞轨道确定方法必须具备较强鲁棒性。同时，在实际观测中，由于各种因素的影响，例如敏感器测量误差、时间计量精度误差等，观测数据中存在一定数量准确度和精确度不高的数据。由于存在轨道动力学建模和观测误差等干扰因素，必须寻找更可靠、更精确的定轨方法，因此航天器伴飞探测轨道确定方法是当前各国航天科研部门重点发展的研究方向之一。

在已发展的轨道确定方法中，最小二乘是最普遍的轨道确定方法之一，它不考虑力学模型误差，在观测数据中含有非线性影响因素的情况下，最小二乘法不再是最优线性估计，它会降低轨道确定的精度。此外很多学者采用EKF或UKF进行轨道确定，在先技术[1](参见LI Heng-nian,LI Ji-sheng,HUANG Yong-xuan.A Dynamic Model Compensation Orbit Determination Method for Maneuvering Satellite[J].JOURNAL OF  ASTRONAUTICS,2010,31(10),2269-2275.)将EKF运用到动力学补偿卫星的轨道确定算法中，实现了卫星在连续推力控制下精确定轨，但系统中观测噪声干扰仅仅假设为高斯白噪声，且对于地球卫星而言，没有复杂的动力学环境，因此如果用于周围引力场分布不均的目标天体动力学模型，且含有其他有色噪声，该轨道确定方法将降低定轨精度，影响伴飞轨道控制效果。

在先技术[2](参见Bhaskaran S,Riedel J E,Synnott S P.Autonomous target tracking of small bodies during flybys[J].Advances in the Astronautical Sciences,2004,119,2079-2096.)，分析了EKF或UKF用法的局限性。为保证精度要求，采用EKF或UKF进行轨道确定时一般假设非主要摄动项、过程噪声、量测噪声均为零均值的高斯白噪声。当过程噪声、量测噪声为有色噪声或存在其他未知概率密度分布的噪声时，EKF或UKF的滤波性能将降低，甚至最终发散。这种情况不适合航天器伴飞探测轨道确定，无法满足任务需求。

在先技术[3](参见Khutorovsky Z N,Samotokhin A S,Alfriend K T.Guaranteed Approach for Orbit Determination with Limited Error Measurements[J].Journal of Guidance,Control,and Dynamics,2013,36(4):1186-1193.)，采用了中心算法求取切比雪夫中心点进行近地卫星轨道确定，降低了有色噪声对定轨滤波的影响，但是如果观测数据中存在数据误差偏大的值，利用中心算法求解轨道参数就会出现大的偏差，同样无法满足伴飞探测高精度的定轨要求。

发明内容

本发明的目的是针对现有空间目标相对伴飞探测轨道确定方法无法保证在有其他噪声信息误差干扰、目标天体周围动力学特性复杂的情况下达到高精度定轨要求的问题，提供一种空间目标相对伴飞轨道确定方法，该方法是通过求取轨道参数的中值点来进行定轨的，能够降低复杂动力学模型及噪声干扰。

该轨道确定方法能够实现空间目标相对伴飞探测轨道确定，保证确定 精度并且能够为航天器定点悬停任务提供定轨保证，以便实现定点悬停控制。该选取方法既可以应用于航天器定点悬停任务的轨道确定，也可以应用于航天器稳定绕飞空间目标任务轨道确定。

本发明的目的是通过下述技术方案实现的。

本发明所述方法是在如下轨道动力学模型基础上进行轨道确定的。

航天器在空间目标固联坐标系下伴飞轨道动力学模型状态方程和观测方程如下：

$\left\{\begin{array}{c}\stackrel{.}{X}=F(X,U,t)+W\\ Y=G(X,t)+\mathrm{\δ}\end{array}\right.---\left(1\right)$

其中为航天器相对于空间目标的三轴位置与三轴速度；U为三轴主动控制加速度；F根据动力学特性定义的关于X,U的函数；W表示动力学建模中的干扰误差总和；Y为观测量向量；G是由X转化为Y的函数；δ为测量误差。

将式(1)在t_k时刻(航天器在轨运行的任意时刻)某一状态X^*处进行Taylor展开得到线性化模型

$\left\{\begin{array}{c}\mathrm{\Δ}\stackrel{.}{X}=A(X,U,t_k)\mathrm{\ΔX}+W\\ \mathrm{\ΔY}=C(X,t_k)\mathrm{\ΔX}+\mathrm{\δ}\end{array}\right.---\left(2\right)$

式(2)中ΔX＝X-X^*，A、C表达式为

$A(X,U,t_k)=\frac{\∂F(X,U,t_k)}{\∂k}{|}_{X_k},C(X,t_k)=\frac{\∂G(X,t_k)}{\∂X}{|}_{X_k}---\left(3\right)$

并且由A能够得到t_k时刻的状态转移矩阵Φ(t_k,t_k+1)，所以线性化离散后的模型为

$\left\{\begin{array}{c}\mathrm{\Δ}{X}_{k+1}=\mathrm{\Φ}(t_k,t_{k+1}){\mathrm{\ΔX}}_k+W\\ {\mathrm{\ΔY}}_k=C(X,t_k){\mathrm{\ΔX}}_k+\mathrm{\δ}\end{array}\right.---\left(4\right)$

从t₀时刻开始进行观测，记录连续n+1时刻的观测量并且预先估计t₀和t_n时刻(进行轨道确定的当前时刻)的初始轨道参数  $X_0={[x_0,y_0,z_0,{\stackrel{.}{x}}_0,{\stackrel{.}{y}}_0,{\stackrel{.}{z}}_0]}^T,X_n={[x_n,y_n,z_n,{\stackrel{.}{x}}_n,{\stackrel{.}{y}}_n,{\stackrel{.}{z}}_n]}^T,$ 可以根据轨道递推得到任意连续时刻状态转移矩阵Φ(t_k,t_k+1)。那么，从t_n时刻起，倒推t_k时刻ΔX_k为

${\mathrm{\ΔX}}_k=\underset{m=k}{\overset{n-1}{\mathrm{\Π}}}{\mathrm{\Φ}}^{-1}(t_k,t_{k+1})\·{\mathrm{\ΔX}}_n={\mathrm{\Φ}}_k\·{\mathrm{\ΔX}}_n---\left(5\right)$

本发明所述的一种空间目标相对伴飞轨道确定方法，对当前时刻t_n的轨道状态X_n进行确定步骤如下：

1)根据轨道递推求出任意时刻X_k，并由此得到预测观测量Y_k，将实际观测值与预测值做差得到ΔY_k。

2)利用式(4)的第二个方程解如下线性方程组

ΔY_k＝C(X,t_k)Φ_k·ΔX_n＝Η·ΔX_n    (6)

考虑观测量向量的维数，利用一个或连续几个时刻的观测量组成观测量集合与轨道状态向量维数相同，确保方程(6)有唯一解；任意k时刻都能得到一组线性方程组解ΔX_n，k，线性方程组的任意一项表示为

$h_{i1}\·{\mathrm{\Δx}}_n+h_{i2}\·{\mathrm{\Δy}}_n+h_{i3}\·{\mathrm{\Δz}}_n+h_{i4}\·\mathrm{\Δ}{\stackrel{.}{x}}_n+h_{i5}\·\mathrm{\Δ}{\stackrel{.}{y}}_n+h_{i6}\·\mathrm{\Δ}{\stackrel{.}{z}}_n={\mathrm{\Δy}}_i---\left(7\right)$

其中i(i≤6)表示测量向量的第i个测量量。

3)对求得的n组ΔX_n,k当中每一项Δx_n,i,k参数按数值大小递增排列，得到新的参数序列剔除相邻两个之差大于测量误差限的点。因为测量误差δ存在上下限，即|δ_i|≤ε_i，对于排在新序列前列的参数，如果存在参数之差满足剔除及之前的参数；同理在新序列后列的参数，若满足剔除 及之后的参数。最终每一项的参数序列为

$\{\mathrm{\Δ}{\tilde{x}}_{n,i,1},\mathrm{\Δ}{\tilde{x}}_{n,i,2},...,\mathrm{\Δ}{\tilde{x}}_{n,i,m}\}(m\≤n).$

4)由3)对2)中求解的参数进行递增排序和剔除超出误差限的参数后得到m(m≤n)组参数。若m为奇数，那么这组序列的中值点为

${\mathrm{\Δx}}_{n,i,c}=\mathrm{\Δ}{\tilde{x}}_{n,i,m+1/2}---\left(8\right)$

若m为偶数，那么序列的中值点为

${\mathrm{\Δx}}_{n,i,c}=0.5(\mathrm{\Δ}{\tilde{x}}_{n,i,m/2}+\mathrm{\Δ}{\tilde{x}}_{n,i,m/2+1})---\left(9\right)$

即有

${[{\mathrm{\Δx}}_{n,c},{\mathrm{\Δy}}_{n,c},{\mathrm{\Δz}}_{n,c},\mathrm{\Δ}{\stackrel{.}{x}}_{n,c},\mathrm{\Δ}{\stackrel{.}{y}}_{n,c},\mathrm{\Δ}{\stackrel{.}{z}}_{n,c]}}^T={[{\mathrm{\Δx}}_{n-1,c},{\mathrm{\Δx}}_{n,2,c},{\mathrm{\Δx}}_{n,3,c},{\mathrm{\Δx}}_{n,4,c},{\mathrm{\Δ}}_{n,5,c},{\mathrm{\Δx}}_{n,6,c}]}^T---\left(10\right)$

这样就得到了t_n时刻的 ${\mathrm{\ΔX}}_n={[{\mathrm{\Δx}}_{n,c},{\mathrm{\Δy}}_{n,c},{\mathrm{\Δz}}_{n,c},\mathrm{\Δ}{\stackrel{.}{x}}_{n,c},\mathrm{\Δ}{\stackrel{.}{y}}_{n,c},\mathrm{\Δ}{\stackrel{.}{z}}_{n,c}]}^T.$

5)当前t_n时刻轨道参数估计值：

${\hat{X}}_n=X_n+{\mathrm{\ΔX}}_n---\left(11\right)$

根据测量误差干扰和动力学模型复杂程度，选取n的值可以调节轨道确定精度。

有益效果

本发明所给出的空间目标相对伴飞探测轨道确定方法，改善了非线性、含其他噪声干扰对系统的影响，提高了轨道确定精度，能够实现多种伴飞与绕飞探测任务轨道确定精度需求。同时，由于在定点悬停任务中，航天器位置与定点位置偏差较小，可以直接选择在定点位置处对伴飞轨道动力学模型进行线性离散化，简化算法、大大降低计算量，更加适用悬停任务的轨道确定，为准确实现定点悬停控制提供保证。对以小行星Eros433为目标天体进行伴飞探测定点悬停数学仿真结果显示，添加有色噪声干扰，在悬停任务初始位置误差超过一公里的情况下，采用该方法可以最终将轨道确定精度控制在100m以内，相比较在相同条件下加权最小二乘误差超过500m，EKF仿真后期轨道确定误差为收敛；并且悬停控制加速度最终趋于稳定，即实现了定点悬停探测任务。

具体实施方式

下面结合实施例对本发明做详细说明。

以探测器定对小天体定点悬停探测为例，详细说明本发明的具体实现方法。

首先给出探测器实现定点悬停的动力学模型，该轨道确定方法是基于此模型实现的。目标小天体固联坐标系下伴飞探测轨道动力学方程

$\left\{\begin{array}{c}\stackrel{..}{x}=2\mathrm{\ω}\stackrel{.}{y}+{\mathrm{\ω}}^{2}x+u_x+\frac{\∂V}{\∂x}+n_x\\ \stackrel{..}{y}=-2\mathrm{\ω}\stackrel{.}{x}+{\mathrm{\ω}}^{2}y+u_y+\frac{\∂V}{\∂y}+n_y\\ \stackrel{..}{z}=u_z+\frac{\∂V}{\∂z}+n_z\end{array}\right.---\left(12\right)$

上式中(x y z)^T是航天器在小天体固联坐标系下位置，ω为小天体自旋 角速度。U＝(u_x u_y u_z)^T为航天器在小天体固联坐标系中的控制加速度；N＝(n_x n_y n_z)^T为干扰加速度；V为小天体非球形引力位函数， 是小天体非球型引力摄动项。

若需要实现航天器的定点悬停控制，即要使航天器相对于目标小天体保持在某一固定位置，相对速度保持为零。与之前一致，将动力学方程改写成状态方程，假设悬停点位置为X₀＝[x₀,y₀,z₀,0,0,0]^T，任意时刻航天器位置悬停点与实际位置之间存在偏差设为ΔX，任意时刻将动力学方程在X₀处Taylor展开，得到线性化方程

$\mathrm{\Δ}\stackrel{..}{X}+\mathrm{P\Δ}\stackrel{.}{X}+\mathrm{Q\ΔX}-\frac{{\∂}^{2}V}{{\∂X}^2}{|}_{X_0}\mathrm{\ΔX}=-Q{X}_0+\frac{\∂V}{\∂X}{|}_{X_0}+U+N---\left(13\right)$

上式中

$P=\left(\begin{array}{ccc}0& -2\mathrm{\ω}& 0\\ 2\mathrm{\ω}& 0& 0\\ 0& 0& 0\end{array}\right),Q=-\left(\begin{array}{ccc}{\mathrm{\ω}}^2& 0& 0\\ 0& {\mathrm{\ω}}^2& 0\\ 0& 0& 0\end{array}\right)$

对式(13)采用滑模变结构控制实现定点悬停，设滑模为其中系数c为正常数。因此经过推倒得到控制加速度U为

$U={\mathrm{QX}}_0-\∂V/\∂X{|}_{X_0}+(P-c)\mathrm{\Δ}\stackrel{.}{X}+(Q-{\∂}^{2}V/{\∂X}^2{|}_{X_0})\mathrm{\ΔX}-\mathrm{kS}+\mathrm{\λsgn}\left(S\right)---\left(14\right)$

上式中系数k为正常数，λ＜-d，d是干扰加速度误差上限。sgn(S)形式如下：

$\mathrm{sgn}\left(s_i\right)=\frac{s_i}{\left|s_i\right|+\mathrm{\δ}},i=x,y,z;$ δ为极小量       (15) 

将控制加速度代U入方程式(11)，状态方程简化为

$\mathrm{\Δ}\stackrel{.}{X}={\left[\begin{array}{cccccc}0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\\ -\mathrm{kc}& 0& 0& -k-c& 0& 0\\ 0& -\mathrm{kc}& 0& 0& -k-c& 0\\ 0& 0& -\mathrm{kc}& 0& 0& -k-c\end{array}\right]}_{X_0}\left(\begin{array}{c}\mathrm{\Δx}\\ \mathrm{\Δy}\\ \mathrm{\Δz}\\ \mathrm{\Δ}\stackrel{.}{x}\\ \mathrm{\Δ}\stackrel{.}{y}\\ \mathrm{\Δ}\stackrel{.}{z}\end{array}\right)+\left(\begin{array}{c}0\\ 0\\ 0\\ \mathrm{\λsgn}(\mathrm{\Δ}\stackrel{.}{x}+\mathrm{c\Δx})\\ \mathrm{\λsgn}(\mathrm{\Δ}\stackrel{.}{y}+\mathrm{c\Δy})\\ \mathrm{\λsgn}(\mathrm{\Δ}\stackrel{.}{z}+\mathrm{c\Δz})\end{array}\right)+\left(\begin{array}{c}0\\ 0\\ 0\\ n_x\\ n_y\\ n_z\end{array}\right)---\left(16\right)$

上式中令

$A\left(X_0\right)={\left[\begin{array}{cccccc}0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0\\ -\mathrm{kc}& 0& 0& -k-c& 0& 0\\ 0& -\mathrm{kc}& 0& 0& -k-c& 0\\ 0& 0& -\mathrm{kc}& 0& 0& -k-c\end{array}\right]}_{X_0},W=\left(\begin{array}{c}0\\ 0\\ 0\\ \mathrm{\λsgn}(\mathrm{\Δ}\stackrel{.}{x}+\mathrm{c\Δx})\\ \mathrm{\λsgn}(\mathrm{\Δ}\stackrel{.}{y}+\mathrm{c\Δy})\\ \mathrm{\λsgn}(\mathrm{\Δ}\stackrel{.}{z}+\mathrm{c\Δz})\end{array}\right)+\left(\begin{array}{c}0\\ 0\\ 0\\ n_x\\ n_y\\ n_z\end{array}\right)$

这样状态方程形式如式(2)所示。

根据小天体与航天器的相对距离ρ、方位角α和俯仰角β，得到观测方程为

$\left\{\begin{array}{c}\mathrm{\ρ}=\sqrt{x^2+y^2+z^2}\\ \mathrm{\α}=\arctan (x/y)\\ \mathrm{\β}=\arcsin (z/\mathrm{\ρ})\end{array}\right.---\left(17\right)$

其中Y＝[ρ,α,β]^T是测量值，实际观测方程中各测量量均带有测量误差干扰δ，且|δ_i|≤ε_i,i＝ρ,α,β。按照式(1)的形式将式(17)改写为向量形式

Y＝G(X,t)+δ     (18) 

同样对观测方程在X₀处线性化，与轨道动力学方程一起离散后得到

$\left\{\begin{array}{c}{\mathrm{\ΔX}}_{k+1}=\mathrm{\Φ}\left(t_k\right){\mathrm{\ΔX}}_k+W\\ {\mathrm{\ΔY}}_k=C\left(X_0\right){\mathrm{\ΔX}}_k+\mathrm{\δ}\end{array}\right.---\left(19\right)$

上式中由A(X₀)能够得到t_k时刻的状态转移矩阵Φ(t_k)，C(X₀)为式(3)中所表示在X₀处的取值的。

在t₀时刻开始进行观测，记录连续n+1时刻的观测量并且已知t₀和t_n时刻的预估轨道参数X₀＝[x₀,y₀,z₀,0,0,0]^T，X_n＝[x₀,y₀,z₀,0,0,0]^T，任意时刻状态对t₀时刻转移矩阵为Φ(t_k)。那么，从t_n时刻起，倒推t_k时刻ΔX_k为

ΔX_k＝Φ(t_k)·(Φ^-1(t_n)ΔX_n)＝Φ_k·ΔX_n     (20) 

本发明所述的相对伴飞轨道确定方法，对探测器定点悬停任务当前时刻t_n的轨道状态X_n进行确定步骤如下：

1)由任意时刻X_k-1,X_k，其中X_k-1＝X₀,X_k＝X₀得到预测观测量Y_k-1,Y_k，将实际观测值与预测值做差得到ΔY_k-1,ΔY_k。

2)利用式(4)第二个方程以及ΔY_k-1和ΔY_k解如下线性方程

$\left[\begin{array}{c}{\mathrm{\ΔY}}_{k-1}\\ {\mathrm{\ΔY}}_k\end{array}\right]=\left[\begin{array}{c}C\left(X_0\right){\mathrm{\Φ}}_{k-1}\\ C\left(X_0\right){\mathrm{\Φ}}_k\end{array}\right]\·{\mathrm{\ΔX}}_n=H\·{\mathrm{\ΔX}}_n---\left(21\right)$

式(21)中ΔY_k-1,ΔY_k都是三维向量，ΔX_n为六维向量。在任意k时刻能够得到一组线性方程组解ΔX_n，k。式(21)展开得

$\left\{\begin{array}{c}{h}_{11}\·{\mathrm{\Δx}}_n+h_{12}\·{\mathrm{\Δy}}_n+h_{13}\·{\mathrm{\Δz}}_n+h_{14}\·\mathrm{\Δ}{\stackrel{.}{x}}_n+h_{15}\·\mathrm{\Δ}{\stackrel{.}{y}}_n+h_{16}\·\mathrm{\Δ}{\stackrel{.}{z}}_n={\mathrm{\Δ\ρ}}_{k-1}\\ h_{21}\·{\mathrm{\Δx}}_n+h_{22}\·{\mathrm{\Δy}}_n+h_{23}\·{\mathrm{\Δz}}_n+h_{24}\·\mathrm{\Δ}{\stackrel{.}{x}}_n+h_{25}\·\mathrm{\Δ}{\stackrel{.}{y}}_n+h_{26}\·\mathrm{\Δ}{\stackrel{.}{z}}_n={\mathrm{\Δ\α}}_{k-1}\\ h_{31}\·{\mathrm{\Δx}}_n+h_{32}\·{\mathrm{\Δy}}_n+h_{33}\·{\mathrm{\Δz}}_n+h_{34}\·\mathrm{\Δ}{\stackrel{.}{x}}_n+h_{35}\·\mathrm{\Δ}{\stackrel{.}{y}}_n+h_{36}\·\mathrm{\Δ}{\stackrel{.}{z}}_n={\mathrm{\Δ\β}}_{k-1}\\ h_{41}\·{\mathrm{\Δx}}_n+h_{42}\·{\mathrm{\Δy}}_n+h_{43}\·{\mathrm{\Δz}}_n+h_{44}\·\mathrm{\Δ}{\stackrel{.}{x}}_n+h_{45}\·\mathrm{\Δ}{\stackrel{.}{y}}_n+h_{46}\·\mathrm{\Δ}{\stackrel{.}{z}}_n={\mathrm{\Δ\ρ}}_k\\ h_{51}\·{\mathrm{\Δx}}_n+h_{52}\·{\mathrm{\Δy}}_n+h_{53}\·{\mathrm{\Δz}}_n+h_{54}\·\mathrm{\Δ}{\stackrel{.}{x}}_n+h_{55}\·\mathrm{\Δ}{\stackrel{.}{y}}_n+h_{56}\·\mathrm{\Δ}{\stackrel{.}{z}}_n={\mathrm{\Δ\α}}_k\\ h_{61}\·{\mathrm{\Δx}}_n+h_{62}\·{\mathrm{\Δy}}_n+h_{63}\·{\mathrm{\Δz}}_n+h_{64}\·\mathrm{\Δ}{\stackrel{.}{x}}_n+h_{65}\·\mathrm{\Δ}{\stackrel{.}{y}}_n+h_{66}\·\mathrm{\Δ}{\stackrel{.}{z}}_n={\mathrm{\Δ\β}}_k\end{array}\right.---\left(22\right)$

其中h_ij是矩阵H的第i行第j列的数值。

3)对求得的n组ΔX_n,k当中每一项Δx_n,i,k参数按数值大小递增排列，得到新的参数序列剔除相邻两个之差大于测量误差限的点。对于排在新序列前列的参数，如果存在参数之差满足 剔除及之前的参数；同理在新序列后列的参数，若满足剔除及之后的参数。最终每一项的参数序列为(m≤n)。

4)由3)对2)中求解的参数进行递增排序和剔除超出误差限的参数后得到m(m≤n)组参数。若m为奇数，那么这组序列的中值点为

${\mathrm{\Δx}}_{n,i,c}=\mathrm{\Δ}{\tilde{x}}_{n,i,m+1/2}---\left(23\right)$

若m为偶数，那么序列的中值点为

${\mathrm{\Δx}}_{n,i,c}=0.5(\mathrm{\Δ}{\tilde{x}}_{n,i,m/2}+\mathrm{\Δ}{\tilde{x}}_{n,i,m/2+1})---\left(24\right)$

即有

${[{\mathrm{\Δx}}_{n,c},{\mathrm{\Δy}}_{n,c},{\mathrm{\Δz}}_{n,c},\mathrm{\Δ}{\stackrel{.}{x}}_{n,c},\mathrm{\Δ}{\stackrel{.}{y}}_{n,c},\mathrm{\Δ}{\stackrel{.}{z}}_{n,c}]}^T={[{\mathrm{\Δx}}_{n,1,c},{\mathrm{\Δx}}_{n,2,c},{\mathrm{\Δx}}_{n,3,c},{\mathrm{\Δx}}_{n,4,c},{\mathrm{\Δx}}_{n,5,c},{\mathrm{\Δx}}_{n,6,c}]}^T---\left(25\right)$

这样就得到了t_n时刻的 ${\mathrm{\ΔX}}_n={[{\mathrm{\Δx}}_{n,c},{\mathrm{\Δy}}_{n,c},{\mathrm{\Δz}}_{n,c},\mathrm{\Δ}{\stackrel{.}{x}}_{n,c},\mathrm{\Δ}{\stackrel{.}{y}}_{n,c},\mathrm{\Δ}{\stackrel{.}{z}}_{n,c}]}^T.$

4)当前t_n时刻轨道参数估计值：

${\hat{X}}_n=X_n+{\mathrm{\ΔX}}_n={[x_0+{\mathrm{\Δx}}_{n,c},y_0+{\mathrm{\Δy}}_{n,c},z_0+{\mathrm{\Δz}}_{n,c},\mathrm{\Δ}{\stackrel{.}{x}}_{n,c},\mathrm{\Δ}{\stackrel{.}{y}}_{n,c},\mathrm{\Δ}{\stackrel{.}{z}}_{n,c}]}^T---\left(26\right)$

以上五步骤即从某一时刻开始进行观测，到当前时刻为止，利用这段 时间的观测数据对当前时刻轨道参数进行估计，得到当前时刻需要的轨道状态。

有以上案例分析易知，在进行悬停任务中，动力学模型的状态方程和观测方程线性离散化后的参数矩阵均为定常值，并未随轨道位置改变而改变，这样大大减少了星载计算机计算量，有利于轨道确定与控制实时性。

Claims (1)

1.一种空间目标相对伴飞轨道确定方法，是在如下轨道动力学模型基础上进行轨道确定的；航天器在空间目标固联坐标系下伴飞轨道动力学模型状态方程和观测方程如下：

$\left\{\begin{array}{c}\stackrel{\·}{X}=F(X,U,t)+W\\ Y=G(X,t)+\mathrm{\δ}\end{array}\right.---\left(1\right)$

其中为航天器相对于空间目标的三轴位置与三轴速度；U为三轴主动控制加速度；F根据动力学特性定义的关于X,U的函数；W表示动力学建模中的干扰误差总和；Y为测量值向量；G是由X转化为Y的函数；δ为测量误差；

将式(1)在t_k时刻(航天器在轨运行的任意时刻)某一状态X^*处进行Taylor展开得到线性化模型

$\left\{\begin{array}{c}\mathrm{\Δ}\stackrel{\·}{X}=A(X,U,t_k)\mathrm{\ΔX}+W\\ \mathrm{\ΔY}=C(X,t_k)\mathrm{\ΔX}+\mathrm{\δ}\end{array}\right.---\left(2\right)$

式(2)中ΔX＝X-X^*，A、C表达式为

$A(X,U,t_k)=\frac{\∂F(X,U,t_k)}{\∂X}{|}_{X_k},C(X,t_k)=\frac{\∂G(X,t_k)}{\∂X}{|}_{X_k}---\left(3\right)$

并且由A能够得到t_k时刻的状态转移矩阵Φ(t_k,t_k+1)，所以线性化离散后的模型为

$\left\{\begin{array}{c}\mathrm{\Δ}{X}_{k+1}=\mathrm{\Φ}(t_k,t_{k+1})\mathrm{\Δ}{X}_k+W\\ \mathrm{\Δ}{Y}_k=C(X,t_k)\mathrm{\Δ}{X}_k+\mathrm{\δ}\end{array}\right.---\left(4\right)$

从t₀时刻开始进行观测，记录连续n+1时刻的观测量并且预先估计t₀和t_n时刻(进行轨道确定的当前时刻)的初始轨道参数 $X_0={[x_0,y_0,z_0,{\stackrel{\·}{x}}_0,{\stackrel{\·}{y}}_0,{\stackrel{\·}{z}}_0]}^T,$ $X_n={[x_n,y_n,z_n,{\stackrel{\·}{x}}_n,{\stackrel{\·}{y}}_n,{\stackrel{\·}{z}}_n]}^T,$ 可以根据轨道递推得到任意连续时刻状态转移矩阵Φ(t_k,t_k+1)；那么，从t_n时刻起，倒推t_k时刻ΔX_k为

$\mathrm{\Δ}{X}_k=\underset{m=k}{\overset{n-1}{\mathrm{\Π}}}{\mathrm{\Φ}}^{-1}(t_k,t_{k+1})\·\mathrm{\Δ}{X}_n={\mathrm{\Φ}}_k\·\mathrm{\Δ}{X}_n---\left(5\right),$

其特征在于，对当前时刻t_n的轨道状态X_n进行确定步骤如下：

1)根据轨道递推求出任意时刻X_k，并由此得到预测观测量Y_k，将实际观测值与预测值做差得到ΔY_k；

2)利用式(4)的第二个方程解如下线性方程组

ΔY_k＝C(X,t_k)Φ_k·ΔX_n＝Η·ΔX_n       (6)

考虑观测量向量的维数，利用一个或连续几个时刻的观测量组成观测量集合与轨道状态向量维数相同，确保方程(6)有唯一解；在任意k时刻都能得到一组线性方程组解ΔX_n,k，线性方程组的任意一项表示为

$h_{i1}\·\mathrm{\Δ}{x}_n+h_{i2}\·\mathrm{\Δ}{y}_n+h_{i3}\·\mathrm{\Δ}{z}_n+h_{i4}\·\mathrm{\Δ}{\stackrel{\·}{x}}_n+h_{i5}\·\mathrm{\Δ}{\stackrel{\·}{y}}_n+h_{i6}\·\mathrm{\Δ}{\stackrel{\·}{z}}_n=\mathrm{\Δ}{y}_i---\left(7\right)$

其中i(i≤6)表示测量向量的第i个测量量；

3)对求得的n组ΔX_n,k当中每一项Δx_n,i,k参数按数值大小递增排列，得到新的参数序列剔除相邻两个之差大于测量误差限的点；因为测量误差δ存在上下限，即|δ_i|≤ε_i，对于排在新序列前列的参数，如果存在参数之差满足剔除及之前的参数；同理在新序列后列的参数，若满足剔除及之后的参数；最终每一项的参数序列为 $\{\mathrm{\Δ}{\tilde{x}}_{n,i,1},\mathrm{\Δ}{\tilde{x}}_{n,i,2},...,\mathrm{\Δ}{\tilde{x}}_{n,i,m}\}(m\≤n);$

4)由3)对2)中求解的参数进行递增排序和剔除超出误差限的参数后得到m(m≤n)组参数；若m为奇数，那么这组序列的中值点为

$\mathrm{\Δ}{x}_{n,i,c}=\mathrm{\Δ}{\tilde{x}}_{n,i,m+1/2}---\left(8\right)$

若m为偶数，那么序列的中值点为

$\mathrm{\Δ}{x}_{n,i,c}=0.5(\mathrm{\Δ}{\tilde{x}}_{n,i,m/2}+\mathrm{\Δ}{\tilde{x}}_{n,i,m/2+1})---\left(9\right)$

即有

${[\mathrm{\Δ}{x}_{n,c},\mathrm{\Δ}{y}_{n,c},\mathrm{\Δ}{z}_{n,c},\mathrm{\Δ}{\stackrel{\·}{x}}_{n,c},\mathrm{\Δ}{\stackrel{\·}{y}}_{n,c},\mathrm{\Δ}{\stackrel{\·}{z}}_{n,c}]}^T={[\mathrm{\Δ}{x}_{n,1,c},\mathrm{\Δ}{x}_{n,2,c},\mathrm{\Δ}{x}_{n,3,c},\mathrm{\Δ}{x}_{n,4,c},\mathrm{\Δ}{x}_{n,5,c},\mathrm{\Δ}{x}_{n,6,c}]}^T---\left(10\right)$

这样就得到了t_n时刻的 ${\mathrm{\Δ}{X}_n=[\mathrm{\Δ}{x}_{n,c},\mathrm{\Δ}{y}_{n,c},\mathrm{\Δ}{z}_{n,c},\mathrm{\Δ}{\stackrel{\·}{x}}_{n,c},\mathrm{\Δ}{\stackrel{\·}{y}}_{n,c},\mathrm{\Δ}{\stackrel{\·}{z}}_{n,c}]}^T;$

5)当前t_n时刻轨道参数估计值：

${\hat{X}}_n=X_n+\mathrm{\Δ}{X}_n---\left(11\right)$

根据测量误差干扰和动力学模型复杂程度，选取n的值调节轨道确定精度。