CN103308940B - 地震剖面的经验模态分解方法 - Google Patents

Abstract

本发明涉及油气地震勘探领域，具体是一种地震剖面的经验模态分解方法，该方法涉及地震剖面二维瞬时属性分析，是一种从地震剖面分解单分量二维固有模态函数的计算方法。在二维经验模态分解过程中，采用时间方向两点邻域的极值检测方法和等式约束的紧框架域稀疏化曲面插值方法，解决地震剖面的极值检测和大规模不规则数据的包络面插值问题。本发明计算结果明显优于现有的二维经验模态分解快速方法，而计算时间大体相当。

Description

地震剖面的经验模态分解方法

技术领域

本发明涉及油气地震勘探领域，更具体地，涉及地震瞬时属性分析中的地震剖面的经验模态分解方法。

背景技术

传统的二维地震剖面复分析技术（A.E.Barnes，Theory of 2D complex seismictrace analysis，Geophysics,Vol.61,No.1，1996，264–272）是对地震剖面直接做二维Hilbert变换（Bidimensional Hilbert Transform，BHT），组成解析信号，进而计算瞬时振幅、瞬时相位和瞬时频率等地震瞬时特征属性，反映了地震波传播过程中的微局部变化情况，可以帮助确定地下岩层的岩性、岩相变化和岩石孔隙中所含流体的性质。但是，这种复分析技术仅适用于单分量数据，而由单道地震信号组成的地震剖面是典型的非线性、非平稳的多分量数据，直接计算地震瞬时属性不符合数学意义和物理意义。

为了得到单分量数据，需要对多分量数据进行二维经验模态分解（BidimensionalEmpirical Mode Decomposition，BEMD），得到一系列具有单一振荡模式分量的二维固有模态函数（Bidimensional Intrinsic Mode Function，BIMF）。由于每个BIMF都满足BHT必要条件，计算瞬时属性变得合理。BEMD计算每个BIMF过程中，计算上下包络面需要进行不规则数据的曲面插值或者拟合，花费大量的运算时间，而且在大多数情况下，还涉及到大规模矩阵和其特征值的计算，计算复杂度高，难以处理大规模数据。为简化计算提高效率，ZhaohuaWu等（Zhaohua Wu,N.E.Huang and Xianyao Chen,The Multi-dimensional EnsembleEmpirical Mode Decomposition Method Advances in Adaptive Data Analysis,Vol.1,No.3，2009，339–372）提出了基于双向一维联合经验模态分解（Ensemble Empirical ModeDecomposition），组合实现二维EEMD（BEEMD），能够快速有效地分解二维数据。这里EEMD是在每次分解前，加入不同的振幅有限的随机白噪声，可以提高分解计算的稳定性。

虽然都是非线性、非平稳的多分量数据，地震剖面不同于一般的图像。一维地震信号是时间方向的地震波形记录，二维地震剖面是由空间方向排列的多道一维地震信号组成，在相邻空间位置的信号波形有一定的错位相似性，表现为连续的起伏波峰和波谷所示的波形同相轴特征。由于地下介质情况的复杂性，地震数据往往不具有高阶光滑性，而且波形同相轴更多表现为时间方向极值的空间连续性。同时，地震数据处理的规模更大，对计算效率有更高的要求。对于地震剖面而言，Zhaohua Wu提出的BEEMD方法，优点是计算量小运算快，但它过于强调这个空间和时间两个方向的极值信息，而忽略这两个方向的差异，丢失很多的波形同相轴信息，不适合地震剖面数据处理。同时，由于分解过程中不断加入白噪声，最后由BIMF计算瞬时属性时，会产生随机噪声背景。

发明内容

本发明的主要目的在于克服了BEEMD方法的不足，提出一种适合地震剖面的快速经验模态分解方法，该分解方法针对地震剖面特点，能够捕捉更多同相轴信息，实现真正的BEMD，简化计算保证效率。

为了实现上述目的，本发明的技术方案为：

一种地震剖面的经验模态分解方法，包括以下步骤：

S1.初始化，r₀=Iori，i=1；给定迭代终止条件中ε的值；Iori表示原始剖面数据，r代表当前残差剖面数据，i为第i个BIMF；

