CN104463193A - 基于深度稀疏ica的极化sar图像分类方法 - Google Patents

Abstract

本发明公开了一种基于深度稀疏ICA的极化SAR图像分类方法，主要克服现有技术中提取的极化SAR图像特征单一且特征层次低造成分类精度不高的问题。其实现步骤为：先对待分类的极化SAR图像进行精致Lee滤波；再提取滤波后图像的散射、纹理、颜色及协方差矩阵特征并将所有特征进行归一化和白化处理，之后选取无标签样本集，训练样本集，测试样本集；利用深度稀疏ICA方法提取训练样本集和测试样本集原始特征的深度稀疏特征表示；最后训练分类器并预测分类，输出结果。本发明相对于已有的方法，图像的特征描述更全面细致，提高了分类精度，可用于极化合成孔径雷达SAR图像目标的检测和目标识别。

Description

基于深度稀疏ICA的极化SAR图像分类方法

技术领域

本发明属于图像处理技术领域，更进一步涉及极化合成孔径雷达SAR图像分类方法，可用于对极化SAR图像的地物分类和目标识别。

背景技术

极化SAR雷达能够获得丰富的目标信息，在农业、林业、军事、地质学、水文学和海洋等方面具有广泛的研究和应用价值，如地物种类的识别、农作物成长监视、产量评估、地物分类、海冰监测、地物沉降监测，目标检测和海洋污染检测等。而极化SAR图像分类又是极化SAR图像处理的中的重要内容，也是极化SAR图像解译的关键技术之一，它是将解译系统中的前端部分单独提取出来作为具体应用的一个典型实例。快速、准确的极化SAR图像分类是实现各种实际应用的前提。

首都师范大学在其专利申请“极化SAR图像分类方法”(专利申请号：201310685303.7，授权公告号：CN103617427A)中提出了一种基于极化SAR图像散射熵H、反熵A和散射角α的极化SAR分类方法。该方法首先提取极化SAR图像的散射熵H、反熵A和散射角α并将其作为第一特征集，然后将所述极化SAR图像分解为两个子孔径图像，分别提取所述两个子孔径图像的特征：散射熵H、反熵A和散射角α从而得到两个子特征集(H₁,A₁,α₁)、(H₂,A₂,α₂)，并将所述的两个子特征集中的各相应特征的值相减，得到所述各相应特征的差值的集合(ΔH,ΔA,Δα)作为第二特征集，最后将所述第一特征集和所述的第二特征集输入到决策树分类模型中，得到所述极化SAR图像的分类结果。该专利技术虽然不仅利用了从极化SAR图像本身提取的特征集合，同时利用了两个子孔径图像提取的特征各相应的差值的集合，即两种不同类型的特征集合，使得所提取的特征相比于以往仅仅提取原极化SAR图像的散射特征更加丰富，但是该技术仍然仅仅只是利用了单一的散射特征，而忽略了更加丰富的图像的底层特征如纹理特征、颜色特征、空间位置关系等，这在一定程度上限制了极化SAR图像分类精度的提高。

西安电子科技大学在其专利申请“基于SDIT和SVM的极化SAR图像分类方法”(专利申请号：201410089692.1，授权公告号：CN103824084A)中提出了一种基于SDIT和SVM的极化SAR图像分类方法。该方法首先对待分类的极化SAR图像进行精致Lee滤波，然后提取图像的散射、偏振及纹理特征，将得到的特征进行组合并归一化，最后用归一化后的特征训练分类器，预测分类并计算分类精度。该专利技术虽然充分利用了极化SAR图像包含的丰富的纹理信息和偏振信息，弥补了仅仅依靠散射特性对极化SAR图像分类的不足，在一定程度上提高了分类的精度，但仍然存在的不足是，该专利技术只是将所提取的散射、偏振和纹理信息进行简单的堆叠，然后输入SVM直接用于分类，这导致输入的特征中包含较多的冗余信息，并且特征维数过高，使得分类效率大大下降。

发明内容

本发明的目的在于针上对述现有技术的不足，提出一种基于深度稀疏ICA的极化SAR图像分类方法。以充分利用极化SAR的散射特性和图像的底层特征，提高分类的精度和分类效率。

本发明的技术方案是：对待分类的极化SAR图像进行精致Lee滤波，提取图像的散射、纹理、颜色及协方差矩阵特征作为原始特征；利用深度稀疏ICA方法提取原始特征的深度稀疏特征表示；用得到的深度稀疏特征表示训练分类器，用训练好的分类器对整幅待分类极化SAR图像进行地物分类，对分类后的极化SAR图像进行上色输出。其实现步骤包括如下：

