CN106338753A - 一种基于地面站/星间链路/gnss联合测量的地球同步轨道星座定轨方法 - Google Patents

Abstract

本发明公开了一种基于地面站/星间链路/GNSS联合测量的地球同步轨道星座定轨方法，具体包括以下几个步骤：步骤一：在集中式结构下，建立GSO星座定轨系统的状态方程；步骤二：在集中式结构下，建立GSO星座定轨系统的测量方程；步骤三：联合优化用于定轨的地面站和GNSS测量组合；步骤四：根据上述建立的基于地面站/星间链路/GNSS的GSO星座联合定轨系统模型，确定实现轨道参数估计的滤波方法；步骤五：在集中式结构下，基于地面站/星间链路/GNSS的GSO星座联合定轨方法的具体实现。本发明基于精度因子法建立了地面站和GNSS测量的联合分布优选策略，在保证定轨精度的同时降低计算量。

Description

一种基于地面站/星间链路/GNSS联合测量的地球同步轨道星 座定轨方法

技术领域

本发明属于地球同步轨道卫星星座的定轨领域，具体来说，是一种通过有效融合利用地面站/星间链路/GNSS三种测量方式分别获取的星-地、星-星以及GNSS的测量信息，实现地球同步轨道卫星星座联合定轨方法。

背景技术

随着航天科技的快速发展，地球同步轨道卫星(Geosynchronous OrbitSatellite，GSO)星座在导航定位、授时、通讯和侦查等领域得到广泛应用。高精度的卫星定轨技术是星座实现整体功能和性能的重要前提保障。地球静止轨道(Geostationary EarthOrbits，GEO)卫星和倾斜地球同步轨道(Inclined Geosynchronous Orbit，IGSO)卫星构成了GSO星座。该星座内的卫星均属于高轨卫星，且GEO具有静地特性；IGSO的星下点轨迹是以赤道为对称轴的“8”字形，“8”字覆盖的区域与轨道倾角有关。与其他轨道类型卫星相比，GSO卫星精密轨道确定存在较大的困难。

采用基于地面站测量技术对GSO卫星定轨存在一定的局限性，(1)GEO的“静地”特性及小倾角的IGSO卫星星下点轨迹区域分布范围小的特点，使之与地面站的相对位置变化不明显，星站间的观测几何变化很小，即卫星相对于地面站的动力学约束信息较弱，此时多普勒跟踪作用不明显，增加观测时间带来的信息量也有限；(2)当仅能选用国境内的地面站用于卫星测控时，由于布设范围相对较小，GSO的高轨特性使卫星与地面站之间构成的观测几何结构强度相当差，将导致在估计卫星轨道参数过程中，观测方程可能会呈现病态，影响估计精度。(3)为了保持GEO卫星对地同步，需对其频繁地实施机动控制,也给GEO卫星精密轨道的确定和预报带来较大麻烦。因此，受到地面站布设地域、设备条件等限制，通过提高测量精度和改善观测几何强度(合理进行地面站选取)，来提高基于地面站的卫星定轨精度能力是有限的。因此，需要融合其他测量技术来保证卫星星座的实时定轨精度，同时降低对地面站的依赖性。

美国已经成功将GPS(Global Positioning System)用于低轨卫星的定轨中。但对于高轨卫星来说，高轨处的GPS信号强度弱，且GPS卫星可见性差。已有研究结果显示，即使将接收机灵敏度提高至28dB-Hz，GSO卫星在相当一部分时间内仍无法观测到定轨所需的至少4颗以上的GPS卫星。显然，GPS几何定轨法将难以在GSO卫星全轨道适用。我国的北斗卫导系统采用GEO+IGSO+MEO的星座构型，目前正在建设完善中，这使高轨卫星对GNSS(GlobalNavigation Satellite System)卫星的可见性得到有效改善。因此，采用基于北斗/GPS的双模GNSS技术可以辅助高轨卫星实时定轨，降低对地面站测量的依赖性。在利用地面站和星载GNSS接收机提供的测量信息有效实现单星绝对定轨的基础上，再融合利用星间链路获取的星间相对测量信息，对于改善星间相对定位的性能及卫星星座的自主定轨能力具有重要作用。

目前，国内关于高轨卫星星座综合利用星-地、星-星及GNSS测量实现定轨的研究还很少。

发明内容

本发明为了解决上述问题，提出一种基于地面站/星间链路/GNSS联合测量的地球同步轨道星座定轨方法，该方法利用地面站/星间链路/GNSS三种测量方式分别获取的星-地、星-星以及GNSS的测量信息，在基于所建立的测量信息优选判别策略实现测量信息有效性判别的基础上，利用扩展卡尔曼滤波(EKF，Extended Kalman Filter)等非线性状态估计法，采用集中式结构融合有效测量信息，实现GSO星座中各个卫星轨道参数的解算。该星座定轨方法首先采用集中式结构建立GSO星座定轨的系统模型，其中状态方程主要基于地心惯性直角坐标系下的卫星轨道动力学模型来建立，测量方程依据地面站/星间链路/GNSS的测量模型来建立；然后基于所建立的测量信息优选判别策略来实现测量信息有效性判别，基于最大单位矢端多面体体积法实现地面站测量信息的优选，基于加权精度因子法同时考虑观测几何结构和测量精度实现星载GNSS测量信息的实时优选；最后，采用EKF等一类非线性状态方法对有效测量信息进行融合实现GSO多星星座定轨。以4颗GSO卫星构成的星座为例验证方法的有效性。

一种基于地面站/星间链路/GNSS联合测量的地球同步轨道星座定轨方法，具体包括以下几个步骤：

步骤一：在集中式结构下，建立GSO星座定轨系统的状态方程；

步骤二：在集中式结构下，建立GSO星座定轨系统的测量方程；

步骤三：定轨测量组合的优化策略；

步骤四：根据上述建立的基于地面站/星间链路/GNSS的GSO星座联合定轨系统模型，确定实现轨道参数估计的滤波方法；

步骤五：在集中式结构下，基于地面站/星间链路/GNSS的GSO星座联合定轨方法的具体实现。

本发明的优点在于：

(1)基于精度因子法建立了地面站分布优选策略，在此基础上优选地面站测量信息，在保证定轨精度的同时降低计算量；

(2)融合了星间链路提供的相对测量信息，改善了星间相对定位性能，在一定程度上降低对地面站测量的依赖性；

(3)融合了利用星载GNSS接收机获得GPS/BD双系统测量提供的辅助信息，提高了定轨精度，并降低星座对地面站的依赖，提高了自主运行能力；

(4)建立GNSS和地面站联合测量信息的优选策略，同时考虑待测卫星与GNSS卫星和地面站的几何构型和测量精度，基于加权精度因子法的实现了GNSS和地面站测量的联合优选，有效利用地面站和星载GNSS接收机的测量信息以确保定轨精度，同时降低系统计算量。

附图说明

图1为基于地面站/星间链路/GNSS的GSO星座定轨系统构成示意图；

图2为基于地面站/星间链路/GNSS联合测量的GSO星座定轨方法结构图；

图3为单位矢端多面体体积示意图；

图4为集中式信息融合算法流程图。

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

本发明是一种基于地面站/星间链路/GNSS的GSO星座定轨方法，基于地面站/星间链路/GNSS联合测量的GSO星座定轨系统构成参见图1，图中显示该星座定轨系统由待测GSO星座、分布在地球表面的地面站以及GNSS星座构成，其中待测GSO星座由m颗高轨同步卫星组成，编号依次为1,2,…,j,…m；n个地面站用于监测GSO星座运行情况，编号依次为1,2,…,i,…n；GSO星座可见的GNSS卫星有l个，编号依次为1,2,…,s,…,l。一种基于地面站/星间链路/GNSS联合测量的地球同步轨道星座定轨方法的基本实现过程参见图2，图中显示该定轨方法的基本实现过程主要包括测量信息获取、星座轨道参数估计。其中获取的测量信息包括地面站的测量信息、星载GNSS接收机提供的可见星测量信息和星间链路提供的星间相对测量信息；在进行星座轨道参数估计时，①建立基于卫星简化轨道动力学模型的状态方程，②根据三种测量方式的测量原理建立测量方程，③利用由状态方程和测量方程构成的系统模型，采用非线性状态估计方法，实现轨道参数的估计；最后，输出估计得到的星座中卫星的位置和速度以及星间相对位置和速度。在地心惯性坐标系下，采用直角坐标描述法(即用位置和速度分量来表示卫星当前的运行状态)，利用EKF等非线性估计方法，结合轨道动力学方程和测量信息建立系统模型，实现GSO星座定轨，具体步骤为：

步骤一：在集中式结构下，建立GSO星座定轨系统的状态方程；

准确的系统模型是保障星座运行状态参数估计精度的主要因素之一。本发明基于卫星轨道动力学模型来建立系统状态方程。

卫星受到地球质心引力和摄动力两类作用力。为了满足计算高效和高精度要求，根据摄动力影响大小来选择主要摄动项建模。由于GSO轨道高度较高，在两天的弧度上，地球非球形对称带来的摄动影响可达10⁴米级，其中3阶以上引力场摄动影响与J2项相比小很多，可以忽略；日月引力摄动影响可达10³米级，太阳光压摄动影响可达10²米级；而潮汐、相对论效应和反照光压摄动等其他摄动的影响在10^-1米级以下；同时高轨处大气密度甚微，大气阻力摄动影响可以忽略。因此，本发明基于地球质心引力以及J2项、日月引力和太阳光压三种主要摄动力构成的轨道动力学模型来建立系统状态方程。

