CN102269819A - 一种用于精密单点定位方法的估计方法 - Google Patents

Abstract

本发明提供了一种用于精密单点定位方法的估计方法，采用kalman滤波估计方式，其特征在于：通过历元间差分计算建立提供速度观测量和钟漂观测量的观测方程，在kalman滤波的动态模型里引入速度状态产生和钟漂状态参数，根据所述观测方程和动态模型进行kalman递推估计。采用本发明所提出利用历元间差分信息进行精密单点定位收敛的方法，不光能加速收敛，对于提高解的精度也有很大的帮助。

Description

一种用于精密单点定位方法的估计方法

技术领域

本发明属于GNSS动态定位技术领域，特别是涉及一种利用历元间差分建立精确的动态状态模型，加速精密单点定位收敛及提高精度的方法。

背景技术

目前，GNSS高精度快速动态定位应用非常广泛。其中精密单点定位(PPP)技术由于只需要一台接收机进行观测即可精确确定观测点的位置，因此应用范围非常广泛。但是制约精密单点定位使用的一个最大的缺陷就是它不能很快的完成初始化，一般需要30分钟才能够使用收敛，这对于需要进行即刻定位的用户来说，是难以忍受的。

精密单点定位之所以要这么久的收敛时间，主要是由于在精密单点定位中，模糊度是不具有整周特性的，因此没有办法采用在传统的双差定位中的搜索的方法来确定模糊度，而只能采用滤波逐历元解算让其逐渐收敛到正确的模糊度值。

发明内容

针对上述问题，本发明提出一种能够加速精密单点定位收敛过程的估计方法。

本发明提供的技术方案是一种用于精密单点定位方法的估计方法，采用kalman滤波估计方式，通过历元间差分计算建立提供速度观测量和钟漂观测量的观测方程，在kalman滤波的动态模型里引入速度状态参数和钟漂状态参数，根据所述观测方程和动态模型进行kalman递推估计；

所述建立观测方程的过程如下，在没有周跳的情况下，同一颗卫星在相邻两个历元的无电离层组合的观测方程分别：

omc₁＝l₁ΔX₁+m₁ΔY₁+n₁ΔZ₁+trop₁+cdt₁+Nλ

omc₂＝l₂ΔX₂+m₂ΔY₂+n₂ΔZ₂+trop₂+cdt₂+Nλ

由于接收机在两个相邻历元间的运动位移很小，以上两个方程是按照同一个近似坐标进行线性化；其中，omc₁、omc₂为相位观测值减去计算值得到的结果，l₁、l₂、m₁、m₂、n₂、n₂为卫星在X、Y、Z三个方向的方向余弦，ΔX₁、ΔX₂、ΔY₁、ΔY₂、ΔZ₁、ΔZ₂为对近似坐标三个方向的改正数，trop₁、trop₂是对流层倾斜路径的延迟量，cdt₁、cdt₂是接收机钟差的距离表示，Nλ是整周模糊度的距离形式，下标1和2用于区别标识两个相邻历元；将以上两式相减得到：

${\mathrm{omc}}_1-{\mathrm{omc}}_2=\frac{l_1+l_2}{2}\mathrm{\Δx}+\frac{m_1+m_2}{2}\mathrm{\Δy}+\frac{n_1+n_2}{2}\mathrm{\Δz}+c({\mathrm{dt}}_1-{\mathrm{dt}}_2)$

当有R颗卫星时，得到的观测方程为：z_k＝H_kx_k+v_k，其中，z_k、H_k、x_k、v_k的含义分别为第k个历元的观测值、观测设计矩阵、状态参数和观测值残差；

其中，l1_k…lR_k、m1_k…mR_k、n1_k…nR_k分别为在第k个历元，第1到第R颗卫星的X、Y、Z三个方向的方向余弦；l1_k-1…lR_k-1、m1_k-1…mR_k-1、n1_k-1…nR_k-1分别为在第k-1个历元，第1到第R颗卫星的X、Y、Z三个方向的方向余弦；map1_k…mapR_k表示在第k个历元，第1到第R颗卫星的对流层投影函数；

