文章快速检索    
  地震地磁观测与研究  2018, Vol. 39 Issue (2): 84-89  DOI: 10.3969/j.issn.1003-3246.2018.02.012
0

引用本文  

谢江涛, 林丽萍, 谌亮, 等. 地震台站台基噪声功率谱概率密度函数Matlab实现[J]. 地震地磁观测与研究, 2018, 39(2): 84-89. DOI: 10.3969/j.issn.1003-3246.2018.02.012.
Xie Jiangtao, Lin Liping, Chen Liang, et al. The program of probability density function of power spectral density curves from seismic noise of a station based on Matlab[J]. Seismological and Geomagnetic Observation and Research, 2018, 39(2): 84-89. DOI: 10.3969/j.issn.1003-3246.2018.02.012.

基金项目

国家科技支撑计划(项目编号:2014BAK03B04);四川省地震局科技专项(项目编号:LY1618,LY1810,LY1505)

作者简介

谢江涛(1986-), 男, 硕士研究生, 工程师, 主要从事数字地震观测、地震预警及背景噪声层析成像研究工作。E-mail:jiangtaoxie@outlook.com

文章历史

本文收到日期:2017-03-09
地震台站台基噪声功率谱概率密度函数Matlab实现
谢江涛 , 林丽萍 , 谌亮 , 赵敏     
中国成都 610041 四川省地震局
摘要:选取2015年四川数字测震台网中筠连和华蓥山地震台记录的垂直分向连续波形数据,利用Matlab软件,计算地震台站台基噪声功率谱概率密度函数,分析地震台站环境噪声特征。结果表明,台站环境噪声功率谱密度概率密度分布对地震事件波形(体波、面波)、人为噪声(台站周围人为活动、车辆及机器噪声等高频干扰)、系统瞬变(数据丢失、地震计小故障)以及仪器标定信号等反映较好。使用台基噪声功率谱概率密度函数方法,有利于监测地震台站数据记录,提高观测数据质量。
关键词台基噪声    仪器响应    功率谱密度    概率密度函数    Matlab    
The program of probability density function of power spectral density curves from seismic noise of a station based on Matlab
Xie Jiangtao, Lin Liping, Chen Liang, Zhao Min     
Sichuan Earthquake Agency, Chengdu 610041, China
Abstract: This paper introduces the design of the program for calculating the probability density function (PDF) of power spectral density (PSD) curves form seismic noise based on the Matlab software. And we computed the PDF of station JLI and HYS BHZ in Sichuan Digital Seismic Network, using PSDs during the period from January to December, 2015. The results manifest the region of the PDF is dominated by high power occurrences of naturally occurring earthquakes (body and surface waves), cultural nosie (source from actions of human being, traffic and machinery energy into the Earth), and recording system transients (data gaps and sensor glitches). And the calibrations are easily identifiable in the PDF. The PSD-PDF method can be helpful to monitor data of seismic stations and improve the quality of seismic data.
Key Words: ambient noise    instrument response    power spectrum density    probability density function    Matlab    
0 引言

地震信号通常是指,由一个天然或人工震源发出并经地下介质传播的瞬态波形,可以用于震源定位、孕震分析以及传播介质结构研究。地震台站的环境噪声水平决定了记录地震信号的能力,对地震噪声的量化是认识噪声水平的第一步。功率谱密度(Power Spectrum Density,PSD)是定量评价地震台站环境噪声水平的常规参数。Peterson(1993)研究全球范围地震台站的环境噪声功率谱密度,得到地球低噪声新模型(New Low Noise Model,NLNM)和地球高噪声新模型(New High Noise Model,NHNM),广泛应用于地震台站环境噪声水平评价。McNamara和Buland(2004)发展了Peterson的地球噪声模型估算方法,通过计算大量功率谱密度曲线的概率密度函数(Probability Density Function,PDF)分布,得到台站噪声水平最大概率分布模型和台网噪声低概率模型。概率密度函数(PDF)计算直接使用连续波形记录,地震数据并未筛选。因此,同时得到地震体波、面波信号、系统瞬变(如数据丢失、地震计小故障)及仪器(如调零、标定)干扰等概率分布。目前,概率密度函数方法被美国地质调查局国家地震信息中心(USGS National Earthquake Information Center)、IRIS数据管理中心以及新西兰地震台网用于地震台站背景噪声水平评价,也被用于美国(McNamara et al,2004)、意大利(Marzorati et al,2006)、新西兰(Rastin et al,2012)、中国华北地区(吴建平等,2012)及伦敦中部地区(Green et al,2017)的地震环境噪声特征分析。

