基于再分析资料ERA5计算对流层延迟精度分析

Accuracy Analysis of Tropospheric Delay Calculation Based on Reanalysis Data ERA5

摘要: 对流层延迟是GNSS导航定位的重要误差源，目前有多种方法可实现天顶对流层延迟(ZTD)的估算。ERA5作为ECMWF发布的最新一代再分析产品，能够提供覆盖全球的高时空分辨率气象资料，据此可估算出高精度的ZTD。然而，当前对ERA5对流层延迟产品精度评价的相关研究仍甚少。本文以高精度的IGS对流层延迟产品为参考，在编程实现ERA5对流层延迟估算的基础上，对其精度进行了评价。试验结果表明：相较于GPT2w-1S模型和GPT3-1S模型，由ERA5气象资料估算的ZTD产品具有更高的精度，均方根误差分别提高了27.53 mm和27.50 mm。

Abstract: Tropospheric delay is an important error source for GNSS navigation and positioning, and there are many methods to estimate the zenith tropospheric delay (ZTD). As the latest generation of reanalysis products released by ECMWF, ERA5 can provide meteorological data with high spatial and temporal resolution covering the whole world, based on which high-precision ZTD can be estimated. However, there are still few related studies on the accuracy evaluation of ERA5 tropospheric delay products. This paper takes the high-precision IGS tropospheric delay product as a reference, and evaluates its accuracy on the basis of programming the ERA5 tropospheric delay estimation. The experimental results show that compared with the GPT2w-1S model and the GPT3-1S model, the ZTD product estimated from the ERA5 meteorological data has higher accuracy, and the root-mean-square deviation is increased by 27.53 mm and 27.50 mm, respectively.

1. 引言

电磁波信号在经过对流层时，会受到大气折射的影响，产生对流层延迟ZTD (Zenith Total Delay) [1]。ZTD作为卫星定位解算过程中的一种重要误差源，需要对其进行改正，从而达到提高GNSS导航定位精度的目的。Hopfield模型和Saastamoinen模型是目前比较常用的对流层延迟改正模型 [2] [3]，但受限于实测气象参数的获取，这两种模型并不适用于多数地区的对流层延迟解算。为摆脱实测气象参数限制，气象再分析资料被应用于对流层延迟改正。当前，欧洲中期天气预报中心ECMWF发布的第五代全球气候再分析资料数据集ERA5在国际上应用较广，与上一代产品ERA-Interim相比，ERA5在时间分辨率上有了较大提升，从6 h提高到了1 h。张永林等 [4] 提出一种基于再分析资料ERA5的天顶对流层延迟计算方法，并与前一代再分析资料ERA-interim的计算结果进行对比分析，结果表明ERA5在精度上较ERA-Interim有明显提升。李瑞辉 [5] 基于气象要素的检核方案，利用ERA5计算气象要素和ZTD的精度，研究表明，利用ERA5提取的气象要素可代替实测值用于对流层研究。

但目前对ERA5对流层延迟产品精度评价的相关研究仍甚少，且多数研究是进行同类模型间的比较。因此，本文在利用ERA5估算对流层延迟ZTD的基础上，与分辨率为 $1˚\times 1˚$ 的GPT2w_1 + Saastamoinen、GPT3_1 + Saastamoinen (以GPT2w-1S、GPT3-1S表示)估算的ZTD进行比较分析，旨在为ERA5资料在对流层延迟改进中的应用提供参考。

2. 对流层延迟的估计方法

2.1. 基于Saastamoinen模型的对流层估计

Saastamoinen模型是基于Shell定律，利用天顶距三角函数展开被积函数进行两层积分所建立 [3]。从地面到对流层10 km为第一层，温度下降率为6.5℃/km；从对流层10 km顶部到70 km处的层顶区间为第二层，温度为一个恒定值 [6]。Saastamoinen模型计算公式如下：

$ZTD=\frac{0.002277}{f\left(\phi ,h\right)}\times \left[P+e\times \left(0.05+\frac{1255}{T}\right)\right]$ (1)

$e=rh6.11\times {10}^{\frac{7.5T-273.15}{T}}$ (2)

$f\left(\phi ,h\right)=1-0.00266\cos 2\phi -0.00028h$ (3)

式中，T为地面温度(K)，P为地面压强(hPa)，e为地面水汽压(hPa)，rh为地面相对湿度。 $f\left(\phi ,h\right)$ 为地球自转引起的重力加速度的修正， $\phi$ 为测站纬度，h为测站高程(km)。本文中，上述气象参数(T、P、e)可由ERA5气象再分析资料提供或利用分辨率为 $1˚\times 1˚$ 的GPT2w/3模型解求。

2.2. ERA5气象资料与Saastamoinen模型协同的对流层延迟估计