其中，状态参数向量形式为：

$x_k=\left[\begin{array}{cccccccccccc}x& y& z& \stackrel{\·}{x}& \stackrel{\·}{y}& \stackrel{\·}{z}& \mathrm{dt}& d\stackrel{\·}{t}& \mathrm{ztd}& \mathrm{amb}1& ......& \mathrm{ambR}\end{array}\right]$

其中，xyz为位置参数，

为引入的速度状态参数，dt为接收机钟差参数，

为引入的接收机的钟漂状态参数，ztd为天顶对流层延迟，amb1…ambR分别为第1到第R颗卫星的模糊度参数；

所述在kalman滤波的动态模型里引入速度状态参数和钟漂状态参数，实现方式如下，kalman滤波的动态方程为：

x_k＝Φ_k，k-1x_k-1+w_k-1

其中，x_k、x_k-1分别为第k和第k-1个历元的状态参数，是未知数向量；Φ_k，k-1为第k-1到第k个历元的状态转移矩阵，w_k-1为过程噪声；引入了历元间差分之后的待求解的状态参数与上述观测模型一致：

$x_k=\left[\begin{array}{cccccccccccc}x& y& z& \stackrel{\·}{x}& \stackrel{\·}{y}& \stackrel{\·}{z}& \mathrm{dt}& d\stackrel{\·}{t}& \mathrm{ztd}& \mathrm{amb}1& ......& \mathrm{ambR}\end{array}\right]$

相应的状态转移矩阵为：

${\mathrm{\Φ}}_{k,k-1}=\left[\begin{array}{cccc}{I}_c& {\mathrm{\Δt}\·I}_c& 0& 0\\ 0& I_c& 0& 0\\ 0& 0& {\mathrm{\Φ}}_{\mathrm{Ck},k-1}& 0\\ 0& 0& 0& I_a\end{array}\right]$

其中，I_c为3×3单位阵，I_a为R×R单位阵；Φ_Ck，k-1为钟差的状态转移矩阵，形式为：

Δt是两个相邻历元的时间间隔；对应的过程噪声矩阵为：

$Q_{k,k-1}=E\[{\mathrm{ww}}^T\]=\left[\begin{array}{cccc}\left(S_P\frac{{\mathrm{\Δt}}^3}{3}\right)\·I_p& \left(S_P\frac{{\mathrm{\Δt}}^2}{2}\right)\·I_p& 0& 0\\ \left(S_P\frac{{\mathrm{\Δt}}^2}{2}\right)\·I_p& \left(S_P\mathrm{\Δt}\right)\·I_p& 0& 0\\ 0& 0& Q_c& 0\\ 0& 0& 0& \left(S_a\mathrm{\Δt}\right)I_a\end{array}\right]$

其中，

表示噪声的方差矩阵，S_P是位置参数的噪声因子，I_n是3×3单位阵，S_a是模糊度的噪声因子，I_a为R×R单位阵，

$Q_c=\left[\begin{array}{cc}{S}_{\mathrm{dt}}\mathrm{\Δt}+S_{d\stackrel{\·}{t}}\frac{{\mathrm{\Δt}}^3}{3}& S_{d\stackrel{\·}{t}}\frac{{\mathrm{\Δt}}^{\text{2}}}{2}\\ S_{d\stackrel{\·}{t}}\frac{\mathrm{\Δ}{t}^2}{2}& S_{d\stackrel{\·}{t}}\mathrm{\Δt}\end{array}\right]$

其中，S_dt是接收机钟漂的噪声因子。

本发明利用同一颗卫星两个相邻历元间的观测方程进行求差，消去对流层，将方向余弦取平均提取出来，得到速度、钟漂观测值；根据引入了速度观测值的动态系统的特性，确定系统的动态模型、状态转移矩阵和过程噪声矩阵等。因此，采用本发明所提出利用历元间差分信息加速PPP收敛的方法，不光能加速收敛，对于提高解的精度也有很大的帮助。

