CN109085614B

CN109085614B - 基于时-频谱值平滑与分割的gnss干扰特征辨识方法及系统

Info

Publication number: CN109085614B
Application number: CN201810731054.3A
Authority: CN
Inventors: 丁梦羽; 许睿; 刘建业; 曾庆化; 孟骞; 熊智
Original assignee: Nanjing University of Aeronautics and Astronautics
Current assignee: Nanjing University of Aeronautics and Astronautics
Priority date: 2018-07-05
Filing date: 2018-07-05
Publication date: 2020-10-20
Anticipated expiration: 2038-07-05
Also published as: CN109085614A

Abstract

本发明结合GNSS信号特点对经典时频方法进行优化，对时频谱值进行平滑，并根据频段能量分布情况对谱值分割处理，从而提出了一种时‑频谱值平滑与分割的GNSS干扰特征辨识方法及系统。该方法可改善对宽带小功率干扰信号特征的辨识能力，提高干扰信号特征辨识效果，为后续干扰的抑制提供必要的条件。

Description

基于时-频谱值平滑与分割的GNSS干扰特征辨识方法及系统

技术领域

本发明属于GNSS干扰信号识别领域，具体涉及一种基于时-频谱值平滑与分割的GNSS干扰特征辨识方法。

背景技术

随着经济科技的飞速发展，全球卫星导航系统(GNSS)逐渐成为当今社会社提供位置服务的重要基础设施。但其可靠性却受到诸多因素的威胁，卫星信号经过长距离的传播，到达地面时的信号功率很低，很容易受到各种无意或有意的无线电干扰信号影响。压制式干扰信号是最恶劣的GNSS干扰类型之一，干扰信号影响卫星信号的跟踪性能，导致伪距误差、定位误差增大，严重降低GNSS系统的可靠性和准确性，而更大功率的干扰会直接淹没可用卫星信号，导致GNSS系统无法提供导航定位服务。因此，提高GPS系统抗干扰能力显得尤为重要。

干扰检测是进行干扰抑制的必要步骤，也是提高GPS系统鲁棒性的重要前提。以跟踪环路相关器为分结点，常见的干扰检测技术可分为后相关检测和前相关检测两大类。后相关干扰检测技术主要使用有效载噪比，载波相位、相关器输出、卫星轨道信息和定位误差等作为干扰检测的主要参数。而前相关干扰检测是检测RF前端或基带模块处的干扰，可以使用AGC模块、A/D转换模块、中频信号、频域信号、时频分析等进行干扰检测。这两种方案相比，后相关干扰检测较慢，干扰信号的灵敏度较低，而前相关干扰检测速度更快，对干扰的灵敏度也更高。

此外，上述大部分干扰检测方法只能检测干扰是否存在，无法对干扰信号的周期特征、频率变化特征等关键参数进行估计，因此无法识别干扰特性。考虑时频分析方法，在检测干扰信号存在的同时，能对干扰信号的参数进行估计，然而，传统时频分析只对大功率干扰信号的参数具有较好的估计效果，而对于较小功率干扰，特别是宽带小功率干扰信号的特征参数估计效果不佳。

发明内容

针对上述问题与难点，本发明结合GNSS信号特点对经典时频方法进行优化，对时频谱值进行平滑，并根据频段能量分布情况对谱值分割处理，从而提出了一种时-频谱值平滑与分割的GNSS干扰特征辨识方法。该方法可改善对宽带小功率干扰信号特征的辨识能力，提高干扰信号特征辨识效果，为后续干扰的抑制提供必要的条件。具体技术方案如下：

本发明公开一种基于时-频谱值平滑与分割的GNSS干扰特征辨识方法，包括如下步骤：

(1)对输入的中频信号进行时频分析，得到信号的二维时-频谱结果；

(2)对时-频谱结果进行数值处理，包括：将二维时-频谱值进行归一化处理，减小数据处理计算量；对归一化数值进行滤波处理，消除GNSS信号和噪声造成的随机误差；对滤波后的数据进行二值化处理，实现数据的分割分类；

