引言

随着地震科技的不断发展，监测仪器的精度和广度不断提升，地震波形中包含的信号也愈加丰富。这使得一些弱震级的地震事件很容易被噪声淹没，难以提取地震信息和分析震源机制，限制了台站密度较小地区的弱震级地震分析。因此，借助数字信号处理技术和计算机技术，运用数学方法对低信噪比的地震数据进行去噪处理，从而提升信号质量，也是地震数据分析的重要的基础性工作之一。在日常地震观测中，地震记录中的背景噪声主要来自地脉动、仪器设备的自带噪声、台风、海浪潮汐等，传统的地震数据去噪方法在进行大量数据去噪时会面临两个问题，一是逐个事件进行噪声特征分析和去噪极大地增加了数据处理的工作量，二是设置统一参数进行去噪导致许多有效数据也被剔除。因此，研究一种设置简便且能够在较强噪声背景下保留波形特征的地震数据去噪方法，具有十分重要的科学意义和实用价值。

地震数据去噪技术一直受到国内外各研究院所、专家学者的高度重视，许多不同领域的方法已经被应用到地震数据的去噪处理中，并在实际应用中取得了较好的效果。例如：f-k （频率-波数）分析法（万光南，2014；牛永效，2017）主要是利用傅里叶变换将原始波形信号从时间-空间域（t-x）转换到频率-波数域（f-k），再根据波形信号中地震数据与噪声信号的速度与频率差异，计算它们在频率-波数域的分布特征，并对噪声集中区域进行剔除，从而实现去噪；小波变换去噪技术（蔡剑华，肖晓，2015；Liu et al，2016；Lu et al，2019）则是利用小波变换将信号按照频率高低分解为不同频率的波形数据，通过去除噪声信号所在的频带，对地震数据进行重构，实现信号去噪；Ｓ变换去噪方法（陈学华等，2008；孙月，2012；Wang et al，2016）主要是同时结合傅里叶变换和小波变换的特点，通过在传统傅里叶变换的基础上加时间窗，计算地震数据与噪声信号的时间频率特征，进行构造阈值从而对地震数据去噪；基于深度学习的去噪方法（韩卫雪等，2018）利用数字图形处理、深度学习、神经网络等技术，分别计算地震数据与噪声信号的多种不同特征，利用特征参数对神经网络进行训练，优化神经网络节点，从而识别和去除噪声信号；基于模态分解技术（empirical mode decomposition，缩写为EMD）的去噪方法（杨凯，刘伟，2012；张杏莉等，2018）是利用EMD技术根据波形信号本身特点进行自适应分解，找出噪声所在本征模函数（intrinsic mode function，缩写为IMF）并予以剔除，从而实现去噪。

上述方法在信号去噪中均取得了较好的效果，但是由于不同台站记录的噪声信号所在的频带、波数、振幅等特征存在差异，在处理不同类型事件时需要不停地调整各种特征参数才能达到较好的去噪效果。鉴于此，本文提出一种基于自适应噪声完全集合经验模态分解（complete ensemble empirical mode decomposition with adaptive noise，缩写为CEEMDAN）算法和Hurst指数（Hurst，1951）的地震数据去噪方法。该方法首先利用CEEMDAN算法的特点，根据地震数据自身频率，自适应地将地震数据分解为一系列IMF波形分量；然后，计算各IMF分量的Hurst指数，并以此作为识别地震数据和噪声信号的特征判据；最后，对地震数据所在的IMF进行重构实现去噪。

1.   相关算法

1.1   本文算法设计

一般来说，一套完整的数据去噪算法主要由三套模块组成：① 波形分解模块，用于将目标数据分解成细小分量；② 成分判定模块，利用已提取出的特征参数对分量类型进行判定；③ 波形重构模块，将识别出的有效数据分量进行重构。本文在设计三个模块时重点关注了以下三点：① 由于不同类型地震事件具有不同的优势频率，简单地根据频段进行分解难以将全部有效数据分解出来；② 考虑到进行大批量数据处理时的便捷性，特征判据应能快速、有效地将地震数据与噪声信号区分开来，避免因事件类型不同需要反复调整参数的情况；③ 从现实应用角度出发，设置尽可能少的参数以确保应用更加便捷。

