基于CRA40产品的对流层延迟计算及对比分析

周要宗，楼益栋，张卫星，梁 宏，施 闯，吴 迪，曹云昌

1. 武汉大学卫星导航定位技术研究中心，湖北 武汉 430079； 2. 中国气象局气象探测中心，北京 100081； 3. 北京航空航天大学卫星导航与移动通信融合技术工业和信息化部重点实验室，北京 100191； 4. 南宁师范大学地理与海洋研究院，广西 南宁 530001

对流层延迟是全球卫星导航系统(GNSS)、甚长基线干涉测量(very long baseline interferometry，VLBI)、卫星激光测距(satellite laser ranging，SLR)、卫星测高等空间大地测量技术实现高精度参数估计的主要误差源之一，在数据处理中可利用外部数据产品估算对流层延迟并予以改正[1-3]。全球大气再分析产品，如欧洲中尺度天气预报中心(European Centre for Medium-range Weather Forecasts，ECMWF)再分析中间产品(ECMWF reanalysis interim，ERA-Interim，简称“ERAI”)和第5代产品(the fifth generation ECMWF reanalysis，ERA5)、美国国家环境预报中心再分析产品(The United States National Centers for environmental prediction reanalysis，NCEP)，通过同化多源历史气象观测资料，具有精度高、时空覆盖完整的优点[4]，被广泛用于对流层天顶延迟(zenith path delay，ZPD)和斜路径延迟(slant path delay，SPD)的计算及建模[5-8]，发布了一系列ZPD经验模型(如SHAO系列模型[9]、IGGtrop系列模型[10]、TropGrid系列模型[11]、ITG模型[12]、GTrop模型[13]和HGPT模型[14])和映射函数模型(如UNB-VMF1[15]、VMF3[16]和GFZ-VMF1[17])。