(1)输入待分类的极化合成孔径雷达SAR图像T1和待分类SAR图像的真实地物标记图像T2；

(2)采用精致Lee滤波方法，对待分类的极化SAR图像T1进行滤波，去除斑点噪声，得到滤波后的极化SAR图像T1’；

(3)提取滤波后的极化SAR图像T1’的散射、纹理、颜色及协方差矩阵特征；

(4)将散射特征、纹理特征、颜色特征及协方差矩阵特征组合成原始样本特征集，并将该原始样本特征集的值归一化到0～1之间，进行白化处理；

(5)选取无标签样本集，训练样本集和测试样本集；

5a)在白化处理过的原始样本特征集中随机选取2％的样本特征作为无标签样本集X_{un_0}；

5b)根据待分类SAR图像的真实地物标记图像T2，在白化处理过的原始样本特征集中选出与该T2中已标记真实地物相同的样本特征集，并在该样本特征集中随机选取20％的样本特征作为训练样本集X_{tr_0}，剩余的80％的样本特征作为测试样本集X_{te_0}；

(6)用深度稀疏ICA方法，分别提取无标签样本集、训练样本集和测试样本集的深度稀疏特征表示；

6a)用均值为0，方差为1的正态分布初始化第一层稀疏ICA的基矩阵W₁；

6b)利用基矩阵W₁，计算无标签样本集的方差衰减值，获取第一层稀疏ICA的最优基矩阵W_{Opt_1}；

6c)将获得的最优基矩阵W_{Opt_1}分别与无标签样本集X_{un_0}，训练样本集X_{tr_0}，测试样本集X_{te_0}相乘，得到无标签样本集，训练样本集，测试样本集的第一层深度稀疏特征表示X_{un_1}，X_{tr_1}，X_{te_1}；

6d)将由第一层稀疏ICA提取到的无标签样本集的第一层深度稀疏特征表示X_{un_1}作为第二层稀疏ICA的输入，重复6a)～6c)获得第二层稀疏ICA的最优基矩阵W_{Opt_2}，进而得到无标签样本集，训练样本集，测试样本的第二层深度稀疏特征表示X_{un_2}，X_{tr_2}，X_{te_2}，

其中，稀疏ICA的层数可以据分类结果进行调整，每加一层，都要重复上述6a)～6c)一次，得到训练样本集和测试样本集最终的深度稀疏特征表示分别是X_tr和X_te；

(7)将训练样本集的最终深度稀疏特征表示X_tr输入到支持向量机分类器中进行训练，得到训练好的分类器；

(8)利用训练好的分类器，对待分类的极化SAR图像T1的每个像素进行分类，得到极化SAR图像T1的每个像素所属地物类别；

(9)对步骤(8)中的分类结果，根据红、绿、蓝三基色原理，用同一种颜色标出同类地物，得到上色后的极化SAR图像T3，输出上色后的极化SAR图像T3。

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

1.本发明充分利用了极化合成孔径雷达SAR图像的散射、纹理、颜色和协方差矩阵的特性，提取的特征更加丰富多样；

2.本发明在提取的原始特征的基础上，利用深度稀疏ICA方法提取原始特征的深度稀疏特征表示，该深度稀疏特征表示相比于原始特征，冗余性小，提高了分类效率；

3.本发明用提取的深度稀疏特征表示分类，比原始的特征直接输入分类器进行分类精度更高。

附图说明

图1是本发明的实现流程图；

图2是本发明实验使用的待分类极化SAR图像；

图3是本发明实验所用的待分类极化SAR图像真实地物标记图像；

图4是用Wishart有监督分类法对待分类极化SAR图像的分类结果图；

图5是用本发明提取的原始特征对待分类极化SAR图像的分类结果；

图6是用本发明提取的深度稀疏特征对待分类极化SAR图像的分类结果。

具体实施方式

参照图1，本发明的具体实现步骤如下：

步骤1，输入图像。

输入两幅极化合成孔径雷达SAR图像：一幅为待分类的极化合成孔径雷达SAR图像T1，另一幅为带有真实地物标记的极化合成孔径雷达SAR图像T2。

步骤2，对待分类图像进行滤波处理。

采用精致Lee滤波方法，对待分类的极化SAR图像T1进行滤波，去除斑点噪声，得到滤波后极化SAR图像T1’。

步骤3，提取散射、纹理、颜色及协方差矩阵特征。

