CN107545573A - 一种基于模态mrf随机场的水下前视声纳图像分割方法 - Google Patents

Abstract

本发明涉及一种基于模态MRF随机场的水下前视声纳图像分割方法。本发明的核心是利用模态Potts网络对图像进行建模和分割。具体步骤包括：初始化Potts态；用EM算法估计类条件参数；每个Potts单元根据热浴动力进行更新,Potts单元的内场满足Markov依赖,Potts的外场是像素的类条件概率的对数；更新Potts单元之间的交互参数，以及Potts单元与外场之间的交互参数；根据参数的收敛程度判断是否退出。本发明采用基于热浴动力的Potts网络进行图像分割，符合大脑的认知机制；本发明在去除斑块噪声的同时能几乎完美的保留目标的边缘信息。

Description

一种基于模态MRF随机场的水下前视声纳图像分割方法

技术领域

本发明属于水下声纳图像处理领域，具体说是一种基于模态Markov随机场的前视声纳图像分割方法。

背景技术

水下机器人是探索水下未知环境、执行水下特定任务的载体。由于光电信号在海水中的作用距离十分有限，声纳就成为了水下机器人感知水下环境的基本必备工具。声纳图像在海洋环境监测、海底资源勘探、水声情报搜集和协助水下工程作业等方面发挥着非常重要的角色。而其中，前视声纳是近距离调查和作业的首选工具。

然而，受到水体的反射和声波的散射等因素的干扰，前视声纳图像中包含大量的斑块噪声。因此，前视声纳图像处理的最大难点在于如何在去除斑块噪声的同时尽可能的保留目标区域的完整轮廓信息。

Markov随机场模型(Markov random field,MRF)和基于均值场近似的期望最大算法(Expectation-Maximization with Mean-Field approximation,EMMF)算法在声纳图像中取得了不错的效果。但是，却存在如下的缺点：(1)MRF模型处理斑块噪声的能力有限，在绝大多数情况下都需要进行形态学后处理去除先验概率较大的噪声。(2)EMMF算法虽然具有较好的平滑效果，但是却存在过平滑的现象，会丢失大量的轮廓信息。

发明内容

本发明目的是提供一种水下前视声纳图像分割方法，在去除斑块噪声的同时尽可能的保留目标区域的轮廓信息。

本发明为实现上述目的所采用的技术方案是：一种基于模态MRF随机场的水下前视声纳图像分割方法，包括以下步骤：

1)采用模糊C均值聚类算法对声纳图像进行聚类得到隶属度矩阵，并用隶属度矩阵初始化当前图像每个Potts单元的激活向量Q为Potts单元的状态数即类的个数，i为像素序号；

2)设置交互参数的当前学习速率η₀为设定的初始学习速率，t为当前迭代次数；

3)采用期望最大算法估计条件参数Θ_k～(π_k；μ_k,∑_k)；

4)根据热浴动力更新Potts单元；

5)根据热平衡原理更新交互参数；返回步骤3)，直到交互参数的变化量Δβ≤ε且Δλ≤ε为止，ε为常数，此时得到的图像为分割后的图像。

所述采用期望最大算法估计条件参数Θ_k～(π_k；μ_k,∑_k)包括以下步骤：

(1)计算每个像素x_{i＝1,...,H×W}属于每个类k＝1,...,Q的概率：

其中，π_k为先验概率，μ_k为均值，∑_k为方差；f(x_i|Θ_k)是第i个Potts单元属于第k个类的概率，服从高斯分布；H和W分别为图像的高度和宽度；

(2)根据似然概率最大求取条件参数：

像素的个数N＝HW。

所述根据热浴动力更新Potts单元通过下式得到：

其中，局部场包括内场和外场

内场和外场分别为

