CN110794458B - 基于时频分析的含气性检测方法、装置及存储介质 - Google Patents
基于时频分析的含气性检测方法、装置及存储介质 Download PDFInfo
- Publication number
- CN110794458B CN110794458B CN201911042135.3A CN201911042135A CN110794458B CN 110794458 B CN110794458 B CN 110794458B CN 201911042135 A CN201911042135 A CN 201911042135A CN 110794458 B CN110794458 B CN 110794458B
- Authority
- CN
- China
- Prior art keywords
- time
- seismic
- frequency
- wavelet
- representing
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000004458 analytical method Methods 0.000 title claims abstract description 23
- 238000001514 detection method Methods 0.000 title claims abstract description 16
- 238000003860 storage Methods 0.000 title claims abstract description 16
- 238000001228 spectrum Methods 0.000 claims abstract description 62
- 238000000034 method Methods 0.000 claims abstract description 50
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 27
- 239000013598 vector Substances 0.000 claims description 22
- 239000011159 matrix material Substances 0.000 claims description 16
- 238000010276 construction Methods 0.000 claims description 8
- 238000004590 computer program Methods 0.000 claims description 5
- 238000013461 design Methods 0.000 claims description 5
- 238000011161 development Methods 0.000 abstract description 5
- 239000007789 gas Substances 0.000 description 40
- 230000006870 function Effects 0.000 description 36
- 230000008569 process Effects 0.000 description 11
- 230000006872 improvement Effects 0.000 description 9
- 238000000354 decomposition reaction Methods 0.000 description 7
- 238000005516 engineering process Methods 0.000 description 7
- 238000010586 diagram Methods 0.000 description 6
- 238000012545 processing Methods 0.000 description 5
- 230000003595 spectral effect Effects 0.000 description 5
- 238000013519 translation Methods 0.000 description 4
- 230000008859 change Effects 0.000 description 3
- 235000021185 dessert Nutrition 0.000 description 3
- 238000005553 drilling Methods 0.000 description 3
- 230000002123 temporal effect Effects 0.000 description 3
- 230000002238 attenuated effect Effects 0.000 description 2
- 235000009508 confectionery Nutrition 0.000 description 2
- 238000009826 distribution Methods 0.000 description 2
- 238000011156 evaluation Methods 0.000 description 2
- 238000004519 manufacturing process Methods 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- VNWKTOKETHGBQD-UHFFFAOYSA-N methane Chemical compound C VNWKTOKETHGBQD-UHFFFAOYSA-N 0.000 description 2
- 230000000750 progressive effect Effects 0.000 description 2
- 238000011160 research Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 230000009466 transformation Effects 0.000 description 2
- OAICVXFJPJFONN-UHFFFAOYSA-N Phosphorus Chemical compound [P] OAICVXFJPJFONN-UHFFFAOYSA-N 0.000 description 1
- 230000005856 abnormality Effects 0.000 description 1
- 238000010521 absorption reaction Methods 0.000 description 1
- 238000013459 approach Methods 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000015572 biosynthetic process Effects 0.000 description 1
- 238000004364 calculation method Methods 0.000 description 1
- 230000001413 cellular effect Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 238000004891 communication Methods 0.000 description 1
- 229910000514 dolomite Inorganic materials 0.000 description 1
- 239000010459 dolomite Substances 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000005284 excitation Effects 0.000 description 1
- 239000012530 fluid Substances 0.000 description 1
- 238000005755 formation reaction Methods 0.000 description 1
- 238000010230 functional analysis Methods 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 230000004807 localization Effects 0.000 description 1
- 230000014759 maintenance of location Effects 0.000 description 1
- 230000005055 memory storage Effects 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 239000003345 natural gas Substances 0.000 description 1
- 230000003287 optical effect Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 239000003209 petroleum derivative Substances 0.000 description 1
- 238000010183 spectrum analysis Methods 0.000 description 1
- 238000006467 substitution reaction Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/282—Application of seismic models, synthetic seismograms
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/62—Physical property of subsurface
- G01V2210/624—Reservoir parameters
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
本说明书实施例提供一种基于时频分析的含气性检测方法、装置及存储介质。所述方法包括：接收目标区域的地震信号；根据地震信号，通过褶积模型构建反演方程；将Lp范数作为稀疏约束项，根据所述反演方程建立目标函数；使用快速迭代软阈值算法求解目标函数，得到目标区域的时频谱；根据所述时频谱检测目标区域的含气性。本说明书实施例提供的含气性检测方法可以获得高分辨、聚焦性好的时频谱属性，能够更为准确的进行地下储层含气性的检测，降低储层勘探和开发的风险。
Description
技术领域
本说明书实施例涉及石油、天然气地震勘探反演和解释领域，特别涉及一种基于时频分析的含气性检测方法、装置及存储介质。
背景技术
随着油气勘探对象逐渐面向构造复杂、介质复杂、地表复杂、深层的非常规油气藏，勘探难度越来越大，取得的油气勘探“亮点”越来越少，在深层和非常规油气藏中寻找油气已经成为现阶段石油工业保持稳产或提高产量的主要任务。此外，现有油气勘探逐渐以薄层探测为主，在这种情况下，常规的谱分解方法已经越来越不能满足地震解释对精度和分辨率的要求，在薄层探测中甚至会得出错误的结论。为降低勘探风险，对勘探区做出更准确、更精细的解释和评价，急需形成一套非常规深层储层含气性预测技术。
传统的时频分析技术有短时傅里叶变换(STFT)、小波变换(CWT)、S变换和广义S变换等，它们已经被广泛地应用到地震勘探的各个方面，但目前存在的这些谱分解方法都会受到海森堡测不准原理的限制，时间分辨率和频率分辨率不能同时提高，很难进行高精度的储层和层位识别，具有一定的局限性，因此有必要发明一种分辨率更高、精确度更高的含气性检测技术。
发明内容
本说明书实施例的目的是提供一种基于时频分析的含气性检测方法、装置及存储介质，以更加准确地进行地下储层含气性的检测，降低储层勘探和开发的风险。
为解决上述问题，本说明书实施例提供一种基于时频分析的含气性检测方法、装置及存储介质是这样实现的。
一种基于时频分析的含气性检测方法，所述方法包括：接收目标区域的地震信号；根据所述地震信号，通过褶积模型构建反演方程；将Lp范数作为稀疏约束项，根据所述反演方程建立目标函数；使用快速迭代软阈值算法求解目标函数，得到目标区域的时频谱；根据所述时频谱检测目标区域的含气性。
一种基于时频分析的含气性检测装置，所述装置包括：接收模块，用于接收目标区域的地震信号；构建模块，用于根据所述地震信号，通过褶积模型构建反演方程；建立模块，用于将Lp范数作为稀疏约束项，根据所述反演方程建立目标函数；求解模块，用于使用快速迭代软阈值算法求解目标函数，得到目标区域的时频谱；检测模块，用于根据所述时频谱检测目标区域的含气性。
一种计算机可读存储介质，其上存储有计算机程序指令，所述计算机程序指令被执行时实现：接收目标区域的地震信号；根据所述地震信号，通过褶积模型构建反演方程；将Lp范数作为稀疏约束项，根据所述反演方程建立目标函数；使用快速迭代软阈值算法求解目标函数，得到目标区域的时频谱；根据所述时频谱检测目标区域的含气性。
由以上本说明书实施例提供的技术方案可见，本说明书实施例提供了一种基于时频分析的含气性检测方法。本说明书实施例可以利用多维褶积模型构建反演方程，使用Lp准则约束条件产生稀疏时频分布，并且利用快速迭代软阈值算法(FISTA)用来求解Lp范数正则化问题，获得高分辨、聚焦性好的时频谱属性，进而通过不同频率谱成分的对比检测低频阴影，能够消除强振幅的影响，从而能够更为准确的进行地下储层含气性的检测，降低储层勘探和开发的风险。
附图说明
为了更清楚地说明本说明书实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本说明书中记载的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本说明书实施例一种基于时频分析的含气性检测方法的流程图；
图2为本说明书实施例含气性检测所在的研究工区图；
图3为本说明书实施例研究工区连井地震剖面图；
图4a为本说明书实施例35Hz频率的幅谱剖面图；
图4b为本说明书实施例30Hz频率的幅谱剖面图；
图4c为本说明书实施例21Hz频率的幅谱剖面图；
图5a为本说明书实施例上段储层原始地震振幅属性图；
图5b为本说明书实施例上段储层原始地震甜点属性图；
图5c为本说明书实施例上段储层30Hz频率成分的振幅谱属性图；
图5d为本说明书实施例上段储层21Hz频率成分的振幅谱属性图；
图6a为本说明书实施例下段储层原始地震振幅属性图；
图6b为本说明书实施例下段储层原始地震甜点属性图；
图6c为本说明书实施例下段储层35Hz频率成分的振幅谱属性图；
图6d为本说明书实施例下段储层30Hz频率成分的振幅谱属性图；
图7为本说明书实施例一种基于时频分析的含气性检测装置的功能模块图。
具体实施方式
下面将结合本说明书实施例中的附图，对本说明书实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本说明书一部分实施例，而不是全部的实施例。基于本说明书中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都应当属于本说明书保护的范围。
传统的时频分析技术有短时傅里叶变换(STFT)、小波变换(CWT)、S变换和广义S变换等，它们已经被广泛地应用到地震勘探的各个方面，但目前存在的这些谱分解方法都会受到海森堡测不准原理的限制，时间分辨率和频率分辨率不能同时提高，很难进行高精度的储层和层位识别，具有一定的局限性，因此有必要发明一种分辨率更高、精确度更高的含气性检测技术。
图1为本说明书实施例一种基于时频分析的含气性检测方法的流程图。在本说明书实施例中，执行所述设备运行状态监测方法的主体可以是具有逻辑运算功能的电子设备，所述电子设备可以是服务器或客户端。所述客户端可以为台式电脑、平板电脑、笔记本电脑、工作站等。当然，客户端并不限于上述具有一定实体的电子设备，其还可以为运行于上述电子设备中的软体，或者，还可以是一种通过程序开发形成的程序软件。该程序软件可以运行于上述电子设备中。
如图1所示，所述基于时频分析的含气性检测方法可以包括以下步骤。
S110：接收目标区域的地震信号。
地震信号采集是油气地震勘探工程中第一道工序，也是最为重要的工序，在这道工序里可以采用地震信号接收和记录系统进行地震信号采集。在一些实施例中，可以把传感地震信号的装置称为地震检波器，把采集和记录地震信号的装置称为地震勘探仪器，它是地震仪器是关键装备。地震检波器和地震勘探仪器可以联合工作来实现完整的地震数据采集功能。
在一些实施例中，服务器可以采用任何方式接收目标区域的地震信号。例如，地震数据采集系统可以直接向服务器发送目标区域的地震信号，服务器可以进行接收；又如除去所述服务器以外的其它电子设备可以向服务器发送目标区域的地震信号，服务器可以进行接收，在本说明书实施例中，对服务器采用何中方式接收目标区域的地震信号不作限定。
S120：根据所述地震信号，通过褶积模型构建反演方程。
在一些实施例中，可以通过以下方法步骤来构建反演方程。
S121：根据所述地震信号，构建由不同频率的地震子波组成的子波库；。
在本说明书实施例中，所述地震子波是一段具有确定的起始时间、能量有限且有一定延续长度的信号，它是地震记录中的基本单元。一般认为，地震震源激发时所产生的地震波仅是一个延续时间极短的尖脉冲，随着尖脉冲在粘弹性介质中传播，尖脉冲的高频成分很快衰减，波形随之增长，便形成了地震子波，一个地震子波一般有2至3个相位的延续长度，大约有90ms左右，然后以地震子波的形式在地下传播。
在一些实施例中，可以根据接收到的地震信号，通过频谱分析，构建不同频率的地震子波组成的子波库。
在一些实施例中，由于不同种类的地震子波具有不同的特点，在反演方程构建过程中，使用不同种类的地震子波，所得到的结果也不同。例如，Ricker子波具有形状简单、收敛块等特点，在本说明书实施例中，可以构建由不同频率的零相位的Ricker子波组成的子波库。在一些实施例中，还可以构建由Morlet小波或Haar小波等地震子波组成的子波库。
S122：根据所述地震信号和子波库，通过褶积模型构建反演方程。
在地震勘探信号中，高频信号部分的时间分辨率高，低频信号的频率分辨率高，如果在处理地震信号高频部分的时候牺牲部分时间分辨率来弥补频率分辨率，在低频的时候牺牲部分频率分辨率来弥补时间分辨率，则可以得到一种综合分辨率较高的时频谱。这种方式可以通过小波变换来实现。
所述小波变换(wavelet transform，WT)是一种新的变换分析方法，它继承和发展了短时傅立叶变换局部化的思想，同时又克服了窗口大小不随频率变化等缺点，能够提供一个随频率改变的“时间-频率”窗口，是进行信号时频分析和处理的理想工具。它的主要特点是通过变换能够充分突出问题某些方面的特征，能对时间(空间)频率的局部化分析，通过伸缩平移运算对信号(函数)逐步进行多尺度细化，最终达到高频处时间细分，低频处频率细分，能自动适应时频信号分析的要求，从而可聚焦到信号的任意细节，从而解决了Fourier变换的困难问题。
在本说明书实施例中，小波逆变换是在小波变换的基础上进行的逆变换。一段地震信号s(t)的小波逆变换通常可以表示为：
其中，t表示时间，a为尺度因子，尺度因子a与频率有关，a>1时子波被拉伸，对应较低的频率，a<1时子波被压缩，对应较高的频率；τ为平移因子，C(τ,a)为小波变换的谱或系数，ψ((t-τ)/a)为母小波通过拉伸和平移构成的子波库。在地震勘探中，母小波通常可以取为Morlet小波、Ricker子波或Haar小波。
对公式(1)中的平移因子的积分用褶积模型表达，得到如下表达式：
其中，*表示褶积运算，a-5/2ψ(t/a)为地震子波组成的子波库或基函数。
对公式(2)中尺度因子的积分离散表示，得到由地震信号和子波库，通过褶积模型构建的反演方程：
其中，fg表示频率，w(t,fg)对应a-5/2ψ(t/a)表示地震子波组成的子波库或基函数，r(t,fg)对应C(t,a)表示某一频率或尺度下的谱或系数，其中，N＝1,2,3...。
S130：将Lp范数作为稀疏约束项，根据所述反演方程建立目标函数。
范数(norm)是数学中的一种基本概念。在泛函分析中，它定义在赋范线性空间中，并满足一定的条件，即非负性、齐次性、三角不等式。在数学上，范数通常可以包括向量范数和矩阵范数，向量范数表征向量空间中向量的大小，矩阵范数表征矩阵引起变化的大小。一种非严密的解释就是，对应向量范数，向量空间中的向量都是有大小的，这个大小如何度量，就是用范数来度量的，不同的范数都可以来度量这个大小，就好比米和尺都可以来度量远近一样；对于矩阵范数，通过运算，可以将向量X变化为B，矩阵范数可以来度量这个变化大小。
在本说明书实施例中，可以将Lp范数作为稀疏约束项，根据所述反演方程建立目标函数。
具体的，可以将上述公式(3)写成矩阵-向量形式，可以得到：
Gm＝s (4)
其中，G＝[W(t,f1)W(t,f2)...W(t,fg)...W(t,fN)]是由N个子波褶积矩阵W(t,fg)组成的一个列大的矩阵，其中，W(t,fg)为w(t,fg)的向量形式；
是由不同频率或尺度下的谱或系数构成的一个长的列向量，其中，r(t,fg)为r(t,fg)的向量形式，g＝1,2,3...N，N＝1,2,3...。
在本说明书实施例中，由于矩阵G的行数普遍远小于列数，所以通过直接反演的方式获得的谱分解结果是不唯一的。由于谱分解的普遍目的是去获得高分辨率或聚焦性好的时频谱，即时频谱m期望是稀疏的。在一些实施例中，可以对谱分解结果m施加先验的稀疏约束去既减小甚至克服谱分解的多解问题同时提供高分辨率或聚焦性好的时频谱。
在一些实施例中，可以采用Lp范数最小来稀疏表示大尺度时频谱的稀疏性，可以建立如下目标函数：
其中，J1为稀疏约束项，用于控制解的稀疏度以及解的位置；J0为数据匹配项，用来衡量评价所求取的解与实际信号的匹配程度；||m||p为m的Lp范数，
为(Gm-s)的L2范数的平方。
在一些实施例中，所述Lp范数p的取值可以为0＜p≤1，以获得高分辨率或聚焦性好的时频谱，也可以取值为其他任意大于零的值，以获得不同分辨率或聚焦性的时频谱。
在一些实施例中地震信号中通常会出现噪声，考虑到噪声的影响，可将公式(5)改写为：
其中，σ为噪声因子，σ的取值范围可以是大于等于零，通常σ的值越大，则表明噪声的影响越大。
将公式(6)进一步转换为正则化问题，可得到如下公式：
其中，λ为稀疏度参数，用于控制稀疏约束在目标函数中所占的比重，即平衡两者之间的一个参数。当λ值相对较大时，稀疏约束所占的权重较大，其解比较稀疏，所得到的时频谱比较稀疏；当λ值相对较小时，其稀疏约束所占的权重较小，其解比较光滑。
基于以上，将Lp范数作为稀疏约束项，根据所述反演方程建立如下目标函数：
J＝J0+λJ1 (8)
其中，J为目标函数。
S140：使用快速迭代软阈值算法求解目标函数，得到目标区域的时频谱。
迭代软阈值算法(ISTA)是在稀疏信号处理中求解目标函数常用的算法，其实际上是软阈值跟梯度算法的结合。迭代软阈值算法在处理线性反演正则化问题上具有简单易行的特点，而且在其求解过程中不需求取矩阵的逆，被广泛地使用。迭代软阈值算法的收敛速度是一阶收敛，在对处理数据量比较大的正则化问题的时候，其需要大量的计算时间，但是其收敛速度较慢。
在本说明书实施例中，所述快速迭代软阈值算法是对迭代软阈值算法进行了优化提升的一种算法，所述快速迭代软阈值算法可以在迭代软阈值算法中引入了一个中间量，在每次迭代中，都可以利用前两次的迭代结果，从而加速了迭代软阈值法的收敛速度，为二阶收敛。所述快速迭代软阈值算法在每次迭代的过程中只需要计算一次梯度，而且其运算简单，且收敛速度明显比迭代软阈值法快。进一步的，通过快速迭代软阈值算法求解目标函数，可以获得高分辨、聚焦性好的时频谱属性。
在一些实施例中，所述快速迭代软阈值算法的第k阶Landweber算子的递推公式可以包括：
其中，
为迭代步长且大于零，“-”表示下降方向，mk为第k次迭代的m值，k＝1,2,3…。公式(9)可以等价写为相对于mk-1的成本函数：
其中，<m-mk-1,▽J0(mk-1)>＝(m-mk-1)T▽J0(mk-1)，将该梯度算法思想应用到目标函数即公式(8)中，其迭代过程可以包括：
忽略公式(11)中的常数项，可以将公式(11)简化为：
其中，公式(12)等价于求解如下成本函数的最优解：
其中，d为由公式(12)得到公式(13)的一个中间量，d与m具有相同的维数，对公式(13)求导并将之等于零可得：
其中，sgn(number)表示符号函数，用来表示输入的变量的正负，number参数是任何有效的数值表达式。参数值如果数字大于0，则sgn返回1，参数值如果数字等于0，则返回0，参数值如果数字小于0，则返回-1。其表达式如下：
在一些实施例中，可以根据以下非线性方程得到阈值
当p＝1时，可以得到L1范数约束的是Lp范数约束的一个特例，如公式(17)所示：
联立公式(10)-(17)可以得出，目标函数的解如下：
其中，其中α的值大于GHG的最大特征值，α与子波相关，子波不变，α的值也不改变；yk为快速迭代阈值算法中引入的一个中间量，yk是一个包含mk和mk-1的综合项，其表达式为：
其中，γk是一个随迭代更新的参数。正是由于yk的代入应用，使快速迭代阈值法的每一次迭代过程中，利用了两次的迭代的解，进而在一定程度上加速了该算法的速度。
在一些实施例中，可以通过设置地震信号的频率范围、以及设置稀疏度参数来获取目标区域的时频谱。
在一些实施例中，所述频率范围可以依据接收到的地震信号进行设定，例如确定地震信号中某一段频率内的地震信号，所述频率范围可以是0-200Hz、50-300Hz等，本说明书实施例对此不作限定。其中，所述频率范围可以是服务器预先设定的，也可以是工作人员进行设定，服务器可以接收工作人员设定的频率范围，服务器还可以通过其他任何方式确定频率范围，本说明书实施例对此不作限定。
在一些实施例中，可以通过稀疏度参数控制稀疏约束在目标函数中所占的比重，从而得到高分辨率或聚焦性好的时频谱。所述稀疏度参数可以是服务器预先设定的，也可以是工作人员进行设定，服务器可以接收工作人员设定的稀疏度参数，服务器还可以通过其他任何方式确定稀疏度参数，本说明书实施例对此不作限定。
S150：根据所述时频谱检测目标区域的含气性。
地震波在地下传播的过程中，部分能量会转化为热能从而引起能量衰减。其中，地震波能量衰减的原因还可以包括波阵面的扩展，通过界面时的反射与折射、通过不均匀介质时的散射以及地下岩层的非完全弹性相关的振幅吸收衰减。地震波在通过含油气层时，高频能量将发生明显的衰减，其中，含流体地层滞性衰减系数越大，地震波的主频越低。
在本说明书实施例中，可以获得高分辨、聚焦性好的时频谱属性，通过不同频率谱成分的对比检测低频阴影，能够消除强振幅的影响，从而能够更为准确的进行地下储层含气性的检测，降低储层勘探和开发的风险。
为了清楚的说明本说明书实施方式的有益效果，通过本说明书实施例提供的含气性检测方法选对某地区深部白云岩储层的的含气性进行了检测，下面结合图2-图6d进行说明。
图2为本说明书实施例含气性检测所在的研究工区图，其中黑色圆圈代表井位，分别用Well1、Well2和Well3来表示；黑线代表穿过三口井的任意线。图3为对应图2中任意线的连井地震剖面，其中T1、T2和T3为解释层位，其覆盖范围为主要储层段，分为上段储层(T1-T2)和下段储层(T2-T3)。钻井结果已证实，Well3井主要在上段储层产气，而Well1井和Well2井主要在下段储层产气，但是原始振幅却看不出任何明显的异常。
图4a、4b和4c展示了利用本说明书实施例提供的含气性检测方法得到的分频振幅谱剖面，其中图4a、4b、4c分别对应高频(35Hz)、中频(30Hz)、低频(21Hz)的振幅谱剖面，可以看到不同频率成分结果之间存在明显差异，并且在储层底部存在明显的低频振幅异常(图4b和4c中的椭圆区域)，可以推测该强振幅异常可能是由于储层含气使得高频能量衰减而低频能量保留所致，从而进一步推测储层含气的位置。
图5a、5b、5c和5d显示了上段储层含气性预测结果。其中，图5a为原始地震振幅(Amp)属性，图5b为原始地震甜点(Sweetness)属性，两者都表现出强的振幅异常且横向展布特征相似。图5c和图5d分别为30Hz和21Hz频率成分的振幅谱属性，可以很清楚识别出低频异常位置(图5c和5d中虚线圈出范围)，这些区域则可能为该段储层的含气区域。该预测结果与Well3钻井结果相吻合，证明了本说明书实施例提供的含气性检测方法的有效性。
同样，图6a、6b、6c和6d显示了下段储层含气性预测结果，其中，图6a为原始地震振幅属性，图6b为原始地震甜点属性，此时两者特征存在较大差异。图6c和图6d分别为35Hz和30Hz频率成分的振幅谱属性，同样可以识别出低频异常位置(图6c和图6d中虚线圈出范围)，这与Well1和Well2钻井结果相吻合。进一步说明了本说明书实施例提供的含气性检测方法可以消除原始地震数据中强振幅异常的干扰，能够更为准确进行地下储层含气性的检测。
本说明书实施例提供的基于时频分析的含气性检测方法，可以通过多维褶积模型构建反演方程，使用Lp准则约束条件产生稀疏时频分布，并且利用快速迭代软阈值算法(FISTA)用来求解Lp范数正则化问题，可以获得高分辨、聚焦性好的时频谱属性，进而可以通过不同频率谱成分的对比检测低频阴影，能够消除强振幅的影响，从而能够更为准确的进行地下储层含气性的检测，降低储层勘探和开发的风险。
本说明书实施例还提供了一种基于时频分析的含气性检测方法的计算机可读存储介质，所述计算机可读存储介质存储有计算机程序指令，在所述计算机程序指令被执行时实现：根据所述地震信号，通过褶积模型构建反演方程；将Lp范数作为稀疏约束项，根据所述反演方程建立目标函数；使用快速迭代软阈值算法求解目标函数，得到目标区域的时频谱；根据所述时频谱检测目标区域的含气性。
在本说明书实施例中，上述存储介质包括但不限于随机存取存储器(RandomAccess Memory,RAM)、只读存储器(Read-Only Memory,ROM)、缓存(Cache)、硬盘(HardDisk Drive,HDD)或者存储卡(Memory Card)。所述存储器可用于存储所述计算机程序和/或模块，所述存储器可主要包括存储程序区和存储数据区，其中，存储程序区可存储操作系统、至少一个功能所需的应用程序(比如声音播放功能、文字转换功能等)等；存储数据区可存储根据用户终端的使用所创建的数据(比如音频数据、文字消息数据等)等。此外，存储器可以包括高速随机存取存储器，还可以包括非易失性存储器。在本实施方式中，该计算机可读存储介质存储的程序指令具体实现的功能和效果，可以与其它实施方式对照解释，在此不再赘述。
请参阅图7，本说明书实施例还提供了一种基于时频分析的含气性检测装置，该装置具体可以包括以下的结构模块。
接收模块710，用于接收目标区域的地震信号；
构建模块720，用于根据所述地震信号，通过褶积模型构建反演方程；
建立模块730，用于将Lp范数作为稀疏约束项，根据所述反演方程建立目标函数；
求解模块740，用于使用快速迭代软阈值算法求解目标函数，得到目标区域的时频谱；
检测模块750，用于根据所述时频谱检测目标区域的含气性。
在一些实施例中，所述构建模块720还包括：第一构建子模块，用于根据所述地震信号，构建由不同频率的地震子波组成的子波库；第二构建子模块，用于根据所述地震信号和子波库，通过褶积模型构建反演方程。
需要说明的是，本说明书中的各个实施例均采用递进的方式描述，各个实施例之间相同或相似的部分互相参见即可，每个实施例重点说明的都是与其它实施例的不同之处。尤其，对于装置实施例和设备实施例而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。
本领域技术人员在阅读本说明书文件之后，可以无需创造性劳动想到将本说明书列举的部分或全部实施例进行任意组合，这些组合也在本说明书公开和保护的范围内。
在20世纪90年代，对于一个技术的改进可以很明显地区分是硬件上的改进(例如，对二极管、晶体管、开关等电路结构的改进)还是软件上的改进(对于方法流程的改进)。然而，随着技术的发展，当今的很多方法流程的改进已经可以视为硬件电路结构的直接改进。设计人员几乎都通过将改进的方法流程编程到硬件电路中来得到相应的硬件电路结构。因此，不能说一个方法流程的改进就不能用硬件实体模块来实现。例如，可编程逻辑器件(Programmable Logic Device,PLD)(例如现场可编程门阵列(Field Programmable GateArray，FPGA))就是这样一种集成电路，其逻辑功能由用户对器件编程来确定。由设计人员自行编程来把一个数字系统“集成”在一片PLD上，而不需要请芯片制造厂商来设计和制作专用的集成电路芯片。而且，如今，取代手工地制作集成电路芯片，这种编程也多半改用“逻辑编译器(logic compiler)”软件来实现，它与程序开发撰写时所用的软件编译器相类似，而要编译之前的原始代码也得用特定的编程语言来撰写，此称之为硬件描述语言(Hardware Description Language，HDL)，而HDL也并非仅有一种，而是有许多种，如ABEL(Advanced Boolean Expression Language)、AHDL(AlteraHardware DescriptionLanguage)、Confluence、CUPL(Cornell University Programming Language)、HDCal、JHDL(Java Hardware Description Language)、Lava、Lola、MyHDL、PALASM、RHDL(RubyHardware Description Language)等，目前最普遍使用的是VHDL(Very-High-SpeedIntegrated Circuit Hardware Description Language)与Verilog2。本领域技术人员也应该清楚，只需要将方法流程用上述几种硬件描述语言稍作逻辑编程并编程到集成电路中，就可以很容易得到实现该逻辑方法流程的硬件电路。
上述实施例阐明的系统、装置、模块或单元，具体可以由计算机芯片或实体实现，或者由具有某种功能的产品来实现。一种典型的实现设备为计算机。具体的，计算机例如可以为个人计算机、膝上型计算机、蜂窝电话、相机电话、智能电话、个人数字助理、媒体播放器、导航设备、电子邮件设备、游戏控制台、平板计算机、可穿戴设备或者这些设备中的任何设备的组合。
通过以上的实施方式的描述可知，本领域的技术人员可以清楚地了解到本说明书可借助软件加必需的通用硬件平台的方式来实现。基于这样的理解，本说明书的技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品可以存储在存储介质中，如ROM/RAM、磁碟、光盘等，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本说明书各个实施例或者实施例的某些部分所述的方法。
本说明书中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于系统实施例而言，由于其基本相似于方法实施例，所以描述的比较简单，相关之处参见方法实施例的部分说明即可。
本说明书可用于众多通用或专用的计算机系统环境或配置中。例如：个人计算机、服务器计算机、手持设备或便携式设备、平板型设备、多处理器系统、基于微处理器的系统、置顶盒、可编程的消费电子设备、网络PC、小型计算机、大型计算机、包括以上任何系统或设备的分布式计算环境等等。
本说明书可以在由计算机执行的计算机可执行指令的一般上下文中描述，例如程序模块。一般地，程序模块包括执行特定任务或实现特定抽象数据类型的例程、程序、对象、组件、数据结构等等。也可以在分布式计算环境中实践本说明书，在这些分布式计算环境中，由通过通信网络而被连接的远程处理设备来执行任务。在分布式计算环境中，程序模块可以位于包括存储设备在内的本地和远程计算机存储介质中。
虽然通过实施例描绘了本说明书，本领域普通技术人员知道，本说明书有许多变形和变化而不脱离本说明书的精神，希望所附的权利要求包括这些变形和变化而不脱离本说明书的精神。
Claims (7)
1.一种基于时频分析的含气性检测方法，其特征在于，所述方法包括：
接收目标区域的地震信号；
根据所述地震信号，通过褶积模型构建反演方程；
将Lp范数作为稀疏约束项，根据所述反演方程建立目标函数；所述目标函数包括：
J＝J0+J1
其中，J表示目标函数，J0表示资料匹配项，
表示Gm-s的L2范数的平方，s表示地震信号s(t)的向量形式，G＝[W(t,f1)W(t,f2)...W(t,fg)...W(t,fN)]是由N个子波褶积矩阵W(t,fg)组成的矩阵，其中，W(t,fg)表示w(t,fg)的向量形式，w(t,fg)表示由地震子波组成的子波库或基函数，t表示时间，fg表示频率；
是由不同频率或尺度下的谱或系数构成的一个长的列向量，其中，r(t,fg)表示r(t,fg)的向量形式，r(t,fg)表示某一频率或尺度下的谱或系数；J1为稀疏约束项，J1＝||m||p，||m||p为m的Lp范数，λ为稀疏度参数，g＝1,2,3...N，N＝1,2,3...，0＜p<1；
使用快速迭代软阈值算法求解目标函数，得到目标区域的时频谱；包括：根据以下非线性方程得到阈值
其中λ表示稀疏度参数，目标函数的解如下：
其中，k表示迭代次数，
表示迭代步长且大于零，α的值大于GHG的最大特征值，GH表示矩阵G的共轭转置矩阵，α与子波相关，子波不变，α的值也不改变，soft表示软阈值函数，sgn表示符号函数；yk为快速迭代阈值算法中引入的一个中间量，yk是一个包含mk和mk-1的综合项，其表达式为：
其中，γ是一个随迭代更新的参数，下标表示迭代次数；
根据所述时频谱检测目标区域的含气性。
2.根据权利要求1所述的方法，其特征在于，所述根据所述地震信号，通过褶积模型构建反演方程包括：
根据所述地震信号，构建由不同频率的地震子波组成的子波库；
根据所述地震信号和子波库，通过褶积模型构建反演方程。
3.根据权利要求2所述的方法，其特征在于，根据所述地震信号和子波库，通过褶积模型构建如下反演方程：
其中，*表示褶积运算，s(t)表示地震信号，t表示时间，fg表示频率，w(t,fg)表示由地震子波组成的子波库或基函数，r(t,fg)表示某一频率或尺度下的谱或系数，其中，N＝1,2,3...。
4.一种基于时频分析的含气性检测装置，其特征在于，所述装置包括：
接收模块，用于接收目标区域的地震信号；
构建模块，用于根据所述地震信号，通过褶积模型构建反演方程；
建立模块，用于将Lp范数作为稀疏约束项，根据所述反演方程建立目标函数；所述目标函数包括：
J＝J0+J1
其中，J表示目标函数，J0表示资料匹配项，
表示Gm-s的L2范数的平方，s表示地震信号s(t)的向量形式，G＝[W(t,f1)W(t,f2)...W(t,fg)...W(t,fN)]是由N个子波褶积矩阵W(t,fg)组成的矩阵，其中，W(t,fg)表示w(t,fg)的向量形式，w(t,fg)表示由地震子波组成的子波库或基函数，t表示时间，fg表示频率；
是由不同频率或尺度下的谱或系数构成的一个长的列向量，其中，r(t,fg)表示r(t,fg)的向量形式，r(t,fg)表示某一频率或尺度下的谱或系数；J1为稀疏约束项，J1＝||m||p，||m||p为m的Lp范数，λ为稀疏度参数，g＝1,2,3...N，N＝1,2,3...，0＜p<1；
求解模块，用于使用快速迭代软阈值算法求解目标函数，得到目标区域的时频谱；包括：根据以下非线性方程得到阈值
其中λ表示稀疏度参数，目标函数的解如下：
其中，k表示迭代次数，
表示迭代步长且大于零，α的值大于GHG的最大特征值，GH表示矩阵G的共轭转置矩阵，α与子波相关，子波不变，α的值也不改变，soft表示软阈值函数，sgn表示符号函数；yk为快速迭代阈值算法中引入的一个中间量，yk是一个包含mk和mk-1的综合项，其表达式为：
其中，γ是一个随迭代更新的参数，下标表示迭代次数；
检测模块，用于根据所述时频谱检测目标区域的含气性。
5.根据权利要求4所述的装置，其特征在于，所述构建模块包括：
第一构建子模块，用于根据所述地震信号，构建由不同频率的地震子波组成的子波库；
第二构建子模块，用于根据所述地震信号和子波库，通过褶积模型构建反演方程。
6.根据权利要求5所述的装置，其特征在于，根据所述地震信号和子波库，通过褶积模型构建如下反演方程：
其中，*表示褶积运算，s(t)表示地震信号，t表示时间，fg表示频率，w(t,fg)表示由地震子波组成的子波库或基函数，r(t,fg)表示某一频率或尺度下的谱或系数，其中，N＝1,2,3...。
7.一种计算机可读存储介质，其上存储有计算机程序指令，所述计算机程序指令被执行时实现：接收目标区域的地震信号；根据所述地震信号，通过褶积模型构建反演方程；将Lp范数作为稀疏约束项，根据所述反演方程建立目标函数；所述目标函数包括：
J＝J0+J1
其中，J表示目标函数，J0表示资料匹配项，
表示Gm-s的L2范数的平方，s表示地震信号s(t)的向量形式，G＝[W(t,f1)W(t,f2)...W(t,fg)...W(t,fN)]是由N个子波褶积矩阵W(t,fg)组成的一个列大的矩阵，其中，W(t,fg)表示w(t,fg)的向量形式，w(t, fg)表示由地震子波组成的子波库或基函数，t表示时间，fg表示频率；
是由不同频率或尺度下的谱或系数构成的一个长的列向量，其中，r(t,fg)表示r(t,fg)的向量形式，r(t,fg)表示某一频率或尺度下的谱或系数；J1为稀疏约束项，J1＝||m||p，||m||p为m的Lp范数，λ为稀疏度参数，g＝1,2,3...N，N＝1,2,3...，0＜p<1；使用快速迭代软阈值算法求解目标函数，得到目标区域的时频谱；包括：根据以下非线性方程得到阈值
其中λ表示稀疏度参数，目标函数的解如下：
其中，k表示迭代次数，
表示迭代步长且大于零，α的值大于GHG的最大特征值，GH表示矩阵G的共轭转置矩阵，α与子波相关，子波不变，α的值也不改变，soft表示软阈值函数，sgn表示符号函数；yk为快速迭代阈值算法中引入的一个中间量，yk是一个包含mk和mk-1的综合项，其表达式为：
其中，γ是一个随迭代更新的参数，下标表示迭代次数；
根据所述时频谱检测目标区域的含气性。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911042135.3A CN110794458B (zh) | 2019-10-30 | 2019-10-30 | 基于时频分析的含气性检测方法、装置及存储介质 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201911042135.3A CN110794458B (zh) | 2019-10-30 | 2019-10-30 | 基于时频分析的含气性检测方法、装置及存储介质 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN110794458A CN110794458A (zh) | 2020-02-14 |
| CN110794458B true CN110794458B (zh) | 2021-05-25 |
Family
ID=69442136
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201911042135.3A Active CN110794458B (zh) | 2019-10-30 | 2019-10-30 | 基于时频分析的含气性检测方法、装置及存储介质 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN110794458B (zh) |
Families Citing this family (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114428329A (zh) * | 2020-09-22 | 2022-05-03 | 中国石油化工股份有限公司 | 地质信息的检测方法、设备和存储介质 |
| CN112213782B (zh) * | 2020-09-29 | 2022-03-04 | 中国石油大学(北京) | 分相位地震数据的处理方法、装置和服务器 |
| CN113777650B (zh) * | 2021-08-12 | 2022-10-25 | 西安交通大学 | 一种基于混合范数和小波变换的稀疏时频谱分解方法、装置、设备及存储介质 |
| CN115267898B (zh) * | 2022-08-11 | 2024-06-14 | 河北地质大学 | 天然地震数据重建方法、装置和电子设备 |
Family Cites Families (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| AU2014200168B2 (en) * | 2013-01-14 | 2018-02-15 | Cgg Services Sa | High-fidelity adaptive curvelet domain primary-multiple separation processing of seismic data |
| WO2016162724A1 (en) * | 2015-04-10 | 2016-10-13 | Cgg Services Sa | Methods and data processing apparatus for deblending seismic data |
| CN106842298A (zh) * | 2015-12-04 | 2017-06-13 | 中国石油化工股份有限公司 | 一种基于匹配追踪的不整合强反射自适应分离方法 |
| CN110174702B (zh) * | 2018-09-30 | 2020-10-16 | 中海油海南能源有限公司 | 一种海上地震数据低频弱信号恢复的方法和系统 |
| CN109946739A (zh) * | 2019-03-15 | 2019-06-28 | 成都理工大学 | 一种基于压缩感知理论的地震剖面增强方法 |
-
2019
- 2019-10-30 CN CN201911042135.3A patent/CN110794458B/zh active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN110794458A (zh) | 2020-02-14 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN110794458B (zh) | 基于时频分析的含气性检测方法、装置及存储介质 | |
| Gebraad et al. | Bayesian elastic full‐waveform inversion using Hamiltonian Monte Carlo | |
| Zhang et al. | Unsplit complex frequency-shifted PML implementation using auxiliary differential equations for seismic wave modeling | |
| Anderson et al. | Time-reversal checkpointing methods for RTM and FWI | |
| Zhu et al. | A Bayesian approach to estimate uncertainty for full-waveform inversion using a priori information from depth migration | |
| Burstedde et al. | Algorithmic strategies for full waveform inversion: 1D experiments | |
| Ji et al. | Source description of the 1999 Hector Mine, California, earthquake, part I: Wavelet domain inversion theory and resolution analysis | |
| US10054714B2 (en) | Fast viscoacoustic and viscoelastic full wavefield inversion | |
| Panning et al. | Importance of crustal corrections in the development of a new global model of radial anisotropy | |
| Gosselin-Cliche et al. | 3D frequency-domain finite-difference viscoelastic-wave modeling using weighted average 27-point operators with optimal coefficients | |
| US10877175B2 (en) | Seismic acquisition geometry full-waveform inversion | |
| CN113945994B (zh) | 使用有限差分模型进行高速多源加载和波场检索的方法 | |
| Vavryčuk | Ray velocity and ray attenuation in homogeneous anisotropic viscoelastic media | |
| Sjögreen et al. | Source estimation by full wave form inversion | |
| Lyu et al. | Efficiency of the spectral element method with very high polynomial degree to solve the elastic wave equation | |
| Li et al. | Pore type identification in carbonate rocks using convolutional neural network based on acoustic logging data | |
| Favorskaya et al. | The use of full-wave numerical simulation for the investigation of fractured zones | |
| Fabien-Ouellet | Seismic modeling and inversion using half-precision floating-point numbers | |
| Yablokov et al. | Uncertainty quantification of multimodal surface wave inversion using artificial neural networks | |
| Ojeda et al. | Effect of windowing on lithosphere elastic thickness estimates obtained via the coherence method: Results from northern South America | |
| WO2021257097A1 (en) | Acoustic dispersion curve identification based on reciprocal condition number | |
| Liu et al. | Quantitative analysis of near‐surface seismological complexity based on the generalized Lippmann–Schwinger matrix equation | |
| US20220390631A1 (en) | Laplace-fourier 1.5d forward modeling using an adaptive sampling technique | |
| Walters et al. | Analytic and numerical solutions to the seismic wave equation in continuous media | |
| Yuan et al. | Finite‐difference modeling and characteristics analysis of Love waves in anisotropic‐viscoelastic media |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |