沉积地层旋回划分方法及装置

摘要

本发明提供一种沉积地层旋回划分方法及装置，其中方法包括：对获取的反映地层原始沉积环境及岩性变化的资料进行预处理，以获得预处理资料；对预处理资料进行经验模态分解，以获取多个本征模态函数；对于每个本征模态函数，对本征模态函数进行变换，以获取与本征模态函数对应的瞬时频率曲线；根据瞬时频率曲线，获取米兰科维奇旋回曲线；根据米兰科维奇旋回曲线，获取可容纳空间变化曲线；根据可容纳空间变化曲线，获取可容纳空间上升半旋回和可容纳空间下降半旋回。上述划分方法计算过程客观，受人为因素影响较小，是一种可操作性较强的沉积地层旋回划分方法。

说明书

技术领域

本发明涉及地质勘探领域，尤其涉及一种沉积地层旋回划分方法及装置。

背景技术

高频旋回识别是地质勘探工作重要组成内容，进行高精度的层序地层划分与对比是进行精细储层预测必要基础。

现有技术将经验模态分解应用于地层旋回识别的方法有：在建立对应目的层的地震或测井模型基础上，通过频率变化对地层旋回进行识别与划分(CN101349764A；CN101349764B；孙兴刚等,希尔伯特-黄变换时频分析在沉积旋回划分中的应用,油气地质与采收率,2012年11月)；直接将不同级别的本征模态函数与不同级次的层序级别相对应(徐敬领等，基于经验模态分解法的层序地层划分及对比研究,石油物探,2010年3月；李秋思,海拉尔盆地贝14区储层描述,中国海洋大学,2013年5月)；通过原曲线及各本征模态函数曲线间相关系数值进行优选，选择与原曲线相关性较好的本征模态函数进行加和形成复合曲线，在复合曲线基础上进行旋回划分(赵妮等,经验模态分解法在测井层序界面识别中的应用,测井技术,2014年8月)；对原始信号及各本征模态函数进行傅里叶变换求取频率谱进行比较，帮助判断不同级次沉积旋回对应的本征模态函数(张翔等,基于测井资料的经验模态分解法的沉积旋回界面划分,石油天然气学报,2010年10月)；目的层深度段的发育时间除以该段各本征模态函数的波数来求取平均周期，周期与地质周期对比判断各本征模态函数对应的地质周期(郭志永,北大西洋U1313站深海沉积物记录的早更新世气候变化,南京大学,2012年5月)。现有技术在识别旋回界面方面：通常是认为某级次的本征模态函数不同频率段间的突变区域代表地层边界的位置(徐敬领等，基于经验模态分解法的层序地层划分及对比研究,石油物探,2010年3月；文政等,基于测井数据Hilbert-Huang变换的地层层序划分,测井技术2007年12月；李秋思,海拉尔盆地贝14区储层描述,中国海洋大学,2013年5月；张翔等,基于测井资料的经验模态分解法的沉积旋回界面划分,石油天然气学报,2010年10月)。

上述方法存在问题有：层序界面的识别不具有与地层成因相联系确切的物理定义，或界面的识别需要人工经验性判识，受人为因素影响，缺乏客观性；无法应用深度上非连续的数据，如地球化学分析数据，进行数学分析。

发明内容

本发明提供一种沉积地层旋回划分方法及装置，用于解决现有技术中对沉积地层进行划分时需要人工经验性判识，受人为因素影响，缺乏客观性的技术问题。

本发明一方面提供一种沉积地层旋回划分方法，包括：

对获取的反映地层原始沉积环境及岩性变化的资料进行预处理，以获得预处理资料，资料包括岩石样品地球化学分析数据、地球物理测井数据和地震数据；

对预处理资料进行经验模态分解，以获取多个本征模态函数；

对于每个本征模态函数，对本征模态函数进行变换，以获取与本征模态函数对应的瞬时频率曲线；

根据瞬时频率曲线，获取米兰科维奇旋回曲线；

