CN104090302B - 工区地下介质频率域异常分析的方法 - Google Patents

Abstract

本发明公开了一种工区地下介质频率域异常分析的方法，工区地下介质频率域异常分析的方法，其特征在于，包括以下步骤：通过测井资料的岩石物理分析，区分工区的地下介质的砂岩地层和泥岩地层；通过时频分析算法，求取工区的地震资料在砂岩地层和泥岩地层的频率信息资料；通过地震波衰减特征提取技术，对频率信息资料进行衰减特征提取，获得地震资料在砂岩地层和泥岩地层的频率衰减特征参数；利用最小二乘法对衰减特征参数进行线性拟合或多项式拟合，得到砂岩地层的衰减趋势线和泥岩地层的衰减趋势线；以砂岩地层的衰减趋势线为主，以泥岩地层的衰减趋势线为附，分析地震资料在砂岩地层和泥岩地层的衰减趋势的差异。

Description

工区地下介质频率域异常分析的方法

技术领域

本发明涉及海上油田工程技术领域，具体是一种工区地下介质频率域异常分析的方法。

背景技术

石油和天然气作为动力能源，控制着国民经济的发展和国家战略命脉，石油天然气的自主勘探开发是维护国家稳定、国际竞争力的重要支柱。准确预测储层和储层含油气性在石油天然气开发中显得尤为重要。从上世纪开始，利用地球物理勘探方法寻找油气储层的技术已经得到发展和应用，进入本世纪，随着简单构造油气藏已经基本被勘探完毕，复杂岩性裂缝型油气藏是当下油气勘探的主要目标，随之也产生了相应的基础理论和技术手段对储层进行预测和评价，形成了一系列技术发展较为成熟、国际业内较为认可、理论基础较为先进的储层预测技术和烃类检测技术。

(1)现阶段，已发展较为成熟的地震储层预测技术主要有以下几个方面。

①地震反演技术：反演包括各种叠后、叠前反演技术，叠后反演主要预测储层厚度，叠前反演在识别储层岩性及流体方面具有重要作用。

叠后反演发展有稀疏脉冲波阻抗反演技术等。稀疏脉冲波阻抗反演过程中以精细的层位解释为基础，建立符合地质规律的波阻抗模型，应用测井约束地震反演技术，反演得到高分辨率的波阻抗数据体，对储层有利砂体的展布规律进行初步的分析。

叠前反演技术主要是基于弹性阻抗的反演和地震波全波形反演等。叠前弹性阻抗反演技术是根据地层反射振幅随炮检距(入射角)的变化规律来计算岩石弹性参数的一项技术，通过选择适当的弹性参数表达式，进行拉梅常数(λρ/μρ)交会分析，利用λρ和μρ剖面比较直观地解释含油气储层。扩展弹性阻抗反演技术是油气勘探领域正在兴起的一项新技术，利用分角度叠加数据及横波、纵波、密度等测井资料，联合反演出与岩性、含油气性相关的多种弹性参数，用于综合判别储层的物性及含油气性。地震波全波形反演技术是一种基于全波场正演模拟，从地震数据反演地球物理参数的技术方法。地震波全波形反演利用叠前地震波场的运动学和动力学信息重建地下速度结构，具有揭示复杂地质背景下构造与岩性细节信息的潜力。由于一次逆时偏移相当于全波形反演的一次迭代，因此全波形反演技术具有常规反演技术无法比拟的优越性：与地震处理技术有机融合，能够在更加准确的地层构造成像和更加真实振幅条件下精确揭示储层岩性与流体特征。

②地震属性分析技术：包括各种常规地震属性分析，地震属性在预测油气储层等方面具有独到之处。

常规地震属性是指振幅类、频率类、相位类及相关类等属性，应用的最为广泛的是亮点技术，在一定地质条件下，利用亮点技术进行储层预测是一种有效可行的方法。理论证明，地下岩层含油气之后岩石波速降低，从而引起地震反射波振幅相对增强，因此利用地震资料上的亮点进行含油气储层预测，具有一定的可行性。在综合分析钻、测井和地震资料基础上，对含油气区域的地震资料进行研究，寻找亮、暗点与含油气的关系，阐明亮、暗点形成的原因，并利用多种技术手段精细刻画亮、暗点，指出有利含油气区分布范围，为下步勘探指明方向。局部频率属性技术是在研究常规地震记录瞬时频率属性基础上提出的方法，该方法综合利用转换波数据进行交互解释，最终得到的局部频率属性反映储层流体的频变特征，具有比瞬时频率属性更高信噪比及较强的解释性，适用于天然气藏的储层预测和流体定性解释。此外还发展有多属性融合技术预测“甜点”，以此来进行储层预测。

