基于功率譜概率密度函數檢測測震臺網JOPENS系統數據異常

安 全 張劍影 翟 浩 蘇日亞 趙鐵鎖

（中國呼和浩特 010010 內蒙古自治區地震局）

0 引言

地震觀測數據是國家防震減災工作的基礎性戰略資源。數據質量直接影響防震減災工作效果。高質量數據是決定地震臺站與地震臺網成為國際一流的先決條件。因此，保障地震觀測數據質量是地震監測工作的核心任務（鄭秀芬等，2017）。Peterson 等（1993）通過對全球正常背景噪聲功率譜（PSD）的研究，給出了全球低噪聲新模型（NLNM）和高噪聲新模型（NHNM）。目前，該模型被廣泛應用于地震臺站噪聲水平評價。McNamara 等（2005）提出了地震噪聲功率譜概率密度函數（Probability Density Function，PDF）方法。該方法可用于臺站噪聲水平和波形質量的測定。PDF 方法直接使用連續波形數據，不對地震數據進行篩選，因此，可同時得到地震事件波形、觀測系統儀器響應參數變化及儀器故障等概率分布。目前PDF 方法已被用于意大利（Marzorati et al，2006）、新西蘭（Rastin et al，2012）、倫敦中部地區（Green et al，2017）的地震環境噪聲特征分析。在國內，應用PSD 方法，孫晴等（2007）分析了河北省地震局數字地震臺網子臺地動噪聲水平，岳秀俠等（2009）分析了天津測震臺網子臺地脈動；應用PDF 方法，廖詩榮等（2008）實現自動處理地震臺站勘選數據，徐嘉雋等（2010）檢測了地震儀器系統工作狀態，林彬華（2015）實現實時監測地噪聲，李小軍等（2019）分析地震計自噪聲水平。王俊等（2013）結合JOPENS 系統，利用MATLAB 平臺研發了運用PDF 方法自動計算臺站背景噪聲的程序，從而實現了對數字地震臺網觀測系統運行狀態以及臺站背景噪聲水平的準實時監控。

造成地震觀測數據異常的因素主要有觀測環境、觀測設備運行狀態、JOPENS 系統MySQL 數據庫儀器響應參數（靈敏度、零極點、歸一化因子等）準確性。在國內相關研究中，關于JOPENS 系統MySQL 數據庫儀器響應參數有誤導致觀測數據異常的文獻較少，筆者基于功率譜概率密度函數，來檢測內蒙古測震臺網JOPENS 系統的數據異常。

1 波形數據質量分析軟件部署

以內蒙古測震臺網JOPENS 系統為基礎，以Continuous Wave Quality Look（CWQL）軟件為支撐，從JOPENS 系統MySQL 數據庫實時讀取48 個臺站儀器響應參數，從流服務獲取48 個臺站實時波形數據，準實時自動計算各臺站各通道PSD 值，自動繪制PDF 圖，對內蒙古測震臺網JOPENS 系統觀測數據質量進行評估。數據處理流程如圖1 所示。

圖1 數據處理流程Fig.1 The data access flowchart

波形數據質量分析CWQL 軟件，結合了JOPENS 系統，利用MATLAB 平臺實現了基于PDF 方法的背景噪聲自動計算，從而實現了對數字地震臺網數據記錄質量、觀測系統運行狀態以及臺站背景噪聲水平的準實時監控（王俊等，2013）。

2 原理

2.1 功率譜密度計算

周期時間序列y(t)的有限范圍傅里葉變換可表示為

式中，Tr為時間序列段長度，f為頻率。

對離散頻率值fk，傅里葉變換定義為

式中，fk=k/(NΔt),其中k=1,2,3，…，N-1，Δt為采樣間隔（0.01 s），N=Tr/Δt為截取時間段的采樣點數。

根據維納—辛欽定理，功率譜密度（PSD）定義為

將速度PSD 值轉換為加速度PSD，公式如下

需要扣除儀器傳遞函數影響，以反映真實地動噪聲物理量值，公式如下

