一种基于无迹递推的空间操作相对轨道控制误差计算方法

摘要

一种基于无迹递推的空间操作相对轨道控制误差计算方法，涉及分析空间操作任务中相对轨道控制误差的分布状况领域；针对初始位置、敏感器和执行机构等设计参数，分析空间操作任务中相对轨道控制误差的分布状况；该方法根据系统的维数进行确定性的采样，再通过无迹变换在时间域内对主要随机变量的均值和误差协方差矩阵进行递推式的处理。这一方法在数据处理量和处理时间方面优于基于蒙特卡洛计算的误差计算方法。同时，由于无迹变换只近似系统状态变量的高斯分布，而不是近似系统的非线性函数，该方法避免了求解解析的雅克比矩阵，较线性协方差计算方法更容易实现。

说明书

技术领域

本发明涉及一种分析空间操作任务中相对轨道控制误差的分布状况领域，特别是一种基于无迹递推的空间操作相对轨道控制误差计算方法。

背景技术

相对轨道控制技术是空间操作任务中一项不可或缺的核心技术。该项技术主要是使航天器安全、有效地抵近空间任务目标，并在目标附近建立并维持稳定的相对轨道运动状态，为开展后继在轨操作提供有力的支撑。相对轨道控制的误差是衡量该技术的一项核心指标。

相对轨道控制的误差主要是通过误差的均值和协方差矩阵进行描述。由于相对轨道控制系统是一个构成复杂的非线性随机系统，系统的动态特性显著。少数几次的数值仿真结果难于反映其控制误差的整体分布状况。通常情况下，该控制系统需要进行成百上千次的蒙特卡洛仿真计算，再根据仿真结果统计控制系统的误差均值和误差协方差矩阵。但是，这种基于蒙特卡洛计算的误差计算方法数据处理量大、等待周期长。线性协方差计算计算方法是近些年来提出的另一种误差计算计算方法(Linear Covariance Techniques for Orbital Rendezvous Analysis and Autonomous Onboard Mission Planning,Journal of Guidance Control and Dynamics,Vol.29,No.6,2006；Multiple Event Triggers in Linear Covariance Analysis for Spacecraft Rendezvous,Journal of Guidance Control and Dynamics,Vol.35,No.2,2012)。该方法在小量假设前提下，通过线性函数近似非线性函数，简化了系统模型，然后在时间域内递推求解系统中主要状态变量的均值和协方差矩阵。该方法较基于蒙特卡洛计算的误差计算计算方法，在数据的处理量和处理时间均大为减少。但是，该方法在仿真分析之前需要对非线性函数求解解析的雅克比矩阵，而相对轨道控制系统中复杂的函数映射关系将使得这一求解工作更为困难。

发明内容

本发明的目的在于克服现有技术的上述不足，提供一种基于无迹递推的空间操作相对轨道控制误差计算方法，在数据处理量和处理时间方面优于基于蒙特卡洛计算的误差计算方法，避免了求解解析的雅克比矩阵，较线性协方差计算方法更容易实现。

本发明的上述目的是通过如下技术方案予以实现的：

一种基于无迹递推的空间操作相对轨道控制误差计算方法，误差计算方法包括如下步骤：

步骤(一)、误差递推计算

设定t_k时刻误差递推状态的均值为误差协方差矩阵为可以表示为：

$X_{100}^{t_k}=f_E\left(r_T^{t_k},v_T^{t_k},r_C^{t_k},v_C^{t_k},{\mathrm{\Delta r}}_{TC,NaviErr}^{t_k},{\mathrm{\Delta v}}_{TC,NaviErr}^{t_k}\right)$

${\mathrm{PX}}_{100}^{t_k}=f_{Cov}(r_T^{t_k},v_T^{t_k},r_C^{t_k},v_C^{t_k},{\mathrm{\Delta r}}_{TC,NaviErr}^{t_k},{\mathrm{\Delta v}}_{TC,NaviErr}^{t_k})$

其中，

为t_k时刻目标航天器在地心赤道惯性坐标系下的位置；

为t_k时刻目标航天器在地心赤道惯性坐标系下的速度；

为t_k时刻追踪航天器在地心赤道惯性坐标系下的位置；

为t_k时刻追踪航天器在地心赤道惯性坐标系下的速度；

为t_k时刻目标航天器在追踪航天器轨道坐标系下的相对位置的导航误差；

为t_k时刻目标航天器在追踪航天器轨道坐标系下的相对速度的导航误差；

f_E()为均值计算函数，按照如下关系进行计算：

$f_E(a_1,a_2,...,a_N)=\left(\begin{array}{c}E\left(a_1\right)\\ E\left(a_2\right)\\ \begin{array}{c}.\\ .\\ .\end{array}\\ E\left(a_N\right)\end{array}\right)$

其中，

a₁、a₂、…、a_N为函数f_E()输入的自变量；

E()为函数f_E()输入自变量a₁、a₂、…、a_N的均值；

N为输入变量的个数；N为正整数；

f_Cov()为误差协方差计算函数，按照如下关系进行计算：

其中，

b₁、b₂、…、b_M分别为函数f_Cov()输入的自变量；

P()为函数f_Cov()输入自变量b₁、b₂、…、b_N的误差协方差矩阵；

P(b_i,b_j)中(i,j∈[1、2、……M]且i≠j)是输入自变量b_i和b_j之间的误差协方差矩阵；

M为输入的自变量的个数；M为正整数；

步骤(二)、相对位置速度计算

根据t_k时刻误差递推状态的均值和误差协方差计算相对位置速度输入变量的均值和误差协方差采用UT变换计算t_k时刻追踪航天器相对位置和相对速度在目标航天器轨道坐标系下的均值t_k时刻追踪航天器相对位置和相对速度在目标航天器轨道坐标系下的误差协方差矩阵

并根据误差递推状态的递推时间t_k以及进行轨道机动的时刻t_Cmd进行判断，当t_k≠t_Cmd时，不进行轨道机动，进入步骤(三)，当t_k＝t_Cmd时，进行轨道机动，进入步骤(四)；

步骤(三)、轨道积分递推计算

当t_k时刻不进行轨道机动时，根据t_k时刻误差递推状态的均值和误差协方差计算t_k时刻轨道积分递推输入变量的均值和误差协方差采用UT变换计算t_k+1时刻轨道积分递推输出变量的均值和误差协方差并将和传递给t_k+1时刻误差递推状态的均值和误差协方差矩阵

步骤(四)、轨道机动递推计算

根据t_k时刻误差递推状态的均值和误差协方差计算t_k时刻轨道机动递推输入变量的均值和误差协方差采用UT变换计算t_k时刻轨道机动递推输出变量的均值和误差协方差并将和传递给t_k时刻误差递推状态的均值和误差协方差矩阵

步骤(五)、按照步骤(一)至步骤(四)，依次递推完成t₁时刻到t₂时刻、t₂时刻到t₃时刻、…、t_N-1时刻到t_N时刻的递推计算，并根据步骤(二)中，计算得到的t₁时刻、t₂时刻、…、t_N时刻追踪航天器相对位置和相对速度在目标航天器轨道坐标系下的均值和误差协方差矩阵通过追踪航天器相对位置和相对速度在目标航天器轨道坐标系下的均值和误差协方差矩阵可以分析t₁时刻、t₂时刻、…、t_N时刻追踪航天器的相对位置和相对速度在目标航天器轨道坐标系下，评价空间操作相对轨道控制的控制精度。

在上述的一种基于无迹递推的空间操作相对轨道控制误差计算方法，所述步骤(二)中，相对位置速度输入变量的均值和误差协方差的计算方法为：

$X_{200}^{t_k}=f_E\left(r_T^{t_k},v_T^{t_k},r_C^{t_k},v_C^{t_k}\right)$

${\mathrm{PX}}_{200}^{t_k}=f_{Cov}\left(r_T^{t_k},v_T^{t_k},r_C^{t_k},v_C^{t_k}\right)$

其中，

为t_k时刻目标航天器在地心赤道惯性坐标系下的位置；

为t_k时刻目标航天器在地心赤道惯性坐标系下的速度；

为t_k时刻追踪航天器在地心赤道惯性坐标系下的位置；

为t_k时刻追踪航天器在地心赤道惯性坐标系下的速度；

f_E()为均值计算函数；

f_Cov()为误差协方差计算函数。

在上述的一种基于无迹递推的空间操作相对轨道控制误差计算方法，所述步骤(二)中，采用UT变换计算t_k时刻追踪航天器在目标航天器轨道坐标系下的相对位置和相对速度的均值和误差协方差矩阵计算方法为：

对t_k时刻输入的追踪航天器相对位置和相对速度在目标航天器轨道坐标系下的均值误差协方差矩阵以及t_k时刻非线性函数映射关系进行计算；得到：

