宋翠玉，呂大煒

山東科技大學(xué)地球科學(xué)與工程學(xué)院，山東青島 266590

0 引言

周期性現(xiàn)象普遍存在于地球的各個角落。在地質(zhì)記錄中，周期性表現(xiàn)為韻律，是地層學(xué)研究最早發(fā)現(xiàn)的現(xiàn)象之一。地層學(xué)里把長而復(fù)雜的周期性節(jié)律稱為旋回[1]。1988 年在意大利Perugia 召開的一次沉積學(xué)會議上，旋回地層學(xué)（Cyclostratigraphy）這一術(shù)語被首次公開使用[2]。對地層記錄的（準）周期性旋回變化進行識別、描述、對比和成因解釋，并應(yīng)用于地質(zhì)年代學(xué)以提高地層年代框架的精度和分辨率，實現(xiàn)高精度地層劃分與對比，是旋回地層學(xué)的基本任務(wù)[3]。沉積地層中的韻律或旋回有外周期型和內(nèi)周期型兩類，前者由周期變化的外力（如氣候、構(gòu)造、海平面變化等）驅(qū)動，后者則產(chǎn)生于沉積系統(tǒng)內(nèi)部，如一次沉積事件造成的沉積搬運等[4]。軌道周期是地球系統(tǒng)中十分常見的外周期。20世紀40年代，塞爾維亞機械工程師米蘭科維奇提出了第四紀冰期成因的天文假說[5]。他認為地球自轉(zhuǎn)和公轉(zhuǎn)的三個重要的軌道參數(shù)：偏心率、黃赤交角或地軸傾斜度和歲差受其它天文因素的影響發(fā)生周期性變化，而這些變化引起了太陽輻射量的變化，地表北緯65 度附近接收到夏季太陽輻射量的變化是驅(qū)動第四紀冰期旋回的主要原因[5]。米蘭科維奇理論的冰期旋回，至今是地球科學(xué)中唯一公認有定量基礎(chǔ)的周期變化[4]。為了突出地球軌道參數(shù)變化對地層沉積作用的影響，徐道一等[6]在旋回地層學(xué)基礎(chǔ)上，基于天文學(xué)的地球軌道三要素的周期性和準周期性變化提出了“天文地層學(xué)”這一概念，指出可通過功率譜分析和數(shù)字濾波等信號處理與運算進行地層劃分與對比[7]。地層中由地球軌道驅(qū)動造成的旋回性記錄稱之為米蘭科維奇旋回[8]，簡稱米氏旋回，是旋回地層學(xué)研究的重點。

目前，米氏旋回很大程度上解釋了許多地質(zhì)學(xué)領(lǐng)域的問題。越來越多研究成果表明，氣候變化只是軌道力驅(qū)動的一個方面，從海水混合到巖漿作用，也都對軌道周期有所響應(yīng)[4,9]。在不同地質(zhì)歷史時期的海相、陸相地層中，米氏旋回正不斷地被揭示出來[8,10]。識別并研究地質(zhì)記錄中的米氏旋回信號，對于提高定年精度[11-15]、構(gòu)建高分辨率年代地層格架[16-18]、估算地層剝蝕量[19-21]、研究古氣候—古環(huán)境演化[22-26]等均具有重要意義。

米氏旋回的研究方法總體上可歸為兩類[8]：一類是巖性直觀識別方法，通過直接觀察巖性和巖相的變化，根據(jù)地層的組合特征、方式及級序結(jié)構(gòu)來判別是否存在旋回信號[12,27]。該類方法要求研究者具有豐富的野外工作經(jīng)驗，能夠準確地把握巖性的識別和讀取，需要有雄厚的地質(zhì)功底。但是隨著研究深入，人們逐漸發(fā)現(xiàn)：在深湖、遠洋、半遠洋等沉積環(huán)境中，有些軌道參數(shù)的變化未必能引起明顯的巖性或巖相變化，故直接識別法不僅難度大，還容易于造成旋回信息的遺漏[8]。第二類是時間序列分析法，利用巖性巖相特征的數(shù)字化值（如特殊巖石的顏色變化等）或者地層的古氣候替代指標等（如地層的測井曲線、地層中某種化學(xué)元素的含量變化等）構(gòu)建包含地層環(huán)境信息的時間（深度）序列，進而利用數(shù)學(xué)變換對這些數(shù)據(jù)序列進行定量分析[2,28]。這種方法需要借助大量的先進技術(shù)手段獲取數(shù)據(jù)序列，比如自然伽馬（Gamma Ray，GR）測 井 設(shè) 備[29-30]、XRF（X Ray Fluorescence）元素掃描儀[31-34]、顏色識別儀[33,35]、磁化率測試設(shè)備[32-33]等，并且也要求相關(guān)人員具有較強的地質(zhì)學(xué)基礎(chǔ)。應(yīng)該說，時間序列分析法是目前最常用的旋回地層分析法，它極大地推進了旋回地層學(xué)的發(fā)展。

