CN105138966B - 基于快速密度峰值聚类的极化sar图像分类方法 - Google Patents

Abstract

本发明公开了一种基于快速密度峰值聚类的极化SAR图像分类方法，主要解决现有无监督极化SAR分类方法分类精度不高的问题。其实现步骤为：1.对待分类极化SAR图像进行滤波；2.对滤波后的极化SAR图像进行Yamaguchi分解，计算每个像素点的四种散射功率；3.提取每个像素点主散射功率，将整幅图像初始划分为四类；4.按照初始类别中所有像素点的主散射功率，将整个极化SAR图像划分为M类；5.以M个类别的中心点作为新像素点进行聚类，并将聚类结果转化为整幅图像的预分类结果；6.对预分类结果进行迭代分类，得到最终分类结果。实验表明，本发明的分类效果更好，可用于对各种极化SAR图像进行无监督分类。

Description

基于快速密度峰值聚类的极化SAR图像分类方法

技术领域

本发明属于图像处理技术领域，涉及极化SAR图像分类方法，可用于目标识别。

背景技术

随着遥感技术在卫星、载人航天、探月工程等领域的发展，极化SAR已成为SAR的发展趋势，极化SAR能够得到更丰富的目标信息，在农业、林业、军事、地质学、水文学和海洋等方面具有广泛的研究和应用价值，如地形测绘、资源勘探、灾害监测及天文研究等不同领域。极化图像分类的目的是利用机载或星载极化传感器获得的极化测量数据，确定每个像素所属的类别。近年来，利用极化SAR测量数据进行的分类在国际遥感领域受到高度重视，已成为SAR图像分类的主要研究方向。

根据处理方法的不同，全极化SAR图像分类方法可以分为监督分类和非监督分类。监督分类方法主要包括基于统计知识、神经网络、小波分析和模糊逻辑的分类方法。非监督分类方法中最多的是基于目标散射机理的分解方法，如H/α目标分解方法和Freeman分解方法。相对于监督分类方法而言，基于散射机理的全极化SAR非监督分类方法不需要已知类别的训练数据，具有更强的适应性。

经典的极化SAR图像分类方法包括：

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

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

发明内容

本发明的目的在于针对已有技术的不足，提出一种基于快速密度峰值聚类的极化SAR图像分类方法，以提高分类正确率。

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

(1)对待分类极化SAR图像进行滤波，去除斑点噪声，得到滤波后的极化SAR图像；

(2)对滤波后的极化SAR图像中每个像素点的相干矩阵T进行Yamaguchi分解，得到每个像素点的体散射功率P_v、二面角散射功率P_d、表面散射功率P_s和螺旋体散射分量P_h；

(3)按照每个像素点四个散射功率P_s,P_d,P_v,P_h的最大值将极化SAR图像初始划分为四类：

如果max(P_s,P_d,P_v,P_h)＝P_s，则将其对应的像素点划分为一类，其中P_s为该类的主散射功率；

如果max(P_s,P_d,P_v,P_h)＝P_d，则将其对应的像素点划分为一类，其中P_d为该类的主散射功率；

如果max(P_s,P_d,P_v,P_h)＝P_v，则将其对应的像素点划分为一类，其中P_v为该类的主散射功率；

如果max(P_s,P_d,P_v,P_h)＝P_h，则将其对应的像素点划分为一类，其中P_h为该类的主散射功率；

其中，max(·)表示最大值；

(4)对每类初始类别所有像素点的主散射功率按由小到大排序，并将每类中每300个像素点再划分为一类，最终将整个极化SAR图像划分为M类；

(5)在已获得的M个类别中，将每一类的中心点作为新像素点，得到M个新像素点A_i，i＝1,…,M，并用每类的新像素点A_i代表该类别内所有像素点；

(6)对上述M个新像素点进行快速密度峰值聚类，将M个新像素点聚为k类；

(7)在M个新像素点的聚类结果中，将由新像素点A_i代表的所有像素点标记为与新像素点A_i相同的类别，完成对整幅图像的预分类；

(8)对整个极化SAR图像的预分类结果进行复Wishart迭代，得到更为准确的分类结果。

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

1.本发明充分利用了Yamaguchi分解得到的四种散射功率P_v,P_d,P_s,P_h在分类中的有效性，使一些不满足反射对称性的复杂区域能更够好的分类。

2.本发明针对现有极化SAR分类中很多经典分类方法都局限于特定的分类类别数的缺陷，采用快速密度峰值聚类，可根据不同极化SAR图像的具体情况进行自适应聚类，自适应的选择类别数目，且分类结果区域一致性划分较好，不同区域划分后的边缘也更加清晰。

附图说明

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

图2是本发明中快速密度峰值聚类的子流程图；

图3是本发明仿真使用的原始的旧金山极化SAR数据图像；

