天文地层学研究程序及其在渤海湾盆地东营凹陷的应用实例分析

彭 军，于乐丹，许天宇，韩浩东，杨一茗

（1.西南石油大学地球科学与技术学院，四川 成都 610500；2.中国地质调查局成都地质调查中心，四川 成都 610081）

旋回地层学（cyclostratigraphy）一词最早由A.G.Fischer在意大利佩鲁贾（Perugia）和法国迪涅（Digne）举办的全球沉积地质大会上正式提出。它是对地层记录的（准）周期性变化进行识别、描述、对比和成因解释，以期提高年代地层格架的精度，进而能更高精度地开展地层高精度划分与对比的一门地层学分支学科［1-3］。天文地层学是由旋回地层学等演变而成的一个新名词，强调了由地球引力场的周期性摄动导致地球表面接收到的太阳日照量发生周期性变化，并引起沉积地层序列相应的周期性变化，且将地层旋回性的最终驱动力归结为地球轨道参数（偏心率、地轴斜率及岁差）的周期性变化［3-5］，从而避免和旋回地层学相混淆。综合来看，天文地层学是研究受地球轨道周期性变化控制形成沉积地层序列的地层学分支学科，是以米兰科维奇理论为基础的一门新兴的地层学分支学科，米兰科维奇旋回序列是天文地层学的研究重点［5］。

天文地层学在解决重大地质事件问题、探究气候和生物的演化以及分析全球碳循环等方面都起到了关键作用，即天文轨道驱动力是解决这些问题的重要驱动力［6-11］。在天文地层学研究过程中，经常被提及的一个问题就是米兰科维奇旋回信号是否存在，而天文地层学研究就是通过技术方法和研究手段将地质信号中的干扰剔除，提取出有效的地层序列替代参数，发现记录在沉积地层中的米兰科维奇旋回信号，即偏心率、斜率及岁差周期等轨道参数［12-13］。天文地层学研究方法主要分为两大类，包括岩性直观识别法和时间序列分析法，前者主要依据露头剖面和钻井岩心资料，受此类资料不连续和取样密度的限制，其运用难度大且易造成旋回信息的遗漏，同时还要求研究者具有很强的地质专业知识和丰富的研究经验；而后者是一种数字处理技术，得到了广泛的应用。该方法主要基于广泛应用的数字信号处理技术，基本步骤包括数据的预处理、频谱分析、时-频分析和相关性研究［14-18］，而这些研究的前提是借助大量的先进测试技术手段获取可用于天文旋回研究的连续数据序列作为沉积地层序列的替代参数（以下简称“替代参数”），即选择能够反映古气候变化的测井参数曲线（主要包括自然伽马、自然电位、声波时差、电阻率等）、碳-氧同位素、磁化率和各种地球化学分析数据以及高分辨率的连续元素录井数据等［19-26］，并应用天文旋回理论来研究这些数据资料中反映的轨道周期。再结合放射性同位素年龄、生物地层和磁性地层年龄，建立更加精细的天文地质年代标尺，在精细时间尺度上理解沉积旋回的地质意义，以期达到精细定年的目的。其后在单井单剖面天文旋回划分基础上，开展高精度的等时地层对比，建立高时间精度的等时地层对比格架，为编制高精度大比例尺等时沉积微相平面分布图、探讨储层的时空展布及演化、预测生储盖组合等奠定地质基础［27-32］。

近几年，天文地层学发展迅速，但由于存在多种替代参数及时频分析方法，导致诸多学者给出的同一研究区天文旋回划分结果存在多解性，致使其划分结果的可靠性存在争议［33-36］。因此，本文调研了近几年国内外米兰科维奇旋回的研究成果和进展，从替代参数的选取、时频分析方法的运用及天文调谐等3个方面入手，阐述各替代参数和分析方法的优缺点，总结天文旋回划分研究中存在的问题及今后的研究方向，并以渤海湾盆地济阳坳陷东营凹陷沙河街组三段和四段天文旋回划分作为具体的研究实例，阐述天文地层学的研究思路和方法、研究流程，优选适用的替代参数和数据处理方法，为其他地区天文旋回研究提供参考和借鉴，推动天文地层学的向前发展。

1 地球轨道参数

太阳系中地球与太阳、月亮和其他行星之间的万有引力作用导致地球绕太阳的公转运动和地球的自转轴发生（准）周期性的摄动，表现为地球轨道参数的偏心率、地轴斜率和岁差发生近似周期性的变化。

1.1 地球轨道偏心率

地球轨道偏心率是指地球公转轨道由于受到太阳之外的其他星体引力作用而偏离正圆轨道的程度，定义为地球绕太阳公转椭圆轨道的赤道半径与极半径之差再与赤道半径之比，其值在0.000 21～0.067 00周期性波动，现在为0.016 70［5］（图1）。偏心率主要周期有405，95，99，124，131，2360和1 000 kyr。偏心率的变化主要受地球和太阳的相对位置影响，会对地球公转的椭圆轨道产生影响，其循环形式造成了日地距离微小的季节性变化，即当偏心率较大时，地球绕太阳运动轨迹越趋近于椭圆形，地球整体接受日照量多，为间冰期，气候暖湿，同时季节变化明显；当偏心率较小时，地球绕太阳运动轨迹越趋近于圆形，地球整体接受日照量少，为冰期，气候干冷，同时四季变化不明显［3-4，37］。

