基于Emardson模型的GNSS-PWV转换系数精化

李黎， 侯晓玲， 张雯雯， 高颖， 赵伟， 陈国栋

1 苏州科技大学地理科学与测绘工程学院， 苏州 215009 2 苏州科技大学北斗导航与环境感知研究中心， 苏州 215009

0 引言

水汽在大气中所占比例很小，仅占大气体积的0.1%～0.3%，但它却是龙卷风、雷暴雨和台风等灾害性天气形成的重要因素之一(宋淑丽等, 2004; 毕研盟等, 2006)，获取高时空分辨率的水汽资料可以提高降水预报的精度，减少灾害性天气给国民带来的经济损失.由于传统水汽观测方法的时空分辨率较低，限制了降水预报精度，实践中使用并不广泛(Baker et al., 2001).自Bevis等(1992)提出GNSS气象学以后，GNSS技术以其全球覆盖、全天候、精度高、垂直分辨率高(≤1 km)和费用低等优点逐渐成为获取大气可降水量 (Precipitable Water Vapour, PWV) 的重要手段之一.由GNSS技术(李黎等, 2012)解算得到的对流层延迟(Zenith Total Delay, ZTD)减去根据Saastamoinen模型计算的对流层干延迟(Zenith Hydrostatic Delay, ZHD)可得对流层湿延迟(Zenith Wet Delay, ZWD)，ZWD再乘以水汽转换系数(K值)可得到GNSS-PWV，可见K值是影响GNSS-PWV反演精度的主要因素之一，因此研究并建立高精度的K值模型非常必要.

近40年来，GNSS气象学的发展开拓了高时空分辨率GNSS-PWV的应用研究(Bevis et al., 1994; Wen et al., 2012).而K值的精度在一定程度上受限于大气加权平均温度(Tm)的精度，不少学者已深入研究Tm并建立区域Tm模型.计算Tm需用到地表温度、气压和高程等信息，这些气象要素的实测精度和实时程度会进一步影响Tm的最终精度.Emardson和Derks(2000)基于欧洲数十个探空站资料，建立一个与位置和年积日相关的K值计算模型(Emardson模型).该模型未引入Tm参数，从而避免气象要素的误差积累，进一步简化了计算过程，在实时反演GNSS-PWV中有着广泛的应用前景.鉴于Emardson模型在不同地区适用性不同，姚朝龙等(2015)建立了一个适用于低纬度地区的K值模型，发现其精度主要取决于纬度和高程，与单站模型或者通过Tm计算的K值精度相当；刘立龙等(2016)选取新疆地区9个无线电探空站连续5年的探空资料，根据积分法获取的K值分析Emardson模型在该地区的时空适用性；刘立龙等(2017)分析了地理位置对于K值的影响，在广西地区建立的本地化K值模型显示顾及高程的Emardson模型精度更好；黄良珂等(2019)利用不同的K值模型反演拉萨站的GNSS-PWV，其结果与探空PWV具有较好的一致性，表明顾及高程的Emardson-H模型适用于缺乏气象参数的高海拔地区；侯晓玲等(2021)根据长三角地区7个探空站2年的探空数据，建立了顾及与不顾及高程的两个精化模型，发现在低海拔地区两种模型得到的GNSS-PWV精度相当.

鉴于Emardson模型具有地区差异性，中国地区缺乏统一的高精度K值预报模型，本文利用2018年中国地区的42个探空站资料，建立顾及与不顾及高程的两套Emardson精化模型，并以2019年K值作为参考，分析以上两个模型的K值预报精度，评估其在中国地区的适用性，以期快速获取GNSS-PWV，提高其在气象预报领域的实时应用效率.

1 数据处理方法与数据来源及分析

1.1 K值计算方法

利用探空数据所提供的气压、水汽压和温度等气象参数和积分法是目前应用最广泛的K值计算方法.

(1)

(2)

式中，Pv为平均水汽压(hPa)，T为绝对温度(K)，Δh为大气的厚度(m).T，Δh均可在探空资料中得到，Tm是指探空站上空水汽压与绝对温度沿天顶方向的积分值，是计算K值的关键参数.大气水汽一般分布在离地面12 km以下，故可以通过数值积分法来逼近Tm真值.在数值积分法中，观测误差与逼近误差的综合影响小于1 K(刘焱雄等, 2000)，式(2)是目前公认精度最高的Tm计算方法(Duan et al., 1996; Ingold et al., 1998).Pv无法直接观测，但能通过露点温度和饱和水汽压公式间接计算,

