引言

合成孔径雷达(Synthetic aperture radar, SAR)图像目标识别的一般过程为：图像预处理、特征提取和分类识别。分类识别的前提和关键在于特征提取。合成孔径雷达图像的特征提取方法主要分成两类：一类是在SAR图像域的特征提取，常用特征有：峰值、分形维数、阴影和纹理等^[1]；另一类是基于数学变换的特征提取：主要是利用数学变换方法将变换域的系数作为图像的特征，常见的方法包括主成分分析、小波变换、离散余弦变换、独立成分分析和线性判决分析等。

主成分分析(Principal component analysis, PCA)也称为离散K-L变换，是一种通过特征线性组合的数据降维方法，由于其通过线性变换后产生的新特征正交，变换后的新特征更稳定、能量更集中，使得该方法被广泛应用于SAR图像特征提取中^[2-5]。主成分分析面向的数组是向量形式，在降维过程中首先需要将图像数据转化为一维向量，导致特征提取要在高维向量空间中进行，很难准确估计协方差矩阵且维数很大，并且丢失了图像本身的空间结构信息。有学者将PCA特征提取和其他特征提取方法相结合来提高识别率和识别速度。李映等^[6]提出一种基于核的奇异值分解(Kernel singular value decomposition，KSVD)与主成分分析相结合的SAR图像目标的组合特征提取方法，首先利用核的奇异值分解得到图像非线性的代数特征，然后进一步经过PCA变换得到图像的最终分类特征，实验结果表明，该方法不仅有效地提高了目标的正确识别率，而且大大降低了对目标方位的敏感度。贺志国等^[7]从实时性的角度出发，采用PCA特征提取和神经网络相结合的SAR目标识别方法，实验结果表明，在维持较高识别率的前提下，该类方法具有内存需求少、运行速度快的优点，可用于实时处理。李勇等^[8]应用在PCA基础发展起来的二维PCA(2DPCA)方法进行SAR特征提取，该方法对SAR图像进行二层小波分解后提取低频子带图像的二维主成分分量作为目标的识别特征，利用支持向量机完成目标分类，取得了不错的识别效果。胡利平等^[9]提出一种两级2DPCA的图像特征提取方法，进一步压缩特征维数，减少识别运算量，实验结果表明，该方法能有效降低特征维数，提高识别率，并且对目标方位角具有较强的鲁棒性。尽管上述特征提取方法获得了较高的识别率或者识别速度，但这些方法都不可避免地破坏了样本的空间结构^[10]。多线性主成分分析(Multilinear principal component analysis, MPCA)在PCA的基础上直接使用张量的形式提取有效特征，避免了“维数灾难”的同时又具有处理张量数据时保留原始数据空间结构信息的优点^[11-13]。为了克服当前SAR图像目标识别中PCA特征提取使得图像空间结构信息丢失的问题，进一步提高目标正确识别率，本文提出一种基于多线性主成分分析和张量分析的SAR图像目标识别方法，该方法对SAR图像构建四阶张量样本，并采用多线性主成分分析提取特征，有效地保留了图像的空间结构信息。

1 多线性主成分分析

多线性主成分分析旨在寻求一组多线性投影矩阵，通过多线性投影矩阵，把原始张量映射到核心张量子空间中，其详细步骤如下^[9]：

(1) 首先对原始张量样本X∈R^{I₁×I₂×…×I_N}进行中心化处理，即每个张量样本减去张量样本集的平均值为

$ {{\mathit{\boldsymbol{\bar X}}}_m} = {\mathit{\boldsymbol{X}}_m} - \mathit{\boldsymbol{\bar X}} $

(1)

$ \mathit{\boldsymbol{\bar X}} = \frac{1}{M}\sum\limits_{m = 1}^M {{\mathit{\boldsymbol{X}}_m}} $

(2)

式中：M为样本数；X为样本均值；X_m为原始样本；X_m为中心化后的样本。

(2) 初始化协方差矩阵

$ {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\left( n \right)}} = \sum\limits_{m = 1}^M {\left( {{{\mathit{\boldsymbol{\bar X}}}_{m\left( n \right)}}{{\mathit{\boldsymbol{\bar X}}}_m}{{\left( n \right)}^{\rm{T}}}} \right)} $

