CN104408326A - 一种对深空探测自主导航滤波算法的评估方法 - Google Patents

Abstract

本发明公开了一种对深空探测自主导航滤波算法的评估方法，先分别采用EKF、UKF、UPF三种导航滤波算法时对获取的导航数据进行处理，分别获取到对应导航滤波算法下的精度、可用性、连续性、实时性、稳定性等性能参数指标的评估值，综合评估模块再结合性能参数指标的评估值对深空探测自主导航滤波算法进行评估，这样综合评估模块不仅能够科学、精确地分配各种性能参数指标间的相对权重，而且还能准确的评估深空探测自主导航滤波算法，适用于当今深空探测技术领域发展的需要。

Description

一种对深空探测自主导航滤波算法的评估方法

技术领域

本发明属于深空探测技术领域，更为具体地讲，涉及一种对深空探测自主导航滤波算法的评估方法。

背景技术

人类在深空探测领域已取得了丰硕成果，我国深空探测的步伐正在加快，目前已制定完成了2030年前深空探测路线图，计划于2018年发射火星探测器、2019年发射小行星探测器。

与近地空间任务相比，深空探测具有飞行距离远、环境位置因素多、飞行程序复杂、器地通信时延大、自主性要求高等特点，其中自主导航能力对任务的完成至关重要。在导航方法和理论的研究和设计过程中，对性能的评估和验证是一个非常重要的环节。

限于成本、安全性等诸多因素的考虑，深空探测相关实验只能通过模拟测试，而很难进行实际测量。常用的方法主要分为以下几类：

1.理论分析：从深空探测的机理上进行定性分析，从飞行原理、飞行器的轨道动力学模型、姿态动力学模型以及运动和控制机理等方面进行相应公式的计算和推导，对其各项性能进行假设和演算，从原理上验证导航方法是否可行。

2.软件仿真：根据系统几何模型，将整个动态过程用程序语言进行实现。常用的程序语言包括MATLAB、C++/C#等，其中MATLAB针对数学计算进行了优化，操作简便，并能够更加直观的将运算结果以图表、曲线等形式进行展示。C++/C#较为灵活，可以根据需要的不同来设置更加复杂的环境条件，并能够更加方便的对导航系统运行过程进行控制。由于所有的运行条件、扰动和各种参数都是人为设置的，所以程序的运行过程与实际环境的接近程度就成为了该种方法的主要限制因素。

3.半物理仿真：也即半实物仿真。针对仿真研究内容，将被仿真对象(导航系统)的一部分以实物(或物理模型)的方式引入仿真回路，被仿真对象系统的其余部分以数学模型描述，并把它转化为仿真计算模型。借助物理效应模型，进行实时数学仿真与物理仿真的联合仿真。半物理仿真的逼真度较高，其逼真度取决于接入的实物部件的多寡、仿真计算机的速度、精度和功能，转台和各目标模拟器的性能等。

4.实地测试：完全采用实际的导航算法、飞行器等，并将被测试的系统放入真实的环境中，通过系统的运行结果进行导航方法的验证。但在大多数情况下，我们无法采用实际的飞行测试，因而这种方法很少被采用。

发明内容

本发明的目的在于克服现有技术的不足，提供一种对深空探测自主导航滤波算法的评估方法，结合精度、可用性、连续性、实时性、稳定性等性能参数评估值对导航滤波算法进行评估，这样能够科学、精确地分配各种指标的相对权重，从而准确的评估深空探测自主导航滤波算法。

为实现上述发明目的，本发明一种对深空探测自主导航滤波算法的评估方法，其特征在于，包括以下步骤：

(1)、深空探测自主导航滤波算法采用EKF滤波时，对获取得到的导航数据进行处理

计算探测器在第i个采样点的位置差值ep(i)：

ex(i)＝Sat_p(x,i)-sat_p(x,i)；

ey(i)＝Sat_p(y,i)-sat_p(y,i)；

ez(i)＝Sat_p(z,i)-sat_p(z,i)；

$\mathrm{ep}\left(i\right)=\sqrt{(\mathrm{ex}{\left(i\right)}^2+\mathrm{ey}{\left(i\right)}^2+\mathrm{ez}{\left(i\right)}^2)};$

其中，i＝1,2,...,k，k表示采样点的个数，Sat_p(x,i)、Sat_p(y,i)、Sat_p(z,i)分别表示探测器在x、y、z三方向上的位置观测值，sat_p(x,i)、sat_p(y,i)、sat_p(z,i)分别表示探测器在x、y、z三个方向上位置的滤波优化值，ex(i)、ey(i)、ez(i)分别表示探测器在x、y、z三个方向上位置的差值；

计算探测器在第i个采样点的速度差值ev(i)：

evx(i)＝Sat_v(x,i)-sat_v(x,i)；

evy(i)＝Sat_v(y,i)-sat_v(y,i)；

ezv(i)＝Sat_v(z,i)-sat_v(z,i)；

$\mathrm{ev}\left(i\right)=\sqrt{(\mathrm{evx}{\left(i\right)}^2+\mathrm{evy}{\left(i\right)}^2+\mathrm{evz}{\left(i\right)}^2)};$

