CN103235339A - 一种时频分解地震流体识别方法 - Google Patents

Abstract

本发明涉及石油勘探技术领域，具体而言涉及一种时频分解地震流体识别方法，根据m(t)＝exp[-β·f²(t-τ)²]exp[i(2πf(t-τ)+φ)]的morlet小波函数建立时频原子库D，由地震道和复地震道方法计算得到morlet小波函数的初始匹配原子；针对地震道匹配分解，在初始匹配原子的邻域内以时频原子库D为约束进行迭代优选最佳的匹配原子，当达到预定的停止条件时停止匹配分解，这样就可以把原始地震道表示为一系列的morlet小波原子的线性组合；将这些最优匹配的morlet原子变换到时频域，从而可获取原始地震道的时频谱分布；在地震资料的时频谱上直接提取目的层层段的地震流体活动性属性；根据流体活动性属性预测气藏的分布范围和空间展布。本发明实施例可以准确的预测气藏的分布范围和空间展布，为天然气勘探的有利目标优选提供技术保障。

Description

一种时频分解地震流体识别方法

技术领域

本发明涉及石油勘探技术领域，具体而言，本发明涉及一种时频分解地震流体识别方法。

背景技术

随着油气勘探工作的不断深入、难度不断加大和构造圈闭的日益匮乏，有目的地寻找隐蔽油气显得越来越迫切，油气勘探开发的目标不再是以构造为主的油气藏，而是更为隐蔽的以岩性圈闭等类型为主的复杂油气藏。

近30年来，地震反演和AVO分析技术已成为储层预测和流体识别的核心技术。叠后地震波阻抗反演参数的单一性，使得它很难有效判断储层中流体的性质。基于反射振幅随偏移距变化来预测岩性、预测油气的AVO分析和叠前多参数反演技术，由于存在计算量大、稳定性差和噪音干扰等诸多问题，方法仍有待完善。利用地震资料的频率属性等相关信息来判断油气由来已久。M.A.Biot于1956年基于油气双相介质地震波传播理论建立了双相介质耗散时的地震波方程，初步分析了地震波的吸收衰减机理。Dilay等人于1995年讨论了储层内部以及储层上下方的频率谱，分析了含油气性对地震波频率的影响，引起了广泛的关注。尤其是随着时频谱分析技术的不断完善，利用地震资料的低频和高频频率信息进行流体识别技术的研究成为新的热点。但是目前在应用地震频率属性进行流体识别还存在如下问题：常规叠后地震资料在处理过程中可能对原始频率成分改造或是滤掉了低频分量，如果应用损失了对油气储层较敏感频率成分的叠后地震资料进行油气检测，就会产生虚假信息；频谱分解技术是地震资料时频属性提取和流体检测的关键，频谱分解方法比较多，并且各种方法对不同的流体检测方法适应性不同；常规的频率属性流体识别因子对油气的敏感性不尽相同，与油气的关系不明确。

在现有技术中，常规的利用地震频率属性预测油气方法是根据高、低频率段的地震波吸收衰减特征。其基本原理是当地震波在地层介质中传播时，受到波前扩散、介质吸收、界面的透射与反射、介质的各向异性、多次反射、反射界面的形态及振幅随偏移距的变化等多种因素的影响，主要表现为振幅和相位的变化。如果地震波在地层介质中传播时传播速度与频率无关，那么就不存在频散现象，地震波的衰减主要表现为振幅的变化；如果地震波速度存在频散时，地震波的衰减同时表现为振幅和相位的变化，其中地震波振幅与传播的距离和品质因子（Q）密切相关。

研究表明，如果储层岩石中含有流体（特别是油气），则储层具有低Q的特征，地震波在聚集了石油、天然气的储层中传播时，地震波会发生非弹性衰减，对高频成分的吸收衰减更强，低频能量相对增强，因此利用高、低频率段的地震波吸收衰减特征可以间接预测油气存在以及分布范围。

常规的地震频率属性是利用高、低频率段的地震波吸收衰减特征可以间接预测油气，但是实际地震信号的高频成分近似与直线段衰减，且吸收系数描述信号从主频开始衰减到结束的整个过程，而且在高频段包含了较多的噪声，严重影响求取吸收衰减属性的结果。因此在地层结构较为稳定、岩性变化不大的情况下，利用地震波高频衰减梯度因子对地层的含油气性进行检测具有比较明显的效果。但是，对于非均质性强的岩石（如火山岩）储层，其地震波场特征的复杂性，使得高频衰减梯度属性预测气层存在较强的多解性，很难实现气层与水层的正确判断，并且地震信号的高频段信噪比低，也导致利用常规的地震频率属性进行油气识别方法不稳定。

另一种现有技术中，地震时频谱分析方法是S变换和广义S变换，它们克服了短时傅里叶变换时窗固定的缺陷，具有多分辨率的能力，且计算效率快。标准S变换是由Stockwell和Mansinha等学者于1999年提出的一种新的时频分析方法，该方法是以Morlet小波为母小波的连续小波变换的延续。S变换的小波基函数是由简谐波和高斯窗函数的乘积构成，简谐波在时域仅做尺度伸缩不发生平移，而高斯窗函数则进行伸缩和平移变换。广义S变换是标准S变换的改进算法，它主要对小波基函数进行适当的调整，来达到相应分辨率，以适应实际地震资料的处理要求。

广义S变换小波基函数的形式如下：

