高分辨率的水汽场在天气预报、数值模拟、人工影响天气等领域都具有重要的应用价值，对水汽总量、变化和分布情况的了解程度直接影响数值天气预报中降水预报的准确性。随着全球定位系统技术的发展，利用GPS探测水汽成为当前研究的热点。目前，地基GPS经过一系列的处理可以估计出高精度的天顶湿延迟(zenith wet delay，ZWD)，ZWD和可降水量(precipitable water vapor，PWV)存在着如下的转换关系：

式中，Π为ZWD与PWV的转换系数；ZWD一般为几厘米至几十厘米。R_v=461.495 J/(kg·K)，为水汽特定的气体常数；ρ_w为液态水的密度；k′ ₂=22.1 K/hPa；k₃=3.739×10⁵K²/hPa；T_m为加权平均温度，是利用式(2)计算转换系数Π的关键参数。

由此可见，为了利用GPS估计的ZWD计算得到PWV，转换系数Π非常关键。在实际应用中，如果精度要求不是很高，可近似地取Π为常数0.15^{[1, 2]}，更复杂的则通过建立加权平均温度模型^{[3, 4]}并利用式(2)计算Π值。如果能对转换系数Π直接建模^{[5, 6, 7]}(通常是对Π^-1=10^-6ρ_wR_v\5[(k₃/T_m)+k′ ₂]进行建模)，就可以省去T_m的计算过程，而直接通过ZWD得到PWV。目前的模型都是局部模型，其适用性受到地域的限制，因而建立一个适用于全球的经验模型是很有必要的。

在缺乏地面气温观测资料的条件下,由于GPS反演水汽的公式与气象观测资料的关系是非常清楚、密切的，只是由于缺乏探空资料且积分计算不易，应用时才通过统计方法构建了T_m的Bevis经验公式，把应用探空气象资料的难度降低为应用地面气象资料，而本文又进一步把依赖地面气象资料的难度降到只需地理及天文资料，虽然拓展了模型的易用性，但牺牲其精度以及与实时气象信息的联系也是必然的。

欧洲中期天气预报中心(European Centre for Medium-Range Weather Forecasts,ECMWF)的再分析资料^{[8, 9]}能够以不高于0.75°×0.75°的分辨率提供每天0:00、6:00、12:00和18:00 UTC的PWV格网数据，GGOS Atmosphere 能够以2°×2.5°的分辨率提供每天0:00、6:00、 12:00和18:00 UTC的ZWD格网数据，这些数据是利用ECMWF的ERA-40再分析资料计算得到的。研究表明，ECMWF再分析资料具有较高的精度^{[10, 11]}，因此本文选取其提供的数据进行建模。本文利用2005~2011年的PWV和ZWD两种格网数据的时间序列，求得各格网点上Π^-1的时间序列。通过分析其时空分布特征，发现Π^-1与测站高程高度相关，并具有年周期和半年周期特性，于是先顾及年周期与半年周期项以及高程改正对Π^-1进行拟合，得到相应的拟合系数，接着借鉴文献[4]对拟合系数进行球谐展开，最终建立了转换系数Π的全球经验模型。

确定转换系数Π的时空分布特征是建立精密的转换系数模型的基础和关键。在最新的 T_m的建模过程中，大都顾及了海拔高程对T_m的影响并予以改正，由于转换系数Π是T_m的函数，我们有理由相信，转换系数Π与海拔高程也是相关的。为了探讨转换系数Π与海拔高程的关系，本文利用2005~2011年的ZWD和PWV格网数据，求得Π^-1的时间序列，再对每个格网点求得Π^-17 a的平均值，并给出其与对应格网点海拔高程之间的关系图(见图 1)。

图 1 2005~2011年全球13 104个格网点的 Π-1的平均值与海拔高程的关系图 Fig. 1 Relation Between the Average Π-1of the Global 13 104 Grid Points and Altitude from 2005 to 2011

图选项

从图 1中可以看出，Π^-1的平均值在2~8之间变化，主要集中在4~6.5之间，如果按文献[2]直接近似取Π为常数0.15，则其误差比较大。同时，格网点大部分分布在较低海拔地区，在海拔超过4 000 m之后点明显较少。而且Π^-1的平均值与海拔高程之间存在一定的线性关系，这种关系在低海拔地区尤为明显。因此对Π建模时，需要对海拔高程的影响进行线性改正。

从文献[5]可以看出，Π^-1具有年周期，为了分析Π^-1值是否有半年周期，本文分别在较高纬度和中低纬度选取了4个格网点2005~2011年转换系数Π^-1的时间序列，通过傅立叶频谱分析给出其功率谱密度(power spectral density，PSD)，如图 2所示。

图 2 选取的4个点在2005~2011年转换系数的时间序列及相应的傅立叶频谱分析结果 Fig. 2 Time Series of the Four Selected Points’ Conversion Factor and Fourier Spectrum Analysis from 2005 to 2011

