Research on BeiDou full-band GNSS-R water level inversion
-
摘要: 北斗三号卫星导航系统(BDS-3)拥有丰富的频段数据，使得该系统在全球卫星导航反射测量(GNSS-R)领域具有广阔的应用前景. 本文系统地开展了北斗全频段GNSS-R水位反演研究，首先对比了频谱分析和非线性拟合两种算法在水位反演中的性能，结果表明非线性拟合的整体精度优于频谱分析；然后针对反演结果中存在的各种误差，引入Savitzky-Golay滤波算法去噪，数值表明去噪前后的结果与实测水位之差的均方根误差(RMSE)平均降低了7.4 cm，且北斗全频段的反演精度与GPS L5相当.
-
关键词:
- GNSS-R /
- 水位反演 /
- 北斗卫星导航系统(BDS) /
- 频谱分析 /
- 非线性拟合 /
- Savitzky-Golay滤波
Abstract: China's BeiDou-3 (BDS-3) has abundant frequency band data, which makes the system promising for application in the field of Global Navigation Satellite System Reflectometry (GNSS-R). In this paper, we systematically carry out the study of BDS-R water level inversion, firstly, we compare the performance of two algorithms of spectral analysis and nonlinear fitting in water level inversion, and the results show that the overall accuracy of nonlinear fitting is better than that of spectral analysis; then, we introduce Savitzky-Golay filtering algorithm to denoise the inversion results for various errors, and the numerical values show that the difference between the results before and after denoising and the measured water level The average RMSE of the difference between the results before and after denoising and the measured water level is reduced by 7.4 cm, and the inversion accuracy of the full frequency band of BeiDou is comparable to that of GPS L5.-
Keywords:
- GNSS-R /
- water level inversion /
- BeiDou system /
- spectral analysis /
- nonlinear fitting /
- Savitzky-Golay filtering
-
0. 引 言
北斗卫星导航系统(BDS)作为中国自主建设、独立运行的卫星导航系统，与其他卫星系统相比，提供了更多的卫星信号，使得该系统具有全球卫星导航反射测量(GNSS-R)应用的巨大潜力. Jin等[1]首次利用BDS三频信号(B1I、B2I和B3I)的信噪比(SNR)数据与三频相位组合反演海平面高度，与验潮站结果相比，相关系数为0.87~0.91，均方根误差(RMSE)小于0.6 m. Wu等[2]研究了BDS中地球静止规道(GEO)卫星的水位反演性能，实验结果表明GEO卫星的时延多普勒和相位测高具有较高的反演精度，适用于监测海平面变化小的区域，且时间分辨率较高. 陈发德等[3]基于MAYG站的B1C、B2a和B2 SNR数据对水位进行了反演，实验结果表明BDS的三频数据反演精度略低于其他系统. Zheng等[4]分析了多频多模GNSS-R反演水位精度，结果表明北斗二号卫星导航系统(BDS-2) 地球同步轨道卫星（GEO）不适用于海岸水位反演，BDS-2中轨道卫星(MEO)的监测精度优于BDS-2倾斜地球同步轨道(IGSO)，北斗三号卫星导航系统(BDS-3)MEO的监测精度与BDS-2 MEO相当. Liu等[5]对北斗系统SNR反演水位高度的潜力进行了研究，但并未对B2b和B2a+B2b等新信号进行水位反演. 上述研究主要集中于如何提高反演精度和不同轨道高度卫星的测高性能，并未对北斗系统全频段的水位反演性能进行全面分析.
本文主要评估北斗系统全频段GNSS-R水位反演的性能. 首先描述GNSS-R水位反演方法，然后对比频谱分析与非线性拟合算法的反演算法精度，在选择合适的反演算法基础上，引入Savitzky-Golay滤波对北斗全频段反演结果进行滤波算法去噪，提升水位观测精度.
1. GNSS-R水位反演方法
1.1 GNSS-R水位反演原理
图1为 GNSS-R反演海平面高度原理图. 图中
$h$ 为天线相位中心(APC)到瞬时海平面的距离，θ为直射信号与瞬时海平面的夹角.为了有效地抑制测站周围反射物引起的多路径效应，直射信号和反射信号的振幅值需要满足如下的关系：
$$ {A_{\rm{d}}} \gt \gt {A_{\rm{m}}}. $$ (1) 干涉合成信号振幅的总体趋势取决于卫星发射的直射信号的振幅
${A_{\rm{d}}}$ ，而多路径引起的反射信号的振幅${A_{\rm{m}}}$ 则表现为低高度角下的周期性波动. SNR与干涉后合成信号振幅之间满足如下关系：$$ {{\rm{SNR}}^2} = A_{\rm{c}}^2 = A_{\rm{d}}^2 + A_{\rm{m}}^2 + {A_{\rm{d}}}{A_{\rm{m}}}\cos \; \psi .$$ (2) 式中：
${A_{\rm{c}}}$ 为GNSS测站天线接收的合成信号振幅；$\cos \; \psi$ 为直射信号与反射信号夹角的余弦值.结合式(1)可知
${A_{\rm{d}}}$ 与${A_{\rm{m}}}$ 在数值上相差较大，通常利用低阶多项式去除SNR的趋势项${A_{\rm{d}}}$ . 把图2(a)中的趋势项移除后，剩余项即可认为是受多路径影响所产生的SNR振荡序列，如图2(b)所示，可利用其反演反射体的环境参数. SNR振荡序列可表达为$$ \delta {\text{SNR}} = A\cos \left(\frac{{4\text{π} h}}{\lambda }\sin\;\theta + \varphi \right) .$$ (3) 式中：
$ A $ 为振幅值；$ h $ 为反演的距离；$ \lambda $ 为波长；$ \theta $ 为卫星高度角；$ \varphi $ 为延迟相位.由式(3)可知，水位反演就是依据
$ \delta {\text{SNR}} $ 观测值求解$h$ . Lomb-Scargle谱分析和非线性拟合是两种常用求解$h$ 的方法，下面分别给予描述.1.1.1 Lomb-Scargle谱分析方法
设
$t = \sin \;\theta$ ，$ f = {{2h} \mathord{\left/ {\vphantom {{2h} \lambda }} \right. } \lambda } $ ，则可将式(3)简化成标准的余弦函数形式：$$ \delta {\rm{SNR}} = A\cos (2 \text{π} ft + \varphi ). $$ (4) 通过式(4)可获得
$f$ ，则$$ h = \frac{{f\lambda }}{2}. $$ (5) 由式(5)中可知，只需要求解出SNR振荡序列的频率
$f$ 即可反算出高度$h$ .SNR振荡序列是以高度角正弦值为变量的不等间距信号，而且接收到的信号中存在噪声干扰. 在应用傅里叶变换时，这些噪声和不均匀分布会对频谱分析产生影响，导致频谱图与SNR信号的实际结果有很大偏差，甚至会出现虚假峰. Lomb-Scargle频谱分析方法不仅可以有效地获得弱周期信号的频谱估计值、解决非均匀采样信号的频谱估计等问题，还可以克服数据缺失和序列长度不足等困难. 下面简要阐述Lomb-Scargle频谱分析方法.
对于时序序列
$X({t}_{j})，j=1,2,3,\cdots,N$ ，其功率谱可定义为关于频率的函数：$$ \begin{split} {{{P}}_x}(f) =& \frac{1}{{2N}}\left\{ \frac{{\left[\displaystyle\sum\limits_{j = 1}^N {X({t_j})} \cos {{(2\text{π} f({t_j} - \tau ))}}\right]^2}}{{\displaystyle\sum\limits_{j = 1}^N {X({t_j})} {{\cos }^2}(2\text{π} f({t_j} - \tau ))}} +\right. \\& \left.\frac{{\left[\displaystyle\sum\limits_{j = 1}^N {X({t_j})} \sin {{(2\text{π} f({t_j} - \tau ))}}\right]}^2}{{\displaystyle\sum\limits_{j = 1}^N {X({t_j})} {{\sin }^2}(2\text{π} f({t_j} - \tau ))}}\right\} . \end{split} $$ (6) 式中：
$P{}_x(f)$ 为频率$f$ 的周期信号功率；$X({t_j})$ 为离散实验数据；${t_j}$ 为离散实验数据的时间；$N$ 为实验数据统计量；$\tau $ 为时间平移不变量.1.1.2 非线性拟合SNR信号方法
Strandberg 等[6]提出利用非线性最小二乘法求解海面高度，将式(3)添加一项衰减因子得到下式:
$$ \delta {\rm{SNR}} = A\cos \left(\frac{{4\text{π} h}}{\lambda }\sin\; \theta + \varphi \right)\exp ( - 4{k^2}\gamma {\sin ^2}\theta ). $$ (7) 式中，
$\gamma $ 与反射体的性质有关.为了保证有稳定的数值解，将式(7)转换为式(8)：
$$ \begin{split} \delta {\rm{SNR}} =& \left({C_1}\sin \left(\frac{{4\text{π} h}}{\lambda }\sin\; \theta \right) + {C_2}\cos \left(\frac{{4\text{π} h}}{\lambda }\sin \;\theta \right)\right) \\& \exp ( - 4{k^2}\gamma {\sin ^2}\theta ) . \\[-13pt]\end{split} $$ (8) 式中：
${C_1}$ 与${C_2}$ 共同决定式(9)中的振幅$A$ 与式(10)的相位$\varphi $ ：$$ A = \sqrt {\mathop C\nolimits_1^2 + \mathop C\nolimits_2^2 }, $$ (9) $$ \varphi = {{\rm{arctan}}}(\mathop C\nolimits_2^{} /\mathop C\nolimits_1^{} ). $$ (10) 同理，基于SNR振荡序列按照式(8)进行非线性拟合，即可求解反射高.
1.2 Savitzky-Golay滤波去噪
受测站环境因素影响，GNSS-R水位反演结果中不可避免会存在一些误差，导致反演精度下降. 本文拟通过Savitzky-Golay滤波算法[7]降低各种误差带来的影响. Savitzky-Golay 滤波是一种时间域范围内基于局部多项式的最小二乘拟合滤波方法，它能实现对数据特征的无衰减平滑.
下面对Savitzky-Golary滤波算法进行描述. 取某点附近的
$2m + 1$ 个连续数据进行$p$ 阶多项式拟合，其中$p \leqslant 2m + 1$ ，记该多项式为$$ {{F}}(t) = \sum\limits_{k = 0}^p {{b_k}} {t^k} = {b_{p0}} + {b_{p1}}t + {b_{p2}}{t^2} + \cdot \cdot \cdot + {b_{pp}}{t^p}. $$ (11) 式中：
$F(t)$ 为拟合后的多项式；bk为多项式系数；${t^p}$ 为多项式自变量的幂次方.最小二乘法通过最小化RMSE之和进行曲线拟合：
$$ S = \sum\limits_{i = - m}^m {(F(i)} - f(i){)^2}.$$ (12) 式中：
$S$ 为RMSE；$f(i)$ 为拟合前离散时间序列在第i点处的值.为了使该误差平方项达到最小值，令误差表达式
$S$ 对各系数求导数并使之为零，即：$$ \begin{gathered} \frac{{\partial {\text{S}}}}{{\partial {b_{pr}}}} = 2\sum\limits_{i = - m}^m \left({\sum\limits_{k = 0}^p {bpr} } {i^k} - f(i)\right){i^r} = 0, \\ r = 0,1,2, \cdot \cdot \cdot ,p. \\ \end{gathered} $$ (13) 从而有：
$$ \sum\limits_{k = 0}^p {bpr} \sum\limits_{i = - m}^m {{i^{k + r}}} = \sum\limits_{i = - m}^m {f(i){i^r}}. $$ (14) 因此，给定需要拟合的点数
$m$ 和多项式阶数$p$ ，求解方程组(14)，就可以求出系数序列$({b_{p0}},{b_{p1}},{b_{p2}}, \cdot \cdot \cdot , {b_{pp}})$ ，从而确定多项式$F(t)$ .使用Savitzky-Golay滤波算法对反演结果序列进行去噪时，根据反演结果序列的特征选择合适的多项式阶数及窗口长度. 本文选取窗口大小为21，曲线拟合阶数为3.
2. 实验及结果分析
2.1 实验数据
实验采用国际GNSS服务(IGS)SAS2站的北斗系统SNR数据进行水位反演，图3(a)为SAS2站点环境，图3(b)为该站方位角80°~360°、高度角5°~30°的菲涅尔区，可看出GNSS接收机距离水面较近，周边存在一些非水面反射. 为客观评价北斗全频段的反演效果，同时选取距站点约8 m的验潮站Sassnitz(http://www.ioc-sealevelmonitoring.org/station.php?code=sass)的实测水位作为验证数据.
2.2 频谱分析与非线性拟合方法比较分析
本实验分别使用频谱分析与非线性拟合对SAS2站2021年年积日1~365天的北斗全频段SNR数据求取反演结果并与该站GPS L5的反演结果作对比，研究表明GPS L5可取得较好的反演精度[8]. 通过将反演结果与验潮站实测水位进行对比分析以验证频谱分析与非线性拟合算法的准确性. 实验选择了RMSE、相关系数、反演点数量(所有SNR振荡信号拟合出的反演高度
$h$ 的数量)作为GNSS-R水位反演性能的指标.上述实验结果相关统计指标如表1所示. 可以看出：频谱分析方法中B2a频段反演效果相对于其他频段较差，B3频段反演效果最优；非线性拟合方法中B1频段反演效果相对于其他频段较差，B3反演效果最优；非线性拟合方法的反演结果都相差不大，其中RMSE最大与最小的差值为7 cm，而频谱分析方法RMSE最大与最小的差值达到了9.6 cm. 北斗全频段中通过非线性拟合的反演结果中，除B1C、B1、B2a的RMSE相对GPS L5的RMSE平均提高了0.87 cm，其余频段的RMSE都优于GPS L5频段. 而频谱分析方法的B2a和B2b的RMSE相对GPS L5的RMSE平均提高了2 cm. 虽然非线性拟合方法得到的B1C、B1、B3的反演精度略低于频谱分析，但其反演时间分辨率均高于频谱分析，且其余频段反演精度均高于频谱分析，尤其B2a频段的RMSE相对于频谱分析提高了17.05%.
表 1 北斗全频段水位反演结果统计信息方法 波段 相关系数 RMSE/cm 反演点数量 频谱分析 L5 0.726 14.7 6235 B1C 0.786 12.5 7568 B1 0.747 13.9 7886 B2a 0.659 17.6 9805 B3 0.889 8.0 11641 B2b 0.708 15.8 11213 B2a+B2b 0.728 14.5 7564 非线性拟合 L5 0.760 13.7 5983 B1C 0.746 13.9 7986 B1 0.717 15.2 8321 B2a 0.750 14.6 9907 B3 0.893 8.2 12940 B2b 0.778 13.3 11877 B2a+B2b 0.820 11.8 7167 2.3 水位反演结果Savitzky-Golay去噪
受测站环境因素的影响，水位反演值中不可避免地会带来误差，本文使用Savitzky-Golay去噪算法降低误差影响. 对2.2节中每个频段的非线性拟合反演结果进行去噪处理，去噪前后的反演结果如图4所示. 图中仅展示了SAS2第1天到第30天滤波前后水位数据与验潮站实测数据的对比，可以看出去噪后的反演结果明显与验潮站实测水位数据更加符合. 图5展示了滤波前后水位数据与实测数据的残差图，可以看出,滤波前的误差图中反演误差变化范围较大；而使用Savitzky-Golay去噪后，反演误差变化范围小于滤波前的误差变化范围.
去噪前后的结果与验潮站实测水位的统计结果如表2所示. 可以看出，经过Savitzky-Golay去噪处理后的水位反演精度均有很大提升，去噪后北斗全频段的RMSE相对去噪前的RMSE平均减少了7.4 cm，其中提升效果最显著的是B2a频段，其RMSE减少了9 cm，相关系数平均提升了16.40%.
表 2 GNSS-R滤波前后精度统计方法 波段 相关系数 RMSE/cm 滤波前 L5 0.760 13.7 B1C 0.746 13.9 B1 0.717 15.2 B2a 0.750 14.6 B3 0.893 8.2 B2b 0.778 13.3 B2a+B2b 0.820 11.8 滤波后 L5 0.936 6.0 B1C 0.934 6.0 B1 0.927 6.6 B2a 0.948 5.6 B3 0.968 4.3 B2b 0.950 5.4 B2a+B2b 0.961 4.8 为提高水位监测的时间分辨率和精度，考虑到北斗不同频率之间的互补性，将北斗系统全频段的反演水位组合在一起. 图6(a)是北斗全频段反演水位值组合，可以看出反演点数量明显增加，且时间分辨率相对于单一频率的时间分辨率更高，其RMSE为5.4 cm，优于B1C、B1以及B2a频段. 图6(b)为反演值与验潮站相关性分析，可以发现相关系数优于B1C、B1以及B2a等频段，与验潮站实测数据有良好的一致性. 总的来说，直接将北斗系统不同频率的反演结果组合，可达到提高精度和时间分辨率的目的.
3. 结束语
本文系统评估了我国北斗全频段SNR数据的GNSS-R水位反演性能. 首先比较了两种常用的反演水位算法的精度，实验表明非线性拟合方法综合效果优于频谱分析方法，虽个别频段反演精度略低于频谱分析方法，但非线性拟合方法的反演点数量多于频谱分析方法. 然后，针对反演结果中存在随机误差较大，本文通过Savitzky-Golay算法对反演的水位结果进行了滤波去噪，数值表明滤波后反演水位的RMSE极大降低，达到厘米级，且提升了与潮位数据的相关性. 整体上，北斗全频段的反演性能与GPS L5相当.
-
表 1 北斗全频段水位反演结果统计信息
方法 波段 相关系数 RMSE/cm 反演点数量 频谱分析 L5 0.726 14.7 6235 B1C 0.786 12.5 7568 B1 0.747 13.9 7886 B2a 0.659 17.6 9805 B3 0.889 8.0 11641 B2b 0.708 15.8 11213 B2a+B2b 0.728 14.5 7564 非线性拟合 L5 0.760 13.7 5983 B1C 0.746 13.9 7986 B1 0.717 15.2 8321 B2a 0.750 14.6 9907 B3 0.893 8.2 12940 B2b 0.778 13.3 11877 B2a+B2b 0.820 11.8 7167 
下载: 导出CSV
表 2 GNSS-R滤波前后精度统计
方法 波段 相关系数 RMSE/cm 滤波前 L5 0.760 13.7 B1C 0.746 13.9 B1 0.717 15.2 B2a 0.750 14.6 B3 0.893 8.2 B2b 0.778 13.3 B2a+B2b 0.820 11.8 滤波后 L5 0.936 6.0 B1C 0.934 6.0 B1 0.927 6.6 B2a 0.948 5.6 B3 0.968 4.3 B2b 0.950 5.4 B2a+B2b 0.961 4.8
下载: 导出CSV
-
[1] JIN S G, QIAN X D, WU X R. Sea level change from BeiDou Navigation Satellite System-Reflectometry (BDS-R): first results and evaluation[J]. Global and planetary change, 2017, 149: 20-25. DOI: 10.1016/j.gloplacha.2016.12.010
[2] WU J M, CHEN Y L, GAO F, et al. Sea surface height estimation by ground-based BDS GEO satellite reflectometry[J]. IEEE journal of selected topics in applied earth observations and remote sensing, 2020(13): 5550-5559. DOI: 10.1109/JSTARS.2020.3024743
[3] 陈发德, 刘立龙, 黄良珂, 等. 基于多模GNSS-MR海平面测高研究[J]. 地球物理学进展, 2018, 33(5): 1767-1772. [4] ZHENG N Q, CHEN P, LI Z. Accuracy analysis of ground-based GNSS-R sea level monitoring based on multi GNSS and multi SNR[J]. Advances in space research, 2021, 68(8): 1789-1801. DOI: 10.1016/j.asr.2021.04.024
[5] LIU Z Z, DU L, ZHOU P Y, et al. BDS/GNSS multipath reflectometry (BDS/GNSS-MR) based altimetry with new signals: initial assessment and comparison[J]. Advances in space research, 2022, 69(2): 282-291. DOI: 10.1016/j.asr.2021.08.025
[6] STRANDBERG J, HOBIGER T, HAAS R. Improving GNSS-R sea level determination through inverse modeling of SNR data[J]. Radio science, 2016, 51(8): 1286-1296. DOI: 10.1002/2016RS006057
[7] 马心宇, 张波, 王峰, 等. GPS L5反射信号测高性能研究[C]//第十届全国信号和智能信息处理与应用学术会议, 2016. [8] 王积丹, 匡翠林. Savitzy-Golay算法改善GNSS-MR土壤湿度反演研究[C]//第十一届中国卫星导航年会论文集——S01 卫星导航行业应用. 2020. -
期刊类型引用(1)
1. 李绪建. 基于北斗的水轮发电机维护远程采集系统研究. 北斗与空间信息应用技术. 2024(05): 25-28 . 
百度学术
其他类型引用(1)
计量
- 文章访问数: 347
- HTML全文浏览量: 92
- PDF下载量: 56
- 被引次数: 2