1 方法原理 1.1 功率谱密度

通过对随机稳态的离散地震数据进行快速傅里叶变换,可估算地震噪声功率谱密度值(McNamara et al,2004吴建平等,2012)。对周期时间序列y(t)的有限范围傅里叶变换可表示为

$ Y\left( {f,{T_{\rm{r}}}} \right) = \int_0^{{T_{\rm{r}}}} {y\left( t \right){e^{ - i2\pi ft}}{\rm{d}}t} $ (1)

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

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

$ {Y_k} = \frac{{Y\left( {{f_k},{T_{\rm{r}}}} \right)}}{{\Delta t}} $ (2)

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

根据维纳—辛钦定理,功率谱密度(PSD)定义为

$ {P_k} = \frac{{2\Delta t}}{N}{\left| {{Y_k}} \right|^2} $ (3)

将地震台站单分向连续时间的地震波形数据分割为1 h长度的时间序列;对波形数据进行去均值、去长周期成分处理;将每小时时间序列分为13段,每段长度3600/13≈276.92 s,设采样率为100 Hz,即每段长度为27 692个采样点,每段之间重叠50%,为了提高傅里叶变换的计算速度,每个数据分段取最接近2的幂次方的点数,即每段长度为215=32 768个采样点;1 h的功率谱密度值由13段功率谱密度值进行平均得到。

速度型地震计记录的连续波形数据为地脉动速度量,将速度PSD值转换为加速度PSD值,采用以下公式

$ {P_{a.k}} = {\left( {2\pi {f_k}} \right)^2}{P_k} $ (4)

式中,Pk为速度功率谱密度,Pa.k为加速度功率谱密度。

地震计记录的地面运动可以表示为

$ y\left( t \right) = g\left( t \right)*h\left( t \right) $ (5)

其中,符号*表示褶积,g(t)为地震计能够观测的真实地面运动信息,h(t)为仪器响应。

在频率域表示为

$ Y\left( f \right) = G\left( f \right)H\left( f \right) $ (6)
$ Y\left( s \right) = G\left( s \right)H\left( s \right) $ (7)

以四川数字测震台网筠连地震台为例,台站采用港震机电技术有限公司生产的EDAS-24IP型数据采集器和Guralp公司生产的CMG-3ESPC型宽频带(60 s—50 Hz)地震计组成观测系统,幅频特性曲线见图 1,仪器响应信息RESP文件使用rdseed 5.3.1软件从.seed记录波形中得到,绘图使用JplotResp软件。在仪器频带范围及附近频带计算噪声功率谱时,是否扣除仪器传递函数随频率变化的影响,对PSD曲线总体形态影响不大(吴建平等,2012);超出频带范围越远,仪器自噪声等因素的影响越大,若未扣除仪器传递函数,PSD曲线将明显失真。

图 1 四川数字测震台网筠连地震台观测仪器(地震计和数据采集器)幅频特性曲线 Fig.1 The amplitude-frequency characteristics curve of instrument (seismic sensor and datalogger) of Junlian Seismic Station in Sichuan Digital Seismic Network

因此,利用地震记录波形数据计算地震噪声功率谱密度,需扣除仪器传递函数影响,以反映真实地面运动特征。

$ {\rm{PS}}{{\rm{D}}_a} = \frac{{{P_{a.k}}}}{{{{\left| {H\left( s \right){|^2}} \right|}_{s = i \cdot 2\pi f}}}} $ (8)

式中,PSDa为真实地面运动加速度功率谱,H (s)为系统传递函数。

1.2 概率密度函数

为了评估地震台站环境噪声的变化特征,通常对每小时的功率谱密度进行计算,产生大量光滑PSD曲线,然后计算台站地震噪声概率密度函数(PDF)(McNamara et al,2004吴建平等,2012)。为了对功率谱密度进行充分采样,按1⁄8倍频程间隔计算整个周期范围内各中心频率的平均值。功率谱密度值在短周期(高频率)Ts与长周期(低频率)Tl = 2Ts之间进行平均,对应的中心周期Tc为倍频程内的几何平均值。

$ {T_{\rm{c}}} = \sqrt {{T_{\rm{s}}} \times {T_{\rm{l}}}} $ (9)

