CN103295226B - 一种基于mrf模型的非监督声纳图像分割方法 - Google Patents

Abstract

本发明涉及一种基于MRF模型的非监督声纳图像分割方法，其特征在于：步骤一：对原始声纳图像进行高斯金字塔预处理，得到预处理后的图像；步骤二：计算出预处理后的声纳图像的灰度直方图；步骤三：根据步骤二获得的灰度直方图，计算声纳图像分类及分类个数；步骤四：根据步骤三中的图像分类个数和判别函数,计算MRF分割模型的初始化参数;将初始化参数代入MRF分割模型对声纳图像进行分割。

Description

一种基于MRF模型的非监督声纳图像分割方法

技术领域

本发明涉及一种基于MRF模型的非监督声纳图像分割方法。

背景技术

随着海洋资源开发越来越受世界各国的重视，对水下资源的勘探、搜索和考察也是如火如荼，由于水下环境的复杂性，与光学摄像、激光探测一类传感器相比，声纳是目前应用于水下探测的最有效的传感器。声纳系统在上世纪50年代末诞生以来，应用于军事领域的声纳设备，主要用于障碍物避碰，以及一些军事目标的发现、跟踪和识别；近年来，随着海洋开发活动的开展，声纳设备的应用已不仅仅局限于军事目的，甚至已应用到商业和民用领域，如海底资源开发、石油勘探、海难救助、海底地形地貌图自动绘制以及鱼群探测等等。

但是对水下声纳图像进行目标分割是非常复杂和困难的，它不仅取决于被分割的不同目标区域，还与海底混响噪声、背景区域等有着紧密的联系。对声纳图像分割的目的就是要从复杂的海底混响区域中提取出目标和阴影，并尽量保留图像原始边缘信息，它是图像分析的关键步骤，如何有效地对水下声纳图像进行分割是国内外研究者们研究的热点与难点。

现有技术中，常用的声纳图像分割模型为MRF（马尔科夫）分割模型，但应用声纳图像分割模型时，现有的理论和算法都需人为地确定声纳图像要分类的类别个数，还没有一种完全自动的声纳图像分割模型，现有算法的复杂度较高，耗费时间长，不利于实时性的声纳图像处理与目标识别的要求。

发明内容

本发明目的在于提供一种基于MRF模型的非监督声纳图像分割方法，能够自动确定图像要分类的类别个数。

实现本发明目的技术方案：

一种基于MRF模型的非监督声纳图像分割方法，其特征在于：

步骤一：对原始声纳图像进行高斯金字塔预处理，得到预处理后的图像；

步骤二：计算出预处理后的声纳图像的灰度直方图；

步骤三：根据步骤二获得的灰度直方图，计算声纳图像分类及分类个数，自动确定声纳图像分类及分类个数的模型为：

