CN115085831A - 一种基于混合时频分析的无人机遥控信号识别系统和方法 - Google Patents

Abstract

本发明公开了一种基于混合时频分析的无人机遥控信号识别系统及方法，该系统包含：无人机、遥控站和监测装置；该无人机遥控信号识别方法包含：数据采集及合并，采用基于零中频接收机架构的数据采集系统，对无人机信号进行采集；对合并后的时域信号做混合时频变换得到混合时频分析；无人机遥控信号提取以及无人机遥控信号识别；摆脱了单独的线性时频分析或非线性时频分析的局限性，以兼顾时频分辨率和克服干扰；通过提取挠率特征将三维时频曲线的相似性转化为二维挠率分布的距离估计，方法简洁有效，评价的三维时频曲线相似性结果符合感官判断，能够准确地反映三维时频曲线的相似程度。

Description

一种基于混合时频分析的无人机遥控信号识别系统和方法

技术领域

本发明属于无人机遥控信号识别技术领域，具体涉及一种基于混合时频分析的无人机遥控信号识别系统和方法。

背景技术

近年来伴随社会经济与无线网络的发展，无人机在世界各地的使用日益增多，非法无人入侵问题越来越受到各国的关注，如何更快、更精确地实现对无人机的识别是无人机反制的重要前提。现有的通过无线电频谱监测识别无人机的方式，理论上可以对非法无人机连续监测，发现目标并采取反制措施。目前多数的无人机遥控信号采用跳频扩频技术，遥控信号常用的2.4GHz等频段已十分拥堵，又因为其频率跳变的高速性、变频方式的多样性和所在频谱环境的复杂性，使得现有的无线电频谱监测方法往往难以实时、准确地检测到无人机遥控信号，非法无人机监测难以取得令人满意的效果。

文献(李光伟.频谱探测技术在无人机探测与反制领域的应用[J])总结了电磁频谱测量技术在无人机探测和反制中的应用和前景，并给出了几种检测方法，但没有具体算法验证。文献(刘丽.民用无人机跳频信号分析与识别技术研究[D])提出了一种跳频信号的分类系统，可以实现对无人机跳频信号进行自动分类，但是该方法需要大量无人机样本数据，并且算法较复杂。

通常无人机的遥控信号是跳频信号，属于非平稳信号，信号的频率随着时间呈现非线性变化，传统的傅里叶变换不再适用于跳频信号的分析。因此，有必要从时域和频域多角度联合分析，全面展现遥控信号频率随时间变化情况。

发明内容

针对现有技术存在的不足，本发明的目的在于提供一种基于混合时频分析的无人机遥控信号识别系统及方法，混合时频分析通过线性时频分析与非线性时频分析的有效结合，精确的描述时间-频率域来表征并且进行处理和分析信号，摆脱了单一的线性时频分析的时频聚集性低、抗噪能力不佳，及非线性时频分析的交叉项等局限性。

为了实现上述目的，本发明采用以下技术方案予以实现：

一种基于混合时频分析的无人机遥控信号识别方法，包括以下步骤：

步骤一、数据采集及合并：采用基于零中频接收机架构的数据采集系统，对无人机信号进行采集；

接收机内部的模拟/数据器件对接收的模拟信号进行模数转换，并通过正交解调、数字下变频将接收信号转换为同相分量I和正交分量Q两路数字基带信号，然后将I、Q两路信号进行合并，即时域信号sig＝I+j·Q；其中，j表示虚数单位，

步骤二、对合并后的时域信号做混合时频变换得到混合时频分析：

采用混合时频变换对合并后的时域信号进行以时间为横轴，频率为纵轴进行分析，通过对采集的Num_采数据点进行处理，每组数据为Num_组点，共分析Num_采/Num_组次，每次代表一个时间点；Num_组点对应Bw带宽内的信号频率，每一点频率为Bw/Num_组；

混合时频分析方法如下：

S1、将Gabor变换结果TFR_G和乔-威廉姆斯(Choi-Williams distribution，CWD)变换结果TFR_CWD进行哈达玛(Hadamard)乘积操作得到自项时频图；

