GPS从入门到放弃（十四） --- 电离层延时

GPS从入门到放弃（十四） — 电离层延时

电离层概念

电离层（Ionosphere）是地球大气的一个电离区域。它是受到太阳高能辐射以及宇宙线的激励而电离的大气高层。50千米以上的整个地球大气层都处于部分电离或完全电离的状态，电离层是部分电离的大气区域，完全电离的大气区域称磁层。

电离层的范围从离地面约50公里开始一直伸展到约1000公里高度的地球高层大气空域。
电离层的主要特性由电子密度、电子温度、碰撞频率、离子密度、离子温度和离子成分等空间分布的基本参数来表示。
电离层的研究对象主要是电子密度随高度的分布。电子密度（或称电子浓度）是指单位体积的自由电子数，随高度的变化与各高度上大气成分、大气密度以及太阳辐射通量等因素有关。电离层内任一点上的电子密度，决定于上述自由电子的产生、消失和迁移三种效应。在不同区域，三者的相对作用和各自的具体作用方式也大有差异。
电离层中存在相当多的自由电子和离子，能使无线电波改变传播速度，发生折射、反射和散射，产生极化面的旋转并受到不同程度的吸收。

电离层对卫星导航的影响

因为卫星都处于大气层外，卫星信号传播到地球上必须要经过大气层中的电离层，于是必然要受到电离层的影响。
电离层是弥散性介质，即介质的介电常数与频率有关，所以在电离层中，不同频率的电磁波具有不同的传播速度。
卫星信号载波的频率是固定的，称之为相速度，vp=λfv_p=\lambda fvp=λf，其中 λ\lambdaλ 是载波波长，fff 是载波频率。对应的相折射率为 npn_pnp。而载波上调制的数据信息可以看做是多个不同频率的波叠加而成的群波，其速度称之为群速度 vgv_gvg。对应的群折射率为 ngn_gng。
由电磁波的理论，我们可以得到群折射率与相折射率之间的关系为：
ng=np+fdnpdfn_g = n_p + f\frac{dn_p}{df} ng=np+fdfdnp

在弥散性介质中，相折射率 npn_pnp 和频率 fff 的关系可以表达为：
np=1+a1f2+a2f3+a3f4+⋯n_p = 1 + \frac{a_1}{f^2}+ \frac{a_2}{f^3} + \frac{a_3}{f^4} + \cdots np=1+f2a1+f3a2+f4a3+⋯

其中的系数 a1,a2,a3,⋯a_1, a_2, a_3, \cdotsa1,a2,a3,⋯ 与自由电子浓度有关，通过实际测量可以得到其数值。
忽略高阶项，有
np=1+a1f2n_p = 1 + \frac{a_1}{f^2} np=1+f2a1

于是可得
ng=1−a1f2n_g = 1 - \frac{a_1}{f^2} ng=1−f2a1

如果信号传播路径上的自由电子浓度用 NeN_eNe 表示，则系数 a1a_1a1 可以表示为：
a1=−40.28Nea_1 = -40.28 N_ea1=−40.28Ne

令路径自由电子总计数 TEC 为
TEC=∫SNe(s)ds\textrm{TEC} = \int_{S}N_e(s)ds TEC=∫SNe(s)ds

其中 SSS 为信号在电离层中的传播路径。则卫星信号在电离层中传播的相延迟量
Ip=∫S(np−1)ds=∫Sa1f2ds=−40.28f2∫SNe(s)ds=−40.28TECf2I_p= \int_{S}(n_p-1)ds = \int_{S}\frac{a_1}{f^2}ds = -\frac{40.28}{f^2}\int_{S}N_e(s)ds = -\frac{40.28\ \textrm{TEC}}{f^2}Ip=∫S(np−1)ds=∫Sf2a1ds=−f240.28∫SNe(s)ds=−f240.28 TEC

类似的，卫星信号在电离层中传播的群延迟量
Ig=∫S(ng−1)ds=−∫Sa1f2ds=40.28f2∫SNe(s)ds=40.28TECf2I_g= \int_{S}(n_g-1)ds = -\int_{S}\frac{a_1}{f^2}ds = \frac{40.28}{f^2}\int_{S}N_e(s)ds = \frac{40.28\ \textrm{TEC}}{f^2}Ig=∫S(ng−1)ds=−∫Sf2a1ds=f240.28∫SNe(s)ds=f240.28 TEC

对于卫星信号来说，对载波相位观测量造成影响的相延迟量，而对伪距观测量造成影响的是群延迟量。

双频接收机消除电离层延时

