第6章 气象上常用小波及其应用

第6章 气象上常用小波及其应用实例(1)

前面五章讲述了小波分析方法的由来和原理，这些基本知识为气象上实际应用奠定了基础。本章将介绍气象上常用的几种小波，特别是Haar 小波和墨西哥帽（Mexihat ）小波，以及小波分析的应用实例。

6.1 二进小波

二进小波的产生基于第4章的“二分法”。它的基本思路是把连续型函数)(t f 及其连续小波变换),(b a W f 离散化，以便于实际应用。作为一种方便和常用的形式，是对小波参数中的放（伸）缩因子a 进行二进制离散。若小波函数系的表达式

{}

Z

m,n n t a

m

∈-- ),(ψ

（??

） 中的放缩因子)(2Z j a j ∈=，则称)(t ψ为二进小波。把经过这种离散化后的二进小波的变换，

称为二进小波变换。

强调说明：在应用时所用的小波函数系（式（??））与前面第4章第3节的式（?）有所不同。比照这两式：

),()(),2(2

)(,2

/,b at a t k t t b a j

j k

j -=

-=ψψψψ

或

（?）

),()(,n t a

t m

n m -=-ψψ

（??

） 可以看出，二者主要的不同点是t 的系数a 指数正负号恰好相反。所以，用式（?）的小波作变换时，随着a （或者j ）的增大，W 曲线变窄；而用式（??）的小波作变换时，随着a （或者j ）的增大，W 曲线变宽。

定义6.1 函数)()(2R L t ∈ψ被称为二进小波，若存在两个常数∞<≤

∑

∈-≤≤

Z

j j B A .)2(?2

ωψ

（6.1）

上述条件式（6.1）称为稳定性条件；若A =B ，则称最稳定条件。而函数序列{}Z

j j f W ∈2称为二进小波变换，其中，

,d 2)(2

1)()(2

2b b t t f t f t f W R

j

j

j

j ?

??

?

??-=

*=ψψ

（6.2）

这里j a 2=是放缩因子，b 是平移因子。 由卷积定理知：

)2(?)(?)(?2

ωψωωj f f W j ?= （6.3） 据此，稳定性条件等价于：对任意)()(2R L t f ∈，有：

∑

∈≤≤

Z

j f

B f

W f

A j 20

20

220

（6.4）

下面定理说明，二进小波一定是允许小波。

定理6.1 设)(t ψ是二进小波，则它一定是允许小波，且

.2ln d )(?,

d )(?2ln 0

2

2

B A ≤-≤

?

?

∞

∞

ωω

ωψ

ωω

ωψ

（6.5）

若A =B ，则上式变成：

.2ln 2d )(?2

A C ==

?

∞

+∞

-ωω

ωψ

ψ （6.6）

其中A 、B 的含义同式（6.1）。

式（6.5）的证明如下：用变量替换，可得

ω

ω

ωψ

ωω

ωψ

d )(?d )2(?1

2

2

2

2

1

2

?

?

+=

j j

j

因此，用ω除式（6.1）中每一项并在区间（1，2）上对ω积分，有

.2ln d )(?2ln 0

2

B A ≤≤

?

∞

ωω

ωψ

同样，用-ω除式（6.1）每一项并在区间（-2，-1）上对ω积分，得

.2ln d )(?2ln 0

2

B A ≤-≤

?

∞

+ωω

ωψ

证毕。 定理6.1表明，当ψ使得式（6.1）成立时，由二进小波变换可完全重构原信号。重构公式是：

.d 2)(2

1)()()( 22

2b b t t f W t t f W t f R j j

j j

j ?

??

?

??-?=

*=ψψ

（6.7）

二进小波介于连续小波和离散小波之间，它只是对尺度参数进行离散化，而在时域上仍

保持平移量（b ）连续变化。因此二进小波变换仍具有连续小波变换的平移不变性。这是它有别于离散小波变换的独特优点，后者不具备平移不变性质。

6.2 正交小波中常用的Haar 小波

正交小波是一种重要的小波型式，正交小波变换已在第4章讲过。在正交小波中，Haar 小波是一种常用的型式。Haar 小波型式简单，使用方便，计算快捷；然而它的缺点也很明显，即缺乏连续性。本节将详细介绍Haar 小波，由此可以加深对正交小波的理解。

图6.1 Haar 小波

(a) 三次样条函数)(t θ； (b) )(t θ的一阶导数，即小波函数t t d /d )(θψ=。

(b) 虚线是小波函数的“简化”，即Haar 小波。

图6.1(a)是作为光滑函数的三次样条函数)(t θ；（b ）是由t d /d θ构成的小波函数，其中

虚线意为实曲线的“简化”，即Haar 小波。

也可以从尺度函数与小波函数关系的角度来理解Haar 小波的构造，也就是该小波是如何构造出来的。这时将信号序列{}4321,,,x x x x 看成单位区间上的一个函数：

),

