CN109917458B

CN109917458B - 油气分析方法及装置

Info

Publication number: CN109917458B
Application number: CN201910303134.3A
Authority: CN
Inventors: 田仁飞; 李晶; 胡江涛
Original assignee: Chengdu Univeristy of Technology
Current assignee: Chengdu Univeristy of Technology
Priority date: 2019-04-16
Filing date: 2019-04-16
Publication date: 2020-03-06
Anticipated expiration: 2039-04-16
Also published as: CN109917458A

Abstract

本实施例提供一种油气分析方法及装置，本实施例所述方案能够提取出地震信号中，对油气灵敏度高的特征，因此使用该方案对油气进行分析能够提高分析的精度。所述方法在获取待分析区域多道地震信号后，针对每道地震信号，对所述地震信号进行预处理获得该地震信号的多个短时功率谱；接着将所述功率谱映射到Bark域，获得所述短时功率谱对应的多段临界频带；再接着针对每段所述临界频带，计算该临界频带的临界频带谱；接着再利用等响度曲线对临界频带谱的幅值进行缩放处理，获得预加重系数；然后根据全极点模型对预加重系数进行变换获得线性预测参数；最后提取每道地震信号的所述线性预测参数的谱峰值并判断待分析区域的油气分布。

Description

油气分析方法及装置

技术领域

本申请涉及油气勘探技术领域，具体而言，涉及一种油气分析方法及装置。

背景技术

在地球物理勘探开发中，利用地震信号分析油气是常用的高效、可行的方法之一。现有技术中，在对油气进行分析时，一般是选取部分对油气响应特征灵敏的参数，然后将这些参数与钻井资料进行对比分析，从而获得油气分析的结果。然而，不同地区的地质构造、储存的油气物理性质差异较大，同一地震属性参数本身也可能属于多个地质构造，仅仅提取部分响应灵敏的参数忽略了响应微弱的信号，因此，这种方法存在分析精度低的问题。

发明内容

为了至少克服现有技术中的上述不足，本申请的目的之一在于提供一种油气分析方法，所述方法包括：

获取待分析区域多道地震信号；

针对每道地震信号，对所述地震信号进行预处理获得该地震信号的多个短时功率谱，所述多个短时功率谱是通过不同的时间窗口获得；

将所述短时功率谱映射到Bark域，获得所述短时功率谱对应的多段临界频带；

针对每段所述临界频带，计算该临界频带的临界频带谱；

利用等响度曲线对所述临界频带谱的幅值进行缩放处理，获得预加重系数；

根据全极点模型对所述预加重系数进行变换，获得线性预测参数；

提取每道地震信号的所述线性预测参数的谱峰值并判断待分析区域的油气分布。

可选地，所述对所述地震信号进行预处理获得该地震信号的多个短时功率谱的步骤包括：

去除所述地震信号中的噪声，获得第一中间信号；

对所述第一中间信号进行加窗处理，获得每个时间窗口内的第二中间信号；

针对每个所述第二中间信号进行离散傅里叶变换，获得该第二中间信号的短时功率谱。

可选地，所述将所述短时功率谱映射到Bark域的计算公式如下：

z₁(w)＝6ln{w/1200π+sqrt[(w/1200π)²+1]}

其中，w代表频率。

可选地，所述利用等响度曲线对所述临界频带谱进行加重处理，获得预加重系数的步骤为：

计算所述临界频带谱与该临界频带谱对应的等响度曲线之积，获得预加重系数；

计算等响度曲线的公式如下：

E₁(w)＝[(w²+56.8*10⁶)w⁴]/[(w²+6.3*10⁶)²*(w²+0.38*10⁹)]

其中，E₁(w)为等响度曲线。

可选地，所述根据全极点模型对所述预加重系数进行变换，获得线性预测参数的步骤前，所述方法还包括：

对所述预加重系数进行压缩处理，获得新的预加重系数；

根据全极点模型对所述新的预加重系数进行变换，获得所述线性预测参数。

本申请的另一目的在于提供一种油气分析装置，所述装置包括第一获取模块、处理模块、映射模块、计算模块、加重模块、第二获取模块和判断模块；