TFR_自项(t,f)＝TFR_G(t,f)⊙TFR_CWD(t,f)；

式中⊙是哈达玛乘积运算，t为时间变量，f为频率变量；

S2、进行截断处理，得到自项时频区域结果；

式中，| |表示求绝对值；ε为阈值，这里ε＝μ·mean(|TFR_自项(t,f)|)；μ为阈值因子，取μ＝0.15；mean()表示求均值；

S3、将自项时频区域结果与CWD变换结果进行哈达玛乘积，得到混合时频变换的时频图：

TFR_混合(t,f)＝TFR_截断(t,f)⊙TFR_CWD(t,f)；

以上共经过Num_采/Num_组次混合时频变换处理，每一次输出皆可得到混合时频谱；再计算每组数据沿着时间轴的混合时频谱最大值，将各组数据的最大值拼接起来，得到全部样本数据的混合时频谱最大值图；

步骤三、无人机遥控信号提取：低秩矩阵恢复先将时频谱矩阵TFR_混合(t,f)表示为低秩矩阵A与稀疏矩阵E之和，再通过求解范数优化问题来恢复低秩矩阵；

即求解下列最优化问题：

式中，min_A,E( )表示针对低秩矩阵A和稀疏矩阵E求最小值，s.t.表示约束条件，rank( )表示求秩，‖ ‖₀表示求0范数，λ为平衡因子，取

利用快速矩阵分解算法求解(1)式得到稀疏矩阵E；

步骤四、无人机遥控信号识别：经过低秩矩阵恢复方法可以进一步抑制图传信号干扰，在稀疏矩阵E中得到清晰的遥控信号时频表示，从而提取出遥控信号的时频特征，进一步与无人机遥控信号特征库中的图样进行比对；

采用基于挠率特征的三维时频曲线相似性搜索方法，对提取的时频谱特征进行识别，得到该无人机为的确定的无人机型号。

本发明还具有以下技术特征：

优选的，所述的步骤三中最优化问题式(1)的求解方法如下：

a)将最优化问题式(1)正则化，便得到优化问题：

式中，‖ ‖₁为1范数，‖ ‖_1.2为(1,2)范数，

为F范数的平方。其中，μ＝||TFR_混合(t,f)||₂/1.5；

为对不满足线性等式约束的惩罚因子，

是背景先验知识诱导核范数，ρ是低秩矩阵A秩的初始预估计值，这里取ρ＝10；

平衡参数一

平衡参数二β＝0.5×||TFR_混合(t,f)||₂；

b)迭代交替更新矩阵A，E：

当A＝A_k时,

当E＝E_k时,

当达到最大迭代次数10或者误差

时，迭代结束；

则得到稀疏矩阵E，或称为采样数据时频序列。

优选的，所述的步骤四中无人机遥控信号识别，包括以下步骤：

S1)根据遥控信号时频谱数据，拟合遥控信号的三维时频曲线并得到拟合方程；

对于无人机遥控信号标准库中选取目标型号无人机的时频谱数据，等间隔抽取N点时频谱数据，用高阶多项式函数拟合三维曲线，这里取二阶多项式，拟合方程为：

tfr₁(t,f)＝a₁+b₁×sin(m₁×π×t×f)+c₁×exp(-(w₁×f)²) (3)

其中，a₁,b₁,m₁,c₁,w₁为方程(3)的系数，t为时间变量，f为频率变量；

同时，得到采样数据的时频三维拟合曲线，拟合方程为：

tfr₂(t,f)＝a₂+b₂×sin(m₂×π×t×f)+c₂×exp(-(w₂×f)²) (4)

其中，a₂,b₂,m₂,c₂,w₂为方程(4)的系数；

S2)计算三维曲线的挠率：

1)在无人机标准库中选取的无人机型号对应的三维曲线S₁上沿曲线方向等弧长取点集S₁＝{s_1,1,s_1,2,…,s_1,N}，i(i＝1,2,…,N)为点集中数据的序号，N表示总数据数；s_1,i表示时点集S₁中的第i个数据；

2)计算点集S₁中每个元素在曲线上的挠率值，则无人机标准曲线S₁的挠率集合为：