在集中式信息融合方法中，状态向量除了考虑所有卫星运行状态之外，还需将星载GNSS接收机时钟的钟差和频差也作为状态变量扩维到状态方程中。由于其数量级很小，直接引入计算很可能被计算机舍入误差湮没，因此，根据钟差和频差与距离之间的关系，将它们转化为等效测距误差δl_tu和距离变化率误差δl_tru。并建立星载GNSS接收机等效时钟误差的状态方程如下：

$\begin{array}{c}\mathrm{\δ}{\stackrel{\·}{l}}_{tu}={\mathrm{\δl}}_{tru}+w_{tu}\\ \mathrm{\δ}{\stackrel{\·}{l}}_{tru}=-\frac{1}{{\mathrm{\τ}}_{tru}}{\mathrm{\δl}}_{tru}+w_{tru}\end{array}---\left(1\right)$

式中，τ_tru为相关时间常数，w_tu和w_tru为零均值高斯白噪声。

在地心直角惯性坐标系下，结合轨道动力学模型和接收机等效钟差模型可建立对星座中待测卫星j进行定轨的状态方程如下：

${\stackrel{\·}{X}}_j\left(t\right)=f_j\[X_j\left(t\right),t\]+W_j\left(t\right)---\left(2\right)$

式中，待测卫星j对应的状态向量为这里[x_j,y_j,z_j]^T和分别为该星的位置矢量和速度矢量，为该星载GNSS接收机的等效钟差和频差。f_j为系统状态函数，W_j(t)为系统噪声，满足均值为零，方差为Q_j(t)的高斯白噪声。式(2)可以进一步详细写为：

$\left[\begin{array}{c}{\stackrel{\·}{x}}_j\\ {\stackrel{\·}{y}}_j\\ {\stackrel{\·}{z}}_j\\ {\stackrel{\·}{v}}_{x_j}\\ {\stackrel{\·}{v}}_{y_j}\\ {\stackrel{\·}{v}}_{z_j}\\ \mathrm{\δ}{\stackrel{\·}{l}}_{tu,j}\\ \mathrm{\δ}{\stackrel{\·}{l}}_{tru,j}\end{array}\right]=\left[\begin{array}{c}{v}_{x_j}\\ v_{y_j}\\ v_{z_j}\\ a_{1x_j}+a_{2x_j}+a_{3x_j}\\ a_{1y_j}+a_{2y_j}+a_{3y_j}\\ a_{1z_j}+a_{2z_j}+a_{3z_j}\\ {\mathrm{\δl}}_{tru,j}\\ -\frac{1}{{\mathrm{\τ}}_{tru}}{\mathrm{\δl}}_{tru,j}\end{array}\right]+\left[\begin{array}{c}{w}_{x_j}\\ w_{y_j}\\ w_{z_j}\\ w_{v_{x_j}}\\ w_{v_{y_j}}\\ w_{v_{z_j}}\\ w_{tu,j}\\ w_{tru,j}\end{array}\right]$

其中，和为速度噪声分量，和是未建模的摄动加速度分量，w_tu,j和w_tru,j分别为等效钟差和频差噪声；设待测卫星j受到的地球引力、日月引力和太阳光压引起的加速度分别为a_1,j、a_2,j和a_3，j，具体表达式如下：

(1)本发明考虑的卫星j的地球引力加速度a_1,j中包含地球形状的J2摄动力，满足：

$a_{1,j}=\left[\begin{array}{c}{a}_{1x_j}\\ a_{1y_j}\\ a_{1z_j}\end{array}\right]=\left[\begin{array}{c}-\frac{{\mathrm{\μx}}_j}{r_{e\_sat,j}^3}\{1+\frac{3J_2}{2}{\left(\frac{R_e}{r_{e\_sat,j}}\right)}^2\[1-5{\left(\frac{z_j}{r_{e\_sat,j}}\right)}^2\]\}\\ -\frac{{\mathrm{\μy}}_j}{r_{e\_sat,j}^3}\{1+\frac{3J_2}{2}{\left(\frac{R_e}{r_{e\_sat,j}}\right)}^2\[1-5{\left(\frac{z_j}{r_{e\_sat,j}}\right)}^2\]\}\\ -\frac{{\mathrm{\μz}}_j}{r_{e\_sat,j}^3}\{1+\frac{3J_2}{2}{\left(\frac{R_e}{r_{e\_sat,j}}\right)}^2\[3-5{\left(\frac{z_j}{r_{e\_sat,j}}\right)}^2\]\}\end{array}\right]---\left(3\right)$

式中，为待测卫星j的地心距，R_e为地球赤道半径，J₂为二阶带形谐系数，μ为地球质量与引力常数G的乘积，即地球引力常数。

(2)日月引力引起的摄动加速度a_2,j满足：

$a_{2,j}=\left[\begin{array}{c}{a}_{2x_j}\\ a_{2y_j}\\ a_{2z_j}\end{array}\right]=\left[\begin{array}{c}{\mathrm{Gm}}_S\left(\frac{x_j}{r_{e\_Sun}^3}\right)+{\mathrm{Gm}}_m\left(\frac{x_j}{r_{e\_m}^3}\right)\\ {\mathrm{Gm}}_S\left(\frac{y_j}{r_{e\_Sun}^3}\right)+{\mathrm{Gm}}_m\left(\frac{y_j}{r_{e\_m}^3}\right)\\ {\mathrm{Gm}}_S\left(\frac{z_j}{r_{e\_Sun}^3}\right)+{\mathrm{Gm}}_m\left(\frac{z_j}{r_{e\_m}^3}\right)\end{array}\right]---\left(4\right)$

式中，r_{e_Sun}和r_{e_m}分别为地心到太阳和月球的距离，m_m和m_S分别为月球和太阳的质量，G为引力常数。

(3)太阳光压引起的摄动加速度a_3,j满足：

$a_{3,j}=\left[\begin{array}{c}{a}_{3x_j}\\ a_{3y_j}\\ a_{3z_j}\end{array}\right]=\left[\begin{array}{c}{\mathrm{PC}}_r{r}_{e\_Sun}^2\left(\frac{S_j}{m_{sat,j}}\right)\frac{(x_{Sun}-x_j)}{{r_{sat,j\_Sun}}^3}\\ {\mathrm{PC}}_r{r}_{e\_Sun}^2\left(\frac{S_j}{m_{sat,j}}\right)\frac{(y_{Sun}-y_j)}{{r_{sat,j\_Sun}}^3}\\ {\mathrm{PC}}_r{r}_{e\_Sun}^2\left(\frac{S_j}{m_{sat,j}}\right)\frac{(z_{Sun}-z_j)}{{r_{sat,j\_Sun}}^3}\end{array}\right]---\left(5\right)$

式中，P为太阳光压强值；C_r＝1+ε为光压系数，ε＝0.2～0.9为卫星制造中典型材料的反射系数；S_j和m_sat,j分别为卫星j的截面积和质量；[x_Sun,y_Sun,z_Sun]^T为在地心惯性系下的太阳位置矢量，r_{sat,j_Sun}为卫星j到太阳的距离，

设GSO星座由m颗GSO卫星组成，根据上述单颗卫星定轨的状态方程(2)，建立集中式结构的星座定轨状态方程如下：

$\stackrel{\·}{X}=F\[X\left(t\right),t\]+W\left(t\right)=\left[\begin{array}{c}{f}_1\left(X_1\right)\\ .\\ .\\ .\\ f_j\left(X_j\right)\\ .\\ .\\ .\\ f_m\left(X_m\right)\end{array}\right]+\left[\begin{array}{c}{W}_1\left(t\right)\\ .\\ .\\ .\\ W_j\left(t\right)\\ .\\ .\\ .\\ W_m\left(t\right)\end{array}\right]---\left(6\right)$

式中，为状态向量，这里X_j(j＝1,...,m)为待测卫星j对应的状态向量；F为整个星座定轨系统的状态函数矢量；W(t)表示系统噪声，W_j(t)表示待测卫星j的系统噪声。

步骤二：在集中式结构下，建立GSO星座定轨系统的测量方程；

本发明采用地面站测量伪距、星载GNSS接收机测量伪距和伪距率、星间链路测量星间相对距离这三种方式来获取观测信息，据此分别建立三种测量方式对应的测量模型。

地面站利用伪码测距技术获得到所测卫星的伪距测量信息。对于待测卫星j，设地面站i在地心惯性系下的位置矢量为它到待测卫星j的距离为那么同时利用n_j个地面站对待测卫星j进行伪码测距的测量方程如下：

z_gr,j＝h_gr,j+v_gr,j (7)

式中，为地面站的联合观测函数向量，这里是地面站i到待测卫星j的距离；为经过电离层、对流层延迟补偿之后残余的伪距观测噪声，其统计特性满足零均值，标准差为σ_ρ,j的白噪声，σ_ρ,j的大小可以反映伪距测量精度。

设卫星j上的星载GNSS接收机可观测到l_j个GNSS卫星，对应的伪距测量方程为：

z_gnss,j＝h_gnss,j+v_gnss,j (8)

式中，为星载GNSS观测函数向量，这里，δl_tu,j为卫星j的星载GNSS接收机等效钟差，为待测星j到第s颗GNSS卫星的距离 位第s颗GNSS卫星在地心直角惯性系下的位置矢量；为观测噪声向量，为经过电离层、对流层延迟补偿后的GNSS卫星s到卫星j伪距观测的残余噪声，设为零均值，标准差为的白噪声，通过的大小反映星载GNSS接收机伪距测量精度。