所述第一获取模块，用于获取待分析区域多道地震信号；

所述处理模块，用于针对每道地震信号，对所述地震信号进行预处理获得该地震信号的多个短时功率谱，所述多个短时功率谱是通过不同的时间窗口获得；

所述映射模块，用于将所述短时功率谱映射到Bark域，获得所述短时功率谱对应的多段临界频带；

所述计算模块，用于针对每段所述临界频带计算该临界频带的临界频带谱；

所述加重模块，用于利用等响度曲线对所述临界频带谱的幅值进行缩放处理，获得预加重系数；

所述第二获取模块，用于根据全极点模型对所述预加重系数进行变换，获得线性预测参数；

所述判断模块，用于提取每道地震信号的所述线性预测参数的谱峰值并判断待分析区域的油气分布。

可选地，所述处理模块具体用于：

去除所述地震信号中的噪声，获得第一中间信号；

可选地，所述映射模块用于将所述短时功率谱映射到Bark域的计算公式如下：

z₁(w)＝6ln{w/1200π+sqrt[(w/1200π)²+1]}

其中，w代表频率。

可选地，所述加重模块具体用于：

计算等响度曲线的公式如下：

E₁(w)＝[(w²+56.8*10⁶)w⁴]/[(w²+6.3*10⁶)²*(w²+0.38*10⁹)]

其中，E₁(w)为等响度曲线。

可选地，所述第二获取模块具体用于：

对所述预加重系数进行压缩处理，获得新的所述预加重系数；

相对于现有技术而言，本申请具有以下有益效果：

本申请实施例中通过将地震信号转化到Bark域上，然后再从转化后的地震信号中，提取高灵敏度的信号，通过提取出的高灵敏度信号来对油气进行分析，由于该方案能够提取出对油气反映灵敏的微弱信号的特征，因此，本实施例所述方案具有较高的分析精度。

附图说明

为了更清楚地说明本申请实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本申请的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。

图1是本申请实施例提供的油气分析设备的结构示意框图；

图2是本申请实施例提供的油气分析方法的流程示意图一；

图3是本申请实施例提供的油气分析方法的流程示意图二；

图4是本申请实施例提供的油气分析方法的流程示意图三；

图5是本申请实施例提供的过井地震剖面示意图；

图6是本申请实施例提供的地震剖面对应地震信号的线性预测参数谱峰值；

图7是本申请实施例提供的过C108、C22-3和C22井单道地震信号计算的频谱；

图8是本申请实施例提供的过C108、C22-3和C22井单道地震信号计算的线性预测参数谱；

图9是本申请实施例提供的沿目标层T2得到线性预测参数谱峰值；

图10是本申请实施例提供的沿目标层T2得到原始振幅属性；

图11是本申请实施例提供的油气分析装置的结构示意框图。

图标：100-油气分析设备；110-油气分析装置；111-第一获取模块；112-处理模块；113-映射模块；114-计算模块；115-加重模块；116-第二获取模块；117-判断模块；120-存储器；130-处理器。

具体实施方式

为使本申请实施例的目的、技术方案和优点更加清楚，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本申请一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本申请实施例的组件可以以各种不同的配置来布置和设计。

因此，以下对在附图中提供的本申请的实施例的详细描述并非旨在限制要求保护的本申请的范围，而是仅仅表示本申请的选定实施例。基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本申请保护的范围。

应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步定义和解释。

对于本领域的普通技术人员而言，可以具体情况理解上述术语在本申请中的具体含义。

在油气地球物理勘探开发中，利用地震信号检测油气是常用方法之一，该方法具有高效、可行的特点。在利用地震信号进行油气检测时，提取出能反映对油气灵敏的地震属性参数非常重要。在现有技术基于地震信号等地震资料进行油气检测的技术，主要是从地震资料中提取各种地震属性参数，然后与钻井资料进行对比，以选出对油气响应特征灵敏的地震属性参数。在这种方法中，虽然每种不同的地震属性参数有其独特的优势，但是由于各个地区储层物理性质差异大，而且地震属性本身也存在多解性，因此，现有的这种方法并不能够检测出所有地下复杂类型的油气藏。尤其是在勘探深度大和油气藏复杂度高的情况下，由于油气藏的地震信号十分微弱，常用的检测方法都很难做到获得对油气检测灵敏度高的特征。