附图说明

图1为本发明实施例的流程图。

具体实施方式

以下根据实施例和附图说明本发明技术方案：

参见附图1，实施例的基本流程为：一方面，考察上一历元和当前新历元，通过历元间差分计算引入历元间差分方程建立提供速度观测量和钟漂观测量的观测方程。另一方面，根据解算信息，在kalman滤波的动态模型里引入速度状态产生和钟漂状态参数。最后根据所述观测方程和动态模型进行kalman递推估计，实现kalman滤波估计方式。

一、Kalman滤波

目前，在精密单点定位中使用的估计方式为kalman滤波，实施例也以该方式为基础。这种估计方法可以使用本历元之前所有历元的观测信息(序贯最小二乘)，同时通过一个动态模型描述两个历元之间的未知量的变化。在使用kalman滤波估计的时候，未知量，也即状态向量的动态特性描述的愈是精确，滤波计算就会愈好。

Kalman滤波的动态模型可以通过如下的动态方程来描述：

x_k＝Φ_k，k-1x_k-1+w_k，k-1，w_k，k-1～N(0，Q_k，k-1)

其中，x_k系统在k历元的状态；Φ_k k-1是系统状态转移矩阵，它把系统当前历元k-1的状态和下一个历元k的状态联系到了一块。w_k，k-1是系统的噪声向量，Q_k，k-1是系统的过程噪声矩阵，它描述了动态方程的不确定性，也即动态方程所没有描述到的部分。N(0，Q_k，k-1)表示w_k，k-1是一个均值为零，方差为Q_k，k-1的正态分布。

Kalman滤波的观测模型为：

z_k＝H_kx_k+v_k，v_k～N(0，R_k)

其中，z_k是观测向量，H_k是设计矩阵，x_k是系统的参数向量。v_k是系统的观测噪声向量，R_k是观测噪声矩阵。N(0，R_k)表示v_k是服从均值为零，方差为R_k的正态分布。

有了动态模型和观测方程，就可以进行kalman递推估计，该部分属于现有技术。为便于实施参考起见，本发明予以简述：其基本的计算过程可以归结为预测、滤波增益和滤波计算三步。

1)预测值的计算

首先，根据前一次滤波值(或初值)计算预测值：

${{\hat{x}}_{k,k-1}=\mathrm{\Φ}}_{k,k-1}{x}_{k-1}$ 根据前一次得到的滤波误差方差矩阵P_k，k-1以及系统的过程噪声矩阵Q_k，k-1计算预测的误差方差矩阵：

P_k，k-1＝Φ_k，k-1P_k，k-1Φ^T _k，k-1+Q_k，k-1

2)计算kalman滤波增益

滤波增益矩阵为：

$K_k=P_{k,k-1}{H}_k^T{\[H_k{P}_{k,k-1}{H}_k^T+R_k\]}^{-1}$

其中，H_k是设计矩阵，R_k是观测噪声矩阵，接下来，计算更新的观测值

$v_k=z_k-H_k{\hat{x}}_{k,k-1}$

3)滤波估计值计算

计算滤波估计：

x_k，k＝x_k，k-1+K_kv_k

计算滤波误差方差矩阵：

P_k，k＝[I-K_kH_k]P_k，k-1

二、历元间差分方法

本发明在以上kalman滤波的基础上，采用历元间差分的方法，在没有周条的情况下，可以消除模糊度参数，而且为系统提供了一个速度观测量。这个速度观测量可以为系统提供更加精确的动态模型描述。可以将普通精密单点定位的滤波过程与历元间差分结合，从而加强原来系统的动态模型的精度。实施例设定采用7颗卫星，即R＝7。

1)动态方程

引入了历元间差分的动态方程为：

