CN103365985B - 一种类别自适应的极化sar分类方法 - Google Patents

Abstract

一种类别自适应的极化SAR分类方法，属于图像处理技术领域。其分类过程为：对输入数据进行Freeman分解得散射功率矩阵P_s，P_d，P_v；根据P_s，P_d，P_v值将极化SAR图像数据初始划分为三类；计算每个像素的同极化比R，选取两个不同的阈值将每类进一步分为三类；计算每类中各像素点的自极化参数δ，按照δ值将每类细分成N类；对细分结果表征类别差异性，获得相异矩阵R_D，用可视化聚类趋势估计法对R_D重排，得到新矩阵R_D ¹；将R_D ¹变换为相异图像Im，对Im黑框识别处理，得到类别数和聚类中心；对所有输入的极化SAR数据用复Wishart迭代分类并上色，得到最终彩色分类结果图。本发明用于极化SAR图像目标识别的分割区域一致性好保留信息完整。

Description

一种类别自适应的极化SAR分类方法

技术领域

本发明属于图像处理领域，具体地说是针对极化SAR图像分类，可应用于对极化SAR的目标检测和目标识别。

背景技术

与传统的合成孔径雷达SAR系统相比，极化SAR可以得到更丰富的目标信息，大大提高对地物的识别能力。因此，极化SAR已成为SAR的发展趋势，并在军事领域、地质与资源勘探、地形测绘和制图、海洋应用与研究、水资源应用、农业和林业应用等多方面有广泛的应用。其中，地物分类是极化SAR图像解译的重要内容。目前经典分类方法有：

1997年，Cloude和Pottier提出的基于H/α极化分解的分类方法，并成为基于散射机制进行极化SAR图像分类的典型代表，见CloudeSR，PottierE.AnentropybasedclassificationschemeforlandapplicationsofpolarimetricSAR[J].IEEETrans.Geosci.RemoteSensing.1997，35(1)：549-557.

1999年，Lee等人提出了基于H/α目标分解和Wishart分类器的H/α-Wishart分类方法，见LeeJS，GrunesMR，AinsworthTL，etal.UnsupervisedclassificationusingpolarimetricdecompositionandthecomplexWishartclassifier[J].IEEEClassificationusingpolarimetricdecompositionandthecomplexWishartclassifier[J].IEEETrans.GeosciRemoteSensing.1999，37(5)：2249-2258.该方法是在原始H/α分类基础上增加Wishart迭代，主要是对H/α划分后的8类利用Wishart分类器对每一个像素进行重新划分，从而有效的提高了分类的精度，但是不能很好的保持各类的极化散射特性。

2011年，西安电子科技大学ShuangWang等利用Freeman分解得到的三种散射功率对图像进行初始大类的划分，并利用同极化比对初始大类进行更细致的划分，最后在初始类别划分的基础上又对整幅图像进行复Wishart迭代，进一步提高了分类精度，见文献ShuangWang，JingjingPei，KunLiu，etal.UnsupervisedclassificationofPOLSARdatabasedonthepolarimetricdecompositionandtheco-polarization.GeoscienceandRemoteSensingSymposium(IGARSS)，2011IEEEInternational，2011：424-427.该算法思想简单，易于理解，具有较高的分类精度，但是该算法仍然存在着一定的局限性，该算法分类类别数目固定不变，一般为9类，因此对于类别多于9类或少于9类的数据来说，该算法的分类效果会有所影响。

发明内容

本发明的目的在于克服上述已有技术的不足，在上述基于Freeman分解和同极化比的极化SAR图像分类方法的基础上，提出了结合另一个自极化参数、可视化聚类趋势估计算法(VAT，VisualAssessmentof(Cluster)Tendency)和黑框识别算法(DBE，darkblockextraction)，实现分类数目的自适应，同时提高了分类效果。

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

(1)对输入数据进行Freeman分解，得到三种散射功率矩阵P_s，P_d，P_v，其中，P_s表示表面散射功率矩阵，P_d表示二面角散射功率矩阵，P_v表示体散射功率矩阵；

