0 引言

与模拟地震仪相比，数字地震仪因其动态范围大、频带宽、分辨率高等优点，容易记录到包括台基噪声在内的诸多复杂噪声信号（张宇等，2013）。台基背景噪声就是其一，是由风、浪、交通、人为活动等引起地球表面的微小震动，容易对地震观测产生干扰，从而影响观测效果（许可等，2015）。由前人研究可知，台基背景噪声取决于地震台站的台基条件和观测环境的好坏，而台站（台网）的地震监测能力不仅受仪器性能的影响，还与台基背景噪声有关（裴晓等，2012；张宇等，2013）。理论上来说，地震台站观测能力较大程度地被台基背景噪声的存在所限制，现今最为灵敏的观测仪器，也不能清晰记录到小于台基背景噪声的微小信号（许可等，2015）。因此，分析数字地震台的噪声水平及特征、掌握数字地震台的观测动态范围，有利于改善并提高地震台站（台网）的监测能力（童汪练，2000；康建红等，2014；杨龙翔等，2015；沙成宁等，2016）。

1 台网现状

1999年北京市测震台网（下文简称北京台网）投入建设，2年后通过验收并投入使用，囊括一个台网中心和18个数字遥测台，2017年5月，接收原北京数字遥测地震台网（中国地震局地球物理研究所）10个无人值守数字遥测台站。目前，北京台网由一个台网中心和28个数字地震台组成，包括17个地面地震计观测站和11个井下地震计观测站（图 1），配备BBVS-120、BBVS-60、BBVS-60DBH和FSS-3DBH地震计，配置EDAS-24GN和EDAS-24IP数据采集器。其中，除天坛观测站采取CDMA传输方式外，其余观测站均采用光纤介质，以SDH/MSTP方式，将观测数据实时传输至测震台网中心。

2 数据获取

为客观真实地了解台站背景噪声现状，并保证地震计系统响应特性稳定（许可等，2017），选取各台站24 h无震、无明显干扰记录的波形数据。北京台网的台站记录速度值采用SEED格式，以各台站单通道24 h实时记录波形为计算样本，将其进行分道并转换为ASC II码格式文本。

由于多种原因，测震台站观测记录中会出现直流偏移现象，计算前需剔除。把整个记录长度求和除以记录总数，即得到记录的平均值，将其作为该记录的直流偏移量，逐点消除直流偏离，采用以下数学公式予以校正。

$ v = {v_i} - \frac{1}{N}\sum\limits_{i = 1}^N {{v_i}} $

(1)

式中，v_i（value）为每个采样点的值；N为采样点个数。

3 背景噪声估算

为评估北京市测震台网运行状态，需计算并分析观测动态范围、噪声水平及特征等。

3.1 观测动态范围计算

使用数字地震仪记录数据计算台站观测动态范围，由于该数据为电压数字值，在实际计算过程中，需转换为地动速度。实测地动速度值计算公式如下

$ V=\frac{N \times F_{\mathrm{s}}}{R \times K \times S} $

(2)

式中，V为实测地动速度值，单位m/s；N为实际记录数字值，单位counts；F_s（full scale）为采集器输入满幅电压，单位V；R为2ⁿ，n为ADC字长；K为数据采集器实际工作增益；S为地震计工作灵敏度，单位V·s/m。

1-20 Hz频带地动噪声均方根值（RMS）是衡量台基背景噪声水平的重要指标，可以用同一标准衡量不同来源噪声在主要地震观测频带的水平。计算1-20 Hz频带内地动噪声均方根RMS值，有以下2种方法。

（1）对原始记录进行1-20 Hz带通滤波，按下式计算RMS值。

$ N_{\mathrm{RMS}}=\sqrt{\frac{1}{n} \sum\limits_{i=1}^{n}\left(V_{i}-\overline{V}\right)^{2}} $

(3)

$ \overline{V}=\frac{1}{n} \sum\limits_{i=1}^{n} V_{i} $

(4)

式中，N_RMS为RMS值；${\bar V} $为速度均值，单位m/s；V_i为某点实测地动速度值，单位m/s。

（2）根据GB/T 19531.1-2004《地震台站观测环境技术要求》，已知功率谱密度PSD，计算频带1-20 Hz内的RMS值。计算公式如下

$ N_{\mathrm{RMS}}=\sqrt{2 \mathrm{PSD}_{v} \cdot f_{\mathrm{c}} \cdot \mathrm{RBW}} $

(5)

式中，PSD_v为速度功率谱密度，单位(m/s)²/Hz；f_c为分度倍频程中心频率；RBW为相对带宽，计算公式如下

$ \mathrm{RBW}=\left(f_{\mathrm{H}}-f_{\mathrm{L}}\right) / f_{\mathrm{c}} $

(6)

式中，f_H为分度倍频程上限频率；f_L为分度倍频程下限频率。则式（5）可转化为

$ \mathrm{RBW}=\sqrt{2 \mathrm{PSD}_{v} \cdot\left(f_{\mathrm{H}}-f_{\mathrm{L}}\right)} $

(7)

台站的观测动态范围，可反映观测仪器自身性能和台基环境噪声背景水平，且有效动态范围大小反映了该台站可以记录到地震信号的最大能力。然而，由于台站环境噪声与地震信号叠加，仪器往往达不到所设计的动态范围，只能达到有效测量动态范围，计算公式为

$ D=20 \lg \frac{F_{\mathrm{s}}}{K \times S \times N_{\mathrm{RMS}} \times \sqrt{2}} $

(8)

由上述公式计算北京市测震台网数字地震台观测动态范围，结果见表 1。

3.2 台基噪声功率谱密度计算

功率谱密度（power spectral density，PSD）是描述随机信号（随机信号只能用统计平均量来表征）统计规律的重要特征参量，其物理意义为随机信号功率在频率上的分布，表现为一条功率密度和频率的关系曲线，计算方法包括直接法、间接法和改进的直接法等。文中采用改进的直接法，即Welch算法计算PSD值。

3.2.1 计算方法

直接法又称周期图法，是把速记序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅里叶变换X(k)，即

$ X(k) = \sum\limits_{n = 0}^{N - 1} x (n){{\rm{e}}^{ - j\frac{{2\pi }}{N}kn}}\quad (k = 0,1,2, \ldots ,N - 1) $

(9)

取X(k)幅值的平方，除以N，作为序列x(n)真实功率谱的估计，即

$ P(k)=\frac{2 \Delta t}{N}|X(k)|^{2} $

(10)

直接法的优点是不需要进行自相关计算，算法较简单，但同时其不满足一致性估计条件，分辨率和方差难以同时取得较高的参数，因此需要改进（邢晓晴等，2018）。Welch算法是对周期图法的一种改进，它包括对信号重叠分段、加窗处理，求出各段功率谱后求平均功率谱, 其具体步骤如下。

（1）各分段间可重叠，以减小方差。设随机序列数据总长度为N，按50%重叠的方式分k段，每段长度为L，则各段数据可以表示为

$ x_{i}(L)=\left[x\left(\frac{L}{2} i\right), \ldots, x\left(\frac{L}{2}(i+2)\right)\right] \quad(i=0,1, \ldots, k-1) $

(11)

（2）选择适当的窗函数W(k)，在周期图计算前直接加窗平滑处理，并计算各段频谱。由此，无论何种窗函数均可使谱估计非负。

$ {X_i}(m) = \sum\limits_{n = 0}^{L - 1} {\left[ {W(k) \cdot {x_i}(k)} \right]} {{\rm{e}}^{ - j\frac{{2\pi }}{N}mn}} $

(12)

（3）各段周期图估计公式如下

$ S_{i}(m)=\frac{1}{U L}\left|X_{i}(m)\right|^{2} $

(13)

其中，$ U=\frac{1}{L} \sum\limits_{n=0}^{L-1}[W(n)]^{2}$。

（4）将k段结果进行平均，得到估计的信号功率谱。公式如下

$ \hat S(m) = \frac{1}{k}\sum\limits_{i = 0}^{k - 1} {{S_i}} (m) $

(14)

Welch算法属于基于FFT的非参数功率谱估计方法，其优点是直观、便于理解和计算，适用于较长的数据记录（王俊等，2009）。

3.2.2 PSD计算

根据实际仪器种类、零极点和增益参数，基于以上算法，使用中国地震局港震机电技术有限公司童汪练研制开发的软件，计算北京市测震台网28个数字地震台地动噪声加速度功率谱密度值，单位(m²/s⁴)/Hz，以分贝表示（表 1）。

由于台基情况较复杂，将台基地噪声按观测频段进行高低分类，而非按全频段划分（任枭等，2004；张宇等，2013；高伟亮等，2015）。按照台基地噪声高低分类标准（GB/T 19531.1-2004《地震台站观测环境技术要求》），将北京台网28个测震台站的台基地动噪声按1-20 Hz观测频段（表 1）划分为5类，其中：Ⅰ类台基噪声水平为Enl＜3.16×10^-8 m/s，台站有刘斌堡、琉璃庙、马坊、密云、西拨子、斋堂、马道峪、四座楼、南山村、周口店、上方山；Ⅱ类台基噪声水平为3.16×10^-8 m/s ≤Enl＜1.00×10^-7 m/s，台站有大灰场、牛口峪、十三陵、西集、喇叭沟、太师屯、龙泉寺、东三旗、大兴；Ⅲ类台基噪声水平为1.00×10^-7 m/s ≤Enl＜3.16×10^-7 m/s，台站有八宝山、次渠台、凤河营、天坛、杨镇、榆垡；Ⅳ类台基噪声水平为3.16×10^-7 m/s ≤Enl＜1.00×10^-6 m/s，台站有金盏、牛栏山；Ⅴ类台基噪声水平为1.00×10^-6 m/s ≤Enl＜3.16×10^-6 m/s。分析发现，使用不同型号的地震计对台基噪声产生的影响不同，本研究中使用FSS-3DBH型地震计的台站，台基噪声水平普遍较大。未来可对台基噪声较高台站进一步落实，详细分析其台基条件以及周围观测环境对台基噪声的影响。

北京市测震台网各台基噪声背景优势频率具有不同特征，部分台站高频干扰明显，部分台站低频干扰增强，且长周期和高频段噪声功率谱曲线差异大。文中仅给出琉璃庙、密云、四座楼、马道峪等4个Ⅰ类台基的地震台PSD曲线，见图 2。由图 2可见，在1-20 Hz范围内条带能量平稳、较低且位于Peterson（1993）给出的全球地震背景噪声模型的新高噪声模型（NHNM）和新低噪声模型（NLNM）之间。

3 结论

通过对北京市测震台网各台基背景噪声数字化记录进行计算和分析，得到台基背景噪声地动速度均方根值（RMS）、有效测量动态范围和噪声功率谱密度（PSD），结果表明：①根据《地震台站观测环境技术要求》，北京市测震台网现有11个Ⅰ类台、9个Ⅱ类台、6个Ⅲ类台、2个Ⅳ类台；②北京市测震台网各台基噪声背景优势频率各有特征，高低频段噪声功率谱曲线差异大；③不同场地类型对台基噪声的影响有所不同，井下台站的台基噪声水平普遍较大，应与井下台站一般选择沉积层厚、噪声较大的平原区建台有关。

应把各个测震台站的环境地噪声水平分析作为日常工作的一部分，随时关注观测环境变化，为台站的环境保护起到一定促进作用。