GPS地表監測數據去噪技術

李衛東

(河北金山礦業有限公司總經理, 河北 邢臺市 054100)

0 引 言

西郝莊礦礦區進行了多年開采，目前地表存有以前露天開采遺留的露天坑及堆存巖土，地面已出現多處塌陷及較大范圍的裂縫。為此，設立西郝莊鐵礦采空區GPS地表沉降監測系統，應用現代化的傳感和網絡通訊技術對地表沉降進行實時監測，及時掌握地表沉降信息，為地壓災害預測提供依據。

由于各種干擾因素的作用，GPS監測數據中的有用信號和噪聲的時域和頻域特性通常是不一樣的，有用信號在時域和頻域上具有局部化特征，表現為低頻特性或是一些比較平穩的信號；而噪聲在時域與頻域空間中的分布具有全局性特征，它在整個數據監測與數據采集的時域內處處存在，在頻域上表現為高頻特性。根據此特性，對其可以通過限定閥值形式對小波系數進行處理，從而降低噪聲部分值，再通過信號重構恢復信號，實現有用信號和噪聲的分離，實現降噪的目標。

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

1 去噪小波的構 造方法

1.1 基本原理

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

其中，當ω=π，F0(eiω)≠0，且|F0(eiω)|在ω=0～2π范圍內的上界≤2p-1。

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

1.2 Daubechies緊支撐正交小波的構造方法

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

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

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

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

其中，當ω=π，F0(eiω)≠0。

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

當ω=π，z=-1，F0(-1)≠0。

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

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

1.3 緊支撐正交小波的構造方法的實現形式

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

總結起來，可得

2 監測數據去噪實現技術

2.1 去噪實現過程

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

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

第1種是硬閥值形式，即把含噪的小波系數的絕對值與所選定的閥值λ進行比較，將小于等于閥值λ的點變為0，大于閥值λ的點保持不變。

硬閥值法具有徹底去噪的特點，同時也使得信號失去了有用成分。

第2種是軟閥值形式，即把含噪的小波系數與所選定的閥值λ進行比較，將大于閥值的點收縮為該點值與閥值的差；將小于閥值相反數的點擴展為該點值與閥值的和；將絕對值小于等于閥值的點變為0。

在處理GPS監測數據時，從第一層到第J層，每一層選擇一個軟閥值。閥值λ的計算公式為：

這樣具有更好的去噪效果。λ對于不同的信號具有自適應性，去噪徹底，同時也不會去除有用信號，N是數據個數；J是分解層數；σ為監測數據的標準方差，計算方法如下：

式中，median為去中值函數，{wJ-1，k，k=1，2，…，2J-1}為小波細節系數。

第3種是強制去噪法，即強制所有分解層下的細節系數置全為0。

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

第4種是Garrote去噪法，即把含噪的小波系數的絕對值與所選定的閥值λ進行比較，將小于等于閥值λ的點變為0，而對大于閥值λ的點按如下公式計算：

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

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

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

2.2 去噪效果分析

假設所收集到的GPS監測數據為序列為:

s1，s2，…，sN

該序列含有噪聲。通過采用小波方法進行去噪，得到去噪后的GPS監測數據序列為

v1，v2，…，vN

則殘差為

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

2.2.1 殘差統計指標

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

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

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

(6) 頻度：mf={r1，r2，…，rN}中出現次數最多的那個值。

(11) 偏態系數：

(12)峰態系數：

(13) 變異系數：

(14)自相關系數：

有偏自相關系數：

無偏自相關系數：

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

2.2.2 殘差直方圖。

將殘差序列{r1，r2，…，rN}分成若干個長度相等的區間，統計每個區間內殘差值出現的次數，生成殘差直方圖。從殘差直方圖能看出殘差的分布。

2.2.3 殘差FFT分析

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

F=fft(r，N)

當k=0，1，…，N-1時，

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

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

相位：α=angle(F)。

3 各監測點監測數據凈變化量去噪

GPS01監測數據去噪過程如圖1～3所示。

圖1 GPS01點z方向監測數據的分解

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

圖2 GP01點z方向監測數據的去噪

4 結 論

本文介紹了GPS監測數據去噪幾種常用技術的原理及其實現過程。欲達到良好的去噪效果，選擇合適的技術參數是關鍵。

圖3 GP01點z方向監測數據的殘值分析