北斗和GPS观测值周跳探测与修复方法、设备及存储设备与流程

本发明涉及卫星导航定位技术领域，尤其涉及北斗和gps观测值周跳探测与修复方法、设备及存储设备。

背景技术：

整周模糊度是gnss(globalnavigationsatellitessystem)载波相位观测值的未知整周部分，只有正确的解算模糊度，才能利用高精度的载波相位观测值进行高精度应用，例如rtk(real-timekinematic)、ppp(precisepointpositioning)等。然而，在实际环境中，gnss接收机在接收卫星信号过程中，卫星信号会受到外部环境的遮挡，从而导致接收机不能连续跟踪信号也称为信号失锁，当重新跟踪到卫星信号时，其载波相关观测值的未知整周部分会发生变化，这种变化即为周跳，在gnss数据处理的过程中，必须对周跳进行处理，否则会引入误差，例如gpsl1频段的载波相位观测值发生1周的周跳，弱不能将其探测到，则会引入0.19米的误差。另外，若对发生的周跳不进行修复，则需要对该卫星进行重新初始化，影响定位的连续性。

北斗和gps静止观测站是基于北斗和gps位置服务中的重要设施之一，例如连续运行参考站网，其是现代化、大众化、集约化、高质量测绘基准体系的重要基础设施之一，也是高精度时空位置服务的重要基础设施之一。参考站网的连续稳定的模糊度解算是实现其高精度时空基准的基础，所以参考站的周跳探测与修复对于整个连续运行参考站网络的有效运行具有十分重要的意义。另外类似的静止观测站应用还包括：北斗和gps变形监测站、北斗和gps气象监测站、北斗和gps时间比对观测站等等。所以研究北斗和gps静止观测值的周跳探测与修复具有重要的理论和应用意义。

目前对实时gnss周跳的探测和修复研究较多，主要可以分为两种：第一种即为分卫星进行单独探测，这种方法主要通过多种观测值或者不同频率间的观测值进行组合，例如blewitt提出的经典turboedit方法，主要就是利用gf组合和mw组合进行周跳探测与修复。这类方法可以对每颗卫星进行探测和修复，流程简单且有效。然而，这类方法需要将多频观测值进行组合探测，或者引入噪声较大的观测值例如伪距、多普勒等，引入其他传感器观测值时会增加硬件成本。

另一种方法则以历元间差分法为代表，顾及卫星的几何结构，利用粗差探测法进行周跳探测，进而采用最小二乘的方法进行周跳修复。这类方法利用了电离层延迟、对流层延迟、硬件延迟等误差的时间相关性，利用历元间差分将其消除，在高采样率(例如1hz)时可以对单频非差非组合数据进行周跳探测与修复，缺点在于探测周跳时需要利用粗差探测的方法，当存在多个周跳时，该方法不易准确探测周跳。

静止观测站多具有采样率高，卫星数多的特点，且ppp、ppp-rtk等技术的兴起也对非差非组合的周跳探测与修复提出了要求，所以使用历元间差分法进行参考站的周跳探测与修复，同时利用参考站处于静止状态的特点，迫切需要构建一种适用于北斗和gps静止观测站的非差非组合周跳修复与探测方法，能够探测观测值中整周周跳和非整周粗差，从而保证静止观测站连续可靠的模糊度解算，实现连续可靠的高精度位置服务。

技术实现要素：

为了解决上述问题，本发明提供了北斗和gps观测值周跳探测与修复方法、设备及存储设备，北斗和gps观测值周跳探测与修复方法，主要包括以下步骤：

请参考图1，图1是本发明实施例中北斗和gps观测值周跳探测与修复方法的流程图，具体包括如下步骤：

s101：获取卫星接收机观测到的北斗和gps卫星的载波相位观测值；根据相邻历元间的载波相位观测值生成相邻历元间的载波相位差分观测值；

s102：按照卫星高度角从高到低选择一颗参考卫星，根据相邻历元间的载波相位差分观测值得到各卫星与参考卫星间的双差观测值；根据双差观测值选择出不存在粗差的卫星观测值；粗差包括整周周跳和非整周的粗差；

s103：剔除不存在粗差的卫星观测值，得到所有存在粗差的卫星观测值，对存在粗差的卫星观测值的模糊度进行设置，即每个存在粗差的卫星观测值的模糊度对应一个模糊度变化参数；

s104：采用阻尼最小二乘法对模糊度变化参数进行求解，得到模糊度变化浮点解和协方差矩阵；模糊度变化浮点解由各卫星观测值的模糊度变化参数对应的浮点解组成；

