CN103364832A - 一种基于自适应最优核时频分布的地震衰减定性估计方法 - Google Patents

Abstract

本发明公开了一种基于自适应最优核时频分布的地震衰减定性估计方法，包括一下步骤，1）采集原始地震数据，2）拾取地震数据目标层的层位，3）确定高频频率f_H和低频频率f_L，4）对叠后地震数据体中每道信号计算自适应最优核时频分布并得到高频分量和低频分量，5）确定修正因子，6）估计地震衰减。该方法首先利用自适应最优核时频分布对地震数据进行时频分解，得到地震数据的高频分量和低频分量，然后利用低频分量和高频分量的差来定性估计地震衰减。与常用的基于时频分布的地震衰减定性估计方法相比，本发明具有较高的时频分辨率，能够更为精细地刻画地震衰减，此外，该技术方案易于实现，可操作性强。

Description

一种基于自适应最优核时频分布的地震衰减定性估计方法

技术领域

本发明属于地球物理勘探领域，涉及一种定性地震衰减估计方法，特别涉及一种基于自适应最优核时频分布的地震衰减定性估计方法。

背景技术

地震波在地下介质中传播，散射和吸收衰减特性等均与频率有关，因此时频分析是地震信号分析的重要手段。Fourier变换已在地震信号处理中取得了很多应用，例如：去噪、频谱分析等，所以也是目前商业处理软件必备的技术模块。加窗傅立叶变换用移动的窗函数来截取信号，然后对分析窗内的信号做傅立叶变换。加窗傅立叶变换高频处和低频处的时间-频率联合分辨率是固定的。连续小波变换在低频处使用大尺度来获得较宽的窗函数，在高频处使用小尺度来获得较窄的窗函数，因此在高低频有不同的时间-频率联合分辨率。连续小波变换利用尺度来控制窗的宽窄，获得的是时间-尺度分布，而尺度和频率并无直接的对应关系。S变换直接利用频率控制窗函数的宽窄，既保留了小波变换的优点，又克服了小波变换尺度和频率没有对应关系的缺点。上述加窗傅立叶变换、小波变换及S变换均为线性变换，其时间-频率联合分辨率均受限于不确定性原理。Wigner-Ville分布具有最高的时间-频率联合分辨率，但是受交叉项干扰严重，因此在实际应用中受限。Cohen类时频分布利用不同的核函数对Wigner-Ville分布进行滤波以抑制交叉项干扰，得到不同的时频分布，因此Cohen类时频分布的结果完全取决于核函数。许多学者设计了固定型核函数，因此只能适合分析特定种类的信号。Baraniuk和Jones提出依赖信号的最优核时频分布，在整个时间持续期内根据信号对核函数进行迭代优化，获得一个适应于信号的时频分布。Jones和Baraniuk对依赖信号的时频分布进行改进，提出自适应最优核时频分布。该分布的核函数是随时间的变化而变化的，故而自适应最优核时频分布具有较高的时间-频率联合分辨率以及良好的交叉项抑制能力。多位学者将该分布用于分析线性调频信号，取得了较好的效果。此外，常用的时频分析方法还有匹配追踪算法及其改进算法等，并得到了广泛的发展和应用。

Taner等人于1979年观察到地震波在穿过含气储层以后视频率会降低，而这种视频率的异常可以用来指示含气储层。为了提取这种异常，需要时频分析的方法。多种时频分析方法已被用来检测这种高低频分量的不同衰减特性，并将其用于直接碳氢储层指示：匹配追踪类方法被广泛用于直接碳氢储层检测以及碳氢引起的低频阴影检测，但是匹配追踪类方法的结果对字典依赖性较强且运算量巨大；陈文超和高静怀等人研究含气储层对地震波的影响，并利用改进的匹配地震子波作为母小波对地震波衰减特性进行分析，在某油田的实际应用取得了明显效果；陈学华等人分析了低频阴影的机理，利用广义S变换对某地区三维实际资料进行频谱分解，能够检测油气储层的低频阴影，刻画油气储层的空间展布。然而稳定的、时间-频率联合分辨率高的时频工具仍是关键。匹配追踪类时频分析方法具有很高的时间-频率联合分辨率，但是计算极为耗时且结果不稳定，此外结果受所选取的字典影响。常用的连续小波变换及广义S变换结果很稳定，但其时间-频率联合分辨率受限于不确定性原理，并且在高频处和低频处的时间-频率联合分辨率不，因此，影响分析结果。

发明内容