图4是用现有H/α方法、H/α-Wishart方法及本发明方法对图3的分类结果图；

图5是本发明仿真使用的原始的弗莱福兰农田极化SAR数据图像；

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

具体实施方式

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

步骤1，对待分类极化SAR图像进行滤波，去除斑点噪声，得到滤波后的极化SAR图像。

对极化SAR图像的滤波通常都是采用现有的精致极化LEE滤波法，滤波窗口的大小为7×7。

步骤2，对滤波后的极化SAR图像中每个像素点的相干矩阵T进行Yamaguchi分解，得到每个像素点的体散射功率P_v、二面角散射功率P_d、表面散射功率P_s和螺旋散射分量P_h。

Yamaguchi分解是一种极化目标分解方法，不同的地物目标具有不同的散射特性，该方法根据不同地物目标的散射特性，将每个像素点的相干矩阵分解成多种散射分量的线性组合，该方法相对其它目标分解方法的优势在于它增加了螺旋体散射分量，使一些不满足反射对称性的复杂区域能更够好的分类，见Yamaguchi Y,Maoriyama T,Ishido M,eta1.Four-component scattering model for polarimetric SARimage decomposition[J].IEEE Trans.Geosci.Remote Sensing.2005,43(8):1699-1706。Yamaguchi分解的具体步骤如下：

(2a)读入滤波后图像的每个像素点，这些像素点为一个3×3的相干矩阵T，根据相干矩阵T得到协方差矩阵C；

其中，U是中间变量，H表示水平极化，V表示垂直极化，S_HH表 示水平向发射和水平向接收的回波数据，S_VV表示垂直向发射和垂直向接收的回波数据，S_HV 表示水平向发射垂直向接收的回波数据，()^*表示这个数据的共轭，〈·〉表示按视数平均；

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

其中，f_s为平面散射分量的分解系数，f_d为二面角散射分量的分解系数，f_v为体散射分量的分解系数，f_h为螺旋散射分量的分解系数，β是水平发射水平接收后向散射反射系数与垂直发射垂直接收后向散射发射系数的比值，α被定义为α＝I_HI_VH/I_VI_VV，I_H和I_V分别表示地表的水平及垂直反射系数，I_VH和I_VV分别表示竖直墙体的水平及垂直反射系数，j表示虚数；

(2c)由协方差矩阵C分解的结果，获得一个具有六个未知数f_s,f_v,f_d,f_h,α,β和五个方程的方程组如下：

(2d)计算像素点协方差矩阵C中的的值，如果则令α＝-1，如果则令β＝1，给定α或β的值后，剩余的4个未知数则根据式3)求解，其中Re(·)表示取实部；

(2e)根据求解出的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+|a|²)

P_v＝f_v

P_h＝f_h。

步骤3，按照每个像素点四个散射功率P_s,P_d,P_v,P_h的最大值将极化SAR图像初始划分为四类：

如果max(P_s,P_d,P_v,P_h)＝P_s，则将其对应的像素点划分为一类，其中P_s为该类的主散射功率；

如果max(P_s,P_d,P_v,P_h)＝P_d，则将其对应的像素点划分为一类，其中P_d为该类的主散射功率；

如果max(P_s,P_d,P_v,P_h)＝P_v，则将其对应的像素点划分为一类，其中P_v为该类的主散射功率；

如果max(P_s,P_d,P_v,P_h)＝P_h，则将其对应的像素点划分为一类，其中P_h为该类的主散射功率；

其中，max(·)表示最大值。

步骤4，对每类初始类别所有像素点的主散射功率按由小到大排序，并将每类中每300个像素点再划分为一类，最终将整个极化SAR图像划分为M类。

步骤5，在已获得的M个类别中，将每一类的中心点作为新像素点，得到M个新像素点A_i，i＝1,…,M，并用每类的新像素点A_i代表该类别内所有像素点。

步骤6，对上述M个新像素点进行快速密度峰值聚类，将M个新像素点聚为k类。

常用的聚类方法主要包括K-means聚类、谱聚类、模糊聚类和快速密度峰值聚类。本实例采用快速密度峰值聚类，该聚类方法是机器学习领域一种新的聚类方法，该方法认为聚类中心由一些具有低局部密度点的数据包围，而这些低密度数据点距离其它高密度数据点距离较大，并以此为根据自适应选择聚类中心和聚类类别，见RodriguezA,andAlessandro L.Clustering by fast search and find of density peaks[J].Science27.2014,344(6191):1492-1496.

参照图2，本步骤的具体实现如下：

6a)计算任意两个新像素点A_i和A_j之间的相互距离d_ij:

d_ij＝Tr((T_i)^-1+(T_j)^-1T_i)-q，

