CN103606164B - 基于高维三重马尔可夫场的sar图像分割方法 - Google Patents

Abstract

本发明公开了一种基于高维三重马尔可夫场的SAR图像分割方法，主要解决现有技术中分割区域一致性不足的问题，其实现步骤是：(1)输入SAR图像，并对其初始化分割获得分割标记场；(2)通过自协方差运算，建立附加标记场；(3)计算像素点分割标记可靠性参数；(4)计算像素点与其高维邻域内像素点的起伏相似性参数；(5)构建高维邻域下分割标记场、附加标记场以及SAR图像三者的联合分布概率；(6)基于联合分布概率，构建后验边缘概率分割模型；(7)最大化后验边缘概率更新分割标记场和附加标记场，获得最终分割结果。本发明提高了同质区域的分割一致性及异质区域的分割精度，可用于SAR图像目标检测与识别。

Description

基于高维三重马尔可夫场的SAR图像分割方法

技术领域

本发明属于图像处理技术领域，更进一步涉及SAR图像分割方法，可用于对目标检测和目标识别。

背景技术

合成孔径雷达(synthetic aperture radar,SAR)系统因其全天时，全天候，具有穿透性等特点而广泛应用于军事和民用邻域。SAR图像分割是SAR图像目标识别与解译技术的重要环节，它可以提供图像整体结构信息，揭示SAR图像本质,为SAR系统的自动目标识别建立基础,并推动SAR的应用。这一领域已逐渐成为近年来国内外的一个研究热点。然而SAR系统成像时发射相干电磁波照射目标会使图像产生大量的斑点噪声，使得传统的光学图像分割算法用于SAR图像分割时无法取得令人满意的结果。

在众多的SAR图像分割算法中，马尔可夫随机场（Markov random field,MRF）模型因其能够充分考虑图像局部相关性而受到广泛关注。近年来，MRF模型也在不断地发展和完善，涌现出多种扩展模型，包括树结构马尔可夫随机场，条件随机场，三重马尔可夫随机场。尽管MRF模型能够有效的描述分割结果的邻域相关性，但该模型并未充分考虑到SAR图像的非平稳特性。针对SAR图像非高斯、非平稳的统计特性，Benboudjema提出了三重马尔可夫场(Triplet Markov Fields，TMF)模型。TMF模型考虑了图像的非平稳特性，并且能够采用多样的统计模型对数据进行精确建模，在处理非平稳SAR图像分割问题时取得了令人满意的结果。

然而，受制于模型复杂度以及参数个数，TMF模型仅仅考虑传统的低维4邻域系统，在构建势能函数的过程中，并未充分考虑5×5以上的高维邻域系统内局部结构的方向信息与各向异性，而是仅局限于利用一阶导数进行统计建模。因此，TMF模型无法精确描述纹理结构丰富、关联性强的SAR图像，进而不能对其进行精确的统计建模，阻碍了分割精度的进一步提高。同时，基于低维4邻域的TMF模型仅考虑一阶邻域像素的相关性，无法克服SAR图像中斑点噪声的影响，阻碍了分割结果中区域一致性的提高。

发明内容

本发明的目的在于针对上述已有问题的缺点，提出一种高维三重马尔可夫场的SAR图像分割方法，以在保证分割准确性的同时，提高分割结果的区域一致性。

实现本发明目的的技术方案是：利用SAR图像像素点间的自协方差分析，获取SAR图像中高维邻域内局部结构信息，并依此重新定义附加标记场的意义，指导分割标记场与附加标记场在高维邻域条件下的联合分布概率的构建，依据贝叶斯概率推导，完成对SAR图像的分割。具体步骤包括如下：

(1)输入SAR图像Y，Y＝{Y_s|s∈Q}，Y_s为像素点s的灰度值，Y_s∈[0,1,...,255]，Q为SAR图像像素点集；

(2)对SAR图像Y进行初始分割，获得分割标记场X，X＝{X_s|s∈Q}，X_s为像素点s的分割标记值，X_s∈[1,2,...,K]，K为分割标记的总类别数，取值为正整数；

(3)对SAR图像Y进行自协方差运算，建立附加标记场U，U＝{U_s|s∈Q}，U_s为像素点s的附加标记值，U_s∈[0,1]，U_s为0表示像素点s处于同质性区域，U_s为1表示像素点s处于异质性区域；

(4)对于像素点s，计算其分割标记值X_s的可靠性参数a_s：

4a)在SAR图像Y中，获取像素点s的灰度值Y_s，获取像素点s高维邻域N_s内所有像素的灰度值计算Y_s和的联合分布概率

4b)将像素点s的灰度值人工设为Y_s ^*，Y_s ^*从集合[0,1,...,255]中取值，按照步骤4a)，计算Y_s ^*和的联合分布概率

