ICEEMDAN在GNSS-MR海平面高度反演中的应用

Application of ICEEMDAN in GNSS-MR Sea Surface Height Inversion

摘要: 针对全球导航卫星系统多路径反射测量(GNSS multipath reflectometry, GNSS-MR)技术存在的信噪比(Signal-to-Noise Ratio, SNR)信号中掺杂噪声的问题，提出利用改进的自适应噪声的完全集合经验模态分解方法对原始信噪比数据进行分解，筛选出不含噪声的有效残差序列，实现去噪处理，再应用于GNSS-MR技术反演海平面高度。以西印度洋非洲东海岸MAYG测站GPS数据进行实验，结果表明该方法相比二次多项式拟合方法，均方根误差降低了29%，相关系数为0.98，验证了该方法应用于海平面测高的有效性。

Abstract: For the problem of doping noise in Signal-to-Noise Ratio(SNR) signals existing in GNSS multipath reflectometry (GNSS-MR) technology, improved complete ensemble empirical mode decomposition with adaptive noise method is proposed to decompose the original Signal-to-noise ratio data, filter out the effective residual sequence without noise, implement denoising and applied to GNSS-MR technology to retrieve sea level height. Experiment with GPS data from MAYG station on the east coast of Africa in the western Indian Ocean, the results show that compared with the quadratic polynomial fitting method, the root mean square error has been reduced by 29%, the correlation coefficient is 0.98, the validity of the method applied to sea level altimetry is verified.

1. 引言

近年来，全球气候逐渐变暖，两极冰川融化，海平面也随之上升。海平面的上升严重影响了沿海地区的生态环境以及沿海居民的生命财产安全。而风暴潮和海啸等自然灾害也频繁发生，因此，监测海平面的高度变化对沿海地区和当地的居民生活具有重大的意义。传统的验潮站监测具有很高的测高精度，但其主要分布在陆海边缘，不利于远海区域的监测。卫星测高具有较高的测高精度，并且还有良好的覆盖能力，可以弥补传统验潮站在远海区域的缺陷。随着GNSS技术的不断发展，利用全球导航卫星系统多路径反射测量(GNSS-MR)技术进行海平面高度反演已逐步成为一种新的遥感方法，该方法具有全天候、实时自动化、成本低等优点，广泛应用于海平面高度反演 [1] [2] [3]，积雪深度反演 [4] 和土壤湿度反演 [5] 等方面。1993年，Martin-Neira等 [6] 首次提出PARIS (passive reflectometry and interferometry system)的概念，其主要思想是利用海面反射的GNSS信号来测量海平面高度，并利用多个卫星的信号来提高测量的精度。从PARIS的概念提出之后，国内外的许多学者采用二次多项式拟合的方法来获取SNR的残差序列，然后利用GNSS-R技术反演海平面高度 [7] [8] [9]。在此基础上，刘立龙等 [1] 从反演所需数据的选择方法的角度进行了研究，依据测站观测环境、有效反射区域反算的卫星高度角可以对反演所需的数据进行选择，结果表明海平面高度反演结果与验潮仪实测数据具有较强的相关性。陈发德等 [3] 为了探究多系统监测海平面高度的精度，利用GPS和BDS的SNR数据对海平面高度进行了反演，并与验潮站实测数据进行对比分析，结果表明GPS、BDS联合反演结果与验潮站实测高度的RMSE达到0.286 m。但上述研究采用的二次多项式拟合方法无法拟合出更为准确的SNR信号趋势，并且因为卫星高度角的增加，天线接收的有效反射信号受到干扰，使海平面高度反演结果出现较大的偏差。为了得到更为准确的SNR信号，王瑞芳 [10] 利用经验模态分解方法对SNR信号数据进行处理，剔除其中的噪声，得到了更接近验潮站实测值的海平面高度。但经验模态分解方法存在模态混叠的问题，在一定程度上影响GNSS-MR海平面测高的精度。针对此问题，Colominas M. A.等人提出改进的自适应噪声的完全集合经验模态分解(improved complete ensemble empirical mode decomposition with adaptive noise, ICEEMDAN) [11] 方法，该方法在原始SNR序列中添加经过EMD方法 [12] 分解之后对应模态阶数的高斯白噪声，能有效地解决EMD方法存在的模态混叠问题，防止其它信号的干扰。因此，本文将探讨利用ICEEMDAN方法代替二次多项式来获取SNR的残差序列，并以位于西印度洋非洲东海岸附近Dzaoudzi市的MAYG测站数据为例进行海平面高度反演，验证本文方法的可行性与有效性。

