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

2022-07-11姜桂婷李水艳

测试技术学报 2022年4期

姜桂婷，李水艳

(河海大学理学院，江苏南京 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).