CN109031443B - 基于米兰科维奇旋回定年的剥蚀量恢复方法 - Google Patents

Abstract

本发明涉及一种基于米兰科维奇旋回定年的剥蚀量恢复方法，涉及地质学盆地分析技术领域，解决了现有技术中存在的局限性较大、主观因素较多而使误差较大以及操作复杂、不容易实现的技术问题。本发明的方法包括确定测量区域内标志层的步骤；识别测量区域的米兰科维奇周期的步骤；获得顶界年龄和底界年龄的步骤；获得不整合下伏残余地层平均沉积速率的步骤；获得被剥蚀地层曾经拥有的时间的步骤；以及获得不整合剥蚀量的步骤。通过采用上述方法，能够获得误差小、可信度高的不整合剥蚀量，且测量的成本低、适用的范围广。

Description

基于米兰科维奇旋回定年的剥蚀量恢复方法

技术领域

本发明涉及地质学盆地分析技术领域，特别地涉及一种基于米兰科维奇旋回定年的剥蚀量恢复方法。

背景技术

地层剥蚀是沉积盆地中普遍存在的现象，如果剥蚀量不大，则对油气生成、运移和聚集的影响可以不必考虑；但如果有过较大的剥蚀，则会对盆地中油气的生成、运移和聚集等产生影响，这时就要恢复剥蚀量。恢复地层剥蚀厚度对古构造演化史分析以及油气生成、运聚和保存史分析，都具有重要意义。

目前采用的剥蚀量恢复的方法有以下几种：

第一种是地层对比法，该方法是将研究区内被剥蚀层段与邻区未被剥蚀层段进行对比，求得被剥蚀岩层的厚度，这时可以考虑厚度递减的原则或采用其它外推法进行校正。由于该方法依据角度不整合计算地层剥蚀量，因此，该方法只能反映这个地区的最小剥蚀量，而对于平行不整合的区域该方法无能为力，因此该方法存在局限性。

第二种方法是沉积速率分析法，该方法需要两个参数，一是地层剥蚀的时限，二是地层剥蚀速率，这两个参数相乘及可以得到地层剥蚀量。该方法将地层剥蚀的时限设置为整个沉积间断时限的一半，而剥蚀速率则用沉积速率来代替。但是这种方法的人为因素较大，计算出来的剥蚀量并不可靠。

图1-a为某区域第一时间段内地层A、B和C之间整合且未变形的状态示意图；图1-b为某区域第二时间段内地层A、B和C之间整合但褶皱隆起的变形的状态示意图；图4为某区域第三时间段内地层A、B和C褶皱变形后被快速侵蚀夷平的状态示意图；图1-c为某区域第四时间段在夷平面之上新沉积的地层A、B和C之间的不整合接触关系的状态示意图。如图1-a至1-d所示，随着时间由第一时间段更新到第四时间段，图1-d显示了某区域内风化壳不发育型不整合的形成过程。其中，1为地层A，2为地层B，3为地层C，4为地层D，5为钻井A，6为钻井B。

如图1-a至1-d所示，图中显示了4套地层，其中地层A、地层B和地层C为构造运动发生之前沉积的层系，地层D为构造运动发生之后沉积的层系。

从第一阶段到第四阶段为时间的更新，随着时间的更新，地层A、地层B和地层C均遭受了剥蚀；很显然，在应用沉积速度法恢复钻井A处地层A的剥蚀量时，将剥蚀时间定为第四时间段中地层A的顶面年龄与地层D的底面年龄之差的一半是不合理的，同样的道理，在应用沉积速度法恢复钻井B处地层B剥蚀量时，将剥蚀时间定为第四阶段中地层B的顶面年龄与地层C的底面年龄之差的一半的是不合理的，因此根据该参数获得的剥蚀量误差较大。

第三种方法是测井曲线法，由于在正常压实情况下，页岩压实与上覆的负荷或埋深有关，孔隙度是页岩压实程度的度量，而声波测井资料直接反映了页岩压实程度的大小。如果泥岩埋藏压实后又被抬升，则正常压实趋势会被破坏，不整合面上下声波会出现跳跃，该方法是根据这一异常压实现象可以推算剥蚀量。但是在很多情况下，不整合面上下地层的声波时差并不存在跳跃异常现象，原因在于上覆地层的厚度如果远远大于剥蚀的厚度，尤其是盆地早期演化过程中形成的剥蚀量，因此这种方法也存在局限性。

