CN103675861B - 一种基于星载gnss多天线的卫星自主定轨方法 - Google Patents

Abstract

一种基于星载GNSS多天线的卫星自主定轨方法，本方法设计了扩展卡尔曼滤波器，充分利用多个GNSS天线的实测伪距观测值对用高精度力学模型的轨道预报值进行实时滤波校正，得到高精度的卫星轨道信息。本发明能够解决复杂姿态机动情况下的高精度轨道确定问题。本发明实现的定轨结果精度高、稳定性好，实时性强，可满足低轨卫星高精度卫星轨道确定需求，可广泛应用于空间站、高分辨率对地观测卫星等航天任务的高精度定轨，具有广阔的推广应用前景。

Description

一种基于星载GNSS多天线的卫星自主定轨方法

技术领域

本发明涉及一种基于星载GNSS多天线的卫星自主定轨方法，能够解决复杂姿态机动情况下的高精度轨道确定问题，属于卫星的自主定轨技术领域。

背景技术

随着卫星应用技术的快速发展，卫星为了完成相应的科学任务，需要具有三轴姿态大角度敏捷机动能力。在敏捷模式下，GNSS接收机应能保证一直定位。然而，在大角度三轴姿态机动下，低轨卫星如果只使用一个GNSS天线，其主瓣不能保证一直指向天顶方向，甚至有可能指向地心，严重影响GNSS接收机正常工作。

另外一方面，整星上载荷、天线较多，如数传天线、中继天线等工作时需要调整天线指向对准中继卫星以便完成数据传输等任务，在调整天线指向时，可能会对GNSS天线产生物理遮挡，影响GNSS接收机收星情况。

综上可以看出，随着卫星应用技术的快速发展，仅利用GNSS单天线进行导航定位时会存在一定的适用局限性，已不能保证GNSS接收机输出高精度的导航定位结果，因此，需要研究适用范围更广的导航定位技术。

发明内容

本发明的技术解决问题是：克服现有技术的不足，提供一种基于星载GNSS多天线的卫星自主定轨方法，本发明能解决复杂姿态机动情况下的高精度轨道确定问题，相比单天线定轨，该方法适用范围更广，具有重大的工程实用价值和军事战略意义。

本发明的技术解决方案是：

一种基于星载GNSS多天线的卫星自主定轨方法包括步骤如下：

（1）自主定轨开始，系统初始化。

（2）获取星载GNSS多天线观测数据和广播星历，利用获取的观测数据进行单点定位测速解算得到卫星位置、速度结果和接收机钟差，进而观测数据历元时间减去接收机钟差得到滤波器时间；

（3）判断自主定轨标记是否为启动，若卫星自主定轨启动，则进入步骤（4）；若没有启动则将步骤（2）得到的卫星位置、速度结果由相应的天线相位中心转到卫星质心并初始化卡尔曼滤波状态量，将卡尔曼滤波状态量定义为数据1，并初始化卡尔曼滤波状态误差协方差阵，将方差阵定义为数据2，并置自主定轨标记启动进入步骤启动（4）；

（4）建立基于高精度轨道动力学模型的卫星运动状态方程，利用四阶龙格-库塔积分器对卫星运动状态方程进行积分，得到滤波器时间的卫星位置速度预报值和状态转移矩阵，分别将其定义为数据3和数据4；若为初次自主定轨，则卫星运动状态方程的积分初值使用步骤（2）中的卫星位置及速度结果；否则使用数据10；

（5）进行卡尔曼滤波的时间更新，计算系统动态噪声协方差阵，将其定义为数据5，利用卡尔曼滤波状态误差协方差阵和步骤（3）中得到的状态转移矩阵计算得到预测状态误差协方差阵，将其定义为数据6；若为初次自主定轨，则卡尔曼滤波状态误差协方差阵使用步骤（3）中的数据2，否则使用数据11；

（6）若步骤（2）中得到的GNSS多天线观测数据为第i个GNSS天线的观测数据，则将数据3从卫星质心处转换到第i个GNSS天线相位中心处，得到数据7；

（7）进行卡尔曼滤波的测量更新，更新步骤如下：

（a）首先计算滤波最优增益阵，建立以步骤（2）中的观测数据为观测量的观测方程，得到观测矩阵，将其定义为数据8，然后采用UD（Uppertriangular matrix-Diagonal matrix）分解滤波原理（一次卡尔曼滤波测量更新只处理一颗GNSS卫星的观测数据），利用观测噪声协方差阵（不同的GPS/GLONASS/BD2（全球定位系统/global position system,格洛纳斯，俄罗斯卫星导航系统/北斗2/beidou2）多系统设置不同的观测噪声协方差阵）、数据8和步骤（5）中得到的数据5和数据6计算滤波最优增益阵，将其定义为数据9；

