一种五系统GNSS超快速轨道并行解算方法

技术领域

本发明属于导航定位技术领域，特别涉及一种五系统GNSS超快速轨道并行解算方法。

背景技术

精密的卫星轨道是GNSS高精度应用的基础，尤其对于基于非差技术的精密单点定位(Precise Point Positioning，PPP)应用。而随着用户对更高时效的GNSS高精度服务的需求加大，事后精密星历已经不能满足用户需要，这就推动了超快速实时轨道产品的发展。

GNSS精密轨道的确定本质上都是基于大规模GNSS网的参数估计问题。由于需要使用复杂的算法流程处理大规模的各类观测数据来求解大量的待估参数，同时具有计算密集型和数据密集型特性，因此GNSS高精度产品生成往往都需要经过较长时间的计算求解过程。而随着全球导航卫星系统(Global Navigation Satellite Systems，GNSS)的蓬勃发展与不断建设，进一步给GNSS数据的高效处理带来了巨大挑战，大量的GNSS卫星与海量的多频观测数据引入了数以万计的待估参数，传统的单机集中式数据处理模式或以单纯增加计算机配置的解决方案已经难以满足实时高精度产品生产的需求。

发明内容

本发明的目的是提供一种五系统GNSS超快速轨道并行解算方法，通过五系统的加入来充分发挥GNSS大网的优势，并利用多种现代高性能计算技术对GNSS高精度实时定轨多处功能模块进行充分加速，可以大幅提高GNSS超快轨道产品的生产效率与产品质量，有效地解决大规模数据量与实时快速解算需求之间的矛盾。

本发明采取的技术方案是一种五系统GNSS超快速轨道并行解算方法，包括以下步骤：

S1，观测数据预处理；

S2，轨道初始化；

S3，轨道浮点解计算；

S4，轨道固定解计算；

S5，轨道产品生成。

进一步地，所述S1的具体步骤包括：

S11，接收系统调度命令，定时启动；

S12，接收机钟跳探测与修复；

S13，使用Turboedit方法进行周跳探测与粗差剔除。

进一步地，所述S13中，Turboedit方法主要包含两部分，即分别组成Melbourne-Wübbena组合和GF组合进行周跳探测和粗差剔除。

进一步地，所述S1中，为提高数据预处理的速度，采用OpenMP技术实现多线程并行处理，具体步骤包括：

⑴通过OpenMP按观测数均匀规划出10个线程；

⑵每个线程各自启动钟跳探测和修复程序；

⑶修复后每个线程各自启动周跳探测程序。

进一步地，所述S2的具体步骤包括：

S21，基于广播星历求出不含动力学参数的轨道状态，然后进行一次轨道积分得到基于该轨道状态参数的卫星轨道；

S22，通过轨道动力学拟合求出顾及轨道力学模型的轨道状态参数；

S23，再进行一次轨道积分得到顾及轨道力学模型的状态转移矩阵初值。

进一步地，轨道积分的具体步骤包括：

⑴变量初始化，并读取轨道状态数据及必要产品；

⑵逐卫星进行轨道积分，直到所有卫星处理完毕；

⑶输出完整积分结果；

由于轨道积分是逐颗卫星循环处理，且各卫星的轨道参数相互独立，因此可进行星间并行轨道积分，基于OpenMP技术实现轨道积分并行化，具体步骤包括：

⑴设置线程数；

⑵定义并行域并指定任务调度方式；

⑶轨道积分。

进一步地，所述S3的具体步骤为：得到轨道状态参数迭代初值和状态转移矩阵初值后，通过序贯最小二乘算法进行参数估计，模糊度参数按照实数进行处理；每完成一次最小二乘计算都要基于更新后的轨道状态参数进行轨道积分，以得到下一次迭代所需要的状态转移矩阵；此外，还需根据参数估计残差按照一定的粗差剔除策略对观测数据进行残差编辑；其中，采用的迭代方式为：进行三次最小二乘迭代，其中前两次均进行残差编辑。

进一步地，为避免对高阶矩阵求逆，序贯最小二乘法在法方程叠加的过程中逐历元引入新的活跃参数并消除非活跃参数，最后再根据需要恢复被消除的相关参数。