第四种方法是流体包裹体法，利用流体包裹体计算地层剥蚀厚度，国外还未成熟，国内尚在开始阶段。其基本原理是由于流体包裹体记载了它们所经历的整个受热地质历史中不同时期沉积物所处的温度、压力等热力学条件的信息，因此在连续沉积过程中，捕获的包裹体温度(或压力)与埋藏深度的对应数值一般呈良好的线性关系。然而，在侵蚀不整合面上下两边的地层中，它们的温度或压力系统往往不同，因此温度与埋藏深度曲线在侵蚀不整合面之处往往表现为曲线明显的温度跃变现象。在温度-深度坐标图上，只需将剥蚀面以下深度、温度(或压力)对应数值的点用回归方法联结成的直线，向上延伸至地表温度处，即为古地表温度。对应于这一温度坐标的标高面就是古地表面，由剥蚀面至古地表面的距离就是地层剥蚀厚度。但是该方法需要古地表温度参数，而该参数较难求取，此外，不整合界面上下地层中流体包裹体温度分布跨度较大，因此该方法获得的剥蚀量误差较大。

第五种方法是磷灰石裂变径迹分析法，磷灰石裂变径迹分析法是近十几年发展起来的恢复沉积盆地热史的一种新方法。该方法主要建立在磷灰石所含的U²³⁸自发裂变产生的径迹，即裂变径迹，在地质历史时间内受温度作用而发生退火行为这一化学动力学原理基础之上。磷灰石裂变径迹分析法虽然可以同时计算出地层剥蚀的时间和剥蚀量，但是该方法仅适合于埋藏不是很大、磷灰石未经历过退火的碎屑岩，因此对于埋藏较深的地层不整合剥蚀量问题，较难应用该方法求取。此外，该方法样品来自沉积盆地的碎屑岩地层，不适用于碳酸盐地层。因此该方法也存在一定的局限性。

第六种方法是沉积盆地波动分析法，该方法是在开展地层古生物、不整合及沉积环境研究的基础上，建立各研究区沉积速率-地质年代直方图，借助滑动窗口对沉积速率直方图进行滤波处理，找到能够代表该区沉积—剥蚀过程的周期以及振幅，最终建立波动方程，然后利用波动方程预测地层缺失时间段内地层的剥蚀量。但是该方法在建立波动方程时有很多人为因素，需要反复调试，费时费力，并需要保证沉积间断处沉积与剥蚀的平衡，因此实际操作复杂，不容易实现。

发明内容

本发明提供一种不同采集参数基于米兰科维奇旋回定年的剥蚀量恢复方法，用于解决现有技术中存在的局限性较大、主观因素较多而使误差较大以及操作复杂、不容易实现的缺陷。

本发明提供一种不同采集参数基于米兰科维奇旋回定年的剥蚀量恢复方法，包括以下步骤：

S10：确定测量区域内可进行区域对比且可标记绝对年龄的标志层；

S20：依据测井曲线识别测量区域的米兰科维奇周期；

S30：以所述标志层为地层时间起算点，以所述米兰科维奇周期为步长进行标记，分别获得测量区域残余地层的顶界年龄T_顶和底界年龄T_底；

S40：测量残余地层的沉积厚度D，获得残余地层的平均沉积速率V，残余地层的平均沉积速率V根据下列定义式进行：

S50：确定与所述残余地层相对应的标准地层的顶界年龄T_标；获得被剥蚀地层曾经拥有的时间T_剥，被剥蚀地层曾经拥有的时间T_剥根据下列定义式进行：

T_剥＝T_顶-T_标；

S60：获得不整合剥蚀量Q，不整合剥蚀量Q根据下列定义式进行：

Q＝T_剥×V。

根据本发明的一个实施例，步骤S50中，通过与相邻区域地层对比以及与国际地质年表对比，获得标准地层的顶界年龄T_标。

根据本发明的一个实施例，所述沉积厚度D为残余地层的顶界深度与底界深度之差。

根据本发明的一个实施例，步骤S10中，通过发育于所述标志层中的古生物资料与国际地质年代表对比，获得所述标志层的绝对年龄。

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

根据本发明的一个实施例，步骤S20中，通过对测井曲线作傅里叶变换或Morlet连续小波变换以识别米兰科维奇周期，并得到米兰科维奇周期曲线。

根据本发明的一个实施例，步骤30中，通过对所述米兰科维奇周期曲线进行标定，获得所述残余地层的顶界年龄T_顶和底界年龄T_底，进行标定时使所述标志层的底界与所述米兰科维奇周期曲线的波谷相对应。

