CN107421550B - 一种基于星间测距的地球-Lagrange联合星座自主定轨方法 - Google Patents

一种基于星间测距的地球-Lagrange联合星座自主定轨方法 Download PDF

Info

Publication number
CN107421550B
CN107421550B CN201710611461.6A CN201710611461A CN107421550B CN 107421550 B CN107421550 B CN 107421550B CN 201710611461 A CN201710611461 A CN 201710611461A CN 107421550 B CN107421550 B CN 107421550B
Authority
CN
China
Prior art keywords
earth
satellite
noise
lagrange
time
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201710611461.6A
Other languages
English (en)
Other versions
CN107421550A (zh
Inventor
杨静
卢帅
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beihang University
Original Assignee
Beihang University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beihang University filed Critical Beihang University
Priority to CN201710611461.6A priority Critical patent/CN107421550B/zh
Publication of CN107421550A publication Critical patent/CN107421550A/zh
Application granted granted Critical
Publication of CN107421550B publication Critical patent/CN107421550B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/24—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 specially adapted for cosmonautical navigation
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
    • G01C21/00—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00
    • G01C21/02—Navigation; Navigational instruments not provided for in groups G01C1/00 - G01C19/00 by astronomical means

Landscapes

  • Engineering & Computer Science (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Astronomy & Astrophysics (AREA)
  • Automation & Control Theory (AREA)
  • General Physics & Mathematics (AREA)
  • Other Investigation Or Analysis Of Materials By Electrical Means (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

本发明公开了一种基于星间测距的地球‑Lagrange联合星座自主定轨方法，包括以下几个步骤：步骤一：建立地球‑Lagrange联合星座自主定轨***的状态方程；步骤二：建立地球‑Lagrange联合星座自主定轨***的测量方程；步骤三：确定实现轨道参数估计的滤波方法；步骤四：基于选定滤波算法的地球‑Lagrange联合星座自主定轨方法的具体实现。本发明通过引入Lagrange卫星，有效解决了仅利用星间测距信息进行自主定轨时存在的“亏秩”问题，减轻了***设备的复杂性；通过自适应非线性滤波算法在线实时估计***噪声的统计特性，对噪声先验信息要求较低，提高了自主定轨滤波算法的稳定性，提高了自主定轨精度。

Description

一种基于星间测距的地球-Lagrange联合星座自主定轨方法
技术领域
本发明属于地球卫星星座的自主定轨领域，具体来说，是一种仅利用地球-拉格朗日(Lagrange)联合星座内的星间测距信息，基于自适应滤波实现地球卫星星座自主定轨的方法。
背景技术
近几十年来，卫星导航在国民经济和军事斗争领域中发挥着越来越重要的作用，提高卫星自主运行能力对于减轻地面测控负担、降低卫星运行费用、提高卫星生存能力和扩展航天器的应用潜力等方面具有重要意义。卫星的自主定轨技术是实现卫星自主控制的前提，对保障其自主运行发挥着至关重要的作用，也成为当今卫星导航与控制技术的发展趋势。
目前星座的自主定轨方法主要分为各卫星独立自主定轨和星座整体自主定轨两类方法。前者是通过卫星上搭载的各种惯性器件、星敏感器和导航接收机等获取卫星相对于其他自然天体或者导航星的距离、方向等信息来对自身位置进行在线估计，***构成简单，运算量小，但其定轨精度受到自然天体表面的不规则程度以及敏感器的测量精度等因素的限制，单独使用难以达到导航级别应用的精度；后者则是充分利用星座内部各卫星间的相对测量信息来实现星座的整网定轨，短弧段定轨精度高，但存在着“亏秩”问题，即仅利用星间相对测量信息无法确定星座在惯性空间中的整体性方向旋转，造成星座卫星的长弧段定轨精度随时间增长而逐步下降。为了解决“亏秩”问题，国外有学者指出当星座或卫星编队所处的引力场高度不对称时，便会由轨道的“唯一性”引入绝对信息，从而解决星座整体旋转的不可观问题。研究发现，在众多非对称摄动引力场中，第三体引力的相对强度最大，特别是在L1和L2两个Lagrange点附近区域，因此通过联合地球卫星星座和Lagrange点星座进行联合自主定轨可以解决“亏秩”问题。利用此种自主定轨方法可以大大降低航天器自主导航***的复杂性，减少导航测量设备，提高自主定轨精度。同时，Lagrange导航星座还可以作为深空探测导航基站，为其他深空探测器提供导航信息，从而有效解决地面深空探测网在深空导航精度、实时性、安全性等方面存在的问题，因此研究地球-Lagrange联合星座的自主定轨具有重大意义。
目前，国内外关于地球-Lagrange联合星座的自主定轨方法中主要以批处理算法、扩展卡尔曼滤波(Extend Kalman Filter，EKF)算法和无迹卡尔曼滤波(Unscented KalmanFilter，UKF)算法等常规算法为基础。但是，批处理算法往往用于事后处理，不适用于实时定轨；常规的EKF和UKF算法虽然采用了便于实时计算的递推形式，但往往要求***噪声和测量噪声的先验信息已知，否则将造成***性能下降甚至引起滤波发散。由于在轨卫星的实际受力情况复杂，特别是Lagrange卫星，难以获得较为精确的噪声统计特性，而且卫星轨道维持的机动等导致的不规律建模误差也会造成统计特性的改变，使得先验信息失去意义，因此研究能够根据测量信息序列自适应调整的自主导航算法对于提高自主定轨***的鲁棒性，保证自主定轨精度有着重要意义。
发明内容
本发明为了解决上述问题，提出一种基于星间测距的地球-Lagrange联合星座自主定轨方法，该方法仅利用星间链路获取的星-星距离测量信息，利用自适应非线性滤波(Adaptive Nonlinear Filter,ANF)算法，采用集中式结构融合有效测量信息，实现地球-Lagrange联合星座中各个卫星位置和速度的解算。该星座定轨方法首先采用集中式结构建立地球-Lagrange联合星座自主定轨的***模型，其中地球导航卫星的状态方程主要基于地心惯性直角坐标系下的受摄二体问题动力学模型来建立，Lagrange导航卫星的状态方程主要基于质心会合坐标系下的圆型限制性三体问题动力学模型来建立，测量方程采用地心惯性直角坐标系下的星间距离测量模型；然后，采用自适应非线性滤波算法对有效测量信息进行融合实现地球-Lagrange联合星座的自主定轨；最后以3颗中地球轨道(MEO)卫星和2颗Lagrange卫星组成的联合星座为例验证方法的有效性。
一种基于星间测距的地球-Lagrange联合星座自主定轨方法，具体包括以下几个步骤：
步骤一：在集中式结构下，建立地球-Lagrange联合星座自主定轨***的状态方程；
步骤二：在集中式结构下，建立地球-Lagrange联合星座自主定轨***的测量方程；
步骤三：在集中式结构下，根据上述建立的地球-Lagrange联合星座自主定轨***模型，确定实现轨道参数估计的滤波方法；
步骤四：在集中式结构下，基于选定滤波算法的地球-Lagrange联合星座自主定轨方法的具体实现。
本发明的优点在于：
(1)通过引入Lagrange卫星为***提供绝对信息，有效解决“亏秩”问题，提高了地球-Lagrange联合星座长期自主定轨精度；
(2)仅利用星间测距信息，减轻了***设备的复杂性；
(3)Lagrange导航星座对地球导航星座的覆盖率高，提高了自主定轨精度；
(4)利用自适应非线性滤波算法在线实时估计***噪声的统计特性，提高了自主定轨的滤波算法精度和稳定性，从而提高自主定轨精度。
附图说明
图1为基于星间测距的地球-Lagrange联合星座自主定轨***构成示意图；
图2为基于星间测距的地球-Lagrange联合星座自主定轨方法结构图；
图3为基于自适应非线性滤波的自主导航算法流程图；
图4为实施例***精度随bQ、bC、bq变化趋势图；
图5为实施例GNSS卫星的位置误差曲线图。
具体实施方式
下面将结合附图和实施例对本发明作进一步的详细说明。
本发明是一种基于星间测距的地球-Lagrange联合星座自主定轨方法，其***构成见图1，图中显示该联合星座定轨***由地球导航星座及其星间链路测距、Lagrange导航星座及其星间链路测距以及两星座卫星间的星间链路测距构成，其中地球导航星座由m颗中/高轨卫星组成，编号依次为1,2,…,i,…m；Lagrange导航星座由n颗运行在L1和L2点附近Halo周期轨道上的卫星组成，编号依次为1,2,…,k,…n。一种基于星间测距的地球-Lagrange联合星座自主定轨方法的基本实现过程参见图2，图中显示该定轨方法的基本实现过程主要包括测量信息获取和星座轨道参数估计。其中获取的测量信息包括地球导航星座内部的星间距离测量信息、Lagrange导航星座内部的星间距离测量信息以及地球导航卫星和Lagrange导航卫星间的距离测量信息三类距离测量信息；在进行星座轨道参数估计时，①基于地球卫星和Lagrange卫星的简化轨道动力学模型建立联合星座定轨***的状态方程，②建立三类距离测量信息的测量方程，③利用由状态方程和测量方程构成的***模型，采用合适的自适应非线性滤波算法，实现轨道参数的估计；最后，输出估计得到的星座中卫星的位置和速度。在地心惯性坐标系和质心会和坐标系下，采用直角坐标描述法(即用位置和速度分量来表示卫星当前的运行状态)，利用自适应非线性滤波估计方法，结合根据轨道动力学方程和测量信息建立的***模型，实现地球-Lagrange联合星座的自主定轨，具体步骤为：
步骤一：在集中式结构下，建立地球-Lagrange联合星座自主定轨***的状态方程；
准确的***模型是保障星座运行状态参数估计精度的主要因素之一。本发明基于卫星轨道动力学模型来建立***状态方程。
(1)地球导航卫星轨道动力学模型
地球导航卫星主要是中高轨卫星，受到地球的质心引力作用远大于其他作用力，一般采用受摄二体问题动力学模型。为了满足计算高效和高精度的要求，根据摄动力影响大小来选择主要摄动项建模。对于中、高轨卫星来说，在诸多摄动项中，地球非球形对称带来的摄动影响可达104米级，其中3阶以上引力场摄动影响与J2项相比小很多，可以忽略；日月引力摄动影响可达103米级；太阳光压摄动影响可达102米级；而潮汐、相对论效应和反照光压摄动等其他摄动的影响在10-1米级以下，可以忽略；同时高轨处大气密度甚微，大气阻力摄动影响可以忽略。因此，本发明主要考虑了地球质心引力以及J2项和日、月引力三种主要摄动力构成的轨道动力学模型来建立***状态方程。
在地心直角惯性坐标系下，结合轨道动力学模型建立对联合星座中待测地球卫星i进行定轨的状态方程如下：
Figure BDA0001359618900000041
式中，待测地球卫星i对应的状态向量为
Figure BDA0001359618900000042
这里[xEi,yEi,zEi]T和
Figure BDA0001359618900000043
分别为该地球卫星的位置矢量和速度矢量。FE为***状态函数，WEi(t)为***噪声，满足均值为qEi(t)，方差为QEi(t)的高斯白噪声。式(1)可以进一步详细写为：
Figure BDA0001359618900000051
其中，
Figure BDA0001359618900000052
和
Figure BDA0001359618900000053
为x、y、z三轴速度噪声分量，
Figure BDA0001359618900000054
和
Figure BDA0001359618900000055
是x、y、z三轴加速度等效噪声分量，包括未建模摄动项和噪声项；
Figure BDA0001359618900000056
和
Figure BDA0001359618900000057
分别为地球卫星i对应的已建模速度矢量和加速度矢量，其中加速度主要考虑地球中心引力加速度a0,Ei、J2摄动项加速度
Figure BDA0001359618900000058
和日月引力摄动加速度aMS,Ei，具体表达式如下：
1)本发明考虑的卫星i的地球中心引力加速度a0,Ei满足：
Figure BDA0001359618900000059
式中，
Figure BDA00013596189000000510
为待测卫星i的地心距，μ为地球质量ME与引力常数G的乘积，即地球引力常数。
2)J2摄动项加速度
Figure BDA00013596189000000511
满足：
Figure BDA00013596189000000512
式中，Re为地球赤道半径，J2为二阶带谐系数。
3)日月引力引起的摄动加速度aMS,Ei满足：
Figure BDA0001359618900000061
式中，(xS,yS,zS)和(xM,yM,zM)分别为太阳和月球在地心直角惯性坐标系下的三维位置坐标，
Figure BDA0001359618900000062
和
Figure BDA0001359618900000063
分别为卫星i到太阳和月球的距离，
Figure BDA0001359618900000064
和
Figure BDA0001359618900000065
分别为地心到太阳和月球的距离，μS和μM分别为太阳引力常数和月球引力常数。
(2)Lagrange导航卫星轨道动力学模型
与地球导航卫星不同，Lagrange卫星运行在地月系Lagrange点附近，受到的地球中心引力作用和月球中心引力作用量级相近，所以在描述其运动时，本发明采用圆型限制性三体模型。
在质心会合坐标系下，结合轨道动力学模型可建立对Lagrange星座中待测月球卫星k进行定轨的状态方程如下：
Figure BDA0001359618900000066
式中，待测Lagrange月球卫星k对应的状态向量为
Figure BDA0001359618900000067
这里[xLk,yLk,zLk]T和
Figure BDA0001359618900000068
分别为该星的无量纲位置矢量和速度矢量，其特征质量、特征长度和特征时间如式(7)所示，其中ME和MM分别为地球质量和月球质量，L为地月距离。FL为***状态函数，WLk(t)为***噪声，满足均值为qLk(t)，方差为QLk(t)的高斯白噪声。
Figure BDA0001359618900000069
式(6)可以进一步详细写为：
Figure BDA0001359618900000071
其中
Figure BDA0001359618900000072
为Lagrange卫星k对应的已建模速度矢量，μ0＝MM/(MM+ME)，ΔE,Lk、ΔM,Lk分别为Lagrange卫星k的地心距和月心距，
Figure BDA0001359618900000073
和
Figure BDA0001359618900000074
为Lagrange卫星k的x、y、z三轴速度噪声分量，
Figure BDA0001359618900000075
和
Figure BDA0001359618900000076
是Lagrange卫星k的x、y、z三轴加速度等效噪声分量，包括未建模摄动项和噪声项。
(3)地球-Lagrange联合星座自主定轨***的状态方程
假设整个自主定轨***包含m颗地球卫星和n颗Lagrange卫星。根据上述单颗卫星定轨的状态方程(1)和(6)，建立集中式结构的星座定轨状态方程如下：
Figure BDA0001359618900000077
即：
Figure BDA0001359618900000078
式中，
Figure BDA0001359618900000079
为状态向量，F为整个星座定轨***的状态函数矢量，W(t)表示均值为q(t)，方差为Q(t)的***噪声。
步骤二：在集中式结构下，建立地球-Lagrange联合星座自主定轨***的测量方程；
地球-Lagrange联合星座自主定轨***所用到的观测信息只有星间距离测量，但由于两类卫星的动力学模型是在不同坐标系中建立的，所以整个星座的星间距离测量模型被分为三类：
(1)地球导航卫星间的星间距离测量模型
对于地球导航卫星，记卫星i和卫星j之间的距离观测量为ρEi,Ej，则有：
ρEi,Ej＝hE(XEi,XEj)+vEi,Ej，i≠j (11)
其中
Figure BDA0001359618900000081
vEi,Ej为卫星i和卫星j之间的距离等效测量噪声，包括建模误差和量测噪声。
(2)Lagrange导航卫星间的星间距离测量模型
记Lagrange卫星i和j之间的距离观测量为ρLi,Lj，有：
ρLi,Lj＝hL(XLi,XLj)+vLi,Lj，i≠j (12)
其中
Figure BDA0001359618900000082
vLi,Lj为Lagrange卫星i和j之间的距离等效测量噪声，包括建模误差和量测噪声。
(3)地球导航卫星和Lagrange卫星间的测量模型
为了表示GNSS卫星i和Lagrange卫星k之间的距离关系，需要将两颗卫星的位置参数放在同一坐标系中表示，则有：
Figure BDA0001359618900000083
其中，
Figure BDA0001359618900000084
是Lagrange卫星在地心惯性坐标系中的位置向量，它是关于Lagrange卫星在质心会合坐标系中的状态向量XLk中的无量纲位置向量[xLk,yLk,zLk]T的函数，其转换关系为：
Figure BDA0001359618900000085
其中，R表示由质心会合坐标系到地心惯性坐标系的旋转矩阵，表达式如下：
Figure BDA0001359618900000091
其中uM＝ωM+θM，iM、ΩM、ωM和θM分别表示月球绕地轨道的轨道倾角、升交点赤经、近地点幅角以及真近点角。
(4)地球-Lagrange联合星座自主定轨***的测量方程
根据上述测量方程(11)、(12)和(13)，建立集中式结构的星座定轨测量方程如下：
Z(t)＝h[X(t),t]+V(t) (16)
即：
Figure BDA0001359618900000092
其中V(t)为零均值，方差为R(t)的等效测量噪声。式(17)给出的是一般情况下的观测方程，在实际应用时，某些距离测量量会受到地球或者月球的遮挡等因素的影响而出现不可测的情况，这时需要从式(17)中将相应的测距信息剔除。
步骤三：根据上述建立的基于星间距离测量的地球-Lagrange联合星座自主定轨***模型，确定实现轨道参数估计的滤波方法；
该非线性***的状态需要采用非线性滤波方法来估计。常规的非线性滤波方法，如EKF、UKF、容积滤波和粒子滤波等，在噪声的先验统计信息精确已知的情况下可以获得较为理想的结果。但是在本***中，等效***噪声Wk除了包含未建模误差之外，还包含滤波算法中的离散化误差(和线性化误差)，这使得噪声具有时变非高斯统计特性，且难以准确获取。为此，一般在应用这些非线性滤波算法时，通常将噪声近似为平稳随机白噪声序列，即认为Qk≡Q，并根据经验来确定Q。为了保证在整个滤波过程中状态估计误差能收敛到稳定值，Q的取值往往偏大，即满足Q≥max{Qk}，这使得***性能难以达到最优。
为了解决这一问题，本发明在常规非线性滤波算法的基础上，结合Sage-Husa自适应算法的思想对***噪声的统计特性进行在线估计和修正，以提高滤波的精度和稳定性。考虑到工程实用性，本发明采用自适应EKF方法或自适应Sigma点卡尔曼滤波(SigmaPoints Kalman Filter,SPKF)方法，实现卫星轨道参数估计。
第一，当***对运算速度要求较高或***存储空间较小时，采用自适应EKF方法对经过离散化后的非线性***进行状态估计。自适应EKF的基本思想是利用泰勒级数展开的方法对连续非线性的状态方程和观测方程进行离散化和线性化处理，然后在经典卡尔曼滤波的基础上，通过时变噪声统计估值器，实时估计和修正***噪声和观测噪声的统计特性，再进行状态估计。具体步骤如下：
(1)时间更新
计算一步预测状态：
Figure BDA0001359618900000101
式中，
Figure BDA0001359618900000102
为k时刻的***噪声均值的估计值，
Figure BDA0001359618900000103
为利用动力学方程进行的状态一步递推，一种方法是通过离散化和线性化后的方程(19)计算得到，其中
Figure BDA0001359618900000104
为k时刻的状态估计，
Figure BDA0001359618900000105
为(k+1)时刻的一步预测状态估计，Δt为积分步长；另一种方法是通过4阶Runge-Kutta(RK4)方法直接对连续状态方程(9)数值积分计算得到。
Figure BDA0001359618900000106
计算状态转移矩阵：
Figure BDA0001359618900000107
式中，Φk+1/k为从k时刻到(k+1)时刻的状态转移矩阵。
计算一步预测误差协方差矩阵：
Figure BDA0001359618900000108
式中，Pk+1/k为一步预测状态的误差协方差矩阵，Pk为k时刻估计状态的误差协方差矩阵，
Figure BDA0001359618900000109
为k时刻***噪声方差阵的估计值。
考虑到滤波初始的一段时间内***状态估计误差尚未收敛且用于估计***噪声的样本较少，此时估计的噪声统计特性波动较大且不够准确，将使滤波收敛较慢，严重时甚至导致滤波发散。因此在滤波初始的一段时间内(0～T)，本算法并不将估计的噪声统计特性应用于时间更新过程，而是采用噪声的先验统计信息，此时的时间更新算法如下：
Figure BDA0001359618900000111
Figure BDA0001359618900000112
(2)测量更新：
计算增益矩阵：
Figure BDA0001359618900000113
式中，Kk+1为(k+1)时刻的增益矩阵，
Figure BDA0001359618900000114
为基于一步预测状态估计
Figure BDA0001359618900000115
计算的观测向量h的Jacobian矩阵，Rk+1为(k+1)时刻测量噪声协方差阵。
计算新息向量：
Figure BDA0001359618900000116
式中υk+1为(k+1)时刻的新息向量，Zk+1为(k+1)时刻的测量向量。
计算状态估计值：
Figure BDA0001359618900000117
计算状态估计误差协方差阵：
Figure BDA0001359618900000118
计算***噪声统计特性的估计值：
Figure BDA0001359618900000119
Figure BDA00013596189000001110
Figure BDA00013596189000001111
其中，
Figure BDA00013596189000001112
Figure BDA0001359618900000121
Figure BDA0001359618900000122
式中0＜bq＜1，0＜bQ＜1，0＜bC＜1分别为计算q、Q、C的遗忘因子。随着k的增大，
Figure BDA0001359618900000123
Figure BDA0001359618900000124
和
Figure BDA0001359618900000125
的值分别趋向于(1-bq)、(1-bQ)和(1-bC)，因此bq、bQ和bC越大，当前信息的比重越小。
Figure BDA0001359618900000126
侧重于跟踪***噪声均值q的变化，取值过大会使
Figure BDA0001359618900000127
偏小，取值过小则会引起滤波发散；
Figure BDA0001359618900000128
侧重于平滑
Figure BDA0001359618900000129
的取值，取值过小平滑效果较差，取值过大则使
Figure BDA00013596189000001210
偏小，增大了滤波发散的风险；
Figure BDA00013596189000001211
用来跟踪方差C的变化，取值过小会使估值的波动变大，取值过大则会使跟踪产生延迟。
在实际应用时，虽然式(29)和(30)的计算具有一定的平滑作用，但仍然会出现
Figure BDA00013596189000001212
负定的情况，为了保证滤波的稳定，在负定的情况下对式(29)中的
Figure BDA00013596189000001213
项的对角线元素取绝对值，非对角线元素取0来近似处理。
第二，当***更注重精度指标，而对运算速度和占用存储空间无特别限制时，采用自适应的Sigma点Kalman滤波(Sigma Points Kalman Filter,SPKF)方法对连续非线性***进行状态估计。自适应SPKF的基本思想是将根据确定性采样策略得到的Sigma采样点经非线性函数直接传递，然后在经典卡尔曼滤波算法框架的基础上，通过时变噪声统计估值器，实时估计和修正***噪声和观测噪声的统计特性，再进行状态估计。具体步骤如下：
1)选择Sigma点采样策略
目前常用的Sigma点采样策略有对称采样、单形采样、中心差分采样和容积采样等多种采样方式。选择一种采样策略，计算SPKF的权系数Wi m、Wi c,其中i＝0,1,…,L，L为由具体的采样策略决定的采样点个数，Wi m用于均值估计，Wi c用于方差估计。
2)时间更新：
根据1)中选择的采样策略以及k时刻的状态估计值
Figure BDA00013596189000001214
和估计误差协方差Pk计算k时刻的Sigma点集χi,k(i＝0,1,…,L)。
计算Sigma点的一步预测值：
χi,k+1/k＝F[χi,k,k]，i＝0,1,…,L (34)
其中，χi,k为第i个Sigma点在k时刻的值，χi,k+1/k为第i个Sigma点在k+1时刻的一步预测值。
计算一步预测状态：
Figure BDA0001359618900000131
式中，
Figure BDA0001359618900000132
为k时刻的***噪声均值的估计值，
Figure BDA0001359618900000133
为(k+1)时刻的一步预测状态估计。
计算状态一步预测误差协方差矩阵：
Figure BDA0001359618900000134
式中，Pk+1/k为一步预测状态的误差协方差矩阵，
Figure BDA0001359618900000135
为k时刻***噪声方差阵的估计值；
在滤波初始的时间0～T内，采用噪声的先验统计信息，时间更新如下：
Figure BDA0001359618900000136
Figure BDA0001359618900000137
3)测量更新：
根据1)中选择的采样策略以及
Figure BDA0001359618900000138
和Pk+1/k更新Sigma点集；
估计测量值：
γi，k+1/k＝h(χi，k+1/k) (39)
Figure BDA0001359618900000139
其中γi，k+1/k为由第i个Sigma点计算的第(k+1)时刻测量向量的估计值，
Figure BDA00013596189000001310
为第(k+1)时刻测量向量的加权估计值。
计算增益矩阵：
Figure BDA00013596189000001311
Figure BDA00013596189000001312
Figure BDA00013596189000001313
其中，Rk+1为(k+1)时刻测量噪声协方差阵，Py,k+1/k为第(k+1)时刻测量向量估计值的协方差，Pxy,k+1/k为第(k+1)时测量向量估计值与状态一步预测值的协方差，Kk+1为(k+1)时刻的增益矩阵。
计算新息向量：
Figure BDA0001359618900000141
式中υk+1为(k+1)时刻的新息向量，Zk+1为(k+1)时刻的测量向量；
计算(k+1)时刻的状态估计值
Figure BDA0001359618900000142
Figure BDA0001359618900000143
计算(k+1)时刻的状态估计误差协方差阵Pk+1：
Figure BDA0001359618900000144
4)计算***噪声统计特性的估计值：
Figure BDA0001359618900000145
Figure BDA0001359618900000146
Figure BDA0001359618900000147
其中，
Figure BDA0001359618900000148
Figure BDA0001359618900000149
Figure BDA00013596189000001410
式中0＜bq＜1，0＜bQ＜1，0＜bC＜1分别为计算q、Q、C的遗忘因子。
当出现
Figure BDA00013596189000001411
负定的情况时，对式(48)中的
Figure BDA00013596189000001412
项的对角线元素取绝对值，非对角线元素取0。
步骤四：在集中式结构下，基于选定滤波方法的地球-Lagrange联合星座自主定轨算法的具体实现；
根据本发明提出的自适应非线性滤波算法实现地球-Lagrange联合星座自主定轨算法的流程图参见图4，具体过程如下：
具体过程如下：
(1)数据初始化；
初始化的参数包括：k＝0，滤波步长h，噪声切换参数T，初始状态
Figure BDA0001359618900000151
初始误差协方差阵P0、初始***噪声均值
Figure BDA0001359618900000152
初始***噪声方差阵
Figure BDA0001359618900000153
测量噪声方差阵R以及遗忘因子bq、bQ、bC；
其中，噪声切换参数T为时间更新算法中噪声统计特性切换的时间节点：当k＜T时，使用噪声的先验统计特性进行时间更新；当k≥T时，使用估计的噪声统计特性进行时间更新；初始状态
Figure BDA0001359618900000154
初始误差协方差阵P0、初始***噪声方差阵
Figure BDA0001359618900000155
根据需要确定；初始***噪声均值
Figure BDA0001359618900000156
取零；测量噪声方差阵R根据测量设备的性能确定；遗忘因子bq、bQ取值介于0.99～1.0之间，bC取值介于0.9～0.99之间；
(2)如果采用自适应EKF算法，直接转到(4)；如果采用自适应SPKF算法，转到(3)；
(3)根据选择的Sigma点采样策略计算权系数Wi m、Wi c；
(4)如果k＜T，采用噪声的先验统计特性进行时间更新，转到(5)；否则采用估计的噪声统计特性进行时间更新，转到(6)；
(5)采用噪声的先验统计特性的时间更新；完成后转到(7)；
利用k时刻的状态估计
Figure BDA0001359618900000157
及误差协方差阵Pk作为k时刻的状态初值及误差协方差阵，依次计算状态一步预测值
Figure BDA0001359618900000158
及预测误差协方差矩阵Pk+1/k,实现自适应Kalman滤波估计方法的时间更新。当采用自适应EKF算法时，计算公式分别为式(22)和式(23)；当采用自适应SPKF算法时，计算公式分别为式(37)和(38)；
(6)采用估计的噪声统计特性的时间更新；
利用k时刻的状态估计
Figure BDA0001359618900000159
及误差协方差阵Pk作为k时刻的状态初值及误差协方差阵，依次计算状态一步预测值
Figure BDA00013596189000001510
及预测误差协方差矩阵Pk+1/k,实现自适应Kalman滤波估计方法的时间更新。当采用自适应EKF算法时，计算公式分别为式(18)和式(21)；当采用自适应SPKF算法时，计算公式分别为式(35)和(36)；
(7)测量更新；
依次计算(k+1)时刻的增益矩阵Kk+1、新息矩阵υk+1、状态估计值
Figure BDA0001359618900000161
及误差协方差矩阵Pk+1，实现自适应Kalman滤波的测量更新。当采用自适应EKF算法时，计算公式分别为式(24)～(27)；当采用自适应SPKF算法时，计算公式分别为式(43)～(46)；
(8)估计***噪声统计特性；
依次计算***噪声的均值估计值
Figure BDA0001359618900000162
和协方差估计值
Figure BDA0001359618900000163
当采用自适应EKF算法时，计算公式分别为式(28)和(29)；当采用自适应SPKF算法时，计算公式分别为式(47)和(48)；
(9)判断滤波过程是否结束；如需继续进行滤波，k＝k+1，并返回到(4)。
本发明首先研究仅利用星间距离测量时存在的“亏秩”问题，引入Lagrange卫星来提供绝对基准；然后在集中式结构下，采用自适应非线性滤波估计方法，实现地球-Lagrange联合星座的自主定轨。
本发明提出的一种基于星间测距的地球-Lagrange联合星座自主定轨方法，该方法首先在集中式结构下建立地球-Lagrange联合星座的状态方程；然后建立地球-Lagrange联合星座的测量方程；最后采用自适应非线性滤波方法实现地球-Lagrange联合的自主定轨。
实施例：
以某地球-Lagrange星座为例，该星座由3颗MEO卫星和2颗Lagrange卫星构成。MEO卫星编号为E1～E3；Lagrange卫星编号为L1～L2。在仿真中发现，MEO星座中各卫星的定轨精度近似，Lagrange星座中各卫星的定轨精度也近似，因此文中仅以编号为E1的MEO卫星和编号为L1的Lagrange卫星为例进行说明。
仿真条件：
采用STK来模拟星座实际运行情况，产生星座运行过程中的标称轨道数据，。
(1)标称数据生成参数设置
仿真时间为2016年1月1日12:00:00—2016年6月28日12:00:00，历时180d，步长T＝1s。地球卫星的轨道动力学模型采用高精度轨道预报模型(The High-precision orbitpropagator,HPOP)，其中地球模型采用WGS84-EGM96模型，考虑21×21阶非球形摄动；其他摄动项包括太阳引力、月球引力、太阳光压(光压系数Cr＝1.0)和大气阻力(大气阻力系数Cd＝2.2)；卫星截面积与卫星质量之比S/m＝0.02m2/kg。Lagrange卫星的轨道动力学模型采用圆型限制性三体问题模型。
(2)基本滤波参数设置
地球卫星的状态方程同样采用受摄二体问题模型，但摄动项仅考虑J2项和日月引力；Lagrange卫星的状态方程采用圆型限制性三体问题模型，滤波周期T＝100s。
在集中式结构下，地球-Lagrange联合星座自主定轨的***噪声方差阵初值为Q0＝diag([QE1，0 QE2，0 QE3,0 QL1,0 QL2,0])，其中
Figure BDA0001359618900000171
Figure BDA0001359618900000178
Figure BDA0001359618900000172
σLx0＝σLy0＝σLz0＝10-10[L]，
Figure BDA0001359618900000179
测量噪声方差阵为
Figure BDA0001359618900000173
σEi,Ej＝σEi,Lk＝σLk,Ll＝5m；初始状态估计误差协方差阵P0＝diag([PE1,0 PE2,0 PE3,0 PL1,0 PL2,0])，
其中
Figure BDA0001359618900000174
σEx,0＝σEy,0＝σEz,0＝1m，
Figure BDA00013596189000001710
σLx,0＝σLy,0＝σLz,0＝1m，
Figure BDA00013596189000001711
(3)评估方法
因为实际***中的***状态和估计误差均为随机量，所以在得到各个仿真时刻的估计误差后，采用统计方法(计算误差数据的均方根误差)来评估星座中卫星的绝对定轨和相对定位的效果。采用滤波估计收敛后的估计状态作为样本来统计误差。均方根误差计算公式如下：
Figure BDA0001359618900000176
其中，
Figure BDA0001359618900000177
为第i时刻的估计值，Xi为第i时刻的真值，n为数据总数。
为了便于量化分析，这里分别以三轴位置误差的平方和和三轴速度误差的平方和作为卫星定轨的距离误差和速率误差，其均方根误差的计算公式如下：
卫星距离估计误差：
Figure BDA0001359618900000181
卫星速率估计误差：
Figure BDA0001359618900000182
实施方式：集中式地球-Lagrange联合星座自主定轨***模型的状态向量由五星的位置、速度构成，采用星间距离测量方式，重点研究自适应EKF对定轨精度的提高。
步骤一：建立地球-Lagrange联合星座自主定轨***的状态方程
在地心直角惯性坐标系下，地球导航卫星i的状态向量为
Figure BDA0001359618900000183
其中[xEi,yEi,zEi]T和
Figure BDA0001359618900000184
分别为该星的位置矢量和速度矢量，结合轨道动力学模型式(1)可建立卫星Ei定轨的状态方程：
Figure BDA0001359618900000185
式中，卫星Ei的地球中心引力、J2摄动力和日月引力引起的加速度分别为a0,Ei、aJ2,Ei和aMS，Ei，具体定义参见式(3)～式(5)。WEi(t)表示卫星Ei的***噪声，包含未建模摄动，满足E[WEi(t)]＝0，E[WEi(t)WEi(τ)]＝QEi(t)δ(t-τ)，δ(t-τ)是狄拉克函数，即
Figure BDA0001359618900000186
Lagrange卫星k对应的状态向量为
Figure BDA0001359618900000187
其中[xLk,yLk,zLk]T和
Figure BDA0001359618900000188
分别为该星的无量纲位置矢量和速度矢量，结合轨道动力学模型式(6)可建立卫星Lk定轨的状态方程：
Figure BDA0001359618900000191
式中，WLk(t)表示卫星Ei的***噪声，包含未建模摄动，满足E[WLk(t)]＝0，E[WLk(t)WLk(τ)]＝QLk(t)δ(t-τ)。
根据上述单颗卫星的状态方程，可以建立集中式结构的五星星座状态方程如下：
Figure BDA0001359618900000192
式中，
Figure BDA0001359618900000193
为包含星座中五颗卫星所有状态的状态向量，W(t)表示***噪声，满足E[W(t)]＝q(t)，E[W(t)W(τ)]＝Q(t)δ(t-τ)。
步骤二：建立地球-Lagrange联合星座自主定轨***的测量方程
在集中式结构下，五星星座采用星间链路测量的测量方程如下：
Figure BDA0001359618900000194
式中，hE为三颗地球卫星间的测量方程，hEL为地球卫星和Lagrange卫星间的测量方程，hL为两颗Lagrange卫星间的测量方程，参见式(11)～式(13)。
V(t)＝[ΔEi,Ej…ΔEi,Lk…ΔLk,Ll…]T为经过电离层、对流层延迟等误差补偿之后残余的伪距观测噪声向量，为零均值、标准差为R(t)的白噪声。
步骤三：在集中式结构下，根据上述建立的地球-Lagrange联合星座自主定轨***模型，确定实现轨道参数估计的滤波方法。
本例要求算法的计算速度快和计算资源需求少，因此采用自适应EKF算法。
步骤四：在集中式结构下，采用自适应EKF估计方法实现地球-Lagrange五星星座联合自主定轨，并仿真分析地球-Lagrange五星星座联合自主定轨的性能。
结合地球-Lagrange五星星座集中式联合自主定轨流程图(图4)和发明内容步骤四，采用自适应EKF估计方法实现地球-Lagrange五星星座联合自主定轨。
(1)数据初始化。初始化的参数包括：k＝0，滤波步长h＝100s，阈值T＝(10天×86400s/天)/100s，初始状态
Figure BDA0001359618900000201
初始误差协方差阵P0、初始***噪声均值
Figure BDA0001359618900000202
初始***噪声方差阵
Figure BDA0001359618900000203
测量噪声方差阵R＝5m以及遗忘因子bq＝0.9999、bQ＝0.9999、bC＝0.90。
注：因为bq取值过小将引起滤波发散，所以这里取一个较大值作为初始值，先调节bQ和bC到最佳值，最后调节bq。
(2)如果k＜T，则采用噪声的先验统计特性进行时间更新，转到(3)；否则采用估计的噪声统计特性进行时间更新，转到(4)。
(3)采用噪声的先验统计特性的时间更新。利用k时刻的状态估计
Figure BDA0001359618900000204
及误差协方差阵Pk作为k时刻的状态初值及误差协方差阵，根据式(22)和式(23)分别计算状态一步预测值
Figure BDA0001359618900000205
及预测误差协方差矩阵Pk+1/k,实现自适应EKF估计方法的时间更新。转到(5)。
(4)采用估计的噪声统计特性的时间更新。利用k时刻的状态估计
Figure BDA0001359618900000206
及误差协方差阵Pk作为k时刻的状态初值及误差协方差阵，根据式(18)和式(21)分别计算状态一步预测值
Figure BDA0001359618900000207
及预测误差协方差矩阵Pk+1/k,实现自适应EKF估计方法的时间更新。
(5)自适应EKF测量更新。根据式(24)～(27)依次计算(k+1)时刻的增益矩阵Kk+1、新息矩阵υk+1、状态估计值
Figure BDA0001359618900000208
及误差协方差矩阵Pk+1，实现EKF的测量更新。
(6)估计***噪声统计特性。根据式(28)和(29)分别计算***噪声的均值估计值
Figure BDA0001359618900000209
和协方差估计值
Figure BDA0001359618900000211
(7)判断滤波过程是否结束。如需继续进行滤波，k＝k+1，并返回到(2)。
步骤五：调节遗忘因子bQ、bC、bq。
结合发明内容步骤三中关于三个遗忘因子的作用以及步骤四中给出的遗忘因子的大致范围，依次调节bQ、bC、bq，重复执行步骤三中滤波过程，使***达到相对最优，具体调节步骤如下：
(1)固定bQ＝0.9999、bq＝0.9999，逐渐增大bC，根据滤波结果确定bC的最佳取值，记为
Figure BDA0001359618900000212
(2)固定
Figure BDA0001359618900000213
bq＝0.9999，逐渐减小bQ，根据滤波结果确定bQ的最佳取值，记为
Figure BDA0001359618900000214
(3)固定
Figure BDA0001359618900000215
逐渐减小bq，根据滤波结果确定bq的最佳取值，记为
Figure BDA0001359618900000216
***精度随bQ、bC、bq变化趋势图如图4所示，据图可确定最优遗忘因子组合为：bq＝0.9993、bQ＝0.9997、bC＝0.96。
本发明在分析仅利用星间距离测量时存在的“亏秩”问题的基础上，引入Lagrange卫星以提供绝对基准；然后在集中式结构下，采用自适应EKF估计方法，实现了地球-Lagrange联合星座的自主定轨。
在星间测距误差为零均值，标准差为5m的情况下，得到以下仿真结果：
由于星座中三颗MEO卫星的定轨精度类似，因此以卫星E1为例，其x、y、z三轴定轨距离误差分别为7.78m，7.60m和10.06m，其误差曲线如图5所示，距离估计误差为14.81m；x、y、z三轴定轨速度误差分别为1.43×10-3m/s，3.92×10-3m/s和4.74×10-3m/s，速度估计误差为6.31×10-3m/s。同样，两颗Lagrange的定轨精度类似，以卫星L1为例，其x、y、z三轴定轨距离误差分别为2.66m，3.60m和3.39m，距离估计误差为5.62m；x、y、z三轴定轨速度误差分别为2.11×10-5m/s，1.64×10-5m/s和1.48×10-5m/s，速度估计误差为3.06×10-5m/s。
验证了本发明方法的有效性。

