CN103839073A - 一种基于极化特征和近邻传播聚类的极化sar图像分类方法 - Google Patents

Abstract

本发明公开了一种基于极化特征和近邻传播聚类的极化SAR图像分类方法，主要解决现有无监督极化SAR分类方法分类精度不高的问题。其实现步骤为：对每个像素点进行特征分解和Freeman分解，提取像素点的三种系数和四种散射功率；根据获得的系数和散射功率对图像进行划分，得到7种类别；对获得的每一类，利用数据分布特征将其分成3个小类；对得到的21个类别计算相似度矩阵，然后用近邻传播聚类对其聚类，得到图像的预分类结果；最后，对预分类得到的图像用Wishart分类器进行迭代分类，得到最终分类结果。本发明与经典分类方法相比，对极化SAR图像的划分更加严谨，分类效果更佳，计算复杂度相对较小，可用于对极化SAR图像的进行地物分类和目标识别。

Description

一种基于极化特征和近邻传播聚类的极化SAR图像分类方法

技术领域

本发明属于图像处理技术领域，涉及极化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，eta1.Unsupervised classificationusing polarimetric decomposition and the complex Wishart classifier[J].IEEE Trans.Geosci.Remote Sensing.1999，37(5)：2249-2258.该方法是在原始H/α分类基础上增加了Wishart迭代，主要是对H/α划分后的8类利用Wishart分类器对每一个像素进行重新划分，从而有效的提高了分类的精度，但是又存在不能很好的保持各类的极化散射特性的不足。

近邻传播聚类是近年出现的一种在数据挖掘领域极具竞争力的聚类算法，相比较于传统聚类算法，近邻传播聚类方法能够在较短时间内完成大规模多类别数据集的聚类，并且该算法能够很好地解决非欧空间问题。然而该算法在应用到图像分割领域时，计算量和存储量是难以承受的，严重阻碍了该算法的性能发挥。

上述图像处理方法未能体现极化SAR数据的分布特性；也没有减少了近邻传播聚类算法的计算量和存储量；且不能根据图像的具体情况对图像自适应的进行聚类。

发明内容

本发明的目的在于针对已有技术的不足，提出一种基于极化特征和近邻传播聚类的极化SAR图像分类方法，可用于对极化SAR图像的地物分类和目标识别分类效果进行提高。

为实现上述目的，本发明技术方案是：一种基于极化特征和近邻传播聚类的极化SAR图像分类方法，其特征是：包括如下步骤：

(1)对待分类的极化SAR图像中每个像素点的相干矩阵T进行特征分解，得到每个像素点的单一散射机制系数ρ_s，双散射机制系数ρ_d和三种散射机制系数ρ_r；

(2)对待分类的极化SAR图像中每个像素点进行Freeman分解，得到每个像素点的表面散射功率P_s，二面角散射功率P_d，体散射功率P_v；

(3)根据每个像素点的单一散射机制系数ρ_s，双散射机制系数ρ_d和三种散射机制系数ρ_r，对待分类的极化SAR图像进行初始划分：根据max(ρ_s,ρ_d,ρ_v)的值的大小，将待分类的极化SAR图像数据初始划分为7类，即当max(ρ_s,ρ_d,ρ_v)＝ρ_s时，认为只有一种主散射机制，将max(P_s,P_d,P_v)＝P_s的对应像素点划分为P_s类，将max(P_s,P_d,P_v)＝P_d的对应像素点划分为P_d类，将max(P_s,P_d,P_v)＝P_v对应像素点划分为P_v类；当max(ρ_s,ρ_d,ρ_v)＝ρ_d时，认为同时存在两种主散射机制，将min(P_s,P_d,P_v)＝P_s的对应像素点划分为P_dP_v类，将min(P_s,P_d,P_v)＝P_d的对应像素点划分为P_sP_v类，将min(P_s,P_d,P_v)＝P_v对应像素点划分为P_sP_d类；当max(ρ_s,ρ_d,ρ_v)＝ρ_v时，认为三种散射同时存在，将对应的像素点划分为P_sP_dP_v类；

(4)对得到的7类初始划分结果，计算每一类中每个像素点的分布特征参数χ_L：

(4a)将每个像素点及其周围的像素点总共九个像素点作为一个小区域，计算该区域的相对峰值RK：

$\mathrm{RK}=\frac{1}{3}(\frac{E\left\{{\left|S_{\mathrm{HH}}\right|}^2\right\}}{E{\left\{\right|S_{\mathrm{HH}}\left|\right\}}^2}+\frac{E\left\{{\left|S_{\mathrm{HV}}\right|}^2\right\}}{E{\left\{\right|S_{\mathrm{HV}}\left|\right\}}^2}+\frac{E\left\{{\left|S_{\mathrm{VV}}\right|}^2\right\}}{E{\left\{\right|S_{\mathrm{VV}}\left|\right\}}^2}),$

其中，S_HH表示水平向发射和水平向接收的回波数据，S_VV表示垂直向发射和垂直向接收的回波数据，S_HV表示水平向发射垂直向接收的回波数据，|·|表示取这个数的模值，E{·}表示取这个数的均值；

(4b)根据相对峰值RK，计算分布特征参数χ_L：

${\mathrm{\χ}}_L=\frac{\mathrm{Ld}+1}{d+1}/(\mathrm{RK}-1),$

其中，L为极化SAR图像的视数，d为向量维数，取d=3；

(5)根据分布特征参数χ_L的值进一步将每一类划分结果划分为三类：如果χ_L＜25，将其对应的像素点划分为一类，如果25＜χ_L＜75，将其对应的像素点划分为一类，如果χ_L＞75，将其对应的像素点划分为一类，从而将整个极化SAR图像划分为21类别；

(6)利用划分后得到的21个类别，计算每一个小类别的平均聚类中心V_u：

$V_u=\frac{\underset{\mathrm{\κ}=1}{\overset{m_u}{\mathrm{\Σ}}}{T}_{\mathrm{\κ}}}{m_u},u=1,...,21,\mathrm{\κ}=\mathrm{1,2},...,m_u,$

其中T_κ是第u个小类别中各像素点的相干矩阵，m_u是属于第u个小类别中像素点的个数；

(7)对极化SAR图像每个像素点的相干矩阵T进行四分量分解，得到整幅图像的特征空间F；

(8)根据每一类得到的平均聚类中心V_u，构造每一类的相似度矩阵W，新相似度矩阵W中的每一个元素W(i,j)为：

