CN102663684A - 基于高斯混合模型参数块迁移聚类的sar图像分割方法 - Google Patents

Abstract

本发明公开了一种基于高斯混合模型参数块迁移聚类的SAR图像分割方法，主要解决现有技术分割结果不稳定、精度和区域一致性不理想的问题。其实现过程是：1)输入图像，确定分割类数；2)提取特征；3)用Normalized Cut算法对待分割的SAR图像进行预分割；4)根据块大小排序选择源域和目标域；5)用EM算法对源域进行参数估计；6)对样本、源域和目标域分别求平均得到新样本、新源域和新目标域；7)根据K近邻法把源域参数直接迁移至新目标域；8)把迁移过来的参数当做基于GMM的EM算法的初始参数，对新样本进行EM聚类得到图像最终分割结果。本发明具有分割速度快、效果稳定良好、区域一致性好的优点，可用于雷达目标检测和识别。

Description

基于高斯混合模型参数块迁移聚类的SAR图像分割方法

技术领域

本发明属于图像处理技术领域，涉及SAR图像分割，可用于雷达目标检测和目标识别。

背景技术

合成孔径雷达SAR具有全天时、全天候的探测与侦察能力。它利用脉冲压缩技术获得高的距离分辨率，利用合成孔径原理提高方位分辨率，从而相比真实孔径雷达在遥感领域具有独特的优势。对SAR图像的理解和解译属于图像处理范畴，还涉及信号处理，模式识别及机器学习等众多学科。由于SAR具有的独特作用，使得SAR图像的理解与解译在国防和民用领域正受到越来越广泛的关注，SAR图像分割作为SAR图像后续解译处理的关键环节之一，也就显得愈加重要。现有的SAR图像分割方法大致可以分为基于区域分割方法和基于边缘的分割方法，基于区域的分割方法有聚类、阈值分割、区域生长法等，基于边缘的分割方法有边界跟踪法和图像滤波法等。

其中基于聚类的SAR图像分割方法，是将SAR图像中具有某方面相似特征的区域划分成一类。经典聚类方法EM算法因简单、容易实现而被广泛使用，但它本身还存在对初始值敏感，容易陷入局部最优的缺陷，造成图像分割的不稳定以及错分现象。此外，随着科学技术的发展，在实际SAR图像分割中，经常遇到大规模、超高维、复杂分布的数据，对于这些数据，现有的样本迁移聚类算法速度慢，分割效果不理想，无法在有限的时间内给出令人满意的分割结果。

发明内容

本发明的目的在于克服上述已有技术的不足，提出了一种基于高斯混合模型参数块迁移聚类的SAR图像分割方法，提高SAR图像分割的速度和性能。

实现本发明的技术思路是：从待分割SAR图像中提取小波特征为待划分的样本，同时将SAR图像划分成若干块，选择少数几块作为源域，剩余的块作为目标域，然后从源域学习参数，通过近邻法将源域参数迁移至目标域作为EM算法的初始参数对样本进行EM聚类，得到图像分割结果。具体实现步骤包括如下：

(1)输入待分割SAR图像，确定图像分割类别数C；

(2)对待分割的SAR图像进行3层平稳小波变换，图像像素点总个数为n，对每个像素点提取出10维子带能量特征，构成大小为n×10的输入数据样本集X＝{x_i}，i＝1，2，...，n：

$x_i=\frac{1}{M\×N}\underset{e=1}{\overset{M}{\mathrm{\Σ}}}\underset{u=1}{\overset{N}{\mathrm{\Σ}}}\left|\mathrm{codf}(e,u)\right|$

其中，M×N为利用滑动窗口确定的子带大小，coef(e，u)为平稳小波子带中第e行第u列的系数值；

(3)根据Shi和Malik提出的Normalized Cut算法把待分割的SAR图像划分成M块，将相应的数据样本集X划分成M组样本：X＝{D_r}，r＝1，2，...，M，其中D_r＝{x_b}，b＝1，2，...，s，s表示第r组样本个数大小；

