基于MATLAB的測震數據分析處理方法的設計與應用

徐文海，丁 勇

（四川省地震局，四川 成都 611730）

在數字化觀測系統中，系統的瞬態變化、儀器毛刺、數據記錄的階躍和尖峰以及由系統故障引起的信號失真和重大環境干擾源的出現等，都是影響數據記錄質量的重要因素。除了真實的地震信號外，還包含了由于系統故障引起的各種失真信號以及臺站周圍的環境噪聲干擾等信息[1]。特別是地震烈度速報與預警工程的大背景下，如何應對地震臺站已存在及未來可能存在的各種類型的干擾，是我們需要長期面對的問題。這也就要求地震臺站或中心站的監測人員需要依托自身熟悉自己臺站或轄區內臺站觀測環境的優勢，及時對可疑的干擾進行分析研究。本文認為在臺站干擾分析排查這項工作中，臺站監測人員應該盡可能地將工作做在前面，才能有效避免潛在的負面影響。為了區分以及消除干擾，首先就需要對地震波進行分析。為此，本文采用了傅里葉變換、希爾伯特-黃變換進行頻譜分析，再結合功率譜對地震波干擾信息進行區分，根據干擾的頻率區間設計相應的濾波器來消除干擾。

1 成都地震監測中心站職責

成都地震監測中心站的主責主業分為3部分：①圍繞地震危險區地震活動態勢及未來地震形勢的嚴峻性，不斷提升區域地震監測預報綜合能力。將日常監測預報工作與重點監視防御區的強化短臨跟蹤工作相結合，開展地震震情監視跟蹤和年度地震趨勢會商工作。②制定針對性監測業務拓展方案。以專業臺與市州局共促區域性地震監測事業發展為基點，與市縣防震減災部門協同合作，共建共享，利用成都地震監測中心站技術力量優勢和地方良好的監測基礎資源，優化區域監測臺網布局，提升區域地震監測能力。③推動地震預警項目開展。組織人員開展預警項目建設實施指導服務和驗收工作，目前已完成轄區60余個站點的驗收工作。待地震速報與預警項目建設全面完成后，成都站將承擔所轄片區內284個預警臺共計305臺預警儀器的運行維護工作。

2 測震觀測系統

成都站是龍門山地震帶和四川中北部地區唯一的多學科綜合型觀測臺站。近幾年，隨著“中國數字地震觀測網絡建設項目”“東半球空間環境地基綜合監測子午鏈項目”“中國地震局背景場項目”“汶川地震災后恢復重建項目”“蘆山地震災后恢復重建項目”以及臺站優化改造、標準化建設項目的相繼完成，成都站已建成一個基礎設施完善、觀測設備現代化、科研技術力量強勁的綜合臺站。成都站測震數據以及地磁數據均參與國際資料交換，服務于全球地震定位及地磁場空間變化研究。成都和松潘測震觀測系統運行穩定，觀測資料記錄連續，震相分析準確，后續震相豐富，其觀測資料精度和分析水平在國內始終處于前列，為中國地震臺網以及地震數據庫提供了大量、可靠的波形資料和震相分析資料，為地震監測、科學研究及防震減災事業作出了應有的貢獻。成都臺作為離5·12汶川地震最近的國家級臺站，擁有一套周期360 s的超寬頻帶地震計，所記錄的寶貴地震資料被廣泛地應用于國內外對地震的科學研究。

3 方法的MATLAB實現

第一步，讀取指定文件名的地震波形數據，確定參數，包括采樣頻率、地震計周期、數采轉換因子。

第二步，調用函數y=fft（x）實現快速傅里葉變換。

第三步，調用函數z=emd（x）對數據進行EMD分解，求得邊際譜和每個IMF分量及最后一個剩余分量residual的圖形。

第四步，用xt=x×C/2 000將count值轉化為速度值，其中C為數采轉化因子。然后調用mean()函數，用公式xt=xt-mean（xt）去均值。

第五步，調用pwelch()函數實現Welch算法，即：

式（1）中：pxx為功率譜估計值；f為得到的頻率點；x為進行功率譜估計的有限長序列，即上述已去均值的xt；window指所用窗函數；noverlap指重疊的點數；Nfft指采用的FFT算法長度；Fs為采樣頻率。

第六步，用函數sys=tf（[100]，[12×0.707×pi/T/2pi×pi/（T×T）/4]）求出頻率響應，再由函數[m p]=bode（sys，2×pi×f）求出幅值m，經過矩陣轉置后，由公式pxx=pxx/abs（mag^2）扣除儀器響應，再用繪出速度功率譜。其中T為地震計周期。

第七步，用db4小波對原始信號進行5層分解并提取小波系數：[c，l]=wavedec（x，5，'db4'），再用appcoef()和detcoef()求得低頻與高頻部分，隨后進行小波消噪，之后再進行高通濾波，本文所選擇的是IIR高通濾波器。

4 應用實例

本文截取了2016-07-15T02：00記錄于成都地震監測中心站JCZ-1T的垂直向波形數據，時間長度為300 s，包含了一個地震事件。

4.1 分析部分

由于原始波形受到嚴重的干擾，地震信息被干擾掩蓋難以分辨，如圖1所示。快速傅里葉變換（FFT）分析后的結果如圖2所示，由圖2可知，波形數據中存在低頻、高頻的干擾，2～4 Hz的頻率區間幅度較大。在Hilbert邊際譜上，同樣存在2～4 Hz頻率區間的較大幅度反應。這是由于2～4 Hz頻率范圍的干擾已成為成都地震臺（現成都地震監測中心站）測震資料中的固有干擾[2]。

為了驗證分析的結果，設計了FIR帶阻濾波器對2～4 Hz的頻率區間進行濾波，毛刺干擾被壓制，如圖3所示，證實了2～4 Hz頻率區間的較大幅度反應是波形中存在毛刺的原因。

4.2 處理部分

db4小波分解的結果如圖4所示，由圖4可知，低頻信號振幅與原信號更接近，而高頻信號比較集中且振幅較小，與原波形相差較大，去噪便是對高頻信號進行處理。將去噪后進行IIR高通濾波，得到最終處理結果，如圖5所示。由處理結果可知，整個處理過程中原始波形中的干擾被有效地壓制。

圖1 原始波形

圖2 FFT幅頻圖與Hilbert邊際譜

圖3 消除2～4 Hz頻率區間后的波形

圖4 db4小波分解后各層低頻、高頻信號

圖5 小波去噪及IIR高通濾波后的波形

5 認識與討論

本文所提出的測震資料分析處理方法在地震臺站的實際應用中，取得了良好的效果，對干擾進行了有效的壓制。文中詳細地介紹了測震數據分析處理方法的設計原理、工作流程，核心是干擾類別的識別。測震工作人員不僅局限于地震震相的劃分，還需要對數據中存在的干擾類型進行判斷并區分是固定干擾還是隨機干擾。

本文所提出的分析處理方法的初衷是加強對干擾類別的判斷與識別，以傅里葉變換、希爾伯特-黃變換以及速度功率譜作為分析判別的理論依據，小波變換和濾波器作為分析結果的驗證手段。未來將對方法進行軟件化，并希望在更多地震臺站進行試用。