小波分析 - 图文

西安石油大学本科毕业设计（论文）

小波多分辨率分析在地震资料处理中的应用

摘 要： 小波分析是近年来发展起来的新的数学理论和方法，在噪声消除方面有

着广泛的应用。小波分析能同时在时频域内对信号进行分析，所以它能有效区分信号中的突变部分和噪声，从而实现对非平稳信号的消噪。小波多分辨率分析是利用小波变换对信号进行多尺度分解分析。利用小波多分辨率分析方法对地震记录进行去噪，可以有效地对噪声和信号进行分离，从而使地质信息丰富清晰，便于解释。利用此方法去噪的过程，可分为如下三个步骤：对记录进行小波多尺度分解，选择一个小波并确定分解的层次，然后进行分解计算；小波分解高频系数的阈值量化，对各个分解尺度下的高频系数选择一个阈值进行阈值量化处理；小波重构，根据小波分解的低频系数和各层高频系数进行小波重构。这三个步骤中，最关键的是如何选择阈值以及进行阈值量化。在某种程度上，它关系到信号消噪的质量。

关键词: 小波变换；多分辨率分析；去噪；阈值；地质信息提取

- I -

西安石油大学本科毕业设计（论文）

The wavelet multi-resolution analysis used in the processing of

seimic data

Abstract: The wavelet analysis is a newly mathematic theory and method that has

developed recently. And specially ,it is widely used in de-noising. The wavelet analysis can simultaneously analyze signals in both time-domain and frequency-domain, so it can distinguish sudden changes of the signals from noises efficiently and be used to remove noises from nonstationary signals. The wavelet multi-resolution analysis uses the wavelet analysis to multi-scale decomposition for signals. It can be used to de-noise seismograms and can separate noises from signals efficiently，therefore the geological information can be made more plentiful and cleaning.The profile can be easily interpretated. There three steps to use this method to de-noise. First, do multi-scale decomposition for the signal. Select a wavelet and determine the scale of the decomposition, then decompose. Second, fix a threshold for high frequency coefficients under each scale and de-noise. Third, do wavelet reconstruction,according to the low frequency coefficients and the de-noised high frequency coefficients. The most important things in the three steps are the selection of thresholds and de-noising using the thresholds.To some degree, it has an effect on the quality of the result of de-noising.

Key words: Wavelet transform； Multi-resolution analysis；De-noising；

Threshold；The extract of geological information

- II -

西安石油大学本科毕业设计（论文）

目 录

1 绪论…………………………………………………………………………………….1 1.1 课题背景及意义…………………………………………………………………1 1.2 国内外研究现状…………………………………………………………………1 1.3 课题研究方法……………………………………………………………………2 2 小波变换基本理论…………………………………………………………………….3

2.1 Fourier变换……………………………………………………………………...3 2.2 短时Fourier变换………………………………………………………..............3 2.3 小波变换…………………………………………………………………………4 2.3.1 连续小波变换…………………………………………………………...5 2.3.2 二进小波变换…………………………………………………………...5 2.3.3 离散小波变换…………………………………………………………...5 3.1 多分辨率分析概念的引入………………………………………………………7 3.2 小波函数与小波空间……………………………………………………………7 3.3 二维多分辨率分析………………………………………………………………8 4.1 小波域阈值滤波算法…………………………………………………………...10 4.2 阈值函数的选取………………………………………………………………...12 4.3 阈值确定方法…………………………………………………………………...12 5.1 一维离散小波变换的Matlab函数……………………………………………..15

5.2 二维离散小波变换的Matlab函数……………………………………………...17 5.3 Matlab信号去噪………………………………………………………………...20 5.4 Matlab图像去噪………………………………………………………………...20

3 多分辨率分析………………………………………………………………………….7

4 小波域阈值滤波………………………………………………………………………10

5 小波域阈值滤波的Matlab实现……………………………………………………...15

6 小波多分辨率分析对地震记录去噪与地质信息提取的Matlab实现……………...22 6.1 程序编写…………………………………………………………………….......22 6.2 模型去噪与地质信息提取……………………………………………………...22 6.3 实际资料去噪与地质信息提取……………………………………...................24 7 结论与认识……………………………………………………………………………48 参考文献…………………………………………………………………………………..49 致谢………………………………………………………………………………………..50 附录1……………………………………………………………………………………...51 附录2……………………………………………………………………………………...56

- III -

西安石油大学本科毕业设计（论文）

1 绪 论

1.1 课题背景及目的

小波变换是20世纪80年代中期发展起来的一种时频分析方法，比离散傅里叶变换的性能更优越，被广泛应用于调和分析、话音处理、图像分割、石油勘探和雷达探测等方面，也被应用于音频、图像和视频的压缩编码。

分析平稳信号的理想工具是傅里叶变换。对于非平稳信号，傅里叶变换不再是有效地工具没因为其无法描述信号的局部频率特性。而小波变换正是分析非平稳信号的有力工具。小波变换时傅里叶变换的新发展，它既保留了傅里叶变换的优点，又弥补了傅里叶变换在信号分析上的一些不足。在传统的基于傅里叶变换的信号处理方法中，要使信号和噪声的频带重叠部分尽可能地小，这样，在频域就可以通过时不变滤波方法将信号同噪声分开。而当它们的频谱重叠时，这种方法就无能为力了。基于小波变换的非线性滤波是完全不同的，在这种方法中，谱可以重叠，但是谱的幅度（而不是谱的位置）要尽可能不同。在小波变换域，可通过对小波系数进行切削、缩小幅度等非线性处理，以达到滤除噪声的目的。小波的多分辨率特性就是在不同尺度上描述信号的局部特征，如边缘、尖峰、断点等，处理的实质在于减小甚至完全剔除由噪声产生的系数，同时最大限度地保留真实信号的系数，最后由经过处理的小波系数重构源信号，得到真实信号的最优估计。

1.2 国内外研究现状

小波变换是近十几年信号处理领域研究的一个热点，许多学者将小波仔理论上的研究成果应用到诸如图像压缩、特征提取、信号滤波和数据融合等方面，而且小波变换的应用领域还在不断发展当中。小波之所以在信号处理领域具有很大的优势，在于小波变换可以获得信号的多分辨率描述。