(3)根据二值化分割的数值结果，提取数据出现梯度变化时频点，获取干扰信号时频特征，构造干扰时-频曲线；通过对干扰时-频曲线进行FFT分析，提取干扰信号的周期；

(4)结合周期结果和干扰时-频曲线的连续性判断辨识干扰曲线是否完整，若完整，直接输出干扰信号的特征参数；若不完整，则转到步骤(5)；

(5)增长时频分析的中频信号数据长度，并对中频信号时-频二维谱值进行按照步骤(3)提取的干扰信号周期分段和叠加，平滑谱图的随机波动，增强干扰信号的特征；

(6)将平滑后的时-频谱值按频段分割成多段，并对各段时-频谱值采用如步骤(2)所述的方法进行数值处理；

(7)将各段数值处理结果重新组合，然后对重组后的时频二值化结果采用如步骤(3)～(4)所述的方法进行再处理，直到辨识出干扰信号的特征参数。

本发明还公开一种基于时-频谱值平滑与分割的GNSS干扰特征辨识系统，包括处理器和存储器，存储器内存储有GNSS干扰特征辨识程序，该程序配置有多个程序模块，包括联合时频域分析模块、综合数值处理模块、多路分块数值处理模块、周期提取模块、边缘完整性判决模块、数据平滑模块、谱值频段分割模块和数据重组模块，该程序在被处理器运行时，执行以下步骤：

(1)联合时频域分析模块对输入的中频信号进行时频分析，得到信号的二维时-频谱结果；

(2)综合数值处理模块对时-频谱结果进行数值处理，包括：将二维时-频谱值进行归一化处理，减小数据处理计算量；对归一化数值进行滤波处理，消除GNSS信号和噪声造成的随机误差；对滤波后的数据进行二值化处理，实现数据的分割分类；

(3)周期提取模块根据二值化分割的数值结果提取数据出现梯度变化时频点，获取干扰信号时频特征，构造干扰时-频曲线，然后通过对干扰时-频曲线进行FFT分析，提取干扰信号的周期；

(4)边缘完整性判决模块结合周期结果和干扰时-频曲线的连续性判断辨识干扰曲线是否完整，若完整，直接输出干扰信号的特征参数，若不完整，则转到步骤(5)；

(5)增长时频分析的中频信号数据长度，数据平滑模块对对中频信号时-频二维谱值进行按照步骤(3)提取的干扰信号周期分段和叠加，平滑谱图的随机波动，增强干扰信号的特征；

(6)谱值频段分割模块将平滑后的时-频谱值按频段分割成多段，并对各段时-频谱值采用如步骤(2)所述的方法进行数值处理；

(7)数据重组模块将各段数值处理结果重新组合，然后对重组后的时频二值化结果采用如步骤(3)～(4)所述的方法进行再处理，直到辨识出干扰信号的特征参数。

有益效果：

(1).采用联合时频域分析的方法分析GNSS信号，能够同时得到信号的频域特征、时域特征和时频域特征，提取的干扰信号特征更全面详细。

(2).采用归一化和二值分割的处理方案对时频谱图进行平滑和处理，比传统的谱图等高线方法提取的干扰信号更准确，误差更小。

(3).根据干扰信号的周期对时-频谱值进行优化，减小谱值的波动情况，强化干扰信号的特征，提高干扰信号提取的灵敏度。

(4).通过二值化分割的数据边缘提取的曲线完整性决定是否需要增长时频分析的数据长度，灵活调整干扰提取所需时间，效率更高。

(5).接收机前端滤波器会衰减干扰信号带外频段的能量，造成干扰信号能量分布不均匀，出现提取的干扰信号(特别是宽带小功率干扰)时-频曲线不完整情况，而本发明提出的频段分割谱图的方案优化图像处理结果，能够有效避免这一问题，提高了宽带小功率干扰信号特征的辨识能力。

(6).本发明基于软件算法完成，具有很大的灵活性和低成本的优势。

