烏加河地震臺DSQ水管傾斜儀觀測的小波及閾值去噪分析

胡 瑋，馮雪東，石 偉，郭 偉，劉 偉

（烏加河地震臺，內蒙古自治區 巴彥淖爾市 015323）

0 引言

數字形變儀器測量會受到多種因素的干擾，包括自然、人為、地球物理事件、儀器故障等因素[1-2]。由于小波分析具有能夠根據分析對象自動調整有關參數的“自適應性”和能夠根據觀測對象自動“調焦”的特性，而廣泛應用于干擾識別及噪聲去除等研究領域[3]。宋治平和倪友忠等[4-5]的研究表明，小波分析可應用在地震震相識別和前兆異常分析等領域；張軍等[6]應用小波分析方法，在流體學科中較好地抑制了隨機噪聲；劉建明等[7]通過運用小波分析將形變觀測數據分解為不同的頻段，提高了識震能力。去除形變觀測中的干擾信號，對提高觀測質量、識別前兆異常有一定的作用。目前對定點形變儀器的小波分析及閾值去噪的研究較少，本文應用db小波分析方法，對烏加河地震臺DSQ水管傾斜儀所受的干擾進行了研究，總結分析其干擾規律和特征，并通過軟閾值去噪方法剔除干擾，取得較好效果。

1 臺站概況

烏加河地震臺（以下簡稱烏加河臺）是國家級綜合觀測臺站，地處狼山前東西向大斷裂帶和河套斷陷盆地北緣陰山緯向構造帶的中西段與狼山弧型構造帶的復合部位[8]。觀測學科有：形變、測震、地磁、地電和重力等。烏加河臺形變山洞于1972年人工開鑿，山洞洞室進深253 m，臺基屬花崗巖，完整性好，平均覆蓋厚度189 m，觀測條件較好（圖1）。烏加河臺DSQ水管儀于2001年進行了數字化改造。

圖1 烏加河臺山洞洞室平面示意圖Fig.1 The cavity plan of the digital deformation instrument at Wujiahe Seismic Station

2 小波分析理論基礎

2.1 小波變換基本原理

小波變換是一種信號的時頻分析方法，具有多尺度（分辨率）分析的特點，能有效區分平穩信號中疊加的非平穩部分，對數字化前兆資料的干擾識別與去除以及對不同頻率的信息識別功能較強[4]。

離散信號的小波變換可表示為：

式中，f（x）為離散的信號序列，Ψa，b（x）為小波函數，a為尺度（伸縮）因子，b為平移因子。

采用Mallat算法，將數字信號f（x），分解成不同頻率成分[9]：

式中，j為分解階數，Aj f（x）是信號f（x）的頻率不超過2-j的部分，稱近似部分，Dj f（x）分是信號f（x）的頻率介于2-j與2-j+1之間的部分，稱細節部分，Aj與Dj分別為二進制離散近似與二進制離散細節，Aj與Dj滿足如下關系：

式中，H為尺度函數對應的低通濾波器，G為小波函數對應的高通濾波器，J為最佳尺度。

設信號f（x）的采樣率為Fs，根據采樣定理，f（x）的最大頻率為奈奎斯特頻率fNyquist=Fs/2，記為fN，f（x）頻帶可表示為［0，fN］，則Mallat塔式算法可圖示為：

由圖2可知，與小波分解階數j對應的近似部分Aj f（x）和細節部分Dj f（x）的頻帶分別為［0，fN/2j］、［fN/2j，fN/2j-1］。每一階的近似部分Aj f（x）相當于對原始數據進行低通濾波，過濾了D1f（x）、D2f（x）…Dj f（x）等細節部分，只保留原始序列中頻帶為［0，fN/2j］的信號主體部分。

圖2 應用Mallat塔式算法對信號進行小波多尺度分解示意圖Fig.2 Schematic diagram of wavelet multi-scale decomposition of signal using Mallat tower algorithm

2.2 小波基的選擇

小波變換與Fourier變換的區別之一是小波變換的小波基函數具備多樣性，選擇適合的小波基函數是利用小波分析識別并去除數字前兆異常信號的前提條件[9]。小波基函數選擇的一般原則為：小波基函數具備緊支性、對稱性、正則性（光滑性）和恰當的消失矩階數等[10-11]。根據尹繼堯、王嘉琦和孫伶俐等[12-14]的研究，應用Daubechies（dbN）小波分析方法可對形變觀測數據中的噪聲和干擾進行有效地識別。選取dbN小波為小波基函數可兼顧正交小波的緊支性和正則性[15]。

