基于測(cè)井?dāng)?shù)據(jù)的米氏旋回分析及浮動(dòng)天文年代標(biāo)尺的建立

徐敬領(lǐng), 霍家慶， 宋連騰， 昝靈

1 中國地質(zhì)大學(xué)(北京)地球物理與信息技術(shù)學(xué)院， 北京 100083 2 中國石油勘探開發(fā)研究院， 北京 100083 3 中國石油化工股份有限公司華東油氣分公司勘探開發(fā)研究院， 南京 210019

0 引言

米氏旋回即米蘭科維奇旋回，是指天文旋回導(dǎo)致氣候(日照量)的周期性變化進(jìn)而形成沉積韻律，最早由南斯拉夫?qū)W者M(jìn)ilankovitch(1941)在研究第四紀(jì)冰期旋回時(shí)提出，后來又被相關(guān)學(xué)者(Shackleton and Opdyke，1973；Hays et al.，1976；Berger, 1976, 1978; Berger and Loutre, 1991)通過大量的數(shù)據(jù)研究來證實(shí).隨著旋回地層學(xué)的不斷發(fā)展，米蘭科維奇旋回理論在地層劃分對(duì)比、地層年代框架建立以及古環(huán)境古氣候的研究中得到了廣泛應(yīng)用(李鳳杰等，2008；王燕等，2019；吳懷春等，2008；Puetz et al.，2016；Westphal et al.，2004；賈澤南等，2019).總結(jié)前人研究米氏旋回的分析方法主要可分為兩類，一類是巖性觀察法(Weedon，2003；龔一鳴等，2004)，主要通過對(duì)鉆井巖心或野外剖面的巖性巖相變化來觀察分析，識(shí)別其中的米氏旋回記錄.該方法識(shí)別米蘭科維奇旋回存在一定困難，需要研究人員有豐富的經(jīng)驗(yàn)與技能.另一類是時(shí)間序列分析方法，Weedon(2003)、鄭民等(2007)、謝灝辰等(2016)采用快速傅里葉變換對(duì)研究區(qū)測(cè)井?dāng)?shù)據(jù)進(jìn)行了頻譜分析，識(shí)別出地層中的米氏沉積旋回，為沉積作用的分析和高頻層序的劃分提供依據(jù).然而快速傅里葉變換的時(shí)域分辨力差，對(duì)于非平穩(wěn)的實(shí)際測(cè)井信號(hào)序列來說，其無法分析米氏旋回沉積厚度在深度域的變化特征.袁偉等(2016)、楊雨等(2021)、劉浩童(2020)將傅里葉變換和小波分析方法相結(jié)合，兩種方法的識(shí)別結(jié)果不僅可以相互對(duì)比驗(yàn)證，進(jìn)而分析出米氏旋回在深度域上的變化特征，并計(jì)算沉積速率，最終精細(xì)地劃分地層.然而該方法的研究重點(diǎn)還是在識(shí)別米氏旋回和劃分地層上，沒有進(jìn)一步研究地層的年代學(xué)信息和建立高精度、高分辨率的時(shí)間地層格架.

為了從時(shí)間尺度上解決地球科學(xué)多個(gè)領(lǐng)域的關(guān)鍵問題，專家學(xué)者們開始尋找引起地質(zhì)沉積周期變化的天文計(jì)時(shí)標(biāo)準(zhǔn)(汪品先，2006；吳懷春等，2011)，于是米氏旋回理論逐漸在浮動(dòng)天文年代標(biāo)尺的建立中得到應(yīng)用.其中，Prokoph等(2001)利用米氏旋回理論，實(shí)現(xiàn)加拿大西部白堊紀(jì)Cenomanian-Turonian期界線附近地層連續(xù)的時(shí)間標(biāo)定，估計(jì)了第二次大洋缺氧事件(OAE2)的持續(xù)時(shí)間約為320 ka.吳懷春等(2008)通過對(duì)松科1井南井青山口組的測(cè)井資料進(jìn)行米氏旋回性分析，建立了青山口組浮動(dòng)天文年代標(biāo)尺，估計(jì)出青山口組1段底部發(fā)生的湖泊缺氧事件的持續(xù)時(shí)間約為250 ka，為后續(xù)研究青山口組沉積期的各類地質(zhì)事件提供了年代學(xué)依據(jù).魏小松等(2018)對(duì)北部灣盆地流沙港組沉積地層的自然伽馬(GR)曲線進(jìn)行米氏旋回性分析，并用鋯石提供年代(時(shí)間錨點(diǎn))，最終建立該地區(qū)的浮動(dòng)天文年代標(biāo)尺，為地層高頻層序劃分及提高油氣勘探精度提供思路.鄧秀芹等(2013)、張輝等(2014)、Wang等(2014)、張文等(2017)主要采用次離子質(zhì)譜(SIMS)和激光剝蝕電感耦合等離子體質(zhì)譜(LA-ICP-MS)等實(shí)驗(yàn)方法對(duì)鄂爾多斯盆地長(zhǎng)7段開展了鋯石定年分析，得到長(zhǎng)7段年齡數(shù)據(jù)誤差約為2.0 Ma；Zhu等(2019)通過鋯石U-Pb ID-TIMS測(cè)年獲得長(zhǎng)7段相對(duì)精確的年齡為241.06±0.12Ma至241.558±0.093 Ma，數(shù)據(jù)誤差范圍約0.1 Ma，Chu等(2020)在此基礎(chǔ)上研究了古氣候、古環(huán)境與旋回，而基于測(cè)井?dāng)?shù)據(jù)研究長(zhǎng)7段米氏旋回及沉積速率存在欠缺.以上天文年代標(biāo)尺的建立有助于對(duì)研究區(qū)地質(zhì)事件、古環(huán)境和古氣候事件的持續(xù)時(shí)間做出精確估計(jì)，但都沒有明確分析方法的原理、方法的可行性及從模擬的角度對(duì)米氏旋回分析的可行性做出探討.

