CN104021400A - 基于分数谱稀疏表示的sar图像目标识别方法 - Google Patents

Abstract

本发明公开了一种基于分数谱稀疏表示的SAR图像目标识别方法，本发明方法包括：用训练图像集构建原子库和对待测图像的分类识别两个过程，具体的在特征提取方面同时提取了分数谱域和图像像素域两类特征；在特征分类方面，利用图像分数谱和像素两类特征的相关性，设计新的联合稀疏表示和稀疏分类方法对待测图像进行分类识别，通过利用分数域时频变换获得了信号的高分辨频谱特征，将其与图像空间像素域特征进行联合稀疏表示和稀疏分类，有效地改善了现有方法在特征提取和特征分类中的问题，提高了SAR目标识别的准确率。

Description

基于分数谱稀疏表示的SAR图像目标识别方法

技术领域

本发明属于图像处理技术领域中，涉及图像分类识别方法，具体涉及运用分数谱稀疏表示来进行合成孔径雷达(Synthetic Aperture Radar，SAR)图像目标识别的方法。

背景技术

合成孔径雷达是一种工作在微波波段的相干成像雷达，以其高分辨和全天候、全天时、大面积的成像探测能力，成为世界各国普遍重视的对地观测手段，具有良好的应用前景。SAR目标识别是SAR图像解译的一个重要方面，目的是实现对图像中感兴趣目标的自动分类识别。SAR目标识别过程通常包括对特征提取和特征分类两个环节。

在特征提取方面，已经提出了如Gabor变换，小波变换，主成分分析，几何投影变换等多种方法；SAR目标具有多变性特点：除了成像系统固有的乘性噪声斑点外，SAR目标受到如目标方位角、姿态变化、背景地表变化、有意或无意的遮挡、杂波干扰等多种因素的影响，这使得其识别任务变得非常困难。因此，基于单一特征的识别方式难以达到高精度、高可靠性的要求，需要考虑一种多特征融合的识别技术和手段。目前，在图像信号分析中，分数域时频分析(Fractional time-frequency analysis，FrTFA)方法逐渐得到关注，该方法利用时频旋转特性及能量聚焦性，获取信号的高分辨率频谱。相对于传统时频分析方法，分数域分析改善了信号的时频分布，更好地反映信号的局部特性。目前分数域时频分析已经被应用于科学研究和工程技术领域，如信号检测与参数估计，其中包括SAR图像目标检测。但至今未见有将分数域时频分析用于SAR目标自动识别的报道。

在特征分类方面，传统采用模式识别、机器学习的方法来训练分类模型。典型的方法有：线性判别分析(Linear Discriminant Analysis，LDA)、支持向量机(Support Vector Machine，SVM)、神经网络等，其中，SVM的使用最为广泛。这些分类方法一般针对低维空间，在训练中存在过拟合问题，且计算复杂度较高。近年来，一种基于稀疏表示的分类器(SparseRepresentation Classifier，SRC)也应用到SAR目标识别中，获得了很好的分类效果，且复杂度较低，从而弥补了传统统计分析理论和机器学习方法的不足。但SRC的存在如下问题：SRC过程是先在训练图像集(或原子库)中选择出一系列原子，用其对测试图像进行稀疏线性表示，然后将表示误差最小的类作为测试图像的最终类别。在这个过程中，原子选择非常关键，但目前的方法趋向于选择一组相关的原子来表示一个测试样本，并不能保证所选原子来源于同一个正确的类，造成相似的局部特征可能会被不同类别的原子组表示，影响了分类结果。

发明内容

本发明的目的是针对现有技术存在的上述缺陷，提出一种基于分数谱稀疏表示的SAR图像目标识别方法。

本发明的技术方案是：一种基于分数谱稀疏表示的SAR图像目标识别方法，具体包括如下步骤：

S1.选择目标区域：在带有类别标签的训练图像中选择目标所在的区域，目标区域内的图像数据作为目标图像；

S2.提取分数谱特征：对目标图像进行分数域时频变换，提取其分数谱特征；

S3.构建原子库：

S31.将每幅训练的目标图像的数据排列为一列，对其进行降维处理，并将得到的低维列向量作为一原子，加入到图像特征原子库，设为D_p,D_p的一列为某训练图像的图像特征数据；

