经验模态分解(EMD)在地球物理资料中的应用(附MATLAB程序)

摘  要
经验模态分解(EMD)是由Huang等人提出的一种新的分析非线性、非平稳信号的方法。本文研究经验模态分解原理及其在地球物理资料中的应用。首先研究经验模态分解的基本原理和算法，对地球物理资料（地震资料，重磁资料）进行EMD分解试验分析，然后研究基于EMD的Hilbert变换原理及其在提取地震属性信息中的应用，对实际地震时间剖面和时间切片进行EMD时频分析试验。
本文的方法研究和数据试验分析表明：经EMD分解变换得到的IMF序列是直接从原始时序数据中分离出来的,事先无需确定分解阶次,能更好反映原始数据固有的物理特性,每阶IMF序列都代表了某种特定意义的频带信息；EMD分解获得的IMF序列具有稳态性,对IMF进行Hilbert变换，就可以得到单个固有模态函数的瞬时振幅、瞬时相位和瞬时频率，这些信息可以清楚的显示信号的时频特征；EMD分析方法用于分解地球物理资料和作时频分析是有效的。
关键词：经验模态分解；地球物理；Hilbert变换；固有模态函数；时频分析
 
ABSTRACT
Empirical Mode Decomposition(EMD), which was developed by huang, is a new method to analyse nonlinear and nonstationary signals. In this paper, we study the theory of EMD and its applications in handling geophysical data. Firstly, we introduce the theory and the Methodology about EMD ,then we will use this method to analyse the geophysical information, including the gravity anomaly data and seism’s data. Based on the EMD, we will study the theory of the Hilbert transform, and then use it to obtain the images,from which we can deal with the seism’s slice by time- frequency analysis in order to distill the seism’s information.
The studying of EMD and the data testing in this paper indicate: intrinsic mode functions(IMF) is comes from the original signal by the EMD, in this course, we need not fix on the Decomposition number and would not influenced by some men’s factors. Every intrinsic mode function stand for some given information and can reflect the intrinsic physical Information very well. Meantime, take the IMF for Hilbert transform, then we can get the IMF’ frequency and the amplitude, which can show the signal’s characteristics very clearly. So it is very useful to use the EMD to study the geophysical data.
Key words: Empirical Mode Decomposition; geophysics; Hilbert transform;
intrinsic mode function; time- frequency analysis

目前，信号处理[1]方法主要有傅立叶变换、小波变换等，这些方法在处理线性、平稳的信号时，具有简单、效率高等优点，然而，在实际的地球物理数据测量中，许多信号是非平稳的或非线性的，如果在处理这些信号时，仍假设数据是平稳或线性的，会使我们得到错误的分析结果，为了正确有效地分解非平稳信号，本文研究一种新的方法，即经验模态分解(EMD) [2]，该方法是由Huang等人提出的一种新的分析非线性、非平稳信号的方法。
经验模态分解（EMD）方法可将信号分解为若干个固有模函数(Instrinsic Mode Function，IMF)，这些IMF是完备的并且正交的。由于EMD方法是依据数据本身的时域信息进行的时域分解，得到的IMF的个数通常是有限的和平稳的，基于这些IMF分量进行的Hilbert变换，得到IMF的瞬时振幅，瞬时频率，以及瞬时相位，从这些频谱图上我们可以得到真实清楚的物理信息，因此，其Hilbert谱也能够准确反映出信号能量、频率在空间或时间尺度上的分布，在非线性和非平稳过程的分析中具有很高的应用价值。
由于经验模态分解法主要用于分析非线性非稳定性数据，所以它在各种资料处理中应用非常广泛。在过去的几年,EMD方法已经被广泛用于非线性的海洋波动数据分析[3]、地震信号和结构分析[4]、位场数据分析[5]、机械故障诊断[6]、医学信号[7]、桥梁和建筑物状况监测的分析、太阳辐射的变化分析等领域的研究。也可以将经验模态分解法扩展至二维，应用于图象处理[8]，将图像分解为一系列不同尺度的信息，为深入处理图象提供足够的不同尺度的信息，在实际应用中取得了很好的效果。
本文研究经验模态分解原理及其在地球物理资料中的应用。本文首先研究经验模态分解的基本原理和算法，对地球物理资料（地震资料，重磁资料）进行EMD分解试验分析，然后研究基于EMD的Hilbert变换原理及其在提取地震属性信息中的应用，对实际地震时间剖面和时间切片进行试验分析。

本文首先研究了EMD信号分析方法的基本原理和算法，对地球物理资料（地震资料，重磁资料）进行EMD分解试验分析，然后研究了基于EMD的Hilbert变换原理及其在提取地震属性信息的应用，对实际地震时间剖面和时间切片做了试验分析。
本文的方法研究和数据试验分析表明：
1、经EMD分解变换得到的IMF序列是直接从原始时序数据中分离出来的,事先无需确定分解阶次,不存在机械分解。因此IMF序列能更好反映原始数据固有的物理特性,其分解是客观的,内在的和自适应的。每阶IMF序列都代表了某种特定意义的频带信息,给实际应用与解释工作带来了方便；
2、原始信号经过EMD分解获得的IMF序列具有稳态性,对IMF进行Hilbert变换，就可以得到单个固有模态函数的瞬时振幅图，瞬时频率图，瞬时相位图，这些图可以清楚的显示信号的时频特征；
3、EMD分析方法用于分解地球物理资料和作时频分析是有效的。







目    录
第一章  绪  论 1
第二章  经验模态分解(EMD) 2
第一节  EMD基本原理 2
第二节  EMD算法流程 3
第三节  EMD示例 5
第三章  EMD分解地球物理信号 8
第一节  EMD分解位场资料 8
一、EMD分解一维重力异常数据 8
二、EMD分解二维重力异常数据 9
第二节  EMD分解地震资料 13
一、EMD分解一维地震数据 13
(毕业设计)
二、EMD分解地震时间剖面 14
三、EMD分解地震时间切片 17
第三章  基于EMD的Hilbert变换与地震信号属性信息 19
第一节  基于EMD的Hilbert变换原理 19
第二节  提取地震信号属性信息 20
一、提取地震时间剖面属性信息 20
二、提取地震时间切片属性信息 22
结  论 23
致  谢 24
参考文献 25
附录：经验模态分解MATLAB程序 26

附录：经验模态分解MATLAB程序
经验模态分解程序
function [imf,ort,nbits] = emd(varargin);

[x,t,sd,sd2,tol,display_sifting,sdt,sd2t,ner,nzr,lx,r,imf,k,nbit,NbIt,MAXITERATIONS,FIXE,FIXE_H,MAXMODES,INTERP,mask] = init(varargin{:}); 

if display_sifting
  figure
end

% maximum number of iterations
% MAXITERATIONS=2000;

%main loop : requires at least 3 extrema to proceed
while ~stop_EMD(r) & (k < MAXMODES+1 | MAXMODES == 0) & ~any(mask)

 % current mode
 m = r;
 

 % mode at previous iteration
 mp = m;

   
 if FIXE
  [stop_sift,moyenne] = stop_sifting_fixe(t,m,INTERP);
 elseif FIXE_H
  stop_count = 0;
  [stop_sift,moyenne,stop_count] = stop_sifting_fixe_h(t,m,INTERP,stop_count,FIXE_H);
  stop_count = 0;
 else
  [stop_sift,moyenne] = stop_sifting(m,t,sd,sd2,tol,INTERP);
 ...