(4)对M组样本数据集X＝{D_r}，r＝1，2，...，M根据组内样本个数大小排序，选择前面的10组为源域X^S＝{D_v}，v＝1，2，...，10，剩余的(M-10)组为目标域X^T＝{D_q}，q＝11，12，...，M；

(5)用步骤(2)提取的样本特征分别对源域X^S＝{D_v}，v＝1，2，...，10的每一组D_v计算均值和协方差，得到10组均值μ_v和协方差σ_v两个类参数，v表示源域组数。

(6)分别对M组样本X＝{D_r}，r＝1，2，...，M中的每一个块D_r＝{x_b}，b＝1，2，...，s求平均得到新样本 $X^{\′}=\left\{{\stackrel{\‾}{D}}_r\right\},r=\mathrm{1,2},...,M,$ 其中 ${\stackrel{\‾}{D}}_r=\underset{b=1}{\overset{s}{\mathrm{\Σ}}}{x}_b/s,$ 则与步骤(4)源域X^S＝{D_v}，v＝1，2，...，10和目标域X^T＝{D_q}，q＝11，12，...，M对应的新源域样本为 $X^{S^{\′}}=\left\{{\stackrel{\‾}{D}}_v\right\},v=\mathrm{1,2},...,10,$ 新目标域样本为 $X^{T^{\′}}=\left\{{\stackrel{\‾}{D}}_q\right\},q=\mathrm{11,12},...,M;$

(7)根据基于欧氏距离的K近邻法对新目标域样本

中的每一个样本

寻找离其最近的新源域样本

并将步骤(5)得到的源域样本类参数均值μ_v和协方差σ_v迁移至新目标域样本

的类参数均值μ_q和协方差σ_q：

μ_q＝μ_v，σ_q＝σ_v；

(8)把迁移过来的那些均值μ_q和协方差σ_q两个类参数当作基于高斯混合模型的EM算法的初始参数，利用该算法对新样本

进行聚类，得到最终图像分割结果。

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

(1)本发明通过基于图论的Normalized Cut算法对原图进行初始划分，实际上是加入了图像像素之间的局部空间信息，从而能有效提高聚类算法的性能。

(2)本发明通过块内所有样本的平均值代替整个块进行参数迁移及聚类，实际上是减少了样本个数，从而比基于样本的参数迁移聚类算法速度快。

仿真实验表明本发明在细节划分和区域一致性上都优于传统EM算法。

附图说明：

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

图2是现有包含山地和平坦地区的SAR图像；

图3是现有包含森林和田地的SAR图像；

图4是本发明与现有方法对图2的仿真分割对比图；

图5是本发明与现有方法对图3的仿真分割对比图。

具体实施方式：

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

步骤1.输入待分割SAR图像，根据图像内容，判断需要识别的主要目标及背景，确定分割类别数C，本实例中的C取值为2。

步骤2.提取待分割SAR图像的特征。

SAR图像不仅数据量大，在成像过程中不同的地物有不同的后向发射和散射特性，从而具有丰富的幅度、相位、极化和纹理信息，并且图像固有的相干斑噪声对于分割性能直接产生影响，因此，有必要在图像分割前对SAR图像进行纹理分析，提取有效的纹理特征进行聚类。

在以上分析基础上，对待分割的SAR图像进行3层平稳小波变换，图像像素点总个数为n，对每个像素点提取出10维子带能量特征，构成大小为n×10的输入数据样本集X＝{x_i}，i＝1，2，...，n：

$x_i=\frac{1}{M\×N}\underset{e=1}{\overset{M}{\mathrm{\Σ}}}\underset{u=1}{\overset{N}{\mathrm{\Σ}}}\left|\mathrm{codf}(e,u)\right|$

其中，M×N为利用滑动窗口确定的子带大小，这里取值为16×16，coef(e，u)为平稳小波子带中第e行第u列的系数值，这样，对于本实例中大小为256×256的SAR图像，就构成65536×10的矩阵作为高斯混合模型块参数迁移聚类方法的输入数据。

步骤3.根据Shi和Malik提出的Normalized Cut算法把待分割的SAR图像划分成M块，将相应的数据样本集X划分成M组样本：X＝{D_r}，r＝1，2，...，M，其中D_r＝{x_b}，b＝1，2，...，s，s表示第r组样本个数大小。