$\left\{\begin{array}{c}\mathrm{\Γ}(y_i,n)=\hat{f}\left(p\left(r_k\right)\right)\\ i=\mathrm{1,2,3}\\ k=\mathrm{0,1},...,L-1\end{array}\right.$

其中，p(r_k)为图像经过高斯金字塔处理后模型的归一化统计直方图;y₁为图像是否含有阴影区的判别函数，y₁∈{0,1}；y₂为图像是否含有背景区的判别函数，y₂∈{0,1}；y₃为图像是否含有目标区的判别函数，y₃∈{0,1};n为图像分类的个数，n∈{1,2,3}；

步骤四：根据步骤三中的图像分类个数和判别函数,计算MRF分割模型的初始化参数;将初始化参数代入MRF分割模型对声纳图像进行分割。

优选地，步骤三中，声纳图像分类及分类个数的模型通过如下方法计算：

定义像素峰值r_pv为图像统计直方图取得极大值时的灰度级；左边像素谷值r_vvl为图像统计直方图从像素峰值r_pv开始向左统计直方图和变化较小时的灰度级；右边像素谷值r_vvr为图像统计直方图从像素峰值r_pv开始向右统计直方图和变化较小时的灰度级；

并定义判别函数：

$q(x,\mathrm{\θ})=\left\{\begin{array}{cc}0,& x\≤\mathrm{\θ}\\ 1,& x>\mathrm{\θ}\end{array}\right.$

（1）设定迭代条件iter，计数个数count＝0，迭代步数l＝0；

（2）计算像素峰值：

$r_{\mathrm{pv}}=\arg \underset{r_k}{\max }p\left(r_k\right),k=\mathrm{0,1},...,L-1$

并计算灰度级为r_pv时的概率：p^l＝p(r_pv)；

（3）计算灰度区间[r_pv-l,r_pv+l]的概率统计和：

p^l+1＝p([r_pv-l,r_pv+l])；若|p^l+1-p^l|＜iter则count++；

若count大于给定一整数m，则转4，否则，

l＝l+1，并转3；

（4）计算左右边像素谷值：

$\left\{\begin{array}{c}{r}_{\mathrm{vvl}}=r_{\mathrm{pv}}-l\\ r_{\mathrm{vvr}}=r_{\mathrm{pv}}+l\end{array}\right.$

若p([r_vvl,r_vvr])＞1-iter，则令 $\left\{\begin{array}{c}{y}_1=0\\ y_2=1\\ y_3=0\\ n=1\end{array}\right.,$ 算法终止；

（5）计算图像的判别函数y_i及类别个数n，其中λ₁∈(5,20)、λ₂∈(0.5,1)为经验参数：

$\left\{\begin{array}{c}{y}_1=q(\underset{k=0}{\overset{r_{\mathrm{vvl}}}{\mathrm{\Σ}}}p\left(r_k\right),\frac{1-p\left(\right[r_{\mathrm{vvl}},r_{\mathrm{vvr}}\left]\right)}{{\mathrm{\λ}}_1})\\ y_2=q(\underset{k=r_{\mathrm{vvl}}}{\overset{r_{\mathrm{vvr}}}{\mathrm{\Σ}}}p\left(r_k\right),{\mathrm{\λ}}_2)\\ y_3=q(\underset{k=r_{\mathrm{vvr}}}{\overset{L-1}{\mathrm{\Σ}}}p\left(r_k\right),\frac{1-p\left(\right[r_{\mathrm{vvl}},r_{\mathrm{vvr}}\left]\right)}{{\mathrm{\λ}}_1})\\ n=\underset{i=1}{\overset{3}{\mathrm{\Σ}}}{y}_i\end{array}\right.$

优选地，步骤四中，MRF分割模型的初始化参数通过如下方法计算：

令(x,y)为某一图像中像素的坐标，令S_xy表示一确定大小的邻域，其中心在(x,y)，r_s,t表示邻域S_xy中一像素点(s,t)∈S_xy处的像素值，|S_xy|表示邻域S_xy中所有像素点的总个数，则在S_xy中像素的平均值、能量和方差的计算公式为，

$m_{S_{\mathrm{xy}}}=\frac{1}{\left|S_{\mathrm{xy}}\right|}\underset{(s,t)\∈S_{\mathrm{xy}}}{\mathrm{\Σ}}{r}_{s,t}$

$E_{S_{\mathrm{xy}}}=\underset{(s,t)\∈S_{\mathrm{xy}}}{\mathrm{\Σ}}{r_{s,t}}^2$

${{\mathrm{\σ}}_{S_{\mathrm{xy}}}}^2=\frac{1}{\left|S_{\mathrm{xy}}\right|*\left(\right|S_{\mathrm{xy}}|-1)}\left(\underset{(s,t)\∈S_{\mathrm{xy}}}{\mathrm{\Σ}}{(r_{s,t}-m_{S_{\mathrm{xy}}})}^2\right);$

根据能量函数的大小、声纳图像分类的个数和判别函数来得到MRF分割模型的初始化参数。

本发明具有的有益效果：

本发明提出了一种能够自动确定声纳图像分类个数的模型，不需人为确定声纳图像要分类的类别个数，有效提高了声纳图像分割速度，如图7所示，本发明能够达到出色的声纳图像分割效果。本发明还提出了一种结合本发明模型确定马尔科夫随机场模型参数的快速有效的方法。本发明算法简单，速度快，可有效地用于声纳图像处理与目标识别系统中。

附图说明

图1为本发明基于MRF模型的非监督声纳图像分割方法流程图；

图2为原始声纳图像；

图3为图2的直方图；

图4为图2经过高斯金字塔预处理后的图像；

图5为图4的直方图；

图6为MRF初始分割后的声纳图像；

图7是图6去除噪声后的分割结果图像。

具体实施方式

如图1所示，本发明基于MRF模型的非监督声纳图像分割方法实施步骤如下：

步骤一：对原始声纳图像（如图2所示）进行高斯金字塔预处理，得到预处理后的图像（如图4所示）。

步骤二：计算出预处理后的声纳图像的灰度直方图（如图5所示）。灰度级为[0,L-1]范围的数字图像的直方图是离散函数h(r_k)＝n_k，这里r_k是第k级灰度，n_k是图像中的灰度级为r_k的像素个数。一个归一化的直方图由P(r_k)＝n_k/n给出，这里k＝0,1,...,L-1。

步骤三：根据步骤二获得的灰度直方图，计算声纳图像分类及分类个数，自动确定声纳图像分类及分类个数的模型为：

$\left\{\begin{array}{c}\mathrm{\Γ}(y_i,n)=\hat{f}\left(p\left(r_k\right)\right)\\ i=\mathrm{1,2,3}\\ k=\mathrm{0,1},...,L-1\end{array}\right.$

其中，p(r_k)为图像经过高斯金字塔处理后模型的归一化统计直方图;y₁为图像是否含有阴影区的判别函数，y₁∈{0,1}；y₂为图像是否含有背景区的判别函数，y₂∈{0,1}；y₃为图像是否含有目标区的判别函数，y₃∈{0,1};n为图像分类的个数，n∈{1,2,3}；

定义像素峰值r_pv为图像统计直方图取得极大值时的灰度级；左边像素谷值r_vvl为图像统计直方图从像素峰值r_pv开始向左统计直方图和变化较小时的灰度级；右边像素谷值r_vvr为图像统计直方图从像素峰值r_pv开始向右统计直方图和变化较小时的灰度级；

并定义判别函数：

$q(x,\mathrm{\θ})=\left\{\begin{array}{cc}0,& x\≤\mathrm{\θ}\\ 1,& x>\mathrm{\θ}\end{array}\right.$

自动分类模型算法的实现：

（1）设定迭代条件iter，计数个数count＝0，迭代步数l＝0；

（2）计算像素峰值：

$r_{\mathrm{pv}}=\arg \underset{r_k}{\max }p\left(r_k\right),k=\mathrm{0,1},...,L-1$

并计算灰度级为r_pv时的概率：p^l＝p(r_pv)；

（3）计算灰度区间[r_pv-l,r_pv+l]的概率统计和：

p^l+1＝p([r_pv-l,r_pv+l])；若|p^l+1-p^l|＜iter则count++；

若count大于给定一整数m，则转4，否则，

l＝l+1，并转3；

（4）计算左右边像素谷值：

$\left\{\begin{array}{c}{r}_{\mathrm{vvl}}=r_{\mathrm{pv}}-l\\ r_{\mathrm{vvr}}=r_{\mathrm{pv}}+l\end{array}\right.$

若p([r_vvl,r_vvr])＞1-iter，则令 $\left\{\begin{array}{c}{y}_1=0\\ y_2=1\\ y_3=0\\ n=1\end{array}\right.,$ 算法终止；

（5）计算图像的判别函数y_i及类别个数n，其中λ₁∈(5,20)、λ₂∈(0.5,1)为经验参数：

$\left\{\begin{array}{c}{y}_1=q(\underset{k=0}{\overset{r_{\mathrm{vvl}}}{\mathrm{\Σ}}}p\left(r_k\right),\frac{1-p\left(\right[r_{\mathrm{vvl}},r_{\mathrm{vvr}}\left]\right)}{{\mathrm{\λ}}_1})\\ y_2=q(\underset{k=r_{\mathrm{vvl}}}{\overset{r_{\mathrm{vvr}}}{\mathrm{\Σ}}}p\left(r_k\right),{\mathrm{\λ}}_2)\\ y_3=q(\underset{k=r_{\mathrm{vvr}}}{\overset{L-1}{\mathrm{\Σ}}}p\left(r_k\right),\frac{1-p\left(\right[r_{\mathrm{vvl}},r_{\mathrm{vvr}}\left]\right)}{{\mathrm{\λ}}_1})\\ n=\underset{i=1}{\overset{3}{\mathrm{\Σ}}}{y}_i\end{array}\right.$

步骤四：根据步骤三中的图像分类个数和判别函数,计算MRF分割模型的初始化参数;将初始化参数代入MRF分割模型对声纳图像进行分割。

MRF分割模型的初始化参数通过如下方法计算：

令(x,y)为某一图像中像素的坐标，令S_xy表示一确定大小的邻域，其中心在(x,y)，r_s,t表示邻域S_xy中一像素点(s,t)∈S_xy处的像素值，|S_xy|表示邻域S_xy中所有像素点的总个数，则在S_xy中像素的平均值、能量和方差的计算公式为，

$m_{S_{\mathrm{xy}}}=\frac{1}{\left|S_{\mathrm{xy}}\right|}\underset{(s,t)\∈S_{\mathrm{xy}}}{\mathrm{\Σ}}{r}_{s,t}$

$E_{S_{\mathrm{xy}}}=\underset{(s,t)\∈S_{\mathrm{xy}}}{\mathrm{\Σ}}{r_{s,t}}^2$

${{\mathrm{\σ}}_{S_{\mathrm{xy}}}}^2=\frac{1}{\left|S_{\mathrm{xy}}\right|*\left(\right|S_{\mathrm{xy}}|-1)}\left(\underset{(s,t)\∈S_{\mathrm{xy}}}{\mathrm{\Σ}}{(r_{s,t}-m_{S_{\mathrm{xy}}})}^2\right);$

局部能量极化也就是把图像分成很多大小相同的区域，然后分别计算源图像中各个区域的均值、能量和方差。根据能量函数的大小以及图像分类的个数和判别函数得到MRF分割模型的初始化参数。

所说的MRF分割模型为现有技术，如下所述：

假定我们观测到的图像数据为F，图像上所有像素点的集合记为S。图像的分割问题是要求解的问题是满足最大后验概率准则，对每个像素的分类标号（标号场），记为ω。这样由Bayes后验概率准则：

$P\left(\mathrm{\ω}\right|F)=\frac{1}{P\left(F\right)}P\left(\mathrm{\ω}\right)P\left(F\right|\mathrm{\ω})$

其中P(F)为观测数据的先验分布，当数据给定后为常数，所以不参与计算过程，可以不予以考虑。

P(ω)是标号场的先验的联合的Gibbs分布，即满足马尔可夫性。假设C表示S所有的集簇，c表示C中的元素，U₂(ω)为能量函数，V_c(ω)是与集簇相关的势函数，那么：

$P\left(\mathrm{\ω}\right)\∝\exp \left(U_2\left(\mathrm{\ω}\right)\right)=\exp (\underset{c\∈C}{\mathrm{\Σ}}-V_c\left({\mathrm{\ω}}_c\right))$

P(F|ω)是似然概率，在很多情况下，假定为各个位置的像素是独立同分布的，即满足：

$P\left(F\right|\mathrm{\ω})=\underset{s\∈S}{\mathrm{\Π}}P\left(F_s\right|{\mathrm{\ω}}_s)$

当我们假定每个P(F_s|ω_s)是高斯分布时，每个类的类参数都是由两个参数唯一确定该分布，即为λ和σ。将Bayes后验概率准则：

$P\left(\mathrm{\ω}\right|F)=\frac{1}{P\left(F\right)}P\left(\mathrm{\ω}\right)P\left(F\right|\mathrm{\ω})$

取对数，得到的目标函数为lnP(ω)+lnP(F|ω)，即我们要求的是使该表达式最大的时候ω的估计，将似然函数和先验Gibbs分布的表达式带入，这时我们可以形成MRF-MAP下的目标函数的最优解问题：

则我们根据声纳图像的性质，选取一阶的邻域系统，并且势函数为ISING模型的势函数：

$V\left({\mathrm{\ω}}_c\right)=V_{\{s,r\}}({\mathrm{\ω}}_c,{\mathrm{\ω}}_r)=\left\{\begin{array}{ccc}-\mathrm{\β}& \mathrm{if}& {\mathrm{\ω}}_c={\mathrm{\ω}}_r\\ +\mathrm{\β}& \mathrm{if}& {\mathrm{\ω}}_c\≠{\mathrm{\ω}}_r\end{array}\right.$

其中β为耦合系数。对每个像素点选取P(F_s|ω_s)服从高斯分布，那么取对数后其函数形式为：

lnP(F|ω)＝U₁(ω,F)

$=\underset{s\∈S}{\mathrm{\Σ}}\left(\ln (\sqrt{2\mathrm{\π}}{\mathrm{\σ}}_{{\mathrm{\ω}}_s}+\frac{{(f_s-{\mathrm{\μ}}_{{\mathrm{\ω}}_s})}^2}{2{{\mathrm{\σ}}_{{\mathrm{\ω}}_s}}^2})\right)$

取lnP(ω)+lnP(F|ω)为目标函数，得到的分割结果为：

$\mathrm{Arg}\left(\underset{\mathrm{\ω}}{\max }(U_1(\mathrm{\ω},F)+U_2\left(\mathrm{\ω}\right))\right)$

其中U₂(ω)为能量函数，且

Claims (3)

1.一种基于MRF模型的非监督声纳图像分割方法，其特征在于：

步骤一：对原始声纳图像进行高斯金字塔预处理，得到预处理后的图像；

步骤二：计算出预处理后的声纳图像的灰度直方图；

步骤三：根据步骤二获得的灰度直方图，计算声纳图像分类及分类个数，自动确定声纳图像分类及分类个数的模型为：

$\left\{\begin{array}{c}\mathrm{\Γ}(y_i,n)=\hat{f}\left(p\right(r_k\left)\right)\\ i=1,2,3\\ k=0,1,...,L-1\end{array}\right.$

其中，p(r_k)为图像经过高斯金字塔预处理后模型的归一化统计直方图；y₁为图像是否含有阴影区的判别函数，y₁∈{0,1}；y₂为图像是否含有背景区的判别函数，y₂∈{0,1}；y₃为图像是否含有目标区的判别函数，y₃∈{0,1}；n为图像分类的个数，n∈{1,2,3}；

步骤四：根据步骤三中的图像分类个数和判别函数,计算MRF分割模型的初始化参数；将初始化参数代入MRF分割模型对声纳图像进行分割。

2.根据权利要求1所述的基于MRF模型的非监督声纳图像分割方法，其特征在于：步骤三中，声纳图像分类及分类个数的模型通过如下方法计算：

定义像素峰值r_pv为图像统计直方图取得极大值时的灰度级；左边像素谷值r_vvl为图像统计直方图从像素峰值r_pv开始向左统计直方图和变化较小时的灰度级；右边像素谷值r_vvr为图像统计直方图从像素峰值r_pv开始向右统计直方图和变化较小时的灰度级；

并定义判别函数：

$q(x,\mathrm{\θ})=\left\{\begin{array}{cc}0,& x\≤\mathrm{\θ}\\ 1,& x>\mathrm{\θ}\end{array}\right.$

(1)设定迭代条件iter，计数个数count＝0，迭代步数l＝0；

(2)计算像素峰值：

k＝0,1,...,L-1

并计算灰度级为r_pv时的概率：p^l＝p(r_pv)；

(3)计算灰度区间[r_pv-l,r_pv+l]的概率统计和：

p^l+1＝p([r_pv-l,r_pv+l])；若|p^l+1-p^l|＜iter则count++；

若count大于给定一整数m，则转步骤(4)，否则，

l＝l+1，重复执行步骤(3)；

(4)计算左右边像素谷值：

$\left\{\begin{array}{c}{r}_{vvl}=r_{pv}-l\\ r_{vvr}=r_{pv}+l\end{array}\right.$

若p([r_vvl,r_vvr])＞1-iter，则令 $\left\{\begin{array}{c}{y}_1=0\\ y_2=1\\ y_3=0\\ n=1\end{array}\right.,$ 算法终止；

(5)计算图像的判别函数y_i及类别个数n，其中λ₁∈(5,20)、λ₂∈(0.5,1)为经验参数：

$\left\{\begin{array}{c}{y}_1=q\left(\underset{k=0}{\overset{r_{vvl}}{\Σ}}p\right(r_k),\frac{1-p(\[r_{vvl},r_{vvr}\])}{{\mathrm{\λ}}_1})\\ y_2=q\left(\underset{k=r_{vvl}}{\overset{r_{vvr}}{\Σ}}p\right(r_k),{\mathrm{\λ}}_2)\\ y_3=q\left(\underset{k=r_{vvr}}{\overset{L-1}{\Σ}}p\right(r_k),\frac{1-p(\[r_{vvl},r_{vvr}\])}{{\mathrm{\λ}}_1})\\ n=\underset{i=1}{\overset{3}{\Σ}}{y}_i\end{array}\right..$

3.根据权利要求2所述的基于MRF模型的非监督声纳图像分割方法，其特征在于：步骤四中，MRF分割模型的初始化参数通过如下方法计算：

令(x,y)为某一图像中像素的坐标，令S_xy表示一确定大小的邻域，其中心在(x,y)，r_s,t表示邻域S_xy中一像素点(s,t)∈S_xy处的像素值，|S_xy|表示邻域S_xy中所有像素点的总个数，则在S_xy中像素的平均值、能量和方差的计算公式为，

$m_{S_{xy}}=\frac{1}{\left|S_{xy}\right|}\underset{(s,t)\∈S_{xy}}{\Σ}{r}_{s,t}$

$E_{S_{xy}}=\underset{(s,t)\∈S_{xy}}{\Σ}{r_{s,t}}^2$

${{\mathrm{\σ}}_{S_{xy}}}^2=\frac{1}{\left|S_{xy}\right|*\left(\right|S_{xy}|-1)}\left(\underset{(s,t)\∈S_{xy}}{\Σ}{(r_{s,t}-m_{S_{xy}})}^2\right);$

根据能量函数的大小、声纳图像分类的个数和判别函数来得到MRF分割模型的初始化参数。