其中，Sat_v(x,i)、Sat_v(y,i)、Sat_v(z,i)分别表示探测器在x、y、z三方向上的速度观测值，sat_v(x,i)、sat_v(y,i)、sat_v(z,i)分别表示探测器在x、y、z三个方向上速度的滤波优化值，evx(i)、evy(i)、evz(i)分别表示探测器在x、y、z三个方向上速度的差值；

(2)、获取EKF滤波时参数指标的评估值

(2.1)、利用前向均方差计算EKF滤波时精度评估值

计算探测器在第i采样点时刻对应EKF滤波的位置精度：

$R_p=\sqrt{\frac{1}{i}\underset{i=1}{\overset{i}{\mathrm{\Σ}}}\mathrm{ep}{\left(i\right)}^2}$

计算探测器在第i采样点时刻对应EKF滤波的速度精度：

$R_v=\sqrt{\frac{1}{i}\underset{i=1}{\overset{i}{\mathrm{\Σ}}}\mathrm{ev}{\left(i\right)}^2}$

(2.2)、计算EKF滤波时可用性评估值

将探测器在第i个采样点的位置差值ep(i)与预设的阈值M₁比较，当ep(i)＜M₁时，则可用次数sum_a加1，如果ep(i)≥M₁，则可用次数sum_a保持不变，直到所有采样点比较完；从而计算出EKF滤波时可用性评估值K_s为：

$K\_s=\frac{\mathrm{sum}\_a}{k};$

(2.3)、计算EKF滤波时连续性评估值

计算探测器在相邻采样点间的位置差值c(i)：

c(i)＝|ep(i)-ep(i-1)|

在第i个采样点时刻，将得到的c(i)与预设的阀值M₂比较，如果c(i)＜M₂，则连续次数sum加1，如果c(i)≥M₂，则连续次数sum保持不变，直到所有采样点比较完；从而计算出EKF滤波时连续性评估值P为：

(2.4)、计算EKF滤波时实时性评估值

EKF导航滤波程序运行时间记为t₀，EKF导航滤波程序运行结束时间记为t₁，则EKF滤波时实时性评估值T为：

T＝t₁-t₀

(2.5)、通过计算EKF滤波的信噪比确定EKF滤波时稳定性评估值

SNR＝S/N

其中，S表示EKF滤波处理外部设备输入的电子信号，N表示原信号中不存在的无规则额外信号；

当对采用EKF滤波时获取得到的导航数据处理完成后，按照步骤(1)和步骤(2)的方法，继续处理采用UKF、UPF滤波时获取得到的导航数据，分别获取到UKF滤波和UPF滤波下参数指标的评估值；

(3)、对深空探测自主导航滤波算法进行综合评估

(3.1)、构造准则层判断矩阵

根据9/9-9/1标度层次分析法，综合评估模块构造准则层判断矩阵A＝a_i'j'，且a_i'j'＞0,a_j'i'＝1/a_i'j'，其中，i',j'∈[1,5]，判断矩阵A是5个参数指标值映射出来的5×5的矩阵，a_i'j'表示映射的第i'个参数指标和映射的第j'个参数指标的重要性的比值，a_i'j'越大，第i'个参数指标较第j'个参数指标更重要；

计算判断矩阵A的特征值和最大特征根：

计算准则层判断矩阵A每一行比值的乘积M_i'

$M_{i^{\′}}=\underset{j^{\′}=1}{\overset{5}{\mathrm{\Π}}}{a}_{i^{\′}{j}^{\′}}$

计算M_i'的5次方根V_i'

$V_{i^{\′}}=5\sqrt{M_{i^{\′}}}$

对V_i'进行归一化，即：

$F_{i^{\′}}=V_{i^{\′}}/{\mathrm{\Σ}}_{i^{\′}=1}^5{V}_{i^{\′}}$

则特征向量F＝(F₁,F₂,...,F₅)^T；

计算准则层判断矩阵A的最大特征根λ_max

${\mathrm{\λ}}_{\max }=\underset{i^{\′}=1}{\overset{5}{\mathrm{\Σ}}}\frac{{(A\·F)}_{i^{\′}}}{{5F}_{i^{\′}}}$

验证准则层判断矩阵A的一致性：

计算准则层判断矩阵A的一致性比率C_R，

其中，C_I为一致性检验指标，R_I＝1.12

当C_R＜0.1时，该准则层判断矩阵A为一致矩阵，否则，返回步骤(3.1)调整a_i'j'的大小，使准则层判断矩阵A最终满足一致矩阵；

(3.2)、将步骤(2)得到的评估值全部输入到综合评估模块，构造方案层判断矩阵B_j

将导航滤波算法EKF、UKF、UPF在每一个相同参数指标下构造一个3×3判断矩阵B_j，j＝1,2,...,5表示参数指标的个数，且其中,判断矩阵B_j是第j个参数指标对应三种导航滤波算法映射出来的3×3的矩阵，表示在该参数指标下映射的第个导航滤波算法的评估值同映射的第个导航滤波算法的评估值的比值，越大，在该参数指标下，映射的第个导航滤波算法的效果越好；