(7).该发明技术通过中频信号分析，不受卫星信号结构的影响，适用于GPS、北斗等各种GNSS系统的干扰信号提取，通用性强，具有重要的工程意义。

附图说明

图1基于时频分析和图像处理的GNSS干扰特征辨识方法示意图。

具体实施方式

实施例中公开一种GNSS干扰特征辨识系统，其主要包括：输入的GNSS中频信号，联合时频域分析模块、综合数值处理模块、多路分块数值处理模块、周期提取模块、边缘完整性判决模块、数据平滑模块、谱值频段分割模块和数据重组模块。其中，综合数值处理模块和多路分块数值处理模块内部结构一致，均由归一化模块、滤波器和二值化数据分割模块构成；周期提取模块包括梯度提取模块和周期估计器；谱值频段分割模块具有多个输出端，数组重组模块具有多个输入端。

GNSS中频信号作为输入信号输入到联合时频域分析模块的输入端，联合时频域分析模块的第一输出端连接数值平滑模块的第一输入端，联合时频域分析模块的第二输出端连接综合数值处理模块的输入端，综合数值模块的输出端连接周期提取模块的第一输入端，周期提取模块的输出端连接边缘完整性判决模块的输入端，边缘完整性判决的第一输出端输出干扰信号特征辨识结果，边缘完整性判决的第二输出端连接数值平滑模块的第二输入端，数值平滑模块的输出端连接谱值频段分割模块的输入端，谱图分割模块的多个输出端对应连接各分块数值处理模块的输入端，各分块数值处理模块的输出端分别连接数值重组的多个输入端，数值重组的输出端连接周期提取模块的第二输入端。

综合、分块数值处理模块内部各模块连接结构：联合时频分析模块的第一输出端和谱值频段分割模块的各输出端连接归一化模块的输入端，归一化模块的输出端连接滤波模块的输入端，滤波模块的输出端连接二值化数据分割模块的输入端，二值化数据分割模块的输出端作为综合、分块数值处理模块的输出；

周期提取模块的内部各模块连接结构：综合数值处理模块的输出端和数值重组模块的输出端连接梯度提取模块的输入端，梯度提取模块的第一输出端直接作为周期提取模块的第一输出结果，梯度提取模块的第二输出端连接周期估计器的输出端，周期估计器的输出端直接作为周期提取模块的第二输出结果。

结合图1所示，基于上述系统，本发明还公开一种GNSS干扰特征辨识方法，可对GNSS中频采样器采集的GNSS中频信号进行处理，具体包括以下步骤：

S1、GNSS中频信号输入到联合时频域分析模块进行时频分析(可采用短时傅里叶、WT小波变换、WVD维格纳变换等方法)，得到信号的二维时-频谱结果。

此处，以采用短时傅里叶时频分析方法为例，对输入信号x(m)进行如下：

其中，m为信号x序列，n为短时傅里叶时间序列值，ω为短时傅里叶频率值，x(m)为输入信号，w(n-m)为窗函数。

进一步，时-频能量密度谱p(n,ω)为：

其中，A为实值标量。

对于单边能量谱密度，计算公式如下：

其中，w(n)表示窗函数，F_s为采样频率，在0频率和奈奎斯特频率处，分子上的因子为1。

S2、对得到的二维时-频谱结果输入到综合数据处理模块进行数值处理，具体包括：将二维时-频谱值进行归一化处理，减小数据处理计算量；然后通过滤波器(可采用均值滤波、中值滤波、高斯滤波等)对归一化数值进行滤波处理，消除GNSS信号和噪声造成的随机误差；最后将滤波后的数据进行二值化处理，实现数据的分割分类。

具体处理方法如下：

将二维时频能量密度谱结果转换到0～1之间，归一化处理得到的数值I₁(n,ω)如下：

其中，min(·)和max(·)为取最小值和最大值，[·]为向下取整的符号；

此处采用均值滤波器对GNSS信号中的伪随机噪声和各种热噪声会造成的噪点进行滤除，选取由I₁(n,ω)其邻近若干数值组成模板S，由模板均值I₂(n,ω)替代原归一化数值：

