CN103344971A - 一种适用于gnss实时数据处理的优化方法 - Google Patents

Abstract

一种适用于GNSS实时数据处理的优化方法，采用实时递推算法，首先利用双频宽巷MW组合及电离层变化速度TECR组合，对周跳进行实时探测并修复，然后利用实时数据解算过程中产生的残差信息判断周跳修复的准确性。采用本发明所提出的GNSS实时数据优化方法，不仅能很好的发现并修复较小周跳，对提高实时解算精度也有很大帮助。

Description

一种适用于GNSS实时数据处理的优化方法

技术领域

本发明涉及一种适用于GNSS实时数据处理的优化方法，属于卫星导航应用技术领域。

背景技术

在GNSS定轨、定位与导航应用中，快速动态、高精度和稳定可靠一直以来是GNSS技术研究的主要内容，如何提高其定位精度是我们关注的重点之一。其中，载波相位观测值观测精度比伪距观测高几个量级，用于大地测量型接收机的载波相位观测精度可达到1-2mm，绝对定位可以达到厘米级，相对定位甚至可以达到毫米级的精度。因此，毫米级的精确定位与科学研究主要是利用载波相位观测资料。

在接收机进行连续的载波相位测量过程中，若由于某种原因而导致整周计数发生错误，就会使载波相位观测值较正确值出现一个整数周的跳跃，但不足一周的部分仍然正常，该跳跃被称为周跳。周跳的发生打断了载波相位跟踪测量的连续性，正确的周跳探测是正确确定整周模糊度参数的关键。GNSS实时数据质量控制是精密数据处理的重要内容，其主要工作为发现并及时处理数据中的周跳或粗差。

目前，国内外研究人员已经对周跳探测进行了深入的研究，提出了许多周跳探测的方法，如高次差法、多项式拟合法、小波分析法等，但大部分都是基于事后处理。另一主要的周跳处理方法为1990年Blewitt提出的TurboEdit方法，该方法主要利用宽巷组合及无几何伪距组合，但该方法对于小周跳不敏感，且不能实时修复周跳。

发明内容

本发明的技术解决问题是：克服现有技术的不足，提供了一种适用于GNSS实时数据处理的优化方法，以解决在精密数据处理过程中收敛缓慢，提升数据解算精度。

本发明的技术解决方案是：

一种适用于GNSS实时数据处理的优化方法，包括以下步骤：

（1）将观测卫星的GNSS双频观测数据组成MW宽巷组合及TECR组合，形成探测窗口；

（2）如果|N_wl(k-1)-N_wl(k)|≤4δ_N并且

则k历元载波不存在周跳，进入步骤（4），否则认为有周跳发生，进入步骤（3）；

ΔN₁(k)-ΔN₂(k)=N_wl(k-1)-N_wl(k)，在k历元L1和L2载波相位观测值Φ₁(k)和Φ₂(k)发生周跳，分别记为ΔN₁(k)和ΔN₂(k)，δ_N为MW宽巷组合的方差，δ_TECR为TECR组合的方差；

为TECR(k)的估值，且

$\mathrm{TE}\tilde{C}R\left(k\right)=\mathrm{TECR}(k-1)+{\mathrm{TECR}}^{\′}(k-1)\·\mathrm{\Δt},$ ${\mathrm{TECR}}^{\′}(k-1)=\frac{\mathrm{TECR}(k-1)-\mathrm{TECR}(k-2)}{\mathrm{\Δt}};$

（3）计算周跳，对载波数据进行修复，并利用

和

更新探测窗口；

（4）重复步骤（1）～（3），直到k历元所有观测卫星执行完毕；

（5）利用实时数据解算过程中产生的残差信息对观测值进行判断，若对于第j+1=k个历元而言,δ_j+1>c.δ_j或|e_j+1|>T成立，则认为该历元平差计算中存在有问题的观测值数据，之后进入步骤（6），否则进入步骤（7）；c和T均为预设常数；e_j+1为第_j+1历元的残差；

