全球ZWD垂直剖面滑动窗口的格网模型

黄良珂，朱 葛，彭 华，陈 华，刘立龙，姜卫平

1. 桂林理工大学测绘地理信息学院，广西 桂林 541004； 2. 广西空间信息与测绘重点实验室，广西 桂林 541004； 3. 武汉大学测绘学院，湖北 武汉 430079； 4. 武汉大学卫星导航定位技术研究中心，湖北 武汉 430079

对流层延迟是影响全球导航卫星系统(GNSS)及其他空间技术高精度应用的关键因素。对流层天顶湿延迟(zenith wet delay,ZWD)与水汽变化密切相关，由于水汽在不同大气高度范围变化趋势不一致，因此将影响对流层延迟改正的精度。近年来，利用大气再分析资料(如NCEP、ERA-Interim、ERA5和MERRA-2等)积分计算的对流层延迟信息通过空间插值进行GNSS精密单点定位(precise point positioning,PPP)和GNSS大气水汽反演获得了广泛关注。其中，高精度的对流层延迟(ZWD/ZHD)初始值信息可加速精密单点定位的收敛时间，提升高程分量估计的精度[1-4]。由于大气再分析资料的格网点高度与GNSS站、探空站等用户所在高度不一致(排除不同数据源之间的高程基准差异)，这种高程差异在高海拔地区更为显著。ZWD在高程上的变化远大于其在水平方向上，因此若要将大气再分析资料格网ZWD数据精确插值到用户位置处，需依赖高精度ZWD垂直剖面模型对ZWD进行垂直插值。ZWD垂直剖面函数也是构建高精度ZWD模型的关键，因此研究实时高精度全球ZWD垂直剖面模型的构建具有重要的现实意义。

常用的ZWD垂直剖面模型主要分为两大类：需要实测气象参数的模型和非气象参数模型。在需要实测气象参数的对流层垂直剖面模型中，主要有Hopfield模型[5]、Saastamoinen模型[6]和Black模型[7]等。这些模型的ZWD垂直剖面函数依赖于气象参数，如果缺乏水汽垂直分布参数，计算得到的ZWD垂直改正精度有限。此外，由于实测气象参数不易实时地获取，限制了这些模型在高精度GNSS实时定位中的应用。针对上述模型的不足，诸多学者建立了不依赖于实测气象参数的对流层垂直剖面模型，如EGNOS模型[8]、UNB系列模型[9-10]和TropGrid系列模型[11]等。前两种模型使用的ZWD垂直剖面是基于物理的方法构建，气象参数的获取仅需查询以15°纬度间隔提供的气象参数表，因此模型的使用较为便捷。但是这些模型的气象参数存在时空分辨率偏低及在低纬度地区未顾及年际参数变化等不足，从而导致其在某些特定区域存在较大的系统误差，尤其在低纬度地区。TropGrid2格网模型由TropGrid模型改进得到，该模型采用负指数函数来表达ZWD垂直剖面，但是其ZWD的高程缩放因子在全球仅采用了一个常数。与此同时，中国学者在ZWD垂直剖面模型构建方面也开展了广泛的研究，并取得了丰富的成果。常用的区域或者全球对流层延迟模型，其垂直剖面模型主要有二次多项式模型[12]、高度格网模型[13-14]和负指数函数模型[15-19]。但是，对流层垂直剖面的模型参数空间分辨率依然较低。为提高对流层垂直剖面模型参数的空间分辨率，诸多学者基于二次多项式、负指数函数、高斯函数等，利用大气再分析资料月均剖面信息建立了高空间分辨率的ZWD垂直剖面格网模型[20-21]，并取得了较好的ZWD垂直改正效果。

尽管现有的全球ZWD垂直剖面模型均表现出了各自的优越性，但仍存在模型构建仅使用单一格网点数据、建模数据使用月均剖面资料等不足。因此，亟需开展实时高精度全球ZWD垂直剖面格网模型的构建研究。基于以上分析，本文拟引入滑动窗口算法，将全球剖分为大小一致的规则窗口，并结合多年全球的MERRA-2(The second Modern-Era Retrospective analysis for Research and Applications)大气再分析资料积分计算的6 h分辨率ZWD格网点分层剖面信息，在分析ZWD高程缩放因子时空特性基础上，构建顾及ZWD高程缩放因子精细季节变化的ZWD垂直剖面模型，最终建立了高精度全球ZWD垂直剖面格网模型。