根据米兰科维奇旋回曲线，获取可容纳空间变化曲线；

根据可容纳空间变化曲线，获取可容纳空间上升半旋回和可容纳空间下降半旋回。

进一步的，对获取的反映地层原始沉积环境及岩性变化的资料进行预处理，以获得预处理资料，具体包括：

对地球物理资料进行归一化与标准化处理，以获得标准化资料，其中，地球物理资料包括地球物理测井数据和地震数据；

对目的层段包含的标准化资料进行曲线拼接，使得曲线拼接处理后获取的预处理资料连续且能够覆盖完整目的层段深度范围；

对岩石样品地球化学分析数据进行二次样条差值处理，使得经过二次样条差值处理后获取的预处理资料连续且能够覆盖完整目的层段深度范围；

预处理资料包括对标准化资料进行曲线拼接后所获得的数据和对岩石样品地球化学分析数据进行二次样条差值处理后的数据。

进一步的，根据瞬时频率曲线，获取米兰科维奇旋回曲线，具体包括：

获取相邻三条瞬时频率曲线，将在目的层段深度范围内超过第一预设阈值的瞬时频率曲线上的点的比值最接近第二预设阈值的三条曲线作为米兰科维奇旋回的曲线；

或者获取瞬时频率曲线中目的层段深度范围内的主频，将主频的比值最接近第二预设阈值的瞬时频率曲线作为米兰科维奇旋回曲线；

或者获取本征模态函数在目的层段深度范围内的周期数量，并根据目的层段深度范围内发育的地质历史时间的时限除以周期数量，以获得瞬时频率曲线的地质历史时间周期，将地质历史时间周期最接近第三预设阈值的三条瞬时频率曲线作为米兰科维奇旋回曲线。

进一步的，根据米兰科维奇旋回曲线，获取可容纳空间变化曲线，具体包括：

从米兰科维奇旋回曲线中任意选择一条曲线，用费希尔图解原理进行分析，以获取可容纳空间变化曲线。

进一步的，根据可容纳空间变化曲线，获取可容纳空间上升半旋回和可容纳空间下降半旋回，具体包括：

对可容纳空间变化曲线求取导数，以获得所述可容纳空间变化曲线的极值点；

根据可容纳空间变化曲线的极值点将目的层段划分成至少一个包含可容纳空间上升半旋回和可容纳空间下降半旋回的部分。

本发明另一方面提供一种沉积地层旋回划分装置，包括：

预处理模块，用于对获取的反映地层原始沉积环境及岩性变化的资料进行预处理，以获得预处理资料，资料包括岩石样品地球化学分析数据、地球物理测井数据和地震数据；

本征模态函数获取模块，用于对预处理资料进行经验模态分解，以获取多个本征模态函数；

瞬时频率曲线获取模块，用于对于每个本征模态函数，对本征模态函数进行变换，以获取与本征模态函数对应的瞬时频率曲线；

米兰科维奇旋回曲线获取模块，用于根据瞬时频率曲线，获取米兰科维奇旋回曲线；

可容纳空间变化曲线获取模块，用于根据米兰科维奇旋回曲线，获取可容纳空间变化曲线；

划分模块，用于根据可容纳空间变化曲线，获取可容纳空间上升半旋回和可容纳空间下降半旋回。

进一步的，预处理模块具体包括：

第一预处理子模块，用于对地球物理资料进行归一化与标准化处理，以获得标准化资料，其中，地球物理资料包括地球物理测井数据和地震数据；

曲线拼接子模块，用于对目的层段包含的标准化资料进行曲线拼接，使得曲线拼接处理后获取的预处理资料连续且能够覆盖完整目的层段深度范围；

二次样条差值处理子模块，用于对岩石样品地球化学分析数据进行二次样条差值处理，使得经过二次样条差值处理后获取的预处理资料连续且能够覆盖完整目的层段深度范围；

