榮成地震臺形變觀測資料小波分析

徐芳芳 崔居全 荊 強 安巴雅爾 孫宗強

1)中國山東 264300 榮成地震臺

2)中國山東 264000 威海市地震局

3)中國云南 657800 水富縣防震減災局

0 引言

小波分析(wavelet analysis)是目前信號處理分析的重要方法，在時頻兩域具有表征信號局部特征的能力，是一種窗口大小固定不變而形狀可改變，時間窗和頻率窗均可改變的時頻局部化分析方法，即在低頻部分具有較低的時間分辨率和較高的頻率分辨率，在高頻部分具有較高的時間分辨率和較低的頻率分辨率，適于分析非平穩信號和提取信號的局部特征。宋治平等（2003）、李杰等（2005）、張燕等（2004，2009）、李艷等（2011）、張軍等（2011）對數字化地殼形變前兆資料進行一系列研究，并取得一定成果，可見小波變換理論對地震前兆觀測資料分析具有廣闊的應用前景。

采用Matlab軟件，運用小波分析方法，對榮成地震臺（以下簡稱榮成臺）CZB-Ⅱ豎直擺鉆孔傾斜儀觀測資料進行處理，研究形變數據的波動、人為標定、同震響應、固體潮等。根據以往研究成果，采用db4作為母小波，對形變觀測資料進行4—7層小波分解，近似部分清晰可見數據趨勢變化，細節部分能清晰顯示包括1/4日波、半日波、日波和半月波等在內的固體潮汐波，并能清晰可見數據中固體潮汐波的周期特征及隨時間變化情況。

1 小波分析理論

1.1 小波變換理論及小波基函數選取