式中，Pa(f)為真實地面運動加速度功率譜，H(f)為儀器傳遞函數。

2.2 平滑處理

為了使PSD 在對數坐標中等間隔分布，采用1/3 頻率積分的平滑處理方法

其中，fl=2-1/6fc，為低頻拐角頻率；fh=21/6fc，為高頻拐角頻率；n為介于二者之間頻率f的個數。由式(6)得到中心頻率fc的Pa(f)平均值Pa(fc)。以Pa(fc)作為fc的加速度功率譜密度，中心頻率fc以1/9 頻程為增加步長，即下一個中心頻率為原中心頻率的21/9倍，重新計算相應的fl和fh，并將新的fl和fh之間的PSD 平均值作為下一個中心頻率fc的PSD 值。這樣，在fc的取值范圍0.02—40 Hz 內，每個記錄段的PSD 值可由在對數坐標系中間隔分布的中心頻率的PSD 值表示。

2.3 概率密度函數計算

每個中心頻fc率的PSD 概率密度函數為

其中，Nfc為fc頻點的記錄段總數，為fc頻點的PSD 值落在某PSD 取值范圍內的記錄段個數。在本研究中，PSD 窗長與步長均取1 dB，變化范圍為-200— -50 dB。以頻率為橫坐標，以PSD 為縱坐標，以顏色深淺表示PPSD（fc），繪制三維平面圖，即為功率譜概率密度函數(PDF)分布圖。不同色塊代表某頻點在一定PSD 窗內的功率譜概率。

臺站環境背景噪聲水平的速度均方根值（RMS）計算公式為

式中，RBW=(fh-fl)/fc為相對寬度。

3 臺站地動速度記錄

地動數據接入測震臺網JOPENS 系統流程如圖2 所示。總體上分為地震計輸出y(t)和數據采集器輸出g(t)，其中地震計輸出y(t)為數采輸入信號。

圖2 數據采集器接入流程Fig.2 The flowchart of data

圖2 中，Y(f)為y(t)的傅里葉變換，X(f)為x(t)的傅里葉變換，H(f) 為地震計傳遞函數，H(f)可表示為

其中，Sd為地震計靈敏度，A0為歸一化因子，(z1，z2，…，zm)為m個零點，(p1，p2，…，pn)為n個極點。

測震臺網JOPENS 系統存儲的數據g(t)單位為counts，將數據還原為速度的頻域計算公式為

其中G(f)為g(t)傅里葉變換，K為數據采集器轉換因子倒數（靈敏度）。JOPENS 系統接入的波形數據為無量綱數值，并非真正的地動速度或地動位移，必須根據公式（10）扣除觀測數據儀器響應（傳遞函數）才能得到地動速度，積分后得到地動位移。而K和H(f)以靈敏度（Sd）、零極點（zm，pn）、歸一化因子（A0）等的形式存儲于JOPENS 系統MySQL 數據庫數中，如圖3 所示。上述參數輸入有誤或者參數的變化均會影響觀測數據質量。上述錯誤或變化以波形瀏覽等方式無法分辨，而應用PDF 方法生成的任意時間段PSD 曲線、PDF 圖則可以直觀反映觀測數據“健康”狀態。

圖3 JOPENS 系統臺站參數樹狀結構Fig.3 Parameter tree of a station in JOPENS system

4 數據計算及分析

通過CWQL 軟件，從內蒙古測震臺網中心JOPENS 系統讀取48 個臺站儀器響應參數和原始波形，并進行預處理。對波形數據進行去均值、去長周期成分處理，將原始波形分段，計算PSD 值并進行平滑處理，計算PDF 值及RMS 值，畫PDF 圖。把計算結果存儲于本地。計算結果可按時間段、臺站進行分項調取、分析。

選取內蒙古測震臺網48 個臺站7 月20 日14 時至2019 年7 月27 日14 時JOPENS 系統內的PDF 值評估樣本（168 h，共計24 192 條）和1—20 Hz 頻段垂直向RMS 值數據，從樣本中挑選典型例子進行分析。

