水汽是地球大气的重要组成部分，在气候与天气现象中扮演着关键角色。传统的水汽探测手段存在空间或时间分辨率不足、精度有限或者适用范围较小等局限性。相比传统的大气水汽探测技术，地基GNSS探测水汽技术具有实时性、全天候、高精度等优势，且单站布设成本较低，一旦布设完毕即可投入连续运行，可实现大范围高密度的水汽监测。

对流层延迟是GNSS定位的重要误差源。天顶对流层延迟包含天顶静力延迟(zenith hydrostatic delay, ZHD)与天顶湿延迟(zenith wet delay, ZWD)两部分。1992年，Bevis等^[1]提出了“GPS meteorology”的理念，即利用GPS ZWD反演GPS天顶可降水量(precipitable water vapor, PWV)，建立了GNSS对流层延迟与大气水汽之间的关系。Tong等^[2]通过各种水汽监测手段证实地面GPS测量是有用的水汽长时间变化的监测技术；宋淑丽等^[3]对利用地基GPS探测大气中水汽的原理和技术做了深入系统分析，并将其应用于气象学中；陈永奇等^[4]探讨了香港实时GPS水汽监测系统的若干关键技术，为GNSS水汽业务化提供了技术支持；于胜杰等^[5]对不同约束条件对GPS水汽层析解算的影响进行了分析，为不同条件下的水汽解算提供相应的约束条件；张双成^[6]通过对影响地基GPS探测水汽精度的因素做深入分析，搭建了实时地基GPS遥感水汽系统; 金双根等^[7]采用GPS解算的PWV和其他数据第一次在全球范围内研究多尺度的水汽变化，发现在局部地区夏季的水汽最大，冬季的水汽最小的季节性规律。

然而，应用GNSS高分辨率数据分析全球及区域性的对流层水汽时空特征等相关研究仍存在不少问题。为了挖掘GNSS高分辨率的水汽时间序列中蕴藏的有用信息，需要采用相关方法对时间序列进行分析。小波分解(wavelet decomposition, WD)算法是一种常用的从时间序列中获取周期信息的方法，本文选用了一种自适应的盲源分解方法-经验模态分解(empirical mode decomposition, EMD)方法。该方法是一种自适应性信号处理滤波方法，可以将原始时间序列分解成有限个频率逐渐降低的本征模态分量(intrinsic mode function, IMF)，对于分析时间序列信号具有重要意义。

为了研究GPS水汽时间序列的时空特征，本文利用EMD和WD方法，在分析两种分解方法的利弊之后，选取合适的方法进行水汽原始时间序列周期性分析，再结合地理环境、地形、纬度等因素探讨水汽变化的规律以及引起变化的物理原因。

1 水汽解算原理

GPS卫星发射的电磁波信号在大气中传播时会受到对流层的影响，导致传播速度减慢，传播路径发生弯曲。因此，假设对流层某处的大气折射系数为n，则电磁波信号在该处的传播速度V=c/n，c为光速。所以当电磁波信号在对流层中传播时，由对流层折射引起的过剩路径长度—对流层延迟ΔL可表示为^[8]：

式中，R为未受对流层折射的卫星到接收机间的直线距离^[9]。地基GPS气象学中通常用天顶总延迟T与天顶静力延迟H的差来间接计算天顶湿延迟W^[10]，即：

利用天顶静力学延迟模型，结合测站气象参数可以计算静力学延迟，并进一步得到ZWD^[11]。地基GPS气象学的主要思想是通过获取高精度的测站天顶方向GPS信号的湿延迟来确定可降水量，由ZWD可以通过式(3)计算测站上空的可降水量P^[12]：

式中，∏是无量纲的水汽转换系数。

2 时间序列分解算法

2.1   EMD算法

EMD方法对原始时间序列进行分解，可以得到有明确物理意义的本征模态分量^[13]。EMD可以将时间序列信号分解为一系列具有不同时间尺度的IMF和一个余项，每个IMF需要符合两个条件：①在全体数据序列中，过零点与极值点的个数必须相同或者最多相差一个；②在任意一点，由局部极小值点组成的下包络线和局部极大值点确定的上包络线的算术平均值为零。

分解过程通过一个称为“筛选”的步骤来完成^[14]，原始时间序列D(t)可表示为：

式中，r_n代表最终余项，表示信号的趋势成分；C₁, C₂, ..., C_n表示频率从高到低排列的n个IMF分量。

2.2   小波分解算法

