CN110727002A - 一种基于稀疏正则化的单频单站动态gnss载波相位信号周跳修复方法 - Google Patents

一种基于稀疏正则化的单频单站动态gnss载波相位信号周跳修复方法 Download PDF

Info

Publication number
CN110727002A
CN110727002A CN201910891666.3A CN201910891666A CN110727002A CN 110727002 A CN110727002 A CN 110727002A CN 201910891666 A CN201910891666 A CN 201910891666A CN 110727002 A CN110727002 A CN 110727002A
Authority
CN
China
Prior art keywords
cycle slip
epoch
observation
carrier phase
error
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.)
Pending
Application number
CN201910891666.3A
Other languages
English (en)
Inventor
常国宾
钱妮佳
陈超
张书毕
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
China University of Mining and Technology CUMT
Original Assignee
China University of Mining and Technology CUMT
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by China University of Mining and Technology CUMT filed Critical China University of Mining and Technology CUMT
Priority to CN201910891666.3A priority Critical patent/CN110727002A/zh
Publication of CN110727002A publication Critical patent/CN110727002A/zh
Pending legal-status Critical Current

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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/00Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
    • G01S19/01Satellite radio beacon positioning systems transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
    • G01S19/13Receivers
    • G01S19/23Testing, monitoring, correcting or calibrating of receiver elements

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

本发明公开了一种基于稀疏正则化的单频单站动态GNSS载波相位信号周跳修复方法,步骤如下:S1:采集任一频点上预设时段内的GNSS载波相位观测量和伪距观测量;S2:构建时间差分载波相位观测量TDCP和时间差分伪距观测量TDPR;S3:在相应历元处,对TDCP和TDPR进行无几何线性组合;S4:根据无几何时间差分观测量,逐历元修复大周跳;S5:将大周跳修复后的无几何时间差分观测量基于稀疏正则化的周跳估计方法估计存在的周跳的浮点解,对浮点解进行取整操作,获取周跳的整数解,进行周跳修复。本发明在采用FISTA数值迭代方法的同时,还充分挖掘了周跳稀疏性这一先验条件,并充分考虑观测量历元间的相关性,从而无需肉眼观测以及其他人为修改,具有高度自动化。

Description