③地震岩石物理技术，随着地球物理勘探开发从宏观向微观的不断深入研究，地震岩石物理技术可以有效地模拟储层的地震响应与传播特征，对于研究非常规油气藏地球物理特征非常重要。

④多波多分量技术，多分量地震资料包含了纵波和横波信息，比常规地震资料的信息更加丰富。岩石实验结果证明，岩石含油气后体积模量降低，而剪切模量基本不变，导致纵波速度明显降低、横波速度基本不变，从而使纵波和转换横波剖面的地震属性参数出现不同的异常。因此，利用多波地震资料联合属性分析技术进行岩性和流体识别的可靠性比单纯利用纵波要高得多，多波地震联合属性分析成为利用多波资料实现岩性、流体预测的关键技术之一。

(2)烃类检测技术目前发展有以下两大类：一是基于AVO的多种属性分析的烃类检测技术，二是基于地震波吸收和衰减两个方面，以各种时频分析为工具，发展的多项烃类检测技术。

①AVO技术是以弹性波理论为基础，利用叠前CDP道集对地震反射振幅随炮检距(或入射角)的变化特征进行研究，分析振幅随炮检距的变化规律，得到反射系数与炮检距之间的关系，并对地下反射界面上覆、下伏介质的岩性特征和物性参数做出分析，达到利用地震反射振幅信息检测油气的目的。该技术通过分析叠前地震信息随偏移距变化特征，目的就是要把AVO信息与岩性和油气联系起来，揭示AVO属性异常和烃类关系，给予AVO属性的地质含义。频率域AVO分析指通过时频转换将地震数据转换到频率域，再对目的层段地震振幅(能量)随频率以及炮检距的变化情况进行研究。由于在频率域内能够反映一些在时间域无法表征的地震信号特征，通过频率域AVO分析可以对不同频率、不同炮检距下的有效储层振幅(能量)变化进行判别，依此进行复杂气藏流体识别。与常规AVO分析不同的是，频率域AVO分析手段不但要考虑炮检距的变化情况，而且还要考虑频率域的变化特征。

②吸收和衰减主要是指地震波在经过烃类储层时，地震记录的低频段能量显著增强，而高频段能量显著减弱，通常称其为“低频共振、高频衰减”。地震信息中同时存在的低频共振和高频衰减特性是指示烃类储层的重要标志，低频共振的大小与可流动油气孔隙度有关，而高频衰减与可流动油气的渗透率有关。水层一般不会出现明显的低频共振、高频衰减特征。Klinversion烃类检测方法的基础是“低频共振、高频衰减”原理，结合钻井揭示的孔隙度、渗透率、含油气饱和度等多参数进行综合分析，精确地识别油气层、油水同层、水层或干层(含束缚水)。流体活动属性分析也是基于烃类对地震波高频能量的吸收作用，表现出的高频能量的衰减现象，其原理与Klinversion的原理基本一致，不过这种方法是通过谱分解的方式来实现不同频率地震波能量的对比。一般来讲，随着频率的增加，能量的损失逐渐加大，这种能量损失随着频率变化而变化的特征是烃类储层的典型特征。地震衰减梯度技术进行烃类检测的理论基础是地层含油气后高频信息衰减，地震波吸收系数增大，含油气储层内衰减属性与振幅属性呈正比关系。利用吸收系数进行储层含油气预测。地震微属性差异烃类检测技术，就是利用地震波的在穿过油气储层时振幅和频率微小的差异，通过分析地震波的变化，来进行差异分析，从而对岩石的性质特别是由于含有油气所造成的变化进行预测分析，进而进行烃类检测。子波分解技术正是利用子波分解重构后的地震频率、振幅及波形等属性开展烃类检测。将目的储层上方的地震反射资料的频谱减去目的储层下方的频谱得到剩余频谱，即为地震波通过目的储层的频谱衰减特征，应用频谱衰减可进行烃类检测，直接预测含油气有利区域。

