1 引言

极化合成孔径雷达(Polarimetric Synthetic Aperture Radar, PolSAR)是一种全天时、全天候、能够大范围获取多维度、高分辨率地面场景信息的遥感对地观测系统。作为极化SAR数据解译的核心内容之一，极化SAR图像分类在军事和民用领域具有广阔的应用前景。协方差矩阵是极化SAR数据最常用的表达形式之一，如何度量协方差矩阵之间的相似性已成为极化SAR图像解译的关键问题之一。通过多年研究，学者们普遍认为在现有的协方差矩阵距离测度中，没有对所有应用均适合的最优测度。在文献[1]中，Yang等人针对不同的极化SAR图像应用详细探讨了多种相似度和距离测度。基于极化协方差矩阵的复Wishart分布模型的对数似然函数，Lee等人^[2]提出Wishart距离并将其应用于极化SAR图像的分类任务中。进一步地，Frery等人^[3]基于复Wishart分布模型讨论了4种常用的概率距离并给出了它们的解析表达式。针对Wishart距离并不严格满足数学意义上所定义的测度性质，Anfinisen等人^[4]又提出了对称修正的Wishart距离。Banerjee等人从指数族函数分布模型的角度出发定义了协方差矩阵的另一种距离度量方法——Bregman散度^[5]。基于假设检验理论，Kersten等人^[6]使用Bartlett距离来度量两个协方差矩阵之间的差异，Song等人^[7]将Bartlett距离应用于极化SAR图像的大尺度谱聚类框架中。考虑到极化协方差矩阵的黎曼流形结构，Song等人^[8]提出了一种结合Jensen-Bregman散度(也称为Stein散度)^[9]与K-means聚类的极化SAR图像非监督分类方法。

稀疏编码理论通过将信号分解为字典中少数原子的线性组合，从而对信号进行编码。目前该理论在极化SAR图像处理任务中得到了广泛应用。鉴于稀疏编码算法主要针对欧式空间中的特征矢量，研究人员将其应用范围从向量值数据扩展到更复杂的数据形式，比如对称正定(Symmetric Positive Definite, SPD)矩阵和厄密特正定(Hermitian Positive Definite, HPD)矩阵。在极化SAR图像分类任务中，考虑到极化协方差矩阵黎曼流形结构的非线性特点，不能直接将向量值数据稀疏编码进行扩展。一种常用的方法是提取极化SAR数据的特征向量，例如Zhang等人^[10]通过提取多种极化特征并对组合的特征向量进行稀疏编码从而实现极化SAR图像监督分类的任务。然而该方法需要人工选择初始特征，不能充分利用极化协方差矩阵的信息。

因此如何充分挖掘极化协方差矩阵所包含的丰富信息已经成为极化SAR图像解译领域的核心问题之一。在计算机视觉领域，Harandi等人^[11]通过将流形嵌入到再生核希尔伯特空间(Reproducing Kernel Hilbert Space, RKHS)以求解黎曼流形上的稀疏编码和字典学习问题。通过利用Stein核进行嵌入操作，文献[12,13]成功地将稀疏表达分类器应用于极化SAR图像地物分类任务。在此基础上，通过进一步讨论极化SAR数据在高维再生核希尔伯特空间中的结构问题，我们提出了一种能够有效融合空间信息的极化SAR图像监督分类方法。和希尔伯特空间嵌入的稀疏编码方法不同，Cherian等人^[14]首先提出了针对SPD矩阵的黎曼稀疏编码方法，我们将其进一步扩展到HPD矩阵，并提出了基于极化协方差矩阵黎曼稀疏诱导相似度^[15]的极化SAR图像非监督分类方法。实验结果表明本文所提出的监督分类与非监督分类方法有效地提高了极化SAR图像的分类精度。

2 极化协方差矩阵和黎曼流形测度

对于单基站极化合成孔径雷达测量目标，每个分辨率单元内的全部极化信息可以通过复散射向量 ${k} = {\left[ {{S_{{\rm{hh}}}},\sqrt 2 {S_{{\rm{hv}}}},{S_{{\rm{vv}}}}} \right]^{\rm{T}}}$ 进行表示，其中h表示水平极化，v则表示垂直极化，而上角标T表示向量的转置运算。通常情况下，极化SAR数据都会经过多视处理以抑制相干斑效应。因此，多视极化SAR数据能够由如下形式的协方差矩阵C所表征：