本申请的发明人经研究发现，地震信号具有与声音信号相同的一些特点，地震信号也可以利用声学模型的波动方程进行表示。因此，也可以提取地震信号中的高灵敏度信号来进行分析。

请参见图1，图1是本申请实施例提供的油气分析设备100的结构示意框图。所述油气分析设备100包括存储器120和处理器130，所述存储器120与所述处理器130通信连接，所述存储器120中存储有可执行程序，所述处理器130用于执行所述可执行程序。

请参见图2，图2是本申请实施例提供的油气分析方法的流程示意图，所述方法可以应用于上述油气分析设备100，以解决至少一个上述问题。所述方法包括步骤S110-步骤S170。

步骤S110，获取待分析区域多道地震信号。

具体地，本实施例用于获取信息采集设备采集的待分析区域内的多道地震信号。以某油田过C108、C22-3和C22连井剖面为例子，该地震剖面共有302道地震数据，那么，本实施例便可以用于获取这302道地震数据，每道地震数据的采样数为221点，采样频率1ms。例如，可以首先提取第1道地震信号为s₁(t)。

步骤S120，针对每道地震信号，对所述地震信号进行预处理获得该地震信号的多个短时功率谱，所述多个短时功率谱是通过不同的时间窗口获得。

请参照图2，具体地，步骤S120中，对所述地震信号进行预处理获得该地震信号的功率谱可以包括步骤S121-步骤S123。

步骤S121，去除所述地震信号中的噪声，获得第一中间信号。

本实施例中，在去除所述地震信号的噪声时，可以采用的滤波器的传递函数为H(z)＝1-az^-1。其中z为经过z变换后的地震信号。也就是说，可以将经过z变换后的地震信号利用H(z)＝1-az^-1进行滤波，以去除信号中的高频噪声，获得滤波后的第一中间信号。

步骤S122，对所述第一中间信号进行加窗处理，获得每个时间窗口内的第二中间信号。

具体地，本实施例用于采用窗函数对第一中间信号进行处理，获得第一中间信号经过窗函数处理后的多个第二中间信号，也就是经过上述步骤处理后的地震信号。

例如，所述窗函数可以是汉宁窗函数，所述汉宁窗函数的时间窗口可以是20ms。

以采样点的数量为221个，时间窗口为20ms为例，则在第一个采样点前补10个频率为0的采样点，在最后一个采样点后补10个频率为0的采样点。这样，补0后共得到241个采样点。对这241个采样点进行加窗处理，使窗函数从241个采样点中的第一个采样点开始每次移动1ms，便可以获得与221个时间窗口分别对应的第二中间信号。

步骤S123，针对每个所述第二中间信号进行离散傅里叶变换，获得该第二中间信号的短时功率谱。

具体地，本实施例中，针对每个第二中间信号，首先对第二中间信号进行离散傅里叶变换，然后再根据傅里叶变换的结果计算第二中间信号的短时功率谱。

仍然以时间窗口为20ms的汉宁窗函数为例，如果经过该汉宁窗函数处理后第二中间信号为x₁(t)，那么，首先对x₁(t)进行离散傅里叶变换，得到由虚部和实部构成的变换后的信号X₁(w)。在获得经过傅里叶变换后的信号X₁(w)后，便可以计算地震信号在一个时间窗口内的功率谱，即第二中间信号的短时功率谱。其中，计算短时功率谱P₁(w)的公式为：

P₁(w)＝[ReX₁(w)]²+[Im(X₁(w))]²

本实施例中，t为接收到信号的时间，w为地震信号中信号的频率。

请继续参照图2，步骤S130，将所述功率谱映射到Bark域，获得所述短时功率谱对应的多段临界频带。

本实施例中的Bark域即声音的心理声学尺度。