1.2 地轴斜率

地轴斜率是指地轴的倾斜度，它是地球公转的轨道面与地球赤道面的夹角，其值在22.5°～24.5°变化，现在斜率值为23.26°。现在斜率的主要周期为41.0 kyr还包括39.6，54.6和29.0 kyr次要周期［3-5］（图1）。斜率大小影响地球上不同纬度日照量分布均匀程度，以及某一固定纬度不同季节气候差异性大小，即地轴斜率增大，太阳直射点在南回归线与北回归线间移动，全球日照量分布均匀，且某一固定纬度一年之中日照量差异明显，高纬度夏季日照量越多，冬季日照量越少，四季分明；若地轴斜率为零，太阳直射赤道，一年之中日照量差异小，不同纬度间气候差异大，但对于某一固定纬度，高纬度夏季日照量越少，冬季日照量越多，几乎无四季之分［3-5，37-40］。

1.3 岁差

岁差是指在固体潮汐作用下导致“地球差异旋转”，即地球自转轴的进动（图1），从而导致回归年与恒星年产生时间差的天文现象，它能够反映偏心率和地轴斜率的综合影响。岁差的主要周期为23.7，22.4，19.0和17.0 kyr［3-5］。岁差本身并不改变地球日照量的分布，其对地球气候产生影响主要受偏心率的调控作用。岁差决定了季节发生的时间，当岁差较小时，对于北半球来说，夏季日照量极大，冬季日照量极小，即北半球冬至日到达远日点，夏至日位于近日点，经历了一个短暂而炎热的夏季和一个漫长而寒冷的冬季，此时，由于整体接受日照量少，气候表现为干冷，但是由于冬夏日照量差异大，表现为季节性强；反之，当岁差较大时，北半球经历了一个漫长而温暖的夏季和一个短暂而温暖的冬季，气候表现为暖湿但季节性弱［3-5，35，37］。

图1 地球轨道参数地球轨道偏心率、地轴斜率和岁差周期示意图（据参考文献［5］修改）Fig.1 Schematic diagrams showing the Earth orbital parameters，i.e.the eccentricity，obliquity and precession period（modified after reference［5］）

2 沉积地层序列替代参数

沉积地层序列替代参数指的就是能够反映过去沉积介质条件变化尤其是气候变化的各类古生物、地球物理和地球化学等参数，由于气候变化与轨道参数之间具有密切联系，因此从理论上来说，各类与气候变化相关联的指标均可作为替代参数用于天文地层学研究［3，41］。通过调研前人研究成果可知，用于天文旋回分析的替代参数包括古生物参数、环境磁学参数、地球物理参数和地球化学参数等，通过对各类参数进行综合分析，从而构建包含地层环境变化信息的时间序列。

2.1 古生物参数

众所周知，生物丰度、生物灭绝速率和微体化石组合特征可作为层序界面、体系域界面的识别标志，从而为层序及体系域的划分提供新的研究思路。前人在旋回识别中多运用虫孔及生物扰动、腕足类等低镁方解石生物壳、牙形石带、孢粉、有孔虫、钙质超微、沟鞭藻等古生物化石，生物丰度、种类及分异度随着海/湖平面的升降表现为周期式的增减，基本反映了海/湖平面升降旋回性变化［42-45］。尤其值得关注的是，湖相沉积中的孢粉能够反映古气候条件及其变化，同时，结合孢粉母体对生长温度及湿度的不同适应性，将孢粉分别划分为喜热组、喜温组、广温组及湿生组、旱生组和中生组，基于孢粉的发育情况对气候条件进行表征［37］。古生物参数在旋回划分中能够发挥重要作用，但是古生物化石的系统取样很难，在细粒沉积岩中尤其困难，且在数量较少的情况下，生物数据又无法很好地说明问题，因此，古生物参数的运用要与其他旋回划分方法相结合，开展综合分析。

2.2 地球物理参数

在用于天文旋回分析的替代参数中，测井数据（或测井曲线）是最理想的参数之一。测井数据不仅分辨率高，达到分米级尺度，而且能够敏感、连续地反映所测地层的岩性及物性特征［46］。不同规模、不同尺度地球轨道参数的周期性变化必定引起气候变化，这种变化会记录在沉积物的结构、构造、岩性及岩相等周期性变化上，并反映在测井数据中，表现为不同幅度、形态及频率特征［3，41］，刻录了地层的沉积旋回信息，故测井数据是分析地层天文旋回的主要资料。不同类型测井数据虽然对地层物理特征敏感程度不同，但所反映的地层沉积旋回是一致的［5］。

调研前人研究成果可知，在天文旋回分析中可以运用单一测井曲线或多种测井曲线相结合来识别旋回，而使用最多的是自然伽马测井数据［5，17，22-23，47-49］，其次是将自然伽马和自然电位测井数据相结合［16］，还有学者采用自然伽马、自然电位和电阻率数据相结合的方法［15，50］，部分学者采用自然伽马能谱测井识别旋回［51-55］，少数学者采用了自然电位测井数据识别旋回［56］。综上可知，多数学者选用自然伽马测井曲线识别米兰科维奇周期，该类曲线对岩性及岩相的变化较为敏感，能够较好地反映其随深度的变化情况，从而能反映沉积环境的变化，因此自然伽马测井数据已经被广泛用于古环境、古气候的重建，在天文旋回地层研究中得到了广泛的应用。

