航空重力測量數(shù)據(jù)的小波濾波處理

王靜波, 熊盛青，羅鋒，王冠鑫

(1.北方工業(yè)大學(xué) 理學(xué)院，北京 100144;2.中國自然資源航空物探遙感中心，北京 100083)

0 引言

在飛行狀態(tài)下進(jìn)行的航空重力測量，必然受到飛機(jī)發(fā)動機(jī)的震動、氣流的顛簸和氣流引起的飛機(jī)高度變化等因素產(chǎn)生擾動加速度的影響[1-9]。擾動加速度的量級很大(可達(dá)106mGal)，并以高頻成分為主，而重力異常信號幅值較小(約102mGal左右)，完全被高頻噪聲干擾所淹沒，通過低通濾波壓制高頻噪聲，從而獲得相對低頻部分有意義的重力異常信號[1-9]。然而，噪聲干擾可能涵蓋整個頻帶，低頻信息既包含重力異常信號，也包含噪聲信號。因此，如何既保留低頻有效信息，又有效去除低頻噪聲，成為高精度航空重力測量需要解決的技術(shù)難題。

小波分析是應(yīng)用數(shù)學(xué)和工程學(xué)科中一個迅速發(fā)展的新領(lǐng)域，經(jīng)過近30年的探索研究，重要的數(shù)學(xué)形式化體系已經(jīng)建立，理論基礎(chǔ)更加扎實(shí)[10-11]。小波變換具有良好的時頻局部化特性，在信號處理中可用于弱信號的提取分離，而對信號進(jìn)行多尺度分解，則可以抑制噪聲(濾波)。國內(nèi)眾多學(xué)者開展了小波濾波技術(shù)研究，并在不同領(lǐng)域得到廣泛應(yīng)用。在大地測量中，柳林濤等構(gòu)造三類連續(xù)小波，用于航空重力測量數(shù)據(jù)濾波處理[12]，孫中苗等初步探討了小波閾值濾波法在航空重力測量數(shù)據(jù)處理中的適用性[13]；在重力勘探領(lǐng)域，羅鋒等選擇Daubechies N小波系和小波閾值濾波法，對引進(jìn)的俄制GT-1A航空重力勘查系統(tǒng)獲取的測量數(shù)據(jù)進(jìn)行濾波試驗(yàn)研究[11]；在基于飛艇的地空電磁探測中，李肅義等采用sym8小波，對地空電磁測量數(shù)據(jù)進(jìn)行了綜合濾波方法技術(shù)研究[14]；在變形監(jiān)測領(lǐng)域，章浙濤等采用小波包分析和基于頻率順序的信息分段的多閾值準(zhǔn)則，對變形監(jiān)測數(shù)據(jù)進(jìn)行了濾波處理研究[15]。盡管上述研究取得了一定的進(jìn)展，但對于航空重力測量數(shù)據(jù)處理來說，與傳統(tǒng)濾波技術(shù)一樣，還需研究、設(shè)計(jì)相應(yīng)的小波低通濾波器，根據(jù)需要方便使用。此外，在小波選取、小波系數(shù)閾值處理方案等方面，也還需進(jìn)行深入、系統(tǒng)的研究。研究分析不同類型小波的不同特性，從理論和應(yīng)用兩方面綜合考慮，選取合適的小波用于航空重力數(shù)據(jù)處理；在通用閾值處理方案的基礎(chǔ)上，應(yīng)根據(jù)航空重力異常信號的頻率特性和噪聲分布特點(diǎn)，研究適合于航空重力測量數(shù)據(jù)處理的小波系數(shù)閾值處理方案[16]。

小波包分析不僅像小波分析一樣對信號的低頻部分進(jìn)行分解，也對信號的高頻部分進(jìn)行分解[15,17-23]。相對而言，小波包分析對信號在頻率域可提供更復(fù)雜和精細(xì)的分解，且通過對各分解層信號的優(yōu)化組合，可獲得接近濾波器濾波頻率范圍的信號?？紤]到航空重力測量信號頻率的分布特征，本文擬采用小波包閾值濾波法研究、設(shè)計(jì)所需的小波低通濾波器，并利用它對引進(jìn)俄制GT-1A航空重力勘查系統(tǒng)獲得的測量數(shù)據(jù)進(jìn)行濾波試驗(yàn)研究。