(2)计算同极化比R，并根据功率矩阵P_s，P_d，P_v和同极化比R对极化SAR图像数据初始化并划分成9类，具体步骤如下：

2a)利用下式计算每个像素点的同极化比R值，得到每类数据的一系列R值：

$R=10*\log \left(\frac{{\left|S_{\mathrm{HH}}\right|}^2}{{\left|S_{\mathrm{VV}}\right|}^2}\right)---\left(1\right)$

其中H表示水平极化，V表示垂直极化，S_HH表示水平向发射和水平向接收的回波数据，S_VV表示垂直向发射和垂直向接收的回波数据，|·|表示取这个数的模值；

2b)根据max(P_s，P_d，P_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对应像素点划分为体散射类；

2c)选取两个不同的阈值threshod1和threshod2，将步骤2b)划分的每类进一步划分为三类，即将R＜threshod1对应的像素点划分为一类，将threshod1＜R＜threshod2对应的像素点划分为一类，将R＞threshod2对应的像素点划分为一类，从而将极化SAR图像数据划分为9类，其中，threshod1＝-2，threshod2＝2；

(3)计算自极化参数δ，利用自极化参数δ对步骤(2)的每一类再进行细分为N类，得到共9N类，具体步骤如下：

3a)利用下式计算步骤2c)每一类中每个像素点的自极化参数δ值：

$\mathrm{\δ}=2*\frac{{\left|S_{\mathrm{HV}}\right|}^2}{{\left|S_{\mathrm{HH}}\right|}^2}---\left(2\right)$

为了计算方便，类比式(1)，将式(2)变为

$\mathrm{\δ}=10*\log \left(\frac{{\left|S_{\mathrm{HV}}\right|}^2}{{\left|S_{\mathrm{HH}}\right|}^2}\right)---\left(3\right)$

其中，S_HV表示水平向发射垂直向接收的回波数据；

3b)对每类按照δ的值将每类数据按间隔均匀细分成N类，得到共9N类，其中N＝30；

(4)对步骤(3)细分结果计算类别差异性，获得相异矩阵R_D，用可视化聚类趋势估计算法对R_D重排序，得到新矩阵R_D ¹；

(5)将矩阵R_D ¹变换为相异图像Im1，对Im1图像进行黑框识别，获得类别数目n和聚类中心V_i，其中i＝1，2，...，n，具体步骤如下：

5a)将R_D ¹转化成0～255的灰度图像Im1，用Otsu自动阈值法将图像Im1二值化处理，然后用形态学腐蚀膨胀去除噪声，得到一组沿着对角线黑框组成的图像Im2；

5b)对图像Im2作距离变换转化得到灰度图像Im3，将Im3图像上所有像素的灰度值投影到主对角线上，构成一维的投影信号H1；

5c)对H1进行平滑，得到信号H2，求H2的一阶导数，根据导数的过零点检测投影信号H2中存在的波峰和波谷，波峰的个数即为数据集的类别数目n，第i个波峰位置所对应的数据作为一个聚类中心V_i，其中i＝1，2，...，n；

(6)利用步骤(5)获得的类别数n和聚类中心V_i，其中i＝1，2，...，n，用复Wishart迭代方法对所有输入的极化SAR数据分类；

(7)用红色R、绿色G、蓝色B三个颜色分量作为三基色，给分类结果上色，得到最终彩色分类结果图。

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

1、本发明由于使用可视化聚类趋势估计算法和黑框识别算法，实现了类别自适应的极化SAR分类；

2、本发明结合了Freeman分解、自极化比R、自极化参数δ和复Wishart迭代分类，充分使用了极化信息，得到更好的分类结果；

3、仿真结果表明，本发明方法较传统的H/α-Wishart方法分类有更好的分类结果效果，改变了ShuangWang的基于Freeman分解和同极化比的极化SAR图像分类方法的类别不变问题，同时获得了较好的结果。

附图说明

图1是本发明的流程示意图；

图2是的本发明方法对Flevoland数据用可视化聚类趋势估计算法和黑框识别算法处理过程的部分过程图；

图3是用本发明与ShuangWang的基于Freeman分解和同极化比的极化SAR图像分类方法、H/α-Wishart分类方法对Flevoland数据的分类仿真结果图；

