基于功率谱概率密度函数检测测震台网JOPENS系统数据异常

安 全 张剑影 翟 浩 苏日亚 赵铁锁

（中国呼和浩特 010010 内蒙古自治区地震局）

0 引言

地震观测数据是国家防震减灾工作的基础性战略资源。数据质量直接影响防震减灾工作效果。高质量数据是决定地震台站与地震台网成为国际一流的先决条件。因此，保障地震观测数据质量是地震监测工作的核心任务（郑秀芬等，2017）。Peterson 等（1993）通过对全球正常背景噪声功率谱（PSD）的研究，给出了全球低噪声新模型（NLNM）和高噪声新模型（NHNM）。目前，该模型被广泛应用于地震台站噪声水平评价。McNamara 等（2005）提出了地震噪声功率谱概率密度函数（Probability Density Function，PDF）方法。该方法可用于台站噪声水平和波形质量的测定。PDF 方法直接使用连续波形数据，不对地震数据进行筛选，因此，可同时得到地震事件波形、观测系统仪器响应参数变化及仪器故障等概率分布。目前PDF 方法已被用于意大利（Marzorati et al，2006）、新西兰（Rastin et al，2012）、伦敦中部地区（Green et al，2017）的地震环境噪声特征分析。在国内，应用PSD 方法，孙晴等（2007）分析了河北省地震局数字地震台网子台地动噪声水平，岳秀侠等（2009）分析了天津测震台网子台地脉动；应用PDF 方法，廖诗荣等（2008）实现自动处理地震台站勘选数据，徐嘉隽等（2010）检测了地震仪器系统工作状态，林彬华（2015）实现实时监测地噪声，李小军等（2019）分析地震计自噪声水平。王俊等（2013）结合JOPENS 系统，利用MATLAB 平台研发了运用PDF 方法自动计算台站背景噪声的程序，从而实现了对数字地震台网观测系统运行状态以及台站背景噪声水平的准实时监控。

造成地震观测数据异常的因素主要有观测环境、观测设备运行状态、JOPENS 系统MySQL 数据库仪器响应参数（灵敏度、零极点、归一化因子等）准确性。在国内相关研究中，关于JOPENS 系统MySQL 数据库仪器响应参数有误导致观测数据异常的文献较少，笔者基于功率谱概率密度函数，来检测内蒙古测震台网JOPENS 系统的数据异常。

1 波形数据质量分析软件部署

以内蒙古测震台网JOPENS 系统为基础，以Continuous Wave Quality Look（CWQL）软件为支撑，从JOPENS 系统MySQL 数据库实时读取48 个台站仪器响应参数，从流服务获取48 个台站实时波形数据，准实时自动计算各台站各通道PSD 值，自动绘制PDF 图，对内蒙古测震台网JOPENS 系统观测数据质量进行评估。数据处理流程如图1 所示。

图1 数据处理流程Fig.1 The data access flowchart

波形数据质量分析CWQL 软件，结合了JOPENS 系统，利用MATLAB 平台实现了基于PDF 方法的背景噪声自动计算，从而实现了对数字地震台网数据记录质量、观测系统运行状态以及台站背景噪声水平的准实时监控（王俊等，2013）。

2 原理

2.1 功率谱密度计算

周期时间序列y(t)的有限范围傅里叶变换可表示为

式中，Tr为时间序列段长度，f为频率。

对离散频率值fk，傅里叶变换定义为

式中，fk=k/(NΔt),其中k=1,2,3，…，N-1，Δt为采样间隔（0.01 s），N=Tr/Δt为截取时间段的采样点数。

根据维纳—辛钦定理，功率谱密度（PSD）定义为

将速度PSD 值转换为加速度PSD，公式如下

需要扣除仪器传递函数影响，以反映真实地动噪声物理量值，公式如下

式中，Pa(f)为真实地面运动加速度功率谱，H(f)为仪器传递函数。

2.2 平滑处理

为了使PSD 在对数坐标中等间隔分布，采用1/3 频率积分的平滑处理方法

其中，fl=2-1/6fc，为低频拐角频率；fh=21/6fc，为高频拐角频率；n为介于二者之间频率f的个数。由式(6)得到中心频率fc的Pa(f)平均值Pa(fc)。以Pa(fc)作为fc的加速度功率谱密度，中心频率fc以1/9 频程为增加步长，即下一个中心频率为原中心频率的21/9倍，重新计算相应的fl和fh，并将新的fl和fh之间的PSD 平均值作为下一个中心频率fc的PSD 值。这样，在fc的取值范围0.02—40 Hz 内，每个记录段的PSD 值可由在对数坐标系中间隔分布的中心频率的PSD 值表示。

2.3 概率密度函数计算

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

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

台站环境背景噪声水平的速度均方根值（RMS）计算公式为

式中，RBW=(fh-fl)/fc为相对宽度。

3 台站地动速度记录

地动数据接入测震台网JOPENS 系统流程如图2 所示。总体上分为地震计输出y(t)和数据采集器输出g(t)，其中地震计输出y(t)为数采输入信号。

图2 数据采集器接入流程Fig.2 The flowchart of data

图2 中，Y(f)为y(t)的傅里叶变换，X(f)为x(t)的傅里叶变换，H(f) 为地震计传递函数，H(f)可表示为

