（一）GNSS定位中的电离层延迟估计

写在前面：

今天开一个板块记录一下自己在GNSS定位（PPP，PPP-RTK）方面的学习过程和心得想法。因为之前是做电离层预报，从大三开始接触算起也两年有余了，或多或少有一些成果。后面想学习一下精密单点定位以及多源融合，可以的话再学学自动驾驶相关的框架和算法。这个板块的内容争取更新频率高一些，更的快说明学的快☺

因为之前和电离层接触的比较多，所以从定位中的电离层延迟开始写起。

GNSS定位中为什么要对电离层建模？

电离层是日-地空间环境的重要组成部分，会对人类的生产和生活产生影响。当全球卫星导航系统（Global Navigation Satellite System，GNSS）发射的电磁波信号在电离层中传播时，会产生几米到几十米的误差，是制约GNSS用户（尤其单频GNSS用户）高精度定位的主要误差源之一。电离层总电子含量（Total Electron Content，TEC）是描述电离层变化特性的关键参量，构建实时电离层TEC模型可为实时导航定位用户提供电离层延迟改正，加快精密单点定位（Precise Point Positioning，PPP）收敛速度，实现对空间天气的精准监测。
与传统电离层探测技术相比，地基GNSS探测技术具有探测精度高、覆盖范围广、可实现24h全天候连续探测并可提供电离层底部到GNSS卫星之间TEC值等优点，已被广泛应用于电离层TEC模型的构建。同时，为满足GNSS用户对高精度电离层延迟的实时改正需求并实现对电离层的实时精准监测，国际GNSS服务组织（International GNSS Service，IGS）电离层分析中心于近些年开始关注实时电离层TEC模型的构建，并取得了一些初步结果。

电离层概况

电离层是地球大气的一个电离区域，通常指从离地面约60 km到约2000 km磁顶层的地球高层大气空域，因该区域存在大量的自由电子和离子而被命名为电离层。电离层中的电子密度(又称电子浓度)是指单位体积的自由电子数。在不同的高度上，伴随着大气成分、大气密度、太阳辐射的变化以及太阳活动、地磁活动等的影响下，电子密度会呈现不同的变化规律。而根据电子密度的大小不同，在高度方向电离层从低往高可以依次分为D层，E层和F层，其中F层又分为F1和F2。电离层分层结构如图1。

D层是电离层最接近地面的一层，大致位于地球表面60-90km处，电子密度大致为10²-10⁴/cm3。D层中电子产生主要是受到太阳辐射的影响导致中性分子和原子电离。该层电子密度在午后达到最大，太阳降落后逐渐降低，在夜间达到最小甚至消失。同时该层的电子密度也会出现季节变化，夏季的电子密度大于冬季。在太阳活动高年，该层的电子密度也明显的高于太阳活动低年。由于在D层中，中性离子和自由电子的碰撞比较剧烈，导致对电磁波的影响较大。
E层属于电离层的中间区域，大约位于地面90-150km处，电子密度大致为10³-10⁵/cm3 。E层主要受到X射线的影响，也存在明显的昼夜变化、季节变化和太阳活动周期等的变化。在日间E层电子密度较大，夜间则会急剧下降，但不会消失。总体而言E层的结构比较稳定，但有时也会出现一个高度约为100km~200km左右的不规则结构，其密度高于上下区域的两倍，称作偶发E层。
F层属于电离层的顶部，大约位于地面150-500km处，是电子密度最大的一层，大致为10⁵-10⁶/cm3。该层的形成主要受到20-90nm的太阳辐射与地球磁场的作用，处于漂移平衡的状态。F层的电子密度为整个电离层最高，在较大尺度的电离层变化中，往往与该层的关系最为密切。F层进一步可以分为F1层与F2层，F1层的高度约为150-210km，F2层的高度约为210km-600km，电子密度的峰值高度一般在250km~400km处。其中F2层的变化最为剧烈，是电离层的主要部分，其变化对整个电离层的影响也最大，也是现如今研究的重点区域。F1层的峰值经常不出现，在夏季中纬度F1层才会出现，并且在夜间F1层通常会消失。
由于电离层受到太阳活动的影响较大，因此电离层在周期上往往会出现与太阳周期类似的变化，如昼夜日周期变化，季节性周期变化，年周期变化，11年周期变化等。除此之外，电离层还受到大气分子浓度，高层大气风速，地磁活动等因素的影响，导致电离层在周期上出现不规则的变化特性。

电离层延迟

