利用聚类法实现PPP 固定解的快速静态定位

罗玉山，张小红

（武汉大学 测绘学院，武汉 430079）

0 引言

文献[1]于 1997 年提出了精密单点定位（precise point positioning,PPP）的概念。随着全球卫星导航系统（global navigation satellite system,GNSS）的逐步建成，可见卫星数量增多，PPP 相关技术发展迅速。PPP 是指利用国际GNSS 服务组织（International GNSS Service,IGS）或者其他机构提供的高精度卫星轨道和钟差产品，通过模型改正、参数估计、随机模型控制等方法精确改正各项误差，并利用单台接收机的伪距和相位观测值进行定位解算，实现高精度定位的方法[2]。

随着PPP 技术的发展，浮点解技术逐步完善，国内外的研究热点也逐步转向固定解技术。固定解技术的核心难题是恢复载波相位测量中的模糊度的整数特性，从而达到更高的定位精度。实现PPP 固定解的方法有相位小数偏差（fractional cycle bias,FCB）方法、整数钟法、钟差去耦法等，都成功地提高了PPP 浮点解的定位精度[3]。针对不同的应用场景，PPP 可分为动态解和静态解。其中，静态PPP 固定解以其毫米级的定位精度被广泛地用于大地测量、工程测量等领域。

PPP 技术的理论、模型、算法已基本成熟，在定位模型研究、数据预处理、参数估计方法等方面均取得了一系列成果，但是在定位解算结果的处理和优化上的相关研究仍然较少[3]。高精度PPP 固定解的出现显著改善了PPP的定位精度；但是，由于受到各误差影响，固定解解算通常需要较长的初始化时间，且在单系统观测条件下，通常需要30 min 甚至更长的初始化时间，极大地延长了外出测量作业的时间，影响了作业效率。尤其是在铁路、公路、输电线路、管线等施工测量中，数据采集时间往往比较短，因此如何利用有限的短时间数据解算出最优的定位结果十分重要。本文提出一种利用位置聚类法对短时间的静态PPP 固定解结果进行优化的方法，并利用全球IGS 跟踪站数据验证该方法的精度和可靠性。

1 理论方法

1.1 静态PPP 原理

本文静态PPP 解算均采用GPS 单系统、双频无电离层组合函数模型。静态观测方程为：

式中：下标1、2、IF 分别代表频率1、频率2 和无电离层组合；PIF、LIF分别为伪距和相位观测值；ρ为卫地距；c为光速；dtr,IF、为接收机钟差和卫星钟差；TIF为对流层延迟；NIF为模糊度；λIF为无电离层组合波长；εP、εL为伪距观测噪声和相位观测噪声；X代表伪距观测值或者相位观测值；f1、f2为卫星载波频率；α、β为频率组合比例系数。

精密单点定位中的主要误差处理方法如表1所示。

表1 误差源及其处理方法

1.2 整数模糊度固定方法

本文采用FCB 方法实现了静态PPP 固定解的解算。无电离层组合模糊度NIF不具有整数特性，无法直接固定，可以将其分解为宽巷（wide lane,WL）与窄巷（narrow lane,NL）模糊度分别进行固定，即可实现固定解，可表示为

式中：NWL为宽巷模糊度；NNL为窄巷模糊度。

FCB 方法从改正或消除卫星和接收机端的小数周偏差的角度出发，首先对宽巷墨尔本-维贝纳（Melbourne-Wübbena,MW）组合平滑得到的宽巷模糊度作星间单差，再利用宽巷FCB 产品进行改正以恢复宽巷整周模糊度；对于相关性较强的窄巷模糊度，则利用窄巷FCB 产品进行改正并通过最小二乘模糊度降相关平差（least-square ambiguity decorrelation adjustment,LAMBDA）方法搜索得到窄巷整周模糊度。当宽巷与窄巷模糊度都成功固定时，即可获得无电离层组合的PPP 整周模糊度[4]。

因此，本文使用了武汉大学测绘学院公开发布的FCB 产品，即相位小数周偏差改正数后，PPP 解算的模糊度参数具有整数特性，可实现固定解。

1.3 模糊度检核指标

PPP 固定解精度提高的前提是模糊度参数被正确固定。得到模糊度固定解后，须对其质量进行检核和验证；因为固定错误得到的模糊度参数会导致分米级甚至米级的误差，严重影响定位精度[5]。

1）模糊度精度因子（ambiguity dilution of precision,ADOP）是一种有效判断整周模糊度是否固定正确的模糊度检核指标。文献[6]于1997 年引入了ADOP 这一指标，用来描述模糊度参数的精度特性，可表示为

式中：det(QNˆ)表示模糊度参数协方差阵的行列式的值；n为模糊度估值的数量。