小波变换(wavelet transform，WT)是一种新的变换分析方法，它继承和发展了短时傅里叶变换局部化的思想，同时又克服了窗口大小不随频率变化等缺点，能够提供一个随频率改变的“时间-频率”窗口，是进行信号时频分析和处理的理想工具。根据小波分析理论，如果二进离散小波函数簇ψ_{j, k}(t), j, k∈Z构成L²(R)中的标准正交基，对于时变信号f(t)∈L²(R)的正交小波分解有：

式中，N为分解层数；d_{j, k}为小波分解系数；c_{j, k}为尺度分解系数；φ_{j, k}为分解后的残余序列。

3 水汽解算算例分析

3.1   数据来源

本文选取了大气可降水服务网(http://120.126.139.90:8080/)提供的水汽解算数据，该网站提供中国台湾及福建离岛共计100个测站2006-2013年时间间隔为1 h的GPS-PWV资料，使用Bernese 5.0软件来处理GPS观测资料^[15]。

3.2   算例分析

文章利用中国台湾高分辨率GPS-PWV资料，由每个GPS站的2006-2013年PWV时间序列计算得到PWV的日平均值，图 1展示了PWV的分布特性，并据此研究水汽变化的影响因素。

图  1  台湾GPS站平均PWV与地理位置、地形、纬度的对应情况

Figure  1.  The Average PWV with the Interaction of Geographical Location, Topography, Latitude

下载: 全尺寸图片 幻灯片

图 1(a)表示平均PWV随地理位置变化的分布情况，图 1(b)表示水汽含量随着高程的变化情况，图 1(c)展示了在GPS站上平均PWV与纬度的关系。因为水汽并不是单一因素控制，而是受到多种因素的影响，这就导致图 1(b)、1(c)中出现偏离红线较远的点。例如，在图 1(b)、1(c)中，蓝色点为同一站点，而该点在图 1(c)中偏离红线较远，这是因为该站点处在高海拔低纬度地区，在高海拔低水汽的影响下，并没有表现出如其他点的低纬度高水汽现象。除去与该点类似的交互影响情况，图 1(b)、1(c)中表现出明显的相关趋势。

从图 1(a)中可以明显看出，GPS站的平均水汽范围为13.0~43.0 mm，平均水汽的分布特点如下。

1) 沿海岸线的各站平均水汽(40 mm)明显大于远离海岸的内陆站(30 mm)，这说明季风对台湾地区可降水量的影响：东亚季风气候区夏季风强，经太平洋补充，带来了充足的水汽；

2) 中国台湾东部沿海站的平均水汽大于西部沿海站。这一现象是中国台湾独特的气候和地形特征所导致的。中国台湾地区地处东南季风地区，为海洋性气候，夏季风从东面而来，中国台湾有南北延伸高耸的山脉，东部湿润的海风被阻于山的东部，由于气流上升而气温下降造成降雨，所以中国台湾东部比西部水汽大。

3) 中国台湾低纬度地区的平均水汽大于高纬度地区的平均水汽。图 1(c)显示，随纬度的递增，平均PWV从43 mm递减为33 mm。而低纬度地区温度高，这表明大气中水汽的分布和气温有关。

4) 水汽含量随着高程的增高降低。在图 1(a)中，平均水汽的最低值都出现在高山地区，这是因为季风带来的水汽受到高山阻挡。而图 1(b)则显示出高程和水汽的负相关性，在高程从0 m递增为3 000 m的时候，平均PWV也从40 mm递减为15 mm。究其原因，高山地区的水汽主要来自地表水蒸发，因此，随高度上升，水汽含量减少。

总而言之，中国台湾地区的大气水汽分布主要受纬度、地形特征和气候条件的控制^[15]。

3.3   EMD分解结果

本文通过EMD算法分解了中国台湾各GPS站在长时间序列(2006-2013)和短时间序列(2013-12-24-2013-12-31)的水汽时间序列，得出了频率依次递减的IMF。这里以KIYN站和SAND站为例，结果如图 2所示，图中横轴为年积日(DOY)或年份。

图  2  KYIN和SAND站长时间和短时间序列以及相应的EMD-IMF

Figure  2.  The Long and Short Time Series and the Corresponding EMD-IMF of KYIN and SAND

下载: 全尺寸图片 幻灯片

在图 2(b)中，KYIN站整天周期振幅为2 mm，且波形较均匀、规整，周期表现明显；半天周期振幅为0.9 mm，但该分量的周期表现不统一，在2013-12-24-2013-12-27以及2013-12-31，半天周期明显，即频率信号均匀统一，而在这一周内的其他天，周期并不为明显的12 h。上述分析说明，EMD并不能很好地分解出0.5 a和0.5 d的周期分量信号。本文将用WD方法分解水汽时间序列，以此探究WD方法是否能分解出明显的半年以及半天的周期分量。