$W(i,j)=\left\{\begin{array}{cc}0& i=j\\ \frac{1}{1+{\mathrm{\Δ}}^{\mathrm{SRW}}(i,j)\×{\mathrm{\Δ}}^F(i,j)\×d_{\mathrm{ij}}}& i\≠j\end{array}\right.,$

其中， ${\mathrm{\Δ}}^{\mathrm{SRW}}(i,j)=\left\{\begin{array}{cc}{d^2}_{\mathrm{SRW}}(V_i,V_j)& d_{\mathrm{SRW}}(V_i,V_j)\≤t_1\\ t_1^2+{2t}_1\×(d_{\mathrm{SRW}}(V_i,V_j)-t_1)& d_{\mathrm{SRW}}(V_i,V_j)>t_1\end{array}\right.$

${\mathrm{\Δ}}^F(i,j)=\left\{\begin{array}{cc}{d^2}_F(F_i,F_j)& d_F(F_i,F_j)\≤t_2\\ t_2^2+{2t}_2\×(d_F(F_i,F_j)-t_2)& d_F(F_i,F_j)>t_2\end{array}\right.$

d_ij是平均聚类中心V_i和V_j坐标的欧式距离，t₁t₂为常数，d_F(F_i,F_j)是特征向量F_i和F_j坐标的欧式距离， $d_{\mathrm{SRW}}(V_i,V_j)=\frac{1}{2}\left|\right\{\ln \left(\right|V_i\left|\right)+\ln \left(\right|V_j\left|\right)+\mathrm{Tr}({V_i}^{-1}{V}_j+{V_j}^{-1}{V}_i)\left\}\right|,$ Tr(·)是矩阵的迹，(·)^-1是矩阵的逆，i＝1,...,21,j＝1,...，21；

(9)根据相似度矩阵W，构造拉普拉斯矩阵L：

$L=D^{-\frac{1}{2}}{\mathrm{WD}}^{-\frac{1}{2}}$

其中D是对角矩阵，对角矩阵D上的每一个元素D_ij为：

$D_{\mathrm{ij}}=\left\{\begin{array}{cc}\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{W}_{\mathrm{ij}}& i=j\\ 0& i\≠j\end{array}\right.;$

(10)对拉普拉斯矩阵L进行特征值分解，得到前k个最大特征值对应的特征向量x₁,x₂,...,x_k，形成特征向量矩阵X＝[x₁,x₂,...,x_k]；

(11)用近邻传播聚类对特征向量矩阵X进行聚类，得到整幅图像的预分类结果；

(12)对预分类得到的整幅图像用K-Wishart分类器进行迭代分类，得到最终的分类结果。

上述发明方法中，步骤(1)所述的对待分类的极化SAR图像中每个像素点的相干矩阵T进行特征分解，按如下步骤进行：

(1a)读入图像的每个像素点，这些像素点为一个含有9个元素的3×3相干矩阵T；

$T=\frac{1}{2}\left[\begin{array}{ccc}\⟨{|S_{\mathrm{HH}}+S_{\mathrm{VV}}|}^2\⟩& \⟨(S_{\mathrm{HH}}+S_{\mathrm{VV}}){(S_{\mathrm{HH}}-S_{\mathrm{VV}})}^{*}\⟩& \⟨2(S_{\mathrm{HH}}+S_{\mathrm{VV}})S_{\mathrm{HV}}^{*}\⟩\\ \⟨(S_{\mathrm{HH}}-S_{\mathrm{VV}}){(S_{\mathrm{HH}}+S_{\mathrm{VV}})}^{*}\⟩& \⟨{|S_{\mathrm{HH}}-S_{\mathrm{VV}}|}^2\⟩& \⟨2(S_{\mathrm{HH}}-S_{\mathrm{VV}})S_{\mathrm{HV}}^{*}\⟩\\ \⟨2S_{\mathrm{HV}}{(S_{\mathrm{HH}}+S_{\mathrm{VV}})}^{*}\⟩& \⟨2S_{\mathrm{HV}}{(S_{\mathrm{HH}}-S_{\mathrm{VV}})}^{*}\⟩& \⟨4{\left|S_{\mathrm{HV}}\right|}^2\⟩\end{array}\right]$

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

(1b)将相干矩阵T进行特征值分解，由于相干矩阵T是一个3×3的矩阵，故分解后得到相干矩阵T的特征值为λ₁,λ₂,λ₃；

(1c)根据获得的特征值λ₁,λ₂,λ₃，计算像素点的单一散射机制系数ρ_s，双散射机制系数ρ_d和三种散射机制系数ρ_r：

ρ_s＝p₁-p₂

ρ_d＝2(p₂-p₃)，

ρ_r＝3p₃

其中， $p_i={\mathrm{\λ}}_i/\left(\underset{j=1}{\overset{3}{\mathrm{\Σ}}}{\mathrm{\λ}}_j\right),i=\mathrm{1,2,3}.$

其中：步骤(2)所述的对每个像素点进行Freeman分解，按如下步骤进行：

(2a)读入图像的每个像素点，该像素点为一个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{HV}}^{*}\⟩& \⟨S_{\mathrm{HH}}{S}_{\mathrm{VV}}\⟩\\ \sqrt{2}\⟨S_{\mathrm{HV}}{S}_{\mathrm{HH}}^{*}\⟩& 2\⟨{\left|S_{\mathrm{HV}}\right|}^2\⟩& \sqrt{2}\⟨S_{\mathrm{HV}}{S}_{\mathrm{VV}}^{*}\⟩\\ \⟨S_{\mathrm{VV}}{S}_{\mathrm{HH}}^{*}\⟩& \sqrt{2}\⟨S_{\mathrm{VV}}{S}_{\mathrm{HV}}^{*}\⟩& \⟨{\left|S_{\mathrm{VV}}\right|}^2\⟩\end{array}\right]---1)$

其中，U是中间变量， $U=\frac{1}{\sqrt{2}}\×\left(\begin{array}{ccc}1& 0& 1\\ 1& 0& -1\\ 0& \sqrt{2}& 0\end{array}\right),$ H表示水平极化，V表示垂直极化，S_HH表示水平向发射和水平向接收的回波数据，S_VV表示垂直向发射和垂直向接收的回波数据，S_HV表示水平向发射垂直向接收的回波数据，()^*表示这个数据的共轭，〈·〉表示按视数平均；

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

$C=F_s\left[\begin{array}{ccc}{\left|\mathrm{\β}\right|}^2& 0& \mathrm{\β}\\ 0& 0& 0\\ {\mathrm{\β}}^{*}& 0& 1\end{array}\right]+F_d\left[\begin{array}{ccc}{\left|\mathrm{\α}\right|}^2& 0& \mathrm{\α}\\ 0& 0& 0\\ {\mathrm{\α}}^{*}& 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]---2)$

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

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

$\left\{\begin{array}{c}\⟨{\left|S_{\mathrm{HH}}\right|}^2\⟩=F_s{\left|\mathrm{\β}\right|}^2+F_d{\left|\mathrm{\α}\right|}^2+F_v\\ \⟨{\left|S_{\mathrm{VV}}\right|}^2\⟩=F_s+F_d+F_v\\ \⟨S_{\mathrm{HH}}{S}_{\mathrm{VV}}^{*}\⟩=F_s\mathrm{\β}+F_d\mathrm{\α}+F_v/3\\ \⟨{\left|S_{\mathrm{HV}}\right|}^2\⟩=F_v/3\end{array}\right.---3);$

(2d)计算像素点协方差矩阵C中的

的值，如果

则令α＝-1，如果

则令β＝1，给定α或β的值后，剩余的4个未知数则根据式3）求解，其中Re(·)表示取实部；

(2e)根据求解出的F_s,F_v,F_d,α,β，求解出体散射功率P_v、二面角散射功率P_d和表面散射功率P_s：

$\begin{array}{c}{P}_v=\frac{8F_v}{3}\\ P_d=F_d(1+{\left|\mathrm{\α}\right|}^2).\\ P_s=F_s(1+{\left|\mathrm{\β}\right|}^2)\end{array}---4)$

其中：步骤(11)所述的用近邻传播聚类对特征向量矩阵X进行聚类，得到整幅图像的预分类结果，按如下步骤进行：

(11a)根据特征向量矩阵X计算新相似矩阵S，

S(i,j)＝-||x_i-x_j||²

(11b)将新相似矩阵S的对角线上的值都设为I，I取值为-1000，将归属度矩阵AS和吸引度矩阵LS初始值设置为0；

(11b)将归属度矩阵AS和吸引度矩阵LS分别赋值给中间变量a和b，即a＝AS，b＝LS；

(11c)根据新相似度矩阵S和归属度矩阵AS，计算吸引度矩阵LS，吸引度矩阵LS中的每一个元素LS_iu：