2. GNSS-MR监测海平面高度变化原理

卫星发射的直射信号和经海平面反射的反射信号同时被GNSS接收机接收，因为直射信号和反射信号两者的相位不同，所以会引起多路径效应。多路径效应主要记录在信噪比中，用来评估观测信号的质量。图1是利用GNSS-MR技术反演海平面高度的示意图，H为GNSS天线相位中心与验潮站高程基准(平均低低潮面，MLLW)间的垂直距离， $h_{\text{GNSS-MR}}$ 为GNSS-MR反演海平面高度， $\theta$ 为卫星高度角。

Figure 1. Retrieval principle of GNSS-MR sea level height

图1. GNSS-MR海平面高度反演原理

GNSS测量型接收机为了接收直射信号从而抑制反射信号，直射信号的振幅A_d与反射信号的振幅A_r存在如下关系

$A_d\gg A_r$ (1)

记合成信号的振幅为A_c， $\cos \phi$ 为直射信号与反射信号夹角的余弦值，其关系为

${\text{SNR}}^2=A_c^2=A_d^2+A_r^2+2A_d{A}_r\cos \phi$ (2)

因为直射信号的振幅A_d远远大于反射信号的振幅A_r的数值，所以合成信号总体的变化趋势由直射信号决定，为了获取低高度角下因多路径影响所生成的SNR残差序列，通常采用二次多项式拟合的方法来消除趋势项A_d。低高度角下因多路径效应所形成的SNR残差序列可表示为

$A_r=A\cos \left(\frac{4\pi h}{\lambda }\sin \theta +\phi \right)$ (3)

式中， $\lambda$ 为载波波长， $\theta$ 为卫星高度角，h为垂直反射距离，对SNR残差序列进行Lomb-Scargle频谱分析得到频率f,由 $h=f\lambda /2$ 可得天线相位中心到海平面的距离h，再把h转换为和验潮站相同基准下的海平面高度，从而实现GNSS-MR对海平面高度的反演。

3. ICEEMDAN的原理与方法

ICEEMDAN [13] [14] 是一种改进的自适应噪声的完全集合经验模态分解方法，该算法基于EEMD [15] [16] 和CEEMDAN [17] 算法做了进一步改进，通过实验得出ICEEMDAN克服了EEMD算法可能存在残余噪声的问题以及CEEMDAN算法可能存在伪模态的问题，能够有效地提高降噪的效果。

定义 $M(•)$ 为产生信号的局部均值， $E_k(•)$ 为运用EMD方法 [18] 分解所得到的第k个模态分量，设 ${\xi }^i$ 为高斯白噪声。ICEEMDAN所产生的第k个模态分量记为 $IM{F}_k$，ICEEMDAN算法的分解步骤如下 [19]：

1) 向原始序列添加一组高斯白噪声 ${\xi }^i$，构造序列 $x^i\left(n\right)=x\left(n\right)+{\epsilon }_0{E}_1\left({\xi }^i\left(n\right)\right)$，得到第一个残差：

$r_1=\left(M\left(x^i\left(n\right)\right)\right)$ (4)

其中 ${\epsilon }_0={\delta }_{0}std\left(x\right)/std\left(E_1\left({\xi }^i\left(n\right)\right)\right)$。

2) 计算可以得出第一个模态分量：

