前兆觀測干擾信號頻譜特征分析

李繼業(yè) 武曉軍 劉長生 任建輝 高研孫鵬宇 孫強(qiáng) 秦麗巖 戰(zhàn)玉明

1）黑龍江省地震局，哈爾濱市南崗區(qū)鴻翔路24號 150090

2）哈爾濱市防震減災(zāi)技術(shù)中心，哈爾濱 150021

0 引言

在數(shù)字信號中，對隨機(jī)信號的處理有重要意義。隨機(jī)信號又可分為平穩(wěn)隨機(jī)信號和非平穩(wěn)隨機(jī)信號，而現(xiàn)實中大多為非平穩(wěn)隨機(jī)信號，所謂非平穩(wěn)隨機(jī)信號，就是其統(tǒng)計特性是時間的函數(shù)。由于計算機(jī)技術(shù)的飛速發(fā)展與普及，非平穩(wěn)隨機(jī)信號分析與處理的應(yīng)用前景廣闊（Cohen，1989）。

過往的地震前兆研究表明，在觀測數(shù)據(jù)中一般都包含著動因復(fù)雜的各種信息，可將其分為正常動態(tài)及異常動態(tài)，而異常動態(tài)中又有與地震活動有關(guān)的震兆異常和與地震活動無關(guān)的干擾信號。地震前兆觀測數(shù)據(jù)由于不同程度地受到觀測站點的場地、自然環(huán)境、人為活動等各種因素的制約與干擾，其干擾信號出現(xiàn)的頻次往往比前兆異常信息要多得多。若干擾信號幅度大于正常觀測背景值時，地震活動前兆異常信息就難以提取。因此，如何識別和排除干擾，是提高地震監(jiān)測研究水平的一項非常重要的基礎(chǔ)性工作（邱永平，2011）。

前兆觀測數(shù)據(jù)可視為不同頻率、周期成分的疊加。數(shù)字化觀測中各種干擾信號與震兆異常信息混雜，有些干擾信號無論是形態(tài)上還是變幅上，都極易與震兆異常信息混淆，從而給數(shù)字化數(shù)據(jù)的使用和地震預(yù)測研究帶來很大的困難。本文選用基于傅里葉變換的傅里葉譜、功率譜以及基于希爾伯特-黃變換的邊際譜方法，將不同周期成分分離，形成不同頻率隨時間變化的多個時間系列，分析頻譜特征，進(jìn)而提取黑龍江數(shù)字化水位和豎直擺傾斜測項常見干擾類型，初步總結(jié)了水位高頻干擾、周期干擾以及豎直擺傾斜的風(fēng)擾、同震應(yīng)變等常見干擾信號的頻譜特征，以期探索出用譜分析方法提取前兆信息的有效途徑。

1 方法簡述

1.1 基于希爾伯特-黃變換的邊際譜方法

希爾伯特-黃變換（Hilbert-Huang Transform，簡稱 HHT）由 Huang等（1998、1999、2003）提出，是一種運用于非線性、非平穩(wěn)信號的處理方法，是對傅里葉變換基本信號和頻率定義做了創(chuàng)造性改進(jìn)，其不再認(rèn)為信號是由一組固定頻率和振幅的正弦信號的加權(quán)和，而是由固有模態(tài)函數(shù)（intrinsic mode function，簡稱IMF）信號所組成。希爾伯特-黃變換方法由經(jīng)驗?zāi)B(tài)分解（empiricalmode decomposition，簡寫為 EMD）（Rilling et al，2003、2007）方法和 Hilbert變換組成。EMD方法即Huang變換，它依據(jù)信號本身時間尺度特征，將信號分解為含有不同時間尺度且滿足以下2個定義條件的一組IMF，每個IMF可以被認(rèn)為是信號中固有的一個模態(tài)函數(shù)：①對于一列數(shù)據(jù)極值點和過零點數(shù)目必須相等或至多相差一點；②在任意點，由局部極大點和極小點構(gòu)成的2條包絡(luò)線的平均值為零。

