基于CEEMDAN的希尔伯特变换在海底天然气水合物地震探测中的应用

夏秋萍，刘怀山,2

(1.中国海洋大学 海底科学与探测技术教育部重点实验室，山东 青岛 266100；2.海洋国家实验室 海洋矿产资源评价与探测技术功能实验室，山东 青岛 266071)

1 引 言

天然气水合物是天然气与水在高压低温条件下形成的类似冰状的结晶物质，目前已被正式确定为新的矿产资源，又因为其绿色清洁没有污染，并且在世界上均有分布，因此多个国家已经开始进行天然气水合物的勘探开发。研究人员应用各种高分辨率地震调查技术进行勘探研究来寻找水合物，目前已经发现了很多标志能够帮助证明水合物的存在，比如BSR(Bottom Simulating Reflector)、极性反转等[1-5]。地震资料的瞬时属性有三个方面，分别是瞬时振幅、相位、频率，以上属性对于处理复杂地层，能更加全面地反映储层的特征，在天然气储层预测中得到了广泛应用。

1998年，Huang等[6]提出了一种新的方法，用来进行时频分析，这个方法就是希尔伯特-黄变换(HHT)。该方法通过对信号进行经验模态分解EMD(Empirical Mode Decomposition)，分解出固有模态函数IMF(Intrinsic Mode Function)分量之后，对照实际信息选取分量进行Hilbert变换来获得瞬时参数[7-11]。这个方法改善了时频分析的局限性，地震资料的分辨率也得到了提高。但是信号中若存在间断或者干扰会引起模态混叠从而导致分解的准确程度有所降低，因此在地震数据中应用EMD进行处理与分析的方法也受到了一定制约[12]。为了解决EMD存在的问题，Wu等[13]提出了集合经验模态分解EEMD(Ensemble Empirical Mode Decomposition)方法，然而该方法无法保证重构信号的准确性，同时会产生很多非必要的迭代过程，大大降低了计算效率[14，15]。

2011年，Torres等[16]提出了自适应噪声的完备总体集合经验模态分解(Complete Ensemble Empirical Mode Decomposition of Adaptive Noise，CEEMDAN)，该方法是建立在经验模态分解的基础之上的，对信号进行经验模态分解，从高频到低频获得一系列IMF分量的过程中，将有限次的自适应白噪声加进每个分解阶段来消除分解过程产生的模态混叠现象，能在较少进行平均的条件下,对信号进行重构且几乎没有误差，能够得到更好的模态频谱分离结果。

本文将CEEMDAN应用于天然气水合物的OBS数据与垂直缆数据中，将经过CEEMDAN分解得到的IMF分量进行Hilbert变换，获得瞬时振幅、相位、频率，将其所形成的瞬时属性图与只进行了Hilbert变换获得的瞬时属性图进行对比，分析得到经过CEEMDAN处理后的瞬时属性图分辨率更高，可以当作地震勘探数据精细解释的工具，为进行储层分析提供有效支撑。

2 方法原理

2.1 EMD算法

Huang将具有以下两个性质的信号定义为IMF，第一是该信号有同等数量的极值点和零交叉点，如有差异，最多相差1个；第二是由这个信号绘制的极值上下包络线均值为零。EMD算法如下[6]：

给定实信号x(t)，首先找出x(t)的所有极值点，利用插值函数拟合所有的局部极大和极小值点，确定上下包络线，求上下包络线的平均值为h1(t),

m1(t)=x(t)-h1(t)

(1)

把m1(t)作为新信号，重复式(1)过程，得到上下包络的均值h2(t)

m2(t)=m1(t)-h2(t)

(2)

重复k次，直到mk(t)满足IMF条件，得到给定信号的第一阶IMF分量

mk(t)=mk-1(t)-hk(t)

(3)

第一阶IMF分量是频率最高的一阶分量，记为c1(t)，从给定信号x(t)中减去c1(t),获得残余量r1(t),

r1(t)=x(t)-c1(t)

(4)

把r1(t)作为新的给定信号，重复上述过程以获得全部的IMF分量，当rn(t)比预定值小或者rn(t)变为单调函数的时候，分解停止。最终得到给定信号:

(5)

2.2 CEEMDAN算法

从本质上来说，自适应噪声的总体集合经验模态分解[17]依然是以EMD为基础进行改进得到的方法[18]。类似于EEMD，对信号进行噪声辅助分析的方法，与EMD和EEMD不同的是，CEEMDAN克服了模态混叠，在提高计算效率的同时精确重构了原始信号，得到了不错的模态分离谱[19]。

首先设Ej(·)为所要分析的信号经过EMD获得的第j阶模态分量；ωi(n)为高斯白噪声，i=1，2，…，I；εk为各阶段的信噪比，k=0，…，K。设x(t)为目标信号， CEEMDAN的算法如下[20]。