计算方案层判断矩阵B₁的特征值和最大特征根：

计算方案层判断矩阵B₁每一行比值的乘积

$N_{\tilde{i}}=\underset{\tilde{j}}{\overset{3}{\mathrm{\Π}}}{b}_{\tilde{i}\tilde{j}}$

计算的3次方根

$H_{\tilde{i}}=3\sqrt{N_{\tilde{i}}}$

对进行归一化处理，即：

$S_{\tilde{i}1}=H_{\tilde{i}}/{\mathrm{\Σ}}_{\tilde{i}=1}^3{H}_{\tilde{i}}$

则特征向量SS₁＝(S₁₁,S₂₁,S₃₁)^T；

计算判断矩阵的最大特征根λ_max1:

${\mathrm{\λ}}_{\max 1}=\underset{\tilde{i}=1}{\overset{3}{\mathrm{\Σ}}}\frac{{(B_1\·{\mathrm{SS}}_1)}_{\tilde{i}}}{{3S}_{\tilde{i}1}}$

验证方案层判断矩阵B₁的一致性：

计算方案层判断矩阵B₁的一致性比率

其中，C_I1为一致性检验指标，R_I1＝0.58

当C_R1＜0.1时，该方案层判断矩阵B₁为一致矩阵，否则，返回步骤(3.2)调整的大小，使准则层判断矩阵B₁最终满足一致矩阵；

同理，可以得到方案层判断矩阵B₂～B₅的特征向量SS₂～SS₅和最大特征根λ_max2～λ_max5；

用方案层判断矩阵B_j的所有特征向量SS₁～SS₅构造矩阵C＝(SS₁,SS₂,...,SS₅)；(3.3)、深空探测自主导航滤波算法的综合评估

将步骤(3.2)求得的特征向量C与步骤(3.1)构造的矩阵F相乘，求出权值W，即：W＝C·F，分别计算出权值W中每一行对应的权值其中，表示滤波算法的个数，是矩阵C中第行、第j列个元素；对所有权值进行总排序，其权值最大项即为对应的导航滤波方法最好。

本发明的发明目的是这样实现的：

本发明一种对深空探测自主导航滤波算法的评估方法，先分别采用EKF、UKF、UPF三种导航滤波算法时对获取的导航数据进行处理，分别获取到对应导航滤波算法下的精度、可用性、连续性、实时性、稳定性等性能参数指标的评估值，综合评估模块再结合性能参数指标的评估值对深空探测自主导航滤波算法进行评估，这样综合评估模块不仅能够科学、精确地分配各种性能参数指标间的相对权重，而且还能准确的评估深空探测自主导航滤波算法，适用于当今深空探测技术领域发展的需要。

同时，本发明一种对深空探测自主导航滤波算法的评估方法还具有以下有益效果：

(1)、本发明对深空探测自主导航滤波算法的评估比较全面且精准。本发明中结合导航滤波算法的精度、可用性、连续性、实时性、稳定性等性能参数对深空探测自主导航滤波算法进行评估，保证了评估的全面性和准确性，从而克服了现有方法的主观性强，随意性大，精度不够，性能参数指标间的相对权重不合理等问题。

(2)、增加了基于层次分析的综合评估方法，这样综合评估模块能够科学、精确地分配各种性能参数间相对权重。

附图说明

图1是本发明对深空探测自主导航滤波算法的评估的流程图；

图2是EKF、UKF、UPF三种导航滤波算法的精度评估结果仿真图；

图3是探测器的位置误差分布图；

图4是导航系统在巡航阶段和捕获阶段的层次结构模型图；

表1是准则层判断矩阵A；

表2是方案层判断矩阵B₁；

表3是权值的单项排序表。

具体实施方式

下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和设计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。

实施例

为了方便描述，先对具体实施方式中出现的相关专业术语进行说明：

EKF(Extended Kalman Filter)：扩展卡尔曼滤波；

UKF(Unscented Kalman Filter)：无损卡尔曼滤波、无迹卡尔曼滤波或者去芳香卡尔曼滤波；

UPF(Unscented particle filter)：无迹粒子滤波；

由于这是深空探测自主导航滤波算法的评估，所以我们不可能得到真实深空导航数据，因此，我们采用STK和MATLAB混合编程得到仿真数据Satellite_p、Satellite_v。STK提供太阳、地球、火星及其卫星、目标恒星和目标小行星的速度和位置，用MATLAB仿真天文测角导航、红移测速导航以及测速测角组合导航过程，我们得到探测器的位置、速度估计值satellite_p、satellite_v，也即我们所需要的导航结果；本实施例中，采用了三种导航滤波算法EKF、UKF、UPF，分别对获取得到的导航数据进行处理，再计算出三种导航滤波算法EKF、UKF、UPF下的五个性能参数指标的评估值，然后输入到综合评估模块，通过综合评估模块进行评估，我们由此来判断导航滤波算法的性能。

在下文中，Satellite_p简写为Sat_p，Satellite_v简写为Sat_v。

在本实施例中，如图1所示，本发明一种对深空探测自主导航滤波算法的评估方法，包括以下步骤：

