CN102508295B - 一种地震地层厚度变化分析方法 - Google Patents

Abstract

本发明公开了一种地震地层厚度变化分析方法，该方法利用基于经验模式分解的高分辨率瞬时频率来分析地层厚度变化。首先对地震信号进行经验模式分解过程，得到多级内蕴模式函数；然后计算内蕴模式函数的瞬时频率，根据该瞬时频率的变化分析地层厚度变化。本发明在经验模式分解过程中采用三次B样条作为插值函数，采用基于连续小波变换的带有阻尼因子的瞬时频率计算方法计算瞬时频率，使得得到的内蕴模式函数和瞬时频率剖面的横向连续性好，假象少，符合实际，且抗噪性能好。本发明有助于提高地震地层厚度变化分析的有效性和可靠性，适用于大规模的实际地震资料处理与解释，为储层预测提供有力的依据，提高油气田勘探开发的精度。

Description

一种地震地层厚度变化分析方法

技术领域：

本发明是一种有关于地震勘探中的地震处理和解释方法，尤其是有关于一种基于经验模式分解和瞬时频率的高分辨率地震地层厚度变化分析方法。

背景技术：

目前用于地震地层厚度变化分析的瞬时属性主要是指基于常规希尔伯特变换的瞬时属性，包括瞬时振幅、瞬时相位以及瞬时频率等。其中瞬时频率可用于地层厚度变化分析以及异常分析。地震同相轴通常是由地下的一组反射界面反射的地震波合成的。当这组反射结构的厚度和岩性逐渐变化时，对应的合成地震波的瞬时频率也随之变化，当有尖灭或者岩性边界存在，瞬时频率将急剧变化。因此可以用瞬时频率来分析地层厚度变化。

现有技术1：

对地震数据的每一道地震信号重复步骤(1)-(3)：

(1)计算原始地震道信号的希尔伯特变换；

(2)构建一个解析信号：将原始地震道信号作为解析信号的实部，将它的希尔伯特变换作为解析信号的虚部；

(3)计算瞬时属性：解析信号的虚部除以实部，然后取反正切，得到的角度称为瞬时相位；对瞬时相位求导数可以得到瞬时频率。

(4)根据得到的地震瞬时频率剖面，分析地层厚度的横向、纵向变化；

现有技术1的缺点：

(1)基于常规希尔伯特变换的瞬时频率要求信号是窄带信号，而实际的地震信号很难满足这一要求，因此基于常规希尔伯特变换的瞬时频率有时没有物理意义，与实际不相符，它更多的表现为一种数学计算；

(2)基于常规希尔伯特变换的瞬时频率分辨率较低，对薄层的厚度变化不敏感，且受噪声影响较大，因此对薄互层结构很难判断其厚度变化。

现有技术2：

对地震道进行经验模式分解和希尔伯特变换，得到地震道的时频谱，根据地震道的时频谱分析地层的厚度变化；

现有技术2的缺点：

(1)针对的是单道地震记录，只能分析地层的纵向厚度变化；

(2)经验模式分解的过程中，采用三次样条作为插值函数存在过冲和欠冲问题，拟合的局部极大值和局部极小值包络受数据采样点波动或小扰动影响极大，可能会得到偏差较大的局部极大值和极小值包络，使得相邻道的经验模式分解结果变化较大，进而导致剖面的横向连续性很差，不能从整个剖面上分析地层的横向厚度变化，很难在实际地震信号处理中得到大量应用。

(3)经验模式分解过程中，采用定义内蕴模式函数的基本条件来判断筛选出的信号是否为内蕴模式函数，导致计算效率较低，计算量较大；内蕴模式函数满足的基本条件是：

①函数的局部极值点个数和零点个数相等或者至多相差1；

②在任一点上，函数的局部极大值包络和局部极小值包络的平均值等于零。

发明内容：

为了解决上述现有技术所存在的问题，根据经验模式分解可以对信号进行多分辨率分解的特性，本发明提供一种高分辨率地震地层厚度变化分析方法，对地震信号进行经验模式分解得到内蕴模式函数，然后对各级内蕴模式函数，计算它的基于连续小波变换的瞬时频率，依据该瞬时频率变化分析地层的横向和纵向厚度变化。该方法采用三次B样条函数直接对局部极大值点和局部极小值点插值拟合平均值，利用标准偏差准则来判断筛选得到的信号是否满足内蕴模式函数的要求。计算瞬时频率时，该方法采用基于连续小波变换的带有阻尼因子的瞬时频率计算公式。