3a)取滤波后的极化SAR图像T1’的每个像素点的协方差矩阵C，该协方差矩阵的大小为3×3，C的表示形式如下：

$C=\left[\begin{array}{ccc}\⟨{\left|S_{\mathrm{HH}}\right|}^2\⟩& \sqrt{2}\⟨S_{\mathrm{HH}}{S_{\mathrm{HV}}}^{*}\⟩& \⟨S_{\mathrm{HH}}{S_{\mathrm{VV}}}^{*}\⟩\\ \sqrt{2}\⟨S_{\mathrm{HV}}{S_{\mathrm{HH}}}^{*}\⟩& 2\⟨{\left|S_{\mathrm{HV}}\right|}^2\⟩& \sqrt{2}\⟨S_{\mathrm{HV}}{S_{\mathrm{VV}}}^{*}\⟩\\ \⟨S_{\mathrm{VV}}{S_{\mathrm{HH}}}^{*}\⟩& \sqrt{2}\⟨S_{\mathrm{VV}}{S_{\mathrm{HV}}}^{*}\⟩& \⟨{\left|S_{\mathrm{VV}}\right|}^2\⟩\end{array}\right],$

其中，H表示水平极化方向，V表示垂直极化方向，S_HH表示水平向接收水平向发射的回波数据，S_HV表示垂直向接收水平向发射的回波数据，S_VV表示垂直向接收垂直向发射的回波数据。〈·〉表示按视数平均，*表示共轭，|·|表示取模；

3b)由协方差矩阵C得到相干矩阵：T＝UCU^T

其中，U表示中间变量 $\frac{1}{\sqrt{2}}\left[\begin{array}{ccc}1& 0& 1\\ 1& 0& -1\\ 0& \sqrt{2}& 0\end{array}\right],$ (·)^T表示转置操作；

3c)提取滤波后极化SAR图像T1’的散射特征：

对滤波后极化SAR图像T1’的每个像素，采用Yamaguchi分解方法得到p_s，p_v，p_d，p_h，f_s，f_v，f_d，f_h共8个散射参数，采用Cloude分解方法得到A，H，α，λ₁，λ₂，λ₃共6个散射参数，采用Huynen分解方法得到|2A₀|，|C-iD|，|H+iG|共3个散射参数，将得到的这些散射参数作为17维的散射特征；

3d)提取滤波后极化SAR图像T1’的纹理特征：

由滤波后极化SAR图像T1’每个像素点的协方差矩阵，得到该像素点的水平向接收水平向发射的回波数据|S_HH|²，垂直向接收水平向发射的回波数据|S_HV|²，垂直向接收垂直向发射的回波数据|S_VV|²，分别由每个像素点的这三个极化值|S_HH|²，|S_HV|²，|S_VV|²形成和滤波后的极化SAR图像T1’大小相同的三幅极化通道图像T4、T5、T6；计算每幅极化通道图像的灰度共生矩阵，根据灰度共生矩阵提取每幅极化通道图像每个像素的能量、逆差距、对比度三个纹理特征；将三幅极化通道图像得到的9个纹理特征，作为滤波后的极化SAR图像的9维的纹理特征。

以极化通道图像T4为例：

按照下式，得到极化通道图像T4的灰度共生矩阵：

G(i,j|θ)＝[(x,y),(x+dx,y+dy)]，

其中，G(i,j|θ)表示极化通道图像T4的灰度共生矩阵，i，j分别是极化通道图像T4的两个像素的像素值，θ表示极化通道图像T4的坐标偏移量的方向，分别取0度，45度，90度，135度共四个方向，x表示极化通道图像T4像素的横坐标，y表示极化通道图像T4像素的纵坐标，dx表示极化通道图像T4像素的横坐标偏移量，dy表示极化通道图像T4像素的纵坐标偏移量；

按照下式，分别得到极化通道图像T4的能量、逆差矩和对比度特征；