(1)吸收衰减技术：地震波在地层中传播时，受到炮检点耦合、地震记录仪器的增益控制、波的球面发散、反透射的损失、偏移距的变化以及地层吸收等的影响，地震波能量是一个持续衰减的过程。地层吸收是引起地震波能量衰减的一个因素，其主要机制是岩石内部流体的相互作用及流动性。通常地层吸收系数是用来描述地层吸收大小的一个重要参数，可以用来研究目的层吸收系数的空间分布，与其它地震、测井和地质信息相结合直接用于圈定油气的分布范围，估算储量。在由固、液、气构成的多相介质中，对吸收性质影响最显著的是气态物质，在岩石孔隙饱和液中渗入少量气态物质，可以明显提高对纵波能量的吸收。相较而言，吸收系数对含气地层特别敏感，含流体地层次之。在73届(2003年)SEG年会上，Eugene Lichman发表了论文《在瞬时地震子波上检测天然气和流体的统一方法》，他结合其他人研究成果认为，在瞬时子波的幅谱上，高频部分的吸收异常往往预示着气藏的存在，而低频部分的吸收异常则常常与岩石孔隙中的流体相关。由于水、油、气本身性质存在的差异，其引起的吸收也存在着差异。当仅仅存在油、水两相介质时，这种差异要大于背景的影响，可以在先验信息下区别二者的边界。现在常用的有基于瞬时子波振幅谱分析法、地震波动力学参数法等相关方法。

Mitchell等人(1996)提出了一种计算地震信号能量衰减的分析方法——E1技术，该项技术的核心是求取信号谱的高频指数衰减系数，指数衰减函数的形式为exp(-a，w)，a为我们感兴趣的衰减系数(或吸收系数)，w为频率。计算是以一系列小视窗对地震道连续作谱分析，并计算得到相应的衰减系数，视窗大小以略大于地震波的周期为准则，连续计算使衰减系数成为时间的函数。该技术重要的思想是对背景(均匀)能量衰减的消除，因为我们只需要关心衰减的异常部分，该技术假定背景的能量衰减变化(在时间轴上)是缓慢的，消除背景后的吸收异常可以更直接地反映岩石的岩性或含油气性。

Eugene Lichman结合前人的研究成果，拓展了Mitchell等人提出的EAA分析技术，他在实验的基础上提出：弹性介质(固体和液体)和塑性介质(气体)对声波传播过程中的能量再分配的形式有着明显的不同，高吸收，特别是振幅谱高频部分的高吸收往往与岩石孔隙含气有关，低频部分的高吸收则可能预示着岩石孔隙的含水(油)饱和度。对于一个多孔介质来说，孔隙中气体对能量吸收的影响机制主要是声波在气体中的传播速度远小于在岩石骨架中的传播速度，液体部分的机制则是液体与岩石骨架之间的摩擦，液体、气体部分对总能量损失的贡献是:

ΔE_gas---气体导致的能量吸收；

ΔE_liquid---液体导致的能量吸收；

其中：E_w---总能量；

a_gas---气体的吸收系数；

a_liquid---液体的吸收系数；

w_ref---参考频率，一般为1Hz。

由上式可以看到，气体导致的能量衰减与频率之间是线性关系，随着频率的升高，其衰减就会越快，其高频段信息变化就将会更剧烈；液体引起的能量衰减与频率之间是对数关系，频率越高吸收反而变慢了。与Mitchell等提出的方法相比，EugeneLichman的研究不仅利用了瞬时子波谱中高频成分解释了气层吸收的特征，而且在实验室分析的基础上，阐明了岩石孔隙中气体、液体的能量吸收机制，着重解释了孔隙流体对瞬时子波中低频成份的影响，总结出利用低频进行含油检测的方法。

基于此种原理，我们引入一个新的比例系数θ，其定义为：

其中ΔS_h为主频右端振幅谱衰减变化，ΔS_l为主频左端振幅谱衰减变化。对于砂岩储层来说，若是θ>1，则说明含气概率更高，θ越大，则含气率可能越高；若是θ<1，则孔隙内含有流体的概率更大，当θ越小，说明地层的有效孔隙度越高。相较于吸收系数，θ对储层孔隙中流体与气体更为敏感。

