基于模糊的三马尔可夫场SAR图像分割方法

摘要

本发明公开了一种基于模糊的三马尔可夫场SAR图像分割方法，其实现步骤为：(1)输入待分割的SAR图像；(2)初始化标号场；(3)模糊化标号场；(4)建立附加场；(5)获得模糊联合先验概率；(6)构建后验边缘的分割模型；(7)最大化后验边缘概率去模糊更新标号场分割；(8)判断标号场的变化率是否大于阈值；(9)输出最终分割结果。本发明既可以保持图像中不同区域的边缘的准确性，又可以提高分割结果的区域一致性，具有计算效率和分割精度高的优点。可应用于合成孔径雷达SAR图像分割和SAR图像目标识别。

说明书

技术领域

本发明属于图像处理技术领域，更进一步涉及目标识别领域中基于模糊理论 的三马尔可夫场合成孔径雷达(SAR)图像分割方法。该方法可应用于目标识别图 像信息的获取、SAR图像目标识别，能够对图像的不同区域准确地进行分割。

背景技术

图像分割是通过对图像的分析，将有意义的或者感兴趣的特征部分提取出 来，是成功进行图像分析、理解和描述的关键技术。合成孔径雷达(SAR)是一种 高分辨率成像雷达，它以全天候、全天时、多波段、多极化、可变侧视角以及穿 透能力强等诸多优点被广泛应用于军事和国民经济领域。SAR图像分割的主要 应用分为两种，一种是将感兴趣的目标提取出来，另一种是将图像中不同的地物 按照不同的类别分割成互不重叠的区域。由于SAR图像具有很严重的相干斑噪 声，而且在缺乏先验的条件下对图像分割非常困难，因此SAR图像的处理经常 建立在统计模型的基础上。

马尔可夫随机场模型方法利用观测图像、统计决策和估计理论把图像分割问 题转化为能量函数的最优化问题。它能很好地描述了当前像素与其邻域中像素之 间的相互关系，又能反映图像的随机性和潜在的结构性质。基于马尔可夫随机场 模型的分割方法能较好地分割噪声图像和纹理复杂的SAR图像。当前流行的基 于马尔可夫模型的分割方法有：基于隐马尔可夫树模型的分割方法和基于空间像 素的分割方法等。例如：

西安电子科技大学在其专利申请“基于第二代Bandelet域隐马尔科夫树模型 的图像分割方法”(专利申请号：200810232009.X，公开号：CN101425184A) 中提出了一种基于第二代Bandelet域隐马尔科夫树模型的图像分割方法。这种方 法首先对图像提取图像块进行多尺度变换，训练参数，然后建立树模型利用最大 似然概率分割，最后融合各尺度的分割结果得到分割结果，虽然充分利用了SAR 图像的多分辨信息，但是仍然存在的不足是，这种方法没有深入挖掘图像像素点 之间的空间相关性，所以区域一致性不高，且分割的过程比较复杂。

西北工业大学在其专利申请“一种无监督马尔可夫随机场图像分割方法”(专 利申请号：200710017875.2，公开号：CN101286227A)中提出了一种无监督马 尔可夫随机场图像分割方法。这种方法利用了SAR图像本身像素点之间的空域 相关性，通过新的势能函数提高了抗噪声性能，虽然步骤简单，但是仍然存在的 不足是，没有充分利用图像信息，会造成分割精度较低。

2007年，D.Benboudjema和W.Pieczynski在文献1“Unsupervised Statistical  Segmentation of Nonstationary Images Using Triplet Markov Fields”(IEEE Trans. Pattern Analysis and Machine Intelligence，vol.29，pp.1367-1378，2007)中提出了 基于三马尔可夫随机场模型的非平稳图像的无监督分割方法。该方法将三马尔可 夫场图像分割方法用于非平稳图像中，在模拟和真实的SAR图像中均得到了较 好的分割结果，但是仍然存在的不足是，能量函数的惩罚形式容易造成边缘像素 点的错分割。

发明内容

本发明的目的在于解决上述现有技术的不足，提出了一种基于模糊理论的三 马尔可夫场SAR图像分割方法，利用模糊马尔可夫场的方法对标号场进行模糊 化，并改进了三马尔可夫场中先验概率中能量函数的形式，在保证分割结果区域 一致性的同时，提高了分割结果区域边缘的准确性。

本发明实现上述目的的思路是：先对待分割SAR图像进行均值漂移分割初 始化标号场，模糊化标号场，再对待分割图像进行K-means聚类提取图像的纹理 特征建立附加场，然后利用标号场、附加场、观测场这三个联合场的马尔可夫性 建模图像的联合分布，最后利用贝叶斯最大后验边缘概率准分割图像，逐点更新 像素点标号，得到最终分割结果。

