基于雷达目标距离像时频特征提取的雷达目标识别方法

摘要

本发明提供了一种基于雷达目标距离像时频特征提取的雷达目标识别方法，该方法首先采用匹配追踪时频分析方法计算得到雷达目标距离像的时频分布矩阵，然后应用非负矩阵分解技术分解时频分布矩阵，得到能够有效反映雷达目标HRRP内部时频相干结构的时频特征；最后，再根据雷达目标距离像的时频特征提取方案进行雷达目标识别；该方法能够通过计算机运行执行，减少雷达目标识别观测人员的工作量，有助于提高雷达目标识别效率，并且该方法中采用非负矩阵分解有效实现了目标距离像时频矩阵的降维，并分解得到雷达目标散射的时相结构，能够进一步提高雷达目标识别的运算性能和准确性。

说明书

技术领域

本发明属于雷达目标识别技术领域，尤其涉及一种基于雷达目标距离像时频特征提取 的雷达目标识别方法。

背景技术

雷达目标识别是国防科技领域中的重要研究课题。利用雷达目标距离像进行雷达目标 识别是主要途径，其中的关键是雷达目标距离像的有效特征提取。

雷达目标识别是国防科技领域中的重要研究课题，而利用雷达目标距离像进行雷达目 标识别又是雷达目标识别的主要实现途径。利用雷达目标距离像进行雷达目标识别的具体 方法通常分为两大步骤：步骤一，确定雷达目标距离像的目标特征提取方案，即选定雷达 目标距离像中一类可识别特征作为目标特征，并建立该目标特征的提取方案；步骤二，将 根据前述方案提取的雷达目标距离像的目标特征作为雷达目标识别特征，进行雷达目标识 别，即采集已知雷达目标的雷达目标距离像，采用步骤一所述的目标特征的提取方案提取 出已知雷达目标的雷达目标距离像对应的目标特征作为已知雷达目标对应的目标特征模 板，然后再采用步骤一所述的目标特征的提取方案对待测雷达目标的雷达目标距离像提取 出待识别的目标特征，而后将该待识别的目标特征与各个目标特征模板进行匹配，若该待 识别的目标特征存在相匹配的目标特征模板，则判定待测雷达目标与所述相匹配的目标特 征模板对应的已知雷达目标属于同类雷达目标，即完成雷达目标识别。通过上述的雷达目 标识别过程可以看出，如何确定雷达目标距离像的目标特征提取方案，是影响雷达目标识 别准确性的一个关键因素，若目标特征自身选择不准确，或者目标特征提取过程使得提取 结果不够准确，都会影响雷达目标识别的准确性。

目前，常见的雷达目标距离像的目标特征提取方法有：(1)采用雷达目标距离像的功率 谱作为目标特征，但功率谱只能反映雷达目标距离像频谱的全局信息，不能描述雷达目标 距离像的局部频率特性，导致雷达目标识别的分类能力有限，限制了雷达目标识别的范围 和准确性。(2)采用Relax方法提取目标距离像的散射中心幅度、位置参数作为目标特征， 但是这种目标特征对噪声和雷达目标闪烁都很敏感，而噪声和雷达目标闪烁则容易影响雷 达目标识别的准确性。(3)采用时频分析提取雷达目标距离像的时频特征作为目标特征，由 于雷达目标距离像反映了目标散射的物理机理，而复杂目标的散射物理机制较为复杂，例 如目标上某些部件的散射并不局限在某个时间区域上，而可能会引起回波在距离向上的扩 展，这些散射现象仅从散射中心来描述是不够的，而时频域分析提取的取雷达目标距离像 的时频特征可以有效地揭示并描述这些复杂的电磁散射特性信息，有助于提高雷达目标识 别的准确性，因此非常适宜作为雷达目标距离像的目标特征。然而，雷达目标距离像的时 频平面分布数据直接作为目标特征用于雷达目标识别的分类较为困难，因为雷达目标距离 像的时频分布数据矩阵的维数很高，直接作为目标特征将导致数据维数太大，数据处理过 于复杂，而且其中能够用于有效识别的特性信息只分布在时频平面上的少数局部位置，而 其他相当多的数据元素对目标识别是无用的，因此还需要对雷达目标距离像的时频平面分 布数据进行进一步的时频特征提取，一方面从中提取能够用于有效识别的特性信息，另一 方面以降低数据维数，从能够更有效的进行雷达目标的识别。但在现有技术中，针对雷达 目标距离像的时频特征提取，通常都由人工从时频图谱上进行观察识别后加以提取，这种 人工观察识别提取方法不仅非常耗时，并且由于不同观测技术人员的认知差别，其识别提 取标准参差不齐，导致雷达目标距离像的时频特征提取效率低、雷达目标识别准确性难以 保证。

发明内容

针对现有技术中存在的上述问题，为了解决现有技术中雷达目标识别效率低、准确性 不高的问题，本发明提供了一种基于雷达目标距离像时频特征提取的雷达目标识别方法， 该雷达目标识别方法能够通过计算机运行执行，有助于减少雷达目标识别观测人员的工作 量，提高识别效率，更重要的是该方法能够有效提高雷达目标识别的准确性。

为实现上述目的，本发明采用了如下技术手段：

基于雷达目标距离像时频特征提取的雷达目标识别方法，包括如下步骤：

A）建立包含有X个基函数的基函数库G_D；其中任意第x个基函数G_x(t)∈G_D，其表 达式为：

$G_x\left(t\right)=\frac{1}{\sqrt{s_x}}g\left(\frac{t-p_x}{s_x}\right)\exp \left[j2(\pi f_{x}t+{\phi }_x)\right];$

其中，t表示时间；s_x、p_x、f_x和φ_x分别表示第x个基函数G_x(t)的比例控制参数、时 域位置控制参数、频率参数和相角参数；为高斯函数g(t)的时移函数；

