基于张量MPCA的极化SAR图像分类方法与流程

本发明属于图像处理技术领域，更进一步涉及目标识别技术领域中的一种基于张量MPCA的极化合成孔径雷达(Synthetic Aperture Radar，SAR)图像分类方法。本发明可用于对雷达图像中不同目标区域进行分类。

背景技术：

雷达是一种可以实现全天候工作的主动探测系统，它可以穿透一定的地表，并且可以改变发射波的频率、强度。极化SAR是用来测量回波信号的新型雷达，可以记录不同极化状态组合回波的相位差信息，所以能够得到更丰富的目标信息，在农业、林业、军事、地质学、水文学和海洋等方面具有广泛的研究和应用价值，如地物种类的识别、灾害监测和评估、农作物成长监视等。现有的极化SAR图像分类方法大致可以分为基于散射特性的方法和基于统计特性方法，如基于目标分解的分类方法和基于散射熵、散射角的分类方法等。

电子科技大学在其申请的专利“一种基于Cloude分解的极化SAR图像分类方法”(专利申请号：201410341457.9，公开号：104123563A)中提出了一种利用Cloud分解的SAR图像分类的方法。该方法首先对图像中的每个像素点进行Cloude分解，得到散射熵H和散射角α；计算散射熵H和散射角α的统计直方图，并对直方图进行分割，得到分割阈值；将得到的分割阈值作为散射熵H和散射角α特征构成的二维平面的划分点，对极化SAR图像进行初始划分为m类；将得到的初始分类的类中心和类别数输入Wishart分类器，得到极化SAR图像的分类结果。该方法比原始的H/α法硬性划分要更为合理，能根据图像本身的特点自适应的调节。但是，该方法的不足之处是，需要将原始数据转换成矢量的形式，破坏了原始数据中的空间结构信息，没有充分利用图像信息，使得该方法需要较多的训练样本来训练分类器。

武汉大学在其申请的专利“一种弱监督的SAR图像分类方法”(专利申请号：201010221605.5，公开号：101894275A)中提出了一种弱监督的SAR图像分类方法。该方法先将SAR图像分割成多个子图像，然后将各子图像格网划分成互不重叠的矩形区域，然后将每个矩形区域采用一个建立在多尺度信息上的局部四叉树描述，然后建立分级的马尔科夫主题模型，通过该模型对图像进行分类。该方法能同时利用图像的数据信息，图像与图像之间的相关性息和图像在多个尺度上的相关信息，因而很好的克服了SAR图像分类中的一些局部不确定性问题。但是，该方法仍然存在的不足之处是，该方法的抗噪性不强，在噪声较大的情况下分类器性能下降比较明显。

技术实现要素：

本发明的目的在于解决上述现有技术的不足，提出了一种基于张量MPCA的极化SAR图像分类方法，相比现有方法将原始数据转换成向量形式，本发明将原始数据转换成张量形式，保留了原始数据三个通道之间的空间结构信息，充分利用了极化SAR图像三个通道的数据，在提高抗噪声性能的同时，降低了对训练样本数目的要求提高了分类的准确度。

本发明实现上述目的的思路是：先将极化SAR图像原始数据转换成张量形式，然后利用多线性主成分张量分析MPCA方法对图像特征进行提取，再训练支撑矢量机SVM分类器，利用训练好的支撑矢量机SVM对待分类数据进行分类，得到分类结果。

本发明的步骤包括如下：

(1)输入数据：

输入极化合成孔径雷达图像SAR三个通道的数据；

(2)数据转换：

(2a)对每个通道的像素取3×3的邻域窗口，得到每个像素对应的特征矩阵；

(2b)将三个通道数据中对应像素的特征矩阵构成3×3×3的三阶张量，作为每个像素的特征张量；

(3)提取图像特征：

(3a)根据极化合成孔径雷达图像SAR图像的实际地物分布，对每一类数据随机选取训练样本；

(3b)利用多线性主成分张量分析MPCA方法，对每一类数据的训练样本进行特征提取，得到每一类数据训练样本的特征向量；

所述的多线性主成分张量分析MPCA方法的具体步骤如下：

第一步，按照下式，对每一类数据的训练样本进行中心化处理：

其中，表示中心化处理后每个像素的特征张量，表示每一类数据的训练样本中每个像素的特征张量，m表示每一类数据的训练样本像素点的序号，m的取值范围为1、2、…M；M表示所选取每一类数据的训练样本的数量Σ表示求和操作；

第二步，按照下式，计算所有像素特征张量的协方差矩阵：

其中,φ⁽ⁿ⁾表示所有像素特征张量的协方差矩阵，Σ表示求和操作，表示中心化处理后每个像素的特征张量的n模展开矩阵，n表示每个像素特征张量的阶数，n的取值范围为1、2、3，m表示每一类数据的训练样本像素点的序号，表示中心化处理后每个像素的特征张量的n模展开矩阵的转置矩阵，m的取值范围为1,2,…M，M表示所选每一类数据的训练样本的数量；

第三步，利用奇异值分解方法，对所有像素特征张量的协方差矩阵进行特征值分解，得到所有像素特征张量的协方差矩阵的特征值及其与特征值相对应的特征向量，将特征向量按照相对应的特征值大小从大到小排序，选择前6个特征值对应的特征向量组成每一类数据的训练样本的投影矩阵U_(n)，其中，U_(n)表示每一类数据的训练样本的投影矩阵，n表示每个像素特征张量的阶数，n的取值范围为：1、2、3；

第四步，设定迭代次数k，其中，迭代次数的取值范围为[10，100]内的正整数；

第五步，按照下式，计算每一类数据的训练样本的变换矩阵：

其中，Φ⁽ⁿ⁾表示每一类数据的训练样本的变换矩阵，Σ表示求和操作，表示中心化处理后每个像素的特征张量的n模展开矩阵，n表示每个像素特征张量的阶数，n的取值范围为1、2、3，m表示每一类数据的训练样本像素点的序号，m的取值范围为1,2…M，M表示所选每一类数据的训练样本的数量，表示中心化处理后每个像素的特征张量的n模展开矩阵的转置矩阵，U_(n)表示每一类数据的训练样本的投影矩阵，表示做克罗内克积操作；

第六步，利用奇异值分解方法，对每一类数据的训练样本的变换矩阵进行特征值分解，得到变换矩阵的特征值及其与特征值相对应的特征向量，将特征向量按照相对应的特征值大小从大到小排序，选择前6个特征值对应的特征向量更新每一类数据的训练样本的投影矩阵；

第七步，判断是否达到设定的迭代次数，如果是，执行第八步，否则，执行第五步；

第八步，按照下式，将每一类数据的训练样本进行张量分解，得到每一类数据训练样本的特征向量：

其中，y_m表示每一类数据的训练样本的特征向量，X_m表示每一类数据的训练样本中每个像素的特征张量，m表示每一类数据的训练样本像素点的序号，m的取值范围为1,2,…M，M表示所选取每一类数据的训练样本的数量，×_n表示做n阶乘积操作，n表示每个像素特征张量的阶数，n的取值范围为：1、2、3，表示每一类数据的训练样本的投影矩阵的转置矩阵；

(4)训练分类器：

利用每一类数据训练样本的特征向量训练支撑矢量机支撑SVM分类器，得到训练好的支撑矢量机SVM分类器；

(5)图像分类：

使用训练好的支撑矢量机SVM分类器，对待分类数据进行分类，得到分类结果。

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

第一，由于本发明将极化合成孔径雷达SAR图像三个通道的原始数据转换成张量形式加以利用，保留了原始数据中的空间结构信息，克服了现有技术没有充分挖掘极化合成孔径雷达SAR图像原始数据三个通道空间结构信息的不足，降低了对训练数据数目的要求，使得本发明可以用更少的训练数据得到更准确的分类结果。

第二，由于本发明采用多线性主成分张量分析MPCA方法提取图像特征，利用了三个通道数据间的相关性，能够保持不同地物间的散射特性，克服了现有技术对相干斑噪声敏感的不足，使得本发明提高了抗噪声的性能，具有较强的鲁棒性。

附图说明

图1为本发明的流程图；

图2为本发明的仿真图。

具体实施方式

下面结合附图，对本发明做进一步的详细描述。

参照附图1，本发明的具体步骤如下。

步骤1.输入数据。

输入极化合成孔径雷达图像SAR三个通道的数据，具体包括水平发射水平接收HH通道、水平发射垂直接收HV通道和垂直发射垂直接收VV通道，三个通道的散射强度值。

步骤2.数据转换。

对每个通道的像素取3×3的邻域窗口，得到每个像素对应的特征矩阵，然后将三个通道数据中对应像素的特征矩阵构成3×3×3的三阶张量，作为每个像素的特征张量。

步骤3.提取图像特征。

首先，根据极化合成孔径雷达图像SAR图像的实际地物分布，对每一类数据选取训练样本，然后，利用多线性主成分张量分析MPCA方法对每一类训练样本进行特征提取，得到每一类训练数据的特征向量，多线性主成分张量分析MPCA方法是一种用来在所有张量模式方向上寻求它们主要成分的多线性算法，具体步骤如下。

第一步，按照下式，对每一类数据的训练样本进行中心化处理：

其中，表示中心化处理后每个像素的特征张量，X_m表示每一类数据的训练样本中每个像素的特征张量，m表示每一类数据的训练样本像素点的序号，m的取值范围为1、2、…M；M表示所选取每一类数据的训练样本的数量，Σ表示求和操作。

第二步，按照下式，计算所有像素特征张量的协方差矩阵：

其中,φ⁽ⁿ⁾表示所有像素特征张量的协方差矩阵，Σ表示求和操作，表示中心化处理后每个像素的特征张量的n模展开矩阵，n表示每个像素特征张量的阶数，n的取值范围为1、2、3，m表示每一类数据的训练样本像素点的序号，表示中心化处理后每个像素的特征张量的n模展开矩阵的转置矩阵，m的取值范围为1,2,…M，M表示所选每一类数据的训练样本的数量。