其中，T_i，T_j分别表示新像素点A_i和A_j的相干矩阵，(T_i)^-1和(T_j)^-1分别表示对矩阵T_i和T_j求逆，q为常数，取值为q＝3，Tr(·)是矩阵的迹；

6b)计算新像素点A_i的局部密度ρ_i：

其中，当d_ij＜d_c时，χ(d_ij-d_c)＝1，否则χ(d_ij-d_c)＝0；M表示新像素点A_i的个数；d_c为常数，其取值为将所有点的相互距离d_ij由小到大排列，将位于该排列2％位置处的值作为d_c取值；

6c)计算新像素点A_i的距离δ_i：

6d)选择局部密度ρ_i和距离δ_i乘积最大的k个新像素点作为聚类中心；

6e)聚类中心确定后，用每个聚类中心代表一类，k个聚类中心共代表k类，比较剩余新像素点与每一类聚类中心的距离，如果该新像素点与第m类聚类中心的距离最近，则将该新像素点划分到第m类,m＝1,…,k。

步骤7，在M个新像素点的聚类结果中，将由新像素点A_i代表的所有像素点标记为与新像素点A_i相同的类别，完成对整幅图像的预分类；

步骤8，对预分类得到的整幅图像用能反映极化SAR分布特性的Wishart分类器进行迭代分类，得到更为准确的分类结果。

Wishart分类器是一种能反映极化SAR数据分布特性的分类器，通过比较像素点与各聚类中心的距离，判断像素点的所属类别，见Lee J S，Grunes M R，Pottier E，eta1.Unsupervised terrain classification preserving polarimetric scatteringcharacteristic[J].IEEE Trans.Geosci.Remote Sensing.2004，42(4)：722-731.其步骤如下：

(8a)对整个极化SAR图像预分类得到的k类划分结果，根据下式求每一类的聚类中心B_c：

其中T_ρ是第c类中各像素点的相干矩阵，n_c是属于第c类的像素点的个数；

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

其中T是像素点的相干矩阵，〈·〉表示按视数平均，[·]表示矩阵的行列式，Tr(·)表示矩阵的迹，表示对聚类中心B_c求逆；

(8c)比较每个像素点与各聚类中心的距离，如果该新像素点与第m类聚类中心的距离最近，则将该像素点划分到第m类,m＝1,…,k，完成对预分类后整幅极化SAR图像类别的重新划分；

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

本发明的效果可以通过以下实验验证:

1、实验条件与方法

硬件平台为：Intel(R)Pentium(R)1CPU2.4GHz；

软件平台为：WindowXPProfessional，MATLAB 7.0.4；

实验方法：分别为现有H/α的方法以及H/α-Wishart方法和本发明，其中现有的这两种方法都是极化SAR数据分类中引用较多的经典方法。

2.实验内容及结果

实验一，将图3所示视数为四的SanFranciscoBay极化SAR图像作为测试图像，用本发明方法与现有H/α方法以及H/α-Wishart方法对图3进行分类仿真，分类结果见图4，其中，图4(a)是H/α方法分类的结果，图4(b)是H/α-Wishart方法分类的结果，图4(c)为本发明方法的分类结果。

由图4(a)可见，图像中的水域部分得到了比较好的划分，但城区和绿地等混淆严重。因此，该方法分类规则过于武断，导致分类效果不佳。

由图4(b)可见，结合H/α和Wishart分类器的H/α-Wishart分类方法，图像区域划分的比较细致，但还有较多区域划分不清楚，甚至出现错误划分，如图4(b)右上角的海洋区域出现的错误划分。

由图4(c)可见，本发明的分类结果从视觉上看效果更好，其中高尔夫球场、跑马场、停车场等区域的区域一致性明显好于前两种方法，不同区域之间分类后的边缘也更加平滑。

实验二，将图5所示视数为四的Flevoland极化SAR图像作为测试图像。用本发明对弗莱福兰省农田极化SAR图像进行分类仿真，分类结果见图6。

由图6可见，本发明对农田的区域划分较细致，边缘保持较好，分类精度高。

综上所述，本发明提出的基于快速密度峰值聚类的极化SAR图像分类方法对极化SAR图像分类能取得更好的分类结果，并可用于对各种极化SAR图像进行分类。

Claims (3)

1.一种基于快速密度峰值聚类的极化SAR图像分类方法，包括如下步骤：

(1)对待分类极化SAR图像进行滤波，去除斑点噪声，得到滤波后的极化SAR图像；

(2)对滤波后的极化SAR图像中每个像素点的相干矩阵T进行Yamaguchi分解，得到每个像素点的体散射功率P_v、二面角散射功率P_d、表面散射功率P_s和螺旋体散射分量P_h；

(3)按照每个像素点四个散射功率P_s,P_d,P_v,P_h的最大值将极化SAR图像初始划分为四类：

