何思源,李貴元,劉華姣,鐘李彬,徐文海
(四川省地震局成都地震基準(zhǔn)臺(tái),四川 成都 611730)
地震觀測(cè)是研究地震和地球內(nèi)部結(jié)構(gòu)的重要手段,也是地震學(xué)理論體系建立和發(fā)展的重要基礎(chǔ),數(shù)字地震記錄已成為現(xiàn)代地震研究的基石(李興泉等,2018)。目前,MSDP軟件廣泛應(yīng)用于臺(tái)站,該軟件提供了波形分析功能,缺少波形干擾處理功能,無(wú)法滿足實(shí)際工作需求。雖然通過(guò)濾波程序可以達(dá)到濾掉部分干擾的效果,但整體缺乏綜合分析處理能力,難以推廣應(yīng)用。地震波形中包含了眾多的信息,其中不乏各種干擾,這些干擾對(duì)數(shù)字地震記錄的質(zhì)量造成了影響,便會(huì)對(duì)于震相的識(shí)別與分析造成干擾,為測(cè)震分析工作增加難度和工作量,因此如何快速有效地完成對(duì)干擾的去除就顯得十分必要。本研究基于MATLAB開(kāi)發(fā)了一款測(cè)震資料干擾分析處理軟件,其內(nèi)部不僅集成了快速傅里葉變換、希爾伯特-黃變換、小波變換、速度功率譜估計(jì)這些分析算法,還包括了不同類型的去噪方法,可以滿足裝配速度型地震計(jì)臺(tái)站觀測(cè)資料的干擾分析與處理需求。
如何去除地震波形中的干擾信息,保留其有效信息,首先需要了解有效信息和干擾的特點(diǎn),才能正確地將二者加以區(qū)分。傅里葉變換、希爾伯特-黃變換、小波變換和速度功率譜估計(jì)則是整個(gè)“波形分析”功能的基礎(chǔ),整體流程如圖1所示?!安ㄐ翁幚怼辈捎玫氖荈IR濾波器和小波去噪兩種方法,F(xiàn)IR濾波器由于具有線性相位的特性而被廣泛使用,小波變換去噪實(shí)質(zhì)是軟閾值去噪法。閾值法去噪的優(yōu)點(diǎn)是噪聲幾乎完全得到有效抑制,且反映原始信號(hào)的特征尖峰點(diǎn)也得到很好的保留,能夠使估計(jì)信號(hào)實(shí)現(xiàn)最大,均方誤差最小化,即去噪后的估計(jì)信號(hào)是原始信號(hào)的近似最優(yōu)估計(jì)。小波變換具有時(shí)頻局部化特性、多分辨特性、去相關(guān)特性和選基靈活性的特點(diǎn)(周宇峰等,2008)。前兩個(gè)性質(zhì)決定了小波變換去噪相比傳統(tǒng)方法具有獨(dú)特的優(yōu)勢(shì),是信號(hào)處理的得力工具。在利用小波變換分析和處理信號(hào)時(shí),不同領(lǐng)域的學(xué)者對(duì)如何選取合適的小波基有不盡相同的看法,而且對(duì)同一信號(hào)選取不同的小波基進(jìn)行處理時(shí)會(huì)出現(xiàn)完全不同的結(jié)果(張華等,2011)。因此,需要根據(jù)信號(hào)的特點(diǎn),采取折衷的辦法進(jìn)行小波基的選取(高玉寶等,2008)。軟件中給出了一些常用類型的小波基供用戶進(jìn)行選擇,如圖2所示。
圖1 波形分析與處理總流程
圖2 軟件所提供的小波基類型
選取成都地震臺(tái)JCZ-1T地震計(jì)記錄的2016年7月27日19時(shí)59分四川平武2.5級(jí)地震垂直分向波形數(shù)據(jù),見(jiàn)如圖3a。選擇“波形分析”功能中的“頻譜分析”功能,可以分別得到FFT頻譜圖和希爾伯特邊際譜,如圖3b所示。“Welch速度功率譜”可以根據(jù)設(shè)置不同的參數(shù)(FFT長(zhǎng)度、窗函數(shù)長(zhǎng)度、重采樣精度)得到相應(yīng)的速度功率譜估計(jì)圖,如圖3c。根據(jù)FFT頻譜圖及希爾伯特邊際譜可知,在2~4 Hz的頻率區(qū)間內(nèi)存在較高的幅值,而速度功率譜估計(jì)也在此區(qū)間出現(xiàn)了高幅值。
圖3 波形分析功能示例
為了研究該2~4 Hz的頻率區(qū)間所表示的內(nèi)容,文中首先采用FIR濾波器對(duì)地震波形中的2~4 Hz頻率范圍進(jìn)行帶阻濾波,以此研究2~4 Hz頻率范圍是否與地震波形中的干擾有關(guān)。選擇軟件中的“帶阻濾波器”功能,通過(guò)設(shè)置通、阻帶頻率,便可以實(shí)現(xiàn)帶阻濾波,如圖4所示。濾波前后的對(duì)比結(jié)果顯示:濾除2~4 Hz頻率范圍的干擾后,干擾被壓制的效果十分明顯,原始波形中的“毛刺”干擾也得到了有效的改善,表明濾除的2~4 Hz頻段波正是干擾所在的頻率范圍。
圖4 帶阻濾波前后波形
選取成都地震臺(tái)JCZ-1T地震計(jì)記錄的2017年12月10日08點(diǎn)14分九寨溝2.9級(jí)地震記錄作為處理實(shí)例。原始波形自身存在明顯的波形“毛刺”,同時(shí)波形中還存在高頻干擾信息,可以采用小波去噪和2~4 Hz帶阻濾波法綜合去除噪聲的辦法進(jìn)行處理,其中2~4 Hz帶阻濾波法是為了消除測(cè)震資料中的日常普遍性干擾(何思源等,2019);再將處理的結(jié)果進(jìn)行仿真,如圖5所示。其中,小波去噪法選擇的小波基為bior2.4小波,因?yàn)閎ior2.4小波基重構(gòu)精度高,去噪效果良好,能夠更加完整地保留了有效信號(hào)(何思源等,2019),其分解層數(shù)為4。原始波形經(jīng)過(guò)處理后, “毛刺”及高頻干擾得到了有效的壓制,仿真后波形的垂直分向上原本被干擾所掩蓋的初至震相變得清晰可辨,信噪比得到了明顯的提升。
圖5 去噪前后波形及仿真結(jié)果
由于經(jīng)濟(jì)發(fā)展、城市規(guī)模擴(kuò)大以及交通建設(shè)等對(duì)部分臺(tái)站的觀測(cè)條件產(chǎn)生了影響。瀘州地震臺(tái)由于地處市區(qū)且長(zhǎng)期受到周邊城建施工的影響,臺(tái)站觀測(cè)數(shù)據(jù)整體呈“晝差夜良”的特點(diǎn),如圖6所示。瀘州地震臺(tái)的觀測(cè)數(shù)據(jù)質(zhì)量保持“晝夜一致”對(duì)提高該地震臺(tái)測(cè)震數(shù)據(jù)質(zhì)量有著積極的意義。分別選取一段2019年3月12日瀘州地震臺(tái)白天和夜晚記錄的波形數(shù)據(jù)進(jìn)行波形分析。通過(guò)“波形分析”功能中的“頻譜分析”可以分別得到白天和夜晚觀測(cè)數(shù)據(jù)的FFT頻譜圖,白天相比夜晚主要在高頻段具有更高的幅值水平(見(jiàn)圖7),再采用小波去噪與濾波器相結(jié)合的方式進(jìn)行處理。結(jié)果表明,白天的觀測(cè)數(shù)據(jù)去噪處理后,波形質(zhì)量得到明顯改善,能夠與夜晚的觀測(cè)數(shù)據(jù)質(zhì)量保持一致,如圖8所示。
圖6 瀘州地震臺(tái)波形數(shù)據(jù)
圖7 瀘州地震臺(tái)白天及夜晚波形數(shù)據(jù)FFT頻譜圖
圖8 瀘州地震臺(tái)去噪前后波形質(zhì)量對(duì)比
此項(xiàng)功能的作用是通過(guò)記錄用戶使用帶阻濾波器濾除的頻率區(qū)間,當(dāng)某一頻率區(qū)間累計(jì)出現(xiàn)了20次,那么軟件會(huì)自動(dòng)將其識(shí)別為“疑似普遍性干擾”,并及時(shí)提醒用戶注意此頻率區(qū)間信息的干擾。該功能的原理是用戶在第一次使用“帶阻濾波器”功能時(shí),會(huì)在D盤(pán)自動(dòng)建立一個(gè)文件夾且命名為“濾波器結(jié)果文件”,里面會(huì)自動(dòng)生成一個(gè)名為“帶阻濾波”的Excel表格并記錄下本次所選擇濾除的頻率區(qū)間。之后用戶每使用一次“帶阻濾波器”功能,軟件都會(huì)對(duì)選擇處理的頻率區(qū)間進(jìn)行記錄并統(tǒng)計(jì)出現(xiàn)的次數(shù)。當(dāng)有一組頻率區(qū)間第20次出現(xiàn)在Excel表格中時(shí),軟件就會(huì)彈出消息提醒框,提示用戶注意該頻率區(qū)間。例如,軟件在進(jìn)行帶阻濾波處理時(shí),用戶有20次選擇濾除2~4 Hz的頻率區(qū)間,那么軟件便會(huì)彈出提示信息,如圖9所示。
圖9 “疑似普遍性干擾頻率”提示功能示例
在測(cè)震資料干擾分析處理軟件的開(kāi)發(fā)和使用過(guò)程中,有以下幾點(diǎn)認(rèn)識(shí):
(1)軟件的核心是“波形分析”與“波形處理”功能?!安ㄐ畏治觥睂⒏道锶~變換、希爾伯特-黃變換、小波變換和速度功率譜估計(jì)作為基礎(chǔ),對(duì)于測(cè)震資料中的干擾能夠?qū)崿F(xiàn)時(shí)域、頻率的分析;“波形處理”功能以FIR濾波器及小波去噪為手段,能夠去一般的干擾信號(hào)。(2)“疑似普遍干擾頻率”提示功能作為附加功能,可以及時(shí)地為用戶反饋測(cè)震資料中可能存在的普遍性干擾,提示用戶對(duì)干擾產(chǎn)生的原因進(jìn)行核實(shí)。(3)在實(shí)際應(yīng)用中分析地震波形中存在的干擾,并能進(jìn)行有效的處理,改善波形數(shù)據(jù)的質(zhì)量。(4)經(jīng)過(guò)分析發(fā)現(xiàn),成都地震臺(tái)測(cè)震資料中普遍存在頻率區(qū)間為2~4 Hz的干擾,需要對(duì)該頻段可疑干擾源進(jìn)行核實(shí)。