1 引言 极化合成孔径雷达(Polarimetric Synthetic Aperture Radar, PolSAR)是一种先进的对地观测技术，不易受天气时间等因素的影响，能够长期获取感兴趣目标的几何结构和物理特性等信息[1]。作为极化SAR数据信息提取和自动化解译的重要步骤，极化SAR图像分类技术得到了深入的研究。目前，极化SAR图像分类已被广泛应用于土地覆盖分类、灾害监测、地质勘探和城市规划等领域[2]。通常，极化SAR图像分类方法按照是否需要标记样本可以分为监督分类、半监督分类和非监督分类。监督分类方法需要根据已有的先验知识挑选大量的训练样本，因此在缺少先验信息的情况下非监督分类方法更加实用。 极化SAR图像非监督分类方法总地来说可以归纳为两大类。第1类方法通过分析极化散射机理，将极化目标分解理论和统计分布模型相结合来完成分类任务。Lee等人[3]利用Freeman极化分解特征初始化经典的Wishart分类器，在保持每类地物散射特性的同时得到了稳定的分类结果。在此基础上，Ferro-Famil等人[4]利用H/A/a-Wishart分类器对多波段全极化SAR图像进行非监督分类，进一步提高了分类精度。第2类方法主要依赖于图像处理和聚类分析技术。Ersahin等人[5]首先将图理论中的谱聚类方法应用于PolSAR图像的非监督分类，而Kersten等人[6]则利用最大期望聚类方法来得到PolSAR的分类结果，这些方法都取得了良好的分类效果。 但是上述基于单个像素的分类方法容易受到相干斑噪声的影响，在分类过程中利用像素之间的空间相关信息能够得到更加鲁棒的分类结果。Yang等人[1]基于极化协方差矩阵的黎曼几何特性，结合黎曼稀疏编码和稀疏相似性进行PolSAR图像分类。而Song等人[7]利用超像素和谱聚类方法完成了大尺度遥感影像的非监督分类任务。Wang等人[8]将张量聚类分析和马尔科夫场相结合，在分割PolSAR图像的过程中有效地融合了边缘信息。这些方法通过在分类过程中引入空间信息，在一定程度上削弱了相干斑噪声的影响，提高了非监督分类的精度。 合适的分类数目对于非监督分类的性能具有重要的影响。研究人员已经提出了多种模型选择方法来求解最优类别数，例如贝叶斯信息准则和最小描述长度准则。Rodriguez等人[9]提出一种密度峰值快速搜索聚类(Density Peaks Clustering, DPC)算法，能够有效地确定聚类数目，实验结果表现出了良好的性能。针对PolSAR图像，Tran等人[10]和Cao等人[11]基于凝聚层次聚类理论提出了自适应类别数的非监督分类方法，伪似然信息和数据的对数似然函数被用于聚类后的验证。在文献[12]基于超像素的分类框架中，Liu等人通过超像素之间的逐对相似性信息来估计类别数，取得了相当高的分类精度。 本文针对PolSAR图像的非监督分类问题，基于混合Wishart模型和密度峰值快速搜索聚类理论提出一种区域级的极化SAR图像非监督分类方法。首先，使用SLIC过分割算法[13]将极化SAR图像分割成多个超像素区域；然后利用混合Wishart模型对每一个超像素区域进行建模，并通过Cauchy-Schwarz (CS)散度来衡量逐对超像素区域之间的距离；最后通过密度峰值快速搜索聚类确定分类数目并得到最终的非监督分类结果。在EMISAR和AIRSAR数据上的实验结果验证了本文方法的有效性。 2 方法介绍 2.1 极化SAR数据和复Wishart分布 在单基站极化SAR测量中，考虑到互易定理，交叉极化分量${S_{{\rm{hv}}}} = {S_{{\rm{vh}}}}$，此时像元内目标的极化信息可以由式(1) 中的复向量进行表示，其中h和v分别表示水平极化和垂直极化，T表示转置运算。 $ \boldsymbol{k} = {\left[{{S_{{\rm{hh}}}}, \sqrt 2 {S_{{\rm{hv}}}}, {S_{{\rm{vv}}}}} \right]^{\rm{T}}} $ (1) 对于经过多视处理的极化SAR数据，每一个数据点都可以由协方差矩阵C来表示，其表达式如下： $ {\boldsymbol{C}} = \frac{1}{N}\mathop \sum \limits_{i = 1}^N {{\boldsymbol{k}}_i}{\boldsymbol{k}}_i^{\rm{H}} $ (2) 其中，N为图像视数，H表示厄密特转置运算。 协方差矩阵C满足复Wishart分布，其概率密度函数为： $ \begin{aligned} & {p{_\boldsymbol{C}}} \left( {{\boldsymbol{C}}{\rm{|}}n, \mathit{\boldsymbol{ \boldsymbol{\varSigma} }}} \right) \\ & \quad \ \;= {n^{qn}}{\left| {\boldsymbol{C}} \right|^{n - q}}\exp \left[{-n \cdot {\rm{Tr}}\left( {{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}^{-1}}{\boldsymbol{C}}} \right)} \right]\\ & \quad \;\; \;\;\;/K\left( {n, q} \right){\left| \mathit{\boldsymbol{ \boldsymbol{\varSigma} }} \right|^n}, \\ & K\left( {n, q} \right) = {{\rm{π}} ^{\frac{{q\left( {q - 1} \right)}}{2}}}{\rm{\Gamma }}\left( n \right) ·\!·\!· {\rm{\Gamma }}\left( {n - q + 1} \right) \end{aligned} $ (3) 其中，q在单站极化数据中为3, n为等效视数，Tr(·)表示矩阵的迹，K为归一化因子，$\mathit{\boldsymbol{ \boldsymbol{\varSigma} }} \!=\! E\left( {\boldsymbol{C}} \right)$表示期望参数。 指数分布族包含高斯分布、Gamma分布和多项式分布等常用的统计分布，其概率密度函数的标准形式为： $ p\left( {x;\lambda } \right) \!=\! {p_F}\left(\! {x;\theta } \right) \!=\! {\rm{exp}}\!\left( \left\langle \right. {\!{t\left( x \right), \theta } \left. \right\rangle \!-\! F\left( \theta \right) \!+\! k\left( x\right)} \right) $ (4) 式中, λ是原参数，t(x)表示充分统计量，θ被称为自然参数，$\left\langle \cdot \right\rangle $指代内积运算，F(·)表示对数归一化因子。复Wishart分布${W_{\boldsymbol{C}}}\left( {n, p, \mathit{\boldsymbol{ \boldsymbol{\varSigma} }}} \right)$也属于指数分布族： $ \begin{aligned} {W_{\boldsymbol{C}}}\left( {n, p, \mathit{\boldsymbol{ \boldsymbol{\varSigma} }}} \right) =&\exp \left( \left\langle \right.{ \!-\! {\boldsymbol{X}}, n{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}^{ - 1}}\left. \right\rangle_F + \left\langle \right.{\rm{log}}\left| {\boldsymbol{X}} \right|, n} \left. \right\rangle \right)\\ & - \left( {{\rm{log}}{{\rm{\Gamma }}_d}\left( n \right) + n{\rm{log}}\left| \mathit{\boldsymbol{ \boldsymbol{\varSigma} }} \right| - dn{\rm{log}}n} \right) \\ & - d{\rm{log}}\left| {\boldsymbol{X}} \right| \end{aligned} $ (5) 对于极化SAR图像区域建模与场景分类问题，参数n是未知常量，只需针对每个区域估计参数$\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}$，而n采用全局估计量。在这种情况下，复Wishart分布${W_{\boldsymbol{C}}}\left( {n, p, \mathit{\boldsymbol{ \boldsymbol{\varSigma} }}} \right)$的参数可以表示为如下形式： $ \theta = \left( {{\theta _s}, {\theta _M}} \right) = \left( {n, n{\mathit{\boldsymbol{ \boldsymbol{\varSigma} }}^{ - 1}}} \right), \\ t\left( x \right) = \left( {\log \left| x \right|, - x} \right), \\ F\left( \theta \right) = \log {{\rm{\Gamma }}_p}\left( n \right) + n{\rm{log}}\left| \mathit{\boldsymbol{ \boldsymbol{\varSigma} }} \right| - pn{\rm{log}}n $ (6) 2.2 混合Wishart模型和Cauchy-Schwarz散度 极化SAR图像的数据通常可以由均值协方差矩阵C来表示，并且服从复Wishart分布C~${W_{\boldsymbol{C}}}\left( {n, p, \mathit{\boldsymbol{ \boldsymbol{\varSigma} }}} \right)$。假定通过SLIC算法已经将PolSAR图像分割成多个超像素区域。对于一个超像素区域，可以通过中心化来表示该区域，例如用超像素m内所有像素的协方差矩阵的平均值Cm来表示。然而这种单模型方式并不适合描述异质性区域，因此需要更加合适的理论模型来对异质性区域进行建模。相较于单一的Wishart模型，由简单模型形成的混合模型更加灵活，且不会涉及到很复杂的参数求解问题。本文使用混合Wishart模型[14]来描述超像素区域： $ m\left( {{\boldsymbol{C}};\mathit{\Phi }} \right) = \mathop \sum \limits_{i = 1}^K {w_i}{p_F}\left( {{\boldsymbol{C}}, {\theta _i}} \right) $ (7) 其中，K表示混合模型的分量个数，$\mathit{\Phi } \!=\! \left( {{w_1}, \!\cdots\!, {w_K}, } \right.$$\left. {{\theta _1}, \cdots, {\theta _K}} \right)$是混合模型的未知参数，${w_1}, \cdots, {w_K}$是每种分量的分布权重，其约束条件为非负数并且和为1。${p_F}\left( {{\boldsymbol{C}};{\theta _i}} \right)$是第i种Wishart分布的概率密度函数，其参数${\theta _i} = \left( {{n_i}, {\mathit{\Sigma}_i}} \right)$。 对于每一个Wishart分布，有$\left( {{w_i}, {n_i}, {\mathit{\Sigma}_i}} \right)$等3个参数需要求解，那么整个超像素区域将会有3K个参数需要训练，因此估计整幅图像的模型参数需要花费大量的时间。本文使用k最大似然估计(k-MLE)算法[15]对混合模型的参数进行快速有效地训练。k-MLE算法将数据迭代分配给最可能的权重部分，对于给定的混合权重使用最大似然估计量更新模型分量。对于独立同分布的观测数据集$x = \left\{ {{x_1}, \cdots, {x_N}} \right\}$ $ =\left\{ {{C_1}, \cdots, {C_N}} \right\}$，将$L\left( {\theta ;x} \right)$定义为似然比函数，$\bar l\left( {\theta ;x} \right)$定义为平均对数似然比函数： $ \begin{aligned} L\left( {\theta ;x} \right) =&\prod\limits_{i = 1}^N {p_F}\left( {{x_i}, \theta } \right) = \prod\limits_{i = 1}^N {\rm{exp}}\left\{ {\left\langle {t\left( {{x_i}} \right), \theta } \right\rangle } \right.\\ & \left. { + k\left( {{x_i}} \right) - F\left( \theta \right)} \right\} \end{aligned} $ (8) $ \ \;\bar l\left( {\theta ;x} \right) = \frac{1}{N}\mathop \sum \limits_{i = 1}^N \left\{\left\langle \right. {t\left( {{x_i}} \right), \theta \left. \right\rangle + k\left( {{x_i}} \right) - F\left( \theta \right)} \right\} $ (9) 对于混合Wishart模型，第i个分量参数的最大似然估计为$\hat \theta = {\rm{argma}}{{\rm{x}}_\theta}\bar l\left( {\theta ;x} \right)$，其中$\theta$满足如下条件： $ \nabla F\left( {\hat \theta } \right) = \frac{1}{N}\mathop \sum \limits_{i = 1}^N t\left( {{x_i}} \right) $ (10) 因此，可以得到最大似然估计$\hat \theta = {\left( {\nabla F} \right)^{ - 1}}$ $ \cdot \left[{\left( {1/N} \right)\sum\nolimits_{i = 1}^N t \left( {{x_i}} \right)} \right]$。令$\eta = \left( {{\eta _S}, {\eta _M}} \right) = E\left( {t\left( x \right)} \right)$为期望参数，那么模型参数可以由式(11) 估计： $ \begin{array}{l} {{\hat \eta }_M} = - {\theta _S}\, \theta _M^{ - 1}, \\ {{\hat \eta }_S} = \mathop \sum \limits_{j = 1}^p \mathit{\Psi }\left( {{\theta _S} - j + 1} \right) - {\rm{log}}\left| {{\theta _S}} \right| \end{array} $ (11) 其中，$\mathit{\Psi }\left( x \right)$是Digamma函数。将式(11) 代入式(6)，可以求解参数${\lambda _i} = \left( {{n_i}, {\mathit{\Sigma}_i}} \right)$。求得Wishart模型的全局参数之后，即可对每一个超像素区域训练混合模型。超像素区域M的权重参数${w_M} = \left( {{w_{{M_1}}}, \cdots, {w_{{M_K}}}} \right)$可以由该区域属于混合模型不同分量的样本数量之和的归一化值来表示。 在使用混合Wishart模型完成PolSAR图像的拟合建模之后，超像素区域之间的距离可以由混合Wishart模型之间的差异度来衡量。衡量分布模型之间的差异可以利用信息论散度，但是针对混合模型的KL散度并没有解析解，因此本文采用Cauchy-Schwarz(CS)散度[16]来衡量混合模型之间的距离： $ {d_{{\rm{CS}}}}\left( {p, q} \right) = - \log \frac{{\int {p(x)q(x){\rm{d}}x} }}{{\sqrt {\int {p{{\left( x \right)}^2}{\rm{d}}x} \int {q{{(x)}^2}{\rm{d}}x} } }} $ (12) 该散度由Cauchy-Schwarz不等式转化而来，只需要考虑$p\left( x \right)q\left( x \right)$的积分。针对指数分布族有限混合模型，${\int {p(x)q(x){\rm d}x} }$的表达式如下： $ \int\!\! {p(x)q(x){\rm d}x} \!\!=\!\! \mathop \sum \limits_{i = 1}^k \mathop \sum \limits_{j = 1}^k\! w_i^pw_j^q{{\rm e}^{F\left( {\theta _i^p + \theta _j^q} \right) - \left( {F\left( {\theta _i^p} \right) + F\left( {\theta _j^q} \right)} \right)}} $ (13) 2.3 密度峰值聚类(DPC)算法 密度峰值快速搜索聚类算法[9, 17]可以给出样本聚类数目的参考值，实现样本集的快速聚类。该算法假设理想的聚类中心包含两个基本特性：(1) 聚类中心被局部密度小于它的邻居点所环绕；(2) 不同聚类中心之间的相对距离较大。该算法通过定义数据量ρi和δi来表征聚类中心的两个基本特性，其中数据量ρi表示数据点i的局部密度，而δi则是数据点i到局部密度大于它的数据点j的最小距离。 相比于计算数据点精确的局部密度值，该算法更加关心局部密度的相对大小，因此局部密度ρi的定义如下： $ {\rho _i} = \mathop \sum \limits_j \chi \left( {{d_{ij}} - {d_c}} \right) $ (14) 其中，dij表示数据点i和j之间的距离，dc称为截断距离。当dij＞dc时，$\chi \left( x \right) = 0$，否则$\chi \left( x \right) = 1$。该定义式说明局部密度值ρi为到数据点i的距离小于截断距离dc的数据点的个数。参考文献[9]中的建议，我们在实验中选择合适的截断距离dc使得平均每个样本的邻域点个数为样本总数的2%。数据点i到其最近邻高局部密度值点j的距离的定义式为： $ {\delta _i} = \mathop {\min }\limits_{j:{\rho _j} > {\rho _i}} {d_{ij}} $ (15) 对于局部密度值最大的数据点i，其${\delta _i} = {\max}\!{_j}{d_{ij}}$。 计算出每个样本点的ρi和δi值后，以ρi为横坐标，δi为纵坐标绘制决策图。真正的聚类中心拥有较大的ρi和δi值，反映在决策图上即为右上方的离散点。根据决策图确定聚类中心之后，剩余的数据点根据每一类的边界域被一次性分配到距其最近的聚类中心。 3 实验及结果分析 本文选取两幅不同的极化SAR图像作为实验数据集。第1组实验数据为丹麦Foulum地区的EMISAR机载L波段全极化SAR图像，其覆盖区域主要包括农田、森林和一些建筑物。该幅图像的大小为300×150像素，其Pauli基(红色|Shh–Svv|，绿色|Shv|，蓝色|Shh+Svv|)合成伪彩图如图 1(a)所示。为抑制相干斑效应，使用窗口大小为5×5的Boxcar滤波器对极化SAR图像进行预处理。SLIC算法的尺度参数NS=10，正则化参数Nm=0.5。NS值越大，超像素越大，Nm越大，超像素越紧凑。在混合Wishart模型的全局训练阶段，当混合模型的分量个数很小时，超像素区域不能被精确地建模；而当模型分量的个数很大时，计算复杂度会相应增加，但是分类性能却没有明显变化。考虑到混合Wishart模型的建模精度和计算复杂度，模型分量的个数K设置为25。本文的对比方法有两组：一组是经典的Wishart-Kmeans聚类算法；另一组是基于超像素之间Bartlett距离[16]的密度峰值快速搜索聚类算法，其中每一个超像素由位于该超像素内所有像素点的协方差矩阵的平均值来表示。 图 1 EMISAR数据的实验结果 Fig.1 The experiment results of EMISAR data 在EMISAR数据上的实验结果如图 1所示。图 1(b)为超像素分割结果，图 1(d)为Bartlett方法的决策图，图 1(e)为本文方法的决策图。从决策图中可以观察到，Bartlett方法将整幅图像分为6类，本文方法将整幅图分为7类，它们的分类结果分别如图 1(f)和图 1(h)所示。比较它们的结果可以看到，两种方法都将各类地物基本区分开来，但是本文方法能够将不同高度的树木区分开来(类3和类4)，而Bartlett方法则将它们分为同一类(类3)。为了更好地分析不同方法的分类性能，我们在图 1(c)中标记了两块异质性比较明显的区域。Bartlett方法将区域2基本分为了同一类(类3)，Wishart方法在区域1和区域2的分类效果都较差，而本文方法不仅将区域2成功区分开来(类3和类4)，而且在标示区域的分类结果也更加平滑。第1组数据可视化的实验结果表明本文方法是有效的。 为了进一步评估本文所提出的方法，我们对另一组实验结果进行了定量分析。第2组实验数据是在1989年获取于荷兰Flevoland地区的AIRSAR机载L波段全极化SAR图像，数据所覆盖的农业区域包含多种典型的农作物。该幅图像的大小为400×400像素，Pauli基伪彩图如图 2(a)所示，图 2(b)为实验区域真实地物类型参考图，其中定义了各种颜色代表的地物种类，包括甜菜、草地、小麦、土豆、大麦、苜蓿、裸地、油菜和豌豆等9类地物。混合Wishart模型分量的个数为K=15。本文选择标记地物的分类精度、全精度(OA)、Kappa系数、F1-score值和纯度(Purity)[18]作为非监督分类结果定量分析的评价指标，其中OA、Kappa系数、F1-score和Purity的数值越大，表明分类精度越高，分类算法的性能越好。 图 2 AIRSAR数据 Fig.2 AIRSAR data AIRSAR数据的实验结果如图 3所示。图 3(a)为超像素分割结果，图 3(b)和图 3(c)分别为Bartlett方法与本文方法的决策图。图 3(b)中有7个明显的类中心，整幅图像被分为7类。而图 3(c)中则有9个类中心，整幅图像被分为9类。它们的分类结果分别如图 3(d)和图 3(f)所示。其中图 3(g)－图 3(i)为只含有标记区域的分类结果图。参考图 2(b)中真实地物图能够发现，Bartlett方法的分类结果中只有小麦、大麦、苜蓿和油菜等地物的分类基本正确，而土豆、豌豆和草地则存在不同程度的分类错误，并且整幅图像的大部分区域都被错分成为一大类，分类效果比较差。而在本文方法的分类结果中，除了部分甜菜区域被错误分类成大麦之外，其他几类地物的绝大部分区域都得到了正确的分类。相比于图 3(d)中的分类结果，在图 3(g)中绝大部分地物的分类精度都得到了提升。与图 3(e)和图 3(h)中Wishart方法的分类结果相比，本文方法的分类结果更加平滑，匀质性区域分为多类的情况明显减少，例如图 3(h)中大麦和油菜中都出现了明显的错分现象。可视化的实验结果表明，本文方法的分类性能比其他两种方法更加良好。 图 3 AIRSAR数据的实验结果 Fig.3 The experimental results of AIRSAR data 为了更加准确地评价不同方法的分类性能，表 1中列出了AIRSAR数据中标记区域的定量分析结果。除了甜菜的分类精度较低之外，本文方法对其他地物的分类精度都达到了90%以上，其中OA值和Kappa系数分别为0.9420和0.9336，而且F1-score值和Purity值也分别达到了0.9337和0.9432，这些评价指标都远高于其他两种对比方法。和可视化的分类结果相似，数字化的定量分析指标再一次表明了本文方法良好的分类性能。 表 1 AIRSAR数据的分类结果 Tab.1 The classification results of AIRSAR data 本文方法主要分为3个步骤，第1步是使用简单线性迭代聚类(SLIC)算法生成超像素，第2步是利用混合Wishart模型和Cauchy-Schwarz散度计算不同超像素区域之间的距离(MixWishart)，第3步是使用密度峰值快速搜索聚类(DPC)得到非监督分类结果。本文实验中EMISAR数据大小为300×150像素，AIRSAR数据大小为400×400像素。表 2为不同分量个数的情况下，混合Wishart模型在EMISAR和AIRSAR数据上求解的计算时间。表 3列出了本文方法各个步骤在EMISAR和AIRSAR数据上所对应的计算时间。相关程序使用Matlab2014进行编程，执行环境为3.4 GHz的单核Intel CPU。表 2和表 3中的结果都是取10次试验的平均值，从其中的计算时间可以观察到，本文方法的运算效率较高，具有较强的实用性。 表 2 混合Wishart模型不同分量个数的计算时间 Tab.2 The computational time of Wishart mixture models with different components 表 3 不同步骤的计算时间 Tab.3 The computational time corresponding to each step 4 结论 本文针对极化SAR图像的非监督分类问题，提出了一种基于混合Wishart模型和密度峰值快速搜索聚类的区域级分类方法。在EMISAR和AIRSAR实验数据上的分类结果表明，该方法能够有效地提高PolSAR图像非监督分类的效果，同时基于区域级的分类也能有效地降低相干斑噪声的影响，得到更加鲁棒的分类结果。在进一步的工作中，我们将探索如何自适应地确定混合Wishart模型的分量个数。