CN111190157A - 一种ipix雷达回波数据时频分析方法及系统 - Google Patents

Abstract

本发明提供了一种IPIX雷达回波数据时频分析方法及系统，其方法包括：首先计算信号的最优旋转阶数；通过STFrFT方法处理信号，在最优旋转阶数下计算信号的短时分数阶傅里叶变换；计算STFrFT时频谱的瞬时频率；然后由SET技术以STFrFT的处理结果为基础构建同步提取算子SEO，对STFrFT时频谱进行同步提取，获得高聚集度的时频分布。本发明的有益效果是：更加准确的估计信号瞬时频率，得到高聚集度的时间‑频率分布，更为精准地获取由目标运动所产生的频率曲线，区分和辨别运动目标和海杂波；算法创新，效果良好，满足科学分析研究的需要，对雷达目标识别研究以及非平稳信号的研究有积极的意义。

Description

一种IPIX雷达回波数据时频分析方法及系统

技术领域

本发明涉及雷达目标识别研究领域及各种非平稳信号的研究领域，尤其涉及一种IPIX雷达回波数据时频分析方法及系统。

背景技术

雷达目标识别技术是集目标，环境和信号处理技术为一体的一项复杂的系统工程，对于提高军队的指挥自动化水平，攻防能力，国土防空反导能力及战略预警能力具有十分重要的作用。雷达目标识别技术作为现代雷达技术的重要发展方向之一，将成为未来武器系统中的一个重要技术组成部分，具有广阔的应用前景和巨大的军事应用价值，美、俄等军事强国均把它作为发展未来智能化武器系统的重点和需要首先突破的关键技术。

在我国，目前还没有人提出过利用一种高聚集度的时频分析算法来研究海杂波背景下IPIX雷达回波数据，也没有人提出过利用SET-STFrFT来研究IPIX雷达回波数据。

发明内容

为了解决上述问题，本发明提供了一种IPIX雷达回波数据时频分析方法及系统，一种IPIX雷达回波数据时频分析方法，主要包括以下步骤：

S101：获取IPIX雷达回波信号x(t)，并计算信号x(t)的最优旋转阶数；

S102：通过STFrFT方法处理信号x(t)，在最优旋转阶数下计算信号x(t)的短时分数阶傅里叶变换，得到信号x(t)的短时分数阶傅里叶变换时频分布；

S103：根据信号x(t)的短时分数阶傅里叶变换时频分布，计算STFrFT时频谱的瞬时频率；

S104：采用SET技术，以STFrFT时频谱的瞬时频率为基础，构建同步提取算子SEO；

S105：对STFrFT时频谱进行同步提取，获得IPIX雷达回波信号x(t)高聚集度的时频分布。

进一步地，步骤S101中，获取IPIX雷达回波信号x(t)，并计算信号x(t)的最优旋转阶数，具体包括：

S201：计算IPIX雷达回波信号x(t)的信号调频率，记为2a；

S202：对于信号x(t)，采用公式(1)计算信号x(t)的旋转角度α_opt：

S203：根据公式(2)计算信号x(t)的最优旋转阶数p_opt：

进一步地，步骤S102中，在最优旋转阶数下计算信号x(t)的短时分数阶傅里叶变换的方法为：

选择窗函数h(τ-t)，对信号x(t)进行处理；其中，所述窗函数为高斯窗，且长度为预设值；

在最优旋转阶数p_opt下，根据公式(3)计算信号x(t)的短时分数阶傅里叶变换STFrFT_p(t,ω)：

上式中，p为旋转阶数，p≠2n，n为任意整数；α为角度参数，即旋转角度，α与旋转阶数p的关系为α＝pπ/2；K_p(t,u)为核函数，其计算公式如公式(4)所示：

其中，

为分数阶傅里叶变换结果的幅度；u为分数阶傅立叶变换域，简称分数域。

进一步地，步骤S103中，根据STFrFT_p(t,ω)，采用公式(5)计算STFrFT时频谱的瞬时频率ω₀(t,ω)：

进一步地，步骤S104中，采用公式(6)建立同步提取算子SEO：

其中，h′表示窗函数的导数；