為選取合適的小波基，本文采用控制變量法：將分解層數均設置為8（對應的細節周期范圍為2～512 min，可分解出固體潮信息），閾值均使用固定閾值，閾值去噪方法均為軟閾值去噪，通過計算同一含噪信號的不同小波尺度均方根誤差，確定最佳分解和重構尺度。

選取含噪信號為2020年11月25日水管儀NS向的數據，表1為計算的均方根誤差結果。

表1 不同小波尺度下的均方根誤差

表1中RMSE為均方根誤差，該指標可體現原始信號與去噪后信號之間的整體偏差信息，數值越小，表明去噪效果越好[16]。其計算公式為：

式f（i）為原始含噪信號，f′（i）為去噪后信號，n為數據長度。

由表1可以得出，隨著小波尺度的增大，RMSE先減小，后增大。當小波尺度選取為5時，RMSE最小。

為驗證db5為合適的小波基，選取2018年6月18日水管儀NS向數據進行小波變換，圖3為應用db5的小波重建信號及其誤差。

圖3 2018年6月18日水管儀NS向小波重建信號及其誤差Fig.3 The wavelet reconstruction signal and its error of NS direction of water-tube tiltmeter on June 18，2018

由圖3可以看出，小波重構后的數據與原始數據有一致的曲線形態，且小波重建誤差非常微小。

為驗證選擇db5為小波基的合理性，表2計算了2018年6月18日水管儀NS向數據1-8階小波分解下的RMSE。從表中可以得出，當階數是5時，RMSE最小。

表2 不同小波尺度下的均方根誤差

綜上，應用db5對烏加河臺水管儀數據進行小波變換是合適的。

3 資料選取

烏加河臺DSQ水管儀記錄精度較高，儀器穩定，固體潮形態清晰，且能夠清晰記錄到各種干擾。通過查看水管儀觀測日志和數據跟蹤分析平臺可知，水管儀的干擾因素主要有：自然環境干擾、地球物理事件干擾和人為干擾等。由于水管儀北南向和東西向觀測干擾類型及其特征相似，故本文只選取了北南向數據進行了分析。本文統計了2018—2020年水管儀北南向受到的各類干擾及其主要特征（表3）。

表3 水管儀北南向干擾特征（2018—2020年）

在各類干擾中，人為干擾由于受觀測人員進洞的影響，其干擾幅度是最大的，干擾形態也是最為明顯的，易于識別。人為干擾時段內的數據是錯誤數據，根據形變學科要求，預處理時應將其做缺數處理，予以直接剔除，不必使用小波分析方法。2018年8月27日至9月1日的降雨事件，在烏加河地區屬罕見天氣，累計降雨量103.8 mm，造成水管儀觀測數據趨勢性變化，2019年2月才基本恢復正常。一般而言，前兆觀測資料所受干擾的周期遠小于固體潮，即有用信號為低頻，而干擾為高頻，通過小波分析，可將其高頻干擾部分逐步分解出來并通過閾值去噪進行剔除，因此降雨干擾用小波分析方法效果不佳。

3.1 地球物理事件干擾

烏加河臺水管儀對全球6級以上地震均能清晰記錄其同震響應，干擾形態主要為突跳。圖3為水管儀記錄到2020年6月26日5時5分新疆和田地區（35.73°N，82.33°E）6.4級地震，震源深度10 km，震中距2316 km，最大變化幅度為6.02 ms。

從原始數據中很難看出受地球物理事件干擾的頻段信息，應用小波分析方法，對水管儀北南向原始數據進行8層分解，分解出了不同頻段的細節信息，d1～d8為該地震的小波分解1～8層細節，周期范圍為2～512 min，對應Mallat塔式算法中的D1f（x）—D8f（x）（圖2），f（x）為水管儀原始數據序列，采樣率為1 min-1，其頻帶為0～0.5 min-1。

由圖4可知，該次地球物理事件的干擾主要體現在細節d1～d3的相對高頻部分，對應周期為2～16 min；在細節d4～d8的相對低頻部分，地球物理事件干擾基本被過濾出去，主要為水管儀記錄到的包括固體潮在內的低頻信息。

圖4 2020年6月25日地球物理事件干擾及其小波分析Fig.4 Geophysical event interference and its wavelet analysis on June 25，2020

3.2 氣壓干擾

氣壓大幅度突變會引起形變儀器觀測曲線的畸變，畸變程度和氣壓變化幅度成正相關。2019年6月9日15時至21時，烏加河地區氣壓發生較大幅度變化，短時內變化幅度達3.2 hPa。從圖5可以看出：氣壓變化對水管儀數據造成干擾，干擾形態主要為噪聲大，和氣壓曲線同步性較好；氣壓干擾主要體現在細節d1～d3，對應周期為2～16 min；在細節d4～d8的相對低頻部分，氣壓干擾基本被過濾出去，主要為水管儀記錄到的包括固體潮在內的低頻信息。

