CN114063122B - 电推进转移轨道航天器星载gnss在轨实时定轨方法 - Google Patents

电推进转移轨道航天器星载gnss在轨实时定轨方法 Download PDF

Info

Publication number
CN114063122B
CN114063122B CN202111456809.1A CN202111456809A CN114063122B CN 114063122 B CN114063122 B CN 114063122B CN 202111456809 A CN202111456809 A CN 202111456809A CN 114063122 B CN114063122 B CN 114063122B
Authority
CN
China
Prior art keywords
orbit
satellite
pseudo
bds
gps
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
CN202111456809.1A
Other languages
English (en)
Other versions
CN114063122A (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.)
Wuhan University WHU
Space Star Technology Co Ltd
Original Assignee
Wuhan University WHU
Space Star Technology Co Ltd
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 Wuhan University WHU, Space Star Technology Co Ltd filed Critical Wuhan University WHU
Priority to CN202111456809.1A priority Critical patent/CN114063122B/zh
Publication of CN114063122A publication Critical patent/CN114063122A/zh
Application granted granted Critical
Publication of CN114063122B publication Critical patent/CN114063122B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01—Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13—Receivers
    • G01S19/35—Constructional details or hardware or software details of the signal processing chain
    • G01S19/37—Hardware or software details of the signal processing chain