1 小波包閾值濾波

小波包分析是小波分析的一種推廣，它為信號分析處理提供了更精細(xì)的信號分解途徑。作為小波包分析在信號分析中的應(yīng)用之一，小波包閾值濾波過程包含4個基本步驟[17,22]：① 選取一種小波，確定小波包分解層級N，對信號進(jìn)行小波包分解；② 優(yōu)化設(shè)計(jì)小波包分解樹；③ 對小波包分解系數(shù)進(jìn)行閾值處理；④ 基于小波包分解原系數(shù)和修改的系數(shù)，進(jìn)行小波包重構(gòu)計(jì)算，恢復(fù)信號。

1.1 信號的小波包分解與重構(gòu)

假設(shè)離散信號序列為s(k)(k=1,2,3,…,N)，則可對信號S={s(k)}進(jìn)行離散小波分解，其二進(jìn)3層小波分解樹如圖1所示，小波分解單元如圖2所示，j為小波分解層次，Aj、Dj為第j層的小波分解系數(shù)，它們分別代表分解信號的低頻分量和高頻分量，亦稱之為低頻系數(shù)和高頻系數(shù)[17]。

圖1 小波包分解樹(3層)Fig.1 Wavelet decomposition tree at level 3

小波包分析不僅對信號的低頻分量進(jìn)行分解，也對高頻分量進(jìn)行分解，信號S={s(k)}的3層小波包分解樹如圖3所示，小波包分解單元如圖4所示，j為小波包分解層次，n為第j層的小波包節(jié)點(diǎn)序號(n可能值：0,1,…,2j-1)，dj,n為小波包節(jié)點(diǎn)(j,n)對應(yīng)的小波包分解系數(shù)[17]。

選用正交或雙正交小波，則小波包系數(shù)分解與重構(gòu)算法如下：

分解問題是：已知系數(shù)dj,n，求系數(shù)dj+1,2n和dj+1,2n+1。分解公式為[15,20-22]

(1)

重構(gòu)問題是：已知系數(shù)dj+1,2n和dj+1,2n+1，求系數(shù)dj,n。重構(gòu)公式為[15,20-22]

(2)

圖3 小波包分解樹(3層)Fig.3 Wavelet packet decomposition tree at level 3

圖4 小波包分解單元Fig.4 Unit of wavelet packet decomposition

1.2 小波包分解系數(shù)的頻率順序

對于信號分析來講，重要的是將小波包系數(shù)按頻率順序排列，而不是按小波包節(jié)點(diǎn)的自然順序排列。由于小波包分解時，每對高頻系數(shù)進(jìn)行一次分解，相應(yīng)分解系數(shù)的頻率排位順序就會“翻轉(zhuǎn)”一次，從而造成頻率順序與節(jié)點(diǎn)自然順序不一致現(xiàn)象[15,21]。按此規(guī)律，可推算出小波包分解系數(shù)頻率順序與對應(yīng)小波包節(jié)點(diǎn)號的對應(yīng)關(guān)系。以3層小波包分解為例，按前所述規(guī)律，可推算出小波包分解系數(shù)頻率由小到大排位與節(jié)點(diǎn)序號的對應(yīng)關(guān)系，具體分析結(jié)果見表1。由表1進(jìn)一步分析可推算小波包分解系數(shù)頻率由小到大對應(yīng)節(jié)點(diǎn)號的順序，如第3分解層的小波包分解系數(shù)頻率由小到大對應(yīng)節(jié)點(diǎn)號的順序?yàn)?3,0)，(3,1)，(3,3)，(3,2)，(3,6)，(3,7)，(3,5)，(3,4)[15,21]。

表1 小波包系數(shù)頻率排位與節(jié)點(diǎn)序號(3層)Table 1 The “frequency” order of wavelet packet coefficients and natural nodes order at level 3

