CN103870714A - 一种基于高阶可观阵的导航系统可观测分析方法 - Google Patents

Abstract

一种基于高阶可观阵的导航系统可观测分析方法，包括根据导航系统的状态转移模型和观测模型，求取状态转移矩阵和观测矩阵；确定导航系统完全可观的最小历元数，构造低阶可观测矩阵，并利用其判断导航系统是否具备完全可观测性，若导航系统具备完全可观测性，再构造了高阶可观测矩阵，并求取不同阶数下的可观测矩阵条件数；最后利用这一系列条件数构成的条件数矢量评价导航系统可观测度。与现有可观测分析方法相比，本发明充分利用了前期导航观测量，提供了更加准确的可观测度信息，并且计算方法简单，便于工程实现。

Description

一种基于高阶可观阵的导航系统可观测分析方法

技术领域

本发明属于航天器自主导航领域，特别涉及一种导航系统可观测分析方法。

背景技术

系统的可观测分析对整个自主导航系统的提高具有重要意义。而系统可观测分析可分为可观测性和可观测度分析。其中，可观测性决定了该导航系统是否能对航天器进行定位，而可观测度则决定了系统对不同导航信息的有效确定程度，也即定位性能“好”与“坏”的问题。对于航天器自主导航系统这种非线性时变系统而言，其可观测分析对精度提升起着至关重要的作用。

目前，分段线性定常系统(Piece-Wise Constant Systems,PWCS)分析法是一种常用的导航系统可观测性分析方法。该方法假设非线性系统的系数矩阵变化量可以忽略不计，并利用条件数或奇异值分解(Singular Value Decomposition,SVD)获得的奇异值分析和评价可观测度。基于微分几何的分析方法利用李导数求解非线性系统的可观测矩阵。基于混合条件数的PWCS分析方法将真实的和近似的可观测阵的特征值条件数作为可观测度的评价准则。这些方法在构造可观性阵时将非线性项忽略或视为干扰。但是，在某些情况下，导航系统非线性有助于提高导航系统可观测度，即强非线性导航系统定位精度高于弱非线性导航系统定位精度。

发明内容

本发明提出了一种导航系统可观测分析方法，旨在为导航系统提供更加准确的可观测度分析结果。

本发明技术方案提供一种基于高阶可观阵的导航系统可观测分析方法，包括以下步骤，

步骤1，根据导航系统的状态转移模型和观测模型，按照式(1)和式(2)分别求取状态转移矩阵F和观测矩阵H如下，

$F=\frac{\∂f(X\left(t\right),t)}{\∂X\left(t\right)}---\left(1\right)$

$H=\frac{\∂h(X\left(t\right),t)}{\∂X\left(t\right)}---\left(2\right)$

其中，f(X(t),t)和h(X(t),t)分别为状态转移模型和观测模型，X(t)为时刻t的状态矢量；

步骤2，根据以下条件确定导航系统完全可观的最小历元数N，N≥1，

N＝argmin(N)s·t rank(O_L)＝6    (3)

$O_L=\underset{m=0}{\overset{N-1}{\mathrm{\Σ}}}{\left({\mathrm{HF}}^m\right)}^T\left({\mathrm{HF}}^m\right)---\left(4\right)$

其中，O_L为低阶可观测矩阵，T表示矢量转置，rank(·)表示矩阵的秩，变量m的取值为0,1,…,N-1；

如果N不存在，则导航系统是不可观测的，结束流程；

反之，该导航系统具备完全可观测性，继续步骤3；

步骤3，构造高阶可观性矩阵O_H(D_n)如下，

其中，D_n≥0，6D_n+6为高阶可观测矩阵阶数，m¹的取值为0,1,…,N-1，m²的取值为N,N+1,…,2N-1，

的取值为(D_n-1)N,(D_n-1)N+1,…,D_nN-1；

步骤4，求取条件数矢量C＝[C₁,C₂,…,C_J]，其中，C_j＝Cond(O_M(j))，j的取值为1,2,…,J，Cond(·)表示矩阵的条件数；通过比较条件数矢量C得到导航系统可观测分析结果。

而且，步骤4中，通过比较条件数矢量C得到导航系统可观测分析结果的实现方式为，对两个导航系统，依次比较C₁,C₂,…,C_J，当出现两个导航系统的某条件数C_j不等时，判断具有较小C_j的导航系统的可观测度更好。

