CN105629272A - 短弧段批处理卫星自主定轨方法及装置 - Google Patents

Abstract

本发明公开了一种短弧段批处理卫星自主定轨方法及装置，该方法包括：采集当前观测弧段的GNSS观测数据，获得整个所述观测弧段的无粗差星间单差观测数据；积分预报所述观测弧段所有观测时刻的卫星状态信息；采用批处理方法迭代估计卫星状态信息；基于迭代估计的所述卫星状态信息计算验后残差和，并基于所述验后残差和判断所述迭代估计的敛散性，如果判断结果为收敛则执行步骤S15；积分预报下一观测弧段的卫星状态信息；根据所述GNSS观测数据计算接收机钟差，统计预报下一观测弧段的接收机钟差；将所述下一观测弧段的接收机钟差置入导航信号跟踪环路，并进入下一观测弧段。本发明提供的方法能够满足中低高卫星的轨道精度需求。

Description

短弧段批处理卫星自主定轨方法及装置

技术领域

本发明涉及卫星测控技术领域，具体涉及一种基于星载全球导航卫星系统(GlobalNavigationSatelliteSystem，GNSS)观测数据的短弧段批处理卫星自主定轨方法及装置。

背景技术

卫星自主轨道确定技术，简称卫星自主定轨，是指卫星不依赖于地面支持，而完全基于自载仪器设备，实时地确定自身位置和速度的技术。目前能够实现自主定轨的系统有三种：一是惯性导航系统。该系统具有隐蔽性强，自主性好等优点，缺点是随着卫星的长时间运行，系统误差不断积累，导致定位精度越来越差；二是天文定轨系统。该系统的缺点在于受限于天体的可见性；三是GNSS。该系统通过测量用户星与导航星之间的伪距或载波相位以及它们的变化率来确定用户星的位置和速度。与前两个系统相比，GNSS具有精度高，成本低，设备轻便等特点，已在低轨道卫星上得到了广泛的应用。

目前，在低轨卫星上广泛应用的卫星自主定轨方法主要有几何学轨道确定方法、基于几何法实时定轨结果的动力学滤波算法和基于伪距观测数据的滤波算法。前两种方法在可视导航卫星数不足四颗的情况下均无法完成自主定轨，第三种方法在可视导航卫星数长期不足四颗时定轨效果也十分不理想。

由于GNSS卫星均分布在中等高度的卫星轨道上，导航信号发射天线向地面方向发射，导致了中高轨卫星仅能接收来自地球对面的导航信号，首先，由于这种情况信号强度低，观测质量差，可视导航卫星数长期不足四颗，因此这些较差的观测条件限制了GNSS技术在中高轨卫星测控领域中的应用。其次，卫星为了完成一定的科学任务，需要具有三轴姿态大角度敏捷机动能力。在机动情况下，GNSS接收机可观测导航卫星数时常不足四颗，不能保证一直定位。最后，卫星工作时需要调整通信天线以对准中继卫星，完成数据传输等任务。在调整天线指向时，可能会对GNSS天线产生物理遮挡，影响GNSS接收机跟踪捕获的导航卫星数。

综上所述，现有的卫星自主定轨方法适用范围有限，不适用于观测条件较差的卫星，无法满足中高轨卫星和快速机动卫星的自主定轨需求。

发明内容

有鉴于此，为克服上述至少一个缺点，并提供下述至少一种优点。本发明公开了一种短弧段批处理卫星自主定轨方法及装置。

为解决上述技术问题，本发明采用以下技术方案：

一种短弧段批处理卫星自主定轨方法，包括：

步骤S11、采集当前观测弧段的GNSS观测数据，并基于所述GNSS观测数据获得整个所述观测弧段的无粗差星间单差观测数据；

步骤S12、以所述观测弧段终点时刻预报的卫星状态信息为先验信息，积分预报所述观测弧段所有观测时刻的卫星状态信息；

步骤S13、基于所述无粗差星间单差观测数据和积分预报的所述卫星状态信息，采用批处理方法迭代估计卫星状态信息；

步骤S14、基于迭代估计的所述卫星状态信息计算验后残差和，并基于所述验后残差和判断所述迭代估计的敛散性，如果判断结果为收敛则执行步骤S15；

步骤S15、基于所述迭代估计的卫星状态信息积分预报下一观测弧段的卫星状态信息；

步骤S16、根据所述GNSS观测数据计算接收机钟差，统计预报下一观测弧段的接收机钟差；

步骤S17、将所述下一观测弧段的接收机钟差置入导航信号跟踪环路，并进入下一观测弧段，执行步骤S11-S17。

如上所述的短弧段批处理卫星自主定轨方法，在一种可能的实现方式中，所述步骤S12中基于式(5)计算所述观测弧段所有观测时刻的卫星状态信息：

$\begin{array}{c}{Z}^i=h(X_0,T^i){|}_{X_0={X_0}^{\mathrm{ref}}}+\frac{\∂h(X_0,T_i)}{{\∂X}_0}{|}_{X_0={X_0}^{\mathrm{ref}}}(X_0-{X_0}^{\mathrm{ref}})+{\mathrm{\ϵ}}^i\\ =g\left(X^i\right){|}_{X^i=\mathrm{Int}({X_0}^{\mathrm{ref}},T^i)}+\frac{\∂g\left(X^i\right)}{{\∂X}^i}\frac{{\∂X}^i}{{\∂X}_0}{|}_{X^i=\mathrm{Int}({X_0}^{\mathrm{ref}},T^i),X_0={X_0}^{\mathrm{ref}}}(X_0-{X_0}^{\mathrm{ref}})+{\mathrm{\ϵ}}^i\end{array}---\left(5\right)$