图4是用本发明与ShuangWang的基于Freeman分解和同极化比的极化SAR图像分类方法、H/α-Wishart分类方法对SanFranciscoBay数据分类的仿真结果图。

具体实施方式

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

步骤一，对输入数据进行Freeman分解，得到三种散射功率矩阵，P_s(表面散射功率矩阵)P_d(二面角散射功率矩阵)P_v(体散射功率矩阵)；

Freeman分解见文献AnthonyFreeman，Athree-componentScatteringModelforPolarimetricSARData.IEEETrans.Geosci.RemoteSensing，36(3)：963-973，May，1998.其主要思想是假设每一个像素都由三种散射类别构成，具体步骤如下：

1a)读入数据的每个像素点为一个含有9个元素的3×3极化协方差矩阵C；

$C=\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]---\left(1\right)$

其中，H表示水平极化，V表示垂直极化，S_HH表示水平向发射和水平向接收的回波数据，S_VV表示垂直向发射和垂直向接收的回波数据，S_HV表示水平向发射垂直向接收的回波数据，^*表示这个数据的共轭，<·>表示按视数平均，|·|表示取这个数的模值；

1b)根据Freeman分解的思想将协方差矩阵C重新表示为表示为：

$C=f_v\left[\begin{array}{ccc}1& 0& 1/3\\ 0& 2/3& 0\\ 1/3& 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}_s\left[\begin{array}{ccc}{\left|\mathrm{\&beta;}\right|}^2& 0& \mathrm{\&beta;}\\ 0& 0& 0\\ {\mathrm{\&beta;}}^{*}& 0& 1\end{array}\right]---\left(2\right)$

其中，f_s为平面散射分量的分解系数，f_d为二面角散射分量的分解系数，f_v为体散射分量的分解系数，β是水平发射水平接收后向散射反射系数与垂直发射垂直接收后向散射发射系数的比值，α被定义为α＝R_ghR_vh/R_gvR_vv，R_gh和R_gv分别表示地表的水平及垂直反射系数，R_vh和R_vv表示竖直墙体的水平及垂直反射系数；

1c)令式令式(2)右边的矩阵的元素与式(1)右边的矩阵的元素相对应，可以获得一个具有五个未知数f_s，f_v，f_d，α，β和四个方程的方程组，其具体表示如下：

$\left\{\begin{array}{c}{<\left|S_{\mathrm{HH}}\right|}^2>=f_s{\left|\mathrm{\&beta;}\right|}^2+f_d{\left|\mathrm{\&alpha;}\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{\&beta;}+f_d\mathrm{\&alpha;}+f_v/3\\ <{\left|S_{\mathrm{HV}}\right|}^2>=f_v/3\end{array}\right.---\left(3\right)$

1d)计算像素点协方差矩阵C中的的值并判断正负，如果令α＝-1，如果则令β＝1，其中Re(·)表示取实部，给定α或β的值后，根据式(3)求解剩余4个未知数.根据求解出的f_s，f_v，f_d，α，β，按照下式求解出散射功率矩阵P_s，P_d，P_v；

$P_v=\frac{8f_v}{3},P_d{=f}_d(1+{\left|\mathrm{\&alpha;}\right|}^2),P_d=f_d(1+{\left|\mathrm{\&alpha;}\right|}^2)---\left(4\right)$

其中，P_v表示体散射功率矩阵，P_d表示二面角散射功率矩阵，P_s表示表面散射功率矩阵。

步骤二，计算同极化比R，并根据功率矩阵P_s，P_d，P_v和同极化比R对极化SAR图像数据进行初始划分为9类：

2a)利用下式计算每个像素点的同极化比R值，得到每类数据的一系列R值：

$R=10*\log \left(\frac{{\left|S_{\mathrm{HH}}\right|}^2}{{\left|S_{\mathrm{VV}}\right|}^2}\right)---\left(5\right)$

其中，H表示水平极化，V表示垂直极化，S_HH表示水平向发射和水平向接收的回波数据，S_VV表示垂直向发射和垂直向接收的回波数据，|·|表示取这个数的模值；