为了提高星间相对定位精度以及星座的自主定轨能力，利用建立的星间链路，采用单向伪码测距方法实现星座内两两卫星间的相对距离测量，建立星座的星间链路伪距观测方程如下：

z_sat＝h_sat+v_sat (9)

式中，h_sat＝[r_1,2 … r_j,p … r_m-1,m]^T(j＝1,…,m；p＝j+1,…,m；且j＜p)为星间链路伪距观测函数向量，其中为卫星j与卫星p间的距离，卫星j和卫星p在地心直角惯性系下位置分别为[x_j,y_j,z_j]^T和[x_p,y_p,z_p]^T；v_sat为星间链路伪码测距的观测噪声向量，是零均值，标准差为σ_sat的白噪声，σ_sat的大小反映星间测距的精度。

在集中式结构下，根据地面站伪距测量模型式(7)、星载GNSS接收机伪距测量模型式(8)和星间链路相对距离测量模型式(9)，建立GSO星座的测量方程如下：

式中，h为星座整体的观测函数向量，v为系统测量噪声。m为GSO星座中卫星的个数，n_j为待测卫星j所选用的地面站数目，l_j为待测卫星j可用的GNSS可见星数目。当参与融合的各测量子系统可用测量数目改变时，测量方程的维数随之自适应调整，如观测方程式(10)的维数可随着星座内各卫星星载GNSS接收机的可见星数目实时变化而发生改变。

步骤三：定轨测量组合的优化策略；

在高轨星座整体运行的各轨道周期内，一方面不同地面站分布提供的可用观测信息的定轨性能不同，另一方面，星座对GNSS导航星在不同时间段的可见情况不同，因此，需要实时判断测量信息的可用性和有效性以保证轨道参数估计精度。本发明对地面站和GNSS这两种测量方式进行优选判别，从而确定有效测量组合。

(1)地面站选择策略

站-星的空间几何分布是影响卫星定轨精度的重要因素之一，因此合理的地面站分布对于提高定轨精度非常重要。精度因子DOP是描述地面站-卫星空间几何分布的特征量，它反映了地面站-卫星空间几何分布造成的测量误差与卫星状态估计误差的比例系数。

在将式(9)中与观测待测卫星j对应的n_j个地面站进行观测的伪距观测方程h_gr,j围绕卫星j的一步预测位置估计进行一阶泰勒级数展开进行线性化处理时，可得到地面站伪距测量的观测函数向量h_gr,j对应的Jacobian矩阵：

$H_{gr,j,\mathrm{\ρ}}=\frac{\∂h_{gr,j}}{\∂X_{\mathrm{\ρ}}}{|}_{{\hat{X}}_{\mathrm{\ρ}}\left(k\right|k-1)}=\left[\begin{array}{ccc}\frac{(x_j-x_g^{\left(1\right)})}{r_{g\_sat,j}^{\left(1\right)}}& \frac{(y_j-y_g^{\left(1\right)})}{r_{g\_sat,j}^{\left(1\right)}}& \frac{(z_j-z_g^{\left(1\right)})}{r_{g\_sat,j}^{\left(1\right)}}\\ & .& \\ & .& \\ & .& \\ \frac{(x_j-x_g^{\left(i\right)})}{r_{g\_sat,j}^{\left(i\right)}}& \frac{(y_j-y_g^{\left(i\right)})}{r_{g\_sat,j}^{\left(i\right)}}& \frac{(z_j-z_g^{\left(i\right)})}{r_{g\_sat,j}^{\left(i\right)}}\\ & .& \\ & .& \\ & .& \\ \frac{(x_j-x_g^{\left(n_j\right)})}{r_{g\_sat,j}^{\left(n_j\right)}}& \frac{(y_j-y_g^{\left(n_j\right)})}{r_{g\_sat,j}^{\left(n_j\right)}}& \frac{(z_j-z_g^{\left(n_j\right)})}{r_{g\_sat,j}^{\left(n_j\right)}}\end{array}\right]---\left(11\right)$

式中，待测卫星j在地心惯性系下的位置矢量[x_j,y_j,z_j]^T，地面站i在地心惯性系下的位置矢量为它到待测卫星j的距离为

当仅根据测量方程采用单点定位方法实现待测卫星j的位置解算时，由H_gr,j,ρ描述的地面站与卫星间的相对几何关系决定了定轨精度。位置估计误差ε_X的方差矩阵为：

$Cov\left({\mathrm{\ϵ}}_X\right)=E\left({\mathrm{\ϵ}}_X{\mathrm{\ϵ}}_X^T\right)={({H_{gr,j,\mathrm{\ρ}}}^T\·H_{gr,j,\mathrm{\ρ}})}^{-1}{\mathrm{\σ}}_{URE}^2---\left(12\right)$

式中，为地面站测量误差方差。令矩阵G＝(H_gr,j,ρ ^T·H_gr,j,ρ)^-1，可见G由星站间的相对位置关系确定。

式(12)清晰地表明了位置估计误差协方差矩阵是由测量误差方差通过权系数阵G放大后得到。可见定轨精度与伪距测量误差和地面站的几何分布有关。为此定义位置精度因子如下：

$DOP=\sqrt{tr{({H_{gr,j,\mathrm{\ρ}}}^T\·H_{gr,j,\mathrm{\ρ}})}^{-1}}=\sqrt{tr\left(G\right)}---\left(13\right)$

由于地面站与卫星所围成多面体几何特性可以直接反应DOP值的大小，通过求解星-站构成的多面体体积即可获得地面站与卫星所围成多面体的几何特性，因此采用最大单位矢端多面体体积法来确定地面站的最佳分布。最大单位矢端多面体体积法首先求得卫星j到n个地面站的单位方向矢量然后计算由单位方向矢量(e₁～e_n)所构成的矢端多面体的体积，从所有组合中选取体积最大的一组作为地面站最佳分布方案。单位矢端多面体体积示意图如图3所示，图中卫星j到地面站1、2、i、n的单位方向矢量分别为e₁、e₂、e_i、e_n，由多个单位方向矢量张成的多面体体积反映了地面站与卫星间的几何关系。

鉴于GSO星座的“静地”、高轨和星座构型相对稳定等特点，站-星间几何关系基本保持不变。因此，可以利用该方法根据注入候选地面站信息在定轨解算之前事先确定最优地面站分布用于导航解算。若需更改地面站分布方案，则应重新注入候选地面站信息，并采用最大单位矢端多面体体积法重新选择地面站。

(2)GNSS选星策略

在可用地面站分布确定基础上，确定GNSS选星策略法。与地面站分布可事前确定不同，星座中每个GSO卫星可见的GNSS导航卫星可能不同，且几何分布和数目均是实时变化的。当可见星数目多时，若都参与导航解算，会增加计算量，而其中几何分布不好的导航星还会影响定轨精度。因此，需要实时进行GNSS选星。

GNSS选星策略有以下两种：

(1)GNSS独立选星

当GNSS可见星数≤4时，利用全部GNSS可见星的测量信息；当可见星数>4时，从候选GNSS可见星中优选出4颗。由于GSO卫星在轨运行期间星载GNSS接收机可见星数>4的时间区间在整个轨道周期中所占百分比较低，因此在大多数时间所有GNSS可见星测量均可使用。而在可见星数>4的时间段内，与低轨载体或地面应用相比，候选GNSS可见星数相对较少，因此选择简单实用的基于精度因子的快速选星方法，可满足算法实时性需求。

(2)基于地面站的GNSS选星

一般地，实际系统中星载GNSS接收机的伪距测量精度比地面站的伪距测量精度低。当星载GNSS接收机测量精度很低时，即便选择了较好构型的GNSS可见星，也不能保证在地面站测量基础上引入星载GNSS测量可以提高定轨精度。因此，需要利用地面站的测量信息辅助GNSS选星。

这里综合考虑地面站和GNSS可见星与待测卫星的几何分布和测量精度，采用加权精度因子法(Weighted DOP，WDOP)进行GNSS选星。

定义加权精度因子：

$WDOP=\sqrt{tr{\left(H_{g\_gnss,j}^T{\mathrm{Weight}}_{g\_gnss,j}^{-1}{H}_{g\_gnss,j}\right)}^{-1}}---\left(14\right)$

式中，Η_{g_gnss,j}为将卫星j的地面站和星载GNSS接收机的测量函数向量在一步预测位置估计处进行线性化后所得到的Jacobian矩阵，即这里，H_gr,j,ρ为地面站伪距观测函数向量h_gr,j的Jacobian矩阵，参见式(11)。H_gnss,j,ρ为GNSS伪距观测函数向量h_gnss,j的Jacobian矩阵，即

$H_{gnss,j,\mathrm{\ρ}}=\frac{\∂h_{gnss,j}}{\∂X_{\mathrm{\ρ}}}{|}_{{\hat{X}}_{\mathrm{\ρ}}\left(k\right|k-1)}=\left[\begin{array}{ccc}\frac{(x_j-x_{gnss}^{\left(1\right)})}{r_{gns\_sat,j}^{\left(1\right)}}& \frac{(y_j-y_{gnss}^{\left(1\right)})}{r_{gnss\_sat,j}^{\left(1\right)}}& \frac{(z_j-z_{gnss}^{\left(1\right)})}{r_{gnss\_sat,j}^{\left(1\right)}}\\ & .& \\ & .& \\ & .& \\ \frac{(x_j-x_{gnss}^{\left(s\right)})}{r_{gns\_sat,j}^{\left(s\right)}}& \frac{(y_j-y_{gnss}^{\left(s\right)})}{r_{gnss\_sat,j}^{\left(s\right)}}& \frac{(z_j-z_{gnss}^{\left(s\right)})}{r_{gnss\_sat,j}^{\left(s\right)}}\\ & .& \\ & .& \\ & .& \\ \frac{(x_j-x_{gnss}^{\left(l_j\right)})}{r_{gns\_sat,j}^{\left(l_j\right)}}& \frac{(y_j-y_{gnss}^{\left(l_j\right)})}{r_{gnss\_sat,j}^{\left(l_j\right)}}& \frac{(z_j-z_{gnss}^{\left(l_j\right)})}{r_{gnss\_sat,j}^{\left(l_j\right)}}\end{array}\right]---\left(15\right)$