本文基于鄂爾多斯盆地B159井的測(cè)井?dāng)?shù)據(jù)，從傅里葉變換與小波深頻分析的原理入手，通過數(shù)值模擬手段研究傅里葉變換和小波深頻分析識(shí)別米氏旋回的可行性，再以B159井長(zhǎng)7段的GR數(shù)據(jù)為研究對(duì)象，綜合傅里葉分析和小波深頻分析方法，識(shí)別出長(zhǎng)7段地層中的米氏旋回，研究浮動(dòng)天文年代標(biāo)尺建立的方法，最終建立長(zhǎng)7段地層的浮動(dòng)天文年代標(biāo)尺，為建立地層的時(shí)間格架提供一種高效的方法與思路.

1 米氏旋回頻譜分析方法

1.1 米氏旋回?cái)?shù)值模擬與頻譜分析的理論基礎(chǔ)

自然界中具有周期性的運(yùn)動(dòng)是普遍存在的.地層沉積的過程就是一個(gè)復(fù)雜的周期運(yùn)動(dòng)，它是由多個(gè)單周期的沉積運(yùn)動(dòng)疊加而成，每種周期的簡(jiǎn)單運(yùn)動(dòng)代表了這段地質(zhì)歷史時(shí)期所發(fā)生的沉積事件，它具有周期性、疊置性特點(diǎn)，表現(xiàn)為沉積物的結(jié)構(gòu)、構(gòu)造、巖性及巖相的周期性變化(鄭民等，2007).測(cè)井信號(hào)本身也是一種由多信號(hào)疊加的復(fù)合信號(hào)，它記錄了地層沉積的物質(zhì)變化過程，能夠敏感、連續(xù)地反映所測(cè)地層的旋回性、周期性等沉積特征，可以通過一定的方法技術(shù)提取地層旋回性特征，找出主要的周期，恢復(fù)地層的疊置韻律特征.

頻譜分析是研究周期性事件的一種常見數(shù)學(xué)方法，其原理是利用傅里葉變換，把一條復(fù)雜的周期信號(hào)分解為一系列不同頻率的正弦波或余弦波，找出其中振幅強(qiáng)度大且數(shù)學(xué)上呈現(xiàn)與米氏旋回周期相近的頻率(胡麗瑩等，2011).在實(shí)際工作過程中，采樣是離散和有限的，因此在實(shí)際數(shù)據(jù)處理過程中都采用離散傅里葉變換，對(duì)于一段N項(xiàng)離散時(shí)間序列Xn，其離散傅里葉變換公式為

(1)

式中，Xk為頻譜值，i為虛數(shù)單位，k=0，1，2,…，N-1.

而快速傅里葉變換(FFT)則是離散傅里葉變換(DFT)的快速算法，能迅速提高DFT的運(yùn)算速度，使運(yùn)算量減少兩個(gè)數(shù)量級(jí).

為了論證快速傅里葉變換識(shí)別米氏旋回的有效性，對(duì)米氏旋回做數(shù)值模擬研究.根據(jù)謝灝辰等(2016)、袁偉等(2016)用Berger法的計(jì)算結(jié)果：晚三疊世長(zhǎng)短偏心率周期為404和123 ka，地軸傾角長(zhǎng)短周期分別為44.3和35.4 ka，歲差長(zhǎng)短周期分別為21.3和17.8 ka，首先模擬出周期分別為6個(gè)天文周期的正弦單信號(hào)曲線(如圖4第二、三、四、五列)，然后將它們疊加，得到復(fù)合信號(hào)，該復(fù)合信號(hào)則可代表周期性沉積疊置地層的測(cè)量信號(hào)，最后對(duì)復(fù)合信號(hào)做FFT，如圖1所示.圖1左側(cè)是6個(gè)天文周期的復(fù)合信號(hào)曲線，中間是對(duì)復(fù)合信號(hào)做FFT的結(jié)果，右側(cè)是對(duì)復(fù)合信號(hào)做多窗口頻譜分析(MTM)的結(jié)果.

圖1 模擬的6個(gè)天文周期復(fù)合信號(hào)與頻譜分析結(jié)果(a) 模擬的6個(gè)天文周期復(fù)合信號(hào)； (b) 復(fù)合信號(hào)的FFT分析結(jié)果； (c) 復(fù)合信號(hào)的MTM分析結(jié)果.Fig.1 The simulated composite signal of six astronomical periods and spectrum analysis results(a) The simulated composite signal of 6 astronomical periods; (b) The FFT analysis result of the composite signal; (c) The MTM analysis result of the composite signal.

