CN107576943B - 基于瑞利熵的自适应时频同步压缩方法 - Google Patents

Abstract

一种基于瑞利熵的自适应时频同步压缩方法，该方法利用瑞利熵衡量信号局部能量的聚集性，估计出适合信号时频分析的最佳窗口参数序列，通过时频同步压缩提高聚集性，实现多分量信号的分离。具体步骤包括：1、采集信号，2、预处理，3、选取适合信号点序列的窗口参数集合，4、计算瑞利熵，5、获取局部时刻的最优窗函数，6、获取平滑的最佳窗口参数序列，7、自适应同步压缩短时傅立叶变换，8、获得各分量信号。本发明克服了传统分离方法适用信号单一、信号形式固定、抗噪声能力差的问题，有利于多分量非平稳信号的分离，更提高了多分量信号分离的鲁棒性。

Description

基于瑞利熵的自适应时频同步压缩方法

技术领域

本发明属于通信技术领域，更进一步涉及信号与信息处理技术领域中的一种基于瑞利熵的自适应时频同步压缩方法。本发明可以用于对雷达侦察接收机接收到的多分量雷达辐射源信号进行分离。

背景技术

近年来，随着雷达技术的发展，以及导航、通信、电子对抗等电子装备的广泛适用，战场电磁环境日益复杂。在雷达电子侦察中，接收机在同一时间内可能会截获到多个雷达辐射源信号，且信号脉冲密度大，波形多变，导致信号时频域上发生重叠，使得传统的时域或频域已无法对信号分量进行有效地分离。目前最常用的雷达辐射源信号分离方法有盲源分离法、经验模式分解法EMD以及时频分析法。盲源分离法主要是利用观测信号所携带的信息与源信号之间相互统计独立的假设下，实现对源信号的分离，无法应用于相关的辐射源信号。经典模式分解法EMD是一种局部化的，完全基于数据自适应信号分析方法，该方法总是将接收到的信号分解为有限多个固有模式函数。但它是一种基于经验的方法，没有严格的数学理论基础，且分解过程中存在模式混叠。时频分析法是将信号在时频域上展开，根据信号时频变换的特点反映信号的局部特性达到分离的目的，但该方法严重受所加窗的大小以及环境中噪声的影响。

北京理工大学在其申请的专利文献“一种同时到达两线性调频信号的分数阶信道化分离方法”(申请号201210189114，公开号CN102685049A)中公开了一种分数阶信道化分离方法。该方法利用分数阶傅里叶变换对非平稳信号能量聚焦特点，对信号进行分析，判断信号的不同组合情形，采用准聚焦分数阶信道化提取的分离处理方式实现分离的目的。虽然，该方法一定程度上减少了分离信号间的相互干扰且达到信号分离的效果，但是，仍然存在不足之处是，该方法只局限于两个线性调频信号，且信号形式固定，当接收的信号多于两个信号时，该方法的分离效果不明显。

上海大学在其申请的专利文献“基于时频分析的LFM信号分离及参数估计方法”(申请号201610094975，公开号CN105717490A)中公布了一种聚合分割法。该方法在小波变换的基础上，通过分段聚合的方式将信号分割，对聚合的分段通过阀值设置，依次进行估计，并利用去重加权的方式进行优化，直到分离出全部信号。虽然，该方法计算量低，一定程度上分离出了每个信号分量，但是，仍然存在不足之处是，该方法受噪声影响很大，且分离过程中阀值的选择没有严格的选取标准。

综上所述，对于雷达辐射源信号的分离，上述已有技术容易受噪声干扰，且适用的信号单一，从而使信号分离出来的分量与实际分量存在很大误差。

发明内容

本发明的目的在于克服上述已有技术的不足，在传统短时傅里叶变换的基础上，提出一种基于瑞利熵的自适应时频同步压缩方法，旨在低信噪比噪声、辐射源信号复杂的情况下实现多分量信号的准确分离。