S1、深空探测自主导航滤波方算法采用EKF滤波时，对获取得到的导航数据进行处理

计算探测器在第i个采样点的位置差值ep(i)：

ex(i)＝Sat_p(x,i)-sat_p(x,i)；

ey(i)＝Sat_p(y,i)-sat_p(y,i)；

ez(i)＝Sat_p(z,i)-sat_p(z,i)；

$\mathrm{ep}\left(i\right)=\sqrt{(\mathrm{ex}{\left(i\right)}^2+\mathrm{ey}{\left(i\right)}^2+\mathrm{ez}{\left(i\right)}^2)};$

其中，i＝1,2,...,k，k表示采样点的个数，Sat_p(x,i)、Sat_p(y,i)、Sat_p(z,i)分别表示探测器在x、y、z三方向上的位置观测值，sat_p(x,i)、sat_p(y,i)、sat_p(z,i)分别表示探测器在x、y、z三个方向上位置的滤波优化值，ex(i)、ey(i)、ez(i)分别表示探测器在x、y、z三个方向上位置的差值；

计算探测器在第i个采样点的速度差值ev(i)：

evx(i)＝Sat_v(x,i)-sat_v(x,i)；

evy(i)＝Sat_v(y,i)-sat_v(y,i)；

ezv(i)＝Sat_v(z,i)-sat_v(z,i)；

$\mathrm{ev}\left(i\right)=\sqrt{(\mathrm{evx}{\left(i\right)}^2+\mathrm{evy}{\left(i\right)}^2+\mathrm{evz}{\left(i\right)}^2)};$

其中，Sat_v(x,i)、Sat_v(y,i)、Sat_v(z,i)分别表示探测器在x、y、z三方向上的速度观测值，sat_v(x,i)、sat_v(y,i)、sat_v(z,i)分别表示探测器在x、y、z三个方向上速度的滤波优化值，evx(i)、evy(i)、evz(i)分别表示探测器在x、y、z三个方向上速度的差值；

本实施例中，还可以利用损失函数计算探测器在第i个采样点的位置差值ep(i)'和速度误差ev(i)'，具体计算如下：

计算探测器在第i个采样点的位置差值ep(i)'：

将ex(i)与预设的位置判断值beta_p比较，如果ex(i)＜beta_p，则ex(i)保持不变，即ex(i)'＝ex(i)；如果ex(i)≥beta_p，则ex(i)'＝2beta_p·|ex(i)|-beta_p²；

同理再分别将ey(i)、ez(i)与预设的位置判断值beta_p比较，得到最终的ey(i)'、ez(i)'；

则探测器在第i个采样点的位置差值ep(i)'为：

$\mathrm{ep}{\left(i\right)}^{\′}=\sqrt{{\left[\mathrm{ex}{\left(i\right)}^{\′}\right]}^2+{\left[\mathrm{ey}{\left(i\right)}^{\′}\right]}^2+{\left[\mathrm{ez}{\left(i\right)}^{\′}\right]}^2}$

计算探测器在第i个采样点的速度误差ev(i)'：

将evx(i)与预设的位置判断值beta_v比较，如果evx(i)＜beta_v，则evx(i)保持不变，即evx(i)'＝evx(i)；如果evx(i)≥beta_v，则evx(i)'＝2beta_v·|evx(i)|-beta_v²；

同理再分别将evy(i)、evz(i)与预设的位置判断值beta_v比较，得到最终的evy(i)'、evz(i)'；

则探测器在第i个采样点的速度误差ev(i)'为：

$\mathrm{ev}{\left(i\right)}^{\′}=\sqrt{{\left[\mathrm{evx}{\left(i\right)}^{\′}\right]}^2+{\left[\mathrm{evy}{\left(i\right)}^{\′}\right]}^2+{\left[\mathrm{evz}{\left(i\right)}^{\′}\right]}^2}$

S2、获取EKF滤波时参数指标的评估值

S2.1)、利用前向均方差计算EKF滤波时精度评估值

计算探测器在第i采样点时刻对应EKF滤波的位置精度：

$R_p=\sqrt{\frac{1}{i}\underset{i=1}{\overset{i}{\mathrm{\Σ}}}\mathrm{ep}{\left(i\right)}^2}$

计算探测器在第i采样点时刻对应EKF滤波的速度精度：

$R_v=\sqrt{\frac{1}{i}\underset{i=1}{\overset{i}{\mathrm{\Σ}}}\mathrm{ev}{\left(i\right)}^2}$

本实施例中，采用EKF、UKF、UPF三种导航滤波算法时，对计算得到的导航参数指标的评估值输入到综合评估模块，其评估结果如图2所示，其中实线、空心线、心形线曲线分别表示EKF、UKF、UPF的评估结果。由图可知，EKF方法的误差较大，UKF的误差其次，UPF的误差最小，具体的讲，UPF曲线有向下的趋势，所以其导航滤波算法是可以稳定的，而UKF、EKF的误差也一直保持在10×10^-3～11×10^-3范围内，故其导航滤波算法也是稳定可行的。

S2.2)、计算EKF滤波时可用性评估值

首先，我们对可用性做介绍，根据可用性的定义可表示为：