表示信号在h′(τ-t)下计算得到的二维矩阵(t,ω)；Δω是离散频率间隔；σ＝Δω/2是判定阈值，数值可根据不同信号在时频表示中的频率间隔设定；Re(g)表示取实部。

进一步地，步骤S105中，利用SEO算子对STFrFT时频谱进行同步提取，提取时频谱中脊线ω＝ω₀位置的时频系数，获得高聚集度的时频分布：

SET_p(t,ω)＝STFrFT_p(t,ω)·SEO(t,ω) (7)。

进一步地，一种IPIX雷达回波数据时频分析系统，其特征在于：包括以下模块：

最优旋转阶数计算模块，用于获取IPIX雷达回波信号x(t)，并计算信号x(t)的最优旋转阶数；

第一处理模块，用于通过STFrFT方法处理信号x(t)，在最优旋转阶数下计算信号x(t)的短时分数阶傅里叶变换，得到信号x(t)的短时分数阶傅里叶变换时频分布；

第二处理模块，用于根据信号x(t)的短时分数阶傅里叶变换时频分布，计算STFrFT时频谱的瞬时频率；

算子提取模块，用于采用SET技术，以STFrFT时频谱的瞬时频率为基础，构建同步提取算子SEO；

时频分布获取模块，用于对STFrFT时频谱进行同步提取，获得IPIX雷达回波信号x(t)高聚集度的时频分布。

进一步地，最优旋转阶数计算模块中，获取IPIX雷达回波信号x(t)，并计算信号x(t)的最优旋转阶数，具体包括以下单元：

信号调频率计算单元，用于计算IPIX雷达回波信号x(t)的调频率，记为2a；

旋转角度计算单元，用于对于信号x(t)，采用公式(8)计算信号x(t)的旋转角度α_opt：

最优旋转阶数计算单元，用于根据公式(9)计算信号x(t)的最优旋转阶数p_opt：

进一步地，第一处理模块中，在最优旋转阶数下计算信号x(t)的短时分数阶傅里叶变换的方法为：

选择窗函数h(τ-t)，对信号x(t)进行处理；其中，所述窗函数为高斯窗，且长度为预设值；

在最优旋转阶数p_opt下，根据公式(10)计算信号x(t)的短时分数阶傅里叶变换STFrFT_p(t,ω)：

上式中，p为旋转阶数，p≠2n，n为任意整数；α为角度参数，即旋转角度，α与旋转阶数p的关系为α＝pπ/2；K_p(t,u)为核函数，其计算公式如公式(11)所示：

其中，

为分数阶傅里叶变换结果的幅度；u为分数阶傅立叶变换域，简称分数域。

进一步地，第二处理模块中，根据STFrFT_p(t,ω)，采用公式(12)计算STFrFT时频谱的瞬时频率ω₀(t,ω)：

进一步地，算子提取模块中，采用公式(13)建立同步提取算子SEO：

其中，h′表示窗函数的导数；

表示信号在h′(τ-t)下计算得到的二维矩阵(t,ω)；Δω是离散频率间隔；σ＝Δω/2是判定阈值，数值可根据不同信号在时频表示中的频率间隔设定；Re(g)表示取实部。

进一步地，时频分布获取模块中，利用SEO算子对STFrFT时频谱进行同步提取，提取时频谱中脊线ω＝ω₀位置的时频系数，获得高聚集度的时频分布：

SET_p(t,ω)＝STFrFT_p(t,ω)·SEO(t,ω) (14)。

本发明提供的技术方案带来的有益效果是：本发明所提出的技术方案可以更加准确的估计信号瞬时频率，得到高聚集度的时间-频率分布，更为精准地获取由目标运动所产生的频率曲线，区分和辨别运动目标和海杂波；算法创新，效果良好，满足科学分析研究的需要，对雷达目标识别研究以及非平稳信号的研究有积极的意义。

附图说明

下面将结合附图及实施例对本发明作进一步说明，附图中：

图1是本发明实施例中一种IPIX雷达回波数据时频分析方法的流程图；

图2是本发明实施例中使用MATLAB程序实现的流程图的示意图；