Pv=6.112exp[17.62td/(243.12+td)],

(3)

式中，td为露点温度(℃).

1.2 PWV计算方法

根据GNSS技术解算，可以求得GNSS-ZTD(精度优于6 mm)，一般时间分辨率为5 min.根据地面纬度、气压、高程等信息，基于Saastamoinen模型可精确计算ZHD(Saastamoinen, 1972)：

(4)

式中，Ps为表面压力(hPa)，φ为纬度，H为大地高度(km).测站位置确定后，ZHD是地表压力的函数，压力每增加1 hPa只会引起大约0.2 mm的ZHD误差(Tregoning and Herring, 2006).

通过GNSS-ZTD减去ZHD，可得到ZWD，ZWD与K相乘可得到GNSS-PWV.

(5)

1.3 数据来源及其分析

1.3.1 数据来源与测站信息

本文的探空资料来源于2018—2019年全国均匀分布的42个无线电探空站每天2次的实测数据，采样间隔为12 h，垂直分辨率为1 km.实测数据可从美国怀俄明大学天气数据网(http:∥weather.uwyo.edu/upperair/sounding.html)获取，包括气压、露点温度、高度、地面温度、相对湿度等气象要素，该机构以世界气象组织制定的规则为依据，已对资料进行质量检查.GNSS数据可从中国地震局GNSS产品服务平台(http:∥www.cgps.ac.cn/index.htm)获取(Zhao et al., 2015; Yu et al., 2019)，其中包括GNSS-PWV及其梯度、气压、温度等参数，时间分辨率为1 h.

本文将用2018年中国地区的42个探空站积分法所得的K值建立两套Emardson模型，之后用2019年的42个探空站和11个GNSS站分别验证两套预报模型的K值精度及其PWV反演精度.图1是中国地区42个探空站和11个GNSS站的分布图，表1和表2分别是探空站和GNSS站的地理位置信息(包括站名、站号、代码、经纬度和高程等).

1.3.2K值与测站位置的相关性分析

我国地处太平洋西岸，横跨多个时区，本文以式(1)计算的2018年42个探空站K值作为真值，分析研究K值与测站经度和纬度之间的相关性，结果见图2.由图2可知，K值与测站经度之间的相关系数为0.592，二者的相关性较小；而K值与测站纬度之间的相关系数为0.735，变化趋势呈线性负相关，即纬度越高，测站K值则越小.因此，在建立K值精化模型时，测站经度的影响可以不予考虑，但需要把纬度作为建模的影响因子之一.

中国地区的地形总体呈现西高东低，且起伏较大，因此有必要分析K值与测站高程之间的线性关系.由图3可知，K值与测站高程呈负相关，K值随着测站高程的增加而逐渐下降，两者之间的相关系数为0.941，尤其高程在0～1000 m区间时，下降趋势更为明显.因此，在中国地区建立K值精化模型时，尤其需要考虑测站高程对K值模型精度的影响，要将高程作为建模的主要影响因子之一.

图1 中国地区探空站及GNSS站地理位置分布图Fig.1 Distribution map of selected radiosondes and GNSS stations in China

表1 中国地区42个探空站地理位置Table 1 Geographical position of 42 selected radiosonde stations in China

表2 中国地区11个GNSS站地理位置Table 2 Geographical position of 11 selected GNSS stations in China

图2 K值与探空站经度和纬度之间的相关性分析Fig.2 Correlation analysis between K value and longitude and latitude of the radiosonde stations

图3 K值与探空站高程的相关性分析Fig.3 Correlation analysis between K value and altitude of radiosonde stations

2 中国地区Emardson模型的建立及其精度分析

2.1 Emardson模型的建立

根据年积日、测站纬度和高程等时空特性与K值之间的相关性，分别建立基于测站纬度和年积日的Emardson-I模型和基于测站纬度、高程和年积日的Emardson-H模型(Emardson and Derks, 2000).

(6)

+α4H,

(7)

式中，φ为大地纬度(°)，doy为年积日，H为高程(m)，α0、α1、α2、α3、α4为模型系数，α0为常数，α1为φ的系数，α4为H的系数.

本文将2018年中国地区42个探空数据计算的K值作为真值，使用多元线性拟合法解算得到模型系数(如表3).

表3 中国地区的Emardson-I和Emardson-H模型系数Table 3 The model coefficients of Emardson-I and Emardson-H in China