綜上所述，米蘭科維奇理論已經(jīng)大大促進了地層學(xué)研究的發(fā)展，但是隨著數(shù)學(xué)方法的不斷革新，米氏旋回研究的方法學(xué)也相應(yīng)隨之改進和提高，為此，本文在近幾年的米氏旋回方法學(xué)研究基礎(chǔ)上，針對時間序列分析法的主要步驟，將分別從地層序列的數(shù)據(jù)類型、天文檢驗及天文調(diào)諧等三方面，簡述現(xiàn)有方法的基本原理，總結(jié)其優(yōu)勢和局限性，討論存在的問題，并展望時間序列分析法的發(fā)展方向，為米氏旋回方法學(xué)研究提供一定的基礎(chǔ)。

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

利用時間序列分析法進行米氏旋回研究的工作方法較為成熟。主要涉及地層序列的數(shù)據(jù)類型選取、數(shù)據(jù)預(yù)處理、頻譜分析、天文驅(qū)動檢驗和天文調(diào)諧等環(huán)節(jié)和內(nèi)容，其一般流程如圖1所示。

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

不同的地層數(shù)據(jù)序列蘊含著不同的地質(zhì)信息，對米氏旋回信號的記錄亦有差異，因此，選擇合適的數(shù)據(jù)序列是米氏旋回分析的首要且關(guān)鍵步驟。數(shù)據(jù)預(yù)處理是指頻譜分析之前對原始數(shù)據(jù)進行的一系列處理。主要包括重采樣、剔除異常值、去趨勢化和預(yù)白化等，其目的是盡量消除地層數(shù)據(jù)中的各種環(huán)境“噪聲”，使結(jié)果更容易解釋[8,32]。頻譜分析將地層序列的信號強度按頻率順序展開，使其成為頻率的函數(shù)，進而識別出信號中（準）周期性的成分[8]。它是檢驗地層中是否包含米氏旋回記錄的前序處理。通過頻譜分析，可以獲得數(shù)據(jù)中蘊含的一系列主周期信號，天文驅(qū)動檢驗就是檢驗這些周期信號中是否包含地球天文軌道驅(qū)動而成的米蘭科維奇周期信號。如果地層中存在米氏旋回信號，則可以通過天文調(diào)諧將從數(shù)據(jù)序列提取出的旋回記錄對比到歲差、斜率和（或）偏心率的天文目標曲線上，使深度域的地層沉積序列變換到時間域，以建立天文年代標尺或浮動天文年代標尺，用于定年或確定地層∕地質(zhì)事件的持續(xù)時間。

目前，已有多種米氏旋回時間序列分析方法，如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 地層序列的數(shù)據(jù)類型

目前，應(yīng)用于米氏旋回研究的數(shù)據(jù)序列可大致歸為四類（表1）。第一類為直觀反映巖石物理性質(zhì)的數(shù)據(jù)，例如顏色[33,35]、粒度、相變與層理韻律[33,41]等；第二類為地球化學(xué)數(shù)據(jù)，如Fe 含量[38]、CaCO3質(zhì)量分數(shù)[36]、碳同位素[42-44]、主量元素含量[24]、微量元素含量[24]、元素比值[31,45]等；第三類為地球物理數(shù)據(jù)，如磁化率[22,24,32-33,45-46]、非磁滯剩磁（Anhysteretic Remanent Magnetization，ARM）[32-33]、GR（Gamma Ray）測 井 數(shù)據(jù)[23,29,33,35,46-49]、電阻率測井數(shù)據(jù)[46]等。最后一類為古生物相關(guān)的數(shù)據(jù)，例如生物豐度、生物滅絕速率[32,50]等。

表1 用于米氏旋回研究的數(shù)據(jù)序列類型列表Table 1 Data series types used in Milankovitch cycle study

由于對古環(huán)境、古氣候變化反映敏感且測試較方便，巖石磁性特征和自然伽馬測井數(shù)據(jù)等地球物理指標被廣泛應(yīng)用于旋回地層學(xué)分析中。磁性礦物是氣候和環(huán)境變化信息的重要載體。表征物質(zhì)磁學(xué)特征的磁化率則反映了巖石中所含磁性礦物的種類、含量和粒徑等信息，是良好的古氣候及古環(huán)境替代指標[22]。GR 測井測量沉積物中伽馬射線的強度，可以敏感反映沉積物中泥質(zhì)含量，進而反映古氣候和古環(huán)境的微小變化，是米氏旋回分析的理想數(shù)據(jù)[29]。一般而言，GR 高值對應(yīng)富含黏土沉積物，GR低值與砂巖、富含碳酸鹽的沉積物有關(guān)[52]。自然伽馬能譜測井除了獲取得出的GR值外，還記錄鉀（K）、鈾（U）和釷（Th）等元素的含量[30,33,52]。

不同數(shù)據(jù)指標對天文軌道驅(qū)動和非天文噪聲的響應(yīng)是有差異的[33,35,53]。例如，在海相沉積地層中，ARM 和Th∕U 可反映早三疊世的內(nèi)陸風(fēng)化作用，GR、K、U、Th、磁化率和非碳酸鹽組分等指示陸源補給，而巖石顏色數(shù)字化后生成的CIE 分量L*、a*以及巖性序列（lithologic rank）則分別指示生產(chǎn)力、氧化還原狀態(tài)及相對海平面[33]。此外，有些數(shù)據(jù)指標受軌道力驅(qū)動，而有些指標可能主要受非軌道噪聲的影響[35]。因此，基于單一沉積序列構(gòu)建的天文年代標尺具有不確定性。

