大气水汽是诸多极端天气系统形成的重要物质基础和主要动力因子，因此研究大气水汽的分布和变化特征对于极端天气的研究和预警预报具有十分重要的价值，然而极端天气的产生和演变过程十分迅速，传统的水汽观测手段(如无线电探空仪)存在时间分辨率过低而无法及时地捕获极端天气中水汽的变化信息等缺陷。微波辐射计虽能获得相对连续的观测数据，但设备较为昂贵，且地面测站数目有限，无法对大范围区域的水汽进行有效检测。自Bevis等^[1]提出用高时间分辨率的连续观测的GPS数据反演大气可降水量(precipitable water vapor, PWV)以来，地基GPS气象学迅速发展，通过地基GPS获取的数据的水汽观测精度与探空测站和微波辐射仪的相当^{[2, 3]}，利用地基连续GNSS观测资料反演大气水汽的技术成为弥补传统观测手段的有效手段。地基GPS水汽反演因其高时间分辨率、高精度、低费用等优势，被广泛用于改善数值天气预报模式^{[4, 5]}以及用于气候和特殊天气的研究^[6-9]。地基GPS气象学基于卫星观测方程估算的对流层天顶总延迟(zenith total delay, ZTD)来反演大气水汽，反演时首先将地面气压观测值代入天顶静力延迟模型求解天顶静力学延迟(zenith hydro‐static delay, ZHD)，而后用ZTD减去ZHD可以得到天顶湿延迟(zenith wet delay, ZWD)，再将ZWD乘以一个转换系数，可以得到PWV。GPS反演PWV需要地面实测的气象参数，如计算ZHD时需要测站处气压，计算加权平均温度时需要地表实测温度。但地基跟踪站并不是专门为气象建立的，因此国内外的众多地基GPS站点并未配备气象传感器。为充分发挥地基测站积累的大量观测资料的作用和价值，有不少学者通过其他方式来获取地面气象参数。刘智敏等^[10]通过GPT2模型^[11]反演出了山东区域连续运行基准站(continuously operating reference station, CORS)网的PWV; Bai等^[12]通过内插澳大利亚气象数据网络的方法获取了测站处的温压；Jade等^[13]利用美国国家环境预报中心及国家大气研究中心(National Centers for Environmen‐tal Prediction/National Center for Atmospheric Re‐search, NCEP/NCAR)再分析资料反演了印度地区的PWV；常亮等^[14]同样针对NCEP/NCAR再分析资料，提出一种GPS测站附近无气象传感器时估算PWV的方法。以上研究都是借助全球气温气压经验模型或第三方资料来确定地面气压和气温，从而反演PWV。曲建光等^[15]首次提出一种无需气象数据，直接由ZTD推算PWV的方法，并证明了其可行性；王勇等^[16]基于ZTD和PWV的高相关系数建立了两者线性回归的转换模型；易正晖等^[17]以线性回归方程建立了ZTD直接推导PWV的季节性和全年性转换模型。

上述以线性回归建立的ZTD和PWV的转换模型从数据相关性方面进行研究，本文则基于以上研究，从数学原理层次提出一种无需气象数据，由ZTD直接推导PWV的新方法。该方法基于最小二乘算法，从GPS反演大气水汽的反演方程出发，建立ZTD推算PWV的模型，即建立自变量为ZTD，因变量为PWV，其他参数或系数为常数的方程。

1 模型推导

地基GPS气象学以非差或双差模型获取ZTD，而后将测站实测气压代入求解ZHD的模型，求得ZHD，其中，计算ZHD的模型有Hopfield、Saas‐tamoinen等模型，本文选用Saastamoinen模型，其表达如下：