${\mathrm{LS}}_{\mathrm{iu}}=S_{\mathrm{iu}}-\underset{\mathrm{\φ}\≠u}{\max }(S_{\mathrm{i\φ}}+{\mathrm{AS}}_{\mathrm{i\φ}}),$

其中，i＝1,...,21，u＝1,...,21，φ＝1,...,21；

(11d)根据新相似度矩阵S和吸引度矩阵LS，计算归属度矩阵AS，归属度矩阵AS中每一个元素AS_iu为；

${\mathrm{AS}}_{\mathrm{iu}}=\left\{\begin{array}{c}\underset{\mathrm{\φ}\≠u}{\mathrm{\Σ}}\max (0,{\mathrm{LS}}_{\mathrm{\φu}}),i=u\\ \min \{0,{\mathrm{LS}}_{\mathrm{uu}}+\underset{\mathrm{\φ}\≠i,u}{\mathrm{\Σ}}\max (0,{\mathrm{LS}}_{\mathrm{\φu}})\},i\≠u\end{array}\right.,$

其中，i＝1,...,21，u＝1,...,21，φ＝1,...,21；

(11e)更新归属度矩阵AS和吸引度矩阵LS：根据已获得的归属度矩阵AS和吸引度矩阵LS，则中间变量g＝(1-λ)×AS+λ×a，中间变量f＝(1-λ)×LS+λ×b，再令AS＝g，LS＝f，更新归属度矩阵AS和吸引度矩阵LS，其中，λ是收敛系数，λ∈[0.5,1)；

(11g)重复步骤(11b)-(11e)，直到迭代次数等于给定的迭代次数n，其中n＝1000，根据LS(u,u)+AS(u,u)的值判断u是否为聚类中心，如果LS(u,u)+AS(u,u)＞0，则u是聚类中心，反之，u不是聚类中心，u＝1,...,21，从而获得预分类结果。

其中：步骤(12)所述的对预分类得到的整幅图像用Wishart分类器进行迭代分类，按如下步骤进行：

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

$B_c=\frac{\underset{\mathrm{\ρ}=1}{\overset{n_c}{\mathrm{\Σ}}}{T}_{\mathrm{\ρ}}}{n_c},c=1,...,k,\mathrm{\ρ}=\mathrm{1,2},...,n_c,$

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

(12b)根据下式计算每个像素点的分布特征参数χ_L：

${\mathrm{\χ}}_L=\frac{\mathrm{qL}+1}{q+1}/(\mathrm{RK}-1),$

其中，L为极化SAR图像的视数，q为向量维数，取q=3，RK是每个像素点及其周围的像素点总共九个像素点组成的区域的相对峰值，RK定义如下：

$\mathrm{RK}=\frac{1}{3}(\frac{E\left\{{\left|S_{\mathrm{HH}}\right|}^2\right\}}{E{\left\{\right|S_{\mathrm{HH}}\left|\right\}}^2}+\frac{E\left\{{\left|S_{\mathrm{HV}}\right|}^2\right\}}{E{\left\{\right|S_{\mathrm{HV}}\left|\right\}}^2}+\frac{E\left\{{\left|S_{\mathrm{VV}}\right|}^2\right\}}{E{\left\{\right|S_{\mathrm{VV}}\left|\right\}}^2}),$

其中，S_HH表示水平向发射和水平向接收的回波数据，S_VV表示垂直向发射和垂直向接收的回波数据，S_HV表示水平向发射垂直向接收的回波数据，|·|表示取这个数的模值，E{·}表示取这个数的均值；

(12c)根据每个像素点的分布特征参数χ_L，计算每个像素点到第i类聚类中心V_i的距离d(〈C〉,V_i)：

$d(\⟨C\⟩,V_i)=\mathrm{Lin}\left[V_i\right]+\ln \left(\mathrm{\Γ}\left({\mathrm{\χ}}_L\right)\right)-\frac{{\mathrm{\χ}}_L+\mathrm{qL}}{2}\ln \left({\mathrm{L\χ}}_L\right)-\frac{{\mathrm{\χ}}_L-\mathrm{qL}}{2}\ln \mathrm{Tr}({V_i}^{-1}\⟨C\⟩)-\ln K_{{\mathrm{\χ}}_L-\mathrm{qL}}\left(2\sqrt{L{\mathrm{\χ}}_L\mathrm{Tr}({V_i}^{-1}\⟨C\⟩)}\right)$

其中C是像素点的协方差矩阵，〈·〉表示按视数平均，[·]表示矩阵的行列式，Tr(·)表示矩阵的迹，V_i ^-1表示对聚类中心V_i求逆，Γ(·)表示Gamma函数，

表示(χ_L-qL)阶的第二类修正贝塞尔函数，L为图像的视数，q为向量维数，取q=3；

(12e)根据每个像素点到第i类聚类中心的距离对极化SAR图像进行重新划分：如果d(〈C〉,B_τ)≤d(〈C〉,B_ψ)，则将该像素点划分为第τ类，如果d(〈C〉,B_τ)＞d(〈C〉,B_ψ)，则将该像素点划分为第ψ类，其中d(〈C〉,B_τ)是像素点到第τ类聚类中心的距离，d(〈C〉,B_ψ)是像素点到第ψ类聚类中心的距离，τ,ψ＝1,...,8,τ≠ψ；

(12f)重复步骤(12a)-(12e)直到迭代次数等于给定的迭代次数μ，其中μ＝3，得到更为准确的分类结果。

本发明技术效果和优点如下：本发明采用对每个像素点进行特征分解和Freeman分解，提取像素点的三种系数和三种散射功率；根据获得的系数和散射功率对图像进行划分，得到7种类别；对获得的每一类，利用数据分布特征将其分成3个小类；对得到的21个类别计算相似度矩阵，然后用近邻传播聚类对其聚类，得到图像的预分类结果；最后，对预分类得到的图像用Wishart分类器进行迭代分类，得到最终分类结果；并且具有以下优点：

1.本发明用改进的Wishart距离代替传统的欧式距离计算像素点的相似性，更能体现极化SAR数据的分布特性；

2.本发明先对极化SAR图像划分成7类，将划分得到的区域，作为近邻传播聚类的输入数据点，减少了近邻传播聚类算法的计算量和存储量；

3.本发明由于采用了近邻传播聚类算法，因此可根据图像的具体情况对图像自适应的进行聚类。

4.对极化SAR图像的划分更加严谨，分类效果更佳，计算复杂度相对较小，可用于对极化SAR图像的进行地物分类和目标识别。

附图说明

下面结合实施例附图对本发明的具体实施方式予以详细描述。

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

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

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

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

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

具体实施方式

实施例1，参照附图1描述，具体实现步骤如下：

一种基于极化特征和近邻传播聚类的极化SAR图像分类方法，包括如下步骤：

(1)对待分类的极化SAR图像中每个像素点的相干矩阵T进行特征分解，得到每个像素点的单一散射机制系数ρ_s，双散射机制系数ρ_d和三种散射机制系数ρ_r；

(2)对待分类的极化SAR图像中每个像素点进行Freeman分解，得到每个像素点的表面散射功率P_s，二面角散射功率P_d，体散射功率P_v；

(3)根据每个像素点的单一散射机制系数ρ_s，双散射机制系数ρ_d和三种散射机制系数ρ_r，对待分类的极化SAR图像进行初始划分：根据max(ρ_s,ρ_d,ρ_v)的值的大小，将待分类的极化SAR图像数据初始划分为7类，即当max(ρ_s,ρ_d,ρ_v)＝ρ_s时，认为只有一种主散射机制，将max(P_s,P_d,P_v)＝P_s的对应像素点划分为P_s类，将max(P_s,P_d,P_v)＝P_d的对应像素点划分为P_d类，将max(P_s,P_d,P_v)＝P_v对应像素点划分为P_v类；当max(ρ_s,ρ_d,ρ_v)＝ρ_d时，认为同时存在两种主散射机制，将min(P_s,P_d,P_v)＝P_s的对应像素点划分为P_dP_v类，将min(P_s,P_d,P_v)＝P_d的对应像素点划分为P_sP_v类，将min(P_s,P_d,P_v)＝P_v对应像素点划分为P_sP_d类；当max(ρ_s,ρ_d,ρ_v)＝ρ_v时，认为三种散射同时存在，将对应的像素点划分为P_sP_dP_v类；

(4)对得到的7类初始划分结果，计算每一类中每个像素点的分布特征参数χ_L：

(4a)将每个像素点及其周围的像素点总共九个像素点作为一个小区域，计算该区域的相对峰值RK：

$\mathrm{RK}=\frac{1}{3}(\frac{E\left\{{\left|S_{\mathrm{HH}}\right|}^2\right\}}{E{\left\{\right|S_{\mathrm{HH}}\left|\right\}}^2}+\frac{E\left\{{\left|S_{\mathrm{HV}}\right|}^2\right\}}{E{\left\{\right|S_{\mathrm{HV}}\left|\right\}}^2}+\frac{E\left\{{\left|S_{\mathrm{VV}}\right|}^2\right\}}{E{\left\{\right|S_{\mathrm{VV}}\left|\right\}}^2}),$

其中，S_HH表示水平向发射和水平向接收的回波数据，S_VV表示垂直向发射和垂直向接收的回波数据，S_HV表示水平向发射垂直向接收的回波数据，|·|表示取这个数的模值，E{·}表示取这个数的均值；

(4b)根据相对峰值RK，计算分布特征参数χ_L：

${\mathrm{\χ}}_L=\frac{\mathrm{Ld}+1}{d+1}/(\mathrm{RK}-1),$

其中，L为极化SAR图像的视数，d为向量维数，取d=3；

(5)根据分布特征参数χ_L的值进一步将每一类划分结果划分为三类：如果χ_L＜25，将其对应的像素点划分为一类，如果25＜χ_L＜75，将其对应的像素点划分为一类，如果χ_L＞75，将其对应的像素点划分为一类，从而将整个极化SAR图像划分为21类别。

(6)利用划分后得到的21个类别，计算每一个小类别的平均聚类中心V_u：

$V_u=\frac{\underset{\mathrm{\κ}=1}{\overset{m_u}{\mathrm{\Σ}}}{T}_{\mathrm{\κ}}}{m_u},u=1,...,21,\mathrm{\κ}=\mathrm{1,2},...,m_u,$

其中T_κ是第u个小类别中各像素点的相干矩阵，m_u是属于第u个小类别中像素点的个数；

(7)对极化SAR图像每个像素点的相干矩阵T进行四分量分解，得到整幅图像的特征空间F；

(8)根据每一类得到的平均聚类中心V_u，构造每一类的相似度矩阵W，新相似度矩阵W中的每一个元素W(i,j)为：

$W(i,j)=\left\{\begin{array}{cc}0& i=j\\ \frac{1}{1+{\mathrm{\Δ}}^{\mathrm{SRW}}(i,j)\×{\mathrm{\Δ}}^F(i,j)\×d_{\mathrm{ij}}}& i\≠j\end{array}\right.,$

其中， ${\mathrm{\Δ}}^{\mathrm{SRW}}(i,j)=\left\{\begin{array}{cc}{d^2}_{\mathrm{SRW}}(V_i,V_j)& d_{\mathrm{SRW}}(V_i,V_j)\≤t_1\\ t_1^2+{2t}_1\×(d_{\mathrm{SRW}}(V_i,V_j)-t_1)& d_{\mathrm{SRW}}(V_i,V_j)>t_1\end{array}\right.$

${\mathrm{\Δ}}^F(i,j)=\left\{\begin{array}{cc}{d^2}_F(F_i,F_j)& d_F(F_i,F_j)\≤t_2\\ t_2^2+{2t}_2\×(d_F(F_i,F_j)-t_2)& d_F(F_i,F_j)>t_2\end{array}\right.$

d_ij是平均聚类中心V_i和V_j坐标的欧式距离，t₁t₂为常数，d_F(F_i,F_j)是特征向量F_i和F_j坐标的欧式距离， $d_{\mathrm{SRW}}(V_i,V_j)=\frac{1}{2}\left|\right\{\ln \left(\right|V_i\left|\right)+\ln \left(\right|V_j\left|\right)+\mathrm{Tr}({V_i}^{-1}{V}_j+{V_j}^{-1}{V}_i)\left\}\right|,$ Tr(·)是矩阵的迹，(·)^-1是矩阵的逆，i＝1,...,21,j＝1,...,21；

(9)根据相似度矩阵W，构造拉普拉斯矩阵L：

$L=D^{-\frac{1}{2}}{\mathrm{WD}}^{-\frac{1}{2}}$

其中D是对角矩阵，对角矩阵D上的每一个元素D_ij为：

$D_{\mathrm{ij}}=\left\{\begin{array}{cc}\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{W}_{\mathrm{ij}}& i=j\\ 0& i\≠j\end{array}\right.;$

(10)对拉普拉斯矩阵L进行特征值分解，得到前k个最大特征值对应的特征向量x₁,x₂,...,x_k，形成特征向量矩阵X＝[x₁,x₂,...,x_k]；

(11)用近邻传播聚类对特征向量矩阵X进行聚类，得到整幅图像的预分类结果；

(12)对预分类得到的整幅图像用K-Wishart分类器进行迭代分类，得到最终的分类结果。

上述图像处理方法体现了极化SAR数据的分布特性；也减少了近邻传播聚类算法的计算量和存储量；且能根据图像的具体情况对图像自适应的进行聚类。

实施例2，结合附图1至5描述。

在实施例1的基础上，对该极化SAR图像分类方法进行细化，其中：

步骤(1)所述的对待分类的极化SAR图像中每个像素点的相干矩阵T进行特征分解，按如下步骤进行：

(1a)读入图像的每个像素点，这些像素点为一个含有9个元素的3×3相干矩阵T；

$T=\frac{1}{2}\left[\begin{array}{ccc}\⟨{|S_{\mathrm{HH}}+S_{\mathrm{VV}}|}^2\⟩& \⟨(S_{\mathrm{HH}}+S_{\mathrm{VV}}){(S_{\mathrm{HH}}-S_{\mathrm{VV}})}^{*}\⟩& \⟨2(S_{\mathrm{HH}}+S_{\mathrm{VV}})S_{\mathrm{HV}}^{*}\⟩\\ \⟨(S_{\mathrm{HH}}-S_{\mathrm{VV}}){(S_{\mathrm{HH}}+S_{\mathrm{VV}})}^{*}\⟩& \⟨{|S_{\mathrm{HH}}-S_{\mathrm{VV}}|}^2\⟩& \⟨2(S_{\mathrm{HH}}-S_{\mathrm{VV}})S_{\mathrm{HV}}^{*}\⟩\\ \⟨2S_{\mathrm{HV}}{(S_{\mathrm{HH}}+S_{\mathrm{VV}})}^{*}\⟩& \⟨2S_{\mathrm{HV}}{(S_{\mathrm{HH}}-S_{\mathrm{VV}})}^{*}\⟩& \⟨4{\left|S_{\mathrm{HV}}\right|}^2\⟩\end{array}\right]$

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

(1b)将相干矩阵T进行特征值分解，由于相干矩阵T是一个3×3的矩阵，故分解后得到相干矩阵T的特征值为λ₁,λ₂,λ₃；

(1c)根据获得的特征值λ₁,λ₂,λ₃，计算像素点的单一散射机制系数ρ_s，双散射机制系数ρ_d和三种散射机制系数ρ_r：

ρ_s＝p₁-p₂

ρ_d＝2(p₂-p₃)，

ρ_r＝3p₃

其中， $p_i={\mathrm{\λ}}_i/\left(\underset{j=1}{\overset{3}{\mathrm{\Σ}}}{\mathrm{\λ}}_j\right),i=\mathrm{1,2,3}.$

其中：步骤(2)所述的对每个像素点进行Freeman分解，按如下步骤进行：

(2a)读入图像的每个像素点，该像素点为一个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{HV}}^{*}\⟩& \⟨S_{\mathrm{HH}}{S}_{\mathrm{VV}}\⟩\\ \sqrt{2}\⟨S_{\mathrm{HV}}{S}_{\mathrm{HH}}^{*}\⟩& 2\⟨{\left|S_{\mathrm{HV}}\right|}^2\⟩& \sqrt{2}\⟨S_{\mathrm{HV}}{S}_{\mathrm{VV}}^{*}\⟩\\ \⟨S_{\mathrm{VV}}{S}_{\mathrm{HH}}^{*}\⟩& \sqrt{2}\⟨S_{\mathrm{VV}}{S}_{\mathrm{HV}}^{*}\⟩& \⟨{\left|S_{\mathrm{VV}}\right|}^2\⟩\end{array}\right]---1)$

其中，U是中间变量， $U=\frac{1}{\sqrt{2}}\×\left(\begin{array}{ccc}1& 0& 1\\ 1& 0& -1\\ 0& \sqrt{2}& 0\end{array}\right),$ H表示水平极化，V表示垂直极化，S_HH表示水平向发射和水平向接收的回波数据，S_VV表示垂直向发射和垂直向接收的回波数据，S_HV表示水平向发射垂直向接收的回波数据，()^*表示这个数据的共轭，〈·〉表示按视数平均；

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

$C=F_s\left[\begin{array}{ccc}{\left|\mathrm{\β}\right|}^2& 0& \mathrm{\β}\\ 0& 0& 0\\ {\mathrm{\β}}^{*}& 0& 1\end{array}\right]+F_d\left[\begin{array}{ccc}{\left|\mathrm{\α}\right|}^2& 0& \mathrm{\α}\\ 0& 0& 0\\ {\mathrm{\α}}^{*}& 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]---2)$

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

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

$\left\{\begin{array}{c}\⟨{\left|S_{\mathrm{HH}}\right|}^2\⟩=F_s{\left|\mathrm{\β}\right|}^2+F_d{\left|\mathrm{\α}\right|}^2+F_v\\ \⟨{\left|S_{\mathrm{VV}}\right|}^2\⟩=F_s+F_d+F_v\\ \⟨S_{\mathrm{HH}}{S}_{\mathrm{VV}}^{*}\⟩=F_s\mathrm{\β}+F_d\mathrm{\α}+F_v/3\\ \⟨{\left|S_{\mathrm{HV}}\right|}^2\⟩=F_v/3\end{array}\right.---3);$

(2d)计算像素点协方差矩阵C中的

的值，如果

则令α＝-1，如果

则令β＝1，给定α或β的值后，剩余的4个未知数则根据式3）求解，其中Re(·)表示取实部；

(2e)根据求解出的F_s,F_v,F_d,α,β，求解出体散射功率P_v、二面角散射功率P_d和表面散射功率P_s：

$\begin{array}{c}{P}_v=\frac{8F_v}{3}\\ P_d=F_d(1+{\left|\mathrm{\α}\right|}^2).\\ P_s=F_s(1+{\left|\mathrm{\β}\right|}^2)\end{array}---4)$