其中，Zⁱ表示第i个观测时刻的观测数据，Tⁱ表示第i个观测时刻与所述当前弧段起点时刻的时间差，X₀表示所述当前观测弧段终点时刻的卫星的状态信息，X₀ ^ref表示参考初始状态信息，g(Xⁱ)及h(X₀,Tⁱ)表示第i个观测时刻根据预报的卫星状态信息计算出的观测数据，εⁱ为第i个观测时刻的观测噪声。

如上所述的短弧段批处理卫星自主定轨方法，在一种可能的实现方式中，

所述步骤S13中基于式(8)迭代估计卫星状态信息：

${{\hat{X}}_0}^{j+1}={\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{H_i}^T{\left(R^i\right)}^{-1}{H}_i\right)}^{-1}\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{H_i}^T{\left(R^i\right)}^{-1}{y}^i\right)+{\left({X_0}^{\mathrm{ref}}\right)}^j---\left(8\right)$

其中，表示第j+1次迭代估计的卫星初始状态信息，表示第j次迭代估计的卫星初始状态信息，E(ε^iTεⁱ)＝Rⁱ， $y^i=Z^i-{\stackrel{\‾}{g}}^i.$

如上所述的短弧段批处理卫星自主定轨方法，在一种可能的实现方式中，

所述步骤S14中基于式(9)计算验后残差和：

$U^{j+1}={({\hat{X}}_0-{\left({X_0}^{\mathrm{ref}}\right)}^j)}^T\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{H_i}^T{\left(R^i\right)}^{-1}{H}_i\right)({\hat{X}}_0-{\left({X_0}^{\mathrm{ref}}\right)}^i)---\left(9\right)$

其中，U^j+1表示第j+1次迭代估计的验后残差和。

如上所述的短弧段批处理卫星自主定轨方法，在一种可能的实现方式中，所述卫星状态信息包括卫星轨道信息，或所述卫星轨道信息和动力学参数，所述动力学参数包括大气阻力参数和太阳光压系数。

如上所述的短弧段批处理卫星自主定轨方法，在一种可能的实现方式中，所述基于所述GNSS观测数据获得整个所述弧段的无粗差的星间单差观测数据包括：

基于所述GNSS观测数据生成整个所述观测弧段的星间单差观测数据；

根据所预报的卫星的轨道信息对所述星间单差观测数据进行粗差探测和剔除，得到无粗差的星间单差观测数据。

一种短弧段批处理卫星自主定轨装置，包括：

采集模块，用于采集当前观测弧段的GNSS观测数据，并基于所述GNSS观测数据获得整个所述观测弧段的无粗差星间单差观测数据；

第一积分预报模块，用于以所述观测弧段终点时刻预报的卫星状态信息为先验信息，积分预报所述观测弧段所有观测时刻的卫星状态信息；

迭代估计模块，用于基于所述无粗差星间单差观测数据和积分预报的所述卫星状态信息，采用批处理方法迭代估计卫星状态信息；

判断模块，用于基于迭代估计的所述卫星状态信息计算验后残差和，并基于所述验后残差和判断所述迭代估计的敛散性；

第二积分预报模块，用于在所述判断模块的判断结果为收敛时基于所述迭代估计的卫星状态信息积分预报下一观测弧段的卫星状态信息；

钟差计算模块，用于根据所述GNSS观测数据计算接收机钟差，统计预报下一观测弧段的接收机钟差；以及

校正模块，用于将所述下一观测弧段的接收机钟差置入导航信号跟踪环路，并进入下一观测弧段。

如上所述的短弧段批处理卫星自主定轨装置，在一种可能的实现方式中，所述第一积分预报模块基于式(5)计算所述观测弧段所有观测时刻的卫星状态信息：

$\begin{array}{c}{Z}^i=h(X_0,T^i){|}_{X_0={X_0}^{\mathrm{ref}}}+\frac{\∂h(X_0,T_i)}{{\∂X}_0}{|}_{X_0={X_0}^{\mathrm{ref}}}(X_0-{X_0}^{\mathrm{ref}})+{\mathrm{\ϵ}}^i\\ =g\left(X^i\right){|}_{X^i=\mathrm{Int}({X_0}^{\mathrm{ref}},T^i)}+\frac{\∂g\left(X^i\right)}{{\∂X}^i}\frac{{\∂X}^i}{{\∂X}_0}{|}_{X^i=\mathrm{Int}({X_0}^{\mathrm{ref}},T^i),X_0={X_0}^{\mathrm{ref}}}(X_0-{X_0}^{\mathrm{ref}})+{\mathrm{\ϵ}}^i\end{array}---\left(5\right)$

其中，Zⁱ表示第i个观测时刻的观测数据，Tⁱ表示第i个观测时刻与所述当前弧段起点时刻的时间差，X₀表示所述当前观测弧段终点时刻的卫星的状态信息，X₀ ^ref表示参考初始状态信息，g(Xⁱ)及h(X₀,Tⁱ)表示第i个观测时刻根据预报的卫星状态信息计算出的观测数据，εⁱ为第i个观测时刻的观测噪声。

如上所述的短弧段批处理卫星自主定轨装置，在一种可能的实现方式中，所述迭代估计模块基于式(8)迭代估计卫星状态信息：