$IM{F}_1\left(n\right)=x\left(n\right)-r_1\left(n\right)$ (5)

3) 继续添加高斯白噪声，利用局部均值分解计算第二个残差 $r_1\left(n\right)+{\epsilon }_1{E}_2\left({\xi }^i\left(n\right)\right)$，并定义第二个模态分量：

$IM{F}_2\left(n\right)=r_1\left(n\right)-r_2\left(n\right)=r_1\left(n\right)-\left(M\left(r_1\left(n\right)+{\epsilon }_1{E}_2\left({\xi }^i\left(n\right)\right)\right)\right)$ (6)

4) 对于剩余阶段 $k=3,\cdots ,K$，计算出第k个残差：

$r_k\left(n\right)=\left(M\left(r_{K-1}\left(n\right)+{\epsilon }_{K-1}{E}_K\left({\xi }^i\left(n\right)\right)\right)\right)$ (7)

计算得到第k个模态分量：

$IM{F}_k\left(n\right)=r_{k-1}\left(n\right)-r_k\left(n\right)$ (8)

5) 重复步骤(4)，直至分解结束，得到所有的模态分量与残余分量。

4. ICEEMDAN在GNSS-MR技术反演海平面高度模型中的应用

本文将ICEEMDAN方法应用于GNSS-MR技术反演模型来反演海平面高度，其思路如下：首先从GNSS观测原数据中提取出SNR、高度角和方位角等数据，并设置合适的高度角和方位角区间筛选出在低高度角范围内的卫星SNR数据；然后通过ICEEMDAN方法分离直射分量和反射分量，获取反射信号的高精度SNR数据；最后，利用Lomb-Scargle谱分析法对SNR残差序列进行频谱分析得到频率，通过频率则可进一步计算出海平面高度。ICEEMDAN在GNSS-MR技术反演海平面高度模型中的应用的流程如图2所示。

5. 实验分析

本文选用MAYG观测站(http://www.igs.org)的观测数据进行实验分析，该测站位于西印度洋非洲东海岸附近的Dzaoudzi市(纬度−12.78˚，经度45.26˚)。图3为MAYG测站周围环境图，GNSS接收机为Trimble NETR9型，天线型号为TRM59800.00，数据采样间隔为30秒。位于距离MAYG测站2 m远的Dzaoudzi验潮站提供实测数据，可从REFMAR官网(http://refmar.shom.fr/)下载时间分辨率为1 min的验潮数据。

Figure 2. GNSS-MR retrieval sea level height process

图2. GNSS-MR反演海平面高度流程

Figure 3. Surrounding environment of MAYG station

图3. MAYG测站周围环境

本实验利用MAYG观测站2021年5月30日至2021年6月8日连续10天的GPS观测数据，以均误差(ME)、均方根误差(RMSE)和相关系数(R)三项精度指标来对比分析ICEEMDAN获得的有效残差序列和二次多项式拟合获得的有效残差序列反演海平面高度的精度：

$\{\begin{array}{l}\text{ME}=\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}\left(X_1\left(i\right)-X_2\left(i\right)\right)}}{n}\hfill \\ \text{RMSE}=\sqrt{\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(X_1\left(i\right)-X_2\left(i\right)\right)}^2}}{n}}\hfill \\ \text{R}=\frac{{\displaystyle \underset{i=1}{\overset{n}{\sum }}\left(X_1\left(i\right)-{\bar{X}}_1\left(i\right)\right)\left(X_2\left(i\right)-{\bar{X}}_2\left(i\right)\right)}}{\sqrt{{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(X_1\left(i\right)-{\bar{X}}_1\left(i\right)\right)}^2}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\left(X_2\left(i\right)-{\bar{X}}_2\left(i\right)\right)}^2}}}\hfill \end{array}$ (9)

其中，n为反演数据总数，X₁(i)、X₂(i)分别为验潮站实测数据和反演海平面高度。