（b）其次更新卡尔曼滤波状态量；利用步骤（a）中的数据9、步骤（6）中的数据7和步骤（2）中的观测数据更新卡尔曼滤波状态量，将其定义为数据10；

（c）最后更新状态协方差矩阵；利用步骤（a）中的观测噪声协方差阵和数据8以及步骤（5）中的数据6更新状态协方差矩阵，将其定义为数据11；

（8）卡尔曼滤波的测量更新过程结束后，将步骤（7）得到的数据10，由第i个天线相位中心处转换到卫星质心处，得到数据12，重复步骤（2）～(8)即可实现星载实时自主定轨。

所述步骤（2）中的观测数据为GPS、GLONASS和BD2观测数据，GPS、GLONASS和BD2观测数据以不等权方式参与测量更新，GPS、GLONASS和BD2导航系统参与测量更新的权重可以作为定轨可调参数通过星载上注方式调节改变，以获得GPS/GLONASS/BD2多系统最优数据融合结果。

所述步骤（6）和步骤（8）中的i为1、2、…N，N为GNSS天线个数。

所述步骤（6）从卫星质心处转换到GNSS天线相位中心充分考虑天线相位中心在卫星本体系下的安装位置和整星姿态，具体实施步骤如下：

（1）获取第i个天线相位中心在卫星本体系中的位置整星姿态和卫星在惯性系中的位置速度，所述位置速度为卫星位置速度预报值数据3；

（2）利用整星姿态计算卫星轨道系到卫星本体系转换矩阵M_bo；

（3）利用卫星在惯性系中的位置速度计算卫星轨道系到惯性系转换矩阵M_io；

（4）利用M^T _bo和M_io计算得到第i个天线相位中心在惯性系中的安装位置矢量

（5）利用卫星质心在惯性系中的位置矢量和矢量计算得到第i个天线相位中心在惯性系中的位置矢量从而完成数据3到数据7的转换。

所述步骤（8）第i个天线相位中心转换到卫星质心处充分考虑天线相位中心在卫星本体系下的安装位置和整星姿态，具体实施步骤如下：

（1）获取第i个天线相位中心在卫星本体系中的位置整星姿态和卫星在惯性系中的位置速度，所述位置速度来源于卡尔曼滤波状态量数据10；

（2）利用整星姿态计算卫星轨道系到卫星本体系转换矩阵M_bo；

（3）利用步骤（2）中的M^T _bo和位置速度计算卫星轨道系到惯性系转换矩阵M_io；

（4）利用M^T _bo和M_io计算得到第i个天线相位中心在惯性系中的安装位置矢量

（5）利用第i个天线相位中心在惯性系中的位置矢量和矢量计算得到卫星质心在惯性系中的位置矢量从而完成数据10到数据12的转换。

本发明与现有技术相比有益效果是：

(1)本发明以位置速度作为滤波状态量，不仅适用于经典开普勒轨道根数，也适用于小偏心率、低倾角、临界倾角轨道相应的无奇点轨道根数。

(2)本发明在测量更新时，充分考虑各GNSS天线的相位中心与卫星质心的转换，转换过程考虑天线相位中心在整星的安装位置和整星的姿态，相比单天线定轨，适用范围更广，能够适用于复杂姿态机动情况下的高精度轨道确定问题。

(3)本发明在测量更新时，GPS/GLONASS/BD2多系统观测数据以不等权方式参与测量更新，各导航系统权重可以通过星载上注方式调节改变，以获得GPS/GLONASS/BD2多系统最优数据融合结果，定轨精度更高。

附图说明

图1为本发明的方法流程图；

图2为本发明的卫星轨道系定义及天线相位偏心示意图；

图3为本发明的接收天线相位中心与卫星质心的相互转换示意图；

图4为采用GNSS信号源对本发明进行仿真验证示意图。

具体实施方式

下面结合附图对本发明的具体实施方式进行进一步的详细描述。

本发明提供的一种基于星载GNSS多天线实时定轨方法，该技术充分使用各天线GNSS观测量进行实时滤波处理，保证GNSS接收机输出高精度的导航定位结果，相比单天线定轨，该方法适用范围更广，具有重大的工程实用价值和军事战略意义。

