基于Freeman分解和同极化比的极化SAR图像分类方法

摘要

本发明公开了一种基于Freeman分解和同极化比的极化SAR图像分类方法，主要解决现有技术计算复杂度较高和分类效果差的问题，其实现步骤为：(1)输入极化SAR数据的协方差矩阵；(2)对输入的矩阵进行Freeman分解，获取平面散射、二面角散射及体散射三种散射功率矩阵；(3)根据三种散射功率矩阵将极化SAR数据进行初始划分；(4)计算每类极化SAR数据各像素点的同极化比；(5)选择阈值依据同极化比将步骤(3)中每类极化SAR数据划分为3类，从而将整个极化SAR数据划分为9类；(6)对整个极化SAR数据的划分结果进行复Wishart迭代并上色，得到最终彩色分类结果图。本发明与经典分类方法相比，对极化SAR数据的划分更加严谨，分类结果明显，计算复杂度相对较小。

说明书

技术领域

本发明属于图像数据处理技术领域，具体地说是一种图像分类方法，该方法可用于对极化SAR数据的分类。

背景技术

随着雷达技术的发展，极化SAR已成为SAR的发展趋势，极化SAR能够得到更丰富的目标信息，在农业、林业、军事、地质学、水文学和海洋等方面具有广泛的研究和应用价值，如地物种类的识别、农作物成长监视、产量评估、地物分类、海冰监测、地面沉降监测，目标检测和海洋污染检测等。极化图像分类的目的是利用机载或星载极化传感器获得的极化测量数据，确定每个像素所属的类别。经典的极化SAR分类方法包括：

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

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

J.S.Lee等基于Freeman分解于提出了一种基于Freeman-Durden分解的多极化图像非监督分类算法，见Lee J S，Grunes M R，Pottier E，et al.Unsupervised terrain classification preserving polarimetric scattering characteristic[J].IEEE Trans.Geosci.Remote Sensing.2004，42(4)：722-731.该方法主要是通过Freeman分解获取表征散射体散射特性的的三个特征：平面散射功率、二面角散射功率和体散射功率，然后根据这三个特征的大小对极化数据进行划分，并对初始划分进行类别合并，最后再利用Wishart分类器对每一个像素进行重新划分。这种算法结合了Freeman散射模型和复Wishart分布，具有保持多极化SAR的主要散射机制纯净性的特性，但是该方法中由于Freeman分解中的多类的划分以及合并，因而计算复杂度较高。

发明内容

本发明的目的在于克服已有技术的不足，在上述基于Freeman-Durden分解的多极化图像非监督分类算法的基础上，提出一种基于Freeman分解和同极化比的极化SAR图像分类方法，以降低计算复杂度并进一步提高分类效果。

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

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

(2)根据功率矩阵P_s，P_d，P_v对极化SAR图像数据进行初始划分：

2a)根据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对应像素点划分为体散射类；

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

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

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

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

(3)对整个极化SAR图像数据的9类划分结果进行复Wishart迭代，得到更为准确的分类结果；

(4)用红色R、绿色G、蓝色B三个颜色分量作为三基色，根据三基色原理给步骤(3)得到更为准确的分类结果上色，得到最终彩色分类结果图。

本发明具有如下优点：

A)本发明以Freeman分解为基础，提取极化SAR图像数据中三种散射功率矩阵P_s，P_d，P_v，并结合同极化比来对极化SAR图像数据进行划分，由于极化SAR图像数据中不同散射体的三种散射功率的分布以及同极化比的大小存在较大的差异，因此结合散射功率和同极化比对极化数据可以进行有效的划分。

B)本发明中将极化SAR图像数据直接划分为9类，避免了Freeman分解中的多类的划分以及合并，在实现上更加简单易行，更加严谨，且降低了计算复杂度。

C)由于极化SAR图像数据没有严格的地物类别数，所以在分类时类别数没有严格的界定，本发明根据具体数据将其划分为9类，由分类的结果可以看出本发明在一些小区域的划分结果上，明显优于现有的极化SAR图像数据分类引用较多的两种经典方法H/α和H/α-Wishart的分类结果，且区域一致性划分更好，不同区域划分后的边缘也更加清晰。

附图说明

图1是本发明的流程图；

图2是本发明使用的两组原始极化SAR数据的图像；

图3是用本发明与现有的H/α及H/α-Wishart分类方法对San Francisco Bay数据的分类仿真结果图；

图4是用本发明与现有的H/α及H/α-Wishart分类方法对Flevoland数据划分的仿真结果图。

具体实施方式

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

步骤1，对输入数据进行Freeman分解，得到散射功率矩阵P_s，P_d，P_v，其中，P_s表示表面散射功率矩阵，P_d表示二面角散射功率矩阵，P_v表示体散射功率矩阵。

Freeman分解见文献Freeman A and Durden S.A three-component scattering model for polarimetric SAR data.IEEE Transactions on Geoscience and Remote Sensing 1998，36(3)：963-973，其具体步骤如下：

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)---1)$

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

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

C＝<C>_s+<C>_d+<C>_v            2)

其中，${<C>}_s=f_s\left(\begin{array}{ccc}{\left|\beta \right|}^2& 0& \beta \\ 0& 0& 0\\ {\beta }^{*}& 0& 1\end{array}\right),$${<C>}_d=f_d\left(\begin{array}{ccc}{\left|\alpha \right|}^2& 0& \alpha \\ 0& 0& 0\\ {\alpha }^{*}& 0& 1\end{array}\right),$${<C>}_V=f_v\left(\begin{array}{ccc}1& 0& 1/3\\ 0& 2/3& 0\\ 1/3& 0& 1\end{array}\right),$