近年来，小波滤波这一概念不断见之于有关信号及图像处理研究的文献中，这标志着一种新的信号滤波思想的出现。在早期的多尺度信号处理工作中，人们就已注意到信号和噪声在不同尺度上有不同的特征表现，并试图有效地利用这些特征，小波变换的出现为这一思想提供了一个自然而完美的工具，使信号图像的多尺度处理技术得到迅速发展。

1989年，Mallat提出了多分辨率分析的概念，从空间的概念上形象的说明了小波的多分辨率分析特性，并将在此以前各种正交小波基的构造方法统一起来，给出了小波变换的快速算法，即Mallat算法。

总而言之，近年来有关小波的多分辨率分析滤波的文献很多，并取得了不少的成果。

- 1 -

西安石油大学本科毕业设计（论文）

1.3 课题研究方法

本课题的研究分三个阶段，第一个阶段为理论分析阶段，学习小波变换的原理，掌握小波多分辨率分析的概念和原理，以及应用；学习Matlab。第二个阶段为程序设计、模型验证以及实际资料处理阶段，利用所掌握的小波多分辨率分析和Matlab知识，编写程序，对所给地震资料模型和实际资料进行处理，验证模型，并对实际资料去噪。第三个阶段是进行对比和总结，得出结论。

- 2 -

西安石油大学本科毕业设计（论文）

2 小波变换基本理论

2.1 Fourier变换

Fourier变换由下列公式定义： 正变换

逆变换

1 f(t)?2??f(?)??f(t)e?j?tdt (2-1)

????????f(?)ej?tdt (2-2)

对于确定信号和平稳随机信号，傅里叶变换时信号分析和信号处理技术的理论基础，有着非凡的意义，起着重大作用。

傅里叶变换把时间域与频率域联系起来，f(?)具有明确的物理含义，通过研究

??f(?)来研究f(t)，许多在时间域内难以看清的问题，在频率中往往表现得非常清楚。

但正是由于傅里叶变换的域变换特性，f(t)与f(?)彼此之间是整体刻画，不能够反映各自在局部区域上的特征，因此不能用于局部分析。作为变换积分核的e?j?t的幅值在任何情况下均为1，即e??j?t?=1，因此，频谱f(?)在任一频率处的值是由时间

?过程f(t)在整个时间域（??，?）上的贡献决定的；反之，过程f(t)在某一时刻的状态也是由f(?)在整个频率域（??，?）上的贡献决定的。如果要知道所分析的信号在突变时刻的频率成分，那么傅里叶变换是无能为力的，因为傅里叶变换的积分作用平滑了非平稳信号的突变成分。

简言之，傅里叶变换能提取出函数在整个频率轴上的频率信息，却不能反映信号在局部时间范围内的特征。然而，对于变频信号如音乐、地震信号、雷达回波等。此时所关心的恰恰是信号在局部时间范围（特别是突变时刻）内的信号特征（一般是频率成分）。例如，在音乐和语音信号中，人们关心的是在什么位置出现什么样的反射波。

对非平稳信号用傅里叶变换进行分析，不能提供完全的信息，也即通过傅里叶变换可以知道信号所含有的频率信息，但无法知道这些频率信息究竟出现在哪些时间段上。可见，若要提取局部时间段（或瞬间）的频率特征信号，傅里叶变换已经不再适用了。

2.2 短时Fourier变换

为了克服傅里叶分析的局限性，使其对非平稳信号也能作较好的分析，通过对信号在时域上加一个窗函数g(t-?),使其对信号f(t)进行乘积运算以实现在?附近的开窗盒平移，再对加窗的信号进行傅里叶分析，这就是短时傅里叶变换，或者称为

- 3 -

西安石油大学本科毕业设计（论文）

加窗傅里叶变换。短时傅里叶变换定义如下：

sf(?,?)??f(t)g(t??)e?j?tdt (2-3)

???其中，窗口函数f(t)一般取为光滑的低通函数，保证g(t??)只在?点附近局部地测量了频率分量?的幅度值，得到信号在t=?时刻附近的频率信息。这是时间频率局部化的一种标准技术。若采用高斯函数作为窗口函数，其相应的傅里叶变换仍旧是高斯函数，从而保证短时傅里叶变换在时域和频域内均有局部化功能。

如果选取的窗口函数在时域和频域都具有良好的局部性质（如呈指数衰减的高斯函数），此时短时傅里叶变换能够同时在频域和时域内提取关于信号的精确信息。

但短时傅里叶变换存在其固有的局限，其时间频率窗口是固定不变的，一旦窗口函数g(t)选定，其时频分辨率也就确定了，并且不随频率?和时间?而变化。也就是说，它对所有的频率都使用同样的窗口。我们若想提高时间分辨率，就要把窗口缩的很窄，但这样势必会降低频率分辨率。由测不准原理可知，不可能在时间和频率上均有任意高的分辨率，因为时间和频率的最高分辨率受下式的制约：

1????t? (2-4)

4?式中?t和??分别代表时间域和频率域的窗口宽度。这表明，任一方分辨率的提高都意味着另一方分辨率的降低。

2.3 小波变换

上节的的分析表明，短时傅里叶变换的问题的症结在于使用了固定的窗口，而对实际时变信号的分析需要时频窗口具有自适应性：对于高频谱的信息，时间间隔要相对地小以给出较高的精度；对于低频谱的信息，时间间隔要相对地宽以给出完全的信息。换句话说，重要的是要有一个灵活可变的时间-频率窗，能够在“高中心频率”时自动变窄，而在“低中心频率”时自动变宽。而小波函数就是为此而设计的。

2（R）设f(t)是平方可积函数，记作f(t)?L，?(t)为基本小波或母小波，它一般是

时域上以t=0为中心的带通函数，在时域和频域都具有局部化（紧支撑），且均值为零，即

?如果?(t)满足容许性条件

????(t)dt?0 (2-5)

?2 C???则

??(?)?d??? (2-6)

0 Wf(a,b)?f,?a,b??f(t)???1a?(t?b)dt (2-7) a称为f(t)的小波变换。其中

- 4 -

西安石油大学本科毕业设计（论文）

?a,b= a为尺度因子，b为位移因子。

小波变换的特点有：

1t?b?(), a>0, b?R (2-8)

aa(1) 时频局域性、多分辨分析、数学显微镜； (2) 自适应窗口滤波：低频宽、高频窄； (3) 适用于去噪、滤波、边缘检测等。