步骤4.对M组样本数据集X＝{D_r}，r＝1，2，...，M根据组内样本个数大小排序，选择前面的10组为源域X^S＝{D_v}，v＝1，2，...，10，剩余的(M-10)组为目标域X^T＝{D_q}，q＝11，12，...，M。

步骤5.用步骤(2)提取的样本特征分别对源域X^S＝{D_v}，v＝1，2，...，10的每一组D_v计算均值和协方差，得到10组均值μ_v和协方差σ_v两个类参数，v表示源域组数。

步骤6.分别对M组样本X＝{D_r}，r＝1，2，...，M中的每一个块D_r＝{x_b}，b＝1，2，...，s求平均得到新样本 $X^{\′}=\left\{{\stackrel{\‾}{D}}_r\right\},r=\mathrm{1,2},...,M,$ 其中 ${\stackrel{\‾}{D}}_r=\underset{b=1}{\overset{s}{\mathrm{\Σ}}}{x}_b/s,$ 则与步骤(4)源域X^S＝{D_v}，v＝1，2，...，10和目标域X^T＝{D_q}，q＝11，12，...，M对应的新源域样本为 $X^{S^{\′}}=\left\{{\stackrel{\‾}{D}}_v\right\},v=\mathrm{1,2},...,10,$ 新目标域样本为 $X^{T^{\′}}=\left\{{\stackrel{\‾}{D}}_q\right\},q=\mathrm{11,12},...,M.$

步骤7.根据基于欧氏距离的K近邻法对新目标域样本

中的每一个样本

寻找离其最近的新源域样本

并将步骤(5)得到的源域样本类参数均值μ_v和协方差σ_v迁移至新目标域样本的类参数均值μ_q和协方差σ_q：

μ_q＝μ_v，σ_q＝σ_v。

步骤8.对新样本进行基于高斯混合模型的EM算法聚类。

现有的基于高斯混合模型的EM算法在类参数初始化时一般都是随机生成，会导致聚类结果不稳定及一定的错分现象。本发明把迁移过来的那些类参数均值μ_q和协方差σ_q当作基于高斯混合模型的EM算法的初始参数，再对新样本

进行聚类，得到最终图像分割结果，具体步骤如下：

(8a)首先针对新样本

定义归属矩阵z_rj：

