经验模态分解法在测井层序界面识别中的应用

赵妮, 李瑞, 陈立军

(1.油气藏地质开发工程国家重点实验室(成都理工大学), 四川 成都 610059; 2.西安工业大学计算机科学与工程学院, 陕西 西安 710032; 3.延长油田股份有限公司勘探开发技术研究中心, 陕西 延安 716000)

0 引 言

层序地层学是研究旋回式、成因上有联系的、以侵蚀面或者与其可以对比的整合面为界的年代地层格架以及沉积层序内部地层、岩相分布模式的地层学分支学科。一般采用露头、岩心、测井和地震资料进行层序地层划分和层序地层界面识别,在露头和岩心资料不足、地震资料分辨率有限的情况下,测井资料在层序界面划分方面就显得尤为重要。随着高分辨率层序地层学理论的提出,小波分析、Hilbert-Huang变换等新方法也应用于层序地层划分。Hilbert-Huang变换是基于Hilbert谱的算法,包括经验模态分解(Empirical Mode Decomposition,EMD)和希尔伯特谱分析2部分。经验模态分解即EMD方法针对非平稳信号提取局部均值曲线,将复杂信号分解为一系列分量函数,以得到具有明显物理意义的瞬时频率,可以很好地处理非平稳、非线性信号[1-2]。本文将EMD方法应用于测井数据处理当中,为测井层序界面识别提供了一种有效的信号处理手段。

1 经验模态分解(EMD)方法

1.1 EMD方法理论基础

EMD分解是基于数据本身的一种自适应变换方法,其本质在于利用信号的特征时间尺度获得信号的固有波动模式。与小波变换和傅里叶变换相比,它可以很好地处理非平稳、非线性信号。该方法建立在以下假设之上:①信号至少有1个极大值和1个极小值点;②特征时间尺度由极值点之间的时间间隔定义;③对于没有极值点但有变形点的情况,可以通过微分获取极值点再通过积分获得分解结果。

1.2 EMD分解流程

EMD分解流程见图1。其中,固有模态函数应满足2个条件(即流程中的IMF条件)[3]:①极值点数目和过0点数目相等或相差1个;②在任意点处,最大值包络线与最小值包络线的均值为0。第1个条件保证了数据的局部最大值为正,极小值为负;第2个条件剔除了瞬时频率中由于波形不对称引起的震荡。

图1 EMD分解流程

从图1可见,EMD分解是将原信号x(t)分解成一系列从高频到低频的固有模态函数IMFi(t)和一个趋势残余项R(t),即

(1)

1.3 EMD方法的优点

傅里叶变换和小波分析是测井信号分析中常用的方法。傅里叶变换是经典信号频谱分析的基础,将原信号分解为若干正弦信号的加权和,其中每1个正弦信号对应1个固定的频率和固定的幅值,适合于分析频率不随时间变化的平稳信号。测井数据是非平稳的空间信号,经傅里叶变换后得到的频率信号在频域上有最大分辨率,但其本身不包含深度定位信息,也就是说可以找出最大频率,但不能确定其深度位置。小波变换是傅里叶变换的发展,可以根据频率高低自动调节窗口大小,具有多分辨分析功能。离散小波分解是通过2个互补的滤波器产生低频的近似信号和高频的细节信号,在应用中通常认为低频部分是最重要的。小波变换的局限性在于,①小波基函数的选取会影响分解信号的精度;②每次提取的分量不是完整的信息,确切的说是信号经过滤波器后的部分信息。而EMD分解是对整个信号进行的分解,每一个分量中都包含了信号完整的信息,反映信号的瞬时频率特征,这也是用EMD分解划分旋回优于小波变换的主要原因。

2 EMD分解在测井层序识别中的应用

2.1 指标的选择与处理

在测井数据采集过程中,各种干扰因素不可避免地影响信号的分辨率。对相关性高的敏感指标进行线性组合是提高测井解释准确性的有效方法。

(1) 针对某一区块选择取心井测井数据,通过多指标的雷达图比较法,并与岩心进行对比,选取能反映岩性变化特征的敏感指标,再对选出的敏感指标求取相关系数,选择相关性较好的指标进行组合。比如,自然电位和自然伽马的组合、声波时差和密度的组合等。对多个线性相关的指标进行组合是对地层特征信息的强化。由于燥声是随机的,叠加信号突出了有用信息,降低了干扰的影响。

(2) 不同测井指标的量纲往往不相同。为了将不同量纲的指标组合在一起,需要进行指标的归一化处理,再将无量纲的数据进行组合,构成新的分析指标。

2.2 地层基准面旋回识别方法

2.2.1 数据处理原则