4c)根据Y_s和的联合分布概率以及Y_s ^*和的联合分布概率计算在条件下Y_s的概率和在条件下Y_s ^*的概率

4d)根据上述得到的以及计算像素点s的分割标记值X_s的可靠性参数a_s：

$a_s=\exp [-\frac{\frac{1}{255}\underset{{Y_s}^{*}\≠Y_s}{\mathrm{\Σ}}p\left({Y_s}^{*}\right|Y_{N_s})}{p\left(Y_s\right|Y_{N_s})}],$ 其中，exp为指数运算符号；

4e)重复步骤4a)-4d)，计算每个像素点的分割标记可靠性参数{a_s|s∈Q}；

(5)对于像素点s及其高维邻域N_s内像素点t，计算两像素点之间的起伏程度相似性参数b_st：

5a)以像素点s为中心，提取第一起伏程度特征向量C_s；以像素点t为中心，提取第二起伏程度特征向量C_t；

5b)计算像素点s与像素点t之间的起伏程度相似性参数b_st＝exp(-dis(C_s,C_t))，其中，dis(·)为对输入的两个行向量进行欧氏距离计算；

5c)重复步骤5a)和5b)，计算每个像素点与其邻域像素点之间的起伏程度相似性参数：{b_st|t∈N_s,s∈Q}；

(6)计算高维邻域条件下分割标记场X和附加标记场U的联合分布先验概率p(X,U)：

6a)计算像素点s的局部高维邻域势能

$\begin{array}{c}{W}_s(X_s,U_s,X_{N_s},U_{N_s})=\\ \underset{t\∈N_s}{\mathrm{\Σ}}{b}_{\mathrm{st}}(1-2\mathrm{\δ}(U_s,U_t))+\left(\underset{t\∈N_s}{\mathrm{\Σ}}{a}_t(1-2\mathrm{\δ}(X_s,X_t))\right)\mathrm{\δ}(U_s,0)\\ +(\underset{t\∈N_s}{\mathrm{\Σ}}(a_t(1-2\mathrm{\δ}(X_s,X_t))\mathrm{\ϵ}(C(s,t)-\mathrm{th})+a_t(2\mathrm{\δ}(X_s,X_t)-1)\mathrm{\ϵ}(\mathrm{th}-C(s,t)))\mathrm{\δ}(U_s,1)\end{array}$

其中，th为自协方差门限，取值为0.5，当s处于同质区域时，δ(U_s,0)＝1，δ(U_s,1)＝0，当s处于异质区域时，δ(U_s,1)＝1，δ(U_s,0)＝0，ε(·)为单位阶跃函数，a_t为像素点t的分割标记可靠性参数，b_st为像素点s与像素点t之间的起伏程度相似性参数；

6b)根据局部高维邻域势能计算分割标记场X以及附加标记场U的整体势能： $W(X,U)=\underset{s\∈Q}{\mathrm{\Σ}}{W}_s(X_s,U_s,X_{N_s},U_{N_s});$

6c)根据整体势能W(X,U)，利用吉布斯随机场概率公式计算分割标记场X以及附加标记场U联合分布的先验概率p(X,U)：

p(X,U)＝γexp(-W(X,U))，其中，γ为归一化常数，

(7)利用先验概率p(X,U)计算分割标记场X、附加标记场U以及SAR图像Y三者的联合分布概率p(X,U,Y)；

(8)利用联合分布概率p(X,U,Y)，计算像素点s的分割标记类别后验边缘概率p(X_s|Y)，以及附加标记类别后验边缘概率p(U_s|Y)；

(9)利用贝叶斯最大后验边缘概率准则，分别确定使像素点s对应后验边缘概率p(X_s|Y)最大时的分割标记值，以及使像素点s对应后验边缘概率p(U_s|Y)最大时的附加标记值，并更新分割标记场和附加标记场；

(10)统计更新前后的分割标记场中类别发生变化的像素点数目，计算变化像素点数目与SAR图像总像素点数目的比值，将该比值作为终止的检测条件，若该比值大于设定的阈值ω=10^-6，返回步骤(6)，否则，输出分割标记场X作为最终的SAR图像Y的分割结果。

本发明与现有技术相比具有以下优点：

第一，本发明在对SAR图像分割的过程中，通过附加标记场对SAR图像中同质区域以及异质区域进行划分，并针对同质区域和异质区域，采用不同的函数形式进行局部高维邻域势能构建，与现有技术相比，本发明即能在异质区域得到较为精细的分割结果，又能在同质区域得到一致性较好的分割结果；

