基于小波分析對(duì)形變擾動(dòng)信號(hào)的觸發(fā)檢測(cè)方法

劉明輝 李 江 沙成寧 周銀興 陳 陽(yáng) 崔仁勝 朱小毅 林 湛

1 中國(guó)地震局地震預(yù)測(cè)研究所地震預(yù)測(cè)重點(diǎn)實(shí)驗(yàn)室,北京市復(fù)興路63號(hào),100036 2 青海省地震局,西寧市興海路1號(hào),810001

震源區(qū)應(yīng)力場(chǎng)在地震前的短臨發(fā)震階段處于變動(dòng)中,越接近發(fā)震,發(fā)震斷層在斷裂前擴(kuò)展具有的脈沖性、間隙性表現(xiàn)越突出,將引起以突跳性、突發(fā)性和幾起幾落為主要特征的前兆現(xiàn)象[1],常被稱為前驅(qū)波、慢地震、靜地震等。這種前兆現(xiàn)象一般在震前7 d至地震發(fā)生前出現(xiàn),周期為幾十秒到幾十分鐘[2-3],可能會(huì)被鉆孔應(yīng)變儀、重力儀、傾斜儀、超長(zhǎng)周期地震儀等儀器記錄下來(lái)[4-7],但儀器同時(shí)也會(huì)記錄到地球固體潮和干擾信號(hào)。地震前驅(qū)波和干擾信號(hào)相對(duì)于地球固體潮來(lái)說(shuō)周期短、能量弱,疊加在固體潮記錄上形成顫抖擾動(dòng),不易直觀提取和判別,給地震研究帶來(lái)較大困難。

鉆孔應(yīng)變儀的記錄頻帶范圍為DC～10 Hz,能夠記錄到清晰完整的固體潮[8]。通過(guò)計(jì)算鉆孔應(yīng)變觀測(cè)的4個(gè)分量之間是否滿足自洽方程,可檢查觀測(cè)數(shù)據(jù)的穩(wěn)定性和可靠性[9]。小波分析是地殼形變資料分析中一種常用的方法[10],STA/LTA觸發(fā)檢測(cè)算法是一種經(jīng)典的大信號(hào)檢測(cè)方法,廣泛應(yīng)用于日常地震監(jiān)測(cè)、地震烈度速報(bào)、地震預(yù)警和地震科學(xué)研究[11-12]。本文以2008年汶川8.0級(jí)地震前玉樹(shù)地震臺(tái)鉆孔應(yīng)變觀測(cè)數(shù)據(jù)為例,用小波分析從鉆孔應(yīng)變觀測(cè)數(shù)據(jù)中提取地震前的擾動(dòng)信號(hào),用STA/LTA觸發(fā)檢測(cè)算法進(jìn)行擾動(dòng)信號(hào)觸發(fā)檢測(cè)。

1 數(shù)據(jù)處理方法

1.1 自洽檢驗(yàn)

鉆孔應(yīng)變儀的4個(gè)觀測(cè)分量(S1、S2、S3和S4)按照順時(shí)針?lè)较蛳嗖?5°依次分布,如圖1所示。其中,S1與S3、S2與S4之間呈90°正交,分別合成分量S13、S24,自洽方程為[9]:

圖1 鉆孔應(yīng)變觀測(cè)分量分布示意圖Fig.1 Observation components distribution sketch of strainmeter

S1+S3=S2+S4

(1)

式(1)成立時(shí),4個(gè)分量的觀測(cè)數(shù)據(jù)可靠,能真實(shí)地反映地層應(yīng)變變化。本文通過(guò)計(jì)算S13與S24之間的相關(guān)系數(shù)進(jìn)行自洽檢驗(yàn)。

1.2 相干分析

利用相干分析得到S13和S24信號(hào)頻率相干程度。2個(gè)信號(hào)x(t)和y(t)間相干函數(shù)為[13]:

(2)

式中,Pxy(ω)為x(t)和y(t)的互功率譜密度,Pxx(ω)、Pyy(ω)分別為x(t)和y(t)的頻率域內(nèi)的自功率譜密度。Cxy(ω)的取值范圍為[0,1],若y(t)與x(t)完全不相關(guān),則Cxy(ω)=0;若y(t)為x(t)的線性響應(yīng),則Cxy(ω)=1。

1.3 小波分析

