多视SAR图像分割方法及装置与流程

本发明涉及图像分割技术领域，具体而言，涉及一种多视SAR图像分割方法及装置。

背景技术：

合成孔径雷达(Synthetic Aperture Radar,SAR)是一种高分辨率成像雷达，通过接收目标散射的电磁波并转换成图像以记录地物形态，其特有的成像机制所导致的固有斑点噪声给图像分割带来了巨大的困难。尽管多视技术能够减少部分噪声，但是，在实际应用中多视SAR图像仍然存在大量斑点噪声，因此，多视SAR图像分割方法的抗噪性及准确性一直是研究的热点问题。

目前，多视SAR图像分割的方法主要有：阈值法、边界法、聚类法及统计模型法等。其中，应用最广泛的为统计法，其通常采用混合模型刻画图像的复杂分布情况。混合模型中最常用的为高斯混合模型(Gaussian Mixture Model,GMM)，其假设图像中像素的灰度值服从高斯分布，但是，传统GMM中的权重系数为用向量表示的仅与聚类相关的单权重，且SAR图像服从Gamma分布，上述问题均会导致传统GMM对SAR图像建模不够准确。

针对上述多视SAR图像分割方法抗噪性能差及分割结果不理想的问题，目前尚未提出有效的解决方案。

技术实现要素：

有鉴于此，本发明的目的在于提供一种多视SAR图像分割方法及装置，能够提高图像分割的抗噪性能且增强分割精度。

为了实现上述目的，本发明实施例采用的技术方案如下：

第一方面，本发明实施例提供了一种多视SAR图像分割方法，包括：

读取待分割多视SAR图像；上述待分割多视SAR图像使用特征场z表示：

z＝{z_i(x_i,y_i):i＝1,…,n}

其中，i为像素索引，n为总像素数，z_i为像素i的强度，(x_i，y_i)∈D为像素i的格点位置，D为图像域；

初始化双权重w：

w_i＝(w_il:l＝1,…,k)

其中，w_il为类属l包含像素i的类属权重，满足k为类数；

重复执行下述步骤：计算Gamma分布尺度参数β：

其中，α为Gamma分布形态参数；

计算上述图像的或然率p(z|w)；

其中p(z_i|w_il)为以双权重Gamma混合模型定义的z_i的概率密度，如下：

计算上述双权重w的分布函数p(w)：

其中，A为归一化系数，η为邻域作用系数，N_i为以(x_i,y_i)为中心的8邻域像素集合，且满足i’∈N_i,i’≠i；

计算品质函数L，上述品质函数L为或然率p(z|w)与分布函数p(w)联合概率分布函数的对数函数；

根据梯度法更新上述双权重w：

w^(t+1)＝w^(t)+ξΔw^(t)

其中，t为迭代次数，ξ为步长，Δw^(t)为梯度，如下：

将更新后的双权重w代入品质函数L；

直到|L^(t+1)-L^(t)|小于预设的阈值ε时，停止执行上述步骤，根据当前的双权重w确定上述图像中各个像素所属的类别；

按照各个像素所属的类别输出分割结果。

结合第一方面，本发明实施例提供了第一方面的第一种可能的实施方式，其中，计算品质函数L包括：

计算或然率p(z|w)与分布函数p(w)联合概率分布函数p(z,w)：

p(z,w)＝p(z|w)p(w)；

对该联合概率分布函数取对数：

L(w)＝log p(z,w)＝log p(z|w)+log p(w)。

结合第一方面，本发明实施例提供了第一方面的第二种可能的实施方式，其中，确定上述像素所属的类别包括：

计算双权重的最大值：

c_i＝arg max{w_il,l＝1,...,k}；

其中，c_i为第i个像素所属类别的标号；

将上述最大值作为上述像素所属的类别。

结合第一方面，本发明实施例提供了第一方面的第三种可能的实施方式，其中，按照各个像素所属的类别输出分割结果包括：用图像中同一类别中的所有像素的强度的均值更新该同一类别中的所有像素的强度，得到均值图像；通过显示装置显示该均值图像。

结合第一方面，本发明实施例提供了第一方面的第四种可能的实施方式，还包括：设置常数，上述常数包括：类数k、Gamma分布形态参数α、归一化常数A和邻域作用系数η。

第二方面，本发明实施例还提供了一种多视SAR图像分割装置，包括：读取模块，用于读取待分割多视SAR图像；上述待分割多视SAR图像使用特征场z表示：

z＝{z_i(x_i,y_i):i＝1,…,n}

其中，i为像素索引，n为总像素数，z_i为像素i的强度，(x_i，y_i)∈D为像素i的格点位置，D为图像域；

双权重初始化模块，用于初始化双权重w：

w_i＝(w_il:l＝1,…,k)

其中，w_il为类属l包含像素i的类属权重，满足k为类数；