当储层中存在油水二相介质时，其吸收差异不是十分明显，但是由于油水本身性质的区别，存在着(特别是低频段区域)一定的差异。油和水存在着差异：①速度差异地层含油后其速度要低于水层和干层，其对相对高频能量吸收更强；②粘滞性差异油的粘滞系数要大大高于水层的粘滞系数，其内部摩擦要大于水层，引起的能量衰减必然要大于水层；③两相介质油层大多是油水两相介质，其内部流体之极的作用要远远大于水层等单相介质地层。综合以上分析，水层和油层引起吸收必定存在着差异，在先验信息下，准确地认识到这种差异，可以有效地帮助识别油层和水层的展布范围。

(2)时频分析技术：时频分析技术已经在地球物理勘探领域受到广泛关注，并且已经应用到层位厚度预测和烃类勘探中。时频分析技术对于地震反演来讲并不是一门新技术，它被用来做时频分析已经有半个多世纪了。由于海森堡原理不确定性原理，使得我们不断提高其在时间频率的准确性。时频分析是将一维的信号变换到二维的时间频率平面上，提供长时间序列的局部信息特征。为了得到二维的时频谱，目前主要方法有：短时傅里叶变换，最大熵算法，魏格纳变换及从魏格纳变换推导出来的变换，小波变换，S变换，匹配追踪算法等等，严格来说这些算法并不是完全独立的，他们都有相互之间的联系，并不能说某一种方法的对或错，每一种方法都有自己的优势与劣势，不同的需求决定使用不同的方法。

短时傅里叶变换和最大熵算法需要用到明确的固定窗口，窗口的特性对最后时频谱的中时间域和频率域都结果有直接影响。短时傅里叶变换主要用来分析宽时窗的频谱信息，这些宽时窗包含很多的同相轴，而且频率成分主要是同相轴出现的频率。最大熵算法由于难以参数化，所以可能产生预想不到的结果。小波变换相当于用窄时窗对地震信号进行滤波，对比短时傅里叶变换，小波变换的优点是进行滤波的母小波是具有频率独立性的，不固定在某一频率；对于小波变换来说最大的缺点就是所用的小波必须是正交的，另外，选取合适的Gaussian窗，短时傅里叶变换基本上可以取得与小波变换一致的效果。S变换是小波变换和短时傅里叶变换的一种组合。S变换基于连续小波变换的思想，采用一个可移动、可变的高斯窗，因此S变换具有一些小波变换不具有的特征：与傅里叶频谱直接相关，而且与频率相关的频率分辨率。匹配追踪算法是最耗时的算法，如果选取了合适的母小波，同时也是在时间频率分辨率上最好的算法。匹配追踪的母小波只要满足容许性条件，匹配就会收敛，最重要的就是这些小波不需要满足正交性，运算的结果是得到一系列的具有时频信息的小波族。另外，还有指数匹配追踪算法，该方法通过相干计算准确的到达时间，指数匹配追踪算法的优势是算法过程中并没有用到时窗，所以也没有相应的频率干扰。下面介绍几种常见的时频分析方法。

①短时傅里叶变换

短时傅里叶变换是由Griffin在1984年提出的算法，为了显示信号的局部特征，更利于对特定时间段的信号的频率信息进行分析。

短时傅里叶的变换公式为：

公式中，STFT(ω,τ)是短时傅里叶变换后的结果，x(t)是要变换的时间信号，w(t-τ)是变换中需要的窗函数。

具体来说，短时傅里叶变换的处理方法是对信号x(t)施加一个滑动窗w(t-τ)(τ反映滑动窗的位置)后，再做傅里叶变换，窗函数w(t-τ)可以是矩形窗或者是高斯窗，因为短时傅里叶变换时频谱存在时间和频率分辨率不能同时提高的问题(海森堡不确定性)，时窗的长短和时窗移动的步长决定了时频谱的分辨率和运算效率。时窗越短时间分辨率越高，相应的频率分辨率越差，反之亦然。步长一般为时窗长度的一半。

短时傅里叶变换的出现使得信号的频谱分析功能得到质变，从单纯的频率谱变换到二维的时频平面上，进而可以分析不同时段的频率信息，不在是以整体信号进行分析，消除了其他部分的信号干扰。以后发展的各种时频分析方法只是在提高时频分辨率上。

②小波变换

