王江濤，張 昆，李 歡

(中煤航測遙感集團(tuán)有限公司，西安 710100)

0 引言

利用GPS定位技術(shù)進(jìn)行變形監(jiān)測時，所獲得的變形數(shù)據(jù)序列包括系統(tǒng)變形、各種觀測誤差以及噪聲，最終反映為變形體的位移、噪聲和粗差。特別是對微變形體進(jìn)行監(jiān)測時，GPS測量誤差成為強(qiáng)噪聲，從受強(qiáng)噪聲干擾的序列觀測數(shù)據(jù)中提取微弱的變形信息，是GPS動態(tài)監(jiān)測應(yīng)解決的一個關(guān)鍵問題[1]。此外，隨著動態(tài)變形監(jiān)測方法的革新，動態(tài)變形監(jiān)測數(shù)據(jù)源的特性發(fā)生了變化，變形信息數(shù)據(jù)的采樣頻率得到大幅提高，數(shù)據(jù)量成倍增長，傳感器采集數(shù)據(jù)的特點又決定了變形信息必然要表現(xiàn)出強(qiáng)污染特征[2]。因此，有必要對動態(tài)變形監(jiān)測數(shù)據(jù)進(jìn)行綜合處理，去除或削弱各種誤差和噪聲的影響，獲取真實的動態(tài)變形信息。

因此本文結(jié)合小波分析具有良好的時頻局部性和多尺度分解能力的特點，提出了一種基于小波分析的動態(tài)變形信息提取多尺度分析方法，實現(xiàn)對信號的濾波去噪、特征提取、多路徑消除等功能，這為信號濾波、信噪分離及強(qiáng)噪聲背景下弱信號特征的提取提供了有效途徑。

1 小波系數(shù)與信號、噪聲的關(guān)系

1.1 動態(tài)變形信息的含噪模型及其小波變換

設(shè)變形監(jiān)測系統(tǒng)所獲得的數(shù)據(jù)有兩部分組成，令s(t)為變形信號，n(t)為隨機(jī)噪聲，即n(t)～N(0,σ2)，受污染的輸出信號為x(t)，則監(jiān)測變形信號的基本模型可以描述為[3]：

x(t)=s(t)+n(t)

(1)

在有用信號s(t)中，既有可能是實際變形信號sd(t)或確定性噪聲sn(t)，也有可能是兩者的混合[4]，即：

s(t)=sd(t)+sn(t)

(2)

但是，這些信號都具有不同的時頻特性。因此，可以借助小波變換的多尺度分析，有效地對不同頻率成分進(jìn)行分離，最終達(dá)到降噪的目的。

小波變換是一種線性變換，變形信號s(t)和噪聲n(t)這兩個信號經(jīng)變換后的離散近似信號和離散細(xì)節(jié)信號分別是各個信號變換后的近似信號和細(xì)節(jié)信號之和，具有如下關(guān)系：

(3)

(4)

=WTs(j,k)+WTn(j,k)

(5)

1.2 小波系數(shù)與信號特征的關(guān)系

小波變換的實質(zhì)是把原始信號不同頻率段的信息經(jīng)過多尺度分解抽取出來，并反映在時間軸上，體現(xiàn)出信號的時頻特性[4]。近似信號多反映在大尺度變換后的低頻部分，細(xì)節(jié)信號多反映在小尺度變換后的高頻部分。

信號的過零點、極值點以及過零間隔等是描述信號的一些重要特征。在變形監(jiān)測系統(tǒng)中，信號的振蕩性使得信號會發(fā)生急劇變化，即信號的奇異性，這是分析變形的關(guān)鍵之處。信號的奇異一般可分為兩種情況：①信號在某一時刻的突變引起信號的不連續(xù)(稱為第一類間斷點)；②高階奇異信號(稱為第二類間斷點)，即信號的外觀很光滑，其幅值沒有突變，但其K階導(dǎo)數(shù)存在突變點。由信號的小波變換的奇異點在多尺度上的綜合表現(xiàn)來表示信號(特變是它的突變或瞬態(tài)特征)是小波變換引人注意的另一應(yīng)用領(lǐng)域[5]。在數(shù)學(xué)上，常采用李氏指數(shù)表示函數(shù)局部光滑和奇異特征，且α越大表明某點光滑度越好，越小表明某點奇異性越大。

1992年Mallat等人在經(jīng)過深入研究給出了小波變換與李氏指數(shù)之間的關(guān)系，參照文獻(xiàn)[6]，給出小波系數(shù)與信號特征之間的關(guān)系如下：

設(shè)0≤α≤1，函數(shù)x(t)在[a,b]上有一致李氏指數(shù)α的充要條件是存在一個常數(shù)ξ>0，使得存在t∈[a,b]，二進(jìn)小波變換滿足：