其中，M表示模板S内的频谱值个数。

对归一化数值进行0-1分割，获得数值整体和局部特征的二值化结果I₃(n,ω)：

其中，η表示二值化分割的阈值。

通过近似一维Means方法确定适当的阈值η，具体包括：

1)随机方法生成初始化阈值η；

2)I₂(n,ω)根据阈值η分为对象数据G₁与背景数据G₂；

3)根据G₁和G₂的平均值m₁和m₂获得新阈值η^′＝(m₁+m₂)/2；

4)重复2)～3)步骤，直至计算出的新阈值等于上一次阈值为止。

S3、梯度提取模块根据二值化分割的数值结果，提取数据出现梯度变化的时频点，获取干扰信号时频特征，构造干扰时-频曲线；然后通过周期估计器对时-频曲线进行FFT分析，提取干扰信号的周期。

具体算法如下：

根据二值化数值的梯度变换构造干扰信号的时频关系ω＝f(n)，干扰信号本身能量较高，因此会在干扰信号处出现梯度变化，即提取二值化结果为1的频率和时序可以拟合出干扰曲线：

ω＝f(n),其中，I₃(n,ω)＝1 (7)

进一步，对干扰信号时频曲线进行周期估计，利用FFT变换互相关原理，周期对齐时，得到的傅里叶结果最大，因此，对信号进行傅里叶变换：

非零频率处FFT结果的最大值对应的频率倒数即为干扰信号周期T：

S4、边缘完整性判决模块结合周期结果和干扰时-频曲线的连续性判断辨识干扰曲线是否完整，若完整，直接输出干扰信号的特征参数，包括类型、频段、周期等；若不完整，转到步骤S5。

具体包括：通过图像边缘提取的干扰信号的时频关系ω＝f(n)进行连续性判断，若每个时间片n都对应有ω，且ω不是随机分布，则认为干扰信号提取完整，直接输出干扰信号时频关系及周期；否则认为干扰信号提取不完整，需将周期信号输入到数值平滑模块进一步处理。

S5、在原本数据长度的基础上，增长时频分析的中频信号数据长度，数据平滑模块对中频信号时-频二维谱值进行按照S3提取的干扰信号周期T分段和叠加，并平滑谱图的随机波动，以增强干扰信号的特征。

其中，数值平滑方法具体包括：

GNSS中频输入信号x(m)：

x(m)＝r(m)+J(m) (10)

其中，r(m)包含GPS可见卫星信号和噪声，J(m)为压制式干扰信号。

对输入信号x(m)进行短时傅里叶变换可得：

因此，时频能量谱密度由正常信号能量p_r、干扰信号能量p_J和交叉项Δ_r,J组成：

由于时频分析的结果本身存在波动，当干扰信号功率较小时，可能会导致干扰信号的能量谱密度淹没在随机波动中，出现无法正常提取干扰信号的情况。而根据已估计的干扰信号的周期T对时频谱值进行叠加平均，可减小谱图随机波动。例如，对数据长度为L的时频谱值按照周期(时序)T分成k段，即实现k＝L/T次叠加，每次能量密度结果均采用不同数据段估计，近似为相互独立的，因此平均时-频能量密度p′(n,ω)的方差var[p′]会减小，进一步，干扰信号的特征得到增强：

其中，i表示第i段时频谱值。

干扰信号作用频段的时-频能量密度谱为p(n,ω_J)＝p_J+p_s+Δ_r,J+ε，ε为计算和噪声引起的随机波动，其标准差为σ(n)。

当p(n,ω_J)满足如下条件时，能够实现干扰信号谱值的正确二值化分割，即正确辨识干扰信号：

Min(p(n,ω_J))>μ(n,ω)+3σ(n) (15)

其中，μ(n,ω)＝p_s(n,ω)为信号在(n,ω),ω≠ω_J处能量密度均值，根据正态分布，误差最大值为3σ，因此式(15)可进一步写为：