能量： $\mathrm{En}=\underset{i=0}{\overset{P-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{P-1}{\mathrm{\Σ}}}{\left(G(i,j)\right)}^2,$

逆差矩： $\mathrm{Idm}=\underset{i=0}{\overset{P-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{P-1}{\mathrm{\Σ}}}\frac{G(i,j)}{1+{(i-j)}^2},$

对比度： $\mathrm{Con}=\underset{n=0}{\overset{P-1}{\mathrm{\Σ}}}{n}^2\left\{\underset{|i-j|=n}{\mathrm{\Σ}}G(i,j)\right\},$

其中： $u_i=\underset{i=0}{\overset{P-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{P-1}{\mathrm{\Σ}}}\mathrm{iG}(i,j),u_j=\underset{i=0}{\overset{P-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{P-1}{\mathrm{\Σ}}}\mathrm{jG}(i,j),S_i^2=\underset{i=0}{\overset{P-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{P-1}{\mathrm{\Σ}}}G(i,j){(i-u_i)}^2,$ $S_j^2=\underset{i=0}{\overset{P-1}{\mathrm{\Σ}}}\underset{j=0}{\overset{P-1}{\mathrm{\Σ}}}G(i,j){(j-u_j)}^2$

对三幅极化通道图像T4、T5、T6中的每一幅图像都要分别在0度，45度，90度，135度共四个方向求相应的能量、逆差距和对比度；然后取四个方向的能量、逆差距和对比度平均值作为相应极化通道图像最终的能量、逆差距和对比度特征，得到三幅极化通道图像T4、T5、T6的纹理特征E_HH，I_HH，C_HH，E_HV，I_HV，C_HV，E_VV，I_VV，C_VV，将得到的这些值作为9维的纹理特征；

3e)提取滤波后极化SAR图像T1’的颜色特征：

根据待分类SAR图像T1及其真实地物标记图像T2，从每类地物中随机选择三个像素，并求这三个像素的平均R，G，B颜色值作为该类地物颜色的R，G，B分量，得到K个RGB颜色值，其中，K为待分类极化SAR图像T1的地物类别数目；令这K个RGB颜色值为聚类中心，对于滤波后的极化SAR图像T1’上的每个像素，计算其9*9的矩形窗口邻域中每个像素的RGB颜色值I＝(r,g,b)到各个聚类中心的欧式距离；根据最短距离原则，分别把窗口内的每个像素归到最小距离对应的地物类别中，再计算窗口内被归到各个地物类别的像素个数占窗口内总像素数目的比值，用得到的K个比值作为该窗口中心像素的K维颜色特征；

3f)提取滤波后极化SAR图像T1’的散射协方差矩阵特征：

对滤波后的极化SAR图像T1’的协方差矩阵，计算下面的各值：

HH-VV相关系数的幅度： ${\mathrm{\ρ}}_{\mathrm{HH}-\mathrm{VV}}=\left|\frac{\⟨S_{\mathrm{HH}}{S_{\mathrm{VV}}}^{*}\⟩}{\sqrt{{\left|S_{\mathrm{HH}}\right|}^2{\left|S_{\mathrm{VV}}\right|}^2}}\right|$

HH-VV相位差：Φ_HH-VV＝arg(〈S_HHS_VV ^*〉)

HV/VV通道比： $10*\log \left(\frac{{\left|S_{\mathrm{HV}}\right|}^2}{{\left|S_{\mathrm{VV}}\right|}^2}\right)$

同极化比： $10*\log \left(\frac{{\left|S_{\mathrm{VV}}\right|}^2}{{\left|S_{\mathrm{HH}}\right|}^2}\right)$

交叉极化： $10*\log \left(\frac{{\left|S_{\mathrm{HV}}\right|}^2}{{\left|S_{\mathrm{HH}}\right|}^2}\right)$

同极化系数比： $\mathrm{\γ}=\frac{\⟨S_{\mathrm{VV}}{S}_{\mathrm{VV}}^{*}\⟩}{\⟨S_{\mathrm{HH}}{S}_{\mathrm{HH}}^{*}\⟩}$

去极化系数比： $\mathrm{\δ}=\frac{\⟨S_{\mathrm{HV}}{S_{\mathrm{HV}}}^{*}\⟩}{\⟨S_{\mathrm{HH}}{S}_{\mathrm{HH}}^{*}\⟩+\⟨S_{\mathrm{VV}}{S}_{\mathrm{VV}}^{*}\⟩}$

C(1,2)的角度： ${\mathrm{\θ}}_1=\arctan \frac{\mathrm{imag}\left(C\left(\mathrm{1,2}\right)\right)}{\mathrm{real}\left(C\left(\mathrm{1,2}\right)\right)}$

C(1,3)的角度： ${\mathrm{\θ}}_2=\arctan \frac{\mathrm{imag}\left(C\left(\mathrm{1,3}\right)\right)}{\mathrm{real}\left(C\left(\mathrm{1,3}\right)\right)}$

C(2,3)的角度： ${\mathrm{\θ}}_3=\arctan \frac{\mathrm{imag}\left(C\left(\mathrm{2,3}\right)\right)}{\mathrm{real}\left(C\left(\mathrm{2,3}\right)\right)}$

C(1,1)的模值：abs(C(1,1))

C(1,2)的模值：abs(C(1,2))

C(1,3)的模值：abs(C(1,3))

C(2,2)的模值：abs(C(2,2))

C(2,3)的模值：abs(C(2,3))

C(3,3)的模值：abs(C(3,3))

其中，H表示水平极化方向，V表示垂直极化方向，C(a,b)分别表示协方差矩阵C第a行第b列的元素。将得到的这16个值作为16维的协方差矩阵特征。

步骤4，特征归一化并白化。

将散射特征、纹理特征、颜色特征和协方差矩阵特征组合成原始样本特征集，并将该原始特征集的值归一化到0～1之间，进行白化处理。

步骤5，选取无标签样本集，训练样本集和测试样本集。

5a)在白化处理过的原始样本特征集中随机选取2％的样本特征作为无标签样本集X_{un_0}；

5b)根据待分类SAR图像的真实地物标记图像T2，在白化处理过的原始样本特征集中选出与该T2中已标记真实地物相同的样本特征集，并在该样本特征集中随机选取20％的样本特征作为训练样本集X_{tr_0}，剩余的80％的样本特征作为测试样本集X_{te_0}。

步骤6，用深度稀疏ICA方法，分别提取无标签样本集、训练样本集和测试样本集的深度稀疏特征表示。

6a)用均值为0，方差为1的正态分布初始化第一层稀疏ICA的基矩阵W₁；