S32.将每幅训练的目标图像经步骤S2处理后的分数谱特征数据排列为一列，对其进行降维处理，并将得到的低维列向量作为一原子，加入到分数谱特征原子库，设为D_f，D_f的一列为某训练图像的分数谱特征数据；

S33.将图像特征原子库与分数谱特征原子库合并为最终的原子库，设为D；

S4.联合稀疏表示：根据步骤S33获得的原子库D，对输入的某测试图像进行联合稀疏表示，其具体过程如下：

S41.将输入的某测试图像，按步骤S1进行目标定位与降维，并将降维后的图像数据矢量化为一列作为待测图像特征数据，设为t_p；

按步骤S2提取分数谱特征，将分数谱特征数据矢量化为一列作为待测试图像的分数谱特征数据，设为t_f；

将两个特征数据联合形成最终的测试数据，设为T，即T＝[t_p,t_f]；

S42.使用原子库D，对测试数据T进行稀疏分解；

S5.联合稀疏识别：根据原子库原子的类别标签，对每一类别，只用稀疏分解得到的该类的原子和对应系数分别对t_p和t_f进行重建，计算两个重建误差的累计；选择累计误差最小的类别为输入测试图像的类别。

进一步的，步骤S2所述的提取其分数谱特征的具体过程如下：对图像的每一行和每一列分别进行分数谱特征的提取，设目标图像的某行(或列)信号为g(t)，则特征提取按如下步骤进行：

S21.寻找分数域时频变换的最优旋转阶次r^*：

$r^{*}=\underset{r}{\arg \min }\{{({W_t}^2\left(g\right)+{W_t}^2\left(h\right))}^{1/2}\·{(W_f^2\left(g\right)+W_f^2\left(h\right))}^{1/2}\}---\left(1\right)$

其中，h为窗函数，这里等号右边的h表示分数阶傅里叶变换时的旋转阶次，在每个阶次时求{}括号内的值，然后比较哪个值最小，当{}括号内的值最小时，即得到优旋转阶次r^*；

W_t(g)为g(t)的时间宽度，其形式为：

$W_t\left(g\right)=\frac{{\left[{\∫(t-{\mathrm{\ζ}}_t)}^2{\left|g\left(t\right)\right|}^2\mathrm{dt}\right]}^{1/2}}{\left|\right|g\left|\right|}---\left(2\right)$

其中，||g||为g(t)的2范数，即g(t)中各个元素平方之和再开根号，ζ_t为

${\mathrm{\ζ}}_t=\frac{\∫t{\left|g\left(t\right)\right|}^2\mathrm{dt}}{{\left|\right|g\left|\right|}^2}---\left(3\right)$

W_f,r(g)为g(t)的频率宽度，其形式为：

$W_{f,r}\left(g\right)=\frac{{[\∫{(f-{\mathrm{\ζ}}_f)}^2{\left|{\tilde{g}}_r\left(f\right)\right|}^2\mathrm{df}]}^{1/2}}{\left|\right|g\left|\right|}---\left(4\right)$

其中，为g(t)的傅里叶变换，ζ_f为

${\mathrm{\ζ}}_f=\frac{\∫f{\left|{\tilde{g}}_r\left(f\right)\right|}^2\mathrm{df}}{{\left|\right|g\left|\right|}^2}---\left(5\right)$

S22.计算该最优旋转阶次r^*下的最优窗函数h^*：

$h^{*}=\mathrm{FT}\left(e^{{-\mathrm{\πt}}^2{W}_{{f,r}^{*}}\left(g\right)/W_t\left(g\right)}\right)---\left(6\right)$

其中，FT(x)表示x的傅里叶变换；

S23.根据谱峰度系数最大准则，求解如下最优化问题，得到最优变换阶次q^*：

$q^{*}=\arg \underset{q}{\max }(\frac{E\left[{{\tilde{g}}_q}^2\left(f\right){\left(\stackrel{\‾}{{\tilde{g}}_q\left(f\right)}\right)}^2\right]}{E{\left[{\tilde{g}}_q\left(f\right)\stackrel{\‾}{{\tilde{g}}_q\left(f\right)}\right]}^2}-2)---\left(7\right)$