针对上述缺陷或不足，本发明提供了一种基于自适应最优核时频分布的地震衰减定性估计方法，目的在于克服上述时频分析方法的缺点，利用高频\低频处具有相同分辨率且时频分辨率能够突破不确定性原理的自适应最优核时频分布作为时频分析工具，根据地震波在穿过含气储层后高频分量和低频分量衰减不同的原理，实现对地震衰减定性进行估计。

为达到以上目的，本发明的技术方案为：

包括以下步骤：

1）采集原始地震数据

采集原始地震数据的道号、线号、以及采样点号时间，对采集到的原始地震数据进行预处理，得到叠后地震数据体，记为S(x,y,t)，其中x表示道号，y表示线号，t表示采样点号时间；

2）拾取地震数据目标层的层位

拾取目标层位及目标层上方的层位，记目标层位为H_T(x,y)，目标层位上面的层位为H_A(x,y)；

3）确定高频频率f_H和低频频率f_L

首先，对目的层邻近的地震数据进行频谱分析，然后选择高频频率fH和低频频率fL；

4）对叠后地震数据体中每道信号计算自适应最优核时频分布并得到高频分量AOKTFR(x,y,t,fH)和低频分量AOKTFR(x,y,t,fL)；

5）确定修正因子α（x，y）

$\mathrm{\α}(x,y)=\frac{\mathrm{AOKTFR}(x,y,H_A(x,y),f_L)}{\mathrm{AOKTFR}(x,y,H_A(x,y),f_H)}---\left(1\right)$

6）估计地震衰减

利用得到的修正因子α（x，y）、高频分量AOKTFR(x,y,t,f_H)以及低频分量AOKTFR(x,y,t,f_L)估计目标层位H_T(x,y)附近的衰减SA（x,y,t）：

SA(x,y,t)=AOKTFR(x,y,t,f_L)-α(x,y)AOKTFR(x,y,t,f_H)     （2）

所述确定高频频率f_H和低频频率f_L具体包括：

首先，对目的层邻近地震数据进行频谱分析，得到振幅谱，然后，选取高频的幅度与低频的幅度接近、高频与低频之间的频率间隔大、并且高频与低频都在地震子波的频带范围之内的高频频率和低频频率为所需高频频率f_H和低频频率f_L。

所述对叠后地震数据体中每道信号计算自适应最优核时频分布并得到高频分量AOKTFR(x,y,t,f_H)和低频分量AOKTFR(x,y,t,f_L)具体为：

定义信号s(t)的对称短时模糊函数SSTAF为：

$\begin{array}{c}\mathrm{SSTAF}(t;\mathrm{\θ},\mathrm{\τ})=\frac{1}{2\mathrm{\π}}{\∫}_{-\∞}^{+\∞}{s}^{*}(u-\mathrm{\τ}/2)w^{*}(u-t-\mathrm{\τ}/2)s(u+\mathrm{\τ}/2)w(u-t+\mathrm{\τ}/2)e^{-\mathrm{j\θu}}\mathrm{du}\\ =\frac{1}{2}{\∫}_{t-T+|\mathrm{\τ}/2|}^{t+T-|\mathrm{\τ}/2|}{s}^{*}(u-\mathrm{\τ}/2)s(u+\mathrm{\τ}/2)e^{-\mathrm{j\θu}}\mathrm{du}\end{array}---\left(3\right)$

其中，w是宽度为2T的对称矩形窗函数,τ为时延，θ为多普勒频率。利用SSTAF和最优化过程可以得到时变最优核函数

其中

在时变最优核函数

得到以后，最优核时频分布如下：

$\mathrm{AOKTFR}(t,\mathrm{\ω})=\frac{1}{2\mathrm{\π}}{\∫}_{-\∞}^{+\∞}{\∫}_{-\∞}^{+\∞}\mathrm{SSTAF}(t;\mathrm{\θ},\mathrm{\τ}){\mathrm{\Φ}}_{\mathrm{opt}}(t;\mathrm{\θ},\mathrm{\τ})e^{\mathrm{j\θt}-\mathrm{j\τ\ω}}\mathrm{d\θd\τ}---\left(4\right)$

记S(x,y,t)的自适应最优核时频分布的结果为AOKTFR(x,y,t,f)，得到的高频分量为AOKTFR(x,y,t,f_H)，低频分量为AOKTFR(x,y,t,f_L)。

本发明提供了一种基于自适应最优核时频分布的地震衰减定性估计方法，首先利用自适应最优核时频分布对地震数据进行时频分解，得到地震数据的高频分量和低频分量，然后利用低频分量和高频分量的差来定性估计地震衰减。与常用的基于时频分布的地震衰减定性估计方法相比，本发明具有较高的时频分辨率，能够更为精细地刻画地震衰减，此外，该技术方案易于实现，可操作性强。

