基于米兰科维奇旋回的古水深定量恢复方法和存储介质与流程

1.本发明涉及地球化学和地球物理中古水深恢复技术领域，具体涉及一种基于米兰科维奇旋回的古水深定量恢复方法和存储介质。


背景技术：

2.现有技术中，在文献(杨万芹,朱德顺,银燕,丁桔红,2015.古水深的地球化学恢复方法及在层序地层划分中的应用.地质论评,61(s1):756
‑
757.)中，根据fe/co、fe/mn和th/u三个元素比值综合，在结合地区平均古水深进行古水深曲线恢复。此方法主要建立对古水深变化反映灵敏的地球化学指标，进行多个指标综合，进行古水深的量化恢复，并应用恢复的古水深曲线。但是此方法无法对无取心或无微量元素测试的井段进行水深的恢复。
3.在另一篇专利文件(申请公布号cn 106019401 a)“一种古水深的定量恢复方法以及装置”中，确定古湖盆的研究层位中滩坝发育区域；确定所述滩坝发育区域中至少一口典型井的单点古水深；确定所述滩坝发育区域中的至少一口所述典型井中的th含量和u含量，并根据所述th含量与所述u含量计算th与u比值；根据所有典型井的单点古水深以及单点古水深所在位置所对应的th与u比值，拟合单点古水深与th/u值的关系方程，恢复所述典型井的定量古水深；根据单点古水深与th/u值的关系方程，恢复所述古湖盆中非滩坝发育区域其它井的定量古水深。通过该方法能够获取古湖盆同一沉积时期水体深度垂向和平面变化的连续定量古水深数据。但是此方法限于滩坝沉积体系特殊性，无法对无滩坝发育区域建立古水深与th/u值的关系方程；此外，表征古水深的关键参数选取及恢复路径相对独立。


技术实现要素：

