基于变分模态分解的米兰科维奇旋回划分厚层砂砾岩地层层序的方法及其应用

丁 恺 赵福海 高莲凤 李丙喜 高晨阳 靳雪彬

（1. 辽宁工程技术大学矿业学院，辽宁 阜新 123000；2. 辽宁省矿产资源绿色开发地质保障重点实验室，辽宁 阜新 123000；3. 中国石油大庆油田有限责任公司勘探事业部，黑龙江 大庆 163453；4. 北京天元云开科技有限公司，北京 100085）

0 引 言

中国东北地区地质条件复杂，大部分构造油气藏也已被钻探发现［1-2］，在这些地区寻找有利的成藏区带，特别是有利地层与岩性圈闭存在很大困难，为了解决这些困难，需要进行细分层序的研究［3-4］。目前，时频方法如小波变换、经验模态分解（Empirical Mode Decomposition，EMD）等方法被广泛应用于层序划分中［5-8］。这些方法可以帮助分离出不同频率和时间尺度上的地质事件，从而识别出不同的地层单元。然而，使用传统的时频方法存在一些局限性，如方法上的模态混叠、易受干扰问题［9］，以及仅用时频分析结果指导划分层序受人为因素影响的问题，这些都会影响层序划分的准确性。

近年来，变分模态分解（Variational Mode Decomposition，VMD）这一新的时频分析技术被提出［10-11］，该方法具有良好的分辨率，并在信号处理、图像处理和地震勘探等领域得到广泛应用，取得了显著的效果。针对前文提出的传统方法的局限性，本文引入VMD 方法应用于地层层序划分，处理模态混叠和信号干扰的问题。而在解决人为因素影响方面，一些学者从天文地质学方面探讨了天体运行轨道周期对沉积盆地层序地层形成的控制作用，并建立了各天文周期与地层基准面旋回之间的相对应关系［12-14］。米兰科维奇旋回作为一种天然存在于地层中的地质时间标尺，可以对VMD 分解结果进行约束，从而进一步提高层序划分的准确性。

类似研究区营城组四段这种发育厚层砂砾岩为主的地层，其测井曲线的形态变化不清晰，根据传统的层序研究方法难以识别各地层的界面。本文将VMD 方法应用于层序划分问题，并结合米兰科维奇旋回对其进行约束约束，旨在可以解决多期次形成的地层界限认识不清的问题。

1 基础原理及VMD优势

1.1 基础原理

变分模态分解方法是由Konstantin Dragomiretskiy 提出的一种模态分解技术，但其原理与EMD 类分解算法完全不同，此方法为全新的信号非递归分解算法［15］。VMD 方法可以根据不同尺度提取出原始信号中实际蕴含的波动，为测井数据的非线性本征分解和地质分析提供有效的信号处理手段。VMD 方法比EMD 方法分解精度更高并且具有很强的鲁棒性。

常规的EMD 方法在对非线性信号分解时采用的是递归方式，而这种方式依靠对极值点的寻求，缺少严谨的数学推导，EMD 方法存在严重的模态混叠问题［16］。VMD 方法不同于EMD 方法，VMD由若干个Wiener 滤波器的组成，是求最优解的变分问题［17］。此外VMD 不同于EMD 类分解算法，对信号分解后得出的IMF（本征模态函数）分量层数不可提前知晓，VMD 在分解信号前需要预先设定分解层数k，因此VMD 方法也被视为一种阶层可调的分解算法，同时也因为分解前给定了分解层数，给定的分解层数恰当也会很好地避免模态混叠现象的产生。相比于EMD 方法，VMD 方法具有更严谨的数学基础，对信号分解采用的是一种完全非递归的方式，模态混叠的问题也得到了有效的抑制。将本征模态函数（IMF）定义为一个信号，其表达式为［15］

式中：k——分解层数；t——时间，s；uk(t)——各模态分量；φk(t)——各模态分量的瞬时相位，rad；Ak(t)——各模态分量的瞬时振幅，dB。

变分模态分解（VMD）方法主要对信号进行变换计算处理：

VMD 是一个约束变分问题，公式为［15］：

式中：j——虚数单位；s——常数；e——自然常数；uk——第k个模态分量；ωk——第k个模态分量的中心频率；δ(t)——Dirichlet 函数；f(t)——输入信号。

