基于连续小波变换及正则化方法的潮位分析

姜桂婷， 李水艳

(河海大学 理学院， 江苏 南京 211100)

0 引 言

潮位的精确检测对于我国海洋事业的发展起着不可或缺的作用. 目前， 有成千上万的全球定位导航系统(Global Navigation Satellite System， GNSS)[1]正在全世界运行， 多路径原理[2]即信号中同时包含直射信号以及周围地表的反射信号信息， 从而利用现成的大量商用GNSS提供其周围环境的几何信息[3-4].

GNSS干涉遥感(GNSS-Interferometr Reflectomety, GNSS-IR)技术基于常规的大地测量型接收机， 利用SNR序列中的干涉振荡特性， 便可完成对站点周围环境的监测[5-7]. Larson等[8]首次将地基(岸基)GPS 的信噪比(Signal-to-Noise Ratio， SNR) 序列用于潮位监测中. 此后， Larson等[9]在 Kachemak 湾利用 PBAY 站的 SNR 信号进行潮位反演， 并且提出改正动态海面偏差的方法. 随后， Löfgren等[10]又进一步研究改进此方法.

目前， GNSS-IR 潮位反演技术已经有一套较为完善的经典潮位反演原理及技术流程： 即利用 LSP(Lomb-Scargle Periodogram)提取 SNR 中的频率信息， 并将该信息转换为天线相位中心与反射面之间的垂直距离(Reflector Height， RH)， 简称为有效高度， 再将其统一到潮位基准上. LSP 方法对于一段输入的 SNR 序列只能输出一个频率值. 海面的时刻变化特性表明对于具有一定时间长度的SNR来说， 与海面高度相关的频率也随之变化， 由此可知LSP方法的不合理性， 且其损失了反演点的分辨率. 针对这个问题， 前人提出Window Lomb-Scargle谱分析(winLSP)方法， 通过滑动窗口， 用LSP方法提取由窗口调整的每个段的频率. 该方法虽然能够输出更多的反演点值， 但其结果精度远不及 LSP方法.

小波变换常作为数据信噪比分析研究的一种时频分析方法. Bilich等[11]利用小波分析提取的SNR序列的光谱含量来确定造成显著多路径错误的因素. Wang 等[12-13]利用连续小波变换来提取多模多频 SNR 序列中的瞬时频率， 以期望提高 GNSS-IR 潮位反演中的数据分辨率， 与此同时证明了该频率在潮位检索中的潜力.

对于多模多频站点潮位分析， 利用连续小波变换可以得到不同SNR序列的小波谱图. 本文基于系统不同质量的SNR序列分析， 给出一定的正则化降噪模型， 从而获得一个高质量的RH序列来进行潮位分析.

1 背景知识

1.1 GNSS-IR分析原理

信号接收器接收来自卫星的直射信号与反射信号. 在只有一次多路径反射时， 与直射信号相比， 反射信号传播了额外的距离D=D1-D2， 如图 1 所示.

图1 一次多路径反射原理图Fig.1 Single multipath reflection schematic diagram

由图 1 中的几何关系容易得到

(1)

根据式(2)计算出直射信号与反射信号之间的相位差φ

(2)

式中：λ为信号的波长.

在静态或准静态的情况下， 即不考虑海面和卫星的运动时， 有

(3)

式中：f为正弦线SNR序列的频率； d表示微分符号.

由此得

(4)

考虑到海面时刻都在变化的特性， 特别是在潮波变化较大的地方， 静态假设变得不合理. 高度速率误差可以通过增加改变海面高度的一阶修正项来修正， 即

修正后的RH时间序列h为

(5)

最终海平面高度为

s=C-h(t),

(6)

式中：s为海平面高度；C为恒定常数， 与站高、 GNSS系统的框架基准和验潮站的高度基准有关.

1.2 连续小波变换

与傅里叶变换一样， 连续小波变换使用内积来度量信号f和分析函数之间的相似性， 然而连续小波变换不再使用三角函数基， 而是使用随时间急速衰减的小波基ψ， 通过平移和伸缩， 得到小波系数矩阵

(7)

式中：a为尺度参数， 控制小波的伸缩；b为平移参数， 控制小波的平移.

连续小波变换同时从时域和频域来分析信号， 其频率信息包含在尺度中， 尺度减小， 中心频率提高， 反之亦然. 通过连续小波变换可以了解信号的频率成分， 与此同时还能确定它在时域上存在的具体位置. Abry在1997年给出了在给定尺度上具有指定采样周期的小波与频率之间的映射.

2 数据处理过程

2.1 连续小波变换提取瞬时频率

对于δS(sine)序列， 其连续小波变换为

W(a,b；δS(sine),ψ(sine))=

(8)

式中： *代表复共轭；δS(sine)为通过二次拟合去除直射信号后的SNR序列.

利用a和b的值， 可以得到小波变换系数矩阵W(a,b). 同时， MATLAB小波工具箱软件提供了2个函数： centfrq 和scal2frq， 给出了指定尺度小波和实际频率的近似关系. 小波分析频谱中， 每个历元对应多个频率(即有效高度)， 本文选择每一历元能量最大值对应的RH值， 认为该 RH 值是该历元的最显著频率对应的RH值.

2.2 频率域加权

基于期望得到的新的RH序列具有一定的连续性， 可通过计算每个RH序列异常点的个数来确定其加权的权重， 异常点越多的序列所占权重越小.

若得到的RH序列分别为ri(i=1,2…)， 称满足条件|r(i+1)-r(i)|>0.5的点为异常点， 其中r(i)代表r序列的第i个元素.

