喻小濤，胡光明，羅水亮，奉偉，李宇 (長江大學地球科學學院，湖北 武漢 430100)

張志鵬 (中國地質(zhì)大學(北京)能源學院，北京 100083)

李志威 (長江大學地球科學學院，湖北 武漢 430100)

1920年南斯拉夫?qū)W者米蘭柯維奇為解釋第四紀冰期成因提出米蘭柯維奇理論，該理論認為北半球高緯夏季太陽輻射變化信號被放大、傳輸進而使全球氣候發(fā)生周期性變化[1]。周期性的氣候變化通過沉積物的結(jié)構(gòu)、構(gòu)造及沉積層的厚度等形式被記錄下來，且米蘭柯維奇旋回的記錄已不限于第四紀，在顯生宙各個時代中均留下了痕跡[2～7]。如何將米蘭柯維奇周期準確提取出來，即如何從地質(zhì)信息中識別出地球公轉(zhuǎn)軌道三要素(偏心率、斜率和歲差)的變化周期，成為精確劃分地層、校準天文地質(zhì)年代和分析古氣候等課題所關(guān)注的問題之一。早先學者們常用傅里葉變換來提取米蘭柯維奇周期的頻率，但傅里葉變換丟掉了時間信息，無法判斷一個特定的信號是什么時間發(fā)生的，且相較于傅里葉變換中用正弦函數(shù)逼近信號，小波變換用不規(guī)則的小波函數(shù)來逼近尖銳信號的效果更好[8]，故小波變換逐漸成為識別米蘭柯維奇周期的常用工具。筆者以鄂爾多斯盆地中二疊統(tǒng)山西組1段～石盒子組8段(以下簡稱“山1段～盒8段”)的測井曲線為例，闡述了如何利用小波變換識別米蘭柯維奇周期。

1 資料選取

在二疊紀，鄂爾多斯盆地北緣的古亞洲洋消亡，大華北陸相沉積盆地雛形形成，鄂爾多斯盆地處于構(gòu)造相對穩(wěn)定的階段[9,10]，南北物源形成的沉積體系在盆地中南部匯聚,形成近東西向展布的淺湖[11]。

構(gòu)造與氣候是陸相沉積體系發(fā)育的兩大主要控制因素，構(gòu)造與古氣候信息都會被記錄在沉積地層中。鄂爾多斯盆地二疊紀構(gòu)造相對穩(wěn)定，意味著沉積記錄中的古氣候信息受到構(gòu)造(相當于噪音)的干擾較少且相對突出，有利于米蘭柯維奇周期的識別，這是選取山1段～盒8段作為研究對象的主要原因。另外，考慮到相對于沖積扇、河流、三角洲而言，處于盆地匯水區(qū)的淺湖地層最為完整，因此選取位于淺湖(匯水區(qū))的A299井、A1116井、A1118井、A1022井的自然伽馬曲線進行分析，確保了氣候信息的完整性，且?guī)r性相對穩(wěn)定，在一定程度上降低了巖性對米蘭柯維奇周期的影響。

2 求取各層頻率

以A299井為例，利用Matlab小波分析提取各層頻率的操作步驟分為兩大步。

1)準備數(shù)據(jù) ①先將A299井山1段～盒8段的自然伽馬(代碼GR)數(shù)據(jù)置于Excel表的一列中，第一行為A299GR，保存為A299GR.xls；②在Matlab主頁中選取“導入數(shù)據(jù)”，選擇A299GR.xls文件并打開；③進入導入窗口，選中自然伽馬整列數(shù)據(jù)，然后在視圖窗口的菜單中選擇“列向量”；④點擊菜單中“導入所選內(nèi)容”；⑤返回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；③根據(jù)a=2n(其中，a為尺度，n為層數(shù))，在系數(shù)圖中右擊選擇相應的層數(shù)(據(jù)選擇的尺度推算)，然后點擊New Coefficients Line，讀出該層的頻率(Frequency)；④重復步驟③選擇下一個層數(shù)并讀出相應的頻率，依次讀出d1～d10層的頻率(見表1)。

表1 各井不同尺度下的相對頻率

表1顯示各井的同層頻率是一致的，說明該地區(qū)各井所處位置的沉積物受到相同外部因素的影響。

為了便于與天文周期進行對比，將各層的頻率比上最小頻率(0.001)得到相應的值(見表2)。由于d9與d10對應的頻率一樣(表1)，說明d9層是軟件所能分析的極限了，d10層的數(shù)據(jù)不采用。

表2 不同尺度下的相對頻率與最小相對頻率的比值

3 求取理論周期與頻率

周期與頻率互為倒數(shù)關(guān)系，得到了周期值，也就得到了頻率值。地質(zhì)歷史中某一時期的偏心率、斜率和歲差周期的比值是相對固定的，如果能在研究區(qū)目的層段的系列周期(或頻率)中找出3個周期，其比值與理論米蘭柯維奇周期(或頻率)一致，即識別出了米蘭柯維奇周期(或頻率)。

偏心率周期405ka被認為是天文上最為穩(wěn)定的地球軌道參數(shù)，在過去幾億年是不變的[12]，因此選取405ka做為偏心率周期。地質(zhì)歷史中斜率和歲差并不是固定的，中生代以來的斜率和歲差可以根據(jù)公式計算[13,14 ]，但是更古老的地層，其中的干擾因素較多，需要通過其他方式求取。