小波分析將信號(hào)分解到尺度域,通過(guò)多分辨率的分解使原始信號(hào)中的弱信號(hào)成分變得突出,具有優(yōu)良的時(shí)頻局部能力。給定一個(gè)基本函數(shù)[14]:

(3)

式中,a、b為常數(shù),a>0,φa,b(t)由基本函數(shù)φ(t)先移位再伸縮后得到,若a、b不斷變化,可得到一組φa,b(t)。信號(hào)x(t)的小波分析為:

(4)

式中,*代表共軛,a、b和t為連續(xù)變量,式(4)的頻率域?yàn)?

WTx(a,b)=

(5)

式中,X(Ω)為x(t)的傅里葉變換,如果φ(t)是幅頻特性比較集中的帶通函數(shù),則小波分析就具有表征待分析信號(hào)在頻域上局部性質(zhì)的能力。當(dāng)a較小時(shí),時(shí)域觀測(cè)范圍較小,在頻域上相當(dāng)于用高頻小波進(jìn)行細(xì)致觀察;當(dāng)a較大時(shí),時(shí)域觀察范圍較大,在頻率域上相當(dāng)于用低頻小波進(jìn)行概貌觀察。這為選擇小波分析的參數(shù)及提取擾動(dòng)信號(hào)的目標(biāo)層提供了參考。

1.4 時(shí)間-頻率-能量分析

對(duì)小波分析提取出來(lái)的擾動(dòng)信號(hào)進(jìn)行時(shí)間-頻率-能量分析,得到時(shí)間-頻率變化關(guān)系及能量變化趨勢(shì)。短時(shí)傅里葉變換以固定的滑動(dòng)窗口對(duì)信號(hào)進(jìn)行分析,隨著窗函數(shù)的滑動(dòng),可表征信號(hào)的局域頻率特性。令g(t)為時(shí)間寬度很窄的窗函數(shù),并沿時(shí)間軸滑動(dòng),信號(hào)z(t)的短時(shí)傅里葉變換(STFT)定義為[15]:

(6)

由于信號(hào)z(t)乘以一個(gè)相當(dāng)窄的窗函數(shù)g(u-t)等價(jià)于取出信號(hào)在分析時(shí)間點(diǎn)t附近的一個(gè)切片,所以STFT(t,f)可以理解為信號(hào)z(t)在分析時(shí)間t附近的局部頻譜。

1.5 STA/LTA檢測(cè)算法

STA/LTA檢測(cè)算法中,STA是短時(shí)間窗口內(nèi)信號(hào)幅度的絕對(duì)均值,目的是捕獲突發(fā)大信號(hào);LTA是長(zhǎng)時(shí)間窗口內(nèi)信號(hào)幅度的絕對(duì)均值,目的是得到突發(fā)大信號(hào)到達(dá)之前的平均背景噪聲水平,作為突發(fā)大信號(hào)觸發(fā)檢測(cè)的參考基礎(chǔ)。STA/LTA算法的步驟為:首先根據(jù)檢測(cè)目標(biāo)的需要選擇適當(dāng)長(zhǎng)度的短時(shí)間和長(zhǎng)時(shí)間窗口,分別計(jì)算對(duì)應(yīng)時(shí)間窗口內(nèi)信號(hào)的絕對(duì)均值STA和LTA;其次計(jì)算這2個(gè)均值的比值,即STA/LTA;最后將比值與啟動(dòng)觸發(fā)閾值進(jìn)行比較。若比值大于啟動(dòng)觸發(fā)閾值,認(rèn)為檢測(cè)到大信號(hào);當(dāng)比值低于停止觸發(fā)閾值時(shí),認(rèn)為大信號(hào)結(jié)束。

2 數(shù)據(jù)處理結(jié)果

2.1 數(shù)據(jù)及其自洽檢查

2008年汶川8.0級(jí)地震前,玉樹(shù)地震臺(tái)鉆孔應(yīng)變儀記錄到的同震效應(yīng)信號(hào)幅度非常大。因此本文以汶川地震前5 d的記錄數(shù)據(jù)為例,用小波分析提取其擾動(dòng)信號(hào),進(jìn)行儀器觸發(fā)檢測(cè)。利用Tsoft軟件對(duì)玉樹(shù)地震臺(tái)鉆孔應(yīng)變記錄到的4個(gè)分量進(jìn)行去零點(diǎn)漂移等預(yù)處理[16],波形如圖2(a)所示,根據(jù)式(1)計(jì)算分量S13和S24,波形如圖2(b)所示。