2.3.1 连续小波变换

对任意连续函数或信号f(t)进行小波变换，如果基函数?a,b(t)的两个参数a和b均为连续变量，则被称为连续变量，则被称为连续小波变换(continuous wavelet transform,CWT)。

连续小波变换的形式可以写为

Wf(a,b)?????f(t)1a?(t?b)dt (2-9) a2.3.2 二进小波变换

为了将小波变换应用于信号分析的实践，需要将小波变换进行离散化，就像将傅里叶变换离散化一样。将小波变换离散化就是对变换参数进行离散化。

本节先讨论尺度因子a的离散化问题。在小波变换中，常令a取离散值 aj?2j (j?Z) (2-10) 称之为二进离散化。

这时，小波基函数的形式为

?2,b(t)?2?j/2?(jt?b) (2-11) j2对应的小波变换为

Wf(2j,b)??f(t)???12j?(t?b)dt (2-12) j2j称为分辨率级别。

2.3.3 离散小波变换

下面再将二进小波变换中的平移因子也离散化。令b=n2j，则可得离散小波变换：

?1t?n2jjj Wf(2,n2)??f(t)?()dt (2-13) j??j22 对小波变换进行离散化的一般情形是小波基函数?a,b的尺度因子a和位移因子

b都只限定在某些离散点上取值。一种最通常的离散方法如下：尺度因子按幂级数进

- 5 -

西安石油大学本科毕业设计（论文）

行离散化，位移因子在尺度内均匀离散化，在尺度间具有幂次关系，即有： a?a0j， b?nb0a0j；a0?1,b0?0, j,n?Z。 (2-14) 此时的小波基函数表示为如下形式：

t?nb0a0j1?j/2?j?()?a?(at?nb0) (2-15) ?j,n(t)?00jaj0a0任意函数f(t)离散小波变换为

Wf(j,n)?f(t),?j,n(t) ?????f(t)?j,n(t)dt

?a?j/20???f(t)?(a?j?0t?nb0)dt (

- 6 -

2-16)

西安石油大学本科毕业设计（论文）

3 多分辨率分析

3.1 多分辨率分析概念的引入

1989年，Mallat提出了多分辨率分析的概念，从空间的概念上形象的说明了小波的多分辨率分析特性。若把尺度理解为照相机的镜头的话，当尺度由大到小变化时，就相当于将照相机镜头由远及近地观察目标。在大尺度空间里，对应远镜头下观察到的目标，只能看到目标大致的概貌；在小尺度空间里，对应近镜头下观察目标，可观察到目标的细节部分。这种由粗及精对事物的分析就称为多分辨率分析。

空间L2(R)的一个多分辨率分析是指满足以下性质的一个闭子空间序列

Vj?L2(R)，j?Z。

（1） 一致单调性：

Vj?Vj?1，对任意j?Z (3-1) （2） 伸缩规则性：

f(t)?Vj?f(2t)?Vj?1 (3-2) （3） 平移不变性：

f(t)?V0?f(t?n)?V0 ，对所有n?Z (3-3) （4） 逼近性：

?V??0?,?Vjj?Zj?Zj?L2(R) (3-4)

（5） 正交基存在性： 存在一积分值非零的函数?(t)?V0，使??(t?n)n?Z?是

V0的标准正交基。

称?(t)为多分辨率分析的尺度函数，Vj为尺度j上的尺度空间。由多分辨率分析的定义知，所有的闭子空间Vj的尺度空间。

??j?Z都是由同一函数?(t)伸缩后的平移系列张成

3.2 小波函数与小波空间

记Vj在Vj?1的补空间为Wj，即

Vj?1?Vj?Wj，Wj?Vj (3-5) 显然，任意子空间Wm和Wn（m?n）是相互正交的，并且，由式（3-1）和式（3-4）可知

?Wj?L2(R) (3-6)

j因此，Wj???j2是的一系列正交子空间。若，显然有f(t)?WL(R)f(2t)?Wj。 0j?Z- 7 -

西安石油大学本科毕业设计（论文）

与尺度函数的定义类似，称?为小波函数，相应地，称Wj是尺度为j的小波空间。显然，小波空间是两个相邻尺度空间的差，也就是说，空间Wj包含了函数f投影到尺度空间Vj与Vj?1间的细节差别，因此小波空间又称为细节空间。

由上的分析可知，多分辨率分析的核心是V0，W0空间的标准正交基?(t?n)以及

?(t?n),n?Z。只要它们已知，分析就可以逐级进行。而要找到这两组标准正交基，其关键是找到合适的尺度函数?(t)和小波函数?(t)。

3.3 二维多分辨率分析

构造多维小波基最简单且常用的方法是张量积方法。应用这种方法可将一维多分辨率分析很容易地推广到二维多分辨率分析。L2(R2)的多分辨率分析是L2(R2)的子空间序列。设Vj2(j?Z)是L2(R2)的一个子空间，信号f(x,y)在分辨率2?j的近似等于信号在向量空间VJ2上的正交投影。类似于一维多分辨率分析，在二维情况下存在着二维尺度函数?(x,y)，使其伸缩和平移??j,m,n(x,y)?构成二维多分辨率分析的标准正交基。其中

?j,m,n(x,y)?2?j?(2?jx?n,2?jy?m),j,m,n?Z (3-7)

设Vj2?Vj?Vj，“?”表示张量积，则Vj2为L2(R2)的多分辨率分析的充要条件是Vjj?Z(j?Z)为L2(R)的一个多分辨率分析。这里只讨论可分离的二维小波，设 ?(x,y)??(x)?(y) (3-8) 设Vj2(j?Z)是L2(R2)的可分离多分辨率分析，?(x,y)??(x)?(y)是相应地二维尺度函数，?(x)是一维多分辨率分析中与?(x)对应的小波函数，则

?1(x,y)??(x)?(y) (3-9) ?2(x,y)??(x)?(y) (3-10) ?3(x,y)??(x)?(y) (3-11) 的平移和伸缩?2?j?i(2?jx?n,2?jy?m)j,m,n?Z,i?1,2,3?构成L2(R2)的标准正交基。记