S2.迭代筛分计算第i个BIMF；

S21.初始化：h₀=r_i-1，j=1；

S22.检测h_j-1的二维局部极小值点和极大值点，构造上下包络曲面E_up和E_down；

S23.计算均值曲面m=(E_up+E_down)/2；

S24.h_j=h_j-1-m；

S25.判断当前h是否满足筛分终止条件，若满足，则得到第i个BIMFb_i=h_j，跳转到步骤S3；否则，j=j+1，跳转到步骤S22，继续筛分迭代；

S3.计算当前残差，r_i=r_i-1-b_i；

S4.判断r_i是否满足BIMF终止条件，若满足，则分解结束，记I=i，r_i是最终残余分量；否则i=i+1，跳转到步骤S2；

BEMD计算结果为：

上述采用时间方向两点邻域的极值检测方法检测h_j-1的极值；

上述采用等式约束的紧框架域稀疏化曲面插值方法获取上下包络曲面。

本发明提出的地震剖面的经验模态分解方法为SPBEMD（Seismic ProfileBidimensional Empirical Mode Decomposition，BEMD），是一种改进的BEMD方法，通过时间方向两点邻域的极值检测方法优化极值检测，通过等式约束的紧框架域稀疏化曲面插值方法优化包络曲面差值。两点邻域极值点一般有两种，分别是横向和纵向的，在本发明中为了更好捕捉地震剖面的同相轴特征，采用了时间方向两点邻域的极值检测方法。为了解决包络面的2D不规则数据插值问题，本发明应用等式约束的紧框架域稀疏化曲面插值方法，显著提升BEMD的运算速度。

附图说明

图1为实际地震示意图。

图2为采用两点邻域极大值点检测示意图。

图3为采用四点邻域极大值点检测示意图。

图4为采用八点邻域极大值点检测示意图。

图5为原始地震数据剖面示意图。

图6为BEEMD计算的第1个BIMF的剖面示意图。

图7为采用本发明计算的第1个BIMF的剖面示意图。

图8为本发明的流程图。

具体实施方式

下面结合附图对本发明做进一步描述，但本发明的实施方式并不限于此。

一维IMF已给出了明确的定义，但是对于BIMF尚未有定义。对于一维的IMF，满足4个条件，局部极大值点和极小值点的个数最多相差一个；具有一种振动模式，即在两个过零点之间仅存在一个极大值或者极小值；在任意时刻，上下包络的均值线为零函数或趋于一个零函数；同一信号分解得到的IMF之间应相互正交。由于前两2个条件在分解过程中难以保证，计算过程中BIMF只要满足一维IMF的后2个条件，就作为二维单分量信号的有效近似。为满足这2个条件，类似一维EMD，BEMD通过不断筛分分解BIMF。

如图8所示，本发明具体包括以下步骤：

S1.初始化：r₀=Iori，i=1；给定迭代终止条件中ε的值，Iori表示原始剖面数据，r代表当前残差剖面数据；

S2.迭代筛分计算第i个BIMF；具体为：

S21.初始化：h₀=r_i-1，j=1；

S22.检测h_j-1的二维局部极小值点和极大值点，构造上下包络曲面E_up和E_down；

S23.计算均值曲面m=(E_up+E_down)/2；

S24.h_j=h_j-1-m；

S25.判断当前h是否满足筛分终止条件，若满足，则得到第i个BIMFb_i=h_j，跳转到步骤S3；否则，j=j+1，跳转到步骤S22，继续筛分迭代；

步骤S25的迭代筛分终止条件是考虑限制前后两次迭代剖面数据的标准偏差小于给定值：

$\underset{m=0}{\overset{M}{\mathrm{\&Sigma;}}}\underset{n=0}{\overset{N}{\mathrm{\&Sigma;}}}\frac{{(h_j(m,n)-h_{j-1}(m,n))}^2}{{\left(h_{j-1}(m,n)\right)}^2}<\mathrm{\&epsiv;}$

