CN106526634A - 一种基于自调整卡尔曼滤波的多普勒辅助载波相位平滑伪距方法 - Google Patents

Abstract

本发明公开了一种基于自调整卡尔曼滤波的多普勒辅助载波相位平滑伪距方法，属于全球卫星导航系统伪距单点定位领域。本发明方案利用多普勒频移对码伪距的初值进行估计，并结合载波相位完成周跳检测，以决定卡尔曼滤波器的平滑方式以及过程噪声方差的大小；之后利用信噪比测量值的变化情况实现了卡尔曼滤波器的观测噪声方差的自调整，最终完成对码伪距的平滑。本发明方案解决了码伪距精度低的问题，消除了周跳对载波相位平滑伪距的影响，能够在复杂环境下得到精度较高的伪距，进而能够提高单点定位精度，且不仅限于静态定位场合，具有广阔的应用前景。

Description

一种基于自调整卡尔曼滤波的多普勒辅助载波相位平滑伪距 方法

技术领域

本发明涉及全球卫星导航系统伪距单点定位领域，特别是一种基于自调整卡尔曼滤波的多普勒辅助载波相位平滑伪距方法。

背景技术

从全球导航定位系统(Global Navigation Satellite System,GNSS)接收机的基带信号中可以获得码伪距、载波相位以及多普勒频移观测量。码伪距在过去一直被视为GNSS接收机最主要的基本距离测量值，在同一时刻利用至少4颗不同可见卫星的伪距测量值，接收机就可以实现三维绝对定位、测速和定时；但由于其包含钟差、电离层延时等各种误差，观测噪声大，多路径干扰对其影响也更大，精度为m级。作为另一个基本距离测量值，载波相位测量值非常平滑，精度很高，达到cm甚至mm级，但会受到周跳的影响，存在整周模糊度问题，它与码伪距具有明显区别，又呈现互补特性。而多普勒频移为载波相位的时间差分值，精度较高，可达到cm级，且不存在周跳。

已有技术提出基于Hatch滤波的载波相位平滑码伪距方法，即用两次测量的载波相位差计算两次测量之间的伪距变化量，再用这个伪距变化量去平滑测量的码伪距。虽然该方法在工程上被广泛应用，但其存在如下问题：

一、如果平滑时使用的码伪距初始值误差较大，那么此后的平滑过程中将很难消除此误差，这样即使码伪距得到了平滑，精度也不高；

二、在“城市峡谷”等复杂的接收环境中，信号容易被遮挡，GNSS接收机载波跟踪环路容易产生信号失锁和相位失周等现象，导致载波相位测量值发生较大变化，从而使该方法的应用受到限制。

三、理论上随着观测历元的增加和平滑过程的推进，平滑后的码伪距精度逐渐接近载波相位的精度，但是由于电离层对码伪距和载波相位的影响是反向的，当平滑处理周期数逐渐增大或电离层发生扰动时，会发生单频伪距平滑值发散的情况。虽然现有技术提出了采用卡尔曼滤波器作为载波相位平滑码伪距的工具，但传统的卡尔曼滤波器的量测噪声方差为常值，在动态GNSS定位中，不能很好地适合环境的变化，从而也达不到理想的平滑效果。

发明内容

本发明的目的是针对现有技术的不足，提供一种基于自调整卡尔曼滤波的多普勒辅助载波相位平滑伪距方法。

本发明解决其技术问题所采用的技术方案，对于每一颗卫星，均包括以下步骤：

(1)设定一个窗口长度N，利用求平均的方法，计算第N-1个历元t_N-1的码伪距估计值

(1-1)接收机第一次接收到码伪距测量值ρ的历元时刻记作t₀；当到达第N-1个历元t_N-1时，通过t₀至t_N-1的多普勒频移测量值对相应时刻的码伪距测量值ρ进行修正；对于t₀至t_N-1，依次依照下式计算出t_N-1的码伪距估计值

其中：t_i为第i个历元，ρ(t_i)为第i个历元的码伪距测量值，D(t_i)为多普勒频移测量值，λ为卫星信号的载波波长；

(1-2)用下式对进行求平均，可以得到利用多普勒频移推算的t_N-1的码伪距值

(2)建立平滑码伪距的卡尔曼滤波模型：

状态方程：

观测方程：

其中：