其中,

分别为对应于三维曲线S₁上点s_1,1,s_1,2,…,s_1,N的挠率值；

S3)估计三维时频曲线的挠率分布的切比雪夫距离，判别三维曲线的相似性：切比雪夫距离可用于评估特征空间两个多维向量的差异性，对于任意两条空间曲线，其形状越相近，切比雪夫距离值越小，相似度越高，反之，则表示相似程度低、差异大；

无人机标准库中选取的无人机型号的三维时频曲线S₁的挠率分布与采样数据三维时频曲线S₂的挠率分布间的归一化欧氏距离(Euclidean Distance，ED)为：

当

时，这里取门限值

此时，采样数据三维时频曲线S₂的挠率分布

与无人机标准库中选取的无人机型号的三维时频曲线S₁的挠率分布

相似，可以判断该机型为无人机标准库中选取的无人机型号。

进一步的，所述的挠率值计算方法如下：

将方程(5)改写为参数方程

γ₁(θ)＝(t(θ),f(θ),tfr(θ))；

其中θ为参数，则无人机标准库中选取的无人机型号的三维曲线S₁的挠率为：

其中，γ′₁＝(t′(θ),f′(θ),tfr′(θ))为一阶导组成的三维向量，γ″₁＝(t″(θ),f″(θ),tfr″(θ))为二阶导组成的三维向量，γ″′₁＝(t″′(θ),f″′(θ),tfr″′(θ))为三阶导组成的三维向量，||为求绝对值；×表示外积运算，·表示内积运算。

同时，得到采样数据三维曲线S₂的挠率为：

其中,

分别为对应于采样数据三维曲线S₂上点s_2,1,s_2,2,…,s_2,N的挠率值。

本发明还保护一种采用如上所述的方法的基于混合时频分析的无人机遥控信号识别系统，包括无人机、遥控站和监测装置；

所述的无人机处于飞行状态时，向地面遥控站传输图传信号，所述的遥控站向无人机发射遥控信号，所述的监测装置收集遥控站发射的遥控信号和无人机发射的图传信号，并采用基于混合时频分析的无人机遥控信号识别方法进行识别。

本发明与现有技术相比，具有如下技术效果：

本发明基于混合时频分析的无人机遥控信号识别系统及方法，摆脱了单独的线性时频分析或非线性时频分析的局限性；混合时频分析采用线性时频分析与非线性时频分析的组合，以兼顾时频分辨率和克服干扰,这里采用线性时频分析Garbo变换与非线性时频分析CWD的混合时频分析，即将Garbo变换结果与CWD变换的结果进行哈达玛乘积，二者都存在信号项时,得到哈达玛乘积结果，而交叉项对应的区域为零，从而实现不含交叉项的同时，保持良好的时频聚焦性；

本发明基于低秩矩阵恢复时频滤波方法，采用混合时频变换将无人机信号变换到时频域上，然后依据遥控信号与图传信号在时频域上能量分布的相关性的差异，通过低秩矩阵恢复方法将遥控信号与图传信号分到稀疏矩阵和低秩矩阵中，从而分离遥控信号与图传信号，降低图传信号对遥控信号的干扰；

本发明基于挠率的归一化欧式距离的三维时频曲线相似性方法，以挠率分布作为描述三维时频曲线的特征，采用挠率分布的切比雪夫距离作为衡量曲线的相似度；通过提取挠率特征将三维时频曲线的相似性转化为二维挠率分布的距离估计，方法简洁有效，评价的三维时频曲线相似性结果符合感官判断，能够准确地反映三维时频曲线的相似程度。

附图说明

图1为本发明的基于时频分析的无人机遥控信号识别方法的流程图；

图2为本发明基于挠率特征的三维时频曲线相似性搜索方法的流程图；

图3为本发明基于时频分析的无人机遥控信号识别系统的示意图；

图4为本发明无人机采集信号的混合时频变换的时频谱；

图5为本发明无人机采集信号的混合时频谱在时频平面的等高线图；

图6为本发明无人机遥控信号标准库中大疆精灵4型的时频三维拟合曲线；