预处理资料包括对标准化资料进行曲线拼接后所获得的数据和对岩石样品地球化学分析数据进行二次样条差值处理后的数据。

进一步的，米兰科维奇旋回曲线获取模块具体用于：

获取相邻三条瞬时频率曲线，将在目的层段深度范围内超过第一预设阈值的瞬时频率曲线上的点的比值最接近第二预设阈值的三条曲线作为米兰科维奇旋回的曲线；

或者获取瞬时频率曲线中目的层段深度范围内的主频，将主频的比值最接近第二预设阈值的瞬时频率曲线作为米兰科维奇旋回曲线；

或者获取本征模态函数在目的层段深度范围内的周期数量，并根据目的层段深度范围内发育的地质历史时间的时限除以周期数量，以获得瞬时频率曲线的地质历史时间周期，将地质历史时间周期最接近第三预设阈值的三条瞬时频率曲线作为米兰科维奇旋回曲线。

进一步的，可容纳空间变化曲线获取模块具体用于：

从米兰科维奇旋回曲线中任意选择一条曲线，用费希尔图解原理进行分析，以获取可容纳空间变化曲线。

进一步的，划分模块具体包括：

求导子模块，用于对可容纳空间变化曲线求取导数，以获得可容纳空间变化曲线的极值点；

划分子模块，用于根据可容纳空间变化曲线的极值点将目的层段划分成至少一个包含可容纳空间上升半旋回和可容纳空间下降半旋回的部分。

本发明提供的沉积地层旋回划分方法及装置，对获取的反映地层原始沉积环境及岩性变化的资料进行预处理，然后对预处理后的资料进行经验模态分解，获得本征模态函数，之后根据本征模态函数获取米兰科维奇旋回曲线，进一步获取可容纳空间变化曲线，最终根据可容纳空间变化曲线，获取可容纳空间上升半旋回和可容纳空间下降半旋回，从而实现对沉积地层的划分，上述划分方法可以利用地球化学分析数据等在深度上不连续的资料进行计算，而不仅仅局限于应用地球物理测井数据和地震数据，计算过程客观，受人为因素影响较小，是一种可操作性较强的地层旋回划分方法，为识别层序地层界面或沉积旋回划分的界面提供一个量化的判识依据。

附图说明

在下文中将基于实施例并参考附图来对本发明进行更详细的描述。其中：

图1为本发明实施例一提供的沉积地层旋回划分方法的流程示意图；

图2为本发明实施例二提供的沉积地层旋回划分方法的流程示意图；

图3为本发明实施例三提供的沉积地层旋回划分装置的结构示意图；

图4为本发明实施例四提供的沉积地层旋回划分装置的结构示意图。

在附图中，相同的部件使用相同的附图标记。附图并未按照实际的比例绘制。

具体实施方式

下面将结合附图对本发明作进一步说明。

实施例一

图1为本发明实施例一提供的沉积地层旋回划分方法的流程示意图；如图1所示，本实施例提供一种沉积地层旋回划分方法，包括步骤101-步骤106。

其中，步骤101，对获取的反映地层原始沉积环境及岩性变化的资料进行预处理，以获得预处理资料，资料包括岩石样品地球化学分析数据、地球物理测井数据和地震数据。

本实施例中选取岩石样品地球化学分析数据、地球物理测井数据和地震数据作为资料。根据以上数据与沉积相、岩相变化的相关性，或对已发表文献记载的能够反映地层原始沉积环境及岩性变化的资料进行选择。其中，地球物理测井数据与地震数据一般为连续资料，岩石样品地球化学分析数据属于不连续资料。

步骤102，对预处理资料进行经验模态分解，以获取多个本征模态函数。

对步骤101中预处理资料进行经验模态分解，根据经验模态分解可以将上述预处理资料分解为有限数量个本征模态函数。

步骤103，对于每个本征模态函数，对本征模态函数进行变换，以获取与本征模态函数对应的瞬时频率曲线。对每个本征模态函数进行希尔伯特变换，每个本征模态函数均会求得一条瞬时频率曲线。