图3是本发明实施例中单分量调频调幅信号x₁(t)的时频分布的示意图；

图4是本发明实施例中信号x₁(t)不同时频分布的局部放大的示意图；

图5是本发明实施例中单分量非线性调频信号x₂(t)的时频分布示意图；

图6是本发明实施例中信号x₂(t)不同时频分布的局部放大示意图；

图7是本发明实施例中第1个距离单元纯海杂波回波的时频分布示意图；

图8是本发明实施例中第8个距离单元含目标回波的时频分布示意图；

图9是本发明实施例中第10个距离单元受影响回波的时频分布示意图；

图10是本发明实施例中每个距离单元数据基于SET-STFrFT的时频分布的瑞利熵示意图；

图11是本发明实施例中一种IPIX雷达回波数据时频分析系统的模块组成示意图。

具体实施方式

为了对本发明的技术特征、目的和效果有更加清楚的理解，现对照附图详细说明本发明的具体实施方式。

本发明的实施例提供了一种IPIX雷达回波数据时频分析方法及系统。

请参考图1，图1是本发明实施例中一种IPIX雷达回波数据时频分析方法的流程图，具体包括如下步骤：

S101：获取IPIX雷达回波信号x(t)，并计算信号x(t)的最优旋转阶数；

S102：通过STFrFT方法处理信号x(t)，在最优旋转阶数下计算信号x(t)的短时分数阶傅里叶变换，得到信号x(t)的短时分数阶傅里叶变换时频分布；

S103：根据信号x(t)的短时分数阶傅里叶变换时频分布，计算STFrFT时频谱的瞬时频率；

S104：采用SET技术，以STFrFT时频谱的瞬时频率为基础，构建同步提取算子SEO；

S105：对STFrFT时频谱进行同步提取，获得IPIX雷达回波信号x(t)高聚集度的时频分布。

步骤S101中，获取IPIX雷达回波信号x(t)，并计算信号x(t)的最优旋转阶数，具体包括：

S201：计算IPIX雷达回波信号x(t)的信号调频率，记为2a；

S202：对于信号x(t)，采用公式(1)计算信号x(t)的旋转角度α_opt：

S203：根据公式(2)计算信号x(t)的最优旋转阶数p_opt：

步骤S102中，在最优旋转阶数下计算信号x(t)的短时分数阶傅里叶变换的方法为：

选择窗函数h(τ-t)，对信号x(t)进行处理；其中，所述窗函数为高斯窗，且长度为预设值，本发明实施例中，高斯窗的长度为401或者501；

在最优旋转阶数p_opt下，根据公式(3)计算信号x(t)的短时分数阶傅里叶变换STFrFT_p(t,ω)：

上式中，p为旋转阶数，p≠2n，n为任意整数；α为角度参数，即旋转角度，α与旋转阶数p的关系为α＝pπ/2；K_p(t,u)为核函数，其计算公式如公式(4)所示：

其中，

为分数阶傅里叶变换结果的幅度；u为分数阶傅立叶变换域，简称分数域。

步骤S103中，根据STFrFT_p(t,ω)，采用公式(5)计算STFrFT时频谱的瞬时频率ω₀(t,ω)：

步骤S104中，采用公式(6)建立同步提取算子SEO：

其中，h′表示窗函数的导数；

表示信号在h′(τ-t)下计算得到的二维矩阵(t,ω)；Δω是离散频率间隔；σ＝Δω/2是判定阈值，数值可根据不同信号在时频表示中的频率间隔设定；Re(g)表示取实部。

步骤S105中，利用SEO算子对STFrFT时频谱进行同步提取，提取时频谱中脊线ω＝ω₀位置的时频系数，获得高聚集度的时频分布：

SET_p(t,ω)＝STFrFT_p(t,ω)·SEO(t,ω) (7)。

步骤S101中，由于在不同旋转阶数下，信号在短时分数阶傅里叶变换域中呈现出不同的能量聚集的特性。当旋转角度参数与信号的调频率相匹配时，信号能量在相应的分数域表现出最为冲激的效果(即能量最为聚集)，此时的旋转角度达到最优值，用α_opt表示，所对应的p为最优旋转阶数p_opt。在匹配的旋转阶数下，可以获得能量高度聚集的分数谱，最终得到的时频分布图也会具有更高的时频聚集度。