本发明的目的是通过以下技术方案来解决的，一种地震地层厚度分析方法：

对地震记录的每一道地震信号重复步骤01-步骤04：

步骤01：对每一道地震信号进行经验模式分解，分解为多级内蕴模式函数分量；

对每一级内蕴模式函数重复步骤02-04：

步骤02：对内蕴模式函数，计算它的连续小波变换；

步骤03：对内蕴模式函数构建对应的解析信号；

步骤04：计算内蕴模式函数对应的瞬时频率；

步骤05：根据得到的各级内蕴模式函数的瞬时频率变化，分析地震的地层厚度变化；

根据本发明的实施例，步骤01所述的经验模式分解过程是一个不断筛选的过程，首先将原始地震信号作为输入信号

筛选过程包括如下步骤：

(A).基于输入信号的局部极大值，极小值点，采用三次B样条函数插值拟合得到平均值；

(B).从输入信号中减去平均值，得到一个新的信号；

(C).步骤(B)得到的新信号如果满足给定的条件，则可以视作一级内蕴模式函数，并执行步骤(D)，否则，将步骤(B)得到的新信号作为输入信号，回到步骤(A)；

(D).提取内蕴模式函数分量，将原始地震信号减去已得各级内蕴模式函数的信号记为残余信号，若其呈现单调趋势或者值很小时，经验模式分解过程结束，否则将其视作为下一次筛选过程的输入信号，回到步骤(A)。

根据本发明的实施例，上述步骤(A)中对输入信号的局部极大值、极小值点采用三次B样条函数插值拟合平均值时，采用四点镜像取反的方法处理端点处的平均值拟合，以输入信号左右两个时间端点为镜像点，对输入信号及其局部极大值、极小值点在左右两端分别进行镜像延拓，并对延拓的信号和局部极大值、极小值点的值取反，然后再将原输入信号的局部极大值、极小值点序列左右各增加四个局部极值点，做为新的局部极大值、极小值点序列，利用该序列拟合平均值；

根据本发明的实施例，上述步骤(C)中新信号满足的给定条件是指新信号和对应输入信号的标准方差小于0.05。

根据本发明的实施例，步骤02中各级内蕴模式函数的连续小波变换由以下关系式决定：

$S_n(b,a)=\frac{1}{a}{\∫}_{T_1}^{T_2}{s}_n\left(t\right)\stackrel{\‾}{g}\left(\frac{t-b}{a}\right)\mathrm{dt},$

式中，s_n(t)表示第n级内蕴模式函数，n＝1，2，3，…N，N为总分解级数，t为时间，t属于实数集合，最小值为T₁，最大值为T₂，尺度因子a＞0，b为时移，b属于实数集合，b的最小值为T₁，最大值为T₂；S_n(b，a)为s_n(t)的连续小波变换，为小波函数

的复共轭， $g\left(\frac{t-b}{a}\right)=e^{i{\mathrm{\ω}}_0\frac{t-b}{a}}{e}^{-\frac{{(t-b)}^2}{2a^2}}$ 是对基本小波函数 $g\left(t\right)=e^{i{\mathrm{\ω}}_{0}t}{e}^{-\frac{t^2}{2}}$ 经过伸缩a和时移b得到的，基本频率ω₀＝6.0， $i=\sqrt{-1}$ 为虚数单位。

根据本发明的实施例，步骤03构建内蕴模式函数s_n(t)的解析信号的方法是将内蕴模式函数s_n(t)作为解析信号的实部，解析信号的虚部

由内蕴模式函数s_n(t)的连续小波变换结果决定：

$s_n^{*}\left(t\right)=\mathrm{Im}\left(\frac{1}{C_g}{\∫}_0^{\∞}{S}_n(t,a)\frac{\mathrm{da}}{a}\right);$

其中，Im(·)表示对括号内的表达式取虚部，

