米兰科维奇旋回时间序列分析法研究进展

宋翠玉，吕大炜

山东科技大学地球科学与工程学院，山东青岛 266590

0 引言

周期性现象普遍存在于地球的各个角落。在地质记录中，周期性表现为韵律，是地层学研究最早发现的现象之一。地层学里把长而复杂的周期性节律称为旋回[1]。1988 年在意大利Perugia 召开的一次沉积学会议上，旋回地层学（Cyclostratigraphy）这一术语被首次公开使用[2]。对地层记录的（准）周期性旋回变化进行识别、描述、对比和成因解释，并应用于地质年代学以提高地层年代框架的精度和分辨率，实现高精度地层划分与对比，是旋回地层学的基本任务[3]。沉积地层中的韵律或旋回有外周期型和内周期型两类，前者由周期变化的外力（如气候、构造、海平面变化等）驱动，后者则产生于沉积系统内部，如一次沉积事件造成的沉积搬运等[4]。轨道周期是地球系统中十分常见的外周期。20世纪40年代，塞尔维亚机械工程师米兰科维奇提出了第四纪冰期成因的天文假说[5]。他认为地球自转和公转的三个重要的轨道参数：偏心率、黄赤交角或地轴倾斜度和岁差受其它天文因素的影响发生周期性变化，而这些变化引起了太阳辐射量的变化，地表北纬65 度附近接收到夏季太阳辐射量的变化是驱动第四纪冰期旋回的主要原因[5]。米兰科维奇理论的冰期旋回，至今是地球科学中唯一公认有定量基础的周期变化[4]。为了突出地球轨道参数变化对地层沉积作用的影响，徐道一等[6]在旋回地层学基础上，基于天文学的地球轨道三要素的周期性和准周期性变化提出了“天文地层学”这一概念，指出可通过功率谱分析和数字滤波等信号处理与运算进行地层划分与对比[7]。地层中由地球轨道驱动造成的旋回性记录称之为米兰科维奇旋回[8]，简称米氏旋回，是旋回地层学研究的重点。

目前，米氏旋回很大程度上解释了许多地质学领域的问题。越来越多研究成果表明，气候变化只是轨道力驱动的一个方面，从海水混合到岩浆作用，也都对轨道周期有所响应[4,9]。在不同地质历史时期的海相、陆相地层中，米氏旋回正不断地被揭示出来[8,10]。识别并研究地质记录中的米氏旋回信号，对于提高定年精度[11-15]、构建高分辨率年代地层格架[16-18]、估算地层剥蚀量[19-21]、研究古气候—古环境演化[22-26]等均具有重要意义。

米氏旋回的研究方法总体上可归为两类[8]：一类是岩性直观识别方法，通过直接观察岩性和岩相的变化，根据地层的组合特征、方式及级序结构来判别是否存在旋回信号[12,27]。该类方法要求研究者具有丰富的野外工作经验，能够准确地把握岩性的识别和读取，需要有雄厚的地质功底。但是随着研究深入，人们逐渐发现：在深湖、远洋、半远洋等沉积环境中，有些轨道参数的变化未必能引起明显的岩性或岩相变化，故直接识别法不仅难度大，还容易于造成旋回信息的遗漏[8]。第二类是时间序列分析法，利用岩性岩相特征的数字化值（如特殊岩石的颜色变化等）或者地层的古气候替代指标等（如地层的测井曲线、地层中某种化学元素的含量变化等）构建包含地层环境信息的时间（深度）序列，进而利用数学变换对这些数据序列进行定量分析[2,28]。这种方法需要借助大量的先进技术手段获取数据序列，比如自然伽马（Gamma Ray，GR）测 井 设 备[29-30]、XRF（X Ray Fluorescence）元素扫描仪[31-34]、颜色识别仪[33,35]、磁化率测试设备[32-33]等，并且也要求相关人员具有较强的地质学基础。应该说，时间序列分析法是目前最常用的旋回地层分析法，它极大地推进了旋回地层学的发展。

综上所述，米兰科维奇理论已经大大促进了地层学研究的发展，但是随着数学方法的不断革新，米氏旋回研究的方法学也相应随之改进和提高，为此，本文在近几年的米氏旋回方法学研究基础上，针对时间序列分析法的主要步骤，将分别从地层序列的数据类型、天文检验及天文调谐等三方面，简述现有方法的基本原理，总结其优势和局限性，讨论存在的问题，并展望时间序列分析法的发展方向，为米氏旋回方法学研究提供一定的基础。

1 基于时间序列分析的米氏旋回研究方法

利用时间序列分析法进行米氏旋回研究的工作方法较为成熟。主要涉及地层序列的数据类型选取、数据预处理、频谱分析、天文驱动检验和天文调谐等环节和内容，其一般流程如图1所示。

