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

1 長安大學地質(zhì)工程與測繪學院，西安市雁塔路126號，710054

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

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

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

1 降噪原理與分界準則

1.1 CEEMD降噪原理

k=2,…,m-1

(1)

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

1.2 平均周期準則

LocalMax[cj(t)]i-1}

(2)

1.3 降噪效果評價指標

引入RMS、半周年振幅、冪律噪聲、速度不確定度等4個指標來評價平均周期準則、CMSE準則和相關(guān)系數(shù)準則的CEEMD方法的降噪效果。

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

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

(3)

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

2 實驗數(shù)據(jù)與預處理

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

圖1 GNSS站分布和觀測時長Fig.1 Distribution and observation durationsof GNSS stations

表1 其余測站的觀測時間

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

3 實驗與分析

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

3.1 季節(jié)信號分析

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

圖3 各測站降噪后時間序列Fig.3 The time sequence of each station after denoising

剔除過度降噪的站點，CEEMD-E方法降噪后的時間序列與原序列更為符合，甚至在許多站點上與原序列完全一致。因此，為進一步對比3種方法的降噪效果，剔除過度降噪的100個測站，對剩余127個測站進行分析。

3.2 RMS改正

降噪前后127個測站的RMS改正率分布見圖4。由圖可得，經(jīng)CEEMD-T降噪后的127個GNSS測站的平均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改正率最高可以達到40.42%，而CEEMD-E、CEEMD-R方法的RMS改正率最大值分別為30.82%、27.49%。以上結(jié)果表明，CEEMD-T方法對優(yōu)化RMS的效果最佳。

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

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

降噪前后時間序列的冪律噪聲值見圖5。由圖可見，CEEMD-T方法對冪律噪聲的剔除效果優(yōu)于其他2種方法，能大幅度剔除所有站點序列中的冪律噪聲，降噪后冪律噪聲值均在5 mm/a0.21以下。CEEMD-R方法效果較差，降噪后大部分冪律噪聲值在1～10 mm/a0.21波動，其中，QHTR站降噪后的冪律噪聲值高達13.09 mm/a0.21。CEEMD-E方法效果最差，大部分測站降噪后冪律噪聲值也在1～10 mm/a0.21波動，其中，SCGU站甚至沒有起到剔除效果。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垂向坐標序列中的冪律噪聲，改正率較其他2種方法分別提高10.56%、7.07% 。

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

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

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

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

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

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

圖7 227個GNSS測站降噪前后時間序列各指標改正率分布Fig.7 The correction rate distribution of each indicator of time sequence of 227 GNSS stations before and after denoising

4 結(jié) 語

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