CN114384560A

CN114384560A - 一种基于tvf-emd-svd的gnss多径信号抑制方法

Info

Publication number: CN114384560A
Application number: CN202210067628.8A
Authority: CN
Inventors: 陈熙源; 于凌宇; 闫晣
Original assignee: Southeast University
Current assignee: Southeast University
Priority date: 2022-01-20
Filing date: 2022-01-20
Publication date: 2022-04-22

Abstract

一种基于TVF－EMD－SVD的GNSS多径信号抑制方法，首先将含有多径效应的GNSS信号进行时变滤波经验模态分解，分解得到m个本征模态函数分量，多径信号往往集中在低频的本征模态函数分量中，利用分解得到的m个本征模态函数分量构造时频矩阵，通过对IMF分量构建的时频矩阵进行奇异值分解，将多径信号的各频率成分分解到不同的时频子空间中，然后，计算各时频子空间自相关函数与原信号自相关函数的相关系数，选择相关系数大于等于0.5的时频子空间进行信号重构，从而实现对多径信号的抑制。本发明能够避免EMD在分解时出现的模态混叠情况，同时对于SVD降噪时不能准确选取重构阶次的问题，引入自相关函数和自相关系数，来帮助SVD进行准确的重构阶次选择。

Description

一种基于TVF-EMD-SVD的GNSS多径信号抑制方法

技术领域

本发明属于卫星导航领域，涉及GNSS多径信号抑制方法，具体地，具体为一种基于TVF-EMD-SVD的GNSS多径信号抑制方法。

背景技术

随着全球导航卫星系统(Global Navigation Satellite System,GNSS)的建成和应用，卫星导航开启了人类导航的新时代，同时人们对其定位精度的要求越来越高。然而，卫星信号在传播的过程中会经过障碍物产生反射或折射，导致接收机天线除了接收到一个从卫星发射后经直线传播的电磁波信号之外，还可能接收到一个或者多个由该电磁波经周围地物反射后的信号，而每个反射信号又可能是经过一次或者多次反射后到达天线的，进而引起多径误差。因此，多径效应成为制约定位精度的主要误差源之一。

目前，针对GNSS多径信号的抑制方法主要有三大类：一是基于接收机天线的设计；二是基于接收机跟踪环路的改进；三是基于数据后处理的多径抑制。前两者在实现上存在一定的困难和缺点：基于天线端的设计，如：采用扼流圈、光子晶体天线等，虽然可以从接收端的源头抑制多径信号，但是其体积大、成本高；基于接收机跟踪环路的多径抑制技术，如：多径估计延迟锁定环(Multipath Estimation Delay Lock Loop,MEDLL)、Strobe相关器、PAC技术等，通过增加相关器的个数来实验多径信号抑制，但是其增加了接收机内部的计算量，使得接收机的成本变高。近年来，随着信号处理技术的不断发展，基于数据处理的各种滤波降噪算法正成为多径抑制的主要方法。经验模态分解(Empirical ModeDecomposition,EMD)是一种数据驱动的自适应时域分解方法，有很强的局部自适应分解特性，适用于非线性、非平稳信号的处理，但是该方法在分解的过程中会引起幅值的改变，通常重构存在误差。针对GNSS多径信号低频的特点，传统的EMD分解后得到的本征模态函数(Intrinsic Mode Function,IMF)在低频段存在少量的模态混叠，导致重构后的信号依旧含有少量多径信号。

发明内容

为解决上述技术问题，本发明提出了一种基于TVF-EMD-SVD的GNSS多径信号抑制方法。

一种基于TVF-EMD-SVD的GNSS多径信号抑制方法，具体步骤如下：

步骤1：选取含有多径效应的GNSS信号作为输入信号s(t)；

步骤2：计算局部截止频率；

首先，通过Hilbert变换得到信号s(t)的瞬时幅值A(t)和瞬时频率

确定A(t)的局部最大值A({t_max})与局部最小值A({t_min})，以及

的局部最大值

与局部最小值

然后进行差值运算得到η₁(t)和η₂(t)，同时，利用时变滤波器对输入信号进行滤波，获得瞬时均值a₁(t)及瞬时包络a₂(t)，进而计算局部截止频率

步骤3：依据调整后的截止频率将信号重构为h(t)，令

步骤4：判断是否满足停止准则，计算截止准则θ(t)

其中

为加权均值瞬时频率；B_L(t)为Loughlin瞬时带宽；

为第i阶分量的瞬时频率；

判断是否满足停止准则，如果θ(t)≤ξ，其中ξ为给定带宽阈值，则判定s(t)为一个IMF，否则，将h(t)的极值点作为节点，将h(t)分成n段，每段步长为m，采用B样条插值对信号s(t)进行逼近，得到逼近结果为m(t)，即令s₁(t)＝s(t)-m(t)，重复执行上述步骤(2)～步骤(4)；