图1 基于时间序列分析的米氏旋回分析流程图Fig.1 Flow chart of Milankovitch cycle studyusing time series analysis

不同的地层数据序列蕴含着不同的地质信息，对米氏旋回信号的记录亦有差异，因此，选择合适的数据序列是米氏旋回分析的首要且关键步骤。数据预处理是指频谱分析之前对原始数据进行的一系列处理。主要包括重采样、剔除异常值、去趋势化和预白化等，其目的是尽量消除地层数据中的各种环境“噪声”，使结果更容易解释[8,32]。频谱分析将地层序列的信号强度按频率顺序展开，使其成为频率的函数，进而识别出信号中（准）周期性的成分[8]。它是检验地层中是否包含米氏旋回记录的前序处理。通过频谱分析，可以获得数据中蕴含的一系列主周期信号，天文驱动检验就是检验这些周期信号中是否包含地球天文轨道驱动而成的米兰科维奇周期信号。如果地层中存在米氏旋回信号，则可以通过天文调谐将从数据序列提取出的旋回记录对比到岁差、斜率和（或）偏心率的天文目标曲线上，使深度域的地层沉积序列变换到时间域，以建立天文年代标尺或浮动天文年代标尺，用于定年或确定地层∕地质事件的持续时间。

目前，已有多种米氏旋回时间序列分析方法，如ASM（Average Spectral Misfit）[36]、eASM（evolutive Average Spectral Misfit）[37]、eCOCO（evolutionary Correlation Coefficient）[38]、TimeOpt（Time scale Optimization）[39]和eTimeOpt（evolutive Time scale Optimization）[40]等。

1.1 地层序列的数据类型

目前，应用于米氏旋回研究的数据序列可大致归为四类（表1）。第一类为直观反映岩石物理性质的数据，例如颜色[33,35]、粒度、相变与层理韵律[33,41]等；第二类为地球化学数据，如Fe 含量[38]、CaCO3质量分数[36]、碳同位素[42-44]、主量元素含量[24]、微量元素含量[24]、元素比值[31,45]等；第三类为地球物理数据，如磁化率[22,24,32-33,45-46]、非磁滞剩磁（Anhysteretic Remanent Magnetization，ARM）[32-33]、GR（Gamma Ray）测 井 数据[23,29,33,35,46-49]、电阻率测井数据[46]等。最后一类为古生物相关的数据，例如生物丰度、生物灭绝速率[32,50]等。

表1 用于米氏旋回研究的数据序列类型列表Table 1 Data series types used in Milankovitch cycle study

由于对古环境、古气候变化反映敏感且测试较方便，岩石磁性特征和自然伽马测井数据等地球物理指标被广泛应用于旋回地层学分析中。磁性矿物是气候和环境变化信息的重要载体。表征物质磁学特征的磁化率则反映了岩石中所含磁性矿物的种类、含量和粒径等信息，是良好的古气候及古环境替代指标[22]。GR 测井测量沉积物中伽马射线的强度，可以敏感反映沉积物中泥质含量，进而反映古气候和古环境的微小变化，是米氏旋回分析的理想数据[29]。一般而言，GR 高值对应富含黏土沉积物，GR低值与砂岩、富含碳酸盐的沉积物有关[52]。自然伽马能谱测井除了获取得出的GR值外，还记录钾（K）、铀（U）和钍（Th）等元素的含量[30,33,52]。

不同数据指标对天文轨道驱动和非天文噪声的响应是有差异的[33,35,53]。例如，在海相沉积地层中，ARM 和Th∕U 可反映早三叠世的内陆风化作用，GR、K、U、Th、磁化率和非碳酸盐组分等指示陆源补给，而岩石颜色数字化后生成的CIE 分量L*、a*以及岩性序列（lithologic rank）则分别指示生产力、氧化还原状态及相对海平面[33]。此外，有些数据指标受轨道力驱动，而有些指标可能主要受非轨道噪声的影响[35]。因此，基于单一沉积序列构建的天文年代标尺具有不确定性。