為了減弱單一數(shù)據(jù)序列中環(huán)境“噪聲”造成的不確定性，同時綜合利用不同數(shù)據(jù)序列中蘊含的地質(zhì)信息，已有大量研究利用多種數(shù)據(jù)指標進行旋回地層學(xué)研究。這些研究中，往往首先將多個數(shù)據(jù)指標逐一進行時間序列分析，然后聯(lián)合解讀它們對米氏旋回的指示作用，進而進行地層劃分或者古環(huán)境研究[30-32,54]。從松遼盆地晚白堊世嫩江組的地球化學(xué)元素比值Rb∕Sr、K∕Ti、Ti∕Al 和Zr∕Rb 中識別的米氏旋回信號，揭示了嫩江組沉積時期氣候由濕潤向半干旱、半濕潤過渡的古氣候響應(yīng)機制[31]。GR 和SP（Spontaneous Potential）測井曲線的頻譜分析結(jié)果表明，北黃海東部坳陷始新統(tǒng)地層很好地保存了米氏旋回，據(jù)此可進一步分析其沉積環(huán)境并進行地層精細劃分[54]。天津薊縣剖面前寒武系洪水莊組—鐵嶺組地層中，磁化率和伽馬能譜數(shù)據(jù)指標中均記錄了完整的米氏旋回[30]。房強[32]利用磁化率、ARM 和地球化學(xué)等多個指標對上寺和渡口剖面茅口組進行旋回地層學(xué)分析，提取了米氏旋回在晚古生代冰期末期的記錄，并討論了其古環(huán)境響應(yīng)。上述研究盡管采用了多種數(shù)據(jù)序列進行旋回分析，但它們均作為單一輸入獨立進行了時間序列分析，并未真正實現(xiàn)多種地層信息的有效整合。近年來，融合多種數(shù)據(jù)類型的米氏旋回分析方法正引起關(guān)注。根據(jù)多種數(shù)據(jù)指標的獨立模型來構(gòu)建綜合的年代模型，被證實是估計年齡模型不確定性的有效方法[35]。此外，對多種數(shù)據(jù)指標先做主成分變換等處理，繼而對某一個或者某幾個分量進行進行米氏旋回分析，也有助于識別完整的米蘭科維奇周期信號[46]。

1.2 天文驅(qū)動檢驗方法進展

米氏旋回研究解決的核心問題之一，就是檢驗地層記錄中是否存在米氏周期信號，即檢驗沉積記錄中的周期信號是否由地球天文軌道驅(qū)動[40,55]。周期的倒數(shù)是頻率，因此，通過頻譜分析把深度域的沉積序列觀測值變換到頻率域，是檢驗天文驅(qū)動信號的關(guān)鍵步驟。

1.2.1 頻譜分析方法研究現(xiàn)狀

旋回地層分析中，常用的頻譜分析方法有快速傅里葉變換（Fast Fourier Transform，F(xiàn)FT）、周期圖法（Periodogram）、Welch 法、多窗譜分析法（Multi-taper Method，MTM）和 最 大 熵 譜 法（Maximum Entropy Method，MEM）等。不同的頻譜分析方法所基于的算法不同，優(yōu)缺點也各不相同，但計算結(jié)果基本一致[8]。

傅里葉變換是最經(jīng)典的頻譜分析方法之一。它把在時間域里看起來很復(fù)雜的多頻率疊加信號分解成不同頻率的正弦分量，將空間域變換到頻率域[56]。由于實際數(shù)據(jù)的長度有限，通常采用離散傅里葉變換（Discrete Fourier Transform，DFT）對數(shù)據(jù)進行譜分析，即對有限長度的離散時間序列進行傅里葉變換?？焖俑道锶~變換（FFT）是DFT的快速算法，將其運算量減少了幾個數(shù)量級。

周期圖法是一種算法簡便且具代表性的功率譜估計方法，其原理是對觀測到的數(shù)據(jù)直接進行傅里葉變換，然后取模的平方獲得功率譜。周期圖是對功率譜密度（Power Spectral Density，PSD）的有偏估計。在信噪比大、數(shù)據(jù)夠長的時候，周期圖是有用的譜估計器。但在數(shù)據(jù)序列不夠長時，周期圖法的分辨率往往受限。Liet al.[38]提出的eCOCO 和Meyers[39]提出的TimeOpt 方法中，頻譜分析法均采用周期圖法。

Welch法是對周期圖法的一種改進，它先將數(shù)據(jù)序列劃分為不同的段（可以有重疊），然后對每段進行加窗處理，再分別計算周期圖后取平均。該方法能改善方差性能和分辨率，是較常用的譜估計方法之一[2]。

