基于高斯函数的ZTD模型系数时空特性分析

张继洪,谢劭峰,张亚博,曾 印,唐友兵,熊 思

(1.桂林理工大学 测绘地理信息学院，广西 桂林 541004；2.湖北科技学院 资源环境科学与工程学院，湖北 咸宁 437100)

0 引言

对流层延迟是电磁波信号在传播过程中的主要误差源，其在天顶方向上的延迟约为2 m，随着卫星高度角的降低，其影响将增至20 m[1-2]。对流层天顶延迟(Zenith Total Delay，ZTD)模型对全球导航卫星系统(Global Navigation Satellite System，GNSS)实时高精度导航定位及GNSS水汽反演具有重要意义。在已有的ZTD模型中，Hopfield模型[3]和Saastamoinen模型[4]虽然可以提供精度较高的ZTD，但其应用时需要输入实测气象参数；UNB3m模型[5]和EGNOS模型[6]虽然无需输入实测气象参数，但其较低的时空分辨率影响了ZTD计算精度。近年来，众多学者利用大气再分析资料计算的ZTD进行ZTD模型的建立，但由于ZTD在垂直方向的变化比水平方向的变化剧烈且复杂，大气再分析资料格网点与用户位置存在高度差，由格网点ZTD插值到用户位置时将产生较大误差。为此，有学者在对流层垂直剖面函数方面开展了广泛的研究[7-8]，通常采用二次多项式[9]、负指数函数[10-15]和高斯型函数[16]进行模型构建。研究发现，利用高斯型函数拟合垂直剖面ZTD效果略优于二次多项式和负指数型函数。以上模型还存在时空分辨率较低、在部分区域(如广西地区)适应性较差等不足。因此亟需构建适用于广西地区高精度、高时空分辨率的ZTD垂直剖面模型。

本文利用广西地区2016—2017年高时空分辨率的ERA5再分析资料积分计算的分层ZTD，经高斯函数对垂直剖面ZTD进行模型化，分析模型系数b与模型系数c的精细时空特性，为广西地区高精度ZTD垂直剖面模型的构建提供参考。

1 数据来源与计算原理

1.1 数据来源

ERA5大气再分析资料(https:∥apps.ecmwf.int/datasets/data/interim-full-daily)是欧洲中期天气预报中心(European Centre for Medium-Range Weather Forecasts,ECMWF)第5代全球气候再分析数据集，该数据集提供逐小时(UTC时)覆盖全球的气压、温度和比湿等气象参数，包括地表数据和分层数据，水平分辨率为0.25°×0.25°(纬度×经度)，垂直方向包括37层，以标准气压层进行分层。ERA5大气再分析资料是目前时空分辨率最高的大气再分析资料，是进行大气科学研究和对流层模型构建的重要数据源。文献[17]利用桂林市GNSS观测站数据解算的ZTD对ERA5大气再分析资料在桂林市计算ZTD的精度进行了评估，发现ERA5大气再分析资料在该地区计算ZTD具有较高的精度。本文选取2016—2017年广西地区(104°～113°E，20°～27°N)的ERA5大气再分析资料。

1.2 计算原理

根据ERA5大气再分析资料提供的地表和分层数据，利用积分法计算广西地区每个格网点的垂直分层ZTD[18]：

(1)

(2)

(3)

式中，e为水汽压，单位hPa；Sh为比湿；P为气压，单位hPa；N为总折射率；T为温度，单位K；H为高程，单位km；hlow为大气积分计算的最底层高度；htop为大气积分计算的最顶层高度；k1，k2，k3均为常系数，其值分别为77.604 K/Pa，64.79 K/Pa，375 463 K2/Pa。

由于大气再分析资料ERA5顶层存在残余大气，且主要存在对流层天顶静力学延迟(Zenith Hydrostatic Delay，ZHD)，而对流层湿延迟(Zenith Wet Delay，ZWD)较小，可忽略不计。为了提高ZTD的积分计算精度，利用Saastamoinen模型计算ERA5再分析资料顶层上大气残余ZHD，将其附加在每一层的积分结果上：

(4)

式中，Ptop为对流层顶层大气压，单位hPa；φ为纬度；ZHDSaas为Saastamoinen模型计算的ERA5层顶残余ZHD。由于ERA5数据是离散的分层数据而非连续数据，因此将式(3)进行离散化：

(5)

式中，Ni为第i层大气的总折射率；ΔHi为第i层大气的厚度；n为积分层数。

鉴于高斯函数对垂直剖面ZTD的拟合效果略优于二次多项式和指数形式，因此，本文选用高斯函数对垂直剖面ZTD进行拟合：

(6)