$z_{\mathrm{rj}}=\left\{\begin{array}{cc}1& {\stackrel{\‾}{D}}_r\∈j\\ 0& {\stackrel{\‾}{D}}_r\∉j\end{array}\right.,r=\mathrm{1,2},...,M;j=\mathrm{1,2},...,C$

其中

是m维向量，M表示样本个数，C表示类别数，若z_rj＝1表示向量

属于第j类，z_rj＝0表示向量

不属于第j类；

(8b)确定初始类参数，包括初始类均值μ_j、类方差σ_j和类混合权重α_j，其中初始类均值μ_j、类方差σ_j为步骤(7)迁移过来的参数均值μ_q和协方差σ_q，即μ_j＝μ_q，σ_j＝σ_q，初始类混合权重α_j为：α_j＝1/C；

(8c)根据上述类参数求得任意一项向量

的概率为：

$p\left({\stackrel{\‾}{D}}_r\right)=\underset{j=1}{\overset{C}{\mathrm{\Σ}}}{\mathrm{\α}}_{j}G({\stackrel{\‾}{D}}_r;{\mathrm{\μ}}_j,{\mathrm{\σ}}_j),i=\mathrm{1,2},...,M$

其中 $G({\stackrel{\‾}{D}}_r,{\mathrm{\μ}}_j,{\mathrm{\σ}}_j)=\frac{1}{{\left(2\mathrm{\π}\right)}^{\frac{m}{2}}{\left|{\mathrm{\σ}}_j\right|}^{\frac{1}{2}}}\exp \{-\frac{1}{2}{({\stackrel{\‾}{D}}_r-{\mathrm{\μ}}_j)}^T{{\mathrm{\σ}}_j}^{-1}({\stackrel{\‾}{D}}_r-{\mathrm{\μ}}_j)\},$ 表示向量

在第j类的高斯分布概率，m为维数；

(8d)利用所求得的概率计算归属矩阵z_rj的期望值E；

$E=E\left(z_{\mathrm{rj}}\right|{\stackrel{\‾}{D}}_r;{\mathrm{\α}}_j,{\mathrm{\μ}}_j,{\mathrm{\σ}}_j)={\mathrm{\α}}_{j}G({\stackrel{\‾}{D}}_r;{\mathrm{\μ}}_j,{\mathrm{\σ}}_j)/p\left({\stackrel{\‾}{D}}_r\right)$

其中表示期望值E是在类混合权重α_j、类均值μ_j和类方差σ_j条件下求得的向量

的归属矩阵z_rj的期望值；

(8e)用产生的期望值E计算新的类混合权重α′_j、类均值μ′_j和类方差σ′_j：

${{\mathrm{\α}}_j}^{\′}=\frac{1}{n}\underset{r=1}{\overset{M}{\mathrm{\Σ}}}E$

${{\mathrm{\μ}}_j}^{\′}=\underset{r=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{Ex}}_{\mathrm{rl}}/\underset{r=1}{\overset{M}{\mathrm{\Σ}}}E,l=\mathrm{1,2},...,m$

${{\mathrm{\σ}}_j}^{\′}=\underset{r=1}{\overset{M}{\mathrm{\Σ}}}E({\stackrel{\‾}{D}}_r-{{\mathrm{\μ}}_j}^{\′}){({\stackrel{\‾}{D}}_r-{{\mathrm{\μ}}_j}^{\′})}^T/\underset{r=1}{\overset{M}{\mathrm{\Σ}}}E;$

(8f)将归属矩阵z_rj的期望值E最大化：

首先，用所述更新后的类均值μ′_j、类方差σ′_j和类混合权重α′_j重新进行步骤(8c)(8e)得到归属矩阵z_rj的新期望E′；

然后，反复迭代步骤(8c)(8e)直至|E′_t+1-E′_t|≤0.000001循环结束，E′_t表示第t次循环得到的z_rj的期望，E′_t+1表示第t+1次循环得到的z_rj的期望；

(8g)利用z_rj的新期望E′更新归属矩阵z_rj，得到样本的最终划分：

z_rj＝1当且仅当 $j=\arg \underset{j=1}{\overset{C}{\max }}\left(E^{\′}\right)$

z_rj＝0当且仅当 $j\≠\arg \underset{j=1}{\overset{C}{\max }}\left(E^{\′}\right).$

本发明的效果可以通过如下仿真实验进行验证：

(一)实验条件设置

选择2幅SAR图像进行验证，图像大小为256×256，划分类别数C取值为2。

(二)仿真内容及结果

A)用现有的EM聚类算法、基于样本的参数迁移聚类算法和本发明方法对图2包含山地和平坦地区的SAR图像进行仿真实验，分割结果如图4所示。其中图4(a)为EM算法分割结果，图4(b)为基于样本的参数迁移聚类算法分割结果，图4(c)为本发明分割结果。

从图4(a)可以看出，EM分割结果中山地部分存在很多误分。

从图4(b)可以看出，基于样本的参数迁移聚类算法比EM算法分割效果好，但区域一致性上表现不理想。

从图4(c)可以看出，本发明方法由于考虑到图像的局部信息，在区域一致性上明显占优势，分割结果明显地呈现山地和平坦地区两部分，这也是符合地质勘探上对地形整体把握需要的。

B)用现有的EM聚类算法、基于样本的参数迁移聚类算法以及本发明方法对图3包含森林和田地的SAR图像进行仿真实验，分割结果如图5所示。其中图5(a)为EM算法分割结果，图5(b)为基于样本的参数迁移聚类算法分割结果，图5(c)为本发明方法分割结果。

从图5(a)可以看出，EM聚类的结果不仅田地中混杂了森林区域，而且轮廓划分也较为粗糙。

从图5(b)可以看出，基于样本的参数迁移聚类算法在田地中的区域一致性表现不理想。

从图5(c)可以看出，本发明方法的分割结果在森林和田地的一致性有很大的改善，同时右下角的细节部分保留较好。