为了减弱单一数据序列中环境“噪声”造成的不确定性，同时综合利用不同数据序列中蕴含的地质信息，已有大量研究利用多种数据指标进行旋回地层学研究。这些研究中，往往首先将多个数据指标逐一进行时间序列分析，然后联合解读它们对米氏旋回的指示作用，进而进行地层划分或者古环境研究[30-32,54]。从松辽盆地晚白垩世嫩江组的地球化学元素比值Rb∕Sr、K∕Ti、Ti∕Al 和Zr∕Rb 中识别的米氏旋回信号，揭示了嫩江组沉积时期气候由湿润向半干旱、半湿润过渡的古气候响应机制[31]。GR 和SP（Spontaneous Potential）测井曲线的频谱分析结果表明，北黄海东部坳陷始新统地层很好地保存了米氏旋回，据此可进一步分析其沉积环境并进行地层精细划分[54]。天津蓟县剖面前寒武系洪水庄组—铁岭组地层中，磁化率和伽马能谱数据指标中均记录了完整的米氏旋回[30]。房强[32]利用磁化率、ARM 和地球化学等多个指标对上寺和渡口剖面茅口组进行旋回地层学分析，提取了米氏旋回在晚古生代冰期末期的记录，并讨论了其古环境响应。上述研究尽管采用了多种数据序列进行旋回分析，但它们均作为单一输入独立进行了时间序列分析，并未真正实现多种地层信息的有效整合。近年来，融合多种数据类型的米氏旋回分析方法正引起关注。根据多种数据指标的独立模型来构建综合的年代模型，被证实是估计年龄模型不确定性的有效方法[35]。此外，对多种数据指标先做主成分变换等处理，继而对某一个或者某几个分量进行进行米氏旋回分析，也有助于识别完整的米兰科维奇周期信号[46]。

1.2 天文驱动检验方法进展

米氏旋回研究解决的核心问题之一，就是检验地层记录中是否存在米氏周期信号，即检验沉积记录中的周期信号是否由地球天文轨道驱动[40,55]。周期的倒数是频率，因此，通过频谱分析把深度域的沉积序列观测值变换到频率域，是检验天文驱动信号的关键步骤。

1.2.1 频谱分析方法研究现状

旋回地层分析中，常用的频谱分析方法有快速傅里叶变换（Fast Fourier Transform，FFT）、周期图法（Periodogram）、Welch 法、多窗谱分析法（Multi-taper Method，MTM）和 最 大 熵 谱 法（Maximum Entropy Method，MEM）等。不同的频谱分析方法所基于的算法不同，优缺点也各不相同，但计算结果基本一致[8]。

傅里叶变换是最经典的频谱分析方法之一。它把在时间域里看起来很复杂的多频率叠加信号分解成不同频率的正弦分量，将空间域变换到频率域[56]。由于实际数据的长度有限，通常采用离散傅里叶变换（Discrete Fourier Transform，DFT）对数据进行谱分析，即对有限长度的离散时间序列进行傅里叶变换。快速傅里叶变换（FFT）是DFT的快速算法，将其运算量减少了几个数量级。

周期图法是一种算法简便且具代表性的功率谱估计方法，其原理是对观测到的数据直接进行傅里叶变换，然后取模的平方获得功率谱。周期图是对功率谱密度（Power Spectral Density，PSD）的有偏估计。在信噪比大、数据够长的时候，周期图是有用的谱估计器。但在数据序列不够长时，周期图法的分辨率往往受限。Liet al.[38]提出的eCOCO 和Meyers[39]提出的TimeOpt 方法中，频谱分析法均采用周期图法。

Welch法是对周期图法的一种改进，它先将数据序列划分为不同的段（可以有重叠），然后对每段进行加窗处理，再分别计算周期图后取平均。该方法能改善方差性能和分辨率，是较常用的谱估计方法之一[2]。

MTM 是由Thompson 提出的一种多窗谱分析和信号重建方法，能提供更优的PSD估计[57]。它使用一组正交、离散的扁平类球体序列（Discrete Prolate Spheroidal Sequences，DPSS，也叫做Slepian 序列）获得最优滤波器，进而计算功率谱估计值。通过调整MTM 中的时间—带宽参数，可以实现估计方差和分辨率之间的平衡。在非常短且有噪声的时间序列中，MTM具有较高的频率分辨率，并且为每个频率提供了统计显著性检验（F-方差比检验）。此外，这种显著性检验与周期的振幅无关，因此，该方法不仅能够识别显著性水平高的低振幅峰值，还能指示哪些高振幅峰值可能不具有统计显著性。旋回地层学研究中，MTM是应用较广的频谱分析方法之一[30,36,58-59]。

MEM 的原理是根据已知数据信息，在不进行任何新的假设的情况下，按信息熵最大准则预测未知延迟离散时间上的相关函数[60]。MEM法的优势在于解决含有低噪声的间隔紧密的正弦短数据信号。与其他谱分析方法（如周期图法）相比较，MEM 具有不受取样长度限制等优点[41,61]。但受噪声、正弦信号初始相位等因素的影响，MEM 谱估计会产生频率偏移。

图2列出了某剖面3 130 m至3 441 m深度GR序列的5种频谱估计结果，采样间隔为5 cm，共有6 219个数据点。可以看出，FFT 与Welch、MTM 与周期图法获得的功率谱形态分别有一定相似性，而MEM分辨率低，谱非常光滑，难以筛选米氏周期信号。与FFT 相比，Welch 获得的谱更光滑，分辨率降低了。MTM 和周期图法获得的谱的分辨率均比其它方法高，但二者功率谱极值对应的频率位置具有差异性。MTM的计算复杂度较高，是计算DPSS的代价。