步骤S103中，因为本申请是在传统短时傅里叶变换的基础上推广到分数域，得到短时分数阶傅里叶变换。在传统的短时傅里叶变换中，对信号s(t)的STFT的结果进行偏导计算得：

因为对于任意(t,ω)，STFT_e(t,ω)≠0，所以STFT谱的瞬时频率

同样的，对于STFrFT，只需将STFT变成STFrFT即可，瞬时频率

步骤S104中，在实际的程序实现过程中，将步骤S103中得到的瞬时频率ω₀(t,ω)代入公式

中即可。与步骤S103类似，本申请是类比STFT的SEO推导结果推导了STFrFT的SEO结果，STFT推导过程如下：

(1)最初的建立同步提取算子的公式是：

(2)然后公式(5)已经计算出

需要代入上式计算SEO，但是公式(5)中的分子

还没算出来；

(3)首先计算

在离散计算中，偏导数通常近似由有限差分算子实现，即

(4)为了更准确地估计参数，

可以计算为：

(5)将上式结果代入公式

得到：

(6)将上式结果代入公式

中：

其中，δ(t)是冲击函数，只在t＝0时即δ(0)＝1，t为其他任何不为0的数时，δ(t)＝0；

(7)为了避免程序中遇到0的情况出现计算错误，并考虑到实际应用的情况，可对ω₀(t,ω)只取实部，则

其中，Δω是离散频率间隔；σ＝Δω/2是判定阈值，数值可根据不同信号在时频表示中的频率间隔设定；Re(g)表示取实部。

请参阅图2，图2是本发明所提出的技术方案使用MATLAB程序实现的流程图。

以下为对本发明实施例所提出的技术方案进行时频聚集度仿真验证：

时频聚集度是指在由时频分析方法获得的时频分布中，估计的频率分量在真实频率位置的聚集程度。时频聚集度越高就表明时频分布对于信号频率分量的局部定位和分辨越准确，从而说明该时频分析方法的处理效果越好。

在本发明实施例汇总，将瑞利熵(Rényi entropy)作为衡量时频分布聚集度的性能评价指标，其计算公式为：

其中，α表示瑞利熵的次数，在本发明实施例中，α＝3；TFR(t,f)表示信号的时频分布。瑞利熵Ep的值越小，表示时频分布的聚集度越高。

SET-STFrFT方法可以绘制信号的时-频分布图，能有效的提高时频聚集度，且抗噪声性能比较优越。本发明实施例通过下面的实验验证了本发明所提出的技术方案的优点，设计了多组仿真实验分别加以证明，证明了本发明所提出的技术方案的广泛性以及有效性：

通过下面几组不同仿真信号来对SET-STFrFT与WVD、SST、STFT、STFrFT等方法来进行对比，从而来验证SET-STFrFT方法在处理不同信号时的优越性。

(1)仿真信号x₁(t)：x₁(t)＝(2+sin(2πt))sin(2π(50t²+100t))

信号时间区间为[0s,2s]，采样频率为800。该信号的频率为f₁(t)＝100t+100，中心频率为100Hz，调频率为100Hz/s，幅度呈正弦趋势变化。信号x₁(t)的时域波形以及由五种不同的时频分析方法处理得到的时频分布结果如图3所示，图3中：(a)时域波形；(b)WVD；(c)SST；(d)STFT；(e)STFrFT；(f)SET-STFrFT。通过这五种时频分析方法处理单分量调频调幅信号得到的结果：(a)为信号x₁(t)的时域波形。对于不含噪声的线性调频调幅信号，(b)的WVD和(c)的SST算法都可以得到聚集度很高的时频分布，但WVD有明显的分布不均。图(d)～(e)是本发明所提出的技术方案的整个过程展现，可以看出图(d)的STFT时频图聚集度较差，频率曲线呈带状分布；在最优阶数下图(e)的STFrFT时频图有明显改善，再进行同步提取运算后得到图(f)SET-STFrFT的结果，时频聚集度进一步提高，达到更为理想的效果。

