一种基于区域划分与Gamma混合模型的多视SAR图像分割方法与流程

本发明属于图像处理技术领域，具体涉及一种基于区域划分与Gamma混合模型的多视SAR图像分割方法。

背景技术：

图像分割是指将图像分成各具同质特性区域的过程，是图像处理的关键环节之一，涉及到工业检测、环境监测和军事等多个领域。合成孔径雷达(Synthetic Aperture Radar，SAR)由于具有覆盖范围广、全天时、全天候等特点，被看作新一代重要遥感信息源。但由于SAR的相干成像特性，使得SAR图像中存在大量斑点噪声，给SAR图像分割造成很大困难。

混合模型是一种统计模型，它属于一种半参数的密度估计方法，不局限于特定的概率密度函数形式，且模型的复杂度仅与所求解相关，与样本大小无关，所以在许多领域有着广泛的应用。高斯混合模型由于其形式简单、计算方便，成为一种比较普遍应用的混合模型。到目前为止，有学者利用高斯混合模型对实际SAR图像的背景杂波进行拟合，得到分布模型后进而完成目标检测；或结合高斯混合模型与马尔科夫随机场(Markov random field，MRF)分割灰度图像。虽然高斯混合模型的应用非常广泛，但由于SAR图像中大量斑点噪声的干扰，使得实际数据呈现出非高斯特性，导致高斯混合模型不能非常有效地对SAR图像进行分割。

根据对SAR图像噪声的统计分析，可以证明在分布式地物覆盖区域内SAR图像中像素强度满足Gamma分布；由此，提出一种将MRF与Gamma混合模型相结合的SAR图像分割方法；此外，Gamma混合模型亦被用于分割覆盖海洋溢油的SAR图像，结果表明该模型能够成功分割海洋溢油和海洋背景这两类目标；但将Gamma混合模型用于SAR图像分割时，多数学者只利用图像像素的灰度信息，因此不能很好地消除SAR图像中斑点噪声的影响。

技术实现要素：

针对现有技术的不足，本发明提出一种基于区域划分与Gamma混合模型的多视SAR图像分割方法，该方法首先将图像域区域化，以区域子块为处理单元，可以有效解决噪声问题，并结合Gamma混合模型建模各划分区域像素强度的分布特性，由此使得分割结果最优化。

一种基于区域划分与Gamma混合模型的多视SAR图像分割方法，包括以下步骤：

步骤1、获得待分割原始图像的图像域上的随机特征场；

步骤2、在图像域上随机产生多个生成点，根据生成点采用Voronoi划分方法将图像域划分成一系列区域子块；

步骤3、建立图像分割模型，具体步骤如下：

步骤3-1、构建每个图像域上每个像素点属于任一像素灰度类别的概率密度函数，即构建Gamma混合模型；

步骤3-2、根据步骤3-1构建的概率密度函数，获得图像的联合概率密度函数；

步骤3-3、通过构建Potts模型的方式描述图像整体空间依赖关系，进而获得某一区域子块内像素取得某类标号时，该多边形与其邻域间的依赖关系，即获得每一个区域子块的先验概率；

步骤3-4、构建以区域子块中所有像素的分布总和最大化为目的的目标函数，该目标函数的待估计量为像素点标号场；

步骤4、采用EM算法，在标号场为均质场的前提下，通过反复迭代更新生成点的方式，更新标号场的对数似然函数，通过比较获得标号场的最大值，获得图像最优分割结果。

步骤3-1所述的构建每个图像域上每个像素点属于任一像素灰度类别的概率密度函数，即构建Gamma混合模型，具体如下：

其中，p(z_i|(x_i，y_i)∈P_j，π_j，β_j)表示区域子块P_j中像素点z_i的概率密度函数，π_j＝(π_j1，...，π_jk)表示区域子块P_j内Gamma混合模型的权重系数矢量，π_jl表示区域子块P_j内第l类Gamma混合模型的权重系数，l＝1，...，k，k表示像素点标号总类数，j表示第j个区域子块；β_j＝(β_j1，...，β_jk)，β_jl表示区域子块P_j内第l类的Gamma混合模型的尺度参数；α表示Gamma混合模型分布的形状参数，即为多视SAR图像的视数，Γ(·)为Gamma函数，Ga(z_i|β_jl)表示区域子块P_j中的像素点z_i属于第l类的Gamma分布。

步骤3-2所述的图像的联合概率密度函数，具体如下：

