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

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

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

0 引言

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

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

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

1 方法簡(jiǎn)述

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

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

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

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

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

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

定義瞬時(shí)頻率

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

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

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

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

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

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

2 干擾頻譜特征

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

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

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

圖2 典型的水位高頻干擾信號(hào)的頻譜特征

表1 干擾信號(hào)的頻譜特征

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

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

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

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

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

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

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

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

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

3 討論與結(jié)論

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

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

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