每个模式的带宽由解析信号的L2的范数决定。引入二次罚项α 和拉格朗日乘子λ，将式（2）的约束问题转化为非约束问题，即

式中：α——惩罚因子；λ——拉格朗日乘子；λ(t)——Lagrange 函数；δ(t)——对时间t的偏导数。

然后，通过交替方向的乘法算子求出式（4）的鞍点，即寻找式（2）最优解并不断更新迭代，，，具体步骤如下［15］：

式中：——输入信号频域的傅里叶等距变换；(ω)——离散模态频域的傅里叶等距变换；(ω)——拉格朗日乘子频域的傅里叶等距变换。

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

（5）重复步骤（2）—步骤（4），对于给定判断精度e = 2.7 ＞0，若则停止迭代。得到具有调幅-调号特征的k个IMF 分量信号。

VMD 方法的分解过程是假定所有分量都集中在各自中心频率，通过寻找最优解，围绕中心频率扰动，不断更迭每个模态函数和中心频率，从而使得最终分解的每个模态函数为具有窄带特性的调幅-调频信号。

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 个含噪信号分解后发现（图2（a）、（b）），含有1 dB 噪声的信号分解出来的IMF 分量波形与含有3 dB 噪声的信号差别不大。而EMD 分解因为不能预设分解的层数，只可以按照高频、次高频的顺序依次分解，因此EMD 分解IMF 的个数不可控制。用EMD 方法对2 个含噪信号进行分解，发现分解出来的IMF 分量受模态混叠影响，波形变化严重（图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 方法相较于其他时频分析方法具有以下几个优点：

第一，VMD 方法具有高效的信号分解能力。通过将信号分解为多个独立的分量，我们能够更加精细地分析信号的频谱特征，从而实现对米兰科维奇旋回等地层特征的高精度分析。

第二，VMD 方法具有分解层数选择的能力。在信号分解的过程中，VMD 方法能够根据试验结果调整分解参数，从而更好地适应信号的特征，提高分解的稳定性和精度，避免过度分解和欠度分解。

第三，VMD 方法具有鲁棒性强的特点。由于能够对噪声和干扰具有较强的鲁棒性，VMD 方法能够实现更为稳定的分解结果，适用于对实际工程信号的处理和分析。

此外，VMD 方法计算效率高，适合处理大数据量的信号分析和处理。因此VMD 方法在米兰科维奇旋回的识别和层序划分等方面具有较大的应用价值。

2 研究方法

本次研究是基于变分模态分解的米氏旋回研究，因此在进行研究之前需要先通过对测井的频谱分析来确定地层中包含米氏周期。确定米氏周期后，对测井曲线数据进行变分模态分解处理，得到多个特征分量，分析不同特征分量的频谱，选取与米兰科维奇旋回周期相对应的特征分量。根据选取的特征分量的曲线特征，结合岩性等地质信息，对地层进行层序划分［18-19］。

2.1 米氏周期的确定

根据米兰科维奇旋回理论，是地球轨道的周期性的改变引起太阳日射的变化，进而影响了海平面的周期性变化［20］，揭示了沉积记录具有周期性规律的根源，当前确定的地球周期轨道参数包括长偏心率、短偏心率、斜率和岁差。长偏心率主要周期为405 ka，短偏心率主要周期为100 ka，斜率主要周期为29～45 ka，岁差主要周期为17～22 ka。米兰科维奇旋回在漫长的地质历史时期中具有相对的稳定性，各周期的比率关系也相对稳定，因此应用频谱分析技术可以验证沉积记录中存在米氏旋回［21-22］。

测井曲线本质上是不同周期地质活动在时间域或深度域上产生的地层旋回信息的叠合，该叠合信息通过频谱分析从深度域或时间域变换为频率域，产生频谱曲线。研究表明，沉积岩的自然放射性的强弱程度可以很好地对应大多数陆地环境下（除了快速堆积的粗相带外）的沉积物泥质含量、粒度以及沉积速率的变化，因此自然伽马数据可以作为分析气候变化的依据。对自然伽马曲线进行频谱分析，频谱能量强度表示地层在沉积时各天文轨道周期的占比，强度最大时为影响地层沉积的主导周期。

以研究区井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。研究发现米兰科维奇各周期比值与厚度比值有着良好的对应关系［18］，证明了营四段中比较完整地保存了米兰科维奇旋回。

