地震臺站噪聲水平和臺網監測能力自動化程序的實現

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

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

0 引言

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

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

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

1 設計思路與流程

為了更好的達到程序的高可移植性、跨計算機運行等各項功能。本程序能夠實現在Linux操作系統下，通過CLI（命令行界面）或者GUI（圖形用戶界面）運行。

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

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

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

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

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

2 噪聲水平和臺網理論監測能力自動化程序計算與實現方法

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

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

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

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

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

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

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

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

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

2.4 理論監測能力的計算

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

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

3 程序運行

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

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

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

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

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

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

圖5 黑龍江地震臺網理論監測能力Fig.5 Theoretical monitoring capability of Heilongjiang Seismic Network

4 結論與討論

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

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