对信号进行TVF-EMD分解后，得到一组IMF分量imf_i(t)，即

其中，IMF指本征模态函数，i为模态分量阶数，imf_i(t)为第i阶模态分量序列，N为分解得到的IMF总数；

步骤5：利用得到的IMF构造时频矩阵P；

步骤6：对时频矩阵进行奇异值分解；

步骤7：计算时频子空间与原始信号的相关函数和相关系数；

步骤8：选择相关系数大于等于0.5的时频子空间进行重构，并对选取的子空间进行奇异值分解逆变换；

步骤9：将重构的信号转化为一维信号，即得到的信号为多径抑制后的信号。

作为本发明进一步改进，所述步骤1的多径信号s(t)应满足：

其中，A为信号幅值，p(t)是值为±1的数据码与伪码的异或和，f是考虑多普勒效应后的载波频率，α_i是反射波的衰减系数，τ_i是反射波i相对于直射波的传播延时，

是信号在各个反射面反射前后的相位变化。

作为本发明进一步改进，所述步骤4中的带宽阈值ξ＝0.08，B样条阶数为20。

作为本发明进一步改进，所述步骤6奇异值分解具体操作步骤如下：

(1)对时频矩阵P进行奇异值分解，得到下式所示的分解形式

P＝UQV^T

其中，U是一个m×m阶的左正交矩阵，V是一个n×n阶的右正交矩阵，Q是一个对角矩阵，具体表示为：

在上式中，σ_i(i＝1,2,3,…,q)为时频矩阵P的奇异值，并且σ₁＞σ₂＞σ₃＞...＞σ_q＞0，q是通过分解获得的奇异值的总数

(2)：根据时频子空间计算得到的相关系数，选择相关系数大于等于0.5的时频子空间，判断出有效奇异值的个数为k。

(3)：保留对角矩阵Ω的前面k个有效奇异值，然后将剩下的其他较小的奇异值做置0处理，以获得新的奇异矩阵Q′，将新的奇异矩阵Q′，左正交矩阵U_m×m和右正交矩阵V_n×n代入，运算得到新的矩阵P′。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：

1.TVF-EMD是完全自适应的，适用于线性和非平稳信号的分析。与EMD相比，所提出的方法能够提高频率分离性能，以及低采样率下的稳定性。局部截止频率是通过充分利用瞬时幅度和频率信息自适应设计的。然后采用非均匀B样条近似作为时变滤波器。为了解决间歇性问题，还引入了截止频率重排算法。为了提高低采样率下的性能，提出了固有模式函数的带宽准则。此外，所提出的方法对噪声干扰具有鲁棒性。

2.传统的SVD降噪算法中，所保留的某些奇异值以噪声贡献为主，不能准确选取奇异值的个数，由此得到的纯信号矩阵估计含有的噪声也较多。本方法中含有多径信号的时频特征信息被分解到一系列的由U和V构成的奇异向量时频子空间中，且时域信息通过V表达，频域信息通过U表达，因此选取左右奇异量信号的有效时域成分来重构矩阵的阶次，此方法较传统的SVD降噪，能够更多的保留有用信号。

附图说明

图1本发明的算法流程图；

图2多径信号波形图；

图3多径信号经TVF-EMD分解后的分解结果；

图4时频子空间自相关函数；

图5抑制后的信号。

具体实施方式

下面结合附图与具体实施方式对本发明作进一步详细描述：

如图1所示，本发明公开了一种基于TVF-EMD-SVD的GNSS多径信号抑制方法，首先将含有多径效应的GNSS信号进行TVF-EMD分解，分解得到m个本征模态函数分量，多径信号往往集中在低频的本征模态函数分量中，利用分解得到的m个本征模态函数分量构造时频矩阵，通过对IMF分量构建的时频矩阵进行奇异值分解，将多径信号的频率信息分解到不同的时频子空间中，然后，计算各时频子空间自相关函数与原信号自相关函数的相关系数，根据相关系数选取合适的时频子空间进行信号重构，从而实现对多径信号的抑制。

步骤1：输入含有低频多径信号的GNSS信号s(t)。

是信号在各个反射面反射前后的相位变化。

步骤2：利用希尔伯特变换计算出多径信号s(t)的瞬时幅值A(t)和瞬时频率

如式(2)和(3)：

式中，

为信号的希尔伯特变换。

步骤3：确定A(t)的局部最大值A({t_max})与局部最小值A({t_min})，同时确定

的局部最大值

与局部最小值

步骤4：对A(t)的极值点分别进行插值，得到曲线分别为β₁(t)和β₂(t)。

步骤5：计算出瞬时均值a₁(t)和瞬时包络a₂(t)，其中：

a₁(t)＝[β₁(t)+β₂(t)]/2 (4)

a₂(t)＝[β₁(t)-β₂(t)]/2 (5)

