联合GNSS/RS多源数据反演三维大气水汽分布研究

张文渊， 张书毕， 郑南山， 张秋昭*， 丁楠

1 中国矿业大学自然资源部国土环境与灾害监测重点实验室， 江苏徐州 221116 2 中国矿业大学环境与测绘学院， 江苏徐州 221116 3 苏黎世联邦理工学院大地测量与摄影测量研究所， 苏黎世 8093 4 江苏师范大学地理测绘与城乡规划学院， 江苏徐州 221116

0 引言

水汽是地球大气的重要组成部分，是引起天气变化的重要因素之一.水汽的相变会吸收或释放巨大热量，进而引发台风、暴雨、干旱、洪水等破坏性的天气现象.因此，准确地获取大气水汽的时空信息对于改善天气预报精度和灾害性气象预警具有重要作用(Zhang et al., 2015; Wang et al., 2018; He et al., 2019).但由于水汽时空变化极为复杂，传统探测手段受观测成本、采样间隔、测站分布等因素影响，难以实现对大气水汽的实时高分辨率监测(Zhang et al., 2015；姚宜斌等，2017b).Flores等(2000)首次提出GPS (Global Positioning System)水汽层析技术，利用夏威夷地区的GPS站网反演出局域三维湿折射率场.依靠其全天候监测、高时空分辨率、低成本等优势，GNSS(Global Navigation Satellite System)层析技术在过去20年得到了快速发展(Champollion et al., 2005; 夏朋飞等, 2015; Benevides et al., 2017; 赵庆志等, 2018; Trzcina and Rohm, 2019; 张文渊等, 2021).三维湿折射率(wet refractivity, WR)和水汽密度(water vapor density, WVD)是层析技术的主要产品，目前已应用于改善数值预报模式(Trzcina and Rohm, 2019)，研究台风周期变化(Zhu et al., 2020)，校正大气延迟误差(Yu et al., 2018)等领域.近年来，随着GPS、BDS (BeiDou Navigation Satellite System)、GLONASS (GLObal NAvigation Satellite System)、Galileo (Galileo Satellite Navigation System)四大系统的发展和现代化，传统GPS层析技术已发展为Multi-GNSS层析技术，利用多系统的GNSS观测信号，可有效增加观测值数量并改善GNSS信号的观测几何缺陷，提高层析结果的精度和可靠性(Benevides et al., 2017；Zhao et al., 2018；夏朋飞等，2015；胡鹏等，2020).

然而，受GNSS卫星和测站空间分布的限制，GNSS信号穿刺的“倒锥形”观测区域与“盒形”层析区域存在明显的不匹配性，尤其是在近地层区域，导致大量的层析体素块无法被GNSS信号线穿过，引起层析方程组的不适定性问题(Rohm, 2013)，这是GNSS水汽层析领域的研究难点和热点(姚宜斌等，2017b).针对该问题，Flores等(2000)添加水平约束和垂直约束来估计空白体素块的水汽参数；Song等(2006)利用高斯加权平均函数确定水平约束中周围体素块的权重，使得层析结果趋于真实的大气水汽分布；何林等(2015)在水平方向采用均值滤波器进行平滑处理，利用二阶Laplacian算子确定周围不同体素块的权重信息；基于此，Zhang 等(2017)利用Helmert方差分量估计建立了自适应Laplacian平滑方法确定最优的权重，改进后的层析结果与探空数据更为吻合.在层析模型优化方面，Chen 和 Liu(2014)利用GNSS站网的分布特点，根据实验确定最佳的层析模型水平分辨率，使得被穿刺的体素块数量达到最大值；随后，Rohm 等(2014)和Zhang 等(2015)共同提出内外双分辨率水汽层析方法，来降低层析区域边界处空白体素块的数量.Haji-Aghajany等(2020)利用WRF(Weather Research and Forecasting)模型获得初始体素块的先验水汽参数，然后构造新的优化体素块，大大地减少了体素块的数量.姚宜斌等(2017a)和Brenot 等(2020)分别提出两种附加虚拟观测信号的精化水汽层析模型，通过增加穿过研究区域的射线条数和射线穿过的体素块数，使层析结果更加逼近真值.

