前言：
本章节代码均在Gitee中开源：
电离层对流层延迟解算
https://gitee.com/Ehundred/navigation-engineering/tree/master/%E5%8D%AB%E6%98%9F%E5%AF%BC%E8%88%AA%E5%8E%9F%E7%90%86/%E7%94%B5%E7%A6%BB%E5%B1%82%E5%AF%B9%E6%B5%81%E5%B1%82%E8%AF%AF%E5%B7%AE%E8%A7%A3%E7%AE%97其中涉及到时空转换和卫星位置解算的部分，在另一章节有详细讲解，如有疑惑可以先阅读上一章节：
卫星位置解算-CSDN博客https://blog.csdn.net/qq_74260823/article/details/139361271?spm=1001.2014.3001.5501因为这章是学校作业，所以稍微正经点.
电离层和对流层
电离层和对流层是什么？
粒子的电离
经历过九年义务教育的都知道，空气中含有氮气氧气还有小部分其他气体，虽然成分复杂，但是他们都是稳定的中性气体，不会发生剧烈的反应(不然我们也不知道再多看一眼天空就会爆炸)。不过，会有一个因素导致他们发生电离——太阳辐射。
高中物理学过，电子接收到能量，会由基态跃迁到激发态，当能量达到一定阈值的时候，电子就会被电离出去。没学过也没关系，只需要知道一件事——一个中性的粒子，在接受到足够多的能量之后，就会变得不再中性，转而变成带电粒子。
那么，地球的主要能量来自于哪里？太阳辐射。粒子在接受到太阳辐射之后，吸收到了足够多的能量，电子就会跃迁，中性粒子也就变成了带电粒子。还不理解可以想想太阳能电力板，电力板吸收了足够的太阳辐射之后，便会产生大量带电粒子形成电流，产生电能。
一层空气的电离强度，实际上是什么？是带电粒子的密度。而且，带电粒子的密度有哪些影响因素？粒子的稀薄程度，和吸收能量的多少。
- 当中性粒子不够多，当所有粒子都变为带电粒子，就算给再多能量，也无法再产生更多的带电粒子
- 当能量不够多，就算有再多的粒子，也只会有固定数量的粒子吸收到能量电离出去，剩余的粒子因为没有足够的能量而保持中性
所以，电离强度的基本影响因素也只有两个：
- 中性粒子的密度
- 太阳辐射的强度
太阳辐射和电离强度
太阳辐射下来，其能量肯定会有一个规律：
越高处的粒子先接触到太阳辐射，先进行能量吸收，而低处的粒子接收到的能量，是高处粒子吸收剩下来的能量。
所以电离强度随高度和时间的变化图一定是这样：
- 白天比夜晚大的原因：白天有太阳，太阳辐射高，总能量大
- 随着高度升高，先升的原因：被吸收的太阳辐射更少，剩下的太阳辐射更大
- 后降的原因：空气的大气变得稀薄，中性粒子减少，已经全部转换为带点粒子
电离层和大气层
OK，以上说了这么多，只想说明一件事：电离强度会在大约50-400km处到达顶峰。
而这一个区间，带电粒子的影响及其剧烈，我们便把他称为电离层
而在电离层以下，50km以下，带电粒子的影响变得很小，但是因为在这个高度，有着平流层对流层，大气中也开始充满水蒸气，这些因素也会造成较大的影响，我们便把这一段称为大气层。
电离层和大气层都会较为明显影响信号波的传输，但是他们的模型却有极大的区别，所以我们把他们区分开来，就有了电离层和大气层。
为什么会有电离层和对流层延迟？
卫星信号在发射到接收机接收的时候，会经过两个阶段：
- 在太空的真空环境传播
- 在大气内的环境传播
而在大气内，又会经过两个阶段：
- 50km以上，电离影响及其剧烈，极大减慢了卫星信号波的运动速度，造成信号传输的延时
- 50km以下，受到对流层的气流和水蒸气影响，造成卫星信号传播路径的弯曲和速度的降低，也导致了信号传输的延时
而这两个影响因素，便分别构成了电离层延迟和对流层延迟。
其中，在电离层中，信号的传播速度与频率有关，而在对流层中则与频率无关。
但是，卫星信号发射到接收机上，肯定不是与电离层垂直的，其可能为下图的情况：
我们实际计算电离层和对流层误差的时候， 就必须把他转换到垂直的方向上计算：
电离层对流层延迟解算
电离层和对流层延迟解算，基本思想都是模型改正。而这些模型大部分都是经验模型，在一般的工程函数里，都会给出一个函数接口：
但是，这里的azel，在广播星历中其实并没有，是由我们人工计算出来的，虽然大部分都没有直接给出计算方法。因为在电离层和对流层都会用到这个azel，所以我们直接写出一个基类，然后让电离层和对流层都去继承这个基类：
class base_data
{
protected:
satellite_position_calculation rs_place;//解算卫星坐标
Position<Space_Coordinates> _rs;//卫星位置
Position<Space_Coordinates> _rr;//接收站位置
double _d;//卫星到接收站的距离
vector<double> _e;//卫星到接收站的单位向量
vector<double> _enu;//卫星在站星坐标系下的位置ENU
double _az;//方位角
double _el;//高度角
double _central_angle;//中心角
};
什么是高度角和方位角？
高度角和方位角都是以接收站为原点的概念。 高度角就是表示卫星高度z轴与水平面的角度，而方位角则是表示卫星在xy平面上的偏向。
而计算高度角和方位角的基本计算思路，便是把卫星转换到站星坐标系上。卫星和接收站的坐标，都是在地心地固坐标系下的坐标，通过这两个坐标，我们可以把卫星的坐标转换到以接收站为原点的坐标系上，方便计算。
1.计算接收机到卫星之间的几何距离
double Distance_from_receiver_to_satellite()
{
double tmp = 0;
for (int i = 0; i < 3; i++)
tmp += pow((_rs.Get__Space_Coordinates()[i] - _rr.Get__Space_Coordinates()[i]), 2);
_d = sqrt(tmp);
return _d;
}2.视距单位向量
vector<double> Unit_vector()
{
_e.resize(3);
for (int i = 0; i < 3; i++)
{
_e[i] = (_rs.Get__Space_Coordinates()[i] - _rr.Get__Space_Coordinates()[i]) / _d;
}
return _e;
}3.地心坐标系转站心坐标系
vector<double> rs_topocentric()
{
_enu.resize(3);
double sinL = sin(_rr.Get__Geodetic_Coordinate_System()[1] * Pi / 180);
double cosL = cos(_rr.Get__Geodetic_Coordinate_System()[1] * Pi / 180);
double sinB = sin(_rr.Get__Geodetic_Coordinate_System()[0] * Pi / 180);
double cosB = cos(_rr.Get__Geodetic_Coordinate_System()[0] * Pi / 180);
auto matrix_multiplication = [this](double factor1, double factor2, double factor3)
{
return factor1 * (_rs.Get__Space_Coordinates()[0] - _rr.Get__Space_Coordinates()[0])
+ factor2 * (_rs.Get__Space_Coordinates()[1] - _rr.Get__Space_Coordinates()[1])
+ factor3 * (_rs.Get__Space_Coordinates()[2] - _rr.Get__Space_Coordinates()[2]);
};
_enu[0] = matrix_multiplication(-sinL, cosL, 0);
_enu[1] = matrix_multiplication(-sinB * cosL, -sinB * sinL, cosB);
_enu[2] = matrix_multiplication(cosB * cosL, cosB * sinL, sinB);
return _enu;
}4.计算卫星高度角和方位角
double azimuth()
{
_az = atan2(_enu[0],_enu[1]);
return _az;
}
double Height_angle()
{
_el = asin(_enu[2] / _d);
return _el;
}5.计算地球中心角
double central_angle()
{
_central_angle = 0.0137 / ((_el / Pi) + 0.11) - 0.022;
return _central_angle;
}代码汇总
对于一组卫星和接收机，因为只是做了坐标系的转换，所以这些基本参数是固定的。与之前的思想一致，我们直接在构造函数中调用这些函数，完成初始化：
class base_data
{
protected:
base_data(int Serial, const vector<double>& rr)
:rs_place(Serial),
_rs(rs_place.xyz[0], rs_place.xyz[1], rs_place.xyz[2]),
/*_rs(12712882.254 ,23247798.196 ,-2637709.427),*/
_rr(rr[0], rr[1], rr[2])
{
Distance_from_receiver_to_satellite();
Unit_vector();
rs_topocentric();
azimuth();
Height_angle();
central_angle();
}
satellite_position_calculation rs_place;//解算卫星坐标
Position<Space_Coordinates> _rs;//卫星位置
Position<Space_Coordinates> _rr;//接收站位置
double _d;//卫星到接收站的距离
vector<double> _e;//卫星到接收站的单位向量
vector<double> _enu;//卫星在站星坐标系下的位置ENU
double _az;//方位角
double _el;//高度角
double _central_angle;//中心角
double Distance_from_receiver_to_satellite()
{
double tmp = 0;
for (int i = 0; i < 3; i++)
tmp += pow((_rs.Get__Space_Coordinates()[i] - _rr.Get__Space_Coordinates()[i]), 2);
_d = sqrt(tmp);
return _d;
}
vector<double> Unit_vector()
{
_e.resize(3);
for (int i = 0; i < 3; i++)
{
_e[i] = (_rs.Get__Space_Coordinates()[i] - _rr.Get__Space_Coordinates()[i]) / _d;
}
return _e;
}
vector<double> rs_topocentric()
{
_enu.resize(3);
double sinL = sin(_rr.Get__Geodetic_Coordinate_System()[1] * Pi / 180);
double cosL = cos(_rr.Get__Geodetic_Coordinate_System()[1] * Pi / 180);
double sinB = sin(_rr.Get__Geodetic_Coordinate_System()[0] * Pi / 180);
double cosB = cos(_rr.Get__Geodetic_Coordinate_System()[0] * Pi / 180);
auto matrix_multiplication = [this](double factor1, double factor2, double factor3)
{
return factor1 * (_rs.Get__Space_Coordinates()[0] - _rr.Get__Space_Coordinates()[0])
+ factor2 * (_rs.Get__Space_Coordinates()[1] - _rr.Get__Space_Coordinates()[1])
+ factor3 * (_rs.Get__Space_Coordinates()[2] - _rr.Get__Space_Coordinates()[2]);
};
_enu[0] = matrix_multiplication(-sinL, cosL, 0);
_enu[1] = matrix_multiplication(-sinB * cosL, -sinB * sinL, cosB);
_enu[2] = matrix_multiplication(cosB * cosL, cosB * sinL, sinB);
return _enu;
}
double azimuth()
{
_az = atan2(_enu[0],_enu[1]);
return _az;
}
double Height_angle()
{
_el = asin(_enu[2] / _d);
return _el;
}
double central_angle()
{
_central_angle = 0.0137 / ((_el / Pi) + 0.11) - 0.022;
return _central_angle;
}
};电离层
因为电离层主要采用的是模型修正的方法，其计算步骤我们就不太去深入。其思想也很简单：因为信号进入大气层的时候，并非是垂直进入的，我们必须将其分为两个向量，只考虑垂直方向上向量的影响，然后求出垂直方向向量进入大气层的经纬度，再采用求出的经纬度进行模型改正，便可以求出电离层延迟。
1.计算电离层穿刺纬度
注意：
在这里的pos是接收机的大地坐标。因为大地坐标的经纬度可以用角度表示，也可以用弧度表示，但是其计算方式是不一样的。如果用角度表示，即0~180度，则这里的Pi应该取180；而如果用弧度表示，则这里的Pi应该取3.14
double FAIi()
{
_FAIi = _rr.Get__Geodetic_Coordinate_System()[0] / 180 + _central_angle * cos(_az);
if (_FAIi > 0.416) _FAIi = 0.416;
if (_FAIi < -0.416) _FAIi = -0.416;
return _FAIi;
}2.计算电离层穿刺经度
同计算纬度时一样，最后一个式子pos是接收机的大地坐标，而Pi的取值看大地坐标是弧度还是角度
double Lambdai()
{
_Lambdai = _rr.Get__Geodetic_Coordinate_System()[1] / 180 +
(_central_angle * sin(_az)) / cos(_FAIi * Pi);
return _Lambdai;
}3.计算电离层穿刺点地理纬度
double FAIm()
{
_FAIm = _FAIi + 0.064 * cos((_Lambdai - 1.617) * Pi);
return _FAIm;
}4.计算电离层穿刺点当地时间
double get_t()
{
double tmp = 43200 * _Lambdai + curtime.GetGPStime()._second;
_t = tmp - floor(tmp / 86400) * 86400;
return _t;
}//curtime 为观测时间5.计算电离层延迟幅度
在一些函数接口中，α和β系列参数都会由数组给出。但是这些参数，其实在广播星历头文件中给出了。用什么方式实现都可以，我这里是用了宏定义来实现：
#define alpha0 0.8382E-8
#define alpha1 -0.7451E-8
#define alpha2 -0.596E-7
#define alpha3 0.596E-7
double Ai()
{
_Ai = alpha0 + _FAIm * (alpha1 + _FAIm * (alpha2 + _FAIm * alpha3));
if (_Ai < 0) _Ai = 0;
return _Ai;
}6.计算电离层延迟周期
#define Beta0 0.8806E5
#define Beta1 -0.3277E5
#define Beta2 -0.1966E6
#define Beta3 0.1966E6
double cyclei()
{
_cyclei = Beta0 + _FAIm * (Beta1 + _FAIm * (Beta2 + _FAIm * Beta3));
if (_cyclei < 72000) _cyclei = 72000;
return _cyclei;
}7.计算电离层延迟相位
double Xi()
{
_Xi = 2 * Pi * (_t - 50400) / _cyclei;
return _Xi;
}8.计算倾斜因子
double obliquity_factor()
{
_F = 1.0 + 16.0 * pow((0.53 - _el / Pi), 3);
return _F;
}9.计算电离层延迟
double Ionospheric_delay()
{
if (fabs(_Xi) >= 1.57)
return 5e-9 * clight * _F;
else
return clight * (5E-9 + _Ai * (1 - 0.5 * _Xi * _Xi + pow(_Xi, 4) / 24.0)) * _F;
}代码汇总：
class iono_calculate:public base_data
{
public:
iono_calculate(int Serial, const vector<double>& rr)
:base_data(Serial,rr)
{}
double iono_delay()
{
return _iono_calculate();
}
private:
double _FAIi;//电离层穿刺点的纬度
double _Lambdai;//电离层穿刺点经度
double _FAIm;//电离层穿刺点的地磁纬度
double _t;//GPS周内秒
double _Ai;//电离层延迟的幅度
double _cyclei;//电离层延迟的周期
double _Xi;//电离层延迟的相位
double _F;//倾斜因子
double FAIi()
{
_FAIi = _rr.Get__Geodetic_Coordinate_System()[0] / 180 + _central_angle * cos(_az);
if (_FAIi > 0.416) _FAIi = 0.416;
if (_FAIi < -0.416) _FAIi = -0.416;
return _FAIi;
}
double Lambdai()
{
_Lambdai = _rr.Get__Geodetic_Coordinate_System()[1] / 180 +
(_central_angle * sin(_az)) / cos(_FAIi * Pi);
return _Lambdai;
}
double FAIm()
{
_FAIm = _FAIi + 0.064 * cos((_Lambdai - 1.617) * Pi);
return _FAIm;
}
double get_t()
{
double tmp = 43200 * _Lambdai + curtime.GetGPStime()._second;
_t = tmp - floor(tmp / 86400) * 86400;
return _t;
}
double Ai()
{
_Ai = alpha0 + _FAIm * (alpha1 + _FAIm * (alpha2 + _FAIm * alpha3));
if (_Ai < 0) _Ai = 0;
return _Ai;
}
double cyclei()
{
_cyclei = Beta0 + _FAIm * (Beta1 + _FAIm * (Beta2 + _FAIm * Beta3));
if (_cyclei < 72000) _cyclei = 72000;
return _cyclei;
}
double Xi()
{
_Xi = 2 * Pi * (_t - 50400) / _cyclei;
return _Xi;
}
double obliquity_factor()
{
_F = 1.0 + 16.0 * pow((0.53 - _el / Pi), 3);
return _F;
}
double Ionospheric_delay()
{
if (fabs(_Xi) >= 1.57)
return 5e-9 * clight * _F;
else
return clight * (5E-9 + _Ai * (1 - 0.5 * _Xi * _Xi + pow(_Xi, 4) / 24.0)) * _F;
}
double _iono_calculate()
{
FAIi();
Lambdai();
FAIm();
get_t();
Ai();
cyclei();
Xi();
obliquity_factor();
return Ionospheric_delay();
}
};调试数据：
同样，为了方便大家调试，给出一个测试用例：
α和β为上面给出的值：
#define alpha0 0.8382E-8
#define alpha1 -0.7451E-8
#define alpha2 -0.596E-7
#define alpha3 0.596E-7
#define Beta0 0.8806E5
#define Beta1 -0.3277E5
#define Beta2 -0.1966E6
#define Beta3 0.1966E6最终计算出的值为：
对流层
对流层则更为简单，直接计算出大气参数，然后分别计算干延迟和湿延迟，最后把干延迟和湿延迟相加，便是对流层延迟的结果。
因为步骤就两步，我就不分开为每一个具体的函数了，我们直接用一个函数：
这里e-5不是表示数学中的自然底数，而是十进制中的E-5，表示10的-5次方
而pos也是表示接收站的大地坐标
void calculate()
{
_T = 15.0 - 6.5E-3 * _rr.Get__Geodetic_Coordinate_System()[2] + 273.16;
_e = 6.108 * humi * exp((17.15 * _T - 4684.0) / (_T - 38.45));
_p = 1013.25 * pow((1 - 2.2557 * 1e-5 * _rr.Get__Geodetic_Coordinate_System()[2]), 5.2568);
double z = Pi / 2 - _el;
_trph = 0.0022768 * _p / (cos(z) * (1.0 - 0.00266 *
cos(2 * _rr.Get__Geodetic_Coordinate_System()[0] * Pi / 180)
- 0.00028 * _rr.Get__Geodetic_Coordinate_System()[2] / 1000));
_trpw = 0.002277 * (1255.0 / _T + 0.05) * _e / cos(z);
_trp = _trph + _trpw;
}代码汇总：
#define humi 0.7//相对湿度
class tropo_calculate:public base_data
{
public:
tropo_calculate(int Serial, const vector<double>& rr)
:base_data(Serial,rr)
{}
double tropo_delay()
{
calculate();
return _trp;
}
private:
double _p;//大气总压强
double _T;//空气绝对温度
double _e;//偏心率
double _trph;//对流层干延迟
double _trpw;//对流层湿延迟
double _trp;//对流层总延迟
void calculate()
{
_T = 15.0 - 6.5E-3 * _rr.Get__Geodetic_Coordinate_System()[2] + 273.16;
_e = 6.108 * humi * exp((17.15 * _T - 4684.0) / (_T - 38.45));
_p = 1013.25 * pow((1 - 2.2557 * 1e-5 * _rr.Get__Geodetic_Coordinate_System()[2]), 5.2568);
double z = Pi / 2 - _el;
_trph = 0.0022768 * _p / (cos(z) * (1.0 - 0.00266 *
cos(2 * _rr.Get__Geodetic_Coordinate_System()[0] * Pi / 180)
- 0.00028 * _rr.Get__Geodetic_Coordinate_System()[2] / 1000));
_trpw = 0.002277 * (1255.0 / _T + 0.05) * _e / cos(z);
_trp = _trph + _trpw;
}
};调试数据：
接收机和卫星坐标还是和电离层一样：
其调试结果为：
更多的测试用例可以见Gitee上的完整代码。
最后，给自己叠个甲。因为自己才是导航工程大二的本科生，有些概念理解可能不到位，而又想用最容易理解的方式表达出来，所以可能正确性会稍微有些偏差。但是对初学者来说，应该不会存在太大的错误，如果可以帮到你，真的荣幸之极。还有