估计参数为k时刻平滑后的码伪距，作为状态向量；观测值y_k为k时刻的码伪距测量值，y作为观测向量；所述k时刻自滤波开始起算，k作为第k步滤波运算；控制输入u_k,k-1由k-1时刻至k时刻的载波相位变化量或多普勒频移测量值决定，u作为输入向量；w、v为满足正态分布的白噪声，均值为0，w_k-1为k-1时刻的过程噪声向量，方差记为Q_k-1；v_k为k时刻的观测噪声向量，方差记为R_k；A为由k-1时刻至k时刻的状态转移矩阵；B为k-1时刻输入向量与状态向量之间的关系矩阵；C为k时刻观测向量与状态向量之间的关系矩阵；为滤波初始时刻的码伪距，即步骤(1)中的码伪距值 的初始方差P₀＝R₀；A、B、C均为常系数；上述向量及矩阵的维数均为1；

优选的，A、B、C的取值均为1。

(3)卡尔曼滤波过程：

在t_N-1时刻开始执行卡尔曼滤波，包括以下两个过程：

①预测过程

其中，为的一步预测值，P^～为的误差的协方差，P为的误差的协方差。

本发明方案中，上式可简化为：

②校正过程

其中：K为权重值。

本发明方案中，上式可简化为：

优选的，上述技术方案的步骤(2)中：

控制输入u_k,k-1的取值方法如下：

利用多普勒频移进行周跳检测，即若

则认为该时刻发生周跳，其中：φ_k为第k时刻的载波相位测量值，D_k为第k时刻的多普勒频移测量值，Δt为相邻两个时刻的时间间隔，CS_Thrd为设定的检测周跳的阈值。

若检测到该时刻发生周跳，此时使用多普勒频移来对码伪距进行平滑，即

若检测到该时刻未发生周跳，将使用测距精度最高的载波相位来对码伪距进行平滑，即

u_k,k-1＝λ(φ_k-φ_k-1)

优选的，考虑到载波相位和多普勒频移的测距精度的不同，根据周跳检测的结果，过程噪声方差Q_k-1设定方法如下：

优选的，考虑到多路径干扰对码伪距的精度影响更大，而接收机接收卫星信号的环境从良好到存在多路径干扰的过程中，卫星信噪比处于急剧的下降状态；且受到多路径干扰的过程中，卫星信噪比处于起伏波动的状态。因此，本方案利用卫星信噪比下降率DropRate(SNR)以及信噪比标准差Std(SNR)，来判断多路径干扰的存在，k时刻的DropRate(SNR)和Std(SNR)的计算方式分别如下：

其中，(SNR)_k为k时刻的信噪比测量值，单位dB；M为统计卫星信噪比下降率使用的连续时刻个数；

其中，M′为统计卫星信噪比标准差使用的连续时刻个数，μ为M′个时刻内所有信噪比测量值的平均值；

则k时刻的卡尔曼滤波器观测噪声的方差R_k的取值方法如下：

R_k＝R_Basic×{(1+k₁×max(DropRate_Thrd-DropRate(SNR)_k,0)+k₂×max(Std(SNR)_k-Std_Thrd,0))}

其中，R_Basic＝(1)²为码伪距的噪声方差的常规值；k₁、k₂分别为设定的卫星信噪比下降率和标准差的权重系数；DropRate_Thrd为接收机接收卫星信号的环境从良好到存在多路径干扰时的卫星信噪比下降率的阈值；Std_Thrd为受到多路径干扰的过程中，卫星信噪比处于起伏波动的状态时的卫星信噪比标准差的阈值。

有益效果：

本发明方案解决了码伪距精度低的问题，消除了周跳对载波相位平滑码伪距的影响，减轻了接收机噪声以及多路径干扰的影响，即使在复杂的环境下也能得到精度较高的伪距，从而提高单点定位精度，且不仅限于静态定位场合；另外，由于本发明方案中的卡尔曼滤波算法采用的是标量的运算形式，计算量较小，因此可以应用于实时单点定位。

附图说明

图1：本发明平滑码伪距算法框图

图2：不同算法的伪距精度对比图

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步的详细说明。应当理解，此处所描述的具体实施方式仅仅用以解释本发明，并不限定本发明的保护范围。

结合图1，给出本发明一种基于自调整卡尔曼滤波的多普勒辅助载波相位平滑伪距方法的最优实施例：使用自主搭建的基站采集数据进行测试，数据时长为24小时，历元间隔为1s，数据采集的时间为2016年5月17日。对于每一颗卫星，均包括以下步骤：

(1)设定一个窗口长度N，利用求平均的方法，计算第N-1个历元t_N-1的码伪距估计值下面以某颗卫星为例具体说明计算方法：

(1-1)接收机第一次接收到码伪距测量值ρ的历元时刻记作t₀；当到达第N-1个历元t_N-1时，通过t₀至t_N-1的多普勒频移测量值对相应时刻的码伪距测量值ρ进行修正；对于t₀至t_N-1，依次按照下式计算出t_N-1的码伪距估计值