图选项

从图 2可以看出，(70°N，40°W)、(24°N，92.5°E)和(84°S，140°E)这三个点处的转换系数时间序列的PSD图有两个明显的峰值，分别在0.5和1位置处，说明这些点的转换系数具有明显的年周期和半年周期。其中，(24°N，92.5°E)和(84°S，140°E)的半年周期尤其明显；而对于点(18°S，72.5°W)来说，则只在1位置处有较明显的峰值，因而该点的转换系数只有较明显的年周期，而没有半年周期或半年周期不明显。由此可以看出，转换系数Π^-1的时间序列不仅具有年周期，而且存在一定的半年周期，在建模时考虑半年周期影响也是有必要的，这也与T_m具有年周期和半年周期完全吻合。而是否有日周期，经过计算发现，加上日周期对结果的改进只有千分之一mm级，增加更多的参数只得到这样的改进显然得不偿失，因此本文不考虑日周期对建模的影响。

通过上文分析，本文在对某个格网点Π^-1建模时顾及了年周期、半年周期以及海拔高程相关的改正，并估计了年周期项和半年周期项初相，最后得到一个新模型：

式中，α₀为平均值；α₁为年周期项振幅；α₂为半年周期项振幅；α₃为高程改正系数；C₁为年周期项初相；C₂为半年周期项初相；doy为年积日；h为海拔高，单位为m。

将式(3)展开并令a₀=α₀，a₁=α₁cos ，则有：

本文在求解模型系数的过程中采用了2005~2011年GGOS Atmosphere提供的ZWD格网数据和ECMWF提供的PWV格网数据时间序列，通过式(1)求得转换系数Π的时间序列并进行拟合，得到的模型系数在全球的分布如图 3所示，其中常系数a₅≈-0.000 433 3。根据式(4)构建误差方程，利用每个格网点的时间序列数据，采用最小二乘原理即可求解出模型系数。

图 3 模型系数全球分布图 Fig. 3 Global Distribution of the Model Coefficients

图选项

从图 3(a)可以看出，在青藏高原附近地区和南美洲西海岸等极个别地区转换系数的倒数平均值a₀比其他地区小了很多(只有3左右)，这可能对建模产生了一定的影响。同时在全球范围内，各拟合系数与经纬度的关系较为密切，不同位置处的系数也表现出一定的差异性。

通过上述建模过程，本文获得了每个格网点上的模型系数，但通常在转换系数Π值应用中，并不都是在格网点上，由此将所求得的格网系数表达成空间位置相关的函数是非常重要的。文献[12]首次使用球谐函数建立全球地表温度和气压模型，取得良好的效果。文献[4]中也是采用球谐函数建立GWMT模型，减少模型参数个数的同时满足GPS气象学反演水汽精度要求，这都说明球谐函数在表征地球物理参数中具有极大的优势和应用前景。于是本文也采用9阶9次的球谐函数将上述a₀,a₁,a₂,a₃,a₄这5个参数进行球谐展开，即：

式中，P_nm为勒让德多项式；φ和λ分别为格网点的纬度和经度；A_nmⁱ和B_nmⁱ为使用最小二乘原理确定的n阶m次球谐函数的系数。

于是根据式(4)、式(5)即可建立Π的全球经验模型，该模型无需站点气象数据，仅用较少的参数就可以通过站点经纬度、年积日和海拔高得到全球任意站点的Π值。

为了分析本文所建立的Π模型的有效性和适用性，验证其在全球范围内的精度和稳定性，本节将采用其他数据源数据对模型的精度进行外部检验。参与检验的数据包括未参与建模的2012年GGOS Atmosphere提供的ZWD格网数据、ECMWF提供的PWV格网数据及2012年无线电探空资料。

由于建模时采用的是2005~2011年的ECMWF提供的PWV格网数据与GGOS Atmosphere提供的ZWD格网数据，本文采用2012年的格网数据对模型进行检验。利用式(5)分别计算出每个格网点的模型系数a₀、a₁、a₂、a₃、a₄，再根据式(4)计算出对应时间上每个格网点的转换系数Π值，将GGOS Atmosphere提供的2012年每天0:00、6:00、12:00、18:00(UTC)等4个时刻的各个格网点上的ZWD通过式(1)计算得到PWV，并分别与ECMWF提供的相应格网点的PWV值进行比较分析，最后统计出模型精度检验结果，其Bias与RMS统计结果如表 1所示。

从表 1可以看出，在13 104个格网上，通过新模型计算转换系数所得的PWV与2012年ECMWF格网PWV之间的Bias的平均值为-0.179 mm，最大值为5.972 mm，最小值为-11.146 mm；而两者之间的RMS的平均值为1.806 mm，最大值为11.825 mm，最小值为0.059 mm。由于平均Bias值较小，故可以认为所建立的转换系数模型没有明显的系统性偏差；而平均RMS值也比较小，则说明模型精度较高。

