許　可，劉瑞瑞，孔繁旭，朱　宏，柳艷麗，李　曄

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

基于matlab的背景噪聲計(jì)算程序的設(shè)計(jì)與應(yīng)用

許可，劉瑞瑞，孔繁旭，朱宏，柳艷麗，李曄

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

摘要：文章介紹了利用軟件平臺(tái)，開發(fā)出采用經(jīng)典譜估計(jì)中算法的背景噪聲計(jì)算程序的設(shè)計(jì)原理及實(shí)際應(yīng)用。該程序的應(yīng)用，解決了天津測(cè)震臺(tái)網(wǎng)各臺(tái)站背景噪聲的統(tǒng)一計(jì)算、批量成圖等問題，從而實(shí)現(xiàn)了對(duì)數(shù)字地震臺(tái)網(wǎng)數(shù)據(jù)記錄質(zhì)量、背景噪聲變化的監(jiān)控以及各頻段噪聲特性的對(duì)比分析，提高了天津臺(tái)網(wǎng)背景噪聲計(jì)算的工作效率。

關(guān)鍵詞：背景噪聲；功率譜密度；天津測(cè)震臺(tái)網(wǎng)

0引言

在數(shù)字化觀測(cè)系統(tǒng)中，系統(tǒng)的瞬態(tài)變化、儀器毛刺、數(shù)據(jù)記錄的階躍和尖峰以及由于系統(tǒng)故障引起的信號(hào)失真和重大環(huán)境干擾源的出現(xiàn)等，都是影響數(shù)據(jù)記錄質(zhì)量的重要因素[1]。隨著信息化社會(huì)的到來，人們的認(rèn)識(shí)水平不斷提高，逐漸關(guān)注實(shí)時(shí)記錄的背景噪聲，希望它能提供科學(xué)有用的信息[2]。在數(shù)字地震觀測(cè)系統(tǒng)記錄中，除真實(shí)的地震信號(hào)外，還包含了由于系統(tǒng)故障引起的各種失真信號(hào)以及臺(tái)站周圍的環(huán)境噪聲干擾等信息[3]。這些干擾信息嚴(yán)重影響觀測(cè)數(shù)據(jù)的質(zhì)量，需要加以區(qū)分以及排除，確保觀測(cè)數(shù)據(jù)質(zhì)量[4]。因此，如何及時(shí)了解數(shù)字地震觀測(cè)系統(tǒng)運(yùn)行的狀態(tài)，保證地震臺(tái)網(wǎng)記錄的數(shù)據(jù)質(zhì)量以及識(shí)別干擾源，就成為當(dāng)前數(shù)字化地震臺(tái)網(wǎng)系統(tǒng)監(jiān)控技術(shù)中的重要研究課題之一。

“十五”項(xiàng)目后，天津測(cè)震臺(tái)網(wǎng)建成1個(gè)臺(tái)網(wǎng)中心和31個(gè)地震臺(tái)站。臺(tái)網(wǎng)孔徑東西約127 km，南北約160 km，地震臺(tái)站平均臺(tái)間距為30.6 km。臺(tái)網(wǎng)中心采用JOPENS對(duì)數(shù)據(jù)進(jìn)行接收、交換存儲(chǔ)及處理。在經(jīng)過“首都圈測(cè)震臺(tái)網(wǎng)的升級(jí)改造”項(xiàng)目后，增添多種不同的儀器類型，導(dǎo)致地震臺(tái)站的背景噪聲計(jì)算步驟繁雜，需花費(fèi)大量的時(shí)間。為此，本文利用matlab平臺(tái)研發(fā)了運(yùn)用經(jīng)典譜估計(jì)中的Welch算法，解決天津測(cè)震臺(tái)網(wǎng)31個(gè)臺(tái)站背景噪聲的統(tǒng)一計(jì)算、批量成圖等問題，從而實(shí)現(xiàn)對(duì)數(shù)字地震臺(tái)網(wǎng)數(shù)據(jù)記錄質(zhì)量、背景噪聲變化的監(jiān)控以及各頻段噪聲特性的對(duì)比分析。

1計(jì)算原理

1.1　去直流偏移

在數(shù)字地震記錄中會(huì)出現(xiàn)直流偏移，這個(gè)信號(hào)不表示真實(shí)的地面運(yùn)動(dòng)，必須消除這一直流分量[5]。在實(shí)際計(jì)算過程中，把整個(gè)記錄長(zhǎng)度求和除以記錄總數(shù)，取記錄的平均值作為記錄的直流偏離量。然后，逐點(diǎn)減掉直流偏離量作為第一步的校正過程，公式為：

式中：fi表示每個(gè)采樣點(diǎn)的值；N為記錄長(zhǎng)度的總采樣個(gè)數(shù)。

1.2　扣除儀器響應(yīng)

數(shù)字地震觀測(cè)系統(tǒng)是由地震計(jì)拾取地面運(yùn)動(dòng)，數(shù)據(jù)采集器對(duì)地震計(jì)拾取的信號(hào)進(jìn)行數(shù)字化采集，通過數(shù)據(jù)專線傳輸?shù)脚_(tái)網(wǎng)中心的數(shù)據(jù)接收系統(tǒng)。地動(dòng)速度變換成數(shù)字記錄的過程就是儀器響應(yīng)的過程。在實(shí)際計(jì)算噪聲功率譜的時(shí)候，要考慮觀測(cè)儀器對(duì)地動(dòng)噪聲的影響，要想得到地動(dòng)噪聲的絕對(duì)量值，還應(yīng)從計(jì)算結(jié)果中扣除儀器的影響，需利用數(shù)字地震儀器的傳遞函數(shù)來扣除儀器響應(yīng)。數(shù)字地震儀器的傳遞函數(shù)[6]為：

式中：Sd代表儀器的靈敏度；A0是歸一化常數(shù)；r1為復(fù)零點(diǎn)；Pm為復(fù)極點(diǎn)；L為零點(diǎn)個(gè)數(shù)；M為極點(diǎn)個(gè)數(shù)。

1.3　功率譜密度的計(jì)算

經(jīng)典譜估計(jì)方法在工程中都是以離散傅立葉變換為基礎(chǔ)的，它隱含著對(duì)無限長(zhǎng)數(shù)據(jù)加窗處理，所以，經(jīng)典譜估計(jì)有著分辨率不高、能量泄漏的固有缺點(diǎn)。為克服這些缺點(diǎn)，經(jīng)過研究，提出了各種算法。目前，經(jīng)典譜估計(jì)算法有周期圖法、Barlett算法、Welch算法、Nattall算法等。Welch算法是由welch提出的對(duì)周期圖的修正算法，是經(jīng)典譜估計(jì)中獲得有效應(yīng)用的一種算法。

功率譜密度的計(jì)算采用經(jīng)典譜估計(jì)中的Welch算法，Welch算法是在周期圖法的基礎(chǔ)上改進(jìn)得到的。Welch算法譜估計(jì)采取數(shù)據(jù)分段加窗處理再求平均的辦法，先分別求出每段的譜估計(jì)，然后進(jìn)行總平均。即：把某一長(zhǎng)度為N的數(shù)據(jù)x(n)分成L段，每一段為M(在分段時(shí)可允許每一段的數(shù)據(jù)有部分的重疊)，分別求每段的功率譜，然后加以平均。每一段的數(shù)據(jù)窗口可使用漢寧窗或哈明窗。這樣可以改善由于矩形窗邊瓣較大所產(chǎn)生的譜失真。用p(ω)表示用Welch算法估計(jì)的功率譜[7]。即：

式中：U為歸一化因子；d(n)是數(shù)據(jù)窗。

2程序的matlab實(shí)現(xiàn)

matlab是一套用于科學(xué)計(jì)算的可視化、高性能的語言和軟件環(huán)境，它集數(shù)值分析、矩陣運(yùn)算、信號(hào)處理和圖形顯示于一體，已經(jīng)成為數(shù)字信號(hào)處理應(yīng)用中分析和仿真的主要工具，工具箱中包含各種經(jīng)典和現(xiàn)代的數(shù)字信號(hào)處理技術(shù)，能實(shí)現(xiàn)各種數(shù)字信號(hào)的處理。