EMD分解可以将原始信号分解成不同频率的内蕴函数IMF(i)和一个趋势残余项R(t),每个IMF(i)为信号在某一特征尺度参数上的模态,其中相对频率较高的分量反映了主要受天文因素控制的短期旋回特征,低频分量则反映了受构造作用控制的较为稳定的中长期旋回特征;R(t)为信号本身的趋势或仪器漂零[4]。利用这一规律,对EMD分解后得到的固有模态函数IMF(i)展开相关性研究。①计算各阶IMF(i)与原信号的相关系数。往往那些与原信号相关性好的IMF(i)也是揭示稳定沉积环境的最佳信号分量。②再分析各阶IMF(i)之间的相关性,将相关性好的IMF(i)进行叠加组合。相关性好说明有一致性的变化规律,将其叠加可以加强反映地质特征的有用信息,削弱干扰因素,从而提高识别地层界面的精确度。

2.2.2 EMD分解响应特征

测井信号经过EMD分解后,可以从不同尺度上提取出原测井信号中的时频特征,反映出不同的沉积旋回级次。若能与各级层序界面建立一定的对应关系,则可以作为测井层序地层划分的依据。

当河流在动力学的作用下达到平衡时,将此时的平衡剖面抽象为一个势能面,称之为基准面[5-6]。地表或沉积界面不断地通过沉积或侵蚀作用向靠近基准面的方向运动以便达到新的平衡,这样的一个旋回升降过程称为成因层序。层序地层的这种旋回性在测井资料上表现为粒度和岩性的韵律性,反映到测井曲线上表现为不同的幅度和频率特征,而曲线突变点往往是层序界面[7-9]。据此对测井指标的EMD分解曲线响应特征在各种沉积环境中的形态表现作镜像图,概括出5种基本形态类型。不同沉积环境的测井曲线形态大都是这5种基本形态的组合或演化(见表1)。

表1 测井指标EMD分解曲线响应特征的沉积环境基本类型

3 应用实例

应用EMD分解识别地层基准面旋回的方法对四川盆地川西坳陷中段须家河组×井层序地层进行了划分。

3.1 测井指标组合分析

自然伽马反映沉积物的粒度、分选及泥质含量,可以用来判断岩性、进行地层对比及估算泥质含量[10]。通常情况下,自然伽马在砂岩中读数为低值,在泥岩中读数为高值,但若砂岩中存在诸如云母、锆石等放射性矿物,会直接影响自然伽马的读值。自然电位常用于砂泥岩剖面中划分岩性、进行地层对比、确定渗透层界面位置、估计地层中的泥质含量、展开沉积相研究等方面[11]。通过分析,这2个测井参数均是适合该地区划分层序的灵敏指标,并且曲线指示具有很好的对应性,可以进行组合。

(1) 分析自然伽马和自然电位指标的相关性,得到它们之间的相关系数值是0.438 3,说明它们的相关性较好。

(2) 分别对自然伽马和自然电位进行归一化处理,选用线性归一化公式,即

y=(x-xmin)/(xmax―xmin)

(2)

式中,x为指标测值;xmax、xmin分别为x的最大值和最小值;y为对x归一化后的值。

(3) 对归一化后的自然伽马和自然电位进行叠加,取均值,组合成新的指标记为GR-SP。

(4) 分析GR-SP指标与自然伽马和自然电位的相关性,得到GR-SP与自然伽马的相关系数为0.936 3;GR-SP与自然电位的相关系数为0.725 9,说明新指标与自然伽马和自然电位都具有较高的相关性,可以替代原始指标,具有应用价值。

3.2 旋回识别与地层划分

对四川盆地川西坳陷中段须家河组×井的GR-SP数据进行EMD分解得到本征模态函数IMF(i)(见图2)。

图2 GR-SP指标EMD分解图

从图2可以看到,IMF(1)至IMF(4),显示出强烈的震荡,但高幅值出现的位置又不尽相同,这说明高频分量中存在噪音信号的影响。随着各阶IMF信号频率的降低,信号中稳定的变化趋势越来越强,尤其是IMF(5)、IMF(6)、IMF(7)、IMF(9)明显具有一致的变化规律,这种一致性反映出地层相对稳定的沉积环境特征。

旋回级次是对地层特征在不同尺度上的划分。EMD各阶分量中的高频分量显示出包含细节的多少,而这些细节分量中包含有用的信息也包含了噪音的干扰。沉积旋回是一个相对稳定的变化过程,更应当注重信号中的中低频成分,那些高频信号的加入只是起到添加剂的作用,这一点与小波分析的观点是一致的。

设Ri为GR-SP与其EMD分解函数IMF(i)的相关系数,利用MATLAB提供的相关函数corrcoef()可求出原信号与固有模态函数IMF(i)的相关系数值(见表2)。

表2 GR-SP与其EMD分解函数IMF(i)的相关系数值

从表2可知,EMD分解函数IMF(9)、IMF(5)、IMF(6)、IMF(7)与GR-SP均有较好的相关性。这一结果同前面EMD分解图得出的结论是一致的。同理,再计算出IMF(i)函数之间的相关系数值(见表3)。

从表3中可以得出相关系数大于0.1的有IMF(6)与IMF(7),为0.118 184;IMF(5)与IMF(9)的相关系数为0.106 921,说明IMF(6)与IMF(7)、IMF(5)与IMF(9)具有较好的相关性。相关系数接近0.1的还有IMF(1)与IMF(2)、IMF(3)与IMF(4)、IMF(5)与IMF(6),可以作为参考组合。将相关性好的EMD分解信号进行组合,目的是加强反映地层特征的固有信息,降低干扰的影响,比较结果说明在所有组合信号中只有那些同时与原信号相关性高的组合才能更好的揭示出地层的旋回特征。