第二，本发明采用局部高维邻域势能进行分割标记场和附加标记场联合势能的构建，与现有技术采用的结构简单的低维4邻域局部势能相比，能够充分考虑高维邻域系统内的结构信息，在对纹理结构丰富且含有斑点噪声的SAR图像进行分割时得到更好的分割结果。

仿真结果表明，本发明与现有的传统三重马尔可夫场分割方法相比，既能在SAR图像异质区域得到更精细的分割结果，又能在SAR图像的同质区域得到一致性更好的分割结果。

附图说明

图1是本发明基于高维三重马尔可夫场的SAR图像分割的流程图；

图2是用本发明和现有三重马尔可夫场方法对机场SAR图像的分割结果图；

图3是用本发明和现有三重马尔可夫场方法对公路SAR图像的分割结果图；

图4是用本发明和现有三重马尔可夫场方法对军事基地SAR图像的分割结果图；

图5是用本发明和现有三重马尔可夫场方法对平原地区SAR图像的分割结果图。

具体实施方式

下面结合附图对本发明的具体实施和效果作进一步描述：

参照图1，本发明的实施步骤如下：

步骤1.输入一幅SAR图像Y。

该SAR图像Y为256灰度图像，每个像素点s的灰度值Y_s为从0到255中取正整数，将该SAR图像中的所有像素点记为集合Q，则该SAR图像Y可表示为Y＝{Y_s|s∈Q}，Y_s∈[0,1,...,255]。

步骤2.对SAR图像Y进行初始分割，获得分割标记场X。

对SAR图像进行初始分割有多种方法，包括聚类方法、分水岭方法、区域生长方法以及直方图阈值方法等，在本实例中，采用K-means聚类对SAR图像Y进行初始分割，即将SAR图像Y中所有像素点依据其灰度值自动划分为指定K类，根据聚类结果得到初始分割标记场X，人工设定分割类别K，则分割标记集合为[1,2,...,K]，对SAR图像Y中像素点s，判断其分割标记值X_s，该标记值从分割标记集合[1,2,...,K]中取值，对所有像素点对应的分割标记值进行判断，获得SAR图像Y对应的分割标记场X，表示为X＝{X_s|s∈Q}，X_s∈[1,2,...,K]。

步骤3.对SAR图像Y进行自协方差运算，得到所有像素点与其各自高维邻域内像素点的自协方差度量。

自协方差是对图像中任意两像素点进行相似性大小判断的度量，该度量越大，则说明两像素点之间的相似性越高，计算自协方差是成熟的现有方法，本实例采用现有的自协方差算法，对像素点s与其高维邻域N_s内像素点之间的自协方差度量C(s,t)进行按如下计算：

3.1)以像素点s为中心，提取7×7窗口内除像素点s以外的所有像素点作为像素点s的高维邻域N_s；在像素点s的高维邻域N_s内选取像素点t，以s为中心，提取大小为9×9的像素矩阵，并按行排列，得到第一行向量Y^s；以t为中心，提取大小为9×9的像素矩阵，并按行排列，得到第二行向量Y^t；

3.2)按照下式计算像素点s和像素点t的自协方差度量C(s,t)：

$C(s,t)=\frac{1}{81}(Y^s-{\mathrm{\μ}}_{Y^s})*{(Y^t-{\mathrm{\μ}}_{Y^t})}^T,$

其中，为第一行向量Y^s的均值，为第二行向量Y^t的均值，T表示对向量进行转置；

3.3)重复步骤3.1)和步骤3.2)得到所有像素点与其各自高维邻域内像素点的自协方差度量{C(s,t)|t∈N_s,s∈Q}。

步骤4.计算每个像素点的局部起伏特征，建立附加标记场U。

像素点的局部起伏特征反映了像素点所处局部区域图像像素灰度值的起伏性，对像素点局部起伏特征的计算可采用局部区域内像素灰度值的方差、局部区域内的灰度共生矩阵以及局部区域内自协方差度量集合的方差等方法。在本实例中，采用局部区域内自协方差度量集合的方差来计算每个像素点的局部起伏特征，并依据该特征，建立附加标记场，具体步骤如下：

4.1)计算像素点s与其高维邻域N_s内像素点t的自协方差度量C(s,t)；将t在N_s中遍历，得到s与其高维邻域N_s内所有像素点的自协方差度量集合F_s＝{C(s,t)|t∈N_s}；

4.2)计算自协方差度量集合F_s的方差，作为像素点s的局部起伏特征，方差越大，则局部起伏程度越高；

4.3)重复步骤4.1)和步骤4.2)，获得SAR图像Y中所有像素点各自对应的局部起伏特征，并对所有像素点的局部起伏特征进行归一化处理；

4.4)根据设定的门限β＝0.5，对所有像素点所处区域进行判定：

