小波分析

一、小波的高分辨率层序分析技术 在碎屑岩沉积地层中，自然伽马曲线与其他测井曲线相比，能最敏感地反映泥质含量变化，利用高分辨率的自然伽马曲线进行高分辨率层序地层划分是最为有效的。鄂尔多斯盆地构造单一，缺少大的构造与断层，没有断层和构造变形引起的地层重复，是利用测井曲线进行小波分析不同级别基准面旋回的理想场所。在此以鄂尔多斯盆地城壕地区城86井自然伽马测井数据进行小波变换分析，求取测井曲线中保存的不同级别的基准面旋回。 以城86井位于鄂尔多斯盆地三叠系延长组华池三角洲沉积体系为例，该地区位于三角洲沉积体系的前缘，该钻穿三角洲沉积体系前缘的钻井测井类型多样、层位齐全(图12-1)。在鄂尔多斯盆地三叠纪延长湖的演化中，自上而下经历了长1至长10共十个油层组发育期，湖泊经历了从发生、发展到消亡的完整过程。长10油层组是湖泊演化的初始阶段，长9和长8油层组以湖泊扩张为主，长7油层组是湖泊面积最大的时期，而长6油层组是湖盆由兴盛走向衰落的转折期，其底部为深灰泥岩，顶部以块状砂岩为主，地层厚度90～120m。长4+长5油层组是在湖退大趋势背景下的次一级湖进-湖退沉积旋回，主要为浅灰色、灰色细砂岩和深灰色、灰色泥岩互层的组合。长3、长2和长1油层组为湖盆大规模收缩至消亡阶段的沉积物，发育大规模的细砂岩、粉砂岩。在湖退的大背景下，长2和长1油层组之间仍有次一级的湖进，在长2油层组顶和长1油层组底有较多的暗色泥岩。通过对城86井长9至长1油层组自然伽马(GR)曲线的小波分析后，可以提供湖泊从湖侵到鼎盛再到衰退、最后消亡的演化细节和高分辨率的地层结构及层序界面特征。 利用一维小波变换对城86井自然伽马曲线进行五次分解，分别获得了五条小波变换曲线，其中一条(D1)代表低频信号旋回、一条(D2)代表了中频信号旋回、三条(G1、G2和G3)代表高频信号旋回曲线(图12-1)。 图12-1城壕地区城86井延长组自然伽马曲线小波分析与层序划分图 小波变换得到的高频信号，充分表现了小波分析“数学放大镜”的功能，由它可以判断相变剧烈程度和相转换面。在高频基准面旋回G3中，在相变剧烈处，曲线表现为剧烈跳跃;图12-1中的1190m、1421m、1675m和1893m处，G3曲线变化剧烈，峰值较大。参照岩屑录井和其他资料，中频信号曲线D2对应于中期基准面旋回，与偏心率周期中气候波动引起的基准面升降和物质供给变化有关，反映同一沉积体系中相似或相邻相序的进积-退积过程，是一套水深变化幅度不大的、彼此间成因联系密切的地层叠加所组成的湖进-湖退沉积序列。高频旋回曲线G1对应于短期基准面旋回，与斜率周期中气候波动引起的基准面升降和A/S值变化有关，反映韵律性沉积旋回的进积-退积过程，是一套具有低幅水深变化的、彼此间成因联系极为密切，或由相似岩性、岩相地层叠加组成的湖进-湖退沉积序列。高频旋回曲线G2超短期基准面旋回，与岁差周期中气候波动引起的基准面升降和A/S值的变化有关，反映单一岩性或数个岩性组合的进积、退积、加积或进积-退积过程，是一套代表最小成因地层单元的单一岩性或相关岩性的叠加变化的特点。对比高频旋回曲线G3上的异常相变点、低频信号曲线D1周期变换界面以及长期基准面旋回的层序界面之间的对应，三者之间存在着良好的对应关系，由此可以看出由小波变换后的不同频率曲线划分的层序与原始GR曲线具有很好的对应关系。 在层序划分中，由于测井曲线是不同级次基准面旋回的叠加，各级次基准面旋回的相互影响和噪声的污染，使基准面旋回的划分难以把握，特别是长周期的基准面变化从原始曲线上难于直观的观察，而通过小波变换分析所得到的低频信号曲线能够直观地反映曲线的变化趋势。其次在原始测井曲线中，层序界面位置难以确定，特别是在测井过程中上下围岩的影响、为消除噪音而进行的滑动平均处理、中长周期对曲线的控制作用都使原始曲线不能准确地反映岩石的界面，而小波分析处理后的高频信号可以准确地反映信号的奇异相变点。高频信号G3曲线中，可以清晰地确定相转换面和剥蚀面的位置。由此可见，用小波分析来进行高分辨率层序地层划分的方法是可行的，也是可信的。 二、延长组基准面旋回与米兰科维奇旋回的对应关系 米兰科维奇旋回是原南斯拉夫学者米兰科维奇在解释第四纪冰期成因时提出的气候与地球轨道定量关系的理论。地球的偏心率、地轴倾角及岁差发生周期变化直接影响气候(Huggett，1991)，气候的周期性波动影响湖平面的变化，被沉积物记录，并保存在地质地层中。米兰科维奇旋回岁差长短周期分别约为23ka、19ka;地轴倾角长短周期分别约为54ka、41ka;偏心率周期为100ka、123ka、413ka，而且在过去的500Ma内，偏心率旋回没有变化，地轴倾角和岁差旋回有随着时间的推移而降低的趋势(Berger，1989)。 鉴别地层中保存的沉积旋回是否属于米兰科维奇旋回，必须具备三个条件: (1)该旋回的年限及其在地层中的重复性必须符合米兰科维奇旋回的周期性; (2)该旋回必须与米兰科维奇提出的气候效应有内在联系，气候效应表现在沉积记录的物理、化学及生物的特征; (3)由于区域性乃至全球性气候制约的结果，该旋回必须同时发育于不同的沉积地点及不同的沉积相区(Kauffman，1988;郝维城等，2000)。 (一)延长组中米兰科维奇旋回发育条件 鄂尔多斯盆地晚三叠世延长组为大型坳陷湖泊发育期，接受陆相河流—三角洲—湖泊沉积，保存了厚达千米的生、储、盖组合(图12-2)。根据岩、电性特征及含油性的差异，自上而下可划分为长1至长10油层组共十个油层组。其中长10至长7油层组为湖盆扩张期，长6至长1油层组为湖盆收缩期。郑荣才等(2002a)研究志丹地区长6油层组的高分辨率层序地层学的超短期、短期和中期基准面旋回的沉积记录时，分析认为这三种类型的沉积旋回在地层中的发育具有规律性的重复性，沉积过程中的韵律性旋回等特征符合米兰科维奇旋回，并具备了米兰科维奇旋回的前两个条件。因此要验证米兰科维奇旋回在鄂尔多斯盆地延长组发育中的一致性，还要研究不同的沉积地点及不同沉积相区是否均有保存。为此以陇东地区庆36井延长组为例(图12-2)，在基准面旋回分析和高分辨率层序地层划分的基础上，按沉积过程的韵律性旋回，结合数学方法对延长组可能的米兰科维奇旋回进行分析，同时与陕北地区(张小会等，2002)、靖安油田(李斌等，2005)延长组米兰科维奇旋回进行对比，从而确定鄂尔多斯盆地晚三叠世延长组中发育的米兰科维奇旋回的一致性。 (二)延长组米兰科维奇旋回分析 频谱分析技术是研究周期性现象中最为常用的一种统计分析方法，其中最常用的是快速傅里叶变换法(FFT)。本文采用该方法，对鄂尔多斯盆地陇东地区延长组自然伽马曲线进行频谱分析，从中提取其所包含的古气候周期变化信息。古气候的周期变化反映在沉积韵律的变化中，这种变化保存在地层中，最终通过钻井资料记录在测井曲线上。测井曲线不但分辨率高，达到米级尺度，而且能够敏感、连续地反映所测地层的特征，沉积物的结构、构造、岩性及岩相等周期性变化必定记录在测井数据中，是高频旋回研究的主要资料。通过数学计算，可以把这些沉积韵律周期定量分析出来。 图12-2鄂尔多斯盆地城壕地区庆36井延长组高分辨率层序地层综合划分图 1.自然伽马曲线的频谱分析 利用第四章所提供的数学方法，首先对区内庆36井延长组长1至长10油层组自然伽马曲线按深度等间距0.2m取值，进行数字离散化，然后对所得离散化数据进行快速傅里叶变换(FFT)的时-频转换，将资料的时间(深度)领域转换到频率领域，得出频谱曲线(图12-3)。图中横坐标代表频率，纵坐标代表相对功率，相对功率表示该频率的沉积旋回在地层中的重要性，功率值越大，表明该周期的沉积旋回在地层中出现的越频繁，因此高点处的频率对应于曲线的主要频率，这样就可以找出频谱曲线中的主要频率值，进而可以求出相应的波长，得出旋回周期。进行快速傅里叶变换，求取地层中保存的主频周期旋回，分析各油层组旋回厚度及其比值结果(表12-1)。从表12-1看出，陇东地区庆36井延长组沉积旋回基本存在8.69～17.56m的长周期，3.22～6.39m的较长周期，2.57～5.04m的中短周期，2.03～4.48m的较短周期，以及1.63～3.02m和1.25～2.63m的短周期。各沉积旋回周期的比值近似为1∶(0.359～0.380)∶(0.279～0.30)∶(0.228～0.267)∶(0.169～0.187)∶(0.142～0.152)。 图12-3庆36井长4+长5油层组自然伽马曲线FFT频谱分析曲线 2.米兰科维奇旋回的确定 根据Berger(1992)计算的地质历史时期米兰科维奇旋回周期的变化，求得晚三叠世时长、短两个地轴倾角周期分别为44.3ka和35.4ka;长、短两个岁差周期分别为21.3ka和17.8ka;按偏心率周期为恒定不变的123ka，偏心率周期与长、短两个地轴倾角周期的比率分别为1∶0.36和1∶0.288，偏心率周期与长、短两个岁差周期的比率分别为1∶0.173和1∶0.145。从表12-1中可以看出，本次计算的测井结果中，B旋回与A旋回的比值变化范围0.359～0.380，与对应的米兰科维奇旋回偏心率周期和地轴倾角长周期之比0.36的误差不超过5.6%;C旋回与A旋回的比值变化范围0.279～0.30，与对应的米兰科维奇旋回偏心率周期和地轴倾角短周期之比0.289的误差不超过3.8%;E旋回与A旋回的比值变化范围0.169～0.187，与对应的米兰科维奇旋回偏心率周期和岁差长周期之比0.173的误差不超过8.1%;F旋回与A旋回的比值变化范围0.142～0.152，与对应的米兰科维奇旋回偏心率和岁差短周期之比0.145的误差不超过4.8%。此外，D旋回与A旋回的比值变化范围0.228～0.267，也有很好的稳定性。上述各比值中，除D旋回与A旋回的比值外，其余均与米兰科维奇旋回周期之比误差不超过8.1%。米兰科维奇旋回主要为偏心率、地轴倾角和岁差三种周期、五种旋回类型。 表12-1鄂尔多斯盆地庆36井延长组地层旋回厚度和厚度比 *代表旋回厚度，**代表比值，/代表该井没有显著周期。 鄂尔多斯盆地陇东地区延长组属于稳定克拉通背景发育的大型陆相坳陷型湖泊三角洲沉积体系，该地层保存完整，泥岩不但发育，而且对气候变化反映敏感。在构造运动稳定区，短周期的基准面旋回变化受控于湖平面高频振荡性波动，湖平面的这种变化又起因于米兰科维奇轨道旋回所驱动的古气候变迁(郑荣才等，2002a)。古气候的周期性变化，势必造成短期基准面上升和下降，与之对应的沉积物粗细旋回变化也被记录在地层中。古气候温暖潮湿和干旱周期性变化，在湖相地层中表现得尤为突出，随着气候的改变，大量泥岩的出现与消失均被忠实地记录在地层中，而这些由气候引起地层旋回的变化又被自然伽马测井曲线很好地反映出来，这也使沉积地层旋回定量化研究成为可能。频谱分析的结果也反映了本区稳定的旋回厚度比值与米兰科维奇旋回比值之间很好的一致性，因此可以认为，米兰科维奇旋回是影响本区地层沉积旋回发育的主要因素，古气候周期性变迁决定着坳陷湖盆湖泊三角洲沉积体系中高频沉积旋回的发生、发展和定格。该区延长组的沉积旋回周期分别为123，44.3，35.4，21.3和17.8ka，其中偏心率周期、地轴倾角周期和岁差周期引起的地层旋回厚度变化范围分别为8.69～17.56m、2.57～6.39m和1.25～3.02m。 三、盆地范围内的延长组米兰科维奇旋回发育的一致性分析 气候对沉积旋回的制约是区域性甚至是全球性的，因此米兰科维奇旋回在鄂尔多斯盆地延长组中不同的沉积地点、不同的沉积相区同样发育，可作为区域地层等时对比的重要依据。 比较鄂尔多斯盆地陇东、靖安和陕北三地延长组地层保存的沉积周期旋回，三者均与米兰科维奇周期旋回之间存在着较好的对应关系，而且同样保存着米兰科维奇三个周期旋回(表12-1，表12-2)。分析三地的沉积环境，陇东、靖安和陕北三地在鄂尔多斯盆地所处的位置不同，在延长组发育期间，各地的沉积相发育不同，如长4+长5油层组，当靖安、陕北处于三角洲平原沉积亚相(张小会等，2002;李斌等，2005)时，陇东地区仍保持水体较深的浅湖-前三角洲沉积环境，但是同样保存着米兰科维奇旋回。这表明在同一时期不同的沉积环境中，在不同的沉积地点受米兰科维奇旋回的影响是一致的，这一研究成果符合文献(Kauffman，1988;郝维城等，2000)提出的鉴定米兰科维奇旋回的第三个条件，由于区域性乃至全球性气候制约的结果，该旋回必须同时发育于不同的沉积地点，不同的沉积相区，从而进一步支持了米兰科维奇理论。 表12-2鄂尔多斯盆地陕15井与XP14井延长组地层米兰科维奇旋回厚度表 注:“无”代表该井没有显著周期。