${{\hat{X}}_0}^{j+1}={\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{H_i}^T{\left(R^i\right)}^{-1}{H}_i\right)}^{-1}\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{H_i}^T{\left(R^i\right)}^{-1}{y}^i\right)+{\left({X_0}^{\mathrm{ref}}\right)}^j---\left(8\right)$

其中，表示第j+1次迭代估计的卫星初始状态信息，表示第j次迭代估计的卫星初始状态信息，E(ε^iTεⁱ)＝Rⁱ， $y^i=Z^i-{\stackrel{\‾}{g}}^i.$

如上所述的短弧段批处理卫星自主定轨装置，在一种可能的实现方式中，所述判断模块基于式(9)计算验后残差和：

$U^{j+1}={({\hat{X}}_0-{\left({X_0}^{\mathrm{ref}}\right)}^j)}^T\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{H_i}^T{\left(R^i\right)}^{-1}{H}_i\right)({\hat{X}}_0-{\left({X_0}^{\mathrm{ref}}\right)}^i)---\left(9\right)$

其中，U^j+1表示第j+1次迭代估计的验后残差和。

通过采用上述技术方案，本发明的所达到的有益效果为：首先，批处理算法能够适应较差观测条件，相比实时滤波算法，具有适用范围广，精度高，算法稳定等优点；本方法采用星间单差GNSS观测数据，在参数估计中消除了接收机钟差参数，待估参数只有卫星位置，速度，大气阻力系数和太阳光压系数等八个参数，能够适应星载处理器的性能要求；本方法采用龙哥库塔算法预报卫星轨道，相比其他批处理算法能够满足对卫星轨道信息的实时性需求；本方法根据卫星预报轨道和原始观测数据信息，统计接收机钟差变化规律，预报未来接收机钟差参数，并置入跟踪环路，保证GNSS接收机稳定输出高精度观测数据；本方法能够适应低轨，中轨和高轨等不同轨道高度卫星自主定轨应用，尤其在中高轨卫星自主定轨中具有明显优势。

附图说明

为了更清楚地说明本发明实施例中的技术方案，下面将对本发明实施例描述中所需要使用的附图作简单的介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据本发明实施例的内容和这些附图获得其他的附图。

图1为本发明一个实施例提供的短弧段批处理卫星自主定轨方法的流程图；

图2为本发明实施例提供的时序示意图；

图3为本发明另一实施例提供的短弧段批处理卫星自主定轨方法的流程图；

图4为本发明实施例提供的仿真实验示意图；

图5为本发明又一实施例提供的短弧段批处理卫星自主定轨装置的结构示意图。

具体实施方式

为使本发明解决的技术问题、采用的技术方案和达到的技术效果更加清楚，下面将结合附图对本发明实施例的技术方案作进一步的详细描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

下面结合附图并通过具体实施方式来进一步说明本发明的技术方案。

本发明人提出了一种基于星载GNSS观测数据的短弧段批处理卫星自主定轨方法及装置，使用整个观测弧段内的所有GNSS星间单差观测数据进行短弧段批处理，结合高精度卫星动力学模型，能够稳定且实时地输出高精度卫星轨道信息。相比传统卫星自主定轨方法，本方法适用于低轨，中轨和高轨等不同轨道高度的卫星自主定轨。

下面首先描述本发明实施例提供的的短弧段批处理卫星自主定轨方法的数学模型。

根据卫星的动力学模型可知，通过数值积分的方法，可以建立任意时刻T卫星质心的位置和速度与某一时刻卫星状态的联系，由式(1)来表示：

Xⁱ＝Int(X₀,T)(1)

其中，Xⁱ表示第i个观测时刻的的卫星状态信息，X₀表示在观测弧段起点时刻的卫星状态信息，T表示第i个观测时刻与观测弧段终点时刻的时间差，这里卫星状态信息指卫星的轨道信息，或者轨道信息以及大气动力学参数。通过该公式可以建立观测弧段内每个观测数据与观测弧段终点时刻卫星轨道信息和时间差的函数关系，由式(2)来表示：

Zⁱ＝g(Xⁱ)+ε＝h(X₀,Tⁱ)+ε(2)

其中，Zⁱ表示第i个观测时刻的观测数据，g(Xⁱ)表示第i个观测时刻根据预报的卫星状态信息，具体为轨道信息计算出的观测数据，h(X₀,T_i)同样表示第i个观测时刻根据预报的卫星轨道信息所计算出的观测数据，Tⁱ表示第i个观测时刻与观测弧段起点时刻的时间差，ε为观测噪声。且满足式(3)所表示的微分关系：

$\frac{\∂h(X_0,T^i)}{{\∂X}_0}=\frac{\∂g\left(X^i\right)}{{\∂X}^i}\frac{{\∂X}^i}{{\∂X}_0}---\left(3\right)$

此时观测方程如式(4)所示，

Zⁱ＝h(X₀,T_i)+εⁱ(4)

其中，εⁱ为观测弧段内第i个观测时刻的观测噪声，其统计特性为E(εⁱ)＝0，E(ε^iTεⁱ)＝Rⁱ，对观测方程(4)在参考初始状态信息处线性化，可得以下式(5)所表示的关系：