MTM 是由Thompson 提出的一種多窗譜分析和信號重建方法，能提供更優(yōu)的PSD估計[57]。它使用一組正交、離散的扁平類球體序列（Discrete Prolate Spheroidal Sequences，DPSS，也叫做Slepian 序列）獲得最優(yōu)濾波器，進而計算功率譜估計值。通過調(diào)整MTM 中的時間—帶寬參數(shù)，可以實現(xiàn)估計方差和分辨率之間的平衡。在非常短且有噪聲的時間序列中，MTM具有較高的頻率分辨率，并且為每個頻率提供了統(tǒng)計顯著性檢驗（F-方差比檢驗）。此外，這種顯著性檢驗與周期的振幅無關(guān)，因此，該方法不僅能夠識別顯著性水平高的低振幅峰值，還能指示哪些高振幅峰值可能不具有統(tǒng)計顯著性。旋回地層學(xué)研究中，MTM是應(yīng)用較廣的頻譜分析方法之一[30,36,58-59]。

MEM 的原理是根據(jù)已知數(shù)據(jù)信息，在不進行任何新的假設(shè)的情況下，按信息熵最大準則預(yù)測未知延遲離散時間上的相關(guān)函數(shù)[60]。MEM法的優(yōu)勢在于解決含有低噪聲的間隔緊密的正弦短數(shù)據(jù)信號。與其他譜分析方法（如周期圖法）相比較，MEM 具有不受取樣長度限制等優(yōu)點[41,61]。但受噪聲、正弦信號初始相位等因素的影響，MEM 譜估計會產(chǎn)生頻率偏移。

圖2列出了某剖面3 130 m至3 441 m深度GR序列的5種頻譜估計結(jié)果，采樣間隔為5 cm，共有6 219個數(shù)據(jù)點?？梢钥闯觯現(xiàn)FT 與Welch、MTM 與周期圖法獲得的功率譜形態(tài)分別有一定相似性，而MEM分辨率低，譜非常光滑，難以篩選米氏周期信號。與FFT 相比，Welch 獲得的譜更光滑，分辨率降低了。MTM 和周期圖法獲得的譜的分辨率均比其它方法高，但二者功率譜極值對應(yīng)的頻率位置具有差異性。MTM的計算復(fù)雜度較高，是計算DPSS的代價。

圖2 基于不同方法的某剖面GR 數(shù)據(jù)頻譜分析結(jié)果Fig.2 Spectral analysis results of GR series in an unknown section using five methods

頻譜分析的作用就是評估能量較高的非隨機頻段是否與米氏旋回的頻率吻合。對沉積序列數(shù)據(jù)做頻譜分析后，功率譜峰值對應(yīng)的頻率即為米氏旋回信號的候選值。然而，這些周期（頻率的倒數(shù)）信號既包含真實的米氏旋回信號，也包括一些無關(guān)的、隨機出現(xiàn)的“噪聲”。這些噪聲的能量可能較低，但或許會出現(xiàn)在所有頻率中。研究表明，沉積記錄中的噪聲在低頻譜段能量較高，隨頻率增大能量逐漸降低，與光類似，故稱為“紅噪”（Red Noise）[62-63]。紅噪的出現(xiàn)表明隨機振蕩變化是氣候的一個普遍特征。即使沒有氣候驅(qū)動參數(shù)（如日照量）的變化，氣候系統(tǒng)中的各組成部分及其相互作用也可以產(chǎn)生強烈的低頻振蕩[40]。因此，需要對噪聲進行估計，并通過顯著性檢驗將這種隨機振蕩與天文驅(qū)動信號區(qū)分出來。

常采用一階自回歸模型（Auto Regressive Lag-1，AR1）估計沉積序列中紅噪聲，生成紅噪聲理論曲線并進行顯著性檢驗，將設(shè)定置信水平（90%、95%或99%）之下的峰值對應(yīng)的頻率剔除。自回歸Lag-1（AR1）隨機噪聲模型是天文年代學(xué)研究中評價旋回顯著性的最常用方法之一，它是一種直接基于自身數(shù)據(jù)的估計模型[2]。圖3為Chaohu剖面不同深度GR數(shù)據(jù)的MTM 頻譜分析結(jié)果，噪聲估計均基于魯棒的紅噪模型[51]。對于10~74 m 的地層而言，置信水平95%以上的峰值對應(yīng)波長有10 m、2.3 m、1 m、0.8 m、0.7 m和0.5 m。

1.2.2 天文檢驗方法研究現(xiàn)狀

在噪聲壓制后的功率譜中，判別峰值對應(yīng)的頻段組合是否由天文軌道驅(qū)動的方法主要有三種，即比值法、調(diào)幅分析法和假設(shè)檢驗法。