步骤104，根据瞬时频率曲线，获取米兰科维奇旋回曲线。

具体的，获取米兰科维奇旋回曲线包括以下三种方法：获取相邻三条瞬时频率曲线，将在目的层段深度范围内超过第一预设阈值的瞬时频率曲线上的点的比值最接近第二预设阈值的三条曲线作为米兰科维奇旋回的曲线；或者获取瞬时频率曲线中目的层段深度范围内的主频，将主频的比值最接近第二预设阈值的瞬时频率曲线作为米兰科维奇旋回曲线；或者获取本征模态函数在目的层段深度范围内的周期数量，并根据目的层段深度范围发育的地质历史时间的时限除以周期数量，以获得各瞬时频率曲线的地质历史时间周期，将地质历史时间周期最接近第三预设阈值的三条瞬时频率曲线作为米兰科维奇旋回曲线。

第一预设阈值、第二预设阈值和第三预设阈值可根据实际情况进行设置，在本实施例中，第一预设阈值设置为50％，第二预设阈值设置为4:2.5:1，第二预设阈值包括三个值，分别为25ka、41ka及100-400ka。上述第一种获取米兰科维奇旋回曲线的方法中，由于在步骤102中对预处理资料进行经验模态分解时，获取的多个本征模态函数会按照高频到低频进行排序，所以根据本征模态函数获取到的相应的瞬时频率曲线也会有顺序，因此，通过对相邻三条瞬时频率曲线进行分析，求取相邻三条瞬时频率曲线的比值，将瞬时频率曲线上超过50％的点的比值最接近4:2.5:1的三条曲线作为米兰科维奇旋回的曲线，也即求取所有相邻三条瞬时频率曲线的比值，以获得瞬时频率曲线上超过50％的点的比值，选取比值最接近4:2.5:1的三条曲线作为米兰科维奇旋回的曲线。

上述第三种获取米兰科维奇旋回曲线的方法中，计算每个本征模态函数在目的层段深度范围内的周期数量，用其目的层段深度范围内发育的地质历史时间的时限除以各个本征模态函数的周期数量得到每个周期的地质历史时间周期，选取地质历史时间周期接近25ka、41Ka及100-400ka的三条曲线作为米兰科维奇旋回的三条曲线。

步骤105，根据米兰科维奇旋回曲线，获取可容纳空间变化曲线。

进一步的，从米兰科维奇旋回曲线中任意选择一条瞬时频率曲线，用费希尔图解原理进行分析，以获取可容纳空间变化曲线。具体的，根据瞬时频率曲线所对应的本征模态函数获得旋回数和平均厚度累积偏差，即目的层深度段内单旋回厚度对于平均旋回厚度的累积偏差，然后以旋回数为横坐标、平均厚度累积偏差为纵坐标绘制费希尔图解曲线，该费希尔图解曲线即可容纳空间变化曲线。

步骤106，根据可容纳空间变化曲线，获取可容纳空间上升半旋回和可容纳空间下降半旋回。

本实施例提供的沉积地层旋回划分方法，对获取的反映地层原始沉积环境及岩性变化的资料进行预处理，然后对预处理后的资料进行经验模态分解，获得本征模态函数，之后根据本征模态函数获取米兰科维奇旋回曲线，进一步获取可容纳空间变化曲线，最终根据可容纳空间变化曲线，获取可容纳空间上升半旋回和可容纳空间下降半旋回，从而实现对沉积地层的划分，上述划分方法可以利用地球化学分析数据等在深度上不连续的资料进行计算，而不仅仅局限于应用地球物理测井数据和地震数据，计算过程客观，受人为因素影响较小，是一种可操作性较强的地层旋回划分方法，为识别层序地层界面或沉积旋回划分的界面提供一个量化的判识依据。

实施例二

本实施例是在上述实施例的基础上进行的补充说明。