此外，融合外部气象数据也是提高层析反演质量的重要手段之一.Xia 等(2013)基于GPS观测值和COSMIC(constellation observing system for meteorology, ionosphere and climate)廓线信息建立了两步重构法的GNSS水汽层析模型.Benevides等(2018)利用大气红外探测仪(Atmospheric Infrared Sounder, AIRS)获取的高精度水汽廓线对GNSS层析模型进行初始化，有效提高了层析解算结果的精度.赵庆志等(2018)和Zhang等(2020)分别利用ECMWF(European Centre for Medium-Range Weather Forecasts)的再分析资料和30年的Radiosonde数据构造比例因子模型，用于精确估计侧边信号的水汽含量，大幅提高观测信号数量和层析结果质量.除此之外，高水平分辨率的遥感气象数据也可以提供丰富的大气水汽信息以及良好的几何观测结构，例如，MODIS(Moderate-Resolution Imaging Spectroradiometer)、MERIS(Medium Resolution Imaging Spectrometer)、以及我国发射的FY-3A卫星上的MERSI(Medium-Resolution Spectral Imager)传感器都可以提供1 km或5 km水平分辨率的PWV观测值(Chang et al., 2015; He and Liu, 2019)，这为融合RS和GNSS观测数据提供了基础.

本文以联合GNSS与RS多源观测的水汽数据为切入点，充分利用两种观测信号的空间互补优势，构造GNSS/RS联合水汽层析观测方程组，建立GNSS/RS联合水汽层析模型，有效解决了传统层析模型中观测几何限制的问题.利用徐州地区的GNSS实测数据和遥感水汽数据进行实验分析，并以研究区域内的无线电探空数据和ERA5再分析资料对层析结果进行精度评定，结果表明GNSS/RS层析模型可提高层析结果的精度，准确获取近地空间大气水汽的三维分布情况.

1 水汽层析原理

水汽层析技术是指利用观测信号在对流层中的湿延迟信息反演研究区域的三维湿折射率场或水汽密度场(Flores et al., 2000).在GNSS水汽层析中，观测信号的斜路径湿延迟(slant wet delay, SWD)可表示为湿折射率沿着信号方向从接收机到对流层顶的积分：

(1)

式中，s表示GNSS信号的传输路径，Nw为湿折射率.为反演研究区域的三维水汽密度分布，将SWD利用(2)式可转换为斜路径水汽含量(slant water vapor, SWV):

(2)

式中，105/RV·(k3/Tm+k′2)表示转换因子，记为Π.其中k′2=16.48 K·hPa-1，k3=3.75×105K2·hPa-1，两者均表示大气折射率常数，RV表示水汽气体常数，值为461.53 J/(kg·K) (Davis et al., 1985).Tm表示大气加权平均温度，本文利用Bevis模型计算Tm的值(Bevis et al., 1992)，该模型是目前应用最为广泛的Tm模型之一.根据水汽层析基本原理，将层析区域划分为均匀分布的体素块，则可建立每个体素块的水汽参数与卫星信号的SWV值的函数关系(Flores et al., 2000)：

(3)

式中，SWV单位为mm；dijk为斜路径射线穿过位于第k层，第i行，第j列的体素块的截距信息，单位为km；xijk表示层析模型的待求参数，即体素块的水汽密度值，单位为g·m-3.由于层析模型中大量体素块未被观测信号穿过，这些体素块所对应的截距信息为0，导致了层析观测方程的不适定性，通常采用水平和垂直经验函数对观测方程进行约束(Flores et al., 2000; Song et al., 2006; Benevides et al., 2017).

