基于变分模态分解的三维沉积旋回模型构建方法

技术领域

本发明属于石油天然气勘探开发领域，尤其涉及基于变分模态分解的三维沉积旋回模型构建方法。

背景技术

目前，在地震勘探领域，层序地层学在沉积岩的划分、对比和分析中发挥着重要作用，层序界面的划分和地层对比在地层评价工作中占有重要地位，正逐渐由定性走向定量、并由宏观走向微观，其在砂岩中有着广阔的应用前景。

常用的时频分析方法，可以有效的将非平稳信号在时间域与频率域间相互转换，因此，被广泛的应用于勘探领域的层序划分中，其是岩性油气藏勘探的重要技术手段。根据沉积分选理论：正旋回时，随着深度的增加，沉积物粒度逐渐变大，地层厚度、速度和密度也相应增大；而反旋回时，随着深度的增加，沉积物粒度逐渐变小，地层厚度、速度和密度也相应减小。相关研究表明，正旋回在时频谱上表现为频率减小、振幅减弱、能量团减弱；而反旋回则表现为频率增加、振幅增加、能量团增加；正—反旋回表现为频率先减小后增加、振幅先减弱后增强、能量团先减弱后增强；反—正旋回表现为频率先增加后减小、振幅先增强后减弱、能量团先增强后减弱。

由于地震资料相对于录井、测井信息而言，其纵向分辨率较低，不能进行精细层序划分。因此，目前主要使用地震资料进行三级层序以内的划分和中长期旋回的识别划分，四、五级层序或更小级别界面划分需要借助录井、测井等资料，但是这些界面在地震剖面上很难进行准确的标定和追踪。且由于常规的三级层序格架难以满足复合区岩性圈闭识别、描述、评价与优选的精度要求，所以，建立高频层序格架的约束尤为重要；同时，高精度层序地层学的研究对岩性油气藏勘探开发意义重大。

常用的时频分析法将地震信号从时间域变换到频率域，并分解成不同的频率成分分量，由于不同频率的地震信号对各种地质构造的响应特征差异很大，因此，其可以通过频率特征的变化差异来有效地区分地质构造和划分层序。

发明内容

本发明目的在于提供一种基于变分模态分解的三维沉积旋回模型构建方法,以解决四、五级层序或更小级别界面划分需要借助录井、测井等资料，但这些层序界面在地震剖面上，很难进行准确的标定和追踪的技术问题。

为实现上述目的，本发明的基于变分模态分解的三维沉积旋回模型构建方法的具体技术方案如下：

一种基于变分模态分解的三维沉积旋回模型构建方法，包括以下步骤：

第一步：输入叠后地震数据，选定计算线道号；

第二步：基于变分模态分解，将信号分解成多个有限带宽的信号分量，并通过计算得到地震数据的瞬时相位；

第三步：对瞬时相位体中的各道数据分别进行S变换，得到每道时域数据对应的频域道集数据；

第四步：调节γ和ρ使得频域道集数据达到频率分辨率最高、时间分辨率最低且横向连续性较好；

第五步：对频域道集数据进行降维处理，计算每道频域数据的振幅最大值对应的频率，就得到频域道集数据对应的峰值频率；

第六步：将所得峰值频率曲线逐道放入对应线道号的地震数据中进行替换，完成沉积旋回体的构建。

进一步，所述第二步中，具体的计算过程如下：

经过VMD分解的IMF信号能够定义为：

其中，A

为了让每个通过模态分解得到的分量，其对应的频带都在有效频带范围内，需要进行以下几个步骤：

①用Hilbert变换求取每个IMF分量对应的单边频谱；

②将各IMF分量的中心频带通过指数函数运算整合到基带；

③通过L2范数梯度的平方，计算每个IMF分量的带宽。

进一步，所述第二步中，为了将约束变分问题转化为非约束变分问题，引入二次项惩罚参数α和增广拉格朗日方程来优化约束结果，如下式：

其中，u

其具体步骤如下：

1.定义u

2.n＝n+1为主循环，k从1到k-1为第一个内循环，更新求解uk，直到k＝K结束；

