航空重力測量數據的小波濾波處理

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

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

0 引言

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

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

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

1 小波包閾值濾波

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

1.1 信號的小波包分解與重構

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

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

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

選用正交或雙正交小波，則小波包系數分解與重構算法如下：

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

(1)

重構問題是：已知系數dj+1,2n和dj+1,2n+1，求系數dj,n。重構公式為[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 小波包分解系數的頻率順序

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

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

1.3 小波包系數的閾值處理

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

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

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

除上述方法和方案外，還有多種閾值估計方法和閾值施加方案。具體采用何種閾值估計方法和閾值施加方案取決于實際的應用，后續將結合航空重力測量數據處理做進一步探討研究。

2 小波低通濾波器設計

基于小波包分析方法，按估算的小波包分解層次對應的信號頻率范圍、低通濾波器的截止頻率和小波包系數頻率節點的排列順序，優化小波包樹，設計小波低通濾波器。小波低通濾波器的設計需重點關注3個方面：小波的選取、小波包樹的優化設計和小波包系數的閾值量化處理。

2.1 小波的選取

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

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

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

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

2.2 小波包樹的優化設計

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

2.3 小波包系數閾值量化處理方案

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

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

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

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

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

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

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

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

3 濾波試驗

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

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

3.1 小波包樹的優化設計

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

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

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

3.3 小波濾波試驗

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

濾波結果如圖8～19所示，表5和表6分別為60 s和100 s濾波結果與GT-1A系統濾波結果的差值統計，表中數據為均方差值，單位為mGal。

3.4 結果分析

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

表3 小波包重構計算結果與GT-1A系統濾波結果的差值統計(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 小波包重構計算結果與GT-1A系統濾波結果的差值統計(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系統60 s濾波航空自由空間重力異常對比Fig.8 Airborne gravity free air anomaly of dmey wavelet and GT-1A 60 s filter

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

表7 小波濾波試驗結果與GT-1A系統濾波結果的差值統計Table 7 The difference between wavelet filtering and GT-1A filtering result and statistics

4 結論

1) 基于小波包系數頻率順序，本文提出的小波包分解樹的優化設計方案是可行的。根據濾波器的截止頻率(濾波周期的倒數)、估算的分解層次對應的信號頻率范圍和小波包系數頻率由小到大節點的排列順序，優化設計小波包分解樹，為小波濾波器的設計奠定基礎。

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

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

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

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