分析圖1中的復(fù)合信號(hào)曲線及FFT與MTM頻譜分析結(jié)果得出：(1)6個(gè)天文周期組合的復(fù)合信號(hào)是比較復(fù)雜的，無法直觀看出其包含的周期性信息；(2)通過對(duì)復(fù)合信號(hào)做FFT與MTM頻譜分析，可以識(shí)別出復(fù)合信號(hào)所包含的有效頻率信息，如：0.006、0.02、0.056、0.07、0.118和0.14個(gè)旋回/m，F(xiàn)FT與MTM對(duì)復(fù)合信號(hào)的分析結(jié)果一致，特別是6個(gè)尖峰頻率是一樣的，沒有差異；(3)FFT頻譜分析中的低頻信號(hào)代表了長(zhǎng)周期沉積的時(shí)間，如頻率值0.006對(duì)應(yīng)長(zhǎng)偏心率時(shí)間404 ka，高頻信號(hào)代表短周期沉積的時(shí)間，如頻率值0.14對(duì)應(yīng)短歲差時(shí)間17.8 ka；(4)該復(fù)合信號(hào)FFT頻譜分析中含有6個(gè)主要頻率信號(hào)，對(duì)應(yīng)頻率分別為0.006、0.02、0.056、0.07、0.118和0.14個(gè)旋回/m，根據(jù)頻率倒數(shù)乘以采樣頻率可得旋回厚度(或周期時(shí)間)的公式，計(jì)算其對(duì)應(yīng)的旋回厚度比為1∶0.3∶0.11∶0.086∶0.051∶0.043，這與6種天文周期的比值1∶0.3∶0.11∶0.088∶0.053∶0.044完全吻合.綜上，通過對(duì)6個(gè)軌道周期的模擬與頻譜分析，其結(jié)果不僅證明了FFT頻譜分析米氏旋回的有效性與精確性，更為頻譜分析識(shí)別劃分地層旋回奠定了理論基礎(chǔ).

1.2 基于實(shí)際測(cè)井?dāng)?shù)據(jù)的米氏旋回頻譜分析方法

基于米氏旋回的數(shù)值模擬與頻譜分析結(jié)果，選擇鄂爾多斯盆地三疊系長(zhǎng)7段作為研究對(duì)象，以B159井長(zhǎng)7段為例，該段厚度為124 m，是典型的陸相深湖—半深湖沉積，其構(gòu)造穩(wěn)定、連續(xù)性好，未發(fā)生沉積間斷，且該段地層的巖性以泥頁巖為主，而泥頁巖中的黏土礦物對(duì)放射性元素吸附能力強(qiáng)，所以自然伽馬曲線能夠反映地層中黏土礦物含量的變化，從而敏感地識(shí)別沉積過程中受天文周期影響的米氏旋回(李鳳杰等，2004；石巨業(yè)等，2017；Boulila et al.，2021).

選取B159井長(zhǎng)7段的自然伽馬曲線，采樣間距Δ為0.125 m，提取該井段地層米氏旋回的頻譜分析方法步驟為：(1)對(duì)GR數(shù)據(jù)進(jìn)行l(wèi)g(x+1)標(biāo)準(zhǔn)化預(yù)處理，再對(duì)預(yù)處理后的數(shù)據(jù)進(jìn)行快速傅里葉變換，得到頻譜圖，如圖2所示；(2)計(jì)算晚三疊世天文周期時(shí)間比為404 ka：123 ka：44.3 ka：35.4 ka：21.3 ka：17.8 ka，即天文周期比為1∶0.3：0.11∶0.088∶0.053∶0.044；(3)計(jì)算圖2中有效頻率所對(duì)應(yīng)周期的比值，再與步驟(2)中的天文周期比作對(duì)比，結(jié)合圖1模擬信號(hào)得出的6個(gè)米氏旋回對(duì)應(yīng)的頻率值，直接可找出與其相近的5個(gè)有效頻率：0.00605、0.02015、0.0564、0.0705和0.1128個(gè)旋回/m，分別命名為A、B、C、D、E5個(gè)頻率，且FFT頻譜分析結(jié)果與MTM頻譜分析結(jié)果一致；(4)通過周期旋回厚度或周期時(shí)間計(jì)算公式d=(1/f)×fs，計(jì)算得到A、B、C、D、E頻率的旋回厚度分別為20.68 m、6.2 m、2.2 m、1.77 m和1.11 m，對(duì)應(yīng)的旋回厚度比為1∶0.3∶0.11∶0.086∶0.054.其中，f為頻率，fs為采樣頻率，值為8；(5)由于該段地層為深湖沉積環(huán)境、沉積穩(wěn)定，其沉積物供給速率相對(duì)穩(wěn)定，則識(shí)別出A頻率的旋回厚度20.68 m對(duì)應(yīng)為404 ka的長(zhǎng)偏心率周期，B頻率的旋回厚度6.2 m對(duì)應(yīng)為123 ka的短偏心率周期，C頻率的旋回厚度2.2 m對(duì)應(yīng)為44.3 ka的地軸傾角長(zhǎng)周期，D頻率的旋回厚度1.77 m對(duì)應(yīng)為35.4 ka的地軸傾角短周期，E頻率的旋回厚度1.11 m對(duì)應(yīng)為21.3 ka的長(zhǎng)歲差周期；(6)根據(jù)上述米氏旋回的識(shí)別結(jié)果，即長(zhǎng)7段包含的天文時(shí)間沉積周期，可用公式(2)計(jì)算出B159井長(zhǎng)7段平均沉積速率為5.05 cm·ka-1.

