林彬華 金星 廖詩榮 李軍 黃玲珠 朱耿青

1)福建省地震局，福州市華鴻路7號(hào) 350003

2)福州大學(xué)，福州市閩侯縣大學(xué)城學(xué)園路2號(hào) 350002

3)中國(guó)地震局工程力學(xué)研究所，哈爾濱 150080

0 引言

地震波形數(shù)據(jù)的質(zhì)量很大程度取決于地震觀測(cè)儀器系統(tǒng)是否正常工作。為了檢測(cè)地震儀器是否發(fā)生異常，通過每天對(duì)其進(jìn)行脈沖標(biāo)定，可以發(fā)現(xiàn)大部分故障。但是這種方法需要人工操作判斷，工作量大、實(shí)時(shí)性差，很難及時(shí)發(fā)現(xiàn)異常。除此之外，對(duì)于要求分秒必爭(zhēng)的地震預(yù)警，要實(shí)時(shí)監(jiān)測(cè)數(shù)據(jù)質(zhì)量狀況，為后續(xù)到來的地震，快速計(jì)算震級(jí)和定位提供客觀、科學(xué)的數(shù)據(jù)。因此目前迫切需要一種新手段來實(shí)時(shí)監(jiān)測(cè)地震記錄系統(tǒng)是否能正常工作，從而保證地震數(shù)據(jù)的最優(yōu)質(zhì)量。

國(guó)際上，McNamara等(2004、2005a)提出用功率譜概率密度函數(shù)(PDF)方法進(jìn)行臺(tái)站地震噪聲水平監(jiān)測(cè)，并在GSN與ANSS等臺(tái)網(wǎng)應(yīng)用于日常儀器工作狀態(tài)檢測(cè)。廖詩榮等(2008)選取臺(tái)站48小時(shí)的噪聲記錄，并以160s為單位對(duì)噪聲記錄進(jìn)行分段，然后采用概率密度函數(shù)(PDF)方法求出臺(tái)址的總體噪聲水平，最后完成了自動(dòng)化處理地震臺(tái)站勘選測(cè)試。Sleeman(2007)在Orfeus數(shù)據(jù)中心通過監(jiān)測(cè)實(shí)時(shí)波形數(shù)據(jù)中地震噪聲加速度功率譜密度(PSD)值的變化，在VEBSN(Virtual European Broadband Seismograph Network)實(shí)現(xiàn)了對(duì)寬頻帶地震臺(tái)網(wǎng)波形數(shù)據(jù)質(zhì)量的自動(dòng)檢測(cè)。徐嘉雋等(2010)研發(fā)了數(shù)據(jù)質(zhì)量檢測(cè)軟件，其主要原理也是利用概率密度函數(shù)方法計(jì)算一天的噪聲記錄，然后針對(duì)自動(dòng)繪制出的PDF圖及RMS值，直觀地判斷觀測(cè)數(shù)據(jù)的質(zhì)量，并已經(jīng)應(yīng)用于日常地震觀測(cè)系統(tǒng)數(shù)據(jù)質(zhì)量的檢測(cè)。為了能夠?qū)崟r(shí)監(jiān)測(cè)地震儀器工作狀態(tài)，本文采用功率譜概率密度函數(shù)(PDF)方法繪制出臺(tái)站的PDF圖，再用網(wǎng)格概率方法確定臺(tái)站高低噪聲參照線，進(jìn)而研究噪聲異常的遴選方法，最后形成地震噪聲實(shí)時(shí)監(jiān)測(cè)系統(tǒng)。

1 地震噪聲概率密度函數(shù)方法概述

可將地震噪聲視為平穩(wěn)隨機(jī)信號(hào)，依據(jù)隨機(jī)過程理論，通常用概率統(tǒng)計(jì)方法來描述隨機(jī)信號(hào)。平穩(wěn)隨機(jī)信號(hào)在時(shí)間上是無限的，其能量也是無限的，但其功率卻是有限的。平穩(wěn)隨機(jī)信號(hào)的功率譜反映了信號(hào)的功率在頻域內(nèi)隨頻率的分布，也稱功率譜密度(簡(jiǎn)稱PSD)。

Peterson等(1993)通過分析全球75個(gè)地震臺(tái)站近2000條噪聲記錄的功率譜密度分布，給出了全球低噪聲模型(NLNM)與高噪聲模型(NHNM)，該模型目前廣泛應(yīng)用于對(duì)臺(tái)址噪聲水平的評(píng)價(jià)、儀器標(biāo)準(zhǔn)定義以及不同背景噪聲水平下地震計(jì)響應(yīng)的預(yù)測(cè)等。需要指出的是，Peterson用來進(jìn)行建立地震噪聲模型的2000條噪聲記錄是從近12000條波形記錄中篩選出來的，其它近10000條包含地震、爆破、標(biāo)定等非地震噪聲的記錄段并未參與最后的處理。

