CN113759424A - 基于频谱分解和机器学习的岩溶储层充填分析方法和系统 - Google Patents

基于频谱分解和机器学习的岩溶储层充填分析方法和系统 Download PDF

Info

Publication number
CN113759424A
CN113759424A CN202111066284.0A CN202111066284A CN113759424A CN 113759424 A CN113759424 A CN 113759424A CN 202111066284 A CN202111066284 A CN 202111066284A CN 113759424 A CN113759424 A CN 113759424A
Authority
CN
China
Prior art keywords
data
seismic
curve
logging
well
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.)
Granted
Application number
CN202111066284.0A
Other languages
English (en)
Other versions
CN113759424B (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.)
Institute of Geology and Geophysics of CAS
Original Assignee
Institute of Geology and Geophysics of CAS
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 Institute of Geology and Geophysics of CAS filed Critical Institute of Geology and Geophysics of CAS
Priority to CN202111066284.0A priority Critical patent/CN113759424B/zh
Publication of CN113759424A publication Critical patent/CN113759424A/zh
Application granted granted Critical
Publication of CN113759424B publication Critical patent/CN113759424B/zh
Priority to US17/697,544 priority patent/US11802985B2/en
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. analysis, for interpretation, for correction
    • 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
    • G01V1/44—Seismology; Seismic or acoustic prospecting or detecting specially adapted for well-logging using generators and receivers in the same well
    • G01V1/48—Processing data
    • G01V1/50—Analysing data
    • 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. analysis, for interpretation, for correction
    • G01V1/30—Analysis
    • G01V1/301—Analysis for determining seismic cross-sections or geostructures
    • G01V1/302—Analysis for determining seismic cross-sections or geostructures in 3D data cubes
    • 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. analysis, for interpretation, for correction
    • G01V1/30—Analysis
    • G01V1/307—Analysis for determining seismic attributes, e.g. amplitude, instantaneous phase or frequency, reflection strength or polarity
    • 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. analysis, for interpretation, for correction
    • G01V1/30—Analysis
    • G01V1/308—Time lapse or 4D effects, e.g. production related effects to the formation
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V5/00—Prospecting or detecting by the use of nuclear radiation, e.g. of natural or induced radioactivity
    • G01V5/04—Prospecting or detecting by the use of nuclear radiation, e.g. of natural or induced radioactivity specially adapted for well-logging
    • G01V5/08—Prospecting or detecting by the use of nuclear radiation, e.g. of natural or induced radioactivity specially adapted for well-logging using primary nuclear radiation sources or X-rays
    • G01V5/12—Prospecting or detecting by the use of nuclear radiation, e.g. of natural or induced radioactivity specially adapted for well-logging using primary nuclear radiation sources or X-rays using gamma or X-ray sources
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N20/00—Machine learning
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N3/00—Computing arrangements based on biological models
    • G06N3/02—Neural networks
    • G06N3/08—Learning methods
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N7/00—Computing arrangements based on specific mathematical models
    • G06N7/01—Probabilistic graphical models, e.g. probabilistic networks
    • 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. analysis, for interpretation, for correction
    • 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/40—Transforming data representation
    • G01V2210/42—Waveform, i.e. using raw or pre-filtered trace data
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00—Details of seismic processing or analysis
    • G01V2210/40—Transforming data representation
    • G01V2210/43—Spectral
    • 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/61—Analysis by combining or comparing a seismic data set with other data
    • G01V2210/614—Synthetically generated data
    • 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/61—Analysis by combining or comparing a seismic data set with other data
    • G01V2210/616—Data from specific type of measurement
    • G01V2210/6163—Electromagnetic
    • 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/61—Analysis by combining or comparing a seismic data set with other data
    • G01V2210/616—Data from specific type of measurement
    • G01V2210/6167—Nuclear
    • 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/61—Analysis by combining or comparing a seismic data set with other data
    • G01V2210/616—Data from specific type of measurement
    • G01V2210/6169—Data from specific type of measurement using 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/63—Seismic attributes, e.g. amplitude, polarity, instant phase
    • 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/63—Seismic attributes, e.g. amplitude, polarity, instant phase
    • G01V2210/632—Amplitude variation versus offset or angle of incidence [AVA, AVO, AVI]
    • 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/64—Geostructures, e.g. in 3D data cubes
    • 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/64—Geostructures, e.g. in 3D data cubes
    • G01V2210/643—Horizon tracking
    • 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/66—Subsurface modeling
    • G01V2210/665—Subsurface modeling using geostatistical modeling

Abstract

本发明属于数据识别、记录载体的处理领域，具体涉及了一种基于频谱分解和机器学习的岩溶储层充填分析方法和系统，旨在解决现有的石油勘探技术无法预测横向变化快的储层、无法识别大范围内复杂盆地碳酸盐岩洞穴型储层的发育特征的问题。本发明包括：获取标准化测井曲线数据；通过混合相位子波反褶积和扩散滤波，获得高精度的三维地震振幅数据体；确定能够反映整体地质情况的最优样本数参量进而确定时深转化关系和对储层敏感的特征参数；通过波阻抗解释结论门槛值计算缝洞型储层结构特征与空间分布规律；通过解释结论门槛值进行交会分析获得古岩溶洞穴充填物三维空间形态特征。本发明达到识别大范围内复杂盆地碳酸盐岩岩溶洞穴型储层的发育特征的效果提高了刻画的精度。

Description