小波变换(wavelet transform)是上世纪80年代后期发展起来的应用数学分支。虽然之前很多数学家物理学家对此理论做出贡献，但是理论上对小波变换进行全面而系统的表述的主要是法国数学家Y.Meyer、地质物理学家J.Morlet和理论物理学家A.Grossman的贡献。而把这一理论引入工程应用，特别是信号处理领域，法国学者I.Daubechies和S.Mallat则起着极为重要的作用。

小波变换的含义就是选取某一函数g(t)，常见有墨西哥帽小波、Hermitian小波、Morlet小波等，作为基本小波(也叫母小波)，在不同时移和伸缩变换后与待分析信号x(t)做内积：

a是尺度因子，b是时移因子，

可见小波变换具有以下功能：

一是具有多分辨率(multi-resolution)，通过不同尺度因子对原始信号的褶积，可以观察信号在长时间段的全局低频信息，也可以分析短时间段的局部高频信息；

二是在频率域可以看成是尺度a所代表的频率的G(ω)的滤波器对原始信号的滤波，由于母小波的尺度因子a与其在频率域的主频成反比，所以尺度因子越大频率越低，反之亦然；

三是通过选择合适的母小波，使得子波在时间域上为有限支撑，频率域上也比较集中，则可以使得变换后的时频谱在时间域、频率域都有表征信号局部特征的能力，有利于检测信号的瞬时特征，更加有利于地层的判断。

③S变换

S变换是由Stockwell等学者在1996年首次提出的一种可逆时频分析方法，是小波变换和短时傅里叶变换的一种组合。S变换基于连续小波变换的思想，采用一个可移动，可变的高斯窗，因此S变换具有一些小波变换不具有的特征：与傅里叶频谱直接相关，而且与频率相关的频率分辨率。

信号h(t)的一维连续S变换定义如下：

公式中：t——时间，f——频率。

S变换是具有代表性的局部时频谱，可以很简单的在时间域相加得到傅里叶变换频谱。也就是公式1-6：

其中，H(f)是h(t)的付氏变换谱，所以也就很轻松的得到S反变换计算公式1-7，可以明显看出与小波变换的差别。

S变换还可以从信号的频谱计算得到：

公式1-8的离散变换可以运用快速傅里叶变换(Fast Fourier Transform)和褶积运算，从而实现S变换的快速运算。运算结果并没有改变运算结果，从S变换得到的实部信息和虚部信息可以获得不同频率、任意时间的相位和振幅。

1.3现有技术的缺点

(1)短时傅里叶变换缺点：STFT虽然在一定程度上弥补了常规傅立叶变换不具有局部分析能力的不足，但也存在着自身不可克服的缺陷，即当窗函数确定后矩形窗的形状就确定了，相应的时间分辨率和频率分辨率也就确定了，不能随着信号频率的变化而变化，不具有自适应能力。因此我们在进行短时傅立叶变换时频分析时必须选择合理的时窗长度，太短的时窗长度会出现高频的假象，而太长的时窗长度会出现不同时刻的频率重叠现象，达不到时频分析的目的。

(2)小波变换缺点：小波变换是一种多尺度方法，对不同频率用不同尺度进行分析，它能给出比较好的时间精度，在低频区有很好的频率精确度，而在高频区频率分辨能力弱，并且小波变换有一个致命的弱点，小波变换是沿着时间方向绝对平移暗含了调制沿着包络方向传播，这样造成了相位信息仅仅是局部的，会失去其物理意义造成难以理解。

(3)S变换缺点：S变换解决了小波变换相位局部化问题，但是它实际上隐含存在一个窗口，因而在频率间断的地方频率刻画不准确，而且它也难以准确刻画振幅值，在峰值(谷值)与过零点地方振幅值相差较大。

频率异常计算的缺点：目前地层衰减异常分析主要采用大地趋势背景，没有考虑不同岩性之间的差异，因此计算出的异常存在很强的岩性影响，在油气检测中也存在着极大的不确定性，如何分岩性进行衰减异常评价，是利用频率衰减进行油气检测的关键。

发明内容

本发明的目的是克服现有技术中的不足，提供一种能够更准确地判断地层流体性质和预测储层的工区地下介质频率域异常分析的方法。

为实现上述目的，本发明所采用的技术方案是：