一种基于稀疏正则化的单频单站动态GNSS载波相位信号周跳 修复方法
技术领域
本发明涉及GNSS卫星导航数据技术领域,尤其涉及一种基于稀疏正则化的单频单站动态GNSS载波相位信号周跳修复方法。
背景技术
GNSS可以提供载波相位和伪距两种重要的观测量,其中前者的精度比后者高两个数量级,所以在很多高精度GNSS应用中载波相位观测是必不可少的。然而,载波相位存在整周模糊性,这使得载波相位的模糊度解算成为一个必要的数据处理步骤。模糊解算的成功实施一般需要在一段时间段内载波相位信号是连续跟踪的,即不存在周跳。然而现实中由于观测环境的不理想,周跳的存在是不可避免的。有鉴于此,周跳的探测和修复则成为GNSS高精度应用中的一个必要的预处理步骤。
发明内容
发明目的:针对在单频、单站、动态的情况下,如何修复GNSS周跳的问题,本发明提出一种基于稀疏正则化的单频单站动态GNSS载波相位信号周跳修复方法。
技术方案:为实现本发明的目的,本发明所采用的技术方案是:
一种基于稀疏正则化的单频单站动态GNSS载波相位信号周跳修复方法,所述周跳修复方法具体包括如下步骤:
S1:采集任一频点上预设时间段内的GNSS载波相位观测量和伪距观测量;
S2:对所述GNSS载波相位观测量和伪距观测量分别进行历元间差分操作,构建时间差分载波相位观测量和时间差分伪距观测量;
S3:在相应历元处,对所述时间差分载波相位观测量和时间差分伪距观测量进行无几何线性组合,构建无几何时间差分观测量;
S4:根据所述无几何时间差分观测量,逐历元进行大周跳修复;
S5:将所述大周跳修复后的无几何时间差分观测量基于稀疏正则化的周跳估计方法估计存在的周跳的浮点解,对所述浮点解进行四舍五入取整操作,获取周跳的整数解,进行周跳修复。
进一步地讲,所述时间差分载波相位观测量和时间差分伪距观测量,具体为:
Figure BDA0002208933720000011
其中:γ=1/λ,
为时间差分载波相位观测量,dk为时间差分伪距观测量,λ为载波波长,xk为第k历元处的几何量增量,yk为第k历元处的周跳,εk为第k历元处伪距的观测误差,εk-1为第k-1历元处伪距的观测误差,ek为第k历元处载波相位的观测误差,ek-1为第k-1历元处载波相位的观测误差。
进一步地讲,所述观测误差具体为:
Figure BDA0002208933720000023
其中:
Figure BDA0002208933720000024
εk-1为第k-1历元处伪距的观测误差,ek-1为第k-1历元处载波相位的观测误差,ek为第k历元处载波相位的观测误差,εk为第k历元处伪距的观测误差。
进一步地讲,所述所述无几何时间差分观测量的观测方程,具体为:
其中:γ=1/λ,
Figure BDA0002208933720000026
wk为无几何时间差分观测量,λ为载波波长,
Figure BDA0002208933720000027
为时间差分载波相位观测量,dk为时间差分伪距观测量,yk为第k历元处的周跳,εk-1为第k-1历元处伪距的观测误差,ek-1为第k-1历元处载波相位的观测误差,ek为第k历元处载波相位的观测误差,εk为第k历元处伪距的观测误差,ηk为观测误差。
进一步地讲,所述步骤S4中逐历元进行大周跳修复,具体如下:
S4.1:根据所述无几何时间差分观测量的观测方程,所述无几何时间差分观测量即为第k历元处周跳的直接观测量,在周跳修复成功率满足固定成功率要求时,对所述周跳的浮点解进行四舍五入取整操作,具体为:
Figure BDA0002208933720000028
其中:
Figure BDA0002208933720000029
为第k历元处周跳的浮点解进行四舍五入取整操作后的值,[]表示四舍五入取整算子,wk为无几何时间差分观测量;
S4.2:根据所述四舍五入取整操作后的周跳浮点解,对所述周跳进行修复,并将所述大周跳修复后的观测量表示为ξk,所述大周跳修复后的观测量对应的观测方程具体为:
Figure BDA0002208933720000031
其中:ξk为大周跳修复后的观测量,wk为无几何时间差分观测量,
Figure BDA0002208933720000032
为第k历元处周跳的浮点解进行四舍五入取整操作后的值,yk为第k历元处的周跳,ηk为观测误差,zk为待求解参数。
进一步地讲,所述步骤S5中获取周跳的整数解,具体如下:
S5.1:根据所述大周跳修复后的观测量对应的观测方程,采用所有所述大周跳修复后的无几何时间差分观测量,构建向量形式的观测方程,具体为:
ξ=z+η
其中:ξ为所有ξk组成的观测向量,ξk为大周跳修复后的观测量,z为所有zk组成的未知数参数向量,zk为待求解参数,η为所有ηk组成的观测误差向量,ηk为观测误差;
S5.2:将所述观测误差向量对应的协方差矩阵标记为Qηη,所述协方差矩阵具体为:
Figure BDA0002208933720000033
其中:
Figure BDA0002208933720000034
γ=1/λ,
Figure BDA0002208933720000035
Qηη为观测误差向量对应的协方差矩阵,
Figure BDA0002208933720000036
为伪距的观测误差εk的方差,为载波相位的观测误差ek的方差,λ为载波波长,ek为第k历元处载波相位的观测误差,εk为第k历元处伪距的观测误差,n表示历元的总个数;
S5.3:预设N个超参数,并从N个所述超参数中选出最优的超参数;
S5.4:将所述最优超参数对应的周跳浮点解作为最终的浮点解,并将所述最终浮点解进行四舍五入取整操作,获取所述周跳的整数解,完成周跳修复。
进一步地讲,所述步骤S5.3中从N个超参数中选出最优的超参数,具体如下:
S5.3.1:预设N个超参数,并从N个所述超参数中任意选取一个超参数;
S5.3.2:根据所述选取的超参数,未知数参数向量通过FISTA算法使稀疏正则化代价函数最小,获取稀疏正则化代价函数最小时对应的浮点解,所述稀疏正则化代价函数,具体为:
Figure BDA0002208933720000041
其中:f为稀疏正则化代价函数,ξ为所有ξk组成的观测向量,ξk为大周跳修复后的观测量,z为所有zk组成的未知数参数向量,zk为待求解参数,Qηη为观测误差向量对应的协方差矩阵,μ为超参数;
S5.3.3:将所述稀疏正则化代价函数最小时对应的浮点解代入广义交叉检核公式中,获取所述选取的超参数对应的GCV值,所述广义交叉检核公式具体为:
Figure BDA0002208933720000042
其中:GCV(μ)为超参数μ对应的GCV值,n表示历元的总个数,n0
Figure BDA0002208933720000043
中非零元素的个数,ξ为所有ξk组成的观测向量,ξk为大周跳修复后的观测量,
Figure BDA0002208933720000044
为稀疏正则化代价函数最小时对应的浮点解,Qηη为观测误差向量对应的协方差矩阵;
S5.3.4:从预设的N个所述超参数中选出一个未被选择过的超参数,并根据所述重新选择的超参数,重复步骤S5.3.2-步骤S5.3.3,直至预设的N个所述超参数中的所有超参数均被选择过,并将每个所述超参数对应的GCV值进行比较,从中选出最小GCV值,所述最小GCV值对应的超参数即为最优超参数。
进一步地讲,所述步骤S5.3.2中未知数参数向量通过FISTA算法使稀疏正则化代价函数最小,具体如下:
S5.3.2.1:将待求解参数进行初始化,具体为:
Figure BDA0002208933720000051
其中:z0为迭代前待求解参数的初步估计值,θ1和s1均为第一次迭代运算引入的辅助变量,ξ为所有ξk组成的观测向量,ξk为大周跳修复后的观测量;
S5.3.2.2:将所述待求解参数的初步估计值、第一次迭代运算引入的辅助变量均进行迭代计算,直至算法收敛,确定出每次迭代过程中所述待求解参数的大小,具体为:
Figure BDA0002208933720000052
其中:
Figure BDA0002208933720000053
zk和zk-1为待求解参数,
Figure BDA0002208933720000054
Figure BDA0002208933720000055
的软阈值函数,θk和sk为第k次迭代运算引入的辅助变量,λmax为Qηη的最大特征值,Qηη为观测误差向量对应的协方差矩阵,λmin为Qηη的最小特征值,ξ为所有ξk组成的观测向量,ξk为大周跳修复后的观测量,sk+1和θk+1为第k+1次迭代运算引入的辅助变量,μ为超参数。
进一步地讲,所述软阈值函数
Figure BDA0002208933720000057
算子中
Figure BDA0002208933720000058
通过高效的三对角矩阵算法TMA进行计算,令则g=Qηηf,同时依次执行如下公式,具体为:
Figure BDA00022089337200000510
Figure BDA00022089337200000511
Figure BDA0002208933720000061
其中:
Figure BDA0002208933720000062
γ=1/λ
Figure BDA0002208933720000063
Figure BDA0002208933720000064
均为辅助变量,fk为向量f的第k个元素,gk为向量g的第k个元素,
Figure BDA0002208933720000065
为伪距的观测误差εk的方差,
Figure BDA0002208933720000066
为载波相位的观测误差ek的方差,λ为载波波长,ek为第k历元处载波相位的观测误差,εk为第k历元处伪距的观测误差,n表示历元的总个数。
进一步地讲,所述软阈值函数
Figure BDA0002208933720000067
的求解过程,具体如下:
第一步:给定向量x,对给定的向量x进行软阈值操作,得到一个同维数的向量,所述同维数向量的第j个元素,具体为:
Figure BDA0002208933720000068
其中:
Figure BDA0002208933720000069
x表示软阈值函数中的自变量,(|xj|-tμ)+为铰链函数,sgn(xj)为符号函数,μ为超参数,λmax为Qηη的最大特征值,Qηη为观测误差向量对应的协方差矩阵,λmin为Qηη的最小特征值,θk为第k次迭代运算引入的辅助变量,ξ为所有ξk组成的观测向量,ξk为大周跳修复后的观测量;
第二步:对所述铰链函数进行求解,具体为:
Figure BDA00022089337200000610
其中:
Figure BDA00022089337200000611
(|xj|-tμ)+为铰链函数,μ为超参数,λmax为Qηη的最大特征值,Qηη为观测误差向量对应的协方差矩阵,λmin为Qηη的最小特征值,xj为x的第j个分量,x表示软阈值函数中的自变量;
第三步:对所述符号函数进行求解,具体为:
Figure BDA0002208933720000071
其中:sgn(xj)为符号函数,xj为x的第j个分量,x表示软阈值函数中的自变量。
有益效果:与现有技术相比,本发明的技术方案具有以下有益技术效果:
本发明的周跳修复方法适用于单频、单站、动态的场合,在采用FISTA数值迭代方法的同时,还充分挖掘了周跳稀疏性这一先验条件,并充分考虑观测量历元间的相关性,从而针对性地采用高效的三对角矩阵算法,进而无需肉眼观测以及其他人为残余,具有高度自动化。
附图说明
图1是本发明周跳修复方法的流程示意图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚,下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述。其中,所描述的实施例是本发明一部分实施例,而不是全部的实施例。因此,以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围,而是仅仅表示本发明的选定实施例。
实施例1
参考图1,本实施例提供了一种基于稀疏正则化的单频单站动态GNSS载波相位信号周跳修复方法,其中单频是指只采用单频观测,即单一频点上的载波相位观测和伪距观测,该单一频点可以是原始观测中的某一频点,也可以是多频情形时多频信号的某一线性组合。单站是指不需要参考站的辅助,可适用于只有流动站的单站绝对定位情况。动态是指适用于流动站处于动态的情形。步骤具体如下:
步骤S1:采集任一频点上预设时间段内的GNSS载波相位观测量和伪距观测量。该GNSS可以是GPS、Beidou、Galileo等任一码分多址卫星导航系统。
其中任一频点可以是单频情况下的观测量,也可以是多频情况下多频观测的任意一种线性组合,譬如:宽巷或窄巷组合。
同样地,预设时间段内可以是任意一颗卫星的连续观测时间,也可以是该连续观测时间中的一部分。
步骤S2:构建时间差分载波相位观测量TDCP即将GNSS载波相位观测量在相邻历元间做差,并将其标记为
Figure BDA0002208933720000081
构建时间差分伪距观测量TDPR即将伪距观测量在相邻历元间做差,并将其标记为dk。在
Figure BDA0002208933720000082
和dk中k=1,2,...,n,其中:k表示历元的序号,n表示历元的总个数。具体地讲,时间差分过程可以消去大部分历元间共模误差或时间慢变误差,譬如:电离层延迟误差和多路径误差。
在本实施例中,时间差分载波相位观测量TDCP和时间差分伪距观测量TDPR,具体为:
Figure BDA0002208933720000083
其中:γ=1/λ,
Figure BDA0002208933720000084
Figure BDA0002208933720000085
为时间差分载波相位观测量,dk为时间差分伪距观测量,λ为载波波长,xk为第k历元处的几何量增量,yk为第k历元处的周跳,εk为第k历元处伪距的观测误差,εk-1为第k-1历元处伪距的观测误差,ek为第k历元处载波相位的观测误差,ek-1为第k-1历元处载波相位的观测误差。
具体地讲,第k历元处的周跳yk即为本实施例中致力于求解的未知数。且该未知数是否存在也是未知的的,当该未知数存在时,则该第k历元处的周跳yk为不为零的整数,当该未知数不存在时,则该第k历元处的周跳yk为零。
步骤S3:在相应历元处,对步骤S2中时间差分载波相位观测量TDCP和时间差分伪距观测量TDPR进行无几何线性组合。无几何线性组合操作可以消除卫星-载体间的几何距离、对流层延迟误差、卫星轨道误差等非弥散,即与频率无关的变量。
同时进行无几何线性组合处理后的时间差分载波相位观测量TDCP可以认为仅包含周跳和观测噪声两部分。其中观测噪声包含所有未建模、未消除的剩余误差。
进行无几何线性组合处理后的时间差分载波相位观测量TDCP和时间差分伪距观测量TDPR,可以按照下式构建无几何时间差分观测量,从而无几何时间差分观测量的观测方程,具体为:
其中:γ=1/λ,
Figure BDA0002208933720000091
wk为无几何时间差分观测量,λ为载波波长,
Figure BDA0002208933720000092
为时间差分载波相位观测量,dk为时间差分伪距观测量,yk为第k历元处的周跳,εk-1为第k-1历元处伪距的观测误差,ek-1为第k-1历元处载波相位的观测误差,ek为第k历元处载波相位的观测误差,εk为第k历元处伪距的观测误差,ηk为观测误差。
具体地讲,观测误差ηk具体为:
Figure BDA0002208933720000093
其中:
Figure BDA0002208933720000094
εk-1为第k-1历元处伪距的观测误差,ek-1为第k-1历元处载波相位的观测误差,ek为第k历元处载波相位的观测误差,εk为第k历元处伪距的观测误差。
步骤S4:根据时间差分载波相位观测量TDCP和时间差分伪距观测量TDPR构建的无几何时间差分观测量,逐历元进行大周跳修复。
在逐历元进行大周跳修复时,需要根据信号波长或频率以及观测量噪声特性,在考虑成功率的前提下,进行大周跳修复,即确定周跳修复成功率大于等于给定值时才可进行修复,否则不进行修复。在本实施例中,给定值选择为99%。
这是由于:处理后的大周跳要么消失、要么变为小周跳,在后续步骤中,周跳越接近于1,则周跳向量的L1范数就越接近于L0范数。同时考虑修复成功率的目的是为了尽可能避免因为错误修复而人为引入周跳,过多地人为引入周跳将破坏周跳稀疏性,并致使步骤S5无法进行。具体过程如下:
步骤S4.1:根据无几何时间差分观测量的观测方程,可以发现,无几何时间差分观测量wk即为要求解的未知数第k历元处的周跳yk的直接观测量,从而可以直接令第k历元处的周跳yk等于无几何时间差分观测量wk,进而获取该周跳的浮点解。当确定周跳修复成功率满足固定成功率要求时,对该周跳的浮点解进行四舍五入取整操作,具体为:
Figure BDA0002208933720000095
其中:为第k历元处周跳的浮点解进行四舍五入取整操作后的值,[]表示四舍五入取整算子,wk为无几何时间差分观测量。
在本实施例中,关于固定成功率要求,譬如:对于GPS L1频点,上述过程中浮点周跳精度则约为5,从而对大于等于5的周跳进行四舍五入取整操作。
步骤S4.2:对第k历元处的周跳yk进行修复,令大周跳修复后的观测量表示为ξk,则大周跳修复后的观测量对应的观测方程,具体为:
Figure BDA0002208933720000101
其中:ξk为大周跳修复后的观测量,wk为无几何时间差分观测量,
Figure BDA0002208933720000102
为第k历元处周跳的浮点解进行四舍五入取整操作后的值,yk为第k历元处的周跳,ηk为观测误差,zk为待求解参数。
具体地讲,待求解参数zk为下一个需要求解的新周跳。从而一旦获取了待求解参数zk的大小,即可得到原始周跳yk的大小。
步骤S5:将大周跳修复后的无几何时间差分观测量ξk基于稀疏正则化的周跳估计方法估计存在的周跳的浮点解,并对该浮点解进行四舍五入取整操作,获取周跳的整数解,从而进行周跳修复。具体如下:
步骤S5.1:根据步骤S4.2中大周跳修复后的观测量对应的观测方程,采用所有大周跳修复后的无几何时间差分观测量ξk,构建如下向量形式的观测方程,具体为:
ξ=z+η
其中:ξ为所有ξk组成的观测向量,ξk为大周跳修复后的观测量,z为所有zk组成的未知数参数向量,zk为待求解参数,η为所有ηk组成的观测误差向量,ηk为观测误差。
步骤S5.2:将所有ηk组成的观测误差向量η对应的协方差矩阵标记为Qηη,该协方差矩阵标记为Qηη具体为:
Figure BDA0002208933720000111
其中:
Figure BDA0002208933720000112
γ=1/λ,
Qηη为观测误差向量对应的协方差矩阵,为伪距的观测误差εk的方差,
Figure BDA0002208933720000115
为载波相位的观测误差ek的方差,λ为载波波长,ek为第k历元处载波相位的观测误差,εk为第k历元处伪距的观测误差,n表示历元的总个数。
具体地讲,当伪距的观测误差εk的方差载波相位的观测误差ek的方差
Figure BDA0002208933720000117
均为常数时,则观测误差向量对应的协方差矩阵Qηη为对称三对角Toplitz矩阵。
步骤S5.3:预设N个超参数μ,并根据如下步骤从N个超参数μ中选出最优的超参数μ。具体如下:
步骤S5.3.1:预设N个超参数μ,并从N个超参数μ中任意选取一个超参数μ。
步骤S5.3.2:根据选出的超参数μ,所有zk组成的未知数参数向量z通过FISTA算法使如下稀疏正则化代价函数最小,获取稀疏正则化代价函数最小时对应的浮点解
Figure BDA0002208933720000118
其中以修复大周跳后的时间差分载波相位和时间差分伪距的无几何组合为观测量,以可能存在的周跳为待估参数,采用L1范数正则化法进行参数的稀疏求解,即采用加权残差平方和与参数L1范数之线性组合作为代价函数,组合系数又称为正则化参数,从而待估参数的求解过程即为此代价函数的最小化过程。
稀疏正则化代价函数,具体为:
Figure BDA0002208933720000121
其中:f为稀疏正则化代价函数,ξ为所有ξk组成的观测向量,ξk为大周跳修复后的观测量,z为所有zk组成的未知数参数向量,zk为待求解参数,Qηη为观测误差向量对应的协方差矩阵,μ为超参数。
在构建稀疏正则化代价函数时,其中的加权残差平方和构建过程中的协方差矩阵考虑了观测量相邻历元间的相关性,该相关性来自于步骤S2中相邻历元间的差分操作,从而其中的协方差矩阵为三对角对称矩阵。
在本实施例中,所有zk组成的未知数参数向量z通过FISTA算法,使上述稀疏正则化代价函数最小,具体如下:
步骤S5.3.2.1:将待求解参数zk进行初始化,具体为:
Figure BDA0002208933720000122
其中:z0为迭代前待求解参数的初步估计值,θ1和s1均为第一次迭代运算引入的辅助变量,ξ为所有ξk组成的观测向量,ξk为大周跳修复后的观测量。
步骤S5.3.2.2:将待求解参数的初步估计值z0、第一次迭代运算引入的辅助变量θ1和s1均进行迭代计算,直至算法收敛,确定出每次迭代过程中待求解参数的大小,具体为:
Figure BDA0002208933720000123
其中:
Figure BDA0002208933720000124
zk和zk-1为待求解参数,
Figure BDA0002208933720000126
的软阈值函数,θk和sk为第k次迭代运算引入的辅助变量,λmax为Qηη的最大特征值,Qηη为观测误差向量对应的协方差矩阵,λmin为Qηη的最小特征值,ξ为所有ξk组成的观测向量,ξk为大周跳修复后的观测量,sk+1和θk+1为第k+1次迭代运算引入的辅助变量,μ为超参数。
在本实施例中,算法收敛即为相邻两次迭代获取的待求解参数的之间差值的绝对值小于预设阈值。
具体地讲,
Figure BDA0002208933720000131
Figure BDA0002208933720000132
通过高效的三对角矩阵算法TMA进行计算,令
Figure BDA0002208933720000133
则g=Qηηf。从而的求解即为求解方程g=Qηηf,求解方程g=Qηηf需要依次执行如下公式,具体为:
Figure BDA0002208933720000135
Figure BDA0002208933720000136
其中:
Figure BDA0002208933720000137
γ=1/λ
Figure BDA0002208933720000138
Figure BDA0002208933720000139
均为辅助变量,fk为向量f的第k个元素,gk为向量g的第k个元素,
Figure BDA00022089337200001310
为伪距的观测误差εk的方差,
Figure BDA00022089337200001311
为载波相位的观测误差ek的方差,λ为载波波长,ek为第k历元处载波相位的观测误差,εk为第k历元处伪距的观测误差,n表示历元的总个数。
在本实施例中,软阈值函数
Figure BDA00022089337200001312
的求解过程,具体如下:
第一步:给定向量x,对给定的向量x进行软阈值操作,得到一个同维数的向量,该同维数向量的第j个元素,具体为:
Figure BDA00022089337200001313
其中:
Figure BDA0002208933720000141
x表示软阈值函数中的自变量,(|xj|-tμ)+为铰链函数,sgn(xj)为符号函数,μ为超参数,λmax为Qηη的最大特征值,Qηη为观测误差向量对应的协方差矩阵,λmin为Qηη的最小特征值,θk为第k次迭代运算引入的辅助变量,ξ为所有ξk组成的观测向量,ξk为大周跳修复后的观测量。
第二步:对铰链函数(|xj|-tμ)+进行求解,具体为:
Figure BDA0002208933720000142
其中:
Figure BDA0002208933720000143
(|xj|-tμ)+为铰链函数,μ为超参数,λmax为Qηη的最大特征值,Qηη为观测误差向量对应的协方差矩阵,λmin为Qηη的最小特征值,xj为x的第j个分量,x表示软阈值函数中的自变量。
第三步:对符号函数sgn(xj)进行求解,具体为:
Figure BDA0002208933720000144
其中:sgn(xj)为符号函数,xj为x的第j个分量,x表示软阈值函数中的自变量。
步骤S5.3.3:将稀疏正则化代价函数最小时对应的浮点解代入广义交叉检核公式中,获取步骤S5.3.1中任意选取的超参数μ对应的GCV值,该广义交叉检核公式具体为:
Figure BDA0002208933720000146
其中:GCV(μ)为超参数μ对应的GCV值,n表示历元的总个数,n0
Figure BDA0002208933720000147
中非零元素的个数,ξ为所有ξk组成的观测向量,ξk为大周跳修复后的观测量,
Figure BDA0002208933720000148
为稀疏正则化代价函数最小时对应的浮点解,Qηη为观测误差向量对应的协方差矩阵。
步骤S5.3.4:从预设的N个超参数μ中选出一个未被选择过的超参数μ,并根据该重新选择的超参数μ,重复步骤S5.3.2-步骤S5.3.3,直至预设的N个超参数μ中的所有超参数均被选择过。
同时将每个超参数μ对应的GCV值进行比较,从中选出最小GCV值,并将该最小GCV值对应的超参数μ作为最优超参数。
步骤S5.4:根据选出的最优超参数μ,将该最优超参数μ对应的周跳浮点解作为最终的浮点解,同时将该最终浮点解进行四舍五入取整操作获取该周跳的整数解,从而完成周跳修复。
以上示意性的对本发明及其实施方式进行了描述,该描述没有限制性,附图中所示的也只是本发明的实施方式之一,实际的结构和方法并不局限于此。所以,如果本领域的普通技术人员受其启示,在不脱离本发明创造宗旨的情况下,不经创造性的设计出与该技术方案相似的结构方式及实施例,均属于本发明的保护范围。