圖2 玉樹(shù)地震臺(tái)鉆孔應(yīng)變儀記錄Fig.2 The records of borehole strainmeter at Yushu seismic station

計(jì)算可得,S13與S24的相關(guān)系數(shù)為0.999,滿足自洽方程,說(shuō)明玉樹(shù)地震臺(tái)的鉆孔應(yīng)變儀工作正常。由于固體潮的周期長(zhǎng)、能量強(qiáng),而擾動(dòng)信號(hào)的周期短、能量弱,從圖2(a)和2(b)中無(wú)法直接分辨出明顯的擾動(dòng)信號(hào)。

2.2 相干分析

根據(jù)式(2)計(jì)算合成分量S13和S24的相干系數(shù),結(jié)果如圖3所示。圖3顯示,S13和S24中不相干的頻率主要分布在120～130 s(2 min左右)和1 000～6 000 s(16～100 min),也是前驅(qū)波、慢地震、靜地震等的頻帶范圍[17-18]。S13和S24在固體潮頻段的信號(hào)是完全相干的,說(shuō)明玉樹(shù)臺(tái)鉆孔應(yīng)變儀工作穩(wěn)定可靠。

2.3 小波分析和時(shí)間-頻率-能量分析

小波分析將原始信號(hào)分成不同的層,使其含有不同頻率的信號(hào)成分[19]。研究表明,地震前兆異常的判據(jù)不僅包括幅度變化,還包括頻率變化等[20]。根據(jù)式(5)選擇的時(shí)域觀察范圍較大,而在頻率域上相當(dāng)于用低頻小波作概貌觀察。利用db小波分析原始數(shù)據(jù),將第1層信號(hào)提取出來(lái),并進(jìn)行時(shí)間-頻率-能量分析,結(jié)果如圖4所示。從圖4(a)中可以看出,50～60 h之間存在一個(gè)擾動(dòng)信號(hào),其在圖4(b)中形成了一個(gè)大峰值,在圖 4(c)中顯示的能量最強(qiáng),4個(gè)分量中尖峰的信號(hào)頻率各有不同,能量也各有不同,S2分量中的頻率多于其他3個(gè)分量。從圖4(a)中無(wú)法直接看到汶川8.0級(jí)地震前的擾動(dòng)信號(hào),但從圖4(b)和4(c)中可以看到地震前的擾動(dòng)信號(hào)。圖4(c)顯示,峰值之后信號(hào)能量逐漸增強(qiáng)。

為進(jìn)一步分析,將汶川地震前24 h的信號(hào)作歸一化處理后再進(jìn)行時(shí)間-頻率-能量分析,結(jié)果如圖5所示。從圖5(a)看出,汶川8.0級(jí)地震發(fā)生前2 h,4個(gè)分量提取的信號(hào)中都有明顯的擾動(dòng)成分,S1與S3的擾動(dòng)信號(hào)形態(tài)大致相同,S2與S4的擾動(dòng)信號(hào)形態(tài)大致相同,4個(gè)分量都具有間歇性的特征。圖5(b)顯示,不同頻率的擾動(dòng)信號(hào)的能量有差異。

圖5 歸一化原始信號(hào)、擾動(dòng)信號(hào)及擾動(dòng)信號(hào)時(shí)間-頻率-能量分布Fig.5 Normalized raw signal, disturbance signal and distribution of time-frequency-energy of disturbance signal

2.4 觸發(fā)檢測(cè)

STA/LTA觸發(fā)檢測(cè)算法有4個(gè)基本觸發(fā)參數(shù):STA和LTA的時(shí)間窗口長(zhǎng)度、啟動(dòng)和停止觸發(fā)閾值。圖4(c)和5(b)顯示,汶川8.0級(jí)地震的擾動(dòng)信號(hào)中能量較高的頻率范圍主要為2 000～6 000 s。將STA和LTA時(shí)間窗口長(zhǎng)度分別設(shè)置為1 000 s和9 600 s,啟動(dòng)和停止觸發(fā)閾值分別設(shè)置為1.6和0.8,對(duì)6個(gè)分量(S1、S2、S3、S4、S13和S24)中利用小波分析提取出來(lái)的擾動(dòng)信號(hào)進(jìn)行觸發(fā)檢測(cè)。結(jié)果表明,6個(gè)分量的檢測(cè)效果都很好。