1.3 小波包系數(shù)的閾值處理

小波包系數(shù)閾值處理的關(guān)鍵是閾值估計(jì)和閾值施加方案。小波包分析中有4種常用的閾值估計(jì)方法和3種簡單的閾值施加方案，具體如下：

4種閾值估計(jì)方法[11,13-15,17,21]：① 自適應(yīng)閾值(Rigrsure)；② 固定形式閾值("sqtwolog")；③ 啟發(fā)式閾值(Heursure)；④ 極小化極大閾值(Minimax)。

3種閾值施加方案[19]：① 硬閾值處理(Hard Thresholding)；② 軟閾值處理(Soft Thresholding)；③ 比例閾值處理(Percentage Thresholding)。

除上述方法和方案外，還有多種閾值估計(jì)方法和閾值施加方案。具體采用何種閾值估計(jì)方法和閾值施加方案取決于實(shí)際的應(yīng)用，后續(xù)將結(jié)合航空重力測量數(shù)據(jù)處理做進(jìn)一步探討研究。

2 小波低通濾波器設(shè)計(jì)

基于小波包分析方法，按估算的小波包分解層次對應(yīng)的信號頻率范圍、低通濾波器的截止頻率和小波包系數(shù)頻率節(jié)點(diǎn)的排列順序，優(yōu)化小波包樹，設(shè)計(jì)小波低通濾波器。小波低通濾波器的設(shè)計(jì)需重點(diǎn)關(guān)注3個方面：小波的選取、小波包樹的優(yōu)化設(shè)計(jì)和小波包系數(shù)的閾值量化處理。

2.1 小波的選取

正交或雙正交小波包分解可將信號按頻率分解到無重疊的子帶上，易于小波包快速算法的實(shí)現(xiàn)。因此，在小波包分析中，小波濾波器通常是選用正交或雙正交小波來實(shí)現(xiàn)的。常用的正交、雙正交小波有Discrete Meyer小波(dmey)、Daubenchies小波系(dbN)、Symlets小波系(symN)、Coiflet小波系(coifN)和Biorthogonal小波系(biorNr.Nd)等[17]。

不同類型的小波系有不同的特性，其主要特性見表2[17]。具有緊支撐集的小波(時域或空域)，局部化能力強(qiáng)，易于算法實(shí)現(xiàn)；對稱性使小波濾波器具有線性相位，避免了信號失真；消失矩的大小決定了用小波逼近光滑函數(shù)的收斂率；正則性是與光滑性相關(guān)聯(lián)的，正則性越好，重構(gòu)信號就越光滑，而光滑性決定了濾波頻率分辨率的高低[10,17,24]。

Daubechies原理表明除Haar小波以外不存在對稱的緊支撐正交小波[18]，且緊支撐性與光滑性二者不可兼得[10]。構(gòu)造一個有正交性、緊支撐集、平滑性和對稱性的小波是困難的[10]，應(yīng)根據(jù)應(yīng)用實(shí)際需求，在各項(xiàng)特性指標(biāo)間做出取舍，選擇適合的小波。

表2 常用正交或雙正交小波系的主要特性Table 2 Wavelet families and main associated properties

2.2 小波包樹的優(yōu)化設(shè)計(jì)

信號的多尺度小波包分解相當(dāng)于帶通濾波(分頻)，對于第j小波包分解層，每個節(jié)點(diǎn)分解系數(shù)的頻帶寬度為2-jHz(采樣頻率1 Hz)[17,20,24]。根據(jù)濾波器的截止頻率和估算的分解層次對應(yīng)的信號頻率范圍，選取分解層次N，并優(yōu)化設(shè)計(jì)小波樹。如周期為60 s的低通濾波器，相對應(yīng)截止頻率fc近似為 0.016 7 Hz(注fc≈2-6+2-10+2-14≈0.016 7 Hz)，取分解層次N=14，按小波包系數(shù)頻率由小到大對應(yīng)的小波包節(jié)點(diǎn)號順序，優(yōu)化設(shè)計(jì)小波包分解樹，如圖5所示。

2.3 小波包系數(shù)閾值量化處理方案