2.3 地球化学参数

天文轨道参数的周期性变化会引起地球表面气候的变化，进而引起地层记录中地球化学的相应变化，因此，以能表征沉积介质条件的地化数据作为替代参数开展天文旋回分析是可行的。目前，用于天文旋回分析的地球化学数据主要包括碳和氧同位素、有机地球化学和元素分析数据等。δ13C和δ18O同位素数据主要应用于碳酸盐岩和深水沉积体系的旋回划分研究［42，57］。在岩心和露头资料缺乏、地震分辨率受限的情况下，利用有机地球化学数据可实现对深水细粒沉积地层精细的层序单元划分［58-59］。元素地球化学数据既可用于碳酸盐岩旋回地层的划分，也可用于陆相细粒沉积地层旋回的识别［21，60］。从理论上来说地球化学参数在天文旋回研究中的适用范围是非常广泛的，但从已有的研究成果来看，实际使用的却较少，这是由于测试技术及成本和样品的限制，地球化学数据的获取具有局限性，很难获得高分辨率且连续性较好的数据系列进行天文旋回分析。需要特别注意的是，在研究混积岩发育的地层时，其同位素或元素含量会由于火山碎屑等物质的注入发生失真，使得参数值发生畸变，故在使用该参数开展天文旋回分析时，要与其他旋回分析方法结合，相互验证。

2.4 环境磁学参数——磁化率和非磁滞剩磁

以环境磁学参数作为替代参数的天文地层学研究认为，地质历史时期的环境解释可以运用岩石磁性参数这一指标，其主要是依据磁性矿物含量以及磁性矿物粒度的大小建立与沉积环境之间的关系，并对包括磁性矿物含量在内的一些数据进行时间序列分析，进而利用环境磁学参数——磁化率（MS）和非磁滞剩磁（ARM）作为替代参数进行天文地层学的研究［18，55，60-64］。由于天文旋回的研究对替代参数数据的连续性要求高，需要大量等间距采集样品开展测试分析，这就要求测试分析方法简单、经济实惠且易于解释，而环境磁学参数的测量较为方便且成本低，正好适用于天文地层学分析，因此高分辨率磁化率数据作为沉积地层序列的替代参数，与自然伽马测井数据一样，应用范围较为广泛［65-69］。但该方法的不足就是测试时可能受到人为因素的干扰，因此在测试时应创造良好的测试环境。

3 天文地层学分析的研究方法

天文地层学研究多选择在地层连续且露头良好的剖面上进行，而常用的研究方法分为岩性直观识别法和时间序列分析法。岩性直观识别法是基于沉积学理论的岩性和岩相韵律性变化的识别和划分，一般用于出露完整、有显著韵律性的野外剖面和连续的长井段岩心，直接用肉眼观察其岩性组合特征，识别出纹层等不同级次的旋回，分析地层中是否记录有天文轨道周期信号。时间序列分析法是一种广泛应用的数字信号处理技术，该方法主要涉及地层序列的数据选取、预处理、频谱分析、时频分析和相关性研究，从而实现多尺度、多分辨率信号分解［18］。各种地层学（放射性同位素年代地层学、生物地层学和磁性地层学等）的约束是天文旋回分析和可靠性验证的重要依据［3］。

3.1 岩性直观识别法

天文地层学研究最直接的方法就是通过观察岩性和岩相的变化来识别旋回界面，从而达到划分旋回的目的。对于沉积岩性差异较大的韵律性地层而言，通过对连续露头或岩心的观察，可以识别出相似的沉积韵律组合和不同级次的韵律组合特征［70］。前人研究表明，可通过纹层/层等的组合特征、方式及级序结构来判别沉积过程是否受到地球轨道参数的影响［3，71］；通过露头及岩性测井资料识别韵律性湖进-湖退沉积充填序列，探究其与米兰科维奇旋回的关系［47，72］；运用岩心资料，结合地球化学和磁化率数据，分析湖相细粒沉积地层（岩相）叠加样式，划分沉积旋回，并探究基于米兰科维奇理论的高频层序定量划分［25，34，73］。从前人的研究成果可以看出，旋回的岩性直观识别法需要与其他数据或方法相结合使用，尤其对于深海相和部分湖相细粒沉积岩岩性的差异往往肉眼难以识别，需要借助其他替代参数和岩性识别手段共同开展旋回分析；另外在缺乏精确定年的绝对年龄数据情况下，也很难确定旋回界面的级次。

在东营凹陷樊页1井沙四上亚段纯上次亚段（深度3 251.00～3 441.00 m）的岩心观察中，可见到由纹层状构造、块状构造或条带状构造组成的韵律层。以3 380.94～3 382.94 m深度段为例，自下而上划分出由厚层灰色纹层状灰质粘土岩（A）和薄层深灰色块状粘土岩（B）组成的8.5个韵律组合（图2）。该韵律组合的识别，为运用时间序列分析法识别旋回提供了沉积学证据，测井时间序列研究得出该韵律组合发育2个22.00 kyr岁差尺度的天文旋回，进一步通过沉积学的精细研究，可以开展更小时间尺度（18.82 kyr和18.69 kyr）天文旋回的识别与划分。