高分辨率层序地层学理论的核心内容是“在基准面旋回变化过程中，由于可容纳空间与沉积物补给通量比值的变化，相同沉积体系域或相域中发生沉积物的体积分配作用和相分异作用所导致的沉积物保存程度、地层堆积样式、相序、相类型及岩石结构和组合类型发生的规律性变化”(Cross，1994)。由于基准面的变化是海平面、构造沉降、沉积物补给、沉积负荷补偿、沉积压实与沉积地形等各种因素变化的综合反映，因此，通过碎屑岩的厚度变化、粒度大小、有机质含量和沉积物类型及岩石的结构构造特征表现出来，而这些地质现象又被高分辨率的测井曲线记录下来，这就为利用数学方法定量分析旋回信息提供了依据。随着计算机技术的发展，数字信号处理、地学信息系统分析使得地质问题的定量化成为地学研究的热点之一。频谱分析和小波分析是进行基准面旋回定量分析的重要技术方法之一，二者均能从复杂的叠加信号提取相对独立的天文周期信息，为基准面旋回的分析提供重要的理论基础和技术支持。 一、小波分析原理 小波是一个衰减的波形，它在有限的区域里存在(不为零)，且其均值为零。小波是尖锐变化而且是不规则的波形，因此用小波能更好地刻画信号的局部特征。小波变换是小波分析的核心，设测井曲线为f(h)，小波变换的算法如下: 0 公式中ωf(i，j)为尺度i下刻度j处的小波系数;a为尺度参数，a=2－j;b=a·j。 采用低通和高通滤波器，求取低频系数CA1和高频系数CD1，然后再分解CA1为CA2和CD2，再分解CA2为CA3和CD3，如此类推(图4-16)，然后再对低频系数和高频系数进行重构高频信号D3、D2、D1和低频信号A3，分别得到不同频率的周期旋回。 图4-16一维小波变换三次分解与重构示意图 二、小波分析在高分辨率层序地层划分中的意义 1.识别转换基准面和剥蚀面 小波变换是傅里叶变换的发展，其实质是引入伸缩、平移思想，对不同频率成分自动地选取时域和取样步长，从而能够聚焦到物体的任意微小细节，因此小波变换被誉为“数学显微镜”。利用小波变换能够有效地检测包含非常重要的信息的奇异点和不规则点的信号。在测井信号中，奇异点常常代表地层信息变化剧烈点，如具有特殊意义的剥蚀面和转换基准面。因此利用小波变换的这一特点，可以很容易地识别测井曲线的剧变点以及剧变程度，分析出对于高分辨率层序划分具有特殊意义的剥蚀面和转换基准面，从而确定层序边界，如图4-17所示为鄂尔多斯盆地A16井延长组与富县组的分界线的位置为不整合界线，经小波变换后检测包含明显信息的奇异点信号，同样在湖泛面的位置也检测出不规则的信号。 图4-17A16井小波分析的延长组与富县组界面和旋回周期性识别 2.识别不完整基准面旋回周期 在地质过程中，存在众多的侵蚀、沉积物过而不沉和欠补偿沉积等作用，导致地质记录和保存的不完整性，测井信号所反映出的周期性、旋回性必然是残缺不全的。而传统的频谱分析方法，如经典的傅里叶变换，是对信号在整个时间域上的周期性进行分析，而不能弥补这些不完整的旋回记录。而小波分析对于地层中发育的单独的上升半旋回、下降半旋回和不对称旋回，可以实现对信号在时间域和频率域局部化分析，又如图4-17所示，在鄂尔多斯盆地A16井1860～2020m段的自然伽马测井曲线中，直观的分析曲线中以正旋回为主，逆旋回难以观察，经小波分析后，可以明显的分析出该井段的测井曲线包含2.5个较高频旋回，由此来看，小波分析为基准面旋回的识别提供了技术支持。 3.识别测井曲线的隐含周期 高分辨率等时地层对比的关键是识别地层记录中这些代表多级次基准面旋回的地层旋回，在基准面旋回的识别中，测井序列是迄今为止所能获得的分辨率最高、连续性最好的地质数据，其中蕴藏着丰富的地质信息。不同的测井数据在不同程度上记录着地质演化的历史，从不同侧面反映着地层形成演化的条件和影响因素，如海平面变化、古环境、古地理、古气候信息及其变化情况等，是了解地质过程的最好工具。但由于地质过程的多时间多尺度特征和各种串级过程，地层沉积序列实际上是各种地质周期沉积响应的叠加，再加上不确定因素和局部因素引起的随机波动的干扰，从测井曲线上很难直观地从测井信号中判读各种隐含周期。 图4-18叠加信号的小波分析分解图 小波分析技术可以把测井信号分解在任意精度的不同频带内(图4-18)，根据感兴趣的信号频带范围，把信号在一定尺度上分解，从而提取相应频带的信息，得到相应的地质周期(余继峰等，2003)，对信号中的低频成分，采用宽的时间窗，得到低的频率分辨率;对信号中的高频成分，采用窄的时间窗，得到高的频率分辨率(Daubechies，1991)。由此可以利用小波分析的这种自适应特征对测井曲线进行多尺度分析，选取信号中代表地质长周期的低频部分来确定大的层序地层旋回，中等频率用来确定中等的地层层序，高频部分则代表短周期的旋回，可以用于小层的精细对比和划分。如图4-18中的A曲线是一个由3个正弦函数f(t)=sin2πt+sin3πt+sin4πt组成的叠加信号，经过小波分析后能够将其分解为3个独立的信号，如图4-18中的曲线B、C和D，由此可见，小波分析可以用于测井曲线中隐含的旋回周期的识别。 三、频谱分析原理 1.频谱分析基本原理 快速傅里叶变换法(FFT)是频谱分析技术研究周期性最为常用的一种统计分析方法，主要通过对一复合波系进行傅里叶变换，将其分解成若干振幅和相位不同的间谐波，并找出其中振幅最大的波，即该复合波中的主要频率。 傅里叶变换函数常以连续函数得出，若变换函数为x(t)，则傅里叶变换由下式: 0 式中:t为时间;f代表频率。 在实际过程中，采样都是离散和有限的，如果是连续信号，在应用计算机处理时也需要进行截断并离散化，因此在实际数据处理时，都采用离散傅里叶变换，其变换公式如下: 假定有一段N项离散时间序列xm，其离散傅里叶变换为 0 其中Xk称为频谱值，k=0，1，2，…，N－1。 快速傅里叶变换(FFT)是计算离散傅里叶变换(DFT)的一种快速算法，能迅速提高DFT的运算速度。FFT的算法种类很多，基于FFT算法中的频率抽选法(Cara，1982)是其中的常用的一种。 当k为偶数时，公式为 0 k=0，1，…，A－1;A=N/2，此处，k已用2k所代替。 当k为奇数时，用2k+1代替k，可以得到: 0 k=0，1，…，A－1;A=N/2。 由于频谱值Xk没有直观意义，习惯上总把其转化成能量，构成直观的能量-频率图，对应于频率fk的谐振能量为 Pk=(real(Xk))2+(imag(Xk))2 其中，k=0，1，2，…，N－1。 2.频谱分析的Matlab实现 对于大多数地质工作者来说，要进行上述复杂的频谱分析数学计算，可能会有一定的难度。由Mathwork公司推出的Matlab软件集数字分析、矩阵运算、信号处理和图像处理、显示于一体，构成了一个方便灵活的、界面友好的用户环境。在这个环境下，对所求解的问题用户只需简单地列出数学表达式，其结果便以数值或图形的方式显示出来。以下是应用Matlab实现复合信号频谱分析进行信号分解的一个实例，以期为Matlab在地质中的应用起到抛砖引玉的作用。 利用Matlab中的快速傅里叶变换函数fft()实现频谱分析的实例如下: 0 对上述信号signal，利用Matlab进行快速傅里叶变换: R=fft(signal，N) 式中:R为测井信号的快速傅里叶变换频谱值序列，fft为快速傅里叶变换函数，N是signal的长度。 频谱值R对应于频率fk的谐振能量为 Pf(k)=Rk*conj(Rk) 式中:Pf(k)为频率fk的谐振能量值，Rk频率fk的快速傅里叶变换(FFT)值，conj(Rk)为Rk的共轭复数。由于所求功率为模数，频谱曲线左右对称，因此仅取其中的一半进行计算，所以k=0，1，2…Round(N/2)，Round(N/2)是对N/2进行取整。 0 运行结果如图4-19所示。 图4-19复合正弦曲线的Matlab频谱分析 上述实例可以看出，使用Matlab进行数据处理十分方便与灵活。Matlab为用户提供了大量的功能函数，可以为研究人员避免大量重复性的数学运算、而把更多的精力集中到专业的方法研究中提供了便利。 对于测井曲线频谱分析，是测井曲线的数字化处理，这些与测井相关的问题正是Matlab很方便解决地。 3.频谱分析在高分辨率层序地层划分中的意义 (1)天文周期的确定。在六个级次的基准面旋回级次划分方案中所强调的主控制因素各不相同，以陆相盆地为例，低频长周期的旋回如超长期和长期旋回主要受构造作用控制，而高频短周期旋回如中期、短期和超短期基准面旋回则分别受天文因素的偏心率长周期、偏心率短周期和岁差周期所控制(郑荣才等，2001)。因此，如何从复杂的地层信息中识别出保存有天文周期的信息是进行高频基准面旋回分析的重要内容，而测井曲线的频谱分析可以从地层中获取包含不同级次的天文周期信息，从而更有效地进行高频周期旋回控制因素的分析。 (2)估算地层沉积速率。通过频谱分析所得到的天文周期旋回，在地质历史时期，组成米兰科维奇旋回的偏心率周期和岁差周期是稳定的，因此在已知沉积周期厚度和周期持续的时间，就可以得到相应的沉积速率。