利用背景噪聲的分形方法監控臺站觀測系統運行狀態

丁莉莎 馬潔美 齊軍偉 謝劍波 廖一帆 盧子晉 葉世山 勞 謙 呂仲杭

1 廣東省地震局，廣州市先烈中路81號，510070 2 中國地震局地球物理研究所，北京市民族大學南路5號，100081 3 湖北省地震局，武漢市洪山側路40號，430071

背景噪聲是地震計在未發生地震時拾取的在地球內部傳播的彈性波，包含地球內部的全部信息，可用于地殼結構成像、火山爆發監測等。關于地球背景噪聲的研究主要集中在2個方面：1)地球噪聲模型、噪聲產生的機制和源。如Peterson[1]提出地球新噪聲模型；Stehly等[2]通過求取噪聲相關函數的包絡，獲得歸一化的噪聲背景能量流(normalized background energy flux, NBEF)；Stutzmann等[3]分析全球氣候對背景噪聲的時間演化關系。2)地球背景噪聲的利用。根據背景噪聲的自互相提取格林函數進行地震成像研究[4]以及地震傳感器其他性能研究，包括估算地震計自噪聲[5]、地震計方位角偏差[6-7]以及數據質量評估[8]，其中針對記錄數據質量評估的研究主要集中在能量-頻率分布特征方面。本文將采用分形新方法對地震傳感器記錄的數據進行質量評價。

1 地震背景噪聲的分形性

1.1 時間序列分形屬性

如果時間序列x(t)具有自相似或自仿射等分形特征，則對于任何τ≥1,t1,t2,…tτ∈T和a>0，x(t)應滿足標度律：

[x(at1),x(at2),…,x(atτ)]=

[anx(t1),anx(t2)…anx(tτ)],n>0

(1)

式中，a為標度指數。由x(t)滿足標度律可知，其功率譜密度S(f)滿足負冪律：

S(f)∝f-β

(2)

式中，β為功率譜密度指數。S(f)滿足負冪律，且1<β≤3為該時間序列分形的必要條件，即PSD法。

在實際數據處理中，實驗對象不可能嚴格服從該規律，但在相當寬度的頻帶范圍內滿足標度律時可解釋為分形[9]。

1.2 地球背景噪聲分形性質

通過對比Peterson[1]基于全球75個臺站記錄得到的新噪聲模型(N噪聲模型)、Khairy等[10]基于埃及國家地震臺網7臺寬頻帶地震計記錄得到的E噪聲模型以及El Fellah等[11]基于摩洛哥國家地震臺網23臺寬頻帶地震計記錄得到的M噪聲模型(圖1)，得出第二類地脈動中0.2～1 Hz頻段內功率譜密度滿足負冪律，因此推斷該頻段背景噪聲具有分形性。基于該性質，下文將重點介紹該頻帶范圍內背景噪聲的分形分析與應用。

圖1 不同地球背景噪聲模型對比Fig.1 Comparison of different Earth backgroundnoise models

2 背景噪聲分形盒維數模型

2.1 時間序列計盒維數法

分形維數(fractal dimension)是分形理論中最核心的概念與內容[12]。其中，分形盒維數表明波形(圖像)在一定標度下的不規則性或復雜性與波形頻譜有關，也表示復雜性隨標度減小而加劇的速率[13]。

對分形曲線y=f(x)，用尺度為δ的網格對其進行覆蓋，所需網格數為Nδ，則Nδ隨δ的減小而增大，且滿足：

Nδ=kδ-D

(3)

式中，D為曲線的分維數，也稱為計盒維數，k為整數。

一維時間序列分形盒維數法采用已知寬度的柵格覆蓋分形曲線(圖2)，柵格尺度(寬度)為δ。如第i個柵格內曲線段的計盒維數為：

(4)

通過構造已知Hurst值的一維分形波形并進行多次比較驗證，結果表明，上述算法能很好地反映理論盒維數的變化情況。

圖2 一維時間序列分形盒維數原理Fig.2 The diagram of fractal box dimension of one-dimensional time sequence

2.2 數據背景

本文選取國家臺網中心及廣東省地震臺網提供的23個在運行臺站2015年的數據及中國地震局地球物理研究所提供的17個實驗設備的記錄數據進行分析，實驗數據有效涵蓋了我國臺站運行環境。所選臺站中，緯度最高的約53°N為漠河臺(MOH)，最低的約16°N為西沙臺(XSA)，實驗范圍跨越中國各個緯度帶；在地域高程方面，都蘭臺(DUL)處于3 500 m以上的高寒地區。

表1為本文使用的臺站和地震計類型，包含我國在運行的4 種主流地震計及17種國內外在運行或即將投入使用的地震計，可使實驗結果具有一般性和普適性。

3 實驗結果與分析

小波的自適應分解能力能有效地對背景噪聲進行尺度(頻段)劃分，滿足在一定尺度內對背景噪聲進行自相似性分析的要求。根據本文所采用的地震儀的采樣率(100 Hz)，選擇哈爾小波為小波基進行8階小波分解，提取子帶6至子帶8的數據(頻帶范圍為0.2～1 Hz)進行盒維數估算。