附图说明：

图1是发明基于自适应最优核时频分布的地震衰减定性估计方法的流程图示意；

图2是时频分析方法的时间-频率联合分辨率比较,其中，（a）为合成信号，(b)为Wigner-Ville分布结果，(c)为64点加窗傅立叶变换结果，(d)为128点加窗傅立叶变换结果，(e)小波变换结果，(f)自适应最优核时频分布结果；

图3是模型及合成信号示意图，其中，(a)为地质模型，(b)为正演得到的地震信号，（c）为正演地震信号的瞬时振幅；

图4是各种不同时频分析方法得到的高频分量及低频分量，其中，(a)为STFT提取的20Hz分频数据，(b)为STFT提取的30Hz分频数据，(c)为CWT_Morlet提取的20Hz分频数据，(d)为CWT_Morlet提取的30Hz分频数据，(e)为CWT_3TP提取的20Hz分频数据，(f)为CWT_3TP提取的30Hz分频数据，(g)为AOKTFR提取的20Hz分频数据，(h)为AOKTFR提取的30Hz分频数据；

图5基于各种不同时频分析方法定性估计合成信号的地震衰减，其中，(a)为基于STFT的地震衰减定性估计结果，(b)为基于CWT_Morlet的地震衰减定性估计结果，(c)为基于CWT_3TP的地震衰减定性估计结果，(d)为基于AOKTFR的地震衰减定性估计结果。

具体实施方式

下面结合附图对本发明做详细描述。

参考图1所示，本发明提供了一种基于自适应最优核时频分布的地震衰减定性估计方法，包括以下步骤：

1）采集原始地震数据

采集原始地震数据的道号、线号、以及采样点号时间，对采集到的原始地震数据进行预处理，得到叠后地震数据体，称叠后地震数据体为观测地震数据体，记为S(x,y,t)，其中x表示道号（x坐标），y表示线号(y坐标)，t表示采样点号(时间)；

2）拾取地震数据目标层的层位

利用商业软件拾取，或者直接从数据采集方获得目标层位及目标层上方的层位，记目标层位为H_T(x,y)，目标层位上面的层位为H_A(x,y)；

3）确定高频频率f_H和低频频率f_L

为了选择合适的高频f_H和低频f_L，首先，对目的层邻近的地震数据进行频谱分析，然后选择高频频率f_H和低频频率f_L，具体的，高频频率f_H和低频频率f_L应该满足三个条件：(a)由目的层邻近的地震数据得到振幅谱以后，高频f_H的幅度应与低频f_L的幅度大致相同；(b)高频f_H与低频f_L之间的频率间隔应该足够大，以保证高频成分和低频成分有足够的衰减差异；(c)高频f_H与低频f_L都应该在地震子波的频带范围之内。

4）对叠后地震数据体中每道信号计算自适应最优核时频分布并得到高频分量AOKTFR(x,y,t,f_H)和低频分量AOKTFR(x,y,t,f_L)；

5）确定修正因子α（x，y）

我们采用低频分量和高频分量的差来刻画层位H_T(x,y)附近的衰减，因此需要在H_A(x,y)层位上将高频和低频的幅度差异预先消除掉，采用如下因子来修正：

$\mathrm{\α}(x,y)=\frac{\mathrm{AOKTFR}(x,y,H_A(x,y),f_L)}{\mathrm{AOKTFR}(x,y,H_A(x,y),f_H)}---\left(1\right)$

6）估计地震衰减

利用得到的修正因子α（x，y）、高频分量AOKTFR(x,y,t,f_H)以及低频分量AOKTFR(x,y,t,f_L)估计目标层位H_T(x,y)附近的衰减SA（x,y,t）：

SA(x,y,t)=AOKTFR(x,y,t,f_L)-α(x,y)AOKTFR(x,y,t,f_H)     （2）

其中，所述对叠后地震数据体中每道信号计算自适应最优核时频分布并得到高频分量AOKTFR(x,y,t,f_H)和低频分量AOKTFR(x,y,t,f_L)具体为：

定义信号s(t)的对称短时模糊函数SSTAF(Symmetrical Short-timeAmbiguity Function)为：

$\begin{array}{c}\mathrm{SSTAF}(t;\mathrm{\θ},\mathrm{\τ})=\frac{1}{2\mathrm{\π}}{\∫}_{-\∞}^{+\∞}{s}^{*}(u-\mathrm{\τ}/2)w^{*}(u-t-\mathrm{\τ}/2)s(u+\mathrm{\τ}/2)w(u-t+\mathrm{\τ}/2)e^{-\mathrm{j\θu}}\mathrm{du}\\ =\frac{1}{2}{\∫}_{t-T+|\mathrm{\τ}/2|}^{t+T-|\mathrm{\τ}/2|}{s}^{*}(u-\mathrm{\τ}/2)s(u+\mathrm{\τ}/2)e^{-\mathrm{j\θu}}\mathrm{du}\end{array}---\left(3\right)$