×100=5.05 cm·ka-1,(2)

式中，Vp為頻譜分析計(jì)算的平均沉積速率，dA、dB、dC、dD、dE分別為A、B、C、D、E頻率的旋回厚度.

1.3 實(shí)際測(cè)井?dāng)?shù)據(jù)識(shí)別米氏旋回與計(jì)算沉積速率方法的應(yīng)用效果分析

根據(jù)前文對(duì)B159井GR數(shù)據(jù)的頻譜分析結(jié)果，分析該方法在識(shí)別米氏旋回與計(jì)算沉積速率的應(yīng)用效果：(1)通過FFT頻譜分析，可以將復(fù)雜的自然伽馬數(shù)據(jù)變換到頻率域，分析其頻率信息，并提取能夠反映米氏旋回周期的有效頻率，從而更加有效地分析地層疊置旋回序列與時(shí)間周期；(2)根據(jù)FFT頻譜分析及天文周期比識(shí)別出B159井長(zhǎng)7段所包含的5個(gè)米氏旋回有效頻率為0.00605、0.02015、0.0564、0.0705和0.1128個(gè)旋回/m(圖2)；(3)計(jì)算5個(gè)有效頻率對(duì)應(yīng)的米氏旋回周期比為1∶0.3∶0.11∶0.086：0.054，其與天文周期比為1∶0.3∶0.11∶0.088∶0.053，非常接近，且誤差不超過2.12%；(4)計(jì)算每個(gè)有效頻率對(duì)應(yīng)天文周期的旋回厚度分別為20.66 m、6.2 m、2.2 m、1.77 m和1.11 m，如表1所示，并可根據(jù)每個(gè)旋回厚度與旋回周期計(jì)算出每個(gè)天文周期的沉積速率，得出B159井長(zhǎng)7段平均沉積速率Vp為5.05 cm·ka-1，這與Zhu等(2019)通過鋯石U-Pb ID-TIMS測(cè)年獲得長(zhǎng)7段泥頁巖沉積速率(5 cm·ka-1)非常吻合.

圖2 Bai159井FFT頻譜分析(a) 測(cè)井曲線； (b) GR的FFT分析結(jié)果； (c) GR的MTM分析結(jié)果.Fig.2 FFT spectrum analysis of Well Bai159(a) The logging curve; (b) The FFT analysis result of GR; (c) The MTM analysis result of GR.

表1 鄂爾多斯盆地延長(zhǎng)組長(zhǎng)7段各天文周期的沉積速率計(jì)算結(jié)果Table 1 Calculation results of deposition rate in each astronomical period of Chang 7 Member of Yanchang Formation in Ordos Basin

2 米氏旋回小波深頻分析方法

2.1 米氏旋回?cái)?shù)值模擬與小波深頻分析的理論基礎(chǔ)

(3)

對(duì)測(cè)井信號(hào)做小波深頻分析提取旋回信息的方法步驟如下：

(1)選取小波基函數(shù)，即母小波，一般選morlet小波，其不但具有非正交性而且還是由Gaussian指數(shù)復(fù)值小波，可以得到連續(xù)平滑的振幅和相位；

(2)選取一定尺度a的母小波與信號(hào)f(x)的起始點(diǎn)做卷積；

(3)計(jì)算得到小波變換后的系數(shù)CWTf(a,b)，它表示尺度a母小波與該段信號(hào)的相關(guān)程度.CWTf(a,b)越大，它們?cè)较嗨疲?/p>

(4)逐步按一定的步長(zhǎng)向右平移小波進(jìn)行計(jì)算，重復(fù)步驟(2)、(3)直到所有信號(hào)被覆蓋，如圖3所示，得到一組該尺度下的小波系數(shù)矩陣；

(5)然后，伸縮變換小波尺度因子a，重復(fù)步驟(2)—(4)；

(6)最后得到所有尺度a對(duì)應(yīng)的小波系數(shù)CWTf(a,b)矩陣，即不同尺度(頻率)的分頻信息與能量旋回結(jié)構(gòu).

圖3 小波深頻分析變換步驟Fig.3 Wavelet deep-frequency analysis and transformation steps

從小波深頻分析變換的步驟中可以看出，小波深頻分析方法可以將一維深度域的測(cè)井?dāng)?shù)據(jù)轉(zhuǎn)換為二維深度-尺度域的小波系數(shù)矩陣Ca×b.尺度a小波與信號(hào)的相似度越高，小波系數(shù)值越大，說明該深度位置很有可能存在尺度為a的沉積旋回.因此，可以通過小波系數(shù)值的大小變化，找出各種尺度的沉積旋回，再從中識(shí)別出尺度比(頻率比)符合天文周期比的米氏旋回(徐敬領(lǐng)等，2009；張?zhí)沟龋?017；閆建平等，2017).

為了論證小波深頻分析識(shí)別米氏旋回的有效性，對(duì)米氏旋回做數(shù)值模擬研究：(1)模擬部分天文周期(長(zhǎng)偏心率、短偏心率、短地軸傾角和短歲差周期)的單信號(hào)曲線，并將它們組合成一個(gè)復(fù)合信號(hào)(圖4，第二—六列)，該復(fù)合信號(hào)則可代表周期性沉積疊置地層的測(cè)量信號(hào)；(2)對(duì)復(fù)合信號(hào)做小波變換，得到小波深頻能量譜圖(圖4，最后一列)；(3)對(duì)不同尺度a的小波系數(shù)做模平均值，構(gòu)建不同尺度a的小波系數(shù)模平均值曲線(圖5).

