CN101526995A - 基于对角子类判决分析的合成孔径雷达目标识别方法 - Google Patents

Abstract

本发明提出一种对角子类判决分析的合成孔径雷达目标识别方法，主要解决现有的合成孔径雷达目标识别性能差的问题。该方法的过程是：对原始图像进行自适应的阈值分割、形态学滤波、几何聚类操作和图像增强的预处理；对预处理后的每类目标采用二维快速全局k均值聚类算法进行最优的子类划分；用对角子类判决分析或者对角子类判决分析和二维子类判决分析找到最优的投影矩阵；将预处理后的训练和测试图像向投影矩阵投影，得到它们的特征矩阵；计算测试目标与每个训练目标的特征矩阵之间的欧氏距离，并采用最近邻准则确定测试目标的类别属性。仿真实验表明，本发明抑制背景杂波效果好、目标图像质量高、特征维数低的优点，可用于遥感系统中。

Description

基于对角子类判决分析的合成孔径雷达目标识别方法

技术领域

本发明属于雷达目标识别技术领域，具体的说是一种合成孔径雷达(SAR)的目标识别方法，用于遥感系统中。

背景技术

由于合成孔径雷达SAR成像不受天气、光照等条件的限制，可对感兴趣的目标进行全天候、全天时的侦察，因而广泛应用于军事和遥感领域。由于合成孔径雷达SAR技术对地面目标尤其是地面静止目标探测的优势，在现代战场感知中具有很好的应用前景，所以如何对合成孔径雷达SAR图像进行快速准确的理解与识别已经越来越引起重视。

在一个给定的场景中，通常既包含感兴趣的目标，又包含大量的背景杂波，如果直接对原始目标图像进行特征提取及识别，背景杂波势必会影响识别性能。因此需要对原始图像进行预处理，将目标从杂波背景中分割出来，以减弱背景杂波对识别性能的影响。

此外，若用图像的全部像素点作为识别的特征，则维数过大，对内存、运算量要求较高，还会出现信息“冗余”，造成识别性能损失。因此，需要提取有效的目标特征，以降低样本维数、提高识别性能。

目前国内外有许多机构都开展了合成孔径雷达SAR图像目标识别技术的研究，这些研究在图像预处理、特征提取和分类器设计等方面各不相同。文献[1，2]是直接在原始图像域内用支持向量机SVM完成对测试目标的识别。文献[3]首先对原始图像作对数变换，再作二维傅里叶变换后取一半的幅频信息，采用主分量分析PCA、核主分量分析KPCA提取特征，最后用支撑向量机SVM对测试目标进行识别。上述这些方法的共同缺点就是没有将目标从杂波背景中提取出来，由于地面目标所处的背景杂波具有多样性，不同的背景杂波特性会影响识别性能。另外，文献[1，2]直接采用原始图像的像素点作为目标的特征，这样会导致特征维数过大，会大大增加内存需求量和计算负担。文献[3]也只采用了主分量分析PCA或者核主分量分析KPCA来提取目标特征，但是这两种特征提取方法在用于图像的特征提取时要将二维图像矩阵转化为一维向量，这会损失图像的二维结构信息，估计得到的训练样本的协方差矩阵维数过大，带来“维数灾难”等问题。另外，主分量分析PCA假设每类目标数据均服从单模分布，例如单高斯分布，对于复杂的合成孔径雷达SAR目标图像数据来说，这一假设是很不合理的。

针对上述背景杂波的问题，我们提出了有效的预处理方法，将目标从杂波背景中分割出来，以减弱背景杂波对识别性能的影响。针对上述特征提取方法的缺陷，我们提出对角子类判决分析DiaCDA的合成孔径雷达SAR图像的特征提取方法来提取目标特征，该方法直接基于二维图像矩阵寻求最优的投影矢量，保留了更多的图像的二维空间结构信息，且估计得到的散布矩阵维数远远小于主分量分析PCA得到的协方差矩阵维数，对其进行特征分解更省时、效率更高。另外，对角子类判决分析DiaCDA假设目标数据服从多模分布，例如混合高斯分布，相对于主分量分析PCA的单模分布假设，多模分布假设更科学更合理。

[1]Zhao Q，Principe J C.Support Vector Machine for SAR Automatic Target Recognition.IEEE Trans on Aerospace and Electronic Systems，2001，37(2)：643-654.

[2]Bryant M，Garber F.SVM Classifier Applied to the MSTAR Public Data Set.SPIE，Orlando，Florida，1999，3721(4)：355-360.

[3]韩萍，吴仁彪.基于KPCA准则的SAR目标特征提取与识别.电子与信息学报，2003，10(25)：1297-1301.

发明内容

本发明的目的在于克服上述已有技术进行目标识别性能差的不足，提出一种基于对角子类判决分析DiaCDA的合成孔径雷达目标识别方法，通过对原始图像进行有效的预处理和目标图像特征提取，抑制背景杂波和降低样本维数，提高目标识别的性能。

为实现上述目的，本发明的目标识别方法有如下两种技术方案：

技术方案1

一种基于对角子类判决分析的合成孔径雷达目标识别方法，包括如下步骤：

(1)对原始的训练和测试图像依次进行阈值分割、形态学滤波、几何聚类操作、图像增强和二维傅里叶变换的预处理；

(2)采用二维快速全局k均值聚类算法对预处理后的每类目标图像进行最优的子类划分，采用对角子类判决分析DiaCDA找到最优子类的最优投影矢量；

(3)将预处理后的训练和测试图像向投影矢量进行投影，得到训练和测试目标的特征矩阵；

(4)计算测试目标与每个训练目标的特征矩阵之间的欧氏距离，并采用最近邻分类器准则确定测试目标的类别属性。