2013年11月，中国气象局启动了我国的全球大气再分析计划，总体目标是建成我国第1代全球大气再分析业务系统，并建成跨度40 a(1979—2018)的中国第1代全球大气/陆面再分析产品(China 's first generation 40 a global atmosphere and land reanalysis，CRA40)[18-19]，质量超过国际第2代，在中国区域接近或达到国际第3代大气再分析水平(http:∥data.cma.cn/analysis/cra40)。2020年12月CRA40产品正式上线并业务化发布[20]。探讨利用CRA40进行对流层延迟计算的方法，并评估CRA40用于对流层延迟计算的精度，尤其是在中国区域的精度，具有重要的意义。文献[21]针对CRA40在中国区域GNSS水汽反演中的应用，对2016年全年CRA40和ERA5的气温、气压、大气水汽加权平均温度和天顶总延迟(zenith total delay，ZTD)进行了较为系统的比较评估，结果表明以GNSS ZTD为参考，CRA40 ZTD精度为1.35 cm，总体较ERA5差。但文献[21]主要讨论的是中国区域的天顶延迟结果，并未开展全球范围内包括斜路径延迟在内的对流层延迟的系统评估。

本文主要对基于CRA40计算的对流层延迟开展系统的评估工作，首先讨论基于CRA40采用射线追踪技术计算对流层延迟的方法，然后分别在231个国际GNSS服务(international GNSS service，IGS)站和213个中国大陆构造环境监测网络(crustal movement observation network of China，CMONOC)站处对ERAI、ERA5和CRA40对流层延迟计算精度(ZPD和5°高度角SPD)进行系统比较评估，并对中国区域CRA40对流层湿延迟计算精度的变化特征进行分析。

1 数据与方法

1.1 CRA40产品及预处理

CRA40(1979—2018)是由中国气象局国家气象信息中心牵头，联合多家单位开发的，整体经历了两个阶段，即准备阶段(2014—2016年)和实施阶段(2016—2019年)[20]。准备阶段主要完成历史气象观测资料的准备和模式同化系统的建设[22-23]，实施阶段主要完成CRA-Interim(2006年7月—2016年12月)产品的制作和评估及CRA40产品的研制[24-26]。CRA40产品是分10段进行研制的，每段4.5 a，段与段之间有半年重叠，以保证产品的连续性。2019年8月全部40 a CRA40再分析产品研制完成，并于2020年12月正式上线业务化[20]。本文采用CRA40气压层产品进行对流层延迟计算，产品模式同化系统为GFS/GSI-3Dvar，时间分辨率为6 h。产品水平格网采用高斯经纬度投影，经向为均匀网格划分，分辨率为0.312 5°，纬向为非均匀网格划分，分辨率在0.312 5°左右。选用参数包括位势高、温度和比湿，位势高和温度参数气压层层数为47，层顶气压为1 hPa(42 km左右)，比湿参数气压层层数为37，层顶气压为100 hPa(16 km左右)。

相关研究表明：同一再分析产品随着水平分辨率的提高，射线追踪的计算效率会显著降低，而对射线追踪ZTD计算精度的影响并不显著[7]，因此综合考虑射线追踪计算效率和精度，本文采用NCEP官网提供的wgrib2工具利用双线性插值方法将再分析资料的水平分辨率统一重采样为1°，来开展后续的再分析资料评估试验。此外，CRA40比湿参数仅有37层，与位势高和温度参数气压层层数不一致，因此需要对比湿参数层顶100 hPa以上气压层(1、2、3、5、7、10、20、30、50、70 hPa)的比湿参数进行处理，得到气压层层数为47的比湿参数。可选的处理方法主要有3种：方法1认为100 hPa(16 km)以上水汽已经相当稀疏，因此，可以简单地将100 hPa以上气压层的比湿参数值设为0 g·kg-1；方法2由现有层顶100 hPa和相邻气压层(如125、150、175 hPa等)的比湿参数来线性外推层顶以上的比湿参数；方法3直接采用100 hPa的比湿参数值作为其他层顶气压层的数值。

为评估以上方法的差异，图1(a)给出了2018年7月19日UTC 15时BZRG测站基于ERA5计算的天顶湿延迟(zenith wet delay，ZWD)随高度的变化曲线，图1(b)给出了2018年7月ERA5月平均产品全球比湿参数均值随气压(37层)的变化曲线。由图1可知，16 km(100 hPa)以上ZWD数值在亚毫米量级，100 hPa以上比湿参数并不为0，几乎处于一个不变值，因此方法1可行但不严密，而方法3相对更为合理。相比之下，方法2利用层顶100 hPa和相邻气压层比湿来线性外推顶层比湿，可能会得到小于0 g·kg-1的比湿。综合考虑计算的可行性和合理性，本文采用方法3处理顶层比湿参数。经预处理后，CRA40产品时间和水平分辨率分别为6 h和1.0°，气压层层数为47，层顶气压为1 hPa(42 km)。

图1 ZWD随高度和比湿随气压的变化曲线Fig.1 Changing curves of ZWD with height and specific humidity with pressure

1.2 ERA5和ERAI产品

ERA5(1950年至今)是ECMWF发布的最新一代全球大气再分析产品，产品模式同化方法为Ensemble of 4D-Var，时间延迟在5 d左右，时间分辨率为1 h，水平分辨率为31 km，气压层产品垂直层数为37[27-28]。ERAI(1979—2019年)是ERA5的上一代全球大气再分析产品，产品模式同化方法为4D-Var，时间分辨率为6 h，水平分辨率为80 km，气压层产品垂直层数为37[29]。本文选用ERA5和ERAI气压层产品进行对流层延迟计算，参数包括位势、温度和比湿，产品时间和水平分辨率与CRA40重采样产品保持一致，分别为6 h和1.0°，层顶气压为1 hPa。

1.3 基于射线追踪技术的对流层延迟计算方法

组合使用再分析产品(0～42 km)和美国标准大气模型(COESA 1976)(42～84 km)来实现参数的中性大气层全覆盖[30]。依据文献[31]给出的方法，采用EGM2008模型[32]将再分析位势高/位势转换成大地高。利用文献[33]给出的方法将再分析比湿转换为水汽压。采用文献[34]给出的大气垂直分层方法(0～2 km：10 m间隔；2～6 km：20 m间隔；6～16 km：50 m间隔；16～36 km：100 m间隔；36～84 km：500 m间隔)来进行垂直大气细化分层。分别使用线性、指数和指数插值方法将温度、气压和水汽压由等压层内插到细化高度层[30]。以温度、气压和水汽压为输入，通过式(1)来计算静力学折射率(nh)、湿折射率(nw)和总折射率(n)[35]

(1)

(2)

分别以6 h和1.0°时空分辨率的CRA40(47层)、ERA5(37层)和ERAI(37层)再分析温度、比湿和位势高/位势为输入，采用前述的对流层延迟射线追踪方法，计算2018年全年全球231个IGS站处(图2(a))和国内213个CMONOC站处(图2(b))的ZPD和5°高度角、16个方位角(0°～337.5°，间隔22.5°)下的SPD。取16个方位角的SPD均值作为5°高度角的SPD。最终得到3种再分析产品在所有选定测站处2018年时间分辨率为6 h的ZPD和5°高度角SPD产品。

图2 IGS和CMONOC站点分布Fig.2 Distribution of selected IGS and CMONOC stations

2 对流层天顶总延迟精度评估

事后GNSS ZTD未被CRA40、ERAI和ERA5同化[26-29]，因此可以将其作为独立参考，来评估基于再分析资料计算的ZTD精度。从IGS官方FTP(gdc.cddis.eosdis.nasa.gov)下载2018年全年选定231个IGS站处时间分辨率为5 min的ZTD产品，产品统计精度可达4 mm[39](https:∥www.igs.org/products#about)。采用PANDA软件和文献[40]所示的数据处理策略对所有选定213个CMONOC测站的GNSS观测数据进行事后精密单点定位(precise point positioning，PPP)处理[41]，得到选定CMONOC测站2018年全年时间分辨率为1 h的ZTD产品，产品统计精度优于5 mm[40]。

以GNSS ZTD为参考，统计2018年全年每个选定测站基于再分资料计算的ZTD误差偏差(Bias)和均方根(root mean square，RMS)。图3给出了全球范围(IGS站网)和中国区域(CMONOC站网)基于3种再分析产品的ZTD Bias地理分布和直方图，图3(a)、(b)和(c)给出全球范围CRA40、ERA5和ERAI的ZTD Bias分布，图3(d)、(e)和(f)给出相应的统计直方图，图3(g)、(h)和(i)给出中国区域CRA40、ERA5和ERAI的ZTD Bias分布，而图3(j)、(k)和(l)则给出相应的统计直方图。在全球范围，ERA5 ZTD Bias平均最小，均值为0.20 cm(图3(b)、(e))。ERAI ZTD Bias所有测站均值为0.34 cm，ZTD偏大的测站主要分布在欧洲、中亚、北美等地区(图3(c)、(f))。相比之下，CRA40 ZTD Bias最大，均值为0.48 cm，主要出现在欧洲、北美及赤道区域(图3(a)、(d))。

图3 ZTD差值Bias分布和直方图Fig.3 Distribution and histogram of ZTD Bias

在中国区域，ERA5 ZTD Bias仍然最小，均值为0.02 cm(图3(h)、(k))。ERAI ZTD Bias次之，均值为0.11 cm(图3(i)、(l))。CRA40 ZTD Bias同样最大，平均为0.28 cm，川滇地区、新疆北部的ZTD正偏差较为明显，负偏差则主要出现在新疆南部区域(图3(g)、(j))。再分析ZTD整体偏差均为正，表明再分析ZTD与GNSS ZTD间可能存在系统性偏差(再分析ZTD数值偏大)，文献[42—43]也给出了类似的结果。

除了平均偏差外，图4也相应地给出了基于3种再分析资料计算的ZTD RMS地理分布和直方图，图4(a)、(b)和(c)分别给出全球范围CRA40、ERA5和ERAI的ZTD RMS分布，图4(d)、(e)和(f)分别给出相应的统计直方图，图4(g)、(h)和(i)分别给出中国区域CRA40、ERA5和ERAI的ZTD RMS分布，而图4(j)、(k)和(l)则给出相应的统计直方图。在全球范围，再分析ZTD RMS与测站纬度相关，赤道及其附近区域RMS明显偏大，这主要是受赤道及其附近区域含量丰富且变化剧烈水汽的影响[44]。ERA5 ZTD RMS整体最小，平均为1.16 cm(图4(b)、(e))。CRA40 ZTD RMS次之，平均为1.39 cm，其中赤道及其附近区域相比ERA5 RMS明显更大(图4(a)、(d))。ERAI ZTD RMS整体最大，平均为1.47 cm，RMS较大的测站主要分布在赤道及其附近区域以及中东、中亚等地区(图4(c)、(f))。

在中国区域，ERA5 ZTD RMS仍然最好，均值为1.16 cm(图4(h)、(k))。ERAI ZTD RMS次之，川滇地区以及新疆、河北、山东、辽宁等地均相对ERA5明显变差，均值为1.39 cm(图4(i)、(l))。CRA40 ZTD RMS相对ERAI略差，主要集中在川滇地区和新疆，均值为1.41 cm(图4(g)、(j))。整体统计上，ERA5 ZTD RMS明显小于CRA40和ERAI，而CRA40和ERAI ZTD RMS比较接近，表明ERA5 ZTD精度显著优于ERAI和CRA40，而ERAI和CRA40 ZTD精度相近。

图4 ZTD差值RMS分布和直方图Fig.4 Distribution and histogram of ZTD RMS

3 对流层斜路径延迟精度评估

由于斜路径延迟缺乏可靠的参考，无法直接对基于再分析资料计算的SPD绝对精度进行评估，而对流层天顶总延迟精度评估结果显示ERA5 ZTD精度最高，理论上ERA5计算的ZPD和SPD精度也最高，因此就以ERA5计算的ZPD和SPD为参考，来评估CRA40和ERAI计算的ZPD和SPD精度。首先计算2018年全年所有选定测站CRA40和ERAI ZPD和SPD同ERA5 ZPD和SPD的差异Bias和RMS，然后统计出差异Bias和RMS的最小值(min)、最大值(max)和均值(mean)，结果见表1。

由表1可知，5°高度角SPD差异Bias和RMS约是相应ZPD差异RMS的10倍，这是因为5°高度角SPD和ZPD比值是在10.2左右[45]。CRA40和ERA5对流层延迟差异主要体现在湿延迟部分。全球范围(IGS站网)CRA40和ERA5对流层延迟(ZPD和SPD)差异RMS分别为1.04(ZTD)、0.13(ZHD)、1.04(ZWD)、10.83(STD)、1.37(SHD)和10.80 cm(SWD)，中国区域(CMONOC站网)为1.18、0.14、1.16、12.30、1.49和12.13 cm。以ERA5对流层延迟为参考，在全球范围CRA40同ERAI的对流层延迟计算精度基本相当，而在中国区域ERAI对流层延迟计算精度略优于CRA40，可能原因是CRA40在中国区域同化了一些国内特有的气象观测数据，导致在中国区域CRA40与ERA5在同化数据源上有更大差异。

表1 再分析ZPD和SPD差异Bias和RMS统计Tab.1 Statistical bias and RMS between reanalysis-based ZPD and SPD cm

4 中国区域CRA40对流层湿延迟精度变化特征

在中国区域CRA40对流层总延迟(ZTD和STD)精度主要取决于湿延迟部分(ZWD和SWD)计算精度，因此重点对中国区域CRA40湿延迟精度的变化特征进行分析。从中国区域选取4个代表性IGS测站，即BJFS(北京房山，温带季风气候)、JFNG(武汉九峰，亚热带季风气候)、LHAZ(拉萨，高原山地气候)和URUM(乌鲁木齐，温带大陆性气候)，来分析2018年CRA40 ZWD和SWD精度变化特征，结果如图5所示。图5中，横轴表示年积日，纵轴(左，蓝色)表示ZWD/SWD差值，蓝色曲线表示CRA40 ZWD/SWD同ERA5 ZWD/SWD的差值时序，纵轴(右，红色)表示ZWD/SWD，红色曲线表示CRA40 ZWD/SWD时序。

由图5可知，同一站点的ZWD精度变化相位同SWD一致、变化振幅约是SWD的1/10，而不同站点间ZWD和SWD精度变化相位和振幅存在明显差异，主要与站点所处气候类型相关。同属季风气候类型的BJFS和JFNG站，ZWD和SWD精度均有明显的季节变化特征：夏季精度最差，冬季精度最好，且温带BJFS站的ZWD和SWD精度变化振幅明显小于亚热带JFNG站，尤其是在冬季(图5(a)、(b)、(c)和(d))。高原山地气候的LHAZ站，ZWD和SWD精度全年较为稳定，基本没有季节变化(图5(e)、(f))。温带大陆性气候的URUM站，ZWD和SWD精度也存在一定的季节变化特征，但变化振幅相对BJFS和JFNG站要小得多(图5(g)、(h))。

图5 2018年BJFS、JFNG、LHAZ和URUM站CRA40 ZWD和SWD差值时序Fig.5 Time series of CRA40 ZWD and SWD differences at BJFS, JFNG, LHAZ and URUM stations in 2018

实际上，CRA40 ZWD和SWD精度直接与对流层水汽的含量及变化幅度相关，对流层水汽含量越丰富、变化幅度越剧烈，ZWD和SWD精度就越差，反之亦然。中国区域季风气候区对流层水汽含量及变化幅度夏季大于冬季，因此季风气候区ZWD和SWD精度夏季差于冬季，高原山地气候区对流层水汽含量较低、变化幅度较小，因此高原山地气候区ZWD和SWD精度全年较为稳定。

5 结论与展望

本文探讨了基于CRA40采用射线追踪技术计算对流层延迟的方法，在全球范围和中国区域对CRA40、ERAI和ERA5对流层延迟(ZPD和5°高度角SPD)的精度进行了系统比较评估，并对中国区域CRA40 ZWD和SWD精度变化规律进行了分析。结果表明：

(1) 以GNSS ZTD为参考，CRA40 ZTD统计Bias在全球范围和中国区域分别为0.48和0.28 cm；统计RMS约为1.40 cm，在全球范围略优于ERAI，在中国区域基本与ERAI相当，整体上较ERA5差。

(2) 以ERA5 STD为参考，CRA40 STD统计Bias在全球范围和中国区域分别为2.84和2.58 cm；统计RMS在全球范围和中国区域分别为10.83和12.30 cm，同ERAI无明显差异。

(3) 中国区域CRA40 ZWD和SWD精度变化规律与气候类型相关，非季风气候区精度全年较为稳定，季风气候区冬季精度明显优于夏季。

结果证实国产自主可控CRA40产品对流层延迟射线追踪精度基本与第3代全球大气再分析ERAI相当，可以作为可靠的数据源来进行对流层延迟射线追踪、建模等相关研究。本文仅对基于CRA40计算的对流层延迟的精度进行了评估，下一步将利用CRA40对流层延迟进行映射函数、水平梯度建模及北斗/GNSS定位应用研究。

致谢：感谢ECMWF提供的ERAI和ERA5数据及中国气象局提供的CRA40数据；感谢IGS提供的IGS ZTD产品和中国地震局提供的GNSS观测数据。