因此，联合层析观测方程、水平约束方程和垂直约束方程，可组成对流层层析方程组：

(4)

式中，SWVO和AO分别表示信号的SWV观测向量和对应的观测矩阵；AH和AV表示水平约束和垂直约束的系数矩阵，X为水汽密度的未知参数向量.

2 GNSS/RS联合水汽层析模型

2.1 GNSS/RS SWV联合观测值

由式(4)可知，精确的SWV观测值是反演大气水汽密度参数X的必要条件，GNSS/RS水汽层析模型的关键在于获取高精度的GNSS SWV和RS SWV观测值.图1给出了GNSS/RS联合水汽层析模型的具体流程，本节将详细介绍这两类关键观测值的计算过程.

图1 GNSS/RS联合水汽层析模型流程图Fig.1 Flow chart of the GNSS/RS combined water vapor tomography model

在GNSS观测中，测站可以接收到来自不同卫星的电磁波信号，利用双差网解法或非差精密单点定位方法可以对GNSS测站处的对流层延迟和梯度信息进行估计(范世杰, 2013)，然后利用(5)式可计算出观测信号的SWD值(Flores et al., 2000)：

SWD=mfw(ele)·ZWD+mfg(ele)

(5)

在遥感观测中，搭载在卫星上的传感器可以直接探测近地层的高分辨率水汽信息.例如，Terra和Aqua卫星上的MODIS传感器，Envisat卫星上的MERIS传感器，以及我国发射的FY-3A卫星上的MERSI传感器都可以提供高分辨率的PWV观测值，并且具有全球观测的优势(Chang et al., 2015; He et al., 2019).如图2所示，与GNSS SWV的计算方法相似，可以利用下式估计RS像元到遥感卫星的信号路径上的RS SWV观测值：

SWVRS(ε,α)=mw(ε)·PWVRS+Π·mfg(ε)

(6)

图2 反演RS SWV观测值的示意图Fig.2 Schematic of the retrieval of RS SWV observations

由式(6)可知，获取可靠的RS SWV观测值的关键是确定精确的RS PWV和水平梯度延迟量.但是，由于传感器的光谱信号易受云的影响，原始的RS PWV水汽产品存在观测误差.将其用于获取RS SWV观测值前，需要先对其进行校正(Li et al., 2003; Chang et al., 2015).另一方面，大气水汽的各向异性分布特点也会对斜路径湿延迟量造成一定影响，本文选用GPT3w(Global Pressure and Temperature 3 wet)模型获取大气水平梯度分布信息(Landskron and Böhm, 2017).接下来将详细论述RS PWV校正和梯度延迟估计的过程.

2.1.1 RS PWV校正

遥感水汽数据可以提供高分辨率的PWV观测影像，以MODIS数据为例，其包括MOD05_L2和MYD05_L2两种水汽产品，分别来自Terra和Aqua卫星，均可以提供1 km水平分辨率的近红外产品和5 km水平分辨率的红外产品，通常近红外数据只能在白天可以获取，而红外数据在白天和夜间均可以获取，因此5 km的PWV数据的时间分辨率要高于1 km的PWV数据.另外，GNSS层析模型的水平分辨率通常在10 km到50 km之间，若将1 km水平分辨率的PWV数据融合层析模型中，则会产生超大量的观测信号.尽管增加了观测值的数量，但是也会为层析方程解算带来一定的负担，而5 km的PWV数据既可以增加一定数量的观测值，又不会加重计算负担.综上，我们选用 5 km水平分辨率的PWV数据.同时，为了保持一致性，需要将1 km水平分辨率的MERSI PWV数据稀释成5 km水平分辨率的PWV产品.但由于遥感水汽产品观测精度不高，通常以观测区域内同步GNSS数据获取的GNSS PWV作为参考值，对RS PWV进行校正.

