基于變分模態(tài)分解的米蘭科維奇旋回劃分厚層砂礫巖地層層序的方法及其應用

丁 愷 趙福海 高蓮鳳 李丙喜 高晨陽 靳雪彬

（1. 遼寧工程技術大學礦業(yè)學院，遼寧 阜新 123000；2. 遼寧省礦產資源綠色開發(fā)地質保障重點實驗室，遼寧 阜新 123000；3. 中國石油大慶油田有限責任公司勘探事業(yè)部，黑龍江 大慶 163453；4. 北京天元云開科技有限公司，北京 100085）

0 引 言

中國東北地區(qū)地質條件復雜，大部分構造油氣藏也已被鉆探發(fā)現［1-2］，在這些地區(qū)尋找有利的成藏區(qū)帶，特別是有利地層與巖性圈閉存在很大困難，為了解決這些困難，需要進行細分層序的研究［3-4］。目前，時頻方法如小波變換、經驗模態(tài)分解（Empirical Mode Decomposition，EMD）等方法被廣泛應用于層序劃分中［5-8］。這些方法可以幫助分離出不同頻率和時間尺度上的地質事件，從而識別出不同的地層單元。然而，使用傳統的時頻方法存在一些局限性，如方法上的模態(tài)混疊、易受干擾問題［9］，以及僅用時頻分析結果指導劃分層序受人為因素影響的問題，這些都會影響層序劃分的準確性。

近年來，變分模態(tài)分解（Variational Mode Decomposition，VMD）這一新的時頻分析技術被提出［10-11］，該方法具有良好的分辨率，并在信號處理、圖像處理和地震勘探等領域得到廣泛應用，取得了顯著的效果。針對前文提出的傳統方法的局限性，本文引入VMD 方法應用于地層層序劃分，處理模態(tài)混疊和信號干擾的問題。而在解決人為因素影響方面，一些學者從天文地質學方面探討了天體運行軌道周期對沉積盆地層序地層形成的控制作用，并建立了各天文周期與地層基準面旋回之間的相對應關系［12-14］。米蘭科維奇旋回作為一種天然存在于地層中的地質時間標尺，可以對VMD 分解結果進行約束，從而進一步提高層序劃分的準確性。

類似研究區(qū)營城組四段這種發(fā)育厚層砂礫巖為主的地層，其測井曲線的形態(tài)變化不清晰，根據傳統的層序研究方法難以識別各地層的界面。本文將VMD 方法應用于層序劃分問題，并結合米蘭科維奇旋回對其進行約束約束，旨在可以解決多期次形成的地層界限認識不清的問題。

1 基礎原理及VMD優(yōu)勢

1.1 基礎原理

變分模態(tài)分解方法是由Konstantin Dragomiretskiy 提出的一種模態(tài)分解技術，但其原理與EMD 類分解算法完全不同，此方法為全新的信號非遞歸分解算法［15］。VMD 方法可以根據不同尺度提取出原始信號中實際蘊含的波動，為測井數據的非線性本征分解和地質分析提供有效的信號處理手段。VMD 方法比EMD 方法分解精度更高并且具有很強的魯棒性。

常規(guī)的EMD 方法在對非線性信號分解時采用的是遞歸方式，而這種方式依靠對極值點的尋求，缺少嚴謹的數學推導，EMD 方法存在嚴重的模態(tài)混疊問題［16］。VMD 方法不同于EMD 方法，VMD由若干個Wiener 濾波器的組成，是求最優(yōu)解的變分問題［17］。此外VMD 不同于EMD 類分解算法，對信號分解后得出的IMF（本征模態(tài)函數）分量層數不可提前知曉，VMD 在分解信號前需要預先設定分解層數k，因此VMD 方法也被視為一種階層可調的分解算法，同時也因為分解前給定了分解層數，給定的分解層數恰當也會很好地避免模態(tài)混疊現象的產生。相比于EMD 方法，VMD 方法具有更嚴謹的數學基礎，對信號分解采用的是一種完全非遞歸的方式，模態(tài)混疊的問題也得到了有效的抑制。將本征模態(tài)函數（IMF）定義為一個信號，其表達式為［15］

式中：k——分解層數；t——時間，s；uk(t)——各模態(tài)分量；φk(t)——各模態(tài)分量的瞬時相位，rad；Ak(t)——各模態(tài)分量的瞬時振幅，dB。

變分模態(tài)分解（VMD）方法主要對信號進行變換計算處理：

VMD 是一個約束變分問題，公式為［15］：