(3)

式中：X_m(n)为张量样本X_m的n模展开后的矩阵，对该协方差矩阵进行特征分解，取最大Pⁿ个特征值对应的Pⁿ个特征向量组成投影矩阵U⁽ⁿ⁾∈R^I_n×P_n，n=1, 2, …, N。

(3) 局部优化

① 利用投影矩阵进行张量子空间映射，计算

$ {\mathit{\boldsymbol{S}}_m} = {{\mathit{\boldsymbol{\bar X}}}_m} \times {}_1{\mathit{\boldsymbol{U}}^{\left( 1 \right){\rm{T}}}} \times {}_2{\mathit{\boldsymbol{U}}^{\left( 2 \right){\rm{T}}}} \times \cdots \times {}_N{\mathit{\boldsymbol{U}}^{\left( N \right){\rm{T}}}} $

(4)

式中：m=1, 2, …, M；×_n, n=1, 2, …, N为n模乘积。

② 计算

$ {\mathit{\Psi }_{{S_0}}} = \sum\limits_{m = 1}^M {{{\left\| {{\mathit{\boldsymbol{S}}_m}} \right\|}^2}} $

(5)

式中||S_m||为张量的范数。

③ 循环以下步骤：

a.对于张量样本的每一阶方向n=1, 2, …, N，求协方差矩阵Φ⁽ⁿ⁾的特征分解，取最大Pⁿ个特征值对应的Pⁿ个特征向量组成投影矩阵U⁽ⁿ⁾∈R^I_n×P_n, n=1, 2, …, N，其中

$ {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\left( n \right)}} = \sum\limits_{m = 1}^M {\left( {{{\mathit{\boldsymbol{\bar X}}}_{m\left( n \right)}}{\mathit{\boldsymbol{U}}_{{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\left( n \right)}}}}\mathit{\boldsymbol{U}}_{{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\left( n \right)}}}^{\rm{T}}\mathit{\boldsymbol{\bar X}}_{m\left( n \right)}^{\rm{T}}} \right)} $

(6)

$ {\mathit{\boldsymbol{U}}_{{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\left( n \right)}}}} = {\mathit{\boldsymbol{U}}^{\left( {n + 1} \right)}} \otimes \cdots \otimes {\mathit{\boldsymbol{U}}^{\left( N \right)}} \otimes {\mathit{\boldsymbol{U}}^{\left( 1 \right)}} \otimes \cdots \otimes {\mathit{\boldsymbol{U}}^{\left( {n - 1} \right)}} $

(7)

式中：⊗为矩阵的克罗内克(Kronecker)乘积。

b.计算S_m和Ψ_Sk

$ {\mathit{\boldsymbol{S}}_m} = {{\mathit{\boldsymbol{\bar X}}}_m} \times {}_1{\mathit{\boldsymbol{U}}^{\left( 1 \right){\rm{T}}}} \times {}_2{\mathit{\boldsymbol{U}}^{\left( 2 \right){\rm{T}}}} \times \cdots \times {}_N{\mathit{\boldsymbol{U}}^{\left( N \right){\rm{T}}}} $

(8)

$ {\mathit{\Psi }_{{S_k}}} = \sum\limits_{m = 1}^M {{{\left\| {{\mathit{\boldsymbol{S}}_m}} \right\|}^2}} $

(9)

c.如果Ψ_Sk－Ψ_Sk－1 < η，退出循环，否则回到步骤a，其中η为自定义的阈值。

(4) 投影计算最终降维后的核心张量

$ {\mathit{\boldsymbol{S}}_m} = {{\mathit{\boldsymbol{\bar X}}}_m} \times {}_1{\mathit{\boldsymbol{U}}^{\left( 1 \right){\rm{T}}}} \times {}_2{\mathit{\boldsymbol{U}}^{\left( 2 \right){\rm{T}}}} \times \cdots \times {}_N{\mathit{\boldsymbol{U}}^{\left( N \right){\rm{T}}}} $

(10)

2 基于MPCA和张量分析的SAR图像目标识别