$\begin{aligned}\!\! {C} = & \frac{1}{N}\mathop \sum \limits_{i = 1}^N {{k}_i}{k}_i^{\rm{H}}\\ = & \left[\!\!\!\! {\begin{array}{*{20}{c}}{\left\langle {{{\left| {{S_{{\rm{hh}}}}} \right|}^2}} \right\rangle } & {\sqrt 2 \left\langle {{S_{{\rm{hh}}}}S_{{\rm{hv}}}^*} \right\rangle } & {\left\langle {{S_{{\rm{hh}}}}S_{{\rm{vv}}}^{\rm{*}}} \right\rangle }\\{\sqrt 2 \left\langle {{S_{{\rm{hv}}}}S_{{\rm{hh}}}^*} \right\rangle } & {\left\langle {{{\left| {{S_{{\rm{hv}}}}} \right|}^2}} \right\rangle } & {\sqrt 2 \left\langle {{S_{{\rm{hh}}}}S_{{\rm{vv}}}^{\rm{*}}} \right\rangle }\\{\left\langle {{S_{{\rm{vv}}}}S_{{\rm{hh}}}^*} \right\rangle } & {\sqrt 2 \left\langle {{S_{{\rm{vv}}}}S_{{\rm{hh}}}^*} \right\rangle } & {\left\langle {{{\left| {{S_{{\rm{vv}}}}} \right|}^2}} \right\rangle }\end{array}} \!\!\!\!\right]\end{aligned}$

(1)

其中，N表示图像视数，H表示厄密特转置运算，*表示共轭操作符。

为了充分利用协方差矩阵所包含的信息，研究人员提出了多种协方差矩阵的距离度量方法。除了前文提及的基于统计分布的距离度量方法^[2–6]，我们还可以从几何的角度来定义协方差之间的距离。经过多视处理的极化协方差矩阵为厄密特正定矩阵，它们所构成的黎曼流形并不满足欧式空间的几何性质，因此简单地将HPD矩阵线性化并使用欧氏距离来度量的方法在实际应用中表现出较差的性能^[6]。考虑到黎曼流形的几何特性，HPD矩阵之间的相似度可以用测地线距离来度量。对于黎曼流形上的两个HPD矩阵X和Y，一种常用的测地线距离为仿射不变黎曼测度^[16](Affine Invariant Riemannian Metric, AIRM)，其定义式为：

${d_{{\rm{AIRM}}}}\left( {X,Y} \right) = {\rm{ }}{\left\| {{\rm{log}}\left( {{X^{ - 1/2}}Y\;{X^{ - 1/2}}} \right)} \right\|_{\rm{F}}}$

(2)

其中， ${\left\| {{\rm{log}}\left( \cdot \right)} \right\|_F}$ 表示Frobenius范数， ${{\rm{log}}\left( \cdot \right)}$ 为矩阵对数运算。对于黎曼流形，AIRM是应用最为广泛的一种距离度量。然而由于其涉及特征值计算，具有较高的计算复杂度，因此对数欧式黎曼测度^[17](Log-Euclidean Riemannian Metric, LERM)常常作为另一种可供选择的替代方案：

${d_{{\rm{LERM}}}}\left( {{X},{Y}} \right) = {\left\| {\log \left( {X} \right) - {\rm{log}}\left( {Y} \right)} \right\|_F}$

(3)

该距离度量通过对数变换将流形数据映射到切平面空间从而利用欧式距离衡量样本间的相似性。由于对数欧式变换可以进行离线计算，因此LERM的计算复杂度相对于AIRM有所降低。但是LERM利用切平面映射平坦化流形也导致距离计算不够精确。为了进一步降低计算复杂度，研究者又提出了Bartlett距离^[6]：

${d_{\rm{B}}}\left( {{X},{Y}} \right) = {\rm{log}}\left| {\frac{{{X} + {Y}}}{2}} \right| - \frac{1}{2}\left| {{XY}} \right|$

(4)

Bartlett距离也称为Jensen-Bregman对数行列式散度或者Stein散度^[9]。作为AIRM的另一种替代选择，Bartlett距离具有良好的计算效率，并且同LERM一样具有与AIRM相似的性质。

