高精度的斜路径水汽延迟(SWV)、精确确定“水汽层层顶”以及合理地划分垂直方向层析网格，是提高地基GNSS反演水汽三维信息精度的3个关键技术^[1]。高精度的SWV主要取决于斜路径湿延迟和湿度转换因子。通常情况下，利用非差或者双差的方法可以获取高质量的信号斜路径湿延迟^[2]，因此实时获得高精确湿度转化因子是获得高质量SWV的关键因素^[3]。湿度转换因子一般可以表示为随大气加权平均温度变化的函数，故求解出高精度大气加权平均温度是获得高质量湿度转换因子的前提。“水汽层层顶”的高度通常可以根据高质量、高垂直分辨率的大气产品(如无线探空、掩星、EAR-Interim等)来精确确定^[3]。然而，如何科学划分垂直方向层析网格目前还没有可靠的准则和依据^[4]。

ERA-Interim再分析资料可通过网站http://apps.ecmwf.int/datasets/data/interim-full-daily/免费下载，其时间分辨率为6 h，按照压强在空间上从1 000~1 hPa的变化分37层，提供了大气压强、大气温度、湿度信息等大气元素。本文利用ERA-Interim再分析资料提供的大气产品，首先建立区域性大气加权平均温度模型；然后根据水汽压的空间和时间上的分布信息，确定“水汽层层顶”的高度和“水汽转折层”的高度，并提出一种新的划分垂直方向层析网格的方法。通过对香港地区2014-06的实测数据进行试算，验证了优化后的层析方法相对于传统层析方法的优势。

1 水汽层析的关键技术 1.1 确定湿度转换因子

获取斜路径水汽含量SWV需要将估算的湿度分量SWD转换成SWV^[1]：

$ {\rm{SWV = }}\mathit{\prod} \cdot {\rm{SWD}} $

(1)

式中，Π为湿度分量转换因子。Π是对流层大气加权平均温度T_m的函数，在不同的地区和季节，其值不同^[2]：

$ \mathit{\prod} {\rm{ = }}\frac{{{{10}^6}}}{{{\rho _w} \cdot \frac{R}{{{m_w}}}.\left[{\frac{{{k_3}}}{{{T_m}}} + {k_2}-\frac{{{m_w}}}{{{m_d}}} \cdot {k_1}} \right]}} $

(2)

$ {T_m} = \frac{{\int_{{h_0}}^\infty {\left( {\frac{{{P_w}}}{T}} \right) \cdot {\rm{d}}h} }}{{\int_{{h_0}}^\infty {\left( {\frac{{{P_w}}}{{{T^2}}}} \right) \cdot {\rm{d}}h} }} = \frac{{\sum {\frac{{\left( {{h_2}-{h_1}} \right){P_w}}}{T}} }}{{\sum {\frac{{\left( {{h_2}-{h_1}} \right){P_w}}}{{{T^2}}}} }} $

(3)

式中，ρ_w为液态水的密度，k₁、k₂、k₃分别为常数，k₁=77.6 K/hPa、k₂=70.4 K/hPa、k₃=3.739×10⁵ K/hPa, m_d、m_w分别为干大气与水汽的摩尔质量, R为普适气体常数, P_w为水汽压(单位hPa), T为温度(单位K)。

由式(3)可知，T_m可以通过高垂直分辨率的大气产品来估算，此外还可以对其建模，然后近实时地估算出T_m。目前较为常用的是Bevis等^[2]建立的T_m模型：

$ {T_m} = 70.2 + 0.72{T_s} $

(4)

式中，T_s为地表温度。Bevis给出的T_m模型适用于全球，但只考虑了地面温度的影响。为了提高T_m的计算精度，本文利用香港地区2011~2013年的ERA-Interim产品，根据式(3)计算出T_m，顾及其与时间变量D的余弦关系和地表温度T_s的线性关系的影响，建立香港地区区域性T_m模型，如式(5)：

$ \begin{array}{l} f\left( {D, {T_s}} \right) = 136.2 + 0.705\;8 \times \\ \cos \left( {\frac{{-3.811\left( {D-30} \right)}}{{365.25}}} \right) + 0.505\;9{T_s} \end{array} $

(5)

式中，D为年积日。

图 1给出了利用ERA-Interim产品建立的T_m模型的结果。