Claims (4)

1.一种基于星间测距的地球-Lagrange联合星座自主定轨方法，包括以下几个步骤：
步骤一：在集中式结构下，建立地球-Lagrange联合星座自主定轨***的状态方程；
具体为：
(1)建立地球导航卫星轨道动力学模型
在地心直角惯性坐标系下，联合星座中待测地球卫星i进行定轨的状态方程如下：
Figure FDA0002569392000000011
式中，待测地球卫星i对应的状态向量为
Figure FDA0002569392000000012
[xEi,yEi,zEi]T和
Figure FDA0002569392000000013
分别为该地球卫星的位置矢量和速度矢量；FE为***状态函数，WEi(t)为***噪声，满足均值为qEi(t)，方差为QEi(t)的高斯白噪声；式可以进一步详细写为：
Figure FDA0002569392000000014
其中，
Figure FDA0002569392000000015
和
Figure FDA0002569392000000016
为x、y、z三轴速度噪声分量，
Figure FDA0002569392000000017
和
Figure FDA0002569392000000018
是x、y、z三轴加速度等效噪声分量，包括未建模摄动项和噪声项；
Figure FDA0002569392000000019
和
Figure FDA00025693920000000110
分别为地球卫星i对应的已建模速度矢量和加速度矢量，其中加速度主要考虑地球中心引力加速度a0,Ei、J2摄动项加速度
Figure FDA00025693920000000111
和日月引力摄动加速度aMS,Ei，具体表达式如下：
1)卫星i的地球中心引力加速度a0,Ei满足：
Figure FDA00025693920000000112
式中，
Figure FDA00025693920000000113
为待测卫星i的地心距，μ为地球质量与引力常数G的乘积，即地球引力常数；
2)J2摄动项加速度
Figure FDA0002569392000000021
满足：
Figure FDA0002569392000000022
式中，Re为地球赤道半径，J2为二阶带谐系数；
日月引力引起的摄动加速度aMS,Ei满足：
Figure FDA0002569392000000023
式中，(xS,yS,zS)和(xM,yM,zM)分别为太阳和月球在地心直角惯性坐标系下的三维位置坐标，
Figure FDA0002569392000000024
和
Figure FDA0002569392000000025
分别为卫星i到太阳和月球的距离，
Figure FDA0002569392000000026
和
Figure FDA0002569392000000027
分别为地心到太阳和月球的距离，μS和μM分别为太阳引力常数和月球引力常数；
(2)Lagrange导航卫星轨道动力学模型
在质心会合坐标系下，建立对Lagrange星座中待测月球卫星k进行定轨的状态方程：
Figure FDA0002569392000000028
式中，待测Lagrange月球卫星k对应的状态向量为
Figure FDA0002569392000000029
[xLk,yLk,zLk]T和
Figure FDA00025693920000000210
分别为该星的无量纲位置矢量和速度矢量，其特征质量、特征长度和特征时间如(7)式所示，其中ME和MM分别为地球质量和月球质量，L为地月距离；FL为***状态函数，WLk(t)为***噪声，满足均值为qLk(t)，方差为QLk(t)的高斯白噪声；
Figure FDA0002569392000000031
式(6)写为：
Figure FDA0002569392000000032
其中，
Figure FDA0002569392000000033
为Lagrange卫星k对应的已建模速度矢量，μ0＝MM/(MM+ME)，ΔE,Lk、ΔM,Lk分别为卫星的地心距和月心距，
Figure FDA0002569392000000034
和
Figure FDA0002569392000000035
为x、y、z三轴速度噪声分量，
Figure FDA0002569392000000036
和
Figure FDA0002569392000000037
是x、y、z三轴加速度等效噪声分量，包括未建模摄动项和噪声项；
(3)地球-Lagrange联合星座自主定轨***的状态方程
假设整个自主导航***包含m颗地球卫星和n颗Lagrange卫星，根据待测地球卫星i进行定轨的状态方程和待测月球卫星k进行定轨的状态方程，建立集中式结构的星座定轨状态方程如下：
Figure FDA0002569392000000038
即：
Figure FDA0002569392000000039
式中，
Figure FDA00025693920000000310
为状态向量，F为整个星座定轨***的状态函数矢量，W(t)表示均值为q(t)，方差为Q(t)的***噪声；
步骤二：建立地球-Lagrange联合星座自主定轨***的测量方程；
步骤三：在集中式结构下，根据上述建立的基于星间距离测量的地球-Lagrange联合星座自主定轨***模型，确定实现轨道参数估计的滤波方法；
步骤四：在集中式结构下，基于选定滤波算法的地球-Lagrange联合星座自主定轨。
2.根据权利要求1所述的一种基于星间测距的地球-Lagrange联合星座自主定轨方法，所述的步骤二，具体为：
(1)地球导航卫星间的星间距离测量模型
对于地球导航卫星，记卫星i和卫星j之间的距离观测量为ρEi,Ej，则有：
ρEi,Ej＝hE(XEi,XEj)+vEi,Ej，i≠j (11)
其中
Figure FDA0002569392000000041
vEi,Ej为卫星i和卫星j之间的等效测量噪声，包括建模误差和量测噪声；
(2)Lagrange导航卫星间的星间距离测量模型
记Lagrange卫星i和j之间的距离观测量为ρLi,Lj，有：
ρLi,Lj＝hL(XLi,XLj)+vLi,Lj，i≠j (12)
其中
Figure FDA0002569392000000042
vLi,Lj为Lagrange卫星i和j之间的等效测量噪声，包括建模误差和量测噪声；
(3)地球导航卫星和Lagrange卫星间的测量模型
将GNSS卫星i和Lagrange卫星k在同一坐标系中表示，则有：
Figure FDA0002569392000000043
其中，
Figure FDA0002569392000000044
是Lagrange卫星在地心惯性坐标系中的坐标，它是关于Lagrange卫星在质心会合坐标系中的无量纲坐标XLk的函数，其转换关系为：
Figure FDA0002569392000000051
其中，R表示由质心会合坐标系到地心惯性坐标系的旋转矩阵，表达式如下：
Figure FDA0002569392000000052
其中uM＝ωM+θM，iM、ΩM、ωM和θM分别表示月球绕地轨道的轨道倾角、升交点赤经、近地点幅角以及真近点角；
(4)地球-Lagrange联合星座自主定轨***的测量方程
根据上述测量方程(11)、(12)和(13)，建立集中式结构的星座定轨测量方程如下：
Z(t)＝h[X(t),t]+V(t) (16)
即：
Figure FDA0002569392000000053
其中V(t)为均值为零，方差为R(t)的等效测量噪声。
3.根据权利要求1所述的一种基于星间测距的地球-Lagrange联合星座自主定轨方法，所述的步骤三，采用EKF方法或者SPKF方法，具体的：
(1)当采用自适应EKF方法对经过离散化后的非线性***进行状态估计，其具体步骤如下：
1)时间更新
计算一步预测状态：
Figure FDA0002569392000000061
式中，
Figure FDA0002569392000000062
为k时刻的***噪声均值的估计值，
Figure FDA0002569392000000063
为利用动力学方程进行的状态一步递推，一种方法是通过离散化和线性化后的方程(19)计算得到，其中
Figure FDA0002569392000000064
为k时刻的状态估计，
Figure FDA0002569392000000065
为(k+1)时刻的一步预测状态估计，Δt为积分步长；另一种方法是通过4阶Runge-Kutta方法直接对连续状态方程(9)数值积分计算得到；
Figure FDA0002569392000000066
计算状态转移矩阵：
Figure FDA0002569392000000067
式中，Φk+1/k为从k时刻到(k+1)时刻的状态转移矩阵；
计算一步预测误差协方差矩阵：
Figure FDA0002569392000000068
式中，Pk+1/k为一步预测状态的误差协方差矩阵，Pk为k时刻估计状态的误差协方差矩阵，
Figure FDA0002569392000000069
为k时刻***噪声方差阵的估计值；
在滤波初始的时间0～T内，采用噪声的先验统计信息，时间更新如下：
Figure FDA00025693920000000610
Figure FDA00025693920000000611
2)测量更新：
计算增益矩阵：
Figure FDA00025693920000000612
式中，Kk+1为(k+1)时刻的增益矩阵，
Figure FDA00025693920000000613
为基于一步预测状态估计
Figure FDA00025693920000000614
计算的观测向量h的Jacobian矩阵，Rk+1为(k+1)时刻测量噪声协方差阵；
计算新息向量：
Figure FDA00025693920000000615
式中υk+1为(k+1)时刻的新息向量，Zk+1为(k+1)时刻的测量向量；
计算状态估计值：
Figure FDA0002569392000000071
计算状态估计误差协方差阵：
Figure FDA0002569392000000072
3)计算***噪声统计特性的估计值：
Figure FDA0002569392000000073
Figure FDA0002569392000000074
Figure FDA0002569392000000075
其中，
Figure FDA0002569392000000076
Figure FDA0002569392000000077
Figure FDA0002569392000000078
式中0＜bq＜1，0＜bQ＜1，0＜bC＜1分别为计算q、Q、C的遗忘因子，随着k的增大，
Figure FDA0002569392000000079
Figure FDA00025693920000000710
和
Figure FDA00025693920000000711
的值分别趋向于(1-bq)、(1-bQ)和(1-bC)；
当出现
Figure FDA00025693920000000712
负定的情况时，对式(29)中的
Figure FDA00025693920000000713
项的对角线元素取绝对值，非对角线元素取0；
(2)当采用自适应SPKF方法对连续非线性***进行状态估计，其具体步骤如下：
1)选择Sigma点采样策略
计算Sigma点Kalman滤波的权系数Wi m、Wi c，其中i＝0,1,…,L，L为采样点个数，Wi m用于均值估计，Wi c用于方差估计；
2)时间更新：
根据1)中选择的采样策略以及k时刻的状态估计值
Figure FDA00025693920000000714
和估计误差协方差Pk计算k时刻的Sigma点集χi,k；
计算Sigma点的一步预测值：
χi,k+1/k＝F[χi,k,k]，i＝0,1,…,L (34)
其中，χi,k为第i个Sigma点在k时刻的值，χi,k+1/k为第i个Sigma点在k+1时刻的一步预测值；
计算一步预测状态：
Figure FDA0002569392000000081
式中，
Figure FDA0002569392000000082
为k时刻的***噪声均值的估计值，
Figure FDA0002569392000000083
为(k+1)时刻的一步预测状态估计；
计算状态一步预测误差协方差矩阵：
Figure FDA0002569392000000084
式中，Pk+1/k为一步预测状态的误差协方差矩阵，
Figure FDA0002569392000000085
为k时刻***噪声方差阵的估计值；
在滤波初始的时间0～T内，采用噪声的先验统计信息，时间更新如下：
Figure FDA0002569392000000086
Figure FDA0002569392000000087
3)测量更新：
根据1)中选择的采样策略以及
Figure FDA0002569392000000088
和Pk+1/k更新Sigma点集；
估计测量值：
γi,k+1/k＝h(χi,k+1/k) (39)
Figure FDA0002569392000000089
其中：γi,k+1/k为由第i个Sigma点计算的第(k+1)时刻测量向量的估计值，
Figure FDA00025693920000000810
为第(k+1)时刻测量向量的加权估计值；
计算增益矩阵：
Figure FDA00025693920000000811
Figure FDA00025693920000000812
Figure FDA0002569392000000091
其中，Rk+1为(k+1)时刻测量噪声协方差阵，Py,k+1/k为第(k+1)时刻测量向量估计值的协方差，Pxy,k+1/k为第(k+1)时测量向量估计值与状态一步预测值的协方差，Kk+1为(k+1)时刻的增益矩阵；
计算新息向量：
Figure FDA0002569392000000092
式中υk+1为(k+1)时刻的新息向量，Zk+1为(k+1)时刻的测量向量；
计算(k+1)时刻的状态估计值
Figure FDA0002569392000000093
Figure FDA0002569392000000094
计算(k+1)时刻的状态估计误差协方差阵Pk+1：
Figure FDA0002569392000000095
4)计算***噪声统计特性的估计值：
Figure FDA0002569392000000096
Figure FDA0002569392000000097
Figure FDA0002569392000000098
其中，
Figure FDA0002569392000000099
Figure FDA00025693920000000910
Figure FDA00025693920000000911
式中0＜bq＜1，0＜bQ＜1，0＜bC＜1分别为计算q、Q、C的遗忘因子；
当出现
Figure FDA00025693920000000912
负定的情况时，对式(48)中的
Figure FDA00025693920000000913
项的对角线元素取绝对值，非对角线元素取0。
4.根据权利要求1所述的一种基于星间测距的地球-Lagrange联合星座自主定轨方法，所述的步骤四，具体为：
(1)数据初始化；
初始化的参数包括：k＝0，滤波步长h，噪声切换参数T，初始状态
Figure FDA0002569392000000101
初始误差协方差阵P0、初始***噪声均值
Figure FDA0002569392000000102
初始***噪声方差阵
Figure FDA0002569392000000103
测量噪声方差阵R以及遗忘因子bq、bQ、bC；
其中，噪声切换参数T为时间更新算法中噪声统计特性切换的时间节点：当k＜T时，使用噪声的先验统计特性进行时间更新；当k≥T时，使用估计的噪声统计特性进行时间更新；初始状态
Figure FDA0002569392000000104
初始误差协方差阵P0、初始***噪声方差阵
Figure FDA0002569392000000105
根据需要确定；初始***噪声均值
Figure FDA0002569392000000106
取零；测量噪声方差阵R根据测量设备的性能确定；遗忘因子bq、bQ取值介于0.99～1.0之间，bC取值介于0.9～0.99之间；
(2)如果采用自适应EKF算法，直接转到(4)；如果采用自适应SPKF算法，转到(3)；
(3)根据选择的Sigma点采样策略计算权系数Wi m、Wi c；
(4)如果k＜T，采用噪声的先验统计特性进行时间更新，转到(5)；否则采用估计的噪声统计特性进行时间更新，转到(6)；
(5)采用噪声的先验统计特性的时间更新；完成后转到(7)；
利用k时刻的状态估计
Figure FDA0002569392000000107
及误差协方差阵Pk作为k时刻的状态初值及误差协方差阵，依次计算状态一步预测值
Figure FDA0002569392000000108
及预测误差协方差矩阵Pk+1/k,实现自适应Kalman滤波估计方法的时间更新；当采用自适应EKF算法时，计算公式分别为式(22)和式(23)；当采用自适应Sigma点Kalman滤波算法时，计算公式分别为式(37)和(38)；
(6)采用估计的噪声统计特性的时间更新；
利用k时刻的状态估计
Figure FDA0002569392000000109
及误差协方差阵Pk作为k时刻的状态初值及误差协方差阵，依次计算状态一步预测值
Figure FDA00025693920000001010
及预测误差协方差矩阵Pk+1/k,实现自适应Kalman滤波估计方法的时间更新；当采用自适应EKF算法时，计算公式分别为式(18)和式(21)；当采用自适应SPKF算法时，计算公式分别为式(35)和(36)；
(7)测量更新；
依次计算(k+1)时刻的增益矩阵Kk+1、新息矩阵υk+1、状态估计值
Figure FDA00025693920000001011
及误差协方差矩阵Pk+1，实现自适应Kalman滤波的测量更新；当采用自适应EKF算法时，计算公式分别为式(24)～(27)；当采用自适应SPKF算法时，计算公式分别为式(43)～(46)；
(8)估计***噪声统计特性；
依次计算***噪声的均值估计值
Figure FDA0002569392000000111
和协方差估计值
Figure FDA0002569392000000112
当采用自适应EKF算法时，计算公式分别为式(28)和(29)；当采用自适应SPKF算法时，计算公式分别为式(47)和(48)；
(9)判断滤波过程是否结束；如需继续进行滤波，k＝k+1，并返回到(4)。
CN201710611461.6A 2017-07-25 2017-07-25 一种基于星间测距的地球-Lagrange联合星座自主定轨方法 Active CN107421550B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201710611461.6A CN107421550B (zh) 2017-07-25 2017-07-25 一种基于星间测距的地球-Lagrange联合星座自主定轨方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201710611461.6A CN107421550B (zh) 2017-07-25 2017-07-25 一种基于星间测距的地球-Lagrange联合星座自主定轨方法