B）提取雷达目标距离像的时频特征；该步骤具体为：

b1）用基函数库G_D中的基函数对雷达目标距离像Θ(t)进行匹配追踪分解，得到其分 解表达式：

$\Theta \left(t\right)=\underset{i=1}{\overset{I}{\Sigma }}{a}_i{G}_i\left(t\right)+R_I;$

其中，I表示匹配追踪分解的级数，且I＜X，i∈{1,2,…,I}；R_I表示雷达目标距离像Θ(t) 在I级匹配追踪分解后的残余信号；G_i(t)表示第i级分解所匹配的基函数，且G_i(t)∈G_D； a_i表示第i级分解的补偿因子，且a_i=|<R_i-1，G_i(t)>|，其中R_i-1雷达目标距离像Θ(t)在第(i-1) 级匹配追踪分解后的残余信号，|<R_i-1，G_i(t)>|表示取R_i-1与G_i(t)进行内积运算后的绝对值；

b2）取雷达目标距离像Θ(t)分解表达式的匹配追踪分解组成部分作为雷达目 标距离像Θ(t)对应的时频函数Γ(t)，即令求取所述时频函数Γ(t)的时频 分布函数Ω(t,f)：

$\Omega (t,f)=\underset{i=1}{\overset{I}{\Sigma }}{\left|a_i\right|}^2\mathrm{WV}{G}_i(t,f)=\underset{i=1}{\overset{I}{\Sigma }}{\left|a_i\right|}^2{\hat{g}}^2\left(s_i(f-f_i)\right)g^2\left(\frac{t-p_i}{s_i}\right);$

其中，f表示频率；WVG_i(t,f)表示基函数G_i(t)的维格纳-威利分布函数；s_i、p_i和f_i分别表示基函数G_i(t)的比例控制参数、时域位置控制参数和频率参数；为高斯 函数g(t)的时移函数；为所述时移函数经傅里叶变换得到的频域函 数；

b3）根据雷达目标距离像Θ(t)的时域采样点总数N以及预设定的频域采样点总数M， N≤M≤2N，求得时频分布函数Ω(t,f)对应的维数为M行×N列的时频分布矩阵V_Ω(n,m)：

$V_{\Omega }(n,m)=\underset{i=1}{\overset{I}{\Sigma }}{\left|a_i\right|}^2\left(\begin{array}{c}{\hat{g}}^2\left(s_i(f_1-f_i)\right)\\ {\hat{g}}^2\left(s_i(f_2-f_i)\right)\\ .\\ .\\ .\\ {\hat{g}}^2\left(s_i(f_m-f_i)\right)\\ .\\ .\\ .\\ {\hat{g}}^2\left(s_i(f_M-f_i)\right)\end{array}\right)·\left(\begin{array}{cccccc}{g}^2\left(\frac{t_1-p_i}{s_i}\right)& g^2\left(\frac{t_2-p_i}{s_i}\right)& ...& g^2\left(\frac{t_n-p_i}{s_i}\right)& ...& g^2\left(\frac{t_N-p_i}{s_i}\right)\end{array}\right);$

其中，t_n表示时频分布函数Ω(t,f)在时域上第n个采样点对应的时间，n∈{1,2,…,N}； f_m表示时频分布函数Ω(t,f)在频域上第m个采样点对应的频率，m∈{1,2,…,M}；

b4）根据预设定的分解维度K和分解门限阈值ε，K＜N，0<ε≤10^-2，对所述时频分布矩 阵V_Ω(n,m)进行非负矩阵分解，获得维数为M行×K列的基元矩阵W_Ω(k,m)和维数为K行×N 列的系数矩阵H_Ω(n,k)，使其满足欧氏距离平方值E(W_Ω,H_Ω)<ε；所述欧氏距离平方值 E(W_Ω,H_Ω)表示时频分布矩阵V_Ω(n,m)与矩阵乘积W_Ω(k,m)·H_Ω(n,k)之间的欧氏距离的平 方，即：

E(W_Ω,H_Ω)=||V_Ω(n,m)-W_Ω(k,m)·H_Ω(n,k)||²；

其中，基元矩阵W_Ω(k,m)和系数矩阵H_Ω(n,k)均为非负矩阵；k∈{1,2,…,K}； b5）获取基元矩阵W_Ω(k,m)的各个列向量以及系数矩阵H_Ω(n,k)的各个行向量 k∈{1,2,…,K}；

${\overrightarrow{w}}_k\left(m\right)=\left(\begin{array}{c}{w}_{1,k}\\ w_{2,k}\\ .\\ .\\ .\\ w_{m,k}\\ .\\ .\\ .\\ w_{M,k}\end{array}\right);$

${\overrightarrow{h}}_k\left(n\right)=\left(\begin{array}{cccccc}{h}_{k,1}& h_{k,2}& ...& h_{k,n}& ...& h_{k,N}\end{array}\right);$

其中，w_m，k表示基元矩阵W_Ω(k,m)中第m行第k列的元素；h_k,n表示系数矩阵H_Ω(n,k) 中第k行第n列的元素；

b6）根据预设定的频域特征阶数Q和时域特征阶数P，计算基元矩阵W_Ω(k,m)各个列 向量对应的频域稀疏特征和各阶频域向量矩以及系数矩阵H_Ω(n,k)各 个行向量对应的时域稀疏特征和各阶时域向量矩

${\mathrm{SF}}_{{\overrightarrow{w}}_k}=\frac{\sqrt{M}-[\left(\underset{m=1}{\overset{M}{\Sigma }}{w}_{m,k}\right)/\sqrt{\underset{m=1}{\overset{M}{\Sigma }}{w_{m,k}}^2}]}{\sqrt{M}-1};$