圖3為不同尺度下盒維數分布圖，子帶6的盒維數集中在1.57～1.59之間，子帶7的盒維數集中在1.49～1.50之間，子帶8的盒維數集中在1.42～1.44之間，表明正常運行的寬頻帶地震計背景噪聲在特定尺度下的分形盒維數穩定收斂。上述結果也表明，由于0.2～1 Hz為21種寬頻帶地震計(包括甚寬頻帶及超寬頻段)的通頻帶，該尺度下分形盒維數與地震計種類及系統參數無關；同時，臺站的地理位置和觀測條件對實驗結果影響較小。

為進一步檢驗分形對儀器運行狀態的監控能力，對部分已發生基本故障的地震計數據及部分地震事件進行分析。圖4為地震事件對實驗結果的影響，圖4(a)為正常地脈動數據疊加近場干擾及地震事件的盒維數結果，其中紫色圈盒維數變化對應圖4(b)。從圖4(a)可以看出，近場地震事件的盒維數結果類似于近場干擾，盒維數變化與能量成反比，且與子帶的頻帶范圍有關，這種異常形態在背景噪聲進行分形處理后較易識別。

表1 地震儀器分類及參數

圖3 不同子帶盒維數頻度圖Fig.3 Frequentness of box dimensionof different subbands

圖5中地震儀處于不穩定工作狀態，儀器異常時盒維數會出現明顯突變。圖5(a)為正常地脈動數據疊加近場干擾及故障的盒維數結果，其中方框盒維數變化對應圖5(b)。圖5(b)為地震計零點漂移直至逐漸靠擺的過程，從圖中可以看出，故障時段的盒維數與近場干擾存在差異，盒維數變化與子帶的頻帶范圍無關。近場干擾表現為盒維數變小的擾動，類似于圖5中近震事件；零漂至靠擺的盒維數表現為整體平行跳變。

圖4 地震事件影響Fig.4 Influence of earthquake event

圖5 故障事件數據分析Fig.5 Data analysis of fault event

結合盒維數結果與頻域分析結果，根據采樣定律可判斷采集的數據頭文件中采樣率是否真實有效，當實際數據的采樣率與本文默認采樣率100 Hz不符合時，所計算的子帶頻域范圍會因不滿足標度律而不具有分形性。

4 監控觀測系統運行狀態的可行性探討

利用噪聲概率密度方法(圖6)和分形盒維數方法(圖7)對NY臺3個月的數據進行處理。由噪聲概率密度方法可知，該時間段內臺基噪聲頻譜相對穩定，存在2個1.2～1.5 Hz頻帶的未知干擾，同時也存在少量異常記錄。由圖7可知，臺基噪聲(0.2～1 Hz)波形在602 h(2015-10-26 07:00)零點漂移逐漸至靠擺；620 h(2015-10-27 01:00)地震觀測系統無輸出；745～749 h(2015-11-01 06:00～10:00)更換設備并進行電標定測試等，異常波形如圖8所示。

圖7 NY臺噪聲分形盒維數Fig.7 Fractal box dimensions of NY station

與功率譜方法相比，分形盒維數方法可定量分析0.2～1 Hz頻段范圍內的異常波形，能準確定位到異常波形的時間。同時，該方法及結果與儀器參數無關，但該方法缺少全頻帶信息。可以看出，噪聲功率譜方法及分形盒維數方法各有優缺點，通過噪聲分析評價波形質量只能用來初步估計儀器的運行狀態。

在測震臺網產出維護中，現有的監測地震儀運行狀態的方法主要是定期使用電流標定線圈激勵法，根據標定數據計算結果來檢查和驗證地震計性能，包括階躍標定、正弦標定以及偽隨機碼標定。與將噪聲功率譜方法作為輔助分析類似，分形盒維數方法可通過對地震計背景噪聲數據進行分形計算，獲取盒維數并定量評價波形質量，進一步反映儀器工作狀態，從而減少信息中斷，遠程監測地震計的運行狀態。

圖8 NY臺異常波形Fig.8 Abnormal waveform of NY station

5 結 語

利用地震計的日常記錄評估地震計運行狀態及數據質量，對提升臺站觀測系統運維能力具有實用意義。功率譜概率密度方法可分析頻域范圍內能量分布的相對變化，但無法在時間域內對地震儀器狀態進行定量檢測；而分形方法提出的時域波形復雜度研究可定量評估特定頻段內波形狀態。因此，相應標度下背景噪聲的自相似性可與地震監測相結合，利用背景噪聲在一定尺度下(頻帶范圍內)的分形性質對地震觀測儀器進行定量分析。在監控地震觀測系統運行狀態方面，該方法具有定量化程度高、抗干擾能力強、不受觀測系統約束、可實現準實時化等優點。本文方法也可用于數據質量評價，同時實現非人為化操作并提高效率，進一步完善我國地震監測和預警系統的運維。