用MATLAB实现地震数据的小波变换

1、文章编号: 10011749 (2002) 02016306用m a tl a b实现地震数据的小波变换李晋平 , 王大庆 , 许新刚(中国矿业大学 地球信息技术系, 江苏徐州221008)摘要: 小波变换借助于时频局部分析特性, 已经成为地震资料处理中的一种重要工具。作者在阐述其基本原理后, 讨论了如何选择小波基 ( 函数) , 并利用m a tla b 语言及小波工具箱, 实 现了 se g2y 格式地震数据文件的小波分析, 取得了令人满意的效果。关键词: 小波变换; 时频局部性; 小波基; 地震数据中图分类号:p 631. 4+ 2文献标识码:aim pl em net ing w a

2、vel et tra nsfo rmo f se ism ic d a ta us ing m a tl a bl i j in 2p ing , w a n g d a2q ing , xu x in 2gang(d ep a rtm en t of th e e a rth i nf orm a tion t ech nolog y , c h ina u n iv ersity of m in ing and t ech nolog y , x uz h ouj iang su221008, c h ina)a bstra ct: w avelet t ran sfo rm , dep

3、end ing upon th e p rop er ty of t im e2f requency loca liza t ion , h a s a lready beenan im po r tan t too l in seism ic da ta p rocessing. t h e au tho r s exp la in th e ba sic th eo ry and d iscu ss th e select ion of w avelet ba sis (funct ion ) , t h en w e m ak e u se of th e m a tla b langu

4、age and it s w avelet too lbox to im p lem en tw avelet ana ly sis of seism ic da ta f ile w ith se g2y fo rm a t and th e resu lt s show th e eff iciency of th e m ethod.key words: w avelet t ran sfo rm ; t im e2f requency loca liza t ion; w avelet ba sis; seism ic da ta引言0小波概念首先是由法国从事石油信号处理的地球物理学家

5、 j. m o r let 在 1984 年提出来的。后来, 数学家 y. m eyer 构造了具有一定衰减性质的光滑小波函数 , 而 i. d aubench ies 构造了具有有限支集的正交小 波 基。 1989 年, s. m a lla t 与 y. m eyer 合作建立了构造小波基的通用方法2多尺度分析 m ra (m u t i2r eso lu t ion a na ly sis) , 并提出了著名的m a lla t 快速算法, 使小波变换成为重要的实用工具。 传统的傅里叶变换是一种纯频域的分析方法, 时频分析比较差。小波变换不仅继承和发展了窗口傅里叶变换的局部化思想, 而且

6、克服了窗口大小不随频率变化、缺乏离散正交基的缺点, 这在理论和实际应用 中都是非常重要的性质。 借助小波变换, 可以进行时频局部分析, 小波变换在地震勘探领域具有很大的应 用潜力。m a tla b 是一种优秀的数值计算编程语言, 主要功能有: 数值分析、矩阵运算、信号处理和图形显 示。m a tla b 还为不同领域的专业人员提供了扩展的工具箱, 从而方便了学习和使用。 作者在本文中通收稿日期: 2001- 10- 19过使用m a tla b 中的小波分析工具箱, 编制了m a tla b 程序, 实现了 se g2y 地震数据格式文件的小波分析。 输入为原始剖面文件, 输出为同样格式的

7、se g2y 小波各尺度分解剖面文件。 从整个实现过程来 看, 我们可以充分体会m a tla b 在解决实际问题上的方便性和高效性。1 小波变换及其分类小波函数的确切定义为: 设 ( t) 为一平方可积函数 ( 即 ( t) l 2 (r ) ) , 若其傅里叶变换 ( ) 满足条件1 ( ) 1 d 2(1)r则称 ( t) 为一个基本小波或小波母函数, 并称式 (1) 为小波函数的允许条件。函数 f表达式为:( t) 的连续小波变换的1 f( t) ( t-b) d tw t f (a, b) = =(2)aa r1t- b 其中 a, b ( t) =) 称为依赖于参数 a、b 的小波