图7为本发明采样数据的时频三维拟合曲线。

具体实施方式

下面将对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

实施例1

本实施例给出一种基于混合时频分析的无人机遥控信号识别系统及方法，如图3示，其中无人机向地面遥控站传输图传信号，遥控站向无人机发射遥控信号，监测装置收集遥控信号和图传信号。

如图1所示，为本发明的基于混合时频分析的无人机遥控信号识别方法的流程图，该识别方法包含：

步骤一、数据采集及合并：采用基于AD9361零中频接收机架构的数据采集系统，对大疆某型号无人机信号进行采集；接收机内部的模拟/数据器件对接收的模拟信号进行模数转换，并通过正交解调、数字下变频将接收信号转换为同相分量I和正交分量Q两路数字基带信号，然后将I、Q两路信号进行合并，得到时域信号sig＝I+j·Q，j表示虚数单位，

步骤二、对合并后的时域信号做混合时频变换得到TFR_混合时频分析：由于无人机遥控信号属于非平稳信号，因此采用混合时频变换对合并后的时域信号进行以时间为横轴，频率为纵轴进行分析，通过对采集的8000×2048数据点进行处理，每组数据为2048点，共分析8000次，每次代表一个时间点；2048点对应40MHz带宽内的信号频率，每一点频率为0.019532625MHz；

以上共经过8000次混合时频变换处理，每一次输出皆可得到TFR_混合时频谱；再计算每组数据沿着时间轴的TFR_混合时频谱最大值，将各组数据的最大值拼接起来，得到全部样本数据的TFR_混合时频谱最大值图，如图4、图5所示，t为时间坐标轴，f为频率坐标轴，tfr为时频谱坐标轴。图4为无人机采集信号的混合时频图、图5为无人机采集信号的混合时频图在时频平面的等高线图。

混合时频分析实施方法如下：

S1、将Gabor变换结果TFR_G和乔-威廉姆斯变换结果TFR_CWD进行哈达玛乘积操作得到自项时频图；

TFR_自项(t,f)＝TFR_G(t,f)⊙TFR_CWD(t,f)

式中⊙是哈达玛乘积运算，t为时间变量，f为频率变量。

S2、进行截断处理，得到自项时频区域结果；

式中，| |表示求绝对值；ε为阈值，这里ε＝μ·mean(|TFR_自项(t,f)|)；μ为阈值因子，取μ＝0.15；mean()表示求均值。

S3、将自项时频区域结果与CWD变换结果进行哈达玛乘积，得到混合时频变换的时频图。

TFR_混合(t,f)＝TFR_截断(t,f)⊙TFR_CWD(t,f)

采用线性时频分析与非线性时频分析的组合,以兼顾时频分辨率和克服干扰,这里采用线性时频分析Garbo变换与非线性时频分析CWD的组合时频分析。即将TFR_截断(t,f)与TFR_CWD(t,f)的结果进行哈达玛乘积,二者都存在信号项时,得到哈达玛乘积结果，而交叉项对应的线性变换区域为零，从而实现不含交叉项的同时，保持良好的时频聚焦性。

步骤三、无人机遥控信号提取：根据在时频域上遥控信号与图传信号相关性存在的差异，利用低秩矩阵恢复的方法将相关性较强部分分配到低秩矩阵中，而稀疏矩阵中存放相关性弱的部分。图传信号在带宽范围内是连续存在的，而遥控信号在时频域上是一条频率随时间变化的阶梯状斜(或折)线，因此在分离两者时，大部分的图传信号分在低秩矩阵中，遥控信号被分在稀疏矩阵中，从而达到分离遥控信号与图传信号的目的。

低秩矩阵恢复先将时频谱矩阵TFR_混合(t,f)表示为低秩矩阵A与稀疏(噪声)矩阵E之和，再通过求解范数优化问题来恢复低秩矩阵。

即求解下列最优化问题：

式中，min_A,E( )表示针对低秩矩阵A和稀疏矩阵E求最小值，s.t.表示约束条件，rank()表示求秩，‖‖₀表示求0范数，λ为平衡因子，取

利用快速矩阵分解算法求解(1)式得到稀疏矩阵E。

