结合K‑S距离和RJMCMC算法的SAR图像分割方法与流程

本发明涉及图像处理技术领域，尤其涉及一种结合k-s(kolmogorov-smirnov)距离和rjmcmc(reversiblejumpmarkovchainmontecarlo，可逆跳马尔科夫链蒙特卡洛)算法的sar图像分割方法。

背景技术：

图像分割是sar(syntheticapertureradar，合成孔径雷达)图像处理的主要任务。目前，已经提出很多相关算法，主要分为：阈值分割、聚类分割和统计分割等，其中，统计分割被认为最为有效的方法。

目前，大部分的统计分割均是以建立图像光谱测度的统计模型为基础，即假设图像光谱测度服从某一统计分布，如假设sar图像强度服从gamma分布、gaussian分布等，建立其强度的统计模型。但这些统计模型仅仅是对sar图像强度特征的一种假设，有时难以反映其本质规律，特别是对于sar图像分割问题；由于sar图像地物目标复杂、多样且斑点噪声明显，导致难以准确建立图像光谱测度的统计模型。

技术实现要素：

针对现有技术不足，本发明提供一种结合k-s距离和rjmcmc算法的sar图像分割方法。

本发明的技术方案如下：

一种结合k-s距离和rjmcmc算法的sar图像分割方法，包括如下步骤：

步骤1：输入待分割图像，并将其定义为图像域上的特征场的一个实现；

步骤2：利用规则划分技术将输入的sar图像的图像域划分为多个规则子块；

步骤3：在划分的图像域上，建立基于规则子块的sar图像分割模型，包括以下步骤：

步骤3.1：在划分的图像域上，建立特征场与标号场的关系模型；

步骤3.2：在划分的图像域上，建立标号场模型；

步骤3.3：在步骤3.1和步骤3.2的基础上，建立基于规则子块的sar图像分割模型；

步骤4：针对已建立的sar图像分割模型，设定迭代次数并设计合理的移动操作，每次迭代中，遍历所有移动操作，利用rjmcmc算法，求解基于规则子块的sar图像分割模型；具体包括以下步骤：

步骤4.1：通过更新标号场中的标号来更新标号场；

步骤4.2：在步骤4.1更新标号场的基础上通过增加或减少规则子块个数来更新划分的图像域及标号场，更新图像域过程中随机选择增加规则子块或减少规则子块的操作；

步骤4.3：将步骤4.2中更新后的图像域上的标号场的实现代入非约束gibbs概率分布函数，得到非约束gibbs概率分布函数的一个函数值；

步骤4.4：按设定的迭代次数重复执行步骤4.1至步骤4.3，得到非约束gibbs概率分布函数的一个函数值集合，集合中函数值最大时所对应的分割结果为基于规则子块的sar图像分割模型的最优解；

步骤5：输出sar图像的分割结果。

优选地，将步骤1中输入的待分割sar图像x＝{xs,s＝1,...,s}定义为图像域p上的特征场x＝{xs,s＝1,...,s}的一个实现，其中，s为像素索引，xs为像素s的强度，s为图像像素总数，xs表示像素s的强度的随机变量；

步骤2利用规则划分技术划分图像域，具体为，利用规则划分技术将图像域p划分成j个规则子块，即，p＝{pj,j＝1,...,j}，其中，pj表示规则子块，j表示规则子块索引，j表示规则子块的总个数，规则子块pj的行数或列数为2的整数倍数，允许最小的规则子块包括2×2个像素点；sar图像的每个同质区域由一个或多个具有相同标号的规则子块拟合而成。

优选地，步骤3中首先在规则划分的图像域上，定义标号场l＝{lj,j＝1,...,j}，其中，lj表示规则子块pj所属标号的随机变量，lj∈{1,...,k}，k为图像的总类别数；在规则划分的图像域中，特征场x定义为x＝{xj,j＝1,...,j}，其中xj＝{xs；s∈pj}为规则子块pj内所有像素为s的强度的随机变量的集合；标号场l的每一个实现l＝{lj,j＝1,...,j}为sar图像x对应的分割结果，lj为规则子块pj的标号，且规则子块pj内所有像素的标号均为lj；

建立基于规则子块的sar图像分割模型的具体方法为：

步骤3.1：在规则划分的图像域上，以同质区域间异质性势能函数之和作为特征场与标号场的关系模型，建立的特征场与标号场的关系模型ux(x,l)表示为：

其中，为异质性势能函数，xj＝{xs；s∈pj}表示规则子块pj内所有像素为s的强度的集合，表示图像域中像素为s、标号为lj的所有像素的强度的集合，其中，ls为图像域中像素s的标号；

