CN103955697A - 基于2dpca的极化sar图像分类方法 - Google Patents

Abstract

本发明公开了一种基于2DPCA的极化SAR图像分类方法，主要解决现有无监督极化SAR分类方法分类精度不高的问题。其实现步骤为：对每个像素点进行Freeman分解，提取像素点的三种散射功率；根据获得的散射功率对图像进行划分，得到3种类别；对获得的每一类，将其用2DKPCA进行自适应降维分类；最后，对预分类得到的图像用Wishart分类器进行迭代分类，得到最终分类结果。本发明与经典分类方法相比，对极化SAR图像的划分更加严谨，分类效果更佳，计算复杂度相对较小，可用于对极化SAR图像的进行地物分类和目标识别。

Description

基于2DPCA的极化SAR图像分类方法

技术领域

本发明属于图像处理技术领域，涉及极化SAR图像地物分类领域的应用，具体的说是一种基于2DPCA的极化SAR图像分类方法，可用于对极化SAR图像的地物分类和目标识别。

背景技术

极化SAR雷达能够得到更丰富的目标信息，在农业、林业、军事、地质学、水文学和海洋等方面具有广泛的研究和应用价值，如地物种类的识别、农作物成长监视、产量评估、地物分类、海冰监测、地面沉降监测，目标检测和海洋污染检测等。极化SAR图像分类的目的是利用机载或星载极化传感器获得的极化测量数据，确定每个像素所属的类别。经典的极化SAR分类方法包括：

Cloude等人提出了基于H/α目标分解的极化SAR图像非监督分类方法，见Cloude S R，Pottier E.An entropy based classification scheme for land applications ofpolarimetric SAR[J].IEEE Trans.Geosci.Remote Sensing.1997，35(1)：549-557.该方法主要是通过Cloude分解提取两个表征极化数据的特征H和α，然后根据H/α平面将极化图像划分为9个区域，其中一个理论上不可能存在的区域，因此最终将图像划分为8类。H/α分类存在的两个缺陷：一个是区域的划分过于武断；另一个是当同一个区域里共存几种不同的地物时，将不能有效区分。

Lee等人提出了基于H/α目标分解和Wishart分类器的H/α-Wishart非监督分类方法，见Lee J S，Grunes M R，Ainsworth T L，et a1.Unsupervised classificationusing polarimetric decomposition and the complex Wishart classifier[J].IEEETrans.Geosci.Remote Sensing.1999，37(5)：2249-2258.该方法是在原始H/α分类基础上增加了Wishart迭代，主要是对H/α划分后的8类利用Wishart分类器对每一个像素进行重新划分，从而有效的提高了分类的精度，但是又存在不能很好的保持各类的极化散射特性的不足。

2DPCA是近年来用于人脸识别种特征提取的，它是在PCA的基础上发展起来的，与PCA方法相比，它有两个优点：一是更易于准确估算方差矩阵；二是计算本征矢量时会更快。

发明内容

本发明的目的在于针对已有技术的不足，提出一种基于2DPCA的极化SAR图像分类方法，以提高分类效果。

为实现上述目的，本发明包括如下步骤：

步骤1，对待分类的极化SAR图像中每个像素点的协方差矩阵C进行Freeman分解，得到每个像素点的体散射功率P_v、二面角散射功率P_d和表面散射功率P_s；

步骤2，根据每个像素点的体散射功率P_v、二面角散射功率P_d和表面散射功率P_s，对极化SAR图像进行初始划分;

步骤3，对待分类的极化SAR图像中每个像素点的协方差矩阵C变换，得到变换矩阵变换的规则为：

其中，是变换后的矩阵，参数代表方向角，χ代表椭圆率，T(ρ)是单位极化变换矩阵：

