目前实时动态PPP的解算策略主要采用序贯最小二乘与Kalman滤波两种算法。李冲等^[1]研究了Kalman滤波模型在精密单点定位中的应用，并针对其在数据异常时易发散的特点，将抗差自适应滤波模型运用到精密单点定位中，得到了较好的结果。许长辉等^[2]研究了基于Kalman滤波残差的抗差方法，并与基于新息的抗差方法进行对比，得到Kalman滤波新息作为抗差等处理的变量不可靠，新息表示的残差作为变量更有效可靠的结论。丁文武等^[3]提出一种充分顾及电离层预报精度的实时快速重新初始化算法。陈泉余等^[4]针对PPP中各类参数和观测值的特性以及观测星座的空间几何结构，提出一种适用于PPP的抗差自适应滤波算法。上述研究主要集中于Kalman抗差自适应算法以及对残差的研究，对于历元间的约束信息涉及较少。

序贯最小二乘法因其算法简单、无需先验信息、不易发散等优点，被广泛应用于实时精密单点定位解算中。但常规序贯最小二乘实时PPP算法将坐标参数视为时变参数，历元间不继承，损失了历元间先验信息的约束。Kalman滤波中常用的CV(常速度)模型、CA(常加速度)模型都适用于运动状态已知且稳定的载体，但运动载体在实际运动状态中很难保证状态稳定，若状态发生改变，先验信息存在偏差都会导致状态估计的发散。基于此，本文探讨一种顾及实时位置预测信息的动态精密单点定位算法，并给出具体算法实现流程。通过实验数据进行验证分析，结果表明，该方法能有效改善高频动态PPP的精度，对无基准站支持下的缓变型形变监测具有一定的应用价值。

1 数学模型

PPP数据处理过程中，通常采用观测值组合、模型改正及参数估计对传播误差进行消除或减弱。常用的PPP模型有无电离层(ionosphere-free, IF)组合模型^[5]、UofC模型、基于原始观测值的非组合模型。为了消除电离层造成的延迟误差，本文构造IF组合观测值，对于和卫星有关的轨道和钟差通过IGS提供的精密星历和钟差进行改正，其他误差通过常规模型改正，简化公式为：

$ \begin{array}{l} {P_{{\rm{IF}}}} = \frac{{f_1^2 \cdot {P_1}-f_2^2 \cdot {P_2}}}{{f_1^2-f_2^2}}\; = \\ \;\;\;\;\;\;\rho + {\rm{cdt}} + {d_{{\rm{trop}}}} + {\varepsilon _{{P_{{\rm{IF}}}}}} \end{array} $

(1)

$ \begin{array}{l} {L_{{\rm{IF}}}} = \frac{{f_1^2 \cdot {L_1}-f_2^2 \cdot {L_2}}}{{\;f_1^2-f_2^2}}\; = \\ \;\;\;\;\;\;\rho + {\rm{cdt}} + {d_{{\rm{trop}}}} + {N_1} + {\varepsilon _{{L_{{\rm{IF}}}}}}\; \end{array} $

(2)

式中，P_IF、L_IF分别表示伪距和载波LC组合观测值，f_i表示频率，ρ表示卫星到测站的距离，c表示光速，d_trop表示未被模型化的对流层湿延迟，N₁为组合整周模糊度，ε_PIF为伪距组合测量值的噪声以及未被模型化的误差改正，ε_LIF为载波组合测量值的噪声以及未被模型化的误差改正。

根据所估计参数的特性，利用序贯最小二乘进行参数解算，得到如下参数估计值：

$ {{\mathit{\boldsymbol{\hat x}}}_k} = \;{{\mathit{\boldsymbol{\hat x}}}_{k{\rm{-}}1}} + {\mathit{\boldsymbol{Q}}_{{{\mathit{\hat x}}_\mathit{k}}}}{\mathit{\boldsymbol{A}}^{\rm{T}}}\;{\mathit{\boldsymbol{P}}_k}(\mathit{\boldsymbol{l}}{\rm{-}}{\mathit{\boldsymbol{A}}_k}{{\mathit{\boldsymbol{\hat x}}}_{k{\rm{-}}1}}) $