基于频谱分解和机器学习的岩溶储层充填分析方法和系统
技术领域
本发明属于数据识别、记录载体的处理领域，具体涉及了一种基于频谱分解和机器学习的岩溶储层充填分析方法和系统。
背景技术
近年来，盆地深层油气资源正逐步成为油气勘探领域的重心。深层碳酸盐岩油气勘探潜力巨大，该类油藏占全球已探明石油储量的52%，产量的60%，是世界上许多盆地中最重要的油气靶区之一。虽然目前我国大部分油气产自碎屑岩储层，但碳酸盐岩储层正在发挥越来越重要的作用。中国海相碳酸盐岩储层油气资源丰富，在塔里木、鄂尔多斯、四川和渤海湾等陆上海相盆地已发现大量的石油资源，所探明石油储量为340×108t，然而，探明率仅为4.56%和13.17%。因此，海相碳酸盐岩储层具有巨大的油气开采潜力，是我国重要的油气勘探接替领域。
经研究发现，岩溶储层的充填特性会对油气产能产生很大的影响，如未充填、半充填的古岩溶洞穴有利于油气运移，能有效储集并在后期调整中良好成藏；而全充填的古岩溶洞穴则难以形成有效储集空间。据学者统计，砂泥质、洞穴垮塌角砾以及化学充填物等占据了古岩溶洞穴储层70%以上的空间。充填物岩性、物性的差异严重影响了油气储层储集空间与渗流作用。
碳酸盐岩古岩溶储层空间分布复杂，相互联系紧密。古岩溶储层研究的关键问题是如何准确认识深埋古岩溶系统的构造，如何利用地球物理资料有效地识别其位置，并且高精度地刻画古岩溶洞穴系统内部充填特征。
常规地球物理探测手段面临储层识别的瓶颈问题。由于沉积、成岩作用及后期岩溶改造，所形成的碳酸盐岩古岩溶储层洞穴、裂缝离散随机，具有极强的非均质性，我们从浅层碎屑岩和多孔碳酸盐岩储层得到的探测经验不能直接用于塔河油田的古岩溶储层勘探。稀疏的测井曲线纵向分辨率高，但是其探测范围有限，远离井筒的储层参数无法准确判断，这些特征决定了传统的测井探测手段量化储集体精细预测难度大。
由于古岩溶洞穴型储层埋藏深，地震波能量吸收衰减作用强烈，导致地震资料主频显著降低，分辨率大打折扣，深层时窗内信噪比过低。而钻井资料和岩芯样品表明，大多数洞穴高度不足5米，最大高度约为15米，原始地震资料分辨率难以达到准确识别效果。此外，与测井曲线相比，地震属性所包含的岩性物性信息含量较低，无法有效识别储层充填特征。
因此，如何充分利用地球物理信息，建立高精度测井解释方法与大范围地震探测数据的联系是改善储层预测分辨率与纵横向探测能力的关键。
井震联合反演可以将地震反射数据转化为岩石性质的定量估计，被广泛应用于储层预测中。基于地质统计学的井震联合反演是利用地质统计学原理，以地质模式为指导，根据随距离变化的变差函数将测井、地震等多类、多尺度数据整合，可实现古岩溶洞穴与致密围岩的区分，在此基础上可表征储集体的岩性、物性参数。由于地质统计思想在储层建模中发挥了地质学优势，解决了物探领域的探测精度与探测范围的理论问题，突破了地震分辨率限制，达到了薄储层的随机模拟的效果。然而，该方法样本井的优选依赖于先验模型，并参考随距离变化的变差函数。地震资料只是起到对随机模拟结果优选的作用，没有考虑地震信息对样本优选的贡献。导致反演结果在高频段部分直接受控于井间插值，横向预测随机性强，无法预测横向变化快的储层。
发明内容
为了解决目前储层识别的技术难题，即现有的石油勘探洞穴型储层填充特征分析方法只将地质资料作为对随机模拟结构优选，并未考虑地震信息对测井样本优选的贡献，导致反演结果在高频段部分直接受控于井间插值，横向预测随机性强，无法预测横向变化快的储层的问题。本发明提供了一种基于频谱分解和机器学习的岩溶储层充填分析方法，所述方法包括：
步骤S100，获取原始地球物理测井资料：通过测井设备获取各样本井的原始测井数据，包括：通过测量电极测量各样本井的自然电位SP，通过自然伽马井下装置和自然伽马地面仪器测量各样本井的自然伽马GR，通过井径臂获取各样本井的井径CAL；通过传统测井设备获取电阻率曲线数据：深侧向测井RLLD、浅侧向测井RLLS和微侧向测井RLLM，物性表征曲线数据：补偿中子CNL、补偿声波曲线AC和密度曲线DEN；
基于成像测井信息、钻井信息、录井信息和岩心信息获得个别深度段的确定岩性信息和物性信息，进而确定目的层位标志层的深度数据；
步骤S200，获取地震数据，通过地震波激发装置与接收装置获取原始地震波反射信号数据，并根据原始地震波反射信号数据的波形获取目的层位标志层的等时三维展布；
步骤S300，原始地球物理测井数据预处理：基于所述样本井的所有原始测井数据，绘制测井曲线数据，并进行异常值处理和标准化处理，获得标准化测井曲线数据；
步骤S400，地震数据预处理：基于所述地震波反射信号数据，通过混合相位子波反褶积和扩散滤波，获得高精度的三维地震振幅数据体；
步骤S500，井震标定与特征参数选取：基于所述标准化测井曲线数据中的补偿声波曲线AC和密度曲线DEN获取样本井的波阻抗曲线，进而计算反射系数曲线，获取雷克子波的优选频率使其与高精度三维地震振幅数据体的主频保持一致，将所述雷克子波与所述反射系数曲线进行褶积运算，得到合成地震记录，并将所述目的层位标志层的深度数据与所述目的层位标志层的等时三维展布相比对进行井震标定，计算所述合成地震记录与井旁地震道波形的相关度，当相关度大于或等于预设的第一阈值时，判定井震标定结果合格，获得测井曲线数据与地震记录的时深转化关系和对储层敏感的特征参数；
所述对储层敏感的特征参数，其获得方法包括：
通过对井旁不同地质体产生的测井参数的响应绘制柱状统计图，当某一标准化测井曲线数据的数值可以将不同测井解释结论预设的第二阈值以上的数据点区分开时，则选定该标准化测井曲线数据为对储层敏感的特征参数；所述对储层敏感的特征参数至少包括波阻抗IMP，还可包括井径CAL、自然伽马GR、自然电位SP，电阻率曲线数据：深侧向测井RLLD、浅侧向测井RLLS、微侧向测井RLLM，物性表征曲线数据：补偿中子CNL、补偿声波曲线AC、密度曲线DEN中的一个或多个；
步骤S600，构建等时格架模型：基于高精度三维地震振幅数据体所反映的沉积地层规律和所述时深转化关系构建等时格架模型；
步骤S700，井间储层参数模拟：基于各样本井的标准化测井曲线数据，确定能够反映整体地质情况的最优样本数参量，选取最优样本数参量条地震波形相关性最高的样本井的对储层敏感的特征参数数据构建初始模型，不断更正初始模型的参数，输出高精度特征值模拟结果数据体，所述高精度特征值模拟结果数据体为与所述对储层敏感的特征参数一一对应的数据体；
步骤S800，刻画井间溶洞系统边界：基于对储层敏感的特征参数中的样本井的波阻抗曲线IMP，绘制波阻抗曲线直方图，确定区分古岩溶洞穴与围岩波阻抗门槛值；
基于所述古岩溶洞穴与围岩波阻抗门槛值，将高精度特征值模拟结果数据体中超过所述区分古岩溶洞穴与围岩波阻抗门槛值的区域定义为碳酸盐岩基岩，将低于所述区分古岩溶洞穴与围岩波阻抗门槛值的区域定义为古岩溶洞穴储层发育位置，获得古岩溶洞穴储层结构特征和空间分布规律；
步骤S900，刻画洞穴内部充填岩性物性边界：基于所述对储层敏感的特征参数构建交会版图，根据所述个别深度段的岩性信息和物性信息在交汇版图中框选数据点作为测井解释样本点，获得充填程度与充填物类型的解释结论门槛值；
通过所述解释结论门槛值对所述高精度特征值模拟结果数据体进行交会分析获得洞穴内部充填程度与充填物组合三维空间形态特征；
步骤S1000，古岩溶洞穴结构与充填描述：基于所述古岩溶洞穴储层结构特征、空间分布规律、基于洞穴内部充填程度和充填物组合三维空间形态特征，采用岩性遮挡技术和三维雕刻技术雕刻出古岩溶洞穴空间展布和内部充填的发育特征。
在一些优选的实施方式中，所述通过测量电极测量各样本井的自然电位SP，通过自然伽马井下装置和自然伽马地面仪器测量各样本井的自然伽马GR，具体包括：
所述通过测量电极测量各样本井的自然电位SP：
将测量电极N设置于地面，测量电机M通过电缆设置于井下；
沿井轴提升测量电极M测量自然电位随井深的变化；
自然电位值的计算方法为：
Figure 681907DEST_PATH_IMAGE001
其中，
Figure 221473DEST_PATH_IMAGE002
为总自然电位，
Figure 602775DEST_PATH_IMAGE003
为扩散电位系数，
Figure 708135DEST_PATH_IMAGE004
为扩散吸附电位系数，
Figure 657636DEST_PATH_IMAGE005
为泥 浆滤液电阻率值，
Figure 789890DEST_PATH_IMAGE006
为地层水电阻率值；
所述通过自然伽马井下装置和自然伽马地面仪器测量各样本井的自然伽马GR：
自然伽马井下装置包括探测器、放大器、高压电源；
通过探测器获取自然伽马射线，并将所述自然伽马射线转化为电脉冲信号，并通过放大器进行放大；
所述地面仪器把每分钟形成的电脉冲计数转化为电位差进行记录。
在一些优选的实施方式中，所述步骤S300，包括：
步骤S310，基于所述原始测井数据绘制原始测井曲线数据；
步骤S320，基于所述原始测井曲线数据，去除离群点获得去除离群点的测井曲线数据；
步骤S330，基于所述去除离群点的测井曲线数据，叠合工区内所有样本井位的单个测井曲线柱状图数据，通过整合阈值，获得标准化测井曲线数据。
在一些优选的实施方式中，所述步骤S400，具体包括：
步骤S410，基于所述地震波反射信号数据，将频率域地震记录褶积模型表示为：
Figure 850250DEST_PATH_IMAGE007
其中，
Figure 810116DEST_PATH_IMAGE008
表示经傅氏变换后的地震记录，
Figure 55153DEST_PATH_IMAGE009
表示经傅氏变换后的子波，
Figure 303731DEST_PATH_IMAGE010
表示经傅氏变换后的反射系数的频谱，
Figure 964520DEST_PATH_IMAGE011
表示角频率；
步骤S420，将所述频率域地震记录褶积模型的等式两边取对数转化为线性系统，获得线性地震记录褶积模型：
Figure 857521DEST_PATH_IMAGE012
步骤S430，将所述线性地震记录褶积模型进行反傅氏变换，获得复赛谱序列：
Figure 148825DEST_PATH_IMAGE013
其中，
Figure 9333DEST_PATH_IMAGE014
表示地震波形记录的复赛谱序列，
Figure 473813DEST_PATH_IMAGE015
表示地震子波的复赛谱序列，
Figure 345954DEST_PATH_IMAGE016
表示地层反射系数的复赛谱序列，
Figure 182060DEST_PATH_IMAGE017
表示地震波形记录时间；
步骤S440，基于所述复赛谱序列，通过低通滤波器进行子波与反射系数分离，提取子波振幅谱；
步骤S450，通过最小二乘法获取模拟地震子波振幅谱：
Figure 201969DEST_PATH_IMAGE018
其中，最小二乘法拟合参数
Figure 407822DEST_PATH_IMAGE019
为常数，
Figure 259104DEST_PATH_IMAGE020
表示子波振幅谱，
Figure 157790DEST_PATH_IMAGE021
和
Figure 664994DEST_PATH_IMAGE022
表示f 的多项式，f表示地震子波的频率；
步骤S460，基于所述模拟地震子波振幅谱，获得子波最大相位分量和最小相位分量；
设子波
Figure 18746DEST_PATH_IMAGE023
的最大相位分量为
Figure 131059DEST_PATH_IMAGE024
、最小相位分量为
Figure 731804DEST_PATH_IMAGE025
，则子波
Figure 54201DEST_PATH_IMAGE023
为：
Figure 336278DEST_PATH_IMAGE026
振幅谱的复赛谱中表示为：
Figure 99835DEST_PATH_IMAGE027
其中，振幅谱的复赛谱
Figure 717154DEST_PATH_IMAGE028
在复赛谱的正、负轴上对称显示，
Figure 871055DEST_PATH_IMAGE029
为地震子波最 大相位分量
Figure 347036DEST_PATH_IMAGE030
所对应的最小相位函数的复赛谱，
Figure 965099DEST_PATH_IMAGE031
为地震子波最小相位分量
Figure 110910DEST_PATH_IMAGE032
所 对应的最大相位函数的复赛谱；
步骤S470，基于所述振幅谱中的复赛谱确定一组具有相同振幅谱的混合相位子波集合，不断调整俞氏子波参数，保持低频、拓展高频和适当提高主频构建期望输出子波形态，在井曲线控制下以信噪比谱作参考寻找提高分辨率与保真度之间的最佳平衡点，获得整形后波形数据；
步骤S480，基于所述整形后波形数据，构建张量扩散模型：
Figure 893052DEST_PATH_IMAGE033
其中，
Figure 579248DEST_PATH_IMAGE034
表示扩散时间，
Figure 520659DEST_PATH_IMAGE035
表示散度算子，D表示扩散滤波器的张量型扩散系数，U 表示扩散滤波结果，
Figure 962005DEST_PATH_IMAGE036
表示
Figure 621656DEST_PATH_IMAGE034
=0时的扩散滤波结果，
Figure 423128DEST_PATH_IMAGE037
表示
Figure 953467DEST_PATH_IMAGE038
时刻的整形后波形数 据，作为张量扩散模型的初始条件，
Figure 503397DEST_PATH_IMAGE039
表示扩散滤波结果的梯度；
基于所述张量扩散模型构建梯度结构张量：
Figure 243820DEST_PATH_IMAGE040
其中，U表示扩散滤波结果，
Figure 740660DEST_PATH_IMAGE041
表示梯度向量张量积；
Figure 187822DEST_PATH_IMAGE042
表示尺度为
Figure 721702DEST_PATH_IMAGE043
的高斯函数为：
Figure 824788DEST_PATH_IMAGE044
其中，r表示计算半径，
结构张量的特征向量为：
Figure 656477DEST_PATH_IMAGE045
其中，
Figure 551621DEST_PATH_IMAGE046
、
Figure 115458DEST_PATH_IMAGE047
和
Figure 768156DEST_PATH_IMAGE048
表示为梯度结构张量的3个特征向量，可视为局部正交坐标系，
Figure 977771DEST_PATH_IMAGE046
指向地震信号的梯度方向，
Figure 71629DEST_PATH_IMAGE047
和
Figure 134262DEST_PATH_IMAGE048
组成的平面平行于地震信号的局部结构特征，
Figure 70994DEST_PATH_IMAGE049
、
Figure 713328DEST_PATH_IMAGE050
和
Figure 989589DEST_PATH_IMAGE051
为分别与
Figure 301753DEST_PATH_IMAGE046
、
Figure 601147DEST_PATH_IMAGE047
和
Figure 578330DEST_PATH_IMAGE048
对应的三个特征值；
步骤S490，基于所述结构张量的特征向量分别计算线状结构置信度量、面状结构置信度量和扩散张量；
所述现状结构置信度量
Figure 36993DEST_PATH_IMAGE052
为：
Figure 379113DEST_PATH_IMAGE053
所述面状结构置信度量
Figure 228120DEST_PATH_IMAGE054
为：
Figure 320579DEST_PATH_IMAGE055
所述扩散张量D为：
Figure 509115DEST_PATH_IMAGE056
其中，
Figure 84453DEST_PATH_IMAGE057
、
Figure 483073DEST_PATH_IMAGE058
和
Figure 5321DEST_PATH_IMAGE059
表示扩散张量的三个非负特征值，他们分别表示扩散滤波器沿
Figure 110681DEST_PATH_IMAGE046
、
Figure 935548DEST_PATH_IMAGE047
和
Figure 431252DEST_PATH_IMAGE048
这三个特征方向的滤波强度；
步骤S4100，重复步骤S380-S390的步骤，直至达到预设的迭代次数，获得扩散滤波结果，即为所述高精度的三维地震振幅数据体。
在一些优选的实施方式中，所述模拟地震子波振幅谱的计算方法为：
定位地震波反射信号数据中振幅谱的极大值和极大值所对应的频率；
将所述地震信号振幅谱的极大值和所述模拟地震子波振幅谱通过最小二乘法拟 合的方式得到参数
Figure 553928DEST_PATH_IMAGE060
和
Figure 576111DEST_PATH_IMAGE061
多项式的系数，进而获得拟合的极大值相应频率振幅值；
将所述地震信号振幅谱的极大值除以所拟合出的相应频率的振幅值，进而用商拟 合多项式
Figure 696514DEST_PATH_IMAGE022
的系数。
在一些优选的实施方式中，所述测井曲线数据与地震记录的时深转化关系的时深转化关系为：
Figure 7409DEST_PATH_IMAGE062
其中，
Figure 248291DEST_PATH_IMAGE063
表示声波测井起始深度对应的地震信号时间，
Figure 347DEST_PATH_IMAGE064
表示声波时差，
Figure 353968DEST_PATH_IMAGE065
表示表 示每个采样点之间的时间运算序号，
Figure 214476DEST_PATH_IMAGE066
表示测井曲线数据采样间隔，
Figure 616639DEST_PATH_IMAGE067
表示地震波双程旅 行时。
在一些优选的实施方式中，所述步骤S700，具体包括步骤S710-步骤S790：
步骤S710，任选一样本井作为参考目标井，设定初始的样本数参量为1；
步骤S720，根据波形相似性原则选取数量为样本数参量的样本井标准化测井曲线的特征参数数据与参考目标井标准化测井曲线的特征参数数据进行相关性分析获得样本数参量-参考目标井曲线特征参数的相关性值；
步骤S730，逐1增加样本数参量，重复步骤S620的方法获得各样本数参量对应的样本数参量-参考目标井曲线特征参数的相关性值，将所有样本数参量-参考目标井曲线特征参数的相关性值连接，获得参考井曲线特征参数相关性随样本数参量变化的相关性曲线；
步骤S740，选取另一样本井作为参考目标井，重复步骤S710-步骤S730的方法，获得多个参考井的曲线特征参数相关性随样本数参量变化的相关性曲线，将所有参考井的特征参数的相关性随样本数参量变化的相关性曲线拟合为整体相关性曲线，选取所述整体相关性曲线中相关性随样本参数增加而升高最终保持平稳处的拐点，确定最优样本数参量；
步骤S750，基于所述高精度的三维地震振幅数据体和等时格架模型，计算待测点位与样本井位的波形相关性，将所述波形相关性由大到小排序，选取最优样本数参量条地震波形相关性最高的样本井的对储层敏感的特征参数数据；基于所述相关性最高的样本井的地震波形特征数据所对应的样本井，通过井间特征参数插值方式构建初始模型；
步骤S760，基于所述初始模型，选取最优样本数参量条地震波形关联度最高的样本井待模拟曲线参数对应的对储层敏感的特征参数作为先验信息；所述待模拟曲线参数为与所述对储层敏感的特征参数一一对应的曲线参数，至少包括波阻抗IMP，还可包括井径CAL、自然伽马GR、自然电位SP，电阻率曲线数据：深侧向测井RLLD、浅侧向测井RLLS、微侧向测井RLLM，物性表征曲线数据：补偿中子CNL、补偿声波AC、密度曲线DEN中的一个或多个；
步骤S770，将所述初始模型与先验信息进行匹配滤波，进而获得最大似然函数；
步骤S780，基于所述最大似然函数和先验信息，在贝叶斯框架下求得后验概率统计分布密度，对所述后验概率统计分布密度进行采样获得目标函数；
步骤S790，以所述目标函数作为所述初始模型的输入，通过马尔科夫链蒙特卡罗方法（Markov chain Monte Carlo method，MCMC）和Metropolis-Hastings抽样准则对后验概率分布抽样，不断优化初始模型的参数，选取目标函数取最大值时的解作为随机实现，取多次随机实现的均值作为期望值输出，将所述期望值输出作为高精度特征值模拟结果数据体；所述高精度特征值模拟结果数据体中的参数与所述对储层敏感的特征参数一一对应。
在一些优选的实施方式中，所述步骤S780，具体包括：
步骤S781，利用白噪声满足高斯分布的规律，将高精度特征值模拟结果数据体的参数表示为：
Figure 285517DEST_PATH_IMAGE068
其中，Y表示测井曲线高精度特征值模拟结果数据体的参数，X表示待求解的地下地层实际特征参数值，N表示随机噪声；
步骤S782，因为
Figure 888668DEST_PATH_IMAGE069
也满足高斯分布，可以确定初始目标函数为：
Figure 846260DEST_PATH_IMAGE070
其中，
Figure 848851DEST_PATH_IMAGE071
表示与后验信息有关的函数，
Figure 965712DEST_PATH_IMAGE072
表示基于最优样本数选取样本井对样本 井的特征曲线进行匹配滤波后，求得后验概率统计分布密度，进而计算得到的特征参数期 望值，
Figure 598818DEST_PATH_IMAGE073
表示白噪声的协方差；
步骤S783，基于所述初始目标函数，通过最大后验估计，在目标函数中引入先验信息，获得稳定的目标函数为：
Figure 840444DEST_PATH_IMAGE074
其中，
Figure 223889DEST_PATH_IMAGE075
表示待模拟的特征参数，
Figure 70623DEST_PATH_IMAGE076
表示与地质、测井资料等先验信息有关的函数，
Figure 671368DEST_PATH_IMAGE077
表示用于协调
Figure 993765DEST_PATH_IMAGE078
和
Figure 275842DEST_PATH_IMAGE076
之间的相互影响的平滑参数。
在一些优选的实施方式中，所述步骤S690，具体步骤为：
步骤S791，设M为目标空间，n为总样本数，m为马尔科夫链趋于平稳时的样本数；
步骤S792，预设一条马尔科夫链，使马尔科夫链收敛至平稳分布；
步骤S793，由M中的某一点
Figure 39399DEST_PATH_IMAGE079
出发，通过马尔科夫链进行抽样模拟，产生点序列：
Figure 889674DEST_PATH_IMAGE080
；
步骤S794，函数
Figure 309154DEST_PATH_IMAGE081
的期望估计为：
Figure 457239DEST_PATH_IMAGE082
其中，n表示产生的总样本数，m表示马尔科夫链达到平稳时的样本数，k表示累加参量；
步骤S795，选取一转移函数
Figure 403198DEST_PATH_IMAGE083
和初始值
Figure 549008DEST_PATH_IMAGE084
，若第i次迭代开始时的参数 值为
Figure 252522DEST_PATH_IMAGE085
，则第i次迭代过程为：
从
Figure 512953DEST_PATH_IMAGE083
中抽取一个备选值
Figure 188785DEST_PATH_IMAGE086
，计算备选值
Figure 302234DEST_PATH_IMAGE086
的接受概率
Figure 820940DEST_PATH_IMAGE087
：
Figure 248510DEST_PATH_IMAGE088
步骤S796，以
Figure 919794DEST_PATH_IMAGE089
，置
Figure 469724DEST_PATH_IMAGE090
，以概率
Figure 85513DEST_PATH_IMAGE091
，置
Figure 441408DEST_PATH_IMAGE092
；
步骤S797，不断扰动所述初始模型的参数，重复步骤S792-步骤S796的方法达到预 设的迭代次数n，获得后验样本
Figure 888570DEST_PATH_IMAGE093
，进而计算后验分布的各阶矩阵获得期望 输出值，将期望值输出作为高精度特征值模拟结果数据体。
在一些优选的实施方式中，所述步骤S900，优选的包括：
步骤S910，以波阻抗值与自然GR值作为对储层敏感的特征参数构建交会图模板；
步骤S920，基于所述交会图模板，根据所述个别深度段岩性信息和物性信息，将样本点划分为全填充、半填充、基岩；
步骤S930，框选储集性能较好的半填充样本点，计算得到二维交会图解释结论门槛值；
步骤S940，将所述二维交会图解释结论门槛值输入三维自然GR特征参数模拟数据体与波阻抗反演数据体，刻画半填充储层的空间结构形态，获得古岩溶洞穴半充填部分三维空间形态特征。
在充填物识别的实施方式中，所述步骤S900，还包括获取其他古岩溶洞穴特征的步骤：
步骤S950，基于所述交会图模板，根据所述个别深度段岩性信息和物性信息，将样本点划分为基岩、沉积充填、化学充填、垮塌充填和混合充填，框选各类填充物样本点，分别计算对应解释结论门槛值，分别刻画对应的空间结构形态，获得古岩溶洞穴充填物三维空间形态特征。
本发明的另一方面，提出了一种基于频谱分解和机器学习的岩溶储层充填分析系统，所述系统包括：原始地球物理测井资料获取模块、地震数据获取模块、原始地球物理测井数据预处理模块、地震数据预处理模块、井震标定与特征参数选取模块、等时格架模型构建模块、井间特征参数模拟模块、井间溶洞系统边界刻画模块、溶洞内部充填岩性物性边界刻画模块、古岩溶洞穴结构与充填特征三维描述模块；
所述原始地球物理测井资料获取模块，配置为通过测井设备获取各样本井的原始测井数据，包括：通过测量电极测量各样本井的自然电位SP，通过自然伽马井下装置和自然伽马地面仪器测量各样本井的自然伽马GR，通过井径臂获取各样本井的井径CAL；通过传统测井设备获取电阻率曲线数据：深侧向测井RLLD、浅侧向测井RLLS和微侧向测井RLLM，物性表征曲线数据：补偿中子CNL、补偿声波曲线AC和密度曲线DEN；
基于成像测井信息、钻井信息、录井信息和岩心信息获得个别深度段的确定岩性信息和物性信息，进而确定目的层位标志层的深度数据；
所述地震数据获取模块，配置为通过地震波激发装置与接收装置获取原始地震波反射信号数据，并根据原始地震波反射信号数据的波形获取目的层位标志层的等时三维展布；
所述原始地球物理测井数据预处理模块，配置为基于所述样本井的原始测井数据，绘制测井曲线数据，并进行异常值处理和标准化处理，获得标准化测井曲线数据；
所述地震数据预处理模块，配置为基于所述地震波反射信号数据，通过混合相位子波反褶积和扩散滤波，获得高精度的三维地震振幅数据体；
所述井震标定与特征参数选取模块，配置为基于所述标准化测井曲线数据中的补偿声波曲线AC和密度曲线DEN获取样本井的波阻抗曲线，进而计算反射系数曲线，获取雷克子波的优选频率使其与高精度三维地震振幅数据体的主频保持一致，将所述雷克子波与所述反射系数曲线进行褶积运算，得到合成地震记录，并将所述目的层位标志层的深度数据与所述目的层位标志层的等时三维展布相比对进行井震标定，计算所述合成地震记录与井旁地震道波形的相关度，当相关度大于或等于预设的第一阈值时，判定井震标定结果合格，获得测井曲线数据与地震记录的时深转化关系和对储层敏感的特征参数；
所述对储层敏感的特征参数，其获得方法包括：
通过对井旁不同地质体产生的测井参数的响应绘制柱状统计图，当某一标准化测井曲线数据的数值可以将不同测井解释结论预设的第二阈值以上的数据点区分开时，则选定该标准化测井曲线数据为对储层敏感的特征参数；所述对储层敏感的特征参数至少包括波阻抗IMP，还可包括井径CAL、自然伽马GR、自然电位SP，电阻率曲线数据：深侧向测井RLLD、浅侧向测井RLLS、微侧向测井RLLM，物性表征曲线数据：补偿中子CNL、补偿声波曲线AC、密度曲线DEN中的一个或多个；
所述等时格架模型构建模块，配置为基于高精度三维地震振幅数据体所反映的沉积地层规律和所述时深转化关系构建等时格架模型；
所述井间特征参数模拟模块，配置为基于各样本井的标准化测井曲线数据，确定能够反映整体地质情况的最优样本数参量，选取最优样本数参量条地震波形相关性最高的样本井的对储层敏感的特征参数数据构建初始模型，不断更正初始模型的参数，输出高精度特征值模拟结果数据体，所述高精度特征值模拟结果数据体为与所述对储层敏感的特征参数一一对应的数据体；
所述井间溶洞系统边界刻画模块，配置为基于样本井对储层敏感的特征参数中的波阻抗曲线IMP，绘制波阻抗曲线直方图，确定区分古岩溶洞穴与围岩波阻抗门槛值；
基于所述古岩溶洞穴与围岩波阻抗门槛值，将高精度特征值模拟结果数据体中超过所述阻抗门槛值的区域定义为碳酸盐岩基岩，将低于所述阻抗门槛值的区域定义为古岩溶洞穴储层发育位置，获得古岩溶洞穴储层结构特征和空间分布规律；
所述溶洞内部充填岩性物性边界刻画模块，配置为刻画洞穴内部充填岩性物性边界：基于所述对储层敏感的特征参数构建交会版图，根据所述个别深度段的岩性信息和物性信息在交汇版图中框选数据点作为测井解释样本点，获得充填程度与充填物类型的解释结论门槛值；
通过所述解释结论门槛值对所述高精度特征值模拟结果数据体进行交会分析获得洞穴内部充填程度与充填物组合三维空间形态特征；
所述古岩溶洞穴结构与充填特征三维描述模块，配置为古岩溶洞穴结构与充填描述：基于所述古岩溶洞穴储层结构特征、空间分布规律、基于洞穴内部充填程度和充填物组合三维空间形态特征，采用岩性遮挡技术和三维雕刻技术雕刻出古岩溶洞穴空间展布和内部充填的发育特征。
本发明的第三方面，提出了一种电子设备，包括：至少一个处理器；以及至少一个所述处理器通信连接的存储器；其中，所述存储器存储有可被所述处理器执行的指令，所述指令用于被所述处理器执行以实现上述的基于频谱分解和机器学习的岩溶储层充填分析方法。
本发明的第四方面，提出了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机指令，所述计算机指令用于被所述计算机执行以实现上述的基于频谱分解和机器学习的岩溶储层充填分析方法。
本发明的有益效果：
本发明通过利用高精度的三维地震波形数据之间的关联性在初始等时地质格架模型基础上井震联合模拟阻抗以及其他储层岩性物性参数，根据测井数据得到古岩溶洞穴型储层波阻抗门槛值，充填物岩性解释门槛值，进而刻画储层成因结构与充填组合类型提高了刻画的精度，能够预测横向变化快的储层。为评价碳酸盐岩古岩溶洞穴储集性能及油气运移提供可靠的理论依据。
本发明在贝叶斯理论指导下的马尔科夫链抽样准则与蒙特卡罗估计法建立地球物理探测方法之间的映射关系，以达到识别大范围内复杂盆地碳酸盐岩古岩溶洞穴型储层的发育特征的效果。
附图说明
通过阅读参照以下附图所作的对非限制性实施例所作的详细描述，本申请的其它特征、目的和优点将会变得更明显：
图1是本发明基于频谱分解和机器学习的岩溶储层充填分析方法实施例的流程示意图；
图2是本发明去除离群点的测井曲线示意图；
图3是本发明进行标准化流程中叠合所有曲线的示意图；
图4是本发明实施例子波振幅谱图；
图5是本发明实施例中单井井震标定图；
图6是本发明实施例中波阻抗曲线柱状统计图；
图7是本发明实施例中以T74标志层以下的等时格架模型图；
图8是本发明实施例中古岩溶洞穴型储层半充填样本框选；
图9是本发明实施例中古岩溶洞穴型储层半充填储层刻画；
图10是本发明实施例中将将所有参考井的特征参数的相关性随样本数参量变化的相关性曲线拟合为整体相关性曲线示意图。
具体实施方式
下面结合附图和实施例对本申请作进一步的详细说明。可以理解的是，此处所描述的具体实施例仅用于解释相关发明，而非对该发明的限定。另外还需要说明的是，为了便于描述，附图中仅示出了与有关发明相关的部分。
需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。下面将参考附图并结合实施例来详细说明本申请。
为了更清晰地对本发明基于频谱分解和机器学习的岩溶储层充填分析方法进行说明，下面结合图1对本发明实施例中各步骤展开详述。
本发明第一实施例的基于频谱分解和机器学习的岩溶储层充填分析方法，包括步骤S100-步骤S1000，各步骤详细描述如下：
经实际地震数据统计分析表明，古岩溶洞穴所产生的地震波形特征通常由多组波谷和波峰组成，经研究发现，这种特征性反射是由地震波干涉造成的。并且反射波的振幅与洞穴充填物组合有关。Xu等（2016）通过物理模型模拟同样可以发现不同形状、厚度、组合关系的古岩溶洞穴型目标体会产生多样的地震波形反射特征，即反射形态特征对于洞穴直径、洞穴形态与分布规律有关。学者们进一步发现在等时地层格架的同一相带内，波形特征的相似性可以代表具有一定相关性的岩性组合。因此，基于波形指示的思路进行等时界面分区反演与特征参数模拟，利用了地震波形的横向变化信息，更好地体现了沉积环境的约束，更符合沉积地质规律。
步骤S100，获取原始地球物理测井资料：通过测井设备获取各样本井的原始测井数据，包括：通过测量电极测量各样本井的自然电位SP，通过自然伽马井下装置和自然伽马地面仪器测量各样本井的自然伽马GR，通过井径臂获取各样本井的井径CAL；通过传统测井设备获取电阻率曲线数据：深侧向测井RLLD、浅侧向测井RLLS和微侧向测井RLLM，物性表征曲线数据：补偿中子CNL、补偿声波曲线AC和密度曲线DEN；
基于成像测井信息、钻井信息、录井信息和岩心信息获得个别深度段的确定岩性信息和物性信息，进而确定目的层位标志层的深度数据；在本实施例中，在工区井眼处探测深度范围为5500~5750m的九条常规测井曲线数据，采样间隔设置为0.01m；运用三维地震勘探方法，借助地震波激发源与地震信号检波器，得到工区面积约27km2的三维地震资料数据，信号记录双程旅行时为4s，采样点时间间隔为1ms，探测深度超过6000m。
在本实施例中，所述通过测量电极测量各样本井的自然电位SP，通过自然伽马井下装置和自然伽马地面仪器测量各样本井的自然伽马GR，具体包括：
所述通过测量电极测量各样本井的自然电位SP：
将测量电极N设置于地面，测量电机M通过电缆设置于井下；
沿井轴提升测量电极M测量自然电位随井深的变化；
自然电位值的计算方法为：
Figure 547085DEST_PATH_IMAGE001
其中，
Figure 712487DEST_PATH_IMAGE002
为总自然电位，
Figure 121340DEST_PATH_IMAGE003
为扩散电位系数，
Figure 360692DEST_PATH_IMAGE004
为扩散吸附电位系数，
Figure 252424DEST_PATH_IMAGE005
为泥 浆滤液电阻率值，
Figure 233019DEST_PATH_IMAGE006
为地层水电阻率值；
所述通过自然伽马井下装置和自然伽马地面仪器测量各样本井的自然伽马GR：
自然伽马井下装置包括探测器、放大器、高压电源；
通过探测器获取自然伽马射线，并将所述自然伽马射线转化为电脉冲信号，并通过放大器进行放大；
所述地面仪器把每分钟形成的电脉冲计数转化为电位差进行记录。
双侧向电极系由主电极、监督电极、环状屏蔽电极与柱状电极组成。主电极
Figure 540503DEST_PATH_IMAGE094
居 中，上下对称分布监督电极
Figure 40886DEST_PATH_IMAGE095
、
Figure 103520DEST_PATH_IMAGE096
和
Figure 915618DEST_PATH_IMAGE097
、
Figure 354690DEST_PATH_IMAGE098
，以及环状屏蔽电极
Figure 958846DEST_PATH_IMAGE099
、
Figure 130065DEST_PATH_IMAGE100
，在
Figure 491776DEST_PATH_IMAGE099
、
Figure 783473DEST_PATH_IMAGE100
的外侧对 称位置上加了两个柱状电极。深侧向电极系中两个柱状电极是屏蔽电极
Figure 117503DEST_PATH_IMAGE101
、
Figure 521939DEST_PATH_IMAGE102
；浅侧向电 极系中两个柱状电极是回路电极
Figure 698843DEST_PATH_IMAGE103
和
Figure 151821DEST_PATH_IMAGE104
。在电极系较远处装有对比电极N和深侧向电极系 的回路电极B；微球形聚焦电极系的主电极
Figure 402673DEST_PATH_IMAGE094
是矩形片状电极，依次向外的矩形框电极是测 量电极
Figure 56640DEST_PATH_IMAGE105
、辅助电极
Figure 330626DEST_PATH_IMAGE099
、监督电极
Figure 977508DEST_PATH_IMAGE106
和
Figure 82868DEST_PATH_IMAGE097
，回路电极B装在仪器外壳上或极板支撑架上。测 量任一监督电极与对比电极之间的电位差变化即反映介质电阻率的变化，其视电阻率表达 式为：
Figure 766790DEST_PATH_IMAGE107
式中，
Figure 324810DEST_PATH_IMAGE108
为电极系系数，深侧向测井、浅侧向测井与微球型聚焦测井各不相同；
Figure 24651DEST_PATH_IMAGE109
为监督电极M的电位；
Figure 922199DEST_PATH_IMAGE110
为主电流。
井下仪器的发射换能器晶体振动，引起周围介质的质点发生振动，产生向井内泥 浆及岩层中传播的声波。在井中可以用接收换能器R、R2先后接收到滑行波，记录时差
Figure 839340DEST_PATH_IMAGE111
，进 而测量地层的声波速度。到达
Figure 478132DEST_PATH_IMAGE112
和
Figure 76603DEST_PATH_IMAGE113
的时刻分别为
Figure 969604DEST_PATH_IMAGE114
和
Figure 323225DEST_PATH_IMAGE115
，那么到达两个接收换能器的时 间差
Figure 793521DEST_PATH_IMAGE111
为:
Figure 320317DEST_PATH_IMAGE116
式中，
Figure 254775DEST_PATH_IMAGE046
为泥浆声速；
Figure 716980DEST_PATH_IMAGE047
为地层的声速
密度测井仪器包括有一个伽马源，两个接收伽马射线的探测器，即长源距探测器和短源距探测器。它们安装在滑板上，测井时被推靠到井壁上。在下井仪器的上方装有辅助电子线路。通常用137Cs作伽马源，它发射的伽马射线具有中等能量(0.661MeV)，用它照射物质只能产生康普顿散射和光电效应。地层的密度不同，则对伽马光子的散射和吸收的能力不同，探测器接收到的伽马光子的计数率也就不同。已知通过距离为L的伽马光子的计数率N为：
Figure 471309DEST_PATH_IMAGE117
若只存在康普顿散射，则
Figure 337151DEST_PATH_IMAGE118
即为康普顿散射吸收系数，变换得：
Figure 63799DEST_PATH_IMAGE119
式中，C为与地层有关的常数，L为接收源距离伽马源的距离。
源距选定后，对仪器进行刻度标定，找到
Figure 493643DEST_PATH_IMAGE120
和N的这种关系，则记录散射伽马光子 计数率N就可以测得地层的密度
Figure 328744DEST_PATH_IMAGE120
补偿中子测井由一个中子发射源与两个接收探测器组成。通过对两个探测器计 数，得到两个计数率
Figure 72709DEST_PATH_IMAGE121
和
Figure 716180DEST_PATH_IMAGE122
，取其比值可以反映地层的氢含量：
Figure 661133DEST_PATH_IMAGE123
式中，
Figure 327738DEST_PATH_IMAGE124
与
Figure 265607DEST_PATH_IMAGE125
为两个探测器到中子源的距离；
Figure 29164DEST_PATH_IMAGE126
为减速长度，反映了地层含氢量。
井径测量时，将仪器下人井预计的深度上，然后通过传统的方式打开井径臂，于是互成90°的四个井径腿便在弹簀力的作用下向外伸张，其末端紧贴井壁。随着仪器的向上提升，井径臂就会由于井径的变化而发生张缩，并带动连杆作上下运动。连杆与一个电位器的滑动端相连，于是井径的变化便可转换成电阻的变化。当给该可动电阻上通以一-定强度的电流时，可动电阻的某一固定端与滑动端之间的电位差将随着其间电阻值的变化而变化。于是测量这一电位差，便可反映井径大小：
Figure 738494DEST_PATH_IMAGE127
式中，
Figure 797454DEST_PATH_IMAGE128
为井径测量仪器电位差；
Figure 945539DEST_PATH_IMAGE129
为初始值；c为仪器常数。
步骤S200，获取地震数据，通过地震波激发装置与接收装置获取原始地震波反射信号数据，并根据原始地震波反射信号数据的波形获取目的层位标志层的等时三维展布；
步骤S300，原始地球物理测井数据预处理：基于所述样本井的原始测井数据，绘制测井曲线数据，并进行异常值处理和标准化处理，获得标准化测井曲线数据；
在本实施例中，所述步骤S300，包括：
步骤S310，基于所述原始测井数据绘制原始测井曲线数据；
步骤S320，基于所述原始测井曲线数据，去除离群点获得去除离群点的测井曲线数据；离群点，即在数据集中存在不合理的值，统计所有井数据的单个曲线参数直方图，合理调整保留区间阈值进行离群点去除，去除离群点的测井曲线直方图如图2所示；
步骤S330，基于所述去除离群点的测井曲线数据，叠合工区内所有样本井位的单个测井曲线柱状图数据，通过整合阈值，获得标准化测井曲线数据；以补偿声波曲线AC标准化处理为例，获得标准化测井曲线数据。由于仪器差异或受到其他因素影响，导致不同的井之间常规测井数据会出现整体偏大或整体偏小，需要将曲线进行标准化，叠合所有曲线的如图3所示。
步骤S400，地震数据预处理：地震数据预处理：基于所述地震波反射信号数据，通过混合相位子波反褶积和扩散滤波，获得高精度的三维地震振幅数据体；
目前25m×25m道间距的三维地震测网技术被广泛应用于石油勘探领域，通常在垂向上以2ms的采样间隔接收地震波反射信号，采样总时长在6s内，以探测不同深度层段的地质特征。本发明常用于2m以上的目标体探测，对地震资料主频要求较高，应在50~60Hz范围内。当古岩溶洞穴发育层段内地震资料振幅数据体主频低于50HZ时，则需要采用混合相位子波反褶积（mixed-phase wavelet estimation and maximum posteriorideconvolution）与扩散滤波（enhancing diffusion filtering）进行三维地震数据的拓频降噪处理，得到高分辨率、高信噪比的三维地震数据体。
混合相位子波反褶积是一种在确保处理后的地震数据具有较高保真度高的前提下，通过拓宽有效频带从而提高地震信号分辨率的数据处理方法，相当于本实施例中的S410-步骤S4100。
在本实施例中，步骤S400，具体包括：
步骤S410，基于所述地震波反射信号数据，将频率域地震记录褶积模型表示为：
Figure 501285DEST_PATH_IMAGE130
其中，
Figure 443833DEST_PATH_IMAGE008
表示经傅氏变换后的地震记录，
Figure 475243DEST_PATH_IMAGE009
表示经傅氏变换后的子波，
Figure 364702DEST_PATH_IMAGE010
表示经傅氏变换后的反射系数的频谱，
Figure 915900DEST_PATH_IMAGE011
表示角频率；
步骤S420，将所述频率域地震记录褶积模型的等式两边取对数转化为线性系统，获得线性地震记录褶积模型：
Figure 294929DEST_PATH_IMAGE012
步骤S430，将所述线性地震记录褶积模型进行反傅氏变换，获得复赛谱序列：
Figure 423422DEST_PATH_IMAGE013
其中，
Figure 241205DEST_PATH_IMAGE014
表示地震波形记录的复赛谱序列，
Figure 833861DEST_PATH_IMAGE015
表示地震子波的复赛谱序列，
Figure 321474DEST_PATH_IMAGE016
表示地层反射系数的复赛谱序列，
Figure 314094DEST_PATH_IMAGE017
表示地震波形记录时间；
步骤S440，基于所述复赛谱序列，通过低通滤波器进行子波与反射系数分离，提取子波振幅谱；子波与反射系数序列平滑程度的差异在复赛谱中易于区分：子波能量出现于原点附近，而反射系数序列则远离原点。利用低通滤波器就可以将复赛谱中的子波与反射系数分离达到子波振幅谱提取的目的。
步骤S450，通过最小二乘法获取模拟地震子波振幅谱：
Figure 607672DEST_PATH_IMAGE018
其中，最小二乘法拟合参数
Figure 992517DEST_PATH_IMAGE019
为常数，
Figure 41244DEST_PATH_IMAGE020
表示子波振幅谱，
Figure 206647DEST_PATH_IMAGE021
和
Figure 976019DEST_PATH_IMAGE022
表示f 的多项式，f表示地震子波的频率；
在本实施例中，所述模拟地震子波振幅谱的计算方法为：
定位地震波反射信号数据中振幅谱的极大值和极大值所对应的频率；
将所述地震信号振幅谱的极大值和所述模拟地震子波振幅谱通过最小二乘法拟 合的方式得到参数
Figure 356316DEST_PATH_IMAGE060
和
Figure 248049DEST_PATH_IMAGE061
多项式的系数，进而获得拟合的极大值相应频率振幅值；
将所述地震信号振幅谱的极大值除以所拟合出的相应频率的振幅值，进而用商拟 合多项式
Figure 572851DEST_PATH_IMAGE022
的系数；
步骤S460，基于所述模拟地震子波振幅谱，获得子波最大相位分量和最小相位分量；
设子波
Figure 536128DEST_PATH_IMAGE023
的最大相位分量为
Figure 957882DEST_PATH_IMAGE024
、最小相位分量为
Figure 223778DEST_PATH_IMAGE025
，则子波
Figure 409778DEST_PATH_IMAGE023
为：
Figure 583270DEST_PATH_IMAGE026
振幅谱的复赛谱中表示为：
Figure 797214DEST_PATH_IMAGE027
其中，振幅谱的复赛谱
Figure 358645DEST_PATH_IMAGE028
在复赛谱的正、负轴上对称显示，
Figure 720356DEST_PATH_IMAGE029
为地震子波最 大相位分量
Figure 900802DEST_PATH_IMAGE030
所对应的最小相位函数的复赛谱，
Figure 110197DEST_PATH_IMAGE031
为地震子波最小相位分量
Figure 514634DEST_PATH_IMAGE032
所 对应的最大相位函数的复赛谱；
步骤S470，基于所述振幅谱中的复赛谱确定一组具有相同振幅谱的混合相位子波集合，不断调整俞氏子波参数，保持低频、拓展高频和适当提高主频构建期望输出子波形态，在井曲线控制下以信噪比谱作参考寻找提高分辨率与保真度之间的最佳平衡点，获得整形后波形数据；子波振幅谱如图4所示；
经过混合相位子波反褶积后，拓展了地震数据的有效频带，高频部分得到了合理加强。在地震波形上表现为同相轴数目增多，更容易反映地震波反射信息的细节变化，并且在振幅、相位、频率方面改善了同一反射波组波形的一致性。在古岩溶洞穴地震响应上，“串珠状”反射特征尤为明显，且串珠内部形态的细节能够清晰地显示，代表了不同结构特征与充填物组合的复杂古岩溶洞穴型储层地震反射，有助于后期精细地质解释。
Fhemers和Hocker于2003年首次将扩散滤波技术应用在地震资料的处理解释中。该技术不仅可以有效压制噪声，而且可以尽可能保留地震数据中的细节：如地质体边缘、断层、不整合面、尖灭等，为后续的地震解释和储层预测工作提供了可靠的地震资料，大大提高油气勘探开发的成功率。
为达到衰减地震噪声，并增强地质结构特征的扩散滤波效果，寻求扩散张量是该方法最关键的步骤：
步骤S480，基于所述整形后波形数据，构建张量扩散模型：
Figure 35745DEST_PATH_IMAGE131
其中，
Figure 144516DEST_PATH_IMAGE034
表示扩散时间，
Figure 395368DEST_PATH_IMAGE035
表示散度算子，D表示扩散滤波器的张量型扩散系数，U 表示扩散滤波结果，
Figure 908389DEST_PATH_IMAGE036
表示
Figure 553347DEST_PATH_IMAGE034
=0时的扩散滤波结果，
Figure 872333DEST_PATH_IMAGE037
表示
Figure 180955DEST_PATH_IMAGE038
时刻的整形后波形数 据，作为张量扩散模型的初始条件，
Figure 989511DEST_PATH_IMAGE039
表示扩散滤波结果的梯度；
基于所述张量扩散模型构建梯度结构张量：
Figure 813110DEST_PATH_IMAGE132
其中，U表示扩散滤波结果，
Figure 607891DEST_PATH_IMAGE041
表示梯度向量张量积；
Figure 911965DEST_PATH_IMAGE042
表示尺度为
Figure 829105DEST_PATH_IMAGE043
的高斯函数为：
Figure 812105DEST_PATH_IMAGE044
其中，r表示计算半径；
结构张量的特征向量为：
Figure 800789DEST_PATH_IMAGE133
其中，
Figure 880741DEST_PATH_IMAGE046
、
Figure 906465DEST_PATH_IMAGE047
和
Figure 16242DEST_PATH_IMAGE048
表示为梯度结构张量的3个特征向量，可视为局部正交坐标系，
Figure 480721DEST_PATH_IMAGE046
指向地震信号的梯度方向，
Figure 352862DEST_PATH_IMAGE047
和
Figure 939701DEST_PATH_IMAGE048
组成的平面平行于地震信号的局部结构特征，
Figure 959610DEST_PATH_IMAGE049
、
Figure 431042DEST_PATH_IMAGE050
和
Figure 33056DEST_PATH_IMAGE051
为分别与
Figure 462901DEST_PATH_IMAGE046
、
Figure 907788DEST_PATH_IMAGE047
和
Figure 41966DEST_PATH_IMAGE048
对应的三个特征值；
步骤S490，基于所述结构张量的特征向量分别计算线状结构置信度量、面状结构置信度量和扩散张量；
所述现状结构置信度量
Figure 685437DEST_PATH_IMAGE052
为：
Figure 489445DEST_PATH_IMAGE134
所述面状结构置信度量
Figure 798460DEST_PATH_IMAGE054
为：
Figure 408433DEST_PATH_IMAGE135
所述扩散张量D为：
Figure 109673DEST_PATH_IMAGE056
其中，
Figure 209216DEST_PATH_IMAGE057
、
Figure 628696DEST_PATH_IMAGE058
和
Figure 776781DEST_PATH_IMAGE059
表示扩散张量的三个非负特征值，他们分别表示扩散滤波器沿
Figure 207893DEST_PATH_IMAGE046
、
Figure 416021DEST_PATH_IMAGE047
和
Figure 322797DEST_PATH_IMAGE048
这三个特征方向的滤波强度；
步骤S4100，重复步骤S480-S490的步骤，直至达到预设的迭代次数，获得扩散滤波结果，即为所述高精度的三维地震振幅数据体。扩散滤波算法保留了“串珠状”反射古岩溶洞穴型储层地质特征，增强了地震数据对目标地质体的成像能力。同时起到了压制噪声、提高同相轴横向连续性与地震信号的信噪比的效果。
步骤S500，井震标定与特征参数选取：基于所述标准化测井曲线数据中的补偿声波曲线AC和密度曲线DEN获取样本井的波阻抗曲线，进而计算反射系数曲线，获取雷克子波的优选频率使其与高精度三维地震振幅数据体的主频保持一致，将所述雷克子波与所述反射系数曲线进行褶积运算，得到合成地震记录，并将所述目的层位标志层的深度数据与所述目的层位标志层的等时三维展布相比对进行井震标定，计算所述合成地震记录与井旁地震道波形的相关度，当相关度大于或等于预设的第一阈值时，判定井震标定结果合格，获得测井曲线数据与地震记录的时深转化关系和对储层敏感的特征参数；所述第一阈值，优选的选取85%；
所述对储层敏感的特征参数，其获得方法包括：
通过对井旁不同地质体产生的测井参数的响应绘制柱状统计图，当某一标准化测井曲线数据的数值可以将不同测井解释结论预设的第二阈值以上的数据点区分开时，则选定该标准化测井曲线数据为对储层敏感的特征参数；所述对储层敏感的特征参数至少包括波阻抗IMP，还可包括井径CAL、自然伽马GR、自然电位SP，电阻率曲线数据：深侧向测井RLLD、浅侧向测井RLLS、微侧向测井RLLM，物性表征曲线数据：补偿中子CNL、补偿声波曲线AC、密度曲线DEN中的一个或多个。所述第二阈值，优选的选取70%。
在寻找碳酸盐岩古岩溶洞穴储层时，高度为0.5~5.0m的小型溶洞储集体具有分布广、数量多的特点，作为我们储层预测的重点对象。地震反射和反演波阻抗特征仅对5m以上的储层有一定响应，且无法确定性识别；同时，虽然常规测井表现为低电阻、低密度、中子和声波时差增大等特征，但其探测范围有限。通过井震标定技术，将测井曲线作为“硬数据”，三维地震波形作为“软数据”，为大范围的储层地震资料解释建立测井约束，可以大大提高此类储层的识别精度。
首先井震标定整合测井与地震资料，根据补偿声波曲线AC与密度曲线DEN计算波阻抗曲线，继而计算反射系数曲线。以目的层段地震主频为依据构建雷克子波，合成地震记录，从标志层位置与地震反射同相轴两方面比对人工合成记录与井旁地震道，通过二者相关系数质控，最终得出时深转化关系。
在本实施例中，使用25HZ（目标层段地震数据主频）的雷克子波，根据补偿声波曲线、密度值曲线计算出的波阻抗信息得到合成地震记录。
通过观察地震波形数据确定鹰山组顶界标志层：与上覆巴楚组泥岩层相比，碳酸盐岩内幕区波形反射不规则，无一定方向，振幅可强可弱，同相轴可长可短连续性差；且具有非系统性同相轴反射终止和分叉现象。
在本实施例中，所述测井曲线数据与地震记录的时深转化关系为：
Figure 336889DEST_PATH_IMAGE136
其中，
Figure 809459DEST_PATH_IMAGE063
表示声波测井起始深度对应的地震信号时间，
Figure 391750DEST_PATH_IMAGE064
表示声波时差，
Figure 894144DEST_PATH_IMAGE066
表示 测井曲线数据采样间隔，T表示地震波双程旅行时。
单井井震标定图如图5所示；
根据钻录井、岩心数据确定目的层段古岩溶洞穴位置，依据充填程度划分为全充填、半充填储层与基岩；依据充填物岩性与组合关系划分出沉积充填、垮塌充填、化学充填及混合充填的储层类型。运用测井曲线直方图、交会图的统计方式，优选对划分储层类型敏感的测井曲线参数。例如波阻抗曲线反映了岩石物理性质的差异，可用于区分古岩溶洞穴与围岩；波阻抗-自然伽马交会图图版利用了地质体岩性、物性特征不仅可以有效区分全充填储层与半充填储层，而且可以划分储层充填物类型。波阻抗曲线柱状统计图如图6所示，波阻抗-自然GR交会图充填程度响应特征如图8所示。
步骤S600，构建等时格架模型：构建等时格架模型：基于高精度三维地震振幅数据体所反映的沉积地层规律和所述时深转化关系构建等时格架模型。确定目标层段位置后，根据地震数据剖面图选取同相轴连续、沉积环境稳定的界面或者标志层作为待预测范围的顶底界面。基于实际地质构造背景和所述时深转化关系构建等时格架模型；以T74标志层以下的等时格架模型为例，如图7所示。
步骤S700，井间储层参数模拟：基于各样本井的标准化测井曲线数据，确定能够反映整体地质情况的最优样本数参量，选取最优样本数参量个地震波形相关性最高的样本井的对储层敏感的特征参数数据构建初始模型，不断更正初始模型的参数，输出高精度特征值模拟结果数据体，所述高精度特征值模拟结果数据体为与所述对储层敏感的特征参数一一对应的数据体；
在本实施例中，所述步骤S700，具体包括具体包括步骤S710-步骤S790：
步骤S710，任选一样本井作为参考目标井，设定初始的样本数参量为1；
步骤S720，根据波形相似性原则选取数量为样本数参量的样本井标准化测井曲线特征参数数据与参考目标井标准化测井曲线特征参数数据进行相关性分析获得样本数参量-参考目标井特征参数的相关性值；
步骤S730，逐1增加样本数参量，重复步骤S720的方法获得各样本数参量对应的样本数参量-参考目标井曲线特征参数的相关性值，将所有样本数参量-参考目标井曲线特征参数的相关性值连接,获得参考井曲线特征参数相关性随样本数参量变化的相关性曲线；
步骤S740，选取另一样本井作为参考目标井，重复步骤S710-步骤S730的方法，获得多个参考井曲线特征参数相关性随样本数参量变化的相关性曲线，将所有参考井的特征参数的相关性随样本数参量变化的相关性曲线拟合为整体相关性曲线，选取所述整体相关性曲线中相关性随样本参数增加而升高最终保持平稳处的拐点，确定最优样本数参量；将所有参考井的特征参数的相关性随样本数参量变化的相关性曲线拟合为整体相关性曲线如图10所示；
本实施例中，在样本井中利用波形相似性和空间距离双变量优选低频结构相似的井作为空间估值样本，能够较好地反映空间结构的低频变化。
地震波形相似的两口井，表明所处沉积环境是相似的，其低频成分具有共性，可以增强反演结果低频段的确定性，同时约束了高频的取值范围，提高了反演与模拟结果的可靠性。
步骤S750，基于所述高精度的三维地震振幅数据体和等时格架模型，计算待测点位与样本井位的波形相关性，将所述波形相关性由大到小排序，选取最优样本数参量条地震波形相关性最高的样本井的对储层敏感的特征参数数据；基于所述相关性最高的样本井的地震波形特征数据，通过井间特征参数插值方式构建初始模型；
步骤S760，基于所述初始模型，选取最优样本数参量条地震波形关联度最高的样本井待模拟曲线参数对应的对储层敏感的特征参数作为先验信息；所述待模拟曲线参数为与所述对储层敏感的特征参数一一对应的曲线参数，至少包括波阻抗IMP，还可包括井径CAL、自然伽马GR、自然电位SP，电阻率曲线数据：深侧向测井RLLD、浅侧向测井RLLS、微侧向测井RLLM，物性表征曲线数据：补偿中子CNL、补偿声波AC、密度曲线DEN中的一个或多个；
步骤S770，将所述初始模型与先验信息进行匹配滤波，进而获得最大似然函数；
步骤S780，基于所述最大似然函数和先验信息，在贝叶斯框架下求得后验概率统计分布密度，对所述后验概率统计分布密度进行采样获得目标函数；
在本实施例中，所述步骤S780，具体包括：
步骤S781，利用白噪声满足高斯分布的规律，将高精度特征值模拟结果数据体的参数表示为：
Figure 384031DEST_PATH_IMAGE068
其中，Y表示测井曲线高精度特征值模拟结果数据体的参数，X表示待求解的地下地层实际特征参数值，N表示随机噪声；
步骤S782，因为
Figure 179949DEST_PATH_IMAGE069
也满足高斯分布，可以确定初始目标函数为：
Figure 526617DEST_PATH_IMAGE070
其中，
Figure 204723DEST_PATH_IMAGE071
表示与后验信息有关的函数，
Figure 701563DEST_PATH_IMAGE072
表示基于最优样本数选取样本井对样本 井的特征曲线进行匹配滤波后，求得后验概率统计分布密度，进而计算得到的特征参数期 望值，
Figure 961774DEST_PATH_IMAGE073
表示白噪声的协方差；
步骤S783，基于所述初始目标函数，通过最大后验估计，在目标函数中引入先验信息，获得稳定的目标函数为：
Figure 682606DEST_PATH_IMAGE137
其中，
Figure 785691DEST_PATH_IMAGE075
表示波阻抗或待模拟的特征参数，
Figure 679698DEST_PATH_IMAGE076
表示与地质、测井资料等先验信息有 关的函数，
Figure 246945DEST_PATH_IMAGE138
表示用于协调
Figure 341940DEST_PATH_IMAGE078
和
Figure 37714DEST_PATH_IMAGE076
之间的相互影响的平滑参数。
步骤S790，以所述目标函数作为所述初始模型的输入，通过马尔科夫链蒙特卡罗方法（Markov chain Monte Carlo method，MCMC）和Metropolis-Hastings抽样准则对后验概率分布抽样，不断更正初始模型的参数，选取目标函数取最大值时的解作为随机实现，取多次随机实现的均值作为期望值输出，将所述期望值输出作为高精度特征值模拟结果数据体；所述高精度特征值模拟结果数据体中的参数与所述对储层敏感的特征参数一一对应。
步骤S791，设M为目标空间，n为总样本数，m为马尔科夫链趋于平稳时的样本数；
步骤S792，预设一条马尔科夫链，使马尔科夫链收敛至平稳分布；
步骤S793，由M中的某一点
Figure 876357DEST_PATH_IMAGE079
出发，通过马尔科夫链进行抽样模拟，产生点序列：
Figure 298111DEST_PATH_IMAGE080
；
步骤S794，函数
Figure 423062DEST_PATH_IMAGE081
的期望估计为：
Figure 500739DEST_PATH_IMAGE139
其中，n表示产生的总样本数，m表示马尔科夫链达到平稳时的样本数，k表示累加参量；
步骤S795，选取一转移函数
Figure 674232DEST_PATH_IMAGE083
和初始值
Figure 763541DEST_PATH_IMAGE084
，若第i次迭代开始时的参数 值为
Figure 997077DEST_PATH_IMAGE085
，则第i次迭代过程为：
从
Figure 562050DEST_PATH_IMAGE083
中抽取一个备选值
Figure 601550DEST_PATH_IMAGE086
，计算备选值
Figure 732317DEST_PATH_IMAGE086
的接受概率
Figure 340016DEST_PATH_IMAGE087
：
Figure 500608DEST_PATH_IMAGE088
步骤S796，以
Figure 281482DEST_PATH_IMAGE089
，置
Figure 470018DEST_PATH_IMAGE090
，以概率
Figure 107673DEST_PATH_IMAGE091
，置
Figure 647239DEST_PATH_IMAGE092
；
步骤S797，不断扰动所述初始模型的参数，重复步骤S792-步骤S796的方法达到预 设的迭代次数n，获得后验样本
Figure 966224DEST_PATH_IMAGE093
，进而计算后验分布的各阶矩阵获得期望 输出值，将期望值输出作为高精度特征值模拟结果数据体。
步骤S800，刻画井间溶洞系统边界：基于所述样本井的对储层敏感的特征参数中的波阻抗曲线，绘制波阻抗曲线直方图，确定古岩溶洞穴与围岩波阻抗门槛值；
基于所述古岩溶洞穴与围岩波阻抗门槛值，将高精度特征值模拟结果数据体中超过所述阻抗门槛值的区域定义为碳酸盐岩基岩，将低于所述阻抗门槛值的区域定义为古岩溶洞穴型储层发育位置，获得古岩溶洞穴型储层结构特征和空间分布规律。可以得到2m以上的古岩溶洞穴型储层结构特征与空间分布规律。波阻抗曲线柱状统计图如图6所示。
步骤S900，刻画洞穴内部充填岩性物性边界：刻画洞穴内部充填岩性物性边界：基于所述对储层敏感的特征参数构建交会版图，根据所述个别深度段的岩性信息和物性信息在交汇版图中框选数据点作为测井解释样本点，获得充填程度与充填物类型的解释结论门槛值；
通过所述解释结论门槛值对所述高精度特征值模拟结果数据体进行交会分析获得洞穴内部充填程度与充填物组合三维空间形态特征；
在本实施例中，所述步骤S900，充填特征识别步骤包括：
步骤S910，以波阻抗值与自然GR值作为对储层敏感的特征参数构建交会图模板；
步骤S920，基于所述交会图模板，根据所述个别深度段岩性信息和物性信息，将样本点划分为全填充、半填充、基岩；
步骤S930，框选储集性能较好的半填充样本点，计算得到二维交会图解释结论门槛值；
步骤S940，将所述二维交会图解释结论门槛值输入三维自然GR特征参数模拟数据体与波阻抗反演数据体，刻画半填充储层的空间结构形态，获得古岩溶洞穴半充填部分三维空间形态特征。
在本实施例中，所述步骤S900，还包括获取其他古岩溶洞穴充填特征的步骤：
步骤S950，基于所述交会图模板，根据所述个别深度段岩性信息和物性信息，将样本点划分为基岩、沉积充填、化学充填、垮塌充填和混合充填，框选各类填充物样本点，分别计算对应解释结论门槛值，分别刻画对应的空间结构形态，获得古岩溶洞穴充填物三维空间形态特征。最终可得到得到古岩溶洞穴空间展布、充填物分布规律、充填程度所对应的三维空间形态。
古岩溶洞穴型储层半充填样本框选如图8所示，古岩溶洞穴型储层半充填储层刻画如图9所示；
步骤S1000，古岩溶洞穴结构与充填描述：基于所述古岩溶洞穴储层结构特征、空间分布规律、基于洞穴内部充填程度和充填物组合三维空间形态特征，采用岩性遮挡技术和三维雕刻技术雕刻出古岩溶洞穴空间展布和内部充填的发育特征。
本发明第二实施例的基于频谱分解和机器学习的岩溶储层充填分析系统，所述系统包括：：原始地球物理测井资料获取模块、地震数据获取模块、原始地球物理测井数据预处理模块、地震数据预处理模块、井震标定与特征参数选取模块、等时格架模型构建模块、井间特征参数模拟模块、井间溶洞系统边界刻画模块、溶洞内部充填岩性物性边界刻画模块、古岩溶洞穴结构与充填特征三维描述模块；
所述原始地球物理测井资料获取模块，配置为通过测井设备获取各样本井的原始测井数据，包括：通过测量电极测量各样本井的自然电位SP，通过自然伽马井下装置和自然伽马地面仪器测量各样本井的自然伽马GR，通过井径臂获取各样本井的井径CAL；通过传统测井设备获取电阻率曲线数据：深侧向测井RLLD、浅侧向测井RLLS和微侧向测井RLLM，物性表征曲线数据：补偿中子CNL、补偿声波曲线AC和密度曲线DEN；
基于成像测井信息、钻井信息、录井信息和岩心信息获得个别深度段的确定岩性信息和物性信息，进而确定目的层位标志层的深度数据；
所述地震数据获取模块，配置为通过地震波激发装置与接收装置获取原始地震波反射信号数据，并根据原始地震波反射信号数据的波形获取目的层位标志层的等时三维展布；
所述原始地球物理测井数据预处理模块，配置为基于所述样本井的原始测井数据，绘制测井曲线数据，并进行异常值处理和标准化处理，获得标准化测井曲线数据；
所述地震数据预处理模块，配置为基于所述地震波反射信号数据，通过混合相位子波反褶积和扩散滤波，获得高精度的三维地震振幅数据体；
所述井震标定与特征参数选取模块，配置为基于所述标准化测井曲线数据中的补偿声波曲线AC和密度曲线DEN获取样本井的波阻抗曲线，进而计算反射系数曲线，获取雷克子波的优选频率使其与高精度三维地震振幅数据体的主频保持一致，将所述雷克子波与所述反射系数曲线进行褶积运算，得到合成地震记录，并将所述目的层位标志层的深度数据与所述目的层位标志层的等时三维展布相比对进行井震标定，计算所述合成地震记录与井旁地震道波形的相关度，当相关度大于或等于预设的第一阈值时，判定井震标定结果合格，获得测井曲线数据与地震记录的时深转化关系和对储层敏感的特征参数；
所述对储层敏感的特征参数，其获得方法包括：
通过对井旁不同地质体产生的测井参数的响应绘制柱状统计图，当某一标准化测井曲线数据的数值可以将不同测井解释结论预设的第二阈值以上的数据点区分开时，则选定该标准化测井曲线数据为对储层敏感的特征参数；所述对储层敏感的特征参数至少包括波阻抗IMP，还可包括井径CAL、自然伽马GR、自然电位SP，电阻率曲线数据：深侧向测井RLLD、浅侧向测井RLLS、微侧向测井RLLM，物性表征曲线数据：补偿中子CNL、补偿声波曲线AC、密度曲线DEN中的一个或多个；
所述等时格架模型构建模块，配置为基于高精度三维地震振幅数据体所反映的沉积地层规律和所述时深转化关系构建等时格架模型；
所述井间特征参数模拟模块，配置为基于各样本井的标准化测井曲线数据，确定能够反映整体地质情况的最优样本数参量，选取最优样本数参量条地震波形相关性最高的样本井的对储层敏感的特征参数数据构建初始模型，不断更正初始模型的参数，输出高精度特征值模拟结果数据体，所述高精度特征值模拟结果数据体为与所述对储层敏感的特征参数一一对应的数据体；
所述井间溶洞系统边界刻画模块，配置为基于样本井对储层敏感的特征参数中的波阻抗曲线IMP，绘制波阻抗曲线直方图，确定区分古岩溶洞穴与围岩波阻抗门槛值；
基于所述古岩溶洞穴与围岩波阻抗门槛值，将高精度特征值模拟结果数据体中超过所述阻抗门槛值的区域定义为碳酸盐岩基岩，将低于所述阻抗门槛值的区域定义为古岩溶洞穴储层发育位置，获得古岩溶洞穴储层结构特征和空间分布规律；
所述溶洞内部充填岩性物性边界刻画模块，配置为刻画洞穴内部充填岩性物性边界：基于所述对储层敏感的特征参数构建交会版图，根据所述个别深度段的岩性信息和物性信息在交汇版图中框选数据点作为测井解释样本点，获得充填程度与充填物类型的解释结论门槛值；
通过所述解释结论门槛值对所述高精度特征值模拟结果数据体进行交会分析获得洞穴内部充填程度与充填物组合三维空间形态特征；
所述古岩溶洞穴结构与充填特征三维描述模块，配置为古岩溶洞穴结构与充填描述：基于所述古岩溶洞穴储层结构特征、空间分布规律、基于洞穴内部充填程度和充填物组合三维空间形态特征，采用岩性遮挡技术和三维雕刻技术雕刻出古岩溶洞穴空间展布和内部充填的发育特征。
所属技术领域的技术人员可以清楚地了解到，为描述的方便和简洁，上述描述的系统的具体工作过程及有关说明，可以参考前述方法实施例中的对应过程，在此不再赘述。
需要说明的是，上述实施例提供的基于频谱分解和机器学习的岩溶储层充填分析系统，仅以上述各功能模块的划分进行举例说明，在实际应用中，可以根据需要而将上述功能分配由不同的功能模块来完成，即将本发明实施例中的模块或者步骤再分解或者组合，例如，上述实施例的模块可以合并为一个模块，也可以进一步拆分成多个子模块，以完成以上描述的全部或者部分功能。对于本发明实施例中涉及的模块、步骤的名称，仅仅是为了区分各个模块或者步骤，不视为对本发明的不当限定。
本发明第三实施例的一种电子设备，包括：
至少一个处理器；以及
与至少一个所述处理器通信连接的存储器；其中，
所述存储器存储有可被所述处理器执行的指令，所述指令用于被所述处理器执行以实现上述的基于频谱分解和机器学习的岩溶储层充填分析方法。
本发明第四实施例的一种计算机可读存储介质，所述计算机可读存储介质存储有计算机指令，所述计算机指令用于被所述计算机执行以实现上述的基于频谱分解和机器学习的岩溶储层充填分析方法。
所属技术领域的技术人员可以清楚地了解到，为描述的方便和简洁，上述描述的存储装置、处理装置的具体工作过程及有关说明，可以参考前述方法实施例中的对应过程，在此不再赘述。
可以以一种或多种程序设计语言或其组合来编写用于执行本申请的操作的计算机程序代码，上述程序设计语言包括面向对象的程序设计语言—诸如Java、Smalltalk、C++，还包括常规的过程式程序设计语言—诸如“C”语言或类似的程序设计语言。程序代码可以完全地在用户计算机上执行、部分地在用户计算机上执行、作为一个独立的软件包执行、部分在用户计算机上部分在远程计算机上执行、或者完全在远程计算机或服务器上执行。在涉及远程计算机的情形中，远程计算机可以通过任意种类的网络——包括局域网(LAN)或广域网(WAN)—连接到用户计算机，或者，可以连接到外部计算机(例如利用因特网服务提供商来通过因特网连接)。
附图中的流程图和框图，图示了按照本申请各种实施例的系统、方法和计算机程序产品的可能实现的体系架构、功能和操作。在这点上，流程图或框图中的每个方框可以代表一个模块、程序段、或代码的一部分，该模块、程序段、或代码的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。也应当注意，在有些作为替换的实现中，方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如，两个接连地表示的方框实际上可以基本并行地执行，它们有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，框图和/或流程图中的每个方框、以及框图和/或流程图中的方框的组合，可以用执行规定的功能或操作的专用的基于硬件的系统来实现，或者可以用专用硬件与计算机指令的组合来实现。
术语“第一”、“第二”等是用于区别类似的对象，而不是用于描述或表示特定的顺序或先后次序。
术语“包括”或者任何其它类似用语旨在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备/装置不仅包括那些要素，而且还包括没有明确列出的其它要素，或者还包括这些过程、方法、物品或者设备/装置所固有的要素。
至此，已经结合附图所示的优选实施方式描述了本发明的技术方案，但是，本领域技术人员容易理解的是，本发明的保护范围显然不局限于这些具体实施方式。在不偏离本发明的原理的前提下，本领域技术人员可以对相关技术特征做出等同的更改或替换，这些更改或替换之后的技术方案都将落入本发明的保护范围之内。