6b)利用基矩阵W₁，计算整体无标签样本集的均方差衰减值，获取第一层稀疏ICA的最优基矩阵W_{Opt_1}；

6b1)根据第一层稀疏ICA，计算第一层深度特征表示的稀疏度：

$P(W_1,X_{\mathrm{un}\_0})=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}\underset{j=1}{\overset{t}{\mathrm{\Σ}}}g\left(W_{1\_j}{x}^{\left(i\right)}\right),$

其中，W₁表示选取的第一层稀疏ICA的基矩阵，其初始值服从均值为0，方差为1的正态分布，W_{1_j}表示基矩阵W₁中的第j个基向量，x⁽ⁱ⁾表示第i个无标签样本，g表示非线性凸函数，t表示基矩阵W₁中基向量的数目，m表示无标签样本集中样本数；

6b2)计算无标签样本集与第一层稀疏ICA的输出W₁ ^TW₁xⁱ之间的误差：

$J_1(W_1,X_{\mathrm{un}\_0})=\frac{\mathrm{\λ}}{m}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\left|\right|{W_1}^T{W}_1{x}^{\left(i\right)}-x^{\left(i\right)}\left|\right|}_2^2,$

其中，λ为平衡参数，用于对误差J₁(W₁,X_{un_0})和稀疏度P(W₁,X_{un_0})进行平衡，W₁x⁽ⁱ⁾为第一层深度稀疏特征表示，W₁ ^TW₁xⁱ为第一层深度稀疏特征表示W₁x⁽ⁱ⁾对无标签样本x⁽ⁱ⁾进行恢复的结果；

6b3)根据步骤6b1)和6b2)的计算结果，计算整体无标签样本集的均方差衰减值：

J(W₁,X_{un_0})＝J₁(W₁,X_{un_0})+P(W₁,X_{un_0})，

6b4)对步骤6b3)得到的整体无标签样本集的均方差衰减值，采用如下的梯度下降公式，获取每一次迭代时的第一层稀疏ICA的基矩阵：

${W_1}^{n+1}={W_1}^n-\mathrm{\β}\frac{\∂}{\∂{W_1}^n}J({W_1}^n,X_{\mathrm{un}\_0}),$

其中，n表示迭代的次数，W₁ ⁿ表示第n次迭代时第一层稀疏ICA的基矩阵，W₁ ⁿ⁺¹表示第n+1次迭代时第一层稀疏ICA的基矩阵，β表示对第一层稀疏ICA基矩阵进行学习的速率，β的取值范围为0<β<1，表示对第n次迭代时第一层稀疏ICA的基矩阵求偏导数；

6b5)将步骤6b4)每一次迭代得到的第一层稀疏ICA的基矩阵值代入步骤6b3)公式中求整体无标签样本的均方差衰减值J(W₁,X_{un_0})，当整体无标签样本的均方差衰减值J(W₁,X_{un_0})达到全局最小值时，停止步骤6b4)的迭代，并将停止迭代时第一层稀疏ICA的基矩阵作为第一层稀疏ICA的最优基矩阵W_{Opt_1}。