其中，为g(t)的q阶傅里叶变换，表示x的复共轭，E[x]表示求x的数学期望。

S24.在最优窗函数h^*和最优变换阶次q^*下，对g(t)进行短时傅里叶变换，得到分数域时频分布FT(t,u)：

$\mathrm{FT}(t,u)={\∫}_{-\∞}^{\∞}g\left(\mathrm{\τ}\right)h^{*}(t-\mathrm{\τ})\frac{\sqrt{1-j\cot (q^{*}\mathrm{\π}/2)}}{2\mathrm{\π}}\exp (\frac{j({\mathrm{\τ}}^2+u^2)}{2\tan (q^{*}\mathrm{\π}/2)}-\frac{\mathrm{j\τu}}{\sin (q^{*}\mathrm{\π}/2)})\mathrm{d\τ}---\left(8\right)$

其中，u为频率，τ为积分变量；

S25.对分数域时频分布FT(t,u)进行局部区域能量聚焦分析，提取能量衰减斜率特征作为图像的分数谱特征F。

进一步的，步骤S42所述的稀疏分解的具体过程如下：

设分解的系数x＝[x_p,x_f]，其中，x_p对应t_p在D_p上的分解系数，x_f对应t_f在D_f上的分解系数，则求解如下最优化问题得到x：

其中，L为设定的稀疏度，Class(z)表示系数类别判定算子，其具体描述为：

·取系数z的非零元素组成矢量设z有L个非零元素，则有z_i≠0,i＝1,2,3,...,L；

·求的元素所对应的原子的类别标签，设所对应原子的类别标签为c_i，则有Class(z)＝[c₁,c₂,c₃,…,c_L]。

这里可以看出，在原子选择上要满足以下3个条件：

(1)仅用D_p中的原子来分解t_p部分；

(2)仅用D_f中的原子来分解t_f部分；

(3)每次原子选择时，在D_p和D_f中选择的原子都来源于同类训练图像。

本发明的有益效果：本发明方法利用分数域时频变换获得了信号的高分辨频谱特征，将其与图像空间像素域特征进行联合稀疏表示和稀疏分类，有效地改善了现有方法在特征提取和特征分类中的问题，提高了SAR目标识别的准确率。

附图说明

图1为本发明方法流程示意图；

图2为MSTAR图像数据库中的T72，BTR70与BMP23种目标及其对应的SAR图像、分数谱特征图像示例图。

具体实施方式

下面结合附图和具体的实施例对本发明作进一步详细的描述：

本发明实施例提供的方法包括：用训练图像集构建原子库和对待测图像的分类识别两个过程，具体的在特征提取方面同时提取了分数谱域和图像像素域两类特征；在特征分类方面，利用图像分数谱和像素两类特征的相关性，设计新的联合稀疏表示和稀疏分类方法对待测图像进行分类识别，从而弥补现有的方法在特征提取和特征分类上的不足，提高了目标识别的准确率。

本发明实施例测试的数据具体为：对MSTAR图像数据库中的T72，BTR70和BMP2等3类目标进行分类识别，其中，T72定为第1类，包含sn-132、sn-812和sn-s7三种型号，BTR70定为第2类，包含sn-71一种型号，BMP2定为第3类，包含sn-c9563、sn-c9566和sn-c21三种型号。每个目标图像成像的俯视角有15度和17度两个角度，方位角在360度内连续变化。

将俯视角为17度的图像作为训练图像，训练图像共698张，俯视角为15度的图像作为测试图像，测试图像共1365张，详细类型与数量如表1所示。

表 1

S1.选择目标区域：在每一训练图像的中心位置，选取大小为64×64的区域为目标区域，设代表该目标区域内的图像(目标图像)；

S2.提取分数谱特征F：对图像的每一行和每一列分别进行分数谱特征的提取，设目标图像的某行(或列)信号为g(t)，则特征提取按如下步骤进行：

S21.根据g(t)寻找分数域时频变换的最优旋转阶次r^*；

S22.根据最优旋转阶次r^*计算该阶次下的最优窗函数h^*；

S23.根据峰度系数最大准则计算最优变换阶次q^*；

