一般情况下,在背景噪声的计算研究过程中,记录波形中的地震信号、观测系统的瞬态变化,仪器的毛刺(如数据记录间跃、限幅、尖峰、脉冲标定等)都会被去除掉或被进行归一化处理(Campillo,Paul,2003; McNamara,Buland,2004)。但McNamara和Buland(2004)在研究美国大陆地区背景噪声水平时,认为记录过程中的这些信号相对于低水平的背景噪声是低概率的事件,在对台站的背景噪声水平进行评价时,可以将它们包含在内,并且认为这些信号是非常有用的,其概率值的大小可以用来评价地震观测系统的运行状态和记录数据的质量。为了实现定量化分析,本文以各个台站单通道每小时的实时记录波形作为计算样本,通过软件中的CTN模块自动从台网中心JOPENS系统的流服务器上读取原始波形数据(图1),并进行必要的预处理:(1)去除记录数据中的线性趋势;(2)将每条计算样本分段,并按50%进行重叠,每段的长度L参照Peterson(1993)的做法,并取为81.92 s(即L=2.∧13)。

功率谱密度的计算采用Matlab中提供的WELCH方法,该方法是改进后的分段平均周期法,属于功率谱密度的非参数估计方法(万永革,2007; 张志勇,2011)。为了最大程度地减小重叠后“频谱泄露”效应、增加频峰的宽度,采用了汉宁窗(Hanning)函数来提高频谱的分辨率,汉宁窗长度取N=L/2,各段的数据表示为

x_i(L)=[x(L/2i),…,x(L/2(i+2))](i=0,1,…,k-1).(1)

加窗平滑处理并计算各段频谱

X_i(m)=∑^L-1_n=0[W(k).x_i(k)]e^-j(2π)/Nmn.(2)

各段周期图估计为

S_i(m)=1/(UL)|X_i(m)|².(3)

其中, U=1/L∑^L-1_n=0[W(n)]².

求k段谱的平均

S^{^}(m)=1/k∑^k-1_i=0S_i(m).(4)

其中,每段信号的功率谱密度是运用FFT算法来实现,取每段信号功率谱的平均值作为整个计算样本的功率谱值。再根据仪器的传递函数在频域里扣除仪器响应,即得到每一条计算样本的速度功率谱密度结果。

可近似地认为数字地震观测系统具有因果性、时不变的特征(中国地震局监测预报司,2007),因此每一条记录信号反映的是记录系统在某一特定时刻的状态,其记录过程具有唯一性,因此信号出现的概率值代表观测系统运行中所呈现出的状态。于是,可以利用信号随频率的概率密度函数分布来对观测系统的运行状态进行评估(McNamara,Boaz,2005),具体步骤如下:

(1)将每一条功率谱密度在全频段内,以1/8倍频为单位间隔滑动地计算平均功率谱,在短拐角周期 T_S和长拐角周期T_L=2×T_S内的平均功率谱值,赋给几何中心周期T_C,

T_C=sqrt(T_S×T_L).(5)

T_S按1/8倍频增加,即T_S=T_S×2.∧0.125,对于采样率为100 Hz的观测系统来说,有效频率范围可以达到90 s～35 Hz。

(2)在某一时间段内,对于某一给定的中心频点,各种记录信号的功率谱概率密度函数(PDF)表示为

P(T_C)=N_pTc/N_Tc.(6)

N_pTc是功率谱值在-200～-80 dB 范围内1 dB单位间隔内的功率谱条数,N_Tc是在中心频点T_C上,单位时间内的功率谱总数。

(3)台站环境背景噪声水平的均方根值(RMS),根据Bormann(2012)提出的下式进行计算,

a_RMS={P×(f_u-f₁)}^1/2=(P×f₀×RBW)^1/2.(7)

式中,P为Matlab计算出来的双边功率谱密度,RBW=(f_u-f₁)/f₀为相对带宽,f_u为分度倍频程上限频率,f_l为分度倍频程下限频率。

为了使软件运行稳定且计算高效,本研究将系统中的各个功能进行分解,按模块进行设计,各个模块之间使用JMS进行消息传递,使用EJB、ODBC技术进行远程数据传递和访问数据库。通过将各个模块添加到Taskcontral定时任务控制器中,实现对各模块的联系、定时控制及管理,从而减少系统中各个功能模块之间的直接联系,使得整个系统中各个功能模块之间的逻辑部署变得简单而清晰,结构框架见图1。

图1 自动准实时背景噪声分析计算软件系统模块结构
Fig.1 The system module structure of automatic quasi-real-time background noise analysis and calculation software

CTN模块:通过Server2Server 服务,从台网中心JOPENS系统的流服务器上获取连续波形数据。

PDT模块:负责将在原始数据进行计算前的各种预处理。

CCN模块:是核心计算部分,并通过ODBC数据库连接技术将计算结果写入指定的数据库进行存储。

DEL模块:定时清理已计算过的数据和内存变量,以保证系统的稳定和持久运行。

GUI模块:负责根据用户的指定条件获取数据里的计算结果、进行显示、编辑等。

各模块的启动和实际运行是通过Taskcontral来控制的,如图2所示,用户将只需将各模块的路径和执行时间添加到任务列表中,控制终端就会自动检测各模块的执行时间,一旦设定的时间与计算机当前时间一致时,就立即启动相应模块,其余时间处于等待状态,设计允许添加的任务最多可达120个。

图2 TaskContral定时任务控制终端
Fig.2 The timing task control terminal by TaskContral

使用GUI客户端模块可以查看数据库中任意时段内某一台站、通道的计算结果,操作界面如图4所示。

