本发明涉及一种定量识别米兰科维奇周期的方法，属于地质地层技术领域。
背景技术
地球轨道几何形态的周期性变化在控制气候变化的同时对沉积物充填也会产生重要的影响，milankovitch对轨道要素(偏心率、斜率和岁差)与气候间联系的定量描述是旋回地层学发展的重要里程碑，偏心率、斜率和岁差对应的周期被统称为米兰科维奇周期。各级米兰科维奇周期间存在相对稳定的比例关系，在稳定沉积环境相对整合的地层中，这一比例应与各级地层旋回厚度的比值一致。这也是米兰科维奇周期研究的基础。稳定沉积地层的岩心、露头及与气候变化相关联的替代性指标均可用于米兰科维奇周期的研究。
常用的米兰科维奇周期识别方法是，通过对替代指标数据进行频谱分析，将数据从深度域转换为频率域，记录下频谱图上可信度之上峰值对应的频率，通过将频率之间的比例关系与理论轨道周期的比例关系进行比对，找到与理论周期比例接近的一组频率组合，该频率组合对应的厚度即与米兰科维奇周期有关。因此，确定与理论轨道周期比例最接近的频率组合，是识别米兰科维奇周期的关键步骤。目前对此关键一环，主要是根据理论周期比例与频率比例近似比对，认为二者接近即可，属于定性分析判断，没有定量确定二者之间的接近程度，也就无法说明是否还有更加接近的比例组合，这样会造成识别出的米兰科维奇周期存在偏差。当分析的目标地层厚度较大时，频谱分析得到的峰值频率较多，相应的频率之间的比例关系组合也较多，通过定性分析找到最有匹配就显得费时费力，异常困难。
技术实现要素：
针对上述问题，本发明的目的是提供一种定量识别米兰科维奇周期的方法，利用该方法识别出的米兰科维奇周期更为准确。
为实现上述目的，本发明采用以下技术方案：一种定量识别米兰科维奇周期的方法，包括以下步骤：
1)选择自然伽马曲线作为米兰科维奇周期识别的替代性指标；
2)对自然伽马数据进行预处理；
3)对理论轨道周期值进行识别，获得理论周期的降序排列集合m；
4)对地层中的米兰科维奇周期进行识别，具体过程如下：
①对目标地层中的自然伽马数据进行整理；
②对选定地层的自然伽马数据进行频谱分析，将数据从深度域转换为频率域，获得数据频谱分布图后，分别记录90％置信度以上的峰值对应的频率f，所有频率f记作升序集合f＝{x1,x2,...,xn}；由于周期与频率呈倒数关系，为了便于比对匹配，再对集合f中的频率取倒数，得到新的频率值倒数降序集合t＝{t1,t2,...,tn}，同时记录下与该频率对应的频谱振幅y，记作降序集合y＝{y1,y2,...,yn}，以反映频率的显著程度；
③实施匹配：从理论周期集合m中随机挑选n个数据，有
从频率值倒数集合t中随机挑选n个数据，有
引入参数σ来定量表征两个比例数据组合l和s的差异，定义
σ(i,j)＝[pi(1)-qj(1)]2+[pi(2)-qj(2)]2+...+[pi(n)-qj(n)]2
其中i为理论周期数组中某行数据的行标，
由于米兰科维奇周期在地质历史时期的显著性，每个频率的显著程度也要考虑，因此引入另一个参数λ表示频率的显著程度：
λ＝z1+z2+...+zn，其中z1，z2，…，zn分别是频率倒数组合s中第j行元素对应的振幅值；
根据σ大小按降序排列，选出σ最小的5种组合，考虑每种频率组合中各频率对应的振幅，按λ降序排列，找到λ最小值对应的理论周期比值组合pi和频率倒数比值组合qj，与使λ最小的理论周期比值组合pi对应的理论周期组合(数组k中的第i行)即为该地层中保存的米兰科维奇周期，使λ最小的频率倒数比值组合qj除以采样间隔即为各米兰科维奇周期对应的旋回厚度。
所述步骤2)中对自然伽马数据进行预处理过程如下：先以井径曲线为参考对自然伽马数据进行环境校正，再对其进行标准化使研究区内各井相同测井曲线具有相同的刻度水平。
所述步骤4)的①中，对目标地层中的自然伽马数据进行整理的过程为：跟进井位分层数据确定分析时窗，将自然伽马数据整理成各自时窗的单个文件，每个文件未两列数据，第一列为等采样间隔的深度数据，第二列是与深度数据对应的自然伽马数据。
本发明由于采取以上技术方案，其具有以下优点：本发明通过频谱分析确定理论周期值和实际地层中的峰值频率值后，引入参数σ表征二者比例之间的差异大小，根据σ值优选出与理论周期比例接近的c(一般小于等于5)个实际地层频率组合，再根据各频率组合的振幅值之和对其降幂排序，从中选出振幅值最大的组合作为与米兰科维奇周期对应的地层旋回，即该旋回厚度受米兰科维奇周期控制。本发明提出了一种定量的识别米兰科维奇周期的方法，并且考虑了振幅的影响，因此所识别的结果更为准确。
附图说明
图1是对某裂谷盆地中-晚中新统地层计算偏心率、斜率与岁差变化的理论值后，进行频谱分析所得到的结果；其中，图(a)是频谱分析识别出的偏心率周期，图(b)是频谱分析识别出的斜率周期，图(c)是频谱分析识别出的岁差周期；
图2是某裂谷盆地中-晚中新统地层gr频谱分析结果曲线。
具体实施方式
下面结合附图和实施例对本发明进行详细的描述。
本发明提出了一种定量识别米兰科维奇周期的方法，包括以下步骤：
1)选择自然伽马曲线作为米兰科维奇周期识别的替代性指标。
稳定沉积地层的岩心、露头及与气候变化相关联的替代性指标均可用于米兰科维奇旋回的研究。替代性指标包括符合采样密度要求的地球化学分析数据(如δ18o、δ13c、δ88sr、caco3含量、磷/钛等)和能反映气候变化的地球物理参数(如磁化率、自然伽马曲线、色率等)。自然伽马测井是在井内测量岩层中的放射性元素原子核衰变过程中放射出的伽马射线的强度。粘土物质和有机质对放射性物质的吸附能力较强，而且泥质物质颗粒小，沉积缓慢，使得放射性元素有足够的时间从溶液中分离出来。因此，自然伽马曲线能够反映沉积物中的泥质和有机质含量变化，进而反映古环境和古气候的变化，能够作为米兰科维奇周期研究的替代性指标。另外，由于几乎所有油田的测井系列中均包含自然伽马曲线，使得该方法也具有易于推广的普遍性。
2)对自然伽马数据进行预处理，具体地，先以井径曲线为参考对自然伽马数据进行环境校正，再对其进行标准化使研究区内各井相同测井曲线具有相同的刻度水平。
在对米兰科维奇周期进行识别之前，为了消除环境因素造成的偏差，在进行旋回地层分析前，以井径曲线为参考曲线对gr数据(自然伽马数据)进行环境校正。另外，本发明使用的是自然伽马绝对值，测井数据除了受到环境因素干扰产生偏差外，还会因测井时间和次数不同、测井仪器类型和刻度不同造成系统偏差。为了增强多井间测井数据的可对比性，保证所用数据的可靠性，在对测井曲线进行环境校正后，再对其进行标准化，使研究区内各井相同测井曲线具有相同的刻度水平。一般情况下，自然伽马数据都是等间距采样的，如果存在非等间距的数据，要对其进行重采样，确保数据采样间隔相等。数据预处理工作均通过forward测井解释软件完成。
3)对理论轨道周期值进行识别，具体识别方法可采用现有的轨道要素计算方法。
地球轨道参数的变化在对地球施加影响的同时，由于受到地月系统围绕太阳公转和地月相互作用的影响也在发生改变。地月引力导致的地表潮汐摩擦力和轨道重力会影响地球自转速率和地球形状，使地球的转动速度变慢，岁差和地轴斜率周期变长。目前的轨道要素计算方法中，具有代表性的有berger、laskar等提出的方案，其中laskar(2004)提出的解决方案la(2004)则综合考虑了上述影响因素。以某裂谷盆地中-晚中新统地层为例，采用la(2004)计算了该层位沉积时间段16ma(百万年)至5.3ma偏心率、斜率与岁差变化的理论值，采样间隔为1ka(千年)。上述对轨道周期理论数据分别进行频谱分析后，识别出3个偏心率周期405ka、125ka、95ka，2个斜率周期53ka、40ka，3个岁差周期24ka、22ka、19ka，共8个理论周期，记作降序集合m＝{405，125，95，53，40，24，22，19}，如图1所示。这些理论轨道周期之间存在的稳定比例关系，可作为确定天文周期的基准。
4)对地层中的米兰科维奇周期进行识别，具体过程如下：
①对目标地层中的自然伽马数据进行整理，跟进井位分层数据确定分析时窗，将自然伽马数据整理成各自时窗的单个文件，每个文件未两列数据，第一列为等采样间隔的深度数据，第二列是与深度数据对应的自然伽马数据。
②对选定地层的自然伽马数据进行频谱分析，将数据从深度域转换为频率域，获得数据频谱分布图后，分别记录90％置信度以上的峰值对应的频率f，所有频率f记作升序集合f＝{x1,x2,...,xn}。由于周期与频率呈倒数关系，为了便于比对匹配，再对集合f中的频率取倒数，得到新的频率值倒数降序集合t＝{t1,t2,...,tn}。同时记录下与该频率对应的频谱振幅y，记作集合y＝{y1,y2,...,yn}，以反映频率的显著程度。
根据比例关系识别米兰科维奇周期应遵循以下几个原则：①频谱频率间的比例关系与理论周期比例的差异应尽可能的小，一般控制在5％以内；②一段地层中识别出的天文周期个数应尽量多，一般不小于3个；③识别出的天文周期与对应的旋回厚度算出的沉积速率应符合研究区沉积规律。
③实施匹配，具体地：
从理论周期集合m中随机挑选n个数据，有
从频率值倒数集合t中随机挑选n个数据，有
引入参数σ来定量表征两个比例数据组合l和s的差异，定义
σ(i,j)＝[pi(1)-qj(1)]2+[pi(2)-qj(2)]2+...+[pi(n)-qj(n)]2
其中i为理论周期数组中某行数据的行标，
由于米兰科维奇周期在地质历史时期的显著性，每个频率的显著程度也应被考虑，本发明引入参数λ表示频率的显著程度：
λ＝z1+z2+...+zn，其中z1，z2，…，zn分别是频率倒数组合s中第j行元素对应的振幅值。
根据σ大小按降序排列，选出σ最小的5种组合，考虑每种频率组合中各频率对应的振幅，按λ降序排列，找到λ最小值对应的理论周期比值组合pi和频率倒数比值组合qj。与使λ最小的理论周期比值组合pi对应的理论周期组合(数组k中的第i行)即为该地层中保存的米兰科维奇周期，使λ最小的频率倒数比值组合qj除以采样间隔即为各米兰科维奇周期对应的旋回厚度。
根据识别出的米兰科维奇周期及其对应的旋回厚度，可以估算出该地层的沉积速率，该结果必须符合区域沉积变化规律，不能偏差过大，一般情况下差异应小于10％。通常与自然伽马曲线对应的地层深度为测量深度，当地层倾角和井斜角较大时，需要将测量地层厚度校正到真实地层厚度，以获得更准确的沉积速率数据。只有经过校正后的真实沉积速率与区域沉积规律匹配，才能最终认为识别出的旋回是受米兰科维奇周期控制的地层旋回。如果差异较大，应重复步骤4)中的③，直到找到符合区域沉积规律的米兰科维奇周期组合及其对应的旋回厚度。
下面以一具体实施例来说明本发明的效果：
应用本发明的专利思路与步骤，对某裂谷盆地中-晚中新统地层开展了米兰科维奇周期研究，进行应用验证。研究区属于陆相扇三角洲前缘沉积，地层主要为砂泥岩互层，受气候因素影响显著。自然伽马数据采样间隔为0.15米，数值分布范围30-200api，地层倾角7度，井斜角22度。
对该地层的gr数据预处理后进行频谱分析，记录下频谱图上的17个峰值对应频率，记作升序集合f，然后对集合f中的频率取倒数，得到新的频率值倒数降序集合t，并记录下各频率对应的振幅值，记作集合y。频率的倒数为对应的旋回测深厚度md，真实旋回厚度tvt等于md乘以地层倾角和井斜角之和的余弦，校正后的真实地层厚度即米兰科维奇周期对应的真实厚度，可以估算比较精确的沉积速率。
以n＝5为例，即从集合f中随机挑选5个频率进行归一化，从理论周期中随机挑选5个周期并归一化处理，采用基于matlab编制的自动匹配程序，对二者进行匹配分析。差异σ最小的5种组合见表1，再根据振幅值y进行优选，获得n＝5时的最优选择为表1中的组合3。因此，该地层中保存的5个米兰科维奇周期分别是125ka、95ka、53ka、40ka和24ka，其对应的旋回厚度分别是20.99m、16.28m、8.67m、6.70m、4.07m，如图2所示。同样的步骤可以对n取其他值的进行分析，其对应的最终组合，再将n不同时得到的组合按照σ排序，以获得最终最为可靠的米兰科维奇周期及其对应的旋回厚度数据。本次研究发现，n＝5时得到的周期组合是该裂谷盆地中-晚中新统地层中的最优米氏周期组合，说明该地层中保留了125ka、95ka、53ka、40ka和24ka五种米兰科维奇周期。据此计算沉积速率在0.17m/ka左右，该结果与调研的区域沉积规律一致，因此认为本次分析结果可靠，准确性高。由此可见，本发明可以自动快速识别出地层中保留的米兰科维奇周期，并计算出相应的沉积速率，对于加深陆相地层沉积控制因素及演化规律具有重要意义。该方法简便快捷，具有极高的推广应用价值。
表1n＝5时自动匹配结果
上述各实施例仅用于说明本发明，其中方法的实施步骤等都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。