图2 东营凹陷樊页1井沙四上亚段纯上次亚段（埋深3 380.94～3 382.94 m）岩心韵律特征Fig.2 The core rhythmcharacteristics of the pure upper sub-interval of the E s4U in Well Fanye 1 in the Dongying Sag（at burial depth of 3 380.94-3 382.94 m）

通过岩性观察来识别天文旋回的级序结构和组合特征，要求研究者具有扎实的地质基础知识和丰富的野外工作与岩心观察经验［41］，要搞清楚露头剖面和岩心中频繁变换的岩性和岩相，从而识别出不同级次的天文旋回［3］。然而，并非所有的轨道参数变化都能够引起岩性或岩相的明显变化，尤其在深湖-半深湖、远洋-半远洋等环境的细粒沉积岩组合中，依靠岩性直观识别方法易于造成旋回信息的遗漏［3，41］。目前随着各种技术手段、仪器设备和实验方法的进步，露头剖面工作和岩心观察不再局限于肉眼，可以借助一些仪器设备和测试分析手段将剖面和岩心地质信息或样品采集回实验室测试转化为连续的数据信息，并对数据进行处理来开展天文旋回分析，这就要借助于一种广泛应用的数字信号处理技术——时间序列分析法。

3.2 时间序列分析法

在天文地层学研究领域，时间序列分析法得到了广泛应用，无论是野外露头剖面还是岩心观察，时间序列分析方法的应用，其重要的第一步就是要获取用于研究的连续数据序列，进而开展时间序列分析。也就是说，运用该方法开展旋回识别的关键就是要把深度域数据向频率域转化，再完成频率域到时间域的转变，通过一系列的时频分析方法获得可以与天文轨道周期参数理论值对比的频谱来划分旋回。

3.2.1 数据取样及预处理

时间序列分析最基础和最关键的工作是确定采样密度和采样间距，采样密度过疏会造成旋回遗漏而无法确保旋回的完整性，采样密度过密则会费时费力且成本高。因此，在采集数据时，多运用奈奎斯特（Nyquist）采样定理，从而保证最短的旋回至少有2个点控制［3，74］。在具体操作中，可根据现有的生物地层、古地磁等确定的初始年代框架估计出沉积速率，进而确定出采样密度。例如，如果一个岁差旋回代表20 kyr，沉积速率是5 cm/kyr，对应的岁差旋回厚度约100 cm，取样间隔至少应控制在50 cm以下，而这种密度是旋回分析数据采集的最低标准，而实际上，一个20 kyr的岁差旋回有3～4个数据点控制是较为合适的，也就是采样间距在5～7 kyr，取样间隔就是25～30 cm，这样既不会造成不必要的时间和资金浪费，也可以保证获得真实的旋回信息［18］。

由于对古气候变化反应敏感的各种沉积序列替代性参数中，会包含一些除了正常信号之外的干扰信号，即“噪音”，因此，要对获取的数据进行预处理。预处理的方法包括插值、消噪和去趋势化［3，74］，插值是为了对用于时间序列分析的不等间距数据进行等间距处理；消噪是为了去除信号中除了轨道周期信号之外的所有因素造成的信息干扰，当噪声过大时会使轨道周期信号难以识别［75］；去趋势化是指消除由于构造和沉积作用影响而出现的假长周期，这种周期的存在会导致真实长周期被压制而漏失。本文借助Matlab软件提供的一维连续小波工具箱，对樊页1井沙三下亚段（深度3 051～3 251 m）自然伽马和电阻率测井数据开展消噪和去趋势化的预处理，通过分析预处理后的数据可知，自然伽马测井数据信号“噪音”少，消噪和去趋势化后的曲线基本与原始曲线一致；而电阻率曲线存在的“噪音”较大，消噪后去除了一些异常值点，消除了一些假周期信号的干扰。

3.2.2 频谱分析

识别地层中米兰科维奇旋回信号的关键一步就是对沉积地层序列进行频谱分析。将替代参数数据连续的深度域数值序列看作时间域，把时间序列看成频率组合，便可对替代参数数据进行属性提取分析，其目的就是要确定出时序信号中周期性或准周期性的成分［3，76］。将优势频率取倒数就得到优势旋回厚度，若旋回厚度的比值与理论地球轨道参数的长偏心率、短偏心率、斜率和岁差周期之比近似1∶2∶5∶20，就可以初步判定该数据的时间序列分析结果记录了米兰科维奇旋回信息［3］。

频谱分析研究常采用包含多种主流频谱分析方法的Past 3.0统计分析软件，选择其Spectral analysis模块中的REDFIT软件包开展分析，与其他方法相比，REDFIT程序的一大优势就是无需插值就可直接应用［16，21，47-48，50，76］。图3为东营凹陷樊页1井沙三下亚段自然伽马测井数据的频谱分析结果图，即在Significance lines窗口中选择红噪曲线（图3中的绿色虚线）和90%置信区间曲线（图3中的黑虚线）。图中横坐标表示频率，地层旋回厚度对应于优势频率的倒数，纵坐标表示频谱能量即频率幅值，为某频率的显著程度。结合频谱分析结果，提取图中的优势频率成分，运用频率与其对应的旋回厚度之间的倒数关系求出相应的旋回厚度，最终通过优势旋回厚度比值与理论轨道周期比值对比判断所研究地层是否受到天文轨道周期的驱动。

