平流层产品及其验证

3 minute read

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的温湿压延伸到地面台站高度。