s105：根据模糊度变化浮点解和协方差矩阵，采用lambda方法对模糊度变化整数解进行搜索，得到模糊度变化整数解；模糊度变化整数解由各卫星观测值的模糊度变化参数对应的整数解组成；

s106：根据模糊度变化整数解，对模糊度变化整数解中未被固定的模糊度参数所对应的卫星观测值进行剔除，将模糊度变化整数解中被固定的模糊度参数对应的卫星观测值作为存在整周周跳的卫星观测值，进行周跳修复。

进一步地，在所述步骤s102的具体步骤中，包括：

s201：根据卫星的高度角，从高到低将所有卫星进行排序；

s202：选择序号为i的卫星为参考卫星，计算得到剩余的每颗卫星与参考卫星的双差观测值；其中第一次选择的参考卫星的序号i＝1；

s203：将各双差观测值的绝对值从小到大进行排序；

s204：若所有的卫星观测值中有n颗不存在粗差，则从小到大依次选择n-1颗卫星对应的双差观测值求均方根误差，均方根误差rmse的计算公式如公式(1)所示：

s205：若均方根误差rmse小于预设值x1，则序号为i的参考卫星和参与均方根误差计算的n-1颗卫星为不存在粗差的n颗卫星，到步骤s208；

s206：若均方根误差rmse大于等于预设值x1，则序号为i的参考卫星为存在粗差的卫星，将序号为i的参考卫星剔除，并判断所有的卫星是否都已经遍历；若是，则到步骤s208，否则，到步骤s207；

s207：将i更新为i+1，返回步骤s202；

s208：结束搜索。

进一步地，在步骤s205中，预设值x1的值为1cm。

进一步地，在步骤s104中，采用阻尼最小二乘法对存在粗差的卫星观测值的模糊度变化参数求解的具体方法为：

对存在粗差的卫星观测值设置模糊度变化参数，如公式(2)所示：

公式(2)中，为载波相位差分观测值；δε为接收机噪声和极小的未消除的误差；e^x，e^y和e^z为三个分量的方向余弦；和分别为相邻历元间的各位置变化分量；δt为相邻历元间的接收机钟差变化值；

根据公式(2)，构建阻尼最小二乘法方程，如公式(3)所示：

其中，为待解算的未知量，代表钟差变化加模糊度变化参数；

代表位置参数；p代表观测值权阵，p的取值根据双差观测值的噪声取值确定；a和b为系数矩阵，根据公式获取；代表载波相位差分观测值矩阵；px为阻尼因子，代表先验坐标的权阵，px的计算如公式(4)：

其中，(σxσyσz)为相邻历元间标准差，为相邻历元间载波相位差分观测值的单位权方差；

求解出公式(3)，即可得到模糊度变化浮点解和协方差矩阵，模糊度变化浮点解为t_n中剔除钟差变化cδt后剩下的部分，协方差矩阵为a^tpa+px的逆。

进一步地，相邻历元间标准差的值为(111)，相邻历元间载波相位差分观测值的单位权方差的值为0.0001m。

进一步地，在步骤s105中，采用lambda方法对模糊度变化整数解进行搜索的具体步骤为：

s301：根据步骤s104中计算得到的协方差矩阵，求得模糊度变化整数解的变化范围和精度，以步骤s104中求出的模糊度变化浮点解为初始解，结合模糊度变化整数解的范围和精度，随机生成若干组模糊度变化整数解，每组模糊度变化整数解对应的各模糊度参数与模糊度变化浮点解中的各模糊度参数对应，将这若干组模糊度变化整数解所组成的组合作为模糊度变化备选组；

s302：根据公式(5)搜索模糊度变化备选组中的最优模糊度变化整数解和次优模糊度变化整数解：

其中，z为降相关矩阵，￠为整数集，为协方差矩阵；和z分别为模糊度变化浮点解和模糊度变化整数解a通过降相关矩阵z转换后的对应量；对计算得到的所对应的标量值从小到大进行排序，最小的标量值对应的模糊度变化整数解为最优模糊度变化整数解，次最小的标量值对应的模糊度变化整数解为次优模糊度变化整数解；

s303：根据最优模糊度变化整数解和次优模糊度变化整数解，计算ratio值，ratio值为次优模糊度变化整数解和最优模糊度变化整数解的残差二次型的比值，计算公式如公式(6)所示:

公式(6)中，为最优模糊度变化整数解，为次优模糊度变化整数解；

s304：根据ratio值，通过阈值t判断是否接受最优模糊度变化整数解，阈值t由载波相位差分观测值和阻尼最小二乘法求取的模糊度变化浮点解精度确定；

