一种利用旋回分析法恢复不整合剥蚀厚度的方法与流程

本发明属于地质勘探技术领域，尤其涉及一种利用旋回分析法恢复不整合剥蚀厚度的方法。

背景技术：

恢复不整合剥蚀厚度是盆地构造演化分析和油气成藏研究的重要内容，其结果的准确性对烃源岩热演化、储层发育与改造、油气运聚与保存条件的评价具有重要的影响。

目前恢复不整合剥蚀厚度的方法较多，包括泥岩声波时差法、镜质体反射率法、包裹体均一温度法，以及基于地震解释的地震趋势分析法、邻层厚度比值法、参考层厚度变化率法等。然而，实践发现这些方法仍存在一些不足：(1)由于上述现有方法均是通过建立一种数学关系来实现对剥蚀厚度的定量恢复，而这种数学关系的建立往往具有不唯一性，从而导致剥蚀厚度值具有多解性，且误差较大；(2)上述现有方法都有各自的适用条件和局限性，例如：镜质体反射率法需要较多的测试数据，且受热事件影响十分明显；基于地震解释的地震趋势分析法则需建立在研究区剥蚀范围较小的基础上，对于剥蚀范围覆盖全区的地区不适用；(3)由于一个不整合可能是一个层系的剥蚀产物，也可能是多个层系的剥蚀产物，当剥蚀程度较大、多个层系剥蚀形成一个不整合时，利用上述现有方法仅能恢复一套不整合的剥蚀厚度，即多层系剥蚀厚度之和，无法对不同层系的剥蚀厚度进行恢复。

因此，如何提供一种准确度高、适用范围广且能够恢复不同剥蚀层系剥蚀厚度的方法是当前急需解决的一项技术难题。

技术实现要素：

本发明针对上述的现有恢复不整合剥蚀厚度的方法误差大、适用范围具有局限性、且无法对不同层系的剥蚀厚度进行恢复的技术问题，提出一种利用旋回分析法恢复不整合剥蚀厚度的方法，能够有效还原研究区在整个地质历史时期的地层沉积-剥蚀状态，可恢复不同剥蚀层系的剥蚀厚度，且该方法准确度高、适用范围广。

为了达到上述目的，本发明采用的技术方案为：

一种利用旋回分析法恢复不整合剥蚀厚度的方法，包括以下步骤：

沉积速率曲线绘制：选取研究区，结合研究区的地质资料，通过去压实校正计算压实系数和压实前地层的原始厚度，根据地层年代表确定每个地层沉积的起始时间和结束时间，计算各沉积期的平均地层沉积速率，绘制沉积速率曲线；

沉积速率曲线的频谱分析：根据沉积速率曲线的旋回规律确定缺失曲线段；根据缺失曲线段前后沉积速率曲线的周期值，推算缺失曲线段对应的周期、缺失曲线段波峰和波谷对应的横坐标；对沉积速率曲线进行分频处理得到对应低频子波和高频子波，根据缺失曲线段前后对应的低频子波和高频子波的波峰和波谷的纵坐标,推算缺失曲线段波峰和波谷对应的纵坐标；

计算不整合剥蚀厚度：根据频谱分析结果恢复完整的沉积速率曲线，计算不整合剥蚀厚度，所述不整合剥蚀厚度的计算公式如下：

其中，H_e为不整合剥蚀厚度，单位：m；F(t)为完整沉积速率曲线函数；H_oc为在沉积时间段内对F(t)积分得到的地层原始厚度；H_r为现今残余厚度，单位：m；C为总压实系数；t₁为剥蚀地层沉积开始时间，单位：Ma；t₂为剥蚀地层沉积结束时间，单位：Ma；

作为优选，在沉积速率曲线绘制步骤中，所述压实前地层的原始厚度的计算公式为：

H₀＝H_s/C_s+H_m/C_m

其中，H₀为地层压实前的原始厚度，单位：m；H_s为砂岩厚度，单位：m；C_s为砂岩压实系数，利用去压实校正方法计算获得；H_m为泥岩厚度，单位：m；C_m为泥岩压实系数，利用去压实校正方法计算获得；

在计算不整合剥蚀厚度时，所述总压实系数的计算公式为：

其中，C为总压实系数；H_s为砂岩厚度，单位：m；C_s为砂岩压实系数，利用去压实校正方法计算获得；H_m为泥岩厚度，单位：m；C_m为泥岩压实系数，利用去压实校正方法计算获得。

作为优选，所述地质资料包括地层分层数据、地层年代表，以及录井资料或孔隙度资料。