是基本小波函数g(t)的实部的傅里叶变换，ω为角频率；基于t和b的一致性，此处将步骤02的S_n(b，a)用S_n(t，a)表示，即将变量b用t代替；对应的解析信号Z_n(t)表示为：

$Z_n\left(t\right)=s_n\left(t\right)+i{s}_n^{*}\left(t\right);$

根据本发明的实施例，步骤04中内蕴模式函数的瞬时频率由以下关系式决定：

$f_n\left(t\right)=\frac{1}{2\mathrm{\π}}\frac{s_n\left(t\right)\frac{{\mathrm{ds}}_n^{*}\left(t\right)}{\mathrm{dt}}-s_n^{*}\left(t\right)\frac{d{s}_n\left(t\right)}{\mathrm{dt}}}{e^2\left(t\right)+\mathrm{\ϵ}{e}_{\max }^2},$

其中，s_n(t)表示第n级内蕴模式函数，f_n(t)表示第n级内蕴模式函数s_n(t)的瞬时频率， $e^2\left(t\right)={s_n}^2\left(t\right)+{\left(s_n^{*}\left(t\right)\right)}^2,$

表示e²(t)的最大值，阻尼因子ε取为0.01或0.005；

根据本发明的实施例，步骤05中分析地层厚度变化的主要依据是：对应瞬时频率值小的地层较厚，对应瞬时频率值较大的地层薄；当瞬时频率由大变小时，说明地层由薄变厚，当瞬时频率由小变大时，说明地层由厚变薄。

本发明在经验模式分解中采用三次B样条函数直接对局部极大值点和局部极小值点插值拟合平均值，能够减少经验模式过程中的过冲和欠冲以及减弱内蕴模式函数受假象的干扰程度，可以得到更为光滑、准确的内蕴模式函数分量，进而地震剖面的横向连续性也变得更好；采用筛选得到的新信号和对应输入信号的标准方差是否小于0.05的准则判定筛选过程的结束与否既能有效的提高计算效率，又不会影响筛选出的内蕴模式函数的精细程度。计算瞬时频率时，采用基于连续小波变换的带有阻尼因子的瞬时频率计算方法，可以有效的减弱瞬时频率受噪声的干扰程度。利用本发明得到的瞬时频率比较稳定、具有较为明确的物理意义，分辨率较高，抗噪性较好，更能有效的反映地层的横向、纵向厚度变化，提高地震地层厚度分析的有效性和可靠性，适用于大规模的实际地震资料处理及应用。

附图说明：

图1为本发明的地震地层厚度变化分析方法流程图；

图2为经验模式分解流程图；

图3为经验模式分解过程中端点处的极大值、极小值点处理方法-四点镜像取反方法示意图；

图4A-图4G为本发明第一实施例单道地震记录模型厚度变化分析：A等幅度反射系数序列；B合成地震记录；C各级内蕴模式函数分量；D各级内蕴模式函数的傅里叶幅度谱；E瞬时频率；F幅度不等的反射系数序列；G瞬时频率；

图5A-图5G为本发明第二实施例的地层厚度变化分析：A二维地震资料；B采用现有技术2得到的第一级内蕴模式函数；C采用现有技术2得到的第二级内蕴模式函数；D采用本专利方法得到的第一级内蕴模式函数；E采用本专利方法得到的第二级内蕴模式函数；F采用现有技术2得到的第一级内蕴模式函数的瞬时频率；G采用本专利方法得到的第一级内蕴模式函数的瞬时频率；

图6A-图6B为本发明第二实施例中基于本专利方法的单道地震记录的地层厚度变化分析：A图5G中W1井处的地震道分析；B图5G中W2井处的地震道分析；

具体实施方式：

下面结合附图和具体实施方式对本发明做进一步详细的说明。

本发明提供一种基于经验模式分解和瞬时频率的高分辨率地震地层厚度分析方法，其是根据各级内蕴模式函数的瞬时频率变化来分析地震地层厚度变化趋势。

如图1所示，本发明所提供的地震地层厚度分析方法，具体通过如下步骤实施：

准备地震记录，地震记录可以为一维、二维或三维的地震记录，对地震记录的每一道地震信号s(t)重复步骤1-4：