权重参数迭代更新模块，用于重复执行下述计算：计算Gamma分布尺度参数β：

其中，α为Gamma分布形态参数；

计算图像的或然率p(z|w)；

其中p(z_i|w_il)为以双权重Gamma混合模型定义的z_i的概率密度，如下：

计算上述双权重w的分布函数p(w)：

其中，A为归一化系数，η为邻域作用系数，N_i为以(x_i,y_i)为中心的8邻域像素集合，且满足i’∈N_i,i’≠i；

计算品质函数L，上述品质函数L为或然率p(z|w)与分布函数p(w)联合概率分布函数的对数函数；

根据梯度法更新上述双权重w：

w^(t+1)＝w^(t)+ξΔw^(t)

其中，t为迭代次数，ξ为步长，Δw^(t)为梯度，如下：

将更新后的双权重w代入品质函数L；

直到|L^(t+1)-L^(t)|小于预设的阈值ε时，停止执行上述计算；

类别确定模块，用于根据当前的双权重w确定图像中各个像素所属的类别；

输出模块，用于按照各个像素所属的类别对图像进行分割，输出分割结果。

结合第二方面，本发明实施例提供了第二方面的第一种可能的实施方式，其中，上述权重参数迭代更新模块包括：第一计算单元，用于计算或然率p(z|w)与分布函数p(w)的联合概率分布函数p(z,w)：

p(z,w)＝p(z|w)p(w)；

第二计算单元，用于对上述联合概率分布函数取对数：

L(w)＝log p(z,w)＝log p(z|w)+log p(w)。

结合第二方面，本发明实施例提供了第二方面的第二种可能的实施方式，其中，类别确定模块包括：最大值计算单元，用于计算双权重的最大值：

c_i＝arg max{w_il,l＝1,...,k}；

其中，c_i为第i个像素所属类别的标号；

类别确定单元，用于将上述最大值作为像素所属的类别。

结合第二方面，本发明实施例提供了第二方面的第三种可能的实施方式，其中，上述输出模块包括：均值单元，用于图像中同一类别中的所有像素的强度的均值更新该同一类别中的所有像素的强度，得到均值图像；显示单元，用于通过显示装置显示该均值图像。

结合第二方面，本发明实施例提供了第二方面的第四种可能的实施方式，还包括：常数设置模块，用于设置常数，上述常数包括：类数k、Gamma分布形态参数α、归一化常数A和邻域作用系数η。

本发明实施例提供的多视SAR图像分割方法及装置，采用基于类属场表示的双权重Gamma混合模型对图像特征场进行建模；为了引入邻域关系，基于马尔科夫随机场模型(Markov Random Filed,MRF)，结合误差平方和理论将中心像素类属权重视为均值，构造其与邻域像素类属权重差的平方和函数以描述邻域窗口内像素的差异性，用误差平方和表征差异能充分利用类属信息，增强了图像的分割精度。同时在估计Gamma混合模型中的双权重时，本发明实施例采用梯度法对双权重进行求解，迭代速度快、不易出现局部最优问题。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，应当理解，以下附图仅示出了本发明的某些实施例，因此不应被看作是对范围的限定，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他相关的附图。

图1示出了本发明实施例所提供的一种多视SAR图像分割装置的流程图；

图2示出了本发明实施例所提供的一种多视SAR图像分割装置的结构示意图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。因此，以下对附图中提供的本发明的实施例的详细描述并非旨在限制要求保护的本发明的范围，而是仅仅表示本发明的选定实施例。基于本发明的实施例，本领域技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

考虑到现有技术中的多视SAR图像分割方法存在的抗噪性能及准确性差的问题，本发明实施例提供了一种多视SAR图像分割方法及装置，该技术可以采用相应的软件和硬件实现。下面通过实施例进行描述。

实施例1

图1示出了本发明实施例提供的多视SAR图像分割方法的流程示意图。下面将对图1所示方法的具体流程进行详细阐述。

步骤S110，读取待分割多视SAR图像。

该待分割多视SAR图像使用特征场z表示：

z＝{z_i(x_i,y_i):i＝1,…,n}

其中，i为像素索引，n为总像素数，z_i为像素i的强度，(x_i，y_i)∈D为像素i的格点位置，D为图像域。

步骤S120，初始化双权重w。

在本实施例中，定义双权重w表示像素与类属之间关系，具体如下：

w_i＝(w_il:l＝1,…,k)

其中，w_il为类属l包含像素i的类属权重，满足k为类数。其中，初始双权重随机生成，行表示像素，列表示类别，矩阵中元素取值范围为0～1，并且满足每一行和为1。

步骤S130，计算Gamma分布尺度参数β。

在本实施例中，定义Gamma分布尺度参数β_l为关于类属权重w_il的函数，具体如下：