根据本发明的一个实施例，步骤S30中，所述米兰科维奇周期为10万年周期、4万年周期或2万年周期。

根据本发明的一个实施例，经过所述Morlet连续小波变换后，获得尺度因子-模值曲线，所述尺度因子-模值曲线的横坐标为尺度因子，纵坐标为模值；

根据尺度因子-模值曲线获得模值的极大值处对应的Morlet连续小波的周期个数，根据国际地质年表获得测量区域的理论时间跨度；

所述米兰科维奇周期为最接近理论时间跨度与周期个数之比的周期。

根据本发明的一个实施例，所述标志层为最大海泛期形成的凝缩层或火山喷发形成的广泛分布的火山凝灰岩层。

与现有技术相比，本发明的优点在于：

(1)通过地层年代学研究找回被剥蚀地层曾经拥有的时间T_剥，即可以恢复地层剥蚀量，因此无需要花费高额的分析测试费用，其成本底、可操作性好。

(2)通过对米兰科维奇周期曲线进行标定，从而获得不整合下伏残余地层平均沉积速率V，并恢复地层剥蚀量，因此在进行分析时不会受到岩石类型的限制，因此其应用没有局限性，应用更广泛。

(3)由于上述获得地层剥蚀量的两个参数被剥蚀地层曾经拥有的时间T_剥和下伏残余地层平均沉积速率V均是经过实际测量和计算获得，因此避免了人为因素的影响，使地层剥蚀量的误差小、可信度高。

附图说明

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

图1-a为某区域第一时间段内地层A、B和C之间整合且未变形的状态示意图；

图1-b为某区域第二时间段内地层A、B和C之间整合但褶皱隆起的变形的状态示意图；

图1-c为某区域第三时间段内地层A、B和C褶皱变形后被快速侵蚀夷平的状态示意图；

图1-d为某区域第四时间段在夷平面之上新沉积的地层A、B和C之间的不整合接触关系的状态示意图；

图2为本发明实施例中某区域的自然伽玛曲线与米兰科维奇周期曲线对比图；

图3为本发明实施例中某区域时间尺度的地层剖面示意图；

图4为本发明实施例中塔里木盆地巴楚地区巴探5井的井下古生界古地理位置图；

图5-a为本发明实施例中巴探5井寒武系自然伽玛原始数据图；

图5-b为本发明实施例中巴探5井寒武系Morlet连续小波变换时频色谱图；

图6为本发明实施例中巴探5井寒武系Morlet连续小波变换的尺度因子-模值曲线图；

图7-a为本发明实施例中巴探5井寒武系自然伽玛曲线图；

图7-b为本发明实施例中巴探5井寒武系2万年米兰科维奇周期曲线图；

图7-c为本发明实施例中巴探5井寒武系4万年米兰科维奇周期曲线图；

图7-d为本发明实施例中巴探5井寒武系10万年米兰科维奇周期曲线图；

图8-a为本发明实施例中巴探5井奥陶系自然伽玛原始数据图；

图8-b为本发明实施例中巴探5井寒武系Morlet连续小波变换时频色谱图；

图9为本发明实施例中巴探5井奥陶系Morlet连续小波变换的尺度因子-模值曲线图；

图10-a为本发明实施例中巴探5井奥陶系自然伽玛曲线图；

图10-b为本发明实施例中巴探5井奥陶系2万年米兰科维奇周期曲线图；

图10-c为本发明实施例中巴探5井奥陶系4万年米兰科维奇周期曲线图；

图10-d为本发明实施例中巴探5井奥陶系10万年米兰科维奇周期曲线图；

图11为巴探5井寒武系-奥陶系鹰山组岩性-深度剖面与Morlet连续小波变换的米兰科维奇周期曲线的对应关系示意图；

图12为巴探5井寒武系-奥陶系鹰山组依据米兰科维奇周期标定的岩性-时间剖面示意图。

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

附图标记：

1-地层A； 2-地层B； 3-地层C；

4-地层D； 5-钻井A； 6-钻井B；

7-不整合间断层； 8-巴探5井。

具体实施方式

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

本发明提供一种基于米兰科维奇旋回定年的剥蚀量恢复方法，依靠识别米兰科维奇旋回从而对测量区域进行标定。