⒊k从1到k-1为第二个内循环，更新求解w

4.对于所有w

其中，τ为噪声容限参数。

5.给定精度ε>0，重复上述步骤，直到满足条件；

因此，每一个IMF分量都可以当作一个单一模态，每个模态的瞬时振幅、瞬时频率、瞬时相位可以通过一下公式计算：

其中，A(t)、F(t)和θ(t)分别为瞬时振幅、瞬时频率和瞬时相位，R(·)和I(·)为信号的实部和虚部，R′和I′为R(·)和I(·)对时间t的导数。

进一步，所述第三步中，改进的广义S变换具体过程为：

1.给定一个信号x(t)，通过傅里叶变换得到它的频率域表达X(ω)：

式中ω为变量t的频率域表达；

傅里叶变换具有可逆性，傅里叶逆变换可以将信号从频率域转换到时间域。令频率域信号函数为X(ω)，则时间域信号函数x(t)的表达式为：

一般情况下，X(ω)为复数，因此，X(ω)还可以表示为：

X(ω)＝A(ω)exp(iφ(ω)) (10)

其中，A(ω)与φ(ω)分别代表振幅谱和相位谱，计算式为：

其中，X

X(ω)＝X

其中：

2.对于信号x(t)在某一时间点τ的频率分布，可以求取以τ为中心的一个时窗I

其中，|I

令L

其中，A代表幅值，r代表角频，

信号x(t)的短时傅里叶变换表达式可表示为：

g(t-τ)为窗函数，在纵向上随着时间τ的变化，对信号做分段时频分析；

将上式做离散化处理，得到式(10)，

其中，x(t)为信号，a为时间；w(a)为窗函数；短时傅里叶变换同样存在缺点，因为时窗长度是固定的，因此一旦时窗确定，不管时间τ怎样变换，截取的信号长度都不会变，所以该方法对信号的处理不够灵活；

⒊若R平面上的平方可积函数空间为

其中：

小波变换在应用中也有一定的限制条件，小波基应满足以下要求：

其中，ψ(ω)是

其中：

当

4.令h(t)∈L

其中，τ是时间，f为频率；ST(τ,f)是h(t)的S变换结果；

S变换中，窗函数可表式为：

Gaussian函数表达式定义为：

通过对τ积分得到信号h(t)的频率的H(f)：

S变换是完全可逆的，表达式为：

5.添加两个控制参数：

则函数h(t)的改进广义S变换表达式为：

公式(27)是在S变换基础上，对Gaussian函数添加时间和频率控制参数改进的广义S变换；其中，f是频率，t是时间，τ是时间上的变化参数，i是虚数符，γ和ρ分别对应对时间和频率窗口的改动。

进一步，所述第六步中，再对沉积旋回体进行横向连续性优化，在优化后的沉积旋回体上进行地层划分、层序界面解释的研究。

本发明的基于变分模态分解的三维沉积旋回模型构建方法具有以下优点：

1、本发明由于采用“一种自适应、完全非递归的模态变分和信号处理的方法”即：VMD的算法，而VMD算法具有自适应的特点，它可以通过非递归的方式，将信号分解为多个准正交固有模态函数；其与模态分解，即：EMD算法所不同的是，VMD算法在分解信号的过程中，有严格的数学推导做支撑；因此，基于VMD算法可以获得更加准确的地震数据的瞬时相位分量。

2、本发明利用广义S变换灵活调节时窗、具有高时频分辨率与聚焦性的优势，其将瞬时相位数据转化为频域道集数据，从频域进行沉积旋回识别和层序界面划分。

3、本发明沉积旋回体剖面，反映了纵横向频率变化，揭示了地层纵向和横向的沉积变化规律；且横向上具有较强的可追踪性，对层序界面的解释具有一定的指导作用。

附图说明

图1为本发明的流程示意图；

图2-1为本发明的单道S变换地震数据示意图(其为屏幕上的实际图形)；

图2-2为本发明的单道S短时傅里叶逆变换示意图(其为屏幕上的实际图形)；

图2-3为本发明的单道S三角滤波变换示意图(其为屏幕上的实际图形)；

