基于matlab的背景噪声计算程序的设计与应用

许　可，刘瑞瑞，孔繁旭，朱　宏，柳艳丽，李　晔

(天津市地震局，天津　300201)

基于matlab的背景噪声计算程序的设计与应用

许可，刘瑞瑞，孔繁旭，朱宏，柳艳丽，李晔

(天津市地震局，天津300201)

摘要：文章介绍了利用软件平台，开发出采用经典谱估计中算法的背景噪声计算程序的设计原理及实际应用。该程序的应用，解决了天津测震台网各台站背景噪声的统一计算、批量成图等问题，从而实现了对数字地震台网数据记录质量、背景噪声变化的监控以及各频段噪声特性的对比分析，提高了天津台网背景噪声计算的工作效率。

关键词：背景噪声；功率谱密度；天津测震台网

0引言

在数字化观测系统中，系统的瞬态变化、仪器毛刺、数据记录的阶跃和尖峰以及由于系统故障引起的信号失真和重大环境干扰源的出现等，都是影响数据记录质量的重要因素[1]。随着信息化社会的到来，人们的认识水平不断提高，逐渐关注实时记录的背景噪声，希望它能提供科学有用的信息[2]。在数字地震观测系统记录中，除真实的地震信号外，还包含了由于系统故障引起的各种失真信号以及台站周围的环境噪声干扰等信息[3]。这些干扰信息严重影响观测数据的质量，需要加以区分以及排除，确保观测数据质量[4]。因此，如何及时了解数字地震观测系统运行的状态，保证地震台网记录的数据质量以及识别干扰源，就成为当前数字化地震台网系统监控技术中的重要研究课题之一。

“十五”项目后，天津测震台网建成1个台网中心和31个地震台站。台网孔径东西约127 km，南北约160 km，地震台站平均台间距为30.6 km。台网中心采用JOPENS对数据进行接收、交换存储及处理。在经过“首都圈测震台网的升级改造”项目后，增添多种不同的仪器类型，导致地震台站的背景噪声计算步骤繁杂，需花费大量的时间。为此，本文利用matlab平台研发了运用经典谱估计中的Welch算法，解决天津测震台网31个台站背景噪声的统一计算、批量成图等问题，从而实现对数字地震台网数据记录质量、背景噪声变化的监控以及各频段噪声特性的对比分析。

1计算原理

1.1　去直流偏移

在数字地震记录中会出现直流偏移，这个信号不表示真实的地面运动，必须消除这一直流分量[5]。在实际计算过程中，把整个记录长度求和除以记录总数，取记录的平均值作为记录的直流偏离量。然后，逐点减掉直流偏离量作为第一步的校正过程，公式为：

式中：fi表示每个采样点的值；N为记录长度的总采样个数。

1.2　扣除仪器响应

数字地震观测系统是由地震计拾取地面运动，数据采集器对地震计拾取的信号进行数字化采集，通过数据专线传输到台网中心的数据接收系统。地动速度变换成数字记录的过程就是仪器响应的过程。在实际计算噪声功率谱的时候，要考虑观测仪器对地动噪声的影响，要想得到地动噪声的绝对量值，还应从计算结果中扣除仪器的影响，需利用数字地震仪器的传递函数来扣除仪器响应。数字地震仪器的传递函数[6]为：

式中：Sd代表仪器的灵敏度；A0是归一化常数；r1为复零点；Pm为复极点；L为零点个数；M为极点个数。

1.3　功率谱密度的计算

经典谱估计方法在工程中都是以离散傅立叶变换为基础的，它隐含着对无限长数据加窗处理，所以，经典谱估计有着分辨率不高、能量泄漏的固有缺点。为克服这些缺点，经过研究，提出了各种算法。目前，经典谱估计算法有周期图法、Barlett算法、Welch算法、Nattall算法等。Welch算法是由welch提出的对周期图的修正算法，是经典谱估计中获得有效应用的一种算法。

功率谱密度的计算采用经典谱估计中的Welch算法，Welch算法是在周期图法的基础上改进得到的。Welch算法谱估计采取数据分段加窗处理再求平均的办法，先分别求出每段的谱估计，然后进行总平均。即：把某一长度为N的数据x(n)分成L段，每一段为M(在分段时可允许每一段的数据有部分的重叠)，分别求每段的功率谱，然后加以平均。每一段的数据窗口可使用汉宁窗或哈明窗。这样可以改善由于矩形窗边瓣较大所产生的谱失真。用p(ω)表示用Welch算法估计的功率谱[7]。即：

式中：U为归一化因子；d(n)是数据窗。

2程序的matlab实现

matlab是一套用于科学计算的可视化、高性能的语言和软件环境，它集数值分析、矩阵运算、信号处理和图形显示于一体，已经成为数字信号处理应用中分析和仿真的主要工具，工具箱中包含各种经典和现代的数字信号处理技术，能实现各种数字信号的处理。