$W{F}_{{\overrightarrow{w}}_k}\left(q\right)=\underset{m=1}{\overset{M}{\Sigma }}{(m-{\mu }_{{\overrightarrow{w}}_k})}^q{w}_{m,k};$

${\mathrm{ST}}_{{\overrightarrow{h}}_k}=\frac{\sqrt{N}-[\left(\underset{n=1}{\overset{N}{\Sigma }}{h}_{k,n}\right)/\sqrt{\underset{n=1}{\overset{N}{\Sigma }}{w_{k,n}}^2}]}{\sqrt{N}-1};$

$M{T}_{{\overrightarrow{h}}_k}\left(p\right)=\underset{n=1}{\overset{N}{\Sigma }}{(n-{\mu }_{{\overrightarrow{h}}_k})}^p{w}_{k,n}$

其中，k∈{1,2,…,K}；表示基元矩阵W_Ω(k,m)中第k列的列向量对应 的第q阶的频域向量矩，q∈{1,2,…,Q}；表示列向量中各个元素的均值，即

表示系数矩阵H_Ω(n，k)中第k行的行向量对应的第p阶的 时域向量矩，p∈{1,2,…,P}；表示行向量中各个元素的均值，即

b7）将基元矩阵W_Ω(k,m)中第k列的列向量对应的频域稀疏特征和各阶频 域向量矩以及系数矩阵H_Ω(n,k)中第k行的行向量对应的时域稀疏特征和各阶时域向量矩的集合作为雷达目标距离像Θ(t)的第k个时频特征矢量ξ_k，即：

${\xi }_k=[{\mathrm{MT}}_{{\overrightarrow{h}}_k}\left(1\right),{\mathrm{MT}}_{{\overrightarrow{h}}_k}\left(2\right),...,{\mathrm{MT}}_{{\overrightarrow{h}}_k}\left(p\right),...,{\mathrm{MT}}_{{\overrightarrow{h}}_k}\left(p\right),$

${\mathrm{MF}}_{{\overrightarrow{w}}_k}\left(1\right),{\mathrm{MF}}_{{\overrightarrow{w}}_k}\left(2\right),...,{\mathrm{MF}}_{{\overrightarrow{w}}_k}\left(q\right),...,{\mathrm{MF}}_{{\overrightarrow{w}}_k}\left(Q\right),;$

${\mathrm{ST}}_{{\overrightarrow{h}}_k},{\mathrm{SF}}_{{\overrightarrow{w}}_k}]$

由此，得到雷达目标距离像Θ(t)的时频特征矢量序列ζ={ξ₁,ξ₂,…,ξ_k,…,ξ_K}，并将所 述时频特征矢量序列ζ作为从雷达目标距离像Θ(t)提取的时频特征；

C）将根据步骤B）所述方案提取的雷达目标距离像的时频特征作为雷达目标识别特 征，进行雷达目标识别。

上述的雷达目标识别方法中，进一步，所述匹配追踪分解的级数I的取值范围为 20≤I≤100。

上述的雷达目标识别方法中，进一步，所述分解维度K的取值范围为0.05N≤K≤0.5N。

上述的雷达目标识别方法中，进一步，所述频域特征阶数Q的取值范围为1≤Q≤6；所 述时域特征阶数P的取值范围为1≤P≤6。

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

1、本发明提供了一种基于雷达目标距离像时频特征提取的雷达目标识别方法，该方 法能够通过计算机运行执行，减少雷达目标识别观测人员的工作量，并且有助于提高雷达 目标识别效率。

2、在本发明基于雷达目标距离像时频特征提取的雷达目标识别方法中，基于稀疏表 示的雷达目标距离像时频矩阵计算具有无交叉项干扰、时频分辨率高等优点。

3、在本发明基于雷达目标距离像时频特征提取的雷达目标识别方法中，采用非负矩 阵分解有效实现了目标距离像时频矩阵的降维，以提高运算性能；另外，分解得到的基矩 阵矢量表征了目标散射的频谱结构，分解得到的系数矩阵矢量表征了雷达目标散射的时相 结构，这样可以有效提取时频矩阵中的显著特征，更有助于时频特征的描述以及雷达目标 的识别，能够进一步提高雷达目标识别的准确性。

附图说明

图1为本发明基于雷达目标距离像时频特征提取的雷达目标识别方法的流程框图；

图2为本发明实施例中五类雷达目标的光学图像；

图3为本发明实施例中一个BMP2雷达目标信号图像的雷达目标距离像序列数据中第 1个雷达目标距离像图；

图4为本发明实施例中一个BMP2雷达目标信号图像的雷达目标距离像序列数据中第 50个雷达目标距离像图；

图5为本发明实施例中一个BMP2雷达目标信号图像的雷达目标距离像序列数据中第 100个雷达目标距离像图；

图6为图3所示的雷达目标距离像的时频分布矩阵图形；

图7为图6所示时频分布矩阵经非负矩阵分解得到的基元矩阵W_Ω(k,m)的10个列向量 图；

图8为图6所示时频分布矩阵经非负矩阵分解得到的系数矩阵H_Ω(n,k)的10个行向量 图；

图9为本发明实施例中五类雷达目标的雷达目标距离像时频特征分布图。

具体实施方式

下面结合附图和实施例对本发明的技术方案作进一步的说明。

本发明提出了一种基于雷达目标距离像时频特征提取的雷达目标识别方法，该方法首 先采用匹配追踪时频分析方法计算得到雷达目标HRRP（High Resolution Range Profile，雷 达目标高分辨一维距离像，通常简称为雷达目标距离像）的时频分布矩阵，然后应用非负 矩阵分解技术分解时频分布矩阵，得到能够有效反映雷达目标HRRP内部时频相干结构的 时频特征；最后，再根据雷达目标HRRP的时频特征提取方案进行雷达目标识别。本发明 基于雷达目标距离像时频特征提取的雷达目标识别方法，该方法的流程如图1所示，具体 包括如下步骤：

