小波分析在层序地层划分中的应用

1.小波分析简介

20世纪80年代后期至今，一种著名的、在各行各业有重要应用价值的数学理论和方法技术在科学技术界得到了广泛的重视和采用，它就是被誉为“数学显微镜”的小波分析(李世雄，1994)。小波分析的主要功能和特点是，它具有多分辨分析或多尺度分析功能，可以把信号分解成各种不同的尺度成分它具有很强的局部分析功能，同时具有时间(或空间)域和频率域的局部分析性质，它可自动地通过伸缩、平移聚焦到信号的任一细节对其加以分析(侯遵泽，1998)

(1)小波分析基本原理。小波(wavelet)，即小区域的波，是一种特殊的长度有限、平均值为0的波形。它有两个特点:一是“小”，即在时域都具有紧支集或近似紧支集二是正负交替的波动性。如果用小波和构成傅里叶分析基础的正弦波做比较的话，傅里叶分析所用的正弦波在时间上没有限制，从负无穷到正无穷，但小波则倾向于不规则与不对称。

傅里叶分析是把信号分解到一组相互正交的正弦波上的，也就是基函数，我们可以把基函数看成是度量信号某些特征的一把“尺子”，傅里叶分析度量的就是信号的频谱特征，但是如果这把“尺子”过于规则，有时候就不能十分精确地表达信号蕴含的信息，而在小波分析中，“尺子”换成了规则程度更低的小波函数，从而可以更加有效地表达信号中信息的成分。

小波变换对不同频率在时域上的取样步长是调节性的，即在低频时小波变换的时间分辨率较差，而频率分辨率较高在高频时小波变换的时间分辨率较高，而频率分辨率较低(图2－13)，这正符合低频信号变化缓慢而高频信号变化迅速的特点(胡昌华，1999)。这就构成了利用小波变换进行信号分析的基础。

图2－13 数字信号的小波变换

(2)一维连续小波变换。小波变换实际上是求取信号在各小波函数上的投影值。每个小波函数均由一个母小波函数经过尺度伸缩与时间平移得来的。信号分析的一般思路就是分解与组合，寻找一组最能代表信号特征的函数形式，将信号用这些量来逼近，或者写成这些量的线性组合的形式。

小波分析的思想来源于伸缩和平移方法:对波形的尺度伸缩就是在时间轴上对信号进行压缩与伸展，而时间平移就是指小波函数在时间轴上的波形平行移动。

(3)离散小波变换。由于连续小波变换的伸缩和平移系数是相互独立的，所以通过伸缩和平移得到的各个小波函数之间有一定的相似性，但由于这两个系数之间的独立，就引入了信息的冗余。在分辨率一定的情况下，描述了多余的信息，使得反映信号特征的一些参数相互重叠，给我们的分析带来不便，但这些特点可以用在本身就有自相似性的信号上，可以让我们更清楚地看到信号自身的自相似性。

此外，由于冗余信息的存在，也使得小波逆变换的重构过程不唯一，也就是说，由同一母小波生成的不同的小波变换函数可能重构成同一个信号。为了减少冗余信息，就引入了离散小波变换的概念，其中的伸缩和平移系数是可数的，重构过程用求和的形式给出。如果伸缩和平移系数满足一定的对应关系，则称为二进小波变换(尺度以2的幂的形式给出)。离散小波变换主要是建立在二进制小波变换的基础上的。

测井曲线数据也恰好是离散数据，符合离散变换的要求。在利用小波分析进行层序地层划分时，主要是对测井曲线进行多尺度分解，得到不同尺度下的小波变换图，利用其表现出来的特征来划分不同级次的层序。

2.利用小波分析进行层序地层划分

利用小波分析方法是层序地层划分方法上的一种新的尝试，其目的是尽量减少层序划分过程中的主观因素，依靠地层自身表现出来的客观特征来识别层序、准层序组以及准层序。在我们研究的沉积岩地层中，沉积物的特征可以反映沉积时水体的特征。随着沉积水深的变化，沉积物中多种特征都会相应的发生变化，如放射性物质含量、有机质含量或其他微量元素的含量等，这种变化就会在相应的测井曲线上反映出来。而沉积水深变化受到了多种因素的影响，有长期和短期的旋回，是多个不同周期旋回的叠加，因此测井曲线应该是沉积地层中某种随水深变化的特征的多种频率变化的响应的叠加。也就是说，测井曲线中包含着沉积水深不同周期变化的信息，是多个沉积水深变化周期相互叠加的响应。而小波分析能够将信号分解为不同频率不同周期的旋回，因此，小波分析的特点恰好可以和测井曲线的特点相对应，利用小波分析的方法可以比较准确地将测井曲线中相互叠加的反映水深变化的不同周期的信息分别识别出来，识别出的这些信息就可以用来进行沉积旋回的划分。