其中，w是宽度为2T的对称矩形窗函数,τ为时延，θ为多普勒频率。利用SSTAF和最优化过程可以得到时变最优核函数

其中

在时变最优核函数

得到以后，最优核时频分布如下：

$\mathrm{AOKTFR}(t,\mathrm{\ω})=\frac{1}{2\mathrm{\π}}{\∫}_{-\∞}^{+\∞}{\∫}_{-\∞}^{+\∞}\mathrm{SSTAF}(t;\mathrm{\θ},\mathrm{\τ}){\mathrm{\Φ}}_{\mathrm{opt}}(t;\mathrm{\θ},\mathrm{\τ})e^{\mathrm{j\θt}-\mathrm{j\τ\ω}}\mathrm{d\θd\τ}---\left(4\right)$

记S(x,y,t)的自适应最优核时频分布的结果为AOKTFR(x,y,t,f)，得到的高频分量为AOKTFR(x,y,t,f_H)，低频分量为AOKTFR(x,y,t,f_L)。

下面为对合成信号效果的分析对比：

首先用一个典型的合成信号来说明自适应最优核时频分布的时间-频率联合分辨率及交叉项抑制能力。合成信号如图2(a)所示，该合成信号包含两个经过调制的高斯包络，两个幂率调频信号。图2(b)为Wigner-Ville分布的结果，可以看出时间-频率联合分辨率极高，但是受交叉项干扰很严重。图2(c)为64点矩形窗加窗傅立叶变换结果，可以看出时间-频率联合分辨率低，且旁瓣比较明显。图2(d)为128点高斯窗加窗傅立叶变换结果，旁瓣效应有所降低，但是时间-频率联合分辨率还是较低。图2(e)为采用Morlet小波作为母小波的一维连续小波变换结果，时间-频率联合分辨率较差，另外高频处和低频处的时间分辨率以及频率分辨率不同。图2(f)为自适应最优核时频分布的结果，虽然时间-频率联合分辨率略低于Wigner-Ville分布(如图2(b))，但是交叉项干扰抑制效果相当明显。

然后用本发明分析正演信号。含有三个地层的地质模型如图3(a)所示。第一层的厚度、纵波速度及Q值(品质因子)分别为2500m、4600m/s以及100；第三层的纵波速度及Q值分别为4600m/s以及100；第二层两端的厚度、纵波速度及Q值分别为100m、4500m/s以及80，第二层中部为储层，其厚度、纵波速度及Q值分别为100m、4200m/s以及20。所使用的震源子波为30Hz的Ricker子波，正演记录(叠后剖面)如图3(b)所示，对应的瞬时振幅如图3(c)所示。

利用加窗傅立叶变换(WFT)、以Morlet小波作为母小波的连续小波变换(CWT_Morlet)，以三参数小波作为母小波的连续小波变换(CWT_3TP)以及自适应最优核时频分布这四种时频分析方法来分析这个正演地震记录，并估计其地震衰减，比较能否反映强衰减区域。经过频谱分析，选取20Hz作为低频fL，30Hz作为高频fH。

以STFT、CWT_Morlet、CWT_3TP和AOKTFR作为时频分析工具提取的20Hz分频数据体如图4(a)、图4(c)、图4(e)及图4(g)所示。以STFT、CWT_Morlet、CWT_3TP和AOKTFR作为时频分析工具提取的30Hz分频数据体如图4(b)、图4(d)、图4(f)及图4(h)所示。可以看出，相比较于STFT、CWT_Morlet以及CWT_3TP的结果，AOKTFR(图4(g)，图4(h))均表现出很高的分辨率。

利用上节给出的流程，我们得到基于STFT、CWT_Morlet、CWT_3TP以及AOKTFR的地震衰减定性估计结果，分别如图5(a)、图5(b)、图5(c)及图5(d)所示。基于STFT的地震衰减定性估计结果，能够给出强衰减区域的指示，但是在第一层反射的位置处出现了虚假的衰减指示；基于CWT_Morlet的地震衰减定性估计结果，没有很好的指示强衰减区域，且在第一层反射的位置处给出了很强的虚假衰减指示，原因是CWT_Morlet在低频时时间分辨率低；基于CWT_3TP的地震衰减定性估计结果，相比于前两种衰减估计方法，能够较好的给出强衰减区域的指示，但在第一层反射的位置处给出了较强的虚假衰减指示。基于AOKTFT的地震衰减定性估计结果，很好的给出强衰减区域的指示，在第一层反射的位置处几乎没有虚假衰减指示。

