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

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

0 引言

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

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

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

1 小波包閾值濾波

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

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

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

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

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

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

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

(1)

重構(gòu)問(wèn)題是：已知系數(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ù)的頻率順序

對(duì)于信號(hào)分析來(lái)講，重要的是將小波包系數(shù)按頻率順序排列，而不是按小波包節(jié)點(diǎn)的自然順序排列。由于小波包分解時(shí)，每對(duì)高頻系數(shù)進(jìn)行一次分解，相應(yīng)分解系數(shù)的頻率排位順序就會(huì)“翻轉(zhuǎn)”一次，從而造成頻率順序與節(jié)點(diǎn)自然順序不一致現(xiàn)象[15,21]。按此規(guī)律，可推算出小波包分解系數(shù)頻率順序與對(duì)應(yīng)小波包節(jié)點(diǎn)號(hào)的對(duì)應(yīng)關(guān)系。以3層小波包分解為例，按前所述規(guī)律，可推算出小波包分解系數(shù)頻率由小到大排位與節(jié)點(diǎn)序號(hào)的對(duì)應(yīng)關(guān)系，具體分析結(jié)果見表1。由表1進(jìn)一步分析可推算小波包分解系數(shù)頻率由小到大對(duì)應(yīng)節(jié)點(diǎn)號(hào)的順序，如第3分解層的小波包分解系數(shù)頻率由小到大對(duì)應(yīng)節(jié)點(diǎn)號(hào)的順序?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)序號(hào)(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種簡(jiǎn)單的閾值施加方案，具體如下：

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é)合航空重力測(cè)量數(shù)據(jù)處理做進(jìn)一步探討研究。

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

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

2.1 小波的選取

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

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

Daubechies原理表明除Haar小波以外不存在對(duì)稱的緊支撐正交小波[18]，且緊支撐性與光滑性二者不可兼得[10]。構(gòu)造一個(gè)有正交性、緊支撐集、平滑性和對(duì)稱性的小波是困難的[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ì)

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

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

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

按1.2小節(jié)所述方法分析可知：小波包低頻系數(shù)頻率由小到大對(duì)應(yīng)的小波包節(jié)點(diǎn)號(hào)順序?yàn)?6,0)，(10,24)，(14,408)，高頻系數(shù)頻率由小到大對(duì)應(yīng)的小波包節(jié)點(diǎn)號(hào)順序?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)行處理。以異常信號(hào)為主的低頻段，閾值要??；而以噪聲為主的高頻段，閾值要大。具體方案如下：

圖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ù)(過(guò)度頻帶)采用通用閾值估計(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ù)全部都置“零”。

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

3 濾波試驗(yàn)

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

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

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

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

100 s與60 s小波低通濾波器的設(shè)計(jì)方法相同，其小波包低頻系數(shù)頻率由小到大對(duì)應(yīng)的小波包節(jié)點(diǎn)號(hào)順序?yàn)?7,0)，(9,6)，(12,60)，高頻系數(shù)頻率由小到大對(duì)應(yīng)的小波包節(jié)點(diǎn)號(hào)順序?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)從正交性、對(duì)稱性、消失矩和正則性等多方面綜合考慮。對(duì)于dbN、symN、coifN和biorNr.Nd小波系，隨著小波階數(shù)N的增大，小波具有更高的正則性[24]，相應(yīng)也具有更高的頻率分辨率。但處于不同頻段的重力異常測(cè)量信號(hào)的信噪比不同，故濾波周期不同的濾波器對(duì)小波頻率分辨率的需求也不盡相同。試驗(yàn)中，60 s濾波器選取dmey、db7、sym7、coif5、bior5.5和bior6.8小波，而100 s濾波器選取dmey、db11、sym10、coif5、bior5.5和bior6.8小波，對(duì)GT-1A系統(tǒng)獲取的航空重力數(shù)據(jù)進(jìn)行濾波試驗(yàn)研究。對(duì)于60 s和100 s濾波器小波的選取， dbN、symN小波系存在差異，而其他小波系則相同(在小波系現(xiàn)有小波中，所選用小波的頻率分辨率已最高)。

