高 鐵，徐英奎，顏 萍，吳燕岡

（沈陽地質(zhì)調(diào)查中心方法室，沈陽 110034）

小波閾值去噪方法在重力資料處理中的應(yīng)用

高 鐵，徐英奎*，顏 萍，吳燕岡

（沈陽地質(zhì)調(diào)查中心方法室，沈陽 110034）

簡要介紹了小波變換的基本理論與小波閾值去噪實(shí)現(xiàn)過程。通過單球體與雙球體模型試驗(yàn)，探討了閾值去噪時(shí)小波基以及分解層數(shù)的選擇問題；最后利用小波閾值去噪方法，對(duì)某實(shí)測(cè)布格重力資料進(jìn)行處理，有效濾除了隨機(jī)噪聲。去噪后計(jì)算的水平方向?qū)?shù)明顯比去噪前直接計(jì)算的水平方向?qū)?shù)效果好。

小波閾值去噪；小波基；分解層數(shù)；重磁資料

0 引言

在重磁勘探中，由于受到各種因素的影響，使得重磁資料中常含有較強(qiáng)的噪聲，這些伴隨在有用信號(hào)中的噪聲是影響重磁資料后續(xù)處理精度的重要因素（例如在各種反演方法中常需要求取異常的水平和垂向高階導(dǎo)數(shù)，但高階導(dǎo)數(shù)方法對(duì)高頻干擾噪聲具有較強(qiáng)的放大作用，會(huì)大大降低反演的精度）。因此在對(duì)重磁資料進(jìn)行各種特殊處理之前，在不降低信號(hào)分辨率的基礎(chǔ)上，對(duì)其進(jìn)行去噪處理是十分必要的。

重磁資料去噪傳統(tǒng)的方法主要是基于傅里葉變換（例如匹配濾波、維納濾波等），但是傅里葉變換是把信號(hào)從時(shí)域完全變換到頻域進(jìn)行分析，計(jì)算過程中丟失了時(shí)域的全部信息，在時(shí)域方面分辨性很差，而在頻域方面不能反映信號(hào)的細(xì)節(jié)［1－4］。

小波變換是20世紀(jì)80年代發(fā)展起來的一種新的方法技術(shù)，由于其具有時(shí)頻分析、多分辨率和去相關(guān)性等特點(diǎn)，在信號(hào)去噪方面得到了廣泛地應(yīng)用。目前，已有許多學(xué)者做了大量的試驗(yàn)，研究利用小波變換對(duì)地震信號(hào)、測(cè)井曲線、探地雷達(dá)信號(hào)等方面的去噪［2－10］，但是在重磁資料的應(yīng)用，研究還比較少［11］。本研究首先介紹了小波變換的基本理論以及小波閾值去噪的原理與過程，然后探討了在重磁資料去噪處理中，閾值函數(shù)、小波基以及空間尺度的選取問題，最后將小波閾值去噪方法應(yīng)用到某實(shí)測(cè)布格重力異常中。

1 方法原理

1．1 小波變換

設(shè)ψ（t）∈L2（R），若其傅里葉變換滿足條件為一個(gè)基本小波或母小波。將母函數(shù)ψ（t）經(jīng)伸縮和平移后得ψa，b（t），稱其為一個(gè)小波序列。

其中：a為伸縮因子；b為平移因子。

對(duì)于任意函數(shù)f（t）∈L2（R）的連續(xù)小波變換為

小波變換的一個(gè)重要特征是有一個(gè)靈活可變的時(shí)間－頻率窗，它在高的“中心頻率”時(shí)自動(dòng)變窄，而在低的“中心頻率”時(shí)自動(dòng)變寬，這樣它對(duì)高頻信號(hào)有較高的分辨率，而對(duì)低頻信號(hào)又能給出完整的信息［］。

1．2 小波閾值去噪

閾值去噪方法的思想就是對(duì)小波分解后的各層系數(shù)中模大于和小于某閾值的系數(shù)分別處理，然后對(duì)處理完的小波系數(shù)再進(jìn)行反變換，重構(gòu)出經(jīng)過去噪后的信號(hào)［12］。具體步驟如下［5］：

1）選擇一個(gè)小波基并確定分解層數(shù)，對(duì)原始含噪數(shù)據(jù)進(jìn)行小波變換，得到含噪小波系數(shù)。