式中：j——虛數單位；s——常數；e——自然常數；uk——第k個模態(tài)分量；ωk——第k個模態(tài)分量的中心頻率；δ(t)——Dirichlet 函數；f(t)——輸入信號。

每個模式的帶寬由解析信號的L2的范數決定。引入二次罰項α 和拉格朗日乘子λ，將式（2）的約束問題轉化為非約束問題，即

式中：α——懲罰因子；λ——拉格朗日乘子；λ(t)——Lagrange 函數；δ(t)——對時間t的偏導數。

然后，通過交替方向的乘法算子求出式（4）的鞍點，即尋找式（2）最優(yōu)解并不斷更新迭代，，，具體步驟如下［15］：

式中：——輸入信號頻域的傅里葉等距變換；(ω)——離散模態(tài)頻域的傅里葉等距變換；(ω)——拉格朗日乘子頻域的傅里葉等距變換。

式中τ——拉格朗日乘子λ更新的步長。

（5）重復步驟（2）—步驟（4），對于給定判斷精度e = 2.7 ＞0，若則停止迭代。得到具有調幅-調號特征的k個IMF 分量信號。

VMD 方法的分解過程是假定所有分量都集中在各自中心頻率，通過尋找最優(yōu)解，圍繞中心頻率擾動，不斷更迭每個模態(tài)函數和中心頻率，從而使得最終分解的每個模態(tài)函數為具有窄帶特性的調幅-調頻信號。

1.2 模型測試

測井曲線本質上是不同周期地質活動在時間域或深度域上產生的地層旋回信息的疊合，本次研究設計了1 個由多個不同頻率的子信號合成的復合信號來模擬實際的地層情況。設計1 個由3 個頻率的信 號 合 成 的 模 擬 信 號X= 4 cos (40πt) +3 sin (60πt) + 5 sin (100πt)（圖1（a）），3 個子信號分別為20 Hz（圖1（b））、30 Hz（圖1（c））和50 Hz（圖1（d））。而在實際的工作中，會不可避免的受到噪聲的干擾，因此在純凈信號中分別加入信噪比為1 dB 和3 dB 的高斯白噪構成2 種含噪信號（圖1（e）、（f）），對信號進行分解。VMD分解結果表明，在已知信號成分的情況下，選取合適的k值，可分解近似得出3 個子頻率信號。利用VMD 分別對2 個含噪信號分解后發(fā)現（圖2（a）、（b）），含有1 dB 噪聲的信號分解出來的IMF 分量波形與含有3 dB 噪聲的信號差別不大。而EMD 分解因為不能預設分解的層數，只可以按照高頻、次高頻的順序依次分解，因此EMD 分解IMF 的個數不可控制。用EMD 方法對2 個含噪信號進行分解，發(fā)現分解出來的IMF 分量受模態(tài)混疊影響，波形變化嚴重（圖2（c）、（d））。

圖1 模擬信號及含噪信號Fig. 1 Analog signal and noisy signal

圖2 模擬信號加噪后分解結果Fig. 2 Decomposition results of noisy signal after adding noise

在本次測試中，采用了1 種基于VMD 方法的時頻分析技術，成功地識別了合成信號中的3 個子信號，并不受噪聲影響將3 個子信號分解出來。VMD 方法相較于其他時頻分析方法具有以下幾個優(yōu)點：

第一，VMD 方法具有高效的信號分解能力。通過將信號分解為多個獨立的分量，我們能夠更加精細地分析信號的頻譜特征，從而實現對米蘭科維奇旋回等地層特征的高精度分析。

第二，VMD 方法具有分解層數選擇的能力。在信號分解的過程中，VMD 方法能夠根據試驗結果調整分解參數，從而更好地適應信號的特征，提高分解的穩(wěn)定性和精度，避免過度分解和欠度分解。

第三，VMD 方法具有魯棒性強的特點。由于能夠對噪聲和干擾具有較強的魯棒性，VMD 方法能夠實現更為穩(wěn)定的分解結果，適用于對實際工程信號的處理和分析。

此外，VMD 方法計算效率高，適合處理大數據量的信號分析和處理。因此VMD 方法在米蘭科維奇旋回的識別和層序劃分等方面具有較大的應用價值。

2 研究方法

本次研究是基于變分模態(tài)分解的米氏旋回研究，因此在進行研究之前需要先通過對測井的頻譜分析來確定地層中包含米氏周期。確定米氏周期后，對測井曲線數據進行變分模態(tài)分解處理，得到多個特征分量，分析不同特征分量的頻譜，選取與米蘭科維奇旋回周期相對應的特征分量。根據選取的特征分量的曲線特征，結合巖性等地質信息，對地層進行層序劃分［18-19］。