$w_{\mathrm{\λ},p}(\mathrm{\τ}-t,f)=\frac{\mathrm{\λ}*{\left|f\right|}^p}{\sqrt{2\mathrm{\π}}}\exp (-\frac{{\mathrm{\λ}}^2{(\mathrm{\τ}-t)}^2{f}^{2p}}{2})*\exp (-i2\mathrm{\πft})---\left(1\right)$

式中，λ、p用于调节小波基函数的时宽和衰减趋势。

设信号h(t)∈L²(R)，则h(t)的广义S变换定义为，

$S_{\mathrm{\λ},p}(\mathrm{\τ},f)={\∫}_{-\∞}^{+\∞}h\left(t\right)\·w_{\mathrm{\λ},p}(\mathrm{\τ}-t,f)\mathrm{dt}---\left(2\right)$

当λ＝1,p=1时，式（2）即为标准S变换，而当p=0时，上式蜕变为短时傅里叶变换。

广义S变换能够根据信号的频率和具体形态自适应的调整窗函数的宽度，但是正是由于采用这种“加窗”的思想，导致它们的分辨率受Heisenberg不确定性准则限制，使得低频处有较高频率分辨率，高频处有较高时间分辨率，以至于产生了时间与频率分辨率“矛盾”。

发明内容

为了解决现有技术中油气的预测不准确的问题，本发明实施例提供了一种时频分解地震流体识别方法，准确的预测气藏的分布范围和空间展布，为天然气勘探的有利目标优选提供技术保障。

本发明实施例提供了一种时频分解地震流体识别方法，包括，

步骤1，根据m(t)＝exp[-β·f²(t-τ)²]exp[i(2πf(t-τ)+φ)]的morlet小波函数建立时频原子库D，其中τ为中心时间；f为主频；φ相位；β为能量衰减因子，所述β用于调节小波的时间延续度和衰减速度；所述时频原子库D表示为

为小波原子，γ_n＝{τ_n,f_n,φ_n,β_n}，τ_n、f_n、φ_n和β_n分别为第n个小波原子的中心时间、主频、相位和衰减因子；

步骤2，将输入的地震道的数据进行希尔伯特黄变换（HHT）得到复地震道，根据所述地震道和复地震道计算初始中心时间τ、主频f和相位φ，得到上述morlet小波函数的初始匹配原子，其中β在计算初始匹配原子时取值β＝4ln2；

步骤3，对所述初始匹配原子以所述时频原子库D为约束进行迭代匹配，得到最优匹配原子；

步骤4，针对输入的地震道数据分时窗进行扫描进行上述步骤2-3，得到多个最优匹配原子，多个最优匹配原子的线性组合就得到了重构的地震数据；

步骤5，根据所述重构的地震道数据，与对应的地震道数据做差，得到残差；

步骤6，当不满足预定的停止条件时，将该残差作为地震道数据重复步骤2-5，直到满足预定停止条件时停止，此时的最优匹配原子为最终匹配原子；

步骤7，将所述最终匹配原子变换到时频域，获取地震道时频谱分布；

步骤8，针对所有地震道数据进行上述步骤2-7；

步骤9，从地震资料的时频谱上直接提取目的层层段的地震流体活动性属性；

步骤10，根据所述流体活动性属性预测气藏的分布范围和空间展布。

通过上述本发明实施例的方法，应用高精度的匹配追踪算法对地震剖面进行频谱分解，有利精确匹配地震信号，提高地震资料的信噪比，可以获取丰富的地震时频属性信息，在此基础上，应用流体活动性频率属性进行储层含油气性检测，流体活动性属性效果显著，流体识别结果完全忠实于地震资料，低频信息稳定，能有效预测气层的存在与否，对于勘探初期少井或无井的地区具有推广应用价值，在气田含气预测中得到了很好的验证。

附图说明

结合以下附图阅读对实施例的详细描述，本发明的上述特征和优点，以及额外的特征和优点，将会更加清楚。

图1给出了根据本发明的一个实施例一种时频分解地震流体识别方法的流程图；

图2所示为本发明实施例一种时频分解地震流体识别方法的具体流程图；

图3a为地震资1500-2000ms频谱图；

图3b为目的层营城组2000-2500ms的频谱图；

图4a为本发明实施例w_1井旁地震道时频谱分布图；

图4b为本发明实施例w_2井旁地震道时频谱分布图；

图4c为本发明实施例w_3井旁地震道时频谱分布图；

图5所示为w_1、w_2和w_3井储层段流体活动属性的示意图。

具体实施方式

下面的描述可以使任何本领域技术人员利用本发明。具体实施例和应用中所提供的描述信息仅为示例。这里所描述的实施例的各种延伸和组合对于本领域的技术人员是显而易见的，在不脱离本发明的实质和范围的情况下，本发明定义的一般原则可以应用到其他实施例和应用中。因此，本发明不只限于所示的实施例，本发明涵盖与本文所示原理和特征相一致的最大范围。

下面的详细说明以流程图、逻辑模块和其他的符号操作表达的形式给出，可以在计算机系统上执行。一个程序、计算机执行步、逻辑块，过程等，在这里被设想为得到所希望的结果的一个或多个步骤或指令的自洽序列。这些步骤是对物理量的物理操作。这些物理量包括电、磁或者无线电信号，它们在计算机系统中被存储、传输、组合、比较以及其他操作。这些信号可是比特、数值、元素、符号、字符、条件、数字等。每个步骤都可以通过硬件、软件、固件或它们的组合执行。