Landscapes

  • Engineering & Computer Science (AREA)
  • Signal Processing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Remote Sensing (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Position Fixing By Use Of Radio Waves (AREA)

Abstract

本发明公开了一种电推进转移轨道航天器星载GNSS在轨实时定轨方法，属于航天器自主轨道测定领域，本发明通过GNSS接收机实时获取GNSS观测数据、卫星姿态以及电推力广播数据，根据电推力器的安装方向，将电推力产生的本体系加速度转换到惯性系，在伪距实时定轨求解卫星运动方程中顾及电推力产生的惯性系加速度的影响，利用伪距观测值的验前残差检验先验轨道信息的准确性，当判别先验轨道信息不准确时，可通过调整补偿加速度的过程噪声来进一步保证自主定轨的稳定输出，不受轨道机动中电推力建模误差的影响。采用本发明可以克服转移轨道航天器受电推进轨道机动导致常规自主定轨滤波发散的问题，定轨计算稳定，提高了自主定轨的可靠性。

Description

电推进转移轨道航天器星载GNSS在轨实时定轨方法
技术领域
本发明涉及航天器自主轨道测定领域，尤其是涉及一种电推进转移轨道航天器星载GNSS在轨实时定轨方法。
背景技术
转移轨道航天器是一类轨道可变的特殊类型航天器，比较典型如地球同步转移轨道(GTO)，月球探测器在发射后的地月转移轨道等。目前完成此类航天器的轨道转移，主要采用的推进系统有化学推进、冷气推进、电推进等，由于离子电推进具有宽范围连续精细调节、低噪声、长寿命和高比冲等特点，经常被用于转移轨道航天器的任务中。
轨道转移航天器高精度轨道确定是保证科学任务顺利完成的关键，因此对其导航系统的精度及可靠性提出了较高的要求。星载GNSS因其具有可连续观测、精度高、成本功耗低、体积小重量轻等优点，是完成航天器高精度轨道确定的主要技术手段。传统的GNSS伪距自主定轨只适用于无轨道机动的正常工况下，对于有离子电推进的轨道机动工况，会出现定轨滤波发散情形，因而不能保证自主定轨的精度、连续性及可靠性。此外，轨道转移航天器因轨道高度较高，星载GNSS信号多为弱信号，对GNSS弱信号跟踪捕获较为困难，难免存在伪距粗差多、粗差量级大的情况，如果不能有效识别伪距粗差，也必将影响定轨精度，严重时还将引起定轨滤波发散。这些都对GNSS伪距自主定轨数据处理也带来了挑战。
综上所述，传统GNSS伪距自主定轨方法无法满足有离子电推进情形的轨道转移航天器全过程高精度高可靠性轨道测控的需求，因而亟需一种顾及电推进影响的转移轨道航天器星载GNSS在轨实时定轨的数据处理方法。
发明内容
本发明的目的是提供一种电推进转移轨道航天器星载GNSS在轨实时定轨方法，可以克服转移轨道航天器受电推进轨道机动导致常规自主定轨滤波发散的问题，定轨计算稳定，提高了自主定轨的可靠性。
为实现上述目的，本发明提供了一种电推进转移轨道航天器星载GNSS在轨实时定轨方法，包括以下具体步骤：
步骤S1：实时获取GNSS观测数据、卫星姿态以及电推力广播数据；
步骤S2：当判断先验轨道信息有效标记有效时，计算并统计验前伪距观测值残差的均值和标准差，当判断先验轨道信息有效标记无效时，跳转到步骤S4；
步骤S3：根据步骤S2中得到的伪距观测值残差的均值和标准差信息判断先验轨道信息是否准确和接收机各通道是否含有伪距粗差，当先验轨道信息不准确时，调整补偿加速度的过程噪声；
步骤S4：完成顾及伪距粗差的伪距单点定轨；
步骤S5：完成顾及电推力影响的定轨时间更新过程；
步骤S6：完成顾及伪距粗差的定轨测量更新过程；
步骤S7：完成先验轨道的预报及先验轨道精度计算，标识先验轨道信息有效；
步骤S8：输出定轨滤波结果和先验轨道预报结果。
进一步的，步骤S2中，利用步骤S7中计算得到的先验轨道预报结果代入GNSS伪距观测方程来计算伪距观测值的验前残差，并根据导航系统类型的不同分别计算每个导航系统的伪距残差序列的均值和标准差，然后计算所有通道的伪距残差序列的均值和标准差。
进一步的，在步骤S3中，判断先验轨道信息是否准确的方法如下：
当连续多次满足伪距验前残差标准差大于伪距观测值的噪声时，则判断先验轨道信息是准确的，否则，不准确；
其中伪距观测值的噪声大小根据转移轨道航天器的轨道高度分段自适应调整，且伪距观测值的噪声阈值设置是通过地面站参数上注的方式进行修改。
进一步的，在步骤S5中根据电推力器的安装方向，使用卫星姿态数据将电推力产生的本体系加速度转换到惯性系下，在伪距实时定轨使用数值积分方法求解卫星变分方程时通过电推力引起的加速度进行修订，所述的电推力引起的加速度的计算公式如下：
Figure GDA0003698627160000031
其中，
Figure GDA0003698627160000032
为卫星本体系到地心惯性系的转换矩阵并由姿态广播数据计算得到；
Figure GDA0003698627160000033
为电推力参考坐标系到卫星本体系的转换矩阵，由离子电推进器在卫星本体系具体安装方向计算得到；
Figure GDA0003698627160000034
为离子电推进器在三个方向上的电推矢量，由电推力广播数据得到，m为卫星质量参数，单位为kg。
因此，本发明采用上述一种电推进转移轨道航天器星载GNSS在轨实时定轨方法，具有以下有益效果：
(1)、通过本发明中的一种顾及多通道伪距粗差的伪距验前残差统计方法，可有效探测识别GNSS弱信号中可多通道伪距粗差，从而提高自主定轨精度和可靠性。
(2)、本发明利用伪距验前残差统计信息检验先验轨道信息的准确性，当判别先验轨道信息不准确时，可通过调整补偿加速度的过程噪声来进一步保证自主定轨的稳定输出，不受轨道机动中电推力建模误差的影响。
(3)、本发明可以克服转移轨道航天器受电推进轨道机动导致常规自主定轨滤波发散的问题，定轨计算稳定，不仅提高了自主定轨的可靠性，还扩展了其适用范围。
下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。
附图说明
图1是本发明实施例提供的一种电推进转移轨道航天器星载GNSS在轨实时定轨方法流程图；
图2是本发明实施例提供的一种顾及多通道伪距粗差的伪距验前残差统计方法流程图；
图3是本发明实施例提供的单个系统残差序列统计信息的方法流程图；
图4是本发明实施例提供的判断先验轨道信息是否准确和接收机各通道是否含有伪距粗差的方法流程图；
图5是本发明实施例提供的半实物仿真测试验证方案示意图。
具体实施方式
实施例
如图1所示，是本发明实施例提供的一种电推进转移轨道航天器星载GNSS在轨实时定轨方法流程图，包括以下步骤：
步骤S1：实时获取GNSS观测数据、卫星姿态以及电推力广播数据。
所述的GNSS观测数据包括GNSS伪距、多普勒等观测值，是由高灵敏GNSS接收机对在轨实收的GNSS信号下变频、跟踪以及捕获得到，所述的卫星姿态、电推力广播数据由高灵敏GNSS接收机实时接收并解码卫星星务系统通过CAN总线播发的姿态、电推力广播信息获取得到。
步骤S2：计算并统计验前伪距观测值残差的均值和标准差。
判断先验轨道信息有效标记是否有效，如果先验轨道信息有效，则利用GNSS伪距观测方程计算伪距观测值的验前残差，并根据导航系统类型的不同分别计算每个导航系统的伪距残差序列的均值和标准差，然后计算所有通道的伪距残差序列的均值和标准差。如果先验轨道信息无效，直接跳至步骤S4。
步骤S3：判断先验轨道信息是否准确和接收机各通道是否含有伪距粗差；根据步骤S2中所述伪距验前残差的均值和标准差信息，判断先验轨道信息是否准确和接收机各通道是否含有伪距粗差，如先验轨道信息不准确，则需调整补偿加速度的过程噪声。所述的先验轨道信息是由步骤S7中计算得到。
步骤S4：完成顾及伪距粗差的伪距单点定轨；
以GPS/BDS卫星为例，星载GPS/BDS伪距观测方程可表示为：
Figure GDA0003698627160000051
式(1)中，
Figure GDA0003698627160000052
分别为经过电离层改正后的星载GPS/BDS的伪距观测值；δtG和δtC分别为GPS/BDS接收机钟差；c为真空中的光速；
Figure GDA0003698627160000053
和
Figure GDA0003698627160000054
分别为第i颗GPS卫星和第j颗BDS卫星的钟差，可分别由GPS/BDS广播星历计算得到；
Figure GDA0003698627160000055
和
Figure GDA0003698627160000056
分别为单频GPS/BDS电离层延迟改正值，
Figure GDA0003698627160000057
和
Figure GDA0003698627160000058
分别为GPS/BDS伪距观测值噪声；
Figure GDA0003698627160000059
为相应伪距观测值对应的站星几何距离，用下式表示：
Figure GDA00036986271600000510
式(2)中，(x,y,z)分别为GNSS接收机在地心地固系内三维位置坐标向量，
Figure GDA00036986271600000511
分别为第i颗GPS卫星和第j颗BDS卫星在地心地固系内三维位置坐标向量，可分别由GPS/BDS广播星历计算得到。
假定GNSS接收机的位置坐标(x,y,z)和接收机钟差(cδtG,cδtC)组成的向量为X＝(x,y,z,cδtG,cδtC)T，其初值和改正值分别为
Figure GDA00036986271600000512
和ΔX＝(dx,dy,dz,dcδtG,dcδtC)T，则GPS/BDS伪距观测值
Figure GDA00036986271600000513
关于X的偏导数为：
Figure GDA0003698627160000061
如果t时刻星载GPS/BDS接收机共观测到m颗GPS和n颗BDS卫星，假设步骤S3中第k颗GPS卫星和第l颗BDS卫星被标识含有伪距粗差，则这两颗卫星不参与计算观测值矩阵A和残差向量b，其计算公式如下：
Figure GDA0003698627160000062
式(4)中，
Figure GDA0003698627160000063
和
Figure GDA0003698627160000064
分别为第i颗GPS卫星和第j颗BDS卫星伪距观测值
Figure GDA0003698627160000065
关于X的偏导数，初始计算时，是将X＝X0带入式(3)；
Figure GDA0003698627160000066
分别为GPS卫星实测的伪距观测值和伪距观测值的计算值；
Figure GDA0003698627160000067
分别为BDS卫星实测的伪距观测值和伪距观测值的计算值；
可采用最小二乘原理，计算星载GNSS接收机初始坐标和接收机钟差的改正值为ΔX＝(ATA)-1ATb，使用ΔX来更新X，即X＝X0+ΔX，再计算新的观测值矩阵A和残差向量b后，重新计算ΔX，直到满足迭代终止条件
Figure GDA0003698627160000068
步骤S5：完成顾及电推力影响的定轨时间更新过程。
在地心惯性系中，卫星(转移轨道航天器)运动方程可以用一阶微分方程组来表示：
Figure GDA0003698627160000069
式(5)中，
Figure GDA0003698627160000071
分别为卫星的位置、速度；
Figure GDA0003698627160000072
为卫星受到的各种摄动加速度之和，可以表示为：
Figure GDA0003698627160000073
式(6)中，
Figure GDA0003698627160000074
为近地航天器所受的惯性系总加速度；
Figure GDA0003698627160000075
为由保守力(包括地球中心引力、非球形引力、N体引力、地球固体潮汐和海洋潮汐摄动力等)引起的加速度；
Figure GDA0003698627160000076
为由非保守力(包括大气阻力、太阳光压力等)引起的加速度；
Figure GDA0003698627160000077
为电推力引起的加速度；
Figure GDA0003698627160000078
为为人为引入的经验补偿加速度，用于补偿无法模型化或错误模型的微小摄动力的影响，采用一阶高斯-马尔可夫随机模型对径向(Radial，R)、切向(Along，A)、法向(Cross，C)3个方向进行动力学模型补偿；T是RTN坐标系到地心惯性系的转换矩阵。
由电推力引起的加速度
Figure GDA0003698627160000079
计算如下式所示：
Figure GDA00036986271600000710
式(7)中，
Figure GDA00036986271600000711
为卫星本体系到地心惯性系的转换矩阵，可由步骤S1实时获取的姿态广播数据计算，例如，姿态广播数据为姿态四元数(q1,q2,q3,q4)，则
Figure GDA00036986271600000712
Figure GDA00036986271600000713
为电推力参考坐标系到卫星本体系的转换矩阵，可由离子电推进器在卫星本体系具体安装方向计算得到；
Figure GDA00036986271600000714
为离子电推进器在三个方向上的电推力矢量，单位为N，由步骤S1实时获取的电推力广播数据得到，一般地，电推力广播数据误差<1％，当没有电推力作用时，电推力广播值为0；m为卫星质量参数，单位为kg，可以通过地面参数上注的方式修改；
考虑到在一般的精密定轨中，除了要考虑位置
Figure GDA0003698627160000081
参数和速度
Figure GDA0003698627160000082
参数之外，还要考虑动力学参数
Figure GDA0003698627160000083
(太阳光压模型参数、大气阻力参数等)，且有
Figure GDA0003698627160000084
因此将式(5)进行扩展，并令
Figure GDA0003698627160000085
则有
Figure GDA0003698627160000086
当设定一初始状态
Figure GDA0003698627160000087
后，通过运动方程积分就可得到卫星的参考状态向量
Figure GDA0003698627160000088
该初值可以从步骤S4中得到星载GNSS接收机位置作为粗略轨道获得。由于上面的状态方程是非线性的，所以对
Figure GDA0003698627160000089
在
Figure GDA00036986271600000810
上展开取一次项，可得：
Figure GDA00036986271600000811
令
Figure GDA00036986271600000812
则有：
Figure GDA00036986271600000813
其解可以写为：
Figure GDA00036986271600000814
则有：
Figure GDA00036986271600000815
式(7)称为变分方程，其中I为单位矩阵，Φ(t,t0)被称为状态转移矩阵。
一般的动力学定轨中，通过数值积分方法，按式(8)和式(11)，可以同时求解卫星的运动方程和变分方程。
伪距定轨卡尔曼滤波的状态方程为：
Xk＝Φk,k-1Xk-1+Wk (12)
式(5)中,Xk＝(r,v,cδtG,cδtC,Cd,Cr,ω)T为伪距定轨卡尔曼滤波状态量，Φk,k-1为离散形式的状态转移矩阵，Wk为系统噪声矩阵。
Figure GDA0003698627160000091
式(13)中，Φrr,Φrv,Φvr,Φvv分别位置和速度分量的状态转移矩阵；
Figure GDA0003698627160000092
位置和速度分别关于大气阻力和太阳光压的状态转移矩阵；Φrw,Φvw分别位置和速度关于RTN方向的补偿加速度的状态转移矩阵，Φww为补偿加速度关于其自身的状态转移矩阵。
伪距定轨卡尔曼滤波的观测方程为：
Zk＝HkXk-1+Vk (14)
式(14)中，Zk为观测矢量，Hk为观测矩阵，Vk为观测噪声矢量，且系统噪声Wk和观测噪声Vk，为零均值白噪声序列，Qk和Rk分别为系统动态噪声协方差阵和观测噪声协方差阵。
伪距定轨时间更新具体过程是，通过数值积分方法，一般采用单步法，比如龙格-库塔(Runge-Kutta，RK)积分法，同时求解卫星的运动方程和变分方程，得到Φk,k-1，再按式-(12)完成对滤波状态量的时间更新，再完成对状态误差协方差矩阵的时间更新；
Figure GDA0003698627160000093
式(15)中，
Figure GDA0003698627160000094
为k-1时刻的状态误差协方差矩阵，
Figure GDA0003698627160000095
为k时刻的状态误差协方差矩阵。
步骤S6：完成顾及伪距粗差的定轨测量更新过程；
式(3)中伪距观测值关于伪距定轨滤波器状态量的偏导数H可以表示为：
Figure GDA0003698627160000101
式(16)中，UT为J2000惯性系到地固系的转换矩阵，
Figure GDA0003698627160000102
为航天器相对于GPS卫星i和BDS卫星j的视线向量。
顾及伪距粗差的定轨测量更新过程为：按式(17)，依次对于步骤S3中没有被标识含有伪距粗差的GPS/BDS伪距观测值进行处理，更新绝对定轨卡尔曼滤波的状态量及状态协方差阵；
Figure GDA0003698627160000103
式(16)中，
Figure GDA0003698627160000104
为当前k时刻第i颗GPS/BDS对应的卡尔曼增益矩阵；
Figure GDA0003698627160000105
为式(16)中的观测矩阵；
Figure GDA0003698627160000106
为观测噪声协方差阵，
Figure GDA0003698627160000107
为第i颗GPS/BDS测量更新后的滤波状态量；
Figure GDA0003698627160000108
为步骤S5中当前k时刻时间更新后的滤波状态量，
Figure GDA0003698627160000109
第i颗GPS/BDS实测的伪距观测值；
Figure GDA00036986271600001010
为式(15)中当前k时刻时间更新后的状态误差协方差矩阵；
Figure GDA00036986271600001011
为当前k时刻测量更新后的状态误差协方差矩阵。
步骤S7：完成先验轨道的预报及先验轨道精度计算，并标识先验轨道信息有效。
在步骤S6中定轨测量更新最后得到的滤波状态量
Figure GDA00036986271600001012
的基础上，通过数值积分方法，一般采用单步法，比如龙格-库塔(Runge-Kutta，RK)积分法，按式(10)求解卫星运动方法，得到下一个时刻的先验轨道PreOrb，先验轨道精度Accu由在步骤S6中状态误差协方差矩阵
Figure GDA0003698627160000111
计算；当步骤S6中参与测量更新的卫星数大于0时，标识先验轨道信息有效，否则，标识先验轨道信息无效。
步骤S8：输出定轨滤波结果和先验轨道预报结果。
将步骤S6中得到的滤波状态量
Figure GDA0003698627160000112
和步骤S7中得到先验轨道PreOrb及其精度Accu作为输出，供下个时刻定轨滤波周期使用。
如图2所示，是本发明实施例提供的一种顾及多通道伪距粗差的伪距验前残差统计方法流程图，包括以下子步骤：
步骤S2-1：由先验轨道信息计算全部的残差序列Val。
首先将步骤S7中得到的先验轨道信息PreOrb代入式(2)中先计算站星几何距离，然后代入伪距观测值方程式(1)中得到每颗GPS/BDS卫星计算的伪距观测值
Figure GDA0003698627160000115
最终得到所有卫星的残差序列
Figure GDA0003698627160000113
其中Pi为第i颗GPS/BDS卫星实测的伪距观测值，m和n分别为当前历元观测到的GPS和BDS的卫星数。
步骤S2-2：按导航系统类型的不同，将全部的残差序列Val分为GPS残差序列GVal和BDS残差序列BVal。
Figure GDA0003698627160000114
步骤S2-3：根据计算单个系统残差序列统计信息的方法，分别计算残差序列GVal和BVal的均值和标准差GMean、GStd和BMean、BStd。所述的计算单个系统残差序列统计信息的方法充分顾及了有多通道伪距粗差的情况，其流程图详见图3。
步骤S2-4：将残差序列GVal和BVal分别减去各自均值GMean和Bmean。
Figure GDA0003698627160000121
步骤S2-5：合并残差序列GVal和BVal得到新的残差序列Val。即将步骤步骤S2-4中减去各自均值后的残差序列GVal和BVal合并得到新的残差序列Val。
步骤S2-6：根据计算单个系统残差序列统计信息的方法，计算新残差序列Val的均值Mean和标准差Std。即将步骤S2-5中得到的新的残差序列Val，当做单个系统，按步骤S2-3中计算单个系统残差序列统计信息的方法，重新计算新残差序列Val的均值Mean和标准差Std。
步骤S2-7：设置新残差序列Val的标准差下限值Std_min。设置规则如下：
Figure GDA0003698627160000122
如图3所示，是本发明实施例提供的单个系统残差序列统计信息的方法流程图，是本发明实施例提供的一种顾及多通道伪距粗差的伪距验前残差统计方法步骤S2-3的主要内容，包括以下子步骤：
步骤S2-3-1：对单个系统所有通道残差序列进行由大到小排序。所述的残差序列排序方法可采用实数冒泡排序法，排序后的残差序列记为Val＝(dP1 ... dPi ... dPn)T,i＝1,...,n。其中有dP1≥...dPi≥...dPn。
步骤S2-3-2：取排序后的残差序列中间位置m的元素作为初始均值Mean0。将设步骤S2-3-1中排序后的残差序列Val维数为n，则m＝n/2，其中m和n均为整数。
步骤S2-3-3：设置初始标准差Std0。可以下列规则设置残差序列Val的初始标准差：
Figure GDA0003698627160000123
式(21)中，Accu为步骤S7中的先验轨道精度。
步骤S2-3-4：以Mean0和Std0计算不满足条件1门限的通道个数k1。所述的条件1为：|dPi-Mean0|＜3*Std0，即对于步骤S2-3-1中残差序列Val的每个元素dPi，统计不能满足条件1的通道个数k1。
步骤S2-3-5：判断步骤S2-3-4中k1是否大于等于门限Thsh1。如果满足，则进行步骤S2-3-6，否则，进行步骤S2-3-7。
所述的门限Thsh1计算方法如下：
Figure GDA0003698627160000131
式(22)中，m为步骤S2-3-2中残差序列Val的中间位置。
步骤S2-3-6：去除残差最大最小的通道后，直接统计其余通道残差的均值Mean和标准差Std。即对去除残差最大最小的通道后新的残差序列Val＝(dP2 ... dPi ... dPn-1)T,i＝2,...,n-1，直接统计其余通道残差的均值Mean和标准差Std，在计算标准差时，是使用步骤S2-3-2中的Mean0作为参考均值的。
步骤S2-3-7：设置变量初值j＝1，如果m-j>0，则进行步骤S2-3-8，否则，进行步骤S2-3-12。
步骤S2-3-8：从残差序列位置m-j处开始，到位置min(n-1,m+j)处结束，分别计算满足步骤S2-3-4中条件1|dPi-Mean0|＜3*Std0的各通道残差的平方和M与各通道残差减去残差均值Mean0的平方和S，并统计参与该计算的通道个数K2。
步骤S2-3-8中所述的残差均值Mean0和Std0初始值分别由步骤S2-3-2和步骤S2-3-3得到，后续循环计算采用的Mean0和Std0由步骤S2-3-10计算得到。
步骤S2-3-9：判断步骤S2-3-8中k2是否大于等于1。如果满足，则进行步骤S2-3-10，否则，进行步骤S2-3-11。
步骤S2-3-10：重新计算残差序列的均值Mean0和标准Std0。即有：
Figure GDA0003698627160000141
式(23)中，M、S、k2由步骤S2-3-8计算得到。
步骤S2-3-11：设置变量j＝j+1，跳转至步骤S2-3-7继续按流程顺序执行其后续步骤。
步骤S2-3-12：判断步骤S2-3-8中k2是否大于等于3。如果满足，则进行步骤S2-3-13，否则，结束流程。
步骤S2-3-13：重新计算残差序列的均值Mean和标准Std。计算方法同步骤S2-3-10。
如图4所示，是本发明实施例提供的判断先验轨道信息是否准确和接收机各通道是否含有伪距粗差的流程图，包括以下步骤：
步骤S3-1：设置静态变量初值ManFlag＝Num1＝Num2＝0。其中ManFlag作为先验轨道是否准确的标记，Num1和Num2分别为连续判断先验轨道是否准确的计数值。
步骤S3-2：判断步骤S2-7中的Std_min是否大于伪距测量噪声CodeNoise，如果是，则进行步骤S3-3，否则，进行步骤S3-7。所述额的距观测值的噪声大小CodeNoise可根据转移轨道航天器的轨道高度分段自适应调整，且伪距观测值的噪声阈值设置是可通过地面站参数上注的方式进行修改。
步骤S3-3：Num1计数值增加1。
步骤S3-4：判断Num1是否大于50，如果是，则进行步骤S3-5，否则，进行步骤S3-7。
步骤S2-5：设置先验轨道标记不准确标记ManFlag＝1。
步骤S3-6：放大补偿加速度的过程噪声到Sigma_M，设置Num1＝0。一般地，可设置Sigma_M＝1×10-5，所述的补偿加速度是式(6)中的
Figure GDA0003698627160000151
Figure GDA0003698627160000152
采用采用一阶高斯-马尔可夫随机模型，Sigma_M为所述的一阶高斯-马尔可夫随机模型的过程噪声值。
步骤S3-7：Num2计数值增加1。
步骤S3-8：判断Num2是否大于50，如果是，则进行步骤S3-9，否则，进行步骤S3-11。
步骤S3-9：设置先验轨道标记不准确标记ManFlag＝0。
步骤S3-10：恢复补偿加速度的过程噪声到Sigma_N,设置Num2＝0。一般地，可设置Sigma_N＝1×10-8。
步骤S3-11：按导航系统类型的不同，分别将不满足条件2：|dPi-Mean|＜3*Std的GPS和BDS的各通道标记为含有伪距粗差。所述的条件2中的dPi为步骤S2-2中算残差序列GVal和BVal的元素；所述的条件2中的Mean是步骤S2-3中计算的残差序列GVal和BVal的均值GMean和Bmean，由按导航系统类型决定取值；所述的条件2中的Std由步骤S2-6中计算得到的。
如图5所示，是本发明实施例提供的半实物仿真测试验证方案示意图，步骤如下：
(1)在GNSS信号模拟器的仿真控制软件中设置好转移轨道航天器的仿真测试场景，仿真控制软件驱动GNSS信号模拟器硬件，实时模拟生成转移轨道航天器在仿真参考轨道上能够接收到的GNSS信号射频信号，同时将NEMA信息输出给地检设备。
(2)地检设备接收GNSS信号模拟器输出的NEMA信息，然后解析NEMA信息中的时间信息，在预先设计并存储在本地磁盘的姿态、电推力文件中查找与仿真测试场景时间同步的姿态、电推力广播数据，然后模拟在轨卫星星务系统，将卫星姿态数据、电推力广播数据通过CAN总线播发给高灵敏度GNSS接收机。
(3)高灵敏度GNSS接收机实际接收GNSS信号模拟器输出GNSS信号射频信号和地检设备播发的卫星姿态数据和电推力广播数据。在高灵敏度GNSS接收机嵌入式处理器中，按本发明方法进行星载GNSS在轨实时定轨计算。
(4)最后，将GNSS接收机输出的实时定轨结果与仿真参考轨道进行对比验证，可完成本发明方法的定轨结果的精度评估。
最后应说明的是：以上实施例仅用以说明本发明的技术方案而非对其进行限制，尽管参照较佳实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对本发明的技术方案进行修改或者等同替换，而这些修改或者等同替换亦不能使修改后的技术方案脱离本发明技术方案的精神和范围。

Claims (3)

1.一种电推进转移轨道航天器星载GNSS在轨实时定轨方法，其特征在于，包括以下具体步骤：
步骤S1：实时获取GNSS观测数据、卫星姿态以及电推力广播数据；
步骤S2：当判断先验轨道信息有效标记有效时，计算并统计验前伪距观测值残差的均值和标准差，当判断先验轨道信息有效标记无效时，跳转到步骤S4；
步骤S3：根据步骤S2中得到的伪距观测值残差的均值和标准差信息判断先验轨道信息是否准确和接收机各通道是否含有伪距粗差，当先验轨道信息不准确时，调整补偿加速度的过程噪声；
步骤S4：完成顾及伪距粗差的伪距单点定轨,
星载GPS/BDS伪距观测方程表示为：
Figure FDA0003698627150000011
式(1)中，
Figure FDA0003698627150000012
分别为经过电离层改正后的星载GPS/BDS的伪距观测值；δtG和δtC分别为GPS/BDS接收机钟差；c为真空中的光速；
Figure FDA0003698627150000013
和
Figure FDA0003698627150000014
分别为第i颗GPS卫星和第j颗BDS卫星的钟差，分别由GPS/BDS广播星历计算得到；
Figure FDA0003698627150000015
和
Figure FDA0003698627150000016
分别为单频GPS/BDS电离层延迟改正值，
Figure FDA0003698627150000017
和
Figure FDA0003698627150000018
分别为GPS/BDS伪距观测值噪声；
Figure FDA0003698627150000019
为相应伪距观测值对应的站星几何距离，用下式表示：
Figure FDA00036986271500000110
式(2)中，(x,y,z)分别为GNSS接收机在地心地固系内三维位置坐标向量，
Figure FDA00036986271500000111
分别为第i颗GPS卫星和第j颗BDS卫星在地心地固系内三维位置坐标向量，分别由GPS/BDS广播星历计算得到；
假定GNSS接收机的位置坐标(x,y,z)和接收机钟差(cδtG,cδtC)组成的向量为X＝(x,y,z,cδtG,cδtC)T，其初值和改正值分别为
Figure FDA0003698627150000021
和ΔX＝(dx,dy,dz,dcδtG,dcδtC)T，则GPS/BDS伪距观测值
Figure FDA0003698627150000022
关于X的偏导数为：
Figure FDA0003698627150000023
如果t时刻星载GPS/BDS接收机共观测到m颗GPS和n颗BDS卫星，假设步骤S3中第k颗GPS卫星和第l颗BDS卫星被标识含有伪距粗差，则这两颗卫星不参与计算观测值矩阵A和残差向量b，其计算公式如下：
Figure FDA0003698627150000024
式(4)中，
Figure FDA0003698627150000025
和
Figure FDA0003698627150000026
分别为第i颗GPS卫星和第j颗BDS卫星伪距观测值
Figure FDA0003698627150000027
关于X的偏导数，初始计算时，是将X＝X0带入式(3)；
Figure FDA0003698627150000028
分别为GPS卫星实测的伪距观测值和伪距观测值的计算值；
Figure FDA0003698627150000029
分别为BDS卫星实测的伪距观测值和伪距观测值的计算值；
步骤S5：完成顾及电推力影响的定轨时间更新过程，在地心惯性系中，卫星运动方程用一阶微分方程组来表示：
Figure FDA00036986271500000210
式(5)中，
Figure FDA00036986271500000211
分别为卫星的位置、速度；
Figure FDA00036986271500000212
为卫星受到的各种摄动加速度之和，可以表示为：
Figure FDA00036986271500000213
式(6)中，
Figure FDA0003698627150000031
为近地航天器所受的惯性系总加速度；
Figure FDA0003698627150000032
为由保守力引起的加速度；
Figure FDA0003698627150000033
为由非保守力引起的加速度；
Figure FDA0003698627150000034
为电推力引起的加速度；
Figure FDA0003698627150000035
为人为引入的经验补偿加速度，用于补偿无法模型化或错误模型的微小摄动力的影响，采用一阶高斯-马尔可夫随机模型对径向、切向、法向3个方向进行动力学模型补偿；T是RTN坐标系到地心惯性系的转换矩阵；
由电推力引起的加速度
Figure FDA0003698627150000036
计算如下式所示：
Figure FDA0003698627150000037
式(7)中，
Figure FDA0003698627150000038
为卫星本体系到地心惯性系的转换矩阵，由步骤S1实时获取的姿态广播数据计算；
Figure FDA0003698627150000039
为电推力参考坐标系到卫星本体系的转换矩阵，由离子电推进器在卫星本体系具体安装方向计算得到；
Figure FDA00036986271500000310
为离子电推进器在三个方向上的电推力矢量，单位为N，由步骤S1实时获取的电推力广播数据得到；在精密定轨中，除了要考虑位置
Figure FDA00036986271500000311
参数和速度
Figure FDA00036986271500000312
参数之外，还要考虑动力学参数
Figure FDA00036986271500000313
且有
Figure FDA00036986271500000314
因此将式(5)进行扩展，并令
Figure FDA00036986271500000315
则有
Figure FDA00036986271500000316
当设定一初始状态
Figure FDA00036986271500000317
后，通过运动方程积分得到卫星的参考状态向量
Figure FDA00036986271500000318
该初值从步骤S4中得到星载GNSS接收机位置作为粗略轨道获得，由于上面的状态方程是非线性的，所以对
Figure FDA00036986271500000319
在
Figure FDA00036986271500000320
上展开取一次项，可得：
Figure FDA00036986271500000321
令
Figure FDA0003698627150000041
则有：
Figure FDA0003698627150000042
其解写为：
Figure FDA0003698627150000043
则有：
Figure FDA0003698627150000044
式(7)称为变分方程，其中I为单位矩阵，Φ(t,t0)被称为状态转移矩阵；
动力学定轨中，通过数值积分方法，按式(8)和式(11)，同时求解卫星的运动方程和变分方程；
伪距定轨卡尔曼滤波的状态方程为：
Xk＝Φk,k-1Xk-1+Wk (12)
式(5)中,Xk＝(r,v,cδtG,cδtC,Cd,Cr,ω)T为伪距定轨卡尔曼滤波状态量，Φk,k-1为离散形式的状态转移矩阵，Wk为系统噪声矩阵；
Figure FDA0003698627150000045
式(13)中，Φrr,Φrv,Φvr,Φvv分别位置和速度分量的状态转移矩阵；
Figure FDA0003698627150000046
位置和速度分别关于大气阻力和太阳光压的状态转移矩阵；Φrw,Φvw分别位置和速度关于RTN方向的补偿加速度的状态转移矩阵，Φww为补偿加速度关于其自身的状态转移矩阵；
伪距定轨卡尔曼滤波的观测方程为：
Zk＝HkXk-1+Vk (14)
式(14)中，Zk为观测矢量，Hk为观测矩阵，Vk为观测噪声矢量，且系统噪声Wk和观测噪声Vk，为零均值白噪声序列，Qk和Rk分别为系统动态噪声协方差阵和观测噪声协方差阵；
伪距定轨时间更新具体过程是，通过数值积分方法，同时求解卫星的运动方程和变分方程，得到Φk,k-1，再按式(12)完成对滤波状态量的时间更新，再完成对状态误差协方差矩阵的时间更新；
Figure FDA0003698627150000051
式(15)中，
Figure FDA0003698627150000052
为k-1时刻的状态误差协方差矩阵，
Figure FDA0003698627150000053
为k时刻的状态误差协方差矩阵；
步骤S6：完成顾及伪距粗差的定轨测量更新过程，式(3)中伪距观测值关于伪距定轨滤波器状态量的偏导数H可以表示为：
Figure FDA0003698627150000054
式(16)中，UT为J2000惯性系到地固系的转换矩阵，
Figure FDA0003698627150000055
为航天器相对于GPS卫星i和BDS卫星j的视线向量；
顾及伪距粗差的定轨测量更新过程为：按式(17)，依次对于步骤S3中没有被标识含有伪距粗差的GPS/BDS伪距观测值进行处理，更新绝对定轨卡尔曼滤波的状态量及状态协方差阵；
Figure FDA0003698627150000056
式(16)中，
Figure FDA0003698627150000061
为当前k时刻第i颗GPS/BDS对应的卡尔曼增益矩阵；
Figure FDA0003698627150000062
为式(16)中的观测矩阵；
Figure FDA0003698627150000063
为观测噪声协方差阵，
Figure FDA0003698627150000064
为第i颗GPS/BDS测量更新后的滤波状态量；
Figure FDA0003698627150000065
为步骤S5中当前k时刻时间更新后的滤波状态量，
Figure FDA0003698627150000066
第i颗GPS/BDS实测的伪距观测值；
Figure FDA0003698627150000067
为式(15)中当前k时刻时间更新后的状态误差协方差矩阵；
Figure FDA0003698627150000068
为当前k时刻测量更新后的状态误差协方差矩阵；
步骤S7：完成先验轨道的预报及先验轨道精度计算，标识先验轨道信息有效，在步骤S6中定轨测量更新最后得到的滤波状态量
Figure FDA0003698627150000069
的基础上，通过数值积分方法，按式(10)求解卫星运动方法，得到下一个时刻的先验轨道PreOrb，先验轨道精度Accu由在步骤S6中状态误差协方差矩阵
Figure FDA00036986271500000610
计算；当步骤S6中参与测量更新的卫星数大于0时，标识先验轨道信息有效，否则，标识先验轨道信息无效；
步骤S8：输出定轨滤波结果和先验轨道预报结果。
2.根据权利要求1所述的一种电推进转移轨道航天器星载GNSS在轨实时定轨方法，其特征在于：步骤S2中，利用步骤S7中计算得到的先验轨道预报结果代入GNSS伪距观测方程来计算伪距观测值的验前残差，并根据导航系统类型的不同分别计算每个导航系统的伪距残差序列的均值和标准差，然后计算所有通道的伪距残差序列的均值和标准差。
3.根据权利要求1所述的一种电推进转移轨道航天器星载GNSS在轨实时定轨方法，其特征在于：在步骤S3中，判断先验轨道信息是否准确的方法如下：
当连续多次满足伪距验前残差标准差大于伪距观测值的噪声时，则判断先验轨道信息是准确的，否则，不准确；
其中伪距观测值的噪声大小根据转移轨道航天器的轨道高度分段自适应调整，且伪距观测值的噪声阈值设置是通过地面站参数上注的方式进行修改。
CN202111456809.1A 2021-12-02 2021-12-02 电推进转移轨道航天器星载gnss在轨实时定轨方法 Active CN114063122B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111456809.1A CN114063122B (zh) 2021-12-02 2021-12-02 电推进转移轨道航天器星载gnss在轨实时定轨方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111456809.1A CN114063122B (zh) 2021-12-02 2021-12-02 电推进转移轨道航天器星载gnss在轨实时定轨方法

Publications (2)

Publication Number Publication Date
CN114063122A CN114063122A (zh) 2022-02-18
CN114063122B true CN114063122B (zh) 2022-08-02

Family

ID=80228280

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111456809.1A Active CN114063122B (zh) 2021-12-02 2021-12-02 电推进转移轨道航天器星载gnss在轨实时定轨方法

Country Status (1)

Country Link
CN (1) CN114063122B (zh)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114861320B (zh) * 2022-05-19 2023-02-10 北京航天飞行控制中心 一种航天器姿控推力建模及定轨解算方法

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1037404A2 (en) * 1999-03-16 2000-09-20 Hitachi, Ltd. Elliptical satellite communication system
EP1076005A2 (en) * 1999-08-13 2001-02-14 Hughes Electronics Corporation Spacecraft orbit control using orbit position feedback
CN101435863A (zh) * 2008-12-25 2009-05-20 武汉大学 一种导航卫星实时精密定轨的方法
CN101872021A (zh) * 2010-05-20 2010-10-27 武汉大学 Gps双频实时星载数据处理方法
CN101893712A (zh) * 2010-07-09 2010-11-24 中国科学院测量与地球物理研究所 用于地球静止卫星精密定轨的选权拟合方法
CN101968361A (zh) * 2009-07-28 2011-02-09 韩春好 基于星光观测的空间绝对定向技术
CN105651516A (zh) * 2014-11-11 2016-06-08 航天恒星科技有限公司 基于gnss观测值的发动机推力标定方法及装置
CN112083452A (zh) * 2020-08-28 2020-12-15 中国人民解放军63921部队 卫星导航接收机的环路跟踪系统及方法

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1037404A2 (en) * 1999-03-16 2000-09-20 Hitachi, Ltd. Elliptical satellite communication system
EP1076005A2 (en) * 1999-08-13 2001-02-14 Hughes Electronics Corporation Spacecraft orbit control using orbit position feedback
CN101435863A (zh) * 2008-12-25 2009-05-20 武汉大学 一种导航卫星实时精密定轨的方法
CN101968361A (zh) * 2009-07-28 2011-02-09 韩春好 基于星光观测的空间绝对定向技术
CN101872021A (zh) * 2010-05-20 2010-10-27 武汉大学 Gps双频实时星载数据处理方法
CN101893712A (zh) * 2010-07-09 2010-11-24 中国科学院测量与地球物理研究所 用于地球静止卫星精密定轨的选权拟合方法
CN105651516A (zh) * 2014-11-11 2016-06-08 航天恒星科技有限公司 基于gnss观测值的发动机推力标定方法及装置
CN112083452A (zh) * 2020-08-28 2020-12-15 中国人民解放军63921部队 卫星导航接收机的环路跟踪系统及方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
"高轨航天器GNSS技术发展";王猛等;《测绘学报》;20200930;第49卷(第9期);第1158-1167页 *

Also Published As

Publication number Publication date
CN114063122A (zh) 2022-02-18

Similar Documents

Publication Publication Date Title
US6735523B1 (en) Process and system of coupled real-time GPS/IMU simulation with differential GPS
CN110764127B (zh) 易于星载在轨实时处理的编队卫星相对定轨方法
CN108594271B (zh) 一种基于复合分层滤波的抗欺骗干扰的组合导航方法
CN111965685B (zh) 一种基于多普勒信息的低轨卫星/惯性组合导航定位方法
US10908296B2 (en) Method for calculating a speed of an aircraft, method for calculating a protection radius, positioning system and associated aircraft
CN101395443A (zh) 混合定位方法和设备
CN110673175A (zh) 一种基于gnss广播星历的高轨卫星高精度自主定轨方法
CN114063122B (zh) 电推进转移轨道航天器星载gnss在轨实时定轨方法
Mahmoud et al. Integrated INS/GPS navigation system
Langel et al. Tightly coupled GPS/INS integration for differential carrier phase navigation systems using decentralized estimation
Biswas et al. Computationally efficient unscented Kalman filtering techniques for launch vehicle navigation using a space-borne GPS receiver
Carpenter et al. Semimajor axis knowledge and GPS orbit determination
CN110307840B (zh) 一种基于多波束测距测速和惯性的着陆段鲁棒融合方法
Vana et al. Benefits of motion constraining for robust, low-cost, dual-frequency GNSS PPP+ MEMS IMU navigation
Hou et al. A dual-satellite GNSS positioning algorithm of high accuracy in incomplete condition
CN116299599A (zh) 一种ins辅助的gnss伪距粗差探测方法
CN114964215A (zh) 多个探测目标主体的轨道确定方法及装置、电子设备
CN111538045A (zh) 一种星载导航接收机在轨精度预先评估方法
CN112799105A (zh) 一种编队leo卫星星间时间同步和评估方法
CN113866732B (zh) 一种单部雷达短弧测轨能力的计算方法
Ismaeel Design of Kalman Filter of Augmenting GPS to INS Systems
Bolandi et al. GPS based onboard orbit determination system providing fault management features for a LEO satellite
Upadhyay et al. Precision relative navigation for automated rendezvous and docking
Langel et al. Cycle ambiguity reacquisition in UAV applications using a novel GPS/INS integration algorithm
CN111238484B (zh) 一种基于球形无迹变换的环火轨道自主导航方法

Legal Events

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