当ratio≥t时，接受最优模糊度变化整数解，将最优模糊度变化整数解中所有的模糊度参数均标记为被固定的模糊度参数，到步骤s307；

当ratio<t时，不接受最优模糊度变化整数解，到步骤s305；

s305：采用部分模糊度固定法，剔除最优模糊度变化整数解和次优模糊度变化整数解中不同的模糊度参数，保留最优模糊度变化整数解和次优模糊度变化整数解中相同的模糊度参数；

判断剔除的模糊度参数对应的模糊度变化浮点解参数个数是否大于预设值x2；若是，则到步骤s306；若否，则将最优模糊度变化整数解和次优模糊度变化整数解中不同的模糊度参数标记为未被固定的模糊度参数，将最优模糊度变化整数解和次优模糊度变化整数解中相同的模糊度参数标记为被固定的模糊度参数，到步骤s307；

s306：对被剔除的模糊度参数再次进行搜索，并将新搜索出来的模糊度参数更新到最优模糊度变化整数解和次优模糊度变化整数解中；到步骤s303；

s307：结束搜索，最终得到的最优模糊度变化整数解即为搜索到的模糊度变化整数解。

进一步地，预设值x2的值为5。

一种存储设备，所述存储设备存储指令及数据用于实现北斗和gps观测值周跳探测与修复方法。

一种北斗和gps观测值周跳探测与修复设备，包括：处理器及所述存储设备；所述处理器加载并执行所述存储设备中的指令及数据用于实现北斗和gps观测值周跳探测与修复方法。

本发明提供的技术方案带来的有益效果是：本发明所提供的技术方案，可以满足北斗和gps静止观测站非差非组合的周跳探测与修复，可单独对参考站各个频率的载波相位观测值进行处理，在高采样率下可以实现稳定可靠的周跳探测与修复；在周跳探测阶段仅需选择出部分未发生周跳的卫星，不需根据残差来进行周跳探测，提高了周跳探测的可靠性，同时简化了探测流程；在周跳修复阶段引入了阻尼最小二乘法和部分模糊度固定方法，可以有效的利用静止观测值坐标变化缓慢的特点，同时可以在修复的过程中准确的选择整周周跳进行修复，剔除非整周粗差的影响。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1是本发明实施例中北斗和gps观测值周跳探测与修复方法的流程图；

图2是本发明实施例中利用双差观测值搜索无周跳或粗差的卫星的流程图示意图；

图3是本发明实施例中基于lambda方法搜索整周周跳及非整周粗差的流程图示意图；

图4是本发明实施例中硬件设备工作的示意图。

具体实施方式

为了对本发明的技术特征、目的和效果有更加清楚的理解，现对照附图详细说明本发明的具体实施方式。

本发明的实施例提供了北斗和gps观测值周跳探测与修复方法、设备及存储设备。

请参考图1，图1是本发明实施例中北斗和gps观测值周跳探测与修复方法的流程图，具体包括如下步骤：

s101：获取卫星接收机观测到的北斗和gps卫星的载波相位观测值；根据相邻历元间的载波相位观测值生成相邻历元间的载波相位差分观测值；

获取北斗和gps原始观测数据，这里原始观测数据既可以是后期解算的rinex文件，也可以是实时数据流进行解码获得的数据，原始观测数据主要包括载波相位观测值和伪距观测值等，以gps的l1频率为例，载波相位观测值计算公式如公式(1)所示：

其中:j代表gps卫星序号，为载波相位观测值，λ是载波的波长，ρ是卫星和测站间的几何距离，c是真空中的光速，tg为接收机钟差，t^j是卫星钟差，i和t分别为电离层和对流层斜路径延迟，n为整周模糊度(即未知的整周计数)，ε包括了未模型化的误差和噪声。通过数据读取或者解码程序，将北斗和gps载波相位观测值读取出来，利用相邻历元的载波相位观测值差分构建相邻历元的载波相位差分观测值计算公式如公式(2)所示：

δ代表差分运算，例如表示历元间的几何距离的差分值，δn^j为整周模糊度变化值，由于对流层误差、电离层误差、轨道误差、卫星钟差、硬件延迟等误差具有很强的时间相关性，即在短时间(例如1s)内的变化很小，那么历元差分可以很大程度上消除这些误差的影响。根据上述分析，历元间的载波相位差分观测值计算公式可以统一写为公式(3)：

其中δε主要为接收机噪声和极小的未消除的误差。

将公式(3)线性化可得载波相位差分观测值的计算公式(4)：

其中e^x，e^y和e^z为三个分量的方向余弦，而和分别为相邻历元的位置变化分量，δt为相邻历元的接收机的钟差变化，δn为整周模糊度变化。本发明没有引入伪距观测值，从而可以避免伪距噪声的影响。

