顾及鲸鱼变分模态奇异谱分析的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基准站坐标时间序列特性研究提供可靠的参考资料。