顧及鯨魚(yú)變分模態(tài)奇異譜分析的GNSS時(shí)序降噪軟件設(shè)計(jì)與實(shí)現(xiàn)

孫喜文, 魯鐵定, 賀小星, 陳紅康, 侯增楠

(1.東華理工大學(xué) 測(cè)繪與空間信息工程學(xué)院,江西 南昌 330013; 2.東華理工大學(xué) 自然資源部環(huán)鄱陽(yáng)湖區(qū)域礦山環(huán)境監(jiān)測(cè)與治理重點(diǎn)實(shí)驗(yàn)室,江西 南昌 330013;3.江西理工大學(xué) 土木與測(cè)繪工程學(xué)院,江西 贛州 341000;4.華東交通大學(xué) 交通工程學(xué)院,江西 南昌 330013)

近幾十年對(duì)全球衛(wèi)星導(dǎo)航系統(tǒng)(GNSS)連續(xù)運(yùn)行參考站(CORS)坐標(biāo)時(shí)間序列的研究為大地測(cè)量與地球動(dòng)力學(xué)提供了豐富的數(shù)據(jù)(姜衛(wèi)平等,2018;賀小星,2013)。隨著對(duì)大地測(cè)量成果所要求的精度越來(lái)越高,精確穩(wěn)定的時(shí)序參數(shù)及其非線性時(shí)變受到更多關(guān)注。從GNSS原始坐標(biāo)時(shí)間序列中提取地球物理相關(guān)信號(hào),成為大地測(cè)量及地球動(dòng)力學(xué)、地殼形變等研究的熱點(diǎn)之一,但GNSS站坐標(biāo)時(shí)間序列成分極為復(fù)雜,含有多種噪聲項(xiàng)(He et al.,2020; 梁沛等,2022)。郭翔(2016)提出對(duì)GNSS坐標(biāo)時(shí)間序列中季節(jié)信號(hào)進(jìn)行降噪,表示經(jīng)驗(yàn)?zāi)B(tài)分解方法(EMD)與互補(bǔ)集合經(jīng)驗(yàn)?zāi)B(tài)分解方法(CEEMD)在GNSS坐標(biāo)時(shí)間序列在降噪過(guò)程中具有重要的作用。EMD是一種針對(duì)非平穩(wěn)非線性信號(hào)的自適應(yīng)信號(hào)分析方法,然而胡愛(ài)軍等(2011)提出EMD方法的理論機(jī)理存在端點(diǎn)效應(yīng)和模態(tài)混疊等問(wèn)題;Dragomiretskiy等(2014)提出了變分模態(tài)分解(VMD)算法,但該方法設(shè)置的模態(tài)分解數(shù)值過(guò)小,會(huì)出現(xiàn)分解不充分的現(xiàn)象;Vautard等(1992)利用一種提取長(zhǎng)時(shí)間序列噪聲信號(hào)等的奇異譜分析(SSA)方法,分解重構(gòu)時(shí)間序列不同成分的信號(hào),如長(zhǎng)期趨勢(shì)信號(hào)、周期信號(hào)、噪聲信號(hào)等,從而對(duì)時(shí)間序列的結(jié)構(gòu)進(jìn)行分析。隨著對(duì)GNSS長(zhǎng)時(shí)間序列的研究,鑒于計(jì)算數(shù)據(jù)量龐大及分析過(guò)程煩瑣,熊常亮等(2021)實(shí)現(xiàn)了三種改進(jìn)的經(jīng)驗(yàn)?zāi)B(tài)分解降噪算法,但未顧及時(shí)間序列中的噪聲模型特性。賀小星等(2020)考慮降噪過(guò)程的復(fù)雜性,開(kāi)發(fā)了GNSS時(shí)間序列降噪軟件,但其僅是基于EMD進(jìn)行改進(jìn)的GNSS時(shí)間序列降噪。綜上所述,設(shè)計(jì)實(shí)現(xiàn)一種既能確定模態(tài)分解數(shù)與懲罰因子的方法,又能與多種降噪方法進(jìn)行對(duì)比的分析軟件,在北斗/GNSS坐標(biāo)序列非線性變化精密建模具有重要意義?；诖?筆者顧及長(zhǎng)時(shí)間序列中存在的噪聲模型特性及融合多種GNSS時(shí)間序列降噪方法,提出鯨魚(yú)變分模態(tài)奇異譜分析(WOA-VMD-SSA)新方法,這是一種融合自適應(yīng)噪聲完備經(jīng)驗(yàn)?zāi)B(tài)分解方法(CEEMDAN)(孫喜文等,2023)、EMD、EEMD(集合經(jīng)驗(yàn)?zāi)B(tài)分解)的降噪方法,并利用MATLAB圖形用戶界面設(shè)計(jì)了顧及WOA-VMD-SSA的GNSS時(shí)序降噪軟件,軟件界面清晰,降噪算法較為全面,無(wú)需外插接口及工具箱,可為GNSS長(zhǎng)時(shí)間序列降噪分析提供了有效參考。