實際工作中數據信號是離散的，需要對信號進行離散小波變換DWT（Discrete Wavelet Transform）。對于離散序列信號f(x)，在小波函數ψ(t)∈L2(R)中，尺度因子(伸縮因子a和平移因子b）（a，b∈R）也需要離散化，則應用小波變換作為不同頻率的信息識別基礎（宋志平等，2003），即

在資料處理中，采用a=2k。隨著k的增加，信號從最高頻向低頻分解。當k=0時，信號為采樣頻率；k=1時將頻率二等分，依此類推。

對于數字信號f(x)，可以近似表示為

小波分解的快速算法稱為馬拉( Mallat) 算法。Mallat 算法（薛志芳等，2012）在小波分析中的地位相當于快速傅里葉變換算法在經典傅里葉分析中的地位。Mallat 算法為正交展開系數 與 分別滿足以下關系。

而Mallat重構算法為

式中，hn與gn分別是尺度函數和小波函數濾波器，且=h-n，=g-n。

小波基不是唯一的，其選取沒有明確標準，選擇合適的小波基進行信號處理需要考慮其正則性和消失矩、所處理信號與小波基的相似性等因素，正交小波變換對此優勢較大。Daubechies小波是離散正交小波，可表示為dbN，其中N是小波的階。dbN不具對稱性（即非線性相位），沒有明確的解析表達式，小波函數ψ與尺度函數φ的有效支撐長度為2N－1，ψ的消失矩為N。根據以往研究結果（宋治平，2003；李杰等，2005），發現db4小波就提取形變資料異常效果較好，故本研究繼續采用db4小波。

1.2 頻率分析

圖1 小波分解近似及細節部分與頻率關系Fig.1 Relationship among approximation part, detail part and frequency of wavelet decomposition

小波變換系數與實際為窗口小波變換。當j固定時，隨著k的變化，與均占滿時間域未重疊；隨著j的變化， 占滿頻率域未重疊。

設信號f(x)在時間域[0，T]采樣，共采N個點，采樣間隔為Δt=T/N。{fi}(i=0，1，…，N－1)的頻帶為[0，Ω]，求出近似部分的系數{}(k=0，1，…，N－1)，其頻帶為 [0，Ω]。頻帶為 [0，Ω]，頻帶為[0，Ω/2]，的頻率為 [Ω/2，Ω]，以此類推，與的頻帶分別為 [0，Ω/2j]與 [Ω/2j，Ω/2j-1](圖 1)。

通過小波變換方法可對不同頻率范圍內的信息進行識別與分離。根據公式(4)及小波分解的近似部分(x)與細節部分(x)與頻率的關系（圖1），對觀測資料進行近似部分(低頻)與細節部分(高頻)信息分離。根據觀測量的物理含義，在分析干擾因素的基礎上，從近似部分中確定趨勢變化信息，從細節部分識別短期異常。

2 小波分析

榮成臺位于山東半島最東端，地處海西頭—俚島斷裂西南15 km，乳山斷裂東80 km，構造涉及膠遼斷塊南部和北黃海凹陷，距離東部海域6 km。臺站安裝3套形變儀器，TJ-Ⅱ體應變儀于2006年9月安裝運行，2012年6月因探頭故障停測；CZB-Ⅱ豎直擺鉆孔傾斜儀位于榮成臺院內西側，2007年7月正式運行；RZB-Ⅱ分量應變儀位于院內西南角，于2014年12月24日安裝運行，因運行時間較短，主要采用豎直擺傾斜儀數據。本研究選擇儀器運行以來正常、固體潮清晰、數據連續可靠的形變測項的月小時或分鐘值及日分鐘值數據進行小波分析，對個別缺數進行線性插值處理，保證資料的連續可靠。

2.1 高頻與低頻信息識別

圖2為傾斜儀北南分量2014年5月—10月整點值曲線及小波變換的4階a4近似部分和d1—d4細節部分。從細節部分看出，1階信號主要為最高頻率信息，包含突跳干擾、噪聲信號等，2階信號可以看出1/4日波、半日波、日波信息；3階、4階信號可以看出半日波、日波、半月波信息，從4階近似部分可以看出數據的變化趨勢。由圖2可見，2014年6月13日因數采面板電壓不穩定，數據出現階變；7月24—25日降大暴雨，降雨量達243 mm，因降雨地面上覆土層荷載增大，數據出現階變。由此可知，細節部分主要反映固體潮汐半日波、日波、半月波的高頻變化信息，近似部分信息反映整體趨勢變化過程，因此可以將小波分析用于地震前兆觀測資料年變趨勢對比分析，尋找長期異常變化?？梢?，應用小波分析可對不同頻率范圍內的信息進行有效識別。

圖2 2014年5—10月數據小波變換近似部分和細節部分Fig.2 Inclinometer NS direction data wavelet approximation part and detail part from May to Oct.，2014

2.2 不同頻率信息分離

借助 Matlab（2013）小波分析軟件，以傾斜儀北南分量2014年4月的分鐘值為例，對形變觀測資料進行近似部分（低頻）與細節部分（高頻）信息的分離。圖3分別給出傾斜儀分鐘值原始信號s、小波變換低頻（a1，a2，…，a5）與高頻（d1，d2，…，d6）曲線。隨著分離層數的增加，高頻信息被剔除，使固體潮信號隨時間尺度的變化趨勢更加明顯、直觀，清晰顯示數據固體潮信息特征。圖3(a)給出不同頻率范圍的近似部分（低頻）信息，反映數據整體變化趨勢。圖3(b)給出不同頻率范圍細節（高頻）部分，反映數據短期變化趨勢。從圖3清晰可見，4月的數據有6次同震效應，包括大震遠震和近中、小地震，經過6階分解，可消除地震同震效應產生的影響。由圖3(a)可見，數據經過6階分解后，干擾消除，顯示出原始變化趨勢，清晰可見4月26—28日因降雨影響數據出現的趨勢轉折，并能準確看出受降雨影響的時間段。

選取圖3中兩種代表性數據進行處理：①同震效應：2014年4月2日智利8.1級地震同震效應分鐘值曲線[圖3(a)，圖4]；②儀器標定：2014年4月30日標定儀器日分鐘值曲線[圖3(b)，圖5]。數字化儀器精度較高，能記錄到同震效應（地震面波）。同震效應是影響數據分析的干擾，需要進行識別和排除。由圖4近似部分可見，數據經過分解，同震波被消除，顯示出原始曲線變化趨勢，從細節部分可以提取同震波。可見，小波變換對同震效應的識別、提取及消除均具有良好效果。經過分解，將2014年4月30日人為標定脈沖干擾去除，反映了原始數據的變化趨勢，細節部分也提取到標定干擾信息，并可見具體時間。由圖3—圖5可知，干擾、脈沖標定和地震同震記錄經6層分離基本消失，6階曲線較平滑，分離效果較好。

圖3 2014年4月傾斜儀北南向數據小波分解（a）近似部分；（b）細節部分Fig.3 Inclinometer NS direction data wavelet decomposition in April, 2014

圖4 2014年4月2日同震效應的小波分解Fig.4 Wavelet decomposition for coseismic effects on April 2, 2014

圖5 2014年4月30日儀器標定的小波分解Fig.5 Wavelet decomposition for instrument calibration on April 30, 2014

2.3 地震前兆異常

以2013年5月18日06∶02∶23黃海海域(37.7°N， 124.7°E)M5.1地震（震中距210 km）為例，分析地震前兆異常。采用多種常規方法，對2013年5月傾斜儀北南向和東西向數據進行分析，發現從2013年5月15日開始出現趨勢異常。對2013年5月整點值數據進行小波變換，得到6階近似部分與細節部分（圖6）。由圖6可見，自2013年5月15日開始，數據出現趨勢異常變化，5月18日發生黃海海域M5.1地震?？梢?，應用小波分析方法，可對地震趨勢異常與短期異常進行有效識別。

圖6 2013年5月1日—31日數據小波分解Fig.6 Wavelet decomposition for the data from May 1 to 31, 2013

3 結束語

通過小波分析方法在榮成臺形變觀測數據處理的應用，驗證該信號處理方法對于數據波動、儀器標定、同震效應及不同頻率信息的識別、分離具有較好效果，一些地震短臨變化也能清晰突顯，使各種長短期震前異常識別變得容易，有利于利用地震前兆數據進行地震預測。威海地區可以嘗試采用小波分析方法進行地震數據處理，并可以推廣到其他地震前兆觀測項及地震臺。

李杰，劉希強，李紅，等.利用小波變換方法分析形變觀測資料的正常背景變化特征[J].地震學報，2005，27（1）：33-41.

李艷，高振強，馮建琴，等.臨汾地震臺體應變觀測資料映震能力分析[J].地震地磁觀測與研究，2011，32（4）：99-104.

薛志芳，于春頌，李非，等.河北昌黎地震臺大地電場數字化觀測資料小波分析[J].山西地震，2012，152（4）：22-40.

張燕，吳云，劉永啟，等.小波分析在地殼形變資料處理中的應用[J].地震學報，2004，26（Supp）：103-109.

張燕，吳云.2008年汶川地震前的形變異常及機理解釋[J].武漢大學學報·信息科學版，2010，35（1）：25-30.

張軍，陳宇衛，劉澤民，等.小波分析在前兆數據處理中的應用[J].地震地磁觀測與研究，2011，26（2）：99-104.

宋治平，武安緒，王梅，等.小波分析方法在形變數字化資料處理中的應用[J].大地測量與地球動力學，2003，23（4）：21-27.