x_k＝Φ_k，k-1x_k-1+w_k-1

x_k、x_k-1分别为第k和第k-1个历元的状态参数，是未知数向量；Φ_k，k-1为第k-1到第k个历元的状态转移矩阵，w_k-1为过程噪声。

 其中，

(因为实施例有七颗卫星)。对应为位置，速度，钟差，钟漂，对流层，模糊度。即x y z为位置参数，

为引入的速度状态参数，dt为接收机钟差参数，dt为引入的接收机的钟漂状态参数，ztd为天顶对流层延迟，amb1…amb7分别为第1到第7颗卫星的模糊度参数。

相应的状态转移矩阵为：

${\mathrm{\Φ}}_{k,k-1}=\left[\begin{array}{cccc}{I}_c& {\mathrm{\Δt}\·I}_c& 0& 0\\ 0& I_c& 0& 0\\ 0& 0& {\mathrm{\Φ}}_{\mathrm{Ck},k-1}& 0\\ 0& 0& 0& I_a\end{array}\right]$

其中，I_c为3×3单位阵，I_a为7×7单位阵。Φ_Ck，k-1为钟差的状态转移矩阵，形式为：

Δt是两个相邻历元的时间间隔。

对应的过程噪声矩阵为：

$Q_{k,k-1}=E\[{\mathrm{ww}}^T\]=\left[\begin{array}{cccc}\left(S_P\frac{{\mathrm{\Δt}}^3}{3}\right)\·I_p& \left(S_P\frac{{\mathrm{\Δt}}^2}{2}\right)\·I_p& 0& 0\\ \left(S_P\frac{{\mathrm{\Δt}}^2}{2}\right)\·I_p& \left(S_P\mathrm{\Δt}\right)\·I_p& 0& 0\\ 0& 0& Q_c& 0\\ 0& 0& 0& \left(S_a\mathrm{\Δt}\right)I_a\end{array}\right]$

其中，

表示噪声的方差矩阵，S_P是位置参数的噪声因子，I_p是3×3单位阵，S_a是模糊度的噪声因子，

$Q_c=\left[\begin{array}{cc}{S}_{\mathrm{dt}}\mathrm{\Δt}+S_{d\stackrel{\·}{t}}\frac{{\mathrm{\Δt}}^3}{3}& S_{d\stackrel{\·}{t}}\frac{{\mathrm{\Δt}}^{\text{2}}}{2}\\ S_{d\stackrel{\·}{t}}\frac{\mathrm{\Δ}{t}^2}{2}& S_{d\stackrel{\·}{t}}\mathrm{\Δt}\end{array}\right]$

S_dt是接收机钟漂的噪声因子。具体实施时，S_P要根据物体的运动特征设置一个合理的值。S_a设置为一个非常小的值，但是不为零，这样允许模糊度吸收一些未模型化的偏差，例如多路径等。

2)观测方程

加入了历元间差分信息的观测方程为：

z_k＝H_kx_k+v_k

其中z_k是观测值减去近似值，即OMC。H_k为观测设计矩阵，表达形式为：

上述矩阵为21×16的矩阵。其中第一至第七行是无电离层的伪距观测值对应的部分；八至第十四行为为电离层的相位对应的部分；第十五至第二十一行为历元间差分相位对应的部分。其中，l1_k…l7_k、m1_k…m7_k、n1_k…n7_k分别为在第k个历元，第1到第R颗卫星的X、Y、Z三个方向的方向余弦；l1_k-1…lR_k-1、m1_k-1…mR_k-1、n1_k-1…nR_k-1分别为在第k一1个历元，第1到第R颗卫星的X、Y、Z三个方向的方向余弦；map1_k…mapR_k表示在第k个历元，第1到第R颗卫星的对流层投影函数。

其中，第十五至第二十一行的方向余弦是本历元方向余弦与上一历元方向余弦的平均。原理如下：

同一颗卫星在两个历元的观测方程分别：