Claims (3)

1.一种基于自适应最优核时频分布的地震衰减定性估计方法，其特征在于，包括以下步骤：

1）采集原始地震数据

采集原始地震数据的道号、线号、以及采样点号时间，对采集到的原始地震数据进行预处理，得到叠后地震数据体，记为S(x,y,t)，其中x表示道号，y表示线号，t表示采样点号时间；

2）拾取地震数据目标层的层位

拾取目标层位及目标层上方的层位，记目标层位为H_T(x,y)，目标层位上面的层位为H_A(x,y)；

3）确定高频频率f_H和低频频率f_L

首先，对目的层邻近的地震数据进行频谱分析，然后选择高频频率f_H和低频频率f_L；

4）对叠后地震数据体中每道信号计算自适应最优核时频分布并得到高频分量AOKTFR(x,y,t,f_H)和低频分量AOKTFR(x,y,t,f_L)；

5）确定修正因子α（x，y）

$\mathrm{\α}(x,y)=\frac{\mathrm{AOKTFR}(x,y,H_A(x,y),f_L)}{\mathrm{AOKTFR}(x,y,H_A(x,y),f_H)}---\left(1\right);$

6）估计地震衰减

利用得到的修正因子α（x，y）、高频分量AOKTFR(x,y,t,f_H)以及低频分量AOKTFR(x,y,t,f_L)估计目标层位H_T(x,y)附近的衰减SA（x,y,t）：

SA(x,y,t)=AOKTFR(x,y,t,f_L)-α(x,y)AOKTFR(x,y,t,f_H)    （2）。

2.根据权利要求1所述的基于自适应最优核时频分布的地震衰减定性估计方法，其特征在于，所述确定高频频率f_H和低频频率f_L具体包括：

首先，对目的层邻近地震数据进行频谱分析，得到振幅谱，然后，选取高频的幅度与低频的幅度接近、高频与低频之间的频率间隔大、并且高频与低频都在地震子波的频带范围之内的高频频率和低频频率为所需高频频率f_H和低频频率f_L。

3.根据权利要求1所述的基于自适应最优核时频分布的地震衰减定性估计方法，其特征在于，所述对叠后地震数据体中每道信号计算自适应最优核时频分布并得到高频分量AOKTFR(x,y,t,f_H)和低频分量AOKTFR(x,y,t,f_L)具体为：

定义信号s(t)的对称短时模糊函数SSTAF为：

$\mathrm{SSTAF}(t:\mathrm{\θ},\mathrm{\τ})=\frac{1}{2\mathrm{\π}}{\∫}_{-\∞}^{+\∞}{s}^{*}(u-\mathrm{\τ}/2)w^{*}(u-t-\mathrm{\τ}/2)s(u+\mathrm{\τ}/2)w(u-t+\mathrm{\τ}/2)e^{-\mathrm{j\θu}}\mathrm{du}$

$=\frac{1}{2\mathrm{\π}}{\∫}_{t-T+|\mathrm{\τ}/2|}^{t+T-|\mathrm{\τ}/2|}{s}^{*}(u-\mathrm{\τ}/2)s(u+\mathrm{\τ}/2)e^{-\mathrm{j\θu}}\mathrm{du}---\left(3\right)$

其中，w是宽度为2T的对称矩形窗函数,τ为时延，θ为多普勒频率，利用SSTAF和最优化过程可以得到时变最优核函数

其中 $r=\sqrt{{\mathrm{\τ}}^2+{\mathrm{\θ}}^2},$

在时变最优核函数得到以后，最优核时频分布如下：

$\mathrm{AOKTFR}(t,\mathrm{\ω})=\frac{1}{2\mathrm{\π}}{\∫}_{-\∞}^{+\∞}{\∫}_{-\∞}^{+\∞}\mathrm{SSTAF}(t;\mathrm{\θ},\mathrm{\τ}){\mathrm{\Φ}}_{\mathrm{opt}}(t;\mathrm{\θ},\mathrm{\τ})e^{\mathrm{j\θt}-\mathrm{j\τ\ω}}\mathrm{d\θd\τ}---\left(4\right)$

记S(x,y,t)的自适应最优核时频分布的结果为AOKTFR(x,y,t,f)，得到的高频分量为AOKTFR(x,y,t,f_H)，低频分量为AOKTFR(x,y,t,f_L)。

Images