$\begin{array}{c}{Z}^i=h(X_0,T^i){|}_{X_0={X_0}^{\mathrm{ref}}}+\frac{\∂h(X_0,T_i)}{{\∂X}_0}{|}_{X_0={X_0}^{\mathrm{ref}}}(X_0-{X_0}^{\mathrm{ref}})+{\mathrm{\ϵ}}^i\\ =g\left(X^i\right){|}_{X^i=\mathrm{Int}({X_0}^{\mathrm{ref}},T^i)}+\frac{\∂g\left(X^i\right)}{{\∂X}^i}\frac{{\∂X}^i}{{\∂X}_0}{|}_{X^i=\mathrm{Int}({X_0}^{\mathrm{ref}},T^i),X_0={X_0}^{\mathrm{ref}}}(X_0-{X_0}^{\mathrm{ref}})+{\mathrm{\ϵ}}^i\end{array}---\left(5\right)$

为了简化公式(5)，令式(6)成立：

$\left\{\begin{array}{c}{\stackrel{\‾}{g}}^i=g\left(X^i\right){|}_{X^i=\mathrm{Int}({X_0}^{\mathrm{ref}},T^i)}\\ {\stackrel{\‾}{G}}^i=\frac{\∂g\left(X^i\right)}{{\∂X}^i}{|}_{X^i=\mathrm{Int}({X_0}^{\mathrm{ref}},T^i)}\\ \mathrm{\Φ}(t^i,t_0)=\frac{{\∂X}^i}{{\∂X}_0}{|}_{X_0={X_0}^{\mathrm{ref}}}\\ H_i={\stackrel{\‾}{G}}^i\mathrm{\Φ}\left(T^i\right)\end{array}\right.---\left(6\right)$

则线性化后的观测方程可简化为式(7)：

$\left[\begin{array}{c}.\\ .\\ .\\ Z^i\\ .\\ .\\ .\end{array}\right]=\left[\begin{array}{c}.\\ .\\ .\\ {\stackrel{\‾}{g}}^i+H_i(X_0-{X_0}^{\mathrm{ref}})+{\mathrm{\ϵ}}^i\\ .\\ .\\ .\end{array}\right]---\left(7\right)$

再令根据最小二乘原理，得到卫星初始状态信息的估计值

${{\hat{X}}_0}^{j+1}={\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{H_i}^T{\left(R^i\right)}^{-1}{H}_i\right)}^{-1}\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{H_i}^T{\left(R^i\right)}^{-1}{y}^i\right)+{\left({X_0}^{\mathrm{ref}}\right)}^j---\left(8\right)$

验后残差和以式(9)来表示：

$U^{j+1}={({\hat{X}}_0-{\left({X_0}^{\mathrm{ref}}\right)}^j)}^T\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{H_i}^T{\left(R^i\right)}^{-1}{H}_i\right)({\hat{X}}_0-{\left({X_0}^{\mathrm{ref}}\right)}^i)---\left(9\right)$

其中，上标j+1表示第j+1次迭代估计，U^j+1表示第j+1次迭代估计的验后残差和，表示第j+1次迭代估计的卫星初始状态信息，(X₀ ^ref)^j表示第j次迭代估计的卫星初始状态信息，第一次迭代估计时可以使用轨道机动前的状态信息初始化。

实施例1

如图1所示，为本发明一个实施例提供的短弧段批处理卫星自主定轨方法的流程图，结合图2所示的时序图，该方法包括如下步骤：

步骤S21、对于当前观测弧段，按照预先设定的观测数据采样频率采集原始GNSS观测数据，并基于GNSS观测数据生成星间单差观测数据，直至生成整个观测弧段的星间单差观测数据；然后对该星间单差观测数据进行粗差探测和剔除，得到不含粗差的星间单差观测数据。

本步骤中，采集原始GNSS观测数据时要保证每个采样点处观测的导航卫星数大于等于2。

步骤S22、以观测弧段终点时刻预报的卫星状态信息为先验信息，积分预报观测弧段内所有观测时刻的卫星轨道信息；

本实施例中，根据卫星的在观测弧段终点时刻的初始位置、初始速度所预报的卫星状态信息作为先验信息，结合高精度卫星动力学模型，基于线性化观测方程(式(5))向观测弧段起点方向进行积分，预报S21步骤中得到的无粗差星间单差观测数据对应的所有观测时刻的卫星状态信息；

$\begin{array}{c}{Z}^i=h(X_0,T^i){|}_{X_0={X_0}^{\mathrm{ref}}}+\frac{\∂h(X_0,T_i)}{{\∂X}_0}{|}_{X_0={X_0}^{\mathrm{ref}}}(X_0-{X_0}^{\mathrm{ref}})+{\mathrm{\ϵ}}^i\\ =g\left(X^i\right){|}_{X^i=\mathrm{Int}({X_0}^{\mathrm{ref}},T^i)}+\frac{\∂g\left(X^i\right)}{{\∂X}^i}\frac{{\∂X}^i}{{\∂X}_0}{|}_{X^i=\mathrm{Int}({X_0}^{\mathrm{ref}},T^i),X_0={X_0}^{\mathrm{ref}}}(X_0-{X_0}^{\mathrm{ref}})+{\mathrm{\ϵ}}^i\end{array}---\left(5\right)$

式中，Zⁱ表示第i个观测时刻的观测数据，h(X₀,T_i)表示根据预报的卫星轨道信息所计算出的第i个观测时刻的观测数据，X₀表示观测弧段起点时刻的卫星状态信息，Xⁱ表示在观测弧段终点时刻所预报的卫星状态信息，(X₀ ^ref)^j表示卫星初始状态信息，X₀ ^ref为参考初始状态信息，εⁱ为弧段内第i个观测数据的噪声，g(Xⁱ)为第i个观测时刻根据预报的卫星轨道信息计算出的观测数据。