ADOP 值极高程度地反映了模糊度的平均精度，因为该值考虑了模糊度对应方差和协方差的所有信息[7]。文献[8]指出，当其小于0.12 周时，模糊度的固定成功率接近100%。

2）模糊度固定成功率（bootstrapping）给出了模糊度正确固定概率的量化信息，体现了数据处理模型的强度[9]，可表示为：

式中：Ni,(l)是Ni,(i−1,…,1)的缩写，表示第i个模糊度，是参考（i– 1）个整数模糊度得到的估值；σ为Ni,(l)的标准差。

成功率指标对模糊度的检核效果好，同时计算方便，可以反映模糊度整数解的全局质量[10]。

3）模糊度比率检验（ratio-test）指标体现了浮点解与最优整数解向量的接近程度，即次优整数解残差2 次型与最优整数解残差2 次型的比值[11]，可表示为

通常利用ratio-test 指标的经验阈值来判断模糊度是否固定正确，不足之处在于不能给出固定解全局质量[12]。

所以，合理地利用ADOP 值、bootstrapping 成功率、ratio-test 指标可以有效地反映出整周模糊度的固定准确度，且整周模糊度固定是否准确将直接影响到固定解的定位精度，因此可把这3 个指标作为PPP 固定解结果的检核指标。

1.4 利用聚类法实现快速静态定位

本文所提出的快速静态定位方法，是指针对短时间的观测数据进行解算和相关算法优化后得到最优定位结果的过程。该方法的主要思想是，通过正反向的卡尔曼（Kalman）滤波解算来增大短时间数据所解算的结果数量，利用误差分布特性去寻找最密集的坐标空间，通过相关的模糊度指标对固定解进行质量检核，最终得到快速静态定位的优化结果。算法流程如图1 所示。

图1 算法流程

算法的关键步骤如下：

1）正反向Kalman 滤波解算。对短时间数据进行事后的正向和反向Kalman 滤波解算，并采用FCB法固定模糊度参数。其中，正反向滤波解算可以增加解算结果的数量，解决短时间数据固定解结果不足的问题。解算得到的信息包括总历元数A、测站的3 维坐标以及每个历元的ADOP 值、bootstrapping成功率、ratio-test 指标、固定解状态等。

2）解算结果预处理。PPP 解算的前15 min 通常为初始化收敛时间，解算结果质量往往较差。为提高优化方法的结果质量，选择去除这段PPP 初始化时间内的解算结果。统计解算结果的固定解总数A1，如果A1≥0.5A，则去除未固定历元的解算结果，仅采用定位精度更高的固定解结果进行结果优化。

3）确定位置聚类空间。统计处理后剩余历元的3 维坐标结果，并找到其最大值和最小值。根据这3 组阈值便可确定位置聚类的3 维空间，且所有历元的3 维坐标(xi,y i,zi)均满足

4）寻找位置聚类空间的最密集区域。对于静态定位而言，测站的3 维坐标是一个确定的坐标，较优的静态解算结果总是较紧密地分布在这个确定值的周围，而较劣的结果则较分散地分布在远离这个确定值的区域。所以为了实现PPP 固定解的快速静态定位，需要找到较优解的空间聚集区域，即位置聚类空间的最密集区域。

首先计算

图2 分割平面示意图

从图可以看出，通过循环迭代地分割坐标点数量最多的区域，可无限逼近最密集区域的中心。同理，将该方法推广到3 维空间，可通过有限次的分割，找到坐标空间的最密集区域。统计各个空间的坐标点数量ni（i=1,2,…,8）以及坐标点数量最大值nmax，选择保留所有ni≥ 0.5nmax（i=1,2,…,8）的坐标点空间。通过数据实验发现，坐标点数量最多的区域并不一定能得到最好的优化结果，所以将数量较多的几个空间都纳入参考，并通过精度指标对参考空间进行评判。

5）寻找最优的坐标空间。计算保留下来的坐标空间的3 维坐标标准差σx、σy、σz，标准差可体现一个空间里的坐标点的离散程度，标准差越小，坐标点分布越紧密，计算公式为

式中：k为坐标空间的坐标点总数量；变量r表示3 维坐标x、y、z；下标i表示坐标空间的第i个点。

式中mADOP,i为第i个坐标点对应的ADOP 值。

同时，计算bootstrapping 成功率、ratio-test 指标的平均值。通过这4 个指标，从坐标点的空间聚集程度、固定解模糊度的固定准确度2 个方面对坐标空间进行评判，即可从多个参考空间中筛选出1 个最优的坐标空间。重复步骤2）、步骤4）、步骤5），不断地分割位置聚类空间，缩小坐标空间范围，逐步逼近位置聚类空间的最密集区域并根据各项指标判定最优坐标空间，最后对最优坐标空间的坐标点（xi、yi、zi）求平均值，得到最终的PPP 固定解快速静态定位结果。