步骤1：对s(t)进行经验模式分解，分解为N级内蕴模式函数s_n(t)，n＝1，2，3，…N，在实际应用中，一般取N＝3；

参照图2，经验模式分解过程包括如下步骤：

首先将原始地震道信号s(t)作为输入信号，令内蕴模式函数级数n＝1，筛选次数k＝1，令h₁₀(t)＝s(t)；

(1)找出输入信号的局部极大值，极小值点，采用三次B样条函数插值拟合得到平均值m_nk(t)；

(2)从输入信号中减去平均值m_nk(t)，得到一个新的信号h_nk(t)＝h_n(k-1)(t)-m_nk(t)，h_n(k-1)(t)为本次筛选过程中(1)的输入信号；

(3)由(2)得到的新信号h_nk(t)如果满足设定的条件，则可以视作一级内蕴模式函数，并执行(4)，否则，将(2)得到的新信号作为输入信号，k＝k+1，即筛选次数增加1，回到(1)；

(4)令s_n(t)＝h_nk(t)，将原始地震信号s(t)减去已得各级内蕴模式函数的信号记为r(t)，若其呈现单调趋势或者值很小时，经验模式分解过程结束，跳出循环；否则将此r(t)视作为下一次筛选过程的输入信号，n＝n+1，即内蕴模式函数的级数增加1，重置筛选次数k＝1，h_n0(t)＝r(t)，回到(1)，重新开始筛选。

上述(1)中对输入信号的局部极大值、极小值点拟合平均值m_nk(t)时，需要考虑信号最小时间和最大时间点处的端点处理问题。在本专利中采用四点镜像取反的方法处理端点处的平均值拟合，从而抑制端点效应，图3所示为四点镜像取反方法的示意图，可以此类推：

①设信号的局部极大值点、局部极小值点序列为{τ₀，τ₁，τ₂，τ₃，...τ_M-3，τ_M-2，τ_M-1，τ_M}，其中，τ₀，τ₁，τ₂，τ₃，...τ_M-3，τ_M-2，τ_M-1，τ_M表示局部极大值或者局部极小值点处的值大小，信号总的极值点数目为M+1，信号的左右端点一般不作为局部极大值点或局部极小值点；

②对输入信号及其局部极大值、极小值点在左右两端分别进行镜像延拓，左边以最小时间点T₁镜像对称延拓，右边以最大时间点T₂镜像对称延拓，并对左右延拓的信号和局部极大值和极小点处的值取反，图3所示为信号右端的延拓示意图，T₂右边为镜像延拓并取反的局部极大、极小值点，虚线为延拓取反的信号部分，左端的延拓方法以此类推；

③将局部极大值点、局部极小值点的值序列在信号左右端点两边各延拓四个点，延拓为{-τ₃，-τ₂，-τ₁，-τ₀，τ₀，τ₁，τ₂，τ₃，...τ_M-3，τ_M-2，τ_M-1，τ_M，-τ_M，-τ_M-1，-τ_M-2，-τ_M-3}，其中最左端-τ₃与τ₃关于最小时间点T₁镜像对称，大小取反，-τ₂与τ₂关于最小时间点T₁镜像对称，大小取反，-τ₁与τ₁关于最小时间点T₁镜像对称，大小取反，-τ₀与τ₀关于最小时间点T₁镜像对称，大小取反；最右端-τ_M-3与τ_M-3关于最大时间点T₂镜像对称，大小取反，-τ_M-2与τ_M-2关于最大时间点T₂镜像对称，大小取反，-τ_M-1与τ_M-1关于最大时间点T₂镜像对称，大小取反，-τ_M与τ_M关于最大时间点T₂镜像对称，大小取反；

④基于③中的局部极大值点、局部极小值点的延拓序列{-τ₃，-τ₂，-τ₁，-τ₀，τ₀，τ₁，τ₂，τ₃，...τ_M-3，τ_M-2，τ_M-1，τ_M，-τ_M，-τ_M-1，-τ_M-2，-τ_M-3}，利用三次B样条函数拟合平均值m_nk(t)。

上述(3)所涉及的新信号满足设定的条件是指新信号和对应输入信号的标准方差SD小于0.05，SD的表达式如下：

