CN102722883B - 一种具有空间自适应性的极化sar图像分割方法 - Google Patents

Abstract

本发明提供一种具有空间自适应性的极化SAR图像分割方法，主要解决现有的分割技术在空间复杂度自适应性上的不足导致分割结果无法体现图像细节信息的问题。该方法首先结合H/α-ML Wishart聚类和四叉树分解得到大小不等、能够自适应场景复杂度的初始分割区域，然后利用复Wishart分布和马尔科夫随机场(MRF)对初始分割区域的大小和形状进行调整，得到最终的分割结果。本发明充分有效利用了极化SAR图像的极化信息，具有良好的空间自适应性，分割结果能够很好的保留极化SAR图像中的细节信息，分割速度较快，结果较为精确。

Description

一种具有空间自适应性的极化SAR图像分割方法

技术领域

本发明涉及图像处理领域，更具体的涉及一种具有空间自适应性的极化合成孔径雷达（SAR）图像分割方法。

背景技术

相比与传统的合成孔径雷达（SAR），极化SAR利用多个通道的散射信息，可以获得对目标更加全面的认识。经过几十年的发展，极化SAR的理论和应用研究都取得了巨大进步，各种极化SAR系统平台获取了大量的地面数据，对这些多极化数据进行分析和解译已经成为近年的研究热点。而图像分割正是其中非常关键的一个步骤，由于严重的相干斑噪声的影响，极化SAR图像的分割十分困难。基于边缘检测的方法和基于像素聚类的方法对极化SAR图像进行分割时往往需要确定经验阈值，极易受到相干斑噪声的影响，缺乏灵活性和鲁棒性，通常不能形成连续的边界，无法获得有效的区域信息。虽然可以通过图像预处理（滤波等）和后处理改善分割效果，但是并不能从根本上解决上述问题。而基于区域的分割方法结合图像中像素的各种特征信息，充分利用像素的空间相关性，能够有效的抑制相干斑噪声的影响。这类方法中最具代表性的是基于Markov随机场的分割方法。

基于MRF的方法是一种能够有效利用像素结构信息的分割方法，这类方法基于图像数据的局部相关性，用一个二维随机场描述图像模型。它假定待分割图像的像素只与其邻域内的像素相关，与邻域外的像素无关。基于该假设可定量计算图像局部的先验结构信息，并参与到最大后验概率(MAP)或最大后验边缘概率(MPM)准则的分割处理，从而能有效地利用像素结构信息，得到较好的分割精度。吴永辉等人利用极化数据的目标协方差矩阵和相干矩阵服从复Wishart分布，将MRF与Wishart分布相结合，提出了一种基于Wishart-MRF的极化SAR图像分割方法。该方法先将极化SAR图像分割为大量互不重叠的矩形子区域，然后将MRF与描述协方差矩阵的Wishart分布结合起来，利用MAP准则和ICM算法调整子区域的形状和大小，从而充分利用极化SAR数据的统计先验知识和相邻像素的空间相关性，使分割效果较好的符合地物分类情况。该方法具有较高的分割精度，分割结果中孤立像素少，图像连通性较好，但是缺点也十分明显。该方法没有考虑图像中场景复杂度的变化，在初始分割中对图像采取统一的矩形分块，在MRF中采用全局的能量平滑系数，因此分割不具备空间的自适应性，对图像中的复杂区域分割效果较差。

发明内容

本发明目的是针对现有的极化SAR图像分割方法在空间自适应性方面的不足之处，以提高极化SAR图像分割的准确性，提供了一种结合了H/α-ML Wishart聚类、四叉树分解和Wishart-MRF的极化SAR图像分割方法。该方法的基本思路是先利用H/α-ML Wishart聚类将极化SAR图像按照不同的极化散射类型分为8类，得到初始聚类结果。在此基础上结合四叉树分解得到大量互不重叠、大小不等，具有空间自适应性的正方形子区域，最后利用Wishart-MRF对子区域的形状和大小进行调整，从而实现极化SAR图像的精确分割。

本发明提供的具有空间自适应性的极化SAR图像分割方法包括以下具体步骤：

步骤1，预处理，即对极化SAR图像做多视处理和Lee滤波，降低图像中的相干斑噪声。

步骤2，初始聚类，对步骤1中得到的滤波后的极化SAR图像做如下处理：

步骤2-1，H/α分类，对滤波图像的极化相干矩阵T做极化目标分解，分别计算极化熵H和平均散射角α，利用H和α将极化SAR图像按照不同的极化散射类型分成8类；

步骤2-2，基于像素的最大似然（ML）分类，将步骤2-1中的H/α分类结果作为训练样本送入ML分类器，分别计算第x类中所有元素的极化相干矩阵的平均值Vx，x={1,2…8}。然后根据如下公式，计算极化SAR图像每一个像素i处的相干矩阵T与Vx的相异性d(T,Vx)，对图像进行重新分类：

$d(T,V_x)=\ln \left|V_x\right|+\mathrm{Tr}\left(V_x^{-1}T\right)$

其中Tr(·)代表矩阵的迹，像素i的类别由最小的d(T,Vx)决定。最后将分类结果再次作为训练样本送入ML分类器进行迭代分类，直到分类结果趋于稳定为止，即迭代前后的分类结果中每一个像素对应的类别标号不再发生改变；

步骤2-3，多数投票，即根据步骤2-2中的分类结果，在每一个像素的邻域内进行多数投票，以决定该像素最终的类别，由此可以消除分类结果中的异常值，得到最终的初始聚类图，也就是极化SAR图像分类结果图。

步骤3，四叉树分解，即对步骤2得到的极化SAR图像分类结果图做四叉树分解，具体做法是将分类图像分为相同大小的四块，判断每一块内所有像素的类别标签是否一致，同时判断检查块的大小是否达到设定阈值。若判断结果为否，则将每一块再次分解为相同大小的子块；若判断结果为是，则停止分解。

步骤4，软分割，主要是对步骤3得到的每一块分割区域进行大小和形状调整，具体做法如下：

步骤4-1，计算同质性协同参数C_H，在步骤2-3得到的初始聚类结果图上的每一点处定义一个矩形窗，该点为矩形中心，分别计算窗口内属于每一类的像素个数占窗口内像素总数的百分比，取其中的最大值作为该点处的参数C_H，C_H用来衡量这一点处矩形窗内的区域同质性，以对下一步MRF模型中的空间平滑系数β进行自适应的调整。

步骤4-2，将复Wishart分布与MRF相结合，利用最大后验概率（MAP）准则和条件迭代模式（ICM）算法对步骤3得到的四叉树分解结果进行迭代调整，使得调整后的分割结果更好地符合地物的实际分布。研究表明多视极化相干矩阵T服从复Wishart分布：