ERA5是欧洲中期天气预报中心ECMWF最新发布的一种再分析产品(用户可从https://cds.climate. copernicus.eu/注册下载)，水平分辨率为 $0.25˚\times 0.25˚$，垂直分辨率为37层，相比于ERA-Interim资料，主要在时间分辨率上有了较大提升，从6 h提高到了1 h。由于ERA5数据顶层高度一般为47 km左右，该位置几乎没有湿延迟影响 [7]，且ERA5只能提供气象数据，不能直接提供对流层延迟结果，所以本文将基于再分析资料ERA5计算对流层延迟分两部分进行。对于顶层以下的ERA5数据，利用积分延迟计算方法求取对流层延迟 $ZT{D}_1$，再引入Saastamoinen干延迟模型计算顶层以上对流层延迟 $ZT{D}_2$ (以ERA5-S表示)，最终得到总对流层延迟，即为两部分对流层延迟之和。具体计算方法如下：

$ZT{D}_1={10}^{-6}{\displaystyle \int N\text{d}h}={10}^{-6}{\displaystyle \underset{i}{\overset{n-1}{\sum }}\left(N_i+N_{i+1}\right)}\times \left(h_{i+1}-h_i\right)/2$ (4)

$ZT{D}_2=\frac{0.002277\times P_{\text{top}}}{f\left(\phi ,h_{\text{top}}\right)}$ (5)

$N=\frac{k_1\left(P-e\right)}{T}+\frac{k_{2}e}{T}+\frac{k_{3}e}{T}$ (6)

$e=\frac{qP}{\left(0.622+0.378q\right)}$ (7)

上式中，k₁ = 77.604 K/hPa，k₂ = 64.79 K/hPa，k₃ = 377,600 K²/hPa，N为大气折射指数，n为测站上方ERA5数据包含的气压总层数，P_top为顶层气压值，h_top为顶层高度，其余气象参数含义与式(1)一致。

ERA5数据按气压分层，各等压面之间的高度存在差距，然而测站点很少出现与ERA5数据格网点重合，或与等压面高程相等的情况，因此在使用ERA5数据时还需对其做以下处理：

1) 时间与高程的单位统一。ERA5数据采用格里高利历作为时间单位，在计算时应将其转换为UTC时。GNSS测站采用大地高作为高程系统，而ERA5数据采用位势高来表示高度，在处理测站处气象参数之前需进行高程统一，具体方法参考文献 [8]。

2) 完成高程单位的统一后，根据测站点坐标找到距离站点最近的四个格网点，通过线性外推或内插方法，得到四个格网点在高程等于测站点高程时对应的气象参数 [9]。

3) 对步骤2)中的气象参数进行双线性内插，得到测站位置处的气象参数。

4) 对测站点高程上方各层等压面进行双线性内插，得到平面位置与测站坐标一致、高度等于各层高度时的气象参数。

5) 将处理后的气象参数代入公式(4)至(7)即可得到最终的对流层延迟ZTD。

上述过程均在MATLAB软件平台上完成。

2.3. 基于GPT2w/3模型的对流层延迟估计

Böhm等在GPT2模型的基础上增加了大气加权平均温度和水汽压垂直梯度两个参数，建立了GPT2w (Global Pressure and Temperature 2 wet)模型 [10]；Landskron等通过改进映射函数VMF3，提出一种新的基于全球网格的经验模型GPT3 [11]。GPT2w/3模型均可提供两种空间分辨率的地表温度、地表气压、水汽压和大气加权平均温度等气象参数，为目前公开的标称精度较高的对流层延迟经验模型。该类模型计算各气象参数的公式为：

$\begin{array}{c}\Phi =A_0+A_1\cos \left(\frac{doy}{365.25}2\pi \right)+B_1\sin \left(\frac{doy}{365.25}2\pi \right)\\ +A_2\cos \left(\frac{doy}{365.25}4\pi \right)+B_2\sin \left(\frac{doy}{365.25}4\pi \right)\end{array}$ (8)

上式中，doy为年积日， $A_0$ 为气象参数的平均值， $A_1$ 和 $B_1$ 为气象参数的年周期振幅， $A_2$ 和 $B_2$ 为气象参数的半年周期振幅。表1为上述2种对流层延迟模型参数类型对比。

Table 1. Comparison of GPT2w/3 model parameters

表1. GPT2w/3模型参数对比

表1中，(Lat, Lon, H)分别为测站经纬度和大地高，dmjd为简化儒略日，year为年份，doy为年积日，ah为干映射函数，aw为湿映射函数，la为水汽下降因子，undu为大地水准面差距，Gn_h为流体静力北梯度，Ge_h为流体静力东梯度，Gn_w为湿北梯度，Ge_w为湿东梯度。

通过测站坐标找到距离最近的四个格网点，对其输出的气象参数进行高程改正，求出测站高程处的格网点气象参数：

