一种识别沉积地层中的米兰科维奇周期的方法及系统

摘要

本发明公开了一种识别沉积地层中的米兰科维奇周期的方法及系统，其中方法包括以下步骤：对能反映泥质量含量的测井曲线作Morlet连续小波变换，得到不同小波尺度对应的小波分解曲线，分别统计各小波尺度对应的小波分解曲线中包含的周期个数，作小波尺度与小波分解曲线周期个数关系图，在图中搜索小波尺度比值和小波分解曲线周期个数比值分别与米兰科维奇周期不同轨道参数周期之间的周期比值和频率比值同时吻合的点组，然后在引入的地质约束条件下，筛选得到的点组对应的小波分解曲线即为相应的米兰科维奇周期。该方法与系统相对应，能够解决当米兰科维奇旋回在反映泥质含量的测井曲线频率分解后的能量谱中不对应能量峰值时的米兰科维奇周期识别问题。

说明书

技术领域

本发明涉及地球物理勘探领域，具体涉及一种识别沉积地层中的米兰科维奇周期的方法及系统。

背景技术

米兰科维奇周期最早由南斯拉夫学者米兰科维奇(Milankovitch)于1920年提出，是已被天文计算和地质纪录证实的地球公转轨道参数(偏心率、黄赤交角和岁差)的变化周期，主要为2±万年的岁差周期、4±万年黄赤交角周期、10±万年的地球轨道短偏心率周期和40±万年的地球轨道长偏心率周期。米兰科维奇周期能在沉积地层中记录的原因是地球轨道参数变化引起地球日照量的变化，进而引起地球气候变化，最终导致盆地沉积变化，因此在地层沉积中记录了与米兰科维奇周期相同的沉积旋回，该沉积旋回被称为米兰科维奇(Milankovitch)旋回。地层记录中常见的米兰科维奇旋回为2±万年的岁差旋回、4±万年黄赤交角旋回和10±万年的短偏心率旋回，而40±万年的长偏心率旋回则较少见。

依据沉积地层中米兰科维奇(Milankovitch)旋回的年代地层意义可对地层进行定年。由于天文轨道周期具有稳定性，从地层定年层段内的沉积记录中识别和提取出米兰科维奇周期，根据计算层段内的米兰科维周期个数和每个周期对应的时间长度可得到分析层段的持续地质时间，并可根据米兰科维奇周期在地层中的分布位置，对地层进行相对定年，建立相对地层年代格架，计算沉积速率，同时，如果已知起算点处的地质年龄，则可计算地层的绝对地质年龄。

识别地层中的米兰科维奇周期主要是依据米兰科维奇周期中各地球公转轨道参数变化周期之间的比例关系，若在地层沉积中识别出这种比例关系的沉积周期，也即识别出了米兰科维奇周期。为寻找地层中这种比例关系的周期，在露头区可以采用微层划分、厚度测量、费舍尔图解的方法；覆盖区则可利用钻井测井资料进行米兰科维奇周期识别。

利用钻井测井资料进行米兰科维奇周期识别的现有的方法主要有基于傅立叶变换的频谱分析法和基于小波变换的小波分析法两大类，其中，小波分析法主要是通过对能反映泥质含量的测井曲线作连续小波变换，对得到的小波系数矩阵使用模平均值、模范数或小波方差等方法进行主周期分析，在分析出的所有主周期中寻找与米兰科维奇各周期之间的比值关系符合的主周期即为米兰科维奇周期。但影响盆地沉积的因素很多，如盆地构造波动。如果盆地沉积演化过程中，米兰科维奇周期控制下的沉积过程被其它因素控制下的沉积过程干扰，则米兰科维奇旋回在频率分解的能量谱中可能不对应能量峰值，此时主周期分析法将无法根据各主周期之间的比值关系有效识别出米兰科维奇周期，如在低纬度地区。

发明内容

针对现有技术中，当米兰科维奇旋回在反映地层泥质含量的测井曲线小波时频分析结果中不表现为能量峰值或主周期时，则基于钻井测井资料小波时频分解的主周期分析法将不能有效识别出米兰科维奇周期的技术问题，本发明提出一种不同于现有主周期分析法识别米兰科维奇周期的方法，解决当主周期分析法不适用时的兰科维奇周期识别问题。

本发明提供一种识别沉积地层中的米兰科维奇周期的方法，包括以下步骤：

S10：获得可以反映泥质含量的测井曲线；

S20：对步骤S10中的测井曲线作连续小波变换，得到不同小波尺度对应的小波分解曲线；

S30：统计步骤S20中各小波尺度对应的小波分解曲线中包含的周期个数，作小波尺度与小波分解曲线周期个数对应关系图；

S40：在步骤S30得到的小波尺度与小波分解曲线周期个数对应关系图中，搜寻同时符合米兰科维奇不同轨道参数周期之间的周期比值和频率比值的点组；

S50：对步骤S40得到的点组对应的小波分解曲线按照米兰科维奇周期的时间意义得到相应的地层年龄；