本实施例用于将功率谱映射到bark频率基z函数上，从而获得多段临界频带。其中，将所述短时功率谱映射到Bark域的计算公式为：

z₁(w)＝6ln{w/1200π+sqrt[(w/1200π)²+1]}

步骤S140，针对每个所述临界频带，计算该临界频带所对应的临界频带谱。

在计算出临界频带后，对于每个时间窗口，分别计算出该时间窗口内的采样点信号映射到Bark域上的每个临界频带对应的临界频带谱系数。

步骤S140的具体步骤为将各个时间窗口对应的临界频带乘以一定的比重后进行叠加，获得该临界频带对应的临界频带谱。

具体地，本实施例用于在计算出临界频带谱系数后，对于不同的时间窗口对应的临界频带中各个频率的幅值乘以临界频带谱系数后的值求和，这样便可以得到每个临界频带对应的临界频带谱z₁(w_i)。

其中，i表示一个时间窗口内，临界频带的序号。

步骤S150，利用等响度曲线对所述临界频带谱的幅值进行缩放处理，获得预加重系数。

计算所述临界频带谱与该临界频带谱对应的等响度曲线之积，获得预加重系数。

本实施例用于在计算出各个临界频带谱后，采用多个等响度曲线分别对各个临界频带谱进行预加重处理。

例如，等响度曲线的表达式为E₁(w_i)，那么，将等响度曲线与所述临界频带谱进行相乘，这样便可以获得预加重系数D₁(w_i)。也就是说，预加重系数的计算公式为：

D₁(w_i)＝E₁(w_i)·z₁(w_i)

其中，等响度曲线的计算方式为：

E₁(w)＝[(w²+56.8*10⁶)w⁴]/[(w²+6.3*10⁶)²*(w²+0.38*10⁹)]

本实施例中，所述预加重系数可以反映地震信号对不同频率的不同敏感性，因此根据该预加重系数可以识别油气信号的特征。

步骤S160，根据全极点模型对所述预加重系数进行变换，获得线性预测参数。

可选地，本实施例中，步骤S160包括步骤S161和步骤S162。

请参照图4，步骤S161，对所述预加重系数进行压缩处理，获得新的预加重系数。

具体地，对所述预加重系数进行压缩处理以使所述预加重系数按照预设条件缩放，获得新的预加重系数。

例如，对预加重系数进行压缩处理可以采用公式：

DF₁(w)＝D₁(w)^0.33

其中，DF₁(w)代表压缩后的预加重系数。

步骤S162，根据全极点模型对所述新的预加重系数进行变换，获得所述线性预测参数。

可选地，步骤S162的具体过程为：对压缩后的预加重系数，即中间加重系数进行z变换，然后对z变换后的新的预加重系数采用全极点模型进行变换。其中，采用全极点模型进行变换的公式为：

其中，G为增益，a_i为数字滤波器系数，Z为Z变换，p为预测阶数。

请继续参照图2，步骤S170，提取每道地震信号的所述线性预测参数的谱峰值并判断待分析区域的油气分布。

本实施例用于在计算出线性预测参数后，对所述线性预测参数进行傅里叶变换，从而获得所述线性预测参数对应的频谱，这样可以获得各个线性预测参数的频谱值。

仍然以上述302道地震数据为例，地震剖面的示意图如图5所示，其中，xLine为平面上沿x轴的测线号，InLine为平面上沿y轴的测线号，Two-WayTime代表接收到信号的时间。x轴、y轴分别为平面直角坐标系的两个轴。其中，在本实施例的附图中，相同的符号代表相同的含义。在进行油气检测的过程中，首先对通过步骤S110获取的各道地震信号进行步骤S120-步骤S160的过程处理，获得各道地震信号的线性预测参数。在对各道地震信号进行上述处理后，对各道地震经过步骤S120-步骤S160处理之后得到的线性预测参数分别进行谱峰值提取，各道地震信号对应的谱峰值的图如图6所示。将图6与图5对比知可知感知线性预测参数谱峰值异常区域主要集中在T2层(图中虚线框所在区域)附近，再对比井所在位置C108、C22-3和C22目标层的感知线性预测参数谱峰值变化情况，可知：C108井线性预测参数谱峰值、C22-3井的线性预测参数谱峰值较小，而C22井线性预测参数谱峰值最大。查阅实际钻探井结果，C22井在目标层T2最高日产油46.2吨，C22-3日产水30方、无油，C108钻井全井无油气显示，测井解释在T2层为干层。以上很好的说明了线性预测参数谱峰值能够很好的刻画储层含油气信息。