进一步地，序贯最小二乘前后历元间具有相关性，不适合进行历元间并行，而各测站在进行法方程组建时，彼此之间是相互独立进行的，具有可并行性，因此选择站间并行的方式并基于OpenMP技术实现组建法方程多线程并行优化，具体步骤包括：

⑴声明线程锁对象；

⑵初始化线程锁；

⑶设置线程数；

⑷定义并行域并指定任务调度方式；

⑸站星对(i，j)观测方程组建；

⑹OpenMP上锁操作；

⑺法方程叠加；

⑻OpenMP解锁操作；

⑼销毁线程锁。

进一步地，当参与定轨的测站数目较多时，每个历元参数消除过程本身将变得非常耗时，因此对其进行相应的优化，参数消除算法优化的具体步骤包括：

⑴将元素级循环消参转换为矩阵/向量运算，并使用高性能线性代数库完成；

⑵在步骤⑴的基础上，基于OpenMP对高性能线性代数库的支持，进一步提升参数消除的效率。

进一步地，所述S4的具体步骤包括：

S41，按照基线选取策略依次确定星间单差基线与星地双差基线；

S42，计算并固定所选取的双差模糊度；

S43，将模糊度固定信息代入最小二乘中求出各类参数的固定解。

进一步地，双差模糊度固定数据处理过程的主要耗时部分包括计算各测站MW组合观测值与构建双差模糊度两部分，其中，

计算各测站MW组合观测值的具体步骤包括：

⑴外部逐GNSS系统循环；

⑵内部逐测站计算MW组合观测值；

⑶外部循环结束后退出此处理流程；

构建双差模糊度的具体步骤包括：

⑴外部逐GNSS系统循环；

⑵内部逐基线计算双差模糊度；

⑶外部循环结束后退出此处理流程。

进一步地，进行MW组合观测值的计算与构建双差模糊度时均存在两层循环，即外部的GNSS系统间循环与内部的测站间或基线间循环，这里的数据处理对于不同GNSS系统、不同测站或者不同基线来说都是彼此独立的，因此基于OpenMP技术对以上两个部分的处理进行多线程并行优化，此外，由于本发明为五系统联合解算，所以外层循环只能进行五次，采用系统间大粒度并行方式时的并行线程数至多为五，其中，

基于OpenMP的MW组合观测值计算站间并行优化的具体步骤包括：

⑴GNSS系统间循环；

⑵设置线程数；

⑶定义并行域并指定任务调度方式；

⑷站间循环；

⑸计算MV组合观测值；

基于OpenMP的双差模糊度构建基线间并行优化的具体步骤包括：

⑴GNSS系统间循环；

⑵设置线程数；

⑶定义并行域并指定任务调度方式；

⑷基线间循环；

⑸构建双差模糊度。

进一步地，所述S5的具体步骤为：得到轨道状态参数固定解后，进行一次轨道积分即可得到相应的高精度卫星轨道，同时，还需要进行轨道精度评估，此外，最终向用户发布轨道产品时，需将软件内部格式的卫星轨道数据转换为SP3格式的文件。

进一步地，轨道产品的精度检核可分为外符合和内符合两部分，其中，内符合精度通过重叠弧段比较来实现；外符合精度主要通过SLR检核轨道产品，与其他机构的实时轨道进行精度对比来表现。

进一步地，生成的实时轨道产品在满足质量标准后，可以将轨道产品转换为SSR改正数，并通过播发软件以RTCM-SSR格式转发到Caster上，用户可以通过接收软件接收Caster上的实时产品。

本发明的有益效果在于：

1、本发明实现了一种五系统GNSS超快速轨道并行解算算法，本发明采用OpenMP、高性能线性代数库以及基于CUDA的GPU并行加速等多种现代高性能计算技术对GNSS高精度实时定轨多处功能模块进行充分加速。

2、相比起传统的GNSS高精度轨道产品生产，本发明基于精密定轨流程中各个阶段的耗时情况与算法特点，对钟跳修复、周跳探测、轨道积分、最小二乘计算、残差编辑以及双差模糊度固定都针对性地进行了优化，最终设计出GNSS高精度轨道产品生成效率提升方案。

附图说明

