基于小波變換的ICESAT-GlAS波形處理

邱 賽，邢艷秋，李立存，王 萌

(東北林業(yè)大學(xué)森林作業(yè)環(huán)境研究中心，哈爾濱150040)

搭載在冰、云和陸地高程衛(wèi)星 (the ice，cloud，and land elevation satellite，ICESat)上的地學(xué)激光測高系統(tǒng) (geoscience laser altimeter system，GLAS)獲取的原始波形數(shù)據(jù)由于受系統(tǒng)噪聲和大氣的影響存在著較大的白噪聲，另外ICESATGLAS在以植被分布為主要特征的激光腳點(diǎn)中不僅包含植被信息，還包括地形陡變處、枯落木、石塊等地物信息，造成其波形數(shù)據(jù)失真、疊加［1］，因此需要對其進(jìn)行降噪處理。小波變換1是一種時(shí)間窗和頻率窗都改變的時(shí)頻局域化分析法，既具備短時(shí)傅里葉變換局部化的思想，又克服了窗口大小不隨頻率變化的缺點(diǎn)，具有自適應(yīng)性［2-5］，因此小波變換在信號去噪中具有很好的應(yīng)用前景。用小波變換法的去噪效果在很大程度上取決于小波基，因?yàn)椴⒉皇撬械男〔ɑ歼m合于激光雷達(dá)數(shù)據(jù)信號的處理，即使是同一信號用不同的小波基進(jìn)行處理也會產(chǎn)生不同的結(jié)果。小波基的選擇目前沒有統(tǒng)一的標(biāo)準(zhǔn)，本文是通過比較分析小波基的5個(gè)主要參數(shù)(正交性、對稱性、正則性、消失矩、支撐長度)來確定選用哪種小波基，并以吉林省汪清林業(yè)局經(jīng)營區(qū)的ICESAT-GLAS波形數(shù)據(jù)為研究對象，進(jìn)行去噪效果的比較分析。

1 數(shù)據(jù)收集及其預(yù)處理

1.1 ICESAT-GLAS 數(shù)據(jù)