A）建立包含有X个基函数的基函数库G_D；其中任意第x个基函数G_x(t)∈G_D，其表 达式为：

$G_x\left(t\right)=\frac{1}{\sqrt{s_x}}g\left(\frac{t-p_x}{s_x}\right)\exp \left[j2(\pi f_{x}t+{\phi }_x)\right];$

其中，t表示时间；s_x、p_x、f_x和φ_x分别表示第x个基函数G_x(t)的比例控制参数、时 域位置控制参数、频率参数和相角参数；为高斯函数g(t)的时移函数。

该步骤用于建立基函数库，而建立基函数库的目的在于，利用基函数库中的基函数， 对雷达目标距离像进行匹配追踪分解，以完成对雷达目标距离像的分解表达。在实际应用 中，基函数库应当是一个完备和冗余的基函数的集合；这里的“完备和冗余”可以理解为， 基函数库G_D中所包含的基函数个数X，要远远的大于能够完成雷达目标距离像的分解表达 所需要的基函数个数；因此，基函数库的建立可以在实际应用中不断积累完成，基函数库 中包含的基函数越多，越有利于进行对雷达目标距离像的分解表达，通常基函数库中应当 包含有数以万计的基函数。而基函数以高斯函数g(t)作为基础，目的在于利用高斯函数的 特性体现时频特征。

B）提取雷达目标距离像的时频特征。该步骤的处理过程是本发明的重要发明点，其 处理过程具体为：

b1）用基函数库G_D中的基函数对雷达目标距离像Θ(t)进行匹配追踪分解，得到其分 解表达式：

$\Theta \left(t\right)=\underset{i=1}{\overset{I}{\Sigma }}{a}_i{G}_i\left(t\right)+R_I;$

其中，I表示匹配追踪分解的级数，且I＜X，i∈{1,2,…,I}；R_I表示雷达目标距离像Θ(t) 在I级匹配追踪分解后的残余信号；G_i(t)表示第i级分解所匹配的基函数，且G_i(t)∈g_D； a_i表示第i级分解的补偿因子，且a_i=|<R_i-1,G_i(t)>|，其中R_i-1雷达目标距离像Θ(t)在第(i-1) 级匹配追踪分解后的残余信号，|<R_i-1，G_i(t)>|表示取R_i-1与G_i(t)进行内积运算后的绝对值。

该步骤是通过匹配追踪分解完成雷达目标距离像Θ(t)的分解表达，所获得分解表达式 包括两部分，即匹配追踪分解组成部分和残余信号部分R_I，只要匹配追踪分解的 级数I足够大，就能够将雷达目标距离像中的主要时频相干结构分解到匹配追踪分解组成 部分中，而残余信号R_I则可以看做是随机的噪声或杂波，其中并没有体现出足够 的时频局部化特点，可以在后续的时频体征提取处理中把残余信号R_I忽略。为了能够取得 这样的分解效果，匹配追踪分解的级数I至少应该为20级，匹配追踪分解的级数越多则分 解效果越好；但匹配追踪分解的级数若过多，例如超过100级以后，分解效果的提升就不 太明显了，而且过多的分解技术会增加数据处理量；因此，匹配追踪分解的级数I的优选 取值范围为20≤I≤100。至于匹配追踪分解方法，是现有的成熟算法，其具体过程本文就不 再赘述了。

b2）取雷达目标距离像Θ(t)分解表达式的匹配追踪分解组成部分作为雷达目 标距离像Θ(t)对应的时频函数Γ(t)，即令求取所述时频函数Γ(t)的时频 分布函数Ω(t,f)：

$\Omega (t,f)=\underset{i=1}{\overset{I}{\Sigma }}{\left|a_i\right|}^2\mathrm{WV}{G}_i(t,f)=\underset{i=1}{\overset{I}{\Sigma }}{\left|a_i\right|}^2{\hat{g}}^2\left(s_i(f-f_i)\right)g^2\left(\frac{t-p_i}{s_i}\right);$

其中，f表示频率；WVG_i(t,f)表示基函数G_i(t)的维格纳-威利分布函数；s_i、p_i和f_i分别表示基函数G_i(t)的比例控制参数、时域位置控制参数和频率参数；为高斯 函数g(t)的时移函数；为所述时移函数经傅里叶变换得到的频域函 数。

该步骤中，由于雷达目标距离像分解表达式的匹配追踪分解组成部分包含了 雷达目标距离像中的主要时频相干结构，因此将该部分作为时频函数Γ(t)并进行维格纳- 威利分布，以获得时频分布函数Ω(t,f)。根据维格纳-威利分布的特性，对于高斯函数g(t) 的维格纳-威利分布WVg(t,f)，是一个二维高斯函数，即其中的 为高斯函数g(t)的傅里叶变换函数；由此可得，基函数G_i(t)的维格纳-威利分布函数 WVG_i(t,f)：

${\mathrm{WVG}}_i(t,f)={\hat{g}}^2\left(s_i(f-f_i)\right)g^2\left(\frac{1-p_i}{s_i}\right);$

进而，由此得到时频分布函数Ω(t,f)的表达式。

b3）根据雷达目标距离像Θ(t)的时域采样点总数N以及预设定的频域采样点总数M， N≤M≤2N，求得时频分布函数Ω(t,f)对应的维数为M行×N列的时频分布矩阵V_Ω(n,m)：