图2-4为本发明的单道广义S变换示意图(其为屏幕上的实际图形)；

图2-5为本发明的单道VMD算法瞬时相位示意图(其为屏幕上的实际图形)；

图3-1为本发明的第一沉积旋回体剖面波形显示示意图(其为屏幕上的实际图形)；

图3-2为本发明的第一沉积旋回体峰值频率剖面示意图(其为屏幕上的实际图形)；

图3-3为本发明的第一沉积旋回体瞬时相位时频谱剖面示意图(其为屏幕上的实际图形)；

图3-4为本发明的第二沉积旋回体剖面波形显示示意图(其为屏幕上的实际图形)；

图3-5为本发明的第二沉积旋回体峰值频率剖面示意图(其为屏幕上的实际图形)；

图3-6为本发明的第二沉积旋回体瞬时相位时频谱剖面示意图(其为屏幕上的实际图形)；

图4为本发明的沉积旋回体剖面密度显示示意图(其为屏幕上的实际图形)。

具体实施方式

为了更好地了解本发明的目的、结构及功能，下面结合附图，对本发明一种基于变分模态分解的三维沉积旋回模型构建方法做进一步详细的描述。

如图1所示，本发明包括以下步骤：

第一步：输入叠后地震数据，选定计算线道号；

第二步：基于变分模态分解(Variational Mode Decomposition)，将信号分解成多个有限带宽的信号分量，并通过计算得到地震数据的瞬时相位，具体的计算过程如下：

经过VMD分解的IMF信号可定义为：

其中，Ak(t)为瞬时振幅，

为了让每个通过模态分解得到的分量，其对应的频带都在有效频带范围内，需要进行以下几个步骤：

(1)用Hilbert变换求取每个IMF分量对应的单边频谱；

(2)将各IMF分量的中心频带通过指数函数运算整合到基带；

(3)通过L2范数梯度的平方，计算每个IMF分量的带宽。

为了将约束变分问题转化为非约束变分问题，引入二次项惩罚参数α和增广拉格朗日方程来优化约束结果，如下式：

其中，u

其具体步骤如下：

(1)定义u

(2)n＝n+1为主循环。k从1到k-1为第一个内循环，更新求解uk，直到k＝K结束。

(3)k从1到k-1为第二个内循环，更新求解w

(4)对于所有w

其中，τ为噪声容限参数。

(5)给定精度ε>0，重复上述步骤，直到满足条件。

因此，每一个IMF分量都可以当作一个单一模态，每个模态的瞬时振幅、瞬时频率、瞬时相位可以通过一下公式计算：

其中，A(t)、F(t)和θ(t)分别为瞬时振幅、瞬时频率和瞬时相位，R(·)和I(·)为信号的实部和虚部，R′和I′为R(·)和I(·)对时间t的导数。

第三步：对瞬时相位体中的各道数据分别进行S变换，得到每道时域数据对应的频域道集数据，改进的广义S变换具体过程为：

1)给定一个信号x(t)，通过傅里叶变换得到它的频率域表达X(ω)：

式中ω为变量t的频率域表达。

傅里叶变换具有可逆性，傅里叶逆变换可以将信号从频率域转换到时间域。令频率域信号函数为X(ω)，则时间域信号函数x(t)的表达式为：

一般情况下，X(ω)为复数，因此，X(ω)还可以表示为：

X(ω)＝A(ω)exp(iφ(ω)) (10)

其中，A(ω)与φ(ω)分别代表振幅谱和相位谱，计算式为：

其中，X

X(ω)＝X

其中：

2)由于傅里叶变换在信号处理时的平均效应，使得它在地震数据处理中无法准确分析纵向的信号特征。为了更好的分析地震信号对应的频率信息，短时傅里叶变换应运而生。短时傅里叶变换是在傅里叶变换的基础上增加矩形窗，在纵向上将地震信号截断，人为的将信号分成若干段，以得到若干时间段的频率信息。

对于信号x(t)在某一时间点τ的频率分布，可以求取以τ为中心的一个时窗I

其中，|I