进一步，求解方法如下：

a)将最优化问题式(1)正则化，便得到优化问题：

式中，‖ ‖₁为1范数，‖ ‖_1.2为(1,2)范数，

为F范数的平方；其中，μ＝||TFR_混合(t,f)||₂/1.5，为对不满足线性等式约束的惩罚因子，

是背景先验知识诱导核范数，ρ是低秩矩阵A秩的初始预估计值，这里取ρ＝10；

平衡参数一

平衡参数二β＝0.5×||TFR_混合(t,f)||₂；

b)迭代交替更新矩阵A，E：

当A＝_k时,

当E＝_k时,

当达到最大迭代次数10或者误差

时，迭代结束。

则得到稀疏矩阵E，或称为采样数据时频曲线,如图7所示。

步骤四、无人机遥控信号识别：经过低秩矩阵恢复方法可以进一步抑制图传信号干扰，在稀疏矩阵E中得到清晰的遥控信号时频表示，从而提取出遥控信号的时频特征，这里如图6采样数据时频序列E的时频图中的三维曲线，进一步与无人机遥控信号标准库中的图样进行比对。

采用基于挠率特征的三维时频曲线相似性搜索方法,对提取的时频谱特征进行识别，如图2所示为基于挠率特征的三维时频曲线相似性搜索方法流程图，得到该无人机为大疆精灵4型。

具体实施方法如下：

S1)根据遥控信号时频谱数据，拟合遥控信号的三维时频曲线并得到拟合方程

如图6为无人机遥控信号标准库中大疆精灵4型的时频三维拟合曲线。从无人机遥控信号标准库大疆精灵4型的时频谱中等间隔抽取214点时频谱数据如图6右下图，用非线性最小平方法拟合三维曲线，拟合结果如图6右上图，图6左图为图6右上图在时间-频率平面的等高线图。

拟合方程为：

tfr₁(t,f)＝0.6386-0.02385sin(0.118π×t×f)-0.1455exp[-(8.78f)²] (3)

类似,得到如图7所示采样数据的时频三维拟合曲线。

拟合方程为：

tfr₂(t,f)＝-0.012-0.0146sin(0.2992π×t×f)+0.0394exp[-(9.348f)²](4)

S2)计算三维时频曲线的挠率

1)在无人机标准库大疆精灵4三维曲线S₁上沿曲线方向等弧长取点集S₁＝{s_1,1,s_1,2,…,s_1,N}，i(i＝1,2,…,N)为点集中数据的序号，N表示总数据数，这里N＝214；s_1,表示时点集S₁中的第i个数据；

2)计算点集S₁中每个元素在曲线上的挠率值，则无人机标准曲线S₁的挠率集合为：

其中,

分别为对应于三维曲线S₁上点s_1,1,s_1,2,…,s_1,N的挠率值。实际上，挠率集合可以视为表征一条曲线特征的多维向量。

进一步，挠率计算方法如下：

将方程(5)改写为参数方程

γ₁(θ)＝(t(θ),f(θ),tfr(θ))

其中θ为参数，则大疆精灵4三维曲线S₁的挠率为：

其中，γ′₁＝(t′(θ),f′(θ),tfr′(θ))为一阶导组成的三维向量，γ″₁＝(t″(θ),f″(θ),tfr″(θ))为二阶导组成的三维向量，γ″₁＝(t″′(θ),f″′(θ),tfr″′(θ))为三阶导组成的三维向量，| |为求绝对值；×表示外积运算，·表示内积运算。

类似的，得到采样数据三维曲线S₂的挠率为：

其中,

分别为对应于采样数据三维曲线S₂上点s_2,1,s_2,2,…,s_2,N的挠率值。

S3)估计三维时频曲线的挠率分布的切比雪夫距离，判别三维曲线的相似性。

切比雪夫距离可用于评估特征空间两个多维向量的差异性，对于任意两条空间曲线，其形状越相近，切比雪夫距离值越小，相似度越高，反之，则表示相似程度低、差异大。

大疆精灵4三维时频曲线S₁的挠率分布

与采样数据三维时频曲线S₂的挠率分布