由前面的公式可知，电离层延时与载波频率有关，所以若有两个不同频率的观测量，我们就可以进行组合消除电离层延时。
对于一个双频接收机来说，设在两个不同频点 f1f_1f1、f2f_2f2 下的伪距观测量分别为 ρ1\rho_1ρ1、ρ2\rho_2ρ2，电离层延迟分别为 I1I_1I1、I2I_2I2，则有
I2−I1=ρ2−ρ1I_2 - I_1 = \rho_2 - \rho_1I2−I1=ρ2−ρ1

I1=f22f12−f22(ρ2−ρ1)I_1 = \frac{f_2^2}{f_1^2-f_2^2}(\rho_2-\rho_1) I1=f12−f22f22(ρ2−ρ1)

I2=f12f12−f22(ρ2−ρ1)I_2 = \frac{f_1^2}{f_1^2-f_2^2}(\rho_2-\rho_1) I2=f12−f22f12(ρ2−ρ1)

所以经过电离层校正后的伪距为
ρ1,2=ρ1−I1=ρ2−I2=f12f12−f22ρ1−f22f12−f22ρ2\rho_{1,2} = \rho_1-I_1 = \rho_2-I_2 = \frac{f_1^2}{f_1^2-f_2^2}\rho_1-\frac{f_2^2}{f_1^2-f_2^2}\rho_2ρ1,2=ρ1−I1=ρ2−I2=f12−f22f12ρ1−f12−f22f22ρ2

单频接收机消除电离层延时

单频接收机因为只有一个接收频点，无法通过双频来消除电离层延时，替代的方案是用电离层数学模型。
最常用的模型，也是 GPS 系统中用的模型，是 Klobuchar 模型。此模型包含8个电离层参数：α0,α1,α2,α3,β0,β1,β2,β3\alpha_0,\ \alpha_1,\ \alpha_2,\ \alpha_3,\ \beta_0, \ \beta_1, \ \beta_2, \ \beta_3α0, α1, α2, α3, β0, β1, β2, β3。在GPS导航电文中，第4子帧第18页广播了这8个参数。

根据这8个参数来计算电离层延时的方法如下：设 EEE、AAA 分别为接收机到卫星的仰角和方位角，ϕu\phi_uϕu、λu\lambda_uλu 分别为接收机的位置在WGS-84坐标系中的纬度和经度，TgpsT_{gps}Tgps为当前GPS系统时间，III 为最终计算出的电离层延时，则
ψ=0.0137E+0.11−0.022(semi-circle)\psi=\frac{0.0137}{E+0.11}-0.022 \qquad (\textrm{semi-circle})ψ=E+0.110.0137−0.022(semi-circle)

ϕi={ϕu+ψcos⁡A,∣ϕi∣≤0.416ifϕi>0.416,thenϕi=0.416ifϕi<−0.416,thenϕi=−0.416}(semi-circle)\phi_i = \left\{ \begin{array}{l} \phi_u + \psi \cos A, \qquad |\phi_i| \le 0.416 \\ \textrm{if} \: \phi_i>0.416, \quad \textrm{then} \: \phi_i = 0.416 \\ \textrm{if} \: \phi_i<-0.416, \: \textrm{then} \: \phi_i = -0.416 \end{array} \right\} \qquad (\textrm{semi-circle})ϕi=⎩⎨⎧ϕu+ψcosA,∣ϕi∣≤0.416ifϕi>0.416,thenϕi=0.416ifϕi<−0.416,thenϕi=−0.416⎭⎬⎫(semi-circle)

λi=λu+ψsin⁡Acos⁡ϕi(semi-circle)\lambda_i=\lambda_u+\frac{\psi\sin A}{\cos\phi_i} \qquad (\textrm{semi-circle})λi=λu+cosϕiψsinA(semi-circle)

ϕm=ϕi+0.064cos⁡(λi−1.617)(semi-circle)\phi_m=\phi_i+0.064\cos(\lambda_i-1.617) \qquad (\textrm{semi-circle})ϕm=ϕi+0.064cos(λi−1.617)(semi-circle)

t=43200λi+Tgps,0≤t<86400(sec)t=43200\lambda_i + T_{gps}, \quad 0\le t < 86400 \qquad (\textrm{sec})t=43200λi+Tgps,0≤t<86400(sec)

F=1.0+16.0×(0.53−E)3(semi-circle)F = 1.0 + 16.0 \times (0.53 - E)^3 \qquad (\textrm{semi-circle})F=1.0+16.0×(0.53−E)3(semi-circle)