EMD分解過程：①找出信號x（t）的所有極大點和極小點，將其用3次樣條函數(shù)分別擬合為原數(shù)據(jù)序列的上、下包絡(luò)線，上、下包絡(luò)線的均值為平均包絡(luò)線m1，將原數(shù)據(jù)序列減去m1可得到一個去掉低頻的新數(shù)據(jù)序列h1。一般h1不是一個平穩(wěn)數(shù)據(jù)序列，為此重復(fù)以上過程n次，使所得到的平均包絡(luò)線趨于零，此時的h1n就是第一個IMF（c1），它表示信號數(shù)據(jù)中的最高頻成分。②用x（t）減去c1得到一個去掉高頻成分的新數(shù)據(jù)序列，重復(fù)步驟①，得到一系列cn和最后一個不可分解的序列rn，rn代表x（t）的均值或趨勢項。那么原序列x（t）可表示為IMF分量和一個殘余項的和

信號經(jīng)分解后得到多個IMF分量，對每個IMF分量進(jìn)行Hilbert變換，即可得到每個IMF分量的瞬時頻率，綜合所有IMF分量的瞬時頻譜得到Hilbert譜。

獲得IMF分量以后，對每一階IMF作 Hilbert變換，設(shè)為y（t）

x（t）和y（t）共同組合為一解析信號z（t），采用極坐標(biāo)形式表示

定義瞬時頻率

由式（6）可看出，ω（t）是時間t的單值函數(shù)，即某一時間對應(yīng)某一頻率，每個IMF序列在每一點的頻率唯一。對每一階IMF作Hilbert變換，并求出相應(yīng)解析函數(shù)的幅值譜和瞬時頻率，從而原始信號可以表示為

上式的頻率ωj（t）和幅值aj（t）都是時間的函數(shù)，可構(gòu)成幅值、頻率、時間的3維時頻圖即Hilbert幅值譜，簡稱Hilbert譜，記為

基于Hilbert譜公式（8），通過對時間的積分，可以定義邊際譜h（ω）為

邊際譜表達(dá)了每個頻率在全局上的幅度，它代表了在統(tǒng)計意義上的全部累加幅度。

1.2 基于傅里葉變換的傅里葉譜和功率譜

快速傅里葉變換（Fast Fourier Transform，簡稱FFT），是離散傅里葉變換的快速算法，也可用于計算離散傅里葉變換的逆變換。它是將時間域的信號變換成頻率域的信號，得到與原信號長度相同的時間序列（萬永革等，2009），因為有些信號在時域上很難看出什么特征，但當(dāng)變換到頻域之后，就完全不同了。FFT變換還可以將一個信號的頻譜提取出來，這常用于頻譜分析。由于隨機(jī)信號是一類持續(xù)時間無限長，具有無限大能量的功率信號，它不滿足傅里葉變換條件，而且也不存在解析表達(dá)式，因此就不能應(yīng)用確定信號的頻譜計算方法去分析隨機(jī)信號的頻譜。這樣就可以用功率譜來描述隨機(jī)信號的頻域特性，即信號能量特征隨頻率的變化。為行文簡潔，傅里葉變換的具體算法、公式不再贅述。

2 干擾頻譜特征

過往的地下水微動態(tài)和地傾斜研究表明，井水位和鉆孔傾斜對地殼應(yīng)力-應(yīng)變的響應(yīng)頻率較寬，可記錄到地震波中的瑞雷波（10～20s）、潮汐波動（?3～?2 d）、氣壓波動（幾小時～幾天）、風(fēng)擾波動（十幾～幾十分鐘）。除此之外，井水位和地傾斜觀測還受降雨效應(yīng)（補給、荷載）、地下水開采和觀測儀器故障等干擾因素的影響（張淑亮等，2005；薄萬舉，2010）。因此，排除干擾信號是前兆觀測和數(shù)據(jù)分析的重要前提。由于篇幅關(guān)系，本文僅例舉具有普遍性、代表性的水位高頻、周期性干擾以及豎直擺傾斜的風(fēng)擾、同震應(yīng)變等干擾信號的頻譜分析結(jié)果。

在實際觀測中，甘南臺水位存在明顯的周期性干擾信號。通過對甘南臺水位2013年2月17日0時1分～8時59分的540個分鐘值干擾信號去趨勢（圖1（a））后進(jìn)行頻譜特征分析發(fā)現(xiàn)，傅里葉譜和功率譜周期存在2.2～10min多個頻譜值（圖1（c）、（d），表1），而邊際譜周期存在9～18min多個頻譜值（圖1（b）、表1）。傅里葉譜和功率譜譜峰分布在整個頻域區(qū)間，并有多條等間距的譜峰出現(xiàn)，而邊際譜卻不明顯，說明影響甘南水位的信號應(yīng)該是一個固定周期的信號。可見，傅里葉譜和功率譜識別水位周期干擾信號效果更好，體現(xiàn)了傅立葉變換在平穩(wěn)信號處理中的優(yōu)勢。