其中：t_i为第i个历元，ρ(t_i)为第i个历元的码伪距测量值，D(t_i)为第i个历元的多普勒频移测量值，λ为卫星信号的载波波长(本实施例中，仅研究单频观测量，λ为L₁载波的波长，即λ＝0.1903m)；

本实施例中，设定窗口长度N＝8，对于t₀至t₇，接收机接收到的码伪距测量值ρ(t_i)和多普勒频移测量值D(t_i)以及代入上式解得的见表1；

表1

(1-2)用下式对进行求平均，可以得到利用多普勒频移推算的t₇的码伪距值

将(1-1)中数值代入上式，即可解得为21757917.464m；

(2)建立平滑码伪距的卡尔曼滤波模型：

状态方程：

观测方程：

其中：

估计参数为k时刻平滑后的码伪距，作为状态向量；观测值y_k为k时刻的码伪距测量值，y作为观测向量；所述k时刻自滤波开始起算，k作为第k步滤波运算；控制输入u_k,k-1由k-1时刻至k时刻的载波相位变化量或多普勒频移测量值决定，u作为输入向量；w、v为满足正态分布的白噪声，均值为0，w_k-1为k-1时刻的过程噪声向量，方差记为Q_k-1；v_k为k-1时刻的观测噪声向量，方差记为R_k-1；A为由k-1时刻至k时刻的状态转移矩阵；B为k-1时刻输入向量与状态向量之间的关系矩阵；C为k时刻观测向量与状态向量之间的关系矩阵；为滤波初始时刻的码伪距，即步骤(1)中的本实施例中已求得为21757917.464m；的初始方差P₀＝R₀；

本实施例中，上述向量及矩阵的维数均为1；且A、B、C均为常系数，其值均为1；

本实施例中，由步骤(1)可知，在t₇时刻开始执行卡尔曼滤波；这里以k＝100时刻为例，来说明卡尔曼滤波的实现：

技术方案的步骤(2)中控制输入u_100,99的取值方法如下：

利用多普勒频移进行周跳检测，即若

则认为该时刻发生周跳，其中：k＝100时刻的载波相位测量值φ₁₀₀为111122494.262周，k＝99时刻的载波相位测量值为111119979.727周，k＝100时刻的多普勒频移测量值D₁₀₀为-2511.758Hz，第k＝99时刻的多普勒频移测量值D₉₉为-2514.281Hz，两个时刻的时间间隔Δt为1s，设定的检测周跳的阈值CS_Thrd为1.1周，代入上式，结果如下：

依此，判断k＝100时刻发生周跳，将使用多普勒频移来对码伪距进行平滑，即

本实施例中，考虑到载波相位和多普勒频移的测距精度的不同，根据周跳检测的结果，，过程噪声方差Q₉₉设定方法如下：

因此，过程噪声方差Q₉₉设定为(0.05)²；

进一步的，本实施例考虑到多路径干扰对码伪距的精度影响更大，而接收机接收卫星信号的环境从良好到存在多路径干扰的过程中，卫星信噪比处于急剧的下降状态；且受到多路径干扰的过程中，卫星信噪比处于起伏波动的状态。因此，本方案利用卫星信噪比下降率DropRate(SNR)以及信噪比标准差Std(SNR)，来判断多路径干扰的存在，k时刻的DropRate(SNR)和Std(SNR)的计算方式分别如下：

其中，(SNR)_k为k时刻的信噪比测量值；M为统计卫星信噪比下降率使用的连续时刻个数(本实施例中为20)；

其中，M′为统计卫星信噪比标准差使用的连续时刻个数(本实施例中为20)，μ为M′个时刻内所有信噪比测量值的平均值；

在k＝100时刻，计算DropRate(SNR)₁₀₀和Std(SNR)₁₀₀所需的信噪比数据见表2；

表2：

k	81	82	83	84	85	86	87	88	89	90
											(SNR)k/dB	50	50	50	49	49	49	49	49	48	48
k	91	92	93	94	95	96	97	98	99	100
											(SNR)k/dB	48	48	48	47	47	47	47	47	47	47

解得

DropRate(SNR)₁₀₀为0.15；

Std(SNR)₁₀₀为1.078dB。

在本实施例中，k＝100时刻的卡尔曼滤波器观测噪声的方差R₁₀₀的取值方法如下：

R₁₀₀＝R_Basic×{(1+k₁×max(DropRate_Thrd-DropRate(SNR)₁₀₀,0)+k₂×max(Std(SNR)₁₀₀-Std_Thrd,0))}