omc₁＝l₁ΔX₁+m₁ΔY₁+n₁ΔZ₁+trop₁+cdt₁+Nλ

omc₂＝l₂ΔX₂+m₂ΔY₂+n₂ΔZ₂+trop₂+cdt₂+Nλ

由于接收机在两个相邻历元间的运动位移很小，以上两个方程是按照同一个近似坐标进行线性化；其中，omc₁、omc₂为相位观测值减去计算值得到的结果，l₁、l₂、m₁、m₂、n₁、n₂为卫星在X、Y、Z三个方向的方向余弦，ΔX₁、ΔX₂、ΔY₁、ΔY₂、ΔZ₁、ΔZ₂为对近似坐标三个方向的改正数，trop₁、trop₂是对流层倾斜路径的延迟量，cdt₁、cdt₂是接收机钟差的距离表示，Nλ是整周模糊度的距离形式，下标1和2用于区别标识两个相邻历元。在两个相邻历元中，如果没有周跳，或者周跳被正确修复，上面两个式子中的模糊度参数是一样的。由于位置变化很小，因此对流层参数也是一样的。这两项在历元间求差的过程中可以消除。又由于在间隔很小的历元间(一般是小于1秒的采样间隔)，卫星位置的变化也很小，因此上面两个式子中的卫星方向余弦是非常接近的，可以把他们取个平均作为一个值。

基于以上分析，可以将上面两式求差得到：

${\mathrm{omc}}_1-{\mathrm{omc}}_2=\frac{l_1+l_2}{2}\mathrm{\Δx}+\frac{m_1+m_2}{2}\mathrm{\Δy}+\frac{n_1+n_2}{2}\mathrm{\Δz}+c({\mathrm{dt}}_1-{\mathrm{dt}}_2)$

根据以上推导得到的动态方程和观测方程，就可以进行求解计算。

Claims (1)

1.一种用于精密单点定位方法的估计方法，采用kalman滤波估计方式，其特征在于：通过历元间差分计算建立提供速度观测量和钟漂观测量的观测方程，在kalman滤波的动态模型里引入速度状态参数和钟漂状态参数，根据所述观测方程和动态模型进行kalman递推估计；

所述建立观测方程的过程如下，

在没有周跳的情况下，同一颗卫星在相邻两个历元的无电离层组合的观测方程分别：

omc₁＝l₁ΔX₁+m₁ΔY₁+n₁ΔZ₁+trop₁+cdt₁+Nλ

omc₂＝l₂ΔX₂+m₂ΔY₂+n₂ΔZ₂+trop₂+cdt₂+Nλ

由于接收机在两个相邻历元间的运动位移很小，以上两个方程是按照同一个近似坐标进行线性化；其中，omc₁、omc₂为相位观测值减去计算值得到的结果，l₁、l₂、m₁、m₂、n₁、n₂为卫星在X、Y、Z三个方向的方向余弦，ΔX₁、ΔX₂、ΔY₁、ΔY₂、ΔZ₁、ΔZ₂为对近似坐标三个方向的改正数，trop₁、trop₂是对流层倾斜路径的延迟量，cdt₁、cdt₂是接收机钟差的距离表示，Nλ是整周模糊度的距离形式，下标1和2用于区别标识两个相邻历元；

将以上两式相减得到：

${\mathrm{omc}}_1-{\mathrm{omc}}_2=\frac{l_1+l_2}{2}\mathrm{\Δx}+\frac{m_1+m_2}{2}\mathrm{\Δy}+\frac{n_1+n_2}{2}\mathrm{\Δz}+c({\mathrm{dt}}_1-{\mathrm{dt}}_2)$

当有R颗卫星时，得到的观测方程为：z_k＝H_kx_k+v_k，其中，z_k、H_k、x_k、v_k的含义分别为第k个历元的观测值、观测设计矩阵、状态参数和观测值残差；