Claims (10)

1.一种基于稀疏正则化的单频单站动态GNSS载波相位信号周跳修复方法,其特征在于,所述周跳修复方法具体包括如下步骤:
S1:采集任一频点上预设时间段内的GNSS载波相位观测量和伪距观测量;
S2:对所述GNSS载波相位观测量和伪距观测量分别进行历元间差分操作,构建时间差分载波相位观测量和时间差分伪距观测量;
S3:在相应历元处,对所述时间差分载波相位观测量和时间差分伪距观测量进行无几何线性组合,构建无几何时间差分观测量;
S4:根据所述无几何时间差分观测量,逐历元进行大周跳修复;
S5:将所述大周跳修复后的无几何时间差分观测量基于稀疏正则化的周跳估计方法估计存在的周跳的浮点解,对所述浮点解进行四舍五入取整操作,获取周跳的整数解,进行周跳修复。
2.根据权利要求1所述的一种基于稀疏正则化的单频单站动态GNSS载波相位信号周跳修复方法,其特征在于,所述时间差分载波相位观测量和时间差分伪距观测量,具体为:
Figure FDA0002208933710000011
其中:γ=1/λ,
Figure FDA0002208933710000012
Figure FDA0002208933710000013
为时间差分载波相位观测量,dk为时间差分伪距观测量,λ为载波波长,xk为第k历元处的几何量增量,yk为第k历元处的周跳,εk为第k历元处伪距的观测误差,εk-1为第k-1历元处伪距的观测误差,ek为第k历元处载波相位的观测误差,ek-1为第k-1历元处载波相位的观测误差。
3.根据权利要求2所述的一种基于稀疏正则化的单频单站动态GNSS载波相位信号周跳修复方法,其特征在于,所述观测误差具体为:
Figure FDA0002208933710000014
其中:
Figure FDA0002208933710000015
εk-1为第k-1历元处伪距的观测误差,ek-1为第k-1历元处载波相位的观测误差,ek为第k历元处载波相位的观测误差,εk为第k历元处伪距的观测误差。
4.根据权利要求1或2或3所述的一种基于稀疏正则化的单频单站动态GNSS载波相位信号周跳修复方法,其特征在于,所述所述无几何时间差分观测量的观测方程,具体为:
Figure FDA0002208933710000021
其中:γ=1/λ,
Figure FDA0002208933710000022
wk为无几何时间差分观测量,λ为载波波长,
Figure FDA0002208933710000023
为时间差分载波相位观测量,dk为时间差分伪距观测量,yk为第k历元处的周跳,εk-1为第k-1历元处伪距的观测误差,ek-1为第k-1历元处载波相位的观测误差,ek为第k历元处载波相位的观测误差,εk为第k历元处伪距的观测误差,ηk为观测误差。
5.根据权利要求4所述的一种基于稀疏正则化的单频单站动态GNSS载波相位信号周跳修复方法,其特征在于,所述步骤S4中逐历元进行大周跳修复,具体如下:
S4.1:根据所述无几何时间差分观测量的观测方程,所述无几何时间差分观测量即为第k历元处周跳的直接观测量,在周跳修复成功率满足固定成功率要求时,对所述周跳的浮点解进行四舍五入取整操作,具体为:
其中:
Figure FDA0002208933710000025
为第k历元处周跳的浮点解进行四舍五入取整操作后的值,[]表示四舍五入取整算子,wk为无几何时间差分观测量;
S4.2:根据所述四舍五入取整操作后的周跳浮点解,对所述周跳进行修复,并将所述大周跳修复后的观测量表示为ξk,所述大周跳修复后的观测量对应的观测方程具体为:
Figure FDA0002208933710000026
其中:ξk为大周跳修复后的观测量,wk为无几何时间差分观测量,为第k历元处周跳的浮点解进行四舍五入取整操作后的值,yk为第k历元处的周跳,ηk为观测误差,zk为待求解参数。
6.根据权利要求5所述的一种基于稀疏正则化的单频单站动态GNSS载波相位信号周跳修复方法,其特征在于,所述步骤S5中获取周跳的整数解,具体如下:
S5.1:根据所述大周跳修复后的观测量对应的观测方程,采用所有所述大周跳修复后的无几何时间差分观测量,构建向量形式的观测方程,具体为:
ξ=z+η
其中:ξ为所有ξk组成的观测向量,ξk为大周跳修复后的观测量,z为所有zk组成的未知数参数向量,zk为待求解参数,η为所有ηk组成的观测误差向量,ηk为观测误差;
S5.2:将所述观测误差向量对应的协方差矩阵标记为Qηη,所述协方差矩阵具体为:
其中:
Figure FDA0002208933710000032
γ=1/λ,
Figure FDA0002208933710000033
Qηη为观测误差向量对应的协方差矩阵,为伪距的观测误差εk的方差,
Figure FDA0002208933710000035
为载波相位的观测误差ek的方差,λ为载波波长,ek为第k历元处载波相位的观测误差,εk为第k历元处伪距的观测误差,n表示历元的总个数;
S5.3:预设N个超参数,并从N个所述超参数中选出最优的超参数;
S5.4:将所述最优超参数对应的周跳浮点解作为最终的浮点解,并将所述最终浮点解进行四舍五入取整操作,获取所述周跳的整数解,完成周跳修复。
7.根据权利要求6所述的一种基于稀疏正则化的单频单站动态GNSS载波相位信号周跳修复方法,其特征在于,所述步骤S5.3中从N个超参数中选出最优的超参数,具体如下:
S5.3.1:预设N个超参数,并从N个所述超参数中任意选取一个超参数;
S5.3.2:根据所述选取的超参数,未知数参数向量通过FISTA算法使稀疏正则化代价函数最小,获取稀疏正则化代价函数最小时对应的浮点解,所述稀疏正则化代价函数,具体为:
其中:f为稀疏正则化代价函数,ξ为所有ξk组成的观测向量,ξk为大周跳修复后的观测量,z为所有zk组成的未知数参数向量,zk为待求解参数,Qηη为观测误差向量对应的协方差矩阵,μ为超参数;
S5.3.3:将所述稀疏正则化代价函数最小时对应的浮点解代入广义交叉检核公式中,获取所述选取的超参数对应的GCV值,所述广义交叉检核公式具体为:
Figure FDA0002208933710000042
其中:GCV(μ)为超参数μ对应的GCV值,n表示历元的总个数,n0
Figure FDA0002208933710000043
中非零元素的个数,ξ为所有ξk组成的观测向量,ξk为大周跳修复后的观测量,为稀疏正则化代价函数最小时对应的浮点解,Qηη为观测误差向量对应的协方差矩阵;
S5.3.4:从预设的N个所述超参数中选出一个未被选择过的超参数,并根据所述重新选择的超参数,重复步骤S5.3.2-步骤S5.3.3,直至预设的N个所述超参数中的所有超参数均被选择过,并将每个所述超参数对应的GCV值进行比较,从中选出最小GCV值,所述最小GCV值对应的超参数即为最优超参数。
8.根据权利要求7所述的一种基于稀疏正则化的单频单站动态GNSS载波相位信号周跳修复方法,其特征在于,所述步骤S5.3.2中未知数参数向量通过FISTA算法使稀疏正则化代价函数最小,具体如下:
S5.3.2.1:将待求解参数进行初始化,具体为:
Figure FDA0002208933710000045
其中:z0为迭代前待求解参数的初步估计值,θ1和s1均为第一次迭代运算引入的辅助变量,ξ为所有ξk组成的观测向量,ξk为大周跳修复后的观测量;
S5.3.2.2:将所述待求解参数的初步估计值、第一次迭代运算引入的辅助变量均进行迭代计算,直至算法收敛,确定出每次迭代过程中所述待求解参数的大小,具体为:
Figure FDA0002208933710000051
其中:
Figure FDA0002208933710000052
zk和zk-1为待求解参数,
Figure FDA0002208933710000053
Figure FDA0002208933710000054
的软阈值函数,θk和sk为第k次迭代运算引入的辅助变量,λmax为Qηη的最大特征值,Qηη为观测误差向量对应的协方差矩阵,λmin为Qηη的最小特征值,ξ为所有ξk组成的观测向量,ξk为大周跳修复后的观测量,sk+1和θk+1为第k+1次迭代运算引入的辅助变量,μ为超参数。
9.根据权利要求8所述的一种基于稀疏正则化的单频单站动态GNSS载波相位信号周跳修复方法,其特征在于,所述软阈值函数
Figure FDA0002208933710000055
算子中
Figure FDA0002208933710000056
通过高效的三对角矩阵算法TMA进行计算,令
Figure FDA0002208933710000057
则g=Qηηf,同时依次执行如下公式,具体为:
Figure FDA0002208933710000058
Figure FDA00022089337100000510
其中:
Figure FDA00022089337100000511
γ=1/λ
Figure FDA00022089337100000512
Figure FDA00022089337100000513
均为辅助变量,fk为向量f的第k个元素,gk为向量g的第k个元素,
Figure FDA0002208933710000061
为伪距的观测误差εk的方差,
Figure FDA0002208933710000062
为载波相位的观测误差ek的方差,λ为载波波长,ek为第k历元处载波相位的观测误差,εk为第k历元处伪距的观测误差,n表示历元的总个数。
10.根据权利要求8所述的一种基于稀疏正则化的单频单站动态GNSS载波相位信号周跳修复方法,其特征在于,所述软阈值函数的求解过程,具体如下:
第一步:给定向量x,对给定的向量x进行软阈值操作,得到一个同维数的向量,所述同维数向量的第j个元素,具体为:
Figure FDA0002208933710000064
其中:
Figure FDA0002208933710000065
x表示软阈值函数中的自变量,(|xj|-tμ)+为铰链函数,sgn(xj)为符号函数,μ为超参数,λmax为Qηη的最大特征值,Qηη为观测误差向量对应的协方差矩阵,λmin为Qηη的最小特征值,θk为第k次迭代运算引入的辅助变量,ξ为所有ξk组成的观测向量,ξk为大周跳修复后的观测量;
第二步:对所述铰链函数进行求解,具体为:
Figure FDA0002208933710000066
其中:
Figure FDA0002208933710000067
(|xj|-tμ)+为铰链函数,μ为超参数,λmax为Qηη的最大特征值,Qηη为观测误差向量对应的协方差矩阵,λmin为Qηη的最小特征值,xj为x的第j个分量,x表示软阈值函数中的自变量;
第三步:对所述符号函数进行求解,具体为:
Figure FDA0002208933710000071
其中:sgn(xj)为符号函数,xj为x的第j个分量,x表示软阈值函数中的自变量。
CN201910891666.3A 2019-09-20 2019-09-20 一种基于稀疏正则化的单频单站动态gnss载波相位信号周跳修复方法 Pending CN110727002A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910891666.3A CN110727002A (zh) 2019-09-20 2019-09-20 一种基于稀疏正则化的单频单站动态gnss载波相位信号周跳修复方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910891666.3A CN110727002A (zh) 2019-09-20 2019-09-20 一种基于稀疏正则化的单频单站动态gnss载波相位信号周跳修复方法