為了減免非地震噪聲記錄段的篩選與截取這個(gè)波形預(yù)處理環(huán)節(jié)，對(duì)包括非地震噪聲記錄段在內(nèi)的原始連續(xù)波形記錄直接處理，McNamara等(2005a，2005b)提出應(yīng)用PSD概率密度函數(shù)(PDF)方法進(jìn)行地震噪聲PSD值的計(jì)算。該方法的主要思路是:將原始波形數(shù)據(jù)分成n個(gè)記錄段，采用與Peterson相同的方法對(duì)每個(gè)記錄段計(jì)算PSD值，使用1/3倍頻程的頻率間隔對(duì)每個(gè)記錄段PSD曲線進(jìn)行平滑;然后計(jì)算PSD值落在某一個(gè)頻點(diǎn)某一功率窗內(nèi)的記錄段數(shù)目，以該記錄段數(shù)目與總記錄段數(shù)目n的比值作為該頻段該功率窗的PSD概率密度函數(shù)的取值。圖1(a)是對(duì)YDXS臺(tái)站194545個(gè)原始噪聲波形記錄段(1年數(shù)據(jù))使用概率密度函數(shù)方法處理得到的PDF圖。在該圖中地震噪聲以高概率取值的形式出現(xiàn)，并未因人工截除的地震體波與面波、儀器尖脈沖干擾、脈沖標(biāo)定等而出現(xiàn)低概率值形式，不會(huì)影響對(duì)高概率水平下環(huán)境地震噪聲水平的評(píng)估。通過PDF圖還可以檢測(cè)記錄系統(tǒng)的故障和全面評(píng)判臺(tái)站的數(shù)據(jù)質(zhì)量。

2 福建地區(qū)地震噪聲觀測(cè)資料與處理

2.1 地震噪聲資料收集與預(yù)處理

福建省測(cè)震臺(tái)網(wǎng)由原先的41個(gè)測(cè)震臺(tái)站和《十一五》新建的44個(gè)測(cè)震臺(tái)站組成，圖2給出了福建省這85個(gè)測(cè)震臺(tái)站的分布。本文所研究的地震噪聲資料是福建省85個(gè)測(cè)震臺(tái)站2012年1～12月噪聲的垂直向記錄，采樣頻率為100Hz(李軍等，2011)，其中大部分儀器所記錄的頻段在1/60～50Hz之間?？紤]到頻段前后會(huì)出現(xiàn)一定的衰減，所以本文設(shè)定PSD值的頻率取值范圍為1/50～40Hz。

圖1 功率譜概率密度函數(shù)(PDF)分布圖。 (a)利用臺(tái)站1年的噪聲數(shù)據(jù)計(jì)算出來的PDF圖，點(diǎn)線為確定的高、低噪聲參照線;(b)對(duì)應(yīng)三維PDF圖中1Hz頻點(diǎn)處的功率譜密度概率分布圖

福建臺(tái)網(wǎng)的地震記錄通常以EVT格式或SEED格式存儲(chǔ)，數(shù)值的單位為counts。使用時(shí)都需要對(duì)數(shù)據(jù)做如下處理，其中地動(dòng)速度(李軍，2007)的計(jì)算公式為

對(duì)一些漂移的波形采用最小二乘法作基線校正，具體公式為

式中v(t)為未校正的地面速度記錄，v*(t)為校正后的地面速度記錄，參數(shù)a、b是對(duì)每條記錄v(t)進(jìn)行最小二乘擬合求得的，其計(jì)算公式為

2.2 地震噪聲功率譜計(jì)算過程

利用噪聲記錄計(jì)算噪聲功率譜的數(shù)據(jù)處理過程主要步驟如下:

(1)記錄段的選取。對(duì)1年期的噪聲記錄以每5分鐘(300s)進(jìn)行分段，且為了避免PSD值的變化過大，對(duì)每段記錄間按50%的疊加率選取。記錄段長(zhǎng)度Tr取決于我們感興趣信號(hào)的最長(zhǎng)周期TL，通常Tr的取值需達(dá)到TL的6倍以上(Bormann，2002)。因此取記錄段長(zhǎng)度300s，可以滿足反映50s低頻地震噪聲PSD值的需要。

(2)速度功率譜密度值估算。最常用的估算地震噪聲PSD值的方法是直接傅里葉變換法。該方法通過計(jì)算有限長(zhǎng)度數(shù)據(jù)序列的FFT變換來計(jì)算PSD值。為了得到以頻率為自變量的FFT變換，可以將數(shù)據(jù)序列表示為