从图 1可以看出，联合时间参数及地表温度参数建立的T_m模型能够较好地拟合出T_m的变化。经统计，该模型的内符合精度优于2.31 K。为了检验该模型外符合精度，选取2014年探空站45004的探空数据，利用式(3)计算T_m序列及其与探空产品计算的T_m的差值。此外，计算Bevis模型及Wang等^[5]给出的香港区域性模型(T_m=113.29+0.586 3T_s)计算的T_m与探空资料估算的T_m的差值(图 2)。

从图 2可以看出，利用本文新模型估算的T_m与真值相比整体偏大。经统计，新模型估算的T_m与真值之差的正偏差最大值优于7 K，负偏差最大值优于-4 K。新模型估算的T_m与Radiosonde资料的计算值之间的差异主要来源于模型本身存在的系统误差。二者之间的差值部分情况下为负值，经过核查知，Radiosonde资料观测值的高度有限(小于10 km)，从观测最高点到对流层层顶这段高度处的气象观测数据未给出。T_m的精度直接影响到湿度转换因子的精度，并最终影响到SWV和PWV产品的质量。此外本文又整体统计了N_Tm、B_Tm、W_Tm的值，如表 1所示(单位K)。

1.2 探测水汽层层顶

地基GNSS层析网格最低高度和最高高度通常分别为地球海平面和对流层层顶，但是实际上水汽密度在垂直方向主要分布在某个高度以下，本文记为“水汽层层顶”。从该处到对流层层顶处的水汽密度非常小，接近为0^[3]。在层析反演中，这一段层析网格内的水汽密度值解算的结果通常为负数。因此，在层析反演中，可以将层析网格的上限确定为“水汽层层顶”的高度，从而降低了层析网格的高度，可以增加穿越层析网格的GNSS射线的数目，同时也减少了未知网格的个数，从而能够提高水汽反演的精度。“水汽层层顶”的高度通常可以通过大气可降水量PWV和水汽密度垂直分布情况确定。图 3给出了2013年香港地区ERA-Interim产品(0.15°×0.15°)计算的不同高度处的PWV年变化量以及水汽密度的年变量。

图 3显示，PWV的垂直变化量及时间变量与水汽密度的垂直变化量及时间变化量趋于一致。根据这个现象，联合PWV及水汽密度的垂直分布信息，给出一种确定水汽层层顶的方法：1)按照压强进行分层，求出每一层高度上PWV占整个高度层PWV的比值RPWV; 2)寻找出RPWV≤0.01点的高度范围；3)寻找水汽密度小于0.1 g/cm³的高度范围；4)将RPWV≤0.01且水汽密度小于0.1 g/cm³的高度定义为水汽层层顶。根据这个定义，给出2013年香港地区水汽层层顶随时间的变化规律，如图 4所示。

从图 4可以看出，水汽层层顶随季节变化的规律很明显。冬季空气比较干燥，水汽层层顶的值偏低，大致的范围为6~8 km；夏季空气比较湿润，水汽层层顶的高度偏高，大致的范围为9~11 km。在干燥的季节，可以取每半个月的平均值作为水汽层层顶的高度；在湿润的季节，可以取每个月的平均值作为水汽层层顶的高度。

1.3 垂直方向层析网格划分

GNSS水汽层析的目的是获取水汽的三维分布信息，故合理地划分垂直方向层析网格是改善层析结果的重要因素之一。Flores等^[6]的研究表明，GPS层析的分辨率的最小垂直厚度约为300 m, 如果小于这个厚度，则噪声会影响层析的结果，但其并未给出层析最佳的垂直分辨率。为此，本文提出一种新的垂直方向划分网格的方法。

在垂直方向，PWV随高度的分布近似呈指数形式变化：

$ {\rm{PWV = PW}}{{\rm{V}}_0} \cdot \exp \left( {\frac{{h-{h_0}}}{H}} \right) $

(6)

由式(6)可得：

$ h = H \cdot \left( {\ln \left( {{\rm{PW}}{{\rm{V}}_T}} \right)-\ln \left( {{\rm{PW}}{{\rm{V}}_0}} \right)} \right) + {h_0} $

(7)

式中，PWV_T为水汽层层顶处的PWV，PWV₀为地表处的PWV。

$ H = \frac{{{h_r}-{h_0}}}{{\ln \left( {{\rm{PW}}{{\rm{V}}_r}} \right)-\ln \left( {{\rm{PW}}{{\rm{V}}_0}} \right)}} $

(8)

由式(7)可知，垂直方向层析网格的划分可以通过PWV比值RPWV=PWV_i/PWV₀(i=1, 2, …)来确定。图 5给出了2013年第30天的一个RPWV的剖面。