利用GNSS ZWD和式(2)中的比例因子∏的乘积，可以获得全天候监测的GNSS PWV观测值.大量实验表明与探空站提供的高精度PWV观测值相比，GNSS PWV的RMSE(root mean square error)为1～2 mm，可用于RS PWV的校正(Li et al., 2003; Chang et al., 2019).由于GNSS测站与RS像元未完全重合，一般利用反距离加权算法插值计算出GNSS测站处的RS PWV观测值，并将其作为待检验值，然后利用同一时刻的GNSS PWV对其进行校正(Chang et al., 2019; He et al., 2019).此外，考虑到大气水汽的快速变化特性，本文采用逐影像校正的方法依次对所选RS PWV观测影像构造线性回归校正模型，如下式所示(Li et al., 2003)：

(7)

2.1.2 大气水平梯度延迟

(8)

式中，A0表示平均值，A1和B1表示年周期系数，A2和B2表示半年周期系数，这5个系数都已事先计算好，以格网形式保存在文本文件中(姚宜斌等, 2015)，doy表示年积日，需要结合RS影像的UTC时刻换算成精确的儒略日时间信息再进行梯度估计.本文基于1°×1°的格网数据利用双线性内插得出RS影像获取时刻每个RS像元的两个梯度延迟量(朱明晨等, 2019).

2.2 GNSS/RS SWV联合层析方程组

基于校正后的RS PWV观测值和GPT3w模型估计的水平梯度延迟量，结合RS像元到卫星的信号的高度角和方位角信息，利用式(6)可以准确估计出RS信号的斜路径水汽含量(RS SWV)，因此融合GNSS SWV和RS SWV的层析观测方程如下所示：

(9)

式中，SWVGNSS和SWVRS分别表示GNSS SWV和RS SWV的观测向量，AGNSS和ARS则表示由GNSS和RS信号在三维层析模型中的截距信息组成的观测矩阵，AH和AV和与式(4)中一致.针对不同类型的观测信号，通常需要对其进行合理定权，参考Zhao 等(2017)和Heublein 等(2019)提出的定权方案，对GNSS和RS信号采用高度角定权，如式(10)所示；对于约束条件采用单位阵定权.

Pi=sin2(ei),

(10)

式中，ei和Pi分别表示第i条观测信号的高度角和权重信息.

代数重构算法(Algebraic Reconstruction Technique, ART)采用逐行迭代的方式对线性方程组进行解算，避免系数矩阵的求逆运算，具有稳定、高效的特点，在稀疏型方程解算中具有优势(Bender et al., 2011; 姚宜斌等, 2014; 张文渊等, 2021)，本文利用该算法对上述层析方程组进行解算.

(11)

3 层析试验分析

3.1 实验数据

本文的实验数据包括GNSS数据、MODIS和MERSI水汽数据、探空数据以及ERA5再分析数据.研究时间为2017年8月，徐州地区正值夏季，多雷暴雨天气.GNSS数据包括徐州Continuously Operating Reference Stations(CORS)网中10个均匀分布的GNSS测站提供的单系统GPS的观测数据，GNSS测站的水平分布如图3所示，其高程范围为31.2 m到60.7 m，高程变化对三维层析模型的影响可忽略不计.利用GAMIT10.7软件处理GNSS观测数据，并引入3个IGS站(BJFS站，SHAO站，LHAZ站)的观测数据，采用双差解算策略进行处理，以获取高精度的ZTD观测值(张文渊等, 2021).同时，采用Saastamoinen模型估计ZHD值(Saastamoinen, 1972)，选用VMF1湿映射函数计算GNSS SWV值(Böehm and Schuh, 2004).利用GNSS数据估计ZTD和梯度信息的时间分辨率为5 min.

