具体实施方式

下面结合附图和实施例对本申请作进一步的详细说明。可以理解的是，此处所描述的具体实施例仅用于解释相关发明，而非对该发明的限定。另外还需要说明的是，为了便于描述，附图中仅示出了与有关发明相关的部分。

需要说明的是，在不冲突的情况下，本申请中的实施例及实施例中的特征可以相互组合。下面将参考附图并结合实施例来详细说明本申请。

为了更清晰地对本发明基于频谱分解和机器学习的岩溶储层充填分析方法进行说明，下面结合图1对本发明实施例中各步骤展开详述。

本发明第一实施例的基于频谱分解和机器学习的岩溶储层充填分析方法，包括步骤S100-步骤S1000，各步骤详细描述如下：

经实际地震数据统计分析表明，古岩溶洞穴所产生的地震波形特征通常由多组波谷和波峰组成，经研究发现，这种特征性反射是由地震波干涉造成的。并且反射波的振幅与洞穴充填物组合有关。Xu等（2016）通过物理模型模拟同样可以发现不同形状、厚度、组合关系的古岩溶洞穴型目标体会产生多样的地震波形反射特征，即反射形态特征对于洞穴直径、洞穴形态与分布规律有关。学者们进一步发现在等时地层格架的同一相带内，波形特征的相似性可以代表具有一定相关性的岩性组合。因此，基于波形指示的思路进行等时界面分区反演与特征参数模拟，利用了地震波形的横向变化信息，更好地体现了沉积环境的约束，更符合沉积地质规律。

步骤S100，获取原始地球物理测井资料：通过测井设备获取各样本井的原始测井数据，包括：通过测量电极测量各样本井的自然电位SP，通过自然伽马井下装置和自然伽马地面仪器测量各样本井的自然伽马GR，通过井径臂获取各样本井的井径CAL；通过传统测井设备获取电阻率曲线数据：深侧向测井RLLD、浅侧向测井RLLS和微侧向测井RLLM，物性表征曲线数据：补偿中子CNL、补偿声波曲线AC和密度曲线DEN；

基于成像测井信息、钻井信息、录井信息和岩心信息获得个别深度段的确定岩性信息和物性信息，进而确定目的层位标志层的深度数据；在本实施例中，在工区井眼处探测深度范围为5500~5750m的九条常规测井曲线数据，采样间隔设置为0.01m；运用三维地震勘探方法，借助地震波激发源与地震信号检波器，得到工区面积约27km²的三维地震资料数据，信号记录双程旅行时为4s，采样点时间间隔为1ms，探测深度超过6000m。

在本实施例中，所述通过测量电极测量各样本井的自然电位SP，通过自然伽马井下装置和自然伽马地面仪器测量各样本井的自然伽马GR，具体包括：

所述通过测量电极测量各样本井的自然电位SP：

将测量电极N设置于地面，测量电机M通过电缆设置于井下；

沿井轴提升测量电极M测量自然电位随井深的变化；

自然电位值的计算方法为：

其中，为总自然电位，为扩散电位系数，为扩散吸附电位系数，为泥 浆滤液电阻率值，为地层水电阻率值；

所述通过自然伽马井下装置和自然伽马地面仪器测量各样本井的自然伽马GR：

自然伽马井下装置包括探测器、放大器、高压电源；

通过探测器获取自然伽马射线，并将所述自然伽马射线转化为电脉冲信号，并通过放大器进行放大；

所述地面仪器把每分钟形成的电脉冲计数转化为电位差进行记录。

双侧向电极系由主电极、监督电极、环状屏蔽电极与柱状电极组成。主电极居 中，上下对称分布监督电极、和、，以及环状屏蔽电极、，在、的外侧对 称位置上加了两个柱状电极。深侧向电极系中两个柱状电极是屏蔽电极、；浅侧向电 极系中两个柱状电极是回路电极和。在电极系较远处装有对比电极N和深侧向电极系 的回路电极B；微球形聚焦电极系的主电极是矩形片状电极，依次向外的矩形框电极是测 量电极、辅助电极、监督电极和，回路电极B装在仪器外壳上或极板支撑架上。测 量任一监督电极与对比电极之间的电位差变化即反映介质电阻率的变化，其视电阻率表达 式为：

