艾文强,范锦彪,王 燕

(仪器科学与动态测试教育部重点实验室, 太原 030051)

初至拾取指对传感器最先接收到的有效波进行拾取。实现高精度定位，在信噪比低的情况下提取可靠的首波到达时间是亟待解决的问题。目前常用的去噪方法是傅里叶变换和小波分析法。采用傅里叶去噪时待分析的信号需要具有一定的线性、周期性或平稳性。对频率随时间变化的非平稳信号，傅里叶变换只能给出整体效果，不能完整地掌握某一时刻信号的本质特征，消噪效果不好[1]。采用小波分析算法快速性、时频性较好，应用广泛。但同时，小波基函数、分解层数和小波变换的阈值选择都对去噪效果有影响，选取原则通常取决于经验，这使得小波变换在分析信号时的自适应性较差。

提出联合互补集合经验模态分解(Complementary Ensemble Empirical Mode Decomposition，CEEMD)算法[2]和改进阈值去噪的方法。CEEMD是近几年来发展起来的一种新型自适应信号时频分析法，根据信号的特点它能自主分解出本征模态函数(Intrinsic Mode Function，IMF)，这是一种适用于分析非线性非平稳信号的方法。常用的阈值去噪方法有硬阈值法、软阈值法。本文提出了改进的阈值去噪法。新方法克服了硬阈值函数不连续及软阈值函数中小波系数估计和原系数之间具有固定偏差的缺点。其基本流程为对信号求其包络，并根据信号的均方差或期望设置合理的包络值，去噪部分为包络小于包络值的部分。对信号进行CEEMD分解后再用新的阈值函数进行改进阈值去噪，既可以保证时频性又可以很好地解决自适应问题[3]。

初至提取算法，可以分为以下几类：① 时域分析法，如时窗平均能量法、分形维数法、自回归模型法、相关法等。时域分析法简单易行，但对信噪比较低且起跳不明显的爆炸引起的地震波效果不好。② 时频分析法，如Fourier、小波变换、S变换，与时域分析法相比，这种方法有很强的的抗噪能力，但受窗长度等参数的影响，一致性不好。③ 综合分析法，如人工神经网络法。综合分析法是目前精度较高、一致性较好的一种方法，但耗时长，实现复杂。④ 本文采用的能量比值算法，这种算法结合初至波识别算法可以快速、准确拾取初至点[4]。

图1为初至点提取总体流程示意图。

图1 初至点提取总体流程示意图

1 EMD算法与CEEMD算法

1.1 EMD算法

经验模态分解(Empirical Mode Decomposition，EMD)可以被理解为能通过信号极值时间特征量来表示的一种时空滤波算法[5]。经验模态分解根据信号局部时间特征尺度将信号从小到大分解，并获得有限个频率从大到小的IMF。每个IMF必须满足如下两个条件：

1) 在整个信号中，极值点的数量和过零点的数量相差不超过1；

2) 在任意时刻，上包络和下包络的平均值为0。

通常，实际信号是不符合上述条件的复杂信号。因此，Huang进行了以下假设[2]：任何信号都是由若干本征模态函数组成的；每个本征模态函数可以是线性的，或非线性的，每个IMF的局部零点和极值点数是相同的，并且上包络线和下包络线关于时间轴是局部对称的；任何情况下，信号都可以包含几个IMF，如果各IMF之间相互混叠，就组成了复合信号。

EMD分解的基本过程为：

1) 找到信号所有的极大值点和极小值点，并通过三次样条函数拟合信号的最大包络线e+(t)和最小包络线e-(t)。上包络线和下包络线的平均值被视为原信号的平均包络，如式(1)所示。

(1)

(2)

(3)

一般情况下，上包络和下包络的平均值无法为零，所以当符合式(3)时，就判定满足要求。其中ε为筛分门限，一般取值在0.2到0.3之间。

(4)

3) 用原信号减去c1(t)，得到新信号r1(t)。

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

(5)

对r1(t)重复得到c1(t)的过程，得到第二个IMF分量c2(t),如此反复进行，直到第n阶IMF分量cn(t)或余量小于预设值；当剩余量是单调函数或常量时，分解停止[6]。

4) 原信号经分解后得到：

(6)

1.2 CEEMD算法

CEEMD分解基于EMD分解，包括以下步骤：

1) 将n组正负白噪声加到原始信号，生成的添加噪声后的信号表示为：

(7)

2) 添加噪声的信号的EMD分解产生一组本征模态函数分量。

3) CEEMD分解得到的第j个IMF分量表示为：

(8)

与EMD相同，原始信号的噪声最多的部分位于低阶IMF中。该算法需要添加辅助噪声幅值k和对数N两个参数，一般当N取100时，k取0.01～0.10。

2 CEEMD分解后的小波阈值去噪及信号重构

在CEEMD将信号分解后，获得频率从高到低的有限个本征模态函数，高阶IMF对应信号的低频部分，低阶IMF为信号高频部分。基于CEEMD去噪的主要思想是，对大多数噪声污染的信号，其主要能量集中在低频带，频带越高，其包含的能量越少。因此，一定存在某个本征模态函数，使得对于该函数值后的IMFk+1中信号主导，其前k个IMF中噪声主导，基于CEEMD滤波去噪的主要目的就是寻找到这个IMFk。

自相关是一个信号于其自身在不同时间点的互相关，非正式地来说，它就是两次观察之间的相似度对它们之间的时间差的函数。其定义式为[7]：

R(s,t)=E(x(s)*x(t))

(9)

