優(yōu)化EWT-NLM的自適應(yīng)GNSS高程時間序列降噪方法

魯鐵定 陶 蕊 程遠明 周子琪 何錦亮

1 東華理工大學(xué)測繪工程學(xué)院，南昌市廣蘭大道418號，330013 2 南昌市城市規(guī)劃設(shè)計研究總院，南昌市春暉路599號，330038

如何有效削弱原始時序信號中的各種噪聲影響一直是GNSS時序分析的研究熱點之一。錢文龍等[1]將改進的經(jīng)驗?zāi)B(tài)分解(empirical mode decomposition, EMD)算法應(yīng)用于GPS信號降噪，效果良好，但會出現(xiàn)一定的端點效應(yīng)及模態(tài)混疊現(xiàn)象；張恒璟等[2]改進了基于連續(xù)均方誤差準(zhǔn)則的EEMD噪聲識別方法，可以正確識別信號與噪聲的分界點，但不具備較高的自適應(yīng)性；Tao等[3]應(yīng)用經(jīng)驗小波變換(empirical wavelet transform, EWT)進行GNSS數(shù)據(jù)降噪，可以較好地反映原始信號的細節(jié)特征，降噪效果明顯，且不會發(fā)生模態(tài)混疊及端點現(xiàn)象，但無法自適應(yīng)選擇閾值去除高頻噪聲。針對以上問題，本文設(shè)計一種引入樣本熵(sample entropy, SE)優(yōu)化的EWT-NLM算法應(yīng)用于GNSS坐標(biāo)時間序列降噪中，并通過模擬數(shù)據(jù)和實測數(shù)據(jù)驗證其有效性。

1 原理與方法

首先對原始信號進行EWT分解，獲取一系列經(jīng)驗小波分量和經(jīng)驗尺度分量，以SE為有效信息分量與混合噪聲分量的分界標(biāo)準(zhǔn)，分別提取低頻信息分量與中高頻混合噪聲分量，再對混合噪聲分量進行NLM濾波，最后重構(gòu)有效信息分量和濾波分量得到最終降噪信號。

1.1 EWT

EWT算法的核心思想是將原始信號進行傅里葉變換，分割所獲取的傅里葉頻譜，對不同的分界頻率構(gòu)建經(jīng)驗小波，確定經(jīng)驗尺度函數(shù)和經(jīng)驗小波函數(shù)后，再進行反傅里葉變換，最后獲取具有緊支撐傅里葉頻譜的調(diào)幅-調(diào)頻(AM-FM)成分，具體原理見文獻[4-5]。其簡要過程如下：1)對原始時序信號進行傅里葉變換，求得在支撐區(qū)間(0，π)的傅里葉頻譜F(ω)；2)對傅里葉頻譜F(ω)作自適應(yīng)分割，分解為N個頻帶與N-1個分界頻率；3)根據(jù)各分界頻率構(gòu)建經(jīng)驗小波φ(ω)，確定經(jīng)驗尺度函數(shù)和經(jīng)驗小波函數(shù)，再對F(ω)φ(ω)進行反傅里葉變換，獲取各模態(tài)分量。

原始GNSS坐標(biāo)信號x(t)進行EWT分解后，分解為1+N個經(jīng)驗?zāi)B(tài)分量，其中包含1個代表信號整體趨勢變化的經(jīng)驗尺度分量和N個代表不同信號頻域特征的經(jīng)驗小波分量：

1.2 SE

EWT算法不能準(zhǔn)確確定以噪聲為主導(dǎo)的經(jīng)驗小波分量界限，而SE可以反映非平穩(wěn)信號的復(fù)雜程度，熵值大表示原始信號中的噪聲信號占比高，即噪聲信號特征明顯；熵值小表示原始信號的趨勢性、周期性規(guī)律較好，即噪聲信號特征不明顯。故引入SE來自適應(yīng)劃分信息分量的閾值，將其作為判別所有分量中信息分量與混合分量、噪聲分量的標(biāo)準(zhǔn)，可以有效識別GNSS高程時間序列信號分量的界限。具體計算過程如下[6]：

1)將時間跨度為L的原始時間序列信號x(i) 構(gòu)建為m維向量X(i)：

X(i)=[x(i)x(i+1) …x(i+m-1)],

i=1,2,…,N-m-1

(2)

0

(3)

4)令m=m+1,重復(fù)上述步驟得Cm+1(r)。

5)因時間跨度L為有限數(shù)，SE值可表示為：

可以看出，參數(shù)m和r的設(shè)定會直接影響熵值。據(jù)已有研究，一般m取2,r取原始信號標(biāo)準(zhǔn)差的0.1～0.25倍，考慮到實測GNSS信號的復(fù)雜性，本文取r=0.25σ(X)[6]。

1.3 NLM濾波

從原始信號中提取出信息分量后，有部分真實信息湮沒在混合經(jīng)驗小波分量及噪聲小波分量之中，利用NLM算法可以有效地對剩余信號進行降噪，從而獲取殘余真實信息。NLM算法利用原始信號的自相似特性進行濾波，其設(shè)置一個檢索跨度窗口，在此跨度內(nèi)盡可能多地搜尋相似值，然后運用加權(quán)平均值獲取實際GNSS殘差噪聲中的真實值，即可以反映原始信號的細節(jié)變化，達到降噪的目的[7]。假設(shè)噪聲信號為x(i)，降噪后信號為(i)，則有：

式中，ω表示兩點之間的相似度，即權(quán)重，其由兩點之間的距離決定，且滿足0<ω<1，∑ω=1：

式中，B是以i點為中心的相似區(qū)域，S為區(qū)域B中點的個數(shù)，λ為濾波器參數(shù)，v為時序值。

1.4 優(yōu)化的EWT-NLM算法

優(yōu)化的EWT-NLM算法的具體步驟如下：1)對原始信號x(i)進行EWT分解，獲取1+N個模態(tài)分量。2)計算各分量的SE值，其中，重構(gòu)維數(shù)m=2，閾值r=0.25σ(x)，σ為標(biāo)準(zhǔn)差。通過經(jīng)驗閾值判定，將經(jīng)驗尺度、小波分量劃分為有效信息分量和混合噪聲分量。3)對非信息分量(混合分量、噪聲分量)進行NLM濾波。經(jīng)過多次濾波實驗，得到合適的濾波經(jīng)驗參數(shù)為：搜索窗口t=50，鄰域窗口f=20，濾波參數(shù)h=250。4)將NLM濾波后的有效信號主導(dǎo)分量與信息分量重構(gòu)為最終的降噪信號。

選取信噪比(Rsn)與均方根誤差(RMSE)為降噪評價指標(biāo)，計算公式為：

式中，xi為原始時間序列信號，i為降噪后信號，L為信號長度。

2 實驗分析

2.1 模擬數(shù)據(jù)降噪分析

模擬信號采樣頻率為1 Hz，采樣點1 287個。不含噪聲的模擬數(shù)據(jù)見式(11)，各分量波形圖見圖1。采用白噪聲與有色噪聲組合模型，見式(12)。包含噪聲的模擬數(shù)據(jù)見式(13)，模擬噪聲信號及頻譜圖見圖2。

圖1 模擬信號各分量波形Fig.1 Waveform of each component of simulated signal

圖2 模擬噪聲信號及其頻譜Fig.2 Simulated noise signal and its spectrum

y1=8，y2=12ti,

(11)

式中，ti為坐標(biāo)歷元時刻，w(ti)為均值為0、信噪比為4的白噪聲，f1和f2均為0.4。則：

x=y1+y2+y3+y4+rt

(13)

式中，rt為組合噪聲。

將模擬噪聲信號進行EWT分解(圖3)。可以看出，原始信號被分解成1個經(jīng)驗尺度分量和11個經(jīng)驗小波分量(將f4～f11分量合并)。計算各分量SE值，結(jié)果見圖4。由圖3和4可見，分解所得的分量頻率逐步增大，各分量SE值也依次增大，表明分量的不規(guī)則度逐漸增大。選取經(jīng)驗閾值0.05為有效信號界限[7]?？梢园l(fā)現(xiàn)，前4個分量的SE值小于0.05，表明前4個分量為信號主導(dǎo)的模態(tài)分量，第5～10個分量為混合信號分量，第11～12個分量大于0.2，為噪聲主導(dǎo)的信號分量。

圖3 EWT分解模擬噪聲信號Fig.3 Simulated noise signal decomposed by EWT

圖4 各分量的SE值Fig.4 SE values of each component

將混合分量和噪聲分量重構(gòu)后進行NLM濾波，結(jié)果見圖5。可以看出，NLM能有效濾除較大噪聲，提取出混合噪聲分量中能量占比較小的真實信號，較好地展示原始信號中的振蕩周期變化。圖6為NLM濾波噪聲與模擬噪聲的對比?？梢园l(fā)現(xiàn)，模擬噪聲和濾波噪聲基本吻合。

圖5 NLM濾波結(jié)果Fig.5 NLM filtering results

圖6 NLM濾波噪聲與模擬噪聲對比Fig.6 Comparison of NLM filtering noise and simulated noise

采用自相關(guān)系數(shù)閾值法結(jié)合EMD、EWT和本文方法分別對模擬信號進行降噪，結(jié)果見表1、圖7。

表1 不同降噪方法評價參數(shù)統(tǒng)計Tab.1 Statistics of evaluation parameters for different noise reduction methods

圖7 不同降噪方法對比Fig.7 Comparison of different noise reduction methods

從表1和圖7可知，使用EMD方法降噪后的信號存在一定的端點效應(yīng)，從而影響降噪效果；EWT-NLM法可以有效提取出原始信號中的細節(jié)趨勢變化，更貼合原始信號，降噪的信號也更加平滑。從模擬實驗結(jié)果的評價指標(biāo)看，EWT-NLM法整體優(yōu)于EMD和EWT法，其模擬降噪結(jié)果的RMSE為2.876 4 mm、Rsn為17.67，相較于另外2種方法，RMSE分別降低10.63%和5.78%，Rsn分別提升22.54%和7.42%，從而驗證了本文方法的有效性。

2.2 GNSS實測高程時間序列降噪分析

使用EWT進行實測GNSS坐標(biāo)時間序列分解時，若分解層數(shù)過小會導(dǎo)致不能最大限度地分解信號的不同特征信息，若分解層數(shù)過大可能會造成信號之間的差異性被掩蓋。同時，EWT分解產(chǎn)生的所有分量中必包含有效信號和混合噪聲模態(tài)分量，且主要的有效信號大多集中于低頻信號中，那么分解的所有經(jīng)驗小波分量可以由低頻有效信號分量、中頻混合信號分量和高頻噪聲信號分量3個部分來概括，故分解層數(shù)與SE值選擇不當(dāng)會直接影響降噪效果。

選取我國12個陸態(tài)網(wǎng)基準(zhǔn)站原始高程坐標(biāo)時間序列進行降噪研究，數(shù)據(jù)來自中國地震局GNSS數(shù)據(jù)產(chǎn)品服務(wù)平臺，時間跨度為2010～2018年，采樣間隔為1/365.25 a，采樣頻率為365.25 Hz。為研究本文方法的自適應(yīng)性及抗差性，僅對原始高程坐標(biāo)時間序列進行粗差剔除和階躍修正，將剔除粗差和修正階躍后的坐標(biāo)時間序列作為原始坐標(biāo)序列。限于篇幅，主要闡述BJFS站的降噪過程，其余站點結(jié)果以圖表形式展示。

首先對BJFS站高程坐標(biāo)時間序列進行EWT分解，經(jīng)多次實驗，確定分解層數(shù)為14，分解信號見圖8?？梢钥闯?，趨勢項f1較好地符合原始時間序列信號的趨勢變化，f3分量為年周期信號，f5分量呈現(xiàn)0.5 a周期變化。再對各分量進行SE解算，熵重構(gòu)維數(shù)m取2，閾值r取0.25σ(X)(σ為計算序列標(biāo)準(zhǔn)差)，各分量SE值見表2和圖9。因GNSS高程坐標(biāo)時間序列具有豐富的組合噪聲特性，無法確定信號分量與噪聲分量的界限，故將所有分量劃分為低頻有效信號分量和中高頻混合噪聲分量，再對混合噪聲分量進行NLM濾波，最后重構(gòu)低頻有效信號分量與濾波信號為降噪信號，以達到最優(yōu)的降噪效果。

圖8 BJFS站分解信號結(jié)果Fig.8 The results of decomposing the signal at BJFS station

圖9 BJFS站各分量SE值Fig.9 SE values of each component at BJFS station

由表2可見，只有前4個分量的SE值小于0.05，故提取重構(gòu)前4個分量為有效信號分量，將其余分量重構(gòu)為混合噪聲分量，再對其進行NLM濾波，濾波效果見圖10。可以看出，分解出的有效信號整體較好地貼合原始時序，但不能呈現(xiàn)原始時序中較小的振蕩周期變化，在局部甚至出現(xiàn)較大的偏離。NLM有較好的濾波效果，可以有效剔除混合噪聲中頻率較大的噪聲，提取出混合噪聲信號中的微小振蕩變化。最后疊加濾波信號和有效信號為最終降噪信號。

表2 BJFS站各分量SE值統(tǒng)計Tab.2 Statistics of SE values of each component at BJFS station

圖10 有效信號分量與濾波效果Fig.10 Effective signal components and filtering effect

使用EMD、EWT和EWT-NLM方法進行對比，各方法降噪效果見圖11，評價參數(shù)統(tǒng)計見表3。

表3 不同降噪方法評價參數(shù)統(tǒng)計Tab.3 Statistics of evaluation parameters for different noise reduction methods

由圖11可見，本文方法和EWT方法相比于EMD方法可以較好地避免端點效應(yīng)造成的影響，在原始信號的首尾兩端有較好的降噪效果；經(jīng)過本文方法降噪的信號更貼合原始信號，可以有效反映局部趨勢運動變化及微小的周期振蕩。從表3可以看出，本文方法降噪效果最優(yōu)，相較于其他2種方法，其RMSE分別降低13.41%、7.13%，Rsn分別提升22.03%、9.72%。

圖11 不同降噪方法效果對比Fig.11 Comparison of different noise reduction methods

3 結(jié) 語

本文提出一種引入SE優(yōu)化的EWT結(jié)合NLM濾波的組合自適應(yīng)降噪方法。相較于EMD方法，該方法可以有效削弱端點效應(yīng)的影響，從而提升降噪效果；相較于EMD與EWT方法，該方法可呈現(xiàn)出原始信號中的局部運動趨勢變化及較小振幅的周期振蕩變化。在模擬數(shù)據(jù)和實測數(shù)據(jù)實驗中，該方法降噪效果均最優(yōu)。