[发明专利]一种应用星间跟踪的航天器星座分散化自主导航方法

摘要

一种应用星间跟踪的航天器星座分散化自主导航方法，它有如下步骤：一、各子滤波器初始化；二、各子滤波器进行本地状态采样；三、各子滤波器进行时间更新；四、各航天器之间建立星间通信链路并保持跟踪；五、已建立星间链路的航天器进行星间跟踪观测；六、经星间链路共享各子滤波器的状态采样信息；七、各子滤波器进行本地相关量测采样；八、各子滤波器进行量测更新；九、各子滤波器进行性能监控，判断其运行是否正常；十、各子滤波器将步骤八的量测更新结果作为本地导航估计输出，返回步骤一，开始执行下一个计算周期；十一、各子滤波器将步骤三的时间更新结果作为本地导航估计输出，返回步骤一，开始执行下一个计算周期。

主权项

一种应用星间跟踪的航天器星座分散化自主导航方法，其特征在于：该方法具体步骤如下：步骤1：各子滤波器初始化；步骤2：各子滤波器进行本地状态采样；步骤3：各子滤波器进行时间更新；步骤4：各航天器之间建立星间通信链路并保持跟踪；对于建立链路且相互跟踪成功航天器，进入步骤5；对于未成功建立任何链路的航天器，进入步骤11；步骤5：已建立星间链路的航天器进行星间跟踪观测；对于成功进行星间观测的航天器，根据可用的星间观测确定本地相关观测模型，进入步骤6；对于未成功进行星间观测的航天器，执行步骤11；步骤6：经星间链路共享各子滤波器的状态采样信息；步骤7：各子滤波器进行本地相关量测采样；步骤8：各子滤波器进行量测更新；步骤9：各子滤波器进行性能监控，判断其运行是否正常；若判断结果为正常，则执行步骤10，否则执行步骤11；步骤10：各子滤波器将步骤8的量测更新结果作为本地导航估计输出，返回步骤1，开始执行下一个计算周期；步骤11：各子滤波器将步骤3的时间更新结果作为本地导航估计输出，返回步骤1，开始执行下一个计算周期；其中，步骤1中所述的各子滤波器初始化，其实现方法为：各子滤波器初始化是指确定各子滤波器在当前计算时刻t_k的本地系统状态估计初值及相应的误差协方差矩阵初值对于整体算法的起始时刻，即t₀时刻，各子滤波器系统状态估计初值包括相应本地航天器在惯性参照坐标系中的位置矢量估计初值和速度矢量估计初值${\hat{X}}_0^{+}={\left[\begin{array}{c}{\hat{r}}_X\\ {\hat{v}}_X\end{array}\right]}_0---\left(1\right)$设t₀时刻本地航天器的系统状态真实值为X₀，则相应的误差协方差矩阵初值按照下式计算：$P_{XX,0}^{+}=E\{\[{\hat{X}}_0^{+}-X_0\]{\[{\hat{X}}_0^{+}-X_0\]}^T\}---\left(2\right)$若缺乏系统状态真实值X₀的必要信息，根据工程经验确定对于t_k(k＝1,2,...)时刻，和则等于上一步计算时刻的估计输出；其中，步骤2中所述的各子滤波器进行本地状态采样，其实现方法为：各子滤波器依据t_k时刻本地状态估计初值及相应的误差协方差矩阵初值并行地使用下面的对称采样算法计算相应的本地状态采样$\left\{\begin{array}{c}{\hat{X}}_k^{\left(0\right)}={\hat{X}}_k^{+}\\ {\hat{X}}_k^{\left(i\right)}={\hat{X}}_k^{+}+\sqrt{n+\mathrm{\τ}}{\left(\sqrt{P_{XX,k}^{+}}\right)}_i^T\\ {\hat{X}}_k^{(n+i)}={\hat{X}}_k^{+}-\sqrt{n+\mathrm{\τ}}{\left(\sqrt{P_{XX,k}^{+}}\right)}_i^T\end{array}\right.,(i=1,...,n)---\left(3\right)$其中，采样矢量共2n+1个，带括号的上标表示采样矢量序号；n为系统状态维数；τ为状态采样系数；当系统状态误差满足高斯分布时，选取n+τ＝3；其中，步骤3中所述的各子滤波器进行时间更新，其实现方法为：首先定义子滤波器状态动力学模型f_X(·)，与公式(1)对应，导航系统状态矢量X包含航天器在相应中心天体惯性系下的位置矢量r以及速度矢量v，导航系统状态动力学模型f_X(·)为：$\stackrel{\·}{X}=\left[\begin{array}{c}\stackrel{\·}{r}\\ \stackrel{\·}{v}\end{array}\right]=f_X(X,t)+w=\left[\begin{array}{c}v\\ a_{cen}+a_{ns}+a_{bd}+a_{srp}\end{array}\right]+w---\left(4\right)$其中，航天器受到中心天体质点引力加速度a_cen、中心天体非球形摄动加速度a_ns、太阳系主要天体质点引力加速度a_bg、太阳光压摄动加速度a_srp以及加速度模型误差w影响；根据航天器轨道动力学理论完成各引力项的计算，模型误差建模为零均值高斯白噪声；根据子滤波器状态动力学模型f_X(·)，建立相应的离散化状态模型F_X(·)：$X_{k+1}=F_X(X_k,t_k)={\∫}_{t_k}^{t_{k+1}}{f}_X(X_k,t_k)dt---\left(5\right)$接下来各子滤波器使用离散化状态模型F_X(·)并行地进行本地状态采样的时间更新，得到t_k+1时刻各自的本地状态一步预测及本地状态误差协方差矩阵一步预测计算公式为：$\left\{\begin{array}{c}{\hat{X}}_{k+1}^{-}=\underset{j=0}{\overset{2n}{\Σ}}{W}^{\left(j\right)}{F}_X({\hat{X}}_k^{\left(j\right)},t_k)\\ P_{XX,k+1}^{-}=\underset{j=0}{\overset{2n}{\Σ}}{W}^{\left(j\right)}({\hat{X}}_{k+1}^{\left(j\right)}-{\hat{X}}_{k+1}^{-}){({\hat{X}}_{k+1}^{\left(j\right)}-{\hat{X}}_{k+1}^{-})}^T+Q_k\end{array}\right.---\left(6\right)$其中j＝0,...,2n；Q_k为系统状态模型噪声对应的协方差阵；W^(j)为状态采样权值，计算公式为：$\left\{\begin{array}{c}{W}^{\left(0\right)}=\frac{\mathrm{\τ}}{n+\mathrm{\τ}}\\ W^{\left(i\right)}=\frac{1}{2(n+\mathrm{\τ})}\\ W^{(n+i)}=W^{\left(i\right)}\end{array}\right.,(i=1,...,n);---\left(7\right)$其中，步骤4中所述的各航天器之间建立星间通信链路并保持跟踪，其实现方法为：星间通信链路的建立和保持跟踪，通过星载空间通信及其链路捕获、跟踪、瞄准系统来完成；首先利用各航天器星载通信发射端机产生星间通信信号，通过天线向满足可视条件的其它航天器发射；后者使用天线和接收端机对星间通信信号进行捕获和确认，然后返回信标到发射端，从而完成初步的链路锁定，建立通信链路；接下来发射端航天器根据目标航天器的估计方位，驱动天线ATP伺服机构完成粗跟踪指向，然后提取通信信号的测角信息，导入信号发射方向微调反馈控制回路，保持通信链路稳定精确指向；其中，步骤5中所述的已建立星间链路的航天器进行星间跟踪观测，其实现方法为：首先定义航天器星间跟踪观测量包括航天器间的相对距离，相对速度和导航计算坐标系中的相对方位；以航天器A对航天器B的测量为例，假设其在惯性参考坐标系，记为i系中的位置矢量分别为和速度矢量分别为和相对视线矢量即相对位置矢量为相对速度矢量为相对距离为ρ_AB，相对速度为相对方位单位矢量为采用伪距式载波相位进行星间相对距离测量，利用无线电信号在空间定速传播的特性，测量其发射时刻与接收时刻的时间差来确定相对距离：ρ_AB＝cΔt_AB    (8)式中，c为电磁波传播速度，即光速；Δt_AB是测量信号的传播时间，由测距设备测定；利用多普勒频移可测定相对速度，测量关系为${\stackrel{\·}{\mathrm{\ρ}}}_{AB}=c\·\frac{{\mathrm{\Φ}}^{2}cos\mathrm{\θ}\±\sqrt{1-{\mathrm{\Φ}}^2{\sin }^2\mathrm{\θ}}}{{\mathrm{\Φ}}^2{\cos }^2\mathrm{\θ}+1}---\left(9\right)$其中，Φ为目标天体辐射观测频率与实际频率之比；θ为星间视线方向与相对速度方向的夹角；设参考恒星视线方向为该数据由星历数据库给出，且在惯性空间中视线位置变化极小；以其作为方位参考进行航天器间方位观测，使用恒星敏感器和星间跟踪测量设备，可测量航天器间的相对视线方向与在惯性空间中的相对角偏差接下来按照下面的公式精确获得航天器间的相对方位在相应中心天体惯性系中的单位矢量；$n_{AB}^i=n_S^i+{\mathrm{\Δn}}_{S,AB}^i---\left(10\right)$综合式(8)～(10)，以航天器A对航天器B的测量为例，星间测量值包括：$Z_{AB}=\left[\begin{array}{c}{\mathrm{\ρ}}_{AB}\\ {\stackrel{\·}{\mathrm{\ρ}}}_{AB}\\ n_{AB}^i\end{array}\right]---\left(11\right)$每一对可进行星间测量的航天器间均得到一组星间测量值；对于某一个航天器，所有与之相关的星间测量值组成其本地相关量测矢量Z_r,k+1；其中，步骤6中所述的经星间链路共享各子滤波器的状态采样信息，其实现方法为：经由星间链路，在各个测量相关的航天器间共享相应各子滤波器在步骤2中产生的状态采样信息；对于每个子滤波器，在将本地状态采样上传至星间链路的同时，获得来自所有与其存在星间测量的子滤波器的外部状态采样信息其中，步骤7中所述的各子滤波器进行本地相关量测采样，其实现方法为：首先定义观测模型，对于某个航天器对应的子滤波器，定义本地相关观测模型h_r(·)包括该航天器和所有与其存在星间测量链路的航天器间的相对距离观测模型、相对速度观测模型和相对方位观测模型；根据步骤5中的变量定义，以航天器A和航天器B为例，每一个星间观测量都至少同时与两个航天器的状态相关，相对距离观测模型为：${\widetilde{\mathrm{\ρ}}}_{AB}=\left|\right|r_A^i-r_B^i\left|\right|+{\mathrm{\ϵ}}_{\mathrm{\ρ},AB}---\left(12\right)$其中“～”标记表示相应变量的测量值，ε_ρ,AB表示航天器A和航天器B的相对距离测量误差，包括测量时延、钟差以及随机误差；相对速度观测模型表示为速度矢量差在位置矢量差方向上的投影：${\widetilde{\stackrel{\·}{\mathrm{\ρ}}}}_{AB}=(v_A^i-v_B^i)\·\frac{r_A^i-r_B^i}{\left|\right|r_A^i-r_B^i\left|\right|}+{\mathrm{\ϵ}}_{\stackrel{\·}{\mathrm{\ρ}},AB}---\left(13\right)$其中，表示航天器A和航天器B的相对速度测量误差；相对方位观测模型则为：${\tilde{n}}_{AB}^m=C_i^m\frac{r_A^i-r_B^i}{\left|\right|r_A^i-r_B^i\left|\right|}+{\mathrm{\ϵ}}_{n,AB}---\left(14\right)$其中，表示i系向星间测量坐标系即m系的姿态转换矩阵，由星载姿态确定系统测定；ε_n,AB表示航天器A和航天器B的相对方位测量误差；记为航天器A和航天器B的星间观测矢量，式(12)～式(14)构成了一组星间观测模型：${\tilde{Z}}_{AB}=h_{AB}(X_A,X_B)+{\mathrm{\ϵ}}_{AB}\stackrel{\mathrm{\Δ}}{=}\left[\begin{array}{c}\left|\right|r_A^i-r_B^i\left|\right|\\ (v_A^i-v_B^i)\·\frac{r_A^i-r_B^i}{\left|\right|r_A^i-r_B^i\left|\right|}\\ C_i^m\frac{r_A^i-r_B^i}{\left|\right|r_A^i-r_B^i\left|\right|}\end{array}\right]+\left[\begin{array}{c}{\mathrm{\ϵ}}_{\mathrm{\ρ},AB}\\ {\mathrm{\ϵ}}_{\stackrel{\·}{\mathrm{\ρ}},AB}\\ {\mathrm{\ϵ}}_{n,AB}\end{array}\right]---\left(15\right)$对于某个航天器，完整的本地相关量测模型h_r(·)包括该航天器和所有与之存在星间测量链路的航天器间的星间观测模型；接下来，各子滤波器并行地使用各自的h_r(·)计算相应的本地相关量测采样矢量${\hat{Z}}_{r,k+1}^{-}=\underset{j=0}{\overset{2n}{\Σ}}{W}^{\left(j\right)}{h}_r({\hat{X}}_k^{\left(j\right)},{\hat{Y}}_{r,k}^{\left(j\right)}),(j=0,...,2n);---\left(16\right)$其中，步骤8中所述的各子滤波器进行量测更新，其实现方法为：各子滤波器首先并行地计算相应的本地状态量测协方差矩阵P_XZr,k+1和本地量测协方差矩阵P_ZrZr,k+1：$\left\{\begin{array}{c}{P}_{XZr,k+1}=\underset{j=0}{\overset{2n}{\Σ}}{W}^{\left(j\right)}({\hat{X}}_{k+1}^{\left(j\right)}-{\hat{X}}_{k+1}^{-}){(Z_{r,k+1}-{\hat{Z}}_{r,k+1}^{-})}^T\\ P_{ZrZr,k+1}=\underset{j=0}{\overset{2n}{\Σ}}{W}^{\left(j\right)}({\hat{Z}}_{r,k+1}^{\left(j\right)}-{\hat{Z}}_{r,k+1}^{-}){({\hat{Z}}_{r,k+1}^{\left(j\right)}-{\hat{Z}}_{r,k+1}^{-})}^T\end{array}\right.---\left(17\right)$进而计算相应的增益矩阵K_k+1：$K_{k+1}=P_{XZr,k+1}{P}_{Zr,k+1}^{-1}---\left(18\right)$然后并行地计算t_k+1时刻各子滤波器相应的本地状态估计和本地状态估计误差协方差阵${\hat{X}}_{k+1}^{+}={\hat{X}}_{k+1}^{-}+K_{k+1}(Z_{r,k+1}-{\hat{Z}}_{r,k+1}^{-})---\left(19\right)$$P_{XX,k+1}^{+}=P_{XX,k+1}^{-}-K_{k+1}{P}_{Zr,k+1}{K}_{k+1}^T;---\left(20\right)$其中，步骤9中所述的各子滤波器进行性能监控，判断滤波器运行是否正常，其实现方法为：针对成员航天器可能出现测量或计算失效而造成算法故障的情况，承袭各成员航天器单独估计自身状态的独立估计方式，每个子滤波器独立检测自身故障；故障检测算法采用基于新息的经验卡方分布分析法，方法步骤如下：首先通过下面的表达式计算t_k+1时刻的新息ε_k+1：${\mathrm{\ϵ}}_{k+1}=Z_{r,k+1}-{\hat{Z}}_{r,k+1}^{-}---\left(21\right)$然后定义下面的等价统计函数：$\mathrm{\γ}=\frac{{\mathrm{\ϵ}}_{k+1}^T{P}_{XZr,k+1}^{-}{\mathrm{\ϵ}}_{k+1}}{l}---\left(22\right)$式中，l为量测量的维数，统计量γ是最小值为零的非负数；理论上，如果滤波器模型准确，且滤波未出现发散现象，γ将是一个标准的卡方分布；本算法中，设定γ的一个上限阈值γ_max作为滤波器发散的判据，当γ≤γ_max，则认为滤波器运行较好，且γ越小滤波性能越好；当γ＞γ_max，则认为滤波器出现故障；阈值的取值需要通过监视运行中的系统来确定，工程上通过仿真实验按经验和需求取定上限阈值γ_max。