本发明与现有技术相比的优点在于：

(1)与现有三种导航系统可观测分析方法相比，本发明充分利用了前期导航观测量信息，提供了更加准确的可观测度。

(2)与基于微分几何的导航系统可观测分析方法相比，本发明无需求解高阶李导数，便于工程计算，节约了系统资源。

附图说明

图1为本发明实施例的流程示意图。

具体实施方式

本发明技术方案可采用计算机软件方式支持自动运行流程。以下结合附图和实施例详细说明本发明技术方案。

实施例采用X射线脉冲星导航系统为例对本发明作进一步说明。X射线脉冲星导航系统中的状态转移模型和测量模型分别为航天器轨道动力学模型和脉冲到达时间转换模型。为便于实施参考起见，本发明先对这两个模型进行说明。

航天器轨道动力学模型表示为

其中，f(X(t),t)为状态转移模型，状态矢量X＝[x y z v_x v_y v_z]^T。x,y,z,v_x,v_y,v_z分别是航天器在三个方向上的位置和速度，X(t)为时刻t的状态矢量，

为X(t)的导数。w(t)＝[0,0,0,w_x,w_y,w_z]为状态处理噪声，w_x，w_y，w_z分别为航天器在三个方向上的加速度噪声，T表示矢量转置。该模型具体表达式如下：

$\left\{\begin{array}{c}\frac{\mathrm{dx}}{\mathrm{dt}}=v_x\\ \frac{\mathrm{dy}}{\mathrm{dt}}=v_y\\ \frac{\mathrm{dz}}{\mathrm{dt}}=v_z\\ \frac{{\mathrm{dv}}_x}{\mathrm{dt}}=-\mathrm{\μ}\frac{x}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-1.5)]+\mathrm{\Δ}{F}_x+w_x\\ \frac{{\mathrm{dv}}_y}{\mathrm{dt}}=-\mathrm{\μ}\frac{y}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-1.5)]+\mathrm{\Δ}{F}_y+w_y\\ \frac{{\mathrm{dv}}_z}{\mathrm{dt}}=-\mathrm{\μ}\frac{z}{r^3}[1-J_2{\left(\frac{R_e}{r}\right)}^2(7.5\frac{z^2}{r^2}-4.5)]+\mathrm{\Δ}{F}_z+w_z\end{array}\right.---\left(6\right)$

其中，μ是地球引力常数，

为航天器与地球质心之间的距离，ΔF_x、ΔF_y、ΔF_z为地球非球形摄动的高阶摄动项，日、月摄动，以及太阳光压摄动和大气摄动等影响航天器位置的摄动力，J₂为二阶带谐项系数，R_e为地球半径。

X射线脉冲到达时间转换模型为：

$c(t_b^i-t_{\mathrm{SC}}^i)=n^i{r}_{\mathrm{SC}}+\frac{1}{2D_0^i}[-{\left|r_{\mathrm{SC}}\right|}^2+{\left(n^i{r}_{\mathrm{SC}}\right)}^2-2{\mathrm{br}}_{\mathrm{SC}}+2\left(n^{i}b\right)\left(n^i{r}_{\mathrm{SC}}\right)]+\frac{2{\mathrm{\μ}}_{\mathrm{Sun}}}{c^2}\ln |\frac{n^i{r}_{\mathrm{SC}}+\left|r_{\mathrm{SC}}\right|}{n^{i}b+\left|b\right|}+1|---\left(7\right)$

$n^i=\left[\begin{array}{c}\cos \left({\mathrm{\δ}}^i\right)\cos \left({\mathrm{\α}}^i\right)\\ \cos \left({\mathrm{\δ}}^i\right)\sin \left({\mathrm{\α}}^i\right)\\ \sin \left({\mathrm{\δ}}^i\right)\end{array}\right]---\left(8\right)$

其中，

为航天器与太阳系质心之间的距离在第i颗脉冲星方向上的投影，nⁱ是第i颗脉冲星的方向矢量，i＝1,2,…I，I为导航的脉冲星数量；αⁱ和δⁱ分别为第i颗脉冲星的赤经和赤纬，

和

分别为第i颗脉冲星的脉冲到达航天器时间与到达太阳系质心时间，c为光速，

