2019, Vol. 39
2. 中国科学院大学, 北京市玉泉路19号甲, 100049
在GNSS导航定位中，可利用模型修正法和气象数值预报资料计算对流层延迟值，对于后者广为使用的有欧洲中尺度预报中心和美国国家环境预报中心(NCEP)提供的再分析资料。ECMWF分层数据计算的中国区域ZHD、ZWD的精度优于其他模型[1]，在亚洲地区计算的ZTD的平均bias和RMS分别为-1.0 cm和2.7 cm[2]。目前，利用ECMWF计算ZTD的研究大多采用高度分层数据[1-4]，采用地表数据的较少。本文利用ECMWF再分析地表格网资料计算中国区域ZTD的精度，通过GPT2w模型得到ECMWF格网点的温度递减率和水汽压递减率，并结合格网点与测站间的高差进行高程改正，分别以中国区域探空站数据求得的ZTD和GNSS测站的对流层产品为参考值，对ECMWF地表资料计算的中国区域的ZTD精度进行统计分析，为今后利用ECMWF再分析地表资料计算中国区域ZTD提供参考。
1 数据利用2010~2015年ECMWF地表数据，结合GPT2w模型提供的温度递减率和水汽压递减率参数进行高度归算，计算ZTDECMWF，并以2010~2015年探空数据计算的ZTDsonde和中国大陆构造环境监测网络(陆态网)GNSS测站2015年ZTD数据作为参考值评定其精度。
1.1 ECMWF再分析地表资料选用中国区域2010~2015年ECMWF地表资料，时间分辨率为6 h(UTC 00：00、06:00、12:00、18:00), 水平分辨率为1°×1°，地面数据包括地表气压P0(单位Pa)、2 m露点温度d2m(单位K)、2 m温度T0(单位K)及地面格网点高度H0(单位m)，研究范围为70°~140°E、15°~55°N。
1.2 探空资料探空资料来自全球站点2010~2015年无线电探空资料数据集(IGRA)。探空数据可提供从地面到约30 km高度范围内的大气压、温度、相对湿度等大气参数，时间分辨率为12 h(UTC 11：00，23:00)，考虑到测站在某些观测历元会有数据缺失，为保证研究结果的可靠性，剔除观测历元数少于4 000的测站，确定75个探空站参与本次研究，其空间分布见图 1中蓝色图标。
|
图 1 探空站与陆态网GNSS测站分布 Fig. 1 Distribution of radiosonde and CMONOC GNSS stations |
由于探空分层资料是通过探空气球在不同高度获取的，当其高度太低时，计算得到的ZTDsonde偏差较大[4]，因此只选取探空气球最大高度大于24 km[5]的观测历元数据进行计算。
1.3 陆态网对流层天顶延迟产品使用GNSS测站2015年对流层产品，时间分辨率为1 h，可从网站(ftp://ftp.cgps.ac.cn/products/troposphere/)获取，包括ZTD值及ZTD精度评定值。测站坐标可从网站ftp://ftp.cgps.ac.cn/logs/中的logs文件获取，237个GNSS观测站空间分布见图 1中红色图标。在下载的GNSS测站对流层产品中，剔除标准差大于15.9 mm[6]的数据，提高结果的可靠性。
1.4 GPT2w模型GPT2w模型将全球地表划分为1°× 1°和5°×5°两种格网[7]，分别记作GPT2_1w和GPT2_5w，为与ECMWF数据分辨率(1°×1°)一致，本文主要采用GPT2_1w提供的温度递减率dT及水汽压递减率λ。
2 计算方法 2.1 探空数据计算ZTD对流层天顶方向的总延迟计算式为：
| $ Z \mathrm{TD}=10^{-6} \int_{H}^{\infty} N \cdot \mathrm{d} z $ | (1) |
式中，ZTD为总延迟值(单位m)，H为测站高(单位m)，N为大气折射率。N的计算式为[8]：
| $ N=k_{1} \frac{P_{d}}{T}+k_{2} \frac{e}{T}+k_{3} \frac{e}{T^{2}} $ | (2) |
式中，Pd为大气干气压(单位hPa)；e为水汽压(单位hPa)，T为温度(单位K)。k1、k2、k3为大气折射常数，k1=77.689 K/hPa，k2=71.295 K/hPa，k3=375 463 K2/hPa。对于探空数据分层部分，通过积分法[1]求天顶方向的对流层延迟值；对于探空气球所能达到的最大高度上方的大气层，以最高层的气象参数作为输入参数，利用Saastamoinen模型计算该部分的ZTD，2部分ZTD之和即为测站的ZTDsonde。
2.2 ECMWF再分析地表资料计算ZTD由于ECMWF格网点的平面位置及高程与测站不同，首先对ECMWF格网点的气象参数进行高度归算和水平内插，得到测站的气象参数，再利用Saastamoinen模型计算测站的ZTD。具体步骤为:
1) 读取ECMWF地表数据中的气压P0(为便于后续计算，将其单位换算为hPa)、2 m露点温度d2m、2 m温度T0、格网点高程及数据观测时间。由2 m露点温度求水汽压e0:
| $ {e_0} = 6.1078\exp \left( {17.269\frac{{{d_{2{\rm{m}}}} - 273.15}}{{{d_{2{\rm{m}}}} - 35.85}}} \right) $ | (3) |
2) 读取测站经纬度及高度H，将观测时间转为简化儒略日dmjd，利用GPT2w模型得到格网点温度递减率dT及水汽压递减率λ。
3) 根据测站坐标确定邻近4个ECMWF格网点的经纬度，进而确定格网点的气象参数值(P0，T0，e0)和高程H0。
4) 将测站高度与格网点高程作差得到两者的高差ΔH，结合温度递减率dT及水汽压递减率λ对GPT2w模型气象参数进行高度改正，得到邻近4个格网点在测站同高度位置上的气象参数(Ppoint，Tpoint，epoint)，具体计算公式为：
| $ {T_{{\rm{point }}}} = {T_0} + {\rm{d}}T\Delta H $ | (4) |
| $ {P_{{\rm{point }}}} = {P_0}\exp ( - c\Delta H) $ | (5) |
| $ {e_{{\rm{point }}}} = {e_0}\left( {\frac{{{P_{{\rm{point }}}}}}{{{P_0}}}} \right){\lambda ^{ + 1}} $ | (6) |
| $ c = {g_m}\frac{{{\rm{d}}{M_{{\rm{tr}}}}}}{{{R_g}{T_v}}} $ | (7) |
| $ \begin{array}{l} {g_m} = 9.784\left( {1 - 2.66 \times {{10}^{ - 3}}\cos (2\varphi ) - } \right.\\ \left. {2.8 \times {{10}^{ - 7}}{H_0}} \right) \end{array} $ | (8) |
| $ {T_v} = {T_0}(1 + 0.6077Q) $ | (9) |
| $ Q = 0.622{e_0}\left( {{p_0} - 0.378{e_0}} \right) $ | (10) |
式中，Q为比湿，gm为重力加速度，dMtr为大气摩尔质量，其值为28.965×10-3 kg/mol，Rg为气体常数，其值为8.314 3 J/(mol·K)，Tv为虚温。
5) 利用双线性内插法[9]，由格网点在测站同高度位置的气象参数得到测站气象参数(P，T，e)。
6) 将上述计算得到的测站气象参数(P，T，e)代入Saastamoinen模型，计算得到测站的对流层天顶延迟ZTDECMWF。
2.3 由ECMWF地表资料计算的ZTD精度评价方法在以探空ZTD为参考值时，因探空数据观测时间为11：00及23:00，为减小时间差对计算结果的影响，选取ECMWF数据观测时间为12：00及00:00的观测数据进行计算。同样为保证与ECMWF地表资料时间一致，选取陆态网GNSS测站及ECMWF数据每天00：00、06:00、12：00及18:00的对流层延迟值进行对比分析。将求得的ZTDECMWF与参考值ZTDreference(包括ZTDsonde、ZTDGNSS)作差，再计算差值的平均值与均方根。
3 ECMWF地表资料计算ZTD精度分析为便于区分，以探空数据计算的ZTD为参考值的精度评定结果记为biassonde、RMSsonde，以陆态网GNSS测站的ZTD为参考值的精度评定结果记为biasGNSS、RMSGNSS。
3.1 以探空ZTD为参考值评定由ECMWF地表资料得到的大气参数和ZTD的精度 3.1.1 气象参数(P，T，e)的精度统计结果以中国区域75个探空站2015年实测地表气象参数为参考值，计算通过上述方法得到的探空站气象参数相对于实测值的bias与RMS。表 1为探空站大气压、温度和水汽压的平均bias和RMS，括号内数据为bias与RMS的最小值和最大值。由表 1可知，气压与水汽压的平均RMS分别为1.76 hPa与1.98 hPa，温度的平均RMS为1.96 K。
|
表 1 ECMWF地表数据结合GPT2w模型所得中国区域探空站气象参数的精度统计结果 Tab. 1 Test results of meteorological parameters derived from ECMWF surface data with GPT2w in Chinese radiosonde stations |
各气象参数RMS的空间分布见图 2。由图可知，大多数测站气象参数的计算精度较高，因此ECMWF地表资料结合GPT2w模型得到的气象参数可用于ZTD估计及可降水量计算等。
|
图 2 ECMWF地表资料结合GPT2w模型得到的气象参数的RMS分布 Fig. 2 RMS distribution of meteorological parameters derived from ECMWF surface data and GPT2w |
将§3.1.1得到的大气参数代入Saastamoinen模型计算ZTDECMWF, 将利用探空站数据计算的ZTDsonde作为真值对ZTDECMWF的精度进行验证。为方便比较，同时对由探空站观测的地表数据代入Saastamoinen模型计算得到的ZTD进行精度分析。统计2010~2015年ECMWF地表资料与探空站地面观测资料计算的ZTD的bias和RMS随纬度的变化情况，结果见图 3。
|
图 3 ECMWF地表资料与探空地面观测资料计算ZTD的bias和RMS随纬度分布情况 Fig. 3 Distribution of bias and RMS of ZTD calculated by ECMWF surface and radiosonde surface observed at different latitudes |
由图 3可知，随着纬度的增加，biassonde变化并不明显, RMSsonde呈逐渐减小趋势。同时，有接近1/2测站的由ECMWF地表资料结合GPT2w模型计算的ZTD的RMSsonde小于探空地面观测资料计算的结果，且这些测站大多分布于低纬度地区(15°~30°N)；对于中高纬度地区(30°~55°N)测站，ECMWF地表资料计算ZTD的RMSsonde大多略大于探空数据计算的结果。造成该现象的原因可能是Saastamoinen模型系数在低纬度地区有偏差[10]，当输入参数(ECMWF地表数据结合GPT2w模型得到的气象参数)有偏差时，两种误差部分抵消，反而减小了ZTD的估计误差。
利用2010~2015年ECMWF地表资料与探空地面观测资料计算的ZTD的bias与RMS的平均值见表 2。
|
|
表 2 ECMWF地表资料与探空地面观测资料计算的ZTD bias与RMS的平均值 Tab. 2 Mean values of bias and RMS of ZTD calculated by ECMWF surface and radiosonde surface observed |
从表 2可以看出，ECMWF地表资料计算的ZTD的平均biassonde为0.07 cm，平均RMSsonde为3.72 cm，略小于探空地面观测资料，说明在没有地面实测气象数据时，ECMWF地表资料可得到与探空地面观测资料精度相当的结果。对比文献[10]发现，ECMWF地表资料计算的ZTD的RMS比利用GPT2w模型计算的ZTD的RMS小，精度有所提高，而利用ECMWF地表资料计算的ZTD与陈钦明等[3]的结果基本一致，精度优于无需实测气象参数的经验模型，弱于利用分层资料计算的结果。其原因是ZWD的计算精度与气象参数密切相关，分层资料相比于地面资料在信号传播路径上提供的气象数据更多，将整个信号传播路径上的对流层延迟进行了积分，故计算的ZTD精度更高。
利用2010~2015年ECMWF地表资料与探空站地面观测资料计算的ZTD的bias及RMS随高度的变化情况见图 4。由图可知，随着高度的增加，biassonde变化并不明显，但RMSsonde整体呈减小的趋势。
|
图 4 ECMWF地表资料与探空地面观测资料计算ZTD的bias和RMS随高度分布情况 Fig. 4 Distribution of bias and RMS of ZTD calculated by ECMWF surface and radiosonde surface observed at different altitudes |
将陆态网237个测站2015年的ZTD数据作为参考值，评定利用ECMWF地表资料结合GPT2w模型求得的ZTD精度，研究分析所有测站ZTD的平均bias和RMS的最值与均值，结果见表 3。由表可知，ECMWF地表资料结合GPT2w模型计算的ZTD的平均RMS为3.41 cm, 与文献[11]中RMS最小的GPT2模型相比，平均RMS较小。
|
|
表 3 ECMWF地表资料计算ZTD的平均bias和RMS的最值与均值 Tab. 3 The most and average values of average bias and RMS about ZTD calculated by ECMWF surface |
利用ECMWF地表资料结合GPT2w模型求得的ZTD的bias和RMS随纬度的分布见图 5。由图可知，ECMWF地表资料计算的ZTD的biasGNSS的绝对值随纬度的增加有减小趋势，RMSGNSS随纬度的增加也有递减的趋势。
|
图 5 ECMWF地表资料计算ZTD的bias和RMS随纬度的分布情况 Fig. 5 Distribution of bias and RMS of ZTD calculated by ECMWF surface at different latitudes |
以陆态网GNSS测站的ZTD为参考值，研究利用ECMWF地表资料计算的ZTD的bias和RMS在中国的分布情况，结果见图 6。由图可知，biasGNSS在西部如新疆、西藏等地区有较大的正偏差，在东南地区有较大的负偏差；RMSGNSS在西部地区相对较小，在东南地区较大，其原因可能是相比于西部地区，东南沿海地区水汽含量高、变化较大，计算的ZTD误差较大。
|
图 6 ECMWF地表资料计算ZTD的bias和RMS分布情况 Fig. 6 Distribution of bias and RMS of ZTD calculated by ECMWF surface |
综上所述，利用实测气象参数计算的ZTD和ECMWF地表数据结合GPT2w模型得到的气象参数计算的ZTD精度基本一致，说明在没有地面气象站或气象站数据质量不佳的情况下，可利用ECMWF地表资料的气象参数计算ZTD。
4 结语本文分别以中国地区75个探空站2010~2015年探空数据计算的ZTD及陆态网237个观测站2015年的ZTD产品为参考值，分析评估利用ECMWF再分析地表资料计算的中国地区ZTD的精度，通过统计分析ZTD的bias和RMS的空间分布，论证利用ECMWF再分析地表资料计算ZTD的可行性，得到以下结论：
1) 以探空地面观测气象参数为参考值，ECMWF地表资料结合GPT2w模型计算得到的气象参数的平均bias分别为1.06 hPa、0.26 K、0.94 hPa，平均RMS分别为1.76 hPa、1.96 K、1.98 hPa，计算精度较高，可用于获取测站的气象参数。
2) 以探空数据计算的ZTD为参考值时，ECMWF地表资料计算的ZTD的平均bias和RMS分别为0.07 cm、3.72 cm，与探空地面观测数据计算的结果精度基本相当，甚至在低纬度地区精度更优，说明ECMWF地表数据计算的ZTD已达较高的精度，适用于无地表气象设备情况下中国地区ZTD的估计。随着测站纬度或高度的增加，平均bias的变化趋势并不明显，平均RMS总体呈递减趋势。
3) 以陆态网ZTD为参考值时，ECMWF地表资料计算的ZTD的平均bias和RMS分别为-0.31 cm、3.41 cm，其中RMS比以探空ZTD为参考值的RMS小2 mm。
利用本文提出的方法估算的ZTD的精度在中国各个区域均为cm级，平均精度优于经验对流层模型，可满足GNSS的m级定位精度中对流层延迟改正的需要，也可用作精密定位中对流层的先验估计值。同时，ECMWF地表资料具有数据量小、方便传输等优势，还可提供测站气压、温度等气象参数估计值，方便GNSS气象学等相关应用。本次数据采用的分辨率为1°×1°，接下来将利用其他分辨率数据进行类似研究，从而进一步探讨分辨率对结果的影响。
致谢: 感谢ECMWF、CMONOC、NCDC提供相关的数据资料。
| [1] |
Chen B Y, Liu Z Z. A Comprehensive Evaluation and Analysis of the Performance of Multiple Tropospheric Models in China Region[J]. IEEE Transactions on Geoscience and Remote Sensing, 2016, 54(2): 663-678 DOI:10.1109/TGRS.2015.2456099
(  0) |
| [2] |
Chen Q M, Song S L, Stefan H, et al. Assessment of ZTD Derived from ECMWF/NCEP Data with GPS ZTD over China[J]. GPS Solutions, 2011, 15(4): 415-425 DOI:10.1007/s10291-010-0200-x
( 0) |
| [3] |
陈钦明, 宋淑丽, 朱文耀. 亚洲地区ECMWF/NCEP资料计算ZTD的精度分析[J]. 地球物理学报, 2012, 55(5): 1 541-1 548 (Chen Qinming, Song Shuli, Zhu Wenyao. An Analysis of the Accuracy of Zenith Tropospheric Delay Calculated from ECMWF/NCEP Data over Asian Area[J]. Chinese Journal of Geophysics, 2012, 55(5): 1 541-1 548)
( 0) |
| [4] |
马志泉, 陈钦明, 高德政. 用中国地区ERA-Interim资料计算ZTD和ZWD的精度分析[J]. 大地测量与地球动力学, 2012, 32(2): 100-104 (Ma Zhiquan, Chen Qinming, Gao Dezheng. Study on Accuracy of ZTD and ZWD Calculated from ERA-Interim Data over China[J]. Journal of Geodesy and Geodynamics, 2012, 32(2): 100-104)
( 0) |
| [5] |
Berrada B H, Golé P, Lavergnat J. A Model for the Tropospheric Excess Path Length of Radio Waves from Surface Meteorological Measurements[J]. Radio Science, 1988, 23(6): 1 023-1 038 DOI:10.1029/RS023i006p01023
( 0) |
| [6] |
Hopfield H S. Tropospheric Effect on Electromagnetically Measured Range: Prediction from Surface Weather Data[J]. Radio Science, 1971, 6(3): 357-367
( 0) |
| [7] |
Böhm J, Möller G, Schindelegger M, et al. Development of an Improved Empirical Model for Slant Delays in the Troposphere(GPT2w)[J]. GPS Solutions, 2014, 19(3): 433-441
( 0) |
| [8] |
Smith E K, Weintraub S. The Constants in the Equation for Atmospheric Refractive Index at Radio Frequencies[J]. Proceedings of the IRE, 1953, 41(8): 1 035-1 037 DOI:10.1109/JRPROC.1953.274297
( 0) |
| [9] |
Leandro R F, Santos M C, Langley R B. A North America Wide Area Neutral Atmosphere Model for GNSS Applications[J]. Navigation, 2009, 56(1): 57-71
( 0) |
| [10] |
Zhou C C, Peng B B, Li W, et al. Establishment of a Site-Specific Tropospheric Model Based on Ground Meteorological Parameters over the China Region[J]. Sensors, 2017, 17(8): 1 722-1 739 DOI:10.3390/s17081722
( 0) |
| [11] |
王君刚.GNSS对流层建模与应用[D].上海: 同济大学, 2016 (Wang Jungang. GNSS Tropospheric Delay Modeling and Application[D]. Shanghai: Tongji University, 2016)
( 0) |
2. University of Chinese Academy of Sciences, A19 Yuquan Road, Beijing 100049, China