$\mathrm{SD}=\underset{t=T_1}{\overset{T_2}{\mathrm{\Σ}}}\left[\frac{{|h_{n(k-1)}\left(t\right)-h_{\mathrm{nk}}\left(t\right)|}^2}{h_{n(k-1)}^2\left(t\right)}\right]$

其中，T₁为最小时刻，T₂为最大时刻，判别标准选取为0.05，既不会过多增加计算量，又能保证得到的内蕴模式函数的精细程度。

对每一级内蕴模式函数s_n(t)重复步骤2-4，n＝1，2，3，…N：

步骤2：对内蕴模式函数s_n(t)，计算它的连续小波变换，连续小波变换定义如下：

$S_n(b,a)=\frac{1}{a}{\∫}_{T_1}^{T_2}{s}_n\left(t\right)\stackrel{\‾}{g}\left(\frac{t-b}{a}\right)\mathrm{dt},$

式中，s_n(t)表示第n级内蕴模式函数，n＝1，2，3，…N，t为时间，t属于实数集合，最小值为T₁，最大值为T₂，尺度因子a＞0，b为时移，b属于实数集合，b的最小值为T₁，最大值为T₂；S_n(b，a)为s_n(t)的连续小波变换，

为小波函数

的复共轭， $g\left(\frac{t-b}{a}\right)=e^{i{\mathrm{\ω}}_0\frac{t-b}{a}}{e}^{-\frac{{(t-b)}^2}{2a^2}}$ 是对基本小波函数 $g\left(t\right)=e^{i{\mathrm{\ω}}_{0}t}{e}^{-\frac{t^2}{2}}$ 经过伸缩a和时移b得到的，基本频率ω₀＝6.0， $i=\sqrt{-1}$ 为虚数单位。在实际计算中，T₁和T₂分别为输入信号的最小和最大时刻。

步骤3：对内蕴模式函数s_n(t)，利用以下关系式，构建对应的解析信号：

$Z_n\left(t\right)=s_n\left(t\right)+i{s}_n^{*}\left(t\right);$

式中，这里是将内蕴模式函数s_n(t)作为解析信号的实部，解析信号的虚部

由内蕴模式函数s_n(t)的连续小波变换结果决定：

$s_n^{*}\left(t\right)=\mathrm{Im}\left(\frac{1}{C_g}{\∫}_0^{\∞}{S}_n(t,a)\frac{\mathrm{da}}{a}\right);$

其中，Im(·)表示取虚部，

是g(t)的实部的傅里叶变换，ω为角频率；基于t和b的一致性，将步骤2的S_n(b，a)用S_n(t，a)表示，即将变量b用t代替；

步骤4：利用下式，计算各级内蕴模式函数对应的瞬时频率f_n(t)：

$f_n\left(t\right)=\frac{1}{2\mathrm{\π}}\frac{s_n\left(t\right)\frac{{\mathrm{ds}}_n^{*}\left(t\right)}{\mathrm{dt}}-s_n^{*}\left(t\right)\frac{d{s}_n\left(t\right)}{\mathrm{dt}}}{e^2\left(t\right)+\mathrm{\ϵ}{e}_{\max }^2},$

其中，s_n(t)表示第n级内蕴模式函数，f_n(t)表示第n级内蕴模式函数s_n(t)的瞬时频率， $e^2\left(t\right)={s_n}^2\left(t\right)+{\left(s_n^{*}\left(t\right)\right)}^2,$

表示e²(t)的最大值，阻尼因子ε取为0.01或0.005；

步骤5：根据得到的各级内蕴模式函数的瞬时频率变化，分析地震的地层厚度变化：对应瞬时频率值小的地层比较厚，对应瞬时频率值较大的地层比较薄；当瞬时频率由大变小时，说明地层由薄变厚，当瞬时频率由小变大时，说明地层由厚变薄。

本发明具有如下有益效果：

1)经验模式分解是一种自适应的局部信号分析方法，它非常适合分析具有非线性非平稳特征的地震信号；

