一种基于稀疏正则化的单频单站动态GNSS载波相位信号周跳
修复方法
技术领域
背景技术
[0002] GNSS可以提供载波相位和伪距两种重要的观测量，其中前者的
精度比后者高两个数量级，所以在很多高精度GNSS应用中载波相位观测是必不可少的。然而，载波相位存在整周模糊性，这使得载波相位的模糊度解算成为一个必要的
数据处理步骤。模糊解算的成功实施一般需要在一段时间段内载波相位信号是连续
跟踪的，即不存在周跳。然而现实中由于观测环境的不理想，周跳的存在是不可避免的。有鉴于此，周跳的探测和修复则成为GNSS高精度应用中的一个必要的预处理步骤。
发明内容
[0003] 发明目的：针对在单频、单站、动态的情况下，如何修复GNSS周跳的问题，本发明提出一种基于稀疏正则化的单频单站动态GNSS载波相位信号周跳修复方法。
[0004] 技术方案：为实现本发明的目的，本发明所采用的技术方案是：
[0005] 一种基于稀疏正则化的单频单站动态GNSS载波相位信号周跳修复方法，所述周跳修复方法具体包括如下步骤：
[0006] S1：采集任一频点上预设时间段内的GNSS载波相位观测量和伪距观测量；
[0007] S2：对所述GNSS载波相位观测量和伪距观测量分别进行历元间差分操作，
构建时间差分载波相位观测量和时间差分伪距观测量；
[0008] S3：在相应历元处，对所述时间差分载波相位观测量和时间差分伪距观测量进行无几何线性组合，构建无几何时间差分观测量；
[0009] S4：根据所述无几何时间差分观测量，逐历元进行大周跳修复；
[0010] S5：将所述大周跳修复后的无几何时间差分观测量基于稀疏正则化的周跳估计方法估计存在的周跳的浮点解，对所述浮点解进行四舍五入取整操作，获取周跳的整数解，进行周跳修复。
[0011] 进一步地讲，所述时间差分载波相位观测量和时间差分伪距观测量，具体为：
[0012]
[0013] 其中：γ＝1/λ，
[0014] 为时间差分载波相位观测量，dk为时间差分伪距观测量，λ为载波
波长，xk为第k历元处的几何量增量，yk为第k历元处的周跳，εk为第k历元处伪距的观测误差，εk-1为第k-1历元处伪距的观测误差，ek为第k历元处载波相位的观测误差，ek-1为第k-1历元处载波相位的观测误差。
[0015] 进一步地讲，所述观测误差具体为：
[0016]
[0017] 其中：
[0018] εk-1为第k-1历元处伪距的观测误差，ek-1为第k-1历元处载波相位的观测误差，ek为第k历元处载波相位的观测误差，εk为第k历元处伪距的观测误差。
[0019] 进一步地讲，所述所述无几何时间差分观测量的观测方程，具体为：
[0020]
[0021] 其中：γ＝1/λ，
[0022] wk为无几何时间差分观测量，λ为载波波长， 为时间差分载波相位观测量，dk为时间差分伪距观测量，yk为第k历元处的周跳，εk-1为第k-1历元处伪距的观测误差，ek-1为第k-1历元处载波相位的观测误差，ek为第k历元处载波相位的观测误差，εk为第k历元处伪距的观测误差，ηk为观测误差。
[0023] 进一步地讲，所述步骤S4中逐历元进行大周跳修复，具体如下：
[0024] S4.1：根据所述无几何时间差分观测量的观测方程，所述无几何时间差分观测量即为第k历元处周跳的直接观测量，在周跳修复成功率满足固定成功率要求时，对所述周跳的浮点解进行四舍五入取整操作，具体为：
[0025]
[0026] 其中： 为第k历元处周跳的浮点解进行四舍五入取整操作后的值，[]表示四舍五入取整算子，wk为无几何时间差分观测量；
[0027] S4.2：根据所述四舍五入取整操作后的周跳浮点解，对所述周跳进行修复，并将所述大周跳修复后的观测量表示为ξk，所述大周跳修复后的观测量对应的观测方程具体为：
[0028]
[0029] 其中：ξk为大周跳修复后的观测量，wk为无几何时间差分观测量， 为第k历元处周跳的浮点解进行四舍五入取整操作后的值，yk为第k历元处的周跳，ηk为观测误差，zk为待求解参数。
[0030] 进一步地讲，所述步骤S5中获取周跳的整数解，具体如下：
[0031] S5.1：根据所述大周跳修复后的观测量对应的观测方程，采用所有所述大周跳修复后的无几何时间差分观测量，构建向量形式的观测方程，具体为：
[0032] ξ＝z+η
[0033] 其中：ξ为所有ξk组成的观测向量，ξk为大周跳修复后的观测量，z为所有zk组成的未知数参数向量，zk为待求解参数，η为所有ηk组成的观测误差向量，ηk为观测误差；
[0034] S5.2：将所述观测误差向量对应的协方差矩阵标记为Qηη，所述协方差矩阵具体为：
[0035]
[0036] 其中： γ＝1/λ，
[0037] Qηη为观测误差向量对应的协方差矩阵， 为伪距的观测误差εk的方差， 为载波相位的观测误差ek的方差，λ为载波波长，ek为第k历元处载波相位的观测误差，εk为第k历元处伪距的观测误差，n表示历元的总个数；
[0038] S5.3：预设N个超参数，并从N个所述超参数中选出最优的超参数；
[0039] S5.4：将所述最优超参数对应的周跳浮点解作为最终的浮点解，并将所述最终浮点解进行四舍五入取整操作，获取所述周跳的整数解，完成周跳修复。
[0040] 进一步地讲，所述步骤S5.3中从N个超参数中选出最优的超参数，具体如下：
[0041] S5.3.1：预设N个超参数，并从N个所述超参数中任意选取一个超参数；
[0042] S5.3.2：根据所述选取的超参数，未知数参数向量通过FISTA
算法使稀疏正则化代价函数最小，获取稀疏正则化代价函数最小时对应的浮点解，所述稀疏正则化代价函数，具体为：
[0043]
[0044] 其中：f为稀疏正则化代价函数，ξ为所有ξk组成的观测向量，ξk为大周跳修复后的观测量，z为所有zk组成的未知数参数向量，zk为待求解参数，Qηη为观测误差向量对应的协方差矩阵，μ为超参数；
[0045] S5.3.3：将所述稀疏正则化代价函数最小时对应的浮点解代入广义交叉检核公式中，获取所述选取的超参数对应的GCV值，所述广义交叉检核公式具体为：
[0046]
[0047] 其中：GCV(μ)为超参数μ对应的GCV值，n表示历元的总个数，n0为 中非零元素的个数，ξ为所有ξk组成的观测向量，ξk为大周跳修复后的观测量， 为稀疏正则化代价函数最小时对应的浮点解，Qηη为观测误差向量对应的协方差矩阵；
[0048] S5.3.4：从预设的N个所述超参数中选出一个未被选择过的超参数，并根据所述重新选择的超参数，重复步骤S5.3.2-步骤S5.3.3，直至预设的N个所述超参数中的所有超参数均被选择过，并将每个所述超参数对应的GCV值进行比较，从中选出最小GCV值，所述最小GCV值对应的超参数即为最优超参数。
[0049] 进一步地讲，所述步骤S5.3.2中未知数参数向量通过FISTA算法使稀疏正则化代价函数最小，具体如下：
[0050] S5.3.2.1：将待求解参数进行初始化，具体为：
[0051]
[0052] 其中：z0为
迭代前待求解参数的初步估计值，θ1和s1均为第一次迭代运算引入的辅助变量，ξ为所有ξk组成的观测向量，ξk为大周跳修复后的观测量；
[0053] S5.3.2.2：将所述待求解参数的初步估计值、第一次迭代运算引入的辅助变量均进行迭代计算，直至算法收敛，确定出每次迭代过程中所述待求解参数的大小，具体为：
[0054]
[0055] 其中：
[0056] zk和zk-1为待求解参数， 为 的软
阈值函数，θk和sk为第k次迭代运算引入的辅助变量，λmax为Qηη的最大特征值，Qηη为观测误差向量对应的协方差矩阵，λmin为Qηη的最小特征值，ξ为所有ξk组成的观测向量，ξk为大周跳修复后的观测量，sk+1和θk+1为第k+1次迭代运算引入的辅助变量，μ为超参数。
[0057] 进一步地讲，所述软阈值函数 算子中 通过高效的三对
角矩阵算法TMA进行计算，令 则g＝Qηηf，同时依次执行如下公
式，具体为：
[0058]
[0059]
[0060]
[0061] 其中： γ＝1/λ
[0062] 和 均为辅助变量，fk为向量f的第k个元素，gk为向量g的第k个元素， 为伪距的观测误差εk的方差， 为载波相位的观测误差ek的方差，λ为载波波长，ek为第k历元处载波相位的观测误差，εk为第k历元处伪距的观测误差，n表示历元的总个数。
[0063] 进一步地讲，所述软阈值函数 的求解过程，具体如下：
[0064] 第一步：给定向量x，对给定的向量x进行软阈值操作，得到一个同维数的向量，所述同维数向量的第j个元素，具体为：
[0065]
[0066] 其中：
[0067] x表示软阈值函数中的自变量，(|xj|-tμ)+为
铰链函数，sgn(xj)为符号函数，μ为超参数，λmax为Qηη的最大特征值，Qηη为观测误差向量对应的协方差矩阵，λmin为Qηη的最小特征值，θk为第k次迭代运算引入的辅助变量，ξ为所有ξk组成的观测向量，ξk为大周跳修复后的观测量；
[0068] 第二步：对所述铰链函数进行求解，具体为：
[0069]
[0070] 其中：
[0071] (|xj|-tμ)+为铰链函数，μ为超参数，λmax为Qηη的最大特征值，Qηη为观测误差向量对应的协方差矩阵，λmin为Qηη的最小特征值，xj为x的第j个分量，x表示软阈值函数中的自变量；
[0072] 第三步：对所述符号函数进行求解，具体为：
[0073]
[0074] 其中：sgn(xj)为符号函数，xj为x的第j个分量，x表示软阈值函数中的自变量。
[0075] 有益效果：与
现有技术相比，本发明的技术方案具有以下有益技术效果：
[0076] 本发明的周跳修复方法适用于单频、单站、动态的场合，在采用FISTA数值迭代方法的同时，还充分挖掘了周跳稀疏性这一先验条件，并充分考虑观测量历元间的相关性，从而针对性地采用高效的三对角矩阵算法，进而无需肉眼观测以及其他人为残余，具有高度自动化。
附图说明
[0077] 图1是本发明周跳修复方法的流程示意图。
具体实施方式
[0078] 为使本发明
实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。其中，所描述的实施例是本发明一部分实施例，而不是全部的实施例。因此，以下对在附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的
选定实施例。
[0079] 实施例1
[0080] 参考图1，本实施例提供了一种基于稀疏正则化的单频单站动态GNSS载波相位信号周跳修复方法，其中单频是指只采用单频观测，即单一频点上的载波相位观测和伪距观测，该单一频点可以是原始观测中的某一频点，也可以是多频情形时多频信号的某一线性组合。单站是指不需要参考站的辅助，可适用于只有流动站的单站绝对
定位情况。动态是指适用于流动站处于动态的情形。步骤具体如下：
[0081] 步骤S1：采集任一频点上预设时间段内的GNSS载波相位观测量和伪距观测量。该GNSS可以是GPS、Beidou、Galileo等任一码分多址卫星
导航系统。
[0082] 其中任一频点可以是单频情况下的观测量，也可以是多频情况下多频观测的任意一种线性组合，譬如：宽巷或窄巷组合。
[0083] 同样地，预设时间段内可以是任意一颗卫星的连续观测时间，也可以是该连续观测时间中的一部分。
[0084] 步骤S2：构建时间差分载波相位观测量TDCP即将GNSS载波相位观测量在相邻历元间做差，并将其标记为 构建时间差分伪距观测量TDPR即将伪距观测量在相邻历元间做差，并将其标记为dk。在 和dk中k＝1,2,...,n，其中：k表示历元的序号，n表示历元的总个数。具体地讲，时间差分过程可以消去大部分历元间共模误差或时间慢变误差，譬如：电离层延迟误差和多路径误差。
[0085] 在本实施例中，时间差分载波相位观测量TDCP和时间差分伪距观测量TDPR，具体为：
[0086]
[0087] 其中：γ＝1/λ，
[0088] 为时间差分载波相位观测量，dk为时间差分伪距观测量，λ为载波波长，xk为第k历元处的几何量增量，yk为第k历元处的周跳，εk为第k历元处伪距的观测误差，εk-1为第k-1历元处伪距的观测误差，ek为第k历元处载波相位的观测误差，ek-1为第k-1历元处载波相位的观测误差。
[0089] 具体地讲，第k历元处的周跳yk即为本实施例中致
力于求解的未知数。且该未知数是否存在也是未知的的，当该未知数存在时，则该第k历元处的周跳yk为不为零的整数，当该未知数不存在时，则该第k历元处的周跳yk为零。
[0090] 步骤S3：在相应历元处，对步骤S2中时间差分载波相位观测量TDCP和时间差分伪距观测量TDPR进行无几何线性组合。无几何线性组合操作可以消除卫星-载体间的几何距离、
对流层延迟误差、卫星轨道误差等非弥散，即与
频率无关的变量。
[0091] 同时进行无几何线性组合处理后的时间差分载波相位观测量TDCP可以认为仅包含周跳和观测噪声两部分。其中观测噪声包含所有未建模、未消除的剩余误差。
[0092] 进行无几何线性组合处理后的时间差分载波相位观测量TDCP和时间差分伪距观测量TDPR，可以按照下式构建无几何时间差分观测量，从而无几何时间差分观测量的观测方程，具体为：
[0093]
[0094] 其中：γ＝1/λ，
[0095] wk为无几何时间差分观测量，λ为载波波长， 为时间差分载波相位观测量，dk为时间差分伪距观测量，yk为第k历元处的周跳，εk-1为第k-1历元处伪距的观测误差，ek-1为第k-1历元处载波相位的观测误差，ek为第k历元处载波相位的观测误差，εk为第k历元处伪距的观测误差，ηk为观测误差。
[0096] 具体地讲，观测误差ηk具体为：
[0097]
[0098] 其中：
[0099] εk-1为第k-1历元处伪距的观测误差，ek-1为第k-1历元处载波相位的观测误差，ek为第k历元处载波相位的观测误差，εk为第k历元处伪距的观测误差。
[0100] 步骤S4：根据时间差分载波相位观测量TDCP和时间差分伪距观测量TDPR构建的无几何时间差分观测量，逐历元进行大周跳修复。
[0101] 在逐历元进行大周跳修复时，需要根据信号波长或频率以及观测量噪声特性，在考虑成功率的前提下，进行大周跳修复，即确定周跳修复成功率大于等于给定值时才可进行修复，否则不进行修复。在本实施例中，给定值选择为99％。
[0102] 这是由于：处理后的大周跳要么消失、要么变为小周跳，在后续步骤中，周跳越接近于1，则周跳向量的L1范数就越接近于L0范数。同时考虑修复成功率的目的是为了尽可能避免因为错误修复而人为引入周跳，过多地人为引入周跳将破坏周跳稀疏性，并致使步骤S5无法进行。具体过程如下：
[0103] 步骤S4.1：根据无几何时间差分观测量的观测方程，可以发现，无几何时间差分观测量wk即为要求解的未知数第k历元处的周跳yk的直接观测量，从而可以直接令第k历元处的周跳yk等于无几何时间差分观测量wk，进而获取该周跳的浮点解。当确定周跳修复成功率满足固定成功率要求时，对该周跳的浮点解进行四舍五入取整操作，具体为：
[0104]
[0105] 其中： 为第k历元处周跳的浮点解进行四舍五入取整操作后的值，[]表示四舍五入取整算子，wk为无几何时间差分观测量。
[0106] 在本实施例中，关于固定成功率要求，譬如：对于GPS L1频点，上述过程中浮点周跳精度则约为5，从而对大于等于5的周跳进行四舍五入取整操作。
[0107] 步骤S4.2：对第k历元处的周跳yk进行修复，令大周跳修复后的观测量表示为ξk，则大周跳修复后的观测量对应的观测方程，具体为：
[0108]
[0109] 其中：ξk为大周跳修复后的观测量，wk为无几何时间差分观测量， 为第k历元处周跳的浮点解进行四舍五入取整操作后的值，yk为第k历元处的周跳，ηk为观测误差，zk为待求解参数。
[0110] 具体地讲，待求解参数zk为下一个需要求解的新周跳。从而一旦获取了待求解参数zk的大小，即可得到原始周跳yk的大小。
[0111] 步骤S5：将大周跳修复后的无几何时间差分观测量ξk基于稀疏正则化的周跳估计方法估计存在的周跳的浮点解，并对该浮点解进行四舍五入取整操作，获取周跳的整数解，从而进行周跳修复。具体如下：
[0112] 步骤S5.1：根据步骤S4.2中大周跳修复后的观测量对应的观测方程，采用所有大周跳修复后的无几何时间差分观测量ξk，构建如下向量形式的观测方程，具体为：
[0113] ξ＝z+η
[0114] 其中：ξ为所有ξk组成的观测向量，ξk为大周跳修复后的观测量，z为所有zk组成的未知数参数向量，zk为待求解参数，η为所有ηk组成的观测误差向量，ηk为观测误差。
[0115] 步骤S5.2：将所有ηk组成的观测误差向量η对应的协方差矩阵标记为Qηη，该协方差矩阵标记为Qηη具体为：
[0116]
[0117] 其中： γ＝1/λ，
[0118] Qηη为观测误差向量对应的协方差矩阵， 为伪距的观测误差εk的方差， 为载波相位的观测误差ek的方差，λ为载波波长，ek为第k历元处载波相位的观测误差，εk为第k历元处伪距的观测误差，n表示历元的总个数。
[0119] 具体地讲，当伪距的观测误差εk的方差 载波相位的观测误差ek的方差 均为常数时，则观测误差向量对应的协方差矩阵Qηη为对称三对角Toplitz矩阵。
[0120] 步骤S5.3：预设N个超参数μ，并根据如下步骤从N个超参数μ中选出最优的超参数μ。具体如下：
[0121] 步骤S5.3.1：预设N个超参数μ，并从N个超参数μ中任意选取一个超参数μ。
[0122] 步骤S5.3.2：根据选出的超参数μ，所有zk组成的未知数参数向量z通过FISTA算法使如下稀疏正则化代价函数最小，获取稀疏正则化代价函数最小时对应的浮点解 其中以修复大周跳后的时间差分载波相位和时间差分伪距的无几何组合为观测量，以可能存在的周跳为待估参数，采用L1范数正则化法进行参数的稀疏求解，即采用加权残差平方和与参数L1范数之线性组合作为代价函数，组合系数又称为正则化参数，从而待估参数的求解过程即为此代价函数的最小化过程。
[0123] 稀疏正则化代价函数，具体为：
[0124]
[0125] 其中：f为稀疏正则化代价函数，ξ为所有ξk组成的观测向量，ξk为大周跳修复后的观测量，z为所有zk组成的未知数参数向量，zk为待求解参数，Qηη为观测误差向量对应的协方差矩阵，μ为超参数。
[0126] 在构建稀疏正则化代价函数时，其中的加权残差平方和构建过程中的协方差矩阵考虑了观测量相邻历元间的相关性，该相关性来自于步骤S2中相邻历元间的差分操作，从而其中的协方差矩阵为三对角对称矩阵。
[0127] 在本实施例中，所有zk组成的未知数参数向量z通过FISTA算法，使上述稀疏正则化代价函数最小，具体如下：
[0128] 步骤S5.3.2.1：将待求解参数zk进行初始化，具体为：
[0129]
[0130] 其中：z0为迭代前待求解参数的初步估计值，θ1和s1均为第一次迭代运算引入的辅助变量，ξ为所有ξk组成的观测向量，ξk为大周跳修复后的观测量。
[0131] 步骤S5.3.2.2：将待求解参数的初步估计值z0、第一次迭代运算引入的辅助变量θ1和s1均进行迭代计算，直至算法收敛，确定出每次迭代过程中待求解参数的大小，具体为：
[0132]
[0133] 其中：
[0134] zk和zk-1为待求解参数， 为 的软阈值函数，θk和sk为第k次迭代运算引入的辅助变量，λmax为Qηη的最大特征值，Qηη为观测误差向量对应的协方差矩阵，λmin为Qηη的最小特征值，ξ为所有ξk组成的观测向量，ξk为大周跳修复后的观测量，sk+1和θk+1为第k+1次迭代运算引入的辅助变量，μ为超参数。
[0135] 在本实施例中，算法收敛即为相邻两次迭代获取的待求解参数的之间差值的绝对值小于预设阈值。
[0136] 具体地讲， 中 通过高效的三对角矩阵算法TMA进行计算，令 则g＝Qηηf。从而 的求解即为求解方程g＝Qηηf，求
解方程g＝Qηηf需要依次执行如下公式，具体为：
[0137]
[0138]
[0139] 其中： γ＝1/λ
[0140] 和 均为辅助变量，fk为向量f的第k个元素，gk为向量g的第k个元素， 为伪距的观测误差εk的方差， 为载波相位的观测误差ek的方差，λ为载波波长，ek为第k历元处载波相位的观测误差，εk为第k历元处伪距的观测误差，n表示历元的总个数。
[0141] 在本实施例中，软阈值函数 的求解过程，具体如下：
[0142] 第一步：给定向量x，对给定的向量x进行软阈值操作，得到一个同维数的向量，该同维数向量的第j个元素，具体为：
[0143]
[0144] 其中：
[0145] x表示软阈值函数中的自变量，(|xj|-tμ)+为铰链函数，sgn(xj)为符号函数，μ为超参数，λmax为Qηη的最大特征值，Qηη为观测误差向量对应的协方差矩阵，λmin为Qηη的最小特征值，θk为第k次迭代运算引入的辅助变量，ξ为所有ξk组成的观测向量，ξk为大周跳修复后的观测量。
[0146] 第二步：对铰链函数(|xj|-tμ)+进行求解，具体为：
[0147]
[0148] 其中：
[0149] (|xj|-tμ)+为铰链函数，μ为超参数，λmax为Qηη的最大特征值，Qηη为观测误差向量对应的协方差矩阵，λmin为Qηη的最小特征值，xj为x的第j个分量，x表示软阈值函数中的自变量。
[0150] 第三步：对符号函数sgn(xj)进行求解，具体为：
[0151]
[0152] 其中：sgn(xj)为符号函数，xj为x的第j个分量，x表示软阈值函数中的自变量。
[0153] 步骤S5.3.3：将稀疏正则化代价函数最小时对应的浮点解 代入广义交叉检核公式中，获取步骤S5.3.1中任意选取的超参数μ对应的GCV值，该广义交叉检核公式具体为：
[0154]
[0155] 其中：GCV(μ)为超参数μ对应的GCV值，n表示历元的总个数，n0为 中非零元素的个数，ξ为所有ξk组成的观测向量，ξk为大周跳修复后的观测量， 为稀疏正则化代价函数最小时对应的浮点解，Qηη为观测误差向量对应的协方差矩阵。
[0156] 步骤S5.3.4：从预设的N个超参数μ中选出一个未被选择过的超参数μ，并根据该重新选择的超参数μ，重复步骤S5.3.2-步骤S5.3.3，直至预设的N个超参数μ中的所有超参数均被选择过。
[0157] 同时将每个超参数μ对应的GCV值进行比较，从中选出最小GCV值，并将该最小GCV值对应的超参数μ作为最优超参数。
[0158] 步骤S5.4：根据选出的最优超参数μ，将该最优超参数μ对应的周跳浮点解作为最终的浮点解，同时将该最终浮点解进行四舍五入取整操作获取该周跳的整数解，从而完成周跳修复。
[0159] 以上示意性的对本发明及其实施方式进行了描述，该描述没有限制性，附图中所示的也只是本发明的实施方式之一，实际的结构和方法并不局限于此。所以，如果本领域的普通技术人员受其启示，在不脱离本发明创造宗旨的情况下，不经创造性的设计出与该技术方案相似的结构方式及实施例，均属于本发明的保护范围。
