说明书基于协同表示和深度学习的极化SAR图像分类方法
技术领域
本发明属于图像处理技术领域，特别涉及极化SAR图像分类方法，可用于地物识别。
背景技术
雷达是一种可以实现全天候工作的主动探测系统，它可以穿透一定的地表，并且可以改变发射波的频率、强度。合成孔径雷达SAR是成像雷达技术的一种，它是利用雷达与目标的相对运动把尺寸较小的真实无线孔径用数据处理的方法合成一个较大的等效天线孔径雷达，具有全天候、全天时、高分辨的优势。而极化SAR是用来测量回波信号的新型雷达，它可以记录不同极化状态组合回波的相位差信息，能对目标进行全极化测量成像，大大提高了对地物的识别能力。极化SAR图像分类是极化SAR图像解译的重要步骤，是边缘提取、目标检测和识别的基础，可广泛应用于军事侦察、地形测绘、农作物生长监测等领域。
目前经典的极化SAR图像分类方法有：
1997年，Cloude等人提出了一种基于H/α极化分解的分类方法，该方法通过Cloude分解得到特征参数散射熵H和散射角α，然后根据两个特征参数不同的值，将目标分成8类。该方法的缺陷是位于类别边界特征相似的像素点会以随机的方式分配给不同的类别并且这两个特征不足以表示所有的极化SAR信息。
1999年，Lee等人提出了一种基于H/α极化分解和复Wishart分类器的H/α-Wishart分类方法，该方法将H/α极化分解方法得到的结果作为复Wishart分类器的初始分类，对划分好的8个类别中的每一个像素进行重新划分，从而提高分类的精度。该方法的缺陷是将分类类别数固定为8类，不能适应不同类别数的地物分类。
2004年，J.S.Lee等人提出了一种基于Freeman-Durden分解的极化SAR图像分类方法，该方法首先通过Freeman分解获取表征散射体散射特征的三个特征：平面散射功率、二面角散射功率和体散射功率，然后根据这三个特征的大小对数据进行初始划分，然后利用Wishart分类器进行进一步精确划分。但是该方法由于Freeman分解中的多类的划分以及合并，计算复杂度较高。
发明内容
本发明的目的在于针对上述现有技术方法的不足，提出了一种基于协同表示和深度学习的极化SAR图像分类方法，以降低极化SAR图像分类的计算复杂度，提高分类精度。
为实现上述目的，本发明的技术方案包括如下步骤：
(1)将极化SAR图像中每个3*3大小像素点的极化相干矩阵T作为输入数据，计算每个3*3大小像素点的极化协方差矩阵C，这两个矩阵T和C中均包括9个元素；用T的对角线上的三个元素T11、T22、T33构成总功率特征参数：S＝T11+T22+T33；
(2)从每个像素点的极化相干矩阵T中通过克劳德Cloude分解方法分解出散射熵H和反熵A两个散射参数；从每个像素点的极化协方差矩阵C中通过弗里曼-德登Freeman-Durden分解方法分解出表面散射功率Ps、二面角散射功率Pd和体散射功率Pv三个功率参数；
(3)用所述参数H、A、Ps、Pd、Pv和极化相干矩阵T的9个元素、极化协方差矩阵C的9个元素、总功率特征参数S，这24个特征作为每个像素点的特征矩阵B；用所有像素点的特征矩阵组成整幅图像的特征矩阵F＝[B1,B2,...,Bk,...,BN],其中Bk表示第k个像素点的特征矩阵，k＝1,2,…,N，N为整幅图像的总像素点数；
(4)根据实际地物分布，从每类像素点所对应的特征矩阵中选取100个像素点的特征矩阵作为训练样本集Y，取整幅图像的特征矩阵F作为测试样本集；
(5)将训练样本集Y作为初始字典，利用K-SVD算法学习得到学习字典D；
(6)用步骤(5)得到的学习字典D协同表示训练样本集Y和测试样本集F，利用最小二乘法求解训练样本集Y的表示系数测试样本集F的表示系数
(7)将步骤(6)得到的训练样本集的表示系数输入到一个两层的稀疏自编码器中训练，得到第一层稀疏自编码器的权值W1和偏置b1，第二层稀疏自编码器的权值W2和偏置b2，然后固定两层稀疏自编码器的参数，将训练样本集的表示系数输入，得到输出值hy；
(8)将测试样本集的表示系数输入到步骤(7)固定的两层稀疏自编码器中，得到测试样本集的表示系数的输出值hf；
(9)将步骤(7)得到的输出值hy输入到libSVM分类器中进行训练，并将步骤(8)得到的输出值hf输入到已训练好的libSVM分类器中，得到最终的分类结果。
本发明与现有技术相比，具有以下优点：
1、本发明结合了协同表示技术，有效地降低了计算复杂度；
2、本发明利用稀疏自编码器对表示系数进行深度学习，得到极化SAR图像特征更本质的表示，提高了分类精度；
3、本发明结合了libSVM分类器，降低了分类消耗的时间，提高了分类精度；
仿真结果表明，本发明方法较经典的H/α极化分解的分类方法和H/α-Wishart分类方法能更有效的对极化SAR图像进行分类。
附图说明
图1是本发明的流程图；
图2是本发明仿真所用的两幅测试图像；
图3是本发明与现有两种方法对San Francisco数据的分类实验结果对比图；
图4为本发明与现有两种方法对Flevoland数据的分类实验结果对比图。
具体实施方式
参照图1，本发明的具体实现步骤如下：
步骤一，计算极化协方差矩阵C、总功率特征参数S。
(1a)输入极化SAR图像每个3*3大小像素点的极化相干矩阵T；
(1b)通过下式计算每个像素点的极化协方差矩阵C：C＝M*T*M’，
式中，M＝[1/sqrt(2)]*m，m＝[101；10-1；0sqrt(2)0]，sqrt(2)表示2的平方根，M’表示M的转置矩阵。
(1c)用T的对角线上的三个元素T11、T22、T33构成总功率特征参数：S＝T11+T22+T33。
步骤二，提取极化特征。
(2a)从每个像素点的极化相干矩阵T中通过克劳德Cloude分解方法分解出散射熵H和反熵A两个散射参数，其公式如下：
Σi=13-Pilog3Pi]]>
A=λ2-λ3λ2+λ3,]]>
式中，H表示散射熵，Pi表示极化相干矩阵T的第i个特征值与所有特征值总和的比值，A表示反熵，λ2表示极化相干矩阵T的第二个特征值，λ3表示极化相干矩阵T的第三个特征值；
(2b)将极化协方差矩阵C按如下公式分解：
C=fs|β|20β000β*01+fs|α|20α000α*01+fv101/302/301/301]]>《1》
式中，fs为平面散射分量的分解系数，fd为二面角散射分量的分解系数，fv为体散射分量的分解系数，β是水平发射水平接收后向散射反射系数与垂直发射垂直接收后向散射反射系数的比值，α＝RghRvhRgvRvv，Rgh和Rgv分别表示地表的水平及垂直反射系数，Rvh和Rvv表示垂直墙体的水平及垂直反射系数，*表示矩阵的共轭，|·|2表示绝对值的平方；
(2c)将由步骤(1b)计算得到的极化协方差矩阵C表示为：
C=〈|SHH|2〉2〈SHHSHH*〉〈SHHSVV*〉2〈SHVSHH*〉2〈|SHV|2〉2〈SHVSVV*〉〈SVVSHH*〉2〈SVVSHV*〉〈|SVV|2〉]]>《2》
式中，H表示水平极化，V表示垂直极化，SHH表示水平发射水平接收的回波数据，SVV表示垂直发射垂直接收的回波数据，SHV表示水平发射垂直接收的回波数据，<·>表示按视数平均；
(2d)将式《1》中矩阵的元素与式《2》中极化协方差矩阵C的元素相对应，获得一个具有五个未知数fs、fv、fd、α、β和四个方程的方程组如下：
〈|SHH|2〉=fs|β|2+fd|α|2+fv〈|SVV|2〉=fs+fd+fv〈|SHHSVV*|2〉=fsβ+fdβ+fv/3〈|SHV|2〉=fv/3]]>《3》
(2e)计算像素点协方差矩阵C中的的值并判断正负，如果 则α＝-1，如果则β＝1，给定α或β的值后，可根据式《3》求解得出5个未知数fs、fv、fd、α、β的值，其中Re(·)表示取实部；
(2f)根据求解出的fs、fv、fd、α、β，按照下式求解体散射功率Pv二面角散射功率Pd、表面散射功率Ps：
Pv=8fv/3Pd=fd(1+|α|2)Ps=fs(1+|β|2).]]>
步骤三，获取整幅图像的特征矩阵F。
(3a)用所述参数H、A、Ps、Pd、Pv和极化相干矩阵T的9个元素、极化协方差矩阵C的9个元素、总功率特征参数S，这24个特征作为每个像素点的特征矩阵B；
(3b)用所有像素点的特征矩阵组成整幅图像的特征矩阵F＝[B1,B2,...,Bk,...,BN],其中Bk表示第k个像素点的特征矩阵，k＝1,2,…,N，N为整幅图像的总像素点数。
步骤四，选取训练样本集和测试样本集。
(4a)根据实际地物分布，从每类像素点所对应的特征矩阵中选取100个像素点的特征矩阵作为训练样本集Y；
(4b)取整幅图像的特征矩阵F作为测试样本集。
步骤五，字典学习。
(5a)将训练样本集Y作为K-SVD算法中的初始字典；
(5b)通过K-SVD算法按照如下公式得到学习字典D：
min{||Y-DX||22Subjectto∀j,||Xj||0≤T0,]]>
式中min||·||表示让·的值达到最小，Subject to表示约束条件，X为系数矩阵，表示任意第j列，j＝1,2,…,K，K为系数矩阵X的总列数，||·||0表示向量的0范数，为矩阵的2范数的平方，T0为稀疏表示中稀疏向量中非零值的个数的上限。
步骤六，求解训练样本集的表示系数和测试样本集的表示系数
(6a)利用步骤五求出的学习字典D构建训练样本集Y和测试样本集F的协同表示模型：
α^y=argminαy{||Y-Dαy||22+λ||αy||22},]]>
α^f=argminαf{||F-Dαf||22+λ||αf||22},]]>
式中，argminαy{||Y-Dαy||22+λ||αy||22}]]>表示使目标函数取最小值时变量αy的值，λ表示正则化参数,为训练样本集Y的表示系数，为测试样本集F的表示系数；
(6b)利用最小二乘法求解步骤(6a)中构建的协同表示模型，得到训练样本集Y的表示系数测试样本集F的表示系数
α^y=(DTD+λI)-1·DT·Y,]]>
α^f=(DTD+λI)-1·DT·F,]]>
式中，DT表示学习字典D的转置，(·)-1表示矩阵·的逆，I表示单位矩阵。
步骤七，通过深度学习训练稀疏自动编码器。
(7a)随机初始化两层稀疏自编码器的权值W1、W2，初始化两层稀疏自编码器的偏置b1＝0，b2＝0；
(7b)将训练样本集的表示系数输入到两层自编码器中进行深度学习，得到两层稀疏自编码器训练后的参数，并对其进行固定；
(7c)将训练样本集的表示系数输入到已固定的稀疏自编码器中，得到输出值hy。
步骤八，将测试样本集的表示系数输入到步骤七固定的两层稀疏自编码器中，得到测试样本的表示系数的输出值hf。
步骤九，将步骤七得到的输出值hy输入到libSVM分类器中进行训练；将步骤八得到的输出值hf输入到已训练好的libSVM分类器中，得到最终的分类结果。
本发明的效果可通过以下仿真进一步说明。
1.实验条件和与方法：
实验仿真环境：MATLAB 2013a，Windows XP Professional。
实验方法：分别为H/α极化分解的分类方法和H/α-Wishart分类方法和本发明，其中前两种方法为极化SAR图像分类的经典方法。
2.实验内容与结果分析：
实验内容：本发明使用的是图2所示的两组极化SAR图像数据，图2(a)为美国San Francisco地区数据，视数为四，图2(b)为荷兰Flevoland地区的数据，视数为四，两组数据都来源于美国宇航局喷气推进实验室的AIRSAR传感器。
仿真1，用本发明以及H/α极化分解的分类方法和H/α-Wishart分类方法对美国San Francisco地区数据进行分类实验，结果如图3所示，其中：
图3(a)为H/α极化分解的分类方法的分类结果，分为9类；
图3(b)为H/α-Wishart分类方法的分类结果，分为9类；
图3(c)为用本发明方法的分类结果，分为3类。
从图3可以看出，H/α极化分解的分类方法分类结果很不理想，各个区域中均出现不同程度的混杂现象，H/α-Wishart分类方法的分类结果优于H/α极化分解的分类方法，区域划分的更加细致，但图像细节保持较差；而本发明的分类结果从视觉上看分类效果更好，其中跑马场、高尔夫球场等区域在分类后的图中，区域一致性好于H/α极化分解的分类方法和H/α-Wishart分类方法，左上角的陆地部分分类清晰。
仿真2，用本发明以及H/α极化分解的分类方法和H/α-Wishart分类方法对荷兰Flevoland地区数据进行分类实验，结果如图4所示，其中：
图4(a)为H/α极化分解的分类方法的分类结果，分为9类；
图4(b)为H/α-Wishart分类方法的分类结果，分为9类；
图4(c)为用本发明方法的分类结果，分为13类。
从图4可以看出，H/α-Wishart分类方法和H/α极化分解的分类方法由于固定了类别数目，对此图不能很精确的划分，很多类别都被归为了一类，而本发明清晰了地分出了各类别的轮廓，分类效果明显好于H/α-Wishart分类方法和H/α极化分解的分类方法，并且边缘清晰细节信息完整。