式中，待测卫星j在地心惯性系下的位置矢量[x_j,y_j,z_j]^T；第s颗可见GNSS导航卫星在地心惯性系下的位置矢量为它到待测卫星j的距离为

为卫星j的地面站和星载GNSS接收机伪距观测误差的加权阵，这里weight_g,j和weight_gnss,j分别为卫星j与地面站和星载GNSS接收机测量对应的对角加权阵，满足：和设各地面站测量精度相当，测量噪声标准差为σ_g，那么和其中i＝1～n_j，s＝1～l_j。

在单独地面站、不同数量候选GNSS可见星的不同组合+地面站确定的WDOP值中，选取使WDOP值最小的GNSS卫星作为最佳GNSS选星结果。

方法(1)中只是简单地针对GNSS可见星≥4的情况，考虑可见GNSS卫星与待测卫星的几何构型，利用最佳精度因子法进行选星，在有无地面站测量的情况下均可以实施；而方法(2)中，同时考虑地面站和可见GNSS卫星与待测卫星的几何构型和测量精度，采用加权精度因子法进行选星，在存在地面站测量的情况下实施。

综上所述，本发明算法在进行GNSS选星时，无地面站测量的情况下采用方法(1)，在有地面站测量的情况下采用方法(2)。

步骤四：根据上述建立的基于地面站/星间链路/GNSS的GSO星座联合定轨系统模型，确定实现轨道参数估计的滤波方法；

该非线性系统的状态估计需要采用非线性滤波方法，如EKF(Extended KalmanFilter)、UKF(Unscented Kalman Filter)、粒子滤波和容积滤波等方法。本发明以工程实用性强的EKF法为例，实现卫星轨道参数估计，其他方法如UKF、粒子滤波和容积滤波等也可用于该系统的状态估计，进一步提高估计精度。

EKF的基本思想是利用泰勒级数展开的方法对非线性的状态方程和观测方程进行线性化处理，为此需要计算其Jacobian矩阵。然后采用经典卡尔曼滤波进行估计。

在利用EKF估计系统状态时，对由式(6)和式(10)所描述的连续系统进行离散化可得：

$\left\{\begin{array}{c}X\left(k\right)=X(k-1)+F\[X(k-1)\]\mathrm{\Δ}t+A\[X(k-1)\]F\[X(k-1)\]{\mathrm{\Δt}}^2/2+W(k-1)\\ Z\left(k\right)=h\[X\left(k\right)\]+v\left(k\right)\end{array}\right.---\left(16\right)$

其中，是状态向量；F是系统状态函数；W(k-1)为系统噪声，是零均值、方差为Q(k-1)的高斯白噪声；Δt为离散时间，为Jacobian矩阵。Z(k)是系统观测向量；h[X(k)]是系统观测函数向量，参见式(11)；v(k)为测量噪声，是零均值、方差为R(k)的高斯白噪声向量。

对经过离散化后的非线性系统采用离散EKF方法进行状态估计。具体的滤波步骤如下：

(1)时间更新

计算一步预测状态：

$\hat{X}\left(k\right|k-1)=\hat{X}(k-1)+F\[\hat{X}(k-1)\]\·\mathrm{\Δ}t+A\[\hat{X}(k-1)\]\·F\[\hat{X}(k-1)\]\frac{{\mathrm{\Δt}}^2}{2}---\left(17\right)$

式中，为(k-1)时刻的状态估计，为k时刻的一步预测状态估计。

计算状态转移矩阵：

$\mathrm{\Φ}\left(k\right|k-1)=I+A\[\hat{X}(k-1)\]\·\mathrm{\Δ}t---\left(18\right)$

式中，Φ(k|k-1)为从(k-1)时刻到k时刻的状态状态转移矩阵。

计算一步预测误差协方差矩阵：

P(k|k-1)＝Φ(k|k-1)P(k-1)Φ^T(k|k-1)+Q(k-1) (19)

式中，P(k|k-1)为一步预测状态的误差协方差矩阵，P(k-1)为(k-1)时刻估计状态的误差协方差矩阵，Q(k-1)为(k-1)时刻系统噪声方差阵。

(2)测量更新：

计算增益矩阵：

K(k)＝P(k|k-1)H^T(k)[H(k)P(k|k-1)H^T(k)+R(k)]^-1 (20)

式中，K(k)为k时刻的增益矩阵；为基于一步预测状态估计计算的观测向量h的Jacobian矩阵；R(k)为k时刻测量噪声协方差阵。

计算状态估计：

$\hat{X}\left(k\right)=\hat{X}\left(k\right|k-1)+K\left(k\right)\{Z\left(k\right)-h\[\hat{X}\left(k\right|k-1)\]\}---\left(21\right)$

式中，为k时刻的状态估计。

计算状态估计误差协方差阵：

P(k)＝[I-K(k)H(k)]P(k|k-1)[I-K(k)H(k)]^T+K(k)R(k)K^T(k) (22)

式中，P(k)为k时刻的状态估计误差的协方差阵。

步骤五：在集中式结构下，基于地面站/星间链路/GNSS的GSO星座联合定轨算法的具体实现；

根据本发明提出的定轨测量组合的优化策略自适应确定有效观测，利用EKF实现基于地面站/星间链路/GNSS的GSO星座联合定轨算法的流程图参见图4，具体过程如下：

(1)数据初始化，初始化的参数包括：k＝1，初始状态及误差协方差阵P(0)和系统噪声方差阵Q、初始地面站分布方案及参数。

(2)时间更新。利用(k-1)时刻的状态估计及误差协方差阵P(k-1)作为k时刻的状态初值及误差协方差阵，根据式(17)和式(19)分别计算状态一步预测值及预测误差协方差矩阵P(k|k-1),实现EKF估计方法的时间更新；

(3)初始化卫星编号j＝1；

(4)获取卫星j星载GNSS接收机观测数据，确定候选GNSS可见星数；

(5)获取卫星j的地面站观测数据，检测是否存在地面站测量，如果不存在则按顺序执行步骤(6)；否则按顺序执行步骤(7)；

(6)如果不存在地面站测量，则直接根据卫星j的候选GNSS可见星个数判断。当候选GNSS可见星数≥4时，从中确定4颗GNSS可见星的所有组合及相应的精度因子，选取最小精度因子所对应的可见星组合作为卫星j的GNSS的优选结果；当可见星数<4颗时，所有GNSS可见卫星均作为卫星j的GNSS的优选结果；跳到(10)；

(7)如果存在地面站测量，根据卫星j的地面站观测数据判断是否需要更新卫星j的地面站信息，如果需要，则注入新的可用地面站信息，利用最大单位矢端多面体体积法，在线确定卫星j的最优地面站分布；否则沿用k-1时刻卫星j选用的地面站分布。

(8)根据卫星j的地面站分布，利用卫星j获得的地面站有效测量信息，构建卫星j的地面站相关的测量模型和测量噪声方差阵。

(9)基于步骤(6)卫星j的地面站分布，根据卫星j的候选GNSS卫星的不同选星个数，计算卫星j在地面站和不同星载GNSS测量组合下的加权精度因子，从中选取最小加权精度因子所对应的可见星组合作为卫星j的GNSS优选结果；

(10)根据卫星j的GNSS可见星优选结果，利用所获得GNSS有效测量信息，构建卫星j的星载接收机相关的测量模型和测量噪声方差阵；

(11)判断是否星座内所有卫星的有效观测组合优选已结束。如果j＝m，继续下一步；否则j＝j+1，并返回到(4)；

(12)获取星间链路星间观测信息，构建测量模型和测量噪声方差阵中与星间链路测量相关的部分。

(13)利用所有可用测量信息对应的测量模型及测量噪声方差阵，组合构成k时刻的合成测量模型和测量噪声方差阵R(k)。

(14)测量更新。根据式(20)～(22)依次计算k时刻的增益矩阵K(k)、状态估计值及误差协方差矩阵P(k)，实现EKF的测量更新。

(15)判断滤波过程是否结束。如需继续进行滤波，k＝k+1，并返回到(2)。

本发明首先研究基于地面站的卫星定轨局限性，为了提高星间相对测量精度、星座自主定轨性能，降低对地面站的依赖，引入了星间链路测量和星载GNSS测量；之后确定定轨测量方式的组合策略，最后在集中式结构下，采用EKF估计方法，实现基于地面站/星间链路/GNSS的GSO星座联合定轨。

本发明提出的一种基于地面站/星间链路/GNSS的地球同步轨道卫星星座定轨方法，该方法首先在集中式结构下建立GSO星座的状态方程；然后在建立基于最大单位矢端多面体体积法的地面站测量信息优选判决策略基础上，确定需要融合的地面站测量信息；在建立基于加权精度因子法的地面站+GNSS组合测量信息优选策略基础上，对需要实时融合的星载GNSS接收机测量信息进行优选判别；根据地面站和星载GNSS接收机的优选结果，建立基于地面站/星间链路/GNSS的测量方程；最后采用EKF方法实现GSO星座联合定轨。

实施例：

以某GSO星座为例，该星座由3个IGSO卫星和1个GEO卫星构成。卫星编号1～3为IGSO，采用小轨道偏心率(接近圆轨道)和小轨道倾角(星下点轨迹区域分布小)的地球同步卫星；卫星4为GEO卫星。在仿真中发现，星座中各卫星的定轨和相对定位精度近似，因此文中仅以编号为1的IGSO为例进行说明。

仿真条件：

采用STK软件中的高精度轨道预报器(HPOP)来产生星座运行过程中的真实轨道数据，以模拟星座实际运行情况。

(1)STK仿真参数设置

仿真时间为2013年1月1日00:00:00—2013年1月2日00:00:00；历时一个轨道周期(24h)，采样间隔T＝1s。地球模型采用JGM-3(Joint Gravity Model)，考虑21×21阶非球形摄动；其他摄动项包括太阳引力、月球引力、太阳光压(光压系数C_r＝1.3)等，不考虑大气阻力；卫星截面积与卫星质量之比

候选地面站选取我国境内现有的9个卫星地面监控站，分别为北京、三亚、拉萨、福州、哈尔滨、呼和浩特、昆明、武汉和乌鲁木齐。

(2)基本滤波参数设置

仿真时长为GSO卫星的1个轨道周期(24h)，滤波周期T＝1s。系统状态方程采用考虑带J2项、日月引力和太阳光压等摄动项的轨道动力学模型。

在集中式结构下，基于地面站/星间链路/GNSS的GSO四星星座联合定轨的系统噪声方差阵为Q＝diag([Q₁(0) Q₂(0) Q₃(0) Q₄(0)])，式中q_1j＝0.01m,q_2j＝0.06m/s,q_3j＝0.001m,q_4j＝0.01m/s；初始状态估计误差协方差阵为：P(0)＝diag([P₁(0) P₂(0) P₃(0) P₄(0)])式中，σ_x,j＝σ_y,j＝σ_z,j＝10m,

(3)评估方法

在得到各个仿真时刻的估计误差后，采用统计方法(计算误差数据的均值和标准差)来评估星座中卫星的绝对定轨和相对定位的效果。统计样本采用滤波估计收敛后的状态估计误差。设mean为均值，std为标准差，具体计算方法如下：

$mean=\frac{1}{N}\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}Error\left(k\right)---\left(23\right)$