?j1?j?j ?12?(2x?n,2y?m) (3-12) (x,y)?j,m,n ?j,m,n(x,y)?2 ?j,m,n(x,y)?232?j?2(2?jx?n,2?jy?m) (3-13)

?j?3(2?jx?n,2?jy?m) (3-14)

对二维图像f(x,y)?L2(R2)，可定义其离散小波变换为

d23 (AJf,(D1jf)1?j?J,(Djf)1?j?J,(Djf)1?j?J) 其中

df?(f(x,y),?j,m,n)m,n?Z AJ- 8 -

西安石油大学本科毕业设计（论文）

1 D1jf?(f(x,y),?j,m,n)m,n?Z (3-15) 2f?(f(x,y),? D2j,m,n)m,n?Z j3 D3jf?(f(x,y),?j,m,n)m,n?Z

(3-16) (3-17)

d23(AJf,(D1jf)1?j?J,(Djf)1?j?J,(Djf)1?j?J)构成了信号f(x,y)的二维正交分解。其

2中系数D1jf给出了f(x,y)垂直方向的小波系数，Djf给出了f(x,y)水平方向的小波系

数，D3jf给出了f(x,y)对角线方向的小波系数。

- 9 -

西安石油大学本科毕业设计（论文）

4 小波域阈值滤波

4.1 小波域阈值滤波算法

小波域阈值滤波是现在研究的一个热点。小波滤波方法对噪声性质有三个基本假设：

（1） 噪声经小波变换后大多数小波系数为零或近似为零； （2） 噪声均匀地分布在所有系数中； （3） 噪声水平不是太高。

假设观测数据

fi?gi??i,i?1,2,...,N(N?2n) (4-1) 由真实信号gi和加性噪声?i组成，其向量形式为f=g+?。在理想情况下，?i为 （1） 正态噪声； （2） 不相关的噪声； （3） 方差为常量。

当然实际情况并非如此，因此每种假设条件都可以适当放宽，以满足实际应用的需要。我们的目标是在观察到f得前提下，对g进行估计。

小波的时频局部化特性和多分辨率特性决定了小波滤波方法与传统方法相比，具有独特的优势—能够在去处噪声的同时，很好的保留信号的突变部分或图像的边缘。

信号和噪声在小波域中有不同的性态表现，它们的小波系数幅值随尺度变化的趋势不同。随着尺度的增加，噪声系数的幅值很快衰减为零，而真实信号系数的幅值基本不变。小波滤波就是利用具体问题的先验知识，根据信号和噪声的小波系数在不同尺度上具有不同性质的机理，构造相应地原则，在小波域采用其他数学方法对含噪信号的小波系数进行处理。处理的实质在于减小甚至完全剔除由噪声产生的系数，同时最大限度的保留真实信号的系数，最后由经过处理的小波系数重构原信号，得到真实信号的最优估计。

小波滤波的特点可概括问哦：首先，它不是平滑。平滑是去除高频信息而保留低频信息；而小波滤波是要试图去除所有噪声，保留所有信号，并不考虑它们的频率范围。其次，它是在小波变换域对小波系数进行非线性处理；第三，滤波过程一般由三个步骤完成：

（1） 小波变换；

（2） 对小波系数进行非线性处理，以滤除噪声； （3） 小波逆变换。

设对感测数据经过小波变换后，得

- 10 -

西安石油大学本科毕业设计（论文）

w= ??? (4-2) W(?)和W?1(?)分别表示小波变换和逆变换算子，前述滤波过程的三个步骤可描述为

?w=W（f）,wt?D(w,t),g?W?1(wt) (4-3)

其中，D(?,?)为非线性滤波算子，它是滤波问题的核心。显然，这种原理性的归纳并没有涉及到算子W(?)或D(?,?)是如何作用到信号上的，也没有涉及到它们的选取方法以及对滤波效果的影响。通过选择不同的W(?)和D(?,?)，可以得到多种不同的小波滤波方法。

小波滤波的核心是在上述第二个步骤中按照一定得准则对小波系数进行修改，以在不损失过多信号的前提下，达到降低或去除噪声的目的。

对式（4-1）所示的观测数据滤波，Donoho提出了两个滤波的前提条件：

?（1） 光滑性：在大概率条件下，g至少跟g有同样的光滑度；

?（2） 适应性：g是最小均方差估计。

由条件（1）可以推出，当N??时，下式几乎接近于1的概率成立：

? gF?C1gF (4-4) 式中C1为一常数。这在小波域中意味着

? ?j,i??j,i (4-5) 成立，式中1?i?N为位置，j为尺度； ?j,i表示真实信号在尺度j上的第i个小波系数。

?1N?2对于条件（2），可以理解为E（?gi?gi）求最小值，着等价于求E???FNi?1??的最小值。Donoho证明，此时g、?必须满足

? Eg?g2F???E???2F???N2lnN (4-6)

当小波变换为正交小波变换时，??1。由式（4-1-6）可知

?21 E???F???2lnN (4-7)

N?由式（4-7）可知，对于任何?j,i???2lnN，取?j,i?0，将满足式（4-6），这相当于认为当?j,i???2lnN时，?j,i由噪声所产生，因此，可取阈值为

t??2lnN (4-8)

由此可得硬阈值法滤波的步骤为： （1） 对信号求小波变换；

（2） 除了最粗尺度信号外，将各细节信号作阈值处理，阈值t取为??2lnN，当

某位置小波变换值大于阈值时，保留原值，否则置零，用公式表示为

- 11 -

西安石油大学本科毕业设计（论文）

?j,i???0,?j,i?t (4-9) ?????j,i,其他（3） 利用小波变换重构，求出信号的滤波值。

而软阈值滤波的步骤为： （1） 对信号求小波变换；

（2） 除了最粗尺度信号外，将各细节信号作阈值处理，阈值t取为??2lnN，当某位置小波变换值大于阈值时，向着减小系数幅值的方向作一个收缩t，否则置零，用公式表示为

? ?j,i?sgn(?j,i)?(?j,i?t) 其中sgn（x）为符号函数；

（3） 进行小波变换重构，求出信号的滤波值。

至此，只要知道噪声的方差，就可以实现整个滤波的算法了。Donoho提出了噪声标准差?的估计。至于软阈值法和硬阈值法，哪种方法的滤波效果更好，应视具体情况而定。

