GNSS浮标海潮高滤波光滑算法

2021-10-15胡亮亮王胜利刘以旭

科学技术与工程 2021年27期

胡亮亮，王进，刘欢，刘毅，王胜利，刘以旭

(1.山东科技大学测绘科学与工程学院，青岛 266510； 2.山东科技大学海洋科学与工程学院，青岛 266510)

自古以来，沿海地区是人类活动最为频繁的地方之一，人们的生活也与海洋的潮汐现象密不可分[1-3]，例如，潮汐能发电、水产养殖、农业灌水、水下设备建设等，因此分析潮汐变化具有重大意义，而全球导航卫星系统(global navigation statellite system,GNSS)技术已经成为获取浮标海潮数据的一种重要手段。

基于地基增强的动态实时定位(real time kinematic, RTK)技术保证了GNSS浮标海上高精度定位需求，但受到无线电通信和远洋环境影响，RTK无法满足远洋海洋测量要求。随着精密钟差和卫星轨道产品的优化，GNSS浮标验潮的精密单点定位(precise point positioning，PPP)与后处理RTK处理模式(post processing kinematic，PPK)技术应用越来越广泛，且PPP验潮的方法拥有无需布设基准站和不受作用距离限制等优越性，在远洋精密海潮数据勘测方面有着较大的应用场景[4-6]。

通过GNSS浮标连续观测和分析获取的海面高，包含周期性潮位，以及多重的复杂高频的涌浪信号和高斯白噪声，其时间序列信号拥有较为显著的非平稳、非线性变化特征[7]。需要通过滤波进行处理，门限滤波需要多次烦琐确定波浪周期，且涌浪频率变化微弱时不易分析[3,8]；当潮汐与涌浪频率相差较大时，可以采用低通滤波，使用窗函数算法，设置低通滤波器的截止频率，但是此方法中的阶数选择、性能比需要根据经验进行多次测试，从而获取最优参数[9-10]；小波分析是进行潮位数据滤波较为广泛的一种方法，能够较好地分离出低频的潮汐信号，但是需要噪声阈值和小波基函数，同样也需要进行大量的测试验证才可以确定[4,11]。上述方法都是应用在线性平稳的时间序列分析上，对于海面高时间序列并不适用。

为了分析非线性非平稳的时间序列信号，Huang等[12]提出经验模态分析方法(empirical mode decomposition，EMD)，分解获得本征模态分量(intrinsic mode function， IMF)包含原信号的不同时间尺度的多局部特征信号，体现信号的物理意义，但是容易出现模态混叠现象。Wu等[13]提出集合经验模态分解法(ensemble empirical mode decomposition，EEMD)通过向原始信号中添加高斯白噪声以达到削弱EMD模态混叠现象，但也会存在尺度丢失以及模态混叠现象。基于EMD和EEMD的基础下，Torres等[14]提出了自适应噪声分析进行的完全集合经验模态分解法(complete ensemble empirical mode decomposition，CEEMD)，通过在每次分解过程中添加一个拟定好的高斯白噪声，计算唯一的残差值，其有效减小了迭代次数，增加重构精度，更适应于非线性信号的分析。CEEMD能够有效削弱EMD分解中产生的模式混叠现象和冗余分量现象，比EMD具有更强的噪声鲁棒性以及更弱的端点效应[15]。变分模态分解(variational mode decomposition，VMD)将待分解的信号转化为非递归、变分模态的分解形式，能够很好地对噪声信号进行分解[16-17]。

海潮高时间序列作为一种非平稳非线性信号，Kuo等[18]提出PPK、PPP与验潮仪对比精度，结果证明GPS浮标使用PPP的RMSE达到8 cm，并且使用EMD算法可以成功地在GPS浮标观测中检测到周期为几秒到一天的高频海平面信号，并将其识别为波浪、气象海啸和潮汐，但是在实验结果中也出现模态混叠现象。文献[19]提出了基于CEEMD方法对于海岸处GNSS站高程进行处理，并得到高程变化的趋势，但此类方法目前并没有应用到潮位滤波降噪上。海潮同时作为一种缓变的自然信号，在考虑降噪的同时，也需要考虑到光滑性，因此引入熵理论[20]。排列熵(permutation entropy，PE)是一种用来检测随机性和动力学突变的方法，抗噪能力强。但是在单一尺度下分析时间序列的不规则性具有一定局限性，潘震等[21]提出多尺度排列熵(multi-scale permutation entropy，MPE),衡量时间序列在不同尺度下随机性和复杂性的方法，鲁棒性优良。

