小波分析法在海拉爾地震臺前兆數(shù)據(jù)處理中的應(yīng)用研究

曲浩鑫，王 怡，堵偉鵬

(內(nèi)蒙古自治區(qū)地震局海拉爾地震臺，內(nèi)蒙古 海拉爾 021000)

0 引言

小波分析是一種包含尺度伸縮和時間平移的雙參數(shù)函數(shù)分析方法，能將一個信號分解成對空間和尺度的獨立結(jié)果, 同時又保留原信號所包含的信息[1]?；谛〔ㄗ儞Q能分離出不同頻帶的信號，目前，小波分析法廣泛應(yīng)用于天然地震震相識別，重力、大面積形變測量和定點形變測量等地震異常分析中[2-7]。宋治平等研究表明，小波分析法對形變數(shù)字化資料中的干擾識別與消除、不同頻率潮汐波的分離、趨勢異常與短期異常的提取等方面都具有較強的功能[8]；張軍等應(yīng)用小波分析法對地下流體資料進行研究，得出較好抑制流體數(shù)據(jù)中隨機噪聲的方法，并通過重構(gòu)小波系數(shù)達到消除降雨干擾的目的[9]；劉建明等研究得出小波分析可作為提取定點形變異常的一種有效方法，通過將定點形變觀測數(shù)據(jù)分解成不同頻段，對提取異常信號有較好的效果，提高了識別地震信息的能力[6]；倪友忠等通過小波變換法對面應(yīng)變?nèi)站颠M行分析，發(fā)現(xiàn)在海域大震前，小波變換第三階細節(jié)部分在震前出現(xiàn)持續(xù)數(shù)月以上的高頻脈沖，且面應(yīng)變有鼓包現(xiàn)象[10]。

因此，在小波變換理論基礎(chǔ)上，運用小波分析法對海拉爾地震臺(以下簡稱海拉爾臺)前兆數(shù)據(jù)進行處理，總結(jié)小波分析在前兆資料中的應(yīng)用經(jīng)驗，以期為識別、提取及去除前兆異常信息提供新思路。

1 小波分析理論基礎(chǔ)

此時稱Ψ(t)為一個基本母小波，將母函數(shù)Ψ(t)經(jīng)伸縮和平移后得到如下公式：

式中：a為伸縮因子；b為平移因子。稱其為一個小波序列。在實際運用中，連續(xù)小波必須加以離散化。任意函數(shù)f(t)的離散小波變換可表示為：

為區(qū)分不同分量信息，引入構(gòu)造正交小波的方法和小波變換的快速算法—Mallat。應(yīng)用Mallat算法不必知道具體的小波函數(shù)，同時，也使小波變換在計算上變得可行。Mallat一維分解算法可表示為：

式中：Cj+1f與Dj+1f是在分辨率為j時的近似部分和細節(jié)部分。對于Mallat一維重構(gòu)算法可表示為：

Cj-1[n]=∑nh(n-2k)Cj[k]+∑ng(n-2k)Dj[k],

式中：Cj[k]為低頻系數(shù)；Dj[k]為高頻系數(shù)；h，g分別為低通和高通濾波器。

2 小波基函數(shù)的理論選取

在Mallat算法中小波基函數(shù)的選擇不是唯一的，所有滿足小波條件的函數(shù)都可作為小波函數(shù)，不同的小波基函數(shù)使信號的分解和重構(gòu)有不同的結(jié)果。因此,選取合適的小波基函數(shù)就顯得較重要。目前，常用的小波基函數(shù)有：Haar小波、Daubechies(dbN)小波系、Biorthogonal(biorNr.Nd)小波系、Coiflet(coifN)小波系、SymletsA(symN)小波系、Morlet(morl)小波、Meyer小波等。常用小波基函數(shù)的主要特點如表1所示，選取的小波基函數(shù)既能突出異常信號又能最大限度地保留前兆資料的固體潮信息。

表1 常用小波基的主要特點Table 1 Main characteristics of commonly used wavelet bases

根據(jù)小波函數(shù)的定義和常用小波基的特點，結(jié)合前兆信號資料的特征及研究目的，綜合考慮小波基函數(shù)的對稱性、正則性、緊支撐性及消失矩的階數(shù)[6,12-16]。理論上，Daubechies(dbN)小波系中的db4小波基函數(shù)能滿足前兆資料處理的要求。db4小波的尺度函數(shù)和小波函數(shù)如圖1所示，較好地兼顧了緊支撐性和正則性，同時，消失矩的階數(shù)也較合適，能在保留固體潮信息的基礎(chǔ)上較好地反應(yīng)異常信號。

為檢驗db4小波在實際前兆數(shù)字化資料處理中的效果，選取2020年1月24日海拉爾臺垂直擺NS分量數(shù)據(jù)進行小波變換，圖2為應(yīng)用db4小波基函數(shù)進行離散小波變換的結(jié)果。可以看出，2020年1月24日15點09分塔吉克斯坦5.6級地震的信號在高頻分