Hinnov和Hilgen計[14]算出了過去250Ma(大致為中生代與古生代的分界)以來斜率和歲差的變化(表3)。在該基礎(chǔ)上，以時間為橫坐標，不同的天文周期為縱坐標，擬合出相應的函數(shù)關(guān)系式(圖1)。

表3 斜率和歲差周期在過去250Ma的變化(據(jù)文獻[14]，有修改)

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

圖1 地質(zhì)歷史時期的斜率周期與歲差周期

根據(jù)2014版的中國地層表和2016版的國際年代底層表，確定山1段～盒8段的地質(zhì)年齡界限為272.3±0.5Ma，取其中間值272.3Ma代入圖1中的各個函數(shù)關(guān)系式中，計算出相應的斜率周期和歲差周期分別為(單位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的偏心率、斜率和歲差周期轉(zhuǎn)換為頻率，計算各頻率與最小頻率的比值，取整數(shù)(表4)。將表4中頻率比值取整數(shù)后，與表2對比可知d9∶d6∶d5=1∶10∶21=e∶α1∶P2，即d9、d6、d5的頻率就是研究區(qū)要提取的天文周期所對應的頻率，分別代表e、α1、P2，據(jù)此可以提取與偏心率、斜率和歲差相對應的曲線。需要特別說明的是，表4中的頻率為真實頻率，而表1中的頻率則為范圍在0～1之間的相對頻率。因此，d9、d6、d5的頻率與e、α1、P2并不相等，只是它們之間的比值是相等的。

表4 目的層(272.3±0.5Ma)對應的周期及其相應的處理數(shù)據(jù)

圖2 A1022井的3個天文周期

所提取的偏心率周期為405ka，斜率周期為38.9ka，歲差周期為19.4ka。歲差的周期：斜率的周期：偏心率的周期=1∶2∶21。因此，對于某一段地層來說，其包含的歲差周期、斜率周期、偏心率周期的個數(shù)的比值應接近21∶10∶1，即若某一段地層中含有1個偏心率周期，則含有的歲差周期個數(shù)為21個，斜率周期個數(shù)為10個。如果該次研究提取的偏心率、斜率、歲差是正確的，那么就應該找到這樣的地層，使得地層中只要含有1個偏心率周期，就有對應的21個歲差周期和10個斜率周期。從圖2可以知道，A1022井基本符合要求，圖中紅線所框選的部分中，從下往上看包含1個偏心率周期，10個斜率周期，18個歲差周期，基本滿足該比值，從而驗證了所提取的天文周期的合理性。

5 結(jié)語

天文周期普遍存在于沉積地層中，利用Matlab中的小波變換工具從自然伽馬曲線中提取米蘭柯維奇旋回對應的周期，首先要將目的層段的自然伽馬曲線轉(zhuǎn)換為“.mat”格式，再使用Continuous Wavelet 1-D將其分解為10層，讀取各層的頻率，并計算各層頻率之間的比值，然后將該比值與該時代的理論頻率比值進行對比，從而識別出d9、d6、d5等對應的頻率和周期為研究區(qū)山1段～盒8段偏心率、斜率和歲差所對應的頻率和周期。前人在應用米蘭柯維奇周期時，往往應用的是偏心率，而該次研究成果是從一段測井曲線中提取相應的3個米蘭柯維奇周期(斜率、歲差、偏心率)，對于精確劃分地層，校準天文地質(zhì)年代具有重要的意義。

[參考文獻]

[1]吳瑞堂，張守信. 現(xiàn)代地層學[M]. 武漢：中國地質(zhì)大學出版社，1991：2～6.

[2]吳智勇. 米蘭柯維奇韻律層及其年代地層意義[J]. 地層學雜志，1995，19(2)：156～160.

[3]金之鈞，范國章，劉國臣. 一種地層精細定年的新方法[J]. 地球科學(中國地質(zhì)大學學報)，1999，24(4)：379～381.

[4]童金南，殷鴻福. 浙江長興煤山剖面Griesbachian期旋回地層研究[J]. 地層學雜志，1999，23(2)：130～134.

[5]郭剛，童金南，張世紅，等. 安徽巢湖早三疊世印度期旋回地層研究[J]. 中國科學D輯(地球科學)，2007，37(12)：1571～1578.

[6]龔一鳴，徐冉，湯中道，等. 廣西上泥盆統(tǒng)軌道旋回地層與牙形石帶的數(shù)字定年[J]. 中國科學D輯(地球科學)，2004，34 (7)：635～643.

[7]蘇德辰，李慶謀，羅光文，等. Fischer圖解及其在旋回層序研究中的應用——以北京西山張夏組為例[J]. 現(xiàn)代地質(zhì)，1995，9(3)：279～283.

[8]邵龍義，何志平，魯靜，等. 環(huán)渤海灣西部石炭系-二疊系層序地層及聚煤作用研究[M]. 北京：地質(zhì)出版社，2008：111～112.

[9]周安朝，賈炳文，馬美玲，等. 華北板塊北緣晚古生代火山事件沉積的全序列及其主要特征[J]. 地質(zhì)論評，2001，47(2)：175～182.

[10]楊俊杰. 鄂爾多斯盆地構(gòu)造演化與尤其分布規(guī)律[M]. 北京：石油工業(yè)出版社，2002：3～5.

[11]田景春，吳琦，王峰，等. 鄂爾多斯盆地下石盒子組盒8 段儲集砂體發(fā)育控制因素及沉積模式研究[J]. 巖石學報，2011，27(8)：2403～2411.

[12]黃春菊. 旋回地層學和天文年代學及其在中生代的研究現(xiàn)狀[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.