j和m均为历元数且j>m；

（6）放弃当前观测卫星的数据，返回步骤（5）重新解算生成残差并对剩余观测卫星的残差进行判断，直至满足步骤（5）中的条件并且进入步骤（7）；

（7）本历元结束，执行下一历元。

所述将观测卫星的GNSS双频码伪距观测数据组成MW宽巷组合及TECR组合具体为：

MW宽巷组合整周模糊度

其中，Φ₁和Φ₂分别为载波L1和L2的载波相位观测值；P₁和P₂为码相位观测值；f₁为L1载波频率，f₂为L2载波频率，λ_wl为MW宽巷组合波长；

电离层总电子含量变化率其中，Δt为相邻历元时间间隔， $\mathrm{TEC}\left(k\right)=\frac{{f_1}^2\left\{\right[{\mathrm{\λ}}_1{\mathrm{\Φ}}_1\left(k\right)-{\mathrm{\λ}}_2{\mathrm{\Φ}}_2\left(k\right)]-[{\mathrm{\λ}}_1{N}_1-{\mathrm{\λ}}_2{N}_2\left]\right\}}{40.3\×{10}^{16}\·(\mathrm{\γ}-1)},$ $\mathrm{\γ}={f_1}^2/{f_2}^2,$ λ₁和λ₂依次为L1和L2的波长，N₁和N₂依次为L1和L2的整周模糊度。

所述计算周跳，对载波数据进行修复具体为：

通过求解方程组