Claims (10)

1.一种基于频谱分解和机器学习的岩溶储层充填分析方法，其特征在于，所述方法包括：
步骤S100，获取原始地球物理测井资料：通过测井设备获取各样本井的原始测井数据，包括：通过测量电极测量各样本井的自然电位SP，通过自然伽马井下装置和自然伽马地面仪器测量各样本井的自然伽马GR，通过井径臂获取各样本井的井径CAL；通过传统测井设备获取电阻率曲线数据：深侧向测井RLLD、浅侧向测井RLLS和微侧向测井RLLM，物性表征曲线数据：补偿中子CNL、补偿声波曲线AC和密度曲线DEN；
基于成像测井信息、钻井信息、录井信息和岩心信息获得个别深度段的确定岩性信息和物性信息，进而确定目的层位标志层的深度数据；
步骤S200，获取地震数据，通过地震波激发装置与接收装置获取原始地震波反射信号数据，并根据原始地震波反射信号数据的波形获取目的层位标志层的等时三维展布；
步骤S300，原始地球物理测井数据预处理：基于所述样本井的所有原始测井数据，绘制测井曲线数据，并进行异常值处理和标准化处理，获得标准化测井曲线数据；
步骤S400，地震数据预处理：基于所述地震波反射信号数据，通过混合相位子波反褶积和扩散滤波，获得高精度的三维地震振幅数据体；
步骤S500，井震标定与特征参数选取：基于所述标准化测井曲线数据中的补偿声波曲线AC和密度曲线DEN获取样本井的波阻抗曲线，进而计算反射系数曲线，获取雷克子波的优选频率使其与高精度三维地震振幅数据体的主频保持一致，将所述雷克子波与所述反射系数曲线进行褶积运算，得到合成地震记录，并将所述目的层位标志层的深度数据与所述目的层位标志层的等时三维展布相比对进行井震标定，计算所述合成地震记录与井旁地震道波形的相关度，当相关度大于或等于预设的第一阈值时，判定井震标定结果合格，获得测井曲线数据与地震记录的时深转化关系和对储层敏感的特征参数；
所述对储层敏感的特征参数，其获得方法包括：
通过对井旁不同地质体产生的测井参数的响应绘制柱状统计图，当某一标准化测井曲线数据的数值可以将不同测井解释结论预设的第二阈值以上的数据点区分开时，则选定该标准化测井曲线数据为对储层敏感的特征参数；所述对储层敏感的特征参数至少包括波阻抗IMP，还可包括井径CAL、自然伽马GR、自然电位SP，电阻率曲线数据：深侧向测井RLLD、浅侧向测井RLLS、微侧向测井RLLM，物性表征曲线数据：补偿中子CNL、补偿声波曲线AC、密度曲线DEN中的一个或多个；
步骤S600，构建等时格架模型：基于高精度三维地震振幅数据体所反映的沉积地层规律和所述时深转化关系构建等时格架模型；
步骤S700，井间储层参数模拟：基于各样本井的标准化测井曲线数据，确定能够反映整体地质情况的最优样本数参量，选取最优样本数参量条地震波形相关性最高的样本井的对储层敏感的特征参数数据构建初始模型，不断更正初始模型的参数，输出高精度特征值模拟结果数据体，所述高精度特征值模拟结果数据体为与所述对储层敏感的特征参数一一对应的数据体；
步骤S800，刻画井间溶洞系统边界：基于对储层敏感的特征参数中的样本井的波阻抗曲线IMP，绘制波阻抗曲线直方图，确定区分古岩溶洞穴与围岩波阻抗门槛值；
基于所述古岩溶洞穴与围岩波阻抗门槛值，将高精度特征值模拟结果数据体中超过所述区分古岩溶洞穴与围岩波阻抗门槛值的区域定义为碳酸盐岩基岩，将低于所述区分古岩溶洞穴与围岩波阻抗门槛值的区域定义为古岩溶洞穴储层发育位置，获得古岩溶洞穴储层结构特征和空间分布规律；
步骤S900，刻画洞穴内部充填岩性物性边界：刻画洞穴内部充填岩性物性边界：基于所述对储层敏感的特征参数构建交会版图，根据所述个别深度段的岩性信息和物性信息在交汇版图中框选数据点作为测井解释样本点，获得充填程度与充填物类型的解释结论门槛值；
通过所述解释结论门槛值对所述高精度特征值模拟结果数据体进行交会分析获得洞穴内部充填程度与充填物组合三维空间形态特征；
步骤S1000，古岩溶洞穴结构与充填描述：基于所述古岩溶洞穴储层结构特征、空间分布规律、基于洞穴内部充填程度和充填物组合三维空间形态特征，采用岩性遮挡技术和三维雕刻技术雕刻出古岩溶洞穴空间展布和内部充填的发育特征。
2.根据权利要求1所述的基于频谱分解和机器学习的岩溶储层充填分析方法，其特征在于，所述通过测量电极测量各样本井的自然电位SP，通过自然伽马井下装置和自然伽马地面仪器测量各样本井的自然伽马GR，具体包括：
所述通过测量电极测量各样本井的自然电位SP：
将测量电极N设置于地面，测量电机M通过电缆设置于井下；
沿井轴提升测量电极M测量自然电位随井深的变化；
自然电位值的计算方法为：
Figure 896450DEST_PATH_IMAGE001
其中，
Figure 548011DEST_PATH_IMAGE002
为总自然电位，
Figure 124486DEST_PATH_IMAGE003
为扩散电位系数，
Figure 254116DEST_PATH_IMAGE004
为扩散吸附电位系数，
Figure 68489DEST_PATH_IMAGE005
为泥浆滤液 电阻率值，
Figure 687689DEST_PATH_IMAGE006
为地层水电阻率值；
所述通过自然伽马井下装置和自然伽马地面仪器测量各样本井的自然伽马GR：
自然伽马井下装置包括探测器、放大器、高压电源；
通过探测器获取自然伽马射线，并将所述自然伽马射线转化为电脉冲信号，并通过放大器进行放大；
所述地面仪器把每分钟形成的电脉冲计数转化为电位差进行记录。
3.根据权利要求1所述的基于频谱分解和机器学习的岩溶储层充填分析方法，其特征在于，所述步骤S300，包括：
步骤S310，基于所述原始测井数据绘制原始测井曲线数据；
步骤S320，基于所述原始测井曲线数据，去除离群点获得去除离群点的测井曲线数据；
步骤S330，基于所述去除离群点的测井曲线数据，叠合工区内所有样本井位的单个测井曲线柱状图数据，通过整合阈值，获得标准化测井曲线数据。
4.根据权利要求1所述的基于频谱分解和机器学习的岩溶储层充填分析方法，其特征在于，所述步骤S400，具体包括：
步骤S410，基于所述地震波反射信号数据，将频率域地震记录褶积模型表示为：
Figure 689143DEST_PATH_IMAGE007
其中，
Figure 684781DEST_PATH_IMAGE008
表示经傅氏变换后的地震记录，
Figure 88080DEST_PATH_IMAGE009
表示经傅氏变换后的子波，
Figure 878182DEST_PATH_IMAGE010
表示 经傅氏变换后的反射系数的频谱，
Figure 632511DEST_PATH_IMAGE011
表示角频率；
步骤S420，将所述频率域地震记录褶积模型的等式两边取对数转化为线性系统，获得线性地震记录褶积模型：
Figure 103944DEST_PATH_IMAGE012
步骤S430，将所述线性地震记录褶积模型进行反傅氏变换，获得复赛谱序列：
Figure 424066DEST_PATH_IMAGE013
其中，
Figure 588332DEST_PATH_IMAGE014
表示地震波形记录的复赛谱序列，
Figure 626695DEST_PATH_IMAGE015
表示地震子波的复赛谱序列，
Figure 636239DEST_PATH_IMAGE016
表 示地层反射系数的复赛谱序列，
Figure 76448DEST_PATH_IMAGE017
表示地震波形记录时间；
步骤S440，基于所述复赛谱序列，通过低通滤波器进行子波与反射系数分离，提取子波振幅谱；
步骤S450，通过最小二乘法获取模拟地震子波振幅谱：
Figure 146035DEST_PATH_IMAGE018
其中，最小二乘法拟合参数
Figure 671694DEST_PATH_IMAGE019
为常数，
Figure 750509DEST_PATH_IMAGE020
表示子波振幅谱，
Figure 982907DEST_PATH_IMAGE021
和
Figure 551291DEST_PATH_IMAGE022
表示f的多 项式，f表示地震子波的频率；
步骤S460，基于所述模拟地震子波振幅谱，获得子波最大相位分量和最小相位分量；
设子波
Figure 236351DEST_PATH_IMAGE023
的最大相位分量为
Figure 915594DEST_PATH_IMAGE024
、最小相位分量为
Figure 2498DEST_PATH_IMAGE025
，则子波
Figure 945047DEST_PATH_IMAGE023
为：
Figure 179719DEST_PATH_IMAGE026
振幅谱的复赛谱中表示为：
Figure 334757DEST_PATH_IMAGE027
其中，振幅谱的复赛谱
Figure 604064DEST_PATH_IMAGE028
在复赛谱的正、负轴上对称显示，
Figure 717513DEST_PATH_IMAGE029
为地震子波最大相 位分量
Figure 377165DEST_PATH_IMAGE030
所对应的最小相位函数的复赛谱，
Figure 398211DEST_PATH_IMAGE031
为地震子波最小相位分量
Figure 459707DEST_PATH_IMAGE032
所对应 的最大相位函数的复赛谱；
步骤S470，基于所述振幅谱中的复赛谱确定一组具有相同振幅谱的混合相位子波集合，不断调整俞氏子波参数，保持低频、拓展高频和适当提高主频构建期望输出子波形态，在井曲线控制下以信噪比谱作参考寻找提高分辨率与保真度之间的最佳平衡点，获得整形后波形数据；
步骤S480，基于所述整形后波形数据，构建张量扩散模型：
Figure 478479DEST_PATH_IMAGE033
其中，
Figure 687743DEST_PATH_IMAGE034
表示扩散时间，
Figure 450163DEST_PATH_IMAGE035
表示散度算子，D表示扩散滤波器的张量型扩散系数，U表示 扩散滤波结果，
Figure 631746DEST_PATH_IMAGE036
表示
Figure 618156DEST_PATH_IMAGE034
=0时的扩散滤波结果，
Figure 252400DEST_PATH_IMAGE037
表示
Figure 615248DEST_PATH_IMAGE038
时刻的整形后波形数据，作 为张量扩散模型的初始条件，
Figure 916917DEST_PATH_IMAGE039
表示扩散滤波结果的梯度；
基于所述张量扩散模型构建梯度结构张量：
Figure 11912DEST_PATH_IMAGE040
其中，U表示扩散滤波结果，
Figure 195768DEST_PATH_IMAGE041
表示梯度向量张量积；
Figure 565570DEST_PATH_IMAGE042
表示尺度为
Figure 190586DEST_PATH_IMAGE043
的高斯函数：
Figure 784378DEST_PATH_IMAGE044
其中，r表示计算半径；
结构张量的特征向量为：
Figure 127635DEST_PATH_IMAGE045
其中，
Figure 301127DEST_PATH_IMAGE046
、
Figure 108546DEST_PATH_IMAGE047
和
Figure 810923DEST_PATH_IMAGE048
表示为梯度结构张量的3个特征向量，可视为局部正交坐标系，
Figure 907055DEST_PATH_IMAGE046
指向 地震信号的梯度方向，
Figure 415397DEST_PATH_IMAGE047
和
Figure 15006DEST_PATH_IMAGE048
组成的平面平行于地震信号的局部结构特征，
Figure 950601DEST_PATH_IMAGE049
、
Figure 268449DEST_PATH_IMAGE050
和
Figure 518165DEST_PATH_IMAGE051
为 分别与
Figure 300176DEST_PATH_IMAGE046
、
Figure 344356DEST_PATH_IMAGE047
和
Figure 211818DEST_PATH_IMAGE048
对应的三个特征值；
步骤S490，基于所述结构张量的特征向量分别计算线状结构置信度量、面状结构置信度量和扩散张量；
所述现状结构置信度量
Figure 265224DEST_PATH_IMAGE052
为：
Figure 839425DEST_PATH_IMAGE053
所述面状结构置信度量
Figure 116823DEST_PATH_IMAGE054
为：
Figure 409264DEST_PATH_IMAGE055
所述扩散张量D为：
Figure 782DEST_PATH_IMAGE056
其中，
Figure 491806DEST_PATH_IMAGE057
、
Figure 877788DEST_PATH_IMAGE058
和
Figure 923105DEST_PATH_IMAGE059
表示扩散张量的三个非负特征值，他们分别表示扩散滤波器沿
Figure 115052DEST_PATH_IMAGE046
、
Figure 398265DEST_PATH_IMAGE047
和
Figure 283045DEST_PATH_IMAGE048
这三个特征方向的滤波强度；
步骤S4100，重复步骤S480-S490的步骤，直至达到预设的迭代次数，获得扩散滤波结果，即为所述高精度的三维地震振幅数据体。
5.根据权利要求4所述的基于频谱分解和机器学习的岩溶储层充填分析方法，其特征在于，所述模拟地震子波振幅谱的计算方法为：
定位地震波反射信号数据中振幅谱的极大值和极大值所对应的频率；
将所述地震信号振幅谱的极大值和所述模拟地震子波振幅谱通过最小二乘法拟合的 方式得到参数
Figure 550078DEST_PATH_IMAGE060
和
Figure 483399DEST_PATH_IMAGE061
多项式的系数，进而获得拟合的极大值相应频率振幅值；
将所述地震信号振幅谱的极大值除以所拟合出的相应频率的振幅值，进而用商拟合多 项式
Figure 683436DEST_PATH_IMAGE022
的系数。
6.根据权利要求1所述的基于频谱分解和机器学习的岩溶储层充填分析方法，其特征在于，所述步骤S700，具体包括步骤S710-步骤S790：
步骤S710，任选一样本井作为参考目标井，设定初始的样本数参量为1；
步骤S720，根据波形相似性原则选取数量为样本数参量的样本井标准化测井曲线特征参数数据与参考目标井标准化测井曲线特征参数数据进行相关性分析获得样本数参量-参考目标井的特征参数相关性值；
步骤S730，逐1增加样本数参量，重复步骤S720的方法获得各样本数参量对应的样本数参量-参考目标井曲线特征参数的相关性值，将所有样本数参量-参考目标井曲线特征参数的相关性值连接,获得参考井曲线特征参数相关性随样本数参量变化的相关性曲线；
步骤S740，选取另一样本井作为参考目标井，重复步骤S710-步骤S730的方法，获得多个参考井曲线特征参数相关性随样本数参量变化的相关性曲线，将所有参考井的特征参数的相关性随样本数参量变化的相关性曲线拟合为整体相关性曲线，选取所述整体相关性曲线中相关性随样本参数增加而升高最终保持平稳处的拐点，确定最优样本数参量；
步骤S750，基于所述高精度的三维地震振幅数据体和等时格架模型，计算待测点位与样本井位的波形相关性，将所述波形相关性由大到小排序，选取最优样本数参量条地震波形相关性最高的样本井的对储层敏感的特征参数数据；基于所述相关性最高的样本井的地震波形特征数据所对应的样本井，通过井间特征参数插值方式构建初始模型；
步骤S760，基于所述初始模型，选取最优样本数参量条地震波形关联度最高的样本井待模拟曲线参数对应的对储层敏感的特征参数作为先验信息；所述待模拟曲线参数为与所述对储层敏感的特征参数一一对应的曲线参数，至少包括波阻抗IMP，还可包括井径CAL、自然伽马GR、自然电位SP，电阻率曲线数据：深侧向测井RLLD、浅侧向测井RLLS、微侧向测井RLLM，物性表征曲线数据：补偿中子CNL、补偿声波AC、密度曲线DEN中的一个或多个；
步骤S770，将所述初始模型与先验信息进行匹配滤波，进而获得最大似然函数；
步骤S780，基于所述最大似然函数和先验信息，在贝叶斯框架下求得后验概率统计分布密度，对所述后验概率统计分布密度进行采样获得目标函数；
步骤S790，以所述目标函数作为所述初始模型的输入，通过马尔科夫链蒙特卡罗方法MCMC和Metropolis-Hastings抽样准则对后验概率分布抽样，不断优化初始模型的参数，选取目标函数取最大值时的解作为随机实现，取多次随机实现的均值作为期望值输出，将所述期望值输出作为高精度特征值模拟结果数据体；所述高精度特征值模拟结果数据体中的参数与所述对储层敏感的特征参数一一对应。
7.根据权利要求6所述的基于频谱分解和机器学习的岩溶储层充填分析方法，其特征在于，所述步骤S780，具体包括：
步骤S781，利用白噪声满足高斯分布的规律，将高精度特征值模拟结果数据体的参数表示为：
Figure 676800DEST_PATH_IMAGE062
其中，Y表示测井曲线高精度特征值模拟结果数据体的参数，X表示待求解的地下地层实际特征参数值，N表示随机噪声；
步骤S782，因为
Figure 431129DEST_PATH_IMAGE063
也满足高斯分布，可以确定初始目标函数为：
Figure 964879DEST_PATH_IMAGE064
其中，
Figure 222685DEST_PATH_IMAGE065
表示与后验信息有关的函数，
Figure 449267DEST_PATH_IMAGE066
表示基于最优样本数选取样本井对样本井的 特征曲线进行匹配滤波后，求得后验概率统计分布密度，进而计算得到的特征参数期望值，
Figure 425313DEST_PATH_IMAGE067
表示白噪声的协方差；
步骤S783，基于所述初始目标函数，通过最大后验估计，在目标函数中引入先验信息，获得稳定的目标函数为：
Figure 700437DEST_PATH_IMAGE068
其中，
Figure 140645DEST_PATH_IMAGE069
表示波阻抗或待模拟的特征参数，
Figure 475812DEST_PATH_IMAGE070
表示与地质、测井资料等先验信息有关的 函数，
Figure 939154DEST_PATH_IMAGE071
表示用于协调
Figure 814706DEST_PATH_IMAGE072
和
Figure 47104DEST_PATH_IMAGE070
之间的相互影响的平滑参数。
8.根据权利要求6所述的基于频谱分解和机器学习的岩溶储层充填分析方法，其特征在于，所述步骤S790，具体步骤为：
步骤S791，设M为目标空间，n为总样本数，m为马尔科夫链趋于平稳时的样本数；
步骤S792，预设一条马尔科夫链，使马尔科夫链收敛至平稳分布；
步骤S793，由M中的某一点
Figure 553172DEST_PATH_IMAGE073
出发，通过马尔科夫链进行抽样模拟，产生点序列：
Figure 566127DEST_PATH_IMAGE074
；
步骤S794，函数
Figure 183053DEST_PATH_IMAGE075
的期望估计为：
Figure 535537DEST_PATH_IMAGE076
其中，n表示产生的总样本数，m表示马尔科夫链达到平稳时的样本数，k表示累加参量；
步骤S795，选取一转移函数
Figure 274823DEST_PATH_IMAGE077
和初始值
Figure 447179DEST_PATH_IMAGE078
，若第i次迭代开始时的参数值为
Figure 664533DEST_PATH_IMAGE079
，则第i次迭代过程为：
从
Figure 137103DEST_PATH_IMAGE077
中抽取一个备选值
Figure 719394DEST_PATH_IMAGE080
，计算备选值
Figure 441362DEST_PATH_IMAGE080
的接受概率
Figure 665670DEST_PATH_IMAGE081
：
Figure 727167DEST_PATH_IMAGE082
步骤S796，以
Figure 808256DEST_PATH_IMAGE083
，置
Figure 955203DEST_PATH_IMAGE084
，以概率
Figure 717623DEST_PATH_IMAGE085
，置
Figure 695943DEST_PATH_IMAGE086
；
步骤S797，不断扰动所述初始模型的参数，重复步骤C692-步骤C696的方法达到预设的 迭代次数n，获得后验样本
Figure 885616DEST_PATH_IMAGE087
，进而计算后验分布的各阶矩阵获得期望输出 值，将期望值输出作为高精度特征值模拟结果数据体。
9.根据权利要求1所述的基于频谱分解和机器学习的岩溶储层充填分析方法，其特征在于，所述步骤S900，包括：
步骤S910，以波阻抗值与自然GR值作为对储层敏感的特征参数构建交会图模板；
步骤S920，基于所述交会图模板，根据所述个别深度段岩性信息和物性信息，将样本点划分为全填充、半填充、基岩；
步骤S930，框选储集性能较好的半填充样本点，计算得到二维交会图解释结论门槛值；
步骤S940，将所述二维交会图解释结论门槛值输入三维自然GR特征参数模拟数据体与波阻抗反演数据体，刻画半填充储层的空间结构形态，获得古岩溶洞穴半充填部分三维空间形态特征；
所述步骤S900，还包括获取其他洞穴特征的步骤：
步骤S950，基于所述交会图模板，根据所述个别深度段岩性信息和物性信息，将样本点划分为基岩、沉积充填、化学充填、垮塌充填和混合充填，框选各类填充物样本点，分别计算对应解释结论门槛值，分别刻画对应的空间结构形态，获得古岩溶洞穴充填物三维空间形态特征。
10.一种基于频谱分解和机器学习的岩溶储层充填分析系统，其特征在于，所述系统包括：原始地球物理测井资料获取模块、地震数据获取模块、原始地球物理测井数据预处理模块、地震数据预处理模块、井震标定与特征参数选取模块、等时格架模型构建模块、井间特征参数模拟模块、井间溶洞系统边界刻画模块、溶洞内部充填岩性物性边界刻画模块、古岩溶洞穴结构与充填特征三维描述模块；
所述原始地球物理测井资料获取模块，配置为通过测井设备获取各样本井的原始测井数据，包括：通过测量电极测量各样本井的自然电位SP，通过自然伽马井下装置和自然伽马地面仪器测量各样本井的自然伽马GR，通过井径臂获取各样本井的井径CAL；通过传统测井设备获取电阻率曲线数据：深侧向测井RLLD、浅侧向测井RLLS和微侧向测井RLLM，物性表征曲线数据：补偿中子CNL、补偿声波曲线AC和密度曲线DEN;
基于成像测井信息、钻井信息、录井信息和岩心信息获得个别深度段的确定岩性信息和物性信息，进而确定目的层位标志层的深度数据；
所述地震数据获取模块，配置为通过地震波激发装置与接收装置获取原始地震波反射信号数据，并根据原始地震波反射信号数据的波形获取目的层位标志层的等时三维展布；
所述原始地球物理测井数据预处理模块，配置为基于所述样本井的原始测井数据，绘制测井曲线数据，并进行异常值处理和标准化处理，获得标准化测井曲线数据；
所述地震数据预处理模块，配置为基于所述地震波反射信号数据，通过混合相位子波反褶积和扩散滤波，获得高精度的三维地震振幅数据体；
所述井震标定与特征参数选取模块，配置为基于所述标准化测井曲线数据中的补偿声波曲线AC和密度曲线DEN获取样本井的波阻抗曲线，进而计算反射系数曲线，获取雷克子波的优选频率使其与高精度三维地震振幅数据体的主频保持一致，将所述雷克子波与所述反射系数曲线进行褶积运算，得到合成地震记录，并将所述目的层位标志层的深度数据与所述目的层位标志层的等时三维展布相比对进行井震标定，计算所述合成地震记录与井旁地震道波形的相关度，当相关度大于或等于预设的第一阈值时，判定井震标定结果合格，获得测井曲线数据与地震记录的时深转化关系和对储层敏感的特征参数；
所述对储层敏感的特征参数，其获得方法包括：
通过对井旁不同地质体产生的测井参数的响应绘制柱状统计图，当某一标准化测井曲线数据的数值可以将不同测井解释结论预设的第二阈值以上的数据点区分开时，则选定该标准化测井曲线数据为对储层敏感的特征参数；所述对储层敏感的特征参数至少包括波阻抗IMP，还可包括井径CAL、自然伽马GR、自然电位SP，电阻率曲线数据：深侧向测井RLLD、浅侧向测井RLLS、微侧向测井RLLM，物性表征曲线数据：补偿中子CNL、补偿声波曲线AC、密度曲线DEN中的一个或多个；
所述等时格架模型构建模块，配置为基于高精度三维地震振幅数据体所反映的沉积地层规律和所述时深转化关系构建等时格架模型；
所述井间特征参数模拟模块，配置为基于各样本井的标准化测井曲线数据，确定能够反映整体地质情况的最优样本数参量，选取最优样本数参量条地震波形相关性最高的样本井的对储层敏感的特征参数数据构建初始模型，不断更正初始模型的参数，输出高精度特征值模拟结果数据体，所述高精度特征值模拟结果数据体为与所述对储层敏感的特征参数一一对应的数据体；
所述井间溶洞系统边界刻画模块，配置为基于样本井对储层敏感的特征参数中的波阻抗曲线IMP，绘制波阻抗曲线直方图，确定区分古岩溶洞穴与围岩波阻抗门槛值；
基于所述古岩溶洞穴与围岩波阻抗门槛值，将高精度特征值模拟结果数据体中超过所述阻抗门槛值的区域定义为碳酸盐岩基岩，将低于所述阻抗门槛值的区域定义为古岩溶洞穴储层发育位置，获得古岩溶洞穴储层结构特征和空间分布规律；
所述溶洞内部充填岩性物性边界刻画模块，配置为刻画洞穴内部充填岩性物性边界：基于所述对储层敏感的特征参数构建交会版图，根据所述个别深度段的岩性信息和物性信息在交汇版图中框选数据点作为测井解释样本点，获得充填程度与充填物类型的解释结论门槛值；
通过所述解释结论门槛值对所述高精度特征值模拟结果数据体进行交会分析获得洞穴内部充填程度与充填物组合三维空间形态特征；
所述古岩溶洞穴结构与充填特征三维描述模块，配置为古岩溶洞穴结构与充填描述：基于所述古岩溶洞穴储层结构特征、空间分布规律、基于洞穴内部充填程度和充填物组合三维空间形态特征，采用岩性遮挡技术和三维雕刻技术雕刻出古岩溶洞穴空间展布和内部充填的发育特征。
CN202111066284.0A 2021-09-13 2021-09-13 基于频谱分解和机器学习的岩溶储层充填分析方法和系统 Active CN113759424B (zh)