2b)根据max(P_s，P_d，P_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对应像素点划分为体散射类；

2c)选取两个不同的阈值threshod1和threshod2，将上一步划分的每类进一步划分为三类，即将R＜threshod1对应的像素点划分为一类，将threshod1＜R＜threshod2对应的像素点划分为一类，将R＞threshod2对应的像素点划分为一类，从而将整个极化SAR图像数据划分为9类，threshod1和threshod2的值可根据不同数据的R值的分布而人为选取，本发明的取值为其中，threshod1＝-2，threshod2＝2。

步骤三，计算自极化参数δ，根据δ值对上述所分的9类进一步细分，每类细分成N类，得到共9N类，具体步骤如下：

3a)计算自极化参数δ，利用下式计算上面每一类中每个像素点的自极化参数δ值：

$\mathrm{\&delta;}=2*\frac{{\left|S_{\mathrm{HV}}\right|}^2}{{\left|S_{\mathrm{HH}}\right|}^2}---\left(6\right)$

为了计算方便，类比式(5)，将式(6)变为

$\mathrm{\&delta;}=10*\log \left(\frac{{\left|S_{\mathrm{HV}}\right|}^2}{{\left|S_{\mathrm{HH}}\right|}^2}\right)---\left(7\right)$

其中，S_HH表示水平向发射和水平向接收的回波数据，S_HV表示水平向发射垂直向接收的回波数据；

3b)对每类按照δ的值将每类细分成N类，得到共9N类，其中，N可以根据数据的多少适当增加或减少，如本文取N＝30。

步骤四，对细分结果计算类别差异性，获得相异矩阵R_D，用可视化聚类趋势估计算法(VAT，VisualAssessmentof(Cluster)Tendency)对R_D重排序，具体步骤如下：

4a)计算每类的最大似然中心，并求每类之间的距离，第i类的最大似然中心V_i为：

$V_i=\frac{1}{N_i}\underset{j=1}{\overset{N_j}{\mathrm{\&Sigma;}}}{C}_j,i=\mathrm{1,2},...,9N$

其中，N_i为第i类的像素个数，C_j为第i类的第j个像素的协方差矩阵，本发明采用的是修订的Wishart距离，见文献BinLiu，HaoHu，HuanyuWang，etal.Superpixel-BasedClassificationWithanAdaptiveNumberofClassesforPolarimetricSARImages，IEEEGeoscienceandRemoteSensing，2013，2(51)：907-924，第i类和第j类的修订的Wishart距离D_ij为：

$D_{\mathrm{ij}}=\frac{1}{2}(\mathrm{Tr}(V_i^{-1}<V_j>)+\mathrm{Tr}(V_j^{-1}<V_i>))-q$

其中，当发射和接收天线是同一个时令q＝3，当发射和接收天线分开时令q＝4，V_i和V_j分别是第i类和第j类的最大似然中心，D_ij组成9N*9N相异矩阵R_D＝{D_ij}，相异矩阵为对称矩阵，对角线元素值为0，得到相异矩阵R_D后，用可视化聚类趋势估计算法处理，参见Bezdek，JamesC，Hathaway，etal.VAT：atoolforvisualassessmentof(cluster)tendency，NeuralNetworks，2002.IJCNN′022002，3：2225-2230，步骤如下：

4b)初始化K＝{1，2，...，m}，I＝J＝Φ，P＝(0，...，0)，其中m为类别数，本文中m＝9N，P为长度为m的全零向量；

4c)选择 $(i,j)=\underset{p,q\&Element;K}{\arg \max R_D\mathrm{pq},}P\left(1\right)=i,I=\left\{1\right\},J=K-\left\{i\right\};$

4d)当r＝2:m，选择 $(i,j)=\underset{p\&Element;I,q\&Element;J}{\arg \max R_D\mathrm{pq}}$

P(r)＝j，I＝I∪{j}，J＝J-{j}；

4e)得新矩阵R_D ¹＝R_D ¹ _ij＝R_DP(i)P(j)，1≤i≤m，1≤j≤m，R_D为9N*9N矩阵。

