小波变换与信号瞬时特征分析

小波变换与信号瞬时特征分析

1 实信号的小波变换及其性质

20()(),)()0()()()()1(,)()()()()2()0i b R R i b R g t b s t L R dt g t b Ra g e g t g t g a

t b S b a e g a s d g s t dt a

g C d ωωωωωπωωω

∞-∞∞-∈∈<-=≠=<∞??信号()()f t L R ∈的Fourier 变换()f ω的定义为

()()i t f f t e dt ωω-∞-∞=?，

反变换公式为

1()()2i t f t f e d ωωωπ∞-∞

=?， (,)L R dt 表示绝对可积函数空间，R 为实数集合。

1.1 小波变换

考虑“解析小波” ()g t （即()0g ω=，当0ω<时），满足

12()(,d )(,)g t L R t L R dt ∈ ，（1）

且

12()(\{0},)(\{0},)||||

d d g L R L R ωωωωω∈ ，（2） 任给一个信号

2()(,)s t L R dt ∈，()s t 相对于()g t 的小波变换定义为 1(,)()()t b S b a g s t dt a a

∞-∞-=?，（3） 这里，t R ∈，0a <，b R ∈，2()L R 为平方可积函数空间，(

)t b g a -表示()t b g a -的复共轭。

(,)S b a 的频率域计算公式为

1

(,)()()()2i b t b S b a e g a s d g a

ωωωωπ∞

-∞-=?。（4） 1.2 利用小波变换计算实信号对应的解析信号

定理1 设()g t 是满足（1）、（2）式的解析小波，()g t 的实部()R g t 为偶函数，()R g t 的Fourier 变换记为()R g ω，()R g ω满足

()

R g g C d ωωω∞=<∞?，0g C ≠，（5）

任给一实函数2()(,)s b L R dt ∈，则 01

(,)()[()]g da S b a s b iH s b C a

∞=+?。（6） 这里，(,)S b a 由（3）式定义，[()]H s b 表示()s b 的Hilbert 变换。

证：由（3）式得

1(,)()()1()[()()]11()()()()(,)(,)

R I R I R I t b S b a s t g dt a a

t b t b s t g ig dt a a a

t b t b s t g dt i s t g dt a a a a

S b a iS b a ∞-∞∞-∞∞∞-∞-∞-=--=---=-=+????（7） （7）式中，()R t b g a -及()I t b g a -依次表示()t b g a -的实部和虚部， 1(,)()()R R t b S b a s t g dt a a

∞-∞-=?，（8a ） 1(,)()()I I t b S b a s t g dt a a

∞-∞-=-?，（8b ） 利用定理1，对任给的尺度因子a （0a >），

(,)(,)(,)

(,)[(,)]

