基于优化Morlet小波的分频扩散滤波断层强化方法与流程

本发明涉及油田开发技术领域，特别是涉及到一种基于优化morlet小波的分频扩散滤波断层强化方法。

背景技术：

在地震数据采集过程中，不可避免地会受到传感器灵敏度、野外噪声干扰和数模转换量化问题等因素的影响，导致地震数据信噪比低等问题，给后续构造解释带来很大困难。针对常规叠后地震数据噪声压制的问题，国内外学者开展了大量研究，提出了中值滤波、高斯滤波、时频域滤波等方法。上述方法在噪声压制方面效果较好，但由于无法识别地震反射的局部方向信息，在噪声压制的同时，会破坏边缘信息。

为此，国内外学者将各向异性扩散滤波方法引入到地震数据滤波处理中。perona等最先提出各项异性扩散方程(即p-m模型)，fehmers等将该方法引入到地震资料的处理和解释中，weickert和lavialle等对p-m模型进行了改进，提出一致性增强扩散滤波(ced滤波)，通过使用扩散张量作为扩散系数，实现了真正意义上的非线性各项异性扩散滤波。此后，众多学者在ced滤波的基础上进行改进，并引入地震相干属性、非平稳相似系数等断层因子，提出诸如双向增强扩散滤波、构造导向扩散滤波等。上述扩散滤波方法的基本思想都是通过地震剖面局部图像的方向特征，判断是否存在断层等异常地质体边界，进而确定是否进行滤波，从而达到去噪的同时保护边界的目的。

需要注意的是，滤波处理在提高信噪比的同时会降低地震数据的分辨率。为了解决上述问题，通过对分频地震数据进行各向异性扩散滤波处理，然后重构得到全频地震数据，可以一定程度提高滤波后地震数据的分辨率。目前，用于分频处理的主要是谱分解法、傅里叶变换、s变换和小波变换等。其中小波变换可以随着频率成分的改变而自动调整窗口形态，能够对信号进行局部化分析，被称为数学显微镜，广泛用于提高地震资料的分辨率。

公开号为cn106054250a的专利申请变频分量扩散滤波融合的地震资料噪声消减方法中提到类似的实现过程：(1)输入三维叠后地震数据体；(2)将三维叠后地震数据体分解为多个变频分量数据体；(3)对变频分量数据体进行地震反射倾角导向扩散滤波处理，得到滤波后的变频分量地震数据体；(4)将滤波后的变频分量地震数据体进行加权重构与融合，获得压制了噪声，同时保留了有效的不连续异常新信息的三维地震数据体；(5)对步骤(4)得到的三维地震数据体求取地震体曲率属性；(6)输入目的层位数据，从上述三维数据体中抽取岩层切片，用于储层解释。

基于分频扩散滤波的地震数据去噪方法的关键是获取高分辨率的分频地震数据，而利用小波变换进行分频处理的关键选取合适的小波母函数，其选取标准是小波母函数的波形和频谱与地震子波的波形与频谱相似，并满足容许条件。目前，morlet小波具有良好的局部性能，且当morlet小波主频大于1hz时近似满足容许条件，被广泛用于地震数据的小波分频中，但传统的morlet小波及前人改进后的morlet子波均存在其波形和频谱无法同时匹配地震子波的波形和频谱的问题，从而导致分频结果分辨率较低。为此我们发明了一种新的基于优化morlet小波的分频扩散滤波重构断层强化处理方法，解决了以上技术问题。

技术实现要素：

本发明的目的是提供一种通过优化morlet小波提高小波分频结果的分辨率及准确度，然后结合各项异性扩散滤波，实现地震数据去噪保边的同时提高其分辨率，从而提高地震数据低序级断层识别能力的基于优化morlet小波的分频扩散滤波断层强化方法。