本程序基于matlab设计，主要包含两部分。一部分为程序的主函数，另一部分为程序的功能函数。程序的主要思想是先通过功能函数计算出功率谱，再通过主函数初始化参数调用功能函数进行绘图操作，实现统一计算、批量成图等功能。功能函数只负责计算出功率谱密度，具体的参数设定、绘图操作等在主函数中完成，这样既可以对计算结果分别成图，又可以对计算结果集中成图，或者可将所关心的某几个计算结果集中在一起绘图，便于对计算结果进行对比分析。数据流程如图1所示。

2.1　主函数的实现

该部分主要实现的是：叠加绘制NLNM和NHNM曲线、地震计周期和数采转换因子的赋值、功能函数的调用以及功率谱密度曲线的绘制。

图1　背景噪声计算程序流程图Fig.1　Flow of program of background noise calculation

(1) 确定全球高噪声模型NHNM和NLNM横坐标的对数值和纵坐标的DB值，用函数semilog()绘制全球噪声模型。

(2) 确定地震计的周期Ti和数采的转换因子Ci，调用功能函数Rsta=noisetest(sta，Ci，Ti)，进行参数传递，其中sta为台站名。

(3) 将功能函数的计算结果传递给semilog(x,y)函数绘制图形，其间还可调用subplot()函数以及colormap()函数进行分区绘图和色彩调整。

2.2　功能函数的实现

这部分主要实现的是创建功能函数Rsta=noisetest(sta，Ci，Ti)，该功能函数的运行结果是通过matlab的pwelch()函数计算所得到功率谱密度和频率。

(1) 用x=dlmread(strcat(sta，′.ch1′))函数读取指定文件名数据，strcat()函数的作用是进行字符串的拼接，用来读取批量的数据。

(2) 用a=x*Ci/2000将count值转化为速度值，其中Ci为数采转换因子，然后调用mean()函数，用公式a=a-mean(a)进行去均值处理。

(3) 确定welch算法的技术指标，matlab中用pwelch()函数来实现Welch算法。即：

[Pxx,f]=pwelch(x,window,noverlap,nfft,fs)，

式中：x为进行功率谱估计的有限长序列，即上述中的a；window指定采用的窗函数；noverlap指重叠点数；nfft设定FFT算法的长度；fs为采样频率；Pxx为输出的功率谱估计值；f为得到的频率点。

(4) 用公式Pxx=Pxx*4*pi^2*f^2将速度功率谱转换为加速度功率谱。然后由函数

sys=tf([100],[14*0.707*pi/Ti4*pi^2/Ti^2])，

求出频率响应H(s)，再由函数[mag phase]=bode(sys,2*pi*f)，求出幅值mag。最后，通过公式Pxx=Pxx/abs(mag^2)对加速度功率谱扣除仪器响应。上述函数tf()中，Ti为地震计周期。

实现功率谱计算的matlab代码：

Fsamp=100;%采样率。

x1=dlmread(strcat(Sta,′.ch1′));%读取台站数据。

a1=x1*C/2000；%count值转化为速度值。

a1=a1-mean(a1)；%去均值。

window=hamming(32768);%定义汉明窗和a1分段序列长度。

noverlap=20000;%重叠点数。

nfft=32768;%FFT算法长度。

[pxx1 f]=pwelch(a1,window,noverlap,nfft,Fsamp);%Welch算法。

pxx1=pxx1.*f.*f*4*pi*pi;%速度功率谱转加速度功率谱。

sys=tf([100],[12*0.707*pi/T/2 pi*pi/(T*T)/4]);%频率响应。

[mag phase]=bode(sys,2*pi*f);%求幅值mag。

mag1(1:nfft/2+1)=mag(1,1,1:nfft/2+1);%矩阵转置。

mag1=mag1.*mag1;

pxx1=pxx1./abs(mag1′);%去除仪器响应。

3应用实例

选择天津台网31个台站1小时的垂直向数据，时间尽量选在夜间2~4时，因为这个时段比较安静，可减少人为干扰，并且要保证24小时无地震且干扰不严重的波形。根据上述matlab实现的背景噪声计算程序，计算天津台网31个台站1小时垂直向的背景噪声功率谱，所得结果如图2所示。

图2　天津台网31个台站垂直向噪声功率谱密度曲线Fig.2　Power spectral densities of noise in vertical components of 31 stations in Tianjin Seismic Network