圖4 模擬的天文周期單信號(hào)、復(fù)合信號(hào)及小波深頻分析能量譜圖(a) 模擬的天文周期單信號(hào)與復(fù)合信號(hào)； (b) 復(fù)合信號(hào)的小波深頻分析能量譜圖.Fig.4 Simulated astronomical period single signal, composite signal and wavelet deep-frequency analysis energy spectrum(a) Simulated astronomical period single signal and composite signal; (b) The wavelet deep-frequency analysis energy spectrogram of the composite signal.

分析圖4中模擬的復(fù)合信號(hào)曲線及小波深頻分析結(jié)果可知：(1)6個(gè)天文周期組合的復(fù)合信號(hào)是比較復(fù)雜的，無法直觀看出其包含的周期性信息；(2)小波深頻能量譜圖中(圖4b)橫坐標(biāo)為尺度，從左到右尺度逐漸增大，小尺度代表了短周期沉積時(shí)間的地層疊置，即短周期沉積地層疊置序列，對(duì)應(yīng)能量團(tuán)個(gè)數(shù)較多，大尺度代表長(zhǎng)周期沉積時(shí)間的地層疊置，即長(zhǎng)周期沉積地層疊置序列，對(duì)應(yīng)能量團(tuán)個(gè)數(shù)較少；(3)小波深頻能量譜圖縱坐標(biāo)指示了能量團(tuán)代表的沉積地層疊置旋回邊界位置及疊置旋回個(gè)數(shù)信息；(4)小波深頻能量譜圖(圖4b)從左到右顯示有4種尺度的能量團(tuán)，分別對(duì)應(yīng)4種天文周期的單信號(hào)信息，從左到右尺度逐漸增大，能量團(tuán)也越大，每列能量團(tuán)大小一樣，對(duì)應(yīng)尺度也一樣.如短地軸傾角單頻信號(hào)曲線的每個(gè)波峰和波谷對(duì)應(yīng)著小波分析能量團(tuán)(圖4b)中從左數(shù)第二列每個(gè)能量團(tuán)的中心，即波峰和波谷分別對(duì)應(yīng)一個(gè)黃紅色的能量團(tuán)，因此兩個(gè)相鄰能量團(tuán)中心的距離對(duì)應(yīng)米氏周期的半旋回，3個(gè)能量團(tuán)為米氏周期的一個(gè)完整周期旋回；(5)根據(jù)波峰、波谷與小波能量團(tuán)的對(duì)應(yīng)關(guān)系，對(duì)小波深頻能量譜圖中能量團(tuán)的個(gè)數(shù)進(jìn)行數(shù)數(shù)，小波深頻能量譜圖中共包含0.5個(gè)長(zhǎng)偏心率周期的旋回，2.5個(gè)短偏心率周期的旋回，8.5個(gè)短地軸傾角周期的旋回，17.5個(gè)短斜率周期的旋回，與模擬的各天文周期旋回信號(hào)的個(gè)數(shù)是完全一致的，說明小波深頻分析不僅可以識(shí)別不同級(jí)次(尺度)的周期疊置旋回，而且可以計(jì)算同一級(jí)次與不同級(jí)次(尺度)的周期旋回個(gè)數(shù).

圖5為小波系數(shù)模平均值曲線，主要是建立尺度與周期旋回級(jí)次的關(guān)系.圖5中黑色圓圈為小波系數(shù)模平均值的局部極大值，其橫坐標(biāo)對(duì)應(yīng)最優(yōu)尺度為134、40和15，比值為1∶0.30∶0.11，與3種天文周期的比值1∶0.30∶0.11完全吻合，即對(duì)該模擬信號(hào)進(jìn)行小波深頻分析時(shí)，134尺度的小波系數(shù)或小波能量譜團(tuán)對(duì)應(yīng)長(zhǎng)偏心率周期，40尺度的小波系數(shù)或小波能量譜團(tuán)對(duì)應(yīng)短偏心率周期，15尺度的小波系數(shù)或小波能量譜團(tuán)對(duì)應(yīng)長(zhǎng)地軸傾角周期，可通過其中某個(gè)尺度的小波系數(shù)或小波能量譜團(tuán)識(shí)別對(duì)應(yīng)的地層周期疊置旋回.綜上，通過對(duì)4個(gè)天文周期的模擬與小波深頻分析，其結(jié)果不僅證明了小波深頻分析米氏旋回的有效性與精確性，更為小波深頻分析識(shí)別劃分地層旋回奠定了理論基礎(chǔ).

圖5 復(fù)合信號(hào)小波分析的最優(yōu)尺度選擇(黑圈是局部極大值、為最優(yōu)尺度的備選點(diǎn))Fig.5 Optimal scale selection for wavelet analysis of composite signals (The black circles are local maxima, which can be used as candidate points for the optimal scale)

2.2 基于實(shí)際測(cè)井?dāng)?shù)據(jù)的米氏旋回小波深頻分析方法