从时频分布的视觉观察上，SST和SET-STFrFT都能表现信号的幅度和频率调制，为了进一步比较两者的时频聚集度，表1列出了不同时频分布图的瑞利熵：

表1信号x₁(t)不同时频分布的瑞利熵

WVD	SST	STFT	STFrFT	SET-STFrFT
					13.8564	12.3040	15.9862	14.6257	12.2121

瑞利熵的值越大表示时频聚集度越低，STFT的瑞利熵值最大，对应在图3中STFT时频图的聚集度最差，而SST和SET-STFrFT的瑞利熵值较小对应两者的时频聚集度较高，也说明了瑞利熵值可以作为评价时频聚集度的性能指标。SET-STFrFT的瑞利熵值为12.2121略小于SST的12.3040，即SET-STFrFT的时频聚集度较高，证明本发明所提出的技术方案在提高聚集度方面有效。

图4给出了WVD、SST和SET-STFrFT时频分布结果的局部放大图像，从细节上也可以看出，WVD、SST在每个时刻的频率分布上仍存在频率模糊现象，而SET-STFrFT有明显的改善。图4中：(a)WVD；(b)SST；(c)SET-STFrFT。

(2)仿真信号x₂(t)：x₂(t)＝sin(-2(t-5)³+200t)

时间区间取为[0s,10s]，采样频率为100。其频率曲线是一个二次函数，表达式为f₂(t)＝(-6(t-5)²+200)/2π，信号的调频率随时间不断变化。图5是基于不同分析方法得出的信号的时频分布结果，对比可以看出SET-STFrFT算法在处理非线性调频信号时的优越性。图5中：(a)时域波形；(b)WVD；(c)SST；(d)STFT；(e)STFrFT；(f)SET-STFrFT。

图5中，(a)为仿真信号x₂(t)的时域波形图。在(b)WVD时频中，虽然信号的频率有较为明显且聚集的轨迹，但在周围有较大区域的模糊，尤其是曲线的顶部，影响频率的识别和分析。SST时频(c)是由STFT结果进行同步压缩处理后得到的，虽然聚集度有所提高，但仍然存在能量表现误差较大的问题。STFT时频(d)中，信号频率的聚集程度不同，在频率曲线的峰值处聚集度较高，此时的调频率较小；但在频率曲线上升和下降的两个阶段聚集度较差，能量较为分散，表现出的能量值(线条颜色)与曲线峰值部分有明显差别。说明当信号的调频率较大，即频率快速变化时STFT的处理效果有待提高。而局部最优的STFrFT可以根据信号在每个窗口中不同的调频率选择不同最优的旋转阶数，使整个信号在不同时段都能得到最高聚集度的时频表示，如(e)所示。本发明所提出的技术方案的处理结果如(f)所示，可以看出，与STFrFT相比，时频聚集度进一步提高，频率曲线轨迹更为理想；并且整段信号的能量表示也较均匀，与信号表达式所定义的幅值恒定的特性相吻合，呈现效果优于(c)SST时频图。

表2中比较了由不同时频分析方法得到的仿真信号x₂(t)的瑞利熵，其中SET-STFrFT算法的瑞利熵值最小，说明其时频聚集度最高，与上述的分析结果一致。

表2信号x₂(t)不同时频分布的瑞利熵

WVD	SST	STFT	STFrFT	SET-STFrFT
					15.4688	11.5378	15.4388	14.4942	11.2388

图6是SST和SET-STFrFT时频分布的局部放大对比，图6中：(a)SST；(b)SET-STFrFT。可以看出，SST中的频率曲线较为模糊，尤其是3s左右曲线展宽严重，且能量值差异较大；而图6的(b)中SET-STFrFT的时频表示在每个时刻的效果相当，能量分布均匀一致，对于信号频率分布的细节描述更为准确清晰。所以，对于非线性调频信号的处理，本发明所提出的技术方案与传统时频分析方法相比，能够得出更好的时频分布图，不仅具有很高的时频聚集度，也能够更为准确理想地反映能量/幅值信息，更适合信号的精准分析。