设其对应异常点个数为si(i=1,2…)， 选取权重为

(9)

根据经验可得， 频域率的加权比时间域的加权更加适合， 因此， 进一步加权从SNR序列中得到的小波系数矩阵Wi(i=1,2，…)， 可以得到

(10)

取W每列小波系数最大处对应的RH值， 得到所需要的新r序列.

2.3 正则化去噪

正则化方法最常用的一个应用领域就是去噪. 正则化过程是在已知r序列的基础上寻找一个更好的估计. 由于希望能在保证连续平滑性的基础上减少噪音， 因此， 得到的最小二乘问题为

(11)

则此优化问题的解为

(12)

3 数值实验与比较

3.1 数据来源

站点选取了能够接收到多模多频GNSS信号的Multi-GNSS Experiment (MGEX)站： MAYG站.

MAYG站位于印度洋附近的马约特岛(12.8°S, 45.3°E)， 安装Trimble TRM59800.00天线(不带天线罩)， 该天线与Trimble NETR9 GNSS 接收机连接， 记录采样间隔30 s的观测数据. 该站可以接收到多模多频信号. 在位于MAYG 站几米处政府间海洋学委员会(Intergovernmental Oceanographic Commission) 建立了Dzaoudzi验潮站， 该潮位站以1 min的时间分辨率收集潮位数据， 以此得到了实测潮位. MAYG站的海域方位角为20°～170°.

本文选取MAYG站2017年第250年积日GPS PRN 07, GLONASS PRN 23, Galileo PRN 09卫星的SNR 数据进行分析. 具体使用的SNR类型如表 1 所示.

表 1 SNR数据类型Tab.1 SNR data type

表 1 中， 对于GPS系统， 使用3种SNR 数据类型： S1C， S2W， S2X； 对于GLONASS系统， 使用3种SNR数据类型： S1C， S1P， S2C； 对于Galileo系统， 使用3种SNR数据类型： S1X， S5X， S7X.

3.2 数据实验

首先对GPS系统3种SNR序列进行实验研究， 3种SNR序列如图 2 所示.

利用连续小波变换提取SNR序列的瞬时频率， 得到如图 3 所示的小波频谱图.

(a) S1C

(b) S2W

(c)S2X图3 3种信号对应的小波频谱图Fig.3 Wavelet spectrum corresponding to three signals

图 3(a) 是由GPS系统S1C类型SNR序列得到的小波时频图， 图 3(b) 是由GPS系统S2W类型SNR序列得到的小波时频图， 图 3(c)是由GPS系统S2X类型SNR序列得到的小波时频图. 小波谱图不仅给出了每个历元不同RH(或频率)的功率(或幅度)， 而且还反映了SNR序列的质量. 质量较好的SNR系列具有边界清晰的功率集中区域和较少的由噪声引起的散乱区域. 根据图中小波分析谱图的表现， 可以看出每个SNR类型的质量等级为： S2X> S2W> S1W.

对于3幅小波时频图选取最大能量处对应的RH值， 得到如图 4 所示的RH序列.

应用式(9)和式(10)得到r序列， 需要先对序列的异常点进行基本预处理

r(i+1)=1.001*r(i).

(a) S1C

(b) S2W

(c) S2X图4 3种信号对应的RH序列Fig.4 RH sequence corresponding to three signals

对于得到的新r序列， 应用模型(11)进行降噪， 得到如图 5 所示的序列， 其中， 正则化参数矩阵对角元素为

图5 正则化降噪后的RH序列Fig.5 RH sequence after regularization and noise reduction

3.3 潮位对比

GPS系统3种不同类型RH序列通过式(5)和式(6)得到的潮位数据与实测潮位对比如图 6 所示.

(a) L1C/A 类型

(b) L2W类型

(c) L2C类型图6 未采用正则化去噪处理得到的潮位数据与实测潮位的对比Fig.6 Comparison between tidal level data obtained withoutregularization denoising and measured tidal level

图 6 中黑色点代表实测的潮位数据， 灰色点代表由GPS不同类型SNR数据推导出的潮位.

降噪处理后推导得到的潮位与实测潮位的对比如图 7 所示.

通过对比图 6 与图 7， 能直观地看到降噪处理后推导出的潮位与实际潮位更加接近. 为了使对比更具统计学的意义， 比较其均方根误差(Root Mean Square Error， RMSE)， 该方法衡量了推导值与记录值之间的误差， 数值越小说明误差越小.

同时本文也用此方法对GLONASS和Galileo系统的不同类型的SNR数据进行了潮位分析， SNR数据类型与表1相对应， 由于篇幅有限， 就不展开叙述， 得到具体结果的均方根误差如表 2 所示.

图7 降噪处理后的潮位与实测潮位对比Fig.7 Comparison between the tide level obtained bynoise reduction and the measured tide level

表 2 所示为GPS系统第7号卫星， GLONASS第23号卫星和Galileo第9号卫星对应3种SNR数据潮位分析结果的RMSE值， 以及加权降噪后的RMSE值， 数值对比发现RMSE值有明显的降低. 由此可知本文方法利用同系统不同类型SNR序列的信息， 在提高数据利用率的同时降低了误差.

表 2 均方根误差比较Tab.2 Comparison of RMSE

4 结 论

本文首先利用连续小波变换来提取SNR序列的瞬时频率， 由于实际卫星接收到的信号噪音很大， 因此得到的潮位数据误差较大. 然后利用同系统不同类型SNR序列的综合信息进行正则化的降噪处理， 与正则化方法处理前的SNR序列进行的潮位数据相比， 至少能提高20%左右的精度， 具有一定的应用和推广意义.