经验模态分解方法(Empirical Mode Decomposition,简称EMD)是由美国国家航空航天局的黄锷博士提出的一种信号分析方法。它依据数据自身的时间尺度特征来进行信号分解,无须预先设定任何基函数。通过EMD分解得到的各个固有模态函数(Instrinsic Mode Function,简称IMF)大多能反映真实的物理过程,大量的实验表明,EMD法是目前提取信号趋势或均值的最好方法(Huang et al,1998)。

EMD分解算法也称为筛选过程,筛选过程主要有两个作用:一是去除叠加波,二是使波形更加对称,其算法本身也是对原始观测信号提取IMF的过程。IMF定义为满足下列两个条件的函数:① 在整个信号中,极值点的数量与过零点的数量相等或至多相差1; ② 信号上任意一点,由局部极大值点确定的包络线和由局部极小值点确定的包络线的均值均为0(张义平等,2008)。

EMD算法的具体步骤如下:首先确定原始信号x(t)的所有极值点,再用三次样条函数将所有的极大值点和极小值点分别进行插值,拟合出原始信号x(t)的上包络线x_max(t)和下包络线x_min(t),然后从原始数据中去除包络线的均值m₁(t):

m₁(t)=(x_max(t)+x_min(t))/2.(1)

从而得到

h₁(t)=x(t)-m₁(t).(2)

重复上述过程R次, 得到第k次筛选的数据h_1k(t):

h_1k(t)=h_1(k-1)(t)-m_1(k-1)(t).(3)

筛选过程终止的标准是通过两个连续处理结果之间的标准差SD的值作为判据的,即

SD=∑^r_t=0|(|h_1(k-1)(t)-h_1k(t)|²)/(h²_1(k-1)(t))|.(4)

在判断筛选过程是否停止时,SD取值一定要谨慎。要避免过于严格的标准,以免导致IMF分量变成纯粹的频率调制信号,造成幅值恒定,同时,也应注意选取过于宽松的标准,导致SD与IMF分量要求相差太远。经验表明,SD值取0.2～0.3为宜,这样既可以保证IMF的线性和稳定性,又可以使IMF具有相应的物理意义。

当h_1k(t)满足SD值的要求时,则h_1k(t)为第一阶IMF,记为c₁(t),即

c₁(t)=h_1k(t).(5)

从原始信号x(t)中减去c₁(t)得剩余信号, 即残差r₁(t):

r₁(t)=x(t)-c₁(t).(6)

再将r₁(t)看作一组新信号重复上述模态分解过程, 经多次运算可得到全部的残差r₁(t):

r₁(t)=r_i-1(t)-c_i(t).(7)

当c_n(t)或r_n(t)小于预定的误差, 或残差r_n(t)成为一单调函数,即不可能再从中提取IMF分量时,就终止模态分析过程。至此,原始信号x(t)可由n阶IMF分量及残差r_n(t)构成:

x(t)=∑ⁿ_i=1c_i(t)+r_n(t).(8)

其中,r_n(t)称为残余函数,代表信号的平均趋势。

通过EMD分解,将原始信号进行平稳化处理,分解出的n个IMF包含了原信号从高到低不同频段的成分,残量则代表了信号的平均趋势。n为原始信号分解出的IMF分量的个数,其取值因原始信号和标准差SD取值有关,在本文的计算中取0.3作为收敛条件。

2010年4月14日青海玉树发生M_S7.1地震。笔者对震中600 km范围内6个台站自2009年1月1日至2010年4月13日的倾斜、应变固体潮观测数据进行EMD分解,提取各台站倾斜、应变的长趋势变化。从图1可以看出,EMD分解得到的残余分量能够非常准确地将倾斜变化拟合出来,说明通过EMD分解提取长趋势变化是可行的,具有明确的物理意义。

图1 德令哈台钻孔倾斜北南分量原始数据与EMD分解长趋势对比时序图(2009年1月1日至2010年4月13日)
Fig.1 Diagram of the observed data and the EMD decomposition residual of bole-hole tiltmeter's NS composition in DeLingha station(Jan.1,2009—Apr.13,2010)

对2009年1月1日至2010年4月13日玉树M_S7.1地震震中600 km以内的6个台站的倾斜、应变固体潮观测数据进行EMD分解,提取各台站、各分量的长趋势变化后,将倾斜南北、东西正交的两分量趋势做倾斜矢量图(图2)。

图2 倾斜矢量空间分布图
Fig.2 Distribution of tilt vector

把近似90°的两个应变分量的长趋势进行相加,可以得到观测点面应变的长趋势变化(张雁滨等,1996)。将玉树地震震中600 km范围内6个台站的面应变长趋势变化从空间上进行展示(图3),面应变数据的起止时间为2009年1月1日至2010年4月13日,图中还给出了玉树—甘孜断裂带。

图3 面应变空间分布图
Fig.3 Distribution of plane strain

(1)利用EMD分解得到固体潮观测的长趋势变化,效果好、数据准确。由长趋势变化得到倾斜矢量图和面应变,其物理意义是明确的。

(2)利用两个正交方向的倾斜潮汐观测长趋势变化,可以获得每个观测点的倾斜矢量图。从玉树地震前距震中600 km以内的6个台站的倾斜矢量变化规律来看,玉树台的倾斜矢量在震前半年出现转折,在“打结”时发震,这与一些历史震例是相符的,比如1975年海城7.3级地震、1988年澜沧、耿马7.6、7.2级地震(张雁滨等,1996,李永振,2010)。

(3)利用两个正交方向的应变潮汐观测长趋势变化,可以获得每个观测点的面应变。在距此次玉树地震最近的三个应变台,即位于甘孜—玉树断裂西北端的格尔木台和中部的玉树台的面应变表现为“受压—平缓—继续受压”,位于东南端的甘孜台的面应变表现为“缓慢张性变化—加速张性变化—缓慢压性变化”这种变化可能反应了这个断裂带震前面应变的变化,但仍有待于和GPS解算结果进行比对。

一个观测点的长趋势潮汐变化能反应多大区域内的潮汐应力变化,是否可以与大尺度的GPS观测或其它断层形变资料进行对比,还需要开展更多的研究来论证。