IMF1(t)表示为分解出的第一个固有模态分量，分解方法等同于EMD，以不同的噪声为基础，将EMD分解过程重复I次，对集合求平均值，得到x(t)的IMF1(t)：

(6)

k=1时，求一阶残差r1(t)得：

r1(t)=x(t)-IMF1(t)

(7)

将通过EMD方法获得的噪声分量E1[ωi(t)]代入到式(7)中后，再次做EMD，将获得的集合平均值定为第二个固有模态分量IMF2(t)：

(8)

k=2，…，K时(K表示模态总数)，求k阶残差rk(t)得：

rk(t)=rk-1(t)-IMFk(t)

(9)

将rk(t)+εkEk[ωi(t)]中的IMF分量提出并求取集合平均值，获得x(t)的第k+1个固有模态分量IMF(k+1)(t)：

(10)

将以上步骤反复至残差无法再进行分解的时候，获得结果就是最终残差R(t):

(11)

可得到目标信号x(t):

(12)

从式(12)可知，CEEMDAN方法能准确地重构原始信号。

设e为希望得到的信号分解相对误差，由相关准则将白噪声[17]确定为：

(13)

式(13)中：ε=σn/σs,表示的是白噪声和原信号的幅值标准差的比；ξ=σh/σs,表示的是高频信息和原信号幅值标准差的比。所以，式(13)可写为

(14)

通过Huang 等[14]总结出的规律计算集合平均次数N，设置好的e值范围则为：

(15)

3 实际数据应用

本文将基于CEEMDAN的希尔伯特变换用于南海北部神狐海域的OBS(海底地震仪)和垂直缆数据中。OBS是一种置于海底接收地震信号的仪器，其接收到的是共接收点记录，因为在工作期间，OBS一直位于海底接收不同炮间距的地震信号，因此一般情况下，OBS不仅能接收反射记录，还对深部地层界面的折射波有所记录。不同于传统的检波器在横向布设以接收信号，垂直缆是一种垂向布设检波器并在底部悬挂海底地震仪的接收装置。理论上垂直缆在海水中是垂直布设的，水听器记录来自震源激发后经反射界面传播回来的地震波。

原始数据(图1)经过CEEMDAN之后会分解出多个IMF分量。如图2和图3所示截取了OBS数据和垂直缆数据的前3个IMF分量，第一阶IMF形成的剖面,在地震数据中是频率最高的一阶，会包含一些细节，也存在高频噪音，所以文中分析不选择IMF1。随着IMF阶数增高,信号幅度会迅速下降,主要成分在前几阶分量上能够得到较好的反应，在IMF3之后的分量形成的剖面都比较模糊，因此本文中选择以数据的第二和第三阶分量为基础来进行分析。

图1 原始数据地震记录Fig.1 Seismic record of original data

图2 OBS数据经CEEMDAN处理的前三阶分量Fig.2 The first three order components of OBS data processed by CEEMDAN

3.1 瞬时振幅

瞬时振幅反映的是振幅包络，一种对于反射强度的度量，反映能量的变化，地层反射的强弱表现得更加明显，在瞬时振幅剖面上，强反射更强，弱反射更弱。一般情况下，来自水合物稳定带底面的反射也大致与海底平行，称为BSR，在地震剖面中表现为与海底轮廓相似的强反射层[21],地层中含有天然气水合物的波阻抗要高于下伏地层中含有游离气地层的波阻抗，在剖面上BSR的特征是振幅较强，含有水合物地层因为有着比较均匀的密度，所以特征是振幅相对较弱，瞬时振幅信息比常规地震剖面更加突出BSR的强反射特征。

图4和图5分别为南海神狐海域天然气水合物的OBS数据的和垂直缆数据的瞬时振幅，其中图4(a)和图5(a)中的为原始数据经过希尔伯特变换得到的瞬时振幅图，图4(b)、图4(c)分别为OBS数据经过CEEMDAN获得的IMF2和IMF3通过希尔伯特变换得到的瞬时振幅图图5 (b)、图5 (c)。从图4和图5 中均可看到BSR，用该方法所成的瞬时振幅图，不论是OBS还是垂直缆数据，分辨率都较图4(a)、图5(a)中有所提高，但是可以看到图4(b)和图4(c)中BSR海底反射不太连续，推测是因为OBS直接沉放在海底的，在一定程度上会受到接收的直达波和海底表层反射得到的反射波的干涉效应；在图5中，垂直缆接收信号的直达波与海底一次反射波能清晰分开，水合物储层BSR效应也很明显，图5(b)和图5(c)显示的BSR清晰，与图5(a)中直接经过希尔伯特变换得到的比较，分辨率更高。