以上的算法是基于正交小波得到的，在实际应用中，可将其推广到非正交小波的情形。信号经非正交小波基分解后，不同尺度上的系数彼此相关，故阈值应与尺度有关。设噪声传播到尺度j上的噪声标准差为?j，则式（4-7）应为

?21E???F??j?2lnN (4-11) N这时尺度j上的阈值为tj??j2lnN。

小波域阈值滤波算法中的两个基本要素是阈值和阈值函数。

(4-10)

4.2 阈值函数的选取

阈值函数体现了对小波系数的不同处理策略，主要分为软阈值函数、硬阈值函

数和半阈值函数。它们的基本思想都是去除小幅值的系数，对幅值较大的系数进行收缩或保留。

硬阈值法往往使得滤波结果具有较大的方差（主要因为滤波后的不连续性），而软阈值滤波结果有较大的偏差（主要因为其对所有大于阈值的系数共同作了收缩）。总的说来，硬阈值法可以很好的保留信号或图像的边缘等局部特征，但滤波结果会出现伪Gibbs现象。为了克服软阈值法和硬阈值法的缺点，Gao Hong-Ye提出了另一种阈值函数，即半软阈值函数，并在此基础上推导出了基于半软阈值法的Minimax阈值。这种方法是软阈值法和硬阈值法的折衷形式。它不仅保留了较大的系数，而且具有连续性。然而在半软阈值法中，需要确定两个阈值，增加了算法的复杂度。

4.3 阈值确定方法

- 12 -

西安石油大学本科毕业设计（论文）

在小波域阈值滤波中，阈值的选取直接影响滤波效果。因此，如何确定阈值是一个非常关键的问题。目前主要的阈值确定方法有通用阈值法、极小化风险阈值法、假设检验法等。

（1） 通用阈值法（universa method）

这是Donoho等在高斯噪声模型下，应用多维独立正态变量决策理论，得出的阈值。他们发现当维数趋于无穷时，噪声系数的幅值大于阈值 t??2lnN的概率趋于零，其中?为噪声标准差，N为信号长度。通用阈值具有渐进的Minimax特性，然而由于侧重考虑滤波结果的平滑性，结果表现出较大的偏差。实际应用时需估计噪声方差。在软阈值法中应用该阈值进行滤波时，往往会过平滑掉真实信号。

（2） 极小化风险阈值法

均方差（mean square error，MSE）函数的期望称为风险，极小化风险阈值即MSE意义下的最优阈值。均方差函数描述了滤波结果与真实信号在均方意义下的偏离程度，但是由于在实际应用中，真实信号是未知的，所以人们提出了许多方法来对MSE函数进行估计，这些方法主要有SURE法、交叉验证算法和广义交叉验证算法等。

① SURE法

Donoho等提出通过极小化SURE准则函数来确定阈值，大数定理保证了SURE函数是均方差损失函数的无偏估计。所以，SURE阈值也是最优阈值的无偏估计。

② 交叉验证（cross-validation，CV）算法

G.P.Nason提出的CV算法是一种基于均方差准则确定最优阈值的统计方法。CV算法可以根据实际问题的需要自适应地选取不同的估计准则，以取得更好的滤波效果。当CV算法用于相关噪声滤波时，即使采用基于尺度自适应的阈值，其滤波效果仍不理想。

CV算法采用函数估计的准则为MSE最小，即

1M(t)?

N?(g?ti?gi) (4-12)

i?1N?t表示阈值，gti表示经阈值处理后的函数，g表示真实函数。然而函数g是未知的，所以必须构造出M的一个估计，应用损失函数代替MSE。传统的CV算法是从估计结构中依次去掉每一点，将预测值与去掉的值相比较。

③ 广义交叉验证（generalized cross-validation，GCV）算法

GCV算法是以SURE为基础的，其性能远优于CV算法。GCV是有偏的，其偏差为?2，但其得到的最优阈值是无偏的，且不需要估计噪声方差，故GCV算法应用很广泛。

需要指出的是，GCV算法不能与硬阈值法相结合使用。 （3） 多假设验证法

- 13 -

西安石油大学本科毕业设计（论文）

阈值处理过程可看作是一个多假设验证过程。Abramovich等引入了错判率（false discovery rate，FDR）的概念，它指的是经过阈值化处理后，错判的非零系数个数与所保留的非零系数总数比值的期望。在满足给定的FDR上界的前提下，使所保留的系数个数达到极大值的阈值，即为所要求取的最优阈值。FDR滤波方法是在给出错判率的前提下，导出阈值。此方法与通用阈值法相比有一定灵活性，而且它可以从本质上解释通用阈值依赖于信号长度这个现象。当然FDR滤波方法也有其局限性，其一是如何给定错判率，其二是错判率只考虑了将噪声 误判为信号的情况，却没有考虑将信号误判为噪声的错判率，也就是说着重考虑了滤波结果的平滑性问题。

若阈值在整个滤波过程中固定不变，即对所有的小波系数均采用相等的阈值进行收缩，可将其称为全局阈值，而自适应阈值是根据小波系数的局部特征来对其进行阈值化处理，最简单的自适应阈值是可以是适应于不同尺度上信号和噪声具体特征的阈值，也就是基于尺度自适应（scale-adaptivity）的阈值。

- 14 -

西安石油大学本科毕业设计（论文）

5

小波域阈值滤波的Matlab实现

5.1 一维离散小波变换的Matlab函数

（1） wavedec函数

【函数功能】多尺度一维小波分解 【语法格式】

[C,L]= wavedec(X,N,’wname’)

X代表输入信号；N代表小波分解层数；wname代表使用的小波名称；C为小波分解向量；L是相应地记录向量。

（2） dwt函数

【函数功能】单尺度一维离散小波变换 【语法格式】

① [cA,cD]=dwt(X,’wname’)

② [cA,cD]=dwt(X,’wname’,’mode’,MODE)

X代表输入信号；wname代表使用的小波名称；cA代表计算的低频系数向量；cD代表计算的高频系数向量。

②式表示使用指定的MODE扩展模式计算小波分解。

DWT扩展模式

‘mode’ ‘sym’或’symh’ ‘symw’ ‘asym’ 或’asymh’ ‘saymw’ ‘zpd’ ‘spd’或’spl’ ‘sp0’ ‘ppd’ （3） idwt函数

【函数功能】单尺度一维小波逆变换 【语法格式】

