地震噪声异常实时监测

林彬华 金星 廖诗荣 李军 黄玲珠 朱耿青

1)福建省地震局，福州市华鸿路7号 350003

2)福州大学，福州市闽侯县大学城学园路2号 350002

3)中国地震局工程力学研究所，哈尔滨 150080

0 引言

地震波形数据的质量很大程度取决于地震观测仪器系统是否正常工作。为了检测地震仪器是否发生异常，通过每天对其进行脉冲标定，可以发现大部分故障。但是这种方法需要人工操作判断，工作量大、实时性差，很难及时发现异常。除此之外，对于要求分秒必争的地震预警，要实时监测数据质量状况，为后续到来的地震，快速计算震级和定位提供客观、科学的数据。因此目前迫切需要一种新手段来实时监测地震记录系统是否能正常工作，从而保证地震数据的最优质量。

国际上，McNamara等(2004、2005a)提出用功率谱概率密度函数(PDF)方法进行台站地震噪声水平监测，并在GSN与ANSS等台网应用于日常仪器工作状态检测。廖诗荣等(2008)选取台站48小时的噪声记录，并以160s为单位对噪声记录进行分段，然后采用概率密度函数(PDF)方法求出台址的总体噪声水平，最后完成了自动化处理地震台站勘选测试。Sleeman(2007)在Orfeus数据中心通过监测实时波形数据中地震噪声加速度功率谱密度(PSD)值的变化，在VEBSN(Virtual European Broadband Seismograph Network)实现了对宽频带地震台网波形数据质量的自动检测。徐嘉隽等(2010)研发了数据质量检测软件，其主要原理也是利用概率密度函数方法计算一天的噪声记录，然后针对自动绘制出的PDF图及RMS值，直观地判断观测数据的质量，并已经应用于日常地震观测系统数据质量的检测。为了能够实时监测地震仪器工作状态，本文采用功率谱概率密度函数(PDF)方法绘制出台站的PDF图，再用网格概率方法确定台站高低噪声参照线，进而研究噪声异常的遴选方法，最后形成地震噪声实时监测系统。

1 地震噪声概率密度函数方法概述

可将地震噪声视为平稳随机信号，依据随机过程理论，通常用概率统计方法来描述随机信号。平稳随机信号在时间上是无限的，其能量也是无限的，但其功率却是有限的。平稳随机信号的功率谱反映了信号的功率在频域内随频率的分布，也称功率谱密度(简称PSD)。

Peterson等(1993)通过分析全球75个地震台站近2000条噪声记录的功率谱密度分布，给出了全球低噪声模型(NLNM)与高噪声模型(NHNM)，该模型目前广泛应用于对台址噪声水平的评价、仪器标准定义以及不同背景噪声水平下地震计响应的预测等。需要指出的是，Peterson用来进行建立地震噪声模型的2000条噪声记录是从近12000条波形记录中筛选出来的，其它近10000条包含地震、爆破、标定等非地震噪声的记录段并未参与最后的处理。

为了减免非地震噪声记录段的筛选与截取这个波形预处理环节，对包括非地震噪声记录段在内的原始连续波形记录直接处理，McNamara等(2005a，2005b)提出应用PSD概率密度函数(PDF)方法进行地震噪声PSD值的计算。该方法的主要思路是:将原始波形数据分成n个记录段，采用与Peterson相同的方法对每个记录段计算PSD值，使用1/3倍频程的频率间隔对每个记录段PSD曲线进行平滑;然后计算PSD值落在某一个频点某一功率窗内的记录段数目，以该记录段数目与总记录段数目n的比值作为该频段该功率窗的PSD概率密度函数的取值。图1(a)是对YDXS台站194545个原始噪声波形记录段(1年数据)使用概率密度函数方法处理得到的PDF图。在该图中地震噪声以高概率取值的形式出现，并未因人工截除的地震体波与面波、仪器尖脉冲干扰、脉冲标定等而出现低概率值形式，不会影响对高概率水平下环境地震噪声水平的评估。通过PDF图还可以检测记录系统的故障和全面评判台站的数据质量。

2 福建地区地震噪声观测资料与处理

2.1 地震噪声资料收集与预处理

福建省测震台网由原先的41个测震台站和《十一五》新建的44个测震台站组成，图2给出了福建省这85个测震台站的分布。本文所研究的地震噪声资料是福建省85个测震台站2012年1～12月噪声的垂直向记录，采样频率为100Hz(李军等，2011)，其中大部分仪器所记录的频段在1/60～50Hz之间。考虑到频段前后会出现一定的衰减，所以本文设定PSD值的频率取值范围为1/50～40Hz。

图1 功率谱概率密度函数(PDF)分布图。 (a)利用台站1年的噪声数据计算出来的PDF图，点线为确定的高、低噪声参照线;(b)对应三维PDF图中1Hz频点处的功率谱密度概率分布图

福建台网的地震记录通常以EVT格式或SEED格式存储，数值的单位为counts。使用时都需要对数据做如下处理，其中地动速度(李军，2007)的计算公式为

对一些漂移的波形采用最小二乘法作基线校正，具体公式为

式中v(t)为未校正的地面速度记录，v*(t)为校正后的地面速度记录，参数a、b是对每条记录v(t)进行最小二乘拟合求得的，其计算公式为

2.2 地震噪声功率谱计算过程

利用噪声记录计算噪声功率谱的数据处理过程主要步骤如下:

(1)记录段的选取。对1年期的噪声记录以每5分钟(300s)进行分段，且为了避免PSD值的变化过大，对每段记录间按50%的叠加率选取。记录段长度Tr取决于我们感兴趣信号的最长周期TL，通常Tr的取值需达到TL的6倍以上(Bormann，2002)。因此取记录段长度300s，可以满足反映50s低频地震噪声PSD值的需要。

(2)速度功率谱密度值估算。最常用的估算地震噪声PSD值的方法是直接傅里叶变换法。该方法通过计算有限长度数据序列的FFT变换来计算PSD值。为了得到以频率为自变量的FFT变换，可以将数据序列表示为

式中，f=k/NΔt(k=0，1，…，N-1);N为采样点个数;Δt为采样时间间隔。y(nΔt)是每段原始噪声记录在时域中的数据系列，Y(f)是傅里叶变换得到的振幅谱。

功率谱反映了信号的能量随频率的分布情况，功率谱与傅里叶振幅谱二者在本质上并无多大不同，只是功率谱的纵坐标大体上相当于傅里叶振幅谱纵坐标的平方，其表达式为

(3)加速度功率谱密度的计算。当频率为f时，速度功率谱密度PSDv与加速度功率谱密度PSDa的转换公式为

其中:PSDv的单位为(m/s)2/Hz;PSDa的单位为(m/s2)2/Hz。

项目建设高度重视资源整合。首先，在项目设计阶段，将水文系统现有水文站网、中小河流监测建设及规划站网纳入山洪灾害防治非工程措施暴雨洪水监测站网，充分利用了现有监测站网资源，避免重复建设。项目实施过程中尽可能避免与气象站重合。二是充分利用防汛抗旱指挥系统建设成果，项目建设要求在完全满足国家防总对山洪灾害非工程措施建设技术要求的同时，从监测数据流程规划、数据采集、传输和存储的技术标准等均保持了与防汛抗旱指挥系统标准一致。三是山洪灾害预警平台与防汛抗旱指挥系统进行了完整的融合，山洪预警作为防汛抗旱决策支持系统的一个模块，极大地方便了各类应用。

(4)平滑处理。为了得到PSD在频率对数坐标中呈等间隔采样，本文采用1/3倍频程积分作为平滑处理

其中，fl=2-1/6fc为低频拐角频率;为fh=21/6fc高频拐角频率;n为介于二者之间频率f的个数。由(7)式得到中心频率fc的PSDa平均值PSDa(fc)，作为fc的加速度功率谱密度的PSD值。中心频率fc以1/9倍频程为增加步长，即下一个中心频率fc=21/9fc，重新计算相应的fl和fh，然后将新的fl与fh之间的PSD值平均值作为下一个中心频率fc的PSD取值。这样在fc的取值范围0.02～40Hz内，每个记录段的PSD值随频率变化情况可由在对数坐标系呈等间隔采样的107个中心频率的PSD值来表示(图3)。

2.3 PDF图上确定台站高、低噪声参照线

每个中心频率fc的PSD概率密度函数为

其中，Nfc为fc频点的记录段总数;Npfc为fc频点的PSD值落在某PSD取值范围内的记录段个数，在本研究中PSD窗长与步长都取1dB，变化范围为-200～-50dB。然后，以频率为横坐标、以PSD为纵坐标、以PPSD(fc)色块颜色深浅绘制三维平面图，得到功率谱概率密度函数(PDF)分布图(图1(a))，不同色块代表某频点在一定PSD窗内功率谱概率数。

从PDF图中可以明显看出台站不同噪声水平所处的区域。如图1(a)所示，颜色较深部分(蓝、绿、黄、红区域)是台站正常噪声水平落入的区域。颜色较浅部分(粉红色区域)是台站异常记录(地震波、外界噪声)所处的区域。为了能描绘出正常噪声所处的区域范围，本文引进了台站高、低噪声参照线来表示该区域的上下限。首先，选取蓝色轮廓作为正常噪声的界限，由图1(a)右边的图例可以看到蓝色轮廓线所对应的网格概率是1.5%。

针对图1(a)中1Hz处的加速度功率谱密度概率分布进行分析，以功率谱密度值为横坐标，网格概率为纵坐标，画出对应1Hz频点的功率谱密度概率图(图1(b))。以最大网格概率点为中心，向左、右两边分别找首个网格概率≤1.5%的点，用于确定低噪声参照点和高噪声参照点。对于其它频点采用同样方法确定高、低噪声参照点，最终将这些点连接起来就得到了台站网格概率≈1.5%的高、低噪声参照线(图1(a))。

图2 福建省85个测震台站分布图

3 噪声异常的分类与遴选方法