图2为本发明实施例二提供的沉积地层旋回划分方法的流程示意图；如图2所示，本实施例提供一种沉积地层旋回划分方法，包括：

其中，步骤1011，对地球物理资料进行归一化与标准化处理，以获得标准化资料，其中，地球物理资料包括地球物理测井数据和地震数据。

以岩石样品地球化学分析数据、地球物理测井数据和地震数据为资料，根据以上数据与沉积相、岩相变化的相关性，或已发表文献记载的能够反映地层原始沉积环境及岩性变化的资料进行选择，其中地球物理测井数据与地震数据一般为连续资料，对不同井或不同次地球物理资料进行归一化与标准化处理。

步骤1012，对目的层段包含的标准化资料进行曲线拼接，使得曲线拼接处理后获取的预处理资料连续且能够覆盖完整目的层段深度范围。

对目的层段可能含有的多次地球物理资料进行曲线拼接，使得曲线拼接处理后获取的预处理资料连续且能够覆盖完整目的层段深度范围。

步骤1013，对岩石样品地球化学分析数据进行二次样条差值处理，使得经过二次样条差值处理后获取的预处理资料连续且能够覆盖完整目的层段深度范围。

岩石样品地球化学分析数据属于不连续资料，需要对其进行二次样条差值处理，以使得经过二次样条差值处理后获取的预处理资料连续且能够覆盖完整目的层段深度范围。

步骤102，对预处理资料进行经验模态分解，以获取多个本征模态函数。

预处理资料包括对标准化资料进行曲线拼接后所获得的数据和对岩石样品地球化学分析数据进行二次样条差值处理后的数据。

步骤103，对于每个本征模态函数，对本征模态函数进行变换，以获取与本征模态函数对应的瞬时频率曲线。对每个本征模态函数进行希尔伯特变换，每个本征模态函数均会求得一条瞬时频率曲线。

步骤104，根据瞬时频率曲线，获取米兰科维奇旋回曲线。

步骤105，根据米兰科维奇旋回曲线，获取可容纳空间变化曲线。

步骤102-步骤105具体可参见实施例一中相应的描述，在此不再赘述。

步骤1061，对可容纳空间变化曲线求取导数，以获得可容纳空间变化曲线的极值点。导数为正的部分为可容纳空间上升半旋回，导数为负的部分为可容纳空间下降半旋回。

对可容纳空间变化曲线求取导数，识别极值点作为沉积地层旋回划分的界面；可容纳空间增加的部分即导数为正的部分做为可容纳空间上升半旋回，可容纳空间减少的部分即导数为负的部分做为可容纳空间下降半旋回。

步骤1062，根据可容纳空间变化曲线的极值点将目的层段划分成至少一个包含可容纳空间上升半旋回和可容纳空间下降半旋回的部分。

本实施例提供的沉积地层旋回划分方法，可以利用地球化学分析数据等在深度上不连续的资料进行计算，而不仅仅局限于应用地球物理测井数据和地震数据，计算过程客观，受人为因素影响较小，是一种可操作性较强的地层旋回划分方法。为识别层序地层界面或沉积旋回划分的界面提供一个量化的判识依据。上述沉积地层旋回界面的识别不需要人工经验性判识，而是直接根据与地层成因相联系的资料，如岩石样品地球化学分析数据、地球物理测井数据和地震数据等计算获得，减少了人为因素的影响。

实施例三

本实施例为装置实施例，用于执行上述实施例一中的方法。

图3为本发明实施例三提供的沉积地层旋回划分装置的结构示意图；如图3所示，本实施例提供一种沉积地层旋回划分装置，包括：预处理模块201、本征模态函数获取模块202、瞬时频率曲线获取模块203、米兰科维奇旋回曲线获取模块204、可容纳空间变化曲线获取模块205和划分模块206。