Publications (2)

Publication Number Publication Date
CN107421550A CN107421550A (zh) 2017-12-01
CN107421550B true CN107421550B (zh) 2020-08-28

Family

ID=60431049

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201710611461.6A Active CN107421550B (zh) 2017-07-25 2017-07-25 一种基于星间测距的地球-Lagrange联合星座自主定轨方法

Country Status (1)

Country Link
CN (1) CN107421550B (zh)

Families Citing this family (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108287334A (zh) * 2018-02-06 2018-07-17 西安四方星途测控技术有限公司 一种基于rcs测量数据的自旋卫星姿态估计方法及***
CN109031349B (zh) * 2018-04-20 2022-04-08 南京航空航天大学 一种geo卫星的智能自主运行***
CN109752006B (zh) * 2018-11-23 2022-09-09 中国西安卫星测控中心 一种非完备外测数据在实时滤波中的使用方法
CN109682383B (zh) * 2018-11-23 2022-11-04 中国西安卫星测控中心 一种使用深空三向测量距离和数据的实时滤波定位方法
CN109933847B (zh) * 2019-01-30 2022-09-16 中国人民解放军战略支援部队信息工程大学 一种改进的主动段弹道估计算法
CN109917431B (zh) * 2019-04-02 2021-03-23 中国科学院空间应用工程与技术中心 一种利用dro轨道和星间测量实现gnss卫星自主导航的方法
CN109975839B (zh) * 2019-04-10 2023-04-07 华砺智行(武汉)科技有限公司 一种车辆卫星定位数据的联合滤波优化方法
CN110793528B (zh) * 2019-09-27 2021-04-13 西安空间无线电技术研究所 一种基于低轨星基锚固的北斗导航星座自主定轨方法
CN111522036B (zh) * 2020-04-30 2022-07-26 中国科学院微小卫星创新研究院 星上可用的北斗卫星集中式星座自主导航***及导航方法
CN111947667B (zh) * 2020-06-24 2022-08-12 火眼位置数智科技服务有限公司 一种基于运动学和动力学组合的低轨卫星实时高精度定轨方法
CN112504281B (zh) * 2020-11-16 2023-06-27 中国人民解放军63921部队 基于北斗星间单向链路的航天器测定轨方法
CN113008243B (zh) * 2021-02-23 2024-03-19 重庆两江卫星移动通信有限公司 一种低轨卫星星座自主定轨方法及***
CN113204917B (zh) * 2021-04-25 2021-12-07 中国科学院国家空间科学中心 针对geo目标的天基光学测角弧段初定轨及关联方法

Family Cites Families (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US9091552B2 (en) * 2011-10-25 2015-07-28 The Boeing Company Combined location and attitude determination system and methods
WO2014172675A1 (en) * 2013-04-18 2014-10-23 California Institute Of Technology Real-time and post-processed orbit determination and positioning
US9939260B2 (en) * 2014-08-28 2018-04-10 The Boeing Company Satellite transfer orbit search methods
CN105716612B (zh) * 2016-02-29 2017-05-10 武汉大学 一种捷联惯导***模拟器的设计方法
CN106338753B (zh) * 2016-09-22 2019-03-12 北京航空航天大学 一种基于地面站/星间链路/gnss联合测量的地球同步轨道星座定轨方法
CN106885577B (zh) * 2017-01-24 2020-01-21 南京航空航天大学 拉格朗日导航卫星自主定轨方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
"导航卫星自主导航关键技术研究";肖寅;《中国博士学位论文全文数据库（电子期刊）》;20160115;正文第24-29页 *

Also Published As

Publication number Publication date
CN107421550A (zh) 2017-12-01

Similar Documents

Publication Publication Date Title
CN107421550B (zh) 一种基于星间测距的地球-Lagrange联合星座自主定轨方法
CN109255096B (zh) 一种基于微分代数的地球同步卫星轨道不确定演化方法
Xiangyu et al. An autonomous optical navigation and guidance for soft landing on asteroids
CN107797130A (zh) 低轨航天器多点多参数轨道上行数据计算方法
Takahashi et al. Spin state and moment of inertia characterization of 4179 Toutatis
CN111552003B (zh) 基于球卫星编队的小行星引力场全自主测量***及方法
CN112161632B (zh) 一种基于相对位置矢量测量的卫星编队初始定位方法
CN111522037A (zh) 星座同轨道面卫星自主导航方法及导航***
Hesar et al. Lunar far side surface navigation using linked autonomous interplanetary satellite orbit navigation (LiAISON)
CN104048664A (zh) 一种导航卫星星座自主定轨的方法
CN104848862A (zh) 一种环火探测器精密同步定位守时方法及***
CN109752005A (zh) 一种基于精确轨道模型的航天器初轨确定方法
Lee et al. Vision-based relative state estimation using the unscented Kalman filter
CN111731513B (zh) 一种基于单脉冲轨控的高精度引力场中回归轨道维持方法
Ke et al. Pico-satellite autonomous navigation with magnetometer and sun sensor data
Qian et al. Sun–Earth–Moon autonomous orbit determination for quasi-periodic orbit about the translunar libration point and its observability analysis
Gui et al. A time delay/star angle integrated navigation method based on the event-triggered implicit unscented Kalman filter
CN113589832B (zh) 对地表固定区域目标稳定观测覆盖的星座快速设计方法
Babcock CubeSat attitude determination via Kalman filtering of magnetometer and solar cell data
Lou et al. A consider unscented particle filter with genetic algorithm for UAV multi-source integrated navigation
Baroni Attitude determination by unscented Kalman filter and solar panels as sun sensor
Liu et al. Guidance and control technology of spacecraft on elliptical orbit
Nordkvist et al. Attitude feedback tracking with optimal attitude state estimation
Liuqing et al. An autonomous navigation study of Walker constellation based on reference satellite and inter-satellite distance measurement
CN111854765B (zh) 一种中轨道导航卫星轨道长期预报方法

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant