引言

我国是世界上岩溶最为发育的国家之一，岩溶面积达363万平方千米，主要分布在黔、桂、滇东、川南、湘西等省区，其中分布最广的是贵州省，约占全省面积的73%，其次是广西壮族自治区，约占全区面积的34% （覃子建等，1996）。近年来随着我国经济的飞速发展，工农业、居民生活用水迅速增加，大量地下水被抽取，导致地下水位不断下降，塌陷事件也逐年增多，岩溶塌陷已成为我国六大地质灾害之一（赵博超等，2015）。然而，天然地震与岩溶塌陷的波形具有较高的相似性，分析人员依靠经验根据波形特征进行二者的识别费时较多，缺乏时效性。因此亟需研究一种能够及时准确地识别天然地震与岩溶塌陷事件的方法，建立准确的事件目录，为快速开展塌陷灾害防治工作提供重要的科技支撑，为地方经济发展、国土资源规划提供理论基础。

近年来国内外研究人员针对地震事件类型的识别开展了广泛而深入的研究，总结了许多时域和频域方面的识别判据。时域的判据主要包括初动方向、震源深度、体波震级与面波震级比、P波初动振幅与P波和S波最大振幅比、P波与S波最大振幅比、P波和S波最大振幅与尾波持续时间比、体波和面波震相特征、P波视速度差异等（吴宝峰，闫峰，2009；Pasyanos，2010；Taylor，2011；刘莎等，2012；王婷婷，边银菊，2015）；而频域的判据主要有能量频谱分布、拐角频率、倒谱曲线、瞬态谱等（张萍等，2005；Yılmaz et al，2013 ；靳玉贞等，2015；Kortström et al，2016 ；叶庆东等，2017）。这些方法各有特色，均取得了很好的效果，对于建立准确的地震事件目录具有重要意义；但是这些方法多是单独使用时域特征或频域特征来进行识别，而且地质条件、仪器参数、波形传播路径等差异使得其普适性不足。因此，许多研究人员使用在时域和频域均具有很高分辨率的小波包变换方法来进行特征判据提取。和雪松等（2006）通过对地震和矿震波形进行小波包变换，再计算二者的奇异值，从而达到识别地震的目的；王婷婷等（2014）使用小波包对地震和爆破波形进行分解，根据P波和S波的发育频段不同，使用不同分解系数中的P波和S波能量进行对比，获得了地震和爆破的识别判据；Beccar-Varela等（2016）通过小波变换对美国亚利桑那州M3.0—3.5地震和爆破事件进行分解，利用分解得出的时频能量特征进行地震事件的识别。尽管上述利用小波包分解提取判据的方法能够有效地进行地震识别，但是这些方法主要是针对地震与爆破或矿震，至于地震与岩溶塌陷事件识别的研究则较少。此外，上述方法尽管提取了地震波形在时频域的特征，但均未涉及地震波形分形维数的特征。

为此，本文拟根据天然地震和岩溶塌陷在不同小波包分解频带内存在的差异，结合分形维数计算方法，对比各小波包分解系数的维数差异，以期提出一种能够体现出事件波形时频域特征和分形维数特征并且能够快速有效地识别天然地震与岩溶塌陷的新方法。

1.   资料

本文从广西地震台网中心2008—2015年记录到的塌陷和地震事件中，选取记录台站较多、波形较为清晰的广西境内的39个事件，包括17次塌陷事件和22次浅源地震事件，震源深度为5—10 km，震级范围为M_L1.0—2.9，所有事件的类型均已由广西地震台网中心确认。从这39个事件中，提取出400条记录清晰的事件波形记录，其中天然地震事件和岩溶塌陷事件的波形各200条，最大、最小震中距分别为237 km和25 km。所选地震、塌陷事件和记录台站的分布如图1所示。

本文所选用事件的波形记录的采样率均为100 sps，波形记录的频率范围为0—50 Hz；所有事件只使用垂直向波形进行计算，事件波形的截取宽度均为120 s （即12 000个点），S波的初动位置尽量位于截取窗口的中部。

图  1  事件和地震台站位置示意图

Figure  1.  Location of earthquakes （green circles），karst collapses （red circles） and seismic stations （triangles） in this study

下载: 全尺寸图片 幻灯片

2.   识别算法设计及基本原理

2.1   识别算法设计

