2024, Vol. 44
2. 桂林理工大学旅游与风景园林学院，桂林市雁山街319号，541006
大气可降水量(PWV)是衡量大气中水汽含量的重要指标之一[1]。随着无线电探空站(RS)、卫星遥感、GNSS及大气再分析资料等多种水汽产品的不断发展，越来越多的学者倾向于利用多源水汽资料获取更全面的水汽信息[2-3]。鉴于使用不同数据源的水汽数据高程存在差异，亟需建立一个高精度的水汽垂直改正模型。
诸多学者对PWV改正模型的构建进行了一系列研究[4-8]，这些研究模型虽表现出较好的性能，但未考虑PWV精细的时变特性。中国地区地形起伏较大，气候复杂多变，现有模型难以描述PWV剖面信息，对获取不同高度的PWV信息造成困难。针对这一问题，本文使用2016~2018年ERA5 PWV剖面数据，采用高斯函数构建顾及高程及精细时间变化的中国地区PWV垂直剖面模型(C-GPWV模型)，利用2019年ERA5 PWV剖面数据及RS PWV对模型进行精度检验。
1 实验数据研究区域范围为15°~55°N、70°~135°E，ERA5数据由欧洲中期天气预报中心(ECMWF)免费提供，本文使用0.25°×0.25°小时剖面数据，数据可在官方网站(https://cds.climate.copernicus.eu/)免费获取。探空站数据资料由美国俄怀明大学(http://weather.uwyo.edu/upperair/sounding.html)免费提供下载。
2 模型建立 2.1 基于GPT3模型计算PWVGPT3模型是在GPT2w模型基础上发展的最新一代模型，可提供多种对流层参量，GPT3模型计算公式详见文献[9]。GPT3模型计算PWV需要与Askne模型[10]联合，计算公式为：
| $ \left\{\begin{array}{l} \mathrm{PWV}=\varPi \cdot \mathrm{ZWD} \\ \varPi=\frac{10^6}{\rho_w R_V\left[\left(k_3 / T_m+k_2^{\prime}\right)\right]} \end{array}\right. $ | (1) |
式中，k′ 2和k3为大气折射率常数，Tm为大气加权平均温度，ρw为液态水密度，Rv为水汽气体常数。由于GPT3模型提供的气象参数有1°×1°和5°×5°两种分辨率，用GPT3-1表示1°分辨率的GPT3模型，GPT3-5表示5°分辨率的GPT3模型。
2.2 PWV时空特性分析及C-GPWV模型建立文献[11]评估ERA5 PWV格网数据的精度表明，ERA5 PWV具有良好的精度和可靠性。因此本文利用2016~2018年ERA5分层数据，通过积分法[12]计算出0~12 km范围的格网PWV，使用高斯函数表示PWV垂直剖面变化：
| $ \mathrm{PWV}_h=a \cdot \mathrm{e}^{\left[-\frac{(h-b)^2}{c^2}\right]} $ | (2) |
式中，a、b、c代表高斯模型的系数，h为高程，PWVh为地表至高度h处的水汽总量。进一步可由式(2)推导出：
| $ \mathrm{PWV}_{h_1}=\operatorname{PWV}_{h_2} \cdot \mathrm{e}^{\left[\frac{\left(h_2-b\right)^2-\left(h_1-b\right)^2}{c^2}\right]} $ | (3) |
式中，h2为起算高度，PWVh2为起算高度PWV；h1为目标高度，PWVh1为目标高度PWV。
为验证高斯函数对PWV剖面的拟合效果，选取研究区域两个格网点：格网点1(20°N，80°E)和格网点2(30°N，120°E)。通过对2019-01-01 UTC 00:00的ERA5格网分层数据进行积分计算获取PWV，并采用高斯函数对其进行拟合，结果如图 1所示。
|
图 1 ERA5格网PWV剖面拟合结果及偏差 Fig. 1 ERA5 grids PWV profile fitting results and deviations |
由图 1可知，格网点1和格网点2的PWV剖面与高斯函数拟合线较吻合。通过残差图可知，格网点1在高程处于0~2.5 km时表现出正偏差，当高程处于2.5~12 km时表现出负偏差；整个高度区间内，最小偏差为-0.31 mm，最大偏差为0.82 mm。格网点2在高程处于0~5.7 km时偏差较大，当高程处于0~1.7 km时，主要表现为负偏差，当高程处于1.7~5.5 km时，主要表现为正偏差；在整个高度区间，最小偏差为-0.58 mm，最大偏差为0.85 mm。
为研究模型系数b和c在研究区域2016~2018年的季节性变化特征，对研究区内所有格网PWV剖面使用高斯函数进行逐小时拟合，以获取模型系数。选择格网点1(20°N，80°E)进行展示，通过计算得到2016年~2018年日平均模型系数b和c的时间序列，进行年周期和半年周期的三角函数拟合，并采用快速傅里叶变换(FFT)探测分析b和c的周期特性，结果如图 2所示。
|
图 2 模型系数b、c的时间序列变化与频谱分析 Fig. 2 Time series variation and spectral analysis of model coefficient b, c |
由图 2可知，模型系数b呈现显著的季节性变化特征，数值范围在-35~1 km之间，三角函数拟合曲线与时序散点趋势相同；频谱分析结果进一步揭示模型系数b在年周期和半年周期上都表现出明显的变化。模型系数c也表现出明显的季节性变化特征，数值范围在0~15 km之间；频谱分析结果进一步表明，模型系数c也具有显著的年周期特征，半年周期特征较弱。
为了更深入地研究模型系数b和c精细的时间变化特征，对格网点1的年均逐小时模型系数b和c进行顾及日周期和半日周期的三角函数拟合，并采用快速傅里叶变换(FFT)分析其周期性，结果如图 3所示。
|
图 3 模型系数b、c的日变化与频谱分析 Fig. 3 Diurnal variation and spectral analysis of model coefficient b, c |
由图 3可知，模型系数b呈现明显的日变化与半日变化特征，主要表现为从UTC 00:00开始出现降-升-降的变化趋势，通过频谱分析进一步说明模型系数b具有明显的日周期和半日周期特征。模型系数c也具有相似的变化特征，具体表现为UTC 00:00开始，模型系数c呈升-降-升的变化趋势；由频谱分析可知，模型系数c也具有显著的日周期及半日周期特征。
综上所述，模型系数b、c存在显著的季节和日变化特征，因此对模型系数b、c进行顾及年周期、半年周期、日周期和半日周期的三角函数拟合，表达式为：
| $ \begin{gathered} \text { cof }^i= \\ A_0^i+A_1^i \cos \left(2 \pi \frac{\text { doy }}{365.25}\right)+A_2^i \sin \left(2 \pi \frac{\text { doy }}{365.25}\right)+ \\ A_3^i \cos \left(4 \pi \frac{\text { doy }}{365.25}\right)+A_4^i \sin \left(4 \pi \frac{\text { doy }}{365.25}\right)+ \\ A_5^i \cos \left(2 \pi \frac{\text { hod }}{24}\right)+A_6^i \sin \left(2 \pi \frac{\text { hod }}{24}\right)+ \\ A_7^i \cos \left(4 \pi \frac{\mathrm{hod}}{24}\right)+A_8^i \sin \left(4 \pi \frac{\mathrm{hod}}{24}\right) \end{gathered} $ | (4) |
式中，i表示第i个格网点；cofi表示第i个格网点的模型系数b或c；Ai0表示第i个格网的年均值，振幅系数分别用(Ai1, Ai2)、(Ai3, Ai4)、(Ai5, Ai6)和(Ai7, Ai8)表示，对应第i个格网点的年、半年、日和半日周期；doy和hod分别表示年积日和UTC时。
利用2016~2018年ERA5 PWV剖面数据获取模型系数b和c，使用最小二乘法求解式(4)中的模型系数值，结果如表 1所示。联立式(3)和式(4)可构建中国区域PWV垂直剖面模型(C-GPWV模型)。使用该模型时，只需要输入年积日和UTC时，利用表 1和式(4)计算出相应时刻的模型系数b和c，获得待插值点的起算PWV，将模型系数b和c代入式(3)，输入待插值点的起算高程和目标高程，即可得到对应目标高度的PWV值。
|
表 1 利用ERA5 PWV剖面计算的C-GPWV模型的系数 Tab. 1 Coefficient of C-GPWV model calculated using ERA5 PWV profiles |
为验证C-GPWV模型精度，采用平均绝对误差(MAE)与均方根误差(RMSE)作为标准评价3个模型的精度，将残差偏离其平均值超过3倍标准偏差的值剔除。以2019年ERA5 PWV剖面数据与2019年探空站地表数据为模型验证值。MAE与RMSE计算公式为：
| $ \left\{\begin{array}{l} \mathrm{MAE}=\frac{1}{n} \sum\limits_{i=1}^n\left|y_i-y_i\right| \\ \mathrm{RMSE}=\left(\frac{1}{n} \sum\limits_{i=1}^n\left(y_i-\hat{y}_i\right)^2\right)^{\frac{1}{2}} \end{array}\right. $ | (5) |
式中，yi为模型计算值，
为验证C-GPWV模型在中国范围内对PWV垂直插值的精度和适用性，选取2019年积分计算得到的0~12 km内PWV剖面数据作为参考值，对C-GPWV、GPT3-1和GPT3-5模型值进行精度评定，结果如表 2和图 4所示。
|
|
表 2 中国区域ERA5 PWV剖面精度验证结果 Tab. 2 Accuracy verification of ERA5 PWV profile in China |
|
图 4 中国区域ERA5 PWV剖面插值精度平均MAE和RMSE空间分布 Fig. 4 The average MAE and RMSE spatial distribution of ERA5 PWV profile interpolation accuracy in China |
由表 2可知，相比于GPT3-1和GPT3-5模型，C-GPWV模型进行ERA5 PWV剖面插值时展现出较好的精度，其MAE在0.06~2.27 mm内变化，RMSE在0.10~3.24 mm内变化，平均MAE和RMSE为0.85 mm和1.32 mm，比GPT3-1模型的平均MAE和RMSE分别提高55.03%和54.33%，比GPT3-5模型的平均MAE和RMSE分别提高55.50%和55.10%。与GPT3-5模型相比，GPT3-1模型的平均MAE和RMSE均较低，这主要是因为GPT3-5模型的空间分辨率低于GPT3-1模型，导致GPT3-5模型的系统误差大于GPT3-1模型，进而使得GPT3-5模型的精度较低。
由图 4可知，相对于GPT3-1和GPT3-5模型，整体上C-GPWV模型在中国地区表现出较小的MAE，在中国南方地区和西北地区最为明显。从RMSE来看，C-GPWV模型的RMSE明显低于GPT3-1和GPT3-5模型，C-GPWV模型的RMSE在大部分地区都处于0~2 mm之间，在中国的青藏高原地区和北方有较小的RMSE，说明C-GPWV模型在高海拔地区有较好的改正精度。相比于GPT3-1和GPT3-5模型，C-GPWV模型在南方地区的精度改善最为明显。由此可见，C-GPWV模型对PWV垂直改正的性能和稳定性均优于GPT3-1和GPT3-5模型。
为评估C-GPWV模型在不同气压层上的垂直插值精度，本文均匀选取5个气压层，使用3个模型将ERA5地表PWV数据插值到5个气压层，并统计平均MAE和RMSE，结果见表 3。
|
|
表 3 ERA5 PWV剖面数据在不同大气压力层下的精度验证结果 Tab. 3 Accuracy verification results of ERA5 PWV profiles data under different atmospheric pressure layers |
由表 3可知，在不同气压下，C-GPWV模型在5个气压层下的平均MAE和RMSE均低于GPT3-1和GPT3-5模型。相比于GPT3-1和GPT3-5模型，C-GPWV模型在950 hPa气压层的精度提升最明显。在950 hPa气压层时，相比于GPT3-1和GPT3-5模型，C-GPWV模型平均MAE、RMSE分别提升84.34%、84.41%和85.01%、85.02%。由于GPT3模型的建模源是月平均数据，其模型仅包含年周期和半年周期项[13]，而PWV随高度变化影响较大，但GPT3模型没有考虑高程改正，这就导致GPT3-1和GPT3-5模型的误差较大。
3.2 利用RS PWV对模型进行精度验证使用C-GPWV模型将探空站周围最近的4个ERA5格网点的PWV值改正到探空站相同高度，再使用反距离加权法将4个格网点的PWV值插值到探空站点上，并与GPT3-1和GPT3-5模型进行精度对比，结果如表 4和图 5所示。
|
|
表 4 中国区域RS PWV精度验证结果 Tab. 4 Accuracy verification results of RS PWV in China region |
|
图 5 中国区域RS PWV插值精度平均MAE和RMSE空间分布 Fig. 5 The average MAE and RMSE spatial distribution of RS PWV interpolation accuracy in China region |
由表 4可知，GPT3-1和GPT3-5模型平均MAE分别为5.05 mm和5.20 mm，C-GPWV模型的平均MAE为1.99 mm，并且C-GPWV模型的最大MAE均低于GPT3-1和GPT3-5模型的平均值。C-GPWV模型的平均RMSE是3个模型中最低的，为2.60 mm，相比于GPT3-1和GPT3-5模型分别减少3.06 mm和3.21 mm，平均RMSE相比于GPT3-1和GPT3-5模型分别减少3.91 mm和4.05 mm，RMSE相比于GPT3-1提高60.06%，相比于GPT3-5提高60.90%。经过C-GPWV模型改正后的PWV精度得到较好的改善，表现出较小的MAE和RMSE。
由图 5可知，3个模型在南方区域均表现出较大的MAE和RMSE，西北地区和北方地区的MAE和RMSE均较低，这主要是因为南部沿海气候湿润，全年降雨较为充沛，水汽通常表现出较强的季节性变化，导致模型难以捕获复杂的水汽特征，影响PWV估计值。对比分析GPT3-1和GPT3-5模型，整体上两个模型的MAE和RMSE相差不大，但在南方地区部分探空站点GPT3-1模型的MAE和RMSE小于GPT3-5模型。C-GPWV模型在整个研究区内大多数站点的MAE和RMSE都远小于GPT3-1和GPT3-5模型，进一步表明C-GPWV模型在中国地区具有较好的PWV垂直改正能力与良好的稳定性。
由于PWV受纬度和高程影响较大，为验证3个模型在不同纬度和高度范围内的插值精度变化情况，将中国探空站根据高度划分为5个区间，根据纬度划分为4个区间，绘制3个模型各高度段及纬度段的MAE和RMSE分布，结果如图 6和7所示。
|
图 6 3种模型在探空站不同高度范围内的精度统计 Fig. 6 Accuracy statistics of three models in different altitude range of sounding stations |
|
图 7 3种模型在探空站不同纬度范围内的精度统计 Fig. 7 Accuracy statistics of three models in different latitude range of sounding stations |
由图 6可知，随着探空站点高程的增加，3个模型的MAE和RMSE均出现降-升-降的变化趋势，且在5个高程范围内，C-GPWV模型的MAE和RMSE均是最低的，说明C-GPWV模型在不同高度区间均有很好的稳定性及精度。在5个高程范围内，当高程区间处于0~500 m时，C-GPWV模型的垂直插值精度较低；当高程区间处于500~1 000 m时，C-GPWV模型的垂直改正精度最高。
由图 7可知，3个模型估计PWV的精度都受纬度的影响，C-GPWV模型在4个纬度带范围内的MAE和RMSE均最低。GPT3-1和GPT3-5模型在中低纬度地区估计的PWV精度都较差，高纬度地区的精度较高。在4个纬度带中，C-GPWV模型的MAE和RMSE均低于GPT3-1和GPT3-5模型，在高纬度带区域表现出最好的插值精度，主要是因为低纬度地区主要处于南方及沿海区域，水汽情况复杂多变，对模型的垂直改正影响较大。
4 结语针对中国地区气候多变、水汽复杂、高精度PWV垂直剖面模型缺乏的情况，本文构建了顾及精细季节性变化及日周期变化的中国区域PWV垂直改正模型(C-GPWV模型)，该模型提供模型系数，使用方便。在使用ERA5 PWV剖面和探空站地表PWV数据进行模型精度验证时，C-GPWV模型均表现出较好的插值精度。下一步将研究该模型不同水汽数据的插值性能。
| [1] |
Lu C X, Li X X, Ge M R, et al. Estimation and Evaluation of Real-Time Precipitable Water Vapor from GLONASS and GPS[J]. GPS Solutions, 2016, 20(4): 703-713 DOI:10.1007/s10291-015-0479-8
(  0) |
| [2] |
谢劭峰, 王义杰, 黄良珂, 等. 中国区域大气加权平均温度垂直递减率格网模型[J]. 测绘学报, 2023, 52(2): 206-217 (Xie Shaofeng, Wang Yijie, Huang Liangke, et al. A Grid Model for the Lapse Rate of Atmospheric Weighted Mean Temperature over China[J]. Acta Geodaetica et Cartographica Sinica, 2023, 52(2): 206-217)
( 0) |
| [3] |
黄良珂, 莫智翔, 刘立龙, 等. 顾及时变递减因子的中国大陆地区大气可降水量垂直改正模型[J]. 测绘学报, 2021, 50(10): 1 320-1 330 (Huang Liangke, Mo Zhixiang, Liu Lilong, et al. An Empirical Model for the Vertical Correction of Precipitable Water Vapor Considering the Time-Varying Lapse Rate for China's Mainland[J]. Acta Geodaetica et Cartographica Sinica, 2021, 50(10): 1 320-1 330)
( 0) |
| [4] |
Huang L K, Mo Z X, Liu L L, et al. Evaluation of Hourly PWV Products Derived from ERA5 and MERRA-2 over the Tibetan Plateau Using Ground-Based GNSS Observations by Two Enhanced Models[J]. Earth and Space Science, 2021, 8(5)
( 0) |
| [5] |
Huang L K, Liu W, Mo Z X, et al. A New Model for Vertical Adjustment of Precipitable Water Vapor with Consideration of the Time-Varying Lapse Rate[J]. GPS Solutions, 2023, 27(4)
( 0) |
| [6] |
Mo Z X, Zeng Z L, Huang L K, et al. Investigation of Antarctic Precipitable Water Vapor Variability and Trend from 18 Year(2001 to 2018) Data of Four Reanalyses Based on Radiosonde and GNSS Observations[J]. Remote Sensing, 2021, 13(19)
( 0) |
| [7] |
Zhao Q Z, Du Z, Yao W Q, et al. Precipitable Water Vapor Fusion Method Based on Artificial Neural Network[J]. Advances in Space Research, 2022, 70(1): 85-95 DOI:10.1016/j.asr.2022.04.031
( 0) |
| [8] |
Zheng Y X, Lu C X, Wu Z L, et al. Machine Learning-Based Model for Real-Time GNSS Precipitable Water Vapor Sensing[J]. Geophysical Research Letters, 2022, 49(3)
( 0) |
| [9] |
Landskron D, Böhm J. VMF3/GPT3: Refined Discrete and Empirical Troposphere Mapping Functions[J]. Journal of Geodesy, 2018, 92(4): 349-360 DOI:10.1007/s00190-017-1066-2
( 0) |
| [10] |
Askne J, Nordius H. Estimation of Tropospheric Delay for Microwaves from Surface Weather Data[J]. Radio Science, 1987, 22(3): 379-386 DOI:10.1029/RS022i003p00379
( 0) |
| [11] |
Zhang Y L, Cai C S, Chen B Y, et al. Consistency Evaluation of Precipitable Water Vapor Derived from ERA5, ERA-Interim, GNSS, and Radiosondes over China[J]. Radio Science, 2019, 54(7): 561-571
( 0) |
| [12] |
Ssenyunzi R C, Oruru B, D'ujanga F M, et al. Performance of ERA5 Data in Retrieving Precipitable Water Vapour over East African Tropical Region[J]. Advances in Space Research, 2020, 65(8): 1 877-1 893
( 0) |
| [13] |
杨飞, 郭际明, 陈明, 等. GNSS大气加权平均温度经验模型精化方法的建立和分析[J]. 测绘学报, 2022, 51(11): 2 339-2 345 (Yang Fei, Guo Jiming, Chen Ming, et al. Establishment and Analysis of a Refinement Method for the GNSS Empirical Weighted Mean Temperature Model[J]. Acta Geodaetica et Cartographica Sinica, 2022, 51(11): 2 339-2 345)
( 0) |
2. College of Tourism and Landscape Architecture, Guilin University of Technology, 319 Yanshan Street, Guilin 541006, China