奇异谱分析在GNSS-MR海平面高度反演中的应用

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

无线电工程 2022年11期

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

(1.桂林理工大学测绘地理信息学院，广西桂林 541006；2.广西空间信息与测绘重点实验室，广西桂林 541006)

0 引言

随着全球气候变暖，南极和北极冰川融化速度加快，海平面逐年上升。为了更好地保护人类活动，特别是沿海地区的环境管理，监测海平面的变化至关重要。GNSS多路径反射测量(GNSS Multipath Reflectometry,GNSS-MR)作为海平面高度反演的一种新方法，尤其对于陆地边缘海平面高度的变化监测更为有利，该方法利用测量型接收机来接收GNSS直射信号和反射信号。直射信号和反射信号的相位不同引起多路径效应，相位的变化随着卫星的位置变化而变化，从而形成干涉图样，干涉图样以振幅的形式记录在信噪比(Signal to Noise Ratio,SNR)中，用于进行参数分析，具有全天候、实时自动化和成本低等优点，因此在海平面高度反演[1-3]、有效波高反演[4]、积雪深度反演[5]、土壤湿度反演[6-7]以及形变监测[8]等方面得到广泛应用。Larson等[9]利用GNSS接收机获取的GPS信号数据进行潮位反演，实验表明GPS-MR的反演结果与验潮站数据具有较好的一致性。吴继忠等[10]利用GPS观测数据对静止的湖面高度进行了反演，反演结果与实测结果比较吻合，GPS-MR反演的静止湖面高度标准偏差达到3 cm。张双成等[11]利用GPS观测数据进行海平面高度反演，结果表明GPS-MR反演结果与验潮站实测值相差约10 cm。刘立龙等[1]对反演所需数据的选择方法进行了研究，依据观测站观测环境、有效反射区域反算的卫星高度角可以对反演所需的数据进行选择，结果表明海平面高度反演结果与验潮仪实测数据具有较强的相关性。陈发德等[3]利用GPS和BDS的SNR数据对海平面高度进行了反演，并与验潮站实测数据进行对比分析，结果表明GPS，BDS联合反演结果与验潮站实测高度的均方根误差(Root Mean Square Error,RMSE)达到0.286 m。但上述对于海平面高度反演的研究均采用二次多项式拟合的方法获取SNR残差序列，无法拟合出更为准确的SNR信号趋势。为了获取更高质量的反射信号，苏晓容等[12]利用小波分析的方法研究GNSS-R技术监测潮位高度变化，并对复杂环境下采集的信号进行小波去噪，有效地提高了反演精度。王瑞芳[13]利用经验模态分解的方法对SNR残差序列进行分解，剔除其中的噪声信号，提取出较为纯净的海面反射信号，在一定程度上提高了反演精度。但是以上方法在处理SNR数据时存在一定的局限性，如经验模态分解的模态混叠问题，导致获取的SNR信号质量相对不高，包含噪声和干扰信息，影响GNSS-MR海平面测高的精度。

奇异谱分析(Singular Spectrum Analysis,SSA)是一种对时间序列进行主成分分析的方法，SSA分解的空间结构与时间尺度密切相关，在细节分量分析方面具有优势，它不受正弦波假定的约束，且无需先验信息即可从含有噪声的数据中提取出大量的可靠信息，并对其进行重构[14-15]。然而，鲜有研究SSA在海平面高度反演的应用。本文通过SSA法以布雷斯特港的BRST测站为例进行海平面高度反演，通过对比来验证本文方法的可行性与有效性。

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

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

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

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

Ad≫Ar,

(1)

记合成信号的振幅为Ac，cosφ为直射信号与反射信号夹角的余弦值，其关系为：

(2)

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

(3)

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

2 SSA的原理与方法

SSA是一种基于主成分分析的时间序列处理方法[16]，主要步骤为：构建轨迹矩阵、奇异值分解、分组和重构。具体过程如下[17-19]：

(1)构建轨迹矩阵

假设混有噪声的SNR序列Li的序列长度为N，对其进行滞后排列处理，并选择嵌入维度M，M为整数并且取值为1≤M≤N/2[20]，得到轨迹矩阵XM×K，其中K=N-M+1，轨迹矩阵X为：

(4)

(2)奇异值分解

X=X1+X2+…+Xd。

(5)

(3)分组

(4)重构

(6)