$\{\begin{array}{l}T=T^{\prime }+dT\times dh\\ P=P^{\prime }\times \exp \left[\frac{-g_m\times dM}{R_g\times T\times \left(1+0.6077q\right)}\times dh\right]\\ e^{\prime }=\frac{qp}{\left(0.622+0.378q\right)}\\ e=e^{\prime }\times {\left(\frac{100P}{P^{\prime }}\right)}^{\lambda +1}\end{array}$ (9)

该式中，( $T^{\prime }$, $P^{\prime }$, $e^{\prime }$ )分别为格网点气温、气压、水汽压；(T, P, e)分别为高程改正后测站点高度的气温、气压、水汽压；dT为温度下降率；dh为站点与格网点之间的高差； $g_m$ 为重力加速度，一般取值为9.80665 m/s²；dM为气体相关常数，取值为28.965 × 10⁻³ kg/mol；Rg为通用气体常数，取值为8.3143 J/(K∙mol)；q为比湿； $\lambda$ 为水汽下降因子。将高程改正后的四个格网点气象参数进行双线性插值可得到测站位置的气象参数，结合Saastamoinen模型，即可输出ZTD。

3. 数据源与精度评价指标

3.1. 数据源

本文选用中国地区5个同属于IGS中心的GNSS测站(URUM、LHAZ、JFNG、HKSL、TWTF) 2021年全年的ERA5格网数据(测站分布如图1所示)，水平格网分辨率为 $0.25˚\times 0.25˚$ 、垂直方向上包含37层等压面、时间分辨率为1 h。

Figure 1. GNSS station distribution

图1. GNSS测站分布

3.2. 精度评价指标

IGS/ZTD产品精度高，不确定性小，是当前国际上公认的精度最高的实测数据 [12]，因此本文以IGS/ZTD日均值作为真值，用于检验模型的精度。采用偏差(bias)和均方根误差(RMSE)作为精度评价指标，计算公式为：

$\{\begin{array}{l}\text{bias}={\displaystyle \underset{i}{\overset{n}{\sum }}\left(ZT{D}_i-ZT{D^{\prime }}_i\right)/n}\\ \text{RMSE}=\sqrt{{\displaystyle \underset{i}{\overset{n}{\sum }}{\left(ZT{D}_i-ZT{D^{\prime }}_i\right)}^2/n}}\end{array}$ (10)

上式中，n为样本个数，此处n = 365； $ZT{D}_i$ 为IGS/ZTD真值； $ZT{D^{\prime }}_i$ 为模型所得ZTD。

4. 实验结果分析

利用前文介绍的计算方法分别解算得到ERA5-S、GPT2w-1S及GPT3-1S 3种延迟模型的ZTD值，并以高精度的IGS对流层延迟产品为参考，其结果如图2所示。

由图2可看出，对流层延迟整体上呈现出先上升后下降的变化趋势，波峰部分主要集中在年积日200左右，该时段正处于夏季，符合降雨较多的实际情况。3种模型中，ERA5-S模型的对流层延迟与IGS/ZTD真值变化趋势一致，整体上反映了对流层延迟随时间的变化特征，GPT2w-1S与GPT3-1S模型的对流层延迟曲线几近重合，表明2种模型的精度大致相当，且该2种对流层延迟模型所得ZTD变化较小，与IGS/ZTD真值存在较大差异。

Figure 2. ZTD values of three tropospheric delay models at each station

图2. 各测站3种对流层延迟模型的ZTD值

由图3可看出，不同测站所得偏差bias与RMSE存在差异，总体而言，大部分测站ERA5-S模型的偏差bias与均方根误差RMSE均小于其它两种模型。为更好地分析各模型解算ZTD的精度，给出各测站bias和RMSE计算结果，如表2所示。

Figure 3. Different model deviations and root mean square errors at each station

图3. 各测站不同模型偏差与均方根误差

Table 2. Three model accuracy evaluation indicators

表2. 3种模型精度评价指标

从表2可看出，ERA5-S模型的偏差bias最大值为5.33 mm，平均值为−11.11 mm。相比于GPT2w-1S模型bias最大值51.37 mm、平均值23.55 mm；GPT3-1S模型bias最大值51.38、平均值23.57 mm，均有较大提升。同时对比3种模型的RMSE可知，ERA5-S模型的RMSE均值为23.08 mm，相较于GPT2w-1S模型RMSE均值50.61 mm、GPT3-1S模型RMSE均值50.58 mm，其RMSE分别提高了27.53 mm及27.50 mm，表明利用ERA5-S模型计算得到的对流层延迟ZTD精度在3种模型中精度最高。

5. 结语

针对当前ERA5对流层延迟产品精度评价的相关研究较少的情况。本文基于2021年全年的高时空分辨率的ERA5气象再分析资料，利用分层积分法结合Saastamoinen模型组成ERA5-S模型，解算对流层延迟ZTD，通过与上述2种模型ZTD进行精度对比，验证了ERA5-S模型具有更高精度，表明再分析资料ERA5可较好应用于对流层延迟改正。同时扩大研究区域，延长时间范围，分析ERA5获取ZTD周期性变化以及高程差异的影响，深入探究ERA5资料改进对流层延迟的效果将在后续工作中开展。

基金项目

国家自然科学基金项目(41901225、52064005)；贵州大学测绘科学与技术研究生创新实践基地建设项目(贵大研CXJD[2014]002)。

NOTES

^*第一作者。

^#通讯作者。