$T\left(\mathrm{\ρ}\right)=\frac{1}{1+{\mathrm{\ρ\ρ}}^{*}}\left[\begin{array}{ccc}1& \sqrt{2}\mathrm{\ρ}& {\mathrm{\ρ}}^2\\ -\sqrt{2}{\mathrm{\ρ}}^{*}& 1-{\mathrm{\ρ\ρ}}^{*}& \sqrt{2}\mathrm{\ρ}\\ {{\mathrm{\ρ}}^{*}}^2& -\sqrt{2}{\mathrm{\ρ}}^{*}& 1\end{array}\right]$

ρ是极化率：

步骤4，在初始划分的基础上用2DPCA对每个像素点变换后的矩阵进行降维，并自适应获得类别数；

步骤5，对初始分类得到的整幅图像用Wishart分类器进行迭代分类，得到最终的分类结果。

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

1.本发明用2DPCA来进行极化SAR的降维分类更能体现极化SAR数据的分布特性；

2.本发明可以不用对极化SAR图像进行滤波处理；

3.本发明利用三种散射功率将极化SAR图像大致划分为3大类别，然后在各大类别内进行2DPCA降维分类，可以很好的保持各类地物的极化散射特性；

4.本发明可根据图像的具体情况对图像自适应的进行聚类。

附图说明

图1是本发明的总流程图；

图2是本发明仿真使用的旧金山极化SAR数据的PauliRGB合成图；

图3是用现有H/α分类方法对图2的分类结果图；

图4是用现有基于Freeman分解的分类方法对图2的分类结果；

图5是用本发明对图2的分类结果。

具体实施方式

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

步骤1，对极化SAR图像中每个像素点的协方差矩阵C进行三种分量分解，得到每个像素点的体散射功率P_v、二面角散射功率P_d和表面散射功率P_s。

(1a)极化SAR图像的像素点用一个3×3的相干矩阵T表示，根据相干矩阵T得到协方差矩阵C；

$C=U^{-1}\mathrm{TU}=\left[\begin{array}{ccc}<{\left|S_{\mathrm{HH}}\right|}^2>& \sqrt{2}<S_{\mathrm{HH}}{S}_{\mathrm{HP}}^{*}>& <S_{\mathrm{HH}}{S}_{\mathrm{PP}}^{*}>\\ \sqrt{2}<S_{\mathrm{HP}}{S}_{\mathrm{HH}}^{*}>& 2<{\left|S_{\mathrm{HP}}\right|}^2>& \sqrt{2}<S_{\mathrm{HP}}{S}_{\mathrm{PP}}^{*}>\\ <S_{\mathrm{PP}}{S}_{\mathrm{HH}}^{*}>& \sqrt{2}<S_{\mathrm{PP}}{S}_{\mathrm{HP}}^{*}>& <{\left|S_{\mathrm{PP}}\right|}^2>\end{array}\right]---1)$

其中，U是中间变量， $U=\frac{1}{\sqrt{2}}\&times;\left(\begin{array}{ccc}1& 0& 1\\ 1& 0& -1\\ 0& \sqrt{2}& 0\end{array}\right),$ U^-1是U矩阵的转置矩阵，H表示水平极化，P表示垂直极化，S_HH表示水平向发射和水平向接收的回波数据，S_PP表示垂直向发射和垂直向接收的回波数据，S_HP表示水平向发射垂直向接收的回波数据，|·|表示系数的模，()^*表示这个数据的共轭，<·>表示按视数平均；

(1b)将协方差矩阵C分解成如下表示：

$C=F_s\left[\begin{array}{ccc}{\left|\mathrm{\&beta;}\right|}^2& 0& \mathrm{\&beta;}\\ 0& 0& 0\\ {\mathrm{\&beta;}}^{*}& 0& 1\end{array}\right]+F_d\left[\begin{array}{ccc}{\left|\mathrm{\&alpha;}\right|}^2& 0& \mathrm{\&alpha;}\\ 0& 0& 0\\ {\mathrm{\&alpha;}}^{*}& 0& 1\end{array}\right]+F_v\left[\begin{array}{ccc}1& 0& 1/3\\ 0& 2/3& 0\\ 1/3& 0& 1\end{array}\right]+F_h\left[\begin{array}{ccc}1& \&PlusMinus;j\sqrt{2}& -1\\ \stackrel{\&OverBar;}{+}j\sqrt{2}& 2& \&PlusMinus;j\sqrt{2}\\ -1& \stackrel{\&OverBar;}{+}j\sqrt{2}& 1\end{array}\right]---2)$