S24.在最优窗函数h^*和最优变换阶次q^*下，对g(t)进行短时傅里叶变换，得到分数域时频分布FT(t,u)；

S25.对分数域时频分布FT(t,u)进行局部区域能量聚焦分析，提取能量衰减斜率特征作为图像的分数谱特征F。

图2为MSTAR图像数据库中的3种目标及其对应的SAR图像和分数谱特征图像的示例图。其中，2(a)为原目标图像，2(b)为SAR图像，2(c)为分数谱特征图像。图中的3列分别对应T72，BTR70和BMP2目标。可以看到，同一目标的SAR图像与对应的分数谱特征图像具有很强的相似性，这是后续将二者进行联合稀疏表示的依据。

S3.构建原子库：

S31.将目标图像Y的数据逐列扫描，排列为一列向量对y进行随机采样降维生成图像特征向量如下：其中，为随机采样矩阵，该矩阵的各行向量元素是由一组归一化的零均值独立正态分布随机变量组成，将作为一原子，加入到图像特征原子库D_p中。

S32.将分数谱特征F的数据逐列扫描，排列为一列向量对f采用与步骤S31相同的随机降维方式，生成分数谱特征向量将作为一原子，加入到分数谱特征原子库D_f中。

S33.将图像特征原子库与分数谱特征原子库合并为最终的原子库原子库共有1396个原子，包含第1类的图像和分数谱特征原子共232×2＝464个；第2类的图像和分数谱特征原子共233×2＝466个；第3类的图像和分数谱特征原子共233×2＝466个。

S4.联合稀疏表示：在原子库D的基础上，对输入的测试图像进行联合稀疏表示，以第1类T72中型号为sn-132的SAR图像HB03335作为测试例子，其具体过程如下：

S41.将图像HB03335按步骤S1进行目标定位与降维，并将降维后的图像数据矢量化为1列作为图像特征数据t_p，其值为[11,12,11,13,16,12,...,11,11,12,17,14,11]^T。

按步骤S2提取分数谱特征，将分数谱特征数据矢量化为1列作为分数谱特征数据t_f,其值为[2,2,2,3,2,2,...,1,1,2,3,2,1]^T。

将两个特征数据联合形成最终的测试数据T如下：

$T={\left[\begin{array}{ccccccccccccc}11,& 12,& 11,& 13,& 16,& 12,& ...,& 11,& 11,& 12,& 17,& 14,& 11\\ 2,& 2,& 2,& 3,& 2,& 2,& ...,& 1,& 1,& 2,& 3,& 2,& 1,\end{array}\right]}^T$

S42.基于原子库D，对测试数据T进行联合稀疏分解。该分解过程通过迭代来完成，

一次迭代选择一对类别相同的原子。

在单次迭代中，首先选择某类别中对应系数平方和最大的一对原子，设为(A_P,A_f)，A_P和A_f分别来自D_p和D_f，然后计算当前残差信号E在所有已选择的原子构成的空间上的投影P，最后更新残差：E＝E-P。设置迭代次数为15，则对T进行稀疏分解得到的15个原子对分别为：(140,200)、(30,13)、(584,527)、(119,30)、(425,281)、(13,107)、(447,276)、(668,608)、(135,31)、(261,358)、(206,53)、(126,35)、(564,634)、(579,584)、(444,316)，其中，每个()中的两数分别表示原子A_P和A_f的编号。这些原子对对应的类别编号分别为：(1,1)、(1,1)、(3,3)、(1,1)、(2,2)、(1,1)、(2,2)、(3,3)、(1,1)、(2,2)、(1,1)、(1,1)、(3,3)、(3,3)、(2,2)。可见，在选择的15个原子对中，属于第1类T72的有7个，属于第2类BTR70的有4个，属于第3类BMP2的有4个。这些原子对对应的稀疏系数分别为(273.71,234.47)、(186.11,31.49)、(62.48,54.77)、(-52.25,57.39)、(-20.30,-16.39)、(47.21,17.44)、(108.91,42.43)、(-30.07,-9.94)、(30.60,80.63)、(-29.42,-39.98)、(-72.56,-21.71)、(73.97,54.00)、(-48.35,-19.73)、(39.64,30.31)、(-33.42,11.41)。