$std=\sqrt{\frac{1}{N}\underset{k=1}{\overset{N}{\mathrm{\&Sigma;}}}{\left(Error\right(k)-mean)}^2}---\left(24\right)$

式中Error(k)表示统计样本中第k时刻的估计误差，N表示统计样本个数。对评估中涉及的各种误差定义如下：

卫星轨道参数估计误差：利用k时刻的卫星轨道参数估计值和卫星轨道参数的标称值X(k)计算k时刻卫星的轨道参数估计误差dX(k)：

$dX\left(k\right)=X\left(k\right)-\hat{X}\left(k\right)---\left(25\right)$

式中，

X(k)＝[x(k),y(k),z(k),v_x(k),v_y(k),v_z(k)]^T

dX(k)＝[dx(k),dy(k),dz(k),dv_x(k),dv_y(k),dv_z(k)]^T

k时刻的卫星距离估计误差dr(k)为：

$dr\left(k\right)=\sqrt{dx{\left(k\right)}^2+dy{\left(k\right)}^2+dz{\left(k\right)}^2}---\left(26\right)$

k时刻卫星速率估计误差dv(k)为：

$dv\left(k\right)=\sqrt{v_x{\left(k\right)}^2+v_y{\left(k\right)}^2+v_z{\left(k\right)}^2}-\sqrt{{\hat{v}}_x{\left(k\right)}^2+{\hat{v}}_y{\left(k\right)}^2+{\hat{v}}_z{\left(k\right)}^2}---\left(27\right)$

实施方式：集中式GSO星座联合定轨系统模型的状态向量由四星的位置、速度及星载接收机钟差和频差构成，采用地面站、星间链路和星载GNSS三种测量方式，重点对比定轨测量方式的组合策略对定轨精度的影响。

步骤一：建立集中式GSO星座定轨系统的状态方程

在地心直角惯性坐标系下，卫星j扩维状态向量为其中[x_j,y_j,z_j]^T和分别为该星的位置矢量和速度矢量，为该星载GNSS接收机的等效钟差和频差。结合轨道动力学模型式(3)～(5)可建立卫星j定轨的状态方程：

${\stackrel{\&CenterDot;}{X}}_j\left(t\right)=f_j\&lsqb;X_j\left(t\right),t\&rsqb;+W_j\left(t\right)\&DoubleRightArrow;\left[\begin{array}{c}{\stackrel{\&CenterDot;}{x}}_j\\ {\stackrel{\&CenterDot;}{y}}_j\\ {\stackrel{\&CenterDot;}{z}}_j\\ {\stackrel{\&CenterDot;}{v}}_{x_j}\\ {\stackrel{\&CenterDot;}{v}}_{y_j}\\ {\stackrel{\&CenterDot;}{v}}_{z_j}\\ \mathrm{\&delta;}{\stackrel{\&CenterDot;}{l}}_{tu,j}\\ \mathrm{\&delta;}{\stackrel{\&CenterDot;}{l}}_{tru,j}\end{array}\right]=\left[\begin{array}{c}{v}_{x_j}\\ v_{y_j}\\ v_{z_j}\\ a_{1x_j}+a_{2x_j}+a_{3x_j}\\ a_{1y_j}+a_{2y_j}+a_{3y_j}\\ a_{1z_j}+a_{2z_j}+a_{3z_j}\\ {\mathrm{\&delta;l}}_{tru,j}\\ -\frac{1}{{\mathrm{\&tau;}}_{tru}}{\mathrm{\&delta;l}}_{tru,j}\end{array}\right]+\left[\begin{array}{c}{w}_{x_j}\\ w_{y_j}\\ w_{z_j}\\ w_{v_{x_j}}\\ w_{v_{y_j}}\\ w_{v_{z_j}}\\ w_{tu,j}\\ w_{tru,j}\end{array}\right]---\left(28\right)$

式中，卫星j的地球引力(包括地球形状的J2摄动力)、日月引力和太阳光压引起的摄动加速度分别为a_1,j、a_2,j和a_3，j，具体定义参见式(3)～式(5)。W_j(t)表示系统噪声，包含未建模摄动，满足E[W_j(t)]＝0，E[W_j(t)W^T _j(τ)]＝Q_j(t)δ(t-τ)，δ(t-τ)是狄拉克函数，即

根据上述单颗卫星的状态方程，可以建立集中式结构的四星星座状态方程如下：

$\stackrel{\&CenterDot;}{X}=F\&lsqb;X\left(t\right),t\&rsqb;+W\left(t\right)=\left[\begin{array}{c}{f}_1\left(X_1\right)\\ f_2\left(X_2\right)\\ f_3\left(X_3\right)\\ f_4\left(X_4\right)\end{array}\right]+\left[\begin{array}{c}{W}_1\left(t\right)\\ W_2\left(t\right)\\ W_3\left(t\right)\\ W_4\left(t\right)\end{array}\right]---\left(29\right)$

式中，为包含星座中四颗卫星所有状态的状态向量，W(t)表示系统噪声。

步骤二：建立集中式GSO星座定轨系统的测量方程

在集中式结构下，四星星座采用地面站、星载GNSS接收机和星间链路测量的测量方程如下：

式中，h_gr,j为卫星j的地面站伪距测量向量，其中表示由地面站i测量待测卫星j的伪距，为经过电离层、对流层延迟补偿之后残余的伪距观测噪声向量，为零均值、标准差为σ_ρ,j的白噪声，参见式(7)。h_gnss,j为卫星j星载GNSS接收机测得的l_i个GNSS可见星的伪距观测函数，其中，为卫星j的星载GNSS接收机测得第s个GNSS可见星的伪距测量函数，为观测噪声向量，为经过电离层、对流层延迟补偿后第s颗GNSS卫星到卫星j的残余伪距观测噪声，设为零均值，标准差为的白噪声，参见式(8)。h_sat＝[r_1,2 r_1,3 r_1,4 r_2,3 r_2,4 r_3,4]^T为四星星座星间链路测量的星间距观测函数，r_j,p(j＝1,…,4；p＝j+1,…,4；且j＜p)为卫星j与卫星p间的距离函数，v_sat为星间链路伪码测距噪声，是零均值、标准差为σ_sat的白噪声，参见式(9)。

步骤三：在集中式结构下，采用EKF估计方法实现基于地面站/星间链路/GNSS的GSO四星星座联合定轨，并仿真分析基于地面站/星间链路/GNSS的GSO四星星座集中式联合定轨性能。

结合GSO四星星座集中式联合定轨流程图(图4)和发明内容步骤五，采用EKF估计方法实现基于地面站/星间链路/GNSS的GSO四星星座联合定轨。

(1)数据初始化，初始化的参数包括：k＝1，m＝4,n₁＝n₂＝n₃＝n₄＝4,初始状态及误差协方差阵P(0)和系统噪声方差阵Q、初始地面站分布参数。

(2)时间更新。利用k-1时刻的状态估计及误差协方差阵P(k-1)作为k时刻的状态初值及误差协方差阵，根据式(17)和式(19)分别计算状态一步预测值及预测误差协方差矩阵P(k|k-1),实现EKF估计方法的时间更新；