1 ZWD高程缩放因子特性分析与GZWD-H模型构建

MERRA-2是由美国航空航天局(National Aeronautics and Space Administration,NASA)提供的最新大气再分析资料(https:∥goldsmr4.gesdisc.eosdis.nasa.gov/data/MERRA2)，其平面分辨率为0.625°×0.5°(经度×纬度)，时间分辨率为6 h。大气再分析资料可以分别提供不同等压层上的分层气象参数，根据这些分层气象数据来计算每层的大气折射率等信息，最终利用积分法即可计算获得每个格网点各层的对流层延迟，但本文未考虑大气湍流和大气波导等情况对格网点各层积分计算ZWD的影响，具体积分计算公式为[22]

e=Sh×P/0.622

(1)

(2)

(3)

由于大气再分析资料在使用前需对其进行精度评估，目前尚无文献对MERRA-2资料在全球计算ZWD的精度予以评估。为此，本文首先利用2015年全球412个探空站的数据来检验MERRA-2大气再分析资料计算ZWD的精度，并统计全球各探空站在0:00和12:00 UTC的ZWD数据与MERRA-2资料积分计算的ZWD的偏差值(bias)和均方根(RMS)误差，统计得到的全球年均偏差值和RMS结果见表1。

表1 2015年全球探空站数据检验MERRA-2资料计算ZWD的精度

由表1可知，MERRA-2资料在全球计算ZWD的平均偏差为4.3 mm，平均RMS为13.5 mm。因此，MERRA-2资料在全球计算ZWD具有较高的精度，可作为全球ZWD垂直剖面格网模型构建的数据源。

1.1 ZWD高程缩放因子时空特性分析

在ZWD的垂直剖面表达方面，诸多学者采用了负指数函数来表达ZWD在高程方向上的变化，且采用该方法对ZWD的垂直改正也得了较好的ZWD垂直改正效果[23-25]。为了进一步验证ZWD在高程上的变化关系，在全球选取了2个具有代表性的MERRA-2再分析资料2017年1月1日0时刻(UTC)的格网点数据，分别积分计算出2个格网点在不同等压层上的ZWD信息，并利用负指数函数来拟合ZWD的垂直剖面，结果如图1所示。

由图1可以看出，利用负指数函数可较好地表达ZWD在高程方向上的变化，从选取的这2个具有代表性格网点来看，ZWD的拟合平均RMS在1 mm以内，表现出较高的拟合精度。因此，本文也采用负指数函数来表示ZWD的垂直剖面，其表达式为

图1 MERRA-2资料计算的2个格网点ZWD值在高程上的变化关系与负指数函数拟合结果及其对应拟合误差Fig.1 The MERRA-2 ZWD profile and cures fitting with the exponential function at two grid points and corresponding fit residual errors of the exponential function

(4)

式中，ZWDr表示其在参考高程Hr处的ZWD值；ZWDt表示其在目标高程Ht处的ZWD值;Hw表示ZWD的高程缩放因子，单位为km。

为了构建更为精细的全球ZWD垂直剖面模型，ZWD高程缩放因子的全球时空特性分析尤为关键。本文首先选取全球具有代表性的6个MERRA-2格网点资料，根据积分法计算出2012—2017年ZWD分层剖面信息，然后根据式(4)计算出ZWD高程缩放因子，最终获取其2012—2017年的日均ZWD高程缩放因子时间序列，并采用年周期和半年周期的余弦函数对Hw进行拟合，结果如图2所示。由图2可知，Hw在全球表现出显著的年周期和半年周期变化。然而，东半球赤道上格网点的值没有其在西半球稳定，可能是受东半球赤道上热带雨林气候和热带季风气候等复杂气候的综合影响。为了分析Hw的年均值、年周期和半年周期振幅的全球分布特性，计算了2012—2017年全球MERRA-2资料每个格网点Hw的年均值、年周期和半年周期振幅，结果如图3所示。由图3可知，Hw在低纬度地区的赤道两端，如南亚、非洲北部、南美洲北部、太平洋南部和印度洋东北部，出现了相对较大的年均值，其在低纬度的太平洋西部、南亚、印度洋北部和非洲南部等地区存在较大的年周期振幅，在低纬度的太平洋中部、非洲北部和西南部、大西洋东部等地区出现了较为显著的半年周期振幅，其原因可能是低纬度的这些区域的水汽变化较其他区域大。因此，为了保证ZWD垂直剖面模型的精度，在模型构建时需同时顾及Hw的年周期和半年周期变化。