如果max(P_s,P_d,P_v,P_h)＝P_s，则将其对应的像素点划分为一类，其中P_s为该类的主散射功率；

如果max(P_s,P_d,P_v,P_h)＝P_d，则将其对应的像素点划分为一类，其中P_d为该类的主散射功率；

如果max(P_s,P_d,P_v,P_h)＝P_v，则将其对应的像素点划分为一类，其中P_v为该类的主散射功率；

如果max(P_s,P_d,P_v,P_h)＝P_h，则将其对应的像素点划分为一类，其中P_h为该类的主散射功率；

其中，max(·)表示最大值；

(4)对每类初始类别所有像素点的主散射功率按由小到大排序，并将每类中每300个像素点再划分为一类，最终将整个极化SAR图像划分为M类；

(5)在已获得的M个类别中，将每一类的中心点作为新像素点，得到M个新像素点A_i，i＝1,…,M，并用每类的新像素点A_i代表该类别内所有像素点；

(6)对上述M个新像素点进行快速密度峰值聚类，将M个新像素点聚为k类：

6a)计算新像素点A_i的局部密度ρ_i：

其中，当d_ij＜d_c时，χ(d_ij-d_c)＝1，否则χ(d_ij-d_c)＝0；M表示新像素点A_i的个数；d_ij表示任意两点间的相互距离：

d_ij＝Tr((T_i)^-1+(T_j)^-1T_i)-q

其中，T_i，T_j分别表示新像素点A_i和A_j的相干矩阵，(T_i)^-1和(T_j)^-1分别表示对矩阵T_i和T_j求逆，q为常数，取值为q＝3，Tr(·)是矩阵的迹；d_c为常数，其取值为将所有点的相互距离d_ij由小到大排列，将位于该排列2％位置处的值作为d_c取值；

6b)计算新像素点A_i的距离δ_i：

6c)选择局部密度ρ_i和距离δ_i乘积最大的k个新像素点作为聚类中心；

6d)聚类中心确定后，每个聚类中心代表一类，k个聚类中心共代表k类，比较剩余新像素点与每一类聚类中心的距离，如果该新像素点与第m类聚类中心的距离最近，则将该新像素点划分到第m类,m＝1,…,k；

(7)在M个新像素点的聚类结果中，将由新像素点A_i代表的所有像素点标记为与新像素点A_i相同的类别，完成对整幅图像的预分类；

(8)对整个极化SAR图像的预分类结果进行复Wishart迭代，得到分类结果。

2.根据权利要求书1所述的基于快速密度峰值聚类的极化SAR图像分类方法，其中步骤(2)所述的对每个像素点的相干矩阵T进行Yamaguchi分解，按如下步骤进行：

(2a)读入滤波后图像的每个像素点，这些像素点为一个3×3的相干矩阵T，根据相干矩阵T得到协方差矩阵C：

其中，U是中间变量，H表示水平极化，V表示垂直极化，S_HH表示水平向发射和水平向接收的回波数据，S_VV表示垂直向发射和垂直向接收的回波数据，S_HV表示水平向发射垂直向接收的回波数据，()^*表示这个数据的共轭，<·>表示按视数平均；

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

其中，f_s为平面散射分量的分解系数，f_d为二面角散射分量的分解系数，f_v为体散射分量的分解系数，f_h为螺旋散射分量的分解系数，β是水平发射水平接收后向散射反射系数与垂直发射垂直接收后向散射发射系数的比值，α被定义为α＝I_HI_VH/I_VI_VV，I_H和I_V分别表示地表的水平及垂直反射系数，I_VH和I_VV分别表示竖直墙体的水平及垂直反射系数，j表示虚数；

(2c)由协方差矩阵C分解的结果，获得一个具有六个未知系数f_s,f_v,f_d,f_h,α,β和五个方程的方程组如下：

(2d)计算像素点协方差矩阵C中的的值，如果则令α＝-1，如果则令β＝1，给定α或β的值后，剩余的4个未知数则根据(2c)中方程组求解，其中Re(·)表示取实部；

(2e)根据求解出的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+|a|²)

P_v＝f_v

P_h＝f_h。

3.根据权利要求书1所述的基于快速密度峰值聚类的极化SAR图像分类方法，其中步骤(8)所述的对整个极化SAR图像的预分类结果进行复Wishart迭代，按如下步骤进行：

(8a)对整个极化SAR图像预分类得到的k类划分结果，根据下式求每一类的聚类中心B_c：

其中T_ρ是第c类中各像素点的相干矩阵，n_c是属于第c类的像素点的个数；

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

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

(8c)比较每个像素点与各聚类中心的距离，如果该新像素点与第m类聚类中心的距离最近，则将该像素点划分到第m类,m＝1,…,k，完成对预分类后整幅极化SAR图像类别的重新划分；

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