2)相比于现有方法所使用的对信号局部极大值点和局部极小值点序列分别利用三次样条函数拟合包络，然后再计算平均值，本发明在经验模式分解过程中，采用三次B样条函数直接对局部极大值点和局部极小值点插值拟合得到平均值，能够有效的减少过冲和欠冲，使得得到的内蕴模式函数和瞬时频率剖面的横向连续性更好，假象较少，比较符合实际情况，可适用于大规模的实际地震资料处理；

3)采用基于连续小波变换的带有阻尼因子的瞬时频率计算方法来估算各级内蕴模式函数的瞬时频率，有效的提高了抗噪性能；

4)基于经验模式分解，可以得到地震信号的具有明确物理意义、多分辨率特性的瞬时频率，经验模式分解对原始信号细节的放大能力使得第一级内蕴模式函数相比于原始信号具有更高的时间分辨率，进而由第一级内蕴模式函数得到的瞬时频率具有更高的分辨率，它比直接由原始信号的希尔伯特变换得到的瞬时频率更能细致反映地层的厚度变化，更有利于地震解释和储层预测。

下面利用本发明提供的地震地层厚度分析方法，对(1)一维模型的地层厚度和(2)二维实际地震资料的气层厚度变化进行了分析估计；对于一维模型，将本发明提供的方法和现有技术1进行了对比，对于二维实际地震资料，主要针对第一、第二级内蕴模式函数和第一级内蕴模式函数的瞬时频率IFP1，将本发明提供的技术方法和现有技术2进行了对比讨论：

图4A-图4G为一维模型的地层厚度变化分析：

如图4A所示，其为反射系数序列，反射系数幅度相等，大小为0.25，正负交替，反射序列间隔由12ms减少到2ms再增加到12ms，模拟地层厚度由厚变薄再变厚，横坐标表示时间，单位为s。地震子波选取为45Hz的零相位Ricker子波，与图4A中的反射序列褶积得到的合成记录如图4B所示。利用三次B样条函数进行插值拟合，对图4B所示合成记录进行经验模式分解得到的前三级内蕴模式函数如图4C所示，由上到下依次为第一级内蕴模式函数IMF1、第二级内蕴模式函数IMF2和第三级内蕴模式函数IMF3，图4D为对应的傅里叶幅度谱；可以看出，第一级内蕴模式函数IMF1包含的高频分量最多，它的时间分辨率比原始信号要高，放大了信号的细节分量，而第三级内蕴模式函数则将向低频段移动，高频分量较少，时间分辨率较低，主要反映信号的概貌。利用本专利计算瞬时频率的方法，对图4C所示的三级内蕴模式函数计算的瞬时频率如图4E所示，其中实线为第一级内蕴模式函数的瞬时频率IFP1，实线加“*”为第二级内蕴模式函数的瞬时频率IFP2，实线加“●”为第三级内蕴模式函数的瞬时频率IFP3；另外，实线加“×”为利用现有技术1，基于常规希尔伯特变换得到的瞬时频率IFP。对比图4A和图4E可以看到，由第一级内蕴模式函数得到的瞬时频率IFP1的变化规律与模拟的地层厚度变化很一致：即地层较厚时，对应的瞬时频率较低，地层较薄时，对应的瞬时频率较高。另外，IFP2、IFP3的变化与地层厚度变化也基本一致，这有助于相互印证。而利用现有技术1得到的瞬时频率IFP在地层变薄的时候出现骤然下降的现象，即地层变薄时，IFP变小，如图4E中箭头所指的位置，这与地层厚度变化规律不相符合。

采用本专利的方法和现有技术1分析幅度不等的反射系数序列模型。如图4F所示，其为反射系数序列，反射系数幅度从0.25减小到0.05，再增加到0.25，正负交替，反射序列间隔由12ms减少到2ms再增加到12ms，模拟地层厚度由厚变薄再变厚，横坐标表示时间，单位为s。地震子波选取为45Hz的零相位Ricker子波，与图4F中的反射序列褶积得到合成记录。对该合成记录利用本专利方法和现有技术1可以得到图4G所示的瞬时频率，其中实线为第一级内蕴模式函数的瞬时频率IFP1，实线加“*”为第二级内蕴模式函数的瞬时频率IFP2，实线加“●”为第三级内蕴模式函数的瞬时频率IFP3；另外，实线加“×”为利用现有技术1，基于常规希尔伯特变换得到的瞬时频率IFP。对比可以看出，IFP1、IFP2和IFP3的变化规律与地层厚度变化规律吻合，而IFP的变化规律与地层厚度变化规律不相符合。因此我们可以利用本专利提供的方法，对信号进行经验模式分解，计算内蕴模式函数的瞬时频率，根据其变化分析地层厚度变化。

