梁 沛 楊志強(qiáng) 楊 兵 田 鎮(zhèn) 陳 祥

1 長(zhǎng)安大學(xué)地質(zhì)工程與測(cè)繪學(xué)院，西安市雁塔路126號(hào)，710054

EMD[1]和CEEMD[2]在GNSS坐標(biāo)時(shí)間序列降噪過(guò)程中發(fā)揮著重要作用。但利用EMD、CEEMD進(jìn)行降噪時(shí)，需要確定噪聲和真實(shí)信號(hào)的界限，即需要確定分界IMF函數(shù)[3]。多數(shù)學(xué)者采用相關(guān)系數(shù)準(zhǔn)則[4]和連續(xù)均方誤差(consecutive mean squared error，CMSE)準(zhǔn)則[3]確定分界IMF分量。然而，相關(guān)系數(shù)準(zhǔn)則確定的模態(tài)分界點(diǎn)有時(shí)并不準(zhǔn)確，且無(wú)法直接給出分界IMF分量[5]；CMSE準(zhǔn)則在連續(xù)的IMF分量的均方誤差比較接近時(shí)難以判定信號(hào)與噪聲的分界點(diǎn)[6]。

GNSS坐標(biāo)時(shí)間序列會(huì)表現(xiàn)出明顯的非線性變化，特別是垂向的季節(jié)性變化更加明顯[7]。這種垂向坐標(biāo)時(shí)間序列的季節(jié)性變化主要由各類(lèi)地球物理效應(yīng)和系統(tǒng)誤差造成。對(duì)季節(jié)性信號(hào)的研究不僅有助于建立與維持動(dòng)態(tài)地球參考框架，同時(shí)對(duì)研究板塊構(gòu)造運(yùn)動(dòng)和地球動(dòng)力學(xué)過(guò)程也具有重要意義[8-9]。通過(guò)計(jì)算IMF分量的平均周期，找到符合季節(jié)性信號(hào)周期的IMF，即可將分解后代表季節(jié)信號(hào)和噪聲的IMF分量進(jìn)行分離。

因此，顧及到GNSS垂向坐標(biāo)時(shí)間序列中的季節(jié)信號(hào)，本文采用IMF分量的平均周期作為CEEMD降噪過(guò)程中篩選噪聲分量的準(zhǔn)則，并與以CMSE和相關(guān)系數(shù)為篩選準(zhǔn)則的CEEMD降噪方法進(jìn)行對(duì)比。

1 降噪原理與分界準(zhǔn)則

1.1 CEEMD降噪原理

k=2,…,m-1

(1)

式中，m為分解得到的IMF個(gè)數(shù)，cj為CEEMD分解得到的第j個(gè)IMF分量，rm為殘差分量。

1.2 平均周期準(zhǔn)則

LocalMax[cj(t)]i-1}

(2)

1.3 降噪效果評(píng)價(jià)指標(biāo)

引入RMS、半周年振幅、冪律噪聲、速度不確定度等4個(gè)指標(biāo)來(lái)評(píng)價(jià)平均周期準(zhǔn)則、CMSE準(zhǔn)則和相關(guān)系數(shù)準(zhǔn)則的CEEMD方法的降噪效果。

研究表明，白噪聲加冪律噪聲模型可作為分析GNSS垂向坐標(biāo)時(shí)間序列的最優(yōu)噪聲模型[10]。因此本文在白噪聲和冪律噪聲組合的噪聲模型下，利用極大似然估計(jì)方法[11]求解降噪前后GNSS垂向坐標(biāo)時(shí)間序列的模型參數(shù)，得到測(cè)站的半周年振幅、冪律噪聲和速度不確定度。

GNSS坐標(biāo)時(shí)間序列函數(shù)模型為[3,8]：

(3)

式中，y(t)為測(cè)站在t時(shí)刻下的位移，tR為參考時(shí)間，a為截距，v為測(cè)站運(yùn)動(dòng)速度，nF為頻率個(gè)數(shù)，通常用來(lái)擬合季節(jié)性變化，sk和φk分別為頻率ωk的振幅和相位，本文只考慮周年項(xiàng)和半周年項(xiàng)(k=1,2)，nJ為出現(xiàn)階躍的次數(shù)，bj為各種原因引起的階躍，H為Heaviside階躍函數(shù)，ε(t)為隨機(jī)過(guò)程。

2 實(shí)驗(yàn)數(shù)據(jù)與預(yù)處理

