一种低轨卫星实时定轨方法与流程

本发明涉及实时定位技术，更具体地说，尤其涉及一种低轨卫星实时定轨方法。

背景技术：

当前，诸如高精度测绘、导航增强等领域都需要低轨卫星的高精度定位数据，随着需求的不断扩大，越来越多的应用对定位的实时性和高精度都提出了要求。目前国内外学者进行的研究大部分都是采用以下两种方法：方法一是以测码伪距或侧相位距作为观测量，使用广播星历计算导航星位置和钟差，再结合简化的动力学方程，使用扩展卡尔曼滤波器进行计算；方法二除了将观测量换成载波相位以外，其它步骤与方法一大体相同。但是它们都有各自的不足之处：方法一定位的精度不理想；方法二由于涉及到模糊度的解算，状态向量维度大大增加，在星载计算机上应用的难度很大。

本发明将有效地解决上述缺陷，面向低轨卫星的实际工程应用，充分考虑精度要求以及星载计算机计算能力的限制，解决兼顾计算精度与计算效率的星载算法关键问题。。

技术实现要素：

本发明的目的在于克服现有技术的不足，提供一种兼顾计算精度与计算效率的低轨卫星实时定轨方法，在对星载gnss接收机观测量进行数据预处理的同时，使用地面提前上注的igs超快速预测星历代替精度不高的广播星历来计算导航卫星的位置和钟差，结合低轨卫星的动力学方程，使用扩展卡尔曼滤波器对目标卫星进行定位。

本发明是根据以下技术方案实现的：

一种低轨卫星实时定轨方法，其特征在于，包括以下步骤：

步骤s1：从星载gnss接收机获取待定位历元的双频测码伪距和载波相位观测量；

步骤s2：对获取的双频测码伪距和载波相位观测量分别进行误差修正，包括：电离层延迟、相对论效应、天线相位中心误差、地球自转效应；

步骤s3：对修正误差以后的载波相位观测量进行周跳和粗差的实时探测与修复；

步骤s4：在使用经过处理的载波相位观测值对测码伪距观测量进行平滑；

步骤s5：使用地面上注的igs超快速星历计算导航星的位置和钟差信息；

步骤s6：将上一历元的定位结果作为初值，使用扩展卡尔曼滤波器计算当前历元的低轨卫星定位结果；

步骤s7:判断定位任务是否结束，若没有结束，则返回步骤s1继续进行计算。

上述技术方案中，步骤s3具体包括：

步骤s301：使用m-w组合计算整周模糊度互差项nδ

其中，n1,n2分别为f1,f2频率上的整周模糊度，l1，l2分别为f1,f2频率上的载波相位观测值，ρ1，ρ2分别为f1,f2频率上的测码伪距观测值；

步骤s302：计算前i个历元nδ的平均值＜nδ＞i及其均方根误差σi，为了保证方法的实时性，采用以下方式计算：

其中，nδi为第i个历元的整周模糊度互差项；

步骤s303：使用上述公式计算历元i的整周模糊度互差项nδi、历元i+1的整周模糊度互差项nδ(i+1)以及前i-1个历元整周模糊度互差项的均值＜nδ＞i-1以后，通过|nδi-＜nδ＞i-1|≥σ和|nδ(i+1)-nδi|＜1探测周跳和粗差，如果这两个判断条件同时满足，认为在历元i-1和i之间已经发生了周跳；如果有|nδi-＜nδ＞i-1|≥σ和|nδ(i+1)-nδi|≥1同时成立，认为历元i时刻的相位突变是粗差造成的，公式中的σ为根据星载gnss接收机的实际情况设定好的阈值；

步骤s304：如果判断出来有周跳存在，进行如下计算：

接着分别在周跳发生前后的两个历元联立nδ和w两个方程，分别解出周跳发生前后历元在f1,f2两个频率上的整周模糊度n1,n2，进而解出周跳发生前后历元整周模糊度的变化量，也就是周跳，接着对周跳发生之后的连续观测历元中对观测量进行修正；

如果判断出有粗差存在，则将粗差对应的观测量进行剔除。

上述技术方案中，所述步骤s4进一步包括：如果出现意外情况星上无法获取超快速星历，则使用广播星历计算导航星位置和钟差。

上述技术方案中，步骤s6进一步包括：滤波器中解算的状态向量包括三维位置向量、三维速度向量以及接收机钟差，以平滑后的测码伪距作为观测量，首先使用上一历元的定位结果或事先设定的初值作为卡尔曼滤波器的初始状态向量；然后使用4阶龙格库塔法对低轨卫星的精密动力学方程进行数值积分得到位置、速度向量的一步预测值，同时把接收机钟差当作白噪声模型来计算一步预测值；最后使用经过载波相位平滑以后的测码伪距观测值对状态向量的一步预测值进行修正。

与现有技术相比，本发明具有如下的有益效果：

本发明可以在星上计算资源消耗较小的前提下，完成较高精度的低轨卫星定位工作，可以适用于低轨卫星星载gnss接收机的实时定位。

附图说明

通过阅读参照以下附图对非限制性实施例所作的详细描述，本发明的其它特征、目的和优点将会变得更明显：

图1为本发明的低轨卫星实时定位方法的流程图。

具体实施方式