$ \left\{ {\begin{array}{*{20}{l}} {{\rm{ZHD}} = \left( {2.2768 \pm 0.0024} \right) \times \frac{{{P_s}}}{{f\left( {\theta , H} \right)}}}\\ {f\left( {\theta , H} \right) = 1 - 0.00266 \times \cos 2\theta + 0.00028H} \end{array}} \right. $

(1)

式中，P_s表示气压；θ表示测站纬度；H表示测站大地高。

用ZTD减去ZHD获得ZWD:

$ {\rm{ZWD}} = {\rm{ZTD}} - {\rm{ZHD}} $

(2)

ZWD乘以一个转换系数Π得到PWV:

$ \left\{ {\begin{array}{*{20}{l}} {{\rm{PWV}} = \mathit{\Pi } \times {\rm{ZWD}}}\\ {\mathit{\Pi } = \frac{{1 \times {{10}^6}}}{{{\rho _W}\left( {{k_3}/{T_m} + {{k'}_2}} \right){P_V}}}}\\ {{T_m} = 70.2 + 0.72{T_s}} \end{array}} \right. $

(3)

式中，ρ_W为液态水密度；T_m为大气加权平均温度；k'₂、k₃为大气折射常数；R_V为水汽的比气体常数；T_s为地表温度。

由式(3)可知，转换系数Π为温度的函数，则由式(1)~式(3)可得：

$ {\rm{PWV = }}\mathit{\Pi }\left( {{T_s}} \right) \times {\rm{ZTD - }}\mathit{\Pi }\left( {{T_s}} \right) \times {\rm{ZHD}}\left( {{P_s}, \theta , H} \right) $

(4)

式中，ZHD为气压、测站纬度和大地高的函数，对其进行泰勒展开，可得：

$ \begin{array}{c} {\rm{PWV}} = {\rm{PW}}{{\rm{V}}_0} + \left( {{\rm{ZTD}} - {\rm{ZHD}}\left( {{P_s}, \theta , H} \right)} \right) \times \mathit{\Pi '}\left( {{T_s}} \right)\left| {_{{T_0}}} \right. \times {\rm{d}}{T_s} - \mathit{\Pi '}\left( {{T_s}} \right) \times {\rm{ZHD'}}\left( {{P_s}, \theta , H} \right)\left| {_{{P_0}}} \right. \times {\rm{d}}{P_s} - \mathit{\Pi }\left( {{T_s}} \right) \times \\ {\rm{ZHD'}}\left( \theta \right)\left| {_{{\theta _0}} \times {\rm{ d}}\theta - \mathit{\Pi }({T_s}) \times {\rm{ZHD'}}(H)\left| {_{{H_0}}} \right. \times {\rm{d}}H} \right. \end{array} $

(5)

式中，P₀、T₀、θ₀、H₀分别为P_s、T_s、θ、H的近似值；PWV₀为将近似值参数代入式(4)后的PWV的近似值。

令

$ \left\{ {\begin{array}{*{20}{l}} {a = ({\rm{ZTD}} - {\rm{ZHD}}({P_s}, \theta , H)) \times \mathit{\Pi '}({T_s})\left| {_{{T_0}}} \right.}\\ {b = - \mathit{\Pi '}({T_s}) \times {\rm{ZHD'}}({P_s}, \theta , H)\left| {_{{P_0}}} \right.}\\ {c = - \mathit{\Pi }({T_s}) \times {\rm{ZHD'}}(\theta )\left| {_{{\theta _0}}} \right.}\\ {d = - \mathit{\Pi }({T_s}) \times {\rm{ZHD'}}(H)\left| {_{{H_0}}} \right.} \end{array}} \right. $

(6)

则有：

$ \begin{array}{c} {\rm{PWV}} = {\rm{PW}}{{\rm{V}}_0} + a \times {\rm{d}}T + b \times {\rm{d}}{P_s} + c \times {\rm{d}}\theta + \\ d \times {\rm{d}}H \end{array} $

(7)

地基测站无实测气象数据时，会把无线探空仪的PWV作为观测值，大气中真实的PWV作为真值，则方程可变为：

$ \mathit{\boldsymbol{v}} = \mathit{\boldsymbol{BX}} - \mathit{\boldsymbol{L}} $

(8)

式中，$\left\{ {\begin{array}{*{20}{l}} {\mathit{\boldsymbol{B}} = \left[ {\begin{array}{*{20}{c}} a&b&c&d \end{array}} \right]}\\ {\mathit{\boldsymbol{X}} = {{\left[ {\begin{array}{*{20}{c}} {{\rm{d}}{T_s}}&{{\rm{d}}{P_s}}&{{\rm{d}}\theta }&{{\rm{d}}H} \end{array}} \right]}^{\rm{T}}}, }\\ {\mathit{\boldsymbol{L}} = {\rm{PWV}} - {\rm{PW}}{{\rm{V}}_0}} \end{array}} \right.$

由最小二乘法可得：

$ \mathit{\boldsymbol{X}} = {\left( {{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{PB}}} \right)^{ - 1}}{\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{PL}} $

(9)

从而确定4个常量参数：

$ \left[ {\begin{array}{*{20}{c}} {{T_s}}\\ {{P_s}}\\ \theta \\ H \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{T_0}}\\ {{P_0}}\\ {{\theta _0}}\\ {{H_0}} \end{array}} \right] + \mathit{\boldsymbol{X}} $

(10)

将式(10)代入式(4)求得PWV，该模型的形式即为将P_s、T_s、θ、H看作满足最小二乘常数的结果。

2 数据获取及预处理

本文模型使用ZTD和探空数据推导，并通过以GPS测站处实测气温气压反演得到的PWV对模型推导的PWV进行精度检验。以国际GNSS服务(International GNSS Service, IGS)中国区域的BJFS测站为例，为获取该站点的ZTD和独立的PWV^[18]，本文引入了4个超500 km的长基线测站(DAEJ、WUHN、SHAO、SUWN)^[19]。本文使用GAMIT采用松弛解模式解算2016年年积日第214~300天的数据，其中，第218天数据缺失，共计86天数据。

BJFS测站并未配备探空设备，因此，本文将怀俄明大学提供的中国境内的ZBAA探空站的探空数据作为BJFS探空数据的补充，BJFS与ZBAA两者直线距离33.2 km，高程相差22.4 m，获取的探空数据按照以下公式计算PWV:

$ \left\{ {\begin{array}{*{20}{l}} {W = \int {\left( {q{\rm{/}}g} \right){\rm{d}}p} }\\ {q = \frac{{622e}}{{p - 0.378e}}}\\ {e = 6.1078 \times {{\rm{e}}^{\frac{{a \times {T_d}}}{{b + {T_d}}}}}} \end{array}} \right. $

(11)

式中，W表示可降水量(单位：mm); q表示比湿(单位：g/kg); g为地球的重力加速度(单位：m/s²); p、e、T_d分别为各高度上的气压(单位：hPa)、水汽压(单位：hPa)和露点温度(单位：℃)。系数a、b的取值为：气温低于-40℃时，a=21.87, b=265.49；气温高于0℃时，a=17.26, b=237.29；若气温居于两种情况之间，系数a、b用线性内插算出。实际计算时以各个标准面等压面上的比湿值进行差分计算以代替积分，即从地面(P=P₀)到大气上界(P=0)的可降水水量计算公式为：

$ W = - \frac{1}{g}\sum\limits_{{P_0}}^0 {q\Delta p} $

(12)

探空测站一天只有0时和12时的数据，故共有172个样本。其中，前半部分用来进行模型推导，后半部分用来进行检验，将GPS反演的PWV(GPS/PWV)与探空测站计算得到的PWV(Radio/PWV)进行比较，比较结果如图 1所示。由图 1(b)可知，GPS/PWV与Radio/PWV的均方根(root mean square, RMS)值为11.2 mm，相关系数仅为0.763 2。从图 1(a)中可以看出，在年积日第284~294天的某个时间弧段存在明显的粗差，这是由于探空测站存在系统误差且GAMIT不能完全检测出GPS观测资料的质量问题，因此，有必要在推导模型前对原始PWV数据进行预处理，检核并剔除掉粗差。

极限误差通常用来做剔除粗差的阈值，其值一般为3倍的中误差σ，本文将GPS/PWV与Radio/PWV较差的标准差(standard deviation, STD)作为中误差，以3σ为检核的极限误差，数据预处理结果见表 1。

由表 1可看出，当设置极限误差为3σ且仅检核一次时，剔除了3个样本。从图 1(a)中可以看出，较差超过20 mm的不只3个样本，以3σ迭代检核时剔除了6个样本，相关系数为0.973，相关性较强且剔除的粗差较为合理。因此，将3σ迭代的结果作为后续的基础数据，3σ迭代的预处理结果如图 2所示，GPS/PWV与Radio/PWV两者曲线基本重合。

3 解算结果及精度检验

将上述预处理后的Radio/PWV取前半部分进行最小二乘解算，解算时将Radio/PWV作为观测值，将经GAMIT解算得到的ZTD作为无误差的真值，以前后两次中误差σ的差值是否小于限差为一重迭代条件，为进一步剔除预处理遗漏的粗差，设置3σ作为筛选数据的二重迭代条件，ZTD直接转换PWV的流程图见图 3。将GAMIT解算得到的ZTD及探空测站计算得到的PWV代入式(10)，可得到BJFS站ZTD转换为PWV的模型：

$ {\rm{PWV}} = - 390.900 + 0.171 \times {\rm{ZTD}} $

(13)

为检验模型的可靠性，本文将模型计算的PWV(Model/PWV)与Radio/PWV进行比较，故取预处理数据的后半部分代入模型计算。Model/PWV与对应Radio/PWV的比较如图 4所示。

由图 4可知，两种曲线的变化趋势基本一致，RMS为5.8 mm，相关系数为0.898 7。为对其精度有一个定量的认识，将Model/PWV和Radio/PWV比较的精度与GPS/PWV和Radio/PWV比较的精度进行对比，其中GPS/PWV与Radio/PWV比较如图 5所示。

对比图 4(a)、图 5(a)可知，相较于Model/PWV, GPS/PWV的曲线更加贴合Radio/PWV。由图 5(b)可知，GPS/PWV与Radio/PWV的RMS为3.3 mm，相关系数为0.935 2。Model/PWV与Radio/PWV的RMS比GPS/PWV与Radio/PWV的大了2.5 mm，相关系数较之小了0.036 5，主要原因可能是：BZAA探空站与BJFS站不并址，存在偏差，如果考虑BZAA探空站的PWV水平方向和垂直方向的改正，则可能优化其精度；2016年年积日214~300日处在夏秋两季之间，未考虑ZTD和PWV的季节性变化，可能对模型造成了干扰。

将Model/PWV与GPS/PWV进行比较，进一步评定模型解算的精度，如图 6所示。

由图 6(a)可知，GPS/PWV始终位于Model/PWV下方，存在一个系统性偏差。这是因为模型推导使用的是探空数据，而非GPS数据，而探空数据与GPS数据本身又存在误差。但GPS/PWV与Model/PWV的变化趋势几乎一致，这是GPS/PWV与Model/PWV的RMS(4.5 mm)虽大于GPS/PWV与Radio/PWV的RMS(3.3 mm)，但是相关系数仍高达0.981 8的原因。Model/PWV、GPS/PWV、Radio/PWV三者的比较结果见表 2。

从表 2可以看出，GPS/PWV与Radio/PWV的较差绝对值仅为0.121 mm，由图 5(b)可知，Radio/PWV的点在直线两侧分布较均匀，故较差值较小；Model/PWV与Radio/PWV的BIAS与GPS/PWV和Radio/PWV较差之间的较差为4.246 mm，与Model/PWV和GPS/PWV较差的较差精度相当，均在4 mm左右。

综上，Model/PWV和Radio/PWV的RMS与GPS/PWV和Radio/PWV的RMS之间的较差优于3 mm，较差之间的较差优于5 mm，这是探空站和GPS测站异址且未考虑季节性变化的缘故；Model/PWV与GPS/PWV的相关系数高达0.982，两者整体存在一个系统性偏差，曲线变化趋势几乎一致，在不研究其数值大小而研究其趋势变化时，可以用Model/PWV替代GPS/PWV，若考虑改正系统偏差则可进一步提高模型精度。

4 结束语

本文提出了一种无需气象数据，由ZTD直接推导PWV的新方法，并对其精度进行了检验，得到以下结论：

1) Model/PWV与Radio/PWV的RMS为5.773 mm，其与GPS/PWV和Radio/PWV的RMS较差优于3 mm，这是因为探空站与GPS测站异址，且未考虑季节性变化，若考虑季节性变化及探空站和GPS测站异址改正，则可能进一步提高模型精度。

2) Model/PWV与GPS/PWV的相关系数高达0.982，两者存在一个微小的系统性偏差，在不研究其数值大小只研究其趋势变化时，可以用本文模型替代GPS反演模型。若本文模型考虑添加一个系统偏差性改正，则可以进一步提高其精度。