基于米氏旋回的數(shù)值模擬與小波深頻分析結(jié)果，優(yōu)選鄂爾多斯盆地三疊系長(zhǎng)7段作為研究對(duì)象，并取B159井長(zhǎng)7段的自然伽馬曲線作為分析數(shù)據(jù)，采樣間隔Δ為0.125 m，提取該井段地層米氏旋回的小波分析方法步驟為：(1)根據(jù)B159井的分層數(shù)據(jù)，選用B159井長(zhǎng)7段的GR數(shù)據(jù)，如圖6第二列所示；(2)對(duì)該段GR數(shù)據(jù)進(jìn)行l(wèi)g(x+1)標(biāo)準(zhǔn)化預(yù)處理，然后選取Morlet小波(中心頻率Fc=0.8125)作為小波基函數(shù)進(jìn)行小波深頻分析，得到小波深頻能量譜圖，并做了鏡像處理，如圖6第三列(圖6b)，橫坐標(biāo)為尺度；(3)計(jì)算每個(gè)尺度下小波系數(shù)的模平均值，構(gòu)建不同尺度的小波系數(shù)模平均值曲線，如圖6第四列(圖6c)，局部極大值處為最優(yōu)尺度；(4)根據(jù)頻譜分析法計(jì)算的長(zhǎng)7段平均沉積速率Vp=5.05 cm·ka-1，以及公式(4)、(5)、(6)得到尺度a=(Fc×Vp×T)/(100×Δ)，將6種天文周期T的值代入即可計(jì)算相應(yīng)尺度值，根據(jù)計(jì)算結(jié)果推測(cè)各個(gè)天文周期對(duì)應(yīng)的最大尺度值不超過134；(5)根據(jù)圖6c小波系數(shù)模平均值曲線得到大小不超過134的四個(gè)最優(yōu)尺度：12、15、40和134，再參考小波能量譜圖篩選出符合天文周期比的尺度134、40，其比值為1∶0.3，因此推斷尺度134和40分別對(duì)應(yīng)長(zhǎng)偏心率和短偏心率周期；(6)通過尺度值即可根據(jù)式(4)、(5)、(6)計(jì)算相應(yīng)天文周期的頻率、沉積厚度及沉積速率:

圖6 B159井長(zhǎng)7段GR曲線及其小波深頻分析結(jié)果(a) B159井長(zhǎng)7段GR曲線; (b) GR的小波深頻分析能量譜圖; (c) GR的小波分析最優(yōu)尺度.Fig.6 GR curve and wavelet deep-frequency analysis results of Chang 7 Member of Well B159(a) GR curve of Chang 7 Member of Well B159; (b) The middle picture is the energy spectrum of GR wavelet deep-frequency analysis; (c) The optimal scale of GR wavelet analysis.

(4)

(5)

(6)

其中，f為該尺度下天文周期的頻率，個(gè)旋回/m；d為該尺度天文周期內(nèi)的地層沉積厚度，m；v為該尺度天文周期內(nèi)地層的沉積速率， cm·ka-1；T為沉積時(shí)間，ka.

最后，根據(jù)公式(4)和(5)得到該尺度天文周期內(nèi)地層沉積厚度：d=(a×Δ)/Fc，計(jì)算出長(zhǎng)偏心率和短偏心率周期的沉積厚度分別為20.62 m和6.2 m；(8)根據(jù)計(jì)算的沉積厚度及長(zhǎng)偏心率和短偏心率的周期，可用公式(6)計(jì)算出B159井長(zhǎng)7段平均沉積速率為5.07 cm·ka-1:

(7)

式中，Vp為小波深度分析計(jì)算的平均沉積速率，dA、dB分別為長(zhǎng)偏心率和短偏心率周期的沉積厚度.

2.3 實(shí)際測(cè)井?dāng)?shù)據(jù)識(shí)別米氏旋回與沉積速率的應(yīng)用效果分析

根據(jù)圖6 B159井GR數(shù)據(jù)的小波深頻分析結(jié)果，分析該方法在識(shí)別米氏旋回與計(jì)算沉積速率的應(yīng)用效果：(1)通過小波分析，可以將復(fù)雜的自然伽馬數(shù)據(jù)變換到二維深度-尺度域的小波系數(shù)矩陣；(2)根據(jù)小波分析能夠識(shí)別出地層中的優(yōu)勢(shì)旋回，即包含米氏周期的旋回，并能提供這些優(yōu)勢(shì)旋回隨深度的變化信息及對(duì)應(yīng)尺度；(3)根據(jù)小波深頻分析能量譜圖及天文周期比識(shí)別出B159井長(zhǎng)7段所包含的2個(gè)米氏旋回最優(yōu)尺度為134和40；(4)計(jì)算2個(gè)最優(yōu)尺度對(duì)應(yīng)的米氏旋回周期比為1∶0.2985，其與天文周期比1∶0.3非常接近，且誤差不超過0.5%，即可推斷出2個(gè)最優(yōu)尺度的米氏旋回分別對(duì)應(yīng)長(zhǎng)、短偏心率；(5)識(shí)別出B159井長(zhǎng)7段含有6個(gè)長(zhǎng)偏心率旋回，20個(gè)短偏心率旋回，且尺度為134的1個(gè)長(zhǎng)偏心率旋回約包含尺度為40的3個(gè)短偏心率旋回；(6)小波深頻分析能量譜圖中，尺度等于134的能量譜團(tuán)旋回厚度與計(jì)算的長(zhǎng)偏心率周期沉積厚度一致，尺度等于40的能量譜團(tuán)旋回厚度與計(jì)算的短偏心率周期沉積厚度也一致.如圖6，在2010～2030.6 m深度段，尺度為134的能量譜團(tuán)指示了該段含有1個(gè)長(zhǎng)偏心率旋回，則一個(gè)長(zhǎng)偏心率旋回的厚度約為20.6 m，與前文公式(5)、(6)計(jì)算的長(zhǎng)偏心率周期沉積厚度20.62 m吻合；(7)最后根據(jù)旋回周期與旋回厚度，計(jì)算B159井長(zhǎng)7段平均沉積速率Vp為5.07 cm·ka-1，這與頻譜分析法計(jì)算的平均沉積速率Vp(5.05 cm·ka-1)比較相近，且與Zhu等(2019)通過鋯石U-Pb ID-TIMS測(cè)年獲得長(zhǎng)7段泥頁巖沉積速率(5 cm·ka-1)非常吻合.