式中，f=k/NΔt(k=0，1，…，N-1);N為采樣點(diǎn)個(gè)數(shù);Δt為采樣時(shí)間間隔。y(nΔt)是每段原始噪聲記錄在時(shí)域中的數(shù)據(jù)系列，Y(f)是傅里葉變換得到的振幅譜。

功率譜反映了信號(hào)的能量隨頻率的分布情況，功率譜與傅里葉振幅譜二者在本質(zhì)上并無多大不同，只是功率譜的縱坐標(biāo)大體上相當(dāng)于傅里葉振幅譜縱坐標(biāo)的平方，其表達(dá)式為

(3)加速度功率譜密度的計(jì)算。當(dāng)頻率為f時(shí)，速度功率譜密度PSDv與加速度功率譜密度PSDa的轉(zhuǎn)換公式為

其中:PSDv的單位為(m/s)2/Hz;PSDa的單位為(m/s2)2/Hz。

項(xiàng)目建設(shè)高度重視資源整合。首先，在項(xiàng)目設(shè)計(jì)階段，將水文系統(tǒng)現(xiàn)有水文站網(wǎng)、中小河流監(jiān)測(cè)建設(shè)及規(guī)劃站網(wǎng)納入山洪災(zāi)害防治非工程措施暴雨洪水監(jiān)測(cè)站網(wǎng)，充分利用了現(xiàn)有監(jiān)測(cè)站網(wǎng)資源，避免重復(fù)建設(shè)。項(xiàng)目實(shí)施過程中盡可能避免與氣象站重合。二是充分利用防汛抗旱指揮系統(tǒng)建設(shè)成果，項(xiàng)目建設(shè)要求在完全滿足國(guó)家防總對(duì)山洪災(zāi)害非工程措施建設(shè)技術(shù)要求的同時(shí)，從監(jiān)測(cè)數(shù)據(jù)流程規(guī)劃、數(shù)據(jù)采集、傳輸和存儲(chǔ)的技術(shù)標(biāo)準(zhǔn)等均保持了與防汛抗旱指揮系統(tǒng)標(biāo)準(zhǔn)一致。三是山洪災(zāi)害預(yù)警平臺(tái)與防汛抗旱指揮系統(tǒng)進(jìn)行了完整的融合，山洪預(yù)警作為防汛抗旱決策支持系統(tǒng)的一個(gè)模塊，極大地方便了各類應(yīng)用。

(4)平滑處理。為了得到PSD在頻率對(duì)數(shù)坐標(biāo)中呈等間隔采樣，本文采用1/3倍頻程積分作為平滑處理

其中，fl=2-1/6fc為低頻拐角頻率;為fh=21/6fc高頻拐角頻率;n為介于二者之間頻率f的個(gè)數(shù)。由(7)式得到中心頻率fc的PSDa平均值PSDa(fc)，作為fc的加速度功率譜密度的PSD值。中心頻率fc以1/9倍頻程為增加步長(zhǎng)，即下一個(gè)中心頻率fc=21/9fc，重新計(jì)算相應(yīng)的fl和fh，然后將新的fl與fh之間的PSD值平均值作為下一個(gè)中心頻率fc的PSD取值。這樣在fc的取值范圍0.02～40Hz內(nèi)，每個(gè)記錄段的PSD值隨頻率變化情況可由在對(duì)數(shù)坐標(biāo)系呈等間隔采樣的107個(gè)中心頻率的PSD值來表示(圖3)。

2.3 PDF圖上確定臺(tái)站高、低噪聲參照線

每個(gè)中心頻率fc的PSD概率密度函數(shù)為

其中，Nfc為fc頻點(diǎn)的記錄段總數(shù);Npfc為fc頻點(diǎn)的PSD值落在某PSD取值范圍內(nèi)的記錄段個(gè)數(shù)，在本研究中PSD窗長(zhǎng)與步長(zhǎng)都取1dB，變化范圍為-200～-50dB。然后，以頻率為橫坐標(biāo)、以PSD為縱坐標(biāo)、以PPSD(fc)色塊顏色深淺繪制三維平面圖，得到功率譜概率密度函數(shù)(PDF)分布圖(圖1(a))，不同色塊代表某頻點(diǎn)在一定PSD窗內(nèi)功率譜概率數(shù)。