图3 东营凹陷樊页1井沙三下亚段GR曲线Redfit频谱分析Fig.3 Redfit spectrumanalysis diagramof GRcurve of the E s3L in Well Fanye 1 in the Dongying Sag

频谱分析方法实现了周期信号从时域到频域的转化，但无法从时频双域对周期信号进行观察分析，也不能反映周期性的深度位置，因而不能从整体上划分出不同级别的沉积旋回界面［3，75，77］。

3.2.3 小波变换和功率谱分析

小波变换能通过不同尺度和频率对信号进行观察处理，是一种时频双域分析方法，即小波变换能够将复杂信号分解为不同频率或者周期的旋回曲线，被分解为各自周期独立的沉积旋回且以不同尺度的形式展示出来，从而考察小波时频能量图局部能量团的变化和多种伸缩尺度的周期性振荡特征［25，77-79］，通过分析旋回性与各级层序界面建立对应关系，这是运用小波变换划分沉积旋回的依据［51，56，80-82］。

图4是借助Matlab软件提供的一维连续小波工具箱，运用Morlet小波作为子波对东营凹陷樊页1井沙三下亚段的GR数据进行小波分析处理的结果，该图从上到下依次为：樊页1井沙三下亚段预处理后得到的GR测井信号；经过Morlet一维连续小波变换得到的时频能谱图，红色代表高能量，蓝色代表低能量；尺度值为512时的小波变换系数曲线；时频能谱图所对应的小波变换系数等值线图，纵坐标代表不同小波尺度。在小波能谱图中，颜色越亮代表能量越强，一个亮色的能量环即代表一个周期［83］，而频率的成分相对单一稳定，结合频谱分析等结果可进一步证明研究层段存在较为明显的米兰科维奇旋回。

图4 东营凹陷樊页1井沙三下亚段一维连续小波处理结果示意图（尺度值=512）Fig.4 Schematic diagram of one-dimensional continuous wavelet processing results of the E s3L in Well Fanye 1 in the Dongying Sag（scale value：512）

小波变换方法在识别不同级别沉积旋回方面具有多分辨率的特点，克服了频谱分析不能反映时域、频域局部性特征的缺点，它可以将不同的频率（周期）成分的变化情况较为精确地显示在时间（深度）域上，进而能更好地识别存在的沉积间断和反映沉积速率的变化［18，41］。但运用小波时频能谱图划分沉积旋回界面时会存在一定误差，往往使界面读取不准确。因此，在小波变换分析的基础上，采用Morlet小波基对测井曲线进行小波多尺度分解，即借鉴一种基于小波变换的功率谱分析提取不同级别旋回的方法，以便能反映地层沉积单元级别从小尺度到大尺度的能量变化［77］。

为了进一步判断这些相对单一稳定的频率成分是否受某一具体天文轨道周期的控制，在小波分析的基础上开展功率谱分析，提取小波系数矩阵以计算模极值，将不同尺度的能谱值累加作相应的模极值图［83-84］。图5是将不同尺度能谱经过矩阵计算得到的东营凹陷樊页1井沙三下亚段模极值图，从图中可看出存在两个明显的模极大值点，从小波能谱图中也可发现，亮色能量环主要集中于尺度值108和331两线附近，说明了这两个尺度值附近的频率为信号中的优势频率，可选择这两个极值尺度用于分析信号的小波周期。

图5 东营凹陷樊页1井沙三下亚段GR曲线小波分析模极值Fig.5 Diagramshowing the module maximum from GRcurve wavelet analysis of the E s3L in Well Fanye 1 in the Dongying Sag

3.2.4 短时傅里叶变换

短时傅里叶变换、小波变换、S变换等均属于线性时频分析方法。短时傅里叶变换是一种固定时间窗的时频分析方法，它用窗函数来截取信号，假定信号在窗内是平稳的，采用傅立叶变换来分析窗内信号，确定该时间存在的频率，然后沿着信号时间移动窗函数，得到信号频率随时间的变化关系［85］。短时傅里叶变换方法中窗口是固定的，时频网格等宽，不能随频率需要使网格变宽或变窄，不利于低频、高频信号的检测；时间和频率的精度无法兼顾，使用长时间滑动窗，在频域获得高精度，但时域精度降低；反之，使用短时间滑动窗，在时域获得高精度，又会使频域精度降低［86］。在运用短时傅里叶变换开展时频分析时，多与小波变换方法相结合，从而更好地识别沉积地层中存在的轨道周期［87-90］。

3.2.5 S变换和广义S变换