$V_{\Omega }(n,m)=\underset{i=1}{\overset{I}{\Sigma }}{\left|a_i\right|}^2\left(\begin{array}{c}{\hat{g}}^2\left(s_i(f_1-f_i)\right)\\ {\hat{g}}^2\left(s_i(f_2-f_i)\right)\\ .\\ .\\ .\\ {\hat{g}}^2\left(s_i(f_m-f_i)\right)\\ .\\ .\\ .\\ {\hat{g}}^2\left(s_i(f_M-f_i)\right)\end{array}\right)·\left(\begin{array}{cccccc}{g}^2\left(\frac{t_1-p_i}{s_i}\right)& g^2\left(\frac{t_2-p_i}{s_i}\right)& ...& g^2\left(\frac{t_n-p_i}{s_i}\right)& ...& g^2\left(\frac{t_N-p_i}{s_i}\right)\end{array}\right);$

其中，t_n表示时频分布函数Ω(t,f)在时域上第n个采样点对应的时间，n∈{1,2,…,N}； f_m表示时频分布函数Ω(t,f)在频域上第m个采样点对应的频率，m∈{1,2,…,M}。

在理论上，雷达目标距离像Θ(t)可以视为以时间t作为变量的连续函数，但通过采用 仪器采集而获得雷达目标距离像时，所获得的雷达目标距离像数据实际上是在时域采样的 离散函数，因此在时域采样频率固定的条件下，一定时间内的雷达目标距离像Θ(t)都有其 一定的时域采样点总数N。相应地，频域采样也是如此，因此需要预设定频域采样点总数 M，要求N≤M≤2N，以获得时频分布函数Ω(t,f)对应的维数为M行×N列的时频分布矩阵 V_Ω(n,m)。当然，如果为了减少时频分布矩阵V_Ω(n,m)的数据量，使得运算更加简便，可以 取频域采样点总数的最小值M＝N。

b4）根据预设定的分解维度K和分解门限阈值ε，K＜N，0<ε≤10^-2，对所述时频分布矩 阵V_Ω(n,m)进行非负矩阵分解，获得维数为M行×K列的基元矩阵W_Ω(k,m)和维数为K行×N 列的系数矩阵H_Ω(n,k)，使其满足欧氏距离平方值E(W_Ω,H_Ω)<ε；所述欧氏距离平方值 E(W_Ω,H_Ω)表示时频分布矩阵V_Ω(n,m)与矩阵乘积W_Ω(k,m)·H_Ω(n,k)之间的欧氏距离的平 方，即：

E(W_Ω,H_Ω)=||V_Ω(n,m)-W_Ω(k,m)·H_Ω(n,k)||²；

其中，基元矩阵W_Ω(k,m)和系数矩阵H_Ω(n,k)均为非负矩阵；k∈{1,2,…,K}。

该步骤采用非负矩阵分解的处理方式，以降低时频分布矩阵V_Ω(n,m)的维数。非负矩 阵分解的理想目标是将时频分布矩阵V_Ω(n,m)分解为两个非负矩阵的乘积，即达到：

V_Ω(n,m)=W′(k,m)·H′(n,k)；

其中，矩阵W′(k,m)为M行×K列的非负矩阵，矩阵H′(n,k)为K行×N列的非负矩 阵，由此则可以将维数以N为基数的数学表达转换成维数以K为基数的数学表达，而K＜N， 从而达到降维的目的，以提高本发明方法的运算性能。同时，该矩阵分解过程也是时频相 干结构的分解过程；分解所得的M行×K列的矩阵W′(k,m)可表达为列向量集合，而K行 ×N列的矩阵H′(n,k)可表达为行向量集合，即：

$W^{\prime \prime }(k,m)=[{\overrightarrow{w}}_1^{\prime }\left(m\right),{\overrightarrow{w}}_2^{\prime }\left(m\right),...,{\overrightarrow{w}}_k^{\prime }\left(m\right),...,{\overrightarrow{w}}_K^{\prime }\left(m\right)];$

$H^{\prime }(n,k)==\left(\begin{array}{c}{\overrightarrow{h}}_1^{\prime }\left(n\right)\\ {\overrightarrow{h}}_2^{\prime }\left(n\right)\\ .\\ .\\ .\\ {\overrightarrow{h}}_k^{\prime }\left(n\right)\\ .\\ .\\ .\\ {\overrightarrow{h}}_K^{\prime }\left(n\right)\end{array}\right);$

其中，表示矩阵W′(k,m)中第k列的列向量，表示矩阵H′(n,k)中第k行的行 向量，k∈{1,2,…,K}；矩阵W′(k,m)的每一个列向量则表示了时频分布矩阵V_Ω(n,m) 中的一个时频成分的频率结构特征，而矩阵H′(n,k)中的每一个行向量则表示了时频 分布矩阵V_Ω(n,m)中的一个时频成分的时间结构特征，因此可将矩阵W′(k,m)称为基元矩 阵，将矩阵H′(n,k)称为系数矩阵；而时频成分被划分为K个，分解维度K的取值越大则时 频成分划分越细，因此时频相干结构的表达也就越细致，更有助于时频特征的描述以及雷 达目标的识别，但分解维度K的取值过大则达不到明显的降维目的；因此综合因素考虑， 分解维度K优选的取值范围为0.05N≤K≤0.5N。然而，在实际应用中，很难找到两个非负矩 阵满足非负矩阵分解的理想目标V_Ω(n,m)=W′(k,m)·H′(n,k)，只能够尽可能求取维数为M 行×K列的基元矩阵W_Ω(k,m)和维数为K行×N列的系数矩阵H_Ω(n,k)，基元矩阵W_Ω(k,m) 和系数矩阵H_Ω(n,k)均为非负矩阵，使之满足V_Ω(n,m)≈W_Ω(k,m)·H_Ω(n,k)，对于满足约等 于的判断标准，本发明采用了欧氏距离平方值E(W_Ω,H_Ω)来判断，所述欧氏距离平方值 E(W_Ω,H_Ω)表示时频分布矩阵V_Ω(n,m)与矩阵乘积W_Ω(k,m)·H_Ω(n,k)之间的欧氏距离的平 方，若欧氏距离平方值E(W_Ω,H_Ω)足够小，以至于小于预设定的分解门限阈值ε，0<ε<10^-2， 即满足E(W_Ω,H_Ω)<ε，则认为此时获得的基元矩阵W_Ω(k,m)和系数矩阵H_Ω(n,k)满足非负 矩阵分解的要求；从而可以得到基元矩阵W_Ω(k,m)的列向量集合表达式以及系数矩阵 H_Ω(n,k)的行向量集合表达式：