<C>_s为表面散射分量的协方差矩阵，<C>_d为二面角散射分量的协方差矩阵，<C>_v为体散射分量的协方差矩阵，f_s为平面散射分量的分解系数，f_d为二面角散射分量的分解系数，f_v为体散射分量的分解系数，β是水平发射水平接收后向散射反射系数与垂直发射垂直接收后向散射发射系数的比值，α被定义为α＝R_ghR_vh/R_gvR_vv，R_gh和R_gv分别表示地表的水平及垂直反射系数，R_vh和R_vv表示竖直墙体的水平及垂直反射系数，因此C可重新表示成：

$C={<C>}_s+{<C>}_d+{<C>}_v$

$=f_s\left(\begin{array}{ccc}{\left|\beta \right|}^2& 0& \beta \\ 0& 0& 0\\ {\beta }^{*}& 0& 1\end{array}\right)+f_d\left(\begin{array}{ccc}{\left|\alpha \right|}^2& 0& \alpha \\ 0& 0& 0\\ {\alpha }^{*}& 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)---3)$

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

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

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

1e)根据求解出的f_s，f_v，f_d，α，β，按照公式5)求解出散射功率矩阵P_s，P_d，P_v。

$P_v=\frac{8f_v}{3},P_d=f_d(1+{\left|\alpha \right|}^2),P_s=f_s(1+{\left|\beta \right|}^2)---5)$

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

步骤2，根据功率矩阵P_s，P_d，P_v对极化SAR图像数据进行初始划分。

2a)根据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对应像素点划分为体散射类；

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

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

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

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

步骤3，对整个极化SAR图像数据的9类划分结果进行复Wishart迭代，得到更为准确分类结果。

复Wishart迭代是1994年由国外学者提出的一种迭代方法，见文献Lee J S，Grunes M R.Classification of multi-look polarimetric SAR imagery based on complex Wishart distribution[J].Int.J.Remote Sensing.1994，15(11)：2299-2311，利用该迭代方法对对整个极化SAR图像数据的9类划分结果进行迭代的步骤如下：

3a)对整个极化SAR图像数据的9类划分结果，根据下式求每一类的聚类中心V_i：

$V_i=\frac{\underset{j=1}{\overset{N_i}{\Sigma }}{C}_j}{N_i},i=\mathrm{1,2},···9,j=\mathrm{1,2},···N_i$

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

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

3b)根据下式计算每个像素点到第i类聚类中心的距离：

$d(<C>,V_i)=\ln \left[V_i\right]+\mathrm{Tr}(V_i^{-1}<C>),i=\mathrm{1,2},···9$

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

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

如果d(<C>，V_i)≤d(<C>，V_j)i，j＝1，2，…9，j≠i，则将该像素点划分为第i类，如果d(<C>，V_i)＞d(<C>，V_j)i，j＝1，2，…9，j≠i，则将该像素点划分为第j类，其中d(<C>，V_j)表示该像素点到第j类聚类中心的距离；

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

步骤4，用红色R、绿色G、蓝色B三个颜色分量作为三基色，根据三基色原理给步骤3中极化SAR图像数据分类结果上色，将分类后的图像数据的标记类别，对不同的标记类别赋以不同的颜色，得到最终彩色分类结果图根据。

本发明的效果通过以下实验仿真进一步说明。

1.实验条件与方法

实验仿真环境为：MATLAB 7.0.4，Intel(R)Pentium(R)1CPU 2.4GHz，Window XP Professional。

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

2.实验内容及结果

实验内容：本发明分别使用图2所示两组极化SAR数据做测试实验，其中图2(a)为第一组是San Francisco Bay数据，视数为四，图2(b)为第二组极化SAR数据，是荷兰Flevoland省的区域，视数为四，来源于AIRSAR。

实验一，用现有的H/α及H/α-Wishart分类方法与本发明对San Francisco Bay数据进行分类仿真，分类结果见图3，其中，图3(a)为H/α方法分类结果，图3(b)是H/α-Wishart分类的结果，图3(c)为本发明分类结果。

从图3(a)可见H/α方法虽然经典，但分类结果很不理想，很多区域都没有区分出来；

从图3(b)可见，结合H/α方法和Wishart分类器的H/α-Wishart分类方法分类结果明显优于原始的H/α方法，区域划分的更加细致，但还有较多区域划分不清楚；

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

实验二，用现有的H/α及H/α-Wishart分类方法与本发明对Flevoland数据进行分类仿真，分类结果见图4，其中，图4(a)是H/α方法分类结果，图4(b)是H/α-Wishart分类的结果，图4(c)是本发明分类的结果。

从图4可以看出，本发明的结果明显好于现有的两种经典方法分类的结果，区域划分更加细致且更加精确。

综上所述，本发明提出的对极化SAR数据的分类方法，通过先对数据进行Freeman分解，提取了表征极化特征的三种散射功率，然后有效地结合同极化比，实现了对极化SAR数据的有效划分，最后再对分类后结果进行复Wishart迭代，进一步改善了每一类的划分结果。而且本发明思想比较简单，计算复杂度相对较小，容易理解与应用。