结合极化特征和分水岭的极化SAR图像分类方法

摘要

本发明公开了一种结合极化特征和分水岭的极化SAR图像分类方法。本发明实现的步骤为：(1)读入极化合成孔径雷达SAR数据；(2)计算像素标记；(3)划分区域；(4)标定区域标记；(5)标定区域大小；(6)更新区域标记；(7)分配颜色，输出最终分类结果。本发明相比现有技术基于区域的阈值划分方法，保证了分类结果中区域一致性，较好的保持了图像中不同区域的边缘信息，解决了极化信息利用不充分和区域划分不完整的问题。本发明可应用于极化合成孔径雷达SAR图像目标检测和目标识别。

说明书

技术领域

本发明属于图像处理技术领域，更进一步涉及极化合成孔径雷达(Synthetic Aperture  Radar,SAR)图像分类技术领域中的一种结合极化特征和分水岭的极化合成孔径雷达 SAR图像分类方法。本发明可应用于对极化合成孔径雷达SAR图像的不同区域准确 地进行分类、目标识别。

背景技术

合成孔径雷达SAR图像分类是通过对极化数据的分析，划分出感兴趣的目标或 区域的过程。合成孔径雷达SAR图像分类是遥感图像分类的主要组成部分，随着雷 达技术的发展，利用极化合成孔径雷达SAR图像分类的研究日益受到人们的重视。 极化合成孔径雷达SAR图像分类技术在地质勘探、森林状态监控、海冰分布监视以 及农作物生长状况评估等多个领域得到广泛的应用。

由于雷达图像在成像过程中不可避免地会产生相干斑，从分类的角度看，相干 斑可被认为是噪声，而基于像素的分类方法对噪声的敏感度很高，所以极化合成孔径 雷达SAR图像的分类方法经常是基于目标或对象的区域处理方式。基于区域处理的 方法主要为基于图像分割的方法，例如：

W.Yonghui等发表的论文“Region based classification of polarimetric SAR images  using Wishart MRF”（IEEE Trans on Geoscience and Remote Sensing Letters,2008, 5(4):668-672）提出了一种基于区域的Wishart MRF（Region-based-WMRF）估计的方 法。该方法先利用MRF过分割算法得到划分的子区域，再利用Wishart分类器对每 一个区域重新标记，从而有效地提高了分类的精度，但是该方法仍然存在的不足是， 该方法属于阈值划分方法，分类结果中明显地出现了块效应，不能很好地保持各类的 区域完整性。

上海交通大学在其专利申请“一种具有空间自适应性的极化SAR图像分割方法” （专利申请号：20121011531.9，公开号：CN102722883A）中提出一种具有空间自适 应性的极化SAR图像分割方法。该方法首先对极化数据进行Cloude分解处理，对 处理结果进行基于像素的最大似然分类，然后对分类结果进行四叉树分解得到初始分 割区域，最后利用Wishart分布和马尔科夫随机场对初始分割区域的大小和形状进行 调整，得到最终的分割结果。但是这种方法存在的不足是，Cloude分解只能将图像 划分为8类且四叉树分解采用阈值来设定区域大小，这样会与真实目标产生偏差。

发明内容

本发明的目的在于克服上述现有技术的不足，提出一种结合极化特征和分水岭 的极化合成孔径雷达SAR图像分类方法。本发明相比基于阈值划分区域的方法提高 了目标保持的准确性，充分利用了极化合成孔径雷达SAR图像的极化特征和分水岭 区域划分算法，在保证图像分类信息的完整性的同时，提高了极化合成孔径雷达SAR 图像分类的质量。

本发明实现上述目的的思路是：先对待分类极化合成孔径雷达SAR图像进行散 射角和功率图的联合似然概率计算像素的最大似然标记，然后对待分类图像进行修正 的分水岭算法划分区域，再对划分的区域标定区域的标记和标定区域的大小，最后采 用两种不同的策略分别对大小区域进行区域标记的更新，得到最终的分类结果。

本发明的步骤包括如下：

(1)读入一幅待分类的极化合成孔径雷达SAR图像中相干矩阵T。

(2)计算像素的标记；

(2a)提取待分类的极化合成孔径雷达SAR图像中相干矩阵T的对角线元素 T₁₁,T₂₂,T₃₃，用0.4，0.2，0.4分别依次对对角线元素T₁₁,T₂₂,T₃₃进行加权处理，将 加权处理后的对角线元素相加，对相加的结果进行归一化处理至0到255之间， 得到功率值；