$Y_{200}^{t_k}=f_E({\mathrm{\Delta r}}_{CT}^{t_k},{\mathrm{\Delta v}}_{CT}^{t_k})$

${\mathrm{PY}}_{200}^{t_k}=f_{Cov}({\mathrm{\Delta r}}_{CT}^{t_k},{\mathrm{\Delta v}}_{CT}^{t_k})$

其中，

f_E()为均值计算函数；

f_Cov()为误差协方差计算函数。

在上述的一种基于无迹递推的空间操作相对轨道控制误差计算方法，所述步骤(三)中，t_k时刻轨道积分递推输入变量的均值和误差协方差的计算方法为：

$\left(\begin{array}{c}{X}_{301}^{t_k}=\left(\begin{array}{c}{X}_{100}^{t_k}\\ f_E\left(a_{T,U\mathrm{mod}el},a_{C,U\mathrm{mod}el},{\xi }_{TC,Mea},r_{C,NaviErr},v_{C,NaviErr}\right)\end{array}\right)\\ =f_E\left(r_T^{t_k},v_T^{t_k},r_C^{t_k},v_C^{t_k},{\mathrm{\Delta r}}_{TC,NaviErr}^{t_k},{\mathrm{\Delta v}}_{TC,NaviErr}^{t_k},a_{T,U\mathrm{mod}el},a_{C,U\mathrm{mod}el},{\xi }_{TC,Mea},r_{C,NaviErr},v_{C,NaviErr}\right)\end{array}\right)$

$\left(\begin{array}{c}{\mathrm{PX}}_{301}^{t_k}=\left(\begin{array}{cccccc}{\mathrm{PX}}_{100}^{t_k}& 0_{18\times 3}& 0_{18\times 3}& 0_{18\times 3}& 0_{18\times 3}& 0_{18\times 3}\\ 0_{3\times 18}& P\left(a_{T,U\mathrm{mod}el}\right)& 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}\\ 0_{3\times 18}& 0_{3\times 3}& P\left(a_{C,U\mathrm{mod}el}\right)& 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}\\ 0_{3\times 18}& 0_{3\times 3}& 0_{3\times 3}& P\left({\xi }_{TC,Mea}\right)& 0_{3\times 3}& 0_{3\times 3}\\ 0_{3\times 18}& 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}& P\left(r_{C,NaviErr}\right)& 0_{3\times 3}\\ 0_{3\times 18}& 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}& P\left(v_{C,NaviErr}\right)\end{array}\right)\\ =f_{Cov}\left(r_T^{t_k},v_T^{t_k},r_C^{t_k},v_C^{t_k},{\mathrm{\Delta r}}_{TC,NaviErr}^{t_k},{\mathrm{\Delta v}}_{TC,NaviErr}^{t_k},a_{T,U\mathrm{mod}el},a_{C,U\mathrm{mod}el},{\xi }_{TC,Mea},r_{C,NaviErr},v_{C,NaviErr}\right)\end{array}\right)$

其中，

为t_k时刻误差递推状态的均值；

为t_k时刻误差递推状态的误差协方差；

a_T,Unmodel为追踪航天器的加速度误差；

a_C,Unmodel为目标航天器的加速度误差；

ξ_TC,Mea为相对测量误差；

r_C,NaviErr为追踪航天器位置的导航误差；

v_C,NaviErr为追踪航天器速度的导航误差；

为t_k时刻目标航天器相对于追踪航天器的位置的导航误差；

为t_k时刻目标航天器相对于追踪航天器的速度的导航误差。

在上述的一种基于无迹递推的空间操作相对轨道控制误差计算方法，所述步骤(三)中，采用UT变换计算t_k+1时刻轨道积分递推输出变量的均值和误差协方差并将和传递给t_k+1时刻误差递推状态的均值和误差协方差矩阵具体方法为：

根据t_k时刻输入的状态变量的均值和误差协方差和t_k时刻非线性函数映射关系

进行计算，得到：

$X_{100}^{t_{k+1}}=Y_{301}^{t_{k+1}}=f_E\left(r_T^{t_{k+1}},v_T^{t_{k+1}},r_C^{t_{k+1}},v_C^{t_{k+1}},{\mathrm{\Delta r}}_{TC,NaviErr}^{t_{k+1}},{\mathrm{\Delta v}}_{TC,NaviErr}^{t_{k+1}}\right)$

${\mathrm{PX}}_{100}^{t_{k+1}}={\mathrm{PY}}_{301}^{t_{k+1}}=f_{Cov}\left(r_T^{t_{k+1}},v_T^{t_{k+1}},r_C^{t_{k+1}},v_C^{t_{k+1}},{\mathrm{\Delta r}}_{TC,NaviErr}^{t_{k+1}},{\mathrm{\Delta v}}_{TC,NaviErr}^{t_{k+1}}\right)$

其中，为t_k+1时刻目标航天器在地心赤道惯性坐标系下的位置；

为t_k+1时刻目标航天器在地心赤道惯性坐标系下的速度；

为t_k+1时刻追踪航天器在地心赤道惯性坐标系下的位置；

为t_k+1时刻追踪航天器在地心赤道惯性坐标系下的速度；

为t_k+1时刻目标航天器相对于追踪航天器轨道坐标系下的相对位置的导航误差；

为t_k+1时刻目标航天器相对于追踪航天器轨道坐标系下的相对速度的导航误差。

在上述的一种基于无迹递推的空间操作相对轨道控制误差计算方法，所述步骤(四)中，t_k时刻轨道机动递推输入变量的均值和误差协方差的计算方法为：

$\left(\begin{array}{c}{X}_{302}^{t_k}=\left(\begin{array}{c}{X}_{100}^{t_k}\\ f_E\left(r_{C,NaviErr},v_{C,NaviErr},{\mathrm{\Delta v}}_{C,CmdErr}\right)\end{array}\right)\\ =f_E\left(r_T^{t_k},v_T^{t_k},r_C^{t_k},v_C^{t_k},{\mathrm{\Delta r}}_{TC,NaviErr}^{t_k},{\mathrm{\Delta v}}_{TC,NaviErr}^{t_k}{r}_{C,NaviErr},v_{C,NaviErr},{\mathrm{\Delta v}}_{C,CmdErr}\right)\end{array}\right)$

$\left(\begin{array}{c}{\mathrm{PX}}_{301}^{t_k}=\left(\begin{array}{cccc}{\mathrm{PX}}_{100}^{t_k}& 0_{18\times 3}& 0_{18\times 3}& 0_{18\times 3}\\ 0_{3\times 18}& P\left(r_{C,NaviErr}\right)& 0_{3\times 3}& 0_{3\times 3}\\ 0_{3\times 18}& 0_{3\times 3}& P\left(v_{C,NaviErr}\right)& 0_{3\times 3}\\ 0_{3\times 18}& 0_{3\times 3}& 0_{3\times 3}& P\left({\mathrm{\Delta v}}_{C,CmdErr}\right)\end{array}\right)\\ =f_{Cov}\left(r_T^{t_k},v_T^{t_k},r_C^{t_k},v_C^{t_k},{\mathrm{\Delta r}}_{TC,NaviErr}^{t_k},{\mathrm{\Delta v}}_{TC,NaviErr}^{t_k}{r}_{C,NaviErr},v_{C,NaviErr},{\mathrm{\Delta v}}_{C,CmdErr}\right)\end{array}\right)$

其中，

为t_k时刻误差递推状态的均值；

为t_k时刻误差递推状态的误差协方差；

r_C,NaviErr为追踪航天器位置的导航误差；

v_C,NaviErr追踪航天器速度的导航误差；

Δv_C,CmdErr为追踪航天器轨道机动的执行误差；

为t_k时刻目标航天器相对于追踪航天器的位置的导航误差；

为t_k时刻目标航天器相对于追踪航天器的速度的导航误差。

在上述的一种基于无迹递推的空间操作相对轨道控制误差计算方法，所述步骤(四)中，采用UT变换计算t_k时刻轨道机动递推输出变量的均值和误差协方差具体计算方法为：

根据和传递给t_k时刻误差递推状态的均值和误差协方差矩阵t_k时刻输入的状态变量的均值和误差协方差和t_k时刻非线性函数映射关系

进行计算，得到：

$X_{100}^{t_k}=Y_{302}^{t_k}=f_E(r_T^{t_k},v_T^{t_k},r_C^{t_k},v_C^{t_k},{\mathrm{\Delta r}}_{TC,NaviErr}^{t_k},{\mathrm{\Delta v}}_{TC,NaviErr}^{t_k})$