一种工区地下介质频率域异常分析的方法，其特征在于，包括以下步骤：A、通过测井资料的岩石物理分析，区分工区的地下介质的砂岩地层和泥岩地层；B、通过时频分析算法，求取工区的地震资料在砂岩地层和泥岩地层的频率信息资料；C、通过地震波衰减特征提取技术，对频率信息资料进行衰减特征提取，获得地震资料在砂岩地层和泥岩地层的频率衰减特征参数，并分别统计砂岩地层的频率衰减值和泥岩地层的频率衰减值；D、利用最小二乘法对砂岩地层的频率衰减值进行线性拟合或多项式拟合，得到砂岩地层的频率衰减趋势线；E、利用砂岩地层的频率衰减趋势线，分析工区地下介质频率域的异常。

所述步骤D中，利用最小二乘法对泥岩地层的频率衰减值进行线性拟合或多项式拟合，得到泥岩地层的频率衰减趋势线；所述步骤E中，以砂岩地层的频率衰减趋势线为主，以泥岩地层的频率衰减趋势线为附，分析工区地下介质频率域的异常。

所述步骤A中，应用测井约束地震反演技术，反演得到高分辨率的波阻抗数据体，对测井资料进行岩石物理分析，分析阻抗对于砂岩、泥岩区分度的敏感性，结合岩石物理分析的结果，确定合适的纵波阻抗的阈值，通过纵波阻抗的阈值，区分工区的地下介质的砂岩地层和泥岩地层。

所述时频分析算法，采用短时傅立叶变换算法、小波变换算法、S变换算法或者匹配追踪时频分析算法。

所述匹配追踪时频分析算法，具体的算法流程如下:

①用Hilbert变换计算复地震道；

②计算复地震道的瞬时包络、瞬时相位和瞬时频率；

③找到包络最大值及其对应的时间位置，得到最优小波振幅的估计值；

④计算相应时间处的瞬时相位和瞬时频率，得到最优小波相位角和主频的估计值；

⑤从原地震道中减去上面确定的小波的实部；

⑥重复①至⑤步，直到剩余的地震道幅度小于迭代误差。

所述衰减趋势线拟合的具体算法如下：

①分别统计砂岩和泥岩的衰减特征，建立衰减特征与深度的散点函数关系：

y_i＝f(x_i)i＝1,2…,m，其中，y_i是衰减值，x_i是深度，m是样本数；

②应用最小二乘法，对散点数据进行数值拟合，分别拟合出砂岩地层和泥岩地层的衰减趋势线，求取多项式使得

满足上式的多项式称为最小二乘拟合多项式，特别的当n＝1时，称为线性拟合。

本发明的有益效果：

(1)提高了利用频率域信息进行流体检测的准确度，消除了泥岩引起的频率异常。基于岩性的衰减趋势线一般是按照砂岩、泥岩进行统计分类，储层的预测、流体性质的分析都是集中在砂岩层段，泥岩引起的频率异常导致的衰减异常并没有统计到砂岩的衰减趋势线中，所以砂岩衰减趋势中的异常主要是由于地层流体引起的，是非岩性的，分析砂岩衰减趋势线中的此类异常可以更准确地判断地层流体性质和预测储层。

(2)提高了目前时频算法商业运行的速度。本发明中应用的匹配追踪算法较传统的匹配追踪算法增加了对匹配时频原子的扩张、平移、调制，增加了匹配的相关度，减少了匹配的次数，降低了运算量；同时，通过Wigner变化求取的信号能量分布和能量密度分布对比传统匹配追踪对匹配时频原子的再次时频分析大大减少了运算量，并且提高了匹配时频原子时频谱的分辨率。

附图说明

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

图1为本发明的流程图。

具体实施方式

如图1所示，工区地下介质频率域异常分析的方法，包括以下步骤：

A、通过测井资料的岩石物理分析，区分工区的地下介质的砂岩地层和泥岩地层。应用测井约束地震反演技术，反演得到高分辨率的波阻抗数据体，对测井资料进行岩石物理分析，分析阻抗对于砂岩、泥岩区分度的敏感性，结合岩石物理分析的结果，确定合适的纵波阻抗的阈值，通过纵波阻抗的阈值，区分工区的地下介质的砂岩地层和泥岩地层。

B、通过时频分析算法，求取工区的地震资料在砂岩地层和泥岩地层的频率信息资料；