基于MPE描述时间序列的规则性优点，考虑海潮时间序列的降噪光滑特点，提出多指标多尺度排列熵(multi-index multi-scale permutation entropy，MMPE)，并将其与CEEMD算法结合用于提取最优海潮时间序列。基于实际GNSS浮标测量与分析获得海图高时间序列，分别利用CEEMD+MPPE和VMD+MMPE算法分解重构最优潮位时间序列，并与小波分析结果对比，验证提出的两种算法的可行性和准确性，且通过比较标准差、MMPE值和时频分析展现出CEEMD+MMPE算法优良的降噪光滑效果，在GNSS浮标验潮数据分析中有着更好应用趋势。

1 GNSS测量海潮高

随着不受远洋环境影响以及高频率，高精度的PPP/PPK的技术，可以提供远海潮高数据，其流程与原理分别如图1、图2所示。

图1 GNSS浮标验潮流程图Fig.1 Processing of GNSS buoy tide survey

图2 GNSS浮标验潮原理Fig.2 Principle of GNSS buoy tide survey

首先，通过PPP/PPK技术分析处理浮标天线的原始观测值，获取到WGS84坐标系下的天线的瞬时大地高H；扣除浮标姿态改正GNSS天线高垂直高度h1并加入静态、动态的吃水改正值Δh，获得海面瞬时大地高；通常是用海图高表示海面潮位，所以需要高程基准转换，将大地高转换到海图高[5]。将高程异常值添加到海面大地高获得海面正常高，再通过线性内插的方法装换到海图高。经过上述步骤可以获取海面瞬时海图高，即

T′=H-h1+Δh-(ξ+h2)

(1)

式(1)中：T′为海面瞬时海图高，即海面高，其中包含潮汐T，涌浪信号以及各种观测噪声和自然信号噪声；ξ为大地高与正常高之间的高程异常值；h2为瞬时正常高转海图高的海图深度基准。

2 时间序列信号分析

基于GNSS浮标测量技术与高程转换获取的海面瞬时海图高，其中包含涌浪和噪声，需要通过滤波获取T。CEEMD与VMD算法都是一种适用于分析非线性非平稳信号的自适应方法，将信号在不偏离时间域的情况下以相对高低频分解成为数个IMF，并基于本文提出的MMPE重构得到降噪光滑的海潮时间序列。

2.1 完全集合经验模态分解法(CEEMD)

CEEMD在EMD和EMMD基础上的改进的递归降噪算法，其减小了筛选迭代次数，进一步削弱了模态混叠现象，对于海潮高时间序列滤波具有很大优势。通过以下步骤完成IMF分层提取。

(1)对原始信号x(t)加入一个高斯白噪声ωi(t)，即为x(t)+ε0ωi(t)并对其进行n次EMD分解，平均计算获取到第一个固有本征模态函数IMF1(t)，即

(2)

式(2)中：E1()表示为x(t)通过EMD方法分解获取的第1个本征模态函数；ωi(t)(i=1,2,…,N)为单位方差为0的高斯白噪声；εk表示CEEMD方法分解出K个IMF1(t)下的信噪比系数。

将原始信号x(t)与上述IMF1(t)相减得到一阶残差r1(t)，即

r1(t)=x(t)-IMF1(t)

(3)

(2)在上述获取的一阶残差r1(t)添加高斯白噪声，信号表示为r1(t)+ε1E1[ωi(t)]，再次进行分解，获取第二个固有本征模态函数IMF2(t)，并同理将r1(t)与IMF2(t)相减得到残差r2(t)，即

(4)

r2(t)=r1(t)-IMF2(t)

(5)

(3)由此类推，继续进行k次分解，在得到信号rk(t)+εkEk[ωi(t)]，计算第k+1个本征模态函数分量以及k+1阶残差为

(6)

r(k+1)(t)=rk(t)-IMF(k+1)(t)

(7)

(4)循环计算获取本征模态分量与残差值，直至残差极值点个数不大于2时，分解结束，最终信号表示为

(8)

式(8)中：K表示在CEEMD方法分解得到的IMF分量个数。

