联合地基GNSS及空基GNSS掩星探测水汽三维分布

吴 昊，鄂盛龙，夏朋飞，许海林，罗颖婷，谭理庆

(1.广东电力科学研究院，广州 510080；2.武汉大学GNSS技术研究中心，武汉 430072)

0 引言

全球导航卫星系统(Global Navigation Satel-lite System，GNSS)能够有效地服务于全球导航与定位，此外，GNSS技术又扩展应用到全球气候监测与极端天气的预报[ 1-3]。目前，利用GNSS技术探测地球大气已成为国内研究的热点之一[ 1-5]。作为一种新的大气探测手段，GNSS技术能够准实时、无需人为干扰、全天候地服务大气科学的研究[ 6-10]。根据观测模式的不同，GNSS探测技术可以分为空基GNSS技术和地基GNSS技术[ 8-10]。空基GNSS掩星技术能够提供高垂直分辨率、高精度、全球覆盖的中性大气层到电离层的多级大气剖面产品。然而，受制于卫星星座的分布及数量，目前空基掩星技术还无法实时地服务于全球数值天气预报模式[10]。

地基GNSS系统探测大气是利用平差的方法从观测数据中估算出天顶对流层延迟(Zenith Tropospheric Delay，ZTD)，将其作为研究对流层大气的背景场。通常ZTD可以分为2个部分，即天顶对流层干延迟(Zenith Hydrostatic Delay，ZHD)和天顶对流层湿度延迟(Zenith Wet Delay，ZWD)。ZWD中含有大气水汽的信息，因此，利用一定的数学模型可以从ZWD中反演出大气水汽的分布[ 8-10]。

为了提高GNSS技术探测大气水汽分布的质量，文中试图利用无线电掩星产品辅助地基GNSS技术反演水汽三维分布信息。对流层天顶干延迟模型、湿度转换因子以及层析模型上边界的高度是影响层析结果的3个关键的参数。首先，根据无线电掩星提供的高垂直分辨率的大气产品标定了对流层干延迟模型；然后，建立新的区域性大气加权平均温度模型，提高了湿度转换因子的精度；最后，根据掩星提供的大气产品确定了水汽层层顶，将其作为层析模型的上边界层。利用香港地区2017年6月份的数据进行试算和分析，反演结果验证了上述方法的可行性。

1 标定对流层干项延迟及湿度转换因子

与对流层湿项延迟相比，对流层干项变化较为平稳，通过建立经验模型可以较为精确的估算。国内外研究者建立了许多全球或区域性的天顶对流层干延迟经验模型。其中，较为常用的主要包括如下3个：Saastamoinen模型、Hopfield模型和Black模型[ 11-13]。然而，经典的对流层干延迟模型存在着系统性的偏差，可以利用高精度的气象产品对其进行修正，如探空和掩星产品。利用高精度气象产品估算ZHD的数学公式如下[9]

(1)

(2)

其中，k1为常数，其经验值通常取为k1=77.6K/mbar；Pd为大气干压强；Pw为大气水汽压；T为温度；Rd和Rw分别为干普适气体常数和湿普适气体常数。

选取2007—2016年中六月份与香港地区共址的COSMIC (The Constellation Observing System for Meteorology, Ionosphere, and Climate)和GRACE (Gravity Recovery and Climate Experiment)掩星事件，将wetPrf产品利用式(1)和式(2)的计算结果分别与三种经典ZHD模型的计算结果进行了比较，其统计结果如表1所示。

表1 掩星产品计算的ZHD与经典ZHD模型计算的ZHD的差Tab.1 Bias of ZHD between radio occultation-derived and ZHD models-derived cm

从表1可以看出，将无线电掩星产品估算的ZHD(ROZHD) 视为参考值，三种经验模型与参考值之间都存在一定的偏差。其中，Sasstamoinen模型估算的ZHD (SAAS) 与参考值之间的偏差优于0.93cm，这个精度远远高于Hopfield模型和Black模型估算的ZHD(HOPE和BLA)。此外，HOPF的精度要优于BLA。为了提高ZHD的估算精度，文中试图利用Hopfield模型的表达形式来修正Sasstamoinen模型。

(3)

这里，a、b和c为区域性标定模型的系数。通过2007—2016年6月份与香港地区共址的COSMIC和GRACE掩星产品来估算这3个参数。精确获取ZHD后，便可从ZTD中分离出ZWD，利用湿度转换因子便可得到大气水汽含量，而湿度转换因子是关于大气加权平均温度的函数[14-15]

(4)

(5)

其中，ρv为水的密度，且ρv=103kg/m3；k1、k2、k3都为常数，取值分别为k1=77.6K/hPa，k2=70.4K/hPa，k3=3.739×105K/hPa；md和mw分别为干和湿大气的摩尔质量，取值分别为md=28.96kg/kmol，mw=18.02kg/kmol；R表示普适气体常数，取值为R= 8314Pa·m3·K-1·kmol-1。

根据所选取的香港地区2007—2016年的COSMIC和GRACE无线电掩星产品(RO)，利用式(5)便可估算出大气加权平均温度Tm，建立二阶拟合关系如图1所示。

图1 大气加权平均温度二阶拟合图Fig.1 Second-order fitting of atmospheric weighted average temperature

从图1可以看出，利用二阶模型拟合的大气加权平均温度模型与掩星产品获得的Tm具有较好的一致性，二者偏差的RMS值优于2.43K。图1拟合出的结果便为相关地区大气加权平均温度的模型

(6)

式中，Ts表示地表处的大气温度。