图2 MERRA-2资料计算全球6个代表性格网点ZWD高程缩放因子的时间序列变化Fig.2 The time series and corresponding seasonal fits for the ZWD height scale factor at six grid points

图3 全球MERRA-2资料计算ZWD高程缩放因子的年均值、年周期振幅和半年周期振幅分布Fig.3Distributions of the ZWD height scale factor:annual mean, annual amplitude, and semiannual amplitude over globe

1.2 GZWD-H模型建立

针对已有ZWD垂直剖面模型建模时仅采用单一格网点数据以及使用月均剖面数据等不足，本文引入滑动窗口算法将全球剖分为大小一致的规则窗口，以解决存在的问题。在此，先对滑动窗口的算法做如下简述。

对全球格网剖分进行划分，得到与MERRA-2格网资料大小相同的水平分辨率格网，即0.625°×0.5°(经度×纬度)。滑动窗口算法的关键需要确定其窗口的大小及步长，同时窗口大小的确定需顾及全球窗口剖分个数的整数性、窗口的连续性及窗口内模型参数的可求解性等原则。基于此，本文以3行3列(1.25°×1°的区域范围)为一个滑动窗口大小来举例说明滑动窗口算法，其流程如图4所示。具体过程为：①利用格网左上角第一个窗口N1内(每个框表示一个滑动窗口大小)的数据求出其窗口内的相关模型参数，并将其作为窗口N1中心格网点(框内黑点)的结果；②将窗口向纬度东向移动2个格网点，求解新窗口N2内的相关模型参数，将其作为窗口N2中心格网点的结果，以此类推，求出这一纬度上所有窗口内的相关模型参数；③窗口移动到下一纬度(向下移动两个格网点)，以同样的方法求出该纬度所有窗口内的相关模型参数，以此类推，直到求出全球所有窗口内的相关模型参数；④通过以上过程获得全球所有窗口内的相关模型参数，并将其作为各自窗口中心格网点的结果，最后将全球所有窗口中心格网点组建新的全球格网，如图4中的黑点和虚线所示。

图4 滑动窗口算法Fig.4 The sliding window algorithm

为了进一步优化ZWD垂直剖面模型参数及并提高建模数据的利用率，结合上述滑动窗口算法，拟确定窗口大小设为1.25°×1°(经度×纬度)。本文将利用全球各窗口内所有格网点数据来建立对应窗口顾及ZWD高程缩放因子精细季节变化的ZWD垂直剖面模型，其表达式为

(5)

式中，i表示窗口的编号。由于Hw在全球表现出明显的年周期和半年周期，因此每个窗口的Hw可表示为

(6)

针对全球每个窗口，利用窗口内9个MERRA-2格网点2012—2016年6 h分辨率的ZWD分层剖面信息通过最小二乘法则能估计出全球每个窗口ZWD垂直剖面模型的系数。全球每个窗口ZWD高程缩放因子的5个系数以平面分辨率为1.25°×1°(经度×纬度)的格网形式存储，最终构建了高精度的GZWD-H模型。GZWD-H模型的使用非常便捷，其使用过程如下：①用户仅需提供年积日和目标点位置信息，根据目标点位置信息查找与目标点最近的模型参数格网点；②根据查询获得的最近的模型参数格网点的模型参数，并分别利用式(5)和式(6)将目标点在参考高程处的ZWD值垂直改正到目标高程处。

2 GZWD-H模型精度验证

探空站提供的大气分层剖面资料均为实测值，是当前精度最为可靠的大气分层资料之一(http:∥weather.uwyo.edu/upperair/sounding.html)，为验证GZWD-H模型在全球的垂直插值精度和适用性，以2017年全球321个探空站数据积分计算的00:00和12:00 UTC的ZWD分层剖面信息为参考值来验证GZWD-H模型的垂直插值精度，并与GPT2w模型进行精度对比。由于GPT2w模型的模型参数以两种平面分辨率(1°×1°和5°×5°)的格网进行存储，便于后续描述，分别将其定义为GPT2w-1和GPT2w-5。