实现本发明的具体思路是：首先确定高斯函数的窗口参数集合，用窗口参数集合对应的窗函数集合，对信号进行加窗傅里叶变换，得到短时傅里叶变换系数。利用瑞利熵公式求最小值，得到所有时刻下最优窗函数对应的窗口参数，组成最佳窗口参数序列。用最佳窗口参数序列对应的窗函数集合对信号进行加窗傅里叶变换，得到自适应短时傅里叶变换。对自适应短时傅里叶变换系数进行叠加，得到自适应时频同步压缩系数，利用维特比Viterbi从自适应时频同步压缩系数中提取信号分量的瞬时频率，用重构公式实现信号分量的分离。

本发明的具体步骤包括如下：

(1)采集信号：

信号采集系统通过雷达侦察接收机设备，实时采集目标空域中的任意一段信号；

(2)获取信号点序列：

根据所所采集的信号的长度，对所采集的信号在大于等于1的采样频率范围内，选取一个采样频率进行采样，得信号点序列N＝Tf_s，其中，N表示信号点序列的长度，T表示所采集的雷达侦察接收机信号的时间长度，f_s表示采集信号的采样频率；

(3)选取信号点序列的窗口参数集合：

(3a)选取高斯函数为窗函数；

(3b)在(0,1]区间，以1/N的步长，选取N个不同的值，组成窗口参数集合；

(4)计算瑞利熵：

(4a)在窗口参数集合中，任意选取一个窗口参数，在其对应的窗函数下，对信号点序列进行离散加窗傅里叶变换，得信号点序列的短时傅里叶变换系数；

(4b)利用瑞利熵公式，得到每个时刻的瑞利熵值；

(4c)遍历窗口参数集合中的所有窗口参数，得到每个时刻对应的瑞利熵值集合；

(5)获取最佳窗口参数序列：

(5a)从瑞利熵值集合中选出最小值，将最小值对应的窗函数作为每个时刻的最优窗函数；

(5b)将短时傅里叶变换所有时刻下，最优窗函数对应的最优窗口参数，组成预处理后的完整信号点序列的最佳窗口参数序列；

(6)获取平滑的最佳窗口参数序列：

用滑动平均法，对最佳窗口参数序列内的点进行累加平均，得到平滑的最佳窗口参数序列；

(7)自适应同步压缩短时傅里叶变换：

(7a)用平滑的最佳窗口参数序列对应的窗函数集合，对信号点序列进行离散加窗傅里叶变换，得到自适应短时傅里叶变换系数；

(7b)用维特比Viterbi法，对自适应短时傅里叶变换系数在时频域进行最小路径搜索，得到信号点序列的每个分量的瞬时频率脊线；

(7c)利用同步压缩公式，根据信号点序列的每个分量的瞬时频率脊线，将每个时刻下所有对应同一频率的自适应短时傅里叶变换系数进行叠加，得到自适应同步压缩系数；

(8)获得各分量信号：

利用重构公式，将每个时刻下对应同一个分量瞬时频率脊线的自适应同步压缩系数进行叠加，得到信号点序列中的每个分量信号。

与现有技术相比，本发明具有以下优点：

第一，本发明是用最佳窗口参数序列对应的窗函数集合，对信号点序列进行离散加窗傅里叶变换，克服了现有技术中只局限于两个线性调频信号，且信号形式固定的缺点，使得本发明对信号具有普遍适用性，且更有利于多分量信号的分离。

第二，本发明是计算瑞利熵值集合，将瑞利熵值集合中最小值对应的窗函数作为局部时刻下的最优窗函数，根据瑞利熵值越小，能量聚集性越高的特点，克服了现有技术受噪声影响大的缺点，使得本发明在短时傅里叶变换时频图中具有很高的能量聚集性，且抗噪声能力强，更准确的表现出信号的局部信息。

附图说明

图1为本发明的流程图；

图2为本发明获取的最佳窗口参数序列图；

图3为本发明的短时傅立叶变换时频图；

图4为本发明的各分量信号的误差仿真对比图。

具体实施方式

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