一般而言，通過頻譜分析獲取的旋回周期之比與天文軌道參數(shù)的周期之比一致，即可認為沉積記錄中包含米氏旋回。例如，地球天文軌道參數(shù)長偏心率、短偏心率、斜率和歲差的周期之比約為20∶5∶2∶1，如果頻譜分析中提取的四個周期（頻率取倒數(shù)）之比與此接近，則可認為存在米氏旋回。該方法稱為比值法，是目前絕大多數(shù)旋回地層研究采用的天文檢驗方法[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 數(shù)據(jù)的2π MTM 的頻譜分析結(jié)果[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

長周期旋回對短周期旋回會有明顯的天文調(diào)制，例如短偏心率對歲差進行天文調(diào)制，長偏心率對短偏心率天文調(diào)制[66]。所以，調(diào)幅分析也是一種判斷觀察到的旋回是否受天文軌道控制的有效方法[67-68]。然而，沉積序列數(shù)據(jù)中類似偏心率的振幅變化可能通過調(diào)諧和數(shù)據(jù)處理等過程人為地引入[67]。因此，振幅調(diào)制方法不能作為一種獨立的方法來測試天文調(diào)諧時間尺度的準確性。

檢驗天文驅(qū)動的另一類常用方法為統(tǒng)計學(xué)中的零假設(shè)檢驗法。該類方法先提出零假設(shè)，即假設(shè)篩選出來的頻率都是與天文驅(qū)動無關(guān)的隨機噪聲，然后借助統(tǒng)計檢驗拒絕該假設(shè)，即可說明序列中含有米氏旋回信號。Meyerset al.[36]提出的ASM 方法，Liet al.[38]提出的eCOCO 方法中，均采用零假設(shè)檢驗天文驅(qū)動信號的存在。但當軌道信號因沉積速率的變化而嚴重失真時，ASM 將無法以高置信水平拒絕零假設(shè)[37]。

1.3 天文調(diào)諧方法研究進展

在地層中識別出完整的米氏旋回信號后，就可以將信號調(diào)諧到理論的天文曲線上[69]。通過年代校準即可建立高分辨率地層層序格架。在年齡控制較差和缺少天文目標曲線的古生代和中生代，調(diào)諧只能建立浮動天文年代標尺，來精確地計算地層或地質(zhì)事件的持續(xù)時間。天文調(diào)諧方法可歸納為三類（表2），即最小調(diào)諧法、統(tǒng)計法和反演法等[38]。

表2 主要天文調(diào)諧方法一覽表Table 2 Main astronomical tuning methods

1.3.1 最小調(diào)諧法

最小調(diào)諧法（Minimal tuning）[70]使用最小數(shù)量的校準點將深度域序列曲線與時間域目標天文曲線聯(lián)系起來完成校正，是最經(jīng)典、也是應(yīng)用最廣泛的調(diào)諧方法[51,55]。該方法先將沉積序列數(shù)據(jù)調(diào)諧到單個天文頻率（如偏心率）上，然后對調(diào)諧后的結(jié)果再進行頻譜分析，如果也能檢驗出其它的天文軌道信號（如斜率、歲差）則調(diào)諧成功[32]。實際操作中，為了更準確地定位校準點，常先將沉積序列曲線做帶通濾波（Band Pass Filter）處理。這種調(diào)諧方法直觀簡單，但存在局限性：首先，該方法假定目標曲線正確，沒有利用統(tǒng)計法進行檢驗[55]；其次，選取的目標曲線年代與真實數(shù)據(jù)年代難以完全對應(yīng)，從而產(chǎn)生相位差。

在缺少天文目標曲線的古生代，頻率域的最小調(diào)諧往往獲得理想效果。首先，利用EHA（Evolutive Harmonic Analysis）[71]或改進的FFT[24,48]在滑動窗口上對地層數(shù)據(jù)序列做頻譜分析，獲得頻譜圖，該圖表征頻率信號在深度上的偏移規(guī)律。然后，在頻譜圖中追蹤單個較連續(xù)的旋回信號（如偏心率周期），獲得沉積速率曲線。最后，根據(jù)沉積速率曲線將沉積記錄調(diào)諧到時間域。圖4a為某剖面GR數(shù)據(jù)經(jīng)EHA處理后的頻譜圖，滑動窗口長度為60 m，圖中白點表示追蹤到的長偏心率周期（405 kyr），據(jù)此計算出沉積速率曲線（圖4b）。由于缺少年齡控制，經(jīng)調(diào)諧建立了該剖面的浮動年代標尺（圖4c）。

圖4 頻率域最小調(diào)諧法示意圖（a）EHA頻譜；（b）沉積速率曲線；（c）天文調(diào)諧結(jié)果Fig.4 Sketch map of Frequency Domain Minimal tuning(a)EHA spectrum;(b)sedimentation rate curve;and(c)tuning result

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

圖5 為樊頁1井磁化率指標的ASM分析結(jié)果[17]，沙三下亞段最優(yōu)沉積速率為10.801 cm∕kyr，零假設(shè)檢驗的顯著性水平低于0.1%，表明沉積速率為10.801 cm∕

1.3.2 統(tǒng)計法

統(tǒng)計法基于一系列可能的沉積速率，利用假設(shè)檢驗和蒙特卡洛模擬來估算最優(yōu)沉積速率或追蹤沉積速率的變化。這些方法不需要精確的年代限制，也能更好地保持數(shù)據(jù)的原始相位。ASM[36]及其改進方法eASM[37]、TimeOpt[39]及eCOCO[38]等方法都是統(tǒng)計法。

給定一系列可能的沉積速率，ASM 對天文軌道驅(qū)動的信號進行零假設(shè)檢驗，同時通過蒙特卡洛模擬來評價識別出來的優(yōu)勢頻率與目標軌道頻率的擬合程度，見式（1），根據(jù)ASM 值極小及低零假設(shè)水平確定最優(yōu)沉積速率的客觀估計值。

式中：f為功率譜峰值對應(yīng)的頻率（單位：cycles∕m）；s為待測試的沉積速率（單位：m∕kyr）；fpred為頻率的理論值（周期取倒數(shù)，單位：cycles∕kyr）；ΔfR是最小分辨率帶寬。kyr 的置信度超過99.9%。然而需要注意的是，ASM方法僅是量化匹配的最優(yōu)沉積速率，并不是平均沉積速率，也不是計算某一段的精確沉積速率[17]。

eASM（evolutive ASM）[37]對ASM做了四點改進：①對式（1）中的αk做了調(diào)整（見式（2）），以更好地考慮譜峰位置的不確定性；②對測試的沉積速率做對數(shù)運算，以更好地評估低沉積速率時經(jīng)常出現(xiàn)的ASM快速變化；③引入了滑動窗口，可以適應(yīng)地層區(qū)間沉積速率的變化；④考慮理論軌道目標頻率（fpred）及其不確定性。

Meyers[39]提出的TimeOpt也是一種統(tǒng)計法。該方法基于一系列可能的沉積速率：①利用概率線性回歸法，根據(jù)公式（3）重建歲差的偏心率調(diào)幅，并計算重建結(jié)果與實際調(diào)幅結(jié)果的相關(guān)系數(shù)r2envelope；

yenvelope=Xeβe+ε（3）

式中：yenvelope是時間校準數(shù)據(jù)序列的歲差振幅包絡(luò)線，Xe是代表偏心率周期的正弦和余弦預(yù)測項矩陣，βe是每個預(yù)測值的回歸系數(shù)向量，ε是誤差項向量；

②再次利用概率線性回歸法，根據(jù)式（4）重建包含歲差和偏心率理論值的功率譜，并計算重建結(jié)果與年代校準數(shù)據(jù)之間的相關(guān)系數(shù)r2spectral；

ydata=Xepβep+ε（4）

式中：ydata是時間校準數(shù)據(jù)序列，Xep是代表偏心率和歲差周期的正弦和余弦預(yù)測項矩陣，βep是每個預(yù)測值的回歸系數(shù)向量，ε是誤差項向量；

③乘積r2opt=r2evenlope r2spectral極大值對應(yīng)的沉積速率為最優(yōu)。

圖6 為ODP Site 1 262 顏色值a*（綠∕紅）的TimeOpt分析結(jié)果[39]。由圖6e可以看出，包絡(luò)回歸模型（公式（3））確定的最優(yōu)沉積速率為1.28 cm∕kyr，譜功率回歸模型（公式（4））和聯(lián)合模型（乘積）確定的最優(yōu)沉積速率均為1.33 cm∕kyr，此時，最大r2opt為0.212（圖6f）。經(jīng)過2 000 次AR1 蒙特卡洛模擬，得出在1.33 cm∕kyr 時觀察到的最大r2opt的p 值為0.005，表明我們可以在99.5%的置信水平下拒絕零假設(shè)（圖6g）。在該地層中，TimeOpt得出的沉積速率為1.33 cm∕kyr，與先前基于天文年代學(xué)、磁性地層學(xué)和生物地層學(xué)的估計值（如1.3 cm∕kyr，1.2 cm∕kyr，1.9 cm∕kyr等）相似，證實了該方法在沉積速率估算中的有效性。

圖6 ODP Site 1 262 顏色值a*（綠∕紅）的TimeOpt 結(jié)果（據(jù)Meyers[39]修改）（a）1 262顏色數(shù)據(jù)；（b）顏色數(shù)據(jù)的周期圖，TimeOpt確定的沉積速率為1.33 cm∕kyr（黑線、深灰線分別為線性、對數(shù)光譜），灰色陰影區(qū)為用于評估歲差振幅包絡(luò)的部分，灰色虛線表示偏心率和歲差的目標周期；（c）帶通歲差信號（黑色）與Hilbert變換獲得的數(shù)據(jù)振幅包絡(luò)（灰色）對比；（d）數(shù)據(jù)振幅包絡(luò)（灰色）與TimeOpt重建的偏心率模型（黑色；據(jù)公式（3））對比；（e）每個待測沉積速率下振幅包絡(luò)擬合的平方皮爾遜相關(guān)系數(shù)（r2envelope；灰色線）與譜功率擬合的相關(guān)系數(shù)（r2spectral；黑色線）；（f）每個待測沉積速率下包絡(luò)線和譜功率的聯(lián)合擬合（r2opt）；（g）基于AR1的2 000次蒙特卡羅模擬結(jié)果（ρAR1=0.907），用于評估最大r2op（t0.212）的顯著性；（h）數(shù)據(jù)振幅包絡(luò)線與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]提出的一種基于演化的相關(guān)系數(shù)的米氏旋回分析方法。對于顯生宙古氣候序列而言，eCOCO 能夠同時評價沉積速率和天文驅(qū)動。它基于一系列待“測試”的沉積速率，1）將古氣候替代成分轉(zhuǎn)換為時間序列，然后計算該時間序列的功率譜與天文解決方案功率譜之間的相關(guān)系數(shù)；2）記錄對相關(guān)系數(shù)有貢獻的天文參數(shù)的數(shù)量；3）使用蒙特卡洛模擬法對天文驅(qū)動信號的零假設(shè)進行檢驗。最合適的沉積速率由高相關(guān)系數(shù)、多旋回信號及低零假設(shè)水平共同決定。由于eCOCO 采用了滑動窗口技術(shù)，可以跟蹤替代成分沿地層序列在沉積速率上的變化。圖7為晚三疊世Newark深度序列（0~2 000 m）的eCOCO分析結(jié)果[38]，其中暗紅色（極值）交集決定了特定深度的最優(yōu)沉積速率（圖7b~d）。

圖7 晚三疊世Newark 深度序列（0~2 000 m）的eCOCO 分析結(jié)果[38]（a）深度序列的高斯低通濾波結(jié)果（截止頻率為0.01 m-1），金色水平條為～270 m旋回；（b）eCOCO相關(guān)系數(shù)圖，黑線和綠線分別為Kentet al.[72]和Liet al.[73]的研究結(jié)果；（c）零假設(shè)檢驗H0顯著性水平圖；（d）有貢獻的天文參數(shù)數(shù)量變化圖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 反演法

反演法廣泛應(yīng)用于地球物理的各個領(lǐng)域，其本質(zhì)為基于一定的數(shù)學(xué)方法將觀測數(shù)據(jù)映射到相應(yīng)的地球物理模型[74]。例如，通過地表觀測的地震波數(shù)據(jù)，得到地球的區(qū)域速度結(jié)構(gòu)，或通過觀測的重力數(shù)據(jù)，獲取地殼的密度分布等信息。

Malinvernoet al.[44]將軌道調(diào)諧視作反演問題，提出了基于貝葉斯反演的天文調(diào)諧法。該方法將變化的沉積速率及相應(yīng)的地層深度作為參數(shù)建立年代模型（Age Model），首先應(yīng)用貝葉斯公式定義軌道周期信號驅(qū)動的沉積速率變化的概率分布函數(shù)（而非單一沉積速率）。然后，通過馬爾科夫鏈蒙特卡洛法（Markov Chain Monte Carlo）對概率分布進行采樣，量化沉積速率的不確定性，這些不確定性來源于調(diào)諧周期的不確定性，以及與軌道周期無關(guān)的沉積信號成分。最優(yōu)沉積速率選取似然度（likelihood）較大的沉積速率。

eASM、eCOCO 等統(tǒng)計法在地層記錄中按順序移動窗口，有助于追蹤地層中變化的沉積速率信息，但無法獲取窗口中發(fā)生的沉積速率的顯著變化。因此，Meyers[40]在TimeOpt 基礎(chǔ)上，引入沉積模板（sedimentation template）（圖8）構(gòu)建復(fù)雜的沉積模型，提出擴展的TimeOpt 反演法（對應(yīng)Astrochron 工具包中的timeOptTemplate 函數(shù)）。在該優(yōu)化反演過程中，初始沉積模板與天文目標和數(shù)據(jù)動態(tài)交互調(diào)整，最終重建沉積速率曲線。

圖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）法。該方法定量地將天文學(xué)理論與地質(zhì)數(shù)據(jù)聯(lián)系起來，旨在克服二者的局限性。TimeOptMCMC 的核心有三部分：1）TimeOpt 法，考慮時間尺度的不確定性，并利用天文信號的多個屬性來提高統(tǒng)計可靠性；2）天文學(xué)理論，將歲差、軌道偏心率周期與太陽系、地—月演化的基本頻率聯(lián)系起來；3）貝葉斯馬爾可夫鏈蒙特卡羅方法，探討數(shù)據(jù)、模型空間及不確定性。

2 存在問題

經(jīng)過眾多研究者幾十年的努力，米氏旋回的研究方法不斷發(fā)展與優(yōu)化，基于時間序列分析的新方法也不斷涌現(xiàn)。天文解決方案精度的不斷提高，也為米氏旋回研究提供了重要基礎(chǔ)。然而，縱觀現(xiàn)有研究，米氏旋回時間序列分析法尚存在如下問題：

（1）數(shù)據(jù)序列中存在的噪聲影響天文檢驗和調(diào)諧的精度。地層數(shù)據(jù)序列既記錄了軌道驅(qū)動的環(huán)境變化信息，也包含了一些“噪聲”。這些噪聲可能來源于地層的古沉積環(huán)境，如沉積速率變化、生物擾動、滑塌作用和后期成巖作用等，也可能是在數(shù)據(jù)準備過程中引入的，如采樣處理、儀器噪聲等。在數(shù)據(jù)預(yù)處理環(huán)節(jié)進行去趨勢化及預(yù)白化等處理、使用對噪聲魯棒的頻譜分析方法、對頻譜分析結(jié)果進行紅噪聲抑制等做法均會在一定程度上去除數(shù)據(jù)中的噪聲，但頻譜分析結(jié)果中仍可能產(chǎn)生大量的譜峰[68]，既包含軌道驅(qū)動的周期信號，也包含非軌道驅(qū)動的周期信號，給后續(xù)處理和研究都帶來很大干擾。

（2）新型的定量化的天文調(diào)諧方法仍需要進一步驗證及改進。米氏旋回研究中，絕大多數(shù)天文調(diào)諧方法采用傳統(tǒng)的最小調(diào)諧法[70,76-77]，其次是借助EHA[31]或改進的FFT 等頻率域最小調(diào)諧法[24,78]，而統(tǒng)計法和反演法應(yīng)用相對較少[42,78]。傳統(tǒng)調(diào)諧法直觀性好，但通常需要交互操作，調(diào)諧過程中濾波主頻段的選取、極小值∕極大值位置、目標天文曲線時間段選擇等難免會受處理者主觀因素的影響，不同處理者或同一處理者在不同時間的調(diào)諧結(jié)果都可能有所不同，導(dǎo)致研究的可重復(fù)性較差。統(tǒng)計法和反演法基于假設(shè)檢驗和概率等數(shù)學(xué)原理，能評價各種噪聲帶來的天文信號的不確定性，為定量、客觀的沉積速率估算提供了重要支持。現(xiàn)有方法不僅能夠獲取地層序列的最優(yōu)沉積速率（如ASM，TimeOpt 等），還能夠追蹤沉積速率的變化信息（如eASM，eCOCO，eTimeOpt）。然而，由于應(yīng)用相對較少，這些新型方法在不同數(shù)據(jù)、不同條件下的天文調(diào)諧精度還有待于進一步驗證。此外，對于統(tǒng)計法和反演法而言，參數(shù)調(diào)整中存在不確定性、缺少中生代和古生代精確的天文目標曲線等因素，也約束著這些方法的適用性及準確性。

3 研究展望

基于時間序列分析的米氏旋回研究方法仍處在發(fā)展階段，目前已經(jīng)出現(xiàn)了一些較為成熟的軟件或工 具 包，如AnalySeries、Astrochron (R package)[40]、Acycle[79]等。地球科學(xué)信息化進程的不斷推進，對研究方法智能化水平和可靠性等方面的要求更高。未來挑戰(zhàn)與機遇并存，作者認為以下兩方面值得開展進一步研究工作。

（1）沉積序列數(shù)據(jù)的選擇與信息優(yōu)化。數(shù)據(jù)是時間序列分析的基礎(chǔ)，其信息有效性對分析結(jié)果起著關(guān)鍵作用。不同的數(shù)據(jù)序列對地層特征及環(huán)境變化的記錄各有側(cè)重，因此，使用單一數(shù)據(jù)類型進行旋回地層分析可能存在局限性[33]。對同一地質(zhì)剖面而言，利用不同數(shù)據(jù)序列提取的米氏旋回頻率組合也可能具有多解性[46]。為了減小數(shù)據(jù)類型差異性帶來的不確定性，需要在以下兩方面繼續(xù)開展工作。

第一，需要建立定量的評價標準，在數(shù)據(jù)可獲取的前提下，選擇更適用于特定環(huán)境下旋回地層分析的數(shù)據(jù)序列。信噪比、層次聚類、能量譜分解等已被分別用于評價單序列的信息質(zhì)量、多沉積序列間的關(guān)系及不同數(shù)據(jù)序列對外來氣候驅(qū)動的敏感程度等[33]。形成穩(wěn)定、可靠的評價體系仍需深入工作。

第二，當同一地質(zhì)剖面有多個數(shù)據(jù)序列可用時，如何綜合使用多個數(shù)據(jù)序列的優(yōu)化信息，亦值得研究。在數(shù)據(jù)處理領(lǐng)域，比值運算、主成分變換、獨立成分變換等是常用的特征提取方法。不同沉積數(shù)據(jù)序列可以視作描述地層環(huán)境的多個特征，因此，在沉積理論的指導(dǎo)下，可借助上述特征提取方法獲取地層環(huán)境的優(yōu)化信息，作為米氏旋回研究的數(shù)據(jù)基礎(chǔ)。

（2）基于新方法的沉積速率客觀、量化估算。準確獲取沉積速率曲線是米氏旋回定量化研究的關(guān)鍵環(huán)節(jié)。天文調(diào)諧的統(tǒng)計法和反演法都是基于沉積速率估算的時間序列分析方法。統(tǒng)計檢驗和蒙特卡洛模擬等數(shù)學(xué)手段已經(jīng)成功應(yīng)用于天文檢驗和調(diào)諧中。不同估算方法的原理存在較大差別，但也不乏關(guān)聯(lián)性。面對逐漸多樣化的時間序列分析法，既需要有效整合現(xiàn)有方法，形成理論更嚴謹、性能更可靠、客觀性更強的旋回地層學(xué)分析方法，也需要引入新方法。

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

4 結(jié)語

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

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