本程序基于matlab設(shè)計(jì)，主要包含兩部分。一部分為程序的主函數(shù)，另一部分為程序的功能函數(shù)。程序的主要思想是先通過功能函數(shù)計(jì)算出功率譜，再通過主函數(shù)初始化參數(shù)調(diào)用功能函數(shù)進(jìn)行繪圖操作，實(shí)現(xiàn)統(tǒng)一計(jì)算、批量成圖等功能。功能函數(shù)只負(fù)責(zé)計(jì)算出功率譜密度，具體的參數(shù)設(shè)定、繪圖操作等在主函數(shù)中完成，這樣既可以對(duì)計(jì)算結(jié)果分別成圖，又可以對(duì)計(jì)算結(jié)果集中成圖，或者可將所關(guān)心的某幾個(gè)計(jì)算結(jié)果集中在一起繪圖，便于對(duì)計(jì)算結(jié)果進(jìn)行對(duì)比分析。數(shù)據(jù)流程如圖1所示。

2.1　主函數(shù)的實(shí)現(xiàn)

該部分主要實(shí)現(xiàn)的是：疊加繪制NLNM和NHNM曲線、地震計(jì)周期和數(shù)采轉(zhuǎn)換因子的賦值、功能函數(shù)的調(diào)用以及功率譜密度曲線的繪制。

圖1　背景噪聲計(jì)算程序流程圖Fig.1　Flow of program of background noise calculation

(1) 確定全球高噪聲模型NHNM和NLNM橫坐標(biāo)的對(duì)數(shù)值和縱坐標(biāo)的DB值，用函數(shù)semilog()繪制全球噪聲模型。

(2) 確定地震計(jì)的周期Ti和數(shù)采的轉(zhuǎn)換因子Ci，調(diào)用功能函數(shù)Rsta=noisetest(sta，Ci，Ti)，進(jìn)行參數(shù)傳遞，其中sta為臺(tái)站名。

(3) 將功能函數(shù)的計(jì)算結(jié)果傳遞給semilog(x,y)函數(shù)繪制圖形，其間還可調(diào)用subplot()函數(shù)以及colormap()函數(shù)進(jìn)行分區(qū)繪圖和色彩調(diào)整。

2.2　功能函數(shù)的實(shí)現(xiàn)

這部分主要實(shí)現(xiàn)的是創(chuàng)建功能函數(shù)Rsta=noisetest(sta，Ci，Ti)，該功能函數(shù)的運(yùn)行結(jié)果是通過matlab的pwelch()函數(shù)計(jì)算所得到功率譜密度和頻率。

(1) 用x=dlmread(strcat(sta，′.ch1′))函數(shù)讀取指定文件名數(shù)據(jù)，strcat()函數(shù)的作用是進(jìn)行字符串的拼接，用來讀取批量的數(shù)據(jù)。

(2) 用a=x*Ci/2000將count值轉(zhuǎn)化為速度值，其中Ci為數(shù)采轉(zhuǎn)換因子，然后調(diào)用mean()函數(shù)，用公式a=a-mean(a)進(jìn)行去均值處理。

(3) 確定welch算法的技術(shù)指標(biāo)，matlab中用pwelch()函數(shù)來實(shí)現(xiàn)Welch算法。即：

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

式中：x為進(jìn)行功率譜估計(jì)的有限長(zhǎng)序列，即上述中的a；window指定采用的窗函數(shù)；noverlap指重疊點(diǎn)數(shù)；nfft設(shè)定FFT算法的長(zhǎng)度；fs為采樣頻率；Pxx為輸出的功率譜估計(jì)值；f為得到的頻率點(diǎn)。

(4) 用公式Pxx=Pxx*4*pi^2*f^2將速度功率譜轉(zhuǎn)換為加速度功率譜。然后由函數(shù)

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

求出頻率響應(yīng)H(s)，再由函數(shù)[mag phase]=bode(sys,2*pi*f)，求出幅值mag。最后，通過公式Pxx=Pxx/abs(mag^2)對(duì)加速度功率譜扣除儀器響應(yīng)。上述函數(shù)tf()中，Ti為地震計(jì)周期。

實(shí)現(xiàn)功率譜計(jì)算的matlab代碼：

Fsamp=100;%采樣率。