其中，K表示可用性，L表示整个导航滤波过程，而l表示其中可用的部分；

以时间比例为基准，导航滤波运行时，导航滤波算法的可用性可以通过以下方法计算，即：

将探测器在第i个采样点的位置差值ep(i)与预设的阈值M₁比较，当ep(i)＜M₁时，则可用次数sum_a加1，如果ep(i)≥M₁，则可用次数sum_a保持不变，直到所有采样点比较完；从而计算出EKF滤波时可用性评估值K_s为：

$K\_s=\frac{\mathrm{sum}\_a}{k};$

在本实施例中，需要将探测器在每个采样点的位置差值与预设的阈值M₁比较，因此在整个导航滤波过程中的取值为采样点的个数k；

S2.3)、计算EKF滤波时连续性评估值

探测器在实际的飞行过程中，其位置、速度都是连续变化的，反应到函数图像上，也就是其位置、速度曲线都是连续、平滑的曲线。因而，如果经过导航算法计算后得到的位置、速度曲线突然出现了较大幅度的变化，则必然是受到了外界扰动的影响，使得其内部状态出现了不连续的情况。所以我们可以用如下的方法来计算导航滤波算法的连续性。

计算探测器相邻采样点间的位置差值c(i)：

c(i)＝|ep(i)-ep(i-1)|

在第i个采样点时刻，将得到的c(i)与预设的阈值M₂比较，如果c(i)＜M₂，则连续次数sum加1，如果c(i)≥M₂，则连续次数sum保持不变，直到所有采样点比较完；从而计算出EKF滤波时连续性评估值P为：

在本实施例中，在步骤S1中，导航滤波算法计算出的位置误差位于一个动态的范围之内，如图3所示，从所有误差的均值朝着上下两个方向扩展，可以得到一个范围区间，这个区间包含了绝大多数(95％)的误差点，则这个范围的大小就是导航滤波算法的连续性指标；

S2.4)、计算EKF滤波时实时性评估值

本实施例中，导航滤波算法必须满足实时性的要求，也即导航滤波程序的运行时间必须控制在一个合理的范围内，因而实时性可以用导航滤波程序的运行时间来进行评估；

EKF导航滤波程序运行时间记为t₀，EKF导航滤波程序运行结束时间记为t₁，则EKF滤波时实时性评估值T为：

T＝t₁-t₀

S2.5)、通过计算EKF滤波的信噪比确定EKF滤波时稳定性评估值

SNR＝S/N

其中，S表示EKF滤波处理外部设备输入的电子信号，N表示原信号中不存在的无规则额外信号；

在本实施例中，额外信号N可以分为内部和外部两种。内部信号主要是由于电路设计、制造工艺等因素，由设备自身产生的，而外部信号是由设备所在的电子环境和物理化学环境(自然环境)所造成的；

当对采用EKF滤波时获取得到的导航数据处理完成后，按照步骤(1)和步骤(2)的方法，继续处理采用UKF、UPF滤波时获取得到的导航数据，分别获取到UKF滤波和UPF滤波下参数指标的评估值；

S3、对深空探测自主导航滤波算法进行综合评估

在本实施例中，以导航的巡航阶段和捕获阶段为例，需要考虑精度、可用性、连续性、稳定性、实时性等参数指标，利用层次分析法最终把系统分析归结为最底层相对于最高层的相对重要性权值的确定或相对好坏次序的排序问题。巡航阶段和捕获阶段的层次结构模型如图4所示，下面进行具体说明，如下：

S3.1)、构造准则层判断矩阵

首先，综合评估模块根据9/9-9/1标度层次分析法构造准则层判断矩阵A，A＝a_i'j'，且a_i'j'＞0,a_j'i'＝1/a_i'j'，其中，i',j'∈[1,5]，判断矩阵A是5个参数指标值映射出来的5×5的矩阵，a_i'j'表示映射的第i'个参数指标和映射的第j'个参数指标的重要性的比值，a_i'j'越大，第i'个参数指标较第j'个参数指标更重要，具体如表1所示；

	精度	可用性	连续性	实时性	稳定性
						精度	1	1.286	1.8	3	3
可用性	0.776	1	1.8	3	3
						连续性	0.555	0.555	1	1.125	1.125
实时性	0.333	0.333	0.888	1	1
						稳定性	0.333	0.333	0.888	1	1

表1

计算判断矩阵A的特征值和最大特征根：

计算准则层判断矩阵A每一行比值的乘积M_i'

$M_{i^{\′}}=\underset{j^{\′}=1}{\overset{5}{\mathrm{\Π}}}{a}_{i^{\′}{j}^{\′}}$

计算M_i'的5次方根V_i'

$V_{i^{\′}}=5\sqrt{M_{i^{\′}}}$

对V_i'进行归一化，即：

$F_{i^{\′}}=V_{i^{\′}}/{\mathrm{\Σ}}_{i^{\′}=1}^5{V}_{i^{\′}}$

则特征向量F＝(F₁,F₂,...,F₅)^T；

计算准则层判断矩阵A的最大特征根λ_max

