2018, Vol. 39
2008年上海测震台网开始部署自动速报系统，系统运行期间触发事件大部分为浙江舟山以及上海长江口附近的爆破事件。为提高自动速报系统信息发布的准确性，将上述区域地震信息发布震级阈值提升为ML 2.5，可基本减少爆破信息发布的可能，但导致该区域震级较小天然地震的漏报。因此，需要研究适用于上海区域的地震和爆破识别算法，并应用于实际工作。
20世纪50年代开始，国内外对天然地震与人工爆破的识别进行了广泛和深入的研究，如：刘希强等(2003, 2005)利用小波变换方法识别天然地震与爆破或塌方或矿震事件；边银菊（2002）将遗传算法（GA）和反向传播算法（BP算法）相结合成为GA-BP算法，建立遗传BP神经网络应用于地震和爆破识别；沈萍等（1999）采用非稳态的WD理论，利用地震核爆和天然地震的瞬态谱特征，瞬时频率随时间的变化等特征来识别天然地震与地下核爆炸；杨选辉等（2005）利用小波包分量比方法识别地震和爆破识别；和雪松等（2005）用小波包识别地震和矿震。卢世军（2009）归纳了基于地震波辐射图形、地震波谱分析和震相特征的地震和爆破识别判据，主要包括：波初动、震源深度、倒谱、体波震级mb与面波震级MS之比、P波初动振幅与P波位移最大振幅比、勒夫波和瑞利波振幅比、P波与S波谱振幅比、小波变换、P波谱振幅与勒夫波谱振幅比、瞬态谱等。
本文对上海测震台网记录的上海周边区域天然构造地震和爆破事件波形进行小波包分解，并提取特征向量，提出用支持向量机（SVM）识别天然构造地震和人工爆破的方法。结果表明，基于SVM算法的向量识别分类方法在区域天然地震和爆破事件识别研究中是可用的，识别准确率预计达85%以上。
1 资料收集选取2008年11月至2016年9月上海市地震监测中心记录的34个浙江舟山爆破事件和66个天然地震事件为研究样本，事件目录详见表 1和表 2。
|
表 1 爆破事件目录 Tab.1 Catalogue of explosion events |
|
|
表 2 地震事件目录 Tab.2 Catalogue of earthquake events |
选取事件原始记录波形为SEED格式或EVT格式，为方便在Mablab中进行研究和计算，对选取的数据文件进行预处理，即将SEED/EVT文件按通道导出为ASCII格式文件。ASCII格式文件是地震台站实际记录的地动数据，由纯数字构成。在Matlab中使用importdata函数，将波形记录导入Matlab，使用save函数，将导入的数据保存为Matlab可直接分析处理的mat文件（王军，2017）。为研究方便，选取垂直向信号进行分析。
因所有数据均来源于上海测震台网，地震台站配置相同，可以自行编制基于Matlab软件环境的批处理软件，对待处理地震事件进行预处理。示例代码如下
importdata SH_DAX_SHZ_2.txt %读入SH_DAX_SHZ_2.txt
DAX=ans %将读入值赋值给数组DAX
save 'DAX.mat' DAX %将工作区的DAX数组保存为DAX.mat
2.2 细节分量提取为得到高信噪比的震相记录，将原始垂直向记录进行2层小波分解并提取细节分量作为研究对象，预处理前后波形记录对比见图 1。
|
图 1 滤波处理前后波形记录对比 Fig.1 Comparison of waveform before and after filtering |
处理过程在Matlab软件中完成，示例代码如下
[c, l]=wavedec(DAX, 2, 'haar'); %将原始垂直向信号进行2层Haar滤波
DAX0=detcoef(c, l, 1); %提取第一层细节分量
2.3 波形检查对于每一个事件，不可能所有地震台站均能记录到清晰震相，在实际处理过程中，需对记录进行检查，剔除事件发生时的故障和无记录地震台。为方便波形检查工作，基于Matlab软件编写波形检查代码，画出滤波后所有地震台站记录波形，查找并剔除无数据和事件记录不正常的地震台。13个地震台波形检查结果见图 2。
|
图 2 波形检查 Fig.2 Check the waveform |
小波包变换能够为信号提供一种更加精细的分析方法。该变换将频带进行多层次划分，对多分辨率分析未细分的高频部分进一步分解，并能够根据被分析信号的特征，自适应的选择相应频带，使之与信号频谱相匹配，从而提高时频分辨率，因此小波包具有广泛的应用价值。
对于小波包分析的理解，此处以一个3层分解进行说明。小波包分析树见图 3，图中：A表示低频分量，D表示高频分量，末尾的序号表示小波包分解的层数。分解具有以下关系：S=AAA3+DAA3+ADA3+DDA3+AAD3+DAD3+ADD3+DDD3。
|
图 3 小波包分析树 Fig.3 Decomposition tree of wavelet packet |
将记录进行4层db1分解基的小波包分解，每层分别重构，重构系数利用最大奇异值分解算法（SVD）计算得到原始特征向量，将原始特征向量进行归一化处理，即得到归一化特征向量。将每个事件有记录台站的归一化特征向量取平均值，得到该事件的特征向量。以某事件有HUH和TMS 2个地震台清晰记录为例，基于Matlab环境编写程序代码，部分代码如下
HUHt=wpdec(HUH0, 4, 'db1', 'shannon'); %4层小波包分解预处理后提取的细节向量
TMSt=wpdec(TMS0, 4, 'db1', 'shannon');
For n=1:16
HUHcg(n, :)=wprcoef(HUHt, [4, n-1]); %按层实现第n层重构
TMScg(n, :)=wprcoef(TMSt, [4, n-1]);
HUHnorm(n)=norm(HUHcg(n, :)); %第n层重构后得到的原始特征值
TMSnorm(n)=norm(TMScg(n, :));
end
HUHnormsum=sum(HUHnorm);
TMSnormsum=sum(TMSnorm);
for n=1:16
HUHnorma(n)=HUHnorm(n)/HUHnormsum; %特征值归一化处理
TMSnorma(n)=TMSnorm(n)/TMSnormsum;
R02(n)=(HUHnorma(n)+TMSnorma(n))/2;%求取平均特征向量
end
将66个天然地震事件的特征向量合成一个66×16的矩阵，将34个爆破事件特征向量合成一个34×16的矩阵。实际小波包分解后各结点频带不按照从小到大顺序排列，频谱最终形成格雷码顺序。编写基于Matlab的软件代码，调整100个事件特征向量各元素之间的顺序，得到爆破特征向量库和地震特征向量库，具体数值见图 4、图 5。
|
图 4 爆破事件特征向量库 Fig.4 Eigenvector library of explosion events |
|
图 5 地震事件特征向量库 Fig.5 Eigenvector library of earthquake events |
研究过程中采用基于支持向量机（SVM）的向量分类识别算法。SVM是一种常见识别方法，在机器学习领域通常用于模式识别、分类以及回归分析。如图 6所示，二维平面上有正常元素和异常元素两种数据元素，分别用五角星和圆圈表示，识别算法就是找到一条直线将2类数据分开，一边对应的全是正常元素，另一边所对应的全是异常元素。
|
图 6 识别算法基本原理 Fig.6 Basic principle of recognition algorithm |
SVM的主要思想可以概括为：①基于结构风险最小化理论，在特征空间中建构最优分割超平面，使得学习器得到全局最优化，并在整个样本空间的期望风险以某个概率满足一定上界；②对于线性不可分情况，通过使用非线性映射算法，将低维输入空间线性不可分样本转化为高维特征空间，使其线性可分，从而使得高维特征空间采用线性算法对样本的非线性特征进行线性分析成为可能，见图 7。
|
图 7 SVM算法高维思想 Fig.7 Higher dimensions in SVM |
SVM识别算法基于Matlab软件环境完成。首先生成学习样本库和样本类型向量。部分代码如下
for n=1:66
all(n, :)=T(n, :); %学习样本库all的1-66行为天然地震特征向量
class(n)=0 %类型向量class前1-66行赋值0，代表天然事件
end
for n=67:100
all(n, :)=B(n-66, :); %学习样本库all的67-100行为非天然事件特征向量
class(n)=1 %类型向量class的67-100行赋值1，代表非天然事件
end
采用svmtrain函数进行支持向量机的自我学习，核函数为线性核函数，学习代码为
svmstruct=svmtrain(all，class)
当得到一个未知类型事件特征向量unknow时，可以使用vmclassify函数，结合预测识别方案，对事件类型进行自动识别。识别代码如下
classunknow=svmclassify(svmstruct，unknow) %unknow为未知事件的特征向量
4.3 自验证使用原始样本库中100个事件的特征向量，对预测识别方案进行自验证，自验证代码如下
for n=1:66
group(n)=svmclassify(svmstruct, T(n, :))-class(n);
end
for n=67:100
group(n)=svmclassify(svmstruct, B(n-66, :))-class(n);
end
plot(group)
title('预测结果')
自验证结果显示，66个天然地震中有7个被识别为爆破，34个爆破事件中有4个被识别为天然地震，可见100个事件的自验证结果准确率为89%。自验证结果证明，基于SVM算法的向量识别分类方法在区域天然地震和爆破的识别研究中是可用的。
5 结束语使用基于SVM算法的向量识别分类方法，对研究样本库事件自验证识别的准确率为89%，在后期使用中，因不断有新资料的补充，实际预测识别准确率预计将达85%以上。支持向量机的算法具有学习功能，新的事件发生后，可以以预测识别方案对事件类型进行识别；当事件类型明确后，可将特征向量加入学习样本库，随着样本库的不断增加，预测识别方案也不断得到修正，预测识别的准确率随之得到提高。
在后期工作中，可将识别算法接入上海测震台网自动速报系统，实现地震数据预处理、小波包分解、特征向量提取、事件类型识别等工作的自动化完成，进一步提高工作效率。
| 边银菊. 遗传BP网络在地震和爆破识别中的应用[J]. 地震学报, 2002, 24(5): 516-524. | |
| 和雪松, 李世愚, 沈萍, 等. 用小波包识别地震和矿震[J]. 中国地震, 2006, 22(4): 425-434. | |
| 刘希强, 杜贻合, 徐波, 等. 区域矿震与地震的定量识别方法及其应用[J]. 中国地震, 2005, 21(1): 50-60. | |
| 刘希强, 沈萍, 张玲, 等. 用小波变换能量线性度方法识别天然地震与爆破或塌方[J]. 西北地震学报, 2003, 25(3): 204-209. | |
| 卢世军. 天然地震与人工爆破波形特征提取与识别算法研究[D]. 广西师范大学, 2009 http://www.wanfangdata.com.cn/details/detail.do?_type=degree&id=D069900 | |
| 沈萍, 郑治真. 瞬态谱在地震与核爆识别中的应用[J]. 地球物理学报, 1999, 42(2): 233-240. | |
| 王军. Mallat算法在数字地震信号压缩中的应用[J]. 地震地磁观测与研究, 2017, 38(5): 133-138. | |
| 杨选辉, 沈萍, 刘希强, 等. 地震与核爆识别的小波包分量比方法[J]. 地球物理学报, 2005, 48(1): 148-156. |