图2 基于不同方法的某剖面GR 数据频谱分析结果Fig.2 Spectral analysis results of GR series in an unknown section using five methods

频谱分析的作用就是评估能量较高的非随机频段是否与米氏旋回的频率吻合。对沉积序列数据做频谱分析后，功率谱峰值对应的频率即为米氏旋回信号的候选值。然而，这些周期（频率的倒数）信号既包含真实的米氏旋回信号，也包括一些无关的、随机出现的“噪声”。这些噪声的能量可能较低，但或许会出现在所有频率中。研究表明，沉积记录中的噪声在低频谱段能量较高，随频率增大能量逐渐降低，与光类似，故称为“红噪”（Red Noise）[62-63]。红噪的出现表明随机振荡变化是气候的一个普遍特征。即使没有气候驱动参数（如日照量）的变化，气候系统中的各组成部分及其相互作用也可以产生强烈的低频振荡[40]。因此，需要对噪声进行估计，并通过显著性检验将这种随机振荡与天文驱动信号区分出来。

常采用一阶自回归模型（Auto Regressive Lag-1，AR1）估计沉积序列中红噪声，生成红噪声理论曲线并进行显著性检验，将设定置信水平（90%、95%或99%）之下的峰值对应的频率剔除。自回归Lag-1（AR1）随机噪声模型是天文年代学研究中评价旋回显著性的最常用方法之一，它是一种直接基于自身数据的估计模型[2]。图3为Chaohu剖面不同深度GR数据的MTM 频谱分析结果，噪声估计均基于鲁棒的红噪模型[51]。对于10~74 m 的地层而言，置信水平95%以上的峰值对应波长有10 m、2.3 m、1 m、0.8 m、0.7 m和0.5 m。

1.2.2 天文检验方法研究现状

在噪声压制后的功率谱中，判别峰值对应的频段组合是否由天文轨道驱动的方法主要有三种，即比值法、调幅分析法和假设检验法。

一般而言，通过频谱分析获取的旋回周期之比与天文轨道参数的周期之比一致，即可认为沉积记录中包含米氏旋回。例如，地球天文轨道参数长偏心率、短偏心率、斜率和岁差的周期之比约为20∶5∶2∶1，如果频谱分析中提取的四个周期（频率取倒数）之比与此接近，则可认为存在米氏旋回。该方法称为比值法，是目前绝大多数旋回地层研究采用的天文检验方法[29,51,55,59,64-65]。图3 所示10~74 m 的地层，周期波长～10 m、～2.3 m、0.7~1.1 m和0.5 m的比值与米兰科维奇频率405 kyr、100 kyr、33 kyr和20 kyr的比值接近，故可以判断米氏旋回存在。然而，由于受到噪声、频谱分析误差等影响，满足比例的频率组合可能不唯一。

图3 Chaohu 剖面GR 数据的2π MTM 的频谱分析结果[51]E:长偏心率；e:短偏心率；o：斜率；p：岁差Fig.3 2π MTM power spectra of GR series from intervalsin the Chaohu Section[51]E:long eccentricity;e:short eccentricity;o:obliquity;p:precession

长周期旋回对短周期旋回会有明显的天文调制，例如短偏心率对岁差进行天文调制，长偏心率对短偏心率天文调制[66]。所以，调幅分析也是一种判断观察到的旋回是否受天文轨道控制的有效方法[67-68]。然而，沉积序列数据中类似偏心率的振幅变化可能通过调谐和数据处理等过程人为地引入[67]。因此，振幅调制方法不能作为一种独立的方法来测试天文调谐时间尺度的准确性。

检验天文驱动的另一类常用方法为统计学中的零假设检验法。该类方法先提出零假设，即假设筛选出来的频率都是与天文驱动无关的随机噪声，然后借助统计检验拒绝该假设，即可说明序列中含有米氏旋回信号。Meyerset al.[36]提出的ASM 方法，Liet al.[38]提出的eCOCO 方法中，均采用零假设检验天文驱动信号的存在。但当轨道信号因沉积速率的变化而严重失真时，ASM 将无法以高置信水平拒绝零假设[37]。

1.3 天文调谐方法研究进展

在地层中识别出完整的米氏旋回信号后，就可以将信号调谐到理论的天文曲线上[69]。通过年代校准即可建立高分辨率地层层序格架。在年龄控制较差和缺少天文目标曲线的古生代和中生代，调谐只能建立浮动天文年代标尺，来精确地计算地层或地质事件的持续时间。天文调谐方法可归纳为三类（表2），即最小调谐法、统计法和反演法等[38]。