通过以上数值仿真实验已经验证了本发明所提出的技术方案的优越性，为进一步增加说服力，以下将进一步验证本发明所提出的技术方案在海杂波背景下IPIX雷达回波数据处理研究中的有效性。本发明所提出的技术方案可以较好的得到该信号的时间-频率的联合分布情况，有助于进一步了解海杂波背景下IPIX雷达回波数据的时频分布变化，为雷达海上目标识别研究提供依据和帮助。

在本发明实施例中，选取1993年实测的54#数据(文件名：19931111_163625_starea.cdf)为处理对象，对各距离单元的回波信号进行时频分析，比较和分析纯海杂波回波与含目标回波的时频特征。IPIX雷达实测数据属于典型的非平稳信号，为了获取更多的信息，对54#文件第1、8和9个距离单元雷达实测I通道HH极化数据进行时频分析，得到如图7～图9所示的时频图，从时频联合域分析海杂波及目标的时频特征。图7中：(a)STFT、(b)SET-STFrFT、(c)SST；图8中：(a)STFT、(b)SET-STFrFT、(c)SST；图9中：(a)STFT、(b)SET-STFrFT、(c)SST。

图7为第1个距离单元的回波数据，只有海杂波存在，其能量分布较为分散，频率范围大致在0～200Hz，频率变化没有明显的规律。在图8目标所在单元回波的时频分布中，明显聚集的频率曲线为目标运动所产生的，但是其频率大小相较于海杂波背景的频率范围相对变化不大，说明漂浮目标处于慢速变化的运动中。由于目标的运动，在第10个距离单元的回波数据中(图9)仍可以分辨出较为清晰的目标运动频率曲线，但相对于第8个距离单元，在目标附近存在较低能量的海杂波分布，可判断为受影响的距离单元。

从不同时频分析方法的处理结果比较来看，传统的STFT方法(图7(a)-9(a))可以区分海杂波和漂浮目标，但其时频分布的聚集度较差，尤其是目标运动所产生的频率变化呈带状分布，细节刻画较为粗糙。SST算法(图7(c)-9(c))虽然可以提高时频聚集度，但由于SST对包括海杂波在内的整个时频面进行了压缩，导致目标信号的能量表示产生失真，此问题已在本文的仿真实验部分得到验证。本文所提出的SET-STFrFT方法(图7(b)～9(b))可以很大程度上提高雷达回波数据时频分布的聚集度，准确反应目标信号的能量分布，更容易区别于海杂波背景。

图10为14个距离单元的数据基于SET-STFrFT所得到的时频分布的瑞利熵变化曲线图，可以看出，纯海杂波时频分布的瑞利熵都高于15，目标单元(第8个距离单元)和受影响单元(第7、9、10个距离单元)的瑞利熵明显减小，且目标单元的瑞利熵达到最小值，因此可通过回波数据时频分布的瑞利熵判断目标所在位置。

请参阅图11，图11是本发明实施例中一种IPIX雷达回波数据时频分析系统的模块组成示意图，包括顺次连接的：最优旋转阶数计算模块11、第一处理模块12、第二处理模块13、算子提取模块14和时频分布获取模块15；

最优旋转阶数计算模块11，用于获取IPIX雷达回波信号x(t)，并计算信号x(t)的最优旋转阶数；

第一处理模块12，用于通过STFrFT方法处理信号x(t)，在最优旋转阶数下计算信号x(t)的短时分数阶傅里叶变换，得到信号x(t)的短时分数阶傅里叶变换时频分布；

第二处理模块13，用于根据信号x(t)的短时分数阶傅里叶变换时频分布，计算STFrFT时频谱的瞬时频率；

算子提取模块14，用于采用SET技术，以STFrFT时频谱的瞬时频率为基础，构建同步提取算子SEO；

时频分布获取模块15，用于对STFrFT时频谱进行同步提取，获得IPIX雷达回波信号x(t)高聚集度的时频分布。

最优旋转阶数计算模块11中，获取IPIX雷达回波信号x(t)，并计算信号x(t)的最优旋转阶数，具体包括以下单元：