1 WOA-VMD-SSA降噪新方法

針對(duì)變分模態(tài)分解中不適合的模態(tài)分解數(shù)和懲罰因子出現(xiàn)過(guò)度分解與分解不足等現(xiàn)象(陶國(guó)強(qiáng),2021),筆者提出一種鯨魚(yú)變分模態(tài)奇異譜分析(WOA-VMD-SSA)新方法。SSA是一種將信號(hào)分解出多個(gè)有序的子序列,通過(guò)SSA對(duì)GNSS時(shí)間序列進(jìn)行分解,從而得到其中所包含的周期信息的高頻時(shí)間序列、趨勢(shì)信息的低頻時(shí)間序列和代表噪聲的時(shí)間序列。

將一段長(zhǎng)度為n的時(shí)間序列S={s1,s2,…,sn}轉(zhuǎn)化為m×k維的軌跡矩陣X,其中,k=n-m+1,矩陣X為(徐洪學(xué)等,2021; 周偉輝等,2021):

(1)

計(jì)算XXT并對(duì)其進(jìn)行奇異值分解獲得m個(gè)特征值,從大至小排列為λ1,λ2,…,λm,對(duì)應(yīng)的特征向量為u1,u2,…,um,代入式(1)可得:

X=E1+E2+…Ed

(2)

將式(2)中的Ei劃分成p個(gè)不相交的子集I1,I2,…,Ip,設(shè)I={i1,i2,…,im},則合成矩陣Xi=Xi1,Xi2,…,Xim,計(jì)算合集I={i1,i2,…,im}的每個(gè)合成矩陣,可得式(3)為:

Xi=XI1,XI2,…,XIp

(3)

(4)

由此可將原始序列S={s1,s2,…,sn}被分解重構(gòu)為:

(5)

當(dāng)式(5)中v=1,u=L時(shí),所得重構(gòu)序列即是原始序列,選取不同的v和u進(jìn)行重構(gòu),即可得到包含周期信息的高頻序列和包含趨勢(shì)信息的低頻序列。

將時(shí)間序列通過(guò)迭代的方式進(jìn)行VMD分解,以求解最優(yōu)中心頻率與帶寬,之后對(duì)信號(hào)重構(gòu)并求解約束變分問(wèn)題,其數(shù)學(xué)模型(Dragomiretskiy et al., 2014)如下:

式中,{uk}={u1,u2,…uk}為k個(gè)模態(tài)分量,{ωk}={ω1,ω2…ωk}為k個(gè)模態(tài)分量對(duì)應(yīng)的頻率中心。通過(guò)引入二階懲罰因子α和拉格朗日乘子λ(t)將約束問(wèn)題變?yōu)闊o(wú)約束變分問(wèn)題,則可得:

L({uk},{ωk},λ)=

(7)