s102：按照卫星高度角从高到低选择一颗参考卫星，根据相邻历元间的载波相位差分观测值得到各卫星与参考卫星间的双差观测值；根据双差观测值选择出不存在粗差的卫星观测值；粗差包括整周周跳和非整周的粗差：选择不存在粗差的卫星观测值的方法如下：

s201：根据卫星的高度角，从高到低将所有卫星进行排序；

s202：选择序号为i的卫星为参考卫星，计算得到剩余的每颗卫星与参考卫星的双差观测值；其中第一次选择的参考卫星的序号i＝1；

s203：将各双差观测值的绝对值从小到大进行排序；

s204：若所有的卫星观测值中有n颗不存在粗差，则从小到大依次选择n-1颗卫星对应的双差观测值求均方根误差，均方根误差rmse的计算公式如公式(5)所示：

s205：若均方根误差rmse小于预设值x1，则序号为i的参考卫星和参与均方根误差计算的n-1颗卫星为不存在粗差的n颗卫星，到步骤s208；

s206：若均方根误差rmse大于等于预设值x1，则序号为i的参考卫星为存在粗差的卫星，将序号为i的参考卫星剔除，并判断所有的卫星是否都已经遍历；若是，则到步骤s208，否则，到步骤s207；

s207：将i更新为i+1，返回步骤s202；

s208：结束搜索。

其中，步骤s205中的预设值x1的值推荐取1，请参考图2，图2是本发明以gps的l1频率为例利用双差观测值搜索不存在周跳的卫星观测值的流程图，由于需要进行周跳探测与修复的卫星观测值的参考站均建立在稳固的平台上，所以参考站短时间的位置变化可以认为为0，即和均为0，将该约束代入公式(4)中可得公式(6)：

根据卫星高度角从大到小，依次选择一颗参考卫星i计算相邻历元间的各卫星与参考卫星间的双差观测值，如公式(7)所示：

公式(7)中，i为参考卫星序号，j为与参考卫星构建双差观测值对应的卫星的序号，若参考卫星没有周跳，则有δn^i,j＝δn^j，代入公式(7)可得公式(8)：

其中，δε^i,j的噪声强度为毫米级，若不存在周跳，则双差观测值残差应在毫米级。在本发明的实施例中，假设至少存在四颗(具体数量应根据具体应用情况确定)没有周跳的卫星，按照卫星高度角从大到小依次选择参考卫星，构建剩余各卫星与参考卫星间的双差观测值，选择双差观测值绝对值最小的三个双差观测值按照式计算均方根误差，δ▽li表示双差观测值，当均方根误差小于预设值x1＝1cm时，即说明参考卫星和构成三个双差观测值的三颗卫星没有周跳，则选择出了不存在周跳的四颗卫星；否则，在所有卫星均遍历后，结束搜索。

s103：剔除不存在粗差的卫星观测值，得到所有存在粗差的卫星观测值，对存在粗差的卫星观测值的模糊度进行设置，即每个存在粗差的卫星观测值的模糊度对应一个模糊度变化参数；

s104：采用阻尼最小二乘法对模糊度变化参数进行求解，得到模糊度变化浮点解和协方差矩阵；模糊度变化浮点解由各卫星观测值的模糊度变化参数对应的浮点解组成；具体求解过程如下：

对存在周跳或粗差的卫星观测值设置模糊度变化参数如公式(9)：

公式(9)中，为载波相位差分观测值；δε为接收机噪声和极小的未消除的误差；e^x，e^y和e^z为三个分量的方向余弦；和分别为相邻历元的位置变化分量；δt为相邻历元的接收机钟差变化；

对公式(9)，直接构建阻尼最小二乘法方程如公式(10)：

根据公式(10)，模糊度变化未知参数的浮点解为t_n中剔除钟差变化cδt后剩下的部分，协方差矩阵为a^tpa+px的逆，解公式(10)即可求得模糊度变化浮点解和协方差矩阵；公式(10)中，为待解算的未知量，代表钟差变化加模糊度变化参数；代表位置参数；p代表观测值权阵，p的取值根据双差观测值的噪声取值确定；a和b为系数矩阵，根据公式获取；代表载波相位差分观测值矩阵；px为阻尼因子，代表先验坐标的权阵，px的计算如公式(11)：

以gps的l1频率为例，公式(11)中，历元间标准差(σxσyσz)均为1cm，历元间载波相位差分观测值的单位权方差为0.0001m。