间的归一化欧氏距离为:

当

时，这里取门限值

此时，采样数据三维时频曲线S₂的挠率分布

与大疆精灵4三维时频曲线S₁的挠率分布

相似，可以判断该机型为大疆精灵4。

尽管本发明的内容已经通过上述优选实施例作了详细介绍，但应当认识到上述的描述不应被认为是对本发明的限制。在本领域技术人员阅读了上述内容后，对于本发明的多种修改和替代都将是显而易见的。因此，本发明的保护范围应由所附的权利要求来限定。

Claims (5)

1.一种基于混合时频分析的无人机遥控信号识别方法，其特征在于，包括以下步骤：

步骤一、数据采集及合并：采用基于零中频接收机架构的数据采集系统，对无人机信号进行采集；

接收机内部的模拟/数据器件对接收的模拟信号进行模数转换，并通过正交解调、数字下变频将接收信号转换为同相分量I和正交分量Q两路数字基带信号，然后将I、Q两路信号进行合并，即时域信号sig＝I+j·Q；其中，j表示虚数单位，

步骤二、对合并后的时域信号做混合时频变换得到混合时频分析：

采用混合时频变换对合并后的时域信号进行以时间为横轴，频率为纵轴进行分析，通过对采集的Num_采数据点进行处理，每组数据为Num_组点，共分析Num_采/Num_组次，每次代表一个时间点；Num_组点对应Bw带宽内的信号频率，每一点频率为Bw/Num_组；

混合时频分析方法如下：

S1、将Gabor变换结果TFR_G和乔-威廉姆斯(Choi-Williams distribution，CWD)变换结果TFR_CWD进行哈达玛(Hadamard)乘积操作得到自项时频图；

TFR_自项(t,f)＝TFR_G(t,f)⊙TFR_CWD(t,f)；

式中⊙是哈达玛乘积运算，t为时间变量，f为频率变量；

S2、进行截断处理，得到自项时频区域结果；

式中，||表示求绝对值；ε为阈值，这里ε＝μ·mean(|TFR_自项(t,f)|)；μ为阈值因子，取μ＝0.15；mean()表示求均值；

S3、将自项时频区域结果与CWD变换结果进行哈达玛乘积，得到混合时频变换的时频图：

TFR_混合(t,f)＝TFR_截断(t,f)⊙TFR_CWD(t,f)；

以上共经过Num_采/Num_组次混合时频变换处理，每一次输出皆可得到混合时频谱；再计算每组数据沿着时间轴的混合时频谱最大值，将各组数据的最大值拼接起来，得到全部样本数据的混合时频谱最大值图；

步骤三、无人机遥控信号提取：低秩矩阵恢复先将时频谱矩阵TFR_混合(t,f)表示为低秩矩阵A与稀疏矩阵E之和，再通过求解范数优化问题来恢复低秩矩阵；

即求解下列最优化问题：

式中，min_A,E( )表示针对低秩矩阵A和稀疏矩阵E求最小值，s.t.表示约束条件，rank()表示求秩，‖ ‖₀表示求0范数，λ为平衡因子，取

利用快速矩阵分解算法求解(1)式得到稀疏矩阵E；

步骤四、无人机遥控信号识别：经过低秩矩阵恢复方法可以进一步抑制图传信号干扰，在稀疏矩阵E中得到清晰的遥控信号时频表示，从而提取出遥控信号的时频特征，进一步与无人机遥控信号特征库中的图样进行比对；

采用基于挠率特征的三维时频曲线相似性搜索方法，对提取的时频谱特征进行识别，得到该无人机为的确定的无人机型号。

2.如权利要求1所述的基于混合时频分析的无人机遥控信号识别方法，其特征在于，所述的步骤三中最优化问题式(1)的求解方法如下：

a)将最优化问题式(1)正则化，便得到优化问题：

式中，‖‖₁为1范数，‖‖_1.2为(1,2)范数，

为F范数的平方；其中，μ＝||TFR_混合(t,f)||₂/1.5；

为对不满足线性等式约束的惩罚因子，