步骤五，将矩阵R_D ¹变换为相异图像Im，对Im图像黑框识别算法(DBE，darkblockextraction)识别，见文献WangLiang，Leckie，ChristopherA，etal.AutomaticallyDeterminingtheNumberofClustersinUnlabeledDataSets，IEEETrans，KnowledgeandDataEngineering，2009，3(21)：335-350，获得类别数目n和聚类中心V_i，其中i＝1，2，...，n，具体步骤如下：

5a)将R_D ¹转化成0～255的灰度图像Im1，用Otsu自动阈值法将Im1图像二值化处理、形态学腐蚀膨胀去除噪声得到图像Im2，是一组沿着对角线黑框组成，其中Im1和Im2都是9N*9N图像；

Otsu自动阈值法是一种简单的算法，见文献N.Otsu.AThresholdSelectionMethodfromGray-levelHistograms，IEEETrans.Systems，Man，andCybernetics，vol.9，no.1，pp.62-66，1979，利用Otsu自动阈值法得到阈值σ，当像素点灰度值大于σ，将该像素点标记为1，反之，当像素点灰度值小于或等于σ，将该像素点标记为0，即将Im1图像二值化了；

图像形态学操作，改变图像的空间结构和目标形状，膨胀和腐蚀是两种最基本的形态学操作，见文献P.Soille，MorphologicalImageAnalysis：PrinciplesandApplications.Springer，1999；

5b)对图像Im2做距离变换(DT，Distancetransform)转化得到灰度图像Im3，将Im3图像上所有像素的灰度值投影到主对角线上，构成一维的投影信号H1；

距离变换是计算并标识空间点(对目标点)距离的过程，二值图像的距离变换是计算二值图像中每个点到最近的非零点的距离，得到距离图像，它最终把二值图像变换为灰度图像，图像Im2做距离变换转化得到9N*9N灰度图像Im3；

将Im3图像上所有像素的灰度值投影到主对角线上，是将所有像素点向主对角线上垂直投影，将在同一条垂线上的像素点的灰度值相加，得到一个数据，这样每条垂线处都会得到一个数据，将每条垂线上的数据按顺序排列，组成一维信号H1，一般H1是有多个波谷和波峰的波形信号，H1的长度为Im3的宽度9N的2倍，即H1的长度为18N；

5c)对H1进行平滑，得到信号H2，求H2的一阶导数，根据导数的过零点检测投影信号H2中存在的波峰和波谷，H2波峰的个数即为数据集的类别数目n，第i个波峰位置所对应的数据作为一个聚类中心V_i，其中i＝1，2，...，n；

对H1进行平滑，采用的方法是均值滤波，即选取均值滤波器与H1信号线性卷积；

H2的一阶导数有多个零值点，求导数时，波峰周围的导数变化趋势是从正值到负值，波谷周围的的导数是从负值到正值，因此，寻找波峰是寻找导数的从正值到负值零点，即为导数下降处的零点，搜寻H2，波峰的个数即为数据集的类别数目n，假设第i个波峰的位置是L，则第i个波峰位置所对应的数据V_i就是步骤4a)中的第L/2类最大似然中心V_L/2，其中i＝1，2，...，n。

步骤六，对利用步骤五获得的类别数n和聚类中心V_i，其中i＝1，2，...，n，用复Wishart迭代方法对所有输入的极化SAR数据分类；

复Wishart迭代是1994年由国外学者提出的一种迭代方法，见文献J.S.Lee，M.R.GrunesandR.Kwok.Classificationofmulti-lookpolarimetricSARimagerybasedoncomplexWishartdistribution[A].Int.J.RemoteSensing[C]，15(11)，1994，15(11)：2299-2311，利用该迭代方法对对整个极化SAR图像数据的n类划分结果进行迭代的步骤如下：

6a)对整个极化SAR图像数根据据的下式计算每个像素点到第i类聚类中心的距离为：

$d_m(<C>,V_i)=\ln \left|V_i\right|+\mathrm{Tr}(V_i^{-1}<C>),i=\mathrm{1,2},...,n$

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