查询结果如图5a所示,即某一时段内所有信号在90 s～35 Hz范围内的功率谱概率密度函数分布,用色标来表示台站记录信号随频率变化的概率值,颜色越红表示台站记录的信号在此频点处的概率越大。一般情况下,尽管低水平的背景噪声是观测系统中的高概率事件,但为了能更加客观、真实地评估台站的环境背景水平,本文中取全频带内的中值(图5a中灰白线所示),来代表

图3 CTN 模块自动从流服务器提取数据的过程(a)和进行背景噪声计算分析的过程(b)
Fig.3 Automatic process of extracting data from the streaming server(a),the process of calculating and analyzing ground noise(b)by CTN

图4 GUI模块的主界面(a)和查询条件设定界面(b)
Fig.4 Main interface(a)and query conditions setting interface(b)of the GUI module

图5 计算结果的查询显示(a),及基于Google Map显示的台基噪声水平等级划分(b)
Fig.5 Query display of calculation results(a),and grade classification of station noise level displayed by the Google Map(b)

该台在此时间段内的背景噪声水平值。为了实现对观测系统运行状况的定量化监控,本文以NHNM和NLNM值(Peterson,1993)作为标准,但考虑到区域数字地震台网主要负责区域内近震的监测,根据近震地震波的频率特性,软件中只以高频段的信号为目标对象,具体的做法:在1～20 Hz范围内的取1.2、3.1、4.5、10.0、19.2 Hz 5个频点的值作为基点,计算其高于NHNM模型和低于NLNM模型的概率值,来评估观测系统的健康状况,如图5a所示,表示前三岛台(QSD)在2012年11月1～30日期间约有6%的记录结果高于NHNM模型值,即处于异常状态。

为了使计算结果更加具有实用性,依据《地震台站观测环境技术要求》(GB/T1953.1-19531.4-2004),自动将各台站的台基噪声水平进行等级划分,并运用GoogleMap ApI 技术进行地图的制作和显示,实现与GoogleMap丰富的全球地图资源进行融合,具有美观的显示效果(Google Maps Javascript API V3,2012),如图5b所示,绿色三角代表Ⅰ类台基、黄色三角代表Ⅱ台基、红色三角代表Ⅲ类及以上台基。

为了方便应用,软件还提供数据及图像存储功能,数据结果以文本格式、电子表格保存便于对数据结果的查看与后续处理。图像可存储为png、jpg、ai、eps、emf 5种格式,便于成果报告撰写及发表文章时不同编辑的要求。

运用此项新技术,可以很方便地实现对台基背景噪声水平季节性变化的监控,以2011年江苏省40个宽频带数字地震台的背景噪声水平进行计算分析为例,图6为以各台站每天垂直向的连续记录为计算样本(共计11 972条)的计算结果,根据所有功率密度谱随频率的概率分布情况,取概率最大处的值来代表江苏地区的背景噪声水平值; 按季节来显示计算的结果:2011年全年内0.1～0.5 Hz和0.03～0.1 Hz范围内噪声能量基本相当,其中0.1～0.5 Hz范围内的噪声源更加稳定,0.03～0.1 Hz范围内的噪声则会随季节变化出约5 dB左右的扰动。

图6 2011年江苏省台站背景噪声的功率谱概率密度
Fig.6 The power spectrum probability density of background noise of stations in Jiangsu Province in 2011

在宽频带地震记录系统中,系统的某些瞬态变化、数据记录中的阶跃、通信链路的短时中断、或系统故障引起的信号失真(如记录中叠加了高频电流信号等),以及特定频率的环境干扰源等,系统维护人员通常情况下难以直观地发现。本研究中的计算方法对记录信号的最小分辨率约为0.012 2 Hz,经计算后可清晰、直观地辨别记录信号中不同频率信号间的细微差异。例如,南京台的BHE分向在2012年11月1～30日(共538条功率谱结果)存在1个频率约为2 Hz的干扰源,如图7a所示,但在此期间的日常地脉动记录波形中,即使在扫描时间为6 s的窗口内也几乎不可能被识别出来,图7b是2012年11月1日10时的实际记录波形,可以看到BHE分向的记录与BHZ、BHN分向的记录无明显区别。

图7d是高邮台2012年8月27日1时三分向记录的连续波形,可以看到三分向的记录基本一样,难以直观地发现BHZ分向的记录中有异常情况,但经过噪声计算分析后,结果显示:2012年8月26～29日期间(74条功率谱结果)的噪声概率密度分布呈斜线形(图7c),与正常的噪声概率密度分布完全不一样,最后经仪器维修人员确定为地震计BHZ分向故障。这些结果表明,通过对实时记录的噪声进行准实时的分析,可及时地 “跟踪”到地震观测系统的健康状态。

图7 南京、高邮地震台故障信号的概率密度函数及相应的连续记录波形
(a)南京台BHE分向的背景噪声概率密度分布;(b)南京台三分向记录的连续波形(红色为BHE分向);
(c)高邮台BHZ分量的背景噪声概率密度分布;(d)高邮台三分向记录的连续波形(红色为BHZ分向)
Fig.7 Probability density function of fault signal at Nanjing,Gaoyou Seismic Stations and their corresponding continuous recording waveform (a)PDF of background noise in BHE component recorded by Nanjing Seismic Station;(b)continuous waveform in three components recorded by Nanjing seismic station(the red line is the BHE component);(c)PDF of background noise in BHE component recorded by Gaoyou Seismic Station;(d)continuous waveform in three components recorded by Gaoyou Seismic Station(the red line is the BHZ component)