步骤S23、基于步骤S21得到的无粗差观测数据，步骤S22积分预报的卫星状态信息，通过式(8)使用批处理方法迭代估计卫星状态信息。其中动力学参数包括大气阻力参数和太阳光压系数。

${{\hat{X}}_0}^{j+1}={\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{H_i}^T{\left(R^i\right)}^{-1}{H}_i\right)}^{-1}\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{H_i}^T{\left(R^i\right)}^{-1}{y}^i\right)+{\left({X_0}^{\mathrm{ref}}\right)}^j---\left(8\right)$

其中，表示第j+1次迭代估计的卫星初始状态信息，(X₀ ^ref)^j表示第j+1次迭代估计的卫星初始状态信息，E(ε^iTεⁱ)＝Rⁱ，

为了方便工程化实现，对卫星初始状态信息估计值公式(8)进一步简化，令

$\left\{\begin{array}{c}{N}^j=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{H}_i^T{\left(R^i\right)}^{-1}{H}_i\\ W^j=\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{H}_i^T{\left(R^i\right)}^{-1}{y}^i\end{array}\right.---\left(10\right)$

则，卫星初始状态信息的估计值，验后残差和可以以式(11)来表示：

$\left\{\begin{array}{c}\widehat{X_0^{j+1}}={\left(N^i\right)}^{-1}{W}^j+{\left(X_0^{\mathrm{ref}}\right)}^j\\ \widehat{P_0^{j+1}}={\left(N^j\right)}^{-1}\\ U^{j+1}={({\hat{X}}_0-{\left(X_0^{\mathrm{ref}}\right)}^j)}^T{N}^j({\hat{X}}_0-{\left(X_0^{\mathrm{ref}}\right)}^i)\end{array}\right.---\left(11\right)$

然后，每当采集到该观测弧段内的新的观测数据时，进行如下叠加操作：

$\left\{\begin{array}{c}{N}^j+=H_i^T{\left(R^i\right)}^{-1}{H}_t\\ W^j+=H_i^T{\left(R^i\right)}^{-1}{y}^i\end{array}\right.---\left(12\right)$

直到弧段内所有的观测数据都叠加完毕，利用公式(8)进行迭代估计，当相邻两次迭代所计算的验后残差和的差小于允许误差时可认为迭代收敛，结束迭代过程，输出初始时刻的卫星状态。

步骤S24、基于迭代估计的卫星状态信息，通过式(9)计算验后残差和，并基于验后残差和判断本次迭代估计的敛散性，如果收敛，则继续进行步骤S25的轨道确定，否则执行步骤S21，重新初始化。

$U^{j+1}={({\hat{X}}_0-{\left({X_0}^{\mathrm{ref}}\right)}^j)}^T\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{H_i}^T{\left(R^i\right)}^{-1}{H}_i\right)({\hat{X}}_0-{\left({X_0}^{\mathrm{ref}}\right)}^i)---\left(9\right)$

其中，U^j+1表示第j+1次迭代估计的验后残差和，表示第j+1次估计的卫星初始轨道信息，(X₀ ^ref)^j表示第j次估计的卫星初始轨道信息，第一次估计时可以使用轨道机动前的状态初始化(X₀ ^ref)⁰。

步骤S25、积分预报下一弧段的卫星轨道信息；

Xⁱ＝In(tX₀,Tⁱ)(1)

Xⁱ表示积分预报的第i个观测时刻的卫星状态信息，X₀表示观测弧段起点时刻的卫星状态信息，Tⁱ表示时刻t与弧段起点时刻的时间差。

步骤S26、根据观测弧段的原始观测值即GNSS观测数据计算每个采样时刻的接收机钟差，统计接收机钟差的一阶线性特征，并预报下一观测弧段起点时刻的接收机钟差。接收机钟差可以按照式(13)计算。

$\mathrm{clkBia}{s}^i=Z^i-{\hat{g}}^i\left(X^i\right){|}_{X^i=\mathrm{Int}({\hat{X}}_0,T^i)}---\left(13\right)$

其中，s为第i个观测时刻的接收机钟差，Zⁱ为观测数据，为最后一次迭代估计的观测弧段终点时刻的卫星状态信息。

步骤S27、基于所预报的接收机钟差对GNSS接收机跟踪环路的时钟进行校正；

本步骤中根据预报的钟差对接收机时钟进行校正，以保证接收机钟差不会漂移过大，保证正常卫星轨道的确定。

进入下一弧段，重复步骤S21-S27。

实施例2

如图3所示，为本发明另一实施例提供的短弧段批处理卫星自主定轨方法的流程图，该方法包括以下步骤：

步骤S31、开始并初始化系统；

步骤S32、采集当前观测弧段的GNSS观测数据，并判断是否到达当前观测弧段的边界，是则执行步骤S33，否则执行步骤S42；

步骤S33、判断当前观测弧段是否已经初始化，是则执行步骤S34，否则执行步骤S35；

步骤S34、以观测弧段终点时刻预报的卫星状态信息作为先验信息，并初始化观测弧段起点时刻的轨道信息，初始化法矩阵和误差矩阵，进入步骤S36；

其中法矩阵N和误差矩阵W即式(10)表示的矩阵。

步骤S35、初始化观测弧段轨道信息，初始化法矩阵和误差矩阵，执行步骤S43；

步骤S36、判断当前观测弧段是否还有未处理的GNSS观测数据，是则执行步骤S37，否则执行步骤S38；

步骤S37、计算当前时刻的观测方程系数矩阵以及O-C矩阵，执行步骤S39；

其中O-C矩阵即式(5)所表示的矩阵。

步骤S38、更新GNSS观测数据，并进行迭代估计，进入步骤S40；

步骤S39、叠加法矩阵和误差矩阵，进入步骤S36；

步骤S40、计算验后残差和，并基于验后残差和判断迭代估计的敛散性，判断结果为发散则执行步骤S41，判断结果为收敛则执行步骤S32；

步骤S41、输出错误状态的提示信息，结束自主定轨；

步骤S42、采集当前观测弧段的GNSS观测数据并存储，执行步骤S32；

步骤S43、判断当前观测弧段是否还有未处理的GNSS观测数据，是则执行步骤S44，否则执行步骤S45；

步骤S44、计算当前时刻的观测方程系数矩阵以及O-C矩阵，执行步骤S46；

步骤S45、更新GNSS观测数据，并进行迭代估计，执行步骤S47；

步骤S46、叠加法矩阵和误差矩阵，执行步骤S43；

步骤S47、计算验后残差和，并基于验后残差和判断迭代估计的敛散性，判断结果为发散则执行步骤S48，判断结果为收敛则执行步骤S32；

步骤S48、判断是否达到迭代极限，是则执行步骤S49，否则执行步骤S35；

步骤S49、输出错误状态的提示信息，结束自主定轨。

下面结合特定工程实施案例，给出对本发明提供的方法的仿真验证方案和验证结果。如图4所示，为了验证本发明提供的方法的有效性，采用GNSS信号源进行验证，验证方案步骤如下：

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

步骤S42、在PC机上，使用基于GNSS观测数据的短弧段批处理卫星自主定轨方法对步骤S41中得到的GNSS接收机输出的广播星历和原始观测数据进行模拟自主定轨数据处理，输出高精度的定轨计算结果；

步骤S43、将步骤S42中得到的模拟自主定轨数据处理的解算结果与GNSS信号仿真器输出的理论轨道数据进行比对，完成用GNSS仿真实测数据对基于星载GNSS多天线的卫星自主定轨方案的精度评估与验证。

验证结果表明：本发明提供的短弧段批处理卫星自主定轨方法能够适用于中高轨卫星和快速机动卫星的自主定轨应用，稳定且实时地输出高精度卫星轨道信息。相比传统卫星自主定轨方法，本方法适用范围更广，适用于中、低以及高轨道卫星，具有重大的工程实用价值和军事战略意义。

实施例3

如图5所示，本发明实施例还提供了一种短弧段批处理卫星自主定轨装置，该装置100包括：采集模块10、第一积分预报模块20、迭代估计模块30、判断模块40、第二积分预报模块50、钟差计算模块60以及校正模块70。其中采集模块10用于采集当前观测弧段的GNSS观测数据，并基于GNSS观测数据获得整个观测弧段的无粗差星间单差观测数据；第一积分预报模块20用于以观测弧段终点时刻预报的卫星状态信息为先验信息，积分预报所述观测弧段所有观测时刻的卫星状态信息；迭代估计模块30用于基于无粗差星间单差观测数据和积分预报的所述卫星状态信息，采用批处理方法迭代估计卫星状态信息；判断模块40用于基于迭代估计的所述卫星状态信息计算验后残差和，并基于所述验后残差和判断迭代估计的敛散性；第二积分预报模块50用于在判断模块40的判断结果为收敛时基于所述迭代估计的卫星状态信息积分预报下一观测弧段的卫星状态信息；钟差计算模块60用于根据所述GNSS观测数据计算接收机钟差，统计预报下一观测弧段的接收机钟差；校正模块70用于将下一观测弧段的接收机钟差置入导航信号跟踪环路，并进入下一观测弧段。

其中第一积分预报模块20基于式(5)计算所述观测弧段所有观测时刻的卫星状态信息：

$\begin{array}{c}{Z}^i=h(X_0,T^i){|}_{X_0={X_0}^{\mathrm{ref}}}+\frac{\∂h(X_0,T_i)}{{\∂X}_0}{|}_{X_0={X_0}^{\mathrm{ref}}}(X_0-{X_0}^{\mathrm{ref}})+{\mathrm{\ϵ}}^i\\ =g\left(X^i\right){|}_{X^i=\mathrm{Int}({X_0}^{\mathrm{ref}},T^i)}+\frac{\∂g\left(X^i\right)}{{\∂X}^i}\frac{{\∂X}^i}{{\∂X}_0}{|}_{X^i=\mathrm{Int}({X_0}^{\mathrm{ref}},T^i),X_0={X_0}^{\mathrm{ref}}}(X_0-{X_0}^{\mathrm{ref}})+{\mathrm{\ϵ}}^i\end{array}---\left(5\right)$

其中，Zⁱ表示第i个观测时刻的观测数据，Tⁱ表示第i个观测时刻与所述当前弧段起点时刻的时间差，X₀表示所述当前观测弧段终点时刻的卫星的状态信息，X₀ ^ref表示参考初始状态信息，g(Xⁱ)及h(X₀,Tⁱ)表示第i个观测时刻根据预报的卫星状态信息计算出的观测数据，εⁱ为第i个观测时刻的观测噪声。

迭代估计模块30基于式(8)迭代估计卫星状态信息：

${{\hat{X}}_0}^{j+1}={\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{H_i}^T{\left(R^i\right)}^{-1}{H}_i\right)}^{-1}\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{H_i}^T{\left(R^i\right)}^{-1}{y}^i\right)+{\left({X_0}^{\mathrm{ref}}\right)}^j---\left(8\right)$