8、基函数。(aa1. 1连续小波变换连续小波变换是指尺度参数 a 和平移参数 b 连续变化。 连续小波在计算机实现时, 必须加以离散化,这一离散化都是针对连续的尺度参数 a 和平移参数 b 进行的。由于有许多好性能的连续小波变换, 诸如平 移不变性、伸缩共变性以及冗余性等, 都非常适合于做信号分析, 所以在连续小波变换中, 仅要求小波函数满足式 (1) 的允许条件即可, 这样我们在处理实际问题时选择小波就有很大的自由度。 连续小波变换的冗 余性在去噪、进行数据恢复及特征提取时, 均可以获得好的效果, 但这需要以牺牲计算量、存储量为代价。1. 2离散小波变换把连续小波变换中的尺度参数 a 和平移参

9、数 b 的离散化形式分别取作 a= aj ,0 b= k aj b0。 在实际0中, 常取 a0 = 2, b0 = 1, 即取二进制小波。为了能够使用著名的m a lla t 快速算法, 离散小波变换多采用正交小波变换的形式, 即要求所选小波是正交小波。选用正交小波进行离散变换没有连续小波的变换系数的冗 余性, 在保证不丢失原始信号信息的情况下, 可以大大地减少计算量。小波基(函数) 的选取2与标准傅里叶变换相比, 小波分析中所用的小波函数不具有唯一性, 即小波函数 ( t) 具有多样性 ( 我们可以选择非正交小波、正交小波、双正交小波、甚至是线性相关的)。在小波分析的工程应用中, 一个很重

10、 要的问题是最优小波基的选取问题, 这是因为应用不同的小波基解决同一个问题会产生不同的结果。 另外, 根据要从信号中提取的信息不同, 也应恰当地选择小波和构造小波函数。 作者在本文中针对地震信号处理的要求, 首先从理论上讨论了小波函数的选取准则, 并利用小波分析后重构信号和原始信号的误差大小来比较判定小波基函数的好坏, 最终选定解决实际问题的最优小波基。在对地震信号进行小波分析时, 我们可以考虑选择和地震子波形状很相近的波形作为小波 (满足允许条件)。 常见的地震子波有零相位子波, 最小相位子波和混合相位子波等, 文献1给出了这些常用地震子 波的模拟公式。 此外, 我们会注意到, 存在的各种小

11、波中, 有很多核小波函数和地震子波的形状很相近, 也可以用来做试验。连续小波变换是一种冗余变换, 子波在空间两点之间的关联增加了分析解释变换结果的困难, 而离散 的正交小波变换则不会出现这种缺陷。选择和构造一个正交小波要求其具有一定的紧支集、平滑性和对称 性。紧支集保证有优良的空间局部性质; 对称性保证子波的滤波特性有线性相移, 不会造成信号的失真; 平 滑性保证频率分辨率的高低。但是上述三点不可能同时得到满足。紧支撑性与平滑性二者不可兼得, 要求 小波具有较高的光滑性, 必然要求增加小波支集的长度; 反之, 为了保证小波分析的局部特性, 利于算法实2 期李晋平等: 用m a tla b 实现

12、地震数据的小波变换165现, 支集的长度要尽量小, 但这又保证不了光滑性。综合考虑, 我们必须采取某种折衷作法, 保证一定的紧支撑性、对称性和平滑性来选择正交小波。我们选用 db4 小波抽取某一标准道地震数据进行离散小波变换, 变换结果如图 1 所示。图 1 误差分析图f ig. 1 e r ro r ana ly sis由图 1 可以看出, 原始信号和重建信号完全一样, 只是有微小的误差, 误差数量级为 10- 7。 这说明我 们选用的小波来分解和重构信号是可行的。下面我们考虑两种情况来对误差结果进行分析:(1) 滤波器长度不同 (选择同一家族小波, 如对于 db 族某一小波 dbn , 滤

13、波器长度为 2n )。(2) 滤波器长度相同 (不同家族小波, 选择 db 族小波和 sym 族小波)。 误差的标准有两种: 最大误差 (一道地震记录中所有误差值的绝对值的最大值) ;平均误差 (一道地震记录中所有误差绝对值的和的平均)。表 1 滤波器长度不同的误差结果比较t ab. 1 com p a r ison of d ifferen t f ilter2leng th er ro r resu lt s表 2 滤波器长度相同的误差结果比较t ab. 2 com p a r ison of sam e f ilter2leng th er ro r resu lt s ( 注: db

14、是指数学家 i. d aubench ie s 提出的一组小波的简称)误差结果分析:对于支集长度较小的 db1、db2、db3 小波来说, 由于小波的光滑度不够, 重构信号的误差都比较大。而对于 db7db9 小波来说, 随着滤波器长度的增加, 紧支集区间也变大, 小波的光滑性得到了保证。但误差逐渐增大说明了由于支集长度的增加, 局部性下降。从 db4 和 db6 小波中可以看出, 在所有 db 家族小波db4sm y4db8sym 8最大误差1. 157910- 75. 096410- 8最大误差2. 691810- 71. 369710- 8平均误差2. 688010- 81. 28801

15、0- 8平均误差7. 034310- 73. 383710- 9db1db2db3db4db5db6db7db8db9最大误差1. 455210- 117. 394910- 86. 503710- 71. 157910- 71. 538510- 78. 011910- 81. 338110- 72. 691810- 62. 896510- 6平均误差3. 117210- 121. 368510- 81. 472010- 72. 668010- 73. 860710- 72. 087010- 83. 372610- 86. 509210- 87. 034310- 7中, 误差较小的有两个, 说明

16、这两种小波能够很好地顾及正交小波的紧支集和平滑性。对于滤波器长度相同 ( 支集长度相同) 的不同小波类型来说, 我们看出在选定的两种滤波器长度下, sym 族小波都优于 db 族小波。 从两种类型的小波函数图形可以看出, sym 族小波比 db 族小波与地震子 波更为相近, 理论分析和误差结果说明了这一点。语言实现小波变换3m a tla b使用上面的误差分析方法, 我们选择比较好的、适合地震数据处理的小波基函数, 运用m a tla b 语言编制小波程序来实现小波分析。3. 1小波变换程序的实现funct ion fo r- cw t (f ilein , f ileou t, t race

17、s, sam p les, sca le, w avelet)连续小波变换%fo r- cw t (f ilein, f ileou t, t races, sam p les, sca le, w avelet)输入f ilein:原始的地震剖面文件 (3 . sgy)f ileou t: 小波分解剖面文件 (3 . sgy)整个剖面文件的总道数t races:sam p les: 每一道地震记录的采样点数小波变换尺度值sca le:w av lelet: 小波类型: m o r l, m o r let 小波; gau s, gau ssian 小波; m exh , 墨西哥帽子小波生成由

18、输入所指定的小波分解剖面文件输出inf ile= ch a r (f ilein)o u tf ile= ch a r (f ileou t) t yp e= ch a r (w avelet)fp 1= fop en ( inf ile, r) ;% 打开一原始的输入文件fp 2= fop en (o u tf ile, w ) ; % 要写入的小波分解剖面文件cc= f read (fp 1, 1800, in teger3 2) ; % 读取 sgy 文件的文件头信息fw r ite (fp 2, cc, in teger32) ;fo r k = 1: t races % 循环处理bb

19、= f read (fp 1, 120, in teger3 2) ; % 读取 sgy 文件的道头信息fw r ite (fp 2, bb, in teger3 2) ;s= f read (fp 1, sam p les, in teger32) ; % 读出一道数据w = cw t ( s, sca le, w avelet) ; % 小波变换后该道的变换系数数据fw r ite (fp 2, w , in teger3 2) ; % 将该道的小波变换数据写入到文件中endfclo se (fp 1) ;fclo se (fp 2) ;funct ion fo r- dw t (f il

20、ein , t races, sam p les, w avelet)% funct ion fo r- dw t (f ilein, t races, sam p les, w avelet)离散小波变换输入原始的地震剖面文件 (3 . sgy)整个剖面文件的总道数 每一道地震记录的采样点数%f ilein:t races:sam p les:w avelet: 小波类型; 可供选择的有d aubech ies 小波、co if let s 小波、sym let s 小波等生成不同分层结构的小波分解剖面文件本程序以 3 层分解为例。w av01. sgy 表示尺度 1 的小波分解剖面文件 (

21、3 . sgy)输出inf ile= ch a r (f ilein)t yp e= ch a r (w avelet)2 期李晋平等:用m a tla b 实现地震数据的小波变换167fp = fop en ( inf ile, r) ;fp 1= fop en (w av01.fp 2= fop en (w av02. fp 3= fop en (w av03.sgy,sgy, sgy,w ) ;w ) ;w ) ;cc= f read (fp ,in teger32) ;% 读取 sgy 文件的文件头信息写入1800,fw r ite (fp 1,fw r ite (fp 2, fw r

22、 ite (fp 3,cc, in teger3cc, in teger3cc, in teger32) ;2) ;2) ;fo r k = 1: t racesbb= f read (fp ,120, in teger3 2) ;% 读取 sgy 文件的道头信息fw r ite (fp 1,fw r ite (fp 2, fw r ite (fp 3,bb, in teger3bb, in teger3bb, in teger32) ;2) ;2) ;s= f read (fp , sam p les, in teger32) ;% 读出一道地震数据% 3 层分解% 对应该道的尺度 1 低频

23、系数重构c , l = w avedec ( s, 3, w avelet) ; a1= w rcoef (a, c , l , w avelet, 1) ; fw r ite (fp 1, a1, in teger3 2) ;a2= w rcoef (a, c , l , w avelet, 2) ;fw r ite (fp 2, a2, in teger3 2) ;a3= w rcoef (a, c , l , w avelet, 3) ;% 对应该道的尺度 2 低频系数重构% 对应该道的尺度 3 低频系数重构fw r ite (fp 3, a3, in teger32) ;endfclo

24、 se (fp ) ; fclo se (fp 1) ; fclo se (fp 2) ; fclo se (fp 3) ;% 小波变换例子t races= 280; sam p les= 700; sca le= 1. 5 fo r- cw t ( 329 stk. sgy, fo r- dw t ( 329 stk. sgy,329cw t. sgy,t races, sam p les, sca le, m o r l)% 连续小波变换t races, sam p les, db4)% 离散小波变换% 329 stk. sgy是一原始地震剖面文件;m o r l和db4分别为所选的小波。

25、图形显示:由于输入和输出数据文件格式为 se g2y 格式文件, 所以这里我们采用方便的地震资料处理软件v ista66 来完成数据的显示。图 2 山东许厂地区某一原始地震剖面f ig. 2 o r ig ina l seism ic p rof ile of shandong xuchang a rea3. 2计算实例图 2 为山东许厂地区某一地震剖面的局部, 目的层 3 煤在 410m s 附近, 在第 90 100 道间发现有断点 存在。 其小波分解剖面 (最佳尺度) 如图 3 所示。图 3 山东许厂地区某一小波分解剖面f ig. 3 w avelet decom po sing p rof ile of shandong xuchang a rea可以