其中：步骤(11)所述的用近邻传播聚类对特征向量矩阵X进行聚类，得到整幅图像的预分类结果，按如下步骤进行：

(11a)根据特征向量矩阵X计算新相似矩阵S，

S(i,j)＝-||x_i-x_j||²

(11b)将新相似矩阵S的对角线上的值都设为I，I取值为-1000，将归属度矩阵AS和吸引度矩阵LS初始值设置为0；

(11b)将归属度矩阵AS和吸引度矩阵LS分别赋值给中间变量a和b，即a＝AS，b＝LS；

(11c)根据新相似度矩阵S和归属度矩阵AS，计算吸引度矩阵LS，吸引度矩阵LS中的每一个元素LS_iu：

${\mathrm{LS}}_{\mathrm{iu}}=S_{\mathrm{iu}}-\underset{\mathrm{\φ}\≠u}{\max }(S_{\mathrm{i\φ}}+{\mathrm{AS}}_{\mathrm{i\φ}}),$

其中，i＝1,...,21，u＝1,...,21，φ＝1,...,21；

(11d)根据新相似度矩阵S和吸引度矩阵LS，计算归属度矩阵AS，归属度矩阵AS中每一个元素AS_iu为；

${\mathrm{AS}}_{\mathrm{iu}}=\left\{\begin{array}{c}\underset{\mathrm{\φ}\≠u}{\mathrm{\Σ}}\max (0,{\mathrm{LS}}_{\mathrm{\φu}}),i=u\\ \min \{0,{\mathrm{LS}}_{\mathrm{uu}}+\underset{\mathrm{\φ}\≠i,u}{\mathrm{\Σ}}\max (0,{\mathrm{LS}}_{\mathrm{\φu}})\},i\≠u\end{array}\right.,$

其中，i＝1,...,21，u＝1,...,21，φ＝1,...,21；

(11e)更新归属度矩阵AS和吸引度矩阵LS：根据已获得的归属度矩阵AS和吸引度矩阵LS，则中间变量g＝(1-λ)×AS+λ×a，中间变量f＝(1-λ)×LS+λ×b，再令AS＝g，LS＝f，更新归属度矩阵AS和吸引度矩阵LS，其中，λ是收敛系数，λ∈[0.5,1)；

(11g)重复步骤(11b)-(11e)，直到迭代次数等于给定的迭代次数n，其中n＝1000，根据LS(u,u)+AS(u,u)的值判断u是否为聚类中心，如果LS(u,u)+AS(u,u)＞0，则u是聚类中心，反之，u不是聚类中心，u＝1,...,21，从而获得预分类结果。

其中：步骤(12)所述的对预分类得到的整幅图像用Wishart分类器进行迭代分类，按如下步骤进行：

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

$B_c=\frac{\underset{\mathrm{\ρ}=1}{\overset{n_c}{\mathrm{\Σ}}}{T}_{\mathrm{\ρ}}}{n_c},c=1,...,k,\mathrm{\ρ}=\mathrm{1,2},...,n_c,$

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

(12b)根据下式计算每个像素点的分布特征参数χ_L：

${\mathrm{\χ}}_L=\frac{\mathrm{qL}+1}{q+1}/(\mathrm{RK}-1),$

其中，L为极化SAR图像的视数，q为向量维数，取q=3，RK是每个像素点及其周围的像素点总共九个像素点组成的区域的相对峰值，RK定义如下：

$\mathrm{RK}=\frac{1}{3}(\frac{E\left\{{\left|S_{\mathrm{HH}}\right|}^2\right\}}{E{\left\{\right|S_{\mathrm{HH}}\left|\right\}}^2}+\frac{E\left\{{\left|S_{\mathrm{HV}}\right|}^2\right\}}{E{\left\{\right|S_{\mathrm{HV}}\left|\right\}}^2}+\frac{E\left\{{\left|S_{\mathrm{VV}}\right|}^2\right\}}{E{\left\{\right|S_{\mathrm{VV}}\left|\right\}}^2}),$

其中，S_HH表示水平向发射和水平向接收的回波数据，S_VV表示垂直向发射和垂直向接收的回波数据，S_HV表示水平向发射垂直向接收的回波数据，|·|表示取这个数的模值，E{·}表示取这个数的均值；

(12c)根据每个像素点的分布特征参数χ_L，计算每个像素点到第i类聚类中心V_i的距离d(〈C〉,V_i)：

$d(\⟨C\⟩,V_i)=\mathrm{Lin}\left[V_i\right]+\ln \left(\mathrm{\Γ}\left({\mathrm{\χ}}_L\right)\right)-\frac{{\mathrm{\χ}}_L+\mathrm{qL}}{2}\ln \left({\mathrm{L\χ}}_L\right)-\frac{{\mathrm{\χ}}_L-\mathrm{qL}}{2}\ln \mathrm{Tr}({V_i}^{-1}\⟨C\⟩)-\ln K_{{\mathrm{\χ}}_L-\mathrm{qL}}\left(2\sqrt{L{\mathrm{\χ}}_L\mathrm{Tr}({V_i}^{-1}\⟨C\⟩)}\right)$

其中C是像素点的协方差矩阵，〈·〉表示按视数平均，[·]表示矩阵的行列式，Tr(·)表示矩阵的迹，V_i ^-1表示对聚类中心V_i求逆，Γ(·)表示Gamma函数，表示(χ_L-qL)阶的第二类修正贝塞尔函数，L为图像的视数，q为向量维数，取q=3；

(12e)根据每个像素点到第i类聚类中心的距离对极化SAR图像进行重新划分：如果d(〈C〉,B_τ)≤d(〈C〉,B_ψ)，则将该像素点划分为第τ类，如果d(〈C〉,B_τ)＞d(〈C〉,B_ψ)，则将该像素点划分为第ψ类，其中d(〈C〉,B_τ)是像素点到第τ类聚类中心的距离，d(〈C〉,B_ψ)是像素点到第ψ类聚类中心的距离，τ,ψ＝1,...,8,τ≠ψ；

(12f)重复步骤(12a)-(12e)直到迭代次数等于给定的迭代次数μ，其中μ=3，得到更为准确的分类结果。

该实施例中对极化SAR图像的划分更加严谨，分类效果更佳，计算复杂度相对较小，可用于对极化SAR图像的进行地物分类和目标识别。凡实施例1中未展开描述的步骤均属于现有技术或公知常识，通常通过网络便可得知。

本发明的极化SAR图像分类方法还可以通过下面的实验仿真进一步说明处理效果：

1、实验条件与方法

硬件平台为：IntelCore2DuoCPUE65502.33GHZ、2GBRAM；

软件平台为：MATLABR2011a；

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

2、仿真内容与结果

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

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

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

仿真二，用现有的基H/α/wishart分类方法对图2进行分类，分类结果见图4。

由图4可见，分类效果与图3相比，类别区分更合理，类别区分也较为合理，但由于类别数的限制，还有较多区域划分不清楚，且分类是对于区域边界的也过于武断。

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

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

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

以上例举仅仅是对本发明的举例说明，并不构成对本发明保护范围的限制，凡是与本发明相同或相似的设计均属于本发明的保护范围之内。

Claims (5)

1.一种基于极化特征和近邻传播聚类的极化SAR图像分类方法，其特征是：包括如下步骤：

(1)对待分类的极化SAR图像中每个像素点的相干矩阵T进行特征分解，得到每个像素点的单一散射机制系数ρ_s，双散射机制系数ρ_d和三种散射机制系数ρ_r；

(2)对待分类的极化SAR图像中每个像素点进行Freeman分解，得到每个像素点的表面散射功率P_s，二面角散射功率P_d，体散射功率P_v；

(3)根据每个像素点的单一散射机制系数ρ_s，双散射机制系数ρ_d和三种散射机制系数ρ_r，对待分类的极化SAR图像进行初始划分：根据max(ρ_s,ρ_d,ρ_v)的值的大小，将待分类的极化SAR图像数据初始划分为7类，即当max(ρ_s,ρ_d,ρ_v)＝ρ_s时，认为只有一种主散射机制，将max(P_s,P_d,P_v)＝P_s的对应像素点划分为P_s类，将max(P_s,P_d,P_v)＝P_d的对应像素点划分为P_d类，将max(P_s,P_d,P_v)＝P_v对应像素点划分为P_v类；当max(ρ_s,ρ_d,ρ_v)＝ρ_d时，认为同时存在两种主散射机制，将min(P_s,P_d,P_v)＝P_s的对应像素点划分为P_dP_v类，将min(P_s,P_d,P_v)＝P_d的对应像素点划分为P_sP_v类，将min(P_s,P_d,P_v)＝P_v对应像素点划分为P_sP_d类；当max(ρ_s,ρ_d,ρ_v)＝ρ_v时，认为三种散射同时存在，将对应的像素点划分为P_sP_dP_v类；

(4)对得到的7类初始划分结果，计算每一类中每个像素点的分布特征参数χ_L：

(4a)将每个像素点及其周围的像素点总共九个像素点作为一个小区域，计算该区域的相对峰值RK：

$\mathrm{RK}=\frac{1}{3}(\frac{E\left\{{\left|S_{\mathrm{HH}}\right|}^2\right\}}{E{\left\{\right|S_{\mathrm{HH}}\left|\right\}}^2}+\frac{E\left\{{\left|S_{\mathrm{HV}}\right|}^2\right\}}{E{\left\{\right|S_{\mathrm{HV}}\left|\right\}}^2}+\frac{E\left\{{\left|S_{\mathrm{VV}}\right|}^2\right\}}{E{\left\{\right|S_{\mathrm{VV}}\left|\right\}}^2}),$

其中，S_HH表示水平向发射和水平向接收的回波数据，S_VV表示垂直向发射和垂直向接收的回波数据，S_HV表示水平向发射垂直向接收的回波数据，|·|表示取这个数的模值，E{·}表示取这个数的均值；

(4b)根据相对峰值RK，计算分布特征参数χ_L：

${\mathrm{\χ}}_L=\frac{\mathrm{Ld}+1}{d+1}/(\mathrm{RK}-1),$

其中，L为极化SAR图像的视数，d为向量维数，取d=3；

(5)根据分布特征参数χ_L的值进一步将每一类划分结果划分为三类：如果χ_L＜25，将其对应的像素点划分为一类，如果25＜χ_L＜75，将其对应的像素点划分为一类，如果χ_L＞75，将其对应的像素点划分为一类，从而将整个极化SAR图像划分为21类别。

(6)利用划分后得到的21个类别，计算每一个小类别的平均聚类中心V_u：

$V_u=\frac{\underset{\mathrm{\κ}=1}{\overset{m_u}{\mathrm{\Σ}}}{T}_{\mathrm{\κ}}}{m_u},u=1,...,21,\mathrm{\κ}=\mathrm{1,2},...,m_u,$

其中T_κ是第u个小类别中各像素点的相干矩阵，m_u是属于第u个小类别中像素点的个数；

(7)对极化SAR图像每个像素点的相干矩阵T进行四分量分解，得到整幅图像的特征空间F；

(8)根据每一类得到的平均聚类中心V_u，构造每一类的相似度矩阵W，相似度矩阵W中的每一个元素W(i,j)为：

$W(i,j)=\left\{\begin{array}{cc}0& i=j\\ \frac{1}{1+{\mathrm{\Δ}}^{\mathrm{SRW}}(i,j)\×{\mathrm{\Δ}}^F(i,j)\×d_{\mathrm{ij}}}& i\≠j\end{array}\right.,$

其中， ${\mathrm{\Δ}}^{\mathrm{SRW}}(i,j)=\left\{\begin{array}{cc}{d^2}_{\mathrm{SRW}}(V_i,V_j)& d_{\mathrm{SRW}}(V_i,V_j)\≤t_1\\ t_1^2+{2t}_1\×(d_{\mathrm{SRW}}(V_i,V_j)-t_1)& d_{\mathrm{SRW}}(V_i,V_j)>t_1\end{array}\right.$

${\mathrm{\Δ}}^F(i,j)=\left\{\begin{array}{cc}{d^2}_F(F_i,F_j)& d_F(F_i,F_j)\≤t_2\\ t_2^2+{2t}_2\×(d_F(F_i,F_j)-t_2)& d_F(F_i,F_j)>t_2\end{array}\right.$

d_ij是平均聚类中心V_i和V_j坐标的欧式距离，t₁t₂为常数，d_F(F_i,F_j)是特征向量F_i和F_j坐标的欧式距离， $d_{\mathrm{SRW}}(V_i,V_j)=\frac{1}{2}\left|\right\{\ln \left(\right|V_i\left|\right)+\ln \left(\right|V_j\left|\right)+\mathrm{Tr}({V_i}^{-1}{V}_j+{V_j}^{-1}{V}_i)\left\}\right|,$ Tr(·)是矩阵的迹，(·)^-1是矩阵的逆，i＝1,...,21,j＝1,...,21；

(9)根据相似度矩阵W，构造拉普拉斯矩阵L：

$L=D^{-\frac{1}{2}}{\mathrm{WD}}^{-\frac{1}{2}}$

其中D是对角矩阵，对角矩阵D上的每一个元素D_ij为：

$D_{\mathrm{ij}}=\left\{\begin{array}{cc}\underset{j=1}{\overset{M}{\mathrm{\Σ}}}{W}_{\mathrm{ij}}& i=j\\ 0& i\≠j\end{array}\right.;$

(10)对拉普拉斯矩阵L进行特征值分解，得到前k个最大特征值对应的特征向量x₁,x₂,...,x_k，形成特征向量矩阵X＝[x₁,x₂,...,x_k]；

(11)用近邻传播聚类对特征向量矩阵X进行聚类，得到整幅图像的预分类结果；

(12)对预分类得到的整幅图像用K-Wishart分类器进行迭代分类，得到最终的分类结果。

2.根据权利要求1所述的基于极化特征和近邻传播聚类的极化SAR图像分类方法，其特征是和：其中步骤(1)所述的对待分类的极化SAR图像中每个像素点的相干矩阵T进行特征分解，按如下步骤进行：

(1a)读入图像的每个像素点，这些像素点为一个含有9个元素的3×3相干矩阵T；

$T=\frac{1}{2}\left[\begin{array}{ccc}\⟨{|S_{\mathrm{HH}}+S_{\mathrm{VV}}|}^2\⟩& \⟨(S_{\mathrm{HH}}+S_{\mathrm{VV}}){(S_{\mathrm{HH}}-S_{\mathrm{VV}})}^{*}\⟩& \⟨2(S_{\mathrm{HH}}+S_{\mathrm{VV}})S_{\mathrm{HV}}^{*}\⟩\\ \⟨(S_{\mathrm{HH}}-S_{\mathrm{VV}}){(S_{\mathrm{HH}}+S_{\mathrm{VV}})}^{*}\⟩& \⟨{|S_{\mathrm{HH}}-S_{\mathrm{VV}}|}^2\⟩& \⟨2(S_{\mathrm{HH}}-S_{\mathrm{VV}})S_{\mathrm{HV}}^{*}\⟩\\ \⟨2S_{\mathrm{HV}}{(S_{\mathrm{HH}}+S_{\mathrm{VV}})}^{*}\⟩& \⟨2S_{\mathrm{HV}}{(S_{\mathrm{HH}}-S_{\mathrm{VV}})}^{*}\⟩& \⟨4{\left|S_{\mathrm{HV}}\right|}^2\⟩\end{array}\right]$

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

(1b)将相干矩阵T进行特征值分解，由于相干矩阵T是一个3×3的矩阵，故分解后得到相干矩阵T的特征值为λ₁,λ₂,λ₃；

(1c)根据获得的特征值λ₁,λ₂,λ₃，计算像素点的单一散射机制系数ρ_s，双散射机制系数ρ_d和三种散射机制系数ρ_r：

ρ_s＝p₁-p₂

ρ_d＝2(p₂-p₃)，

ρ_r＝3p₃

其中， $p_i={\mathrm{\λ}}_i/\left(\underset{j=1}{\overset{3}{\mathrm{\Σ}}}{\mathrm{\λ}}_j\right),i=\mathrm{1,2,3}.$

3.根据权利要求1所述的基于极化特征和近邻传播聚类的极化SAR图像分类方法，其特征是：其中步骤(2)所述的对每个像素点进行Freeman分解，按如下步骤进行：

(2a)读入图像的每个像素点，该像素点为一个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{HV}}^{*}\⟩& \⟨S_{\mathrm{HH}}{S}_{\mathrm{VV}}\⟩\\ \sqrt{2}\⟨S_{\mathrm{HV}}{S}_{\mathrm{HH}}^{*}\⟩& 2\⟨{\left|S_{\mathrm{HV}}\right|}^2\⟩& \sqrt{2}\⟨S_{\mathrm{HV}}{S}_{\mathrm{VV}}^{*}\⟩\\ \⟨S_{\mathrm{VV}}{S}_{\mathrm{HH}}^{*}\⟩& \sqrt{2}\⟨S_{\mathrm{VV}}{S}_{\mathrm{HV}}^{*}\⟩& \⟨{\left|S_{\mathrm{VV}}\right|}^2\⟩\end{array}\right]---1)$

其中，U是中间变量， $U=\frac{1}{\sqrt{2}}\×\left(\begin{array}{ccc}1& 0& 1\\ 1& 0& -1\\ 0& \sqrt{2}& 0\end{array}\right),$ H表示水平极化，V表示垂直极化，S_HH表示水平向发射和水平向接收的回波数据，S_VV表示垂直向发射和垂直向接收的回波数据，S_HV表示水平向发射垂直向接收的回波数据，()^*表示这个数据的共轭，〈·〉表示按视数平均；

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

$C=F_s\left[\begin{array}{ccc}{\left|\mathrm{\β}\right|}^2& 0& \mathrm{\β}\\ 0& 0& 0\\ {\mathrm{\β}}^{*}& 0& 1\end{array}\right]+F_d\left[\begin{array}{ccc}{\left|\mathrm{\α}\right|}^2& 0& \mathrm{\α}\\ 0& 0& 0\\ {\mathrm{\α}}^{*}& 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]---2)$

其中，F_s为平面散射分量的分解系数，F_d为二面角散射分量的分解系数，F_v为体散射分量的分解系数，β是水平发射水平接收后向散射反射系数与垂直发射垂直接收后向散射发射系数的比值，α被定义为α＝I_HI_VH/I_VI_VV，I_H和I_P分别表示地表的水平及垂直反射系数，I_VH和I_VV分别表示竖直墙体的水平及垂直反射系数；

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

$\left\{\begin{array}{c}\⟨{\left|S_{\mathrm{HH}}\right|}^2\⟩=F_s{\left|\mathrm{\β}\right|}^2+F_d{\left|\mathrm{\α}\right|}^2+F_v\\ \⟨{\left|S_{\mathrm{VV}}\right|}^2\⟩=F_s+F_d+F_v\\ \⟨S_{\mathrm{HH}}{S}_{\mathrm{VV}}^{*}\⟩=F_s\mathrm{\β}+F_d\mathrm{\α}+F_v/3\\ \⟨{\left|S_{\mathrm{HV}}\right|}^2\⟩=F_v/3\end{array}\right.---3);$

(2d)计算像素点协方差矩阵C中的

的值，如果

则令α＝-1，如果

则令β＝1，给定α或β的值后，剩余的4个未知数则根据式3）求解，其中Re(·)表示取实部；

(2e)根据求解出的F_s,F_v,F_d,α,β，求解出体散射功率P_v、二面角散射功率P_d和表面散射功率P_s：

$\begin{array}{c}{P}_v=\frac{8F_v}{3}\\ P_d=F_d(1+{\left|\mathrm{\α}\right|}^2)\\ P_s=F_s(1+{\left|\mathrm{\β}\right|}^2)\end{array}---4).$