步骤2所述的“采用对角子类判决分析DiaCDA找到最优子类的最优投影矢量”，按如下步骤进行：

(a)将每类预处理后的合成孔径雷达训练目标图像采用二维全局k均值聚类算法找到最优的子类划分结果；

(b)将预处理后的合成孔径雷达训练图像沿对角方向将相应的图像像素取出来，重新排列后转化为相应的对角图像；

(c)根据上述的子类划分结果和相应的对角图像，定义如下公式：

$\left\{\begin{array}{c}{G}_B=\underset{j=1}{\overset{c-1}{\mathrm{\Σ}}}\underset{l=j+1}{\overset{c}{\mathrm{\Σ}}}\underset{i=1}{\overset{d_j}{\mathrm{\Σ}}}\underset{h=1}{\overset{d_l}{\mathrm{\Σ}}}{p}_{j,i}{p}_{l,h}{({\stackrel{\‾}{D}}^{j,i}-{\stackrel{\‾}{D}}^{l,h})}^T({\stackrel{\‾}{D}}^{j,i}-{\stackrel{\‾}{D}}^{l,h})\∈R^{n\×n}\\ G_W=\underset{j=1}{\overset{c}{\mathrm{\Σ}}}\underset{i=1}{\overset{d_j}{\mathrm{\Σ}}}\underset{k=1}{\overset{N_{j,i}}{\mathrm{\Σ}}}{(D_k^{j,i}-{\stackrel{\‾}{D}}^{j,i})}^T(D_k^{j,i}-{\stackrel{\‾}{D}}^{j,i})\∈R^{n\×n}\end{array}\right.$

其中，G_B和G_W分别表示所有训练图像对应对角图像的子类类间和子类类内散布矩阵，c为训练样本的类别数，d_j表示第j(j＝1，…，c)类目标包含的子类个数，d_l表示第l(l＝1，…，c)类目标包含的子类个数，D^j，i表示第j类的第i个子类的均值，D^l，h表示第l类的第h个子类的均值，D_k ^j，i表示I_k ^j，i的对角图像，I_k ^j，i表示第j类的第i个子类中的第k个样本，p_j，i＝N^j，i/M为第j类的第i个子类的先验概率，p_l，h＝N^l，h/M为第l类的第h个子类的先验概率，N^j，i为第j类的第i个子类中的样本个数，N^l，h为第l类的第h个子类中的样本个数，M为所有训练样本的个数，R^n×n表示子类类间和子类类内散布矩阵的大小，T表示矩阵转置；

(d)对矩阵G_W求逆后与G_B相乘，得到矩阵G_W ^-1G_B，再对它进行特征分解，得到相应的特征值和特征向量，并取前d个较大特征值对应的特征向量w₁，…，w_d组成投影矩阵W＝[w₁，…，w_d]∈R^n×d。

步骤3所述的“将预处理后的训练和测试图像向投影矢量进行投影，得到训练和测试目标的特征矩阵”，按如下步骤进行：

(a)将第i个训练样本I_i∈R^m×n向投影矩阵投影，得到它的特征矩阵B_i：

$B_i=[y_1^{\left(i\right)},\·\·\·,y_d^{\left(i\right)}]=(I_i-\stackrel{\‾}{I})[w_1,\·\·\·,w_d]=(I_i-\stackrel{\‾}{I})W\∈R^{m\×d},$ i＝1，2，…，M，

其中，y_k ⁽ⁱ⁾表示I_i∈R^m×n向w_k(k＝1，2，…，d)投影得到的特征分量，I为所有训练图像的均值图像；

(b)将测试样本I∈R^m×n向投影矩阵投影，得到它的特征矩阵B为：

B＝[y₁，…，y_d]＝(I-I)[w₁，…，w_d]＝(I-I)W∈R^m×d，

其中，y_k表示I_i∈R^m×n向w_k(k＝1，2，…，d)投影得到的特征分量。

技术方案2

一种基于对角子类判决分析和二维子类判决分析DiaCDA+2DCDA的合成孔径雷达目标识别方法，包括如下步骤：

(A)对原始的训练和测试图像依次进行阈值分割、形态学滤波、几何聚类操作、图像增强和二维傅里叶变换的预处理；

(B)采用二维快速全局k均值聚类算法对预处理后的每类目标图像进行最优的子类划分，采用对角子类判决分析DiaCDA和二维子类判决分析2DCDA找到最优子类的两个投影矩阵；

(C)将预处理后的训练和测试图像向投影矩阵进行投影，得到训练和测试目标的特征矩阵；

(D)计算测试目标与每个训练目标的特征矩阵之间的欧氏距离，并采用最近邻分类器准则确定测试目标的类别属性。

步骤B所述的“采用对角子类判决分析DiaCDA和二维子类判决分析2DCDA找到最优子类的两个投影矩阵”，按如下步骤进行：

(B1)采用对角子类判决分析DiaCDA得到投影矩阵W∈R^n×d；

(B2)采用二维子类判决分析2DCDA得到另一个投影矩阵V∈R^m×q。

步骤C所述的“将预处理后的训练和测试图像向投影矩阵进行投影，得到训练和测试目标的特征矩阵”，按如下步骤进行：

(C1)将第i个训练样本I_i∈R^m×n向投影矩阵W∈R^n×d和V∈R^m×q投影，得到它的特征矩阵C_i为：

C_i＝V^T(I_i-I)W∈R^q×d，

其中，T表示矩阵转置，I为所有训练图像的均值图像，R^q×d表示特征矩阵C_i的大小；

(C2)将测试样本I∈R^m×n向投影矩阵W∈R^n×d和V∈R^m×q投影，得到它的特征矩阵C为：

C＝V^T(I-I)W∈R^q×d

其中，T表示矩阵转置，I为所有训练图像的均值图像，R^q×d表示特征矩阵C的大小。

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

(1)现有的文献中大多都没有对原始图像做任何预处理或者只是做了简单的对数变换和傅里叶变换，由于原始图像中既包含我们感兴趣的目标又包含大量的背景杂波，如果不去除这些背景杂波，识别性能就会受到影响。在本发明中，由于对原始图像先用自适应阈值分割方法将目标从杂波背景中大致分割出来，再采用形态学滤波和几何聚类操作去除掉分割结果中出现的假目标的处理，因而减弱了背景杂波对识别性能的影响；同时由于对上述处理后图像采用基于幂变换的空域灰度变换方法增强，因而增强了目标图像质量。

(2)现有的许多文献均采用原始图像的全部像素点作为后续识别阶段的特征，由于特征维数过大，对内存需求量及运算量要求比较高，还会出现信息冗余，造成识别性能损失，因此需要对目标图像进行特征提取。在本发明中，由于采用对角子类判决分析或者对角子类判决分析和二维子类判决分析提取目标图像特征，因而提取得到了有效的目标特征，并且降低了特征维数。

(3)传统的基于向量的子空间特征提取方法，如主分量分析、线性判决分析，在用于图像的特征提取时，要将二维图像矩阵转化为一维图像向量，这会带来两个方面的问题：1.损失图像的二维空间结构信息，造成最终的识别性能损失；2.导致特征提取要在高维向量空间中进行，而在高维向量空间中，很难准确估计训练图像的协方差矩阵、类内、类间散布矩阵等，并且得到的这些矩阵维数过大，对它们进行特征分解时需要过多的内存需求和运算负担。本发明中提出的对角子类判决分析或者对角子类判决分析和二维子类判决分析的特征提取方法直接基于二维对角图像矩阵，不需要将二维图像矩阵转化为一维图像向量，其好处是：1.保留了更多的图像二维空间结构信息；2.估计得到的协方差或者散布矩阵维数远远小于基于向量的子空间特征提取方法得到的矩阵维数，更准确有效，对它们进行特征分解时效率更高。

(4)传统的特征提取方法大多都假设每类目标数据均服从单模分布，即每类目标具有相同的协方差矩阵而均值不同。然而，对于复杂的SAR目标图像数据来说，这一假设往往是不合理的。在本发明中，假设每类目标数据服从多模分布，由于采用二维快速全局k均值聚类算法找到每类目标图像的最优的子类划分，因而这一假设更科学、更合理，使得最终的目标识别性能更好。

附图说明

图1是本发明方案1的目标识别流程图；

图2是本发明预处理前后的图像对比；

图3是本发明将原始图像转化为对角图像示意图；

图4是本发明方案2的目标识别流程图；

图5是本发明预处理中幂变换的幂次对图像质量的影响；

图6是本发明不同子类个数下对角子类判决分析的识别性能；

图7是本发明与现有方法的识别性能对比图；

图8是本发明不同子类个数下对角子类判决分析和二维子类判决分析的识别性能。

具体实施方式

参照图1，本发明技术方案1的具体步骤如下。

步骤1，对原始合成孔径雷达SAR图像进行预处理。

在一个给定的场景中，通常既包含感兴趣的目标，又包含背景杂波，如图2(a)所示，如果直接对数据库中的原始图像进行特征提取及识别，背景杂波势必会影响识别性能。因此需要对其进行预处理，将目标从杂波背景中分割出来，以减弱背景杂波对识别性能的影响。此外，合成孔径雷达SAR图像中往往包含有较强的相干斑，需要进行适当处理以降低其对识别性能的影响。为此，给出了合成孔径雷达SAR图像预处理方法，其实现步骤如下：

(1)作对数变换

合成孔径雷达SAR图像通常是非高斯的，且存在大量乘性相干斑。对数变换可将相干斑由乘性模型变为加性模型，使图像更适合用高斯分布描述。设F(x，y)是原始图像的幅度矩阵F的第(x，y)个像素点，对其作对数变换，得

G(x，y)＝10lg[F(x，y)+0.001]+30.                    (1)

由于对数变换在零点没有定义，所以在变换前对原图像加一个小常数(0.001)，为了保证变换后图像灰度值仍非负，在最后加一个相应的常数(30)。

(2)自适应阈值分割

估计当前的对数图像G的均值μ和方差σ。假定目标强度基本大于背景强度，对G的每个像素点(x，y)，有

其中，T_ar(x，y)、B_ac(x，y)分别表示目标和背景的第(x，y)个像素点，c由训练样本先验地统计得到(0.9)。

(3)形态学滤波和几何聚类操作

由于相干斑的存在，阈值分割后的图像中不可避免地会有许多伪目标，在目标内部也可能会有一些“小洞”，如图2(b)所示。为了得到平滑干净的目标图像，对T_ar进行形态学滤波和几何聚类操作。

形态学滤波的目的是平滑边界，去除尖凸，切断细长的搭界，接合窄的缺口，填充凹处，去除小洞等。

一般情况下，上述处理后的图像中还可能包含一些非目标区域，与目标相比，它们小很多，如图2(c)所示。几何聚类的目的是去除这些区域，其方法是：首先检测出所有独立的连通区域，对它们进行标号；然后统计每个区域的面积，面积最大的区域即是我们感兴趣的目标，于是得到平滑干净的目标图像，如图2(d)所示。

将上述处理后的T_ar与对数图像G作像素点乘运算，得到的图像即为去掉杂波后的目标图像H。

(4)图像增强和归一化处理

本发明采用基于幂变换的灰度变换增强方法，对图像的每个像素点作幂变换，即

K(x，y)＝[H(x，y)]^α.                                (3)

其中，H(x，y)、K(x，y)分别为增强前后的目标图像的第(x，y)个像素点，α是常数。

为了去除目标相对于雷达的距离不同产生的影响，对目标图像作能量归一化处理，即

$J(x,y)=K(x,y)/\sqrt{\underset{x}{\mathrm{\Σ}}\underset{y}{\mathrm{\Σ}}{\left|K(x,y)\right|}^2}.---\left(4\right)$

其中，K(x，y)、J(x，y)分别为归一化前后的目标图像的第(x，y)个像素点。

由于目标在场景中具有位置不确定性，取傅里叶变换的幅频，以保证目标具有平移不变性。由于傅里叶变换的幅频信息是对称的，只需用一半的幅频信息作为特征提取的输入即可，这样可降低样本维数，减小运算量。

步骤2，对预处理后的合成孔径雷达SAR图像采用对角子类判决分析DiaCDA进行特征提取。

特征提取是目标识别的关键步骤之一。若用图像的全部像素点作为特征，则维数过大，对内存、运算量要求较高，还会出现信息冗余造成识别性能损失。因此，本发明方案1采用子类判决分析DiaCDA的特征提取方法提取目标特征。

假定有c类M个预处理后的训练样本为I₁，I₂，…，I_M，其中，I_i∈R^m×n，i＝1，2，…，M。对于每类预处理后的训练样本，采用二维快速全局k均值聚类算法得到最优的子类划分，如下：第j(j＝1，…，c)类目标包含有d_j个子类，且第j类的第i个子类中有N^j，i个样本，I_k ^j，i表示第j类的第i个子类中的第k个样本，I^j，i为第j类的第i个子类的均值，I为所有训练图像的均值。

对于每幅预处理后训练图像I₁，I₂，…，I_M，采用如下定义方式将其变为对角图像：

(1)当图像的高度m大于宽度n时，采用图3(a)将I_i转为其对角图像D_i；

(2)当图像的高度m小于等于宽度n时，采用图3(b)将I_i转为其对角图像D_i。

不妨假设图像高度m大于宽度n，采用图3(a)的方法将每幅预处理后的训练图像I_i转化为对角图像D_i(i＝1，…，M)。基于这些对角图像，分别定义其子类类间散布矩阵G_B和子类类内散布矩阵G_W如下：

$\left\{\begin{array}{c}{G}_B=\underset{j=1}{\overset{c-1}{\mathrm{\Σ}}}\underset{l=j+1}{\overset{c}{\mathrm{\Σ}}}\underset{i=1}{\overset{d_j}{\mathrm{\Σ}}}\underset{h=1}{\overset{d_l}{\mathrm{\Σ}}}{p}_{j,i}{p}_{l,h}{({\stackrel{\‾}{D}}^{j,i}-{\stackrel{\‾}{D}}^{l,h})}^T({\stackrel{\‾}{D}}^{j,i}-{\stackrel{\‾}{D}}^{l,h})\∈R^{n\×n}\\ G_W=\underset{j=1}{\overset{c}{\mathrm{\Σ}}}\underset{i=1}{\overset{d_j}{\mathrm{\Σ}}}\underset{k=1}{\overset{N_{j,i}}{\mathrm{\Σ}}}{(D_k^{j,i}-{\stackrel{\‾}{D}}^{j,i})}^T(D_k^{j,i}-{\stackrel{\‾}{D}}^{j,i})\∈R^{n\×n}\end{array}\right..---\left(5\right)$

其中，d_j表示第j(j＝1，…，c)类目标包含的子类个数，d_l表示第l(l＝1，…，c)类目标包含的子类个数，D^j，i表示第j类的第i个子类的均值，D^l，h表示第l类的第h个子类的均值，D_k ^j，i表示I_k ^j，i的对角图像，I_k ^j，i表示第j类的第i个子类中的第k个样本，p_j，i＝N^j，i/M为第j类的第i个子类的先验概率，p_l，h＝N^l，h/M为第l类的第h个子类的先验概率，N^j，i为第j类的第i个子类中的样本个数，N^l，h为第l类的第h个子类中的样本个数，M为所有训练样本的个数，R^n×n表示子类类间和子类类内散布矩阵的大小，T表示矩阵转置。

根据上面的散布矩阵的定义，可以看出，该方法的目的是最大化同类目标的不同子类之间的距离，同时最小化子类的散度。

则最优的投影矢量为矩阵G_W ^-1G_B的前d个较大特征值对应的特征向量w₁，…，w_d，形成投影矩阵为W＝[w₁，…，w_d]∈R^n×d。

步骤3，将预处理后的训练和测试样本向投影矩阵投影，得到其特征矩阵。

将第i个预处理后的训练样本I_i∈R^m×n向投影矩阵W＝[w₁，…，w_d]∈R^n×d进行投影，得到它的特征矩阵B_i为：

$B_i=[y_1^{\left(i\right)},\·\·\·,y_d^{\left(i\right)}]=(I_i-\stackrel{\‾}{I})[w_1,\·\·\·,w_d]=(I_i-\stackrel{\‾}{I})W\∈R^{m\×d},$ i＝1，2，…，M.            (6)

其中，y_k ⁽ⁱ⁾表示I_i∈R^m×n向w_k(k＝1，2，…，d)投影得到的特征分量，I为所有预处理后的训练图像的均值图像，R^m×d表示特征矩阵的维数。

将预处理后的测试样本I∈R^m×n向投影矩阵W＝[w₁，…，w_d]∈R^n×d进行投影，得到它的特征矩阵B为：

B＝[y₁，…，y_d]＝(I-I)[w₁，…，w_d]＝(I-I)W∈R^m×d.            (7)

其中，y_k表示I∈R^m×n向w_k(k＝1，2，…，d)投影得到的特征分量，I为所有预处理后训练图像的均值图像，R^m×d表示特征矩阵的维数。

步骤4，计算测试样本与每个训练样本的特征矩阵之间的欧氏距离，并采用最近邻准则确定测试目标的类别属性。

首先，计算预处理后的测试样本I∈R^m×n与第i(i＝1，2，…，M)个预处理后的训练样本I_i∈R^m×n的特征矩阵B和B_i之间的欧式距离为

$d(B,B_i)=\underset{k=1}{\overset{d}{\mathrm{\Σ}}}{\left|\right|y_k-y_k^{\left(i\right)}\left|\right|}_2.---\left(8\right)$

其中，向量y_k和y_k ⁽ⁱ⁾分别表示测试图像I∈R^m×n和训练I_i∈R^m×n的特征矩阵B和B_i的第k列，d表示特征矩阵B和B_i的列数，||y_k-y_k ⁽ⁱ⁾||₂表示向量y_k与y_k ⁽ⁱ⁾的2-范数。

然后，对于测试样本I∈R^m×n，经上述计算后，得到了该测试样本和M个训练样本的特征矩阵之间的距离d(B，B₁)，d(B，B₂)，…，d(B，B_M)，将使“距离最小”的训练样本所属的类别作为该测试样本的类别。

参照图4，本发明技术方案2的具体步骤如下：

步骤1，与技术方案1的步骤相同。

步骤2，对预处理后的SAR图像对角子类判决分析DiaCDA和二维子类判决分析2DCDA进行特征提取。

首先，采用上述的对角子类判决分析DiaCDA得到投影矩阵为W＝[w₁，…，w_d]∈R^n×d。

然后，基于上述的预处理后的训练样本，采用二维子类判决分析2DCDA找到另外一个投影矩阵V＝[v₁，…，v_q]，方法如下：

(1)当图像的高度m等于宽度n时，投影矩阵V＝[v₁，…，v_q]可由矩阵C_W ^-1C_B的前q个较大特征值对应的特征向量组成，并且，C_B和C_W分别为预处理后的训练图像的子类类间、子类类内散布矩阵，定义为：

$\left\{\begin{array}{c}{C}_B=\underset{j=1}{\overset{c-1}{\mathrm{\Σ}}}\underset{l=j+1}{\overset{c}{\mathrm{\Σ}}}\underset{i=1}{\overset{d_j}{\mathrm{\Σ}}}\underset{h=1}{\overset{d_l}{\mathrm{\Σ}}}{p}_{j,i}{p}_{l,h}{({\stackrel{\‾}{I}}^{j,i}-{\stackrel{\‾}{I}}^{l,h})}^T({\stackrel{\‾}{I}}^{j,i}-{\stackrel{\‾}{I}}^{l,h})\∈R^{n\×n}\\ C_W=\underset{j=1}{\overset{c}{\mathrm{\Σ}}}\underset{i=1}{\overset{d_j}{\mathrm{\Σ}}}\underset{k=1}{\overset{N_{j,i}}{\mathrm{\Σ}}}{(I_k^{j,i}-{\stackrel{\‾}{I}}^{j,i})}^T(I_k^{j,i}-{\stackrel{\‾}{I}}^{j,i})\∈R^{n\×n}\end{array}\right..---\left(9\right)$

其中，d_j表示第j(j＝1，…，c)类目标包含的子类个数，d_l表示第□(l＝1，…，c)类目标包含的子类个数，I_k ^j，i表示第j类的第i个子类中的第k个样本，I^j，i表示第j类的第i个子类的均值，I^l，h为第l类的第h个子类的均值，p_j，i＝N^j，i/M为第j类的第i个子类的先验概率，p_l，h＝N^l，h/M为第l类的第h个子类的先验概率，N^j，i为第j类的第i个子类中的样本个数，N^l，h为第l类的第h个子类中的样本个数，M为所有训练样本的个数，R^n×n表示子类类间和子类类内散布矩阵C_B和C_W的大小，T表示矩阵转置。

(2)当图像高度m不等于宽度n时，投影矩阵V＝[v₁，…，v_q]可由矩阵S_W ^-1S_B的前q个较大特征值对应的特征向量组成，并且，S_B和S_W分别为训练图像的子类类间、子类类内散布矩阵，定义为：

$\left\{\begin{array}{c}{S}_B=\underset{j=1}{\overset{c-1}{\mathrm{\Σ}}}\underset{l=j+1}{\overset{c}{\mathrm{\Σ}}}\underset{i=1}{\overset{d_j}{\mathrm{\Σ}}}\underset{h=1}{\overset{d_l}{\mathrm{\Σ}}}{p}_{j,i}{p}_{l,h}({\stackrel{\‾}{I}}^{j,i}-{\stackrel{\‾}{I}}^{l,h}){({\stackrel{\‾}{I}}^{j,i}-{\stackrel{\‾}{I}}^{l,h})}^T\∈R^{m\×m}\\ S_W=\underset{j=1}{\overset{c}{\mathrm{\Σ}}}\underset{i=1}{\overset{d_j}{\mathrm{\Σ}}}\underset{k=1}{\overset{N_{j,i}}{\mathrm{\Σ}}}(I_k^{j,i}-{\stackrel{\‾}{I}}^{j,i}){(I_k^{j,i}-{\stackrel{\‾}{I}}^{j,i})}^T\∈R^{m\×m}\end{array}\right..---\left(10\right)$

其中，d_j表示第j(j＝1，…，c)类目标包含的子类个数，d_l表示第l(l＝1，…，c)类目标包含的子类个数，I_k ^j，i表示第j类的第i个子类中的第k个样本，I^j，i表示第j类的第i个子类的均值，I^l，h为第l类的第h个子类的均值，p_j，i＝N^j，i/M为第j类的第i个子类的先验概率，p_l，h＝N^l，h/M为第l类的第h个子类的先验概率，N^j，i为第j类的第i个子类中的样本个数，N^l，h为第l类的第h个子类中的样本个数，M为所有训练样本的个数，R^m×m表示子类类间和子类类内散布矩阵S_B和S_W的大小，T表示矩阵转置。

步骤3，将预处理后的训练和测试样本向两个投影矩阵进行投影，得到其特征矩阵。

将第i个预处理后的训练样本I_i∈R^m×n向投影矩阵W∈R^n×d和V∈R^m×q投影，得到它的特征矩阵C_i为：

C_i＝V^T(I_i-I)W∈R^q×d.                            (11)

其中，I_i∈R^m×n为第i个预处理后的训练样本，I为所有预处理后的训练图像的均值图像，T表示矩阵转置，R^q×d表示特征矩阵C_i的大小；

将预处理后的测试样本I∈R^m×n向投影矩阵W∈R^n×d和V∈R^m×q投影，得到它的特征矩阵C为：

C＝V^T(I-I)W∈R^q×d.                            (12)

其中，I∈R^m×n为预处理后的测试图像，I为所有预处理后的训练图像的均值图像，T表示矩阵转置，R^q×d表示特征矩阵C的大小。

步骤4，计算测试样本与每个训练样本的特征矩阵之间的欧氏距离，并采用最近邻准则确定测试目标的类别属性。

首先，按如下三种方式定义预处理后的测试样本I∈R^m×n与第i(i＝1，2，…，M)个预处理后的训练样本I_i∈R^m×n的特征矩阵C和C_i之间的欧式距离为：

(1)沿行定义距离：不妨将特征矩阵C和C_i写成如下形式：C＝[x₁，…，x_q]^T， $C_i={[x_i^{\left(i\right)},\·\·\·,x_q^{\left(i\right)}]}^T.$ 定义两特征矩阵之间的距离为

$d_1(C,C_i)=\underset{k_1=1}{\overset{q}{\mathrm{\Σ}}}{\left|\right|x_{k_1}-x_{k_1}^{\left(i\right)}\left|\right|}_2.---\left(13\right)$

其中，

和

分别表示测试样本I∈R^m×n和第i个训练样本I_i∈R^m×n的特征矩阵C和C_i的第k₁个行向量，q表示特征矩阵C和C_i的行数。

(2)沿列定义距离：将特征矩阵C和C_i写成如下形式：C＝[y₁，…，y_d]， $C_i=[y_1^{\left(i\right)},\·\·\·,y_d^{\left(i\right)}].$ 定义两特征矩阵之间的距离为

$d_2(C,C_i)=\underset{k_2=1}{\overset{d}{\mathrm{\Σ}}}{\left|\right|y_{k_2}-y_{k_2}^{\left(i\right)}\left|\right|}_2.---\left(14\right)$

其中，和分别表示测试样本I∈R^m×n和第i个训练样本I_i∈R^m×n的特征矩阵C和C_i的第k₂个列向量，d表示特征矩阵C和C_i的列数。

(3)沿行和列定义距离：综合上述两种距离定义形式，则待测图像与第i个训练图像I∈R^m×n和I_i∈R^m×n的特征矩阵C和C_i之间距离为：

d(C，C_i)＝d₁(C，C_i)+d₂(C，C_i).                    (15)

然后，对于测试样本I∈R^m×n，按行定义形式计算距离，得到了该测试样本和M个训练样本的特征矩阵之间的距离d₁(C，C₁)，d₁(C，C₂)，…，d₁(C，C_M)，将使“距离最小”的训练样本所属的类别作为该测试样本的类别。

或者，对于测试样本I∈R^m×n，按列定义形式计算距离，得到了该测试样本和M个训练样本的特征矩阵之间的距离d₂(C，C₁)，d₂(C，C₂)，…，d₂(C，C_M)，将使“距离最小”的训练样本所属的类别作为该测试样本的类别。

或者，对于测试样本I∈R^m×n，按行和列定义距离形式计算，得到了该测试样本和M个训练样本的特征矩阵之间的距离d(C，C₁)，d(C，C₂)，…，d(C，C_M)，将使“距离最小”的训练样本所属的类别作为该测试样本的类别。

在本发明的两个方案的特征提取方法中，假定每类目标的子类结构是已知的，但实际上是不知道的。因此，需要采用聚类算法以找到每类目标的子类结构。在本发明的两个方案中，采用的聚类算法均是二维快速全局k均值聚类算法，它是将快速全局k均值聚类算法推广到二维空间中。在二维快速全局k均值聚类算法中，需要确定每类目标可能包含的子类个数，为了方便起见，假设每类目标具有相同的子类数目k_max。

仿真实验

为了验证本发明提出的合成孔径雷达SAR图像预处理方法和特征提取方法，即对角子类判决分析DiaCDA、对角子类判决分析和二维子类判决分析DiaCDA+2DCDA的有效性，做了如下实验。

表1  实验采用的训练和测试数据

实验采用的是美国运动和静止目标获取与识别MSTAR计划录取的实测合成孔径雷达SAR地面静止目标数据，如表1所示，其中“-”表示没有用到。训练样本是合成孔径雷达SAR在俯仰角17°时对地面目标的成像数据，包括BMP2sn_9563(233个)、BTR70sn_c71(233个)、T72sn_132(232个)。测试样本是合成孔径雷达SAR在俯仰角15°时对地面目标的成像数据，包括BMP2sn_c21(196个)、BMP2sn_9563(195个)、BMP2sn_9566(196个)、BTR70sn_c71(196个)、T72sn_132(196个)、T72sn_812(195个)、T72sn_s7(191个)。BMP2和T72都有不同的变体，不同的变体间在配置和序列号上各不相同。所有样本的方位覆盖范围是0°～360°，所有原始图像大小为128×128。预处理后图像大小变为128×64。

以T72为例分析预处理中对数变换的作用及图像增强时幂次不同对目标图像质量的影响。图5(a)～(h)分别给出原始图像、对数变换后图像、二值图像、图5(c)与图5(b)掩模后的图像、幂次为2时图5(d)增强后的图、幂次为3时图5(d)增强后的图像、幂次为4时图5(d)增强后的图像和幂次为5时图5(d)增强后的图像。由于原始的合成孔径雷达SAR图像中存在乘性相干斑，且整体灰度值较低，目标的许多细节不可见，如图5(a)所示。对数变换，一方面将相干斑由乘性变为加性，使图像更适合用高斯分布描述便于阈值分割，另一方面扩展目标的灰度值，显示出更多的细节。但目标内部的对比度却下降了，如图5(d)所示，因此对比度增强是必需的，这可由幂次α＞1的幂变换实现。从图5(e)～(h)看出，α从2增大到4，图像对比度明显增强，但当α再增大时，目标中有些地方太暗，从而丢失一些细节。因此得出，对比度与可辨别细节的最好增强效果在α＝4左右。由于幂次α的选择对识别性能的影响很难从数学上进行准确描述，在实际应用中可基于训练和测试数据进行预设α的值(如不指明，以下均取α＝3.5)。

图6给出了不同子类个数下，对角子类判决分析DiaCDA的最高识别率。从该图可以看出，子类数目从1到10变化时，识别率都能达到97％左右，并且当子类数目为2时，识别性能最优。

为了进一步验证本发明的预处理方法和对角子类判决分析DiaCDA的特征提取方法的有效性，图7给出了不同的特征参数d，即不同的特征维数为128×d下，结合本发明提出的预处理方法，对角子类判决分析DiaCDA与现有的特征提取方法，如二维主分量分析2DPCA、二维线性判决分析2DLDA、二维子类判决分析2DCDA的识别率对比。从该图不难看出，二维子类判决分析2DCDA和对角子类判决分析DiaCDA的识别性能优于二维主分量分析2DPCA和二维线性判决分析2DLDA，这是因为：后者假设每类目标是线性可分的，而前者充分利用了每类目标中包含的子类信息，因此后者可以缓解前者的线性不可分问题。另外，由于对角子类判决分析DiaCDA通过对角图像构造子类类间、子类类内散布矩阵，而现有的二维子类判决分析2DCDA特征提取方法仅仅利用了图像的行信息，因此，对角子类判决分析DiaCDA的识别性能稍优于二维子类判决分析2DCDA。

表2  不同特征提取方法的最高识别率对比

表2也给出了本发明提出的特征提取方法与现有的常用的特征提取方法的最高识别率及特征维数的对比。从该表可以看出，对角子类判决分析DiaCDA的性能优于主分量分析PCA、二维主分量分析2DPCA、二维线性判决分析2DFLD和二维子类判决分析2DCDA。

另外，从表2也可以看出，尽管二维主分量分析2DPCA、二维线性判决分析2DLDA、二维子类判决分析2DCDA和对角子类判决分析DiaCDA均取得了相当好的识别性能，但是它们得到的特征维数较大，这会导致内存需求量较大、识别速度较慢。

为了降低特征维数、改善识别性能，将对角子类判决分析DiaCDA和二维子类判决分析2DCDA结合起来给出DiaCDA+2DCDA的特征提取方法。为了验证DiaCDA+2DCDA的有效性，图8给出了在不同子类个数下，DiaCDA+2DCDA的最高识别率。从该图不难看出，在三种不同的距离定义形式下，DiaCDA+2DCDA取得了相当好的性能，且最优识别率高达98.01％。表2也证明了DiaCDA+2DCDA在大大降低了特征维数的同时，又改善了识别性能。

另外，由表2还看出，尽管现有的主分量分析PCA的特征维数最小，但是它在用于图像的特征提取时，需要将二维图像矩阵转化为一维向量，这会导致两个方面的问题：1.损失图像的二维空间结构信息；2.特征提取要在高维向量空间中进行，带来“维数灾难”等问题。

Claims (9)

1.一种基于对角子类判决分析的合成孔径雷达目标识别方法，包括如下步骤：

(1)对原始的训练和测试图像依次进行阈值分割、形态学滤波、几何聚类操作、图像增强和二维傅里叶变换的预处理；

(2)采用二维快速全局k均值聚类算法对预处理后的每类目标图像进行最优的子类划分，采用对角子类判决分析找到最优子类的最优投影矢量；

(3)将预处理后的训练和测试图像向投影矢量进行投影，得到训练和测试目标的特征矩阵；

(4)计算测试目标与每个训练目标的特征矩阵之间的欧氏距离，并采用最近邻分类器准则确定测试目标的类别属性。

2.根据权利要求1所述的目标识别方法，其中步骤1所述的“形态学滤波”主要采用了matlab工具箱中bwmorph函数，去除阈值分割后图像中孤立的像素点、尖凸、小洞、接合缺口和平滑边界。

3.根据权利要求1所述的目标识别方法，其中步骤1所述的“几何聚类操作”按如下步骤进行：

(3a)检测出上述形态学滤波后的图像中所有独立的连通区域，并对它们进行标号；

(3b)统计每个区域的面积；

(3c)将面积最大的区域作为感兴趣的目标。

4.根据权利要求1所述的目标识别方法，其中步骤1所述的“图像增强”是对图像的每个像素点采用如下基于幂变换的空域灰度变换公式进行：

K(x，y)＝[H(x，y)]^α

其中，H(x，y)、K(x，y)分别为增强前后的目标图像的第(x，y)个像素点，α是常数。

5.根据权利要求1所述的目标识别方法，其中步骤2所述的“采用对角子类判决分析找到最优子类的最优投影矢量”，按如下步骤进行：

(5a)将每类预处理后的合成孔径雷达训练目标图像采用二维全局k均值聚类算法找到最优的子类划分结果；

(5b)将预处理后的合成孔径雷达训练图像沿对角方向将相应的图像像素取出来，重新排列后转化为相应的对角图像；

(5c)根据上述的子类划分结果和相应的对角图像，定义如下公式：

$\left\{\begin{array}{c}{G}_B=\underset{j=1}{\overset{c-1}{\mathrm{\Σ}}}\underset{l=j+1}{\overset{c}{\mathrm{\Σ}}}\underset{i=1}{\overset{d_j}{\mathrm{\Σ}}}\underset{h=1}{\overset{d_I}{\mathrm{\Σ}}}{p}_{j,i}{p}_{l,h}{({\stackrel{\‾}{D}}^{j,i}-{\stackrel{\‾}{D}}^{l,h})}^T({\stackrel{\‾}{D}}^{j,i}-{\stackrel{\‾}{D}}^{l,h})\∈R^{n\×n}\\ G_W=\underset{j=1}{\overset{c}{\mathrm{\Σ}}}\underset{i=1}{\overset{d_j}{\mathrm{\Σ}}}\underset{k=1}{\overset{N_{j,i}}{\mathrm{\Σ}}}{(D_k^{j,i}-{\stackrel{\‾}{D}}^{j,i})}^T(D_k^{j,i}-{\stackrel{\‾}{D}}^{j,i})\∈R^{n\×n}\end{array}\right.$

其中，G_B和G_W分别表示所有训练图像对应对角图像的子类类间和子类类内散布矩阵，c为训练样本的类别数，d_j表示第j(j＝1，…，c)类目标包含的子类个数，d_l表示第l(l＝1，…，c)类目标包含的子类个数，D^j，i表示第j类的第i个子类的均值，D^l，h表示第l类的第h个子类的均值，K_k ^j，i表示I_k ^j，i的对角图像，I_k ^j，i表示第j类的第i个子类中的第k个样本，p_j，i＝N^j，i/M为第j类的第i个子类的先验概率，p_l，h＝N^l，h/M为第l类的第h个子类的先验概率，N^j，i为第j类的第i个子类中的样本个数，N^l，h为第l类的第h个子类中的样本个数，M为所有训练样本的个数，R^n×n表示子类类间和子类类内散布矩阵的大小，T表示矩阵转置；

(5d)对矩阵G_W求逆后与G_B相乘，得到矩阵G_W ^-1G_B，再对它进行特征分解，得到相应的特征值和特征向量，并取前d个较大特征值对应的特征向量w₁，…，w_d组成投影矩阵W＝[w₁，…，w_d]∈R^n×d。

6.根据权利要求1所述的目标识别方法，其中步骤3所述的“将预处理后的训练和测试图像向投影矢量进行投影，得到训练和测试目标的特征矩阵”，按如下步骤进行：

(6a)将第i个训练样本I_i∈R^m×n向投影矩阵投影，得到它的特征矩阵B_i为：

$B_i=[y_1^{\left(i\right)},\·\·\·,y_d^{\left(i\right)}]=(I_i-\stackrel{\‾}{I})[w_1,\·\·\·,w_d]=(I_i-\stackrel{\‾}{I})W\∈R^{m\×d},i=\mathrm{1,2},\·\·\·,M,$

其中，y_k ⁽ⁱ⁾表示I_i∈R^m×n向w_k(k＝1，2，…，d)投影得到的特征分量，I为所有训练图像的均值图像；

(6b)将测试样本I∈R^m×n向投影矩阵投影，得到它的特征矩阵B为：

B＝[y₁，…，y_d]＝(I-I)[w₁，…，w_d]＝(I-I)W∈R^m×d，

其中，y_k表示I_i∈R^m×n向w_k(k＝1，2，…，d)投影得到的特征分量。

7.一种基于对角子类判决分析和二维子类判决分析的合成孔径雷达目标识别方法，包括如下步骤：

(A)对原始的训练和测试图像依次进行阈值分割、形态学滤波、几何聚类操作、图像增强和二维傅里叶变换的预处理；

(B)采用二维快速全局k均值聚类算法对预处理后的每类目标图像进行最优的子类划分，采用对角子类判决分析和二维子类判决分析找到最优子类的两个投影矩阵；

(C)将预处理后的训练和测试图像向投影矩阵进行投影，得到训练和测试目标的特征矩阵；

(D)计算测试目标与每个训练目标的特征矩阵之间的欧氏距离，并采用最近邻分类器准则确定测试目标的类别属性。

8.根据权利要求7所述的目标识别方法，其中步骤B所述的“采用对角子类判决分析和二维子类判决分析找到最优子类的两个投影矩阵”，按如下步骤进行：

(B1)采用对角子类判决分析得到投影矩阵W∈R^n×d；

(B2)采用二维子类判决分析得到另一个投影矩阵V∈R^m×q。

9.根据权利要求7所述的目标识别方法，其中步骤C所述的“将预处理后的训练和测试图像向投影矩阵进行投影，得到训练和测试目标的特征矩阵”，按如下步骤进行：

(C1)将第i个训练样本Iⁱ∈R^m×n向投影矩阵W∈R^n×d和V∈R^m×q投影，得到它的特征矩阵C_i为：

C_i＝V^T(I_i-I)W∈R^q×d，

其中，T表示矩阵转置，I为所有训练图像的均值图像，R^q×d表示特征矩阵C_i的大小；

(C2)将测试样本I∈R^m×n向投影矩阵W∈R^n×d和V∈R^m×q投影，得到它的特征矩阵C为：

C＝V^T(I-I)W∈R^q×d，

其中，T表示矩阵转置，I为所有训练图像的均值图像，R^q×d表示特征矩阵C的大小。

Images