其中，表示第j+1次迭代估计的卫星初始状态信息，表示第j次迭代估计的卫星初始状态信息，E(ε^iTεⁱ)＝Rⁱ， $y^i=Z^i-{\stackrel{\‾}{g}}^i.$

判断模块40基于式(9)计算验后残差和：

$U^{j+1}={({\hat{X}}_0-{\left({X_0}^{\mathrm{ref}}\right)}^j)}^T\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{H_i}^T{\left(R^i\right)}^{-1}{H}_i\right)({\hat{X}}_0-{\left({X_0}^{\mathrm{ref}}\right)}^i)---\left(9\right)$

其中，U^j+1表示第j+1次迭代估计的验后残差和。

本发明提供的短弧段批处理卫星自主定轨方法及装置，与现有技术相比，有益效果包括以下几个方面：首先，批处理算法能够适应较差观测条件，相比实时滤波算法，具有适用范围广，精度高，算法稳定等优点；本方法采用星间单差GNSS观测数据，在参数估计中消除了接收机钟差参数，待估参数只有卫星位置，速度，大气阻力系数和太阳光压系数等八个参数，能够适应星载处理器的性能要求；本方法采用龙哥库塔算法预报卫星轨道，相比其他批处理算法能够满足对卫星轨道信息的实时性需求；本方法根据卫星预报轨道和原始观测数据信息，统计接收机钟差变化规律，预报未来接收机钟差参数，并置入跟踪环路，保证GNSS接收机稳定输出高精度观测数据；本方法能够适应低轨，中轨和高轨等不同轨道高度卫星自主定轨应用，尤其在中高轨卫星自主定轨中具有明显优势。