${\mathrm{PX}}_{100}^{t_k}={\mathrm{PY}}_{302}^{t_k}=f_{Cov}\left(r_T^{t_k},v_T^{t_k},r_C^{t_k},v_C^{t_k},{\mathrm{\Delta r}}_{TC,NaviErr}^{t_k},{\mathrm{\Delta v}}_{TC,NaviErr}^{t_k}\right)$

其中，

为用于替换的t_k时刻目标航天器在地心赤道惯性坐标系下的位置；

为用于替换的t_k时刻目标航天器在地心赤道惯性坐标系下的速度；

为用于替换的t_k时刻追踪航天器在地心赤道惯性坐标系下的位置；

为用于替换的t_k时刻追踪航天器在地心赤道惯性坐标系下的速度；

为用于替换的t_k时刻目标航天器在追踪航天器轨道坐标系下的相对位置的导航误差；

为用于替换的t_k时刻目标航天器在追踪航天器轨道坐标系下的相对速度的导航误差。

本发明与现有技术相比具有如下优点：

(1)本发明将滤波估计中通常使用的无迹变换方法，应用于航天器相对轨道控制误差分布状况的计算分析之中，根据系统的维数进行确定性的采样，再通过无迹变换在时间域内对主要随机变量的均值和协方差进行递推式地处理；

(2)本发明所述方法根据系统维数进行确定性的采样计算，相比于基于蒙特卡洛仿真计算的误差计算方法，减少了数据处理量，缩短了计算时间；

(3)本发明所述方法只近似系统状态变量的高斯分布而不是近似系统的非线性函数，相比于线性协方差计算方法，无需求解解析的雅克比矩阵，更容易实现计算。

附图说明

图1为基于无迹递推的空间操作相对轨道控制误差计算方法的计算流程示意图；

图2为采用基于无迹递推的空间操作相对轨道控制误差计算方法得到的相对位置误差分布算例结果；

图3为基于蒙特卡洛计算的误差计算方法得到的相对位置误差分布算例结果。

具体实施方式

下面结合附图和具体实施例对本发明作进一步详细的描述：

如图1所示为基于无迹递推的空间操作相对轨道控制误差计算方法的计算流程示意图，由图可知，一种基于无迹递推的空间操作相对轨道控制误差计算方法，其特征在于：误差计算方法包括如下步骤：

步骤(一)、误差递推计算

设定t_k时刻误差递推状态的均值为误差协方差矩阵为可以表示为：

$X_{100}^{t_k}=f_E(r_T^{t_k},v_T^{t_k},r_C^{t_k},v_C^{t_k},{\mathrm{\Delta r}}_{TC,NaviErr}^{t_k},{\mathrm{\Delta v}}_{TC,NaviErr}^{t_k})$

${\mathrm{PX}}_{100}^{t_k}=f_{Cov}(r_T^{t_k},v_T^{t_k},r_C^{t_k},v_C^{t_k},{\mathrm{\Delta r}}_{TC,NaviErr}^{t_k},{\mathrm{\Delta v}}_{TC,NaviErr}^{t_k})$

其中，

为t_k时刻目标航天器在地心赤道惯性坐标系下的位置，m；

为t_k时刻目标航天器在地心赤道惯性坐标系下的速度，m/s；

为t_k时刻追踪航天器在地心赤道惯性坐标系下的位置，m；

为t_k时刻追踪航天器在地心赤道惯性坐标系下的速度，m/s；

为t_k时刻目标航天器在追踪航天器轨道坐标系下的相对位置的导航误差，m；

为t_k时刻目标航天器在追踪航天器轨道坐标系下的相对速度的导航误差，m/s；

f_E()为均值计算函数，按照如下关系进行计算：

$f_E(a_1,a_2,...,a_N)=\left(\begin{array}{c}E\left(a_1\right)\\ E\left(a_2\right)\\ \begin{array}{c}.\\ .\\ .\end{array}\\ E\left(a_N\right)\end{array}\right)$

其中，

a₁、a₂、…、a_N为函数f_E()输入的自变量；

E()为函数f_E()输入自变量a₁、a₂、…、a_N的均值；

N为输入变量的个数；N为正整数；

f_Cov()为误差协方差计算函数，按照如下关系进行计算：

其中，

b₁、b₂、…、b_M分别为函数f_Cov()输入的自变量；

P()为函数f_Cov()输入自变量b₁、b₂、…、b_N的误差协方差矩阵；

P(b_i,b_j)中(i,j∈[1、2、……M]且i≠j)是输入自变量b_i和b_j之间的误差协方差矩阵；

M为输入的自变量的个数；M为正整数；

步骤(二)、相对位置速度计算

根据t_k时刻误差递推状态的均值和误差协方差计算相对位置速度输入变量的均值和误差协方差采用UT变换计算t_k时刻追踪航天器相对位置和相对速度在目标航天器轨道坐标系下的均值t_k时刻追踪航天器相对位置和相对速度在目标航天器轨道坐标系下的误差协方差矩阵

并根据误差递推状态的递推时间t_k以及进行轨道机动的时刻t_Cmd进行判断，当t_k≠t_Cmd时，不进行轨道机动，进入步骤(三)，当t_k＝t_Cmd时，进行轨道机动，进入步骤(四)；

步骤(三)、轨道积分递推计算

当t_k时刻不进行轨道机动时，根据t_k时刻误差递推状态的均值和误差协方差计算t_k时刻轨道积分递推输入变量的均值和误差协方差采用UT变换计算t_k+1时刻轨道积分递推输出变量的均值和误差协方差并将和传递给t_k+1时刻误差递推状态的均值和误差协方差矩阵

步骤(四)、轨道机动递推计算

根据t_k时刻误差递推状态的均值和误差协方差计算t_k时刻轨道机动递推输入变量的均值和误差协方差采用UT变换计算t_k时刻轨道机动递推输出变量的均值和误差协方差并将和传递给t_k时刻误差递推状态的均值和误差协方差矩阵

步骤(五)、按照步骤(一)至步骤(四)，依次递推完成t₁时刻到t₂时刻、t₂时刻到t₃时刻、…、t_N-1时刻到t_N时刻的递推计算，并根据步骤(二)中，计算得到的t₁时刻、t₂时刻、…、t_N时刻追踪航天器相对位置和相对速度在目标航天器轨道坐标系下的均值和误差协方差矩阵通过追踪航天器相对位置和相对速度在目标航天器轨道坐标系下的均值和误差协方差矩阵可以分析t₁时刻、t₂时刻、…、t_N时刻追踪航天器的相对位置和相对速度在目标航天器轨道坐标系下，评价空间操作相对轨道控制的控制精度。

所述步骤(二)中，相对位置速度输入变量的均值和误差协方差的计算方法为：

$X_{200}^{t_k}=f_E\left(r_T^{t_k},v_T^{t_k},r_C^{t_k},v_C^{t_k}\right)$

${\mathrm{PX}}_{200}^{t_k}=f_{Cov}\left(r_T^{t_k},v_T^{t_k},r_C^{t_k},v_C^{t_k}\right)$

其中，

为t_k时刻目标航天器在地心赤道惯性坐标系下的位置，m；

为t_k时刻目标航天器在地心赤道惯性坐标系下的速度，m/s；

为t_k时刻追踪航天器在地心赤道惯性坐标系下的位置，m；

为t_k时刻追踪航天器在地心赤道惯性坐标系下的速度，m/s；

f_E()为均值计算函数；

f_Cov()为误差协方差计算函数。

采用UT变换计算t_k时刻追踪航天器在目标航天器轨道坐标系下的相对位置和相对速度的均值和误差协方差矩阵计算方法为：

对t_k时刻输入的追踪航天器相对位置和相对速度在目标航天器轨道坐标系下的均值误差协方差矩阵以及t_k时刻非线性函数映射关系进行计算；得到：

$Y_{200}^{t_k}=f_E({\mathrm{\Delta r}}_{CT}^{t_k},{\mathrm{\Delta v}}_{CT}^{t_k})$

${\mathrm{PY}}_{200}^{t_k}=f_{Cov}({\mathrm{\Delta r}}_{CT}^{t_k},{\mathrm{\Delta v}}_{CT}^{t_k})$

其中，

f_E()为均值计算函数；

f_Cov()为误差协方差计算函数。

所述步骤(三)中，t_k时刻轨道积分递推输入变量的均值和误差协方差的计算方法为：

$\left(\begin{array}{c}{X}_{301}^{t_k}=\left(\begin{array}{c}{X}_{100}^{t_k}\\ f_E\left(a_{T,U\mathrm{mod}el},a_{C,U\mathrm{mod}el},{\xi }_{TC,Mea},r_{C,NaviErr},v_{C,NaviErr}\right)\end{array}\right)\\ =f_E\left(r_T^{t_k},v_T^{t_k},r_C^{t_k},v_C^{t_k},{\mathrm{\Delta r}}_{TC,NaviErr}^{t_k},{\mathrm{\Delta v}}_{TC,NaviErr}^{t_k},a_{T,U\mathrm{mod}el},a_{C,U\mathrm{mod}el},{\xi }_{TC,Mea},r_{C,NaviErr},v_{C,NaviErr}\right)\end{array}\right)$

