[GAMP学习笔记]计算STEC程序中遇到的一些问题小结

计算STEC程序中遇到的一些问题小结

1. 相关知识背景

1.1 穿刺点及电离层投影函数原理

1.2 根据接收机和卫星xyz坐标（方法一）计算穿刺点：

首先将接收机和卫星xyz坐标转化成经纬度，再根据以下函数计算卫星高度角，利用GAMP中计算穿刺点，再插值得到VTEC，经过投影函数计算得到STEC。

//readsat_xyzpos、readrec_xyzpos分别为卫星、接收机xyz坐标
ecef2pos(readsat_xyzpos, readsat_llhpos);
ecef2pos(readrec_xyzpos, readrec_llhpos);//xyz转llh

geodist(readsat_xyzpos, readrec_xyzpos, e);//compute geometric distance and receiver-to-satellite unit vector

satazel(readrec_llhpos, e, azel);//satellite azimuth/elevation angle
ele = azel[1] * R2D;//卫星高度角

得到卫星高度角后，再结合接收机位置（llh）、地球半径、电离层单层高度（标准投影函数通常取450km），利用ionppp计算穿刺点位置（llh）

mf = ionppp(pos, azel, nav->tec->rb, hion, posp);//返回值为标准投影函数值

其中，ionppp函数的各项参数说明如下：

/* ionospheric pierce point position -------------------------------------------
* compute ionospheric pierce point (ipp) position and slant factor
* args   : double *pos      I   receiver position {lat,lon,h} (rad,m)
*          double *azel     I   azimuth/elevation angle {az,el} (rad)
*          double re        I   earth radius (km)
*          double hion      I   altitude of ionosphere (km)
*          double *posp     O   pierce point position {lat,lon,h} (rad,m)
* return : slant factor
* notes  : see ref [2], only valid on the earth surface
*          fixing bug on ref [2] A.4.4.10.1 A-22,23
*-----------------------------------------------------------------------------*/
extern double ionppp(const double *pos, const double *azel, double re,
                     double hion, double *posp)

根据当前时刻与读取的每幅GIM图时刻进行比对，找出对应时刻的GIM图位置：

for (itt = 0; itt < nav->nt; itt++) 
        {
            if (timediff(nav->tec[itt].time, time) > 0.0) break;
        }// itt判断当前用哪张GIM map
        
        if (itt == 0 || itt >= nav->nt)
        {
            printf("%s: tec grid out of period\n", time_str(time, 0));
            return 0;
        }

        if ((tt = timediff(nav->tec[itt].time, nav->tec[itt - 1].time)) == 0.0)
        {
            printf("tec grid time interval error\n");
            return 0;
        }

利用找出的GIM对应的tec、穿刺点位置posp可以计算出vtec，此处分别计算当前GIM图+下一张GIM图，根据时间进行线性插值的方法，得到最终vtec：

stat[0] = interptec(nav->tec + itt - 1, 0, posp21, vtec2, rms2);
stat[1] = interptec(nav->tec + itt, 0, posp22, vtec2 + 1, rms2 + 1);
      
if (stat[0] && stat[1]) { /* linear interpolation by time */
	a = timediff(time, nav->tec[itt - 1].time) / tt;
	vtec22 = vtec2[0] * (1.0 - a) + vtec2[1] * a;

	}
else if (stat[0]) { /* nearest-neighbour extrapolation by time */
	vtec22 = vtec2[0];
	}
else {
	vtec22 = vtec2[1];
}

ionppp函数中，计算穿刺点的方法在rtklib手册上有详细介绍：

$$\alpha =z-z^{{}^{\prime }}$$

$${\varphi }_{ipp}=arcsin(cos\alpha sin{\varphi }_r+sin\alpha cos{\varphi }_{r}cosA{z}_r^s)$$

$$({\varphi }_r>70^{°}andtan\alpha cosA{z}_r^s>tan(\pi /2-{\varphi }_r))or({\varphi }_r<-70^{°}and-tan\alpha cosA{z}_r^s>tan(\pi /2+{\varphi }_r)))$$

$${\lambda }_{ipp}={\lambda }_r+\pi -arcsin\frac{sin\alpha sinA{z}_r^s}{cos{\varphi }_{ipp}}$$

$$(otherwise)$$

$${\lambda }_{ipp}={\lambda }_r+\pi +arcsin\frac{sin\alpha sinA{z}_r^s}{cos{\varphi }_{ipp}}$$

1.3 根据三角关系（方法二）计算穿刺点

直接根据接收机和卫星的xyz坐标，根据三角关系计算接收机-穿刺点的距离，求出穿刺点zyz坐标，再将其转到llh坐标。其原理主要为：

$$sta=\sqrt{{X_R}^2+{Y_R}^2+{Z_R}^2}$$

$$es=\sqrt{{X_S}^2+{Y_S}^2+{Z_S}^2}$$

$$ss=\sqrt{({X_R-X_S)}^2+({Y_R-Y_S)}^2+({Z_R-Z_S)}^2}$$

$$scale=\frac{sta\ast sin(\pi /2-ele-z^{{}^{\prime }})/sin{z}^{{}^{\prime }}}{ss}$$（正弦定理）

对应高度角计算方法为：

$$cosele=\frac{s{s}^2+st{a}^2-e{s}^2}{2\ast ss\ast sta}$$

程序中对应部分为：