式中，为电极系系数，深侧向测井、浅侧向测井与微球型聚焦测井各不相同； 为监督电极M的电位；为主电流。

井下仪器的发射换能器晶体振动，引起周围介质的质点发生振动，产生向井内泥 浆及岩层中传播的声波。在井中可以用接收换能器R、R2先后接收到滑行波，记录时差，进 而测量地层的声波速度。到达和的时刻分别为和，那么到达两个接收换能器的时 间差为:

式中，为泥浆声速；为地层的声速

密度测井仪器包括有一个伽马源，两个接收伽马射线的探测器，即长源距探测器和短源距探测器。它们安装在滑板上，测井时被推靠到井壁上。在下井仪器的上方装有辅助电子线路。通常用137Cs作伽马源，它发射的伽马射线具有中等能量(0.661MeV)，用它照射物质只能产生康普顿散射和光电效应。地层的密度不同，则对伽马光子的散射和吸收的能力不同，探测器接收到的伽马光子的计数率也就不同。已知通过距离为L的伽马光子的计数率N为：

若只存在康普顿散射，则即为康普顿散射吸收系数，变换得：

式中，C为与地层有关的常数，L为接收源距离伽马源的距离。

源距选定后，对仪器进行刻度标定，找到和N的这种关系，则记录散射伽马光子 计数率N就可以测得地层的密度

补偿中子测井由一个中子发射源与两个接收探测器组成。通过对两个探测器计 数，得到两个计数率和，取其比值可以反映地层的氢含量：

式中，与为两个探测器到中子源的距离；为减速长度，反映了地层含氢量。

井径测量时，将仪器下人井预计的深度上，然后通过传统的方式打开井径臂，于是互成90°的四个井径腿便在弹簀力的作用下向外伸张，其末端紧贴井壁。随着仪器的向上提升，井径臂就会由于井径的变化而发生张缩，并带动连杆作上下运动。连杆与一个电位器的滑动端相连，于是井径的变化便可转换成电阻的变化。当给该可动电阻上通以一-定强度的电流时，可动电阻的某一固定端与滑动端之间的电位差将随着其间电阻值的变化而变化。于是测量这一电位差，便可反映井径大小：

式中，为井径测量仪器电位差；为初始值；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，基于所述地震波反射信号数据，将频率域地震记录褶积模型表示为：

其中，表示经傅氏变换后的地震记录，表示经傅氏变换后的子波， 表示经傅氏变换后的反射系数的频谱，表示角频率；

步骤S420，将所述频率域地震记录褶积模型的等式两边取对数转化为线性系统，获得线性地震记录褶积模型：

步骤S430，将所述线性地震记录褶积模型进行反傅氏变换，获得复赛谱序列：

其中，表示地震波形记录的复赛谱序列，表示地震子波的复赛谱序列，表示地层反射系数的复赛谱序列，表示地震波形记录时间；

步骤S440，基于所述复赛谱序列，通过低通滤波器进行子波与反射系数分离，提取子波振幅谱；子波与反射系数序列平滑程度的差异在复赛谱中易于区分：子波能量出现于原点附近，而反射系数序列则远离原点。利用低通滤波器就可以将复赛谱中的子波与反射系数分离达到子波振幅谱提取的目的。

步骤S450，通过最小二乘法获取模拟地震子波振幅谱：

其中，最小二乘法拟合参数为常数，表示子波振幅谱，和表示f 的多项式，f表示地震子波的频率；

在本实施例中，所述模拟地震子波振幅谱的计算方法为：

定位地震波反射信号数据中振幅谱的极大值和极大值所对应的频率；