Ts以1/8倍频程增加用于计算下一间隔的功率谱密度平均值,Ts = Ts×20.125。重新计算TlTs的值,并计算下一个中心周期TcTlTs周期范围内的功率谱密度平均值。

对于每一个给定的中心周期Tc,概率密度函数表示为

$ P\left( {{T_{\rm{c}}}} \right) = \frac{{{N_{PT{\rm{c}}}}}}{{{N_{T{\rm{c}}}}}} $ (10)

其中,NPTc是功率谱密度值落在某个1 dB间隔范围的数量,P的范围为-200—-40 dB,NTc是中心周期Tc在功率谱密度值范围内估计值总数。

2 Matlab实现

地震台站台基噪声功率谱概率密度函数程序基于Matlab软件设计,主要计算功率谱密度和概率密度函数2部分。利用Matlab强大的矩阵处理能力,计算单个台站1小时长度的时间序列PSD曲线,然后计算TsTl间的平均功率值,Ts按1/8倍频程增加,得到平滑PSD曲线,存入矩阵;读取矩阵中存储的平滑PSD曲线,根据式(10)计算某个1 dB区间的概率值,得到PSD值的概率密度分布图。

功率谱密度计算采用Welch(1967)方法,该方法是一种改进的周期图谱估计(periodogram spectrum estimating)方法和Bartlett法,对时间序列分段时允许每段数据有重叠,可以采用非矩形窗对数据进行加窗处理。使用Matlab程序的pwelch()函数实现,即[Pxx, F] = pwelch(X, WINDOW, NOVERLAP, NFFT, Fs),其中,X为进行功率谱估计的有限序列长度输入,WINDOW用于指定窗函数(boxcar,hamming,blackman),NOVERLAP为重叠点数,NFFT为FFT算法的设定长度,Fs为采样率,Pxx为功率谱估计输出值,F为所得频率点。

(1)实现功率谱密度计算,部分Matlab代码为

Sta.sampF = 100;%采样率Hz

BHZ = BHZ-mean(BHZ); %去均值,BHZ为1小时长度的随机离散地震数据

window = hamming(32768); %定义汉明窗及长度

noverlap = 16384; %重叠点数

nfft = 32768; %分段数据长度

[pbhz(:, k) f] = pwelch(BHZ, window, noverlap, nfft, Sta.sampF); %Welch算法计算功率谱密度

pbhz(:, k)=pbhz(:, k).*f.*f*4*pi*pi; %速度功率谱转换为加速度功率谱

sys=zpk(PZ.Z, PZ.P, PZ.K); %地震计和数据采集器总的传递函数

[mag phase] = bode(sys, 2*pi*f);

mag1(1:nfft/2+1)=mag(1, 1, 1:nfft/2+1);

mag1=mag1.*mag1;