表3 EMD分解本征模态函数IMF(i)之间的相关系数值

已知×井信号取样间隔为0.1 m,IMF(6)+IMF(7)的瞬时频率范围为0.003 4～0.516 1,进一步算出对应可分辨的深度为0.2～29.4 m;同理,可求出IMF(5)+IMF(9)的瞬时频率范围为(9.931 1×10-4)～0.503 4,其对应可分辨的深度为0.2～100.69 m。不同级次的旋回通常以分布时限来划分,体现到层次实体规模上,一般短期旋回厚度在数米级至数十米级,中期旋回厚度在数十米至近百米,长期旋回厚度在近百米至数百米。据此,可将IMF(6)+IMF(7)作为短旋回的参考依据;IMF(5)+IMF(9)作为中长期的参考依据。具体旋回划分结果见图3。

图3 EMD分解法识别的×井短期和中长期旋回层序地层柱状图

×井位于四川盆地川西坳陷中段孝泉-丰谷构造带。图3分析的是须家河组二段地层的基准面旋回特征,岩性以细、中粒岩屑砂岩、页岩为主,属于三角洲前缘亚相。依据表1指示的曲线形态特征对EMD分解后的不同频率组合信号进行分析,可以明显识别短期和中长期旋回特征。现对沉积微相分析如下。①水下分支河道:岩性为灰色细-中粗粒岩屑砂岩,表现为正粒序旋回特征;发育交错层理、斜层理、块状层理,可见冲刷-充填构造。②河口坝:岩性为分选好的细砂或粉砂组成,发育交错层理,具有反粒序旋回特征。③支流间湾:岩性以黏土为主,含少量的粉砂及细砂,具有水平层理。④席状砂:岩性为分选好的细砂或粉砂组成,发育交错层理,具有反粒序旋回特征。⑤远砂坝:位于河口坝前缘较远部位,岩性主要为粉砂,砂体内部发育交错层理,具有反粒序旋回特征。

4 结 论

(1) 测井指标从不同角度反映了地层的岩性、构造、层序及储层的物性特征,针对地区优选灵敏指标加以组合可以提高测井解释的准确性。

(2) 自然电位和自然伽马是测井资料中能反映地质构造和岩性的灵敏指标,采用数学方法对其进行组合,新指标是对地层固有特征的强化,可以提高分层效果。

(3) EMD分解是一种全新的时频分析方法,不仅保证了信号分解后的非平稳特征,还避免了加窗傅里叶变换受时窗限制以及小波变换受基函数影响的不足,可以更准确地描述信号时变特征。

(4) 对多个测井组合指标进行EMD分解,探讨各阶本征模态函数及其与原指标的相关性,寻找有价值的组合规律,创建IMF(i)组合信号,有利于更好地提取异常,改善识别地层基准面旋回的效果。

参考文献:

[1] 李彦军, 胡祥云. 何展翔HHT变换在地球物理中的应用现状及前景 [J]. 工程地球物理学报, 2010, 7(5): 537-543

[2] Huang N E, Zheng Shen, Long S R, et al. The Empirical Mode Decomposition and the Hilbert Spectrum for Nonlinear Non-stationary Time Series Analysis [C]∥Proc R Soc London. 1998, 45(4): 903-995.

[3] 郑天翔, 杨力华. 经验模式分解算法的探讨和改进 [J]. 中山大学学报: 自然科学版, 2007, 46(1): 1-6.

[4] 张婷, 聂春燕. 基于EMD时频分析的声波测井油层的信息提取 [J]. 长春大学学报, 2011, 21(8): 15-18.

[5] 郑荣才, 尹世民, 彭军. 基准面旋回结构与叠加样式的沉积动力学分 [J]. 沉积学报, 2001, 19(2): 259-255.

[6] 郑荣才, 彭军, 吴朝容. 陆相盆地基准面旋回的级次划分和研究意义 [J]. 沉积学报, 2000, 18(3): 249-255.

[7] 郭小燕, 王成, 于雷. 测井资料识别层序的方法及问题讨论 [J]. 西部探矿工程, 2011, 2: 67-70.

[8] 周祺, 郑荣才, 李凤杰, 等. 测井曲线在陆相层序地层界面识别中的应用 [J]. 大庆石油地质与开发, 2008, 27(4): 135-138.

[9] 李新虎. 测井曲线拐点在测井层序地层分析中的应用研究 [J]. 天然气地球科学, 2006, 17(6): 815-819.

[10] 曾婷, 桂志先, 何加成, 等. 储层预测中有关测井参数的分析及应用 [J]. 工程地球物理学报, 2010, 7(3): 296-299.

[11] 李君文, 陈洪德, 田景春. 鄂尔多斯盆地石千峰组高分辨率层序地层学特征 [J]. 大庆石油地质与开发, 2007, 26(1): 45-49.