以上实施例提供的技术方案中的全部或部分内容可以通过软件编程实现，其软件程序存储在可读取的存储介质中，存储介质例如：计算机中的硬盘、光盘或软盘。

注意，上述仅为本发明的较佳实施例及所运用技术原理。本领域技术人员会理解，本发明不限于这里所述的特定实施例，对本领域技术人员来说能够进行各种明显的变化、重新调整和替代而不会脱离本发明的保护范围。因此，虽然通过以上实施例对本发明进行了较为详细的说明，但是本发明不仅仅限于以上实施例，在不脱离本发明构思的情况下，还可以包括更多其他等效实施例，而本发明的范围由所附的权利要求范围决定。

Claims (10)

1.一种短弧段批处理卫星自主定轨方法，其特征在于，包括：

步骤S11、采集当前观测弧段的GNSS观测数据，并基于所述GNSS观测数据获得整个所述观测弧段的无粗差星间单差观测数据；

步骤S12、以所述观测弧段终点时刻预报的卫星状态信息为先验信息，积分预报所述观测弧段所有观测时刻的卫星状态信息；

步骤S13、基于所述无粗差星间单差观测数据和积分预报的所述卫星状态信息，采用批处理方法迭代估计卫星状态信息；

步骤S14、基于迭代估计的所述卫星状态信息计算验后残差和，并基于所述验后残差和判断所述迭代估计的敛散性，如果判断结果为收敛则执行步骤S15；

步骤S15、基于所述迭代估计的卫星状态信息积分预报下一观测弧段的卫星状态信息；

步骤S16、根据所述GNSS观测数据计算接收机钟差，统计预报下一观测弧段的接收机钟差；

步骤S17、将所述下一观测弧段的接收机钟差置入导航信号跟踪环路，并进入下一观测弧段，执行步骤S11-S17。

2.如权利要求1所述的短弧段批处理卫星自主定轨方法，其特征在于，所述步骤S12中基于式(5)计算所述观测弧段所有观测时刻的卫星状态信息：

$\begin{array}{c}{Z}^i=h(X_0,T^i){|}_{X_0={X_0}^{\mathrm{ref}}}+\frac{\∂h(X_0,T_i)}{{\∂X}_0}{|}_{X_0={X_0}^{\mathrm{ref}}}(X_0-{X_0}^{\mathrm{ref}})+{\mathrm{\ϵ}}^i\\ =g\left(X^i\right){|}_{X^i=\mathrm{Int}({X_0}^{\mathrm{ref}},T^i)}+\frac{\∂g\left(X^i\right)}{{\∂X}^i}\frac{{\∂X}^i}{{\∂X}_0}{|}_{X^i=\mathrm{Int}({X_0}^{\mathrm{ref}},T^i),X_0={X_0}^{\mathrm{ref}}}(X_0-{X_0}^{\mathrm{ref}})+{\mathrm{\ϵ}}^i\end{array}---\left(5\right)$

其中，Zⁱ表示第i个观测时刻的观测数据，Tⁱ表示第i个观测时刻与所述当前弧段起点时刻的时间差，X₀表示所述当前观测弧段终点时刻的卫星的状态信息，X₀ ^ref表示参考初始状态信息，g(Xⁱ)及h(X₀,Tⁱ)表示第i个观测时刻根据预报的卫星状态信息计算出的观测数据，εⁱ为第i个观测时刻的观测噪声。

3.如权利要求2所述的短弧段批处理卫星自主定轨方法，其特征在于，所述步骤S13中基于式(8)迭代估计卫星状态信息：

${{\hat{X}}_0}^{j+1}={\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{H_i}^T{\left(R^i\right)}^{-1}\right)}^{-1}\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{H_i}^T{\left(R^i\right)}^{-1}{y}^i\right)+{\left({X_0}^{\mathrm{ref}}\right)}^j---\left(8\right)$