(2b)对待分类的极化合成孔径雷达SAR图像中相干矩阵T，进行Cloude分 解，得到散射角α；

(2c)在待分类的极化合成孔径雷达SAR图像中每一类地物上，按照SAR图 像总像素点的0.5%随机取像素点，分别提取每一类地物选取像素点所对应的功 率值和散射角，将提取的功率值和散射角作为训练数据；

(2d)计算待分类的极化合成孔径雷达SAR图像中每一个像素点的联合似然 概率，得到各像素点的最大似然标记。

(3)划分区域；

(3a)将像素的最大似然标记与待分类的极化合成孔径雷达SAR图像中功率 值的梯度图进行叠加，得到修正后的梯度图；

(3b)对修正后的梯度图进行分水岭划分，得到划分的区域。

(4)标定区域标记；

(4a)如果划分区域中存在像素标记相同的像素个数达到二分之一以上，将该 像素标记标定为划分区域的标记，若划分区域中不存在像素标记相同的像素个数 达到二分之一以上，则执行步骤(4b)，否则，执行步骤(4c)；

(4b)找出划分区域中相邻区域的标记最多的类别，将该类别标定为划分区 域的标记；

(4c)判断待分类的极化合成孔径雷达SAR图像划分区域的标记是否全部完 成，若否，重复步骤(4a)、步骤(4b)，若是，执行步骤(5)。

(5)标定区域大小；

(5a)依次读入待分类的极化合成孔径雷达SAR图像中划分区域，判断读入 的划分区域的像素个数是否大于100，若大于，执行步骤(5b)，否则，执行步骤 (5c)；

(5b)将读入的划分区域标定为大区域，

(5c)将读入的划分区域标定为小区域；

(5d)判断待分类的极化合成孔径雷达SAR图像中划分区域的大小标记是否 全部完成，若否，重复步骤(5a)、步骤(5b)、步骤(5c)，若是，执行步骤(6)。

(6)更新区域标记；

(6a)用聚类中心公式，计算每一类的聚类中心；

(6b)依次读取极化合成孔径雷达SAR图像标定的区域大小，若为小区域，执 行步骤(6c)，否则，执行步骤(6d)；

(6c)对小区域，用以下最大后验概率公式，得到更新的区域标记：

${\hat{x}}_r=\arg \underset{x_r}{\max }\left\{P\left(y_r\right|x_r)P\left(x_r\right|{\eta }_{x_r})\right\}$

其中，表示小区域r更新后的区域标记，arg表示求取满足条件的参数操 作，max表示求解参数最大值的操作，x_r表示小区域r当前的区域标记，x_r的 取值范围为1,2,...K，K表示极化合成孔径雷达SAR图像所划分的类别数，r表 示极化合成孔径雷达SAR图像中的小区域，y_r表示小区域r的观测数据， P(y_r|x_r)表示小区域的最大后验概率中的似然概率，表示小区域的最 大后验概率中的背景先验概率，表示小区域r的相邻区域；

(6d)对大区域，用以下最大似然概率公式，得到更新的区域标记：

其中，表示大区域更新的区域标记，arg表示求取满足条件的参数操作， max表示求解参数最大值的操作，表示大区域当前的区域标记，的取值 范围为1,2,...K，K表示极化合成孔径雷达SAR图像所划分的类别数，表示 极化合成孔径雷达SAR图像中大区域，表示大区域的最大似然概率中 的似然概率，表示大区域的观测数据；

(6e)判断读取的极化合成孔径雷达SAR图像区域是否完成，若否，重复步 骤(6b)、步骤(6c)、步骤(6d)，若是，执行步骤(6f)；

(6f)判断更新的区域标记迭代次数是否达到20次，若否，重复步骤(6a)、步 骤(6b)、步骤(6c)、步骤(6d)、步骤(6e)，若是，执行步骤(7)。

(7)用红色R、绿色G、蓝色B三个颜色分量作为三基色，根据三基色原理给步 骤6中的区域标记结果分配颜色，得到彩色分类结果图。

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

第一，由于本发明采用散射角和功率图的联合似然概率计算像素的最大似然标 记，克服了现有技术极化目标分解方法所带来的像素标记结果中像素类别不准确的不 足，使得本发明所标定的像素标记能够对极化协方差数据进行有效的信息提取，保证 了像素标记的准确性和极化散射特性的完整性。

第二，由于本发明采用修正的分水岭算法划分区域，克服了现有技术阈值方法所 带来的划分结果中区域边缘不准确的不足，使得本发明所划分的区域更符合实际地物 的边缘分布。