2.2 变分模态分解

VMD是将信号分解为给定数量的有限带宽模态，是一种自动准正交的分解方法，其中获取的每个IMF分量的频带限制在各自估计的中心频率。对于每个IMF分量，需就行Hilbert变换得到单边频谱，通过指数函数将每个分量的中心频带移至基带，以L2范数梯度的平方估计IMF分量的带宽。

由此上述转化为一个约束变分问题，即

(9)

式(9)中：{uk}={u1,u2,…,uk}和{ωk}={ω1,ω2,…,ωk}分别表示的是所有IMF分量和中心频率集合；f是分解的信号；δ(t)为狄拉克函数。

采用具有良好收敛性的二次惩罚因子和严格执行约束条件的拉格朗日乘数λ(t)，重构问题变为非约束变分问题，其表达式为

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

(10)

式(10)中：α称为惩罚因子，表示每一个模态初始中心约束强度。

而后采用交替方向算法(alternate direction method of multipliers，ADAM)迭代发现更新的IMF分量和中心频率，公式为

(11)

综上所述，实际应用的VMD算法步骤如下。

(1)初始化u1、ω1、λ1。

(2)根据式(12)更新uk、ωk。

(3)通过式(12)，更新λ,即

(12)

(4)重复步骤②进行迭代计算，以满足条件

(13)

式(13)中：ε为判别精度。当满足判别精度要求时，得到了K个模态分量；否则继续返回步骤(2)继续计算，直到满足步骤(3)。

2.3 多指标多尺度排列熵

MMPE的基本思想是在MPE的基础上引入降噪度和相关系数，通过MPE反映重构信号在不同尺度下的排列熵，即规则性，并考虑与真值的降噪度和相关系数，获取最优的MMPE值，重构信号即为最优的降噪光滑时间序列。

(14)

(15)

(16)

(17)

当Pr=1/m!时，Hp(n)达到最大值ln(m!)。

多尺度排列熵Hp只能反映时间序列的光滑程度，在此基础下引入算法降噪度和相关系数，分别反映重构信号的精准度与相关性，构成多指标多尺度排列熵，即

(18)

式(18)中：α、β、1-α-β分别为多指标影响因子；0

3 试验及其分析

3.1 实验概况

本次实验为2020年11月20日在日照市码头附近进行的GNSS浮标验潮对比测试。实验过程中，在岸边设立GNSS基准站(见图3，中海达AT-53501CR GNSS扼流圈天线，中海达VNet8 主机)，在距离验潮仪附近设立GNSS浮标(见图4，中海达HZZACF-S808 GNSS扼流圈天线，NovAtel OEM6接收机)，可接收到四系统双频信号，并包括IMU测姿系统。GNSS浮标与基准站距离11 km，与验潮仪相距300 m，基准站、浮标和验潮仪的位置及采样频率见表1。

表1 GNSS 基准站、浮标和验潮站坐标Table 1 The coordinates of GNSS base station, buoy and tide gauge station

图3 GNSS基准站Fig.3 GNSS base station

图4 GNSS浮标Fig.4 GNSS buoy

GNSS浮标的数据分析主要包括了PPP处理和PPK处理。并通过图1所示的方法，将连同验潮仪测量的潮位值归算到海图基准。最终得到将3个结果进行对比(图5)，可以看出GNSS浮标测量计算得到潮位变化与验潮仪测量结果在波动趋势上大致一致，但GNSS技术得到的海潮数据包含大量的噪声值，需要剔除。

图5 GNSS浮标分析与验潮仪潮位变化Fig.5 The change of GNSSbuoy bnalysis and bide gauge

验潮仪常作为高精度的传统验潮方法，拥有广泛的认可。将验潮仪数据拟定为基准值，经统计分析得到PPP解算的精度为7.25 cm，PPK解算的精度为6.78 cm，适用于远洋的海潮测量。与传统验潮仪相比，PPP解算结果95%置信区间为±7.67 cm(图6)，PPK解算结果中95%置信区间为±6.39 cm(图7)。

图6 PPP结果与验潮仪对比误差Fig.6 Comparison error between PPP results and tide gauge

图7 PPK结果与验潮仪对比误差Fig.7 Comparison error between PPK results and tide gauge

3.2 时间序列分解

