基于多参数自适应VMD 的GNSS 形变监测序列分解

高 旺，龚舒宁，潘树国，倪江生，李慧生

（1.东南大学 仪器科学与工程学院，南京 210096；2.深圳市北斗云信息技术有限公司，深圳 518057）

大型构筑物（高层楼房、桥梁、大坝等）易受到各种因素影响而产生振动与变形，进而影响构筑物的结构健康状况。结构健康监测（Structural Health Monitoring）技术[1]通过周期性记录构筑物的形变信息，并对其变形特征进行分析，进而判断结构的健康状态。基于全球导航卫星系统（Global Navigation Satellite System,GNSS）的形变监测系统凭借其较高的监测精度以及监测高频结构振动的能力，在大多数构筑物的结构健康监测系统中成为必要的组成部分[2]。

GNSS 形变监测序列的特征提取与分析是构筑物结构健康监测技术中的关键问题，其中主要包括：①对构筑物的变形趋势提取；②对构筑物的振动特性提取，如对其自振频率的提取以及各种外界影响因素导致的振动提取。由于监测序列中序列趋势、各类振动分量与观测噪声相互混叠，仅对原始监测序列进行分析无法实现对各类变形特征的区分与定量监测，需要对序列进行分解，从而提取序列中不同变形特征。

现有的序列分解方法主要有小波分解[3]（Wavelet Decomposition,WD）、经验模态分解[4（]Empirical Mode Decomposition,EMD）与奇异谱分析[5]（Singular Spectrum Analysis,SSA）等。小波分解通过将信号分解为不同频率的小波基函数来分析信号特性，但其缺点是对小波基函数的选取敏感度较高，对于成分复杂多变、难以预知信号特征的GNSS 时间序列，选择错误的小波基会使分解效果大大下降。经验模态分解主要利用循环筛选的方法对序列进行递归分解，但该算法不可避免地存在模态混叠与端部效应问题。奇异谱分析基于对时间序列的动力重构进行序列分解，但存在严重的端部效应与过分解问题。

变分模态分解[6（]Variational Mode Decomposition,VMD）算法主要基于维纳滤波、希尔伯特变换和混频的外差解调，通过构造并求解一个约束变分问题，将原信号分解为指定数量的本征模态函数（Intrinsic Mode Function,IMF）分量。该算法在分解效果方面显著优于小波分解与经验模态分解等方法，但需要对分解模态数、惩罚因子、拉格朗日惩罚算子等多组分解参数进行设定。若参数设定不当，则会显著影响分解结果的精度与可靠性。

针对VMD 算法中分解参数设定的问题，现有的部分研究仅对分解模态数这一参数进行优化。如文献[7]采用能量差参数确定分解模态数，文献[8]利用样本熵、中心频率比以及相关系数确定分解模态数，文献[9]采用峰度和包络谱峰度确定分解模态数等；这类优化方法忽略了惩罚因子等参数对分解效果的影响，无法获得最优的分解效果。

关于惩罚因子与分解模态数的综合优化，现有研究往往通过设定目标函数，并引入多参数优化算法以确定参数组合。如文献[10]设定最小包络熵为目标函数，采用灰狼优化算法进行参数优化；文献[11]采用排列熵与互信息设定目标函数；文献[12]采用融合影响指数为目标函数，引入人工蜂群算法进行参数寻优；文献[13]基于粒子群优化算法与最大熵模型进行参数优化。

此类设定目标函数的参数优化方法往往仅适用于某种特殊频域特性的时间序列，针对频域特性多变的GNSS 时间序列，目标函数的鲁棒性不足，无法保证良好的分解效果。同时，VMD 算法中初始中心频率参数的优化以及及拉格朗日乘子参数的设定往往被忽略，没有针对这两组参数的相关研究。