11()()[()()R I R R R R S b a S b a iS b a S b a iH S b a t b t b s t g dt iH s t g dt a a a a

∞∞-∞-∞=+=+--=+??，（9） （8）式中两边同乘以

1g C a ，然后对a 在区间(0,)∞作积分，得 2002020(,)1

11()()1

1()()()1

11[()]()2R R g g

R g i t R g S b a t b da s t g dtda C a C a a t b s b s t g dtda C a a

b t s e d g dtda C a a

ωωωπ∞∞∞-∞∞∞-∞∞∞∞

-∞-∞-=-=-=???????? 上式中，先对t 积分，作变量替换并利用g C 的定义，有

00(,)()1

11()21()()

2i b R R g g

i b S b a g y da dy s e d C a C y s e d s b ωωωωπωωπ∞∞∞-∞∞

-∞===????

即2011()()()R g t b s b s t g dtda C a a

∞∞-∞-=??。（10） 上式中y R ∈，（9）式两边同乘以1g C a

，然后对变量a 在区间(0,)∞作积分，注意到Hilbert 变换为线性变换，再利用（10）式，即可得（6）式。

1.3 瞬时参数的定义

瞬时参数的定义与通常由解析信号求瞬时参数定义同。下面列出瞬时频率、瞬时相位及瞬时振幅等参数的计算公式。

1.3.1 不同尺度下的瞬时参数

根据定理1，对任给的尺度因子a ，可求出该尺度下的瞬时参数。

(,)e t a =（11a ）

(,)arctan[(,)/(,)]I R t a S t a S t a θ=，

（11b ） (,)arctan[(,)/(,)]I R d t a S t a S t a dt

ω=，（11c ） 式中，(,)e t a 、(,)t a θ及(,)t a ω依次表示尺度为a 时信号的瞬时振幅、瞬时相位及瞬时频率。

1.3.2 实信号()s t 对应的瞬时参数

()e t =（12a ）

*()()arctan ()

s t t s t θ=，（12b ） *()()[arctan ]()

d s t t dt s t ω=，（12c ） 上式中，*()s t 为()s t 对应的解析信号虚部，

*01()[()]Im (,)g da s t H s t S t a C a

∞==?，（13） ()e t 、()t θ及()t ω依次表示()s t 的瞬时振幅、瞬时相位及瞬时频率。

实际中常用的还有加权平均频率及阻尼瞬时频率。阻尼瞬时频率定义为

**222max

()()()()11()()22()ds t ds t s t s t dt dt f t t e t e ωππε-==+?，（14）

这里，()e t 为()s t 的瞬时振幅，22max max(())e e t =，（15）

ε为小于1的正实数。

2 不同尺度瞬时参数的含义

地震资料分析中一个常用模型是褶积模型。尽管实际地震记录满足“变子波模型”， 但在一定条件下可用褶积模型近似。

用()t ω、()R t 及()s t 依次记地震波、反射系数序列及地震记录。根据根据地震记录的褶积模型。

()()()()()s t R t d R t d ωττττωττ∞∞

-∞-∞=-=-??，（16） 把（16）式代入（3）中得

1(,)()()1[()()]()()()t b S b a s t g dt a a

t b R t d g dt a a S b R d ωτωτττττ∞-∞∞∞-∞-∞∞

-∞-=

-=-=-????，（17） 式中，1(,)()()1()()()(,)(,)

R I t b S b a t g dt a a t b t g dt a a

S b a iS b a ωωωτωττωττ∞-∞∞-∞--=

---==-+-??，（18） 根据定义，(,)S b a ωτ-为地震子波()t ω的小波变换。（18）式中，

1()(,)()()R R t b S b a t g dt a a

ωττω∞-∞---=

?，（19） 1()(,)()()I I t b S b a t g dt a a ωττω∞-∞---=?，（20） 根据定理1有

(,)[(,)]I R S b a H S b a ωωττ-=-，

设()g t 满足定理2，则()()R R g t g t -=，（19）式可写成

1()(,)()()1()(,)R R R b t S b a t g dt a a

g b a a ωττωωτ∞-∞---==*-?，（21）

式中，R g ω*表示ω与R g 的褶积。由（21）式知，给定a ，(,)R S b a ωτ-实际上是

用()R g t 作为滤波器对地震子波作带通滤波的结果。如果把(,)S a ωτ看成新的震源

子波，根据（17）式，(,)R S b a 就是对应新震源子波的地震记录，(,)I S b a 是(,)R S b a 的Hilbert 变换。(,)S b a 是与新地震记录对应的解析信号。利用(,)S b a 计算出的瞬时参数，描述了新地震记录的特征。

根据小波变换的性质，当a 变化时，()(

)R b t g a τ--的频带宽度及中心频率都将发生变化。A 减小时，()()R b t g a

τ--中心频率升高，频带变宽，由（21）式知，(,)S a ωτ的分辨能力变高；当a 增大时，()()R b t g a

τ--中心频率变低，频带变窄，由（21）式知，(,)R S a ωτ的分辨能力变低。

综上所述，当a 变化时，(,)R S b a 为一系列具有不同分辨率的地震记，(,)S b a 是与其对应的解析信号。由(,)S b a 计算的瞬时参数，描述了具有不同分辨率的地震记录特征。这就是不同尺度下计算出的瞬时参数的物理意义。

3 基于（6）式计算瞬时参数

设20.005ε=，Hilbert 滤波因子取801个点。

图1 无噪声时两种方法计算结果比较

图2 有噪声时Hilbert变换方法的计算结果

当不考虑噪声时，只要Hilbert 滤波因子取得足够长，用Hilbert 变换方法与由（6）式计算的解析信号及瞬时参数吻合得很好。图1a 是一道合成地震记录。子波采用Ricker 子波。图1b 中示出了两种方法计算结果，图1c 为两种方法计算结果相减取绝对值。均方误差为0.007。

噪声对Hilbert 变换方法计算结果影响很大。例如，图2a 是图1a 中加上白噪声的结果。图2b 是用Hilbert 变换方法求出的瞬时频率。可见，有噪声时直接利用Hilbert 变换方法计算出的瞬时参数，不能很好地反映有效信号的特征。

当信号中含有噪声时，如何最大限度地压制噪声，而又使有效信号畸变最小，这一直是重要的研究课题。一种常用方法是：用双边滤波器对信号作滤波以压制噪声，然后用Hilbert 变换对滤波结果求瞬时参数。图2c 是对图2a 作带通滤波的结果。图2d 是用Hilbert 变换方法求出的图2c 的瞬时频率，比较图2d 与2b ，可见经滤波后求出的瞬时频率，精度明显地提高。

对信号()s t 作小波分解时，当所选小波形状与有效信号比较接近时，有效信号能量分布在时间一尺度域一个小的闭子空间V 中，而干扰波及随机噪声能量分布在时间一尺度域的另一个大的闭子空间V1。V 与V1可能完全分离，也可能部分重叠，当我们在V 空间讨论问题时，噪声得到了一定的压制，当V 与V1有重叠时，V 中各分量的信噪比可能是不同的。

在用（6）式求瞬时参数时，我们可以选实部与有效信号最接近的解析小波，然后把信号()s t 分解在时间一尺度域，确定有效信号的能量分布空间V 。在V 中，如果噪声出现在某个频道可以先对此频道作处理，然后利用（6）式在处理后的V 空间中计算瞬时参数。这种方法在压噪时，引起有效信号的失真最小。

下面用上述小波方法计算图2a 中信号的瞬时频率。图3a 是对图2a 作时间-尺度域分析的结果。所用小波函数如（10）式取4τ=，6m =，图3b 是在时间-尺度域压噪后，利用（6）式计算的瞬时频率。

比较图3b 与图2d ，小波变换方法更接近无噪声时的结果。小波变换方法的均方误差为7.20，而滤波加Hilbert 变换方法的均方误差为14.15。可见前者优于后者。

图3 用小波方法计算图2a 的瞬时频率

图4 复杂模型计算结果比较

图4a为一地质模型，它由一水平层和一楔形模型构成。图4b为该模型的自激自收合成记录，其中水平层和楔形模形中地震子波依次为50和30Hz的Ricker 子波。图4c是加上白噪声后的结果。图4d、e依次是用带通滤波加Hilbert变换及小波方法计算出楔形模型平分线上的瞬时频率。两图中的细实线是无噪声时的瞬时频率。由图可见，小波方法计算结果明显地优于Hilbert变换加滤波方法。

4 结论

本文讨论了用小波变换分析信号的瞬时特征及求取瞬时参数问题。文中给出的定理2揭示了实信号相对于解析小波的小波变换与它的Hilbert变换之间的联系。模型举例说明了基于定理2计算瞬时参数的有效性。由于小波变换有良好的时间一尺度域局部化性质，因此，根据有效信号与干扰波时间一尺度域特征差别，在实际问题中，可望构造出抗干扰能力强、计算精度高的瞬时参数计算方法。