式中α為懲罰因子;λ(t)為拉格朗日乘子。通過(guò)式(7)更新各分量和中心頻率得到無(wú)約束模型的鞍點(diǎn),即最佳解。選用WOA優(yōu)化算法來(lái)確定VMD分解中最佳模態(tài)數(shù)k與懲罰因子,選用包絡(luò)熵作為WOA優(yōu)化的適應(yīng)度函數(shù),包絡(luò)熵代表原始信號(hào)的稀疏特性,當(dāng)內(nèi)涵模態(tài)分量中噪聲較多時(shí)包絡(luò)熵值較大,反之,則包絡(luò)熵值較小,包絡(luò)熵?cái)?shù)學(xué)模型(Mirjalili et al., 2016)為:

(8)

式中,N為信號(hào)的采樣點(diǎn)數(shù),Pj是α(j)的歸一化形式,α(j)為信號(hào)x(j)經(jīng)Hilbert解調(diào)后得到的包絡(luò)信號(hào)。

2 軟件結(jié)構(gòu)設(shè)計(jì)與模塊功能

2.1 軟件總體框架設(shè)計(jì)

顧及長(zhǎng)時(shí)間序列中存在的噪聲模型特性及融合多種GNSS時(shí)間序列降噪方法,利用本文提出的鯨魚(yú)變分模態(tài)奇異譜分析新方法,采用MATLAB圖形用戶界面設(shè)計(jì),結(jié)合MATLAB圖形用戶界面實(shí)現(xiàn)了WOA-VMD-SSA、EMD、EEMD、CEEMDAN四種GNSS時(shí)間序列降噪方法。軟件設(shè)計(jì)包括仿真數(shù)據(jù)模塊、時(shí)間序列降噪模塊、結(jié)果評(píng)價(jià)模塊。軟件總體設(shè)計(jì)結(jié)構(gòu)如圖1所示。

圖1 軟件總體結(jié)構(gòu)

軟件各個(gè)模塊設(shè)置可視化界面,可獨(dú)立運(yùn)行,對(duì)mom格式文件可進(jìn)行批處理,相關(guān)源程序及其數(shù)據(jù)輸出結(jié)果在對(duì)應(yīng)文件夾中,運(yùn)行輸出結(jié)果清晰可靠,各降噪模型模塊功能完整,操作簡(jiǎn)便,能有效完成GNSS時(shí)間序列降噪的數(shù)據(jù)處理工作。

2.2 仿真數(shù)據(jù)模塊

GNSS時(shí)間序列中最初存在的噪聲被假定為白噪聲,通過(guò)對(duì)長(zhǎng)時(shí)間序列的研究發(fā)現(xiàn),只假定白噪聲會(huì)導(dǎo)致參數(shù)估計(jì)偏差(魯鐵定等,2020),對(duì)GNSS站速度及其不確定度估計(jì)不準(zhǔn)確等。基于此,本軟件設(shè)計(jì)采用白噪聲(WN)、白噪聲+閃爍噪聲(WN+FN)、白噪聲+冪律噪聲(WN+ PL)、白噪聲+隨機(jī)游走噪聲+閃爍噪聲(WN+FN+RW)、白噪聲+高斯馬爾科夫(WN +GGM)五種噪聲模型,采用經(jīng)典的Trajectory Model模型對(duì)GNSS坐標(biāo)時(shí)間序列坐標(biāo)序列進(jìn)行仿真,其表達(dá)式如式(9)所示,并分析GNNS時(shí)間序列中存在的周期信號(hào)與時(shí)變信號(hào),在進(jìn)行噪聲模型估計(jì)時(shí),周期信號(hào)分析與時(shí)變分析時(shí)參數(shù)根據(jù)實(shí)際需求進(jìn)行設(shè)置(戴國(guó)強(qiáng)等,2022;張恒璟等,2014)。

y(ti)=a+bti+csin (2πti)+dcos (2πti)+

Tpost)/τj)H(ti-Tpost)+vi

(9)