在未知同质区域内像素光谱测度的统计分布模型中，k-s距离能以函数间最大垂直距离作为两种不同统计分布间的相似性测度。因此，以k-s距离为图像分割准则，用k-s距离定义异质性势能函数其表示为：

其中，dks代表k-s距离，即直方图和的最大间距，和分别代表两个数据集合xj＝{xs；s∈pj}和的采样分布函数，和分别表示为：

其中，n1和n2分别为两个数据集合xj和中元素的个数，h为强度值索引，若对于n位图像，h∈[0,2ⁿ-1]；

步骤3.2：在划分的图像域上，利用势能函数定义标号场模型，标号场模型ul(l)表示为：

其中，β为邻域子块的空间作用参数，npj为规则子块pj的八邻域规则子块的集合，j′为规则子块pj的邻域规则子块pj′的索引，lj′为规则子块pj′的标号；若lj＝lj′，则δ(lj,lj′)＝1；若lj≠lj′，则δ(lj,lj′)＝0；

步骤3.3：结合步骤3.1中特征场与标号场的关系模型和步骤3.2中的标号场模型，定义图像分割的全局势能函数u(x,l)，其表示为：

利用非约束gibbs概率分布函数刻画全局势能函数u(x,l)，得到基于规则子块的sar图像分割模型g(x,l)，其表示为：

优选地，步骤4.1中更新标号场中标号的操作方法为：从图像域p＝{pj,j＝1,...,j}中随机抽取一个规则子块pj，其对应标号为lj；从图像的总类别数集合{1,…,k}中随机抽取规则子块pj的候选标号且未被抽取的规则子块所属标号不变；利用非约束gibbs概率分布函数g(x,l)，得到更新标号场中标号lj为的接受率为：

其中，

其中，

从0～1中产生一个随机数，判断该随机数与接受率的大小，当所述接受率大于该随机数时，接受本次更新标号场中标号的操作，否则放弃本次更新标号场中标号的操作。

优选地，步骤4.2中增加规则子块个数是通过规则子块的分裂而实现的，具体实现步骤为：

步骤4.2.1：从通过步骤4.1更新标号场后的由j个规则子块拟合而成的图像域p中随机抽取一个规则子块pj，其对应标号为lj；

步骤4.2.2：判断所选规则子块pj能否实现分裂操作；如果pj的像素数大于4且其行数或列数为2的整数倍数，则进行分裂pj操作，执行步骤4.2.3，增加图像域中规则子块个数；否则不分裂规则子块pj，不增加图像域中规则子块个数，结束本次分裂操作；

步骤4.2.3：将图像域p中能实现分裂操作的规则子块pj分裂成两个新的规则子块pj1和pj2，新的规则子块pj1的对应标号为新的规则子块pj2对应标号为且图像域中的规则子块数增加一个，图像域p更新为图像域j^*表示更新后的规则子块的总个数；未分裂的规则子块不变，对应标号也不变；图像域p更新为图像域p^*的接受率为：

as(p,p^*)＝min{1,rs}

其中，

其中，为更新后的图像域p^*的图像的一个实现，为更新后的标号场的实现，为更新后的标号场的标号；

步骤4.2.4：从0～1中产生一个随机数，判断该随机数与接受率as(p,p^*)的大小，当所述接受率大于该随机数时，接受本次通过增加规则子块来更新图像域及标号场的操作，否则放弃本次通过增加规则子块来更新图像域及标号场的操作。

优选地，步骤4.2中减少规则子块个数通过图像域p中随机抽取的一个规则子块和其任一个邻域规则子块合并成一个新的规则子块来实现，该操作为分裂规则子块的对偶操作，因此，通过减少规则子块个数来更新图像域的接受率为：

am(p,p^*)＝min{1,1/rs}

从0～1中产生一个随机数，判断该随机数与接受率am(p,p^*)的大小，当所述接受率大于该随机数时，接受本次通过减少规则子块来更新图像域及标号场的操作，否则放弃本次通过减少规则子块来更新图像域及标号场的操作。

由上述技术方案可知，本发明的有益效果在于：本发明提供一种结合k-s距离和rjmcmc算法的sar图像分割方法，在未知sar图像模型的情况下，结合k-s距离和规则划分技术，建立的基于规则划分的图像分割模型，克服了sar图像固有斑点噪声对图像分割的影响，提高了采用本发明方法所分割的图像的同质区域的区域一致性；针对建立的基于规则划分的图像分割模型，在rjmcmc算法中设计更新标号场标号和增加或减少规则子块个数的移动操作，克服了图像分割时所划分的规则子块跨越同质区域边界而降低了分割边缘精度的问题，提高了分割结果的边缘准确性。