6b)根据每个像素点到第i类聚类中心的距离对极化SAR图像数据进行重新划分：

如果d(<C>，V_i)≤d(<C>，V_j)，i，j＝1，2，...，n，i≠j，则将该像素点划分为第i类，其中C为像素点的协方差矩阵，V_i和V_j分别为第i类和第j类的聚类中心，d(<C>，V_i)和表示该像素点到第i类聚类中心的距离，d(<C>，V_j)表示该像素点到第j类聚类中心的距离；

6c)重新计算每一类的聚类中心，第i类的聚类中心V_i为：

$V_i=\frac{1}{N_i}\underset{j=1}{\overset{N_j}{\mathrm{\&Sigma;}}}{C}_j,i=\mathrm{1,2},...,n$

其中，N_i为第i类的像素个数，C_j为第i类的第j个像素的协方差矩阵；

该步骤中的V_i为初始聚类中心，Wishart迭代中每一次迭代聚类中心都会发生变化，直至达到规定的迭代次数或其他的迭代终止条件，V_i才作为最终的聚类中心；

6d)重复步骤6a)-6c)直到迭代次数等于给定的迭代次数cishu，其中cishu＝5，这里是本发明规定的迭代终止条件，为规定迭代次数，其它的终止条件也可以是前后两次聚类中心浮动，达到人为规定范围，或者是分类结果的变化小于一个百分比，像素点归属稳定。

步骤七，用红色R、绿色G、蓝色B三个颜色分量作为三基色，给分类结果上色，得到最终彩色分类结果图。

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

1仿真内容：应用本发明方法和传统的H/α-Wishart方法，ShuangWang的基于Freeman分解和同极化比的极化SAR图像分类分别对二幅SAR图像进行分类实验，并从分类结果的区域一致性、错分情况、边缘保持、分类类别数等方面评价。

2仿真实验结果

A、三类分类方法对Flevoland数据的分类仿真：

用本发明方法和传统的H/α-Wishart方法，ShuangWang的基于Freeman分解和同极化比的极化SAR图像分类方法对Flevoland数据分类仿真。

本发明方法对Flevoland数据用可视化聚类趋势估计算法和黑框识别算法处理过程的部分过程结果，见图2，其中图2(a)为相异矩阵R_D可视化聚类趋势估计算法重排序后的结果图，图2(b)为Im3图像上所有像素的灰度值投影到主对角线上，构成一维的投影信号H1和平滑信号H2，上方为H1，下方为H2，图2(c)为H2的梯度图。从图2(a)中可以看出沿对角线的黑块，从图2(c)的H2梯度图可以识别出8类聚类中心，即画“o”处的个数为8，每个“o”处对应的位置对应的数据为聚类中心。

三种方法对Flevoland数据分类结果见图3，其中，图3(a)为原始的Flevoland数据，图3(b)为H/α-Wishart分类的结果，图3(c)为的基于Freeman分解和同极化比的极化SAR图像分类，图3(d)为本发明分类结果。从图3可见，前两种方法类别基本固定，H/α-Wishart方法分为8类，有较多区域划分不清楚，区域的一致性不好，差别较小的类别未能分开，而ShuangWang的方法分为9类，划分的较细致，但同质区域的一致性也不是太好。本发明实现类别自适应，对本图分类别数为8类，解决了ShuangWang的方法不能改变分类数目固定为9的问题，且从效果上看，本发明方法在同质区域的一致性比前两种方法更好，对差别较小的类别也可以分开，边缘保持也较好。

B、三类分类方法对SanFranciscoBay数据的分类仿真：

用本发明方法和传统的H/α-Wishart方法，ShuangWang的基于Freeman分解和同极化比的极化SAR图像分类方法对SanFranciscoBay数据分类仿真，分类结果见图4，其中，图4(a)表示原始的SanFranciscoBay数据，图4(b)为H/α-Wishart分类的结果，图4(c)为ShuangWang的基于Freeman分解和同极化比的极化SAR图像分类，图4(d)为本发明分类结果。从图4(b)可见，H/α-Wishart分类方法的分类结果，分为8类，区域划分的较细致，但还有较多区域划分不清楚；从图4(c)可见，ShuangWang的基于Freeman分解和同极化比的极化SAR图像分类方法分为9类，分类结果从视觉上看效果更好，其中区域分类一致性明显好于前种方法，分类后的边缘保持也较好；从图4(d)可见，本方法为10类，分类结果比ShuangWang的基于Freeman分解和同极化比的极化SAR图像分类结果好，本方法区域分的更细致，结构更清晰，同质区域和边缘分的较好。