在本次的模型垂直插值检验中，从探空剖面的地表层开始，依次对探空剖面的相邻两层ZWD信息进行垂直插值(即以其中一层为参考层，则另一层为目标层)，直至插值到探空剖面的顶层。最后对全球每个探空站ZWD的垂直插值的偏差值和RMS进行统计，结果见表2和图5所示。

表2 不同模型对全球2017年各探空站在相邻两层的ZWD垂直插值精度统计

由表2可以看出，所有模型进行全球ZWD垂直插值时均表现为负偏差。在RMS方面，GZWD-H模型在全球的ZWD垂直插值中表现为最小值，相对于GPT2w-1和GPT2w-5模型，其RMS值分别减少了4%和7%。由此表明，GZWD-H模型在全球探空站的ZWD垂直插值中表现出了最优的精度和稳定性。由图5可知，所有模型在全球大部分区域体现为负偏差,而GPT2w-1和GPT2w-5模型在中国东南部和太平洋西部的部分探空站表现为正偏差值。相对于GPT2w模型，GZWD-H模型在全球表现更为稳定，其在南亚和太平洋西部的部分探空站仍具有一定的精度改善。在RMS方面，所有模型在全球中、高纬度地区表现为相对较小的值，而在低纬度地区表现为相对较大的值，尤其在南亚地区，主要原因是在中、高纬度地区的ZWD值远小于低纬度地区，且其变化比低纬度地区更为稳定，在低纬度地区(尤其在南亚)的水汽含量大且变化较为剧烈，因此所有模型在低纬度地区的ZWD垂直插值效果不如其在中、高纬度地区。然而，GZWD-H模型顾及了ZWD高程缩放因子的精细季节变化，其相对于GPT2w模型在南亚和太平洋地区具有一定的精度改善。

图5 不同模型对2017年各探空站剖面在相邻两层的ZWD垂直插值全球精度分布Fig.5 Distributions of bias and RMS of the ZWD layered vertical interpolation at two adjacent levels using different models compared with radiosonde data in 2017 over globe

为了分析ZWD垂直插值偏差值和RMS的季节变化情况，选取了全球3个探空站，统计了其2017年日均偏差值和RMS的时间序列，结果如图6所示。

由图6可知，对北半球高纬度地区的22113探空站，所有模型在夏季存在相对较显著的负偏差值和相对较大的RMS；对于南半球高纬度的89611探空站，所有模型在全年的偏差值和RMS均较小，其偏差值的变化在夏季相对较大，而RMS无明显季节变化；对位于低纬度地区的91413探空站，GPT2w-1和GPT2w-5模型在全年的大部分时间均表现出相对较大且较为剧烈的日均偏差值和RMS，并未发现其明显的季节特性，而GZWD-H模型在全年的时间均具有较小和稳定的日均偏差值和RMS。由此说明在ZWD垂直插值的偏差值和RMS季节变化方面，相对于GPT2w模型，GZWD-H模型在低纬度地区表现出了良好的季节性能。

图6 不同模型在3个探空站ZWD剖面垂直插值的日均偏差和RMS时序变化Fig.6 Variation of average daily bias and RMS timing sequence of the ZWD layered vertical interpolation of the radiosonde with different models

由于ZWD的垂直分布与纬度密切相关，为了分析ZWD垂直插值的偏差值和RMS在纬度上的分布情况，对全球所有探空站按照15°纬度间隔对偏差值和RMS进行分类统计，由于南半球纬度大于60°以上探空站数量较少，为此将该范围的探空站化为一个纬度区间，结果如图7所示。由图7可知，所有模型在大部分纬度区间范围内存在相对较小的负偏差值，在0°S—15°S范围内表现为相对较大的偏差值。所有模型的RMS变化均表现为以赤道向两极逐渐减小的特性，这个特性跟ZWD在全球纬度上的分布特性较为一致。GZWD-H模型在15°N—15°S范围内相对GPT2w模型仍具有一定的精度改善。

图7 不同模型对ZWD垂直插值的偏差和RMS在纬度上的分布Fig.7 Distributions of bias and RMS of the ZWD layered vertical interpolation in different latitude bands using different models

3 利用探空资料检验GZWD-H模型在ZWD空间插值中的精度

