GNSS观测资料来自CMONOC并使用GAMIT/GLOBK（10.4版本）软件进行处理^[16]。中国地壳运动观测网络建于1997—2000年，为GNSS气象学等领域提供高精度的观测资料。在此基础上CMONOC也于2011年投入观测，包含大量卫星激光测距（satellite laser ranging，SLR）、甚长基线干涉测量（very long baseline interferometry，VLBI）、重力及GNSS测站，其中GNSS测站中包含260个基准站以及约2 000个非连续观测站^[17]。本文中使用的ZTD数据为1 h分辨率，由于时间序列的残缺，最终选择黄土高原地区33个GNSS站点处的2014—2017年数据进行使用。

ERA-Interim与ERA5数据分别为ECMWF中的第4代与第5代数据集，包含1979年至今2~3个月的数据集，原始数据可通过https://www.ecmwf.int/免费下载，ECMWF的这两类产品均提供了不同时空分辨率的气压、气温与PWV等格网气象数据，其中，ERA-Interim数据集时间分辨率为6 h^[18]，ERA5数据集时间分辨率为1 h。

本文选择无线电探空数据来验证ERA5中的P和T与GNSS PWV的精度。无线电探空仪数据来自集成的全球无线电探空仪存档版本2（IGRA2）数据集。IGRA2是由美国国家气候数据中心于2016年生成，它比之前的数据集（IGRA1）包含更多的无线电探空站。该数据集包括气压、气温、相对湿度和其他参数，时间分辨率为每天2次或4次（ftp：//ftp.ncdc.noaa.gov/pub/data/igra）。本实验选取黄土高原地区的8个探空站点，黄土高原地区GNSS站与无线电探空（radio sonde，RS）站分布如图 1所示。

图  1  黄土高原地区GNSS和RS站点分布图

Figure 1.  Distribution of GNSS and RS Stations in Loess Plateau

ERA5的地形高度为位势，首先，参照Wang等^[19]和Zhao等^[20]的方法将ERA5的位势转换为位势高；其次，将位势高转换为正高；然后，对ERA5格网的气压与气温数据利用双线性插值获取GNSS站点的气压、气温数据，并对其进行高程改正^[21]：

式中，P为站点的气压数据；$P_{\mathrm{0}}$为格网点的气压数据；h为站点高；$h_{\mathrm{0}}$为格网高；dT为站点温度与格网温度的差值；dh为站点高与格网高的差值。

利用高程改正后的GNSS站点处的气压，依据Saastamoinen模型计算出GNSS站点的天顶静力延迟（zenith hydrostatic delay，ZHD）^[22]，利用ZTD与ZHD可获得天顶湿延迟（zenith wet delay，ZWD）。在获得ZWD后，进行PWV的计算：