3 极化协方差矩阵稀疏编码

将稀疏编码方法扩展到极化协方差矩阵需要求解黎曼流形上的稀疏编码问题。黎曼流形上的稀疏编码问题可以表述为：将需要编码的流形上的数据点表示为给定的字典元素的稀疏组合。与欧氏空间不同，流形结构本身的非线性导致对应的距离测度的非线性，在其上稀疏组合涉及到非线性运算，这使得直接定义在流形上的相关优化问题的求解变得相对困难。

3.1 核空间稀疏编码方法

希尔伯特空间是欧式空间的推广，它是一个完备的内积空间(可能为无限维)，其中定义了由内积导出的范数。RKHS是一种特殊的希尔伯特空间。对于一个定义在非空集合 ${\cal X}$ 上的RKHS( ${\cal H},{\left\langle \cdot \right\rangle _{\cal H}}$ )以及函数 $\phi :{\cal X} \to {\cal H}$ , ${\cal H}$ 上的内积 ${\left\langle \cdot \right\rangle _{\cal H}}$ 可由二元函数，也即核函数 $k:{\cal X} \times {\cal X} \to {\mathbb{R}}$ 来定义： ${\left\langle {\phi \left( x \right),\phi \left( y \right)} \right\rangle _{\cal K}}= $ k(x, y)，其中 $x,y \in {\cal X}$ 。该性质使得无需知道数据在 ${\cal H}$ 中的具体形式就可以直接计算其在 ${\cal H}$ 上的内积。将低维数据嵌入到高维希尔伯特空间通常都是针对欧式空间的数据，然而Jayasumana等人^[18]的工作却表明这一概念同样适用于流形值数据，例如极化SAR数据的HPD矩阵。具体地说，对于黎曼流形 ${\cal M}$ 上的极化协方差矩阵X和Y，可以通过函数 $\phi :{\cal M} \to {\cal H}$ 映射到高维的希尔伯特空间 ${\cal H}$ ，内积运算由核函数 $k\left( \cdot \right)$ 来确定： ${\left\langle {\phi \left( {X} \right),\phi \left( {Y} \right)} \right\rangle _{\cal K}} = k\left( {{X},{Y}} \right)$ 。

对于极化SAR协方差矩阵数据，其希尔伯特空间嵌入可以通过基于黎曼流形测度的高斯径向基核函数 ${k_{\rm{G}}}\left( \cdot \right)$ 来实现：

${k_{\rm{G}}}\left( {{X},{Y}} \right) = {\rm{exp}}\left( { - \beta d\left( {{X},{Y}} \right)} \right)$

(5)

在高斯径向基核形式下，LERM和Bartlett距离均能获得正定核。由LERM可推导出对数欧式核函数：

${k_{{\rm{LERM}}}}\left( {{X},{Y}} \right) = {\rm{exp}}\left( { - \beta {{\left\| {\log \left( {X} \right) - {\rm{log}}\left( {Y} \right)} \right\|}_F}} \right)$

(6)

由Bartlett距离(Stein散度)可推导出Stein核：

${k_{{\rm{Stein}}}}\left( {{X},{Y}} \right) = {2^{\beta \gamma }}\frac{{\sqrt {{{\left| {X} \right|}^\beta }{{\left| {Y} \right|}^\beta }} }}{{{{\left| {{X} + {Y}} \right|}^\beta }}}$

(7)

其中 $\gamma$ 表示维度，对于单站极化SAR协方差矩阵数据而言， $\gamma = 3$ 。

通过核函数将流形值数据嵌入到希尔伯特空间，极化协方差矩阵的稀疏编码问题便可以转化为希尔伯特空间上的线性组合问题。特别地，给定一个由黎曼流形上数据组成的字典 ${M} = \left\{ {{{D}_1},{{D}_2}, ·\!·\!· ,{{D}_K}} \right\}$ ${{D}_i} \in {\cal M}$ 和一个嵌入核函数 $\phi :{\cal M} \to {\cal H}$ ，对于流形上的极化协方差矩阵X，我们寻找一个稀疏的系数矢量 ${v} \in {{\mathbb{R}}^K}$ 使得 $\phi \left( {X} \right)$ 表示为以v为系数的 $\left\{ {\phi \left( {{{D}_1}} \right), ·\!·\!· ,\phi \left( {{{D}_K}} \right)} \right\}$ 的线性组合。该问题可以表示为核化的Lasso问题^[19]形式：