本发明的步骤包括如下：

(1)输入待分割的SAR图像；

(2)初始化标号场；

(3)建立附加场，利用K-means聚类工具，将待分割的SAR图像聚为输入的 类别数，聚类后的矩阵作为附加场；

(4)利用隶属度函数模糊化标号场；

(5)获得联合模糊先验概率：

5a)利用下式计算能量函数：

$W(x,u)=\underset{(s,t)\in C}{\Sigma }{\alpha }^1\left(2\right||x_s-x_t||-1)-\left({\alpha }_a^2\delta (u_s,u_t,a)\right)+{\alpha }_b^2\delta (u_s,u_t,b))\left(\right||x_s-x_t|\left|\right)$

其中，W(x，u)为标号场x和附加场u的能量函数，∑为求和符号，(s，t)是 像素集合中的一对相邻像素点，C为像素点所属的集合，a¹为约束标号场能量的 参数，||x_s-x_t||为标号场的惩罚函数，x_s为当前像素点的标号值，x_t为当前像素 邻域像素点的标号值，为约束附加场能量的参数，a，b为附加场中包含的 两种纹理类别，δ(u_s，u_t，a)、δ(u_s，u_t，b)为附加场的惩罚函数，u_s为当前像素点的 类别，u_t为当前像素邻域像素点的类别；

5b)利用下式获得联合模糊先验概率；

p(x，u)＝γexp[-W(x，u)]

其中，p(x，u)为标号场x和附加场u的联合模糊先验概率，γ为归一化常数， W(x，u)为标号场x和附加场u的能量函数；

(6)构建贝叶斯后验边缘概率的分割模型：

6a)利用下式计算图像中各像素点的似然概率；

$p\left(y_s\right|x_s)=\frac{1}{\Gamma \left(L\right)}\left(\frac{L}{\sigma }\right){y_s}^{L-1}\exp (-\frac{L{y}_s}{\sigma }),$y_s≥0

其中，p(y_s|x_s)为似然概率，y_s为像素点的灰度值，x_s为像素点的标号，s 为像素点，Γ为伽马函数，L为SAR图像的等效视数，σ为各类灰度的均值；

6b)利用统计概率公式计算三马尔可夫场联合概率分布；

6c)利用贝叶斯后验概率方法计算各像素点的后验边缘概率；

(7)利用贝叶斯最大后验边缘概率准则确定每个像素点新的标号，逐点更新 标号场中各像素点的标号，分割图像；

(8)将更新前后标号场中变化的像素点个数和标号场像素点总数的比率作为 终止条件，如果比率大于输入的阈值，返回步骤(4)，否则执行下一步骤；

(9)输出最终分割结果。

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

第一，由于本发明采用均值漂移分割方法初始化图像的标号场，克服了现有 技术多尺度方法所带来的分割结果中图像分割区域边缘不准确的不足，使得本发 明所分割的图像可以保持边缘准确性。

第二，由于本发明采用模糊理论将标号场模糊化的同时，修改了联合分布中 能量函数的形式，削弱了初始标号场错分割对最终分割结果的影响，克服了现有 技术惩罚函数过重造成错分割的不足，使得本发明提高了分割结果的边缘准确 性。

第三，由于本发明提取图像的纹理特征建立图像的附加场，克服了现有技术 利用图像信息不充分的不足，使得本发明所分割的图像提高了图像中同质区域的 区域一致性。

附图说明

图1为本发明的流程图；

图2为本发明使用的三幅原始测试图像和在测试图像上分割两类的结果图；

图3为本发明使用的三幅原始测试图像和在测试图像上分割三类的结果图。

具体实施方式

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

步骤1.输入待分割的SAR图像。

步骤2.初始化标号场。对待分割的SAR图像进行均值漂移分割，输入窗口 宽度为5，输出分割结果图，将分割结果图中灰度值相同的像素标记为一类，得 到的标号矩阵作为初始标号场。

步骤3.建立附加场。

利用K-means聚类工具，将待分割的SAR图像聚为输入的类别数，聚类后 的矩阵作为附加场，本方法将图像聚为两类，表示SAR图像中含有两类不同的 纹理。

步骤4.模糊化标号场。按照下式计算标号场中每一个像素点的隶属度向量；

${\mu }_{\mathrm{sm}}=\frac{N_m}{\underset{m}{\Sigma }{N}_m}$