(3)初始化卫星编号j＝1；

(4)获取卫星j星载GNSS接收机观测数据，确定候选GNSS可见星数；

(5)获取卫星j的地面站观测数据，检测是否存在地面站测量，如果不存在则按顺序执行步骤(6)；否则按顺序执行步骤(7)；

(6)如果不存在地面站测量，则直接根据卫星j的候选GNSS可见星个数判断。当候选GNSS可见星数≥4时，从中确定4颗GNSS可见星的所有组合及相应的精度因子，选取最小精度因子所对应的可见星组合作为卫星j的GNSS的优选结果；当可见星数<4颗时，所有GNSS可见卫星均作为卫星j的GNSS的优选结果；跳到(10)；

(7)如果存在地面站测量，根据卫星j的地面站观测数据判断是否需要更新卫星j的地面站信息，如果需要，则注入新的可用地面站信息，利用最大单位矢端多面体体积法，在线确定卫星j的最优地面站分布；否则沿用k-1时刻卫星j选用的地面站分布。

(8)根据卫星j的地面站分布，利用卫星j获得的地面站有效测量信息，构建卫星j的地面站相关的测量模型和测量噪声方差阵。

(9)基于步骤(6)卫星j的地面站分布，根据卫星j的候选GNSS卫星的不同选星个数，计算卫星j在地面站和不同星载GNSS测量组合下的加权精度因子，从中选取最小加权精度因子所对应的可见星组合作为卫星j的GNSS优选结果；

(10)根据卫星j的GNSS可见星优选结果，利用所获得GNSS有效测量信息，构建卫星j的星载接收机相关的测量模型和测量噪声方差阵；

(11)判断是否星座内所有卫星的有效观测组合优选已结束。如果j＝m，继续下一步；否则j＝j+1，并返回到(4)；

(12)获取星间链路星间观测信息，构建测量模型和测量噪声方差阵中与星间链路测量相关的部分。

(13)利用所有可用测量信息对应的测量模型及测量噪声方差阵，组合构成k时刻的合成测量模型和测量噪声方差阵R(k)。

(14)测量更新。根据式(20)～(22)依次计算k时刻的增益矩阵K(k)、状态估计值及误差协方差矩阵P(k)，实现EKF的测量更新。

(15)判断滤波过程是否结束。如需继续进行滤波，k＝k+1，并返回到(2)。

本发明在分析基于地面站的卫星定轨的局限性基础上，从提高星座内卫星的绝对定位和相对定位精度、星座自主运行能力、降低对地面站的依赖性出发，融入了星间链路和星载GNSS接收机的测量信息，制定了地面站和GNSS测量信息的组合策略；最后，在集中式结构下，采用以EKF为例的非线性滤波方法实现了基于地面站/星间链路/GNSS的GSO星座联合定轨。

在地面站测量误差均值为零，标准差为0.1m；星间测距误差均值为零，标准差为0.01m；GNSS测量误差均值为零，标准差为10m的情况下，有以下2种仿真结果：

(1)采用最佳精度因子法对可见GNSS卫星进行选星(参见步骤四的方法(1))。由于星座中四颗卫星的定轨精度类似，因此以卫星1为例，其定轨距离误差均值为0.3442m，标准差为0.1893m，定速精度均值为1.0471×10^-5m/s，标准差为0.0092m/s。

(2)利用加权精度因子法对可见GNSS卫星进行选星(参见步骤四的方法(2))，卫星1的定轨距离误差均值为0.3397m，标准差为0.1883m，定速精度均值为1.9756×10^-6m/s，标准差为0.0090m/s。

对比上述两种情况的仿真结果可知，当与地面站测量相比，星载GNSS接收机的伪距测量精度差别较大时，结合地面站测量，采用加权精度因子法对GNSS可见卫星进行选星的定轨精度有明显提高。验证了本发明方法的有效性。

Claims (1)

1.一种基于地面站/星间链路/GNSS联合测量的地球同步轨道星座定轨方法，具体包括以下几个步骤：

步骤一：在集中式结构下，建立GSO星座定轨系统的状态方程；

基于地球质心引力以及J2项、日月引力和太阳光压三种摄动力构成的轨道动力学模型建立系统状态方程；

根据钟差和频差与距离之间的关系，将其转化为等效测距误差δl_tu和距离变化率误差δl_tru，建立星载GNSS接收机等效时钟误差的状态方程如下：

$\begin{array}{c}\mathrm{\&delta;}{\stackrel{\&CenterDot;}{l}}_{tu}={\mathrm{\&delta;l}}_{tru}+w_{tu}\\ \mathrm{\&delta;}{\stackrel{\&CenterDot;}{l}}_{tru}=-\frac{1}{{\mathrm{\&tau;}}_{tru}}{\mathrm{\&delta;l}}_{tru}+w_{tru}\end{array}---\left(1\right)$

式中，τ_tru为相关时间常数，w_tu和w_tru为零均值高斯白噪声；

在地心直角惯性坐标系下，结合轨道动力学模型和接收机等效钟差模型，建立对星座中待测卫星j进行定轨的状态方程如下：

${\stackrel{\&CenterDot;}{X}}_j\left(t\right)=f_j\&lsqb;X_j\left(t\right),t\&rsqb;+W_j\left(t\right)---\left(2\right)$

式中，待测卫星j对应的状态向量为[x_j,y_j,z_j]^T和分别为该星的位置矢量和速度矢量，为该星载GNSS接收机的等效钟差和频差；f_j为系统状态函数，W_j(t)为系统噪声，满足均值为零，方差为Q_j(t)的高斯白噪声；式(2)进一步详细写为：

$\left[\begin{array}{c}{\stackrel{\&CenterDot;}{x}}_j\\ {\stackrel{\&CenterDot;}{y}}_j\\ {\stackrel{\&CenterDot;}{z}}_j\\ {\stackrel{\&CenterDot;}{v}}_{x_j}\\ {\stackrel{\&CenterDot;}{v}}_{y_j}\\ {\stackrel{\&CenterDot;}{v}}_{z_j}\\ \mathrm{\&delta;}{\stackrel{\&CenterDot;}{l}}_{tu,j}\\ \mathrm{\&delta;}{\stackrel{\&CenterDot;}{l}}_{tru,j}\end{array}\right]=\left[\begin{array}{c}{v}_{x_j}\\ v_{y_j}\\ v_{z_j}\\ a_{1x_j}+a_{2x_j}+a_{3x_j}\\ a_{1y_j}+a_{2y_j}+a_{3y_j}\\ a_{1z_j}+a_{2z_j}+a_{3z_j}\\ {\mathrm{\&delta;l}}_{tru,j}\\ -\frac{1}{{\mathrm{\&tau;}}_{tru}}{\mathrm{\&delta;l}}_{tru,j}\end{array}\right]+\left[\begin{array}{c}{w}_{x_j}\\ w_{y_j}\\ w_{z_j}\\ w_{v_{x_j}}\\ w_{v_{y_j}}\\ w_{v_{z_j}}\\ w_{tu,j}\\ w_{tru,j}\end{array}\right]$

其中，和为速度噪声分量，和是未建模的摄动加速度分量，w_tu,j和w_tru,j分别为等效钟差和频差噪声；设待测卫星j受到的地球引力、日月引力和太阳光压引起的加速度分别为a_1,j、a_2,j和a_3，j，具体表达式如下：

(1)卫星j的地球引力加速度a_1,j中包含地球形状的J2摄动力，满足：

$a_{1,j}=\left[\begin{array}{c}{a}_{1x_j}\\ a_{1y_j}\\ a_{1z_j}\end{array}\right]=\left[\begin{array}{c}-\frac{{\mathrm{\&mu;x}}_j}{r_{e\_sat,j}^3}\left\{1+\frac{3J_2}{2}{\left(\frac{R_e}{r_{e\_sat,j}}\right)}^2\&lsqb;1-5{\left(\frac{z_j}{r_{e\_sat,j}}\right)}^2\&rsqb;\right\}\\ -\frac{{\mathrm{\&mu;y}}_j}{r_{e\_sat,j}^3}\left\{1+\frac{3J_2}{2}{\left(\frac{R_e}{r_{e\_sat,j}}\right)}^2\&lsqb;1-5{\left(\frac{z_j}{r_{e\_sat,j}}\right)}^2\&rsqb;\right\}\\ -\frac{{\mathrm{\&mu;z}}_j}{r_{e\_sat,j}^3}\left\{1+\frac{3J_2}{2}{\left(\frac{R_e}{r_{e\_sat,j}}\right)}^2\&lsqb;1-5{\left(\frac{z_j}{r_{e\_sat,j}}\right)}^2\&rsqb;\right\}\end{array}\right]---\left(3\right)$

式中，为待测卫星j的地心距，R_e为地球赤道半径，J₂为二阶带形谐系数，μ为地球质量与引力常数G的乘积，即地球引力常数；

(2)日月引力引起的摄动加速度a_2,j满足：

$a_{2,j}=\left[\begin{array}{c}{a}_{2x_j}\\ a_{2y_j}\\ a_{2z_j}\end{array}\right]=\left[\begin{array}{c}{\mathrm{Gm}}_S\left(\frac{x_j}{r_{e\_Sun}^3}\right)+{\mathrm{Gm}}_m\left(\frac{x_j}{r_{e\_m}^3}\right)\\ {\mathrm{Gm}}_S\left(\frac{y_j}{r_{e\_Sun}^3}\right)+{\mathrm{Gm}}_m\left(\frac{y_j}{r_{e\_m}^3}\right)\\ {\mathrm{Gm}}_S\left(\frac{z_j}{r_{e\_Sun}^3}\right)+{\mathrm{Gm}}_m\left(\frac{z_j}{r_{e\_m}^3}\right)\end{array}\right]---\left(4\right)$