2）選取一種去噪準(zhǔn)則，確定噪聲的閾值，給出小波系數(shù)估值。

3）由小波系數(shù)估值進(jìn)行小波反變換，重構(gòu)信號(hào)得到信號(hào)估值。

常用的閾值函數(shù)有硬閾值函數(shù)和軟閾值函數(shù)：

1）硬閾值函數(shù)。它可以很好地保留信號(hào)邊緣等局部特征，其表達(dá)式為

2）軟閾值函數(shù)。它的處理結(jié)果會(huì)相對(duì)圓滑，但是會(huì)造成邊緣模糊等失真現(xiàn)象，其表達(dá)式為

3）半軟閾值函數(shù)。它可兼顧軟閾值和半軟閾值方法的優(yōu)點(diǎn)，其表達(dá)式為

2 模型試驗(yàn)

重力異常是由地下淺部與深部所有地質(zhì)體的綜合反映，所以等值線會(huì)比較圓滑，不會(huì)出現(xiàn)跳變點(diǎn)，因此選用軟閾值比硬閾值更適合重力資料的去噪。

小波去噪時(shí)需要考慮的另外兩個(gè)參數(shù)是小波基和分解層數(shù)，選擇不同的小波解決同一個(gè)問題可能會(huì)產(chǎn)生不同的結(jié)果，因此去噪處理時(shí)必須考慮選擇最優(yōu)小波基，而隨著分解層數(shù)的增加，噪聲的幅值越來越小，有用信號(hào)的幅值逐漸增大，因此要選擇合適的分解層數(shù)以達(dá)到既保持原有信號(hào)的有用部分，又能濾去干擾部分的效果。

小波基和分解層數(shù)的選擇對(duì)去噪的結(jié)果有很大的影響，前人已經(jīng)做過許多理論分析和實(shí)驗(yàn)：通過對(duì)地震信號(hào)進(jìn)行小波閾值去噪后的信號(hào)與原始信號(hào)的誤差大小進(jìn)行分析，得出sym族小波優(yōu)于db族小波，而sym N小波族中，N取4～6的中等數(shù)，得到的誤差最小，而地震信號(hào)中的噪聲只存在于前三尺度的分解結(jié)果中，故采用3次分解去噪即可消除地震中的噪聲［6］；db小波比較適合探地雷達(dá)數(shù)據(jù)處理［8］，分解層數(shù)取“3”既能保證去噪的效果，又能盡量減少處理的工作量［9］。

本研究通過對(duì)疊加隨機(jī)噪聲的單球體模型和雙球體模型重力異常，選用不同的小波基以及分解層數(shù)進(jìn)行去噪試驗(yàn)，利用去噪后的信號(hào)與原始信號(hào)的標(biāo)準(zhǔn)偏差作為誤差分析的標(biāo)準(zhǔn)，來探討重磁資料小波閾值去噪處理中小波基和分解層數(shù)的選擇問題。

2．1 單球體模型

選取單個(gè)球體為模型，其參數(shù)為：球體中心坐標(biāo)為（0，0），埋深50 m，半徑100 m，剩余密度為1．0 g／cm3。其重力異常等值線如圖1（a）所示，最大值為3．5 g．u．。疊加幅度為（－0．1，0．1）的隨機(jī)噪音后，重力異常等值線如圖1（b）所示。表1為不同小波基以及分解層數(shù)的誤差結(jié)果。當(dāng)選取的小波基為db6，分解層數(shù)為4層時(shí)的小波閾值去噪后的重力異常等值線見圖1（c）。

2．2 雙球體模型

在單球體模型的基礎(chǔ)上再疊加一個(gè)球體，其參數(shù)為：球體中心坐標(biāo)為（50，40），埋深20 m，半徑15 m。此球體埋深相對(duì)淺，規(guī)模也小，產(chǎn)生的重力異常相當(dāng)于淺部的高頻有效信號(hào)。其重力異常等值線如圖2（a）所示，最大值為4．42 g．u．，同樣疊加幅度為（－0．1，0．1）的隨機(jī)噪音后，重力異常等值線如圖2（b）所示。表2為其不同小波基以及分解層數(shù)的誤差結(jié)果。當(dāng)選取的小波基為db6，分解層數(shù)為3層時(shí)的小波閾值去噪后的重力異常等值線見圖2（c）。