图4 OBS数据瞬时振幅Fig.4 Instantaneous amplitude of OBS data

图5 垂直缆数据瞬时振幅Fig.5 Instantaneous amplitude of vertical cable data

图6 OBS数据瞬时相位Fig.6 Instantaneous phase of OBS data

图7 垂直缆数据瞬时相位Fig.7 Instantaneous phase of vertical cable data

3.2 瞬时相位

瞬时相位是同相轴连续性的度量，一种不考虑振幅强度的相位剖面，地震波穿过有岩性差异的地层时，瞬时剖面上能反映出地震波的相位变化。图6和图7分别为南海神狐海域天然气水合物的OBS和垂直缆数据的瞬时相位图，其中图6(a)、图7(a)为原始数据经过希尔伯特变换得到的剖面，图6(b)、图6(c)和图7 (b)、图7(c)分别为原始数据经过CEEMDAN得到的IMF2和IMF3分量经过希尔伯特变换得到的剖面。从图6和图7中可以看到，原始数据经过希尔伯特变换得到的图6(a)和图7(a)中，相位属性信息复杂，难以分辨有效信息，而经过CEEMDAN处理后的进行希尔伯特变换得到的图6(b)和图6(c)中，也没有明显的追踪到BSR，在瞬时相位属性图中，如果BSR不连续，是能在图中找到，对后续解释有很大帮助；而如果BSR界面平行于地层或相交的角度比较小时，在图中效果不明显。因此，判断在本文的属性图中没有有效识别到BSR的原因，应是本文使用的数据得到的BSR与地层大致平行或交角较小的缘故。

3.3 瞬时频率

瞬时频率表示瞬时相位的时间变化率。地层是否含有气可以通过反射波某一时刻的即时频率即瞬时频率来反映，在AVO属性中，游离气的赋存状态主要通过振幅随偏移距的变化情况来描述，而瞬时频率则通过地震波传播到不同地层的不同频率响应来描述，含游离气的地层在瞬时频率剖面中表现为低频特征(图8、图9)。

图8 OBS数据瞬时频率Fig.8 Transient frequency of OBS data

图9 垂直缆数据瞬时频率 Fig.9 Instantaneous frequency of vertical cable data

如图8、图9所示，分别为南海神狐海域OBS数据和垂直缆数据的瞬时频率图，由于希尔伯特变换是全通滤波，极其容易受到高频随机噪声的干扰，如果接收到的地震信号含有高频噪声，噪声会把经过希尔伯特变换获得的瞬时特点掩盖。在图8(a)和图9(a)中，不论是OBS数据还是垂直缆数据以希尔伯特变换为基础得到的瞬时频率图，无法清楚辨别出高低频信息，难以分辨各层的位置，且无法在图8(a)、图9(a)中识别到BSR，而图8(b)、图8(c)和图9(b)、图9(c)是经过CEEMDAN的IMF2和IMF3再做希尔伯特变换得到的瞬时频率图。相较于图8(a)和图9(a)，虽然依然含有噪音，但是都能识别到BSR。如图8和9所指示的反射能量较强的地方，也就是水合物和游离气共存产生 BSR 相应的位置。 BSR以上表示水合物储集区(该区水合物储层厚度在 30 m 左右)。同时垂直缆与OBS相比，由于垂直缆接收数据具有较高的信噪比和相对较宽的频带，对比图8与图9可知，图9中垂直缆数据的剖面分辨率更高。

4 结 论

通过对天然气水合物实际数据应用的分析可知，瞬时振幅、频率、相位可以很好地识别出天然气水合物的赋存特征，基于CEEMDAN的希尔伯特变换可以用在海洋地震勘探数据OBS和垂直缆数据的分析中，能获得更加清楚的层间信息，为减少地质解释的多解性提供了可参考的依据。

1)瞬时振幅中，垂直缆数据经过处理后分辨率提高，而OBS数据分辨率虽有提高，但BSR有些间断，是因为OBS中包含了纵波分量和横波分量，横波地震记录中水合物的BSR效应不明显，而垂直缆竖直投放在海水中，环境噪音相对较小，因此垂直缆的分辨率要高于OBS；

2)瞬时相位中，无论是只经过希尔伯特变换处理得到的瞬时属性图还是经过CEEMDAN后再进行希尔伯特变换的瞬时属性图，都没有有效追踪到BSR，判断BSR与地层大致平行或交角较小；

3)瞬时频率中，经过处理后的剖面比只经过希尔伯特变换获得的剖面分辨率更高，能识别到BSR的存在；

4)在获得的地震瞬时属性中，对比垂直缆数据和OBS数据，垂直缆数据能获得分辨率更高的信息。