${\mathrm{\λ}}_{\max }=\underset{i^{\′}=1}{\overset{5}{\mathrm{\Σ}}}\frac{{(A\·F)}_{i^{\′}}}{{5F}_{i^{\′}}}$

验证准则层判断矩阵A的一致性：

计算准则层判断矩阵A的一致性比率C_R，

其中，C_I为一致性检验指标，R_I＝1.12

当C_R＜0.1时，该准则层判断矩阵A为一致矩阵，否则，返回步骤S3.1)调整a_i'j'的大小，使准则层判断矩阵A最终满足一致矩阵；

S3.2)、将步骤S2得到的评估值全部输入到综合评估模块，构造方案层判断矩阵B_j

将导航滤波算法EKF、UKF、UPF在每一个相同参数指标下构造一个3×3判断矩阵B_j，j＝1,2,...,5表示参数指标的个数，且其中,判断矩阵B_j是第j个参数指标对应三种导航滤波算法映射出来的3×3的矩阵，表示在该参数指标下映射的第个导航滤波算法的评估值同映射的第个导航滤波算法的评估值的比值，越大，在该参数指标下，映射的第个导航滤波算法的效果越好；

本实施例中，导航滤波算法采用EKF、UKF、UPF三种滤波，接下来利用三种滤波算法得到的五个参数指标((精度、可用性、实时性、稳定性、连续性)的评估值分别来构造5个方案层判断矩阵。

以精度的评估值为例，构造方案层判断矩阵B₁，如表2所示，若UKF的精度评估值是EKF的精度评估值3倍，则取b₁₂＝3，而 的其他取值同理，从而构造出方案层判断矩阵B₁；

精度	EKF	UKF	UPF
				EKF	b11	b12	b13
UKF	b21	b22	b23
				UPF	b31	b32	b33

表2

同理，利用其他四个参数指标((可用性、实时性、稳定性、连续性)的评估值来构造方案层判断矩阵B₂～B₅；

本实施例中，继续以方案层判断矩阵B₁为例，计算出B₁的特征值和最大特征根：

计算方案层判断矩阵B₁每一行比值的乘积

$N_{\tilde{i}}=\underset{\tilde{j}}{\overset{3}{\mathrm{\Π}}}{b}_{\tilde{i}\tilde{j}}$

计算的3次方根

$H_{\tilde{i}}=3\sqrt{N_{\tilde{i}}}$

对进行归一化处理，即：

$S_{\tilde{i}1}=H_{\tilde{i}}/{\mathrm{\Σ}}_{\tilde{i}=1}^3{H}_{\tilde{i}}$

则特征向量SS₁＝(S₁₁，S₂₁，S₃₁)^T；

计算判断矩阵的最大特征根λ_max1：

${\mathrm{\λ}}_{\max 1}=\underset{\tilde{i}=1}{\overset{3}{\mathrm{\Σ}}}\frac{{(B_1\·{\mathrm{SS}}_1)}_{\tilde{i}}}{{3S}_{\tilde{i}1}}$

验证方案层判断矩阵B₁的一致性：

计算方案层判断矩阵B₁的一致性比率C_R1，

其中，C_I1为一致性检验指标，R_I1＝0.58

当C_R1＜0.1时，该方案层判断矩阵B₁为一致矩阵，否则，返回步骤S3.2)调整的大小，使准则层判断矩阵B₁最终满足一致矩阵；

同理，可以得到其他的方案层判断矩阵B₂～B₅的特征向量SS₂～SS₅和最大特征根λ_max2～λ_max5；

用方案层判断矩阵B_j的所有特征向量SS₁～SS₅构造矩阵C＝(SS₁,SS₂,...,SS₅)；

S3.3)、深空探测自主导航滤波算法的综合评估

将步骤(3.2)求得的特征向量C与步骤(3.1)构造的矩阵F相乘，求出权值W，即：W＝C·F，分别计算出权值W中每一行对应的权值，其中，表示滤波算法的个数，对所有权值进行总排序，其权值最大项即为对应的导航滤波方法最好。

本实施例中，如表3所示，将精度判断矩阵的特征向量SS₁＝(S₁₁,S₂₁,S₃₁)^T作为元素S₁的层次单排序权值，将可用性判断矩阵的特征向量SS₂＝(S₁₂,S₂₂,S₃₂)^T作为对于元素S₂的层次单排序权值，……，将稳定性判断矩阵的特征向量SS₅＝(S₁₅,S₂₅,S₃₅)^T作为对于元素S₅的层次单排序权值；分别计算出权值W中每一行对应的权值 是矩阵C中第行、第j列个元素。

表3

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

Claims (2)

1.一种对深空探测自主导航滤波算法的评估方法，其特征在于，包括以下步骤：

(1)、深空探测自主导航滤波算法采用EKF滤波时，对获取得到的导航数据进行处理

计算探测器在第i个采样点的位置差值ep(i)：

ex(i)＝Sat_p(x,i)-sat_p(x,i)；

ey(i)＝Sat_p(y,i)-sat_p(y,i)；

ez(i)＝Sat_p(z,i)-sat_p(z,i)；