$\left(\begin{array}{c}{\mathrm{PX}}_{301}^{t_k}=\left(\begin{array}{cccccc}{\mathrm{PX}}_{100}^{t_k}& 0_{18\times 3}& 0_{18\times 3}& 0_{18\times 3}& 0_{18\times 3}& 0_{18\times 3}\\ 0_{3\times 18}& P\left(a_{T,U\mathrm{mod}el}\right)& 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}\\ 0_{3\times 18}& 0_{3\times 3}& P\left(a_{C,U\mathrm{mod}el}\right)& 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}\\ 0_{3\times 18}& 0_{3\times 3}& 0_{3\times 3}& P\left({\xi }_{TC,Mea}\right)& 0_{3\times 3}& 0_{3\times 3}\\ 0_{3\times 18}& 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}& P\left(r_{C,NaviErr}\right)& 0_{3\times 3}\\ 0_{3\times 18}& 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}& 0_{3\times 3}& P\left(v_{C,NaviErr}\right)\end{array}\right)\\ =f_{Cov}\left(r_T^{t_k},v_T^{t_k},r_C^{t_k},v_C^{t_k},{\mathrm{\Delta r}}_{TC,NaviErr}^{t_k},{\mathrm{\Delta v}}_{TC,NaviErr}^{t_k},a_{T,U\mathrm{mod}el},a_{C,U\mathrm{mod}el},{\xi }_{TC,Mea},r_{C,NaviErr},v_{C,NaviErr}\right)\end{array}\right)$

其中，

为t_k时刻误差递推状态的均值；

为t_k时刻误差递推状态的误差协方差；

a_T,Unmodel为追踪航天器的加速度误差；

a_C,Unmodel为目标航天器的加速度误差；

ξ_TC,Mea为相对测量误差；

r_C,NaviErr为追踪航天器位置的导航误差；

v_C,NaviErr为追踪航天器速度的导航误差；

为t_k时刻目标航天器相对于追踪航天器的位置的导航误差；

为t_k时刻目标航天器相对于追踪航天器的速度的导航误差。

采用UT变换计算t_k+1时刻轨道积分递推输出变量的均值和误差协方差并将和传递给t_k+1时刻误差递推状态的均值和误差协方差矩阵具体方法为：

根据t_k时刻输入的状态变量的均值和误差协方差和t_k时刻非线性函数映射关系

进行计算，得到：

$X_{100}^{t_{k+1}}=Y_{301}^{t_{k+1}}=f_E\left(r_T^{t_{k+1}},v_T^{t_{k+1}},r_C^{t_{k+1}},v_C^{t_{k+1}},{\mathrm{\Delta r}}_{TC,NaviErr}^{t_{k+1}},{\mathrm{\Delta v}}_{TC,NaviErr}^{t_{k+1}}\right)$

${\mathrm{PX}}_{100}^{t_{k+1}}={\mathrm{PY}}_{301}^{t_{k+1}}=f_{Cov}\left(r_T^{t_{k+1}},v_T^{t_{k+1}},r_C^{t_{k+1}},v_C^{t_{k+1}},{\mathrm{\Delta r}}_{TC,NaviErr}^{t_{k+1}},{\mathrm{\Delta v}}_{TC,NaviErr}^{t_{k+1}}\right)$

其中，为t_k+1时刻目标航天器在地心赤道惯性坐标系下的位置，m；

为t_k+1时刻目标航天器在地心赤道惯性坐标系下的速度，m/s；

为t_k+1时刻追踪航天器在地心赤道惯性坐标系下的位置，m；

为t_k+1时刻追踪航天器在地心赤道惯性坐标系下的速度，m/s；

为t_k+1时刻目标航天器相对于追踪航天器轨道坐标系下的相对位置的导航误差，m；

为t_k+1时刻目标航天器相对于追踪航天器轨道坐标系下的相对速度的导航误差，m/s。

所述步骤(四)中，t_k时刻轨道机动递推输入变量的均值和误差协方差的计算方法为：

$\left(\begin{array}{c}{X}_{302}^{t_k}=\left(\begin{array}{c}{X}_{100}^{t_k}\\ f_E\left(r_{C,NaviErr},v_{C,NaviErr},{\mathrm{\Delta v}}_{C,CmdErr}\right)\end{array}\right)\\ =f_E\left(r_T^{t_k},v_T^{t_k},r_C^{t_k},v_C^{t_k},{\mathrm{\Delta r}}_{TC,NaviErr}^{t_k},{\mathrm{\Delta v}}_{TC,NaviErr}^{t_k}{r}_{C,NaviErr},v_{C,NaviErr},{\mathrm{\Delta v}}_{C,CmdErr}\right)\end{array}\right)$

$\left(\begin{array}{c}{\mathrm{PX}}_{301}^{t_k}=\left(\begin{array}{cccc}{\mathrm{PX}}_{100}^{t_k}& 0_{18\times 3}& 0_{18\times 3}& 0_{18\times 3}\\ 0_{3\times 18}& P\left(r_{C,NaviErr}\right)& 0_{3\times 3}& 0_{3\times 3}\\ 0_{3\times 18}& 0_{3\times 3}& P\left(v_{C,NaviErr}\right)& 0_{3\times 3}\\ 0_{3\times 18}& 0_{3\times 3}& 0_{3\times 3}& P\left({\mathrm{\Delta v}}_{C,CmdErr}\right)\end{array}\right)\\ =f_{Cov}\left(r_T^{t_k},v_T^{t_k},r_C^{t_k},v_C^{t_k},{\mathrm{\Delta r}}_{TC,NaviErr}^{t_k},{\mathrm{\Delta v}}_{TC,NaviErr}^{t_k}{r}_{C,NaviErr},v_{C,NaviErr},{\mathrm{\Delta v}}_{C,CmdErr}\right)\end{array}\right)$

其中，

为t_k时刻误差递推状态的均值；

为t_k时刻误差递推状态的误差协方差；

r_C,NaviErr为追踪航天器位置的导航误差，m；

v_C,NaviErr追踪航天器速度的导航误差，m/s；

Δv_C,CmdErr为追踪航天器轨道机动的执行误差；

为t_k时刻目标航天器相对于追踪航天器的位置的导航误差，m；

为t_k时刻目标航天器相对于追踪航天器的速度的导航误差，m/s。

采用UT变换计算t_k时刻轨道机动递推输出变量的均值和误差协方差具体计算方法为：

根据和传递给t_k时刻误差递推状态的均值和误差协方差矩阵t_k时刻输入的状态变量的均值和误差协方差和t_k时刻非线性函数映射关系

进行计算，得到：

$X_{100}^{t_k}=Y_{302}^{t_k}=f_E(r_T^{t_k},v_T^{t_k},r_C^{t_k},v_C^{t_k},{\mathrm{\Delta r}}_{TC,NaviErr}^{t_k},{\mathrm{\Delta v}}_{TC,NaviErr}^{t_k})$

${\mathrm{PX}}_{100}^{t_k}={\mathrm{PY}}_{302}^{t_k}=f_{Cov}\left(r_T^{t_k},v_T^{t_k},r_C^{t_k},v_C^{t_k},{\mathrm{\Delta r}}_{TC,NaviErr}^{t_k},{\mathrm{\Delta v}}_{TC,NaviErr}^{t_k}\right)$

其中，

为用于替换的t_k时刻目标航天器在地心赤道惯性坐标系下的位置，m；

为用于替换的t_k时刻目标航天器在地心赤道惯性坐标系下的速度，m/s；

为用于替换的t_k时刻追踪航天器在地心赤道惯性坐标系下的位置，m；

为用于替换的t_k时刻追踪航天器在地心赤道惯性坐标系下的速度，m/s；

为用于替换的t_k时刻目标航天器在追踪航天器轨道坐标系下的相对位置的导航误差，m；

为用于替换的t_k时刻目标航天器在追踪航天器轨道坐标系下的相对速度的导航误差m/s。

所述步骤(二)中相对位置速度计算中，t_k时刻非线性的函数映射关系可以表示为

$\left(\begin{array}{c}{\mathrm{\Delta r}}_{CT}^{t_k}\\ {\mathrm{\Delta v}}_{CT}^{t_k}\end{array}\right)=\left(\begin{array}{cc}{L}_{OI}\left(r_T^{t_k},v_T^{t_k}\right)& 0_{3\times 3}\\ -L_{OI}\left(r_T^{t_k},v_T^{t_k}\right){\omega }^{\times }\left(r_T^{t_k},v_T^{t_k}\right)& L_{OI}\left(r_T^{t_k},v_T^{t_k}\right)\end{array}\right)\left(\begin{array}{c}{r}_C^{t_k}-r_T^{t_k}\\ v_C^{t_k}-v_T^{t_k}\end{array}\right)$

