利用CEEMDAN进行GNSS-MR雪深反演

2022-01-11薛张芳刘立龙吴昊舰刘睿国

大地测量与地球动力学 2022年1期

薛张芳刘立龙吴昊舰张志刘睿国

1 桂林理工大学测绘地理信息学院，桂林市雁山街319号，5410062 广西空间信息与测绘重点实验室，桂林市雁山街319号，541006

基于多路径效应的全球导航卫星系统多路径反射测量(GNSS multipath reflectometry，GNSS-MR)技术具有时间分辨率高、全天候、实时自动化、覆盖面积广等优点，在积雪深度探测方面具有较大潜力。张双成等[1]利用实测GPS数据对基于信噪比(signal-to-noise, SNR)的GNSS-MR探测雪深的算法进行了初步验证。Ozeki等[2]提出无几何距离的线性组L4观测值反演积雪深度，与SNR观测量探测的结果较为一致。但上述研究均采用二次多项式拟合方法获取残差序列，无法准确拟合SNR信号趋势，同时由于卫星高度角的增加，使天线接收的反射信号来自不同地表和位置，干扰有效反射信号，导致反演结果极易出现跳变现象。本文利用自适应噪声的完全集合经验模态分解[3](complete ensemble empirical mode decomposition with adaptive noise, CEEMDAN)代替二次多项式获取残差序列，并以美国科罗拉多州NWOT测站数据为例进行雪深反演，验证本文方法的有效性。

1 GNSS-MR反演雪深原理

GNSS天线接收机可以同时接收卫星发射的直射信号Ad和经地表面反射的反射信号Ar。图1是基于信噪比观测值的GNSS-MR雪深探测示意图。

图1 GNSS-MR雪深探测示意图Fig.1 Diagram of GNSS-MR for snow depth

直射信号与反射信号的振幅存在如下关系：

Ad≫Ar

(1)

记Ac为合成信号的振幅,φ为直射信号与经地表积雪面反射的信号的夹角，其关系为[4]：

(2)

由于Ad≫Ar，通常采用低阶多项式拟合方法来去除趋势项，获得低卫星高度角下SNR残差序列[1]。多路径效应下的反射信号振幅可表示为：

(3)

式中，λ为载波波长，θ为卫星高度角，h为垂直反射距离。取f=2h/λ，对SNR残差序列进行Lomb-Scargle频谱分析，可得到GPS 多路径反射信号的频率f，从而得到h；再结合图1，可获取积雪厚度。

2 CEEMDAN

CEEMDAN是对经验模态分解(empirical mode decomposition, EMD)改进后得到的一种变换形式[5],其每个本征模态函数(intrinsic mode function, IMF)对应一个频率，具有各自的物理意义。因此，可通过CEEMDAN将原始SNR信号分解成各阶IMF并重构，从而获取高质量的残差序列。

设原始SNR信号为x(t)，Ej(·)是对信号进行EMD分解后的第j阶IMF分量，ωi(t)为第i次添加的高斯白噪声(i=1，2，…，I)，εk为幅值(其中k=0，1，…，K，K是模态总体数量)。CEEMDAN分解步骤如下：

1)与EEMD[6]算法相同，在x(t)中加入白噪声εkωi(t)，对其进行EMD分解得到多个IMF分量并求集合平均值。第1阶最终分量IMF1(t)为：

(4)

2)当k=1时，计算1阶残差r1(t):

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

(5)

3)对1阶残差r1(t)添加经EMD分解后的噪声分量，得到一个新信号，即r1(t)+ε1E1[ωi(t)]，再对其分解，得到的集合平均值为该阶的IMF最终分量：

IMF2(t)=

(6)

4)当k=2，…，K时,计算第k阶残差rk(t)：

首先，要多主体参与教学评价环节，学习者、管理者、企业、教师多角度、多方面进行教学评价，对评价信息及时反馈，实行教学改进，引导教学贴近大众多样化需求；

rk(t)=rk-1(t)-IMFk(t)

(7)

5)对k阶残差rk(t)添加噪声分量εkEk[ωi(t)]，获得一个新信号，计算集合平均值，得第k+1阶的最终IMF为：