鉴于此，本文采用CEEMDAN算法与Hurst指数相结合的方式进行数据去噪处理。CEEMDAN算法是在EMD算法基础上改进而来，该方法能够根据波形数据自身的频率特征，将事件按照频率高低逐步分解，从而确保数据分解的精确性和有效性。Hurst指数是一种用于描述信号长期持续性的重要指标（Zhang et al，2015），由于在地震观测中地震数据与噪声信号的持续性存在明显差异，因此使用该指数作为特征判据时不需要进行反复修改调整。图1给出了本文算法流程框图。

图  1  本文算法流程框图

Figure  1.  A flowchart of the method proposed in this paper

下载: 全尺寸图片 幻灯片

1.2   CEEMDAN算法

由于传统的EMD分解方法存在模式混叠问题，即不同频段和长度的波形叠加共存于一段数据中，会造成分解后同一IMF分量中存在不同频率的波形混叠。该问题不仅降低了分解的效率，还降低了滤波性能。为了解决该问题并提高精度，本文采用由Torres等（2011）提出的CEEMDAN方法。CEEMDAN方法的计算过程如下：

计算IMF1分量。在原始波形信号x（t）中分别加入正负成对的高斯白噪声$E_1 （ n^{ （ i ） } ） $ （i＝1，2，···，N）构成信号集合{x（t）＋E₁（n^（i））}，利用EMD分别对每一个信号进行分解，用分解得到的第一个IMF分量组成集合{$h_1^{ （ i ） }$}，对该集合求平均得到IMF1的分量h₁ ［ 式（1） ］ ，并计算IMF1的残余r₁ ［ 式（2） ］ 。

同理计算IMF2分量，得到IMF2分量h₂，并计算残余r₂。

依次类推，可计算得到原始信号的IMFm分量h_m和残余r_m分别为

重复上述步骤，直到残余信号为单调函数，完成分解。此时原始信号x（t）的分解结果为

1.3   Hurst指数

本文使用聚合方差法计算Hurst指数。具体步骤为：对于原始信号序列x（t），t＝1，2，···，N，用长度为$ M $的时间窗，将x（t）划分为X_k（t）（k＝1，2，···，N/M）个区间，分别计算每个区间$ X_k （ t ） $的平均值。

式中N/M不为整数时向上取整，各区间相互不重叠，最后一个时间窗的长度以实际的数据量为准。

计算各区间$ X_k （ t ） $平均值的样本方差，即

分别计算M＝1，2，···，N/4时对应的${\rm{Var}}_M （ t ） $，并在双对数坐标系上用最小二乘法拟合$ ［ M， {\rm{Var}}_M （ t ） ］ $，即

拟合后直线的斜率A即为Hurst指数$ h $。

2.   波形去噪

2.1   基于CEEMDAN的信号分解

为了验证本文算法的有效性，从近年来广西地震台的地震记录中挑选四种不同类型且背景噪声较大的事件波形进行去噪处理，事件具体信息列于表1。四个事件的波形均由宽频带地震计记录，采样率为100 sps，事件长度为120 s。本文选择垂直向波形进行分析，图2给出了四个事件的垂直向波形。从图中可以看出，由于几个事件的震级小、距离远，而背景噪声又比较强，事件的震相较难识别。

表  1  本文所选事件的信息

Table  1.  The information of events selected in this paper

事件序号	事件类型	发震日期	震中位置	记录台站	震中距/km	ML
a	天然地震	2021-04-28	广西百色	FUN	73	1.3
b	天然地震	2017-11-25	广西灵川	YOF	148	1.9
c	人工爆破	2021-04-04	广西南宁	DAQ	88	2.0
d	山体坍塌	2020-06-23	广西桂林	GUL	36	2.3

下载: 导出CSV 

| 显示表格

图  2  事件垂直向原始波形

Figure  2.  Vertical original waveforms of the events

下载: 全尺寸图片 幻灯片

在使用CEEMDAN算法时，为了避免添加过高的噪声对原始信号产生干扰，且过大的运算量会导致运行时间过长，将CEEMDAN算法的附加噪声振幅定为0.01，最大筛选迭代次数设为300。图3给出了四个事件的CEEMDAN分解结果，可以看出，地震波形主要集中在IMF1−IMF5，其余的则是背景噪声信号。

图  3  CEEMDAN算法分解结果

（a） 百色地震事件波形；（b） 灵川地震事件波形

Figure  3.  Decomposition results of CEEMDAN

（a） Baise earthquake waveform；（b） Lingchuan earthquake waveform

下载: 全尺寸图片 幻灯片

图  3  CEEMDAN算法分解结果

（c） 南宁爆破事件波形；（d） 桂林坍塌事件波形

Figure  3.  Decomposition results of CEEMDAN