信号调频率计算单元，用于计算IPIX雷达回波信号x(t)的信号调频率，记为2a；

旋转角度计算单元，用于对于信号x(t)，采用公式(8)计算信号x(t)的旋转角度α_opt：

最优旋转阶数计算单元，用于根据公式(9)计算信号x(t)的最优旋转阶数p_opt：

第一处理模块12中，在最优旋转阶数下计算信号x(t)的短时分数阶傅里叶变换的方法为：

选择窗函数h(τ-t)，对信号x(t)进行处理；其中，所述窗函数为高斯窗，且长度为预设值；

在最优旋转阶数p_opt下，根据公式(10)计算信号x(t)的短时分数阶傅里叶变换STFrFT_p(t,ω)：

上式中，p为旋转阶数，p≠2n，n为任意整数；α为角度参数，即旋转角度，α与旋转阶数p的关系为α＝pπ/2；K_p(t,u)为核函数，其计算公式如公式(11)所示：

其中，

为分数阶傅里叶变换结果的幅度；u为分数阶傅立叶变换域，简称分数域。

第二处理模块13中，根据STFrFT_p(t,ω)，采用公式(12)计算STFrFT时频谱的瞬时频率ω₀(t,ω)：

算子提取模块14中，采用公式(13)建立同步提取算子SEO：

其中，h′表示窗函数的导数；

表示信号在h′(τ-t)下计算得到的二维矩阵(t,ω)；Δω是离散频率间隔；σ＝Δω/2是判定阈值，数值可根据不同信号在时频表示中的频率间隔设定；Re(g)表示取实部。

时频分布获取模块15中，利用SEO算子对STFrFT时频谱进行同步提取，提取时频谱中脊线ω＝ω₀位置的时频系数，获得高聚集度的时频分布：

SET_p(t,ω)＝STFrFT_p(t,ω)·SEO(t,ω) (14)。

本发明的有益效果是：本发明所提出的技术方案进行海杂波背景下IPIX雷达回波数据的时频联合分析，可分析海杂波以及海上漂浮目标的时频特征，更为精准地获取由目标运动所产生的频率曲线，从而有效地进行运动目标和海杂波的区分和辨别。雷达目标识别技术作为现代雷达技术的重要发展方向之一，目前在整个军事雷达研究领域还没有深入的更多研究，对于战略预警能力具有十分重要的作用和意义。介于目前雷达目标识别研究领域没有一个精确的时频分析手段，对于新的雷达信号数据研究很有必要去寻找一个更好的方法。

以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种IPIX雷达回波数据时频分析方法，其特征在于：包括以下步骤：

S101：获取IPIX雷达回波信号x(t)，并计算信号x(t)的最优旋转阶数；

S102：通过STFrFT方法处理信号x(t)，在最优旋转阶数下计算信号x(t)的短时分数阶傅里叶变换，得到信号x(t)的短时分数阶傅里叶变换时频分布；

S103：根据信号x(t)的短时分数阶傅里叶变换时频分布，计算STFrFT时频谱的瞬时频率；

S104：采用SET技术，以STFrFT时频谱的瞬时频率为基础，构建同步提取算子SEO；

S105：对STFrFT时频谱进行同步提取，获得IPIX雷达回波信号x(t)高聚集度的时频分布。

2.如权利要求1所述的一种IPIX雷达回波数据时频分析方法，其特征在于：步骤S101中，获取IPIX雷达回波信号x(t)，并计算信号x(t)的最优旋转阶数，具体包括：

S201：计算IPIX雷达回波信号x(t)的信号调频率，记为2a；

S202：对于信号x(t)，采用公式(1)计算信号x(t)的旋转角度α_opt：

S203：根据公式(2)计算信号x(t)的最优旋转阶数p_opt：

3.如权利要求2所述的一种IPIX雷达回波数据时频分析方法，其特征在于：步骤S102中，在最优旋转阶数下计算信号x(t)的短时分数阶傅里叶变换的方法为：

选择窗函数h(τ-t)，对信号x(t)进行处理；其中，所述窗函数为高斯窗，且长度为预设值；