作为优选，所述去压实校正方法为经验公式法或孔隙度-岩性方法；当所述地质资料包括录井资料时，所述去压实校正方法采用经验公式法；当所述地质资料包括孔隙度资料时，所述去压实校正方法采用孔隙度-岩性方法。

作为优选，在计算不整合剥蚀厚度步骤后，还包括利用地层分层数据或现有恢复不整合剥蚀厚度的方法验证所得不整合剥蚀厚度准确性的步骤。

作为优选，利用地层分层数据进行准确性验证的具体步骤为：利用恢复的完整沉积速率曲线，通过对未发生剥蚀的地层对应的沉积速率曲线段进行积分得到未剥蚀地层的原始厚度，通过压实计算得到未剥蚀地层的压实后地层厚度，与地层分层数据中的现今地层厚度进行对比，验证准确性。

作为优选，所述现有恢复不整合剥蚀厚度的方法选自泥岩声波时差法、镜质体反射率法、磷灰石裂变径迹法、物质平衡法、包裹体均一温度法、沉积速率法、砂岩孔隙度法、波动分析法、地震趋势分析法、邻层厚度比值法和参考层厚度变化率法中的任一种。

与现有技术相比，本发明的优点和积极效果在于：

1、本发明能够有效还原研究区在整个地质历史时期的地层沉积-剥蚀状态，恢复任何一个时间范围内的地层沉积原始厚度，从而可恢复不同剥蚀层系的剥蚀厚度，可以为构造演化、油气在不同历史时期的运移方向等研究奠定重要基础，有利于降低油气勘探风险；

2、本发明仅利用地层分层数据、地层年代表、录井资料和孔隙度资料等基础数据，对数据的要求简单且需要的数据量较少，排除了因实验数据数量少或机器误差、选样差异等带来的误差，提高了准确性；

3、本发明对于研究区的剥蚀范围或地层压力等均无特殊限制与要求，适用范围更为广泛；

4、本发明不依赖于实验，不需要特殊实验设备的支撑，操作简单且成本更低。

附图说明

图1为本发明实施例所提供的利用旋回分析法恢复不整合剥蚀厚度的方法的流程图；

图2为本发明实施例1所提供的B202井恢复前的沉积速率折线图；

图3为本发明实施例1所提供的B202井恢复前的沉积速率曲线图；

图4为本发明实施例1所提供的B202井恢复前沉积速率曲线对应的低频子波曲线图；

图5为本发明实施例1所提供的B202井恢复前沉积速率曲线对应的高频子波曲线图；

图6为本发明实施例1所提供的B202井完整的沉积速率折线图；

图7为本发明实施例1所提供的B202井完整的沉积速率曲线图；

图8为本发明实施例1所提供的B202井完整沉积速率曲线对应的低频子波曲线图；

图9为本发明实施例1所提供的B202井完整沉积速率曲线对应的高频子波曲线图。

具体实施方式

下面将对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

本发明实施例提供了一种利用旋回分析法恢复不整合剥蚀厚度的方法，其流程图如图1所示，包括以下步骤：

S1沉积速率曲线绘制：选取研究区，结合研究区的地质资料，通过去压实校正计算压实系数和压实前地层的原始厚度，根据地层年代表确定每个地层沉积的起始时间和结束时间，计算各沉积期的平均地层沉积速率，绘制沉积速率曲线。

在本步骤中，需要说明的是，由于在后沉积的地层对在先沉积的地层具有压实作用，通过去压实校正可获得压实前地层的原始厚度，采用此原始厚度计算得到的平均地层沉积速率更接近于该沉积期真实的地层沉积速率。所述平均地层沉积速率的计算公式为：

V＝H₀/(t₄-t₃) (1)

式(1)中，V为各沉积期的平均地层沉积速率，单位：m/Ma；H₀为地层压实前的原始厚度，单位：m；t₃为各沉积期地层沉积开始时间，单位：Ma；t₄为各沉积期地层沉积结束时间，单位：Ma。

此外，还需要说明的是，沉积速率曲线的绘制可利用Matlab软件进行，以时间为横坐标、平均地层沉积速率为纵坐标绘制沉积速率折线图，进一步通过三次样条插值方法得到沉积速率曲线图。

S2沉积速率曲线的频谱分析：根据沉积速率曲线的旋回规律确定缺失曲线段；根据缺失曲线段前后沉积速率曲线的周期值，推算缺失曲线段对应的周期、缺失曲线段波峰和波谷对应的横坐标；对沉积速率曲线进行分频处理，得到沉积速率曲线对应的低频子波和高频子波，根据低频子波和高频子波推算缺失曲线段波峰和波谷对应的纵坐标。