针对上述问题，本文从原始序列以及分解结果的频域特性出发，自适应调整分解模态数、惩罚因子、初始中心频率及拉格朗日乘子四组参数，建立多参数自适应 VMD（Multi-parameter Adaptive-VMD,MA-VMD）算法，并通过仿真序列与实测GNSS 序列对算法的分解精度和可靠性进行验证。

1 VMD 算法特性分析

1.1 VMD 算法原理

算法分解得到的IMF 为一个有带宽限制的调幅-调频函数uk(t)，表达式为：

其中，k为IMF 的序号；t为历元标识；Ak(t)为IMF的幅值函数；ϕk(t)为频率函数。

算法首先基于分解模态数M与预估各IMF 的中心频率{ωk}，利用Hilbert 变换将各IMF 的频谱调制到相应的基频带，并求其梯度平方的范数，则可建立受约束的变分问题，如式(2)所示。

其中，{uk}代表分解得到的M个IMF；{ωk}表示各IMF 对应的中心频率；∂t表示对时间t求偏导；δ(t)表示单位冲激函数；j表示单位虚数；f(t)表示原时间序列。

为了求解该变分问题，引入惩罚因子α和拉格朗日乘子λ(t)，将约束性变分问题转化为非约束性变分问题。增广拉格朗日的表达式为：

利用交替方向乘子迭代算法（Alternate Direction Method of Multipliers,ADMM）对式(3)进行求解：

①根据设定的分解模态数M、初始中心频率及初始拉格朗日乘法算子1对各个IMF 进行初始化，得到

其中，τ为拉格朗日乘子的更新倍率。

④ 重复②③，直至满足精度判据。

其中，ϵ 为精度阈值。

1.2 VMD 算法特性分析

在VMD 算法中，可以对分解模态数M、惩罚因子α、初始中心频率与拉格朗日乘子λ(t)进行调整，参数变化会直接影响序列分解结果的精度与准确性。

1）拉格朗日乘子λ(t)

由式(4)与式(6)可见，拉格朗日乘子λ(t)的主要作用是保证分解出的各个IMF之和与原序列的差值尽量接近0。该参数主要应用于信噪比较高的序列，对于低信噪比序列，λ(t)的引入会使VMD 算法的迭代次数大大增加，且分解出的IMF 中包含过多噪声。

2）惩罚因子α

由式(4)可见，IMF 在频域中偏离中心频率{ωk}的部分会受到削弱，而惩罚因子α在算法中负责控制削弱作用的强度，即α越大，整体的削弱作用越强，IMF的带宽越低。在其他参数不变的情况下，若α过低，则各个IMF 中会包含过多的噪声；若α过高，则提取出的各个IMF 带宽与幅值会出现缩减，导致部分有效信息残留在分解残差中无法被有效提取。

3）分解模态数M

分解模态数M控制算法最终分解得到的IMF 数量。当分解模态数M过小时，由于分解出的IMF 数量小于序列实际包含的信号分量数量，故需要采用远小于最佳值的惩罚因子α。通过扩大IMF 带宽的方式使得序列中的有效信息被完全提取，但由于IMF 数量过小，故分解结果中会出现模态混叠现象（即某个IMF中同时含有多个不同频率的信号分量），同时某些IMF基本仅由随机噪声组成。

当分解模态数M过大时，由于分解出的IMF 数量大于序列实际包含的信号分量数量，故分解结果中会出现过分解现象（即某一个信号分量被分解成多个IMF）。

2 MA-VMD 算法

基于VMD 算法的上述特性，为了达到最佳的分解效果，首先需要确定最佳的分解模态数M与惩罚因子α，同时还需确定合适的初始中心频率，避免算法出现模态混叠、过分解等问题。因此，本文提出了一种根据时间序列的频域特性自适应确定分解参数的MA-VMD 算法。通过对VMD 分解结果与残差进行频域分析，自适应调节惩罚因子α、分解模态数M与初始中心频率，使算法提取出尽可能多的有效信息，同时使分解残差基本只含随机误差等噪声。