（c） Nanning blasting waveform；（d） Guilin collapse waveform

下载: 全尺寸图片 幻灯片

2.2   基于Hurst指数的IMF选择

Hurst指数作为一种用于描绘时间序列波形长期记忆性的重要参数，可以反映事件波形震动的持续性特征。根据式（10）计算得到的$ h $值在表示波形的持续性方面具有以下特征：① 当h＝0.5时，说明该波形具有随机性，表现出随机震动；② 当h＜0.5时，说明该波形具有反持续性；③ 当h＞0.5时，说明该波形具有长期记忆性。

在日常的地震监测中可以发现，地震记录的背景噪声是持续存在的，且占据了大部分的波形记录，表现出长期记忆性，而地震波相较于背景噪声则具有频域丰富、持续时间短的特点。因此，Hurst指数可以作为检验分解后的IMF分量是否为地震数据的重要指标。

根据Hurst指数的定义，本文将h＜0.5的分量作为地震数据的IMF分量，将h≥0.5的分量作为噪声IMF分量。表2给出了四个事件IMF分量Hurst指数的计算结果。根据Hurst指数不同，将h＜0.5的IMF分量作为地震数据进行重构。图4展示了四个事件的去噪效果。结果表明，利用本文方法去噪后，波形震相更加清晰且保留了波形特征。

表  2  表1中四个事件的IMF分量的Hurst指数

Table  2.  The Hurst exponents of different IMF components for four events in Table 1

分量序号	事件a	事件b	事件c	事件d
IMF1	0.302 5	0.294 6	0.258 9	0.359 4
IMF2	0.322 6	0.407 5	0.425 5	0.325 4
IMF3	0.307 0	0.287 5	0.356 1	0.466 3
IMF4	0.418 9	0.333 7	0.306 3	0.487 1
IMF5	0.483 8	0.384 6	0.515 9	0.361 1
IMF6	0.655 8	0.640 8	0.612 2	0.647 8
IMF7	0.686 0	0.686 2	0.658 9	0.683 1
IMF8	0.763 8	0.792 9	0.796 6	0.653 0
IMF9	0.806 2	0.877 7	0.870 9	0.864 5
IMF10	0.969 3	0.916 5	0.992 9	0.908 0
IMF11	0.973 3	0.966 8	1.007 0	0.955 1
IMF12	0.985 1		1.009 2	1.004 0
IMF13	0.994 9

下载: 导出CSV 

| 显示表格

图  4  使用本文方法所得四个事件的去噪结果

左侧为原始波形，右侧为去噪后波形。（a） 百色地震波形；（b） 灵川地震波形；（c） 南宁爆破事件波形；（d） 桂林坍塌事件波形

Figure  4.  Denoising results of the four events by the method proposed in this paper

The left column stands for original waveforms and the right column stands for denoising waveforms. （a） Baise earthquake waveform；（b） Lingchuan earthquake waveform；（c） Nanning blasting waveform；（d） Guilin collapse waveform

下载: 全尺寸图片 幻灯片

3.   去噪效果分析

为了验证本文去噪方法的有效性，将文本方法结果与传统的带通滤波器、小波包变换及EMD等去噪方法进行对比分析。

图5为使用上述三种方法及本文方法分别对灵川地震事件和桂林坍塌事件进行去噪后的结果，其中带通滤波器去噪方法采用的是巴特沃斯（Butterworth）带通滤波器，小波包变换去噪方法采用db4基函数进行四层分解去噪。从图5可以看出：虽然传统的方法能够较好地滤除背景噪声，但是在去噪的同时也会降低该频段地震数据（如短周期面波）滤除，导致地震波产生变形；而本文的方法在去除噪声的同时，还能够较好地保存地震数据不同频段的波形特征。

图  5  不同方法对灵川地震事件（a）和桂林坍塌事件（b）的去噪结果

Figure  5.  Denoising results of Lingchuan earthquake （a） and Guilin collapse （b） by different methods

下载: 全尺寸图片 幻灯片

为了能够对几种方法的去噪效果进行量化评价，本文采用质量因子Q值（Chen，Sacchi，2015）和均方根误差（RMSE）作为评价去噪效果优劣的指标。

式中，s_o（i）为原始波形信号，s_d（i）为去噪后的波形信号。Q值越高而RMSE越低，说明去噪的效果越好。表3给出了上述三种方法与本文方法去噪效果的质量因子Q值，可见：与带通滤波器和小波包变换等传统去噪方法相比，对于低信噪比的波形（表1事件b），本文方法的Q值比传统方法Q值高32%，而RMSE降低3%；对于高信噪比的波形（表1事件d），本文方法的Q值则高出648.3%，而RMSE降低33%。这表明相对于传统方法，本文方法的去噪性能更好。