在本步骤中，需要说明的是，由于地层的沉积具有一定的旋回性特征，具体表现为：地层的沉积速率呈近似正弦波动，其周期、波峰和波谷数值均呈现一定的渐变规律。然而，地层的剥蚀会造成沉积速率数据缺失，进而影响沉积速率曲线的这种旋回规律，因此，可确定在沉积速率曲线中与前后旋回规律不一致的一段曲线段即为缺失曲线段。此外，结合区域地质资料及地球物理响应(即测井曲线突变、地震不整合识别等)，可验证缺失曲线段对应时间区间的准确性。

同时，沉积速率曲线的周期、波峰和波谷数值也呈现一定的线性渐变规律，即沉积速率曲线的周期、波峰和波谷数值呈递减或递增规律。而且，缺失曲线段对应的时间必须介于前后曲线之间，即时间约束原则。因此，利用这种线性渐变规律和时间约束原则，根据缺失曲线段前后沉积速率曲线的周期值，可推算出缺失曲线段缺失的周期数，进一步可推算出缺失曲线段对应的周期、缺失曲线段波峰和波谷对应的横坐标。

此外，虽然缺失曲线段波峰和波谷对应的纵坐标(即沉积速率)数值也可同理进行推算，但利用分解重构思想对沉积速率曲线进行分频处理，得到的高频子波的波动变化可反映沉积的旋回规律，低频子波的波动变化可反映沉积的大致趋势，由此推算得到的数据更为精确性，有利于提高频谱分析的准确性。具体的，对沉积速率曲线进行分频处理得到低频子波和高频子波，分别读取两个子波对应的波峰值，二者叠加即可得到总的波峰数值Yc，同理可得波谷数值Yt。通过此方法分别计算出缺失曲线段前后沉积速率曲线对应的波峰和波谷数值，再利用线性渐变规律，即可推算出缺失曲线段波峰和波谷对应的纵坐标。

还需要说明的是，分频处理可以在Matlab软件中完成，具体步骤为：通过(i-N+1)次N点取平均方法获得沉积速率曲线的低频子波，将沉积速率曲线减去所述低频子波得到高频子波。其中，N点取平均方法具体为：对于任意一维数据x，选取固定时窗对于数据内各个采样点x_i，逐步滑动时窗求取平均值，求取过程可表示为：

式(2)中，i为数据点个数；t_i为第i个数据点所对应的时间点；N为数据采样点个数；j为循环变量；Δt为时间采样间隔。

S3计算不整合剥蚀厚度：根据频谱分析结果恢复完整的沉积速率曲线，计算不整合剥蚀厚度，所述不整合剥蚀厚度的计算公式如下：

在本步骤中，根据频谱分析获得的缺失曲线段对应的周期、缺失曲线段波峰和波谷对应的横坐标和纵坐标，即可恢复完整的沉积速率曲线，其中波峰处为沉积速率最大的时期，波谷处为沉积速率最小的时期，该曲线可真实反映沉积速率的旋回变化和沉积演化过程。进一步，通过对完整沉积速率曲线积分，即可得到未受压实作用的地层原始厚度，利用之前计算获得的总压实系数对该厚度进行压实计算，即可得到地层的原始厚度，而剥蚀厚度即为地层的原始厚度与现今残余厚度的差值。因此，采用式(3)即可计算不整合剥蚀厚度。需要说明的是，式(3)中的参数t₁和t₂是针对某一具体剥蚀地层的沉积开始和结束时间，因而本实施例提供的利用旋回分析法恢复不整合剥蚀厚度的方法可以针对不同剥蚀地层进行单独计算，进而可恢复不同剥蚀层系的剥蚀厚度。

在一优选实施例中，在沉积速率曲线绘制步骤中，所述压实前地层的原始厚度的计算公式为：

H₀＝H_s/C_s+H_m/C_m (4)

在计算不整合剥蚀厚度时，所述总压实系数的计算公式为：

在本实施例中，进一步限定了压实前地层的原始厚度的计算公式和总压实系数的计算公式，需要说明的是，由于地层中包含砂岩和/或泥岩，而砂岩与泥岩的压实系数不同，因此，需要分别对砂岩和泥岩进行去压实校正。

