CN107153209B - 一种短弧段低轨导航卫星实时精密定轨方法 - Google Patents

Abstract

本发明公开了一种短弧段低轨导航卫星实时精密定轨方法，该方法利用四台或四台以上地面接收机同时跟踪若干颗GNSS(全球导航卫星系统)导航卫星和LEO(低轨)导航卫星的导航信号，并实时发送到数据中心，包括这些步骤：步骤一，观测数据的获取与数据预处理；步骤二，各个接收机间精密时间同步；步骤三，利用天顶方向对流层延迟和投影函数计算LEO视线方向的对流层延迟；步骤四，修正接收机硬件延迟，形成LEO定轨几何观测方程；步骤五，联合LEO的动力学模型和几何观测方程求解LEO精密轨道和精密钟差；步骤六，LEO精密轨道和精密钟差的实时分发。采用本发明的方法能够解决低轨卫星精密轨道无法实时计算、需事后处理计算的问题。

Description

一种短弧段低轨导航卫星实时精密定轨方法

技术领域

本发明属于卫星导航技术领域，具体涉及一种短弧段低轨导航卫星实时精密定轨方法。

背景技术

低轨卫星和微小卫星以其低廉的发射和制造成本，获得了越来越多的青睐。低轨卫星在通信、遥感、资源管理以及空间基准维持等方面做出了巨大的贡献，其应用也从传统的对地观测扩展到导航领域，因而越来越多的低轨卫星发射计划被提上日程。著名的铱星计划通过66颗低轨卫星完成了基于低轨卫星的全球卫星通信系统的组网，并于2016年宣布其卫星授时和位置服务(STL)的计划，宣告导航卫星进入了低轨卫星时代。国内方面，武汉大学正在研制的低轨卫星珞珈一号搭载了天基导航载荷，也开始了低轨导航卫星的探索性研究。

实时获取精确的卫星轨道是实现精密卫星导航定位的关键技术。目前低轨卫星轨道确定方法主要包括三种：地面监测站定轨、星上自主定轨和事后精密定轨。目前GNSS系统，包括GPS，GLONASS，BDS等都使用少数的地面跟踪站实现轨道全弧段观测，然后计算和预报卫星轨道，上注到导航卫星。该方法对中高轨卫星较理想，但是低轨卫星轨道低，地面可见区域小，要实现全弧段跟踪通常需要全球均匀分布的数十个甚至上百个跟踪站，受地理条件限制，难以实现。星上自主定轨是依靠星上的GNSS接收机或惯性测量单元进行轨道确定，其中基于GNSS观测值的星上自主定轨方法能够提供实时，自主连续的LEO轨道，因而应用最广泛。然而受GNSS导航卫星广播星历的轨道误差和钟差的影响，基于GNSS观测值的LEO星上自主定轨通常只能获得数米量级的定轨精度。事后精密定轨是利用低轨卫星下传的星载GNSS观测数据和精密的GNSS轨道和钟差信息在地面求解LEO的精密轨道，也可以加入地面激光测卫(SLR)的数据联合解算。后处理LEO精密定轨精度可达2-5cm，但是由于星上观测数据下传窗口有限，难以做到实时下传。另外后处理通常使用批处理最小二乘方式确定轨道，该算法精度较高，但仅适用于后处理或近实时处理。鉴于以上原因，后处理LEO精密定轨无法满足实时定轨需求。

综上所述，现有的低轨卫星定轨方法无法提供实时的高精度LEO轨道信息，无法满足低轨卫星导航实时获取高精度轨道的需求。

发明内容

针对现有技术中存在的问题，为了提供实时的低轨导航卫星的精密轨道和时钟信息，用以支持低轨导航卫星的精密定位应用，本发明公开了一种短弧段低轨导航卫星实时精密定轨方法。

一种短弧段低轨导航卫星实时精密定轨方法，该方法利用四台或四台以上地面接收机同时跟踪若干颗GNSS(全球导航卫星系统)导航卫星和LEO(低轨)导航卫星的导航信号，并实时发送到数据中心，包括以下步骤：

步骤一，观测数据的获取与数据预处理；

步骤二，各个接收机间精密时间同步；

步骤三，利用天顶方向对流层延迟和投影函数计算LEO视线方向的对流层延迟；

步骤四，修正接收机硬件延迟，形成LEO定轨几何观测方程；

步骤五，联合LEO的动力学模型和几何观测方程求解LEO精密轨道和精密钟差；

步骤六，LEO精密轨道和精密钟差的实时分发。

进一步地，所述导航信号包括GPS，GLONASS，BDS，Galileo，IRNSS，QZSS等卫星导航系统，导航信号内容包括伪距码观测值、载波相位观测值和多普勒观测值。

进一步地，所述时间同步采用绝对时间同步，其方法是数据中心利用各接收机接收到的导航卫星信号和从外部数据流解码获得的精密卫星轨道和钟差信息，约束地面接收机坐标，采用滤波的方式解算精密的接收机钟差、天顶方向对流层延迟和模糊度参数。

