GNSS/C/C++——SPP中的大气延迟
作者：互联网
大气延迟
对流层延迟
在GNSS领域，对流层延迟是定位误差来源之一。卫星导航定位中的对流层延迟通常是泛指电磁波信号在通过高度为50km以下的未被电离的中性大气层时所产生的信号延迟。在研究信号延迟的过程中，我们不再将该大气层细分为对流层和平流层（如大气科学中那样），也不再顾及两者之间性质上的差别。
由于80%的延迟发生在对流层，所以我们将发生在该中性大气层中的信号延迟通称为对流层延迟。 对流层是大气层较低的部分，对于直到高达15GHz的频率来说它是非色散的。在这种介质中，与L1和L2上GPS载波和信号信息（PRN码和导航数据）相关联的相速和群速，相当于自由空间传播被同等地延迟了。这种延迟随对流层折射率而变，其折射率取决于当地的温度、压力和相对湿度。
如果不补偿，这种延迟的等效距离可从卫星在天顶和用户在海平面上的2.4m左右到卫星在约5°仰角上的25m左右。
电离层延迟
电离层概念
电离层（Ionosphere）是地球大气的一个电离区域。它是受到太阳高能辐射以及宇宙线的激励而电离的大气高层。50千米以上的整个地球大气层都处于部分电离或完全电离的状态，电离层是部分电离的大气区域，完全电离的大气区域称磁层。
- 电离层的范围从离地面约50公里开始一直伸展到约1000公里高度的地球高层大气空域。
- 电离层的主要特性由电子密度、电子温度、碰撞频率、离子密度、离子温度和离子成分等空间分布的基本参数来表示。
- 电离层的研究对象主要是电子密度随高度的分布。电子密度（或称电子浓度）是指单位体积的自由电子数，随高度的变化与各高度上大气成分、大气密度以及太阳辐射通量等因素有关。电离层内任一点上的电子密度，决定于上述自由电子的产生、消失和迁移三种效应。在不同区域，三者的相对作用和各自的具体作用方式也大有差异。
- 电离层中存在相当多的自由电子和离子，能使无线电波改变传播速度，发生折射、反射和散射，产生极化面的旋转并受到不同程度的吸收。
电离层延迟误差
电离层延迟误差（ionospheric delay error）亦称电离层折射误差(ionospheric refraction error)，是由电离层效应引起的观测值误差。
与其他电磁波一样，当导航卫星信号通过电离层时，导航卫星信号的路径也要发生弯曲，传播速度也会发生变化。电离层延迟对距离测量的影响，在天顶方向(仰角0度)最大可达50m，在近地平方向时(仰角20度)可达150m，因此必须仔细地加以改正，否则将严重影响伪距观测精度。
电离层折射可用三种方法来消除：
一是利用电离层模型加以改正；
二是利用同步观测值求差，这种方法对于短基线的效果尤为明显；
三是利用双频观测值，即通过不同频率的观测值组合来对电离层的延迟进行改正。
SPP中模型改正方法
声明
/*-----------------------------解码模块;-----------------------*/
/*星历结构体;*/
struct GPSEPHEM
{
unsigned long PRN; //Satallite PRN number;
GPSTIME TOC;
GPSTIME TOE;
/*第一数据块(第一子帧）;*/
unsigned long week; //GPS reference week number;
double URN; //User Range Accurancy variance,m2;
unsigned long health; //Health status;
double tgd; //Estimated group delay difference,seconds;
unsigned long idoc; //Issue of data clock;
double a_f0; //Clock aging paramenter,seconds(s);
double a_f1; //Clock aging paramenter(s/s);
double a_f2; //Clock aging paramenter(s/s/s);
/*第二数据块（第二、三子帧）;*/
//星历参考时刻toe时刻的轨道根数;
double M0; //Mean anomaly of reference time,radians;
double w; //Argument of perigee,radians;近地点角距;
double I0; //Inclination angle at reference time,radian;
double w0; //Right ascension,radians;
double A; //Semi-major axis,metres;半长轴;
double ecc; //Eccentricity,dimensionless-quantity defined for a conic section;
double tow; //Time stamp of subframe;子帧的时间标记;
//9个轨道摄动参数;
double delta_N; //Mean motion difference,radians/second;
double V_I; //Rate of inclination angle,radians/second;
double V_w0; //Rate of Right ascension,radians;
double cuc; //Argument of latitude(amplitude of cosine,radians);
double cus; //Argument of latitude(amplitude of sine,radians);
double crc; //Orbit radius(amplitude of cosine,metres);
double crs; //Orbit radius(amplitude of sine,metres);
double cic; //Inclination(amplitude of cosine,radian);
double cis; //Inclination(amplitude of sine,radian);
//包含的参数;
unsigned long IODE1;
unsigned long IODE2;
unsigned long z_week; //Z count week number;
double toe; //Reference time for ephemeris ,seconds;
double toc; //SV clcok correction term,seconds;
int AS; //0=FALSE;1=TRUE;
double N; //Corrected mean motion,radians/seconds;
};
/*Ion与UTC解码结构体;*/
/*Message_ID==8 IONUTC Ionospheric and UTC Data*/
struct IONOSPHERIC
{
double a0; //Alpha parameter constant term;
double a1; //Alpha parameter 1st order term;
double a2; //Alpha parameter 2nd order term;
double a3; //Alpha parameter 3rd order term;
double b0; //Beta parameter constant term;
double b1; //Beta parameter 1st order term;
double b2; //Beta parameter 2nd order term;
double b3; //Beta parameter 3rd order term;
IONOSPHERIC()
{
a0 = a1 = a2 = a3 = b0 = b1 = b2 = b3 = 0;
}
};
struct UTC
{
unsigned long utc_wn; //UTC reference week number;
unsigned long tot; //Reference time of UTC parameters;
double A0; //UTC constant term of polynomial;
double A1; //UTC 1st order term of polynomial;
unsigned long wn_lsf; //Future week number;
unsigned long dn; //Day number(1=Sunday and 7=Saturday);
long deltat_ls; //Delta time due to leap second;
long deltat_lsf; //Future delta time due to leap second;
unsigned long deltat_utc;//Time difference;
};
/*接收机观测值结构体;*/
/*Message_ID==43 RANGE Satellite Range Information*/
struct RANGE
{
long obsnum; //Number of oberservations with information to follow;
unsigned short PRN_slo; //GPS:1-32,SBAS:120-138,GLONASS:38-61;
unsigned short glofreq; //GLONASS Frequency+7;
double psr; //Pseudorange measurement(m);
float psr_std; //Pseudorange measurement standard deviation(m);
double adr; //Carrier phase,in cycles(accumulated Doppler range);
float adr_std; //Estimated carrier Doppler frequency(Hz);
float dopp; //Instantaneous carrier Doppler frequency;
float C_No; //Carrier to nosie density radio;C/No=10[log10(S/No)](dB-Hz);
float locktime; //#of second of continuous tracking(no cycle slipping);
unsigned long ch_tr_status; //Tracking status;
};
/*存储接收机自己解算的位置结构体;*/
/*Message_ID==496 PDPPOS PDP filter position*/
struct PDPPOS
{
int sol_status; //Solution status;
int pos_type; //Position type ;
double lat; //Latitude;
double lon; //Longitude;
double hgt; //Height above mean sea level;海拔高 ellipsoid(m);
float undulation; //Undulation-the relationship between the geoid and the WGS84 ellipsoid;
int datum_id; //Datum ID number;
float lat_delta; //Latitude standard deviation;
float lon_delta; //Longitude standard deviation;
float hgt_delta; //Height standard deviation;
float diff_age; //Differential age in seconds;
float sol_age; //Solution age in seconds;
unsigned char sats; //Number of satellite vehicles tracked;
unsigned char sats_soln; //Number of satellites in the solution;
};
/*每个历元的观测值数据;*/
struct EPOCHOBS
{
GPSTIME Time; //GPS time;
short SATnum[CHANNELNUM]; //which SAT
long Obsnums; //Number of observation with information to follow;
RANGE RangeObs[CHANNELNUM]; //Number of Channels;
unsigned long ReceiverStatus; //Receiver Status;
EPOCHOBS()
{
Obsnums = 0;
ReceiverStatus = -1;
}
};
/*存储所有数据，包括星历，观测值等的结构体;*/
struct RAW
{
EPOCHOBS EpochObs; //每个历元的观测数据;
GPSEPHEM EphemObs[SATNUM]; //卫星星历数据;
IONOSPHERIC IonObs; //电离层数据;
UTC UTCObs; //时间数据;
PDPPOS PDPpos;
};
/*对流层气象参数结构体，包含在SATPOS结构体中;*/
struct TROPOSPHERE
{
double T; //测站上的干温;
double p; //测站上的气压;
double RH; //测站上的相对湿度;
double H; //测站高度;
TROPOSPHERE()
{
T = 0.0; //测站上的干温;
p = 0.0; //测站上的气压;
RH = 0.0; //测站上的相对湿度;
H = 0.0; //测站高度;
}
};
/*存储接收机解算结果的结构体;*/
struct RECIVIERPOS
{
GPSTIME Time; //什么时刻的位置;
XYZ XYZPOS; //空间直角坐标系下的位置;
BLH BLHPOS; //大地坐标系下位置;
double Velocity[3]; //接收机速度;
double delta_t; //接收机钟差;
double Vdelta_t;//接收机钟速;
TROPOSPHERE TropPara;
double PDOP;
double DOP;
double sigmaP;
double sigmaV;
RECIVIERPOS()
{
DOP = 0;
PDOP = 0;
delta_t=0; //接收机钟差;
Vdelta_t=0;//接收机钟速;
Velocity[0] = Velocity[1] = Velocity[2] = 0.0;
sigmaV = sigmaP = 0;
}
};
/*解算卫星星历过程中得到的卫星轨道参数;*/
struct SATORBPAR
{
double n0;
double n;
double tk;
double Mk;
double Ek;
double vk;
double PHIk;
double delta_uk;
double delta_rk;
double delta_ik;
double uk;
double rk;
double ik;
double xk1;
double yk1;
double wk;
double xk;
double yk;
double zk;
double e;
double A;
};
/*存储卫星星历解算结果的结构体;*/
struct SATPOS
{
short PRN;
GPSTIME Time;
XYZ XYZPOS;
BLH BLHPOS;
EAR O2S_EARPOS;
double delta_Ion;
double delta_Trop;
double SignalTime;
double Velocity[3];
double delta_t;
double Vdelta_t;
SATPOS()
{
Velocity[0] = Velocity[1] = Velocity[2] = 0.0;
delta_t = 0;
Vdelta_t = 0;
SignalTime = 0;
delta_Ion = 0;
delta_Trop = 0;
}
};
实现
/*-----------------------电离层对流层误差解算结果;-------------------------*/
//模块功能：用以解算对流层与电离层对定位造成的误差干扰，并给出改正数;
//函数包含：;
// Klobuchar:电离层改正函数;
// Hopefiled:对流层改正函数;
// 公式来源：《GPS测量与数据处理》中关于电离层与对流层改正部分，王老师所给PPT;
/*------------------------------------------------------------------------*/
#include "Fuction.h"
/*-------------------------电离层改正模型;---------------------------------*/
//函数功能：计算电离层对定位的影响;
//函数输入：;
// RawObs:接收机输出数据;
// SatPos:卫星位置结构体;
// RecivierPos:接收机位置结构体;
// delta_Ion:电离层延迟;
//公式来源：《GPS测量与数据处理》中关于电离层与对流层改正部分，王老师所给PPT;
/*--------------------------------------------------------------------------*/
bool Klobuchar(RAW*RawObs, SATPOS *SatPos,RECIVIERPOS *RecivierPos, double *delta_Ion)
{
ENU ENUPos;
XYZ2BLH(RecivierPos->XYZPOS, &RecivierPos->BLHPOS);
XYZ2ENU(RecivierPos->XYZPOS, SatPos->XYZPOS, &ENUPos);
ENU2EAR(ENUPos, &SatPos->O2S_EARPOS);
double psi, phi_I, lambda_I, phi_m, t = 0;
double A_I,P_I,F_I,X_I=0;
//计算测站点在地心的夹角;
psi = 0.0137 / (SatPos->O2S_EARPOS.E/PI + 0.11) - 0.022;
phi_I = RecivierPos->BLHPOS.B / PI + psi*cos(SatPos->O2S_EARPOS.A);
if (phi_I > 0.461)phi_I = 0.461;
if (phi_I < -0.461)phi_I =- 0.461;
//升交点经度;
lambda_I = RecivierPos->BLHPOS.L / PI + (psi*sin(SatPos->O2S_EARPOS.A)) / cos(phi_I*PI);
//地磁维度;
phi_m = phi_I + 0.064*cos((lambda_I - 1.617)*PI);
//计算地方时;
t = 43200 * lambda_I + RawObs->EpochObs.Time.SecOfWeek;
t = t - (int)(t / 86400.0) * 86400;
while (t<0 || t>86400)
{
if (t >= 86400)t = t - 86400;
if (t < 0)t = t + 86400;
}
//振幅;
A_I = RawObs->IonObs.a0 + RawObs->IonObs.a1*phi_m + RawObs->IonObs.a2*phi_m*phi_m +
RawObs->IonObs.a3*phi_m*phi_m*phi_m;
if (A_I < 0)A_I = 0.0;
//周期;
P_I = RawObs->IonObs.b0 + RawObs->IonObs.b1*phi_m + RawObs->IonObs.b2*phi_m*phi_m +
RawObs->IonObs.b3*phi_m*phi_m*phi_m;
if (P_I < 72000)P_I = 72000;
//相位延迟;
X_I = 2 * PI*(t - 50400) / P_I;
//倾斜因子;
F_I = 1.0 + 16.0*(0.53 - SatPos->O2S_EARPOS.E / PI)*(0.53 - SatPos->O2S_EARPOS.E / PI)
*(0.53 - SatPos->O2S_EARPOS.E / PI);
//计算电离层延迟;
if (fabs(X_I) <= 1.57)
{
*delta_Ion = (5e-9 + A_I*cos(X_I))*F_I;
}
else
{
*delta_Ion = 5e-9*F_I;
}
*delta_Ion = *delta_Ion*WGS84_V_Light;
return true;
}
/*-----------------------------Hopefield 模型;----------
函数功能：用于计算对流层改正参数;
函数输入：;
SatPos:卫星位置结构体;
ReciverPos:接收机位置结构体;
Trop:接收机位置气象参考;
delta_Trop:对流层改正;
公式来源：《GPS测量与数据处理》中关于电离层与对流层改正部分，王老师所给PPT;
----------------------------------------------------*/
bool Hopefield(SATPOS *SatPos, RECIVIERPOS ReciverPos, TROPOSPHERE Trop,double *delta_Trop)
{
ENU Obs2Sat;
BLH BLHObsPos;
XYZ2ENU(ReciverPos.XYZPOS, SatPos->XYZPOS, &Obs2Sat);
ENU2EAR(Obs2Sat, &SatPos->O2S_EARPOS);
XYZ2BLH(ReciverPos.XYZPOS, &BLHObsPos);
if (SatPos->O2S_EARPOS.E < 0){ *delta_Trop = 0; return false; }
//标准气象元素;
double H0, T0, p0, RH0 = 0;
H0 = 0;
T0 = 15 + 273.15;//单位是开尔文;
p0 = 1013.25;
RH0 = 0.5;
double RH, p, T, e, hw, hd, Kw, Kd = 0;
Trop.H = BLHObsPos.H;
Trop.RH = RH0*exp(-0.0006396*(Trop.H - H0));
Trop.p = p0*pow(1 - 0.0000226*(Trop.H - H0)*1.0, 5.225);
Trop.T = T0 - 0.0065*(Trop.H - H0);
e = Trop.RH*exp(-37.2465 + 0.213166*Trop.T - 0.000256908*Trop.T*Trop.T);
hw = 11000.0;
hd = 40136 + 148.72*(T0 - 273.16);
Kw = (155.2e-7)*4810.0*e*(hw - Trop.H) / (Trop.T*Trop.T);
Kd = (155.2e-7)*Trop.p*(hd - Trop.H) / Trop.T;
double delta_d, delta_w = 0;
delta_d = Kd / sin(sqrt(pow(SatPos->O2S_EARPOS.E, 2) + 6.25*D2R*D2R ));
delta_w = Kw / sin(sqrt(pow(SatPos->O2S_EARPOS.E, 2) + 2.25*D2R*D2R));
*delta_Trop = delta_w + delta_d;
return true;
}
标签：phi,SPP,double,电离层,C++,unsigned,GNSS,delta,Trop 来源： https://blog.csdn.net/qq_32109917/article/details/113423687