第三，由于本发明采用两种不同的策略分别对大小区域进行区域标记的更新，克 服了现有技术对相干斑噪声敏感的不足，使得本发明提高了抗噪声的性能，有效地保 持了同质区域，也得到了理想的分类边缘。

附图说明

图1为本发明的流程图；

图2为本发明和现有的Region-based-WMRF分类方法，在Flevoland极化合成 孔径雷达SAR测试图像上的仿真示意图；

图3为本发明和现有的Region-based-WMRF分类方法，在SanFranciscobay极 化合成孔径雷达SAR测试上的图像仿真示意图。

具体实施方式

下面结合附图1对本发明的步骤做进一步的详细描述。

步骤1.读入一幅待分类极化合成孔径雷达SAR图像的相干矩阵T。

读入极化合成孔径雷达SAR图像的每个像素点，每个像素点为一个含有9个元 素的极化相干矩阵T：

$T=\left(\begin{array}{ccc}⟨\frac{1}{2}{|S_{\mathrm{EE}}+S_{\mathrm{VV}}|}^2⟩& ⟨\frac{1}{2}(S_{\mathrm{EE}}+S_{\mathrm{VV}}){(S_{\mathrm{EE}}-S_{\mathrm{VV}})}^{*}⟩& ⟨{(S_{\mathrm{EE}}+S_{\mathrm{VV}})S_{\mathrm{EV}}}^{*}⟩\\ ⟨\frac{1}{2}{(S_{\mathrm{EE}}+S_{\mathrm{VV}})}^{*}(S_{\mathrm{EE}}+S_{\mathrm{VV}})⟩& ⟨\frac{1}{2}{|S_{\mathrm{EE}}-S_{\mathrm{VV}}|}^2⟩& ⟨{(S_{\mathrm{EE}}-S_{\mathrm{VV}})S_{\mathrm{EV}}}^{*}⟩\\ ⟨{(S_{\mathrm{EE}}+S_{\mathrm{VV}})}^{*}{S}_{\mathrm{EV}}⟩& ⟨{(S_{\mathrm{EE}}-S_{\mathrm{VV}})}^{*}{S}_{\mathrm{EV}}⟩& ⟨2{\left|S_{\mathrm{EV}}\right|}^2⟩\end{array}\right)$

其中，S_EE表示水平向发射和水平向接收的回波数据，E表示水平极化，S_VV表示 垂直向发射和垂直向接收的回波数据，V表示垂直极化，S_EV表示水平向发射和垂直 向接收的回波数据，*表示共轭转置操作，〈·〉表示按雷达的多视参数求取平均值的 操作，|·|表示求取模值操作。

步骤2.计算像素的标记。

提取待分类极化合成孔径雷达SAR图像中相干矩阵T的对角线元素T₁₁,T₂₂,T₃₃， 用0.4，0.2，0.4分别依次对对角线元素T₁₁,T₂₂,T₃₃进行加权处理，将加权处理后的对 角线元素相加，对相加的结果进行归一化处理至0到255之间，得到功率图，并对待 分类极化合成孔径雷达SAR图像的相干矩阵T进行Cloude分解得到散射角α。

所述Cloude分解步骤如下：

按照下式，对相干矩阵T进行特征分解：

$T=\left[U_3\right]\left(\begin{array}{ccc}{\lambda }_1& 0& 0\\ 0& {\lambda }_2& 0\\ 0& 0& {\lambda }_3\end{array}\right){\left[U_3\right]}^{*}$

其中，T表示相干矩阵，[U₃]表示酉矩阵，$\left[U_3\right]=\left(\begin{array}{ccc}{u}_1& u_2& u_3\end{array}\right),$u₁,u₂,u₃分别表 示相干矩阵T的正交特征向量，λ₁,λ₂,λ₃表示相干矩阵T的特征值，*表示共轭转置 操作。

计算每个像素点的散射随机程度概率：

$G_{\phi }=\frac{{\lambda }_{\phi }}{\underset{j=1}{\overset{3}{\Sigma }}{\lambda }_j}$

其中，G_φ表示每个像素点的散射随机程度概率，φ表示相干矩阵T的特征值的个 数，φ＝1,2,3，λ_φ和λ_j均表示相干矩阵T的特征值。

计算每个像素点的散射熵：

$H=-\underset{\phi =1}{\overset{3}{\Sigma }}{G}_{\phi }{\log }_3\left(G_{\phi }\right)$