图1为本发明的总体流程图；

图2为本发明采用OpenMP技术实现多线程并行处理的流程图；

图3为本发明轨道积分的流程图；

图4为本发明计算各测站MW组合观测值的流程图；

图5为本发明构建双差模糊度的流程图。

具体实施方式

为了使本发明的目的、技术方案和优点更加清楚、明白，下面将结合附图对本发明的技术方案作进一步详细描述。

本发明所提出的一种五系统GNSS超快速轨道并行解算方式本质上是充分利用现代计算资源，发挥各类现代高性能计算技术的优势，提升超快速轨道的计算效率。本发明采用OpenMP、高性能线性代数库以及基于CUDA的GPU并行加速等现代高性能计算技术针对GNSS高精度实时轨道产品生成效率的提升方法按照一定的技术路线进行了详细研究并给出了相应的最优加速方案设计。

本发明在对原始观测数据进行钟跳探测与修复之后，使用TurboEdit方法进行周跳探测与粗差剔除；随后，基于广播星历通过轨道动力学拟合求出顾及轨道力学模型的轨道状态参数作为最小二乘迭代的初值；然后，经过三次序贯最小二乘迭代对浮点解轨道状态参数进行参数估计更新；最后，计算并固定所选取的双差模糊度，将模糊度固定信息代入最小二乘中求出各类参数的固定解，再进行一次轨道积分即可得到相应的高精度卫星轨道产品。

本发明公开一种五系统GNSS超快速轨道并行解算方法，其中，五系统GNSS为全球卫星导航系统，包括美国的GPS系统、俄罗斯的GLONASS系统、中国的北斗导航卫星系统(BeiDou Navigation Satellite System，BDS)、欧洲的伽利略(Galileo)系统以及日本的QZSS系统，超快速轨道为实时预报的精密轨道，并行解算为多个处理单元来对一个问题进行同步求解，如图1所示，所述方法具体包括以下步骤：

S1，观测数据预处理。

由于测站接收机的钟跳会影响到周跳探测的准确性，因此首先要对原始观测数据进行钟跳探测与修复，之后在此基础上使用TurboEdit方法进行周跳探测与粗差剔除，具体步骤包括：

S11，接收系统调度命令，定时启动。

S12，接收机钟跳探测与修复。

原始观测文件和广播星历中不可避免存在钟跳，钟跳对周跳探测极为不利，它将导致现有的部分周跳探测方法失效，因此在数据处理前需要进行钟跳修复。

S13，使用Turboedit方法进行周跳探测与粗差剔除。

Turboedit方法主要包含两部分，即分别组成Melbourne-Wübbena组合和GF组合进行周跳探测和粗差剔除。

不考虑硬件延迟情况下，Melbourne-Wübbena组合观测值为：

式中，λ

宽巷观测有较长的波长，使粗差的检测更加容易，也使在较短的几个历元之内就能利用递推公式求出宽巷模糊度的估计值和方差，具体公式如公式(2)、公式(3)：

对于当前历元n，若

GF(无几何距离)组合观测值消除历元间观测几何图形的影响，并且消除了接收机钟差、卫星钟差及对流层等所有与频率无关的误差的影响，仅包含电离层影响和整周模糊度项及有频率相关的观测噪声：

L

由于在未发生周跳的情况下，整周模糊度保持不变，且电离层影响变化缓慢，因此，此组合观测值尤为适合粗差的剔除、周跳的探测。在实际应用中，一般利用多项式Q拟合L

