0 引言

开展强震前临震微震动信号的拾取和特征机制分析，对于开展区域地震活动性和趋势判定具有重要作用。国际地震学家做了很多地震预测的研究，许多与地震有关的临震破裂前兆现象被人工分析识别出来(Obara et al，2004；Ito et al，2006)。相关研究表明，宽频带仪器地震记录中叠加于正常脉动基础上的“干扰事件”可能来源于临震震源预活动、微破裂、微震动或与地震孕育进入临震阶段相关联的活动构造微活动、微破裂、裂隙扩展等过程(杨立明等，2018a)；对汶川、玉树等历史强震开展强震临震微波动的研究结果显示，强震在震前12~14天左右距震中一定范围的台站在11~16Hz频带的噪声会有明显增加(杨立明等，2018a、2018b)。国内一些学者试图通过分析连续波形查找汶川地震震前是否存在低频扰动，通过使用快速傅里叶变化分析震前12天的连续波形，发现一些台站在0.1~0.3Hz频带内的振幅在临震前快速持续增大(杨立明等，2015)；通过对一些台站震前数十小时连续波形进行频谱分析，发现0.1~0.18Hz的低频扰动信号在地震发生前迅速增强(胡小刚等，2008)。通过分析低频事件的活动，即地震孕育间断期整个应力增长和释放的物理过程，可以提高长期预测地震活动的能力(Kanamori，2008)。

地震是地壳运动的结果，存在一个由缓慢变形向快速发展的过程(马瑾，2016)，相关研究采用的方法主要在观测数据0.1Hz以上高频带开展(郝春月等，2012)。同时，随着我国地震监测台网的不断加密，在开展临震预测研究提供了更丰富的数据的同时，海量数据的处理也对处理效率提出了更高的要求。杨立明等(2017)在兰州搭建地脉动实时跟踪技术系统，通过6台服务器处理200多个台站数据。我国已启动的“国家地震烈度速报与预警工程”将建设总数超过10000个台站的地震烈度速报与预警观测网络(蒋长胜等，2016)，如何对海量地震监测数据开展高性能计算和分析显得尤为重要。

因此，本文提出一种对宽频带地震连续观测数据进行多频点背景噪声异常变化分析的方法，并针对海量数据计算场景采用分布式大数据技术提升计算效率，采用列式存储数据库Hbase提升数据检索效率，通过计算结果频点蕴含的能量变化来对强震前是否孕育相关信号特征进行识别，并以2013年芦山7.0级地震为例进行验证分析。

1 数据处理与分析方法 1.1 噪声功率谱计算方法

台站地动噪声功率谱密度(PSD)是描述台站环境地动噪声的重要参数，简称噪声功率谱。本文以小时为单位对记录的原始数据采用噪声功率谱算法，计算连续波形数据的地脉动噪声，从50Hz周期开始至200s周期提取多频点噪声数据，且不对原始观测数据进行任何滤波处理。本文采用IRIS DMC的标准算法对测震台网的连续波形记录数据计算功率谱密度(Peterson，1993)，不进行任何事件或者干扰的去除工作，对于离散频率值f_k=k/(NΔt)，k=1，2，…，N-1，傅里叶分量定义为

$ Y_k=\frac{Y\left(f_k, T_r\right)}{\Delta t} $

(1)

其中，Δt为采样间隔，N为该时间序列段中的采样点数，即N=T_r/Δt，Y_k为快速傅里叶变换算法(FFT)给出的离散傅里叶变换(DFT)分量。总功率谱密度估计值定义为

$ P_k=\frac{2 \Delta t}{N}\left|Y_k\right|^2 $

(2)

对每一段数据，采用如下步骤计算其功率谱密度：

(1) 以1h为单位对进行预处理，从时间序列中去掉均值和趋势；

(2) 应用汉宁(Hann)窗以减少谱泄漏；

(3) 计算傅里叶变换，见式(1)；

(4) 去掉仪器响应，得到以地动物理单位表示的频谱；

(5) 应用1/8倍频程带宽高斯滤波器作用于傅里叶变换结果，按照式(2)计算PSD。

1.2 分布式计算和存储

