利用全球导航卫星系统干涉反射(GNSS-IR)技术来监测海面高度变化是近年来新兴的监测方法^[1-3]。Löfgren等^[4]采用载波相位法反演海面高度，其精度可达到cm级；Jin等^[5]利用岸基BDS信号反演海面高度，其结果与验潮仪结果呈正相关；胡媛等^[6]使用信噪比数据进行海面高度的反演实验，实验结果与验潮仪测量结果相关度为86%。在GNSS-IR技术模型中，SNR数据处理方法主要为多项式拟合法去趋势和Lomb-Scargle谱分析法(LSP)提频相结合。本文提出一种基于小波变换和改进Burg算法的GNSS-IR海平面测高方法。小波变换中的小波函数可任意伸缩和平移，能较好地识别信号的局部特征，并且小波变换在对低频信号进行分离时不会丢失信号原有的高频信息，因此应用小波变换能有效去除信号中的趋势项。改进的Burg算法对短序列数据进行分析时不仅分辨率高，而且能有效抑制谱峰偏移或谱线分裂现象，提高谱估计精度。

1 GNSS-IR海面测高原理 1.1 经典GNSS-IR海面测高模型

如图 1所示，接收机可同时接收卫星直射信号和来自海面的反射信号，接收机天线距水面的反射高度为h，卫星仰角为θ。由几何关系可得，反射信号相对于直射信号存在额外路径，根据额外路径和仰角θ即可求得反射高度h^[7]。再结合反射信号的波长及相位，可将反射信号的额外路径转换为额外相位：

$ \varphi _{\rm{r}}^s = \frac{{4{\rm{ \mathsf{ π} }}h}}{\lambda }{\rm{sin}}\theta $

(1)

式中，λ为载波波长。接收机接收到的SNR信号可表示为：

$ \begin{array}{l} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{SNR}} = \left( {{A^S}} \right) = \\ \sqrt {{{(A_{\rm{d}}^S)}^2} + {{(A_{\rm{r}}^S)}^2} + 2A_{\rm{d}}^SA_{\rm{r}}^S{\rm{cos}}\left( {\frac{{4{\rm{ \mathsf{ π} }}h}}{\lambda }{\rm{sin}}\theta } \right)} \end{array} $

(2)

式中，(A_d^S)²+(A_r^S)²为导致信号向上升的趋势项(干扰项)。去除趋势项后的新信噪比值SNR_rv为干涉引起的SNR振荡项，可表示为：

$ \begin{array}{*{20}{l}} {{\rm{SN}}{{\rm{R}}_{{\rm{rv}}}} = 2A_{\rm{d}}^SA_{\rm{r}}^S{\rm{cos}}(\varphi _{\rm{r}}^S) = }\\ {\;\;\;2A_{\rm{d}}^SA_{\rm{r}}^S{\rm{cos}}\left( {\frac{{4{\rm{\pi }}h}}{\lambda }{\rm{sin}}\left( \theta \right)} \right)} \end{array} $

(3)

由式(3)可知，反射高度h与SNR_rv角频率有关，即

$ {\omega _{{\rm{SNR}}}} = \frac{{4{\rm{ \mathsf{ π} }}h}}{\lambda } = 2{\rm{ \mathsf{ π} }}{f_{{\rm{SNR}}}} $

(4)

式中，f_SNR为SNR振荡项的振荡频率，则海面高度可表示为：

$ h = {\rm{ }}\frac{{\lambda \cdot{f_{{\rm{SNR}}}}}}{{2}} $

(5)

1.2 小波变换原理

小波变换^[8]是利用一系列小波函数对原始信号进行逼近，分离出低频部分和高频部分。小波函数可任意伸缩和平移，能较好地识别信号的局部特征，并且小波变换在对低频信号进行分离时不会丢失信号原有的高频信息，因此应用小波变换能有效去除信号中的趋势项。由于小波函数以及小波分解层数的选择会决定SNR数据去趋势的效果，本文选择sym3小波函数和8层分解层数。

1.3 改进Burg算法

Burg算法^[9]是建立在线性预测原理基础上求解AR系数的有效算法，其特点是在Levison-Durbin递归约束条件下，求出前、后向预测误差功率之和最小时的反射系数，并据此计算出所有阶次的AR参数。

根据线性预测原理对信号x(n)某一时刻n前后的p个数值进行预测，得出前、后向预测误差e^f(n)、e^b(n)以及前、后向预测误差功率ρ^f、ρ^b：