$\hat{ {v}} = {\rm{argmin}}\left( {{{\left\| {\phi \left( {X} \right) \!-\! \mathop \sum \limits_{i = 1}^K {v_i}\phi \left( {{{D}_i}} \right)} \right\|}^2} \!+\! \lambda {{\left\| {v} \right\|}_1}} \right)$

(8)

由核函数的性质可知，式(8)中的第1项可以表达为如下形式：

$\begin{aligned} & \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! {\left\| {\phi \left( {X} \right) - \mathop \sum \limits_{i = 1}^K {v_i}\phi \left( {{{D}_i}} \right)} \right\|^2} \\ \;\;\;\; = & k\left( {{X},{X}} \right) - 2\mathop \sum \limits_{i = 1}^K {v_i}k\left( {{X},{{D}_i}} \right) \\&+ \mathop \sum \limits_{i,j = 1}^K {v_i}{v_j}k\left( {{{D}_i},{{D}_j}} \right)\end{aligned}$

(9)

如前所述，核函数 $k\left( { \cdot , \cdot } \right)$ 可以采用对数欧式核函数或者Stein核函数。核函数映射将黎曼流形上稀疏编码问题转化为希尔伯特空间的稀疏编码问题，使用经典的交叉方向乘子法^[20](Alternative Direction Method of Multiplier, ADMM)求解该凸优化问题即可得到流形数据在希尔伯特空间的稀疏编码向量v。

3.2 黎曼稀疏编码方法

和希尔伯特空间方法不同，黎曼稀疏编码^[15,21]不需要通过核函数将流形数据映射到高维希尔伯特空间。黎曼稀疏编码可以描述为将给定的极化协方差矩阵分解为多个HPD矩阵原子的稀疏线性组合，并使得重建误差最小。具体来讲，对于给定的极化协方差矩阵X和编码字典 ${M} = \left\{ {{{D}_1},{{D}_2}, ·\!·\!· ,{{D}_K}} \right\},$ ${{D}_i} \in {\cal M}$ ，黎曼稀疏编码的目的是找到一组非负的稀疏系数 ${{α}} = {\left[ {{\alpha _1},{\alpha _2}, ·\!·\!· ,{\alpha _K}} \right]^{\rm{T}}}$ ，使得线性组合 $\sum\nolimits_{i = 1}^K {{\alpha _i}} {D_i}$ 与HPD矩阵X的距离最小。该问题可以表示为：

$\min \phi \left( {{α}} \right): = {d^2}\left( {\mathop \sum \limits_{i = 1}^K {\alpha _i}{{D}_i},{X}} \right) + {\cal L}\left( {{α}} \right)$

(10)

其中， ${{\alpha _i}}$ 表示编码字典中第i个原子D_i的编码系数，d表示距离度量， ${\cal L}\left( {{α}} \right)$ 为损失函数。考虑到极化协方差矩阵的黎曼流形结构，式(10)中距离度量d可以使用AIRM来代替。文献[21]中已经证明，函数 $\phi \left( \alpha \right) = d_{{\rm{AIRM}}}^{\;2}\left( {\sum\nolimits_{i = 1}^K {{\alpha _i}} {D_i},X} \right)$ 在集合 ${\cal A}$ 上为凸函数：

${\cal A}: = \left\{ {{α} \Biggr|\mathop \sum \limits_{i = 1}^K {\alpha _i}{{D}_i} \underline \prec {X},{\alpha _i} \ge 0} \right\}$

(11)

损失函数 ${\cal L}\left( {{α}} \right)$ 可以写为 ${\cal L}\left( {{α}} \right) = \lambda \parallel {{{α}} \parallel_1}$ ，其中 $\lambda $ 为正则化参数。考虑到 ${\alpha _i} \ge 0$ ，我们用 $\lambda \sum\nolimits_i {{\alpha _i}} $ 替换损失函数，黎曼稀疏编码问题可以表示为如下的优化问题：