其中，p(z|G，π，β)表示图像z的联合概率密度函数，G表示生成点集，m表示区域子块总个数，p(z_i|(x_i，y_i)∈P_j，π_j，β_j)表示区域子块P_j中像素点z_i的概率密度函数，π_j＝(π_j1，...，π_jk)表示区域子块P_j内Gamma混合模型的权重系数矢量，π_jl表示区域子块P_j内第l类Gamma混合模型的权重系数，l＝1，...，k，k表示像素点标号总类数，j表示第j个区域子块；β＝{β_j：j＝1，...，m}，β_j＝(β_j1，...，β_jk)，β_jl表示区域子块P_j内第l类的Gamma混合模型的尺度参数；Ga(z_i|β_jl)表示区域子块P_j中的像素点z_i属于第l类的Gamma分布。

步骤3-3所述的通过构建Potts模型的方式描述图像整体空间依赖关系，具体如下：

其中，p(L)表示随机标号场L的联合概率，j表示第j个区域子块，m表示区域子块总个数，p(L_j|L_j′，j′∈NP_j)表示区域子块的先验概率，NP_j表示区域子块P_j的邻域多边形集合，l＝1，...，k，k表示像素点标号总类数，η表示区域子块之间的空间作用参数；t(a，b)为指示函数，当a＝b时，t(a，b)＝0，否则，t(a，b)＝1；

所述的进而获得某一区域子块内像素取得某类标号时，该多边形与其邻域间的依赖关系，即获得每一个区域子块的先验概率，具体如下：

其中，π_jl表示当区域子块所属类别确定时其先验概率，即第j个区域子块内像素点取得第l类标号时，该区域子块与其邻域间区域子块的依赖关系。

步骤3-4所述的构建以区域子块中所有像素的分布总和最大化为目的的目标函数，具体如下：

其中，表示标号场L的估计值，l＝1，...，k，k表示像素点标号总类数，π_jl表示当区域子块所属类别确定时其先验概率，即第j个区域子块内像素点取得第l类标号时，该区域子块与其邻域间区域子块的依赖关系，Ga(z_i|β_jl)表示区域子块P_j中的像素点z_i属于第l类的Gamma分布。

步骤4所述的采用EM算法，在标号场为均质场的前提下，通过反复迭代更新生成点的方式，更新标号场的对数似然函数，通过比较获得标号场的最大值，获得图像最优分割结果，具体步骤如下：

步骤4-1、构建图像所对应标号场的对数似然函数；

步骤4-2、使标号场的对数似然函数的一阶导数等于0，估计对数似然函数中的尺度参数集；

步骤4-3、对生成点进行位移操作，获得对数似然函数中新的生成点集，返回执行步骤4-1获得更新后的标号场的对数似然函数；

步骤4-4、判断更新后的标号场的对数似然函数是否大于更新前的标号场的对数似然函数，若是，则保留当前生成点集，并返回执行步骤4-3，否则，保留之前生成点集，并返回执行步骤4-3；

步骤4-5、当达到最大迭代次数时，获得最优生成点集，即获得最大标号场的对数似然函数，获得图像最优分割结果。

步骤4-1所述的构建图像所对应标号场的对数似然函数，具体如下

其中，L(β，G)表示图像所对应标号场的对数似然函数，G表示生成点集，尺度参数集β＝{β_j：j＝1，...，m}，β_j＝(β_j1，...，β_jk)，β_jl表示区域子块P_j内第l类的Gamma混合模型的尺度参数；j表示第j个区域子块，m表示区域子块总个数，π_j＝(π_j1，...，π_jk)表示区域子块P_j内Gamma混合模型的权重系数矢量，π_jl表示区域子块P_j内第l类Gamma混合模型的权重系数，l＝1，...，k，k表示像素点标号总类数；Ga(z_i|β_jl)表示区域子块P_j中的像素点z_i属于第l类的Gamma分布。

本发明优点：

本发明提出一种基于区域划分与Gamma混合模型的多视SAR图像分割方法，该方法通过Gamma混合模型来模拟多视SAR图像的分布，从而达到分割SAR图像的目的，并使用Voronoi多边形来代替像素作为处理单元，能够有效消除SAR图像本身存在的斑点噪声影响，而且在Gamma混合模型的基础上利用Potts模型引入了邻域关系，同时约束了多边形的类属性，进而提高SAR图像的分割精度。

附图说明

图1为本发明一种实施例的基于区域划分与Gamma混合模型的多视SAR图像分割方法流程图；

图2为本发明一种实施例的利用Voronoi划分方法对图像域进行划分的示意图，其中，空心圆圈为生成点，实心圆圈为像素格点；

图3为本发明一种实施例的位移生成点的示意图；

图4为本发明一种实施例的待分割原始SAR图像；

图5为本发明一种实施例的分割结果；

图6为本发明一种实施例的对应图5的定性评价结果；

图7为本发明一种实施例的分割结果轮廓线与区域划分结果的叠加。

具体实施方式

下面结合附图对本发明一种实施例做进一步说明。