(3)

$ {\mathit{\boldsymbol{Q}}_{{{\mathit{\hat x}}_\mathit{k}}}} = {(\mathit{\boldsymbol{Q}}_{{{\mathit{\hat x}}_{k{\rm{-}}1}}}^{-1} + \mathit{\boldsymbol{A}}_k^{\rm{T}}{\mathit{\boldsymbol{P}}_k}{\mathit{\boldsymbol{A}}_k})^{-1}} $

(4)

式中，${\mathit{\boldsymbol{\hat x}}_\mathit{k}}$为当前历元状态参数，k代表历元数，${\mathit{\boldsymbol{Q}}_{{{\mathit{\hat x}}_\mathit{k}}}}$为参数协因数阵，A为观测方程线性化后的系数矩阵，P_k为观测值的权矩阵，l为单位阵。

实时动态精密单点定位解算过程中，对于时变参数通常是将其进行约化，因此历元间的坐标关系无法继承，而时不变参数则利用序贯最小二乘算法递推求解。静态序贯最小二乘算法认为坐标为时不变参数，考虑了历元间的坐标约束，但只适用于静态观测数据，对于缓变型载体并不适用；动态序贯最小二乘算法认为坐标为时变参数，没有考虑历元间的坐标约束，历元间坐标信息无法利用。

因此本文探索一种自适应的顾及位置预测信息的实时动态精密单点定位算法。首先对未知参数进行分类，进行如下设置：$\mathit{\boldsymbol{\hat x}}$ = [x₁ x₂]^T，A = [A₁ A₂，x₁]为时变参数接收机钟差，x₂= [x y z ZTD N]^T，将坐标信息等作为时不变参数进行处理，其中坐标参数可通过逐历元附加不同先验约束进行时变特征转换。

然后利用前序历元(考虑速度预测的滞后性，设置固定窗口为前序两个历元)的坐标信息求解平均速度参数，并利用粗差探测技术进行数据异常剔除，进而结合当前历元坐标信息预测下一历元坐标参数值以及坐标参数的协因数阵：

$ {{\mathit{\boldsymbol{\bar X}}}_k} = \;{{\mathit{\boldsymbol{\hat x}}}_{k{\rm{-}}1}} + \Delta t*{\mathit{\boldsymbol{V}}_{k-1}} $

(4)

$ {\mathit{\boldsymbol{Q}}_{{{\mathit{\bar X}}_\mathit{k}}}} = {\mathit{\boldsymbol{Q}}_{{{\mathit{\hat x}}_{\mathit{k-1}}}}} + \Delta {t^2}*{\mathit{\boldsymbol{Q}}_{{V_{k-1}}}} $

(5)

式中，Q_Vk－1为前一历元速度协因数阵。

协因数阵为本历元内的相对数字标准^[6]，通过单位权方差进行绝对化，才能表征精度的绝对指标。对协因数阵在进行继承之前进行历元间的标准化统一：

$ {Q_{{{\mathit{\hat x}}_\mathit{k}}}} = \frac{{{\sigma _{k{\rm{序}}}}}}{{{\sigma _{k-1}}}}{Q_{{{\hat x}_{k-1}}}} $

(6)

式中，σ_k序为当前历元动态序贯求解所得单位权中误差，σ_k－1为前一历元求解所得单位权中误差。

对于各坐标分量的自适应因子确定，首先通过动态序贯最小二乘法进行当前历元定位求解(前序坐标参数不继承，仅继承前序模糊度和对流层参数信息)，得到当前历元的非继承坐标参数值 X_k序和其协因数信息 Q_xk序。利用当前历元坐标值与前序历元预测坐标值的差值和当前历元各坐标分量方差值(${\sigma _{k序}}\sqrt {{\mathit{\boldsymbol{Q}}_{{\mathit{x}_{\mathit{k}序}}}}} $ )进行对比确定各坐标分量自适应因子α，具体确定公式为：

