GPS地表監(jiān)測(cè)數(shù)據(jù)去噪技術(shù)

李衛(wèi)東

(河北金山礦業(yè)有限公司總經(jīng)理, 河北 邢臺(tái)市 054100)

0 引 言

西郝莊礦礦區(qū)進(jìn)行了多年開(kāi)采，目前地表存有以前露天開(kāi)采遺留的露天坑及堆存巖土，地面已出現(xiàn)多處塌陷及較大范圍的裂縫。為此，設(shè)立西郝莊鐵礦采空區(qū)GPS地表沉降監(jiān)測(cè)系統(tǒng)，應(yīng)用現(xiàn)代化的傳感和網(wǎng)絡(luò)通訊技術(shù)對(duì)地表沉降進(jìn)行實(shí)時(shí)監(jiān)測(cè)，及時(shí)掌握地表沉降信息，為地壓災(zāi)害預(yù)測(cè)提供依據(jù)。

由于各種干擾因素的作用，GPS監(jiān)測(cè)數(shù)據(jù)中的有用信號(hào)和噪聲的時(shí)域和頻域特性通常是不一樣的，有用信號(hào)在時(shí)域和頻域上具有局部化特征，表現(xiàn)為低頻特性或是一些比較平穩(wěn)的信號(hào)；而噪聲在時(shí)域與頻域空間中的分布具有全局性特征，它在整個(gè)數(shù)據(jù)監(jiān)測(cè)與數(shù)據(jù)采集的時(shí)域內(nèi)處處存在，在頻域上表現(xiàn)為高頻特性。根據(jù)此特性，對(duì)其可以通過(guò)限定閥值形式對(duì)小波系數(shù)進(jìn)行處理，從而降低噪聲部分值，再通過(guò)信號(hào)重構(gòu)恢復(fù)信號(hào)，實(shí)現(xiàn)有用信號(hào)和噪聲的分離，實(shí)現(xiàn)降噪的目標(biāo)。

由于Daubechies小波具有很好的正交性、緊支性(緊支撐長(zhǎng)度為2N-1)、正則性、近似對(duì)稱性和適用Mallat快速算法等特點(diǎn)，其時(shí)域分析優(yōu)良，在監(jiān)測(cè)信號(hào)的奇異性中具有很大優(yōu)勢(shì)，適用于地表監(jiān)測(cè)數(shù)據(jù)的處理，且在許多工程領(lǐng)域中應(yīng)用比較成熟。dbN中N越大，其濾波的長(zhǎng)度越長(zhǎng)，其中db4具有最短的時(shí)窗和最高的時(shí)間分辨率，去噪效果良好。下面著重介紹Daubechies去噪小波的基本原理。

1 去噪小波的構(gòu) 造方法

1.1 基本原理

構(gòu)造去噪小波的充要條件是：對(duì)于離散濾波器h={h1，h2，…，hN}，兩尺度方程

其中，當(dāng)ω=π，F(xiàn)0(eiω)≠0，且|F0(eiω)|在ω=0～2π范圍內(nèi)的上界≤2p-1。

若{h1，h2，…，hN}滿足上述條件(1)～(3)，則φ(t)只存在一個(gè)解，并且φ(t)是緊支撐正交尺度函數(shù)，其支撐區(qū)間[0，N]。若令gn=(-1)nh1-n，則

1.2 Daubechies緊支撐正交小波的構(gòu)造方法

若ψ(t)有p階消失矩，則對(duì)于任何p-1次多項(xiàng)式x(t)=a0+a1t+…+ap-1tp-1，都有〈ψ(t),x(t)〉=0，即ψ(t)與任何p-1次多項(xiàng)式都是正交的。由于

將f(t+2-jk)在2-jk處展開(kāi)，得

由于ψ(t)有p階消失矩，則

(1) 小波函數(shù)ψ(t)具有p階消失矩；

其中，當(dāng)ω=π，F(xiàn)0(eiω)≠0。