基于MPCA和张量分析的SAR图像目标识别方法，首先对SAR训练图像构建四阶张量训练样本，用MPCA对训练张量降维，得到投影矩阵，然后将原始样本向张量子空间进行映射得到核心张量，再通过线性判别分析(Linear discriminant analysis，LDA)获得分类函数权向量，最后进行测试样本的分类识别。

2.1 构建四阶张量训练样本

首先对原始的SAR图像进行预处理，将图像统一调整为大小为128像素×128像素的幅度图像P，并对幅度数据进行归一化处理，处理后图像幅度值的均值为0，方差为1，归一化公式为

$ P = \left( {P - {\rm{mean}}P} \right)/{\rm{std}}P $

(11)

式中:meanP为图像幅度数据的均值，stdP为标准差。

按照图像空间x轴、图像空间y轴、方位角和样本类别，构建一个四阶的张量训练样本，该四阶张量可以表示为X_m∈R^{I₁×I₂×I₃×I₄}，I₁, I₂, I₃, I₄分别代表图像空间x轴、图像空间y轴、方位角和样本类别四阶张量的维度，四阶张量构造示意图如图 1所示。

2.2 多线性主成分分析获得多线性投影矩阵

多线性主成分分析首先对原始的四阶张量训练样本进行中心化，即每个张量样本减去张量样本集的均值为

$ {{\mathit{\boldsymbol{\bar X}}}_m} = {\mathit{\boldsymbol{X}}_m} - \mathit{\boldsymbol{\bar X}} $

(12)

$ \mathit{\boldsymbol{\bar X}} = \frac{1}{M}\sum\limits_{m = 1}^M {{\mathit{\boldsymbol{X}}_m}} $

(13)

然后对中心化后的四阶张量样本沿各阶模展开到一组矩阵X_n－MODE，该过程对组成张量的所有阶按交错次序采样，在整个采样过程中对不同阶的特征值进行混合交错采样，从而实现了张量不同阶特征值之间的传递和融合。

最后对四阶张量展开后的矩阵进行多线性主成分分析，获得各阶上的投影矩阵U⁽ⁿ⁾, n=1, 2, 3, 4。

2.3 构造核心张量 2.3.1 核心张量的计算

与奇异值分解类似，U⁽ⁿ⁾, n=1, 2, 3, 4可以被看成是一组正交变换基对，而核心张量S_m是X_m在该变换基对上的投影，中心化后的训练张量样本可以表示为

$ {{\mathit{\boldsymbol{\bar X}}}_m} = {\mathit{\boldsymbol{S}}_m} \times {}_1{\mathit{\boldsymbol{U}}^{\left( 1 \right)}} \times {}_2{\mathit{\boldsymbol{U}}^{\left( 2 \right)}} \times {}_3{\mathit{\boldsymbol{U}}^{\left( 3 \right)}} \times {}_4{\mathit{\boldsymbol{U}}^{\left( 4 \right)}} $

(14)

将中心化后的样本通过投影矩阵向核心张量子空间进行映射得到核心张量S_m，表达式为

$ {\mathit{\boldsymbol{S}}_m} = {{\mathit{\boldsymbol{\bar X}}}_m} \times {}_1{\mathit{\boldsymbol{U}}^{\left( 1 \right){\rm{T}}}} \times {}_2{\mathit{\boldsymbol{U}}^{\left( 2 \right){\rm{T}}}} \times {}_3{\mathit{\boldsymbol{U}}^{\left( 3 \right){\rm{T}}}} \times {}_4{\mathit{\boldsymbol{U}}^{\left( 4 \right){\rm{T}}}} $

(15)

S_m可以很好地表征中心化后的张量训练样本X_m的特性，从而可以用该核心张量取代原始的张量训练样本进行后续的训练与识别。

2.3.2 核心张量特征的选择

将映射后的核心张量特征转化为一维数组，计算每一个特征分量的类内离散度和类间离散度以及类间离散度和类内离散度的比值

$ \begin{array}{*{20}{c}} {{{\mathit{\boldsymbol{\bar V}}}_k} = \frac{1}{M}\sum\limits_{m = 1}^{{M_k}} {{\mathit{\boldsymbol{V}}_m}} }&{k = 1,2, \cdots ,K} \end{array} $

(16)

$ {{\mathit{\boldsymbol{\bar V}}}_A} = \frac{1}{K}\sum\limits_{k = 1}^K {{{\mathit{\boldsymbol{\bar V}}}_k}} $