其中，为t_k时刻由地心赤道惯性坐标系到目标航天器的轨道坐标系的转移矩阵，具体计算公式为

$L_{OI}\left(r_T^{t_k},v_T^{t_k}\right)={\left(\begin{array}{ccc}\frac{\left(r_T^{t_k}\times v_T^{t_k}\right)\times r_T^{t_k}}{|\left(r_T^{t_k}\times v_T^{t_k}\right)\times r_T^{t_k}|}& -\frac{r_T^{t_k}\times v_T^{t_k}}{|r_T^{t_k}\times v_T^{t_k}|}& -\frac{r_T^{t_k}}{\left|r_T^{t_k}\right|}\end{array}\right)}^T$

为t_k时刻地心赤道惯性坐标系下目标航天器的轨道角速度，具体计算公式为

$\omega \left(r_T^{t_k},v_T^{t_k}\right)=\frac{r_T^{t_k}\times v_T^{t_k}}{\left|r_T^{t_k}\right|}$

为的叉乘矩阵。

所述步骤(三)轨道积分递推计算中，t_k时刻非线性的函数映射关系是指按照如下顺序关系(a轨道积分计算；b生成测量信息；c滤波估计计算；d相对导航误差计算)进行计算：

a.轨道积分计算

以t_k时刻目标航天器和追踪航天器在地心赤道惯性坐标系下的位置和速度作为初值，以t_k+1-t_k为积分步长，采用4阶RungeKutta方法(数学手册，高等教育出版社，1979.5，P687)对下述公式

$\frac{d^2{r}_T^{t_k}}{{\mathrm{dt}}^2}=-\frac{{\mathrm{\mu r}}_T^{t_k}}{|r_T^{t_k}{|}^3}+a_{J234}\left(r_T^{t_k}\right)+a_{T,Un\mathrm{mod}el}$

$\frac{d^2{r}_C^{t_k}}{{\mathrm{dt}}^2}=-\frac{{\mathrm{\mu r}}_C^{t_k}}{|r_C^{t_k}{|}^3}+a_{J234}\left(r_C^{t_k}\right)+a_{C,Un\mathrm{mod}el}$

进行积分，得到t_k+1时刻目标航天器和追踪航天器在地心赤道惯性坐标系下的其中，a_T,Unmodel和a_C,Unmodel分别为追踪航天器和目标航天器未建模的加速度误差，满足均值为E(a_T,Unmodel)和E(a_C,Unmodel)、协方差为P(a_T,Unmodel)和P(a_C,Unmodel)的白噪声分布；分别为的模值，分别为地球J2、J3、J4摄动项，分别按照如下函数计算计算

$\left|r_T^{t_k}\right|=\sqrt{{\left(x_T^{t_k}\right)}^2+{\left(y_T^{t_k}\right)}^2+{\left(z_T^{t_k}\right)}^2}$

$\left|r_C^{t_k}\right|=\sqrt{{\left(x_C^{t_k}\right)}^2+{\left(y_C^{t_k}\right)}^2+{\left(z_C^{t_k}\right)}^2}$

$\left(\begin{array}{c}{a}_{J234}\left(r_T^{t_k}\right)=\frac{{\mu }_E{J}_2{R}_E^2}{|r_T^{t_k}{|}^5}\left(\begin{array}{c}{x}_T^{t_k}\left(7.5{\left(z_T^{t_k}/\left|r_T^{t_k}\right|\right)}^2-1.5\right)\\ y_T^{t_k}\left(7.5{\left(z_T^{t_k}/\left|r_T^{t_k}\right|\right)}^2-1.5\right)\\ z_T^{t_k}\left(7.5{\left(z_T^{t_k}/\left|r_T^{t_k}\right|\right)}^2-4.5\right)\end{array}\right)\\ +\frac{{\mu }_E{J}_3{R}_E^3}{|r_T^{t_k}{|}^6}\left(\begin{array}{c}{x}_T^{t_k}\left(17.5{\left(z_T^{t_k}/\left|r_T^{t_k}\right|\right)}^3-7.5z_T^{t_k}/\left|r_T^{t_k}\right|\right)\\ y_T^{t_k}\left(17.5{\left(z_T^{t_k}/\left|r_T^{t_k}\right|\right)}^3-7.5z_T^{t_k}/\left|r_T^{t_k}\right|\right)\\ \left|r_T^{t_k}\right|\left(17.5{\left(z_T^{t_k}/\left|r_T^{t_k}\right|\right)}^4-15{\left(y_T^{t_k}/\left|r_T^{t_k}\right|\right)}^2+1.5\right)\end{array}\right)\\ +\frac{{\mu }_E{J}_4{R}_E^4}{|r_T^{t_k}{|}^7}\left(\begin{array}{c}{x}_T^{t_k}\left(39.375{\left(z_T^{t_k}/\left|r_T^{t_k}\right|\right)}^4-26.25{\left(z_T^{t_k}/\left|r_T^{t_k}\right|\right)}^2+1.875\right)\\ y_T^{t_k}\left(39.375{\left(z_T^{t_k}/\left|r_T^{t_k}\right|\right)}^4-26.25{\left(z_T^{t_k}/\left|r_T^{t_k}\right|\right)}^2+1.875\right)\\ \left|r_T^{t_k}\right|\left(39.375{\left(z_T^{t_k}/\left|r_T^{t_k}\right|\right)}^5-43.75{\left(z_T^{t_k}/\left|r_T^{t_k}\right|\right)}^3+9.375\left(z_T^{t_k}/\left|r_T^{t_k}\right|\right)\right)\end{array}\right)\end{array}\right)$

$\left(\begin{array}{c}{a}_{J234}\left(r_C^{t_k}\right)=\frac{{\mu }_E{J}_2{R}_E^2}{|r_C^{t_k}{|}^5}\left(\begin{array}{c}{x}_C^{t_k}\left(7.5{\left(z_C^{t_k}/\left|r_C^{t_k}\right|\right)}^2-1.5\right)\\ y_C^{t_k}\left(7.5{\left(z_C^{t_k}/\left|r_C^{t_k}\right|\right)}^2-1.5\right)\\ z_C^{t_k}\left(7.5{\left(z_C^{t_k}/\left|r_C^{t_k}\right|\right)}^2-4.5\right)\end{array}\right)\\ +\frac{{\mu }_E{J}_3{R}_E^3}{|r_C^{t_k}{|}^6}\left(\begin{array}{c}{x}_C^{t_k}\left(17.5{\left(z_C^{t_k}/\left|r_C^{t_k}\right|\right)}^3-7.5z_C^{t_k}/\left|r_C^{t_k}\right|\right)\\ y_C^{t_k}\left(17.5{\left(z_C^{t_k}/\left|r_C^{t_k}\right|\right)}^3-7.5z_C^{t_k}/\left|r_C^{t_k}\right|\right)\\ \left|r_C^{t_k}\right|\left(17.5{\left(z_C^{t_k}/\left|r_C^{t_k}\right|\right)}^4-15{\left(y_C^{t_k}/\left|r_C^{t_k}\right|\right)}^2+1.5\right)\end{array}\right)\\ +\frac{{\mu }_E{J}_4{R}_E^4}{|r_C^{t_k}{|}^7}\left(\begin{array}{c}{x}_C^{t_k}\left(39.375{\left(z_C^{t_k}/\left|r_C^{t_k}\right|\right)}^4-26.25{\left(z_C^{t_k}/\left|r_C^{t_k}\right|\right)}^2+1.875\right)\\ y_C^{t_k}\left(39.375{\left(z_C^{t_k}/\left|r_C^{t_k}\right|\right)}^4-26.25{\left(z_C^{t_k}/\left|r_C^{t_k}\right|\right)}^2+1.875\right)\\ \left|r_C^{t_k}\right|\left(39.375{\left(z_C^{t_k}/\left|r_C^{t_k}\right|\right)}^5-43.75{\left(z_C^{t_k}/\left|r_C^{t_k}\right|\right)}^3+9.375\left(z_C^{t_k}/\left|r_C^{t_k}\right|\right)\right)\end{array}\right)\end{array}\right)$

式中分别为矢量在地心赤道惯性坐标系的分量，分别为矢量在地心赤道惯性坐标系的分量。

b.生成测量信息

根据积分后得到的t_k+1时刻追踪航天器和目标航天器在地心赤道惯性坐标系下的位置按照如下公式计算t_k+1时刻目标航天器在追踪航天器轨道坐标系下的相对位置