在最优旋转阶数p_opt下，根据公式(3)计算信号x(t)的短时分数阶傅里叶变换STFrFT_p(t,ω)：

上式中，p为旋转阶数，p≠2n，n为任意整数；α为角度参数，即旋转角度，α与旋转阶数p的关系为α＝pπ/2；K_p(t,u)为核函数，其计算公式如公式(4)所示：

其中，

为分数阶傅里叶变换结果的幅度；u为分数阶傅立叶变换域，简称分数域。

4.如权利要求3所述的一种IPIX雷达回波数据时频分析方法，其特征在于：步骤S103中，根据STFrFT_p(t,ω)，采用公式(5)计算STFrFT时频谱的瞬时频率ω₀(t,ω)：

5.如权利要求4所述的一种IPIX雷达回波数据时频分析方法，其特征在于：步骤S104中，采用公式(6)建立同步提取算子SEO：

上式中，h′表示窗函数的导数；STFrFT_p ^h′(t,ω)表示信号在h′(τ-t)下计算得到的二维矩阵(t,ω)；Δω是离散频率间隔；σ＝Δω/2是判定阈值，根据不同信号在时频表示中的频率间隔设定；Re(g)表示取实部。

6.如权利要求5所述的一种IPIX雷达回波数据时频分析方法，其特征在于：步骤S105中，利用SEO算子对STFrFT时频谱进行同步提取，提取时频谱中脊线ω＝ω₀位置的时频系数，获得高聚集度的时频分布：

SET_p(t,ω)＝STFrFT_p(t,ω)·SEO(t,ω) (7)。

7.一种IPIX雷达回波数据时频分析系统，其特征在于：包括以下模块：

最优旋转阶数计算模块，用于获取IPIX雷达回波信号x(t)，并计算信号x(t)的最优旋转阶数；

第一处理模块，用于通过STFrFT方法处理信号x(t)，在最优旋转阶数下计算信号x(t)的短时分数阶傅里叶变换，得到信号x(t)的短时分数阶傅里叶变换时频分布；

第二处理模块，用于根据信号x(t)的短时分数阶傅里叶变换时频分布，计算STFrFT时频谱的瞬时频率；

算子提取模块，用于采用SET技术，以STFrFT时频谱的瞬时频率为基础，构建同步提取算子SEO；

时频分布获取模块，用于对STFrFT时频谱进行同步提取，获得IPIX雷达回波信号x(t)高聚集度的时频分布。

8.如权利要求7所述的一种IPIX雷达回波数据时频分析方法，其特征在于：最优旋转阶数计算模块中，获取IPIX雷达回波信号x(t)，并计算信号x(t)的最优旋转阶数，具体包括以下单元：

信号调频率计算单元，计算IPIX雷达回波信号x(t)的信号调频率，记为2a；

旋转角度计算单元，用于对于信号x(t)，采用公式(8)计算信号x(t)的旋转角度α_opt：

最优旋转阶数计算单元，用于根据公式(9)计算信号x(t)的最优旋转阶数p_opt：

9.如权利要求8所述的一种IPIX雷达回波数据时频分析方法，其特征在于：第一处理模块中，在最优旋转阶数下计算信号x(t)的短时分数阶傅里叶变换的方法为：

选择窗函数h(τ-t)，对信号x(t)进行处理；其中，所述窗函数为高斯窗，且长度为预设值；

在最优旋转阶数p_opt下，根据公式(10)计算信号x(t)的短时分数阶傅里叶变换STFrFT_p(t,ω)：

上式中，p为旋转阶数，p≠2n，n为任意整数；α为角度参数，即旋转角度，α与旋转阶数p的关系为α＝pπ/2；K_p(t,u)为核函数，其计算公式如公式(11)所示：

其中，

为分数阶傅里叶变换结果的幅度；u为分数阶傅立叶变换域，简称分数域。

10.如权利要求9所述的一种IPIX雷达回波数据时频分析方法，其特征在于：第二处理模块中，根据STFrFT_p(t,ω)，采用公式(12)计算STFrFT时频谱的瞬时频率ω₀(t,ω)：

Images