以2017年全球321个探空站资料积分计算的12 h分辨率地表ZWD信息为参考值，检验GZWD-H模型在GGOS大气格网ZWD中的空间插值精度。本文使用的GGOS大气格网ZWD产品的水平分辨率为2.5°×2°(经度×纬度)，GGOS大气格网产品可免费下载获取(https:∥vmf.geo.tuwien.ac.at)。由于GGOS大气格网产品的高程是大地高，而探空站的高程是海拔高，两者之间存在高程基准差异，因此在进行ZWD空间插值前，需对两者的高程基准进行统一，本文采用EGM2008模型来实现两者高程基准统一[26]。各模型用于GGOS大气格网ZWD空间插值的精度统计如表3和图8所示。

表3 利用探空资料验证GZWD-H、GPT2w-1和GPT2w-5模型进行GGOS大气格网ZWD空间插值的精度对比

由表3可以看出，针对GGOS大气格网ZWD的空间插值，GPT2w-5具有最大的平均RMS为27.7 mm，相对于GPT2w-5，GPT2w-1的空间插值精度提升了6.1 mm；GZWD-H模型相对于GPT2w-1和GPT2w-5，其RMS分别减少了3.6 mm(17%)和9.7 mm(35%)，这些误差的减少对应在GNSS数据处理中可对高程分量估计的精度分别提升7.2 mm和19.4 mm[24]。图8表明，所有模型在全球高纬度地区均具有良好的ZWD垂直改正精度，GPT2w-5在全球南北半球的中纬度地区出现了相对较大的RMS；此外，GPT2w-1和GPT2w-5在低纬度地区的东南亚、太平洋西部、南美洲中部和北美洲南部以及欧洲中部的部分探空站均出现了较大的偏差值和RMS，尤其在东南亚地区，部分探空站出现了显著的偏差值和RMS，可能是受该地区热带雨林气候和热带季风气候的影响，导致该地区温度、水汽压等气象因子相对于中高纬度地区较为活跃，难以对气象因子精确模型化，进而导致GPT2w-1和GPT2w-5的精度较低。而GZWD-H模型充分利用了多个格网点的ZWD垂直剖面信息进行建模且顾及了ZWD高程缩放因子的精细季节变化，因此在全球范围内均保持了较优和稳定的ZWD垂直改正性能，尤其在低纬度地区，GZWD-H模型表现出显著的优势。

图8 利用探空资料验证GZWD-H、GPT2w-1和GPT2w-5模型进行GGOS大气格网ZWD空间插值的偏差和RMS全球分布Fig.8 Distribution of bias and RMS of the spatial interpolation for GGOS Atmosphere gridded ZWD compared with radiosonde data

4 结 论

ZWD垂直剖面模型是进行ZWD垂直改正的关键，也是高精度对流层延迟模型构建的基础。本文基于滑动窗口算法建立了全球ZWD垂直剖面格网模型(GZWD-H模型)。同时，以2017年全球321个探空站资料作为参考值，检验了GZWD-H模型在全球的垂直插值精度及其在ZWD空间插值中的应用，并与全球精度优异的GPT2w模型进行了精度对比。结果表明：以全球探空站资料计算的ZWD分层剖面信息为参考值，在ZWD的垂直插值检验中，相对于GPT2w-1和GPT2w-5模型，GZWD-H模型在全球范围内均表现出了最优的精度和稳定性，尤其在南亚和太平洋地区，其体现出了显著的精度改善。以全球探空数据计算的地表ZWD信息为参考值，GZWD-H模型在全球GGOS大气格网ZWD空间插值中均呈现出了最优的精度和稳定性，尤其在低纬度地区，其相对于GPT2w模型具有显著的优势。相对于GPT2w-1模型计算ZWD垂直改正时所需的模型参数，GZWD-H模型的模型参数减少了3.75倍。GZWD-H模型在全球进行ZWD垂直改正时，可在一定程度上提升模型的计算效率。由于GZWD-H模型在全球大气格网对流层产品的空间插值中表现出优异的性能，从而促进其在高精度GNSS大气探测和GNSS精密定位中的应用。本文在构建GZWD-H模型时只考虑了固定的窗口大小和滑动步长，后续研究将采用高分辨率的ERA5资料分析不同窗口大小和步长的选取对全球ZWD垂直剖面建模的影响及模型在实时PPP中的改正精度。

致谢：感谢美国NASA提供的MERRA-2格网资料、美国怀俄明大学提供的探空数据和GGOS Atmosphere提供的ZWD格网数据。