第三步，利用奇异值分解方法，对所有像素特征张量的协方差矩阵进行特征值分解，得到所有像素特征张量的协方差矩阵的特征值及其与特征值相对应的特征向量，将特征向量按照相对应的特征值大小从大到小排序，一般情况下，前10％甚至1％的特征值的和就占了全部的特征值总和的99％以上，因此我们可以用前10％甚至1％的特征值来近似的描述所有像素特征张量的协方差矩阵，我们按照保留所有像素特征张量的协方差矩阵总信息量的97％来选取特征值个数，选择最大的6个特征值对应的特征向量组成每一类数据的训练样本的投影矩阵U_(n)，其中，U_(n)表示每一类数据的训练样本的投影矩阵，n表示每个像素特征张量的阶数，n的取值范围为：1、2、3。

第四步，设定迭代次数k，其中，迭代次数的取值范围为[10，100]内的正整数。

第五步，按照下式，计算每一类数据的训练样本的变换矩阵：

其中，Φ⁽ⁿ⁾表示每一类数据的训练样本的变换矩阵，Σ表示求和操作，表示中心化处理后每个像素的特征张量的n模展开矩阵，n表示每个像素特征张量的阶数，n的取值范围为1、2、3，m表示每一类数据的训练样本像素点的序号，m的取值范围为1,2…M，M表示所选每一类数据的训练样本的数量，表示中心化处理后每个像素的特征张量的n模展开矩阵的转置矩阵，U_(n)表示每一类数据的训练样本的投影矩阵，表示做克罗内克积操作。

第六步，利用奇异值分解方法，对每一类数据的训练样本的变换矩阵进行特征值分解，得到变换矩阵的特征值及其与特征值相对应的特征向量，将特征向量按照相对应的特征值大小从大到小排序，一般情况下，前10％甚至1％的特征值的和就占了全部的特征值总和的99％以上，因此我们可以用前10％甚至1％的特征值来近似的描述所有像素特征张量的协方差矩阵，我们按照保留训练样本的变换矩阵总信息量的97％来选取特征值个数，这里选择最大的6个特征值对应的特征向量更新每一类数据的训练样本的投影矩阵。

第七步，判断是否达到设定的迭代次数，如果是，执行第八步，否则，执行第五步。

第八步，按照下式，将每一类数据的训练样本进行张量分解，得到每一类数据训练样本的特征向量：

其中，y_m表示每一类数据的训练样本的特征向量，X_m表示每一类数据的训练样本中每个像素的特征张量，m表示每一类数据的训练样本像素点的序号，m的取值范围为1,2,…M，M表示所选取每一类数据的训练样本的数量，×_n表示做n阶乘积操作，n表示每个像素特征张量的阶数，n的取值范围为：1、2、3，表示每一类数据的训练样本的投影矩阵的转置矩阵。

步骤4.训练分类器。

利用训练样本的特征向量训练支撑矢量机支撑SVM分类器，得到训练好的支撑矢量机SVM分类器。

通过求解下面的优化问题，得到最优分类超平面：

其中，max表示取最大值操作，Σ表示求和操作，a_i表示拉格朗日乘子，M表示所选每一类数据训练样本的数量，y_i,y_j表示每一类数据训练样本的类别，x_i,x_j表示每一类数据的训练样本，i，j表示每一类训练样本像素点的序号，i，j的取值范围为1,2,…M，<·，·>表示做欧式内积操作，s.t.表示约束条件。

步骤5.图像分类。

使用训练好的支撑矢量机SVM分类器，对待分类数据进行分类，得到分类结果。

本发明的效果可通过以下仿真进一步说明

1.仿真条件：

本发明的仿真实验在主频2.5GHZ的Pentium Dual_Core CPU E5200、内存4.00GB的硬件环境和MATLAB R2014a的软件环境下进行的。

2.仿真内容：

本发明仿真使用的实验数据为美国旧金山地区的L波段极化SAR数据，图2(a)为本发明仿真使用的旧金山地区L波段极化SAR数据的RGB合成图。图2(b)为使用现有技术的H/α分类方法，对本发明仿真实验数据的分类结果图。图2(c)为使用现有技术的主成分分析方法，对本发明仿真实验数据的分类结果图。图2(d)为利用本发明对实验数据的分类结果图。

3.仿真效果分析：

由图2(b)可见，使用现有技术的H/α分类方法，地物可以大致得到划分，但是还有较多地域划分不准确，而且对于边界的划分过于武断。由图2(c)可见，使用现有技术的主成分分析方法，分类结果杂点较多，左上角山体和水域无法准确区分。由图2(d)可见，本发明的分类结果与图2(b)、图2(c)相比，分类结果中杂点数量明显降低，地物信息体现得更加精细，山脉，植被和城区也变得更易区分，其中左上方的山峰、跑马场均可辨识出来。