|WTj,kx(t)|≤ξ2jα

(6)

兩邊取對數(shù)后變?yōu)椋?/p>

log2|WTj,kx(t)|≤log2ξ+jα

(7)

jα這一項給出了小波變換的對數(shù)值隨尺度j或α的變化規(guī)律，同時也把小波變換的尺度特征j與李氏指數(shù)α聯(lián)系起來了?？梢缘弥?，如果信號的s(t)的李氏指數(shù)α>0，則該信號的小波變換的極大值將隨尺度j的增大而增大；反之，若α<0，則隨尺度j的增大而減小[7]。因此，可以通過小波系數(shù)的極大值與尺度間的關(guān)系來描述信號的特征。

2 小波分析降噪過程

一般地，對一維動態(tài)變形信號進(jìn)行降噪的過程，可以歸納為以下幾步：

①原始信號的小波分解：選擇并確定一個小波基及其分解層數(shù)N，這是信號分解好壞的關(guān)鍵，不同小波基適用情況有所差異，分解層數(shù)也影響著提取精度。

②進(jìn)行閾值量化：對小波分解得到的細(xì)節(jié)系數(shù)進(jìn)行閾值(如硬閾值、軟閾值)量化處理；

③一維信號的小波重建過程：將降噪處理后的系數(shù)通過小波進(jìn)行重構(gòu)，得到消噪后的變形數(shù)據(jù)序列估計值。

上述②中閾值量化處理，是小波降噪的核心，直接影響信號降噪的質(zhì)量。在小波變換中，對各層系數(shù)降噪所需的閾值一般是根據(jù)原信號的信號噪聲比來選取的，閾值量化主要分硬閾值和軟閾值兩種方式。實驗表明，硬閾值處理后的信號比軟閾值處理后的信號更粗糙。

3 降噪模型性能分析

為了便于定量比較不同小波基和不同閾值量化方式的降噪效果，定義原信號為x(t)和去噪后的估計信號為x(t)，則濾波后噪聲部分的均方根誤差(NRMSE)為：

(8)

信噪比(RSN)是測量信號中噪聲強(qiáng)度的傳統(tǒng)方法，將其定義為：

(9)

其中n為離散信號的長度。

能量成分(EC)是指消噪后的信號在原信號中的能量比，將其定義為：

EC=norm(x(t))/norm(x(t))

(10)

上述中，均方誤差NRMSE越大，信噪比RSN越高，則表明估計的信號越接近原始信號，降噪效果越好，但前提是盡可能多的保留原始信號的有用成分。能量成分EC越大表明損失原信號的能量越大，但對于強(qiáng)噪聲信號而言，損失的能量主要是強(qiáng)噪聲引起的，因此這種損失可認(rèn)為是有益的。

4 計算工具及仿真算例分析

4.1 強(qiáng)噪聲背景下弱信號特征提取算例

設(shè)x(t)為一附加強(qiáng)噪聲的弱信號，噪聲方差σ=20，表達(dá)式如下：

x(t) =s(t)+n(t)=-8(t+30)/2 047.5+

4sin(2πft)+N(0,20)2

(11)

信號的頻率f=0.05Hz。對該信號進(jìn)行時間間隔為0.05 s的采樣，即采樣頻率為20Hz，采樣點數(shù)為211=2 048個點。含強(qiáng)噪聲的信號如圖1。

圖1 附加強(qiáng)噪聲的弱信號Figure 1 Weak signal with strong noise added

由已知條件和多尺度分析，在頻率帶寬為10Hz進(jìn)行小波分解，各個尺度上對應(yīng)的頻率帶寬分別為：

原始信號20：[0,10Hz]；

21尺度水平：低頻段a1[0,5Hz]+高頻段d1[5 Hz,10Hz]；

22尺度水平：低頻段a2[0,2.5Hz]+高頻段d2[2.5Hz,5Hz]；

23尺度水平：低頻段a3[0,1.25Hz]+高頻段d3[1.25Hz,2.5Hz]；

24尺度水平：低頻段a4[0,0.625Hz]+高頻段d4[0.625Hz,1.25Hz]；

…… ……

27尺度水平：低頻段a7[0, 0.078 125Hz] + 高頻段d7[0.078 125 Hz, 0.156 25Hz]。

由已知采樣條件可知，弱信號的頻率0.05Hz落在第7分解尺度[0, 0.078 125Hz]上，因此確定小波分解尺度為7層。

對不同的數(shù)據(jù)來源和處理目的，應(yīng)選擇不同的小波基函數(shù)，只有這樣才能使處理結(jié)果效果達(dá)到最佳[8]。根據(jù)原信號的性質(zhì)，選擇了正則性和擴(kuò)展性都比較好的dbN小波基函數(shù)，為了確定N的取值，分別取N=4，5，6，8，10時，用dbN小波對信x(t)號進(jìn)行7層分解，分解結(jié)果如圖2所示。

為了評價各種小波基分解的效果，采用上文的三個指標(biāo)進(jìn)行評價，如表1：

表1 不同小波基分解性能比較

圖2 dbN小波分解結(jié)果Figure 2 Results of dbN wavelet decomposition

從上表可以看出，db8小波分解效果較其它小波基較好，因此選擇了db8進(jìn)行小波分解。為了在提取特征信息的同時，有效的消除強(qiáng)噪聲，通常的方法有強(qiáng)制消噪法和給定閾值消噪法。強(qiáng)制消噪法是指將小波分解結(jié)構(gòu)中的細(xì)節(jié)系數(shù)強(qiáng)制置為零，然后對信號進(jìn)行小波重構(gòu)；而給定閾值消噪法是通過經(jīng)驗公式給出閾值進(jìn)行消噪并在此基礎(chǔ)上重構(gòu)信號。本文選擇強(qiáng)制消噪方法，選取db8小波7層濾波后重構(gòu)了信號，重構(gòu)信號與原信號如3圖所示。

圖3 db8做7層濾波信號Figure 3 Wavelet db8 7 layered signal filtering

從圖3可以看出，消噪后的信號能基本實現(xiàn)原信號的恢復(fù)，很好地反映了原信號的趨勢，且周期性明顯。消除的噪聲的標(biāo)準(zhǔn)差為19.75左右，這跟所加噪聲20非常接近，說明小波濾波效果很好。

4.2 信號奇異性檢驗

奇異性可分為兩種，即第一類間斷點和第二類間斷點，本文將對第一類型的間斷點進(jìn)行檢測分析。如圖4所示，給定的信號由三個分段的余弦信號組成，各個信號的頻率不同，目的是為了檢驗信號的幅值在局部時間上的突變，即信號的奇異性。分段余弦信號表達(dá)式為：

(12)

由于實際測量得到的信號往往受到測量過程和不可避免的噪聲影響而被模糊化，為了討論小波變換能在噪聲影響下能否檢測到信號奇異點，在原信號基礎(chǔ)上添加了一幅值為[-2,2]的隨機(jī)白噪聲，形成一個附加噪聲的幅值突變信號x(t)。噪聲信號如5圖，含噪聲后的突變余弦信號如6圖。

分別采用Daubechies小波db1，db2，db3，db5，db8對無噪信號s(t)和含噪信號x(t)的奇異性進(jìn)行比較檢測。 圖7為db2和db3對兩個信號的小波分解結(jié)果，檢測結(jié)果如表2所示。

圖4 幅值突變的余弦信號Figure 4 Amplitude mutational cosine signal

圖5 噪聲信號Figure 5 Noise signal

圖6 附加噪聲后的突變余弦信號Figure 6 Noise added mutational cosine signal

表2 奇異性檢測結(jié)果

(注：為識別出奇異點的最佳細(xì)節(jié)級別)

由表2可知，在無噪聲的情況下，db1～db8小波都能較好地在d1，d2細(xì)節(jié)上檢測到突變點的瞬間位置，而且奇異點定位性能非常好。且由圖7顯示，db2小波定位精度較db3更高，消失矩階數(shù)N越小，定位性能越好。在有噪聲情況下，較低細(xì)節(jié)上主要是噪聲信號，低尺度上基本只能檢測到第二個頻率突變點(因為第二個突變點頻率變化較大，特征明顯)，信號的所有突變點只能在較高尺度上檢測到，如尺度d4，d5上可以看到所有突變情況。雖然具有隨機(jī)噪聲的影響，但db1～db8都能清楚地檢測到信號的突變點，說明小波的去噪能力還是不錯的。但需注意的是，隨著小波消失矩階數(shù)N的增大，檢測到所有突變點的尺度也越來越大，奇異點定位能力逐漸下降。這說明，選擇較小消失矩的的小波在檢測此類突變點時效果更好，有利于精確定位奇異點。

5.試驗結(jié)果及分析

本試驗將采用于山東鮑店煤礦副井樓頂所測GPS-RTK動態(tài)變形監(jiān)測數(shù)據(jù)為例進(jìn)行降噪研究。數(shù)據(jù)于2015年3月13日09時53分到16時02分觀測，采樣間隔為1 s，信號頻率分析后為0.005Hz。數(shù)據(jù)用單歷元解算法解算得到變形量。分析中選取4 096個歷元，由于樓頂變形量較小，對所選取數(shù)據(jù)進(jìn)行適當(dāng)放大處理。X方向變形量如圖8所示。

數(shù)據(jù)采集時環(huán)境為：RTK基準(zhǔn)站架設(shè)在煤礦工業(yè)廣場上，移動站通過強(qiáng)制對中安置在觀測墩上，副井工作室有升降機(jī)在每隔幾分鐘進(jìn)行一次升降，致使副井出現(xiàn)不同強(qiáng)度的振動，總體振動幅度還是比較小的?；鶞?zhǔn)站和移動站間距大概250 m。

選擇dbN(N=4,5,6,8,10)小波基對X方向變形量進(jìn)行小波分解，分解尺度為7，分解結(jié)果如圖9所示，性能評價如表3所示。

表3 不同小波基分解性能比較

圖 7 奇異檢測小波分解結(jié)果Figure 7 Singularity test wavelet decomposed results

圖8 變形數(shù)據(jù)X方向變形量Figure 8 Deformation data X axis deformation amount

圖9 dbN(N=3，4，5，6，8，10)小波分解結(jié)果Figure 9 Wavelet decomposed results of dbN (N=3, 4, 5, 6, 8 and 10)

圖10 db10強(qiáng)制消噪重構(gòu)信號Figure 10 Reconstructed signals of db10 forced noise reduction

比較不同評價指標(biāo)，選擇db10小波基作為小波分解函數(shù)。通過強(qiáng)制消噪后重構(gòu)的信號能基本體現(xiàn)井筒在X方向的變形趨勢，如圖10所示。

由圖9、圖10和表3可以看出：①對同一信號，N取值不同對信號降噪效果差異不是很大，但從更多保留真實信號的角度來說，db10在對該信號消噪上效果更好；②由于副井樓頂動態(tài)變形是微弱的，信號中沒有大的突變點，因此在dbN小波分解中沒有明顯發(fā)現(xiàn)奇異點出現(xiàn)的時刻；③通過對原信號采用強(qiáng)制消噪后，重構(gòu)的信號能很好的反映副井振動變形趨勢。

6 結(jié)論

①本文針對GPS變形監(jiān)測獲得的受強(qiáng)干擾的變形數(shù)據(jù)序列提取微變形這一關(guān)鍵難題，在深入研究小波分析在變形信號降噪中的作用后，提出了基于小波分析的動態(tài)變形信息提取多尺度分析方法，由于小波分析具有良好的時頻局部性和多尺度分解的能力，這為信號濾波、信噪分離及強(qiáng)噪聲背景下弱信號特征的提取提供了有效途徑。

②通過仿真分析和實測數(shù)據(jù)，研究了如何在強(qiáng)噪聲背景下提取弱信號特征和對含噪信號進(jìn)行奇異性檢查，這為研究動態(tài)變形趨勢、特征提取和周期性分析提供指導(dǎo)作用。

③通過引用能量成分、濾波后噪聲部分的均方誤差和信噪比等指標(biāo)來評價各種小波降噪性能。分析了如何將小波分析應(yīng)用到對強(qiáng)噪聲背景下弱信號的提取和對信號的奇異性檢驗中。發(fā)現(xiàn)選擇消失矩階數(shù)較大的基本小波能更有效地提取強(qiáng)噪聲背景下弱信號的特征。在對信號進(jìn)行奇異性檢驗時，選擇較小消失矩的的小波在檢測突變點時效果更好，更能精確定位奇異點。

參考文獻(xiàn):

[1]楊玉龍.GPS變形監(jiān)測技術(shù)的現(xiàn)狀及發(fā)展趨勢[J].測繪科學(xué), 2014(2): 319-319.

[2]黃清,唐詩華,許虹偉.基于小波去噪的變形監(jiān)測數(shù)據(jù)處理的方法研究[J].北京測繪,2015(6):31-33.

[3]譚超.基于小波去噪—灰色預(yù)測混合模型在變形監(jiān)測中的應(yīng)用研究[D].陜西:西北大學(xué),2016.

[4]沙成滿,韓合新,楊冬梅.基于小波去噪的改進(jìn)灰色自適應(yīng)變形預(yù)測[J].東北大學(xué)學(xué)報(自然科學(xué)版),2011, 32 (8):1195-1197.

[5]侯瑞.基于小波變換降噪方法的穩(wěn)定性比較模型[D].中國地質(zhì)大學(xué)(北京),2012.

[6]Donoho D L. No-linear Wavelet Methods for Recovery of Signal, Densities and Spectra-form Indirect and Noisy Data[C]. Proceeding of Symposia in Applied Mathematics and Stanford Report,1993.

[7]鄧敏,吳光強(qiáng).基于變異系數(shù)的小波降噪綜合評價法[M].機(jī)電一體化, 2016,3(6):30-35.

[8]向東,貢建兵.變形序列小波去噪最佳分解尺度量化指標(biāo)的確定[J].武漢大學(xué)學(xué)報(信息科學(xué)版), 2014,39(4):468-470.