在S变换时频分析中，首先将原始信号从时间域经过S变换变到时间-频率域，在对信号进行分析后，又将时频信号变换到时间域，体现了S变换的无损可逆性［91］。S变换构建的高斯窗函数包含与频率有关的因子，频率的倒数决定了窗函数的大小，而当频率变大时，其窗函数的宽度可以自适应变小，一般来说，高频信号时窗较窄，时域分辨率高，低频信号时窗较宽，时域分辨率低，当频率变大时，S变换中所构建的高斯窗函数宽度就变窄，对应的信号幅值增大［91］。而广义S变换比S变换具有更好的时频分辨率，并且具有多分辨率特性。小波变换因具有更多的小波类型作为信号分析的选择，因此小波变换和广义S变换对数据处理具有更好的灵活性，同时，该类时频分析方法除了单独使用外，还常与小波变换时频能谱图结合使用识别旋回［92-96］。

3.2.6 旋回划分结果检验方法

在对各种地球物理参数和地球化学参数等旋回分析的替代参数开展时频分析时，要特别注意验证旋回分析结果是否受到天文轨道周期的控制。天文旋回划分结果的检验常用的是理论轨道周期比值法，还有学者采用假设检验法和天文调制法。

在距今500 Myr中，地轴斜率和岁差的周期值都在不断的变化，而地球轨道偏心率周期值，尤其是振幅最大的404 kyr周期值，受金星和木星万有引力的影响，在地质历史时期非常稳定，变化较少［5］。因此，在运用理论轨道周期比值法时，要先研究地层所处的地质时代，运用Laskar具有代表性的轨道要素计算方法［97］，计算得出用于研究层段地质时代的理论轨道周期。各级天文周期间存在相对稳定的比例关系，这种比例关系可保存于稳定沉积环境中相对整合的沉积地层里，也就是说各级天文周期的比值响应于相对整合的各级旋回厚度的比值。一般而言，天文轨道的长偏心率、短偏心率、斜率和岁差周期的比值近似为20∶5∶2∶1，故在频谱分析获得优势旋回厚度之后，将优势旋回厚度的比值与理论轨道周期比值作对比，如果比值一致或者相近，就可以认为沉积地层受到天文旋回的控制，该方法目前运用最多、范围最广［35，47，52-53，98-100］。

根据前期研究可知，东营凹陷沙四上亚段至沙三中亚段地层年龄在50～30 Ma［74］，故本次采用Laskar的方案计算了北纬38°夏至日（6月21日）50～30 Ma期间轨道周期变化的理论值，最终计算出偏心率、斜率和岁差周期之间的比例关系（表1）。而由图3的频谱分析结果可知，主要峰值的旋回厚度分别为38.95，12.98，4.10，2.40和1.95 m等，旋回厚度比值为38.95 m∶12.98 m∶4.10 m∶2.40 m∶1.95 m=19.974∶6.656∶2.103∶1.231∶1，该比值与古近纪米兰科维奇旋回理论轨道周期405 kyr∶124.22 kyr∶39.76 kyr∶23.28 kyr∶22 kyr的 比 值18.409∶5.646∶1.807∶1.058∶1很接近。因此进一步确定樊页1井沙三下亚段受米兰科维奇旋回控制。

表1 50～30 Ma地球轨道参数周期及其比值关系Table 1 Earth orbital parameter periodsand their ratios during the 30—50 Ma

除理论轨道周期比值法之外，少数学者还使用天文调制法［101］和假设检验法［102］验证旋回识别的结果。但是天文调制法由于调制过程中人为因素干扰大而不可单独使用；假设检验法则是要先提出信号与天文旋回无关的假设，然后再拒绝该假设，但该方法对信号的稳定性要求很高，若信号失真则会产生很大误差。综上可知，旋回识别的检验方法，首选的是理论轨道周期比值法。

3.2.7 调谐方法

天文调谐就是将替代参数的数据序列对比到偏心率、斜率或岁差理论目标曲线上，即分析由替代参数获得的米兰科维奇旋回与理论轨道周期曲线之间的关系，然后将放射性同位素测年数据、古地磁数据、生物化石记录确定的年龄数据作为“锚点”，将天文轨道周期值与放射性同位素年龄、生物地层年龄、磁性地层年龄相互校正，最终建立精确的天文年代标尺［3，18，54，103-106］，从而为地质年代提供更精确的时间约束。

滤波是对信号中特定波段频率滤除的操作，在天文地层学天文调谐研究中应用广泛。一般来说，数字滤波的方法有低通滤波、高通滤波以及帯通滤波［3］。高通滤波和低通滤波可滤除时间序列中与米兰科维奇旋回无关的超高频和超低频信息，带通滤波能提取出目标频率信号（如岁差、斜率和偏心率旋回信号），将其与原始时间序列信号叠加在一起，以显示目标信号在时间（深度）域上的变化特征［3，107-108］。天文地层学研究中一般使用带通滤波的方法来提取符合地球轨道周期的滤波频率信号，并且将其与原始的深度域（时间域）数据序列叠置在一起，以突出目标周期信号（偏心率、斜率、岁差）在深度域（时间域）上的变化特征，以及这些目标周期信号之间的相位和振幅变化关系［103-105］。

例如，尺度值为108时，对应124.22 kyr的短偏心率周期。尺度值为331时，则对应405 kyr的长偏心率周期。在小波分析基础上分别提取出尺度值为108和331所对应的小波系数曲线（图6），用以代表东营凹陷樊页1井沙三下亚段的周期旋回曲线［74］，识别出研究层段存在的轨道周期。