采用正交或雙正交小波濾波器對數(shù)據(jù)進(jìn)行濾波處理，可將數(shù)據(jù)按頻率劃分為若干子帶。按小波包系數(shù)頻率大小，分別采用不同的閾值準(zhǔn)則和處理方案，對小波包系數(shù)進(jìn)行量化處理?，F(xiàn)以60 s低通濾波器為例，加以具體說明。

按1.2小節(jié)所述方法分析可知：小波包低頻系數(shù)頻率由小到大對應(yīng)的小波包節(jié)點(diǎn)號順序?yàn)?6,0)，(10,24)，(14,408)，高頻系數(shù)頻率由小到大對應(yīng)的小波包節(jié)點(diǎn)號順序?yàn)?14,409)，(13,205)，(12,103)，(11,50)，(9,13)，(8,7)，(7,2)，(5,1)，(4,1)，(3,1)，(2,1)，(1,1)。按照小波包系數(shù)頻率大小，分別采用不同的閾值方案進(jìn)行處理。以異常信號為主的低頻段，閾值要??；而以噪聲為主的高頻段，閾值要大。具體方案如下：

圖5 小波低通濾波器的小波包分解樹(濾波周期：60 s)Fig.5 Wavelet packet decomposition tree of the low pass filter (filtering period: 60 s)

方案1：保留全部低頻系數(shù)，高頻系數(shù)全部都置“零”；

方案2：節(jié)點(diǎn)(6,0)低頻系數(shù)進(jìn)行平滑處理，節(jié)點(diǎn)(10,24)和(14,408)低頻系數(shù)按估算比例壓縮或采用通用閾值估計(jì)進(jìn)行量化處理，各層高頻系數(shù)全部都置“零”；

方案3：在方案2基礎(chǔ)上，節(jié)點(diǎn)(14,409)高頻系數(shù)(過度頻帶)采用通用閾值估計(jì)進(jìn)行量化處理，其余各節(jié)點(diǎn)高頻系數(shù)全部都置“零”；

方案4：在方案3基礎(chǔ)上，節(jié)點(diǎn)(13,205)，(12,103)，(11,50)，(9,13)高頻系數(shù)采用通用閾值估計(jì)進(jìn)行量化處理，其余各節(jié)點(diǎn)高頻系數(shù)全部都置“零”。

對于航空重力異常測量信號來講，由于小波包分解低頻系數(shù)同時包含信號和噪聲兩部分，且不同節(jié)點(diǎn)小波包系數(shù)的信噪比也不相同，故通用閾值處理方案有其局限性。在通用閾值處理方案基礎(chǔ)上，本文提出的小波包系數(shù)閾值處理方案，針對不同節(jié)點(diǎn)的小波包低頻系數(shù)采用平滑或按估算比例壓縮的方法壓制高頻(相對)噪聲，故能顯著提高濾波精度。

3 濾波試驗(yàn)

圖6為GT-1A系統(tǒng)3010測線原始未濾波航空自由空間重力異常測量數(shù)據(jù)，異常信號被幅度大致在±5 000 mGal的噪聲信號所淹沒，橫軸Fid為測線基準(zhǔn)點(diǎn)號，間距平均為30 m[2,7-9 ]。選用正交或雙正交小波，采用本文提出的閾值處理方案和設(shè)計(jì)的小波低通濾波器，對圖6原始未濾波數(shù)據(jù)進(jìn)行了低通濾波試驗(yàn)研究，并與GT-1A系統(tǒng)軟件濾波結(jié)果進(jìn)行對比分析。

圖6 GT-1A系統(tǒng)原始未濾波航空自由空間重力異常Fig.6 GT-1A raw unfiltered airborne gravity free air anomaly

3.1 小波包樹的優(yōu)化設(shè)計(jì)

在航空重力測量中，低通濾波器的截止頻率fc一般在0.005～0.016 6 Hz之間，也就是濾波周期為200～60 s[5]。試驗(yàn)中，小波低通濾波器的濾波周期選取具有代表性的60 s和100 s(對應(yīng)的截止頻率fc分別為0.016 7 Hz 和0.01 Hz)，優(yōu)化設(shè)計(jì)的小波包樹分別如圖5和圖7所示。