X=idwt(cA,cD,’wname’)

X代表重构信号；wname代表使用的小波名称；cA代表计算的低频系数向量；

- 15 -

DWT扩展模式 对称延拓模式：对称复制边缘值。默认模式 对称延拓模式：对称复制边缘值 反对称延拓模式：反对称复制边缘值 反对称延拓模式：反对称复制边缘值 补零 1阶平滑模式（边缘处进行1阶导数插补） 0阶平滑模式（边缘处进行常数插补） 周期延拓模式（边缘处进行周期延拓） 西安石油大学本科毕业设计（论文）

cD代表计算的高频系数向量。

（4） waverec函数

【函数功能】多尺度一维小波重构 【语法格式】

X=waverec(C,L,’wname’)

X代表重构信号；wname代表使用的小波名称；C为小波分解向量；L是相应地记录向量。

（5） wrcoef函数

【函数功能】由一维小波系数进行单支重构 【语法格式】

①X=wrcoef(‘type’, C,L,’’wname’,N) ②X=wrcoef(‘type’, C,L,’wname’)

X代表重构信号；wname代表使用的小波名称；C为小波分解向量；L是相应地记录向量。N代表在N层计算重构系数向量。

（6） wmaxlev函数

【函数功能】计算小波分解的最大尺度 【语法格式】

L=wmaxlev(S, ‘wname’)

S表示一维信号或者二维图像；wname代表使用的小波名称；L代表函数返回值，是信号或图像的最大分解尺度。使用它可以避免分解时超过这个值。一般情况下，小波分解尺度都大于理论上的最大值。一般一维分解尺度为5，二维分解尺度为3。

（7） detcoef函数

【函数功能】提取一维小波高频系数 【语法格式】 ①D=detcoef(C,L,N) ②D=detcoef(C,L)

D代表提取的一维小波高频系数；C为小波分解向量；L是相应地记录向量。N代表提取N层高频系数。①式表示提取第N层的高频系数，②式表示提取最后一层（即最大层）的高频系数。

（8） appcoef函数

【函数功能】提取一维小波低频系数 【语法格式】

① A=appcoef(C,L, ‘wname’，N) ② A=appcoef(C,L, ‘wname’)

A代表提取的一维小波低频系数；C为小波分解向量；L是相应地记录向量。N

- 16 -

西安石油大学本科毕业设计（论文）

代表提取N层低频系数。①式表示提取第N层的低频系数，②式表示提取最后一层（即最大层）的低频系数。

（9） dwtmode函数

【函数功能】离散小波变换扩展模式 【语法格式】 ① ST=dwtmode ② dwtmode(‘mode’)

dwtmode命令用来设置信号或者图像进行离散小波和小波包变换的扩展模式。扩展模式表示信号或者图像分析时边界问题的处理方法。

dwtmode或者dwtmode(‘status’)显示当前模式。 ST返回当前模式。 ‘mode’为DWT扩展模式。 （10）upwlev函数

【函数功能】一维单尺度小波重构 【语法格式】

[NC,NL,cA]= upwlev(C,L,’wname’)

上格式用于对小波分解结构[C,L]进行单尺度重构，返回新的结构[NC,NL]，并提取最后一层的低频系数向量cA。

C为小波分解向量；L是相应地记录向量；wname代表使用的小波名称。 （11）upcoef函数

【函数功能】一维小波系数直接重构 【语法格式】

① Y=upcoef(O,X, ‘wname’,N) ② Y=upcoef(O,X, ‘wname’,N,L) ③ Y=upcoef(O,X, ‘wname’)

格式①计算向量X向上N步的重构系数。wname代表使用的小波名称，N必须为严格的正整数。如果O=‘a’,则重构低频系数；如果O=‘d’,则重构高频系数。

格式②计算向量X向上N步的重构系数，并提取结果中长度为L的中间部分。 Y=upcoef(O,X, ‘wname’)等同于Y=upcoef(O,X, ‘wname’，1)。

5.2 二维离散小波变换的Matlab函数

（1） upwlev2函数

【函数功能】二维小波分解的单尺度重构 【语法格式】

[NC,NS,cA]= upwlev2(C,S, ‘wname’)

- 17 -

西安石油大学本科毕业设计（论文）

上格式用于对小波分解结构[C,S]进行单尺度重构，返回新的结构[NC,NS]，并提取最后一层的低频系数向量cA。wname代表使用的小波名称。

（2） wenergy2函数

【函数功能】计算二维小波分解能量 【语法格式】

①[Ea,Eh,Ev,Ed]=wenergy2(C,S) ②[Ea,Edetail]=wenergy2(C,S)

格式①返回的Ea是低频部分能量的百分比；Eh,Ev,Ed分别是高频部分水平、垂直、对角方向能量百分比的向量。

格式②返回Ea和Edetail，后者是向量Eh,Ev,Ed的和。 （3） wavedec2函数

【函数功能】二维小波多尺度分解 【语法格式】

[C,S]= wavedec2(X,N, ‘wname’)

使用小波‘wname’返回矩阵X尺度为N时的小波分解。输出是分解向量C和相应地记录矩阵S，N必须是严格的正整数。

（4） dwt2函数

【函数功能】二维小波单尺度分解 【语法格式】

[cA,cH,cV,cD]= dwt2(X,’wname’)

根据输入矩阵X进行小波分解，计算低频系数矩阵cA和高频系数矩阵cH（水平方向）,cV（垂直方向）,cD（对角方向）。

（5） appcoef2函数

【函数功能】提取二维小波分解的低频系数 【语法格式】

①A= appcoef2(C,S, ‘wname’,N) ②A= appcoef2(C,S, ‘wname’)

格式①使用小波分解结构[C,S]，计算N层的低频系数。wname代表使用的小波名称，N必须为严格的正整数。

格式②提取最后一层的低频系数。 （6） wrcoef2函数

【函数功能】由二维小波系数重构单支 【语法格式】

①X=wrcoef2(‘type’, C,S,’wname’,N) ②X=wrcoef2(‘type’, C,S,’wname’)

- 18 -

西安石油大学本科毕业设计（论文）

格式①基于小波分解结构[C,S]，计算尺度为N时的重构系数矩阵。wname代表使用的小波名称，N必须为严格的正整数。如果’type’=’a’，则重构低频系数；’type’=’h’、’v’、’d’时，分别重构高频的水平、垂直、对角系数。