附图说明

图1为本发明实施例提供的结合k-s距离和rjmcmc算法的sar图像分割方法流程图；

图2为本发明实施例提供的三个待分割图像；其中，(a)为第一个待分割图像；(b)为第二个待分割图像；(c)为第三个待分割图像；

图3为本发明实施例提供的一种规则划分示意图；

图4为本发明实施例提供的更新标号场的操作流程图；

图5为本发明实施例提供的增加规则子块数的操作流程图；

图6为图2中三个待分割图像的分割结果图像。

具体实施方式

下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述。以下实施例用于说明本发明，但不用来限制本发明的范围。

一种结合k-s距离和rjmcmc算法的sar图像分割方法，如图1所示，包括的具体步骤如下。

步骤1：输入待分割图像，并将其定义为图像域上的特征场；

图2(a)、图2(b)和图2(c)为本实施例采用的3幅待分割的sar图像。作如下定义：sar图像x＝{xs,s＝1,...,s}定义为图像域p上的特征场x＝{xs,s＝1,...,s}的一个实现，其中，s为像素索引，xs为像素s的强度，s为图像像素总数，xs表示像素s的强度的随机变量。

步骤2：利用规则划分技术划分图像域；

利用规则划分技术将图像域p划分成j个规则子块，即，p＝{pj,j＝1,...,j}，其中，pj表示规则子块，j表示规则子块索引，j表示图像域p中规则子块总个数。每个同质区域由一个或多个规则子块拟合而成。图3为将图像域利用规则划分技术划分为10个长方形规则子块，即，p＝{pj,j＝1,...,10}，得到规则子块p1-p10，规则子块p1-p10由多个像素点构成，每个规则子块的行数或列数为2的整数倍数，允许最小的规则子块包括2×2个像素点。

步骤3：在划分的图像域上，建立基于规则子块的sar图像分割模型；

在规则划分的图像域上，定义标号场l＝{lj,j＝1,...,j}，其中，lj表示规则子块pj所属标号的随机变量，lj∈{1,...,k}，k为图像的总类别数；在规则划分的图像域中，特征场x定义为x＝{xj,j＝1,...,j}，其中xj＝{xs；s∈pj}为规则子块pj内所有像素的强度随机变量的集合；标号场l的每一个实现l＝{lj,j＝1,...,j}为sar图像x对应的分割结果，lj为规则子块pj的标号，且规则子块pj内所有像素的标号均为lj；

建立基于规则子块的sar图像分割模型的具体方法为：

步骤3.1：在划分的图像域上，建立特征场与标号场的关系模型；

以同质区域间异质性势能函数之和作为特征场与标号场的关系模型，特征场与标号场的关系模型ux(x,l)表示为：

其中，为异质性势能函数，xj＝{xs；s∈pj}表示规则子块pj内所有像素为s的强度的集合，表示图像域中像素为s、标号为lj的所有像素的强度的集合，其中，ls为划分的图像域中像素s的标号。

在未知同质区域内像素光谱测度的统计分布模型中，k-s距离能以函数间最大垂直距离为两种不同统计分布间的相似性测度。因此，本实施例以k-s距离为图像分割准则，用k-s距离定义该异质性势能函数其表示为：

其中，dks代表k-s距离，即直方图和的最大间距，和分别代表两个数据集合xj＝{xs；s∈pj}和的采样分布函数，和分别表示为：

其中，n1和n2分别为两个数据集合xj和中元素的个数，h为强度值索引，若对于n位图像，h∈[0,2ⁿ-1]；

步骤3.2：在划分的图像域上，建立标号场模型；

利用势能函数定义标号场模型，标号场模型ul(l)表示为：

其中，β为邻域子块的空间作用参数，npj为规则子块pj的八邻域规则子块的集合，j′为规则子块pj的邻域规则子块pj′的索引，lj′为规则子块pj′的标号；若lj＝lj′，则δ(lj,lj′)＝1；若lj≠lj′，则δ(lj,lj′)＝0；

步骤3.3：建立基于规则子块的sar图像分割模型；

结合步骤3.1中特征场与标号场的关系模型和步骤3.2中的标号场模型，定义图像分割的全局势能函数u(x,l)，其表示为：

利用非约束gibbs概率分布刻画上述全局势能函数u(x,l)，得到基于规则子块的sar图像分割模型g(x,l)，其表示为：

步骤4：针对已建立的sar图像分割模型，设定迭代次数并设计合理的移动操作，每次迭代中，遍历所有移动操作，利用rjmcmc算法，求解基于规则子块的sar图像分割模型。本实施例中，设定迭代次数为10000次。具体包括如下步骤：

步骤4.1：通过更新标号场中的标号更新标号场；

