小波分析方法在鉆孔應(yīng)變數(shù)字化資料處理中的應(yīng)用

倪友忠，吳珊珊，教聰聰，方 韜，葉 青

(上海佘山基準(zhǔn)地震臺(tái)，上海市 201602)

0 引言

利用小波變換方法可對(duì)觀測(cè)資料不同頻率段內(nèi)的信息進(jìn)行分解,從而提取到觀測(cè)資料不同頻率段的信息，為形變數(shù)字化資料的分析處理提供一種新的行之有效的分析方法。文章應(yīng)用小波變換理論，對(duì)佘山臺(tái)歷史記錄到的鉆孔應(yīng)變資料，結(jié)合宋治平等[1]、張燕等[2]提到的方法和邱澤華[3]文中通過(guò)頻率響應(yīng)和小波分解來(lái)研究氣象變化影響中采用的計(jì)算方法，求解得到小波分解細(xì)節(jié)部分d1-d4的頻率段信息以及氣象和應(yīng)變的影響系數(shù)，分析及識(shí)別觀測(cè)資料中存在的各種干擾現(xiàn)象。比如，系統(tǒng)故障、氣象、地震(同震效應(yīng))干擾等。只有通過(guò)識(shí)別及排除這些干擾因素，才可以進(jìn)一步對(duì)資料進(jìn)行計(jì)算分析，從中識(shí)別和提取地震趨勢(shì)與短臨異常。在如何識(shí)別和提取地震趨勢(shì)與短臨異常的方法上，李杰等[4-5]在小波分解方法分析形變觀測(cè)資料的正常背景變化特征和分析跨斷層形變異常等文中，都有具體提及。張燕等[6]在大姚6級(jí)雙震前地震前兆特征一文中，也具體闡述了在地震前震中附近的形變臺(tái)站，觀測(cè)資料在小波分解第三層都記錄到了相同頻段的異常信號(hào)。

1 小波分析的基本原理

小波分析是傅里葉方法的發(fā)展和延伸，分析信號(hào)頻率隨時(shí)間不斷變化。其主要作用是能在信號(hào)的低頻段具有較高的頻率分辨率，在高頻信號(hào)段具有較高的時(shí)間分辨率。能在觀測(cè)信號(hào)中，探測(cè)到正常信號(hào)中夾帶的瞬態(tài)間的反常信息，并將其反常的成分展示出來(lái)。觀測(cè)臺(tái)站記錄到的形變資料f(x)組成的序列并不是平穩(wěn)的，而是隨時(shí)間不斷變化的，其中夾雜著許多高頻干擾。Fourier變換并不能提取f(x)中的奇異性和突變的信息。小波變換理論的出現(xiàn)，對(duì)于解決和展示信號(hào)的突變和反常成分，具有很強(qiáng)的實(shí)用性?？梢岳秒x散小波變換DWT(Discrete Wavelet Transform)對(duì)觀測(cè)資料的不同頻率段的信號(hào)進(jìn)行識(shí)別并加以分析。

運(yùn)用公式：

(1)

計(jì)算中，采用a=2k。隨著k的增加，信號(hào)從最高頻向低頻分解。當(dāng)k=0時(shí)，信號(hào)為采樣頻率；k=1時(shí)將頻率二等分，依此類(lèi)推。對(duì)于數(shù)字信號(hào)f(x)，可以近似地表示為：

(2)

(3)

Mallat的重構(gòu)算法為：

(4)

2 資料處理結(jié)果分析

2.1 小波變換對(duì)噪聲的識(shí)別和排除

數(shù)字化前兆資料的特征是信息量大且干擾源多，所以在分析和研究數(shù)字化前兆資料時(shí)，必須先對(duì)干擾源進(jìn)行有效的識(shí)別并加以排除。小波變換理論具備較強(qiáng)的干擾識(shí)別和排查能力，所以近年來(lái)在數(shù)字化前兆領(lǐng)域，出現(xiàn)了不少關(guān)于小波變換方面的研究。文章在分析處理鉆孔應(yīng)變前兆信號(hào)時(shí)，初步采用Morlet小波進(jìn)行細(xì)節(jié)分析。

