一种北斗卫星导航系统载波相位周跳探测与修复方法与流程

本发明涉及一种北斗卫星导航系统载波相位周跳探测与修复方法，属于北斗定位精度控制技术领域。

背景技术：

北斗卫星导航系统载波相位观测值中1周的周跳即可对定位结果造成分米级的偏差，因此周跳的探测与修复在全球卫星导航系统高质量定位中至关重要。由于北斗单频观测数据不能构成周跳探测精度较高的mw组合及电离层残差组合，其微小周跳难以探测与修复，因此有效探测与修复北斗单频观测数据中1周左右的微小周跳成为了周跳探测与修复领域的一个研究热点。许多学者将周跳信号看作信号中的奇异值，采用emd方法进行单频周跳探测与修复，但emd方法存在过包络、欠包络、模态混叠及端点效应等问题，会影响周跳探测精度。与emd方法相比，itd方法在分解速度上有着明显优势，通过对itd进行改进可以削弱端点效应的影响。

技术实现要素：

本发明的目的在于提供一种北斗卫星导航系统载波相位周跳探测与修复方法，用于对载波相位中发生的小周跳信号进行探测与修复。

本发明的技术方案是：本发明首先利用北斗载波相位观测量与伪距观测量构造周跳信号，然后用改进的固有时间尺度分解方法即改进的itd方法对周跳信号进行分解，得到若干个相互独立的固有旋转分量信号即pr分量，筛选出含有周跳的pr分量，最后对含有周跳的pr分量信号进行hilbert谱分析，探测出周跳发生的历元；接着对含有周跳的分量信号进行修复，将含有周跳的pr分量信号作为训练样本，采用粒子群优化算法对最小二乘支持向量机即ls-svm进行参数优化，利用优化后的ls-svm进行回归预测，通过计算实测值与预测值的差值来确定周跳的大小并修复周跳。

本发明的具体步骤如下：

step1、北斗卫星导航系统的原始观测量由伪距、载波相位组成，其观测方程表示如下所示：

式中：为卫星e至接收机u的几何距离；c为光速；dtu、δt^e分别为接收机钟差、卫星钟差；分别为对流层延迟、电离层延迟；分别为伪距多路径误差、载波相位多路径误差；分别为伪距观测噪声、载波相位观测噪声；为整周模糊度与波长的乘积；f＝(1,2,3)表示北斗卫星导航系统的三个频率；

公式(1)减去(2)可以消除卫星到接收机的距离、接收机的时钟差、卫星的时钟差、对流层延迟量，省略角标，得到公式(3)，如下所示：

再经过历元间一次差分消除整周模糊度得到公式(4)，如下所示：

式中：t为历元，δi为电离层延迟差分量，δe为观测噪声差分量，δm为多路径差分量；

若t+1历元发生周跳，公式(4)可以表示成公式(5)，如下所示：

式中，w表示周跳，单频周跳检测量表示为公式(6)，如下所示：

由于公式(6)中只包含电离层延迟差分量、观测噪声差分量及多路径差分量，将公式(6)等号后面的公式作为周跳检测量来探测周跳；

step2、改进的itd方法对周跳信号进行分解，步骤如下：

step2.1、确定信号xm所有局部极值点xk及其对应的时刻τk，其中，m为信号的数据点，k∈[1,2,…,m]，m为极值点总数，由于τ0不是极值点，定义τ0＝0，在连续的极值点间隔[τk,τk+1]上定义分段线性基线提取算子l，如公式(7)、(8)所示：

其中，α是线性缩放因子，用于控制固有旋转分量的幅值；lm为被提取算子l提取出来的基线信号。

通过公式(8)求出l2，，，lm-1，用三次b样条插值函数拟合代替公式(7)拟合l2，，，lm-1，得到拟合后的l2，，，lm-1；

step2.2、由于公式(8)中lk+1的值是从l2到lm-1，应用于周跳探测时，若周跳发生在端点处，将淹没周跳，导致无法探测出周跳，造成端点效应，故要估计端点l1和lm的值，在左端点处取拟合后的值l2,l3,l4，采用多项式拟合法拟合得到l1，采用同样的方法，对右端点拟合得lm，得到了拟合后的l1,l2，，，lm-1值；

step2.3、利用step2.2拟合后的l1,l2，，，lm-1值和公式h1(m)＝xm-lm计算h1(m)，其中h1(m)为固有旋转分量；

step2.4、将基线信号lm作为下一次分解的输入信号，重复步骤step2.1～step2.3，分解终止条件为基线信号lk为单调时，最终得到若干个pr分量信号和一个具有单调性的函数；

step3、从步骤2.4得到的若干个pr分量信号中筛选出含有周跳的分量信号，步骤如下：

step3.1、先求取周跳检测量的标准差s；

step3.2、计算各个分量信号的最大振幅a；

step3.3、若a大于2s则认为该分量信号包含周跳，保留该分量信号进行hilbert谱分析，根据hilbert谱中模极大值点出现的位置可以准确探测出周跳发生的历元；

