胡宝慧，张 浩，教智博，雷凯悦

（黑龙江省地震局，黑龙江 哈尔滨 150090）

0 引言

地震台站的噪声研究受到国内外专家广泛关注。Peterson[1]通过分析全球75 个地震台站近2000 条噪声记录的功率谱密度分布，给出全球低噪声模型（NLNM）与高噪声模型（NHNM），目前广泛应用于台址噪声水平评价与不同背景噪声水平下地震计响应的预测。自2004 年以来，McNamara[2]认为，噪声功率谱密度（PSD）估计能较好确定台站噪声水平，起到实时监控和及时发现仪器故障的作用，国内许多专家[3-4]对以上方法进行了应用，取得了较好的效果。但现有噪声水平计算程序存在安装设置困难、数据格式受限、参数配置复杂等诸多问题，长期缺乏一套具备自动数据接入、处理运算和绘图功能的数据处理软件。大部分软件程序还需要手动收集处理观测数据的情况，存在速度慢、效率低、开发后难维护等问题，远远满足不了日常业务工作的需要。

地震台网监测能力，是地震台站布局及单台监测能力的综合体现。地震台网监测能力除了与台站的布局情况有关之外，还与台站的背景噪声水平、观测系统响应灵敏度、仪器动态范围有关。本程序基于地震噪声水平和震级衰减关系的理论监测能力评估方法[5]，以台网内所有台站噪声水平计算结果作为输入，依据近震震级计算公式，建立震级大小和震中距的对应关系，按照0.1°×0.1°划分网格得到整个台网的理论监测能力。其结果可用于分析区域内测震台网地震监测能力变化。

本次开发全部程序采用解释性脚本语言，极大简化了“开发、部署、测试和调试”的周期，为自动化单台噪声水平计算和台网的理论监测能力的快速估算提供了可能。

1 设计思路与流程

为了更好的达到程序的高可移植性、跨计算机运行等各项功能。本程序能够实现在Linux操作系统下，通过CLI（命令行界面）或者GUI（图形用户界面）运行。

程序主要包含文件处理、数据处理和绘图三个模块。其中文件处理模块完成数据的接入、小时波形数据格式的转换、文件合并、文件重命名等工作；数据处理模块完成仪器信息文件的解析、噪声水平的计算和理论监测能力的计算；绘图模块完成各类图件的快速绘制。模块间采用函数调用方式实现。

图1 程序设计思路Fig.1 Program design ideas

图2 程序计算流程Fig.2 Program calculation flow

程序编制过程中，根据噪声水平和台网监测能力的设计思路和计算流程，通过JOPENS系统中AWS 服务获取所需台站该时间段内小时波形数据服务，然后采用Perl 脚本程序调用SAC 对波形数据进行自动批量处理。首先将SEED 格式的数据转换成SAC 格式，然后将可能因多种因素出现间断的同一通道的SAC 数据合并起来完成小时波形的检查和数据合并，并对文件进行统一重命名。该流程将得到各台站SAC 格式的小时波形文件和对应的文本格式的仪器信息文件，其中通过解析仪器信息文件获得的数据采集器转换因子、地震计灵敏度、零级点参数等信息，可用于建立传递函数模型。

SAC 提供的命令可以实现地震数据的预处理，但无法实现所有的数据分析功能。因此在数据处理阶段，本程序使用SAC 中自带的I/O接口程序ReadSac.m 读写SAC 文件的MATLAB脚本。通过MATLAB 脚本程序实现SAC 格式文件固定长度的头段区和非固定长度的数据区的读取，进而开展数据的计算[6]。

2 噪声水平和台网理论监测能力自动化程序计算与实现方法

2.1 噪声功率谱（PSD）的计算

在噪声功率谱（PSD）估计方面，本程序应用基于FFT 的非参数功率谱估计方法Welch平均周期法通过分段选取数据进行加窗求功率，再进行平均计算。

程序编制过程中，选用Hamming 窗函数，数据置信水平为95%。首先根据地震计传递函数，去除仪器响应，得到地面运动的速度记录；然后将1 小时的数据段分成14 小段，每小段数据的长度为1000 S，数据段重合率为80%，对于每一小段数据去均值与线性趋势；计算所有14 小段的速度功率谱密度，取平均值，得到该小时数据的速度功率谱密度随频率的分布；再次对得到的速度功率谱密度进行1/8 倍频程滤波，得到平滑且在对数坐标上均匀分布的速度功率谱密度；最后将速度功率谱密度转化为加速度功率谱密度，并与NLNM、NHNM 进行对比分析，分不同的频段研究加速度功率谱密度的特征。离散噪声信号的功率谱密度可以表示为[7-8]：

同时可将噪声的单位用dB 表示为：

2.2 功率谱概率密度（PDF）的计算