式中,ti(i=1,2,…,N)為觀測(cè)時(shí)間(a),a、b分別為觀測(cè)站坐標(biāo)時(shí)間序列起始位置(mm)和速度(mm/a),c、d表示測(cè)站周年項(xiàng)運(yùn)動(dòng)(mm/a),e、f表示測(cè)站半周年項(xiàng)運(yùn)動(dòng)(mm/a),H表示階躍函數(shù),gj代表地震發(fā)生時(shí)刻Teq發(fā)生的突變值(mm),ng表示跳變個(gè)數(shù)(nh、nk是一樣的含義),hj表示地震發(fā)生之后測(cè)站速度的變化率,kj表示震后Tpost時(shí)刻松弛或滑移位移大小(mm),τj表示測(cè)站震后松弛時(shí)間常數(shù),vi表示觀測(cè)誤差(mm)。

2.3 時(shí)間序列降噪模塊

時(shí)間序列降噪模塊主要對(duì)GNSS時(shí)間序列的四種降噪方法逐一運(yùn)行實(shí)現(xiàn)。其中鯨魚(yú)變分模態(tài)奇異譜分析(WOA-VMD-SSA)降噪方法是本軟件的核心降噪方法,針對(duì)變分模態(tài)分解中不適合的模態(tài)分解數(shù)和懲罰因子出現(xiàn)過(guò)度分解與分解不足等現(xiàn)象,利用本次設(shè)計(jì)的軟件界面,用戶只需要選擇降噪方法,即可對(duì)載入的數(shù)據(jù)進(jìn)行降噪處理,運(yùn)行結(jié)束后軟件自行繪制圖形,其他三種降噪方法與WOA-VMD-SSA方法的界面操作一致,在數(shù)據(jù)處理完成后,可對(duì)四種方法的降噪結(jié)果進(jìn)行對(duì)比分析。

2.4 結(jié)果評(píng)價(jià)模塊

GNSS時(shí)間序列結(jié)果評(píng)價(jià)采用均方根誤差(RMSE)、信噪比(SNR)、相關(guān)系數(shù)(CORR)等進(jìn)行評(píng)價(jià),相關(guān)系數(shù)越接近1、均方根誤差越小,信噪比值越大,則說(shuō)明GNSS時(shí)間序列降噪效果越好。本文采用均方根誤差、信噪比、相關(guān)系數(shù)公式(賀小星,2020;付杰,2022; 惠振陽(yáng),2020),分別為式(10)、(11)、(12)。

(10)

(11)

(12)

3 方法算例實(shí)現(xiàn)與分析

3.1 軟件界面

為檢驗(yàn)本軟件的可靠性實(shí)用性,利用武漢大學(xué)IGS數(shù)據(jù)中心及http://garner.ucsd.edu/pub/measuresESESES_products/Timeseries/網(wǎng)站發(fā)布的1 000組站點(diǎn)數(shù)據(jù)進(jìn)行批處理測(cè)試,限于篇幅,本次列舉5個(gè)GNSS站的計(jì)算結(jié)果。因數(shù)據(jù)量較大,本文所用設(shè)備是處理器為i7-12700H,CPU為2.70 GHz,機(jī)帶RAM為32 GB的64位Windows11系統(tǒng)。經(jīng)過(guò)測(cè)試,軟件運(yùn)行環(huán)境良好,精度結(jié)果可靠,軟件主界面如圖2所示。

圖2 軟件主界面

3.2 算例分析

通過(guò)仿真數(shù)據(jù)對(duì)周期信號(hào)和時(shí)變信號(hào)進(jìn)行噪聲模型估計(jì),分別采用WN、WN+FN、WN+PL、WN+FN+RW、WN+GGM五種噪聲模型進(jìn)行估計(jì)。根據(jù)實(shí)際需求,在圖2軟件界面點(diǎn)擊噪聲模型模塊并設(shè)置周期與時(shí)變信號(hào)的參數(shù),輸出結(jié)果如圖3所示。