在短时傅里叶变换中，窗口长度的大小影响了时频分辨率，大时窗对应着更高的频率分辨率，但是却会导致时间分辨率的降低；同样的，当采用小时窗时信号所对应的时间分辨率会提高，频率分辨率会降低。

令L

其中，A代表幅值，r代表角频，

信号x(t)的短时傅里叶变换表达式可表示为：

g(t-τ)为窗函数，在纵向上随着时间τ的变化，对信号做分段时频分析。

将上式做离散化处理，得到式(10)，

其中x(t)为信号，a为时间；w(a)为窗函数。短时傅里叶变换同样存在缺点，因为时窗长度是固定的，因此一旦时窗确定，不管时间τ怎样变换，截取的信号长度都不会变。所以该方法对信号的处理不够灵活。

3)在对非平稳信号的处理中，短时傅里叶变换存在明显不足，其原因在于窗函数固定，为了更好的分析信号的频率域信息，能变换不同基小波的Gabor变换被提出。虽然小波变换也采用加窗的思想，但Gabor变换改进了短时傅里叶变换的不足之处，窗函数可以根据不同的信号特征改变，表现出了很好的自适应性

若R平面上的平方可积函数空间为

其中：

小波变换在应用中也有一定的限制条件，小波基应满足以下要求：

其中，ψ(ω)是

其中：

为了使小波变化处理后的时频结果达到最佳的效果，可以针对不同的信号选取不同的小波基，并且对小波基函数的时窗长度做灵活调整。当

4)虽然，小波变换中引入了尺度控制函数以及时间控制函数，解决了短时傅里叶变换时窗固定不变的缺点，让算法变得灵活。但是，小波变换不是严格意义上的时间-频率谱，其本质反映的是时间与尺度之间的关系。同时，因为算法中需要设定基小波，所以，该算法始终还是会受时窗的影响。针对小波变换所存在的问题，Stockwell(1996)等人提出了S变换。在S变换中，简谐波控制时间尺度上的伸缩，Gaussian函数，同时控制伸缩和平移。这些改变让S变换可以根据频率的变化自适应变化，并且，不受容许性条件约束。

令h(t)∈L

其中，τ是时间，f为频率；ST(τ,f)是h(t)的S变换结果。

S变换中，窗函数可表式为：

Gaussian函数表达式定义为：

通过对τ积分得到信号h(t)的频率的H(f)：

S变换是完全可逆的，表达式为：

5)在S变换中，简谐波函数和Gaussian函数相乘得到了基小波，虽然，在尺度和时间上可以任意变化，但是，它的基函数形状是固定不变的。使得S变换在实际应用中也会受限。针对S变换的这一缺点，在Gaussian函数中做进一步改进，添加两个控制参数：

则函数h(t)的改进广义S变换表达式为：

公式(27)是在S变换基础上，对Gaussian函数添加时间和频率控制参数改进的广义S变换。其中，f是频率，t是时间，τ是时间上的变化参数，i是虚数符，γ和ρ分别对应对时间和频率窗口的改动。

第四步：调节γ和ρ使得频域道集数据达到频率分辨率最高、时间分辨率最低且横向连续性较好。

第五步：如图3-1至图3-6所示，对频域道集数据进行降维处理，计算每道频域数据的振幅最大值对应的频率，就得到频域道集数据对应的峰值频率。

第六步：如图4所示，将所得峰值频率曲线逐道放入对应线道号的地震数据中进行替换，完成沉积旋回体的构建。然后对沉积旋回体进行横向连续性优化，在优化后的沉积旋回体上进行地层划分、层序界面解释等研究。

上述未作说明的技术为现有技术，故不再赘述。

可以理解，本发明是通过一些实施例进行描述的，本领域技术人员知悉的，在不脱离本发明的精神和范围的情况下，可以对这些特征和实施例进行各种改变或等效替换。另外，在本发明的教导下，可以对这些特征和实施例进行修改以适应具体的情况及材料而不会脱离本发明的精神和范围。因此，本发明不受此处所公开的具体实施例的限制，所有落入本申请的权利要求范围内的实施例都属于本发明所保护的范围内。