关于遥感水汽数据，本文采用Aqua、Terra两个卫星提供的MODIS数据和FY-3A卫星提供的MERSI数据.由于近几年已不提供MERIS数据，所以文中没有使用该数据(He et al., 2019).根据层析区域的观测条件确定可用的MODIS和MERSI水汽数据，本文挑选置信度大于95%的RS PWV影像参与到层析实验中(Li et al., 2003)，表1列出了符合要求的40幅影像的观测卫星和观测时间.

图3 徐州地区GNSS测站分布示意图Fig.3 Distribution of GNSS stations in Xuzhou

表1 2017年8月徐州地区40幅RS PWV影像的观测卫星及获取时间Table 1 Observation satellites and acquisition time of the 40 RS PWV maps in Xuzhou in August 2017

此外，本文采用探空数据和ERA5数据对层析结果进行精度验证.探空数据由美国俄怀明大学(University of Wyoming)提供，其包括不同高度层的气压、温度、相对湿度、比湿等气象数据，可以精确地计算出不同高度处的水汽密度含量，常用于验证GNSS反演结果的精度(Li et al., 2003).ERA5数据是第五代大气再分析资料，可以提供0.25°×0.25°水平分辨率的37个气压层逐小时分辨率三维气象参数.相比于传统的ERA-Interim产品，ERA5数据具有更高的精度和时间分辨率(Cong et al., 2018).因此，我们分别采用高精度的探空垂直廓线和ERA5三维水汽分布作为验证数据来计算GNSS/RS层析模型的反演精度.

3.2 层析策略

在本文的层析实验中，层析区域为：经度116.52—117.96°E，纬度33.90—34.62°N.参考现有研究中有关层析格网划分的设定(Zhang et al., 2015, 2021; Benevides et al., 2018; Heublein et al., 2019; Haji-Aghajany et al., 2020)，考虑到徐州地区GNSS站网的稀疏分布特性，本文在经度方向和纬度方向的水平分辨率均设置为0.18°，对应的体素块个数为8个和4个.垂直方向采用非均匀划分的形式，共划分15层，1到10层每层间隔为0.4 km，11到15层的间隔分别为1.0 km、1.0 km、1.5 km、1.5 km、2.0 km，层析顶层高度为11 km，根据徐州地区10年的探空历史数据确定(张文渊等, 2021).单个层析历元的时长为30 min，并且假设在一个层析历元内体素块的水汽密度不发生变化(Benevides et al., 2018).本文层析时间段的设置主要根据遥感影像的观测时间而定，使得RS PWV影像观测时间恰好位于层析时段的中间时刻(例如，Image02的观测时间为UTC 18∶30，则该影像对应的层析历元为UTC 18∶15—18∶45)，保证GNSS数据与RS数据的时间一致性.由于两种遥感数据探测徐州地区的过境时间不一致，本文没有将两种遥感数据整合到一个层析历元中，这也是遥感观测的主要缺陷，其时间分辨率对于实时水汽反演有一定的限制.

本文设计如下两种实验方案来进行对比分析，由于Terra、Aqua和FY-3A卫星在徐州地区的过境时间并非每天的UTC 00∶00和UTC 12∶00，本文需要采用三次样条插值对探空数据进行时间域的插值(Chang et al., 2019)，以获得层析时刻的探空水汽廓线信息.然后利用探空插值数据和ERA5再分析资料分别对两种层析结果进行精度评定，计算层析结果的RMSE、STD(standard deviation)和MAE(mean absolute error)等精度指标，下边给出精度指标的计算公式：

方案一：只采用GNSS观测信号进行水汽反演，以“GNSS”表示；

方案二：采用GNSS与RS观测信号联合进行水汽反演，以“GNSS/RS”表示.

(12)

(13)

(14)

3.3 RS PWV校正分析

根据2.1.1节介绍，选用逐影像的校正方法依次利用GNSS PWV观测值对每幅遥感影像的PWV反演值进行校正，表2列出了校正前后多项精度指标的统计结果，可以看出校正后的RS PWV的精度有了明显提高，其平均RMSE值从8.06 mm降低到1.05 mm，与GNSS PWV的相关系数由0.56提高到0.93，可以满足水汽层析的要求(Heublein et al., 2019).