首先，采用GT-1A系統(tǒng)數(shù)據(jù)處理軟件對(duì)圖6的原始未濾波數(shù)據(jù)進(jìn)行60 s和100 s低通濾波處理，并將濾波結(jié)果作為近似標(biāo)準(zhǔn)；其次，采用本文設(shè)計(jì)的60 s和100 s小波包分解樹，對(duì)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小波低通濾波器和本文提出的閾值量化處理方案，對(duì)圖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é)果誤差會(huì)變小，小波濾波重構(gòu)精度會(huì)提高。而濾波試驗(yàn)結(jié)果，隨著頻率由低到高重構(gòu)包含小波包系數(shù)增加，重構(gòu)計(jì)算結(jié)果誤差出現(xiàn)局部波動(dòng)，總體趨勢(shì)和理想情況分析吻合。③由①和②進(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濾波航空自由空間重力異常對(duì)比Fig.8 Airborne gravity free air anomaly of dmey wavelet and GT-1A 60 s filter

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

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

圖11 coif5小波濾波器與GT-1A系統(tǒng)60 s濾波航空自由空間重力異常對(duì)比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濾波航空自由空間重力異常對(duì)比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濾波航空自由空間重力異常對(duì)比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濾波航空自由空間重力異常對(duì)比Fig.14 Airborne gravity free air anomaly of dmey wavelet and GT-1A 100 s filter

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

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

圖17 coif5小波濾波器與GT-1A系統(tǒng)100 s濾波航空自由空間重力異常對(duì)比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濾波航空自由空間重力異常對(duì)比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濾波航空自由空間重力異常對(duì)比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)

對(duì)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分析可知：重力異常信號(hào)主要集中在小于 0.02 Hz低頻段，在0～0.01 Hz頻段的信噪比大于處于0.01～0.02 Hz頻段的信噪比。

由表3或表4可以看出處于不同頻帶的小波包節(jié)點(diǎn)系數(shù)對(duì)濾波結(jié)果的影響，以60 s濾波器為例，對(duì)本文采用的閾值量化處理方案做進(jìn)一步說(shuō)明。方案1：保留低頻系數(shù)，高頻系數(shù)置“零”；方案2：在方案1基礎(chǔ)上，對(duì)低頻系數(shù)細(xì)化處理；節(jié)點(diǎn)(6,0)低頻系數(shù)，以異常信號(hào)為主，對(duì)其進(jìn)行平滑處理，抑制高頻(相對(duì))噪聲干擾。節(jié)點(diǎn)(10,24)和(14,408)低頻系數(shù)，噪聲幅值遠(yuǎn)大于異常信號(hào)幅值，按估算比例壓縮小波包系數(shù)幅值，降低高頻(相對(duì))噪聲的影響，亦可采用通用閾值估計(jì)進(jìn)行量化處理。方案3：在方案2基礎(chǔ)上，對(duì)節(jié)點(diǎn)(14,409)高頻(過(guò)度頻帶)系數(shù)，采用通用閾值估計(jì)進(jìn)行量化處理；因小波包分解存在頻率分辨率的問(wèn)題，有可能造成頻率失真，通過(guò)方案3彌補(bǔ)方案2的不足。因節(jié)點(diǎn)(14,409)小波包系數(shù)對(duì)應(yīng)的頻帶寬度較“窄”，方案3對(duì)濾波效果的改善作用影響不大。方案4：在方案3基礎(chǔ)上，對(duì)節(jié)點(diǎn)(13,205)，(12,103)，(11,50)，(9,13)部分高頻系數(shù)，采用通用閾值估計(jì)進(jìn)行量化處理。對(duì)于高頻系數(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é)果。通過(guò)對(duì)比分析可知：① 采用方案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)效果與選用小波的頻率分辨率和對(duì)稱性有關(guān)，根據(jù)所選小波頻率分辨率高低，再考慮小波的對(duì)稱性，綜合判斷小波濾波試驗(yàn)效果。在濾波試驗(yàn)選用的小波中，db7和db11小波無(wú)對(duì)稱性，sym7、sym10和coif5小波均具有準(zhǔn)對(duì)稱性，而dmey、bior5.5和bior6.8小波才具有完全對(duì)稱性。由表5和表6中方案1的濾波效果可判斷出所選用小波頻率分辨率的高低，亦可按小波頻率分辨率定義計(jì)算比較[24]。60 s濾波試驗(yàn)中，選用的coif5和 dmey小波頻率分辨率較高；100 s濾波試驗(yàn)中，選用的coif5、dmey、sym10和db11小波頻率分辨率較高。綜合小波頻率分辨率和對(duì)稱性二者因素，可判斷在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濾波航空自由空間重力異常對(duì)比Fig.21 Airborne gravity free air anomaly of wavelets and GT-1A 60 s filter

圖22 小波濾波器與GT-1A系統(tǒng)100 s濾波航空自由空間重力異常對(duì)比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ù))、估算的分解層次對(duì)應(yīng)的信號(hào)頻率范圍和小波包系數(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的濾波效果最佳(相對(duì)GT-1A系統(tǒng)濾波結(jié)果)。

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

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