本文选取2021年5月31日GPS卫星09号L2波段、卫星高度角为5˚~20˚、方位角为0˚~180˚和330˚~360˚的原始信噪比序列。经ICEEMDAN分解后结果如图4所示，横轴表示高度角正弦值，纵轴表示振幅。其中，IMF1为高频噪声信号，IMF2和IMF3为平稳变化的SNR信号，残余分量为SNR信号的趋势项，可以看出，IMF2和IMF3更为光滑，没有高频振荡部分。

Figure 4. Decomposition of SNR sequences by ICEEMDAN

图4. ICEEMDAN分解SNR序列

经过大量的实验，最终选取IMF2 + IMF3作为有效残差序列，并对其进行L-S谱分析得到频率f。图5是通过ICEEMDAN得到的信噪比残差序列的Lomb-Scargle谱分析图，横轴表示天线相位中心到海平面的距离，纵轴表示SNR反射信号频谱振幅，其中振幅的最高峰值所对应的就是天线相位中心到海平面的垂直反射距离。

本实验统计了GPS卫星在2021年5月30日至2021年6月8日的观测数据，如图6所示，验潮站的实测海平面高度变化因受到潮汐的影响具有很明显的日周期性；从图中可以看出，ICEEMDAN获得的有效残差序列和二次多项式拟合获得的有效残差序列两者的反演结果和验潮站的实测结果基本一致，验证了这两种方法反演海平面高度变化的可行性。图6为ICEEMDAN获得的有效残差序列和二次多项式拟合获得的有效残差序列反演的海平面高度结果对比图。

Figure 5. Lomb-Scargle spectral analysis

图5. Lomb-Scargle谱分析

Figure 6. Two methods to invert sea level and measured sea level

图6. 两种方法反演海平面高度与实测海平面高度

图7为ICEEMDAN获得的有效残差序列和二次多项式拟合获得的有效残差序列反演的海平面高度结果与验潮站实测结果的相关性分析，表1为ICEEMDAN获得的有效残差序列和二次多项式拟合获得的有效残差序列反演的海平面高度结果与验潮站实测结果的精度统计，结果表明利用ICEEMDAN方法很好地对原始SNR序列进行了去噪处理，避免了模态混叠，有效地防止了其他信号的干扰，进而得到与验潮站监测值更为接近的反演海平面高度。结合图6、图7、表1可知，通过ICEEMDAN获得的有效残差序列反演的海平面高度的均方根误差比二次多项式拟合降低了29%；其相关系数为0.98，比二次多项式拟合提高了0.01，表明ICEEMDAN方法获得的有效残差序列用于海平面高度反演具有良好的精度。

Figure 7. Correlation analysis between the inversion results of the two methods and the measured results

图7. 两种方法反演结果与实测结果相关性分析

Table 1. Accuracy statistics of sea level height and sea level height at tide gauge stations retrieved by quadratic polynomial fitting and ICEEMDAN

表1. 二次多项式拟合和ICEEMDAN反演的海平面高度与验潮站海平面高度的精度统计

6. 结语

本文利用MAYG测站的数据进行GNSS-MR海平面高度反演，通过利用ICEEMDAN对原始SNR序列进行分析，将有用的信号与噪声分离，提取出高质量的残差序列，进而利用GNSS-MR技术反演海平面高度。结果表明：

1) ICEEMDAN得到的有效残差序列能够准确地反映原始SNR信号的变化趋势，有效去除了原始SNR信号的噪声，从而得到精度更高的海平面高度。

2) ICEEMDAN方法能够更加精确地提取SNR直反信号的干涉信号，以验潮站实测高度为参考，ICEEMDAN方法获得的有效残差序列反演的海平面高度的RMSE为0.15 m，比二次多项式拟合方法降低了29%，验证了ICEEMDAN方法应用于海平面测高的有效性。

致谢

感谢国家自然科学基金提供资助，感谢刘老师和吴晗对论文提供指导，感谢文章里引用文献的所有者。

基金项目

国家自然科学基金42064002。

NOTES

^*通讯作者。