进一步地，所述绝对时间同步的解算方法如下：

导航卫星的伪距和载波相位观测值表示为：

式(1)中：P_i和φ_i分别表示第i个频点的伪距和载波相位观测值(以米为单位)，ρ表示卫星和接收机间的几何距离，δt^S和δt^R分别表示卫星和接收机的时钟偏差(以秒为单位)，I_i和δ_trop分别表示电离层和对流层延迟影响，N_i表示第i个频点的整周模糊度，ε_Pi和ε_φi分别表示伪距和载波相位的接收机观测噪声，c表示真空中光速；λ_i表示第i个频点载波波长(以米为单位)；

时间同步的观测模型表示为：

式(2)中：E(·)为GNSS观测值向量，C为钟差参数的系数矩阵，T为一个包含对流层投影函数的列向量作为天顶对流层延迟参数的系数矩阵，Λ为模糊度参数的系数矩阵；

考虑到实时估计的需求和参数变化的非线性特征，参数估计使用扩展的卡尔曼滤波(EKF)实现；

考虑到接收机钟差较好的短期稳定度和长期漂移特性，滤波过程中接收机钟差参数使用随机游走过程模拟，根据接收机钟的质量和稳定度施加合适的过程噪声；

采用随机游走过程对对流层天顶方向湿分量影响进行模拟，对流层延迟量与导航信号传播路径有关，在参数化时利用对流层投影函数将各个接收机卫星视线方向的对流层延迟投影至测站天顶方向，记作

δ_trop＝δ_trop,z*MF (3)

式(3)中：δ_trop表示信号视线方向对流层延迟，δ_trop,z表示天顶方向对流层延迟，MF表示对流层投影函数，选择一种投影函数实现对流层延迟从导航卫星信号视线方向到测站天顶方向的投影；

模糊度参数在接收机载波环正常锁定的情况下当作常数，在载波环抖动或者失锁的情况下，模糊度参数重置；

滤波完成后，获得各个接收机钟面时和GNSS系统时之间的偏差以及各个接收机天顶方向的对流层延迟。

进一步地，所述对流层投影函数MF包括Marini函数、Chao函数、Herring函数、Neil投影函数(NMF)、Global投影函数(GMF)、Vienna投影函数(VMF)，还包括自定义投影函数和根据区域跟踪网计算的投影函数。

进一步地，利用已知的各个接收机天线相位中心坐标和近似的LEO卫星轨道计算LEO卫星在各个接收机处的站心地平极坐标系坐标，即各个接收机观测到的LEO卫星信号的视线方向；利用测站坐标、信号发射时间以及LEO卫星信号的视线方向，计算对流层投影函数，并利用该投影函数将上述估计得到的天顶方向对流层延迟投影到LEO卫星信号的视线方向，获得对流层延迟对LEO导航信号的延迟量；使用实验室标定、在线标定或等效时间基准的方式处理硬件延迟引起的时间系统偏差，将时间同步的参考时间基准从GNSS系统时转换到LEO系统时。

进一步地，利用各接收机接收到的LEO导航信号列出LEO轨道确定的几何观测方程，LEO导航信号观测值表示为：

式(4)中：LEO导航信号的卫星端钟差作为参数估计，LEO导航信号的接收机端钟差通过参数估值和硬件延迟校正计算得到，电离层延迟误差I_i通过两个LEO导航频点组成无电离层组合消除，几何距离项ρ按照泰勒级数一阶展开，可得:

式(5)中：ρ₀为站星几何距离近似值，dx、dy、dz表示近似LEO轨道的改正数，x^S,y^S,z^S表示卫星三维坐标，x^R,y^R,z^R表示接收机天线相位中心三维坐标，ε是泰勒展开的高阶项。

进一步地，联合LEO动力学模型和几何观测求解轨道钟差的过程如下：

考虑到低轨卫星在空间运动，其受力情况表示为：

式(6)中：r，分别表示卫星位置、速度和加速度矢量，GM为地球重力场常数，f_p为扰动加速度，Q₁,…,Q_d为d个待估动力学参数，动力学参数包括大气阻力系数C_d、太阳光压系数C_r和伪随机脉冲参数等；

通过线性化的方法求解出加速度对各个参数的偏导数矩阵，得到状态转移矩阵Φ(t,t₀)，使用状态转移矩阵进行时间更新：

式(7)中：P_t ^-分别为t时刻时间更新后的状态量及其方差协方差矩阵，分别为(t-1)时刻轨道状态滤波解及其方差协方差矩阵，Q(t)为过程噪声矩阵。

进一步地，联合动力学信息和几何观测信息进行轨道确定的卡尔曼滤波解可以表示为

上式中：K_t为t时刻滤波增益矩阵，R_t为几何观测值方差协方差矩阵。y为根据LEO导航信号计算得到的几何观测信息向量。