为了与常规的频谱参数进行对比，我们在一道地震数据中，提取了C108、C22-3和C22目标层的频谱如图7所示，其中，图7中横坐标代表样点数，纵坐标代表振幅。过C108、C22-3和C22井沿T开20ms时窗计算得到频谱图，然后提取该频谱的最大值，即可得到该地震道的线性预测参数的频谱图，请参见图8。从图7和图8中可见，频谱及其频谱峰值有一定差异，同时也可明显的看出干井C108井频谱也存在异常，且三口井频谱峰值差异不太大；而图8中，干井C108线性预测参数的频谱几乎无异常，而产油井异常特征更为明显，特别是峰值更大。进一步，我们提取频谱峰值和线性预测参数谱峰值，并计算油层/干层，油层/水层和水层/干层的频谱峰值比，其结果见表1。

表1

项目	油井/干井	油井/水井	水井/干井
				SP比	2.6	1.4	1.9
PLPSP比	29.5	2.7	11.0

从表1中，可以明显的看出油层/干层线性预测参数谱峰值比为29.5；而油层/干层谱峰值比仅2.6。即使油层/水层线性预测参数谱峰值比也有2.6。所以线性预测参数更能够有效的区别油层，也能够区分水层，其灵敏度比常规的频谱属性具有更好的识别效率。为了刻画目标层横向的油气分布，针对整个检测区地震数据，重复步骤S120—步骤S160，可计算所有地震信号的线性预测参数，然后便可以根据各个线性预测参数计算对应地震信号的感知线性预测参数谱峰值，沿T层提取获得感知线性预测参数谱峰值如图9。上述区域的原始地震振幅属性如图10。从图9中，也能够明显看出C108、C22-3和C22附近线性预测参数谱峰值特征具有明显差异，含油层具有更大的线性预测参数谱峰值(图中黑色异常所示)。相比原始地震振幅属性如图5中，则较难区分C22油层和C22-3水层。因此，图9中黑色异常区分布范围可表示T2油层分布特征和范围。

请参见图11，本申请的另一目的在于提供一种油气分析装置110，所述装置包括第一获取模块111、处理模块112、映射模块113、计算模块114、加重模块115、第二获取模块116和判断模块117。所述油气分析装置110包括一个可以软件或固件的形式存储与所述存储器120中或者固化在所述油气分析设备100的操作系统(operating system，OS)中的软件功能模块。

所述第一获取模块111，用于获取待分析区域多道地震信号。

本实施例中，所述第一获取模块111用于执行图2所示的步骤S110，关于第一获取模块111的具体描述可参照对所述步骤S110的描述。

所述处理模块112，用于针对每道地震信号，对所述地震信号进行预处理获得该地震信号的多个短时功率谱，所述多个短时功率谱是通过不同的时间窗口获得。

可选地，对所述地震信号进行预处理获得该地震信号的功率谱的包括步骤S121-步骤S123。

步骤S121，去除所述地震信号的噪声，获得第一中间信号。

步骤S123，对所述第二中间信号进行离散傅里叶变换，获得所述地震信号的短时功率谱。

本实施例中，所述处理模块112用于执行图2所示的步骤S120，关于处理模块112的具体描述可参照对所述步骤S120的描述。

所述映射模块113，用于将所述短时功率谱映射到Bark域，获得所述短时功率谱对应的多段临界频带。

本实施例中，所述映射模块113用于执行图2所示的步骤S130，关于映射模块113的具体描述可参照对所述步骤S130的描述。

所述计算模块114，用于针对每个所述临界频带计算该临界频带的临界频带谱。