其中，λ是Potts单元的外场的权重，β是Potts单元之间的交互权重；表示第i个Potts单元位于状态k时的概率，表示第j个Potts单元位于状态l时的概率；c_i表示第i个Potts单元的邻域，j∈c_i表示第i个Potts单元的邻域中的第j个Potts单元；f(x_i|Θ_k)是第i个Potts单元属于第k个类的概率，δ_kl是Kronecker函数，且

所述根据热平衡原理更新交互参数包括以下步骤：

分别求取Potts单元之间的交互权重β和Potts单元的外场的权重λ的变化量：

将两个变化量分别与上次的β、λ相加得到当前的交互参数。

本发明具有以下有益效果及优点：

1.本发明采用基于热浴动力的Potts网络进行图像分割，符合大脑的认知机制。

2.本发明指出外场和类条件概率的对数成正比，并根据Markov链采样的原理得出交互参数的学习公式。从图5可以看出算法能够快速的收敛。

3.本发明在去除斑块噪声的同时能几乎完美的保留目标的边缘信息，如图6所示。这表明，本发明对于水下前视声纳图像具有较好的分割效果。当然，本发明也可用于分割其它包含强斑块噪声的图像。

附图说明

图1为本发明的前视声纳图像；

图2为FCM分割的结果图像；

图3为Potts网络的内场和外场示意图；

图4为系统的平衡态结果图；

图5为交互参数学习的过程曲线图；

图6为本发明分割的结果图像；

图7为MRF分割的结果图像；

图8为EMMF分割的结果图像；

图9为本发明的方法流程图。

具体实施方式

下面结合附图对本发明做进一步的详细说明。

如图9所示，一种基于模态Markov随机场的前视声纳图像分割方法，包括如下步骤：

算法输入：前视声纳图像X＝{x_i}_{i＝1,...,H×W}(如图1所示)，其中x_i表示第i个像素的灰度值，H和W分别为图像的高度和宽度。

算法输出：类别矩阵L＝{l_i}_{i＝1,...,H×W}；l_i表示第i个像素的类别。

步骤一：设定Potts单元的状态数Q，并用变量q表示任意状态，即q＝1,...,Q。对于前视声纳图像，一般将Q设定为2或3。

步骤二：Potts网络初始化。利用模糊C均值聚类(Fuzzy C-Means，FCM)算法对声纳图像进行聚类，其中聚类的个数即为Potts单元的状态数。聚类的结果是一个维度为3的隶属度矩阵其中表示第i个像素对应的隶属度向量，且向量的长度为Q。对于每个像素i，用隶属度向量初始化对应的Potts单元的激活向量为了可视化隶属度矩阵，可以将其转化为一个二维矩阵。即，对每个像素所对应的隶属度向量取最大值，然后将最大值所对应的类作为灰度(参考步骤八)，这样可以得到一个灰度级为Q的二维矩阵(即初步的分割结果)，如图2所示。

步骤三：设置交互权重(λ和β)的当前学习速率其中，t表示第t次迭代，η₀为初始学习速率，可设为1。λ是Potts单元的外场的权重，β是Potts单元内场的交互权重，二者使用相同的学习速率。

步骤四：条件参数估计。声纳图像每个类的像素都近似服从Gauss分布；第k个Gauss分布的参数集为Θ_k＝(π_k；μ_k,Σ_k)。用f(x_i|Θ_k)表示第i个像素属于第k个Gauss分布的条件概率，

本步骤估计每个类的先验概率π_k、均值μ_k和方差Σ_k。采用期望最大算法(Expectation-Maxmization,EM)估计条件参数。

(1)期望步：计算每个像素x_{i＝1,...,H×W}属于每个类k＝1,...,Q的概率；

(2)最大化：根据似然最大(Maximum Likelihood,ML)求取条件参数。

N＝H*W为图像中像素的个数。

步骤五：Potts态更新。根据热浴动力更新每个Potts网络的状态：

局部场包括内场和外场k＝1,...,Q，

内场和外场之间的关系如图3所示，且分别为：