(17)

$ \begin{array}{*{20}{c}} {S{w_k} = \sum\limits_{m = 1}^{{M_k}} {{{\left( {{\mathit{\boldsymbol{V}}_m} - {{\mathit{\boldsymbol{\bar V}}}_k}} \right)}^2}} }&{k = 1,2, \cdots ,K} \end{array} $

(18)

$ \begin{array}{*{20}{c}} {S{b_k} = {{\left( {{{\mathit{\boldsymbol{\bar V}}}_k} - {{\mathit{\boldsymbol{\bar V}}}_A}} \right)}^2}}&{k = 1,2, \cdots ,K} \end{array} $

(19)

$ R = \frac{{\sum\limits_{k = 1}^K {S{b_k}} }}{{\sum\limits_{k = 1}^K {S{w_k}} }} $

(20)

式中：K为类别数；M_k为每类样本的样本数；V_m为当前第k类的样本；V_k为每类样本的均值；V_A为全体样本的均值；Sw_k为每类样本的类内离散度；Sb_k为每类样本的类间离散度；R为类间离散度和类内离散度的比值。

对核心张量向量化后特征的每一特征分量按照R大小进行降序排序，选取R最大的若干特征分量进行后续的训练和识别，根据式(21)计算选取的特征分量数占总特征分量数的比例为

$ 特征分量比例 = \frac{{选取特征分量数}}{{总特征分量数}} $

(21)

2.4 线性判别分析

核心张量可表征原始训练样本，因此可以通过核心张量训练线性判别分类器；对于K类的线性分类问题要寻找K个线性判别函数，每个线性判别函数的求解目的是寻找一个权向量，使得训练样本X的错分最小，因而，这里需要求得K个权向量W_k，k=1, 2, …, K，K为类别数；定义误差向量为

$ {\mathit{\boldsymbol{e}}_k} = \mathit{\boldsymbol{X}}{\mathit{\boldsymbol{W}}_k} - {b_k} $

(22)

式中：b_k为已知的样本类别。

将目标函数定义为平方误差的形式为

$ J\left( {{\mathit{\boldsymbol{W}}_k}} \right) = {\left\| {{\mathit{\boldsymbol{e}}_k}} \right\|^2} = {\left\| {\mathit{\boldsymbol{X}}{\mathit{\boldsymbol{W}}_k} - {b_k}} \right\|^2} $

(23)

W_k的优化目标为使得J(W_k)最小，即求J(W_k)的梯度为0，表达式为

$ \nabla J\left( {{\mathit{\boldsymbol{W}}_k}} \right) = \sum\limits_{i = 1}^M {2\left( {\mathit{\boldsymbol{W}}_k^{\rm{T}}{\mathit{\boldsymbol{X}}_i} - {b_{k,i}}} \right){x_i}} = 2{\mathit{\boldsymbol{X}}^{\rm{T}}}\left( {\mathit{\boldsymbol{X}}{\mathit{\boldsymbol{W}}_k} - {b_k}} \right) = 0 $

(24)

从而得到权向量为

$ {\mathit{\boldsymbol{W}}_k} = {\left( {{\mathit{\boldsymbol{X}}^{\rm{T}}}\mathit{\boldsymbol{X}}} \right)^{ - 1}}{\mathit{\boldsymbol{X}}^{\rm{T}}}{b_k} $

(25)

2.5 测试样本分类

对于给定的测试样本X_test，首先将测试样本通过投影矩阵U⁽ⁿ⁾, n=1, 2, 3, 4映射到张量子空间中，得到测试样本核心张量为

$ {\mathit{\boldsymbol{S}}_{{\rm{test}}}} = {\mathit{\boldsymbol{X}}_{{\rm{test}}}} \times {}_1{\mathit{\boldsymbol{U}}^{\left( 1 \right){\rm{T}}}} \times {}_2{\mathit{\boldsymbol{U}}^{\left( 2 \right){\rm{T}}}} \times {}_3{\mathit{\boldsymbol{U}}^{\left( 3 \right){\rm{T}}}} \times {}_4{\mathit{\boldsymbol{U}}^{\left( 4 \right){\rm{T}}}} $