20世纪初，南斯拉夫气象学家米兰科维奇(Milankovitch)为解释第四纪大冰期的成因，研究了天文学近60万年来轨道三要素变化规律的资料，结合一系列数学公式计算出在过去60万年内不同纬度地区由于地球轨道三要素(偏心率，倾角或黄赤交角，岁差)的变化引起夏季和冬季所在上下半年太阳辐射量曲线的变化，并将这种变化曲线与第四纪气候变迁的经典模式进行了比较，发现二者存在许多相似性。因而米兰科维奇认为，太阳对地球辐射量的变化导致了第四纪气候的变迁，即太阳对地球辐射量越大，地球上的温度就会越高，同时也认为，北纬65度地区天文辐射夏半年总量的多少决定了全球冰期的消长，这一假设就是米兰科维奇理论。该假说很长一段时间没有被地质学家接受，直到20世纪70年代，海斯(Hays)等依据更新世深海钻孔中有孔虫壳的氧同位素数据的频率分析，指示出了全球冰量的变化，同时期海洋表面的温度也根据来自南半球放射虫的组合被估算出来，这二者很好的证实了晚更新世气候变化与轨道旋回有相同的变化周期，从而证实了米兰科维奇的理论，并引起了天文古气候学家和地质学家两大领域对米兰科维奇理论研究的重视。

米兰柯维奇周期主要有3个，即10万年的地球公转偏心率变化周期、4.1万年的地轴倾斜度变化周期以及2万年左右的岁差变化周期，他们之间存在着一定的比例关系，如10万年周期中包含了2.5个4万年的周期和5个2万年的周期。如果在地层旋回中识别出这种比例关系的旋回，也即识别出了米兰柯维奇周期。为寻找地层中这种比例关系的旋回，可利用测井数据进行数学变换实现米兰科维奇周期的寻找。目前用于识别米兰科维奇旋回的数学变换方法主要有两种，即傅立叶变换方法和小波变换方法。

80、90年代前期人们采用序列的谱分析方法，一般采用傅立叶变换或最大嫡谱法求取功率谱或振幅谱，通过分析谱峰频率的倒数来识别米兰科维奇旋回。但这种方法对数据的频域分析不足，任一频率分量都是对整个剖面定义区间上的积分，只能反映整个处理剖面的平均波谱，不能反映波谱随深度的变化细节，即不具备位置和频率的“定位”功能，因此应用效果也有其局限性。

小波变换(wavelet transform，WT)是一种新的变换分析方法，它继承和发展了短时傅立叶变换局部化的思想，同时又克服了窗口大小不随频率变化等缺点，能够提供一个随频率改变的“时间-频率”窗口，是进行信号时频分析和处理的理想工具。它的主要特点是通过变换能够充分突出问题某些方面的特征，能对时间(空间)频率的局部化分析，通过伸缩平移运算对信号(函数)逐步进行多尺度细化，最终达到高频处时间细分，低频处频率细分，能自动适应时频信号分析的要求，从而可聚焦到信号的任意细节，解决了傅立叶变换的困难问题，成为继傅立叶变换以来在科学方法上的重大突破。本发明即是以该方法的基础上，进一步通过插值处理将地层-深度剖面转化为地层-时间剖面。

具体来说，本发明提供的基于米兰科维奇旋回定年的剥蚀量恢复方法包括以下步骤：第一步，确定测量区域内可进行区域对比且可标记绝对年龄的标志层。这是最基础的前期工作，该标志层必须在测量区域内广泛分布，例如可选择最大海泛期形成的凝缩层以及火山喷发形成的广泛分布的火山凝灰岩层，火山岩类岩石由于可以同位素定年，是最好的标志层。沉积岩标志层的绝对年龄可通过发育于其中的古生物资料与国际地质年代表对比而获得。

第二步，依据测井曲线识别测量区域的米兰科维奇周期。

进一步地，所采用的测井曲线为自然伽玛曲线。之所以可以应用自然伽玛来识别米兰科维奇周期，是因为在气候温暖潮湿时期，沉积岩标志层中泥质成分会增加，相应地自然伽玛高值就会在测井曲线上反映出来。

进一步地，本发明的一个实施例是通过对自然伽玛曲线作Morlet连续小波变换以识别米兰科维奇周期，并得到米兰科维奇周期曲线。

第三步，以标志层为地层时间起算点，以米兰科维奇周期为步长进行标记，分别获得测量区域残余地层的顶界年龄T_顶和底界年龄T_底。

其中，米兰科维奇周期为10万年周期、4万年周期或2万年周期。

例如以米兰科维奇10万年周期为步长，依据米兰科维奇周期曲线分别向上或向下进行标定，从而获得测量区域残余地层的顶界年龄T_顶和底界年龄T_底。其中，向上为沿靠近地层的方向，向下为沿远离地层的方向。