其中，μ_Sm为像素点s属于标号m类的隶属度，取值范围为[0，1]，N_m为 像素点s的邻域中属于m类的像素点个数，为像素点s的邻域中像素点总数 的求和符号，然后将像素点s在隶属度向量中的最大值作为当前像素点的标号值 x_s，其中，隶属度向量为：{μ_s1，μ_s2，…，μ_sk}，μ_s1为像素点s属于标号为1 的隶属度，μ_s2为像素点s属于标号为2的隶属度，μ_sk为像素点s属于标号为 k类的隶属度，μ为隶属度的值，k为标号场的总类别数。

步骤5.获得联合模糊先验概率。

首先，利用下式计算能量函数：

$W(x,u)=\underset{(s,t)\in C}{\Sigma }{\alpha }^1\left(2\right||x_s-x_t||-1)-\left({\alpha }_a^2\delta (u_s,u_t,a)\right)+{\alpha }_b^2\delta (u_s,u_t,b))\left(\right||x_s-x_t|\left|\right)$

其中，W(x，u)为能量函数，∑为求和符号，(s，t)是像素集合中的一对相邻 像素点，C为像素点所属的集合，a¹为约束标号场能量的参数，本发明取值为 a¹＝2，||x_s-x_t||为标号场的惩罚函数，x_s为当前像素点的标号值，x_t为当前像 素邻域像素点的标号值，为约束附加场能量的参数，本发明取值为 a，b为附加场中包含的两种纹理类别，δ(u_s，u_t，a)、δ(u_s，u_t，b) 为附加场的惩罚函数，u_s为当前像素点的类别，u_t为当前像素邻域像素点的类别； 标号场的惩罚函数按下式计算：

$\left|\right|x_s-x_t\left|\right|=\underset{m=1}{\overset{k}{\Sigma }}|{\mu }_{\mathrm{sm}}-{\mu }_{\mathrm{tm}}|$

其中，||x_s-x_t||为标号场的惩罚函数，x_s为当前像素点的标号值，x_t为当前 像素邻域像素点的标号值，为标号场的标号m取值从标号1到标号k的求 和，k为标号场的总类别数，||为绝对值符号，μ_Sm为像素点s属于标号m类的 隶属度，μ_tm为像素点t属于标号m类的隶属度；

附加场的惩罚函数按下式计算：

其中，δ(u_s，u_t，a)为附加场类别为a的惩罚函数，u_s为当前像素点的类别，u_t为当前像素邻域像素点的类别；

其中，δ(u_s，u_t，b)为附加场类别为b的惩罚函数，u_s为当前像素点的类别，u_t为当前像素邻域像素点的类别。

然后，利用下式获得联合模糊先验概率；

p(x，u)＝γexp[-W(x，u)]

其中，p(x，u)为标号场x和附加场u的联合模糊先验概率，γ为归一化常数， W(x，u)为模糊标号场x和附加场u的能量函数；

步骤6.构建贝叶斯后验边缘概率的分割模型。

利用下式计算图像中各像素点的似然概率：

$p\left(y_s\right|x_s)=\frac{1}{\Gamma \left(L\right)}\left(\frac{L}{\sigma }\right){y_s}^{L-1}\exp (-\frac{L{y}_s}{\sigma }),$y_s≥0

其中，p(y_s|x_s)为似然概率，y_s为像素点的灰度值，x_s为像素点的标号，s 为像素点，Γ为伽马函数，L为SAR图像的等效视数，σ为待估计的各类灰度 的均值，本发明使用最大似然法估计，然后利用下式计算三马尔可夫场联合概率 分布：

p(x，u，y)＝p(x，u)·p(y|x，u)＝γexp[-W(x，u)+∑_s∈Slogp(y_s|x_s)]

其中，p(x，u，y)为标号场x和附加场u和灰度值场y的三马尔可夫场联合概 率分布，p(x，u)为联合先验概率，p(y|x，u)为似然概率，γ为归一化常数，W(x，u) 为联合先验概率中标号场x和附加场u的能量函数，s为像素点，S为像素点所 属的集合，∑为求和符号，p(y_s|x_s)为各像素点的似然概率，y_s为像素点的 灰度值，x_s为像素点的标号；

最后，利用贝叶斯后验概率方法计算各像素点的后验边缘概率：首先，对三 马尔可夫场进行联合概率分布进行吉布斯采样，在像素的8个邻域像素点计算局 部条件概率，依概率随机取值赋给当前像素点，采样20-50次标号场和附加场的 样本，选出样本中每个像素点出现次数最多的标号值，计算各像素点的后验概率 p(x_s，u_s|y)；然后，利用下式计算各像素点的后验边缘概率。

$p\left(x_s\right|y)={\Sigma }_{u_s\in \Lambda }p(x_s,u_s|y)$

其中，p(x_s|y)为后验边缘概率，x_s为像素点的标号，y为灰度值场，s为像 素点，∑为求和符号，u_s为附加场中像素点的类别，Λ为附加场类别所属的集 合。