針對(duì)GNSS時(shí)間序列降噪方法多樣化,本軟件設(shè)計(jì)的時(shí)序降噪模塊為核心模塊,通過(guò)WOA-VMD-SSA、EMD、EEMD、CEEMDAN四種方法進(jìn)行降噪,分別選用6個(gè)GNSS站2010—2021年的數(shù)據(jù)進(jìn)行降噪,再與原始數(shù)據(jù)對(duì)比分析。EMD、EEMD、CEEMDAN三種方法降噪結(jié)果如圖4所示。

圖4 EMD、EEMD、CEEMDAN降噪結(jié)果圖

在進(jìn)行WOA-VMD-SSA降噪時(shí),先輸入所需參數(shù),采樣數(shù)為3 652個(gè),采樣頻率為1 Hz(李萌等,2020),振幅周期信號(hào)如圖5所示,最優(yōu)解收斂變換如圖6所示, VMD模態(tài)分解頻譜圖如圖7所示,WOA-VMD-SSA降噪前后結(jié)果對(duì)比如圖8所示。

圖5 振幅周期信號(hào)

圖6 最優(yōu)解收斂變換

圖7 VMD模態(tài)分解頻譜圖

圖8 WOA-VMD-SSA降噪前后結(jié)果對(duì)比

軟件通過(guò)時(shí)間序列處理模塊的解算,實(shí)現(xiàn)了多方法融合降噪,驗(yàn)證了不同的降噪方法將GNSS時(shí)間序列進(jìn)行降噪的可行性,事后對(duì)降噪結(jié)果進(jìn)行對(duì)比分析,降噪結(jié)果可靠。算法解算完成后,對(duì)4種方法進(jìn)行結(jié)果評(píng)價(jià),此時(shí),在結(jié)果評(píng)價(jià)模塊,直接點(diǎn)擊相關(guān)評(píng)價(jià)模式,即可運(yùn)行并繪制相關(guān)結(jié)果圖。

均方根誤差值越小、相關(guān)性系數(shù)值越接近1、噪比越大則降噪效果越好(姚文偉等,2012),本例中選用AC07、AV08、AB17、AHID、BBDM 、FALK站等6個(gè)GNSS站,對(duì)比分析四種降噪方法精度結(jié)果評(píng)價(jià)如表1所示。

表1 GNSS站精度結(jié)果評(píng)價(jià)

由表1可知,利用WOA-VMD-SSA方法降噪的各項(xiàng)精度評(píng)價(jià)結(jié)果最優(yōu)。WOA-VMD-SSA方法估計(jì)的RMSE值約為EEMD、EMD估計(jì)值的0.5～1.0,相關(guān)系數(shù)更接近于1.0,WOA-VMD-SSA方法估計(jì)的CORR值比EEMD、EMD估計(jì)值更精準(zhǔn),表明了該方法能夠較好地改善GNSS降噪效果。同時(shí),本軟件也實(shí)現(xiàn)了多方法對(duì)比的GNSS數(shù)據(jù)處理與降噪分析的目標(biāo)。

4 結(jié)論

(1)軟件可視化界面清晰,操作簡(jiǎn)易便捷,具有良好的穩(wěn)定性。

(2)軟件實(shí)現(xiàn)了WOA-VMD-SSA方法融合EMD、EEMD、CEEMDAN四種GNSS時(shí)間序列降噪方法,實(shí)現(xiàn)了WOA-VMD-SSA新方法的降噪運(yùn)行,可根據(jù)實(shí)際需要與其他三種降噪方法對(duì)比分析降噪結(jié)果。

(3)WOA-VMD-SSA新方法能獲取更高精度的GNSS時(shí)間序列降噪結(jié)果,能較好改善GNSS降噪效果,并以可視化圖形輸出,為進(jìn)一步探索GNSS基準(zhǔn)站坐標(biāo)時(shí)間序列特性研究提供可靠的參考資料。