6c)将获得的最优基矩阵W_{Opt_1}分别与无标签样本集X_{un_0}，训练样本集X_{tr_0}，测试样本X_{te_0}相乘，得到无标签样本集，训练样本集，测试样本的第一层深度稀疏特征表示X_{un_1}，X_{tr_1}，X_{te_1}，即：

无标签样本集的第一层深度稀疏特征表示：X_{un_1}＝W_{Opt_1}*X_{un_0}，

训练样本集的第一层深度稀疏特征表示：X_{tr_1}＝W_{Opt_1}*X_{tr_0}，

测试样本集的第一层深度稀疏特征表示：X_{te_1}＝W_{Opt_1}*X_{te_0}。

6d)将由第一层稀疏ICA提取到的无标签样本集的第一层深度稀疏特征表示X_{un_1}作为第二层稀疏ICA的输入，重复步骤6a)～6c)获得第二层稀疏ICA的最优基矩阵W_{Opt_2}，进而得到无标签样本集，训练样本集，测试样本的第二层深度稀疏特征表示X_{un_2}，X_{tr_2}，X_{te_2}，即

无标签样本集的第二层深度稀疏特征表示：X_{un_2}＝W_{Opt_2}*X_{un_1}，

训练样本集的第二层深度稀疏特征表示：X_{tr_2}＝W_{Opt_2}*X_{tr_1}，

测试样本集的第二层深度稀疏特征表示：X_{te_2}＝W_{Opt_2}*X_{te_1}，

其中，稀疏ICA的层数可以据分类结果进行调整，每加一层，都要重复上述步骤6a)～6c)一次，得到训练样本集和测试样本集最终的深度稀疏特征表示分别是X_tr和X_te。

步骤7，将训练样本集的最终深度稀疏特征表示X_tr输入到支持向量机分类器中进行训练，得到训练好的分类器。

步骤8，利用训练好的分类器，对待分类的极化SAR图像的每个像素进行分类，得到极化SAR图像的每个像素所属地物类别。

步骤9，输出结果。

对步骤(8)中的分类结果，将红色，绿色，蓝色三个颜色作为三基色，根据三基色上色法，用同一种颜色标出同一类地物，得到上色后的极化SAR图像T3，输出上色后的极化SAR图像T3。

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

1、仿真条件

硬件平台为：Intel(R)Xeon(R)CPU E56062.13GHZ、8GB RAM；

软件平台为：MATLAB R2013a；

仿真所使用的输入图像如图2所示，图像尺寸大小为750×1024。该图为NASA—JPL记载L波段AIRSAR系统于1991年获取的荷兰Flevoland农田区域极化合成孔径雷达图像，该图像分辨率是12.1m×6.7m。该极化合成孔径雷达SAR图像中区域的地表覆盖物为各种农作物和水域，按15类进行分类。

2、仿真内容与结果

仿真1，用现有的Wishart有监督分类方法对图2进行分类，结果如图4。

仿真2，用本发明提取的原始特征直接输入SVM对图2进行分类，结果如图5。

仿真3，用本发明提取的深度稀疏特征表示输入SVM对图2进行分类，结果如图6。

从图4、图5和图6的对比可以看出，图6所示的本发明提取的深度稀疏特征表示输入SVM分类得到的分类结果中，各区域均划分的较为细致，不同区域之间分类后的边缘也较为平滑，水域和农田的交接处被清晰的划分出来，边缘完整且准确。

以图3中的待分类极化SAR真实地物标记图像作为精度评价标准，对Wishart有监督分类方法、本发明提取的原始特征直接输入SVM分类方法及本发明提取的深度稀疏特征表示输入SVM的分类方法的分类精度进行统计，结果见表1。其中，Alg1表示Wishart有监督分类方法，Alg2表示直接将本发明提取的原始特征输入SVM分类的方法，Alg3表示将本发明提取的深度稀疏特征表示输入SVM分类的方法。

表1.三种方法在仿真中得到的分类精度

从表1中可以看出，用本发明提取的深度稀疏特征表示输入SVM分类相比于另外两种对比实验，精度有较大的提高，这主要是因为本发明提取的深度稀疏特征表示相比与本发明提取的原始特征和极化SAR图像的相干矩阵本身，包含了更加丰富的地物信息，且特征冗余性小，有利于分类，分类精度更高。

Claims (6)

1.一种基于深度稀疏ICA的极化SAR图像分类方法，包括如下步骤：

(1)输入待分类的极化合成孔径雷达SAR图像T1和待分类SAR图像的真实地物标记图像T2；