格式②是重构最后一层的系数。 （7） upcoef2函数

【函数功能】由二维小波系数的直接重构 【语法格式】

①Y=upcoef2(O,X, ‘ wname’,N) ②Y=upcoef2(O,X, ‘wname’,N,S) ③Y=upcoef2(O,X, ‘wname’)

格式①计算向量X向上N步的重构系数。wname代表使用的小波名称，N必须为严格的正整数。如果O=’a’,则重构低频系数；如果O=’h’,则重构水平方向高频系数；如果O=’v’,则重构垂直方向高频系数；如果O=’d’,则重构对角方向高频系数；

格式②计算向量X向上N步的重构系数，并提取结果中长度为S的中间部分。 Y=upcoef2(O,X, ‘wname’)等同于Y=upcoef2(O,X, ‘wname’，1)。 （8） waverec2函数

【函数功能】多尺度二维小波重构 【语法格式】

X=waverec2(C,S,’wname’)

X代表重构图像；wname代表使用的小波名称；C为小波分解向量；S是相应地记录向量。

（9） idwt2函数

【函数功能】二维小波单尺度逆变换 【语法格式】

①X=idwt2(cA,cH,cV,cD, ‘wname’) ②X=idwt2(cA,cH,cV,cD, ‘wname’,S)

③X=idwt2(cA,cH,cV,cD, ‘wname’,’mode’,MODE)

X代表重构图像；wname代表使用的小波名称；cA表示低频系数矩阵；cH表示水平方向高频系数矩阵,cV表示垂直方向高频系数矩阵,cD表示对角方向高频系数矩阵。S代表尺度，’mode’表示延拓模式。

（10） detcoef2函数

【函数功能】提取二维小波高频系数 【语法格式】

D= detcoef2(O,C,S,N)

该函数由小波分解结构[C,S]提取尺度为N时的高频系数，如果O=’h’,则提取水

- 19 -

西安石油大学本科毕业设计（论文）

平方向高频系数；如果O=’v’,则提取垂直方向高频系数；如果O=’d’,则提取对角方向高频系数。

5.3 Matlab信号去噪

Matlab小波工具箱中用于信号去噪的一维小波函数是wden.m的wdencmp.m。 进行去噪处理一般有下述3中方法：

（1） 默认阈值消噪处理。该方法利用函数ddencmp生成信号的默认阈值，然后利用函数wdencmp进行消噪处理。

（2） 给定阈值消噪处理。在实际的消噪处理过程中，阈值往往可通过经验公式获得，且这种比默认阈值的可信度高。在进行阈值量化处理时可利用函数wthresh。

（3） 强制消噪处理。该方法是将小波分解结构中的高频系数全部置为0，即滤掉所有高频部分，然后对信号进行小波重构。这种方法比较简单，且消噪后的信号比较平滑，但是容易丢失信号中的有用成分。

默认阈值消噪处理和给定阈值消噪处理在实际中应用的更为广泛。 在Matlab小波工具箱中，设置软或硬阈值的函数为wthresh.m。

该函数根据参数sorh的值计算分解系数的软阈值或硬阈值。其中，硬阈值对应于最简单的处理方法，而软阈值具有很好的数学特性，并且所得到的理论结果是可用的。

阈值的选取有以下4个规则，其中每一条规则对应于函数thselect中输入参数tptr的一个选项。

（1）选项tptr=’rigrsure’,是一种基于Stein无偏似然估计原理的自适应阈值选择。给定一个阈值t，得到它的似然估计，再将非似然最小化，就可以得到所选的阈值。这是一种软阈值估计器。

（2）选项tptr=’sqtwolog’,是一种固定的阈值形式，它所产生的阈值为sprt(2*log(length(X)))。

（3）选项tptr=’heursure’,是前两种阈值的综合，所选择的是最优预测变量阈值。如果信噪比很小，而SURE估计有很大的噪声，此时就需要采用这种固定的阈值形式。

（4）选项tptr=’minimaxi’,也是一种固定的阈值选择形式，它所产生的是一个最小均方差的极值，而不是无误差。

5.4 Matlab图像去噪

在Matlab小波工具箱中，图像的小波去噪与信号去噪类似，步骤也大致相同，所用函数也一致，而所不同的是，处理时用二维小波分析代替了一维小波分析。

二维小波分析用于图像去噪的步骤如下：

- 20 -

西安石油大学本科毕业设计（论文）

（1）二维图像信号的小波分解。在这一步，应当选择合适的小波和恰当的分解层次（记为N），然后对待分析的二维图像信号X进行N层分解计算。

（2）对分解后的高频系数进行阈值量化。对于分解的每一层，选择一个恰当的阈值，并对该层高频系数进行软阈值量化处理。在此，阈值选择规则同前面的信号处理部分。

（3）二维小波的重构图像信号。同样的，根据小波分解后的第N层近似（低频系数）和经过阈值量化处理后的各层细节（高频系数），来计算二维信号的小波重构。

- 21 -

西安石油大学本科毕业设计（论文）

6 小波多分辨率分析对地震记录去噪与地质信息提取的

Matlab实现

6.1 程序编写

本论文中所用到的程序全是基于Matlab的小波处理工具箱中的函数。其算法如下：

读入地震记录 对记录进行二维小波多尺度分解 提取各层高频系数 由thselect函数生成各层高频系数的阈值，进行阈值去噪 二维小波单尺度重构 利用wextend函数进行信号延拓 输出处理后的地震记录

在对地震记录进行小波多尺度分解时，小波的道数和采样点成倍减少，例如，尺度为1时的重构记录为100×1001，尺度为2时的重构记录变为50×502，尺度为3时的重构记录为25×252。如果不对重构后的信号进行处理，将无法从程序中输出，即得不到处理结果。在本论文中，使用信号延拓方法，解决了这个问题。

6.2 模型去噪与地质信息提取

地震记录模型共有60道，每道1200个采样点，采样率为0.5ms。

从图6-1-1可以看出，模型中有两个水平层1和2，背景噪声很强，只有上部的水平反射层1比较明显，下部的反射层2淹没在噪声中，几乎无法分辨，