3 浮動(dòng)天文年代標(biāo)尺的建立

為了研究地層沉積、地質(zhì)事件的持續(xù)時(shí)間，以及古氣候、古環(huán)境與地質(zhì)事件的關(guān)系，需要建立高精度且具有相對(duì)時(shí)間概念的浮動(dòng)天文年代標(biāo)尺(吳懷春等，2011).由于三疊紀(jì)時(shí)期火星和地球軌道的諧振作用發(fā)生了變化，導(dǎo)致短偏心率周期產(chǎn)生一定的改變，而長(zhǎng)偏心率周期在該地質(zhì)歷史時(shí)期則具有較強(qiáng)的穩(wěn)定性(Laskar，2004).因此，研究對(duì)象為鄂爾多斯盆地三疊系長(zhǎng)7段時(shí)，以長(zhǎng)偏心率周期為主、以其它天文周期為輔來重構(gòu)各尺度下的測(cè)井曲線，作為建立浮動(dòng)天文年代標(biāo)尺的依據(jù).

三疊紀(jì)長(zhǎng)7段浮動(dòng)天文年代標(biāo)尺建立過程為：(1)根據(jù)B159井長(zhǎng)7段的小波分析結(jié)果，得到反映米氏旋回周期的最優(yōu)尺度為134和40，分別對(duì)應(yīng)長(zhǎng)偏心率周期和短偏心率周期；(2)提取出尺度134和40對(duì)應(yīng)的小波系數(shù)，小波系數(shù)隨深度的變化曲線即為該尺度下的重構(gòu)測(cè)井曲線，曲線中兩個(gè)相鄰極大值點(diǎn)的距離表示該尺度旋回的一個(gè)周期；(3)編寫算法剔除重構(gòu)曲線中具有干擾性的無效極大值點(diǎn)，僅標(biāo)記有效極大值點(diǎn)并計(jì)算旋回的個(gè)數(shù)，如圖7所示.

圖7 小波深頻分析重構(gòu)米氏旋回曲線(紅圈為有效極大值點(diǎn)，兩個(gè)相鄰紅圈間距為一個(gè)周期)Fig.7 Reconstructed Milankovitch cycle curves of wavelet deep-frequency analysis (The red circle is the effective maximum point, and the distance between two adjacent red circles is one period)

(4)Zhu等(2019)采用鋯石U-Pb ID-TIMS測(cè)年技術(shù)獲得長(zhǎng)7段頂部烴源巖年齡約241.558 Ma，根據(jù)測(cè)得頂部巖石年齡和尺度134的長(zhǎng)偏心率重構(gòu)曲線得到的該地層沉積時(shí)間，可得出地層的時(shí)間序列，即地質(zhì)歷史時(shí)間，如圖8第八列所示；(5)以尺度134的長(zhǎng)偏心率重構(gòu)曲線為依據(jù)，繪制長(zhǎng)偏心率周期旋回的地層厚度及沉積速率隨時(shí)間的變化曲線，如圖8最后一列，圖中曲線上標(biāo)注的數(shù)字為深度點(diǎn)，兩個(gè)深度點(diǎn)的距離為一個(gè)長(zhǎng)偏心率周期時(shí)間，曲線斜率為沉積速率；(6)綜合測(cè)井曲線、小波深頻色譜圖(圖8第五列)、不同尺度對(duì)應(yīng)的重構(gòu)曲線(圖8第六/七列)、地層歷史時(shí)間(圖8第八列)、旋回地層厚度及沉積速率變化(圖8第九列)，得出綜合分析圖.綜上，通過小波深頻分析，從B159井長(zhǎng)7段GR數(shù)據(jù)中提取出404 ka的長(zhǎng)偏心率周期和123 ka的短偏心率周期兩條曲線，將這兩條曲線作為調(diào)諧曲線建立了長(zhǎng)7段的浮動(dòng)天文年代標(biāo)尺.

圖8 B159井長(zhǎng)7段浮動(dòng)天文年代標(biāo)尺第五列為小波深頻分析能量譜圖;第六列為重構(gòu)的143尺度對(duì)應(yīng)長(zhǎng)偏心率周期曲線；第七列為重構(gòu)的40尺度對(duì)應(yīng)短偏心率周期曲線；第八列為確定的年代；第九列黑色橫杠為長(zhǎng)偏心率旋回地層的深度點(diǎn)與年代點(diǎn)，紅色斜線為沉積速率曲線，斜率越大、沉積速率越大.Fig.8 Floating astronomical data scale of Chang 7 member of Well B159The fifth column is the wavelet deep frequency analysis energy spectrum; The sixth column is the reconstructed 143-scale long eccentricity period curve; The seventh column is the reconstructed 40-scale short eccentricity period curve; The 8th column is the determined age; The black horizontal bar in the ninth column is the depth point and age point of the formation with long eccentricity period, and the red oblique line is the deposition rate curve, the greater the slope, the greater the deposition rate.