3.4   WD分解结果

在有关信号分解获取周期信号的分析中，常用的算法是小波分解算法。为了对比分析EMD方法与小波分解方法的优劣，也为了进一步验证EMD分解结果的可信度，本文采用Mallat快速算法及Danbechies4小波算法对台湾KYIN站和SAND站的原始信号进行8层分解，结果如图 3所示。

图  3  KYIN站和SAND站长时间序列和短时间序列小波分解结果

Figure  3.  The Long and Short Time Series and the Corresponding the WD-IMF of KYIN and SAND

下载: 全尺寸图片 幻灯片

图 3(a)、3(c)显示了小波分解长时间原始序列及其分解结果，蓝色线为原始序列，红色线为趋势项，橙色线为周年振荡的IMF，绿色线为半周年振荡的IMF，中国台湾其他各GPS站都表现出了周年和半周年振荡的IMF。图 3(b)、图 3(d)显示了小波分解短时间序列及其分解结果，蓝色线为原始序列，红色线为趋势项，紫色线为一天周期振荡的IMF，粉色线为半天周期振荡的IMF，台湾其他各GPS站都表现出了一天和半天(12 h)振荡的IMF。

在图 3(a)中，KYIN站年周期振幅表现与图 2(a)中的表现一致；但半年周期振幅表现比EMD的分解结果好，周期明显、振幅均匀，振幅为3 mm。在图 3(b)中，KYIN站一天周期振幅表现与图 2(a)中的表现一致；但半天周期振幅表现比EMD分解结果好，周期明显均匀，可以看出是规律的半天周期分量。

由图 3、图 4可以看出，EMD算法和小波分解算法都可以分解得到周年、半周年、一天、半天的周期分量，甚至在振幅、波峰的表现都如出一辙。

图  4  周年、半周年、一天、半天周期振荡IMF振幅分布情况(白色点为异常高值)

Figure  4.  Amplitude of Annual, Semiannual, Diurnal, Semidiurnal Oscillation

下载: 全尺寸图片 幻灯片

1) 图 2(a)、图 2(c)和图 3(a)、图 3(c)中，橙色线对应的年周期振荡分量的最大峰值集中在夏季，即7~8月；相反，水汽的最小值出现在冬季。水汽最大值在夏季的明显季节规律是由于受到夏季湿季风的影响，最小值出现在冬季是受低温的影响。图 2(a)、图 2(c)和图 3(a)、图 3(c)中，绿色线对应的半周年振荡分量的最大波峰在2~4月和8~10月。

2) 图 2(b)、图 2(d)和图 3(b)、图 3(d)中，紫色线对应的1 d周期振荡分量变化的波峰和波谷分别为正午(12时)和午夜(24时)。图 3(b)、图 3(d)中，粉色线对应的半天振荡分量的波峰在早晨6点(6时)和晚上6点(18时)。

但是两种算法在周期信号分解的应用中各有利弊。经小波分解后，获得的半周年分量和半天分量(图 3)的周期性明显比EMD算法获得的半周年分量(图 2)的周期性更加明显且稳定。这说明小波分解在长时间大量数据的信息处理中，对细节信息更加敏感。但小波分解受到分解层数和小波基的限制，在无经验模态分解提供的先验信息(将长时间序列分解至8层左右可以获得半周年和周年周期信号，短时间序列分解至4层左右可以获得一天和半天周期信号)的情况下，不易获得准确周期信号, 且利用小波分析并不能获得原始信号的趋势项。

本文认为先利用EMD方法对水汽时间序列进行盲源分析，即在没有任何先验信息的情况下进行分解，在得到初步周期信息和分解层数等信息后，再用WD方法进行分解，以此得到明显、均匀的周期分解分量。本文对中国台湾地区每个站的水汽时间序列都做了分解，并统计周期、振幅信息，来研究这种周期现象产生的原因。

3.5   物理原因分析

为了研究PWV水汽周年、半周年、一天、半天的振荡曲线的振幅和地形、地理环境的对应关系，将PWV原始序列按照§3.4所述的方法进行分解，得到的周年、半周年、一天、半天周期对应的IMF振幅的分布情况显示在图 4(a)~图 4(d)中。

图 4水汽分布图中，所表现出的周年、半周年、一天、半天周期振荡的特征及其出现的物理原因分析如下。

1) 水汽的季节性变化