2．3 參數(shù)選擇

由表1和表2可知，Haar方波的誤差很大，而當(dāng)N和分解層數(shù)任意時(shí)，coif族小波的誤差變化范圍比sym族小波及db族小波大，所以sym族小波與db族小波優(yōu)于haar小波及coif族小波。sym N和db N小波族中，N不能取“1”，其他取值對(duì)結(jié)果影響很小，最好取3～6；分解層數(shù)取3層最優(yōu)。

圖1 單球體模型小波閾值去噪Fig．1 The application of wavelet threshold de－noising method to single sphere model

圖2 雙球體模型小波閾值去噪Fig．2 The application of wavelet threshold de－noising method to double sphere model

表1 不同小波基以及分解層數(shù)的誤差結(jié)果比較（單球體模型）Tab．1 The error by different wavelet basis function and decomposition order（single model）

3 實(shí)測(cè)數(shù)據(jù)處理

圖3是某區(qū)實(shí)測(cè)的布格重力異常等值線圖（單位毫伽：10－5m／s2），從圖3可以看出，等值線總體呈東北方向展布，西北高東南低。由于測(cè)量和各項(xiàng)改正過程中各種誤差的存在，導(dǎo)致等值線雜亂，甚至在某些地方出現(xiàn)明顯的沿南北測(cè)線方向的條帶狀虛假異常，如果在處理之前不對(duì)這些干擾去除，將嚴(yán)重影響后續(xù)各種反演與解釋的精度。采用二維小波閾值去噪對(duì)實(shí)測(cè)資料進(jìn)行處理，采用的小波基為sym7，尺度為3，得到結(jié)果如圖4所示。從圖4可以看出，去噪后的等值線完好地保持了原有的整體形態(tài)，并且在沒有產(chǎn)生虛假異常的基礎(chǔ)上，等值線變得相對(duì)圓滑，去除了各種隨機(jī)干擾以及南北向條帶狀虛假異常。

表2 不同小波基以及分解層數(shù)的誤差結(jié)果比較（雙球體模型）Tab．2 The error by different wavelet basis function and decomposition order（double model）

水平方向?qū)?shù)可以突出某個(gè)特定方向的異常，是重磁數(shù)據(jù)處理中的一種經(jīng)典方法，在此引入用以說明去噪的重要性。圖5是直接用實(shí)測(cè)的布格重力異常利用頻率域方法得到的水平方向?qū)?shù)，由于噪音的存在，求導(dǎo)過程中將高頻成分放大，導(dǎo)致結(jié)果非常雜亂，肉眼很難辨別出規(guī)律。圖6是先做小波閾值去噪處理然后再求導(dǎo)得到的水平方向?qū)?shù)，四個(gè)方向的導(dǎo)數(shù)等值線都比較光滑，很容易辨別出東北向的異常最強(qiáng)。

4 結(jié)論

圖3 實(shí)測(cè)布格重力異常等值線圖Fig．3 The practical gravity data contours

圖5 原始布格重力異常水平方向?qū)?shù)Fig．5 Directional derivatives calculated from the original gravity data

由于小波變換具有時(shí)頻分析、多分辨率和去相關(guān)性等特點(diǎn)，在信號(hào)去噪方面得到了廣泛地應(yīng)用。通過建立簡單的球體模型，用小波閾值方法對(duì)其疊加了隨機(jī)噪聲后的數(shù)據(jù)進(jìn)行去噪處理，探討了在重磁資料去噪處理中，軟硬閾值函數(shù)、小波基以及分解層數(shù)的選取問題。結(jié)論如下：

1）軟閾值函數(shù)比較適合于重磁數(shù)據(jù)去噪。

2）通過誤差分析比較，各類小波都可以很好地濾除噪聲，但db N族小波和sym N族小波優(yōu)于其他小波基，N不能取1，其他取值對(duì)結(jié)果影響很小，最好取3～6。

3）分解層數(shù)取3層去噪效果最佳。