利用小波分析基本原理中上述公式(3)和(4)計(jì)算佘山臺(tái)鉆孔應(yīng)變面應(yīng)變Sa數(shù)據(jù)，在小波分解細(xì)節(jié)部分尺度為1-4階等高頻率段的信息。在圖1中分別列出5張小圖，分別是Sa面應(yīng)變數(shù)據(jù)、d1-d4細(xì)節(jié)1階至4階高頻信號(hào)數(shù)據(jù)。d1反映的是其中第1張小圖Sa面應(yīng)變，年變特征清晰，曲線趨勢(shì)穩(wěn)定。圖中有2處臺(tái)階突跳，均為儀器電源故障導(dǎo)致，從圖中很難識(shí)別出其他干擾信號(hào)。利用小波變換方法分解出不同頻率域的信息中，發(fā)現(xiàn)在2017年10月30日后的Sa數(shù)據(jù)的d1-d3頻率變大，特別是d1階的頻率和原來(lái)相比有很大幅度的變化，表明干擾源主要來(lái)自高頻。d1-d3主要反映的還是高頻信號(hào)段，這種影響一直持續(xù)到2018年12月12日仍未消失。在d4階時(shí)，高頻干擾信號(hào)有所減弱，固體潮信號(hào)段開(kāi)始顯示出來(lái)，表明小波分析可以清晰地分解出不同頻率域的信息，對(duì)數(shù)據(jù)資料處理有較強(qiáng)的識(shí)別能力。

圖1 鉆孔應(yīng)變面應(yīng)變分鐘變化及小波d1-d4細(xì)節(jié)Fig.1 Minute strain variation of borehole strain surface and details of wavelet d1-d4

放大2017年10月30日前后5 d的觀測(cè)數(shù)據(jù)(見(jiàn)第24頁(yè)圖2)發(fā)現(xiàn)，佘山臺(tái)四元件鉆孔應(yīng)變儀2017年10月30日電源故障導(dǎo)致斷電，觀測(cè)系統(tǒng)故障，數(shù)據(jù)造成較大幅度的漂移。儀器維修恢復(fù)后，元件3的原始觀測(cè)曲線出現(xiàn)毛刺增多，曲線變粗。所以，在利用小波變化分析時(shí)，d1-d3階頻率段即2～16 min周期段的帶通信號(hào)存在明顯的高頻異?，F(xiàn)象。之后聯(lián)系廠家了解情況，得知頻繁的斷電故障及之前使用UPS不間斷電源供電，對(duì)鉆孔元件3的傳感器電容元件造成不可修復(fù)的損傷，導(dǎo)致元件3的傳感器輸出的觀測(cè)信號(hào)曲線毛刺增多，曲線變粗，觀測(cè)精度有所下降。

2.2 小波變換對(duì)氣象變化影響的識(shí)別和排除

氣象變化對(duì)鉆孔應(yīng)變觀測(cè)的影響主要包括氣壓、氣溫及強(qiáng)降雨變化等的影響。文章利用小波變換，主要識(shí)別和排除氣壓對(duì)鉆孔應(yīng)變觀測(cè)的影響。用小波方法分解氣壓和面應(yīng)變兩條觀測(cè)曲線的觀測(cè)時(shí)間序列，可以對(duì)比分析兩條曲線d1-d8不同頻率上氣壓和觀測(cè)應(yīng)變的變化情況及其對(duì)應(yīng)關(guān)系。氣象變化在已知研究領(lǐng)域中表明，對(duì)水平剪應(yīng)變幾乎沒(méi)有影響，只對(duì)面應(yīng)變產(chǎn)生影響。所以，只提取鉆孔面應(yīng)變某時(shí)段觀測(cè)的時(shí)間序列，利用小波變換方法來(lái)加以對(duì)比，識(shí)別和排除氣壓對(duì)面應(yīng)變的干擾。

第24頁(yè)圖3的小波分解將Sa面應(yīng)變信號(hào)與氣壓信號(hào)分解至第9層，主要考察離散小波變換各層的細(xì)節(jié)系數(shù)，可以看成原始信號(hào)高頻部分不同周期段的帶通濾波信號(hào)。鉆孔應(yīng)變觀測(cè)的采樣間隔是1 min，第1層細(xì)節(jié)(d1)的信號(hào)實(shí)際上相當(dāng)于周期為2～4 min的帶通信號(hào)，第2層細(xì)節(jié)(d2)相當(dāng)于周期為4～8 min的帶通信號(hào)，依此類(lèi)推，一直到第九層(d9)。圖中，在面應(yīng)變Sa第1層至第3層細(xì)節(jié)(d1-d3)中，可以一一對(duì)應(yīng)地找到氣壓變化對(duì)面應(yīng)變的高頻干擾信號(hào)。第4層至第9層的細(xì)節(jié)(d4-d9)中，就顯得不明顯。第25頁(yè)圖4進(jìn)行小波分解后，再分別對(duì)面應(yīng)變和氣壓進(jìn)行回歸分析，得到各頻帶的氣壓影響系數(shù)。利用得到的一系列隨周期變化的氣壓影響系數(shù)(kp)，繪制曲線圖4b。