图3 井1营四段自然伽马曲线和IMF频谱图Fig. 3 GR curves and IMF spectrum of Ying-4 Member in Well 1

2.2 层序地层的划分

对自然伽马曲线进行变分模态分解，得到多个特征分量IMF（图4），其地质意义：分解得到的多个IMF 曲线的振荡位置既是突变点，又是一种层序界面的响应，可与层序界面建立一定的关系，能够作为层序分析的依据。对分解出来的各个IMF 进行时频分析，发现IMF 曲线与自然伽马曲线中包含的各周期天文轨道旋回信息有很好的对应关系，各IMF 曲线中的每个波峰都代表了该曲线所代表周期的一个旋回，应选取合适级别的IMF 划分层序界面。

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

本文将井1 的自然伽马曲线（图4）通过VMD方法进行模态分解，进而得到多个IMF，通过对多个IMF 的频谱分析，发现原始的自然伽马曲线的频谱分析结果（图3（a））与IMF1（图3（b））、IMF2（图3（c））、IMF3（图3（d））有着很好的对应关系，可以看作是对原始自然伽马曲线的天文周期进行滤波。从图3 可以看出，IMF1 对应100 ka 短偏心率，IMF2 对应41 ka 短斜率，IMF3对应17 岁差。

3 应用实例

3.1 研究区概况

徐家围子断陷位于松辽盆地北部，近南北向展布，整体走向NNW 向，早白垩世晚期断陷开始向坳陷转化［23］，营四段为陆相断陷湖盆沉积。研究区位于徐家围子断陷中部的徐西凹陷，西临中央古隆起，东接徐东坳陷。研究区目的层为营城组四段，营四段沉积时期断裂对古地形的影响作用已非常弱，发育一套以砂砾岩为主的地层，是断陷末期较为平缓的古地貌背景之上发育的一套厚度巨大的粗碎屑岩沉积建造，营四段沉积环境为湖相、冲积扇和辫状河三角洲相［24-25］，后两者属于典型的高能沉积环境。营四段的沉积背景受到了区域构造背景和全球气候变化的影响，其沉积特点呈现出明显的周期性变化。

3.2 实例分析

前人将研究区营城组四段按岩性组合及曲线特征划分了SQ1 和SQ2 两个三级层序，SQ1 岩性以泥岩为主，储层物性差，SQ2 发育大套砂砾岩，气藏开发潜力巨大［26］。但是营四段地层厚度大，上部层序（SQ2）存在的大套砂砾岩地层层序标志难以识别，导致出现四级层序认识不清，不利于认识储集砂体的分布。本文根据米兰科维奇旋回理论，识别出地层不同级次的旋回，并以米兰科维奇旋回作为高频层序划分的标尺，为层序划分提供了不受人为影响的标准，保证层序划分方案的统一性，可以大幅提高其准确度。

前文对研究区井1 的自然伽马曲线进行频谱分析后发现了100、41 和13 ka 的周期。根据层序地层学理论，四级层序受405 ka 周期的长偏心率控制，五级层序受100 ka 周期的短偏心率控制。然而，频谱分析并未发现对应四级层序的405 ka 周期的长偏心率，因此基于100 ka 周期的短偏心率将营四段细分为五级层序。对自然伽马曲线进行VMD 处理，得到10 个IMF（图4）。

通过频谱分析对比发现（图3），IMF1 分量对应100 ka 周期的短偏心率，IMF1 中每个波峰都代表着一个100 ka 周期。基于IMF1 在营四段识别出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 为营四段最底部层序，电阻率曲线突然剧烈变化，且营一段的灰白色凝灰岩与营四段岩性明显不同，判断底部发生沉积间断，该层序由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 的岩性由泥岩转砂砾岩，所划分的层序与岩性变化一致，符合区域地质认识。

4 结 论

（1）VMD 方法相比与其他时频方法具有更好的鲁棒性和更高的精确性，相比于EMD 方法克服了模态混叠的问题。对自然伽马曲线进行VMD 处理，可获得不同周期的米氏旋回。

（2）基于不同周期米氏旋回对高频沉积旋回的控制作用，对松辽盆地徐家围子断陷井1 进行细分层序划分，所划分的细分层序与岩性变化一致，同时可将厚层砂砾岩的沉积旋回特征展现的更清晰，符合区域地质认识。