利用小波变换从测井曲线中提取米兰柯维奇周期

喻小涛，胡光明，罗水亮，奉伟，李宇 (长江大学地球科学学院，湖北 武汉 430100)

张志鹏 (中国地质大学(北京)能源学院，北京 100083)

李志威 (长江大学地球科学学院，湖北 武汉 430100)

1920年南斯拉夫学者米兰柯维奇为解释第四纪冰期成因提出米兰柯维奇理论，该理论认为北半球高纬夏季太阳辐射变化信号被放大、传输进而使全球气候发生周期性变化[1]。周期性的气候变化通过沉积物的结构、构造及沉积层的厚度等形式被记录下来，且米兰柯维奇旋回的记录已不限于第四纪，在显生宙各个时代中均留下了痕迹[2～7]。如何将米兰柯维奇周期准确提取出来，即如何从地质信息中识别出地球公转轨道三要素(偏心率、斜率和岁差)的变化周期，成为精确划分地层、校准天文地质年代和分析古气候等课题所关注的问题之一。早先学者们常用傅里叶变换来提取米兰柯维奇周期的频率，但傅里叶变换丢掉了时间信息，无法判断一个特定的信号是什么时间发生的，且相较于傅里叶变换中用正弦函数逼近信号，小波变换用不规则的小波函数来逼近尖锐信号的效果更好[8]，故小波变换逐渐成为识别米兰柯维奇周期的常用工具。笔者以鄂尔多斯盆地中二叠统山西组1段～石盒子组8段(以下简称“山1段～盒8段”)的测井曲线为例，阐述了如何利用小波变换识别米兰柯维奇周期。

1 资料选取

在二叠纪，鄂尔多斯盆地北缘的古亚洲洋消亡，大华北陆相沉积盆地雏形形成，鄂尔多斯盆地处于构造相对稳定的阶段[9,10]，南北物源形成的沉积体系在盆地中南部汇聚,形成近东西向展布的浅湖[11]。

构造与气候是陆相沉积体系发育的两大主要控制因素，构造与古气候信息都会被记录在沉积地层中。鄂尔多斯盆地二叠纪构造相对稳定，意味着沉积记录中的古气候信息受到构造(相当于噪音)的干扰较少且相对突出，有利于米兰柯维奇周期的识别，这是选取山1段～盒8段作为研究对象的主要原因。另外，考虑到相对于冲积扇、河流、三角洲而言，处于盆地汇水区的浅湖地层最为完整，因此选取位于浅湖(汇水区)的A299井、A1116井、A1118井、A1022井的自然伽马曲线进行分析，确保了气候信息的完整性，且岩性相对稳定，在一定程度上降低了岩性对米兰柯维奇周期的影响。

2 求取各层频率

以A299井为例，利用Matlab小波分析提取各层频率的操作步骤分为两大步。

1)准备数据 ①先将A299井山1段～盒8段的自然伽马(代码GR)数据置于Excel表的一列中，第一行为A299GR，保存为A299GR.xls；②在Matlab主页中选取“导入数据”，选择A299GR.xls文件并打开；③进入导入窗口，选中自然伽马整列数据，然后在视图窗口的菜单中选择“列向量”；④点击菜单中“导入所选内容”；⑤返回Matlab主菜单，在工作空间中找到名为A299GR的向量，右击另存为A299GR.mat文件。

2)求取各层频率 ①运行Wavemenu，开启小波分析窗口，点击Continuous Wavelet1-D，在Continuous Wavelet1-D窗口的File下拉菜单中选择Load Signal，打开A299GR.mat文件；②Wavelet选择db5，Scales setting选择Power 2 Mode，Power选择10，点击Analyze；③根据a=2n(其中，a为尺度，n为层数)，在系数图中右击选择相应的层数(据选择的尺度推算)，然后点击New Coefficients Line，读出该层的频率(Frequency)；④重复步骤③选择下一个层数并读出相应的频率，依次读出d1～d10层的频率(见表1)。

表1 各井不同尺度下的相对频率

表1显示各井的同层频率是一致的，说明该地区各井所处位置的沉积物受到相同外部因素的影响。

为了便于与天文周期进行对比，将各层的频率比上最小频率(0.001)得到相应的值(见表2)。由于d9与d10对应的频率一样(表1)，说明d9层是软件所能分析的极限了，d10层的数据不采用。

表2 不同尺度下的相对频率与最小相对频率的比值

3 求取理论周期与频率

周期与频率互为倒数关系，得到了周期值，也就得到了频率值。地质历史中某一时期的偏心率、斜率和岁差周期的比值是相对固定的，如果能在研究区目的层段的系列周期(或频率)中找出3个周期，其比值与理论米兰柯维奇周期(或频率)一致，即识别出了米兰柯维奇周期(或频率)。

偏心率周期405ka被认为是天文上最为稳定的地球轨道参数，在过去几亿年是不变的[12]，因此选取405ka做为偏心率周期。地质历史中斜率和岁差并不是固定的，中生代以来的斜率和岁差可以根据公式计算[13,14 ]，但是更古老的地层，其中的干扰因素较多，需要通过其他方式求取。

Hinnov和Hilgen计[14]算出了过去250Ma(大致为中生代与古生代的分界)以来斜率和岁差的变化(表3)。在该基础上，以时间为横坐标，不同的天文周期为纵坐标，拟合出相应的函数关系式(图1)。

表3 斜率和岁差周期在过去250Ma的变化(据文献[14]，有修改)

注：α1、α2、α3、α4、α5分别代表斜率1、斜率2、斜率3、斜率4、斜率5；P1、P2、P3、P4、P5分别代表岁差1、岁差2、岁差3、岁差4、岁差5。

