两种常用的水文气象谱分析方法的对比研究

两种常用的水文气象谱分析方法的对比研究

0水文过程谱分析一般分为三大环光谱分析是一种通用的数据分析方法，具有很强的适用性。它通常用于分析数据系统的物理属性、周期成分或过滤噪声。频谱分析大体经历了调和分析(傅立叶级数展开,也称谐波分析)、功率谱分析(自相关函数的傅立叶变换)和最大熵谱分析(简称MEM)等发展过程。一个连续的时间序列在某种意义上构成了一种波,地下水水文过程线可以认为由若干个不同频率、振幅的简单波形组成。频谱分析是把降水量等复杂的组分(或波形),分解成单纯的成分(波形),从而了解其固有的性质(如周期性、滞后性等),并给出定量的结果。频谱分析是在滞后域和频率域上研究序列的性质,这一点不同于在时间域上的分析。降水预测分析中经常遇到确定降水周期的问题,本文分别采用功率谱分析和最大熵谱分析方法对挠力河流域菜咀子站1964~2003年的降水量资料进行周期分析,并比较了两种方法的优缺点。1,[n评分设某一实测气象要素时间序列为Xt(t=1,2,…,n),其距平序列为xt。其中xt=Xt-ˉXtxt=Xt−X¯¯¯t。采用功率谱分析和最大熵谱分析对xt进行周期分析。1.1功率谱的确定功率谱分析是以傅立叶变换为基础的频域分析方法,将时间序列的总能量分解为不同频率上的分量,根据不同频率波的方差贡献判断序列隐含的显著周期。功率谱是应用极为广泛的一种周期分析方法。根据功率谱密度与自相关函数互为傅立叶变换的性质,可以通过自相关函数间接进行功率谱估计。(1)最大时滞m的确定对距平序列xt,最大时滞为m的自相关系数r(j)(j=0,1,2,…,m)为r(j)=ln-jn-jΣt=1(xt-ˉxS)(xt+j-ˉxS)(1)r(j)=ln−jΣt=1n−j(xt−x¯S)(xt+j−x¯S)(1)式中:ˉxx¯、S分别为序列的均值和标准差,其表达式如下:ˉx=1nnΣt=1xt,S=√1nnΣt=1(xt-ˉx)2计算自相关系数首先要确定最大时滞m。如果m取值过大,谱的峰值就多,会影响自相关函数和谱估计的稳定性。当m太小时,谱估计过于平滑,不易出现峰值,难以确定主要周期,而且得不到低频部分(长周期)的谱估计。因此m的选取十分重要,一般取为n/10~n/3。(2)计算粗谱估计值不同波数k(k=0,1,2,…,m)的粗谱估计值为:ˆˆSk=1m[r(0)+2m-1Σj=1r(j)coskπjm+r(m)coskπ](2)(3)工艺平滑谱值的确定由式(2)得到的谱估计与真实谱存在一定的误差,因此需要进行平滑处理,消除粗谱估计的小波动。常用Hanning平滑系数进行平滑,即{S0=0.5ˆS0+0.5ˆS1Sk=0.25ˆSk-1+0.5ˆSk+0.25ˆSk+1Sm=0.5ˆSm-1+0.5ˆSm(3)另外,在平滑过程中,由于端点S0及Sm仅是两点等权重平滑,为保持权重的一致,对(3)式的平滑谱值在端点上乘以1/2。(4)k轴级谱图以波数k为横轴,平滑功率谱估计值Sk为纵轴作图。在谱图横轴上,同时标出对应的周期或频率。它们与波数k的关系为ωk=πk/m,Τk=2m/k(5)标准谱和白噪音检验为了确定谱值在哪一个波段最突出并了解该谱值的统计意义,需要求出一个标准过程谱以便比较。标准谱有两种情况:①红噪音标准谱:S0k=ˉS[1-r(1)21+r(1)2+2r(1)cos(πk/m)](4)其中,ˉS为m+1个谱估计的均值。②白噪音标准谱:S0k=ˉS(5)用红噪音还是白噪音标准谱进行检验,Rao等人给出一个判据:r0=-1+tα√n-2n-1(6)在信度α=0.05下,tɑ取1.645。若r(1)>r0,则认为序列存在明显的持续性,可使用红噪音检验;否则,就用白噪音检验。计算S′0k=S0k(X2α/v)(7)其中,X2α为给定显著性水平α条件时自由度为v的X2分布,v=(2n-m2)/m(8)若谱估计值Sk>S′0k,则表明波数k对应的周期波动是显著的。1.2最大熵谱的计算由功率谱分析的计算过程可以看出,它不能兼顾高频和低频段的需要;某些短周期振动易在一些周期长度为它们整数倍的长周期中表现出来,又混在长周期中;所取样本较短时,得出的周期可能与实际有偏离。因此,1967年Burg提出了最大熵谱估计方法,以克服上述困难。该方法具有突出的分辨能力,峰值偏离小,提取的主次周期更符合实际。最大熵谱分析主要用Burg算法,其思路是,建立适当阶数的自回归模型,计算最大熵谱。设变量x的自回归模型为:xt=a1xt-1+a2xt-2+⋯+akxt-k+εt假设均值为0,k阶预报误差的方差为σ2k(k=1,2,…,n-1),相应的系数为ak,1,ak,2,…,ak,k。它们可由(9)式的推导公式求得。{ak,k=2nΣt=k+1(xt-k-1Σj=1ak-1,jxt-j)(xt-k-k-1Σj=1ak-1,jxt-k+j)nΣt=k+1[(xt-k-1Σj=1ak-1,jxt-j)2+(xt-k-k-1Σj=1ak-1,jxt-k+j)2]ak,j=ak-1,j-ak,kak-1,k-j,σ2k=(1-a2k,k)σ2k-1(9)其中,a1,1=2n-1Σt=1xtxt+1/n-1Σt=1(x2t+x2t+1),σ20=1nnΣt=1x2t。由上述的递推过程可以看出,Burg算法的优点在于直接用时间序列来计算谱估计中的参数,不必提前算出自相关系数,克服了功率谱估计中自相关系数最大时滞主观选择的不足。最大熵谱为SE(l)=σ2k0[1-k0Σk=1ak0,kcos(πlkm)]2+[k0Σk=1ak0,ksin(πlkm)]2(10)其中,k0为自回归模型的阶数,可由最终预测误差(FPE)准则确定。即使FΡE(k)=n+kn-kσ2k(11)最小时的k值就是自回归模型的最佳阶数。2年降水量资料对计算过程编制VB程序,以挠力河流域菜咀子站的年降水量资料(图1)为例,来说明功率谱分析及最大熵谱分析的应用。对实测数据取距平值,并进行计算。2.1标准曲线的绘制取m=n/3=13,计算得到的自相关系数、粗谱估计值及平滑谱估计值见表1,功率谱图见图2。由式(6)得r0=0.234>r(1),因此用白色噪音检验。由式(5)得S0k=0.071,由式(8)得v=5.65,取显著性水平α=0.05,查X2分布表可得X2α=12.06,由式(7)得S′0k=0.152。在图2中画出白噪音标准谱直线。由图2可以看出,k=8时的Sk值为一峰值,且明显超过标准谱,故第一个显著性周期为3.25年。另外,k=3时的Sk值为第二个峰值,比较接近标准谱,故另一个显著性周期为8.67年。因此,菜咀子站年降水量存在3年左右和8~9年的周期振荡。2.2最大熵谱图取m=n/2=20,计算得到的年降水量距平序列的最大熵谱见表2。最大熵谱图见图3。由图3可以看出,图上有两个明显的峰值点,最高峰值对应的周期为3.6年,次峰值对应的周期为8.0年。因此菜咀子站年降水量主要存在3.6年和8.0年的周期。2.3线准三年周期的配方优化由功率谱及最大熵谱估计可知,挠力河流域降水量存在3~4年和8~9年的水文周期。3~4年的周期与徐国昌、李栋梁等得到的降水周期一致,它的形成可能与南亚高压脊线准三年周期的南北振荡有关。8~9年的周期符合我国北方的降水周期规律。对比图2和图3可知,功率谱分析得到的谱图峰值不太明显,k=2、3时的Sk值很接近,不容易区分显著周期;最大熵谱分析得到的谱图峰值非常明显,分辨率很高。这是因为功率谱估计的谱分布在很大程度上取决于自相关系数的最大时滞,同时它假定序列的自相关系数在最大时滞以外为零,导致了功率谱估计的偏移,而最大熵谱分析没有这些不符合实际情况的假定,因此在分辨率、谱偏移和对短序列的适应性等方面都优于功率谱估计。3最大熵谱的特点(1)频谱分析是研究时间序列周期性的重要手段,它能给出定量分析结果,揭示时间序列的内部