$W_{\Omega }(k,m)=[{\overrightarrow{w}}_1\left(m\right),{\overrightarrow{w}}_2\left(m\right),...,{\overrightarrow{w}}_k\left(m\right),...,{\overrightarrow{w}}_K\left(m\right)];$

$H_{\Omega }(n,k)==\left(\begin{array}{c}{\overrightarrow{h}}_1^{\prime }\left(n\right)\\ {\overrightarrow{h}}_2^{\prime }\left(n\right)\\ .\\ .\\ .\\ {\overrightarrow{h}}_k^{\prime }\left(n\right)\\ .\\ .\\ .\\ {\overrightarrow{h}}_K^{\prime }\left(n\right)\end{array}\right).$

其中，表示基元矩阵W_Ω(k,m)中第k列的列向量，表示系数矩阵H_Ω(n,k)中 第k行的行向量，k∈{1,2,…,K}。至于基元矩阵W_Ω(k,m)和系数矩阵H_Ω(n,k)的非负矩阵 分解求取方法，通常可随机设定基元矩阵W_Ω(k,m)和系数矩阵H_Ω(n,k)的初始值进行迭代 运算（也可根据数据特点自定义设定初始值，以提高迭代效率），迭代过程可采用Lin提出 的投影梯度边界约束优化算法（可参考文献“Lin C J.Projected gradient methods for  non-negative matrix factorization[J].Neural Computing,2007,19(10):2756-2779.”），每次迭代 均判断是否满足E(W_Ω,H_Ω)<ε，若不满足则继续迭代，直至满足时终止迭代，得到基元矩 阵W_Ω(k,m)和系数矩阵H_Ω(n,k)的分解结果；介于非负矩阵分解为成熟的现有技术，本文 不再对非负矩阵分解的具体运算过程及其迭代过程加以详述。

b5）获取基元矩阵W_Ω(k,m)的各个列向量以及系数矩阵H_Ω(n,k)的各个行向量 k∈{1,2,…,K}；

${\overrightarrow{w}}_k\left(m\right)=\left(\begin{array}{c}{w}_{1,k}\\ w_{2,k}\\ .\\ .\\ .\\ w_{m,k}\\ .\\ .\\ .\\ w_{M,k}\end{array}\right);$

${\overrightarrow{h}}_k\left(n\right)=\left(\begin{array}{cccccc}{h}_{k,1}& h_{k,2}& ...& h_{k,n}& ...& h_{k,N}\end{array}\right);$

其中，w_m，k表示基元矩阵W_Ω(k,m)中第m行第k列的元素；h_k,n表示系数矩阵H_Ω(n,k) 中第k行第n列的元素。

该步骤中，获取基元矩阵W_Ω(k,m)的各个列向量以及系数矩阵H_Ω(n,k)的各个 行向量即相当于获取了时频分布矩阵V_Ω(n,m)中各个时频成分的频率结构特征和时 间结构特征，因此可以将基元矩阵W_Ω(k,m)中的每个列向量看做一个频域特征向量， 将系数矩阵H_Ω(n,k)中的每个行向量看做一个时域特征向量。而频域特征向量和时域 特征向量的向量矩则含有描述联合时频平面矩的信息，频域特征向量和时域特征向量的稀 疏特征则可用于鉴别雷达目标距离像中的瞬态成分和连续成分。因此，接下来求取基元矩 阵W_Ω(k,m)各个列向量以及系数矩阵H_Ω(n,k)各个行向量对应的稀疏特征和 多阶向量矩，并以此作为从雷达目标距离像Θ(t)中提取得到的时频特征。

b6）根据预设定的频域特征阶数Q和时域特征阶数P，计算基元矩阵W_Ω(k,m)各个列 向量对应的频域稀疏特征和各阶频域向量矩以及系数矩阵H_Ω(n,k)各 个行向量对应的时域稀疏特征和各阶时域向量矩

${\mathrm{SF}}_{{\overrightarrow{w}}_k}=\frac{\sqrt{M}-[\left(\underset{m=1}{\overset{M}{\Sigma }}{w}_{m,k}\right)/\sqrt{\underset{m=1}{\overset{M}{\Sigma }}{w_{m,k}}^2}]}{\sqrt{M}-1};$

$M{F}_{{\overrightarrow{w}}_k}\left(q\right)=\underset{m=1}{\overset{M}{\Sigma }}{(m-{\mu }_{{\overrightarrow{w}}_k})}^q{w}_{m,k};$

${\mathrm{ST}}_{{\overrightarrow{h}}_k}=\frac{\sqrt{N}-[\left(\underset{n=1}{\overset{N}{\Sigma }}{h}_{k,n}\right)/\sqrt{\underset{n=1}{\overset{N}{\Sigma }}{w_{k,n}}^2}]}{\sqrt{N}-1};$

$M{T}_{{\overrightarrow{h}}_k}\left(p\right)=\underset{n=1}{\overset{N}{\Sigma }}{(n-{\mu }_{{\overrightarrow{h}}_k})}^p{h}_{k,n}$