图1给出了根据本发明的一个实施例一种时频分解地震流体识别方法的流程图。

包括步骤101，根据m(t)＝exp[-β·f²(t-τ)²]exp[i(2πf(t-τ)+φ)]的morlet小波函数建立时频原子库D，其中τ为中心时间；f为主频；φ相位；β为能量衰减因子，所述β用于调节小波的时间延续度和衰减速度；所述时频原子库D表示为

γ_n＝{τ_n,f_n,φ_n,β_n}，τ_n、f_n、φ_n和β_n分别为第n个小波原子的中心时间、主频、相位和衰减因子。其中，为了节省运行效率和内存空间，可以基于3个参数f、φ和β构成理论的原子库D，在实际地震道匹配分解中，初始的中心时间主要通过复地震道技术得到。

步骤102，将输入的地震道的数据进行希尔伯特黄变换（HHT）得到复地震道，根据所述地震道和复地震道计算初始中心时间τ、主频f和相位φ，得到上述morlet小波函数的初始匹配原子，其中β在计算初始匹配原子时取值β＝4ln2。

作为可以理解的，上述步骤101和步骤102不限定执行的先后。

步骤103，对所述初始匹配原子以所述时频原子库D为约束进行迭代匹配，得到最优匹配原子。

步骤104，针对输入的地震道数据分时窗进行扫描进行上述步骤102-103，得到多个最优匹配原子，多个最优匹配原子的线性组合就得到了重构的地震记录；

步骤105，根据所述重构的地震道数据，并与对应的地震道数据做差，得到残差；

步骤106，当不满足预定的停止条件时，将该残差作为地震道数据重复步骤102-105，直到满足预定停止条件时停止，此时的最优匹配原子为最终匹配原子；

步骤107，将所述最终匹配原子变换到时频域，获取地震道时频谱分布。

其中可以利用Wigner-Villa分布（WVD）方式将所述最优匹配原子变换到时频域。

步骤108，针对所有地震道数据进行上述步骤102-107；

步骤108，从地震资料的时频谱上直接提取目的层层段的地震流体活动性属性，流体活动性属性的表达式如下：

其中，M为流体活动性属性，F为流体函数，v为饱和流体多孔介质速度，ρ为饱和流体的岩石密度，κ为储层的渗透率，η为流体粘度，R为反射系数，w为反射频率。

步骤109，根据所述流体活动性属性预测气藏的分布范围和空间展布。

通过上述实施例提高计算效率，适合定量分析地震资料的频谱变化特征，并提高了流体识别的精度，利用本发明实施例中的方法可以提高对气藏预测的准确率。

如图2所示为本发明实施例一种时频分解地震流体识别方法的具体流程图。

包括步骤201，构建时频原子库。

morlet小波随频率具有衰减性，且与地震子波的波形相似，适合分析地震反射信号的频谱特征。但是标准的morlet小波函数是以固定的趋势随频率变化，因此需要对其加以改造。

首先定义小波函数：

m(t)＝exp[-β·f²(t-τ)²]exp[i(2πf(t-τ)+φ)]（3）

式中，τ为中心时间；f为主频；φ为相位；β为能量衰减因子，用于调节小波的时间延续度和衰减速度，当β＝4ln2时式（3）即为标准的morlet小波。

时频原子库的创建是基于主频、相位和能量衰减因子三参数的morlet小波，每个参数的范围足够大（如频率5～120Hz，相位-180°～180°，能量衰减因子1～10）。

所述时频原子库D可以表示为：

$D=\{m_{{\mathrm{\γ}}_n},{\mathrm{\γ}}_n\∈\mathrm{\Γ}\}$

其中，γ_n＝{τ_n,f_n,φ_n,β_n}，τ_n、f_n、φ_n和β_n分别为第n个小波原子的中心时间、主频、相位和衰减因子。

步骤202，输入地震道数据。

步骤203，将所述地震道数据进行希尔伯特黄变换（HHT）得到复地震道，并根据输入的地震道数据计算得到morlet小波函数的初始匹配原子。

首先利用希尔伯特黄变换（HHT）构建复地震道，设输入的地震道数据为s(t)，则该地震道的Hilbert变换为

$s^{\′}\left(t\right)=\frac{1}{\mathrm{\πt}}*s\left(t\right)---\left(4\right)$

则复地震道为

$y\left(t\right)=s\left(t\right)+{\mathrm{js}}^{\′}\left(t\right)=s\left(t\right)+\frac{j}{\mathrm{\πt}}*s\left(t\right)---\left(5\right)$

则瞬时振幅（或包络）、瞬时频率和瞬时相位分别为

$a\left(t\right)=\sqrt{s^2\left(t\right)+s^{\′2}\left(t\right)}---\left(6\right)$

$f\left(t\right)=\frac{1}{2\mathrm{\π}}\frac{\mathrm{d\φ}\left(t\right)}{\mathrm{dt}}---\left(7\right)$

$\mathrm{\φ}\left(t\right)={\tan }^{-1}\left(\frac{s^{\′}\left(t\right)}{s\left(t\right)}\right)---\left(8\right)$

因此复地震道y(t)的最大振幅包络a(t)所对应的时间为匹配小波原子

的中心时间τ_n，瞬时频率f(t)为

的主频f_n，瞬时相位φ(t)为