参照图1，本发明的具体实施步骤如下。

步骤1，采集信号。

信号采集系统通过雷达侦察接收机设备，实时采集目标空域中的任意一段信号x(t)，信号x(t)包含多个分量，具体模型可表示为：

其中，∑表示求和符号，x_k(t)表示信号x(t)的第k个信号分量，A_k(t)表示信号x(t)的第k个信号分量的幅度，表示信号x(t)的第k个信号分量的相位，K表示信号x(t)的信号分量个数，T表示所采集的目标空域中信号的时间长度。

步骤2，获得信号点序列。

根据所采集信号的长度，对所采集的信号在大于等于1的采样频率范围内，选取一个采样频率进行采样，得信号点序列N＝Tf_s，其中，N表示信号点序列的长度，T表示所采集的雷达侦察接收机信号的时间长度，f_s表示采集信号的采样频率。

步骤3，选取信号点序列的窗口参数。

选取高斯函数为窗函数。信号点序列的短时傅里叶变换时频分辨特性依赖于窗函数的选择，由不确定性准则，窗函数为高斯函数时，信号点序列的时频局部化能力更强，我们选择为高斯窗函数时域形式，其中，e表示指数符号，σ为控制着高斯窗形状变化的窗口参数。

在(0,1]区间，以1/N的步长，选取N个不同的值，作为窗函数的窗口参数，形成窗口参数集合{σ₁,σ₂,…,σ_N}，以及相对应的窗函数集合{g₁,g₂,…,g_N}。

步骤4，计算瑞利熵。

在窗口参数集合{σ₁,σ₂,…,σ_N}中，任意选取一个窗口参数σ_p∈{σ₁,σ₂,…,σ_N}，在其对应的窗函数g_p∈{g₁,g₂,…,g_N}下，对信号点序列进行离散加窗傅里叶变换，得到信号点序列的短时傅里叶变换系数。

按照下式瑞利熵公式，计算得到局部时刻下的瑞利熵值：

其中，R表示瑞利熵值，α表示瑞利熵系数，α≥2.5，log₂表示以2为底的对数操作，∫∫表示二重积分操作，|·|表示取绝对值操作，V_s(t,f)表示信号点序列s在短时傅里叶变换中第t个时刻第f个频率点对应的短时傅里叶变换系数。

遍历窗函数集合{g₁,g₂,…,g_N}中的所有窗口参数，得到局部时刻下对应的瑞利熵值集合

步骤5，获取最佳窗口参数序列。

从瑞利熵集合R′_α中选出最小值，将最小值对应的窗函数作为局部时刻下的最优窗函数。

将短时傅里叶变换所有时刻{t₁,t₂,…,t_N}下，最优窗函数对应的最优窗口参数，组成信号点序列的最佳窗口参数序列{β_t,t＝t₁,t₂,...,t_N}，即不同时刻所加的窗函数是不一样的。

步骤6，获取平滑的最佳窗口参数序列。

用滑动平均法，对最佳窗口参数序列{β_t,t＝t₁,t₂,...,t_N}内的每个点进行累加平均，得到平滑后的最佳窗口参数序列{σ_t,t＝t₁,t₂,...,t_N}。

步骤7，自适应同步压缩短时傅里叶变换。

用平滑的最佳窗口参数序列{σ_t,t＝t₁,t₂,...,t_N}对应的窗函数集合{g_t,t＝t₁,t₂,...,t_N}，对信号点序列进行离散加窗傅傅里叶变换，得到自适应短时傅里叶变换系数。

用维特比Viterbi法，对自适应短时傅里叶变换系数在时频域进行最小路径搜索，得到信号点序列的每个分量的瞬时频率脊线。

首先把所有时刻点的自适应短时傅里叶变换系数按从大到小的顺序排列，计算每个时频点对应的顺序值，确定每个时间点的所有频率点与其前一个时间点的所有频率点的局部最优路径并记录，选取代价函数最小的路径为两个时间点之间的最优路径。