選取中國(guó)地殼運(yùn)動(dòng)觀測(cè)網(wǎng)絡(luò)227個(gè)(圖1)GNSS測(cè)站的垂向時(shí)間序列進(jìn)行降噪分析，其中，171個(gè)測(cè)站觀測(cè)時(shí)間為2010.33～2020.99年，25個(gè)測(cè)站觀測(cè)時(shí)間為2011.01～2020.99年，19個(gè)測(cè)站觀測(cè)時(shí)間為2012.00～2019.75年，4個(gè)測(cè)站觀測(cè)時(shí)間為2012.31～2020.99年，其余測(cè)站的觀測(cè)時(shí)間見(jiàn)表1。所有測(cè)站觀測(cè)數(shù)據(jù)已由GAMIT/GLOCK軟件進(jìn)行基線解算和約束平差處理。

圖1 GNSS站分布和觀測(cè)時(shí)長(zhǎng)Fig.1 Distribution and observation durationsof GNSS stations

表1 其余測(cè)站的觀測(cè)時(shí)間

上述坐標(biāo)時(shí)間序列中仍存在由環(huán)境、設(shè)備、數(shù)據(jù)處理策略等因素引起的粗差和數(shù)據(jù)缺失，在進(jìn)行降噪前，需要進(jìn)行數(shù)據(jù)預(yù)處理[11]。首先基于四分位距統(tǒng)計(jì)量識(shí)別和剔除粗差，然后采用分段三次Hermite插值法進(jìn)行插值[12]，補(bǔ)全序列中缺失的數(shù)據(jù)。

3 實(shí)驗(yàn)與分析

圖2 FJWY站降噪后時(shí)間序列Fig.2 Denoised time sequence at FJWY station

3.1 季節(jié)信號(hào)分析

GNSS垂向坐標(biāo)時(shí)間序列中的季節(jié)信號(hào)主要包括周年和半周年季節(jié)信號(hào)，反映地表水文變化、大氣等地球物理效應(yīng)對(duì)地表形變監(jiān)測(cè)的影響，因此降噪過(guò)程中剔除季節(jié)信號(hào)的現(xiàn)象，即被認(rèn)為是過(guò)度降噪。本文通過(guò)計(jì)算降噪前后GNSS垂向坐標(biāo)時(shí)間序列的半周年振幅，判斷3種方法是否會(huì)造成過(guò)度降噪，結(jié)果見(jiàn)圖3。由圖可見(jiàn)，CEEMD-T方法降噪后的序列與原序列的半周年振幅值較為符合，半周年振幅值損失較小，差值穩(wěn)定在一定范圍，表明CEEMD-T方法對(duì)原序列中的半周年信號(hào)保留較好，造成半周年振幅損失的部分原因是噪聲的扣除提高了其精度。而經(jīng)CEEMD-E和CEEMD-R方法降噪后，分別存在24個(gè)和79個(gè)測(cè)站的半周年振幅的改正率高于90%，表明這2種方法分別有10.57%和34.80%的測(cè)站出現(xiàn)過(guò)度降噪現(xiàn)象。

圖3 各測(cè)站降噪后時(shí)間序列Fig.3 The time sequence of each station after denoising

剔除過(guò)度降噪的站點(diǎn)，CEEMD-E方法降噪后的時(shí)間序列與原序列更為符合，甚至在許多站點(diǎn)上與原序列完全一致。因此，為進(jìn)一步對(duì)比3種方法的降噪效果，剔除過(guò)度降噪的100個(gè)測(cè)站，對(duì)剩余127個(gè)測(cè)站進(jìn)行分析。

3.2 RMS改正

降噪前后127個(gè)測(cè)站的RMS改正率分布見(jiàn)圖4。由圖可得，經(jīng)CEEMD-T降噪后的127個(gè)GNSS測(cè)站的平均RMS改正率為19.13%，比CEEMD-E、CEEMD-R的結(jié)果分別提高4.72%、3.36%。CEEMD-T、CEEMD-E、CEEMD-R三種方法降噪前后RMS改正率介于0%～20%的比例分別為49.61%、80.32%、64.57%，介于20%～40%的比例分別為49.60%、19.68%、35.43%。CEEMD-T方法的RMS改正率最高可以達(dá)到40.42%，而CEEMD-E、CEEMD-R方法的RMS改正率最大值分別為30.82%、27.49%。以上結(jié)果表明，CEEMD-T方法對(duì)優(yōu)化RMS的效果最佳。

圖4 127個(gè)GNSS測(cè)站的RMS改正率分布Fig.4 Distribution of RMS change rate of 127 GNSS stations

3.3 冪律噪聲與速度不確定度

降噪前后時(shí)間序列的冪律噪聲值見(jiàn)圖5。由圖可見(jiàn)，CEEMD-T方法對(duì)冪律噪聲的剔除效果優(yōu)于其他2種方法，能大幅度剔除所有站點(diǎn)序列中的冪律噪聲，降噪后冪律噪聲值均在5 mm/a0.21以下。CEEMD-R方法效果較差，降噪后大部分冪律噪聲值在1～10 mm/a0.21波動(dòng)，其中，QHTR站降噪后的冪律噪聲值高達(dá)13.09 mm/a0.21。CEEMD-E方法效果最差，大部分測(cè)站降噪后冪律噪聲值也在1～10 mm/a0.21波動(dòng)，其中，SCGU站甚至沒(méi)有起到剔除效果。CEEM-T、CEEMD-E、CEEMD-R三種方法降噪后冪律噪聲平均值由18.91 mm/a0.21分別降至2.21 mm/a0.21、4.27 mm/a0.21、3.55 mm/a0.21，其平均改正率分別為88.29%、77.73%、81.22%。以上結(jié)果表明，采用CEEMD-T方法能有效、穩(wěn)定地剔除GNSS垂向坐標(biāo)序列中的冪律噪聲，改正率較其他2種方法分別提高10.56%、7.07% 。

圖5 127個(gè)GNSS測(cè)站降噪前后時(shí)間序列的冪律噪聲值Fig.5 Power law noise values of the time sequence before and after denoising of 127 GNSS stations

圖6為127個(gè)GNSS測(cè)站降噪前后時(shí)間序列的速度不確定度。由圖可見(jiàn)，相比于其他2種方法，CEEMD-T方法能更好地改正所有站點(diǎn)序列的速度不確定度，降噪后速度不確定度均降低至0.2 mm/a以下，平均速度不確定度由0.60 mm/a降至0.08 mm/a，平均改正率為86.46% 。CEEMD-R方法除在QHTR站上速度不確定度為0.46 mm/a外，其余測(cè)站降噪后速度不確定度在0～0.3 mm/a范圍，平均改正率為78.13%。CEEMD-E方法對(duì)速度不確定度的改善效果最差，測(cè)站降噪后速度不確定度在0～0.3 mm/a波動(dòng)，但在大多數(shù)測(cè)站上的改正效果不如其他2種方法，且在SCGU站上同樣沒(méi)有起到改正效果，其速度不確定度平均改正率為73.72%。綜上，CEEMD-T方法改正率較其他2種方法分別提高8.33%、12.74%，可以穩(wěn)定、顯著地降低時(shí)間序列的速度不確定度。

圖6 127個(gè)GNSS測(cè)站降噪前后序列的速度不確定度Fig.6 The velocity uncertainty of the time sequence before and after denoising of 127 GNSS stations

通過(guò)上述分析可知，CEEMD-T方法降噪效果優(yōu)于CEEMD-E和CEEMD-R方法。CEEMD-E和CEEMD-R方法降噪不完全且不穩(wěn)定，獲得的降噪信號(hào)更接近帶噪聲的原始信號(hào)。這能在一定程度上解釋去除過(guò)度降噪站點(diǎn)后，CEEMD-E方法降噪后的半周年振幅與原序列更吻合、在部分測(cè)站上與原序列完全一致的現(xiàn)象。

3.4 CEEMD-T結(jié)果分析

對(duì)CEEMD-T方法得到的RMS、冪律噪聲、速度不確定度的改正率分布進(jìn)行分析，結(jié)果見(jiàn)圖7。由圖可見(jiàn)，在227個(gè)測(cè)站上，RMS改正率在0%～20%范圍的測(cè)站比例為65.49%，在20%～40%范圍的比例為34.07%，平均改正率為15.63%。87.61%的測(cè)站的冪律噪聲改正率介于80%～90%之間，同時(shí)其平均改正率達(dá)到87.41%。速度不確定度的平均改正率為85.57%，且超過(guò)95%的測(cè)站的速度不確定度改正率高于80%。綜上，CEEMD-T方法在所有測(cè)站均適用，對(duì)GNSS垂向坐標(biāo)時(shí)間序列的RMS、冪律噪聲、速度不確定度均有較好的改正效果。

圖7 227個(gè)GNSS測(cè)站降噪前后時(shí)間序列各指標(biāo)改正率分布Fig.7 The correction rate distribution of each indicator of time sequence of 227 GNSS stations before and after denoising

4 結(jié) 語(yǔ)

本文提出一種顧及GNSS坐標(biāo)時(shí)間序列中季節(jié)信號(hào)的CEEMD降噪方法CEEMD-T，并通過(guò)降噪前后的RMS值、冪律噪聲值、速度不確定度等指標(biāo)進(jìn)行降噪效果分析。結(jié)果表明，CEEMD-T方法能夠很好地保留GNSS測(cè)站垂向坐標(biāo)時(shí)間序列中的半周年振幅，不會(huì)出現(xiàn)過(guò)度降噪的現(xiàn)象，具有良好的穩(wěn)健性。