圖1 典型的水位周期性干擾頻譜特征

圖2 典型的水位高頻干擾信號的頻譜特征

表1 干擾信號的頻譜特征

在實際觀測中，孫吳臺水位存在明顯的漏電高頻干擾信號。通過對孫吳臺水位2009年10月1日0時1分～8時59分的540個分鐘值干擾信號去趨勢（圖2（a））后進(jìn)行頻譜特征分析發(fā)現(xiàn)，傅里葉譜和功率譜周期存在2～20min多個頻譜值（圖2（c）、（d），表1），而邊際譜周期存在13～20min多個頻譜值（圖2（b）、表1）。傅里葉譜和功率譜譜峰分布在整個頻域區(qū)間，而邊際譜卻不明顯，說明影響孫吳臺水位的高頻干擾信號亦是一個有固定周期的信號，屬于平穩(wěn)信號。可見，傅里葉譜和功率譜識別水位高頻干擾信號效果更好，體現(xiàn)了傅立葉變換在平穩(wěn)信號處理中的優(yōu)勢。

形變觀測受多種因素影響，包括氣象變化、人為影響、儀器故障等因素，使排除干擾成為識別異常必不可少的工作。在目前認(rèn)識水平下，干擾形態(tài)與地震前兆形態(tài)極為相似，給干擾與前兆的判斷帶來很大困難。通河臺豎直擺傾斜存在明顯的季節(jié)性大風(fēng)干擾，通過對通河臺豎直擺傾斜2012年8月29日14時1分～23時59分NS、EW2個分量各600個分鐘值干擾信號去趨勢（圖3a、4a）后進(jìn)行頻譜特征分析發(fā)現(xiàn)，傅里葉譜和功率譜對風(fēng)擾信號頻譜特征提取效果較好，而邊際譜提取效果并不明顯。通河臺豎直擺傾斜NS測項傅里葉譜和功率譜周期存在 11～26min多個頻譜值（圖 3（c）、（d），表 1），而邊際譜頻譜周期存在 8～27min多個譜值（圖3（b）、表1）。EW測項傅里葉譜和功率譜周期存在 3.7～18min多個頻譜值（圖 4（c）、（d），表 1），邊際譜周期存在 4～26min多個頻譜值（圖 4（b）、表 1）。功率譜和傅里葉譜出現(xiàn)的頻譜周期均小于27min，邊際譜存在77min的獨立周期?？梢?，傅里葉譜和功率譜識別傾斜風(fēng)擾信號效果更好，體現(xiàn)了傅立葉變換在平穩(wěn)信號處理中的優(yōu)勢。

圖3 典型的豎直擺傾斜風(fēng)擾信號的頻譜特征（NS）

圖4 典型的豎直擺傾斜風(fēng)擾信號的頻譜特征（EW）

在傾斜儀記錄中，由于地震引起地表變形使觀測曲線偏離震前的位置，這種突變稱傾斜同震應(yīng)變。當(dāng)?shù)卣鸢l(fā)生時，孕震體釋放能量，震源區(qū)應(yīng)變場發(fā)生調(diào)整，由于地震動導(dǎo)致的應(yīng)力-應(yīng)變傳遞，使得各地應(yīng)變場也隨之調(diào)整，在震源區(qū)外的儀器能夠記錄到這種直接調(diào)整，這就是為什么在上千千米范圍內(nèi)的儀器都能觀測到同震應(yīng)變的原因（陳德福，1993）。前兆觀測儀器記錄到的同震應(yīng)變對于正常背景觀測數(shù)據(jù)而言也可以看作是一種干擾信號。

2013年 4月 6日 12時 42分，印度尼西亞（138.5°E、3.5°S）發(fā)生M7.0強(qiáng)震，震源深度70km，在震中距5560km的通河臺豎直擺傾斜儀記錄到了明顯的同震應(yīng)變波。通過對通河豎直擺傾斜2013年4月6日12時1分～16時59分的NS、EW2個分量各300個分鐘值同震應(yīng)變信號去趨勢（圖5（a）、6（a））后進(jìn)行頻譜特征分析發(fā)現(xiàn)，傅里葉譜和功率譜識別同震應(yīng)變的頻譜特征較為一致，其中NS測項傅里葉譜和功率譜周期存在2.9～12min多個頻譜值（圖 5（c）、（d），表 1），而邊際譜周期存在 2.8～14min多個頻譜值（圖 5（b）、表 1）。EW測項傅里葉譜和功率譜周期存在 2.1～9min多個頻譜值（圖6（c）、（d），表1），邊際譜周期存在2.2～14min多個頻譜值（圖6（b），表1）?？梢?，邊際譜較傅里葉譜和功率譜識別同震應(yīng)變信號效果更好，體現(xiàn)了邊際譜在非平穩(wěn)信號處理中的優(yōu)勢。