GNSS浮标测量并处理得到海面高时间序列，其包含了大量的涌浪、自然和观测噪声，其精度并不完全符合验潮要求。基于本文中提出的CEEMD和VMD两种时间序列分解算法，以PPP处理下的海面高时间序列为例进行分析，并通过傅里叶变换得到每一个IMF分量的频谱。如图8所示，对PPP技术获得海面高信号中分别加入标准差为0.03 m的白噪声，通过100次的迭代，CEEMD将信号都分解为14个IMF分量，其中包括一个residual序列，以及各个分量对应的频谱图，从图8中看出，CEEMD方法分解出来的每一个分量频带宽度大致相等，频带出依次递减且无重叠，直至获取到低频的信号，很好地避免了分解过程中容易出现的模态混叠现象。图9显示的是VMD算法分别对PPP的海面高序列分解结果,获取到IMF1～IMF6以及一个residual序列共计7个分量,每一个分量图右边对应的是对数形式的频谱图，VMD方法分解每一个分量信号都成功消除了模态混叠影响，信号频率从高依次递减，并且获取的低频信号IMF6即为海潮时间序列信息。能够看出两种方法都可以实现信号频域内各个分量的自适应分割，可以有效克服分解中产生的模式混叠现象和冗余分量。

图8 CEEMD对PPP海潮高的分解结果Fig.8 CEEMD decomposition result of PPP sea tide height

图9 VMD对PPP海潮高分解结果Fig.9 VMD decomposition result of PPP sea tide height

3.3 最优降噪光滑海潮高

通过CEEMD和VMD算法分解PPP方法下海潮高时间序列，只是获取信号本身的IMF分量，其过程并没有达到降噪的作用，需要通过重构获取海潮时间序列。基于式(18)提出的MMPE算法，分析不同IMF分量重构的时间序列值，得到HMMEP。通过图10和图11中MMPE值看出，在CEEMD方法下的最优重构信号就是residual本身，而在VMD方法下最优重构信号是IMF6分量本身。

图10 CEEMD下MMPEFig.10 MMPE of CEEMD

图11 VMD下MMPEFig.11 MMPE of VMD

在GNSS浮标验潮的降噪过程中，通常使用的是小波分析对GNSS浮标数据进行再处理。将常用的小波滤波分析的海潮时间序列与本文基于MMPE下的CEEMD或VMD分解重构的结果进行对比，同时本文也使用的PPK处理的海潮时间序列数据进行上述分析，如表2所示。

表2 3种滤波在PPP和PPK数据结果对比Table 2 Comparison of three kinds of filtering results in PPP and PPK data

基于CEEMD+MMPE方法和VMD+MMPE方法的分解重构得到的海潮时间序列精度均优于小波分析结果，通过对比PPP和PPK的分解潮位数据，均有较好的吻合性，在整体趋势上有较高的一致性。CEEMD+MMPE方法下统计结果表明：在PPK模式下的标准差为1.39 cm，相关系数为0.999 8，MMPE值为0.973 7；PPP模式下的标准差为3.0 cm，相关系数为0.998 3，MMPE值为0.965 1。通过图12和图13可以看出3种方法下得到海潮时间序列与验潮仪测量的潮位变化趋势相同，且CEEMD+MMPE算法分析的海潮序列最大误差也在5 cm，符合海洋调差规范。由实验结果得知基于GNSS浮标获取海潮高时间序列，并通过CEEMD+MMPE算法得到潮位变化信息，在降噪、相关性和光滑都与验潮仪吻合较好，符合海潮验潮的精度要求，同时也保证数据的可信度。

图12 PPP方法数据解算结果对比Fig.12 Comparison of data solution results of PPP method

图13 PPP方法数据解算结果差值对比Fig.13 Comparison of difference between data solution results of PPP method

在潮汐时间序列中，瞬时频率可以作为一种区分信号类型的有效方法。因此本文分别对小波分析、CEEMD和VMD分解提取的海潮时间序列进行Hilbert变换获得时间序列的瞬时频率。如图14所示，可以看出3种方法下的重构的海潮时间序列效果较好，整体趋势光滑，且都未出现明显的模态混叠现象，拥有优良的时频聚集性，突出了海潮时间序列的时频信息，说明CEEMD与VMD在潮位时间序列提取上有着很好的应用可行性和准确性。