如图1所示，本发明的具体步骤如下：

一种基于星载GNSS多天线的卫星自主定轨方法包括步骤如下：

（1）自主定轨开始，系统初始化（设定全局变量的初值）。

（2）获取星载GNSS多天线观测数据和广播星历，利用获取的观测数据进行单点定位测速解算得到卫星位置、速度结果和接收机钟差，进而观测数据历元时间减去接收机钟差得到滤波器时间；

（3）判断自主定轨标记是否为启动，若卫星自主定轨标记为启动，则进入步骤（4）；若没有启动则将步骤（2）得到的卫星位置、速度结果由相应的天线相位中心转到卫星质心并初始化卡尔曼滤波状态量，将卡尔曼滤波状态量定义为数据1，并初始化卡尔曼滤波状态误差协方差阵，将方差阵定义为数据2，并置自主定轨标记启动进入步骤启动（4）；

（4）建立基于高精度轨道动力学模型的卫星运动状态方程，利用四阶龙格-库塔积分器对卫星运动状态方程进行积分，得到滤波器时间的卫星位置速度预报值和状态转移矩阵，分别将其定义为数据3和数据4；若为初次自主定轨，则卫星运动状态方程的积分初值使用步骤（2）中的卫星位置及速度结果；否则使用数据10；

在地心惯性系，卫星运动状态方程是一个二阶微分方程，如下所示：

$\stackrel{\→}{\stackrel{\·\·}{r}}=-\frac{\mathrm{GM}\stackrel{\→}{r}}{r^3}+\stackrel{\→}{R}$

其中，为卫星在惯性系下的位置矢量，G为地球万有引力常数，M为地球质量，表示卫星受到的各种摄动加速度之和，包括保守力和非保守力。为了便于数值求解，可将其转化为一阶微分方程组来表示：