其中，k∈{1,2,…,K}；表示基元矩阵W_Ω(k,m)中第k列的列向量对应 的第q阶的频域向量矩，q∈{1,2,…,Q}；表示列向量中各个元素的均值，即 表示系数矩阵H_Ω(n,k)中第k行的行向量对应的第p阶的 时域向量矩，p∈{1,2,…,P}；表示行向量中各个元素的均值，即

该步骤中设置频域特征阶数Q和时域特征阶数P，主要是为了在一定程度上利用多阶 向量矩提高时频特征识别成分，因此频域特征阶数Q和时域特征阶数P的取值都不宜过大， 以避免盲目增多时频特征数据维数；通常情况下，频域特征阶数Q的优选取值范围为 1≤Q≤6，时域特征阶数P的优选取值范围为1≤P≤6，二者可以各自独立取值。

b7）将基元矩阵W_Ω(k,m)中第k列的列向量对应的频域稀疏特征和各阶频 域向量矩以及系数矩阵H_Ω(n,k)中第k行的行向量对应的时域稀疏特征和各阶时域向量矩的集合作为雷达目标距离像Θ(t)的第k个时频特征矢量ξ_k，即：

${\xi }_k=[{\mathrm{MT}}_{{\overrightarrow{h}}_k}\left(1\right),{\mathrm{MT}}_{{\overrightarrow{h}}_k}\left(2\right),...,{\mathrm{MT}}_{{\overrightarrow{h}}_k}\left(p\right),...{,\mathrm{MT}}_{{\overrightarrow{h}}_k}\left(p\right),$

${\mathrm{MF}}_{{\overrightarrow{w}}_k}\left(1\right),{\mathrm{MF}}_{{\overrightarrow{w}}_k}\left(2\right),...,{\mathrm{MF}}_{{\overrightarrow{w}}_k}\left(q\right),...,{\mathrm{MF}}_{{\overrightarrow{w}}_k}\left(Q\right),;$

${\mathrm{ST}}_{{\overrightarrow{h}}_k}{\mathrm{SF}}_{{\overrightarrow{w}}_k}]$

由此，得到雷达目标距离像Θ(t)的时频特征矢量序列ζ={ξ₁,ξ₂,…,ξ_k,…,ξ_K}，并将所 述时频特征矢量序列ζ作为从雷达目标距离像Θ(t)提取的时频特征。

至此，步骤B）的处理过程即得以完成，由此确定了雷达目标距离像的时频特征提取 方案，也即相当于建立了雷达目标距离像对应的时频特征模型。

C）将根据步骤B）所述方案提取的雷达目标距离像的时频特征作为雷达目标识别特 征，进行雷达目标识别。具体而言，该步骤通常的处理流程为，采集已知雷达目标的雷达 目标距离像，采用步骤B）所述方案提取已知雷达目标的雷达目标距离像的时频特征并以 此作为已知雷达目标对应的时频特征模板；然后对于待测雷达目标的雷达目标距离像，也 采用采用步骤B）所述方案提取待测雷达目标的雷达目标距离像的待识别时频特征，而后 将该待识别时频特征与各个时频特征模板进行匹配，若该待识别时频特征存在相匹配的时 频特征模板，则判定待测雷达目标与所述相匹配的时频特征模板对应的已知雷达目标属于 同类雷达目标，即完成雷达目标识别。在计算机匹配识别过程可采用的现有成熟算法则非 常多，本文在此不再一一列举。

由此，本发明提供了一种基于雷达目标距离像时频特征提取的雷达目标识别方法，该 方法能够通过计算机运行执行，减少雷达目标识别观测人员的工作量，并且有助于提高雷 达目标识别效率。

下面通过实施例对本发明的技术方案作进一步的说明。

实施例：

本实施例采用本发明的雷达目标识别方法，对MSTAR公司公开发布的五类雷达目标 进行识别。这五类雷达目标均为地面军事车辆，其雷达目标代号分别为BMP2_SN9563、 BTR70_SNC71、T72_SN132、D7、ZSU23_4（后文中为叙述方便，分别将这雷达目标代号 简称为BMP2、BTR70、T72、D7、ZSU23），五类雷达目标的光学图像分别如图2中的（2a） ～（2e）所示，其中BMP2为步兵战车，BTR70为装甲运兵车，T72为坦克，D7为推土机， ZSU23为防空火炮。雷达目标信号是一种包含了幅度和相位的复图像数据，经过一系列的 滤波变换运算处理能够转化为雷达目标距离像序列数据，每个雷达目标距离像序列数据中 包含有按方位角大小的顺序均匀排列的若干个雷达目标距离像，所包含雷达目标距离像的 数量根据滤波变换运算处理过程确定，每个雷达目标距离像表征了雷达目标在一个方位角 位置上的距离特征。例如，一个BMP2雷达目标信号图像处理为包含100个方位角的雷达 目标距离像序列数据，该雷达目标距离像序列数据中第1个、第50个和第100个雷达目 标距离像分别如图3～图5所示。对这些雷达目标距离像进行特征提取，即可根据提取的特 征对雷达目标进行识别。下面通过本实施例，对本发明雷达目标识别方法的具体应用过程 以及技术效果进行进一步的说明。

A）建立包含有X个基函数的基函数库G_D。

本实施例中构建了包含66612个基函数的基函数库G_D（即X=66612），使其构成一个 完备和冗余的基函数库。

B）提取雷达目标距离像的时频特征。

