胡寶慧，張 浩，教智博，雷凱悅

（黑龍江省地震局，黑龍江 哈爾濱 150090）

0 引言

地震臺站的噪聲研究受到國內(nèi)外專家廣泛關注。Peterson[1]通過分析全球75 個地震臺站近2000 條噪聲記錄的功率譜密度分布，給出全球低噪聲模型（NLNM）與高噪聲模型（NHNM），目前廣泛應用于臺址噪聲水平評價與不同背景噪聲水平下地震計響應的預測。自2004 年以來，McNamara[2]認為，噪聲功率譜密度（PSD）估計能較好確定臺站噪聲水平，起到實時監(jiān)控和及時發(fā)現(xiàn)儀器故障的作用，國內(nèi)許多專家[3-4]對以上方法進行了應用，取得了較好的效果。但現(xiàn)有噪聲水平計算程序存在安裝設置困難、數(shù)據(jù)格式受限、參數(shù)配置復雜等諸多問題，長期缺乏一套具備自動數(shù)據(jù)接入、處理運算和繪圖功能的數(shù)據(jù)處理軟件。大部分軟件程序還需要手動收集處理觀測數(shù)據(jù)的情況，存在速度慢、效率低、開發(fā)后難維護等問題，遠遠滿足不了日常業(yè)務工作的需要。

地震臺網(wǎng)監(jiān)測能力，是地震臺站布局及單臺監(jiān)測能力的綜合體現(xiàn)。地震臺網(wǎng)監(jiān)測能力除了與臺站的布局情況有關之外，還與臺站的背景噪聲水平、觀測系統(tǒng)響應靈敏度、儀器動態(tài)范圍有關。本程序基于地震噪聲水平和震級衰減關系的理論監(jiān)測能力評估方法[5]，以臺網(wǎng)內(nèi)所有臺站噪聲水平計算結果作為輸入，依據(jù)近震震級計算公式，建立震級大小和震中距的對應關系，按照0.1°×0.1°劃分網(wǎng)格得到整個臺網(wǎng)的理論監(jiān)測能力。其結果可用于分析區(qū)域內(nèi)測震臺網(wǎng)地震監(jiān)測能力變化。

本次開發(fā)全部程序采用解釋性腳本語言，極大簡化了“開發(fā)、部署、測試和調(diào)試”的周期，為自動化單臺噪聲水平計算和臺網(wǎng)的理論監(jiān)測能力的快速估算提供了可能。

1 設計思路與流程

為了更好的達到程序的高可移植性、跨計算機運行等各項功能。本程序能夠?qū)崿F(xiàn)在Linux操作系統(tǒng)下，通過CLI（命令行界面）或者GUI（圖形用戶界面）運行。

程序主要包含文件處理、數(shù)據(jù)處理和繪圖三個模塊。其中文件處理模塊完成數(shù)據(jù)的接入、小時波形數(shù)據(jù)格式的轉換、文件合并、文件重命名等工作；數(shù)據(jù)處理模塊完成儀器信息文件的解析、噪聲水平的計算和理論監(jiān)測能力的計算；繪圖模塊完成各類圖件的快速繪制。模塊間采用函數(shù)調(diào)用方式實現(xiàn)。

圖1 程序設計思路Fig.1 Program design ideas

圖2 程序計算流程Fig.2 Program calculation flow

程序編制過程中，根據(jù)噪聲水平和臺網(wǎng)監(jiān)測能力的設計思路和計算流程，通過JOPENS系統(tǒng)中AWS 服務獲取所需臺站該時間段內(nèi)小時波形數(shù)據(jù)服務，然后采用Perl 腳本程序調(diào)用SAC 對波形數(shù)據(jù)進行自動批量處理。首先將SEED 格式的數(shù)據(jù)轉換成SAC 格式，然后將可能因多種因素出現(xiàn)間斷的同一通道的SAC 數(shù)據(jù)合并起來完成小時波形的檢查和數(shù)據(jù)合并，并對文件進行統(tǒng)一重命名。該流程將得到各臺站SAC 格式的小時波形文件和對應的文本格式的儀器信息文件，其中通過解析儀器信息文件獲得的數(shù)據(jù)采集器轉換因子、地震計靈敏度、零級點參數(shù)等信息，可用于建立傳遞函數(shù)模型。

SAC 提供的命令可以實現(xiàn)地震數(shù)據(jù)的預處理，但無法實現(xiàn)所有的數(shù)據(jù)分析功能。因此在數(shù)據(jù)處理階段，本程序使用SAC 中自帶的I/O接口程序ReadSac.m 讀寫SAC 文件的MATLAB腳本。通過MATLAB 腳本程序?qū)崿F(xiàn)SAC 格式文件固定長度的頭段區(qū)和非固定長度的數(shù)據(jù)區(qū)的讀取，進而開展數(shù)據(jù)的計算[6]。