C)将现有的EM聚类算法、基于样本的参数迁移方法，简称STEM和本发明，简称GTEM三种分割方法在图2和图3上的分割时间进行对比，结果如表1所示：

表1.SAR图像三种算法运行时间对比(单位：秒)

图像	EM	STEM	GTEM
				图2	49.4756	43.7900	14.3288
图3	13.0276	61.7988	14.4243

从表1中可以看出，本发明比基于样本的参数迁移方法速度快，这正是由于块内所有样本的平均值替代整个样本进行迁移聚类，使得样本大幅度减少，从而节省了一定的时间。

Claims (2)

1.一种基于高斯混合模型参数块迁移聚类的SAR图像分割方法，包括如下步骤：

(1)输入待分割SAR图像，确定图像分割类别数C；

(2)对待分割的SAR图像进行3层平稳小波变换，图像像素点总个数为n，对每个像素点提取出10维子带能量特征，构成大小为n×10的输入数据样本集X＝{x_i}，i＝1，2，...，n：

$x_i=\frac{1}{M\×N}\underset{e=1}{\overset{M}{\mathrm{\Σ}}}\underset{u=1}{\overset{N}{\mathrm{\Σ}}}\left|\mathrm{codf}(e,u)\right|$

其中，M×N为利用滑动窗口确定的子带大小，coef(e，u)为平稳小波子带中第e行第u列的系数值；

(3)根据Shi和Malik提出的Normalized Cut算法把待分割的SAR图像划分成M块，将相应的数据样本集X划分成M组样本：X＝{D_r}，r＝1，2，...，M，其中D_r＝{x_b}，b＝1，2，...，s，s表示第r组样本个数大小；

(4)对M组样本数据集X＝{D_r}，r＝1，2，...，M根据组内样本个数大小排序，选择前面的10组为源域X^S＝{D_v}，v＝1，2，...，10，剩余的(M-10)组为目标域X^T＝{D_q}，q＝11，12，...，M；

(5)用步骤(2)提取的样本特征分别对源域X^S＝{D_v}，v＝1，2，...，10的每一组D_v计算均值和协方差，得到10组均值μ_v和协方差σ_v两个类参数，v表示源域组数。

(6)分别对M组样本X＝{D_r}，r＝1，2，...，M中的每一个块D_r＝{x_b}，b＝1，2，...，s求平均得到新样本 $X^{\′}=\left\{{\stackrel{\‾}{D}}_r\right\},r=\mathrm{1,2},...,M,$ 其中 ${\stackrel{\‾}{D}}_r=\underset{b=1}{\overset{s}{\mathrm{\Σ}}}{x}_b/s,$ 则与步骤(4)源域X^S＝{D_v}，v＝1，2，...，10和目标域X^T＝{D_q}，q＝11，12，...，M对应的新源域样本为 $X^{S^{\′}}=\left\{{\stackrel{\‾}{D}}_v\right\},v=\mathrm{1,2},...,10,$ 新目标域样本为 $X^{T^{\′}}=\left\{{\stackrel{\‾}{D}}_q\right\},q=\mathrm{11,12},...,M;$

(7)根据基于欧氏距离的K近邻法对新目标域样本

中的每一个样本

寻找离其最近的新源域样本

并将步骤(5)得到的源域样本类参数均值μ_v和协方差σ_v迁移至新目标域样本

的类参数均值μ_q和协方差σ_q：

μ_q＝μ_v，σ_q＝σ_v；

(8)把迁移过来的那些均值μ_q和协方差σ_q两个类参数当作基于高斯混合模型的EM算法的初始参数，利用该算法对新样本

进行聚类，得到最终图像分割结果。

2.根据权利要求1所述的SAR图像分割方法，其中步骤(8)按照如下步骤进行：

(8a)首先针对新样本

定义归属矩阵z_rj如下：

