CN113759425A - 井震联合评价深层古岩溶储层充填特征的方法与系统 - Google Patents

井震联合评价深层古岩溶储层充填特征的方法与系统 Download PDF

Info

Publication number
CN113759425A
CN113759425A CN202111066355.7A CN202111066355A CN113759425A CN 113759425 A CN113759425 A CN 113759425A CN 202111066355 A CN202111066355 A CN 202111066355A CN 113759425 A CN113759425 A CN 113759425A
Authority
CN
China
Prior art keywords
data
seismic
well
logging
group
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
CN202111066355.7A
Other languages
English (en)
Other versions
CN113759425B (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 CN202111066355.7A priority Critical patent/CN113759425B/zh
Publication of CN113759425A publication Critical patent/CN113759425A/zh
Priority to US17/697,562 priority patent/US11500117B1/en
Application granted granted Critical
Publication of CN113759425B publication Critical patent/CN113759425B/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/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/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
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/282—Application of seismic models, synthetic seismograms
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • 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/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. for interpretation or for event detection
    • G01V1/30—Analysis
    • G01V1/306—Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/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. for interpretation or for event detection
    • G01V1/32—Transforming one recording into another or one representation into another
    • G01V1/325—Transforming one representation into another
    • 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/36—Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • G01V1/364—Seismic filtering
    • 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/46—Data acquisition
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V11/00—Prospecting or detecting by methods combining techniques covered by two or more of main groups G01V1/00 - G01V9/00
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V5/00—Prospecting or detecting by the use of ionising radiation, e.g. of natural or induced radioactivity
    • G01V5/04—Prospecting or detecting by the use of ionising radiation, e.g. of natural or induced radioactivity specially adapted for well-logging
    • G01V5/08—Prospecting or detecting by the use of ionising 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 ionising 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
    • 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/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/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
    • 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/64—Geostructures, e.g. in 3D data cubes
    • G01V2210/641—Continuity of geobodies
    • 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

Landscapes

  • Physics & Mathematics (AREA)
  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Geophysics (AREA)
  • Acoustics & Sound (AREA)
  • Environmental & Geological Engineering (AREA)
  • Geology (AREA)
  • High Energy & Nuclear Physics (AREA)
  • Geophysics And Detection Of Objects (AREA)

Abstract

本发明属于数据识别、记录载体的处理领域，具体涉及了一种井震联合评价深层古岩溶储层充填特征的方法与系统，旨在解决现有的石油勘探技术无法预测横向变化快的储层、无法识别大范围内复杂盆地碳酸盐岩洞穴型储层的发育特征的问题。本发明包括：获取标准化测井曲线数据；通过混合相位子波反褶积和扩散滤波，获得高精度的三维地震振幅数据体；通过PCA方法分析对所述对储层敏感的特征参数进行降维，获得第一组PCA数据；计算缝洞型储层结构特征；通过解释结论门槛值进行交会分析获得古岩溶洞穴空间展布和不同填充物类型的洞穴内部充填的发育特征。本发明达到识别大范围内复杂盆地碳酸盐岩岩溶洞穴型储层的发育特征的效果提高了刻画的精度。

Description

井震联合评价深层古岩溶储层充填特征的方法与系统
技术领域
本发明属于数据识别、记录载体的处理领域，具体涉及了一种井震联合评价深层古岩溶储层充填特征的方法与系统。
背景技术
含油气盆地深层（埋深超过4500m）蕴含丰富的油气资源，已经成为石油工业界勘探的重要领域。碳酸盐岩油藏占全球已探明石油储量的52%，产量的60%，是深层油气勘探的主要目标。
由于碳酸盐岩复杂的岩性、易变的成岩作用和随机的裂缝，古岩溶储层具有更强的非均质性。稀疏的测井曲线纵向分辨率高，仅能反映井筒附近有限深度的地层信息（一般小于10m），远离井筒的古岩溶储层仍存在很大变化。所以，那些浅层碎屑岩和多孔碳酸盐岩储层得到的探测技术和经验不能直接用于深层古岩溶储层的探测。然而，这往往被本领域技术人员所忽略，导致深层古岩溶储层的评价效果较差。
当前古岩溶储层充填特征评价存在以下问题：由于古岩溶洞穴型储层埋藏深，地震波能量吸收衰减作用强烈，导致地震资料主频显著降低（小于25Hz），分辨率大打折扣（一般为10m），深层时窗内信噪比过低。而钻井、测井和岩心样品表明，大多数洞穴高度不超过5米，现有的三维地震资料分辨率难以达到准确识别效果。此外，与测井曲线相比，地震属性所包含的岩性物性信息含量较低，无法有效识别古岩溶储层充填特征。
古岩溶储层研究的关键问题是如何准确认识深埋古岩溶系统的构造，如何利用地球物理资料有效地识别其位置，并且高精度地刻画古岩溶洞穴系统内部充填特征。
发明内容
为了解决目前深层古岩溶储层充填特征评价的技术难题，即现有的石油勘探洞穴型储层充填特征分析方法只将地质资料作为对随机模拟结论的约束，并未考虑地震信息对测井样本优选的贡献，并且采用全部的测井特征参数对储层进行描述数据过于繁杂，导致反演结果在高频段部分直接受控于井间插值，横向预测随机性强，对横向充填特征快速变化的深层古岩溶储层预测精度差的问题。本发明提供了一种井震联合评价深层古岩溶储层充填特征的方法，所述方法包括：
步骤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，获取第一组PCA数据：从所述对储层敏感的特征参数中优选预设条数的标准化测井曲线数据，并通过PCA方法对所述对储层敏感的特征参数分析，获得对储层敏感的第一组PCA数据；
步骤S700，构建等时格架模型：基于高精度三维地震振幅数据体所反映的沉积地层规律和所述时深转化关系构建等时格架模型；
步骤S800，井间储层参数模拟：基于各样本井的第一组PCA数据，确定能够反映整体地质情况的最优样本数参量，选取最优样本数参量个地震波形相关性最高的样本井的第一组PCA数据构建初始模型，不断更正初始模型的参数，输出第一组高精度特征值模拟结果数据体，所述第一组高精度特征值模拟结果数据体为与所述第一组PCA数据一一对应的数据体；
步骤S900，刻画井间溶洞系统边界：基于样本井的第一组PCA数据，进行交会分析，获得第一组各类地质体离散数据点分布图；
基于所述岩性信息和物性信息获取第一组各类地质体离散数据点分布图的解释结论，划分岩溶洞穴与围岩，进而构建第一组PCA数据交会版图，获取岩溶洞穴划分所需的第一组PCA数据门槛值；
基于所述第一组PCA数据门槛值，对所述第一组高精度特征值模拟结果数据体进行交会分析，获得岩溶洞穴三维空间形态特征；
步骤S1000，刻画洞穴内部充填岩性物性边界：将所述岩溶洞穴数据点所对应的对储层敏感的特征参数通过PCA方法进行分析，获得对充填物响应敏感的第二组PCA数据；
基于所述第二组PCA数据，通过如S800所述方法生成第二组高精度特征值模拟结果数据体；
基于所述对充填物响应敏感的第二组PCA数据，进行交会分析，获得第二组各类地质体离散数据点分布图；
基于所述岩性信息和物性信息获取第二组各类地质体离散数据点分布图的解释结论，区分填充物类型，进而构建第二组PCA数据交会版图，获取岩溶洞内部填充物类型划分所需的第二组PCA数据门槛值；
基于所述第二组PCA数据门槛值，对所述第二组高精度特征值模拟结果数据体进行交会分析，获得岩溶洞穴内部不同填充物类型的三维空间形态特征；
步骤S1100，古岩溶洞穴结构与充填描述：基于所述岩溶洞穴三维空间形态特征和岩溶洞穴内部不同填充物类型的三维空间形态特征，采用岩性遮挡技术和三维雕刻技术雕刻出古岩溶洞穴空间展布和不同填充物类型的内部充填的发育特征。
在一些优选的实施方式中，所述通过测量电极测量各样本井的自然电位SP，通过自然伽马井下装置和自然伽马地面仪器测量各样本井的自然伽马GR，具体包括：
所述通过测量电极测量各样本井的自然电位SP：
将测量电极N设置于地面，测量电机M通过电缆设置于井下；
沿井轴提升测量电极M测量自然电位随井深的变化；
自然电位值的计算方法为：
Figure DEST_PATH_IMAGE001
其中，
Figure 923962DEST_PATH_IMAGE002
为总自然电位，
Figure 496894DEST_PATH_IMAGE003
为扩散电位系数，
Figure 73369DEST_PATH_IMAGE004
为扩散吸附电位系数，
Figure 875103DEST_PATH_IMAGE005
为泥 浆滤液电阻率值，
Figure 423896DEST_PATH_IMAGE006
为地层水电阻率值；
所述通过自然伽马井下装置和自然伽马地面仪器测量各样本井的自然伽马GR：
自然伽马井下装置包括探测器、放大器、高压电源；
通过探测器获取自然伽马射线，并将所述自然伽马射线转化为电脉冲信号，并通过放大器进行放大；
所述地面仪器把每分钟形成的电脉冲计数转化为电位差进行记录。
在一些优选的实施方式中，所述PCA分析具体包括：
步骤B100，计算样本均值：
Figure 777517DEST_PATH_IMAGE007
其中，
Figure 703272DEST_PATH_IMAGE008
表示曲线总样本数，i表示样本计数参量，
Figure 433331DEST_PATH_IMAGE009
表示第i个样本特征参数数 据，
Figure 774314DEST_PATH_IMAGE010
表示样本均值；
步骤B200，对每条特征参数曲线进行归一化处理，获得归一化样本特征参数数据
Figure 298836DEST_PATH_IMAGE011
：
Figure 708957DEST_PATH_IMAGE012
其中，x表示样本特征参数采样点数据；
步骤B300，基于所述归一化样本特征参数数据
Figure 977128DEST_PATH_IMAGE011
，计算样本特征参数协方差
Figure 907038DEST_PATH_IMAGE013
：
Figure 71303DEST_PATH_IMAGE014
其中，
Figure 844087DEST_PATH_IMAGE015
表示第c个X样本特征参数数据，
Figure 775002DEST_PATH_IMAGE016
表示第c个Y样本特征参数数据；
步骤B400，基于所述样本特征参数协方差，构建样本特征参数协方差矩阵：
Figure DEST_PATH_IMAGE017
Z表示样本特征参数数据；
包括p种样本特征参数数据的协方差矩阵为：
Figure 621736DEST_PATH_IMAGE018
其中，
Figure 691323DEST_PATH_IMAGE019
表示多个样本特征参数数据，
Figure 951403DEST_PATH_IMAGE020
表示第1个样本特征参数数据与第p种样本特征参数数据的协方差；
步骤B500，基于所述包括p种样本特征参数数据的协方差矩阵，获取该协方差矩阵 的特征值
Figure DEST_PATH_IMAGE021
和p个特征向量
Figure 482747DEST_PATH_IMAGE022
，
Figure 387249DEST_PATH_IMAGE023
，
Figure 893317DEST_PATH_IMAGE024
；
基于所述协方差矩阵的特征值和特征向量计算包括p种特征值的PCA数据：
Figure 375114DEST_PATH_IMAGE025
其中，
Figure 910482DEST_PATH_IMAGE026
表示第j个PCA数据向量，
Figure 59704DEST_PATH_IMAGE027
表示第j个样本特征参数数据向量，
Figure 408776DEST_PATH_IMAGE028
表示第j 个样本特征参数数据的特征向量，j
Figure 581132DEST_PATH_IMAGE029
1,2,
Figure 532907DEST_PATH_IMAGE030
p,为计数参量；
步骤B600，基于包括p种特征值的PCA数据，将对应的特征值从大到小进行排序，得 到从大到小排序的特征值
Figure 926848DEST_PATH_IMAGE031
，…，
Figure 571456DEST_PATH_IMAGE032
计算方差累计贡献量：
Figure 903212DEST_PATH_IMAGE033
其中，
Figure 861940DEST_PATH_IMAGE034
表示第j个从大到小排序的特征值，
Figure 454596DEST_PATH_IMAGE035
为累计前k个协方差累计贡献量，
Figure 660318DEST_PATH_IMAGE036
；
当k从1开始不断增加时，选取
Figure 604003DEST_PATH_IMAGE035
首次大于预设的第三阈值时对应的前k个协方差 对应的PCA数据向量，即为本轮PCA分析获得的PCA数据。
在一些优选的实施方式中，所述步骤S400，具体包括：
步骤S410，基于所述地震波反射信号数据，将频率域地震记录褶积模型表示为：
Figure 366423DEST_PATH_IMAGE037
其中，
Figure 954530DEST_PATH_IMAGE038
表示经傅氏变换后的地震记录，
Figure DEST_PATH_IMAGE039
表示经傅氏变换后的子波，
Figure 596733DEST_PATH_IMAGE040
表示经傅氏变换后的反射系数的频谱，
Figure 230977DEST_PATH_IMAGE041
表示角频率；
步骤S420，将所述频率域地震记录褶积模型的等式两边取对数转化为线性系统，获得线性地震记录褶积模型：
Figure 328246DEST_PATH_IMAGE042
步骤S430，将所述线性地震记录褶积模型进行反傅氏变换，获得复赛谱序列：
Figure 302018DEST_PATH_IMAGE043
其中，
Figure 397013DEST_PATH_IMAGE044
表示地震波形记录的复赛谱序列，
Figure 315290DEST_PATH_IMAGE045
表示地震子波的复赛谱序列，
Figure 343814DEST_PATH_IMAGE046
表示地层反射系数的复赛谱序列，
Figure 31147DEST_PATH_IMAGE047
表示地震波形记录时间；
步骤S440，基于所述复赛谱序列，通过低通滤波器进行子波与反射系数分离，提取子波振幅谱；
步骤S450，通过最小二乘法获取模拟地震子波振幅谱：
Figure 562622DEST_PATH_IMAGE048
其中，最小二乘法拟合参数
Figure 577983DEST_PATH_IMAGE049
为常数，
Figure 548213DEST_PATH_IMAGE050
表示子波振幅谱，
Figure 214687DEST_PATH_IMAGE051
和
Figure 917063DEST_PATH_IMAGE052
表示f 的多项式，f表示地震子波的频率；
步骤S460，基于所述模拟地震子波振幅谱，获得子波最大相位分量和最小相位分量；
设子波
Figure 544354DEST_PATH_IMAGE053
的最大相位分量为
Figure 928062DEST_PATH_IMAGE054
、最小相位分量为
Figure 324408DEST_PATH_IMAGE055
，则子波
Figure 197686DEST_PATH_IMAGE053
为：
Figure 702486DEST_PATH_IMAGE056
振幅谱的复赛谱中表示为：
Figure 483360DEST_PATH_IMAGE057
其中，振幅谱的复赛谱
Figure 203054DEST_PATH_IMAGE058
在复赛谱的正、负轴上对称显示，
Figure 919337DEST_PATH_IMAGE059
为地震子波最 大相位分量
Figure 786799DEST_PATH_IMAGE060
所对应的最小相位函数的复赛谱，
Figure 495998DEST_PATH_IMAGE061
为地震子波最小相位分量
Figure 70199DEST_PATH_IMAGE062
所 对应的最大相位函数的复赛谱；
步骤S470，基于所述振幅谱中的复赛谱确定一组具有相同振幅谱的混合相位子波集合，不断调整俞氏子波参数，保持低频、拓展高频和适当提高主频构建期望输出子波形态，在井曲线控制下以信噪比谱作参考寻找提高分辨率与保真度之间的最佳平衡点，获得整形后波形数据；
步骤S480，基于所述整形后波形数据，构建张量扩散模型：
Figure 82017DEST_PATH_IMAGE063
其中，
Figure 46562DEST_PATH_IMAGE064
表示扩散时间，
Figure 434818DEST_PATH_IMAGE065
表示散度算子，D表示扩散滤波器的张量型扩散系数，U 表示扩散滤波结果，
Figure 863526DEST_PATH_IMAGE066
表示
Figure 191387DEST_PATH_IMAGE064
=0时的扩散滤波结果，
Figure DEST_PATH_IMAGE067
表示
Figure 439966DEST_PATH_IMAGE068
时刻的整形后波形数 据，作为张量扩散模型的初始条件，
Figure DEST_PATH_IMAGE069
表示扩散滤波结果的梯度；
基于所述张量扩散模型构建梯度结构张量：
Figure 100754DEST_PATH_IMAGE070
其中，U表示扩散滤波结果，
Figure 570919DEST_PATH_IMAGE071
表示梯度向量张量积；
Figure 393381DEST_PATH_IMAGE072
表示尺度为
Figure 191573DEST_PATH_IMAGE073
的高斯函数为：
Figure 796998DEST_PATH_IMAGE074
其中，r表示计算半径，
结构张量的特征向量为：
Figure 997035DEST_PATH_IMAGE075
其中，
Figure 990399DEST_PATH_IMAGE076
、
Figure 400520DEST_PATH_IMAGE077
和
Figure 934270DEST_PATH_IMAGE078
表示为梯度结构张量的3个特征向量，可视为局部正交坐标系，
Figure 864180DEST_PATH_IMAGE076
指向地震信号的梯度方向，
Figure 28445DEST_PATH_IMAGE077
和
Figure 535649DEST_PATH_IMAGE078
组成的平面平行于地震信号的局部结构特征，
Figure 997724DEST_PATH_IMAGE079
、
Figure 375615DEST_PATH_IMAGE080
和
Figure 241940DEST_PATH_IMAGE081
为分别与
Figure 642966DEST_PATH_IMAGE076
、
Figure 252939DEST_PATH_IMAGE077
和
Figure 485337DEST_PATH_IMAGE078
对应的三个特征值；
步骤S490，基于所述结构张量的特征向量分别计算线状结构置信度量、面状结构置信度量和扩散张量；
所述现状结构置信度量
Figure 650126DEST_PATH_IMAGE082
为：
Figure DEST_PATH_IMAGE083
所述面状结构置信度量
Figure 69606DEST_PATH_IMAGE084
为：
Figure 686532DEST_PATH_IMAGE085
所述扩散张量D为：
Figure 835754DEST_PATH_IMAGE086
其中，
Figure 434095DEST_PATH_IMAGE087
、
Figure 606450DEST_PATH_IMAGE088
和
Figure 558225DEST_PATH_IMAGE089
表示扩散张量的三个非负特征值，他们分别表示扩散滤波器沿
Figure 702899DEST_PATH_IMAGE076
、
Figure 737720DEST_PATH_IMAGE077
和
Figure 194109DEST_PATH_IMAGE078
这三个特征方向的滤波强度；
步骤S4100，重复步骤S480-S490的步骤，直至达到预设的迭代次数，获得扩散滤波结果，即为所述高精度的三维地震振幅数据体。
在一些优选的实施方式中，所述模拟地震子波振幅谱的计算方法为：
定位地震波反射信号数据中振幅谱的极大值和极大值所对应的频率；
将所述地震信号振幅谱的极大值和所述模拟地震子波振幅谱通过最小二乘法拟 合的方式得到参数
Figure 152838DEST_PATH_IMAGE090
和
Figure 152018DEST_PATH_IMAGE091
多项式的系数，进而获得拟合的极大值相应频率振幅值；
将所述地震信号振幅谱的极大值除以所拟合出的相应频率的振幅值，进而用商拟 合多项式
Figure 967527DEST_PATH_IMAGE052
的系数。
在一些优选的实施方式中，所述测井曲线数据与地震记录的时深转化关系的时深转化关系为：
Figure 114475DEST_PATH_IMAGE092
其中，
Figure 532686DEST_PATH_IMAGE093
表示声波测井起始深度对应的地震信号时间，
Figure 511007DEST_PATH_IMAGE094
表示声波时差，
Figure 372784DEST_PATH_IMAGE095
表示表 示每个采样点之间的时间运算序号，
Figure 7027DEST_PATH_IMAGE096
表示测井曲线数据采样间隔，
Figure 104296DEST_PATH_IMAGE097
表示地震波双程旅 行时。
在一些优选的实施方式中，所述步骤S800，具体包括，具体包括步骤S810-步骤S890：
步骤S810，任选一样本井作为参考目标井，设定初始的样本数参量为1；
步骤S820，根据波形相似性原则选取数量为样本数参量的样本井第一组PCA数据与参考目标井标第一组PCA数据进行相关性分析获得样本数参量-参考目标井的第一组PCA数据相关性值；
步骤S830，逐1增加样本数参量，重复步骤S720的方法获得各样本数参量对应的样本数参量-参考目标井的第一组PCA数据相关性值，将所有样本数参量-参考目标井的第一组PCA数据相关性值连接，获得参考井第一组PCA数据相关性随样本数参量变化的相关性曲线；
步骤S840，选取另一样本井作为参考目标井，重复步骤S810-步骤S830的方法，获得多个参考井第一组PCA数据相关性随样本数参量变化的相关性曲线，将所有参考井的第一组PCA数据波形相关性随样本数参量变化的相关性曲线拟合为整体相关性曲线，选取所述整体相关性曲线中相关性随样本数参量增加而升高最终保持平稳处的拐点，确定最优样本数参量；
步骤S850，基于所述高精度的三维地震振幅数据体和等时格架模型，计算待测点位与样本井位的波形相关性，将所述波形相关性由大到小排序，选取最优样本数参量条地震波形相关性最高的样本井的第一组PCA数据；基于所述相关性最高的样本井的地震波形特征数据所对应的样本井，通过井间特征参数插值方式构建初始模型；
步骤S860，基于所述初始模型，选取最优样本数参量条地震波形关联度最高的样本井第一组PCA数据作为先验信息；
步骤S870，将所述初始模型与先验信息进行匹配滤波，进而获得最大似然函数；
步骤S880，基于所述最大似然函数和先验信息，在贝叶斯框架下求得后验概率统计分布密度，对所述后验概率统计分布密度进行采样获得目标函数；
步骤S890，以所述目标函数作为所述初始模型的输入，通过马尔科夫链蒙特卡罗方法（Markov chain Monte Carlo method，MCMC）和Metropolis-Hastings抽样准则对后验概率分布抽样，不断优化初始模型的参数，选取目标函数取最大值时的解作为随机实现，取多次随机实现的均值作为期望值输出，将所述期望值输出作为高精度特征值模拟结果数据体；所述高精度特征值模拟结果数据体中的参数与所述第一组PCA数据一一对应。
在一些优选的实施方式中，所述步骤S880，具体包括：
步骤S881，利用白噪声满足高斯分布的规律，将高精度特征值模拟结果数据体的参数表示为：
Figure 58827DEST_PATH_IMAGE098
其中，Y表示测井曲线高精度特征值模拟结果数据体的参数，X表示待求解的地下地层实际特征参数值，N表示随机噪声；
步骤S882，因为
Figure DEST_PATH_IMAGE099
也满足高斯分布，可以确定初始目标函数为：
Figure 950560DEST_PATH_IMAGE100
其中，
Figure 9783DEST_PATH_IMAGE101
表示与后验信息有关的函数，
Figure 645163DEST_PATH_IMAGE102
表示基于最优样本数选取样本井对样本 井的特征曲线进行匹配滤波后，求得后验概率统计分布密度，进而计算得到的特征参数期 望值，
Figure 535759DEST_PATH_IMAGE103
表示白噪声的协方差；
步骤S883，基于所述初始目标函数，通过最大后验估计，在目标函数中引入先验信息，获得稳定的目标函数为：
Figure 988606DEST_PATH_IMAGE104
其中，
Figure 394179DEST_PATH_IMAGE105
表示待模拟的特征参数，
Figure 974196DEST_PATH_IMAGE106
表示与地质、测井资料等先验信息有关的函数，
Figure 719299DEST_PATH_IMAGE107
表示用于协调
Figure 218413DEST_PATH_IMAGE108
和
Figure 235916DEST_PATH_IMAGE106
之间的相互影响的平滑参数。
在一些优选的实施方式中，所述步骤S890，具体步骤为：
步骤S891，设M为目标空间，n为总样本数，m为马尔科夫链趋于平稳时的样本数；
步骤S892，预设一条马尔科夫链，使马尔科夫链收敛至平稳分布；
步骤S893，由M中的某一点
Figure 416362DEST_PATH_IMAGE109
出发，通过马尔科夫链进行抽样模拟，产生点序列：
Figure 78288DEST_PATH_IMAGE110
；
步骤S894，函数
Figure DEST_PATH_IMAGE111
的期望估计为：
Figure 420407DEST_PATH_IMAGE112
其中，n表示产生的总样本数，m表示马尔科夫链达到平稳时的样本数，k表示累加参量；
步骤S895，选取一转移函数
Figure 738256DEST_PATH_IMAGE113
和初始值
Figure 643764DEST_PATH_IMAGE114
，若第i次迭代开始时的参数 值为
Figure 425775DEST_PATH_IMAGE115
，则第i次迭代过程为：
从
Figure 142058DEST_PATH_IMAGE113
中抽取一个备选值
Figure 947203DEST_PATH_IMAGE116
，计算备选值
Figure 797348DEST_PATH_IMAGE116
的接受概率
Figure 295850DEST_PATH_IMAGE117
：
Figure 510930DEST_PATH_IMAGE118
步骤S896，以
Figure 334530DEST_PATH_IMAGE119
，置
Figure 863731DEST_PATH_IMAGE120
，以概率
Figure DEST_PATH_IMAGE121
，置
Figure 823597DEST_PATH_IMAGE122
；
步骤S897，不断扰动所述初始模型的参数，重复步骤S892-步骤S896的方法达到预 设的迭代次数n，获得后验样本
Figure 396530DEST_PATH_IMAGE123
，进而计算后验分布的各阶矩阵获得期望 输出值，将期望值输出作为第一组高精度特征值模拟结果数据体。
在一些优选的实施方式中，所述步骤S300，包括：
步骤S310，基于所述原始测井数据绘制原始测井曲线数据；
步骤S320，基于所述原始测井曲线数据，去除离群点获得去除离群点的测井曲线数据；
步骤S330，基于所述去除离群点的测井曲线数据，叠合工区内所有样本井位的单个测井曲线柱状图数据，通过整合阈值，获得标准化测井曲线数据。
本发明的另一方面，提出了一种井震联合评价深层古岩溶储层充填特征的系统，所述系统包括：原始地球物理测井资料获取模块、地震数据获取模块、原始地球物理测井数据预处理模块、地震数据预处理模块、井震标定与特征参数选取模块、第一组PCA数据获取模块、等时格架模型构建模块、井间储层参数模拟模块、刻画井间溶洞系统边界模块、刻画洞穴内部充填岩性物性边界模块和古岩溶洞穴结构与充填描述模块；
所述原始地球物理测井资料获取模块，配置为通过测井设备获取各样本井的原始测井数据，包括：通过测量电极测量各样本井的自然电位SP，通过自然伽马井下装置和自然伽马地面仪器测量各样本井的自然伽马GR，通过井径臂获取各样本井的井径CAL；通过传统测井设备获取电阻率曲线数据：深侧向测井RLLD、浅侧向测井RLLS和微侧向测井RLLM，物性表征曲线数据：补偿中子CNL、补偿声波曲线AC和密度曲线DEN；
基于成像测井信息、钻井信息、录井信息和岩心信息获得个别深度段的确定岩性信息和物性信息，进而确定目的层位标志层的深度数据；
所述地震数据获取模块，配置为通过地震波激发装置与接收装置获取原始地震波反射信号数据，并根据原始地震波反射信号数据的波形获取目的层位标志层的等时三维展布；
所述原始地球物理测井数据预处理模块，配置为基于所述样本井的原始测井数据，绘制测井曲线数据，并进行异常值处理和标准化处理，获得标准化测井曲线数据；
所述地震数据预处理模块，配置为基于所述地震波反射信号数据，通过混合相位子波反褶积和扩散滤波，获得高精度三维地震振幅数据体；
所述井震标定与特征参数选取模块，配置为井震标定与特征参数选取：基于所述标准化测井曲线数据中的补偿声波曲线和密度曲线DEN获取样本井的波阻抗曲线，进而计算反射系数曲线，获取雷克子波的优选频率使其与高精度三维地震振幅数据体的主频保持一致，将所述雷克子波与所述反射系数曲线进行褶积运算，得到合成地震记录，并将所述目的层位标志层的深度数据与所述目的层位标志层的等时三维展布相比对进行井震标定，计算所述合成地震记录与井旁地震道波形的相关度，当相关度大于或等于预设的第一阈值时，判定井震标定结果合格，获得测井曲线数据与地震记录的时深转化关系和对储层敏感的特征参数；
所述对储层敏感的特征参数，其获得方法包括：
通过对井旁不同地质体产生的测井参数的响应绘制柱状统计图，当某一标准化测井曲线数据可以将不同测井解释结论预设的第二阈值以上的数据点区分开时，则选定该标准化测井曲线数据为对储层敏感的特征参数；所述对储层敏感的特征参数至少包括波阻抗IMP，还可包括井径CAL、自然伽马GR、自然电位SP，电阻率曲线数据：深侧向测井RLLD、浅侧向测井RLLS、微侧向测井RLLM，物性表征曲线数据：补偿中子CNL、补偿声波曲线AC、密度曲线DEN中的一个或多个；
所述第一组PCA数据获取模块，配置为从所述对储层敏感的特征参数中优选预设条数的的标准化测井曲线数据，并通过PCA分析对所述对储层敏感的特征参数进行降维，获得第一组PCA数据；
所述等时格架模型构建模块，配置为基于高精度三维地震振幅数据体所反映的沉积地层规律和所述时深转化关系构建等时格架模型；
所述井间储层参数模拟模块，配置为基于各样本井的第一组PCA数据，确定能够反映整体地质情况的最优样本数参量，选取最优样本数参量个地震波形相关性最高的样本井的第一组PCA数据构建初始模型，不断更正初始模型的参数，输出第一组高精度特征值模拟结果数据体，所述第一组高精度特征值模拟结果数据体为与所述第一组PCA数据一一对应的数据体；
所述井间溶洞系统边界刻画模块，配置为基于样本井的第一组PCA数据，进行交会分析，获得第一组各类地质体离散数据点分布图；
基于所述岩性信息和物性信息获取第一组各类地质体离散数据点分布图的解释结论，划分岩溶洞穴与围岩，进而构建第一组PCA数据交会版图，获取岩溶洞穴划分所需的第一组PCA数据门槛值；
基于所述第一组PCA数据门槛值，对所述第一组高精度特征值模拟结果数据体进行交会分析，获得岩溶洞穴三维空间形态特征；
所述洞穴内部充填岩性物性边界刻画模块，配置为刻画洞穴内部充填岩性物性边界：所述岩溶洞穴数据点，将所述岩溶洞穴数据点所对应的对储层敏感的特征参数通过PCA方法进行分析，获得对充填物响应敏感的第二组PCA数据；
基于所述第二组PCA数据，通过如井间储层参数模拟模块所述功能生成第二组高精度特征值模拟结果数据体；
基于所述对充填物响应敏感的第二组PCA数据，进行交会分析，获得第二组各类地质体离散数据点分布图；
基于所述岩性信息和物性信息获取第二组各类地质体离散数据点分布图的解释结论，区分填充物类型，进而构建第二组PCA数据交会版图，获取岩溶洞内部填充物类型划分所需的第二组PCA数据门槛值；
基于所述第二组PCA数据门槛值，对所述第二组高精度特征值模拟结果数据体进行交会分析，获得岩溶洞穴内部不同填充物类型的三维空间形态特征；
所述古岩溶洞穴结构与充填描述模块，配置为：基于所述岩溶洞穴三维空间形态特征和岩溶洞穴内部不同填充物类型的三维空间形态特征，采用岩性遮挡技术和三维雕刻技术雕刻出古岩溶洞穴空间展布和不同填充物类型的内部充填的发育特征。
本发明的第三方面，提出了一种电子设备，包括：
至少一个处理器；
以及至少一个所述处理器通信连接的存储器；其中，
所述存储器存储有可被所述处理器执行的指令，所述指令用于被所述处理器执行以实现上述的井震联合评价深层古岩溶储层充填特征的方法。
本发明的第四方面，提出了一种计算机可读存储介质，所述计算机可读存储介质存储有计算机指令，所述计算机指令用于被所述计算机执行以实现上述的井震联合评价深层古岩溶储层充填特征的方法。
本发明的有益效果：
本发明充分利用测井和地震资料，采用机器学习算法，以波形指示模拟方法为核心，利用主成分分析算法(PCA算法)建立敏感参数，实现测井资料对储层结构与充填的精确识别；然后利用敏感参数和三维地震数据建立三维地震模拟数据体，实现高分辨率的测井解释结果外推至井间三维大范围区域，最终评价深层古岩溶储层的充填特征。
本发明通过利用高精度的三维地震波形数据之间的关联性在初始等时地质格架模型基础上井震联合进行三维参数模拟，并运用两次PCA分析（PCA）对多条经预处理后的测井曲线变换坐标系，将多条常规测井曲线数据所反映的岩性、物性信息降维分析，使得少量的PCA数据具有多条测井曲线对储层的表征能力。
本发明还通过基于PCA数据所构建的二维交会图很好地区分岩溶洞穴及其内部充填的数据点，进而获得岩溶储层识别交会图版，将该图版所体现的储层划分门槛值信息运用至三维参数模拟资料中，得到PCA交会分析约束下的岩溶洞穴纵横向结构及其内部充填特征刻画结论。提高了刻画的精度，能够预测横向变化快的储层。为大范围评价碳酸盐岩古岩溶洞穴储集性能及油气运移提供可靠的理论依据。
本发明在贝叶斯理论指导下的马尔科夫链抽样准则与蒙特卡罗估计法建立地球物理探测方法之间的映射关系，以达到识别大范围内复杂盆地碳酸盐岩古岩溶洞穴型储层的发育特征的效果。
附图说明
通过阅读参照以下附图所作的对非限制性实施例所作的详细描述，本申请的其它特征、目的和优点将会变得更明显：
图1是井震联合评价深层古岩溶储层充填特征的方法与系统实施例的流程示意图；
图2是本发明去除离群点的测井曲线示意图；
图3是本发明进行标准化流程中叠合所有曲线的示意图；
图4是本发明实施例中单井井震标定图；
图5是本发明实施例中以T74标志层以下的等时格架模型图；
图6是本发明实施例中将所有参考井的波形相关性随样本数参量变化的相关性曲线拟合为整体相关性曲线示意图PCA；
图7是本发明实施例中第一组PCA数据交会分析抓取岩溶洞穴数据点并框选门槛值示意图；
图8是本发明实施例中古岩溶洞穴三维分布与结构刻画剖面图；
图9是本发明实施例中第二组PCA数据交会分析划分岩溶洞穴内部充填物并框选门槛值示意图；
图10是本发明实施例中古岩溶洞穴内部充填物发育特征刻画剖面。
具体实施方式
下面结合附图和实施例对本申请作进一步的详细说明。可以理解的是，此处所描述的具体实施例仅用于解释相关发明，而非对该发明的限定。另外还需要说明的是，为了便于描述，附图中仅示出了与有关发明相关的部分。
需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。下面将参考附图并结合实施例来详细说明本申请。
本发明充分利用测井和地震资料，采用机器学习算法，以波形指示模拟方法为核心，利用主成分分析算法(PCA算法)建立敏感参数，实现测井资料对储层结构与充填的精确识别；然后利用敏感参数和三维地震数据建立三维地震模拟数据体，实现高分辨率的测井解释结果外推至井间三维大范围区域，最终评价深层古岩溶储层的充填特征。
为了更清晰地对本发明古岩溶洞穴储层充填特征智能划分方法进行说明，下面结合图1对本发明实施例中各步骤展开详述。
本发明第一实施例的井震联合评价深层古岩溶储层充填特征的方法与系统，包括步骤S100-步骤S1100，各步骤详细描述如下：
经实际地震数据统计分析表明，古岩溶洞穴所产生的地震波形特征通常由多组波谷和波峰组成，经研究发现，这种特征性反射是由地震波干涉造成的。并且反射波的振幅与洞穴充填物组合有关。Xu等（2016）通过物理模型模拟同样可以发现不同形状、厚度、组合关系的古岩溶洞穴型目标体会产生多样的地震波形反射特征，即反射形态特征对于洞穴直径、洞穴形态与分布规律有关。学者们进一步发现在等时地层格架的同一相带内，波形特征的相似性可以代表具有一定相关性的岩性组合。因此，基于波形指示的思路进行等时界面分区反演与特征参数模拟，利用了地震波形的横向变化信息，更好地体现了沉积环境的约束，更符合沉积地质规律。
步骤S100，获取原始地球物理测井资料：通过测井设备获取各样本井的原始测井数据，包括：通过测量电极测量各样本井的自然电位SP，通过自然伽马井下装置和自然伽马地面仪器测量各样本井的自然伽马GR，通过井径臂获取各样本井的井径CAL；通过传统测井设备获取电阻率曲线数据：深侧向测井RLLD、浅侧向测井RLLS和微侧向测井RLLM，物性表征曲线数据：补偿中子CNL、补偿声波曲线AC和密度曲线DEN；
基于成像测井信息、钻井信息、录井信息和岩心信息获得个别深度段的确定岩性信息和物性信息，进而确定目的层位标志层的深度数据；
所述通过测量电极测量各样本井的自然电位SP，通过自然伽马井下装置和自然伽马地面仪器测量各样本井的自然伽马GR，具体包括：
所述通过测量电极测量各样本井的自然电位SP：
将测量电极N设置于地面，测量电机M通过电缆设置于井下；
沿井轴提升测量电极M测量自然电位随井深的变化；
自然电位值的计算方法为：
Figure 973005DEST_PATH_IMAGE124
其中，
Figure 774739DEST_PATH_IMAGE125
为总自然电位，
Figure 57952DEST_PATH_IMAGE126
为扩散电位系数，
Figure 677153DEST_PATH_IMAGE127
为扩散吸附电位系数，
Figure 865557DEST_PATH_IMAGE128
为泥 浆滤液电阻率值，
Figure 330037DEST_PATH_IMAGE129
为地层水电阻率值；
所述通过自然伽马井下装置和自然伽马地面仪器测量各样本井的自然伽马GR：
自然伽马井下装置包括探测器、放大器、高压电源；
通过探测器获取自然伽马射线，并将所述自然伽马射线转化为电脉冲信号，并通过放大器进行放大；
所述地面仪器把每分钟形成的电脉冲计数转化为电位差进行记录。
双侧向电极系由主电极、监督电极、环状屏蔽电极与柱状电极组成。主电极
Figure 733336DEST_PATH_IMAGE130
居 中，上下对称分布监督电极
Figure 398804DEST_PATH_IMAGE131
、
Figure 684292DEST_PATH_IMAGE132
和
Figure 342675DEST_PATH_IMAGE133
、
Figure 600481DEST_PATH_IMAGE134
，以及环状屏蔽电极
Figure DEST_PATH_IMAGE135
、
Figure 968008DEST_PATH_IMAGE136
，在
Figure 944055DEST_PATH_IMAGE135
、
Figure 15916DEST_PATH_IMAGE136
的外侧对 称位置上加了两个柱状电极。深侧向电极系中两个柱状电极是屏蔽电极
Figure 312249DEST_PATH_IMAGE137
、
Figure 444153DEST_PATH_IMAGE138
；浅侧向电 极系中两个柱状电极是回路电极
Figure 641917DEST_PATH_IMAGE139
和
Figure 392835DEST_PATH_IMAGE140
。在电极系较远处装有对比电极N和深侧向电极系 的回路电极B；微球形聚焦电极系的主电极
Figure 687550DEST_PATH_IMAGE130
是矩形片状电极，依次向外的矩形框电极是测 量电极
Figure 928038DEST_PATH_IMAGE141
、辅助电极
Figure 800048DEST_PATH_IMAGE135
、监督电极
Figure 479291DEST_PATH_IMAGE142
和
Figure 238300DEST_PATH_IMAGE133
，回路电极B装在仪器外壳上或极板支撑架上。测 量任一监督电极与对比电极之间的电位差变化即反映介质电阻率的变化，其视电阻率表达 式为：
Figure 649690DEST_PATH_IMAGE143
式中，
Figure 618783DEST_PATH_IMAGE019
为电极系系数，深侧向测井、浅侧向测井与微球型聚焦测井各不相同；
Figure 226351DEST_PATH_IMAGE144
为监督电极M的电位；
Figure 964499DEST_PATH_IMAGE145
为主电流。
井下仪器的发射换能器晶体振动，引起周围介质的质点发生振动，产生向井内泥 浆及岩层中传播的声波。在井中可以用接收换能器R、R2先后接收到滑行波，记录时差
Figure 812370DEST_PATH_IMAGE146
，进 而测量地层的声波速度。到达
Figure 144125DEST_PATH_IMAGE147
和
Figure 165171DEST_PATH_IMAGE148
的时刻分别为
Figure 148039DEST_PATH_IMAGE149
和
Figure 166811DEST_PATH_IMAGE150
，那么到达两个接收换能器的时 间差
Figure 844917DEST_PATH_IMAGE146
为:
Figure 545019DEST_PATH_IMAGE151
式中，
Figure 461023DEST_PATH_IMAGE076
为泥浆声速；
Figure 447433DEST_PATH_IMAGE077
为地层的声速
密度测井仪器包括有一个伽马源，两个接收伽马射线的探测器，即长源距探测器和短源距探测器。它们安装在滑板上，测井时被推靠到井壁上。在下井仪器的上方装有辅助电子线路。通常用137Cs作伽马源，它发射的伽马射线具有中等能量(0.661MeV)，用它照射物质只能产生康普顿散射和光电效应。地层的密度不同，则对伽马光子的散射和吸收的能力不同，探测器接收到的伽马光子的计数率也就不同。已知通过距离为L的伽马光子的计数率N为：
Figure 5978DEST_PATH_IMAGE152
若只存在康普顿散射，则
Figure 368826DEST_PATH_IMAGE153
即为康普顿散射吸收系数，变换得：
Figure 811440DEST_PATH_IMAGE154
式中，C为与地层有关的常数，L为接收源距离伽马源的距离。
源距选定后，对仪器进行刻度标定，找到
Figure 172014DEST_PATH_IMAGE155
和N的这种关系，则记录散射伽马光子 计数率N就可以测得地层的密度
Figure 90292DEST_PATH_IMAGE155
补偿中子测井由一个中子发射源与两个接收探测器组成。通过对两个探测器计 数，得到两个计数率
Figure 115885DEST_PATH_IMAGE156
和
Figure 6481DEST_PATH_IMAGE157
，取其比值可以反映地层的氢含量：
Figure 334694DEST_PATH_IMAGE158
式中，
Figure 615634DEST_PATH_IMAGE159
与
Figure 54706DEST_PATH_IMAGE160
为两个探测器到中子源的距离；
Figure 799808DEST_PATH_IMAGE161
为减速长度，反映了地层含氢量。
井径测量时，将仪器下人井预计的深度上，然后通过传统的方式打开井径臂，于是互成90°的四个井径腿便在弹簀力的作用下向外伸张，其末端紧贴井壁。随着仪器的向上提升，井径臂就会由于井径的变化而发生张缩，并带动连杆作上下运动。连杆与一个电位器的滑动端相连，于是井径的变化便可转换成电阻的变化。当给该可动电阻上通以一-定强度的电流时，可动电阻的某一固定端与滑动端之间的电位差将随着其间电阻值的变化而变化。于是测量这一电位差，便可反映井径大小：
Figure 689135DEST_PATH_IMAGE162
式中，
Figure 316426DEST_PATH_IMAGE163
为井径测量仪器电位差；
Figure 762450DEST_PATH_IMAGE164
为初始值；c为仪器常数。
步骤S200，获取地震数据，通过地震波激发装置与接收装置获取原始地震波反射信号数据，并根据原始地震波反射信号数据的波形获取目的层位标志层的等时三维展布；
在本实施例中，在工区井眼处探测深度范围为5500~5750m的九条常规测井曲线数据，采样间隔设置为0.01m；运用三维地震勘探方法，借助地震波激发源与地震信号检波器，得到工区面积约27km2的三维地震资料数据，信号记录双程旅行时为4s，采样点时间间隔为1ms，探测深度超过6000m。
步骤S300，原始地球物理测井数据预处理：基于所述样本井的所有原始测井数据，绘制测井曲线数据，并进行异常值处理和标准化处理，获得标准化测井曲线数据；
在本实施例中，所述步骤S300，包括：
步骤S310，基于所述原始测井数据绘制原始测井曲线数据；
步骤S320，基于所述原始测井曲线数据，去除离群点获得去除离群点的测井曲线数据；离群点，即在数据集中存在不合理的值，统计所有井数据的单个曲线参数直方图，合理调整保留区间阈值进行离群点去除，在本实施例中优选的去除与中位数偏离最大的前百分之5的数据，去除离群点的测井曲线直方图如图2所示；
步骤S330，基于所述去除离群点的测井曲线数据，叠合工区内所有样本井位的单个测井曲线柱状图数据，通过整合阈值，获得标准化测井曲线数据；以补偿声波曲线AC标准化处理为例，获得标准化测井曲线数据。由于仪器差异或受到其他因素影响，导致不同的井之间常规测井数据会出现整体偏大或整体偏小，需要将曲线进行标准化，叠合所有曲线的如图3所示。
步骤S300，地震数据预处理：基于所述地震波反射信号数据，通过混合相位子波反褶积和扩散滤波，获得高精度三维地震振幅数据体；
目前25m×25m道间距的三维地震测网技术被广泛应用于石油勘探领域，通常在垂向上以2ms的采样间隔接收地震波反射信号，采样总时长在6s内，以探测不同深度层段的地质特征。本发明常用于2m以上的目标体探测，对地震资料主频要求较高，应在50~60Hz范围内。当古岩溶洞穴发育层段内地震资料振幅数据体主频低于50HZ时，则需要采用混合相位子波反褶积（mixed-phase wavelet estimation and maximum posteriorideconvolution）与扩散滤波（enhancing diffusion filtering）进行三维地震数据的拓频降噪处理，得到高分辨率、高信噪比的三维地震数据体。
混合相位子波反褶积是一种在确保处理后的地震数据具有较高保真度高的前提下，通过拓宽有效频带从而提高地震信号分辨率的数据处理方法，相当于本实施例中的S410-步骤S1100。
在本实施例中，步骤S400，具体包括：
步骤S410，基于所述地震波反射信号数据，将频率域地震记录褶积模型表示为：
Figure 299742DEST_PATH_IMAGE165
其中，
Figure 235337DEST_PATH_IMAGE038
表示经傅氏变换后的地震记录，
Figure 208978DEST_PATH_IMAGE039
表示经傅氏变换后的子波，
Figure 458694DEST_PATH_IMAGE040
表示经傅氏变换后的反射系数的频谱，
Figure 975126DEST_PATH_IMAGE041
表示角频率；
步骤S420，将所述频率域地震记录褶积模型的等式两边取对数转化为线性系统，获得线性地震记录褶积模型：
Figure 956988DEST_PATH_IMAGE042
步骤S430，将所述线性地震记录褶积模型进行反傅氏变换，获得复赛谱序列：
Figure 496554DEST_PATH_IMAGE043
其中，
Figure 81119DEST_PATH_IMAGE044
表示地震波形记录的复赛谱序列，
Figure 573762DEST_PATH_IMAGE045
表示地震子波的复赛谱序列，
Figure 585580DEST_PATH_IMAGE046
表示地层反射系数的复赛谱序列，
Figure 815704DEST_PATH_IMAGE047
表示地震波形记录时间；
步骤S440，基于所述复赛谱序列，通过低通滤波器进行子波与反射系数分离，提取子波振幅谱；子波与反射系数序列平滑程度的差异在复赛谱中易于区分：子波能量出现于原点附近，而反射系数序列则远离原点。利用低通滤波器就可以将复赛谱中的子波与反射系数分离达到子波振幅谱提取的目的。
步骤S450，通过最小二乘法获取模拟地震子波振幅谱：
Figure 141644DEST_PATH_IMAGE048
其中，最小二乘法拟合参数
Figure 632668DEST_PATH_IMAGE049
为常数，
Figure 205600DEST_PATH_IMAGE050
表示子波振幅谱，
Figure 516496DEST_PATH_IMAGE051
和
Figure 318230DEST_PATH_IMAGE052
表示f 的多项式，f表示地震子波的频率；
在本实施例中，所述模拟地震子波振幅谱的计算方法为：
定位地震波反射信号数据中振幅谱的极大值和极大值所对应的频率；
将所述地震信号振幅谱的极大值和所述模拟地震子波振幅谱通过最小二乘法拟 合的方式得到参数
Figure 867023DEST_PATH_IMAGE090
和
Figure 220644DEST_PATH_IMAGE091
多项式的系数，进而获得拟合的极大值相应频率振幅值；
将所述地震信号振幅谱的极大值除以所拟合出的相应频率的振幅值，进而用商拟 合多项式
Figure 143469DEST_PATH_IMAGE052
的系数；
步骤S460，基于所述模拟地震子波振幅谱，获得子波最大相位分量和最小相位分量；
设子波
Figure 139107DEST_PATH_IMAGE053
的最大相位分量为
Figure 214511DEST_PATH_IMAGE054
、最小相位分量为
Figure 942295DEST_PATH_IMAGE055
，则子波
Figure 227783DEST_PATH_IMAGE053
为：
Figure 151746DEST_PATH_IMAGE166
振幅谱的复赛谱中表示为：
Figure 206289DEST_PATH_IMAGE057
其中，振幅谱的复赛谱
Figure 104975DEST_PATH_IMAGE058
在复赛谱的正、负轴上对称显示，
Figure 753125DEST_PATH_IMAGE059
为地震子波最 大相位分量
Figure 90566DEST_PATH_IMAGE060
所对应的最小相位函数的复赛谱，
Figure 127179DEST_PATH_IMAGE061
为地震子波最小相位分量
Figure 462346DEST_PATH_IMAGE062
所 对应的最大相位函数的复赛谱；
步骤S470，基于所述振幅谱中的复赛谱确定一组具有相同振幅谱的混合相位子波集合，不断调整俞氏子波参数，保持低频、拓展高频和适当提高主频构建期望输出子波形态，在井曲线控制下以信噪比谱作参考寻找提高分辨率与保真度之间的最佳平衡点，获得整形后波形数据；子波振幅谱如图4所示；
经过混合相位子波反褶积后，拓展了地震数据的有效频带，高频部分得到了合理加强。在地震波形上表现为同相轴数目增多，更容易反映地震波反射信息的细节变化，并且在振幅、相位、频率方面改善了同一反射波组波形的一致性。在古岩溶洞穴地震响应上，“串珠状”反射特征尤为明显，且串珠内部形态的细节能够清晰地显示，代表了不同结构特征与充填物组合的复杂古岩溶洞穴型储层地震反射，有助于后期精细地质解释。
Fhemers和Hocker于2003年首次将扩散滤波技术应用在地震资料的处理解释中。该技术不仅可以有效压制噪声，而且可以尽可能保留地震数据中的细节：如地质体边缘、断层、不整合面、尖灭等，为后续的地震解释和储层预测工作提供了可靠的地震资料，大大提高油气勘探开发的成功率。
为达到衰减地震噪声，并增强地质结构特征的扩散滤波效果，寻求扩散张量是该方法最关键的步骤：
步骤S480，基于所述整形后波形数据，构建张量扩散模型：
Figure 456847DEST_PATH_IMAGE167
其中，
Figure 473344DEST_PATH_IMAGE064
表示扩散时间，
Figure 440163DEST_PATH_IMAGE065
表示散度算子，D表示扩散滤波器的张量型扩散系数，U 表示扩散滤波结果，
Figure 742968DEST_PATH_IMAGE066
表示
Figure 614978DEST_PATH_IMAGE064
=0时的扩散滤波结果，
Figure 28642DEST_PATH_IMAGE067
表示
Figure 115547DEST_PATH_IMAGE068
时刻的整形后波形数 据，作为张量扩散模型的初始条件，
Figure 464620DEST_PATH_IMAGE069
表示扩散滤波结果的梯度；
基于所述张量扩散模型构建梯度结构张量：
Figure 699292DEST_PATH_IMAGE168
其中，U表示扩散滤波结果，
Figure 41281DEST_PATH_IMAGE071
表示梯度向量张量积；
Figure 982692DEST_PATH_IMAGE072
表示尺度为
Figure 892879DEST_PATH_IMAGE073
的高斯函数为：
Figure 224634DEST_PATH_IMAGE074
其中，r表示计算半径；
结构张量的特征向量为：
Figure 183363DEST_PATH_IMAGE169
其中，
Figure 41598DEST_PATH_IMAGE076
、
Figure 981741DEST_PATH_IMAGE077
和
Figure 925426DEST_PATH_IMAGE078
表示为梯度结构张量的3个特征向量，可视为局部正交坐标系，
Figure 687845DEST_PATH_IMAGE076
指向地震信号的梯度方向，
Figure 275953DEST_PATH_IMAGE077
和
Figure 527943DEST_PATH_IMAGE078
组成的平面平行于地震信号的局部结构特征，
Figure 838486DEST_PATH_IMAGE079
、
Figure 139018DEST_PATH_IMAGE080
和
Figure 971844DEST_PATH_IMAGE081
为分别与
Figure 270102DEST_PATH_IMAGE076
、
Figure 922800DEST_PATH_IMAGE077
和
Figure 27022DEST_PATH_IMAGE078
对应的三个特征值；
步骤S490，基于所述结构张量的特征向量分别计算线状结构置信度量、面状结构置信度量和扩散张量；
所述现状结构置信度量
Figure 104568DEST_PATH_IMAGE082
为：
Figure 432782DEST_PATH_IMAGE170
所述面状结构置信度量
Figure 448142DEST_PATH_IMAGE084
为：
Figure 90476DEST_PATH_IMAGE171
所述扩散张量D为：
Figure 897895DEST_PATH_IMAGE172
其中，
Figure 521643DEST_PATH_IMAGE087
、
Figure 352196DEST_PATH_IMAGE088
和
Figure 594959DEST_PATH_IMAGE089
表示扩散张量的三个非负特征值，他们分别表示扩散滤波器沿
Figure 132250DEST_PATH_IMAGE076
、
Figure 802266DEST_PATH_IMAGE077
和
Figure 120115DEST_PATH_IMAGE078
这三个特征方向的滤波强度；
步骤S4100，重复步骤S480-S490的步骤，直至达到预设的迭代次数，获得扩散滤波结果，即为所述高精度的三维地震振幅数据体。扩散滤波算法保留了“串珠状”反射古岩溶洞穴型储层地质特征，增强了地震数据对目标地质体的成像能力。同时起到了压制噪声、提高同相轴横向连续性与地震信号的信噪比的效果。
步骤S500，井震标定与特征参数选取：基于所述标准化测井曲线数据中的补偿声波曲线和密度曲线DEN获取样本井的波阻抗曲线，进而计算反射系数曲线，获取雷克子波的优选频率使其与高精度三维地震振幅数据体的主频保持一致，将所述雷克子波与所述反射系数曲线进行褶积运算，得到合成地震记录，并将所述目的层位标志层的深度数据与所述目的层位标志层的等时三维展布相比对进行井震标定，计算所述合成地震记录与井旁地震道波形的相关度，当相关度大于或等于预设的第一阈值时，判定井震标定结果合格，获得测井曲线数据与地震记录的时深转化关系和对储层敏感的特征参数；所述第一阈值，优选的选取85%；
所述对储层敏感的特征参数，其获得方法包括：
通过对井旁不同地质体产生的测井参数的响应绘制柱状统计图，当某一标准化测井曲线数据的数值可以将不同测井解释结论预设的第二阈值以上的数据点区分开时，则选定该标准化测井曲线数据为对储层敏感的特征参数；所述对储层敏感的特征参数至少包括波阻抗IMP，还可包括井径CAL、自然伽马GR、自然电位SP，电阻率曲线数据：深侧向测井RLLD、浅侧向测井RLLS、微侧向测井RLLM，物性表征曲线数据：补偿中子CNL、补偿声波曲线AC、密度曲线DEN中的一个或多个。所述第二阈值，优选的选取70%。
在寻找碳酸盐岩古岩溶洞穴储层时，高度为0.5~5.0m的小型溶洞储集体具有分布广、数量多的特点，作为我们储层预测的重点对象。地震反射和反演波阻抗特征仅对5m以上的储层有一定响应，且无法确定性识别；同时，虽然常规测井表现为低电阻、低密度、中子和声波时差增大等特征，但其探测范围有限。通过井震标定技术，将测井曲线作为“硬数据”，三维地震波形作为“软数据”，为大范围的储层地震资料解释建立测井约束，可以大大提高此类储层的识别精度。
首先井震标定整合测井与地震资料，根据补偿声波曲线AC与密度曲线DEN计算波阻抗曲线，继而计算反射系数曲线。以目的层段地震主频为依据构建雷克子波，合成地震记录，从标志层位置与地震反射同相轴两方面比对人工合成记录与井旁地震道，通过二者相关系数质控，最终得出时深转化关系。
在本实施例中，使用25HZ（目标层段地震数据主频）的雷克子波，根据补偿声波曲线、密度值曲线计算出的波阻抗信息得到合成地震记录。
通过观察地震波形数据确定鹰山组顶界标志层：与上覆巴楚组泥岩层相比，碳酸盐岩内幕区波形反射不规则，无一定方向，振幅可强可弱，同相轴可长可短连续性差；且具有非系统性同相轴反射终止和分叉现象。
在本实施例中，所述测井曲线数据与地震记录的时深转化关系为：
Figure 291202DEST_PATH_IMAGE173
其中，
Figure 807634DEST_PATH_IMAGE093
表示声波测井起始深度对应的地震信号时间，
Figure 851813DEST_PATH_IMAGE094
表示声波时差，
Figure 329062DEST_PATH_IMAGE096
表示 测井曲线数据采样间隔，T表示地震波双程旅行时。
单井井震标定图如图5所示；
根据钻录井、岩心数据确定目的层段古岩溶洞穴位置，依据充填物岩性与组合关系划分出沉积充填、垮塌充填、化学充填及混合充填的储层类型。运用测井曲线直方图、交会图的统计方式，优选对划分储层类型敏感的测井曲线参数。例如波阻抗曲线反映了岩石物理性质的差异，可用于区分古岩溶洞穴与围岩；自然伽马曲线利用了地质体放射性特征可以有效识别沉积砂泥岩储层充填物。
步骤S600，获取第一组PCA数据：从所述对储层敏感的特征参数中优选预设条数的标准化测井曲线数据，并通过PCA方法对所述对储层敏感的特征参数分析，获得对储层敏感的第一组PCA数据；
在本实施例中，所述PCA分析具体包括：
本实施例中优选波阻抗曲线、补偿中子曲线、深侧向电阻率曲线、浅侧向电阻率曲线、自然伽马曲线为5条对储层分类敏感的测井特征参数进行PCA分析。
步骤B100，计算样本均值：
Figure 179207DEST_PATH_IMAGE007
其中，
Figure 677709DEST_PATH_IMAGE008
表示曲线总样本数，i表示样本计数参量，
Figure 892789DEST_PATH_IMAGE009
表示第i个样本特征参数数 据，
Figure 981968DEST_PATH_IMAGE010
表示样本均值；
不同测井曲线反映了地质体不同物理学参数，由于各种测井曲线的探测方法有别，其量纲对数据分析的影响无法忽略。
步骤B200，对每条参数曲线进行归一化处理，获得归一化样本特征参数数据
Figure 245590DEST_PATH_IMAGE011
：
Figure 471035DEST_PATH_IMAGE012
其中，x表示样本特征参数采样点数据；
步骤B300，基于所述归一化样本特征参数数据
Figure 122596DEST_PATH_IMAGE011
，计算样本特征参数协方差
Figure 558126DEST_PATH_IMAGE013
：
Figure 750073DEST_PATH_IMAGE174
其中，
Figure 970970DEST_PATH_IMAGE015
表示第c个X样本特征参数数据，
Figure 527853DEST_PATH_IMAGE016
表示第c个Y样本特征参数数据；
步骤B400，基于所述样本特征参数协方差，构建协方差矩阵：
Figure 326045DEST_PATH_IMAGE017
Z表示样本特征参数数据；
包括p种样本特征参数数据的协方差矩阵为：
Figure 446316DEST_PATH_IMAGE018
其中，
Figure 849616DEST_PATH_IMAGE019
表示多个样本特征参数数据，
Figure 374138DEST_PATH_IMAGE020
表示第1个样本特征参数数据与第p个样本特征参数数据的协方差；K并没有具体的数值含义；
步骤B500，基于所述包括p种样本特征参数数据的协方差矩阵，获取该协方差矩阵 的特征值
Figure 68742DEST_PATH_IMAGE176
和p个特征向量
Figure 326548DEST_PATH_IMAGE177
，
Figure 412184DEST_PATH_IMAGE178
，
Figure 450547DEST_PATH_IMAGE179
；
基于所述协方差矩阵的特征值和特征向量计算包括p种特征值的PCA数据：
Figure 132195DEST_PATH_IMAGE025
其中，
Figure 510087DEST_PATH_IMAGE180
表示第j个PCA数据向量，
Figure 641991DEST_PATH_IMAGE027
表示第j个样本特征参数数据向量，
Figure 23775DEST_PATH_IMAGE028
表示第j 个样本特征参数数据的特征向量；j
Figure 633748DEST_PATH_IMAGE029
1,2,
Figure 866147DEST_PATH_IMAGE030
p,为计数参量；
步骤B600，基于包括p种特征值的PCA数据，将对应的特征值从大到小进行排序，得 到从大到小排序的特征值
Figure 309897DEST_PATH_IMAGE181
，…，
Figure 791694DEST_PATH_IMAGE182
计算方差累计贡献量：
Figure 408620DEST_PATH_IMAGE033
其中，
Figure 416897DEST_PATH_IMAGE034
表示由大到小排序后的第j个特征值，
Figure DEST_PATH_IMAGE183
为累计前k个协方差累计贡献量，
Figure 562707DEST_PATH_IMAGE184
；
当k从1开始不断增加时，选取
Figure 735062DEST_PATH_IMAGE183
首次大于预设的第三阈值时对应的前k个协方差 对应的PCA数据向量，即为本轮PCA分析获得的PCA数据。在本实施例中，所述第三阈值优选 的取85%。前k个特征参数对应的PCA数据反映了表征储层各项参数的
Figure 686838DEST_PATH_IMAGE183
部分信息，优选PCA 数据可用于进一步利用PCA数据对储层岩性、物性进行分析。
本实施例优选对5个测井样本特征参数通过PCA方法分析，得到5个PCA数据。PC1与PC2的方差贡献率达到了85%以上，PC1占48.72%，PC2占36.96%，可以代表储层的各类测井响应特征。
在本实施例中，经过第一轮PCA方法获得的PC1和PC2为：
Figure 80779DEST_PATH_IMAGE185
Figure 725387DEST_PATH_IMAGE186
其中，PC1表示第一次对所有数据点特征参数曲线进行PCA分析的第一组PCA数据 中的第一数据；PC2表示第一次对所有数据点特征参数曲线进行PCA分析的第一组PCA数据 中的第二数据；
Figure 385038DEST_PATH_IMAGE187
为所有数据点的自然伽马值；
Figure 281450DEST_PATH_IMAGE188
为所有数据点的深侧向电阻率值；
Figure 874106DEST_PATH_IMAGE189
为所有数据点的浅侧向电阻率值；
Figure 814249DEST_PATH_IMAGE190
为所有数据点的补偿中子值；
Figure 961196DEST_PATH_IMAGE191
为所有数 据点的波阻抗值；
步骤S700，构建等时格架模型：基于高精度三维地震振幅数据体所反映的沉积地层规律和所述时深转化关系构建等时格架模型。确定目标层段位置后，根据地震数据剖面图选取同相轴连续、沉积环境稳定的界面或者标志层作为待预测范围的顶底界面。基于实际地质构造背景和所述时深转化关系构建等时格架模型。
步骤S800，井间储层参数模拟：基于各样本井的第一组PCA数据，确定能够反映整体地质情况的最优样本数参量，选取最优样本数参量个地震波形相关性最高的样本井的第一组PCA数据构建初始模型，不断更正初始模型的参数，输出第一组高精度特征值模拟结果数据体，所述第一组高精度特征值模拟结果数据体为与所述第一组PCA数据一一对应的数据体；
在本实施例中，所述步骤S800，具体包括步骤S810-步骤S890：
步骤S810，步骤S810，任选一样本井作为参考目标井，设定初始的样本数参量为1；
步骤S820，根据波形相似性原则选取数量为样本数参量的样本井第一组PCA数据与参考目标井标第一组PCA数据进行相关性分析获得样本数参量-参考目标井的第一组PCA数据相关性值；
步骤S830，逐1增加样本数参量，重复步骤S820的方法获得各样本数参量对应的样本数参量-参考目标井的第一组PCA数据相关性值，将所有样本数参量-参考目标井的第一组PCA数据相关性值连接,获得参考井第一组PCA数据相关性随样本数参量变化的相关性曲线；
步骤S840，选取另一样本井作为参考目标井，重复步骤S810-步骤S830的方法，获得多个参考井第一组PCA数据相关性随样本数参量变化的相关性曲线，将所有参考井的第一组PCA数据波形相关性随样本数参量变化的相关性曲线拟合为整体相关性曲线，选取所述整体相关性曲线中相关性随样本数参量增加而升高最终保持平稳处的拐点，确定最优样本数参量；将所有参考井的第一组PCA数据波形相关性随样本数参量变化的相关性曲线拟合为整体相关性曲线的效果如图6所示；
本实施例中，在样本井中利用波形相似性和空间距离双变量优选低频结构相似的井作为空间估值样本，能够较好地反映空间结构的低频变化。
地震波形相似的两口井，表明所处沉积环境是相似的，其低频成分具有共性，可以增强反演结果低频段的确定性，同时约束了高频的取值范围，提高了反演与模拟结果的可靠性。
步骤S850，基于所述高精度的三维地震振幅数据体和等时格架模型，计算待测点位与样本井位的波形相关性，将所述波形相关性由大到小排序，选取最优样本数参量个地震波形相关性最高的样本井的第一组PCA数据；基于所述相关性最高的样本井的地震波形特征数据所对应的样本井，通过井间特征参数插值方式构建初始模型；
步骤S860，基于所述初始模型，选取最优样本数参量条地震波形关联度最高的样本井第一组PCA数据作为先验信息；
步骤S870，将所述初始模型与先验信息进行匹配滤波，进而获得最大似然函数；
步骤S880，基于所述最大似然函数和先验信息，在贝叶斯框架下求得后验概率统计分布密度，对所述后验概率统计分布密度进行采样获得目标函数；
在本实施例中，所述步骤S880，具体包括：
步骤S881，利用白噪声满足高斯分布的规律，将高精度特征值模拟结果数据体的参数表示为：
Figure 785933DEST_PATH_IMAGE098
其中，Y表示测井曲线高精度特征值模拟结果数据体的参数，X表示待求解的地下地层实际特征参数值，N表示随机噪声；
步骤S882，因为
Figure 374040DEST_PATH_IMAGE192
也满足高斯分布，可以确定初始目标函数为：
Figure DEST_PATH_IMAGE193
其中，
Figure 829292DEST_PATH_IMAGE101
表示与后验信息有关的函数，
Figure 653416DEST_PATH_IMAGE194
表示基于最优样本数选取样本井对样本 井的特征曲线进行匹配滤波后，求得后验概率统计分布密度，进而计算得到的特征参数期 望值，
Figure 750685DEST_PATH_IMAGE103
表示白噪声的协方差；
步骤S883，基于所述初始目标函数，通过最大后验估计，在目标函数中引入先验信息，获得稳定的目标函数为：
Figure 786774DEST_PATH_IMAGE104
其中，
Figure 553873DEST_PATH_IMAGE105
表示待模拟的特征参数，
Figure 737730DEST_PATH_IMAGE106
表示与地质、测井资料等先验信息有关的函数，
Figure 841952DEST_PATH_IMAGE107
表示用于协调
Figure 388340DEST_PATH_IMAGE108
和
Figure 982132DEST_PATH_IMAGE106
之间的相互影响的平滑参数。
步骤S890，以所述目标函数作为所述初始模型的输入，通过马尔科夫链蒙特卡罗方法MCMC和Metropolis-Hastings抽样准则对后验概率分布抽样，不断优化初始模型的参数，选取目标函数取最大值时的解作为随机实现，取多次随机实现的均值作为期望值输出，将所述期望值输出作为高精度特征值模拟结果数据体；所述高精度特征值模拟结果数据体中的参数对应于所述第一组PCA数据一一对应的特征参数。
步骤S900，刻画井间溶洞系统边界：基于样本井的第一组PCA数据，进行交会分析，获得第一组各类地质体离散数据点分布图；
基于所述岩性信息和物性信息获取第一组类地质体离散数据点分布图的解释结论，划分岩溶洞穴与围岩，进而构建第一组PCA数据交会版图，获取岩溶洞穴划分所需的第一组PCA数据门槛值；
在本实施例中，将PCA数据PC1与PC2进行交会分析，得到各类地质体离散数据点分布。基于所述成像测井、钻录井、岩心数据所得到的个别深度段的确定的岩性、物性信息得到优选PCA数据交会图解释结论，区分岩溶洞穴与围岩，构建优选PCA数据交会图版，获取岩溶洞穴划分所需的PC1、PC2数据门槛值。第一组PCA数据交会分析抓取岩溶洞穴数据点并框选门槛值示意图如图7所示；
基于所述第一组PCA数据门槛值，对所述第一组高精度特征值模拟结果数据体进行交会分析，获得岩溶洞穴三维空间形态特征。可以得到2m以上的古岩溶洞穴型储层结构特征与空间分布规律。三维空间形态与结构刻画剖面效果如图8所示。
步骤S1000，刻画洞穴内部充填岩性物性边界：所述岩溶洞穴数据点，将所述岩溶洞穴数据点所对应的对储层敏感的特征参数通过PCA方法进行分析，获得对充填物响应敏感的第二组PCA数据；
基于所述第二组PCA数据，通过如S700所述方法生成第二组高精度特征值模拟结果数据体；
基于所述对充填物响应敏感的第二组PCA数据，进行交会分析，获得第二组各类地质体离散数据点分布图；
基于所述岩性信息和物性信息获取第二组类地质体离散数据点分布图的解释结论，区分填充物类型，进而构建第二组PCA数据交会版图，获取岩溶洞内部填充物类型划分所需的第二组PCA数据门槛值；
基于所述第二组PCA数据门槛值，对所述第二组高精度特征值模拟结果数据体进行交会分析，获得岩溶洞穴内部不同填充物类型的三维空间形态特征；第二组PCA数据交会分析划分岩溶洞穴内部充填物并框选门槛值示意图如图9所示；
在本实施例中，提取PCA数据PC1与PC2交会图版得到的岩溶洞穴解释数据点，将其所对应的五条特征参数曲线再次进行PCA分析，获得对充填物响应敏感的第二组PCA数据F1、F2。将PCA数据F1与F2进行交会分析，得到各类地质体离散数据点分布。基于所述成像测井、钻录井、岩心数据所得到的个别深度段的确定的岩性、物性信息得到第二组PCA数据交会图解释结论，区分沉积充填与垮塌充填，构建第二组PCA数据交会图版，获取岩溶洞穴内部沉积充填与垮塌充填划分所需的F1、F2数据门槛值。
在本实施例中，经过第二轮PCA方法获得的F1和F2数据为：
Figure 997493DEST_PATH_IMAGE195
Figure 905406DEST_PATH_IMAGE196
其中，F1表示第二此对洞穴表征数据点特征参数曲线进行PCA分析的第二组PCA数 据中的第一PCA数据；F2为第二次对洞穴表征数据点特征参数曲线进行PCA分析的第二组 PCA数据中的第二PCA数据；
Figure 447246DEST_PATH_IMAGE197
为洞穴表征数据点的自然伽马值；
Figure 336573DEST_PATH_IMAGE198
为洞穴表征数据点 的深侧向电阻率值；
Figure 963864DEST_PATH_IMAGE199
为洞穴表征数据点的浅侧向电阻率值；
Figure 409889DEST_PATH_IMAGE200
为洞穴表征数据点的 补偿中子值；
Figure 681601DEST_PATH_IMAGE201
为洞穴表征数据点的波阻抗值
通过所述岩溶洞穴内部沉积充填与垮塌充填划分的解释结论门槛值对所述第二组PCA数据模拟结果数据体进行交会分析获得岩溶洞穴内部沉积充填与垮塌充填的三维空间形态特征；
步骤S1100，古岩溶洞穴结构与充填描述：基于所述岩溶洞穴三维空间形态特征和岩溶洞穴内部不同填充物类型的三维空间形态特征，采用岩性遮挡技术和三维雕刻技术雕刻出古岩溶洞穴空间展布和不同填充物类型的内部充填的发育特征。古岩溶洞穴内部充填物发育特征刻画剖面如图10所示；
本发明第二实施例的井震联合评价深层古岩溶储层充填特征的系统，所述系统包括：原始地球物理测井资料获取模块、地震数据获取模块、原始地球物理测井数据预处理模块、地震数据预处理模块、井震标定与特征参数选取模块、第一组PCA数据获取模块、等时格架模型构建模块、井间储层参数模拟模块、刻画井间溶洞系统边界模块、刻画洞穴内部充填岩性物性边界模块和古岩溶洞穴结构与充填描述模块；
所述原始地球物理测井资料获取模块，配置为通过测井设备获取各样本井的原始测井数据，包括：通过测量电极测量各样本井的自然电位SP，通过自然伽马井下装置和自然伽马地面仪器测量各样本井的自然伽马GR，通过井径臂获取各样本井的井径CAL；通过传统测井设备获取电阻率曲线数据：深侧向测井RLLD、浅侧向测井RLLS和微侧向测井RLLM，物性表征曲线数据：补偿中子CNL、补偿声波曲线AC和密度曲线DEN；
基于成像测井信息、钻井信息、录井信息和岩心信息获得个别深度段的确定岩性信息和物性信息，进而确定目的层位标志层的深度数据；
所述原始地球物理测井数据预处理模块，配置为基于所述样本井的原始测井数据，绘制测井曲线数据，并进行异常值处理和标准化处理，获得标准化测井曲线数据；
所述地震数据预处理模块，配置为基于所述地震波反射信号数据，通过混合相位子波反褶积和扩散滤波，获得高精度三维地震振幅数据体；
在本实施例中，基于井眼坐标与地震信号检波器坐标的相关关系，确定井眼与井旁地震道地震波形的一一对应关系，进而确定9条测井曲线的平面位置与其所对应的井旁地震道地震波形。
所述井震标定与特征参数选取模块，配置为井震标定与特征参数选取：基于所述标准化测井曲线数据中的补偿声波曲线和密度曲线DEN获取样本井的波阻抗曲线，进而计算反射系数曲线，获取雷克子波的优选频率使其与高精度三维地震振幅数据体的主频保持一致，将所述雷克子波与所述反射系数曲线进行褶积运算，得到合成地震记录，并将所述目的层位标志层的深度数据与所述目的层位标志层的等时三维展布相比对进行井震标定，计算所述合成地震记录与井旁地震道波形的相关度，当相关度大于或等于预设的第一阈值时，判定井震标定结果合格，获得测井曲线数据与地震记录的时深转化关系和对储层敏感的特征参数；
所述对储层敏感的特征参数，其获得方法包括：
通过对井旁不同地质体产生的测井参数的响应绘制柱状统计图，当某一标准化测井曲线数据的数值可以将不同测井解释结论预设的第二阈值以上的数据点区分开时，则选定该标准化测井曲线数据为对储层敏感的特征参数；所述对储层敏感的特征参数至少包括波阻抗IMP，还可包括井径CAL、自然伽马GR、自然电位SP，电阻率曲线数据：深侧向测井RLLD、浅侧向测井RLLS、微侧向测井RLLM，物性表征曲线数据：补偿中子CNL、补偿声波曲线AC、密度曲线DEN中的一个或多个；
所述第一组PCA数据获取模块，配置为从所述对储层敏感的特征参数中优选预设条数的标准化测井曲线数据，并通过PCA分析对所述对储层敏感的特征参数进行降维，获得第一组PCA数据；
所述等时格架模型构建模块，配置为基于高精度三维地震振幅数据体所反映的沉积地层规律和所述时深转化关系构建等时格架模型；
所述井间储层参数模拟模块，配置为基于各样本井的第一组PCA数据，确定能够反映整体地质情况的最优样本数参量，选取最优样本数参量个地震波形相关性最高的样本井的第一组PCA数据构建初始模型，不断更正初始模型的参数，输出第一组高精度特征值模拟结果数据体，所述第一组高精度特征值模拟结果数据体为与所述第一组PCA数据一一对应的数据体；
所述井间溶洞系统边界刻画模块，配置为基于样本井的第一组PCA数据，进行交会分析，获得第一组类地质体离散数据点分布图；
基于所述岩性信息和物性信息获取第一组各类地质体离散数据点分布图的解释结论，划分岩溶洞穴与围岩，进而构建第一组PCA数据交会版图，获取岩溶洞穴划分所需的第一组PCA数据门槛值；
基于所述第一组PCA数据门槛值，对所述第一组高精度特征值模拟结果数据体进行交会分析，获得岩溶洞穴三维空间形态特征；
所述洞穴内部充填岩性物性边界刻画模块，配置为刻画洞穴内部充填岩性物性边界：所述岩溶洞穴数据点，将所述岩溶洞穴数据点所对应的对储层敏感的特征参数通过PCA方法进行分析，获得对充填物响应敏感的第二组PCA数据；
基于所述第二组PCA数据，通过如井间储层参数模拟模块所述功能生成第二组高精度特征值模拟结果数据体；
基于所述对充填物响应敏感的第二组PCA数据，进行交会分析，获得第二组类地质体离散数据点分布图；
基于所述岩性信息和物性信息获取第二组类地质体离散数据点分布图的解释结论，区分填充物类型，进而构建第二组PCA数据交会版图，获取岩溶洞内部填充物类型划分所需的第二组PCA数据门槛值；
基于所述第二组PCA数据门槛值，对所述第二组高精度特征值模拟结果数据体进行交会分析，获得岩溶洞穴内部不同填充物类型的三维空间形态特征；
所述古岩溶洞穴结构与充填描述模块，配置为：基于所述岩溶洞穴三维空间形态特征和岩溶洞穴内部不同填充物类型的三维空间形态特征，采用岩性遮挡技术和三维雕刻技术雕刻出古岩溶洞穴空间展布和不同填充物类型的内部充填的发育特征。
所属技术领域的技术人员可以清楚地了解到，为描述的方便和简洁，上述描述的系统的具体工作过程及有关说明，可以参考前述方法实施例中的对应过程，在此不再赘述。
需要说明的是，上述实施例提供的古岩溶洞穴储层充填特征智能划分系统，仅以上述各功能模块的划分进行举例说明，在实际应用中，可以根据需要而将上述功能分配由不同的功能模块来完成，即将本发明实施例中的模块或者步骤再分解或者组合，例如，上述实施例的模块可以合并为一个模块，也可以进一步拆分成多个子模块，以完成以上描述的全部或者部分功能。对于本发明实施例中涉及的模块、步骤的名称，仅仅是为了区分各个模块或者步骤，不视为对本发明的不当限定。
本发明第三实施例的一种电子设备，包括：至少一个处理器；以及与至少一个所述处理器通信连接的存储器；其中，所述存储器存储有可被所述处理器执行的指令，所述指令用于被所述处理器执行以实现上述的井震联合评价深层古岩溶储层充填特征的方法。
本发明第四实施例的一种计算机可读存储介质，所述计算机可读存储介质存储有计算机指令，所述计算机指令用于被所述计算机执行以实现上述的井震联合评价深层古岩溶储层充填特征的方法。
所属技术领域的技术人员可以清楚地了解到，为描述的方便和简洁，上述描述的存储装置、处理装置的具体工作过程及有关说明，可以参考前述方法实施例中的对应过程，在此不再赘述。
可以以一种或多种程序设计语言或其组合来编写用于执行本申请的操作的计算机程序代码，上述程序设计语言包括面向对象的程序设计语言—诸如Java、Smalltalk、C++，还包括常规的过程式程序设计语言—诸如“C”语言或类似的程序设计语言。程序代码可以完全地在用户计算机上执行、部分地在用户计算机上执行、作为一个独立的软件包执行、部分在用户计算机上部分在远程计算机上执行、或者完全在远程计算机或服务器上执行。在涉及远程计算机的情形中，远程计算机可以通过任意种类的网络——包括局域网(LAN)或广域网(WAN)—连接到用户计算机，或者，可以连接到外部计算机(例如利用因特网服务提供商来通过因特网连接)。
附图中的流程图和框图，图示了按照本申请各种实施例的系统、方法和计算机程序产品的可能实现的体系架构、功能和操作。在这点上，流程图或框图中的每个方框可以代表一个模块、程序段、或代码的一部分，该模块、程序段、或代码的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。也应当注意，在有些作为替换的实现中，方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如，两个接连地表示的方框实际上可以基本并行地执行，它们有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，框图和/或流程图中的每个方框、以及框图和/或流程图中的方框的组合，可以用执行规定的功能或操作的专用的基于硬件的系统来实现，或者可以用专用硬件与计算机指令的组合来实现。
术语“第一”、“第二”等是用于区别类似的对象，而不是用于描述或表示特定的顺序或先后次序。
术语“包括”或者任何其它类似用语旨在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备/装置不仅包括那些要素，而且还包括没有明确列出的其它要素，或者还包括这些过程、方法、物品或者设备/装置所固有的要素。
至此，已经结合附图所示的优选实施方式描述了本发明的技术方案，但是，本领域技术人员容易理解的是，本发明的保护范围显然不局限于这些具体实施方式。在不偏离本发明的原理的前提下，本领域技术人员可以对相关技术特征做出等同的更改或替换，这些更改或替换之后的技术方案都将落入本发明的保护范围之内。

Claims (10)

1.一种井震联合评价深层古岩溶储层充填特征的方法，其特征在于，所述方法包括：
步骤S100，获取原始地球物理测井资料：通过测井设备获取各样本井的原始测井数据，包括：通过测量电极测量各样本井的自然电位SP，通过自然伽马井下装置和自然伽马地面仪器测量各样本井的自然伽马GR，通过井径臂获取各样本井的井径CAL；通过传统测井设备获取电阻率曲线数据：深侧向测井RLLD、浅侧向测井RLLS和微侧向测井RLLM，物性表征曲线数据：补偿中子CNL、补偿声波曲线AC和密度曲线DEN；；
基于成像测井信息、钻井信息、录井信息和岩心信息获得个别深度段的确定岩性信息和物性信息，进而确定目的层位标志层的深度数据；
步骤S200，获取地震数据，通过地震波激发装置与接收装置获取原始地震波反射信号数据，并根据原始地震波反射信号数据的波形获取目的层位标志层的等时三维展布；
步骤S300，原始地球物理测井数据预处理：基于所述样本井的所有原始测井数据，绘制测井曲线数据，并进行异常值处理和标准化处理，获得标准化测井曲线数据；
步骤S400，地震数据预处理：基于所述地震波反射信号数据，通过混合相位子波反褶积和扩散滤波，获得高精度三维地震振幅数据体；
步骤S500，井震标定与特征参数选取：基于所述标准化测井曲线数据中的补偿声波曲线和密度曲线DEN获取样本井的波阻抗曲线，进而计算反射系数曲线，获取雷克子波的优选频率使其与高精度三维地震振幅数据体的主频保持一致，将所述雷克子波与所述反射系数曲线进行褶积运算，得到合成地震记录，并将所述目的层位标志层的深度数据与所述目的层位标志层的等时三维展布相比对进行井震标定，计算所述合成地震记录与井旁地震道波形的相关度，当相关度大于或等于预设的第一阈值时，判定井震标定结果合格，获得测井曲线数据与地震记录的时深转化关系和对储层敏感的特征参数
所述对储层敏感的特征参数，其获得方法包括：
通过对井旁不同地质体产生的测井参数的响应绘制柱状统计图，当某一标准化测井曲线数据的数值可以将不同测井解释结论预设的第二阈值以上的数据点区分开时，则选定该标准化测井曲线数据为对储层敏感的特征参数；所述对储层敏感的特征参数至少包括波阻抗IMP，还可包括井径CAL、自然伽马GR、自然电位SP，电阻率曲线数据：深侧向测井RLLD、浅侧向测井RLLS、微侧向测井RLLM，物性表征曲线数据：补偿中子CNL、补偿声波曲线AC、密度曲线DEN中的一个或多个；
步骤S600，获取第一组PCA数据：从所述对储层敏感的特征参数中优选预设条数的标准化测井曲线数据，并通过PCA方法对所述对储层敏感的特征参数进行降维，获得对储层敏感的第一组PCA数据；
步骤S700，构建等时格架模型：基于高精度三维地震振幅数据体所反映的沉积地层规律和所述时深转化关系构建等时格架模型；
步骤S800，井间储层参数模拟：基于各样本井的第一组PCA数据，确定能够反映整体地质情况的最优样本数参量，选取最优样本数参量个地震波形相关性最高的样本井的第一组PCA数据构建初始模型，不断更正初始模型的参数，输出第一组高精度特征值模拟结果数据体，所述第一组高精度特征值模拟结果数据体为与所述第一组PCA数据一一对应的数据体；
步骤S900，刻画井间溶洞系统边界：基于样本井的第一组PCA数据，进行交会分析，获得第一组各类地质体离散数据点分布图；
基于所述岩性信息和物性信息获取第一组各类地质体离散数据点分布图的解释结论，划分岩溶洞穴与围岩，进而构建第一组PCA数据交会版图，获取岩溶洞穴划分所需的第一组PCA数据门槛值；
基于所述第一组PCA数据门槛值，对所述第一组高精度特征值模拟结果数据体进行交会分析，获得岩溶洞穴三维空间形态特征；
步骤S1000，刻画洞穴内部充填岩性物性边界：基于第一组PCA数据交会版图所述岩溶洞穴数据点，将所述岩溶洞穴数据点所对应的对储层敏感的特征参数通过PCA方法进行分析，获得对充填物响应敏感的第二组PCA数据；
基于所述第二组PCA数据，通过如S800所述方法生成第二组高精度特征值模拟结果数据体；
基于所述对充填物响应敏感的第二组PCA数据，进行交会分析，获得第二组各类地质体离散数据点分布图；
基于所述岩性信息和物性信息获取第二组各类地质体离散数据点分布图的解释结论，区分填充物类型，进而构建第二组PCA数据交会版图，获取古岩溶洞穴内部填充物类型划分所需的第二组PCA数据门槛值；
基于所述第二组PCA数据门槛值，对所述第二组高精度特征值模拟结果数据体进行交会分析，获得岩溶洞穴内部不同填充物类型的三维空间形态特征；
步骤S1100，古岩溶洞穴结构与充填描述：基于所述岩溶洞穴三维空间形态特征和岩溶洞穴内部不同填充物类型的三维空间形态特征，采用岩性遮挡技术和三维雕刻技术雕刻出古岩溶洞穴空间展布和不同填充物类型的内部充填的发育特征。
2.根据权利要求1所述的井震联合评价深层古岩溶储层充填特征的方法，其特征在于，所述通过测量电极测量各样本井的自然电位SP，通过自然伽马井下装置和自然伽马地面仪器测量各样本井的自然伽马GR，具体包括：
所述通过测量电极测量各样本井的自然电位SP：
将测量电极N设置于地面，测量电机M通过电缆设置于井下；
沿井轴提升测量电极M测量自然电位随井深的变化；
自然电位值的计算方法为：
Figure 605133DEST_PATH_IMAGE001
其中，
Figure 404462DEST_PATH_IMAGE002
为总自然电位，
Figure 193426DEST_PATH_IMAGE003
为扩散电位系数，
Figure 888850DEST_PATH_IMAGE004
为扩散吸附电位系数，
Figure 661634DEST_PATH_IMAGE005
为泥浆滤液 电阻率值，
Figure 467916DEST_PATH_IMAGE006
为地层水电阻率值；
所述通过自然伽马井下装置和自然伽马地面仪器测量各样本井的自然伽马GR：
自然伽马井下装置包括探测器、放大器、高压电源；
通过探测器获取自然伽马射线，并将所述自然伽马射线转化为电脉冲信号，并通过放大器进行放大；
所述地面仪器把每分钟形成的电脉冲计数转化为电位差进行记录。
3.根据权利要求1所述的井震联合评价深层古岩溶储层充填特征的方法，其特征在于，所述PCA分析具体包括：
步骤B100，计算样本均值：
Figure 376966DEST_PATH_IMAGE007
其中，
Figure 240361DEST_PATH_IMAGE008
表示曲线总样本数，i表示样本计数参量，
Figure 234862DEST_PATH_IMAGE009
表示第i个样本特征参数数据，
Figure 844835DEST_PATH_IMAGE010
表 示样本均值；
步骤B200，对每条参数曲线进行归一化处理，获得归一化样本特征参数数据
Figure 608391DEST_PATH_IMAGE011
：
Figure 645617DEST_PATH_IMAGE012
其中，x表示样本特征参数采样点数据；
步骤B300，基于所述归一化样本特征参数数据
Figure 127414DEST_PATH_IMAGE011
，计算样本特征参数协方差
Figure 275499DEST_PATH_IMAGE013
：
Figure 893562DEST_PATH_IMAGE014
其中，
Figure 367269DEST_PATH_IMAGE015
表示第c个X样本特征参数数据，
Figure 70782DEST_PATH_IMAGE016
表示第c个Y样本特征参数数据；
步骤B400，基于所述样本特征参数协方差，构建样本特征参数协方差矩阵：
Figure 288137DEST_PATH_IMAGE017
Z表示样本特征参数数据；
包括p种样本特征参数数据的协方差矩阵为：
Figure 760707DEST_PATH_IMAGE018
其中，
Figure 405315DEST_PATH_IMAGE019
表示多个样本特征参数数据，
Figure 596125DEST_PATH_IMAGE020
表示第1个样本特征参数数据与第p个 样本特征参数数据的协方差；
步骤B500，基于所述包括p种样本特征参数数据的协方差矩阵，获取该协方差矩阵的特 征值
Figure 351591DEST_PATH_IMAGE021
和p个特征向量
Figure 209826DEST_PATH_IMAGE022
，
Figure 494177DEST_PATH_IMAGE023
，
Figure 437862DEST_PATH_IMAGE024
；
基于所述协方差矩阵的特征值和特征向量计算包括p种特征值的PCA数据：
Figure 731440DEST_PATH_IMAGE025
其中，
Figure 444181DEST_PATH_IMAGE026
表示第j个PCA数据向量，
Figure 165012DEST_PATH_IMAGE027
表示第j个样本特征参数数据向量，
Figure 333344DEST_PATH_IMAGE028
表示第j个样 本特征参数数据的特征向量；j
Figure 165034DEST_PATH_IMAGE029
1,2,
Figure 997861DEST_PATH_IMAGE030
p,为计数参量；
步骤B600，基于p种特征值的PCA数据，将对应的特征值从大到小进行排序，得到从大到 小排序的特征值
Figure 889593DEST_PATH_IMAGE031
，…，
Figure 542291DEST_PATH_IMAGE032
计算方差累计贡献量：
Figure 177672DEST_PATH_IMAGE033
其中，
Figure 599426DEST_PATH_IMAGE034
表示从大到小排序后的第j个特征值，
Figure 927639DEST_PATH_IMAGE035
为累计前k个协方差累计贡献量，
Figure 802054DEST_PATH_IMAGE036
；
当k从1开始不断增加时，选取
Figure 241126DEST_PATH_IMAGE035
首次大于预设的第三阈值时对应的前k个协方差对应 的PCA数据向量，即为本轮PCA分析获得的PCA数据。
4.根据权利要求1所述的井震联合评价深层古岩溶储层充填特征的方法，其特征在于，所述步骤S400，具体包括：
步骤S410，基于所述地震波反射信号数据，将频率域地震记录褶积模型表示为：
Figure 517387DEST_PATH_IMAGE037
其中，
Figure 16501DEST_PATH_IMAGE038
表示经傅氏变换后的地震记录，
Figure 378212DEST_PATH_IMAGE039
表示经傅氏变换后的子波，
Figure 620975DEST_PATH_IMAGE040
表示 经傅氏变换后的反射系数的频谱，
Figure 751742DEST_PATH_IMAGE041
表示角频率；
步骤S420，将所述频率域地震记录褶积模型的等式两边取对数转化为线性系统，获得线性地震记录褶积模型：
Figure 421758DEST_PATH_IMAGE042
步骤S430，将所述线性地震记录褶积模型进行反傅氏变换，获得复赛谱序列：
Figure 536344DEST_PATH_IMAGE043
其中，
Figure 317218DEST_PATH_IMAGE044
表示地震波形记录的复赛谱序列，
Figure 568071DEST_PATH_IMAGE045
表示地震子波的复赛谱序列，
Figure 143409DEST_PATH_IMAGE046
表 示地层反射系数的复赛谱序列，
Figure 745292DEST_PATH_IMAGE047
表示地震波形记录时间；
步骤S440，基于所述复赛谱序列，通过低通滤波器进行子波与反射系数分离，提取子波振幅谱；
步骤S450，通过最小二乘法获取模拟地震子波振幅谱：
Figure 64277DEST_PATH_IMAGE048
其中，最小二乘法拟合参数
Figure 432286DEST_PATH_IMAGE049
为常数，
Figure 178525DEST_PATH_IMAGE050
表示子波振幅谱，
Figure 267704DEST_PATH_IMAGE051
和
Figure 124802DEST_PATH_IMAGE052
表示f的多 项式，f表示地震子波的频率；
步骤S460，基于所述模拟地震子波振幅谱，获得子波最大相位分量和最小相位分量；
设子波
Figure 350247DEST_PATH_IMAGE053
的最大相位分量为
Figure 532966DEST_PATH_IMAGE054
、最小相位分量为
Figure 843862DEST_PATH_IMAGE055
，则子波
Figure 504650DEST_PATH_IMAGE053
为：
Figure 850181DEST_PATH_IMAGE056
振幅谱的复赛谱中表示为：
Figure 938223DEST_PATH_IMAGE057
其中，振幅谱的复赛谱
Figure 736414DEST_PATH_IMAGE058
在复赛谱的正、负轴上对称显示，
Figure 200894DEST_PATH_IMAGE059
为地震子波最大相 位分量
Figure 400931DEST_PATH_IMAGE060
所对应的最小相位函数的复赛谱，
Figure 659874DEST_PATH_IMAGE061
为地震子波最小相位分量
Figure 945362DEST_PATH_IMAGE062
所对应 的最大相位函数的复赛谱；
步骤S470，基于所述振幅谱中的复赛谱确定一组具有相同振幅谱的混合相位子波集合，不断调整俞氏子波参数，保持低频、拓展高频和适当提高主频构建期望输出子波形态，在井曲线控制下以信噪比谱作参考寻找提高分辨率与保真度之间的最佳平衡点，获得整形后波形数据；
步骤S480，基于所述整形后波形数据，构建张量扩散模型：
Figure 947953DEST_PATH_IMAGE063
其中，
Figure 2497DEST_PATH_IMAGE064
表示扩散时间，
Figure 697920DEST_PATH_IMAGE065
表示散度算子，D表示扩散滤波器的张量型扩散系数，U表示 扩散滤波结果，
Figure 205125DEST_PATH_IMAGE066
表示
Figure 11407DEST_PATH_IMAGE064
=0时的扩散滤波结果，
Figure 920457DEST_PATH_IMAGE067
表示
Figure 789712DEST_PATH_IMAGE068
时刻的整形后波形数据，作 为张量扩散模型的初始条件，
Figure 784212DEST_PATH_IMAGE069
表示扩散滤波结果的梯度；
基于所述张量扩散模型构建梯度结构张量：
Figure 394185DEST_PATH_IMAGE070
其中，U表示扩散滤波结果，
Figure 157742DEST_PATH_IMAGE071
表示梯度向量张量积；
Figure 460547DEST_PATH_IMAGE072
表示尺度为
Figure 676765DEST_PATH_IMAGE073
的高斯函数：
Figure 90429DEST_PATH_IMAGE074
其中，r表示计算半径；
结构张量的特征向量为：
Figure 708492DEST_PATH_IMAGE075
其中，
Figure 182199DEST_PATH_IMAGE076
、
Figure 885712DEST_PATH_IMAGE077
和
Figure 837488DEST_PATH_IMAGE078
表示为梯度结构张量的3个特征向量，可视为局部正交坐标系，
Figure 575637DEST_PATH_IMAGE076
指向 地震信号的梯度方向，
Figure 220245DEST_PATH_IMAGE077
和
Figure 411055DEST_PATH_IMAGE078
组成的平面平行于地震信号的局部结构特征，
Figure 900942DEST_PATH_IMAGE079
、
Figure 493597DEST_PATH_IMAGE080
和
Figure 309106DEST_PATH_IMAGE081
为 分别与
Figure 987212DEST_PATH_IMAGE076
、
Figure 546370DEST_PATH_IMAGE077
和
Figure 993532DEST_PATH_IMAGE078
对应的三个特征值；
步骤S490，基于所述结构张量的特征向量分别计算线状结构置信度量、面状结构置信度量和扩散张量；
所述现状结构置信度量
Figure 448784DEST_PATH_IMAGE082
为：
Figure 879765DEST_PATH_IMAGE083
所述面状结构置信度量
Figure 708525DEST_PATH_IMAGE084
为：
Figure 541352DEST_PATH_IMAGE085
所述扩散张量D为：
Figure 167505DEST_PATH_IMAGE086
其中，
Figure 85783DEST_PATH_IMAGE087
、
Figure 986743DEST_PATH_IMAGE088
和
Figure 142918DEST_PATH_IMAGE089
表示扩散张量的三个非负特征值，他们分别表示扩散滤波器沿
Figure 471131DEST_PATH_IMAGE076
、
Figure 345546DEST_PATH_IMAGE077
和
Figure 50197DEST_PATH_IMAGE078
这三个特征方向的滤波强度；
步骤S4100，重复步骤S480-S490的步骤，直至达到预设的迭代次数，获得扩散滤波结果，即为所述高精度的三维地震振幅数据体。
5.根据权利要求4所述的井震联合评价深层古岩溶储层充填特征的方法，其特征在于，所述模拟地震子波振幅谱的计算方法为：
定位地震波反射信号数据中振幅谱的极大值和极大值所对应的频率；
将所述地震信号振幅谱的极大值和所述模拟地震子波振幅谱通过最小二乘法拟合的 方式得到参数
Figure 326457DEST_PATH_IMAGE090
和
Figure 825572DEST_PATH_IMAGE091
多项式的系数，进而获得拟合的极大值相应频率振幅值；
将所述地震信号振幅谱的极大值除以所拟合出的相应频率的振幅值，进而用商拟合多 项式
Figure 187283DEST_PATH_IMAGE052
的系数。
6.根据权利要求1所述的井震联合评价深层古岩溶储层充填特征的方法，其特征在于，所述步骤S800，具体包括：
步骤S810，任选一样本井作为参考目标井，设定初始的样本数参量为1；
步骤S820，根据波形相似性原则选取数量为样本数参量的样本井第一组PCA数据与参考目标井标第一组PCA数据进行相关性分析获得样本数参量-参考目标井的第一组PCA数据相关性值；
步骤S830，逐1增加样本数参量，重复步骤S720的方法获得各样本数参量对应的样本数参量-参考目标井的第一组PCA数据相关性值，将所有样本数参量-参考目标井的第一组PCA数据相关性值连接,获得参考井第一组PCA数据相关性随样本数参量变化的相关性曲线；
步骤S840，选取另一样本井作为参考目标井，重复步骤S810-步骤S830的方法，获得多个参考井第一组PCA数据相关性随样本数参量变化的相关性曲线，将所有参考井的第一组PCA数据波形相关性随样本数参量变化的相关性曲线拟合为整体相关性曲线，选取所述整体相关性曲线中相关性随样本数参量增加而升高最终保持平稳处的拐点，确定最优样本数参量；
步骤S850，基于所述高精度的三维地震振幅数据体和等时格架模型，计算待测点位与样本井位的波形相关性，将所述波形相关性由大到小排序；基于所述相关性最高的样本井的地震波形特征数据所对应的样本井，通过井间特征参数插值方式构建初始模型；
步骤S860，基于所述初始模型，选取最优样本数参量条地震波形关联度最高的样本井第一组PCA数据作为先验信息；
步骤S870，将所述初始模型与先验信息进行匹配滤波，进而获得最大似然函数；
步骤S880，基于所述最大似然函数和先验信息，在贝叶斯框架下求得后验概率统计分布密度，对所述后验概率统计分布密度进行采样获得目标函数；
步骤S890，以所述目标函数作为所述初始模型的输入，通过马尔科夫链蒙特卡罗方法MCMC和Metropolis-Hastings抽样准则对后验概率分布抽样，不断优化初始模型的参数，选取目标函数取最大值时的解作为随机实现，取多次随机实现的均值作为期望值输出，将所述期望值输出作为高精度特征值模拟结果数据体；所述高精度特征值模拟结果数据体中的参数对应于所述第一组PCA数据一一对应的特征参数。
7.根据权利要求6所述的井震联合评价深层古岩溶储层充填特征的方法，其特征在于，所述步骤S880，具体包括：
步骤S881，利用白噪声满足高斯分布的规律，将高精度特征值模拟结果数据体的参数表示为：
Figure 164466DEST_PATH_IMAGE092
其中，Y表示测井曲线高精度特征值模拟结果数据体的参数，X表示待求解的地下地层实际特征参数值，N表示随机噪声；
步骤S882，因为
Figure 560812DEST_PATH_IMAGE093
也满足高斯分布，可以确定初始目标函数为：
Figure 965249DEST_PATH_IMAGE094
其中，
Figure 79835DEST_PATH_IMAGE095
表示与后验信息有关的函数，
Figure 595130DEST_PATH_IMAGE096
表示基于最优样本数选取样本井对样本井的 特征曲线进行匹配滤波后，求得后验概率统计分布密度，进而计算得到的特征参数期望值，
Figure 111562DEST_PATH_IMAGE097
表示白噪声的协方差；
步骤S883，基于所述初始目标函数，通过最大后验估计，在目标函数中引入先验信息，获得稳定的目标函数为：
Figure 686900DEST_PATH_IMAGE098
其中，
Figure 288783DEST_PATH_IMAGE099
表示待模拟的特征参数，
Figure 607769DEST_PATH_IMAGE100
表示与地质、测井资料等先验信息有关的函数，
Figure 713128DEST_PATH_IMAGE101
表 示用于协调
Figure 727876DEST_PATH_IMAGE102
和
Figure 285896DEST_PATH_IMAGE100
之间的相互影响的平滑参数。
8.根据权利要求6所述的井震联合评价深层古岩溶储层充填特征的方法，其特征在于，所述步骤S890，具体步骤为：
步骤S891，设M为目标空间，n为总样本数，m为马尔科夫链趋于平稳时的样本数；
步骤S892，预设一条马尔科夫链，使马尔科夫链收敛至平稳分布；
步骤S893，由M中的某一点
Figure 674152DEST_PATH_IMAGE103
出发，通过马尔科夫链进行抽样模拟，产生点序列：
Figure 634018DEST_PATH_IMAGE104
；
步骤S894，函数
Figure 816738DEST_PATH_IMAGE105
的期望估计为：
Figure 127633DEST_PATH_IMAGE106
其中，n表示产生的总样本数，m表示马尔科夫链达到平稳时的样本数，k表示累加参量；
步骤S895，选取一转移函数
Figure 788422DEST_PATH_IMAGE107
和初始值
Figure 868373DEST_PATH_IMAGE108
，若第i次迭代开始时的参数值为
Figure 221994DEST_PATH_IMAGE109
，则第i次迭代过程为：
从
Figure 20186DEST_PATH_IMAGE107
中抽取一个备选值
Figure 484665DEST_PATH_IMAGE110
，计算备选值
Figure 419123DEST_PATH_IMAGE110
的接受概率
Figure 943646DEST_PATH_IMAGE111
：
Figure 229133DEST_PATH_IMAGE112
步骤S896，以
Figure 231725DEST_PATH_IMAGE113
，置
Figure 286268DEST_PATH_IMAGE114
，以概率
Figure 981692DEST_PATH_IMAGE115
，置
Figure 488896DEST_PATH_IMAGE116
；
步骤S897，不断扰动所述初始模型的参数，重复步骤C792-步骤C796的方法达到预设的 迭代次数n，获得后验样本
Figure 295178DEST_PATH_IMAGE117
，进而计算后延分布的各阶矩阵获得期望输出 值，将期望值输出作为第一组高精度特征值模拟结果数据体。
9.根据权利要求1所述的井震联合评价深层古岩溶储层充填特征的方法，其特征在于，所述步骤S300，包括：
步骤S310，基于所述原始测井数据绘制原始测井曲线数据；
步骤S320，基于所述原始测井曲线数据，去除离群点获得去除离群点的测井曲线数据；
步骤S330，基于所述去除离群点的测井曲线数据，叠合工区内所有样本井位的单个测井曲线柱状图数据，通过整合阈值，获得标准化测井曲线数据。
10.一种井震联合评价深层古岩溶储层充填特征的系统，其特征在于，所述系统包括：原始地球物理测井资料获取模块、地震数据获取模块、原始地球物理测井数据预处理模块、地震数据预处理模块、井震标定与特征参数选取模块、第一组PCA数据获取模块、等时格架模型构建模块、井间储层参数模拟模块、刻画井间溶洞系统边界模块、刻画洞穴内部充填岩性物性边界模块和古岩溶洞穴结构与充填描述模块；
所述原始地球物理测井资料获取模块，配置为通过测井设备获取各样本井的原始测井数据，包括：通过测量电极测量各样本井的自然电位SP，通过自然伽马井下装置和自然伽马地面仪器测量各样本井的自然伽马GR，通过井径臂获取各样本井的井径CAL；通过传统测井设备获取电阻率曲线数据：深侧向测井RLLD、浅侧向测井RLLS和微侧向测井RLLM，物性表征曲线数据：补偿中子CNL、补偿声波曲线AC和密度曲线DEN；
基于成像测井信息、钻井信息、录井信息和岩心信息获得个别深度段的确定岩性信息和物性信息，进而确定目的层位标志层的深度数据；
所述地震数据获取模块，配置为通过地震波激发装置与接收装置获取原始地震波反射信号数据，并根据原始地震波反射信号数据的波形获取目的层位标志层的等时三维展布；
所述原始地球物理测井数据预处理模块，配置为基于所述样本井的原始测井数据，绘制测井曲线数据，并进行异常值处理和标准化处理，获得标准化测井曲线数据；
所述地震数据预处理模块，配置为基于所述地震波反射信号数据，通过混合相位子波反褶积和扩散滤波，获得高精度三维地震振幅数据体；
所述井震标定与特征参数选取模块，配置为井震标定与特征参数选取：基于所述标准化测井曲线数据中的补偿声波曲线和密度曲线DEN获取样本井的波阻抗曲线，进而计算反射系数曲线，获取雷克子波的优选频率使其与高精度三维地震振幅数据体的主频保持一致，将所述雷克子波与所述反射系数曲线进行褶积运算，得到合成地震记录，并将所述目的层位标志层的深度数据与所述目的层位标志层的等时三维展布相比对进行井震标定，计算所述合成地震记录与井旁地震道波形的相关度，当相关度大于或等于预设的第一阈值时，判定井震标定结果合格，获得测井曲线数据与地震记录的时深转化关系和对储层敏感的特征参数；
所述对储层敏感的特征参数，其获得方法包括：
通过对井旁不同地质体产生的测井参数的响应绘制柱状统计图，当某一标准化测井曲线数据的数值可以将不同测井解释结论预设的第二阈值以上的数据点区分开时，则选定该标准化测井曲线数据为对储层敏感的特征参数；所述对储层敏感的特征参数至少包括波阻抗IMP，还可包括井径CAL、自然伽马GR、自然电位SP，电阻率曲线数据：深侧向测井RLLD、浅侧向测井RLLS、微侧向测井RLLM，物性表征曲线数据：补偿中子CNL、补偿声波曲线AC、密度曲线DEN中的一个或多个；
所述第一组PCA数据获取模块，配置为从所述对储层敏感的特征参数中优选预设条数的标准化测井曲线数据，并通过PCA分析对所述对储层敏感的特征参数进行降维，获得第一组PCA数据；
所述等时格架模型构建模块，配置为基于高精度三维地震振幅数据体所反映的沉积地层规律和所述时深转化关系构建等时格架模型；
所述井间储层参数模拟模块，配置为基于各样本井的第一组PCA数据，确定能够反映整体地质情况的最优样本数参量，选取最优样本数参量个地震波形相关性最高的样本井的第一组PCA数据构建初始模型，不断更正初始模型的参数，输出第一组高精度特征值模拟结果数据体，所述第一组高精度特征值模拟结果数据体为与所述第一组PCA数据一一对应的数据体；
所述井间溶洞系统边界刻画模块，配置为基于样本井的第一组PCA数据，进行交会分析，获得第一组类地质体离散数据点分布图；
基于所述岩性信息和物性信息获取第一组各类地质体离散数据点分布图的解释结论，划分岩溶洞穴与围岩，进而构建第一组PCA数据交会版图，获取岩溶洞穴划分所需的第一组组PCA数据门槛值；
基于所述第一组PCA数据门槛值，对所述第一组高精度特征值模拟结果数据体进行交会分析，获得岩溶洞穴三维空间形态特征；
所述洞穴内部充填岩性物性边界刻画模块，配置为刻画洞穴内部充填岩性物性边界：所述岩溶洞穴数据点，将所述岩溶洞穴数据点所对应的对储层敏感的特征参数通过PCA方法进行分析，获得对充填物响应敏感的第二组PCA数据；
基于所述第二组PCA数据，通过如井间储层参数模拟模块所述功能生成第二组高精度特征值模拟结果数据体；
基于所述对充填物响应敏感的第二组PCA数据，进行交会分析，获得第二组类地质体离散数据点分布图；
基于所述岩性信息和物性信息获取第二组类地质体离散数据点分布图的解释结论，区分填充物类型，进而构建第二组PCA数据交会版图，获取岩溶洞内部填充物类型划分所需的第二组PCA数据门槛值；
基于所述第二组PCA数据门槛值，对所述第二组高精度特征值模拟结果数据体进行交会分析，获得岩溶洞穴内部不同填充物类型的三维空间形态特征；
所述古岩溶洞穴结构与充填描述模块，配置为：基于所述岩溶洞穴三维空间形态特征和岩溶洞穴内部不同填充物类型的三维空间形态特征，采用岩性遮挡技术和三维雕刻技术雕刻出古岩溶洞穴空间展布和不同填充物类型的内部充填的发育特征。
CN202111066355.7A 2021-09-13 2021-09-13 井震联合评价深层古岩溶储层充填特征的方法与系统 Active CN113759425B (zh)

Priority Applications (2)

Application Number Priority Date Filing Date Title
CN202111066355.7A CN113759425B (zh) 2021-09-13 2021-09-13 井震联合评价深层古岩溶储层充填特征的方法与系统
US17/697,562 US11500117B1 (en) 2021-09-13 2022-03-17 Method and system for evaluating filling characteristics of deep paleokarst reservoir through well-to-seismic integration

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111066355.7A CN113759425B (zh) 2021-09-13 2021-09-13 井震联合评价深层古岩溶储层充填特征的方法与系统

Publications (2)

Publication Number Publication Date
CN113759425A true CN113759425A (zh) 2021-12-07
CN113759425B CN113759425B (zh) 2022-04-01

Family

ID=78795030

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111066355.7A Active CN113759425B (zh) 2021-09-13 2021-09-13 井震联合评价深层古岩溶储层充填特征的方法与系统

Country Status (2)

Country Link
US (1) US11500117B1 (zh)
CN (1) CN113759425B (zh)

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114994759A (zh) * 2022-08-02 2022-09-02 中国科学院地质与地球物理研究所 基于gan网络的碳封存箱智能识别方法和系统
CN114994758A (zh) * 2022-08-02 2022-09-02 北京京鲁聚源能源科技有限公司 碳酸盐岩断控储层的波阻抗提取与结构表征方法和系统
CN115357994A (zh) * 2022-10-20 2022-11-18 中国地质大学（北京） 软岩隧道围岩参数空间随机场建模方法、装置及设备
CN116933196A (zh) * 2023-09-19 2023-10-24 中国科学院地质与地球物理研究所 多维测井数据异常值智能剔除的方法与系统
CN116976705A (zh) * 2023-09-19 2023-10-31 中国科学院地质与地球物理研究所 深地油气精准导航砂泥岩地层物性评价方法与系统
CN117150929A (zh) * 2023-10-19 2023-12-01 中国科学院地质与地球物理研究所 深地油气精准导航井震结合更新碎屑岩岩性模型系统

Families Citing this family (14)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115828072B (zh) * 2022-12-07 2023-11-17 中国科学技术大学 不可通约测井曲线缺失段重构方法及存储介质
CN115932968B (zh) * 2023-01-09 2023-07-18 西南石油大学 一种基于地震振幅比值属性的碳酸盐岩薄储层预测方法
CN116256801B (zh) * 2023-05-16 2023-07-21 中国科学院地质与地球物理研究所 基于图像融合的深地油气精准导航断层表征方法与系统
CN116340731B (zh) * 2023-05-29 2023-08-22 北京科技大学 一种岩石声发射信号的降噪办法及系统
CN116402957B (zh) * 2023-06-09 2023-08-04 武汉中旺亿能科技发展有限公司 基于全直径岩心ct扫描的储层构型相控智能建模方法
CN116840916B (zh) * 2023-07-04 2024-03-26 成都理工大学 一种地震速度信号和加速度信号联合子波提取方法
CN116592951B (zh) * 2023-07-17 2023-09-08 陕西西特电缆有限公司 一种电缆数据智能采集方法及系统
CN116680545B (zh) * 2023-08-02 2023-10-20 西安核音智言科技有限公司 一种基于马尔科夫随机场的煤矿出井人员预测方法
CN116755414B (zh) * 2023-08-22 2023-11-07 山东新巨龙能源有限责任公司 一种基于物联网的矿采设备监管系统
CN116881636B (zh) * 2023-09-08 2023-12-08 山东远盾网络技术股份有限公司 基于人工智能的频域信号分析检测方法
CN116957363B (zh) * 2023-09-19 2023-12-01 中国科学院地质与地球物理研究所 深地油气精准导航砂泥岩地层构造评价方法与系统
CN117075212B (zh) * 2023-10-16 2024-01-26 吉林大学 一种隧道磁共振裂隙结构成像方法
CN117665933B (zh) * 2023-11-27 2024-06-11 西南石油大学 一种低分辨率地震资料限制下稀井网区的沉积微相刻画方法
CN117784244B (zh) * 2024-02-28 2024-05-10 中国石油大学(华东) 基于变纵波速度的细粒混积岩孔隙压力预测方法及系统

Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20060055403A1 (en) * 2004-04-30 2006-03-16 Schlumberger Technology Corporation Method for determining characteristics of earth formations
US20090125239A1 (en) * 2004-03-18 2009-05-14 Baker Hughes Incorporated Rock and Fluid Properties Prediction From Downhole Measurements Using Linear and Nonlinear Regression
CN103529475A (zh) * 2013-04-19 2014-01-22 中国石油大学(华东) 一种识别和解释碳酸盐岩古岩溶储层三维结构的方法
FR3053125A1 (fr) * 2016-06-23 2017-12-29 Storengy Procede de caracterisation du sous-sol d'une region utilisant des signaux sismiques passifs, et systeme correspondant
JP6494147B1 (ja) * 2018-07-24 2019-04-03 株式会社地震科学探査機構 地体区分システム、及び地体区分方法
CN109711429A (zh) * 2018-11-22 2019-05-03 中国石油天然气股份有限公司 一种储层评价分类方法及装置
CN110674841A (zh) * 2019-08-22 2020-01-10 中国石油天然气集团有限公司 一种基于聚类算法的测井曲线识别方法
US20200033495A1 (en) * 2018-07-30 2020-01-30 Cgg Services Sas Method and apparatus for obtaining statistical measures of geological properties values related to lateral wells using seismic-derived maps
CN112213797A (zh) * 2020-10-28 2021-01-12 中国石油大学(华东) 一种利用测井相联合叠后地震属性预测储层产能的方法
CN113009559A (zh) * 2021-03-02 2021-06-22 中国科学院地质与地球物理研究所 基于多类型地球物理数据的地震评价方法

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
RU2760102C2 (ru) * 2016-09-07 2021-11-22 Чайна Петролеум Энд Кемикал Корпорейшн Способ и система автоматического распознавания центра залежи в карстовой пещере
US11143776B2 (en) * 2020-02-17 2021-10-12 China Petroleum & Chemical Corporation Computer-implemented method and system for small cave recognition using seismic reflection data
US20210264262A1 (en) * 2020-02-21 2021-08-26 Saudi Arabian Oil Company Physics-constrained deep learning joint inversion

Patent Citations (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20090125239A1 (en) * 2004-03-18 2009-05-14 Baker Hughes Incorporated Rock and Fluid Properties Prediction From Downhole Measurements Using Linear and Nonlinear Regression
US20060055403A1 (en) * 2004-04-30 2006-03-16 Schlumberger Technology Corporation Method for determining characteristics of earth formations
CN103529475A (zh) * 2013-04-19 2014-01-22 中国石油大学(华东) 一种识别和解释碳酸盐岩古岩溶储层三维结构的方法
FR3053125A1 (fr) * 2016-06-23 2017-12-29 Storengy Procede de caracterisation du sous-sol d'une region utilisant des signaux sismiques passifs, et systeme correspondant
JP6494147B1 (ja) * 2018-07-24 2019-04-03 株式会社地震科学探査機構 地体区分システム、及び地体区分方法
US20200033495A1 (en) * 2018-07-30 2020-01-30 Cgg Services Sas Method and apparatus for obtaining statistical measures of geological properties values related to lateral wells using seismic-derived maps
CN109711429A (zh) * 2018-11-22 2019-05-03 中国石油天然气股份有限公司 一种储层评价分类方法及装置
CN110674841A (zh) * 2019-08-22 2020-01-10 中国石油天然气集团有限公司 一种基于聚类算法的测井曲线识别方法
CN112213797A (zh) * 2020-10-28 2021-01-12 中国石油大学(华东) 一种利用测井相联合叠后地震属性预测储层产能的方法
CN113009559A (zh) * 2021-03-02 2021-06-22 中国科学院地质与地球物理研究所 基于多类型地球物理数据的地震评价方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
WENHAO ZHENG，等: "Electrofacies classification of deeply buried carbonate strata using machine learning methods: A case study on ordovician paleokarst reservoirs in Tarim Basin", 《MARINE AND PETROLEUM GEOLOGY》 *
伍家和，等: "缝洞型碳酸盐岩油藏溶洞描述技术研究", 《石油地质与工程》 *
夏竹，等: "井震联合薄储层沉积微相表征实例研究", 《石油地球物理勘探》 *

Cited By (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114994759A (zh) * 2022-08-02 2022-09-02 中国科学院地质与地球物理研究所 基于gan网络的碳封存箱智能识别方法和系统
CN114994758A (zh) * 2022-08-02 2022-09-02 北京京鲁聚源能源科技有限公司 碳酸盐岩断控储层的波阻抗提取与结构表征方法和系统
CN114994758B (zh) * 2022-08-02 2022-10-28 北京京鲁聚源能源科技有限公司 碳酸盐岩断控储层的波阻抗提取与结构表征方法和系统
CN115357994A (zh) * 2022-10-20 2022-11-18 中国地质大学（北京） 软岩隧道围岩参数空间随机场建模方法、装置及设备
CN116933196A (zh) * 2023-09-19 2023-10-24 中国科学院地质与地球物理研究所 多维测井数据异常值智能剔除的方法与系统
CN116976705A (zh) * 2023-09-19 2023-10-31 中国科学院地质与地球物理研究所 深地油气精准导航砂泥岩地层物性评价方法与系统
CN116976705B (zh) * 2023-09-19 2023-12-22 中国科学院地质与地球物理研究所 深地油气精准导航砂泥岩地层物性评价方法与系统
CN116933196B (zh) * 2023-09-19 2023-12-26 中国科学院地质与地球物理研究所 多维测井数据异常值智能剔除的方法与系统
CN117150929A (zh) * 2023-10-19 2023-12-01 中国科学院地质与地球物理研究所 深地油气精准导航井震结合更新碎屑岩岩性模型系统
CN117150929B (zh) * 2023-10-19 2023-12-26 中国科学院地质与地球物理研究所 深地油气精准导航井震结合更新碎屑岩岩性模型系统

Also Published As

Publication number Publication date
US11500117B1 (en) 2022-11-15
CN113759425B (zh) 2022-04-01

Similar Documents

Publication Publication Date Title
CN113759425B (zh) 井震联合评价深层古岩溶储层充填特征的方法与系统
CN113759424B (zh) 基于频谱分解和机器学习的岩溶储层充填分析方法和系统
CN108931814B (zh) 一种基于多属性融合的基岩裂缝预测的方法
US5831935A (en) Method for geophysical processing and interpretation using seismic trace difference for analysis and display
US7133779B2 (en) Automated borehole geology and petrophysics interpretation using image logs
WO2020123084A1 (en) Machine learning-augmented geophysical inversion
CN114966856B (zh) 基于多频带地震资料的碳封存场址优选方法、系统和设备
CN107783187B (zh) 一种将测井速度和地震速度结合建立三维速度场的方法
CN114994758B (zh) 碳酸盐岩断控储层的波阻抗提取与结构表征方法和系统
CN111506861B (zh) 一种目的层有利区域裂缝强度计算方法
Kianoush et al. Compressional and shear interval velocity modeling to determine formation pressures in an oilfield of SW Iran
WO2022159698A1 (en) Method and system for image-based reservoir property estimation using machine learning
CN109283577B (zh) 一种地震层位标定方法
CN115857047B (zh) 一种地震储层综合预测方法
CN111077578B (zh) 岩层分布预测方法和装置
CN107679614B (zh) 一种基于粒子群优化的声波时差实时提取方法
CN113514884A (zh) 一种致密砂岩储层预测方法
CN113806674A (zh) 古河道纵向尺度的量化方法、装置、电子设备及存储介质
CN114462703A (zh) 声波参数曲线预测、测井曲线预测方法及电子设备
CN113848593A (zh) 一种定量预测含煤地层中岩浆岩侵蚀区的方法
CN112147676A (zh) 一种煤层及夹矸厚度预测方法
CN113589365A (zh) 基于时频域信息的储层尖灭线描述方法
CN114609666B (zh) 一种页岩薄储层的预测方法、装置、设备及存储介质
CN108931813B (zh) 一种识别礁盖、礁核、礁基的方法
Obilo et al. Use of seismic inversion attributes in field development planning

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