从图3、图4可见，本发明可以实现类别自适应的极化SAR分类，可以根据极化SAR数据自动选择合适的分类数目，本发明改变了现有的两种方法分类数目固定不变的问题，且分类结果也明显较好，区域划分更加细致且更加精确。

综上所述，本发明提出的对极化SAR数据的分类方法，通过先对数据进行Freeman分解，提取了表征极化特征的三种散射功率，然后有效地结合同极化比和自极化参数，利用可视化聚类趋势估计算法和黑框识别算法，获得自适应的类别数目和聚类中心，最后再对分类进行复Wishart迭代分类，进一步改善了分类的结果。

Claims (4)

1.一种类别自适应的极化SAR分类方法，包括如下步骤：

(1)对输入数据进行Freeman分解，得到三种散射功率矩阵P_s，P_d，P_v，其中，P_s表示表面散射功率矩阵，P_d表示二面角散射功率矩阵，P_v表示体散射功率矩阵；

(2)计算同极化比R，并根据功率矩阵P_s，P_d，P_v和同极化比R对极化SAR图像数据初始化并划分成9类，具体步骤如下：

2a)利用下式计算每个像素点的同极化比R值，得到每类数据的一系列R值：

其中H表示水平极化，V表示垂直极化，S_HH表示水平向发射和水平向接收的回波数据，S_VV表示垂直向发射和垂直向接收的回波数据，|·|表示取这个数的模值；

2b)根据max(P_s，P_d，P_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对应像素点划分为体散射类；

2c)选取两个不同的阈值threshod1和threshod2，将步骤2b)划分的每类进一步划分为三类，即将R＜threshod1对应的像素点划分为一类，将threshod1＜R＜threshod2对应的像素点划分为一类，将R＞threshod2对应的像素点划分为一类，从而将极化SAR图像数据划分为9类，其中，threshod1＝-2，threshod2＝2；

(3)计算自极化参数δ，利用自极化参数δ对步骤(2)的每一类再进行细分为N类，具体步骤如下：

3a)利用下式计算步骤2c)每一类中每个像素点的自极化参数δ值：

为了计算方便，类比式(1)，将式(2)变为

其中，S_HV表示水平向发射垂直向接收的回波数据；

3b)对每类按照δ的值将每类数据按间隔均匀细分成N类，得到共9N类，其中N＝30；

(4)对步骤(3)细分结果计算类别差异性，获得相异矩阵R_D，用可视化聚类趋势估计算法对R_D重排序，得到新矩阵R_D ¹；

(5)将矩阵R_D ¹变换为相异图像Im1，对Im1图像进行黑框识别，获得类别数目n和聚类中心V_i，其中i＝1，2，...，n，具体步骤如下：

5a)将R_D ¹转化成0～255的灰度图像Im1，用Otsu自动阈值法将图像Im1二值化处理，然后用形态学腐蚀膨胀去除噪声，得到一组沿着对角线黑框组成的图像Im2；

5b)对图像Im2作距离变换转化得到灰度图像Im3，将Im3图像上所有像素的灰度值投影到主对角线上，构成一维的投影信号H1；

5c)对H1进行平滑，得到信号H2，求H2的一阶导数，根据导数的过零点检测投影信号H2中存在的波峰和波谷，波峰的个数即为数据集的类别数目n，第i个波峰位置所对应的数据作为一个聚类中心V_i，其中i＝1，2，...，n；

(6)利用步骤(5)获得的类别数n和聚类中心V_i，其中i＝1，2，...，n，用复Wishart迭代方法对所有输入的极化SAR数据分类；