式中，r_{e_Sun}和r_{e_m}分别为地心到太阳和月球的距离，m_m和m_S分别为月球和太阳的质量，G为引力常数；

(3)太阳光压引起的摄动加速度a_3,j满足：

$a_{3,j}=\left[\begin{array}{c}{a}_{3x_j}\\ a_{3y_j}\\ a_{3z_j}\end{array}\right]=\left[\begin{array}{c}{\mathrm{PC}}_r{r}_{e\_Sun}^2\left(\frac{S_j}{m_{sat,j}}\right)\frac{(x_{Sun}-x_j)}{{r_{sat,j\_Sun}}^3}\\ {\mathrm{PC}}_r{r}_{e\_Sun}^2\left(\frac{S_j}{m_{sat,j}}\right)\frac{(y_{Sun}-y_j)}{{r_{sat,j\_Sun}}^3}\\ {\mathrm{PC}}_r{r}_{e\_Sun}^2\left(\frac{S_j}{m_{sat,j}}\right)\frac{(z_{Sun}-z_j)}{{r_{sat,j\_Sun}}^3}\end{array}\right]---\left(5\right)$

式中，P为太阳光压强值；C_r＝1+ε为光压系数，ε＝0.2～0.9为卫星制造中典型材料的反射系数；S_j和m_sat,j分别为卫星j的截面积和质量；[x_Sun,y_Sun,z_Sun]^T为在地心惯性系下的太阳位置矢量，r_{sat,j_Sun}为卫星j到太阳的距离，

设GSO星座由m颗GSO卫星组成，根据上述单颗卫星定轨的状态方程(2)，建立集中式结构的星座定轨状态方程如下：

$\stackrel{\&CenterDot;}{X}=F\&lsqb;X\left(t\right),t\&rsqb;+W\left(t\right)=\left[\begin{array}{c}{f}_1\left(X_1\right)\\ .\\ .\\ .\\ f_j\left(X_j\right)\\ .\\ .\\ .\\ f_m\left(X_m\right)\end{array}\right]+\left[\begin{array}{c}{W}_1\left(t\right)\\ .\\ .\\ .\\ W_j\left(t\right)\\ .\\ .\\ .\\ W_m\left(t\right)\end{array}\right]---\left(6\right)$

式中，为状态向量，这里X_j为待测卫星j对应的状态向量，j＝1,...,m；F为整个星座定轨系统的状态函数矢量；W(t)表示系统噪声，W_j(t)表示待测卫星j的系统噪声；

步骤二：在集中式结构下，建立GSO星座定轨系统的测量方程；

本发明采用地面站测量伪距、星载GNSS接收机测量伪距和伪距率、星间链路测量星间相对距离这三种方式来获取观测信息，据此分别建立三种测量方式对应的测量模型；

利用地面站利用伪码测距技术获得到所测卫星的伪距测量信息为：

对于待测卫星j，设地面站i在地心惯性系下的位置矢量为它到待测卫星j的距离为同时利用n_j个地面站对待测卫星j进行伪码测距的测量方程如下：

z_gr,j＝h_gr,j+v_gr,j (7)

式中，为地面站的联合观测函数向量，是地面站i到待测卫星j的距离；为经过电离层、对流层延迟补偿之后残余的伪距观测噪声，其统计特性满足零均值，标准差为σ_ρ,j的白噪声；

利用星载GNSS接收机测量伪距为：

设卫星j上的星载GNSS接收机可观测到l_j个GNSS卫星，对应的伪距测量方程为：

z_gnss,j＝h_gnss,j+v_gnss,j (8)

式中，为星载GNSS观测函数向量，s＝1,…,l_j，δl_tu,j为卫星j的星载GNSS接收机等效钟差，为待测星j到第s颗GNSS卫星的距离 为第s颗GNSS卫星在地心直角惯性系下的位置矢量；为观测噪声向量，为经过电离层、对流层延迟补偿后的GNSS卫星s到卫星j伪距观测的残余噪声，设为零均值，标准差为的白噪声；

利用建立的星间链路，采用单向伪码测距方法实现星座内两两卫星间的相对距离测量，建立星座的星间链路伪距观测方程如下：

z_sat＝h_sat+v_sat (9)

式中，h_sat＝[r_1,2 … r_j,p … r_m-1,m]^T为星间链路伪距观测函数向量，j＝1,…,m；p＝j+1,…,m；且j＜p，其中为卫星j与卫星p间的距离，卫星j和卫星p在地心直角惯性系下位置分别为[x_j,y_j,z_j]^T和[x_p,y_p,z_p]^T；v_sat为星间链路伪码测距的观测噪声向量，是零均值，标准差为σ_sat的白噪声；

在集中式结构下，根据地面站伪距测量模型式(7)、星载GNSS接收机伪距测量模型式(8)和星间链路相对距离测量模型式(9)，建立GSO星座的测量方程如下：

式中，h为星座整体的观测函数向量，v为系统测量噪声；m为GSO星座中卫星的个数，n_j为待测卫星j所选用的地面站数目，l_j为待测卫星j可用的GNSS可见星数目；当参与融合的各测量子系统可用测量数目改变时，测量方程的维数随之自适应调整；

步骤三：定轨测量组合的优化策略；

对地面站和GNSS两种测量方式进行优选判别，从而确定有效测量组合；

第一，地面站选择策略

将式(9)中与观测待测卫星j对应的n_j个地面站进行观测的伪距观测方程h_gr,j围绕卫星j的一步预测位置估计进行一阶泰勒级数展开进行线性化处理时，得到地面站伪距测量的观测函数向量h_gr,j对应的Jacobian矩阵：

$H_{gr,j,\mathrm{\&rho;}}=\frac{\&part;h_{gr,j}}{\&part;X_{\mathrm{\&rho;}}}{|}_{{\hat{X}}_{\mathrm{\&rho;}}\left(k\right|k-1)}=\left[\begin{array}{ccc}\frac{(x_j-x_g^{\left(1\right)})}{r_{g\_sat,j}^{\left(1\right)}}& \frac{(y_j-y_g^{\left(1\right)})}{r_{g\_sat,j}^{\left(1\right)}}& \frac{(z_j-z_g^{\left(1\right)})}{r_{g\_sat,j}^{\left(1\right)}}\\ & .& \\ & .& \\ & .& \\ \frac{(x_j-x_g^{\left(i\right)})}{r_{g\_sat,j}^{\left(i\right)}}& \frac{(y_j-x_g^{\left(i\right)})}{r_{g\_sat,j}^{\left(i\right)}}& \frac{(z_j-z_g^{\left(i\right)})}{r_{g\_sat,j}^{\left(i\right)}}\\ & .& \\ & .& \\ & .& \\ \frac{(x_j-x_g^{\left(n_j\right)})}{r_{g\_sat,j}^{\left(n_j\right)}}& \frac{(y_j-y_g^{\left(n_j\right)})}{r_{g\_sat,j}^{\left(n_j\right)}}& \frac{(z_j-z_g^{\left(n_j\right)})}{r_{g\_sat,j}^{\left(n_j\right)}}\end{array}\right]---\left(11\right)$

式中，待测卫星j在地心惯性系下的位置矢量[x_j,y_j,z_j]^T，地面站i在地心惯性系下的位置矢量为它到待测卫星j的距离为

位置估计误差ε_X的协方差矩阵为：

$Cov\left({\mathrm{\&epsiv;}}_X\right)=E\left({\mathrm{\&epsiv;}}_X{\mathrm{\&epsiv;}}_X^T\right)={({H_{gr,j,\mathrm{\&rho;}}}^T\&CenterDot;H_{gr,j,\mathrm{\&rho;}})}^{-1}{\mathrm{\&sigma;}}_{URE}^2---\left(12\right)$

式中，为地面站测量误差方差；令矩阵G＝(H_gr,j,ρ ^T·H_gr,j,ρ)^-1，则G由星站间的相对位置关系确定；

设位置精度因子如下：

$DOP=\sqrt{tr{({H_{gr,k,\mathrm{\&rho;}}}^T\&CenterDot;H_{gr,j,\mathrm{\&rho;}})}^{-1}}=\sqrt{tr\left(G\right)}---\left(13\right)$

采用最大单位矢端多面体体积法来确定地面站的最佳分布，首先求得卫星j到n个地面站的单位方向矢量然后计算由单位方向矢量(e₁～e_n)所构成的矢端多面体的体积，从所有组合中选取体积最大的一组作为地面站最佳分布方案；

利用该方法根据注入候选地面站信息在定轨解算之前确定最优地面站分布用于导航解算；若需更改地面站分布方案，则应重新注入候选地面站信息，并采用最大单位矢端多面体体积法重新选择地面站；

第二，GNSS选星策略

在进行GNSS选星时，无地面站测量的情况下采用方法(1)，在有地面站测量的情况下采用方法(2)，具体为：

(1)GNSS独立选星

当GNSS可见星数≤4时，利用全部GNSS可见星的测量信息；当可见星数>4时，从候选GNSS可见星中优选出4颗；

(2)基于地面站的GNSS选星

利用地面站的测量信息辅助GNSS选星，采用加权精度因子法进行GNSS选星；

定义加权精度因子：

$WDOP=\sqrt{tr{\left(H_{g\_gnss,j}^T{\mathrm{Weight}}_{g\_gnss,j}^{-1}{H}_{g\_gnss,j}\right)}^{-1}}---\left(14\right)$