由计算结果可以看出，天津数字地震台网各台站的垂直向背景噪声功率谱密度曲线在0.125~1 Hz频段内，有一个明显的峰值，该峰值是一个具有小振幅、长周期特性的“单频峰”。另外，由于天津东部临近渤海海域，不可避免地会受到海洋的影响。因此，这个峰值可能是浅层地壳的海浪转化为地震能量，而引起的垂直向压力的变化或直接作用在海岸上所致[8]。在1~3 Hz频段内，大部分台站存在一个明显的峰值，且噪声水平与3~20 Hz频段内的高噪声水平相当，而在3~20 Hz频段内的差异较小，仅有一些分散频率的峰值，这可能是由人为活动引起的。因为周期小于1 s是人为噪声的卓越周期，人为噪声源大多数是分散或移动的，由来自各个方向的波叠加而形成的相当复杂且近似稳定的随机噪声场，它以高频面波的方式传播，随距离和深度的增加而迅速衰减。从整个区域的噪声水平来看，位于天津东南部的滨海新区、中南部的中心城区以及西南部的静海等地区，其噪声水平平均达到了-100 db，为噪声水平的高值区。因为这3个地区是天津经济和工业的发达地区，人口稠密，人类活动对台站背景噪声的影响更明显。尤其天津东南部滨海新区，是天津的沿海地区，受渤海的影响比较大，噪声水平较高。相对于这3个地区的高值，位于天津北部的蓟县，中北部的宝坻、宁河、武清等地区，是背景噪声的低值区。这些地区多为经济不发达地区，台站分布稀疏，受干扰较小。其中蓟县背景噪声值最低，达到-120 db，这是由于蓟县台是唯一的基岩观测台，本身比其他浮土覆盖层记录的背景噪声值低。另外，蓟县地震台位于山洞内，不会受到人为干扰的影响。天津地区背景噪声水平如图3所示。

图3　天津地区背景噪声水平Fig.3　Background noise level in Tianjin area

4认识与讨论

随着计算机和微电子技术的迅速发展，数字信号处理的理论、算法以及实现手段都得到了全面的提高。而matlab是一个强大的数值计算软件，程序设计自由度大，程序的可移植性好。作为国际控制界的标准计算软件，具有强大的分析和解决实际问题的能力，它简明的程序代码书写格式和丰富的命令式程序，可以有效降低程序编写难度。但也需要研究者对研究对象的物理模型、数学模型充分融合才能发挥其作用。

功率谱估计应用范围很广，日益受到各学科和应用领域的极大重视，是用有限长的数据来估计信号的功率谱，是数字信号处理的重要研究内容之一。该研究利用matlab，采用功率谱估计的Welch算法，设计出背景噪声计算程序，解决了天津测震台网背景噪声

的统一计算、批量成图、对比分析等问题，提高了工作效率。通过应用，软件还需进一步的完善，将背景噪声计算程序与JOPENS系统的流服务对接，实现背景噪声的实时计算，进而实现对数字地震台网的数据记录质量、观测系统的运行状态以及背景噪声源变化的准实时监控。

参考文献：

[1]Bomann P.(2012)[2012-09-24].New Manual of Seismological Observatory Practice[M/OL].GFZ:IASPEI.http://www.iaspei.org/.

[2]许可，刘瑞瑞，孔繁旭，等.基于绘制天津乡镇级地震速报分布图[J].地震地磁观测与研究，2014，35(1/2):276-279.

[3]王俊，徐戈，孙业君.江苏省区域地表背景噪声特性的分析[J].地震研究，2009，32(2):155-161.

[4]王建国，栗连弟，崔晓峰，等.数字化地震前兆台网日常工作管理软件[J].地震研究,2009，32(1):79-83.

[5]刘瑞丰，陈培善，党京平，等.宽频带数字地震记录仿真的应用[J].地震地磁观测与研究,1997，25(1):23-28.

[6]任枭，刘瑞丰，梁建宏，等.国家数字地震台网台站地动噪声功率谱分析[J].地震地磁观测与研究,2004，25(1):23-28.

[7]胡广书.数字信号处理—理论、算法与实现[M].北京:清华大学出版社，1997.

[8]Hasselmann K.A statistical analysis of the generation ofmicrosisms[J].Rev Geophys,1963(1):177-209.

Design and Application of Program of Background Noise Calculation Based on Matlab Software Platform

XU Ke, LIU Rui-rui, KONG Fan-xu, ZHU Hong, LIU Yan-li, LI Ye

(Earthquake Administration of Tianjin, Tianjin 300201, China)

Abstract:The program of background noise calculation using the algorithm of classical spectral estimation is developed on Matlab software platform. Design principle and application of the program are introduced in this paper. The application of the program solves the problems of background noise calculation of each station in Tianjin Seismic Network and mapping in batches and can monitor the data record quality, change of background noise and analyze noise character of each frequency band comparatively. So work efficiency of background noise calculation in Tianjin Seismic Network is improved.

Key words:Background Noise; Power spectral density; Seismic network of Tianjin

作者简介：第一许可(1982—)，男，辽宁省沈阳市人。2006年毕业于沈阳师范大学，工程师。

基金项目：天津市地震局科研课题(20141016)。

收稿日期：2014-10-31

中图分类号：P315.99

文献标志码：A

文章编号：1000-6265(2015)02-0042-04