米蘭科維奇旋回時間序列分析法研究進展

宋翠玉，呂大煒

山東科技大學地球科學與工程學院，山東青島 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個數據點?？梢钥闯觯現FT 與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

理論上，對于深度域的地層序列，如果獲得其對應的沉積速率曲線，則可建立浮動年代標尺，根據已知地質年代錨點，即可確定沉積序列的絕對年齡?，F有的米氏旋回研究中，基于沉積速率定量估算的調諧方法主要有統計法和反演法等。

圖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）基于新方法的沉積速率客觀、量化估算。準確獲取沉積速率曲線是米氏旋回定量化研究的關鍵環節。天文調諧的統計法和反演法都是基于沉積速率估算的時間序列分析方法。統計檢驗和蒙特卡洛模擬等數學手段已經成功應用于天文檢驗和調諧中。不同估算方法的原理存在較大差別，但也不乏關聯性。面對逐漸多樣化的時間序列分析法，既需要有效整合現有方法，形成理論更嚴謹、性能更可靠、客觀性更強的旋回地層學分析方法，也需要引入新方法。

基于一系列待測試的沉積速率或初始速率模板將沉積序列轉換為時間序列，然后根據該時間序列與目標天文年代曲線間的擬合度或相似性確定最優速率（曲線），是沉積速率估算的最核心過程?，F有方法，如eCOCO、TimeOpt中，采用的衡量變量間相似性的指標主要為Person相關系數和Spearman相關系數等。結合地層序列數據的分布規律，亦可引入并驗證其它評價匹配度的指標。此外，隨著全球范圍內旋回地層學研究的不斷完善和資料共享，相關數據和研究成果將繼續積累，勢必為未知地層的研究提供足量的學習樣本，機器學習有望在米氏旋回研究領域發揮重要作用。對于暫時缺少精確天文目標曲線的中生代和古生代而言，海量已知地層序列信息的出現，無疑為該領域旋回地層學的深入研究提供了新契機。

4 結語

旋回地層學的發展對理解和解決地球科學領域的眾多科學問題做出了重要貢獻。時間序列分析法促進了米氏旋回分析向客觀性、定量化方向發展，為地層學的進一步發展提供了可靠依據。然而，沉積序列數據中存在的“非天文”噪聲影響著天文檢驗和調諧的精度，新型的天文檢驗和調諧方法仍需要進一步驗證及改進。隨著地質學理論和測試技術的不斷革新，時間序列分析法仍將發揮其地質研究的輔助作用。未來研究中，如何綜合利用多種數據類型蘊含的地質信息、如何有效整合現有方法、以及如何將機器學習等思想引入米氏旋回研究等問題，值得科技工作者繼續開展研究。

致謝 論文評審過程中，收獲了匿名評審專家詳盡、寶貴的修改建議；修改過程中，得到了編輯部工作人員耐心、細致的指導，在此一并致以真誠的謝意！