采用本发明方法提取雷达目标距离像的时频特征。以图3所示的雷达目标距离像为例， 提取其时频特征的过程如下：按步骤b1），用基函数库G_D中的基函数对图3所示的雷达目 标距离像进行匹配追踪分解，分解级数为100级（即取I＝100），得到其分解表达式；然后 按步骤b2），取该雷达目标距离像分解表达式的匹配追踪分解组成部分作为该雷达目标距 离像对应的时频函数，并求取该时频函数的时频分布函数Ω(t,f)；由于图3所示的雷达目 标距离像的采样点数为100个（即N＝100），本实施例中预设定频域采样点总数M＝100，因 此按步骤b3），求得上述时频分布函数Ω(t,f)对应的维数为100行×100列的时频分布矩阵 V_Ω(n,m)，其时频分布矩阵V_Ω(n,m)的图形如图6所示；本实施例中预设定分解维度K＝10 和分解门限阈值ε=10^-2，因此再按步骤b4），对所述时频分布矩阵V_Ω(n,m)进行非负矩阵分 解，获得维数为100行×10列的基元矩阵W_Ω(k,m)和维数为10行×100列的系数矩阵 H_Ω(n,k)，并按步骤b5）获取基元矩阵W_Ω(k,m)的10个列向量以及系数矩阵H_Ω(n,k) 的10个行向量其中基元矩阵W_Ω(k,m)的10个列向量如图7所示，系数矩阵H_Ω(n,k) 的10个行向量图如图8所示；本实施例中预设定频域特征阶数Q=6和时域特征阶数P=6， 因此按步骤b6）计算基元矩阵W_Ω(k,m)各个列向量对应的频域稀疏特征和6阶 频域向量矩以及系数矩阵H_Ω(n,k)各个行向量对应的时域稀疏特征和6 阶时域向量矩从而按步骤b7）得到图3所示的雷达目标距离像的时频特征矢量 序列，并将之作为从图3所示的雷达目标距离像提取的时频特征。

根据本发明中步骤B）提取雷达目标距离像的时频特征的方法，分别提取本实施例中 BMP2、BTR70、T72、D7、ZSU23五类雷达目标的雷达目标距离像时频特征，最终获得 该五类雷达目标的雷达目标距离像时频特征分布情况如图9所示；从图9中可以看到，根 据本发明中步骤B）提取的时频特征，其同类目标的特征聚集性以及不同类目标的特征分 离性都非常好，能够很好的用于雷达目标的识别。

C）将根据步骤B）所述方案提取的雷达目标距离像的时频特征作为雷达目标识别特 征，进行雷达目标识别。

本实施例在识别过程中，针对BMP2、BTR70、T72、D7、ZSU23五类雷达目标，分 别采集了233个BMP2训练样本、233个BTR70训练样本、232个T72训练样本、233个 D7训练样本、234ZSU23个训练样本作为五类已知雷达目标并运算获得其雷达目标距离 像，并采用本发明步骤B）的方案提取了这些已知雷达目标距离像的时频特征并以此作为 已知雷达目标对应的时频特征模板；然后，对于上述五类雷达目标的每一类采集了196个 测试样本作为待测雷达目标并运算获得其雷达目标距离像，采用本发明步骤B）的方案提 取待测雷达目标的雷达目标距离像的待识别时频特征；而后，本实施例采用隐马尔科夫模 型分类识别算法，分别建立各个时频特征模板以及待识别时频特征的隐马尔科夫模型，然 后对待识别时频特征与各个时频特征模板进行匹配分类识别，每个隐马尔科夫模型中隐状 态数为120，高斯混合数为6；如果将BMP2、BTR70、T72、D7、ZSU23五类雷达目标分 别称之为第1类、第2类、第3类、第4类、第5类雷达目标，令d_i，j表示将第i类测试样 本识别为第j类雷达目标的个数，i∈{1,2,3,4,5}，j∈{1,2,3,4,5}；令P_ci表示第i类雷达目 标的正确识别率，P_ci=d_ii/(d_i1+d_i2+…+d_i5)×100%；另表示平均识别率， 由此，得到各类雷达目标测试样本的识别结果如表1 所示：

表1

从表1所示的识别结果可以看到，采用本发明方法结合隐马尔科夫模型分类识别算法 对雷达目标加以识别，正确识别率最佳达到了98.47%，平均识别率可以达到97.55%，识 别效果非常好。

为了加以对比，本实施例还采用了另一种雷达目标识别方法对上述五种雷达目标各自 的196个测试样本进行识别，在该作为对比的雷达目标识别方法中，采用主成分分析算法 （主成分分析算法为现有技术中的陈述算法，可参考文献“Behnaz G.,Sridhar K. Quantifcation and localization of features in time-frequency plane[C].Proceedings of IEEE  Canadian Conference on Electrical and Computer Engineering,May 4-7,2008.”）提取雷达目标 距离像的时频特征，然后再同样采用隐马尔科夫模型分类识别算法根据提取的时频特征进 行雷达目标识别，每个隐马尔科夫模型中隐状态数同样为120，高斯混合数同样为6；该 作为对比的雷达目标识别方法，其时频特征提取方案与本发明不同，而分类识别算法与本 实施例中上述采用本发明方案进行雷达目标识别的实施例相同，由此比较不同频特征提取 方案对雷达目标识别结果的影响。同样分别采用上述233个BMP2训练样本、233个BTR70 训练样本、232个T72训练样本、233个D7训练样本、234ZSU23个训练样本作为五类已 知雷达目标并运算获得其雷达目标距离像，对五种雷达目标各自的196个测试样本进行识 别，该作为对比的雷达目标识别方法的识别结果如表2所示：

表2

通过表2与表1进行对比可以看到，在分类识别算法相同的条件下，该作为对比的雷 达目标识别方法的正确识别率最佳为86.22%，与本发明方法相比低了12.25%；而平均识 别率为83.77%，与本发明方法相比低了13.78%。这是因为时频平面的主成分分析特征仅 利用了整体信息，而基于本发明方法提取得到的时频特征则充分挖掘了信号的时间局部结 构和频谱局部结构特征信息，因此可以进一步提高雷达目标识别的准确性。

最后说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实 施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方 案进行修改或者等同替换，而不脱离本发明技术方案的宗旨和范围，其均应涵盖在本发明 的权利要求范围当中。