step4、利用最小二乘支持向量机即ls-svm修复周跳信号，具体步骤如下：

step4.1、选取径向基函数作为ls-svm的核函数，核函数如公式(9)所示：

k(u,g)＝exp(-||u-g||²/2σ²)(9)

其中：||u-g||²表示空间中任意一点u到某一中心点g的欧式距离，σ＞0为高斯核函数的带宽；

ls-svm回归预测的优化目标为：

其中，c为惩罚参数，φ为不敏感损失函数，d为划分超平面的法向量，b为位移项，f(xv)-yv为回归预测值与真实值之差，v∈[1,2,…,o]表示采样点，o为正整数；

step4.2、采用粒子群优化算法对ls-svm的惩罚参数c以及核函数参数σ进行参数寻优；

step4.3、选择训练样本：对step3筛选出的j个pr分量信号进行信号截取，截取范围为信号的起始历元到周跳发生的前一个历元，截取的信号表示为qj＝{oj(i)}，i＝1,2,…,z-1，z是周跳发生时的历元，i表示历元z之前的历元，利用qj构建训练样本集tqj＝{xj(i),yj(i)},i＝1,2,…,z-n，n为数据滑动窗口，并对该训练样本集进行训练，其中训练输入为xj(i)＝[oj(i),oj(i+1),…,oj(i+n-1)]，输出为yj(i)＝oj(i+n)，得到训练模型；

step4.4、对于z历元处，将[oj(z-n),oj(z-n+1),…,oj(z-1)]作为预测模型的输入，其输出值为在z历元处的预测值将在历元z处的预测值及实测值oj(z)带入公式(9)，求出nj，nj表示发生在不同的pr分量信号上的周跳：

将nj值累加求和并取整即可得到历元z处的周跳值n，将z历元之后每个历元上的载波相位观测量均减去n值，即完成周跳的修复。

所述step2.1中α的取值为0.5。

所述step4.2中采用粒子群优化算法对ls-svm的惩罚参数c以及核函数参数σ进行参数寻优，其中参数及其初始值分别为：

粒子数r，r取30。

最大迭代次数t，t取300。

惯性权重系数r，r取0.5。

加速常数c1和c2，c1和c2均取2。

所述step4.3中n为数据滑动窗口，n取10。

所述step2.4中基线信号lk为单调时，其中单调表现为单调递增或单调递减。

本发明的有益效果是：

(1)本发明通过itd分解周跳信号可以克服emd方法存在的过包络、欠包络、模态混叠及端点效应等问题，而且加快分解速度。

(2)本发明通过对pr分量信号利用ls-svm预测，通过比较实测值与预测值的大小，可以同时修复周跳。

附图说明

图1为本发明的方法流程图；

图2为本发明中的周跳信号示意图；

图3为pr分量信号示意图；

图4为pr分量信号的hilbert谱示意图。

具体实施方式

实施例1：如图1-4所示，一种北斗卫星导航系统载波相位周跳探测与修复方法方法，首先利用北斗载波相位观测量与伪距观测量构造周跳信号，然后用改进的固有时间尺度分解方法即改进的itd方法对周跳信号进行分解，得到若干个相互独立的固有旋转分量信号即pr分量，筛选出含有周跳的pr分量，最后对含有周跳的pr分量信号进行hilbert谱分析，探测出周跳发生的历元；接着对含有周跳的分量信号进行修复，将含有周跳的pr分量信号作为训练样本，采用粒子群优化算法对最小二乘支持向量机即ls-svm进行参数优化，利用优化后的ls-svm进行回归预测，通过计算实测值与预测值的差值来确定周跳的大小并修复周跳。具体步骤如下：

step1、首先利用北斗星通ur370-cors三系统七频接收机在云南省昆明市昆明理工大学信自学院顶楼采集北斗载波相位观测量与伪距观测量，为了接收到较多的卫星信号，同时降低多路径误差，设置卫星截止高度角为10度，实验时接收机采样频率为1hz。选用北斗b3频段的伪距观测量作为公式(1)表示的数据，选择北斗b3频段的载波相位观测量作为公式(2)表示的数据，如下所示：

式中：为卫星e至接收机u的几何距离；c为光速；dtu、δt^e分别为接收机钟差、卫星钟差；分别为对流层延迟、电离层延迟；分别为伪距多路径误差、载波相位多路径误差；分别为伪距观测噪声、载波相位观测噪声；为整周模糊度与波长的乘积；f＝(1,2,3)表示北斗卫星导航系统的三个频率；

公式(1)减去(2)可以消除卫星到接收机的距离、接收机的时钟差、卫星的时钟差、对流层延迟量，省略角标，得到公式(3)，如下所示：

再经过历元间一次差分消除整周模糊度得到公式(4)，如下所示：

式中：t为历元，δi为电离层延迟差分量，δe为观测噪声差分量，δm为多路径差分量；

若t+1历元发生周跳，公式(4)可以表示成公式(5)，如下所示：

式中，w表示周跳，单频周跳检测量表示为公式(6)，如下所示：