步骤7.利用下式贝叶斯最大后验边缘概率准则确定每个像素点新的标号， 逐点更新标号场中各像素点的标号，分割图像：

${\hat{x}}_s=\arg \max p_{\omega \in \Omega }(x_s=\omega |y)$

其中，为本次分割后的标号，x_s为逐个选取标号集合中的标号，s为像 素点，arg max为求最大值符号，Pω∈Ω(x_s＝ω|y)为后验边缘概率，ω为标号值， Ω为标号值所属的集合，Ω＝{ω₁，ω₂，…，ω_k}，k为标号场的类别总数。

步骤8.将更新前后标号场中变化的像素点个数和标号场像素点总数的比率 作为终止条件，如果比率大于输入的阈值，返回步骤4，否则执行下一步骤。其 中，阈值为一个较小的常数，取值范围是[10^-6，10^-8]，因为当更新前后标号场中 变化的像素点很少或者没有变化的像素点时，停止迭代，认为已经达到最优，所 以本方法的阈值为10^-8。

步骤9.输出最终分割结果。

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

1仿真条件

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

2仿真内容

图2为本发明仿真试验中使用三幅原始测试SAR图像和现有技术与本发明 在三幅测试图像上分割两类的结果对比图，大小均为256×256，其中，图2(a)、 图2(b)和图2(c)为两类目标的原始测试SAR图像。图2(d)为利用文献1的方法 对图2(a)测试图像的分割结果图，图2(e)为利用文献1的方法对图2(b)测试图 像的分割结果图，图2(f)为利用文献1的方法对图2(c)测试图像的分割结果图。 图2(g)为利用本发明的方法对图2(a)测试图像的分割结果图，图2(h)为利用本 发明的方法对图2(b)测试图像的分割结果图，图2(i)为利用本发明的方法对图 2(c)测试图像的分割结果图。

图3为本发明仿真试验中使用三幅原始测试SAR图像和现有技术与本发明 在三幅测试图像上分割三类的结果对比图，大小均为256×256，其中，图3(a)， 图3(b)和图3(c)为三幅三类目标的原始测试SAR图像。图3(d)为利用文献1的 方法对图3(a)测试图像的分割结果图，图3(e)为利用文献1的方法对图3(b)测 试图像的分割结果图，图3(f)为利用文献1的方法对图3(c)测试图像的分割结 果图。图3(g)为利用本发明的方法对图3(a)测试图像的分割结果图，图3(h)为 利用本发明的方法对图3(b)测试图像的分割结果图，图3(i)为利用本发明的方 法对图3(c)测试图像的分割结果图。

3仿真效果分析

通过上述两类和三类目标的原始测试SAR图像的仿真实验以及分割结果对 比图，可以看出本发明方法在不同SAR图像中都可以进行准确的分割。

图2是三幅对水域和陆地两类目标SAR图像的分割。图2(a)纹理较简单，边 缘较清晰，图2(b)和图2(c)原图中陆地纹理复杂，水陆交接处边缘复杂。由图2(g)、 图2(h)、图2(i)本发明方法的仿真结果可以看到，由于本发明采用模糊随机场的 方法建模能量函数使得图像分割精度得到了明显的提高。在图2(g)中水域从图像 中清晰的分割出来，边缘完整准确，区域一致性好。在图2(h)中，在准确的保持 区域边缘的同时，提高了陆地区域的区域一致性。在图2(i)中，延伸至水域中的 港口部分的边缘和区域一致比图2(f)结果中保持的很好。本方法的三个结果图 2(g)、图2(h)、图2(i)中消除了区域中的杂点，而在图2(d)、图2(e)、图2(f) 文献1方法的结果中同质区域中还存在错分割的杂点。

图3是对三类目标SAR图像的分割。图3(a)是对水域、农田、城区的分割， 图3(b)是对机场跑道与草地、建筑物的分割。图3(c)是对水域、草地和灌木丛的 分割。由图3(g)、图3(h)、图3(i)仿真结果可以看到，虽然对纹理复杂的城市 区域很难分割出来，但是由于本发明采用模糊随机场的方法建模能量函数的同 时，采用Gamma统计分布的概率密度计算三马尔可夫场中的似然概率，提高了 分割结果的区域一致性和抗噪性。在图3(g)和图3(h)中，白色代表的城市区域从 图像中清晰的分割出来，并且区域一致性很高，而在图3(d)和图3(e)中城市区域 并没有完整的分割出来。在图3(i)中，对具有复杂纹理的灌木丛也可以保持边缘 的准确性和区域一致性，比图3(f)分割结果中的区域一致性好很多。