在一优选实施例中，所述地质资料包括地层分层数据、地层年代表，以及录井资料或孔隙度资料。本实施例进一步限定了最优的地质资料组合，可以理解的是，本实施例并不局限于上述所列举的地质资料种类，还可以为本领域技术人员所熟知的其它常用地质资料。

在一优选实施例中，所述去压实校正方法为经验公式法或孔隙度-岩性方法；当所述地质资料包括录井资料时，所述去压实校正方法采用经验公式法；当所述地质资料包括孔隙度资料时，所述去压实校正方法采用孔隙度-岩性方法。本实施例进一步限定了去压实校正方法，需要说明的是，本领域技术人员可根据地质资料种类和地层岩性特征具体选择合适的去压实校正方法。而且，本实施例并不局限于上述所列举的去压实校正方法，还可以为本领域技术人员所熟知的其它常用去压实校正方法。

此外，还需要说明的是，利用所述经验公式法计算压实系数的公式为：

C'＝b-alnH (6)

式(6)中，当计算砂岩压实系数时，C’即为C_s；当计算泥岩压实系数时，C’即为C_m；a、b为常数，对于砂岩取a＝0.08、b＝1.37，对于泥岩取a＝0.1、b＝1.46；H为埋藏深度，单位：m。

利用所述孔隙度-岩性方法，即利用不同岩性孔隙度与深度的关系计算压实系数的公式为：

式(7)中，当计算砂岩压实系数时，C”即为C_s；当计算泥岩压实系数时，C”即为C_m；H为埋藏深度，单位：m；为深度为H处的孔隙度；为地表初始孔隙度。

在一优选实施例中，在计算不整合剥蚀厚度步骤后，还包括利用地层分层数据或现有恢复不整合剥蚀厚度的方法验证所得不整合剥蚀厚度准确性的步骤，通过该验证步骤可以保证所得不整合剥蚀厚度的准确性。

在一优选实施例中，利用地层分层数据进行准确性验证的具体步骤为：利用恢复的完整沉积速率曲线，通过对未发生剥蚀的地层对应的沉积速率曲线段进行积分得到未剥蚀地层的原始厚度，通过压实计算得到未剥蚀地层的压实后地层厚度，与地层分层数据中的现今地层厚度进行对比，验证准确性。在本实施例中，通过未剥蚀地层的压实后地层厚度与现今地层厚度的误差大小，即可验证所得不整合剥蚀厚度的准确性，这种验证方法不需要其他方法和数据的支撑，操作简单、工作量小且成本低。

在一优选实施例中，所述现有恢复不整合剥蚀厚度的方法选自泥岩声波时差法、镜质体反射率法、磷灰石裂变径迹法、物质平衡法、包裹体均一温度法、沉积速率法、砂岩孔隙度法、波动分析法、地震趋势分析法、邻层厚度比值法和参考层厚度变化率法中的任一种。在本实施例中，需要说明的是，本领域技术人员可根据地质资料和现有恢复不整合剥蚀厚度方法的适用条件，具体选择合适的现有恢复不整合剥蚀厚度方法进行验证。而且，本实施例采用的验证方法并不局限于上述所列举的恢复不整合剥蚀厚度方法，还可以为本领域技术人员所熟知的其它常用方法。

为了更清楚详细地介绍本发明实施例所提供的利用旋回分析法恢复不整合剥蚀厚度的方法，下面将结合具体实施例进行描述。

实施例1

在松辽盆地长岭断陷龙凤山地区，选取一口典型井B202，该井位于龙凤山地区的东南缓坡带，钻井钻到的白垩纪地层从上至下依次为明水组、四方台组、嫩江组、姚家组、青山口组、泉头组、登娄库组、营城组、沙河子组(未钻穿)。

对B202井的不整合剥蚀厚度进行恢复，步骤如下：

S1沉积速率曲线绘制：

结合B202井的地层分层数据、地层年代表和录井资料，通过经验公式法进行去压实校正，利用式(6)计算压实系数，利用式(4)计算压实前地层的原始厚度，根据地层年代表确定每个地层沉积的起始时间和结束时间，利用式(1)计算得到该井各地层沉积时期的平均地层沉积速率，数据如表1所示。

表1 B202井各地层沉积时期的平均地层沉积速率统计

表1(续)B202井各地层沉积时期的平均地层沉积速率统计

根据表1的计算结果，利用Matlab软件以时间为横坐标、平均地层沉积速率为纵坐标，绘制B202井的沉积速率折线图，如图2所示。通过三次样条插值方法得到沉积速率曲线图，如图3所示。

S2沉积速率曲线的频谱分析：