AMP={∑n=03αnϕmn,AMP≥0if AMP<0,AMP=0}(sec)\textrm{AMP} = \left\{ \begin{array}{l} \displaystyle \sum_{n=0}^3\alpha^n \phi_m^n, \qquad \textrm{AMP} \ge 0 \\ \textrm{if AMP}<0, \qquad \textrm{AMP} = 0 \end{array} \right\} \qquad (\textrm{sec})AMP=⎩⎪⎨⎪⎧n=0∑3αnϕmn,AMP≥0if AMP<0,AMP=0⎭⎪⎬⎪⎫(sec)

PER={∑n=03βnϕmn,PER≥72000if PER<72000,PER=72000}(sec)\textrm{PER} = \left\{ \begin{array}{l} \displaystyle \sum_{n=0}^3\beta^n \phi_m^n, \qquad\qquad \textrm{PER} \ge 72000 \\ \textrm{if PER}<72000, \qquad \textrm{PER} = 72000 \end{array} \right\} \qquad (\textrm{sec})PER=⎩⎪⎨⎪⎧n=0∑3βnϕmn,PER≥72000if PER<72000,PER=72000⎭⎪⎬⎪⎫(sec)

x=2π(t−50400)PER(rad)x = \frac{2\pi(t-50400)}{\textrm{PER}} \qquad (\textrm{rad})x=PER2π(t−50400)(rad)

I={F×[5×10−9+A(1−x22+x424)],∣x∣<π2F×5×10−9,∣x∣≥π2}(sec)I = \left\{ \begin{array}{l} \displaystyle F\times\left[5\times 10^{-9}+A(1-\frac{x^2}{2}+\frac{x^4}{24}) \right], \qquad |x|<\frac{\pi}{2} \\ \displaystyle F\times 5\times 10^{-9}, \qquad\qquad\qquad\qquad\qquad\quad |x|\ge\frac{\pi}{2} \end{array} \right\} \qquad (\textrm{sec})I=⎩⎪⎨⎪⎧F×[5×10−9+A(1−2x2+24x4)],∣x∣<2πF×5×10−9,∣x∣≥2π⎭⎪⎬⎪⎫(sec)

以上每条公式中最后的括号内的表示运算单位，semi-circle是半圆，其与弧度的换算关系是 π\piπ 弧度等于一个半圆。

下面摘录了 RTKLIB 中用 Klobuchar 模型计算电离层延时的函数代码，可以跟以上公式进行一一对应。需要注意的是计算时的单位换算，azel[0], azel1分别是方位角和仰角，pos[0], pos1分别是接收机纬度和经度，它们的单位都是弧度，所以换算成半圆时需要除以π\piπ。

double ionmodel(gtime_t t, const double *ion, const double *pos,const double *azel)
{const double ion_default[]={ /* 2004/1/1 */0.1118E-07,-0.7451E-08,-0.5961E-07, 0.1192E-06,0.1167E+06,-0.2294E+06,-0.1311E+06, 0.1049E+07};double tt,f,psi,phi,lam,amp,per,x;int week;if (pos[2]<-1E3||azel[1]<=0) return 0.0;if (norm(ion,8)<=0.0) ion=ion_default; //若没有电离层参数，用默认参数/* earth centered angle (semi-circle) */psi=0.0137/(azel[1]/PI+0.11)-0.022;/* subionospheric latitude/longitude (semi-circle) */phi=pos[0]/PI+psi*cos(azel[0]);if      (phi> 0.416) phi= 0.416;else if (phi<-0.416) phi=-0.416;lam=pos[1]/PI+psi*sin(azel[0])/cos(phi*PI);/* geomagnetic latitude (semi-circle) */phi+=0.064*cos((lam-1.617)*PI);/* local time (s) */tt=43200.0*lam+time2gpst(t,&week);tt-=floor(tt/86400.0)*86400.0; /* 0<=tt<86400 *//* slant factor */f=1.0+16.0*pow(0.53-azel[1]/PI,3.0);/* ionospheric delay */amp=ion[0]+phi*(ion[1]+phi*(ion[2]+phi*ion[3]));per=ion[4]+phi*(ion[5]+phi*(ion[6]+phi*ion[7]));amp=amp<    0.0?    0.0:amp;per=per<72000.0?72000.0:per;x=2.0*PI*(tt-50400.0)/per;return CLIGHT*f*(fabs(x)<1.57?5E-9+amp*(1.0+x*x*(-0.5+x*x/24.0)):5E-9);
}

单频接收机用 Klobuchar 模型计算电离层延时大致能消除电离层延时的50%，而双频的方法能消除绝大部分，精度达到1m以下。