$ \alpha = \frac{{\left| {{{\mathit{\boldsymbol{\bar X}}}_k}-{\mathit{\boldsymbol{X}}_{k序}}} \right|\;}}{{{\sigma _{k{\rm{序}}}}\mathit{\boldsymbol{\;}}\sqrt {{\mathit{\boldsymbol{Q}}_{{x_{k序}}}}} }} $

(7)

依据当前历元实时确定的自适应因子α，既可以有效利用前序历元的信息，又不会因约束过强而对当前结果造成影响。

最后将X_k直接替换${\mathit{\boldsymbol{\hat x}}}$_k－1代入式(3)解算，而 ${\mathit{\boldsymbol{Q}}_{{{\mathit{\bar x}}_\mathit{k}}}}$通过自适应因子α进行调整，即${\mathit{\boldsymbol{Q}}_{{{\mathit{\bar x}}_\mathit{k}}}} = \alpha {\mathit{\boldsymbol{Q}}_{{{\mathit{\bar x}}_\mathit{k}}}}$，得到顾及位置预测信息的精密单点定位估计值：

$ {{\mathit{\boldsymbol{\hat x}}}_k} = \;{{\mathit{\boldsymbol{\bar X}}}_k} + {\mathit{\boldsymbol{Q}}_{{{\mathit{\bar X}}_\mathit{k}}}}\mathit{\boldsymbol{A}}{\;^{\rm{T}}}{\mathit{\boldsymbol{P}}_k}({\mathit{\boldsymbol{l}}_k}{\rm{-}}{\mathit{\boldsymbol{A}}_k}{{\mathit{\boldsymbol{\bar X}}}_k})\; $

(8)

$ {\mathit{\boldsymbol{Q}}_{{{\mathit{\hat x}}_\mathit{k}}}} = {(\mathit{\boldsymbol{Q}}_{{{\mathit{\bar x}}_k}}^{-1} + A\;_k^{\rm{T}}{\mathit{\boldsymbol{P}}_k}{\mathit{\boldsymbol{A}}_k})^{{\rm{-}}1}} $

(9)

依据上述原理所得算法流程和处理策略如图 1和表 1所示。

2 算例分析

为验证该方法在单站动态形变监测中的有效性，于2016-08-27架设两台接收机(型号:和芯星通UR380)获取采样间隔1 s的观测数据，其中一台位于强制对中桩，作为基准站；另一台位于移动平台上，该平台通过相互垂直的导轨可以在东西与南北方向精确移动以此模拟缓变型形变。首先静态观测一段时间采集静态数据，然后利用实验平台进行动态模拟。模拟从第7 800个历元开始，过程如表 2所示。

实验分两个阶段：第一个阶段为未移动天线的静态观测阶段；第二个阶段为移动天线的动态观测阶段；两个阶段都采用实时动态PPP解算模式进行求解。

验证该算法的精度和有效性，设置如下4种解算方案：方案1，RTK算法；方案2，顾及位置预测信息PPP的算法；方案3，常规PPP序贯最小二乘算法；方案4，自适应滤波算法。

静态观测阶段流动站的坐标真值为事前30 min的静态基线解坐标值；动态阶段坐标无真值，仅以短距离(小于30 m)RTK定位结果为基准，对比方案2和方案3的符合精度。在§2.1中讨论其他3种方案与RTK的解算对比，在§2.2中详细分析方案2与方案3的结果。

2.1 1 s采样数据结果统计

图 2、图 3、图 4分别为E、N、U方向4种方案的残差序列，并依据接收机天线是否模拟位移将序列划分为静态与动态。表 3(单位m)为3种不同方案1 Hz采样率的精度统计结果。