以像素点s为例，进行如下判定：

当像素点s的局部起伏特征大于等于β时，像素点s被判定为处于异质性区域；

当像素点s的局部起伏特征小于β时，像素点s被判定为处于同质性区域；

用该判定结果初始化附加标记场U，U_s为像素点s的附加标记值，U_s为0表示像素点s处于同质性区域，U_s为1表示像素点s处于异质性区域；

对SAR图像像素集合Q中所有像素点进行判定，获得附加标记场，该附加标记场可表示为U＝{U_s|s∈Q}，U_s∈[0,1]。

步骤5.计算SAR图像像素点集Q中每个像素点的分割标记可靠性参数{a_s|s∈Q}。

5.1)在SAR图像Y中，分别获取像素点s的灰度值Y_s和像素点s高维邻域N_s内所有像素的灰度值计算Y_s和的联合分布概率

计算Y_s和的联合分布概率可采用多种现有方法，例如基于主成分分析的图像块建模方法、基于专家乘积的图像块建模方法以及基于Gabor滤波器的图像块建模方法。本实例采用基于Gabor滤波器的建模方法来计算Y_s和的联合分布概率其计算步骤如下：

首先，将像素点s的灰度值Y_s与像素点s高维邻域N_s内所有像素的灰度值组成大小为7×7的像素矩阵Y^(s)，采用大小为7×7、方向依次是0,30,60,90,120,150度的二维Gabor滤波器组J，对该像素矩阵Y^(s)进行滤波；

其次，根据滤波器组J中每一个滤波器的输出结果，按照下式计算所述Y_s和的联合分布概率

$p(Y_s,Y_{N_s})=\frac{1}{Z}\underset{i=1}{\overset{6}{\mathrm{\Π}}}(\underset{V_D}{\mathrm{\Σ}}{\mathrm{\λ}}_{V_D}^i*{\mathrm{\δ}}^{*}(V_D,D(\mathrm{row}\left(J_i\right)*\mathrm{col}\left(Y^{\left(s\right)}\right)))),$

其中，J_i为滤波器组J中的第i个滤波器，row(·)为将第i个滤波器J_i按行排列的行向量，col(·)将Y^(s)按列排列的列向量，D(·)为将滤波器与像素矩阵的滤波结果离散化成间隔为1的整数值，V_D为整数值，取值范围为[-200，200]，为第i个滤波器J_i对Y^(s)滤波结果为V_D的概率，可由该滤波器对SAR图像Y的滤波结果直方图计算得到，δ^*(·)为指示函数，当函数中两变量相等时，函数值为1，其它为0，Z为归一化常数；

5.2)将像素点s的灰度值人工设为Y_s ^*，Y_s ^*从集合[0,1,...,255]中取值，按照步骤5.1），计算Y_s ^*和的联合分布概率

5.3)根据Y_s和的联合分布概率以及Y_s ^*和的联合分布概率按照下式计算在条件下Y_s的概率和在条件下Y_s ^*的概率

$p\left(Y_s\right|Y_{N_s})=\frac{p(Y_s,Y_{N_s})}{\underset{{Y_s}^{*}=0}{\overset{255}{\mathrm{\Σ}}}p({Y_s}^{*},Y_{N_s})},$

$p\left({Y_s}^{*}\right|Y_{N_s})=\frac{p({Y_s}^{*},Y_{N_s})}{\underset{{Y_s}^{*}=0}{\overset{255}{\mathrm{\Σ}}}p({Y_s}^{*},Y_{N_s})};$

5.4)根据上述得到的以及按照下式计算像素点s的分割标记值X_s的可靠性参数a_s：

$a_s=\exp [-\frac{\frac{1}{255}\underset{{Y_s}^{*}\≠Y_s}{\mathrm{\Σ}}p\left({Y_s}^{*}\right|Y_{N_s})}{p\left(Y_s\right|Y_{N_s})}],$ 其中，exp为指数运算符号；

5.5)重复步骤5.1)-步骤5.4)，获得SAR图像像素点集Q中每个像素点的分割标记可靠性参数{a_s|s∈Q}。

步骤6.计算SAR图像像素点集Q中每个像素点与其高维邻域内像素点的起伏程度相似性参数{b_st|t∈N_s,s∈Q}。

6.1)在SAR图像Y中，以像素点s为中心，提取大小为7×7的第一像素矩阵Y^(s)；以像素点t为中心，提取大小为7×7的第二像素矩阵Y^(t)；

6.2)按行计算第一像素矩阵Y^(s)中每一像素点与第一像素矩阵Y^(s)的中心像素点s的自协方差度量，并将每个度量值按行排列，得到49维第一行向量C_s；