进一步地，计算出滤波解后，还要对几何观测值的验后残差进行检验。滤波方程的验后残差估值可表示为：

根据验后残差估值可构造χ²检验量，服从χ²(n-t)分布，其中(n-t)表示自由度；若则认为验后残差通过检验，否则认为验后残差检验拒绝，式中是显著性水平为α的阈值；如果验后残差检验未通过，则使用轨道预报解作为精密定轨输出；如果验后残差检验通过，则使用滤波解作为精密定轨输出。

进一步地，将上述计算得到的精密轨道和精密钟差编码成实时数据流格式，如RTCM，通过互联网或者无线电台的方式播发给用户，用户就可以使用收到的精密轨道和钟差来计算；考虑到定轨计算和信息传播的延迟，需要对收到的精密轨道和钟差外推到用户信号接收时刻，再参与地面定位解算。

与现有的方法相比，本方法具有以下优点：

1.能够实时获取低轨卫星短弧段的精密轨道和钟差，产品的延迟时间主要是网络数据传输造成的，仅为数秒，而该数据延迟造成的误差影响可以通过数据外推有效的补偿。

2.地面跟踪设备价格低廉，仅需数台定制化的GNSS接收机，无需精密卫星跟踪设备即可获得高精度的卫星轨道。

3.监测网布设成本低，无需使用光纤、专用电缆等方法实现多接收机精密时间同步。

4.使用递推滤波算法，滤波状态参数维度低，计算效率高。

5.该定轨方法不受观测条件、气象等因素影响，有较强的环境适应能力。

6.该定轨方法同时适用于一颗LEO定轨和多颗LEO同时定轨的情况。

7.该定轨方法不要求LEO卫星播发的导航信号与GNSS导航信号精确同步，也不要求多颗LEO卫星间导航信号时间系统精确同步。

附图说明

图1为本发明提出的短弧段低轨导航卫星实时精密定轨方法的原理示意图。

图2为本发明提出的短弧段低轨导航卫星实时精密定轨方法的流程图。

图3为本发明实施例提供的短弧段低轨导航卫星实时精密定轨方法中滤波算法流程图。

具体实施方式

下面将结合附图以及具体实施例来详细说明本发明，其中的具体实施例以及说明仅用来解释本发明，但并不作为对本发明的限定。

本发明提出的短弧段低轨导航卫星实时精密定轨方法原理如图1所示。图1中的MEO，GEO，IGSO分别指中轨道、地球同步轨道和倾斜同步轨道的导航卫星，包括但不限于GPS，GLONASS，BDS，Galileo，IRNSS，QZSS等全球或区域导航卫星系统的卫星；LEO指具有导航测距信号发射功能的低轨道卫星。为了保证低轨道卫星能够播发稳定的、高精度的测距信号，低轨卫星应配备驯服后的高稳晶振或原子钟。图1中的接收机是指能够同时接收导航卫星和LEO卫星导航信号、并能够实时输出高精度伪距和载波相位的测量型接收机。接收机应保证将LEO卫星导航信号与GNSS导航信号同步采样。接收机应具有网络通信功能或连接至具有网络通信功能的数据转发设备。图1中的数据中心指至少1台具有网络通信功能的高性能服务器，并安装了使用本发明论述的短弧段低轨导航卫星实时精密定轨方法实现的软件。为了保证低轨卫星的轨道和钟差可估，要求地面必须配备4台或4台以上的接收机同步观测导航卫星和低轨卫星的导航信号，形成地面监测网络。考虑到低轨卫星导航信号的地面覆盖范围直径只有数千公里，地面接收机之间的间距建议选在1000km左右来获得LEO轨道较好的几何观测条件和较长的LEO公共可见弧段。各个接收机对应天线的相位中心坐标需要事先精确测定，并且坐标应保证和GNSS精密轨道在同一参考框架下。考虑到低轨卫星发射的导航信号频点与导航卫星不同，接收机内部对不同频率信号的硬件延迟也不同。对于单个可见的LEO来说，接收机射频元器件的硬件延迟偏差需要预先标定。各接收机应保证在LEO卫星过境前、中、后连续稳定的工作，并且实时地将导航卫星、LEO卫星的双频伪距和载波相位观测值实时地、连续地发送至数据中心。图1所示的各卫星(包括GEO/IGSO/MEO/LEO)发射导航或测距信号，地面接收机负责导航信号的接收，数据中心负责轨道计算和分发。本发明涉及的计算方法均在数据中心执行。

本发明提出的短弧段低轨导航卫星实时精密定轨方法的详细流程如图2所示，下面对该流程图进行详细论述。

步骤101：观测数据获取与数据预处理；数据中心对收到的各个接收机的发过来的观测数据进行解码、校验，确认无误后逐个接收机进行数据预处理。数据预处理方法包括使用无几何距离组合等方式进行数据扫描、检测和标记异常数据及可能存在的载波相位周跳。

