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

(上海佘山基準地震臺，上海市 201602)

0 引言

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

1 小波分析的基本原理

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

運用公式：

(1)

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

(2)

(3)

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

(4)

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

2.1 小波變換對噪聲的識別和排除

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

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

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

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

2.2 小波變換對氣象變化影響的識別和排除

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

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

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

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

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

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

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

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

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

2.3 小波變換對地震(同震效應)的識別和排除

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

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

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

2.4 小波變換對鉆孔應變前兆資料趨勢與短臨異常的識別

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

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

3 結(jié)論與討論

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

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

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

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