6.3)按行计算第二像素矩阵Y^(t)中每一像素点与第二像素矩阵Y^(t)的中心像素点t的自协方差度量，并将每个度量值按行排列，得到49维第二行向量C_t；

6.4)按下式计算像素点s与像素点t之间的起伏程度相似性参数b_st：

b_st＝exp(-dis(C_s,C_t))，其中，dis(·)为对输入的两个行向量进行欧氏距离计算；

6.5)重复步骤6.1)-步骤6.4)，计算Q中每个像素点与其高维邻域内像素点之间的起伏程度相似性参数：{b_st|t∈N_s,s∈Q}。

步骤7.计算高维邻域条件下分割标记场X和附加标记场U的联合分布先验概率p(X,U)。

7.1)按照下式计算像素点s的局部高维邻域势能

$\begin{array}{c}{W}_s(X_s,U_s,X_{N_s},U_{N_s})=\\ \underset{t\∈N_s}{\mathrm{\Σ}}{b}_{\mathrm{st}}(1-2\mathrm{\δ}(U_s,U_t))+\left(\underset{t\∈N_s}{\mathrm{\Σ}}{a}_t(1-2\mathrm{\δ}(X_s,X_t))\right)\mathrm{\δ}(U_s,0)\\ +(\underset{t\∈N_s}{\mathrm{\Σ}}(a_t(1-2\mathrm{\δ}(X_s,X_t))\mathrm{\ϵ}(C(s,t)-\mathrm{th})+a_t(2\mathrm{\δ}(X_s,X_t)-1)\mathrm{\ϵ}(\mathrm{th}-C(s,t)))\mathrm{\δ}(U_s,1)\end{array}$

其中，th为自协方差门限，取值为0.5，当s处于同质区域时，δ(U_s,0)＝1，δ(U_s,1)＝0，当s处于异质区域时，δ(U_s,1)＝1，δ(U_s,0)＝0，ε(·)为单位阶跃函数；

7.2)根据局部高维邻域势能按照下式计算分割标记场X以及附加标记场U的整体势能： $W(X,U)=\underset{s\∈Q}{\mathrm{\Σ}}{W}_s(X_s,U_s,X_{N_s},U_{N_s});$

7.3)根据整体势能W(X,U)，利用吉布斯随机场概率公式计算分割标记场X以及附加标记场U联合分布的先验概率p(X,U)：

p(X,U)＝γexp(-W(X,U))，其中，γ为归一化常数，

步骤8.在获得先验概率p(X,U)的基础上，计算分割标记场X、附加标记场U以及SAR图像Y三者的联合分布概率p(X,U,Y)。

8.1)按照下式计算像素点s的灰度值Y_s在其对应分割类别X_s条件下的似然概率p(Y_s|X_s)：

$p\left(Y_s\right|X_s)=\frac{1}{\mathrm{\Γ}\left(L\right)}\left(\frac{L}{{\mathrm{\σ}}_{X_s}}\right){Y_s}^{L-1}\exp (-\frac{{\mathrm{LY}}_s}{{\mathrm{\σ}}_{X_s}}),$

其中，Γ(·)为伽马函数，L为SAR图像的等效视数，为第X_s类的灰度均值；

8.2)根据整体势能W(X,U)以及似然概率p(Y_s|X_s)，按照下式计算分割标记场X、附加标记场U以及SAR图像Y三者的联合分布概率p(X,U,Y)：

$p(X,U,Y)=\mathrm{\γexp}[-W(X,U)-\underset{s\∈Q}{\mathrm{\Σ}}\log p\left(Y_s\right|X_s)].$

步骤9.利用联合分布概率p(X,U,Y)，计算SAR图像像素点集Q中所有像素点的分割标记类别后验边缘概率{p(X_s|Y)|s∈Q}，以及附加标记类别后验边缘概率{p(U_s|Y)|s∈Q}.

9.1)对联合概率p(X,U,Y)进行20次吉布斯采样运算，得到SAR图像像素点集Q中所有像素点的分割标记类别和附加标记类别的联合后验边缘概率{p(X_s,U_s|Y)|s∈Q}；

9.2)根据联合后验边缘概率p(X_s,U_s|Y)，分别按照下式计算所有像素点的分割标记类别后验边缘概率p(X_s|Y)，以及附加标记类别后验边缘概率p(U_s|Y)|：

$p\left(X_s\right|Y)=\underset{U_s}{\mathrm{\Σ}}p(X_s,U_s|Y),$

$p\left(U_s\right|Y)=\underset{X_s}{\mathrm{\Σ}}p(X_s,U_s|Y),s\∈Q.$

步骤10.更新分割标记场和附加标记场。

