CN114994758B - 碳酸盐岩断控储层的波阻抗提取与结构表征方法和系统 - Google Patents

碳酸盐岩断控储层的波阻抗提取与结构表征方法和系统 Download PDF

Info

Publication number
CN114994758B
CN114994758B CN202210919668.0A CN202210919668A CN114994758B CN 114994758 B CN114994758 B CN 114994758B CN 202210919668 A CN202210919668 A CN 202210919668A CN 114994758 B CN114994758 B CN 114994758B
Authority
CN
China
Prior art keywords
wave impedance
seismic
model
data
stratum
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
Application number
CN202210919668.0A
Other languages
English (en)
Other versions
CN114994758A (zh
Inventor
赵艳萍
崔建兵
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beijing Jinglu Juyuan Energy Technology Co ltd
Original Assignee
Beijing Jinglu Juyuan Energy Technology Co ltd
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Beijing Jinglu Juyuan Energy Technology Co ltd filed Critical Beijing Jinglu Juyuan Energy Technology Co ltd
Priority to CN202210919668.0A priority Critical patent/CN114994758B/zh
Publication of CN114994758A publication Critical patent/CN114994758A/zh
Application granted granted Critical
Publication of CN114994758B publication Critical patent/CN114994758B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

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/40—Seismology; Seismic or acoustic prospecting or detecting specially adapted for well-logging
    • 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/622—Velocity, density or impedance
    • G01V2210/6226—Impedance
    • 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

  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Remote Sensing (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

本发明属于地质分析领域，具体涉及了一种碳酸盐岩断控储层的波阻抗提取与结构表征方法和系统，旨在解决常规波阻抗反演方法难以体现垂向断裂带储层发育的地质规律，地球物理反演难度高，多解性强，储层识别结果不准确的问题。本发明包括：根据地震资料和测井数据构建目标地层的等时地层格架模型，再获取波阻抗反演数据体；根据野外地质露头资料和时深转化关系计算背景地层波阻抗模型；将背景地层波阻抗模型与波阻抗反演数据体进行做差运算，获得波阻抗异常数据体；根据波阻抗异常数据体进行部分数据剔除，解释获得碳酸盐岩断控储层解释模型。本发明通过将所述背景地层波阻抗模型与波阻抗反演数据体进行做差运算的方式提高了储层识别的准确性。

Description

碳酸盐岩断控储层的波阻抗提取与结构表征方法和系统
技术领域
本发明属于地质分析领域，具体涉及了一种碳酸盐岩断控储层的波阻抗提取与结构表征方法和系统。
背景技术
在长期的碳酸盐岩储层勘探过程中，风化壳岩溶控制的缝洞体作为主要的油气开采目标。但会面临油井寿命短、综合递减大、钻探成功率低等生产问题。在不断地勘探实践过程中，发现深大断裂带不仅是油源运输的主要通道，而且与构造变形、断控岩溶作用耦合可以形成有效的油气圈闭，储层连通性较好，成为近期油气开采的重要目标。
目前的断控岩溶储层识别领域中，叠前或叠后地震属性数据对断控岩溶储层的结构框定起到了关键性作用。其中包括相干属性、倾角属性、AFE属性、结构张量属性、蚂蚁体属性等，地震属性间的融合与建模在一定程度上使断裂带解释更加准确，符合地质规律。然而断控储层内部储层类型与组合关系是影响油气储集性能的关键要素。在使用常规波阻抗反演技术时，研究人员发现，由于断溶体受到断裂、岩溶作用的控制，储层具有断裂样式多样、裂缝与溶蚀孔洞组合关系复杂、纵横向非均质性极强的特点。再加上地震资料横向地层界面反射干扰性大，常规波阻抗反演方法难以体现垂向断裂带储层发育的地质规律，因此常规地球物理反演难度高，多解性强，储层识别结果不准确。
碳酸盐岩断裂带垂向发育与横向稳定沉积的地质规律存在显著差别，并体现在地震波形的横向变化中。本发明利用了该地质与地球物理规律，将沉积相特征与断控岩溶储层发育模式的地质因素考虑在内，对断裂带区域进行针对性的反演研究是解决复杂的断控岩溶储层识别与划分的重要手段。
发明内容
为了解决现有技术中的上述问题，常规波阻抗反演方法难以体现垂向断裂带储层发育的地质规律，因此常规地球物理反演难度高，多解性强，储层识别结果不准确的问题，本发明提供了一种碳酸盐岩断控储层的波阻抗提取与结构表征方法和系统，所述方法包括：
获取已知井位的原始地震资料和测井资料；
基于所述原始地震资料和测井资料，通过谱整形和扩散滤波的方式，获得提频降噪后的地震资料，并进行井震标定获得时深转化关系；
基于所述提频降噪后的地震资料，构建目标层位的等时地层格架模型，并计算断溶体几何结构和空间分布；
基于所述测井资料和提频降噪后的地震资料，在所述等时地层格架模型提供的地层趋势信息基础上，进行波形指示地震波阻抗反演获得波阻抗反演数据体；
基于野外地质露头资料，确定目标层位的地层密度序列；
基于所述时深转化关系，确认所述地层密度序列与等时地层格架模型的纵向对应关系，进而根据测井资料构建地层密度模型和地层声波时差模型；
基于所述地层密度模型和地层声波时差模型，计算背景地层波阻抗模型；
将所述背景地层波阻抗模型与波阻抗反演数据体进行做差运算，获得波阻抗异常数据体；
基于断溶体几何结构和空间分布，在所述波阻抗异常数据体中仅保留断溶体轮廓内的波阻抗异常数据体中的数据点，获得断溶体几何结构-波阻抗模型；
对比所述断溶体几何结构-波阻抗模型与已知井位中的断溶体井位的测井解释结果，划定孔洞储层特征值区间、断裂带特征值区间和围岩特征值区间，获得碳酸盐岩断控储层解释模型。
在一些优选的实施方式中，所述构建目标层位的等时地层格架模型，其获取方法为：
基于所述提频降噪后的地震资料，对标志层位地震波反射同相轴进行追踪显示，获得标志层层位解释连续面；
基于所述标志层层位解释连续面，预设时窗、切片方向和最大深度，构建目标层位的等时地层格架模型。
时深转换关系指示了地震波反射界面随地震波到时的变化关系；
在一些优选的实施方式中，所述波阻抗反演数据体，其获取方式包括：
基于所述提频降噪后的地震资料，计算待判别地震道波形与已知井的合成地震记录的波形相关性，根据波形相关性最高的井对应的波阻抗曲线建立初始模型；
利用白噪声满足高斯分布的规律，将测井资料中的波阻抗曲线表示为：
Figure 753796DEST_PATH_IMAGE001
其中，
Figure 405357DEST_PATH_IMAGE002
表示测井波阻抗曲线，
Figure 919515DEST_PATH_IMAGE003
表示待求解的地下地层实际波阻抗值，N表示 随机噪声；
根据中心极限定理，
Figure 908200DEST_PATH_IMAGE004
也满足高斯分布，确定初始目标函数为：
Figure 456993DEST_PATH_IMAGE005
其中，
Figure 13876DEST_PATH_IMAGE006
表示与后验信息有关的函数，
Figure 139964DEST_PATH_IMAGE007
表示基于样本井的波阻抗曲线，进行匹配 滤波后，求得后验概率统计分布密度，进而计算得到的波阻抗期望值，
Figure 73285DEST_PATH_IMAGE008
表示白噪声的协方 差；
基于所述初始目标函数，通过最大后验估计，在目标函数中引入先验信息，获得稳定的目标函数为：
Figure 476584DEST_PATH_IMAGE009
其中，
Figure 63423DEST_PATH_IMAGE010
表示待模拟的特征参数，
Figure 552173DEST_PATH_IMAGE011
表示与地质和测井资料先验信息有关的函 数，
Figure 23606DEST_PATH_IMAGE012
表示用于协调
Figure 281412DEST_PATH_IMAGE013
和
Figure 570311DEST_PATH_IMAGE011
之间的相互影响的平滑参数；
以所述稳定的目标函数作为初始模型的输入，通过马尔科夫链蒙特卡罗方法MCMC和Metropolis-Hastings抽样准则对后验概率分布抽样，不断优化初始模型的参数，选取目标函数取最大值时的解作为随机实现，取多次随机实现的均值作为期望值输出，将所述期望值输出作为波阻抗反演数据体。
在一些优选的实施方式中，根据稳定的目标函数获取波阻抗反演数据体，具体包括：
设M为目标空间，n为总样本数，m为马尔科夫链趋于达到平稳时的样本数；
预设一条马尔科夫链，使马尔科夫链收敛至达到平稳；
由M中的某一点
Figure 280778DEST_PATH_IMAGE014
出发，通过马尔科夫链进行抽样模拟，产生点序列：
Figure 290322DEST_PATH_IMAGE015
；
函数
Figure 792848DEST_PATH_IMAGE016
的期望估计为：
Figure 862435DEST_PATH_IMAGE017
其中，n为总样本数，m表示马尔科夫链达到平稳时的样本数，k表示累加参量；
选取一转移函数
Figure 60198DEST_PATH_IMAGE018
和初始值
Figure 998067DEST_PATH_IMAGE019
，若第i次迭代开始时的参数值为
Figure 230465DEST_PATH_IMAGE020
， 则第h次迭代过程为：
从
Figure 736533DEST_PATH_IMAGE018
中抽取一个备选值
Figure 421592DEST_PATH_IMAGE021
，计算备选值
Figure 163152DEST_PATH_IMAGE021
的接受概率
Figure 250057DEST_PATH_IMAGE022
：
Figure 927026DEST_PATH_IMAGE023
以
Figure 99381DEST_PATH_IMAGE024
，置
Figure 116403DEST_PATH_IMAGE025
；以概率
Figure 323394DEST_PATH_IMAGE026
，置
Figure 171264DEST_PATH_IMAGE027
；
不断扰动所述初始模型的参数，重复选取M中的某一点
Figure 955549DEST_PATH_IMAGE014
求解后验样本
Figure 914278DEST_PATH_IMAGE028
，直至 达到预设的迭代次数
Figure 710196DEST_PATH_IMAGE029
，获得后验样本
Figure 728967DEST_PATH_IMAGE030
，进而计算后验分布的期望输 出值，将期望值输出作为波阻抗反演数据体。
在一些优选的实施方式中，所述地层密度模型和地层声波时差模型，其获取方法为：
基于所述时深转化关系，建立地层密度序列与地震波到时Tsm的关系；依据地震波到时匹配地层密度序列与等时地层格架模型的纵向对应关系，得到等时地层格架模型内部的密度分布，即地层密度模型；
基于所述时深转化关系，建立测井资料中声波时差测井曲线与地震波到时Tsm的关系；依据地震波到时匹配地层声波时差与等时地层格架模型的纵向对应关系，得到等时地层格架模型内部的声波时差分布，即地层声波时差模型。
在一些优选的实施方式中，所述背景地层波阻抗模型为：
Figure 734969DEST_PATH_IMAGE031
其中，IMP为波阻抗值，DEN为密度值，vp为纵波波速。
在一些优选的实施方式中，所述断溶体几何结构-波阻抗模型，其获方法为：
基于所述提频降噪后的地震资料计算地震信号方差属性值数据体，剔除地震信号方差值小于预设的第二阈值的区域得到断溶体轮廓；保留断溶体轮廓内的波阻抗异常数据体中的数据点，获得断溶体几何结构-波阻抗模型。
对比所述断溶体几何结构-波阻抗模型与已知井位中的断溶体井位的测井解释结果，划定孔洞储层特征值区间、断裂带特征值区间和围岩特征值区间，获得碳酸盐断岩控储层解释模型。
本发明的另一方面，提出了一种碳酸盐岩断控储层的波阻抗提取与结构表征系统，其特征在于，所述提取系统包括：资料获取模块、地震提频和井震标定模块、地层结构获取模块、波阻抗反演数据体获取模块、地层密度序列确定模块、地层声波时差模型获取模块、背景地层波阻抗模型获取模块、波阻抗异常数据体获取模块、断溶体几何结构-波阻抗模型和碳酸盐岩断控储层解释模型获取模块；
所述资料获取模块，配置为获取已知井位的原始地震资料和测井资料；
所述地震提频和井震标定模块，配置为基于所述原始地震资料和测井资料，通过谱整形和扩散滤波的方式，获得提频降噪后的地震资料，并进行井震标定获得时深转化关系；
所述地层结构获取模块，配置为基于所述提频降噪后的地震资料，构建目标层位的等时地层格架模型，并计算断溶体几何结构和空间分布；
所述波阻抗反演数据体获取模块，配置为基于所述测井资料和提频降噪后的地震资料，在所述等时地层格架模型提供的地层趋势信息基础上，进行波形指示地震波阻抗反演获得波阻抗反演数据体；
所述地层密度序列确定模块，配置为基于野外地质露头资料，确定目标层位的地层密度序列；
所述地层声波时差模型获取模块，配置为基于所述时深转化关系，确认所述地层密度序列与等时地层格架模型的纵向对应关系，进而根据测井资料构建地层密度模型和地层声波时差模型；
所述背景地层波阻抗模型获取模块，配置为基于所述地层密度模型和地层声波时差模型，计算背景地层波阻抗模型；
所述波阻抗异常数据体获取模块，配置为将所述背景地层波阻抗模型与波阻抗反演数据体进行做差运算，获得波阻抗异常数据体；
所述断溶体几何结构-波阻抗模型，配置为基于断溶体几何结构和空间分布，在所述波阻抗异常数据体中仅保留断溶体几何轮廓对应区域的数据，获得断溶体几何结构-波阻抗模型；
所述碳酸盐岩断控储层解释模型获取模块，配置为对比所述断溶体几何结构-波阻抗模型与已知井位中的断溶体井位的测井解释结果，划定孔洞储层特征值区间、断裂带特征值区间和围岩特征值区间，获得碳酸盐岩断控储层解释模型。
本发明的有益效果：
本发明通过波形指示反演算法，挖掘了地震波形的横向变化的细节信息，并转化为波阻抗数据，突出了断溶体横向上由岩溶洞穴到裂缝、孔洞的渐变特征；此外，利用高精度的测井资料计算井旁波阻抗信息用于模拟稳定沉积的波阻抗背景，与实际反演结果对比，突出了断溶体纵向发育特征，避免了传统反演方法对于垂向发育储层识别难度高的问题。
附图说明
通过阅读参照以下附图所作的对非限制性实施例所作的详细描述，本申请的其它特征、目的和优点将会变得更明显：
图1是本发明实施例中碳酸盐岩断控储层的波阻抗提取与结构表征方法的流程示意图；
图2是断溶体几何结构和空间分布；
图3是目标层位的等时地层格架；
图4是波阻抗反演数据体；
图5是碳酸盐岩断控储层解释模型。
具体实施方式
下面结合附图和实施例对本申请作进一步的详细说明。可以理解的是，此处所描述的具体实施例仅用于解释相关发明，而非对该发明的限定。另外还需要说明的是，为了便于描述，附图中仅示出了与有关发明相关的部分。
需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。下面将参考附图并结合实施例来详细说明本申请。
本发明提供一种碳酸盐岩断控储层的波阻抗提取与结构表征方法，通过将所述背景地层波阻抗模型与波阻抗反演数据体进行做差运算的方式，避免了因为地震资料横向地层界面反射干扰性大，常规波阻抗反演方法难以体现垂向断裂带储层发育的地质规律带来的地球物理反演难度高，多解性强，储层识别结果不准确的问题。
中国深层海相碳酸盐岩储层油气资源丰富，在塔里木、鄂尔多斯、四川和渤海湾等陆上海相盆地已发现大量的石油地质资源，所探明石油储量为340×108t，然而，探明率仅为13.17%。
塔里木盆地奥陶系超深层缝洞型碳酸盐岩大油区作为重要代表，呈现连片含油趋势。盆地超深层的石油资源量为34.5×108t，天然气资源量为5.98×1012m3，分别占盆地石油和天然气资源量的46%和51%，是中国最重要的海相碳酸盐岩储层发育地之一。
为了更清晰地对本发明系统进行说明，下面结合图1对本发明实施例中各步骤展开详述。
本发明第一实施例的碳酸盐岩断控储层的波阻抗提取与结构表征方法，包括：
S100，获取已知井位的原始地震资料和测井资料；
S200，基于所述原始地震资料和测井资料，通过谱整形和扩散滤波的方式，获得提频降噪后的地震资料，并进行井震标定获得时深转化关系；
在本实施例中，原始地震资料以15m×15m的道间距设计三维测网，地震信号由2ms采样间隔的地震波振幅数据点组成，总采样时长为6ms。由于本发明用于断控岩溶储层的精细识别，对于原始地震资料主频要求在50~60Hz范围内。当目的层段原始地震资料主频低于50Hz时，需要采用混合相位子波反褶积（mixed-phase wavelet estimation and maximumposteriori deconvolution）与扩散滤波（enhancing diffusion filtering）技术提高地震资料主频、扩展地震资料的有效频带范围，获得提频降噪后的地震资料；
根据标志层层位解释连续面进行三维地质格架的构建，即设定时窗为0.25ms，平行于层位顶设计地层切片，最大深度设置为距层位顶面以下1000m，据此构建等时地层格架模型；测井资料是描述不同深度段仪器探测参数的变化，而地震资料是描述地震波反射能量随地震波双程旅行时的变化关系；须通过井震标定的方式将这两种地球物理数据结合。
在本实施例中，所述提频降噪后的地震资料，其获取方法为：
地震资料提频包括：谱整形和扩散滤波；
谱整形：
步骤A100，将原始地震资料的地震记录褶积模型在频率域的表示为：
Figure 497389DEST_PATH_IMAGE032
其中，
Figure 413392DEST_PATH_IMAGE033
表示傅氏变换后的地震记录频谱，
Figure 462120DEST_PATH_IMAGE034
表示傅氏变换后的地震子波频 谱，
Figure 96364DEST_PATH_IMAGE035
表示傅氏变换后的反射系数频谱，
Figure 396895DEST_PATH_IMAGE036
表示角频率；
步骤A200，将所述地震记录褶积模型在频率域的表示转化为地震记录线性系统频域的表示：
Figure 432984DEST_PATH_IMAGE037
其中，
Figure 652613DEST_PATH_IMAGE038
为地震记录线性系统频域表示，
Figure 508573DEST_PATH_IMAGE039
为地震子波线性系统频域表 示，
Figure 612796DEST_PATH_IMAGE040
为反射系数线性系统频域表示；
步骤A300，将所述地震记录线性系统频域的表示进行反傅氏变换，获得地震记录频谱的复赛谱序列：
Figure 362446DEST_PATH_IMAGE041
其中，
Figure 893921DEST_PATH_IMAGE042
表示地震记录频谱的复赛谱序列，
Figure 237178DEST_PATH_IMAGE043
表示地震子波频谱的复赛谱序 列，
Figure 145091DEST_PATH_IMAGE044
表示反射系数频谱的复赛谱序列，t表示地震波到时；
步骤A400，通过低通滤波器将复赛谱中的子波复赛谱序列和反射系数复赛谱序列进行分离，提取子波振幅；本步骤根据子波与反射稀疏列平滑程度的差异在复赛谱中易于区分的特性：子波能量出现于原点附近，而反射系数序列则远离原点。利用低通滤波器就可以将复赛谱中的子波与反射系数分离达到子波振幅谱提取的目的；
步骤A500，基于所述子波振幅，通过最小二乘法模拟地震子波振幅谱：
Figure 14827DEST_PATH_IMAGE045
其中，
Figure 451625DEST_PATH_IMAGE046
表示大于等于0的常数，
Figure 282177DEST_PATH_IMAGE047
表示地震子波频谱，由
Figure 852836DEST_PATH_IMAGE048
经傅氏变换而来，
Figure 452444DEST_PATH_IMAGE049
和
Figure 325723DEST_PATH_IMAGE050
为待求的关于频率
Figure 643571DEST_PATH_IMAGE051
的多项式用于拟合地震子波振幅谱；
步骤A600，基于所述模拟地震子波振幅谱，获得子波最大相位分量和最小相位分量；
设子波
Figure 585350DEST_PATH_IMAGE052
的最大相位分量为
Figure 39465DEST_PATH_IMAGE053
、最小相位分量为
Figure 83644DEST_PATH_IMAGE054
，则子波
Figure 13423DEST_PATH_IMAGE052
为：
Figure 801250DEST_PATH_IMAGE055
振幅谱的复赛谱中表示为：
Figure 375451DEST_PATH_IMAGE056
其中，振幅谱的复赛谱
Figure 590532DEST_PATH_IMAGE057
在复赛谱的正、负轴上对称显示，
Figure 742027DEST_PATH_IMAGE058
为地震 子波最大相位分量
Figure 333546DEST_PATH_IMAGE059
所对应的最小相位函数的复赛谱，
Figure 762253DEST_PATH_IMAGE060
为地震子波最小相 位分量
Figure 7290DEST_PATH_IMAGE061
所对应的最大相位函数的复赛谱；
步骤A700，基于所述振幅谱的复赛谱确定一组具有相同振幅谱的混合相位子波集合，调整俞氏子波参数，在保证地震子波主频完整的前提下，提升有效带宽至预设的带宽阈值，获得整形后波形数据；优选带宽阈值为50Hz；
经过谱整形后，拓展了地震数据的有效频带，高频部分得到了合理加强。在地震波形上表现为同相轴数目增多，更容易反映地震波反射信息的细节变化，并且在振幅、相位、频率方面改善了同一反射波组波形的一致性。在缝洞体地震响应上，“串珠状”反射特征尤为明显，且串珠内部形态的细节能够清晰地显示，代表了不同结构特征与充填物组合的复杂缝洞型储层地震反射，有助于后期精细地质解释。
扩散滤波：
Fhemers和Hocker与2003年首次将扩散滤波技术应用在地震资料的处理解释中。该技术不仅可以有效压制噪声，而且可以尽可能保留地震数据中的细节：如地质体边缘、断层、不整合面、尖灭等，为后续的地震解释和储层预测工作提供了可靠的地震资料，大大提高油气勘探开发的成功率。为达到衰减地震噪声，并增强地质结构特征的扩散滤波效果，寻求扩散张量是该方法最关键的步骤。
步骤A800，基于所述整形后波形数据，构建张量扩散模型：
Figure 787027DEST_PATH_IMAGE062
其中，
Figure 916657DEST_PATH_IMAGE063
表示扩散时间，
Figure 199871DEST_PATH_IMAGE064
表示散度算子，D表示扩散张量，U表示扩散滤波结果，
Figure 146967DEST_PATH_IMAGE065
表示
Figure 148421DEST_PATH_IMAGE063
=0时的扩散滤波结果，
Figure 81742DEST_PATH_IMAGE066
表示
Figure 78517DEST_PATH_IMAGE067
时刻的整形后波形数据，作为张量扩散模型 的初始条件，
Figure 71880DEST_PATH_IMAGE068
表示扩散滤波结果的梯度；
基于所述张量扩散模型构建梯度结构张量：
Figure 560631DEST_PATH_IMAGE069
其中，U表示扩散滤波结果，
Figure 32063DEST_PATH_IMAGE070
表示梯度向量张量积；
Figure 414503DEST_PATH_IMAGE071
表示尺度为
Figure 578768DEST_PATH_IMAGE072
的高斯函数：
Figure 289235DEST_PATH_IMAGE073
其中，r表示计算半径；
结构张量的特征向量为：
Figure 564359DEST_PATH_IMAGE074
其中，
Figure 801305DEST_PATH_IMAGE075
、
Figure 870892DEST_PATH_IMAGE076
和
Figure 334235DEST_PATH_IMAGE077
表示为梯度结构张量的3个特征向量，可视为局部正交坐标系，
Figure 275033DEST_PATH_IMAGE075
指向地震信号的梯度方向，
Figure 507431DEST_PATH_IMAGE076
和
Figure 482341DEST_PATH_IMAGE077
组成的平面平行于地震信号的局部结构特征，
Figure 557613DEST_PATH_IMAGE078
、
Figure 174539DEST_PATH_IMAGE079
和
Figure 261444DEST_PATH_IMAGE080
为分别与
Figure 797467DEST_PATH_IMAGE075
、
Figure 969823DEST_PATH_IMAGE076
和
Figure 124861DEST_PATH_IMAGE077
对应的三个特征值；
步骤A900，基于所述结构张量的特征向量分别计算线状结构置信度量、面状结构置信度量和扩散张量；
所述线状结构置信度量
Figure 331851DEST_PATH_IMAGE081
为：
Figure 38776DEST_PATH_IMAGE082
所述面状结构置信度量
Figure 698427DEST_PATH_IMAGE083
为：
Figure 391577DEST_PATH_IMAGE084
所述扩散张量D为：
Figure 577707DEST_PATH_IMAGE086
其中，
Figure 596479DEST_PATH_IMAGE087
、
Figure 743427DEST_PATH_IMAGE088
和
Figure 240267DEST_PATH_IMAGE089
表示扩散张量的三个非负特征值，他们分别表示扩散滤波器沿
Figure 280904DEST_PATH_IMAGE075
、
Figure 204998DEST_PATH_IMAGE076
和
Figure 839242DEST_PATH_IMAGE077
这三个特征方向的滤波强度；
扩散滤波算法保留了“串珠状”反射缝洞型储层地质特征，增强了地震数据对目标地质体的成像能力。同时起到了压制噪声、提高同相轴横向连续性与地震信号的信噪比的效果。
步骤A1000，重复步骤A100-步骤A900，直至达到预设的迭代次数，获得扩散滤波结果，即为提频降噪后的地震资料。
深层碳酸盐岩储层与上覆地层的岩性差异大，根据这一地质规律将碳酸盐岩岩层顶面作为标志层。
以提频降噪后的地震资料为基础，对标志层位地震波反射同相轴进行追踪显示，获得标志层层位解释连续面。
在本实施例中，所述时深转化关系，其获取方法为：
基于每个已知井位的测井资料中的声波时差曲线和密度曲线做乘积运算获取波阻抗曲线，进而计算反射系数曲线；
以目的层段地震主频为依据构建雷克子波，将雷克子波与反射系数曲线褶积计算后，得到合成地震记录；
根据原始地震资料，通过地震波激发装置与接收装置获取原始地震波反射信号数据，并根据原始地震波反射信号数据的波形获取目的层位标志层的等时三维展布模型；
将每个钻井井位井眼处的标志层深度数据与标志层三维展布模型对应，计算合成地震记录与井旁地震道提频降噪后的地震资料的相关性，当波形相关性高于第一相关阈值时，第一相关阈值优选85%，井震标定完成，最终得出测井深度与地震反射波双程旅行时之间的时深转化关系：
Figure 264407DEST_PATH_IMAGE090
其中，
Figure 34917DEST_PATH_IMAGE091
表示声波测井标志层深度对应的地震资料的双程旅行时；
Figure 395491DEST_PATH_IMAGE092
为声波时差；
Figure 517031DEST_PATH_IMAGE093
为测井曲线数据采样间隔；
Figure 480307DEST_PATH_IMAGE094
为地震波双程旅行时。
本实施例中还包括将提频后的地震资料进行多尺度分解获取地震资料多尺度Bimf分量，并用以生成精确时深转化关系的步骤，具体包括：
将提频降噪后的地震资料沿“道集”方向切片，获得若干个叠前道集矩阵；
Figure 370903DEST_PATH_IMAGE095
。
通过EMD分解的方式获取Bimf分量：
设置初始的叠前道集矩阵：
Figure 636799DEST_PATH_IMAGE096
其中，
Figure 101760DEST_PATH_IMAGE097
为矩阵行坐标，
Figure 9673DEST_PATH_IMAGE098
为矩阵列坐标，
Figure 754775DEST_PATH_IMAGE099
为叠后三维地震资料数据的值，
Figure 457152DEST_PATH_IMAGE100
为矩 阵元素值；
设置变量
Figure 412339DEST_PATH_IMAGE101
为分解的Bimf层数，
Figure 327205DEST_PATH_IMAGE101
的初始值为1，定义初始化矩阵变量
Figure 192393DEST_PATH_IMAGE102
，获得变量矩 阵
Figure 190305DEST_PATH_IMAGE103
为：
Figure 508153DEST_PATH_IMAGE104
计算关于
Figure 492290DEST_PATH_IMAGE103
的局部最大值矩阵和局部最小值矩阵：
在变量矩阵
Figure 211984DEST_PATH_IMAGE103
中以预设尺寸的观测矩阵进行观测，提取观测矩阵中的元素 值，并获取观测矩阵的最大值；所述预设尺寸可选用
Figure 380798DEST_PATH_IMAGE105
4的小矩阵进行观测；
滑动所述观测矩阵，直至观测矩阵中心遍历所有矩阵变量h，将测得的观测矩阵的 最大值赋值于对应观测矩阵的中心元素位置，获得局部最大值矩阵
Figure 185942DEST_PATH_IMAGE106
；
通过设置观测矩阵获取观测矩阵中的最小值，进而获得局部最小值矩阵
Figure 708191DEST_PATH_IMAGE107
；
通过样条插值的方法获取最大值包络面和最小值包络面；
插值公式为：
Figure 407025DEST_PATH_IMAGE108
其中，A、B和C表示插值过渡矩阵：
Figure 622106DEST_PATH_IMAGE109
Figure 648968DEST_PATH_IMAGE111
Figure 240486DEST_PATH_IMAGE112
其中，u和v为预设的大于0且小于1的值，
Figure 793827DEST_PATH_IMAGE113
为扩频模拟子波；
求得
Figure 179809DEST_PATH_IMAGE114
为最大值包络面，将过渡公式B中的局部最大值矩阵
Figure 959546DEST_PATH_IMAGE106
替换为局 部最小值矩阵
Figure 89176DEST_PATH_IMAGE107
，获得
Figure 231445DEST_PATH_IMAGE115
为最小值包络面；
计算最大值包络面和最小值包络面的均值，获得中间过渡量：
Figure 53907DEST_PATH_IMAGE116
将变量矩阵
Figure 55361DEST_PATH_IMAGE103
减去中间过渡量，获得地震资料对应的第j层Bimf分量：
Figure 847737DEST_PATH_IMAGE117
获取余量
Figure 251036DEST_PATH_IMAGE118
作为计算下一层Bimf分量的输入量：
Figure 244400DEST_PATH_IMAGE119
Figure 467571DEST_PATH_IMAGE120
表示第l层的余量，
Figure 66567DEST_PATH_IMAGE118
表示第
Figure 324373DEST_PATH_IMAGE121
层的余量，第
Figure 223059DEST_PATH_IMAGE121
层的余量作为第l层Bimf分 量的输入量；
迭代计算Bimf分量，直至满足终止条件：
Figure 323739DEST_PATH_IMAGE122
获得地震资料多尺度Bimf分量；其中，r为预设的终止阈值；r通常取值为0.2~0.3，能够保证Bimf的数目和质量，保证其能够更好地反映波形的细节。
获取的Bimf分量的具体用法为：
将所述地震资料Bimf分量以尺度由大到小的顺序，逐步添加Bimf分量获得更新后的地震资料；
每添加一个Bimf分量获得更新后的有效地震资料计算一次更新后的合成地震记录与有效地震资料的第二相关性；
随着有效地震资料逐步添加Bimf分量，第二相关性起初呈现上升趋势，当第二相关性出现下降趋势时，取第二相关性峰值时的更新后的合成地震记录与更新后的地震，计算精确时深转化关系。通过削减小尺度Bimf分量的方式进一步降低了噪音的影响，进一步突出了断溶体纵向发育的特征。
S300，基于所述提频降噪后的地震资料，构建目标层位的等时地层格架模型，并计算断溶体几何结构和空间分布；等时地层格架模型指示地震波反射界面随地震波到时的变化关系；所述断溶体几何结构和空间分布如图2所示；
在本实施例中，所述构建目标层位的等时地层格架模型，其获取方法为：
基于所述提频降噪后的地震资料，对标志层位地震波反射同相轴进行追踪显示，获得标志层层位解释连续面；
基于所述标志层层位解释连续面，预设时窗、切片方向和最大深度，构建目标层位的等时地层格架模型。构建的等时地层格架模型如图3所示；
S400，基于所述测井资料和提频降噪后的地震资料，在所述等时地层格架模型的提供的地层趋势信息基础上，进行波形指示地震波阻抗反演获得波阻抗反演数据体；
在本实施例中，所述波阻抗反演数据体，其获取方式包括：
基于所述提频降噪后的地震资料，计算待判别地震道波形与已知井的合成地震记录的波形相关性，根据波形相关性最高的井对应的波阻抗曲线建立初始模型；在本实施例中，建立初始模型中选用的波形相关性最高的井可以为预设个数的样本井，也可以为最优样本数个样本井；
所述最优样本数，其获得方法为：
任选一样本井作为参考目标井，设定初始的样本数参量为1；
根据波形相似性原则选取数量为样本数参量的样本井标准化测井曲线的特征参数数据与参考目标井标准化测井曲线的特征参数数据进行相关性分析获得样本数参量-参考目标井曲线特征参数的相关性值；
逐1增加样本数参量，重复获取各样本数参量对应的样本数参量-参考目标井曲线特征参数的相关性值，将所有样本数参量-参考目标井曲线特征参数的相关性值连接，获得参考井曲线特征参数相关性随样本数参量变化的相关性曲线；
选取另一样本井作为参考目标井，获得多个参考井的曲线特征参数相关性随样本数参量变化的相关性曲线，将所有参考井的特征参数的相关性随样本数参量变化的相关性曲线拟合为整体相关性曲线，选取所述整体相关性曲线中相关性随样本参数增加而升高最终保持平稳处的拐点，确定最优样本数参量；
基于所述高精度的三维地震振幅数据体和等时地层格架模型，计算待测点位与样本井位的波形相关性，将所述波形相关性由大到小排序，选取最优样本数参量条地震波形相关性最高的样本井的对储层敏感的特征参数数据；基于所述相关性最高的样本井的地震波形特征数据所对应的样本井，通过井间特征参数插值方式构建初始模型；
利用白噪声满足高斯分布的规律，将测井资料中的波阻抗曲线表示为：
Figure 598863DEST_PATH_IMAGE123
其中，
Figure 711175DEST_PATH_IMAGE002
表示测井波阻抗曲线，
Figure 46342DEST_PATH_IMAGE003
表示待求解的地下地层实际波阻抗值，N表示 随机噪声；
根据中心极限定理，
Figure 368738DEST_PATH_IMAGE004
也满足高斯分布，确定初始目标函数为：
Figure 181974DEST_PATH_IMAGE124
其中，
Figure 414372DEST_PATH_IMAGE006
表示与后验信息有关的函数，
Figure 779494DEST_PATH_IMAGE007
表示基于最优样本数选取样本井对样本 井的波阻抗曲线，进行匹配滤波后，求得后验概率统计分布密度，进而计算得到的波阻抗期 望值，
Figure 730133DEST_PATH_IMAGE008
表示白噪声的协方差；
基于所述初始目标函数，通过最大后验估计，在目标函数中引入先验信息，获得稳定的目标函数为：
Figure 347059DEST_PATH_IMAGE125
其中，
Figure 433963DEST_PATH_IMAGE010
表示待模拟的特征参数，
Figure 969987DEST_PATH_IMAGE011
表示与地质和测井资料先验信息有关的函 数，
Figure 142342DEST_PATH_IMAGE012
表示用于协调
Figure 297380DEST_PATH_IMAGE013
和
Figure 363425DEST_PATH_IMAGE011
之间的相互影响的平滑参数；
以所述稳定的目标函数作为初始模型的输入，通过马尔科夫链蒙特卡罗方法MCMC和Metropolis-Hastings抽样准则对后验概率分布抽样，不断优化初始模型的参数，选取目标函数取最大值时的解作为随机实现，取多次随机实现的均值作为期望值输出，将所述期望值输出作为波阻抗反演数据体。获得的波阻抗反演数据体如图4所示；
在本实施例中，根据稳定的目标函数获取波阻抗反演数据体，具体包括：
设M为目标空间，n为总样本数，m为马尔科夫链趋于达到平稳时的样本数；
预设一条马尔科夫链，使马尔科夫链收敛至达到平稳；
由M中的某一点
Figure 211295DEST_PATH_IMAGE014
出发，通过马尔科夫链进行抽样模拟，产生点序列：
Figure 870947DEST_PATH_IMAGE015
；
函数
Figure 688730DEST_PATH_IMAGE016
的期望估计为：
Figure 750227DEST_PATH_IMAGE017
其中，n为总样本数，m表示马尔科夫链达到平稳时的样本数，k表示累加参量；
选取一转移函数
Figure 768999DEST_PATH_IMAGE018
和初始值
Figure 650367DEST_PATH_IMAGE019
，若第i次迭代开始时的参数值为
Figure 537421DEST_PATH_IMAGE020
， 则第h次迭代过程为：
从
Figure 453424DEST_PATH_IMAGE018
中抽取一个备选值
Figure 643097DEST_PATH_IMAGE021
，计算备选值
Figure 277340DEST_PATH_IMAGE021
的接受概率
Figure 433997DEST_PATH_IMAGE022
：
Figure 204507DEST_PATH_IMAGE023
以
Figure 565081DEST_PATH_IMAGE024
，置
Figure 545675DEST_PATH_IMAGE025
；以概率
Figure 649897DEST_PATH_IMAGE026
，置
Figure 274914DEST_PATH_IMAGE027
；
不断扰动所述初始模型的参数，重复选取M中的某一点
Figure 931023DEST_PATH_IMAGE014
求解后验样本
Figure 274280DEST_PATH_IMAGE028
，直至 达到预设的迭代次数
Figure 916614DEST_PATH_IMAGE029
，获得后验样本
Figure 661716DEST_PATH_IMAGE030
，进而计算后验分布的各阶矩 阵获得期望输出值，将期望值输出作为波阻抗反演数据体。
S500，基于野外地质露头资料，确定目标层位的地层密度序列；具体为：利用野外地质露头资料，确定目标层位的岩性序列，并对岩性序列进行物理性质（密度）分析，生成随深度变化的地层密度序列。
S600，基于所述时深转化关系，确认所述地层密度序列与等时地层格架模型的纵向对应关系，进而根据测井资料构建地层密度模型和地层声波时差模型；
在本实施例中，所述地层密度模型和地层声波时差模型，其获取方法为：
基于所述时深转化关系，建立地层密度序列与地震波到时Tsm的关系；依据地震波到时匹配地层密度序列与等时地层格架模型的纵向对应关系，得到等时地层格架模型内部的密度分布，即地层密度模型；
基于所述时深转化关系，建立测井资料中声波时差测井曲线与地震波到时Tsm的关系；依据地震波到时匹配地层声波时差与等时地层格架模型的纵向对应关系，得到等时地层格架模型内部的声波时差分布，即地层声波时差模型。
S700，基于所述地层密度模型和地层声波时差模型，计算背景地层波阻抗模型；
在本实施例中，所述背景地层波阻抗模型为：
Figure 488726DEST_PATH_IMAGE126
其中，IMP为波阻抗值，DEN为密度值，vp为纵波波速。
S800，将所述背景地层波阻抗模型与波阻抗反演数据体进行做差运算，获得波阻抗异常数据体；
S900，基于断溶体几何结构和空间分布，在所述波阻抗异常数据体中仅保留断溶体几何轮廓对应区域的数据，获得断溶体几何结构-波阻抗模型；
在本实施例中，所述断溶体几何结构-波阻抗模型，其获方法为：
基于所述提频降噪后的地震资料计算地震信号方差属性值数据体，剔除地震信号方差值小于预设的第二阈值的区域得到断溶体轮廓；剔除所述断溶体波阻抗反演数据体在断溶体轮廓以外区域的数据，获得断溶体几何结构-波阻抗模型。
在本实施例中，所述地震信号方差属性值数据体，其计算方法为：
设提频降噪后的地震资料中的各采样点数据为
Figure 319279DEST_PATH_IMAGE127
，其中
Figure 499725DEST_PATH_IMAGE128
表示地震测网线号，
Figure 489546DEST_PATH_IMAGE129
代表地震测网道号，k代表1ms采样的地震记录采样点序号；
计算预设采样区域的采样点数据均方差：
Figure 362824DEST_PATH_IMAGE130
将采样区域进行纵向和横向平移，遍历计算所有采样区域的数据均方差，获得三维方差属性体；
对所述三维方差属性体进行切片，获取方差属性数据在平面上的分布特征，将方差属性值低于总体能量均值的区域作为沉积稳定区。
本实施例中，对断溶体波阻抗数据体的能量进行统计；因为这些区域代表沉积稳定区的显示，较高特征值代表由断控溶蚀作用引起的地震资料扰动。因此对特征值低于0.62的数据进行切除。由于断溶体波阻抗数据体低能量区域剔除后还无法确定断裂带边界，需要通过地震波形能量属性进一步刻画。由提频降噪后的地震资料计算地震信号能量属性值数据体，剔除地震信号能量值小于0.8的区域得到断溶体轮廓。
S1000，对比所述断溶体几何结构-波阻抗模型与已知井位中的断溶体井位的测井解释结果，划定孔洞储层特征值区间、断裂带特征值区间和围岩特征值区间，获得碳酸盐岩断控储层解释模型。碳酸盐岩断控储层解释模型如图5所示；
比对测井解释结果与断溶体结构-特征值模型的特征值能量数据，划定特征值大于0.82的区域为裂缝储层，划定特征值处于0.63~0.82的区域为孔洞储层，划定特征值处于0.31~0.63的区域为断裂带，划定特征值小于0.31的区域为围岩，以此作为最终的碳酸盐岩断控储层解释模型。
本发明第二实施例的碳酸盐岩断控储层的波阻抗提取与结构表征系统，包括：资料获取模块、地震提频和井震标定模块、地层结构获取模块、波阻抗反演数据体获取模块、地层密度序列确定模块、地层声波时差模型获取模块、背景地层波阻抗模型获取模块、波阻抗异常数据体获取模块、断溶体几何结构-波阻抗模型和碳酸盐岩断控储层解释模型获取模块；
所述资料获取模块，配置为获取已知井位的原始地震资料和测井资料；
所述地震提频和井震标定模块，配置为基于所述原始地震资料和测井资料，通过谱整形和扩散滤波的方式，获得提频降噪后的地震资料，并进行井震标定获得时深转化关系；
所述地层结构获取模块，配置为基于所述提频降噪后的地震资料，构建目标层位的等时地层格架模型，并计算断溶体几何结构和空间分布；
所述波阻抗反演数据体获取模块，配置为基于所述测井资料和提频降噪后的地震资料，在所述等时地层格架模型的提供的地层趋势信息基础上，进行波形指示地震波阻抗反演获得波阻抗反演数据体；
所述地层密度序列确定模块，配置为基于野外地质露头资料，确定目标层位的地层密度序列；
所述地层声波时差模型获取模块，配置为基于所述时深转化关系，确认所述地层密度序列与等时地层格架模型的纵向对应关系，进而根据测井资料构建地层密度模型和地层声波时差模型；
所述背景地层波阻抗模型获取模块，配置为基于所述地层密度模型和地层声波时差模型，计算背景地层波阻抗模型；
所述波阻抗异常数据体获取模块，配置为将所述背景地层波阻抗模型与波阻抗反演数据体进行做差运算，获得波阻抗异常数据体；
所述断溶体几何结构-波阻抗模型，配置为在所述波阻抗异常数据体中仅保留断溶体几何结构对应区域的数据，获得断溶体几何结构-波阻抗模型；
所述碳酸盐岩断控储层解释模型获取模块，配置为对比所述断溶体几何结构-波阻抗模型与已知井位中的断溶体井位的测井解释结果，划定孔洞储层特征值区间、断裂带特征值区间和围岩特征值区间，获得碳酸盐岩断控储层解释模型。
所属技术领域的技术人员可以清楚地了解到，为描述的方便和简洁，上述描述的存储装置、处理装置的具体工作过程及有关说明，可以参考前述方法实施例中的对应过程，在此不再赘述。
术语“第一”、 “第二”等是用于区别类似的对象，而不是用于描述或表示特定的顺序或先后次序。
术语“包括”或者任何其它类似用语旨在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备/装置不仅包括那些要素，而且还包括没有明确列出的其它要素，或者还包括这些过程、方法、物品或者设备/装置所固有的要素。
至此，已经结合附图所示的优选实施方式描述了本发明的技术方案，但是，本领域技术人员容易理解的是，本发明的保护范围显然不局限于这些具体实施方式。在不偏离本发明的原理的前提下，本领域技术人员可以对相关技术特征做出等同的更改或替换，这些更改或替换之后的技术方案都将落入本发明的保护范围之内。

Claims (10)

1.一种碳酸盐岩断控储层的波阻抗提取与结构表征方法，其特征在于，所述方法包括：
获取已知井位的原始地震资料和测井资料；
基于所述原始地震资料和测井资料，通过谱整形和扩散滤波的方式，获得提频降噪后的地震资料，并进行井震标定获得时深转化关系；
基于所述提频降噪后的地震资料，构建目标层位的等时地层格架模型，并计算断溶体几何结构和空间分布；
基于所述测井资料和提频降噪后的地震资料，在所述等时地层格架模型的提供的地层趋势信息基础上，进行波形指示地震波阻抗反演获得波阻抗反演数据体；
基于野外地质露头资料，确定目标层位的地层密度序列；
基于所述时深转化关系，确认所述地层密度序列与等时地层格架模型的纵向对应关系，进而根据测井资料构建地层密度模型和地层声波时差模型；
基于所述地层密度模型和地层声波时差模型，计算背景地层波阻抗模型；
将所述背景地层波阻抗模型与波阻抗反演数据体进行做差运算，获得波阻抗异常数据体；
基于断溶体几何结构和空间分布，在所述波阻抗异常数据体中仅保留断溶体轮廓内的波阻抗异常数据体中的数据点，获得断溶体几何结构-波阻抗模型；
对比所述断溶体几何结构-波阻抗模型与已知井位中的断溶体井位的测井解释结果，划定孔洞储层特征值区间、断裂带特征值区间和围岩特征值区间，获得碳酸盐岩断控储层解释模型。
2.根据权利要求1所述的碳酸盐岩断控储层的波阻抗提取与结构表征方法，其特征在于，所述提频降噪后的地震资料，其获取方法为：
步骤A100，将原始地震资料的地震记录褶积模型在频率域的表示为：
Figure DEST_PATH_IMAGE001
其中，
Figure 168388DEST_PATH_IMAGE002
表示傅氏变换后的地震记录频谱，
Figure DEST_PATH_IMAGE003
表示傅氏变换后的地震子波频谱，
Figure 679266DEST_PATH_IMAGE004
表示傅氏变换后的反射系数频谱，
Figure DEST_PATH_IMAGE005
表示角频率；
步骤A200，将所述地震记录褶积模型在频率域的表示转化为地震记录线性系统频域的表示：
Figure 618403DEST_PATH_IMAGE006
其中，
Figure DEST_PATH_IMAGE007
为地震记录线性系统频域表示，
Figure 269833DEST_PATH_IMAGE008
为地震子波线性系统频域表示，
Figure DEST_PATH_IMAGE009
为反射系数线性系统频域表示；
步骤A300，将所述地震记录线性系统频域的表示进行反傅氏变换，获得地震记录频谱的复赛谱序列：
Figure 610816DEST_PATH_IMAGE010
其中，
Figure DEST_PATH_IMAGE011
表示地震记录频谱的复赛谱序列，
Figure 338600DEST_PATH_IMAGE012
表示地震子波频谱的复赛谱序列，
Figure DEST_PATH_IMAGE013
表示反射系数频谱的复赛谱序列，
Figure 781345DEST_PATH_IMAGE014
表示地震波到时；
步骤A400，通过低通滤波器将复赛谱中的子波复赛谱序列和反射系数复赛谱序列进行分离，提取子波振幅；
步骤A500，基于所述子波振幅，通过最小二乘法模拟地震子波振幅谱：
Figure DEST_PATH_IMAGE015
其中，
Figure 456040DEST_PATH_IMAGE016
表示大于等于0的常数，
Figure DEST_PATH_IMAGE017
表示地震子波频谱，由
Figure 166376DEST_PATH_IMAGE018
经傅氏变换而来，
Figure DEST_PATH_IMAGE019
和
Figure 533903DEST_PATH_IMAGE020
为待求的关于频率
Figure DEST_PATH_IMAGE021
的多项式用于拟合地震子波振幅谱；
步骤A600，基于所述模拟地震子波振幅谱，获得子波最大相位分量和最小相位分量；
设子波
Figure 198365DEST_PATH_IMAGE022
的最大相位分量为
Figure DEST_PATH_IMAGE023
、最小相位分量为
Figure 411172DEST_PATH_IMAGE024
，则子波
Figure DEST_PATH_IMAGE025
为：
Figure 38331DEST_PATH_IMAGE026
振幅谱的复赛谱中表示为：
Figure DEST_PATH_IMAGE027
其中，振幅谱的复赛谱
Figure 45601DEST_PATH_IMAGE028
在复赛谱的正、负轴上对称显示，
Figure DEST_PATH_IMAGE029
为地震子波最小相位分量
Figure 508944DEST_PATH_IMAGE030
所对应的最大相位函数的复赛谱，
Figure DEST_PATH_IMAGE031
为地震子波最大相位分量
Figure 276174DEST_PATH_IMAGE032
所对应的最小相位函数的复赛谱；
步骤A700，基于所述振幅谱的复赛谱确定一组具有相同振幅谱的混合相位子波集合，调整俞氏子波参数，在保证地震子波主频完整的前提下，提升有效带宽至预设的带宽阈值，获得整形后波形数据；
步骤A800，基于所述整形后波形数据，构建张量扩散模型：
Figure DEST_PATH_IMAGE033
其中，
Figure 711834DEST_PATH_IMAGE034
表示扩散时间，
Figure DEST_PATH_IMAGE035
表示散度算子，D表示扩散张量，U表示扩散滤波结果，
Figure 873694DEST_PATH_IMAGE036
表示
Figure DEST_PATH_IMAGE037
时的扩散滤波结果，
Figure 981590DEST_PATH_IMAGE038
表示
Figure 332937DEST_PATH_IMAGE037
时刻的整形后波形数据，作为张量扩散模型的初始条件，
Figure DEST_PATH_IMAGE039
表示扩散滤波结果的梯度；
基于所述张量扩散模型构建梯度结构张量：
Figure 872371DEST_PATH_IMAGE040
其中，
Figure DEST_PATH_IMAGE041
表示扩散滤波结果，
Figure 752603DEST_PATH_IMAGE042
表示梯度向量张量积；
Figure DEST_PATH_IMAGE043
表示尺度为
Figure 890672DEST_PATH_IMAGE044
的高斯函数：
Figure DEST_PATH_IMAGE045
其中，
Figure 248972DEST_PATH_IMAGE046
表示计算半径；
结构张量的特征向量为：
Figure DEST_PATH_IMAGE047
其中，
Figure 908492DEST_PATH_IMAGE048
、
Figure DEST_PATH_IMAGE049
和
Figure 756362DEST_PATH_IMAGE050
表示为梯度结构张量的3个特征向量，可视为局部正交坐标系，
Figure 150435DEST_PATH_IMAGE048
指向地震信号的梯度方向，
Figure DEST_PATH_IMAGE051
和
Figure 63158DEST_PATH_IMAGE052
组成的平面平行于地震信号的局部结构特征，
Figure DEST_PATH_IMAGE053
、
Figure 249289DEST_PATH_IMAGE054
和
Figure DEST_PATH_IMAGE055
为分别与
Figure 471323DEST_PATH_IMAGE056
、
Figure DEST_PATH_IMAGE057
和
Figure 572265DEST_PATH_IMAGE052
对应的三个特征值；
步骤A900，基于所述结构张量的特征向量分别计算线状结构置信度量、面状结构置信度量和扩散张量；
所述线状结构置信度量
Figure 69105DEST_PATH_IMAGE058
为：
Figure DEST_PATH_IMAGE059
所述面状结构置信度量
Figure 188371DEST_PATH_IMAGE060
为：
Figure DEST_PATH_IMAGE061
所述扩散张量D为：
Figure 299415DEST_PATH_IMAGE062
其中，
Figure DEST_PATH_IMAGE063
、
Figure 136921DEST_PATH_IMAGE064
和
Figure DEST_PATH_IMAGE065
表示扩散张量的三个非负特征值，他们分别表示扩散滤波器沿
Figure 657027DEST_PATH_IMAGE066
、
Figure DEST_PATH_IMAGE067
和
Figure 630799DEST_PATH_IMAGE068
这三个特征方向的滤波强度；
步骤A1000，重复步骤A100-步骤A900，直至达到预设的迭代次数，获得扩散滤波结果，即为提频降噪后的地震资料。
3.根据权利要求2所述的碳酸盐岩断控储层的波阻抗提取与结构表征方法，其特征在于，所述构建目标层位的等时地层格架模型，其获取方法为：
基于所述提频降噪后的地震资料，对标志层位地震波反射同相轴进行追踪显示，获得标志层层位解释连续面；
基于所述标志层层位解释连续面，预设时窗、切片方向和最大深度，构建目标层位的等时地层格架模型。
4.根据权利要求1所述的碳酸盐岩断控储层的波阻抗提取与结构表征方法，其特征在于，所述时深转化关系，其获取方法为：
基于每个已知井位的测井资料中的声波时差曲线和密度曲线做乘积运算获取波阻抗曲线，进而计算反射系数曲线；
以目的层段地震主频为依据构建雷克子波，将雷克子波与反射系数曲线褶积计算后，得到合成地震记录；
将每个钻井井位井眼处的标志层深度数据与标志层三维展布模型对应，计算合成地震记录与井旁地震道提频降噪后的地震资料的相关性，当波形相关性高于第一相关阈值时，井震标定完成，最终得出测井深度与地震反射波双程旅行时之间的时深转化关系：
Figure DEST_PATH_IMAGE069
其中，
Figure 443903DEST_PATH_IMAGE070
表示声波测井标志层深度对应的地震资料的双程旅行时；
Figure DEST_PATH_IMAGE071
为声波时差；
Figure 503126DEST_PATH_IMAGE072
为测井曲线数据采样间隔；
Figure DEST_PATH_IMAGE073
为地震波双程旅行时。
5.根据权利要求4所述的碳酸盐岩断控储层的波阻抗提取与结构表征方法，其特征在于，所述波阻抗反演数据体，其获取方式包括：
基于所述提频降噪后的地震资料，计算待判别地震道波形与已知井的合成地震记录的波形相关性，根据波形相关性高于预设阈值的井对应的波阻抗曲线建立初始模型；
利用白噪声满足高斯分布的规律，将测井资料中的波阻抗曲线表示为：
Figure 561343DEST_PATH_IMAGE074
其中，
Figure DEST_PATH_IMAGE075
表示测井波阻抗曲线，
Figure 655201DEST_PATH_IMAGE076
表示待求解的地下地层实际波阻抗值，
Figure DEST_PATH_IMAGE077
表示随机噪声；
根据中心极限定理，
Figure 842469DEST_PATH_IMAGE078
也满足高斯分布，确定初始目标函数为：
Figure DEST_PATH_IMAGE079
其中，
Figure 405299DEST_PATH_IMAGE080
表示与后验信息有关的函数，
Figure DEST_PATH_IMAGE081
表示基于样本井的波阻抗曲线，进行匹配滤波后，求得后验概率统计分布密度，进而计算得到的波阻抗期望值，
Figure 719737DEST_PATH_IMAGE082
表示白噪声的协方差；
基于所述初始目标函数，通过最大后验估计，在目标函数中引入先验信息，获得稳定的目标函数为：
Figure DEST_PATH_IMAGE083
其中，
Figure 386211DEST_PATH_IMAGE084
表示待模拟的特征参数，
Figure DEST_PATH_IMAGE085
表示与地质和测井资料先验信息有关的函数，
Figure 308161DEST_PATH_IMAGE086
表示用于协调
Figure DEST_PATH_IMAGE087
和
Figure 545239DEST_PATH_IMAGE088
之间的相互影响的平滑参数；
以所述稳定的目标函数作为初始模型的输入，通过马尔科夫链蒙特卡罗方法MCMC和Metropolis-Hastings抽样准则对后验概率分布抽样，不断优化初始模型的参数，选取目标函数取最大值时的解作为随机实现，取多次随机实现的均值作为期望值输出，将所述期望值输出作为波阻抗反演数据体。
6.根据权利要求5所述的碳酸盐岩断控储层的波阻抗提取与结构表征方法，其特征在于，根据稳定的目标函数获取波阻抗反演数据体，具体包括：
设M为目标空间，n为总样本数，m为马尔科夫链达到平稳时的样本数；
预设一条马尔科夫链，使马尔科夫链收敛至达到平稳；
由M中的某一点
Figure DEST_PATH_IMAGE089
出发，通过马尔科夫链进行抽样模拟，产生点序列：
Figure 240531DEST_PATH_IMAGE090
；
函数
Figure DEST_PATH_IMAGE091
的期望估计为：
Figure 43402DEST_PATH_IMAGE092
其中，n为总样本数，m表示马尔科夫链达到平稳时的样本数，k表示累加参量；
选取一转移函数
Figure DEST_PATH_IMAGE093
和初始值
Figure 870675DEST_PATH_IMAGE094
，若第h次迭代开始时的参数值为
Figure DEST_PATH_IMAGE095
，则第h次迭代过程为：
从
Figure 313158DEST_PATH_IMAGE096
中抽取一个备选值
Figure DEST_PATH_IMAGE097
，计算备选值
Figure 766136DEST_PATH_IMAGE097
的接受概率
Figure 220251DEST_PATH_IMAGE098
：
Figure DEST_PATH_IMAGE099
以
Figure 206706DEST_PATH_IMAGE100
，置
Figure DEST_PATH_IMAGE101
；以概率
Figure 949534DEST_PATH_IMAGE102
，置
Figure DEST_PATH_IMAGE103
；
不断扰动所述初始模型的参数，重复选取M中的某一点
Figure 189892DEST_PATH_IMAGE104
求解后验样本
Figure DEST_PATH_IMAGE105
，直至达到预设的迭代次数
Figure 967355DEST_PATH_IMAGE106
，获得后验样本
Figure DEST_PATH_IMAGE107
，进而计算后验分布的各阶矩阵获得期望输出值，将期望值输出作为波阻抗反演数据体。
7.根据权利要求1所述的碳酸盐岩断控储层的波阻抗提取与结构表征方法，其特征在于，所述地层密度模型和地层声波时差模型，其获取方法为：
基于所述时深转化关系，建立地层密度序列与地震波到时的关系；依据地震波到时匹配地层密度序列与等时地层格架模型的纵向对应关系，得到等时地层格架模型内部的密度分布，即地层密度模型；
基于所述时深转化关系，建立测井资料中声波时差测井曲线与地震波到时的关系；依据地震波到时匹配地层声波时差与等时地层格架模型的纵向对应关系，得到等时地层格架模型内部的声波时差分布，即地层声波时差模型。
8.根据权利要求7所述的碳酸盐岩断控储层的波阻抗提取与结构表征方法，其特征在于，所述背景地层波阻抗模型为：
Figure 933168DEST_PATH_IMAGE108
其中，IMP为波阻抗值，DEN为密度值，vp为纵波波速。
9.根据权利要求1所述的碳酸盐岩断控储层的波阻抗提取与结构表征方法，其特征在于，所述断溶体几何结构-波阻抗模型，其获方法为：
基于所述提频降噪后的地震资料计算地震信号方差属性值数据体，剔除地震信号方差值小于预设的第二阈值的区域得到断溶体轮廓；保留断溶体轮廓内的波阻抗异常数据体中的数据点，获得断溶体几何结构-波阻抗模型。
10.一种碳酸盐岩断控储层的波阻抗提取与结构表征系统，其特征在于，所述提取与结构表征系统包括：资料获取模块、地震提频和井震标定模块、地层结构获取模块、波阻抗反演数据体获取模块、地层密度序列确定模块、地层声波时差模型获取模块、背景地层波阻抗模型获取模块、波阻抗异常数据体获取模块、断溶体几何结构-波阻抗模型和碳酸盐岩断控储层解释模型获取模块；
所述资料获取模块，配置为获取已知井位的原始地震资料和测井资料；
所述地震提频和井震标定模块，配置为基于所述原始地震资料和测井资料，通过谱整形和扩散滤波的方式，获得提频降噪后的地震资料，并进行井震标定获得时深转化关系；
所述地层结构获取模块，配置为基于所述提频降噪后的地震资料，构建目标层位的等时地层格架模型，并计算断溶体几何结构和空间分布；
所述波阻抗反演数据体获取模块，配置为基于所述测井资料和提频降噪后的地震资料，在所述等时地层格架模型提供的地层趋势信息基础上，进行波形指示地震波阻抗反演获得波阻抗反演数据体；
所述地层密度序列确定模块，配置为基于野外地质露头资料，确定目标层位的地层密度序列；
所述地层声波时差模型获取模块，配置为基于所述时深转化关系，确认所述地层密度序列与等时地层格架模型的纵向对应关系，进而根据测井资料构建地层密度模型和地层声波时差模型；
所述背景地层波阻抗模型获取模块，配置为基于所述地层密度模型和地层声波时差模型，计算背景地层波阻抗模型；
所述波阻抗异常数据体获取模块，配置为将所述背景地层波阻抗模型与波阻抗反演数据体进行做差运算，获得波阻抗异常数据体；
所述断溶体几何结构-波阻抗模型，配置为基于断溶体几何结构和空间分布，在所述波阻抗异常数据体中仅保留断溶体轮廓内的波阻抗异常数据体中的数据点，获得断溶体几何结构-波阻抗模型；
所述碳酸盐岩断控储层解释模型获取模块，配置为对比所述断溶体几何结构-波阻抗模型与已知井位中的断溶体井位的测井解释结果，划定孔洞储层特征值区间、断裂带特征值区间和围岩特征值区间，获得碳酸盐岩断控储层解释模型。
CN202210919668.0A 2022-08-02 2022-08-02 碳酸盐岩断控储层的波阻抗提取与结构表征方法和系统 Active CN114994758B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202210919668.0A CN114994758B (zh) 2022-08-02 2022-08-02 碳酸盐岩断控储层的波阻抗提取与结构表征方法和系统

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202210919668.0A CN114994758B (zh) 2022-08-02 2022-08-02 碳酸盐岩断控储层的波阻抗提取与结构表征方法和系统

Publications (2)

Publication Number Publication Date
CN114994758A CN114994758A (zh) 2022-09-02
CN114994758B true CN114994758B (zh) 2022-10-28

Family

ID=83021165

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202210919668.0A Active CN114994758B (zh) 2022-08-02 2022-08-02 碳酸盐岩断控储层的波阻抗提取与结构表征方法和系统

Country Status (1)

Country Link
CN (1) CN114994758B (zh)

Families Citing this family (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116299701B (zh) * 2023-02-23 2024-02-13 北京阳光吉澳能源技术有限公司 断溶体储层识别方法、装置及电子设备
CN117094233B (zh) * 2023-10-19 2024-01-30 中国科学院地质与地球物理研究所 深地油气精准导航三维速度场高精度建模方法与系统

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101206264A (zh) * 2007-11-08 2008-06-25 符力耘 高分辨率非线性地震波阻抗反演方法
CN102650701A (zh) * 2012-04-27 2012-08-29 中国石油天然气股份有限公司 一种基于偶极子波进行储层波阻抗预测的方法
CN103439740A (zh) * 2013-07-24 2013-12-11 中国石油天然气股份有限公司 基于偶极地震子波多重积分的相对阻抗预测的方法及装置
CN108535775A (zh) * 2018-03-30 2018-09-14 中国石油大学(北京) 非平稳地震资料声波阻抗反演方法及装置
CN113759425A (zh) * 2021-09-13 2021-12-07 中国科学院地质与地球物理研究所 井震联合评价深层古岩溶储层充填特征的方法与系统

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US5583825A (en) * 1994-09-02 1996-12-10 Exxon Production Research Company Method for deriving reservoir lithology and fluid content from pre-stack inversion of seismic data
US8055449B2 (en) * 2007-10-08 2011-11-08 Westerngeco L.L.C. Determining fault transmissivity in a subterranean reservoir
CN103487835B (zh) * 2012-06-12 2016-03-23 中国石油化工股份有限公司 一种基于模型约束的多分辨率波阻抗反演方法
CN108680955B (zh) * 2017-12-12 2019-09-17 中国地质大学（武汉） 一种岩性识别方法及识别系统

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101206264A (zh) * 2007-11-08 2008-06-25 符力耘 高分辨率非线性地震波阻抗反演方法
CN102650701A (zh) * 2012-04-27 2012-08-29 中国石油天然气股份有限公司 一种基于偶极子波进行储层波阻抗预测的方法
CN103439740A (zh) * 2013-07-24 2013-12-11 中国石油天然气股份有限公司 基于偶极地震子波多重积分的相对阻抗预测的方法及装置
CN108535775A (zh) * 2018-03-30 2018-09-14 中国石油大学(北京) 非平稳地震资料声波阻抗反演方法及装置
CN113759425A (zh) * 2021-09-13 2021-12-07 中国科学院地质与地球物理研究所 井震联合评价深层古岩溶储层充填特征的方法与系统

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
基于井间地震资料的储层精细描述方法;李建华等;《石油地球物理勘探》;20080215;第43卷(第01期);第41-47页 *
塔河油田奥陶系碳酸盐岩缝洞型储层预测技术;韩革华等;《石油与天然气地质》;20061230;第27卷(第06期);第860-870页 *

Also Published As

Publication number Publication date
CN114994758A (zh) 2022-09-02

Similar Documents

Publication Publication Date Title
CN113759425B (zh) 井震联合评价深层古岩溶储层充填特征的方法与系统
CN113759424B (zh) 基于频谱分解和机器学习的岩溶储层充填分析方法和系统
CN109425896B (zh) 白云岩油气储层分布预测方法及装置
CN108802812B (zh) 一种井震融合的地层岩性反演方法
CN114994758B (zh) 碳酸盐岩断控储层的波阻抗提取与结构表征方法和系统
CN111596364B (zh) 基于高精度层序地层格架的地震沉积微相组合分析方法
Yanhu et al. A method of seismic meme inversion and its application
CA2293627A1 (en) Method for characterizing subsurface petrophysical properties using linear shape attributes
US11852771B1 (en) Method and system for optimally selecting carbon storage site based on multi-frequency band seismic data and equipment
CN109541685B (zh) 一种河道砂体识别方法
CN111522063B (zh) 基于分频联合反演的叠前高分辨率流体因子反演方法
CN111506861B (zh) 一种目的层有利区域裂缝强度计算方法
CN111856566A (zh) 湖相滩坝砂体中薄储层预测方法及装置
CN111897011A (zh) 储层孔隙特征确定方法、装置及设备
Eikrem et al. Bayesian estimation of reservoir properties—effects of uncertainty quantification of 4D seismic data
CN111077578B (zh) 岩层分布预测方法和装置
CN113109875B (zh) 一种全波形速度场约束下的盐下碳酸盐岩储层反演方法
CN109283577B (zh) 一种地震层位标定方法
CN114152985B (zh) 一种确定地下古河道边界及其内部薄砂体厚度的方法
CN115857047A (zh) 一种地震储层综合预测方法
CN113589365B (zh) 基于时频域信息的储层尖灭线描述方法
CN115629417A (zh) 基于地震沉积学的多尺度融合及相控颗粒滩刻画方法
CN113806674A (zh) 古河道纵向尺度的量化方法、装置、电子设备及存储介质
CN110244383B (zh) 基于近地表数据的地质岩性综合模型创建方法
CN113514884A (zh) 一种致密砂岩储层预测方法

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