图6 东营凹陷樊页1井沙三下亚段旋回周期曲线Fig.6 Cyclicity curve of the E s3L in Well Fanye 1 in the Dongying Sag

对于新生代的地层，滤波曲线可以直接与Laskar等理论目标曲线进行对比，进而通过天文调谐来建立天文地质年代标尺，即在确定了研究层段某一深度的绝对年龄后，将其作为年龄“锚点”［106］，进而可根据识别出的米兰科维奇旋回个数建立绝对天文年代标尺，计算得出地层持续时间及沉积速率等；若研究层段有生物地层、磁性地层或放射性同位素测年得到的顶底绝对年龄数据［18，106］，则可以确定出研究层段持续时间，同时检验米兰科维奇旋回个数的正确与否，建立高精度的天文年代标尺，也能更加确切地计算出某一地质事件发生的时间及持续时间。但对于中生代或者更老的古生代地层来说，只能通过天文调谐来建立“浮动”天文年代标尺［106］。

除天文调谐之外，还有其他的一些调谐方法，包括平均频谱误差分析（Averge Spetral Misfit，ASM）和滑动相关系数（Evolution Correlation Coeffient，ECOCO）等，它们都是基于沉积速率估算的调谐方法。ASM法是通过一系列综合测试估算受轨道驱动影响的地层间隔中的最佳沉积速率，并且评估输出的数据是否符合天文轨道信号，该方法的原理是利用假设检验和蒙特卡洛模拟来估算最优沉积速率，从而得出由数据的比例和相对应的由天文轨道信号频率加权后降噪的函数，最终将地层序列转化为时间序列［21，98］。ECOCO法是一个强大的工具，该方法是估算天文旋回时间序列分析的功率谱与替代参数序列之间的相关关系，将替代参数转换为一系列的时间序列用于估算沉积速率，同时还考虑了对沉积速率有贡献的天文轨道周期的数量，并使用蒙特卡洛模拟法对天文驱动信号的零假设进行检验，最终认为，最合适的沉积速率由高相关系数、多旋回信号及低零假设水平共同决定［22，98，102］。通过调谐，即结合放射性同位素测年、古生物或磁性地层提供的初始年龄估算平均沉积速率，建立高精度年代地层格架，并进一步使用ASM法或ECOCO法来校验沉积速率。根据建立的高精度年代地层格架，可以精确计算各种地质、气候和生物事件发生的时间，探讨不同天文轨道周期信号与一系列地质事件的耦合机制，为认识地球轨道力对气候和环境变化的影响做出贡献。

总之，天文旋回调谐的方法较多，而天文调谐方法是将运用时间序列方法得到的天文旋回结果与理论轨道周期和滤波方法结合起来完成的，该方法操作简单、直观且应用最为广泛，但是该方法的准确程度稍低，其前提是理论轨道周期曲线是正确的；ASM法是基于沉积速率估算的一种方法，它不受地质年代的限制，但是该方法对处理中的参数选取要求较高；ECOCO法的优势是有助于识别沉积速率变化和可能存在的沉积间断［18，98］，但是该方法处理过程较为繁琐，需要在固定的窗口长度内进行频谱分析，还要以一定的步长将频谱分析窗口滑动至时域内其他区域进行频谱分析［109］。

4 天文地层学研究流程

天文地层学的研究要按照沉积地层序列替代参数选取到时间序列分析方法优选，再到旋回结果检验的流程进行（图7），其中最为关键的就是获取用于旋回研究的沉积地层序列替代性指标，并选取最优的时间序列分析方法，从而保证旋回分析的准确性和可靠性。利用时间序列法开展米兰科维奇旋回研究的方法涉及数据的选取、预处理、频谱分析、时-频分析、天文调谐等一系列的研究过程，常用的时间序列分析方法有频谱分析、快速傅里叶变换、小波变换和滑动窗口频谱分析法等，通过这些方法检验沉积记录中的周期信号是否受到天文轨道驱动的影响，若存在轨道驱动，则需要将信号调谐到理论轨道参数曲线上，从而建立高精度年代地层格架。

图7 天文地层学研究的主要工作流程Fig.7 Workflow for astrostratigraphy study

5 存在的问题及下步研究方向

5.1 存在的问题

1）不同时-频分析方法的优选问题。目前旋回划分研究中，各种时-频分析方法均有涉及，但不同学者针对同一研究区的旋回划分结果稍有差异，其根本原因就在于所运用的时-频分析方法不同，而即使是运用同一时-频分析方法，不同学者在参数值的选取上也存在差异，也会导致旋回划分结果的异同。因此，同一研究区域不同时-频分析方法的对比研究及参数优选尤为重要，在旋回划分时要将多种时-频分析方法综合起来使用，从而优选最为准确的时-频分析方法及参数值。

2）不同替代参数的优选问题。多数学者在开展旋回划分工作时，所选用的替代参数多为测井数据和磁化率，以及少量元素地化分析数据。因资金及测试分析手段有限，成像测井数据及元素地化分析数据的使用受到限制。成像测井可以通过提取基质电导率分析旋回，元素地化分析数据则需要更广泛、更全面的使用。针对同一研究区，可以将不同替代参数的旋回划分结果进行对比，从而优选适合的替代性参数，保证旋回分析的准确性。