电磁波在电离层中传播的相速度Vp=cnp{V_{\rm{p}}} = {c \over {{n_p}}}Vp=npc，其中c为光速，相折射率
np=1−K1Nef−2−K2Ne(H0cos⁡θ)f−3−K3Ne3f−4{n_p} = 1 - {K_1}Ne{f^{ - 2}} - {K_2}Ne({H_0}\cos \theta ){f^{ - 3}} - {K_3}N{e^3}{f^{ - 4}}np=1−K1Nef−2−K2Ne(H0cosθ)f−3−K3Ne3f−4其中K1=e28πε0m{K_1} = {{{e^2}} \over {8\pi {\varepsilon _0}m}}K1=8πε0me2,K2=μ0e316π3ε0m2{K_2} = {{{\mu _0}{e^3}} \over {16{\pi ^3}{\varepsilon _0}{m^2}}}K2=16π3ε0m2μ0e3,K3=e4128π4ε03m2{K_3} = {{{e^4}} \over {128{\pi ^4}\varepsilon _0^3{m^2}}}K3=128π4ε03m2e4，后两项的值非常小，均小于10^-9，可以忽略不计。H0{H_0}H0为地磁场的场强；Ne为电子密度，m为电子的质量，e为电子所带的电荷值，ε0{\varepsilon _0}ε0为真空中的介电系数，θ\thetaθ为地磁场方向和电磁波信号传播方向的夹角，f为电磁波信号的频率。
因此， np≈1−K1Nef−2=1−40.3Nef2{n_p} \approx 1 - {K_1}Ne{f^{ - 2}} = 1 - 40.3{{Ne} \over {{f^2}}}np≈1−K1Nef−2=1−40.3f2Ne
代入可得 Vp=cnp=c1−40.3Nef2=c(1+40.3Nef2){V_p} = {c \over {{n_p}}} = {c \over {1 - 40.3{{Ne} \over {{f^2}}}}} = c(1 + 40.3{{Ne} \over {{f^2}}})Vp=npc=1−40.3f2Nec=c(1+40.3f2Ne)对于电磁波的群速度VG=cnG=c1+40.3Nef2=c(1−40.3Nef2){V_G} = {c \over {{n_G}}} = {c \over {1 + 40.3{{Ne} \over {{f^2}}}}} = c(1 - 40.3{{Ne} \over {{f^2}}})VG=nGc=1+40.3f2Nec=c(1−40.3f2Ne)对于这两个速度这里简单提一下，相速度可以理解为等相位面的传播速度，群速度则是等振幅面的传播速度。所以测距码的速度是群速度，而载波相位的速度是相速度。相速度大于光速切勿理解为波的速度超过光速，物理上没有速度能够超过光速，这里只是波的包络面形成的相位面的传播速度，而不是实际的粒子波传播速度。
对于卫星到接收机的几何距离有
ρ=∫Δt′VGdt=∫Δt′(c−c⋅40.3Nef2)dt=c⋅Δt′−40.3f2∫Δt′c⋅Nedt\rho = \int\limits_{\Delta t'} {{V_G}dt = } \int\limits_{\Delta t'} {(c - c \cdot 40.3{{Ne} \over {{f^2}}})} dt = c \cdot \Delta t' - {{40.3} \over {{f^2}}}\int\limits_{\Delta t'} {c \cdot Nedt} ρ=Δt′∫VGdt=Δt′∫(c−c⋅40.3f2Ne)dt=c⋅Δt′−f240.3Δt′∫c⋅Nedtρ=ρ′−40.3f2∫sNeds\rho = \rho ' - {{40.3} \over {{f^2}}}\int\limits_s {Neds} ρ=ρ′−f240.3s∫Neds所以使用测距码的电离层延迟改正项为(Vion)G=−40.3f2∫sNeds{({V_{ion}})_G} = - {{40.3} \over {{f^2}}}\int\limits_s {Neds} (Vion)G=−f240.3s∫Neds载波相位的电离层改正项也同理，所以有一个重要的结论：测码伪距观测值和载波相位观测值的电离层延迟改正大小相同，符号相反。

总电子含量TEC

由上述介绍电离层概况可知，电子密度Ne和大气高度H有一定的关系。所以，引入总电子含量TEC（Total Electron Content）TEC=∫sNedsTEC = \int\limits_s {Neds} TEC=s∫Neds单位为1TECU=10¹⁶e/m².
因此有−(Vion)G=+(Vion)P=40.3Nef2TEC- {({V_{ion}})_G} = + {({V_{ion}})_P} = 40.3{{Ne} \over {{f^2}}}TEC−(Vion)G=+(Vion)P=40.3f2NeTEC

双频改正模型