因台站的噪声水平随着外界环境变化，功率谱密度也随时间不断变化，为了评估地震台站的环境变化特征，在PSD 值计算结束后，采用McNamara 的方法对1/8 倍频程滤波后的功率谱密度进行统计。以1 dB 为间隔，因为大多数功率谱密度都在-200～80 dB 范围内，因此本文以-200～80 dB 为统计范围，然后计算中心频率f 范围内的所有功率谱密度的数量Nf，Ndf为f处功率谱密度在（d～d+1）dB 范围内的个数，则（f，d）处的概率密度为：P（f，d）=Ndf/Nf。然后计算值落在中心点频率记录段的数目，统计各周期PSD 在不同时间内取某一数值的概率，得到功率谱密度在-200～80 dB 上的统计分布，即为该时段对应的PDF。

2.3 地脉动噪声均方根值（RMS）的计算

地脉动噪声均方根值RMS 可以衡量地震台台基背景噪声水平[6]，是求解台网监测能力的基础。在计算RMS 时，对中心频点在1～20 Hz频带范围内所有PSD 值的振幅值取平方，再求其在该时段内的平均，然后求其平方根。对于离散信号，均方根计算公式为：

其中，n 是样本数，xi是第i 个样本的幅值。

2.4 理论监测能力的计算

理论监测能力方面，以台站为中心、最大记录距离为半径画圆，求3 个以上台站的交集，该交集即为某一震级的监测区。在程序编制过程中，把台网监测区划分为N×M 的网格，各网格点所在坐标为（Xij，Yij），（i=1，…，N；j=1，…，M）。假设网格点为震中，把台网的全部台站按下式计算出震级ML，然后将震级ML由小到大排序，当网格点处发生第四号震级的地震时，则至少有3 个台站可以记录到，定义第四号震级为该网格点的理论最小监测震级。在2.3 节计算获得频带1～20Hz 均方根值RMS的基础上，根据目前地震台网用速度型记录测定近震体波的公式计算得到各台可监测到的最小震级：

式中，Vs为台站记录的体波最大速度，f对应最大速度频率值。R（Δ）为量规函数，即震级起算函数[9]。

3 程序运行

选取黑龙江省及周边区域74 个台站222 个通道的连续波形数据，采用上述噪声计算方法，获得各台站的小时噪声计算结果。在程序实际运行中，首先对单台48 小时连续数据进行功率谱概率密度函数（PDF）的计算，以低噪声模型（NLNM）与高噪声模型（NHNM）参照线作为衡量该台站噪声水平的标准。然后通过任意小时波形数据与48 小时参照数据的比对判定该小时数据偏差，以衡量该小时波形的质量。

选取BWS 和HEH 两个台站（图3）的数据计算结果进行分析，其中BWS 台为深井台站，仪器型号为BBVS-60，HEH 台为地面台站，仪器型号为JCZ-1。两个台站的NS 分量加速度功率谱概率密度函数分布如图4 所示。

图3 小时波形噪声功率谱（PSD）检验结果Fig.3 Hourly waveform noise power spectrum inspection results

图4 连续48 小时概率密度函数PDFFig.4 PDF of seismic station for continuous 48 hours

对于固定台站来说，以上两个台站仪器参数正常，台基环境噪声水平在频域的分布特征是相对稳定的。BWS 台作为深井台本应观测环境优良，但图3（a）显示长周期噪声偏差明显，且48 小时均超出高噪声模型（NHNM），需要查明该台站是何种原因导致的长周期频段异常。HEH 台虽出现该小时微震频段噪声水平异常，但明显从图3（b）可以看出是因为发生地震所致，且48 小时PDF 显示该台站仪器整体稳定，噪声水平在合理区间内。

从理论监测能力运行结果看，除黑龙江西北部局部地区外，全省测震台网基本已达ML2.5 以上地震监测能力，且中东部地区总体监测能力较强，这与监控黑龙江内最大断裂—依舒断裂的目标较为一致。从图5 看出，监测能力不仅与台站密度增加有关，更与该区域台站记录质量好坏密不可分，HEH 区域因台站稀疏导致监测能力较弱，而BWS 区域台站稀疏与记录质量不佳更导致了监测能力的下降，其他如大庆及哈尔滨区域虽然台站密度较高，但数据质量不高。因此，地震监测能力计算对优化测震台网布局，合理布置台站，避免盲目增加地震台站数量起到重要作用[10]。

图5 黑龙江地震台网理论监测能力Fig.5 Theoretical monitoring capability of Heilongjiang Seismic Network

4 结论与讨论

地震台站噪声水平和台网监测能力自动化程序完成自动分析地震台站波形数据质量的工作，产出了相应结果的图件。实现了区域内台站波形的实时接收和处理、自动监控波形质量的目的，为评估台站的地震噪声水平和台网监测能力提供了充足且可靠的依据。经过长时间试运行表明，该套程序能够满足目前我省台站数据质量监控的绝大部分需求，自动产出结果能为后期的台网规划和台站效能评估服务。

对于地震台网的地震监测能力评估，本程序采用的震级估算法虽然能够较为客观地确定区域内最小控震能力，但仍受噪声质量的影响。今后将通过程序改进，增加依靠台网实际产出为依据的诸如PMC 等计算监测能力的方法，实现地震台网监测能力更科学的评估。而在噪声水平分析方面，如何实现与地震目录、仪器实时参数配合分析噪声功率谱，更精确表述仪器的工作状态也是该程序需要逐步改进的部分。