表2 主要天文调谐方法一览表Table 2 Main astronomical tuning methods

1.3.1 最小调谐法

最小调谐法（Minimal tuning）[70]使用最小数量的校准点将深度域序列曲线与时间域目标天文曲线联系起来完成校正，是最经典、也是应用最广泛的调谐方法[51,55]。该方法先将沉积序列数据调谐到单个天文频率（如偏心率）上，然后对调谐后的结果再进行频谱分析，如果也能检验出其它的天文轨道信号（如斜率、岁差）则调谐成功[32]。实际操作中，为了更准确地定位校准点，常先将沉积序列曲线做带通滤波（Band Pass Filter）处理。这种调谐方法直观简单，但存在局限性：首先，该方法假定目标曲线正确，没有利用统计法进行检验[55]；其次，选取的目标曲线年代与真实数据年代难以完全对应，从而产生相位差。

在缺少天文目标曲线的古生代，频率域的最小调谐往往获得理想效果。首先，利用EHA（Evolutive Harmonic Analysis）[71]或改进的FFT[24,48]在滑动窗口上对地层数据序列做频谱分析，获得频谱图，该图表征频率信号在深度上的偏移规律。然后，在频谱图中追踪单个较连续的旋回信号（如偏心率周期），获得沉积速率曲线。最后，根据沉积速率曲线将沉积记录调谐到时间域。图4a为某剖面GR数据经EHA处理后的频谱图，滑动窗口长度为60 m，图中白点表示追踪到的长偏心率周期（405 kyr），据此计算出沉积速率曲线（图4b）。由于缺少年龄控制，经调谐建立了该剖面的浮动年代标尺（图4c）。

图4 频率域最小调谐法示意图（a）EHA频谱；（b）沉积速率曲线；（c）天文调谐结果Fig.4 Sketch map of Frequency Domain Minimal tuning(a)EHA spectrum;(b)sedimentation rate curve;and(c)tuning result

理论上，对于深度域的地层序列，如果获得其对应的沉积速率曲线，则可建立浮动年代标尺，根据已知地质年代锚点，即可确定沉积序列的绝对年龄。现有的米氏旋回研究中，基于沉积速率定量估算的调谐方法主要有统计法和反演法等。

图5 为樊页1井磁化率指标的ASM分析结果[17]，沙三下亚段最优沉积速率为10.801 cm∕kyr，零假设检验的显著性水平低于0.1%，表明沉积速率为10.801 cm∕

1.3.2 统计法

统计法基于一系列可能的沉积速率，利用假设检验和蒙特卡洛模拟来估算最优沉积速率或追踪沉积速率的变化。这些方法不需要精确的年代限制，也能更好地保持数据的原始相位。ASM[36]及其改进方法eASM[37]、TimeOpt[39]及eCOCO[38]等方法都是统计法。

给定一系列可能的沉积速率，ASM 对天文轨道驱动的信号进行零假设检验，同时通过蒙特卡洛模拟来评价识别出来的优势频率与目标轨道频率的拟合程度，见式（1），根据ASM 值极小及低零假设水平确定最优沉积速率的客观估计值。

式中：f为功率谱峰值对应的频率（单位：cycles∕m）；s为待测试的沉积速率（单位：m∕kyr）；fpred为频率的理论值（周期取倒数，单位：cycles∕kyr）；ΔfR是最小分辨率带宽。kyr 的置信度超过99.9%。然而需要注意的是，ASM方法仅是量化匹配的最优沉积速率，并不是平均沉积速率，也不是计算某一段的精确沉积速率[17]。

eASM（evolutive ASM）[37]对ASM做了四点改进：①对式（1）中的αk做了调整（见式（2）），以更好地考虑谱峰位置的不确定性；②对测试的沉积速率做对数运算，以更好地评估低沉积速率时经常出现的ASM快速变化；③引入了滑动窗口，可以适应地层区间沉积速率的变化；④考虑理论轨道目标频率（fpred）及其不确定性。

Meyers[39]提出的TimeOpt也是一种统计法。该方法基于一系列可能的沉积速率：①利用概率线性回归法，根据公式（3）重建岁差的偏心率调幅，并计算重建结果与实际调幅结果的相关系数r2envelope；

yenvelope=Xeβe+ε（3）

式中：yenvelope是时间校准数据序列的岁差振幅包络线，Xe是代表偏心率周期的正弦和余弦预测项矩阵，βe是每个预测值的回归系数向量，ε是误差项向量；

②再次利用概率线性回归法，根据式（4）重建包含岁差和偏心率理论值的功率谱，并计算重建结果与年代校准数据之间的相关系数r2spectral；

ydata=Xepβep+ε（4）

式中：ydata是时间校准数据序列，Xep是代表偏心率和岁差周期的正弦和余弦预测项矩阵，βep是每个预测值的回归系数向量，ε是误差项向量；