这里M和N分别表示地震剖面在空间和时间方向的采样点数，而h_j-1和h_j表示前后两次迭代的剖面数据。一般标准偏差小于0.05时，则筛分终止。为了避免筛分迭代发散，同样需要设定一个最大筛分次数Max_Iter作为辅助筛分迭代终止的条件，即如果已达到筛分迭代次数的上限，则停止筛分迭代，输出当前BIMF。）

S3.计算当前残差，r_i=r_i-1-b_i；

S4.判断r_i是否满足BIMF终止条件，即，极值点小于2或近似于平面（或分解出预定最大个数I个BIMF）。若满足，则分解结束，记I=i，r_i是最终残余分量；否则i=i+1，转步骤S2。

最终BEMD计算结果为：

在本实施例中，采用时间方向两点邻域的极值检测方法检测局部极值点，采用等式约束的紧框架域稀疏化曲面插值方法插值计算上下包络曲面。针对地震剖面数据，本实施例提出了以下计算方法具体如下：

1.时间方向两点邻域的极值检测方法

最简单和最直观的局部极值点是一个局部窗口中的最大值点或者最小值点。由于波场离散采样之后变成纵横排列的矩阵，每个点对应一个波场值。所谓的局部窗口一般形如表1所示，选取其中的一个内点相邻的两个点、四个点或八个点与该点的波场值作比较，这就是极值检测的局部窗口法。显然，选取的窗口越小，局部邻域就越小，检测得到的相应的极值点则越多。两点邻域极值点一般有两种，分别是横向和纵向的。为了更好体现地震剖面的同相轴特征，这里选择所定义的两点邻域极值点是时间方向的两点邻域极值点。

根据邻域的概念，给出一幅两点邻域、四点领域和八点领域的极值点定义。

两点邻域极值点，如表1：若(x₀,y₀)处波场值大于或等于其上、下两个相邻位置的波场值，则(x₀,y₀)为剖面的极大值；相应地，若(x₀,y₀)处波场值小于或等于其上、下两个相邻位置的波场值，则(x₀,y₀)处波场值为剖面的极小值。

表1

（i-1，j）
	（i，j）
（i+1，j）

四点邻域极值点，如表2：若(x₀,y₀)处波场值大于或等于其上、下、左、右四个相邻位置的波场值，则(x₀,y₀)为剖面的极大值；相应地，若(x₀,y₀)处波场值小于或等于其上、下、左、右四个相邻位置的波场值，则(x₀,y₀)处波场值为剖面的极小值。

表2

	（i-1，j）
			（i，j-1）	（i，j）	（i，j+1）
	（i+1，j）

八点邻域极值点，如表3：若(x₀,y₀)处波场值大于或等于其上、下、左、右及四个对角相邻位置的波场值，则(x₀,y₀)为剖面的极大值；相应地，若(x₀,y₀)处波场值小于或等于其上、下、左、右及四个对角相邻位置的波场值，则(x₀,y₀)处波场值为剖面的极小值。

表3

（i-1，j-1）	（i-1，j）	（i-1，j+1）
			（i，j-1）	（i，j）	（i，j+1）
（i+1，j-1）	（i+1，j）	（i+1，j+1）

下面分析一个实际地震剖面例子，如图1-4，比较这三种邻域极值检测的效果。图1是一个实际地震图像，分别采用两点、四点和八点邻域进行极大值点检测，结果为图2、图3和图4。从图上可见，满足八点邻域极值条件的点同时也是满足四点邻域极值条件的点，同样，满足四点邻域极值条件的点也是满足两点邻域极值条件的点。两点邻域极值点检测相对于四点和八点邻域，找到的极值点更多，而且更能刻画剖面中波峰和波谷等地震同相轴特征，这对于进一步得到合理准确的包络曲面有很大的帮助。这个实例表明，按时间方向进行的两点邻域极值点检测是适合地震数据剖面的极值检测方法。

2.等式约束的紧框架域稀疏化曲面插值方法