(26)

通过线性判别函数的权向量W_k计算，即

$ gn{d_k}\left( {{\mathit{\boldsymbol{S}}_{{\rm{test}}}}} \right) = \mathit{\boldsymbol{W}}_k^{\rm{T}}{\mathit{\boldsymbol{S}}_{{\rm{test}}}} $

(27)

比较函数值gnd_k的大小，使gnd_k最大的线性判别函数权向量所对应的类别则是测试样本的类别。

3 实验结果与分析 3.1 实验数据集

采用MSTAR数据集中的实测SAR地面静止目标数据来测试本文方法，并比较PCA算法、2DPCA算法和MPCA算法在SAR图像特征提取上的性能。MSTAR是当前SAR目标识别性能评估的公开数据库，该数据库来自于美国国防预研计划署和空军研究实验室(DARPA/AFRL)共同资助的运动与静止目标的获取与识别计划。实验中选用SAR在俯视角为17°和15°的7类目标图像数据作为目标的训练和测试样本，该7类目标分别为：BTR70_c71, D7, ZSU_23/4, BRDM_2, T72_132, BTR_60和2S1。图 2是7类目标在俯视角17°及不同方位角φ下的SAR原始图像，其训练和测试样本情况见表 1。

3.2 实验结果与分析

本文采用识别率作为SAR目标识别效果的衡量指标，它是正确识别样本数和总样本数的比值。图 3给出了本文提出的MPCA特征提取和LDA分类方法(MPCA+LDA)选取不同数量的特征分量时，分别采用3类、5类和7类目标作为训练和测试样本集得到的识别率。图 3中，特征分量比例从10%增加到100%，3条识别率曲线首先随着特征分量比例的增加而增加，然后随着特征分量比例的进一步增加而趋于稳定。由图 3可知，采用3类目标进行训练测试，当特征分量比例为30%时，识别率最高为96.77%。采用5类目标进行训练测试，当特征分量比例为90%时，识别率最高为96.46%。采用7类目标进行训练测试，当特征分量比例为100%时，识别率最高为92.34%。

分别取3类、5类和7类目标作为训练和测试样本集时的最高识别率，得到的混淆矩阵如表 2—4所示。从表 2—4可见，BTR70目标的类间离散度较高，该类正确分类的几率高，被误判为其他类别的几率很少；而D7和ZSU_23/4的特征相似度较高，在特征空间中距离很近，导致这两类被相互误判的概率较高。从表 3和表 4中数据可知，BRDM_2, BTR_60和2S1在识别的过程中和其他类别目标都发生了误判的情况，说明这3类在特征空间中分布较为广泛。

将本文提出的MPCA特征提取和LDA分类方法(MPCA+LDA)得到的识别率与用PCA特征提取和LDA分类方法(PCA+LDA)以及2DPCA特征提取和LDA分类方法(2DPCA+LDA)得到的识别率进行比较，比较结果如表 5—7所示。由表 5—7可知，在分别采用3类、5类和7类目标作为训练和测试样本集的3组实验中，本文提出的MPCA+LDA方法的平均识别率均要高于PCA+LDA和2DPCA+LDA方法的识别率，这是由于本文方法对图像样本构建四阶张量样本，并采用多线性主成分分析提取特征，有效地保留了图像空间结构信息。因此，目标识别率得到了明显的提高。

将MPCA+LDA方法所需的计算时间与PCA+LDA方法以及2DPCA+LDA方法所需的计算时间进行比较，比较结果如表 8所示。由于在实际工程应用中，训练过程都是离线进行的，因而本文只对识别单幅SAR目标图像所需的计算时间进行比较。由表 8可见，尽管MPCA+LDA方法的计算时间高于PCA+LDA和2DPCA+LDA方法，但MPCA+LDA方法的计算时间还是在毫秒级，完全可满足实际工程应用的实时性需求。

4 结束语

本文对SAR图像构建四阶张量训练样本，并通过多线性主成分分析对SAR图像张量样本进行特征提取，然后对降维后的特征进行线性判别分析，进行目标的分类识别。实验结果表明，该方法可有效保留图像的空间结构信息，较PCA，2DPCA特征提取方法可明显提高SAR目标的正确识别率。