是背景先验知识诱导核范数，ρ是低秩矩阵A秩的初始预估计值，这里取ρ＝10；

平衡参数二β＝0.5×||TFR_混合(t,f)||₂；

b)迭代交替更新矩阵A，E：

当A＝A_k时,

当E＝E_k时,

当达到最大迭代次数10或者误差

时，迭代结束；

则得到稀疏矩阵E，或称为采样数据时频序列。

3.如权利要求1所述的基于混合时频分析的无人机遥控信号识别方法，其特征在于，所述的步骤四中无人机遥控信号识别，包括以下步骤：

S1)根据遥控信号时频谱数据，拟合遥控信号的三维时频曲线并得到拟合方程

于无人机遥控信号标准库中选取目标型号无人机的时频谱数据，等间隔抽取N点时频谱数据，用高阶多项式函数拟合三维曲线，这里取二阶多项式，拟合方程为：

tfr₁(t,f)＝a₁+b₁×sin(m₁×π×t×f)+c₁×exp(-(w₁×f²) (3)

其中，a₁,b₁,m₁,c₁,w₁为方程(3)的系数，t为时间变量，f为频率变量；

同时，得到采样数据的时频三维拟合曲线，拟合方程为：

tfr₂(t,f)＝a₂+b₂×sin(m₂×π×t×f)+c₂×exp(-(w₂×f)²) (4)

其中，a₂,b₂,m₂,c₂,w₂为方程(4)的系数；

S2)计算三维曲线的挠率：

1)在无人机标准库中选取的无人机型号对应的三维曲线S₁上沿曲线方向等弧长取点集S₁＝{s_1,1,s_1,2,…,s_1,N}，i(i＝1,2,…,N)为点集中数据的序号，N表示总数据数；s_1,i表示时点集S₁中的第i个数据；

2)计算点集S₁中每个元素在曲线上的挠率值，则无人机标准曲线S₁的挠率集合为：

其中,

分别为对应于三维曲线S₁上点s_1,1,s_12,,…,s_1,N的挠率值；

S3)估计三维时频曲线的挠率分布的切比雪夫距离，判别三维曲线的相似性：切比雪夫距离可用于评估特征空间两个多维向量的差异性，对于任意两条空间曲线，其形状越相近，切比雪夫距离值越小，相似度越高，反之，则表示相似程度低、差异大；

无人机标准库中选取的无人机型号的三维时频曲线S₁的挠率分布与采样数据三维时频曲线S₂的挠率分布间的归一化欧氏距离(Euclidean Distance，ED)为:

当

时，这里取门限值

此时，采样数据三维时频曲线S₂的挠率分布

与无人机标准库中选取的无人机型号的三维时频曲线S₁的挠率分布

相似，可以判断该机型为无人机标准库中选取的无人机型号。

4.如权利要求3所述的基于混合时频分析的无人机遥控信号识别方法，其特征在于，所述的挠率值计算方法如下：

将方程(5)改写为参数方程

γ₁(θ)＝(t(θ),f(θ),tfr(θ))

其中θ为参数，则无人机标准库中选取的无人机型号的三维曲线S₁的挠率为：

其中，γ′₁＝(t′(θ),f′(θ),tfr′(θ))为一阶导组成的三维向量，γ₁″＝(t″(θ),f″(θ),tfr″(θ))为二阶导组成的三维向量，γ₁″′＝(t″′(θ),f″′(θ),tfr″′(θ))为三阶导组成的三维向量，| |为求绝对值；×表示外积运算，·表示内积运算。

同时，得到采样数据三维曲线S₂的挠率为：

其中,

分别为对应于采样数据三维曲线S₂上点s_2,1,s_2,2,…,s_2,N的挠率值。

5.一种采用权利要求1至4中任一项所述的方法的基于混合时频分析的无人机遥控信号识别系统，其特征在于，包括无人机、遥控站和监测装置；

所述的无人机处于飞行状态时，向地面遥控站传输图传信号，所述的遥控站向无人机发射遥控信号，所述的监测装置收集遥控站发射的遥控信号和无人机发射的图传信号，并采用基于混合时频分析的无人机遥控信号识别方法进行识别。

Images