${\mathrm{\Delta r}}_{TC}^{t_{k+1}}=L_{OI}(r_C^{t_{k+1}},v_C^{t_{k+1}})(r_T^{t_{k+1}}-r_C^{t_{k+1}})$

其中，为t_k+1时刻由地心赤道惯性坐标系到目标航天器的轨道坐标系的转移矩阵，具体计算公式为

$L_{OI}\left(r_C^{t_{k+1}},v_C^{t_{k+1}}\right)={\left(\begin{array}{ccc}\frac{\left(r_C^{t_{k+1}}\times v_C^{t_{k+1}}\right)\times r_C^{t_{k+1}}}{|\left(r_C^{t_{k+1}}\times v_C^{t_{k+1}}\right)\times r_C^{t_{k+1}}|}& -\frac{r_C^{t_{k+1}}\times v_C^{t_{k+1}}}{|r_C^{t_{k+1}}\times v_C^{t_{k+1}}|}& -\frac{r_C^{t_{k+1}}}{\left|r_C^{t_{k+1}}\right|}\end{array}\right)}^T$

根据t_k+1时刻目标航天器在追踪航天器轨道坐标系下的相对位置计算t_k+1时刻追踪航天器上的相对导航敏感器测量目标航天器的测量信息包括相对位置视线仰角和方位角

$Y_{TC,Mea}^{t_{k+1}}=\left(\begin{array}{c}{\rho }_{TC,Mea}^{t_{k+1}}\\ {\alpha }_{TC,Mea}^{t_{k+1}}\\ {\beta }_{TC,Mea}^{t_{k+1}}\end{array}\right)=\left(\begin{array}{c}\sqrt{{\left({\mathrm{\Delta x}}_{TC}^{t_{k+1}}\right)}^2+{\left({\mathrm{\Delta y}}_{TC}^{t_{k+1}}\right)}^2+{\left({\mathrm{\Delta z}}_{TC}^{t_{k+1}}\right)}^2}\\ \arctan 2\left(-{\mathrm{\Delta z}}_{TC}^{t_{k+1}},{\mathrm{\Delta x}}_{TC}^{t_{k+1}}\right)\\ \arctan \left(\frac{{\mathrm{\Delta y}}_{TC}^{t_{k+1}}}{\sqrt{{\left({\mathrm{\Delta x}}_{TC}^{t_{k+1}}\right)}^2+{\left({\mathrm{\Delta z}}_{TC}^{t_{k+1}}\right)}^2}}\right)\end{array}\right)+{\xi }_{TC,Mea}$

式中分别为矢量在轨道坐标系下的分量；ξ_TC,Mea为相对导航敏感器的测量误差，满足均值为E(ξ_TC,Mea)、协方差为P(ξ_TC,Mea)的白噪声分布。

根据t_k时刻追踪航天器在地心赤道惯性坐标系下的位置和速度按照如下关系计算追踪航天器位置和速度的导航信息

$r_{C,Navi}^{t_k}=r_C^{t_k}+r_{C,NaviErr}$

$v_{C,Navi}^{t_k}=v_C^{t_k}+v_{C,NaviErr}$

其中，r_C,NaviErr、v_C,NaviErr分别为追踪航天器位置速度的导航误差，满足均值为E(r_C,NaviErr)和E(v_C,NaviErr)、协方差为P(r_C,NaviErr)和P(v_C,NaviErr)的白噪声分布。

c.滤波估计计算

根据t_k时刻目标航天器和追踪航天器在地心赤道惯性坐标系下的位置和速度和t_k时刻目标航天器在追踪航天器轨道坐标系下的相对位置和相对速度的导航误差和按照如下关系计算t_k时刻相对导航滤波器的系统状态值和误差协方差矩阵

$X_{TC,Navi}^{t_k}=\left(\begin{array}{c}{r}_T^{t_k}-r_C^{t_k}+{\mathrm{\Delta r}}_{TC,NaviErr}^{t_k}\\ v_T^{t_k}-v_C^{t_k}+{\mathrm{\Delta v}}_{TC,NaviErr}^{t_k}\end{array}\right)$

$P_{TC,Navi}^{t_k}=P({\mathrm{\Delta r}}_{TC,NaviErr}^{t_k},{\mathrm{\Delta v}}_{TC,NaviErr}^{t_k})$

根据t_k时刻相对导航滤波器的系统状态值和误差协方差矩阵按照下述公式计算，得到一步递推的相对导航滤波器的系统状态值和误差协方差矩阵

$X_{TC,Navi}^{t_{k+1}/t_k}=\Phi \left(r_{C,Navi}^{t_k}\right)X_{TC,Navi}^{t_k}$

$P_{TC,Navi}^{t_{k+1}/t_k}=\Phi \left(r_{C,Navi}^{t_k}\right)P_{TC,Navi}^{t_k}\Phi {\left(r_{C,Navi}^{t_k}\right)}^T+\left(t_{k+1}-t_k\right)Q_{TC,Navi}$

其中，Q_TC，Navi为相对导航滤波器中的系统误差方差阵，具体表示为

$Q_{TC,Navi}=\left(\begin{array}{cc}{0}_{3\times 3}& 0_{3\times 3}\\ 0_{3\times 3}& {0.003}^2\times I_{3\times 3}\end{array}\right)$

为t_k时刻进行滤波估计一步递推计算的状态转移矩阵，具体为

$\Phi \left(r_{TC,Navi}^{t_k}\right)=\left(\begin{array}{cc}{I}_{3\times 3}& I_{3\times 3}\left(t_{k+1}-t_k\right)\\ \left(\frac{3\mu }{|r_{C,Navi}^{t_k}{|}^5}{r}_{C,Navi}^{t_k}{\left(r_{C,Navi}^{t_k}\right)}^T-\frac{\mu }{|r_{C,Navi}^{t_k}{|}^3}{I}_{3\times 3}\right)\left(t_{k+1}-t_k\right)& I_{3\times 3}\end{array}\right)$

μ为地球引力常数，取为3.986005×10¹⁴m³/s²；

根据一步递推计算的结果完成测量更新得到t_k+1时刻相对导航滤波器的系统状态值和误差协方差矩阵

$K_{TC,Navi}^{t_{k+1}/t_k}=P_{TC,Navi}^{t_{k+1}/t_k}{\left(H_{TC,Navi}^{t_{k+1}/t_k}\right)}^T{(H_{TC,Navi}^{t_{k+1}/t_k}{P}_{TC,Navi}^{t_{k+1}/t_k}{\left(H_{TC,Navi}^{t_{k+1}/t_k}\right)}^T+R_{TC,Navi})}^{-1}$

$X_{TC,Navi}^{t_{k+1}}=X_{TC,Navi}^{t_{k+1}/t_k}+K_{TC,Navi}^{t_{k+1}/t_k}(Y_{TC,Mea}^{t_{k+1}}-Y_{TC,Navi}^{t_{k+1}/t_k})$

$\left(\begin{array}{c}{P}_{TC,Navi}^{t_{k+1}}=\left(I_{6\times 6}-K_{TC,Navi}^{t_{k+1}/t_k}{H}_{TC,Navi}^{t_{k+1}/t_k}\right)P_{TC,Navi}^{t_{k+1}/t_k}{\left(I_{6\times 6}-K_{TC,Navi}^{t_{k+1}/t_k}{H}_{TC,Navi}^{t_{k+1}/t_k}\right)}^T\\ +K_{TC,Navi}^{t_{k+1}/t_k}{R}_{TC,Navi}^{t_{k+1}/t_k}{K}_{TC,Navi}^{t_{k+1}/t_k}\end{array}\right)$

其中，为相对导航的滤波增益矩阵；R_TC，Navi为相对导航滤波器中的测量误差方差阵，具体表示为

$R_{TC,Navi}=\left(\begin{array}{ccc}{10}^2& 0& 0\\ 0& {(0.01/180\pi )}^2& 0\\ 0& 0& {(0.01/180\pi )}^2\end{array}\right)$

为根据一步递推的状态计算得到的测量状态量，可以表示为

$Y_{TC,Navi}^{t_{k+1}/t_k}=\left(\begin{array}{c}{\rho }_{TC,Navi}^{t_{k+1}/t_k}\\ {\alpha }_{TC,Navi}^{t_{k+1}/t_k}\\ {\beta }_{TC,Navi}^{t_{k+1}/t_k}\end{array}\right)=\left(\begin{array}{c}\sqrt{{\left({\mathrm{\Delta x}}_{TC,Navi}^{t_{k+1}/t_k}\right)}^2+{\left({\mathrm{\Delta y}}_{TC,Navi}^{t_{k+1}/t_k}\right)}^2+{\left({\mathrm{\Delta z}}_{TC,Navi}^{t_{k+1}/t_k}\right)}^2}\\ \arctan 2\left(-{\mathrm{\Delta z}}_{TC,Navi}^{t_{k+1}/t_k},{\mathrm{\Delta x}}_{TC,Navi}^{t_{k+1}/t_k}\right)\\ \arctan \left(\frac{{\mathrm{\Delta y}}_{TC,Navi}^{t_{k+1}/t_k}}{\sqrt{{\left({\mathrm{\Delta x}}_{TC,Navi}^{t_{k+1}/t_k}\right)}^2+{\left({\mathrm{\Delta z}}_{TC,Navi}^{t_{k+1}/t_k}\right)}^2}}\right)\end{array}\right)$

