小波分析整理

青海湖地区近50年降水量周期变化分析杨沈斌，张弥，吕开龙1．问题青海湖位于青海省东北部的青海湖盆地内，既是中国昀大的内陆湖泊，也是中国昀大的咸水湖。近年来有多篇关于青海湖面积受气候变化影响的研究报告，认为温度升高，降水减少和蒸发量大是造成面积下降的几个重要原因。为此，本实验拟采用连续小波分析方法对青海湖地区年降水量周期变化进行分析，探讨降水量变化与青海湖面积变化的关系。2.资料考虑到刚察气象观测站正好处于青海湖边上，所得数据对反映青海湖周边降水变化具有较好的代表性。因此，以青海湖西北位置的刚察气象台站1961-2007年年降水量资料为例，对该站年降水资料进行小波分析，获取其周期变化特征。同时，获取了青海湖1961-2007年的面积变化资料。从资料发现，青海湖面积近50年呈现下降趋势，平均下降3km2/a。图1和图2分别显示了1961-2007年刚察气象站年降水量和青海湖面积的变化，以及对应的趋势线。图11961-2007年刚察气象台站年降水量及趋势线图21961-2007年青海湖面积及趋势线3.小波分析方法小波变换具有多分辨率分析的特点，并且在时频两域都具有表征信号局部特征的能力。小波变换通过将时间系列分解到时间频率域内，从而得出时间系列的显著的波动模式，即周期变化动态，以及周期变化动态的时间格局（TorrenceandCompo,1998）。小波（Wavelet），即小区域的波，是一种特殊的、长度有限，平均值为零的波形。它有两个特点：一是“小”，二是具有正负交替的“波动性”，即直流分量为零。小波分析是时间（空间）频率的局部化分析，它通过伸缩平移运算对信号（函数）逐步进行多尺度细化，能自动适应时频信号分析的要求，可聚焦到信号的任意细节。小波分析将信号分解成一系列小波函数的叠加，而这些小波函数都是由一个母小波（motherwavelet）函数经过平移与尺度伸缩得来的。用这种不规则的小波函数可以逼近那些非稳态信号中尖锐变化的部分，也可以去逼近离散不连续具有局部特性的信号，从而更为真实的反映原信号在某一时间尺度上的变化。小波分析这种局部分析的特性使其成为对非稳态、不连续时间序列进行量化的一个有效工具（Stoyetal.,2005）。小波是一个具有零均值且可以在频率域与时间域内进行局部化的数学函数（Grinstedetal.,2004）。一个小波被称为母小波（motherwavelet），母小波可沿着时间指数经过平移与尺度伸缩得到一系列子小波。子小波可以通过尺度（s，频率的反函数）函数和时间（n）位置或平移来描述。利用一系列子小波，一个信号可以在不同的时间尺度上进行计算并显示出详细的特征尺度。拉伸更大的小波窗口，使其宽度更大便可以分析时间系列中波动较大的部分并捕捉大尺度（低频）事件的特征。相反，压缩较小的窗口将包含小尺度（高频）的事件信息。当信号被子小波相乘，被s与n唯一的表达，我们可以计算出信号在时间频率域一个具体位置的系数。如果信号在时间n上的谱成分可以与小波s比较，那么计算的小波系数具有相对较大的值。在其它n与s的组合（如其它的子小波）上都进行这样的计算，那么将会产生一系列系数（小波变换）来表达信号在时间频率域内的分解。通过这样的变化便可得到时间系列的波动模式（周期变化模式）以及这些模式随时间的变化（Furonetal.,2008;Jevrejevaetal.,2003）。小波变换可以分为连续小波变换（theContinuousWaveletTransform,CWT）与离散小波变换（DiscreteWaveletTransform,DWT）。离散小波变换DWT是数据的紧凑表示，常用于降噪与数据压缩。连续小波变换CWT更适合于信号特征的提取（Grinstedetal.,2004）。CWT作为时间系列间歇式波动特征提取的工具被广泛的应用的地球物理学研究中（Grinstedetal.,2004;Furonetal.,2008）。（1）连续小波变换CWT可以将具有等时间步长δt的离散时间系列xn(n=1,…,N)的连续小波变换定义为小波函数ψ0尺度化以及转换下的xn的卷积：10'*'')()(NnnXnstnnxstsW（1）式中*表示共轭复数，N是时间系列的总数据个数，(δt/s)1/2是一个用于小波函数标准化的因子从而使得小波函数在每个小波尺度s上具有单位能量。通过转换小波尺度s并沿着时间指数n进行局部化，昀终可得到一幅展示时间系列在某一尺度上波动特征及其随时间变化的图谱，即小波功率谱（TorrenceandCompo,1998;TorrenceandWebster,1999;Grinstedetal.,2004）。对一个时间系列进行小波转换时，母小波的选择显得尤为重要，Farge（1992）曾经讨论过母小波选择时需要考虑的因素，例如正交与非正交、负值与实值、母小波的宽度与图形等等。正交小波函数一般用于离散小波变换，非正交小波函数即可用于离散小波变换也可用于连续小波变换（TorrenceandCompo,1998）。通常在对时间系列进行分析时，希望能够得到平滑连续的小波振幅，因此非正交小波函数较为合适。此外，要得到时间系列振幅和相位两方面的信息，就要选择复值小波，因为复值小波具有虚部，可以对相位进行很好的表达（TorrenceandCompo,1998）。Morlet小波不但具有非正交性而且还是由Gaussian调节的指数复值小波。2/4/1020)(ttieet（2）式中t为时间，ω0是无量纲频率。当ω0=6，小波尺度s与傅里叶周期（period）基本相等（λ,λ=1.03s）（TorrenceandWebster,1999），所以尺度项与周期项可以相互替代。由此可见，Morlet小波在时间与频率的局部化之间有着很好的平衡（Grinstedetal.,2004）。此外，Morlet小波中还包含着更多的振动信息，小波功率可以将正、负峰值包含在一个宽峰之中（TorrenceandCompo,1998）。（2）小波功率谱为使计算更为快捷，公式（1）的卷积在傅里叶域内执行（TorrenceandCompo,1998;Grinstedetal.,2004）。2)(sWXn定义为小波功率谱(waveletpowerspectrum)，该功率谱表达了时间系列在给定小波尺度和时间域内的波动量级（LafrenièreandSharp,2003）。由于我们采用的Morlet母小波为复值小波，因此)(sWXx也为复数，其复值部分可以解释为局部相位（TorrenceandCompo,1998）。将小波功率谱在某一周期上进行时间平均，我们可以得到小波全谱（globalwaveletspectrum），122)(1)(NonnsWNsW（3）小波全谱能够表明时间系列真实功率谱的无偏、一致估计（TorrenceandCompo,1998）。由于小波全谱可以显示出背景谱量度，所以局部小波谱的峰值可以得到验证。因为该特性，通过小波全谱图中可以清晰的辨别时间系列的周期波动特征及其强度。（3）小波功率谱边缘效应及影响锥由于小波变换假设数据是循环的，所以当我们处理有限长度的时间系列时，在小波功率谱中会出现边缘效应，即在功率谱的起始及末端部分出现误差。由于该原因，需要我们在时间系列的末尾补零从而使得分析的时间系列的总长度N大于2m而小于2m+1。但是，当我们采取这样的措施时会在小波功率图谱边缘引起端点不连续以及谱振幅下降的现象。在这种情况下，需要明确一个概念，即影响锥（ConeofInfluence,COI），影响锥COI表示小波谱区域以及相应的边缘效应。在COI的边缘小波谱值会下降e-2（TorrenceandCompo,1998;Grinstedetal.,2004;Furonetal.,2008）。（4）小波功率谱的显著性检验小波功率谱的统计显著性可以对照一个原假设进行评价，该原假设为假设信号由一个给定背景功率谱（Pk）的稳定过程产生，通常背景功率谱为白噪声或红噪声（TorrenceandCompo,1998;LafrenièreandSharp,2003）。由于许多地球物理时间系列具有红噪声特征（即方差随着尺度的增加或频率的下降而增加），所以常采用红噪声作为背景谱对小波谱进行检验。红噪声过程可以很好的由一阶自回归过程（AR1）来模拟（TorrenceandCompo,1998;Grinstedetal.,2004）。一个由lag-1自相关α处理的AR1的傅里叶功率谱可以定义为：22211kikeP（4）式中k为傅里叶频率指数。通常在研究中，每个尺度上用COI以外的值以5%的显著水平进行估计。（5）尺度选择在进行小波变换时，还需要选择一系列尺度s。本研究使用非正交小波变换，我们可以使用任意一组的尺度来构建较完整的图像。一系列尺度可以用2的分数幂来表达：jjjss20，j=0,1,…,J（5）)/(log021stNjJ（6）式中，s0为可分辨的昀小尺度，J为确定的昀大尺度。s0应该被选择恰当以便使相等的傅里叶周期近似于2δt。一个足够小的δj的选择依赖于小波方程谱空间的宽度。4.小波分析步骤选用开源小波分析软件（），以47年年降水数据为例，操作步骤如下（本实验中使用的数据和部分小波分析代码见附录A）：（1）数据预处理。进行小波分析的年降水量时序数据必须是连续等时间步长的。另外，需要进行标准化处理。标准化处理可选择SPSS软件，对降水量时序数据进行Z值标准化，使其成为无量纲时序数据；（2）母小波选取。由于小波分析能够很好的反应序列的尺度变化、周期变化和趋势变化，将其应用在降水序列的分析中，不仅能反应降水序列的尺度变化，还可以显示出降水序列的周期变化以及变化的时间位置。可选择Mexicanhat小波，或Morlet小波。通常在对时间系列进行分析时，希望能够得到平滑连续的小波振幅，因此非正交小波函数较为合适。本实验中选择Morlet小波，它比Fourier分析更能反映出序列的局部特征。（3）尺度选择。由于年降水时间序列为47年，时间系列长度N=47，为了减小功率谱的边缘效应，在进行交互小波变换时选择26个数据。时间步长dt=1，即一年一个数据。δj可选择0.125。（4）显著性检验。由于许多地球物理时间系列具有红噪声特征（即方差随着尺度的增加或频率的下降而增加），所以常采用红噪声作为背景谱对小波谱进行检验（在程序中lag1=0.72）。在计算中，每个尺度上用COI以外的值以5%的显著水平进行估计。（5）小波分析及绘图。在MATLAB软件中运行附录A中提供的小波分析代码即可输出小波系数、小波功率、对应的周期和尺度以及显著性检验结果。为了更好的表现小波分析的结果，将这些数据整理后导入Surfer绘图软件中进行小波图绘制。本实验主要绘制了小波系数实部等值线图，小波功率谱等值线图和功率谱显著性检验。5.年降水量的小波分析结果图3为标准化后的刚察气象站47年年降水量序列。标准化公式为年降水量序列减去整个序列降水量平均值，再除上降水量时序的标准差。从图中可以看出，标准化后的序列平均值为0。值域在±2.5之间。图4为年降水量数据在小波变换后得到的小波系数实部等值线图。图中横坐标为年份，纵坐标为周期。当小波系数实部值为正时，代表降水增多期；为负时，代表降水减少期。从图中可以清楚看出，刚察地区年降水量变化过程存在多时间尺度特征。总的来说，存在15-30年、6-10年和2-5年三类尺度的周期变化规律。其中15-30年尺度上出现降水增多和减少交替的两次震荡，在整个分析时段表现非常稳定；在6-10年尺度上出现准6次震荡，在1970-1990年间表现稳定。上述两个尺度上的降水量周期变化特征具有全域性。而2-5年尺度的周期变化，在1980-2000年之间表现较为稳定。图5为年降水量小波功率谱。图中小波功率越大，等值线越密集