本文采用的ICESAT-GLAS數(shù)據(jù)是從美國國家冰雪數(shù)據(jù)中心 (NSIDC)(http://nsidc.org/data/icesat)下載的。ICESAT-GLAS數(shù)據(jù)產(chǎn)品可分為3個(gè)級別:0、1和2級，共15種產(chǎn)品，本文采用的是1A級數(shù)據(jù)中的GLA01和GLA14兩種數(shù)據(jù)產(chǎn)品。其中GLA01中包括光斑號、光斑幀號、時(shí)間、脈沖強(qiáng)度 (光電壓值)、經(jīng)緯度等信息，GLA14記錄了地表的高度信息，包括高程、大地經(jīng)緯度、回波振幅、面積和標(biāo)準(zhǔn)差。本次研究利用的數(shù)據(jù)坐標(biāo)是WGS84坐標(biāo)。

1.2 波形數(shù)據(jù)預(yù)處理

(1)ICESAT-GLAS原始波形數(shù)據(jù)的數(shù)據(jù)格式為整數(shù)型二進(jìn)制格式，在MATLAB和IDL平臺上將其轉(zhuǎn)化為十進(jìn)制的ASCII(American Standard Code for Information Interchange，ASCII)格式 (*.rasc)［6］。

(2)根據(jù)激光光斑腳點(diǎn)的經(jīng)緯度坐標(biāo)將GLA01數(shù)據(jù)與GLA14的相關(guān)波形信息合成，包括光電壓值、高程、接收與發(fā)送的UTC時(shí)間等波形信息與參數(shù)。

(3)將數(shù)據(jù)分別保存到文本文檔中，再加載到MATLAB文件中，并轉(zhuǎn)存為為.mat文件的格式，然后對其進(jìn)行去噪處理。

2 波形數(shù)據(jù)去噪

2.1 小波變換原理

小波，是一種小區(qū)域的波，小波變換的實(shí)質(zhì)是將L2(R)空間中的任意函數(shù)f(t)表示成為其在小波核函數(shù)Ψa，b(t)之上的投影的疊加，即小波變換是信號f(t)和Ψa，b(t)的內(nèi)積，公式為:

公式 (1)中a，b，t均為連續(xù)變量。

Ψa，b(t)是由給定的小波基Ψ (t)在伸縮因子a和平移因子b下進(jìn)行伸縮和平移得到函數(shù)族:

ab∈R;a〉0 其中a為伸縮因子，b為平移因子。

2.2 小波基的的選擇

小波變換中采用的小波基函數(shù)Ψ(t)具有多樣性，因此用不同的小波基函數(shù)分析同一問題會有不同的結(jié)果。本文通過對小波基的主要參數(shù)［7］(正交性、對稱性、正則性、消失矩、支撐長度)進(jìn)行比較分析初步判定兩種小波基 (Daubechies、Symlets)對波形數(shù)據(jù)進(jìn)行處理是較好的。本文選取的是db1小波基和sym7小波基進(jìn)行比較。表1為常用小波基所對應(yīng)重要參數(shù)特點(diǎn)。

表1 常用小波參數(shù)特性Tab.1 Common wavelet basis parameters'characteristics

2.3 小波去噪過程

本文假設(shè)激光雷達(dá)波形信號f(t)為離散采樣后的得到的544點(diǎn)的離散信號f(n)，n=1，2，……544。對信號去噪的具體過程如下:

(1)一維小波分解

本文選擇分解為4層，將信號分解成不同頻率通道成分，并將每個(gè)頻率通道成分按相位進(jìn)行了分解，獲取小波分解高頻系數(shù)WΨf(a，b)，如公式(3)所示:

{aN-2k}和 {bN-2k}是分解序列，N為采樣點(diǎn)數(shù)，j表示分解的層數(shù)，Cj=

(2)對小波高頻系數(shù)進(jìn)行閾值化處理

D.L.Dohono提出閾值化處理方法包括硬閾值法和軟閾值法，一般，硬閾值比軟閾值處理后的信號粗糙，所以本文采用軟閾值法處理。具體公式如 (4)所示:

Wδ為量化后的小波系數(shù)值，w為分解后的小波系數(shù)，sgn(·)是符號函數(shù)

用閾值法去噪最關(guān)鍵的是確定閾值δ，本文采用的是Birge-Massart懲罰函數(shù)方法來由小波系數(shù)選擇規(guī)則得到［8］。

(3)一維小波重構(gòu)

根據(jù)小波分解后的最底層低頻系數(shù)和各層經(jīng)過閾值量化處理后的高頻系數(shù)進(jìn)行一維小波重構(gòu)。具體如公式 (5)所示:

{pk-2N}和 {qk-2N}是重建序列，N為采樣點(diǎn)數(shù)，j表示分解的層數(shù)，Cj=Cjk。

3 去噪效果評價(jià)

本文采用信噪比 (SNR)和均方根誤差(RMSE)做為評價(jià)指標(biāo)，對同一波形數(shù)據(jù)去噪效果進(jìn)行評價(jià)［9］。具體如公式 (6)、 (7)所示:s(i)為含噪信號，f(i)為去噪后信號，L為信號長度。

表2 小波變換波形數(shù)據(jù)處理結(jié)果Tab.2 Results of wavelet transform

圖1 分別用db1和Sym7對波形進(jìn)行去噪處理Fig.1 Denoised wave by db1 and sym7

圖1中 (a)為用db1對數(shù)據(jù)去噪，(b)為用sym7對數(shù)據(jù)去噪，用此兩種小波對GLA01_22587635:28去噪結(jié)果表明經(jīng)過兩種小波基去噪后波形都比原始波形平滑了，但經(jīng)db1小波基去噪后的波形在波峰前后尚含有一定噪聲，而經(jīng)Sym7小波基去噪后的波形效果就比db1小波基效果好，波形更加平滑;同時(shí)表2中顯示經(jīng)db1小波基處理后的信噪比為 30.129 26，均方根誤差為 1.234 002，經(jīng)Sym7小波基去噪后的信噪比為32，760 89，均方根誤差為0.912 325。觀察比較表2中的所有數(shù)據(jù)發(fā)現(xiàn)對同一組數(shù)據(jù)Sym7小波基的信噪比比db1小波基的要高，而均方根誤差比db1小波基的要低。信噪比越高，均方根誤差越低說明信號的去噪效果越好。所以在處理大光斑激光雷達(dá)波形數(shù)據(jù)方面Sym7小波基比db1小波基更合適。

4 結(jié)論

本文對吉林長白山汪清林業(yè)局的經(jīng)營區(qū)的大光斑激光雷達(dá)波形數(shù)據(jù)，通過比較分析發(fā)現(xiàn)db1小波基和Sym7小波基的去噪效果比db1小波基的好，更適合對大光斑激光雷達(dá)波形數(shù)據(jù)進(jìn)行去噪處理。原始波形數(shù)據(jù)中的細(xì)節(jié)信號是高頻突變信號，屬于噪聲，逐漸逼近于信號的低頻部分，而小波變換恰恰可以起到低通濾波的作用，因此既可以有效保留信號頻譜成分，又去除了噪聲，克服了波形“重疊”現(xiàn)象。

本實(shí)驗(yàn)只是在數(shù)據(jù)中隨機(jī)選取了5組數(shù)據(jù)采用里兩種小波基進(jìn)行研究，不排除偶然性，在今后研究中應(yīng)該盡量對多組數(shù)據(jù)進(jìn)行處理，降低偶然性誤差，提高準(zhǔn)確性和可信性;小波基有很多種，本文只是比較了其中兩種小波，不排除其他小波基能更好的對數(shù)據(jù)進(jìn)行去噪，比如高斯小波。

［1］劉經(jīng)南，張小紅.利用激光強(qiáng)度信息分類激光掃描測高數(shù)據(jù)［J］.武漢大學(xué)學(xué)報(bào)(信息科學(xué)版)，2005，34(6):696-700.

［2］Lázaro J C.Influence of thresholding procedures in ultrasonic grain noise reduction using wavelets［J］.Ultrasonics，2002，40(8):263-267.

［3］Xu Y S，Weaver JB，Hedly D M.Wavelet transform domain filters:A spatially selective noise filtration technique［J］.IEEE Transactions on Image Processing，1994，3(6):747-758.

［4］Da silva E A D，Ghanbari M.On the performance of linear phase wavelet transforms in low bit-rate image coding［J］.IEEE Transactions on Image Processing，1996，5(5):689-705.

［5］孫仁山，李文彬，徐凱宏.基于小波變換的林業(yè)圖像處理研究［J］.森林工程，2005，21(1):4-6.

［6］邢艷秋，王立海.基于ICESAT-GLAS完整波形的坡地森林冠層高度反演研究——以吉林長白山林區(qū)為例［J］.武漢大學(xué)學(xué)報(bào)(信息科學(xué)版).2009，34(6):696-700.

［7］Antonini M，Barlaud M，Mathieu P，et al.Image codingusing wavelet transform［J］.IEEE Transactions on Image Processing，1992，1(2):205-220.

［8］胡 波，陳 懇，徐建瑜.一種基于新型小波包閾值的圖像去噪方法［J］.寧 波 大 學(xué) 學(xué) 報(bào)(理工版).2009，22(4):454-458.

［9］Labat D.Recent advances in wavelet analyses:Part 1.A review of concepts［J］.Journal of Hydrology，2005，314(1-4):275-288.