s105：根据模糊度变化浮点解和协方差矩阵，采用lambda方法对模糊度变化整数解进行搜索，得到模糊度变化整数解；模糊度变化整数解由各卫星观测值的模糊度变化参数对应的整数解组成；请参考图3，图3是本发明实施例中以gps的l1频率为例，基于lambda方法对模糊度变化整数解进行搜索的示意图。具体搜索步骤如下：

s301：根据步骤s104中得到的协方差矩阵，求得模糊度变化整数解的范围和精度，以步骤s104中求出的模糊度变化浮点解为初始解，结合模糊度变化整数解的范围和精度，随机生成若干组模糊度变化整数解，每组模糊度变化整数解对应的各模糊度参数与模糊度变化浮点解中的各模糊度参数对应，将这若干组模糊度变化整数解所组成的组合称为模糊度变化备选组；

s302：根据公式(12)搜索模糊度变化备选组中的最优模糊度变化整数解和次优模糊度变化整数解，公式(12)如下所示：

公式(12)中z为降相关矩阵，￠为整数集；和z为模糊度变化浮点解和模糊度变化整数解a通过降相关矩阵z转换后的对应量，转换的目的是为了加快搜索速度；根据公式(12)中计算的所对应的标量值从小到大进行排序，最小的标量值对应的模糊度变化整数解为最优模糊度变化整数解，次最小的标量值对应的模糊度变化整数解为次优模糊度变化整数解；

s303：计算ratio值，ratio值由次优模糊度变化整数解和最优模糊度变化整数解残差二次型的比值确定，如公式(13):

公式(13)中，为最优模糊度变化整数解，为次优模糊度变化整数解；

s304：根据ratio值，通过阈值t判断是否接受最优模糊度变化整数解；

当ratio≥t时，接受最优模糊度变化整数解，将最优模糊度变化整数解中所有的模糊度参数均标记为被固定的模糊度参数，到步骤s307；

当ratio<t时，不接受最优模糊度变化整数解，到步骤s305；

其中，阈值t为经验值，由载波相位差分观测值和阻尼最小二乘法求取的模糊度变化浮点解的精度确定；

s305：采用部分模糊度固定法，剔除最优模糊度变化整数解和次优模糊度变化整数解中不同的模糊度参数，保留最优模糊度变化整数解和次优模糊度变化整数解中相同的模糊度参数；

判断剔除的模糊度参数对应的模糊度浮点解参数个数是否大于预设值x2；若是，则到步骤s306；若否，则将最优模糊度变化整数解和次优模糊度变化整数解中不同的模糊度参数标记为未被固定的模糊度参数，将最优模糊度变化整数解和次优模糊度变化整数解中相同的模糊度参数标记为被固定的模糊度参数，到步骤s307；

s306：对被剔除的模糊度参数再次进行搜索，并将新搜索出来的模糊度参数更新到最优模糊度变化整数解和次优模糊度变化整数解中；转到步骤s303；

s307：结束搜索，最终得到的最优模糊度变化整数解即为搜索到的模糊度变化整数解。

其中，步骤s305中的预设值x2的值推荐取5。

s106：根据模糊度变化整数解，对模糊度变化整数解中未被固定的模糊度参数所对应的卫星观测值进行剔除，将模糊度变化整数解中被固定的模糊度参数对应的卫星观测值作为存在整周周跳的卫星观测值，进行周跳修复。

请参见图4，图4是本发明实施例的硬件设备工作示意图，所述硬件设备具体包括：北斗和gps观测值周跳探测与修复设备401、处理器402及存储设备403。

北斗和gps观测值周跳探测与修复设备401：所述北斗和gps观测值周跳探测与修复设备401实现所述北斗和gps观测值周跳探测与修复方法。

处理器402：所述处理器402加载并执行所述存储设备403中的指令及数据用于实现所述北斗和gps观测值周跳探测与修复方法。

存储设备403：所述存储设备403存储指令及数据；所述存储设备403用于实现所述北斗和gps观测值周跳探测与修复方法。

本发明的有益效果是：本发明所提供的技术方案，可以满足北斗和gps静止观测站非差非组合的周跳探测与修复，可单独对参考站各个频率的载波相位观测值进行处理，在高采样率下可以实现稳定可靠的周跳探测与修复；在周跳探测阶段仅需选择出部分不存在周跳的卫星，不需根据残差来进行周跳探测，提高了周跳探测的可靠性，同时简化了探测流程；在周跳修复阶段引入了阻尼最小二乘法和部分模糊度固定方法，可以有效的利用静止观测值坐标变化缓慢的特点，同时可以在修复的过程中准确的选择整周周跳进行修复，剔除非整周粗差的影响。

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