③乘积r2opt=r2evenlope r2spectral极大值对应的沉积速率为最优。

图6 为ODP Site 1 262 颜色值a*（绿∕红）的TimeOpt分析结果[39]。由图6e可以看出，包络回归模型（公式（3））确定的最优沉积速率为1.28 cm∕kyr，谱功率回归模型（公式（4））和联合模型（乘积）确定的最优沉积速率均为1.33 cm∕kyr，此时，最大r2opt为0.212（图6f）。经过2 000 次AR1 蒙特卡洛模拟，得出在1.33 cm∕kyr 时观察到的最大r2opt的p 值为0.005，表明我们可以在99.5%的置信水平下拒绝零假设（图6g）。在该地层中，TimeOpt得出的沉积速率为1.33 cm∕kyr，与先前基于天文年代学、磁性地层学和生物地层学的估计值（如1.3 cm∕kyr，1.2 cm∕kyr，1.9 cm∕kyr等）相似，证实了该方法在沉积速率估算中的有效性。

图6 ODP Site 1 262 颜色值a*（绿∕红）的TimeOpt 结果（据Meyers[39]修改）（a）1 262颜色数据；（b）颜色数据的周期图，TimeOpt确定的沉积速率为1.33 cm∕kyr（黑线、深灰线分别为线性、对数光谱），灰色阴影区为用于评估岁差振幅包络的部分，灰色虚线表示偏心率和岁差的目标周期；（c）带通岁差信号（黑色）与Hilbert变换获得的数据振幅包络（灰色）对比；（d）数据振幅包络（灰色）与TimeOpt重建的偏心率模型（黑色；据公式（3））对比；（e）每个待测沉积速率下振幅包络拟合的平方皮尔逊相关系数（r2envelope；灰色线）与谱功率拟合的相关系数（r2spectral；黑色线）；（f）每个待测沉积速率下包络线和谱功率的联合拟合（r2opt）；（g）基于AR1的2 000次蒙特卡罗模拟结果（ρAR1=0.907），用于评估最大r2op（t0.212）的显著性；（h）数据振幅包络线与TimeOpt重建的偏心率模型在平面“d”中的交会图；灰色虚线为1∶1线Fig.6 TimeOpt analysis of color data (a*, representing the red∕green ratio) from ODP Site 1 262(modified from Meyers[39])(a) The 1262 color data; (b) Periodogram for the Site 1262 color data, given the TimeOpt derived sedimentation rate of 1.33 cm∕kyr (black line = linear spectrum; dark gray line=log spectrum).Gray shaded region indicates the portion of the spectrum bandpassed for evaluation of the precession amplitude envelope.Dashed gray lines indicate the eccentricity and precession target periods; (c) Comparison of the band-passed precession signal (black) and the data amplitude envelope (gray) determined via Hilbert transform;(d)Comparison of the data amplitude envelope (gray)and the TimeOpt-reconstructed eccentricity model (black;derived using equation(3));(e)Squared Pearson correlation coefficient for the amplitude envelope fit (r2envelope;gray line)and the spectral power fit(r2spectral;black line)at each evaluated sedimentation rate;(f)Combined envelope and spectral power fit (r2opt) at each evaluated sedimentation rate; (g) Summary of 2 000 Monte Carlo simulations with AR1 surrogates (ρAR1= 0.907), used to evaluate the significance of the maximum observedr2optof 0.212; and (h) Cross plot of the data amplitude envelope and the TimeOpt-reconstructed eccentricity model in panel“d”; dashed gray line is the 1∶1 line

与ASM类似，TimeOpt生成的也是恒定沉积速率模型，无法获取沉积速率的变化信息。eTimeOpt（evolutive TimeOpt）引入滑动窗口，在地层记录上顺序移动窗口进行TimeOpt 分析，可以追踪由盆地演化、相对海平面变化等引起的沉积速率的渐进长期变化[40]。

eCOCO是2018年Liet al.[38]提出的一种基于演化的相关系数的米氏旋回分析方法。对于显生宙古气候序列而言，eCOCO 能够同时评价沉积速率和天文驱动。它基于一系列待“测试”的沉积速率，1）将古气候替代成分转换为时间序列，然后计算该时间序列的功率谱与天文解决方案功率谱之间的相关系数；2）记录对相关系数有贡献的天文参数的数量；3）使用蒙特卡洛模拟法对天文驱动信号的零假设进行检验。最合适的沉积速率由高相关系数、多旋回信号及低零假设水平共同决定。由于eCOCO 采用了滑动窗口技术，可以跟踪替代成分沿地层序列在沉积速率上的变化。图7为晚三叠世Newark深度序列（0~2 000 m）的eCOCO分析结果[38]，其中暗红色（极值）交集决定了特定深度的最优沉积速率（图7b~d）。