表2 校正前后RS PWV与GNSS PWV的相关系数以及精度统计量(mm)Table 2 Correlation coefficient and precision statistics of MODIS PWV and GNSS PWV before and after correction (mm)

3.4 层析模型观测几何分析

在传统GNSS水汽层析模型中，受GNSS卫星星座和测站分布的限制，位于层析模型低层区域的大量体素块无法被GNSS信号穿刺.但是，不同于GNSS信号的倒锥形观测形状，RS信号的观测几何为正锥形，两者恰好可以相互弥补，达到最优观测结构.图4对比了2017年8月22日UTC 03∶10—03∶40层析历元下，传统GNSS层析模型和GNSS/RS联合层析模型中观测信号的三维空间分布.此外，图5对比了该历元下两种层析模型中三维体素块的穿刺情况.

图4 2017年8月22日UTC 03∶10—03∶40传统GNSS层析模型(a)和GNSS/RS联合层析模型(b)信号分布图，蓝色线代表GNSS信号，红色线代表RS信号Fig.4 Signal distribution of traditional GNSS tomography model (a) and GNSS/RS combing tomography model (b) at UTC 03∶10—03∶40 on August 22, 2017. Blue lines represent GNSS signals while red lines represent RS signals

图5 传统GNSS层析模型(a)和GNSS/RS联合层析模型(b)的三维体素块穿刺对比Fig.5 Comparison of 3D crossed voxels between the traditional GNSS tomography model (a) and GNSS/RS tomography model (b)

由图5可知，在传统GNSS层析模型中，低层存在大量的空白体素块，而RS信号可以通过几乎所有的体素块.经统计可知，GNSS层析模型中1到5层的穿刺体素块数目为56个，占前5层总体素块数量的35%，GNSS/RS模型的穿刺体素块数目为160，所有体素块均被穿刺.进一步地，图6比较了各个层析历元内两种层析模型的观测信号数量以及穿刺体素块数量，其中平均观测信号数量由GNSS模型的479条增加到GNSS/RS模型的663条；相应地，平均穿刺体素块数量由297个提升到462个.综上可以看出加入RS信号可以充分穿刺层析模型的体素块，有效改善传统层析模型中观测信号与层析区域不匹配的缺陷.

图6 不同层析历元下两种模型的观测信号数量(a)以及穿刺体素块数量(b)对比Fig.6 Comparison of the number of observed signals (a) and the number of crossed voxels (b) for the two models during different tomographic epochs

3.5 层析结果精度分析

本节比较了GNSS层析模型和GNSS/RS层析模型反演的水汽分布情况，以探空水汽廓线和ERA5再分析资料为参考值，对两种模型获取的层析结果进行了精度评定.如上所述，由于时间的不一致性，本文引入三次样条插值算法对探空数据进行插值，用于确定层析时刻的探空站参考数据.尽管该算法已广泛应用于水汽产品的时间域插值处理，但不同的天气条件对探空插值结果的可靠性有一定的影响.参考Heublein等 (2019)中关于探空插值数据可靠性的检验方法，当满足以下两个条件时，探空插值廓线可以被视为可靠参考值：(1)在连续两个探空观测时刻上，GNSS PWV和探空PWV有很好的一致性；(2)探空站在层析时刻的插值PWV与同一时刻的GNSS PWV数据非常接近.在这两个条件中，以2 mm作为PWV误差的临界值.图7分别展示了2017年8月GNSS PWV与原始探空PWV的差值以及逐个层析时刻GNSS PWV与探空插值PWV的误差，可以看出历元1、3、4、12、14、15、16、17、22、25、26、27、28、29、35、36这16个层析时刻的探空插值数据符合要求.