從PDF圖中可以明顯看出臺(tái)站不同噪聲水平所處的區(qū)域。如圖1(a)所示，顏色較深部分(藍(lán)、綠、黃、紅區(qū)域)是臺(tái)站正常噪聲水平落入的區(qū)域。顏色較淺部分(粉紅色區(qū)域)是臺(tái)站異常記錄(地震波、外界噪聲)所處的區(qū)域。為了能描繪出正常噪聲所處的區(qū)域范圍，本文引進(jìn)了臺(tái)站高、低噪聲參照線來表示該區(qū)域的上下限。首先，選取藍(lán)色輪廓作為正常噪聲的界限，由圖1(a)右邊的圖例可以看到藍(lán)色輪廓線所對(duì)應(yīng)的網(wǎng)格概率是1.5%。

針對(duì)圖1(a)中1Hz處的加速度功率譜密度概率分布進(jìn)行分析，以功率譜密度值為橫坐標(biāo)，網(wǎng)格概率為縱坐標(biāo)，畫出對(duì)應(yīng)1Hz頻點(diǎn)的功率譜密度概率圖(圖1(b))。以最大網(wǎng)格概率點(diǎn)為中心，向左、右兩邊分別找首個(gè)網(wǎng)格概率≤1.5%的點(diǎn)，用于確定低噪聲參照點(diǎn)和高噪聲參照點(diǎn)。對(duì)于其它頻點(diǎn)采用同樣方法確定高、低噪聲參照點(diǎn)，最終將這些點(diǎn)連接起來就得到了臺(tái)站網(wǎng)格概率≈1.5%的高、低噪聲參照線(圖1(a))。

圖2 福建省85個(gè)測(cè)震臺(tái)站分布圖

3 噪聲異常的分類與遴選方法

本文對(duì)噪聲異常的遴選思路是利用臺(tái)站過去1年的歷史記錄繪出臺(tái)站PDF模型，由該模型確定出臺(tái)站的高、低噪聲參照線，然后將實(shí)時(shí)計(jì)算的PSD線與高、低噪聲參照線進(jìn)行比較，實(shí)現(xiàn)噪聲異常自動(dòng)化遴選。

根據(jù)PDF圖可以直觀的判斷噪聲是否發(fā)生異常。為了便于研究噪聲異常的自動(dòng)化挑選，本文將噪聲異常大致進(jìn)行了分類，分類規(guī)則如下:PDF圖上顯示的異常分布在低噪聲參照線附近及以下，屬于低噪處異常;PDF圖上顯示的異常分布在高噪聲參照線附近及以上，屬于高噪處異常;PDF圖上顯示介于高、低噪聲參照線之間的異常，則屬于中噪處異常;除此之外，缺數(shù)異常也是常見的一種噪聲異常，一般都沒有數(shù)據(jù)，在PDF圖上無法顯示，所以本文直接在時(shí)域上對(duì)其進(jìn)行挑選。挑選方法:若所有點(diǎn)都小于1.0×10-10m/s(由于正常噪聲速度記錄一般為1.0×10-6m/s左右)，則判定為缺數(shù)異常。由于噪聲異常具有不定性和多樣性，那么只用一種挑選方法是無法完成對(duì)噪聲異常的挑選，通過上述對(duì)噪聲異常的分類，然后研究針對(duì)每個(gè)類型的挑選方法，就可以很好地解決這一難題。

3.1 低噪處異常的挑選

3.2 高噪處異常的挑選

從PDF圖中也可以很容易判斷出高噪處異常，對(duì)造成高噪聲異常的某條波形，并繪出該波形的PSD線(圖3(b)，實(shí)線)，從圖3(b)中可以看出，該類異常的特征為PSD值大于臺(tái)站高噪聲參照線。同樣考慮在確定臺(tái)站高噪聲參照線時(shí)存在誤差影響，則由臺(tái)站高噪聲參照線全體向上移動(dòng)3dB而得到虛線。挑選方法:若PSD線(實(shí)線)超出虛線的點(diǎn)數(shù)占總點(diǎn)數(shù)的百分比大于0.45，則判定為高噪處異常。

圖3 噪聲異常挑選方法示意圖。(a)低噪處異常中5分鐘一小段噪聲記錄計(jì)算出來的PSD線與臺(tái)站低噪聲參照線的比較;(b)高噪處異常中5分鐘一小段噪聲記錄計(jì)算出來的PSD線與臺(tái)站高噪聲參照線的比較;(c)中噪處異常中5分鐘一小段噪聲記錄計(jì)算出來的PSD線的特征

3.3 中噪處異常的挑選