图2为本发明实施例中某区域的自然伽玛曲线与米兰科维奇周期曲线对比图；如图2所示，标志层位于岩性柱的顶部，其底界时间为65百万年(65Ma)。最右侧的曲线为经过小波变换处理获得的米兰科维奇周期10万年周期曲线，依据该周期曲线可以划分出M1、M2、M3、M4和M5共计5个米氏旋回层。标志层的底界与米兰科维奇周期的波谷对应，从标志层底界年龄65Ma开始，可向下依次将每一个米氏旋回层的起始时间标记出来，如图6中M1底界年龄标记为65.1Ma，直至将M5底界年龄为65.5Ma标记出来。

进一步地，经过所述Morlet连续小波变换后，获得尺度因子-模值曲线，所述尺度因子-模值曲线的横坐标为尺度因子，纵坐标为模值；

根据尺度因子-模值曲线获得模值的极大值处对应的Morlet连续小波的周期个数，根据国际地质年表获得测量区域的理论时间跨度；

米兰科维奇周期为最接近理论时间跨度与周期个数之比的周期。

第四步，测量残余地层的沉积厚度D，获得残余地层的平均沉积速率V，残余地层的平均沉积速率V根据下列定义式进行：

进一步地，可测量获得残余地层的顶界深度与底界深度，沉积厚度即为二者之差。

第五步，确定与残余地层相对应的标准地层的顶界年龄T_标；获得被剥蚀地层曾经拥有的时间T_剥，被剥蚀地层曾经拥有的时间T_剥根据下列定义式进行：

T_剥＝T_顶-T_标；

进一步地，标准地层的顶界年龄T_标可通过与相邻区域地层对比以及与国际地质年表对比而获得。

图3为本发明实施例中某区域时间尺度的地层剖面示意图；如图3所示，图中自下而上，显示了下伏残余地层(即地层A、B和C)、不整合间断层3和上覆地层(即地层D)，及其所经历的时间，即T_标＝T_顶-T_剥，因此可计算得到T_剥。

第六步，获得不整合剥蚀量Q，不整合剥蚀量Q根据下列定义式进行：

Q＝T_剥×V。

下面以塔里木盆地巴楚地区巴探5井奥陶系中奥陶统鹰山组(O2ys)为例，对本发明的方法进行说明。

图4为本发明实施例中塔里木盆地巴楚地区巴探5井的井下古生界古地理位置图；如图4所示，序号8为巴探5井。塔里木盆地是一多旋回演化的叠合盆地，早古生代处于盆地演化的最早期旋回。从寒武纪到中奥陶世，塔里木盆内发育了伸展背景下的克拉通内碳酸盐台地和斜坡-深水拗拉槽，周边处于大洋扩张和被动大陆边缘环境，盆地内部古地理格局近南北向展布。在此伸展构造背景下，寒武纪自下而上沉积了玉尔吐斯组(∈1y)、肖尔布拉克组(∈1x)、吾松格尔组(∈1w)、沙依里克组(∈2s)、阿瓦塔格组(∈2a)、丘里塔格群下亚群(∈3ql)。继寒武纪之后同样处于伸展期的早中奥陶世，自下而上沉积了下奥陶统的蓬莱坝组(O1p)、中奥陶统的鹰山组(O2ys)和中奥陶统一间房组(O2yj)。

中奥陶世末区域大地构造背景下的伸展构造体制向挤压构造体制的转换，在中奥陶统一间房组(O2yj)之上沉积了上奥陶统的恰尔巴克组(O3q)、良里塔格组(O3l)和桑塔木组(O3s)。这套上奥陶统层系在坳陷区保存完整，但是在隆起区因剥蚀而导致地层保存不完整。

处于隆起区的巴探5井不仅缺失整个上奥陶统的恰尔巴克组(O3q)和中奥陶统一间房组(O2yj)，还缺失了部分鹰山组(O2ys)，在地层剖面上良里塔格组(O3l)与鹰山组(O2ys)直接接触，鹰山组顶界深度和底界深度分别为3210m和3760m。

因此本实施例针对鹰山组剥蚀量的恢复过程如下：

第一步，确定巴探5井可区域对比且可标记绝对年龄的标志层。

在塔里木盆地内，寒武系的沙依里克组灰岩分布广泛，为整个盆地全区对比的标志层，依据国际地质年表，该标志层的底界年龄为510Ma。

第二步，依据自然伽玛曲线识别巴探5井寒武系-奥陶系米兰科维奇周期。