式中，a为模型系数，单位m；b，c也为模型系数，单位km；e为自然常数；h为高程，单位km；ZTD(h)为高程h处的ZTD值，单位m。为了验证高斯函数对垂直剖面ZTD的拟合效果，选取广西地区具有代表性的4个格网点，通过积分法计算2017年1月1日00:00(UTC时)的分层ZTD，采用高斯型函数对其进行拟合，结果如图1和表1所示。

(a)26°N,111°E

表1 ZTD在垂直剖面的高斯函数拟合偏差和均方根误差

由图1可知，每个格网点不同高程的ZTD散点与高斯函数拟合线均具有很高的吻合度，4个格网点在10 km以下均表现出较大的偏差，5 km以下表现出正偏差，5～10 km表现出负偏差，4个格网点在20 km以上均表现出较小的偏差，说明高斯函数在较大的高度时拟合垂直剖面ZTD效果更佳。由表1可以看出，4个格网点的平均bias均为负值，说明高斯型函数拟合垂直剖面ZTD具有一定的系统误差，但仅约-0.6 mm。4个格网点的均方根误差分别为11.439，8.864，7.573，7.592 mm，总体来看，采用高斯函数对垂直剖面ZTD进行拟合具有较高的精度和稳定性。

为了使式(6)能够对任意高度的ZTD进行高程改正，将参考高程的ZTD插值到目标高程，因此，对式(6)进行推导[16]：

(7)

式中，e为自然常数；hr和ht分别为参考高程和目标高程，单位km；b和c为模型系数，单位km；ZTDr为参考高程hr处的ZTD值，单位m；ZTDt为目标高程ht处的ZTD值，单位m。利用式(7)即可将参考高程的ZTD插值到目标高程，更有利于实际应用。

2 ZTD模型系数的时空特性分析

为了探究模型系数b和模型系数c的周期性，分析二者的精细时间变化特性，根据广西地区2016—2017年ERA5大气再分析资料积分计算的分层ZTD，采用高斯函数拟合出该研究区域内所有格网点逐小时(UTC时)的模型系数b与模型系数c。选取具有代表性的一个格网点(25.75°N，107.5°E)，计算该格网点2016—2017年的日均模型系数b与模型系数c，对二者进行年周期和半年周期的三角函数拟合，并采用快速傅里叶变换分析周期性，结果如图2所示。

由图2(a)可知，模型系数b的值为-100～-40 km，具有显著的季节变化特征，具体表现为春冬季节较大，夏季较小，且具有明显的周期性；由图2(c)的频谱分析可以进一步看出，模型系数b具有显著的年周期和半年周期。由图2(b)可知，模型系数c的值为25～40 km，模型系数c同样具有显著的季节变化特征，具体表现为春冬季节较小，而夏季较大，同样具有明显的周期性；由图2(d)的频谱分析可以进一步看出，模型系数c也具有显著的年周期和半年周期。

(a)b的日均时间序列及其拟合线

为了分析模型系数b与模型系数c更精细的时间变化特征，同样选取上述格网点，计算年均UTC00:00—23:00的逐小时模型系数b与模型系数c，对其进行日周期和半日周期的三角函数拟合，并采用快速傅里叶变换分析周期性，结果如图3所示。

(a)b的年均逐小时时间序列及其拟合线

由图3(a)可知，模型系数b具有显著的日变化特征，具体表现为从00:00(UTC时)开始，先增大后减小最后又增大的变化趋势，其值在夜晚较大，白天较小,且具有明显的周期性;由图3(c)的频谱分析可以进一步看出，模型系数b具有显著的日周期和半日周期。由图3(b)可知，模型系数c同样具有显著的日变化特征，具体表现为从00:00(UTC时)开始，先减小后增大最后又减小的变化趋势，其值在夜晚较小，白天较大,同样具有明显的周期性;由图3(d)的频谱分析可以进一步看出，模型系数b也具有显著的日周期和半日周期。因此，在采用高斯型函数构建垂直剖面ZTD模型时，需要考虑模型系数b与模型系数c的年周期、半年周期、日周期以及半日周期。

为了分析模型系数b与模型系数c的年均值、年周期振幅、半年周期振幅、日周期振幅及半日周期振幅在广西地区的区域分布特征，对每个格网点的模型系数b与模型系数c进行顾及年周期、半年周期、日周期及半日周期拟合，并计算每个格网点模型系数b与模型系数c各自的年均值、年周期振幅、半年周期振幅、日周期振幅及半日周期振幅。对离散格网点上的数值进行水平插值，结果如图4和图5所示。