x1=dlmread(strcat(Sta,′.ch1′));%讀取臺(tái)站數(shù)據(jù)。

a1=x1*C/2000；%count值轉(zhuǎn)化為速度值。

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

window=hamming(32768);%定義漢明窗和a1分段序列長(zhǎng)度。

noverlap=20000;%重疊點(diǎn)數(shù)。

nfft=32768;%FFT算法長(zhǎng)度。

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

pxx1=pxx1.*f.*f*4*pi*pi;%速度功率譜轉(zhuǎn)加速度功率譜。

sys=tf([100],[12*0.707*pi/T/2 pi*pi/(T*T)/4]);%頻率響應(yīng)。

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

mag1(1:nfft/2+1)=mag(1,1,1:nfft/2+1);%矩陣轉(zhuǎn)置。

mag1=mag1.*mag1;

pxx1=pxx1./abs(mag1′);%去除儀器響應(yīng)。

3應(yīng)用實(shí)例

選擇天津臺(tái)網(wǎng)31個(gè)臺(tái)站1小時(shí)的垂直向數(shù)據(jù)，時(shí)間盡量選在夜間2~4時(shí)，因?yàn)檫@個(gè)時(shí)段比較安靜，可減少人為干擾，并且要保證24小時(shí)無地震且干擾不嚴(yán)重的波形。根據(jù)上述matlab實(shí)現(xiàn)的背景噪聲計(jì)算程序，計(jì)算天津臺(tái)網(wǎng)31個(gè)臺(tái)站1小時(shí)垂直向的背景噪聲功率譜，所得結(jié)果如圖2所示。

圖2　天津臺(tái)網(wǎng)31個(gè)臺(tái)站垂直向噪聲功率譜密度曲線Fig.2　Power spectral densities of noise in vertical components of 31 stations in Tianjin Seismic Network

由計(jì)算結(jié)果可以看出，天津數(shù)字地震臺(tái)網(wǎng)各臺(tái)站的垂直向背景噪聲功率譜密度曲線在0.125~1 Hz頻段內(nèi)，有一個(gè)明顯的峰值，該峰值是一個(gè)具有小振幅、長(zhǎng)周期特性的“單頻峰”。另外，由于天津東部臨近渤海海域，不可避免地會(huì)受到海洋的影響。因此，這個(gè)峰值可能是淺層地殼的海浪轉(zhuǎn)化為地震能量，而引起的垂直向壓力的變化或直接作用在海岸上所致[8]。在1~3 Hz頻段內(nèi)，大部分臺(tái)站存在一個(gè)明顯的峰值，且噪聲水平與3~20 Hz頻段內(nèi)的高噪聲水平相當(dāng)，而在3~20 Hz頻段內(nèi)的差異較小，僅有一些分散頻率的峰值，這可能是由人為活動(dòng)引起的。因?yàn)橹芷谛∮? s是人為噪聲的卓越周期，人為噪聲源大多數(shù)是分散或移動(dòng)的，由來自各個(gè)方向的波疊加而形成的相當(dāng)復(fù)雜且近似穩(wěn)定的隨機(jī)噪聲場(chǎng)，它以高頻面波的方式傳播，隨距離和深度的增加而迅速衰減。從整個(gè)區(qū)域的噪聲水平來看，位于天津東南部的濱海新區(qū)、中南部的中心城區(qū)以及西南部的靜海等地區(qū)，其噪聲水平平均達(dá)到了-100 db，為噪聲水平的高值區(qū)。因?yàn)檫@3個(gè)地區(qū)是天津經(jīng)濟(jì)和工業(yè)的發(fā)達(dá)地區(qū)，人口稠密，人類活動(dòng)對(duì)臺(tái)站背景噪聲的影響更明顯。尤其天津東南部濱海新區(qū)，是天津的沿海地區(qū)，受渤海的影響比較大，噪聲水平較高。相對(duì)于這3個(gè)地區(qū)的高值，位于天津北部的薊縣，中北部的寶坻、寧河、武清等地區(qū)，是背景噪聲的低值區(qū)。這些地區(qū)多為經(jīng)濟(jì)不發(fā)達(dá)地區(qū)，臺(tái)站分布稀疏，受干擾較小。其中薊縣背景噪聲值最低，達(dá)到-120 db，這是由于薊縣臺(tái)是唯一的基巖觀測(cè)臺(tái)，本身比其他浮土覆蓋層記錄的背景噪聲值低。另外，薊縣地震臺(tái)位于山洞內(nèi)，不會(huì)受到人為干擾的影響。天津地區(qū)背景噪聲水平如圖3所示。