本发明的目的可通过如下技术措施来实现：基于优化morlet小波的分频扩散滤波断层强化方法，该基于优化morlet小波的分频扩散滤波重构断层强化处理方法包括：步骤1，将叠后地震数据转换为matlab可识别的数据格式；步骤2，对地震数据进行频谱分析并提取子波；步骤3，参照地震子波的波形和频谱特征，改进morlet子波母函数，构建morlet子波函数组合；步骤4，利用优化morlet子波进行小波分频处理；步骤5，计算叠后地震数据的结构张量，进而构建基于结构张量迹的新断层因子；步骤6，在结构张量迹的保护下，对分频地震数据进行各项异性扩散滤波；步骤7，对分频扩散滤波结果进行小波反变换，并重构得到提高了信噪比和分辨率，并保护了断层信息的地震数据。

本发明的目的还可通过如下技术措施来实现：

在步骤2中，利用快速傅里叶变换对地震数据进行频谱分析，并进行地震子波提取。

在步骤3中，改进的morlet子波母函数公式为：

式中：fc为小波母函数的主频；t为时间。

在改进morlet小波母函数的基础上，构建不同尺度的morlet小波函数组合：

式中：fc为小波母函数的主频，a为尺度参数，该参数与分频数据主频fa的关系为：a＝fs*fc/fa，其中fs为地震数据采样频率。

在步骤4中，利用优化morlet小波开展叠后地震数据的小波分频处理，其小波变换公式为：

式中：wψ(a，b)为小波分解后的小波系数；ψ(t)为优化morlet小波母函数，a为尺度因子，b为时移因子。

在步骤5中，结构张量的计算公式如下所示：

式中：fδ＝gδ*f，gδ为尺度为δ的高斯核函数，*为卷积运算符，f为地震数据；为梯度算子；t为矩阵转置符号；算子sδ是在尺度δ下的结构张量，包含着与梯度一样的方向信息；

在此基础上，对上述sδ进一步利用尺度为ρ的高斯核函数gρ进行滤波，得到另一尺度下的结构张量：

sρ＝gρ*sδ

对于平稳变化的区域，其不同尺度滤波后的结构张量反映的方向信息应该一致，而对于非平稳变化的区域，受高斯核函数尺度参数的影响，其滤波后反映的方向信息可能发生变化；因此，基于结构张量特征值的变化构建断层因子：

式中：sρ和sδ是经过两个不同尺度滤波的结构张量；tr表示矩阵的迹，其值为矩阵特征值的和；由ε的计算公式可知，ε取值范围为0<ε<1，ε取值较大时，反映图像平稳变化区域，而ε取值较小，则反映图像边界区域。

在步骤6中，在地层连续区域，进行扩散滤波，在断层附近，减弱或停止扩散滤波，从而保护断层信息；

针对不同分频结果，利用步骤3构建的结构张量构建扩散张量；对sρ进行特征值分解得出特征值μ1≥μ2和对应的特征向量(v1，v2)；扩散张量的特征向量取结构张量的特征向量，使得扩散张量能够反映图像梯度变化，扩散张量的特征值则由研究人员赋值及其结构张量的特征值计算得到：

式中：α用来控制沿图像梯度变化较大方向的扩散强度，为了保护图像的结构信息，α通常设为接近0的正数，c为阈值，通常取值为1，m取值为1；

在引入结构张量迹后，扩散张量d构建如下：

式中：a、b、c是扩散张量为矩阵表示形式时的参数。

最终，基于结构张量迹保护的各向异性扩散滤波计算公式如下：

式中：div为散度计算符号，▽为梯度算子，f为地震数据，分别为地震数据在x方向、对角方向和y方向上的二阶偏导。

在步骤7中，小波反变换公式为：

式中：ψ(t)为优化morlet小波母函数，t为时间，ω为圆频率，wψ(a，b)为小波分解后的小波系数；a为尺度因子，b为时移因子。

该基于优化morlet小波的分频扩散滤波断层强化方法还包括，在步骤7之后，在完成断层强化处理后，将处理后数据输出为标准sgy格式的数据。