根据图3中沉积速率曲线的旋回规律可以明显看出，在第7个点和第8个点之间缺失数据，确定缺失曲线段对应的时间区间为109.835Ma-118Ma。根据图2和图3的频谱分析可知，缺失曲线段之前(即94.25Ma-109.835Ma时间内)的曲线为2个周期，其平均周期T₁＝7.7925Ma；缺失曲线段之后(即118Ma-135.1429Ma时间内)的曲线为2.75个周期，其平均周期T₂＝6.234Ma。由此可得，沉积速率曲线的周期随着年代从新到老呈递减规律。另外，根据图3中沉积速率曲线的旋回规律可见，缺失曲线段应为1.25个周期，即缺失一个波峰和一个波谷，因此，缺失曲线段的平均周期为6.532Ma，根据线性渐变和时间约束原则，缺失曲线段波峰对应的横坐标(即时间值)为113.0605Ma，波谷对应的横坐标为116.286Ma。

对图3中的沉积速率曲线进行分频处理，利用式(2)提取低频子波，如图4所示，将沉积速率曲线减去所述低频子波得到高频子波，如图5所示。从图4和图5中分别读取缺失曲线段之前低频子波和高频子波对应的波峰值，二者叠加得到缺失曲线段之前的波峰值Yc₁＝81.79，同理得到缺失曲线段之前的波谷值Yt₁＝43.87；同理，得到缺失曲线段之后的波峰值Yc₂＝148.68、波谷值Yt₂＝39.01。利用沉积速率曲线的波峰和波谷值具有线性渐变规律，可推算出缺失曲线段波峰对应的纵坐标(即沉积速率)为115.235m/Ma，波谷对应的纵坐标为41.44m/Ma。

S3计算不整合剥蚀厚度：

根据以上频谱分析结果，将(113.0605，115.235)和(116.286，41.44)两个数据点补入到图2的数据中，绘制的完整沉积速率折线图如图6所示，对应的完整沉积速率曲线图如图7所示，对应的低频子波和高频子波曲线图分别如图8、图9所示。B202井发生剥蚀的地层为登娄库组(112Ma-124Ma)，利用式(5)计算得到登娄库组对应的总压实系数C为0.7207，将C＝0.7207、t₁＝112Ma、t₂＝124Ma、H_r＝318.8m代入式(3)中，计算得到登娄库组的不整合剥蚀厚度H_e为379.69m。

S4方法验证：

利用地层分层数据进行验证，利用恢复的完整沉积速率曲线，通过对未发生剥蚀的地层对应的沉积速率曲线段进行积分得到未剥蚀地层的原始厚度，通过压实计算得到未剥蚀地层的压实后地层厚度，数据如表2所示。

表2利用B202井地层分层数据验证结果统计表

将未剥蚀地层的压实后地层厚度与从地层分层数据中获取的现今地层厚度进行对比，如表2所示，除青山口组地层外，其余地层的误差率均小于10％。根据松辽盆地构造演化史的相关文献可知，青山口组地层发生过地层剥蚀，该地层的误差正是地层剥蚀导致的。由此可见，本实施例提供的利用旋回分析法恢复不整合剥蚀厚度的方法是可靠的、有效的。

实施例2

与实施例1的不同是，本实施例采用泥岩声波时差法验证所得不整合剥蚀厚度的准确性，除方法验证步骤外，其余步骤及其所得数据结果同实施例1。由于泥岩声波时差法为现有常用的恢复不整合剥蚀厚度的方法，此处不再赘述。本实施例计算所得的B202井登娄库组的不整合剥蚀厚度为379.69m，采用泥岩声波时差法计算结果为408.19m，二者误差率为6.98％，小于10％，由此可见，本实施例提供的利用旋回分析法恢复不整合剥蚀厚度的方法是可靠的、有效的。

实施例3

对松辽盆地长岭断陷龙凤山地区B204井的不整合剥蚀厚度进行恢复，除方法验证步骤外，其余步骤同实施例1，计算得到B204井的不整合剥蚀厚度为418.09m，本实施例采用包裹体均一温度法验证所得不整合剥蚀厚度的准确性。由于包裹体均一温度法为现有常用的恢复不整合剥蚀厚度的方法，此处不再赘述。采用第一类包裹体均一温度法计算的结果为388.58m，误差率为7.60％；采用第二类包裹体均一温度法计算的结果为446.15m，误差率为6.29％。由此可见，本实施例提供的利用旋回分析法恢复不整合剥蚀厚度的方法是可靠的、有效的。