|[L

|[L

则认为i历元发生周跳。

进一步地，为了提高数据预处理的速度，本发明采用OpenMP(Open Multi-Processing，共享存储并行编程)技术来实现多线程并行处理(一般设置为10个线程)，如图2所示，具体步骤包括：

⑴通过OpenMP按观测数均匀规划出10个线程。

⑵每个线程各自启动钟跳探测和修复程序。

⑶修复后每个线程各自启动周跳探测程序。

在数据预处理结束后，将生成每个测站的周跳文件，该文件记录了周跳与粗差探测中被剔除的地面站与卫星观测值。

S2，轨道初始化，具体步骤包括：

S21，为了获得精密定轨最小二乘迭代的初值，首先基于广播星历求出不含动力学参数的轨道状态，然后进行一次轨道积分得到基于该轨道状态参数的卫星轨道。

S22，由于此时尚未顾及各种轨道力学模型，积分结果与广播星历之间存在较大差异，因此还需要通过轨道动力学拟合求出顾及轨道力学模型的轨道状态参数，作为最小二乘迭代的初值。

S23，最后再进行一次轨道积分即可得到顾及轨道力学模型的状态转移矩阵初值，用于构建最小二乘观测方程中的设计矩阵。

具体地，如图3所示，轨道积分的具体步骤包括：

⑴变量初始化，并读取轨道状态数据及必要产品。

⑵逐卫星进行轨道积分，直到所有卫星处理完毕。

⑶输出完整积分结果。

由图3可以看出，轨道积分是逐颗卫星循环处理的，由于各卫星的轨道参数相互独立，因此可进行星间并行轨道积分。本发明基于OpenMP技术实现轨道积分并行化，具体步骤包括：

⑴设置线程数。

⑵定义并行域并指定任务调度方式。

⑶轨道积分。

S3，轨道浮点解计算，具体步骤为：得到轨道状态参数迭代初值和状态转移矩阵初值后，即可通过序贯最小二乘算法进行参数估计，模糊度参数按照实数进行处理。每完成一次最小二乘计算都要基于更新后的轨道状态参数进行轨道积分，以得到下一次迭代所需要的状态转移矩阵。此外，还需根据参数估计残差按照一定的粗差剔除策略对观测数据进行必要的质量控制，本发明将该过程称为残差编辑。由于精密定轨观测方程的非线性较强，通常需要进行多次迭代以消除非线性误差影响。本发明所采用的迭代方式为：进行三次最小二乘迭代，其中前两次均进行残差编辑。

进一步地，为避免对高阶矩阵求逆，节省计算资源，提高数据处理效率，序贯最小二乘法在法方程叠加的过程中逐历元引入新的活跃参数并消除非活跃参数，最后再根据需要恢复被消除的相关参数。

具体地，在历元k，可将待估参数向量X

A

式中，P

相应的法方程可以表示为：

在历元k+1，假设出现新的活跃参数X

V

式中，Φ

使用参数消除算法将历元k的非活跃参数X

消除历元k的非活跃参数X

由于序贯最小二乘前后历元间具有相关性，因此不适合进行历元间并行，而各测站在进行法方程组建时，彼此之间是相互独立进行的，具有可并行性。因此，本发明选择站间并行的方式并基于OpenMP技术实现组建法方程多线程并行优化。

在站间并行组建法方程的具体实现中，考虑到OpenMP采用共享内存的方式进行多线程并行，每个历元的法方程矩阵所对应的变量是线程间共享的，因此需要注意在获得单个测站所有观测方程后不能直接进行法方程的叠加操作，否则可能会因为多个线程同时修改法方程矩阵从而导致访问冲突等异常情况。为保证程序的稳健性与正确性，在进行法方程叠加时需利用OpenMP的线程锁机制来确保同一运行时刻最多只有一个线程能够访问法方程矩阵变量。

基于OpenMP的站间并行组建法方程的具体步骤包括：

⑴声明线程锁对象。

⑵初始化线程锁。

⑶设置线程数。

⑷定义并行域并指定任务调度方式。

⑸站星对(i，j)观测方程组建。

⑹OpenMP上锁操作。

⑺法方程叠加。

⑻OpenMP解锁操作。

⑼销毁线程锁。

同时，当参与定轨的测站数目较多时，每个历元参数消除过程本身将变得非常耗时，因此有必要对其进行相应的优化。参数消除算法优化的具体步骤包括：

⑴将元素级循环消参转换为矩阵/向量运算，并使用高性能线性代数库完成。

⑵在步骤⑴的基础上，基于OpenMP对高性能线性代数库的支持，进一步提升参数消除的效率。

将上述的法方程矩阵分别写成如下所示的向量形式可得：

式中，p和q为1×3的行向量；w和v为3×1的列向量；X

因此可以按照上式以法方程系数阵中的行向量为基本处理单元来逐行进行等效消参，并利用高性能线性代数库来实现程序的高效运行。本发明主要采用高性能计算领域中常用的Eigen库和OpenBLAS库来加速矩阵/向量运算，并通过开启Eigen库与OpenBLAS库的OpenMP支持选项，进一步提高参数消除的执行效率。

S4，轨道固定解计算。

完成轨道浮点解计算并获得较高精度的参数估计后，为了进一步提高定轨精度，需将固定模糊度为整数。

轨道固定解计算的具体步骤包括：

S41，按照基线选取策略依次确定星间单差基线与星地双差基线。

S42，计算并固定所选取的双差模糊度。

S43，将模糊度固定信息代入最小二乘中求出各类参数的固定解。

具体地，双差模糊度固定数据处理过程的主要耗时部分包括计算各测站MW组合观测值与构建双差模糊度两部分。

如图4所示，计算各测站MW组合观测值的具体步骤包括：

⑴外部逐GNSS系统循环。

⑵内部逐测站计算MW组合观测值。

⑶外部循环结束后退出此处理流程。

如图5所示，构建双差模糊度的具体步骤包括：

⑴外部逐GNSS系统循环。

⑵内部逐基线计算双差模糊度。

⑶外部循环结束后退出此处理流程。

从图中可以看出，在进行MW组合观测值的计算与构建双差模糊度时均存在两层循环，即外部的GNSS系统间循环与内部的测站间或基线间循环。显然，这里的数据处理对于不同GNSS系统、不同测站或者不同基线来说都是彼此独立的，因此，本发明基于OpenMP技术对以上两个部分的处理进行多线程并行优化。

此外，由于这里是五系统联合解算，所以外层循环只能进行五次，采用系统间大粒度并行方式时的并行线程数至多为五，若使用更高性能的八核处理器，则存在三个CPU核心始终处于闲置状态，就会造成计算资源的严重浪费。

综上所述，对于构建双差模糊度的并行优化，采用基线间小粒度并行的方法可以达到最好的加速效果，能够充分利用现代计算资源，发挥多核处理器的性能。计算MW组合观测值部分的并行优化与之类似，对内层站间循环进行小粒度的多线程并行优化为最佳优化方式。

基于OpenMP的MW组合观测值计算站间并行优化的具体步骤包括：

⑴GNSS系统间循环。

⑵设置线程数。

⑶定义并行域并指定任务调度方式。

⑷站间循环。

⑸计算MV组合观测值。

基于OpenMP的双差模糊度构建基线间并行优化的具体步骤包括：

⑴GNSS系统间循环。

⑵设置线程数。

⑶定义并行域并指定任务调度方式。

⑷基线间循环。

⑸构建双差模糊度。

S5，轨道产品生成。

得到轨道状态参数固定解后，进行一次轨道积分即可得到相应的高精度卫星轨道，同时，还需要进行轨道精度评估，包括与IGS分析中心精密产品进行一致性比较、轨道重复性检验以及激光测卫检核等，此外，最终向用户发布轨道产品时，需将软件内部格式的卫星轨道数据转换为SP3格式的文件。

具体地，轨道产品的精度检核可分为外符合和内符合两部分。其中，内符合精度通过重叠弧段比较来实现；外符合精度主要通过SLR检核轨道产品，与其他机构的实时轨道进行精度对比来表现。

⑴重叠弧段比较：不同时间产成的实时轨道产品存在公共重叠部分，这部分重叠误差可以作为产品精度的内符合标准。

⑵与其他机构实时轨道比较：目前国际上有多个机构和组织可以提供高精度GNSS数据、产品和服务，其中最为主要的就是IGS/MGEX、iGMAS等。将产品与CODE、GFZ等分析中心的实时精密轨道产品进行比较，可以直观得到本系统实时轨道的外符合精度。

⑶SLR轨道检核：卫星激光测距(Satellite Laser Ranging，SLR)是目前唯一能直接给出无模糊度亚厘米级站星距离观测值的空间大地测量技术，一直是检核卫星精密轨道外符合精度的最可靠手段。

生成的实时轨道产品在满足质量标准后，可以将轨道产品转换为SSR改正数，并通过播发软件以RTCM-SSR格式转发到Caster上，用户可以通过接收软件接收Caster上的实时产品。

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