其中，R_Basic＝(1)²为码伪距的基本噪声方差；k₁、k₂分别为设定的卫星信噪比下降率和标准差的权重系数(本实施例中分别为10/5)；DropRate_Thrd为接收机接收卫星信号的环境从良好到存在多路径干扰时的卫星信噪比下降率的阈值(本实施例中为0.1)；Std_Thrd为受到多路径干扰的过程中，卫星信噪比处于起伏波动的状态时的卫星信噪比标准差的阈值(本实施例中为2.0dB)；

代入DropRate(SNR)₁₀₀和Std(SNR)₁₀₀，解得R₁₀₀＝1.5。

(3)卡尔曼滤波过程：

本实施例中，在t₇时刻开始执行卡尔曼滤波，包括以下两个过程：

①预测过程

其中，为的一步预测值，P^～为的误差的协方差，P为的误差的协方差。

本发明方案中，上式可简化为：

②校正过程

其中：K为权重值。

本发明方案中，上式可简化为：

由此可见，只须提供P₀，Q₀，R₀，u_1,0，y₁(实施例中已举例给出各参数的获取方法)，即可展开上述卡尔曼滤波过程。

对本实施例所用数据进行算法仿真，截取其中的30个历元进行比较，结果如图2所示，其中X坐标表示历元数，Y坐标表示两个相邻时刻的伪距差(由于伪距的数值太大，故使用相邻时刻的伪距差进行比较平滑的性能)。可以看出，本发明所提的这种平滑伪距算法是十分有效的，其平滑精度优于传统载波相位平滑伪距算法，且平滑后的伪距精度明显优于波动程度较大的未经过平滑的码伪距；通过本发明方案的这种平滑伪距算法，能够得到精度较高的伪距，从而提高单点定位精度，且不仅限于静态定位场合。

以上所述实施例仅表达了本发明最优的一种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。

Claims (5)

1.一种基于自调整卡尔曼滤波的多普勒辅助载波相位平滑伪距方法，其特征在于对于每一颗卫星，均包括以下步骤：

(1)设定一个窗口长度N，利用求平均的方法，计算第N-1个历元t_N-1的码伪距估计值

(1-1)接收机第一次接收到码伪距测量值ρ的历元时刻记作t₀；当到达第N-1个历元t_N-1时，通过t₀至t_N-1的多普勒频移测量值对相应时刻的码伪距测量值ρ进行修正；对于t₀至t_N-1，依次依照下式计算出t_N-1的码伪距估计值

$\widehat{\mathrm{\ρ}}{\left(t_{N-1}\right)}_i=\mathrm{\ρ}\left(t_i\right)-\mathrm{\λ}{\∫}_{t_i}^{t_{N-1}}D\left(t\right)dt,0\≤i\≤N-1$

其中：t_i为第i个历元，ρ(t_i)为第i个历元的码伪距测量值，D(t_i)为多普勒频移测量值，λ为卫星信号的载波波长；

(1-2)用下式对进行求平均，可以得到利用多普勒频移推算的t_N-1的码伪距值

$\stackrel{\‾}{\mathrm{\ρ}}\left(t_{N-1}\right)=\frac{1}{N}\underset{i=0}{\overset{N-1}{\Σ}}\widehat{\mathrm{\ρ}}{\left(t_{N-1}\right)}_i$

(2)建立平滑码伪距的卡尔曼滤波模型：

状态方程：

观测方程：

其中：

估计参数为k时刻平滑后的码伪距，作为状态向量；观测值y_k为k时刻的码伪距测量值，y作为观测向量；所述k时刻自滤波开始起算，k作为第k步滤波运算；控制输入u_k,k-1由k-1时刻至k时刻的载波相位变化量或多普勒频移测量值决定，u作为输入向量；w、v为满足正态分布的白噪声，均值为0，w_k-1为k-1时刻的过程噪声向量，方差记为Q_k-1；v_k为k时刻的观测噪声向量，方差记为R_k；A为由k-1时刻至k时刻的状态转移矩阵；B为k-1时刻输入向量与状态向量之间的关系矩阵；C为k时刻观测向量与状态向量之间的关系矩阵；为滤波初始时刻的码伪距，即步骤(1)中的码伪距值 的初始方差P₀＝R₀；A、B、C均为常系数；上述向量及矩阵的维数均为1；

(3)卡尔曼滤波过程：

在t_N-1时刻开始执行卡尔曼滤波，包括以下两个过程：

①预测过程

${\hat{x}}_k^{~}=A{\hat{x}}_{k-1}+{\mathrm{Bu}}_{k,k-1}$