$p^{\left(L\right)}\left(T\right)=\frac{L^{\mathrm{qL}}{\left|T\right|}^{L-q}\exp [-L\·\mathrm{Tr}\left(V_k^{-1}T\right)}{R(L,q){\left|V_k\right|}^L}$

R(L,q)=π^q(q-1)/2Γ(L)…Γ(L-q+1)

其中L是视数，q为极化通道数目，V_k=E{T}是多视极化相干矩阵的空间统计平均，R(L,q)为归一化常数，Г(·)代表Gamma函数：

当x为正整数时，Г(x)=(x-1)!。

此外，图像分割可视为已知观测图像时像素类别标号的估计问题，令Y={y_s，s∈S}为观测图像，X={x_s，s∈S}为图像中的像素类别标号，则根据MAP准则和贝叶斯公式：

$p\left(x_s\right|y_s)=\frac{p\left(y_s\right|x_s)p\left(x_s\right)}{p\left(y_s\right)}$

x_s的取值应使后验概率最大，即：

${\hat{x}}_s=\arg \underset{x\∈\{\mathrm{1,2}...K\}}{\max }p\left(x_s\right|y_s)=\arg \underset{x\∈\{\mathrm{1,2}...K\}}{\max }\left\{p\left(y_s\right|x_s)p\left(x_s\right)\right\}$

式中，K为类别数，p(x_s|y_s)为在给定观测数据时类别标号的后验概率，p(y_s|x_s)为在给定类别标号时观测数据的条件概率，p(x_s)为在给定类别标号时观测数据的条件概率，由于p(y_s)与x_s无关，在最大化过程中可以省略。由于在整幅图像内求解该式的全局最优解需要的计算量较大，在此我们采用ICM迭代算法求取局部最优解代替全局最优解，即利用坐标为s的像素的观测值y_s和邻域标号x_ηs确定该像素的标号，因此上式可改写为：

${\hat{x}}_s=\arg \underset{x\∈\{\mathrm{1,2}...K\}}{\max }p\left(x_s\right|y_s)=\arg \underset{x\∈\{\mathrm{1,2}...K\}}{\max }\left\{p\left(y_s\right|x_s)p\left(x_s\right|x_{\mathrm{\ηs}})\right\}$

根据MRF的定义，易知X为马尔科夫随机场，可以用吉布斯分布描述为：

$p\left(x_s\right|x_{\mathrm{\ηs}})=\frac{\exp \{\mathrm{\β}\·u\left(x_s\right)\}}{\underset{x_s=1}{\overset{K}{\mathrm{\Σ}}}\exp \{\mathrm{\β}\·u\left(x_s\right)\}}$

其中u(x_s)等于坐标为s的像素邻域内具有相同类别标号的像素个数，β>0为空间平滑参数，其值越大，分割结果越平滑。

假设步骤3中的四叉树分解最终将整张图像划分为K个子区域，由此结合复Wishart分布与MRF，根据MAP准则和区域同质性协同参数C_H，在每一次的ICM迭代调整中，第i个像素所属分割区域的标号为：

${\hat{x}}_i=\arg \underset{x\∈\{\mathrm{1,2}...K\}}{\max }\left\{\frac{L^{\mathrm{qL}}{\left|T\right|}^{L-q}\exp [C_H\·\mathrm{\βu}\left(x\right)-L\·\mathrm{Tr}\left(V_k^{-1}T\right)]}{R(L,q){\left|V_k\right|}^L{\mathrm{\Σ}}_{x=1}^K\exp \{C_H\·\mathrm{\βu}\left(x\right)\}}\right\}$

R(L,q)=π^q(q-1)/2Γ(L)…Γ(L-q+1)

其中L是视数，q为极化通道数目，β是MRF模型中的空间平滑系数，一般等于1.4。C_H·β是当前像素处经过自适应调整后的平滑系数，u(x)为与像素i具有相同的分割区域标号的像素个数，V_k是第k个子区域内的相干矩阵的平均值，T为像素i处的相干矩阵，R(L,q)为归一化常数，Г(·)代表Gamma函数：

当x为正整数时，Г(x)=(x-1)!。

值得一提的是，在现有的MRF模型中，每一个像素点处的空间平滑系数β是固定不变的，而在本发明中通过计算C_H*β。对每个像素点处的平滑系数进行了自适应调整。

步骤5，后处理，在步骤4的区域调整过程中会产生一些面积过小、毫无意义的分割块，需要将这种分割块与其邻域内散射特性最为接近的分割块进行融合，得到最终的极化SAR图像分割结果。

本发明与现有的极化SAR图像分割方法相比，充分有效利用了极化SAR图像的极化信息，并且考虑了图像中场景复杂度的变化，分割结果具有良好的空间自适应性，能够很好的保留极化SAR图像中的细节信息，分割速度较快，结果较为精确。

附图说明

图1为本发明的实现流程图；

图2-图5分别为本发明实验采用的极化SAR图像，地表真值图，训练样本图及其地物类型说明图；

图6和图7分别为本发明实验的分割结果图和相应的监督分类结果图；

图8和图9分别为Wishart-MRF方法分割结果图和相应的监督分类结果图。

为让本发明的上述内容能更明显易懂，下文特举较佳实施例，并配合所示附图，作详细说明如下。

具体实施方式

下面结合附图与实施例对本发明做进一步说明。

实验数据是NASA-JPL机载L波段AIRSAR系统于1991年获取的荷兰Flevoland地区多视全极化SAR图像，如图2。该图像被广泛运用于极化SAR图像分割、分类实验之中。该图像尺寸大小为320*400，主要是由6类农作物区域组成，其余为未知农作物区域。对该图像按照本发明图1所示流程进行处理，具体步骤如下：

步骤1如图1中S11，对该图像做Lee滤波，降低相干斑噪声。滤波窗口为3*3大小的矩形窗；

步骤2，初始聚类，对步骤1中得到的滤波图像做如下处理：

步骤2-1如图1中S121，H/α分类，对滤波后极化SAR图像的相干矩阵T做特征分解，在此基础上分别计算极化熵H和平均散射角α，根据Cloude-Pottier目标分解分类器利用H和α将极化SAR图像按照不同的极化散射类型分成8类，即有序二面角、有序植被、有序地表面、二面角、植被、地表面、随机植被和随机地表面；

步骤2-2如图1中S122，基于像素的最大似然（ML）分类，将步骤2-1中的H/α分类结果作为训练样本送入ML分类器，分别计算第x类中所有元素的极化相干矩阵的平均值Vx，x={1,2…8}。然后根据如下公式，计算极化SAR图像每一个像素i处的相干矩阵T与Vx的相异性d(T,Vx)，对图像进行重新分类：

$d(T,V_x)=\ln \left|V_x\right|+\mathrm{Tr}\left(V_x^{-1}T\right)$

其中Tr(·)代表矩阵的迹，像素i的类别由最小的d(T,Vx)决定。最后将分类结果再次作为训练样本送入ML分类器进行迭代分类，迭代的终止条件为：

(1)迭代前与迭代后，图像中每一个像素所属的类别标号不再发生改变；

(2)最大迭代次数为30次；

步骤2-3如图1中S123，多数投票，即根据步骤2-2中的分类结果图，在每一个像素的邻域内进行多数投票，以占最大优势的类别标号决定该像素最终的类别，由此可以消除分类结果图中的异常值，得到最终的初始聚类结果。

步骤3如图1中S13，四叉树分解，即对步骤2得到的极化SAR图像分类结果图做四叉树分解，具体做法是设定分解阈值，即最小分解块大小为4*4，然后将分类图像分为相同大小的四块，判断每一块内所有像素处的类别标号是否一致，同时判断检查块的大小是否达到设定阈值。若判断结果为否，则将每一块再次分解为相同大小的子块；若判断结果为是，则停止分解。分解结束后，即可得到具有自适应空间复杂度的初始分割图像；

步骤4，软分割，主要是对步骤3得到的每一块分割区域进行大小和形状调整，具体做法如下：

步骤4-1如图1中S141，计算同质性协同参数C_H，在步骤2-3得到的极化SAR图像分类结果图上的每一点处定义一个矩形窗，该点为矩形中心，分别计算窗口内属于每一类的像素个数N_x与窗口内像素总数N_tot的比值，取其中的最大值作为该点处的参数C_H，C_H用来衡量这一点处矩形窗内的区域同质性，以对下一步MRF模型中的空间平滑系数β进行自适应的调整；

步骤4-2如图1中S142，将复Wishart分布与MRF相结合，利用最大后验概率（MAP）准则和条件迭代模式（ICM）算法对步骤3中四叉树分解得到的分割区域进行迭代调整，使得调整后的分割结果更好地符合地物的实际分布。实验中四叉树分解最终将整张图像划分为K个子区域（K=4985），则在每一次的调整中，第i个像素所属分割区域的标号为：

${\hat{x}}_i=\arg \underset{x\∈\{\mathrm{1,2}...K\}}{\max }\left\{\frac{L^{\mathrm{qL}}{\left|T\right|}^{L-q}\exp [C_H\·\mathrm{\βu}\left(x\right)-L\·\mathrm{Tr}\left(V_k^{-1}T\right)]}{R(L,q){\left|V_k\right|}^L{\mathrm{\Σ}}_{x=1}^K\exp \{C_H\·\mathrm{\βu}\left(x\right)\}}\right\}$

R(L,q)=π^q(q-1)/2Γ(L)…Γ(L-q+1)

其中L是视数，q为极化通道数目，对于全极化数据q=3。β是MRF模型中的空间平滑系数，一般等于1.4。C_H·β是当前像素处经过自适应调整后的平滑系数。u(x)为与像素i具有相同的分割区域标号的像素个数，V_k是第k个子区域内的所有像素处相干矩阵的平均值，T为像素i处的相干矩阵，R(L,q)为归一化常数，Г(·)代表Gamma函数：

当x为正整数时，Г(x)=(x-1)!。

ICM迭代次数为20次。

步骤5如图1中S15，后处理，在步骤4的区域调整过程中会产生一些面积过小（一般小于5个像素）、毫无意义的分割块，需要将这种分割块与其邻域内散射特性最为接近的分割块进行融合，得到最终的极化SAR图像分割结果。但是这种融合同时会影响图像细节的分割效果，因此在实验中并没有进行这一步处理。

最终的分割结果如图6。然后根据图3的地表真值图，从每一类中选取16*16大小的区域作为训练样本，如图4。用ML-Wishart分类器对每一块分割区域进行分类：

$d(V_n,T_m)=\ln \left|T_m\right|+\mathrm{Tr}\left(T_m^{-1}{V}_n\right)$

实验中的训练样本一共有6类，分割区域总数为N，Tm为第m类训练样本中所有元素的极化相干矩阵的平均值，Vn为第n个分割区域中所有元素的相干矩阵的平均值，最终得到极化SAR图像监督分类结果如图7。同时采用吴永辉的分割方法对极化SAR图像进行分割，选择默认参数：初始分割区域大小为5*5，迭代次数20次，分割结果如图8。然后使用同样的训练样本和分类器对分割图像进行监督分类，监督分类结果如图9。

本发明步骤2中所述的初始聚类是先利用H/α-ML Wishart将地物按照不同的极化散射类型分类，这是一个无监督的分类过程，分类结果是在一定程度上真实反映了地物的复杂程度，然后以此为依据进行自适应复杂度的分割。而图7所示的分类结果属于监督分类，是在已知部分地表真值的情况下选取合适的训练样本进行分类，能够体现本发明在保持图像细节方面的优越性。

从分类结果中可以看出，本发明提供的分割方法与吴永辉的方法相比，具有良好的空间自适应性，在图像中的同质性区域分割块较大，块数较少。而在图像中的复杂区域分割块较小，块数较多。在总分割块数相差不大的情况下，有效的保留了图像中的细节信息，如图像中纵向的细长地带和图像上方中部连续的六块散射较强的小区域，确实改善了极化SAR图像分割方法对于场景复杂度的自适应性，具有良好的分割效果。

虽然本发明已以较佳实施例揭露如上，然其并非用以限定本发明，任何熟习此技艺者，在不脱离本发明的精神和范围内，当可作些许的更动与润饰，因此本发明的保护范围当视申请专利范围所界定者为准。

Claims (4)

1.一种具有空间自适应性的极化SAR图像分割方法，其特征在于，包括以下步骤：

步骤1，预处理操作；

步骤2，初始聚类操作：对步骤1中得到的滤波后的极化SAR图像做初始聚类处理，具体包括如下步骤：

步骤2-1，H/α分类，得到H/α分类结果；

步骤2-2，基于像素的最大似然分类，直到分类结果趋于稳定为止；

步骤2-3，多数投票，根据步骤2-2中的分类结果，在每一个像素的邻域内进行多数投票，以决定该像素最终的类别，由此可以消除分类结果中的异常值，得到极化SAR图像分类结果图；

步骤3，四叉树分解操作：对步骤2得到的极化SAR图像分类结果图做四叉树分解，具体步骤为：将分类图像分为相同大小的四块，判断每一块内所有像素的类别标签是否一致，同时判断检查块的大小是否达到设定阈值，若判断结果为否，则将每一块再次分解为相同大小的子块，若判断结果为是，则停止分解；

步骤4，软分割操作，对步骤3得到的每一块分割区域进行大小和形状调整；

步骤5，后处理操作，将不符合预先设定要求的分割块与其邻域内散射特性最为接近的分割块进行融合，得到最终的极化SAR图像分割结果；

其中，步骤4的软分割处理具体包括如下步骤：

步骤4-1，计算同质性协同参数C_H，在步骤2得到的极化SAR图像分类结果图上的每一点处定义一个矩形窗，该点为矩形中心，分别计算窗口内属于每一类的像素个数占窗口内像素总数的百分比，取其中的最大值作为该点处的参数C_H，C_H用来衡量这一点处矩形窗内的区域同质性；

步骤4-2，将复Wishart分布与MRF相结合，利用最大后验概率准则和条件迭代模式算法对步骤3得到的四叉树分解结果进行迭代调整，使得调整后的分割结果更好地符合地物的实际分布，四叉树分解最终将整张图像划分为K个子区域，结合复Wishart分布与MRF，根据最大后验概率准则和区域同质性协同参数C_H，在每一次的条件迭代调整中，第i个像素所属分割区域的标号为：

${\hat{x}}_i=\arg \underset{x\∈\{\mathrm{1,2}...K\}}{\max }\left\{\frac{L^{\mathrm{qL}}{\left|T\right|}^{L-q}\exp [C_H\·\mathrm{\βu}\left(x\right)-L\·\mathrm{Tr}\left(V_k^{-1}T\right)]}{R(L,q){\left|V_r\right|}^L{\mathrm{\Σ}}_{x=1}^K\exp \{C_H\·\mathrm{\βu}\left(x\right)\}}\right\}$

R(L,q)＝π^q(q-1)/2Γ(L)…Γ(L-q+1)

其中L是视数，q为极化通道数目，β是MRF模型中的空间平滑系数，等于1.4，C_H﹒β是当前像素处经过自适应调整后的平滑系数，u(x)为与像素i具有相同的分割区域标号的像素个数，V_k是第k个子区域内的相干矩阵的平均值，T为像素i处的相干矩阵，R(L,q)为归一化常数，Г(·)代表Gamma函数。

2.根据权利要求1所述的一种具有空间自适应性的极化SAR图像分割方法，其特征在于，步骤1所述的预处理操作，是对极化SAR图像做多视处理和Lee滤波，降低图像中的相干斑噪声。

3.根据权利要求1所述的一种具有空间自适应性的极化SAR图像分割方法，其特征在于，步骤2-1所述的H/α分类，是对滤波图像的极化相干矩阵T做极化目标分解，分别计算极化熵H和平均散射角α，利用H和α将极化SAR图像按照不同的极化散射类型分成8类。

4.根据权利要求1所述的一种具有空间自适应性的极化SAR图像分割方法，其特征在于，所述的基于像素的最大似然分类，是将步骤2-1中的H/α分类结果作为训练样本送入ML分类器，分别计算第y类中所有元素的极化相干矩阵的平均值Vy，y＝{1,2…8}，然后根据如下公式，计算极化SAR图像每一个像素i处的相干矩阵T与Vy的相异性d(T,Vy)，对图像进行重新分类：

$d(T,V_y)=\ln \left|V_y\right|+\mathrm{Tr}\left(V_y^{-1}T\right)$

其中Tr(·)代表矩阵的迹，像素i的类别由最小的d(T,Vy)决定，最后将分类结果再次作为训练样本送入ML分类器进行迭代分类，直到分类结果趋于稳定为止。