本发明中的时频分析算法采用匹配追踪时频算法，通过匹配追踪时频算法中的公式的扩张、平移和调制，获取时频原子库中的最佳匹配时频原子，对最佳匹配时频原子利用Wigner分布获取新的时频能量分布，进行迭代计算，完成地震信号的匹配追踪，对获取的所有匹配时频原子获得的时频能量分布进行叠加，获取原始地震资料的高分辨率时频能量分布。

具体的匹配追踪时频算法如下：

通过研究匹配追踪分解的特性，对每个时频原子的Wigner分布求和，可以得到一个新型时频能量分布，这也是本项目研究的基础——高分辨率的时频谱，具体公式推导如下。

因为时频原子字典是完整的,匹配追踪把每个满足f(t)∈L²(R)的公式分解为

此处γ_n＝(s_n,u_n,ξ_n)且

这些时频原子被选出作为f余项的最佳匹配。

匹配追踪算法主要依赖于选择公式，该公式在每次迭代时从所有矢量中选出满足式(2-2)条件的矢量g_γn。可以通过扩张、平移和调制协变的匹配追踪定义选择公式。用和表示所选用来分别分解f⁰(t)与f′(t)的时频原子库，其中和γ_n′＝(s′_n,u′_n,ξ_n′)，存在一类选择公式如下：

当且仅当所有n≥0时

且

公式的扩张、平移和调制作为所选原子下标的简单改进。通过扩张、平移和调制得到的协方差对于进行信号分析非常重要，所有变化都被计入其中。由于加上了每一个所选原子的Wigner分布，所以从任何一个内含时频字典的f(t)的分解，可得出一个新的时频能量分布。f(t)与h(t)两个公式的相交Wigner分布由下式定义

f(t)式的Wigner分布是Wf(t,ω)＝W[f,f](t,ω)。因为Wigner分布是二次方程，可以f(t)的原子分解(2-1)式得出：

双核对应着Wigner分布的相交项。因此只保留第一个合并定义

为了定义时频平面上该能量的分布，关于时频原子的一个类似的分解算法可由Qian和Chen独立得出。

从Wigner分布的著名扩张和平移特性和时频原子的表达式(2-2)中，可知对于γ＝(s,ξ,u)有

因此

同时Wigner分布也满足

所以能量守恒公式表明

因此可将Ef(t,ω)理解为在时频平面(t,ω)上f的能量密度。与Wigner和Cohen类型分布不同，其中并不包括相交项。在当g(t)为Gaussian情况时，若W_g(t,ω)为正的话它也同时保持为正。另一方面，与某Cohen类分布相反，能量密度Ef(t,ω)不满足临界特性，这些临界特性对于信号处理的重要性无法说清。

当信号f(t)为实数时，为了用实延展系数得到一个分解结果，就必须用到实时频原子字典。对任意γ＝(s,ξ,u)，ξ≠0时，且任意相φ∈[0,2π]，有：

调整常量k(γ,φ)使||g_(γ,φ)||＝1。相φ被隐藏在复数内，现在作为实原子的参数清晰的显现出来。实时频原子的字典由D＝(g_(γ,φ))_(γ,φ)∈Γ×[0,2π]，定义Γ＝R⁺×R²。按此字典运行的匹配追踪可将任意实信号f(t)分解为：

此处的下标(γ_n,φ_n)＝(s_n,u_n,ξ_n,φ_n)被选出用于对f余项进行最佳匹配。对于任意γ＝(s,ξ,u)，实原子由下式与复原子相关

此处γ^-＝(s,-ξ,u)。f(t)的时频能量分布可从其匹配追踪分解中得到，通过加上下划线的复原子的Wigner分布有：

将式(2-9)代入此表达式中可得

此分布也满足能量密度特性。

C、通过地震波衰减特征提取技术，对频率信息资料进行衰减特征提取，获得地震资料在砂岩地层和泥岩地层的频率衰减特征参数，并分别统计砂岩地层的频率衰减值和泥岩地层的频率衰减值；

D、利用最小二乘法对衰减特征参数进行线性拟合或多项式拟合，得到砂岩地层的衰减趋势线和泥岩地层的衰减趋势线；

所述衰减趋势线拟合的具体算法如下：