从图 2可以看出，方案2比方案3、方案4的结果更加平稳，在动态缓变阶段，方案2与方案1结果符合精度更高。从图 3可以看出，N方向的定位精度显然优于E方向，因此方案2对N方向精度改善并不明显。从图 4可看出，U方向方案2的结果略优于方案3、方案4。表 4中，方案2的STD值优于方案3、方案4结果，且在动态阶段方案2平面RMS精度明显优于方案4结果，为cm级，可满足一般突发性形变监测的精度要求。

2.2 不同采样间隔数据结果统计

为了验证本方法计算不同采样间隔数据的差异性，将1 s数据降采样为5 s、10 s、15 s和30 s采样间隔进行解算，并将结果统计如下。

图 5中的静态阶段两种方案RMS精度相当；动态缓变阶段，对于高采样率数据方案2的解算结果优于方案3，随着采样间隔增大，二者精度相当。

U方向由于没有模拟动态，天线在其他方向上的移动都认定为动态阶段。可以看出，方案2在该方向的动态结果与E方向和N方向的效果相当，采样率越高，精度改善越明显，随着采样间隔增大精度趋于方案3。

表 4~表 6(单位m)给出了两种方案利用不同采样率数据定位得到的E、N、U 3个方向的STD和RMS值。

从表 4可以看出，随着采样间隔变大，方案2的静态STD值从4 cm增大到9 cm，当采样间隔为30 s时，两种方案精度值趋于一致；在静态阶段，方案2和方案3的RMS较为接近，方案2的结果略优。对于模拟动态阶段，方案2在高采样率时明显优于方案3结果，在采样间隔达到30 s时，两种方案结果趋近。从表 5也可以看出，数据采样率的提高并不能改进常规PPP序贯算法(方案3)的精度，而本文算法(方案2)随着采样率的提高定位精度有显著提升。

从表 5可以看出，N方向在静态阶段，STD精度优于5.5 cm，RMS优于8 cm；模拟动态阶段两种方案的RMS都小于7 cm，且对于高频数据方案2的精度比方案3提高约10%。

表 6中两种方案STD精度相当，方案2 RMS则要优于方案3，尤其是对高采样率数据的提高明显，约10%。

本文假定短时间内速度不变。在静态阶段，受测量噪声等影响速度预测值并不等于零，本文算法对精度改善并不明显；在动态缓变阶段，高频数据的解算精度显著提高，随着采样间隔的变大，坐标位置信息的预测值与实际值偏差较大，导致自适应因子α变大，解算结果等同于动态序贯最小二乘中将坐标参数视为时变参数，因此趋近于动态序贯最小二乘解算结果。利用观测信息自适应的调整预测信息协因数阵对当前历元解算的贡献，既符合实际情况，又有效地利用了历元间的坐标信息。本文方法对高频数据的解算具有较好的适用性，并对无基准站支持下的缓变型形变监测的适用性较强。随着采样间隔的加大，对定位精度的提高减弱，这与理论推导相符。

3 结语

本文探讨了一种顾及实时位置预测信息的动态精密单点定位算法，通过位置预测信息和自适应因子，兼顾了历元间坐标信息的约束，且不会因约束不合理而导致解算发散，提高了历元间参数信息的利用率，改善了实时动态定位精度。算例表明，载体在缓变形变发生阶段，所提算法E方向解算精度有较大提高，RMS从常规算法的10.6 cm提高到4 cm，高达60%；N方向的RMS在采样间隔为1 s时也有10%的提高；静态阶段该方法在U方向的RMS在采样间隔为1 s时也有10%左右的提高。随着采样间隔变大，位置预测信息误差增大，进而精度的改进效果减弱。因此，本文方法能有效改善高频动态数据处理的精度，对无基准站支持下的缓变型形变监测具有一定的实用价值。