一般来说，一套完整的自动识别算法主要由两个模块组成：其一为特征提取模块，用于从目标数据中提取出识别所需的特征参数；其二为识别判定模块，利用已提取出的特征参数对目标类型进行判定。本文在设计两个模块时考虑到：① 本文所使用的特征提取算法应能体现出天然地震和岩溶塌陷的震源机制在时频域的特点，同时还要能描述出两种事件在不同频带下波形的复杂性；② 由于测震仪器、地质构造、传播路径等差异，每一组事件所提取的特征数据可能存在一定差异，因此识别算法需要具有一定的鲁棒性和容错性，同时从现实应用的角度出发，算法还需具有较快的运算速度。

鉴于此，本文使用小波包分形方法来提取天然地震和岩溶塌陷的特征，再通过径向基函数(radial basis function，简写为RBF)神经网络进行事件分类的识别判定。小波包分形方法利用了小波包在时频域下高分辨的特点，并通过分形理论对不同频带下波形信号的复杂度进行描述，从而达到兼顾地震信号的细节分析与整体描述合理性的目的（曹茂森等，2005；谢全民等，2011；朱权洁等，2012）。而RBF神经网络除了神经网络普遍具有的自适应性强、能够实现各种非线性映射、有较高的鲁棒性和容错性、运算速度快等特点外，还具备结构简单、收敛快等优势。所以本文使用小波包分形作为特征提取模块，RBF神经网络作为识别判定模块。首先，对事件波形进行归一化处理，然后利用小波包分形提取出事件波形特征；其次，利用已知事件组成训练集完成对RBF神经网络的训练，并将训练好的神经网络用于对未知事件的类型进行识别。图2为本文算法的流程框图。

图  2  本文算法流程图

Figure  2.  The flow chart of algorithm in this study

下载: 全尺寸图片 幻灯片

2.2   小波包变换

小波变换是一种重要的时频域波形分析方法，在数字信号处理、图像处理、声纹识别等方面均有广泛应用。而小波包变换是对小波变换的进一步扩展，弥补了小波变换对高频细节特征提取的不足。小波包变换方法是通过多组低通、高通正交滤波器将信号的各频段不断细分，以此来提高信号在整个频带内的分辨能力，进而获得信号中所包含的更多细节。也正因为此，小波包变换也被称作信号时频域的显微镜。

2.3   分形维数

分形维数是分形理论中用于描述曲线或者图形不规则程度的重要特征量，但是，在实际使用中往往难以直接计算得到不规则分形体的维数，而只能计算其近似值。因此，本文从快速计算的角度出发，选用计盒维数来计算维数。

计盒维数D_B（F）定义为（张济忠，1995）：对于需要计算维数的目标图形F，构造足够多的边长均为δ的正方形盒子将其完全覆盖（图3），所用到的盒子数量为N_δ（F），那么F的计盒维数D_B（F）则为当δ趋于0时N_δ（F）增加的对数速率，可以表示为

图  3  波形分形盒维数示意图

Figure  3.  Schematic diagram of fractal box dimension

下载: 全尺寸图片 幻灯片

3.   特征提取

由于地震与岩溶塌陷的震源性质和传播路径有所不同，台站记录到两种事件的波形在时频谱中呈现不同的特征。因此，通过时频域分析来进行地震和塌陷事件的特征提取具有可行性。图4为归一化处理后的地震和塌陷事件的时频能量分布图，所示的两事件是由贞丰地震台记录到的地震事件（2010-01-25，M_L2.9，震中距Δ=128 km）和塌陷事件（2015-09-09，M_L1.8，Δ=124 km），两个事件的震中均位于河池市天峨县。可以看出，地震的优势频率处于0—20 Hz范围内，而塌陷的优势频率处于0—10 Hz范围内。

图  4  地震（a）、岩溶塌陷（b）事件的波形（上）及其时频谱（下）

Figure  4.  Waveforms （upper） and spectrograms （lower） of the earthquake （a） and the karst collapse （b）

下载: 全尺寸图片 幻灯片

3.1   事件波形的小波包分解

使用小波包进行变换前首先需要选取小波包的基函数。常用的基函数有Haar，Daubechies，Mexican Hat，Morlet，Meyer，Symlet，Coiflet，双正交小波等，不同的基函数性质不同。考虑到基函数需具备快速计算的能力，以便增强算法的时效性，同时需要具备精确重构的性质，以便能够完整地反映原始波形的频率特性，为此，本文选取具有紧支撑正交性的Daubechies作为基函数。考虑到需要避免因支撑区间太大而产生边界问题或者因区间太小而使得能量无法集中，选择区间长度为7的4阶Daubechies小波，简写为db4。

