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

摘要

一种基于星载GNSS多天线的卫星自主定轨方法，本方法设计了扩展卡尔曼滤波器，充分利用多个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（Upper  triangular 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；

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

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

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

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

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

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

$\left(\begin{array}{c}\stackrel{·}{\Phi }(t,t_0)=\mathrm{A\Phi }(t,t_0)\\ \Phi (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\times 6}& {\mathrm{TQ}}_{\mathrm{rw}}\\ {\mathrm{TQ}}_{\mathrm{vr}}{T}^T& {\mathrm{TQ}}_{\mathrm{vv}}{T}^T& {\left[0\right]}_{6\times 6}& {\mathrm{TQ}}_{\mathrm{vw}}\\ {\left[0\right]}_{6\times 6}& {\left[0\right]}_{6\times 6}& Q_{\mathrm{other}}& {\left[0\right]}_{6\times 6}\\ Q_{\mathrm{wr}}{T}^T& Q_{\mathrm{wv}}{T}^T& {\left[0\right]}_{6\times 6}& Q_{\mathrm{ww}}\end{array}\right)$

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

$Q_{\mathrm{other}}=\left(\begin{array}{cccccc}\frac{S_g\Delta t^3}{3}+S_f\mathrm{\Delta t}& \frac{S_g\Delta t^2}{2}& 0& 0& 0& 0\\ 0& \frac{S_g\Delta t^3}{3}+S_f\mathrm{\Delta t}& \frac{S_g\Delta t^2}{2}& 0& 0& 0\\ 0& 0& \frac{S_g\Delta t^3}{3}+S_f\mathrm{\Delta t}& \frac{S_g\Delta t^2}{2}& 0& 0\\ 0& 0& \frac{S_g\Delta t^2}{2}& S_g\mathrm{\Delta t}& 0& 0\\ 0& 0& 0& 0& {\sigma }_{C_D}^2\mathrm{\Delta t}& 0\\ 0& 0& 0& 0& 0& {\sigma }_{C_R}^2\mathrm{\Delta t}\end{array}\right)$

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

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

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

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

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

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

$S=\left(\begin{array}{ccc}{\sigma }_R^2& 0& 0\\ 0& {\sigma }_T^2& 0\\ 0& 0& {\sigma }_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（全球定位系统/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)即可实现星载实时自主定轨。

（一）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计算如下：

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

其中，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{\overrightarrow{r}}{\left|\overrightarrow{r}\right|},e_y=\frac{-\overrightarrow{r}\times \overrightarrow{v}}{|\overrightarrow{r}\times \overrightarrow{v|}},e_x=e_y\times e_z$

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

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

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

${\overrightarrow{r}}_M={\overrightarrow{r}}_{p_i}-\overrightarrow{A}$

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

${\overrightarrow{r}}_{p_i}={\overrightarrow{r}}_M+\overrightarrow{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/GLONAS S/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{\delta }_{\mathrm{kj}}\\ E\left(V_k{V}_j^T\right)=R_k{\delta }_{\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}{{\rho }_k},\frac{Y_k-Y_g}{{\rho }_k},\frac{Z_k-Z_g}{{\rho }_k},0_{1\times 3},T_{1\times 4},0,0_{1\times 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}^{-}={\Phi }_{k,k-1}{P}_{k-1}^{+}{\Phi }_{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}}}{\delta },$其中，$\delta =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多天线的卫星自主定轨方案的精度评估与验证。验证结果表明：所 提出的方法，能够适用于复杂姿态机动情况，并且可以满足工程应用高精度 定轨需求。

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