步骤102：接收机间精密时间同步；受到接收机钟的频率准确度和稳定度的限制，各个接收机的时钟不是精确同步的，因此需要首先计算各个接收机钟的钟面时差异来进行时间同步。同步的方式有两种，相对时间同步和绝对时间同步。相对时间同步是将所有接收机的钟面时同步到某一个钟的钟面时。绝对时间同步是将所有接收机的钟面时同步到外部的系统时间，如GPS时或BDS时。由于使用相对时间同步会导致最终确定LEO钟差与外部系统时之间存在偏差，建议使用绝对时间同步。绝对时间同步方法是指数据中心利用已知的各个接收机天线相位中心坐标，从外部数据流解码获得导航卫星精密的卫星轨道和卫星钟差信息，采用滤波的方式解算精密的接收机钟差，天顶方向对流层延迟和模糊度参数。解算方法如下：

导航卫星的伪距和载波相位观测值可以表示为：

P_i＝ρ+δ_orb+c(δt^S-δt^R)+I_i+δ_trop+ε_Pi

φ_i＝ρ+δ_orb+c(δt^S-δt^R)-I_i+δ_trop+λ_iN_i+ε_φi

上式中：P_i和φ_i分别表示第i个频点的伪距和载波相位观测值(以米为单位)；表示卫星和接收机间的几何距离，其中x^S,y^S,z^S表示卫星三维坐标，x^R,y^R,z^R表示接收机天线相位中心三维坐标；δt^S和δt^R分别表示卫星和接收机的时钟偏差(以秒为单位)；I_i和δ_trop分别表示电离层和对流层延迟的影响；N_i表示第i个频点的整周模糊度；ε_Pi和ε_φi分别表示伪距和载波相位的接收机观测噪声；c表示真空中光速；λ_i表示第i个频点的载波波长(以米为单位)。

在步骤102中，接收机的天线相位中心坐标x^R,y^R,z^R可以事先精确测定，卫星坐标x^S,y^S,z^S可以通过精密星历，卫星天线偏差文件精确计算，目前超快速GPS卫星轨道误差优于5cm。在数据处理中，几何距离项ρ可以视为已知值，也可以考虑接收机天线相位中心坐标的不确定性，将接收机坐标作为未知参数进行估计，但同时根据其先验的精度施加一个较强的约束。对于精密星历，轨道误差δ_orb可忽略。δt^S可通过实时估计的精密钟差数据流改正。电离层延迟I_i可通过无电离层组合消除，剩下的误差项，包括接收机钟差δt^R，对流层延迟δ_trop和模糊度参数N_i通过参数估计的方法处理。时间同步的观测模型可表示为：

上式中：E(·)为GNSS观测值向量，C为钟差参数的系数矩阵，矩阵列数与参与计算的卫星系统数相等。T为一个包含对流层投影函数的列向量作为天顶对流层延迟参数的系数矩阵。Λ为模糊度参数的系数矩阵，对于载波相位观测值，该矩阵为对角阵，每个对角线元素为对应载波相位观测值的波长；对于伪距观测值，该矩阵为零矩阵。如果载波相位观测值不参与计算，则不考虑载波相位参数及其系数矩阵。

考虑到实时估计的需求和参数变化的非线性特征，参数估计使用扩展的卡尔曼滤波(EKF)实现。考虑到接收机钟差较好的短期稳定度和长期漂移特性，滤波过程中接收机钟差参数使用随机游走过程模拟，根据接收机钟的质量和稳定度施加合适的过程噪声。常见的GNSS接收机钟包括(但不限于)温补石英晶振(TCXO)、恒温石英晶振(OCXO)、原子钟，对于OCXO和原子钟长期漂移缓慢，在一定时间内也可以当作随机常数估计。由于接收机硬件延迟与信号频率有关，对不同频率的信号接收机端硬件延迟并不相同。如果接收机能够同时跟踪多个GNSS系统的信号，应考虑不同GNSS系统信号间的硬件延迟不同。因此，每个GNSS系统应单独估计一个接收机钟差参数。以GPS/BDS双模接收机为例，滤波模型中接收机钟差向量应表示为由于GPS和BDS信号均采用同一个晶振，在接收机工作环境不变的情况下，和的变化趋势一致。因此，接收机钟差的另一种估计算法将接收机钟差向量表示为该模型中仍采用随机游走模型估计，可采用随机常数模型估计，该项也被称为系统间偏差(ISB)。

根据中性大气理论，对流层对电磁波延迟影响可以分为干分量和湿分量。干分量约占延迟的90％，并且可以用经验模型很好的计算，湿分量和气压，大气水汽分布有关，难以模型化，通常需要通过参数化的方式进行估计。对流层延迟量与导航信号传播路径有关，在参数化时可以利用对流层投影函数将各个接收机卫星视线方向的对流层延迟投影至测站天顶方向，记作