通過對黑龍江數(shù)字化水位周期、高頻干擾以及豎直擺傾斜等的風(fēng)擾、同震應(yīng)變干擾進(jìn)行頻譜分析發(fā)現(xiàn)，甘南臺水位和孫吳臺水位存在固定周期的干擾信號疊加在正常背景觀測信號中，且分布在整個頻域內(nèi)；而通河臺豎直擺傾斜風(fēng)擾和同震應(yīng)變干擾頻譜周期主要集中在10min以下。傅里葉譜和功率譜的頻譜特征基本一致，這與它們均基于傅立葉變換有關(guān)，且功率譜識別干擾信號效果更好。由于邊際譜體現(xiàn)了某一段頻率總的幅值累積特征，邊際譜譜值較傅里葉譜、功率譜譜值要高出幾個數(shù)量級；同時邊際譜存在77min的獨立周期，這可能是希爾伯特-黃變換的虛假波動現(xiàn)象。

圖5 豎直擺傾斜同震應(yīng)變的頻譜特征（NS）

圖6 豎直擺傾斜同震應(yīng)變的頻譜特征（EW）

3 討論與結(jié)論

前兆觀測及研究表明，與現(xiàn)有的儀器觀測精度和干擾水平相比，多數(shù)的地震前兆信息屬于弱信息，接近噪聲水平。為了更好地、真實地、連續(xù)可靠地提取地震前兆信息，對干擾信號進(jìn)行合理的提取和分析，是非常關(guān)鍵和重要的環(huán)節(jié)。

通過采用基于傅里葉變換的功率譜、傅里葉譜以及基于希爾伯特-黃變換的邊際譜方法，對黑龍江數(shù)字化水位和豎直擺傾斜典型干擾數(shù)據(jù)進(jìn)行分析，提取水位周期干擾、高頻干擾以及豎直擺傾斜的風(fēng)擾、同震應(yīng)變等常見干擾信號的頻譜特征。通過對水位、豎直擺傾斜典型干擾信號進(jìn)行頻譜分析，發(fā)現(xiàn)干擾信號優(yōu)勢頻率在1.7×10-3×10-3Hz以上，周期一般在10min以下。其中，水位周期性干擾信號優(yōu)勢頻率為1.72×10-3～7.61×10-3Hz，周期一般為2～10min，高頻干擾信號優(yōu)勢頻率為 1.89×10-3～7.01×10-3Hz，周期一般為 2～9min；豎直擺傾斜大風(fēng)干擾信號優(yōu)勢頻率為1.51×10-3～4.23×10-3Hz，周期一般為 4～11min，豎直擺傾斜同震應(yīng)變干擾信號優(yōu)勢頻率為2.08×10-3～7.51×10-3Hz，周期一般為2～12min。

頻譜變換方法在前兆觀測資料干擾信號提取中各有優(yōu)勢。傅里葉變換可以將任何信號分解為正弦信號的加權(quán)和，而每一個正弦信號對應(yīng)著一個固定的頻率和固定的幅值。希爾伯特-黃變換對于處理非線性、非平穩(wěn)信號有清晰的物理意義，能夠得到信號的振幅-時間-頻率分布特征，能使信號分解具有唯一性，又能在時域和頻域中同時具有良好的局部化性質(zhì)。由于希爾伯特-黃變換突破傅里葉變換僅對平穩(wěn)信號有效的不足之處，使其能夠更好的應(yīng)用于非線性、非平穩(wěn)信號的處理。研究發(fā)現(xiàn)，功率譜和傅里葉譜方法對提取干擾和噪聲等線性、平穩(wěn)信號更加有效，而邊際譜對提取同震應(yīng)變等非線性、非平穩(wěn)信號更加顯著，豐富了頻譜分析方法在前兆觀測信息提取中的應(yīng)用。