图5A-5G、图6A-6B为实际二维地震资料的地层厚度估计分析：

(1)如图5A所示，其为某油田的二维地震资料，共1435道，采样率1ms。已知含油气的目的层位于1300ms附近，故截取时间范围为1100ms至1500ms。图5B-5C为采用现有技术2，利用三次样条函数作为插值函数，对其进行经验模式分解，得到的前两级内蕴模式函数，图5B为第一级内蕴模式函数，图5C为第二级内蕴模式函数。可以看出，图5B和图5C中目的层的同相轴连续性很差，导致了很多假象断层，例如图中黑色箭头所指的位置；进一步由图5B中第一级内蕴模式函数得到的瞬时频率IFP1如图5F所示，白色表示IFP1值较小，颜色越黑表示IFP1值较大。可以看出，相邻两道地震道的IFP1值变化很大，IFP1剖面的横向连续性也很差，假象较多，这将导致对地层厚度横向变化的错误分析。

(2)图5D-5E为采用本发明提供的方法，利用三次B样条函数插值拟合平均值，对图5A的二维地震资料进行经验模式分解，得到的前两级内蕴模式函数，图5D为第一级内蕴模式函数，图5E为第二级内蕴模式函数。对图5D所示的第一级内蕴模式函数计算得到的瞬时频率IFP1如图5G所示，用灰度色标来表示IFP1值的大小，白色对应的IFP1值较小，黑色对应的IFP1值较大，W1、W2为两口已知井，图中给出了两口井的位置；可以看出IFP1剖面的横向连续性较好，横向变化较为稳定，与实际情况相符。进一步分析目的层处的地层厚度变化：从图5G可以看出，黑色箭头所指处的IFP1值较小，我们可以估计此处对应的地层较厚；向地震道号增大的方向，目的层的IFP1值逐渐变大，故可以推断横向上地层逐渐变薄，白色箭头所指处的IFP1值较大，我们可以估计此处对应的地层较薄。而在图5F上，根据IFP1值的大小无法分析出目的层地层厚度的横向变化。

为了验证基于本专利方法的厚度分析结果，我们进一步分析了W1、W2井对应的测井曲线及用本专利方法得到的瞬时频率IFP1。如图6A所示，其为W1井处的测井速度曲线和IFP1。可以看出，在1310ms处有一个较厚的气层，对应的IFP1的值较小。如图6B所示，其为W2井处的测井速度曲线和IFP1。可以看出，在1320ms处有一个较薄的气层，对应的IFP1值较大。这说明IFP1的变化与地层厚度变化相一致。

以上的模型和实际资料算例中，利用本发明的方法计算的瞬时频率与地层厚度变化有很好的对应关系，利用该瞬时频率变化可以有效的分析估计地层厚度变化，进一步用于储层预测和油气勘探。

最后需要说明的是，以上模型和实际资料算例对本发明的目的，技术方案以及有益效果提供了进一步的验证，这仅属于本发明的具体实施算例，并不用于限定本发明的保护范围，在本发明的精神和原则之内，所做的任何修改，改进或等同替换等，均应在本发明的保护范围内。

Claims (1)

1.一种地震地层厚度变化分析方法，其特征在于，该方法包括以下步骤：

对每一道地震信号重复步骤01-步骤04：

步骤01：对每一道地震信号进行经验模式分解过程，分解为n级内蕴模式函数分量，n=1,2,3,"N，N为总分解级数；

对每一级内蕴模式函数重复步骤02-步骤04：

所述的经验模式分解过程是一个不断筛选的过程，首先将原始地震信号作为输入信号，经验模式分解过程包括如下步骤：

(A).基于输入信号的局部极大值和极小值点，采用三次B样条函数插值拟合得到平均值；