电离层延迟和频率平方成反比。卫星发射两种频率，可以认为两个频率的电磁波传播的是同一个路径，如果精确确定这两个频率信号到达接收机的时间差，那么就可以反推出各自受到的电离层延迟。
GPS卫星采用的L1和L2载波频率为：
f1=1575.42MHz{f_1} = 1575.42MHzf1=1575.42MHz, f2=1227.60MHz{f_2} = 1227.60MHzf2=1227.60MHz
令A=−40.3∫sNedsA = - 40.3\int\limits_s {Neds} A=−40.3s∫Neds于是有ρ=ρ′1+Af12\rho = {{\rho '}_1} + {A \over {f_1^2}}ρ=ρ′1+f12Aρ=ρ′2+Af22\rho = {{\rho '}_2} + {A \over {f_2^2}}ρ=ρ′2+f22A两式相减有ρ=ρ1′−ρ2′=c⋅t=Af22−Af12=(Vion)1[f12f22−1]\rho = {\rho '_1} - {\rho '_2} = c \cdot t = {A \over {f_2^2}} - {A \over {f_1^2}} = {({V_{ion}})_1}\left[ {{{f_1^2} \over {f_2^2}} - 1} \right]ρ=ρ1′−ρ2′=c⋅t=f22A−f12A=(Vion)1[f22f12−1]由此就可以解出两个频率的电离层延迟V_ion
同样的，可以得到无电离层组合Ion-Free组合ρ=f12f12−f22ρ1′+−f22f12−f22ρ2′=2.54573ρ1′−1.54573ρ2′\rho = {{f_1^2} \over {f_1^2 - f_2^2}}{\rho '_1} + {{ - f_2^2} \over {f_1^2 - f_2^2}}{\rho '_2} = 2.54573{\rho '_1} - 1.54573{\rho '_2}ρ=f12−f22f12ρ1′+f12−f22−f22ρ2′=2.54573ρ1′−1.54573ρ2′载波相位的推导和这个类似。

双频观测值建立VTEC模型

双频用户可以利用双频来消去电离层，而对于单频接收机，使用机构提供的VTEC模型是一个比较好的办法。
原理
与上述公式类似，只不过是解出A的值，因此有TEC=9.52437(ρ1′−ρ2′)(1TECU=1016e/m2)TEC = 9.52437({\rho '_1} - {\rho '_2}){\rm{ }}(1TECU = {10^{16}}e/{m^2})TEC=9.52437(ρ1′−ρ2′)(1TECU=1016e/m2)这里的TEC是卫星到接收机的路径上的电离层，也称为斜延迟，利用投影函数可以把斜延迟转化为天顶方向的垂直电离层延迟VTEC=TEC⋅M=f(B,L,t)VTEC = TEC \cdot M = f(B,L,t)VTEC=TEC⋅M=f(B,L,t)其中 B为穿刺点（将电离层市委薄层模型，卫星的信号穿过电离层模型的点）的大地纬度、L为大地经度，t为时间。如果某时段有w个观测历元，每个历元从y个监测站对z个卫星进行双频观测，那么就可以得到w×y×z个穿刺点的VTEC值。如考虑这个时间段电离层不变化，选择模型即可拟合出一个VTEC模型f(B,L,t)f(B,L,t)f(B,L,t)。
全球模型

IGS提供的VTEC格网图
从1998年开始，提供时段长度为2h，经差5°，纬差2.5°的VTEC格网图，标准IONEX数据格式。任一时间和位置通过格网内插获得。
CODE球谐函数模型
采用15阶15次的球谐函数建立全球的VTEC模型，提供球谐系数文件，使用者只要带入经纬度即可算出VTEC值。VTEC=∑n=0nmax⁡∑m=0nPˉnm(sin⁡φ)(Cˉnmcos⁡ms+Sˉnmsin⁡ms)VTEC = \sum\limits_{n = 0}^{{n_{\max }}} {\sum\limits_{m = 0}^n {{{\bar P}_{nm}}(\sin \varphi )} } ({\bar C_{nm}}\cos ms + {\bar S_{nm}}\sin ms)VTEC=n=0∑nmaxm=0∑nPˉnm(sinφ)(Cˉnmcosms+Sˉnmsinms)式中， n, m 分别表示对应球谐函数的阶数和次数，max表示球谐展开的最大阶数，Cˉnm{\bar C_{nm}}Cˉnm和Sˉnm{\bar S_{nm}}Sˉnm表示待估的球谐函数系数值。在全球建模时，N通常取15，即每次构建全球电离层模型有256(16*16)个球谐系数。得到左式的是IPP处的VTEC(vertical TEC)。

区域性模型
较多采用曲面拟合模型，VTEC看成纬差和太阳角时差的函数，表达式为：VTEC=∑i=0n∑j=0mEij(φ−φ0)i(S−S0)jVTEC = \sum\limits_{i = 0}^n {\sum\limits_{j = 0}^m {{E_{ij}}} } {(\varphi - {\varphi _0})^i}{(S - {S_0})^j}VTEC=i=0∑nj=0∑mEij(φ−φ0)i(S−S0)j