烏加河地震臺(tái)DSQ水管傾斜儀觀測(cè)的小波及閾值去噪分析

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

（烏加河地震臺(tái)，內(nèi)蒙古自治區(qū) 巴彥淖爾市 015323）

0 引言

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

1 臺(tái)站概況

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

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

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

2.1 小波變換基本原理

小波變換是一種信號(hào)的時(shí)頻分析方法，具有多尺度（分辨率）分析的特點(diǎn)，能有效區(qū)分平穩(wěn)信號(hào)中疊加的非平穩(wěn)部分，對(duì)數(shù)字化前兆資料的干擾識(shí)別與去除以及對(duì)不同頻率的信息識(shí)別功能較強(qiáng)[4]。

離散信號(hào)的小波變換可表示為：

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

采用Mallat算法，將數(shù)字信號(hào)f（x），分解成不同頻率成分[9]：

式中，j為分解階數(shù)，Aj f（x）是信號(hào)f（x）的頻率不超過(guò)2-j的部分，稱近似部分，Dj f（x）分是信號(hào)f（x）的頻率介于2-j與2-j+1之間的部分，稱細(xì)節(jié)部分，Aj與Dj分別為二進(jìn)制離散近似與二進(jìn)制離散細(xì)節(jié)，Aj與Dj滿足如下關(guān)系：

式中，H為尺度函數(shù)對(duì)應(yīng)的低通濾波器，G為小波函數(shù)對(duì)應(yīng)的高通濾波器，J為最佳尺度。

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

由圖2可知，與小波分解階數(shù)j對(duì)應(yīng)的近似部分Aj f（x）和細(xì)節(jié)部分Dj f（x）的頻帶分別為［0，fN/2j］、［fN/2j，fN/2j-1］。每一階的近似部分Aj f（x）相當(dāng)于對(duì)原始數(shù)據(jù)進(jìn)行低通濾波，過(guò)濾了D1f（x）、D2f（x）…Dj f（x）等細(xì)節(jié)部分，只保留原始序列中頻帶為［0，fN/2j］的信號(hào)主體部分。

圖2 應(yīng)用Mallat塔式算法對(duì)信號(hào)進(jìn)行小波多尺度分解示意圖Fig.2 Schematic diagram of wavelet multi-scale decomposition of signal using Mallat tower algorithm

2.2 小波基的選擇

小波變換與Fourier變換的區(qū)別之一是小波變換的小波基函數(shù)具備多樣性，選擇適合的小波基函數(shù)是利用小波分析識(shí)別并去除數(shù)字前兆異常信號(hào)的前提條件[9]。小波基函數(shù)選擇的一般原則為：小波基函數(shù)具備緊支性、對(duì)稱性、正則性（光滑性）和恰當(dāng)?shù)南Ь仉A數(shù)等[10-11]。根據(jù)尹繼堯、王嘉琦和孫伶俐等[12-14]的研究，應(yīng)用Daubechies（dbN）小波分析方法可對(duì)形變觀測(cè)數(shù)據(jù)中的噪聲和干擾進(jìn)行有效地識(shí)別。選取dbN小波為小波基函數(shù)可兼顧正交小波的緊支性和正則性[15]。

為選取合適的小波基，本文采用控制變量法：將分解層數(shù)均設(shè)置為8（對(duì)應(yīng)的細(xì)節(jié)周期范圍為2～512 min，可分解出固體潮信息），閾值均使用固定閾值，閾值去噪方法均為軟閾值去噪，通過(guò)計(jì)算同一含噪信號(hào)的不同小波尺度均方根誤差，確定最佳分解和重構(gòu)尺度。

選取含噪信號(hào)為2020年11月25日水管儀NS向的數(shù)據(jù)，表1為計(jì)算的均方根誤差結(jié)果。

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

表1中RMSE為均方根誤差，該指標(biāo)可體現(xiàn)原始信號(hào)與去噪后信號(hào)之間的整體偏差信息，數(shù)值越小，表明去噪效果越好[16]。其計(jì)算公式為：

式f（i）為原始含噪信號(hào)，f′（i）為去噪后信號(hào)，n為數(shù)據(jù)長(zhǎng)度。

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