2 噪聲水平和臺網(wǎng)理論監(jiān)測能力自動化程序計算與實現(xiàn)方法

2.1 噪聲功率譜（PSD）的計算

在噪聲功率譜（PSD）估計方面，本程序應用基于FFT 的非參數(shù)功率譜估計方法Welch平均周期法通過分段選取數(shù)據(jù)進行加窗求功率，再進行平均計算。

程序編制過程中，選用Hamming 窗函數(shù)，數(shù)據(jù)置信水平為95%。首先根據(jù)地震計傳遞函數(shù)，去除儀器響應，得到地面運動的速度記錄；然后將1 小時的數(shù)據(jù)段分成14 小段，每小段數(shù)據(jù)的長度為1000 S，數(shù)據(jù)段重合率為80%，對于每一小段數(shù)據(jù)去均值與線性趨勢；計算所有14 小段的速度功率譜密度，取平均值，得到該小時數(shù)據(jù)的速度功率譜密度隨頻率的分布；再次對得到的速度功率譜密度進行1/8 倍頻程濾波，得到平滑且在對數(shù)坐標上均勻分布的速度功率譜密度；最后將速度功率譜密度轉化為加速度功率譜密度，并與NLNM、NHNM 進行對比分析，分不同的頻段研究加速度功率譜密度的特征。離散噪聲信號的功率譜密度可以表示為[7-8]：

同時可將噪聲的單位用dB 表示為：

2.2 功率譜概率密度（PDF）的計算

因臺站的噪聲水平隨著外界環(huán)境變化，功率譜密度也隨時間不斷變化，為了評估地震臺站的環(huán)境變化特征，在PSD 值計算結束后，采用McNamara 的方法對1/8 倍頻程濾波后的功率譜密度進行統(tǒng)計。以1 dB 為間隔，因為大多數(shù)功率譜密度都在-200～80 dB 范圍內(nèi)，因此本文以-200～80 dB 為統(tǒng)計范圍，然后計算中心頻率f 范圍內(nèi)的所有功率譜密度的數(shù)量Nf，Ndf為f處功率譜密度在（d～d+1）dB 范圍內(nèi)的個數(shù)，則（f，d）處的概率密度為：P（f，d）=Ndf/Nf。然后計算值落在中心點頻率記錄段的數(shù)目，統(tǒng)計各周期PSD 在不同時間內(nèi)取某一數(shù)值的概率，得到功率譜密度在-200～80 dB 上的統(tǒng)計分布，即為該時段對應的PDF。

2.3 地脈動噪聲均方根值（RMS）的計算

地脈動噪聲均方根值RMS 可以衡量地震臺臺基背景噪聲水平[6]，是求解臺網(wǎng)監(jiān)測能力的基礎。在計算RMS 時，對中心頻點在1～20 Hz頻帶范圍內(nèi)所有PSD 值的振幅值取平方，再求其在該時段內(nèi)的平均，然后求其平方根。對于離散信號，均方根計算公式為：

其中，n 是樣本數(shù)，xi是第i 個樣本的幅值。

2.4 理論監(jiān)測能力的計算

理論監(jiān)測能力方面，以臺站為中心、最大記錄距離為半徑畫圓，求3 個以上臺站的交集，該交集即為某一震級的監(jiān)測區(qū)。在程序編制過程中，把臺網(wǎng)監(jiān)測區(qū)劃分為N×M 的網(wǎng)格，各網(wǎng)格點所在坐標為（Xij，Yij），（i=1，…，N；j=1，…，M）。假設網(wǎng)格點為震中，把臺網(wǎng)的全部臺站按下式計算出震級ML，然后將震級ML由小到大排序，當網(wǎng)格點處發(fā)生第四號震級的地震時，則至少有3 個臺站可以記錄到，定義第四號震級為該網(wǎng)格點的理論最小監(jiān)測震級。在2.3 節(jié)計算獲得頻帶1～20Hz 均方根值RMS的基礎上，根據(jù)目前地震臺網(wǎng)用速度型記錄測定近震體波的公式計算得到各臺可監(jiān)測到的最小震級：

式中，Vs為臺站記錄的體波最大速度，f對應最大速度頻率值。R（Δ）為量規(guī)函數(shù)，即震級起算函數(shù)[9]。

3 程序運行