S60：引入地质约束条件，将步骤S50中得到的地层年龄与所述地质约束条件对比，确定米兰科维奇周期。

根据本发明的一个实施例，步骤S10中采用的测井曲线为自然伽马曲线。

根据本发明的一个实施例，步骤S20中所述连续小波变换为Morlet连续小波变换，其中Morlet连续小波定义式如下：

其中：Ψ₀为Morlet小波，t为变量。

根据本发明的一个实施例，步骤S20中所述的连续小波变换，根据下式进行：

其中：

WT_f为f的连续小波变换，

Ψ为小波函数，

t为变量，

a为尺度因子，

b为平移因子。

根据本发明的一个实施例，步骤S30包括以下子步骤：

S31：对步骤S20中的小波分解曲线作光滑去噪处理，剔除曲线中的扰动噪音；

S32：应用极值点法计算步骤S31中的小波分解曲线中包含的周期个数；

S33：以小波尺度为横坐标，周期个数为纵坐标，将步骤S32中得到的不同小波尺度对应分解曲线包含的周期个数投在图中，得到小波尺度与小波分解曲线周期个数对应关系曲线。

根据本发明的一个实施例，步骤S40中，在所述小波尺度与小波分解曲线周期个数对应关系图上，以小波尺度比值与所述周期比率对应，周期个数比值与所述频率比率对应。

根据本发明的一个实施例，步骤S50中的地层年龄获得方法为以所述周期对应的周期个数和所述周期对应的时间相乘得到地层年龄。

根据本发明的一个实施例，步骤S60中的所述的地质约束条件为地层年代格架或井底地层地质年代范围。

本发明的另一方面，还公开一种识别沉积地层中的米兰科维奇周期的系统，包括：

曲线获得模块，所述曲线获得模块用于获得可以反映泥质含量的测井曲线；

小波变换模块，所述小波变换模块用于对曲线获得模块中得到的测井曲线作连续小波变换，得到不同小波尺度对应的小波分解曲线；

周期统计模块，所述周期统计模块用于对小波变换模块得到的各小波尺度对应的小波分解曲线中包含的周期个数，作小波尺度与小波分解曲线周期个数对应关系图；

点组搜索模块，所述点组搜索模块用于在周期统计模块得到的小波尺度与小波分解曲线周期个数对应关系图中，搜寻同时符合米兰科维奇不同轨道参数周期之间的周期比值和频率比值的点组；

地层年龄计算模块，所述地层年龄计算模块用于对点组搜索模块得到的点组对应的小波分解曲线按照米兰科维奇周期的时间意义进行定年计算相应的地层年龄；和

米兰科维奇周期确定模块，所述米兰科维奇周期确定模块用于引入地质约束条件，将地层年龄计算模块中得到的地层年龄与所述地质约束条件对比，确定米兰科维奇周期。

本发明可以从反映泥质含量的钻井测井曲线中识别出米兰科维奇周期，特别是当米兰科维奇周期在小波分析时频分解的能量谱中不表现为能量峰值而无法通过主周期分析法进行识别的情况，如：在赤道等低纬度地区分布的沉积盆地中，米兰科维奇周期的部分周期值在地层自然伽马测井曲线经Morlet连续小波时频分解后的能量谱中，往往不表现为能量峰值，此时使用主周期分析法将无法有效识别米兰科维奇周期，但使用本发明提出的方法可以克服上述局限，有效识别出米兰科维奇周期。

附图说明

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

图1是本发明实施例一中的一种识别沉积地层中的米兰科维奇周期的方法流程图；

图2是本发明实施例一中的步骤S30中的小波尺度与小波分解曲线所含周期个数对应关系曲线图；

图3-a是本发明实施例一中的2±万年的岁差周期的小波分解曲线；

图3-b是本发明实施例一中的4±万年的黄赤交角周期的小波分解曲线；

图3-c是本发明实施例一中的10±万年的地球轨道偏心率周期的小波分解曲线；

图4是本发明实施例二中的识别沉积地层中的米兰科维奇周期的系统的结构图。

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

具体实施方式

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

实施例一

图1所示为本实施例中的一种识别沉积地层中的米兰科维奇周期的方法的流程图。如图1所示，该方法具体包括以下步骤：

S10：获得可以反映泥质含量的测井曲线；

S20：对步骤S10中的测井曲线作连续小波变换，得到不同小波尺度对应的小波分解曲线；

S30：统计步骤S20中各小波尺度对应的小波分解曲线中包含的周期个数，作小波尺度与小波分解曲线周期个数对应关系图；

S40：在步骤S30得到的小波尺度与小波分解曲线周期个数对应关系图中，搜寻同时符合米兰科维奇不同轨道参数周期之间的周期比值和频率比值的点组；

S50：对步骤S40得到的点组对应的小波分解曲线按照米兰科维奇周期的时间意义得到相应的地层年龄；

S60：引入地质约束条件，将步骤S50中得到的地层年龄与所述地质约束条件对比，确定米兰科维奇周期。