S5.联合稀疏识别：根据原子库原子的类别标签，对每一类别，只用前分解得到的该类的原子和其对应系数对t_p和t_f进行重建，计算这两个重建误差求平方和作为该类的总误差，选择总误差最小的类别作为输入图像的类别。即，求解如下最优化问题，得到输入图像的类别c：

$c=\arg \underset{c}{\min }({\left|\right|t_p-D_p{\mathrm{\δ}}^c\left(x_p\right)\left|\right|}_2^2+{\left|\right|t_f-D_f{\mathrm{\δ}}^c\left(x_f\right)\left|\right|}_2^2)---\left(10\right)$

其中，δ^c(z)表示保留系数z中第c类对应的稀疏系数，而将其他类对应的稀疏系数置0。

实施中得到三个类别的总误差分别为25173，361475，347438，第1个类的总误差最小，故判断本测试图像HB03335属于第1类。

本实施方法按照图1的流程对测试集的图像逐一进行识别。表2显示了最终的分类结果。该结果以识别率来表示，其中识别率定义为正确识别的目标数目与总目标数目之比。表中同时列出了经典稀疏分类方法SRC以及2010年《电子学报》上的基于非负矩阵分解的SAR图像目标识别方法的识别结果，以便比较。可以看到，相对于其他两种方法，本实施方法的平均识别率最高，分别比SRC高2.49％，比电子学报方法高1.32％。

表 2

测试样本类别	电子学报方法	SRC	本实施方法
				T72	93.13	95.02	95.70
BTR70	96.43	100.00	100.00

BMP2	91.48	85.69	90.80
				平均识别率	92.89	91.72	94.21

Claims (3)

1.一种基于分数谱稀疏表示的SAR图像目标识别方法，具体包括如下步骤：

S1.选择目标区域：在带有类别标签的训练图像中选择目标所在的区域，目标区域内的图像数据作为目标图像；

S2.提取分数谱特征：对目标图像进行分数域时频变换，提取其分数谱特征；

S3.构建原子库，具体过程如下：

S31.将每幅训练的目标图像的数据排列为一列，对其进行降维处理，并将得到的低维列向量作为一原子，加入到图像特征原子库，设为D_p,D_p的一列为某训练图像的图像特征数据；

S32.将每幅训练的目标图像经步骤S2处理后的分数谱特征数据排列为一列，对其进行降维处理，并将得到的低维列向量作为一原子，加入到分数谱特征原子库，设为D_f，D_f的一列为某训练图像的分数谱特征数据；

S33.将图像特征原子库与分数谱特征原子库合并为最终的原子库，设为D；

S4.联合稀疏表示：根据步骤S33获得的原子库D，对输入的某测试图像进行联合稀疏表示，具体过程如下：

S41.将输入的某测试图像，按步骤S1进行目标定位与降维，并将降维后的图像数据矢量化为一列作为待测图像特征数据，设为t_p；

按步骤S2提取分数谱特征，将分数谱特征数据矢量化为一列作为待测试图像的分数谱特征数据，设为t_f；

将两个特征数据联合形成最终的测试数据，设为T，即T＝[t_p,t_f]；

S42.使用原子库D，对测试数据T进行稀疏分解；

S5.联合稀疏识别：根据原子库原子的类别标签，对每一类别，只用稀疏分解得到的该类的原子和对应系数分别对t_p和t_f进行重建，计算两个重建误差的累计；选择累计误差最小的类别为输入测试图像的类别。

2.根据权利要求1所述的基于分数谱稀疏表示的SAR图像目标识别方法，其特征在于，步骤S2所述的提取其分数谱特征的具体过程如下：对图像的每一行和每一列分别进行分数谱特征的提取，设目标图像的某行(或列)信号为g(t)，则特征提取按如下步骤进行：

S21.寻找分数域时频变换的最优旋转阶次r^*：

$r^{*}=\underset{r}{\arg \min }\{{({W_t}^2\left(g\right)+{W_t}^2\left(h\right))}^{1/2}\·{(W_f^2\left(g\right)+W_f^2\left(h\right))}^{1/2}\}---\left(1\right)$