同时，小波分析方法还可以帮助解决传统研究方法所不能解决的一些难题，如大段单一岩性地层中的沉积旋回识别。大段单一岩性尤其是大段泥岩、页岩，并不是一个小的沉积旋回里沉积的产物，相反，应是一个相当长时期沉积下来的，但是通过传统的岩性划分方法却很难将其划分开，这就给准层序甚至准层序组的划分造成了困难。小波分析方法可以较好地解决这一问题，利用这种方法可以从测井曲线的细微变化中识别沉积间断和沉积旋回。

(1)测井曲线的选择。不同的曲线具有不同的地质含义，进行相同的变换可能会得到不同的结果。但在研究中通过对GR、AC、COND、电阻率等多条曲线进行小波变换后对比发现，不同测井曲线所得出的变换结果尽管形态上不完全一样，但在旋回的划分上却比较一致(图2－14)。图中曲线a是COND测井曲线经过db5小波变换后的结果，曲线b是同一井段AC曲线变换后的结果。出现这个结果是由于虽然不同的曲线代表着不同的地层响应，会呈现出不同的特征，但地层中各种参数的变化主要受沉积环境的影响，会随着沉积环境的旋回变化呈现出基本一致的旋回特征。这也从一个方面反映了小波变换在沉积旋回划分中的客观性。因此，可以选择目标井的测量精度较高、质量较好的曲线来进行小波变换，进而进行沉积旋回的划分。

图2－14 对COND和AC曲线进行小波变换结果对比

(2)小波的选择。同傅里叶分析不同，小波分析的基(小波函数)不是唯一存在的，所有满足小波条件的函数都可以作为小波函数，那么小波函数的选取就成了十分重要的问题，实际选取小波的标准主要有以下三种。

1)自相似性原则:对二进小波变换(因为在正交小波变换中，取样的方式就是按照小波函数取样的，所以不存在这个问题)，如果选择的小波对信号有一定相似性，则变换后的能量就比较集中，可以有效减少计算量。

2)判别函数:针对某类问题，找出一些关键性的技术指标，得到一个判别函数，将各种小波函数代入其中，得到一个最优原则。

3)支集长度:大部分应用选支集长度在5～9之间的小波。因为支集太长会产生边界问题，支集太短不利于信号能量的集中。

但在实际应用中，因为大部分信号的信息量太大，很难找到相应的模式，因此主要是依靠经验来选取。根据前人研究经验及作者对不同函数所做结果的对比，本书采用的是Daubechies小波，阶数为5。

Daubechies小波是由著名小波学者Ingrid Daubechies所创造，她发明的紧支集正交小波是小波领域的里程碑，使得小波的研究由理论转到可行。这一系列的小波简写成dbN，其中N表示阶数。

(3)工作流程。测井曲线能比较准确地反映井旁地层的电性、物性等特征，但往往会受到测井仪器、钻井泥浆等其他非地层因素的干扰，且不同频率的旋回相互叠加，对正确识别和划分沉积旋回造成一定的影响。小波分析能真正消除干扰信号，放大真实信息，按不同频率反映出测井曲线中包含的真正旋回特征，以此划分不同级别层序单元，同时还可以在划分高精度沉积旋回的基础上，与Fischer图解相结合划分出体系域。

MATLAB软件的小波分析工具箱是一种比较常用的工具。MATLAB是Math works公司于1982年推出的一套高性能的数值计算和可视化软件。MATLAB的推出得到了广大专家学者的广泛关注，其强大的扩展功能为各个领域的应用提供了基础。各个领域的专家学者相继推出了MATLAB工具箱，包括信号处理、神经网络、图像处理、小波分析等。其中小波分析工具箱可以满足对测井曲线进行小波变换的需要。

图2－15 小波分析流程图