圖1 db4小波尺度函數(shù)及小波函數(shù)Fig.1 db4 wavelet scaling function and wavelet function

解中變得突出，同時，應(yīng)用近似5階小波重構(gòu)消除了該地震波對垂直擺正常波形的影響，垂直擺的固體潮信息被較好地保留。由此可見，應(yīng)用db4小波基函數(shù)進行離散小波變換對海拉爾臺前兆異常信號的識別、提取及消除是適用的。

圖2 垂直擺NS分量曲線Fig.2 NS component curve of vertical pendulum

3 小波分析剔除干擾

海拉爾臺目前共有5套數(shù)字化前兆觀測儀器，前兆資料信息量大且不同儀器會受到不同程度干擾，在原始曲線中不易識別和區(qū)分。小波變換理論具有較強的干擾識別和排查能力，因此，選取海拉爾臺受氣壓干擾和場地環(huán)境影響的典型事件，采用小波分析法進行識別和消除干擾。

3.1 氣象干擾

第22頁圖3a為2020年8月2日至12日伸縮儀NS分量原始數(shù)據(jù)曲線，第22頁圖4a為2020年7月1日至3日水管儀NS分量原始數(shù)據(jù)曲線，從圖中難以識別出具體的干擾時段。應(yīng)用小波變換法對圖3a、圖4a原始數(shù)據(jù)曲線進行細節(jié)8階分解，分解細節(jié)部分如圖3d、圖4d所示?？梢钥闯?，d1-d2未見明顯異常，從d3開始曲線出現(xiàn)明顯異常(圖中方框部分)，上下波動幅度超過2倍中誤差。通過與同時段氣壓變化曲線進行對比發(fā)現(xiàn)，d3-d6明顯的上下階躍部分與氣壓異常波動時段一一對應(yīng)，因此，可以判定伸縮儀和水管儀的異常變化由氣壓突變引起。

在小波分解d1-d8中，d1-d4主要是高頻信號段部分，d5-d8為低頻信號段部分，氣壓影響在d3-d6部分清晰顯示，d8頻帶部分明顯過濾了異常變化，顯示的是固體潮信息。綜上所述，海拉爾臺形變儀受氣壓干擾在小波分解d3-d6部分較明顯。

圖3c為通過小波重構(gòu)剔除受氣壓影響后的曲線，與原始曲線圖3a相比，在保留原始信息的情況下，較好地去除了氣壓影響造成的固體潮畸變，證明小波分析法對排除氣壓干擾具有較好效果。

圖3 伸縮儀2020年8月2日至12日NS分量曲線Fig.3 NS component curve of extensometer from August 2 to 12, 2020

圖4 水管儀2020年7月1日至3日NS分量曲線Fig.4 NS component curve of water tube instrument from July 1 to 3, 2020

3.2 場地環(huán)境干擾

由第23頁圖5a看出,從2018年10月25日10時30分左右至27日8時30分左右，伸縮儀受干擾出現(xiàn)固體潮畸變，經(jīng)核查為體應(yīng)變鉆新井工程影響所致，屬于典型的場地環(huán)境干擾。通過小波細節(jié)8階分解顯示有3個異常時間段(見圖5c)，分別為25日10時至16時、26日 13時至22時、27日7時至8時，三個時間段與三次鉆井時間相吻合。從圖5a原始曲線中較難判斷出具體受干擾的時段，小波分解d3-d6部分壓制了其余信息，將場地環(huán)境干擾凸顯出來，其中細節(jié)5階最為明顯。將場地環(huán)境干擾信息提取之后對原始曲線進行干擾剔除，剔除結(jié)果如圖5b所示。與原始曲線相比，曲線光滑度高，高頻干擾信號基本被抑制。

圖5 伸縮儀2018年10月25日至27日NS分量曲線Fig.5 NS component curve of extensometer from October 25 to 27, 2018

4 小波分析地震與爆破

4.1 小波分析近震與遠震

海拉爾臺不同前兆儀器識別地震的能力不同，多種手段互補觀測可以清晰地記錄到遠震、近震以及地方震的同震效應(yīng)[8，10]。對于前兆分析來講，這種同震效應(yīng)幅度大，夾雜在觀測曲線中，過大的振動幅度會壓制前兆儀器正常的固體潮信息，也屬于一種干擾。因此，采用小波分析法對同震效應(yīng)造成的突跳進行識別和處理是必要的。