本发明中的基于优化morlet小波的分频扩散滤波重构断层强化处理方法，在前人基础上进一步改进morlet小波，使优化morlet小波的波形和频谱与地震子波的波形和频谱同时达到最佳匹配，提高分频结果的分辨率和准确度。在此基础上，结合各向异性扩散滤波，实现低序级断层强化处理。与现有的地震数据分频去噪手段相比，本发明主要有三点优势：

(1)基于优化morlet小波的小波分频相比于matlab自带小波分频，在地震数据带宽范围内大幅度提高了分频结果的分辨率和准确度；

(2)基于优化morlet小波的分频扩散滤波断层强化方法相比于现有软件提供的扩散滤波方法，在保证去噪保边效果的同时提高了滤波处理后地震数据的分辨率；

(3)基于地震图像结构张量的计算构建基于结构张量迹的新断层因子，相比于地震相干等断层属性，该断层因子在程序运行中计算得到，避免了地震数据的多次输入输出及属性提取，提高了断层强化处理的效率。

附图说明

图1为本发明的基于优化morlet小波的分频扩散滤波断层强化方法的一具体实施例的流程图；

图2为本发明的一具体实施例中叠后地震数据的地震子波波形图和频谱图；

图3为本发明的一具体实施例中不同morlet小波母函数的波形及频谱图；

图4为本发明的一具体实施例中优化morlet小波母函数波形和频谱与地震数据地震子波和频谱对比图；

图5为本发明的一具体实施例中matlab自带morlet小波与本发明的优化morlet小波对单道地震信号分频结果的对比图；

图6为本发明的一具体实施例中已有软件的扩散滤波结果与本发明的扩散滤波结果对比图。

具体实施方式

为使本发明的上述和其他目的、特征和优点能更明显易懂，下文特举出较佳实施例，并配合附图所示，作详细说明如下。

如图1所示，图1为本发明的基于优化morlet小波的分频扩散滤波重构断层强化处理方法的流程图。

步骤1：将叠后地震数据转换为matlab可识别的数据格式；

步骤2：利用快速傅里叶变换对地震数据进行频谱分析，利用geoframe软件进行地震子波提取；图2展示了实际地震数据提取的地震子波及频谱。

步骤3：参考地震数据的频谱特征及地震子波的波形特征，改进morlet小波母函数，并构建优化morlet小波函数组合；

传统的morlet小波母函数公式为：

其中：fc为小波母函数的主频。传统morlet小波母函数当主频fc较小时(matlab自带morlet小波的主频为0.8126hz)，其波形旁瓣较少，与地震子波较为相似，而当fc与地震数据主频一致时，旁瓣很多，无法满足小波分频的要求。

前人改进的morlet小波母函数公式为：

式中：fc为小波母函数的主频，c为调节参数。通过改变调节参数，该morlet小波母函数的波形与地震子波波形相似度很高，且其频谱特征可以匹配地震数据的频谱，但此时该小波母函数不满足容许条件：

参考前人思路，本文改进的morlet小波母函数公式为：

式中：fc为小波母函数的主频。将前人调节参数c替换为主频fc，同时将公式中2替换为优化后的morlet小波母函数其波形和频谱特征与地震子波波形和地震数据频谱特征匹配度均很高，且满足容许条件。

由图3可以看出，相对于传统的morlet小波及前人改进的morlet小波，本文改进的morlet小波其波形和频谱特征更合理，且满足小波反变换的容许条件。图4则进一步证明本文改进的morlet小波波形及频谱与实际地震数据的子波波形及频谱匹配度更高。

在改进morlet小波母函数的基础上，构建不同尺度的morlet小波函数组合：

式中：fc为小波母函数的主频，a为尺度参数，该参数与分频数据主频fa的关系为：a＝fs*fc/fa，其中fs为地震数据采样频率。该小波函数组合可以用于后续小波分频处理中，得到不同频率的分频结果。