在对测井曲线进行小波变换时，首先需要对所研究层段的顶底界面进行准确的标定，然后将需要变换的该深度段的测井曲线数值建立单独的文本文件作为原始文件。将原始文件导入后保存成.m格式的信号文件。选择MATLAB软件的wavelet(小波分析)工具箱进行离散一维小波变换，小波类型选择db，阶数为5，最大级数定为12，将上述参数选好后进行分析，即可得到一组12条不同级次的db5小波变换曲线(图2－15)。此外对其进行连续一维小波变换，可以得到小波的频谱分析图，选择合适的最大显示值，根据频谱图上图形的闭合方向可以区分出层序的界面和层序单元(图2－6，图2－7)。

(4)单井分析实例。牛100井位于牛庄洼陷西部，地层以砂泥岩互层为主，岩性变化较快(图2－16)。利用小波分析方法对AC、R25两条测井曲线进行了一维连续小波变换，分别得到其小波变换谱系图，对AC曲线进行了一维离散变换，得到不同阶数的小波，根据与地震、测井及录井岩性资料的对比，选用d11，d9，d7三个层的小波变换曲线分别进行层序、准层序组和准层序的划分。

将传统划分方法所得的结果与小波分析方法所得结果进行对比可以比较明显的看出，在层序和准层序组的划分上两种方法划分的层序单元基本一致，可以相互验证。在准层序级别上的划分，小波分析方法的优势就比较明显地体现了出来，划分的旋回数较多，精度也有提高。这也正是小波分析作为“数学显微镜”的特点所决定的。

从图2－16中小波分析得到的d11曲线可以看出，这一段地层可以划分为两个大的旋回，对应两个层序，谱系图上的特征也比较明显。其中每个大的旋回又可以分为三个次一级的旋回，在d9及谱系图上可以找到相关界面，相当于每个层序划分出三个准层序组，每个准层序组在测井曲线及录井资料上也有较明显的反旋回特征。在进行准层组的划分时，小波分析方法可以划分出肉眼不易识别的旋回，从而提高了划分精度。整段地层一共可以划分为21个准层序，代表不同的沉积旋回。其旋回特征在d7曲线上有较好体现，从谱系图上也可以找到各界面的标志。从测井曲线和岩性上看，基本上每一个准层序都是一个反旋回，代表着一期的水体变换，这也完全符合层序地层学的基本原理。

图2－16 牛100井小波分析资料的层序地层划分

王62井位于牛庄洼陷东部，与牛100井相比，划分出的各层序单元的厚度发生了明显的变化，但数目基本一致，这也证明了小波分析划分层序地层的结果是比较可靠的。通过对AC曲线的小波变换得到AC曲线的小波变换谱系图和小波变换曲线，如图2－17所示。从谱系图和d11曲线上可以将整段地层划分为两个大的旋回，分别对应层序Ⅲ和层序Ⅳ。其中每个层序又可以划分为3个准层序组，在d9曲线上可以看到相应的旋回出现，谱系图上可以找到界面的标志(图2－17)。王62井这一段地层一共可以划分成20个准层序，缺失第一个准层序。各准层序在岩石类型、颜色和测井曲线上基本上可以看出反旋回特征，符合层序地层划分方法。

通过牛100井、王62井的划分可以看出，小波分析方法在砂泥岩互层的地层中有较好的应用效果，可以提高层序划分的精度和准确性。在层序划分中有比较好的可重复性，使得全区的划分结果比较客观和统一，减少了人为因素造成的干扰。

首先将你的数据txt文档放在matlab的搜索路径中，确保matlab能找得到，使用load函数可以将改txt文档中的数据读到一个数组中，比如load data.txt就可以将数据读入一个名为data的数组中，当然，data也是一个两列的数组（时间和幅值），使用t = data(:, 1) 和m = data(:, 2)可以提取出data数组中的第一列和第二列的数据构成向量t和m，在函数中就可以使用这两个向量了。

按照你上面的代码，a3 dd3 dd2 dd1是应该都为行向量，它们的长度a3和dd3长度相等，dd2约等于其一倍，dd1又是dd2的约一倍左右，c1=[a3 dd3 dd2 dd1]是4个行向量按顺序依次组成总小波系数数组，这代码没啥问题，waverec和wavedec函数都是依靠CL组构的，它对C处理没啥难以理解的，L数组存放的就是a3 dd3 dd2 dd1长度值，靠他在C中取出最高阶的a和其他各阶的d，a和d的小波系数本来就不等长，这样才能实现DWT理论中各层分解中小波函数伸长的效果。就是因为它们不等长，用元胞数组 *** 作太麻烦，所以才把系数混排成一行，再加一个长度数组L共同 *** 作来访问数据。所以不知道你对这有啥疑问？

---