取上式的z變換，z=eiω，即

當(dāng)ω=π，z=-1，F(xiàn)0(-1)≠0。

構(gòu)造Daobechies緊支撐正交小波濾波器的關(guān)鍵是：合理選擇F0(eiω)或F0(z)，可選方案為：

Daobechies小波的支撐[-p+1，p+1]，尺度函數(shù)的支撐[0，2p-1]。

1.3 緊支撐正交小波的構(gòu)造方法的實(shí)現(xiàn)形式

h0+(-1)h1e-iω+(-1)2h2+…+(-1)NhN=0

(-1)(-i)kh1+(-1)2(-2i)kh2+…+(-1)N(-Ni)khN=0

h1-2kh2+…+(-1)N-1NkhN=0，k=1，2，…，p-1

總結(jié)起來(lái)，可得

2 監(jiān)測(cè)數(shù)據(jù)去噪實(shí)現(xiàn)技術(shù)

2.1 去噪實(shí)現(xiàn)過(guò)程

(1) 選擇小波函數(shù)和小波分解層數(shù)，計(jì)算監(jiān)測(cè)數(shù)據(jù)到第N層的小波分解。小波分解層數(shù)的選擇也同樣影響著數(shù)據(jù)去噪的效果，若分解層次較少,則可能去噪不徹底，若分解層數(shù)太多，則影響計(jì)算速度，同時(shí)使得部分有用信號(hào)消除，且當(dāng)分解到一定程度后，濾波后的信號(hào)幾乎保持不變。

(2) 高頻系數(shù)的閥值選擇。閥值選取分為以下6種形式。

第1種是硬閥值形式，即把含噪的小波系數(shù)的絕對(duì)值與所選定的閥值λ進(jìn)行比較，將小于等于閥值λ的點(diǎn)變?yōu)?，大于閥值λ的點(diǎn)保持不變。

硬閥值法具有徹底去噪的特點(diǎn)，同時(shí)也使得信號(hào)失去了有用成分。

第2種是軟閥值形式，即把含噪的小波系數(shù)與所選定的閥值λ進(jìn)行比較，將大于閥值的點(diǎn)收縮為該點(diǎn)值與閥值的差；將小于閥值相反數(shù)的點(diǎn)擴(kuò)展為該點(diǎn)值與閥值的和；將絕對(duì)值小于等于閥值的點(diǎn)變?yōu)?。

在處理GPS監(jiān)測(cè)數(shù)據(jù)時(shí)，從第一層到第J層，每一層選擇一個(gè)軟閥值。閥值λ的計(jì)算公式為：

這樣具有更好的去噪效果。λ對(duì)于不同的信號(hào)具有自適應(yīng)性，去噪徹底，同時(shí)也不會(huì)去除有用信號(hào)，N是數(shù)據(jù)個(gè)數(shù)；J是分解層數(shù)；σ為監(jiān)測(cè)數(shù)據(jù)的標(biāo)準(zhǔn)方差，計(jì)算方法如下：

式中，median為去中值函數(shù)，{wJ-1，k，k=1，2，…，2J-1}為小波細(xì)節(jié)系數(shù)。

第3種是強(qiáng)制去噪法，即強(qiáng)制所有分解層下的細(xì)節(jié)系數(shù)置全為0。

{wj-1，k，j=1，2，…，J；k=1，2，…，2j-1}

第4種是Garrote去噪法，即把含噪的小波系數(shù)的絕對(duì)值與所選定的閥值λ進(jìn)行比較，將小于等于閥值λ的點(diǎn)變?yōu)?，而對(duì)大于閥值λ的點(diǎn)按如下公式計(jì)算：

第5種是Birge-Massart閾值去噪法，給定一個(gè)指定的分解層數(shù)j，對(duì)j+1以及更高層所有系數(shù)保留；對(duì)i層(1≤i≤j)，保留絕對(duì)值最大的ni個(gè)系數(shù)并確定：