步骤4：利用优化morlet小波开展叠后地震数据的小波分频处理，其小波变换公式为：

式中：wψ(a，b)为小波分解后的小波系数。ψ(t)为优化morlet小波母函数，a为尺度因子，b为时移因子。

由图5可以看出，不同主频条件下，基于优化morlet小波的分频结果精度明显高于matlab自带的morlet小波，且对比分频结果与地震数据，分频后地震数据的分辨率得到明显提高。

步骤5：构建叠后地震数据局部结构张量，并对结构张量开展不同尺度的高斯滤波，构建基于结构张量迹的断层因子；

结构张量是目前获取图像局部方向信息较为稳定的方法，其计算公式如下所示：

式中：fδ＝gδ*f(gδ为尺度为δ的高斯核函数，*为卷积运算符，f为地震数据)；为梯度算子；t为矩阵转置符号。算子sδ是在尺度δ下的结构张量，包含着与梯度一样的方向信息。

在此基础上，对上述sδ进一步利用尺度为ρ的高斯核函数gρ进行滤波，得到另一尺度下的结构张量：

sρ＝gρ*sδ

对于平稳变化的区域，其不同尺度滤波后的结构张量反映的方向信息应该一致，而对于非平稳变化的区域，受高斯核函数尺度参数的影响，其滤波后反映的方向信息可能发生变化。因此，基于结构张量特征值的变化构建断层因子：

式中：sρ和sδ是经过两个不同尺度滤波的结构张量；tr表示矩阵的迹，其值为矩阵特征值的和。由ε的计算公式可知，ε取值范围为0<ε<1，ε取值较大时，反映图像平稳变化区域，而ε取值较小，则反映图像边界区域。

步骤6：在结构张量迹的保护下，开展分频结果的各向异性扩散滤波。在地层连续区域，进行扩散滤波，在断层附近，减弱或停止扩散滤波，从而保护断层信息。

针对不同分频结果，利用步骤3构建的结构张量构建扩散张量。对sρ进行特征值分解得出特征值(μ1≥μ2)和对应的特征向量(v1，v2)。扩散张量的特征向量取结构张量的特征向量，使得扩散张量能够反映图像梯度变化，扩散张量的特征值则由研究人员赋值及其结构张量的特征值计算得到：

式中：α用来控制沿图像梯度变化较大方向的扩散强度，为了保护图像的结构信息，α通常设为接近0的正数，c为阈值，通常取值为1，m取值为1。

在引入结构张量迹后，扩散张量d构建如下：

式中：a、b、c是扩散张量为矩阵表示形式时的参数。

最终，基于结构张量迹保护的各向异性扩散滤波计算公式如下：

式中：div为散度计算符号，▽为梯度算子，f为地震数据，分别为地震数据在x方向、对角方向和y方向上的二阶偏导。

步骤7：对上述扩散滤波后的分频结果进行小波反变换，并重构各分频数据的反变换结果得到断层强化后的地震数据。其中小波反变换公式为：

式中：ψ(t)为优化morlet小波母函数，t为时间，ω为圆频率，wψ(a，b)为小波分解后的小波系数；a为尺度因子，b为时移因子。

步骤8：在完成断层强化处理后，将处理后数据输出为标准sgy格式的数据。图6显示，基于本文提出的断层强化处理方法得到滤波后地震剖面，其信噪比和分辨率明显高于现有软件的处理结果，证明该方法可以有效提高复杂断块区地震数据的品质，有助于低序级断层的精细解释。

在本发明的一实施例中，通过构建结构张量迹作为断层因子，来约束扩散滤波，通过提取地震相干属性、体曲率属性等断层属性作为断层因子，也可实现同样的目的。

本发明的基于优化morlet小波的分频扩散滤波断层强化方法通过对传统morlet小波的改进，提出了优化morlet小波母函数；将基于优化morlet小波的小波分频方法与各向异性扩散滤波方法相结合，在去噪保边的同时提高了地震数据的分辨率。