首先，根据序列本身的频域特性自适应确定分解阈值T，认为原序列在频域中大于分解阈值的部分包含有效信息，小于分解阈值的部分为噪声。然后，利用当前参数对序列进行VMD 分解，并对分解结果的频域特性进行分析，判断是否分解成功、是否产生过分解等，由此更新下次分解的各个参数，直至求得最优分解结果。

图1 为MA-VMD 算法的流程图，具体步骤为：

图1 MA-VMD 算法流程图Fig.1 Flow chart of MA-VMD algorithm

1）初始化

对于GNSS 时间序列，由于序列中较高频的部分基本仅含随机噪声，故可计算序列高频部分（对于1 Hz采样率的序列，则为0.4～0.5 Hz 部分）在频域的均值μ与标准差σ，将分解阈值设置为T=μ+5σ。

将初始的分解模态数M设置为1，采用原序列的峰值频率ωpeak（频谱最大值所在频率）作为算法的初始中心频率。由于GNSS 时间序列为含噪序列，故将λ(t)设置为0。

2）VMD 分解

利用当前参数，对原序列进行VMD 分解，分解得到各个IMF 序列与残差序列。

3）分解结果分析

①峰值频率计算：首先对分解得到的各个IMF进行快速傅里叶变换（Fast Fourier Transform,FFT），并计算各IMF 的峰值频率{ωpeak}。

② 过分解判定：若两个IMF 序列的峰值频率{ωpeak}的比值过大或差值过小，则认为出现了过分解，跳转至步骤4）；

③残差频谱分析：若此时分解残差在频域的幅值均小于分解阈值，则认为残差中仅含噪声，此次分解成功，跳转至步骤4）；

④ 寻找新峰值频率：若残差频谱中存在幅值大于分解阈值的部分，则找出幅值最大且与各IMF 峰值频率不同的新频率ωnew，并将各IMF峰值频率{ωpeak}与ωnew一起组成下次VMD分解的初始中心频率，将模态数M+1，重新进行步骤2）；若找不到新频率ωnew，则认为产生了过分解现象，跳转至步骤4）。

4）二分法迭代优化惩罚因子

①若产生了过分解，则记录产生过分解的最小惩罚因子αF，并将惩罚因子向成功分解方向移动αnew=（αG+α）2；若分解成功，则记录成功分解的惩罚因子αG、分解模态数M与初始中心频率，并将惩罚因子向过分解方向移动

② 将分解模态数M重置为1，用新的惩罚因子αnew开始步骤2）。若惩罚因子的改变步长小于最小步长 Δαmin，则最后一次成功分解的惩罚因子、分解模态数与初始中心频率即为最终确定的分解参数。

3 仿真与实验验证

为了验证本文算法的有效性，分别采用一组仿真序列与桥梁GNSS 测站的实测序列对分解效果进行分析，以验证MA-VMD 算法的多参数优化效果以及算法的序列分解精度。实验中MA-VMD 算法的初始惩罚因子αstart=5000，惩罚因子的最小步长Δαmin=1000。

由于构成仿真序列的各个真实信号分量已知，故可采用互相关系数（Correlation）与均方根误差（Root Mean Square Error,RMSE）作为信号分解精度的评价指标。

其中，ρ为互相关系数；σ为均方根误差；fk(t)为第k个真实信号分量分别为真实序列分量与分解结果分量的均值；N为序列长度。

在仿真实验中，采用MA-VMD 算法对仿真序列进行分解，并对MA-VMD 算法确定的各组分解参数进行偏移，利用偏移后的参数进行VMD 分解，从而验证MA-VMD 算法的多参数优化效果。同时采用EMD、SSA 算法与文献[10]提出的改进变分模态分解算法（Improved Variational Mode Decomposition,IVMD）对同一组仿真序列进行分解，从而验证MA-VMD 算法的优越性。