其中，表示第j+1次迭代估计的卫星初始状态信息，表示第j次迭代估计的卫星初始状态信息，E(ε^iTεⁱ)＝Rⁱ， $y^i=Z^i-{\stackrel{-}{g}}^i.$

4.如权利要求3所述的短弧段批处理卫星自主定轨方法，其特征在于，所述步骤S14中基于式(9)计算验后残差和：

$U^{j+1}={({\hat{X}}_0-{\left({X_0}^{\mathrm{ref}}\right)}^j)}^T\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{H_i}^T{\left(R^i\right)}^{-1}{H}_i\right)({\hat{X}}_0-{\left({X_0}^{\mathrm{ref}}\right)}^i)---\left(9\right)$

其中，U^j+1表示第j+1次迭代估计的验后残差和。

5.如权利要求1所述的短弧段批处理卫星自主定轨方法，其特征在于，所述卫星状态信息包括卫星轨道信息，或所述卫星轨道信息和动力学参数，所述动力学参数包括大气阻力参数和太阳光压系数。

6.如权利要求1所述的短弧段批处理卫星自主定轨方法，其特征在于，所述基于GNSS观测数据获得整个所述观测弧段的无粗差的星间单差观测数据包括：

基于所述GNSS观测数据生成整个所述观测弧段的星间单差观测数据；

根据所预报的卫星的轨道信息对所述星间单差观测数据进行粗差探测和剔除，得到无粗差的星间单差观测数据。

7.一种短弧段批处理卫星自主定轨装置，其特征在于，包括：

采集模块，用于采集当前观测弧段的GNSS观测数据，并基于所述GNSS观测数据获得整个所述观测弧段的无粗差星间单差观测数据；

第一积分预报模块，用于以所述观测弧段终点时刻预报的卫星状态信息为先验信息，积分预报所述观测弧段所有观测时刻的卫星状态信息；

迭代估计模块，用于基于所述无粗差星间单差观测数据和积分预报的所述卫星状态信息，采用批处理方法迭代估计卫星状态信息；

判断模块，用于基于迭代估计的所述卫星状态信息计算验后残差和，并基于所述验后残差和判断所述迭代估计的敛散性；

第二积分预报模块，用于在所述判断模块的判断结果为收敛时基于所述迭代估计的卫星状态信息积分预报下一观测弧段的卫星状态信息；

钟差计算模块，用于根据所述GNSS观测数据计算接收机钟差，统计预报下一观测弧段的接收机钟差；以及

校正模块，用于将所述下一观测弧段的接收机钟差置入导航信号跟踪环路，并进入下一观测弧段。

8.如权利要求7所述的短弧段批处理卫星自主定轨装置，其特征在于，所述第一积分预报模块基于式(5)计算所述观测弧段所有观测时刻的卫星状态信息：

$\begin{array}{c}{Z}^i=h(X_0,T^i){|}_{X_0={X_0}^{\mathrm{ref}}}+\frac{\∂h(X_0,T_i)}{{\∂X}_0}{|}_{X_0={X_0}^{\mathrm{ref}}}(X_0-{X_0}^{\mathrm{ref}})+{\mathrm{\ϵ}}^i\\ =g\left(X^i\right){|}_{X^i=\mathrm{Int}({X_0}^{\mathrm{ref}},T^i)}+\frac{\∂g\left(X^i\right)}{{\∂X}^i}\frac{{\∂X}^i}{{\∂X}_0}{|}_{X^i=\mathrm{Int}({X_0}^{\mathrm{ref}},T^i),X_0={X_0}^{\mathrm{ref}}}(X_0-{X_0}^{\mathrm{ref}})+{\mathrm{\ϵ}}^i\end{array}---\left(5\right)$

其中，Zⁱ表示第i个观测时刻的观测数据，Tⁱ表示第i个观测时刻与所述当前弧段起点时刻的时间差，X₀表示所述当前观测弧段终点时刻的卫星的状态信息，X₀ ^ref表示参考初始状态信息，g(Xⁱ)及h(X₀,Tⁱ)表示第i个观测时刻根据预报的卫星状态信息计算出的观测数据，εⁱ为第i个观测时刻的观测噪声。

9.如权利要求7所述的短弧段批处理卫星自主定轨装置，其特征在于，所述迭代估计模块基于式(8)迭代估计卫星状态信息：

${{\hat{X}}_0}^{j+1}={\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{H_i}^T{\left(R^i\right)}^{-1}\right)}^{-1}\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{H_i}^T{\left(R^i\right)}^{-1}{y}^i\right)+{\left({X_0}^{\mathrm{ref}}\right)}^j---\left(8\right)$

其中，表示第j+1次迭代估计的卫星初始状态信息，表示第j次迭代估计的卫星初始状态信息，E(ε^iTεⁱ)＝Rⁱ， $y^i=Z^i-{\stackrel{-}{g}}^i.$

10.如权利要求7所述的短弧段批处理卫星自主定轨装置，其特征在于，所述判断模块基于式(9)计算验后残差和：

$U^{j+1}={({\hat{X}}_0-{\left({X_0}^{\mathrm{ref}}\right)}^j)}^T\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{H_i}^T{\left(R^i\right)}^{-1}{H}_i\right)({\hat{X}}_0-{\left({X_0}^{\mathrm{ref}}\right)}^i)---\left(9\right)$

其中，U^j+1表示第j+1次迭代估计的验后残差和。