k次叠加平均后，能量密度波动的方差变为原来的1/k，因此：

可见，时频谱图的随机波动被大幅平滑，对干扰频点的能量密度与交叉项之和的要求减小了

最终适应小功率的干扰信号的特征参数辨识。

S6、谱值频段分割模块对平滑后的时-频二维谱值根据不同频段的能量分布情况将谱值分割成多块，每块时-频谱值分别输入到分块数值处理模块，重复步骤S2所述的方法。

谱图频段分割的方法具体如下：

GNSS前端滤波器会对GNSS信号波段外的其他射频信号进行滤波，因此，经过滤波后的GNSS中频信号的能量主要集中在滤波器带宽内，而带宽外的信号能量水平则很低。

前端滤波器对带内干扰的能量没有损耗，但对于宽带干扰而言，带宽外的干扰信号能量被降低到很低的水平，而带宽内的干扰信号被保留，能量水平很高。在GNSS受宽带小功率干扰信号的影响时，带外干扰能量水平低于带内GNSS信号能量水平，对GNSS中频信号谱值归一化、滤波和二值化数据分割处理只能辨识出干扰信号带内部分的特征。

根据奈奎斯特采样定理，信号频域变化的最大值为采样频率F_s的一半，即F_s/2。在极小一段频段范围，信号能量谱密度近似保持不变，将时-频谱值按频段分割成N段，N的大小受短时傅里叶频率分辨率的限制，短时傅里叶的频率分辨间隔为：

其中，NFFT为傅里叶变换点数，N为(NFFT/2+1)的约数。

取N＝NFFT/2+1，对时-频谱图按频段分割成：

S7、各分块数值处理模块的输出结果输入到数值重组模块，将多个频段的二值化数据结果重新组合，然后重复步骤S3～S4所述的方法，最终实现干扰信号特征的辨识。

其中，数组重组方式是分将分割后的各部分数据处理结果组合，形成完整的时频二值化结果：

由此可见，本发明可改善对宽带小功率干扰信号特征的辨识能力，提高干扰信号特征辨识效果，为后续干扰的抑制提供必要的条件，具有重要工程意义。

尽管以上结合附图对本发明的实施方案进行了描述，但本发明并不局限于上述的具体实施方案和应用领域，上述的具体实施方案仅仅是示意性的、指导性的，而不是限制性的。本领域的普通技术人员在本说明书的启示下，在不脱离本发明权利要求所保护的范围的情况下，还可以做出很多种的形式，这些均属于本发明保护之列。

Claims

1.一种基于时-频谱值平滑与分割的GNSS干扰特征辨识方法，其特征在于，包括以下步骤：

2.如权利要求1所述的GNSS干扰特征辨识方法，其特征在于，步骤(1)中，时频分析方法包括但不限于短时傅里叶分析方法、WT小波变换分析方法、WVD维格纳变换分析方法。

3.如权利要求1所述的GNSS干扰特征辨识方法，其特征在于，步骤(2)中，所述滤波方法包括但不限于均值滤波、中值滤波、高斯滤波。

4.如权利要求1所述的GNSS干扰特征辨识方法，其特征在于，步骤(2)中，归一化处理具体包括：

将二维时频能量密度谱结果转换到0～1之间，归一化处理得到的数值I₁(n，ω)如下：

其中，n为短时傅里叶时间序列值，ω为短时傅里叶频率值，p(n，ω)为时-频能量密度谱，min(·)和max(·)为取最小值和最大值，[·]为向下取整的符号；

对GNSS信号中的伪随机噪声和各种热噪声会造成的噪点进行滤除，选取由I₁(n，ω)其邻近若干数值组成模板S，由模板均值I₂(n，ω)替代原归一化数值：

其中，M表示模板S内的频谱值个数；

对归一化数值进行0-1分割，获得数值整体和局部特征的二值化结果I₃(n，ω)：

其中，η表示二值化分割的阈值。

5.如权利要求4所述的GNSS干扰特征辨识方法，其特征在于，通过近似一维Means方法确定适当的阈值η，具体包括：