其中，H表示每个像素点的散射熵，G_φ表示目标散射类型参量，φ表示特征值 的个数。

计算每个像素点的目标散射类型参量：

α_φ＝acos(abs(u_φ))×180/π

其中，α_φ表示每个像素点的目标散射类型参量，u_φ表示相干矩阵T的正交特征 向量，φ表示特征值的个数，φ＝1,2,3。

计算每个像素点的散射角：

$\alpha =\underset{\phi =1}{\overset{3}{\Sigma }}{\Gamma }_{\phi }{\alpha }_{\phi }$

其中，α表示每个像素点表征目标散射机制的散射角，G_φ表示散射随机程度概 率，α_φ表示目标散射类型参量，φ表示特征值的个数，φ＝1,2,3。

在待分类极化合成孔径雷达SAR图像的每一类地物的上随机取像素点，其 数目是图像总像素点的0.5%，将每一类地物所取像素点所对应的功率图和散射 角分别提取出来，作为训练数据，利用训练数据计算散射角和功率图的联合似然 概率，得到各像素点的最大似然标记，联合似然概率公式如下：

$P(N,G|x_s=m)=\frac{1}{\sqrt{{\sigma }_{1m}·{\sigma }_{2m}}}\exp \{-\frac{{(N-{\mu }_{1m})}^2}{{2\sigma }_{1m}^2}-\frac{{(G-{\mu }_{2m})}^2}{{2\sigma }_{2m}^2}\}$

其中，P(·)表示待分类的极化合成孔径雷达SAR图像中像素的联合似然概率，N 表示像素s的功率值，G表示像素s的散射角，x_s表示像素s的类别标记，s表示极 化合成孔径雷达SAR图像的像素点，m表示待分类的极化合成孔径雷达SAR图像的 类别，m的取值范围为1,2,...,K，K表示待分类的极化合成孔径雷达SAR图像的总 类别数，σ_1m表示第m类的功率值方差，σ_2m表示第m类的散射角方差，exp{·}表示 求取指数操作，μ_1m表示第m类的功率值均值，μ_2m表示第m类的散射角均值。

步骤3.划分区域。

将像素的最大似然标记与待分类的极化合成孔径雷达SAR图像中功率值的梯度 图进行叠加，得到修正后的梯度图，对修正后的梯度图进行分水岭划分，得到划分的 区域。

步骤4.标定区域标记。

第一步，如果划分区域中存在像素标记相同的像素个数达到二分之一以上，将该 像素标记标定为划分区域的标记，若划分区域中不存在像素标记相同的像素个数达到 二分之一以上，则执行步骤4的第二步，否则，执行步骤4的第三步。

第二步，找出划分区域中相邻区域的标记最多的类别，将该类别标定为划分区域 的标记。

第三步，判断待分类的极化合成孔径雷达SAR图像划分区域的标记是否全部完 成，若否，重复步骤4的第一步、步骤4的第二步，若是，执行步骤5。

步骤5.标定区域大小。

第一步，依次读入待分类的极化合成孔径雷达SAR图像中划分区域，判断读入的 划分区域的像素个数是否大于100，若大于，执行步骤5的第二步，否则，执行步骤 5的第三步。

第二步，将读入的划分区域标定为大区域。

第三步，将读入的划分区域标定为小区域。

第四步，判断待分类的极化合成孔径雷达SAR图像中划分区域的大小标记是否全 部完成，若否，重复步骤5的第一步、步骤5的第二步、步骤5的第三步，若是，执 行步骤6。

步骤6.更新区域标记。

第一步，用聚类中心公式，计算每一类的聚类中心：

${\Sigma }_i=\frac{1}{N_i}\underset{f=1}{\overset{N_i}{\Sigma }}{T}_f$

其中，Σ_i表示极化合成孔径雷达SAR图像中第i类的聚类中心，N_i表示极化合 成孔径雷达SAR图像中属于第i类像素点的总个数，T_f表示极化合成孔径雷达SAR图 像中属于第i类像素点的相干矩阵，i表示极化合成孔径雷达SAR图像分类类别的个 数。

第二步，依次读取极化合成孔径雷达SAR图像标定的区域大小，若为小区域，执 行步骤6的第三步，否则，执行步骤6的第四步。

第三步，对小区域，用以下最大后验概率公式，得到更新的区域标记：

