GNSS精密单点定位（PPP）基本原理（进阶篇）

上节介绍了精密单点定位的基本原理，本文继续在精密单点定位的基础上进行更深层次的介绍。

一、精密单点定位的函数模型

上节说过，在精密单点定位之前，也有一种绝对定位技术，那就是伪距单点定位。伪距单点定位靠的伪距进行单点定位，但是伪距的精度较差，主要原因是其噪声较大，且多径效应明显。而精密单点定位的精密就精密在其一，轨道的精密、钟差的精密。即使用的是IGS精密星历和精密钟差。二、就是观测值的精密。这里的观测值就是比伪距精度更高的载波相位观测值。三、函数模型的精密。伪距单点定位中的误差模型都已经在书本上介绍了，如大气误差模型。但是那种模型的精度不能满足精密单点定位的需求，因此在精密单点定位中，需要对大气误差延迟进行估计和处理。如对流层延迟一般都是赋予参数进行估计，而电离层延迟则需进行组合消除，如上节介绍的无电离层组合。对于非组合模型，也需要对电离层进行估计。因此精密单点定位的函数模型可以表述如下：
[ P 1 − ρ 1 − D P 1 L 1 − ρ 1 − D L 1 . . . P n − ρ n − D p n L n − ρ n − D L n ] = [ x 0 − x 1 ρ 0 y 0 − y 1 ρ 0 z 0 − z 1 ρ 0 1 M w e t 1 0 . . . 0 x 0 − x 1 ρ 0 y 0 − y 1 ρ 0 z 0 − z 1 ρ 0 1 M w e t 1 1 . . . 0 . . . x 0 − x n ρ 0 y 0 − y n ρ 0 z 0 − z n ρ 0 1 M w e t n 0 . . . 0 x 0 − x n ρ 0 y 0 − y n ρ 0 z 0 − z n ρ 0 1 M w e t n 0 . . . 1 ] [ d x d y d z c ⋅ δ t d Z T D w B 1 . . . B n ] \begin{bmatrix} P_1-\rho_1-D_{P1} \\ L_1-\rho_1-D_{L1} \\ ...\\ P_n-\rho_n-D_{pn} \\ L_n -\rho_n-D_{Ln}\end{bmatrix}=\begin{bmatrix} \frac{x_0-x^1}{\rho_0} & \frac{y_0-y^1}{\rho_0} & \frac{z_0-z^1}{\rho_0} & 1 & M^1_{wet} & 0 & ... & 0 \\ \frac{x_0-x^1}{\rho_0} & \frac{y_0-y^1}{\rho_0} & \frac{z_0-z^1}{\rho_0} & 1 & M^1_{wet} & 1 & ... & 0 \\ ...\\ \frac{x_0-x^n}{\rho_0} & \frac{y_0-y^n}{\rho_0} & \frac{z_0-z^n}{\rho_0} & 1 & M^n_{wet} & 0 & ... & 0 \\ \frac{x_0-x^n}{\rho_0} & \frac{y_0-y^n}{\rho_0} & \frac{z_0-z^n}{\rho_0} & 1 & M^n_{wet} & 0 & ... & 1 \end{bmatrix} \begin{bmatrix} dx \\ dy \\ dz \\ c·\delta t \\ dZTD_w \\ B_1 \\ ... \\ B_n\end{bmatrix} ⎣⎢⎢⎢⎢⎡P1−ρ1−DP1L1−ρ1−DL1...Pn−ρn−DpnLn−ρn−DLn⎦⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎡ρ0x0−x1ρ0x0−x1...ρ0x0−xnρ0x0−xnρ0y0−y1ρ0y0−y1ρ0y0−ynρ0y0−ynρ0z0−z1ρ0z0−z1ρ0z0−znρ0z0−zn1111Mwet1Mwet1MwetnMwetn0100............0001⎦⎥⎥⎥⎥⎥⎥⎤⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡dxdydzc⋅δtdZTDwB1...Bn⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤
上式就是精密单点定位的基本函数模型。即为n个卫星的伪距和载波相位观测方程组成的 y − b = A X y-b=AX y−b=AX的形式。这种形式也是测绘领域最小二乘的基本观测方程式。

二、精密单点定位的随机模型