4.根据权利要求1所述的基于极化特征和近邻传播聚类的极化SAR图像分类方法，其特征是：其中步骤(11)所述的用近邻传播聚类对特征向量矩阵X进行聚类，得到整幅图像的预分类结果，按如下步骤进行：

(11a)根据特征向量矩阵X计算新相似矩阵S，

S(i,j)＝-||x_i-x_j||²

(11b)将新相似矩阵S的对角线上的值都设为I，I取值为-1000，将归属度矩阵AS和吸引度矩阵LS初始值设置为0；

(11b)将归属度矩阵AS和吸引度矩阵LS分别赋值给中间变量a和b，即a＝AS，b＝LS；

(11c)根据新相似度矩阵S和归属度矩阵AS，计算吸引度矩阵LS，吸引度矩阵LS中的每一个元素LS_iu：

${\mathrm{LS}}_{\mathrm{iu}}=S_{\mathrm{iu}}-\underset{\mathrm{\φ}\≠u}{\max }(S_{\mathrm{i\φ}}+{\mathrm{AS}}_{\mathrm{i\φ}}),$

其中，i＝1,...,21，u＝1,...,21，φ＝1,...,21；

(11d)根据新相似度矩阵S和吸引度矩阵LS，计算归属度矩阵AS，归属度矩阵AS中每一个元素AS_iu为；