(2)采用精致Lee滤波方法，对待分类的极化SAR图像T1进行滤波，去除斑点噪声，得到滤波后的极化SAR图像T1’；

(3)提取滤波后的极化SAR图像T1’的散射、纹理、颜色及协方差矩阵特征；

(4)将散射特征、纹理特征、颜色特征及协方差矩阵特征组合成原始样本特征集，并将该原始样本特征集的值归一化到0～1之间，进行白化处理；

(5)选取无标签样本集，训练样本集和测试样本集；

5a)在白化处理过的原始样本特征集中随机选取2％的样本特征作为无标签样本集X_{un_0}；

5b)根据待分类SAR图像的真实地物标记图像T2，在白化处理过的原始样本特征集中选出与该T2中已标记真实地物相同的样本特征集，并在该样本特征集中随机选取20％的样本特征作为训练样本集X_{tr_0}，剩余的80％的样本特征作为测试样本集X_{te_0}；

(6)用深度稀疏ICA方法，分别提取无标签样本集、训练样本集和测试样本集的深度稀疏特征表示；

6a)用均值为0，方差为1的正态分布初始化第一层稀疏ICA的基矩阵W₁；

6b)利用基矩阵W₁，计算整体无标签样本集的方差衰减值，获取第一层稀疏ICA的最优基矩阵W_{Opt_1}；

6c)将获得的最优基矩阵W_{Opt_1}分别与无标签样本集X_{un_0}，训练样本集X_{tr_0}，测试样本集X_{te_0}相乘，得到无标签样本集，训练样本集，测试样本集的第一层深度稀疏特征表示X_{un_1}，X_{tr_1}，X_{te_1}；

6d)将由第一层稀疏ICA提取到的无标签样本集的第一层深度稀疏特征表示X_{un_1}作为第二层稀疏ICA的输入，重复步骤6a)～6c)获得第二层稀疏ICA的最优基矩阵W_{Opt_2}，进而得到无标签样本集，训练样本集，测试样本的第二层深度稀疏特征表示X_{un_2}，X_{tr_2}，X_{te_2}，

其中，稀疏ICA的层数可以据分类结果进行调整，每加一层，都要重复步骤6a)～6c)一次，得到训练样本集和测试样本集最终的深度稀疏特征表示分别是X_tr和X_te；

(7)将训练样本集的最终深度稀疏特征表示X_tr输入到支持向量机分类器中进行训练，得到训练好的分类器；

(8)利用训练好的分类器，对待分类的极化SAR图像T1的每个像素进行分类，得到极化SAR图像T1的每个像素所属地物类别；

(9)对步骤(8)中的分类结果，根据红、绿、蓝三基色原理，用同一种颜色标出同类地物，得到上色后的极化SAR图像T3，输出上色后的极化SAR图像T3。

2.根据权利要求1所述的方法，其中所述步骤(3)中提取滤波后的极化SAR图像的散射特征，是对滤波后的极化SAR图像的每个像素，采用Yamaguchi分解方法得到p_s，p_v，p_d，p_h，f_s，f_v，f_d，f_h共8个散射参数，采用Cloude分解方法得到A，H，α，λ₁，λ₂，λ₃共6个散射参数，采用Huynen分解方法得到|2A₀|，|C-iD|，|H+iG|共3个散射参数，将得到的这些散射参数作为17维的散射特征。

3.根据权利要求1所述的方法，其中所述步骤(3)中提取滤波后的极化SAR图像的纹理特征，按如下步骤进行：

首先，由滤波后的极化SAR图像每个像素点的协方差矩阵，得到该像素点的水平向接收水平向发射的回波数据|S_HH|²、垂直向接收水平向发射的回波数据|S_HV|²、垂直向接收垂直向发射的回波数据|S_VV|²，由每个像素点的这三个极化值|S_HH|²，|S_HV|²，|S_VV|²分别形成与滤波后的极化SAR图像大小相同的三幅极化通道图像；

其次，计算每幅极化通道图像的灰度共生矩阵，根据灰度共生矩阵提取每幅极化通道图像每个像素的能量、逆差距和对比度三个纹理特征，

然后，将三幅极化通道图像得到的9个纹理特征，作为滤波后的极化SAR图像的9维的纹理特征。

4.根据权利要求1所述的方法，其中所述步骤(3)中提取滤波后的极化SAR图像的颜色特征，是利用滤波后的极化SAR图像及其真实地物标记图像进行：

首先，从每类地物中随机选择三个像素，并求这三个像素的平均R，G，B颜色值作为该类地物颜色的R，G，B分量，得到K个RGB颜色值，其中，K为待分类极化SAR图像的地物类别数目；