将所述地震信号振幅谱的极大值和所述模拟地震子波振幅谱通过最小二乘法拟 合的方式得到参数和多项式的系数，进而获得拟合的极大值相应频率振幅值；

将所述地震信号振幅谱的极大值除以所拟合出的相应频率的振幅值，进而用商拟 合多项式的系数；

步骤S460，基于所述模拟地震子波振幅谱，获得子波最大相位分量和最小相位分量；

设子波的最大相位分量为、最小相位分量为，则子波为：

振幅谱的复赛谱中表示为：

其中，振幅谱的复赛谱在复赛谱的正、负轴上对称显示，为地震子波最 大相位分量所对应的最小相位函数的复赛谱，为地震子波最小相位分量所 对应的最大相位函数的复赛谱；

步骤S470，基于所述振幅谱中的复赛谱确定一组具有相同振幅谱的混合相位子波集合，不断调整俞氏子波参数，保持低频、拓展高频和适当提高主频构建期望输出子波形态，在井曲线控制下以信噪比谱作参考寻找提高分辨率与保真度之间的最佳平衡点，获得整形后波形数据；子波振幅谱如图4所示；

经过混合相位子波反褶积后，拓展了地震数据的有效频带，高频部分得到了合理加强。在地震波形上表现为同相轴数目增多，更容易反映地震波反射信息的细节变化，并且在振幅、相位、频率方面改善了同一反射波组波形的一致性。在古岩溶洞穴地震响应上，“串珠状”反射特征尤为明显，且串珠内部形态的细节能够清晰地显示，代表了不同结构特征与充填物组合的复杂古岩溶洞穴型储层地震反射，有助于后期精细地质解释。

Fhemers和Hocker于2003年首次将扩散滤波技术应用在地震资料的处理解释中。该技术不仅可以有效压制噪声，而且可以尽可能保留地震数据中的细节：如地质体边缘、断层、不整合面、尖灭等，为后续的地震解释和储层预测工作提供了可靠的地震资料，大大提高油气勘探开发的成功率。

为达到衰减地震噪声，并增强地质结构特征的扩散滤波效果，寻求扩散张量是该方法最关键的步骤：

步骤S480，基于所述整形后波形数据，构建张量扩散模型：

其中，表示扩散时间，表示散度算子，D表示扩散滤波器的张量型扩散系数，U 表示扩散滤波结果，表示=0时的扩散滤波结果，表示时刻的整形后波形数 据，作为张量扩散模型的初始条件，表示扩散滤波结果的梯度；

基于所述张量扩散模型构建梯度结构张量：

其中，U表示扩散滤波结果，表示梯度向量张量积；

表示尺度为的高斯函数为：

其中，r表示计算半径；

结构张量的特征向量为：

其中，、和表示为梯度结构张量的3个特征向量，可视为局部正交坐标系， 指向地震信号的梯度方向，和组成的平面平行于地震信号的局部结构特征，、和为分别与、和对应的三个特征值；

步骤S490，基于所述结构张量的特征向量分别计算线状结构置信度量、面状结构置信度量和扩散张量；

所述现状结构置信度量为：

所述面状结构置信度量为：

所述扩散张量D为：

其中，、和表示扩散张量的三个非负特征值，他们分别表示扩散滤波器沿 、和这三个特征方向的滤波强度；

步骤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（目标层段地震数据主频）的雷克子波，根据补偿声波曲线、密度值曲线计算出的波阻抗信息得到合成地震记录。

通过观察地震波形数据确定鹰山组顶界标志层：与上覆巴楚组泥岩层相比，碳酸盐岩内幕区波形反射不规则，无一定方向，振幅可强可弱，同相轴可长可短连续性差；且具有非系统性同相轴反射终止和分叉现象。

在本实施例中，所述测井曲线数据与地震记录的时深转化关系为：

其中，表示声波测井起始深度对应的地震信号时间，表示声波时差，表示 测井曲线数据采样间隔，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，利用白噪声满足高斯分布的规律，将高精度特征值模拟结果数据体的参数表示为：

