地震噪聲異常實時監(jiān)測

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

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

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

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

0 引言

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

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

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

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

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

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

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

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

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

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

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

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

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

2.2 地震噪聲功率譜計算過程

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

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

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

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

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

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

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

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

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

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

2.3 PDF圖上確定臺站高、低噪聲參照線

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

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

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

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

圖2 福建省85個測震臺站分布圖

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

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

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

3.1 低噪處異常的挑選

3.2 高噪處異常的挑選

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

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

3.3 中噪處異常的挑選

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

3.4 挑選結(jié)果

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

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

圖5 選取2013年7月份噪聲記錄對福建省85個臺站進行測試的結(jié)果

4 結(jié)論與討論

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

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

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

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

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