图 4(a)展示中国台湾GPS站观测值年周期振荡的PWV的振幅(单位为mm)。较高振幅在西沿海地区，约为27±5 mm，较低振幅分布在内部高山地区，约为13±5 mm。这是由于台湾冬季受大陆冷气团影响，东北季风盛行，雨影效应引起西部水汽值小于东部；而当夏季时，受太平洋副热带高气压影响，西南季风盛行，来自太平洋的东南季风受到中央山脉阻挡，西部水汽值增多。即使得西部水汽含量在冬季达到最小值，夏季达到最大，故认为季风在台湾特殊地形作用下影响了年周期变化的PWV。

图 4(b)中，半周年周期振荡的振幅比周年振荡的振幅弱，在2~12 mm左右，振幅高值出现在沿海区域，振幅低值出现在内部高山区域，但是并没有出现显著的与地理环境间的对应关系。台湾冬季受大陆冷气团影响，东北季风盛行，大陆南下的冷气团经东海带来丰富的水汽；夏季受太平洋副热带高气压影响，西南季风盛行，各地气温普遍升高，空气中水汽增多，对流作用强烈，带来大量降水。因此可认为半周年的振荡主要是由于夏季风与冬季风对台湾岛控制作用的转换引起的。

2) 水汽的昼夜变化

如图 4(c)所示，台湾所有的GPS站的日PWV变化的范围为0.3~5 mm，而且振幅较大的站集中在沿海和内陆高山站，日PWV变化的高峰期为正午和午夜。这主要是由于中国台湾独特的地理位置和地形特征所引起的以1 d为周期转换的两种风系——海陆风和山谷风所导致的。沿海站的高振幅值主要是由于海陆风环流引起的。白天，陆地增温比海洋快，陆地上的气温比海上高，因而形成局地环流，下层风由海洋吹向陆地，使沿海站水汽增加；夜间，陆地降温快，地面冷却，海面降温慢，海面气温高于陆地，产生与白天相反的热力环流，下层风自陆地吹向海洋，使沿海站水汽值减小。而高山站的高振幅主要是由于海陆风-山谷风叠加环流引起的。在山区，白天日出后，山坡受热，其上空气增温快，而同一高度的山谷上空的空气距地面较远，增温较慢。此时由于海陆风，山坡聚集了大量从海上吹来的水汽，暖湿空气沿山坡上升，风由山谷吹向山坡，导致高山站水汽值增加。夜间辐射冷却，气温迅速降低，而同一高度山谷上空的空气冷却较慢，山坡上的空气下滑，下层风由山坡吹向山谷，又因为海陆风环流，进一步吹向大海，形成与白天相反的热力环流。

半天(12 h)PWV变化(图 5(d))比日PWV变化小很多，振幅小于4 mm。半天的PWV振荡与地理环境并没有显著对应关系。通常，半天水汽振荡的峰值发生在早晨6点和下午6点，这分别是一天太阳升落的时刻，也是温度急剧变化的时刻，地面受热不均导致空气作垂直运动(受热上升，冷却下降)，即由于太阳辐射加热引起的局地热对流导致了PWV的半天周期变化。此外，低层大气水分聚集和降雨、地面水汽输送以及其他因素也会对PWV天振荡产生影响。

4 结语

本文以台湾GPS水汽时间序列为研究对象，分析对流层水汽与地理环境的对应关系，并利用经验模态分解方法和小波分解方法将水汽时间序列进行分解，以此来研究水汽时间序列的周期性振荡特性，得到以下结论。

1) 台湾地区的大气水汽分布主要受纬度、地形特性和气候条件控制。

2) EMD算法和小波分解算法都可以分解得到周年、半周年、一天、半天的周期分量，甚至在振幅、波峰的表现都如出一辙。但是小波分解对细节信息更加敏感，且受到分解层数和小波基的限制。据此提出先利用EMD方法对水汽时间序列进行盲源分析，得到初步周期信息和分解层数等信息后，再用小波分解方法进行分解。

3) 在PWV年周期变化的振幅表现中，较高振幅分布在西部沿海地区，较低振幅在内陆高山地区。季风在台湾特定的地形条件下作用导致了PWV的年周期变化。半周年的振幅相比周年振幅弱一些，且并没有出现显著的与地理环境的对应关系。半周年的振荡主要是由于夏季风与冬季风对台湾地区的交替控制导致的。水汽的天变化中，振幅较大的站集中在沿海和内陆高山站，分别是由于海陆风环流、海陆风-山谷风叠加环流引起的。半天的PWV变化振幅较小。半天的PWV振荡与地理环境并没有显著对应关系。半天PWV周期主要是因为受到了太阳辐射加热引起的局地热对流的影响。