100 s與60 s小波低通濾波器的設(shè)計(jì)方法相同，其小波包低頻系數(shù)頻率由小到大對應(yīng)的小波包節(jié)點(diǎn)號順序?yàn)?7,0)，(9,6)，(12,60)，高頻系數(shù)頻率由小到大對應(yīng)的小波包節(jié)點(diǎn)號順序?yàn)?12,61)，(11,31)，(10,14)，(8,2)，(6,1)，(5,1)，(4,1)，(3,1)，(2,1)，(1,1)。

圖7 小波低通濾波器的小波包分解樹(濾波周期：100 s)Fig.7 Wavelet packet decomposition tree of the low pass filter (filtering period: 100 s)

3.2 小波的選取

不同類型的小波，具有不同的特性(表2)，濾波器小波的選取應(yīng)從正交性、對稱性、消失矩和正則性等多方面綜合考慮。對于dbN、symN、coifN和biorNr.Nd小波系，隨著小波階數(shù)N的增大，小波具有更高的正則性[24]，相應(yīng)也具有更高的頻率分辨率。但處于不同頻段的重力異常測量信號的信噪比不同，故濾波周期不同的濾波器對小波頻率分辨率的需求也不盡相同。試驗(yàn)中，60 s濾波器選取dmey、db7、sym7、coif5、bior5.5和bior6.8小波，而100 s濾波器選取dmey、db11、sym10、coif5、bior5.5和bior6.8小波，對GT-1A系統(tǒng)獲取的航空重力數(shù)據(jù)進(jìn)行濾波試驗(yàn)研究。對于60 s和100 s濾波器小波的選取， dbN、symN小波系存在差異，而其他小波系則相同(在小波系現(xiàn)有小波中，所選用小波的頻率分辨率已最高)。

首先，采用GT-1A系統(tǒng)數(shù)據(jù)處理軟件對圖6的原始未濾波數(shù)據(jù)進(jìn)行60 s和100 s低通濾波處理，并將濾波結(jié)果作為近似標(biāo)準(zhǔn)；其次，采用本文設(shè)計(jì)的60 s和100 s小波包分解樹，對GT-1A系統(tǒng)濾波結(jié)果分別進(jìn)行小波包分解計(jì)算；最后，分別將分解的部分小波包系數(shù)(表3或表4：本行以下各行小波包節(jié)點(diǎn)系數(shù))置“零”，再用剩余部分小波包系數(shù)(表3或表4：本行及以上各行小波包節(jié)點(diǎn)系數(shù))進(jìn)行重構(gòu)計(jì)算，重構(gòu)計(jì)算結(jié)果與GT-1A系統(tǒng)濾波結(jié)果的差值統(tǒng)計(jì)見表3和表4，表中數(shù)據(jù)為均方差值，單位為mGal。

3.3 小波濾波試驗(yàn)

選用不同的正交或雙正交小波(同前)，采用本文研究、設(shè)計(jì)的60 s和100 s小波低通濾波器和本文提出的閾值量化處理方案，對圖6所示的原始未濾波數(shù)據(jù)進(jìn)行了低通濾波試驗(yàn)研究。

濾波結(jié)果如圖8～19所示，表5和表6分別為60 s和100 s濾波結(jié)果與GT-1A系統(tǒng)濾波結(jié)果的差值統(tǒng)計(jì)，表中數(shù)據(jù)為均方差值，單位為mGal。

3.4 結(jié)果分析

由表3、表4分析可知：① 試驗(yàn)中，所選取的正交或雙正交小波(見表3和表4)，按優(yōu)化設(shè)計(jì)的小波包樹分解后，能夠精確完整重構(gòu)；② 所選用小波，因頻率分辨率和相位失真等因素的影響，引起重構(gòu)計(jì)算結(jié)果誤差；理想情況下，隨著頻率由低到高重構(gòu)包含小波包系數(shù)的增加，重構(gòu)計(jì)算結(jié)果誤差會變小，小波濾波重構(gòu)精度會提高。而濾波試驗(yàn)結(jié)果，隨著頻率由低到高重構(gòu)包含小波包系數(shù)增加，重構(gòu)計(jì)算結(jié)果誤差出現(xiàn)局部波動，總體趨勢和理想情況分析吻合。③由①和②進(jìn)一步分析可知：本文小波包分解樹的優(yōu)化設(shè)計(jì)方案是可行的。