确定每个时间点的所有频率点与后一个时间点的所有频率点的局部最优路径。以两个时间点之间的局部最优路径为基准，向前、向后进行搜索，选取最小代价函数值对应的路径为最优路径，即第一个分量的瞬时频率脊线。

把第一个分量的瞬时频率脊线从自适应短时傅立叶变换系数中清除，再从剩下的自适应短时傅立叶变换系数中重复上述步骤进行搜索，知道找到所有分量的瞬时频率脊线。

按照下式同步压缩公式，根据信号点序列的每个分量的瞬时频率脊线，将每个时刻下所有对应同一频率的自适应短时傅里叶变换系数进行叠加，计算得到自适应同步压缩系数。

其中，T_s(m,n)表示信号点序列s在自适应短时傅里叶变换中的第m个时刻第n个频率点对应的自适应同步压缩系数，π表示圆周率，g表示高斯窗函数，g(0)表示在0时刻的高斯窗函数值，∫表示积分操作符号，V_s(m,n,σ)表示信号点序列s在最佳窗口参数序列σ下，自适应短时傅里叶变换中第m个时刻第n个频率点对应的短时傅里叶变换系数，δ表示单位冲击函数，ω_s(m,n)表示信号点序列s的每个分量在短时傅里叶变换中第m个时刻第n个频率点对应的瞬时频率。

步骤8，获得各分量信号。

按照下式重构公式，将每个时刻下对应同一个瞬时频率脊线的自适应同步压缩系数进行叠加，计算得到信号点序列中的每个分量信号。

其中，s_k表示信号点序列s中的第k个分量信号，≈表示近似等于符号，表示第k个分量信号在自适应短时傅里叶变换中第m个时刻下的瞬时频率值，＜表示小于符号，χ表示门限阈值，1＜χ＜5，T_s(m,n)表示信号点序列s在自适应短时傅里叶变换中的第m个时刻第n个频率点对应的自适应同步压缩系数。

下面结合仿真实验对本发明的效果做进一步的描述。

1、仿真条件：

本发明的仿真信号为两个正余弦调制信号

s(t)＝s₁(t)+s₂(t)