下面结合具体实施例对本发明进行详细说明。以下实施例将有助于本领域的技术人员进一步理解本发明，但不以任何形式限制本发明。应当指出的是，对本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变化和改进。这些都属于本发明的保护范围。

如图1所示，一种低轨卫星实时定轨方法，它包括以下步骤：

步骤s1：从星载gnss接收机获取待定位历元的双频测码伪距和载波相位观测量；

步骤s2：对获取的双频测码伪距和载波相位观测量分别进行误差修正，包括：电离层延迟、相对论效应、天线相位中心误差、地球自转效应；

步骤s3：对修正误差以后的载波相位观测量进行周跳和粗差的实时探测与修复；具体步骤为：

步骤s3具体包括：

步骤s301：使用m-w组合计算整周模糊度互差项nδ

其中，n1,n2分别为f1,f2频率上的整周模糊度，l1，l2分别为f1,f2频率上的载波相位观测值，ρ1，ρ2分别为f1,f2频率上的测码伪距观测值；

步骤s302：计算前i个历元nδ的平均值＜nδ＞i及其均方根误差σi，为了保证方法的实时性，采用以下方式计算：

其中，nδi为第i个历元的整周模糊度互差项；

步骤s303：使用上述公式计算历元i的整周模糊度互差项nδi、历元i+1的整周模糊度互差项nδ(i+1)以及前i-1个历元整周模糊度互差项的均值＜nδ＞i-1以后，通过|nδi-＜nδ＞i-1|≥σ和|nδ(i+1)-nδi|＜1探测周跳和粗差，如果这两个判断条件同时满足，认为在历元i-1和i之间已经发生了周跳；如果有|nδi-＜nδ＞i-1|≥σ和|nδ(i+1)-nδi|≥1同时成立，认为历元i时刻的相位突变是粗差造成的，公式中的σ为根据星载gnss接收机的实际情况设定好的阈值；

步骤s304：如果判断出来有周跳存在，进行如下计算：

接着分别在周跳发生前后的两个历元联立nδ和w两个方程，分别解出周跳发生前后历元在f1,f2两个频率上的整周模糊度n1,n2，进而解出周跳发生前后历元整周模糊度的变化量，也就是周跳，接着对周跳发生之后的连续观测历元中对观测量进行修正；

如果判断出有粗差存在，则将粗差对应的观测量进行剔除。

步骤s4：使用经过处理的载波相位观测值对测码伪距观测量进行平滑；

考虑到实时计算的需求，使用以下公式来进行处理：

初始条件为：

步骤s5：使用地面上注的igs超快速星历计算导航星的位置和钟差信息；

特别需要说明的是，以gps系统为例，igs发布的精密星历有以下几种：事后精密星历以及钟差的精度最高，但是不能实时获得；快速星历及钟差同样不能实时获取；超快速星历分为观测和预报两种，其中观测的同样不能实时获得，预报的可以提供未来24h的星历和钟差信息，适用于本发明。

由于精密星历和钟差相邻两个点的时间间隔为15min，因此需要使用九阶拉格朗日插值法计算导航卫星信号发射时刻的位置和钟差。

步骤s6：将上一历元的定位结果作为初值，使用扩展卡尔曼滤波器计算当前历元的低轨卫星定位结果；其中，如果是初始历元则改为预先设定的初值。

状态向量为三维位置向量、三维速度向量以及接收机钟差。

观测向量为经过载波相位平滑以后的测码伪距观测量。

扩展卡尔曼滤波器的状态方程为低轨卫星的动力学方程，方程中考虑的几种摄动力模型如表1所示:

表1摄动力模型以及参数设置

在计算一步预测值时，使用4阶龙格库塔法对动力学方程进行数值积分。

为了减小计算量，关于位置和速度的状态转移矩阵为：

其中

对于钟差，则认为其在相邻历元间是不变的，因此对应到转移矩阵中的元素为1，并且与速度和位置无关。

滤波器中的观测方程为：

其中，xi,yi,zi分别对应导航卫星在信号发射时刻在ecef坐标系下的坐标；x,y,z为接收机在ecef坐标系下的坐标；cδt是接收机钟差(导航星钟差已经被修正)。

观测矩阵为：

其中

为接收机与gps星之间的距离。

状态噪声协方差矩阵设置为：

其中qy表示位置和速度的噪声协方差矩阵，在这里认为不同轴上的分量之间相互无关，但是同一轴上的位置和速度是相关的。同一轴上的位置和速度的噪声协方差矩阵形式为：

在对观测噪声r阵进行设置时，认为不同类型的观测值和不同gps星之间的观测值是互不相关的，所以r阵表示为：

步骤s7:判断定位任务是否结束，若没有结束，则返回步骤s1继续进行计算。

如果在当前历元计算完成以后还要接着进行定位，则需要等待一个设定好的历元间隔才可以继续计算，这个历元间隔的设定必须保证星载计算机在这段时间里可以完成步骤s1～s7的所有工作，因此需要针对不同的硬件环境作出调整。

以上对本发明的具体实施例进行了描述。需要理解的是，本发明并不局限于上述特定实施方式，本领域技术人员可以在权利要求的范围内做出各种变化或修改，这并不影响本发明的实质内容。在不冲突的情况下，本申请的实施例和实施例中的特征可以任意相互组合。