为第i颗脉冲星到太阳系质心的距离；b为太阳系质心相对于太阳的位置矢量，|b|为位置矢量b的长度；μ_Sun为太阳引力常数。r_SC是航天器相对于太阳系质心的位置矢量，|r_SC|为位置矢量r_SC的长度。利用标准星历表提供的地球位置r_E，可将r_SC转化为航天器相对于地球的位置矢量r。

r＝r_SC-r_E    (9)

假设X射线脉冲星导航观测量Z为：

$Z=\left[\begin{array}{c}c(t_b^1-t_{\mathrm{SC}}^1)\\ c(t_b^2-t_{\mathrm{SC}}^2)\\ .\\ .\\ .\\ c(t_b^I-t_{\mathrm{SC}}^I)\end{array}\right]---\left(10\right)$

其对应的量测噪声为V。X射线脉冲星导航观测模型可表示为：

Z＝h(X(t),t)+V(t)    (11)

其中，V(t)为时刻t的量测噪声，观测模型h(X(t),t)如下，

$h(X\left(t\right),t)=\left[\begin{array}{c}{h}_1(X\left(t\right),t)\\ h_2(X\left(t\right),t)\\ .\\ .\\ .\\ h_i(X\left(t\right),t)\\ .\\ .\\ .\\ h_I(X\left(t\right),t)\end{array}\right]---\left(12\right)$

其中，第i颗脉冲星的相应项h_i(X(t),t)的表达式如下：

$h_i(X\left(t\right),t)=n^i{r}_{\mathrm{SC}}+\frac{1}{2D_0^i}[-{\left|r_{\mathrm{SC}}\right|}^2+{\left(n^i{r}_{\mathrm{SC}}\right)}^2-2{\mathrm{br}}_{\mathrm{SC}}+2\left(n^{i}b\right)\left(n^i{r}_{\mathrm{SC}}\right)]+\frac{2{\mathrm{\μ}}_{\mathrm{Sun}}}{c^2}\ln |\frac{n^i{r}_{\mathrm{SC}}+\left|r_{\mathrm{SC}}\right|}{n^{i}b+\left|b\right|}+1|---\left(13\right)$

本发明实施例采用的导航脉冲星为B0531+21，B1821-24，B1937+21。其方位参数如表1所示。

表1脉冲星方位参数

参见图1，下面给出本发明实施例的具体实施步骤：

步骤1：状态转移矩阵F和观测矩阵H。根据导航系统的状态转移模型和观测模型，按照式(14)和式(15)分别求取状态转移矩阵和观测矩阵。

$F=\frac{\∂f(X\left(t\right),t)}{\∂X\left(t\right)}---\left(14\right)$

$H=\frac{\∂h(X\left(t\right),t)}{\∂X\left(t\right)}---\left(15\right)$

其中，f(X(t),t)和h(X(t),t)分别为状态转移模型和观测模型，X(t)为时刻t的状态矢量。

实施例求取观测矩阵H和状态转移矩阵F如下，

$F=\frac{\∂f(X\left(t\right),t)}{\∂X\left(t\right)}=\left[\begin{array}{cc}{0}_{3\×3}& I_{3\×3}\\ S_{3\×3}& 0_{3\×3}\end{array}\right]---\left(16\right)$

其中，0_3×3和I_3×3分别为3×3的零矩阵和单位矩阵，忽略地球非球形摄动以及其他摄动力的影响，S_3×3可近似地表示为：

$S_{3\×3}\≈\left[\begin{array}{ccc}\frac{\mathrm{\μ}(3x^2-r^2)}{r^5}& \frac{3\mathrm{\μxy}}{r^5}& \frac{3\mathrm{\μxz}}{r^5}\\ \frac{3\mathrm{\μxy}}{r^5}& \frac{\mathrm{\μ}(3y^2-r^2)}{r^5}& \frac{3\mathrm{\μyz}}{r^5}\\ \frac{3\mathrm{\μxz}}{r^5}& \frac{3\mathrm{\μyz}}{r^5}& \frac{\mathrm{\μ}(3z^2-r^2)}{r^5}\end{array}\right]---\left(17\right)$

由于相对于一阶量，广义相对论高阶量的影响相对较小，求观测矩阵时该影响可以忽略不计。量测矩阵H可表示如下：

$H=\frac{\∂h(X\left(t\right),t)}{\∂X\left(t\right)}=\left[\begin{array}{cc}{n}^T& 0_{3\×3}\end{array}\right]---\left(18\right)$

其中，脉冲星方向矩阵n＝[n¹,n²,…,n^I]。