${\hat{x}}_r=\arg \underset{x_r}{\max }\left\{P\left(y_r\right|x_r)P\left(x_r\right|{\eta }_{x_r})\right\}=\arg \underset{x_r}{\max }\left\{\frac{\exp \{\mathrm{\beta u}\left(x_r\right)-\mathrm{LTr}\left({\Sigma }_{x_r}^{-1}{T}_r\right)\}}{{\left|{\Sigma }_{x_r}\right|}^L}\right\}$

其中，表示小区域r更新的区域标记，arg表示求取满足条件的参数操作，max 表示求解参数最大值的操作，x_r表示小区域r当前的区域标记，x_r的取值范围为 1,2,...K，K表示极化合成孔径雷达SAR图像所划分的类别数，P(y_r|x_r)表示小区 域的最大后验概率中的似然概率，y_r表示小区域r的观测数据，表示小区 域的最大后验概率中的背景先验概率，表示小区域r的相邻区域，β表示空间平 滑参数，u(x_r)表示阶跃函数，L为雷达的多视参数，Tr(·) 表示矩阵的迹，为第x_r类的聚类中心，T_r为区域r的相干矩阵。

第四步，对大区域，用以下最大似然概率公式，得到更新的区域标记：

其中，表示大区域更新的区域标记，arg表示求取满足条件的参数操作，max 表示求解参数最大值的操作，表示大区域当前的区域标记，的取值范围为 1,2,...K，K表示极化合成孔径雷达SAR图像所划分的类别数，表示极化合成孔 径雷达SAR图像中大区域，表示大区域的最大似然概率中的似然概率，表示大区域的观测数据。

第五步，判断读取的极化合成孔径雷达SAR图像区域是否完成，若否，重复步 骤6的第二步、步骤6的第三步、步骤6的第四步，若是，若是，执行步骤6的第六 步。

第六步，判断更新的区域标记迭代次数是否达到20次，若否，重复步骤6的第 一步、步骤6的第二步、步骤6的第三步、步骤6的第四步、步骤6的第五步，若是， 执行步骤7。

步骤7.用红色R、绿色G、蓝色B三个颜色分量作为三基色，根据三基色原理 给步骤6中的区域标记结果分配颜色，得到彩色分类结果图。

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

1仿真条件：

本发明的仿真是在主频2.5GHZ的Pentium Dual_Core CPU E5200、内存1.98GB 的硬件环境和MATLAB R2010a的软件环境下进行的。

2仿真内容：

图2为本发明和现有的Region-based-WMRF分类方法，在Flevoland极化合成 孔径雷达SAR测试图像上的仿真示意图，其中，图2(a)为仿真中使用的原始Flevoland 极化合成孔径雷达SAR测试图像，于1989年8月16日由NASA/JPL AIRSAR雷达 系统获取，大小为750×102，4雷达的多视参数L为4。图2(b)为用现有的 Region-based-WMRF分类方法对图2(a)测试图像的分类仿真结果图，图2(c)为用本 发明对图2(a)测试图像的分类仿真结果图。

图3为本发明和现有的Region-based-WMRF分类方法，在SanFranciscobay极 化合成孔径雷达SAR测试上的图像仿真示意图，其中，图3(a)为仿真中使用的原始 SanFrancisco bay极化合成孔径雷达SAR测试图像，于1994年10月3日由NASA/JPL ARISAR雷达系统获取，大小为900×1024，雷达的多视参数L为4。图3(b)为用现有 的Region-based-WMRF分类方法对图3(a)测试图像的分类仿真结果图，图3(c)为用 本发明对图3(a)测试图像的分类仿真结果图。

3仿真效果分析：

图2(b)为用现有的Region-based-WMRF分类方法对图2(a)测试图像的分类仿 真结果图。由图2(b)的仿真结果可以看到，对于同质区域保持有一定的改善，但块效 应明显，孤立小区域多，边缘边界的清晰性和细节完整性都比较差。图2(c)为用本 发明对图2(a)测试图像的分类仿真结果图。由图2(c)的仿真结果可以看到，本发明的 分类结果从视觉上具有更好地分类性能，能有效区分各类地物，区域一致性明显改善， 不同区域之间分类后的边缘更加平滑清晰，细节信息完整。

图3(b)为用现有的Region-based-WMRF分类方法对图3(a)测试图像的分类仿 真结果图，由图3(b)的仿真结果可以看到，有明显的小区域存在，并且地物边缘不平 滑；图3(c)为用本发明对图3(a)测试图像的分类仿真结果图，由图3(c)的仿真结果可 以看到，本发明的结果明显好于现有的Region-based-WMRF分类的结果，分类图 中地物细节保留完整，图像视觉效果较好，区域的保持更加完整且更加精确。