基于分数阶伪Wigner-Ville分布的卫星导航接收机干扰检测方法

摘要

本发明公开了一种基于分数阶伪Wigner‑Ville分布的卫星导航接收机干扰检测方法，用于全球导航卫星系统接收机的干扰检测。通过实验对发明的干扰检测方法进行了测试，测试的信号是受干扰的北斗卫星导航系统BDS‑B1I信号，实验结果表明本发明提出的全球导航卫星系统干扰检测方法有效地消除了双线性时频分布中存在的交叉项，同时显著改善了时频能量分布聚集特性，极大提升了全球导航卫星系统接收机干扰检测性能。

说明书

技术领域

本发明涉及全球导航卫星系统领域，具体来说是一种基于分数阶伪Wigner-Ville分布(Fractional Pseudo Wigner-Ville Distribution,FrPWVD)的卫星导航接收机干扰检测方法。

背景技术

全球导航卫星系统接收机的安全问题(例如民航和飞机着陆)引起了人们广泛的关注。全球导航卫星系统(Global Navigation Satellite System，GNSS)接收机接收信号中电磁干扰将导致卫星导航接收机导航与定位性能下降，甚至导致无法正常工作。大多数卫星导航系统中使用扩频(Direct Sequence Spread Spectrum，DSSS)技术，其原理是将接收到的GNSS信号功率扩展到更宽的带宽上，确保GNSS接收机中解扩增益，从而减少不良干扰信号造成的损害。DSSS技术使卫星导航系统具有一定的抗干扰能力，但是由于GNSS接收机接收的导航信号功率非常低，即使十分微弱的电磁干扰信号也将导致GNSS接收机的导航与定位性能严重下降。

当前，GNSS干扰检测和缓解技术(Interference Detection and Mitigation，ID&M)成为GNSS应用中非常重要的组成部分。GNSS干扰检测中现已经采用不同的时间-频率表示形式(Time Frequency Reference，TFR)的方法，例如频谱图和韦格纳分布(Wigner-Ville Distribution，WVD)。频谱图方法根据海森堡测不准原理存在时频(TimeFrequency，TF)分析分辨率的权衡问题，呈现较差的TF定位特性，无法用于干扰信号的瞬时频率估计。为了解决时频分辨率的折衷问题，在GNSS接收机的干扰检测中可采用WVD方法。虽然WVD具有许多良好的性能，并在所有时频分布中能够提供几乎最佳的分辨率，但由于不同频率分量的相互作用，它在时频平面中会存在严重的交叉干扰项，这会导致GNSS接收机干扰信号瞬时频率估计产生严重的误差。

为了减少WVD中存在的交叉干扰项，一种合理的解决方法是在时域中引入窗口函数，因此，提出了伪Wigner-Ville分布(Pseudo Wigner-Ville Distribution，PWVD)的概念。PWVD中窗口函数的引入可以一定程度抑制部分交叉干扰项；这种滤波窗口的缺点是使时频分析能量分布聚集特性衰减，因此时频定位精度降低，还可以观察到PWVD中交叉干扰项平行于扫频方向振荡。

在信号处理领域中，傅里叶变换作为一种成熟的数学工具被广泛使用，傅里叶变换是一种将时间信号从时间轴旋转π/2转移到频域轴的线性算子，而分数阶傅里叶变换(Fractional Fourier Transform,FrFT)是可以旋转任意角度的算子，其保留了傅里叶变换的原有性质又具备新的技术优势，被认为是一种广义的傅里叶变换。FrFT通过从0到1分数阶数的连续变换，可充分体现信号从时域到频域的变换特性。FrFT可作为一种有效的时频分析工具，可用于信号估计瞬时频率估计和相位信息恢复。

随着全球卫星导航系统的现代化进程和快速发展，GNSS干扰检测与缓解技术的要求越来越高，为此需构建一种更高效的基于分数阶伪Wigner-Ville分布的卫星导航接收机干扰检测方法。

发明内容

为了克服现有技术的不足，本发明提供一种基于分数阶伪Wigner-Ville分布的卫星导航接收机干扰检测方法，采用FrFT和PWVD相结合对包含扫频干扰的BDS-B1I信号进行分析，提取GNSS接收信号中扫频干扰的特征；FrPWVD方法可有效消除双线性时频WVD中存在的交叉项对GNSS干扰检测性能的影响，并可较好保留GNSS干扰信号的自项成分，显著改善GNSS干扰信号的时频聚集特性，从而有效提升GNSS接收机干扰检测性能。

为了便于描述本发明的内容，首先对GNSS接收信号进行说明，每次处理的卫星导航接收信号采样频率为100MHz，数字中频为14.58MHz。GNSS有效信号中加入的干扰信号为线性调频干扰。本发明采用的技术方案是：

基于分数阶伪Wigner-Ville分布的卫星导航接收机干扰检测方法，其特征在于，包括以下步骤：

步骤1：射频信号接收：通过接收机天线接收BDS-B1I射频信号r