2 实验与结果分析

本文选取了IGS 数据中心的DARW、JFNG、KIRU 等6 个测站2018-01-01的数据进行解算实验，观测数据采样率为30 s，精密数据产品采用IGS 提供的精密星历和钟差产品，并采用武汉大学测绘学院公开发布的FCB 产品以实现PPP 固定解。为实现短时间数据的PPP 静态快速固定解，本文的观测数据均为1 h，即1 d的数据分为24 段，每段数据长度为1 h，并分别进行单系统的静态PPP固定解解算，其中单系统指美国全球定位系统（global positioning system,GPS）。

为体现优化方法对于短时间静态PPP 固定解精度的提升效果，本文从东（E）、北（N）、天顶（U）3 个方向进行对比，分析快速静态定位优化结果与浮点解、固定解的精度差异。分别计算每个测站各自的24 段数据，可以得到每1 h的浮点解、固定解均值，以及PPP 固定解快速静态定位结果（以下简称优化结果），其中：浮点解指1 h的浮点解收敛结果，即最后1 个历元的解；固定解均值指该段静态数据所解得所有收敛点后的固定解的平均值。为更加直观、便捷地进行误差的对比分析，每一个测站的浮点解、固定解均值和优化结果均取24 组对应结果的平均值。

图3～图5 分别为浮点解与优化结果基于真值的E、N、U 方向的偏差对比。从图中可以看出：在E 方向上，优化结果的偏差均优于浮点解，且偏差均小于3 cm，KIRU、CEDU、CHPG、AREG测站的偏差低于 1.5 cm，而浮点解偏差均值为3.5 cm，因此在E 方向上，优化结果的精度更优；在N 方向上，优化结果的偏差均优于浮点解，且优化结果的偏差均值为1.02 cm，而浮点解的偏差均值为1.38 cm，有毫米级的提升，所以在N 方向上，优化结果的精度略优；在U 方向上，优化结果的偏差均小于浮点解，且均值为2.64 cm，而浮点解偏差均值为4.39 cm，大多数测站有厘米级的精度提升，所以优化结果在U 方向上的精度也更优。

图3 浮点解与优化结果东（E）方向偏差对比

图4 浮点解与优化结果北（N）方向偏差对比

图5 浮点解与优化结果天顶（U）方向偏差对比

图6～图8 分别为固定解均值与优化结果在E、N、U 方向的偏差对比。从图中可以看出：在E 方向上，大多数测站的固定解均值和优化结果的偏差相差为1 cm 以上，优化结果在E 方向上的精度更高；在N 方向上，固定解均值偏差均值为1.37 cm，优化结果偏差均值为1.02 cm，优化结果的精度略优；在U 方向上，优化结果偏差均小于固定解均值偏差，且固定解均值偏差的平均值为3.97 cm，优化结果偏差均值为2.64 cm，所以优化结果在U方向上的精度也更优。

图6 固定解均值与优化结果东（E）方向偏差对比

图7 固定解均值与优化结果北（N）方向偏差对比

图8 固定解均值与优化结果天顶（U）方向偏差对比

综上所述，相对于浮点解以及固定解，优化结果的精度在E、N、U 方向上都有所提升。为进一步证明优化结果定位精度的稳定性，本文分别计算了这3 种解算结果在E、N、U 方向上的均方根误差（root mean square error，RMSE），其结果如表2 所示。

表2 浮点解、固定解均值及优化结果RMSE 统计 m

从表中可以看出：相对于浮点解，优化结果的精度在E方向提升了3.58 cm，N方向提升了0.7 cm，U 方向提升了2.24 cm；相对于固定解均值，优化结果的精度在E 方向提升了2.77 cm，N 方向提升了0.56 cm，U 方向提升了1.18 cm。这些数据表明，利用位置聚类方法优化的静态固定解结果的精度是最优的。

为了解该算法的计算效率，本文统计了DARW 测站24 组数据的内存消耗和运行时间的平均值，最终结果表明程序内存消耗低，运行时间短，算法效率高。

3 结束语

本文基于静态 PPP 短时作业时精度不高的应用现状，提出了利用位置聚类来优化静态PPP固定解的方法。通过模糊度精度因子、坐标标准差等指标寻找位置域最密集空间，得到优化结果，并通过数据解算分析，验证了该方法所得优化结果的精度相较于浮点解以及固定解均有提升，有较高的精度和稳定性，且算法的计算效率较高，可为短时的快速静态定位作业的实际应用提供参考。