$\begin{aligned}{\rm{min}}\;\;\phi \left( {{α}} \right): =& \frac{1}{2}\left\| {{\rm{log}}\mathop \sum \limits_{i = 1}^K {\alpha _i}{{X}^{ - 1/2}}{{D}_i}{{X}^{ - 1/2}}} \right\|_F^2 \\ & + {\rm{\lambda }}\mathop \sum \limits_{i = 1}^K {\alpha _i}\end{aligned}$

(12)

通过使用AIRM来度量重建误差有效地提高了模型效果，并且该模型可以通过现有的凸优化工具进行求解。给定HPD矩阵B, C和X，对AIRM距离函数 $f\left( x \right) = d_{{\rm{AIRM}}}^{\ 2}\left( {x{B} + {C},{X}} \right)$ 求导可得：

$\begin{aligned}f'\left( x \right) = & 2{\rm{Tr}} \Big({{\rm{log}}\left( {{S} \left( {x{B} + {C}} \right){S}} \right)}\\& \left. \cdot{{{{{S}^{ - 1}}\left( {x{B} + {C}} \right)}^{ - 1}}{BS}} \right),{S} = {{X}^{ -1/2}}\end{aligned}$

(13)

将 $M\left( x \right) \!=\! x{B} \!+\! {C}$ 改写为 $M\left( {{\alpha _p}} \right) \!=\! {\alpha _p}{{D}_p} \!+\! \displaystyle\mathop \sum \nolimits_{i \ne p}$ ${\alpha _i}{{D}_i}$ ，黎曼稀疏编码问题优化函数的梯度 $\nabla \phi \left( {{α}} \right)$ 可以表示为：

$\begin{aligned}\frac{{\partial \phi \left( {α} \right)}}{{\partial {\alpha _p}}} \!=\!\! & {\rm{Tr}}\!\Biggr(\!\! {\log \!\left( \!\!\!{{{X}^{ - 1/2}}\!\!\left( \!\!{{\alpha _p}{{D}_p} \!\!+\!\! \mathop \sum \limits_{i \ne p} \!\!{\alpha _i}{{D}_i}} \!\!\right)\!\!{{X}^{ - 1/2}}} \!\!\right)\!\!{{X}^{1/2}}} \\& \left.\!\! { \cdot \!\!{{\left( {{\alpha _p}{{D}_p} \!+\! \mathop \sum \limits_{i \ne p} {\alpha _i}{{D}_i}} \right)}^{ - 1}}\!\!{{D}_p}{{X}^{ - 1/2}}} \right)\! + \!\lambda \end{aligned}$

(14)

利用谱投影梯度(Spectral Projected Gradient, SPG)算法^[22]求解该优化问题可以得到HPD矩阵X在黎曼稀疏编码框架下的非负稀疏编码向量 ${{α}}$ 。

4 基于极化协方差矩阵稀疏编码的极化SAR 图像分类方法 4.1 基于核空间稀疏编码的极化SAR图像监督分类方法

稀疏表达分类器最早由Wright等人^[23]在人脸识别任务中提出，其核心思想是来自于同一类别的训练样本位于相同的子空间。测试样本在与其不同类别的训练样本空间的投影系数接近于0，因此整个投影系数表现为具有稀疏性。Yang等人^[12]通过核函数稀疏编码方法成功地将稀疏表达分类器应用于极化SAR图像分类。

在核稀疏表达分类器框架下，字典 ${M}$ 由训练样本组成。假设训练样本来自于M个类别，当求得测试样本X的稀疏编码系数 $\hat{ {v}}$ 之后，可以求出由不同类别的字典元素所重建的误差，则测试样本X的类别标号Label(X)由重建误差最小的类别所确定：

${\rm{Label}}\left( {X} \right) = \mathop {{\rm{argmin}}}\limits_{m = 1, \cdots ,M} {\left\| {\phi \left( {X} \right) - \phi \left( {{{M}_m}} \right){{\hat{ {v}}}_m}} \right\|_F}$

(15)

其中， ${\phi \left( {{{M}_m}} \right)}$ 表示第m类的字典原子通过核函数 $\phi \left( \cdot \right)$ 在希尔伯特空间中的映射， ${\hat{ {v}}_m}$ 为第m类的字典原子在希尔伯特空间的稀疏编码系数。