步骤2：采用干扰机生成干扰信号η

y

步骤3：使用射频电缆将射频信号传输至信号采集器进行处理，将处理后的信号用USB电缆传输到GNSS接收机作进一步处理；

步骤4：GNSS接收机读取长度为n的数据流，可用下式表示：

y＝[y

步骤5：对GNSS接收机读取的数据流y进行希尔伯特变换获得

公式(3)中，j为虚根单位。

步骤6：在GNSS接收机中取分数阶数p＝0计算解析信号y

公式(4)中，F

公式(5)中，

公式(5)中，α＝pπ/2，代表时频面的旋转角度；

公式(5)中，当α＝2nπ或α＝(2n±1)π时，K

通过FrFT，得到1×n矩阵型数据F

步骤7：在PWVD中采用高斯窗函数，用PWVD来处理分数阶傅立叶变换后的信号F

公式(6)中，

公式(6)中，h(t)为时域窗函数，满足广义时间带宽积准则，其表达式如下：

公式(7)中，

公式(7)中，最优窗函数可以描述为窗函数

公式(7)中，

最终采用FrPWVD方法得到相应的n×n矩阵型变换域数据FrPWVD(u，v)

步骤8：运用最小二乘法方法检测出矩阵型变换域(u，v)数据FrPWVD(u，v)

步骤9：计算峰值点连线l与分数阶变换域u坐标轴的夹角α，若夹角α与90°角的偏差E大于ε，则回到步骤6将FrFT的分数阶数p增加Δp，重复步骤6至步骤9；偏差E的计算公式为：

E＝|α-90|/90 (8)

步骤10：若夹角α与90°角的偏差小于ε，则可认为FrPWVD(u，v)

步骤11：通过最优分数阶数p和交点坐标m即可计算扫频干扰的调频率和初始频率。

GNSS接收机信号是由接收机单天线采集。

干扰机生成的GNSS干扰信号η

η

公式(9)中，包括两个部分，分别是非平稳射频干扰项j

进一步地，步骤2的中GNSS干扰机生成非稳定线性调频干扰，其具体过程如下：

GNSS干扰机生成线性调频干扰，表达式如下：

公式(10)中，A为线性调频干扰的幅度，f

步骤3中计算GNSS接收机处理的数据流长度，所运用的公式如下：

n＝f

公式(11)中，f

步骤6和步骤7中解析信号的FrFT结果为一维矩阵数据，FrPWVD结果FrPWVD(u，v)为对称阵。

本发明的有益效果是：

本发明的基于分数阶伪Wigner-Ville分布的卫星导航接收机干扰检测方法，有效地消除了传统的双线性时频分布WVD中存在的交叉项的影响，保留了GNSS干扰信号的自项能量，明显改善了GNSS干扰信号的时频能量聚集特性；在复杂电磁干扰环境下，GNSS接收机可以准确识别GNSS干扰类型，并可精确估算其相应的特征参数。

附图说明

图1是本发明实施例提供的基于分数阶伪Wigner-Ville分布的卫星导航干扰检测方法的原理框图。

图2是在载噪比为46dB-Hz，干噪比为-2dB条件下本发明实例对含有扫频干扰的BDS-B1I信号进行检测的结果图。

图3是在载噪比为46dB-Hz，干噪比为-8dB条件下本发明实例对含有扫频干扰的BDS-B1I信号进行检测的结果图。

图4是在不同干噪比下，对BDS-B1I信号中存在的扫频干扰进行初始频率估计获得的均方根误差分析结果。

图5是在不同干噪比下，对BDS-B1I信号中存在的扫频干扰进行调频斜率估计获得的均方根误差分析结果。

图6是在不同干噪比下，对BDS-B1I信号中存在的扫频干扰进行干扰检测获得的峰均比分析结果。

具体实施方式

下面结合实施例和附图对本发明做进一步详细描述。

实施例1：

本发明实施例所述的基于分数阶伪Wigner-Ville分布的卫星导航接收机干扰检测方法的原理框图如图1所示。

步骤1：运用单天线采集BDS-B1I信号r

步骤2：采用干扰机生成干扰信号η

y

软件干扰器生成的干扰η

η

公式(2)中，包括两个部分，分别是非平稳的干扰信号j

步骤2.1：GNSS干扰机生成线性调频干扰信号，表达式如下：

公式(3)中，A为线性调频干扰的幅度，f

步骤3：使用射频电缆将GNSS信号传输至信号采集器，再AD转换之后的信号用USB电缆传输至GNSS接收机。

步骤3.1：计算GNSS接收机处理的数据流长度，运用公式如下：

n＝f

公式(4)中，f

步骤4：GNSS接收机读取长度为n＝2000的数据流，读取的数据表示公式为：

y(t)＝[y

步骤5：对GNSS接收机读取的数据流y进行希尔伯特变换获得

公式(6)中，y

步骤6：以分数阶数p＝0计算y

公式(7)中，F

公式(8)中，

公式(8)中，α＝pπ/2，代表时频面的旋转角度。

公式(8)中，当α＝2nπ或α＝(2n±1)π时，K

通过FrFT变换得到1×n矩阵型数据F

步骤7：在PWVD中采用高斯窗函数，对FrFT的信号F

公式(9)中，

公式(9)中h(t)是一个时域上的窗函数，其满足广义时间带宽积准则，其表达式如下：

公式(10)中，

公式(10)中，最优窗函数h(t)可描述为

公式(10)中，

采用FrPWVD得到对应的n×n矩阵型变换域数据FrPWVD(u，v)

步骤8：运用最小二乘法方法对矩阵型变换域(u，v)数据FrPWVD(u，v)

步骤9：计算峰值点拟合直线l与变换域u坐标轴的夹角α，若夹角α与90°角的偏差E大于ε，本实验取ε＝0.01，则回到步骤6将FrFT阶数p增加Δp，本实验取Δp＝0.01，重复步骤6至步骤9。偏差E的计算公式为：

E＝|α-90|/90 (11)

步骤10：若夹角α与90°角的偏差小于ε，则可认为FrPWVD(u，v)

步骤11：通过最优分数阶数p和交点坐标m即可计算扫频干扰的调频斜率和初始频率。

图2(a)为干噪比为-2dB下本实施对含有线性调频干扰的BDS-B1I信号的检测结果。由于BDS-B1I有用信号的频谱几乎是平坦的，因此FrPWVD分析使得干扰信号能量分布的差异性易于在TF平面上显现并可作出区分。这可在图2(b)中进一步验证，显示了本实施方式的TF分布的轮廓。针对线性调频干扰，在取最优分数阶数条件下下，FrPWVD在变换域平面上出现一条垂直于u坐标轴的切面，而在其他区域，其值几乎保持为零。本发明实例的FrPWVD方法有效抑制了双线性TF分布中存在的交叉项干扰，并且使干扰信号中TF能量聚集特性得到了改善。因此，在干扰环境下，该方法可有效应用于GNSS接收机干扰检测，能够精确识别电磁干扰类型，并可精确估计干扰的特征参数。

图3(a)为干噪比为-8dB条件下本实施方式对含有线性调频干扰的BDS-B1I信号进行检测的结果。与图2线性调频信号干扰的影响类似，在取最优分数阶数条件下，扫频干扰信号的FrPWVD在变换域平面上出现一条垂直于u坐标轴的切面，相较于图2，在排除扫频干扰分量以外的其它区域，平行于干扰分量的能量更明显，对应于导航信号本身和高斯白噪声信号，由于这些分量远小于干扰信号的贡献，对GNSS干扰检测的影响很小，因此，在存在干扰情况下，FrPWVD方法能够有效检测、识别GNSS干扰信号，并可精确估算其特征参数。

为了全面评估所提出的基于分数阶伪Wigner-Ville分布的卫星导航接收机干扰检测方法的性能，将从统计学角度使用均方根误差(Root Mean Square Error，RMSE)来进行验证分析。当BDS-B1I接收信号出现扫频干扰时，不同干噪比(Jamming Noise Ratio，JNR)条件下扫频干扰初始频率和调频斜率估计的均方根误差不同。

图4为本实施方式的不同时频分析方法针对扫频干扰初始频率估计的定量对比结果，在不同JNR条件下，计算BDS-B1I接收信号中线性调频干扰的初始频率的RMSE结果。其中，采用TF分布方法计算出的RMSE值为JNR的函数，干扰初始频率的RMSE随JNR值的增加而减小。在图4中，WVD和PWVD方法显示出非常差的频率估计精度。其中，WVD对干扰初始频率估计的RMSE随着JNR值的增加几乎不变，维持在1×10

图5为针对干扰调频斜率估计，本实施方式的FrPWVD和传统时频分析方法进行对比分析的结果。在不同JNR条件下，采用不同的方法对BDS-B1I接收信号中线性调频干扰的调频斜率进行估计，并对相应的RMSE进行对比分析。其中，不同TF分布计算出的干扰调频斜率估计的RMSE值是以JNR的函数形式出现，从整体上看出扫频干扰调频斜率估计的RMSE随JNR值的增加而减小。其中，WVD对干扰调频斜率估计的RMSE随着JNR值的增加几乎不变，维持在2.1×10

为了更全面地比较FrFT、FrWVD和本发明提出的FrPWVD方法的性能，将使用峰均比(Peak-to-Average Ratio，PAR)来表征GNSS干扰检测灵敏度。PAR是输出信号的峰值与平均值之比，它影响功率放大器的动态范围。PAR的计算公式如下：

公式(12)中，max{|s(t)|

图6显示了在不同JNR条件下FrFT和FrWVD以及本发明所提FrPWVD方法的峰均比对比结果。可以看出，三种方法的PAR均随JNR增加而提高。分析结果表明，随着JNR值的增加，三种方法的检测能力均获得提高。FrFT的PAR结果变化范围大约在1.5×10