其次，令这K个RGB颜色值为聚类中心，对于滤波后的SAR图像上的每个像素，计算其9*9的矩形窗口邻域中每个像素的RGB颜色值I＝(r,g,b)到各个聚类中心的欧式距离；根据最短距离原则，分别把窗口内的每个像素归到最小距离对应的地物类别中，再计算窗口内被归到各个地物类别的像素个数占窗口内总像素数目的比值，用得到的K个比值作为该窗口中心像素的K维颜色特征。

5.根据权利要求1所述的方法，其中所述步骤(3)中提取滤波后的极化SAR图像的协方差矩阵的特征，是从滤波后的极化SAR图像的协方差矩阵C中提取的HH-VV相关系数的幅度、HH-VV相位差、HH/VV通道比、同极化比、交叉极化、同极化系数比、去极化系数比、C(1,2)的角度、C(1,3)的角度、C(2,3)的角度、C(1,1)的模值、C(1,2)的模值、C(1,3)的模值、C(2,2)的模值、C(2,3)的模值、C(3,3)的模值，得到的16维的协方差矩阵特征，其中，H表示水平极化方向，V表示垂直极化方向，C(a,b)表示协方差矩阵C第a行第b列元素值。

6.根据权利要求1所述的方法，其中所述步骤6b)中利用基矩阵W₁，计算整体无标签样本集的方差衰减值，获取第一层稀疏ICA的最优基矩阵W_{Opt_1}，按如下步骤进行：

6b1)根据第一层稀疏ICA，计算第一层深度特征表示的稀疏度：

$P(W_1,X_{\mathrm{un}\_0})=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}\underset{j=1}{\overset{t}{\mathrm{\&Sigma;}}}g\left(W_{1\_j}{x}^{\left(i\right)}\right)$

其中，W₁表示选取的第一层稀疏ICA的基矩阵，其初始值服从均值为0，方差为1的正态分布，W_{1_j}表示基矩阵W₁中的第j个基向量，x⁽ⁱ⁾表示第i个无标签样本，g表示非线性凸函数，t表示基矩阵W₁中基向量的数目，m表示无标签样本集中样本数。

6b2)计算无标签样本集与第一层稀疏ICA的输出W₁ ^TW₁xⁱ之间的误差：

$J_1(W_1,X_{\mathrm{un}\_0})=\frac{\mathrm{\&lambda;}}{m}\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\left|\right|{W_1}^T{W}_1{x}^{\left(i\right)}-x^{\left(i\right)}\left|\right|}_2^2,$

其中，λ为平衡参数，用于对误差J₁(W₁,X_{un_0})和稀疏度P(W₁,X_{un_0})进行平衡，W₁x⁽ⁱ⁾为第一层深度稀疏特征表示，W₁ ^TW₁xⁱ为第一层深度稀疏特征表示W₁x⁽ⁱ⁾对无标签样本x⁽ⁱ⁾进行恢复的结果；

6b3)根据步骤6b1)和6b2)的计算结果，计算整体无标签样本集的均方差衰减值：

J(W₁,X_{un_0})＝J₁(W₁,X_{un_0})+P(W₁,X_{un_0})，

6b4)对步骤6b3)得到的整体无标签样本集的均方差衰减值，采用如下的梯度下降公式，获取每一次迭代时的第一层稀疏ICA的基矩阵：

${W_1}^{n+1}={W_1}^n-\mathrm{\&alpha;}\frac{\&PartialD;}{\&PartialD;{W_1}^n}J({W_1}^n,X_{\mathrm{un}\_0}),$

其中，n表示迭代的次数，W₁ ⁿ表示第n次迭代时第一层稀疏ICA的基矩阵，W₁ ⁿ⁺¹表示第n+1次迭代时第一层稀疏ICA的基矩阵，β表示对第一层稀疏ICA基矩阵进行学习的速率，β的取值范围为0<β<1，表示对第n次迭代时第一层稀疏ICA的基矩阵求偏导数。

6b5)将步骤6b4)每一次迭代得到的第一层稀疏ICA的基矩阵值代入步骤6b3)公式中求整体无标签样本的均方差衰减值J(W₁,X_{un_0})，当整体无标签样本的均方差衰减值J(W₁,X_{un_0})达到全局最小值时，停止步骤6b4)的迭代，并将停止迭代时第一层稀疏ICA的基矩阵作为第一层稀疏ICA的最优基矩阵W_{Opt_1}。