结合表3、式(6)和式(7)可知，α1代表K值与纬度之间的关系，负值表示K值随探空站纬度的升高而减小，纬度每升高1°，K值约减小0.00023；α4代表K值与高程之间的关系，负值表示K值随探空站处高程的增加而减小，高程每增加1 m，K值则减小0.0000012.因此，在地形平缓的低海拔地区，K值随高程变化的影响较小；而在地形起伏较大时，尤其是我国西部3000 m以上的高海拔地区，高程对K值的影响不可忽视.

2.2 K值预报精度分析

利用2.1节得到的Emardson-I和Emardson-H模型预报2019年42个探空站处的K值，以2019年积分法K值作为参考真值，评估两套K值模型的预报精度.限于篇幅，图4给出了19个探空站基于两套模型的预报K值与积分法K值(当作真值)的时序对比.

由图4可知，在低海拔地区(<600 m，上海、海口、北京、章丘、百色、怀化、汉中)和高海拔地区(>1500 m，西昌、昆明、西宁、甘孜、拉萨、玉树、那曲)，两套模型的预报K值相差较大，最大达到0.004.而在中海拔地区(600～1500 m，伊宁、太原、乌鲁木齐、呼和浩特、贵阳)，两套模型的预报K值偏差很小，尤其是太原、贵阳和呼和浩特等站，均仅相差0.0001，可见在我国中海拔地区，高程对Emardson精化模型影响很小.

图4 Emardson-I和Emardson-H模型的2019年预报K值时序变化Fig.4 Time series variation of predicted K value based on Emardson-I and Emardson-H models in 2019

就两套K值预报模型的时序变化而言，积分法K值与Emardson-H模型的一致性更好，并且海拔越高，Emardson-H模型预报的K值精度也会越高.就中海拔地区而言，Emardson-H模型与Emardson-I模型预报K值与积分法K值相差很小.

表4统计了高中低不同海拔19个探空站的K值预报精度.从表中可知，中国地区Emardson-H模型的适用性相对较好，所得结论与图4类似.低海拔地区的Emardson-H模型的K值预报精度较Emardson-I模型略有提高，其预报K值的平均绝对误差(MAE)均值和均方根(RMS)均值分别为0.00144和0.00173，均优于Emardson-I预报模型的0.00183和0.00218；在中海拔地区，Emardson-H模型相对于Emardson-I模型的K值预报精度提升很小，RMS均值从0.00208提高到0.00202；而在高海拔地区，Emardson-H模型相对于Emardson-I模型的预报精度得到了显著的提升，RMS均值从0.00312提高到0.00181.从表4中还可以看出，用Emardson-I模型预报的K值RMS变化较大，在0.001～0.006之间，海拔越高，偏差越大；而利用Emardson-H模型预报的K值RMS集中在0.002左右，稳定性比前者有大幅度改善.

2.3 GNSS-PWV反演精度分析

2.3.1 GNSS-PWV时序变化及其偏差

由表1和表2可知，安庆、银川和拉萨等处的探空站和GNSS站距离较近，高程相差最多38 m，可认为探空站与GNSS站所受到的对流层延迟非常接近，可用探空PWV作为真值来检验基于Emardson-I和Emardson-H模型的GNSS-PWV精度和可靠性.

限于篇幅，本文仅选择6个不同海拔的探空站和GNSS并址站，分析利用Emardson-I和Emardson-H模型的预报K值反演的2019年GNSS-PWV和探空PWV时序变化.如图5所示，两组GNSS-PWV与探空PWV在变化趋势上表现一致，PWV波峰出现在夏季，PWV波谷出现在冬季，峰值和谷值基本吻合，符合中国夏季多雨、冬季干燥的气候特点，说明两套模型反演的GNSS-PWV具有较高的一致性.

由两套模型计算的GNSS-PWV与探空PWV之间的偏差图可看出，基于Emardson-H模型的GNSS-PWV偏差均小于Emardson-I模型，且夏季的GNSS-PWV偏差较大.在高海拔地区(昌都、拉萨)，GNSS-PWV偏差范围相对较小(-0.5～0.5 mm)；在中海拔地区(银川、若羌)，基于Emardson-I和Emardson-H模型的GNSS-PWV偏差范围相对较大(-1.0～0.5 mm)；而在低海拔地区(安庆、青岛)，PWV偏差范围相对更大(-2.0～0.5 mm)，但基于Emardson-H模型的GNSS-PWV更接近探空PWV，比Emardson-I模型有更为明显的改善.

