GAMP源码阅读：PPP中的模型改正：对流层延迟、电离层延迟

原始 Markdown文档、Visio流程图、XMind思维导图见：https://github.com/LiZhengXiao99/Navigation-Learning

在这里插入图片描述

文章目录

- 一、对流层延迟改正
- - 1、原理
  - 2、model_trop()：对流层改正入口函数
  - 3、tropmodel()：Saastamoinen 模型改正计算延迟
  - - 1. Saastamoinen 模型
    - 2. 标准大气模型
    - 3. GPT 模型
    - 4. tropmodel() 代码
  - 4、trop_model_prec()：精密对流层模型
  - 5、tropmapf()：计算干湿延迟投影系数
  - 6、tropmapf_nmf()：Neill 投影函数
- 二、电离层延迟改正
- - 1、原理
  - 2、model_iono()：电离层改正入口函数
  - 3、ionmodel()：克罗布歇模型计算电离层改正量
  - 4、iontec()：TEC 模型计算电离层改正量
  - - 1. readtec()：电离层 IONEX 文件读取
    - 2. iontec()：TEC 格网改正入口函数
    - 3. iondelay()：计算指定时间电离层延迟
    - 4. ionppp()：计算电离层穿刺点
    - 5. dataindex()：获取 TEC 格网数据下标
    - 6. interptec()：插值计算穿刺点处 TEC

一、对流层延迟改正

1、原理

对流层一般指距离地面 50km 内的大气层，是大气层质量的主要部分。当导航信号穿过对流层时，由于传播介质密度的增加，信号传播路径和传播速度会发生改变，由此引起的 GNSS 观测值误差称为对流层延迟。对流层延迟一般可分为干延迟和湿延迟，对于载波相位和伪距完全相同，一般在米级大小，可通过模型改正和参数估计的方法来削弱其影响。修正模型如下：
$T=M_{d r y} T_{d r y}+M_{w e t} T_{w e t}$
式中， $T_{d r y}, T_{w e t}$ 分别表示接收机天顶对流层的干延迟和湿延迟； $M_{d r y}, M_{w e t}$ 分别表示干延迟和湿延迟的投影函数。对流层干延迟比较稳定，主要与测站高度、大气温度和大气压相关，可通过模型改正，常用模型有 Saastamoninen 模型、Hopfield 模型等。湿延迟不同于干延迟，变化较大，主要与水汽含量相关，一般估计天顶对流层湿延迟，通过投影函数计算各卫星的电离层湿延迟，常用的投影函数有全球投影函数（Global Mapping Function，GMF）、Niell 投影函数（NMF）和 Vienna投影函数（Vienna Mapping Function，VMF）等。本文使用 Saastamoninen 模型改正对流层干延迟，使用 GMF 投影函数估计接收机天顶对流层湿延迟。

2、model_trop()：对流层改正入口函数

Saastamoinen 模型：调用 tropmodel() 计算对流层延迟，方差设为 0.3*0.3。
估计对流层：从 x 中取估计的 ZTD，调用 trop_model_prec() 改正，方差在其中计算(0.01*0.01)。

计算得到的对流层改正量，作为参数 dtrp 传出

static int model_trop(gtime_t time, const double *pos, const double *azel,
	const prcopt_t *opt, const double *x, double *dtdx,
	const nav_t *nav, double *dtrp, double *shd, double *var)
{
	double trp[3]={0};
	double zwd;

	// Saastamoinen 模型改正计算延迟，方差设为 0.3*0.3
	if (opt->tropopt==TROPOPT_SAAS) {
		*dtrp=tropmodel(time,pos,azel,REL_HUMI,&zwd,0);
		*dtrp+=zwd;
		*var=SQR(ERR_SAAS);
		return 1;
	}
	// 估计对流层模式，从 x 中取估计的 ZTD，调用 trop_model_prec() 改正，方差在其中计算
	if (opt->tropopt==TROPOPT_EST||opt->tropopt==TROPOPT_ESTG) {
		matcpy(trp,x+IT(opt),opt->tropopt==TROPOPT_EST?1:3,1);
		*dtrp=trop_model_prec(time,pos,azel,opt,trp,dtdx,shd,var);
		return 1;
	}
	return 0;
}

3、tropmodel()：Saastamoinen 模型改正计算延迟

1. Saastamoinen 模型