δ_trop＝δ_trop,z*MF

上式中：δ_trop表示信号视线方向对流层延迟，δ_trop,z表示天顶方向对流层延迟，MF表示对流层投影函数，常见的投影函数包括但不限于Marini函数、Chao函数、Herring函数、Neil投影函数(NMF)、Global投影函数(GMF)、Vienna投影函数(VMF)。选择一种投影函数实现对流层延迟从导航卫星信号视线方向到测站天顶方向的投影。

在步骤102中采用随机游走过程对对流层天顶方向湿分量影响进行模拟。模糊度参数在接收机载波环正常锁定的情况下可以当作常数，在滤波器中也作为随机常数估计。在载波环抖动或者失锁的情况下，模糊度参数需要重置。滤波完成后，可获得各个接收机钟面时和GNSS系统时(如GPS时，BDS时)之间的偏差，以及各个接收机天顶方向的对流层延迟。

步骤103：计算LEO视线方向对流层延迟；利用精确已知的各个接收机天线相位中心坐标和近似的LEO卫星轨道计算LEO卫星在各个接收机处的站心地平极坐标系坐标，即各个接收机观测到的LEO卫星信号的视线方向。利用测站坐标，信号发射时间以及LEO卫星信号的视线方向，计算对流层投影函数，并利用该投影函数将步骤102估计得到的天顶方向对流层延迟投影到LEO卫星信号的视线方向，获得对流层延迟对LEO导航信号的延迟量。

步骤104：形成LEO定轨几何观测方程；步骤102中时间同步是将各接收机的钟面时同步到参考GNSS系统(如GPS系统)的系统时，对于LEO轨道计算需要使用LEO系统时作为参考。对于接收机而言，LEO系统时和参考GNSS系统的系统时之间差异由系统时间差和接收机硬件延迟两部分组成。系统时间差可根据两时间系统的系统定义计算，硬件延迟偏差是由于不同导航信号之间频率不同导致的，该偏差与信号频率，接收机射频元器件，接收机硬件设计，固件，接收机工作环境有关。接收机硬件延迟是可以使用三种方法处理：实验室标定，在线标定和等效时间基准。实验室标定指在实验室内测定接收机射频元器件和硬件链路对不同频率输入的响应时延。在线标定是使用其他方式获得的精密LEO轨道和钟差反算接收机硬件延迟。等效时间基准是对整个监测网使用同种型号接收机时，认为各接收机引起的硬件延迟相同。根据基准变换(S变换)原理，各接收机相同的硬件延迟被自动吸收到估计的LEO卫星钟差中，不会对LEO轨道确定精度产生影响。标定参考GNSS系统与LEO系统的时间差后，再根据步骤102估计的GNSS信号的接收机钟差来推算相对于LEO导航信号的接收机钟差，以单GPS接收机为例，推算公式为：

上式中：可从步骤102的估值中获取，为预先在实验室标定的两频点不同引起的硬件延迟偏差。为LEO轨道确定需要使用的接收机钟差。

若接收机为多模GNSS接收机，如GPS/BDS接收机，则需要选择某一卫星导航系统(如GPS系统)为主卫星导航系统，并在实验室内测定所有接收机LEO导航信号与该主卫星导航系统频率之间的硬件延迟偏差。再根据步骤102估计得到的该主卫星导航系统的接收机钟差来确定接收机相对于LEO导航信号的接收机钟差。

利用各接收机接收到的LEO导航信号列出LEO轨道确定的几何观测方程，LEO导航信号观测值可表示为：

由于轨道参数是待估参数，因此不再将轨道误差列入观测误差源。式中的几何距离项ρ按照泰勒级数一阶展开，可得:

上式中：ρ₀为站星几何距离近似值，可利用预先测定的测站坐标和通过广播星历计算的近似LEO轨道计算。dx，dy，dz表示近似LEO轨道的改正数，列入轨道估计滤波器的状态向量。ε是泰勒展开的高阶项，计算中通常忽略不计。LEO导航信号的卫星端钟差也作为参数估计。LEO导航信号的接收机端钟差可通过步骤102中的参数估值和硬件延迟校正计算得到。电离层延迟误差I_i可通过两个LEO导航频点组成无电离层组合消除，或者使用电离层约束模型进行估计。考虑到大部分LEO的轨道高度低于电离层总厚度，因此LEO卫星播发的导航信号受到电离层影响小于MEO，GEO和IGSO导航卫星。在状态向量中每一个估计一个电离层参数，并且在随机模型中给该电离层参数设置一个合理的先验方差信息，这样可以有效的加快滤波收敛速度。考虑到所有地面接收机的观测值，上式可表示为线性系统

y＝Ax

其中位置参数向量x＝[dx,dy,dz,dt]，y为几何观测值，A为几何观测值线性化后的设计矩阵。

类似的，也可以利用多普勒观测值将卫星的速度矢量和接收机钟变化率表达成同样的形式。