为根据一步递推的状态计算得到的测量状态量的误差转移矩阵，可以表示为

$H_{TC,Navi}^{t_{k+1}/t_k}=\left(\begin{array}{cc}{H}_{11}\left({\mathrm{\Delta r}}_{TC,Navi}^{t_{k+1}/t_k}\right)L_{OI}\left(r_{C,Navi}^{t_k},v_{C,Navi}^{t_k}\right)& 0_{1\times 3}\\ H_{21}\left({\mathrm{\Delta r}}_{TC,Navi}^{t_{k+1}/t_k}\right)L_{OI}\left(r_{C,Navi}^{t_k},v_{C,Navi}^{t_k}\right)& 0_{1\times 3}\\ H_{31}\left({\mathrm{\Delta r}}_{TC,Navi}^{t_{k+1}/t_k}\right)L_{OI}\left(r_{C,Navi}^{t_k},v_{C,Navi}^{t_k}\right)& 0_{1\times 3}\end{array}\right)$

分别为误差转移矩阵的子矩阵，为根据t_k时刻追踪航天器在地心赤道惯性坐标系下的位置和速度的导航信息计算得到的由地心赤道惯性坐标系到追踪航天器轨道系的坐标转移矩阵，分别按照如下关系进行计算：

$H_{11}\left({\mathrm{\Delta r}}_{TC,Navi}^{t_{k+1}/t_k}\right)=\frac{\left(\begin{array}{ccc}{\mathrm{\Delta x}}_{TC,Navi}^{t_{k+1}/t_k}& {\mathrm{\Delta y}}_{TC,Navi}^{t_{k+1}/t_k}& {\mathrm{\Delta z}}_{TC,Navi}^{t_{k+1}/t_k}\end{array}\right)}{\left|{\mathrm{\Delta r}}_{TC,Navi}^{t_{k+1}/t_k}\right|}$

$H_{21}\left({\mathrm{\Delta r}}_{TC,Navi}^{t_{k+1}/t_k}\right)=\frac{\left(\begin{array}{ccc}{\mathrm{\Delta z}}_{TC,Navi}^{t_{k+1}/t_k}& 0& -{\mathrm{\Delta x}}_{TC,Navi}^{t_{k+1}/t_k}\end{array}\right)}{{\left({\mathrm{\Delta x}}_{TC,Navi}^{t_{k+1}/t_k}\right)}^2+{\left({\mathrm{\Delta z}}_{TC,Navi}^{t_{k+1}/t_k}\right)}^2}$

$H_{31}\left({\mathrm{\Delta r}}_{TC,Navi}^{t_{k+1}/t_k}\right)=\frac{\left(\begin{array}{ccc}-{\mathrm{\Delta x}}_{TC,Navi}^{t_{k+1}/t_k}{\mathrm{\Delta y}}_{TC,Navi}^{t_{k+1}/t_k}& {\left({\mathrm{\Delta x}}_{TC,Navi}^{t_{k+1}/t_k}\right)}^2+{\left({\mathrm{\Delta z}}_{TC,Navi}^{t_{k+1}/t_k}\right)}^2& -{\mathrm{\Delta y}}_{TC,Navi}^{t_{k+1}/t_k}{\mathrm{\Delta z}}_{TC,Navi}^{t_{k+1}/t_k}\end{array}\right)}{|{\mathrm{\Delta r}}_{TC,Navi}^{t_{k+1}/t_k}{|}^2\sqrt{{\left({\mathrm{\Delta x}}_{TC,Navi}^{t_{k+1}/t_k}\right)}^2+{\left({\mathrm{\Delta z}}_{TC,Navi}^{t_{k+1}/t_k}\right)}^2}}$

$L_{OI}\left(r_{C,Navi}^{t_k},v_{C,Navi}^{t_k}\right)={\left(\begin{array}{ccc}\frac{\left(r_{C,Navi}^{t_k}\times v_{C,Navi}^{t_k}\right)\times r_{C,Navi}^{t_k}}{|\left(r_{C,Navi}^{t_k}\times v_{C,Navi}^{t_k}\right)\times r_{C,Navi}^{t_k}|}& -\frac{r_{C,Navi}^{t_k}\times v_{C,Navi}^{t_k}}{|r_{C,Navi}^{t_k}\times v_{C,Navi}^{t_k}|}& -\frac{r_{C,Navi}^{t_k}}{\left|r_{C,Navi}^{t_k}\right|}\end{array}\right)}^T$

分别为根据一步递推的状态解算的目标航天器在追踪航天器轨道坐标系的相对位置，按照如下关系式进行计算

$\left(\begin{array}{c}{\mathrm{\Delta x}}_{TC,Navi}^{t_{k+1}/t_k}\\ {\mathrm{\Delta y}}_{TC,Navi}^{t_{k+1}/t_k}\\ {\mathrm{\Delta z}}_{TC,Navi}^{t_{k+1}/t_k}\end{array}\right)=\left(\begin{array}{cc}{L}_{OI}\left(r_{C,Navi}^{t_k},v_{C,Navi}^{t_k}\right)& 0_{3\times 3}\end{array}\right)X_{TC,Navi}^{t_{k+1}/t_k}$

d.相对导航误差计算

根据t_k+1时刻目标航天器和追踪航天器在地心赤道惯性坐标系下的以及t_k+1时刻相对导航滤波器的系统状态值计算t_k+1时刻目标航天器在追踪航天器轨道坐标系下的相对位置和相对速度的导航误差和

$\left(\begin{array}{c}{\mathrm{\Delta r}}_{TC,NaviErr}^{t_{k+1}}\\ {\mathrm{\Delta v}}_{TC,NaviErr}^{t_{k+1}}\end{array}\right)=X_{TC,Navi}^{t_{k+1}}-\left(\begin{array}{c}{r}_T^{t_{k+1}}-r_C^{t_{k+1}}\\ v_T^{t_{k+1}}-v_C^{t_{k+1}}\end{array}\right)$

所述步骤(四)轨道机动递推计算中，t_k时刻非线性的函数映射关系是指按照如下顺序关系(a导航信息计算；b轨道机动指令计算；c更新轨道信息；d相对导航误差计算)进行计算：

a.导航信息计算

根据t_k时刻追踪航天器在地心赤道惯性坐标系下的位置和速度按照如下关系计算t_k时刻追踪航天器在地心赤道惯性坐标系下位置和速度的导航信息

$r_{C,Navi}^{t_k}=r_C^{t_k}+r_{C,NaviErr}$

$v_{C,Navi}^{t_k}=v_C^{t_k}+v_{C,NaviErr}$

其中，r_C,NaviErr、v_C,NaviErr分别为追踪航天器位置速度的导航误差，满足均值为E(r_C,NaviErr)和E(v_C,NaviErr)、协方差为P(r_C,NaviErr)和P(v_C,NaviErr)的白噪声分布。

根据t_k时刻目标航天器和追踪航天器在地心赤道惯性坐标系下的位置和速度和t_k时刻目标航天器在追踪航天器轨道坐标系下的相对位置和相对速度的导航误差和按照如下关系计算t_k时刻相对导航滤波器的系统状态值

$X_{TC,Navi}^{t_k}=\left(\begin{array}{c}{r}_T^{t_k}-r_C^{t_k}+{\mathrm{\Delta r}}_{TC,NaviErr}^{t_k}\\ v_T^{t_k}-v_C^{t_k}+{\mathrm{\Delta v}}_{TC,NaviErr}^{t_k}\end{array}\right)$

b.轨道机动指令计算

根据t_k时刻相对导航滤波器的系统状态值以及t_k时刻追踪航天器在地心赤道惯性坐标系下位置和速度的导航信息按照如下公式计算t_k时刻目标航天器在追踪航天器轨道系下的相对位置和相对速度的导航信息

$\left(\begin{array}{c}{\mathrm{\Delta r}}_{TC,Navi}^{t_k}\\ {\mathrm{\Delta v}}_{TC,Navi}^{t_k}\end{array}\right)=\left(\begin{array}{cc}{L}_{OI}\left(r_{C,Navi}^{t_k},v_{C,Navi}^{t_k}\right)& 0_{3\times 3}\\ -L_{OI}\left(r_{C,Navi}^{t_k},v_{C,Navi}^{t_k}\right){\omega }^{\times }\left(r_{C,Navi}^{t_k},v_{C,Navi}^{t_k}\right)& L_{OI}\left(r_{C,Navi}^{t_k},v_{C,Navi}^{t_k}\right)\end{array}\right)X_{TC,Navi}^{t_k}$

