水文要素周期性分析

水文时间系列的周期分析方法探讨

1 前言

水文现象由于受多种气候因素的影响，具有随机性，其时间系列所隐含的震荡周期的识别、判定是一个重要但又较为困难的问题。解决方法从最初的离散周期图、方差周期图、方差分析过渡到谱分析、小波分析，经历了从时间域到时频域的发展过程，特别是1965年傅立叶变换的出现以及上世纪80年代计算机技术的飞速发展，使时频域分析走向实用并迅速拓展。

2 方法原理

2.1 最大熵谱分析方法

谱分析方法把时间系列看成是多种不同频率的规则波（正弦波或余弦波）迭加而成，比较不同频率波的方差大小，从而找出波动的主要周期，对一时间系列的谱分析，主要有功率谱、最大熵谱分析方法等。其中最大熵谱具有突出的分辨能力，峰值偏离小，提取的主次周期更符合实际。

对一组系列值，熵定义为：

ωωd S H ⎰+∞

∞

-=)(ln （1）

式中：ω为频率，)(ωs 为谱密度。

由于功率谱和自相关函数)(j r 互为傅立叶变换，即：

ωωωd e s j r j i ⎰+∞

∞

-=)()( （2）

因为熵值极大输出功率也达到极大，现在要解决的关键问题就是如何利用)(j r 去估计谱密度)(ωs ，使其满足式（2），且式（1）中的熵谱为最大。利用泛函分析中拉格朗日乘子法可以证明，满足以上条件的熵谱密度为：

2

1

)(2

00

1)(∑=--=

k k k

i k k k H e

a S ωσω （3）

式中0k 为自回归的阶数；)

(0k k

a 为自回归系数；2

0k σ为预报误差方差估计。

将计算得到的不同波长（可换算成周期）的最大熵谱密度绘成图，如果谱密度有尖锐的峰点，其对应的周期就是序列存在的显著周期。 2.2 小波分析方法

小波分析是一种时、频多分辨率分析方法，具有时频局部化功能，可以对函数和信号系列进行多尺度细化分析，以分析不同尺度（周期）随时间的演变情况。

小波分析不仅能将水文时间系列的频率特征在时间域上展现出来，分析出其主要周期，而且能清晰地给出各种时间尺度（周期）的强弱和分布情况以及旱涝变化趋势和突变点。方法原理简述如下。

小波函数是指具有震荡特性、能够迅速衰减到零的一类函数，即：

0)(=⎰

-∞

∞

+dt t φ （4）

目前有很多函数可以选用，本文采用“墨西哥帽”小波，)(t φ通过伸缩和平移构成一族函

数系：

)(

)(2

1

,a

b

t a

t b a -=-φφ （5） 式中：R b ∈，R a ∈，0≠a 。)(,t b a φ为连续小波，称为母小波，对于时间系列函数)(t f ，其小波变换的连续形式为：

dt a

b

t t f a

b a W f )(

)(),(2

1

-=⎰

+∞

∞

--φ （6） 其中)(t φ为)(t φ的共轭函数，),(b a W f 称小波变换系数。在实际工作中，时间系列常常是离散的，如t ∆为取样时间间隔，n 为样本量，则式（6）的离散形式为：

)(

)(),(1

2

1

a

b

t i t i f t a

b a W n

i f -∆∆∆=∑=-φ （7） ),(b a W f 随参数a 和b 变化，作出以b 为横坐标，a 为纵坐标的关于),(b a W f 的二维

等值线图，正的小波系数对应于偏多期，负的小波系数对应于偏少期，小波系数为零对应着突变点；小波系数绝对值越大，表明该时间尺度变化越显著。通过小波系数的分析，可识别降水量时间系列多时件尺度演变特性和突变特征。

3 实例分析

3.1 最大熵谱分析方法

根据上述原理，运用最大熵谱研究1956-2000年共45年降水量时间序列，判断1x ，

2x ，…，n x 显著周期。

(1) 建立适当阶数的自回归模型，并根据最终预报误差（FPE ）最小原则，确定模型的最佳阶数0k ，最终预报误差计算见式（8）

2

)(k k

n k n k FPE σ-+=

（8）

(2) 用递推模型计算最终的自回归系数。

∑∑=-=++=n

t t t n

t t t x x

x

x a 2

2122

1

11

)(2

（9）

∑∑∑∑∑∑+===+-----+===+-----++⎥⎦

⎤⎢⎣⎡-+-⎥

⎦

⎤

⎢⎣⎡--=n k t k j k j j k t kj k t j t kj t n k t k j k j j k t kj k t j t kj t k k x a x x a x x a x x a x a 2112

1122111111

)()())((2 （10）

(3) 利用（3）式计算最大熵谱，由于降水系列为离散值，只能采用离散形式计算熵谱值，因此计算出的最大熵谱的离散形式：

2

1

)(2

1)(20

00

0)sin()cos(1)(⎥⎦⎤⎢⎣⎡+⎥⎦⎤⎢⎣⎡-=

∑∑==k k k k

k k k k k H m lk a m lk a l S ππσ （11）

上式中，频率取m

l

l 22πω= ( l =0，1，2，…，m)，m 为选取的最大波数，在序列样本量不大时，m 通常取

2

n

。 根据上述原理，最终计算出最佳阶数为6以及不同周期的最大熵谱值，如表1，并绘制最大熵谱图，见图1。

表1 年降水最大熵谱表