平流层产品及其验证
Published:
前言
在GNSS技术及其误差中提到GNSS信号在电离层中传播时，受到对流层/电离层/等离子层/磁层的影响，具体如下：
Thayer(1974)给出：
\(N=k_1(P_d/T)Z^{-1}_d+k_2(P_w/T)Z^{-1}_w+k_3(P_w/T^2)/Z^{-1}_w+k_4\frac{n_e}{f^2}\) ，$P_d, P_w$分别是干空气和水汽的分压强；其中$Z^{-1}_d, Z^{-1}_w$分别为干空气和水汽的可压缩系数，
\[Z^{-1}_d=1+P_d[57.97 \times 10^{-8}(1+0.52/T)-9.4611\times 10^{-4}t/T^2]\approx 1\] \[Z^{-1}_w=1+1650(P_w/T^3)(1-0.01317 t+1.75 \times 10^{-5} t^2 +1.44 \times 10^{-6} t^3)\approx 1\]；$T, t$为绝对温度（K）和摄氏温度；或引入状态方程为
\[N=k_1 R_d \rho+k'_2(P_w/T)Z^{-1}_w+k_3(P_w/T^2)/Z^{-1}_w+k_4\frac{n_e}{f^2}\]，其中$k’_2=k_2-k_1(R_d/R_w)=16.48 K\times hPa^{-1}$即
\[N(s)=N_h(s)+N_w(s)+N_e(s)\]令$\Delta L_z, \Delta L_{zh}, \Delta L_{zw}$分别表示天顶总延迟，天顶静力学延迟和天顶湿延迟，则有
\[\Delta L_z=10^{-6}\times\int_{X_j}^{X_k}N(s))ds=\Delta L_{zh}+\Delta L_{zw}+\Delta L_{ze}\]\[\Delta L_{zh}=10^{-6}\times\int_{X_j}^{X_k}N_h(s)ds=10^{-6}\times\int_{X_j}^{X_k} k_1 R_d \rho ds\]
\[\Delta L_{zw}=10^{-6}\times\int_{X_j}^{X_k}N_w(s)ds=10^{-6}\times\int_{X_j}^{X_k} (k'_2\frac{P_w}{T}Z^{-1}_w+k_3\frac{P_w}{T^2}Z^{-1}_w) ds\]
\[\Delta L_{ze}=10^{-6}\times\int_{X_j}^{X_k}N_e(s)ds=10^{-6}\times\int_{X_j}^{X_k} k_4 \frac{n_e}{f^2} ds\]
其中$R_d=287(J\cdot kg^{-1}\cdot K^{-1})$；$R_w=461.5(J\cdot kg^{-1}\cdot K^{-1})$；$\rho=\rho_d+\rho_w$，即干空气及湿水气密度，$k_1=77.604\pm0.014 K \cdot hPa^{-1}$；$k_2=64.79\pm0.08 K \cdot hPa^{-1}$； $k_3=3.776\pm0.014K^2 \cdot hPa^{-1}$,$k_4=4.028\times 10^7$。
饱和水汽压计算公式
在这里引入三个饱和压计算公式
void ptrh2huref_guo(double p, double t, double rh, double*h, double*rd,double*rw,double*re){
double T0 = 273.16, e_sat,e0=6.1121;
// guopeng,doc,
e_sat = e0*pow(10.,7.5*(t-T0)/t);
*h = rh / 100 * e_sat;
//[折射率与温湿压](#ref_eq)
//N = k1(Pd / T)Z^-1_d + k2(Pw / T)Z^-1_w + k3(Pw / T ^ 2)Z^-1_w + k4 ne / f ^ 2, Z^-1_d = Z^-1_w = 1
// k1 = 77.604, k2 = 64.79, k3 = 3.776e5, k4 = 4.03e7
*rd = 77.604*((p-*h) / t);
*rw = 64.79*(*h / t) + 3.776e5*(*h / t / t);
*re=0.0;
}
void ptrh2huref_ucar(double p, double t, double rh, double*h, double*rd,double*rw,double*re){
double T0 = 273.16, Td = 243.5, vaps = 0, vapr = 0, dp, e0 = 611.21;
// ucar,perl
if (t <= 0)
vaps = -999;
else
vaps = (6.112*exp(17.67*(t - T0) / (t - T0 + 243.5)));
if (vaps == -999||rh<=0)dp = -999;
else {
vapr = vaps*rh / 100;
dp = 243.5*(log(6.112) - log(vapr)) / (log(vapr) - log(6.112) - 17.67);
}
if (dp == -999)*h = 0;
else if (dp < -240)*h = -999;
else *h = (6.112*exp(17.67*dp / (dp + 243.5)));
//[折射率与温湿压](#ref_eq)
//N = k1(Pd / T)Z^-1_d + k2(Pw / T)Z^-1_w + k3(Pw / T ^ 2)Z^-1_w + k4 ne / f ^ 2, Z^-1_d = Z^-1_w = 1
// k1 = 77.604, k2 = 64.79, k3 = 3.776e5, k4 = 4.03e7
*rd = 77.604*((p-*h) / t);
*rw = 64.79*(*h / t) + 3.776e5*(*h / t / t);
*re=0.0;
}
void ptrh2huref_ifs(double p, double t, double rh, double*h, double*rd,double*rw,double*re){
double T0 = 273.16, T1 = 250.16, e_sati, e_satw, e_sat, alpha;
// IFSCY46R1, 332:e_sat(T) = a_1 exp(a_3(T - T0) / (T - a4)),
// ice : a1 = 611.21Pa, a3 = 22.587, a4 = -0.7k, T0 = 273.16;
// water:a1 = 611.21Pa, a3 = 17.502, a4 = 32.19k, T0 = 273.16;
e_sati = 6.1121*exp(22.587*(t - T0) / (t + 0.7));
e_satw = 6.1121*exp(17.502*(t - T0) / (t - 32.19));
alpha = t <= T1 ? 0 : t <= T0 ? (t - T0)*(t - T0) / (T1 - T0) / (T1 - T0) : 1.0;
e_sat = alpha*e_satw + (1 - alpha)*e_sati;
*h = rh / 100 * e_sat;
//[折射率与温湿压](#ref_eq)
//N = k1(Pd / T)Z^-1_d + k2(Pw / T)Z^-1_w + k3(Pw / T ^ 2)Z^-1_w + k4 ne / f ^ 2, Z^-1_d = Z^-1_w = 1
// k1 = 77.604, k2 = 64.79, k3 = 3.776e5, k4 = 4.03e7
*rd = 77.604*((p-*h) / t);
*rw = 64.79*(*h / t) + 3.776e5*(*h / t / t);
*re=0.0;
}
ECMWF数据
从ECMWF-ERA5数据中，我们可以获得1 hPa～1000 hPa大气分布数据，具体包括温湿压高度等，我们可以得到湿压与干压。 由于高度是通过地势计算的，所以不一定处于同一海拔高，也可能与地面陆地的高度也不一致。地势与海拔高关系如下：
\[\frac{(r+h)^2}{r^2} d \frac{\phi}{g_1}=d h\]其中，$g_1=(1 - a \cdot cos(2 \phi))*g_0, a=0.00259, g_0= 9.80616$或
\[d \phi=g(\phi,h)d h\]而$g(\phi,h)=g_1(\phi) / (1 + h/r_0)^2,r_0=6371229$。 具体C代码如下：
void z2h(double lat,double*z, double*h, int n){
const double g0 = 9.80616, r0 = 6371229;
double clat = cos(lat*3.1415926 / 180), g1 = (1 - 0.00259*cos(2 * lat*3.1415926 / 180))*g0, g2 = 0;
int i = 0;
h[n - 1] = z[n - 1] * r0 / (g1*(r0 - z[n - 1] / g1));
for (i = n - 2; i >= 0; i--){
g2 = g1 / (1 + h[i + 1] / r0) / (1 + h[i + 1] / r0);
h[i] = h[i + 1] + (z[i] - z[i + 1]) / g2;
}
for (i = 0; i < n; i++)h[i] /= 1e3;
}
我们可以得到一般性的GeoID数据，然后将ECMWF的温湿压延伸到地面台站高度。