图1 地质历史时期的斜率周期与岁差周期

根据2014版的中国地层表和2016版的国际年代底层表，确定山1段～盒8段的地质年龄界限为272.3±0.5Ma，取其中间值272.3Ma代入图1中的各个函数关系式中，计算出相应的斜率周期和岁差周期分别为(单位ka)：

α1=38.95295，α2=32.03703，

α3=31.12670，α4=24.81591，

α5=24.14775;P1=20.40117，

P2=19.44416,P3=16.99083，

P4=16.82506,P5=14.82297。

4 确定米兰柯维奇周期

将272.3±0.5Ma的偏心率、斜率和岁差周期转换为频率，计算各频率与最小频率的比值，取整数(表4)。将表4中频率比值取整数后，与表2对比可知d9∶d6∶d5=1∶10∶21=e∶α1∶P2，即d9、d6、d5的频率就是研究区要提取的天文周期所对应的频率，分别代表e、α1、P2，据此可以提取与偏心率、斜率和岁差相对应的曲线。需要特别说明的是，表4中的频率为真实频率，而表1中的频率则为范围在0～1之间的相对频率。因此，d9、d6、d5的频率与e、α1、P2并不相等，只是它们之间的比值是相等的。

表4 目的层(272.3±0.5Ma)对应的周期及其相应的处理数据

图2 A1022井的3个天文周期

所提取的偏心率周期为405ka，斜率周期为38.9ka，岁差周期为19.4ka。岁差的周期：斜率的周期：偏心率的周期=1∶2∶21。因此，对于某一段地层来说，其包含的岁差周期、斜率周期、偏心率周期的个数的比值应接近21∶10∶1，即若某一段地层中含有1个偏心率周期，则含有的岁差周期个数为21个，斜率周期个数为10个。如果该次研究提取的偏心率、斜率、岁差是正确的，那么就应该找到这样的地层，使得地层中只要含有1个偏心率周期，就有对应的21个岁差周期和10个斜率周期。从图2可以知道，A1022井基本符合要求，图中红线所框选的部分中，从下往上看包含1个偏心率周期，10个斜率周期，18个岁差周期，基本满足该比值，从而验证了所提取的天文周期的合理性。

5 结语

天文周期普遍存在于沉积地层中，利用Matlab中的小波变换工具从自然伽马曲线中提取米兰柯维奇旋回对应的周期，首先要将目的层段的自然伽马曲线转换为“.mat”格式，再使用Continuous Wavelet 1-D将其分解为10层，读取各层的频率，并计算各层频率之间的比值，然后将该比值与该时代的理论频率比值进行对比，从而识别出d9、d6、d5等对应的频率和周期为研究区山1段～盒8段偏心率、斜率和岁差所对应的频率和周期。前人在应用米兰柯维奇周期时，往往应用的是偏心率，而该次研究成果是从一段测井曲线中提取相应的3个米兰柯维奇周期(斜率、岁差、偏心率)，对于精确划分地层，校准天文地质年代具有重要的意义。

[参考文献]

[1]吴瑞堂，张守信. 现代地层学[M]. 武汉：中国地质大学出版社，1991：2～6.

[2]吴智勇. 米兰柯维奇韵律层及其年代地层意义[J]. 地层学杂志，1995，19(2)：156～160.

[3]金之钧，范国章，刘国臣. 一种地层精细定年的新方法[J]. 地球科学(中国地质大学学报)，1999，24(4)：379～381.

[4]童金南，殷鸿福. 浙江长兴煤山剖面Griesbachian期旋回地层研究[J]. 地层学杂志，1999，23(2)：130～134.

[5]郭刚，童金南，张世红，等. 安徽巢湖早三叠世印度期旋回地层研究[J]. 中国科学D辑(地球科学)，2007，37(12)：1571～1578.

[6]龚一鸣，徐冉，汤中道，等. 广西上泥盆统轨道旋回地层与牙形石带的数字定年[J]. 中国科学D辑(地球科学)，2004，34 (7)：635～643.

[7]苏德辰，李庆谋，罗光文，等. Fischer图解及其在旋回层序研究中的应用——以北京西山张夏组为例[J]. 现代地质，1995，9(3)：279～283.

[8]邵龙义，何志平，鲁静，等. 环渤海湾西部石炭系-二叠系层序地层及聚煤作用研究[M]. 北京：地质出版社，2008：111～112.

[9]周安朝，贾炳文，马美玲，等. 华北板块北缘晚古生代火山事件沉积的全序列及其主要特征[J]. 地质论评，2001，47(2)：175～182.

[10]杨俊杰. 鄂尔多斯盆地构造演化与尤其分布规律[M]. 北京：石油工业出版社，2002：3～5.

[11]田景春，吴琦，王峰，等. 鄂尔多斯盆地下石盒子组盒8 段储集砂体发育控制因素及沉积模式研究[J]. 岩石学报，2011，27(8)：2403～2411.

[12]黄春菊. 旋回地层学和天文年代学及其在中生代的研究现状[J]. 地学前缘，2014，21(2)：48～62.

[13]Laskar J，Robutel P，Joutel F，et al. A long-term numerical solution for the insolation quantities of the Earth [J]. Astronomy and Astrophysics，2004，428(1)：261～285.

[14]Hinnov L A，Hilgen F J. Charpter 4：Cyclostratigraphy and astrochronologgy [A].Gradstein F M，Ogg J G，Schmitz M D，et al. The Geologic Time Scale 2012[C]. Amsterdam：Elsevier，2012：63～83.