表3 小波包重構(gòu)計(jì)算結(jié)果與GT-1A系統(tǒng)濾波結(jié)果的差值統(tǒng)計(jì)(60 s低通濾波器)Table 3 Difference statistics between the wavelet packet reconstruction’s computing result and GT-1A filtering result (60 s low-pass filter)

表4 小波包重構(gòu)計(jì)算結(jié)果與GT-1A系統(tǒng)濾波結(jié)果的差值統(tǒng)計(jì)(100 s低通濾波器)Table 4 Difference statistics between the wavelet packet reconstruction’s computing result and GT-1A filtering result (100 s low-pass filter)

圖8 dmey小波濾波器與GT-1A系統(tǒng)60 s濾波航空自由空間重力異常對比Fig.8 Airborne gravity free air anomaly of dmey wavelet and GT-1A 60 s filter

圖9 db7小波濾波器與GT-1A系統(tǒng)60 s濾波航空自由空間重力異常對比Fig.9 Airborne gravity free air anomaly of db7 wavelet and GT-1A 60 s filter

圖10 sym7小波濾波器與GT-1A系統(tǒng)60 s濾波航空自由空間重力異常對比Fig.10 Airborne gravity free air anomaly of sym7 wavelet and GT-1A 60 s filter

圖11 coif5小波濾波器與GT-1A系統(tǒng)60 s濾波航空自由空間重力異常對比Fig.11 Airborne gravity free air anomaly of coif5 wavelet and GT-1A 60 s filter

圖12 bior5.5小波濾波器與GT-1A系統(tǒng)60 s濾波航空自由空間重力異常對比Fig.12 Airborne gravity free air anomaly of bior5.5 wavelet and GT-1A 60 s filter

圖13 bior6.8小波濾波器與GT-1A系統(tǒng)60 s濾波航空自由空間重力異常對比Fig.13 Airborne gravity free air anomaly of bior6.8 wavelet and GT-1A 60 s filter

圖14 dmey小波濾波器與GT-1A系統(tǒng)100 s濾波航空自由空間重力異常對比Fig.14 Airborne gravity free air anomaly of dmey wavelet and GT-1A 100 s filter

圖15 db11小波濾波器與GT-1A系統(tǒng)100 s濾波航空自由空間重力異常對比Fig.15 Airborne gravity free air anomaly of db11 wavelet and GT-1A 100 s filter

圖16 sym10小波濾波器與GT-1A系統(tǒng)100 s濾波航空自由空間重力異常對比Fig.16 Airborne gravity free air anomaly of sym10 wavelet and GT-1A 100 s filter

圖17 coif5小波濾波器與GT-1A系統(tǒng)100 s濾波航空自由空間重力異常對比Fig.17 Airborne gravity free air anomaly of coif5 wavelet and GT-1A 100 s filter

圖18 bior5.5小波濾波器與GT-1A系統(tǒng)100 s濾波航空自由空間重力異常對比Fig.18 Airborne gravity free air anomaly of bior5.5 wavelet and GT-1A 100 s filter

圖19 bior6.8小波濾波器與GT-1A系統(tǒng)100 s濾波航空自由空間重力異常對比Fig.19 Airborne gravity free air anomaly of bior6.8 wavelet and GT-1A 100 s filter

表5 小波濾波結(jié)果與GT-1A系統(tǒng)濾波結(jié)果的差值統(tǒng)計(jì)(濾波周期：60 s)Table 5 Difference statistics between the wavelet filtering result and GT-1A filtering result (filtering period: 60 s)