根據(jù)圖8各列信息的對(duì)比，分析該方法在建立浮動(dòng)天文年代標(biāo)尺的應(yīng)用效果：(1)404 ka長(zhǎng)偏心率曲線和123 ka短偏心率曲線的波峰與波谷位置都與小波深頻色譜圖中的能量團(tuán)有著很好的對(duì)應(yīng)；(2)研究出尺度與典型米氏周期(長(zhǎng)偏心率、短偏心率)的關(guān)系，即尺度134和40分別對(duì)應(yīng)長(zhǎng)、短偏心率，根據(jù)關(guān)系與小波系數(shù)，計(jì)算出B159井長(zhǎng)7段共記錄了6個(gè)長(zhǎng)偏心率旋回，20個(gè)短偏心率旋回；(3)鄂爾多斯盆地三疊系長(zhǎng)7段長(zhǎng)偏心率比較穩(wěn)定(Laskar，2004)，以404 ka長(zhǎng)偏心率曲線為基準(zhǔn)，計(jì)算出B159井長(zhǎng)7段的沉積時(shí)限約為6×404 ka=2424 ka，平均沉積速率V=124 m/2424 ka=5.1 cm·ka-1，這與頻譜分析、小波分析計(jì)算結(jié)果十分吻合，并與Zhu等(2019)前人計(jì)算結(jié)果也完全吻合，論證了以長(zhǎng)偏心率曲線精細(xì)標(biāo)定地層相對(duì)時(shí)間的可靠性；(4)長(zhǎng)7段地層短偏心率周期雖然有不斷變化，但長(zhǎng)偏心率周期與短偏心率周期的旋回個(gè)數(shù)比為20∶6=1∶0.3，與天文周期比完全吻合，說明短偏心率仍能夠作為精細(xì)標(biāo)定長(zhǎng)7段地層相對(duì)時(shí)間的參考依據(jù)；(5)圖8最后一列黑色橫杠除了代表深度點(diǎn)，同時(shí)也表示了疊置層的沉積時(shí)間點(diǎn)，其長(zhǎng)短代表了地層的相對(duì)歷史時(shí)間，其長(zhǎng)度越短表示地層歷史時(shí)間相對(duì)越新，深度由深到淺對(duì)應(yīng)地層時(shí)間由老到新，相鄰黑色橫杠的垂直距離為旋回的沉積厚度，因此紅色線條的斜率越大，相應(yīng)層段的沉積速率越大；(6)將第二—四列測(cè)井曲線與最后一列沉積速率作對(duì)比，發(fā)現(xiàn)沉積速率低的旋回地層，對(duì)應(yīng)測(cè)井響應(yīng)特征為泥巖的地層，沉積速率高的旋回地層，對(duì)應(yīng)測(cè)井響應(yīng)特征為砂巖的地層，即泥巖層的沉積速率小，砂巖層的沉積速率相對(duì)大.

4 結(jié)論

(1)基于頻譜分析與小波深頻分析的理論研究，通過對(duì)米氏周期進(jìn)行數(shù)值模擬并進(jìn)行頻譜分析與小波深頻分析，論證出兩種方法均可識(shí)別信號(hào)中的周期性信息，即可提取地層周期性、旋回性信息，特別是米氏旋回，用于沉積時(shí)間與沉積速率的研究；小波深頻分析還直觀體現(xiàn)出能量譜團(tuán)的變化、疊置旋回的個(gè)數(shù)及邊界，但由于沉積地層存在多變性，兩種分析方法結(jié)合起來使用效果更好，相互印證.

(2)建立了頻譜分析與小波深頻分析米氏旋回、沉積時(shí)間及速率的方法步驟與公式，并對(duì)鄂爾多斯盆地B159井長(zhǎng)7段自然伽馬數(shù)據(jù)分別做頻譜分析和小波深頻分析，識(shí)別出沉積厚度符合天文周期比的米氏旋回，表明長(zhǎng)7段沉積過程受到天文周期的影響，廣泛保存了米氏旋回.

(3)采用頻譜分析和小波深頻分析方法分別計(jì)算出鄂爾多斯盆地B159井長(zhǎng)7段的平均沉積速率為5.05 cm·ka-1和5.07 cm·ka-1，與前人研究結(jié)果吻合.

(4)鄂爾多斯盆地長(zhǎng)7段地層旋回主要受404 ka的長(zhǎng)偏心率周期控制，該地層共識(shí)別出6個(gè)長(zhǎng)偏心率旋回和20個(gè)短偏心率旋回，計(jì)算出長(zhǎng)7段124 m地層沉積時(shí)限約為2.4 Ma.

(5)形成了天文年代標(biāo)尺的計(jì)算方法步驟，并建立了鄂爾多斯盆地長(zhǎng)7段浮動(dòng)天文年代標(biāo)尺，有助于對(duì)地質(zhì)事件的持續(xù)時(shí)間做出精確的估計(jì).