Saastamoinen 模型将对流层分为两层进行积分，第 1 层为从地表到 $\mathrm{~km}$ 高度的对流层顶，该层的温度变化率为；第 2 层为从 $\mathrm{~km}$ 到 $\mathrm{~km}$ 高度的平流层顶，该层的温度视为常数。 Saastamoinen 模型首次将被积函数按照天顶距三角函数展开逐项进行积分, 并把对流层天顶延迟分为对流层干延迟和湿延迟两个分量之和，两个分量的表达式为：
$\left\{\begin{array}{l}Z H D=\frac{0.002277 P}{f(\varphi, h)} \\ Z W D=\frac{0.002277 e}{f(\varphi, h)}\left(0.05+\frac{1255}{T}\right) \\ f(\varphi, h)=1-0.00266 \cos (2 \varphi)-0.00028 h\end{array}\right.$
或者：
$T_{r}^{s}=\frac{0.002277}{\cos z}\left\{p+\left(\frac{1255}{T}+0.05\right) e-\tan ^{2} z\right\}$
其中， $\varphi$ 和 $h$ 分别为地表气压 $(\mathrm{hPa})$ 、地表温度 $(\mathrm{K})$ 、水汽压 $(\mathrm{hPa})$ 、测站纬度 $(\mathrm{rad})$ 和高程 $(\mathrm{km})$ ，可以由标准大气模型基于经验计算，也可由 GPT 模型提供。

2. 标准大气模型

根据经验模型计算求大气压 P、温度 T、大气水汽压力 e：
$\begin{array}{l}p=1013.25 \times\left(1-2.2557 \times 10^{-5} h\right)^{5.2568} \\ T=15.0-6.5 \times 10^{-3} h+273.15 \\ e=6.108 \times \exp \left\{\frac{17.15 T-4684.0}{T-38.45}\right\} \times \frac{h_{r e l}}{100}\end{array}$

其中 $z$ 是天顶角，与高度角互余： $z=\pi / 2-E l_{r}^{s}$

3. GPT 模型

GPT 系列模型是 Boehm 等利用欧洲中尺度天气预报中心 (European Centre for Medium-Range Weather Forecasts, ECMWF) 长期的再分析气象资料建立的全球气象参数经验模型, 仅需知道测站地理位置信息与年积日便可以获得地表温度、大气压力和水汽压等气象参数，在全球范围内得到广泛应用。

作者用的 ftp://tai.bipm.org/iers/convupdt/chapter9/GPT.F 文件对应的模型，把 IGS 中心提供的函数改成 c 语言版本的 getGPT() 函数，这个 GPT.F 可以通过资源管理器下载，用 notepade++、或者 VScode 查看：

5f5cec497ffb8ae62cb16a7ae2519cff

4. tropmodel() 代码

用标准大气模型或 GPT 模型计算一些标准大气值，包括总压 p、大气温度 T、水汽压力 e；然后计算了静力学延迟(干延迟) trph计算了湿延迟 trpw。

humi 传 0，就只计算干延迟，湿延迟作为参数估计

extern double tropmodel(gtime_t time, const double *pos, const double *azel,
                        double humi, double *zwd, int atmodel)
{
    const double temp0=15.0;   /* temparature at sea level */
    double hgt,pres,temp,e,z,trph,trpw;
	int b1,b2;
	mjd_t mjd;
	double dmjd,undo,d1,d2;
    
    if (pos[2]<-100.0||pos[2]>1E6||azel[1]<=0) {
        //���㿪ʼ�ĵ�һ����Ԫ�߳�������-100.0����������������Ϣ
        b1=(PPP_Glo.iEpoch==1)&&(PPP_Glo.revs==0);
        b2=(PPP_Glo.nEpoch-1==PPP_Glo.iEpoch)&&(PPP_Glo.revs==1);

        //����deltaEp����
        if (!b1&&!b2&&PPP_Glo.delEp<20) {
            sprintf(PPP_Glo.chMsg,"*** WARNING: tropmodel: height=%7.3f elev=%4.1f\n",pos[2],azel[1]*R2D);
            outDebug(OUTWIN,OUTFIL,OUTTIM);
        }

        return 0.0;
    }
    if (pos[2]>=1.0/2.2557E-5) return 0.0;

    hgt=pos[2]<0.0?0.0:pos[2];

    //�̹߳���ᵼ�¼��������ֵ���������ӳ���Ϊ�����
    if (hgt>15000.0) {
        //���㿪ʼ�ĵ�һ����Ԫ�߳�������-100.0����������������Ϣ
        b1=(PPP_Glo.iEpoch)==1&&(PPP_Glo.revs==0)&&(PPP_Glo.prcOpt_Ex.solType==0||
			PPP_Glo.prcOpt_Ex.solType==3||PPP_Glo.prcOpt_Ex.solType==4);
        b2=(PPP_Glo.nEpoch-1==PPP_Glo.iEpoch)&&(PPP_Glo.revs==1)&&(PPP_Glo.prcOpt_Ex.solType==1
			||PPP_Glo.prcOpt_Ex.solType==2);

        if (!b1&&!b2&&PPP_Glo.delEp<20) {
			sprintf(PPP_Glo.chMsg,"*** WARNING: tropmodel: height=%7.3f\n",hgt);
            outDebug(OUTWIN,OUTFIL,OUTTIM);
        }
        hgt=15000.0;
    }

    // 计算大气参数，标准大气模型 or GPT 模型
    /* standard atmosphere */
    if (atmodel!=1) {
        pres=1013.25*pow(1.0-2.2557E-5*hgt,5.2568);         // 求大气压P (E.5.1)
        temp=temp0-6.5E-3*hgt+273.16;                       // 求温度temp (E.5.2)
        e=6.108*humi*exp((17.15*temp-4684.0)/(temp-38.45)); // 求大气水汽压力e (E.5.3)
    }
    else {
        time2mjd(time,&mjd);
        dmjd=mjd.day+(mjd.ds.sn+mjd.ds.tos)/86400.0;
        
        undo=0.0;
        getGPT(pos,dmjd,&pres,&temp,&undo);

        d1=1013.25*pow(1.0-2.2557E-5*hgt,5.2568);
        d2=15.0-6.5E-3*hgt;

        temp+=273.16;
        e=6.108*humi*exp((17.15*temp-4684.0)/(temp-38.45));
    }
    
    // 计算对流层延迟
    /* saastamoninen model */
    z=PI/2.0-azel[1];                   // 求天顶角z 卫星高度角azel[1]的余角
    trph=0.0022768*pres/(1.0-0.00266*cos(2.0*pos[0])-0.00028*hgt/1E3)/cos(z);   // 干延迟
    trpw=0.002277*(1255.0/temp+0.05)*e/cos(z);  // 湿延迟

    *zwd=trpw;

    return trph;
}

4、trop_model_prec()：精密对流层模型

先调用 tropmodel() 相对湿度传 0 计算干延迟 zhd。
调用 tropmapf() NMF 投影函数计算干延迟投影系数 m_h、湿延迟投影系数 m_w。
格网梯度模型计算梯度系数
$\begin{array}{l}m\left(E l_{r}^{s}\right)=m_{W}\left(E l_{r}^{s}\right)\left\{1+\cot E l_{r}^{s}\left(G_{N, r} \cos A z_{r}^{s}+G_{E, r} \sin A z_{r}^{s}\right)\right\} \\ T_{r}^{s}=m_{H}\left(E l_{r}^{s}\right) Z_{H, r}+m\left(E l_{r}^{s}\right)\left(Z_{T, r}-Z_{H, r}\right)\end{array}$
方差设为 0.01*0.01。
根据干湿延迟和投影系数计算对流层延迟改正量，返回。
$T=M_{d r y} T_{d r y}+M_{w e t} T_{w e t}$

static double trop_model_prec(gtime_t time, const double *pos, const double *azel, 
	const prcopt_t *opt, const double *x, double *dtdx, double *shd, double *var)
{
	const double zazel[]={0.0,PI/2.0};
	double zhd,zwd,m_h,m_w,cotz,grad_n,grad_e;

	// 利用对流层误差模型进行改正，然后将残余误差当作一个未知参数进行估计
	/* zenith hydrostatic delay */
	zhd=tropmodel(time,pos,zazel,0.0,&zwd,1);
    PPP_Glo.zhd=zhd;

	// 调用 tropmapf() NMF 投影函数计算干延迟投影系数 m_h、湿延迟投影系数 m_w
	/* mapping function */
	m_h=tropmapf(time,pos,azel,&m_w);

	*shd=m_h*zhd;

	// 格网梯度模型计算梯度系数
	if (opt->tropopt>=TROPOPT_ESTG&&azel[1]>0.0) {
		/* m_w=m_0+m_0*cot(el)*(Gn*cos(az)+Ge*sin(az)): ref [6]  (E.5.5) */
		cotz=1.0/tan(azel[1]);
		grad_n=m_w*cotz*cos(azel[0]);
		grad_e=m_w*cotz*sin(azel[0]);
		m_w+=grad_n*x[1]+grad_e*x[2];
		dtdx[1]=grad_n*(x[0]);			// (E.5.6)
		dtdx[2]=grad_e*(x[0]);
	}
	dtdx[0]=m_w;
	*var=SQR(0.01);	// 方差设为 0.01*0.01

	//return m_h*zhd+m_w*(x[0]-zhd);
	return m_h*zhd+m_w*(x[0]);
}

5、tropmapf()：计算干湿延迟投影系数

干投影函数是通过返回值获得的，而湿投影是通过输入/输出参数mapfw获得的，有两种投影函数的计算方法，分别是 GMF 和 NMF ，默认使用的是 NMF 方法，RTKLIB 可以通过定义IERS_MODEL宏来使用 GMF 方法；GAMP 作者好像实现了一个 tropmapf_gmf() 函数，但是没有调用。

extern double tropmapf(gtime_t time, const double pos[], const double azel[],
                       double *mapfw)
{
#ifdef IERS_MODEL
    const double ep[]={2000,1,1,12,0,0};
    double mjd,lat,lon,hgt,zd,gmfh,gmfw;
#endif
    
    if (pos[2]<-1000.0||pos[2]>20000.0) {
        if (mapfw) *mapfw=0.0;
        return 0.0;
    }
#ifdef IERS_MODEL
    mjd=51544.5+(timediff(time,epoch2time(ep)))/86400.0;
    lat=pos[0];
    lon=pos[1];
    hgt=pos[2]-geoidh(pos); /* height in m (mean sea level) */
    zd =PI/2.0-azel[1];
    
    /* call GMF */
    gmf_(&mjd,&lat,&lon,&hgt,&zd,&gmfh,&gmfw);
    
    if (mapfw) *mapfw=gmfw;
    return gmfh;
#else
    return tropmapf_nmf(time,pos,azel,mapfw); /* NMF */
#endif
}

6、tropmapf_nmf()：Neill 投影函数

$\begin{array}{c}M_{d y}=\frac{1+\frac{a_{d r v}}{1+\frac{b_{d r y}}{1+c_{d r y}}}}{\sin E+\frac{a_{d r y}}{\sin E+\frac{b_{d r y}}{\sin E+c_{d r y}}}}+\left(\frac{1}{\sin E}-\frac{1+\frac{a_{h t}}{1+\frac{b_{h t}}{1+c_{h t}}}}{\sin E+\frac{a_{h t}}{\sin E+\frac{b_{h t}}{\sin E+c_{h t}}}}\right) \\ M_{\text {wet }}=\frac{1+\frac{a_{\text {wet }}}{1+\frac{b_{w e t}}{1+c_{\text {wet }}}}}{\sin E+\frac{a_{w e t}}{\sin E+\frac{b_{\text {wet }}}{\sin E+c_{w e t}}}} \\ a_{d r y}=a_{c v g}(\varphi)+a_{u m p}(\varphi) \cos \left(2 \pi \frac{D O Y-28}{365.25}\right)\end{array}$

其中， $E$ 为卫星高度角， $h$ 为测站高程，系数 $a_{h t}=2.53 e-5, b_{h t}=5.49 e-3$ ， $c_{h t}=1.14 e-3, a_{a v g}(\varphi), a_{a m p}(\varphi)$ 查表内插的方法获得， $D O Y$ 为年积日，如果测站在南半球时 $D O Y - 180$

extern double tropmapf_nmf(gtime_t time, const double pos[], const double azel[],
                  double *mapfw)
{
    /* ref [5] table 3 */
    /* hydro-ave-a,b,c, hydro-amp-a,b,c, wet-a,b,c at latitude 15,30,45,60,75 */
    const double coef[][5]={
        { 1.2769934E-3, 1.2683230E-3, 1.2465397E-3, 1.2196049E-3, 1.2045996E-3},
        { 2.9153695E-3, 2.9152299E-3, 2.9288445E-3, 2.9022565E-3, 2.9024912E-3},
        { 62.610505E-3, 62.837393E-3, 63.721774E-3, 63.824265E-3, 64.258455E-3},
        
        { 0.0000000E-0, 1.2709626E-5, 2.6523662E-5, 3.4000452E-5, 4.1202191E-5},
        { 0.0000000E-0, 2.1414979E-5, 3.0160779E-5, 7.2562722E-5, 11.723375E-5},
        { 0.0000000E-0, 9.0128400E-5, 4.3497037E-5, 84.795348E-5, 170.37206E-5},
        
        { 5.8021897E-4, 5.6794847E-4, 5.8118019E-4, 5.9727542E-4, 6.1641693E-4},
        { 1.4275268E-3, 1.5138625E-3, 1.4572752E-3, 1.5007428E-3, 1.7599082E-3},
        { 4.3472961E-2, 4.6729510E-2, 4.3908931E-2, 4.4626982E-2, 5.4736038E-2}
    };
    const double aht[]={ 2.53E-5, 5.49E-3, 1.14E-3}; /* height correction */
    
    double y,cosy,ah[3],aw[3],dm,el=azel[1],lat=pos[0]*R2D,hgt=pos[2];
    int i;
    
    if (el<=0.0) {
        if (mapfw) *mapfw=0.0;
        return 0.0;
    }
    /* year from doy 28, added half a year for southern latitudes */
    y=(time2doy(time)-28.0)/365.25+(lat<0.0?0.5:0.0);
    
    cosy=cos(2.0*PI*y);
    lat=fabs(lat);
    
    for (i=0;i<3;i++) {
        ah[i]=interpc(coef[i  ],lat)-interpc(coef[i+3],lat)*cosy;
        aw[i]=interpc(coef[i+6],lat);
    }
    /* ellipsoidal height is used instead of height above sea level */
    dm=(1.0/sin(el)-mapf(el,aht[0],aht[1],aht[2]))*hgt/1E3;
    
    if (mapfw) *mapfw=mapf(el,aw[0],aw[1],aw[2]);   // mapfw 湿延迟函数
    
    return mapf(el,ah[0],ah[1],ah[2])+dm;   // 返回干延迟投影函数
}

二、电离层延迟改正

1、原理

受太阳辐射的影响，距地面 60km 以上的大气层处于部分电离或完全电离状态，该区域被称为电离层。当电磁波信号通过电离层时，传播速度和传播路径会发生改变，给 GNSS 观测值带来误差，即电离层延迟。电离层延迟大小由电子密度和信号频率决定，影响可达数十米。由于 GNSS 信号的特性，同一频率载波相位和伪距观测值电离层延迟大小相等，符号相反。在一些条件下，可使用经验模型改正或约束观测值中的电离层延迟，常用的电离层延迟经验模型有 Klobuchar 模型、Bent 模型和电离层格网模型等。本文对电离层延迟的处理使用最常用的两种方法：使用非差非组合模型估计各卫星第一频率倾斜电离层延迟；使用无电离层组合消除电离层延迟一阶项，忽略二阶及以上项。

2、model_iono()：电离层改正入口函数

TEC 格网模型改正：调用 iontec() 计算改正量，方差也在 iontec() 中计算。
广播星历改正：调用 ionmodel() 克罗布歇模型计算改正量，方差设为 0.5 乘以电离层改正量再平方。
电离层估计：从参数向量 x 中取改正量，方差设为 0。
L1/L2消电离层组合：电离层改正量和方差设为 0。

计算得到的 L1 频率的电离层改正量，作为参数 dion 传出，当使用其它频率信号时，依据所用信号频组中第一个频率的波长与 L1 波长的比例关系，对上一步得到的电离层延时进行修正，不考虑模糊度情况下改正公式为：
$I=\frac{\Phi_{2}-\Phi_{1}}{1-\left(f_{1} / f_{2}\right)^{2}}$

3、ionmodel()：克罗布歇模型计算电离层改正量

使用克罗布歇模型计算 L1 的电离层改正量，将晚间的电离层时延视为常数，取值为 5ns，把白天的时延看成是余弦函数中正的部分。于是天顶方向调制在 L1 载波上的测距码的电离层时延可表示为：
$T_{g}=5 \times 10^{-9}+A \cos \frac{2 \pi}{P}\left(t-14^{h}\right)$
振幅 $A$ 和周期 $P$ 是模型中需要计算的部分，分别为：
$\begin{array}{l} A=\sum_{i=0}^{3} \alpha_{i}\left(\varphi_{m}\right)^{i} \\ P=\sum_{i=0}^{3} \beta_{i}\left(\varphi_{m}\right)^{i} \end{array}$
全球定位系统向单频接收机用户提供的电离层延迟改正时就采用上述模型。其中 $\alpha_i$ 和 $\beta_i$ 是地面控制系统根据该天为一年中的第几天（将一年分为 37 个区间）以及前 5 天太阳的平均辐射流量（共分10档）从 370 组常数中选取的，然后编入星导航电文播发给用户。卫星在其播发的导航电文中提供了这 8 个电离层延迟参数：
$p_{\text {ion }}=\left(\alpha_{0}, \alpha_{1}, \alpha_{2}, \alpha_{3}, \beta_{0}, \beta_{1}, \beta_{2}, \beta_{3}\right)^{T}$
根据根据参数 $\alpha_0，\alpha_1，\alpha_2，\alpha_3$ 确定振幅 $A$ ，根据根据参数 $\beta_0，\beta_1，\beta_2，\beta_3$ 确定周期 $T$ ，再给定一个以秒为单位的当地时间 $t$ ，就能算出天顶电离层延迟。再由天顶对流层延迟根据倾斜率转为卫星方向对流层延迟。

电离层分布在离地面 60-1000km 的区域内。当卫星不在测站的天顶时，信号传播路径上每点的地方时和纬度均不相同，为了简化计算，我们将整个电离层压缩为一个单层，将整个电离层中的自由电子都集中在该单层上，用它来代替整个电离层。这个电离层就称为中心电离层。中心电离层离地面的高度通常取 350km。需要计算卫星言号传播路径与中心电离层的交点 $P^{\prime}$ 的时角 $t$ 和地磁纬度 $\varphi_{m}$ ，因为只有 $P^{\prime}$ 才能反映卫星信号所受到的电离层延迟的总的情况。

综上，已知大地经度、大地纬度、卫星的高度角和卫星测站的方位角，电离层延迟计算方法如下：

计算测站 $P$ 和 $P^{\prime}$ 在地心的夹角：
$\psi=0.0137 /(E l+0.11)-0.022$
计算交点 $P^{\prime}$ 的地心纬度：
$\varphi_{i}=\varphi+\psi \cos A z$

$\left\{\begin{array}{ll}\varphi_{l}>+0.416 & \varphi_{l}=+0.416 \\ \varphi{l}>-0.416 & \varphi_{l}=-0.416\end{array}\right.$

计算交点 $P^{\prime}$ 的地心经度：
$\lambda_{i}=\lambda+\psi \sin A z / \cos \varphi_{i}$
计算地磁纬度：
$\varphi_{m}=\varphi_{i}+0.064 \cos \left(\lambda_{i}-1.617\right)$
计算观测瞬间交点 $P^{\prime}$ 处的地方时：
$\times 10^{4} \lambda_{i}+t$

$\left\{\begin{array}{ll}t>86400 & t=t-86400 \\ t<0 & t=t+86400\end{array}\right.$

计算倾斜因子：
$\times(0.53-E l)^{3}$
最后计算电离层时间延迟：
$\begin{array}{l}x=2 \pi(t-50400) / \sum_{n=0}^{3} \beta_{n} \varphi_{m}{ }^{n} \\ I_{r}^{s}=\left\{\begin{array}{cc}F \times 5 \times 10^{-9} \\ F \times\left(5 \times 10^{-9}+\sum_{n=1}^{4} \alpha_{n} \varphi_{m}{ }^{n} \times\left(1-\frac{x^{2}}{2}+\frac{x^{4}}{24}\right)\right) & (|x|>1.57)\end{array}\right.\end{array}$

extern 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;         // 计算地心角(E.5.6)
    
    /* subionospheric latitude/longitude (semi-circle) */   
    phi=pos[0]/PI+psi*cos(azel[0]);                 // 计算穿刺点地理纬度(E.5.7)
    if      (phi> 0.416) phi= 0.416;        // phi不超出(-0.416,0.416)范围
    else if (phi<-0.416) phi=-0.416;
    lam=pos[1]/PI+psi*sin(azel[0])/cos(phi*PI);     // 计算穿刺点地理经度(E.5.8)
    
    /* geomagnetic latitude (semi-circle) */
    phi+=0.064*cos((lam-1.617)*PI);                 // 计算穿刺点地磁纬度(E.5.9)
    
    /* local time (s) */
    tt=43200.0*lam+time2gpst(t,&week);              // 计算穿刺点地方时(E.5.10)
    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);            // 计算投影系数(E.5.11)
    
    /* 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;                      // (E.5.12)
    
    return CLIGHT*f*(fabs(x)<1.57?5E-9+amp*(1.0+x*x*(-0.5+x*x/24.0)):5E-9);     //(E.5.13)
}

4、iontec()：TEC 模型计算电离层改正量

格网改正模型电离层格网模型文件 IONEX，通过内插获得穿刺点位置，并结合当天电离层格网数据求出穿刺点的垂直电子含量，获得电离层延迟误差。

1. readtec()：电离层 IONEX 文件读取

IONEX 文件分四大部分：文件头结束于END OF HEADER。多组总电子含量，每组以START OF TEC MAP ，结束于END OF TEC MAP 。多组电子含量均方根误差 ，与总电子含量对应，开始于START OF RMS MAP，结束于END OF RMS MAP，DCB数据块开始于START OF AUS DATA，结束于END OF AUS DATA，也称辅助数据块（Auxiliary Data Blocks)。

readtec() 执行流程：

开辟内存空间
扩展file的*到efile[] ，遍历efile[] ：
fopen()以读的方式打开
调用readionexh()读 ionex 文件头
调用readionexb()读 ionex 文件体，其中会调用 addtec() 将 tec 格网数据存到nav->tec[]
读取完之后，调用combtec()合并 tec 格网数据
存 DCB 参数到nav->cbias

extern void readtec(const char *file, nav_t *nav, int opt)
{
    FILE *fp;
    double lats[3]={0},lons[3]={0},hgts[3]={0},rb=0.0,nexp=-1.0;
    double dcb[MAXSAT]={0},rms[MAXSAT]={0};
    int i,n;
    char *efiles[MAXEXFILE];
    
    trace(3,"readtec : file=%s\n",file);
    
    /* clear of tec grid data option */
    if (!opt) { //如果没有opt，释放nav->tec，nav->ntmax置0
        free(nav->tec); nav->tec=NULL; nav->nt=nav->ntmax=0;
    }
    for (i=0;i<MAXEXFILE;i++) {     //为efile[]开辟空间
        if (!(efiles[i]=(char *)malloc(1024))) {
            for (i--;i>=0;i--) free(efiles[i]);
            return;
        }
    }
    /* expand wild card in file path */
    n=expath(file,efiles,MAXEXFILE);    //扩展file的*到efile[]
    
    //遍历efile[]
    for (i=0;i<n;i++) {
        if (!(fp=fopen(efiles[i],"r"))) {       //以读的方式打开
            trace(2,"ionex file open error %s\n",efiles[i]);
            continue;
        }
        /* read ionex header */     //调用readionexh()读ionex文件头
        if (readionexh(fp,lats,lons,hgts,&rb,&nexp,dcb,rms)<=0.0) {     
            trace(2,"ionex file format error %s\n",efiles[i]);
            continue;
        }
        /* read ionex body */       //调用readionexb()读ionex文件体
        readionexb(fp,lats,lons,hgts,rb,nexp,nav);
        
        fclose(fp);
    }
    for (i=0;i<MAXEXFILE;i++) free(efiles[i]);
    
    /* combine tec grid data */
    if (nav->nt>0) combtec(nav);    //调用combtec()合并tec格网数据
    
    /* P1-P2 dcb */
    for (i=0;i<MAXSAT;i++) {
        nav->cbias[i][0]=CLIGHT*dcb[i]*1E-9; /* ns->m */	//存DCB到nav->cbias
    }
}

调用函数：

readionexh()：读取文件头，循环读取每一行，根据注释读取前面内容，如果遇到START OF AUX DATA ，调用 readionexdcb() 读取 DCB 参数。
readionexdcb()：循环读取 DCB 参数和对应的均方根误差，直到END OF AUS DATA
readionexb()：循环读取 TEC 格网数据和均方根误差，type为 1 是 TEC 格网，type为2是均分根误差，调用addtec()将tec格网数据存到nav->tec[]
combtec()：合并nav->tec[]中时间相同的项。

2. iontec()：TEC 格网改正入口函数

由所属时间段两端端点的 TEC 网格数据时间插值计算出电离层延迟 (L1) (m)

检测高度角和接收机高度是否大于阈值
从 nav_t.tec中找出第一个tec[i].time>time信号接收时间
确保 time是在所给出的nav_t.tec包含的时间段之中，通过确认所找到的时间段的右端点减去左端点，来确保时间间隔不为0
调用iondelay()来计算所属时间段两端端点的电离层延时
由两端的延时，插值计算出观测时间点处的值

extern int iontec(gtime_t time, const nav_t *nav, const double *pos,
                  const double *azel, int opt, double *delay, double *var)
{
    double dels[2],vars[2],a,tt;
    int i,stat[2];

    // 检测高度角和接收机高度是否大于阈值
    if (azel[1]<MIN_EL||pos[2]<MIN_HGT) {
        *delay=0.0;
        *var=VAR_NOTEC;
        return 1;
    }
    // 从 nav_t.tec 中找出第一个 tec[i].time>time 信号接收时间
    for (i=0;i<nav->nt;i++) {
        if (timediff(nav->tec[i].time,time)>0.0) break;
    }
    // 确保 time 是在所给出的 nav_t.tec 包含的时间段之中
    if (i==0||i>=nav->nt) {
        printf("%s: tec grid out of period\n",time_str(time,0));
        return 0;
    }
    // 通过确认所找到的时间段的右端点减去左端点，来确保时间间隔 != 0
    if ((tt=timediff(nav->tec[i].time,nav->tec[i-1].time))==0.0) {
        printf("tec grid time interval error\n");
        return 0;
    }
    // 调用 iondelay 来计算所属时间段两端端点的电离层延时
    /* ionospheric delay by tec grid data */
    stat[0]=iondelay(time,nav->tec+i-1,pos,azel,opt,dels  ,vars  );
    stat[1]=iondelay(time,nav->tec+i  ,pos,azel,opt,dels+1,vars+1);
    
    // 由两端的延时，插值计算出观测时间点处的值
    if (!stat[0]&&!stat[1]) {   // 两个端点都计算出错，输出错误信息，返回 0
        printf("%s: tec grid out of area pos=%6.2f %7.2f azel=%6.1f %5.1f\n",
              time_str(time,0),pos[0]*R2D,pos[1]*R2D,azel[0]*R2D,azel[1]*R2D);
        return 0;
    }
    // 两个端点都有值，线性插值出观测时间点的值，返回 1
    if (stat[0]&&stat[1]) { /* linear interpolation by time */
        a=timediff(time,nav->tec[i-1].time)/tt;
        *delay=dels[0]*(1.0-a)+dels[1]*a;
        *var  =vars[0]*(1.0-a)+vars[1]*a;
    }
    // 只有一个端点有值，将其结果作为观测时间处的值，返回 1
    else if (stat[0]) { /* nearest-neighbour extrapolation by time */
        *delay=dels[0];
        *var  =vars[0];
    }
    else {
        *delay=dels[1];
        *var  =vars[1];
    }

    return 1;
}

3. iondelay()：计算指定时间电离层延迟

while大循环tec->ndata[2]次：
调用ionppp()函数，计算当前电离层高度，穿刺点的位置 {lat,lon,h} (rad,m)和倾斜率
按照M-SLM映射函数重新计算倾斜率
在日固系中考虑地球自转，重新计算穿刺点经度
调用interptec()格网插值获取vtec
*delay+=fact*fs*vtec ,*var+=fact*fact*fs*fs*rms*rms

static int iondelay(gtime_t time, const tec_t *tec, const double *pos,
                    const double *azel, int opt, double *delay, double *var)
{
    const double fact=40.30E16/FREQ1/FREQ1; /* tecu->L1 iono (m) */
    double fs,posp[3]={0},vtec,rms,hion,rp;
    int i;
    
    *delay=*var=0.0;
    
    for (i=0;i<tec->ndata[2];i++) { /* for a layer */
        hion=tec->hgts[0]+tec->hgts[2]*i;

        // 调用ionppp()函数，计算当前电离层高度，穿刺点的位置 {lat,lon,h} (rad,m)和倾斜率( )
        /* ionospheric pierce point position */
        fs=ionppp(pos,azel,tec->rb,hion,posp);
        
        // 按照 M-SLM 映射函数重新计算倾斜率
        /*if (opt&2) {*/
        /* modified single layer mapping function (M-SLM) ref [2] */
        rp=tec->rb/(tec->rb+hion)*sin(0.9782*(PI/2.0-azel[1]));
        fs=1.0/sqrt(1.0-rp*rp);

        // 在日固系中考虑地球自转，重新计算穿刺点经度
        //} 
        //if (opt&1) {
        //    /* earth rotation correction (sun-fixed coordinate) */
        //    posp[1]+=2.0*PI*timediff(time,tec->time)/86400.0;
        //}

        // 调用 interptec() 格网插值获取 vtec
        /* interpolate tec grid data */
        if (!interptec(tec,i,posp,&vtec,&rms)) return 0;
        
        *delay+=fact*fs*vtec;
        *var+=fact*fact*fs*fs*rms*rms;
    }
    
    return 1;
}

4. ionppp()：计算电离层穿刺点

位置 {lat,lon,h} (rad,m)和倾斜率做返回值
$\begin{array}{l}z=\pi / 2-E l_{r}^{s} \\ z^{\prime}=\arcsin \left(\frac{R_{E}}{R_{E}+H} \sin z\right) \\ \alpha=z-z^{\prime} \\ \phi_{I P P}=\arcsin \left(\cos \alpha \sin \phi_{r}+\sin \alpha \cos \phi_{r} \cos A z_{r}^{s}\right)\end{array}$

extern double ionppp(const double *pos, const double *azel, double re,
                     double hion, double *posp)
{
    double cosaz,rp,ap,sinap,tanap;
    // z 并不是仰角 azel[1]，而是仰角关于的补角，所以程序中在计算 rp 是采用的是 cos(azel[1]) 的写法
    rp=re/(re+hion)*cos(azel[1]);
    ap=PI/2.0-azel[1]-asin(rp);     // (E.5.14) (E.5.16)
    sinap=sin(ap);
    tanap=tan(ap);
    cosaz=cos(azel[0]);
    posp[0]=asin(sin(pos[0])*cos(ap)+cos(pos[0])*sinap*cosaz);      // (E.5.17)
    
    if ((pos[0]> 70.0*D2R&& tanap*cosaz>tan(PI/2.0-pos[0]))||
        (pos[0]<-70.0*D2R&&-tanap*cosaz>tan(PI/2.0+pos[0]))) {
        posp[1]=pos[1]+PI-asin(sinap*sin(azel[0])/cos(posp[0]));    // (E.5.18a)
    }
    else {
        posp[1]=pos[1]+asin(sinap*sin(azel[0])/cos(posp[0]));       // (E.5.18b)
    }
    return 1.0/sqrt(1.0-rp*rp); // 返回倾斜率

    // 可能因为后面再从 TEC 网格数据中插值时，并不需要高度信息，所以这里穿刺点位置 posp[2] 没有赋值
}

5. dataindex()：获取 TEC 格网数据下标

先判断点位是否在格网中，之后获取网格点的 tec 数据在 tec.data 中的下标

static int dataindex(int i, int j, int k, const int *ndata) //(i:lat,j:lon,k:hgt)
{
    if (i<0||ndata[0]<=i||j<0||ndata[1]<=j||k<0||ndata[2]<=k) return -1;
    return i+ndata[0]*(j+ndata[1]*k);
}

6. interptec()：插值计算穿刺点处 TEC

通过在经纬度网格点上进行双线性插值，计算第 k 个高度层时穿刺点处的电子数总量 TEC
$\operatorname{TEC}\left(t, \phi_{I P P}, \lambda_{I P P}\right)=\frac{\left(t-t_{i}\right) T E C\left(t_{i}, \phi_{I P P}, \lambda_{I P P}+\omega\left(t-t_{i}\right)\right)+\left(t_{i+1}-t\right) T E C\left(t_{i+1}, \phi_{I P P}, \lambda_{I P P}+\omega\left(t-t_{i+1}\right)\right)}{t_{i+1}-t_{i}}$

static int interptec(const tec_t *tec, int k, const double *posp, double *value,
                     double *rms)
{
    double dlat,dlon,a,b,d[4]={0},r[4]={0};
    int i,j,n,index;
    
    *value=*rms=0.0;    // 将 value和 rms所指向的值置为 0
    
    if (tec->lats[2]==0.0||tec->lons[2]==0.0) return 0;
    
    // 将穿刺点的经纬度分别减去网格点的起始经纬度，再除以网格点间距，对结果进行取整，
    // 得到穿刺点所在网格的序号和穿刺点所在网格的位置(比例) i,j 
    dlat=posp[0]*R2D-tec->lats[0];
    dlon=posp[1]*R2D-tec->lons[0];
    if (tec->lons[2]>0.0) dlon-=floor( dlon/360)*360.0; /*  0<=dlon<360 */
    else                  dlon+=floor(-dlon/360)*360.0; /* -360<dlon<=0 */
    
    a=dlat/tec->lats[2];
    b=dlon/tec->lons[2];
    i=(int)floor(a); a-=i;
    j=(int)floor(b); b-=j;
    
    // 调用 dataindex() 函数分别计算这些网格点的 tec 数据在 tec.data中的下标，
    // 按从左下到右上的顺序
    // 从而得到这些网格点处的 TEC 值和相应误差的标准差
    /* get gridded tec data */
    for (n=0;n<4;n++) {
        if ((index=dataindex(i+(n%2),j+(n<2?0:1),k,tec->ndata))<0) continue;
        d[n]=tec->data[index];
        r[n]=tec->rms [index];
    }
    if (d[0]>0.0&&d[1]>0.0&&d[2]>0.0&&d[3]>0.0) {
        // 穿刺点位于网格内，使用双线性插值计算出穿刺点的 TEC 值
        /* bilinear interpolation (inside of grid) */
        *value=(1.0-a)*(1.0-b)*d[0]+a*(1.0-b)*d[1]+(1.0-a)*b*d[2]+a*b*d[3];
        *rms  =(1.0-a)*(1.0-b)*r[0]+a*(1.0-b)*r[1]+(1.0-a)*b*r[2]+a*b*r[3];
    }
    // 穿刺点不位于网格内，使用最邻近的网格点值作为穿刺点的 TEC 值，不过前提是网格点的 TEC>0
    /* nearest-neighbour extrapolation (outside of grid) */
    else if (a<=0.5&&b<=0.5&&d[0]>0.0) {*value=d[0]; *rms=r[0];}
    else if (a> 0.5&&b<=0.5&&d[1]>0.0) {*value=d[1]; *rms=r[1];}
    else if (a<=0.5&&b> 0.5&&d[2]>0.0) {*value=d[2]; *rms=r[2];}
    else if (a> 0.5&&b> 0.5&&d[3]>0.0) {*value=d[3]; *rms=r[3];}
    // 否则，选用四个网格点中 >0 的值的平均值作为穿刺点的 TEC 值
    else {
        i=0;
        for (n=0;n<4;n++) if (d[n]>0.0) {i++; *value+=d[n]; *rms+=r[n];}
        if(i==0) return 0;
        *value/=i; *rms/=i;
    }
    return 1;
}