圖3a、3b為面應(yīng)變及d1-d8小波分解第1至第8階細(xì)節(jié)部分展示。d1-d4主要是高頻信號(hào)段展示，d5-d8為低頻信號(hào)段展示。圖3c、3d為氣壓及d1-d8小波分解第1至第8階細(xì)節(jié)部分展示，d1-d4主要是高頻信號(hào)段展示，d5-d8為低頻信號(hào)段展示。在d1-d4高頻信號(hào)區(qū)域，可以清晰地找到氣壓在面應(yīng)變高頻段的干擾信號(hào)，d1-d2的細(xì)節(jié)部分信號(hào)最為清晰，d5細(xì)節(jié)部分固體潮汐信號(hào)開(kāi)始顯現(xiàn)，干擾信號(hào)逐漸被排除，直至d8幾乎看不到氣壓的干擾信號(hào)。由此可知，氣壓變化對(duì)鉆孔應(yīng)變觀測(cè)的影響主要集中在高頻信號(hào)段。為了更清晰準(zhǔn)確地了解佘山臺(tái)鉆孔應(yīng)變和氣壓的對(duì)應(yīng)關(guān)系，利用小波分解及回歸分析，求得應(yīng)變氣壓的影響系數(shù)(kp)及曲線。

圖2 鉆孔應(yīng)變四元件分鐘值曲線Fig.2 Minute value curve of borehole strain four elements

圖3 面應(yīng)變、氣壓分鐘值小波變化d1-d9細(xì)節(jié)部分展示Fig.3 Details of strain and pressure minute wavelet transform d1-d9

圖4a為2016年1月1日至2017年9月30日面應(yīng)變Sa及應(yīng)變氣壓影響系數(shù)變化圖，可見(jiàn)佘山臺(tái)面應(yīng)變的變化趨勢(shì)大致與氣壓趨勢(shì)同步，有明顯年變規(guī)律；圖4b的影響系數(shù)顯示，佘山地區(qū)氣壓對(duì)四元件鉆孔面應(yīng)變Sa的影響。由此可見(jiàn)，應(yīng)變氣壓影響系數(shù)周期在1 min，數(shù)值接近0.3 ns/hPa的時(shí)候，開(kāi)始逐漸上升；周期10 min時(shí)，數(shù)值達(dá)到峰值1.5 ns/hPa；周期大于10 min后，應(yīng)變氣壓影響系數(shù)反轉(zhuǎn)開(kāi)始下降；周期達(dá)到50 min時(shí)，應(yīng)變氣壓影響系數(shù)趨向于0。該應(yīng)變氣壓影響系數(shù)曲線同樣清晰表明，在佘山臺(tái)鉆孔應(yīng)變觀測(cè)中，氣壓只在短周期段對(duì)鉆孔面應(yīng)變有影響，最大影響頻率為10 min，對(duì)長(zhǎng)周期數(shù)據(jù)無(wú)影響。

圖4 面應(yīng)變、氣壓分鐘值及應(yīng)變氣壓影響系數(shù)圖Fig.4 Variation of strain, pressure minute value and strain pressure influence coefficient

圖5為氣壓差分?jǐn)?shù)據(jù)(氣壓變化影響kp×DP曲線)與佘山臺(tái)面應(yīng)變差分?jǐn)?shù)據(jù)進(jìn)行對(duì)比。從圖5清晰找出氣壓對(duì)鉆孔應(yīng)變的影響，表明氣壓對(duì)佘山臺(tái)面應(yīng)變存在影響，影響周期主要集中在0～10 min的頻率域范圍內(nèi)，小波分解細(xì)節(jié)d1-d3周期段。

圖5 面應(yīng)變分鐘值差分、氣壓變化影響(kp×DP)曲線Fig.5 The strain minute difference and the effect of pressure change (kp×DP) curve

2.3 小波變換對(duì)地震(同震效應(yīng))的識(shí)別和排除

鉆孔應(yīng)變觀測(cè)可以清晰地記錄到遠(yuǎn)震、大震或近震、中小地震的同震效應(yīng)(即地震面波)。但這種同震效應(yīng)幅度大，夾雜在觀測(cè)曲線中。對(duì)于前兆分析來(lái)說(shuō)，也屬于一種干擾，這種干擾在趨勢(shì)異常分析時(shí)，應(yīng)予以識(shí)別和排除。小波分解方法是排除此種干擾最直接有效的方法之一(見(jiàn)圖6)。

