阅读说明：本技术 易于星载在轨实时处理的编队卫星相对定轨方法 (Relative orbit determination method for formation satellite easy for satellite-borne on-orbit real-time processing ) 是由 张万威 王甫红 龚学文 郭磊 于 2019-10-08 设计创作，主要内容包括：本发明公开了一种易于星载在轨实时处理的编队卫星相对定轨方法,属于卫星编队相对导航领域,该方法使用GPS/BDS广播星历和简化的相对动力学模型,对编队卫星GPS/BDS双频观测数据进行在轨实时处理,考虑星载处理器计算资源有限,设计了动力学轨道预报配合轨道内插的运行模式,并将核心复杂数据处理过程通过分时多步实现来降低单历元的计算耗时,可输出频率可调的编队卫星高精度绝对与相对轨道结果。本发明不仅适用于较大星间距离范围的星间编队,在GPS/BDS信号中断情况下,仍可连续输出满足一定精度要求较为平滑的导航结果,而且具有定轨精度高、稳定性好及对星载处理器性能要求较低,易于工程化实现等特点。(The invention discloses a relative orbit determination method for formation satellites easy to carry out satellite-borne on-orbit real-time processing, which belongs to the field of satellite formation relative navigation. The invention is not only suitable for the inter-satellite formation with a larger inter-satellite distance range, but also can continuously output a navigation result which meets a certain precision requirement and is smoother under the condition of GPS/BDS signal interruption, and has the characteristics of high orbit determination precision, good stability, lower requirement on the performance of a satellite-borne processor, easy engineering realization and the like.)

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

技术领域

本发明属于卫星编队相对导航技术领域，更具体地，涉及一种易于星载在轨实时处理的编队卫星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星星载接收机钟差参数(δt_A,δt_B)，即步骤(7.2)中绝对定轨时间更新间隔或

步骤(7.4)中相对定轨时间更新间隔

或