式中，α為經(jīng)驗(yàn)系數(shù)，在信號(hào)降噪情況下取α=3；缺省情況下取S=L(1)，即第1層分解后系數(shù)的長(zhǎng)度，一般情況下，S滿足L(1)≤S≤2L(1)。設(shè)對(duì)信號(hào)進(jìn)行i層分解，每層保留的小波系數(shù)個(gè)數(shù)最多為n，由于每層保留的小波系數(shù)個(gè)數(shù)不同，系數(shù)較少的層數(shù)采取補(bǔ)零的方法使每層系數(shù)個(gè)數(shù)相同，從而形成小波變換值矩陣。

第6種是Semisoft闡值函數(shù)，即把含噪的小波系數(shù)的絕對(duì)值與所選定的閥值λ1和λ2進(jìn)行比較，按如下公式計(jì)算小波系數(shù)：

2.2 去噪效果分析

假設(shè)所收集到的GPS監(jiān)測(cè)數(shù)據(jù)為序列為:

s1，s2，…，sN

該序列含有噪聲。通過(guò)采用小波方法進(jìn)行去噪，得到去噪后的GPS監(jiān)測(cè)數(shù)據(jù)序列為

v1，v2，…，vN

則殘差為

ri=si-vii=1，2，…，N

2.2.1 殘差統(tǒng)計(jì)指標(biāo)

(1) 最大值:max=max{r1，r2，…，rN}。

(2) 最小值：min=min{r1，r2，…，rN}。

(3) 極差：range=max-min。

(6) 頻度：mf={r1，r2，…，rN}中出現(xiàn)次數(shù)最多的那個(gè)值。

(11) 偏態(tài)系數(shù)：

(12)峰態(tài)系數(shù)：

(13) 變異系數(shù)：

(14)自相關(guān)系數(shù)：

有偏自相關(guān)系數(shù)：

無(wú)偏自相關(guān)系數(shù)：

計(jì)算上述自相關(guān)系數(shù)的函數(shù)是autocorr(r)，若用xcov(r)計(jì)算，則所得結(jié)果是對(duì)稱的，即ACF-j=ACFjACF-N+1,ACF-N+2,…,ACF-1,ACF0,ACF1,…,ACFN-1

2.2.2 殘差直方圖。

將殘差序列{r1，r2，…，rN}分成若干個(gè)長(zhǎng)度相等的區(qū)間，統(tǒng)計(jì)每個(gè)區(qū)間內(nèi)殘差值出現(xiàn)的次數(shù)，生成殘差直方圖。從殘差直方圖能看出殘差的分布。

2.2.3 殘差FFT分析

對(duì)殘差序列r={r1，r2，…，rN}進(jìn)行FFT分析的方法是：

F=fft(r，N)

當(dāng)k=0，1，…，N-1時(shí)，

頻率：f=k/(NΔt)，Δt為取樣時(shí)間間隔；

幅值：A=2|F|/N=2abs(F)/N；

相位：α=angle(F)。

3 各監(jiān)測(cè)點(diǎn)監(jiān)測(cè)數(shù)據(jù)凈變化量去噪

GPS01監(jiān)測(cè)數(shù)據(jù)去噪過(guò)程如圖1～3所示。

圖1 GPS01點(diǎn)z方向監(jiān)測(cè)數(shù)據(jù)的分解

其他方向的去噪方法一樣，在此不再贅述。

圖2 GP01點(diǎn)z方向監(jiān)測(cè)數(shù)據(jù)的去噪

4 結(jié) 論

本文介紹了GPS監(jiān)測(cè)數(shù)據(jù)去噪幾種常用技術(shù)的原理及其實(shí)現(xiàn)過(guò)程。欲達(dá)到良好的去噪效果，選擇合適的技術(shù)參數(shù)是關(guān)鍵。

圖3 GP01點(diǎn)z方向監(jiān)測(cè)數(shù)據(jù)的殘值分析