的相位φ_n，上述通过输入的地震道数据可以计算得到morlet小波函数的三个参数，并且在初始计算中β＝4ln2，这样就可以得到morlet小波函数的初始匹配原子。

步骤204，利用步骤203得到的初始匹配原子，以步骤201中的时频原子库D为约束，对地震道数据进行分时窗动态扫描得到最优匹配原子，其中分时窗是指在所述初始匹配原子的邻域内进行迭代匹配，邻域是指所述初始匹配原子附近的值。这是一个反复迭代的过程，每一次迭代提取最匹配的小波原子（n为迭代次数）。经过N次迭代后，地震数据s(t)可以表示如下：

$s\left(t\right)=\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}{a}_n{m}_{{\mathrm{\γ}}_n}+R^{\left(N\right)}f---\left(9\right)$

式中，a_n是

的振幅；R^(N)f是残差，并且R⁽⁰⁾f＝s(t)，当n=0时，所述a_n=0。

在第n次迭代中，为了快速准确求取最匹配的小波原子，采用全局预测和局部优化相结合的方法，进行分时窗动态扫描最优匹配原子。

根据式（9），第n次迭代的残差如下

$R^{\left(n\right)}f=a_n{m}_{{\mathrm{\γ}}_n}+R^{(n+1)}f$

或 $R^{(n+1)}f=R^{\left(N\right)}f-a_n{m}_{{\mathrm{\γ}}_n}---\left(10\right)$

由上式可知，在求取最优匹配原子

及其振幅a_n的过程中，需要使||R⁽ⁿ⁺¹⁾f||²→min。为此，利用如下优化方程在局部范围开窗[γ_n-Δγ_n,γ_n+Δγ_n]动态扫描确定的参数，：

$m_{{\mathrm{\&gamma;}}_n}\left(t\right)=\arg \underset{m_{{\mathrm{\&gamma;}}_n}\&Element;D}{\max }\frac{|<R^{\left(n\right)}f,m_{{\mathrm{\&gamma;}}_n}>|}{\left|\right|m_{{\mathrm{\&gamma;}}_n}\left|\right|}---\left(11\right)$

其中，γ_n＝{τ_n,f_n,φ_n,β_n}，Δγ_n为邻域，邻域Δγ_n的选择为频率间隔f_n为2；相位间隔φ_n为5°；能量衰减因子间隔β_n为1，τ_n主要是基于复地震道技术求得，在实验中发现通过在τ_n邻域里搜索的匹配原子变化不大，因此在实际操作时，省略了该参数，这样可以进一步的提高计算效率，全局预测体现在首先利用复地震道技术确定初始的匹配原子；局部优化体现在基于初始匹配原子，在其邻域范围内优选最佳的匹配原子。

局部开窗动态扫描确定最优匹配原子是一个反复迭代的过程，直到残差的能量最小就得到了最优匹配原子。

步骤205，根据所述最优匹配原子重构地震道数据，并与对应的地震道数据做差，得到残差。其中，所述重构地震记录是指，将所有的匹配原子进行线性组合从而得到重构的地震道数据。由于迭代匹配n次后，残差基本上为噪音，所以相比于原始地震道，重构地震记录的信噪比要高。

步骤206，在本实施例中采用残差能量值与预定阀值进行比较作为停止条件，判断所述残差是否大于预定的阀值（某个常数），如果大于则进入步骤207，否则进入步骤208。

步骤207，当所述残差大于预定的阀值时，将该残差作为地震道数据重复步骤203-206。

步骤208，直到残差小于预定的阀值时停止，此时的最优匹配原子为最终匹配原子。

所述时频分解就是把原始地震道表示为一系列的最优匹配原子的组合

判定时频分解完成的标准有很多，例如，可设定时频分解的迭代次数作为时频分解结束的条件，如果时频分解次数达到阀值，则时频分解停止；亦可以依据匹配时频原子的振幅参数的大小来判断，如果第n个匹配小波的振幅小于设定的阈值，则可判定时频分解完成；还可根据如上述步骤206-步骤208，以残差信号的能量值作为结束条件，假如残差能量小于某个阈值，则时频分解完成，此时的一系列最优匹配原子才是最终匹配原子，利用这些最终匹配原子来获取地震道的时频谱分布。

在得到最终匹配原子

的主频、相位和能量衰减因子后，可以利用下式计算

的振幅a_n

$a_n=\frac{|<R^{\left(n\right)}f,m_{{\mathrm{\&gamma;}}_n}>|}{{\left|\right|m_{{\mathrm{\&gamma;}}_n}\left|\right|}^2}---\left(12\right)$

步骤209，将步骤208中得到的最终匹配原子变换到时频域，在本例中可以利用最终匹配原子的Wigner-Villa分布（WVD）获取地震道时频谱分布，本发明实施例只是举出了一种变换算法，应当理解为还可以采用其它类似的变换算法，在此不再赘述。

WVD属于非线性的谱分解方法，在分析多种频率成分的信号时存在交叉项的干扰，但对于morlet小波原子而言，WVD具有良好的时频聚集性，WVD能在时间域和频率域同时保持高分辨率，精细表征morlet小波的时频谱特征。

对于时频分解得到的最终匹配原子

其WVD表示如下