其中，h为窗函数，这里等号右边的h表示分数阶傅里叶变换时的旋转阶次，在每个阶次时求{}括号内的值，然后比较哪个值最小，当{}括号内的值最小时，即得到优旋转阶次r^*；

W_t(g)为g(t)的时间宽度，其形式为：

$W_t\left(g\right)=\frac{{\left[{\∫(t-{\mathrm{\ζ}}_t)}^2{\left|g\left(t\right)\right|}^2\mathrm{dt}\right]}^{1/2}}{\left|\right|g\left|\right|}---\left(2\right)$

其中，||g||为g(t)的2范数，即g(t)中各个元素平方之和再开根号，ζ_t为

${\mathrm{\ζ}}_t=\frac{\∫t{\left|g\left(t\right)\right|}^2\mathrm{dt}}{{\left|\right|g\left|\right|}^2}---\left(3\right)$

W_f，r(g)为g(t)的频率宽度，其形式为：

$W_{f,r}\left(g\right)=\frac{{[\∫{(f-{\mathrm{\ζ}}_f)}^2{\left|{\tilde{g}}_r\left(f\right)\right|}^2\mathrm{df}]}^{1/2}}{\left|\right|g\left|\right|}---\left(4\right)$

其中，为g(t)的傅里叶变换，ζ_f为

${\mathrm{\ζ}}_f=\frac{\∫f{\left|{\tilde{g}}_r\left(f\right)\right|}^2\mathrm{df}}{{\left|\right|g\left|\right|}^2}---\left(5\right)$

S22.计算该最优旋转阶次r^*下的最优窗函数h^*：

$h^{*}=\mathrm{FT}\left(e^{{-\mathrm{\πt}}^2{W}_{{f,r}^{*}}\left(g\right)/W_t\left(g\right)}\right)---\left(6\right)$

其中，FT(x)表示x的傅里叶变换；

S23.根据谱峰度系数最大准则，求解如下最优化问题，得到最优变换阶次q^*：

$q^{*}=\arg \underset{q}{\max }(\frac{E\left[{{\tilde{g}}_q}^2\left(f\right){\left(\stackrel{\‾}{{\tilde{g}}_q\left(f\right)}\right)}^2\right]}{E{\left[{\tilde{g}}_q\left(f\right)\stackrel{\‾}{{\tilde{g}}_q\left(f\right)}\right]}^2}-2)---\left(7\right)$

其中，为g(t)的q阶傅里叶变换，表示x的复共轭，E[x]表示求x的数学期望。

S24.在最优窗函数h^*和最优变换阶次q^*下，对g(t)进行短时傅里叶变换，得到分数域时频分布FT(t,u)：

$\mathrm{FT}(t,u)={\∫}_{-\∞}^{\∞}g\left(\mathrm{\τ}\right)h^{*}(t-\mathrm{\τ})\frac{\sqrt{1-j\cot (q^{*}\mathrm{\π}/2)}}{2\mathrm{\π}}\exp (\frac{j({\mathrm{\τ}}^2+u^2)}{2\tan (q^{*}\mathrm{\π}/2)}-\frac{\mathrm{j\τu}}{\sin (q^{*}\mathrm{\π}/2)})\mathrm{d\τ}---\left(8\right)$

其中，u为频率，τ为积分变量；

S25.对分数域时频分布FT(t,u)进行局部区域能量聚焦分析，提取能量衰减斜率特征作为图像的分数谱特征F。

3.根据权利要求1或2所述的基于分数谱稀疏表示的SAR图像目标识别方法，其特征在于，步骤S42所述的稀疏分解的具体过程如下：

设分解的系数x＝[x_p,x_f]，其中，x_p对应t_p在D_p上的分解系数，x_f对应t_f在D_f上的分解系数，则求解如下最优化问题得到x：

其中，L为设定的稀疏度，Class(z)表示系数类别判定算子，其具体描述为：

·取系数z的非零元素组成矢量设z有L个非零元素，则有z_i≠0,i＝1,2,3,...,L；

·求的元素所对应的原子的类别标签，设所对应原子的类别标签为c_i，则有Class(z)＝[c₁,c₂,c₃,…,c_L]。