对本文所选的400条事件波形记录进行db4小波包的4层分解，得到16个分解系数波形，对应的频带宽度分别为0—3.125 Hz，3.125—6.250 Hz，…，46.875—50 Hz。图5和图6分别给出了河池市天峨县的地震和塌陷事件的小波包分解结果。可以看出，虽然地震与塌陷事件的优势频率不同，但是在各频段内均有震动波形的记录。因此，为了获取事件在不同频段内的波形特征，提取特征时需对各频段的数据进行计算。

图  5  2010年1月25日M_L2.9地震事件的第四层小波包分解

Figure  5.  The fourth-layer wavelet packet decomposition of the M_L2.9 earthquake on 25 January 2010

下载: 全尺寸图片 幻灯片

图  6  2015年9月9日M_L1.8岩溶塌陷事件的第四层小波包分解

Figure  6.  The fourth-layer wavelet packet decomposition of the M_L1.8 karst collapse on 9 September 2015

下载: 全尺寸图片 幻灯片

3.2   小波包系数的分形维数计算

本文采用计盒维数法计算分形维数（张济忠，1995），分别计算各个事件的小波包分解系数，最终获得由16个维数值组成的一维向量，此即为本文所要提取的地震和岩溶塌陷事件的特征。下面将以这16个维数值为输入向量，进行RBF人工神经网络训练和事件类型识别。图7给出了地震和岩溶塌陷事件小波包分解系数的分形维数计算结果，可以看出，岩溶坍塌事件在不同频段内波形的维数高于地震事件波形的维数。

图  7  地震与岩溶坍塌事件的小波包分形维数计算结果对比

Figure  7.  Comparison of wavelet packet fractal dimension of earthquake with that of collapse

下载: 全尺寸图片 幻灯片

4.   RBF神经网络

RBF神经网络是一种经典的前馈3层神经网络，根据人脑皮层神经网络结构提出，具备以任意精度逼近任意一个连续函数的能力，由输入层、隐含层、输出层组成，其中输入层将输入数据与网络连接起来，隐含层对输入层传递过来的数据进行非线性变换，输出层则负责对隐含层进行加权后输出（Haykin，2001）。图8给出了单输出的RBF神经网络结构图，由图可知整个网络的输出y_n可以表示为

图  8  RBF神经网络结构图

Figure  8.  Structure of RBF neural networks

下载: 全尺寸图片 幻灯片

式中，n为样本输入数据的数量，m为隐含层节点数，g_i（x）为神经网络的径向基函数。

本文所使用的RBF神经网络结构为：① 输入层：神经网络的输入为小波包分形得出的16个特征值，所以输入层节点数为16，节点的权值均为1；② 隐含层：径向基函数选择应用广泛的高斯函数，其表达式为

This page contains the following errors:

Below is a rendering of the page up to the first error.

RBF神经网络是否优秀，关键在于隐含层的训练，即基函数中心的选取、节点数量的选择以及节点权值的训练。本文使用matlab神经网络工具箱来进行隐含层的训练，目标均方差设为0.01，高斯函数扩展步长为2.0。训练的具体步骤为：① 建立一个0神经元的网络；② 寻找最大方差样本；③ 增加一个神经元，并使用最大方差的输入样本作为基函数中心；④ 重新设计隐含层权值，使用能够实现最小输出误差的值作为隐含层权值；⑤ 检查输出的均方差是否小于所设定的目标值，若满足该条件则网络训练完毕，否则返回第二步，从剩余样本中继续寻找最大方差样本。

5.   识别能力检验

神经网络的训练和测试由训练集和测试集完成。本文按照发震时间的先后顺序，将前10次地震事件的100条波形数据和前9个岩溶塌陷事件的100条波形数据组成训练集，用于神经网络训练；测试集则是由余下12次地震事件的100条波形数据和8次岩溶塌陷事件的100条波形数据组成，用于验证分类能力。之后分别计算测试集和训练集的地震和岩溶塌陷事件的小波包分形维数，计算结果列于表1。

表  1  地震和岩溶塌陷事件的小波包分形维数均值和方差

Table  1.  Mean and variance of fractal dimension of earthquake and collapse events

小波包 分解系数	训练集（地震）			训练集（塌陷）			测试集（地震）			测试集（塌陷）
	维数均值	维数方差		维数均值	维数方差		维数均值	维数方差		维数均值	维数方差