圖6 小波變換識(shí)別面應(yīng)變同震效應(yīng)d1-d4細(xì)節(jié)部分Fig.6 Wavelet transform identifies details of surface strain coseismic effect d1-d4

2018年9月28日18:02:44印度尼西亞發(fā)生7.4級(jí)地震，上海市佘山臺(tái)鉆孔應(yīng)變四個(gè)元件測(cè)項(xiàng)，都清晰地記錄到了同震效應(yīng)。圖6展示了鉆孔應(yīng)變儀元件1測(cè)項(xiàng)記錄到該地震的同震效應(yīng),利用小波變換方法提取了第1階至第4階細(xì)節(jié)的高頻信號(hào)。說(shuō)明小波方法對(duì)于同震效應(yīng)的識(shí)別、提取效果較好。

2.4 小波變換對(duì)鉆孔應(yīng)變前兆資料趨勢(shì)與短臨異常的識(shí)別

選取2013至2018年鉆孔應(yīng)變面應(yīng)變Sa的日均值進(jìn)行小波分析。第26頁(yè)圖7a、7b、7c分別是面應(yīng)變Sa未剔除干擾的日均值變化圖、Sa剔除干擾后的日均值變化圖、Sa小波分解d3階細(xì)節(jié)部分高頻信號(hào)。Sa未剔除干擾日均值變化圖中的干擾包含地震同震效應(yīng)、儀器故障缺數(shù)、氣象干擾等。在日均值圖中表現(xiàn)為脈沖突跳，變化幅度比較大。利用matlab差分函數(shù)進(jìn)行差分計(jì)算后，剔除部分干擾，展現(xiàn)出圖7b面應(yīng)變Sa的日均值長(zhǎng)期變化趨勢(shì)。從Sa的日均值長(zhǎng)趨勢(shì)變化圖中，了解到佘山地區(qū)的面應(yīng)變長(zhǎng)期受壓向下。圖7c是Sa的日均值進(jìn)行小波分解后第3階細(xì)節(jié)部分的高頻信息圖。從圖中可以看出，在2014年3月3日琉球M6.7地震前，Sa的日均值有明顯鼓包現(xiàn)象，小波分解第三階細(xì)節(jié)部分的高頻信號(hào)幅度超出平時(shí)的2倍以上，連續(xù)且密集，持續(xù)了近半年左右，震后該高頻脈沖信號(hào)幅度明顯降低；2015年11月14日東海海域M7.2地震發(fā)生前，Sa的日均值鼓包不太明顯，但Sa的日均值小波分解d3細(xì)節(jié)部分出現(xiàn)了持續(xù)4個(gè)月左右的密集，幅度超過(guò)平靜日2倍以上的高頻脈沖信號(hào)。2016年1至6月，Sa的日均值小波分解d3細(xì)節(jié)部分也曾陸續(xù)出現(xiàn)幅度較大的高頻脈沖信號(hào)，但由于沒(méi)有出現(xiàn)Sa的日均值鼓包，且高頻信號(hào)不連續(xù)，所以該處的d3細(xì)節(jié)部分反映出來(lái)的高頻信號(hào)不是地震前兆異常。

圖7 面應(yīng)變?nèi)站底兓瘓D、小波分解第3階細(xì)節(jié)部分Fig.7 Daily mean strain variation diagram, third-order detail of wavelet decomposition

3 結(jié)論與討論

綜上所述，得出如下結(jié)論。

(1) 中小波變換在數(shù)字化前兆資料分析中可以展現(xiàn)d1-d4階細(xì)節(jié)部分的高頻信號(hào)，能清晰地識(shí)別出地震(同震效應(yīng))、氣象干擾(氣壓等)、儀器故障等干擾，具有較強(qiáng)的干擾識(shí)別與排除功能。

(2) 小波變換方法對(duì)鉆孔應(yīng)變長(zhǎng)趨勢(shì)的識(shí)別、提取和短臨趨勢(shì)異常也有較強(qiáng)的識(shí)別能力。圖7在琉球群島M6.7和東海海域M7.2地震前，在d3階細(xì)節(jié)部分都明顯地提取到了幅度較大，周期大概持續(xù)5～6個(gè)月的高頻信號(hào)；且Sa的日均值有鼓包現(xiàn)象，地震發(fā)生在鼓包回落恢復(fù)的后期，震后該現(xiàn)象明顯消失，故判斷很可能為地震前兆異常。

(3) 2018年7月后，Sa的日均值第3階細(xì)節(jié)部分開(kāi)始出現(xiàn)明顯的高頻脈沖，Sa的日均值也有明顯鼓包，這一異?，F(xiàn)象是否為海域地震前兆異常，值得進(jìn)一步關(guān)注研究。