其中，α为Gamma分布形态参数，对于多视SAR图像，α等于其视数。基于双权重，结合Gamma混合模型定义特征场概率分布，其中，混合模型中分布尺度参数与分布形态参数的乘积为均值。

步骤S140，计算上述图像的或然率p(z|w)。

假设图像中各像素特征值相互独立，定义上述或然率为：

其中p(z_i|w_il)为以双权重Gamma混合模型定义的z_i的概率密度，如下：

步骤S150，计算上述双权重w的分布函数p(w)。

具体地，以误差平方和理论为基础，结合MRF定义中心像素与其邻域像素差异性，进而得到双权重w的分布函数p(w)为：

其中，A为归一化系数，控制聚类尺度；η为邻域作用系数，表征邻域影响强度；N_i为以(x_i,y_i)为中心的8邻域像素集合，且满足i’∈N_i,i’≠i。

步骤S160，计算品质函数L，该品质函数L为或然率p(z|w)与分布函数p(w)联合概率分布函数的对数函数。

其中，计算品质函数L具体包括：

计算或然率p(z|w)与分布函数p(w)的联合概率分布函数p(z,w)：

p(z,w)＝p(z|w)p(w)；

品质函数L定义为p(z,w)的对数函数，即对上式取对数：

L(w)＝log p(z,w)＝log p(z|w)+log p(w)。

步骤S170，根据梯度法更新上述双权重w。

w^(t+1)＝w^(t)+ξΔw^(t)

其中，t为迭代次数，ξ为步长，Δw^(t)为梯度，该梯度定义为上述品质函数L的导数如下：

步骤S180，将更新后的双权重w代入品质函数L。

步骤S190，重复执行上述步骤S130-S180，直到|L^(t+1)-L^(t)|小于预设的阈值ε时，停止执行上述步骤，根据当前的双权重w确定图像中各个像素所属的类别。

确定上述各个像素所属的类别具体包括以下步骤：

(1)计算上述双权重的最大值：

c_i＝arg max{w_il,l＝1,...,k}；

其中，c_i为第i个像素所属类别的标号；

(2)将该最大值作为该像素所属的类别。

由上述步骤S120中双权重的定义可知，其是像素与类属的关系矩阵，每一行中最大值对应的列即为像素所属类别。

步骤S200，按照各个像素所属的类别输出分割结果。

具体地，通过上述步骤确定待分割图像的每个像素属于哪个类别后，将同一类中所有像素的强度(即像素的灰度值)求均值，将该均值作为此类的强度值即得到分割结果。用上述图像中同一类别中的所有像素的强度的均值更新该同一类别中的所有像素的强度，得到均值图像；通过显示装置显示该均值图像。

在本实施例的方法实际实施时，在输入待分割的图像后执行上述步骤前，还包括设置常数的步骤，该常数具体包括：类数k、Gamma分布形态参数α、归一化常数A和邻域作用系数η。

本发明实施例提供的多视SAR图像分割方法，通过采用基于类属场表示的双权重Gamma混合模型对图像特征场进行建模；为了引入邻域关系，基于马尔科夫随机场模型(Markov Random Filed,MRF)，结合误差平方和理论将中心像素类属权重视为均值，构造其与邻域像素类属权重差的平方和函数以描述邻域窗口内像素的差异性，用误差平方和表征差异能充分利用类属信息，增强了图像的分割精度；同时在估计Gamma混合模型中的双权重时，本发明实施例采用梯度法对双权重进行求解。具体如下：

(1)利用Gamma混合模型刻画SAR图像特征场的概率分布，以类属场表征像素与类属间关系，该关系用矩阵表示，并将该类属矩阵作为Gamma混合模型的双权重，从像素与聚类两个角度共同作用于Gamma分布，较传统混合模型中以向量表示的仅与聚类相关的单权重能够更好地刻画特征场；

(2)基于误差平方和及马尔科夫随机场(Markov Random Field,MRF)理论，以中心像素与其8邻域像素类属权重差的平方和刻画邻域窗口内特征差异程度，进而定义双权重的分布函数，当平方和越大时，说明窗口内像素特征差异越大，即其属于不同类别的概率越大；

(3)利用贝叶斯定理定义特征场和双权重的联合概率分布函数，以其对数函数作为品质函数，并利用梯度法求解双权重，以获得最小化品质函数对应的最佳估计值。

本发明实施例的效果可通过如下仿真实验进一步说明：

(1)仿真实验条件

本实施例在CPU为Core(TM)i5-3470 3.20GHz的Windows 7旗舰版系统上使用MATLAB 2011a软件编程实现仿真。

仿真数据1为模拟SAR图像，包含3个同质区域，其由标准模板图像添加形态参数为4，尺度参数分别为2、10、20的Gamma分布随机数得到。仿真数据2为真实SAR影像，包含2个同质区域，由于真实遥感图像无标准模板，因此将手绘模板视为标准模板。以上仿真数据图像大小均为128×128，图像总像素数n＝16384。