＝cos(2π(80t+sin(12πt))+cos(2π(140t+sin(12πt)))

其参数为：时间T＝1s，采样频率为f_s＝1024Hz。仿真软件环境为Intel(R)Core(TM)i3-4160CPU@3.60GHz双核，Windows 7 32bit操作系统下的Matlab R2012b。

本发明采用重构信号与原始信号的平均相对误差RMSE来评价本发明多分量信号分离的性能，假设Monte Carlo仿真实验的次数为M，原始信号表示为f，重构信号表示为则各分量信号的平均相对误差为

本发明中所加噪声均为加性高斯白噪声，且用信噪比SNR＝10lg(p_s/p_n)来度量信号和噪声能量的大小，其中，p_s表示信号的能量，p_n表示噪声的能量。

2、仿真内容与结果分析

仿真实验1：最佳窗口参数序列图。

本发明是在局部时刻熵值最小的条件下，选取的窗口参数，其对应的窗函数在时频图中充分体现了信号点序列的局部信息。为了表明所有时刻下，窗口参数的变化情况，图2给出信号点序列的最佳窗口参数序列图。

图2中的横坐标表示时间，纵坐标表示窗口参数，图中的曲线表示，窗口参数是随时间不断变化的，即在不同时刻，对信号序列所加的窗函数是不同的。与所有时刻用同一固定窗函数的现有技术相比，本发明更能表现出信号序列的局部信息，同时也提高了信号序列在时频图中的时频分辨率。

仿真实验2：短时傅立叶变换时频图。

本发明得到的最佳窗口参数是随时间变化的，为了凸显本发明方法中变化的窗口参数对信号点序列的分离效果，图3(a)给出最佳窗口序列下，信号点序列的短时傅里叶变化时频图，图3(b)给出最佳窗口参数序列下，信号序列的同步压缩短时傅里叶变换时频图，图3(c)给出固定窗口参数0.006下，信号序列的短时傅里叶变换时频图。

图3(a)中的横坐标表示时间，纵坐标表示频率，图中两条线表示两个分量信号的瞬时频率脊线。由时频图可以看出，信号序列包含两个正余弦调制的信号分量，且在最佳窗口参数序列下，两个分量信号得到很好地分离。

图3(b)中的横坐标表示时间，纵坐标表示频率，图中两条线表示经过同步压缩后，两个分量信号的瞬时频率脊线。由时频图可以看出，本发明得到的两个分量信号的瞬时频率脊线清晰明显，能量聚集性高，信号序列包含的分量完全得到分离。

图3(c)中的横坐标表示时间，纵坐标表示频率。由图看出两分量信号的瞬时频率脊线能量发生严重扩散，导致信号的两分量出现交叉混叠现象，导致预处理后的信号序列包含的分量很难被识别和分离。与图3(a)和图3(b)相比，本发明方法在时频图中，能更好的分离开预处理后信号序列的两个分量。

仿真实验3：各分量信号的误差对比图。

本发明的仿真实验3是当改变信噪比SNR，在信噪比从0db到30db变化时，得到信号点序列的各分量信号与原始信号的平均相对误差。本发明的仿真实验3分别采用本发明的方法与现有技术的经典模式分解EMD法，传统的短时傅里叶变换方法得到的平均相对误差做对比，经过50次Monte Carlo仿真实验，其仿真结果如图4(a)，图4(b)所示。

在图4(a)中的横坐标表示信噪比SNR，纵坐标表示平均相对误差。以“五角星”标示的曲线是经典模式分解法EMD，得到的第一个分量信号的误差曲线，以“正方形”标示的曲线是传统短时傅里叶变换固定窗口参数0.006，得到的第一个信号分量的误差曲线，以“三角形”标示的曲线是同步压缩传统短时傅里叶变换固定窗口参数0.006，得到的第一个信号分量的误差曲线，以“*”标示的曲线是最佳窗口参数自适应短时傅立叶变换，得到的第一个分量信号的误差曲线，以“o”标示的曲线是本发明方法，分离得到的第一个分量信号的误差曲线。

在图4(b)中，横坐标表示信噪比SNR，纵坐标表示平均相对误差。以“五角星”标示的曲线是经典模式分解EMD，得到的第二个分量信号的误差曲线，以“正方形”标示的曲线是传统短时傅里叶变换固定窗口参数0.006下，得到的第二个信号分量的误差曲线，以“三角形”标示的曲线是同步压缩传统短时傅里叶变换固定窗口参数0.006下，得到的第二个信号分量的误差曲线，以“*”标示的曲线是最佳窗口参数自适应短时傅立叶变换下，得到的第二个分量信号的误差曲线，以“o”标示的曲线是本发明方法，得到的第二个分量信号的误差曲线。

由图4(a)和4(b)看出，随着信噪比SNR的不断增加，每个分量信号的平均相对误差都在不断减小，最终达到平衡。但是，同步压缩传统短时傅里叶变换下分离得到的分量信号，要比传统短时傅里叶变换下分离得到的分量信号误差小，本发明方法分离得到的分量信号，要比短时傅里叶变换固定窗口参数和经典模式分解EMD下分离得到的分量信号的误差小。本发明的最佳窗口参数自适应同步压缩短时傅立叶变换下，得到的误差曲线最低，误差最小，分离效果最好的。

以上仿真表明，采用本发明，更能有效的对多分量信号进行分离，且具有很强的鲁棒性。

Claims (4)

1.一种基于瑞利熵的自适应时频同步压缩方法，包括如下步骤：

(1)采集信号：

信号采集系统通过雷达侦察接收机设备，实时采集目标空域中的任意一段信号；

(2)获取信号点序列：

根据所采集信号的长度，对所采集的信号在大于等于1的采样频率范围内，选取一个采样频率进行采样，得信号点序列N＝Tf_s，其中，N表示信号点序列的长度，T表示所采集的雷达侦察接收机信号的时间长度，f_s表示采集信号的采样频率；

(3)选取信号点序列的窗口参数集合：

(3a)选取高斯函数为窗函数；

(3b)在(0,1]区间，以1/N的步长，选取N个不同的值，组成窗口参数集合；

(4)计算瑞利熵：

(4a)在窗口参数集合中，任意选取一个窗口参数，在其对应的窗函数下，对信号点序列进行离散加窗傅里叶变换，得到信号点序列的短时傅里叶变换系数；

(4b)利用瑞利熵公式，得到每个时刻的瑞利熵值；

(4c)遍历窗口参数集合中的所有窗口参数，得到每个时刻对应的瑞利熵值集合；

(5)获取最佳窗口参数序列：

(5a)从瑞利熵值集合中选出最小值，将所选最小值对应的窗函数作为每个时刻的最优窗函数；

(5b)遍历每个时刻，得到所有时刻的最优窗函数，将所有时刻下最优窗函数对应的窗口参数，组成信号点序列的最佳窗口参数序列；

(6)获取平滑的最佳窗口参数序列：

用滑动平均法，对最佳窗口参数序列内的点进行累加平均，得到平滑的最佳窗口参数序列；

(7)自适应同步压缩短时傅里叶变换：

(7a)用平滑的最佳窗口参数序列对应的窗函数集合，对信号点序列进行离散加窗傅里叶变换，得到自适应短时傅里叶变换系数；

(7b)用维特比Viterbi法，对自适应短时傅里叶变换系数在时频域进行最小路径搜索，得到信号点序列的每个分量的瞬时频率脊线；

(7c)利用同步压缩公式，根据信号点序列的每个分量的瞬时频率脊线，将每个时刻下所有对应同一频率的自适应短时傅里叶变换系数进行叠加，得到自适应同步压缩系数；

(8)获得各分量信号：

利用重构公式，将每个时刻下对应信号点序列中同一个分量瞬时频率脊线的自适应同步压缩系数进行叠加，得到信号点序列中的每个分量信号。

2.根据权利要求1所述的基于瑞利熵的自适应时频同步压缩方法，其特征在于：步骤(4b)所述的瑞利熵公式如下：

其中，R表示瑞利熵值，α表示瑞利熵系数，α≥2.5，log₂表示以2为底的对数操作，∫∫表示二重积分操作，|·|表示取绝对值操作，V_s(t,f)表示信号点序列s在短时傅里叶变换中第t个时刻第f个频率点对应的短时傅里叶变换系数。

3.根据权利要求1所述的基于瑞利熵的自适应时频同步压缩方法，其特征在于：步骤(7c)所述的同步压缩公式如下：

其中，T_s(m,n)表示信号点序列s在自适应短时傅里叶变换中的第m个时刻第n个频率点对应的自适应同步压缩系数，π表示圆周率，g表示高斯窗函数，g(0)表示在0时刻的高斯窗函数值，∫表示积分操作，V_s(m,n,σ)表示信号点序列s在平滑后的最佳窗口参数序列σ对应的窗函数集合下，自适应短时傅里叶变换中第m个时刻第n个频率点对应的自适应短时傅里叶变换系数，δ表示单位冲击函数，ω_s(m,n)表示信号点序列s的每个分量在自适应短时傅里叶变换中第m个时刻第n个频率点对应的瞬时频率。

4.根据权利要求1所述的基于瑞利熵的自适应时频同步压缩方法，其特征在于：步骤(8)所述的重构公式如下：

其中，s_k表示信号点序列s中的第k个分量信号，≈表示近似等于符号，表示信号点序列中第k个分量信号在自适应短时傅里叶变换中第m个时刻下的瞬时频率值，＜表示小于符号，χ表示门限阈值，1＜χ＜5，T_s(m,n)表示信号点序列s在自适应短时傅里叶变换中的第m个时刻第n个频率点对应的自适应同步压缩系数。