Publications (1)

Publication Number Publication Date
CN110727002A true CN110727002A (zh) 2020-01-24

Family

ID=69219304

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910891666.3A Pending CN110727002A (zh) 2019-09-20 2019-09-20 一种基于稀疏正则化的单频单站动态gnss载波相位信号周跳修复方法

Country Status (1)

Country Link
CN (1) CN110727002A (zh)

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US11662478B2 (en) 2020-12-17 2023-05-30 Swift Navigation, Inc. System and method for fusing dead reckoning and GNSS data streams
US11860287B2 (en) 2022-03-01 2024-01-02 Swift Navigation, Inc. System and method for detecting outliers in GNSS observations
US11906640B2 (en) 2022-03-01 2024-02-20 Swift Navigation, Inc. System and method for fusing sensor and satellite measurements for positioning determination
JP7464790B2 (ja) 2021-03-24 2024-04-09 三菱電機株式会社 測位装置、測位プログラムおよび測位方法
US12019163B2 (en) 2022-09-12 2024-06-25 Swift Navigation, Inc. System and method for GNSS correction transmission

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130271318A1 (en) * 2012-04-12 2013-10-17 Trimble Navigation Limited Advanced global navigation satellite systems (gnss) positioning using precise satellite information
CN105676243A (zh) * 2016-01-11 2016-06-15 昆明理工大学 一种基于无几何相位和电离层残差法的北斗三频周跳探测方法
CN105933008A (zh) * 2016-04-15 2016-09-07 哈尔滨工业大学 基于聚集稀疏正则化正交匹配追踪算法的多频带信号重构方法
CN107728181A (zh) * 2017-10-10 2018-02-23 北京无线电计量测试研究所 一种实时周跳探测修复方法
CN108169774A (zh) * 2017-12-26 2018-06-15 北方信息控制研究院集团有限公司 支持rtppp和rtk的多模gnss单频周跳探测与修复方法
CN108196281A (zh) * 2017-11-22 2018-06-22 同济大学 一种基于位置域曲线约束的单频动态周跳探测与修复方法
CN110208836A (zh) * 2019-05-30 2019-09-06 东南大学 基于卡尔曼滤波的gnss高适应性周跳探测与修复方法

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20130271318A1 (en) * 2012-04-12 2013-10-17 Trimble Navigation Limited Advanced global navigation satellite systems (gnss) positioning using precise satellite information
CN105676243A (zh) * 2016-01-11 2016-06-15 昆明理工大学 一种基于无几何相位和电离层残差法的北斗三频周跳探测方法
CN105933008A (zh) * 2016-04-15 2016-09-07 哈尔滨工业大学 基于聚集稀疏正则化正交匹配追踪算法的多频带信号重构方法
CN107728181A (zh) * 2017-10-10 2018-02-23 北京无线电计量测试研究所 一种实时周跳探测修复方法
CN108196281A (zh) * 2017-11-22 2018-06-22 同济大学 一种基于位置域曲线约束的单频动态周跳探测与修复方法
CN108169774A (zh) * 2017-12-26 2018-06-15 北方信息控制研究院集团有限公司 支持rtppp和rtk的多模gnss单频周跳探测与修复方法
CN110208836A (zh) * 2019-05-30 2019-09-06 东南大学 基于卡尔曼滤波的gnss高适应性周跳探测与修复方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
GUOBIN CHANG 等: "Ionospheric delay prediction based on online polynomial modeling for real-time cycle slip repair of undifferenced triple-frequency GNSS signals", 《MEASUREMENT》 *
高书亮等: "适用于单频GPS用户的周跳探测方法", 《北京航空航天大学学报》 *

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US11662478B2 (en) 2020-12-17 2023-05-30 Swift Navigation, Inc. System and method for fusing dead reckoning and GNSS data streams
US12013472B2 (en) 2020-12-17 2024-06-18 Swift Navigation, Inc. System and method for fusing dead reckoning and GNSS data streams
JP7464790B2 (ja) 2021-03-24 2024-04-09 三菱電機株式会社 測位装置、測位プログラムおよび測位方法
US11860287B2 (en) 2022-03-01 2024-01-02 Swift Navigation, Inc. System and method for detecting outliers in GNSS observations
US11906640B2 (en) 2022-03-01 2024-02-20 Swift Navigation, Inc. System and method for fusing sensor and satellite measurements for positioning determination
US12019163B2 (en) 2022-09-12 2024-06-25 Swift Navigation, Inc. System and method for GNSS correction transmission