10.1)利用贝叶斯最大后验边缘概率准则，确定使SAR图像像素点集Q中每个像素点对应后验边缘概率p(X_s|Y),s∈Q最大时的分割标记值，更新分割标记场；

10.2)利用贝叶斯最大后验边缘概率准则，确定使SAR图像像素点集Q中每个像素点对应后验边缘概率p(U_s|Y),s∈Q最大时的附加标记值，更新附加标记场。

步骤11.统计更新前后的分割标记场中类别发生变化的像素点数目，计算变化像素点数目与SAR图像总像素点数目的比值，将该比值作为终止的检测条件，若该比值大于设定的阈值ω=10^-6，返回步骤7，否则，输出分割标记场X作为最终的SAR图像Y的分割结果。

本发明效果可以通过以下仿真进一步说明：

1.仿真条件

本发明的仿真是在主频2.5GHz Intel(R)Pentium(R)2CPU的硬件环境以及MATLABR2009b，Window XP Professional的软件环境下进行的。

2.仿真内容与结果分析

仿真1，采用本发明的方法和现有三重马尔可夫场方法对实测机场SAR图像进行分割，仿真结果如图2所示。其中图2(a)为待分割的机场SAR图像，大小为256×256，图2(b)为利用现有三重马尔可夫场方法对图2(a)的分割结果，图2(c)为利用本发明方法对图2(a)的分割结果。

图2(a)中跑道部分纹理较为简单，包含有大量的同质区域，在建筑群部分，图像细节信息增多，纹理信息较为复杂。从仿真1结果可得，在图2(c)中，SAR图像跑道部分的分割结果相较于图2(b)具有更好的区域一致性，同时，建筑群部分的分割相较于图2(b)更加精细，不同走向的白色线状目标能够准确的分割出来，而在图2(b)中，白色线状目标被分割成连通起来的区域。

仿真2，用本发明的方法和现有三重马尔可夫场方法对实测公路SAR图像进行分割，仿真结果如图3所示。其中图3(a)为待分割的公路SAR图像，大小为415×438，图3(b)为利用现有三重马尔可夫场方法对图3(a)的分割结果，图3(c)为利用本发明方法对图3(a)的分割结果。

仿真3，采用本发明的方法和现有三重马尔可夫场方法对实测军事基地SAR图像进行分割，仿真结果如图4所示。其中图4(a)为待分割的军事基地SAR图像，大小为317×324，图4(b)为利用现有三重马尔可夫场方法对图4(a)的分割结果，图4(c)为利用本发明方法对图4(a)的分割结果。

仿真4，采用本发明的方法和现有三重马尔可夫场方法对实测平原地区SAR图像进行分割，仿真结果如图5所示。其中图5(a)为待分割的平原地区SAR图像，大小为256×256，图5(b)为利用现有三重马尔可夫场方法对图5(a)的分割结果，图5(c)为利用本发明方法对图5(a)的分割结果。

以上四个仿真表明：本发明方法由于考虑了SAR图像高维邻域内结构信息，从而能够从区域同质性或异质性的不同角度出发，构造恰当的高维势能，同时提高了对SAR图像同质区域分割的一致性，以及对SAR图像异质区域分割的精细度。

Claims (9)

1.一种基于高维三重马尔可夫场的SAR图像分割方法，包括如下步骤：

(1)输入SAR图像Y，Y＝{Y_s|s∈Q}，Y_s为像素点s的灰度值，Y_s∈[0,1,...,255]，Q为SAR图像像素点集；

(2)对SAR图像Y进行初始分割，获得分割标记场X，X＝{X_s|s∈Q}，X_s为像素点s的分割标记值，X_s∈[1,2,...,K]，K为分割标记的总类别数，取值为正整数；

(3)对SAR图像Y进行自协方差运算，建立附加标记场U，U＝{U_s|s∈Q}，U_s为像素点s的附加标记值，U_s∈[0,1]，U_s为0表示像素点s处于同质性区域，U_s为1表示像素点s处于异质性区域；

(4)对于像素点s，计算其分割标记值X_s的可靠性参数a_s：

4a)在SAR图像Y中，获取像素点s的灰度值Y_s，获取像素点s高维邻域N_s内所有像素的灰度值计算Y_s和的联合分布概率

4b)将像素点s的灰度值人工设为Y_s ^*，Y_s ^*从集合[0,1,...,255]中取值，按照步骤4a)，计算Y_s ^*和的联合分布概率

4c)根据Y_s和的联合分布概率以及Y_s ^*和的联合分布概率计算在条件下Y_s的概率和在条件下Y_s ^*的概率

4d)根据上述得到的以及计算像素点s的分割标记值X_s的可靠性参数a_s：

其中，exp为指数运算符号；