步骤105：联合LEO动力学模型和几何观测求解轨道钟差；考虑到低轨卫星在空间运动受力情况可表示为：

上式中：r，分别表示卫星位置、速度和加速度矢量，GM为地球重力场常数，f_p为扰动加速度，Q₁,…,Q_d为d个待估动力学参数。动力学待估参数通常包括大气阻力系数C_d，太阳光压系数C_r和伪随机脉冲参数。r为卫星质心到地球质心间几何距离。通过线性化的方法求解出加速度对各个参数的偏导数矩阵，称为状态转移矩阵Φ(t,t₀)。考虑到计算复杂度，状态转移矩阵计算通常做一定的近似，常用的状态转移矩阵计算方法有Goodyear法，Markley法。求解出状态转移矩阵后，即可使用状态转移矩阵进行时间更新：

上式中：P_t ^-分别为t时刻时间更新后的状态量及其方差协方差矩阵，分别为(t-1)时刻轨道状态滤波解及其方差协方差矩阵。Q(t)为过程噪声矩阵。考虑到线性化后的几何观测方程，轨道确定的卡尔曼滤波解可以表示为

K_t＝P_t ^-A(AP_t ^-A^T+R_t)^-1

P_t ⁺＝(I-K_tA)P_t ^-

上式中：K_t为t时刻滤波增益矩阵，R_t为几何观测值方差协方差矩阵。计算出滤波解后，还需要对几何观测值的验后残差进行检验，如果检验未通过，则仅使用时间更新的解作为最终定轨结果；如果检验通过，则使用滤波解作为最终定轨结果。然后开始下一历元轨道计算。

步骤106：LEO精密轨道和钟差的实时分发；将步骤105计算得到的精密轨道和精密钟差编码成实时数据流格式，如RTCM，通过互联网或者无线电台的方式播发给用户，用户就可以使用收到的精密轨道和钟差来计算。考虑到定轨计算和信息传播的延迟，通常需要对收到的精密轨道和钟差外推到用户信号接收时刻，再参与地面定位解算。

本发明与现有后处理精密定轨的显著区别之一就是本发明使用卡尔曼滤波算法确定实时轨道，而不是批处理最小二乘算法。滤波流程已在步骤105中有所涉及，下面对滤波算法进行更详细的介绍。本发明实施例提供的短弧段低轨导航卫星实时精密定轨方法中滤波算法流程如图3所示。

步骤201:初始化轨道状态向量。初始化轨道状态向量指滤波开始时刻的近似的LEO卫星位置矢量r和速度矢量卫星初始的轨道参数可以使用LEO广播的轨道。如果LEO播发的信号不包括广播星历，也可以使用步骤101～步骤104计算得到的信息使用逆向的单点定位计算求解初始卫星轨道。卫星初始轨道通常只需要数米级的精度，因此普通的星上自主定轨即可满足要求。由于本发明使用的是Kalman滤波计算，轨道计算会有一个收敛过程，而收敛过程时间长短与近似坐标准确程度有关。

步骤202：轨道积分；根据LEO卫星的初始轨道状态，分析卫星的受力情况，再利用卫星的受力与运动的关系通过积分的方式预报卫星速度和位置。自由飞行的LEO在轨卫星受到力主要包括重力、多体引力、太阳辐射压、大气阻力、地球辐射压、地球潮汐、未模型化的力。其中重力可分解为正常重力和重力摄动项。正常重力可以使用已知的地球常量直接计算，地球的摄动项可以展开表示成一定阶次的球谐函数来逼近。目前已有多种重力场模型给出了各自的球谐函数系数，如JGM3，EGM96，EGM2008等。应用这些重力场模型，选择合适的阶次，即可计算出相应精度的重力场模型。米级精密定轨通常选用30～50阶次的重力场模型，厘米级精密定轨则选用90～120阶次的重力场模型。模型阶次越高，计算量越大，计算耗时越长。多体引力主要是使用星表，如JPL DE200，JPL DE405星表，计算其他星体惯性系坐标，再根据相应星体的天文常数确定该星体对卫星的吸引力。米级定轨通常只考虑日月引力，精密定轨还考虑太阳系其他行星的吸引力。太阳辐射压和卫星表面材质，卫星与太阳光线的夹角等有关。米级定轨使用经验的反射系数，精密定轨需要估计卫星三个面的反射系数C_R。卫星进入地影区不受太阳辐射压影响。大气阻力主要和卫星轨道处的大气密度，卫星表面阻尼系数，卫星面质比和卫星速度有关。其中大气密度通常使用经验模型，如Harris-Priester模型、Jacchia-Lineberry模型计算，卫星三个方向的面质比可根据卫星参数确定，表面阻尼系数C_D在米级定轨中使用经验参数，在精密定轨中作为参数估计。地球辐射压、地球潮汐可根据模型计算，未模型化的力在米级定轨中忽略，在精密定轨中使用经验力或者伪随机脉冲吸收。