1	1.379 6	0.010 0		1.463 0	0.009 2		1.389 0	0.010 4		1.440 8	0.008 5
2	1.341 8	0.008 2		1.412 2	0.010 0		1.337 6	0.007 7		1.438 0	0.006 5
3	1.354 4	0.013 3		1.435 1	0.014 0		1.358 0	0.010 1		1.466 3	0.010 2
4	1.344 9	0.009 7		1.400 4	0.013 2		1.333 9	0.008 1		1.444 3	0.009 1
5	1.416 3	0.017 9		1.483 0	0.018 5		1.426 5	0.017 2		1.505 5	0.011 0
6	1.388 8	0.016 0		1.438 4	0.016 0		1.395 4	0.012 8		1.478 4	0.010 2
7	1.368 3	0.014 0		1.450 2	0.019 9		1.364 6	0.012 0		1.471 0	0.008 8
8	1.370 0	0.014 5		1.451 4	0.016 6		1.370 1	0.012 6		1.462 8	0.010 3
9	1.411 5	0.011 8		1.539 9	0.015 9		1.414 9	0.013 4		1.557 0	0.009 4
10	1.438 5	0.017 3		1.508 4	0.010 2		1.443 9	0.015 7		1.524 9	0.011 6
11	1.423 4	0.017 0		1.493 3	0.010 0		1.433 3	0.015 5		1.518 4	0.007 7
12	1.429 8	0.017 4		1.486 3	0.010 5		1.440 5	0.016 2		1.513 1	0.009 6
13	1.422 9	0.019 2		1.530 1	0.009 8		1.438 5	0.017 2		1.539 2	0.007 7
14	1.422 4	0.018 8		1.520 0	0.009 5		1.433 2	0.017 8		1.536 9	0.006 0
15	1.421 1	0.018 4		1.497 6	0.010 0		1.425 9	0.016 6		1.516 1	0.008 6
16	1.420 3	0.017 6		1.501 0	0.016 2		1.429 7	0.017 4		1.524 8	0.007 1

下载: 导出CSV 

| 显示表格

从表1可以看出，岩溶塌陷事件各小波包分解系数计算出的分形维数值均明显大于地震事件对应分解系数的维数值。使用训练集中200条波形小波包分形特征向量进行RBF神经网络的训练，训练完成后，分别使用训练集和测试集的波形进行分类检验。训练集和测试集的识别结果如表2所示，可以看出：对于训练集，本文的识别方法可以达到100%的识别率；而对于测试集，地震和岩溶塌陷的识别率可分别达到87%和92%，能够有效地识别出天然地震与岩溶塌陷事件。由表2还可以看出，本文错误识别的共有9个地震事件13条波形记录和5个岩溶塌陷事件的8条波形，其相关信息列于表3。通过分析错误识别的事件波形，本文认为导致错误识别的原因在于：地脉动等背景噪声幅值较大，导致事件波形信号不明显；事件波形存在零漂现象；事件波形中叠加了少量干扰，出现了多个波形信号；事件波形的毛刺多，信号不够平滑；事件波形中有断记，部分波形为直线。

表  2  小波包分形维数判据的识别结果

Table  2.  The recognition result of the criteria based on wavelet packet fractal

事件类型	测试样本数	正确识别样本数	错误识别样本数	识别率
地震（训练集）	100	100	0	100%
塌陷（训练集）	100	100	0	100%
地震（测试集）	100	87	13	87%
塌陷（测试集）	100	92	8	92%

下载: 导出CSV 

| 显示表格

表  3  错误识别事件的相关信息

Table  3.  The directory of false recognition