汶川地震前5 d及12 h的S1分量觸發(fā)檢測(cè)效果分別如圖6(a)和6(b)所示,圖中灰色為提取的擾動(dòng)信號(hào),紅色為觸發(fā)啟動(dòng)標(biāo)識(shí),綠色為觸發(fā)停止標(biāo)識(shí)。圖6顯示,STA/LTA檢測(cè)算法檢測(cè)出了S1分量中利用小波分析提取的擾動(dòng)信號(hào)的開(kāi)始時(shí)間和觸發(fā)時(shí)間。如果用計(jì)算機(jī)進(jìn)行實(shí)時(shí)在線觸發(fā)檢測(cè),一旦檢測(cè)到擾動(dòng)信號(hào)就可自動(dòng)啟動(dòng)聲、光和圖形報(bào)警,并立即通過(guò)電話、微信等方式通知值班人員進(jìn)行分析處理[21]。

圖6 S1分量觸發(fā)檢測(cè)效果Fig.6 Trigger detection results of S1 component

3 結(jié) 語(yǔ)

地形變觀測(cè)過(guò)程中往往存在干擾信號(hào)和地震前兆信號(hào)疊加在固體潮上形成顫抖擾動(dòng)的情況,給地震研究帶來(lái)較大的困難。本文首先利用鉆孔應(yīng)變自洽方程檢驗(yàn)玉樹(shù)臺(tái)鉆孔應(yīng)變儀的工作狀態(tài)。結(jié)果表明,汶川8.0級(jí)地震前玉樹(shù)臺(tái)鉆孔應(yīng)變儀的觀測(cè)數(shù)據(jù)滿足自洽方程,儀器工作正常,長(zhǎng)周期信號(hào)記錄穩(wěn)定可靠。然后利用互相干分析得到玉樹(shù)臺(tái)鉆孔應(yīng)變儀S13和S24兩個(gè)合成分量中信號(hào)的相干頻率,再通過(guò)小波分析分別從S1、S2、S3和S4分量中提取擾動(dòng)信號(hào),并進(jìn)行時(shí)間-頻率-能量分析,得到地震前擾動(dòng)信號(hào)的時(shí)間-頻率-能量分布。最后利用STA/LTA檢測(cè)算法對(duì)提取出來(lái)的擾動(dòng)信號(hào)進(jìn)行觸發(fā)檢測(cè)。結(jié)果表明,選擇合適的時(shí)間窗口長(zhǎng)度和觸發(fā)閾值能夠有效檢測(cè)出擾動(dòng)信號(hào)的開(kāi)始時(shí)間和結(jié)束時(shí)間。汶川8.0級(jí)地震前2 h,玉樹(shù)臺(tái)鉆孔應(yīng)變儀4個(gè)分量都檢測(cè)到擾動(dòng)信號(hào)。

合理設(shè)置STA/LTA觸發(fā)檢測(cè)參數(shù)是檢測(cè)擾動(dòng)信號(hào)的關(guān)鍵。由于觸發(fā)檢測(cè)算法本身的局限性,會(huì)不可避免地造成誤觸發(fā)和漏觸發(fā),需要運(yùn)行維護(hù)人員掌握觸發(fā)檢測(cè)算法原理及相應(yīng)參數(shù)的意義,并結(jié)合鉆孔應(yīng)變觀測(cè)臺(tái)站的實(shí)際觀測(cè)情況不斷調(diào)整STA/LTA觸發(fā)檢測(cè)參數(shù),以達(dá)到最優(yōu)的檢測(cè)效果。

不同形變臺(tái)站觀測(cè)到的擾動(dòng)信號(hào)的分布形態(tài)、持續(xù)時(shí)間、振動(dòng)頻率和能量分布都可能存在差異,需要提取小波分析的不同層來(lái)進(jìn)行觸發(fā)檢測(cè),也可選擇多層分別進(jìn)行觸發(fā)檢測(cè)。如果多個(gè)臺(tái)站觀測(cè)到類似的擾動(dòng)信號(hào),在排除自然界和人工干擾后,可判斷為前驅(qū)波。

致謝：感謝童汪練研究員提供幫助!