利用特征值变化率选取有用信号和噪声，将贡献率最大的特征分量作为其趋势项；特征值变化率突变拐点之后的有限个特征分量作为数据处理所需的残差序列；剩下的特征分量即为噪声。特征值变化率为:

r=(λk-λk+2)/2,k=1,2,…,M-2，

(7)

一般认为r≤0.001时符合实际情况，此时p=k-1。

3 实验分析

本文选用BRST观测站(http:∥www.igs.org)的观测数据进行实验分析，该测站位于法国西海岸的布雷斯特港(48°22′49.78″N，4°29′47.75″W)，该站点属于法国国家地理研究所。BRST测站周围环境如图2所示，GNSS接收机为天宝NETR9型，天线型号为TRM57971.00，数据采样间隔为30 s。位于BRST测站以北500 m的布雷斯特港验潮站提供实测数据，可从REFMAR官网下载时间分辨率为1 min的验潮数据。

图2 BRST测站周围环境

本实验利用BRST观测站2021年6月1日—15日连续15 d的GPS观测数据，以均误差(ME)、RMSE和相关系数(R)三项精度指标来对比分析SSA法获得的有效残差序列和二次多项式拟合获得的有效残差序列反演海平面高度的精度：

(8)

式中，n为反演数据总数；X1(i)，X2(i)分别为验潮站实测数据和反演海平面高度。

本文选取2021年6月15日GPS卫星08号L2波段、卫星高度角在5°～20°范围内、方位角在130°～270°范围内的原始SNR序列。奇异值贡献率如图3所示，当特征值个数为2时，奇异值贡献率变化速率突然降低，出现了明显的拐点，根据特征值变化率以及大量的实验，最终选择第2～17个特征分量重构有效残差序列。经SSA法分解和重构的结果如图4所示，横轴表示高度角正弦值，纵轴表示振幅。其中IMF1为选取的贡献率最大的特征分量所构成的趋势项；IMF2为通过特征值变化率选取的第2～17个特征分量所构成的有效残差序列；IMF3为选取的第18～25个特征分量所构成的噪声。

图3 奇异值贡献率

(a)原始SNR序列

将IMF2作为有效残差序列，并对其进行L-S谱分析得到频率f。图5是当天用SSA得到的残差序列。图6是通过SSA得到的SNR残差序列的Lomb-Scargle谱分析图，横轴表示天线相位中心到海平面的距离，纵轴表示SNR反射信号频谱振幅，其中,振幅的最高峰值所对应频率为149.530 4 Hz，和实际的天线相位中心到海平面的垂直反射距离18.26 m相符合。

图5 SSA获得的残差序列

图6 Lomb-Scargle 谱分析

本实验统计了GPS卫星在2021年6月1日—15日的观测数据，图7为SSA法获得的有效残差序列和二次多项式拟合获得的有效残差序列反演的海平面高度结果对比图。

图7 2种方法反演海平面高度与实测海平面高度

验潮站的实测海平面高度的变化有很明显的日周期性，主要是因为受潮汐影响。从图7中可以看出，SSA法获得的有效残差序列和二次多项式拟合获得的有效残差序列二者的反演结果和验潮站的实测结果基本一致，验证了这2种方法反演海平面高度变化的可行性，但仍有一些反演的海平面高度和验潮站的高度存在较大偏差，这主要是由波浪在某一时间段发生急剧变化所导致，同时这也是利用GNSS-MR技术反演海平面高度的结果与验潮站监测结果存在偏差的主要原因。为了进一步验证SSA法的有效性，实验选取在连续暴雨天气和外界干扰较严重的环境下，GPS卫星在2021年6月17日—21日的观测数据进行反演，反演结果如图8所示。

图8 连续暴雨天气下，2种方法反演海平面高度与实测海平面高度

表1为SSA法获得的有效残差序列和二次多项式拟合获得的有效残差序列反演的海平面高度结果与验潮站实测结果的精度统计。结果表明，利用SSA法很好地对原始SNR序列进行了去噪处理，进而得到与验潮站监测值更为接近的反演海平面高度。结合图7、图8和表1可知，通过SSA法获得的有效残差序列反演的海平面高度的RMSE比二次多项式拟合降低了6.5%，其相关系数为0.97，比二次多项式拟合提高了0.01；而在连续暴雨天气和外界干扰较严重的环境下，通过SSA法获得的有效残差序列反演的海平面高度的RMSE比二次多项式拟合降低了14%，其相关系数为0.98，表明SSA法获得的有效残差序列用于海平面高度反演具有良好的精度。