紧框架是冗余的正交基，能够实现重构的系数不唯一。这种冗余性，可以实现稳定的重构表示，即使部分数据缺失时，其不利影响也能够容忍，冗余分解系数可以恢复原缺失信息。所以，紧框架这种冗余性，使得它非常适合进行2D不规则数据插值。考虑到BEMD方法中，检测的极值分布往往是不规则的，相应包络面插值正是2D不规则数据插值问题，可以应用等式约束的紧框架域稀疏化曲面插值方法，即下述分裂的Bregman迭代方法。

设x表示要插值的包络曲面E逐列表示的列向量，而P_Λ表示投影算子，其对角线为1或0值的对角阵，这里1代表此处图像有值，对上包络是那些检测到的不规则极大值点组成的集合，对下包络是那些检测到的不规则极小值点组成的集合，不妨把极大值或极小值都记为f，而把极大值点或极小值点集合都记为Λ。为了简化计算提高效率，考虑等式约束的紧框架域稀疏化曲面插值问题

$\underset{x\&Element;R^{M*N}}{\min }{\left|\right|\mathrm{Bx}\left|\right|}_{1}s.t.P_{\mathrm{\&Lambda;}}x=f$

这里B表示Framelet变换，满足B^TB=I；R^M*N是M*N行的实数值列向量空间；‖·‖₁是定义在R^M*N空间上的l₁范数，对向量x∈R^M*N，而s.t.是subject to的缩写，表示“服从……条件”。由此，上式是将插值问题，等价表示为一个等价的优化问题，即：求解满足插值条件PΛx=f的使得‖Bx‖₁最小的列向量x。

应用分裂的Bregman迭代算法，求解该优化问题，得到迭代关系式