具体地，本实施例中，在步骤S10中，采用的测井曲线为自然伽马曲线(GR)。测井曲线的具体制作方法采用现有技术中常用的方法对即可，在此不再赘述。

具体地，本实施例中，在步骤S20中，对能反映泥质含量的测井曲线(如自然伽马(GR)测井曲线)作Morlet连续小波变换，得到不同小波尺度对应的小波分解曲线。

其中，Morlet连续小波的定义形式如式(1)：

式(1)中的Ψ₀为Morlet小波，t为变量，

具体的，t相当于函数变量，给定t一系列任意值，根据式(1)可计算出其对应的Ψ₀值，Ψ₀描绘的曲线即为Morlet小波；

其中Morlet连续小波变换采用式(2)进行：

式(2)中，WT_f为f的连续小波变换，t为变量，Ψ为小波函数，a为尺度因子，b为平移因子。

具体地，本实施例中，在步骤S30中，还包括以下子步骤，统计小波分解曲线中包含的周期个数前需要对小波分解曲线作光滑去噪处理，剔除曲线中的扰动噪音，然后应用极值点法计算曲线中包含的周期个数；以小波尺度为横坐标，周期个数为纵坐标，将不同小波尺度对应分解曲线包含的周期个数投在坐标图中得到小波尺度与小波分解曲线周期个数对应关系曲线图，具体的示例如图2所示。

具体地，本实施例中，在步骤S40中，在所述小波尺度与小波分解曲线周期个数对应关系图上，以小波尺度比值与所述周期比率对应，周期个数比值与所述频率比率对应。以表1所示米兰科维奇周期中不同轨道参数周期之间的周期比率和频率比率为参考，在周期个数曲线图上以小波尺度比值与周期比率对应，周期个数比值与频率比率对应，在周期个数曲线上搜索同时满足表1中的米兰科维奇周期各轨道参数之间的周期比率和频率比率，(在本申请中，例如在表1中Ka为时间单位千年，即1Ka＝1千年)。搜索过程中需要设置合理的误差范围，只要在误差范围内，则可认为该点是有效的。本发明实施例在步骤S40得到的点组见表2所示。

表1

表2

具体地，本实施例中，在步骤S50中依据步骤S40得到的点组对应的10万年偏心率周期包含的周期个数计算地层年龄，计算方法为10万年周期与周期个数相乘，该实施例的全井段(井口——井底)10万年周期个数为114个，则经计算井底地层年龄为11.4Ma(本发明中Ma为时间单位，百万年)；中新统顶界之上井段(井口——中新统顶界)10万年周期个数为52个，则经计算中新统顶界处的地层年龄为5.2Ma(百万年)。

具体地，本实施例中，本实施例在步骤S60中引入钻井分层界面和井底地层所处的地质年代范围为约束条件，与步骤S50计算得到的井底地层年龄和钻井分层界面处的年龄对比。本实施例计算得到的钻井分层中新统顶界处的年龄为5.2Ma，井底地层年龄为11.4Ma，与国际地层年代表(2008)中的中新统顶界年龄5.332Ma，及井底地层位于上中新统底-中中新统顶部(地质年龄5.3Ma～11.6Ma)基本符合。因此，该点组对应的小波分解曲线(如图3-a、图3-b和图3-c所示)即为该井地层中的米兰科维奇旋回对应的米兰科维奇周期曲线，依据该曲线可以在地层中划分出米兰科维奇旋回。

图3-a对应的为2±万年的岁差周期的小波分解曲线，图3-b对应的是4±万年的黄赤交角周期的小波分解曲线，图3-c对应的为10±万年的地球轨道偏心率周期的小波分解曲线。

实施例二

如图4所示，为本实施例中的识别沉积地层中米兰科维奇周期的系统的结构图，由图所示可以看出，该系统包括：

曲线获得模块10，用于获得可以反映泥质含量的测井曲线；

小波变换模块20，用于对曲线获得模块中得到的测井曲线作连续小波变换，得到不同小波尺度对应的小波分解曲线；

周期统计模块30，用于对小波变换模块得到的各小波尺度对应的小波分解曲线中包含的周期个数，作小波尺度与小波分解曲线周期个数对应关系图；

点组搜索模块40，用于在周期统计模块得到的小波尺度与小波分解曲线周期个数对应关系图中，搜寻同时符合米兰科维奇不同轨道参数周期之间的周期比值和频率比值的点组；

地层年龄计算模块50，用于对点组搜索模块得到的点组对应的小波分解曲线按照米兰科维奇周期的时间意义得到相应的地层年龄；和

米兰科维奇周期确定模块60，用于引入地质约束条件，将地层年龄计算模块中得到的地层年龄与所述地质约束条件对比，确定米兰科维奇周期。

在一些实施例中，与所述方法对应地设置相应的处理模块，本领域技术人员员根据本发明中的披露的方法可以对系统进行改进和进一步的设计，在此不再赘述。

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