2.1 米氏周期的確定

根據米蘭科維奇旋回理論，是地球軌道的周期性的改變引起太陽日射的變化，進而影響了海平面的周期性變化［20］，揭示了沉積記錄具有周期性規(guī)律的根源，當前確定的地球周期軌道參數包括長偏心率、短偏心率、斜率和歲差。長偏心率主要周期為405 ka，短偏心率主要周期為100 ka，斜率主要周期為29～45 ka，歲差主要周期為17～22 ka。米蘭科維奇旋回在漫長的地質歷史時期中具有相對的穩(wěn)定性，各周期的比率關系也相對穩(wěn)定，因此應用頻譜分析技術可以驗證沉積記錄中存在米氏旋回［21-22］。

測井曲線本質上是不同周期地質活動在時間域或深度域上產生的地層旋回信息的疊合，該疊合信息通過頻譜分析從深度域或時間域變換為頻率域，產生頻譜曲線。研究表明，沉積巖的自然放射性的強弱程度可以很好地對應大多數陸地環(huán)境下（除了快速堆積的粗相帶外）的沉積物泥質含量、粒度以及沉積速率的變化，因此自然伽馬數據可以作為分析氣候變化的依據。對自然伽馬曲線進行頻譜分析，頻譜能量強度表示地層在沉積時各天文軌道周期的占比，強度最大時為影響地層沉積的主導周期。

以研究區(qū)井1 為例，使用時頻分析軟件，以0.05 的顯著性水平為參考依據進行分析，得到自然伽馬曲線的頻譜分析圖（圖3（a））。根據頻譜能量強度假設主要頻率點A、B 和C，3 個點的頻率分別為0.055 56、0.013 542 和0.031 597 Hz（圖3（b）—（d））。波長為頻率的倒數，轉化為對應的波長為17.998 56、7.384 434 和3.164 857 m，波長本質上是地層中的沉積厚度，其厚度比值為1∶0.410 279∶0.175 839。本次選取的米蘭科維奇軌道參數為：周期100 ka 的短偏心率、周期41 ka 的短軸斜率和周期17 ka 的歲差，各周期的比值為1∶0.41∶0.17。研究發(fā)現米蘭科維奇各周期比值與厚度比值有著良好的對應關系［18］，證明了營四段中比較完整地保存了米蘭科維奇旋回。

圖3 井1營四段自然伽馬曲線和IMF頻譜圖Fig. 3 GR curves and IMF spectrum of Ying-4 Member in Well 1

2.2 層序地層的劃分

對自然伽馬曲線進行變分模態(tài)分解，得到多個特征分量IMF（圖4），其地質意義：分解得到的多個IMF 曲線的振蕩位置既是突變點，又是一種層序界面的響應，可與層序界面建立一定的關系，能夠作為層序分析的依據。對分解出來的各個IMF 進行時頻分析，發(fā)現IMF 曲線與自然伽馬曲線中包含的各周期天文軌道旋回信息有很好的對應關系，各IMF 曲線中的每個波峰都代表了該曲線所代表周期的一個旋回，應選取合適級別的IMF 劃分層序界面。

圖4 井1營四段IMF1—10Fig. 4 IMF1-10 of Ying-4 Member in Well 1

本文將井1 的自然伽馬曲線（圖4）通過VMD方法進行模態(tài)分解，進而得到多個IMF，通過對多個IMF 的頻譜分析，發(fā)現原始的自然伽馬曲線的頻譜分析結果（圖3（a））與IMF1（圖3（b））、IMF2（圖3（c））、IMF3（圖3（d））有著很好的對應關系，可以看作是對原始自然伽馬曲線的天文周期進行濾波。從圖3 可以看出，IMF1 對應100 ka 短偏心率，IMF2 對應41 ka 短斜率，IMF3對應17 歲差。

3 應用實例

3.1 研究區(qū)概況

徐家圍子斷陷位于松遼盆地北部，近南北向展布，整體走向NNW 向，早白堊世晚期斷陷開始向坳陷轉化［23］，營四段為陸相斷陷湖盆沉積。研究區(qū)位于徐家圍子斷陷中部的徐西凹陷，西臨中央古隆起，東接徐東坳陷。研究區(qū)目的層為營城組四段，營四段沉積時期斷裂對古地形的影響作用已非常弱，發(fā)育一套以砂礫巖為主的地層，是斷陷末期較為平緩的古地貌背景之上發(fā)育的一套厚度巨大的粗碎屑巖沉積建造，營四段沉積環(huán)境為湖相、沖積扇和辮狀河三角洲相［24-25］，后兩者屬于典型的高能沉積環(huán)境。營四段的沉積背景受到了區(qū)域構造背景和全球氣候變化的影響，其沉積特點呈現出明顯的周期性變化。