4e)重复步骤4a)-4d)，计算每个像素点的分割标记值的可靠性参数{a_s|s∈Q}；

(5)对于像素点s及其高维邻域N_s内像素点t，计算两像素点之间的起伏程度相似性 参数b_st：

5a)以像素点s为中心，提取第一起伏程度特征向量C_s；以像素点t为中心，提取第二起伏程度特征向量C_t；

5b)计算像素点s与像素点t之间的起伏程度相似性参数b_st＝exp(-dis(C_s,C_t))，其中，dis(·)为对输入的两个特征向量进行欧氏距离计算；

5c)重复步骤5a)和5b)，计算每个像素点与其邻域像素点之间的起伏程度相似性参数：{b_st|t∈N_s,s∈Q}；

(6)计算高维邻域条件下分割标记场X和附加标记场U的联合分布先验概率p(X,U)：

6a)计算像素点s的局部高维邻域势能

其中，th为自协方差门限，取值为0.5，当s处于同质区域时，δ(U_s,0)＝1，δ(U_s,1)＝0，当s处于异质区域时，δ(U_s,1)＝1，δ(U_s,0)＝0，ε(·)为单位阶跃函数，a_t为像素点t的分割标记值的可靠性参数，b_st为像素点s与像素点t之间的起伏程度相似性参数，C(s,t)为像素点s与像素点t的自协方差度量；

6b)根据局部高维邻域势能计算分割标记场X以及附加标记场U的整体势能：

6c)根据整体势能W(X,U)，利用吉布斯随机场概率公式计算分割标记场X以及附加标记场U联合分布的先验概率p(X,U)：

p(X,U)＝γexp(-W(X,U))，其中，γ为归一化常数，(7)利用先验概率p(X,U)计算分割标记场X、附加标记场U以及SAR图像Y三者的 联合分布概率p(X,U,Y)；

(8)利用联合分布概率p(X,U,Y)，计算像素点s的分割标记类别后验边缘概率p(X_s|Y)，以及附加标记类别后验边缘概率p(U_s|Y)；

(9)利用贝叶斯最大后验边缘概率准则，分别确定使像素点s对应后验边缘概率p(X_s|Y)最大时的分割标记值，以及使像素点s对应后验边缘概率p(U_s|Y)最大时的附加标记值，并更新分割标记场和附加标记场；

(10)统计更新前后的分割标记场中类别发生变化的像素点数目，计算变化像素点数目与SAR图像总像素点数目的比值，将该比值作为终止的检测条件，若该比值大于设定的阈值ω＝10^-6，返回步骤(6)，否则，输出分割标记场X作为最终的SAR图像Y的分割结果。

2.根据权利要求1所述的基于高维三重马尔可夫场的SAR图像分割方法，其特征在于：步骤(2)所述的初始分割，是指对SAR图像Y进行K-means聚类分割，将SAR图像所有像素点依据其灰度特征自动划分为指定K类，根据聚类结果得到初始分割标记场X，X＝{X_s|s∈Q}，X_s为像素点s的分割标记值，X_s∈[1,2,...,K]，K为分割标记的总类别数，取值为正整数。

3.根据权利要求1所述的基于高维三重马尔可夫场的SAR图像分割方法，其特征在于：步骤(3)所述的对SAR图像Y进行自协方差运算，按如下步骤进行：

3a)以像素点s为中心，提取7×7窗口内除像素点s以外的所有像素点作为像素点s的高维邻域N_s；在像素点s的高维邻域N_s内选取像素点t，以s为中心，提取大小为9×9的像素矩阵，并按行排列，得到第一行向量Y^s；以t为中心，提取大小为9×9的像素矩阵，并按行排列，得到第二行向量Y^t；

3b)计算像素点s和像素点t的自协方差度量C(s,t)：

其中，为第一行向量Y^s的均值，为第二行向量Y^t的均值，T表示对向量进行 转置；

3c)重复步骤3a)和3b)，计算所有像素点与其各自高维邻域内像素点的自协方差度量{C(s,t)|t∈N_s,s∈Q}，完成对SAR图像Y进行自协方差运算。

4.根据权利要求1所述的基于高维三重马尔可夫场的SAR图像分割方法，其特征在于：步骤(3)所述的建立附加标记场U，按如下步骤进行：

3d)在SAR图像Y中，以像素点s为中心，提取7×7窗口内除像素点s以外的所有像素点作为像素点s的高维邻域N_s，计算像素点s与其高维邻域N_s内像素点t的自协方差度量C(s,t)；将t在N_s中遍历，得到s与其高维邻域N_s内所有像素点的自协方差度量集合F_s＝{C(s,t)|t∈N_s}；