圖6a為海拉爾臺2017年5月12日重力潮汐觀測值原始曲線，可以看出，呼倫貝爾市鄂倫春旗3.7級地震明顯，中美洲沿岸遠海6.2級地震難以辨別。從圖6c中看出，d1-d8中有兩處明顯異常，是由19時11分的中美洲沿岸遠海6.2級遠震和22時19分的呼倫貝爾市鄂倫春旗3.7級近震造成的，在原始曲線上難以辨別的中美洲沿岸遠海地震經(jīng)過小波變換后清晰可見。第24頁圖7a為海拉爾臺2019年9月30日垂直擺EW分量原始曲線，其中，智利中部沿岸近海地震為遠震，呼倫貝爾市新巴爾虎右旗地震為近震，近震在原始曲線中不明顯，經(jīng)過小波變換后可輕易識別。

由圖6c、圖7c看出,近震和遠震在小波變換細節(jié)分解中具有完全不同的特征，遠震在d3-d5階中異常明顯；近震在d1-d3、d6-d8階中較明顯，d1-d3階異常幅度逐漸減弱，d6-d8階異常范圍的寬度逐漸增加。圖6b、圖7b為去除干擾后的重力潮汐觀測曲線，比原始數(shù)據(jù)曲線清晰，日變形態(tài)正常顯現(xiàn)。通過以上分析可以得出，小波分析對同震效應(yīng)的識別、提取及干擾消除均具有較好的效果。

4.2 小波分析爆破

第24頁圖8a為海拉爾臺垂直擺NS分量2019年6月10日原始數(shù)據(jù)曲線，可以看出，11時左右北南分量數(shù)據(jù)形態(tài)表現(xiàn)為上下階躍，難以判斷此階躍是由地震還是爆破引起。因此，對曲線進行小波細節(jié)8階分解，結(jié)果如圖8b所示。對比圖8b、8c發(fā)現(xiàn)，爆破與地震兩種信號在小波細節(jié)分解上有明顯區(qū)別。爆破信號在小波分解d1-d8階中均清晰可見，異常的寬度和幅度變化不大；地震信號在d1-d8階中的某幾部分異常明顯，同時，地震類型的不同，異常的寬度和幅度變化呈現(xiàn)不同的態(tài)勢。這是由于地震信號比爆破信號復(fù)雜，信號的頻帶范圍更寬所致。第25頁圖9為海拉爾臺2020年4月7日垂直擺EW分量記錄到的爆破，通過做小波分解發(fā)現(xiàn)與圖8具有相同的規(guī)律，即在小波分解d1-d8階中爆破信號的寬度和幅度變化不大。由此可見，利用小波細節(jié)分解辨別原始信號異常是由地震還是爆破引起具有較好的效果。

圖6 2017年5月12日重力儀曲線Fig.6 Gravimeter curve of May 12, 2017

圖7 2019年9月30日垂直擺曲線Fig.7 Vertical pendulum curve on September 30, 2019

5 結(jié)論與討論

通過以上分析，得出如下認識：

(1) 小波變換可將不同頻帶的信號較好地分離，實現(xiàn)對各頻段小波系數(shù)的逐一分析，從而達到提取前兆異常、去除特定頻段干擾的目的。

(2) 通過對不同小波基函數(shù)特征的歸納分析,并對前兆資料進行處理，驗證了db4小波是海拉爾臺前兆數(shù)字化資料處理較適合的小波基函數(shù)，能滿足在保留固體潮信息的基礎(chǔ)上檢測前兆異常。

圖8 垂直擺2019年6月10日NS分量曲線Fig.8 NS component curve of vertical pendulum on June 10, 2019

圖9 2020年4月7日垂直擺EW分量曲線Fig.9 EW component curve of vertical pendulum on April 7, 2020

(3) 由海拉爾臺形變儀器小波分析結(jié)果可知，小波分析法對氣壓、場地環(huán)境干擾的識別和剔除具有較好效果。氣壓、場地環(huán)境干擾均在小波細節(jié)分解d3～d6部分異常顯著，其中d5階最明顯，d8階表征的基本是固體潮信息。應(yīng)用小波變換重構(gòu)信號基本去掉了氣壓、場地環(huán)境對原始信號的影響，去噪效果顯著，信號曲線光滑度高，基本抑制了噪聲信號。

(4) 小波變換法在海拉爾臺同震效應(yīng)的識別、提取和去除具有較好效果。通過對比遠震和近震小波變換結(jié)果得出，遠震在小波細節(jié)分解d3-d5部分異常顯著，近震在小波細節(jié)分解d1-d3、d6-d8部分明顯，且d1-d3部分異常幅度逐漸減弱，d6-d8部分異常范圍的寬度逐漸增加。

(5) 利用小波細節(jié)分解法可對原始異常相似的地震和爆破信號進行辨別。地震信號只在小波細節(jié)分解d1-d8階中的某幾部分異常明顯，同時，地震類型的不同，異常的寬度和幅度變化呈現(xiàn)不同的態(tài)勢；爆破信號在小波細節(jié)分解d1-d8階中均清晰可見，異常的寬度和幅度變化不大。