其中，预处理模块201，用于对获取的反映地层原始沉积环境及岩性变化的资料进行预处理，以获得预处理资料，所述资料包括岩石样品地球化学分析数据、地球物理测井数据和地震数据；

本征模态函数获取模块202，用于对所述预处理资料进行经验模态分解，以获取多个本征模态函数；

瞬时频率曲线获取模块203，用于对于每个本征模态函数，对所述本征模态函数进行变换，以获取与所述本征模态函数对应的瞬时频率曲线；

米兰科维奇旋回曲线获取模块204，用于根据所述瞬时频率曲线，获取米兰科维奇旋回曲线；

可容纳空间变化曲线获取模块205，用于根据所述米兰科维奇旋回曲线，获取可容纳空间变化曲线；

划分模块206，用于根据所述可容纳空间变化曲线，获取可容纳空间上升半旋回和可容纳空间下降半旋回。

本实施例是与方法实施例一对应的装置实施例，具体可参见实施例一中的描述，在此不再赘述。

实施例四

本实施例是在实施例三的基础上进行的补充说明，用于执行上述实施例二中的方法。

图4为本发明实施例四提供的沉积地层旋回划分装置，如图4所示，预处理模块201具体包括：第一预处理子模块2011、曲线拼接子模块2012和二次样条差值处理子模块2013。

其中，第一预处理子模块2011，用于对地球物理资料进行归一化与标准化处理，以获得标准化资料，其中，地球物理资料包括地球物理测井数据和地震数据；

曲线拼接子模块2012，用于对目的层段包含的标准化资料进行曲线拼接，使得曲线拼接处理后获取的预处理资料连续且能够覆盖完整目的层段深度范围；

二次样条差值处理子模块2013，用于对岩石样品地球化学分析数据进行二次样条差值处理，使得经过二次样条差值处理后获取的预处理资料连续且能够覆盖完整目的层段深度范围；所述预处理资料包括对所述标准化资料进行曲线拼接后所获得的数据和对所述岩石样品地球化学分析数据进行二次样条差值处理后的数据。

进一步的，米兰科维奇旋回曲线获取模块204具体用于：

获取相邻三条瞬时频率曲线，将在目的层段深度范围内超过第一预设阈值的瞬时频率曲线上的点的比值最接近第二预设阈值的三条曲线作为米兰科维奇旋回的曲线；

或者获取瞬时频率曲线中目的层段深度范围内的主频，将主频的比值最接近第二预设阈值的瞬时频率曲线作为米兰科维奇旋回曲线；

或者获取本征模态函数在目的层段深度范围内的周期数量，并根据目的层段深度范围发育的地质历史时间的时限除以周期数量，以获得瞬时频率曲线的地质历史时间周期，将地质历史时间周期最接近第三预设阈值的三条瞬时频率曲线作为米兰科维奇旋回曲线。

进一步的，可容纳空间变化曲线获取模块205具体用于：

从米兰科维奇旋回曲线中任意选择一条曲线，用费希尔图解原理进行分析，以获取可容纳空间变化曲线。

进一步的，划分模块206具体包括求导子模块2061和划分子模块2062。

其中，求导子模块2061，用于对可容纳空间变化曲线求取导数，以获得所述可容纳空间变化曲线的极值点。导数为正的部分为可容纳空间上升半旋回，导数为负的部分为可容纳空间下降半旋回；

划分子模块2062，用于根据可容纳空间变化曲线的极值点将目的层段划分成至少一个包含可容纳空间上升半旋回和可容纳空间下降半旋回的部分

本实施例是与方法实施例二对应的装置实施例，具体可参见实施例二中的描述，在此不再赘述。

虽然已经参考优选实施例对本发明进行了描述，但在不脱离本发明的范围的情况下，可以对其进行各种改进并且可以用等效物替换其中的部件。尤其是，只要不存在结构冲突，各个实施例中所提到的各项技术特征均可以任意方式组合起来。本发明并不局限于文中公开的特定实施例，而是包括落入权利要求的范围内的所有技术方案。