(2)仿真实验结果

为了证明本实施例算法的有效性，分别以模板图像为标准，对本实施例算法及对比算法分割结果生成混淆矩阵，并计算其用户精度、产品精度及Kappa值，以对本实施例方法进行定量分析(如表1所示)，其中“-”表示图像中无此区域。由表1可知，本实施例的总精度及Kappa值均高于对比算法，从数字角度精确地验证了本实施例方法的有效性。本实施例的方法抗噪性能较好，分割结果误分现象较少，分割边界拟合精确；而对比算法不能有效克服复杂噪声，分割结果中误分像素较多，导致视觉效果极差。

表1

实施例2

结合前述实施例，本实施例提供了一种多视SAR图像分割装置，参见图2所示的多视SAR图像分割装置的结构示意图，该装置包括：读取模块301、双权重初始化模块302、权重参数迭代更新模块303、类别确定模块304和显示模块305。

具体介绍如下：

读取模块301，用于读取待分割多视SAR图像。

该待分割多视SAR图像使用特征场z表示：

z＝{z_i(x_i,y_i):i＝1,…,n}

其中，i为像素索引，n为总像素数，z_i为像素i的强度，(x_i，y_i)∈D为像素i的格点位置，D为图像域。

双权重初始化模块302，用于初始化双权重w。

在本实施例中，定义双权重w表示像素与类属之间关系，具体如下：

w_i＝(w_il:l＝1,…,k)

其中，w_il为类属l包含像素i的类属权重，满足k为类数。其中，初始双权重随机生成，行表示像素，列表示类别，矩阵中元素取值范围为0～1，并且满足每一行和为1。

权重参数迭代更新模块303，用于重复执行下述计算：

计算Gamma分布尺度参数β。

定义Gamma分布尺度参数β_l为关于类属权重w_il的函数，具体如下：

其中，α为Gamma分布形态参数，对于多视SAR图像，α等于其视数。基于双权重，结合Gamma混合模型定义特征场概率分布，其中，混合模型中分布尺度参数与分布形态参数的乘积为均值。

计算图像的或然率p(z|w)。

假设图像中各像素特征值相互独立，定义上述或然率为：

其中p(z_i|w_il)为以双权重Gamma混合模型定义的z_i的概率密度，如下：

计算上述双权重w的分布函数p(w)。

具体地，以误差平方和理论为基础，结合MRF定义中心像素与其邻域像素差异性，进而得到双权重w的分布函数p(w)为：

其中，A为归一化系数，控制聚类尺度；η为邻域作用系数，表征邻域影响强度；N_i为以(x_i,y_i)为中心的8邻域像素集合，且满足i’∈N_i,i’≠i。

计算品质函数L，上述品质函数L为或然率p(z|w)与分布函数p(w)联合概率分布函数的对数函数。

上述权重参数迭代更新模块303包括：第一计算单元，用于计算或然率p(z|w)与分布函数p(w)的联合概率分布函数p(z,w)：

p(z,w)＝p(z|w)p(w)；

第二计算单元，用于对上式取对数：

L(w)＝log p(z,w)＝log p(z|w)+log p(w)。

根据梯度法更新上述双权重w。

w^(t+1)＝w^(t)+ξΔw^(t)

其中，t为迭代次数，ξ为步长，Δw^(t)为梯度，该梯度定义为上述品质函数L的导数如下：

将更新后的双权重w代入品质函数L；

直到|L^(t+1)-L^(t)|小于预设的阈值ε时，停止执行上述步骤。

类别确定模块304，用于根据当前的双权重w确定图像中各个像素所属的类别。

上述类别确定模块304包括：最大值计算单元，用于计算双权重的最大值：

c_i＝arg max{w_il,l＝1,...,k}

其中，c_i为第i个像素所属类别的标号；

类别确定单元，用于将上述最大值作为像素所属的类别。

由上述双权重的定义可知，其是像素与类属的关系矩阵，每一行中最大值对应的列即为像素所属类别。

输出模块305，用于按照各个像素所属的类别输出分割结果。

具体地，上述输出模块305包括：

均值单元，用于图像中同一类别中的所有像素的强度的均值更新该同一类别中的所有像素的强度，得到均值图像；

显示单元，用于通过显示装置显示该均值图像。

在本实施例的装置实际实施时，还包括：常数设置模块，用于设置常数，上述常数包括：类数k、Gamma分布形态参数α、归一化常数A和邻域作用系数η。

本实施例所提供的多视SAR图像分割装置的实现原理及产生的技术效果和前述实施例相同，为简要描述，本实施例部分未提及之处，可参考前述实施例中相应内容。

上述功能如果以软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以权利要求的保护范围为准。