(7)用红色R、绿色G、蓝色B三个颜色分量作为三基色，给分类结果上色，得到最终彩色分类结果图。

2.根据权利要求1所述的一种类别自适应的极化SAR分类方法，步骤(1)所述的对输入数据进行Freeman分解，按如下步骤进行：

2.1)读入数据的每个像素点为一个含有9个元素的3×3极化协方差矩阵C；

其中，*表示这个数据的共轭，<·>表示按视数平均；

2.2)根据Freeman分解的方法将协方差矩阵C重新表示为：

其中，f_s为平面散射分量的分解系数，f_d为二面角散射分量的分解系数，f_v为体散射分量的分解系数，β是水平发射水平接收后向散射反射系数与垂直发射垂直接收后向散射发射系数的比值，α被定义为α＝R_ghR_vh/R_gvR_vv，R_gh和R_gv分别表示地表的水平及垂直反射系数，R_vh和R_vv表示竖直墙体的水平及垂直反射系数；

2.3)令式(5)等号右边的矩阵的元素与式(4)等号右边的矩阵的元素相对应，由此获得一个具有五个未知数f_s，f_v，f_d，α，β和四个方程的方程组，其具体表示如下：

2.4)计算像素点协方差矩阵C中的的值并判断正负，如果令α＝-1，如果则令β＝1，其中Re(·)表示取实部，给定α或β的值后，根据式(6)求解剩余4个未知数，根据求解出的f_s，f_v，f_d，α，β，按照下式求解出散射功率矩阵P_s，P_d，P_v；

其中，P_v表示体散射功率矩阵，P_d表示二面角散射功率矩阵，P_s表示表面散射功率矩阵。

3.根据权利要求1所述的一种类别自适应的极化SAR分类方法，步骤(4)所述的对步骤(3)细分结果表征类别差异性，获得相异矩阵P_D，用可视化聚类趋势估计算法对R_D重排序，得到新矩阵R_D ¹的方法按如下步骤进行：

3.1)计算每类的最大似然中心，并求每类之间的距离，第i类的最大似然中心V_i为：

其中，N_i为第i类的像素个数，C_j为第i类的第j个像素的协方差矩阵，第i类和第j类的修订的Wishart距离D_ij为：

其中，当发射和接收天线是同一个时令q＝3，当发射和接收天线分开时令q＝4，V_i和V_j分别是第i类和第j类的最大似然中心，D_ij组成相异矩阵R_D＝{D_ij}，相异矩阵为对称矩阵，对角线元素值为0，得到相异矩阵R_D后，用可视化聚类趋势估计算法处理，步骤如下：

3.2)初始化K＝{1，2，...，m}，I＝J＝Φ，P＝(0，...，0)，其中m为类别数，本文中m＝9N，P为长度为m的全零向量；

3.3)选择

3.4)当r＝2：m，选择

P(r)＝j，I＝I∪{j}，J＝J-{j}；

3.5)得新矩阵R_D ¹＝R_D ¹ _ij＝R_DP(i)P(j)，1≤i≤m，1≤j≤m。

4.根据权利要求1所述的一种类别自适应的极化SAR分类方法，所述步骤(6)获得的类别数n和聚类中心V_i，其中i＝1，2，...，n，用复Wishart迭代方法对所有输入的极化SAR数据分类，按如下步骤进行：

4.1)利用下式计算每个像素点到第i类聚类中心的距离：

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

4.2)根据每个像素点到第i类聚类中心的距离对极化SAR图像数据进行重新划分：

如果d(<C>，V_i)≤d(<C>，V_j)，i，j＝1，2，...，n，i≠j，则将该像素点划分为第i类，其中C为像素点的协方差矩阵，V_i和V_j分别为第i类和第j类的聚类中心，d(<C>，V_i)和d(<C>，V_j)表示该像素点到第i类和第j类聚类中心的距离；

4.3)重新计算每一类的聚类中心，第i类的聚类中心V_i为：

其中，N_i为第i类的像素个数，C_j为第i类的第j个像素的协方差矩阵；

4.4)重复步骤4.1)-4.3)直到迭代次数等于给定的迭代次数cishu，其中cishu＝5。