()()()()()1,4/3[4)4/3,2/1[3)2/1,4/1[2)4/1,0[1t X x t X x t X x t X x t f +++= （6.8）

其中)(),[t X b a 表示区间),[b a 上的特征函数。)(t f 是一个分段常数函数，如图6.2所示。

).

4/3()(),2/1()(),4/1()()4/1,0[)1,4/3[)4/1,0[)4/3,2/1[)4/1,0[)2/1,4/1[-=-=-=t X t X t X t X t X t X

图6.2 序列{}4321,,,x x x x 对应的分段常数函数

现分析)()4/1,0[t X 、)()2/1,4/1[t X 、)()4/3,2/1[t X 和)()1,4/3[t X 之间的关系。由平移关系

知：

)4/1()()4/1,0[)2/1,4/1[-=t X t X ,

)2/1()()4/1,0[)4/3,2/1[-=t X t X ,

)4/3()()4/1,0[)1,4/3[-=t X t X .

另一方面，)()4/1,0[t X 可以看成单位区间上的特征函数的伸缩，即

)2()(2

)1,0[)4/1,0[t X t X =

于是，)()2/1,4/1[t X 、)()4/3,2/1[t X 和)()1,4/3[t X 都可看成)()1,0[t X 的伸缩和平移。

引入记号： ),()()1,0[t X t =?

（6.9）

定义：

)12

,,1,0(),2()(,-=-=j

j k

j k k t t ??

（6.10）

则有 )()(0,0t t ??=

??

?<≤==.

0,1/2,t 0 ,1)2()(0,1其它t t ??

??

?<≤=-=.

0,1,t 1/2 ,1)12()(1,1其它t t ??

这表明，)()(110,1t t ??和可以由)(0,0t ?的伸缩与平移得到。如图6.3所示。

图6.3 尺度函数的图形

称)(t ?为Haar 尺度函数。于是，式（6.8）定义的函数)(t f 可以重新写成如下形式：

).()()()()(3,242,231,220,21t x t x t x t x t f ????+++= （6.11）

即函数f (t )可以用一个尺度函数)(t ?的伸缩和平移的线性组合表示。

k j ,?

在其上不为零的区间叫做该函数的支撑。可以验证，k j ,?

的支撑为??

??

??+j j

k k 21,

2

，支

撑的宽度j

2

1随j 的增加而减小。j 称为分辨率，

j

2

1称为尺度。

下面讨论用{}4321,,,x x x x 的小波变换{}1,10,10,00.0,,,d d d a 来表达式（6.11）定义的函数)

(t f 。把式（6.11）看成)(t f 最高分辨级的分解。求平均如同把)(t f 进行较低分辨级的表示，

这等价于用0,1?和1,1?来表示)(t f 。在这种表示中，0,1?的系数是0,1a ，它是式（6.11）中)(t f 的前两个系数的平均值；1,1?的系数是1,1a ，它是)(t f 的后两个系数的平均值。于是得到如下的函数：

)()()(1,11,10,10,11t a t a t g ??+= （6.12）

一般地，21x x ≠，所以由1)8/1(x f =及2/)()8/1(210,11x x a g +==知，1g 不等于1f ，即用1g 来表示f 时丢失了一些细节信息，因而需要一个函数来描述这种信息。该函数称为小波函数。 基本小波函数定义如下：

??

?

??<≤<≤=-=

. 0,0,1/2 1,-1/2,

0 ,1 )()()()1,2/1[)2/1,0[其它t t t X t X t ψ （6.13）

于是，)12()2()(--=t t t ??ψ。)(t ψ称为Haar 小波，如图6.4所示。 与尺度函数)(t ψ类似，定义小波函数的伸缩和平移如下：

.12

,,1,0),2()(,-=-=j

j k

j k k t t ψψ

（6.14）

于是有

),()(0,0t t ψψ=

??

?

??<≤<≤==

. 0,1/2,1/4 1,-1/4,

0 ,1 )2()(0,1其它t t t t ψψ

??

?

??<≤<≤=-=

. 0,1,3/4 1,-3/4,1/2 ,1 )12()(1,1其它t t t t ψψ

图6.4 Haar 小波)(t ψ

图6.5 )(0,1t ψ和)(1,1t ψ的图形表示

)(0,1t ψ和)(1,1t ψ的图形如图6.5所示。

现在重新考虑利用{}1,10,10,00.0,,,d d d a 、尺度函数{}1,10,1,??和小波函数{}1,10,1,ψψ表示)

