基于功率譜概率密度函數(shù)檢測測震臺網(wǎng)JOPENS系統(tǒng)數(shù)據(jù)異常

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

（中國呼和浩特 010010 內(nèi)蒙古自治區(qū)地震局）

0 引言

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

造成地震觀測數(shù)據(jù)異常的因素主要有觀測環(huán)境、觀測設(shè)備運(yùn)行狀態(tài)、JOPENS 系統(tǒng)MySQL 數(shù)據(jù)庫儀器響應(yīng)參數(shù)（靈敏度、零極點(diǎn)、歸一化因子等）準(zhǔn)確性。在國內(nèi)相關(guān)研究中，關(guān)于JOPENS 系統(tǒng)MySQL 數(shù)據(jù)庫儀器響應(yīng)參數(shù)有誤導(dǎo)致觀測數(shù)據(jù)異常的文獻(xiàn)較少，筆者基于功率譜概率密度函數(shù)，來檢測內(nèi)蒙古測震臺網(wǎng)JOPENS 系統(tǒng)的數(shù)據(jù)異常。

1 波形數(shù)據(jù)質(zhì)量分析軟件部署

以內(nèi)蒙古測震臺網(wǎng)JOPENS 系統(tǒng)為基礎(chǔ)，以Continuous Wave Quality Look（CWQL）軟件為支撐，從JOPENS 系統(tǒng)MySQL 數(shù)據(jù)庫實(shí)時(shí)讀取48 個(gè)臺站儀器響應(yīng)參數(shù)，從流服務(wù)獲取48 個(gè)臺站實(shí)時(shí)波形數(shù)據(jù)，準(zhǔn)實(shí)時(shí)自動(dòng)計(jì)算各臺站各通道PSD 值，自動(dòng)繪制PDF 圖，對內(nèi)蒙古測震臺網(wǎng)JOPENS 系統(tǒng)觀測數(shù)據(jù)質(zhì)量進(jìn)行評估。數(shù)據(jù)處理流程如圖1 所示。

圖1 數(shù)據(jù)處理流程Fig.1 The data access flowchart

波形數(shù)據(jù)質(zhì)量分析CWQL 軟件，結(jié)合了JOPENS 系統(tǒng)，利用MATLAB 平臺實(shí)現(xiàn)了基于PDF 方法的背景噪聲自動(dòng)計(jì)算，從而實(shí)現(xiàn)了對數(shù)字地震臺網(wǎng)數(shù)據(jù)記錄質(zhì)量、觀測系統(tǒng)運(yùn)行狀態(tài)以及臺站背景噪聲水平的準(zhǔn)實(shí)時(shí)監(jiān)控（王俊等，2013）。

2 原理

2.1 功率譜密度計(jì)算

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

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

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

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

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

將速度PSD 值轉(zhuǎn)換為加速度PSD，公式如下

需要扣除儀器傳遞函數(shù)影響，以反映真實(shí)地動(dòng)噪聲物理量值，公式如下

式中，Pa(f)為真實(shí)地面運(yùn)動(dòng)加速度功率譜，H(f)為儀器傳遞函數(shù)。

2.2 平滑處理

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

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

2.3 概率密度函數(shù)計(jì)算

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

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

臺站環(huán)境背景噪聲水平的速度均方根值（RMS）計(jì)算公式為

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

3 臺站地動(dòng)速度記錄

地動(dòng)數(shù)據(jù)接入測震臺網(wǎng)JOPENS 系統(tǒng)流程如圖2 所示?？傮w上分為地震計(jì)輸出y(t)和數(shù)據(jù)采集器輸出g(t)，其中地震計(jì)輸出y(t)為數(shù)采輸入信號。

圖2 數(shù)據(jù)采集器接入流程Fig.2 The flowchart of data

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

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

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

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

圖3 JOPENS 系統(tǒng)臺站參數(shù)樹狀結(jié)構(gòu)Fig.3 Parameter tree of a station in JOPENS system

4 數(shù)據(jù)計(jì)算及分析

通過CWQL 軟件，從內(nèi)蒙古測震臺網(wǎng)中心JOPENS 系統(tǒng)讀取48 個(gè)臺站儀器響應(yīng)參數(shù)和原始波形，并進(jìn)行預(yù)處理。對波形數(shù)據(jù)進(jìn)行去均值、去長周期成分處理，將原始波形分段，計(jì)算PSD 值并進(jìn)行平滑處理，計(jì)算PDF 值及RMS 值，畫PDF 圖。把計(jì)算結(jié)果存儲于本地。計(jì)算結(jié)果可按時(shí)間段、臺站進(jìn)行分項(xiàng)調(diào)取、分析。

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

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

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

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

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

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

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

圖6 霍林河臺東西向和呼和浩特臺南北向功率譜概率密度函數(shù)（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）中，下面的曲線是因?yàn)閿?shù)據(jù)采集器未正常連接地震計(jì)，沒有采集到地面運(yùn)動(dòng)而形成的；上面的曲線是因?yàn)楣ぷ魅藛T連接地震計(jì)時(shí)忘記調(diào)地震計(jì)機(jī)械零點(diǎn)，地震計(jì)處于靠擺臨界狀態(tài)，地震計(jì)記錄地面運(yùn)動(dòng)的動(dòng)態(tài)范圍縮小而形成的。圖6（b）中，在1—10 Hz 范圍內(nèi)，PSD 曲線明顯高出是由于呼和浩特臺距市區(qū)較近，受車輛、人為干擾較大。只有臺站遷址才能有效降低該頻段噪聲水平。

5 結(jié)束語

利用觀測數(shù)據(jù)功率譜概率密度函數(shù)，檢測測震臺網(wǎng)JOPENS 系統(tǒng)數(shù)據(jù)異常，可以得出以下結(jié)果：①可以監(jiān)控臺站運(yùn)行質(zhì)量，包括地震計(jì)和數(shù)采運(yùn)行狀態(tài)；②能檢測JOPENS 系統(tǒng)MySQL 數(shù)據(jù)庫儀器響應(yīng)參數(shù)準(zhǔn)確性；③能夠?qū)崟r(shí)評估臺站背景噪聲水平。

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