式中，Η_{g_gnss,j}为将卫星j的地面站和星载GNSS接收机的测量函数向量在一步预测位置估计处进行线性化后所得到的Jacobian矩阵，即这里，H_gr,j,ρ为地面站伪距观测函数向量h_gr,j的Jacobian矩阵，参见式(11)；H_gnss,j,ρ为GNSS伪距观测函数向量h_gnss,j的Jacobian矩阵，即

$H_{gnss,j,\mathrm{\&rho;}}=\frac{\&part;h_{gnss,j}}{\&part;X_{\mathrm{\&rho;}}}{|}_{{\hat{X}}_{\mathrm{\&rho;}}\left(k\right|k-1)}=\left[\begin{array}{ccc}\frac{(x_j-x_{gnss}^{\left(1\right)})}{r_{gns\_sat,j}^{\left(1\right)}}& \frac{(y_j-y_{gnss}^{\left(1\right)})}{r_{gnss\_sat,j}^{\left(1\right)}}& \frac{(z_j-z_{gnss}^{\left(1\right)})}{r_{gnss\_sat,j}^{\left(1\right)}}\\ & .& \\ & .& \\ & .& \\ \frac{(x_j-x_{gnss}^{\left(s\right)})}{r_{gns\_sat,j}^{\left(s\right)}}& \frac{(y_j-y_{gnss}^{\left(s\right)})}{r_{gnss\_sat,j}^{\left(s\right)}}& \frac{(z_j-z_{gnss}^{\left(s\right)})}{r_{gnss\_sat,j}^{\left(s\right)}}\\ & .& \\ & .& \\ & .& \\ \frac{(x_j-x_{gnss}^{\left(l_j\right)})}{r_{gns\_sat,j}^{\left(l_j\right)}}& \frac{(y_j-y_{gnss}^{\left(l_j\right)})}{r_{gnss\_sat,j}^{\left(l_j\right)}}& \frac{(z_j-z_{gnss}^{\left(l_j\right)})}{r_{gnss\_sat,j}^{\left(l_j\right)}}\end{array}\right]---\left(15\right)$

式中，待测卫星j在地心惯性系下的位置矢量[x_j,y_j,z_j]^T；第s颗可见GNSS导航卫星在地心惯性系下的位置矢量为它到待测卫星j的距离为 $r_{gnss\_sat,j}^{\left(s\right)}=\sqrt{{(x_j-x_{gnss}^{\left(s\right)})}^2+{(y_j-y_{gnss}^{\left(s\right)})}^2+{(z_j-z_{gnss}^{\left(s\right)})}^2}.$

为卫星j的地面站和星载GNSS接收机伪距观测误差的加权阵，weight_g,j和weight_gnss,j分别为卫星j与地面站和星载GNSS接收机测量对应的对角加权阵，满足：

和设各地面站测量精度相当，测量噪声标准差为σ_g，那么和其中i＝1～n_j，s＝1～l_j；

在单独地面站、不同数量候选GNSS可见星的不同组合+地面站确定的WDOP值中，选取使WDOP值最小的GNSS卫星作为最佳GNSS选星结果；

步骤四：根据上述建立的基于地面站/星间链路/GNSS的GSO星座联合定轨系统模型，确定实现轨道参数估计的滤波方法；

采用EKF实现卫星轨道参数估计，对由式(6)和式(10)所描述的连续系统进行离散化可得：

$\left\{\begin{array}{c}X\left(k\right)=X(k-1)+F\&lsqb;X(k-1)\&rsqb;\mathrm{\&Delta;}t+A\&lsqb;X(k-1)\&rsqb;F\&lsqb;X(k-1)\&rsqb;{\mathrm{\&Delta;t}}^2/2+W(k-1)\\ Z\left(k\right)=h\&lsqb;X\left(k\right)\&rsqb;+v\left(k\right)\end{array}\right.---\left(16\right)$

其中，是状态向量；F是系统状态函数；W(k-1)为系统噪声，是零均值、方差为Q(k-1)的高斯白噪声；Δt为离散时间，为Jacobian矩阵；Z(k)是系统观测向量；h[X(k)]是系统观测函数向量；v(k)为测量噪声，是零均值、方差为R(k)的高斯白噪声向量；

对经过离散化后的非线性系统采用离散EKF方法进行状态估计，具体的滤波步骤如下：

(1)时间更新

计算一步预测状态：

$\hat{X}\left(k\right|k-1)=\hat{X}(k-1)+F\&lsqb;\hat{X}(k-1)\&rsqb;\&CenterDot;\mathrm{\&Delta;}t+A\&lsqb;\hat{X}(k+1)\&rsqb;\&CenterDot;F\&lsqb;\hat{X}(k-1)\&rsqb;\frac{{\mathrm{\&Delta;t}}^2}{2}---\left(17\right)$

式中，为(k-1)时刻的状态估计，为k时刻的一步预测状态估计；

计算状态转移矩阵：

$\mathrm{\&Phi;}\left(k\right|k-1)=I+A\&lsqb;\hat{X}(k-1)\&rsqb;\&CenterDot;\mathrm{\&Delta;}t---\left(18\right)$

式中，Φ(k|k-1)为从(k-1)时刻到k时刻的状态状态转移矩阵；

计算一步预测误差协方差矩阵：

P(k|k-1)＝Φ(k|k-1)P(k-1)Φ^T(k|k-1)+Q(k-1) (19)

式中，P(k|k-1)为一步预测状态的误差协方差矩阵，P(k-1)为(k-1)时刻估计状态的误差协方差矩阵，Q(k-1)为(k-1)时刻系统噪声方差阵；

(2)测量更新：

计算增益矩阵：

K(k)＝P(k|k-1)H^T(k)[H(k)P(k|k-1)H^T(k)+R(k)]^-1 (20)

式中，K(k)为k时刻的增益矩阵；为基于一步预测状态估计计算的观测向量h的Jacobian矩阵；R(k)为k时刻测量噪声协方差阵；

计算状态估计：

$\hat{X}\left(k\right)=\hat{X}\left(k\right|k-1)+K\left(k\right)\{Z\left(k\right)-h\&lsqb;\hat{X}\left(k\right|k-1)\&rsqb;\}---\left(21\right)$

式中，为k时刻的状态估计；

计算状态估计误差协方差阵：

P(k)＝[I-K(k)H(k)]P(k|k-1)[I-K(k)H(k)]^T+K(k)R(k)K^T(k) (22)

式中，P(k)为k时刻的状态估计误差的协方差阵；

步骤五：在集中式结构下，基于地面站/星间链路/GNSS的GSO星座联合定轨方法的具体实现；

根据定轨测量组合的优化策略自适应确定有效观测，利用EKF实现基于地面站/星间链路/GNSS的GSO星座联合定轨算法的具体过程如下：

(1)数据初始化，初始化的参数包括：k＝1，初始状态及误差协方差阵P(0)和系统噪声方差阵Q、初始地面站分布方案及参数；

(2)时间更新；利用(k-1)时刻的状态估计及误差协方差阵P(k-1)作为k时刻的状态初值及误差协方差阵，根据式(17)和式(19)分别计算状态一步预测值及预测误差协方差矩阵P(k|k-1),实现EKF估计方法的时间更新；

(3)初始化卫星编号j＝1；

(4)获取卫星j星载GNSS接收机观测数据，确定候选GNSS可见星数；

(5)获取卫星j的地面站观测数据，检测是否存在地面站测量，如果不存在则按顺序执行步骤(6)；否则按顺序执行步骤(7)；

(6)如果不存在地面站测量，则直接根据卫星j的候选GNSS可见星个数判断；当候选GNSS可见星数≥4时，从中确定4颗GNSS可见星的所有组合及相应的精度因子，选取最小精度因子所对应的可见星组合作为卫星j的GNSS的优选结果；当可见星数<4颗时，所有GNSS可见卫星均作为卫星j的GNSS的优选结果；跳到(10)；

(7)如果存在地面站测量，根据卫星j的地面站观测数据判断是否需要更新卫星j的地面站信息，如果需要，则注入新的可用地面站信息，利用最大单位矢端多面体体积法，在线确定卫星j的最优地面站分布；否则沿用k-1时刻卫星j选用的地面站分布；

(8)根据卫星j的地面站分布，利用卫星j获得的地面站有效测量信息，构建卫星j的地面站相关的测量模型和测量噪声方差阵；

(9)基于步骤(6)卫星j的地面站分布，根据卫星j的候选GNSS卫星的不同选星个数，计算卫星j在地面站和不同星载GNSS测量组合下的加权精度因子，从中选取最小加权精度因子所对应的可见星组合作为卫星j的GNSS优选结果；

(10)根据卫星j的GNSS可见星优选结果，利用所获得GNSS有效测量信息，构建卫星j的星载接收机相关的测量模型和测量噪声方差阵；

(11)判断是否星座内所有卫星的有效观测组合优选已结束；如果j＝m，继续下一步；否则j＝j+1，并返回到(4)；

(12)获取星间链路星间观测信息，构建测量模型和测量噪声方差阵中与星间链路测量相关的部分；

(13)利用所有可用测量信息对应的测量模型及测量噪声方差阵，组合构成k时刻的合成测量模型和测量噪声方差阵R(k)；

(14)测量更新；根据式(20)～(22)依次计算k时刻的增益矩阵K(k)、状态估计值及误差协方差矩阵P(k)，实现EKF的测量更新；

(15)判断滤波过程是否结束；如需继续进行滤波，k＝k+1，并返回到(2)。