步骤6：对

A²({t_min})和

A²({t_max})分别进行插值，得到η₁(t)和η₂(t)，并根据式(6)和(7)计算得到

和

步骤7：计算局部截止频率

步骤8：根据调整后的截止频率重构信号，如式(9)所示：

将h(t)的极值点作为节点，将h(t)分成n段，每段步长为m，其中n称为B样条的阶数。使用式(10)进行B样条插值逼近，逼近结果记为m(t)，代表曲线的均值。

式中：[·]为下采样操作，

为预滤波器且

βⁿ(t)为B-样条函数；n为B-样条阶次，m为节点，t为时间，*代表卷积运算。

步骤9：计算截止频率θ(t)，给定带宽阈值ξ，如果θ(t)≤ξ，则信号可认为是窄带信号，即判定s(t)为一个IMF，否则，令s₁(t)＝s(t)-m(t)，重复执行步骤2至步骤9。加权平均瞬时频率

与B_L瞬时带宽B_L的计算如式(13)～(15)；

经过上述步骤分解后，原信号可表达如式(16):

s(t)＝∑s_i(t) (16)

式中，s_i(t)为分解后的子序列。

假设一维多径信号s(t)的采样点数为n，信号经TVF-EMD分解后得到m个IMF分量，各个IMF分量的数据点为a_ik(i＝1...m；k＝1...n)，则时频矩阵P可表示为：

步骤10：对时频矩阵P进行奇异值分解，得到下式(18)所示的分解形式

P＝UQV^T (18)

在上式中，σ_i(i＝1,2,3,...,q)为时频矩阵P的奇异值，并且σ₁＞σ₂＞σ₃＞...＞σ_q＞0，q是通过分解获得的奇异值的总数。

由于矩阵P包含了信号的时频特征信息，那么矩阵P经奇异值分解后得

到的一系列子矩阵P_i也包含了时频特征信息，

步骤11：计算时频子空间的各个自相关函数R_s1,R_s2,...,R_sj和原信号的自相关函数R_x，自相关函数如式(20)所示：

式中，s(i)为信号某一时刻的状态，N为信号的采样点数。对自相关函数做归一化处理，求R_s1,R_s2,...,R_sj的相关系数ρ(j)，计算公式如式(21)所示：

式中，j为第j个分量。

相关系数的绝对值越大，相关性越强；相关系数越接近于1或-1，相关度越强，相关系数越接近于0，相关度越弱。本方法中相关系数的阈值取值为0.5。

步骤12：根据时频子空间计算得到的相关系数，判断出有效奇异值的个数为k。

步骤13：保留对角矩阵Ω的前面k个有效奇异值，然后将剩下的其他较小的奇异值做置0处理，以获得新的奇异矩阵Q′，将新的奇异矩阵Q′，左正交矩阵U和右正交矩阵V代入上式(19)，运算得到新的矩阵P′。

步骤14：将矩阵A′转换成一维信号s′(t)，s′(t)即为降噪后的信号。

实施例1

本发明提供了一种基于TVF-EMD-SVD的GNSS多径信号抑制方法，采用仿真的GNSS多径信号进行抑制，相关参数设置及步骤如下：

多径信号的采样率为1000hz，采样时间为2s，本实施例中，直达信号的幅值为1，频率为50hz，没有相位延时；第一条多径信号的幅值为0.8，频率为3.5hz，相位延时为

第二条多径信号的幅值为0.62，频率为15hz，相位延时为

同时还加入均值为0，方差σ²＝1的高斯白噪声，多径信号的波形图如图2所示。

经步骤2-9，上述信号分解出9个IMF分量，其中带宽阈值ξ＝0.08，B样条阶数为20，如图3所示。

经步骤10计算自相关函数图如图4所示，自相关系数如表1所示，其中自相关系数大于0.5的分量为S7和S9。

表1 S1-S9与原信号之间的相关系数

经步骤11-14后，抑制前后的信号如图5所示。通过计算抑制前后的信噪比和均方根误差发现，抑制前的信噪比为-0.4939dB，均方根误差为0.7466；抑制后的信噪比为1.7610，均方根误差为0.5775。经过对比，发现此方法可以有效地抑制多径信号，提高信号的信噪比。

本技术领域技术人员可以理解的是，除非另外定义，这里使用的所有术语(包括技术术语和科学术语)具有与本发明所属领域中的普通技术人员的一般理解相同的意义。还应该理解的是，诸如通用字典中定义的那些术语应该被理解为具有与现有技术的上下文中的意义一致的意义，并且除非像这里一样定义，不会用理想化或过于正式的含义来解释。

以上所述，仅是本发明的较佳实施例而已，并非是对本发明作任何其他形式的限制，而依据本发明的技术实质所作的任何修改或等同变化，仍属于本发明所要求保护的范围。