表  3  不同去噪方法的质量因子Q值和均方根误差RMSE

Table  3.  Quality factor Q and root-mean-square error RMSE of different denoising methods

去噪方法	灵川地震事件			桂林坍塌事件
	Q值	RMSE		Q值	RMSE
带通滤波器	0.723 9	0.380 1		1.360 1	0.076 1
小波包变换	0.718 5	0.375 6		1.193 9	0.073 3
EMD方法	0.854 9	0.374 7		7.344 2	0.040 0
本文方法	0.955 9	0.373 9		9.606 9	0.025 7

下载: 导出CSV 

| 显示表格

4.   本文方法在地磁观测数据去噪中的应用

近年来，随着城市的快速扩张，特别是高铁、地铁、高压输电线路、高速公路的扩建，各地的地磁观测数据无可避免地经常受到各种类型的外加干扰，导致记录的数据量暴增，而有效信号的比率却不断下降，这对震前地磁异常特征提取的研究造成了较为严重的干扰。因此，地磁观测数据的去噪研究具有非常重要的意义。

由于定点地磁观测与地震观测都具有长期性、连续性的特征，因此，本文的方法同样可以应用于地磁观测数据的噪声去除。图6给出了邕宁地磁台2021年1月1日00时至2021年1月5日24时Z分向的地磁原始观测数据（分钟值）。可以看出，由于受到地铁干扰，波形中出现了大量高频噪声。

图  6  邕宁地磁台2021年1月1日0时至5日24时观测数据波形

Figure  6.  Observation waveform at Yongning geomagnetic station from 00:00 on January 1 to 24:00 on January 5，2021

下载: 全尺寸图片 幻灯片

图7和表4给出了使用本文方法对该段数据进行分解的结果以及各IMF分量的Hurst指数。可以看出，地铁的噪声信号主要集中在IMF1−IMF5分量，而IMF6之后的分量是地磁每日周期性变化数据。

图  7  地磁数据CEEMDAN分解结果

Figure  7.  Decomposition results of geomagnetic data by CEEMDAN

下载: 全尺寸图片 幻灯片

表  4  地磁数据IMF分量的Hurst指数

Table  4.  The Hurst exponent of different IMFs from geomagnetic data

序号	Hurst指数		序号	Hurst指数		序号	Hurst指数
IMF1	0.284 6		IMF5	0.497 9		IMF9	0.954 7
IMF2	0.323 2		IMF6	0.583 0		IMF10	0.989 9
IMF3	0.208 0		IMF7	0.702 9		IMF11	0.997 1
IMF4	0.352 4		IMF8	0.861 3		IMF12	1.006 3

下载: 导出CSV 

| 显示表格

地磁观测与地震观测不同，地磁观测主要是记录磁场的长期变化规律，而短时的干扰则是需要去除的噪声信号。因此，使用Hurst指数作为判据时，类型选择与地震数据相反，h＞0.5的分量是有效信号，而h≤0.5的分量则作为噪声信号。图8给出了去噪后的地磁波形和噪声波形，通过与原始波形的对比可以看出，地磁数据波形的周期性变化得到较完好的保留。由此可以得出，本文的方法在应用于地磁数据去噪时，可以有效地去除噪声数据，并能够较好地保留波形的变化特征。

图  8  地磁数据去噪结果

Figure  8.  Denoising results of geomagnetic data

下载: 全尺寸图片 幻灯片

5.   讨论与结论

本文提出了一种基于CEEMDAN算法和hurst指数的地震信号去噪方法，该方法已在不同类型地震波形上进行了测试。相较于传统方法，本文方法对于低信噪比波形和高信噪比波形的去噪性能分别提高了33%和6倍。同时，本文方法能够根据信号数据自身频率和振幅变化的特点，自适应地进行噪声信号识别。这意味着在处理大量地震数据时，不再需要逐类型事件分析噪声特征，这将有效减少分析人员处理数据的工作量。而对于地磁信号的滤波结果表明，本文的方法同样能够较为完整地将地铁噪声从地磁数据中识别滤除。在未来的工作中，我们将收集更多类型的地震数据和地磁数据进行验证，同时还将把该方法推广应用于连续重力和形变等观测数据的噪声处理研究中。