顧及鯨魚變分模態奇異譜分析的GNSS時序降噪軟件設計與實現

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

(1.東華理工大學 測繪與空間信息工程學院,江西 南昌 330013; 2.東華理工大學 自然資源部環鄱陽湖區域礦山環境監測與治理重點實驗室,江西 南昌 330013;3.江西理工大學 土木與測繪工程學院,江西 贛州 341000;4.華東交通大學 交通工程學院,江西 南昌 330013)

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

1 WOA-VMD-SSA降噪新方法

針對變分模態分解中不適合的模態分解數和懲罰因子出現過度分解與分解不足等現象(陶國強,2021),筆者提出一種鯨魚變分模態奇異譜分析(WOA-VMD-SSA)新方法。SSA是一種將信號分解出多個有序的子序列,通過SSA對GNSS時間序列進行分解,從而得到其中所包含的周期信息的高頻時間序列、趨勢信息的低頻時間序列和代表噪聲的時間序列。

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

(1)

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

X=E1+E2+…Ed

(2)

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

Xi=XI1,XI2,…,XIp

(3)

(4)

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

(5)

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

將時間序列通過迭代的方式進行VMD分解,以求解最優中心頻率與帶寬,之后對信號重構并求解約束變分問題,其數學模型(Dragomiretskiy et al., 2014)如下:

式中,{uk}={u1,u2,…uk}為k個模態分量,{ωk}={ω1,ω2…ωk}為k個模態分量對應的頻率中心。通過引入二階懲罰因子α和拉格朗日乘子λ(t)將約束問題變為無約束變分問題,則可得:

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

(7)

式中α為懲罰因子;λ(t)為拉格朗日乘子。通過式(7)更新各分量和中心頻率得到無約束模型的鞍點,即最佳解。選用WOA優化算法來確定VMD分解中最佳模態數k與懲罰因子,選用包絡熵作為WOA優化的適應度函數,包絡熵代表原始信號的稀疏特性,當內涵模態分量中噪聲較多時包絡熵值較大,反之,則包絡熵值較小,包絡熵數學模型(Mirjalili et al., 2016)為:

(8)

式中,N為信號的采樣點數,Pj是α(j)的歸一化形式,α(j)為信號x(j)經Hilbert解調后得到的包絡信號。

2 軟件結構設計與模塊功能

2.1 軟件總體框架設計

顧及長時間序列中存在的噪聲模型特性及融合多種GNSS時間序列降噪方法,利用本文提出的鯨魚變分模態奇異譜分析新方法,采用MATLAB圖形用戶界面設計,結合MATLAB圖形用戶界面實現了WOA-VMD-SSA、EMD、EEMD、CEEMDAN四種GNSS時間序列降噪方法。軟件設計包括仿真數據模塊、時間序列降噪模塊、結果評價模塊。軟件總體設計結構如圖1所示。

圖1 軟件總體結構

軟件各個模塊設置可視化界面,可獨立運行,對mom格式文件可進行批處理,相關源程序及其數據輸出結果在對應文件夾中,運行輸出結果清晰可靠,各降噪模型模塊功能完整,操作簡便,能有效完成GNSS時間序列降噪的數據處理工作。

2.2 仿真數據模塊

GNSS時間序列中最初存在的噪聲被假定為白噪聲,通過對長時間序列的研究發現,只假定白噪聲會導致參數估計偏差(魯鐵定等,2020),對GNSS站速度及其不確定度估計不準確等?；诖?本軟件設計采用白噪聲(WN)、白噪聲+閃爍噪聲(WN+FN)、白噪聲+冪律噪聲(WN+ PL)、白噪聲+隨機游走噪聲+閃爍噪聲(WN+FN+RW)、白噪聲+高斯馬爾科夫(WN +GGM)五種噪聲模型,采用經典的Trajectory Model模型對GNSS坐標時間序列坐標序列進行仿真,其表達式如式(9)所示,并分析GNNS時間序列中存在的周期信號與時變信號,在進行噪聲模型估計時,周期信號分析與時變分析時參數根據實際需求進行設置(戴國強等,2022;張恒璟等,2014)。

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

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

(9)

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

2.3 時間序列降噪模塊

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

2.4 結果評價模塊

GNSS時間序列結果評價采用均方根誤差(RMSE)、信噪比(SNR)、相關系數(CORR)等進行評價,相關系數越接近1、均方根誤差越小,信噪比值越大,則說明GNSS時間序列降噪效果越好。本文采用均方根誤差、信噪比、相關系數公式(賀小星,2020;付杰,2022; 惠振陽,2020),分別為式(10)、(11)、(12)。

(10)

(11)

(12)

3 方法算例實現與分析

3.1 軟件界面

為檢驗本軟件的可靠性實用性,利用武漢大學IGS數據中心及http://garner.ucsd.edu/pub/measuresESESES_products/Timeseries/網站發布的1 000組站點數據進行批處理測試,限于篇幅,本次列舉5個GNSS站的計算結果。因數據量較大,本文所用設備是處理器為i7-12700H,CPU為2.70 GHz,機帶RAM為32 GB的64位Windows11系統。經過測試,軟件運行環境良好,精度結果可靠,軟件主界面如圖2所示。

圖2 軟件主界面

3.2 算例分析

通過仿真數據對周期信號和時變信號進行噪聲模型估計,分別采用WN、WN+FN、WN+PL、WN+FN+RW、WN+GGM五種噪聲模型進行估計。根據實際需求,在圖2軟件界面點擊噪聲模型模塊并設置周期與時變信號的參數,輸出結果如圖3所示。

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

圖4 EMD、EEMD、CEEMDAN降噪結果圖

在進行WOA-VMD-SSA降噪時,先輸入所需參數,采樣數為3 652個,采樣頻率為1 Hz(李萌等,2020),振幅周期信號如圖5所示,最優解收斂變換如圖6所示, VMD模態分解頻譜圖如圖7所示,WOA-VMD-SSA降噪前后結果對比如圖8所示。

圖5 振幅周期信號

圖6 最優解收斂變換

圖7 VMD模態分解頻譜圖

圖8 WOA-VMD-SSA降噪前后結果對比

軟件通過時間序列處理模塊的解算,實現了多方法融合降噪,驗證了不同的降噪方法將GNSS時間序列進行降噪的可行性,事后對降噪結果進行對比分析,降噪結果可靠。算法解算完成后,對4種方法進行結果評價,此時,在結果評價模塊,直接點擊相關評價模式,即可運行并繪制相關結果圖。

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

表1 GNSS站精度結果評價

由表1可知,利用WOA-VMD-SSA方法降噪的各項精度評價結果最優。WOA-VMD-SSA方法估計的RMSE值約為EEMD、EMD估計值的0.5～1.0,相關系數更接近于1.0,WOA-VMD-SSA方法估計的CORR值比EEMD、EMD估計值更精準,表明了該方法能夠較好地改善GNSS降噪效果。同時,本軟件也實現了多方法對比的GNSS數據處理與降噪分析的目標。

4 結論

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

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

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