3）天文旋回划分结果的校准问题。旋回划分可运用不同替代参数采用不同时-频分析方法，但运用不同参数的不同方法、不同参数值的选取均会得出不同的旋回划分结果，导致旋回划分的多解性强，从而导致旋回划分结果的可靠性存在争议。主要存在两个方面的问题：①同一口井、同一层位不同学者划分出了不同的旋回；②不同的测井资料时-频分析方法与不同尺度天文周期、沉积单元之间的匹配关系研究还不够深入。因为天文轨道周期与不同级次的高频旋回之间有一定的对应关系［110-111］，而诸多学者对于同一口井、同一层位，运用时-频分析方法划分的不同尺度天文旋回与运用沉积学等方法划分的高频旋回结果存在差异，既包括旋回个数的差异也包括研究层段持续时间的差异［5，34-36，112］。但正常来说，天文轨道参数与高频旋回之间有一定的对应关系，故其同一层段计算得出的持续时间应该是一致或者极其接近的。因此，在天文旋回研究中，缺少一个“标杆”，而这个“标杆”就是在绝对年龄测年数据约束下利用沉积学原理和经典层序地层学理论对露头剖面和岩心开展研究所建立的旋回划分方案，有了这个“标杆”后，就可用于校准使用其他方法得出的旋回划分结果，也就可以优选替代参数、优选时频分析方法和参数值。

4）天文年代标尺建立的问题。在调研天文旋回研究现状时，发现针对同一研究对象建立的年代标尺不一致，分析其原因主要是在建立年代标尺时，各个学者针对同一层位的同一深度所选择的“锚点”年龄值不一致，没有统一的年龄值标准［36，113-115］。比如本文研究的实例区，大部分学者一般沿用了1992年梁鸿德或2002年姚益民公开发表著作中的年龄值［115-116］，其主要是通过古地磁和火山岩同位素测得的年龄值。众所周知，不同的测试机构和测试方法测试的同一深度地层的绝对年龄值是不一样的，且本身也存在一定幅度的误差。

5）地层分层数据（深度）不统一的问题。在旋回划分研究成果调研中，发现导致旋回划分结果不一致的另一个原因是地层分层数据的不统一，对同一口井、同一研究层段的顶-底深度使用不一致［5，17，19-20，111］。

5.2 下步研究方向

1）重视天文旋回分析技术手段和替代参数的优选研究。目前天文旋回分析的技术手段已经相对完善，替代参数也丰富多样，包括古生物参数、地球物理参数、地球化学参数和环境磁学参数等，通过对各类替代参数的处理，构建包含地层环境变化信息的时间序列。随着研究的不断深入和技术的进步，后期需要进一步加强对不同时-频分析方法、不同替代参数及不同参数值优选的研究，现阶段缺少各种方法和各种替代参数的对比研究，通过对比分析分别优选出适用于湖相细粒沉积泥页岩地层、碳酸盐岩地层和海、陆相碎屑岩地层的旋回划分方法和技术流程，按照沉积地层序列替代参数选取、时间序列分析方法优选、旋回结果检验及调谐的流程展开研究。

2）加强元素参数在天文旋回分析中的使用。目前使用较多的是单一元素分析数据，后期要综合考虑单一元素分析数据和元素比值的旋回划分效果，同时还要分别考虑表征古气候、古水深和古盐度等沉积环境各项参数的适用性，全面系统地将表征沉积环境要素的各项地化数据应用于旋回时间序列分析中。

3）重视建立用于旋回划分校准的“标杆”。时间序列分析方法划分天文旋回，由于时频分析方法、替代参数及参数值选取的差异会导致旋回划分结果的多解性，因此后期在开展旋回划分时需要建立一个与之对比的“标杆”用于校准旋回，也可用于方法、参数及参数的优选。以湖相细粒沉积旋回划分为例，首先从沉积学和经典层序地层学角度开展湖相细粒沉积地层高频层序划分，将划分结果与利用不同测井数据的不同时-频分析方法划分出的天文周期进行对比分析，最终优选出天文旋回划分的最优方法、参数和参数值。

4）加强天文轨道周期对沉积作用控制的研究。前人已经针对轨道周期对古生物、沉积环境条件、沉积物质组分的响应特征开展了不同程度的研究，但总体研究程度都比较低，特别是对沉积作用尤其是事件性沉积作用的控制研究更是相当薄弱。因此需要重视并加强轨道周期对沉积作用的控制研究，利用天文旋回的研究成果，建立千年—万年级的天文年代浮动标尺，计算地层沉积速率，分析沉积作用及沉积方式；在建立的高分辨率等时地层格架中探讨轨道周期与事件性沉积的耦合机制，为事件性沉积作用的研究提供新的思路和方法。

5）拓宽天文地层学在油气勘探开发中的应用领域。随着非常规泥页岩油气勘探与开发的不断深入，细粒沉积的研究越发受到重视，诸多研究也证实了中生代-新生代陆相地层较好地记录了天文地层周期。同时，针对厚层泥页岩天文地层学研究也如火如荼的展开，这既有助于更好地解决泥页岩油气勘探开发过程中的地层划分问题并提高地质分层的时间分辨率，也能加强天文地层学在深湖或深海页岩油气勘探与开发中的应用。