2. 中国矿业大学环境与测绘学院, 江苏 徐州 221116;
3. 江苏师范大学地理测绘与城乡规划学院, 江苏 徐州 221116
2. School of Environment Science and Spatial Informatics, China University of Mining and Technology, Xuzhou 221116, China;
3. School of Geography, Geomatics and Planning, Jiangsu Normal University, Xuzhou 221116, China
水汽是大气的重要组成成分，具有极为复杂的时空变化特性，在全球气候变化以及极端气象灾害中都扮演着重要角色。地基GNSS水汽层析技术凭借其全天候监测、高时空分辨率、不受天气影响以及监测成本低等优势，已成为探测大气水汽高时空分辨率变化特性的重要手段之一[1]。
GNSS水汽层析技术是指利用区域地基CORS(continuously operating reference stations)网进行局域对流层三维水汽场的反演，由文献[2]首次提出，并重构出夏威夷地区的三维湿折射率分布。在传统GNSS水汽层析模型中，由于“盒形”层析区域与“倒锥形”GNSS信号簇的空间几何形态的不匹配性，导致三维层析模型中大量的体素块无法被GNSS信号线穿刺，进而引起层析观测方程的秩亏性，这是GNSS水汽层析领域的研究难点和热点[1]。针对该问题，众多学者提出了不同的优化模型与算法[2-19]。文献[2]充分考虑水汽的空间分布特性，构造了水平约束和垂直约束条件来对空白体素块进行约束处理。文献[3—4]分别基于高斯距离加权函数与图像平滑原理对水平约束条件优化，使得三维水汽分布更符合真实状态。文献[5—9]基于GPS、BDS、GLONASS、Galileo多系统的卫星观测信号用于三维水汽层析研究，通过增加GNSS观测信号，大大降低了空白体素块的数量，并提高层析结果的精度。文献[10—11]基于虚拟思想，构造出虚拟GNSS观测站和观测信号来增加观测信号数量，以此来提高层析模型体素块的穿刺率。文献[12—14]对GNSS层析模型的空间分辨率进行深入研究，通过优化水平分辨率与垂直分辨率使GNSS信号穿刺尽可能多的体素块，改善层析结果的质量。此外，融合外部气象数据也是改善层析秩亏问题的重要手段之一。文献[15]利用探空数据和大气红外探测仪(atmospheric infrared sounder, AIRS)获取的高精度水汽廓线作为先验信息来估计三维水汽分布。文献[16]根据COSMIC(constellation observing system for meteorology, ionosphere and climate)提供的水汽信息构造了基于两步重构法的GNSS层析模型，结合先验水汽观测值有效提高了层析解算结果的精度。文献[17—18]利用ECMWF(European Centre for Medium-Range Weather Forecasts)的再分析资料和多年的Radiosonde数据构造比例因子模型，用于精确估计侧边信号的水汽含量，大幅提高观测信号数量和层析结果质量。
近年来，随着各类遥感卫星的升空，大量的水汽遥感数据凭借其监测范围广，空间分辨率高的优势已被广泛用于大范围对流层水汽变化监测研究。文献[20]利用大范围的InSAR数据获取高分辨率的可降水量(precipitable water vapor, PWV)差分观测值，并通过试验证明InSAR数据用于大气监测的有效性与可靠性，这为高分辨率遥感数据用于GNSS水汽层析模型研究提供了可能。文献[21—22]利用上述的PWV差分观测值构造新型GNSS层析模型约束条件，提高GNSS层析结果的精度。但由于InSAR数据提供的是差分PWV信息，在与GNSS信号融合过程中需要将差分PWV转换为完整的PWV数据，一般借助于数值天气预报模型(numerical weather prediction, NWP)提供的全球水汽背景场进行转换[23]，这会降低高分辨率PWV观测值的精度。与InSAR数据类似的有MODIS、MERIS、FY-3A等水汽遥感数据，该类遥感数据可直接提供高分辨率的完整PWV观测值[23-24]，可被视为融合高分辨率遥感数据的GNSS水汽层析技术研究的理想数据源。
本文基于高分辨率MODIS PWV观测信息，充分挖掘其空间几何结构优势，构造出GNSS/MODIS紧耦合水汽层析算法。该算法首先利用GNSS PWV对MODIS PWV进行检验和校正，并分析了传统GNSS层析模型利用MODIS信号的缺点，提出基于体素节点参数化构造GNSS/MODIS紧耦合层析观测方程组的方法，可充分利用MODIS观测信号，并有效缓解GNSS层析模型的病态性问题。此外，利用徐州地区的探空数据对层析结果进行了精度评估，结果表明本文提出的紧耦合算法可以增强GNSS层析模型的稳定性，保证重构三维水汽场的准确性。
1 GNSS水汽层析原理GNSS信号的斜路径水汽总量(slant water vapor, SWV)分为各向同性与各向异性两部分。利用湿映射函数VFM1[25]和水平湿延迟梯度可计算出信号的斜路径湿延迟(slant wet delay, SWD)，再通过式(2)将其转换为SWV
(1)
(2)
式中，ZWD为天顶对流层湿延迟，可通过对流层总延迟(zenith total delay, ZTD)和天顶干延迟(zenith hydrostatic delay, ZHD)相减获得；
(3)
式中，ρw=1×103 kg/m3表示液态水密度；Rv=461.495 J/kg表示水汽气体常数；k3=3.776×105 K2/hPa, k′2=16.52 K/hPa, 两者均表示大气折射率常数；Tm表示大气加权平均温度[27]。
在GNSS层析模型中，通过对三维层析区域进行离散化可以获得均匀分布的体素块，一般用体素块中心位置处的水汽密度值代表该体素块的水汽密度，并且假设在一个层析历元内体素块的水汽密度为常数。因此卫星信号路径上的SWV可由式(4)表达
(4)
式中，SWV表示卫星信号路径上的水汽含量值，层析区域一共分成t层，每一层内按一定的经度差和纬度差均匀划分为n×m个体素块；dijk表示卫星信号穿过位于第k层，第i行，第j列的体素块的截距；xijk表示对应体素块的水汽密度。
2 MODIS水汽产品检验及校正中分辨率成像光谱仪(moderate resolution imaging spectroradiometer, MODIS)是美国发射的两颗环境遥感卫星Terra和Aqua上搭载的主要探测仪器之一，可提供36个光谱波段的地球观测信息，为自然灾害与生态环境监测、全球环境和气候变化研究提供强有力的数据支持。这两颗卫星两天可对全球进行一次观测，可满足突发性、快速变化的环境监测[28]。此外，凭借MODIS数据高时空分辨率的优势，其提供的高分辨率PWV数据已被广泛应用于InSAR数据的大气延迟改正[28]。但由于MODIS水汽产品具有一定的滞后期，且时间分辨率明显低于GNSS的全天候监测，其在实时大气水汽反演中存在一定的局限性。
2.1 MODIS水汽产品简介MODIS数据包括Terra卫星平台提供的MOD05_L2以及Aqua卫星提供的MYD05_L2两种水汽产品，包含每个像元的经纬度、扫描时间、太阳天顶角和方位角、传感器天顶角和方位角、近红外水汽值、红外水汽值、云掩膜产品等13个参数[24]。本文采用来自Terra平台的MOD05_L2水汽产品进行层析试验。
MOD05_L2产品数据的空间分辨率有两种：1 km×1 km和5 km×5 km。在InSAR大气延迟改正中，通常选用与SAR影像接近的1 km×1 km高分辨率水汽数据。在本文中，一方面考虑到层析模型的水平分辨率一般大于10 km×10 km[2-19]；另一方面，1 km×1 km分辨率的MODIS数据会成倍地增加GNSS层析模型的计算量，因此选用5 km×5 km分辨率的水汽产品用于层析试验研究。图 1(a)和(b)分别显示了2016年7月1日14:00时我国东部地区的MODIS水汽数据以及徐州地区的局部放大图，图 1(a)中的灰色区域表示被云朵污染的部分。
|
| 图 1 我国东部地区及徐州地区MODIS水汽分布图 Fig. 1 MODIS water vapor map of eastern China and Xuzhou, Jiangsu province |
2.2 MODIS PWV检验与校正
如图 1(a)所示，MODIS数据可以提供每个像素点的可降水汽含量，简记为MODIS PWV。由于MODIS信号易受云朵的影响，需要利用MODIS的云掩膜产品将受污染的MODIS PWV剔除，此外，还需利用高精度GNSS PWV对有效的MODIS PWV观测值进行检验和校正[29]。本文选用徐州地区的5个GNSS测站数据和对应时间上通过该地区的MODIS PWV数据进行比较分析。GNSS测站分布如图 1(b)所示，ZHJI、LGUO、BNTG、SGOU、SANP代表 5个测站的位置，RS 58027对应的是探空站的位置分布。同时，选取徐州地区2016年7月中15幅观测质量较好的MODIS影像进行云朵影响剔除，然后参与层析试验。表 1列出了15幅MODIS影像的获取时间和对应的层析时段范围。
| MODIS影像 | 获取日期 | 层析时段 |
| Image01 | 7-01 UTC 14:00 | 13:45—14:15 |
| Image02 | 7-06 UTC 14:20 | 14:05—14:35 |
| Image03 | 7-10 UTC 13:55 | 13:40—14:10 |
| Image04 | 7-11 UTC 14:35 | 14:20—14:50 |
| Image05 | 7-17 UTC 14:00 | 13:45—14:15 |
| Image06 | 7-21 UTC 02:35 | 02:20—02:50 |
| Image07 | 7-22 UTC 03:15 | 03:00—03:30 |
| Image08 | 7-23 UTC 02:20 | 02:05—02:35 |
| Image09 | 7-23 UTC 15:00 | 14:45—15:15 |
| Image10 | 7-24 UTC 03:05 | 02:50—03:20 |
| Image11 | 7-24 UTC 14:05 | 13:50—14:20 |
| Image12 | 7-25 UTC 14:50 | 14:35—15:05 |
| Image13 | 7-26 UTC 02:55 | 02:40—03:10 |
| Image14 | 7-27 UTC 14:35 | 14:20—14:50 |
| Image15 | 7-28 UTC 13:40 | 13:25—13:55 |
基于GNSS全天候运行的优势，可以通过处理GNSS数据获取与MODIS影像同一时刻的水汽信息，利用转换系数可以将式(1)中的ZWD转换为PWV，如式(5)所示
(5)
式中，ZWDGNSS为GNSS ZWD, ZWDGNSS和Ⅱ分别表示GNSS测站的天顶湿延迟与转换系数，同式(2)一致。
由图 1(b)可知，通常情况下，MODIS像素点与GNSS测站的位置并非完全重合，因此在对MODIS PWV进行检验以及校正的过程中，需要确定距离GNSS站点最接近的MODIS像元，作为待检验的MODIS PWV。图 2为徐州地区5个测站的GNSS PWV与MODIS PWV对比的散点图。
|
| 图 2 校正前后MODIS PWV与GNSS PWV对比 Fig. 2 Comparison between MODIS PWV and GNSS PWV before and after correction |
由图 2可以看出，GNSS PWV与MODIS PWV之间具有较好的一致性，其相关系数为0.736 4，这为利用GNSS PWV构造MODIS PWV校正模型提供了基础。根据文献[29]所述，本文所选用的MODIS影像范围较小，时间尺度较短，可利用GNSS PWV和MODIS PWV的实测数据，基于线性回归分析构造MODIS PWV的校正模型[29]，其表达式如下
(6)
式中，0.406为模型系数，31.77为常数项。表 2统计校正前后MODIS PWV与GNSS PWV的相关系数以及精度指标。
| 校正前后 | 统计量 | 相关系数 | RMS | STD | Bias |
| 校正前 | 平均值 | 7.24 | 3.14 | 2.74 | |
| 最小值 | 0.7364 | 3.30 | 0.45 | -18.77 | |
| 最大值 | 19.54 | 6.34 | 8.88 | ||
| 校正后 | 平均值 | 1.35 | 1.27 | -0.002 | |
| 最小值 | 0.9658 | 0.55 | 0.18 | -4.30 | |
| 最大值 | 3.58 | 2.57 | 3.25 | ||
由表 2可知，校正后的相关系数为0.965 8。RMS从7.24 mm降低到1.35 mm，表明校正效果较好, 可满足层析试验要求[22]。
3 GNSS/MODIS信号紧耦合层析算法校正后的高分辨率MODIS PWV观测值可被添加到GNSS水汽层析模型中，但在传统层析模型中通常假设一个体素块内水汽密度处处相等，受这一不合理假设的影响，MODIS PWV观测值在传统层析模型中只能作为观测约束，无法充分利用MODIS观测信号。本节基于体素节点模型对GNSS与MODIS两种观测信号进行参数化，提出GNSS/MODIS信号紧耦合层析算法，弥补传统层析模型在融合MODIS信号方面的不足，充分发挥MODIS观测信号的空间几何优势。
3.1 融合MODIS信号的传统层析模型一般情况下，PWV观测值可以看作沿垂直方向上的水汽密度积分值，同理，PWVMODIS可看作MODIS像素点位置处垂直方向上的水汽密度积分值，如式(7)所示
(7)
式中，ρ(s)表示垂直方向路径上的水汽密度值，单位为g/m3，s为垂直方向路径。在传统层析模型中(图 3)，式(7)可离散化为式(8)
|
| 图 3 传统层析模型融合MODIS信号 Fig. 3 Traditional tomographic model fusion with MODIS signal |
(8)
式中，xk表示第k层体素块的水汽密度值，单位为g/m3；dhk表示该体素块的高度，即层析模型第k层的厚度。
在传统层析模型中，受体素块内的水汽密度处处相等这一假设的影响，一列体素块的水汽密度在垂直方向上的积分对应于一个PWV观测值。由图 1可知，MODIS PWV数据的空间分辨率为5 km×5 km，明显高于GNSS层析模型的水平分辨率，因此会导致多个MODIS PWV观测值集中分布在一列体素块中，如图 3所示。若要将这些观测值融合到GNSS层析模型中，需要根据一列体素块内的所有MODIS PWV值计算出体素块中心位置处的PWVcenter，作为该列体素块的PWV观测信息。考虑到大气水汽分布具有一定的空间结构性，本文采用Kriging插值来估计PWVcenter值，利用空间随机场的结构性信息降低估计值的不准确性[30]，计算如式(9)
(9)
式中，λi表示加权系数，即Kriging系数，选用指数模型的变异函数来计算该系数[30]。PWVMODISi表示第i条MODIS信号的PWV观测值。利用PWVcenter可以构造出其所在列体素块的层析观测方程
(10)
式中，AGNSS、AMODIS、AH以及AV分别表示GNSS信号观测方程、MODIS信号观测方程、水平约束方程以及垂直约束方程的系数矩阵，其中AGNSS、AMODIS表示两种信号在体素块层析模型中的截距信息；X为水汽密度的未知参数向量。
分析式(10)可知，在MODIS观测条件良好情况下，MODIS信号所对应的层析观测方程的数量为m×n(m、n分别表示GNSS层析模型经度方向和纬度方向上的体素块个数)。但是，由图 1可知，MODIS观测信号极易受云朵影响，导致部分PWVMODIS观测信号不可用。因此，在传统体素块层析模型中无法发挥MODIS信号空间几何优势，造成有效观测信息的浪费。
3.2 融合MODIS信号的体素节点模型在传统层析模型中，受不合理假设影响，MODIS PWV信号的利用率大幅地降低。为有效利用MODIS信号，缓解层析方程的病态性，本文对每一条MODIS信号构造层析观测方程，以凸显其几何优势。考虑到同一体素块内的多条MODIS信号到体素块8个顶点的距离不同，采用文献[19]提出的Trilinear parameterization方法，可以为每条MODIS PWV观测值构造层析观测方程[19]。
在该方法中，将体素块8个顶点所处位置的水汽密度设为未知参数，并用8个顶点的水汽密度的平均值表示体素块中心的水汽密度值，以此来对传统层析模型进行优化。在本文中，考虑到大气水汽在垂直方向分布极不均匀这一特性，垂直方向采用指数加权方法来确定中心处的水汽密度值，记为改进Trilinear parameterization方法。如图 4所示，体素块中心O0处的水汽密度ρ(O0)可由式(11)计算出
|
| 图 4 改进的trilinear parameterization方法 Fig. 4 The improved tilinear parameterization method |
(11)
式中，ρ(O1)和ρ(O2)分别表示上平面中心O1和下平面中心O2的水汽密度值；H为水汽标高，取值为2 km[26]；dhk与式(8)一致。则体素中心O0的水汽密度值可用式(12)表达
(12)
式中，ρ(Pi)表示图 4中第i个顶点处的水汽密度。
本文基于改进参数化方法分别对GNSS SWV和MODIS PWV观测信号进行参数化。首先介绍基于体素节点模型GNSS信号的参数代过程，传统GNSS层析观测方程组的参数向量X可用新的参数向量Z代替
(13)
式中，m、n、k分别表示传统层析模型中经度、纬度、高度3个方向上体素块的数量；M为转换矩阵，其详细计算方式如式(14)所示
(14)
式中，mij为转换矩阵M中第i行，第j列的元素，其中i=a×b×c，a、b、c分别表示传统层析模型中体素块所处的行、列和层数。基于改进Trilinear parameterization方法，GNSS观测信号的层析方程可写成
(15)
式中，SWVq, 1表示GNSS信号斜路径水汽值组成的列向量。下面结合图 5介绍MODIS信号的节点参数化过程。
|
| 图 5 基于体素节点模型的MODIS信号参数化过程 Fig. 5 The parameterization process of MODIS signal based on the voxel node tomography model |
图 5以1条MODIS PWV观测信号，右侧对应的是其中一个体素块的局部放大图。由图 1可知，该信号未穿过体素块的中心位置，其PWV值用式(16)表达
(16)
式中，ρ(s)和s与式(7)一致，ρ(sk)表示第k层的垂直积分路径sk上的水汽密度值。对位于第k层的分段MODIS信号而言，如图 5所示，其垂直路径的PWV积分值可利用数值积分中的梯形公式(17)来计算[31]
(17)
式中，ρ(S0)、ρ(S1)、ρ(S2)分别表示S0、S1、S2 3点处的水汽密度值。因此，要构造MODIS PWV观测值所对应的层析观测方程，则需计算出S0、S1、S2 3点处水汽密度值与8个顶点水汽密度的数量关系。
对于位于上、下层面的S1、S2来说，设其仅与同一层平面内的4个节点的水汽密度有关，以S1为例，根据大气水平的空间分布特性，采用高斯加权函数[3]可建立如下数量关系
(18)
式中，
对于S0而言，可以用式(11)计算其与上、下层面的S1、S2处的水汽密度的数量关系，然后利用式(18)获得其与周围8个顶点的空间位置关系。因此，可构造第k层的体素块中MODIS PWV与8个顶点的水汽密度的函数关系为
(19)
根据式(19)可构造一条完整的MODIS PWV所对应的层析观测方程，与传统层析模型的融合过程相比，每一条MODIS信号可构造一条唯一的层析观测方程。因此，有多条MODIS信号构成的MODIS层析观测方程组如式(20)所示。此外，对于受云朵影响的MODIS PWV观测值可直接舍去，并不会对MODIS层析观测方程组造成较大的影响。
(20)
式中，
(21)
式中，BGNSS、BMODIS表示GNSS信号与MODIS信号所对应的系数矩阵，不同于式(10)中的AGNSS和AMODIS，这两个系数矩阵表示的是观测信号与节点之间的空间位置关系。BH和BV分别为体素节点模型中的水平约束[3]和垂直约束[2]的系数矩阵，其构造方法与传统体素模型相似，根据节点在水平和垂直方向的位置关系分别构造约束矩阵，Z为所有体素节点的水汽密度参数。
4 层析试验与分析 4.1 层析方案设置本文选取徐州地区2016年7月共15幅的MODIS影像以及相应时段的GNSS观测数据进行试验。利用GAMIT/GLOBK10.6软件对30 s采样率的GNSS数据进行处理，ZTD和梯度项的时间分辨率为5 min。利用Saastamoinen模型[24]估计ZHD值，选用VMF1湿映射函数[25]计算SWV值。如图 1(b)所示，层析试验区域范围为经度116.94°E—117.72°E，纬度33.94°N—34.64°N。在水平分辨率方面，经度方向和纬度方向的分辨率为0.13°×0.14°；垂直分辨率方面，采用下密上疏的非等间距划分方法[13]，共划分为15层。并且根据徐州地区10年(2009—2018年)每年7月探空数据提供的水汽廓线，确定最优对流层顶高度为11 km，考虑到不同天气条件的影响，图 6展示了2009—2018年徐州地区每年7月的水汽垂直分布特征。此外，本文采用代数重构算法解算层析观测方程组，利用层析时段前3天的探空垂直廓线信息作为迭代初值进行层析试验[4]。
|
| 图 6 徐州地区2009—2018年每年7月水汽垂直分布图 Fig. 6 Vertical distribution of water vapor in Xuzhou during each July from 2009 to 2018 |
为了验证本文提出的紧耦合算法，共设计3种试验方案来进行对比分析。
方案1：采用传统层析模型，仅利用GNSS数据作为观测信息。
方案2：采用传统层析模型，并利用GNSS和MODIS数据构造层析观测方程进行解算。
方案3：采用体素节点模型，利用改进Trilinear parameterization方法分别对GNSS信号与MODIS信号进行参数化，并利用两者构成的层析观测方程组进行层析解算。
4.2 MODIS PWV信号贡献性分析在体素节点模型中加入MODIS信号后，部分未被GNSS信号线穿刺的体素块将会被MODIS信号线穿过，可以有效弥补传统层析模型几何结构的缺陷。图 7展示了加入MODIS信号前后，层析模型中第1层、第5层、第10层、第15层的体素块所穿过信号的平均数量对比。
|
| 图 7 加入MODIS信号前后穿过体素块的信号数量对比 Fig. 7 The comparison of the number of observation signals passing through the voxels before and after adding MODIS signals |
由图 7可知，一方面，在0~2 km的近地层，传统层析模型中有大量的体素块未被GNSS信号穿过，尤其第1层中只有GNSS测站所在的体素块被观测信号穿刺，加入MODIS信号后，这些空白体素块均有信号线通过，较好地改善了层析模型的观测几何结构。另一方面，在4 km以上的高层中，加入MODIS信号之后同样使得每个体素块所穿过的信号数量增加4~6条，为层析模型引入更多的有用观测信息。图 8统计了每个层析时段内有效的GNSS信号和MODIS信号的数量，将MODIS信号添加到层析系统后，平均有效观测信号数量从366条增加到491条，提高了34.15%，这表明引入MODIS数据可以提高有效观测值数量。
|
| 图 8 层析时段内不同类型观测信号的数量对比 Fig. 8 Comparison of the number of different observation signals during the tomography period |
4.3 层析结果精度分析
在层析结果精度评定中，通常以探空站在UTC 00:00和UTC 12:00两个时刻提供的高精度垂直水汽廓线信息为参考值，但本文中的层析时段无法与探空观测时间完全一致，因此选用距层析时段最近的探空数据对水汽反演结果进行检验[21-22]。图 9(a)和图 9(b)分别给出了层析时段内3种方案的层析结果的均方根误差(root mean squared error，RMSE)和平均绝对误差(mean square error, MAE)的对比情况。由图 9可知，在大多数层析历元内，方案3的RMSE和MAE值都要优于方案1和方案2，方案2的统计量略优于方案1。此外，表 3列出了15个层析时段内不同方案的层析结果与探空数据对比的统计结果，可知3种方案的层析结果的平均RMSE分别为2.43、2.24、1.82 g/m3；平均MAE分别为1.96、1.72、1.43 g/m3。相对于传统层析模型(方案1)，本文提出的算法的反演精度(RMS和MAE)平均提高了26.07%，表明融合MODIS信号的体素节点模型具有更好的反演能力。
|
| 图 9 层析时段内3种方案的层析结果的RMSE和MAE对比 Fig. 9 Comparison of the RMSE and MAE of the tomographic results derived from three schemes during the tomography periods |
| 试验方案 | RMS | MAE | STD | bias |
| 方案1 | 2.43 | 1.96 | 2.41 | -0.30 |
| 方案2 | 2.24 | 1.72 | 2.23 | -0.21 |
| 方案3 | 1.82 | 1.43 | 1.80 | -0.14 |
为了进一步比较不同算法反演水汽的时空分布情况，对两种天气条件下(晴天和雨天)反演的水汽廓线进行对比分析。其中，7-01、7-11、7-21、7-25为晴天；7-17、7-26、7-27、7-28对应的是雨天。图 11给出了两种天气条件下3种方案反演得到的水汽廓线与探空站提供的水汽廓线的对比图。对比来看，方案1和方案2得到的水汽廓线较为接近，说明将MODIS信号融合到体素块层析模型中并不能起到有效作用，但方案3的水汽廓线明显优于两者，并与探空站反演的水汽廓线相一致。
|
| 图 10 3种方案在不同天气条件下的层析水汽廓线对比 Fig. 10 Comparison of tomographic water vapor profiles derived from three schemes under different weather conditions |
此外，根据水汽廓线比对情况，本节对层析时段内不同高度层的水汽密度进行精度评定，图 11给出了3种方案的层析结果与探空站所在列体素块的水汽密度的误差对比图。可以看出，在4 km以上的高度层中，优化算法的层析结果质量优于传统算法；在2~4 km的高度层中，方案3的反演精度略优于前两种方案；0~2 km的近地层中，富含60%的大气水汽含量，是造成信号延迟的主要区域。通过定量计算可得3种算法在0~2 km的层析结果平均RMSE分别为3.87 g/m3、3.75 g/m3、2.55 g/m3，可以看出优化算法的反演精度明显高于传统算法，进一步表明融合MODIS观测信号的可改善近地层三维水汽场的重构质量。
|
| 图 11 层析时段内3种方案层析结果与探空参考值的误差对比 Fig. 11 The difference between tomographic results and radiosonde along the radiosonde vertical columns |
5 结论
本文基于高分辨率MODIS PWV数据，提出GNSS/MODIS信号紧耦合水汽层析算法。凭借体素节点模型的建模优势，分别建立基于GNSS信号和MODIS信号的层析观测方程组。利用2016年7月徐州地区15幅MODIS影像和同一时刻下的GNSS实测数据，系统地比较分析了紧耦合算法与传统算法的层析结果的质量及精度。主要结论如下：
(1) 引入MODIS观测信号后，平均有效观测信号数量提高了34.15%，且所有体素块中均有观测信号穿过，可有效改善传统层析模型的GNSS信号几何结构缺陷。
(2) 与传统算法的层析结果相比，融合MODIS信号的传统算法和所提算法层析结果的平均RMSE值由2.43 g/m3分别降低为2.24 g/m3和1.82 g/m3，平均反演精度分别提高了7.81%和25.10%；与融合MODIS信号的传统算法相比，紧耦合算法的反演精度提高18.75%。
(3) 在0~2 km的近地层，紧耦合算法反演的三维水汽场质量明显优于其他地方两种层析结果，该高度层的平均RMSE值由3.87 g/m3(方案1)和3.75 g/m3(方案2)降低为2.55 g/m3(方案3)，反演精度分别提高了34.11%和32%，这表明紧耦合算法可最大化利用MODIS观测信号，有效改善层析结果质量。
致谢: 特别感谢徐州市自然资源和规划局提供的XZCORS网监测数据，美国国家航空航天局提供的MODIS影像，美国怀俄明大学提供的长期探空数据，麻省理工学院提供的GAMIT/GLOBK软件。
| [1] |
姚宜斌, 张顺, 孔建. GNSS空间环境学研究进展和展望[J]. 测绘学报, 2017, 46(10): 1408-1420. YAO Yibin, ZHANG Shun, KONG Jian. Research progress and prospect of GNSS space environment science[J]. Acta Geodaetica et Cartographica Sinica, 2017, 46(10): 1408-1420. DOI:10.11947/j.AGCS.2017.20170333 |
| [2] |
FLORES A, RUFFINI G, RIUS A. 4D tropospheric tomography using GPS slant wet delays[J]. Annales Geophysicae, 2000, 18(2): 223-234. DOI:10.1007/s00585-000-0223-7 |
| [3] |
SONG Shuli, ZHU Wenyao, DING Jincai, et al. 3D water-vapor tomography with Shanghai GPS network to improve forecasted moisture field[J]. Chinese Science Bulletin, 2006, 51(5): 607-614. DOI:10.1007/s11434-006-0607-5 |
| [4] |
何林, 柳林涛, 苏晓庆, 等. 水汽层析代数重构算法[J]. 测绘学报, 2015, 44(1): 32-38. HE Lin, LIU Lintao, SU Xiaoqing, et al. Algebraic reconstruction algorithm of vapor tomography[J]. Acta Geodaetica et Cartographica Sinica, 2015, 44(1): 32-38. DOI:10.11947/j.AGCS.2015.20130308 |
| [5] |
ZHAO Qingzhi, YAO Yibin, CAO Xinyun, et al. An optimal tropospheric tomography method based on the multi-GNSS observations[J]. Remote Sensing, 2018, 10(2): 234. DOI:10.3390/rs10020234 |
| [6] |
DONG Zhounan, JIN Shuanggen. 3-D water vapor tomography in Wuhan from GPS, BDS and GLONASS Observations[J]. Remote Sensing, 2018, 10(1): 62. |
| [7] |
BENEVIDES P, NICO G, CATALÄO J, et al. Analysis of Galileo and GPS integration for GNSS tomography[J]. IEEE Transactions on Geoscience and Remote Sensing, 2017, 55(4): 1936-1943. DOI:10.1109/TGRS.2016.2631449 |
| [8] |
夏朋飞, 叶世榕, 江鹏. GPS/GLONASS组合精密单点定位技术在三维水汽层析中的应用[J]. 大地测量与地球动力学, 2015, 35(1): 72-76. XIA Pengfei, YE Shirong, JIANG Peng. Research on three-dimensional water vapor tomography using GPS/GLONASS PPP method[J]. Journal of Geodesy and Geodynamics, 2015, 35(1): 72-76. |
| [9] |
胡鹏, 黄观文, 张勤, 等. 顾及边界入射信号的多模水汽层析方法[J]. 测绘学报, 2020, 49(5): 557-568. HU Peng, HUANG Guanwen, ZHANG Qin, et al. A multi-GNSS water vapor tomography method considering boundary incident signals[J]. Acta Geodaetica et Cartographica Sinica, 2020, 49(5): 557-568. DOI:10.11947/j.AGCS.2020.20190113 |
| [10] |
姚宜斌, 赵庆志, 罗亦泳. 附加虚拟信号精化水汽层析模型的方法[J]. 武汉大学学报(信息科学版), 2017, 42(11): 1658-1664. YAO Yibin, ZHAO Qingzhi, LUO Yiyong. An approach of imposing virtual signals to sophisticate water vapor tomographic model[J]. Geomatics and Information Science of Wuhan University, 2017, 42(11): 1658-1664. |
| [11] |
ADAVI Z, MASHHADI-HOSSAINALI M. 4D tomographic reconstruction of the tropospheric wet refractivity using the concept of virtual reference station, case study: northwest of Iran[J]. Meteorology and Atmospheric Physics, 2014, 126(3-4): 193-205. DOI:10.1007/s00703-014-0342-4 |
| [12] |
CHEN Biyan, LIU Zhizhao. Voxel-optimized regional water vapor tomography and comparison with radiosonde and numerical weather model[J]. Journal of Geodesy, 2014, 88(7): 691-703. DOI:10.1007/s00190-014-0715-y |
| [13] |
XIA Pengfei, YE Shirong, JIANG Peng, et al. Assessing water vapor tomography in Hong Kong with improved vertical and horizontal constraints[J]. Annales Geophysicae Discussions, 2018, 7(9): 1-28. |
| [14] |
ROHM W, ZHANG K, BOSY J. Limited constraint, robust Kalman filtering for GNSS troposphere tomography[J]. Atmospheric Measurement Techniques, 2014, 7(5): 1475-1486. DOI:10.5194/amt-7-1475-2014 |
| [15] |
BENEVIDES P, CATALAO J, NICO G, et al. 4D wet refractivity estimation in the atmosphere using GNSS tomography initialized by radiosonde and AIRS measurements: results from a 1-week intensive campaign[J]. GPS Solutions, 2018, 22(4): 91. DOI:10.1007/s10291-018-0755-5 |
| [16] |
XIA Pengfei, CAI Changsheng, LIU Zhizhao, et al. GNSS troposphere tomography based on two-step reconstructions using GPS observations and COSMIC profiles[J]. Annales Geophysicae, 2013, 31(10): 1805-1815. DOI:10.5194/angeo-31-1805-2013 |
| [17] |
赵庆志, 姚宜斌, 姚顽强, 等. 利用ECMWF改善射线利用率的三维水汽层析算法[J]. 测绘学报, 2018, 47(9): 1179-1187. ZHAO Qingzhi, YAO Yibin, YAO Wanqiang, et al. A method to improve the utilization rate of satellite rays for three-dimensional water vapor tomography using the ECMWF data[J]. Acta Geodaetica et Cartographica Sinica, 2018, 47(9): 1179-1187. |
| [18] |
ZHANG Wenyuan, ZHANG Shubi, DING Nan, et al. A tropospheric tomography method with a novel height factor model including two parts: isotropic and anisotropic height factors[J]. Remote Sensing, 2020, 12(11): 1848. DOI:10.3390/rs12111848 |
| [19] |
PERLER D, GEIGER A, HURTER F. 4D GPS water vapor tomography: new parameterized approaches[J]. Journal of Geodesy, 2011, 85(8): 539-550. DOI:10.1007/s00190-011-0454-2 |
| [20] |
MIRANDA P M A, MATEUS P, NICO G, et al. InSAR meteorology: High-resolution geodetic data can increase atmospheric predictability[J]. Geophysical Research Letters, 2019, 46(5): 2949-2955. DOI:10.1029/2018GL081336 |
| [21] |
BENEVIDES P, NICO G, CATALÃO J, et al. Bridging InSAR and GPS tomography: a new differential geometrical constraint[J]. IEEE Transactions on Geoscience and Remote Sensing, 2016, 54(2): 697-702. DOI:10.1109/TGRS.2015.2463263 |
| [22] |
HEUBLEIN M, ALSHAWAF F, ERDNVß B, et al. Compressive sensing reconstruction of 3D wet refractivity based on GNSS and InSAR observations[J]. Journal of Geodesy, 2019, 93(2): 197-217. DOI:10.1007/s00190-018-1152-0 |
| [23] |
TRZCINA E, ROHM W. Estimation of 3D wet refractivity by tomography, combining GNSS and NWP data: First results from assimilation of wet refractivity into NWP[J]. Quarterly Journal of the Royal Meteorological Society, 2019, 145(720): 1034-1051. DOI:10.1002/qj.3475 |
| [24] |
SAASTAMOINEN J. Contributions to the theory of atmospheric refraction[J]. Bulletin Géodésique, 1972, 105(1): 279-298. DOI:10.1007%2FBF02522083 |
| [25] |
BOEHM J, SCHUH H. Vienna mapping functions in VLBI analyses[J]. Geophysical Research Letters, 2004, 31(1): L01603. |
| [26] |
BEVIS M, BUSINGE S, HERRING T A, et al. GPS meteorology: remote sensing of atmospheric water vapor using the global positioning system[J]. Journal of Geophysical Research: Atmospheres, 1992, 97(D14): 15787-15801. DOI:10.1029/92JD01517 |
| [27] |
YAO Yibin, SUN Zhangyu, XU Chaoqian. Applicability of bevis formula at different height levels and global weighted mean temperature model based on near-earth atmospheric temperature[J]. Journal of Geodesy and Geoinformation Science, 2020, 3(1): 1-11. DOI:10.11947/j.JGGS.2020.0101 |
| [28] |
杨成生. 差分干涉雷达测量技术中水汽延迟改正方法研究[D]. 西安: 长安大学, 2011. YANG Chengsheng. Study on atmosphere delay correction methods in differential interferometric synthetic aperture radar[D]. Xi'an: Chang'an University, 2011. |
| [29] |
刘备, 王勇, 娄泽生, 等. CMONOC观测约束下的中国大陆地区MODIS PWV校正[J]. 测绘学报, 2019, 48(10): 1207-1215. LIU Bei, WANG Yong, LOU Zesheng, et al. The MODIS PWV correction based on CMONOC in Chinese mainland[J]. Acta Geodaetica et Cartographica Sinica, 2019, 48(10): 1207-1215. DOI:10.11947/j.AGCS.2019.20180386 |
| [30] |
何林. 地基GNSS水汽探测及其应用研究[D]. 北京: 中国科学院大学, 2013. HE Lin. Study on the ground GNSS vapor dentection and its application[D]. Beijing: University of Chinese Academy of Sciences, 2013. |
| [31] |
胡建华, 陈兴同, 曹德欣. 数值计算方法[M]. 徐州: 中国矿业大学出版社, 2008. HU Jianhua, CHEN Xingtong, CAO Dexin. Numerical calculation method[M]. Xuzhou: China University of Mining and Technology Press, 2008. |