(B).从输入信号中减去平均值，得到一个新的信号；

(C).步骤(B)得到的新信号如果满足给定的条件，则可以将其视作一级内蕴模式函数，并执行步骤(D)，否则，将步骤(B)得到的新信号作为输入信号，回到步骤(A)；

（D）.提取内蕴模式函数分量，将原始地震信号减去已得各级内蕴模式函数的信号记为残余信号，若其呈现单调趋势或者值很小时，经验模式分解过程结束，否则将其视作为下一次筛选过程的输入信号，回到步骤(A)；

所述步骤(A)中，对输入信号的局部极大值、极小值点采用三次B样条函数插值拟合平均值时，采用四点镜像取反的方法处理端点处的平均值拟合，以输入信号左右两个时间端点为镜像点，对输入信号及其局部极大值、极小值点在左右两端分别进行镜像延拓，并对延拓的信号和局部极大值、极小值点的值取反，然后再将原输入信号的局部极大值、极小值点序列左右各增加四个局部极值点，做为新的局部极大值、极小值点序列，利用该序列拟合平均值；

所述步骤(C)中，新信号满足的给定条件是指新信号和对应输入信号的标准方差小于0.05；

步骤02：对内蕴模式函数，计算它的连续小波变换；所述内蕴模式函数的连续小波变换由以下关系式决定：

$S_n(b,a)=\frac{1}{a}{\∫}_{T_1}^{T_2}{s}_n\left(t\right)\stackrel{\‾}{g}\left(\frac{t-b}{a}\right)\mathrm{dt},$

式中，s_n(t)表示第n级内蕴模式函数，n=1,2,3,…N，t为时间，t属于实数集合，最小值为T₁，最大值为T₂，尺度因子a>0，b为时移，b属于实数集合，b的最小值为T₁，最大值为T₂；S_n(b,a)为s_n(t)的连续小波变换，

为小波函数的复共轭， $g\left(\frac{t-b}{a}\right)=e^{{\mathrm{i\ω}}_0\frac{t-b}{a}}{e}^{-\frac{{(t-b)}^2}{2a^2}}$ 是对基本小波函数

经过伸缩a和时移b得到的，基本频率ω₀=6.0，

为虚数单位；

步骤03：对内蕴模式函数构建对应的解析信号；所述构建内蕴模式函数的解析信号的方法是将内蕴模式函数作为解析信号的实部，解析信号的虚部

由内蕴模式函数s_n(t)的连续小波变换结果决定：

$s_n^{*}\left(t\right)=\mathrm{Im}\left(\frac{1}{C_g}{\∫}_0^{\∞}{S}_n(t,a)\frac{\mathrm{da}}{a}\right);$

其中，Im(·)表示对括号内的表达式取虚部，

是g(t)的实部的傅里叶变换，ω为角频率；基于t和b的一致性，将S_n(b,a)用S_n(t,a)表示，即将变量b用t代替；对应的解析信号Z_n(t)表示为：

$Z_n\left(t\right)=s_n\left(t\right)+i{s}_n^{*}\left(t\right);$

步骤04：计算内蕴模式函数对应的瞬时频率；所述内蕴模式函数的瞬时频率由以下关系式决定：

$f_n\left(t\right)=\frac{1}{2\mathrm{\π}}\frac{s_n\left(t\right)\frac{d{s}_n^{*}\left(t\right)}{\mathrm{dt}}-s_n^{*}\left(t\right)\frac{d{s}_n\left(t\right)}{\mathrm{dt}}}{e^2\left(t\right)+\mathrm{\ϵ}{e}_{\max }^2},$

其中，s_n(t)表示第n级内蕴模式函数，f_n(t)表示第n级内蕴模式函数s_n(t)的瞬时频率，

表示e²(t)的最大值，阻尼因子ε取为0.01或0.005；

步骤05：根据得到的各级内蕴模式函数的瞬时频率变化，分析地震的地层厚度变化值，所述利用瞬时频率分析地层厚度变化的依据是：对应瞬时频率值小的地层为厚层，对应瞬时频率值大的地层为薄层；当瞬时频率由大变小时，说明地层由薄变厚，当瞬时频率由小变大时，说明地层由厚变薄。

Images