Priority Applications (2)

Application Number Priority Date Filing Date Title
CN202111066284.0A CN113759424B (zh) 2021-09-13 2021-09-13 基于频谱分解和机器学习的岩溶储层充填分析方法和系统
US17/697,544 US11802985B2 (en) 2021-09-13 2022-03-17 Method and system for analyzing filling for karst reservoir based on spectrum decomposition and machine learning

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111066284.0A CN113759424B (zh) 2021-09-13 2021-09-13 基于频谱分解和机器学习的岩溶储层充填分析方法和系统

Publications (2)

Publication Number Publication Date
CN113759424A true CN113759424A (zh) 2021-12-07
CN113759424B CN113759424B (zh) 2022-03-08

Family

ID=78795056

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111066284.0A Active CN113759424B (zh) 2021-09-13 2021-09-13 基于频谱分解和机器学习的岩溶储层充填分析方法和系统

Country Status (2)

Country Link
US (1) US11802985B2 (zh)
CN (1) CN113759424B (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114966856A (zh) * 2022-08-02 2022-08-30 中国科学院地质与地球物理研究所 基于多频带地震资料的碳封存场址优选方法、系统和设备
CN115616665A (zh) * 2022-09-30 2023-01-17 中国科学院地质与地球物理研究所 卷积神经网络的处理方法、装置和电子设备
CN116122802A (zh) * 2022-12-23 2023-05-16 中国科学院地质与地球物理研究所 基于Unet双通道输出的钻测录井特征提取方法与系统
WO2023178632A1 (en) * 2022-03-25 2023-09-28 Saudi Arabian Oil Company Rock facies identification method based on seismic attribute classification using a machine learning network

Families Citing this family (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116227308B (zh) * 2023-05-09 2023-07-18 广东石油化工学院 一种浅层测井自然电场的数值模拟方法及系统
CN116522792A (zh) * 2023-05-17 2023-08-01 中国地质大学（武汉） 一种预测古岩溶暗河充填程度的方法
CN116699690A (zh) * 2023-06-14 2023-09-05 成都理工大学 一种基于裂缝精细刻画的钻井轨迹优化设计技术
CN116595396B (zh) * 2023-07-19 2023-09-15 广州海洋地质调查局三亚南海地质研究所 一种基于多窗口锚点的测井曲线标准化方法及装置
CN116976705B (zh) * 2023-09-19 2023-12-22 中国科学院地质与地球物理研究所 深地油气精准导航砂泥岩地层物性评价方法与系统
CN116975987B (zh) * 2023-09-20 2024-03-08 中国地质大学（北京） 基于声学特征的深水浅层岩土工程参数预测方法及装置
CN117077068B (zh) * 2023-10-18 2024-03-08 中国科学院地质与地球物理研究所 深地油气精准导航随钻声波测井数据实时标定方法与系统

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103529475A (zh) * 2013-04-19 2014-01-22 中国石油大学(华东) 一种识别和解释碳酸盐岩古岩溶储层三维结构的方法
US20190293818A1 (en) * 2017-04-28 2019-09-26 Pioneer Natural Resources Usa, Inc. High resolution seismic data derived from pre-stack inversion and machine learning
CN110308488A (zh) * 2018-03-27 2019-10-08 中国石油化工股份有限公司 确定洞穴充填程度的方法及系统
CN112698392A (zh) * 2020-11-20 2021-04-23 中国石油天然气股份有限公司 一种多层系储集体立体刻画与油气空间分布及定容方法
CN113050157A (zh) * 2020-10-15 2021-06-29 中国石油天然气股份有限公司 一种基于露头资料的碳酸盐岩地震储层反演方法及系统

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
RU2760102C2 (ru) * 2016-09-07 2021-11-22 Чайна Петролеум Энд Кемикал Корпорейшн Способ и система автоматического распознавания центра залежи в карстовой пещере

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103529475A (zh) * 2013-04-19 2014-01-22 中国石油大学(华东) 一种识别和解释碳酸盐岩古岩溶储层三维结构的方法
US20190293818A1 (en) * 2017-04-28 2019-09-26 Pioneer Natural Resources Usa, Inc. High resolution seismic data derived from pre-stack inversion and machine learning
CN110308488A (zh) * 2018-03-27 2019-10-08 中国石油化工股份有限公司 确定洞穴充填程度的方法及系统
CN113050157A (zh) * 2020-10-15 2021-06-29 中国石油天然气股份有限公司 一种基于露头资料的碳酸盐岩地震储层反演方法及系统
CN112698392A (zh) * 2020-11-20 2021-04-23 中国石油天然气股份有限公司 一种多层系储集体立体刻画与油气空间分布及定容方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
田飞: "塔河油田奥陶系缝洞型储层小型缝洞及其充填物测井识别", 《石油与天然气地质》 *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2023178632A1 (en) * 2022-03-25 2023-09-28 Saudi Arabian Oil Company Rock facies identification method based on seismic attribute classification using a machine learning network
CN114966856A (zh) * 2022-08-02 2022-08-30 中国科学院地质与地球物理研究所 基于多频带地震资料的碳封存场址优选方法、系统和设备
CN115616665A (zh) * 2022-09-30 2023-01-17 中国科学院地质与地球物理研究所 卷积神经网络的处理方法、装置和电子设备
CN116122802A (zh) * 2022-12-23 2023-05-16 中国科学院地质与地球物理研究所 基于Unet双通道输出的钻测录井特征提取方法与系统
CN116122802B (zh) * 2022-12-23 2023-07-14 中国科学院地质与地球物理研究所 基于Unet双通道输出的钻测录井特征提取方法与系统

Also Published As

Publication number Publication date
US20230083651A1 (en) 2023-03-16
CN113759424B (zh) 2022-03-08
US11802985B2 (en) 2023-10-31

Similar Documents

Publication Publication Date Title
CN113759425B (zh) 井震联合评价深层古岩溶储层充填特征的方法与系统
CN113759424B (zh) 基于频谱分解和机器学习的岩溶储层充填分析方法和系统
CN108931814B (zh) 一种基于多属性融合的基岩裂缝预测的方法
WO2020123084A1 (en) Machine learning-augmented geophysical inversion
US7133779B2 (en) Automated borehole geology and petrophysics interpretation using image logs
WO2017024702A1 (zh) 一种射线弹性参数的反演系统
CA2818790C (en) Seismic trace attribute
WO2008070596A1 (en) Identification of fracture clusters in rock formations
CA2437418C (en) Petrophysical property estimation using an acoustic calibration relationship
CN114994758B (zh) 碳酸盐岩断控储层的波阻抗提取与结构表征方法和系统
CN104502996A (zh) 一种密度曲线校正方法及系统
CN110646850B (zh) 隔夹层地震预测方法及装置
EA030770B1 (ru) Система и способ адаптивной сейсмической оптики
CN111077578B (zh) 岩层分布预测方法和装置
Mari et al. Contribution of seismic and acoustic methods to reservoir model building
CN115857047A (zh) 一种地震储层综合预测方法
CN113806674A (zh) 古河道纵向尺度的量化方法、装置、电子设备及存储介质
CN109283577A (zh) 一种地震层位标定方法
CN113514884A (zh) 一种致密砂岩储层预测方法
CN113848593A (zh) 一种定量预测含煤地层中岩浆岩侵蚀区的方法
CN112147676A (zh) 一种煤层及夹矸厚度预测方法
CN117250658B (zh) 建立研究区的地震数据集的方法
de Oliveira et al. Heterogeneity Index from Acoustic Image Logs and Its Application in Core Samples Representativeness: A Case Study in the Brazilian Pre-Salt Carbonates
CN115685344A (zh) 储层的确定方法、装置、存储介质及电子设备
CN117908155A (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