图7 2017年8月GNSS PWV与探空PWV的误差(a)以及不同层析历元GNSS PWV与探空插值PWV的误差(b)Fig.7 (a) Difference between GNSS PWV and radiosonde PWV in August 2017 and (b) error of GNSS PWV and radiosonde interpolation PWV in different tomographic epochs

为了比较两种水汽层析模型的层析结果的质量，图8首先展示了不同天气条件下(晴天：2017年8月15日UTC 04∶40—05∶10，雨天：2017年8月18日UTC 22∶45—23∶15)两种层析模型反演的大气水汽廓线以及探空插值廓线(对应上述两个时间段).需要注意的是，由于探空站仅能提供其所在位置上空的垂直水汽廓线信息，因此这里只将探空站所处那一列的体素块的层析水汽廓线与探空水汽廓线进行对比.尽管通过验证分析，探空插值数据仍会存在一定的误差，不利用准确分析层析结果的精度，因此引入独立的ERA5数据进行精度验证，该数据还可以验证层析结果的整体三维分布质量，图9 比较了相同时间段内层析结果和ERA5水汽产品的三维分布情况.综合图8和9可以发现，在5 km(第11层)以上，由于两种模型均有大量观测信号穿过且分布较为均匀，对应的两种层析廓线较为接近，且与探空水汽廓线变化趋势一致.但在5 km以下，GNSS层析模型反演的水汽廓线明显低于GNSS/RS模型的水汽廓线，并且与探空水汽廓线也存在较大的差异.这可能是由于在GNSS层析模型中，“倒锥形”GNSS信号簇与“盒形”层析区域的空间几何形态的不匹配性导致GNSS层析模型中大量的体素块无法被GNSS信号线穿过.因此，在GNSS层析方程的解算过程中，这些空白体素块的水汽信息只能通过约束方程来进行估计，导致低层体素块的水汽结果低于GNSS/RS模型和ERA5数据.但在GNSS/RS层析模型中，通过引入“正锥形”RS信号观测值，使得低层的绝大多数体素块可以被信号穿过，保证三维层析模型的结构稳定性，进而提高了层析结果的质量.就层析模型的三维反演精度而言，同样可以看出GNSS/RS层析模型得到的水汽三维分布与ERA5资料提供的水汽三维分布更加接近，进一步表明遥感水汽信号对于提高三维水汽分布的反演精度有重要作用.

表3 两种方法的层析结果的统计量对比 (单位：g·m-3)Table 3 Statistics of the tomography results for two methods (unit: g·m-3)

图8 两种层析模型在不同天气条件下的水汽廓线对比(a) Rainless天气； (b) Rainy天气.Fig.8 Comparison of water vapor profiles for two tomography models under different weather conditions(a) Rainless weather; (b) Rainy weather.

图9 两种层析模型在不同天气条件下的三维水汽分布对比Fig.9 Comparison of 3D water vapor distribution for two tomography models under different weather conditions

为了更好的比较不同层析历元下两种模型的反演精度，图10给出了层析时段内两种模型的层析结果与探空数据和ERA5数据的RMSE对比.由图可知，在绝大多数层析历元内，GNSS/RS层析模型的RMSE都明显优于GNSS模型.此外，表3列出了不同模型的层析结果的精度统计量，可以看出以探空站为参考值，GNSS模型和GNSS/RS模型的平均RMSE分别为2.19 g·m-3、1.54 g·m-3，平均MAE分别为1.69 g·m-3、1.34 g·m-3，两者分别提高了29.68%和20.71%；以ERA5数据为参考值评估40个历元的层析结果可得两个模型的平均RMSE分别为2.56 g·m-3、2.04 g·m-3，平均MAE分别为2.03 g·m-3、1.63 g·m-3，分别提高了20.31%和21.18%.综合来看，相对于传统GNSS层析模型，本文提出的GNSS/RS层析模型的RMSE和MAE分别提高了25.00%和20.95%，这表明融合RS信号可以有效提高层析结果的精度.