在实测数据实验中，由于无法求得各项信号分量的真值，故将MA-VMD 算法的分解结果与原序列在时域、频域两方面进行对比，并对不同时段、不同方向上的GNSS 坐标序列进行分解，以验证算法鲁棒性。

3.1 仿真序列实验与分析

首先建立一组仿真GNSS 坐标时间序列。仿真序列在考虑GNSS 时间序列中可能存在的多路径效应等周期性误差、观测噪声及粗差的基础上采用GNSS 变形分析中常用的坐标时间序列模型[14]。模型可表示为式(10)。

其中，x(ti)为仿真时间序列；F(ti)为序列中的趋势部分，采用随机样条曲线；m0为谐波个数；am、bm为信号中多路径误差等周期项的振幅；fm为周期项频率；为随机噪声；err为粗差。

构成仿真序列的周期分量共有8 个，频率分别为0.005 Hz、0.01 Hz、0.1 Hz、0.12 Hz、0.2 Hz、0.3 Hz、0.35 Hz；周期分量的振幅均为随机生成。采用式(10)所示的坐标时间序列模型随机生成100 组仿真时间序列，并用MA-VMD 算法对仿真序列进行分解。以其中一组仿真序列为例进行分析，如图2 所示。图2(a)为该仿真序列的时域图，图2(b)为构成该仿真序列的真实分量。

图2 仿真实验中的时间序列Fig.2 Simulated sequence in simulation experiment

图3 为该仿真序列的MA-VMD 分解结果。图3(a)为MA-VMD 算法分解得到的9 个IMF 分量的时域序列与频率谱，其中前8 个IMF 为不含趋势的周期项，第9 个IMF 为序列的趋势分量。MA-VMD 算法在进行 16 次迭代后确定最优参数：最优分解模态数M=9，初始中心频率为{0,0.005,0.01,0.1,0.12,0.15,0.2,0.3,0.35} Hz，分解模态数及各分量的初始中心频率均与真实构成序列的各分量相同；最优惩罚因子为α=2144780，惩罚因子迭代过程如图3(b)所示。将分解得到的趋势分量与原序列对比，如图3(c)所示，MA-VMD 算法分解得到的趋势分量与真实趋势分量的吻合程度较好，仅在序列起始处受到周期项波动的影响，表明MA-VMD 算法可有效提取序列趋势项。从图3(a)中分解残差的频谱可见，分解残差中仅由随机噪声构成，不含有效信息；同时，算法有效提取了8 个周期分量，在频率与幅值方面均与图2(b)所示的真实分量相匹配，且没有出现模态混叠与过分解现象。计算提取出的各个分量与真实分量的互相关系数与RMSE，互相关系数的均值为0.9877，RMSE 均值为0.1365，说明分解结果与真实分量的吻合程度较好。

图3 仿真序列的MA-VMD 分解结果Fig.3 MA-VMD decomposition result of simulated sequence

为了探讨算法将拉格朗日乘子λ(t)固定为0 的必要性，在其他参数不变的情况下，设置不同的λ(t)进行VMD 分解，解算结果如表1 所示。由表1 可知：当λ(t)=0时，VMD 算法仅迭代9 次就得到了较好的分解效果；当λ(t) ≠ 0时，算法迭代次数随着λ(t)的增大而大幅增加，且由于各IMF 包含过多噪声，导致分解的准确程度下降。由此可见，对于序列特性类似于GNSS 时间序列的含噪序列，将λ(t)设置为0 有助于提升算法的效率与分解精度。

表1 不同拉格朗日乘子对算法分解效果的影响Tab.1 Effect of different Lagrange multipliers to decomposition

图4 无初始中心频率约束时的VMD 分解结果Fig.4 VMD decomposition result without initial center frequency constraint

为进一步验证MA-VMD 算法对各分解参数的优化效果，对生成的100 组时间序列都进行分解。分解结果显示，100 组序列均正确确定了分解模态数M=9，且初始中心频率均与真实分量相同。