3e)计算自协方差度量集合F_s的方差，作为像素点s的局部起伏特征，方差越大，则局部起伏程度越高；

3f)重复步骤3d)以及3e)，获得SAR图像Y中所有像素点各自对应的局部起伏特征，并对局部起伏特征进行归一化处理；

3g)根据设定的门限β＝0.5，对任意像素点s所处区域为同质性或异质性进行判定：

当像素点s的局部起伏特征大于等于β时，像素点s被判定为处于异质性区域，当像素点s的局部起伏特征小于β时，像素点s被判定为处于同质性区域；

3h)用该判定结果初始化附加标记场U＝{U_s|s∈Q}，

其中，U_s为像素点s的附加标记值，U_s∈[0,1]，U_s为0表示像素点s处于同质性区域，U_s为1表示像素点s处于异质性区域。

5.根据权利要求1所述的基于高维三重马尔可夫场的SAR图像分割方法，其特征在于：所述步骤4a)中计算Y_s和的联合分布概率按如下步骤进行：

4a1)将像素点s的灰度值Y_s与像素点s高维邻域N_s内所有像素的灰度值组成大小为7×7的像素矩阵Y^(s)，采用大小为7×7、方向依次是0,30,60,90,120,150度的二维Gabor 滤波器组J，对该像素矩阵Y^(s)进行滤波；

4a2)根据滤波器组J中每一个滤波器的输出结果，计算所述Y_s和的联合分布概率

其中，J_i为滤波器组J中的第i个滤波器，row(J_i)为第i个滤波器J_i按行排列的行向量，col(Y^(s))为将Y^(s)按列排列的列向量，D(·)为将滤波器与像素矩阵的滤波结果离散化成间隔为1的整数值，V_D为整数值，取值范围为[-200，200]，为第i个滤波器J_i对Y^(s)滤波结果为V_D的概率，δ^*(·)为指示函数，当函数中两变量相等时，函数值为1，其它为0，Z为归一化常数。

6.根据权利要求1所述的基于高维三重马尔可夫场的SAR图像分割方法，其特征在于：步骤4c)所述的计算在条件下Y_s的概率和在条件下Y_s ^*的概率 通过如下公式计算：

其中，p(|)为条件概率，为Y_s与的联合概率，为Y_s ^*与的联合概率。

7.根据权利要求1所述的基于高维三重马尔可夫场的SAR图像分割方法，其特征在于：步骤5a)所述的以像素点s为中心，提取第一起伏程度特征向量C_s；以像素点t为中心，提取第二起伏程度特征向量C_t，按如下步骤提取：

5a1)在SAR图像Y中，以像素点s为中心，提取大小为7×7的第一像素矩阵Y^(s)；以像素点t为中心，提取大小为7×7的第二像素矩阵Y^(t)；

5a2)按行计算第一像素矩阵Y^(s)中每一像素点与第一像素矩阵Y^(s)的中心像素点s的自协方差度量，并将每个度量值按行排列，得到49维第一行向量C_s；

5a3)按行计算第二像素矩阵Y^(t)中每一像素点与第二像素矩阵Y^(t)的中心像素点t的自协方差度量，并将每个度量值按行排列，得到49维第二行向量C_t。

8.根据权利要求1所述的基于高维三重马尔可夫场的SAR图像分割方法，其特征在于：所述步骤(7)中计算分割标记场X、附加标记场U以及SAR图像Y三者的联合分布的概率p(X,U,Y)，按如下步骤进行：

7a)计算像素点s的灰度值Y_s在其对应分割类别X_s条件下的似然概率p(Y_s|X_s)：

其中，Γ(·)为伽马函数，L为SAR图像的等效视数，为第X_s类的灰度均值；

7b)根据整体势能W(X,U)以及似然概率p(Y_s|X_s)，计算分割标记场X、附加标记场U以及SAR图像Y三者的联合分布概率p(X,U,Y)：

其中，W(X,U)为分割标记场X以及附加标记场U的整体势能，Q为SAR图像像素点集，exp为指数运算符号，γ为归一化常数，

9.根据权利要求1所述的基于高维三重马尔可夫场的SAR图像分割方法，其特征在于：步骤(8)所述的计算像素点s的分割标记类别后验边缘概率p(X_s|Y)，以及附加标记类别后验边缘概率p(U_s|Y)，按如下步骤进行：

8a)对联合概率p(X,U,Y)进行20次吉布斯采样运算，得到像素点s的分割标记类别 和附加标记类别的联合后验边缘概率p(X_s,U_s|Y)；

8b)根据联合后验边缘概率p(X_s,U_s|Y)，分别计算像素点s的分割标记类别后验边缘概率p(X_s|Y)，以及附加标记类别后验边缘概率p(U_s|Y)：