$W_{m_{{\mathrm{\&gamma;}}_n}}(t,f)=\frac{1}{2\mathrm{\&pi;}}{\&Integral;}_{-\&infin;}^{\&infin;}{m}_{{\mathrm{\&gamma;}}_n}(t+\frac{\mathrm{\&tau;}}{2})\&CenterDot;m_{{\mathrm{\&gamma;}}_n}(t-\frac{\mathrm{\&tau;}}{2})\&CenterDot;\exp (-i2\mathrm{\&pi;f\&tau;})\mathrm{d\&tau;}$

$\left(13\right)$

$=\frac{1}{\sqrt{2\mathrm{\&pi;\&beta;}}{f}_n}\&CenterDot;\exp (-\frac{{2\mathrm{\&pi;}}^2}{\mathrm{\&beta;}}\frac{{(f-f_n)}^2}{f_n^2})\exp (-2\mathrm{\&beta;}{f}_n^2{(t-{\mathrm{\&tau;}}_n)}^2)$

则地震信号的地震资料时频谱分布为

$\mathrm{Af}(t,f)=\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}\frac{a_n}{\left|\right|m_{{\mathrm{\&gamma;}}_n}\left|\right|}\sqrt{W_{m_{{\mathrm{\&gamma;}}_n}}(t,f)}---\left(14\right)$

将地震信号精确分解为一组morlet小波原子的线性组合后，再通过匹配原子的WVD映射，即可以避免直接求原始信号WVD的交叉项干扰，又能得到高精度的时频谱分布。

步骤210，针对所有地震道数据进行上述步骤202-209，得到多个地震道时频谱分布。

步骤211，利用上述的高精度频谱分解方法把叠前地震道数据的目的层段分解到时间-频率联合域，得到高精度的地震道数据时频谱。在此基础上，从地震资料的时频谱上直接提取目的层层段的地震流体活动性属性，

流体活动性属性的表达式如下：

$M\&ap;F(v,\mathrm{\&rho;},\mathrm{\&kappa;},\mathrm{\&eta;})\&CenterDot;{\left(\frac{\&PartialD;R}{\&PartialD;w}\right)}^2*w---\left(15\right)$

式中F为流体函数，无量纲，M为流体活动性属性，v为饱和流体多孔介质速度，ρ为饱和流体的岩石密度，κ为储层的渗透率，η为流体粘度，R为反射系数，w为反射频率。

上式表明，低频域饱和流体多孔介质储集层中，流体的活动性与地震反射系数对反射频率的偏导的绝对值成正比，而实际地震记录为反射振幅，因此流体的活动性可以近似认为与地震反射振幅对反射频率的偏导的绝对值成正比。

式（15）的具体推导过程如下所示：

（1）饱和流体多孔介质中地震波特征

在饱和流体多孔介质的研究中通过引入渗流理论，Silin D B推导了一个弹性波方程，该方程既可以描述地震波动力学特征，又能把储层的密度、渗透率与流体的粘度等油藏参数直接与地震传播特征联系起来。Silin D B推导的弹性波方程为