更新标号场中的标号的操作流程如图4所示，具体操作为：从图像域p＝{pj,j＝1,...,j}中随机抽取一个规则子块pj，其对应标号为lj；从图像的总类别数集合{1,…,k}中随机抽取规则子块pj的候选标号且利用非约束gibbs概率分布函数g(x,l)，得到更新标号场中标号lj为的接受率为：

其中，

从0～1中产生一个随机数，判断该随机数与接受率的大小，当该接受率大于该随机数时，接受本次更新标号场中标号的操作，否则放弃本次更新标号场中标号的操作。

步骤4.2：在步骤4.1更新标号场的基础上通过增加或减少规则子块个数来更新划分的图像域及标号场，更新图像域过程中随机选择增加规则子块或减少规则子块的操作；

增加规则子块个数是通过规则子块的分裂而实现的，具体操作过程如图5所示，具体步骤如下：

步骤4.2.1：从通过步骤4.1更新标号场后的由j个规则子块拟合而成的图像域p中随机抽取一个规则子块pj，其对应标号为lj；

步骤4.2.2：判断所选规则子块pj能否实现分裂操作；如果pj的像素数大于4且其行数或列数为2的整数倍数，则实现分裂pj操作，执行步骤4.2.3，增加图像域中规则子块个数；否则不分裂规则子块pj，不增加图像域中规则子块个数，结束本次分裂操作；

步骤4.2.3：将图像域p中能实现分裂操作的规则子块pj分裂成两个新的规则子块pj1和pj2，新的规则子块pj1的对应标号为新的规则子块pj2对应标号为且图像域中的规则子块数增加一个，图像域p更新为图像域j^*表示更新后的规则子块的总个数；未分裂的规则子块不变，对应标号也不变；图像域p更新为图像域p^*的接受率为：

as(p,p^*)＝min{1,rs}

其中，

其中，为更新后的图像域p^*的图像的一个实现，为更新后的标号场的实现，为更新后的标号场的标号；

步骤4.2.4：从0～1中产生一个随机数，判断该随机数与该接受率as(p,p^*)的大小，当该接受率大于该随机数时，接受本次通过增加规则子块来更新图像域及标号场的操作，否则放弃本次通过增加规则子块来更新图像域及标号场的操作。

减少规则子块个数是通过图像域p中随机抽取的一个规则子块和其任一个邻域规则子块合并成一个新的规则子块来实现的，该操作为分裂规则子块的对偶操作，因此，通过减少规则子块实现更新图像域的接受率为：

am(p,p^*)＝min{1,1/rs}

从0～1中产生一个随机数，判断该随机数与接受率am(p,p^*)的大小，当所述接受率大于该随机数时，接受本次通过减少规则子块来更新图像域及标号场的操作，否则放弃本次通过减少规则子块来更新图像域及标号场的操作。

步骤4.3：将步骤4.2中更新后的标号场的实现l代入非约束gibbs概率分布函数g(x,l)，得到非约束gibbs概率分布函数g(x,l)的一个函数值。

具体实施中，更新标号场的过程也可以先执行步骤4.2，后执行步骤4.1，将后执行的步骤4.1中更新后的标号场的实现l代入非约束gibbs概率分布函数g(x,l)。

步骤4.4：按设定的迭代次数重复执行步骤4.1至步骤4.3，得到非约束gibbs概率分布函数g(x,l)的一个函数值集合，集合中g(x,l)函数值最大时所对应的分割结果为基于规则子块的sar图像分割模型的最优解。

已知图像x求解l的过程，可直接通过最大化非约束gibbs概率分布函数得到。当非约束gibbs概率分布函数取得最大值时，基于规则子块的sar图像分割模型可得到最优解。

步骤5：输出分割结果图像。

采用本实施例所述的结合k-s距离和rjmcmc算法的sar图像分割方法分别对图2(a)、图2(b)和图2(c)所示的待分割图像进行图像分割，结果分别如图6(a)、图6(b)和图6(c)所示。图6(a)中，能将白色代表的溢油区域和黑色代表的水从sar图像中清晰的分割出来；图6(b)中，能将灰色代表的河岸、黑色代表的坚冰及白色代表的融冰从sar图像中清晰的分割出来；图6(c)中，能将灰色代表的水、黑色代表的坚冰及白色代表的融冰从sar图像中清晰的分割出来；这些分割结果表明分割模型区域一致性很高，并且保持了边缘的准确性；同时从待分割图像可以看出原图像斑点噪声明显，不易分割；从结果图可看出本发明采用k-s距离和规则划分技术相结合的分割模型，可以有效克制斑点噪声。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明权利要求所限定的范围。