本文对噪声异常的遴选思路是利用台站过去1年的历史记录绘出台站PDF模型，由该模型确定出台站的高、低噪声参照线，然后将实时计算的PSD线与高、低噪声参照线进行比较，实现噪声异常自动化遴选。

根据PDF图可以直观的判断噪声是否发生异常。为了便于研究噪声异常的自动化挑选，本文将噪声异常大致进行了分类，分类规则如下:PDF图上显示的异常分布在低噪声参照线附近及以下，属于低噪处异常;PDF图上显示的异常分布在高噪声参照线附近及以上，属于高噪处异常;PDF图上显示介于高、低噪声参照线之间的异常，则属于中噪处异常;除此之外，缺数异常也是常见的一种噪声异常，一般都没有数据，在PDF图上无法显示，所以本文直接在时域上对其进行挑选。挑选方法:若所有点都小于1.0×10-10m/s(由于正常噪声速度记录一般为1.0×10-6m/s左右)，则判定为缺数异常。由于噪声异常具有不定性和多样性，那么只用一种挑选方法是无法完成对噪声异常的挑选，通过上述对噪声异常的分类，然后研究针对每个类型的挑选方法，就可以很好地解决这一难题。

3.1 低噪处异常的挑选

3.2 高噪处异常的挑选

从PDF图中也可以很容易判断出高噪处异常，对造成高噪声异常的某条波形，并绘出该波形的PSD线(图3(b)，实线)，从图3(b)中可以看出，该类异常的特征为PSD值大于台站高噪声参照线。同样考虑在确定台站高噪声参照线时存在误差影响，则由台站高噪声参照线全体向上移动3dB而得到虚线。挑选方法:若PSD线(实线)超出虚线的点数占总点数的百分比大于0.45，则判定为高噪处异常。

图3 噪声异常挑选方法示意图。(a)低噪处异常中5分钟一小段噪声记录计算出来的PSD线与台站低噪声参照线的比较;(b)高噪处异常中5分钟一小段噪声记录计算出来的PSD线与台站高噪声参照线的比较;(c)中噪处异常中5分钟一小段噪声记录计算出来的PSD线的特征

3.3 中噪处异常的挑选

同样可利用PDF图判断出中噪处异常，选取造成中噪声异常的某条波形，绘出该波形的PSD线(图3(c)，实线)，从图3(c)中可以看出，该类异常的特征是PSD线起伏不大，几乎成一条直线。由于正常噪声的PSD线起伏较大，那么它所对应的方差也较大，这样可以通过比较两者方差的不同来挑选出该类异常。挑选方法:先用最小二乘法对PSD线作倾斜校正，再求校正后的PSD线的方差。若方差＜3，则判定为中噪处异常。

3.4 挑选结果

若要将上述方法运用到实际监测中就需将其按一定的次序整合起来，本文给出了地震噪声实时监测系统流程图(图4)，根据该流程图用Matlab编写了对应的挑选程序。选取福建省85个测震台站2013年7月份的噪声记录数据进行验证，依据图4所示的系统流程分别对85个台站的噪声数据进行在线运行。运行结果如图5所示，横坐标表示福建省85个台站，纵坐标表示每个台站应用地震噪声实时监测系统所挑选出来属于异常波形的个数与所挑选出来的波形总数的百分比，图中每个点表示一个台站所识别出来异常的正确率。从图5中可以看出85个台站识别出来的异常正确率均在90%以上，多数台站达到100%。这是因为噪声一旦发生异常，那么其对应的功率谱密度PSD值也是异常的，这样利用台站高、低噪声参照线作为参考标准就可以很容易的将异常挑选出来。

图4 地震噪声实时监测系统流程图

图5 选取2013年7月份噪声记录对福建省85个台站进行测试的结果

4 结论与讨论

本文利用台站1年期的历史记录绘出台站PDF模型，然后由该模型确定出台站的高低噪声参照线，再根据这两条线与实时计算的PSD线进行比较，实现噪声异常自动化挑选。该方法用于噪声实时监测可以取得很好的效果，大大减少了台网仪器系统异常检测的工作量，且能及时发现异常，并提醒台网工作者及时解决故障，保证记录到的地震数据质量优良。本文研究可以得到如下结论:

(1)采用网格概率的方法来确定台站高低参照线并作为噪声异常遴选的标准，其中考虑了台站台基和仪器稳定性等复杂因素的影响，对研究台站本身噪声异常挑选具有很强的适用性。

(2)地震噪声实时监测系统对噪声异常的挑选正确率很高，且普遍适用于福建省85个台站。

(3)通过噪声异常在PDF图上所分布的区域不同，将噪声异常分成缺数(断记)异常、低噪处异常、中噪处异常、高噪处异常等4类，其中囊括了绝大部分异常。高噪处异常较为复杂，包括了仪器异常、地震波、外界噪声、系统瞬变等异常，这些异常有些是非仪器异常引起的，需要进一步研究对它们的区分方法。

(4)通过大量数据统计，给出高、低噪声异常挑选方法的阀值，这为以后新建台站的异常挑选提供了可行的参照标准。通过该参照标准可以很好解决新建台站历史数据不足的问题。