利用Matlab，使用db1小波，对原始模型记录进行3尺度分解，按上述流程进行处理后，再分别对尺度为2、尺度为3的单尺度重构记录进行成图。

- 22 -

西安石油大学本科毕业设计（论文）

图6-1-2和图6-1-3均为使用附录[1]中的程序对原始模型记录处理后的结果。 下图为一地震记录模型：

CDP Time/ms

图6-1-1 原始模型记录 CDP Time/ms

图6-1-2 消噪后尺度为2时重构记录

从图中可以看出，噪声得到了有效压制，背景较清晰，上部水平地层1更加明显，

- 23 -

西安石油大学本科毕业设计（论文）

下部水平地层2比上图要清晰些，基本上从噪声中突出出来了。 CDP Time/ms

图6-1-3 消噪后尺度为3时重构记录

与图6-1-2相比，噪声更少，背景清晰，水平层1和2均非常明显，突出于背景之上。滤波效果较好。

经处理可以看出，模型由噪声和两个水平地层组成。尺度为3时的处理效果最好。

6.3 实际资料去噪与地质信息提取

利用模型对论文中的处理方法进行验证之后，仍需观察该方法对实际资料的处理效果。

应用Matlab对实际资料的处理方法与流程基本与对模型的处理时一致的，所用的程序见附录[2]。

图6-2-1是一张实际地震资料的图形。

实际地震资料共有100道，每道采样点为1001个，采样间隔为1ms。 由图6-2-1可看出，在区域1处可能有断层，但在图中不是很明显。而且，只有区域2和3处的地层反射较明显，可作为标准反射层组，其余部分都被噪声淹没，不能很好的反映出地层形态。

为了较好的提取地质信息，需要对实际资料先进行多尺度分解，再对高频成分进行去噪处理。

- 24 -

西安石油大学本科毕业设计（论文）

CDP Time/ms

图6-2-1 实际资料原始图

在小波域滤波方法中，为了得出最佳处理结果，分别对一下参数进行了对比： （1）由于小波是任意选的，而且分解尺度数也可任意选择，只要不要超过最大分解尺度数（对于一维信号，最大尺度数为5；对于二维图像，最大尺度数为3），所以在本论文中，使用各种小波对资料进行了处理，以便于对比，得出适用于该资料的最佳小波。

（2）由于最大分解尺度为3，尺度为4时图像失真严重，故本论文中所用处理均采用3尺度分解，并分别对各尺度重构图像进行对比。

（3）由于阈值函数分为硬阈值函数和软阈值函数，为了找出哪种阈值函数更适合该资料，对资料分别用两种函数进行了处理。

（4）阈值的选择原则有4种，分别为’rigrsure’’sqtwolog’ ’heursure’ ’minimaxi’，本论文也应用它们分别对资料进行了去噪处理。

对比所用的程序大致相同，只要根据对比需要改变程序中的参数即可。但由于使用的参数不同，由信号所分解得到的低频和高频成分矩阵可能会发生变化，这是需要改变程序中信号延拓的行列数，这需根据实际情况而定。

下面为经过处理后的图像。 当使用小波db1处理时： （1）使用硬阈值函数：

①阈值的选择原则为：heursure，处理结果为：

- 25 -

西安石油大学本科毕业设计（论文）

CDP Time/ms

图6-2-2 去噪处理后尺度为1时的图像

CDP Time/ms

图6-2-3 去噪处理后尺度为2时的图像

- 26 -

西安石油大学本科毕业设计（论文）

CDP Time/ms

图6-2-4 去噪处理后尺度为3时的图像

②阈值的选择原则为：minimaxi，处理结果为：

CDP Time/ms

图6-2-5 去噪处理后尺度为1时的图像

- 27 -

西安石油大学本科毕业设计（论文）

CDP Time/ms

图6-2-6 去噪处理后尺度为2时的图像

CDP Time/ms

图6-2-7 去噪处理后尺度为3时的图像

- 28 -

西安石油大学本科毕业设计（论文）

③阈值的选择原则为：rigrsure，处理结果为：

CDP Time/ms

图6-2-8 去噪处理后尺度为1时的图像

CDP Time/ms

图6-2-9 去噪处理后尺度为2时的图像

- 29 -

西安石油大学本科毕业设计（论文）

CDP Time/ms

图6-2-10 去噪处理后尺度为3时的图像

④阈值的选择原则为：sqtwolog，处理结果为：

CDP Time/ms

图6-2-11 去噪处理后尺度为1时的图像

- 30 -

西安石油大学本科毕业设计（论文）

CDP Time/ms

图6-2-12 去噪处理后尺度为2时的图像

CDP Time/ms

图6-2-13 去噪处理后尺度为3时的图像

- 31 -

西安石油大学本科毕业设计（论文）

（2） 使用软阈值函数

①阈值的选择原则为：heursure，处理结果为：

CDP Time/ms

图6-2-14 去噪处理后尺度为1时的图像

CDP

Time/ms

图6-2-15 去噪处理后尺度为2时的图像

- 32 -

西安石油大学本科毕业设计（论文）

CDP Time/ms

图6-2-16 去噪处理后尺度为3时的图像

②阈值的选择原则为：minimaxi，处理结果为：

CDP Time/ms

图6-2-17 去噪处理后尺度为1时的图像

- 33 -

西安石油大学本科毕业设计（论文）

CDP Time/ms

图6-2-18 去噪处理后尺度为2时的图像

CDP Time/ms

图6-2-19 去噪处理后尺度为3时的图像

- 34 -

西安石油大学本科毕业设计（论文）

③阈值的选择原则为：rigrsure，处理结果为：

CDP Time/ms

图6-2-20 去噪处理后尺度为1时的图像

CDP Time/ms

图6-2-21 去噪处理后尺度为2时的图像

- 35 -

西安石油大学本科毕业设计（论文）

CDP Time/ms

图6-2-22 去噪处理后尺度为3时的图像

④阈值的选择原则为：sqtwolog，处理结果为

CDP Time/ms

图6-2-23 去噪处理后尺度为1时的图像

- 36 -

西安石油大学本科毕业设计（论文）

CDP Time/ms

图6-2-24 去噪处理后尺度为2时的图像

CDP Time/ms

图6-2-25 去噪处理后尺度为3时的图像

- 37 -

本文来源：https://www.bwwdw.com/article/dhhg.html