其中，Y表示测井曲线高精度特征值模拟结果数据体的参数，X表示待求解的地下地层实际特征参数值，N表示随机噪声；

步骤S782，因为也满足高斯分布，可以确定初始目标函数为：

其中，表示与后验信息有关的函数，表示基于最优样本数选取样本井对样本 井的特征曲线进行匹配滤波后，求得后验概率统计分布密度，进而计算得到的特征参数期 望值，表示白噪声的协方差；

步骤S783，基于所述初始目标函数，通过最大后验估计，在目标函数中引入先验信息，获得稳定的目标函数为：

其中，表示波阻抗或待模拟的特征参数，表示与地质、测井资料等先验信息有 关的函数，表示用于协调和之间的相互影响的平滑参数。

步骤S790，以所述目标函数作为所述初始模型的输入，通过马尔科夫链蒙特卡罗方法（Markov chain Monte Carlo method，MCMC）和Metropolis-Hastings抽样准则对后验概率分布抽样，不断更正初始模型的参数，选取目标函数取最大值时的解作为随机实现，取多次随机实现的均值作为期望值输出，将所述期望值输出作为高精度特征值模拟结果数据体；所述高精度特征值模拟结果数据体中的参数与所述对储层敏感的特征参数一一对应。

步骤S791，设M为目标空间，n为总样本数，m为马尔科夫链趋于平稳时的样本数；

步骤S792，预设一条马尔科夫链，使马尔科夫链收敛至平稳分布；

步骤S793，由M中的某一点出发，通过马尔科夫链进行抽样模拟，产生点序列：；

步骤S794，函数的期望估计为：

其中，n表示产生的总样本数，m表示马尔科夫链达到平稳时的样本数，k表示累加参量；

步骤S795，选取一转移函数和初始值，若第i次迭代开始时的参数 值为，则第i次迭代过程为：

从中抽取一个备选值，计算备选值的接受概率：

步骤S796，以，置，以概率，置；

步骤S797，不断扰动所述初始模型的参数，重复步骤S792-步骤S796的方法达到预 设的迭代次数n，获得后验样本，进而计算后验分布的各阶矩阵获得期望 输出值，将期望值输出作为高精度特征值模拟结果数据体。

步骤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)—连接到用户计算机，或者，可以连接到外部计算机(例如利用因特网服务提供商来通过因特网连接)。

附图中的流程图和框图，图示了按照本申请各种实施例的系统、方法和计算机程序产品的可能实现的体系架构、功能和操作。在这点上，流程图或框图中的每个方框可以代表一个模块、程序段、或代码的一部分，该模块、程序段、或代码的一部分包含一个或多个用于实现规定的逻辑功能的可执行指令。也应当注意，在有些作为替换的实现中，方框中所标注的功能也可以以不同于附图中所标注的顺序发生。例如，两个接连地表示的方框实际上可以基本并行地执行，它们有时也可以按相反的顺序执行，这依所涉及的功能而定。也要注意的是，框图和/或流程图中的每个方框、以及框图和/或流程图中的方框的组合，可以用执行规定的功能或操作的专用的基于硬件的系统来实现，或者可以用专用硬件与计算机指令的组合来实现。

术语“第一”、“第二”等是用于区别类似的对象，而不是用于描述或表示特定的顺序或先后次序。

术语“包括”或者任何其它类似用语旨在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者设备/装置不仅包括那些要素，而且还包括没有明确列出的其它要素，或者还包括这些过程、方法、物品或者设备/装置所固有的要素。

至此，已经结合附图所示的优选实施方式描述了本发明的技术方案，但是，本领域技术人员容易理解的是，本发明的保护范围显然不局限于这些具体实施方式。在不偏离本发明的原理的前提下，本领域技术人员可以对相关技术特征做出等同的更改或替换，这些更改或替换之后的技术方案都将落入本发明的保护范围之内。