按式(9)分别计算噪声信号和正弦信号的自相关函数。从图2、图3可以看出，噪声信号的自相关函数在零点处具有自相关函数的最大值，且自相关函数值在其他点处迅速衰减到较小的值；正弦信号，它的最大值也出现在零点处，但由于信号之间的相关性，在其他点处自相关函数值不会快速衰减到很小，而是随时间差而变化，两种信号的自相关函数有着明显的区别。根据该特征，每个IMF自相关处理可以找到噪声主导和信号主导的分界本征模态函数[8]。

图2 噪声信号及其自相关函数

图3 正弦信号及其自相关函数

软硬阈值方法虽然得到了广泛的应用，在一些去噪方法中也取得了较好的效果，但硬阈值函数存在不连续点，容易出现附加振荡影响结果光滑性；软阈值使信号连续性很好，比硬阈值函数结果更光滑，但因为与原始信号存在恒定的偏差，严重影响重构信号精度。

本文提出了一种新的阈值函数：

(10)

2) |wij|≥λ时，函数是高阶可导的，方便进行数学处理。

3 能量比值算法原理及实现

能量比方法是基于地震记录中地震波的瞬时特征拾取初至的。地震信号在到达时间之前和之后有明显差异，即信号的振幅和能量在到达时间之前保持相对稳定；当地震信号到达时，信号的振幅和能量比到达之前出现了不同水平的突增[9]。对于单一通道信号X(i)，前时窗和后时窗长度为M，总窗长为2M，前时窗和后时窗的能量值分别计算，则后窗能量与前窗能量的比为：

(11)

为了提高结果的稳定和提高可靠性，在式(11)的基础上增加一个稳定因子S，S的表达式为：

(12)

式中，L为通道信号的长度。则式(11)可改写为：

(13)

当信号没有到达时，R(i)的值接近于1；当信号到达时，R(i)的值会突然有很大变化。根据实际情况设定一个阈值P，当R(i)超过这个阈值P时就认为该时窗内是初至波，从而确定初至点。通常还需设定另外一个振幅阈值A，有且仅有R(i)和振幅都超过预设阈值时，才认为识别到初至波。

4 模型建立及测试

1) 模型建立

地震信号在地下传播时部分鞥量会被吸收，地震信号衰减，本文采用下式合成衰减地震信号。

(14)

其中:S(ω)是的Fourier变换，ω为频率，k为对应波数，r为传播距离。

子波选用零相位雷克子波，零相位子波的表达式见式(15)。

w(t)=e-(2πfm/γ)2t2cos2πfmt

(15)

其中:fm代表子波的中心频率，γ代表子波宽度[10]。

将零相位子波按式(13)可合成一维地震衰减记录如图所示。从图4合成的一维地震记录观察，可以发现初至点很明确，但是实际地震记录往往还有很多噪声，因此为了更真实的验证本文方法的有效性，对合成记录加噪，得到真实地震记录的模拟记录，如图5。

图4 合成的一维地震记录

2) CEEMD分解及阈值去噪处理

对图5的模拟地震记录进行CEEMD分解，以便后边进一步处理。其分解图如图6所示，共分解为8个IMF分量和一个趋势。其中阶数小的IMF对应于信号的高频成分，阶数大的IMF对应于信号的低频成分。

图5 模拟地震记录

为了找到噪声起主导作用IMF和信号起主导作用IMF的分界点IMFk，求解出各IMF的自相关函数(如图7、图8)。

图6 模拟地震记录的CEEMD分解图

图7 前4个IMF自相关函数图

图8 后4个IMF自相关函数图

根据信号与噪声的自相关特点不能直接确定要处理的分量。为了更精确的确定要处理的IMF分量，分别求信号与自相关结果的能量谱 中心频率，通过对信号的能量谱中心频率与IMF分量的能量谱中心频率进行对比，准确得出要确定的IMF分量[11]。

如图9，可以直观的判定前3个IMF为噪声起主导作用，对IMF1～IMF3分量进行改进阈值去噪。图10为最终去噪后的结果。

图9 自相关结果对应频域图

图10 重构信号

3) 初至拾取

对原始信号和经过去噪处理的信号分别求初至，求解结果如表1所示。从表1中可以得出，基于CEEMD和小波阈值去噪的方法对于初至的精确拾取有明显的改善，该方法能比较精确拾取低信噪比信号的初至。

表1 两种方法拾取时间

5 实际应用及效果分析

为验证算法实际情况的可行性，进行了人工地下爆破模拟地震实验。将3 kg TNT埋设在地下5 m起爆，地震波采集装置的采样频率为20 kHz，采样总时间为10 s。图11为地下起爆点起爆后，传感器获取的局部震动时域信号图。

图11 时域信号图

为方便计算只对前5 s信号进行CEEMD分解，分解图如图12，信号被分解为9个IMF分量及一个趋势项。

图12 实测信号CEEMD分解图

对前四个IMF求解自相关结果如图13。根据自相关结果，不能很直观的得出需要对哪几个分量进行改进阈值去噪。为了更准确的判断对自相关结果，做傅里叶变换，通过能量谱中心频率进一步对需要处理的IMF分量进行判断，如图14。通过计算主频的方法可以精确的确定需要对IMF1、IMF2、IMF3进行改进阈值去噪处理。

图13 前4个IMF分量自相关图

图14 自相关结果对应频域图

去噪后的实测信号重构图如图15，表2所列数据为直接通过能量比值法拾取初至时间与去噪后拾取初至时间的结果。

图15 实测信号重构图

算法时间/ms直接采用能量比值算法3 059.1本文方法3 047.2

6 结论

1) 联合CEEMD与改进阈值去噪再利用改进能量比值法可精确提取初至时间。

2) 通过仿真信号和试验采集实际信号，验证了算法的稳定性、有效性和高精度性，可以对爆炸产生的地震信号进行初至拾取工作。