CN110764127B - 易于星载在轨实时处理的编队卫星相对定轨方法 - Google Patents

易于星载在轨实时处理的编队卫星相对定轨方法 Download PDF

Info

Publication number
CN110764127B
CN110764127B CN201910949615.1A CN201910949615A CN110764127B CN 110764127 B CN110764127 B CN 110764127B CN 201910949615 A CN201910949615 A CN 201910949615A CN 110764127 B CN110764127 B CN 110764127B
Authority
CN
China
Prior art keywords
satellite
orbit
relative
orbit determination
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
CN201910949615.1A
Other languages
English (en)
Other versions
CN110764127A (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
Original Assignee
Wuhan University WHU
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 filed Critical Wuhan University WHU
Priority to CN201910949615.1A priority Critical patent/CN110764127B/zh
Publication of CN110764127A publication Critical patent/CN110764127A/zh
Application granted granted Critical
Publication of CN110764127B publication Critical patent/CN110764127B/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/38—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
    • G01S19/39—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/42—Determining position
    • G01S19/51—Relative positioning
    • 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/38—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
    • G01S19/39—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/42—Determining position
    • G01S19/50—Determining position whereby the position solution is constrained to lie upon a particular curve or surface, e.g. for locomotives on railway tracks

Landscapes

  • Engineering & Computer Science (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

本发明公开了一种易于星载在轨实时处理的编队卫星相对定轨方法，属于卫星编队相对导航领域，该方法使用GPS/BDS广播星历和简化的相对动力学模型，对编队卫星GPS/BDS双频观测数据进行在轨实时处理，考虑星载处理器计算资源有限，设计了动力学轨道预报配合轨道内插的运行模式，并将核心复杂数据处理过程通过分时多步实现来降低单历元的计算耗时，可输出频率可调的编队卫星高精度绝对与相对轨道结果。本发明不仅适用于较大星间距离范围的星间编队，在GPS/BDS信号中断情况下，仍可连续输出满足一定精度要求较为平滑的导航结果，而且具有定轨精度高、稳定性好及对星载处理器性能要求较低，易于工程化实现等特点。

Description

易于星载在轨实时处理的编队卫星相对定轨方法
技术领域
本发明属于卫星编队相对导航技术领域，更具体地，涉及一种易于星载在轨实时处理的编队卫星GPS/BDS相对定轨方法。
背景技术
卫星编队飞行是近年来国内外航天界一个非常热点的研究领域，已广泛应用于分布式雷达干涉、时变地球重力场测量、海洋环境监测及航天军事侦察等科学任务中，比较具有代表的应用包括国外的GRACE-A/B、TerraSAR/TanDEM、SWARM-A/B/C和国内的实践九号A/B、探索四号A/B、环境一号A/B/C等。
卫星及其编队卫星的高精度绝对轨道及相对基线确定，是完成编队卫星飞行任务的关键，星载GNSS因其具有可连续观测、精度高、成本功耗低及体积小重量轻等优点，已成为卫星及其编队绝对定轨与相对定位的主要技术手段。通过地面事后处理的方式，低轨卫星单星绝对定轨精度已达到cm级，卫星编队星间基线解算精度已达到mm级。随着卫星应用技术的快速发展，各类科学任务与应用在精度、实时性和可靠性方面，对卫星及其编队卫星的轨道提出了越来越高的需求，因此通过星载在轨实时处理完成卫星轨道及基线确定已成为未来发展趋势。
然而，卫星上搭载的嵌入式处理器无法与地面PC机相比拟，特别是微小卫星，计算能力和资源相对有限，难以完成在轨实时相对定轨庞大的计算量，因此，迫切地需要一种较易适用于星载在轨实时处理的方法。
发明内容
针对现有技术的以上缺陷或改进需求，本发明提出了一种易于星载在轨实时处理的编队卫星GPS/BDS相对定轨方法，由此解决现有适用于在轨工程化实际应用的星载实时编队卫星相对定轨方式无法同时兼顾高精度和实时性的技术问题。
本发明提供了一种易于星载在轨实时处理的编队卫星相对定轨方法，包括：
(1)在当前历元时刻的星载GPS/BDS实时相对定轨开始后，判断是否需要进行初始化，若需要进行初始化，则进行初始化操作，若不需要进行初始化，则进入步骤(2)；
(2)获取定轨实时数据，根据所述定轨实时数据得到双频无电离层伪距和相位LC组合观测值、电离层残差GF组合观测值及MW组合观测值；
(3)根据所述电离层残差GF组合观测值及所述MW组合观测值对非差相位数据进行实时周跳探测，并标记存在周跳的非差相位观测值；
(4)判断定轨滤波器启动标记是否为已启动，若所述启动标记为已启动，则进入步骤(5)，若所述启动标记为未启动，则需完成定轨滤波的启动，并将所述启动标记设置为已启动，并结束当前历元任务；
(5)由AB星轨道预报得到的轨道插值首尾端点信息，使用Hermit插值分别得到AB星在当前历元整秒时刻的单星绝对位置速度结果及相对位置速度结果；
(6)由当前观测历元时间与插值首端点信息时间的间隔判断是否达到定轨滤波时刻，若已达到定轨滤波时刻，则进入步骤(7)，若未达到定轨滤波时刻，则结束当前历元任务；
(7)基于不存在周跳的非差/单差相位观测值、所述AB星在当前历元整秒时刻的单星绝对位置速度结果及相对位置速度结果进行定轨滤波解算；
(8)重复步骤(1)～步骤(7)，通过在轨实时完成编队卫星GPS/BDS相对定轨。
优选地，步骤(4)中的定轨滤波启动包括如下子步骤：
(4.1)使用编队AB卫星星载GPS/BDS双频伪距及多普勒观测数据分别进行标准单点定位与单点测速，分别得到AB卫星天线相位中心在地心地固系下的位置和速度、接收机钟差和钟漂参数及定位测速中伪距多普勒观测值残差的RMS值；
(4.2)组成AB星站间伪距单差观测值，进行伪距单差相对定位，得到AB卫星天线相位中心之间基线向量在地心地固系下的相对位置和相对钟差；
(4.3)若不满足定轨滤波启动条件，则直接结束当前历元任务，若满足启动条件，则执行步骤(4.4)；
(4.4)将AB卫星天线相位中心在地心地固系下的位置及AB卫星天线相位中心之间基线向量在地心地固系下的相对位置，从地心地固系转换到地心惯性系，再初始化绝对与相对定轨滤波器的滤波状态量及状态误差协方差阵；
(4.5)考虑AB卫星姿态，分别将AB星在地心惯性系下的天线相位中心转到AB卫星质心，以AB卫星质心为轨道积分预报初值，使用AB星动力学模型参数，分别对A星和B星进行动力学轨道预报，得到轨道插值首尾端点信息，并将定轨滤波器启动标记设置为已启动，结束当前历元任务。
优选地，步骤(7)包括：
(7.1)使用星载GPS/BDS双频伪距和多普勒观测数据对编队AB卫星分别进行标准单点定位，得到AB卫星天线相位中心在地心地固系下的位置和AB星载接收机钟差参数；
(7.2)根据AB卫星天线相位中心在地心地固系下的位置和AB星载接收机钟差参数，完成A星绝对定轨时间更新及A星绝对定轨测量更新，依次对每颗没有标记粗差的GPS/BDS非差载波相位观测值进行处理，更新卡尔曼滤波的状态及协方差参数；
(7.3)生成AB星站间单差数据，并根据AB星站间单差数据得到单差数据的电离层残差GF组合观测值和MW组合观测值，根据单差数据的电离层残差GF组合观测值和MW组合观测值完成对单差数据的周跳探测，标记存在周跳的相位单差观测值，其中，存在周跳的相位单差观测值不参与后续的相对定轨测量更新；
(7.4)完成AB星相对定轨时间更新，采用基于相对动力学先验轨道来探测单差相位观测值的验前粗差，标记出存在验前粗差的单差数据，其中，存在验前粗差的单差数据不参与后续的相对定轨测量更新；
(7.5)完成AB星相对定轨测量更新，依次对每颗没有标记粗差的GPS/BDS单差载波相位观测值进行处理，更新卡尔曼滤波的状态及协方差参数；
(7.6)考虑AB卫星姿态，分别将AB星测量更新的结果转到AB卫星质心，以AB卫星质心结果为轨道积分预报初值，分别对A星和B星进行动力学轨道预报，得到轨道插值首尾端点信息。
优选地，在步骤(7.2)中，所述完成A星绝对定轨时间更新包括：
考虑A星卫星姿态，先将上个历元A星测量更新的结果转到A卫星质心，再使用A星动力学模型参数，经动力学轨道积分得到的A星卫星状态及状态转移矩阵，然后再考虑A星卫星姿态，将A星质心的结果转到A星天线相位中心，最后将绝对定轨滤波的状态和状态误差协方差矩阵更新到当前测量更新时刻。
优选地，在步骤(7.4)中，所述完成AB星相对定轨时间更新包括：
考虑AB星卫星姿态，利用上一历元A星绝对定轨测量更后的结果和相对定轨测量更新结果，得到B星绝对位置速度结果，将B星天线相位中心转到B卫星的质心处，再使用B星动力学模型参数，经动力学轨道积分得到的B星卫星状态及状态转移矩阵，再将B星质心的结果转到B星天线相位中心，然后减去A星时间更新后的结果，得到AB星相对轨道状态及相对状态转移矩阵，最后将相对定轨滤波的状态和状态误差协方差矩阵更新到当前测量更新时刻。
优选地，在步骤(7.2)和步骤(7.4)中的动力学轨道积分过程，均使用经典4阶龙格库塔法对卫星运动方程和状态转移矩阵微分方程进行数值积分。
优选地，在步骤(4.5)和步骤(7.6)中的AB星动力学轨道积分预报中使用的重力场模型阶数均为30阶，在步骤(7.2)中的A星绝对定轨时间更新和步骤(7.4)中的AB星相对定轨时间更新中使用的重力场模型阶数均为70阶。
总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列有益效果：
(1)在步骤(4)和步骤(7)中完成AB星轨道预报，生成轨道插值首尾端点信息，然后配合步骤(5)，通过轨道内插即可输出频率可调的编队卫星高精度绝对与相对轨道结果。同时上述模式下，将步骤(4)和步骤(7)中的多个子步骤的计算量，通过分时多步计算来完成，以降低步骤(4)和步骤(7)单历元的运算耗时，使算法更易于在星载嵌入式处理器上实时处理。
(2)在动力学轨道积分预报(步骤(4.5)和步骤(7.6))和时间更新(步骤(7.2)和步骤(7.4))中，为了提高计算精度，均采用经典4阶龙格－库塔法RK4(或RKF5)数值积分方法分别计算卫星轨道状态和状态转移矩阵，其中轨道积分预报时长h＝45s(可通过上注方式修改)，时间更新间隔step＝30s(可通过上注方式修改)，且准确的时间更新间隔考虑了AB星星载接收机钟差参数(δtA,δtB)，即步骤(7.2)中绝对定轨时间更新间隔
Figure BDA0002225353780000051
或
Figure BDA0002225353780000052
步骤(7.4)中相对定轨时间更新间隔
Figure BDA0002225353780000053
或
Figure BDA0002225353780000054
(3)本发明适用于较大星间距离范围的星间编队，相对于传统RTK方法，在GPS/BDS信号中断情况下，仍可连续输出高精度导航结果；
(4)本发明具有定轨精度高、稳定性好、对星载处理器要求不高，易于工程化实现等特点；
(5)本发明适用于GPS L1/L2和北斗二代B1I/B2I或B1I/B3I双频信号，同时可推广到北斗三代B1c/B2a或GLONASS/GALILEO等。
附图说明
图1是本发明实施例提供的一种易于星载在轨实时处理的编队卫星GPS/BDS相对定轨方法的流程示意图；
图2是本发明实施例提供的一种短基线双差固定解方法的流程示意图；
图3是本发明实施例提供的一种GRACE编队卫星相对定轨结果图；
图4是本发明实施例提供的一种轨道预报配合轨道内插模式下的定轨滤波启动和定轨滤波解算过程分时分步计算的流程示意图；
图5是本发明实施例提供的一种1s计算10步定轨滤波解算过程的计算时序图。
具体实施方式
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，但并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
如图1所示是本发明实施例提供的一种易于星载在轨实时处理的编队卫星GPS/BDS相对定轨方法的流程示意图，包括以下步骤：
步骤S1：开始星载GPS/BDS实时相对定轨当前历元时刻计算任务。
步骤S2：判断是否需要进行系统初始化。
若不需要则直接进行步骤S3，若需要则进行系统初始化，系统初始化包括：时间更新时间间隔(step)、轨道预报间隔(h)、卫星星体参数(整星质量、接收机天线相位中心偏差PCO参数)、地球重力场、地球定向参数(EOP)、跳秒值(GPST-UTC)、定轨相关阈值门限参数、定轨相关的全局变量及结构体等的初始化；
步骤S3：定轨输入数据获取及数据预处理。
定轨输入数据包括编队卫星AB(一般地，A为主星，B为从星，下同)的星载GPS/BDS双频观测数据(伪距、载波相位、多普勒)、广播星历数据及卫星姿态数据等。数据预处理包括计算双频无电离层伪距和相位无电离层LC组合观测值、电离层残差GF组合观测值及MW组合观测值；
其中，无电离层相位和伪距LC组合观测值可表示为：
Figure BDA0002225353780000071
电离层残差GF组合观测值可表示为：
Lgf＝L1-L2 (2)
MW组合观测值可表示为：
Figure BDA0002225353780000072
其中，f1、f2为GPS/BDS双频观测数据对应的信号频点，LC和PC分别表示GPS/BDS双频无电离层相位和伪距的组合观测值，L1和L2分别表示GPS L1/BDS B1和GPS L2/BDS B2(或B3)频点的相位观测值，P1和P2分别表示GPS L1/BDS B1和GPS L2/BDS B2(或B3)频点的伪距观测值，Lgf表示电离层残差GF组合观测值，Lmw表示MW组合观测值。
步骤S4：星载GPS/BDS双频非差数据周跳探测。
根据GF和MW组合观测值对非差相位数据进行实时周跳探测，标记存在周跳的非差相位观测值，不参与后续的绝对定轨测量更新。
具体方法是：首先根据式(4)计算历元间GF和MW组合观测值变化量dLgf和dLmw，当同时满足dLgf＜GFThreshold和dLmw＜MWThreshold时，将该非差相位数据标记为1，即认为不存在周跳，否则标记为-2，即认为存在周跳，后续不参与绝对定轨滤波解算。GFThreshold和MWThreshold为周跳探测的具体阈值门限，可通过上注方式修改。
Figure BDA0002225353780000081
其中，Lgf(t)、Lgf(t-1)分别表示当前/上一时刻的电离层残差GF组合观测值，Lmw(t)、Lmw(t-1)分别表示当前/上一时刻的MW组合观测值。
步骤S5：判断定轨滤波器是否已经启动。
即判断定轨滤波器启动标记是否为已启动，如果该标记为已启动，则直接进入步骤S7，否则进入步骤S6。
步骤S6：执行完成定轨滤波的启动；
步骤S6-1：AB星单点定位测速
使用编队AB卫星星载GPS/BDS双频伪距、多普勒观测数据分别进行标准单点定位(SPP)与单点测速(SPV)，得到AB卫星天线相位中心在地心地固系(WGS84)下的位置
Figure BDA0002225353780000082
速度
Figure BDA0002225353780000083
和接收机钟差
Figure BDA0002225353780000084
和
Figure BDA0002225353780000085
钟漂
Figure BDA0002225353780000086
参数和定位测速中伪距多普勒残差的RMS值
Figure BDA0002225353780000087
其中
Figure BDA0002225353780000088
和
Figure BDA0002225353780000089
分别为AB卫星星载GPS/BDS接收机关于GPS和BDS导航系统的钟差。
步骤S6-2：AB星站间伪距单差相对定位
组成AB星站间伪距单差观测值，进行伪距单差相对定位，得到AB卫星天线相位中心之间基线向量在地心地固系(WGS84)下的相对位置
Figure BDA0002225353780000091
和相对钟差
Figure BDA0002225353780000092
其中，
Figure BDA0002225353780000093
基线向量的相对速度
Figure BDA0002225353780000094
步骤S6-3：判断是否满足定轨滤波启动条件
如果不满足定轨滤波启动条件，则直接结束当前历元所有计算任务(步骤S10)，如果满足启动条件，则进行步骤S6-4。启动条件为同时满足
Figure BDA0002225353780000095
其中σp和σv是参考阈值门限，可通过上注方式修改。
步骤S6-4：绝对与相对定轨滤波器初始化
将S6-1和S6-2计算的结果
Figure BDA0002225353780000096
和
Figure BDA0002225353780000097
从地心地固系(WGS84)转换到地心惯性系(J2000)，得到在地心惯性系下结果
Figure BDA0002225353780000098
和
Figure BDA0002225353780000099
再初始化绝对与相对定轨滤波器的滤波状态量(XA,XAB)及状态误差协方差阵(PA,PAB)等，其中：
Figure BDA00022253537800000910
滤波状态量XA,XAB初始值设置如下：
Figure BDA00022253537800000911
Cd、dCd和Cr、dCr分别为待估的绝对与相对大气阻力系数和太阳光压系数，W和dW分别为RTN三个方向待估的动力学模型补偿加速度参数，NG和NC分别为GPS和BDS对应12通道的非差相位模糊度参数(单位为m)，P和L分别为由式(1)计算的双频无电离层伪距和相位LC组合观测值。
Figure BDA00022253537800000912
和
Figure BDA00022253537800000913
分别为GPS和BDS对应12通道的单差相位模糊度参数(单位为m)，PSD和LSD分别为AB星由对应的共视卫星P和L计算的站间单差伪距和相位观测值，上标G和C分别表示对应GPS和BDS系统的各变量。
状态误差协方差阵PA,PAB初始值可设置为：
PA＝PAB＝diag(1E4,1E2,1E4,1E4,10,10,0.01,SigN2,SigN2)
其中，SigN2为相位模糊度的初始方差，可通过上注方式修改。
步骤S6-5：AB星轨道预报
考虑AB卫星姿态，分别将AB星初始轨道状态
Figure BDA0002225353780000101
和
Figure BDA0002225353780000102
由天线相位中心转到AB卫星质心
Figure BDA0002225353780000103
和
Figure BDA0002225353780000104
其中
Figure BDA0002225353780000105
然后以AB卫星质心为轨道积分预报初值，使用AB星动力学模型参数，分别对A、B星进行动力学轨道预报，得到AB星轨道插值所需的首尾端点信息
Figure BDA0002225353780000106
和
Figure BDA0002225353780000107
其中，t,r,v,a分别为时间、位置、速度、加速度。
所述的动力学轨道预报过程，是使用经典4阶龙格－库塔法(RK4)，对考虑了地球引力、日月引力、固体潮、大气阻力及太阳光压等摄动力影响的卫星运动方程进行数值积分，该过程积分步长为h，使用的地球重力场模型阶数为30阶。
步骤S6-6：置标记为滤波器已启动
将定轨滤波器启动标记设置为已启动，结束当前历元所有计算任务(第10步)。
步骤S7：轨道内插输出。
即由定轨滤波启动步骤S6或定轨滤波解算步骤S9中AB星轨道预报得到的轨道插值首尾端点信息，使用5阶Hermit插值分别计算AB星在当前历元整秒时刻的单星绝对位置速度结果，再计算相对位置速度结果；
具体方法为：对于任一时刻t∈(t0,t1)的卫星状态(r(t),v(t))，可用一个五阶的Hermite多项式近似表示如下：
Figure BDA0002225353780000111
所述的式(7)中的位置、速度系数项计算方法如下所示：
Figure BDA0002225353780000112
由步骤S6-5或步骤S9-8计算得到的AB卫星轨道插值首尾端点信息结果
Figure BDA0002225353780000113
和
Figure BDA0002225353780000114
经式(7)计算AB星在当前历元整秒时刻t的单星绝对位置速度结果(rA(t),vA(t))和(rB(t),vB(t))，再计算AB星相对位置速度结果(rAB(t),vAB(t))，其中rAB(t)＝rB(t)-rA(t),vAB(t)＝vB(t)-vA(t)。
步骤S8：判断是否达到定轨滤波时刻。
即计算当前观测历元时间与插值前首端点信息时间的间隔dt，如果dt取整后大于等于卡尔曼滤波时间间隔step，则已达到定轨滤波时刻，需要进入步骤S9，否则进入步骤S10。
步骤S9：执行完成定轨滤波解算任务。
步骤S9-1：AB星单点定位
使用星载GPS/BDS双频伪距、多普勒观测数据对编队AB卫星分别进行标准单点定位(SPP)，得到AB卫星天线相位中心在地心地固系(WGS84)下的位置和AB星载接收机钟差参数
Figure BDA0002225353780000121
和
Figure BDA0002225353780000122
步骤S9-2：A星绝对定轨时间更新
设绝对定轨卡尔曼滤波的状态方程为：
Figure BDA0002225353780000123
式(8)中
Figure BDA0002225353780000124
为绝对定轨卡尔曼滤波状态量，详见步骤S6-4，
Figure BDA0002225353780000125
为状态转移矩阵，Wk为系统噪声矩阵。
Figure BDA0002225353780000126
Φrr,Φrv,Φvr,Φvv分别为位置和速度分量的状态转移矩阵，
Figure BDA0002225353780000127
表示位置和速度分别关于大气阻力和太阳光压的状态转移矩阵，Φrw,Φvw分别位置和速度关于RTN方向的补偿加速度的状态转移矩阵，Φww为补偿加速度关于其自身的状态转移矩阵，
Figure BDA0002225353780000128
为对应GPS/BDS通道非差模糊度关于其自身的状态转移矩阵，下表k表示当前时刻，k-1表示上一时刻。
绝对定轨卡尔曼滤波的观测方程为：
Figure BDA0002225353780000129
其中，
Figure BDA0002225353780000131
为观测矢量，
Figure BDA0002225353780000132
为观测矩阵，Vk为观测噪声矢量，且系统噪声Wk和观测噪声Vk，为零均值白噪声序列，即对所有的k,j，有
Figure BDA0002225353780000133
式中，Qk和Rk分别为系统动态噪声协方差阵和观测噪声协方差阵，δkj为Dirac-δ函数。
A星绝对定轨时间更新具体过程是，考虑A星卫星姿态，先将上个历元A星测量更新的结果
Figure BDA0002225353780000134
转到A卫星质心
Figure BDA0002225353780000135
以
Figure BDA0002225353780000136
为初始值，使用A星动力学模型参数，经动力学轨道积分得到的滤波测量更新时刻的A星卫星状态
Figure BDA0002225353780000137
及状态转移矩阵
Figure BDA0002225353780000138
然后再考虑A星卫星姿态，将积分结果A星质心的结果
Figure BDA0002225353780000139
转到A星天线相位中心
Figure BDA00022253537800001310
最后使用式(8)和式(12)将绝对定轨滤波的状态量和状态误差协方差矩阵更新到当前测量更新时刻；
Figure BDA00022253537800001311
其中，
Figure BDA00022253537800001312
表示当前时刻的时间更新后的状态误差协方差矩阵，
Figure BDA00022253537800001313
表示上一时刻的时间更新后的状态误差协方差矩阵。
所述的动力学轨道积分过程，是使用经典4阶龙格－库塔法(RK4)，同时对卫星运动方程和状态转移矩阵微分方程进行数值积分，该过程实际积分步长为
Figure BDA00022253537800001314
或
Figure BDA00022253537800001315
使用的地球重力场模型阶数为70阶。
步骤S9-3：A星绝对定轨测量更新
双频无电离层LC组合的非差相位观测值可以表示为：
Figure BDA0002225353780000141
其中，ρj表示GPS/BDS发射与接收天线相位中心之间的几何距离，c表示真空中的光速，δtr、
Figure BDA0002225353780000142
分别表示接收机和第j颗导航卫星的钟差，N表示不具有整数特性的非差模糊度，单位米，εP表示相位测量噪声。
式(13)中非差相位观测值关于绝对定轨滤波器状态量的偏导数H可以表示为：
Figure BDA0002225353780000143
式(14)中：UT为J2000惯性系到地固系的转换矩阵，
Figure BDA0002225353780000144
为A星相对于GPS卫星i和BDS卫星j的视线向量，I1×12中第i颗GPS或第j颗BDS对应的通道为1，其余通道为0，下标G和C分别表示GPS和BDS系统。
按式(15)，依次对每颗没有标记粗差的GPS/BDS非差载波相位观测值进行处理，更新绝对定轨卡尔曼滤波的状态量及状态协方差阵；
Figure BDA0002225353780000145
其中，
Figure BDA0002225353780000146
表示当前时刻卡尔曼增益矩阵，
Figure BDA0002225353780000147
表示式(14)中的第i颗对应的观测矩阵，
Figure BDA0002225353780000148
表示式(11)中的观测噪声协方差阵，
Figure BDA0002225353780000149
表示测量更新后的滤波状态量，
Figure BDA00022253537800001410
表示时间更新后的滤波状态量，
Figure BDA00022253537800001411
表示式(10)中的观测值，
Figure BDA00022253537800001412
表示当前时刻测量更新后的状态误差协方差矩阵。
步骤S9-4：生成AB星站间单差观测数据与观测值域周跳探测
首先按照式(16)生成AB星站间单差伪距和相位数据，计算单差数据的电离层残差GF组合和MW组合观测值。
Figure BDA0002225353780000151
式(16)中，
Figure BDA0002225353780000152
分别为AB星GPS/BDS伪距和相位无电离层LC组合观测值，P1 B,P1 A,
Figure BDA0002225353780000153
分别为AB星GPS/BDS导航系统对应的双频非差伪距和相位观测值，f1和f2为GPS/BDS双频信号频点，
Figure BDA0002225353780000154
和
Figure BDA0002225353780000155
分别表示AB星GPS/BDS双频伪距和相位无电离层组合计算的单差相位观测值。
根据单差数据的GF和MW组合观测值对单差数据进行实时周跳探测，标记存在周跳的观测值。即根据(17)式计算历元间单差数据的GF和MW组合观测值变化量
Figure BDA0002225353780000156
和
Figure BDA0002225353780000157
当同时满足
Figure BDA0002225353780000158
和
Figure BDA0002225353780000159
时，将该单差数据标记为1，即认为不存在周跳，否则标记为-2，即认为存在周跳，后续不参与相对定轨滤波解算。
Figure BDA00022253537800001510
其中，
Figure BDA00022253537800001511
分别表示当前/上一时刻的单差数据的GF组合观测值
Figure BDA00022253537800001512
分别表示当前/上一时刻的单差数据的MW组合观测值。
步骤S9-5：AB星相对定轨时间更新
设相对定轨卡尔曼滤波的状态方程和观测方程为：
Figure BDA00022253537800001513
式(18)中，
Figure BDA0002225353780000161
为相对定轨卡尔曼滤波状态量，详见步骤S6-4，此外可推导得出
Figure BDA0002225353780000162
Figure BDA0002225353780000163
表示GPS/BDS单差观测值，
Figure BDA0002225353780000164
表示观测矩阵。
因此，相对定轨时间更新可参考绝对定轨时间更新过程，具体为：考虑AB星卫星姿态，利用上一历元A星绝对定轨测量更新的结果
Figure BDA0002225353780000165
相对定轨测量更新结果
Figure BDA0002225353780000166
得到B星绝对位置速度结果
Figure BDA0002225353780000167
将该结果从天线相位中心转到B卫星的质心处
Figure BDA0002225353780000168
再使用B星动力学模型参数，经动力学轨道积分得到滤波更新时刻的B星卫星状态
Figure BDA0002225353780000169
及状态转移矩阵
Figure BDA00022253537800001610
再将B星质心的结果转到B星天线相位中心
Figure BDA00022253537800001611
然后减去A星时间更新后的结果
Figure BDA00022253537800001612
得到AB星相对轨道状态
Figure BDA00022253537800001613
最后式(18)和式(12)将相对定轨滤波的状态量和状态误差协方差矩阵更新到当前测量更新时刻，相对定轨时间更新过程实际积分步长为
Figure BDA00022253537800001614
或
Figure BDA00022253537800001615
使用的地球重力场模型阶数为70阶。
步骤S9-6：基于动力学先验轨道的单差验前粗差探测；
具体过程为：以时间更新的滤波状态量结果为先验轨道，首先计算单差相位数据的残差向量Val，使用Grubbs准则计算Val中元素的均值Mean和标准差Std，对于每颗卫星，当满足|Val-Mean|>3*Std时，标记该颗卫星为验前粗差，并不参与后续的相对定轨滤波解算。
步骤S9-7：相对定轨测量更新；
单差相位观测值可表示为：
LAB j＝LB j-LA j＝ρAB j+c(δtrAB-δtsAB j)+NSD+εPAB (19)
式(19)中的单差相位观测值，对于中长基线，使用的是式(16)中AB星GPS/BDS双频相位无电离层组合计算的单差相位观测值
Figure BDA0002225353780000171
对于短基线，使用的是式(16)中AB星GPSL1/BDS B1单频非差相位观测值计算的单差相位观测值
Figure BDA0002225353780000172
LAB j表示AB星GPS/BDS单差相位观测值，LA j、LB jAB星第j颗GPS/BDS相位观测值，ρAB j表示AB星站间单差几何距离，c表示真空中的光速，δtrAB、δtsAB j表示AB星载接收机的相对钟差和第j颗GPS/BDS的相对钟差，NSD表示单差模糊度，单位为米，εPAB表示单差相位测量噪声。
将单差几何距离ρAB j在概略点ρAB0 j处线性化展开可得到：
Figure BDA0002225353780000173
将式(20)代入式(19)可得到：
Figure BDA0002225353780000174
根据式(21)，单差相位观测值关于相对定轨滤波器状态量的偏导数H可以表示为：
Figure BDA0002225353780000175
其中：UT为J2000惯性系到地固系的转换矩阵，
Figure BDA0002225353780000176
为B星相对于GPS卫星i和BDS卫星j的视线向量，I1×12中第i颗GPS或第j颗BDS对应的通道为1，其余通道为0，
Figure BDA0002225353780000177
表示AB星基线向量的误差。
按式(15)，依次对每颗没有标记粗差的GPS/BDS单差载波相位观测值进行处理，更新相对定轨卡尔曼滤波的状态量及状态协方差阵。
步骤S9-8：AB星轨道预报
根据绝对与相对测量更新后的绝对与相对定轨滤波状态量结果，计算B星星绝对位置速度，分别考虑AB卫星姿态，将AB星天线相位中心的结果转到AB卫星质心，以AB卫星质心结果为轨道积分预报初值，分别对A、B星进行动力学轨道预报，得到轨道插值首尾端点信息。
步骤S10：结束当前历元所有计算任务。
依次不断重复上诉步骤S1～S10，即可通过在轨实时(星上处理)完成编队卫星GPS/BDS相对定轨，得到AB卫星高精度的绝对与相对位置速度等结果。
请参考图2，本发明的一种易于星载在轨实时处理的编队卫星GPS/BDS相对定轨方法，还可以在主从星之间的星间距离较短的情况下，通过传统RTK方式求解双差固定解，然后以双差固定解约束相对定轨滤波器状态量中的待估单差模糊度参数，来提高短基线相对定轨精度，具体的方法是，在前述相对定轨方法的步骤S9-6与步骤S9-7之间，增加如下a～e等步骤：
步骤a，选取双差参考星。
为建立双差观测值，需要在每个导航系统的单差观测值(式(16)中GPS L1/BDS B1单频非差相位观测值计算的单差相位观测值)中，选取一个卫星作为参考星，其它同系统卫星的单差观测值都与其相应的参考星单差观测值求差，形成双差观测数据。参考星的选取应具有以下要求：(1)单差观测值不存在周跳或粗差数据。(2)单差观测值的卫星不是新升起的卫星。(3)上一个历元的参考星如果高度角大于30°，当前历元可以继续被选取为参考星。(4)如果上述3条要求都不能满足，在没有周跳且非新升起的卫星中，选取高度角最大的卫星作为参考星。(5)如果上述条件都不能满足，不选取参考星，不组成双差观测，也不进行载波相位双差相对定位。
步骤b，传递已固定的模糊度参数。
经过步骤S9-4周跳检查，当前历元某卫星单差数据没有发生周跳，将上一历元该卫星已固定的双差固定解模糊度传递到当前历元，即检查GPS和BDS上一个历元与当前历元的双差参考星是否变化，如果有变化，需要对上一历元的固定解模糊度进行组合，才能传递给当前历元。如果没有变化，直接传递到当前历元，并将其不再作为待估双差模糊度参数。最后设置剩下需待估的双差模糊度参数。
步骤c，双差最小二乘相对定位浮点解。
对于主从星间的短基线而言，不考虑电离层延迟和对流层延迟误差的影响，主从星GPS/BDS相对定位的双差伪距与载波相位观测方程可表示为：
Figure BDA0002225353780000191
式(23)中，
Figure BDA0002225353780000192
分别为GPS/BDS卫星P1伪距和L1载波相位的双差观测值，
Figure BDA0002225353780000193
为站星几何距离的双差值，
Figure BDA0002225353780000194
为双差载波相位模糊度，以米为单位，εP、εL分别表示双差伪距和相位测量噪声。
假定主星A的绝对位置以通过绝对定轨滤波得到，主从星的基线向量的初值
Figure BDA0002225353780000195
可通过相对动力学预报得到，可对式(23)进行线性化，可得：
Figure BDA0002225353780000196
式(24)中，A为
Figure BDA0002225353780000197
关于从星B位置向量的偏导数，WP和WL分别是双差伪距和双差相位观测值残差，δx为基线向量的改正量，为待估参数，PP和PL为双差观测值的权矩阵，PP:PL＝1:10000。
如果主从星的基线向量的初值
Figure BDA0002225353780000198
的相对动力学预报精度很高时，可以增加先验基线向量约束：
Figure BDA0002225353780000199
通过最小二乘求解式(25)，可得到双差模糊度浮点解
Figure BDA00022253537800001910
及其协方差矩阵
Figure BDA0002225353780000201
Figure BDA0002225353780000202
表示浮点解对应的基线向量值，
Figure BDA0002225353780000203
表示由相对动力学预报的先验基线向量
Figure BDA0002225353780000204
的精度，δx为基线向量的改正量(相对于基线向量的初值
Figure BDA0002225353780000205
)。
步骤d，双差模糊度固定解及确认。
以步骤c中的双差模糊度浮点解
Figure BDA0002225353780000206
及其协方差矩阵
Figure BDA0002225353780000207
为输入，首先采用最小二乘去相关算法LAMBDA进行搜索得到双差模糊度的整数解，再根据式(26)计算Ratio、RMS、δdR等值。
Figure BDA0002225353780000208
式(26)中，RS,RB分别为LAMBDA方法输出的次优整周模糊度和最优整周模糊度对应的残差量，V为由最优整周模糊度计算的观测值残差，P为观测值权阵，n为观测值总数，dRekf,dRrtk分别为由先验相对动力学轨道预报的基线长和由最优整周模糊度候选组合计算的基线长，δdR为两者基线长度差值。
记
Figure BDA0002225353780000209
为条件a，其中RatioLim、RMSLim、δdRLim为相应的限差门限值，可通过上注的方式修改，同时记满足条件a的累计次数为M，确认最优整周模糊度的候选组合是否为正确固定的模糊度规则为：(1)在上一个历元模糊度没有固定情况下，当M≥2时，将最优整周模糊度的候选组合设置都设置为固定状态，否则认为还没有固定。(2)当上一个历元模糊度已固定时，需要满足条件a的检验，将步骤b中待估的模糊度固定状态设置都设定为固定。(3)如果RMS没有通过条件a中的限差检验或者双差卫星数少于等于3颗，将已固定的模糊度状态都设置为非固定，即下个滤波周期全部重新搜索模糊度。(4)其它情况，如Ratio或δdR未通过条件a中的相应限差检验，保持现有的模糊度固定状态情况不变。
步骤e，双差模糊度固定解约束。
具体方法为：根据双差和单差观测方程之间的关系，双差相位模糊度观测方程可表示为：
Figure BDA0002225353780000211
根据式(27)，双差相位模糊度L关于相对定轨滤波器状态量的偏导数H可以表示为：
Figure BDA0002225353780000212
其中：I1×12中第i颗GPS/BDS对应的通道值为1，第j颗为GPS/BDS双差参考星对应的通道值为-1，其余通道为0，
Figure BDA0002225353780000213
表示第i颗GPS/BDS卫星对应的双差模糊度，
Figure BDA0002225353780000214
表示第i颗GPS/BDS卫星对应的单差模糊度，
Figure BDA0002225353780000215
表示第j颗GPS/BDS双差参考星对应的单差模糊度。
最后按式(15)，将单差模糊度参数约束为双差固定解，即将经过步骤d检验后的双差相位模糊度固定解NDD作为观测量，依次对每颗GPS/BDS进行处理，更新相对定轨卡尔曼滤波的状态量及状态协方差阵。
如图3所示，为本发明实施例的GRACE编队卫星相对定轨结果图。是将GRACE-A/B卫星在2005年DOY343-345共三天的在轨实测数据作为输入数据流，在地面PC机上仿真模拟在轨实时处理模式，经过本发明的S1-S10步骤进行连续定轨解算，输出相对定轨结果。其中，中长基线(大于15Km)采用双频无电离层组合浮点解进行相对定轨，短基线(小于15Km)采用单频双差固定解进行相对定轨(即还经过了本发明实施例中的步骤a～e进行数据处理)，最后以JPL、ESA等机构公布的后处理精密轨道(精度约为3cm)为参考，对两者进行求差处理，统计分析实时相对定轨结果的精度状况。经实测统计，整个过程相对定轨精度：R方向为2～3cm，T方向为2～3cm，N方向为2cm，3d约为3～5cm，在基线长度BL<15km区间内，R：9mm；T：12mm；N：13mm；3d：21mm，采用本发明的方法，短基线实时相对定轨精度可优于1cm。
如图4所示，为使本发明中的轨道预报配合轨道内插的运行模式更易于理解，进一步的通过时序图对本发明实施例进行详细说明。
在本发明实施例中，在步骤S6和S9中完成AB星轨道预报，生成轨道插值首尾端点信息，然后配合步骤S7，通过轨道内插，即可输出频率可调的编队卫星高精度绝对与相对轨道结果。同时上述模式下，为保证运算资源的合理分配，将步骤S6和S9中的多个子步骤的计算量，通过分时多步计算来完成，以降低步骤S6和S9单历元的运算耗时，使算法更易于在星载嵌入式处理器上实时处理。
具体为，作为一种具体的实施例，一般地，可根据各子步骤的实际计算量，将步骤S6计算任务分为3步完成，将步骤S9(含步骤a～e)计算任务分为24步完成，步骤S6和步骤S9(含步骤a～e)的分步计算规划如下表1和表2所示：
表1步骤S6分步计算规划
步数 计算内容
1 S6-1～S6-4
2 S6-5中A轨道预报
3 S6-5中B轨道预报、S6-6
表2步骤S9(含步骤a～e)分步计算规划
步数 计算内容 步数 计算内容 步数 计算内容
1 S9-1 9 S9-3中部分计算 17 步骤e中部分计算
2 S9-2中部分计算 10 S9-4、S9-5中部分 18 步骤e中部分计算
3 S9-2中部分计算 11 S9-5中部分计算 19 S9-7中部分计算
4 S9-2中部分计算 12 S9-5中部分计算 20 S9-7中部分计算
5 S9-2中部分计算 13 S9-5中部分计算 21 S9-7中部分计算
6 S9-3中部分计算 14 S9-6步骤a～d 22 S9-7中部分计算
7 S9-3中部分计算 15 步骤e中部分计算 23 S9-8中A轨道预报
8 S9-3中部分计算 16 步骤e中部分计算 24 S9-8中B轨道预报
如图4和表3所示，设定轨道滤波启动时刻为T0，经过T0～T0+3时段完成表1中分步规划的3步(步骤S6计算任务)，从而计算出了(T0,T0+h)轨道插值首尾端点，记为轨道预报周期1，那么在时间区间(T0+4,T0+step1+step2)内的绝对和相对定轨结果，可以根据轨道预报周期1的结果通过步骤S7轨道内插获得，在T0+step1时刻开始执行步骤S6定轨滤波解算，经过T0+step1～T0+step1+step2时段完成表2中分步规划的24步(步骤S9计算任务)，从而计算出了(T0+step1,T0+step1+h)轨道插值首尾端点，记为轨道预报周期2，那么在时间区间(T0+step1+step2+1,T0+2*step1+step2)内的绝对和相对定轨结果，可以根据轨道预报周期2的结果通过步骤7轨道内插获得，依次类推。
在本发明实施例中，T0表示轨道滤波启动时刻，step1表示定轨测量更新间隔，step2表示定轨滤波解算分步数，h表示轨道预报步长，一般地取，step1＝30，step2＝24，h＝55。
表3轨道插值输出所用到的插值端点信息
Figure BDA0002225353780000231
如图5所示，为本发明另一实施例提供的1s计算10步定轨滤波解算过程的计算时序图。
在本发明实施例中，GPS/BDS接收机中的星载处理器除了有在轨实时相对定轨数据处理任务外，还包括通过定时中断的方式(如1秒10次中断)，完成对GPS/BDS导航信号进行捕获、跟踪、导航电文解译及GPS/BDS原始观测量生成等任务。假设该任务运算耗时在10ms以内，则可以充分利用星载处理器剩下的90ms空闲时间完成相对定轨的计算任务，即可将步骤S6和步骤S9(含步骤a～e)的各步，以1s计算10步的方式进行数据处理。那么上述规划为24步的步骤S9(含步骤a～e)实际只需2.4s即可全部计算完成，从而可缩小步骤S6-5和步骤S9-8中轨道预报步长h，进而可以提高最终的绝对与相对定轨精度。
进一步的，本发明的另一实施例，是在上述1s计算10步定轨滤波解算过程的基础上，还可以将与导航卫星数有关的步骤(如步骤S9-3、步骤e、步骤S9-7等)再更多的分步，以进一步减小单历元的耗时，从而降低对星载处理器性能要求。同时，还可以在步骤S6-5和S9-8完成AB星轨道预报之后，通过增加一步或多步的轨道预报(预报步长不超过60s)计算任务，以增大轨道插值尾端点的时间跨度，再经过步骤S7轨道内插，可预报距离当前历元几分钟之后的未来时刻的高精度的绝对与相对轨道结果。
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (6)

1.一种易于星载在轨实时处理的编队卫星相对定轨方法，其特征在于，包括：
(1)在当前历元时刻的星载GPS/BDS实时相对定轨开始后，判断是否需要进行初始化，若需要进行初始化，则进行初始化操作，若不需要进行初始化，则进入步骤(2)；
(2)获取定轨实时数据，根据所述定轨实时数据得到双频无电离层伪距和相位LC组合观测值、电离层残差GF组合观测值及MW组合观测值；
(3)根据所述电离层残差GF组合观测值及所述MW组合观测值对非差相位数据进行实时周跳探测，并标记存在周跳的非差相位观测值；
(4)判断定轨滤波器启动标记是否为已启动，若所述启动标记为已启动，则进入步骤(5)，若所述启动标记为未启动，则需完成定轨滤波的启动，并将所述启动标记设置为已启动，并结束当前历元任务；
其中，步骤(4)中的定轨滤波启动包括如下子步骤：
(4.1)使用编队AB卫星星载GPS/BDS双频伪距及多普勒观测数据分别进行标准单点定位与单点测速，分别得到AB卫星天线相位中心在地心地固系下的位置和速度、接收机钟差和钟漂参数及定位测速中伪距多普勒观测值残差的RMS值；
(4.2)组成AB星站间伪距单差观测值，进行伪距单差相对定位，得到AB卫星天线相位中心之间基线向量在地心地固系下的相对位置和相对钟差；
(4.3)若不满足定轨滤波启动条件，则直接结束当前历元任务，若满足启动条件，则执行步骤(4.4)；
(4.4)将AB卫星天线相位中心在地心地固系下的位置及AB卫星天线相位中心之间基线向量在地心地固系下的相对位置，从地心地固系转换到地心惯性系，再初始化绝对与相对定轨滤波器的滤波状态量及状态误差协方差阵；
(4.5)考虑AB卫星姿态，分别将AB星在地心惯性系下的天线相位中心转到AB卫星质心，以AB卫星质心为轨道积分预报初值，使用AB星动力学模型参数，分别对A星和B星进行动力学轨道预报，得到轨道插值首尾端点信息，并将定轨滤波器启动标记设置为已启动，结束当前历元任务；
(5)由AB星轨道预报得到的轨道插值首尾端点信息，使用Hermit插值分别得到AB星在当前历元整秒时刻的单星绝对位置速度结果及相对位置速度结果；
(6)由当前观测历元时间与插值首端点信息时间的间隔判断是否达到定轨滤波时刻，若已达到定轨滤波时刻，则进入步骤(7)，若未达到定轨滤波时刻，则结束当前历元任务；
(7)基于不存在周跳的非差/单差相位观测值、所述AB星在当前历元整秒时刻的单星绝对位置速度结果及相对位置速度结果进行定轨滤波解算；
(8)重复步骤(1)～步骤(7)，通过在轨实时完成编队卫星GPS/BDS相对定轨。
2.根据权利要求1所述的方法，其特征在于，步骤(7)包括：
(7.1)使用星载GPS/BDS双频伪距和多普勒观测数据对编队AB卫星分别进行标准单点定位，得到AB卫星天线相位中心在地心地固系下的位置和AB星载接收机钟差参数；
(7.2)根据AB卫星天线相位中心在地心地固系下的位置和AB星载接收机钟差参数，完成A星绝对定轨时间更新及A星绝对定轨测量更新，依次对每颗没有标记粗差的GPS/BDS非差载波相位观测值进行处理，更新卡尔曼滤波的状态及协方差参数；
(7.3)生成AB星站间单差数据，并根据AB星站间单差数据得到单差数据的电离层残差GF组合观测值和MW组合观测值，根据单差数据的电离层残差GF组合观测值和MW组合观测值完成对单差数据的周跳探测，标记存在周跳的相位单差观测值，其中，存在周跳的相位单差观测值不参与后续的相对定轨测量更新；
(7.4)完成AB星相对定轨时间更新，采用基于相对动力学先验轨道来探测单差相位观测值的验前粗差，标记出存在验前粗差的单差数据，其中，存在验前粗差的单差数据不参与后续的相对定轨测量更新；
(7.5)完成AB星相对定轨测量更新，依次对每颗没有标记粗差的GPS/BDS单差载波相位观测值进行处理，更新卡尔曼滤波的状态及协方差参数；
(7.6)考虑AB卫星姿态，分别将AB星测量更新的结果转到AB卫星质心，以AB卫星质心结果为轨道积分预报初值，分别对A星和B星进行动力学轨道预报，得到轨道插值首尾端点信息。
3.根据权利要求2所述的方法，其特征在于，在步骤(7.2)中，所述完成A星绝对定轨时间更新包括：
考虑A星卫星姿态，先将上个历元A星测量更新的结果转到A卫星质心，再使用A星动力学模型参数，经动力学轨道积分得到的A星卫星状态及状态转移矩阵，然后再考虑A星卫星姿态，将A星质心的结果转到A星天线相位中心，最后将绝对定轨滤波的状态和状态误差协方差矩阵更新到当前测量更新时刻。
4.根据权利要求2或3所述的方法，其特征在于，在步骤(7.4)中，所述完成AB星相对定轨时间更新包括：
考虑AB星卫星姿态，利用上一历元A星绝对定轨测量更后的结果和相对定轨测量更新结果，得到B星绝对位置速度结果，将B星天线相位中心转到B卫星的质心处，再使用B星动力学模型参数，经动力学轨道积分得到的B星卫星状态及状态转移矩阵，再将B星质心的结果转到B星天线相位中心，然后减去A星时间更新后的结果，得到AB星相对轨道状态及相对状态转移矩阵，最后将相对定轨滤波的状态和状态误差协方差矩阵更新到当前测量更新时刻。
5.根据权利要求4所述的方法，其特征在于，在步骤(7.2)和步骤(7.4)中的动力学轨道积分过程，均使用经典4阶龙格库塔法对卫星运动方程和状态转移矩阵微分方程进行数值积分。
6.根据权利要求2所述的方法，其特征在于，在步骤(4.5)和步骤(7.6)中的AB星动力学轨道积分预报中使用的重力场模型阶数均为30阶，在步骤(7.2)中的A星绝对定轨时间更新和步骤(7.4)中的AB星相对定轨时间更新中使用的重力场模型阶数均为70阶。
CN201910949615.1A 2019-10-08 2019-10-08 易于星载在轨实时处理的编队卫星相对定轨方法 Active CN110764127B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910949615.1A CN110764127B (zh) 2019-10-08 2019-10-08 易于星载在轨实时处理的编队卫星相对定轨方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910949615.1A CN110764127B (zh) 2019-10-08 2019-10-08 易于星载在轨实时处理的编队卫星相对定轨方法

Publications (2)

Publication Number Publication Date
CN110764127A CN110764127A (zh) 2020-02-07
CN110764127B true CN110764127B (zh) 2021-07-06

Family

ID=69331240

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910949615.1A Active CN110764127B (zh) 2019-10-08 2019-10-08 易于星载在轨实时处理的编队卫星相对定轨方法

Country Status (1)

Country Link
CN (1) CN110764127B (zh)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111487660B (zh) * 2020-04-24 2022-07-26 北京航空航天大学 一种高精度实时微纳卫星集群导航方法
CN111953401B (zh) * 2020-07-28 2022-06-07 中国西安卫星测控中心 一种微小卫星自主请求式轨道服务系统
CN112649828B (zh) * 2020-11-30 2024-03-01 中国科学院国家天文台 倾斜高圆轨道通信卫星的定轨方法、系统及设备
CN114779285A (zh) * 2022-04-18 2022-07-22 浙江大学 基于微小型低功耗星载双模四频gnss接收机的精密定轨方法
CN114740541B (zh) * 2022-06-09 2022-09-13 武汉大学 基于主从星测速模式的小行星重力场反演方法及系统
CN115598673B (zh) * 2022-09-29 2023-10-24 同济大学 Igs gnss卫星钟差和轨道单天相邻产品边界处偏差计算方法

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2004210032A (ja) * 2002-12-27 2004-07-29 Mitsubishi Electric Corp 編隊飛行衛星
CN102230969A (zh) * 2011-03-22 2011-11-02 航天恒星科技有限公司 一种卫星星座星间链路的长时间自主维持方法
CN103675861A (zh) * 2013-11-18 2014-03-26 航天恒星科技有限公司 一种基于星载gnss多天线的卫星自主定轨方法
CN104142686A (zh) * 2014-07-16 2014-11-12 北京控制工程研究所 一种卫星自主编队飞行控制方法
CN105372691A (zh) * 2015-08-18 2016-03-02 中国人民解放军国防科学技术大学 一种模糊度固定的长基线卫星编队gnss相对定位方法
CN105468882A (zh) * 2014-07-28 2016-04-06 航天恒星科技有限公司 卫星自主定轨方法及系统
CN105652308A (zh) * 2014-11-27 2016-06-08 航天恒星科技有限公司 飞行器相对测量方法及系统

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6859690B2 (en) * 2002-03-13 2005-02-22 The Johns Hopkins University Method for using GPS and crosslink signals to correct ionospheric errors in space navigation solutions

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2004210032A (ja) * 2002-12-27 2004-07-29 Mitsubishi Electric Corp 編隊飛行衛星
CN102230969A (zh) * 2011-03-22 2011-11-02 航天恒星科技有限公司 一种卫星星座星间链路的长时间自主维持方法
CN103675861A (zh) * 2013-11-18 2014-03-26 航天恒星科技有限公司 一种基于星载gnss多天线的卫星自主定轨方法
CN104142686A (zh) * 2014-07-16 2014-11-12 北京控制工程研究所 一种卫星自主编队飞行控制方法
CN105468882A (zh) * 2014-07-28 2016-04-06 航天恒星科技有限公司 卫星自主定轨方法及系统
CN105652308A (zh) * 2014-11-27 2016-06-08 航天恒星科技有限公司 飞行器相对测量方法及系统
CN105372691A (zh) * 2015-08-18 2016-03-02 中国人民解放军国防科学技术大学 一种模糊度固定的长基线卫星编队gnss相对定位方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
《GPS 单差观测量的运动学相对定轨研究》;刘荣芳，等;《载人航天》;20150731;第21卷(第4期);正文第2节及附图1 *
《星载GPS自主定轨理论及其软件实现》;王甫红;《中国博士学位论文全文数据库 基础科学辑》;20081115;正文第2.5.3-2.5.4，4.2.3节 *

Also Published As

Publication number Publication date
CN110764127A (zh) 2020-02-07

Similar Documents

Publication Publication Date Title
CN110764127B (zh) 易于星载在轨实时处理的编队卫星相对定轨方法
CN108120994B (zh) 一种基于星载gnss的geo卫星实时定轨方法
Loyer et al. Zero-difference GPS ambiguity resolution at CNES–CLS IGS Analysis Center
CN107728180B (zh) 一种基于多维粒子滤波偏差估计的gnss精密定位方法
CN104714244A (zh) 一种基于抗差自适应Kalman滤波的多系统动态PPP解算方法
CN108196284B (zh) 一种进行星间单差模糊度固定的gnss网数据处理方法
CN111239787A (zh) 一种集群自主协同中的gnss动态卡尔曼滤波方法
CN107966722B (zh) 一种gnss钟差解算方法
Hwang et al. GPS‐Based Orbit Determination for KOMPSAT‐5 Satellite
CN113624243B (zh) 一种近地轨道卫星的星上实时轨道预报方法
CN111487660B (zh) 一种高精度实时微纳卫星集群导航方法
CN107561565A (zh) 低轨航天器星载gnss差分相对导航换星的处理方法
CN110673175A (zh) 一种基于gnss广播星历的高轨卫星高精度自主定轨方法
WO2023167899A1 (en) System and method for fusing sensor and satellite measurements for positioning determination
Sun et al. Precise real-time navigation of LEO satellites using a single-frequency GPS receiver and ultra-rapid ephemerides
CN108205151B (zh) 一种低成本gps单天线姿态测量方法
Wu et al. Real-time sub-cm differential orbit determination of two Low-Earth Orbiters with GPS bias fixing
CN116594046B (zh) 基于低轨卫星信号多普勒误差补偿的运动目标定位方法
CN110988932B (zh) 一种提高星载gps接收机实时钟差解算精度的方法
CN114063122B (zh) 电推进转移轨道航天器星载gnss在轨实时定轨方法
CN112731504B (zh) 对月球探测器自主定轨的方法及装置
CN113671551B (zh) Rtk定位解算方法
Lee Introduction to navigation systems
Jwo Complementary Kalman filter as a baseline vector estimator for GPS-based attitude determination
CN115859560B (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