$\left\{\begin{array}{c}\mathrm{\&rho;}\frac{{\&PartialD;}^{2}u}{{\&PartialD;t}^2}+{\mathrm{\&rho;}}_f\frac{\&PartialD;W}{\&PartialD;t}=\frac{1}{\mathrm{\&beta;}}\frac{{\&PartialD;}^{2}u}{{\&PartialD;x}^2}-\frac{\&PartialD;P}{\&PartialD;x}\\ W+\mathrm{\&tau;}\frac{\&PartialD;W}{\&PartialD;t}=\frac{\mathrm{\&kappa;}}{\mathrm{\&eta;}}\frac{\&PartialD;P}{\&PartialD;x}-{\mathrm{\&rho;}}_f\frac{\mathrm{\&kappa;}}{\mathrm{\&eta;}}\frac{{\&PartialD;}^{2}u}{{\&PartialD;t}^2}\\ \frac{{\&PartialD;}^{2}u}{\&PartialD;x\&PartialD;t}+{\mathrm{\&phi;\&rho;}}_f\frac{\&PartialD;P}{\&PartialD;t}=-\frac{\&PartialD;W}{\&PartialD;x}\end{array}\right.---\left(16\right)$

其中，W为饱和介质中流体的达西速度，表示单位时间内通过单位面积的流体流量；P为流体压力；u为固体骨架的位移；η和κ分别为流体粘度与储层的渗透率；β和β_f分别为固体骨架与流体的拉梅系数；τ为弛豫时间，它是孔隙空间形态、流体粘度和流体拉梅系数的函数。假设饱和流体多孔介质的孔隙度为φ，则饱和流体的岩石密度ρ可用流体密度ρ_f和岩石基质密度ρ_g表示为

ρ＝(1-φ)ρ_g+φρ_f

方程（16）两边同时除以密度ρ，可以得到一种简化形式，如下：

$\left\{\begin{array}{c}\frac{{\&PartialD;}^{2}u}{{\&PartialD;t}^2}+\frac{{\mathrm{\&rho;}}_f}{\mathrm{\&rho;}}\frac{\&PartialD;W}{\&PartialD;t}=v_b^2\frac{{\&PartialD;}^{2}u}{{\&PartialD;x}^2}-v_f^2\frac{\&PartialD;P}{\&PartialD;x}\\ {\mathrm{\&lambda;}}_f\frac{{\&PartialD;}^{2}u}{{\&PartialD;t}^2}+W+\mathrm{\&tau;}\frac{\&PartialD;W}{\&PartialD;t}=-D\frac{\&PartialD;P}{\&PartialD;x}\\ \frac{{\&PartialD;}^{2}u}{\&PartialD;x\&PartialD;t}+\frac{\&PartialD;P}{\&PartialD;t}=-\frac{\&PartialD;W}{\&PartialD;x}\end{array}\right.---\left(17\right)$

其中， $v_b^2=\frac{1}{\mathrm{\&beta;\&rho;}},$ $v_f^2=\frac{1}{{\mathrm{\&phi;\&beta;}}_f\mathrm{\&rho;}};$ 压力扩散系数 $D=\frac{\mathrm{\&kappa;}}{{\mathrm{\&phi;\&beta;}}_f\mathrm{\&eta;}};$ ${\mathrm{\&lambda;}}_f={\mathrm{\&rho;}}_f\frac{\mathrm{\&kappa;}}{\mathrm{\&eta;}}.$

以分析沿X方向传播的平面简谐纵波为例来说明流体饱和多孔介质中地震波场的特征。设固体骨架位移、流体的达西速度和压力分别如下：

$u=U_s{e}^{i(\mathrm{wt}-\mathrm{kx})},$ $\stackrel{\&OverBar;}{u}={\stackrel{\&OverBar;}{U}}_f{e}^{i(\mathrm{wt}-\mathrm{kx})},$ $P=P_0{e}^{i(\mathrm{wt}-\mathrm{kx})}---\left(18\right)$

将上式中的平面纵波和流体压力表达式代入简化方程（18）中，整理可得快纵波和慢纵波的低频域波场特征，其形式分别如下：

$\left\{\begin{array}{c}{\stackrel{\&OverBar;}{U}}_f^{\mathrm{fast}}\&ap;\frac{-\mathrm{\&epsiv;w}(1-{\mathrm{\&gamma;}}_{\mathrm{\&rho;}}{\mathrm{\&gamma;}}_v-{\mathrm{\&gamma;}}_{\mathrm{\&rho;}})}{1+{\mathrm{\&gamma;}}_{\mathrm{\&rho;}}}{U}_s^{\mathrm{fast}}+\&CenterDot;\&CenterDot;\&CenterDot;\\ {\stackrel{\&OverBar;}{U}}_f^{\mathrm{slow}}\&ap;-\mathrm{iw}(1+{\mathrm{\&gamma;}}_v)U_s^{\mathrm{slow}}+\&CenterDot;\&CenterDot;\&CenterDot;\end{array}\right.---\left(19\right)$

其中，

（一般地震频率小于1kHz，则ε的量纲小于10^-3）。

（2）饱和流体多孔介质分界面的反射系数推导

对于饱和流体多孔介质，如果平面纵波入射到介质的分界面，则地震波应该满足以下边界条件：①位移、应力连续；②由于上覆的弹性介质无渗透性，因此在边界处流体的达西速度为零。

边界条件可以写成如下形式：

$\left\{\begin{array}{c}{u}_1{|}_{X=0}=u_2{|}_{X=0}\\ -\frac{1}{{\mathrm{\&beta;}}_1}\frac{{\&PartialD;u}_1}{\&PartialD;x}{|}_{x=0}=-\frac{1}{{\mathrm{\&beta;}}_2}\frac{{\&PartialD;u}_2}{\&PartialD;x}{|}_{x=0}+\mathrm{\&phi;p}{|}_{x=0}\\ {\stackrel{\&OverBar;}{u}}_f{|}_{x=0}=0\end{array}\right.---\left(20\right)$

其中，u₁和u₂分别为M1与M2的固相位移。

为了方便起见，仅讨论垂直入射的情况，则介质M1中的垂直入射和反射的平面波可以表示为

$u_1=U_1{e}^{i(\mathrm{wt}-k_{1}x)}+{\mathrm{RU}}_1{e}^{i(\mathrm{wt}+k_{1}x)}---\left(21\right)$

介质M2中的位移和压力分别为：

$\left\{\begin{array}{c}p=\frac{1}{{\mathrm{\&phi;\&beta;}}_f}{P}_0^s{e}^{i(\mathrm{wt}-k_{s}x)}+\frac{1}{{\mathrm{\&phi;\&beta;}}_f}{P}_0^f{e}^{i(\mathrm{wt}-k_{f}x)}\\ u_2=U_2^s{e}^{i(\mathrm{wt}-k_{s}x)}+U_2^f{e}^{i(\mathrm{wt}-k_{f}x)}\end{array}\right.---\left(22\right)$

将式（19）、（21）和（22）代入边界条件（20），整理后可以得到如下方程：

$\left\{\begin{array}{c}(1+R)U_1=U_2^s+U_2^f\\ \frac{{\mathrm{ik}}_1}{{\mathrm{\&beta;}}_1}(1-R)U_1=\frac{{\mathrm{ik}}_2^s}{{\mathrm{\&beta;}}_2}{U}_2^s+\frac{{\mathrm{ik}}_2^f}{{\mathrm{\&beta;}}_2}{U}_2^f+\frac{P_0^f+P_0^s}{{\mathrm{\&beta;}}_f}\\ 0=\mathrm{iw}(1+{\mathrm{\&gamma;}}_v)U_2^s+\frac{\mathrm{\&epsiv;w}(1-{\mathrm{\&gamma;}}_{\mathrm{\&rho;}}{\mathrm{\&gamma;}}_v-{\mathrm{\&gamma;}}_{\mathrm{\&rho;}})}{1+{\mathrm{\&gamma;}}_{\mathrm{\&rho;}}}{U}_2^f\end{array}\right.---\left(23\right)$

式中R为界面的反射系数。

分析可知，在低频端（），饱和流体多孔介质分界面处的快纵波和慢纵波的位移近似成比例，并且慢纵波的位移趋近于零，说明慢纵波并不向前运动，主要在界面处产生反射，对界面的反射系数影响较大。对式（23）进行化简，可以得到低频域反射系数的近似表达式，其形式为：

$R\&ap;\frac{A_2-A_1}{A_2+A_1}+B(v_1,v_2,{\mathrm{\&rho;}}_2,\mathrm{\&kappa;},\mathrm{\&eta;})\sqrt{\frac{\mathrm{\&kappa;\&rho;}}{\mathrm{\&eta;}}w}---\left(24\right)$

式中：A1和A2分别为上下介质的波阻抗，B是饱和流体多孔介质速度、密度、流体粘度等的函数，与频率无关。

（3）基于低频域反射系数推导流体活动属性

反射系数与地震频率有关，与储层的渗透性、密度及流体粘度等油藏参数有关。令流体活动属性

代入式（24）并对其进行角频率求导，便可以得到流体活动属性的表达式，形式如下：

$M\&ap;F(v,\mathrm{\&rho;},\mathrm{\&kappa;},\mathrm{\&eta;})\&CenterDot;{\left(\frac{\&PartialD;R}{\&PartialD;w}\right)}^2*w---\left(25\right)$

式中F为流体函数，无量纲。

对研究区的地震道进行时频谱分解，分别得到了干层和储层的频谱，相比于致密干层，含油气储岩的低频能量相对增强，其频谱变化率出现正异常。因此，可以利用地震资料中含油气储集层与致密层频谱的变化率就可以获得流体活动性的变化量，进而开展储集层储集性能和地层流体变化研究。

步骤212，在利用高精度的频谱分解技术把叠前地震资料分解到时间-频率联合域以及提取地震流体活动性属性的基础上，与测井、地质等已知资料紧密结合完善流体识别方法综合，预测气藏的分布范围和空间展布，为天然气勘探的有利目标优选提供技术保障。

以下作为本发明实施例的一个方案，对本发明中的时频分解地震流体识别方法进行详细的说明和举例。

实验气田主体区位于松辽盆地南部长岭断陷中央隆起带，是达尔罕断裂上升盘形成的断裂构造与火成岩体叠合所形成的复合圈闭，发育两个局部构造高点，区内发育北东及北北西向两组断裂，其中北东向断裂规模大、延伸远。

实验地区营城组广泛发育火山岩，是深层气的主要来源之一，并且完钻井（w_1井、w_2井和w_3井）均在营城组获得工业气流，表明深层火山岩体储层具备较好的储集性能，具有一定的开发前景。由于火山岩储层埋藏深以及其强屏蔽作用，导致地震资料信噪比低、反射波连续性差以及与地层接触关系复杂多变等，致使火山岩储层及含气性的预测十分困难。

针对实验地区天然气勘探开发现状，本实施例主要从地震频率属性分析入手，将丰富的叠前地震信息与先进的预测技术结合来预测含气储层的分布范围。

地震资料的频谱分析是开展火山岩流体识别的基础工作。图3a为地震资1500-2000ms频谱图，图3b为目的层营城组2000-2500ms的频谱图，其中图3b中高频成分吸收衰减，频带变窄，主频仅28Hz左右。

图4a-图4c是w_1井、w_2井和w_3井旁地震道时频谱分布图，分析表明，1-7°、8-14°和15-21°角度叠加数据在储层处均出现明显低频强振、高频衰减特的特征，地震波能量相对移动到中低频（20-40Hz），与围岩差异明显。

地震波传播过程中流体对地震子波能量具有吸收特性，研究表明，在固、液、气构成的多相介质中，对地震波吸收性质影响最显著的是气态物质。因此，利用地震波的频率信息（流体活动属性）对火山岩储层含气性预测比较有效。

基于1-7°地震数据提取10-25Hz的流体活动属性，如图5所示为w_1、w_2和w_3井储层段流体活动属性的示意图，在该图5中均见流体活动属性的异常高值，与钻井解释的气层符合较好。在附图中的箭头是指在流体活动属性剖面上的异常高值正好与钻进曲线上解释的气层段吻合、对应。

通过上述本发明实施例的方法，应用高精度的匹配追踪算法对地震剖面进行频谱分解，有利精确匹配地震信号，提高地震资料的信噪比，可以获取丰富的地震时频属性信息，在此基础上，应用流体活动性频率属性进行储层含油气性检测，流体活动性属性效果显著，流体识别结果完全忠实于地震资料，低频信息稳定，能有效预测气层的存在与否，对于勘探初期少井或无井的地区具有推广应用价值，在气田含气预测中得到了很好的验证。

本发明可以以任何适当的形式实现，包括硬件、软件、固件或它们的任意组合。本发明可以根据情况有选择的部分实现，比如计算机软件执行于一个或多个数据处理器以及数字信号处理器。本文的每个实施例的元素和组件可以在物理上、功能上、逻辑上以任何适当的方式实现。事实上，一个功能可以在独立单元中、在一组单元中、或作为其他功能单元的一部分来实现。因此，该系统和方法既可以在独立单元中实现，也可以在物理上和功能上分布于不同的单元和处理器之间。

在相关领域中的技术人员将会认识到，本发明的实施例有许多可能的修改和组合，虽然形式略有不同，仍采用相同的基本机制和方法。为了解释的目的，前述描述参考了几个特定的实施例。然而，上述的说明性讨论不旨在穷举或限制本文所发明的精确形式。前文所示，许多修改和变化是可能的。所选和所描述的实施例，用以解释本发明的原理及其实际应用，用以使本领域技术人员能够最好地利用本发明和各个实施例的针对特定应用的修改、变形。

Claims (7)

1.一种时频分解地震流体识别方法，其特征在于包括，

步骤1，根据m(t)＝exp[-β·f²(t-τ)²]exp[i(2πf(t-τ)+φ)]的morlet小波函数建立时频原子库D，其中τ为中心时间；f为主频；φ相位；β为能量衰减因子，所述β用于调节小波的时间延续度和衰减速度；所述时频原子库D表示为

为小波原子，γ_n＝{τ_n,f_n,φ_n,β_n}，τ_n、f_n、φ_n和β_n分别为第n个小波原子的中心时间、主频、相位和衰减因子；

步骤2，将输入的地震道的数据进行希尔伯特黄变换（HHT）得到复地震道，根据所述地震道和复地震道计算初始中心时间τ、主频f和相位φ，得到上述morlet小波函数的初始匹配原子，其中β在计算初始匹配原子时取值β＝4ln2；

步骤3，对所述初始匹配原子以所述时频原子库D为约束进行迭代匹配，得到最优匹配原子；

步骤4，针对输入的地震道数据分时窗进行扫描进行上述步骤2-3，得到多个最优匹配原子，多个最优匹配原子的线性组合得到重构的地震道数据；

步骤5，根据所述重构的地震道数据，与对应的地震道数据做差，得到残差；

步骤6，当不满足预定的停止条件时，将该残差作为地震道数据重复步骤2-5，直到满足预定停止条件时停止，此时的最优匹配原子为最终匹配原子；

步骤7，将所述最终匹配原子变换到时频域，获取地震道时频谱分布；

步骤8，针对所有地震道数据进行上述步骤2-7；

步骤9，从地震资料的时频谱上直接提取目的层层段的地震流体活动性属性；

步骤10，根据所述流体活动性属性预测气藏的分布范围和空间展布。

2.根据权利要求1所述的一种时频分解地震流体识别方法，其特征在于，利用Wigner-Villa分布方式将所述最终匹配原子变换到时频域。

3.根据权利要求1所述的一种时频分解地震流体识别方法，其特征在于，在步骤5中还包括，地震数据s(t)表示如下：

其中，a_n是小波原子

的振幅；R^(N)f是残差，n为迭代次数，并且R⁽⁰⁾f＝s(t)，当n=0时，所述a_n=0。

4.根据权利要求3所述的一种时频分解地震流体识别方法，其特征在于，所述步骤5中的预定的停止条件为：

根据残差信号的能量值作为停止条件，如果残差能量小于给定的阈值，则时频分解完成。

5.根据权利要求1所述的一种时频分解地震流体识别方法，其特征在于，所述流体活动性属性的表达式如下：

$M\&ap;F(v,\mathrm{\&rho;},\mathrm{\&kappa;},\mathrm{\&eta;})\&CenterDot;{\left(\frac{\&PartialD;R}{\&PartialD;w}\right)}^2*w$

式中F为流体函数，无量纲，M为流体活动性属性，v为饱和流体多孔介质速度，ρ为饱和流体的岩石密度，κ为储层的渗透率，η为流体粘度，R为反射系数，w为反射频率。

6.根据权利要求1所述的一种时频分解地震流体识别方法，其特征在于，步骤3中所述初始匹配原子以所述时频原子库D为约束进行迭代匹配，得到最优匹配原子中进一步包括，以所述时频原子库D为约束，在所述初始匹配原子的邻域内进行迭代匹配。

7.根据权利要求6所述的一种时频分解地震流体识别方法，其特征在于，在所述初始匹配原子的邻域内进行迭代匹配进一步包括，

每一次迭代提取最匹配的小波原子

其中，n为迭代次数，经过N次迭代后，地震数据s(t)表示如下：

$s\left(t\right)=\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{a}_n{m}_{{\mathrm{\&gamma;}}_n}+R^{\left(N\right)}f---\left(9\right)$

式中，a_n是

的振幅；R^(N)f是残差，并且R⁽⁰⁾f＝s(t)，当n=0时，所述a_n=0；

在第n次迭代中，根据式（9），第n次迭代的残差如下：

$R^{\left(n\right)}f=a_n{m}_{{\mathrm{\&gamma;}}_n}+R^{(n+1)}f$

或 $R^{(n+1)}f=R^{\left(N\right)}f-a_n{m}_{{\mathrm{\&gamma;}}_n}---\left(10\right)$

利用如下优化方程（11）在局部范围开窗[γ_n-Δγ_n,γ_n+Δγ_n]动态扫描确定的参数，其中Δγ_n为γ_n的邻域，

$m_{{\mathrm{\&gamma;}}_n}\left(t\right)=\arg \underset{m_{{\mathrm{\&gamma;}}_n}\&Element;D}{\max }\frac{|<R^{\left(n\right)}f,m_{{\mathrm{\&gamma;}}_n}>|}{\left|\right|m_{{\mathrm{\&gamma;}}_n}\left|\right|}---\left(11\right)$

直到残差的能量最小，就得到了最优匹配原子。

Images