在上述过程中，编码与分类都是逐像素进行的。然而，极化SAR图像本身存在相干斑噪声，逐像素分类会影响算法的分类性能，因此需要在分类过程中考虑空间信息。我们在文献[13]中将权重联合稀疏表达分类^[24]方法引入核稀疏表达分类器框架，通过考虑邻域像素在稀疏编码中的不同贡献引入空间信息，从而有效地提高了极化SAR图像的分类精度。

假定中心像素点X及其邻域 $\cal{{R}}$ 内的像素点X_S的极化协方差矩阵集合为 $S = \left\{ {X,{X_1},{\rm{ }} \cdots ,{X_S}} \right\}$ ，核稀疏表达分类框架下的权重联合稀疏编码问题可以表示为：

$\min {\left\| {\phi \left( {S} \right){W} - \phi \left( {M} \right){{V}\!\!_w}} \right\|_F} + \lambda {\left\| {{{V}\!\!_w}} \right\|_F}$

(16)

其中，权重矩阵W的值表示邻域像素在稀疏编码过程中的不同贡献，V_w为稀疏编码系数矩阵。在权重联合稀疏表达分类框架下，测试样本X的类别标签仍然由重建误差最小的类别所确定：

${\rm{Label}}\left( {X} \right) \!\! =\!\! \mathop {{\rm{argmin}}}\limits_{m = 1, \cdots ,M} {\left\| {\phi \left( {S} \right){W} \!\!-\!\! \phi \left( {{{M}_m}} \right){V}_w^m} \right\|_F}$

(17)

权重矩阵W由非局部空间信息^[25]所确定。假定中心像素点X的邻域为R×R的矩形窗口，中心像素与其邻域窗口内像素X^r的相似性可以定义为：

$w\left( {{X},{{X}^{\rm{r}}}} \right) = {\rm{exp}}\left( { - \Delta \left( {{X},{{X}^{\rm{r}}}} \right)/h} \right)$

(18)

参数h为权重因子， ${\Delta \left( {{X},{{X}^{\rm{r}}}} \right)}$ 由中心位于X和X^r，窗口大小为t×t的图像块的相似性所确定。对于极化SAR图像， ${\Delta \left( {{X},{{X}^{\rm{r}}}} \right)}$ 为两个图像块之间对应的极化协方差矩阵的Bartlett距离之和所确定：

$\Delta \left( {{X},{{X}^{\rm{r}}}} \right) = \mathop \sum \limits_{i = 1}^{{t^2}} {d_{\rm B}}\left( {{{X}_i},{X}_i^{\rm{r}}} \right)$

(19)

本文所提出的监督分类方法(记为“KSC”)的基本流程包括：(1)通过Stein核函数将极化协方差矩阵映射到希尔伯特再生核空间得到高维特征向量；(2)使用非局部均值滤波算法生成包含空间信息的权重矩阵；(3)结合权重矩阵对再生希尔伯特核空间的高维特征向量进行权重联合稀疏编码；(4)利用训练样本所构建的稀疏表达分类器对极化SAR图像进行监督分类。

4.2 基于黎曼稀疏编码的极化SAR图像非监督分类方法

Cheng等人^[26]提出稀疏诱导相似度并将其应用于向量值数据，本文通过黎曼稀疏编码将其扩展至极化协方差矩阵，提出了黎曼稀疏诱导相似度^[15](Riemannian Sparse Induced Similarity, RSIS)，并将其应用于极化SAR图像的非监督分类任务。

假定HPD矩阵集合 ${\cal F} = \left\{ {{{F}_1},{{F}_2}, ·\!·\!· ,{{F}_N}} \right\}$ 由极化协方差矩阵组成，其中 ${{F}_k} \in {\cal F}$ 为给定的待编码矩阵。首先将除了F_k以外的HPD矩阵作为编码字典，即编码字典满足 ${{M}_k} = {\cal F}\backslash {{F}_k} = \left\{ {{{F}_1},{{F}_2}, ·\!·\!· ,} \right. $ $ \left. {{{F}_{k - 1}},{{F}_{k + 1}}, ·\!·\!· ,{{F}_N}} \right\}$ 。然后使用黎曼稀疏编码算法将HPD矩阵F_k分解成编码字典 ${{M}_k}$ 中的原子的线性组合，得到非负编码向量 $${\alpha _k} = \left[ {{\alpha _1},{\alpha _2},{\rm{ }} \cdots ,} \right.{\left. {{\alpha _{k - 1}},{\alpha _{k + 1}},{\rm{ }} \cdots ,{\alpha _N}} \right]^{\rm{T}}}$$ 。那么HPD矩阵F_i与F_k之间的黎曼稀疏诱导相似度可以定义为：