图10 层析时段内两种模型的层析结果分别与有效探空插值数据(a)和ERA5数据(b)对比的RMSE值Fig.10 RMSE values of the tomographic results from the two models are compared with effective radiosonde interpolation data (a) and ERA5 data (b) respectively in the tomographic periods

为了进一步比较不同层析方法在垂直方向上的反演精度，分别以探空插值数据和ERA5数据为参考值，对不同高度层上的层析结果进行对比分析.图11给出了不同高度层上，两种层析水汽结果与探空数据和ERA5数据对比的均方根误差(RMSE)、平均偏差(mean deviation, MD)、标准偏差(STD)以及相对误差(relative error, RE)四个统计指标信息.可以发现，相对于两种参考数据，GNSS/RS联合模型的层析结果质量在绝大多数高度层上都优于传统GNSS模型，尤其是在0～2 km的近地范围内，改善效果最为明显.该范围内的层析结果的平均RMSE值由2.95 g·m-3降低到2.04 g·m-3，平均偏差和标准偏差分别减小了1.11 g·m-3和0.75 g·m-3，相对误差由17.45%改善为10.84%.这进一步证明了GNSS/RS层析方法在反演高精度三维水汽分布方面的优势.

图11 层析时段内两种模型的层析结果分别与探空插值数据(a—d)和ERA5数据(e—h)对比的分层精度信息Fig.11 Layered statistics of the tomographic results from the two models compared with radiosonde interpolation data (a—d) and ERA5 data (e—h) in the tomographic periods

4 结论

本文基于高水平分辨率RS信号的观测优势，首次创建了联合GNSS与RS多源观测的水汽层析模型，并详细论述了GNSS/RS层析模型的构造流程.利用2017年8月徐州地区的GNSS实测数据和MODIS、MERSI水汽遥感数据进行水汽层析实验，系统地比较分析了传统GNSS模型与GNSS/RS模型的观测几何结构和层析结果质量，主要结论如下：

(1)RS PWV作为获得准确RS SWV观测值的前提条件，在进行层析之前必须对其进行校正.利用同步的GNSS PWV对其进行校正后，RS PWV的平均RMSE由8.06 mm降低为1.05 mm，可以满足层析试验的要求.

(2)在层析模型的观测几何结构方面，GNSS/RS层析模型的平均观测信号数量和穿刺体素块数量分别为663条和462个，明显优于传统GNSS模型的479条和297个，两者分别提高了38.41%和55.56%.在GNSS/RS模型中，位于低层的多数体素块均能被RS信号较好地穿刺，有效改善传统层析模型的GNSS信号几何结构缺陷.

(3)分别以探空水汽廓线和ERA5三维水汽分布信息为参考值对两种模型的层析结果进行评估.

发现层析水汽廓线的平均RMSE从GNSS模型的2.19 g·m-3降低到GNSS/RS模型的1.54 g·m-3，精度提高了29.68%；就层析三维水汽分布情况而言，传统模型的RMSE为2.56 g·m-3，多源观测模型为2.04 g·m-3，三维分布精度改善了20.31%.此外，在不同高度层上，GNSS/RS联合模型的层析结果的精度统计量也优于传统的GNSS层析模型.结合本文可知高分辨率遥感数据在改善层析反演结果方面展现出较大的潜力，受其时间分辨率的影响，目前对于实时水汽反演存在一定的局限性.但是，遥感数据的大范围、高空间分辨率的优势，有助于进行大尺度水汽层析实验和三维水汽场精细化建模，这将是我们下一步的重点研究工作.

致谢作者特别感谢徐州市自然资源和规划局提供的XZCORS网监测数据，美国国家航空航天局提供的MODIS PWV影像，中国风云卫星遥感数据服务网提供的MERSI PWV影像，美国俄怀明大学提供的长期探空数据，麻省理工学院提供的GAMIT/GLOBK软件.