序号	发震时间 年-月-日	事件类型	ML	地点	东经/°	北纬/°	台站	震中距/km
1	2010-08-26	地震	2.2	钦州	22.24	109.47	LNS	26
2	2010-09-29	地震	2.7	河池	25.24	106.99	TE	36
3	2010-09-29	地震	2.7	河池	25.24	106.99	YIZ	188
4	2011-02-28	地震	2.9	百色	24.42	108.11	YAX	101
5	2011-02-28	地震	2.9	百色	24.42	108.11	PGX	122
6	2011-02-28	地震	2.9	百色	24.42	108.11	GAL	136
7	2011-10-03	地震	2.9	百色	24.61	105.40	YTT	143
8	2012-03-28	地震	2.5	来宾	23.64	108.90	PGX	136
9	2013-03-28	地震	2.2	崇左	21.88	107.12	PGX	128
10	2013-06-17	地震	2.3	崇左	23.42	107.50	SHX	115
11	2014-02-06	地震	2.5	百色	23.38	107.51	SHX	118
12	2014-02-06	地震	2.5	百色	23.38	107.51	QIP	122
13	2015-01-09	地震	2.3	钦州	22.25	108.04	LNS	112
14	2010-06-01	塌陷	1.9	防城港	23.71	109.23	LIB	34
15	2015-01-14	塌陷	1.0	百色	23.14	106.39	JXX	25
16	2015-09-08	塌陷	2.4	来宾	23.88	108.89	XCT	35
17	2015-09-08	塌陷	2.1	来宾	23.85	108.89	DHX	82
18	2015-09-08	塌陷	2.1	来宾	23.85	108.89	TAP	95
19	2015-09-08	塌陷	2.1	来宾	23.85	108.89	JL	200
20	2015-09-08	塌陷	2.1	来宾	23.88	108.89	TAP	95
21	2015-09-09	塌陷	1.6	河池	25.24	107.00	GAL	90

下载: 导出CSV 

| 显示表格

天然地震事件和岩溶塌陷事件的波形差异的本质是两种事件产生的震源机制不同。天然地震从本质上来讲地壳介质在构造应力作用下的破裂过程，通常震源深度较深，破裂过程较长，同时地震波经地下介质多次反射、折射、散射、叠加，导致地震仪所记录到的地震波频率范围更宽，尾波持续时间更久；而岩溶塌陷则发生在地表或者地下几十米处，是由于塌陷体从周边岩体上脱落并撞击底部岩体而产生的瞬时压缩源，能量释放时间短，波的传播路径较简单，未经过多层介质干扰，频带范围窄，尾波持续时间短（靳玉贞等，2015）。

为了能够得到量化的地震波形小波包分解的特征，本文使用分形维数来计算特征参数。地震波形分形维数的大小能够在一定程度上反映出岩石所受加载应力的差异。当岩石受到不同大小应力作用时，所产生的瞬态弹性波（称为声发射技术）的分形维数变化规律也验证了这一点。尹贤刚等（2009）运用分形理论方法研究了不同加载应力比阶段岩石声发射的分形维数，结果表明加载应力达到岩石峰值应力的40%左右时，分形维数即可达到最大值，之后逐步减小直至发生破裂；张黎明等（2015）对大理岩卸围压破坏过程的研究结果表明，应力达到峰值前，声发射分形维数会大幅度降低。由此可见，由于岩溶塌陷事件震动所产生的加载力并未达到岩石破裂峰值，而天然地震事件的加载力已经达到破裂峰值，所以岩溶塌陷事件波形的分形维数要大于地震事件的分形维数，这可以作为识别两种事件的判据。

6.   讨论与结论

本文利用基于小波包的分形维数和RBF神经网络识别技术，对2008—2015年广西境内的17次岩溶塌陷事件和22次浅源地震事件进行了识别研究和识别效果的检验，得到以下几点认识：

1） 小波包分形维数识别方法对于训练集的识别能力能够达到100%，对于测试集的平均识别能力达89.5%，说明该方法能够有效地识别天然地震与岩溶塌陷事件。而且，一个地震事件通常会被多个台站记录到，综合多个台站的识别结果，可提高事件识别的准确率。

2） 时频谱图和小波包分解系数对比显示，天然地震波形的优势频段为0—20 Hz，岩溶塌陷的优势频段为0—10 Hz。尽管两种事件的优势频率不同，但是二者所产生的波形在各个频段上均有记录，因此进行时频分析时不能只计算优势频率，应将全频段考虑在内。

3） 地震事件的各小波包分形维数均小于塌陷事件的小波包分形维数，这可能与岩石受到的应力大小有关，岩石受力越接近临界状态，所产生震动的分形维数越小。

在实际工作中，可将本文获得的量化指标作为快速自动化识别地震与岩溶塌陷事件的判定依据之一。鉴于地震破裂过程、传播路径等因素的复杂性，地震事件性质的判定还需结合近台记录P波初动符号、P波与S波振幅比、现场地震地质调查、地震宏观现象等资料综合判定。

不同地区的地质背景、塌陷方式、传播路径等方面差异较大，波形记录的特征也有所不同。因此，本研究获得的指标仅具有一定参考价值，不一定适用于其它地区。本文研究样本为已知类型事件，下一步将收集更多的事件样本，进行更广泛的识别训练和研究，不断地完善和验证判定指标的普适性。