pbhz(:, k)=pbhz(:, k)./abs(mag1'); %去除仪器响应

(2)实现PSD曲线平滑,部分Matlab代码为

Ts = 0.02; %定义起始周期

lengthF = length(f);

numP =1;

for n=1:100

    Tl = 2*Ts;

    Tc = sqrt(Ts*Tl);

    sum = 0.0;

    F(n)=1/Tc;

    numF=0;

    for l=1:lengthF

      if (1/Ts) > = f(l) & & f(l) > = (1/Tl)

        sum = pbhz(l, k)+sum;

        numF=numF+1;

      end

    end

    psdT(numP, k) = sum/numF; %计算Ts和Tl周期范围内的平均功率值,存储于psdT矩阵

    Ts=Ts*2^0.125; %Ts以1⁄8倍频程增加,用于计算下一个间隔的平均功率

    numP=numP+1;

end

(3)实现PSD概率密度函数分布计算,部分Matlab代码为

[pdflinenum pdfcolnum]=size(psdT);

for i=1:100

    k=1;

    for dB=-40:-1:-200

        Num.pdf(k, i) = length(find((10*log10(psdT(i, :)) < =dB) & (10*log10(psdT(i, :)) > dB-1))); %10*log10(psdT(i, :) %将计算的功率谱估计值转换为以dB为单位的数值

        pdfT(k, i) = Num.pdf(k, i)/pdfcolnum; %计算在某1 dB范围内的PSD概率值

        k=k+1;

    end

end

3 应用实例

选取四川数字测震台网的筠连、华蓥山地震台2015年记录的垂直分向连续波形数据[包含地震事件引起的地震波信号(体波、面波)及系统产生的瞬变信号(如数据丢失、地震计故障)],采用上述功率谱密度概率密度函数计算程序,分别得到8 760条和8 712条不同时间的PSD曲线,在此基础上获得台站环境噪声的概率密度函数分布,见图 2,不同频率PSD值日变化,见图 3

图 2 2015年筠连和华蓥山地震台垂直分向环境噪声功率谱密度概率密度函数(PDF)分布 (a)筠连台;(b)华蓥山 Fig.2 Ambient seismic noise PDF distribution obtained form continuous observation data in 2015 at BHZ of Junlian and Huayingshan seismic stations
图 3 2015年筠连和华蓥山地震台垂直分向环境噪声加速度PSD值日变化 (a)筠连台;(b)华蓥山 Fig.3 Diurnal variations in acceleration PSD for BHZ of Junlian and Huayingshan seismic stations from 2015

(1)筠连地震台台基噪声水平较低,低于NHNM和NLNM的平均值;高频段(> 1 Hz)噪声功率谱密度有2条较大的概率分布条带[图 2(a)],且昼夜变化明显[图 3(a)],白天高频噪声明显增加,主要由台站附近人为活动、车辆以及区域地震体波等引起;在周期1 s以上,台站环境噪声水平明显降低,第二类地脉动(double-frequency microseism)(约5 s)环境噪声PSD值变化较小[图 3(a)]。图 2(a)中-100—-60 dB范围内出现一条明显的浅粉色线,由每月地震计脉冲标定信号引起。

(2)华蓥山地震台台基噪声水平日变化特征与筠连台类似,高频段噪声水平随昼夜变化明显,白天高频噪声明显增加,噪声水平PSD值略高于NHNM和NLNM的平均值,见图 2(b)图 3(b);在周期60 s以上,超出仪器频带范围越远,仪器自噪声等其他干扰因素的影响越大,获得的结果就存在明显失真,见图 2(b),可见-100—-60 dB范围出现1条明显的浅粉色线,由每月地震计脉冲标定信号引起;-200—-150 dB范围对应1条明显的浅紫色线,与仪器故障有关。图 3中2个地震台高频段噪声水平开始增加的间略有不同,与台站所处区域位置有关(筠连台位于县城周边,华蓥山台位于乡镇附近),进一步说明高频段噪声与人为活动相关。

4 结束语

利用Matlab软件,实现地震台站台基噪声功率谱密度计算,得到不同时间的PSD曲线,按McNamara和Buland(2004)提出的概率密度函数方法,得到台站环境噪声PSD的概率密度分布。地震台站连续观测记录的功率谱密度概率分布,对人为活动和车辆等高频干扰及地震事件、仪器标定、地震计故障等反映较好,有利于监测地震台站数据记录质量。通过计算四川数字测震台网所辖地震台站连续记录的功率谱密度概率分布,得到最低概率分布噪声模型,比地球低噪声新模型(NLNM)更能真实反映该台网背景低噪声。

参考文献
吴建平, 欧阳飚, 王未来, 等. 华北地区地震环境噪声特征研究[J]. 地震学报, 2012, 34(6): 818-829.
Green D N, Bastow I D, Dashwood B, et al. Characterizing broadband seismic noise in Central London[J]. Seismological Research Letters, 2017, 88(1): 113-124. DOI:10.1785/0220160128
Marzorati S, Bindi D. Ambient noise levels in north central Italy[J]. Geochemistry, Geophysics, Geosystems, 2006, 7(9): 1-14.
McNamara D E, Buland R P. Ambient noise levels in the continental United States[J]. Bull Seismol Soc Am, 2004, 94(4): 1517-1527. DOI:10.1785/012003001
Peterson J. Observations and modeling of seismic background noise[J]. USGS Open File Report, 1993: 93-322.
Rastin S J, Unsworth C P, Gledhill K R, et al. A detailed noise characterization and sensor evaluation of the North Island of New Zealand using the PQLX data quality control system[J]. Bull Seismol Soc Am, 2012, 102(1): 98-113. DOI:10.1785/0120110064
Welch P. The use of fast Fourier transform for the estimation of power spectra:a method based on time averaging over short, modified periodograms[J]. IEEE Transactions on audio and electroacoustics, 1967, 15(2): 70-73. DOI:10.1109/TAU.1967.1161901