Similar Documents

Publication Publication Date Title
CN110727002A (zh) 一种基于稀疏正则化的单频单站动态gnss载波相位信号周跳修复方法
Watson et al. Robust navigation in GNSS degraded environment using graph optimization
Wang et al. GPS and GLONASS integration: modeling and ambiguity resolution issues
Nan et al. Integrated method for instantaneous ambiguity resolution using new generation GPS receivers
CN104714244A (zh) 一种基于抗差自适应Kalman滤波的多系统动态PPP解算方法
CN107966722B (zh) 一种gnss钟差解算方法
Watson et al. Evaluation of kinematic precise point positioning convergence with an incremental graph optimizer
CN109738926A (zh) 一种基于bp神经网络技术的gnss多路径效应改正方法
EP3430431A1 (en) Navigation satellite wide-lane bias determination and over-range adjustment system and method
CN110988948A (zh) 一种基于动对动相对定位场景中完好性分析方法
CN107422342A (zh) Gnss卫星钟差实时估计质量控制方法
Wang et al. Adaptive Kalman filtering based on posteriori variance-covariance components estimation
Lehmann Observation error model selection by information criteria vs. normality testing
Sarychev et al. Optimal regressors search subjected to vector autoregression of unevenly spaced TLE series
CN110703284B (zh) 一种基于稀疏核学习的单站gnss瞬时速度和加速度构建方法
Hu et al. Cycle slip detection and repair using an array of receivers with known geometry for RTK positioning
Joerger et al. Sequential residual-based RAIM
Wen et al. Mitigation of multiple outliers using consistency checking for GNSS standard point positioning in urban areas
Cellmer Single-epoch precise positioning using modified ambiguity function approach
Cellmer On-the-fly ambiguity resolution using an estimator of the modified ambiguity covariance matrix for the GNSS positioning model based on phase data
Murata et al. Experimental evaluation of dynamics-free PPP in space using Japanese ALOS2 dataset
Baroni et al. Evaluation of two integer ambiguity resolution methods for real time GPS positioning
Do-yoon et al. Development & performance analysis of Korean WADGPS positioning algorithm
Biswas et al. Application of a fast unscented Kalman Filtering method to satellite position estimation using a space-borne multi-GNSS receiver
Petovello et al. Kalman filter reliability analysis using different update strategies

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
AD01 Patent right deemed abandoned
AD01 Patent right deemed abandoned

Effective date of abandoning: 20231229