由于公式(6)中只包含电离层延迟差分量、观测噪声差分量及多路径差分量，将公式(6)等号后面的公式作为周跳检测量来探测周跳，利用公式(6)构造周跳检测量，周跳检测量如图2；

步骤2、利用改进的itd方法分解周跳检测量；

步骤2.1、求出周跳检测量的所有局部极值点及对应的时刻，通过计算公式(8)，求出l2，，，lm-1，用三次b样条插值函数拟合代替公式(7)拟合l2，，，lm-1，得到拟合后的l2，，，lm-1，公式(7)、公式(8)如下所示；

其中，α是线性缩放因子，用于控制固有旋转分量的幅值，取值为0.5；lm为被提取算子l提取出来的基线信号；

步骤2.2、估计端点l1和lm的值，在左端点处取拟合后的值l2,l3,l4，采用多项式拟合法拟合得到l1，采用同样的方法，对右端点拟合得lm，得到了拟合后的l1,l2，，，lm-1值；

步骤2.3、利用步骤2.2拟合后的l1,l2，，，lm-1值和公式h1(m)＝xm-lm计算h1(m)，其中h1(m)为固有旋转分量；

步骤2.4、将基线信号lm作为下一次分解的输入信号，重复步骤2.1～2.3，分解终止条件为基线信号lk为单调时，最终得到3个pr分量信号和一个具有单调性的函数，分量信号如图3所示：

步骤3、从步骤2.4得到的3个pr分量信号中筛选出含有周跳的分量信号，步骤如下：

步骤3.1、先求取周跳检测量的标准差s；

步骤3.2、计算3个分量信号的最大振幅a；

步骤3.3、3个分量信号中的前两个分量信号的a值大于2s，认为前两个分量信号包含周跳，保留前两个分量信号进行hilbert谱分析，hilbert谱如图4所示，根据图4中hilbert谱中模极大值点出现的位置可以准确探测出：周跳发生在100历元、450历元；

步骤4、对前两个pr分量信号利用ls-svm预测，通过比较实测值与预测值的大小来修复周跳。

步骤4.1、选取径向基函数作为核函数，核函数如公式(9)所示：

k(u,g)＝exp(-||u-g||²/2σ²)(9)

其中：||u-g||²表示空间中任意一点u到某一中心点g的欧式距离，σ＞0为高斯核函数的带宽；

ls-svm回归预测的优化目标为：

其中，c为惩罚参数，φ为不敏感损失函数，d为划分超平面的法向量，b为位移项，f(xv)-yv为回归预测值与真实值之差，v∈[1,2,…,o]表示采样点，o为正整数；

步骤4.2、采用粒子群优化算法对ls-svm的惩罚参数c以及核函数参数σ进行参数寻优，其中参数及其初始值分别为：

粒子数r，r取30。

最大迭代次数t，t取300。

惯性权重系数r，r取0.5。

加速常数c1和c2，c1和c2均取2。

步骤4.3、选择训练样本：对步骤3筛选出的2个pr分量信号进行信号截取，在历元区间[1,99]与[351,449]之间截取信号。修复第100历元的周跳时，依次选择[1,9],[2-10],…,[90,98]历元的数据作为训练模型的输入，第10,11,…,99历元的数据作为训练模型的输出，来训练模型得到预测模型1。修复第450历元的周跳时，依次选择[351,359],[352,360],…,[440,448]历元的数据作为训练模型的输入，第369,361,…,449历元的数据作为训练模型的输出，来训练模型，得到预测模型2：

步骤4.4、对于z历元处，将[oj(z-n),oj(z-n+1),…,oj(z-1)]作为预测模型的输入，其输出值为在z历元处的预测值将在历元z处的预测值及实测值oj(z)带入公式(9)，求出nj，nj表示发生在不同的pr分量信号上的周跳：

将nj值累加求和并取整即可得到历元z处的周跳值n，将z历元之后每个历元上的载波相位观测量均减去n值，即完成周跳的修复，本实施例具体如下：

修复100历元的周跳时，将两个分量信号区间[91,99]历元数据作为预测模型1的输入，在100历元处第1个pr分量的预测值为0.1021，第2个pr分量的预测值为0.0363，相加得到最终预测值为0.1384。在100历元真实值减去预测值为1.7344，取整数为2，则在100历元发生2周周跳。修复450历元的周跳时，将两个分量信号区间[341,449]历元数据作为预测模型2的输入，在450历元处第1个pr分量的预测值为0.1372，第2个pr分量的预测值为0.0275，相加得到最终预测值为0.1647。在450历元真实值减去预测值为1.021，取整数为1，则在450历元发生1周周跳。在b3频率载波相位观测量的100、450历元及其以后历元的载波相位上分别减去1,2即完成了周跳的修复。

上面结合附图对本发明的具体实施例作了详细说明，但是本发明并不限于上述实施例，在本领域普通技术人员所具备的知识范围内，还可以在不脱离本发明宗旨的前提下作出各种变化。