将MA-VMD 算法确定的各序列的惩罚因子α进行偏移，在其他参数不变的情况下，利用偏移后的惩罚因子对100 组序列进行VMD 分解。分别计算MA-VMD 分解结果以及参数偏移后的分解结果与真实分量之间的互相关系数均值与RMSE 均值，结果如图5 所示。由图5 可见，MA-VMD 算法确定的惩罚因子具有最高的互相关系数与最低的RMSE；随着惩罚因子偏移程度增大，分解准确程度下降明显。

图5 采取不同惩罚因子时的互相关系数与RMSEFig.5 Correlation and RMSE of different penalty factors

综合上述实验可见，对于随机生成的多组仿真序列的分解，MA-VMD 算法均能确定正确地分解模态数M，且算法确定的初始中心频率与惩罚因子α接近全局最优。

为了对比验证MA-VMD 算法相对于其他现有算法的优越性，分别采用EMD、SSA、IVMD 算法对相同的仿真序列进行分解，结果如图6 所示。其中，图6(a)为MA-VMD 算法分解结果。图6(b)为IVMD 算法的分解结果，该算法最终确定的分解模态数M=8，惩罚因子为α=316685，分解模态数与序列实际分量数不同；第6 个与第7 个分量的中心频率几乎相同，出现过分解现象；分解残差中仍含有0.1 Hz 与0.2 Hz 的振动分量未被提取。图6(c)为EMD 算法的分解结果，可见分解结果出现了明显的模态混叠与过分解现象，且分解残差中仍存在明显趋势分量，分解效果不佳。图6(d)为SSA 算法的分解结果，可见第2 个与第3 个分量、第4 个与第5 个分量、第6 个与第7 个分量的中心频率几乎相同，出现了明显的过分解现象，并且分解残差中仍存在多个频率的振动分量。

图6 现有算法的仿真序列分解结果与MA-VMD 分解结果的对比Fig.6 Comparison of decomposition result of simulated sequence between existing algorithms and MA-VMD

综合上述实验，MA-VMD 算法能够将分解模态数M、惩罚因子α以及初始中心频率收敛至接近全局最优，且分解效果优于多种现有信号分解算法。

3.2 实测数据实验与分析

为了检验算法分解效果，采用苏通大桥数据对算法分解能力进行分析。测试采用的时间序列为苏通大桥上的GNSS 测站在2022 年11 月18 日共32400 历元（9 h）的观测数据，采样频率为1 Hz；根据Nyquist定律，序列可用于分析0～0.5 Hz 的频域特征。由于交通荷载、温度变化及长江潮汐的影响，桥梁各个部分均会发生变形，且桥梁本身也存在一定的结构振动。

首先截取其中1～3600 历元的监测序列，并对Z方向的序列进行分解实验。图7(a)为该序列的时域图像，图7(b)为序列频谱图。由图7(a)可见，该监测序列存在大量观测噪声；由图7(b)可见，序列在频域中的主要能量集中在其趋势部分（频域中接近0 Hz 的部分），桥梁本身的振动特性在频域中被淹没，无法有效进行分析。

图7 苏通大桥测站Z 方向变形监测数据Fig.7 Deformation monitoring data of Sutong Bridge measurement station in Z direction

利用MA-VMD 算法对该序列进行分解，结果如图8 所示。MA-VMD 算法经过6 次迭代后确定最优惩罚因子α为750，最优分解模态数M为6。图8(a)为算法分解得到的6 个IMF 以及分解残差的时域序列与频谱，由分解残差的频谱可知：残差低频部分的幅值基本与高频部分的随机噪声幅值相同，说明分解残差中已不含有效信息。在分解结果中，IMF 6 为序列的趋势分量，由其频谱可见趋势分量的能量主要集中在0 频率附近，同时可见测站Z方向在1 h 内呈逐渐下降趋势，共下降约160 mm。