4.本发明解决的一个主要问题是现有古水深量化恢复方法无法更直观地恢复沉积时期的全序列的水体深度，从而难以建立沉积时期古地貌和古环境的问题。
5.根据本发明的一个方面，本发明提供一种基于米兰科维奇旋回的古水深定量恢复方法，包括：
6.选取标定井；
7.基于所述标定井选取岩样，并获取所述岩样的岩相敏感测井数据；
8.对所述岩样进行微常量元素分析化验，并利用微常量元素丰度计算获得单点古水深；
9.对所述岩相敏感测井数据进行降噪，获得降噪测井数据；
10.将所述降噪测井数据分别进行multitaper method spectral analysis时间域和深度域频谱分析和bandpass滤波分析，确认所述降噪测井数据是否符合米兰科维奇天文信号；
11.所述降噪测井数据符合所述米兰科维奇天文信号后，通过evolutive harmo nic analysis&evolutive power spectral analysis演化谐波功率谱分析，并通过蒙特卡洛模拟计算演化平均光谱失配evolutive average spectral misfit，获得所述降噪测井数据
的深度域谱图；
12.基于米兰科维奇天文旋回理论，通过easmtrack程序追踪easm程序中的asm零假设显著性水平最小值，获取所述深度域谱图中的沉积速率；
13.将获取的所述沉积速率与所述单点古水深等深度对应，建立并拟合沉积速率
‑
古水深公式；
14.对所述沉积速率
‑
古水深公式进行校验，获得校验后的沉积速率
‑
古水深公式；
15.使用校验后的所述沉积速率
‑
古水深公式计算所述标定井的全序列古水深，恢复目标区域的古水深时空分布。
16.进一步地，所述对所述沉积速率
‑
古水深公式进行校验，获得校验后的沉积速率
‑
古水深公式还包括：
17.使用微常量元素丰度计算获得的所述单点古水深与利用米兰科维奇天文旋回公式计算的古水深进行比对，验证所述沉积速率
‑
古水深公式的准确性。
18.进一步地，所述对所述岩相敏感测井数据进行降噪还包括：
19.利用wavelet analysis软件中的wavelet coefficients selection 1
‑
d程序，用dmey函数将所述岩相敏感测井数据分解为9层，去除分层后的所述岩相敏感测井数据中的背景a9和最高频d1；
20.并将去除背景a9和最高频d1后所述岩相敏感测井数据保存为含有深度和表头的深度域数据以及不含所述深度和所述表头的时间域数据。
21.进一步地，将所述降噪测井数据分别进行时间域和深度域频谱分析以及滤波分析，确认所述降噪测井数据是否符合米兰科维奇天文信号还包括：
22.导入时间域数据，进行频谱分析验证特征峰频率是否满足米氏旋回周期的反比；
23.导入深度域数据，使用“astrochron”软件包中的mtm代码和bandpass代码验证所述深度域数据是否符合米兰科维奇天文周期。
24.进一步地，所述深度域谱图还包括沉积速率、天文周期和零假设显著性水平信息。
25.进一步地，所述微常量元素为co元素。
26.进一步地，所述使用校验后的所述沉积速率
‑
古水深公式计算所述标定井的全序列古水深包括：
27.利用所述沉积速率
‑
古水深公式计算所有单井的全序列古水深，进而形成井网外推全目标区域的古水深水深分布，并基于所述全目标区域的古水深水深分布，定量恢复全目标区域的各时期的古水深分布。
28.进一步地，所述建立并拟合沉积速率
‑
古水深公式包括：
29.将微常量元素丰度计算获得的所述单点古水深与对应深度的所述沉积速率整合到一起，利用origin等线性回归软件拟合两者的函数关系式。
30.根据本发明的另一个方面，还公开一种存储介质，所述存储介质为计算机可读存储介质，所述计算机可读存储介质上存储有如前任一一项所述的一种基于米兰科维奇旋回的古水深定量恢复方法。
31.沉积速率可以反映可容空间及水深变化，通过建立古水深和沉积速率之间的定量关系，实现古水深定量恢复。因此本发明采用一种结合钴元素定量计算沉积速率和米兰科维奇天文旋回理论全井段相对沉积速率来建立多点全井段沉积速率进而定量研究古水深
的方法，更直观、科学地恢复沉积时期的水体深度，从而建立沉积时期古地貌和古环境，为层序界面划分、沉积中心演化、沉降中心迁移、沉积相校验等研究提供指导。
附图说明
32.本发明构成说明书的一部分附图描述了本发明的实施例，并且连同说明书一起用于解释本发明的原理。
33.图1为本发明一种基于米兰科维奇天文旋回的全序列沉积速率
‑
古水深定量恢复方法步骤流程图；
34.图2为本发明实施例中测井数据降噪流程图；
35.图3为本发明实施例中用于识别米兰科维奇天文信号的时间域频谱分析谱图；
36.图4为本发明实施例中用于识别米兰科维奇天文信号的深度域频谱分析谱图；
37.图5为本发明实施例中用于识别405kyr长偏心率周期的滤波分析组图；
38.图6为本发明实施例中演化功率谱图、演化谐波图、eha标准化振幅图；
39.图7为本发明实施例中沉积速率能谱图以及拾取的沉积速率变化图；
40.图8为本发明实施例中验证井沉积速率、古水深投点函数关系图；
41.图9为本发明实施例中利用沉积速率
‑
古水深函数关系建立的单井古水深综合柱状图；
42.图10为本发明实施例中某层段古水深等值线分布图。
具体实施方式
43.下面将结合附图来详细描述本发明的各种示例性实施例。应注意到：除非另外具体说明，否则在这些实施例中阐述的部件和步骤的相对布置、数字表达式和数值不限制本发明的范围。
44.同时，应当明白，为了便于描述，附图中所示出的各个部分的尺寸并不是按照实际的比例关系绘制的。
45.以下对至少一个示例性实施例的描述实际上仅仅是说明性的，决不作为对本发明及其应用或使用的任何限制。
46.为使本发明的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本发明进一步详细说明。
47.对于相关领域普通技术人员已知的技术、方法和设备可能不作详细讨论，但在适当情况下，所述技术、方法和设备应当被视为授权说明书的一部分。
48.在这里示出和讨论的所有示例中，任何具体值应被解释为仅仅是示例性的，而不是作为限制。因此，示例性实施例的其它示例可以具有不同的值。
49.应注意到：相似的标号和字母在下面的附图中表示类似项，因此，一旦某一项在一个附图中被定义，则在随后的附图中不需要对其进行进一步讨论。
50.实施例一，如图1所示，本实施例为一种基于米兰科维奇天文旋回的古水深恢复方法的流程图，该方法主要包括如下步骤：
51.步骤一，选取岩心或岩屑资料丰富的井，并选择合适的岩相敏感曲线(如gr曲线、ac曲线等)进行预处理。如图2所示，为对测井曲线进行降噪的流程图，降噪后对岩心及岩屑
进行微常量元素分析化验并计算古水深。
52.步骤二，如图3
‑
5所示，将降噪后的测井数据分别进行时间域和深度域频谱分析(mtm，multitaper method spectral analysis)和滤波分析(bandpass)，识别米兰科维奇天文信号，验证该方法在该地区目的层段的可行性；
53.步骤三，确定数据符合米兰科维奇天文信号后，如图6所示，通过演化谐波演化功率谱分析(eha，evolutive harmonic analysis&evolutive power spectral analysis)和蒙特卡洛模拟计算演化平均光谱失配(easm，evolutive average spectral misfit)得到该组数据的深度域谱图，其中包含了沉积速率、天文周期、零假设显著性水平等信息，再通过easmtrack程序追踪easm程序中的asm零假设显著性水平最小值，获取如图7所示的沉积速率置信图谱(沉积速率能谱图)及沉积速率数值,；
54.步骤四，如图8所示，将天文旋回计算所得沉积速率与co元素丰度计算的古水深进行等深度对应，建立并拟合二者间函数关系式；
55.步骤五，对公式进行校验，即使用验证井岩心co元素丰度计算古水深与米兰科维奇天文旋回公式计算古水深进行比对，如图9所示，为利用沉积速率
‑
古水深函数关系建立的单井古水深综合柱状图；
56.步骤六，使用最终所得公式计算研究所有井全序列古水深，实现全区古水深时空分布恢复工作，如图10所示，为恢复后的某层段古水深等值线分布图。
57.在一些实施例中，优选地，步骤一中选取岩心或岩屑资料丰富的井以及合适的岩相敏感曲线(gr曲线、ac曲线等)进行预处理，其中对测井曲线进行降噪，对岩心及岩屑进行微常量元素分析化验并计算古水深，其具体方法为：
58.使用下述公式计算样品沉积时的速率vs：
[0059][0060]
其中，v0为近似条件现代湖泊沉积速率，n
co
为正常湖泊沉积物中钴的丰度，s
co
为样品中钴的丰度，t
co
为陆源碎屑岩中钴的丰度，t为样品中镧的含量与陆源碎屑岩中镧的平均丰度的比值；
[0061]
利用第二公式计算所述单点古水深，第二公式为：
[0062][0063]
其中，h为单点古水深的值。
[0064]
得到对应深度的点水深后，resform中的等间距数据导出(或导出后等间距处理)到excel中，新建表格去掉数据的深度和表头，导入matlab软件，使用wavelet analysis软件中的wavelet coefficients selection 1
‑
d程序进行降噪处理，用dmey函数将数据分解为9层，去除背景a9和最高频d1，处理过后导回excel并保存为含有深度和表头的深度域版本和没有深度和表头的时间域版本。
[0065]
在一些实施例中，优选地，步骤二中将降噪后的数据分别进行时间域和深度域频谱分析(mtm)和滤波分析(bandpass)，识别米兰科维奇天文信号，验证该方法在该地区目的层段的可行性，具体方法为：
[0066]
将时间域数据导入到past3.0(或其他带有时间域频谱分析功能的软件，如matlab、redfit等)中，使用timeseries栏中的redfit功能，选择90,95,99置信度计算，进行频谱分析验证特征峰频率是否满足米氏旋回周期的反比，即是否符合米兰科维奇理论；然后将深度域数据导入到r语言(或rstudio软件中)使用“astrochron”软件包中的mtm代码和bandpass代码验证是否符合米兰科维奇天文周期。
[0067]
在一些实施例中，优选地，步骤三中确定数据符合米兰科维奇天文信号后，通过演化谐波演化功率谱分析(eha)和蒙特卡洛模拟计算演化平均光谱失配(easm)得到该组数据的深度域谱图，其中包含了沉积速率、天文周期、零假设显著性水平等信息，再通过easmtrack程序追踪easm程序中的asm零假设显著性水平最小值，得到更为直观的沉积速率图并读值，具体方法为：继续在r语言(或rstudio软件)中使用上一步中的深度域数据运行eha代码、easm代码和easmtrack代码，得到相应的演化功率谱图、演化谐波图、eha标准化振幅图、沉积速率能谱图以及拾取的沉积速率变化图，并初步读值。
[0068]
在一些实施例中，优选地，步骤四中将米兰科维奇天文旋回计算所得沉积速率与钴(co)元素丰度计算的古水深进行等深度对应，建立并拟合二者间函数关系式，具体方法为：先利用吴智平等人提出、张才利等人优化的钴元素古水深公式计算古水深数据，然后将古水深数据与对应深度的米兰科维奇天文旋回计算沉积速率整合到一起，利用origin等线性回归软件拟合两者的函数关系式，针对案例研究区及近似背景研究区适用的古水深恢复公式：
[0069]
v＝36.69/h
0.714
ꢀꢀ
(3)
[0070]
其中，h为利用co元素测试的古水深，v为对应深度下利用米兰科维奇天文旋回计算的沉积速率；
[0071]
古水深恢复公式的幂函数相关系数r2＝0.96。
[0072]
在一些实施例中，优选地，步骤五中对公式进行校验，使用验证井岩心co元素丰度计算古水深与米兰科维奇天文旋回公式计算古水深进行比对，具体方法为：
[0073]
使用已建立好的沉积速率
‑
古水深公式对验证井全序列古水深进行计算，将之与co元素丰度计算的古水深进行比对，验证公式准确性。
[0074]
在一些实施例中，优选地，步骤六中使用最终所得公式计算研究所有井全序列古水深，实现全区古水深时空分布恢复工作，具体方法为：
[0075]
调整好参数后的代码理论上已适用该地区该层段的所有井位，后续将验证后的沉积速率
‑
古水深公式运用到其他单井中计算全序列古水深，进而形成井网外推全区古水深水深分布，最后恢复建立全区的各时期古水深分布。
[0076]
本发明实施例的工作原理：沉积速率可以反映可容空间及水深变化，通过建立古水深和沉积速率之间的定量关系，实现古水深定量恢复。其中建立公式所用单点古水深是通过co元素计算所得，沉积速率计算是基于米兰科维奇天文旋回理论的r语言“astrochron”天文信号处理软件包。
[0077]
本发明实施例可以针对无取心井段古水深进行计算，从而实现有测井曲线覆盖井段全序列古水深恢复，解决了特征元素法对于无取心井段无法计算古水深的问题。根据井上的信息定量计算古水深，能更直观、科学的表征不同构造
‑
沉积单元在不同时期的沉积环境以及演化，指示可容空间变化并模拟沉积过程，建立沉积水深的演化序列，并在一定程度
上对古地貌的恢复具有指导和校验意义。
[0078]
古水深的变化直接反应湖(海)平面变化，与旋回地层学具有高度关联，且研究精度可以达到数十个千年(ky)级别，为各构造单元的高频率层序格架定量构建提供了新的思路。各构造单元中使用该方法处理后井位还可以制作沉积速率
‑
古水深变化曲线，更直观、清晰地显示出湖(海)平面的变化，从而指示有利砂体分布预测及目标优选排序。并为利用地震资料和地球物理资料进行层序格架划分和有利烃源岩、储集体预测提供了一种可操作性强、直观的方法，解决了油气勘探中所涉及到的微量元素样品数不足的问题，具有很广泛的用途和广阔的市场前景。
[0079]
地球物理资料和米兰科维奇天文旋回理论相结合的是本世纪地质学家和古气候学家的最新成果和学科领域的新思路，具有可操作性好，针对性强等优点，本发明为运用测井资料恢复沉积区古水深提供了一种可操作性强的制图流程，并详细介绍了工作思路、原理、方法及流程，体现了米兰科维奇理论在古环境恢复中应用的新思路、新方法，具备先进性、勘探实用性、广泛应用性特点，可以广泛运用于层序地层学、沉积学研究和石油地质勘探中。
[0080]
以上所述仅为本发明的较佳实施例而已，并不用于限制本发明，凡在本发明的精神和原则范围之内所作的任何修改、等同替换以及改进等，均应包含在本发明的保护范围之内。
[0081]
还需要说明的是，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、商品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、商品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个
……”
限定的要素，并不排除在包括要素的过程、方法、商品或者设备中还存在另外的相同要素。