3.2 實例分析

前人將研究區(qū)營城組四段按巖性組合及曲線特征劃分了SQ1 和SQ2 兩個三級層序，SQ1 巖性以泥巖為主，儲層物性差，SQ2 發(fā)育大套砂礫巖，氣藏開發(fā)潛力巨大［26］。但是營四段地層厚度大，上部層序（SQ2）存在的大套砂礫巖地層層序標志難以識別，導致出現四級層序認識不清，不利于認識儲集砂體的分布。本文根據米蘭科維奇旋回理論，識別出地層不同級次的旋回，并以米蘭科維奇旋回作為高頻層序劃分的標尺，為層序劃分提供了不受人為影響的標準，保證層序劃分方案的統一性，可以大幅提高其準確度。

前文對研究區(qū)井1 的自然伽馬曲線進行頻譜分析后發(fā)現了100、41 和13 ka 的周期。根據層序地層學理論，四級層序受405 ka 周期的長偏心率控制，五級層序受100 ka 周期的短偏心率控制。然而，頻譜分析并未發(fā)現對應四級層序的405 ka 周期的長偏心率，因此基于100 ka 周期的短偏心率將營四段細分為五級層序。對自然伽馬曲線進行VMD 處理，得到10 個IMF（圖4）。

通過頻譜分析對比發(fā)現（圖3），IMF1 分量對應100 ka 周期的短偏心率，IMF1 中每個波峰都代表著一個100 ka 周期?；贗MF1 在營四段識別出21 個五級層序。根據長短偏心率的近似比值4∶1，將4 個五級層序合成1 個四級層序，并結合巖性將營四段上部細分為SQ2-4、SQ2-3、SQ2-2 和SQ2-1共4 個四級層序（圖5），其中SQ2-4 為營四段頂部層序，井1 在該時期處于構造高部位，判斷頂部沉積被剝蝕，該四級層序由2 個五級層序組成，巖石組合為含礫砂巖夾灰白色泥巖，IMF1 曲線呈現周期性下降趨勢，海平面下降導致沉積物侵蝕，為退積型沉積類型，即水退旋回。SQ2-3、SQ2-2 以及SQ2-1 均由4 個五級層序組成，巖性組合為厚層砂礫巖夾泥巖。SQ2-1 和SQ2-3 的 IMF1 曲線呈現周期性上升趨勢，海平面上升導致沉積物堆積，判斷為進積型沉積類型，即水進旋回，SQ2-2 對應的IMF1 曲線為周期性下降趨勢，該層為水退旋回。營四段下部可劃分為SQ1-1 和SQ1-2 共2 個四級層序，其中SQ1-1 為營四段最底部層序，電阻率曲線突然劇烈變化，且營一段的灰白色凝灰?guī)r與營四段巖性明顯不同，判斷底部發(fā)生沉積間斷，該層序由3 個五級層序組成，巖性組合為淺灰色深灰色泥巖，SQ1-2 由4 個五級層序組成，巖性組合為砂礫巖夾泥巖，IMF1 曲線為周期性下降趨勢，兩個都為水退旋回，最終完成四級層序的劃分。

圖5 井1營四段層序地層劃分Fig. 5 Sequence stratigraphic division of Ying-4 Member in Well 1

整個營城組四段共劃分為6 個四級層序，如圖5 所示（黑色箭頭為水退旋回，紅色箭頭為水進旋回）。其中SQ2 砂礫巖地層劃分出3 個層序界面，由圖5 中巖心照片可以看出，3 673 m 的巖性由大套灰色砂礫巖轉為雜色砂礫巖夾泥巖，3 743 m 的巖性由粗礫砂礫巖轉為細礫砂礫巖，3 814 m 的巖性由細礫砂礫巖轉為粗礫砂礫巖，3 878 m 的巖性由泥巖轉砂礫巖，所劃分的層序與巖性變化一致，符合區(qū)域地質認識。

4 結 論

（1）VMD 方法相比與其他時頻方法具有更好的魯棒性和更高的精確性，相比于EMD 方法克服了模態(tài)混疊的問題。對自然伽馬曲線進行VMD 處理，可獲得不同周期的米氏旋回。

（2）基于不同周期米氏旋回對高頻沉積旋回的控制作用，對松遼盆地徐家圍子斷陷井1 進行細分層序劃分，所劃分的細分層序與巖性變化一致，同時可將厚層砂礫巖的沉積旋回特征展現的更清晰，符合區(qū)域地質認識。