其中，Sd为地震计灵敏度，A0为归一化因子，(z1，z2，…，zm)为m个零点，(p1，p2，…，pn)为n个极点。

测震台网JOPENS 系统存储的数据g(t)单位为counts，将数据还原为速度的频域计算公式为

其中G(f)为g(t)傅里叶变换，K为数据采集器转换因子倒数（灵敏度）。JOPENS 系统接入的波形数据为无量纲数值，并非真正的地动速度或地动位移，必须根据公式（10）扣除观测数据仪器响应（传递函数）才能得到地动速度，积分后得到地动位移。而K和H(f)以灵敏度（Sd）、零极点（zm，pn）、归一化因子（A0）等的形式存储于JOPENS 系统MySQL 数据库数中，如图3 所示。上述参数输入有误或者参数的变化均会影响观测数据质量。上述错误或变化以波形浏览等方式无法分辨，而应用PDF 方法生成的任意时间段PSD 曲线、PDF 图则可以直观反映观测数据“健康”状态。

图3 JOPENS 系统台站参数树状结构Fig.3 Parameter tree of a station in JOPENS system

4 数据计算及分析

通过CWQL 软件，从内蒙古测震台网中心JOPENS 系统读取48 个台站仪器响应参数和原始波形，并进行预处理。对波形数据进行去均值、去长周期成分处理，将原始波形分段，计算PSD 值并进行平滑处理，计算PDF 值及RMS 值，画PDF 图。把计算结果存储于本地。计算结果可按时间段、台站进行分项调取、分析。

选取内蒙古测震台网48 个台站7 月20 日14 时至2019 年7 月27 日14 时JOPENS 系统内的PDF 值评估样本（168 h，共计24 192 条）和1—20 Hz 频段垂直向RMS 值数据，从样本中挑选典型例子进行分析。

参照《地震台站观测环境技术要求第1 部分：测震》（GB/T 19531.1—2004），以48个测震台站垂直向RMS 值绘制台站噪声分级空间分布图（图4）。由图4 可知，内蒙古测震台网48 个台站背景噪声水平属于Ⅰ类的有45 个，属于Ⅱ类噪声水平的台站有3 个，整体噪声水平较低。Ⅱ类噪声水平台站均位于市区，噪声源主要为人为活动、交通等。

图4 48 个台站背景噪声分级空间分布Fig.4 Background noise classification interspace distribution of 48 stations

图5 为以乌加河台NS 和EW 向连续记录为样本的计算结果，可见两分向PDF 呈明显差异。在图5（a）中，NS 向PSD 曲线在整个频段呈上下2 条曲线。在图5（b）中，低于0.06 Hz 频段PSD 曲线呈2 条曲线。

图5（a）出现上下2 条曲线，是由于JOPENS 系统MySQL 数据库中的地震计传递函数H(f)的归一化因子A0比仪器实际归一化因子小1 个数量级。源数据扣除仪器响应时，降低了H(f)的值，从而高估了PSD 值。下面曲线为A0更正后正常PSD 曲线。在图5（b）中，低于0.06 Hz 频段PSD 曲线呈2 条曲线是由于MySQL 数据库中H(f)零极点（zm，pn）与实际值不一致。扣除仪器传递函数影响时，计算地震计频带宽度有误。参数更正后，低于0.06 Hz 频段PSD 曲线（下）值明显降低，回归正常水平。

图5 乌加河台水平向功率谱概率密度函数（PDF）分布Fig.5 Probability density function（PDF）distribution of horizontal components of the Wujiahe station

图6 为以霍林河台EW 向和呼和浩特台NS 向连续记录为样本的计算结果。图6（a）中PSD 曲线呈上下2 条曲线，且不在全球低噪声新模型（NLNM）和全球高噪声新模型（NHNM）范围内，曲线形状也与NLNM 和NHNM 不一致。图6（b）中呼和浩特台NS 向PSD 曲线基本在NLNM 和NHNM 范围内，形状基本一致，但在1—10 Hz 范围内PSD曲线明显凸出。

图6 霍林河台东西向和呼和浩特台南北向功率谱概率密度函数（PDF）分布Fig.6 Probability density function（PDF）distribution of the EW component of the Huolinhe station and the NS component of the Wujiahe station

图6（a）中，下面的曲线是因为数据采集器未正常连接地震计，没有采集到地面运动而形成的；上面的曲线是因为工作人员连接地震计时忘记调地震计机械零点，地震计处于靠摆临界状态，地震计记录地面运动的动态范围缩小而形成的。图6（b）中，在1—10 Hz 范围内，PSD 曲线明显高出是由于呼和浩特台距市区较近，受车辆、人为干扰较大。只有台站迁址才能有效降低该频段噪声水平。

5 结束语

利用观测数据功率谱概率密度函数，检测测震台网JOPENS 系统数据异常，可以得出以下结果：①可以监控台站运行质量，包括地震计和数采运行状态；②能检测JOPENS 系统MySQL 数据库仪器响应参数准确性；③能够实时评估台站背景噪声水平。

经过地震烈度速报与预警工程项目的建设，内蒙古地区台站数量不断增加，如何高效维护台站、实时检测、及时发现仪器设备故障是目前面临的难题。使用功率谱概率密度函数方法评估台站背景噪声水平，可以帮助运维人员快速发现观测系统中存在的异常信息，从而缩短运维周期，在一定程度上解决目前台站数量增多、运维能力不足的问题。