为了更好地完成对巴探5井寒武系-奥陶系米兰科维奇周期的识别，对巴探5井的寒武系和奥陶系的自然伽玛曲线分别作Morlet连续小波变换。

首先，依据自然伽玛识别巴探5井寒武系米兰科维奇周期。

寒武系最下部的肖尔布拉克组的底界深度5805m，寒武系最上部的下丘里塔格组的顶界深度为4214m。巴探5井完钻井深为5822m，进入前寒武系。在对寒武系层段的自然伽玛曲线进行Morlet连续小波变换处理时，向上和向下分别延长到井深4196.25m和井深5811.875m，图5-a为本发明实施例中巴探5井寒武系自然伽玛原始数据图，如图5-a所示，对这井深4196.25m和井深5811.875m之间共计12926个自然伽玛数据进行连续Morlet连续小波变换。图5-b为本发明实施例中巴探5井寒武系Morlet连续小波变换时频色谱图，如图5-b所示，在进行连续小波变换时，取小波的最大长度为100，最小为0.1，滑动步长为0.1。在图5-b中，横坐标为自然伽玛深度对应的数据；纵坐标为连续小波变换的不同尺度，其尺度从最小0.1到100；其颜色代表连续小波变换的系数大小，冷色为低值，暖色为高值，图5-b显示的色谱图为1000×12926二维数据矩阵直观图形显示，通过该色谱图可以粗略地看出自然伽玛主周期以及关键变化点。

图6为本发明实施例中巴探5井寒武系Morlet连续小波变换的尺度因子-模值曲线图，为了更精确地找到自然伽玛的主周期，可以对不同尺度的Morlet连续小波变换数据进行分析，即计算每一个尺度下的Morlet连续小波的模值，1000×12926小波系数矩阵中每一行即代表了每一尺度的小波系数，以尺度因子为横坐标，以计算得到的1000个模值为纵坐标，可以获得如图6所示的尺度因子-模值曲线图。

在图6中，尺度因子为333处的模值为一极大值，它所对应的Morlet连续小波反映有292个周期；而数据处理的寒武系大体的时间跨度为从肖尔布拉克期开始，至下丘里塔格期结束，与国际地质年表对比，其历经的时间为从521Ma开始，至488.3Ma结束，累计时间跨度为32.7Ma，即理论时间跨度为32.7Ma。因此理论时间跨度与周期个数之比为11.2万年，因此与该比值最为接近的是米兰科维奇周期为10万年周期。

即巴探5井寒武系存在10万年的米兰科维奇周期。

另外一方面，在图6中小于333的尺度对应的模值不存在峰值，说明米兰科维奇周期中的4万年和2万年周期表现得并不明显。尽管识别4万年和2万年米兰科维奇周期存在困难，但是大体可以找到这两个周期所对应的尺度的大体位置，即4万年和2万年米兰科维奇周期对应的尺度因子分别为131和64。

因为4万年周期要比10万年周期个数多出2.5倍，2万年周期的个数要比10万年周期的个数多出5倍，因此分别利用10万年周期计算出4万年的理想的周期为730个(2.5×292＝730)，2万年的周期为1460个(5×292＝1460)。

而尺度因子为131和尺度64所对应的Morlet连续小波分别存在730个和1458个周期，其与利用4万年万年周期和2万年周期计算的理论周期个数，即420个周期和840个周期非常接近，因此尺度131和尺度64的选择非常合理，这也从另一方面验证了米兰科维奇周期为10万年周期的合理性。

图7-a为本发明实施例中巴探5井寒武系自然伽玛曲线图，图7-b为本发明实施例中巴探5井寒武系2万年米兰科维奇周期曲线图；图7-c为本发明实施例中巴探5井寒武系4万年米兰科维奇周期曲线图；图7-d为本发明实施例中巴探5井寒武系10万年米兰科维奇周期曲线图，其中，图7-d所示的10万年米兰科维奇周期曲线由292个周期组成，可用于地层年代的标定。

其次，依据自然伽玛识别巴探5井奥陶系米兰科维奇周期。