表4 2019年基于Emardson-I和Emardson-H模型的K值预报精度Table 4 Prediction accuracy of K value based on Emardson-I and Emardson-H models in 2019

2.3.2 GNSS-PWV精度

由表5可以得到，低海拔地区(青岛、安庆)基于Emardson-I模型的GNSS-PWV的MAE均值为0.629 mm，RMS均值为0.752 mm，而基于Emardson-H模型的GNSS-PWV的MAE和RMS均值分别为0.433 mm和0.530 mm，RMS降低了0.222 mm；中海拔地区(若羌、银川、贵阳)基于Emardson-I模型的GNSS-PWV的MAE均值为0.274 mm，RMS均值为0.342 mm，而基于Emardson-H模型的GNSS-PWV的MAE和RMS均值分别为0.266 mm和0.333 mm，RMS仅降低了0.009 mm；而高海拔地区(和田、都兰、昌都、甘孜、拉萨)基于Emardson-I模型的GNSS-PWV的MAE均值为0.159 mm，RMS均值为0.213 mm，而基于Emardson-H模型的GNSS-PWV的MAE和RMS均值分别为0.119 mm和0.154 mm，两种模型的RMS相差在0.1 mm以内.

综上所述，在中国地区，基于Emardson-H模型的GNSS-PWV反演精度明显高于Emardson-I模型，其MAE和RMS均值均优于Emardson-I模型.在中高海拔地区，基于Emardson-H模型的GNSS-PWV相对于Emardson-I模型，精度有所提高，但提高有限；相比而言，Emardson-H模型在低海拔地区的改正效果更为明显，这是因为低海拔地区的ZTD和PWV值均较大，高精度K值更加有利于提高GNSS-PWV精度，说明在低海拔地区考虑高程因子同样重要.

3 讨论与结论

水汽作为重要的气象要素之一，是灾害性天气形成的重要因子，在降水监测和天气预报中扮演着关键角色.中国幅员辽阔、地形起伏大、探空站时空分辨率低，导致大部分区域的水汽资料比较匮乏.地基GNSS技术能探测高时空分辨率PWV，而K值是从GNSS-ZTD转换到GNSS-PWV的重要参数.因此，本文基于Emardson模型建立了两套适合中国地区的高精度K值模型，重要结论如下：

(续图5)

图5 2019年基于Emardson-I和Emardson-H模型的GNSS-PWV时序变化及其偏差Fig.5 Time series variations and their deviations of GNSS-PWV based on Emardson-I and Emardson-H models in 2019

表5 2019年基于Emardson-I和Emardson-H模型的GNSS-PWV精度统计Table 5 The accuracy statistics of GNSS-PWV based on Emardson-I and Emardson-H models in 2019

(1)K值建模需要考虑纬度和高程这两个关键参数.K值与测站纬度和高程之间呈负相关，其相关系数分别为0.735和0.941，而与测站经度的相关性较小.

(2)Emardson-H模型的K值预报精度优于Emardson-I模型，在高海拔地区尤为明显.Emardson-H模型预报的K值MAE和RMS均值分别为0.00150和0.00182，均优于Emardson-I模型的0.00221和0.00255.

(3)基于Emardson-H模型的GNSS-PWV精度也优于Emardson-I模型，尤其在低海拔和高海拔地区.基于Emardson-H模型的GNSS-PWV的MAE和RMS均值分别为0.226 mm和0.283 mm，均优于基于Emardson-I模型的0.288 mm和0.360 mm.

综上所述，就K值而言，Emardson-H模型预报的K值精度优于Emardson-I模型，尤其在高海拔地区；而就GNSS-PWV来说，不仅在高海拔地区，并且在低海拔地区，Emardson-H模型反演的GNSS-PWV精度均优于Emardson-I模型.这是因为低海拔地区的ZTD值较大，高精度K值更有利于提高PWV反演精度，说明本文所建立的全国K值Emardson模型适用于大部分地区.但分析图4可以发现，某些探空站(比如海口、伊宁、乌鲁木齐、西昌等)在春冬季节的模型预报K值偏差较大，意味着K值模型不仅仅要考虑年积日、纬度和高程，可能也要考虑当地的季节变化和气候特征；此外，本文仅考虑了影响因子的线性特征，而没有考虑各参数的非线性变化特点，后续我们还可基于以上两点进一步完善该模型.

致谢感谢美国怀俄明大学天气数据网和中国地震局GNSS产品服务平台为本文提供了探空站资料和GNSS数据.