最后將小波閾值去噪方法應(yīng)用到某實(shí)測(cè)布格重力異常中，有效地濾除了隨機(jī)噪聲，去噪后計(jì)算的水平方向?qū)?shù)明顯比去噪前直接計(jì)算的水平方向?qū)?shù)效果好。經(jīng)理論和實(shí)測(cè)數(shù)據(jù)的處理結(jié)果表明，小波閾值去噪方法能夠有效地去除重磁資料中的隨機(jī)噪聲，是一種切實(shí)可行的重磁信號(hào)去噪新方法。

圖4 小波閾值去噪后的重力等值線圖Fig．4 The contours after de－noising

圖6 去噪后的布格重力異常水平方向?qū)?shù)Fig．6 Directional derivatives calculated from the denoised gravity data

［1］ 方興．小波分析及其去噪性能研究 ［J］．浙江師范大學(xué)學(xué)報(bào)：自然科學(xué)版，2003，26（4）：345－349．

［2］ 中國生，徐國元，江文武．基于小波變換的爆破地震信號(hào)去噪的應(yīng)用［J］．中南大學(xué)學(xué)報(bào)：自然科學(xué)版，2006，37（1）：155－159．

［3］ 李英，張淑貞，許康生．小波降噪方法在地震信號(hào)處理中的應(yīng)用［J］．西北地震學(xué)報(bào)，2006，28（2）：159－162．

［4］ 張旭東，詹毅，馬永琴．小波變換與傅立葉變換在地震資料去噪中的對(duì)比［J］．內(nèi)蒙古石油化工，2007，7：29－44．

［5］ 姜弢，劉慶普，胡留軍．地震信號(hào)去噪的小波方法研究［J］．哈爾濱工程大學(xué)學(xué)報(bào)，2002，23（4）：86－90．

［6］ 柳建新，韓世禮，馬捷．小波分析在地震資料去噪中的應(yīng)用［J］．地球物理學(xué)進(jìn)展，2006，21（2）：541－545．

［7］ 鄒長春，楊欣德，潘令枝．一種基于小波變換的測(cè)井曲線去噪新方法［J］．物探與化探，1999，23（6）：462－466．

［8］ 鄒海林，寧書年，林捷．小波理論在探地雷達(dá)信號(hào)處理中的應(yīng)用［J］．地球物理學(xué)進(jìn)展，2004，19（2）：268－275．

［9］ 劉秀娟，鄧世坤，徐保林．探地雷達(dá)信號(hào)小波變換去噪［J］．工程勘察，2006，10：66－71．

［10］何蘭芳，王緒本，何展翔，等．MT時(shí)間序列的小波去噪分析［J］．地震地質(zhì)，2001，23（2）：222－226．

［11］謝仕民，王子靜，鐘鈞宇．小波降噪技術(shù)在地磁數(shù)據(jù)庫預(yù)處理中的應(yīng)用研究［J］．戰(zhàn)術(shù)導(dǎo)彈控制技術(shù)，2008，30（2）：15－22．

［12］周偉．MATLAB小波分析高級(jí)技術(shù)［M］．西安：西安電子科技大學(xué)出版社，2006．

Application of wavelet threshold de－noising method in gravity data processing

GAO Tie，XU Ying－kui，YAN Ping，WU Yan-gang
（Shenyang Center China Geological Survey，Shenyang 110034，China）

In this paper，the basic theory of wavelet transform and the process of wavelet threshold denoising are introduced．Then，we discussed the criterion for wavelet bases and scales to be selected by establishing single sphere and double spheres models and doing various experiments on the gravity of the models．Finally，the wavelet threshold de－noising method is applied to practical gravity data，and the result is completely denoised．Directional derivatives calculated from the denoised data are more effective than that calculated from the original data

wavelet threshold de－noising；wavelet bases；wavelet scales；gravity and magnetic data

P 631．1

A

10．3969／j．issn．1001-1749．2014．05．09

1001-1749（2014）05-0566-05

2014-03-03 改回日期：2014-07-23

沈陽地質(zhì)礦產(chǎn)研究所自設(shè)項(xiàng)目（2012007）

高鐵（1983-），男，博士，主要研究方向?yàn)榈厍蛭锢砭C合方法，E-mail：39085496＠qq．com。

*通訊作者：徐英奎（1963-），男，工程師，研究方向?yàn)闃?gòu)造地質(zhì)，E-mail：997182861＠qq．com。