同樣可利用PDF圖判斷出中噪處異常，選取造成中噪聲異常的某條波形，繪出該波形的PSD線(圖3(c)，實(shí)線)，從圖3(c)中可以看出，該類異常的特征是PSD線起伏不大，幾乎成一條直線。由于正常噪聲的PSD線起伏較大，那么它所對(duì)應(yīng)的方差也較大，這樣可以通過比較兩者方差的不同來挑選出該類異常。挑選方法:先用最小二乘法對(duì)PSD線作傾斜校正，再求校正后的PSD線的方差。若方差＜3，則判定為中噪處異常。

3.4 挑選結(jié)果

若要將上述方法運(yùn)用到實(shí)際監(jiān)測(cè)中就需將其按一定的次序整合起來，本文給出了地震噪聲實(shí)時(shí)監(jiān)測(cè)系統(tǒng)流程圖(圖4)，根據(jù)該流程圖用Matlab編寫了對(duì)應(yīng)的挑選程序。選取福建省85個(gè)測(cè)震臺(tái)站2013年7月份的噪聲記錄數(shù)據(jù)進(jìn)行驗(yàn)證，依據(jù)圖4所示的系統(tǒng)流程分別對(duì)85個(gè)臺(tái)站的噪聲數(shù)據(jù)進(jìn)行在線運(yùn)行。運(yùn)行結(jié)果如圖5所示，橫坐標(biāo)表示福建省85個(gè)臺(tái)站，縱坐標(biāo)表示每個(gè)臺(tái)站應(yīng)用地震噪聲實(shí)時(shí)監(jiān)測(cè)系統(tǒng)所挑選出來屬于異常波形的個(gè)數(shù)與所挑選出來的波形總數(shù)的百分比，圖中每個(gè)點(diǎn)表示一個(gè)臺(tái)站所識(shí)別出來異常的正確率。從圖5中可以看出85個(gè)臺(tái)站識(shí)別出來的異常正確率均在90%以上，多數(shù)臺(tái)站達(dá)到100%。這是因?yàn)樵肼曇坏┌l(fā)生異常，那么其對(duì)應(yīng)的功率譜密度PSD值也是異常的，這樣利用臺(tái)站高、低噪聲參照線作為參考標(biāo)準(zhǔn)就可以很容易的將異常挑選出來。

圖4 地震噪聲實(shí)時(shí)監(jiān)測(cè)系統(tǒng)流程圖

圖5 選取2013年7月份噪聲記錄對(duì)福建省85個(gè)臺(tái)站進(jìn)行測(cè)試的結(jié)果

4 結(jié)論與討論

本文利用臺(tái)站1年期的歷史記錄繪出臺(tái)站PDF模型，然后由該模型確定出臺(tái)站的高低噪聲參照線，再根據(jù)這兩條線與實(shí)時(shí)計(jì)算的PSD線進(jìn)行比較，實(shí)現(xiàn)噪聲異常自動(dòng)化挑選。該方法用于噪聲實(shí)時(shí)監(jiān)測(cè)可以取得很好的效果，大大減少了臺(tái)網(wǎng)儀器系統(tǒng)異常檢測(cè)的工作量，且能及時(shí)發(fā)現(xiàn)異常，并提醒臺(tái)網(wǎng)工作者及時(shí)解決故障，保證記錄到的地震數(shù)據(jù)質(zhì)量?jī)?yōu)良。本文研究可以得到如下結(jié)論:

(1)采用網(wǎng)格概率的方法來確定臺(tái)站高低參照線并作為噪聲異常遴選的標(biāo)準(zhǔn)，其中考慮了臺(tái)站臺(tái)基和儀器穩(wěn)定性等復(fù)雜因素的影響，對(duì)研究臺(tái)站本身噪聲異常挑選具有很強(qiáng)的適用性。

(2)地震噪聲實(shí)時(shí)監(jiān)測(cè)系統(tǒng)對(duì)噪聲異常的挑選正確率很高，且普遍適用于福建省85個(gè)臺(tái)站。

(3)通過噪聲異常在PDF圖上所分布的區(qū)域不同，將噪聲異常分成缺數(shù)(斷記)異常、低噪處異常、中噪處異常、高噪處異常等4類，其中囊括了絕大部分異常。高噪處異常較為復(fù)雜，包括了儀器異常、地震波、外界噪聲、系統(tǒng)瞬變等異常，這些異常有些是非儀器異常引起的，需要進(jìn)一步研究對(duì)它們的區(qū)分方法。

(4)通過大量數(shù)據(jù)統(tǒng)計(jì)，給出高、低噪聲異常挑選方法的閥值，這為以后新建臺(tái)站的異常挑選提供了可行的參照標(biāo)準(zhǔn)。通過該參照標(biāo)準(zhǔn)可以很好解決新建臺(tái)站歷史數(shù)據(jù)不足的問題。