參照《地震臺站觀測環境技術要求第1 部分：測震》（GB/T 19531.1—2004），以48個測震臺站垂直向RMS 值繪制臺站噪聲分級空間分布圖（圖4）。由圖4 可知，內蒙古測震臺網48 個臺站背景噪聲水平屬于Ⅰ類的有45 個，屬于Ⅱ類噪聲水平的臺站有3 個，整體噪聲水平較低。Ⅱ類噪聲水平臺站均位于市區，噪聲源主要為人為活動、交通等。

圖4 48 個臺站背景噪聲分級空間分布Fig.4 Background noise classification interspace distribution of 48 stations

圖5 為以烏加河臺NS 和EW 向連續記錄為樣本的計算結果，可見兩分向PDF 呈明顯差異。在圖5（a）中，NS 向PSD 曲線在整個頻段呈上下2 條曲線。在圖5（b）中，低于0.06 Hz 頻段PSD 曲線呈2 條曲線。

圖5（a）出現上下2 條曲線，是由于JOPENS 系統MySQL 數據庫中的地震計傳遞函數H(f)的歸一化因子A0比儀器實際歸一化因子小1 個數量級。源數據扣除儀器響應時，降低了H(f)的值，從而高估了PSD 值。下面曲線為A0更正后正常PSD 曲線。在圖5（b）中，低于0.06 Hz 頻段PSD 曲線呈2 條曲線是由于MySQL 數據庫中H(f)零極點（zm，pn）與實際值不一致。扣除儀器傳遞函數影響時，計算地震計頻帶寬度有誤。參數更正后，低于0.06 Hz 頻段PSD 曲線（下）值明顯降低，回歸正常水平。

圖5 烏加河臺水平向功率譜概率密度函數（PDF）分布Fig.5 Probability density function（PDF）distribution of horizontal components of the Wujiahe station

圖6 為以霍林河臺EW 向和呼和浩特臺NS 向連續記錄為樣本的計算結果。圖6（a）中PSD 曲線呈上下2 條曲線，且不在全球低噪聲新模型（NLNM）和全球高噪聲新模型（NHNM）范圍內，曲線形狀也與NLNM 和NHNM 不一致。圖6（b）中呼和浩特臺NS 向PSD 曲線基本在NLNM 和NHNM 范圍內，形狀基本一致，但在1—10 Hz 范圍內PSD曲線明顯凸出。

圖6 霍林河臺東西向和呼和浩特臺南北向功率譜概率密度函數（PDF）分布Fig.6 Probability density function（PDF）distribution of the EW component of the Huolinhe station and the NS component of the Wujiahe station

圖6（a）中，下面的曲線是因為數據采集器未正常連接地震計，沒有采集到地面運動而形成的；上面的曲線是因為工作人員連接地震計時忘記調地震計機械零點，地震計處于靠擺臨界狀態，地震計記錄地面運動的動態范圍縮小而形成的。圖6（b）中，在1—10 Hz 范圍內，PSD 曲線明顯高出是由于呼和浩特臺距市區較近，受車輛、人為干擾較大。只有臺站遷址才能有效降低該頻段噪聲水平。

5 結束語

利用觀測數據功率譜概率密度函數，檢測測震臺網JOPENS 系統數據異常，可以得出以下結果：①可以監控臺站運行質量，包括地震計和數采運行狀態；②能檢測JOPENS 系統MySQL 數據庫儀器響應參數準確性；③能夠實時評估臺站背景噪聲水平。

經過地震烈度速報與預警工程項目的建設，內蒙古地區臺站數量不斷增加，如何高效維護臺站、實時檢測、及時發現儀器設備故障是目前面臨的難題。使用功率譜概率密度函數方法評估臺站背景噪聲水平，可以幫助運維人員快速發現觀測系統中存在的異常信息，從而縮短運維周期，在一定程度上解決目前臺站數量增多、運維能力不足的問題。