圖3　天津地區(qū)背景噪聲水平Fig.3　Background noise level in Tianjin area

4認(rèn)識(shí)與討論

隨著計(jì)算機(jī)和微電子技術(shù)的迅速發(fā)展，數(shù)字信號(hào)處理的理論、算法以及實(shí)現(xiàn)手段都得到了全面的提高。而matlab是一個(gè)強(qiáng)大的數(shù)值計(jì)算軟件，程序設(shè)計(jì)自由度大，程序的可移植性好。作為國際控制界的標(biāo)準(zhǔn)計(jì)算軟件，具有強(qiáng)大的分析和解決實(shí)際問題的能力，它簡(jiǎn)明的程序代碼書寫格式和豐富的命令式程序，可以有效降低程序編寫難度。但也需要研究者對(duì)研究對(duì)象的物理模型、數(shù)學(xué)模型充分融合才能發(fā)揮其作用。

功率譜估計(jì)應(yīng)用范圍很廣，日益受到各學(xué)科和應(yīng)用領(lǐng)域的極大重視，是用有限長(zhǎng)的數(shù)據(jù)來估計(jì)信號(hào)的功率譜，是數(shù)字信號(hào)處理的重要研究?jī)?nèi)容之一。該研究利用matlab，采用功率譜估計(jì)的Welch算法，設(shè)計(jì)出背景噪聲計(jì)算程序，解決了天津測(cè)震臺(tái)網(wǎng)背景噪聲

的統(tǒng)一計(jì)算、批量成圖、對(duì)比分析等問題，提高了工作效率。通過應(yīng)用，軟件還需進(jìn)一步的完善，將背景噪聲計(jì)算程序與JOPENS系統(tǒng)的流服務(wù)對(duì)接，實(shí)現(xiàn)背景噪聲的實(shí)時(shí)計(jì)算，進(jìn)而實(shí)現(xiàn)對(duì)數(shù)字地震臺(tái)網(wǎng)的數(shù)據(jù)記錄質(zhì)量、觀測(cè)系統(tǒng)的運(yùn)行狀態(tài)以及背景噪聲源變化的準(zhǔn)實(shí)時(shí)監(jiān)控。

參考文獻(xiàn)：

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

[2]許可，劉瑞瑞，孔繁旭，等.基于繪制天津鄉(xiāng)鎮(zhèn)級(jí)地震速報(bào)分布圖[J].地震地磁觀測(cè)與研究，2014，35(1/2):276-279.

[3]王俊，徐戈，孫業(yè)君.江蘇省區(qū)域地表背景噪聲特性的分析[J].地震研究，2009，32(2):155-161.

[4]王建國，栗連弟，崔曉峰，等.數(shù)字化地震前兆臺(tái)網(wǎng)日常工作管理軟件[J].地震研究,2009，32(1):79-83.

[5]劉瑞豐，陳培善，黨京平，等.寬頻帶數(shù)字地震記錄仿真的應(yīng)用[J].地震地磁觀測(cè)與研究,1997，25(1):23-28.

[6]任梟，劉瑞豐，梁建宏，等.國家數(shù)字地震臺(tái)網(wǎng)臺(tái)站地動(dòng)噪聲功率譜分析[J].地震地磁觀測(cè)與研究,2004，25(1):23-28.

[7]胡廣書.數(shù)字信號(hào)處理—理論、算法與實(shí)現(xiàn)[M].北京:清華大學(xué)出版社，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

作者簡(jiǎn)介：第一許可(1982—)，男，遼寧省沈陽市人。2006年畢業(yè)于沈陽師范大學(xué)，工程師。

基金項(xiàng)目：天津市地震局科研課題(20141016)。

收稿日期：2014-10-31

中圖分類號(hào)：P315.99

文獻(xiàn)標(biāo)志碼：A

文章編號(hào)：1000-6265(2015)02-0042-04