近年来大数据技术不断发展，先后推出了MapReduce、Spark、Flink、Storm等分布式处理技术(Chintapalli et al，2016)，并在地震行业开展了一定应用(郭凯等，2017)。MapReduce在数据处理过程中需要在硬盘进行读取和写入(Rizvandi et al，2012)，而基于内存计算的Spark相对MapReduce在计算性能上体现了相当大的优势(Mavridis et al，2017)。开展强震前临震信号的特征提取和分析时，为提高处理效率，需要对海量地震连续观测数据进行高并发的计算与分析，基于PSD计算方法，建立一套具有高可靠、高处理性能、可在线弹性伸缩、可不间断接收任务的处理模型，采用弹性数据集Spark RDD将PSD计算任务自动分配到计算节点，解析存储在HDFS中的测震波形数据，计算结果采用RowKey方式放入分布式数据库HBase中(郭凯等，2021)，数据处理流程如图 1所示。

Spark主节点在接收到计算任务需要处理的参数后(日期、台站以及通道)，将任务分发到各个计算节点，由各个节点完成计算过程并将结果反馈给主节点。当处理层在接收到PSD计算请求时，算法调度模块根据请求参数中的时间范围、台网以及台站等参数查询出需要计算的每个文件的具体位置，同时将文件位置以及任务ID推送至计算模块的sever端；计算模块在接收到计算请求后，通过参数中的路径直接从HDFS上获取需要处理的miniseed文件，并对每个miniseed文件进行计算。由于计算结果数据量较大，为了便于数据存储和分析，数据结果存储格式上采用NoSQL列式存储，如表 1所示。

2 芦山7.0级地震震例分析

芦山地震震中台站分布情况如图 2所示，按照台站分布选取6个台站进行分析，其中震中50km范围内2个台站，震中50~100km范围内2个台站，震中100~150km范围内2个台站，台站信息如表 2所示。

对表 2中台站2013年3月1日—4月19日的数据，按照小时为单位进行计算背景噪声值计算，数据分析周期为0.1~150s，为减少毛刺、突变对结果的影响，对噪声数据做了适当平滑和毛刺过滤，选取0.1s、30s、100s、150s水平向BHN和垂直向BHZ通道部分计算结果进行分析(图 3、图 4)。

在进行地震背景噪声分析时发现，经常会出现单点噪声值的突然升高且出现频率呈一定周期性现象，如BAX台在3月1日、3月13日、3月14日、3月28日、4月4日和4月18日出现该现象，噪声值范围在-60dB~-70dB之间，通过分析对应时序数据，如图 5所示为BAX台对应时间点出现毛刺时的波形数据，可以看出台站可以正常记录到数据，“毛刺”为仪器标定对计算结果造成的影响，而非计算方法或者数据缺失造成的影响。

除标定影响外，还存在非周期性的一些噪声突然升高，在地震台站垂直向背景噪声时序变化图(图 4)中，1s和30s周期在4月16日出现较明显的噪声变化，以MDS台站为例在3月9日、4月16日出现噪声值范围在-110dB左右的突然升高，然后迅速回落，影响在2~3h左右，结合地震目录和波形数据，分析为天然地震产生的影响，分别为2013年3月9日四川汶川4.5级和2013年4月16日西藏墨脱4.8级地震。

通过对多频点计算结果的分析，在芦山7.0级地震震中附近发现一种相对VLF事件更加低频且持续多天的信号(Ghosh et al，2015)，震中50km范围内MDS和BAX两个台站在4月5—18日水平向出现较明显的噪声值变化，周期在30~150s，而50km外的台站并未出现同步的明显变化；在变化幅度上，背景噪声每天有一定周期性变化，即白天噪声较大，晚上噪声较小。图 6为未经过平滑和过滤异常的水平向芦山地震前背景噪声100s频点时序变化，BAX台背景噪声变化每日低值提升了约20dB，每日高值增加了约10dB，MDS台背景噪声每日低值提升了约20dB，每日高值增加了约20dB。由于地震、爆破等能量主要集中于高频，基本不会对30~150s的频带产生持续数天的影响，因此也可以排除这两个干扰因素。

3 结论

本文采用大数据技术，基于分布式计算框架Spark实现了对海量波形数据噪声功率谱的分布式计算，通过对2013年4月20日四川芦山7.0级地震震前多频点背景噪声变化趋势进行分析，发现震前50km范围内的MDS台和BAX台出现相对VLF信号频带更低的长周期异常信号，表现为震前约半个月，在水平方向30~150s周期背景噪声上出现10~20dB增强且持续多天的情况，垂直向未发现明显变化，同时超出50km范围内的台站也未有明显变化。

从相关研究来看，这种长周期信号的异常变化很难通过传统的地震动信号拾取发现，本文提出的基于背景噪声的多频点分析方法，可以从高频50Hz到长周期200s的范围开展分析，通过大数据技术应用可以快速对海量数据完成计算分析工作，对于发现台站记录数据的长周期异常变化，以及进一步开展强震临震前的异常分析发挥一定作用。