其中，λ是Potts单元的外场的权重，β是Potts单元之间的交互权重，即Markov随机场中的热温度；表示第i个Potts单元位于状态k时的概率，表示第j个Potts单元位于状态l时的概率；c_i表示第i个Potts单元的邻域，j∈c_i表示第i个Potts单元的邻域中的第j个Potts单元；f(x_i|Θ_k)是第i个Potts单元属于第k个类的概率。δ_kl是Kronecker函数，且

步骤六：交互参数学习。上一步等价于对外场进行更新。在新的外场下Potts系统要想实现平衡，需要根据平衡态对系统的交互参数λ和β作出调整。为了学习交互参数，本发明设计了一个目标函数即Potts单元的状态更新之后与更新之前的能量差的平方。也就是说，本发明中的热平衡态定义为系统能量不再变化时Potts网络的状态(如图4所示)。根据目标函数最小化和梯度下降的原理可以推导出内场交互和外场交互参数的变化量为：

步骤七：决定算法是否满足终止条件。如果内场交互参数变化量Δβ≤ε且外场交互参数变化量Δλ≤ε，则退出循环，转步骤八；否则转步骤三。其中，ε是某个非常小的常数，可设0～0.1，如ε＝0.05。

步骤八：计算每个像素的类别，

输出最终的类别矩阵

图5和图6分别是交互参数的学习过程和最后的分割结果。作为对本发明的优势的说明和支持，图7和图8给出了Markov随机场和基于均值场的EM算法的分割结果。

Claims (4)

1.一种基于模态MRF随机场的水下前视声纳图像分割方法，其特征在于包括以下步骤：

1)采用模糊C均值聚类算法对声纳图像进行聚类得到隶属度矩阵，并用隶属度矩阵初始化当前图像每个Potts单元的激活向量Q为Potts单元的状态数即类的个数，i为像素序号；

2)设置交互参数的当前学习速率η₀为设定的初始学习速率，t为当前迭代次数；

3)采用期望最大算法估计条件参数Θ_k～(π_k；μ_k，∑_k)；

4)根据热浴动力更新Potts单元；

5)根据热平衡原理更新交互参数；返回步骤3)，直到交互参数的变化量Δβ≤ε且Δλ≤ε为止，ε为常数，此时得到的图像为分割后的图像。

2.根据权利要求1所述的一种基于模态MRF随机场的水下前视声纳图像分割方法，其特征在于所述采用期望最大算法估计条件参数Θ_k～(π_k；μ_k，∑_k)包括以下步骤：

(1)计算每个像素x_{i＝1，...，H×W}属于每个类k＝1，...，Q的概率：

<mrow> <msubsup> <mi>&amp;alpha;</mi> <mi>i</mi> <mi>k</mi> </msubsup> <mo>=</mo> <mfrac> <mrow> <msub> <mi>&amp;pi;</mi> <mi>k</mi> </msub> <mi>f</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>|</mo> <msub> <mi>&amp;mu;</mi> <mi>k</mi> </msub> <mo>,</mo> <msub> <mi>&amp;Sigma;</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> </mrow> <mrow> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>q</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>Q</mi> </munderover> <msub> <mi>&amp;pi;</mi> <mi>q</mi> </msub> <mi>f</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>|</mo> <msub> <mi>&amp;mu;</mi> <mi>q</mi> </msub> <mo>,</mo> <msub> <mi>&amp;Sigma;</mi> <mi>q</mi> </msub> <mo>)</mo> </mrow> </mrow> </mfrac> </mrow>

其中，π_k为先验概率，μ_k为均值，∑_k为方差；f(x_i|Θ_k)是第i个Potts单元属于第k个类的概率，服从高斯分布；H和W分别为图像的高度和宽度；

(2)根据似然概率最大求取条件参数：

<mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>&amp;pi;</mi> <mi>k</mi> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <mi>N</mi> </mfrac> <munder> <mi>&amp;Sigma;</mi> <mi>i</mi> </munder> <msubsup> <mi>&amp;alpha;</mi> <mi>i</mi> <mi>k</mi> </msubsup> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>&amp;mu;</mi> <mi>k</mi> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <msub> <mi>N&amp;pi;</mi> <mi>k</mi> </msub> </mrow> </mfrac> <munder> <mi>&amp;Sigma;</mi> <mi>i</mi> </munder> <msubsup> <mi>&amp;alpha;</mi> <mi>i</mi> <mi>k</mi> </msubsup> <msub> <mi>x</mi> <mi>i</mi> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>&amp;Sigma;</mi> <mi>k</mi> </msub> <mo>=</mo> <mfrac> <mn>1</mn> <mrow> <msub> <mi>N&amp;pi;</mi> <mi>k</mi> </msub> </mrow> </mfrac> <munder> <mi>&amp;Sigma;</mi> <mi>i</mi> </munder> <msubsup> <mi>&amp;alpha;</mi> <mi>i</mi> <mi>k</mi> </msubsup> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>-</mo> <msub> <mi>&amp;mu;</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> <msup> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>-</mo> <msub> <mi>&amp;mu;</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> <mi>T</mi> </msup> </mrow> </mtd> </mtr> </mtable> </mfenced>

像素的个数N＝HW。

3.根据权利要求1所述的一种基于模态MRF随机场的水下前视声纳图像分割方法，其特征在于所述根据热浴动力更新Potts单元通过下式得到：

<mrow> <msubsup> <mi>&amp;sigma;</mi> <mi>i</mi> <mi>k</mi> </msubsup> <mo>=</mo> <mfrac> <msup> <mi>e</mi> <msubsup> <mi>h</mi> <mi>i</mi> <mi>k</mi> </msubsup> </msup> <mrow> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>Q</mi> </munderover> <msup> <mi>e</mi> <msubsup> <mi>h</mi> <mi>i</mi> <mi>l</mi> </msubsup> </msup> </mrow> </mfrac> </mrow>

其中，局部场包括内场和外场

<mrow> <msubsup> <mi>h</mi> <mi>i</mi> <mi>k</mi> </msubsup> <mo>=</mo> <msubsup> <mi>I</mi> <mi>i</mi> <mi>k</mi> </msubsup> <mo>+</mo> <msubsup> <mi>&amp;lambda;O</mi> <mi>i</mi> <mi>k</mi> </msubsup> </mrow>

内场和外场分别为

<mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <msubsup> <mi>I</mi> <mi>i</mi> <mi>k</mi> </msubsup> <mo>=</mo> <mi>&amp;beta;</mi> <mstyle> <munder> <mo>&amp;Sigma;</mo> <mrow> <mi>j</mi> <mo>&amp;Element;</mo> <msub> <mi>c</mi> <mi>i</mi> </msub> </mrow> </munder> </mstyle> <mstyle> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>Q</mi> </munderover> </mstyle> <msub> <mi>&amp;delta;</mi> <mrow> <mi>k</mi> <mi>l</mi> </mrow> </msub> <msubsup> <mi>&amp;sigma;</mi> <mi>i</mi> <mi>k</mi> </msubsup> <msubsup> <mi>&amp;sigma;</mi> <mi>j</mi> <mi>l</mi> </msubsup> </mtd> </mtr> <mtr> <mtd> <mrow> <msubsup> <mi>O</mi> <mi>i</mi> <mi>k</mi> </msubsup> <mo>=</mo> <mi>l</mi> <mi>o</mi> <mi>g</mi> <mi> </mi> <mi>f</mi> <mrow> <mo>(</mo> <msub> <mi>x</mi> <mi>i</mi> </msub> <mo>|</mo> <msub> <mi>&amp;Theta;</mi> <mi>k</mi> </msub> <mo>)</mo> </mrow> </mrow> </mtd> </mtr> </mtable> </mfenced>