图8 苏通大桥测站Z 方向序列的MA-VMD 分解结果Fig.8 MA-VMD decomposition result of the time series from Sutong Bridge measurement station in Z direction

图8(b)为提取出的趋势分量与原序列的对比图，可见MA-VMD 算法提取出的趋势分量与原序列吻合度较高。

IMF 1～5 为序列中的五个周期分量，代表测站在不同频率范围的振动信息。由于原序列的趋势分量能量过大，影响对周期分量提取效果的分析，故首先将原序列减去趋势分量（IMF 6），得到去趋势的时间序列并求其频谱，算法提取出的各个周期分量的频谱与去趋势序列频谱的对比如图8(c)所示。由图8(c)可知：MA-VMD 算法成功提取了序列在不同频率范围的主要振动信息，各IMF 的峰值频率与原序列频谱的局部峰值相吻合。

从各周期分量的频谱与时域序列可知，IMF 2～5主要为桥梁在Z方向上由于车辆载荷等作用而产生的较宽频振动，主要分布在0.01～0.25 Hz 区间。将该频段的各个IMF 相加得到总体振动序列，该序列的平均振幅为94.34 m。除此之外，IMF 1 的峰值频率为0.309 Hz，平均振幅为21.76 mm，说明测站在0.3 Hz附近也有频率范围较窄的显著振动，应为桥梁的自振部分。

为了验证MA-VMD 算法分解的稳定性，对测站X-ECEF 与Y-ECEF 方向上相同时刻的变形监测序列进行分解，并对比去趋势序列的频谱与提取出的各个周期分量的频谱，结果如图9 所示，其中图9(a)为X方向，图9(b)为Y方向。可见X方向与Y方向上各个周期分量的频域构成与Z方向类似：算法均在0.01～0.25 Hz 区间提取出因载荷等作用而产生的宽频振动分量；在0.3 Hz 附近提取出存在较窄频带的振动分量，验证了MA-VMD 进行序列分解的稳定性。

图9 X 方向与Y 方向序列的MA-VMD 分解结果Fig.9 MA-VMD decomposition result of time series from X-direction and Y-direction

为了监测测站Z方向的长期变形振动情况，对Z方向连续9 h（共32400 历元）的监测数据按小时进行分割，并对九段序列分别进行MA-VMD 分解。分解完成后，分别统计峰值频率在0.01～0.25 Hz 频段以及0.25～0.35 Hz 频段的IMF 分量，并计算其平均振幅，计算结果如图10 所示。由图10 可知，0.01～0.25 Hz频段（因载荷产生的振动）的振幅在93 mm 左右波动，0.25～0.35 Hz 频段（自振频率）的振幅在20 mm 左右波动。两种不同频段的振幅变化趋势大致相同，第5 h的振幅均为最低。由此可见，本文算法可对形变监测序列的长期运动趋势及各类振动特性进行定量分析与持续监测。

图10 苏通大桥测站Z 方向各时段的振动特性Fig.10 Vibration characteristics of Sutong Bridge measurement station in Z direction

4 结论

本文首先对变分模态分解算法中多项参数对分解效果的影响进行了综合分析，并从原始序列以及分解结果的频域特性出发，自适应调整分解模态数、惩罚因子、初始中心频率及拉格朗日乘子四组参数，建立MA-VMD 算法。利用MA-VMD 算法对多组仿真序列及真实形变监测序列进行分解。仿真序列实验表明，MA-VMD 算法的分解结果与真实值之间的互相关系数为98.77%、均方根误差为0.1365 mm，均接近全局最优，并优于经验模态分解、奇异谱分析、改进变分模态分解等算法；同时算法对真实GNSS 形变监测序列具有优良的分解效果，能够避免模态混叠与过分解问题，正确提取序列趋势及各类振动相关频率特征。MA-VMD 算法可有效应用于构筑物变形趋势与振动特性的精细化定量分析。