${s_{ki}} = {\alpha _i}{\rm{\biggr/}}\!\!\sum\limits_{j = 1,j \ne k}^N {{\alpha _j}} $

(20)

和文献[26]中的计算方法不同，黎曼稀疏编码系数是非负的，因此不需要在系数值和0之间使用取最大值函数。对HPD矩阵集合中的所有元素进行同样的相似度计算，得到相似度矩阵 ${S} = \left[ {{{s}_1},{{s}_2}, ·\!·\!· ,{{s}_N}} \right]$ 。考虑到对称性，最终的黎曼稀疏诱导相似度矩阵定义为：

${s}_{ij}^{\rm{R}} = \left\{ {\begin{aligned}& {\frac{{{{s}_{ij}} + {{s}_{ji}}}}{2},\;\;\;{i \ne j}}\\& {1,\;\;\;\;\;\;\;\;\;\;\;\;\; \;\; {i = j} }\end{aligned}} \right.$

(21)

在得到黎曼稀疏诱导相似度矩阵S^R后，可以通过谱聚类算法^[27]得到极化SAR图像非监督分类的结果。考虑到谱聚类算法的计算复杂度，本文首先使用简单线性迭代聚类^[28](Simple Linear Iterative Clustering, SLIC)算法将极化SAR图像分割成多个超像素区域，再使用超像素内所有像素的协方差矩阵的平均值来表示该超像素。基于超像素区域的分类不仅有效降低了黎曼稀疏编码过程中矩阵样本的数量和谱聚类的计算量，同时也引入了部分空间约束信息。

本文所提出的非监督分类方法(记为“RSC”)可以表示为如下步骤：(1)通过SLIC算法对极化SAR图像进行分割，生成多个超像素子区域；(2)使用K近邻算法对每一个超像素区域生成黎曼稀疏编码字典；(3)对每一个超像素区域进行黎曼稀疏编码，计算超像素区域之间的黎曼稀疏诱导相似度矩阵；(4)在黎曼稀疏诱导相似度矩阵上执行谱聚类算法，对极化SAR图像进行非监督分类。

5 实验结果与分析 5.1 监督分类实验结果

图1(a)为实验数据的Pauli基伪彩图(红色： $\left| {{S_{{\rm{hh}}}} - {S_{{\rm{vv}}}}} \right|$ ，绿色： $\left| {{S_{{\rm{hv}}}}} \right|$ ，蓝色： $\left| {{S_{{\rm{hh}}}} + {S_{{\rm{vv}}}}} \right|$ )，图像大小为260×260像素。该数据为EMISAR全极化SAR图像，于1998年获取自丹麦Foulum地区的农田，主要由油菜、黑麦、燕麦、小麦和针叶林等5类地物组成。该区域地物的真实分布如图1(e)所示。我们将基于Wishart距离的最大似然分类器^[29]和无空间信息的核函数稀疏表达分类器^[12]作为对比实验，并使用总体分类精度(Overall Accuracy, OA)和Kappa系数^[30]来衡量不同实验方法的分类性能。本文选用窗口大小为3×3的Boxcar滤波器对极化SAR图像进行预处理以抑制相干斑效应。稀疏编码步骤使用Stein核函数，每类地物随机选择100个样本点作为字典原子。在权重联合编码过程中使用15×15的矩形搜索框，h=1, t=3。

EMISAR数据监督分类的结果如图1所示。比较图1(b)～图1(d)中不同方法的分类结果，我们可以发现图1(d)中的KSC方法分类性能最优。由于在分类过程中通过权重联合编码引入了空间信息，图1(d)中的分类效果比图1(b)～图1(c)中的要好，在同一区域的分类结果更加一致，例如图1(d)中小麦和黑麦区域的分类错误更少，分类结果更加平滑，而另外两种对比方法则存在较多的错误分类。为了进一步说明不同分类方法的性能，图1(f)～图1(h)中给出了真实地物分布区域所对应的分类结果。从图1(f)～图1(h)的分类结果中，可以观察到本文方法对这些类地物的分类几乎完全正确，特别是针叶林、油菜和燕麦等3类地物基本没有错分。尽管小麦和黑麦存在极少量的错分现象，但是相比其他两种对比方法，KSC方法在这两类地物的分类精度有显著提高。

表1中列出了不同实验方法在EMISAR数据标记区域中每类地物的分类精度，并计算出相应的OA值和Kappa系数。从表1中的数据可以看到，KSC方法对针叶林、小麦、油菜、燕麦和黑麦等多种地物的分类精度都较对比方法有所提高，而且OA值为0.9981, Kappa值为0.9975，也远高于另外两种方法。可视化的分类结果和数字化的定量分析指标都充分表明了本文所提出的KSC方法的有效性。

5.2 非监督分类实验结果

图2(a)为实验数据的Pauli基伪彩图(红色： $\left| {{S_{{\rm{hh}}}} - {S_{{\rm{vv}}}}} \right|$ ，绿色： $\left| {{S_{{\rm{hv}}}}} \right|$ ，蓝色： $\left| {{S_{{\rm{hh}}}} + {S_{{\rm{vv}}}}} \right|$ )，图像大小为400×400像素。该数据为AIRSAR全极化SAR图像，于1989年获取自荷兰Flevoland地区，其所覆盖的区域是一片典型的农业区。该区域的地物类型主要为甜菜、苜蓿、裸地、草地、小麦、土豆、油菜、豌豆和大麦。该区域的真实地物分布如图2(e)所示。我们选择基于Wishart距离的K-means聚类方法和基于Bartlett距离的谱聚类方法作为对比实验，并使用OA, F1-score^[31], Purity和Entropy^[32]等4个因子作为衡量不同方法分类性能的指标。其中OA为平均分类精度。F-score是衡量分类性能的有效指标，F1-score就是将准确率和召回率按照等权重合并为一个值，其值越大，表明分类性能越好。Purity是衡量每个簇由单个类的对象组成程度的指标，其值越大，表明簇中含有的类别数越少，分类性能也越好。Entropy是衡量簇类中数据来源分布均衡程度的指标，其值越小，表明数据的来源分布越均衡，分类性能也越好。本文选用窗口大小为3×3的Boxcar滤波器对极化SAR图像进行预处理以抑制相干斑效应。SLIC算法中参数N_s=10, N_m=0.1。黎曼稀疏编码中正则化因子 ${\rm{\lambda }} = 0.1$ ，字典原子数为30。谱聚类算法中聚类数目k=9。

AIRSAR数据非监督分类的结果如图2所示。图2(b)为基于Wishart距离的K-means聚类方法的结果，图2(c)为基于Bartlett距离的谱聚类方法的结果，图2(d)为RSC方法的结果，可以看到RSC方法的分类结果比另外两种方法要好。前两种对比方法的错误分类较多，甚至同一类别的相邻区域也被错分为了不同的地物类别。由于黎曼稀疏诱导相似度有效地融合了超像素的邻域信息，RSC方法的分类结果更加平滑，匀质域的分类错误更少，同一类别地物的分类结果也更加一致。参考图2(e)中的真实地物类别分布图，从图2(f)～图2(h)中可以看到本文RSC方法对几乎所有地物的分类精度都有所提高。

为了定量比较不同方法的分类性能，表2中列出了不同方法的分类性能。RSC方法的OA, F1-score和Purity值分别达到了0.8485, 0.8633和0.9047，都明显高于另外两种对比方法。与此同时，RSC方法的Entropy值为0.1344，是所有方法中最低的。可视化的分类结果图和定量评价指标均表明了本文所提方法的有效性。

6 结论

本文针对PolSAR图像的分类问题，基于黎曼流形充分挖掘极化协方差矩阵所包含的丰富信息，提出了基于极化协方差矩阵稀疏编码的PolSAR图像监督分类与非监督分类方法。在不同实验数据上的结果表明考虑黎曼流形结构特性能够有效地提高PolSAR图像的分类精度。下一步的工作是将极化协方差矩阵稀疏表达结合深度学习方法以进一步提高PolSAR图像的解译效果。