仅仅依赖上式是无法将待求参数估计出来的，因为在上式中不同观测方程的精度并不一致，因此需要对观测方程定权。因为在最小二乘系统中是需要有权矩阵的。观测方程噪声中主要由各个模型的误差组成：
σ 2 = σ U E R E 2 = σ e p h 2 + σ c l k 2 + σ i o n 2 + σ t r o p 2 + σ m p 2 + σ n o i s e 2 \sigma^2=\sigma^2_{UERE}=\sigma^2_{eph}+\sigma^2_{clk}+\sigma^2_{ion}+\sigma^2_{trop}+\sigma^2_{mp}+\sigma^2_{noise} σ2=σUERE2=σeph2+σclk2+σion2+σtrop2+σmp2+σnoise2
其中 σ 2 \sigma^2 σ2即为观测方程误差， σ U E R E 2 \sigma^2_{UERE} σUERE2即为用户等效距离误差，后边各项的含义在下标中已经给出。其中 σ i o n 2 \sigma^2_{ion} σion2在使用无电离层组合时可不用考虑； σ t r o p 2 \sigma^2_{trop} σtrop2在估计对流层湿延迟参数时不用考虑。
一般而言，会将随机模型分为与高度角有关和与高度角无关的部分：
σ 2 = σ S I S U R E 2 + σ e l e 2 = σ e p h 2 + σ c l k 2 + σ 2 ( e l e ) \sigma^2=\sigma^2_{SISURE}+\sigma^2_{ele}=\sigma^2_{eph}+\sigma^2_{clk}+\sigma^2(ele) σ2=σSISURE2+σele2=σeph2+σclk2+σ2(ele)
其中 σ S I S U R E 2 \sigma^2_{SISURE} σSISURE2为GNSS空间信号精度，主要由后面的星历误差和种误差组成，其余的即为与高度角相关的部分，一般就是按照高度角定权。即
{ σ ( e l e ) = σ 0 , e l e > 30 σ ( e l e ) = σ 0 2 s i n ( e l e ) , e l s e \begin{cases} \sigma(ele)=\sigma_0, ele > 30 \\ \sigma(ele)=\frac{\sigma_0}{2sin(ele)}, else \end{cases} {σ(ele)=σ0,ele>30σ(ele)=2sin(ele)σ0,else
而权矩阵即为
W = R − 1 = [ 1 σ 1 2 ⋱ 1 σ n 2 ] W=R^{-1}=\begin{bmatrix} \frac{1}{\sigma^2_1} \\ & \ddots & \\ & & \frac{1}{\sigma^2_n} \end{bmatrix} W=R−1=⎣⎢⎡σ121⋱σn21⎦⎥⎤

三、精密单点定位的参数估计方法

对于单站来说，精密单点定位的参数估计方法一般有两种，即最小二乘和卡尔曼滤波方法。最小二乘就是将所有的观测方程形成法方程求解，而卡尔曼滤波估计当前时刻的参数时仅需知道上一时刻的状态参数即可，因此卡尔曼滤波相比最小二乘来说，有较小的耗存。尤其是在动态PPP求解的时候，最小二乘就需要采用序贯平差，但是序贯平差的法方程矩阵会越来越大，因此实际求解的时候一般会采用参数消去的序贯平差求解。而卡尔曼滤波在动态系统有天然的优势性，不管是在定轨、动态定位还是组合导航中都应用广泛。且针对不同的系统，又包括扩展卡尔曼滤波、自适应卡尔曼滤波、UD分解滤波等等。本节不对其展开，仅对这两种方法的最简单形式叙述。

（1）最小二乘
根据上节介绍的PPP的函数模型和随机模型，建立的GNSS定位的数学模型如下：
y = G x + ϵ ， R y=Gx+\epsilon，R y=Gx+ϵ，R
ϵ \epsilon ϵ满足随机模型。即 E [ ϵ ] = 0 E[\epsilon]=0 E[ϵ]=0及 R = E [ ϵ ϵ T ] R=E[\epsilon \epsilon^T] R=E[ϵϵT].
接下来就是我们熟知的组成法方程、法方程求逆、求解待估参数。这里直接给出最后的结果。
x ^ = ( G T R − 1 G ) − 1 G T R − 1 y P = ( G T R − 1 G ) − 1 \hat x=(G^TR^{-1}G)^{-1}G^TR^{-1}y \\ P=(G^TR^{-1}G)^{-1} x^=(GTR−1G)−1GTR−1yP=(GTR−1G)−1
在动态系统中，序贯平差会增加观测方程的个数，也就导致了法方程矩阵的维数越来越多，求逆就变成了一个难题。

对于维数较大的法方程矩阵，一般可以采用Cholesky分解法又称为平方根法求逆。这种方法可由读者自行学习。

那么序贯平差是在做一件什么事呢？

假如有两组观测方程，他们有全部相同的未知参数，也可以有部分相同的未知参数。
y 1 = G 1 x + ϵ 1 , R 1 y 2 = G 2 x + ϵ 2 , R 2 y_1=G_1 x+\epsilon_1,R_1 \\ y_2 = G_2 x + \epsilon_2,R_2 y1=G1x+ϵ1,R1y2=G2x+ϵ2,R2
对这两组观测方程也可以组成法方程矩阵求解，结果为
x ^ = [ G 1 T R 1 − 1 G 1 + G 2 T R 2 − 1 G 2 ] − 1 [ G 1 T R 1 − 1 y 1 + G 2 T R 2 − 1 y 2 ] P = [ G 1 T R 1 − 1 G 1 + G 2 T R 2 − 1 G 2 ] − 1 \hat x= [G^T_1R_1^{-1}G_1+G_2^TR^{-1}_2G_2]^{-1}[G^T_1R^{-1}_1y_1+G^T_2R^{-1}_2y_2] \\ P = [G^T_1R^{-1}_1G_1+G^T_2R^{-1}_2G_2]^{-1} x^=[G1TR1−1G1+G2TR2−1G2]−1[G1TR1−1y1+G2TR2−1y2]P=[G1TR1−1G1+G2TR2−1G2]−1
上式为两组观测方程的解。可以对比一下与经典最小二乘的解。
序贯最小二乘的精髓在于迭代求解，因为是针对与历元解，所以需要对其进行迭代，那么观测方程会越来越多，即会有第三组、第四组…迭代就是每增加一组观测方程，就在前面几组计算的基础上求解。
x ^ 1 = P 1 ⋅ [ G 1 T R 1 − 1 y 1 ] P 1 = [ G 1 T R 1 − 1 G 1 ] − 1 x ^ 2 = P 2 ⋅ [ G 1 T R 1 − 1 y 1 + G T ] P 2 = [ P 1 − 1 + G 2 T R 2 − 1 G 2 ] − 1 \hat x_1=P_1·[G^T_1R^{-1}_1y_1] \\ P_1=[G^T_1R^{-1}_1G_1]^{-1} \\ \hat x_2=P_2 ·[G^T_1R^{-1}_1y_1+G^T] \\ P_2 = [P^{-1}_1+G^T_2R^{-1}_2G_2]^{-1} x^1=P1⋅[G1TR1−1y1]P1=[G1TR1−1G1]−1x^2=P2⋅[G1TR1−1y1+GT]P2=[P1−1+G2TR2−1G2]−1
当增加第三组的时候，那么第i组的解的法方程矩阵就在第i-1组的基础上在增加一个 G i T R i − 1 G i G^T_iR^{-1}_iG_i GiTRi−1Gi。
对于静态系统，也可以采用序贯平差计算，但是此时序贯平差在每个历元的解都不同，只有最后一个历元的解才是真解。
序贯平差的优势是在迭代的过程中可以消除不再感兴趣的参数。如接收机钟差等。在动态定位中，上个历元的位置参数就不再感兴趣，可以消去。当发生周跳时，上个历元的模糊度参数就要消去。
（2）卡尔曼滤波
卡尔曼滤波在GNSS定位定轨中非常常见，因此想要学习GNSS定位定轨，卡尔曼滤波的知识也不可或缺。其实卡尔曼滤波和序贯平差很相似。
对于n-1个历元和n历元，参数和协方差的预测模型如下：
x ^ − ( n ) = Φ ( n − 1 ) x ^ − ( n − 1 ) P − x ( n ) ^ = Φ ( n − 1 ) P − x ( n − 1 ) ^ Φ T ( n − 1 ) + Q ( n − 1 ) \hat x^-(n)=\Phi (n-1)\hat x^-(n-1) \\ {P^-}_{\hat{x(n)}}=\Phi (n-1){P^-}_{\hat{x(n-1)}}\Phi^T(n-1)+Q(n-1) x^−(n)=Φ(n−1)x^−(n−1)P−x(n)^=Φ(n−1)P−x(n−1)^ΦT(n−1)+Q(n−1)
更新方程如下：
P x ( n ) ^ = [ I − K ( n ) G ( n ) ] P − x ( n ) ^ x ^ ( n ) = x ^ − ( n ) + K ( n ) [ y ( n ) − G ( n ) x ^ − ( n ) ] {P}_{\hat{x(n)}}=[I-K(n)G(n)]{P^-}_{\hat{x(n)}} \\ \hat x(n)=\hat x ^-(n)+K(n)[y(n)-G(n)\hat x^-(n)] Px(n)^=[I−K(n)G(n)]P−x(n)^x^(n)=x^−(n)+K(n)[y(n)−G(n)x^−(n)]
K(n)称为卡尔曼增益：
K ( n ) = P − x ( n ) ^ G T ( n ) [ G ( n ) P − x ( n ) ^ G ( n ) T + R ( n ) ] − 1 K(n)={P^-}_{\hat{x(n)}}G^T(n)[G(n){P^-}_{\hat{x(n)}}G(n)^T+R(n)]^{-1} K(n)=P−x(n)^GT(n)[G(n)P−x(n)^G(n)T+R(n)]−1
上面5个公式就是卡尔曼的经典5公式。
对于静态定位，坐标参数不变，钟差参数为随机噪声模型，对流层参数一般为随机游走模型，未发生周跳时，其模糊度参数不变，故其状态转移矩阵 Φ \Phi Φ包含坐标、接收机、对流层延迟和模糊度参数的状态转移。
Φ = [ 1 1 1 0 1 1 ] Q = [ 0 0 0 σ δ t 2 σ t r o p 2 0 ] \Phi=\begin{bmatrix} 1 & \\ & 1 & \\ & & 1 & \\ & & & 0 & \\ & & & & 1 &\\ & & & & & 1 \end{bmatrix} \\ Q=\begin{bmatrix} 0 & \\ & 0 & \\ & & 0 & \\ & & & \sigma^2_{\delta t} & \\ & & & & \sigma^2_{trop} \\ & & & & & 0 \end{bmatrix} Φ=⎣⎢⎢⎢⎢⎢⎢⎡111011⎦⎥⎥⎥⎥⎥⎥⎤Q=⎣⎢⎢⎢⎢⎢⎢⎡000σδt2σtrop20⎦⎥⎥⎥⎥⎥⎥⎤
当发生周跳时， Φ \Phi Φ中右下角为0.

对于动态定位，又分为运动速度已知和未知的情况。若速度未知，则
Φ = [ 0 0 0 0 1 1 ] Q = [ σ d x 2 σ d y 2 σ d z 2 σ δ t 2 σ t r o p 2 0 ] \Phi=\begin{bmatrix} 0 & \\ & 0 & \\ & & 0 & \\ & & & 0 & \\ & & & & 1 &\\ & & & & & 1 \end{bmatrix} \\ Q=\begin{bmatrix} \sigma^2_{dx} & \\ & \sigma^2_{dy} & \\ & & \sigma^2_{dz} & \\ & & & \sigma^2_{\delta t} & \\ & & & & \sigma^2_{trop} \\ & & & & & 0 \end{bmatrix} Φ=⎣⎢⎢⎢⎢⎢⎢⎡000011⎦⎥⎥⎥⎥⎥⎥⎤Q=⎣⎢⎢⎢⎢⎢⎢⎡σdx2σdy2σdz2σδt2σtrop20⎦⎥⎥⎥⎥⎥⎥⎤
其中 [ σ d x 2 σ d y 2 σ d z 2 ] \begin{bmatrix} \sigma^2_{dx} & \sigma^2_{dy} & \sigma^2_{dz} \end{bmatrix} [σdx2σdy2σdz2]为坐标噪声，当高速运动时，可设置为10Km。

若速度已知，坐标参数可变为随机游走。
Φ = [ 1 1 1 0 1 1 ] Q = [ Q d x ′ Δ t Q d y ′ Δ t Q d z ′ Δ t σ δ t 2 σ t r o p 2 0 ] \Phi=\begin{bmatrix} 1 & \\ & 1 & \\ & & 1 & \\ & & & 0 & \\ & & & & 1 &\\ & & & & & 1 \end{bmatrix} \\ Q=\begin{bmatrix} Q^{'}_{dx}\Delta t & \\ &Q^{'}_{dy}\Delta t& \\ & &Q^{'}_{dz}\Delta t & \\ & & & \sigma^2_{\delta t} & \\ & & & & \sigma^2_{trop} \\ & & & & & 0 \end{bmatrix} Φ=⎣⎢⎢⎢⎢⎢⎢⎡111011⎦⎥⎥⎥⎥⎥⎥⎤Q=⎣⎢⎢⎢⎢⎢⎢⎡Qdx′ΔtQdy′ΔtQdz′Δtσδt2σtrop20⎦⎥⎥⎥⎥⎥⎥⎤
其中 [ Q d x ′ Q d y ′ Q d z ′ ] \begin{bmatrix} Q^{'}_{dx} & Q^{'}_{dy} & Q^{'}_{dz} \end{bmatrix} [Qdx′Qdy′Qdz′]为速度噪声。

不管静态还是动态定位，其先验信息都可以设置为
Q 0 = [ 1 2 k m 2 1 2 k m 2 1 2 k m 2 30 0 2 k m 2 0. 1 2 m 2 2 0 2 m 2 ] Q_0=\begin{bmatrix} 1^2km^2& \\ &1^2km^2& \\ & &1^2km^2 & \\ & & & 300^2km^2& \\ & & & & 0.1^2m^2 \\ & & & & & 20^2m^2 \end{bmatrix} Q0=⎣⎢⎢⎢⎢⎢⎢⎡12km212km212km23002km20.12m2202m2⎦⎥⎥⎥⎥⎥⎥⎤

本节介绍了PPP的函数模型、随机模型以及其两种参数估计方法。希望通过本篇，能对初入GNSS的小白提供一点研究思路。计划下节介绍一些多系统PPP的基本原理。

参考文献：GNSS高精度定位原理（董大南）强推此书，希望研究GNSS的小白们人手一本。