本发明实施例中，基于区域划分与Gamma混合模型的多视SAR图像分割方法，方法流程图如图1所示，包括以下步骤：

步骤1、获得待分割原始图像的图像域上的随机特征场；

本发明实施例中，输入待分割原始图像，z＝{z_i(x_i，y_i)：i＝1，...，n}，其中，i为像素索引，n为总像素数，z_i为像素i的强度值，(x_i，y_i)∈D是像素i在图像域D上的像素位置坐标；z可看作定义在图像域D上的随机特征场(也称特征场)Z＝{Z_i(x_i，y_i)：i＝1，...，n}的实现，Z_i(x_i，y_i)表征定义在像素位置(x_i，y_i)上的随机变量，而z_i可看作Z_i的实现；

本发明实施例中，以大小为128×128的256级灰度图像图4中图(a)为例，其总像素数n＝16384，z＝{40，37，...，148}；

步骤2、在图像域上随机产生多个生成点，根据生成点采用Voronoi划分方法将图像域划分成一系列区域子块；

步骤2-1、在图像域D上随机产生m个生成点，生成点集为G＝{(u_j，v_j)：j＝1，...，m}，其中，j为生成点索引，m为生成点个数，(u_j，v_j)∈D为生成点j的位置坐标；

本发明实施例中，设置m＝280，G＝{(1，8)，(1，10)，...，(128，77)}；

步骤2-2、如图2所示，利用Voronoi划分方法将图像域划分成一系列区域子块，可以表示为P＝{P_j：j＝1，...，m}，其中，P_j表示第j个Voronoi多边形(区域子块)；

本发明实施例中，Voronoi划分方法即通过计算各像素点与生成点之间的欧式距离，将该像素点划分到与之距离最近的生成点所在多边形中，第j个多边形的划分方法如下：

其中，δ(·，·)计算两点之间的欧式距离；

本发明实例中，共生成280个多边形，其中第一个多边形P₁包含14个像素点，P₁＝{(1，5)，(1，6)，(1，7)，(1，8)，(1，9)，(2，5)，(2，6)，(2，7)，(2，8)，(2，9)，(3，6)，(3，7)，(3，8)，(4，6)}；

步骤3、建立图像分割模型，具体步骤如下：

步骤3-1、构建每个图像域上每个像素点属于任一像素灰度类别的概率密度函数，即构建Gamma混合模型；

本发明实施例中，已知像素i的强度值z_i，当图像有k类时(图像区域间灰度差值大于阈值50，则认为区域所属不同类别)，则定义在(x_i，y_i)∈P_j上的概率密度函数即可被认为是Gamma混合模型；本发明实施例中，对3幅真实的多视SAR图像(见图4中图(a)、图(b)和图(c)所示，)进行分割实验，其中，图4中图(a)和图(b)的大小均为128×128像素，图4中图(c)大小为256×256像素，依据其灰度分布依次分为2、3和2类；

具体如下：

其中，p(z_i|(x_i，y_i)∈P_j，π_j，β_j)表示区域子块P_j中像素点z_i的概率密度函数，π_j＝(π_j1，...，π_jk)表示区域子块P_j内Gamma混合模型的权重系数矢量，π_jl表示区域子块P_j内第l类Gamma混合模型的权重系数，l＝1，...，k，k表示像素点标号总类数，j表示第j个区域子块；β_j＝(β_j1，...，β_jk)，β_jl表示区域子块P_j内第l类的Gamma混合模型的尺度参数；α表示Gamma混合模型分布的形状参数，即为多视SAR图像的视数，本发明实施例中，视数取值为4；Γ(·)为Gamma函数，Ga(z_i|β_jl)表示区域子块P_j中的像素点z_i属于第l类的Gamma分布；

步骤3-2、假设所有像素强度的概率分布均为独立的，根据步骤3-1构建的概率密度函数，获得图像的联合概率密度函数；具体如下：

其中，p(z|G，π，β)表示图像z的联合概率密度函数，G表示生成点集，m表示区域子块总个数；

步骤3-3、通过构建Potts模型的方式描述图像整体空间依赖关系，进而获得某一区域子块内像素取得某类标号时，该多边形与其邻域间的依赖关系，即获得每一个区域子块的先验概率；

本发明实施例中，假设同一多边形中的所有像素均属于同一类别，如多边形P_j内的所有像素均为标号L_j，L_j为隐含随机变量，在{1，...，k}中任意取值；于是所有L_j(j＝1，...，m)能够在划分区域上建模为随机标号场L，L＝{L_j：j＝1，...，m}，L的实现即对应图像分割结果；为了使邻域多边形具有空间依赖关系，上述标号场被假定为MRF，并利用Potts模型来刻画该MRF，p(L)即为利用Potts模型来描述的图像整体具有的空间依赖关系：