$\mathrm{ep}\left(i\right)=\sqrt{(\mathrm{ex}{\left(i\right)}^2+\mathrm{ey}{\left(i\right)}^2+\mathrm{ez}{\left(i\right)}^2)};$

其中，i＝1,2,…,k，k表示采样点的个数，Sat_p(x,i)、Sat_p(y,i)、Sat_p(z,i)分别表示探测器在x、y、z三方向上的位置观测值，sat_p(x,i)、sat_p(y,i)、sat_p(z,i)分别表示探测器在x、y、z三个方向上，位置的滤波优化值，ex(i)、ey(i)、ez(i)分别表示探测器在x、y、z三个方向上位置的差值；

计算探测器在第i个采样点的速度差值ev(i)：

evx(i)＝Sat_v(x,i)-sat_v(x,i)；

evy(i)＝Sat_v(y,i)-sat_v(y,i)；

ezv(i)＝Sat_v(z,i)-sat_v(z,i)；

$\mathrm{ev}\left(i\right)=\sqrt{(\mathrm{evx}{\left(i\right)}^2+\mathrm{evy}{\left(i\right)}^2+\mathrm{evz}{\left(i\right)}^2)};$

其中，Sat_v(x,i)、Sat_v(y,i)、Sat_v(z,i)分别表示探测器在x、y、z三方向上的速度观测值，sat_v(x,i)、sat_v(y,i)、sat_v(z,i)分别表示探测器在x、y、z三个方向上，速度的滤波优化值，evx(i)、evy(i)、evz(i)分别表示探测器在x、y、z三个方向上速度的差值；

(2)、获取EKF滤波时参数指标的评估值

(2.1)、利用前向均方差计算EKF滤波时精度评估值

计算探测器在第i采样点时刻对应EKF滤波的位置精度：

$R_p=\sqrt{\frac{1}{i}\underset{i=1}{\overset{i}{\mathrm{\Σ}}}\mathrm{ep}{\left(i\right)}^2}$

计算探测器在第i采样点时刻对应EKF滤波的速度精度：

$R_v=\sqrt{\frac{1}{i}\underset{i=1}{\overset{i}{\mathrm{\Σ}}}\mathrm{ev}\left(i\right)}$

(2.2)、计算EKF滤波时可用性评估值

将探测器在第i个采样点的位置差值ep(i)与预设的阈值M₁比较，当ep(i)＜M₁时，则可用次数sum_a加1，如果ep(i)≥M₁，则可用次数sum_a保持不变，直到所有采样点比较完；从而计算出EKF滤波时可用性评估值K_s为：

$K\_s=\frac{\mathrm{sum}\_a}{k};$

(2.3)、计算EKF滤波时连续性评估值

计算探测器在相邻采样点间的位置差值c(i)：

c(i)＝|ep(i)-ep(i-1)|

在第i个采样点时刻，将得到的c(i)与预设的阀值M₂比较，如果c(i)＜M₂，则连续次数sum加1，如果c(i)≥M₂，则连续次数sum保持不变，直到所有采样点比较完；从而计算出EKF滤波时连续性评估值P为：

(2.4)、计算EKF滤波时实时性评估值

EKF导航滤波程序运行时间记为t₀，EKF导航滤波程序运行结束时间记为t₁，则EKF滤波时实时性评估值T为：

T＝t₁-t₀

(2.5)、通过计算EKF滤波的信噪比确定EKF滤波时稳定性评估值

SNR＝S/N

其中，S表示EKF滤波处理外部设备输入的电子信号，N表示原信号中不存在的无规则额外信号；

当对采用EKF滤波时获取得到的导航数据处理完成后，按照步骤(1)和步骤(2)的方法，继续处理采用UKF、UPF滤波时获取得到的导航数据，分别获取到UKF滤波和UPF滤波下参数指标的评估值；

(3)、对深空探测自主导航滤波算法进行综合评估

(3.1)、构造准则层判断矩阵

根据9/9-9/1标度层次分析法，综合评估模块构构造准则层判断矩阵A＝a_i'j'，且a_i'j'＞0,a_j'i'＝1/a_i'j'，其中，i',j'∈[1,5]，判断矩阵A是5个参数指标值映射出来的5×5的矩阵，a_i'j'表示映射的第i'个参数指标和映射的第j'个参数指标的重要性的比值，a_i'j'越大，第i'个参数指标较第j'个参数指标更重要；

计算判断矩阵A的特征值和最大特征根：

计算准则层判断矩阵A每一行比值的乘积M_i'

$M_{i^{\′}}=\underset{j^{\′}=1}{\overset{5}{\mathrm{\Π}}}{a}_{i^{\′}{j}^{\′}}$

计算M_i'的5次方根V_i'

$V_{i^{\′}}=\sqrt[5]{M_{i^{\′}}}$

对V_i'进行归一化，即：

$F_{i^{\′}}=V_{i^{\′}}/{\mathrm{\Σ}}_{i^{\′}=1}^5{V}_{i^{\′}}$

则特征向量F＝(F₁,F₂,…,F₅)^T；

计算准则层判断矩阵A的最大特征根λ_max

