徐文海，丁 勇

（四川省地震局，四川 成都 611730）

在数字化观测系统中，系统的瞬态变化、仪器毛刺、数据记录的阶跃和尖峰以及由系统故障引起的信号失真和重大环境干扰源的出现等，都是影响数据记录质量的重要因素。除了真实的地震信号外，还包含了由于系统故障引起的各种失真信号以及台站周围的环境噪声干扰等信息[1]。特别是地震烈度速报与预警工程的大背景下，如何应对地震台站已存在及未来可能存在的各种类型的干扰，是我们需要长期面对的问题。这也就要求地震台站或中心站的监测人员需要依托自身熟悉自己台站或辖区内台站观测环境的优势，及时对可疑的干扰进行分析研究。本文认为在台站干扰分析排查这项工作中，台站监测人员应该尽可能地将工作做在前面，才能有效避免潜在的负面影响。为了区分以及消除干扰，首先就需要对地震波进行分析。为此，本文采用了傅里叶变换、希尔伯特-黄变换进行频谱分析，再结合功率谱对地震波干扰信息进行区分，根据干扰的频率区间设计相应的滤波器来消除干扰。

1 成都地震监测中心站职责

成都地震监测中心站的主责主业分为3部分：①围绕地震危险区地震活动态势及未来地震形势的严峻性，不断提升区域地震监测预报综合能力。将日常监测预报工作与重点监视防御区的强化短临跟踪工作相结合，开展地震震情监视跟踪和年度地震趋势会商工作。②制定针对性监测业务拓展方案。以专业台与市州局共促区域性地震监测事业发展为基点，与市县防震减灾部门协同合作，共建共享，利用成都地震监测中心站技术力量优势和地方良好的监测基础资源，优化区域监测台网布局，提升区域地震监测能力。③推动地震预警项目开展。组织人员开展预警项目建设实施指导服务和验收工作，目前已完成辖区60余个站点的验收工作。待地震速报与预警项目建设全面完成后，成都站将承担所辖片区内284个预警台共计305台预警仪器的运行维护工作。

2 测震观测系统

成都站是龙门山地震带和四川中北部地区唯一的多学科综合型观测台站。近几年，随着“中国数字地震观测网络建设项目”“东半球空间环境地基综合监测子午链项目”“中国地震局背景场项目”“汶川地震灾后恢复重建项目”“芦山地震灾后恢复重建项目”以及台站优化改造、标准化建设项目的相继完成，成都站已建成一个基础设施完善、观测设备现代化、科研技术力量强劲的综合台站。成都站测震数据以及地磁数据均参与国际资料交换，服务于全球地震定位及地磁场空间变化研究。成都和松潘测震观测系统运行稳定，观测资料记录连续，震相分析准确，后续震相丰富，其观测资料精度和分析水平在国内始终处于前列，为中国地震台网以及地震数据库提供了大量、可靠的波形资料和震相分析资料，为地震监测、科学研究及防震减灾事业作出了应有的贡献。成都台作为离5·12汶川地震最近的国家级台站，拥有一套周期360 s的超宽频带地震计，所记录的宝贵地震资料被广泛地应用于国内外对地震的科学研究。

3 方法的MATLAB实现

第一步，读取指定文件名的地震波形数据，确定参数，包括采样频率、地震计周期、数采转换因子。

第二步，调用函数y=fft（x）实现快速傅里叶变换。

第三步，调用函数z=emd（x）对数据进行EMD分解，求得边际谱和每个IMF分量及最后一个剩余分量residual的图形。

第四步，用xt=x×C/2 000将count值转化为速度值，其中C为数采转化因子。然后调用mean()函数，用公式xt=xt-mean（xt）去均值。

第五步，调用pwelch()函数实现Welch算法，即：

式（1）中：pxx为功率谱估计值；f为得到的频率点；x为进行功率谱估计的有限长序列，即上述已去均值的xt；window指所用窗函数；noverlap指重叠的点数；Nfft指采用的FFT算法长度；Fs为采样频率。

第六步，用函数sys=tf（[100]，[12×0.707×pi/T/2pi×pi/（T×T）/4]）求出频率响应，再由函数[m p]=bode（sys，2×pi×f）求出幅值m，经过矩阵转置后，由公式pxx=pxx/abs（mag^2）扣除仪器响应，再用绘出速度功率谱。其中T为地震计周期。

第七步，用db4小波对原始信号进行5层分解并提取小波系数：[c，l]=wavedec（x，5，'db4'），再用appcoef()和detcoef()求得低频与高频部分，随后进行小波消噪，之后再进行高通滤波，本文所选择的是IIR高通滤波器。

4 应用实例

本文截取了2016-07-15T02：00记录于成都地震监测中心站JCZ-1T的垂直向波形数据，时间长度为300 s，包含了一个地震事件。

4.1 分析部分

由于原始波形受到严重的干扰，地震信息被干扰掩盖难以分辨，如图1所示。快速傅里叶变换（FFT）分析后的结果如图2所示，由图2可知，波形数据中存在低频、高频的干扰，2～4 Hz的频率区间幅度较大。在Hilbert边际谱上，同样存在2～4 Hz频率区间的较大幅度反应。这是由于2～4 Hz频率范围的干扰已成为成都地震台（现成都地震监测中心站）测震资料中的固有干扰[2]。

为了验证分析的结果，设计了FIR带阻滤波器对2～4 Hz的频率区间进行滤波，毛刺干扰被压制，如图3所示，证实了2～4 Hz频率区间的较大幅度反应是波形中存在毛刺的原因。

4.2 处理部分

db4小波分解的结果如图4所示，由图4可知，低频信号振幅与原信号更接近，而高频信号比较集中且振幅较小，与原波形相差较大，去噪便是对高频信号进行处理。将去噪后进行IIR高通滤波，得到最终处理结果，如图5所示。由处理结果可知，整个处理过程中原始波形中的干扰被有效地压制。

图1 原始波形

图2 FFT幅频图与Hilbert边际谱

图3 消除2～4 Hz频率区间后的波形

图4 db4小波分解后各层低频、高频信号

图5 小波去噪及IIR高通滤波后的波形

5 认识与讨论

本文所提出的测震资料分析处理方法在地震台站的实际应用中，取得了良好的效果，对干扰进行了有效的压制。文中详细地介绍了测震数据分析处理方法的设计原理、工作流程，核心是干扰类别的识别。测震工作人员不仅局限于地震震相的划分，还需要对数据中存在的干扰类型进行判断并区分是固定干扰还是随机干扰。

本文所提出的分析处理方法的初衷是加强对干扰类别的判断与识别，以傅里叶变换、希尔伯特-黄变换以及速度功率谱作为分析判别的理论依据，小波变换和滤波器作为分析结果的验证手段。未来将对方法进行软件化，并希望在更多地震台站进行试用。