$ \left\{ \begin{array}{l} {e^{\rm{f}}}\left( n \right) = x\left( n \right) - {{\hat x}^{\rm{f}}}\left( n \right)\\ {\rho ^{\rm{f}}} = E\left[ {|{e^{\rm{f}}}\left( n \right){|^2}} \right]\\ {e^{\rm{b}}}\left( n \right) = x\left( {n{\rm{ - }}p} \right){\rm{ - }}{{\hat x}^{\rm{b}}}\left( {n{\rm{ - }}p} \right)\\ {\rho ^{\rm{b}}} = E\left[ {|{e^{\rm{b}}}\left( n \right){|^2}} \right] \end{array} \right. $

(6)

根据Burg算法的特点，令前、后向预测误差功率之和最小，即

$ {\rho ^{{\rm{fb}}}} = \frac{1}{2}[{\rho ^{\rm{f}}} + {\rho ^{\rm{b}}}] $

(7)

ρ^f和ρ^b的求和范围不是0~(N-1+p)，而是p~(N-1)，即等效于e^f(n)和e^b(n)前后均不加窗，此时有：

$ \left\{ \begin{array}{l} \rho _p^{\rm{f}} = \frac{1}{{N - p}}\sum\limits_{n = p}^{N - 1} {} |e_p^{\rm{f}}\left( n \right){|^2}\\ \rho _p^{\rm{b}} = \frac{1}{{N - p}}\sum\limits_{n = p}^{N - 1} {} |e_p^{\rm{b}}\left( n \right){|^2} \end{array} \right. $

(8)

式中，当阶次m=1, 2, …, p时，e^f(n)和e^b(n)的递推关系可表示为：

$ \left\{ \begin{array}{l} e_m^{\rm{f}}\left( n \right) = e_{m - 1}^{\rm{f}}\left( n \right) + {k_m}e_{m - 1}^{\rm{b}}\left( {n{\rm{ - }}1} \right)\\ e_m^{\rm{b}}\left( n \right) = e_{m - 1}^{\rm{b}}\left( {n{\rm{ - }}1} \right) + {k_m}e_{m - 1}^{\rm{f}}\left( n \right)\\ e_0^{\rm{f}}\left( n \right) = e_0^{\rm{b}}\left( n \right) = x\left( n \right) \end{array} \right. $

(9)

可得使前、后向预测误差功率之和ρ^fb最小的反射系数$ {{\hat k}_m}$为：

$ {{\hat k}_m} = \frac{{ - 2\sum\limits_{n = m}^{N - 1} {} {\rm{ }}e_{m - 1}^{\rm{f}}\left( n \right)e_{m - 1}^{\rm{b}}\left( {n{\rm{ - }}1} \right)}}{{\sum\limits_{n = m}^{N - 1} {} |e_{m - 1}^{\rm{f}}\left( n \right){|^2} + \sum\limits_{n = m}^{N - 1} {} |e_{m - 1}^{\rm{b}}\left( {n - 1} \right){|^2}}} $

(10)

式中，m=1, 2, …, p。根据上式估计的$ {{\hat k}_m}$仍满足|$ {{\hat k}_m}$|＜1。按式(10)估计出$ {{\hat k}_m}$后，通过Levinson-Durbin算法递归获得m阶AR模型系数：

$ \left\{ \begin{array}{l} {{\hat a}_m}\left( k \right) = {{\hat a}_{m{\rm{ - }}1}}\left( k \right) + {{\hat k}_m}{{\hat a}_{m{\rm{ - }}1}}\left( {m{\rm{ - }}k} \right), \\ \;\;\;k = 1, 2, \cdots , m{\rm{ - }}1\\ {{\hat a}_m}\left( m \right) = {k_m}\\ {{\hat \rho }_m} = {\rm{ }}\left( {1{\rm{ - }}|{{\hat k}_m}{|^2}} \right){{\hat \rho }_{m{\rm{ - }}1}} \end{array} \right. $

(11)

由Burg算法的m阶AR模型系数可知，高阶系数由一阶系数递归获得，存在误差累积，即一阶系数的估计误差会对高阶系数的估计精度造成一定影响^[10]。对于给定的信号x(n)，由Yule-Walker方程和反射系数$ {{\hat k}_m}$可计算得到一阶模型参数a₁(1)的真实值和估计值：

$ \left\{ \begin{array}{l} {a_1}\left( 1 \right) = {\rm{ - cos}}\omega \\ {a_1}\left( 1 \right)^{\prime} = {\rm{ - cos}}\omega + \\ \;\;\;\;\;\;\frac{{{\rm{cos}}\left( {\omega n + \theta } \right){\rm{sin}}\left( {N\omega } \right){\rm{si}}{{\rm{n}}^2}\left( \omega \right)}}{{\left( {N - 1} \right) - {\rm{sin}}\left( {2\omega } \right){\rm{cos}}\left( {N\omega + 2\theta } \right){\rm{sin}}\left( {N\omega } \right)/2}} \end{array} \right. $

(12)

从上式可以看出，与a₁(1)真实值相比，其估计值存在估计误差项。因此，鉴于其误差源，需要对一阶系数误差项进行修正，即改进Burg算法。改进的算法是在预测误差功率最小的条件下直接求解模型的高阶系数。同时，为了不增加计算量，可以先计算二阶预测误差功率^[11]。对于给定的信号x(n)，可以先求得二阶误差平均功率ρ₂，再计算二阶模型系数a₂(1)、a₂(2)，最后得到a₁(1)真实值为：

$ {a_1}\left( 1 \right) = {\rm{ }}\frac{{{a_2}\left( 1 \right)}}{{1 + {a_2}\left( 2 \right)}}{\rm{ }} = - {\rm{cos}}\omega $

(13)

可以看出，通过二阶系数的推导间接校正一阶系数的估计误差，可避免误差项的影响，从而极大提高谱估计精度。

2 GNSS-IR海面测高模型验证实验

本文使用瑞典Onsala空间天文台GTGU站(57.392 954 9°N、11.913 488 6°E)的GNSS数据(https://zenodo.org/record/2924309)进行实验^[12]。GTGU站安装在平均海平面上方约4 m处，接收机采样频率为1 s。在距接收机约1 km处安装验潮仪测量实际海面高度变化，并与反演结果进行对比分析。本次研究主要针对SNR去趋势项以及功率谱分析，并对GNSS-IR海面测高模型进行实验验证。

实验所选的GTGU站因靠近北极圈，其潮涨潮落周期无规律。图 2为验潮仪记录的该站点2016-01~2016-06海面高度变化情况，虽然整体波动较小，但十分不规则，尤其是局部突变较多。

2.1 基于小波变换的GNSS-IR海面测高模型精度分析

本文使用DOY32~38的SNR观测数据，结合LSP频谱分析法分别应用小波变换算法和多项式拟合法对海面高度进行反演，结果见表 1。可以看出，基于小波变换的海面测高方法反演误差为4.96 cm，与验潮仪的相关系数达到0.98，且反演的有效点数为344个，能够较好地实现海面高度反演。小波变换算法的3项指标均优于多项式拟合法，在精度上可提高约20%，表明优化后的GNSS-IR海面测高模型具有可行性，同时说明新模型可提高GNSS数据的利用率。

2.2 基于小波变换和改进Burg算法的GNSS-IR海面测高模型精度分析

为了对比LSP和改进的Burg算法对反演结果的影响，选取2016-05的SNR观测值，基于多项式拟合法对海面高度进行反演，结果见表 2。可以看出，基于改进Burg算法的反演结果在精度和有效反演点数上均优于LSP，说明改进的Burg算法在提取SNR振荡频率时具有一定优势。

本文在小波变换基础上，利用改进的Burg算法进行频谱分析来反演海平面高度。为分析仰角范围对反演结果的影响，设置5°~20°、5°~25°及5°~30°三个仰角范围进行实验，结果见表 3。图 3为5°~30°卫星仰角范围的反演结果。可以看出，随着仰角升高，反演误差逐渐降低，反演的有效点数逐渐增加，同时相关系数也有所升高，表明新模型可提高GNSS数据的利用率，并且具有良好的稳定性。

为进一步验证小波变换对反演结果的改善效果，表 4为结合LSP谱分析进行反演得到的最终结果。可以看出，改进的Burg算法反演海面高度误差最小，同时相关度最高。一般情况下将误差作为首要考虑因素，而本文提出的改进Burg算法反演海面高度误差最小，与验潮仪数据的相关度也最高，相比于LSP误差至少降低20%。本文还应用GTGU站2016-01~2016-06(其中DOY85~88数据缺失)的数据验证新模型在长时间海面高度反演中的效果，结果见表 5。从表 4和表 5可见，Burg算法的反演结果差于LSP(长期结果相差不大)，其原因是Burg算法存在谱峰偏移和谱线分裂等异常现象，使谱估计结果出现较大偏差。而改进的Burg算法在改进后可去除一阶参数产生的累积误差，使得反演结果的精度优于其他方法，且反演误差至少降低20%，相关系数能达到0.95(短期结果为0.98)，反演的有效点数也相对较多。因此本文提出的小波变换算法和改进Burg算法的谱分析法能有效反演海面高度，而且适用于类似GTGU站所处的海平面变化无规律性的海域，具有精度高、相关度好、有效点数多等优点。

3 结语

本文提出结合小波变换和改进的Burg算法反演海面高度，构建一种新的GNSS-IR海面测高模型，并在GTGU站进行海面高度反演实验，同时与多项式拟合法、LSP等频谱分析方法的反演结果进行对比。结果表明，短期时间序列的反演结果与验潮仪站测量值的相关度达到98%，反演精度为4.48 cm；长期时间序列反演结果的相关度也达到95%，精度为5.45 cm，表明新模型在海面测高方面具有可靠性和有效性，相比于传统测高模型精度可提高约20%，在稳定性和GNSS数据利用率上也有较大提高。同时，本文提出的新型GNSS-IR海面测高模型在反演海面高度变化无规律的海域时也具有独特优势。