图7 晚三叠世Newark 深度序列（0~2 000 m）的eCOCO 分析结果[38]（a）深度序列的高斯低通滤波结果（截止频率为0.01 m-1），金色水平条为～270 m旋回；（b）eCOCO相关系数图，黑线和绿线分别为Kentet al.[72]和Liet al.[73]的研究结果；（c）零假设检验H0显著性水平图；（d）有贡献的天文参数数量变化图Fig.7 eCOCO analysis of the Late Triassic Newark depth-rank series (0~2 000 m)[38](a)Gaussian low-pass filter output of the depth rank series(cutoff frequency is 0.01 m-1).Gold horizontal bars indicate ～270 m cycles;(b)Evolutionary correlation coefficient map shown with published sedimentation rate curves by Kent et al.(black line)[72]and Li et al.(green line)[73];(c)Evolutionary H0significance level map;and(d)An evolutionary map of the number of contributing astronomical parameters

1.3.3 反演法

反演法广泛应用于地球物理的各个领域，其本质为基于一定的数学方法将观测数据映射到相应的地球物理模型[74]。例如，通过地表观测的地震波数据，得到地球的区域速度结构，或通过观测的重力数据，获取地壳的密度分布等信息。

Malinvernoet al.[44]将轨道调谐视作反演问题，提出了基于贝叶斯反演的天文调谐法。该方法将变化的沉积速率及相应的地层深度作为参数建立年代模型（Age Model），首先应用贝叶斯公式定义轨道周期信号驱动的沉积速率变化的概率分布函数（而非单一沉积速率）。然后，通过马尔科夫链蒙特卡洛法（Markov Chain Monte Carlo）对概率分布进行采样，量化沉积速率的不确定性，这些不确定性来源于调谐周期的不确定性，以及与轨道周期无关的沉积信号成分。最优沉积速率选取似然度（likelihood）较大的沉积速率。

eASM、eCOCO 等统计法在地层记录中按顺序移动窗口，有助于追踪地层中变化的沉积速率信息，但无法获取窗口中发生的沉积速率的显著变化。因此，Meyers[40]在TimeOpt 基础上，引入沉积模板（sedimentation template）（图8）构建复杂的沉积模型，提出扩展的TimeOpt 反演法（对应Astrochron 工具包中的timeOptTemplate 函数）。在该优化反演过程中，初始沉积模板与天文目标和数据动态交互调整，最终重建沉积速率曲线。

图8 扩展后TimeOpt 中的沉积模板示例（阶跃变化模板）[40]Fig.8 Step change sedimentation template of the extended TimeOpt[40]

2018年，Meyerset al.[75]将TimeOpt与贝叶斯反演法融合，提出了TimeOptMCMC（TimeOpt Markov Chain Monte Carlo）法。该方法定量地将天文学理论与地质数据联系起来，旨在克服二者的局限性。TimeOptMCMC 的核心有三部分：1）TimeOpt 法，考虑时间尺度的不确定性，并利用天文信号的多个属性来提高统计可靠性；2）天文学理论，将岁差、轨道偏心率周期与太阳系、地—月演化的基本频率联系起来；3）贝叶斯马尔可夫链蒙特卡罗方法，探讨数据、模型空间及不确定性。

2 存在问题

经过众多研究者几十年的努力，米氏旋回的研究方法不断发展与优化，基于时间序列分析的新方法也不断涌现。天文解决方案精度的不断提高，也为米氏旋回研究提供了重要基础。然而，纵观现有研究，米氏旋回时间序列分析法尚存在如下问题：

（1）数据序列中存在的噪声影响天文检验和调谐的精度。地层数据序列既记录了轨道驱动的环境变化信息，也包含了一些“噪声”。这些噪声可能来源于地层的古沉积环境，如沉积速率变化、生物扰动、滑塌作用和后期成岩作用等，也可能是在数据准备过程中引入的，如采样处理、仪器噪声等。在数据预处理环节进行去趋势化及预白化等处理、使用对噪声鲁棒的频谱分析方法、对频谱分析结果进行红噪声抑制等做法均会在一定程度上去除数据中的噪声，但频谱分析结果中仍可能产生大量的谱峰[68]，既包含轨道驱动的周期信号，也包含非轨道驱动的周期信号，给后续处理和研究都带来很大干扰。