步骤2：确定参数N的值，即使导航系统完全可观的最小历元数N。对于航天器自主导航系统而言，该N(N≥1)应满足以下条件：

N＝argmin(N)s·t rank(O_L)＝6    (19)

$O_L=\underset{m=0}{\overset{N-1}{\mathrm{\Σ}}}{\left({\mathrm{HF}}^m\right)}^T\left({\mathrm{HF}}^m\right)---\left(20\right)$

其中，O_L为低阶可观测矩阵，T表示矢量转置。rank(·)表示矩阵的秩，变量m的取值为0,1,…,N-1。

如果N不存在，则导航系统是不可观测的，结束流程；

反之，该导航系统具备完全可观测性，继续步骤3。

实施例中N＝2。导航系统具备完全可观测性。进行第3步。

步骤3：构造高阶可观性矩阵O_H(D_n)。

其中，变量D_n≥0。6D_n+6为高阶可观测矩阵阶数。变量m¹的取值为0,1,…,N-1，变量m²的取值为N,N+1,…,2N-1，变量

的取值为(D_n-1)N,(D_n-1)N+1,…,D_nN-1。

实施例中，将N的值、H和F的表达式代入上式。可得到高阶可观性矩阵O_H(D_n)表达式为：

步骤4：求取条件数矢量C＝[C₁,C₂,…,C_J]。其中，C_j＝Cond(O_M(j))，变量j的取值为1,2,…,J，Cond(·)表示矩阵的条件数。条件数矢量C中的各条件数越小，导航系统可观测度就越好，即导航定位精度越高。

其中，O_M(j)就是用变量j替代O_H(D_n)中D_n的结果，C_j表示矩阵O_M(j)的条件数。6j+6为高阶可观测矩阵的阶数。6J+6为实施例中的最高阶数。

可对两个导航系统，依次比较C₁,C₂,…,C_J，当出现两个导航系统的某条件数C_j不等时，判断具有较小C_j的导航系统的可观测度更好。实施例将两个系统对比进行时，先比较C₁。具有较小C₁的系统可观测度更好，即精度更高。若两个系统的C₁相等，再比较C₂。依此类推。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (2)

1.一种基于高阶可观阵的导航系统可观测分析方法，其特征在于：包括以下步骤，

步骤1，根据导航系统的状态转移模型和观测模型，按照式(1)和式(2)分别求取状态转移矩阵F和观测矩阵H如下，

$F=\frac{\∂f(X\left(t\right),t)}{\∂X\left(t\right)}---\left(1\right)$

$H=\frac{\∂h(X\left(t\right),t)}{\∂X\left(t\right)}---\left(2\right)$

其中，f(X(t),t)和h(X(t),t)分别为状态转移模型和观测模型，X(t)为时刻t的状态矢量；

步骤2，根据以下条件确定导航系统完全可观的最小历元数N，N≥1，

N＝argmin(N)s·t rank(O_L)＝6    (3)

$O_L=\underset{m=0}{\overset{N-1}{\mathrm{\Σ}}}{\left({\mathrm{HF}}^m\right)}^T\left({\mathrm{HF}}^m\right)---\left(4\right)$

其中，O_L为低阶可观测矩阵，T表示矢量转置，rank(·)表示矩阵的秩，m的取值为0,1,…,N-1；

如果N不存在，则导航系统是不可观测的，结束流程；

反之，该导航系统具备完全可观测性，继续步骤3；

步骤3，构造高阶可观性矩阵O_H(D_n)如下，

其中，D_n≥0，6D_n+6为高阶可观测矩阵阶数，m¹的取值为0,1,…,N-1，m²的取值为N,N+1,…,2N-1，

的取值为(D_n-1)N,(D_n-1)N+1,…,D_nN-1；

步骤4，求取条件数矢量C＝[C₁,C₂,…,C_J]，其中，条件数C_j＝Cond(O_M(j))，j的取值为1,2,…,J，Cond(·)表示矩阵的条件数；通过比较条件数矢量C得到导航系统可观测分析结果。

2.根据权利要求1所述基于高阶可观阵的导航系统可观测分析方法，其特征在于：步骤4中，通过比较条件数矢量C得到导航系统可观测分析结果的实现方式为，对两个导航系统，依次比较C₁,C₂,…,C_J，当出现两个导航系统的某条件数C_j不等时，判断具有较小C_j的导航系统的可观测度更好。

Images