其中，λ是Potts单元的外场的权重，β是Potts单元之间的交互权重；表示第i个Potts单元位于状态k时的概率，表示第j个Potts单元位于状态l时的概率；c_i表示第i个Potts单元的邻域，j∈c_i表示第i个Potts单元的邻域中的第j个Potts单元；f(x_i|Θ_k)是第i个Potts单元属于第k个类的概率，δ_kl是Kronecker函数，且

<mrow> <msub> <mi>&amp;delta;</mi> <mrow> <mi>k</mi> <mi>l</mi> </mrow> </msub> <mo>=</mo> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <mn>1</mn> <mo>,</mo> </mrow> </mtd> <mtd> <mrow> <mi>i</mi> <mi>f</mi> </mrow> </mtd> <mtd> <mrow> <mi>k</mi> <mo>=</mo> <mi>l</mi> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mn>0</mn> <mo>,</mo> </mrow> </mtd> <mtd> <mrow> <mi>i</mi> <mi>f</mi> </mrow> </mtd> <mtd> <mrow> <mi>k</mi> <mo>&amp;NotEqual;</mo> <mi>l</mi> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>.</mo> </mrow>

4.根据权利要求1所述的一种基于模态MRF随机场的水下前视声纳图像分割方法，其特征在于所述根据热平衡原理更新交互参数包括以下步骤：

分别求取Potts单元之间的交互权重β和Potts单元的外场的权重λ的变化量：

<mrow> <mi>&amp;Delta;</mi> <mi>&amp;beta;</mi> <mo>=</mo> <mi>&amp;eta;</mi> <munder> <mo>&amp;Sigma;</mo> <mi>i</mi> </munder> <munder> <mo>&amp;Sigma;</mo> <mi>k</mi> </munder> <mo>&amp;lsqb;</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>i</mi> <mi>k</mi> </msubsup> <mrow> <mo>(</mo> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>-</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>i</mi> <mi>k</mi> </msubsup> <mo>&amp;rsqb;</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>i</mi> <mi>k</mi> </msubsup> <mrow> <mo>(</mo> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>&amp;lsqb;</mo> <msubsup> <mi>I</mi> <mi>i</mi> <mi>k</mi> </msubsup> <mo>-</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>Q</mi> </munderover> <msubsup> <mi>I</mi> <mi>i</mi> <mi>l</mi> </msubsup> <msubsup> <mi>&amp;sigma;</mi> <mi>i</mi> <mi>l</mi> </msubsup> <mrow> <mo>(</mo> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow>

<mrow> <mi>&amp;Delta;</mi> <mi>&amp;lambda;</mi> <mo>=</mo> <mi>&amp;eta;</mi> <munder> <mo>&amp;Sigma;</mo> <mi>i</mi> </munder> <munder> <mo>&amp;Sigma;</mo> <mi>k</mi> </munder> <mo>&amp;lsqb;</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>i</mi> <mi>k</mi> </msubsup> <mrow> <mo>(</mo> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>-</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>i</mi> <mi>k</mi> </msubsup> <mo>&amp;rsqb;</mo> <msubsup> <mi>&amp;sigma;</mi> <mi>i</mi> <mi>k</mi> </msubsup> <mrow> <mo>(</mo> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>&amp;lsqb;</mo> <msubsup> <mi>O</mi> <mi>i</mi> <mi>k</mi> </msubsup> <mo>-</mo> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>Q</mi> </munderover> <msubsup> <mi>O</mi> <mi>i</mi> <mi>l</mi> </msubsup> <msubsup> <mi>&amp;sigma;</mi> <mi>i</mi> <mi>l</mi> </msubsup> <mrow> <mo>(</mo> <mi>t</mi> <mo>+</mo> <mn>1</mn> <mo>)</mo> </mrow> <mo>&amp;rsqb;</mo> </mrow>

将两个变化量分别与上次的β、λ相加得到当前的交互参数。