$z_{\mathrm{rj}}=\left\{\begin{array}{cc}1& {\stackrel{\‾}{D}}_r\∈j\\ 0& {\stackrel{\‾}{D}}_r\∉j\end{array}\right.,r=\mathrm{1,2},...,M;j=\mathrm{1,2},...,C$

其中

是m维向量，M表示样本个数，C表示类别数，若z_rj＝1表示向量

属于第j类，z_rj＝0表示向量

不属于第j类；

(8b)确定初始类参数，包括初始类均值μ_j、类方差σ_j和类混合权重α_j，其中初始类均值μ_j、类方差σ_j为步骤(7)迁移过来的参数均值μ_q和协方差σ_q，即μ_j＝μ_q，σ_j＝σ_q，初始类混合权重α_j为：α_j＝1/C；

(8c)根据上述类参数求得任意一项向量

的概率为：

$p\left({\stackrel{\‾}{D}}_r\right)=\underset{j=1}{\overset{C}{\mathrm{\Σ}}}{\mathrm{\α}}_{j}G({\stackrel{\‾}{D}}_r;{\mathrm{\μ}}_j,{\mathrm{\σ}}_j),i=\mathrm{1,2},...,M$

其中 $G({\stackrel{\‾}{D}}_r,{\mathrm{\μ}}_j,{\mathrm{\σ}}_j)=\frac{1}{{\left(2\mathrm{\π}\right)}^{\frac{m}{2}}{\left|{\mathrm{\σ}}_j\right|}^{\frac{1}{2}}}\exp \{-\frac{1}{2}{({\stackrel{\‾}{D}}_r-{\mathrm{\μ}}_j)}^T{{\mathrm{\σ}}_j}^{-1}({\stackrel{\‾}{D}}_r-{\mathrm{\μ}}_j)\},$ 表示向量

在第j类的高斯分布概率，m为维数；

(8d)利用所求得的概率计算归属矩阵z_rj的期望值E；

$E=E\left(z_{\mathrm{rj}}\right|{\stackrel{\‾}{D}}_r;{\mathrm{\α}}_j,{\mathrm{\μ}}_j,{\mathrm{\σ}}_j)={\mathrm{\α}}_{j}G({\stackrel{\‾}{D}}_r;{\mathrm{\μ}}_j,{\mathrm{\σ}}_j)/p\left({\stackrel{\‾}{D}}_r\right)$

其中

表示期望值E是在类混合权重α_j、类均值μ_j和类方差σ_j条件下求得的向量

的归属矩阵z_rj的期望值；

(8e)用产生的期望值E计算新的类混合权重α′_j、类均值μ′_j和类方差σ′_j：

${{\mathrm{\α}}_j}^{\′}=\frac{1}{n}\underset{r=1}{\overset{M}{\mathrm{\Σ}}}E$

${{\mathrm{\μ}}_j}^{\′}=\underset{r=1}{\overset{M}{\mathrm{\Σ}}}{\mathrm{Ex}}_{\mathrm{rl}}/\underset{r=1}{\overset{M}{\mathrm{\Σ}}}E,l=\mathrm{1,2},...,m$

${{\mathrm{\σ}}_j}^{\′}=\underset{r=1}{\overset{M}{\mathrm{\Σ}}}E({\stackrel{\‾}{D}}_r-{{\mathrm{\μ}}_j}^{\′}){({\stackrel{\‾}{D}}_r-{{\mathrm{\μ}}_j}^{\′})}^T/\underset{r=1}{\overset{M}{\mathrm{\Σ}}}E$

(8f)将归属矩阵z_rj的期望值E最大化：

首先，用所述更新后的类均值μ′_j、类方差σ′_j和类混合权重α′_j重新进行步骤(8c)-(8e)得到归属矩阵z_rj的新期望E′；

然后，反复迭代步骤(8c)-(8e)直至|E′_t+1-E′_t|≤0.000001循环结束，E′_t表示第t次循环得到的z_rj的期望，E′_t+1表示第t+1次循环得到的z_rj的期望；

(8g)利用z_rj的期望E′更新归属矩阵z_rj，得到样本的最终划分：

z_rj＝1当且仅当 $j=\arg \underset{j=1}{\overset{C}{\max }}\left(E^{\′}\right)$

z_rj＝0当且仅当 $j\≠\arg \underset{j=1}{\overset{C}{\max }}\left(E^{\′}\right).$

Images