選取黑龍江省及周邊區(qū)域74 個臺站222 個通道的連續(xù)波形數(shù)據(jù)，采用上述噪聲計算方法，獲得各臺站的小時噪聲計算結果。在程序?qū)嶋H運行中，首先對單臺48 小時連續(xù)數(shù)據(jù)進行功率譜概率密度函數(shù)（PDF）的計算，以低噪聲模型（NLNM）與高噪聲模型（NHNM）參照線作為衡量該臺站噪聲水平的標準。然后通過任意小時波形數(shù)據(jù)與48 小時參照數(shù)據(jù)的比對判定該小時數(shù)據(jù)偏差，以衡量該小時波形的質(zhì)量。

選取BWS 和HEH 兩個臺站（圖3）的數(shù)據(jù)計算結果進行分析，其中BWS 臺為深井臺站，儀器型號為BBVS-60，HEH 臺為地面臺站，儀器型號為JCZ-1。兩個臺站的NS 分量加速度功率譜概率密度函數(shù)分布如圖4 所示。

圖3 小時波形噪聲功率譜（PSD）檢驗結果Fig.3 Hourly waveform noise power spectrum inspection results

圖4 連續(xù)48 小時概率密度函數(shù)PDFFig.4 PDF of seismic station for continuous 48 hours

對于固定臺站來說，以上兩個臺站儀器參數(shù)正常，臺基環(huán)境噪聲水平在頻域的分布特征是相對穩(wěn)定的。BWS 臺作為深井臺本應觀測環(huán)境優(yōu)良，但圖3（a）顯示長周期噪聲偏差明顯，且48 小時均超出高噪聲模型（NHNM），需要查明該臺站是何種原因?qū)е碌拈L周期頻段異常。HEH 臺雖出現(xiàn)該小時微震頻段噪聲水平異常，但明顯從圖3（b）可以看出是因為發(fā)生地震所致，且48 小時PDF 顯示該臺站儀器整體穩(wěn)定，噪聲水平在合理區(qū)間內(nèi)。

從理論監(jiān)測能力運行結果看，除黑龍江西北部局部地區(qū)外，全省測震臺網(wǎng)基本已達ML2.5 以上地震監(jiān)測能力，且中東部地區(qū)總體監(jiān)測能力較強，這與監(jiān)控黑龍江內(nèi)最大斷裂—依舒斷裂的目標較為一致。從圖5 看出，監(jiān)測能力不僅與臺站密度增加有關，更與該區(qū)域臺站記錄質(zhì)量好壞密不可分，HEH 區(qū)域因臺站稀疏導致監(jiān)測能力較弱，而BWS 區(qū)域臺站稀疏與記錄質(zhì)量不佳更導致了監(jiān)測能力的下降，其他如大慶及哈爾濱區(qū)域雖然臺站密度較高，但數(shù)據(jù)質(zhì)量不高。因此，地震監(jiān)測能力計算對優(yōu)化測震臺網(wǎng)布局，合理布置臺站，避免盲目增加地震臺站數(shù)量起到重要作用[10]。

圖5 黑龍江地震臺網(wǎng)理論監(jiān)測能力Fig.5 Theoretical monitoring capability of Heilongjiang Seismic Network

4 結論與討論

地震臺站噪聲水平和臺網(wǎng)監(jiān)測能力自動化程序完成自動分析地震臺站波形數(shù)據(jù)質(zhì)量的工作，產(chǎn)出了相應結果的圖件。實現(xiàn)了區(qū)域內(nèi)臺站波形的實時接收和處理、自動監(jiān)控波形質(zhì)量的目的，為評估臺站的地震噪聲水平和臺網(wǎng)監(jiān)測能力提供了充足且可靠的依據(jù)。經(jīng)過長時間試運行表明，該套程序能夠滿足目前我省臺站數(shù)據(jù)質(zhì)量監(jiān)控的絕大部分需求，自動產(chǎn)出結果能為后期的臺網(wǎng)規(guī)劃和臺站效能評估服務。

對于地震臺網(wǎng)的地震監(jiān)測能力評估，本程序采用的震級估算法雖然能夠較為客觀地確定區(qū)域內(nèi)最小控震能力，但仍受噪聲質(zhì)量的影響。今后將通過程序改進，增加依靠臺網(wǎng)實際產(chǎn)出為依據(jù)的諸如PMC 等計算監(jiān)測能力的方法，實現(xiàn)地震臺網(wǎng)監(jiān)測能力更科學的評估。而在噪聲水平分析方面，如何實現(xiàn)與地震目錄、儀器實時參數(shù)配合分析噪聲功率譜，更精確表述儀器的工作狀態(tài)也是該程序需要逐步改進的部分。