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

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

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

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

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

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

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

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

0引言

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

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

1計算原理

1.1　去直流偏移

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

式中：fi表示每個采樣點的值；N為記錄長度的總采樣個數(shù)。

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

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

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

1.3　功率譜密度的計算

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(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)對加速度功率譜扣除儀器響應(yīng)。上述函數(shù)tf()中，Ti為地震計周期。

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

Fsamp=100;%采樣率。

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

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

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

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

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

nfft=32768;%FFT算法長度。

[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)用實例

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

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

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

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

4認識與討論

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

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

的統(tǒng)一計算、批量成圖、對比分析等問題，提高了工作效率。通過應(yīng)用，軟件還需進一步的完善，將背景噪聲計算程序與JOPENS系統(tǒng)的流服務(wù)對接，實現(xiàn)背景噪聲的實時計算，進而實現(xiàn)對數(shù)字地震臺網(wǎng)的數(shù)據(jù)記錄質(zhì)量、觀測系統(tǒng)的運行狀態(tài)以及背景噪聲源變化的準實時監(jiā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].地震地磁觀測與研究，2014，35(1/2):276-279.

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

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

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

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

[7]胡廣書.數(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

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

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

收稿日期：2014-10-31

中圖分類號：P315.99

文獻標志碼：A

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