(t f 的问题。由式（6.8）知，)(t f 的最高分辨率表示是在长度为1/4的子区间上。下面就

在每个子区间)1,4/3[)4/3,2/1[),2/1,4/1[),4/1,0[和上来看一下f g 和1的差别，以便利用小波函数捕捉丢失的细节信息。 在区间)4/1,0[上，);()()(0,10,110,10,111d a x d a x t g t f +=?=-=- 在区间)2/1,4/1[上，)()()(0,10,120,10,121d a x d a x t g t f -=?-=-=-。 这表明，在区间)2/1,0[上，)()()(0,10,10,10,1t d t a t f ψ?+=。 类似地，在区间)1,2/1[上，有)()()(1,11,11,11,1t d t a t f ψ?+=。

将区间)2/1,0[及)1,2/1[上两个)(t f 的表达式写在一起，可以得到区间)1,0[上)(t f 的一

个新的表达

形式：

).()()()()(1,11,10,10,11,11,10,10,1t d t d t a t a t f ψψ??+++= （6.15）

)

(t f

被表示成支撑宽度为1/2的尺度函数与小波函数的线性和。其中，)(0,10,1t d ψ和

)(1,11,1t d ψ用于捕捉用)(1t g 逼近)(t f 时丢失的细节信息。

现在根据对{}1,10.1,a a 求平均的原理，引入区间)1,0[上的平均函数)(0,00,0t a ?代替式（6.15）中的前两项，而细节项保持不变。令

),()()()(1,11,10,10,10,00,00t d t d t a t g ψψ?++=

考虑用)()(0t f t g 逼近时，如何捕捉丢失的细节信息。这时，

在区间)4/1,0[上，;)()(0,00,10,010d d a x t g t f =--=- 在区间)2/1,4/1[上，0,00,10,020)()(d d a x t g t f =+-=-。 类似地，在区间)1,2/1[上，有0,00)()(d t g t f -=-。

因此，在区间)1,0[上，

).()()()()(1,11,10,10,10,00,00,00,0t d t d t d t a t f ψψψ?+++= （6.16）

式（6.16）表明，函数)(t f 可用序列{}4321,,,x x x x 的小波变换、尺度函数和小波函数表示。

综上所述，给出了序列{}4321,,,x x x x 对应的函数表示，以及该函数的如下多分辨表示：

).

()()()( )()()()( )()()()( )()()()()(1,11,10,10,10

,00,00,00,01,11,10,10,11,11,10,10,13,242,231,220,21)1,4/3[4)4/3,2/1[3)2/1,4/1[2)4/1,0[1t d t d t d t a t d t d t a t a t x t x t x t x t X x t X x t X x t X x t f ψψψ

?ψψ??????+++=+++=+++=+++= （6.17）

其中，)()(),()(0,00,0t t t t ψψ??==。而)()(t t ψ?和满足以下方程：

),()()(1,10,1t t t ???+= （6.18） ).()()(1,10,1t t t ??ψ-= （6.19）

称式（6.18）为双尺度方程；称式（6.19）为小波方程。 把本节的内容和第3章中的“二分法”内容加以对照，可以加深对尺度函数和小波函数的理解。

[例] 将信号{}3,1,2,4-=f 看成区间??

?

??

?

41,0上取值4，??

?

???21,41上取值-2，??

?

???43,21上取值1，

??

????1,43上取值3的分段常数函数，即2V f ∈，且)(3)()(2)(4)(3,22,21,20,2t t t t t f ????++-=。

问f 的小波变换是什么？函数)(t f 的多分辨表示是什么？在同一坐标系中画出)(t f 的多分辨表示的图形。

用塔式算法（见本节后面的《附录》）可以算出，{}3,1,2,4-=f 的小波变换为

?

?

????--1,3,21,23，于是，函数)(t f =的多分辨表示为：

).