①分别统计砂岩和泥岩的衰减特征，建立衰减特征与深度的散点函数关系：

y_i＝f(x_i)i＝1,2…,m，其中，y_i是衰减值，x_i是深度，m是样本数；

②应用最小二乘法，对散点数据进行数值拟合，分别拟合出砂岩地层和泥岩地层的衰减趋势线，求取多项式使得

满足上式的多项式称为最小二乘拟合多项式，特别的当n＝1时，称为线性拟合。

(a)线性拟合

所谓线性拟合，即将原有散点数据按照一定规律拟合成一条直线，即将公式2-18中的散点数据拟合成如公式2-20截距-斜率公式

P(x)＝a₀+a₁x (2-20)

通过最小二乘法，可以保证所有的样本都参与采样，可以得到显示样本变化趋势的线性最优解，样本中的异常点可以不统计入内，线性拟合参数求取算法见公式2-21，

通过上式可以得到直线公式中的斜率a₁和截距a₀，达到线性拟合的目的，进而得到衰减趋势的线性趋势线。

(b)多项式拟合

当n＞1时，即为多项式拟合，多项式拟合对比线性拟合来说，对于有曲线变化趋势的散点数据具有更好的表征作用。拟合的次数决定运算量，同时，拟合次数越高对已有的样本数据具有更好的表征，但是并不是拟合次数越高越能体现曲线的变化趋势，应该根据需求，选取合适的拟合次数进行拟合。样本中的异常数据应当剔除在拟合之外。多项式拟合参数求取算法就是，根据样本数据建立深度矩阵A、衰减特征矩阵Y，然后建立公式2-22，通过对矩阵的求解，求取多项式拟合的系数α，得到衰减趋势的多项式拟合，其中T是矩阵的转置。

A^TAα＝A^TY (2-22)

其中，

α＝[a₀ a₁ … a_n]^T

E、以砂岩地层的衰减趋势线为主，以泥岩地层的衰减趋势线为附，分析地震资料在砂岩地层和泥岩地层的衰减趋势的差异。

以上所述是本发明的优选实施方式而已，当然不能以此来限定本发明之权利范围，应当指出，对于本技术领域的普通技术人员来说，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的保护范围。

Claims (1)

1.一种工区地下介质频率域异常分析的方法，其特征在于，包括以下步骤：A、通过测井资料的岩石物理分析，区分工区的地下介质的砂岩地层和泥岩地层；B、通过时频分析算法，求取工区的地震资料在砂岩地层和泥岩地层的频率信息资料；C、通过地震波衰减特征提取技术，对频率信息资料进行衰减特征提取，获得地震资料在砂岩地层和泥岩地层的频率衰减特征参数，并分别统计砂岩地层的频率衰减值和泥岩地层的频率衰减值；D、利用最小二乘法对砂岩地层的频率衰减值进行线性拟合或多项式拟合，得到砂岩地层的频率衰减趋势线；E、利用砂岩地层的频率衰减趋势线，分析工区地下介质频率域的异常；所述步骤D中，利用最小二乘法对泥岩地层的频率衰减值进行线性拟合或多项式拟合，得到泥岩地层的频率衰减趋势线；所述步骤E中，以砂岩地层的频率衰减趋势线为主，以泥岩地层的频率衰减趋势线为附，分析工区地下介质频率域的异常；所述步骤A中，应用测井约束地震反演技术，反演得到高分辨率的波阻抗数据体，对测井资料进行岩石物理分析，分析阻抗对于砂岩、泥岩区分度的敏感性，结合岩石物理分析的结果，确定合适的纵波阻抗的阈值，通过纵波阻抗的阈值，区分工区的地下介质的砂岩地层和泥岩地层；所述时频分析算法，采用短时傅立叶变换算法、小波变换算法、S变换算法或者匹配追踪时频分析算法；所述匹配追踪时频分析算法，具体的算法流程如下:

①用Hilbert 变换计算复地震道；

②计算复地震道的瞬时包络、瞬时相位和瞬时频率；

③找到包络最大值及其对应的时间位置，得到最优小波振幅的估计值；

④计算相应时间处的瞬时相位和瞬时频率，得到最优小波相位角和主频的估计值；

⑤从原地震道中减去上面确定的小波的实部；

⑥重复①至⑤步，直到剩余的地震道幅度小于迭代误差。