由图 5看出，PWV比值随高度呈现出指数形式的变化。在水汽层层顶处，PWV接近于0；在接近6 km处，PWV出现了转折，定义为水汽转折层。确定水汽转折层的标准为：首先确定PWV比值随高度呈现出指数变化，并估算该指数模型的斜率K，斜率随高度从大到小变化；然后求出地面点与水汽层层顶的连线与X轴的夹角的正切值K_R，当某个高度处的K与K_R无限接近时，则认为该高度为水汽转折层的高度。据此，利用ERA-Interim产品给出2013年香港地区水汽转折层的高度随时间的变化量，如图 6所示。

图 6中水汽转折层的变化规律与图 5中水汽层层顶的高度变化规律相似。同理可得，在干燥季节把每半个月的平均值作为水汽转折层；在湿润季节把每个月的平均值作为水汽转折层。

根据水汽层层顶的高度及水汽转折层的高度，定义垂直方向网格的划分方法。从地面到水汽转折层的高度上，层析网格划分的要求为：网格的高度不小于300 m，且不大于600 m，同时还要满足每一层独立的网格层底和层顶的PWV不能超过12%。在水汽转折层到水汽层层顶的高度上，层析网格划分要求为：网格的高度不小于600 m，且不大于1 000 m。每个月求一个平均的水汽层层顶和水汽转折层的高度。

2 实验结果

选取香港地区12个CORS站2014-06的观测数据及气象观测数据，利用Bernese软件进行处理，获得对流层天顶总延迟ZTD，然后选用Saastamoinen模型及地面气象数据估算出对流层天顶静水延迟ZHD。在水平方向层析网格经度范围为113.87°~114.35°E，纬度范围为22.19°~22.54°N。水平网格在经度方向为0.06°，纬度方向为0.05°，故水平方向共划分7×8=56个网格。在垂直方向，使用上文方法，利用2010~2013年香港地区ERA-Interim数据计算出6月的水汽层层顶为9.2 km, 水汽转折层层顶为6 km。从地面到6 km分13层，其中第1层为400 m，第2层到第9层每层高度为450 m，第10层到第13层每层高度为500 m。6~9.2 km分4层，每层的高度为800 m。然后根据Flores等^[6]提供的GNSS水汽层析的方法，每0.5 h输出一次反演结果。为了检验本文优化后层析结果的优越性，与传统的层析方法进行比较。传统的层析方法为：湿度转换因子取0.152，水汽层层顶选为10 km，层析网格选择平均分层的方法，即每层的高度为500 m。限于篇幅，本文仅给出2014年第157天的反演结果，如图 7所示。

由图 7可以看出，2014年第157天的水汽主要集中在4 km以下，随着时间的变化，0.4~ 3.1 km处的水汽密度发生较为明显的变化。为了检验近实时估算水汽密度的精度，将本文估算的水汽密度值记为TomoWV1，传统层析方法的解算结果记为TomoWV2，然后将两者分别与探空产品进行比较。选择第45004探空站2014-06的产品，并将获得的水汽密度值记为RadWV。图 8给出了2014-06 RadWV与TomoWV1、TomoWV2的差值。

由图 8可以看出，由于水汽“逆增层”现象的存在，TomoWV1与RadWV的最大偏差为6 g/m³，TomoWV2与RadWV的最大偏差达到8 g/m³。在6 km以下，TomoWV1、TomoWV2与RadWV的偏差较大；而在6 km以上，TomoWV1、TomoWV2与RadWV的偏差相对较小。分别分段统计6 km以下和以上的层析水汽值与探空水汽值的差值，如表 2所示。

从表 2看出，将探空产品作为真值，优化后的层析结果与传统的层析结果相比，6 km以下精度提高了12%，6 km以上提高了17%。

3 结语

利用ERA-Interim再分析资料，首先对大气加权平均温度进行精密建模，获得实时湿度转换因子；然后根据ERA-Interim产品提供的水汽值空间分布信息，研究水汽层层顶的变化规律，为水汽层析确定最优的上边界层；最后提出了一种新的划分垂直方向层析网格水汽层析方法。为了检验这些优化方法对GNSS层析反演的影响，选取2014-06的数据进行反演，并与传统层析反演结果进行比较。选择无线探空产品作为真值，将层析结果与探空产品提供的水汽密度进行比较。结果表明，6 km高度以下优化后的水汽层析方法相对于传统层析方法获得的水汽密度精度提高了12%，在6 km以上提高了17%。