$P_k^{~}={\mathrm{AP}}_{k-1}{A}^T+Q_{k-1}$

其中，为的一步预测值，P^～为的误差的协方差，P为的误差的协方差；

②校正过程

$K_k=P_k^{~}{C}^T{({\mathrm{CP}}_k^{~}{C}^T+R_k)}^{-1}$

$P_k=(1-K_{k}C)P_k^{~}$

${\hat{x}}_k={\hat{x}}_k^{~}+K_k(y_k-C{\hat{x}}_k^{~})$

其中：K为权重值。

2.根据权利要求1所述的一种基于自调整卡尔曼滤波的多普勒辅助载波相位平滑伪距方法，其特征在于A、B、C的取值均为1。

3.根据权利要求1所述的一种基于自调整卡尔曼滤波的多普勒辅助载波相位平滑伪距方法，其特征在于步骤(2)中：

控制输入u_k,k-1的取值方法如下：

利用多普勒频移进行周跳检测，即若

$|({\mathrm{\φ}}_k-{\mathrm{\φ}}_{k-1})+(\frac{1}{2}\·\mathrm{\Δ}t\·(D_k+D_{k-1}\left)\right)|>{\mathrm{CS}}_{Thrd}$

则认为该时刻发生周跳，其中：φ_k为第k时刻的载波相位测量值，D_k为第k时刻的多普勒频移测量值，Δt为相邻两个时刻的时间间隔，CS_Thrd为设定的检测周跳的阈值；

若检测到该时刻发生周跳，此时使用多普勒频移来对码伪距进行平滑，即

$u_{k,k-1}=-\mathrm{\λ}(\frac{1}{2}\·\mathrm{\Δ}t\·(D_k+D_{k-1}\left)\right)$

若检测到该时刻未发生周跳，将使用载波相位来对码伪距进行平滑，即

u_k,k-1＝λ(φ_k-φ_k-1)。

4.根据权利要求3所述的一种基于自调整卡尔曼滤波的多普勒辅助载波相位平滑伪距方法，其特征在于根据周跳检测的结果，过程噪声方差Q_k-1设定方法如下：

$Q_{k-1}=\left\{\begin{array}{cc}{\left(0.01\right)}^2& |({\mathrm{\φ}}_k-{\mathrm{\φ}}_{k-1})+(\frac{1}{2}\·\mathrm{\Δ}t\·(D_k+D_{k-1}\left)\right)|\≤{\mathrm{CS}}_{Thrd}\\ {\left(0.05\right)}^2& |({\mathrm{\φ}}_k-{\mathrm{\φ}}_{k-1})+(\frac{1}{2}\·\mathrm{\Δ}t\·(D_k+D_{k-1}\left)\right)|>{\mathrm{CS}}_{Thrd}\end{array}\right..$

5.根据权利要求4所述的一种基于自调整卡尔曼滤波的多普勒辅助载波相位平滑伪距方法，其特征在于利用卫星信噪比下降率DropRate(SNR)以及信噪比标准差Std(SNR)，来判断多路径干扰的存在，k时刻的DropRate(SNR)和Std(SNR)的计算方式分别如下：

$DropRate{\left(SNR\right)}_k=\frac{{\left(SNR\right)}_k-{\left(SNR\right)}_{k-M+1}}{Mt},k\≥M-1$

其中，(SNR)_k为k时刻的信噪比测量值，单位dB；M为统计卫星信噪比下降率使用的连续时刻个数；

$Std{\left(SNR\right)}_k=\sqrt{\frac{1}{M^{\′}}\underset{i=k-M^{\′}+1}{\overset{k}{\Σ}}{({\left(SNR\right)}_i-\mathrm{\μ})}^2},k\≥M^{\′}-1$

其中，M′为统计卫星信噪比标准差使用的连续时刻个数，μ为M′个时刻内所有信噪比测量值的平均值；

则k时刻的卡尔曼滤波器观测噪声的方差R_k的取值方法如下：

R_k＝R_Basic×{(1+k₁×max(DropRate_Thrd-DropRate(SNR)_k,0)+k₂×max(Std(SNR)_k-Std_Thrd,0))}

其中，R_Basic＝(1)²为码伪距的噪声方差的常规值；k₁、k₂分别为设定的卫星信噪比下降率和标准差的权重系数；DropRate_Thrd为接收机接收卫星信号的环境从良好到存在多路径干扰时的卫星信噪比下降率的阈值；Std_Thrd为受到多路径干扰的过程中，卫星信噪比处于起伏波动的状态时的卫星信噪比标准差的阈值。