楊 兵,楊志強(qiáng),田 鎮(zhèn),陳 祥
長安大學(xué)地質(zhì)工程與測繪學(xué)院,陜西 西安 710054
長期積累的GNSS基準(zhǔn)站坐標(biāo)時(shí)間序列廣泛應(yīng)用于大地測量學(xué)、地球物理科學(xué)和環(huán)境學(xué)等領(lǐng)域,為高精度全球參考框架、地震等自然災(zāi)害的監(jiān)測與預(yù)警,以及地球深部動力學(xué)研究提供了重要的基礎(chǔ)數(shù)據(jù)[1-4]。然而,由于受到數(shù)據(jù)處理策略、模型誤差和接收機(jī)噪聲等GNSS技術(shù)誤差,以及包括環(huán)境負(fù)荷、測站熱膨脹效應(yīng)在內(nèi)的地球物理因素的影響,GNSS基準(zhǔn)站坐標(biāo)時(shí)間序列中除反映構(gòu)造運(yùn)動的趨勢項(xiàng)外,還包括了各種非線性信號(主要包括周年項(xiàng)、半周年項(xiàng))和噪聲,嚴(yán)重制約了其應(yīng)用[5-8]。因此,GNSS測站坐標(biāo)時(shí)間序列的降噪分析是目前的研究熱點(diǎn)。
經(jīng)驗(yàn)?zāi)B(tài)分解(empirical mode decomposition,EMD)作為一種非平穩(wěn)、非線性信號的處理方法[9],能夠較好地識別并剔除地震信號[10]、GNSS坐標(biāo)時(shí)間序列[11-13]和高精度測量儀器的數(shù)據(jù)[14-16]中的噪聲。EMD方法將GNSS坐標(biāo)時(shí)間序列自適應(yīng)地分解成一系列從高頻到低頻的本征模態(tài)函數(shù)(intrinsic mode function,IMF),再基于一定的篩選準(zhǔn)則識別噪聲和信號分量,從而達(dá)到降噪效果[9,16]。為了避免過度降噪和降噪不足,選取合適的篩選準(zhǔn)則至關(guān)重要。此外,EMD方法分解得到的IMF分量中可能存在尺度互異或相似的信號,即所謂的模態(tài)混疊效應(yīng),使得信號分量中混有部分噪聲,嚴(yán)重影響降噪效果[9,17-18]。
國內(nèi)外學(xué)者對EMD方法的篩選準(zhǔn)則進(jìn)行了大量研究。文獻(xiàn)[10]利用基于連續(xù)均方誤差(consecutive mean square error,CMSE)的EMD濾波方法[19]有效壓制了微震信號中的隨機(jī)噪聲。文獻(xiàn)[11,13]選擇能量密度E和平均周期T的乘積ET[20]作為分界指標(biāo),驗(yàn)證了該指標(biāo)識別GNSS坐標(biāo)時(shí)間序列噪聲的可行性。文獻(xiàn)[12]在各IMF分量與原始信號的相關(guān)系數(shù)[11,21]發(fā)生明顯變化處開始重構(gòu)信號,顯著提高了GNSS坐標(biāo)時(shí)間序列的信噪比。但是,文獻(xiàn)[12]的方法需要將各IMF分量與原始信號的相關(guān)系數(shù)和預(yù)先設(shè)置的閾值進(jìn)行比較,某些情況下信號的相關(guān)性過強(qiáng)或過弱,會造成噪聲分量的錯(cuò)誤識別。文獻(xiàn)[10—13]均是根據(jù)篩選準(zhǔn)則第一次顯著變化判斷分界分量,具有一定的主觀性,使降噪過程不穩(wěn)定。文獻(xiàn)[11—13]選用的測站僅包括我國區(qū)域內(nèi)的IGS站,不能較全面地反映EMD方法對我國GNSS測站的降噪效果。
針對EMD的模態(tài)混疊效應(yīng),文獻(xiàn)[17—18]研究了EMD分解過程的特點(diǎn),表明連續(xù)數(shù)據(jù)的EMD分解過程中也存在模態(tài)混疊效應(yīng),并分別提出了集合經(jīng)驗(yàn)?zāi)B(tài)分解和互補(bǔ)集合經(jīng)驗(yàn)?zāi)B(tài)分解算法以改善該現(xiàn)象。雖然這兩種改進(jìn)算法在一定程度上抑制了模態(tài)混疊效應(yīng),但需要預(yù)先加入噪聲,造成計(jì)算量增加,影響計(jì)算效率的同時(shí)也可能導(dǎo)致降噪不足[22]。
事實(shí)上,數(shù)據(jù)分布形狀的概率密度函數(shù)在一定程度上可以反映兩種信號的差異,那么度量概率相似性的Hausdorff距離(Hausdorff distance,HD)便可以界定噪聲分量與信號分量[16]。文獻(xiàn)[16]通過空氣動力學(xué)數(shù)據(jù)證明了HD方法在降噪方面的有效性和優(yōu)勢。因此,本文將該篩選準(zhǔn)則引入EMD方法對GNSS坐標(biāo)時(shí)間序列的降噪處理中。而另一方面,小波分解(wavelet decomposition,WD)作為一種時(shí)頻分析方法,能夠基于給定的小波基函數(shù)和分解層數(shù)對時(shí)域信號進(jìn)行多尺度分解,將得到的高頻分量和低頻分量的小波系數(shù)進(jìn)行閾值濾波、重構(gòu),從而實(shí)現(xiàn)含噪信號的降噪處理[23-25]。
基于以上分析,本文利用小波分解對產(chǎn)生模態(tài)混疊效應(yīng)的IMF分量進(jìn)行處理,探索其對EMD方法降噪和模態(tài)混疊效應(yīng)的改善效果,提出將Hausdorff距離與EMD相結(jié)合并聯(lián)合小波分解的降噪算法EMD-HD&WD,對我國大陸構(gòu)造環(huán)境監(jiān)測網(wǎng)絡(luò)(crustal movement observation network of China,CMONOC)的GNSS垂向坐標(biāo)時(shí)間序列進(jìn)行降噪處理。本文在引入新篩選準(zhǔn)則的同時(shí)改善模態(tài)混疊效應(yīng),以拓展EMD方法在GNSS數(shù)據(jù)處理中的應(yīng)用,并提出將GNSS測站的運(yùn)動模型的精度作為降噪效果對比的指標(biāo),以研究降噪方法對GNSS運(yùn)動模型的優(yōu)化效果。
HD通過計(jì)算兩個(gè)點(diǎn)集之間的概率密度函數(shù)來度量相似性。其他度量距離方法多以點(diǎn)集中兩個(gè)點(diǎn)的距離表述空間目標(biāo)的距離,而HD則兼顧目標(biāo)的整體波形和點(diǎn)集中的異常值[16]。其原理如下:
有限集合A={a1,a2,…,am}和B={b1,b2,…,bn}分別是GNSS測站的坐標(biāo)時(shí)間序列和EMD算法分解得到的IMF分量所對應(yīng)的概率密度函數(shù)值(probability density function,PDF)。集合A中元素a到有限集合B之間的距離定義為
(1)
那么,定義A和B之間的單向HD為
(2)
式中,max為最大值操作。
同理,B和A之間的單向HD可以定義為
(3)
單向HD具有非對稱性,這是極大極小函數(shù)的特性。絕對HD為2個(gè)單向HD的最大值,即
HD(A,B)=max{h(A,B),h(B,A)}
(4)
為了利用HD準(zhǔn)則篩選IMF分量,首先利用EMD算法將GNSS測站的坐標(biāo)時(shí)間序列x(t)分解為n個(gè)IMF分量,然后估計(jì)x(t)及各IMF分量的PDF,分別記為PDF(x)與PDF(IMF)={PDF(IMF(1)),PDF(IMF(2)),…,PDF(IMF(i)),…,PDF(IMF(n))},最后根據(jù)式(2)—式(4)計(jì)算x(t)與各IMF分量之間PDF的HD值,記為HD(i)
HD(i)=HD{PDF(x),PDF(IMF(i))}
(5)
式中,i=1,2,…,n。
HD的第1個(gè)全局極大值點(diǎn)所對應(yīng)的階數(shù)記作m,則第一階至第m階的IMF分量為噪聲分量,其余為信號分量,即
(6)
式中,arg max{·}為全局極大值操作。
小波分解具有良好的時(shí)頻特性和多尺度分解能力,可以對產(chǎn)生模態(tài)混疊效應(yīng)的IMF分量進(jìn)行多分辨率分析,剔除其中的噪聲,而HD根據(jù)GNSS坐標(biāo)時(shí)間序列與EMD分解得到的各IMF分量的概率密度函數(shù)進(jìn)行相似性度量,判斷噪聲與信號分量。本文結(jié)合小波分解和HD的特點(diǎn)和優(yōu)勢,提出了EMD-HD&WD算法。該算法一方面引入了HD篩選準(zhǔn)則,可以避免出現(xiàn)其他指標(biāo)的主觀性或陷入局部最小值的情況,減少某些情況下的誤判,另一方面聯(lián)合了小波分解,能夠在不需要先驗(yàn)信息的情況下改善模態(tài)混疊效應(yīng),提高EMD的降噪能力。
EMD-HD&WD算法的具體步驟如下:
(1) 依據(jù)EMD算法對GNSS坐標(biāo)時(shí)間序列進(jìn)行分解,得到n階IMF分量
(7)
式中,IMFi(t)代表從原始信號中分離出來的頻率從高到低的IMF分量;R(t)代表趨勢項(xiàng)分量。
(2) 基于核密度函數(shù)估計(jì)原始信號及各IMF分量的PDF。
(3) 利用HD比較IMF分量與原始信號的PDF的相似性,確定m。
(4) 將m階之前的各IMF分量重構(gòu),得到初始噪聲N0(t)
(8)
(5) 選取coif5小波基和8階分解層數(shù)[25],對m之后的IMF分量(不包括趨勢項(xiàng)分量)進(jìn)行小波分解降噪,得到第j階IMF分量的降噪后信號S1j(t)和噪聲N1j(t)
IMFj(t)=S1j(t)+N1j(t)
(j=m+1,m+2,…,n)
(9)
(6) 將初始噪聲與小波分解得到的噪聲重構(gòu)得到最終噪聲N1(t),最終信號S(t)則為S1(t)與趨勢分量之和
(10)
1.3.1 復(fù)合指標(biāo)T值
傳統(tǒng)降噪方法的質(zhì)量評價(jià)指標(biāo)主要包括均方根誤差(RMSE)、信噪比(SNR)、相關(guān)系數(shù)(ρ)及平滑度(r)。試驗(yàn)表明,單純依靠某一種或兩種指標(biāo),可能會出現(xiàn)不一致的現(xiàn)象[26]。文獻(xiàn)[26]提出通過變異系數(shù)定權(quán)將歸一化后的均方根誤差與平滑度線性組合的復(fù)合指標(biāo),以評價(jià)小波降噪的效果。本文在文獻(xiàn)[26]的基礎(chǔ)上,采用變異系數(shù)定權(quán)方法將上述4種傳統(tǒng)評價(jià)指標(biāo)進(jìn)行融合。4種傳統(tǒng)評價(jià)指標(biāo)的計(jì)算方法可以參考文獻(xiàn)[27],本文重點(diǎn)介紹復(fù)合指標(biāo)構(gòu)建過程中各指標(biāo)的歸一化和定權(quán)方法。
各指標(biāo)的歸一化處理如下
(11)
式中,RMSE、SNR、ρ和r分別表示均方根誤差、信噪比、相關(guān)系數(shù)及平滑度;PRMSE(i)、PSNR(i)、Pρ(i)和Pr(i)分別為測站i的各指標(biāo)經(jīng)歸一化處理的結(jié)果;max(·)和min(·)分別為取最大值和最小值操作。
各指標(biāo)的變異系數(shù)定權(quán)公式為
(12)
式中,CV=σ/u表示各傳統(tǒng)指標(biāo)的變異系數(shù),σ和u分別表示各歸一化指標(biāo)的方差和均值;W表示各歸一化指標(biāo)在復(fù)合指標(biāo)中的權(quán)重。則復(fù)合指標(biāo)T可表示為
T=PRMSE×WPRMSE+PSNR×WPSNR+Pρ×WPρ+Pr×WPr
(13)
根據(jù)4個(gè)傳統(tǒng)指標(biāo)在數(shù)據(jù)降噪過程中的意義與特征以及歸一化和變異系數(shù)定權(quán)的計(jì)算過程可知,在對比不同降噪方法時(shí),復(fù)合指標(biāo)T越小,表示降噪效果越好。
1.3.2 GNSS坐標(biāo)時(shí)間序列函數(shù)模型精度
GNSS坐標(biāo)時(shí)間序列的函數(shù)模型[28-29]可以表示為
(14)
式中,y(t)是測站在時(shí)刻t下的位移;tR是參考時(shí)間;a為截距;v是測站運(yùn)動速度;nF是頻率個(gè)數(shù),通常用來擬合季節(jié)性變化;sk和φk分別是頻率ωk的振幅和相位,本文只考慮周年項(xiàng)(k=1)和半周年項(xiàng)(k=2);nJ是出現(xiàn)階躍的次數(shù);bj是階躍時(shí)刻tj時(shí)測站在地心笛卡兒坐標(biāo)中的瞬時(shí)位移向量;H是Heaviside階躍函數(shù);ε(t)為隨機(jī)過程。
研究表明,白噪聲加閃爍噪聲的噪聲組合模型[6,24]是中國區(qū)域測站垂向坐標(biāo)時(shí)間序列的最優(yōu)噪聲模型。此外,為了直觀地比較降噪前后GNSS時(shí)間序列的噪聲振幅變化和降噪算法的效果,降噪前后的噪聲模型保持一致[6,24]。因此,本文在白噪聲加閃爍噪聲的模型下利用最大似然估計(jì)法[29]解算降噪前后GNSS垂向坐標(biāo)時(shí)間序列的模型參數(shù),得到測站的速度不確定度和閃爍噪聲振幅。通過對比降噪前后速度不確定度和閃爍噪聲振幅判斷模型精度,屬于內(nèi)符合精度。速度不確定度和閃爍噪聲振幅越小,表明降噪效果越好[7,28-29]。GNSS垂向坐標(biāo)時(shí)間序列函數(shù)模型的計(jì)算過程及相關(guān)誤差分析的公式參見文獻(xiàn)[29]的式(7)、式(15—31)及附錄。
以GNSS為主的CMONOC為我國大陸地殼運(yùn)動監(jiān)測提供了高精度、高時(shí)空分辨率的基礎(chǔ)數(shù)據(jù)。考慮到觀測時(shí)長和數(shù)據(jù)缺失情況,本文選取了CMONOC中1999.16—2020.22年的149個(gè)GNSS測站坐標(biāo)時(shí)間序列,其中1999.16—2020.22年共10個(gè)測站;2010.32—2020.22年共104個(gè)測站;2012.19—2019.75年共16個(gè)測站;其余測站的觀測時(shí)間跨度詳見表1。本文采用測站的空間分布和觀測時(shí)長的詳細(xì)信息如圖1和表1所示。本文所選取的GNSS測站的垂向坐標(biāo)時(shí)間序列是由中國地震局GNSS數(shù)據(jù)產(chǎn)品服務(wù)平臺(http:∥www.cgps.ac.cn)所提供。GNSS觀測數(shù)據(jù)經(jīng)GAMIT基線解算處理、GLOBK約束平差,得到GNSS坐標(biāo)時(shí)間序列[4,30]。
表1 部分GNSS基準(zhǔn)站的觀測時(shí)間Tab.1 Observation times of GNSS reference stations
圖1 選取的GNSS基準(zhǔn)站的分布和觀測時(shí)長Fig.1 Locations and observation durations of GNSS continuous stations used in this study
由于GNSS測站設(shè)備更換、觀測條件和數(shù)據(jù)處理策略等因素的影響,致使GNSS坐標(biāo)時(shí)間序列存在粗差和數(shù)據(jù)缺失情況。為了后續(xù)測站運(yùn)動特征分析,在進(jìn)行降噪前,需要對原始GNSS垂向坐標(biāo)時(shí)間序列進(jìn)行預(yù)處理[29,31]:首先基于四分位距統(tǒng)計(jì)量識別和剔除粗差,即當(dāng)GNSS垂向坐標(biāo)時(shí)間序列中的數(shù)值大于3倍的GNSS坐標(biāo)時(shí)間序列的四分位距統(tǒng)計(jì)量時(shí),可認(rèn)為是粗差,予以剔除;然后采用分段三次Hermite插值法對GNSS坐標(biāo)時(shí)間序列中缺失的數(shù)據(jù)進(jìn)行補(bǔ)全。
本文對CMONOC的149個(gè)測站垂向坐標(biāo)時(shí)間序列采用EMD-HD&WD算法和文獻(xiàn)[12—14,16]的方法進(jìn)行降噪處理,計(jì)算復(fù)合指標(biāo)T值、降噪前后GNSS垂向坐標(biāo)時(shí)間序列的運(yùn)動參數(shù)及閃爍噪聲振幅。為方便表述,文獻(xiàn)[12—14,16]中以CMSE、ET、相關(guān)系數(shù)和HD為篩選準(zhǔn)則的EMD降噪方法及處理結(jié)果分別記為EMD-CMSE、EMD-ET、EMD-CORR、EMD-HD,降噪前和本文算法的結(jié)果記為raw和EMD-HD&WD。
GNSS坐標(biāo)時(shí)間序列降噪的主要目的是獲取更真實(shí)的構(gòu)造與季節(jié)性運(yùn)動信號(趨勢項(xiàng)和周期項(xiàng)),因此本文在對比各降噪方法的效果之前,首先分析了降噪后的構(gòu)造及季節(jié)性信號是否受到影響。從信號降噪的原理出發(fā),數(shù)據(jù)處理過程并未涉及測站的趨勢項(xiàng),因而本文的方法不會影響穩(wěn)態(tài)的構(gòu)造運(yùn)動信號。周年項(xiàng)和半周年項(xiàng)的振幅是衡量季節(jié)性運(yùn)動的重要指標(biāo)[4,28]。為進(jìn)一步分析降噪過程是否會對季節(jié)性運(yùn)動信號產(chǎn)生影響,本文計(jì)算了降噪前后149個(gè)GNSS時(shí)間序列的周年項(xiàng)和半周年項(xiàng)的平均振幅(表2)。EMD-CMSE、EMD-CORR、EMD-ET、EMD-HD和EMD-HD&WD所得到的周年、半周年項(xiàng)的振幅均略小于降噪前的結(jié)果,但相差不大,因此本文所采用的降噪方法并沒有平滑掉測站的地球物理信號。圖2為不同降噪方法對HEZJ測站的降噪結(jié)果。
圖2 不同降噪方法得到的噪聲與信號(以HEZJ站為例)Fig.2 The time series of noise and signal by different methods for HEZJ station
表2 GNSS坐標(biāo)時(shí)間序列降噪前后GNSS時(shí)間序列的周年項(xiàng)和半周年項(xiàng)的平均振幅Tab.2 The average amplitudes of annual and semi-annual signal on GNSS coordinate time series with different strategies mm
2.2.1 最優(yōu)篩選準(zhǔn)則的確定
圖3是不同降噪方法所確定的m和T值。本文根據(jù)圖3(a)計(jì)算了不同篩選準(zhǔn)則得到的m所占總階數(shù)的比重。結(jié)果表明,EMD-HD得到的比重平均值為0.52,大于0.48(EMD-CMSE)、0.44(EMD-CORR)、0.43(EMD-ET),說明文獻(xiàn)[13]的降噪結(jié)果與文獻(xiàn)[14]相近,EMD-HD所確定的m較其他準(zhǔn)則整體偏大。表3為圖3(b)中各種降噪方法所得到的T值的平均值。從表3可知,EMD-CMSE、EMD-CORR、EMD-HD和EMD-ET的149個(gè)T值的平均值分別是0.61、0.85、0.52、0.59。結(jié)合圖3(a)和表3可以發(fā)現(xiàn),EMD-CORR與EMD-ET的降噪結(jié)果相近,但是T值卻相差0.24,原因是因?yàn)榍罢咝枰A(yù)先設(shè)置閾值,但該閾值僅對部分測站比較有效[12];EMD-ET確定為噪聲的IMFs分量整體小于EMD-CMSE,但兩者的T值卻接近,這是由于比較的T值是149個(gè)測站數(shù)據(jù)的平均值,而EMD-ET對于部分測站具有較好的作用[11,13],改善了該方法的T值平均值;從m的比重來看,EMD-CMSE、EMD-CORR和EMD-ET所識別的噪聲的IMFs分量小于EMD-HD的結(jié)果,與T值的結(jié)果一致,說明前3種方法降噪不足。以上分析表明:EMD-CORR和EMD-ET僅適用于部分測站;HD優(yōu)于其他篩選準(zhǔn)則。
表3 基于不同準(zhǔn)則得到的T值平均值Tab.3 Averages of T obtained by different criteria
圖3 不同降噪方法所確定的m和T值Fig.3 m and T obtained by different methods
文獻(xiàn)[23,25]表明,白噪聲可以被各類方法很好地識別并剔除,但對閃爍噪聲的剔除效果卻不一致。圖4是降噪前后GNSS垂向坐標(biāo)時(shí)間序列的閃爍噪聲及相對于降噪前的增益率,并給出149個(gè)GNSS測站時(shí)序在降噪前后的閃爍噪聲的平均值,見表4。結(jié)果表明,與降噪前的GNSS坐標(biāo)時(shí)間序列相比,經(jīng)不同方法降噪處理后的GNSS坐標(biāo)時(shí)間序列的閃爍噪聲振幅得到顯著改善,均小于5 mm·a-0.25。其中,EMD-HD的平均改正率為79.7%,增益率大于85%的測站數(shù)占比為54.4%。表4的結(jié)果與圖3(a)相一致,表明HD的增益效果均高于其他準(zhǔn)則的結(jié)果。
圖4 降噪前后GNSS坐標(biāo)時(shí)間序列的閃爍噪聲振幅及相對于降噪前的增益率Fig.4 Amplitudes of flicker noise of GNSS coordinate time series before and after noise reduction and its improvement rates gain rate relative to before noise reduction
圖5給出了降噪前后GNSS垂向坐標(biāo)時(shí)間序列的速度不確定度及相對于降噪前的增益率,以進(jìn)一步分析本文算法對GNSS測站的運(yùn)動模型的優(yōu)化效果。表4是降噪前后149個(gè)GNSS測站的速度不確定度的平均值。由表4可知,速度不確定度從降噪前的0.54 mm·a-1減少到0.200 mm·a-1以下。由圖5可以得出,EMD-HD的平均改正率為79.7%,增益效果均高于其他準(zhǔn)則。通過對比圖5和圖4的結(jié)果可知,不同降噪方法對閃爍噪聲振幅與速度不確定度的改善效果表現(xiàn)出較強(qiáng)的一致性,說明降噪處理不僅可以剔除GNSS坐標(biāo)時(shí)間序列中的噪聲,而且可以提高GNSS測站運(yùn)動速率的估值精度,兩者具有統(tǒng)一性。
表4 降噪前后的GNSS坐標(biāo)時(shí)間序列的閃爍噪聲振幅和速度不確定度的平均值Tab.4 Averages of flicker noise amplitude and velocity uncertainty of GNSS coordinate time series before and after noise reduction
圖5 降噪前后GNSS垂向坐標(biāo)時(shí)間序列的速度不確定度及相對于降噪前的增益率Fig.5 Amplitudes of velocity uncertainty of GNSS coordinate time series before and after noise reduction and its improvement rates gain rate relative to before noise reduction
從以上分析可以看出,T值與噪聲振幅和速度不確定度的結(jié)果雖然并非完全一致,但均表明,HD得到的結(jié)果好于其他準(zhǔn)則的結(jié)果,說明HD能夠更準(zhǔn)確地確定噪聲分量。
2.2.2 模態(tài)混疊效應(yīng)的改善
圖6是EMD-HD&WD和EMD-HD的閃爍噪聲振幅和速度不確定度估值及EMD-HD&WD相對于EMD-HD的增益率。從表4和圖6可知,EMD-HD&WD方法降噪后的速度不確定度和閃爍噪聲振幅的平均值分別是0.065 mm·a-1和1.710 mm·a-0.25,相對于降噪前平均改正率均為88.4%;相對于EMD-HD方法,速度不確定度方面,96%的測站得到增益,平均增益率為40.4%;僅有6個(gè)測站有平均12.8%的負(fù)增益。閃爍噪聲振幅方面,96.6%的測站得到增益,平均增益率為40.1%;僅有5個(gè)測站有平均14.8%的負(fù)增益。結(jié)果表明,EMD-HD&WD算法改善模態(tài)混疊效應(yīng)的同時(shí)進(jìn)一步提高了GNSS坐標(biāo)時(shí)間序列的精度。
圖6 EMD-HD&WD和EMD-HD的運(yùn)動模型參數(shù)及EMD-HD&WD相對于EMD-HD的增益率Fig.6 The motion model parameters of EMD-HD&WD and EMD-HD and the gain rate of EMD-HD&WD compared to EMD-HD
本文針對基于EMD的GNSS基準(zhǔn)站坐標(biāo)時(shí)間序列的噪聲識別與剔除過程中存在的問題,提出了EMD-HD&WD算法,并通過對CMONOC的149個(gè)GNSS垂向坐標(biāo)時(shí)間序列的降噪處理驗(yàn)證了該算法的有效性及效果,得出以下結(jié)論:
(1) HD優(yōu)于現(xiàn)有的其他篩選準(zhǔn)則。相對于CMSE等準(zhǔn)則,EMD-HD和EMD-HD&WD更好地識別和剔除了GNSS垂向坐標(biāo)時(shí)間序列中的噪聲,并且降噪處理后測站的速度不確定度和噪聲振幅顯著小于其他方法,證實(shí)了HD作為篩選準(zhǔn)則的可靠性和普適性。
(2) EMD-HD&WD方法能較好地改善EMD的模態(tài)混疊效應(yīng)。EMD-HD&WD方法有效地將低頻IMF分量中混入的高頻信號分離出來,進(jìn)一步降低了測站的速度不確定度和閃爍噪聲,在保留坐標(biāo)時(shí)間序列中季節(jié)變化特征的同時(shí)具有較好的光滑性,達(dá)到更好的降噪效果,一定程度改善了EMD方法的模態(tài)混疊效應(yīng)。
(3) 結(jié)合GNSS坐標(biāo)時(shí)間序列分析的理論,本文提出將速度不確定度和噪聲振幅作為降噪評價(jià)指標(biāo),效果良好。經(jīng)EMD-HD&WD算法對GNSS坐標(biāo)時(shí)間序列進(jìn)行處理后,噪聲振幅和速度不確定度的平均值從降噪前的14.35 mm·a-0.25和0.54 mm·a-1降低至1.71 mm·a-0.25和0.065 mm·a-1,噪聲振幅和速度不確定度的平均改正率均為88.4%。本文將傳統(tǒng)的降噪方法的評價(jià)指標(biāo)和GNSS測站運(yùn)動參數(shù)進(jìn)行綜合評價(jià)降噪效果,為降噪方法的可靠性評價(jià)提供新的依據(jù)。
(4) 本文所提的降噪算法EMD-HD&WD能夠有效識別并剔除GNSS坐標(biāo)時(shí)間序列的噪聲,提高GNSS坐標(biāo)時(shí)間序列和垂向速度場的精度,不但有助于提取GNSS坐標(biāo)時(shí)間序列中季節(jié)性變化并解釋其物理機(jī)制,而且還為基于GNSS技術(shù)的地球參考框架的建立與精化,以及板塊構(gòu)造運(yùn)動、地球質(zhì)量變遷等地球動力學(xué)過程等研究提供了高精度數(shù)據(jù)。