其中，F_s为表面散射分量的分解系数，F_d为二面角散射分量的分解系数，F_v为体散射分量的分解系数，F_h为螺旋散射分量的分解系数，β是水平发射水平接收后向散射反射系数与垂直发射垂直接收后向散射发射系数的比值，α被定义为α＝I_HI_PH/I_PI_PP，I_H和I_P分别表示地表的水平及垂直反射系数，I_PH和I_PP分别表示竖直墙体的水平及垂直反射系数，j为虚部；

(1c)将式2）中矩阵的元素与式1）中协方差矩阵C的元素相对应，获得一个具有六个未知数F_s,F_v,F_d,F_h,α,β和五个方程的方程组如下：

$\left\{\begin{array}{c}<{\left|S_{\mathrm{HH}}\right|}^2>=F_s{\left|\mathrm{\&beta;}\right|}^2+F_d{\left|\mathrm{\&alpha;}\right|}^2+\frac{1}{15}{F}_v+\frac{F_h}{4}\\ <{\left|S_{\mathrm{HP}}\right|}^2>=\frac{2}{15}{F}_v+\frac{F_h}{4}\\ <{\left|S_{\mathrm{PP}}\right|}^2>=F_s+F_d+\frac{3}{15}{F}_v+\frac{F_h}{4}\\ <S_{\mathrm{HH}}{S}_{\mathrm{PP}}^{*}>=F_s\mathrm{\&beta;}+F_d\mathrm{\&alpha;}+\frac{2}{15}{F}_v-\frac{F_h}{4}\\ \frac{1}{2}\mathrm{Im}(<S_{\mathrm{HH}}{S}_{\mathrm{HP}}^{*}>+<S_{\mathrm{HV}}{S}_{\mathrm{PP}}^{*}>)=\frac{F_h}{4}\end{array}\right.---3);$

其中，Im(·)表示取虚部；

(1d)求解方程组3)，得到F_s,F_v,F_d,F_h,α,β的值：

取像素点的协方差矩阵C中的的实部的值并与零相比较，如果则令β＝1，反之，则令α＝-1，根据得到α或β的值，解方程组3），得到其余的未知数F_s,F_v,F_d,F_h的值，其中，Re(·)表示取实部；

(1e)根据求解出的F_s,F_v,F_d,F_h,α,β，得到体散射功率P_v、二面角散射功率P_d、表面散射功率P_s和螺旋散射功率P_h：

P_s＝F_s(1+|β|²)

P_d＝F_d(1+|α|²)

P_v＝F_v   。4)

步骤2，根据每个像素点的体散射功率P_v、二面角散射功率P_d和表面散射功率P_s，对极化SAR图像进行初始划分。

根据max(P_sP_dP_v)的值，将极化SAR图像数据初始划分为三类，即将max(P_s,P_d,P_v)＝P_s的对应像素点划分为平面散射类，将max(P_s,P_d,P_v)＝P_d的对应像素点划分为二面角散射类，将当max(P_s,P_d,P_v)＝P_v对应像素点划分为体散射类。

步骤3，对极化SAR图像中每个像素点的协方差矩阵C变换，得到变换矩阵

3a)变换的规则为：

其中是变换后的矩阵，参数代表方向角，χ代表椭圆率。T(ρ)是单位极化变换矩阵：