其中，为采用t_k时刻追踪航天器在地心赤道惯性坐标系下位置和速度的导航信息计算得到的t_k时刻由地心赤道惯性坐标系到目标航天器的轨道坐标系的转移矩阵，具体计算公式为

$L_{OI}\left(r_{C,Navi}^{t_k},v_{C,Navi}^{t_k}\right)={\left(\begin{array}{ccc}\frac{\left(r_{C,Navi}^{t_k}\times v_{C,Navi}^{t_k}\right)\times r_{C,Navi}^{t_k}}{|\left(r_{C,Navi}^{t_k}\times v_{C,Navi}^{t_k}\right)\times r_{C,Navi}^{t_k}|}& -\frac{r_{C,Navi}^{t_k}\times v_{C,Navi}^{t_k}}{|r_{C,Navi}^{t_k}\times v_{C,Navi}^{t_k}|}& -\frac{r_{C,Navi}^{t_k}}{\left|r_{C,Navi}^{t_k}\right|}\end{array}\right)}^T$

为采用t_k时刻追踪航天器在地心赤道惯性坐标系下位置和速度的导航信息计算得到的t_k时刻地心赤道惯性坐标系下目标航天器的轨道角速度，具体计算公式为

$\omega (r_{C,Navi}^{t_k},v_{C,Navi}^{t_k})=\frac{r_{C,Navi}^{t_k}\times v_{C,Navi}^{t_k}}{\left|r_{C,Navi}^{t_k}\right|}$

为的叉乘矩阵。

根据t_k时刻目标航天器在追踪航天器轨道系下的相对位置和相对速度的导航信息按照如下CW制导律，计算t_k时刻施加给追踪航天器的速度脉冲指令

${\mathrm{\Delta v}}_{C,Cmd}^{t_k}=-{\left(B_{CW}\left({\mathrm{\Delta T}}_{Trans}\right)\right)}^{-1}\left(-{\mathrm{\Delta r}}_{CT,Des}-A_{CW}\left({\mathrm{\Delta T}}_{Trans}\right){\mathrm{\Delta r}}_{TC,Navi}^{t_k}\right)-{\mathrm{\Delta v}}_{TC,Navi}^{t_k}$

其中，Δr_CT,Des为进行CW制导律计算中追踪航天器在目标航天器轨道坐标系下的目标位置；A_CW(ΔT_Trans)、B_CW(ΔT_Trans)为CW制导律计算中的系数矩阵，是轨道机动转移时间ΔT_Trans的函数：

$A_{CW}\left({\mathrm{\Delta T}}_{Trans}\right)=\left(\begin{array}{ccc}1& 0& 6(n_{CW}{\mathrm{\Delta T}}_{Trans}-sin(n_{CW}{\mathrm{\Delta T}}_{Trans}\left)\right)\\ 0& cos\left(n_{CW}{\mathrm{\Delta T}}_{Trans}\right)& 0\\ 0& 0& 4-3cos\left(n_{CW}{\mathrm{\Delta T}}_{Trans}\right)\end{array}\right)$

$B_{CW}\left({\mathrm{\Delta T}}_{Trans}\right)=\left(\begin{array}{ccc}\frac{4\sin \left(n_{CW}{\mathrm{\Delta T}}_{Trans}\right)}{n_{CW}}-3{\mathrm{\Delta T}}_{Trans}& 0& \frac{2\left(1-\cos \left(n_{CW}{\mathrm{\Delta T}}_{Trans}\right)\right)}{n_{CW}}\\ 0& \frac{\sin \left(n_{CW}{\mathrm{\Delta T}}_{Trans}\right)}{n_{CW}}& 0\\ \frac{2\left(-1+\cos \left(n_{CW}{\mathrm{\Delta T}}_{Trans}\right)\right)}{n_{CW}}& 0& \frac{\sin \left(n_{CW}{\mathrm{\Delta T}}_{Trans}\right)}{n_{CW}}\end{array}\right)$

n_CW为CW制导律计算中使用的追踪航天器的轨道平均角速度。

c.更新轨道信息

根据t_k时刻目标航天器和追踪航天器在地心赤道惯性坐标系下的位置和速度和t_k时刻施加给追踪航天器的速度脉冲指令以及追踪航天器施加速度脉冲机动的误差Δv_C,CmdErr，按照如下公式更新t_k时刻目标航天器和追踪航天器在地心赤道惯性坐标系下的位置和速度和

$\left(\begin{array}{c}{r}_T^{\prime t_k}\\ v_T^{\prime t_k}\end{array}\right)=\left(\begin{array}{c}{r}_T^{t_k}\\ v_T^{t_k}\end{array}\right)$

$\left(\begin{array}{c}{r}_C^{\prime t_k}\\ v_C^{\prime t_k}\end{array}\right)=\left(\begin{array}{c}{r}_C^{t_k}\\ v_C^{t_k}\end{array}\right)+\left(\begin{array}{c}{0}_{3\times 1}\\ {\left(L_{OI}\right(r_C^{t_k},v_C^{t_k}\left)\right)}^T({\mathrm{\Delta v}}_{C,Cmd}^{t_k}+{\mathrm{\Delta v}}_{C,CmdErr})\end{array}\right)$

其中，Δv_C,CmdErr追踪航天器施加速度脉冲机动的误差Δv_C,CmdErr，满足均值为E(Δv_C,CmdErr)、协方差为P(Δv_C,CmdErr)的白噪声分布；为t_k时刻由地心赤道惯性坐标系到目标航天器的轨道坐标系的转移矩阵，具体计算公式为

$L_{OI}\left(r_C^{t_k},v_C^{t_k}\right)={\left(\begin{array}{ccc}\frac{\left(r_C^{t_k}\times v_C^{t_k}\right)\times r_C^{t_k}}{|\left(r_C^{t_k}\times v_C^{t_k}\right)\times r_C^{t_k}|}& -\frac{r_C^{t_k}\times v_C^{t_k}}{|r_C^{t_k}\times v_C^{t_k}|}& -\frac{r_C^{t_k}}{\left|r_C^{t_k}\right|}\end{array}\right)}^T$

根据t_k时刻相对导航滤波器的系统状态值t_k时刻施加给追踪航天器的速度脉冲指令以及采用t_k时刻追踪航天器在地心赤道惯性坐标系下位置和速度的导航信息计算得到的t_k时刻由地心赤道惯性坐标系到目标航天器的轨道坐标系的转移矩阵计算更新后的t_k时刻相对导航滤波器的系统状态值

$X_{TC,Navi}^{\prime t_k}=X_{TC,Navi}^{t_k}+\left(\begin{array}{c}{0}_{3\times 1}\\ -{\left(L_{OI}\right(r_{C,Navi}^{t_k},v_{C,Navi}^{t_k}\left)\right)}^T{\mathrm{\Delta v}}_{C,Cmd}^{t_k}\end{array}\right)$

d.相对导航误差计算

根据替换后的t_k时刻目标航天器和追踪航天器在地心赤道惯性坐标系下的以及替换后的t_k时刻相对导航滤波器的系统状态值计算替换后的t_k时刻目标航天器在追踪航天器轨道坐标系下的相对位置和相对速度的导航误差和

$\left(\begin{array}{c}{\mathrm{\Delta r}}_{TC,NaviErr}^{\prime t_k}\\ {\mathrm{\Delta v}}_{TC,NaviErr}^{\prime t_k}\end{array}\right)=X_{TC,Navi}^{\prime t_k}-\left(\begin{array}{c}{r}_T^{\prime t_k}-r_C^{\prime t_k}\\ v_T^{\prime t_k}-v_C^{\prime t_k}\end{array}\right).$

如图2所示为采用基于无迹递推的空间操作相对轨道控制误差计算方法得到的相对位置误差分布算例结果，由图可知，其中黑色椭圆为在任意选取的若干时刻，根据分析结果绘制的误差椭圆。如图3所示为基于蒙特卡洛计算的误差计算方法得到的相对位置误差分布算例结果，由图可知，其中浅灰色线条为对应的300次蒙特卡洛结果，黑色椭圆为根据300次蒙特卡洛仿真的统计结果绘制对应时刻的误差椭圆。

图2和图3表明，本发明所述的基于无迹变换的递推式误差计算方法可以得到与蒙特卡洛仿真方法相近的分析结果，但是本发明所述的方法仅需进行2n+1次非线性函数计算(n为系统变量的维数)。

本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。