本实施例中，所述计算模块114用于执行图2所示的步骤S140，关于计算模块114的具体描述可参照对所述步骤S140的描述。

所述加重模块115，用于利用等响度曲线对所述临界频带谱的幅值进行缩放处理，获得预加重系数。

具体地，所述计算模块114可以通过计算所述临界频带谱与等响度曲线之积，获得预加重系数。

本实施例中，所述加重模块115用于执行图2所示的步骤S150，关于加重模块115的具体描述可参照对所述步骤S150的描述。

所述第二获取模块116，用于根据全极点模型对所述预加重系数进行变换，获得线性预测参数

具体地，第二获取模块116具体用于执行步骤S161和步骤S162。

步骤S161，对所述预加重系数进行压缩处理以使所述预加重系数按照预设条件减小，获得中间加重系数。

步骤S162，根据全极点模型对所述中间加重系数转换，获得所述线性预测参数。

本实施例中，所述第二获取模块116用于执行图2所示的步骤S160，关于第二获取模块116的具体描述可参照对所述步骤S160的描述。

所述判断模块117，用于提取所述线性预测参数的谱峰值判断待分析区域的油气分布。

本实施例中，所述判断模块117用于执行图2所示的步骤S70，关于判断模块117的具体描述可参照对所述步骤S170的描述。

综上所述，本申请实施例中通过将地震信号转化到Bark域上，然后再从转化后的地震信号中，提取高灵敏度的信号，通过提取出的高灵敏度信号来对油气进行分析，由于该本实施所述方案能够提取出对油气反映灵敏的微弱信号的特征，因此，本实施例所述方案具有较高的分析精度。

本实施例中，以上所述，仅为本申请的各种实施方式，但本申请的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本申请揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本申请的保护范围之内。因此，本申请的保护范围应所述以权利要求的保护范围为准。

Claims

1.一种油气分析方法，其特征在于，所述方法包括：

获取待分析区域多道地震信号；

针对每段所述临界频带，计算该临界频带的临界频带谱；

2.根据权利要求1所述的方法，其特征在于，所述对所述地震信号进行预处理获得该地震信号的多个短时功率谱的步骤包括：

去除所述地震信号中的噪声，获得第一中间信号；

3.根据权利要求1所述的方法，其特征在于，所述将所述短时功率谱映射到Bark域的计算公式如下：

z₁(w)＝6ln{w/1200π+sqrt[(w/1200π)²+1]}

其中，w代表频率。

4.根据权利要求1所述的方法，其特征在于，所述利用等响度曲线对所述临界频带谱进行加重处理，获得预加重系数的步骤为：

计算等响度曲线的公式如下：

E₁(w)＝[(w²+56.8*10⁶)w⁴]/[(w²+6.3*10⁶)²*(w²+0.38*10⁹)]

其中，E₁(w)为等响度曲线，w代表频率。

5.根据权利要求1所述的方法，其特征在于，所述根据全极点模型对所述预加重系数进行变换，获得线性预测参数的步骤前，所述方法还包括：

对所述预加重系数进行压缩处理，获得新的预加重系数；

6.一种油气分析装置，其特征在于，所述装置包括第一获取模块、处理模块、映射模块、计算模块、加重模块、第二获取模块和判断模块；

所述第一获取模块，用于获取待分析区域多道地震信号；

7.根据权利要求6所述的装置，其特征在于，所述处理模块具体用于：

去除所述地震信号中的噪声，获得第一中间信号；

8.根据权利要求6所述的装置，其特征在于，所述映射模块用于将所述短时功率谱映射到Bark域的计算公式如下：

z₁(w)＝6ln{w/1200π+sqrt[(w/1200π)²+1]}

其中，w代表频率。

9.根据权利要求6所述的装置，其特征在于，所述加重模块具体用于：

计算等响度曲线的公式如下：

E₁(w)＝[(w²+56.8*10⁶)w⁴]/[(w²+6.3*10⁶)²*(w²+0.38*10⁹)]

其中，E₁(w)为等响度曲线，w代表频率。

10.根据权利要求6所述的装置，其特征在于，所述第二获取模块具体用于：