表6 小波濾波結(jié)果與GT-1A系統(tǒng)濾波結(jié)果的差值統(tǒng)計(jì)(濾波周期：100 s)Table 6 Difference statistics between the wavelet filtering result and GT-1A filtering result (filtering period: 100 s)

對GT-1A系統(tǒng)原始未濾波數(shù)據(jù)(圖6)、及其經(jīng)GT-1A系統(tǒng)軟件60 s和100 s濾波獲得的數(shù)據(jù)分別進(jìn)行功率譜分析，結(jié)果如圖20所示(黑色、藍(lán)色和紅色曲線分別為原始未濾波、及其GT-1A系統(tǒng)60 s和100 s濾波航空自由空間重力異常功率譜)。由圖20分析可知：重力異常信號主要集中在小于 0.02 Hz低頻段，在0～0.01 Hz頻段的信噪比大于處于0.01～0.02 Hz頻段的信噪比。

由表3或表4可以看出處于不同頻帶的小波包節(jié)點(diǎn)系數(shù)對濾波結(jié)果的影響，以60 s濾波器為例，對本文采用的閾值量化處理方案做進(jìn)一步說明。方案1：保留低頻系數(shù)，高頻系數(shù)置“零”；方案2：在方案1基礎(chǔ)上，對低頻系數(shù)細(xì)化處理；節(jié)點(diǎn)(6,0)低頻系數(shù)，以異常信號為主，對其進(jìn)行平滑處理，抑制高頻(相對)噪聲干擾。節(jié)點(diǎn)(10,24)和(14,408)低頻系數(shù)，噪聲幅值遠(yuǎn)大于異常信號幅值，按估算比例壓縮小波包系數(shù)幅值，降低高頻(相對)噪聲的影響，亦可采用通用閾值估計(jì)進(jìn)行量化處理。方案3：在方案2基礎(chǔ)上，對節(jié)點(diǎn)(14,409)高頻(過度頻帶)系數(shù)，采用通用閾值估計(jì)進(jìn)行量化處理；因小波包分解存在頻率分辨率的問題，有可能造成頻率失真，通過方案3彌補(bǔ)方案2的不足。因節(jié)點(diǎn)(14,409)小波包系數(shù)對應(yīng)的頻帶寬度較“窄”，方案3對濾波效果的改善作用影響不大。方案4：在方案3基礎(chǔ)上，對節(jié)點(diǎn)(13,205)，(12,103)，(11,50)，(9,13)部分高頻系數(shù)，采用通用閾值估計(jì)進(jìn)行量化處理。對于高頻系數(shù)，因信噪比極低，故閾值量化處理效果有限。綜合以上分析可知：① 四種方案中，方案2或方案3濾波效果最佳；② 因在低頻段，不同頻帶信噪比存在差異，頻率越低，信噪比越大，故100 s濾波器的濾波效果好于60 s濾波器的濾波效果。

圖20 GT-1A系統(tǒng)原始未濾波、60 s和100 s濾波航空自由空間重力異常功率譜Fig.20 Power spectrums of GT-1A raw unfiltered, 60 s and 100 s filtered airborne gravity free air anomaly

濾波試驗(yàn)結(jié)果圖8～19、表5和表6也可進(jìn)一步驗(yàn)證上述分析結(jié)果。通過對比分析可知：① 采用方案2～4的濾波效果好于采用方案1的濾波效果；② 采用方案2～4的濾波效果基本相同，其中濾波效果最佳的為方案2或方案3；③ 總體而言，100 s濾波器的濾波效果好于60 s濾波器的濾波效果。

將濾波試驗(yàn)結(jié)果(最佳效果)按選用的小波歸納重新整理，不同小波60 s、100 s濾波結(jié)果分別如圖21、圖22所示，濾波結(jié)果與GT-1A系統(tǒng)濾波結(jié)果差值統(tǒng)計(jì)見表7。由圖21、圖22和表7分析可知：① 60 s小波濾波器：按濾波試驗(yàn)效果，選用小波的排列順序依次為bior5.5、bior6.8、sym7、db7、coif5和 dmey小波，dmey小波濾波效果最佳(均方差值約為0.28 mGal)；② 100 s小波濾波器：按濾波試驗(yàn)效果，選用小波的排列順序依次為bior5.5、bior6.8、sym10、coif5、db11和 dmey小波，dmey小波濾波效果最佳(均方差值約為0.22 mGal)。

