CN109143298A - 北斗和gps观测值周跳探测与修复方法、设备及存储设备 - Google Patents
北斗和gps观测值周跳探测与修复方法、设备及存储设备 Download PDFInfo
- Publication number
- CN109143298A CN109143298A CN201810937216.9A CN201810937216A CN109143298A CN 109143298 A CN109143298 A CN 109143298A CN 201810937216 A CN201810937216 A CN 201810937216A CN 109143298 A CN109143298 A CN 109143298A
- Authority
- CN
- China
- Prior art keywords
- ambiguity
- satellite
- change
- value
- ambiguity change
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 75
- 238000001514 detection method Methods 0.000 title claims abstract description 57
- 230000008439 repair process Effects 0.000 title claims abstract description 20
- 230000008859 change Effects 0.000 claims abstract description 185
- 239000011159 matrix material Substances 0.000 claims description 29
- 238000007667 floating Methods 0.000 claims description 27
- 238000004364 calculation method Methods 0.000 claims description 13
- 238000012897 Levenberg–Marquardt algorithm Methods 0.000 claims description 10
- 238000005067 remediation Methods 0.000 claims description 10
- 238000013016 damping Methods 0.000 claims description 6
- 238000012163 sequencing technique Methods 0.000 claims description 6
- DMBHHRLKUKUOEG-UHFFFAOYSA-N diphenylamine Chemical compound C=1C=CC=CC=1NC1=CC=CC=C1 DMBHHRLKUKUOEG-UHFFFAOYSA-N 0.000 claims description 3
- FKZLAQOTFPPENP-UHFFFAOYSA-N n-(2,3-dihydroxypropyl)-n-(2-hydroxypropyl)nitrous amide Chemical compound CC(O)CN(N=O)CC(O)CO FKZLAQOTFPPENP-UHFFFAOYSA-N 0.000 claims description 3
- 230000003247 decreasing effect Effects 0.000 claims description 2
- 230000008569 process Effects 0.000 abstract description 11
- 230000009286 beneficial effect Effects 0.000 abstract description 3
- 238000012544 monitoring process Methods 0.000 abstract description 3
- 230000008901 benefit Effects 0.000 abstract 1
- 230000007246 mechanism Effects 0.000 abstract 1
- 238000005457 optimization Methods 0.000 abstract 1
- 230000003068 static effect Effects 0.000 description 10
- 238000010586 diagram Methods 0.000 description 5
- 238000005070 sampling Methods 0.000 description 4
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 3
- 238000006243 chemical reaction Methods 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000004458 analytical method Methods 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000001934 delay Effects 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000013507 mapping Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 239000011541 reaction mixture Substances 0.000 description 1
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
- G01S19/38—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
- G01S19/39—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
- G01S19/42—Determining position
- G01S19/43—Determining position using carrier phase measurements, e.g. kinematic positioning; using long or short baseline interferometry
- G01S19/44—Carrier phase ambiguity resolution; Floating ambiguity; LAMBDA [Least-squares AMBiguity Decorrelation Adjustment] method
Landscapes
- Engineering & Computer Science (AREA)
- Radar, Positioning & Navigation (AREA)
- Remote Sensing (AREA)
- Computer Networks & Wireless Communication (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Position Fixing By Use Of Radio Waves (AREA)
Abstract
本发明提供了北斗和GPS观测值周跳探测与修复方法、设备及存储设备，北斗和GPS观测值周跳探测与修复方法仅使用相邻历元的非差非组合观测值进行周跳探测与修复，可实时运行，适用于网络RTK、单基站RTK、变形监测等观测站位置变化缓慢的应用场景，且对接收机的频率数没有要求。本发明的有益效果是：本发明所提供的技术方案有效的利用了多系统下卫星数量增加、卫星几何机构优化的优势，适应了当前北斗和GPS从站间双差组合观测处理向非差非组合发展的趋势，简化了周跳探测的流程，提高了探测的成功率和准确率，在周跳修复的过程中使用了部分模糊度固定的方法对非整周粗差进行剔除，可有效的提高北斗和GPS的定位精度。
Description
技术领域
本发明涉及卫星导航定位技术领域，尤其涉及北斗和GPS观测值周跳探测与修复方法、设备及存储设备。
背景技术
整周模糊度是GNSS(Global Navigation Satellites System)载波相位观测值的未知整周部分，只有正确的解算模糊度，才能利用高精度的载波相位观测值进行高精度应用，例如RTK(Real-time Kinematic)、PPP(Precise Point Positioning)等。然而，在实际环境中，GNSS接收机在接收卫星信号过程中，卫星信号会受到外部环境的遮挡，从而导致接收机不能连续跟踪信号也称为信号失锁，当重新跟踪到卫星信号时，其载波相关观测值的未知整周部分会发生变化，这种变化即为周跳，在GNSS数据处理的过程中，必须对周跳进行处理，否则会引入误差，例如GPS L1频段的载波相位观测值发生1周的周跳，弱不能将其探测到，则会引入0.19米的误差。另外，若对发生的周跳不进行修复，则需要对该卫星进行重新初始化，影响定位的连续性。
北斗和GPS静止观测站是基于北斗和GPS位置服务中的重要设施之一，例如连续运行参考站网，其是现代化、大众化、集约化、高质量测绘基准体系的重要基础设施之一，也是高精度时空位置服务的重要基础设施之一。参考站网的连续稳定的模糊度解算是实现其高精度时空基准的基础，所以参考站的周跳探测与修复对于整个连续运行参考站网络的有效运行具有十分重要的意义。另外类似的静止观测站应用还包括：北斗和GPS变形监测站、北斗和GPS气象监测站、北斗和GPS时间比对观测站等等。所以研究北斗和GPS静止观测值的周跳探测与修复具有重要的理论和应用意义。
目前对实时GNSS周跳的探测和修复研究较多，主要可以分为两种：第一种即为分卫星进行单独探测，这种方法主要通过多种观测值或者不同频率间的观测值进行组合，例如Blewitt提出的经典Turbo Edit方法，主要就是利用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)中，为载波相位差分观测值；Δε为接收机噪声和极小的未消除的误差；ex，ey和ez为三个分量的方向余弦；和分别为相邻历元间的各位置变化分量；Δt为相邻历元间的接收机钟差变化值；
根据公式(2)，构建阻尼最小二乘法方程，如公式(3)所示：
其中，为待解算的未知量，代表钟差变化加模糊度变化参数；
代表位置参数；P代表观测值权阵，P的取值根据双差观测值的噪声取值确定；A和B为系数矩阵，根据公式获取；代表载波相位差分观测值矩阵；PX为阻尼因子，代表先验坐标的权阵，PX的计算如公式(4)：
其中，(σx σy σz)为相邻历元间标准差，为相邻历元间载波相位差分观测值的单位权方差；
求解出公式(3)，即可得到模糊度变化浮点解和协方差矩阵，模糊度变化浮点解为t_N中剔除钟差变化cΔt后剩下的部分，协方差矩阵为ATPA+PX的逆。
进一步地，相邻历元间标准差的值为(1 1 1)，相邻历元间载波相位差分观测值的单位权方差的值为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为接收机钟差，tj是卫星钟差，I和T分别为电离层和对流层斜路径延迟，N为整周模糊度(即未知的整周计数)，ε包括了未模型化的误差和噪声。通过数据读取或者解码程序，将北斗和GPS载波相位观测值读取出来，利用相邻历元的载波相位观测值差分构建相邻历元的载波相位差分观测值计算公式如公式(2)所示：
Δ代表差分运算，例如表示历元间的几何距离的差分值，ΔNj为整周模糊度变化值，由于对流层误差、电离层误差、轨道误差、卫星钟差、硬件延迟等误差具有很强的时间相关性，即在短时间(例如1s)内的变化很小，那么历元差分可以很大程度上消除这些误差的影响。根据上述分析，历元间的载波相位差分观测值计算公式可以统一写为公式(3)：
其中Δε主要为接收机噪声和极小的未消除的误差。
将公式(3)线性化可得载波相位差分观测值的计算公式(4)：
其中ex，ey和ez为三个分量的方向余弦，而和分别为相邻历元的位置变化分量，Δ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为与参考卫星构建双差观测值对应的卫星的序号，若参考卫星没有周跳，则有ΔNi,j＝ΔNj，代入公式(7)可得公式(8)：
其中，Δεi,j的噪声强度为毫米级，若不存在周跳，则双差观测值残差应在毫米级。在本发明的实施例中，假设至少存在四颗(具体数量应根据具体应用情况确定)没有周跳的卫星，按照卫星高度角从大到小依次选择参考卫星，构建剩余各卫星与参考卫星间的双差观测值，选择双差观测值绝对值最小的三个双差观测值按照式计算均方根误差，Δ▽Li表示双差观测值，当均方根误差小于预设值X1＝1cm时，即说明参考卫星和构成三个双差观测值的三颗卫星没有周跳，则选择出了不存在周跳的四颗卫星；否则，在所有卫星均遍历后，结束搜索。
S103：剔除不存在粗差的卫星观测值，得到所有存在粗差的卫星观测值，对存在粗差的卫星观测值的模糊度进行设置，即每个存在粗差的卫星观测值的模糊度对应一个模糊度变化参数；
S104：采用阻尼最小二乘法对模糊度变化参数进行求解，得到模糊度变化浮点解和协方差矩阵；模糊度变化浮点解由各卫星观测值的模糊度变化参数对应的浮点解组成；具体求解过程如下：
对存在周跳或粗差的卫星观测值设置模糊度变化参数如公式(9)：
公式(9)中，为载波相位差分观测值；Δε为接收机噪声和极小的未消除的误差；ex，ey和ez为三个分量的方向余弦；和分别为相邻历元的位置变化分量；Δt为相邻历元的接收机钟差变化；
对公式(9)，直接构建阻尼最小二乘法方程如公式(10)：
根据公式(10)，模糊度变化未知参数的浮点解为t_N中剔除钟差变化cΔt后剩下的部分，协方差矩阵为ATPA+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静止观测站非差非组合的周跳探测与修复，可单独对参考站各个频率的载波相位观测值进行处理，在高采样率下可以实现稳定可靠的周跳探测与修复；在周跳探测阶段仅需选择出部分不存在周跳的卫星，不需根据残差来进行周跳探测，提高了周跳探测的可靠性，同时简化了探测流程；在周跳修复阶段引入了阻尼最小二乘法和部分模糊度固定方法，可以有效的利用静止观测值坐标变化缓慢的特点，同时可以在修复的过程中准确的选择整周周跳进行修复，剔除非整周粗差的影响。
以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
Claims (9)
1.北斗和GPS观测值周跳探测与修复方法，其特征在于：包括以下步骤：
S101：获取卫星接收机观测到的北斗和GPS卫星的载波相位观测值；根据相邻历元间的载波相位观测值生成相邻历元间的载波相位差分观测值；
S102：按照卫星高度角从高到低选择一颗参考卫星，根据相邻历元间的载波相位差分观测值得到各卫星与参考卫星间的双差观测值；根据双差观测值选择出不存在粗差的卫星观测值；粗差包括整周周跳和非整周的粗差；
S103：剔除不存在粗差的卫星观测值，得到所有存在粗差的卫星观测值，对存在粗差的卫星观测值的模糊度进行设置，即每个存在粗差的卫星观测值的模糊度对应一个模糊度变化参数；
S104：采用阻尼最小二乘法对模糊度变化参数进行求解，得到模糊度变化浮点解和协方差矩阵；模糊度变化浮点解由各卫星观测值的模糊度变化参数对应的浮点解组成；
S105：根据模糊度变化浮点解和协方差矩阵，采用LAMBDA方法对模糊度变化整数解进行搜索，得到模糊度变化整数解；模糊度变化整数解由各卫星观测值的模糊度变化参数对应的整数解组成；
S106：根据模糊度变化整数解，对模糊度变化整数解中未被固定的模糊度参数所对应的卫星观测值进行剔除，将模糊度变化整数解中被固定的模糊度参数对应的卫星观测值作为存在整周周跳的卫星观测值，进行周跳修复。
2.如权利要求1所述的北斗和GPS观测值周跳探测与修复方法，其特征在于：所述步骤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：结束搜索。
3.如权利要求2所述的北斗和GPS观测值周跳探测与修复方法，其特征在于：在步骤S205中，预设值X1的值为1cm。
4.如权利要求1所述的北斗和GPS观测值周跳探测与修复方法，其特征在于：在步骤S104中，采用阻尼最小二乘法对存在粗差的卫星观测值的模糊度变化参数求解的具体方法为：
对存在粗差的卫星观测值设置模糊度变化参数，如公式(2)所示：
公式(2)中，为载波相位差分观测值；Δε为接收机噪声和极小的未消除的误差；ex，ey和ez为三个分量的方向余弦；和分别为相邻历元间的各位置变化分量；Δt为相邻历元间的接收机钟差变化值；
根据公式(2)，构建阻尼最小二乘法方程，如公式(3)所示：
其中，为待解算的未知量，代表钟差变化加模糊度变化参数；代表位置参数；P代表观测值权阵，P的取值根据双差观测值的噪声取值确定；A和B为系数矩阵，根据公式获取；代表载波相位差分观测值矩阵；PX为阻尼因子，代表先验坐标的权阵，PX的计算如公式(4)：
其中，(σx σy σz)为相邻历元间标准差，为相邻历元间载波相位差分观测值的单位权方差；
求解出公式(3)，即可得到模糊度变化浮点解和协方差矩阵，模糊度变化浮点解为t_N中剔除钟差变化cΔt后剩下的部分，协方差矩阵为ATPA+PX的逆。
5.如权利要求4所述的北斗和GPS观测值周跳探测与修复方法，其特征在于：相邻历元间标准差的值为(1 1 1)，相邻历元间载波相位差分观测值的单位权方差的值为0.0001m。
6.如权利要求1所述的北斗和GPS观测值周跳探测与修复方法，其特征在于：在步骤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：结束搜索，最终得到的最优模糊度变化整数解即为搜索到的模糊度变化整数解。
7.如权利要求6所述的北斗和GPS观测值周跳探测与修复方法，其特征在于：在步骤S305中，预设值X2的值为5。
8.一种存储设备，其特征在于：所述存储设备存储指令及数据用于实现权利要求1～7所述的任意一种北斗和GPS观测值周跳探测与修复方法。
9.北斗和GPS观测值周跳探测与修复设备，其特征在于：包括：处理器及权利要求8所述的存储设备；所述处理器加载并执行所述存储设备中的指令及数据用于实现权利要求1～7所述的任意一种北斗和GPS观测值周跳探测与修复方法。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810937216.9A CN109143298B (zh) | 2018-08-16 | 2018-08-16 | 北斗和gps观测值周跳探测与修复方法、设备及存储设备 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201810937216.9A CN109143298B (zh) | 2018-08-16 | 2018-08-16 | 北斗和gps观测值周跳探测与修复方法、设备及存储设备 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN109143298A true CN109143298A (zh) | 2019-01-04 |
| CN109143298B CN109143298B (zh) | 2020-08-07 |
Family
ID=64789941
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201810937216.9A Active CN109143298B (zh) | 2018-08-16 | 2018-08-16 | 北斗和gps观测值周跳探测与修复方法、设备及存储设备 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN109143298B (zh) |
Cited By (8)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109917356A (zh) * | 2019-03-13 | 2019-06-21 | 武汉际上导航科技有限公司 | 一种机载激光扫描系统误差标定方法 |
| CN110398758A (zh) * | 2019-07-24 | 2019-11-01 | 广州中海达卫星导航技术股份有限公司 | 实时钟差估计中的粗差探测方法、装置、设备及存储介质 |
| CN111045062A (zh) * | 2019-11-29 | 2020-04-21 | 航天恒星科技有限公司 | 一种基于电磁星的星基电离层反演方法 |
| CN112083464A (zh) * | 2019-06-14 | 2020-12-15 | 北京合众思壮科技股份有限公司 | 一种部分模糊度的固定方法及装置 |
| CN112381309A (zh) * | 2020-11-23 | 2021-02-19 | 珠江水利委员会珠江水利科学研究院 | 水库大坝安全监测预警方法、装置、系统及存储介质 |
| CN112526573A (zh) * | 2021-02-07 | 2021-03-19 | 腾讯科技（深圳）有限公司 | 对象定位方法和装置、存储介质及电子设备 |
| CN113093236A (zh) * | 2019-12-23 | 2021-07-09 | 中国石油天然气集团有限公司 | 动态后处理的方法和装置 |
| CN114594502A (zh) * | 2022-02-18 | 2022-06-07 | 真点科技(北京)有限公司 | 一种实现周跳处理的方法、装置、计算机存储介质及终端 |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JPH09119972A (ja) * | 1995-10-26 | 1997-05-06 | Furuno Electric Co Ltd | 相対測位装置および相対測位方法 |
| CN102288978A (zh) * | 2011-07-20 | 2011-12-21 | 东南大学 | 一种cors基站周跳探测与修复方法 |
| CN103529462A (zh) * | 2013-10-21 | 2014-01-22 | 西南交通大学 | 一种用于全球导航卫星系统的动态周跳探测与修复方法 |
| CN104101890A (zh) * | 2014-07-01 | 2014-10-15 | 同济大学 | 一种对流层长距离亚厘米级实时动态卫星导航定位方法 |
| CN105549056A (zh) * | 2015-12-03 | 2016-05-04 | 中国电子科技集团公司第二十研究所 | 一种相对定位装置及其载波整周模糊度解算方法 |
| CN107102346A (zh) * | 2017-06-08 | 2017-08-29 | 中国电子科技集团公司第五十四研究所 | 一种基于北斗系统的多天线测姿方法 |
-
2018
- 2018-08-16 CN CN201810937216.9A patent/CN109143298B/zh active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| JPH09119972A (ja) * | 1995-10-26 | 1997-05-06 | Furuno Electric Co Ltd | 相対測位装置および相対測位方法 |
| CN102288978A (zh) * | 2011-07-20 | 2011-12-21 | 东南大学 | 一种cors基站周跳探测与修复方法 |
| CN103529462A (zh) * | 2013-10-21 | 2014-01-22 | 西南交通大学 | 一种用于全球导航卫星系统的动态周跳探测与修复方法 |
| CN104101890A (zh) * | 2014-07-01 | 2014-10-15 | 同济大学 | 一种对流层长距离亚厘米级实时动态卫星导航定位方法 |
| CN105549056A (zh) * | 2015-12-03 | 2016-05-04 | 中国电子科技集团公司第二十研究所 | 一种相对定位装置及其载波整周模糊度解算方法 |
| CN107102346A (zh) * | 2017-06-08 | 2017-08-29 | 中国电子科技集团公司第五十四研究所 | 一种基于北斗系统的多天线测姿方法 |
Cited By (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN109917356A (zh) * | 2019-03-13 | 2019-06-21 | 武汉际上导航科技有限公司 | 一种机载激光扫描系统误差标定方法 |
| CN109917356B (zh) * | 2019-03-13 | 2022-10-28 | 武汉际上导航科技有限公司 | 一种机载激光扫描系统误差标定方法 |
| CN112083464A (zh) * | 2019-06-14 | 2020-12-15 | 北京合众思壮科技股份有限公司 | 一种部分模糊度的固定方法及装置 |
| CN112083464B (zh) * | 2019-06-14 | 2023-12-26 | 北京合众思壮科技股份有限公司 | 一种部分模糊度的固定方法及装置 |
| CN110398758A (zh) * | 2019-07-24 | 2019-11-01 | 广州中海达卫星导航技术股份有限公司 | 实时钟差估计中的粗差探测方法、装置、设备及存储介质 |
| CN111045062A (zh) * | 2019-11-29 | 2020-04-21 | 航天恒星科技有限公司 | 一种基于电磁星的星基电离层反演方法 |
| CN113093236A (zh) * | 2019-12-23 | 2021-07-09 | 中国石油天然气集团有限公司 | 动态后处理的方法和装置 |
| CN112381309A (zh) * | 2020-11-23 | 2021-02-19 | 珠江水利委员会珠江水利科学研究院 | 水库大坝安全监测预警方法、装置、系统及存储介质 |
| CN112381309B (zh) * | 2020-11-23 | 2022-04-12 | 珠江水利委员会珠江水利科学研究院 | 水库大坝安全监测预警方法、装置、系统及存储介质 |
| CN112526573A (zh) * | 2021-02-07 | 2021-03-19 | 腾讯科技（深圳）有限公司 | 对象定位方法和装置、存储介质及电子设备 |
| CN114594502A (zh) * | 2022-02-18 | 2022-06-07 | 真点科技(北京)有限公司 | 一种实现周跳处理的方法、装置、计算机存储介质及终端 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN109143298B (zh) | 2020-08-07 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109143298B (zh) | 北斗和gps观测值周跳探测与修复方法、设备及存储设备 | |
| RU2479855C2 (ru) | Зависящее от расстояния уменьшение ошибки при определении местоположения в режиме кинематики реального времени | |
| US6127968A (en) | On-the-fly RTK positioning system with single frequency receiver | |
| CN108169774B (zh) | 支持rtppp和rtk的多模gnss单频周跳探测与修复方法 | |
| CN111308528B (zh) | 一种北斗/gps紧组合虚拟参考站定位方法 | |
| RU2354991C2 (ru) | Способ использования трех частот gps для разрешения целочисленных неоднозначностей фазы несущей | |
| Li | Cycle slip detection and ambiguity resolution algorithms for dual-frequency GPS data processing | |
| CN109765589B (zh) | 一种基于无电离层组合的三频gnss实时周跳固定技术 | |
| Lee et al. | The performance of RTK-GPS for surveying under challenging environmental conditions | |
| CN105699999B (zh) | 一种固定北斗地基增强系统基准站窄巷模糊度的方法 | |
| CN103837879A (zh) | 基于北斗系统民用载波相位组合实现高精度定位的方法 | |
| CN104483688A (zh) | 基于北斗卫星导航系统的高精度基线解算方法 | |
| CN110346823B (zh) | 可用于北斗精密单点定位的三频模糊度解算方法 | |
| CN105738926A (zh) | 一种glonass系统接收机间相位频间偏差标定方法 | |
| CN110687556A (zh) | 一种适用于laas的多径误差模型化方法 | |
| CN115373005A (zh) | 卫星导航信号间高精度产品转化方法 | |
| CN109597105A (zh) | 一种顾及载波系统间偏差的gps/glonass紧组合定位方法 | |
| CN112731496A (zh) | 一种面向智能终端的gnss精密单点定位数据质量控制方法 | |
| CN115933356A (zh) | 一种虚拟原子钟的高精度时间同步系统和方法 | |
| Zhou et al. | Multi-frequency BDS-3 real-time positioning performance assessment using new PPP-B2b augmentation service | |
| Li et al. | Toward Wide-Area and High-Precision Positioning With LEO Constellation Augmented PPP-RTK | |
| Zeng et al. | GPS triple-frequency undifferenced and uncombined precise orbit determination with the consideration of receiver time-variant bias | |
| Wang et al. | Comparison of three widely used multi‐GNSS real‐time single‐frequency precise point positioning models using the International GNSS Service real‐time service | |
| CN112924992B (zh) | 一种geo轨道精度评估方法、装置、电子设备及存储介质 | |
| CN106707311B (zh) | 基于gps增强的glonass rtk定位方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |