自適應(yīng)噪聲整體經(jīng)驗(yàn)?zāi)B(tài)分解法在CORS站數(shù)據(jù)處理中的應(yīng)用

李 婷

(浙江省國(guó)土勘測(cè)規(guī)劃有限公司 浙江杭州 310030)

作為重要的空間數(shù)據(jù)基礎(chǔ)設(shè)施，連續(xù)運(yùn)行參考站(CORS)為社會(huì)各行業(yè)提供穩(wěn)定的位置信息服務(wù)。目前，國(guó)內(nèi)外許多國(guó)家和地區(qū)已經(jīng)開展或建成了CORS系統(tǒng)，建立全天候、全覆蓋、高精度和實(shí)時(shí)動(dòng)態(tài)定位的衛(wèi)星導(dǎo)航系統(tǒng)是國(guó)際大地測(cè)量的一個(gè)發(fā)展趨勢(shì)。對(duì)CORS站進(jìn)行觀測(cè)時(shí)，環(huán)境因素的影響使得其觀測(cè)數(shù)據(jù)在高程方向存在較大噪聲。因此，通過(guò)剔除觀測(cè)數(shù)據(jù)中的噪聲，獲取CORS站高程方向真實(shí)值具有重要的意義。

對(duì)于CORS站高程時(shí)間序列的非線性、非平穩(wěn)信號(hào)，通常利用事先給定的周期函數(shù)模型擬合趨勢(shì)項(xiàng)和周期項(xiàng)去除非平穩(wěn)性，但是，這種方法得到的分析結(jié)果不準(zhǔn)確，甚至?xí)x序列真實(shí)的運(yùn)動(dòng)情況。近些年來(lái)，整體經(jīng)驗(yàn)?zāi)B(tài)分解(Ensemble Empirical Mode Decomposition，EEMD)方法被廣泛應(yīng)用到時(shí)間序列分析中，經(jīng)EEMD法分解后的信號(hào)可以得到若干含有特征時(shí)間尺度的本征模態(tài)函數(shù)(Intrinsic Mode Function，IMF)分量，從而對(duì)信號(hào)進(jìn)行多分辨率分析[1]。吉長(zhǎng)東等[2]將EEMD方法應(yīng)用到CORS站時(shí)間序列分析中，有效地改善了信號(hào)分解過(guò)程中產(chǎn)生的模態(tài)分解問(wèn)題，但還會(huì)存在整體平均不能完全消除噪聲的影響。本文引入一種新的噪聲自適應(yīng)整體經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition with Adaptive Noise，EEMDAN)方法，與EEMD方法不同的是，EEMDAN方法在經(jīng)驗(yàn)?zāi)B(tài)分解處理時(shí)添加白噪聲分解的模態(tài)分量，可以改進(jìn)CORS站高程時(shí)間序列信號(hào)分解的模態(tài)混疊。最后利用某市高程CORS站點(diǎn)數(shù)據(jù)，驗(yàn)證EEMDAN方法處理CORS站高程時(shí)間序列的有效性。

1 EEMD及EEMDAN

1.1 EEMD

EEMD可以對(duì)待處理信號(hào)的極大值特征尺度進(jìn)行信號(hào)分解，通過(guò)不斷篩選，將待分解信號(hào)分解為頻率各不相同的IMF，根據(jù)頻率的不同來(lái)反映周期長(zhǎng)度，體現(xiàn)了對(duì)信號(hào)多分辨率分析的濾波過(guò)程。EEMD方法的濾波過(guò)程是通過(guò)對(duì)信號(hào)進(jìn)行EEMD分解產(chǎn)生IMF分量，根據(jù)IMF分量特征，構(gòu)造時(shí)空濾波器。

分解后的IMF可以反映信號(hào)不同范圍的特征尺度，通過(guò)特征可以對(duì)信號(hào)進(jìn)行濾波。以往是通過(guò)傅里葉變換在頻域上實(shí)現(xiàn)濾波，其缺點(diǎn)是會(huì)在頻域產(chǎn)生諧波信號(hào)，很難將信號(hào)的頻譜有效分離。

對(duì)監(jiān)測(cè)體進(jìn)行變形監(jiān)測(cè)采集的GNSS數(shù)據(jù)，噪聲與有用信號(hào)的時(shí)頻特性是有區(qū)別的。GNSS時(shí)間序列中通常存在由于多路徑影響產(chǎn)生的頻率低于結(jié)構(gòu)特性的低頻噪聲，也會(huì)存在頻率高于結(jié)構(gòu)特性的高頻噪聲。結(jié)合GNSS的多路徑特點(diǎn)，對(duì)GNSS監(jiān)測(cè)數(shù)據(jù)經(jīng)EEMD分解產(chǎn)生的不同頻率特征相關(guān)性進(jìn)行分析，從而消除多路徑誤差對(duì)觀測(cè)信號(hào)的影響。

1.2 EEMDAN

EEMDAN是對(duì)EEMD方法的改進(jìn)。與EEMD分解所有 IMF分量進(jìn)行整體平均不同的是，EEMDAN分解得到第一階分量后即進(jìn)行整體平均，并將第一階整體平均IMF剔除后得到除第一階整體平均IMF的剩余分量，再對(duì)剩余分量進(jìn)行整體平均，得到第二階、第三階直至最后一階IMF分量[3]。EEMDAN濾波過(guò)程的主要步驟如下[4]：

第一步：將白噪聲添加到信號(hào)中做經(jīng)驗(yàn)?zāi)B(tài)分解(Empirical Mode Decomposition，EMD)，對(duì)分解的IMF分量進(jìn)行平均化處理，獲取第一階IMF分量。

第二步：剔除第一階分量后，得到剩余信號(hào)，再分解后得到含有噪聲的模態(tài)分量集；在剩余信號(hào)中加入該模態(tài)分量集，組合成分解信號(hào)，繼續(xù)進(jìn)行EMD分解，得到第二階分量。

第三步：剔除第二階分量后，重復(fù)第二步，直至剔除最后一階模態(tài)分量后得到剩余信號(hào)。

EEMDAN利用白噪聲頻譜分布均勻這一特性，對(duì)CORS站點(diǎn)的時(shí)間序列進(jìn)行投影，以達(dá)到消除由分解帶來(lái)的模態(tài)混疊。EEMDAN可以將上一次分解后的剩余信號(hào)作為下一次分解信號(hào)，減小了重構(gòu)誤差和整體平均噪聲在沒(méi)有消除的情況下對(duì)信號(hào)分解造成的影響。

2 CORS站點(diǎn)時(shí)間序列分解評(píng)價(jià)標(biāo)準(zhǔn)

對(duì)CORS站點(diǎn)時(shí)間序列進(jìn)行分解，可以通過(guò)兩個(gè)指標(biāo)評(píng)價(jià)其分解效果，一是分解后的IMF分量是否存在模態(tài)混疊；二是參考正交性指標(biāo)(Orthogona-lity Index，IO)，分解精度隨著IO的減小而提高。IO可表示為

(1)

將CORS站時(shí)間序列分解為一系列不同的IMF分量，若將Tj和CMSEj定義為第j個(gè)IMF分量的平均周期與能量密度，則可以通過(guò)參考分量的能量密度與平均周期的乘積得到信號(hào)與噪聲的臨界位置(const值)，通過(guò)判別臨界位置來(lái)剔除包含噪聲的分量，完成對(duì)信號(hào)的重構(gòu)。通過(guò)計(jì)算重構(gòu)后的相關(guān)系數(shù)R和均方根誤差(RMSE)來(lái)量化重構(gòu)效果。相關(guān)系數(shù)R越大，重構(gòu)效果就越好；RMSE值越小，重構(gòu)信號(hào)的精度與質(zhì)量越高[5-7]。

R和RMSE可分別表示為

(2)

(3)

3 試驗(yàn)分析

本文選取某市CORS站BTRD站點(diǎn)2016年1月1日—2019年12月31日的GNSS高程觀測(cè)數(shù)據(jù)為試驗(yàn)數(shù)據(jù)，對(duì)高程觀測(cè)數(shù)據(jù)進(jìn)行均值化處理，圖1所示為高程觀測(cè)數(shù)據(jù)時(shí)間序列。

圖1 BTRD站點(diǎn)2016—2019年高程方向上時(shí)間序列

3.1 分解精度分析

為了便于對(duì)比，將原始信號(hào)標(biāo)準(zhǔn)差與噪聲的比值設(shè)置為10，經(jīng)EEMD和EEMDAN兩種方法分解得到部分低頻IMF分量如圖2所示(實(shí)線表示的是EEMD分解結(jié)果，虛線表示的是EEMDAN分解結(jié)果)。

圖2 BTRD站點(diǎn)高程時(shí)間序列分解得到的部分IMF分量

本文使用IMF分量的方差貢獻(xiàn)率獲取CORS站點(diǎn)時(shí)間序列運(yùn)動(dòng)的周期形式。IMF分量在整個(gè)時(shí)間序列中的比重可以用該分量的方差貢獻(xiàn)率反映，時(shí)間序列中的主要周期運(yùn)動(dòng)貢獻(xiàn)項(xiàng)由方差貢獻(xiàn)率高的IMF分量組成。BRTD站點(diǎn)高程時(shí)間序列經(jīng)EEMD和EEMDAN兩種方法分解的低頻分量的方差貢獻(xiàn)率如表1所示。

表1 不同IMF分量的方差貢獻(xiàn)率 單位:%Tab.1 Variance Contribution Rate in Different IMF Components 方法c4c5c6c7c8c9EEMD21.428.64.11.25.70.8EEMDAN24.826.33.93.10.91.4

由表1可知，BRTD站點(diǎn)高程時(shí)間序列經(jīng)EEMD分解后，IMFc4與IMFc5方差貢獻(xiàn)率所占比重最高，分別為21.4%與28.6%，為時(shí)間序列中的主要周期貢獻(xiàn)項(xiàng)。BRTD站點(diǎn)高程時(shí)間序列經(jīng)EEMDAN分解后，同樣是IMFc4與IMFc5的方差貢獻(xiàn)率占比最高，分別為24.8%與26.3%。

對(duì)于BRTD站點(diǎn)高程時(shí)間序列來(lái)說(shuō)，可以通過(guò)IMFc4與IMFc5兩個(gè)分量表現(xiàn)出周期形式。為了評(píng)價(jià)EEMD和EEMDAN兩種方法的分解精度，計(jì)算方差貢獻(xiàn)率較高的IMF分量之間的正交指數(shù)，結(jié)果如表2所示。

由表2可知，取200次整體平均的分解精度比取100次整體平均的分解精度低，也符合文獻(xiàn)[8]中的推薦值。通過(guò)計(jì)算EEMD和EEMDAN兩種方法分解后的IMF分量的IO值，可知EEMDAN方法的分解精度比EEMD方法的分解精度提高了53.7%。

表2 方差貢獻(xiàn)率較高的IMF分量之間的正交指數(shù)Tab.2 Orthogonal Index of IMF Components in Higher Variance Contribution Rate方法N=100N=200EEMD0.054 10.068 7EEMDAN0.035 20.041 6

3.2 重構(gòu)信號(hào)分析

BRTD站點(diǎn)高程時(shí)間序列經(jīng)EEMD和EEMDAN兩種分解方法分解后得到不同頻率的IMF分量，其中IMFc1～I(xiàn)MFc3是高頻分量。為了有效剔除噪聲并重構(gòu)得到有用的信號(hào)，通過(guò)IMF分量的能量密度與平均周期的乘積指標(biāo)來(lái)識(shí)別和剔除噪聲，其計(jì)算結(jié)果見表3。

表3 BRTD站點(diǎn)分解分量的能量密度與平均周期乘積Tab.3 Product of Energy Density and Average Period of Decomposition Component of BRTD Station方 法指 標(biāo)c1c2c3c4c5c6c7EEMDT/d3.4067.25814.32631.43360.221177.274384.214CMSE/(cm2/d)0.3570.2140.1550.1640.0980.2860.179T×CMSE/cm21.2161.5532.2205.1555.90250.70068.774EEMDANT/d3.3577.16615.82328.16344.581159.548181.532CMSE/(cm2/d)0.2440.1950.0970.1470.1260.3270.218T×CMSE/cm20.8191.3971.5354.1405.61752.17239.574

陳仁祥等[9]根據(jù)CMSE值是否發(fā)生明顯變化來(lái)判斷信號(hào)與噪聲的臨界點(diǎn)，但是，這種判別方式在進(jìn)行噪聲剔除時(shí)，會(huì)將雙月與近似月的分量剔除。故本文引入const值作為噪聲與信號(hào)的判別標(biāo)準(zhǔn)[10-12]。

對(duì)于BRTD站高程時(shí)間序列，經(jīng)EEMDAN處理后的前3個(gè)分量的const均值為1.250，前3個(gè)分量const值與均值的最大偏差不超過(guò)34.5%，而IMFc4的const值與均值的偏差超過(guò)2倍；經(jīng)EEMD處理后的前3個(gè)分量的const均值為1.663，前3個(gè)分量const值與均值的最大偏差不超過(guò)36.8%，最小偏差為7.1%，而IMFc4的const值與均值的偏差超過(guò)2倍。因此，EEMD與EEMDAN分解得到的結(jié)果相同。

const值對(duì)噪聲與信號(hào)的區(qū)別判別指標(biāo)主要是針對(duì)信號(hào)中高斯白噪聲，對(duì)有色噪聲的處理效果并不理想。通過(guò)噪聲識(shí)別對(duì)提取的信號(hào)進(jìn)行重構(gòu)，其結(jié)果如圖3所示，重構(gòu)信號(hào)的均方根誤差與相關(guān)系數(shù)見表4。

表4 信號(hào)重構(gòu)后的均方根誤差與相關(guān)系數(shù)Tab.4 RMSE and Correlation Coefficient after Signal Reconstruction方法RMSE/mm相關(guān)系數(shù)/%EEMD1.57596.7EEMDAN1.29698.8

由圖3與表4可知，經(jīng)EEMD與EEMDAN兩種方法重構(gòu)后的時(shí)間序列均能有效提取信號(hào)中的有用信息，但是經(jīng)EEMDAN處理后的信號(hào)均方根誤差比經(jīng)EEMD處理后的均方根誤差減少21.5%，相關(guān)系數(shù)增加2.1%，表明經(jīng)EEMDAN剔除噪聲，重構(gòu)信號(hào)的結(jié)果要比經(jīng)EEMD方法得到的結(jié)果更好。

圖3 兩種方法重構(gòu)后時(shí)間序列

這是因?yàn)镋EMDAN處理CORS站高程時(shí)間序列時(shí)，較好地利用了白噪聲頻譜均勻分布的特性，使信號(hào)能夠在遍布整個(gè)時(shí)頻空間且分布一致的白噪聲背景上，按照合適的尺度進(jìn)行投影，以此來(lái)削弱分解產(chǎn)生的模態(tài)混疊。

4 結(jié) 語(yǔ)

以EEMD為基礎(chǔ)的EENDAN方法改進(jìn)了對(duì)整體平均的處理，通過(guò)計(jì)算分量能量密度與平均周期的乘積、均方根誤差、相關(guān)系數(shù)和正交性指標(biāo)的值，驗(yàn)證了EEMDAN方法的分解精度與重構(gòu)后信號(hào)的效果都要比EEMD方法好。

基于信號(hào)特有的運(yùn)動(dòng)特征，EEMDAN方法可以將原始信號(hào)的有用信息提取出來(lái)，為CORS站點(diǎn)時(shí)間序列的處理提供有益的參考。但是本文對(duì)于如何剔除信號(hào)中有色噪聲沒(méi)有做進(jìn)一步的研究。