（2）新型的定量化的天文调谐方法仍需要进一步验证及改进。米氏旋回研究中，绝大多数天文调谐方法采用传统的最小调谐法[70,76-77]，其次是借助EHA[31]或改进的FFT 等频率域最小调谐法[24,78]，而统计法和反演法应用相对较少[42,78]。传统调谐法直观性好，但通常需要交互操作，调谐过程中滤波主频段的选取、极小值∕极大值位置、目标天文曲线时间段选择等难免会受处理者主观因素的影响，不同处理者或同一处理者在不同时间的调谐结果都可能有所不同，导致研究的可重复性较差。统计法和反演法基于假设检验和概率等数学原理，能评价各种噪声带来的天文信号的不确定性，为定量、客观的沉积速率估算提供了重要支持。现有方法不仅能够获取地层序列的最优沉积速率（如ASM，TimeOpt 等），还能够追踪沉积速率的变化信息（如eASM，eCOCO，eTimeOpt）。然而，由于应用相对较少，这些新型方法在不同数据、不同条件下的天文调谐精度还有待于进一步验证。此外，对于统计法和反演法而言，参数调整中存在不确定性、缺少中生代和古生代精确的天文目标曲线等因素，也约束着这些方法的适用性及准确性。

3 研究展望

基于时间序列分析的米氏旋回研究方法仍处在发展阶段，目前已经出现了一些较为成熟的软件或工 具 包，如AnalySeries、Astrochron (R package)[40]、Acycle[79]等。地球科学信息化进程的不断推进，对研究方法智能化水平和可靠性等方面的要求更高。未来挑战与机遇并存，作者认为以下两方面值得开展进一步研究工作。

（1）沉积序列数据的选择与信息优化。数据是时间序列分析的基础，其信息有效性对分析结果起着关键作用。不同的数据序列对地层特征及环境变化的记录各有侧重，因此，使用单一数据类型进行旋回地层分析可能存在局限性[33]。对同一地质剖面而言，利用不同数据序列提取的米氏旋回频率组合也可能具有多解性[46]。为了减小数据类型差异性带来的不确定性，需要在以下两方面继续开展工作。

第一，需要建立定量的评价标准，在数据可获取的前提下，选择更适用于特定环境下旋回地层分析的数据序列。信噪比、层次聚类、能量谱分解等已被分别用于评价单序列的信息质量、多沉积序列间的关系及不同数据序列对外来气候驱动的敏感程度等[33]。形成稳定、可靠的评价体系仍需深入工作。

第二，当同一地质剖面有多个数据序列可用时，如何综合使用多个数据序列的优化信息，亦值得研究。在数据处理领域，比值运算、主成分变换、独立成分变换等是常用的特征提取方法。不同沉积数据序列可以视作描述地层环境的多个特征，因此，在沉积理论的指导下，可借助上述特征提取方法获取地层环境的优化信息，作为米氏旋回研究的数据基础。

（2）基于新方法的沉积速率客观、量化估算。准确获取沉积速率曲线是米氏旋回定量化研究的关键环节。天文调谐的统计法和反演法都是基于沉积速率估算的时间序列分析方法。统计检验和蒙特卡洛模拟等数学手段已经成功应用于天文检验和调谐中。不同估算方法的原理存在较大差别，但也不乏关联性。面对逐渐多样化的时间序列分析法，既需要有效整合现有方法，形成理论更严谨、性能更可靠、客观性更强的旋回地层学分析方法，也需要引入新方法。

基于一系列待测试的沉积速率或初始速率模板将沉积序列转换为时间序列，然后根据该时间序列与目标天文年代曲线间的拟合度或相似性确定最优速率（曲线），是沉积速率估算的最核心过程。现有方法，如eCOCO、TimeOpt中，采用的衡量变量间相似性的指标主要为Person相关系数和Spearman相关系数等。结合地层序列数据的分布规律，亦可引入并验证其它评价匹配度的指标。此外，随着全球范围内旋回地层学研究的不断完善和资料共享，相关数据和研究成果将继续积累，势必为未知地层的研究提供足量的学习样本，机器学习有望在米氏旋回研究领域发挥重要作用。对于暂时缺少精确天文目标曲线的中生代和古生代而言，海量已知地层序列信息的出现，无疑为该领域旋回地层学的深入研究提供了新契机。

4 结语

旋回地层学的发展对理解和解决地球科学领域的众多科学问题做出了重要贡献。时间序列分析法促进了米氏旋回分析向客观性、定量化方向发展，为地层学的进一步发展提供了可靠依据。然而，沉积序列数据中存在的“非天文”噪声影响着天文检验和调谐的精度，新型的天文检验和调谐方法仍需要进一步验证及改进。随着地质学理论和测试技术的不断革新，时间序列分析法仍将发挥其地质研究的辅助作用。未来研究中，如何综合利用多种数据类型蕴含的地质信息、如何有效整合现有方法、以及如何将机器学习等思想引入米氏旋回研究等问题，值得科技工作者继续开展研究。

致谢 论文评审过程中，收获了匿名评审专家详尽、宝贵的修改建议；修改过程中，得到了编辑部工作人员耐心、细致的指导，在此一并致以真诚的谢意！