其中，l1_k…lR_k、m1_k…mR_k、n1_k…nR_k分别为在第k个历元，第1到第R颗卫星的X、Y、Z三个方向的方向余弦；l1_k-1…lR_k-1、m1_k-1…mR_k-1、n1_k-1…nR_k-1分别为在第k-1个历元，第1到第R颗卫星的X、Y、Z三个方向的方向余弦；map1_k…mapR_k表示在第k个历元，第1到第R颗卫星的对流层投影函数；

其中，状态参数向量形式为：

$x_k=\left[\begin{array}{cccccccccccc}x& y& z& \stackrel{\·}{x}& \stackrel{\·}{y}& \stackrel{\·}{z}& \mathrm{dt}& d\stackrel{\·}{t}& \mathrm{ztd}& \mathrm{amb}1& ......& \mathrm{ambR}\end{array}\right]$

其中，xyz为位置参数，

为引入的速度状态参数，dt为接收机钟差参数，

为引入的接收机的钟漂状态参数，ztd为天顶对流层延迟，amb1…ambR分别为第1到第R颗卫星的模糊度参数；

所述在kalman滤波的动态模型里引入速度状态参数和钟漂状态参数，实现方式如下，kalman滤波的动态方程为：

x_k＝Φ_k，k-1x_k-1+w_k-1

其中，x_k、x_k-1分别为第k和第k-1个历元的状态参数，是未知数向量；Φ_k，k-1为第k-1到第k个历元的状态转移矩阵，w_k-1为过程噪声；引入了历元间差分之后的待求解的状态参数与上述观测模型一致：

$x_k=\left[\begin{array}{cccccccccccc}x& y& z& \stackrel{\·}{x}& \stackrel{\·}{y}& \stackrel{\·}{z}& \mathrm{dt}& d\stackrel{\·}{t}& \mathrm{ztd}& \mathrm{amb}1& ......& \mathrm{ambR}\end{array}\right]$ 相应的状态转移矩阵为：

${\mathrm{\Φ}}_{k,k-1}=\left[\begin{array}{cccc}{I}_c& {\mathrm{\Δt}\·I}_c& 0& 0\\ 0& I_c& 0& 0\\ 0& 0& {\mathrm{\Φ}}_{\mathrm{Ck},k-1}& 0\\ 0& 0& 0& I_a\end{array}\right]$

其中，I_c为3×3单位阵，I_a为R×R单位阵；

Φ_Ck，k-1为钟差的状态转移矩阵，形式为：

Δt是两个相邻历元的时间间隔；对应的过程噪声矩阵为：

$Q_{k,k-1}=E\[{\mathrm{ww}}^T\]=\left[\begin{array}{cccc}\left(S_P\frac{{\mathrm{\Δt}}^3}{3}\right)\·I_p& \left(S_P\frac{{\mathrm{\Δt}}^2}{2}\right)\·I_p& 0& 0\\ \left(S_P\frac{{\mathrm{\Δt}}^2}{2}\right)\·I_p& \left(S_P\mathrm{\Δt}\right)\·I_p& 0& 0\\ 0& 0& Q_c& 0\\ 0& 0& 0& \left(S_a\mathrm{\Δt}\right)I_a\end{array}\right]$

其中，

表示噪声的方差矩阵，S_P是位置参数的噪声因子，I_p是3×3单位阵，S_a是模糊度的噪声因子，I_a为R×R单位阵，

$Q_c=\left[\begin{array}{cc}{S}_{\mathrm{dt}}\mathrm{\Δt}+S_{d\stackrel{\·}{t}}\frac{{\mathrm{\Δt}}^3}{3}& S_{d\stackrel{\·}{t}}\frac{{\mathrm{\Δt}}^{\text{2}}}{2}\\ S_{d\stackrel{\·}{t}}\frac{\mathrm{\Δ}{t}^2}{2}& S_{d\stackrel{\·}{t}}\mathrm{\Δt}\end{array}\right]$

其中，S_dt是接收机钟漂的噪声因子。

Images