LEO卫星某一时刻受到的力等于以上各个力的合力。根据牛顿运动定律，可根据卫星的受力计算卫星的加速度，进而积分得到卫星的速度。由于卫星的受力非常复杂，难以使用解析的方法进行积分计算，通常都使用数值积分的方法逼近，常见的单步数值积分方法有4阶Runge-Kutta-Fehlberg方法(RKF4)，4阶Runge-Kutta方法(RK4)等，也有使用多步积分、变阶次、变步长等更复杂的数值积分方法，如Stoermer-Cowell积分法。通过数值积分的方法，就可以根据初始轨道x₀预报其他时刻的卫星轨道

步骤203：求解状态转移矩阵；状态转移矩阵是当前的卫星状态和卫星加速度对所有参数计算一阶导数。重力场加速度(含摄动项)对卫星轨道向量[x，y，z]的偏导数可用Cunningham方法计算。多体引力、太阳辐射压、大气阻力等对滤波状态向量的偏导数均可使用公式直接计算。类似地，状态转移矩阵也需要先计算某一时刻的状态转移矩阵后，再使用数值积分的方法进行递推到指定时刻。由于状态转移矩阵计算牵涉导数计算，通常在计算中做一定程度的近似，因此状态转移矩阵递推没有步骤202中论述的数值积分精度高，通常用来方差协方差矩阵的时间更新计算。

步骤204：处理几何观测值；几何观测值处理主要包括利用已知的地面站坐标和近似LEO轨道计算近似站星几何距离，改正地面接收机钟差，对流层延迟，消除电离层延迟等方面，详细的计算方法在步骤102和步骤103中阐述。

步骤205：列几何观测方程；LEO卫星的几何观测量可使用泰勒级数展开，取一阶项，近似转化为线性系统，具体的列几何观测方程的方法在步骤104中阐述。

步骤206：滤波求解轨道；联合几何观测和动力学模型，使用扩展的Kalman滤波(EKF)方法求解LEO的精密轨道、速度和钟差。具体的滤波方法在步骤105中阐述。

步骤207：验后残差检验；滤波方程的验后残差估值可表示为：

根据验后残差估值可构造χ²检验量，服从χ²(n-t)分布，其中(n-t)表示自由度。若则认为验后残差通过检验，否则认为验后残差检验拒绝。式中是显著性水平为α的阈值。如果验后残差检验未通过，则使用轨道预报解作为精密定轨输出；如果验后残差检验通过，则使用滤波解作为精密定轨输出。

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

Claims (10)

1.一种短弧段低轨导航卫星实时精密定轨方法，该方法利用四台或四台以上地面接收机同时跟踪若干颗GNSS(全球导航卫星系统)导航卫星和LEO(低轨)导航卫星的导航信号，并实时发送到数据中心，其特征在于，包括以下步骤：

步骤一，观测数据的获取与数据预处理；

步骤二，各个接收机间精密时间同步；

步骤三，利用天顶方向对流层延迟和投影函数计算LEO视线方向的对流层延迟；

步骤四，修正接收机硬件延迟，形成LEO定轨几何观测方程；

步骤五，联合LEO的动力学模型和几何观测方程求解LEO精密轨道和精密钟差；

步骤六，LEO精密轨道和精密钟差的实时分发。

2.根据权利要求1所述的短弧段低轨导航卫星实时精密定轨方法，其特征在于：所述导航信号包括GPS，GLONASS，BDS，Galileo，IRNSS，QZSS卫星导航系统的导航信号，导航信号内容包括伪距码观测值、载波相位观测值和多普勒观测值。

3.根据权利要求2所述的短弧段低轨导航卫星实时精密定轨方法，其特征在于：所述时间同步采用绝对时间同步，其方法是数据中心利用各接收机接收到的导航卫星信号和从外部数据流解码获得的精密卫星轨道和钟差信息，约束地面接收机坐标，采用滤波的方式解算精密的接收机钟差、天顶方向对流层延迟和模糊度参数。

4.根据权利要求3所述的短弧段低轨导航卫星实时精密定轨方法，其特征在于：所述绝对时间同步的解算方法如下：

导航卫星的伪距和载波相位观测值表示为：

式(1)中：P_i和φ_i分别表示第i个频点的伪距和载波相位观测值，以米为单位，ρ表示卫星和接收机间的几何距离，δ_orb表示轨道误差，δt^S和δt^R分别表示卫星和接收机的时钟偏差，以秒为单位，I_i和δ_trop分别表示电离层和对流层延迟影响，N_i表示第i个频点的整周模糊度，ε_Pi和ε_φi分别表示伪距和载波相位的接收机观测噪声，c表示真空中光速；λ_i表示第i个频点载波波长，以米为单位；

时间同步的观测模型表示为：