(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组合观测值可表示为：

电离层残差GF组合观测值可表示为：

L_gf＝L₁-L₂ (2)

MW组合观测值可表示为：

其中，f₁、f₂为GPS/BDS双频观测数据对应的信号频点，L_C和P_C分别表示GPS/BDS双频无电离层相位和伪距的组合观测值，L₁和L₂分别表示GPS L1/BDS B1和GPS L2/BDS B2(或B3)频点的相位观测值，P₁和P₂分别表示GPS L1/BDS B1和GPS L2/BDS B2(或B3)频点的伪距观测值，L_gf表示电离层残差GF组合观测值，L_mw表示MW组合观测值。

步骤S4：星载GPS/BDS双频非差数据周跳探测。

根据GF和MW组合观测值对非差相位数据进行实时周跳探测，标记存在周跳的非差相位观测值，不参与后续的绝对定轨测量更新。

具体方法是：首先根据式(4)计算历元间GF和MW组合观测值变化量dL_gf和dL_mw，当同时满足dL_gf＜GFThreshold和dL_mw＜MWThreshold时，将该非差相位数据标记为1，即认为不存在周跳，否则标记为-2，即认为存在周跳，后续不参与绝对定轨滤波解算。GFThreshold和MWThreshold为周跳探测的具体阈值门限，可通过上注方式修改。

其中，L_gf(t)、L_gf(t-1)分别表示当前/上一时刻的电离层残差GF组合观测值，L_mw(t)、L_mw(t-1)分别表示当前/上一时刻的MW组合观测值。

步骤S5：判断定轨滤波器是否已经启动。

即判断定轨滤波器启动标记是否为已启动，如果该标记为已启动，则直接进入步骤S7，否则进入步骤S6。

步骤S6：执行完成定轨滤波的启动；

步骤S6-1：AB星单点定位测速

使用编队AB卫星星载GPS/BDS双频伪距、多普勒观测数据分别进行标准单点定位(SPP)与单点测速(SPV)，得到AB卫星天线相位中心在地心地固系(WGS84)下的位置速度和接收机钟差

和

钟漂

参数和定位测速中伪距多普勒残差的RMS值

其中

和分别为AB卫星星载GPS/BDS接收机关于GPS和BDS导航系统的钟差。

步骤S6-2：AB星站间伪距单差相对定位

组成AB星站间伪距单差观测值，进行伪距单差相对定位，得到AB卫星天线相位中心之间基线向量在地心地固系(WGS84)下的相对位置

和相对钟差

其中，

基线向量的相对速度

步骤S6-3：判断是否满足定轨滤波启动条件

如果不满足定轨滤波启动条件，则直接结束当前历元所有计算任务(步骤S10)，如果满足启动条件，则进行步骤S6-4。启动条件为同时满足

其中σ_p和σ_v是参考阈值门限，可通过上注方式修改。

步骤S6-4：绝对与相对定轨滤波器初始化

将S6-1和S6-2计算的结果

和

从地心地固系(WGS84)转换到地心惯性系(J2000)，得到在地心惯性系下结果

和

再初始化绝对与相对定轨滤波器的滤波状态量(X_A,X_AB)及状态误差协方差阵(P_A,P_AB)等，其中：

滤波状态量X_A,X_AB初始值设置如下：

C_d、dC_d和C_r、dC_r分别为待估的绝对与相对大气阻力系数和太阳光压系数，W和dW分别为RTN三个方向待估的动力学模型补偿加速度参数，N^G和N^C分别为GPS和BDS对应12通道的非差相位模糊度参数(单位为m)，P和L分别为由式(1)计算的双频无电离层伪距和相位LC组合观测值。

和

分别为GPS和BDS对应12通道的单差相位模糊度参数(单位为m)，P_SD和L_SD分别为AB星由对应的共视卫星P和L计算的站间单差伪距和相位观测值，上标G和C分别表示对应GPS和BDS系统的各变量。

状态误差协方差阵P_A,P_AB初始值可设置为：

P_A＝P_AB＝diag(1E4,1E2,1E4,1E4,10,10,0.01,SigN2,SigN2)

其中，SigN2为相位模糊度的初始方差，可通过上注方式修改。

步骤S6-5：AB星轨道预报

考虑AB卫星姿态，分别将AB星初始轨道状态

和

由天线相位中心转到AB卫星质心

和

其中

然后以AB卫星质心为轨道积分预报初值，使用AB星动力学模型参数，分别对A、B星进行动力学轨道预报，得到AB星轨道插值所需的首尾端点信息

和

其中，t,r,v,a分别为时间、位置、速度、加速度。

所述的动力学轨道预报过程，是使用经典4阶龙格－库塔法(RK4)，对考虑了地球引力、日月引力、固体潮、大气阻力及太阳光压等摄动力影响的卫星运动方程进行数值积分，该过程积分步长为h，使用的地球重力场模型阶数为30阶。

步骤S6-6：置标记为滤波器已启动

将定轨滤波器启动标记设置为已启动，结束当前历元所有计算任务(第10步)。

步骤S7：轨道内插输出。

即由定轨滤波启动步骤S6或定轨滤波解算步骤S9中AB星轨道预报得到的轨道插值首尾端点信息，使用5阶Hermit插值分别计算AB星在当前历元整秒时刻的单星绝对位置速度结果，再计算相对位置速度结果；

具体方法为：对于任一时刻t∈(t₀,t₁)的卫星状态(r(t),v(t))，可用一个五阶的Hermite多项式近似表示如下：

所述的式(7)中的位置、速度系数项计算方法如下所示：

由步骤S6-5或步骤S9-8计算得到的AB卫星轨道插值首尾端点信息结果和

经式(7)计算AB星在当前历元整秒时刻t的单星绝对位置速度结果(r_A(t),v_A(t))和(r_B(t),v_B(t))，再计算AB星相对位置速度结果(r_AB(t),v_AB(t))，其中r_AB(t)＝r_B(t)-r_A(t),v_AB(t)＝v_B(t)-v_A(t)。

步骤S8：判断是否达到定轨滤波时刻。

即计算当前观测历元时间与插值前首端点信息时间的间隔dt，如果dt取整后大于等于卡尔曼滤波时间间隔step，则已达到定轨滤波时刻，需要进入步骤S9，否则进入步骤S10。

步骤S9：执行完成定轨滤波解算任务。

步骤S9-1：AB星单点定位

使用星载GPS/BDS双频伪距、多普勒观测数据对编队AB卫星分别进行标准单点定位(SPP)，得到AB卫星天线相位中心在地心地固系(WGS84)下的位置和AB星载接收机钟差参数和

步骤S9-2：A星绝对定轨时间更新

设绝对定轨卡尔曼滤波的状态方程为：

式(8)中

为绝对定轨卡尔曼滤波状态量，详见步骤S6-4，

为状态转移矩阵，W_k为系统噪声矩阵。

Φ_rr,Φ_rv,Φ_vr,Φ_vv分别为位置和速度分量的状态转移矩阵，

表示位置和速度分别关于大气阻力和太阳光压的状态转移矩阵，Φ_rw,Φ_vw分别位置和速度关于RTN方向的补偿加速度的状态转移矩阵，Φ_ww为补偿加速度关于其自身的状态转移矩阵，为对应GPS/BDS通道非差模糊度关于其自身的状态转移矩阵，下表k表示当前时刻，k-1表示上一时刻。

绝对定轨卡尔曼滤波的观测方程为：

其中，

为观测矢量，

为观测矩阵，V_k为观测噪声矢量，且系统噪声W_k和观测噪声V_k，为零均值白噪声序列，即对所有的k,j，有

式中，Q_k和R_k分别为系统动态噪声协方差阵和观测噪声协方差阵，δ_kj为Dirac-δ函数。

A星绝对定轨时间更新具体过程是，考虑A星卫星姿态，先将上个历元A星测量更新的结果转到A卫星质心

以

为初始值，使用A星动力学模型参数，经动力学轨道积分得到的滤波测量更新时刻的A星卫星状态

及状态转移矩阵

然后再考虑A星卫星姿态，将积分结果A星质心的结果

转到A星天线相位中心

最后使用式(8)和式(12)将绝对定轨滤波的状态量和状态误差协方差矩阵更新到当前测量更新时刻；

其中，

表示当前时刻的时间更新后的状态误差协方差矩阵，

表示上一时刻的时间更新后的状态误差协方差矩阵。

所述的动力学轨道积分过程，是使用经典4阶龙格－库塔法(RK4)，同时对卫星运动方程和状态转移矩阵微分方程进行数值积分，该过程实际积分步长为

或

使用的地球重力场模型阶数为70阶。

步骤S9-3：A星绝对定轨测量更新

双频无电离层LC组合的非差相位观测值可以表示为：

其中，ρ^j表示GPS/BDS发射与接收天线相位中心之间的几何距离，c表示真空中的光速，δt_r、

分别表示接收机和第j颗导航卫星的钟差，N表示不具有整数特性的非差模糊度，单位米，ε_P表示相位测量噪声。

式(13)中非差相位观测值关于绝对定轨滤波器状态量的偏导数H可以表示为：

式(14)中：U^T为J2000惯性系到地固系的转换矩阵，

为A星相对于GPS卫星i和BDS卫星j的视线向量，I_1×12中第i颗GPS或第j颗BDS对应的通道为1，其余通道为0，下标G和C分别表示GPS和BDS系统。

按式(15)，依次对每颗没有标记粗差的GPS/BDS非差载波相位观测值进行处理，更新绝对定轨卡尔曼滤波的状态量及状态协方差阵；

其中，

表示当前时刻卡尔曼增益矩阵，

表示式(14)中的第i颗对应的观测矩阵，

表示式(11)中的观测噪声协方差阵，表示测量更新后的滤波状态量，

表示时间更新后的滤波状态量，

表示式(10)中的观测值，

表示当前时刻测量更新后的状态误差协方差矩阵。

步骤S9-4：生成AB星站间单差观测数据与观测值域周跳探测

首先按照式(16)生成AB星站间单差伪距和相位数据，计算单差数据的电离层残差GF组合和MW组合观测值。

式(16)中，

分别为AB星GPS/BDS伪距和相位无电离层LC组合观测值，P₁ ^B,P₁ ^A,

分别为AB星GPS/BDS导航系统对应的双频非差伪距和相位观测值，f₁和f₂为GPS/BDS双频信号频点，

和分别表示AB星GPS/BDS双频伪距和相位无电离层组合计算的单差相位观测值。

根据单差数据的GF和MW组合观测值对单差数据进行实时周跳探测，标记存在周跳的观测值。即根据(17)式计算历元间单差数据的GF和MW组合观测值变化量

和

当同时满足

和

时，将该单差数据标记为1，即认为不存在周跳，否则标记为-2，即认为存在周跳，后续不参与相对定轨滤波解算。

其中，分别表示当前/上一时刻的单差数据的GF组合观测值

分别表示当前/上一时刻的单差数据的MW组合观测值。

步骤S9-5：AB星相对定轨时间更新

设相对定轨卡尔曼滤波的状态方程和观测方程为：

式(18)中，

为相对定轨卡尔曼滤波状态量，详见步骤S6-4，此外可推导得出

表示GPS/BDS单差观测值，

表示观测矩阵。

因此，相对定轨时间更新可参考绝对定轨时间更新过程，具体为：考虑AB星卫星姿态，利用上一历元A星绝对定轨测量更新的结果

相对定轨测量更新结果

得到B星绝对位置速度结果将该结果从天线相位中心转到B卫星的质心处

再使用B星动力学模型参数，经动力学轨道积分得到滤波更新时刻的B星卫星状态及状态转移矩阵

再将B星质心的结果转到B星天线相位中心

然后减去A星时间更新后的结果得到AB星相对轨道状态

最后式(18)和式(12)将相对定轨滤波的状态量和状态误差协方差矩阵更新到当前测量更新时刻，相对定轨时间更新过程实际积分步长为

或

使用的地球重力场模型阶数为70阶。

步骤S9-6：基于动力学先验轨道的单差验前粗差探测；

具体过程为：以时间更新的滤波状态量结果为先验轨道，首先计算单差相位数据的残差向量Val，使用Grubbs准则计算Val中元素的均值Mean和标准差Std，对于每颗卫星，当满足|Val-Mean|>3*Std时，标记该颗卫星为验前粗差，并不参与后续的相对定轨滤波解算。

步骤S9-7：相对定轨测量更新；

单差相位观测值可表示为：

L_AB ^j＝L_B ^j-L_A ^j＝ρ_AB ^j+c(δt_rAB-δt_sAB ^j)+N_SD+ε_PAB (19)

式(19)中的单差相位观测值，对于中长基线，使用的是式(16)中AB星GPS/BDS双频相位无电离层组合计算的单差相位观测值

对于短基线，使用的是式(16)中AB星GPSL1/BDS B1单频非差相位观测值计算的单差相位观测值

L_AB ^j表示AB星GPS/BDS单差相位观测值，L_A ^j、L_B ^jAB星第j颗GPS/BDS相位观测值，ρ_AB ^j表示AB星站间单差几何距离，c表示真空中的光速，δt_rAB、δt_sAB ^j表示AB星载接收机的相对钟差和第j颗GPS/BDS的相对钟差，N_SD表示单差模糊度，单位为米，ε_PAB表示单差相位测量噪声。

将单差几何距离ρ_AB ^j在概略点ρ_AB0 ^j处线性化展开可得到：

将式(20)代入式(19)可得到：

根据式(21)，单差相位观测值关于相对定轨滤波器状态量的偏导数H可以表示为：

其中：U^T为J2000惯性系到地固系的转换矩阵，

为B星相对于GPS卫星i和BDS卫星j的视线向量，I_1×12中第i颗GPS或第j颗BDS对应的通道为1，其余通道为0，

表示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相对定位的双差伪距与载波相位观测方程可表示为：

式(23)中，

分别为GPS/BDS卫星P1伪距和L1载波相位的双差观测值，

为站星几何距离的双差值，

为双差载波相位模糊度，以米为单位，ε_P、ε_L分别表示双差伪距和相位测量噪声。

假定主星A的绝对位置以通过绝对定轨滤波得到，主从星的基线向量的初值

可通过相对动力学预报得到，可对式(23)进行线性化，可得：

式(24)中，A为

关于从星B位置向量的偏导数，W_P和W_L分别是双差伪距和双差相位观测值残差，δx为基线向量的改正量，为待估参数，P_P和P_L为双差观测值的权矩阵，P_P:P_L＝1:10000。

如果主从星的基线向量的初值

的相对动力学预报精度很高时，可以增加先验基线向量约束：

通过最小二乘求解式(25)，可得到双差模糊度浮点解

及其协方差矩阵

表示浮点解对应的基线向量值，

表示由相对动力学预报的先验基线向量的精度，δx为基线向量的改正量(相对于基线向量的初值

)。

步骤d，双差模糊度固定解及确认。

以步骤c中的双差模糊度浮点解

及其协方差矩阵

为输入，首先采用最小二乘去相关算法LAMBDA进行搜索得到双差模糊度的整数解，再根据式(26)计算Ratio、RMS、δdR等值。

式(26)中，R_S,R_B分别为LAMBDA方法输出的次优整周模糊度和最优整周模糊度对应的残差量，V为由最优整周模糊度计算的观测值残差，P为观测值权阵，n为观测值总数，dR_ekf,dR_rtk分别为由先验相对动力学轨道预报的基线长和由最优整周模糊度候选组合计算的基线长，δdR为两者基线长度差值。

记

为条件a，其中RatioLim、RMSLim、δdRLim为相应的限差门限值，可通过上注的方式修改，同时记满足条件a的累计次数为M，确认最优整周模糊度的候选组合是否为正确固定的模糊度规则为：(1)在上一个历元模糊度没有固定情况下，当M≥2时，将最优整周模糊度的候选组合设置都设置为固定状态，否则认为还没有固定。(2)当上一个历元模糊度已固定时，需要满足条件a的检验，将步骤b中待估的模糊度固定状态设置都设定为固定。(3)如果RMS没有通过条件a中的限差检验或者双差卫星数少于等于3颗，将已固定的模糊度状态都设置为非固定，即下个滤波周期全部重新搜索模糊度。(4)其它情况，如Ratio或δdR未通过条件a中的相应限差检验，保持现有的模糊度固定状态情况不变。

步骤e，双差模糊度固定解约束。

具体方法为：根据双差和单差观测方程之间的关系，双差相位模糊度观测方程可表示为：

根据式(27)，双差相位模糊度L关于相对定轨滤波器状态量的偏导数H可以表示为：

其中：I_1×12中第i颗GPS/BDS对应的通道值为1，第j颗为GPS/BDS双差参考星对应的通道值为-1，其余通道为0，

表示第i颗GPS/BDS卫星对应的双差模糊度，

表示第i颗GPS/BDS卫星对应的单差模糊度，

表示第j颗GPS/BDS双差参考星对应的单差模糊度。

最后按式(15)，将单差模糊度参数约束为双差固定解，即将经过步骤d检验后的双差相位模糊度固定解N_DD作为观测量，依次对每颗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轨道插值输出所用到的插值端点信息

如图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轨道内插，可预报距离当前历元几分钟之后的未来时刻的高精度的绝对与相对轨道结果。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。