/* ---------------calculate ionospheric pierce point position -----------------
* calculate ionospheric  pierce point position
* args   : double   *recpos_xyz       I   receiver position in xyz
*                                   
*          double  *satpos_xyz        I  satellite position in xyz
* 
*          double   SHELL_H           I  ionosphere effective height(m)
*                                 
*          double    *ipp_llh         O    ipp position{lat,lon, height}
*
*          double    *zen_ele         O   zenith and (elevation+2/pi)
* return : none
* notes  : from PANDA-FORTRAN 2022-05-05 by hr
*-----------------------------------------------------------------------------*/
extern double cal_ipp( double* recpos_xyz, double* satpos_xyz, double SHELL_H, double* ipp, double* zen_ele)
{
    double ssLength, esLength, staHeight;
    double cosE, elev, sinZ, zenith;
    double dx, dy, dz;
    double scale, staLayerLength;
    double ipp_xyz[3] = {0};

    dx = satpos_xyz[0] - recpos_xyz[0];
    dy = satpos_xyz[1] - recpos_xyz[1];
    dz = satpos_xyz[2] - recpos_xyz[2];

    ssLength = sqrt(dx * dx + dy * dy + dz * dz);

    esLength = sqrt(satpos_xyz[0] * satpos_xyz[0] + satpos_xyz[1] * satpos_xyz[1] + satpos_xyz[2] * satpos_xyz[2]);
    
    staHeight = sqrt(recpos_xyz[0] * recpos_xyz[0] + recpos_xyz[1] * recpos_xyz[1] + recpos_xyz[2] * recpos_xyz[2]);

    //COMPUTE OF COS OF ELEVATION IS GEOCCENTRIC FRAME
    cosE = (ssLength * ssLength + staHeight * staHeight - esLength * esLength) / (2.0 * ssLength * staHeight);

    //ELEV SHOULD BETWEEN(0, PI)
    elev = acos(cosE);

    //GEOCCENTRIC IPP ZENITH
    double sinele = sin(elev);
    double factor = staHeight / (staHeight + SHELL_H);
    sinZ = sin(elev) * staHeight / (staHeight + SHELL_H);
    zenith = asin(sinZ);

    staLayerLength = staHeight * sin(PI - elev - zenith) / sinZ;

    scale = staLayerLength / ssLength;

    dx = (satpos_xyz[0] - recpos_xyz[0]) * scale;
    dy = (satpos_xyz[1] - recpos_xyz[1]) * scale;
    dz = (satpos_xyz[2] - recpos_xyz[2]) * scale;

    ipp_xyz[0] = recpos_xyz[0] + dx;
    ipp_xyz[1] = recpos_xyz[1] + dy;
    ipp_xyz[2] = recpos_xyz[2] + dz;

    ecef2pos(ipp_xyz, ipp);

    zen_ele[0] = zenith;
    zen_ele[1] = elev;

 /*   double sum = (zenith + elev) * R2D;
    double ipp_lat = ipp[0] * R2D;
    double ipp_lon = ipp[1] * R2D;*/


    return 0;

}

再根据前述方法插值计算vtec

2. 关于两种不同方法计算穿刺点得到的高度角、穿刺点经纬度的差异

以auck站2014年002天10号卫星，自19950-33150s弧段内数据为例：

此处，方法一计算结果记为方法1，方法二计算结果记为方法2。

可以看出两种方法计算高度角差异rms在0.1°左右，经纬度计算差异rms值在0.01°以下，认为两种方法计算差异不大。

3. 计算VTEC考虑经度改正对结果的影响

3.1 根据穿刺点通过插值计算GIM 中VTEC方法

大部分程序中采用的方法为：分别计算当前时刻对应GIM和下一张GIM中VTEC，再通过时间长度进行线性插值：

$$E(\beta ,\lambda ,t)=\frac{T_{i+1}-t}{T_{i+1}-Ti}{E}_i(\beta ,\lambda )+\frac{t-T_i}{T_{i+1}-Ti}{E}_{i+1}(\beta ,\lambda )$$

其中，$\beta ,\lambda$分别为穿刺点的纬度、经度，$t$为当前插值时刻，$E_i(\beta ,\lambda )$表示第$i$幅GIM图中对应位置计算出的VTEC：$T_i\le t\le T_{i+1}$。

考虑到当前插值时间在两幅GIM图时间间隔内，在计算VTEC时，也需要考虑其和太阳的相对位置变化，即需要对插值计算时穿刺点的经度进行改正：

$$E(\beta ,\lambda ,t)=\frac{T_{i+1}-t}{T_{i+1}-Ti}{E}_i(\beta ,{\lambda }_i^{{}^{\prime }})+\frac{t-T_i}{T_{i+1}-Ti}{E}_{i+1}(\beta ,{\lambda }_{i+1}^{{}^{\prime }})$$

其中，${\lambda }_i^{{}^{\prime }}=\lambda +(t-T_i)$。

3.2 两种方法计算VTEC的造成的差异

此处，仍以auck站2014年002天10号卫星，自19950-33150s弧段内数据为例，不考虑和考虑旋转改正分别记为NoRotaCorr和WithRotaCorr。

可以发现当考虑经度改正时，在一个连续弧段内计算得到的VTEC出现了不连续的情况，如STEC差异图中的小图所示，不考虑经度改正时，计算得到的STEC连续，而考虑经度改正后在21990时刻出现约0.6TECu大小的跳动。经分析主要原因如下（以10号卫星在21960、21990、22020三个历元为例）：

在21990时刻，$\lambda =178.41$，考虑经度改正后${\lambda }_i^{{}^{\prime }}$和${\lambda }_{i+1}^{{}^{\prime }}$分别为180.04和176.29，此时${\lambda }_i^{{}^{\prime }}$已经超出了该穿刺点对应格网的范围，即采用经度改正后的穿刺点计算的VTEC为相邻格网的数据，故此时计算出的VTEC最终结果与前一历元相比出现较大跳动。