GNSS技术及其误差
Published:
简介
位置服务，是GNSS项目的主要任务；然而经过多年的发展，GNSS技术已经广泛应用到了空间科学与地球科学的多个方面。
原理
GNSS卫星及LEO卫星以及地面台站，通过测量信号传播的时间及相位，可以获得信号源（比如GNSS卫星）与接收机（比如地面台站）间的传播时间与相位，进而知道信号源与接收机间距离的伪距与相位。 伪距，指并不是真是的发射机与接收机间的几何距离，而是由于信号传播过程中会受到多种因素影响，其传播时间对应的距离（乘以光速而获得）；相位则指通过记录某一初始时间的初相，记录累计相位观测数，以获得更高精度的距离信息。所以，两种观测数据，分别记为$P_i$与$L_i$，$i$表示不同频率对应的观测值，可以得到以下公式：
\[P^i_{j,k}=\rho_{j,k}+n^{P,i}_{j,k}\]
\[\phi^i L^i_{j,k}=\rho_{j,k}+\phi^i N^{\phi,i}_{j,k}+\phi^i n^{\phi,i}_{j,k}\] \[\rho_{j,k}=\sqrt{(X_j-X_k)^2+(Y_j-Y_k)^2+(Z_j-Z_k)^2}\]
其中 $\phi^i$ 表示频率为$f^i$的信号的波长，即$\frac{C}{f^i}$, $C$ 为光速; $\rho$, 即真实几何距离；$j, k$表示不同的发射机与接收机；$X_j, Y_j, Z_j, X_k, Y_k, Z_k$ 表示不同接收机与发射机的坐标;关于伪距与相位中的$n_{j,k}^{P,i}$、$N_{j,k}^{\phi,i}$、$n_{j,k}^{\phi,i}$ 分别表示伪距观测误差，相位观测误差及相位整周模糊度。
然而信号从等离子层穿过电离层（一般在 60km ～20000 km，即GNSS卫星（特别是GPS卫星）的卫星高度），中性大气（一般在2km～40km），最后直接或者通过反射传递到接收机天线——这一过程受到诸多因素的影响，为定位精度提高带了很多影响，同时也为误差所属研究方向提供了重要的空间观测信息。
注：这里特别指出30km～80km属于电离层与热层耦合区域，特别是电离层D层。
由于定位过程对于时间极度敏感，所以对与钟的稳定性要求极高。由于GNSS卫星相对用户数量较少，可以提供高精度的卫星钟以保证卫星信号在发射机传输时起算时间的准确性；而对于庞大的用户群体，比如全球及各国卫星监测站、个人手机及车辆，IoT产品以及无人驾驶车辆，为每个用户提供一个高精度的钟是不现实的。为此在一般定位过程中，一般假定卫星坐标与卫星时间是准确的，对地基用户或者数目庞大的用户，都是将信号接收时间当做未知量，联合未知的坐标，作为待估参数解算，从而得到准确的时间与坐标——这也是人们常说的，至少需要四颗卫星才可以实现定位，以及GNSS系统可以授时的原因。当然，卫星坐标虽然十分准确，但是卫星轨道需要地面观测站及注入站进行控制，卫星时间也受到卫星钟的稳定性影响，也会导致定位误差。卫星信号中的有实时卫星星历（其中包含钟差及钟漂信息，用以卫星钟校正），一般用于实时定位；同时各大数据分析中心以及IGS服务中心也有事后更加准确的卫星星历（具体信息可以在IGS中心查阅）以及其他产品。
当然，从另一面，如果知道了某个位置的坐标甚至时间，通过GNSS信号可以做什么呢？ 那就是获取“更精确”的$n_{j,k}^{P,i}$、$N_{j,k}^{\phi,i}$、$n_{j,k}^{\phi,i}$来获得中性大气及电离层信息。在考虑多种误差的情况下，对这三个量进行展开如下：
\[n^{P,i}_{j,k}=-C \delta t_j +C \delta t_k - \delta_{iono}^i +\delta_{trop} - C d^i_j + C d^i_k +m+\varepsilon\] \[\phi^i(N^{\phi,i}_{j,k}+n^{\phi,i}_{j,k})=-C \delta t_j +C \delta t_k + \delta_{iono}^i +\delta_{trop} +m+\varepsilon\]
其中$\delta t_j, \delta t_k$与原时对齐下的信号时间与接收机/发射机钟的差值，一般叫卫星钟差与接收机钟差；$d^i_j, d^i_k$ 分别是卫星及地面台站的硬件延迟；$m$ 是多路径误差；$\varepsilon$是其他误差，包括地面钟狭义相对论影响、相对论效应、地球自转、固体潮、天线相位偏差、区域形变、参考框架长时间项，观测噪音等等。由以上可以知道，卫星钟差$\delta t_j, \delta t_k$只与卫星/台站有关，$d^i_j, d^i_k$与卫星/台站和信号频率有关，$\delta_{iono}^i , \delta_{trop}, m$ 与卫星、台站及频率都有关系。通过地面参考框架，可以建立十分准确的大气及卫星模型以进行相关研究与应用，特别是改进精密定位在任和场景的实时解算。
当观测卫星，观测频段，观测台站，及观测时段足够多，便可以建立起对大地参考框架与区域性变形的监测以实现参考框架及地质形变的分析。不过夲节专讲PPP解算，主要讨论其中的误差、消除方式、及现在技术所达到的精度。
对流层
Thayer(1974)给出：
\[N=k_1(P_d/T)Z^{-1}_d+k_2(P_w/T)Z^{-1}_w+k_3(P_w/T^2)/Z^{-1}_w+k_4\frac{n_e}{f^2}\]，$P_d, P_w$分别是干空气和水汽的分压强；其中$Z^{-1}_d, Z^{-1}_w$分别为干空气和水汽的可压缩系数，
\[Z^{-1}_d=1+P_d[57.97 \times 10^{-8}(1+0.52/T)-9.4611\times 10^{-4}t/T^2]\approx 1\] \[Z^{-1}_w=1+1650(P_w/T^3)(1-0.01317 t+1.75 \times 10^{-5} t^2 +1.44 \times 10^{-6} t^3)\approx 1\]；$T, t$为绝对温度（K）和摄氏温度；或引入状态方程为
\[N=k_1 R_d \rho+k'_2(P_w/T)Z^{-1}_w+k_3(P_w/T^2)/Z^{-1}_w+k_4\frac{n_e}{f^2}\]，其中$k’_2=k_2-k_1(R_d/R_w)=16.48 K\times hPa^{-1}$即
\[N(s)=N_h(s)+N_w(s)+N_e(s)\]令$\Delta L_z, \Delta L_{zh}, \Delta L_{zw}$分别表示天顶总延迟，天顶静力学延迟和天顶湿延迟，则有
\[\Delta L_z=10^{-6}\times\int_{X_j}^{X_k}N(s))ds=\Delta L_{zh}+\Delta L_{zw}+\Delta L_{ze}\]\[\Delta L_{zh}=10^{-6}\times\int_{X_j}^{X_k}N_h(s)ds=10^{-6}\times\int_{X_j}^{X_k} k_1 R_d \rho ds\]
\[\Delta L_{zw}=10^{-6}\times\int_{X_j}^{X_k}N_w(s)ds=10^{-6}\times\int_{X_j}^{X_k} (k'_2\frac{P_w}{T}Z^{-1}_w+k_3\frac{P_w}{T^2}Z^{-1}_w) ds\]
\[\Delta L_{ze}=10^{-6}\times\int_{X_j}^{X_k}N_e(s)ds=10^{-6}\times\int_{X_j}^{X_k} k_4 \frac{n_e}{f^2} ds\]
其中$R_d=287(J\cdot kg^{-1}\cdot K^{-1})$；$R_w=461.5(J\cdot kg^{-1}\cdot K^{-1})$；$\rho=\rho_d+\rho_w$，即干空气及湿水气密度，$k_1=77.604\pm0.014 K \cdot hPa^{-1}$；$k_2=64.79\pm0.08 K \cdot hPa^{-1}$； $k_3=3.776\pm0.014K^2 \cdot hPa^{-1}$,$k_4=4.028\times 10^7$。
静力学延迟
Elgered(1991)通过静力学得到$d\rho=-\rho g d z$，即通过静力学延迟延迟积分得到：
\[\Delta L_{zh}=10^{-6}\times k_1 R_d\int_{X_j}^{X_k}\frac{ds}{g}=[0.0022768\pm0.0000024]\frac{P_0}{f(\varphi, H)}\] \[f(\varphi, H)=1-0.0026 cos 2\varphi -0.00028 H\]式中$P_0, \varphi， H$为GPS接收机处高度的气压(hPa)（应该是EC的湿压与干压之和），地理纬度，海拔高度（km），$f(\varphi, H)\times 10$为垂直大气 质量中心加速度，只与海拔与高度有关。由于地面气压可达0.5hPa，上式可以达到1mm精度。
天顶湿延迟
天顶湿延迟约0~50cm，由湿延迟公式可以得到：
\[\Delta L_{zw}=10^{-6}\times Z^{-1}_w\int_{X_j}^{X_k}[k'_2+k_3\frac{P_w/T^2}{P_w/T}]\frac{P_w}{T}dr=10^{-6}\times Z^{-1}_w[k'_2+k_3/T_m]\int_{X_j}^{X_k}\frac{P_w}{T}dr\] \[T_m=\frac{\int_{X_j}^{X_k}P_w/Tdr}{\int_{X_j}^{X_k}P_w/T^2dr}\]，一般，$T_m$可用地面温度估计，写成$T_m=a+b T_g$（美国，a=70.2,b=0.72；中国，a=44.6，b=0.81），$T_g$为地面温度。再次使用状态方程有：
\[\Delta L_{zw}=10^{-6}\times Z^{-1}_w[k'_2+k_3/T_m]R_w\int_{X_j}^{X_k}\rho_w dr=10^{-6}\times Z^{-1}_w[k'_2+k_3/T_m]R_w\times IWV\]令
\[\pi=\frac{10^6}{R_w[(k_3/T_m)+k'_2]}\] \[IWV=\pi\Delta L_{zw}\]可得可降水量$PWV$定义：
\[PWV=\frac{1}{\rho_{water}}\int_{X_j}^{X_k} \rho_w dr=\frac{IWV}{\rho_{water}}=\frac{\pi}{\rho_{water}}\Delta L_{zw}\] \[\frac{\pi}{\rho_{water}}=\frac{10^6}{\rho_{water}R_w[(k_3/T_m)+k'_2]}\approx 0.15\]其中$\rho_{water}$为水密度，最终结果单位与$L_{zw}$一致。
GNSS 解算方法
\[SPD=mf_h\times ZHD+mf_m \times ZWD + mf_g\times(G_N \cos(\alpha)+G_E \sin(\alpha)) + mf_i(\theta)\times STEC\] \[ZHD=\frac{0.0022768P_0}{1-0.0026 \cos{2\phi}-0.00028 h}\] \[ZWD=0.002277(\frac{1255}{T}+0.05)e\]其中$P_0, e, T, \phi, h$是压湿温，以及该处的纬度与高度，可以通过metRnx文件获得。
电离层
著名的Appleton-Hartree公式描述了电离层折射率指数：
\[n^2=1-\frac{X(U-X)}{U(U-X)-\frac{1}{2}Y^2sin^2\theta \pm\sqrt{\frac{1}{4}Y^2sin^4\theta+Y^2cos^2\theta(U-X)}}\]
其中：
\[X=\frac{\omega^2_p}{\omega^2}, X=\frac{\omega_c}{\omega}, Z=\frac{\nu_e}{\omega}, U=1-iZ\]式中：$\omega^2_p=e^2n_e/(m \varepsilon_0)$:等离子体角频率；$\omega_c=eB/m$:电子磁旋转角频率；$\nu_e$:电子有小碰撞频率；$\omega$:电磁波角频率；$n$:电离层折射率指数；$\theta$:地磁场与电磁波传播方向夹角；$e$:电子电荷；$n_e$:电子密度；$m$:电子质量；$\varepsilon$:自由空间介电常数；$B$:电地磁场强度。
对A-H方程进行级数展开
\[n \approx 1-\frac{1}{2}X\pm\frac{1}{2}XY|cos \theta|-\frac{1}{8}X^2-\frac{1}{4}XY^2(1-cos^2\theta)-i\frac{1}{2}XZ\]
由于GPS工作频率为$f_1$=1575.42 Mz，$f_2$=1227.60 Mz，$X«1$，$Y«1$，$Z«1$，电离层影响的一阶近似可以简化为
\[n=1-40.28\frac{n_e}{f^{i,2}}\]
上式也可以通过折射率方程得到。
同时由于电离层是一种色散介质，电磁波在电离层汇总传播时，各频率分量以各自不同的速度传播，最终导致电磁波波包能量传播速度（即群速度）小于光速，而载波相位的传播速度（即相速度大于光速），因而伪距的电离层延迟与载波相位的电离层延迟刚好符号相反，即：
\[n_p=1-40.28\frac{n_e}{f^2}, n_g=1+40.28\frac{n_e}{f^2}\] \[v_p=C/n_p=C(1+40.28\frac{n_e}{f^2}), v_g=C/n_g=C(1-40.28\frac{n_e}{f^2})\]
进而，GPS在电离层传播过程中，$\delta_{iono}$有以下关系
\[{\delta_{iono}}_{j,k}^{P,i}=\int_{\Delta t}\nu_gdt-\rho_{j,k}=C \Delta t-\rho_{j,k}-\frac{40.28}{f^{i,2}}\int_{\Delta t} n_e dt=n^P-\frac{40.28}{f^{i,2}}\int_{\Delta t} n_e dt\] \[{\delta_{iono}}_{j,k}^{\phi,i}=\int_{\Delta t}\nu_pdt-\rho_{j,k}=C \Delta t-\rho_{j,k}+\frac{40.28}{f^{i,2}}\int_{\Delta t} n_e dt=n^\phi+\frac{40.28}{f^{i,2}}\int_{\Delta t} n_e dt\]
令
\[TEC_{j,k}=\int^{X_j}_{X_k} n_e ds\]
则由电离层对载波及相位的影响通过观察电离层误差方程，可以通过双频线性组合的方式获得：
\[TEC^P_{j,k}=\frac{f^{1,2}f^{2,2}}{40.28(f^{1,2}-f^{2,2})}[n^{P,2}_{j,k}-n^{P,1}_{j,k} + C \delta d_j - C \delta d_k - \delta m_{j,k} -\delta \varepsilon_{j,k} ]\] \[TEC^\phi_{j,k}=\frac{f^{1,2}f^{2,2}}{40.28(f^{1,2}-f^{2,2})}[\phi^1(N^{\phi,1}_{j,k}+n^{\phi,1}_{j,k})-\phi^2(N^{\phi,2}_{j,k}-n^{\phi,2}_{j,k}) + \delta m_{j,k} + \delta \varepsilon_{j,k}]\] \[\delta d_j=d^1_j-d^2_j\] \[\delta d_k=d^1_k-d^2_k\] \[\delta m_{j,k}=m^1_{j,k}-m^2_{j,k}\] \[\delta \varepsilon_{j,k}=\varepsilon^1_{j,k}-\varepsilon^2_{j,k}\]
即
\[TEC^P_{j,k}=\frac{f^{1,2}f^{2,2}}{40.28(f^{1,2}-f^{2,2})}[P^2_{j,k}-P^2_{j,k} + C \delta d_j - C \delta d_k - \delta m_{j,k} -\delta \varepsilon_{j,k} ]\] \[TEC^\phi_{j,k}=\frac{f^{1,2}f^{2,2}}{40.28(f^{1,2}-f^{2,2})}[\phi^1 L^1_{j,k}-\phi^2 L^2_{j,k} -D^{\phi}_{j,k} + \delta m_{j,k} + \delta \varepsilon_{j,k}]\] \[D^{\phi}_{j,k}=\phi^1 N^{\phi,1}_{j,k}-\phi^2 N^{\phi,2}_{j,k}\]一般，\(D^{\phi}_{j,k}\)为模糊度；$\delta d_j$，$\delta d_k$分别称为接收机及卫星硬件码偏差；$\delta m_{j,k}$称为多路径效应残差，可通过半天球模型降低；$\varepsilon_{j,k}$为观测误差及其他误差。
由于相位只能获得相对TEC，但是精度高（百分之几的TECU）；而伪距可以获得绝对TEC，但是精度差。通过两者联合求解，可以得到高精度的绝对TEC。具体方法如下：
- 计算一次连续弧段中伪距与相位TEC的偏差为：
- 校正相位TEC：
其他误差
相对论效应
地面钟由于地球自转会产生钟频变化
\[\delta f_R=-\frac{V^2_R}{2c^2}\times f\]。由于$\delta f_R$极小，所以可以把地球当做半径为$R$的球：
\[V_R=V_0\cdot cos \phi， V_0=464m/s\] \[\delta f_{R_0} \approx 0.012(f\cdot 10^{-10})\]，其结果极小。
而引力延迟引起的误差
\(\Delta D_g=\frac{2\nu}{c^2}ln \frac{\rho_j+\rho_k+\rho_{j,k}}{\rho_j+\rho_k-\rho_{j,k}}\)c
， 其中
\[\rho_j=\sqrt{X_j^2+Y_j^2+Z_j^2},\rho_k=\sqrt{X_k^2+Y_k^2+Z_k^2}\]。当站间距离达到1000~3000km时，引力延迟误差在1.3～3.6mm之间。
地球自转
固体潮
由于地球非刚体，在日月引力作用下，地球会产生周期性的形变，可能使地面店在垂直方向产生80km的变化。
设日月二阶引潮位为$V_2$；$h$为第一勒夫数，0.6090；$l$为第二勒夫数，0.0852，则固体潮对测站($R, \theta, \lambda$，分别是向径、余纬及经度)产生的变化为：
\[\begin{cases}\delta R=h\frac{V_2}{g}=\frac{h}{2g}r\frac{V_2}{\alpha R}\\\delta \theta=\frac{l}{g}\frac{\alpha V_2}{\alpha \theta}\\\delta \lambda=\frac{l}{g}\frac{\alpha V_2}{sin \theta \alpha \lambda}\end{cases}\]即
\[u=\begin{bmatrix}\delta R\\\delta \theta \\\delta \lambda\end{bmatrix}=\frac{m_w}{M}\frac{\rho_w^4}{\rho_j^3}\{3l(\rho_w^0\cdot \rho_j^0)\rho_w^0+[3(\frac{h}{2}-l)(\rho_w^0\cdot \rho_j^0)^2-\frac{h}{2}]\rho_j^0\}\]
其中$\rho_w^0$为摄动天体到地心单位矢矩，$\rho_j^0$为地面站到地心单位矢矩，$M$为地球质量，$m_w$是摄动体质量。形变与作用力在时间上滞后$\varphi$（一般取2.5°），所以需要乘以一个滞后因子$L$：
\[L=\begin{bmatrix}cos \varphi&sin \varphi& 0\\-sin \varphi &cos \varphi &0\\0&0&1 \end{bmatrix}\]从测站中的形变转换到协议坐标中，还需要乘以$R$:
\[R=\begin{bmatrix}cos \theta cos \lambda & cos \theta sin \lambda & -sin \theta\\-sin \lambda & cos \theta &0\\sin \theta cos \lambda &sin \theta sin \lambda & cos \theta\end{bmatrix}\]即
\[u_{cmt}=R\cdot L\cdot u\]如果不对测站坐标改正，而对距离观测值$\rho_{j,k}$改正，则
\[\delta \rho=\rho^0\cdot u_{cmt}=\frac{r_s-R}{|r_s-R|} u_{cmt}\]卫星及接收机天线相位偏差
卫星天线相位偏差会影响到卫星钟差的解算。如果已知偏差向量$\Delta$，若星固坐标为$(e_x, e_y, e_z)$，则偏差向量为
\[\Delta_{cmt}=\rho^0_{j,k} \cdot \Delta_{cmt}=\frac{r_s-R}{|r_s-R|}\cdot\Delta_{cmt}\]。
而接收机天线相位中心变化一般通过改正实现，由于GNSS接收机天线的平均相位中心与天线的几何中心一般不重合，其偏差向量在高精度单点定位是需要采用归心改正。而接收机的瞬时相位中心与平均相位中心间的偏差，随着卫星高度角变化，主要影响高程。