利用小波變換從測井曲線中提取米蘭柯維奇周期

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

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

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

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

1 資料選取

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

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

2 求取各層頻率

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

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

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

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

為了便于與天文周期進行對比，將各層的頻率比上最小頻率(0.001)得到相應(yīng)的值(見表2)。由于d9與d10對應(yīng)的頻率一樣(表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ǔ)上，以時間為橫坐標(biāo)，不同的天文周期為縱坐標(biāo)，擬合出相應(yīng)的函數(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)系式中，計算出相應(yīng)的斜率周期和歲差周期分別為(單位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ū)要提取的天文周期所對應(yīng)的頻率，分別代表e、α1、P2，據(jù)此可以提取與偏心率、斜率和歲差相對應(yīng)的曲線。需要特別說明的是，表4中的頻率為真實頻率，而表1中的頻率則為范圍在0～1之間的相對頻率。因此，d9、d6、d5的頻率與e、α1、P2并不相等，只是它們之間的比值是相等的。

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

圖2 A1022井的3個天文周期

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

5 結(jié)語

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

[參考文獻]

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

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

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

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

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

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

[7]蘇德辰，李慶謀，羅光文，等. Fischer圖解及其在旋回層序研究中的應(yīng)用——以北京西山張夏組為例[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]. 巖石學(xué)報，2011，27(8)：2403～2411.

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