$\left\{\begin{array}{c}{x}^k=f+\mathrm{\&alpha;}(I-P_{\mathrm{\&Lambda;}})B^T(y^{k-1}-q^{k-1})\\ y^k=\mathrm{Pro}{x}_{{\mathrm{\&alpha;}\left|\right|\&CenterDot;\left|\right|}_1}({\mathrm{Bx}}^k+q^{k-1})\\ p^k=p^{k-1}+(x^k-f)\\ q^k=q^{k-1}+({\mathrm{Bx}}^k-y^k)\end{array}\right.$

这里x^k,y^k,p^k,q^k是第k次迭代过程中R^M*N空间列向量，它们的初始值都是R^M*N空间零向量；I是单位算子；B和B^T分别表示Framelet分解和重构；Prox表示逐点计算的临近算子

${\mathrm{Prox}}_{{\mathrm{\&alpha;}\left|\right|\&CenterDot;\left|\right|}_1}\left(x\right)=[{\mathrm{Prox}}_{\mathrm{\&alpha;}|\&CenterDot;|}{x}_1,\&CenterDot;\&CenterDot;\&CenterDot;,{\mathrm{Prox}}_{\mathrm{\&alpha;}|\&CenterDot;|}{x}_{M*N}{]}^t,$

其中上标t表示行向量转置为列向量；Proxx_α|·|x_i就是所谓的软门阈算子

${\mathrm{Proxx}}_{\mathrm{\&alpha;}|\&CenterDot;|}{x}_i=\left\{\begin{array}{cc}{x}_i-\mathrm{\&alpha;},& x_i>\mathrm{\&alpha;}\\ 0,& \left|x_i\right|\&le;\mathrm{\&alpha;}\\ x_i+\mathrm{\&alpha;},& x_i<-\mathrm{\&alpha;}\end{array}\right.$

|·|表示绝对值，而参数α>0是迭代插值算法的收敛条件。实践中根据地震数据包络面的光滑程度，该参数取固定经验值，一般取0.001。至于迭代终止条件，类似于筛分终止条件，考虑限制前后两次插值剖面数据的标准偏差小于给定值：

$\underset{i=0}{\overset{M*N}{\mathrm{\&Sigma;}}}\frac{{(x_i^k-x_{i-1}^k)}^2}{{\left(x_{i-1}^k\right)}^2}<\mathrm{\&epsiv;}$

一般这里标准偏差小于0.05时，则插值迭代终止，将当前迭代列向量x重排，即为插值后的包络曲面E。

上述插值迭代算法，Framelet分解和重构都有快速算法，而临近算子是线性运算，所以总体上计算量非常小，并且只需要存储几个列向量，内容占用也很很少，非常适合大规模地震剖面不规则数据的快速高效包络面插值。

这里，称使用时间方向两点邻域的极值检测方法和等式约束的紧框架域稀疏化曲面插值算法的地震剖面经验模态分解算法为SPBEMD（Seismic Profile BidimensionalEmpirical Mode Decomposition，BEMD）。

3.数值试验

数值试验选用一个共偏移距地震剖面数据，应用SPBEMD，给出第1个BIMF。作为比较，也给出了BEEMD计算的第1个BIMF结果。

图5为原始地震数据剖面，图6为BEEMD计算的第1个BIMF，而图7为SPBEMD计算的第1个BIMF。明显的，图6中BIMF的横向波形连续性比较差，似乎处处都存在波形能量的震荡，个别地方波形能量出现明显破碎，由其相邻层位的信息出现混叠的现象比较严重，比较图6-7中黑色圆圈部分，而且背景中有噪声干扰。而图7中，上述不利现象，都得到了明显的改善，尤其横向波形特征连续性好，波形能量也对比明显，分解结果更理想。

最后，比较采用这两种BEMD方法的计算效率。这里共偏移距地震剖面的规模为6000×320的数据点，计算机CPU是Intel(R)Core(TM)2 Quad CPU，主频2.66GHz，内存4.00GB，编程环境是Microsoft Visual C++.Net 2003，实验中共分解得到5个IMF图像，比较它们的运行时间，BEEMD方法用时154.193s，而SPBEMD方法用时162.238s，基本相当。因此，对于大规模的地震剖面数据的BEMD，这里发明的SPBEMD方法，更具优势。

以上所述的本发明的实施方式，并不构成对本发明保护范围的限定。任何在本发明的精神原则之内所作出的修改、等同替换和改进等，均应包含在本发明的权利要求保护范围之内。

Claims (6)

1.一种地震剖面的经验模态分解方法，其特征在于，包括以下步骤：

S1.初始化，r₀＝Iori，i＝1；给定迭代终止条件中ε1的值；Iori表示原始剖面数据，r代表当前残差剖面数据，i为第i个BIMF；

S2.迭代筛分计算第i个BIMF；

S21.初始化：h₀＝r_i-1，j＝1；

S22.检测h_j-1的二维局部极小值点和极大值点，构造上下包络曲面E_up和E_down；

S23.计算均值曲面m＝(E_up+E_down)/2；

S24.h_j＝h_j-1-m；

S25.判断当前h是否满足筛分终止条件，若满足，则得到第i个BIMFb_i＝h_j，跳转到步骤S3；否则，j＝j+1，跳转到步骤S22，继续筛分迭代；

S3.计算当前残差，r_i＝r_i-1-b_i；

S4.判断r_i是否满足BIMF终止条件，若满足，则分解结束，记I＝i，r_I是最终残余分量；否则i＝i+1，跳转到步骤S2；

BEMD计算结果为：

上述采用时间方向两点邻域的极值检测方法检测h_j-1的极值；

上述采用等式约束的紧框架域稀疏化曲面插值方法获取上下包络曲面。

2.根据权利要求1所述的地震剖面的经验模态分解方法，其特征在于，所述采用时间方向两点邻域的极值检测方法具体为：若某点(x₀,y₀)处的波场值大于或等于其上、下两个相邻位置的波场值，则点(x₀,y₀)为剖面的极大值；相应地，若某点(x₀,y₀)处的波场值小于或等于其上、下两个相邻位置的波场值，则点(x₀,y₀)处波场值为剖面的极小值。

3.根据权利要求1或2所述的地震剖面的经验模态分解方法，其特征在于，所述采用等式约束的紧框架域稀疏化曲面插值方法的具体过程为：

设x表示要插值的包络曲面E逐列表示的列向量，P_Λ表示投影算子，P_Λ是对角线为1或0值的对角阵，1代表此处图像有值，对上包络是那些检测到的不规则极大值点组成的集合，对下包络是那些检测到的不规则极小值点组成的集合，将极大值或极小值都记为f，极大值点或极小值点集合都记为Λ，采用等式约束的紧框架域稀疏化曲面插值方法：

$\begin{array}{ccc}\underset{x\&Element;R^{M*N}}{min}\left|\right|Bx|{|}_1& s.t.& P_{\mathrm{\&Lambda;}}x=f\end{array}$

这里B表示Framelet变换，满足B^TB＝I₁；R^M*N是M*N行的实数值列向量空间；||·||₁是定义在R^M*N空间上的l₁范数，对向量x∈R^M*N，s.t.是subject to的缩写，表示“服从……条件”，上式将插值问题，等价表示为一个等价的优化问题，即：求解满足插值条件P_Λx＝f的使得||Bx||₁最小的列向量x；

应用分裂的Bregman迭代算法，求解该优化问题，得到迭代关系式：

$\left\{\begin{array}{c}{x}^k=f+\mathrm{\&alpha;}(I-P_{\mathrm{\&Lambda;}})B^T(y^{k-1}-q^{k-1})\\ y^k=\mathrm{Pr}o{x}_{\mathrm{\&alpha;}\left|\right|\&CenterDot;|{|}_1}(B{x}^k+q^{k-1})\\ p^k=p^{k-1}+(x^k-f)\\ q^k=q^{k-1}+(B{x}^k-y^k)\end{array}\right.$

式中x^k,y^k,p^k,q^k是第k次迭代过程中R^M*N空间列向量，其初始值均为R^M*N空间零向量；I₁是单位算子；B和B^T分别表示Framelet分解和重构；Prox表示逐点计算的临近算子

${\mathrm{Prox}}_{\mathrm{\&alpha;}\left|\right|\&CenterDot;|{|}_1}\left(x\right)={\&lsqb;{\mathrm{Prox}}_{\mathrm{\&alpha;}|\&CenterDot;|}{x}_1,\mathrm{...},{\mathrm{Prox}}_{\mathrm{\&alpha;}|\&CenterDot;|}{x}_{M*N}\&rsqb;}^t,$

其中上标t表示行向量转置为列向量；Proxx_α|·|x_i就是所谓的软门阈算子

${\mathrm{Proxx}}_{\mathrm{\&alpha;}|\&CenterDot;|}{x}_i=\left\{\begin{array}{cc}{x}_i-\mathrm{\&alpha;},& x_i>\mathrm{\&alpha;}\\ 0,& \left|x_i\right|\&le;\mathrm{\&alpha;}\\ x_i+\mathrm{\&alpha;},& x_i<-\mathrm{\&alpha;}\end{array}\right.$

|·|表示绝对值，参数α＞0是迭代插值算法的收敛条件，达到插值迭代终止条件，则插值迭代终止，将当前迭代列向量x重排，即为插值后的包络曲面E。

4.根据权利要求3所述的地震剖面的经验模态分解方法，其特征在于，所述插值迭代终止条件为：前后两次插值剖面数据的标准偏差小于给定值ε2：

$\underset{i=0}{\overset{M*N}{\&Sigma;}}\frac{{(x_i^k-x_{i-1}^k)}^2}{{\left(x_{i-1}^k\right)}^2}<\mathrm{\&epsiv;}2.$

5.根据权利要求4所述的地震剖面的经验模态分解方法，其特征在于，所述S25的筛分终止条件具体为：

前后两次迭代残差剖面数据的标准偏差小于给定值ε3：

$\underset{m=0}{\overset{M}{\&Sigma;}}\underset{n=0}{\overset{N}{\&Sigma;}}\frac{{\left(h_j\right(m,n)-h_{j-1}(m,n\left)\right)}^2}{{\left(h_{j-1}\right(m,n\left)\right)}^2}<\mathrm{\&epsiv;}3$

其中M和N分别表示地震剖面在空间和时间方向的采样点数，h_j-1和h_j分别表示前后两次迭代的残差剖面数据；

或筛分迭代次数达到所设定的最大筛分次数Max_Iter。

6.根据权利要求5所述的地震剖面的经验模态分解方法，其特征在于，所述步骤S4的BEMD终止条件是：极值点小于2或近似于平面或分解出预定最大个数I个BIMF。