${\mathrm{\λ}}_{\max }=\underset{i^{\′}=1}{\overset{5}{\mathrm{\Σ}}}\frac{{(A\•F)}_{i^{\′}}}{{5F}_{i^{\′}}}$

验证准则层判断矩阵A的一致性：

计算准则层判断矩阵A的一致性比率C_R，

其中，C_I为一致性检验指标，R_I＝1.12

当C_R＜0.1时，该准则层判断矩阵A为一致矩阵，否则，返回步骤(3.1)调整a_i'j'的大小，使准则层判断矩阵A最终满足一致矩阵；

(3.2)、将步骤(2)得到的评估值全部输入到综合评估模块，构造方案层判断矩阵B_j

将导航滤波算法EKF、UKF、UPF在每一个相同参数指标下构造一个3×3判断矩阵B_j，j＝1,2,…,5表示参数指标的个数，且其中,判断矩阵B_j是第j个参数指标对应三种导航滤波算法映射出来的3×3的矩阵，表示在该参数指标下映射的第个导航滤波算法的评估值同映射的第个导航滤波算法的评估值的比值，越大，在该参数指标下，映射的第个导航滤波算法的效果越好；

计算方案层判断矩阵B₁的特征值和最大特征根：

计算方案层判断矩阵B₁每一行比值的乘积

$N_{\tilde{i}}=\underset{\tilde{j}}{\overset{3}{\mathrm{\Π}}}{b}_{\tilde{i}\tilde{j}}$

计算的3次方根

$H_{\tilde{i}}=\sqrt[3]{N_{\tilde{i}}}$

对进行归一化处理，即：

$S_{\tilde{i}1}=H_{\tilde{i}}/{\mathrm{\Σ}}_{\tilde{i}=1}^3{H}_{\tilde{i}}$

则特征向量SS₁＝(S₁₁,S₂₁,S₃₁)^T；

计算判断矩阵的最大特征根λ_max1:

${\mathrm{\λ}}_{\max 1}=\underset{\tilde{i}=1}{\overset{3}{\mathrm{\Σ}}}\frac{{(B_1\•{\mathrm{SS}}_1)}_{\tilde{i}}}{{3S}_{\widetilde{i1}}}$

验证方案层判断矩阵B₁的一致性：

计算方案层判断矩阵B₁的一致性比率C_R1，

其中，C_I1为一致性检验指标，R_I1＝0.58

当C_R1＜0.1时，该方案层判断矩阵B₁为一致矩阵，否则，返回步骤(3.2)调整的大小，使准则层判断矩阵B₁最终满足一致矩阵；

同理，可以得到其他的方案层判断矩阵B₂～B₅的特征向量SS₂～SS₅和最大特征根λ_max2～λ_max5；

用方案层判断矩阵B_j的所有特征向量SS₁～SS₅构造矩阵C＝(SS₁,SS₂,…,SS₅)。

(3.3)、深空探测自主导航滤波算法的综合评估

将步骤(3.2)求得的特征向量C与步骤(3.1)构造的矩阵F相乘，求出权值W，即：W＝C·F，分别计算出权值W中每一行对应的权值其中，表示滤波算法的个数，是矩阵C中第行、第j列个元素；对所有权值进行总排序，其权值最大项即为对应的导航滤波方法最好。

2.根据权利要求1所述的一种对深空探测自主导航系统的评估方法，其特征在于，还包括：利用损失函数计算探测器在第i个采样点的位置差值ep(i)'和速度误差ev(i)'；

计算探测器在第i个采样点的位置差值ep(i)'：

将ex(i)与预设的位置判断值beta_p比较，如果ex(i)＜beta_p，则ex(i)保持不变，即ex(i)'＝ex(i)；如果ex(i)≥beta_p，则ex(i)'＝2beta_p·|ex(i)|-beta_p²；

同理再分别将ey(i)、ez(i)与预设的位置判断值beta_p比较，得到最终的ey(i)'、ez(i)'；

则探测器在第i个采样点的位置差值ep(i)'为：

$\mathrm{ep}{\left(i\right)}^{\′}=\sqrt{{\left[\mathrm{ex}{\left(i\right)}^{\′}\right]}^2+{\left[\mathrm{ey}{\left(i\right)}^{\′}\right]}^2+{\left[\mathrm{ez}{\left(i\right)}^{\′}\right]}^2}$

计算探测器在第i个采样点的速度误差ev(i)'：

将evx(i)与预设的位置判断值beta_v比较，如果evx(i)＜beta_v，则evx(i)保持不变，即evx(i)'＝evx(i)；如果evx(i)≥beta_v，则evx(i)'＝2beta_v·|evx(i)|-beta_v²；

同理再分别将evy(i)、evz(i)与预设的位置判断值beta_v比较，得到最终的evy(i)'、evz(i)'；

则探测器在第i个采样点的速度误差ev(i)'为：

$\mathrm{ev}{\left(i\right)}^{\′}=\sqrt{{\left[\mathrm{evx}{\left(i\right)}^{\′}\right]}^2+{\left[\mathrm{evy}{\left(i\right)}^{\′}\right]}^2+{\left[\mathrm{evz}{\left(i\right)}^{\′}\right]}^2}.$