$\left\{\begin{array}{c}\stackrel{\→}{\stackrel{\·}{r}}=\stackrel{\→}{v}\\ \stackrel{\→}{\stackrel{\·}{v}}=-\frac{\mathrm{GM}\stackrel{\→}{r}}{r^3}+\stackrel{\→}{R}\end{array}\right.$

式中，为卫星运动速度矢量；已知t₀时刻的初始状态向量通过数值积分可得到下一时刻t的状态向量

与卫星运动状态方程相似，状态转移矩阵也可以用类似的方程给出

$\left\{\begin{array}{c}\stackrel{\·}{\mathrm{\Φ}}(t,t_0)=\mathrm{A\Φ}(t,t_0)\\ \mathrm{\Φ}(t_0,t_0)=I\end{array}\right.$

其中，A为摄动力模型的相关偏导数，I为单位矩阵，为状态转移矩阵的导数，状态转移矩阵和卫星运动状态方程一样，也是一个具有初值的常微分方程，为了能比较精确计算这些分量的状态转移矩阵，通常将其与卫星状态参数（包括卫星位置和速度共6维）一起，通过数值积分进行求解。

（5）进行卡尔曼滤波的时间更新，计算系统动态噪声协方差阵，将其定义为数据5，利用卡尔曼滤波状态误差协方差阵和步骤（3）中得到的状态转移矩阵计算得到预测状态误差协方差阵，将其定义为数据6；若为初次自主定轨，则卡尔曼滤波状态误差协方差阵使用步骤（3）中的数据2，否则使用数据11；

根据卡尔曼滤波的状态方程和协方差传播理论，状态方程的系统动态噪声协方差矩阵（亦称为过程噪声协方差矩阵）可以表示为

$Q=\left[\begin{array}{cccc}{\mathrm{TQ}}_{\mathrm{rr}}{T}^T& {\mathrm{TQ}}_{\mathrm{rv}}{T}^T& {\left[0\right]}_{6\×6}& {\mathrm{TQ}}_{\mathrm{rw}}\\ {\mathrm{TQ}}_{\mathrm{vr}}{T}^T& {\mathrm{TQ}}_{\mathrm{vv}}{T}^T& {\left[0\right]}_{6\×6}& {\mathrm{TQ}}_{\mathrm{vw}}\\ {\left[0\right]}_{6\×6}& {\left[0\right]}_{6\×6}& Q_{\mathrm{other}}& {\left[0\right]}_{6\×6}\\ Q_{\mathrm{wr}}{T}^T& Q_{\mathrm{wv}}{T}^T& {\left[0\right]}_{6\×6}& Q_{\mathrm{ww}}\end{array}\right]$

其中，T为卫星轨道径向、切向和法向坐标系(RTN坐标系)到空间直角坐标系的转换矩阵。

$Q_{\mathrm{other}}=\left[\begin{array}{cccccc}\frac{S_g\mathrm{\Δ}{t}^3}{3}+S_f\mathrm{\Δt}& \frac{S_g\mathrm{\Δ}{t}^2}{2}& 0& 0& 0& 0\\ 0& \frac{S_g\mathrm{\Δ}{t}^3}{3}+S_f\mathrm{\Δt}& \frac{S_g\mathrm{\Δ}{t}^2}{2}& 0& 0& 0\\ 0& 0& \frac{S_g\mathrm{\Δ}{t}^3}{3}+S_f\mathrm{\Δt}& \frac{S_g\mathrm{\Δ}{t}^2}{2}& 0& 0\\ 0& 0& \frac{S_g\mathrm{\Δ}{t}^2}{2}& S_g\mathrm{\Δt}& 0& 0\\ 0& 0& 0& 0& {\mathrm{\σ}}_{C_D}^2\mathrm{\Δt}& 0\\ 0& 0& 0& 0& 0& {\mathrm{\σ}}_{C_R}^2\mathrm{\Δt}\end{array}\right]$

$Q_{\mathrm{rr}}=[\frac{\mathrm{\Δ}{t}^3}{3}-\mathrm{\Δ}{t}^2\mathrm{\τ}+\frac{1-e^{-2\mathrm{\Δt}/\mathrm{\τ}}}{2}{\mathrm{\τ}}^3+\mathrm{\Δt}(1-e^{-\mathrm{\Δt}/\mathrm{\τ}}){\mathrm{\τ}}^2]{\mathrm{\τ}}^2\·S$

$Q_{\mathrm{rv}}=Q_{\mathrm{vr}}=[\frac{\mathrm{\Δ}{t}^2}{2}+(\frac{1+e^{-2\mathrm{\Δt}/\mathrm{\τ}}}{2}-e^{-\mathrm{\Δt}/\mathrm{\τ}}){\mathrm{\τ}}^2-\mathrm{\Δt}(1-e^{-\mathrm{\Δt}/\mathrm{\τ}})\mathrm{\τ}]{\mathrm{\τ}}^2\·S$

$Q_{\mathrm{vv}}=(\mathrm{\Δt}+\frac{4e^{-\mathrm{\Δt}/\mathrm{\τ}}-e^{-2\mathrm{\Δt}/\mathrm{\τ}}-3}{2}\mathrm{\τ}){\mathrm{\τ}}^2\·S$

$Q_{\mathrm{rw}}=Q_{\mathrm{wr}}=(\frac{1-e^{-2\mathrm{\Δt}/\mathrm{\τ}}}{2}\mathrm{\τ}-\mathrm{\Δt}{e}^{-\mathrm{\Δt}/\mathrm{\τ}}){\mathrm{\τ}}^2\·S$

$Q_{\mathrm{vw}}=Q_{\mathrm{wv}}=(\frac{1+e^{-2}\mathrm{\Δt}/\mathrm{\τ}}{2}-e^{-\mathrm{\Δt}/\mathrm{\τ}}){\mathrm{\τ}}^2\·S$

$Q_{\mathrm{ww}}=\frac{1-e^{-2\mathrm{\Δt}/\mathrm{\τ}}}{2}\mathrm{\τ}\·S$

$S=\left[\begin{array}{ccc}{\mathrm{\σ}}_R^2& 0& 0\\ 0& {\mathrm{\σ}}_T^2& 0\\ 0& 0& {\mathrm{\σ}}_N^2\end{array}\right]$

其中，S_g、S_f为接收机钟差模型参数，Δt为时间步长，τ为一阶高斯马尔可夫模型的相关时间，σ_R,σ_T,σ_N分别为卫星轨道在径向、切向和法向的位置方差。

（6）若步骤（2）中得到的GNSS多天线观测数据为第i个GNSS天线的观测数据，则将数据3从卫星质心处转换到第i个GNSS天线相位中心处，得到数据7；

（7）进行卡尔曼滤波的测量更新，更新步骤如下：

（a）首先计算滤波最优增益阵，建立以步骤（2）中的观测数据为观测量的观测方程，得到观测矩阵，将其定义为数据8，然后采用UD分解滤波原理（一次卡尔曼滤波测量更新只处理一颗GNSS卫星的观测数据），利用观测噪声协方差阵（不同的GPS/GLONASS/BD2（全球定位系统/globalposition system,格洛纳斯，俄罗斯卫星导航系统/北斗2/beidou2）多系统设置不同的观测噪声协方差阵）、数据8和步骤（5）中得到的数据5和数据6计算滤波最优增益阵，将其定义为数据9；

（b）其次更新卡尔曼滤波状态量；利用步骤（a）中的数据9、步骤（6）中的数据7和步骤（2）中的观测数据更新卡尔曼滤波状态量，将其定义为数据10；

（c）最后更新状态协方差矩阵；利用步骤（a）中的观测噪声协方差阵和数据8以及步骤（5）中的数据6更新状态协方差矩阵，将其定义为数据11；

（8）卡尔曼滤波的测量更新过程结束后，将步骤（7）得到的数据10，由第i个天线相位中心处转换到卫星质心处，得到数据12，重复步骤（2）～(8)即可实现星载实时自主定轨。

（一）GNSS天线相位中心到卫星质心的相互转换

GNSS天线相位中心到卫星质心的相互转换需要充分考虑天线相位中心在卫星本体系下的安装位置和整星姿态保证转换的精确度。

（a）从卫星质心处转换到GNSS天线相位中心充分考虑天线相位中心在星固系下的安装位置和整星姿态，具体实施步骤如下：

（1）获取第i个天线相位中心在卫星本体系中的位置整星姿态和卫星在惯性系中的位置速度，所述位置速度为卫星位置速度预报值数据3；

（2）利用整星姿态计算卫星轨道系到卫星本体系转换矩阵M_bo；

（3）利用卫星在惯性系中的位置速度计算卫星轨道系到惯性系转换矩阵M_io；

（4）利用M_bo和M_io计算得到第i个天线相位中心在惯性系中的安装位置矢量

（5）利用卫星质心在惯性系中的位置矢量和矢量计算得到第i个天线相位中心在惯性系中的位置矢量从而完成数据3到数据7的转换。

（b）第i个天线相位中心转换到卫星质心处充分考虑天线相位中心在星固系下的安装位置和整星姿态，具体实施步骤如下：

（1）获取第i个天线相位中心在卫星本体系中的位置整星姿态和卫星在惯性系中的位置速度，所述位置速度来源于卡尔曼滤波状态量数据10；

（2）利用整星姿态计算卫星轨道系到卫星本体系转换矩阵M_bo；

（3）利用步骤（2）中的M_ob和位置速度计算卫星轨道系到惯性系转换矩阵M_io；

（4）利用M_bo和M_io计算得到第i个天线相位中心在惯性系中的安装位置矢量

（5）利用第i个天线相位中心在惯性系中的位置矢量和矢量计算得到卫星质心在惯性系中的位置矢量从而完成数据10到数据12的转换。

（c）下面以一个实施例具体说明相互转换过程

如图2所示，卫星的轨道系定义：坐标系原点为卫星质心M，Z轴指向地球质心O，Y轴指向卫星轨道负法向，X轴与Y、Z轴成右手坐标系。

设M_bo为卫星坐标系到轨道系转换矩阵，则其可以由整星的轨道系欧拉姿态角及姿态转序即可求得，例如，假设轨道系欧拉姿态角为α、β、γ，姿态转序为123，则M_bo计算如下：

$\begin{array}{c}{M}_{\mathrm{bo}}=R_3\left(\mathrm{\γ}\right)\·R_2\left(\mathrm{\β}\right)\·R_1\left(\mathrm{\α}\right)\\ =\left(\begin{array}{ccc}\cos \mathrm{\γ}& \sin \mathrm{\γ}& 0\\ -\sin & \cos \mathrm{\γ}& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{ccc}\cos \mathrm{\β}& 0& -\sin \mathrm{\β}\\ 0& 1& 0\\ \sin \mathrm{\β}& 0& \cos \mathrm{\β}\end{array}\right)\left(\begin{array}{ccc}1& 0& 0\\ 0& \cos \mathrm{\α}& \sin \mathrm{\α}\\ 0& -\sin \mathrm{\α}& \cos \mathrm{\α}\end{array}\right)\end{array}$

其中，R₃(γ)为绕z轴的旋转矩阵，R₂(β)为绕y轴的旋转矩阵，R₁(α)为为绕x轴的旋转矩阵；

假设卫星质心M点的位置和速度分别为和M_io为轨道系到地心惯性系转换矩阵，M_io=[e_x;e_y;e_z]^T，e_x,e_y,e_z是卫星轨道系三个轴方向在地心惯性系中的单位矢量。根据卫星轨道系的定义，该单位矢量可表示为：

$e_z=-\frac{\stackrel{\→}{r}}{\left|\stackrel{\→}{r}\right|},e_y=\frac{-\stackrel{\→}{r}\×\stackrel{\→}{v}}{|\stackrel{\→}{r}\×\stackrel{\→}{v|}},e_x=e_y\×e_z$

如图3所示，设星载GNSS接收机第i个天线相位中心P_i在卫星本体坐标系内的坐标为则它在地心惯性系O-XiYiZi中的坐标表示为：

$\stackrel{\→}{A}=[a_1,a_2,a_3]\·{M^T}_{\mathrm{bo}}\·M_{\mathrm{io}}$

如果给定卫星第i个天线相位中心的位置则卫星质量中心的位置可表示为

${\stackrel{\→}{r}}_M={\stackrel{\→}{r}}_{p_i}-\stackrel{\→}{A}$

如果给定卫星的质量中心的位置则卫星第i个天线相位中心的位置可表示为

${\stackrel{\→}{r}}_{p_i}={\stackrel{\→}{r}}_M+\stackrel{\→}{A}$

（二）卡尔曼滤波器设计

下面以一个实施例具体说明卡尔曼滤波器滤波过程

（a）卡尔曼滤波状态方程

设计的自主定轨扩展卡尔曼滤波的状态方程为

X_k=Φ_k,k-1X_k-1+W_k

式中，X_k为滤波器状态量（被估计量）；Φ_k,k-1为状态转移矩阵；W_k为系统噪声矩阵；

选择自主定轨滤波状态量为：

$X={[x,y,z,\stackrel{\·}{x},\stackrel{\·}{y},\stackrel{\·}{z},b_G,b_R,b_C,{\stackrel{\·}{b}}_u,C_R,w_R,w_T,w_N]}^T$

其中，为J2000惯性系下卫星三轴位置速度；

b_G=cδt_G,b_R=cδt_R,b_C=cδt_C,δt_G,δt_R,δt_C 分别为星载接收机GPS/GLONASS/BD2的钟差和频差，c为真空中的光速；C_D为大气阻力系数；C_R为太阳光压系数；w_R,w_T,w_N为卫星沿径向、切向、法向三个方向的补偿加速度分量。

（b）卡尔曼滤波观测方程

设计的自主定轨扩展卡尔曼滤波的状态方程为

Z_k=H_kX_k+V_k

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

$\left.\begin{array}{c}E\left(W_k\right)=0\\ E\left(V_k\right)=0\\ E\left(W_k{W}_j^T\right)=Q_k{\mathrm{\δ}}_{\mathrm{kj}}\\ E\left(V_k{V}_j^T\right)=R_k{\mathrm{\δ}}_{\mathrm{kj}}\\ E\left(W_k{V}_j^T\right)=0\end{array}\right\}$

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

选取GNSS观测数据为伪距观测值，有：

$H=[\frac{X_k-X_g}{{\mathrm{\ρ}}_k},\frac{Y_k-Y_g}{{\mathrm{\ρ}}_k},\frac{Z_k-Z_g}{{\mathrm{\ρ}}_k},0_{1\×3},T_{1\×4},0,0_{1\×3}]$

其中，为用户卫星的地固系下三维坐标矢量，由J2000惯性系下滤波状态量卫星三轴位置矢量经坐标转换得到，为信号发射时刻的GNSS卫星在地固系下三维坐标矢量，为卫星与GNSS卫星的距离。当为GPS卫星时，T_1×4=[1,0,0,0]，当为GLONASS卫星时，T_1×4=[0,1,0,0]，当为BD2卫星时，T_1×4=[0,0,1,0]。

（c）卡尔曼滤波时间更新

$P_{k|k-1}^{-}={\mathrm{\Φ}}_{k,k-1}{P}_{k-1}^{+}{\mathrm{\Φ}}_{k,k-1}^T+Q_k$

式中，Φ_k,k-1为状态转移矩阵，Φ_k,k-1为Φ_k,k-1的转置矩阵，Q_k为系统动态噪声协方差阵，和分别为时间更新前后的状态误差协方差矩阵，若为初次自主定轨，则使用步骤（3）中的数据2，否则使用数据11，即为后面对所有GNSS卫星进行测量更新完成后的状态误差协方差矩阵

（d）卡尔曼滤波测量更新

卡尔曼滤波的测量更新采用类似于UD分解滤波原理，因此测量更新针对标量来进行计算。对于星载GNSS伪距观测数据来说，一次卡尔曼滤波测量更新只处理一颗GNSS卫星的观测数据，即假设用户卫星共总有N个GNSS天线，第i个GNSS天线中有Mi颗GNSS卫星，对于第i个GNSS天线中的第j颗GNSS卫星，（i=1、2、…N，j=1、2、…M_i），则其测量更新过程如下：

1）计算滤波最优增益阵

$K_k^j=\frac{P_k^{j-1}{H}_k^{\mathrm{jT}}}{\mathrm{\δ}},$ 其中， $\mathrm{\δ}=H_k^j{P}_k^{j-1}{H}_k^{\mathrm{jT}}+{R_k}^2$

式中，为第j颗GNSS卫星对应的卡尔曼滤波最优增益阵，为第j颗GNSS卫星对应的观测矩阵；Rk为观测噪声协方差阵。第j颗GNSS卫星测量更新前的状态误差协方差矩阵；

对于GPS/GLONASS/BD2不同的导航系统，伪距观测值的精度指标也不同，其对应的测噪声协方差阵可设为不同，设其权重比为R_gps:R_glo:R_bd2，可以作为定轨可调参数通过星载上注方式调节改变，以获得GPS/GLONASS/BD2多系统最优数据融合结果。

2）更新卡尔曼滤波状态量

$X_k^{j+}=X_{k|}^{j-}+k_k^j\·y_k^j,$ 其中， $y_k^i={Z_0}^j-{Z_k}^j={Z_0}^j-H_k^j{X}_k^{j-}$

式中，为第j颗GNSS卫星对应的卡尔曼滤波最优增益阵，为观测向量，为第j颗GNSS卫星的实测伪距观测量值与计算的伪距观测量值之差，即新息，Z₀ ^j为第j颗GNSS卫星的接收机实测伪距观测量值，Z_k ^j为第j颗GNSS卫星的计算的伪距观测量值，和为第j颗GNSS卫星测量更新前后的卡尔曼滤波状态向量；

3）更新状态协方差矩阵

$P_k^{j+}=(I-K_k^j{H}_k^j)P_k^{j-}$

式中，I为单位阵，为第j颗GNSS卫星对应的卡尔曼滤波最优增益阵，和为第j颗GNSS卫星测量更新前后的状态误差协方差矩阵；

（三）仿真验证方案设计

如图4所示，为了验证本发明方法的正确性，采用GNSS信号源进行验证，验证方案步骤如下：

（a）GNSS接收机接收GNSS仿真源输出的GNSS射频信号，输出GNSS接收机的广播星历和原始观测数据，并通过GNSS地检设备将广播星历和原始观测数据存到PC计算机上；

（b）在PC机上，使用基于星载GNSS多天线的卫星自主定轨方案对步骤（1）中得到的GNSS接收机输出的广播星历和原始观测数据进行模拟自主定轨数据处理，输出高精度的定轨计算结果；

（c）将步骤（2）中得到的模拟自主定轨数据处理的解算结果与GNSS信号仿真器输出的理论轨道进行比对，完成用GNSS仿真实测数据对基于星载GNSS多天线的卫星自主定轨方案的精度评估与验证。验证结果表明：所提出的方法，能够适用于复杂姿态机动情况，并且可以满足工程应用高精度定轨需求。

本发明说明书中未作详细描述的内容属于本领域技术人员的公知技术。

Claims (5)

1.一种基于星载GNSS多天线的卫星自主定轨方法，其特征在于步骤如下：

(1)自主定轨开始，系统初始化；

(2)获取星载GNSS多天线观测数据和广播星历，利用获取的观测数据进行单点定位测速解算得到卫星位置、速度结果和接收机钟差，进而利用观测数据中的历元时间减去接收机钟差得到滤波器时间；

(3)判断自主定轨标记是否为启动，若卫星自主定轨标记为启动，则进入步骤(4)；若没有启动则卫星为初次定轨将步骤(2)得到的卫星位置、速度结果由相应的天线相位中心转到卫星质心并初始化卡尔曼滤波状态量，将卡尔曼滤波状态量定义为数据1，并初始化卡尔曼滤波状态误差协方差阵，将方差阵定义为数据2，并置自主定轨标记启动进入步骤启动(4)；

(4)建立基于高精度轨道动力学模型的卫星运动状态方程，利用四阶龙格-库塔积分器对卫星运动状态方程进行积分，得到滤波器时间的卫星位置速度预报值和状态转移矩阵，分别定义为数据3和数据4；若为初次自主定轨，则卫星运动状态方程的积分初值使用步骤(2)中的卫星位置及速度结果；否则使用数据10；

(5)进行卡尔曼滤波的时间更新，计算系统动态噪声协方差阵，将系统动态噪声协方差阵定义为数据5，利用卡尔曼滤波状态误差协方差阵和步骤(4)中得到的状态转移矩阵计算得到预测状态误差协方差阵，将预测状态误差协方差阵定义为数据6；若为初次自主定轨，则卡尔曼滤波状态误差协方差阵使用步骤(3)中的数据2，否则使用数据11；

(6)若步骤(2)中得到的GNSS多天线观测数据为第i个GNSS天线的观测数据，则将数据3从卫星质心处转换到第i个GNSS天线相位中心处，得到数据7；

(7)进行卡尔曼滤波的测量更新，更新步骤如下：

(a)首先计算滤波最优增益阵，建立以步骤(2)中的观测数据为观测量的观测方程，得到观测矩阵，将观测矩阵定义为数据8，然后采用UD分解滤波原理，利用观测噪声协方差阵、数据8和步骤(5)中得到的数据5和数据6计算滤波最优增益阵，将滤波最优增益阵定义为数据9；

(b)其次更新卡尔曼滤波状态量；利用步骤(a)中的数据9、步骤(6)中的数据7和步骤(2)中的观测数据更新卡尔曼滤波状态量，将卡尔曼滤波状态量定义为数据10；

(c)最后更新状态协方差矩阵；利用步骤(a)中的观测噪声协方差阵和数据8以及步骤(5)中的数据6更新卡尔曼滤波状态误差协方差阵，将状态协方差矩阵定义为数据11；

(8)卡尔曼滤波的测量更新过程结束后，将步骤(7)得到的数据10，由第i个天线相位中心处转换到卫星质心处，得到数据12，重复步骤(2)～(8)即可实现星载实时自主定轨。

2.根据权利要求1所述一种基于星载GNSS多天线的卫星自主定轨方法，其特征在于：所述步骤(2)中的观测数据为GPS、GLONASS和BD2观测数据，GPS、GLONASS和BD2观测数据以不等权方式参与测量更新，GPS、GLONASS和BD2导航系统参与测量更新的权重可以作为定轨可调参数通过星载上注方式调节改变，以获得GPS/GLONASS/BD2多系统最优数据融合结果。

3.根据权利要求1所述一种基于星载GNSS多天线的卫星自主定轨方法，其特征在于：所述步骤(6)和步骤(8)中的i为1、2、…N，N为GNSS天线个数。

4.根据权利要求1所述一种基于星载GNSS多天线的卫星自主定轨方法，其特征在于：所述步骤(6)从卫星质心处转换到GNSS天线相位中心充分考虑天线相位中心在卫星本体系下的安装位置和整星姿态，具体实施步骤如下：

(1)获取第i个天线相位中心在卫星本体系中的位置整星姿态和卫星在惯性系中的位置速度，所述位置速度为卫星位置速度预报值数据3；

(2)利用整星姿态计算卫星轨道系到卫星本体系转换矩阵M_bo；

(3)利用卫星在惯性系中的位置速度计算卫星轨道系到惯性系转换矩阵M_io；

(4)利用M_bo和M_io计算得到第i个天线相位中心在惯性系中的安装位置矢量

(5)利用卫星质心在惯性系中的位置矢量和矢量计算得到第i个天线相位中心在惯性系中的位置矢量从而完成数据3到数据7的转换。

5.根据权利要求1所述一种基于星载GNSS多天线的卫星自主定轨方法，其特征在于：所述步骤(8)第i个天线相位中心转换到卫星质心处充分考虑天线相位中心在卫星本体系下的安装位置和整星姿态，具体实施步骤如下：

(1)获取第i个天线相位中心在卫星本体系中的位置整星姿态和卫星在惯性系中的位置速度，所述位置速度来源于卡尔曼滤波状态量数据10；

(2)利用整星姿态计算卫星轨道系到卫星本体系转换矩阵M_bo；

(3)利用步骤(2)中的M^T _bo和位置速度计算卫星轨道系到惯性系转换矩阵M_io；

(4)利用M^T _bo和M_io计算得到第i个天线相位中心在惯性系中的安装位置矢量

(5)利用第i个天线相位中心在惯性系中的位置矢量和矢量计算得到卫星质心在惯性系中的位置矢量从而完成数据10到数据12的转换。