${\mathrm{AS}}_{\mathrm{iu}}=\left\{\begin{array}{c}\underset{\mathrm{\φ}\≠u}{\mathrm{\Σ}}\max (0,{\mathrm{LS}}_{\mathrm{\φu}}),i=u\\ \min \{0,{\mathrm{LS}}_{\mathrm{uu}}+\underset{\mathrm{\φ}\≠i,u}{\mathrm{\Σ}}\max (0,{\mathrm{LS}}_{\mathrm{\φu}})\},i\≠u\end{array}\right.,$

其中，i＝1,...,21，u＝1,...,21，φ＝1,...,21；

(11e)更新归属度矩阵AS和吸引度矩阵LS：根据已获得的归属度矩阵AS和吸引度矩阵LS，则中间变量g＝(1-λ)×AS+λ×a，中间变量f＝(1-λ)×LS+λ×b，再令AS＝g，LS＝f，更新归属度矩阵AS和吸引度矩阵LS，其中，λ是收敛系数，λ∈[0.5,1)；

(11g)重复步骤(11b)-(11e)，直到迭代次数等于给定的迭代次数n，其中n＝1000，根据LS(u,u)+AS(u,u)的值判断u是否为聚类中心，如果LS(u,u)+AS(u,u)＞0，则u是聚类中心，反之，u不是聚类中心，u＝1,...,21，从而获得预分类结果。

5.根据权利要求1所述的基于极化特征和近邻传播聚类的极化SAR图像分类方法，其特征是：其中步骤(12)所述的对预分类得到的整幅图像用Wishart分类器进行迭代分类，按如下步骤进行：

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

$B_c=\frac{\underset{\mathrm{\ρ}=1}{\overset{n_c}{\mathrm{\Σ}}}{T}_{\mathrm{\ρ}}}{n_c},c=1,...,k,\mathrm{\ρ}=\mathrm{1,2},...,n_c,$

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

(12b)根据下式计算每个像素点的分布特征参数χ_L：

${\mathrm{\χ}}_L=\frac{\mathrm{qL}+1}{q+1}/(\mathrm{RK}-1),$

其中，L为极化SAR图像的视数，q为向量维数，取q=3，RK是每个像素点及其周围的像素点总共九个像素点组成的区域的相对峰值，RK定义如下：

$\mathrm{RK}=\frac{1}{3}(\frac{E\left\{{\left|S_{\mathrm{HH}}\right|}^2\right\}}{E{\left\{\right|S_{\mathrm{HH}}\left|\right\}}^2}+\frac{E\left\{{\left|S_{\mathrm{HV}}\right|}^2\right\}}{E{\left\{\right|S_{\mathrm{HV}}\left|\right\}}^2}+\frac{E\left\{{\left|S_{\mathrm{VV}}\right|}^2\right\}}{E{\left\{\right|S_{\mathrm{VV}}\left|\right\}}^2}),$

其中，S_HH表示水平向发射和水平向接收的回波数据，S_VV表示垂直向发射和垂直向接收的回波数据，S_HV表示水平向发射垂直向接收的回波数据，|·|表示取这个数的模值，E{·}表示取这个数的均值；

(12c)根据每个像素点的分布特征参数χ_L，计算每个像素点到第i类聚类中心V_i的距离d(〈C〉,V_i)：

$d(\⟨C\⟩,V_i)=\mathrm{Lin}\left[V_i\right]+\ln \left(\mathrm{\Γ}\left({\mathrm{\χ}}_L\right)\right)-\frac{{\mathrm{\χ}}_L+\mathrm{qL}}{2}\ln \left({\mathrm{L\χ}}_L\right)-\frac{{\mathrm{\χ}}_L-\mathrm{qL}}{2}\ln \mathrm{Tr}({V_i}^{-1}\⟨C\⟩)-\ln K_{{\mathrm{\χ}}_L-\mathrm{qL}}\left(2\sqrt{L{\mathrm{\χ}}_L\mathrm{Tr}({V_i}^{-1}\⟨C\⟩)}\right)$

其中C是像素点的协方差矩阵，〈·〉表示按视数平均，[·]表示矩阵的行列式，Tr(·)表示矩阵的迹，V_i ^-1表示对聚类中心V_i求逆，Γ(·)表示Gamma函数，表示(χ_L-qL)阶的第二类修正贝塞尔函数，L为图像的视数，q为向量维数，取q=3；

(12e)根据每个像素点到第i类聚类中心的距离对极化SAR图像进行重新划分：如果d(〈C〉,B_τ)≤d(〈C〉,B_ψ)，则将该像素点划分为第τ类，如果d(〈C〉,B_τ)＞d(〈C〉,B_ψ)，则将该像素点划分为第ψ类，其中d(〈C〉,B_τ)是像素点到第τ类聚类中心的距离，d(〈C〉,B_ψ)是像素点到第ψ类聚类中心的距离，τ,ψ＝1,...,8,τ≠ψ；

(12f)重复步骤(12a)-(12e)直到迭代次数等于给定的迭代次数μ，其中μ=3，得到更为准确的分类结果。

Images