IMFk+1(t)=

(8)

6)重复第4步，直到获得的残差不可再被分解。残差最后满足：

(9)

x(t)最终被分解为K个IMF分量和1个残余分量：

(10)

3 实验分析

本文观测数据来自美国科罗拉多州NWOT观测站(https:∥www.unavco.org)，地理坐标为(40°03′19.4″N,105°35′25.8″W) ，接收机类型为大地测量型Trimble NetRS，天线类型为TRM41249，测站高3 m，数据采集间隔为5 s。实测数据来自U-C气象站(40°01′48″N,105°34′48″W)，其距NWOT测站2.2 km，由美国国家气候中心提供(https:∥gis.ncdc.noaa.gov/maps/ncei/summaries/daily )。

本实验收集了NWOT测站2016-11～2017-03连续125 d的GPS观测原始数据。为对比CEEMDAN算法和二次多项式拟合算法反演雪深的精度差异，以平均误差(ME)、均方根误差(RMSE)和相关系数(R)作为精度评价指标：

(11)

限于篇幅，仅以2016-12-16 GPS 32号卫星原始信噪比序列为例说明，截取高度角为5°～25°，波段为L2，利用CEEMDAN对其分解，结果如图2所示。分解获得的残余分量就是SNR信号的趋势项，作为自适应结果，克服了二次多项式拟合带来的振幅异常情况。

图2 CEEMDAN分解SNR序列Fig.2 Decomposition of SNR sequence by CEEMDAN

考虑到IMF分量与原信噪比序列的相关性[7]和振幅变化趋势，表1列出了部分IMF分量与原始信号的相关系数。可以看出，对于不同年积日和不同卫星，对应的IMF个数及与原始信号的相关性也不同。经过对比，在2016-12-16将IMF6和IMF7组合作为有效残差序列，随后对其组合进行Lomb-Scargle谱分析得到频率f。图3是在2016-12-16用2种方法获取的残差序列，可以看出，二次多项式拟合方法在高度角正弦值为0.17～0.19之间出现高频振荡时没能准确地拟合出SNR变化趋势；而CEEMDAN方法提取的残差序列更为光滑，无高频振荡部分。另外，该方法避免了模态混叠现象，能有效防止其他信号的干扰。图4表示2种方法获取信噪比残差序列的Lomb-Scargle谱分析图，横轴表示天线相位中心到积雪表面的距离，纵轴表示SNR反射信号频谱振幅，位于振幅最高峰值对应的距离就是垂直反射距离。由图可见，采用二次多项式方法时，最高峰值低于准确峰值，使得反演结果出现负值；雪深实测值为0.61 m，本文方法反演值为0.597 m，更接近实测雪深值。

表1 不同IMF分量与原始信号的相关系数

图3 2种方法获得的残差序列Fig.3 Residual sequences obtained by two methods

图4 SNR残差序列Lomb-Scargle谱分析图Fig.4 Lomb-Scargle spectrum analysis diagramof SNR residual sequence

统计3颗GPS卫星在2016-11-16～2017-03-21的观测数据，卫星运行周期为11 h 58 min，1 d内至少有2次经过同一测站上空，因此会出现多个SNR数据，以均值作为当天有效雪深值，2种方法的雪深反演结果如图5所示。由图5(a)看出，PRN28卫星在2016-12-06、2016-12-13～14和2016-12-24～27雪深值均为负值，尤其在12-06雪深探测值为-0.18 m，与其他反演值相差太大，且在其他卫星也出现该异常现象，应剔除该异常值，将其余雪深值作为有效雪深探测值。相比之下，图5(b)反演结果更佳，有效避免了异常跳变的情况。图6表示2种方法在多星融合下反演的雪深变化，可以看出，利用CEEMDAN进行GNSS-MR雪深反演的雪深值与实测雪深有较高的吻合度。结合表2，相比二次多项式拟合方法，CEEMDAN均方根误差降低了30.7%，与实测值相关系数为0.965，该结果进一步说明本文方法用于雪深探测具有良好的精度。