式中，$P$为地表气压；$\phi$为纬度；$H$为大地高；$\rho$为水汽的密度；$K_{\mathrm{3}}$=$(\mathrm{3.776}\pm \mathrm{0.014})\times {\mathrm{10}}^{\mathrm{5}}$ ${\mathrm{K}}^{\mathrm{2}}/\mathrm{h}\mathrm{P}{\mathrm{a}}^{}$；$K_{\mathrm{2}}^{\mathrm{\text{'}}}$=16.48 $\mathrm{K}/\mathrm{h}\mathrm{P}{\mathrm{a}}^{}$；$R_v$=461 $\mathrm{J}/(\mathrm{k}\mathrm{g}\cdot \mathrm{K})$；T_m是大气加权平均温度，本文依据Huang等^[23]改进的中国T_m模型IGPT2w（improved global pressure and temperature 2 wet）获得。

为了验证基于GAMIT/GLOBK（version 10.4）软件得到的GNSS ZTD与其他软件处理的一致性，本文还与基于PANDA软件采用精密单点定位技术计算得到的ZTD^[24-25]进行了对比。PANDA软件自2001年起由武汉大学卫星导航定位技术研究中心自主开发，并被国际多家著名研究机构采用，其目标是实现对GNSS/SLR/DORIS/VLBI等各类观测数据的后处理和实时处理分析^[26]。本文中基于PANDA软件得到的ZTD数据采用Zhang等^[27]依据PANDA软件计算得到的中国区域内的ZTD产品，其通过与IGS的ZTD对比发现，RMS为3.90 mm，最小值和最大值分别为2.70 mm和6.40 mm。本文选择黄土高原地区33个GNSS站点处的两种软件获得的ZTD进行对比，将ZTD时间序列中大于3倍中误差的数据进行剔除。图 2给出了使用GAMIT/GLOBK和PANDA软件计算的2014—2017年GNSS的ZTD差值的RMS和Bias分布。从图 2可以看出，两种软件所获得的GNSS的ZTD相差较小，据统计33个站点处的平均RMS和Bias分别为4.05 mm和-0.46 mm。

图  2  2014-2017年两类软件解算的GNSS ZTD差值的RMS和Bias分布

Figure 2.  RMS and Bias Distribution of GNSS ZTD Difference Calculated by Two Types of Softwares from 2014 to 2017

为分析ERA5数据集中气压与气温的精度，本文对比了ERA5与ERA-Interim格网间的气压与气温，由于两种数据集均属于再分析数据资料，因此本节主要讨论二者之间的差异。图 3给出了ERA5与ERA-Interim格网间气压与气温的标准差（standard deviation，STD）与Bias。从图 3可以看出，两种数据集的气压与气温的STD相对较小，气压的Bias较大，说明两种数据集中的气压相差较大，出现这一现象的可能原因是由于该地区高程差异及两类数据集的格网高不同所导致。据统计，两种数据集的气压差值的STD为0.51 hPa，气温差值的STD为1.29 K；气压差值的Bias为-2.12 hPa，气温差值的Bias为0.08 K。由于ERA5的时间分辨率高于ERA-Interim，因此后续使用ERA5进行研究。

图  3  ERA5与ERA-Interim格网间气压与气温的STD与Bias分布图

Figure 3.  STD and Bias Distribution of Pressure and Temperature Between ERA5 and ERA-Interim

依据3倍中误差法对数据进行粗差剔除并保证站点的数据利用率大于85%，选择黄土高原地区8个RS站点处的气压与气温数据对基于ERA5获得的气压、气温数据进行精度评定。

8个站点所对应的站点经纬度及高程信息见表 1。图 4给出了黄土高原地区8个站点处气压与气温的RMS和Bias分布图。从图 4可以看出，气压与气温的对比结果整体较好，气压的结果差于气温，导致这一现象的主要原因是高程变化对气压的影响高于气温。在部分站点处气压与气温的RMS值和Bias值均较大，出现这一现象的主要原因是这些站点对应的原始格网数据质量相对较差。经统计，8个站点处的气压的平均RMS为3.36 hPa，气温的平均RMS为1.97 K；气压的平均Bias为-0.01 hPa，气温的平均Bias为0.04 K。由此可知，基于ERA5的气压、气温数据的精度相对较好，可应用于1 ‍h分辨率PWV的计算。

图  4  ERA5与RS站点间气压与气温的RMS和Bias分布图

Figure 4.  RMS and Bias Distribution of Pressure and Temperature Between ERA5 and RS Stations

由PWV的计算公式可知PWV精度受到多个因子的影响，通过误差传播定律，依据各项影响因子的精度值，可推算出PWV的理论精度值^[14]。下面分别给出ZHD、转换系数Q、PWV的不确定度的计算方式。

1）ZHD数据的不确定度。ZHD的精度主要受气压和常数c的不确定性影响，1 hPa的地表气压约对ZHD产生0.20 mm的误差影响^[28]，ZHD的不确定度计算如下^[14]：

式中，P₀为站点处的气压值（单位：hPa）；${\sigma }_{P_{\mathrm{0}}}$为气压的误差；${\sigma }_c=\mathrm{0.002}\mathrm{4}$，常数c的不确定度参照Zhao等^[15]研究中所使用的数值；λ和H分别是站点纬度和大地高。

2）转换系数Q的不确定度。转换系数Q与Q的不确定度的计算公式如下：

式中，${\rho }_w$为液态水的密度；$k_{\mathrm{2}}^{\mathrm{\text{'}}}$=16.48 $\mathrm{K}/\mathrm{h}\mathrm{P}\mathrm{a}$；$k_{\mathrm{3}}$=$(\mathrm{3.776}\pm \mathrm{0.014})\times {\mathrm{10}}^{\mathrm{5}}$ $\mathrm{K}²/\mathrm{h}\mathrm{P}\mathrm{a}$；$R_w$=461 $\mathrm{J}/(\mathrm{k}\mathrm{g}\cdot \mathrm{K})$；$T_m$是大气加权平均温度；${\sigma }_Q$、${\sigma }_{k_{\mathrm{3}}}$、${\sigma }_{k_{\mathrm{2}}^{\mathrm{\text{'}}}}$和${\sigma }_{T_m}$分别为转化系数Q、常数$k_{\mathrm{3}}$、$k_{\mathrm{2}}^{\mathrm{\text{'}}}$和大气加权平均温度$T_m$的误差。其中，液态水的密度${\rho }_w$的不确定度^[29]为0.002 kg/m³，$R_w$的不确定度为0.008 $\mathrm{J}/(\mathrm{k}\mathrm{g}\cdot \mathrm{K})$（ http://physics.nist.gov/cuu/Constants/index.html），由于来自这两个参数的影响较小（小于Q总不确定度的0.1%），因此Q的不确定度取决于$T_m$、$k_{\mathrm{3}}$和$k_{\mathrm{2}}^{\mathrm{\text{'}}}$的不确定度^[14]。

3）PWV的不确定度。根据误差传播定律，PWV的不确定度^[14，28]计算如下：

式中，${\sigma }_v$、${\sigma }_{\mathrm{Z}\mathrm{T}\mathrm{D}}$、${\sigma }_{P_{\mathrm{0}}}$、${\sigma }_Q$分别代表PWV、ZTD、气压、转换系数Q的误差；V代表PWV；${\sigma }_c=\mathrm{0.002}\mathrm{4}$。

依据上述各项因子的不确定度计算式可获得各项因子对PWV所引起的误差影响。通过计算可知，转换因子Q的误差约为2.26，P的误差约为3.36 hPa，ZHD的误差约为2.79 mm，ZTD的误差约为4.05 mm。其中，ZTD的误差由Zhang等^[27]利用PANDA软件获得的ZTD进行对比得到，ZHD与Q的误差分别由式（3）和式（4）计算获得。

由上述计算结果可知，ZTD对PWV的理论误差影响最大，其次是ZHD。而造成ZHD出现误差的主要原因是气压的误差相对较大，造成Q出现误差的主要原因是大气加权平均温度存在误差。将上述各因子的误差值代入式（5）可得，黄土高原地区PWV的理论精度值为1.51 mm。

为确定黄土高原地区GNSS PWV数据集的精度，依据RS站的PWV数据对其进行精度评定。对基于ERA5气象参数与CMONOC ZTD获得的2014—2017年的GNSS PWV与RS PWV利用3倍中误差对数据进行粗差剔除，并保证每个站点的数据利用率均在85%以上。通过对数据进行粗差处理后，选择黄土高原地区的5个GNSS与RS并址站进行对比。图 5给出了2014—2017年GNSS PWV与RS PWV的RMS和Bias分布图。从图 5可以看出，除XNIN站以外，其余站点处的RMS值均小于3 mm，出现这一现象的主要原因是该站点位于青藏地区，高程差异较大，气压与气温的数据质量较差，因此导致该站点的PWV精度相对较低。据统计，图 5中5个站点处的GNSS PWV与RS PWV的平均RMS和Bias分别为1.94 mm和0.72 mm，由此可知，PWV的实际精度低于理论精度，主要原因是由于ERA5为再分析资料与实际值有一定差距，并且RS PWV本身也存在着一定误差，因此最终导致PWV的实际精度值低于理论精度值。为了清晰地描述GNSS PWV与RS PWV的长时序变化，图 6给出了GSPL站与SNYA站2014—2017年GNSS PWV与RS PWV的长时序变化。

图  5  2014—2017年GNSS PWV与RS PWV的RMS与Bias分布图

Figure 5.  RMS and Bias Distribution of GNSS PWV and RS PWV from 2014 to 2017

图  6  2014—2017年期间测站GSPL与SNYA位置处的RS PWV与GNSS PWV长时序图

Figure 6.  RS PWV and GNSS PWV Long Time Series of Stations GSPL and SNYA from 2014 to 2017

由图 6可知，GNSS PWV和RS PWV的变化趋势一致，故基于ERA5气象参数与CMONOC ZTD所获得的GNSS PWV具有较高精度。

为进一步确定GNSS PWV的精度，利用狄洛尼三角网法将2016年GNSS站点的PWV数据插值为整个黄土高原地区的面状PWV，与ERA5的格网PWV数据进行对比。图 7给出了2016年GNSS PWV与ERA5 PWV及二者差值的分布图。分别按照春、夏、秋、冬4个季节进行分析，图 7（a）为4个季节的ERA5 PWV分布，图 7（b）为4个季节的GNSS PWV分布，图 7（c）为4个季节ERA5 PWV与GNSS PWV的差值分布。

图  7  2016年GNSS PWV与ERA5 PWV对比分布图

Figure 7.  Comparative Distribution of GNSS PWV and ERA5 PWV in 2016

从图 7可以看出，PWV的分布在春季到冬季之间也有很大差异。比较GNSS PWV和ERA5 PWV可以看出，其分布总体上是一致的，在每个季节ERA5和GNSS的PWV数据都有一定的差值，但是在各个季节大多数地区的差值均在5 ‍mm以下，因此基于ERA5气象参数和CMONOC ZTD获得的黄土高原地区的GNSS PWV数据具有较好的精度。此外，在个别区域可以发现ERA5 PWV与GNSS PWV的差值较大，这可能是由于GNSS站不均匀所导致。

表 2给出了GNSS PWV与ERA5 PWV及二者差值分别在4个季节的均值，同样可知二者相差较小，进而说明GNSS PWV的精度相对较高。

大气水汽的快速变化对区域水汽循环和气候变化具有重要意义，因此，对于PWV的影响因素的分析及精度评估变得尤为重要。本文重点研究黄土高原地区GNSS PWV的理论精度与实际精度，分析了影响PWV精度的各因子的不确定度。理论误差的计算主要依据误差传播定理进行研究，实际误差的计算是以RS数据资料为标准来验证GNSS PWV的精度。

通过本文的研究最终得到如下结论：（1）ERA5和ERA-Interim两类数据集的气压在个别区域相差较大，气温整体相差较小。通过利用RS数据集中的气压与气温对ERA5气压与气温的评定发现，气压的平均RMS为3.36 hPa，气温的平均RMS为1.97 K；气压的平均Bias为-0.01 hPa，气温的Bias为0.04 K。（2）影响PWV理论精度的最大影响因素为ZTD，其次为ZHD与转换系数Q，通过PWV的不确定度公式计算得到黄土高原地区PWV的理论精度为1.51 mm，低于PWV的实际精度1.94 mm，出现该现象的主要原因是ERA5为再分析资料，与真实值存在一定差距，RS数据资料本身也存在一定误差。（3）分析2016年期间4个季节的GNSS PWV和ERA5 PWV发现，二者整体相差较小，在个别区域相差较大，主要原因可能是由于GNSS站点不均匀所导致。

基于ERA5气象参数和CMONOC ZTD获得的GNSS PWV精度较好，且具有较高的时间分辨率，这对于研究黄土高原区域水汽变化与分布具有重要意义。