巴探5井奥陶系最下部的蓬莱坝组的底界深度为4214m，奥陶系最上部的鹰山组的顶界深度为3211m。在对奥陶系层段自然伽玛进行Morlet连续小波变换处理时，向上和向下分别延长到井深3196m和井深4235m，图8-a为本发明实施例中巴探5井奥陶系自然伽玛原始数据图；如图8-a所示，对这两个深度之间共计8314个自然伽玛数据进行连续Morlet连续小波变换。图8-b为本发明实施例中巴探5井寒武系Morlet连续小波变换时频色谱图，如图8-b所示，在进行连续小波变换时，取小波的最大长度为100，最小为0.1，滑动步长为0.1。在图8-b中，横坐标为自然伽玛的深度对应的数据；纵坐标为连续小波变换的不同尺度，其尺度从最小0.1到100；其颜色代表连续小波变换的系数大小，冷色为低值，暖色为高值，图8-b显示的色谱图是1000×8314二维数据矩阵直观图形显示，通过该色谱图可以粗略地看出自然伽玛主周期以及关键变化点。

图9为本发明实施例中巴探5井奥陶系Morlet连续小波变换的尺度因子-模值曲线图，为了更精确地找到自然伽玛的主周期，可以对不同尺度的Morlet连续小波变换数据进行分析，即是计算每一个尺度下的Morlet连续小波的模值，1000×8314小波系数矩阵中每一行即代表了每一尺度的小波系数，以尺度因子为横坐标，以计算得到的1000个模值为纵坐标，可以获得如图9所示的尺度因子-模值曲线图。

在图9中，尺度因子为387处的模值为一极大值，它所对应的Morlet连续小波反映有168个周期；数据处理的奥陶系大体的时间跨度为从奥陶纪蓬莱坝期开始，至鹰山期结束，与国际地质年表对比，历经的时间为从488.3Ma开始，至468.1Ma结束，累计时间跨度为20.2Ma，即理论时间跨度为20.2Ma。

理论时间跨度与周期个数之比为12万年，与该比值最接近的是米兰科维奇周期为10万年周期。

即巴探5井奥陶系存在10万年的米兰科维奇周期。

另外一方面，在图9中小于387尺度对应的模值不存在峰值，说明米兰科维奇周期中的4万年和2万年周期表现得并不明显。尽管识别4万年和2万年米兰科维奇周期存在困难，但是大体可以找到这两个周期所对应的尺度的大体位置，即4万年和2万年米兰科维奇周期对应的尺度因子分别为151和71。

因为4万年周期要比10万年周期个数多出2.5倍，2万年周期的个数要比10万年周期的个数多出5倍，因此分别利用10万年周期计算出4万年的理想的周期为420个(2.5×168＝420)，2万年的周期为840个(5×168＝840)。

而尺度因子为151和71所对应的Morlet连续小波分别存在423个和843个周期，这两个值分别与利用4万年万年周期和2万年周期计算的理论周期个数，即420个周期和840个周期非常接近，因此尺度151和尺度71的选择非常合理，这也从另一方面验证了米兰科维奇周期为10万年周期的合理性。

图10-a为本发明实施例中巴探5井奥陶系原始自然伽玛曲线图，图10-b为本发明实施例中巴探5井奥陶系2万年米兰科维奇周期曲线图；图10-c为本发明实施例中巴探5井奥陶系4万年米兰科维奇周期曲线图；图10-d为本发明实施例中巴探5井奥陶系10万年米兰科维奇周期曲线图，其中，图10-d所示的10万年米兰科维奇周期曲线由168个周期组成，可用于地层年代的标定。

以上是分别依据自然伽玛识别巴探5井寒武系以及奥陶系中米兰科维奇周期的过程，经过分析可知巴探5井寒武系以及奥陶系中米兰科维奇周期为10万年周期。

第三步，获得巴探5井奥陶系鹰山组残余层顶界年龄T_顶与底界底年龄T_底。

图11为巴探5井寒武系-奥陶系鹰山组岩性-深度剖面与Morlet连续小波变换的米兰科维奇周期曲线的对应关系示意图。如图11所示，沙依里克组的底界深度为5480m。

沙依里克组的底界深度为5480m，底界年龄为510Ma，以这个年龄值为起始点，依据米兰科维奇10万年周期为步长，对巴探5井鹰山组顶底界面进行年代标定。图12为巴探5井寒武系-奥陶系鹰山组依据米兰科维奇周期标定的岩性-时间剖面示意图，如图12所示，奥陶系鹰山组残余地层顶界年龄T_顶＝471.13Ma、底界年龄T_底＝479.8Ma。

第四步，获得巴探5井奥陶系鹰山组残余地层的平均沉积速度。

巴探5井奥陶系鹰山组顶界深度和底界深度分别为3210m和3760m，因此沉积厚度D＝550m，则沉积速度V按照下列定义式计算：

可得V＝63.4m/Ma。

第五步，获得巴探5井奥陶系鹰山组曾经拥有的时间T_剥。