圖5 2019年6月9日至10日氣壓干擾及其小波分析Fig.5 Air Pressure Interference and its wavelet analysis from June 9 to 10，2019

3.3 大風干擾

烏加河臺位于內蒙古西部地區，每年春秋兩季大風天氣較多，是形變數字化觀測的主要干擾之一，尤其對水管儀等長基線儀器的干擾明顯。2019年10月3日烏加河地區大幅降溫，伴有6～7級大風天氣，導致水管儀觀測曲線出現明顯擾動，表現為觀測曲線上疊加有毛刺干擾信號，數據噪聲變大（圖6）。

圖6 2019年10月2日至4日大風干擾及其小波分析Fig.6 Wind Interference and its wavelet analysis from October 2 to 4，2019

從圖6可知，大風干擾主要體現在細節d1～d5，對 應 周 期 為：2～64 min；細 節d6～d8為水管儀記錄到的包括固體潮在內的低頻背景信息。

大風干擾與氣壓干擾在干擾頻次、幅度和形態上均表現出相似性（表1），但利用小波分析可看出二者的不同：氣壓擾動持續時間較短，擾動頻段為1/16～1/2 min-1，擾動主要體現在1/8～1/2 min-1范圍內；大風干擾持續時間更長，分布在更寬的頻帶（1/64～1/2 min-1）上，能量主要集中在1/16～1/2 min-1，在1/64～1/16 min-1頻段上仍存在其相對低頻部分，而在這一頻段上，氣壓干擾對水管儀幾乎沒有產生影響。

4 閾值去噪

小波去噪的方法主要包括：屏蔽去噪法、模極大值檢測法以及閾值去噪法等，其中閾值去噪法為最常見的去噪手段[17]。由前面的分析可知，噪聲主要為相對高頻的細節部分，信號主體在低頻的近似部分。閾值去噪先計算出有用信號與噪聲的小波系數之間的閾值，再將小于閾值的小波系數置零，實現去噪目的[16]。

閾值降噪應用最廣泛的是Donoho提出的軟、硬閾值降噪法。硬閾值去噪是將待處理的信號中絕對值小于閾值的小波系數置零，軟閾值去噪則是在硬閾值的基礎上，將待處理信號絕對值大等于閾值的小波系數再減去閾值，使數據處理邊界連續光滑[18]。因此，本文采用軟閾值降噪法處理數據。

本次研究針對每一個小波細節選取一個固定閾值，以此避免有用信號的小波系數損失過重，而噪聲的小波系數沒有得到有效壓制[19]。

固定閾值Thr的計算公式[11]為：

式中，σ為標準偏差，n為信號長度。由上述公式，可計算出各種干擾類型的閾值（表4）。

表4 DSQ水管儀各類干擾的閾值

根據表4中的閾值，對水管儀地球物理事件、氣壓干擾和大風干擾數據進行去噪處理，在去除高頻干擾的同時，較好地保留了原始數據中固體等有用信號，明顯提高了數據的信噪比，去噪效果良好（圖7）。在接下來的工作中，作者擬應用小波閾值去噪方法提取前震異常信號，開展深入研究。

圖7 地球物理事件、氣壓干擾和大風干擾原始信號及小波閾值去噪信號對比Fig.7 Comparison of original signal and wavelet threshold denoising signal of Geophysical Event Interference,Air Pressure Interference and Wind Interference

5 結論

通過以上分析，可以得到以下結論：（1）烏加河臺水管傾斜觀測主要受人為干擾、降雨干擾、地球物理事件干擾、氣壓干擾和大風干擾的影響；（2）利用小波分解不同頻段的信號可知，地球物理事件干擾主要分布在細節d1～d3，對應的周期為2～16 min；氣壓干擾主要分布在細節d1～d3，對應的周期為2～16 min；大風干擾主要分布在細節d1～d5，對應的周期為2～64 min；（3）應用db5對烏加河臺水管儀數據進行小波分析和軟閾值去噪是適用的，基本去除了各種類型對原始數據的噪聲干擾，去噪效果明顯，去噪后的數據曲線光滑度高，同時保留了原始固體潮信息；（4）運用小波分析和閾值去噪方法可以較好去除干擾信號，能清楚的觀測出固體潮等有用信息，對定點形變觀測異常識別有一定的作用，可供形變臺站參考。