1)随机方法生成初始化阈值η；

2)I₂(n，ω)根据阈值η分为对象数据G₁与背景数据G₂；

3)根据G₁和G₂的平均值m₁和m₂获得新阈值η′＝(m₁+m₂)/2；

4)重复步骤2)～3)，直至计算出的新阈值等于上一次阈值为止。

6.如权利要求5所述的GNSS干扰特征辨识方法，其特征在于，步骤(3)的具体方法如下：

根据二值化数值的梯度变换构造干扰信号的时频关系ω＝f(n)，干扰信号处会出现梯度变化，提取二值化结果为1的频率和时序，并拟合出干扰曲线：

ω＝f(n)，其中，I₃(n，ω)＝1 (7)

对干扰信号时频曲线进行周期估计，利用FFT变换互相关原理，周期对齐时得到的傅里叶结果最大，因此，对信号进行傅里叶变换如下：

7.如权利要求6所述的GNSS干扰特征辨识方法，其特征在于，步骤(4)的具体方法如下：

通过图像边缘提取的干扰信号的时频关系ω＝f(n)进行连续性判断，若每个时间片n都对应有ω，且ω不是随机分布，则认为干扰信号提取完整，直接输出干扰信号时频关系及周期；否则认为干扰信号提取不完整，需要转入步骤(5)进一步处理。

8.如权利要求7所述的GNSS干扰特征辨识方法，其特征在于，步骤(5)中，数值平滑的方法具体包括：

根据提取的干扰信号的周期T，对时频谱值进行k次叠加平均，其中，k＝L/T，L表示数据总长度；

干扰信号作用频段的能量密度谱为p(n，ω_J)＝p_J+p_s+Δ_r，J+ε，其中，p_r表示正常信号能量、p_J表示干扰信号能量、Δ_r，J为交叉项，ε为计算和噪声引起的随机波动，其标准差为σ(n)；

当p(n，ω_J)满足如下条件时，能够实现干扰信号谱值的正确二值化分割，即正确辨识干扰信号：

Min(_p(n，ω_J))＞μ(n，ω)+3σ(n) (15)

其中，μ(n，ω)＝p_s(n，ω)为信号在(n，ω)，ω≠ω_J处能量密度均值，根据正态分布，误差最大值为3σ，式(15)进一步表示为：

k次叠加平均，能量密度波动的方差变为原来的1/k，因此：

时频谱图的随机波动被大幅平滑，对干扰频点的能量密度与交叉项之和的要求减小了

最终适应小功率的干扰信号的特征参数辨识。

9.如权利要求8所述的GNSS干扰特征辨识方法，其特征在于，步骤(6)中，谱图频段分割的方法如下：

根据奈奎斯特采样定理，信号频域变化的最大值为采样频率F_s的一半，即F_s/2，在极小一段频段范围，信号能量谱密度近似保持不变，将时-频谱值按频段分割成N段，N的大小受短时傅里叶频率分辨率的限制，短时傅里叶的频率分辨间隔为：

其中，NFFT为傅里叶变换点数，N为(NFFT/2+1)的约数；

取N＝NFFT/2+1，对时-频谱图按频段分割成：

步骤(7)中，数组重组方式如下：

将分割后的各部分数据处理结果组合，形成完整的时频二值化结果I₃′(n，ω)：

10.一种基于时-频谱值平滑与分割的GNSS干扰特征辨识系统，其特征在于，包括处理器和存储器，存储器内存储有GNSS干扰特征辨识程序，该程序配置有多个程序模块，包括联合时频域分析模块、综合数值处理模块、多路分块数值处理模块、周期提取模块、边缘完整性判决模块、数据平滑模块、谱值频段分割模块和数据重组模块，该程序在被处理器运行时，执行以下步骤：

(5)增长时频分析的中频信号数据长度，数据平滑模块对中频信号时-频二维谱值进行按照步骤(3)提取的干扰信号周期分段和叠加，平滑谱图的随机波动，增强干扰信号的特征；