由图4可以看出显著的区域性，其中，模型系数b年均值在广西西北地区较大，在广西东南地区较小，且呈现由西北向东南逐渐减小的变化趋势，该变化可能与广西西北地区海拔高、东南地区海拔低有关；年周期振幅由北向南逐渐减小，该变化可能由广西地区的南部临海气候和北部内陆气候差异所致。而半年周期振幅、日周期振幅和半日周期振幅没有明显的区域变化规律，半年周期振幅在广西西南地区较大，在东南地区较小；日周期振幅在西北和西南地区较大，在中部较小；半日周期振幅在西南地区较大，在东北地区较小。由图5可看出，模型系数c的年均值及各振幅的区域分布特征与模型系数b的区域分布特征有极高的相似性。模型系数c的年均值与年周期振幅同样具有显著的区域变化规律，而半年周期振幅、日周期振幅和半日周期振幅也没有明显的区域变化规律。由于模型系数b与模型系数c的年均值与各振幅在不同地理位置存在较大的差异，因此，在采用高斯函数构建垂直剖面ZTD模型时，每个格网点均需要考虑模型系数b与模型系数c的年周期、半年周期、日周期以及半日周期。

(a)年均值

(a)年均值

为了分析模型系数b和模型系数c与经度的相关性，选取23°N与25°N两个纬度上2017年的模型系数b和模型系数c，计算每个格网点的年均b值与年均c值，结果如图6所示。

(a)23°N

为了分析模型系数b和模型系数c与纬度的相关性，选取108°E与111°E两个经度上2017年的模型系数b和模型系数c，计算每个格网点的年均b值与年均c值，结果如图7所示。

(a)108°E

图6(a)和(b)分别为模型系数b在23°N和25°N两个纬度上随经度变化的散点分布及其拟合线，从图中可看出，模型系数b与经度呈现出较为显著的负相关性；图6(c)和(d)分别为模型系数c在23°N和25°N随经度变化的散点分布及其拟合线，从图中可看出，模型系数c与经度呈现出较为显著的正相关性。图7(a)和(b)分别为模型系数b在108°E和111°E两个经度上随纬度变化的散点分布及其拟合线，从图中可看出，模型系数b与纬度呈现出较为显著的正相关性；图7(c)和(d)分别为模型系数c在108°E和111°E两个经度上随纬度变化的散点分布及其拟合线，从图中可看出，模型系数c与纬度呈现出较为显著的负相关性。对比图6和图7可知，模型系数b、模型系数c与纬度的相关性更高。

为了进一步分析模型系数b、模型系数c与经纬度的相关性，利用皮尔逊相关性分析法计算b值、c值与经纬度的相关系数，结果如表2所示。

表2 模型系数b/c与经纬度的相关系数

在皮尔逊相关性分析法中，2组数据的相关系数越接近1或-1，说明2组数据相关性越高，相关系数为正时，说明2组数据呈正相关，相关系数为负时，说明2组数据呈负相关。由表2可知，模型系数b与经度的相关系数为-0.824，二者呈负相关关系，相关性较高；模型系数b与纬度的相关系数为0.969，二者呈正相关关系，相关性较高，相比经度，模型系数b与纬度的相关性更高。模型系数c与经度的相关系数为0.822，二者呈正相关关系，相关性较高；模型系数c与纬度的相关系数为-0.968，二者呈负相关关系，相关性较高，相比经度，同样可以看出模型系数c与纬度的相关性更高。由于模型系数b、模型系数c与经纬度均存在较高的相关性，因此，在采用高斯型函数构建垂直剖面ZTD模型时，均需要考虑模型系数b、模型系数c与经纬度的相关性。

3 结束语

针对广西地区地形起伏较大、气候复杂的特点，已有的低时空分辨率ZTD垂直剖面模型难以满足该地区的需求，根据2016—2017年广西地区ERA5大气再分析资料，分析基于高斯函数的ZTD垂直剖面模型的模型系数b与模型系数c的时空分布特征，结果表明：

① 模型系数b与模型系数c均存在季节变化特征，二者季节变化差异较大，模型系数b的数值表现出春冬季节较大而夏季较小的特征，模型系数c的数值表现出春冬季节较小而夏季较大的特征，二者均存在显著的年周期与半年周期。

② 模型系数b与模型系数c在日变化方面表现出一定的差异性，模型系数b的数值在夜晚较大而在白天较小，模型系数c的数值与模型系数b的情况相反，二者均存在显著的日周期和半日周期。

③ 模型系数b与模型系数c的年均值、年周期振幅、半年周期振幅、日周期振幅及半日周期振幅在广西地区存在显著的区域性特征，其中，b值、c值的年均值与年周期振幅存在沿某一方向逐渐变化的规律。

④ 模型系数b、模型系数c与经纬度均存在较高的相关性，二者与纬度的相关性高于与经度的相关性。

综上可知，在广西地区构建基于高斯函数的ZTD垂直剖面模型时需考虑模型系数b与模型系数c的年周期、半年周期、日周期、半日周期，以及二者与经纬度的相关性。研究结果可为构建基于高斯函数的ZTD垂直剖面模型提供参考。