$T\left(\mathrm{\&rho;}\right)=\frac{1}{1+{\mathrm{\&rho;\&rho;}}^{*}}\left[\begin{array}{ccc}1& \sqrt{2}\mathrm{\&rho;}& {\mathrm{\&rho;}}^2\\ -\sqrt{2}{\mathrm{\&rho;}}^{*}& 1-{\mathrm{\&rho;\&rho;}}^{*}& \sqrt{2}\mathrm{\&rho;}\\ {{\mathrm{\&rho;}}^{*}}^2& -\sqrt{2}{\mathrm{\&rho;}}^{*}& 1\end{array}\right],$

ρ是极化率：

3b)考虑只有线极化的情况，即χ＝0，此时C矩阵随着方向角的变化而变化，取得到每个像素点的变换后矩阵矩阵是3*3维的，i，j表示矩阵的坐标。

步骤4，在初始划分的基础上用2DPCA对每个像素点变换后的矩阵进行降维，并自适应获得类别数。

4a)对每个像素点的算其协方差矩阵G_ij；

4b)接着计算G_ij的特征向量和特征值，然后把特征向量按照相对应的特征值大小从大到小排序，定义前d个特征向量x₁,x₂,...x_d为每个像素点的主元。

d的取值可以采用如下的方式确定：先将G_ij的特征值从大到小排列；对G_ij的所有特征值求和；从大到小叠加上述特征值，若叠加到第d个时，叠加的特征值的和>所有特征值的和的D%，就取d个。优选地，D=90。

4c)对每个像素点设定统一的阀值，出于每个像素点包含的极化散射特征不一样，满足统一的阀值的主元个数也就会有所区别。

4d)最后以确定的主元数目来自适应的确定分类类别。

步骤5，对预分类得到的整幅图像用Wishart分类器进行迭代分类，得到最终的分类结果。

(5a)对整幅图像预分类得到的结果，计算每一类的聚类中心B_c：

$B_c=\frac{\underset{\mathrm{\&rho;}=1}{\overset{n_c}{\mathrm{\&Sigma;}}}{T}_{\mathrm{\&rho;}}}{n_c},c=1,...,k,\mathrm{\&rho;}=\mathrm{1,2},...,n_c,$

其中T_ρ是第c类中各像素点的相干矩阵，n_c是属于第c类的像素点的个数，k为预分类得到的类别数；

(5b)根据每一类的聚类中心B_c，计算每个像素点到第c类聚类中心的距离：

$d(<T>,B_c)=\ln \left[B_c\right]+\mathrm{Tr}(B_c^{-1}<T>),c=1,...k,$

其中T是像素点的协方差矩阵，ln[·]表示求对数，<·>表示按视数平均，[·]表示矩阵的行列式，Tr(·)表示矩阵的迹，表示对聚类中心B_c求逆；

(5c)根据每个像素点到各聚类中心的距离，对预分类后得到整幅图像的类别重新进行划分：

如果d(<T>,B_τ)≤d(<T>,B_ψ)，则将该像素点划分为第τ类，如果d(<T>,B_τ)＞d(<T>,B_ψ)，则将该像素点划分为第ψ类，

其中，d(<T>,B_τ)是像素点到第τ类聚类中心的距离，d(<T>,B_ψ)是像素点到第ψ类聚类中心的距离，τ,ψ＝1,...,k,τ≠ψ；

(5d)重复步骤(5a)-(5c)，直到迭代次数等于给定的迭代次数μ，得到最终分类结果。

优选地，μ＝3。

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

1、实验条件与方法

硬件平台为：Intel Core2Duo CPU E6550@2.33GHZ、2GB RAM；

软件平台为：MATLAB R2012b；

实验方法：分别为本发明和现有的H/α方法以及基于Freeman分解的方法，其中现有的这两种方法都是极化SAR图像分类中引用较多的经典方法。

2、仿真内容与结果

本发明将图2所示的旧金山极化SAR图像作为测试图像，大小为900×700。

仿真一，用现有的H/α分类方法对图2进行分类，分类结果见图3。