為驗(yàn)證db5為合適的小波基，選取2018年6月18日水管儀NS向數(shù)據(jù)進(jìn)行小波變換，圖3為應(yīng)用db5的小波重建信號(hào)及其誤差。

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

由圖3可以看出，小波重構(gòu)后的數(shù)據(jù)與原始數(shù)據(jù)有一致的曲線形態(tài)，且小波重建誤差非常微小。

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

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

綜上，應(yīng)用db5對(duì)烏加河臺(tái)水管儀數(shù)據(jù)進(jìn)行小波變換是合適的。

3 資料選取

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

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

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

3.1 地球物理事件干擾

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

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

由圖4可知，該次地球物理事件的干擾主要體現(xiàn)在細(xì)節(jié)d1～d3的相對(duì)高頻部分，對(duì)應(yīng)周期為2～16 min；在細(xì)節(jié)d4～d8的相對(duì)低頻部分，地球物理事件干擾基本被過(guò)濾出去，主要為水管儀記錄到的包括固體潮在內(nèi)的低頻信息。

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

3.2 氣壓干擾

氣壓大幅度突變會(huì)引起形變儀器觀測(cè)曲線的畸變，畸變程度和氣壓變化幅度成正相關(guān)。2019年6月9日15時(shí)至21時(shí)，烏加河地區(qū)氣壓發(fā)生較大幅度變化，短時(shí)內(nèi)變化幅度達(dá)3.2 hPa。從圖5可以看出：氣壓變化對(duì)水管儀數(shù)據(jù)造成干擾，干擾形態(tài)主要為噪聲大，和氣壓曲線同步性較好；氣壓干擾主要體現(xiàn)在細(xì)節(jié)d1～d3，對(duì)應(yīng)周期為2～16 min；在細(xì)節(jié)d4～d8的相對(duì)低頻部分，氣壓干擾基本被過(guò)濾出去，主要為水管儀記錄到的包括固體潮在內(nèi)的低頻信息。

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

3.3 大風(fēng)干擾

烏加河臺(tái)位于內(nèi)蒙古西部地區(qū)，每年春秋兩季大風(fēng)天氣較多，是形變數(shù)字化觀測(cè)的主要干擾之一，尤其對(duì)水管儀等長(zhǎng)基線儀器的干擾明顯。2019年10月3日烏加河地區(qū)大幅降溫，伴有6～7級(jí)大風(fēng)天氣，導(dǎo)致水管儀觀測(cè)曲線出現(xiàn)明顯擾動(dòng)，表現(xiàn)為觀測(cè)曲線上疊加有毛刺干擾信號(hào)，數(shù)據(jù)噪聲變大（圖6）。

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

從圖6可知，大風(fēng)干擾主要體現(xiàn)在細(xì)節(jié)d1～d5，對(duì) 應(yīng) 周 期 為：2～64 min；細(xì) 節(jié)d6～d8為水管儀記錄到的包括固體潮在內(nèi)的低頻背景信息。

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

4 閾值去噪

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

閾值降噪應(yīng)用最廣泛的是Donoho提出的軟、硬閾值降噪法。硬閾值去噪是將待處理的信號(hào)中絕對(duì)值小于閾值的小波系數(shù)置零，軟閾值去噪則是在硬閾值的基礎(chǔ)上，將待處理信號(hào)絕對(duì)值大等于閾值的小波系數(shù)再減去閾值，使數(shù)據(jù)處理邊界連續(xù)光滑[18]。因此，本文采用軟閾值降噪法處理數(shù)據(jù)。

本次研究針對(duì)每一個(gè)小波細(xì)節(jié)選取一個(gè)固定閾值，以此避免有用信號(hào)的小波系數(shù)損失過(guò)重，而噪聲的小波系數(shù)沒(méi)有得到有效壓制[19]。

固定閾值Thr的計(jì)算公式[11]為：

式中，σ為標(biāo)準(zhǔn)偏差，n為信號(hào)長(zhǎng)度。由上述公式，可計(jì)算出各種干擾類(lèi)型的閾值（表4）。

表4 DSQ水管儀各類(lèi)干擾的閾值

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

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

5 結(jié)論

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