具体如下：

其中，p(L)表示随机标号场L的联合概率，p(L_j|L_j′，j′∈NP_j)表示区域子块的先验概率，NP_j表示区域子块P_j的邻域多边形集合，η表示区域子块之间的空间作用参数；t(a，b)为指示函数，当a＝b时，t(a，b)＝0，否则，t(a，b)＝1；

获得第j个区域子块内像素取得第1类标号时，该多边形与其邻域间的依赖关系，即获得每一个区域子块的先验概率，具体如下：

其中，π_jl表示当区域子块所属类别确定时其先验概率，即第j个区域子块内像素点取得第l类标号时，该区域子块与其邻域间区域子块的依赖关系；

步骤3-4、构建以区域子块中所有像素的分布总和最大化为目的的目标函数，该目标函数的待估计量为像素点标号场；

本发明实施例中，当在一个Voronoi多边形中的所有像素的分布总和达到最大化时，可获得标号场L的估计值，以多边形P_j为例：

其中，表示标号场L的估计值；

步骤4、采用EM算法，在标号场为均质场的前提下，通过反复迭代更新生成点的方式，更新标号场的对数似然函数，通过比较获得标号场的最大值，获得图像最优分割结果；

本发明实施例中，为估计已知图像z条件下的分割模型参数，通过EM算法分别来改变参数β和G，以达到最大化L的目的，具体步骤如下：

将每次EM迭代更新第j个多边形P_j的先验概率π_jl和标号L_j时，均将L看作是均质场，且π_jl和L_j是通过公式(4)和公式(5)依次得到的；

步骤4-1、构建图像所对应标号场的对数似然函数；具体如下

其中，L(β，G)表示图像所对应标号场的对数似然函数；

步骤4-2、使标号场的对数似然函数的一阶导数等于0，估计对数似然函数中的尺度参数集；

本发明实施例中，令公式(6)中的一阶导数等于0，以此估计β：

令

则

步骤4-3、对生成点进行位移操作，获得对数似然函数中新的生成点集，返回执行步骤4-1获得更新后的标号场的对数似然函数；

本发明实施例中，如图3所示，位移生成点操作：若当前Voronoi划分后的生成点集为G^(r)＝{(u₁^(r)，v₁^(r))，...，(u_j^(r)，v_j^(r))，...，(u_m^(r)，v_m^(r))}，G^(r)表示第r次迭代的生成点集合，以1/m的概率随机选取一生成点(u_j^(r)，v_j^(r))∈P_j，移动该生成点到(u_j^*，v_j^*)∈P_j，则G^*＝{(u₁^(r)，v₁^(r))，...，(u_j^*，v_j^*)，...，(u_m^(r)，v_m^(r))}，相应的Voronoi划分也同样变为P^*＝{P_j^*：j＝1，...，m}；

步骤4-4、判断更新后的标号场的对数似然函数L(β^(r)，G^*)是否大于更新前的标号场的对数似然函数L(β^(r)，G^(r))，若是，则操作被接受，即G(^r+1)＝G^*且P^(r+1)＝P^*，并返回执行步骤4-3，否则，操作被拒绝，即G^(r+1)＝G^(r)且P^(r+1)＝P^(r)，并返回执行步骤4-3；

本发明实施例中，此次迭代完成之后，通过改变生成点集G(亦为图像划分结果)，能得到优化的目标函数，由此进一步确定参数β(执行步骤4-2)，以进行下一次迭代(返回执行步骤4-3)。

步骤4-5、当达到最大迭代次数时，获得最优生成点集，即获得最大标号场的对数似然函数，获得图像最优分割结果。

图5为对应的分割结果，通过图5中图(a)、图(b)和图(c)可以看出，本发明能够正确地将图像类别分割出来，并有效消除了SAR图像中噪声的影响，且对于大小不同的多视SAR图像都能够得到较好的分割结果。

为进一步证明本发明的可行性与有效性，通过定性评价，即将分割结果轮廓线叠加在原始图像上来直观判断分割结果是否准确，定性评价结果如图6所示。通过图6中图(a)、图(b)和图(c)可以看出，本发明能够将图像较准确地分割出来。

本发明实施例在分割过程中均将多边形数目设置为m＝280，从图7中图(a)、图(b)和图(c)能够看出，图7背景为区域划分结果，然后将分割结果轮廓线与其进行叠加，能够看出哪些多边形拟合了图像中的同一类别；由此可以看出本发明能够正确分割多视SAR图像。

以上所述具体实例，对本发明的目的、技术方案和有益效果进行了进一步详细说明。所应理解的是，以上所述仅为本发明的具体实例而已，并不用于限制本发明，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。