由图3可见，地物类别可以大致得到划分，但由于分类类别数的限制，还有较多区域划分不清楚，且对于区域边界的划分过于武断。

仿真二，用现有的基于Freeman分解的分类方法对图2进行分类，分类结果见图4。

由图4可见，分类效果与图3相比，类别区分更合理，并且分类类别数不受限制，但由于Freeman分解并不适用于对城区进行划分，因此分类结果中很多的区域边缘模糊不清晰。

仿真三，用本发明对图2进行分类，分类结果见图5。

由图5可见，本发明得到的分类结果与图3和图4相比，原本不清晰的山脉、植被和城区更易区分，地物细节信息体现得更加精细，其中高尔夫球场、跑马场、左上方的山峰等均可辨识出来。

综上所述，本发明的方法对极化SAR图像的无监督分类能取得更好的分类结果。

Claims (10)

1.一种基于2DPCA的极化SAR图像分类方法，包括如下步骤： 

步骤1，对待分类的极化SAR图像中每个像素点的协方差矩阵C进行Freeman分解，得到每个像素点的体散射功率P_v、二面角散射功率P_d和表面散射功率P_s； 

步骤2，根据每个像素点的体散射功率P_v、二面角散射功率P_d和表面散射功率P_s，对极化SAR图像进行初始划分； 

步骤3，对待分类的极化SAR图像中每个像素点的协方差矩阵C变换，得到变换矩阵变换的规则为： 

其中，是变换后的矩阵，参数代表方向角，χ代表椭圆率，T(ρ)是单位极化变换矩阵： 

ρ是极化率： 

步骤4，在初始划分的基础上用2DPCA对每个像素点变换后的矩阵进行降维，并自适应获得类别数； 

步骤5，对初始分类得到的整幅图像用Wishart分类器进行迭代分类，得到最终的分类结果。 

2.根据权利要求1所述的基于2DPCA的极化SAR图像分类方法，其中，对待分类的极化SAR图像中每个像素点的协方差矩阵C进行Freeman分解，得到每个像素点的体散射功率P_v、二面角散射功率P_d和表面散射功率P_s包括： 

(1a)极化SAR图像的像素点用一个3×3的相干矩阵T表示，根据相干矩阵T得 到协方差矩阵C； 

其中，U是中间变量，U^-1是U矩阵的转置矩阵，H表示水平极化，P表示垂直极化，S_HH表示水平向发射和水平向接收的回波数据，S_PP表示垂直向发射和垂直向接收的回波数据，S_HP表示水平向发射垂直向接收的回波数据，|·|表示系数的模，()^*表示这个数据的共轭，<·>表示按视数平均； 

(1b)将协方差矩阵C分解成如下表示： 

其中，F_s为表面散射分量的分解系数，F_d为二面角散射分量的分解系数，F_v为体散射分量的分解系数，F_h为螺旋散射分量的分解系数，β是水平发射水平接收后向散射反射系数与垂直发射垂直接收后向散射发射系数的比值，α被定义为α＝I_HI_PH/I_PI_PP，I_H和I_P分别表示地表的水平及垂直反射系数，I_PH和I_PP分别表示竖直墙体的水平及垂直反射系数，j为虚部； 

(1c)将式2)中矩阵的元素与式1)中协方差矩阵C的元素相对应，获得一个具有六个未知数F_s,F_v,F_d,F_h,α,β和五个方程的方程组如下： 

其中，Im(·)表示取虚部； 

(1d)求解方程组3)，得到F_s,F_v,F_d,F_h,α,β的值： 

取像素点的协方差矩阵C中的的实部的值并与零相比较，如果 则令β＝1，反之，则令α＝-1，根据得到α或β的值，解方程组3)，得到其余的未知数F_s,F_v,F_d,F_h的值，其中，Re(·)表示取实部； 

(1e)根据求解出的F_s,F_v,F_d,F_h,α,β，得到体散射功率P_v、二面角散射功率P_d、表面散射功率P_s和螺旋散射功率P_h： 