()(3)(2

1

)(23 )(3)()(2)(4)(1,10,10,00,03,22,21,20,2t t t t t t t t t f ψψψ?????-+-=

++-= （6.20）

在同一坐标系中

)(t f 的多分辨表示的图形如图6.6所示。

其中，细实线表示原信号)(t f ，

粗实线表示总平均值)(2

3

0,0t ?，点划线表示低分辨率细节)

(2

10

,0t ψ

-

，虚线表示高分辨率细

节)()(31,10,1t t ψψ-。

图6.6 同一坐标系中)(t f 的多分辨表示的图形

《附录》 平均与细节

以下通过求平均与细节这两个线性变换，介绍信号序列的多分辨表示、小波变换及塔式算法，同时说明将信号的小波变换用于信号压缩的思想。

设{}21,x x 是一个由两元素组成的信号，定义这两个元素的平均与细节为：

,

2/)(,2/)(2121x x d x x a -=+=

则可将{}d a ,作为原信号的另一种表示，而且原信号{}21,x x 可由{}d a ,恢复如下：

.

,21d a x d a x -=+=

信号{}d a ,是{}21,x x 的线性变换的结果。可以看出，当21x x 与非常接近时，d 会很小，

这时，{}21,x x 可近似地用 {}a 来表示，由此可实现信号压缩。重构信号为{}a a ,，误差信号为{}{}d d a x a x ,,21=--。特别地，当21x x =时，0=d ，这时可实现信号的完全重构。这表明，平均值a 可以看成原信号的整体信息，而d 可看成原信号用a 表示时丢失的细节信息。用{}a 近似地表示{}21,x x 可实现信号压缩，而且丢失很小的细节信息对最终信号的重构不会造成大的影响。

对于较多元素的信号{}4321,,,x x x x ，通过如下的平均和细节运算：

平均：

2

/)(2/)(431,1210,1x x a x x a +=+=

（6.21）

细节：

2

/)(2/)(431,1210,1x x d x x d -=-=

（6.22）

可得到信号的另一表示{}1,10,11,10,1,,,d d a a 。若细节信号1,10,1,d d 都很小，则丢失细节信号，可得压缩后的信号为{}1,10,1,a a 。 用类似的方法可对压缩后的信号{}1,10,1,a a 进行压缩：

:

:细节平均

?????-=+=2

/)(2/)(1,10,10,01,10,10,0a a d a a a

（6.23）

如果0,01,10,1,,d d d 都非常小，则可用{}0,0a 代替原信号{}4321,,,x x x x 。事实上，

,

4/)( 2/)2/)(2/)((2/)(432143211,10,10,0x x x x x x x x a a a +++=+++=+= （6.24）

即0,0a 为整个信号所有元素的总平均值，它保留了信号最基本的信息。

另一方面，由{}0,00,0,d a 可重构{}1,10,1,a a ，由{}1,10,1,a a 和{}1,10,1,d d 可重构{}4321,,,x x x x ，于是，由{}1,10,10,00,0,,,d d d a 可恢复{}4321,,,x x x x ，这就找到了原信号序列的另一种表示：

{}1,10,10,00,0,,,d d d a （6.25）

式（6.25）所表示的序列称为原序列{}4321,,,x x x x 的小波变换。这种表示为数据压缩提供了更多的选择。

为体现小波分析中多分辨率表示的思想，称原信号{}4321,,,x x x x 为最高分辨率信息，而根据求平均与细节的不同层次将平均信息（也称低频信息）与细节信息（也称高频信息）与分辨率联系起来，具体如下：

{}4321,,,x x x x ——最高分辨率信息；

{}1,10,1,a a ——次高分辨率低频信息； {}1,10,1,d d ——次高分辨率细节信息； {}0,0a ——最低分辨率低频信息； {}0,0d ——最低分辨率细节信息。

在这个意义上，{}4321,,,x x x x 的小波变换{}1,10,10,00,0,,,d d d a 由信号的整体平均及两个不同

分辨率级的细节信息组成，因而是{}4321,,,x x x x 的一个多分辨率表示。

信号的小波变换过程可以用塔式算法表示。以{}4,2,1,3-为例，它的塔式算法如图6.7所示。

图6.7 小波变换的塔式算法

在所有正交小波中，Haar 小波具有最短的支集，且是奇对称的。因为它仅有一阶消失矩，所以不太适合用来逼近光滑函数。归结起来： Haar 小波为：

??

?

??<≤-<≤=. ,01,t 1/2 ,1,

2/10 ,1 )(其它t t ψ

（6.26）

Haar 小波的频域函数为：

).

4/(sin

4

)(?2

2

/ωω

ωψ

ωi e

i -= （6.27）

Harr 小波的图形见图6.4。

Haar 小波的频域函数图见图6.8。

图6.8 Harr 小波的频域函数图（取模）

6.3 Morlet 小波

Morlet 小波为： ..)(0

2

2/t i t e e t ωψ-=

（6.28）

Morlet 小波的频域函数为： .2)(?2/)(20

ωωπωψ--=e

（6.29）

Morlet 小波的图形见图6.9.

Morlet 小波的频域函数图见图6.10。

Morlet 小波是一个常用的小波，因为它的时-频局部性能比较好（尽管在严格意义上不是紧支撑的）。此外，该函数不满足允许条件，但当50≥ω时，允许条件近似成立。Morlet 小波不存在尺度函数。

图6.9 Morlet 小波

图6.10 Morlet 小波的频域函数的图形

6.4 Gauss 小波

这是Gauss 函数的一阶导数，在信号与图像边缘提取中有着重要的应用。 Gauss 小波函数为：

.21)(2

/2

t te

t --

=π

ψ （6.30）

Gauss 小波的频域函数为：

.)(?2

/2

ω

ωωψ

-=e i （6.31）

可见，Gauss 小波函数及其傅氏变换（频域函数）具有同一函数形式（相差一个常数）。

Gauss 小波及其傅氏变换（频域函数）的图形见图6.11。

图6.11 Gauss 小波及其频域函数的图形

6.5 Marr 小波（也称Mexihat 小波）

它是Gauss 函数的二阶导数（但差个负号）形成的小波。

Mexihat 小波函数为：

.)1(32)(2

/2

2

t e

t t --=

π

ψ （6.32）

图6.12 Mexihat 小波及其频域函数的图形

Mexihat 小波的频域函数为：

.3

22)(?2

/242

ω

ωπ

ωψ

-=e

（6.33）

Mexihat 小波及其频域函数的图形见图6.12。小波函数的形状像墨西哥帽，因而得名。

Mexihat 小波在时域和频域具有很好的局部化，并且满足允许条件。Mexihat 小波和Haar 小波都在气象上常用。Mexihat 小波在计算机视觉领域有重要的应用。使用时，可将式（6.32）和式（6.33）前面的常系数加以简化，如式（6.32）的系数改为π

21

，并常取33≤≤-t 。

随着小波分析方法的推广应用，将会用到更多的小波和小波变换。应用于气象观测所获

得的离散有限时间序列（序列长度为N ）二进型小波变换公式是：

,2)(2

11),(1∑

=??

?

??-?=

N

n j

j

b

b n n f

c b j W ψ 其中 .,,2,1;,,1,0,)(?22

1

n b m j c b ===∑

∞

=ω

ωψ

πω （6.34）

小波逆变换（重构）公式是：

()

∑

∑==??

?

??-?=

n

j j

n

b j

b

b n b j W a

c t f 1

1

2

/3.2),(2)

0.2log(1)(ψ （6.35） 利用式（6.35），取不同的j 值（j =0,1,…,m , m 为某个正的定值），可在不同尺度j 2下进行平滑，以及重构出较原序列)(t f 平滑的（派生）序列)(2t S j

。

需要强调指出：使用Haar 小波变换结果作原序列突变点分析时，注意它的极值点：极大值点对应原序列自下向上的突变点；极小值点则对应自上向下的突变点。使用Mexihat 小波变换时注意它的过零点：自上向下过零点对应原序列自上向下的突变点；反之亦然。 以下介绍小波分析在气象时间序列中的应用实例。