濾波試驗(yàn)效果與選用小波的頻率分辨率和對稱性有關(guān)，根據(jù)所選小波頻率分辨率高低，再考慮小波的對稱性，綜合判斷小波濾波試驗(yàn)效果。在濾波試驗(yàn)選用的小波中，db7和db11小波無對稱性，sym7、sym10和coif5小波均具有準(zhǔn)對稱性，而dmey、bior5.5和bior6.8小波才具有完全對稱性。由表5和表6中方案1的濾波效果可判斷出所選用小波頻率分辨率的高低，亦可按小波頻率分辨率定義計(jì)算比較[24]。60 s濾波試驗(yàn)中，選用的coif5和 dmey小波頻率分辨率較高；100 s濾波試驗(yàn)中，選用的coif5、dmey、sym10和db11小波頻率分辨率較高。綜合小波頻率分辨率和對稱性二者因素，可判斷在60 s和100 s濾波試驗(yàn)中，選用dmey小波的濾波效果為最佳，這與試驗(yàn)結(jié)果與GT-1A結(jié)果差值統(tǒng)計(jì)表5和表6數(shù)據(jù)相吻合，也與相關(guān)文獻(xiàn)的研究結(jié)果一致[24]。

圖21 小波濾波器與GT-1A系統(tǒng)60 s濾波航空自由空間重力異常對比Fig.21 Airborne gravity free air anomaly of wavelets and GT-1A 60 s filter

圖22 小波濾波器與GT-1A系統(tǒng)100 s濾波航空自由空間重力異常對比Fig.22 Airborne gravity free air anomaly of wavelets and GT-1A 100 s filter

表7 小波濾波試驗(yàn)結(jié)果與GT-1A系統(tǒng)濾波結(jié)果的差值統(tǒng)計(jì)Table 7 The difference between wavelet filtering and GT-1A filtering result and statistics

4 結(jié)論

1) 基于小波包系數(shù)頻率順序，本文提出的小波包分解樹的優(yōu)化設(shè)計(jì)方案是可行的。根據(jù)濾波器的截止頻率(濾波周期的倒數(shù))、估算的分解層次對應(yīng)的信號頻率范圍和小波包系數(shù)頻率由小到大節(jié)點(diǎn)的排列順序，優(yōu)化設(shè)計(jì)小波包分解樹，為小波濾波器的設(shè)計(jì)奠定基礎(chǔ)。

2) 選用正交或雙正交小波，本文研究、設(shè)計(jì)的小波低通濾波器，以及提出的閾值處理方案是可行的。濾波試驗(yàn)：60 s和100 s濾波結(jié)果與GT-1A系統(tǒng)濾波結(jié)果的均方差值分別達(dá)到約0.28 mGal和0.22 mGal，獲得與GT-1A系統(tǒng)幾乎同樣滿意的濾波效果。

3) 在低頻段，處于不同頻帶(節(jié)點(diǎn))的小波包系數(shù)信噪比不同，應(yīng)采用不同的方案進(jìn)行量化處理。在本文提出的4種閾值量化處理方案中，方案2或方案3的濾波效果最佳(相對GT-1A系統(tǒng)濾波結(jié)果)。

4) 濾波效果與選用小波的頻率分辨率和對稱性有關(guān)。在濾波試驗(yàn)選中的不同類型的小波中，選用dmey小波的濾波器濾波效果最佳(相對GT-1A系統(tǒng))；選用bior5.5小波的濾波器，方案2相較于方案1的濾波效果改善幅度相對最大(見表5)。

與離散Meyer小波不同，B樣條小波是雙正交的、局部緊支撐、且同樣也具有良好的對稱性和平滑性。如何構(gòu)造出比dmey小波具有更高頻率分辨率的B樣條小波，進(jìn)一步提高小波濾波器的濾波效果，有待于后續(xù)深入研究。