1、青海湖地区近 50 年降水量周期变化分析 杨沈斌，张弥，吕开龙 1 问题 青海湖位于 青海省东北 部的青海湖 盆地内，既 是中国最大 的内陆湖泊 ，也是中国 最 大 的咸水湖。 近年来有多篇关于青海湖面积受气候变化影响的研究报告， 认为温度升高， 降水 减少和蒸发量大是造成面积下降的几个重要原因。 为此， 本实验拟采用连续小波分析方法对 青海湖地区年降水量周期变化进行分析，探讨降水量变化与青海湖面积变化的关系。 2. 资料 考虑到刚察气象观测站正好处于青海湖边上， 所得数据对反映青海湖周边降水变化具有 较好的代表性。因此，以青海湖西北位置的刚察气象台站 1961-2007 年 年降水量资料为例
2、， 对该站年降水资料进行小波分析， 获取其周期变化特征。 同时， 获取了青海湖 1961-2007 年 的面积变化资料。 从资料发现， 青海湖面积近 50 年呈现下降 趋势， 平均下降 3 km 2 /a 。图 1 和图 2 分别 显示了 1961-2007 年刚察 气象站年降水量和青海湖面积的变化， 以及对应的趋势 线。 图 1 1961-2007 年刚察气象台站年降水量及趋势线 图 2 1961-2007 年青海湖面积及趋势线 3. 小 波分析 方法 小波变换具有多分辨率分析的特点，并且在时频两域都具有表征信号局部特征的能力。小波变换通过将时间系列分解到时间频率域内， 从而得出时间系列的显
3、著的波动模式， 即周 期变化动态， 以及周期变化动态的时间格局 （Torrence and Compo, 1998 ） 。小 波（Wavelet ） ， 即小区域的波， 是一种特殊的、 长度有限， 平均值为零的波形。 它有两个特点： 一是 “小” ， 二是具有正 负交替的“ 波动性” ，即 直流分量为 零。小波分 析是时间（ 空间）频率 的局部化 分析， 它通过伸缩平移运算对信号 （函数） 逐步 进行多尺度细化， 能自动适应时频信号分析 的要求， 可聚焦到信号的任意细节。 小波分析将信号分解成一系列小波函数的叠加， 而这些 小波函 数都 是由一 个母 小波（mother wavelet）函 数
4、经过 平移 与尺度 伸缩 得来的 。用 这种不 规则的小波函数可以逼近那些非稳态信号中尖锐变化的部分， 也可以去逼近离散不连续具有 局部特性的信号， 从而更为真实的反映原信号在某一时间尺度上的变化。 小波分析这种局部 分析的特性使其成为对非稳态、 不连续时间序列进行量化的一个有效工具 （Stoy et al., 2005 ） 。 小波是一个具有零均值且可以在频率域与时间域内进行局部化的数学函数（Grinsted et al., 2004 ） 。 一个小波被称为母小波 （mother wavelet ） ， 母小 波可沿着时间指数经过平移与尺 度伸缩得到一系列子小波。 子小波可以通过尺度 （s
5、， 频率的反函数） 函数和时间 （n）位 置 或平移来描述。 利用一系列子小波， 一个信号可以在不同的时间尺度上进行计算并显示出详 细的特征尺度。 拉伸更大的小波窗口， 使其宽度更大便可以分析时间系列中波动较大的部分 并捕捉大尺度 （低频） 事件的特征。 相反， 压缩较小的窗口将包含小尺度 （高频） 的事件信 息。 当信号被子小波相乘， 被 s 与 n 唯一的表达， 我们可以计算出信号在时间频率域一个具 体位置的系数。 如果信号 在时间 n 上 的谱成分可以与小波 s 比较， 那么计 算的小波系数具有 相对较大的值。 在其它 n 与 s 的组合 （如其它的子小波） 上都进行这样的计算， 那么将
6、会产 生一系列系数 （小波变换） 来表达信号在时间频率域内的分解。 通过这样的变化便可得到时 间系列的波动模式 （周期变化模式） 以及这些模式随时间的变化 （Furon et al., 2008; Jevrejeva et al., 2003 ） 。 小波变换可以分为连续小波变换 （the Continuous Wavelet Transform, CWT ） 与离散小波 变换（Discrete Wavelet Transform, DWT ） 。离散小 波变换 DWT 是数据的紧凑表示，常用于 降噪与数据压缩。 连续小波变换 CWT 更适合于信号特征的提取 （Grinsted et al.,
7、 2004 ） 。 CWT 作为时间系列间歇式波动特征提取的工具被广泛的应用的地球物理学研究中 （Grinsted et al., 2004; Furon et al., 2008 ） 。 （1 ）连 续小 波变换 CWT 可以将具有等时间步长 t 的离散时间 系列 x n(n=1, N) 的连续小波变换定义为小波函 数 0 尺度化 以及转换下的 x n 的卷积 ： 1 0 * ) ( ) ( N n n X n s t n n x s t s W （1） 式中* 表示共 轭复数，N 是时间系列的总数据个数，( t/s) 1/2 是一个用于小波函数标准化的因 子从而使得小波函数在每个小波尺度
8、s 上具有单位能量。 通过转换小波尺度 s 并沿着时间指 数 n 进行局 部化， 最终可得到一幅展示时间系列在某一尺度上波动特征及其随时间变化的图 谱，即小波功率谱（Torrence and Compo, 1998; Torrence and Webster, 1999; Grinsted et al., 2004 ） 。 对一个时间系列进行小波转换时，母小波的选择显得尤为重要，Farge （1992 ）曾经 讨 论过母小波选择时需要考虑的因素， 例如正交与非正交、 负值与实值、 母小波的宽度与图形 等等。 正交小波函数一般用于离散小波变换， 非正交小波函数即可用于离散小波变换也可用 于连续小
9、波变换（Torrence and Compo, 1998 ） 。通 常在对时间系列进行分析时，希望能够得 到平滑连续的小波振幅， 因此非正交小波函数较为合适。 此外， 要得到时间系列振幅和相位 两方面的信息，就要选择复值小波，因为复值小波具有虚部，可以对相位进行很好的表达 （Torrence and Compo, 1998 ） 。Morlet 小波不但具有非正交性而且还是由 Gaussian 调 节的指 数复值小波。 2 / 4 / 1 0 2 0 ) ( t t i e e t （2） 式中 t 为时间 ， 0 是无量 纲频率。 当 0 =6 ， 小波尺 度 s 与傅里叶周期 （period
10、 ） 基本相等 （ , = 1.03s ） （Torrence and Webster, 1999 ） ，所以尺度 项与周期项可以相互替代。由此可见， Morlet 小波 在时间与频率的局部化之间有着很好的平衡 （Grinsted et al., 2004 ） 。 此外 ， Morlet 小波中还包含着更多的振动信息， 小波功率可以将正、 负峰值包含在一个宽峰之中 （Torrence and Compo, 1998 ） 。 （2 ）小 波功 率谱 为使计 算更 为快捷 ，公 式（1 ） 的卷 积在傅 里叶 域内执 行（Torrence and Compo, 1998; Grinsted et
11、al., 2004 ） 。 2 ) (s W X n 定义为小波功率谱(wavelet power spectrum) ，该功率 谱表达 了时间系列在给定小波尺度和时间域内的波动量级（Lafrenire and Sharp, 2003 ） 。由 于我们 采用的 Morlet 母小波为复 值小波， 因此 ) (s W X x 也为复数， 其复值部分可以解释为局部相位 （Torrence and Compo, 1998 ） 。将小波 功率谱在某一周期上进行时间平均，我们可以得到小 波全谱（global wavelet spectrum ） ， 1 2 2 ) ( 1 ) ( N o n n s W
12、 N s W （3 ） 小波全谱能够表明时间系列真实功率谱的无偏、 一致估计 （Torrence and Compo, 1998 ） 。 由于小波全谱可以显示出背景谱量度，所以局部小波谱的峰值可以得到验证。因为该特性， 通过小波全谱图中可以清晰的辨别时间系列的周期波动特征及其强度。 （3 ）小 波功 率谱边缘 效 应及影响 锥 由于小波变换假设数据是循环的， 所以当我们处理有限长度的时间系列时， 在小波功率 谱中会出现边缘效应， 即在功率谱的起始及末端部分出现误差。 由于该原因， 需要我们在时 间系列的末尾补零从而使得分析的时间系列的总长度 N 大于 2 m 而小于 2 m+1 。 但是， 当
13、我们 采取这样的措施时会在小波功率图谱边缘引起端点不连续以及谱振幅下降的现象。 在这种情况下，需要明确一个概念，即影响锥（Cone of Influence, COI ） ，影响锥 COI 表示小波谱区域以及相应的边缘效应。 在 COI 的 边缘小波谱值会下降 e -2 （Torrence and Compo, 1998; Grinsted et al., 2004; Furon et al., 2008 ） 。 （4 ）小 波功 率谱的显 著 性检验 小波功率谱的统计显著性可以对照一个原假设进行评价， 该原假设为假设信号由一个给 定背景功率谱 （P k ） 的稳定过程产生， 通常背景功率谱为
14、白噪声或红噪声 （Torrence and Compo, 1998; Lafrenire and Sharp, 2003 ） 。 由于 许多地球物理时间系列具有红噪声特征 （即方差随着 尺度的增加 或频率的下 降而增加） ， 所以常采用 红噪声作为 背景谱对小 波谱进行检 验。红噪 声过程可以很好的由一阶自回归过程 （AR1 ） 来 模拟 （Torrence and Compo, 1998; Grinsted et al., 2004 ） 。 一个由 lag-1 自相关 处理 的 AR1 的 傅里叶功率谱可以定义为： 2 2 2 1 1 k i k e P （4 ） 式中 k 为傅 里叶频率指
15、数。 通常在研究中，每个尺度上用 COI 以外的值以 5% 的显著水平进行估计。 （5 ）尺 度选 择 在进行小波变换时， 还需要选择一系列尺度 s 。 本研究使用非正交小波变换， 我们可以 使用任意一组的尺度来构建较完整的图像。一系列尺度可以用 2 的分 数幂来表达： j j j s s 2 0 ，j = 0, 1, J （5 ） ) / ( log 0 2 1 s t N j J （6 ） 式中，s 0 为可分辨的最小尺度，J 为确定的最大尺度。s 0 应该被选择恰当以便使相等的 傅里叶周期近似于 2 t 。 一个足够小 的 j 的选择 依赖于小波方程谱空间的宽度。 4. 小 波分析 步骤
16、 选用开源小波分析软件 （http:/paos.colorado.edu/research/wavelets/），以 47 年年降水数 据为例，操作步骤如下（本实验中使用的数据和部分小波分析代码见附录 A ）： （1 ） 数 据预处 理 。 进行 小波分析的年降水量时序数据必须是连续等时间步长的。 另外， 需要进行标准化处理。标准化处理可选择 SPSS 软件，对降水量时序数据进行 Z 值 标准化， 使其成为无量纲时序数据； （2 ）母小波 选取 。由于小波分析能 够很好的反 应序列的尺 度变化、周 期变化和趋 势变 化， 将其应用在降水序列的分析中， 不仅能反应降水序列的尺度变化， 还可以显示
17、出降水序 列的周期变化以及变化的时间位置。 可选择 Mexican hat 小波 ， 或 Morlet 小波。 通常在对时 间系列进行分析时， 希望能够得到平滑连续的小波振幅， 因此非正交小波函数较为合适。 本 实验中选择 Morlet 小波 ，它比 Fourier 分析更能反映出序列的局部特征。 （3 ）尺度选 择。由于年降水时间序列为 47 年， 时间系列长度 N=47 ，为 了减小功率谱 的边缘效应， 在进行交互小波变换时选择 2 6 个 数据。 时间步长 dt=1 ， 即一年一个数据。 j 可选择 0.125 。 （4 ）显著性 检验 。由于许多地球物 理时间系列 具有红噪声 特征（即
18、方 差随着尺度 的增加或频率的下降而增加），所以常采用红噪声作为背景谱对小波谱进行检验（在程序中 lag1=0.72 ） 。在计算中，每个尺度上用 COI 以外的值以 5% 的显著水平进行估计。 （5 ）小波分 析及绘图。在 MATLAB 软件中运行附录 A 中提 供的小波分析代码即可输 出小波系数、 小波功率、 对应的周期和尺度以及显著性检验结果。 为了更好的表现小波分析 的结果， 将这些数据整理后导入 Surfer 绘图软件中进行小波图绘制。 本实验主要绘制了小波 系数实部等值线图，小波功率谱等值线图和功率谱显著性检验。 5. 年 降水量 的小波分 析 结果 图 3 为标准 化后的刚察气象
19、站 47 年 年降水量序列。标准化公式为年降水量序列减去整 个序列降水量平均值， 再除上降水量时序的标准差。 从图中可以看出， 标准化后的序列平均 值为 0 。值域 在2.5 之间 。 图 4 为年降 水量数据在小波变换后得到的小波系数实部等值线图。图中横坐标为年份， 纵坐标为周期。 当小波系数实部值为正时， 代表降水增多期； 为负时， 代表降水减少期。 从 图中可以清楚看出， 刚察地区年降水量变化过程存在多时间尺度特征。 总的来说， 存在 15-30 年、6-10 年和 2-5 年三类 尺度的周期变化规律。 其中 15-30 年 尺度上出现降水增多和减少交 替的两次震荡，在整个分析时段表现非
20、常稳定；在 6-10 年尺度上出现准 6 次震荡，在 1970-1990 年 间表现稳定。上述两个尺度上的降水量周期变化特征具有全域性。而 2-5 年尺 度的周期变化，在 1980-2000 年之间 表现较为稳定。 图 5 为年降 水量小波功率谱。 图中小波功率越大， 等值线越密集。 倒锥形线为影响锥 （图 5 中短划线） 。 该锥线以下表明该部分小波功率谱图受到边缘效应的影响， 表现出的周期特 征存在较大不确定性。 在图中， 还有三个闭合的等值线 （图 5 中粗线 ） 。 其中最上 部分的两 个等值线环抱的区域对应周期特征为 2-4 年，分 别出现在 1965-1975 年和 1982-20
21、00 年。 表 明， 这两个时间段内， 年降水量的波动以 2-4 年周期 为主。 在周期为 6-10 年， 对应 1970-1984 年之间同样出现等值线环抱的区域。表明该时间段内年降水量波动的平均周期约为 8 年。 图 6 为对小 波谱进行显著性检验的小波全谱图。 从短划线与小波功率谱曲线的关系可以 看出，当短划线小于小波功率谱曲线时，表明该区段对应的周期特征达到了 0.05 水平的显 著性检验。因此，从图中可以看出，通过显著性检验的周期为 2-4 年。 通过小波分析可以看出， 刚察地区年降水量存在 2-4 年为主 的周期变化特征， 总体呈现 年降水量逐年增加的趋势（见图 1 ） 。而从青海
22、湖面积 47 年 的变化特征来看，青海湖面积 几乎呈现单调的减少过程， 没有明显的周期特征 （见图 2 ） 。 为此， 本实验认为， 刚察气象 站观测到的年降水量变化没有反映出降水的年际变化对青海湖面积变化的影响， 与单点气象 资料不能够充分反映整个青海湖地区总体气候变化特征有关， 也表明青海湖面积的变化不是 由某一个因子的独立变化决定的，而是气候变化和人类活动（填湖）两方面共同引起的。 最后， 本实验以降水量的周期分析为例， 阐述小波分析方法在气候变化影响研究中的应 用，同时探讨了年降水量的周期变化特征与青海湖面积变化的关系。 图 3 标准化后的年降水量时序数据 图 4 年降水量数据的小波系
23、数实部等值线图 周期 ( 年) 图 5 年降水量数据的小波功率谱图 图 6 年降水量数据小波全谱图（实线）及 0.05 显著性水平线（短划线） 参考文献 1 Farge M. Wavelet transforms and their applications to turbulence. Annual Review of Fluid Mechanics. 1992. 24: 395-457. 2 Furon A, Wagner Riddle C, Smith C R, et al. Wavelet analysis of wintertime and spring thaw CO2 and N
24、2O fluxes from agricultural fields. Agricultural and Forest Meteorology. 2008. 148, 1305-1317. 3 Grinsted A, Jevrejeva S, Moore J. Application of the cross wavelet transform and wavelet coherence to geophysical time series. Nonlinear Processes in Geophysics. 2004. 11: 561-566. 4 Jevrejeva S, Moore J
25、 C, Grinsted A. Influence of the Arctic oscillation and El Nino -southern oscillation (ENSO) on ice conditions in the Baltic Sea: the wavelet approach. Journal of Geophysical Research. 2003. 108 (D21), 4677. 5 Lafrenire M, Sharp M. Wavelet analysis of inter-annual variability in the runoff regimes o
26、f glacial and nival stream catchments, Bow Lake, Alberta. Hydrological Process. 2003. 17, 1093-1118. 6 Stoy P C, Katual G G, Siqueira M B S, et al. Variablity in net ecosystem exchange from hourly to inter-annual time scale at adjacent pine and hardwood forests: a wavelet analysis. Tree Physiology.
27、2005. 25: 887-902. 7 Torrence C, Compo G P. A practical guide to wavelet analysis. Bulletin of the American Meteorological Society. 1998. 79(1): 61-78. 8 Torrence C, Webster P J. Interdecadal Changes in the ENSO-Monsoon system. 1999. Journal of climate. 1999. 12: 2679-2690. 附录 A 1. 小波分析数据 年份 降水量 （mm
28、 ） 青海湖面积 （km 2 ） 年份 降水量 （mm ） 青海湖面积 （km 2 ） 年份 降水量 （mm ） 青海湖面积 （km 2 ） 1961 762 4505 1981 788 4317 2001 802 4263 1962 766 4490 1982 770 4318 2002 798 4269 1963 778 4476 1983 758 4326 2003 790 4257 1964 767 4475 1984 784 4327 2004 802 4245 1965 773 4472 1985 780 4308 2005 787 4250 1966 771 4457 1986
29、777 4305 2006 803 4284 1967 740 4473 1987 793 4295 2007 787 4291 1968 766 4503 1988 778 4282 1969 802 4485 1989 759 4317 1970 762 4462 1990 798 4339 1971 783 4448 1991 783 4316 1972 784 4451 1992 761 4292 1973 791 4433 1993 762 4292 1974 766 4410 1994 791 4289 1975 758 4410 1995 770 4281 1976 761 44
30、14 1996 782 4271 1977 764 4402 1997 787 4265 1978 787 4392 1998 814 4255 1979 794 4364 1999 790 4259 1980 796 4335 2000 790 4260 2. Matlab 程序 ：wavetest.m %- 参数 定义- s s t = pre cip ita tio n; % p rec ip itat ion 为保存在 MATLAB workspace 中的年降水量矢量数据 varian ce = s td (s s t) 2; % 计算方差 sst = (sst - mean(sst
31、)/sqrt(variance) ; % 数据 标准化。如果导入的降水量数据已经标准化，此步骤可忽略 n = lengt h(s s t); % 降水量数据的矢量长度 dt = 1 ; % 设置时间步长为 1 年 time = 0:n-1*dt+ 1961 ; % 设置时间矢量，从 1961 到 2007 xlim = 1961,2007; % 绘图时 x 轴范围 pad = 1; % 填补模式为 1。本时间序列长度为 47，不足 2 6 ，故选择 0 填补 dj = 0.125 ; % this will do 4 sub-octaves per octave s0 = 2*dt; % 可分
32、辨的最小尺度为 2 年 j1 = 6/dj ; % 为 2 6 中的 6，即以 2 的 幂指数，根据时间序列长度确定 lag1 = 0 .72; % 红噪声背景谱的自相关检验参数 m other = M orlet ; % 母小波为“Morlet ” %- 小波转换- % 计算小波系数、周期、尺度，请参考软件中该函数说明 wave,period,scale,coi = wavelet(sst,dt,pad,dj,s0,j1,mother); pow er = (abs (w a ve). 2 ; % 计算小波功率谱 %- 显著性检验 (variance=1 for the normalized
33、 SST)- % 计算不同尺度的显著性水平，请参考软件中该函数说明 signif,fft_theor = wave_signif(1.0,dt,scale,0,lag1,-1,-1,mother); sig95 = (signif)*(ones(1,n); % expand signif (J+1)x(N) array sig95 = power ./ sig95; % 当比值大于 1，表明强度值显著 %- 小波全谱和显著性检验- global_ws = variance*(sum(power)/n); % 所有时间尺度上的平均 dof = n - s cal e; % 边缘效应的尺度订正 g
34、lobal_signif = wave_signif(variance,dt,scale,1,lag1,-1,dof,mother); %- 保存小波分析结果- save period.dat period -ascii; save scale.dat scale -ascii; save power.dat power -ascii; save global_ws.dat global_ws -ascii; save global_signif.dat global_signif -ascii; %- 绘图- %- 绘制时间序列图- subplot(position,0.1 0.75 0.6
35、5 0.2) plot(time,sst) set(gca,XLim,xlim(:) xlabel(Time (day) ylabel(S(l/s) title(a) 标准化 后的年降水时序 (seasonal) hold off %- 绘制小波功率谱等值线图- subplot(position,0.1 0.37 0.65 0.28) levels = 0.0625,0.125,0.25,0.5,1,2,4,8,16 ; Yticks = 2.(fix(log2(min(period):fix(log2(max(period); contour(time,log2(period),log2(p
36、ower),log2(levels); %* or use contourfill %imagesc(time,log2(period),log2(power); %* uncomment for image plot xlabel( 年份) ylabel( 周期 ( 年) title(b) 年降水 小波功率谱) set(gca,XLim,xlim(:) set(gca,YLim,log2(min(period),max(period), . YDir,reverse, . YTick,log2(Yticks(:), . YTickLabel,Yticks) % Contourfill ima
37、gesc(time,log2(period),log2(power); xlabel( 年份) ylabel( 周期 ( 年) title(b) 年降水 小波功率谱) set(gca,XLim,xlim(:) set(gca,YLim,log2(min(period),max(period), . YDir,reverse, . YTick,log2(Yticks(:), . YTickLabel,Yticks) % 95% significance contour, levels at -99 (fake) and 1 (95% signif) hold on contour(time,lo
38、g2(period),sig95,-99,1,k); hold on % 影响锥线，锥线一下的周期特征存在不确定性 plot(time,log2(coi),k) hold off %- 绘制小波全谱- subplot(position,0.77 0.37 0.2 0.28) plot(global_ws,log2(period) hold on plot(global_signif,log2(period),-) hold off xlabel( 功率 (l/s2) title(c) 小波全 谱) set(gca,YLim,log2(min(period),max(period), . YDi
39、r,reverse, . YTick,log2(Yticks(:), . YTickLabel,) set(gca,XLim,0,1.25*max(global_ws) %- 补充绘制小波全谱图（反向坐标轴）- subplot(position,0.1 0.1 0.2 0.15) plot(log2(period),global_ws) hold on plot(log2(period),global_signif,-) hold off xlabel( 周期 ( 年) ylabel( 功率(1/s2) title(d) 小波全 谱) set(gca,XLim,log2(min(period),max(period),XTick,log2(Yticks(:), . XTickLabel,Yticks)