为了检验新的Tm模型的精度，文中将45004探空测站2017年的探空产品利用式(5)计算获得的Tm作为真值。此外，文中又分别评估了Bevis模型[16](Tm=70.2+0.72Ts)及Wang[17]等香港区域性模型(Tm=113.29+0.5863Ts)计算的Tm与探空资料估算Tm的差值，如表2所示。

表2 模型估算的Tm与真值的差值Tab.2 Bias of Tm between models-derived and radiosonde-derived K

表2显示了新的模型(NTm)估算的结果要明显优于Bevis模型结果(BTm)和Wang模型结果(WTm)。新的模型与真值之间的偏差的RMS值优于2.94K，偏差的变化幅度为[-3.70, 7.19]。

2 新的确定水汽层层顶的方法

层析网格划分的上下限的高度是从地面到对流层顶。实际上，从某个高度到对流层顶的区域中水汽密度接近于0[9]。由于其值太小，在实际的反演过程中得到的反演结果通常出现负值或异常值。为了提高对流层层析的精度，有学者定义了某个低于对流层顶的高度作为层析的上限高，并命名为水汽层层顶[9-10]，将其作为水汽反演的上限高度。为了精确地确定水汽层层顶的高度，本文提出了一种新的方法确定水汽层层顶。

在某个高度h上的大气可降水量PWVh可以表示为水汽密度的积分

(7)

其中，ρw表示液态水密度；ρk表示水汽密度。

则在垂直方向上，PWV呈现近似高斯指数函数分布

(8)

其中，PWV0为地表处的PWV；H为大气层等效高度。

通过对式(8)求导，便可得到不同高度处的斜率Kh

Kh=tan(α)

(9)

理论上，当Kh接近为0时，便对应为水汽层层顶的高度。为了确定水汽层层顶的高度，需要找到合适的Kh。本文定义当α=179.9°时对应的高度为水汽层层顶。通过香港地区2007—2016年6月份的掩星反演产品计算出香港地区水汽层层顶高度的变化，如图2所示。

图2 利用掩星产品定义的香港地区水汽层层顶高度随时间变化的规律Fig.2 Determining the water vapor top using radio occultation products in Hong Kong

从图2可以看出，香港地区6月份水汽层层顶变化范围为9.0km～10.6km。为了方便计算，文中求取水汽层层顶的平均值9.5km作为香港地区6月份的水汽层层顶。

3 试验分析

选取香港地区12个连续气象监测站2017年6月份的观测数据进行试算。利用Bernese 5.2软件进行处理获取天顶对流层延迟ZTD，基于式(3)对Sasstamoinen模型进行改正后得到ZHD，然后从ZTD中分离出ZWD。根据式(6)估算的大气加权平均温度代入式(4)中来获取湿度转换因子。选用Niell投影函数将ZWD投影到斜路径方向上SWD后，利用湿度转化因子将SWD转换为斜路径水汽含量SWV。把研究的香港地区在水平方向上划分5×5=25个网格，如图3所示。在垂直方向上，0.8km以下每400m划分一层，共2层；在0.8km～8.6km以上每600m划分一层，共13层；8.6km～9.5km划分一层，故垂直网格共划分成16个网格，则总网格为16×5×5=400个。选取Flores提出的经典的对流层层析方法，便可反演出水汽的三维分布信息。

图3 水平方向层析网格划分及测站的分布Fig.3 Horizontal distribution of GNSS stations and the division of voxels

根据上述的数据处理方式，每0.5h输出一次结果。限于篇幅，图4仅仅给出了第160d的处理结果。

图4 水汽的三维分布信息Fig.4 3D distribution of water vapor

由图4可以看出，2017年第160d的水汽主要集中在5km以下，而在5km以上水汽含量较少。随着时间的变化，0.8km～5km处的水汽密度发生较为明显的变化。为了检验优化后的层析结果，文中又利用传统的层析方法反演出对流层水汽密度值。将45004探测站的观测产品作为真值分别检验优化后的层析结果(Opti)和传统的层析结果(Trad)，如图5所示。

图5 层析结果与探空产品比较Fig.5 Statistical results between Opti & Trad and Radisonde

由图5可以看出，层析反演的水汽密度值与探空产品的变化趋于一致。当水汽密度在空间上存在逆增层的情况下，层析水汽技术不能很好地反演出这种情况；Opti结果要稍好于Trad。此外，从图5中还可以看出，Opti、Trad与Radiosonde在5km以上符合性较好，在5km以下差距较大。为了统计层析技术估算的水汽密度与探空水汽密度值之间的差距，文中分段分别统计了2017年6月份5km以下和5km以上的层析水汽值与探空水汽值的差值，如表3所示。

表3 Opti和Trad与Radiosonde偏差统计图Tab.3 The statistical results between Opti & Trad and Radisonde

由表5可以看出，在5km以下由于存在逆增层现象，故Opti、Trad与Radiosonde的偏差较大，而在5km以上偏差较小。Opti与Trad相比，在5km以下精度提高了1.24g/m3，在5km以上提高了0.22g/m3。

4 结论

利用无线电掩星产品辅助地基GNSS探测大气水汽的空间分布。针对对流层层析过程中3个关键技术进行优化：1)利用无线电掩星历史资料标定传统的Sasstamoinen干延迟模型，并建立了区域性大气加权平均温度模型，提高湿度转换因子的精度；2)精确地估算了水汽层层顶的高度，将其作为层析网格的上边界层。选取了香港地区12个监测站2017年6月份的观测数据进行试算，将反演的结果与探空产品进行比较和统计，验证了优化后层析方法的可行性和有效性。

致谢非常感谢香港特别行政区地政署测绘局提供的GNSS观测数据及气象观测数据。非常感谢UCAR提供的COSMIC和GRACE掩星观测产品。