$\left\{\begin{array}{c}{\mathrm{\ΔN}}_1\left(k\right)-{\mathrm{\ΔN}}_2\left(k\right)=N_{\mathrm{wl}}(k-1)-N_{\mathrm{wl}}\left(k\right)\\ {\mathrm{\λ}}_1{\mathrm{\ΔN}}_1\left(k\right)-{\mathrm{\λ}}_2{\mathrm{\ΔN}}_2\left(k\right)=\frac{40.3\×{10}^{16}(\mathrm{\γ}-1)\mathrm{\Δt}\·\mathrm{TECR}\left(k\right)}{{f_1}^2}-{\mathrm{\λ}}_1[{\mathrm{\Φ}}_1\left(k\right)-{\mathrm{\Φ}}_1(k-1)]+{\mathrm{\λ}}_2[{\mathrm{\Φ}}_2\left(k\right)-{\mathrm{\Φ}}_2(k-1)]\end{array}\right.$

计算周跳ΔN₁(k)和ΔN₂(k)，通过公式 $\left\{\begin{array}{c}\widehat{{\mathrm{\Φ}}_1}\left(k\right)={\mathrm{\Φ}}_1\left(k\right)+\mathrm{\Δ}{N}_1\left(k\right)\\ \widehat{{\mathrm{\Φ}}_2}\left(k\right)={\mathrm{\Φ}}_2\left(k\right)+\mathrm{\Δ}{N}_2\left(k\right)\end{array}\right.$ 对载波数据进行修复，

和

为修复周跳后的载波相位观测值。

本发明与现有技术相比的有益效果是：

（1）在传统GNSS数据预处理基础上，重新挖掘MW宽巷组合特性，并引入TECR组合，实现了GNSS载波周跳的实时探测与修复。实验结果证明，该过程实时修复精度较高，可准确探测并修复1周小周跳，且周跳修复结果较稳定，完全适用于高精度GNSS精密数据处理中。

（2）为提高实时周跳修复的准确定，提升数据解算精度，该优化方法中引入了数据解算过程中生成的残差信息，利用残差信息对小周跳或粗差敏感特性，设计了检验条件，提高了GNSS数据解算的准确性。

（3）依据如上多种线性组合及残差信息，设计了适用于GNSS实时数据处理的优化流程。

附图说明

图1为本发明流程图。

具体实施方式

下面结合附图对本发明的具体实施方式进行进一步的详细描述。

本发明提供的适合于GNSS实时数据处理的优化方法，首先利用观测卫星双频观测值形成Melbourne-Wübbena（MW）宽巷组合及电离层总电子含量变化率（TECR）组合，组成探测窗口，然后利用该窗口对周跳进行实时探测并对其进行修复，最后利用解算过程中产生的残差信息进行判断，提升解算精度。其流程见图1，具体包括以下几个步骤：

（1）将观测卫星的GNSS双频观测数据组成MW宽巷组合及TECR组合，形成探测窗口；

假设前k-1历元没有周跳，在k历元载波L1和L2载波相位观测值Φ₁(k)和Φ₂(k)发生周跳，分别记为ΔN₁(k)和ΔN₂(k)，则在k历元正确的载波相位观测值应表示为：

$\widehat{{\mathrm{\Φ}}_1}\left(k\right)={\mathrm{\Φ}}_1\left(k\right)+{\mathrm{\ΔN}}_1\left(k\right)---\left(1\right)$

$\widehat{{\mathrm{\Φ}}_2}\left(k\right)={\mathrm{\Φ}}_2\left(k\right)+{\mathrm{\ΔN}}_2\left(k\right)---\left(2\right)$

其中，

和

为修复周跳后的载波相位观测值。

MW宽巷组合：

$N_{\mathrm{wl}}=N_1-N_2={\mathrm{\Φ}}_1-{\mathrm{\Φ}}_2-\frac{f_1\·P_1+f_2\·P_2}{{\mathrm{\λ}}_{\mathrm{wl}}(f_1+f_2)}---\left(3\right)$

其中λ_wl、N_wl分别称为宽巷波长及宽巷整周模糊度。由式（1）-（3）可得：

ΔN₁(k)-ΔN₂(k)=N_wl(k-1)-N_wl(k)   （4）

TECR组合：

根据电离层总电子含量（TEC）定义及载波相位观测值的原理性观测方程可推导出TEC的数学公式，如下：

$\mathrm{TEC}\left(k\right)=\frac{{f_1}^2\left\{\right[{\mathrm{\λ}}_1{\mathrm{\Φ}}_1\left(k\right)-{\mathrm{\λ}}_2{\mathrm{\Φ}}_2\left(k\right)]-[{\mathrm{\λ}}_1{N}_1-{\mathrm{\λ}}_2{N}_2\left]\right\}}{40.3\×{10}^{16}\·(\mathrm{\γ}-1)}---\left(5\right)$

由式（6）可推算出相邻历元的电离层总电子含量变化率(TECR)：

$\mathrm{TECR}\left(k\right)=\frac{\mathrm{TEC}\left(k\right)-\mathrm{TEC}(k-1)}{\mathrm{\Δt}}---\left(6\right)$

其中Δt为相邻历元时间间隔。整理式（5）和（6），则k历元的周跳可估计为：

${\mathrm{\λ}}_1{\mathrm{\ΔN}}_1\left(k\right)-{\mathrm{\λ}}_2{\mathrm{\ΔN}}_2\left(k\right)=\frac{40.3\×{10}^{16}(\mathrm{\γ}-1)\mathrm{\Δt}\·\mathrm{TECR}\left(k\right)}{{f_1}^2}----\left(7\right)$

${\mathrm{\λ}}_1[{\mathrm{\Φ}}_1\left(k\right)-{\mathrm{\Φ}}_1(k-1)]+{\mathrm{\λ}}_2[{\mathrm{\Φ}}_2\left(k\right)-{\mathrm{\Φ}}_2(k-1)]$

由于要在k历元探测与修复周跳，且不能确定k历元是否存在周跳，因此不能直接使用式（6）估计TECR(k)，因此采用递推方法。例如在k历元，利用k-1历元和k-2历元的观测值来获得TECR(k)。具体步骤为：首先，在k-1历元计算该历元TECR的变化率TECR′(k-1)，如式（8）；然后计算TECR的估值

如式（9）。为提高精度，可对其进行拟合，滤去噪声影响。

${\mathrm{TECR}}^{\′}(k-1)=\frac{\mathrm{TECR}(k-1)-\mathrm{TECR}(k-2)}{\mathrm{\Δt}}---\left(8\right)$

$\mathrm{TE}\tilde{C}R\left(k\right)=\mathrm{TECR}(k-1)+{\mathrm{TECR}}^{\′}(k-1)\·\mathrm{\Δt}---\left(9\right)$

（2）如果|N_wl(k-1)-N_wl(k)|≤4δ_N并且则k历元载波不存在周跳，进入步骤（4），否则认为有周跳发生，进入步骤（3）；

条件一：利用MW宽巷组合探测周跳。利用式（10）可递推出k历元宽巷模糊度的均值E(N_k)和方差

如果通过式（4）计算得到的|N_wl(k-1)-N_wl(k)|≤4δ_N，则k历元不存在周跳。否则，就认为有周跳发生。

$\left\{\begin{array}{c}E\left(N_k\right)=E\left(N_{k-1}\right)+\{N_k-E\left(N_{k-1}\right)\}/k\\ E\left(N_k^2\right)=E\left(N_{k-1}^2\right)+\{N_k^2-E\left(N_{k-1}^2\right)\}/k\\ {\mathrm{\δ}}_{N\left(k\right)}^2=E\left(N_k^2\right)-{\left[E\left(N_k\right)\right]}^2\end{array}---\left(10\right)\right.$

条件二：利用TECR组合探测周跳。先于k历元的TECR数据的均值和方差可以通过与式（10）相似的方法递推出来。认为TECR在较短的时间内是一个常量，把由式（6）计算出来的k历元处的TECR与式（9）估计出的k历元处的TECR进行比较。如果则k历元不存在周跳。否则，就认为有周跳发生并用式（7）来确定[λ₁ΔN₁(k)-λ₂ΔN₂(k)]的大小。

（3）计算周跳，对载波数据进行修复，并利用和

更新探测窗口；

只要满足上述条件中的一个，就认为有周跳发生。联立方程（4）和（7）

$\left\{\begin{array}{c}{\mathrm{\ΔN}}_1\left(k\right)-{\mathrm{\ΔN}}_2\left(k\right)=N_{\mathrm{wl}}(k-1)-N_{\mathrm{wl}}\left(k\right)\\ {\mathrm{\λ}}_1{\mathrm{\ΔN}}_1\left(k\right)-{\mathrm{\λ}}_2{\mathrm{\ΔN}}_2\left(k\right)=\frac{40.3\×{10}^{16}(\mathrm{\γ}-1)\mathrm{\Δt}\·\mathrm{TECR}\left(k\right)}{{f_1}^2}-{\mathrm{\λ}}_1[{\mathrm{\Φ}}_1\left(k\right)-{\mathrm{\Φ}}_1(k-1)]+{\mathrm{\λ}}_2[{\mathrm{\Φ}}_2\left(k\right)-{\mathrm{\Φ}}_2(k-1)]\end{array}\right.$

四舍五入则可解算出ΔN₁(k)、ΔN₂(k)。在方程联立解算过程中，顾及到宽巷组合受噪声影响较大，首先对方程（4）进行取整，滤去噪声对前后历元的共性影响，然后与方程（7）联立取整。实际数据验证该策略获得了较高的解算精度。然后通过公式 $\left\{\begin{array}{c}\widehat{{\mathrm{\Φ}}_1}\left(k\right)={\mathrm{\Φ}}_1\left(k\right)+\mathrm{\Δ}{N}_1\left(k\right)\\ \widehat{{\mathrm{\Φ}}_2}\left(k\right)={\mathrm{\Φ}}_2\left(k\right)+\mathrm{\Δ}{N}_2\left(k\right)\end{array}\right.$ 对载波数据进行修复更新，并利用

和

及双频码观测数据重新形成MW宽巷组合及TECR组合，更新探测窗口。

（4）重复步骤（1）～（3），直到k历元所有观测卫星执行完毕；

（5）利用实时数据解算过程中产生的残差信息对观测值进行判断，若对于第j+1=k个历元而言，δ_j+1>c.δ_j或|e_j+1|>T成立，则认为该历元平差计算中存在有问题的观测值数据，之后进入步骤（6），否则进入步骤（7）；

为提高数据解算精度，引入实时数据解算过程中产生的残差信息。观测值验后残差随着时间不断发生变化，其变化是有规律且平滑的，特别是相邻历元间的变化更加平稳，主要出现偶然误差的特性，变化极小，然而周跳或粗差的出现将破坏其平稳的规律性。据此，可以设计根据残差信息进行数据质量控制的探测条件。

假设从第m历元开始某颗卫星开始参与平差计算，在j（j>m）历元，可通过解算，获得该卫星从m到j历元的解算残差信息矩阵

由下公式可获得第j历元的残差中误差：

${\mathrm{\δ}}_j=\sqrt{\frac{\underset{i=m}{\overset{j}{\mathrm{\Σ}}}{e}_i\·e_i}{j-m+1}}---\left(11\right)$

那么对于第j+1=k个历元而言，若满足以下条件：

δ_j+1>c.δ_j或|e_j+1|>T   （12）

即，该卫星第j+1历元残差的标准差达到第j历元标准差的c倍，或者第j+1历元残差的绝对值大于某一值T（其中c及T为经验预设常数，例如取c=1.4，T=6δ_j），则认为该历元平差计算中存在不“干净”的观测值数据。

（6）放弃当前观测卫星的数据，返回步骤（5）重新解算生成残差并对剩余观测卫星的残差进行判断，直至满足步骤（5）中的条件并且进入步骤（7）；

（7）本历元结束，执行下一历元。

为了测试上述流程的特性及对周跳的敏感性，选取IGS ALGO站2010年4月2日采集的采样率为1s的GPS数据进行测试，选取2号卫星作为测试卫星（截至高度角为10°）。

表1给出了模拟一系列的周跳对MW宽巷组合、TECR组合及SRIF模块验后残差的变化，从表中可以看出对于高频数据而言（如1Hz），当发生小周跳时，MW组合观测值变化很小，然而TECR值变化很剧烈，因此TECR组合可用于探测高频小周跳。当GPS发生f₁/f₂=77/60比例的周跳时，TECR组合探测失效，此时可以使用MW进行周跳探测。因此综合使用MW组合观测值和TECR，可有效地进行周跳探测与修复。

表1：GPS周跳产生对MW款项组合、TECR组合及验后残差变化

表中： $C={f_1}^2/[40.3\×{10}^{16}(\mathrm{\γ}-1)].$

表中模拟了小周跳对验后残差的影响，目的是检测残差信息对小周跳的敏感性（其中将观测数据每500个历元分为一观测弧段，求解其δ_j）。在上文中分析到，在周跳解算过程中可能由于其它一系列误差的影响，导致根据修复的周跳值可能与正确值相差1～2周。从表中可以看出，验后残差信息可有效地反映周跳修复错误历元残差的变化，由此可准确的确定周跳修复错误的卫星及历元，提高实时处理的准确性。

Claims (3)

1.一种适用于GNSS实时数据处理的优化方法，其特征在于包括以下步骤：

（1）将观测卫星的GNSS双频观测数据组成MW宽巷组合及TECR组合，形成探测窗口；

（2）如果|N_wl(k-1)-N_wl(k)|≤4δ_N并且

则k历元载波不存在周跳，进入步骤（4），否则认为有周跳发生，进入步骤（3）；

ΔN₁(k)-ΔN₂(k)=N_wl(k-1)-N_wl(k)，在k历元L1和L2载波相位观测值Φ₁(k)和Φ₂(k)发生周跳，分别记为ΔN₁(k)和ΔN₂(k)，δ_N为MW宽巷组合的方差，δ_TECR为TECR组合的方差；

为TECR(k)的估值，且

$\mathrm{TE}\tilde{C}R\left(k\right)=\mathrm{TECR}(k-1)+{\mathrm{TECR}}^{\′}(k-1)\·\mathrm{\Δt},$ ${\mathrm{TECR}}^{\′}(k-1)=\frac{\mathrm{TECR}(k-1)-\mathrm{TECR}(k-2)}{\mathrm{\Δt}};$

（3）计算周跳，对载波数据进行修复，并利用

和

更新探测窗口；

（4）重复步骤（1）～（3），直到k历元所有观测卫星执行完毕；

（5）利用实时数据解算过程中产生的残差信息对观测值进行判断，若对于第j+1=k个历元而言,δ_j+1>c.δ_j或|e_j+1|>T成立，则认为该历元平差计算中存在有问题的观测值数据，之后进入步骤（6），否则进入步骤（7）；c和T均为预设常数；e_j+1为第_j+1历元的残差；

，j和m均为历元数且j>m；

（6）放弃当前观测卫星的数据，返回步骤（5）重新解算生成残差并对剩余观测卫星的残差进行判断，直至满足步骤（5）中的条件并且进入步骤（7）；

（7）本历元结束，执行下一历元。

2.根据权利要求1所述的一种适用于GNSS实时数据处理的优化方法，其特征在于：所述将观测卫星的GNSS双频码伪距观测数据组成MW宽巷组合及TECR组合具体为：

MW宽巷组合整周模糊度

Φ₁和Φ₂分别为载波L1和L2的载波相位观测值；P₁和P₂为码相位观测值；f₁为L1载波频率，f₂为L2载波频率，λ_wl为MW宽巷组合波长；

电离层总电子含量变化率

其中，Δt为相邻历元时间间隔， $\mathrm{TEC}\left(k\right)=\frac{{f_1}^2\left\{\right[{\mathrm{\λ}}_1{\mathrm{\Φ}}_1\left(k\right)-{\mathrm{\λ}}_2{\mathrm{\Φ}}_2\left(k\right)]-[{\mathrm{\λ}}_1{N}_1-{\mathrm{\λ}}_2{N}_2\left]\right\}}{40.3\×{10}^{16}\·(\mathrm{\γ}-1)},$ $\mathrm{\γ}={f_1}^2/{f_2}^2,$ λ₁和λ₂依次为L1和L2的波长，N₁和N₂依次为L1和L2的整周模糊度。

3.根据权利要求1所述的一种适用于GNSS实时数据处理的优化方法，其特征在于：所述计算周跳，对载波数据进行修复具体为：

通过求解方程组

$\left\{\begin{array}{c}{\mathrm{\ΔN}}_1\left(k\right)-{\mathrm{\ΔN}}_2\left(k\right)=N_{\mathrm{wl}}(k-1)-N_{\mathrm{wl}}\left(k\right)\\ {\mathrm{\λ}}_1{\mathrm{\ΔN}}_1\left(k\right)-{\mathrm{\λ}}_2{\mathrm{\ΔN}}_2\left(k\right)=\frac{40.3\×{10}^{16}(\mathrm{\γ}-1)\mathrm{\Δt}\·\mathrm{TECR}\left(k\right)}{{f_1}^2}-{\mathrm{\λ}}_1[{\mathrm{\Φ}}_1\left(k\right)-{\mathrm{\Φ}}_1(k-1)]+{\mathrm{\λ}}_2[{\mathrm{\Φ}}_2\left(k\right)-{\mathrm{\Φ}}_2(k-1)]\end{array}\right.$

计算周跳ΔN₁(k)和ΔN₂(k)，通过公式 $\left\{\begin{array}{c}\widehat{{\mathrm{\Φ}}_1}\left(k\right)={\mathrm{\Φ}}_1\left(k\right)+\mathrm{\Δ}{N}_1\left(k\right)\\ \widehat{{\mathrm{\Φ}}_2}\left(k\right)={\mathrm{\Φ}}_2\left(k\right)+\mathrm{\Δ}{N}_2\left(k\right)\end{array}\right.$ 对载波数据进行修复，

和

为修复周跳后的载波相位观测值。

Images