P_s＝F_s(1+|β|²) 

P_d＝F_d(1+|α|²) 

P_v＝F_v                4)。 

3.根据权利要求1所述的基于2DPCA的极化SAR图像分类方法，其中，根据每个像素点的体散射功率P_v、二面角散射功率P_d和表面散射功率P_s，对极化SAR图像进行初始划分包括：根据max(P_v,P_d,P_s)的值，将极化SAR图像数据初始划分为三类，即平面散射类、二面角散射类和体散射类。 

4.根据权利要求1所述的基于2DPCA的极化SAR图像分类方法，其中，根据max(P_v,P_d,P_s)的值，将极化SAR图像数据初始划分为三类包括： 

将max(P_s,P_d,P_v)＝P_s的对应像素点划分为平面散射类，将max(P_s,P_d,P_v)＝P_d的对应像素点划分为二面角散射类，将当max(P_s,P_d,P_v)＝P_v对应像素点划分为体散射类。 

5.根据权利要求1所述的基于2DPCA的极化SAR图像分类方法，其中，在只有线极化的情况下，即χ＝0，此时C矩阵随着方向角的变化而变化，取 得到每个像素点的变换后矩阵i，j表示矩阵的坐标。 

6.根据权利要求1所述的基于2DPCA的极化SAR图像分类方法，其中，在初始 划分的基础上用2DPCA对每个像素点变换后的矩阵进行降维，并自适应获得类别数具体包括： 

4a)对每个像素点的算其协方差矩阵Gi_j； 

4b)接着计算G_ij的特征向量和特征值，然后把特征向量按照相对应的特征值大小从大到小排序，定义前d个特征向量x₁,x₂,...x_d为每个像素点的主元； 

4c)对每个像素点设定统一的阀值； 

4d)根据阈值确定主元数目； 

4e)以确定的主元数目来自适应的确定分类类别。 

7.根据权利要求6所述的基于2DPCA的极化SAR图像分类方法，其中，d的取值可以采用如下的方式确定： 

将G_ij的特征值从大到小排列； 

对G_ij的所有特征值求和； 

从大到小叠加上述特征值，若叠加到第d个时，叠加的特征值的和大于所有特征值的和的D％，就取d个。 

8.根据权利要求7所述的基于2DPCA的极化SAR图像分类方法，其中，D＝90。 

9.根据权利要求1所述的基于2DPCA的极化SAR图像分类方法，其中，对初始分类得到的整幅图像用Wishart分类器进行迭代分类，得到最终的分类结果包括： 

(5a)对整幅图像预分类得到的结果，计算每一类的聚类中心B_c： 

其中T_ρ是第c类中各像素点的相干矩阵，n_c是属于第c类的像素点的个数，k为预分类得到的类别数； 

(5b)根据每一类的聚类中心B_c，计算每个像素点到第c类聚类中心的距离： 

其中T是像素点的协方差矩阵，ln[·]表示求对数，<·>表示按视数平均，[·]表示 矩阵的行列式，Tr(·)表示矩阵的迹，表示对聚类中心B_c求逆； 

(5c)根据每个像素点到各聚类中心的距离，对预分类后得到整幅图像的类别重新进行划分： 

如果d(<T>,B_τ)≤d(<T>,B_ψ)，则将该像素点划分为第τ类，如果d(<T>,B_τ)＞d(<T>,B_ψ)，则将该像素点划分为第ψ类， 

其中，d(<T>,B_τ)是像素点到第τ类聚类中心的距离，d(<T>,B_ψ)是像素点到第ψ类聚类中心的距离，τ,ψ＝1,...,k,τ≠ψ； 

(5d)重复步骤(5a)-(5c)，直到迭代次数等于给定的迭代次数μ，得到最终分类结果。 

10.根据权利要求1所述的基于2DPCA的极化SAR图像分类方法，其中，μ＝3。