图 4、图 5分别给出了Bias和RMS的整体区间分布情况以及全球分布图。

图 4 2012年格网数据检验Bias与RMS直方图和饼图 Fig. 4 Histograms and Pie Charts of Bias and RMS Differences Validated by Gridded Data in 2012

图选项

图 5 2012年格网数据检验Bias与RMS全球分布图 Fig. 5 Bias and RMS’s Global Distribution of Differences Validated by Gridded Data in 2012

图选项

从图 4中的 Bias和RMS的直方图可以看出，两者整体符合统计学规律。其中，Bias高度集中在0附近，主要在-2~2 mm范围内，而且随着Bias逐渐增大，所在区间的个数逐渐减少，这说明模型系统性偏差较小；RMS值则主要集中在0~4 mm，而大于10 mm的点只有两个。饼图更加直观地显示出在全球13 104个格网点中，Bias值在-2.5~2.5 mm范围内占总数的99%；而RMS值在0~4.5 mm范围内个数占总数的99% 以上，在0~3.5 mm范围内超过94%。这也说明模型与格网数据整体符合较好，在全球范围内整体系统误差较小，精度和稳定性较高。

从图 5可以看出Bias和RMS在全球大部分地区都比较小，其中Bias主要在-2~2 mm，RMS主要在4 mm以内。在青藏高原附近地区和南美洲西海岸等极个别地区较大，这与图 3(a)表现一致，验证了之前的想法，说明这些地区气候较异常而导致Π的值较其他地区大得多，从而拟合时会出现一定的误差，这也可能是由于ECMWF提供的水汽产品在这些地区精度本身有限而导致的。而青藏高原气候异常的原因则是由于海拔较高，大气环流较剧烈引起的，南美洲西海岸则是受到厄尔尼诺现象的影响。整体上看来，模型基本消除了经纬度对转换系数Π的影响，具有较高的精度。

本文在全球范围内选取了661个无线电探空站2012年的每天0:00和12:00(UTC)探空资料，在利用无线电探空仪的观测数据时，剔除了那些观测层数过少或观测高度过低(在中纬度地区，观测高度一般应大于12 km)的数据，以免产生较大的离散误差或忽略了大气较高处的水汽信息，分别按照式(6)、式(7)计算出各点的ZWD^[13] 以及PWV^[14]，并以此为真值，进行检验与分析。对于不同测站全年的检验结果而言，Bias与RMS统计结果如表 2所示。

式中，为第i层的湿折射率；n为观测层数；P_wⁱ为第i层大气的平均水汽压(hPa)；T_i为第i层大气的平均温度(K)；Δh为每一层的高度。

从表 2中可以看出，通过本文模型计算转换系数所得的PWV与无线电探空资料直接计算的PWV之间的Bias的平均值为0.465 mm，最大值为5.547 mm，最小值为-2.756 mm；而两者之间的RMS的平均值为0.789 mm，最大值为0.102 mm，最小值为5.629 mm。由于平均Bias值和RMS值都比较小，同样可以认为模型没有明显的系统性偏差，精度较高。该检验从整体上再次证明本文所建模型在全球范围内具有较高的精度与稳定性。

为了进一步分析模型在全球范围内的精度变化情况，本文给出了这661个探空站点的Bias与RMS全球分布图，如图 6所示。

图 6 2012年无线电探空数据检验Bias与RMS全球分布图 Fig. 6 Bias and RMS’s Global Distribution of Differences Validated by Radiosonde Data in 2012

图选项

从图 6中可以看出，在全球范围内模型所得PWV与探空站PWV两者的Bias和RMS都比较小，Bias集中在-2~2 mm范围内，RMS则主要集中在0~2.5 mm。这也说明模型已基本消除经纬度的影响，全球范围内具有较强的适用性。仅在极少部分地区RMS超过4 mm，而这些地区与之前用格网数据检验完全一致，再次说明这些地区的数据较为异常，从而导致对建模有一定的影响。

本文利用2005~2011年的GGOS Atmosphere提供的2.5°×2°的ZWD格网数据和ECMWF提供的2.5°×2°的PWV格网数据建立了两者之间转换系数的全球经验模型，该模型顾及了年周期、半年周期以及海拔高程相关的改正，仅需知道站点的经纬度、高程以及年积日即可得到转换系数Π。经过检验，采用格网数据Bias的平均值为-0.179 mm，RMS的平均值为1.806 mm；采用无线电探空资料平均Bias为0.465 mm，平均RMS为0.789 mm。这说明新建立的模型在全球范围内反演PWV系统性偏差较小，精度较高。与先计算T_m再利用式(2)得到转换系数相比，本文模型更加直接地求得转换系数，计算更为方便且更直观地显示ZWD与PWV两者之间的关系，可用于无气象数据条件下求得转换系数，从而进行PWV反演。

致谢: 感谢ECMWF提供全球PWV格网数据和GGOS Atmosphere提供全球ZWD格网数据。