式(2)中：E(·)为GNSS观测值向量，C为钟差参数的系数矩阵，T为一个包含对流层投影函数的列向量作为天顶对流层延迟参数的系数矩阵，Λ为模糊度参数的系数矩阵；

考虑到实时估计的需求和参数变化的非线性特征，参数估计使用扩展的卡尔曼滤波(EKF)实现；

考虑到接收机钟差较好的短期稳定度和长期漂移特性，滤波过程中接收机钟差参数使用随机游走过程模拟，根据接收机钟的质量和稳定度施加合适的过程噪声；

采用随机游走过程对对流层天顶方向湿分量影响进行模拟，对流层延迟量与导航信号传播路径有关，在参数化时利用对流层投影函数将各个接收机卫星视线方向的对流层延迟投影至测站天顶方向，记作

δ_trop＝δ_trop,z*MF (3)

式(3)中：δ_trop表示信号视线方向对流层延迟，δ_trop,z表示天顶方向对流层延迟，MF表示对流层投影函数，选择一种投影函数实现对流层延迟从导航卫星信号视线方向到测站天顶方向的投影；

模糊度参数在接收机载波环正常锁定的情况下当作常数，在载波环抖动或者失锁的情况下，模糊度参数重置；

滤波完成后，获得各个接收机钟面时和GNSS系统时之间的偏差以及各个接收机天顶方向的对流层延迟。

5.根据权利要求4所述的短弧段低轨导航卫星实时精密定轨方法，其特征在于：所述对流层投影函数MF包括Marini函数、Chao函数、Herring函数、Neil投影函数(NMF)、Global投影函数(GMF)、Vienna投影函数(VMF)，还包括自定义投影函数和区域投影函数。

6.根据权利要求4所述的短弧段低轨导航卫星实时精密定轨方法，其特征在于：利用已知的各个接收机天线相位中心坐标和近似的LEO卫星轨道计算LEO卫星在各个接收机处的站心地平极坐标系坐标，即各个接收机到LEO卫星的视线方向；利用测站坐标、信号发射时间以及LEO卫星信号的视线方向，计算对流层投影函数，并利用该投影函数将上述估计得到的天顶方向对流层延迟投影到LEO卫星信号的视线方向，获得对流层延迟对LEO导航信号的延迟量；使用实验室标定、在线标定或等效时间基准的方式处理硬件延迟引起的时间系统偏差，将时间同步的参考时间基准从GNSS系统时转换到LEO系统时。

7.根据权利要求6所述的短弧段低轨导航卫星实时精密定轨方法，其特征在于：利用各接收机接收到的LEO导航信号列出LEO轨道确定的几何观测方程，LEO导航信号观测值表示为：

式(4)中：LEO导航信号的卫星端钟差作为参数估计，LEO导航信号的接收机端钟差通过参数估值和硬件延迟校正计算得到，电离层延迟误差I_i通过两个LEO导航频点组成无电离层组合消除，几何距离项ρ按照泰勒级数一阶展开，可得:

式(5)中：ρ₀为站星几何距离近似值，dx、dy、dz表示近似LEO轨道的改正数，x^S,y^S,z^S表示卫星三维坐标，x^R,y^R,z^R表示接收机天线相位中心三维坐标，ε是泰勒展开的高阶项。

8.根据权利要求7所述的短弧段低轨导航卫星实时精密定轨方法，其特征在于：联合LEO动力学模型和几何观测求解轨道钟差的过程如下：

考虑到低轨卫星在空间运动，其受力情况表示为：

式(6)中：r，分别表示卫星位置、速度和加速度矢量，GM为地球重力场常数，f_p为扰动加速度，Q₁,…,Q_d为d个待估动力学参数，动力学参数包括大气阻力系数C_d、太阳光压系数C_r和伪随机脉冲参数；

通过线性化的方法求解出加速度对各个参数的偏导数矩阵，得到状态转移矩阵Φ(t,t₀)，使用状态转移矩阵进行时间更新：

式(7)中：P_t-分别为t时刻时间更新后的状态量及其方差协方差矩阵，分别为(t-1)时刻轨道状态滤波解及其方差协方差矩阵，Q(t)为过程噪声矩阵。

9.根据权利要求8所述的短弧段低轨导航卫星实时精密定轨方法，其特征在于：计算出滤波解后，还要对几何观测值的验后残差进行检验；滤波方程的验后残差估值表示为：

根据验后残差估值构造χ²检验量；服从χ²(n-t)分布，其中(n-t)表示自由度；若则认为验后残差通过检验，式中是显著性水平为α的阈值；如果验后残差检验未通过，则使用轨道预报解作为精密定轨输出；如果验后残差检验通过，则使用滤波解作为精密定轨输出。

10.根据权利要求9所述的短弧段低轨导航卫星实时精密定轨方法，其特征在于：将上述计算得到的精密轨道和精密钟差编码成实时数据流格式，通过互联网或者无线电台的方式播发给用户。