理论上鹰山组保存完整时的顶面年龄为国际地质年表中大坪阶的顶面年龄468.1Ma，即T_标＝468.1Ma，而残余的奥陶系鹰山组的顶面年龄T_顶＝471.13Ma，因此，巴探5井奥陶系鹰山组曾经拥有的时间T_剥＝3.03Ma。

第六步，获得巴探5井奥陶系鹰山组不整合剥蚀量。

不整合剥蚀量Q＝T_剥×V＝3.03×63.4＝192.2m。

因此，巴探5井奥陶系鹰山组剥蚀量为192.2m。

综上所述，本发明提供的基于米兰科维奇旋回定年的剥蚀量恢复方法，是通过对识别米兰科维奇周期，并由此对测量区域的参数就进行标定，从而分别获得残余地层的平均沉积速率V与被剥蚀地层曾经拥有的时间T_剥这两个参数，并通过这两个参数获得不整合剥蚀量。因此本发明提供的方法是通过地层年代学研究找回被剥蚀地层曾经拥有的时间，即可以恢复地层剥蚀量，无需要花费高额的分析测试费用，其成本低、可操作性好，此外，本发明提供的方法不会受到岩石类型的限制，因此其应用没有局限性，应用更广泛。

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

Claims (8)

1.一种基于米兰科维奇旋回定年的剥蚀量恢复方法，包括以下步骤：

S10：确定测量区域内可进行区域对比且可标记绝对年龄的标志层；

S20：通过对测井曲线作Morlet连续小波变换以识别测量区域的米兰科维奇周期；

S30：以所述标志层为地层时间起算点，以所述米兰科维奇周期为步长进行标记，分别获得测量区域残余地层的顶界年龄T_顶和底界年龄T_底；

S40：测量残余地层的沉积厚度D，获得残余地层的平均沉积速率V，残余地层的平均沉积速率V根据下列定义式进行：

S50：确定与所述残余地层相对应的标准地层的顶界年龄T_标；获得被剥蚀地层曾经拥有的时间T_剥，被剥蚀地层曾经拥有的时间T_剥根据下列定义式进行：

T_剥＝T_顶-T_标；

S60：获得不整合剥蚀量Q，不整合剥蚀量Q根据下列定义式进行：

Q＝T_剥×V；

所述步骤S20中，对测井曲线作Morlet连续小波变换后，对不同尺度的Morlet连续小波变换数据进行分析，计算每一个尺度下的Morlet连续小波的模值，获得尺度因子-模值曲线，所述尺度因子-模值曲线的横坐标为尺度因子，纵坐标为模值；根据尺度因子-模值曲线获得模值的极大值处对应的Morlet连续小波的周期个数，根据国际地质年表获得测量区域的理论时间跨度；所述米兰科维奇周期为最接近理论时间跨度与周期个数之比的周期。

2.根据权利要求1所述的基于米兰科维奇旋回定年的剥蚀量恢复方法，其特征在于，步骤S50中，通过与相邻区域地层对比以及与国际地质年表对比，获得标准地层的顶界年龄T_标。

3.根据权利要求1或2所述的基于米兰科维奇旋回定年的剥蚀量恢复方法，其特征在于，所述沉积厚度D为残余地层的顶界深度与底界深度之差。

4.根据权利要求1或2所述的基于米兰科维奇旋回定年的剥蚀量恢复方法，其特征在于，步骤S10中，通过发育于所述标志层中的古生物资料与国际地质年代表对比，获得所述标志层的绝对年龄。

5.根据权利要求1或2所述的基于米兰科维奇旋回定年的剥蚀量恢复方法，其特征在于，步骤S20中测井曲线为自然伽玛曲线。

6.根据权利要求5所述的基于米兰科维奇旋回定年的剥蚀量恢复方法，其特征在于，步骤30中，通过对所述米兰科维奇周期曲线进行标定，获得所述残余地层的顶界年龄T_顶和底界年龄T_底，进行标定时使所述标志层的底界与所述米兰科维奇周期曲线的波谷相对应。

7.根据权利要求1或2所述的基于米兰科维奇旋回定年的剥蚀量恢复方法，其特征在于，步骤S30中，所述米兰科维奇周期为10万年周期、4万年周期或2万年周期。

8.根据权利要求1或2所述的基于米兰科维奇旋回定年的剥蚀量恢复方法，其特征在于，所述标志层为最大海泛期形成的凝缩层或火山喷发形成的广泛分布的火山凝灰岩层。

Images