CN103914846B - 基于三重判别随机场的sar图像分割方法 - Google Patents

Abstract

本发明公开一种基于三重判别随机场模型的SAR图像分割方法，克服了现有技术采用三重马尔可夫随机场TMRF模型在捕获观测数据的统计相关性时，其作用距离必须依赖于所采用的邻域系统的缺点，本发明包括如下步骤：1.输入待分割SAR图像；2.提取特征和初始化；3.获得联合概率分布；4.获得后验边缘概率分布；5.分割SAR图像；6.估计模型参数；7.输出分割图像。本发明通过引入SAR图像的非平稳信息和像素点观测值之间的相关信息，保持了分割区域一致性，提高了分割结果边缘的准确性。

Description

基于三重判别随机场的SAR图像分割方法

技术领域

本发明属于图像处理技术领域，更进一步涉及雷达图像处理技术领域中的一种基于三重判别随机场的SAR图像分割方法。本发明可用于对合成孔径雷达SAR图像进行分割处理。

背景技术

合成孔径雷达SAR是一种高分辨率成像雷达。其在民用和军用中的使用需要合成孔径雷达SAR图像解译技术作为支撑，合成孔径雷达SAR图像分割是合成孔径雷达SAR图像解译技术的重要环节，它可以提供合成孔径雷达SAR图像的整体结构信息，揭示合成孔径雷达SAR图像的本质。随机场模型被认为是处理合成孔径雷达SAR图像分割问题的强有力工具，其优势在于能够在图像分类过程中充分考虑各像素间的空域相关性。目前，基于随机场的图像分割效果较好的方法有三重马尔可夫随机场TMRF的方法和条件随机场的方法。

西安电子科技大学拥有的专利技术“基于模糊的三马尔可夫场SAR图像分割方法”（申请号：201110356011X，授权公告号：CN102496142B）提出利用三重马尔可夫随机场TMRF模型进行图像分割的方法。该专利技术通过辅助场的引入描述了图像的非平稳特性，克服了隐马尔可夫随机场HMRF模型不具备描述图像非平稳特性的能力，实现对非平稳图像空域结构的精确建模。辅助场的不同取值描述了图像在不同平稳状态下的分布。三重马尔可夫随机场TMRF模型明确地考虑了图像的非平稳性，对合成孔径雷达SAR图像的非平稳性实现了准确的建模，因此三重马尔可夫随机场TMRF模型适合处理合成孔径雷达SAR图像分割问题。但是该专利技术存在的不足是，作为生成模型，三重马尔可夫随机场TMRF模型在捕获观测数据的统计相关性时，其作用距离必须依赖于所采用的邻域系统。

Takahiro Toyoda等人发表的论文“Random Field Model for Integration ofLocal Information and Global Information”（Trans.IEEE PAMI,2008,8:1483-1489）提出利用条件随机场进行图像分割的方法。该论文引用的条件随机场模型隶属于判别模型，其利用图像的统计特征直接对后验概率进行建模，充分利用观测数据的统计相关特性来提高模型精度，且条件随机场模型构建的后验概率分布满足条件马尔可夫性。条件随机场模型具备对图像特征有效统计建模的能力，能够精确地描述图像信息。但是该论文存在的不足是，此类模型需要训练数据学习模型参数，只限于处理有监督的图像分割问题，而且此类模型并未明确地从空域结构的角度上考虑图像的非平稳性。

发明内容

本发明的目的在于克服现有技术的不足，提出了一种基于三重判别随机场的SAR图像分割方法，利用判别随机场对待分割合成孔径雷达SAR图像进行建模，并引入待分割合成孔径雷达SAR图像的非平稳信息，改进了三重马尔可夫场，没有引入待分割合成孔径雷达SAR图像中像素点观测值之间的相关信息的问题，在保证分割区域一致性的同时，提高了分割结果边缘的准确性。

实现本发明目的的具体思路是，采用判别随机场对用于判断分割类别的联合概率分布进行建模，在建模过程中引入图像的非平稳信息，实现三重判别随机场模型的建立，用该三重判别随机场模型计算用于分割的后验边缘概率分布，利用该后验边缘概率分布分割待分割的合成孔径雷达SAR图像，得到分割结果图像。

本发明的具体步骤包括如下：

(1)输入待分割SAR图像；

(2)提取特征和初始化：

(2a)将半径为7的窗口在待分割的合成孔径雷达SAR图像内滑动，采用统计方法，计算滑窗内待分割的合成孔径雷达SAR图像像素点灰度值的均值、方差和熵三个图像统计特征，利用半方差方法，提取滑窗内待分割的合成孔径雷达SAR图像的四个半方差特征；

(2b)将半径为15的窗口在待分割的合成孔径雷达SAR图像内滑动，采用几何矩估计方法，计算滑窗内待分割合成孔径雷达SAR图像的非平稳各向异性高斯核参数；

(2c)采用直方图阈值方法，对待分割的合成孔径雷达SAR图像进行初始分割，得到待分割的合成孔径雷达SAR图像的初始分割图像；

(2d)采用模糊C均值聚类方法，对待分割合成孔径雷达SAR图像的平稳态进行初始分类，得到待分割的合成孔径雷达SAR图像的初始平稳态分类；

(2e)利用最小二乘法，根据待分割的合成孔径雷达SAR图像的初始分割图像和待分割的合成孔径雷达SAR图像的初始平稳态分类，估计用于判定分割类别的邻域内像素对的相互势能参数；

(2f)将用于判定分割类别的联合势能函数的参数全部初始化为1；

(2g)将当前的图像分割迭代次数设为1，将图像分割的最大迭代次数设为20；

(3)获得联合概率分布：

(3a)按照下式，构建用于判定分割类别的联合势能：

$A=\underset{n}{\mathrm{\Σ}}\underset{m}{\mathrm{\Σ}}\left(\mathrm{\δ}(x_s,n)\mathrm{\δ}(u_s,m)\right)\log G$

其中，A表示用于判定分割类别的联合势能；Σ表示求和操作；n表示待分割合成孔径雷达SAR图像的分割值；m表示待分割合成孔径雷达SAR图像的平稳态；δ(x_s,n)表示一个常数，当分割图像中s坐标位置的像素点和待分割合成孔径雷达SAR图像的分割值n相同时，δ(x_s,n)＝1，当分割图像中s坐标位置的像素点和待分割合成孔径雷达SAR图像的分割值n不同时，δ(x_s,n)＝0；x_s表示分割图像s坐标位置的像素点；s表示像素点在待分割合成孔径雷达SAR图像中的坐标位置；δ(u_s,m)表示一个常数，当待分割合成孔径雷达SAR图像中s坐标位置的像素点的平稳态和待分割合成孔径雷达SAR图像的平稳态m相同时，δ(u_s,m)＝1，当待分割合成孔径雷达SAR图像中s坐标位置的像素点的平稳态和待分割合成孔径雷达SAR图像的平稳态m不同时，δ(u_s,m)＝0；u_s表示待分割合成孔径雷达SAR图像中s坐标位置的像素点的平稳态；log表示求对数操作；G表示待分割图像中s坐标位置的像素点为分割值n，且待分割合成孔径雷达SAR图像中s坐标位置的像素点的平稳态为平稳态m的概率；

(3b)按照下式，构建用于判定分割类别的邻域内像素对的相互势能：

$I=\underset{s}{\mathrm{\Σ}}\underset{t\∈N_s}{\mathrm{\Σ}}(2\mathrm{\δ}(x_s,x_t)-1)\×(\mathrm{\α}(1-\mathrm{\δ}(u_s,u_t))+\underset{m}{\mathrm{\Σ}}{\mathrm{\α}}_m{\mathrm{\δ}}^{*}(u_s,u_t,m))\mathrm{\μ}$

其中，I表示用于判定分割类别的邻域内像素对的相互势能；s表示像素点在待分割合成孔径雷达SAR图像中的坐标位置；Σ表示求和操作；t表示在待分割合成孔径雷达SAR图像中坐标位置为s的像素点的邻域系统内的像素点；N_s表示待分割合成孔径雷达SAR图像中s坐标位置的像素点的邻域系统；δ(x_s,x_t)表示一个常数，当分割图像中s坐标位置的像素点和分割图像中t坐标位置的像素点相同时，δ(x_s,x_t)＝1，当分割图像中s坐标位置的像素点和分割图像中t坐标位置的像素点不同时，δ(x_s,x_t)＝1；x_s表示分割图像s坐标位置的像素点；x_t表示分割图像t坐标位置的像素点；α表示平稳态不同时，待分割合成孔径雷达SAR图像中s坐标位置的像素点和待分割合成孔径雷达SAR图像中t坐标位置的像素点之间相关性的参数；δ(u_s,u_t)表示一个常数，当待分割合成孔径雷达SAR图像中s坐标位置的像素点的平稳态和待分割合成孔径雷达SAR图像中t坐标位置的像素点的平稳态相同时，δ(u_s,u_t)＝1，当待分割合成孔径雷达SAR图像中s坐标位置的像素点的平稳态和待分割合成孔径雷达SAR图像中t坐标位置的像素点的平稳态不同时，δ(u_s,u_t)＝0；u_s表示待分割合成孔径雷达SAR图像中s坐标位置的像素点的平稳态；u_t表示待分割合成孔径雷达SAR图像中t坐标位置的像素点的平稳态；m表示待分割合成孔径雷达SAR图像的平稳态；α_m表示在平稳态为m时，待分割合成孔径雷达SAR图像中s坐标位置的像素点和待分割合成孔径雷达SAR图像中t坐标位置的像素点的相关性的参数；δ^*(u_s,u_t,m)表示一个常数，当待分割合成孔径雷达SAR图像中s坐标位置的像素点的平稳态、待分割合成孔径雷达SAR图像中t坐标位置的像素点的平稳态和平稳态m均相同时，δ^*(u_s,u_t,m)＝1，当待分割合成孔径雷达SAR图像中s坐标位置的像素点的平稳态、待分割合成孔径雷达SAR图像中t坐标位置的像素点的平稳态和平稳态m不同时，δ^*(u_s,u_t,m)＝0；μ表示待分割合成孔径雷达SAR图像中s坐标位置的像素点的特征和待分割合成孔径雷达SAR图像中t坐标位置的像素点的特征，组成的特征向量；

(3c)按照下式，构建用于判定分割类别的似然分布：

$f=\underset{s}{\mathrm{\Π}}\frac{\mathrm{\β}{y_s}^{\mathrm{\β\λ}-1}}{{\mathrm{\α}}^{\mathrm{\β\λ}}\mathrm{\Γ}}\exp (-{\left(\frac{y_s}{\mathrm{\α}}\right)}^{\mathrm{\β}})$

其中，f表示用于判定分割类别的似然分布；Π表示求乘积操作；s表示像素点在待分割合成孔径雷达SAR图像中的坐标位置；β表示广义伽马Gamma分布的形状参数；y_s表示待分割合成孔径雷达SAR图像中s坐标位置的像素点；λ表示广义伽马Gamma分布的指示形状参数；α表示广义伽马Gamma分布的尺度参数；Γ表示伽马Gamma函数；

(3d)按照下式，构建用于判定分割类别的联合概率分布：

$p=\frac{1}{Z}\exp (\underset{s}{\mathrm{\Σ}}A+\underset{s}{\mathrm{\Σ}}\underset{t\∈N_s}{\mathrm{\Σ}}I+\underset{s}{\mathrm{\Σ}}\log f)$

其中，p表示用于判定分割类别的联合概率分布；Z表示归一化因子；exp表示求指数操作；Σ表示求和操作；s表示像素点在待分割合成孔径雷达SAR图像中的坐标位置；A表示用于判定分割类别的联合势能；t表示在待分割合成孔径雷达SAR图像中s坐标位置的像素点的邻域系统内像素点的坐标位置；N_s表示待分割合成孔径雷达SAR图像中s坐标位置的像素点的邻域系统；I表示用于判定分割类别的邻域内像素对的相互势能；f表示用于判定分割类别的似然分布；

(4)获得后验边缘概率分布：

对用于判定分割类别的联合概率分布，利用吉布斯Gibbs采样方法，得到用于判定分割类别的后验边缘概率分布；

(5)分割SAR图像：

利用贝叶斯最大后验边缘概率准则，确定每个像素点新的标号，逐点更新待分割合成孔径雷达SAR图像的标号，得到分割图像；

(6)估计模型参数：

(6a)在待分割合成孔径雷达SAR图像的分割图像上，利用对比发散方法，估计用于判定分割类别的联合势能参数；

(6b)利用最小二乘法，根据待分割合成孔径雷达SAR图像的分割图像，估计用于判定分割类别的邻域内像素对的相互势能参数；

(7)判断是否达到最大迭代次数：

将当前的分割迭代次数加1作为新的分割迭代次数，判定新的迭代次数是否达到图像分割所设定的分割最大迭代次数20，如果没有达到，返回步骤(3)；否则，执行步骤(8)；

(8)输出分割图像。

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

第一，由于本发明采用判别随机场对用于判断分割类别的联合概率进行建模，克服了现有技术采用三马尔可夫随机场TMRF模型在捕获观测数据的统计相关性时，其作用距离必须依赖于所采用的邻域系统的缺点。使得本发明在考虑图像观测信息的整体性的基础上，充分利用了观测信息，提高了分割结果的准确性和对噪声的鲁棒性。

第二，由于本发明在建立三重判别随机场时引入了图像的非平稳信息，克服了条件随机场未明确地从空域结构的角度上考虑图像的非平稳性的缺点，使得本发明建立的模型与待分割合成孔径雷达SAR图像的模型拟合性得到改善，提高了分割结果的一致性和边缘定位的准确性。

附图说明

图1为本发明的流程图；

图2为本发明的仿真图。

具体实施方式

下面结合附图对本发明做进一步的描述。

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

步骤1，输入待分割SAR图像。

步骤2，提取特征和初始化。

提取特征主要是为了充分利用图像信息。

选取大小为7×7的滑动窗口,计算待分割图像中滑动窗口内灰度值的均值h¹、方差h²、熵h³等图像统计特征，以及半方差特征h⁴,h⁵,h⁶,h⁷，其中半方差特征是由下式提取得到：

其中，v表示待分割的合成孔径雷达SAR图像灰度值的半方差特征；N_d(L)表示在待分割的合成孔径雷达SAR图像中所选取的滑动窗口内，d方向上间隔距离为L的象素对的个数；d表示待分割的合成孔径雷达SAR图像中滑动窗口内，像素对坐标的相对方向；L表示待分割的合成孔径雷达SAR图像中滑动窗口内，像素对之间的距离；Σ表示求和操作；表示待分割的合成孔径雷达SAR图像中滑动窗口内像素对的序号；y_s表示待分割的合成孔径雷达SAR图像中s坐标位置的象素点；s表示像素点在待分割合成孔径雷达SAR图像中的坐标位置；表示待分割的合成孔径雷达SAR图像中坐标位置的象素点；表示在待分割合成孔径雷达SAR图像中所选取滑动窗口内，第个像素对所对应的坐标。

将半径为15的窗口在待分割的合成孔径雷达SAR图像内滑动，采用几何矩估计方法，计算滑窗内待分割合成孔径雷达SAR图像的非平稳各向异性高斯核参数h⁸,h⁹。

采用直方图阈值方法，对待分割的合成孔径雷达SAR图像进行初始分割得到待分割的合成孔径雷达SAR图像的初始分割图像。

采用模糊C均值聚类方法，对待分割合成孔径雷达SAR图像的平稳态进行初始分类得到待分割的合成孔径雷达SAR图像的初始平稳态分类结果。

利用最小二乘法，估计用于判定分割类别的邻域内像素对的相互势能函数的参数。

将用于判定分割类别的联合势能函数的参数全部初始化为1。

将当前的图像分割迭代次数设为1，将图像分割的最大迭代次数设为20。

步骤3，获得联合概率分布。

按照下式，计算用于判定分割类别的联合势能：

$A=\underset{n}{\mathrm{\Σ}}\underset{m}{\mathrm{\Σ}}\left(\mathrm{\δ}(x_s,n)\mathrm{\δ}(u_s,m)\right)\log G$

其中，A表示用于判定分割类别的联合势能；Σ表示求和操作；log表示求对数操作；n表示待分割合成孔径雷达SAR图像的分割标记；m表示待分割合成孔径雷达SAR图像的辅助标记；δ(x_s,n)表示待分割合成孔径雷达SAR图像中坐标位置为s的像素点的分割标记和待分割合成孔径雷达SAR图像的分割标记n的一致性的判断，取值为1时代表分割标记相同，取值为0时代表分割标记不相同；x_s表示待分割合成孔径雷达SAR图像中坐标位置为s的像素点的分割标记；s表示像素点在待分割合成孔径雷达SAR图像中的位置坐标；δ(u_s,m)表示待分割合成孔径雷达SAR图像中坐标位置为s的像素点的平稳态标记和待分割合成孔径雷达SAR图像的平稳态标记m的一致性的判断，取值为1时代表平稳态标记相同，取值为0时代表平稳态标记不相同；u_s表示待分割合成孔径雷达SAR图像中坐标位置为s的像素点的辅助标记；G表示待分割合成孔径雷达SAR图像中坐标位置为s的像素点的分割标记为分割标记n，且待分割合成孔径雷达SAR图像中坐标位置为s的像素点的辅助标记为辅助标记m的概率大小，由下式计算得到：

$G=\left\{\begin{array}{cc}\frac{\exp \left({\mathrm{\ω}}_{\mathrm{nm}}^{T}h\right)}{1+\underset{l=1}{\overset{N-1}{\mathrm{\Σ}}}\underset{k=1}{\overset{M}{\mathrm{\Σ}}}\exp \left({\mathrm{\ω}}_{\mathrm{lk}}^{T}h\right)}& \mathrm{nm}\≠\mathrm{NM}\\ \frac{1}{1+\underset{l=1}{\overset{N-1}{\mathrm{\Σ}}}\underset{k=1}{\overset{M}{\mathrm{\Σ}}}\exp \left({\mathrm{\ω}}_{\mathrm{lk}}^T{h}_s\left(y\right)\right)}& n=N,m=M\end{array}\right.$

其中，G表示待分割合成孔径雷达SAR图像中坐标位置为s的像素点的分割标记为分割标记n，且待分割合成孔径雷达SAR图像中坐标位置为s的像素点的辅助标记为辅助标记m的概率大小；exp表示求指数操作；ω_nm表示在分割标记为n、辅助标记为m时的联合势能的参数；n表示待分割合成孔径雷达SAR图像的分割标记；m表示待分割合成孔径雷达SAR图像的辅助标记；T表示转置操作；h＝[1,h¹,...,h⁹]^T表示待分割合成孔径雷达SAR图像的特征向量；Σ表示求和操作；l表示分割标记；k表示辅助标记；N表示分割标记的总个数；M表示辅助标记的总个数；ω_lk表示在分割标记为l、辅助标记为k时的联合势能的参数。

按照下式，计算用于判定分割类别的邻域内像素对的相互势能：

$I=\underset{s}{\mathrm{\Σ}}\underset{t\∈N_s}{\mathrm{\Σ}}(2\mathrm{\δ}(x_s,x_t)-1)\×(\mathrm{\α}(1-\mathrm{\δ}(u_s,u_t))+\underset{m}{\mathrm{\Σ}}{\mathrm{\α}}_m{\mathrm{\δ}}^{*}(u_s,u_t,m))\mathrm{\μ})$

其中，I表示用于判定分割类别的邻域内像素对的相互势能；s表示像素点在待分割合成孔径雷达SAR图像中的位置坐标；Σ表示求和操作；t表示在待分割合成孔径雷达SAR图像中位置坐标为s的像素点的邻域系统内的像素点；N_s表示待分割合成孔径雷达SAR图像中坐标位置为s的像素点的邻域系统；δ(x_s,x_t)表示待分割合成孔径雷达SAR图像中坐标位置为s的像素点的分割标记和待分割合成孔径雷达SAR图像中坐标位置为t的像素点的分割标记的一致性的判断，取值为1时代表分割标记相同，取值为0时代表分割标记不相同；x_s表示待分割合成孔径雷达SAR图像中坐标位置为s的像素点的分割标记；x_t表示待分割合成孔径雷达SAR图像中坐标位置为t的像素点的分割标记；α表示用于表示辅助信息重要性的相互势能参数；α表示平稳态标记不同时，待分割合成孔径雷达SAR图像中坐标位置为s的像素点和待分割合成孔径雷达SAR图像中坐标位置为t的像素点之间相关性的参数；δ(u_s,u_t)表示待分割合成孔径雷达SAR图像中坐标位置为s的像素点的平稳态标记和待分割合成孔径雷达SAR图像中坐标位置为t的像素点的平稳态标记的一致性的判断，取值为1时代表平稳态标记相同，取值为0时代表平稳态标记不相同；u_s表示待分割合成孔径雷达SAR图像中坐标位置为s的像素点的辅助标记；u_t表示待分割合成孔径雷达SAR图像中坐标位置为t的像素点的辅助标记；m表示待分割合成孔径雷达SAR图像的辅助标记；α_m表示在平稳态标记为m时，待分割合成孔径雷达SAR图像中坐标位置为s的像素点和待分割合成孔径雷达SAR图像中坐标位置为t的像素点的相关性的参数；δ^*(u_s,u_t,m)表示待分割合成孔径雷达SAR图像中坐标位置为s的像素点的平稳态标记、待分割合成孔径雷达SAR图像中坐标位置为t的像素点的平稳态标记和平稳态标记m的一致性的判断，取值为1时代表平稳态标记相同，取值为0时代表平稳态标记不相同；μ表示待分割合成孔径雷达SAR图像中坐标位置为s的像素点的特征和待分割合成孔径雷达SAR图像中坐标位置为t的像素点的特征组成的特征向量。

按照下式，计算用于判定分割类别的似然分布：

$f=\underset{s}{\mathrm{\Π}}\frac{\mathrm{\β}{y_s}^{\mathrm{\β\λ}-1}}{{\mathrm{\α}}^{\mathrm{\β\λ}}\mathrm{\Γ}}\exp (-{\left(\frac{y_s}{\mathrm{\α}}\right)}^{\mathrm{\β}})$

其中，f表示用于判定分割类别的似然分布；Π表示求乘积操作；s表示像素点在待分割合成孔径雷达SAR图像中的位置坐标；β表示广义伽马Gamma分布的形状参数；y_s表示待分割合成孔径雷达SAR图像中坐标为s的像素点的灰度值；λ表示广义伽马Gamma分布的指示形状参数；α表示广义伽马Gamma分布的尺度参数；Γ表示伽马Gamma函数；exp表示求指数操作。

按照下式，计算用于判定分割类别的联合概率分布：

$p=\frac{1}{Z}\exp (\underset{s}{\mathrm{\Σ}}A+\underset{s}{\mathrm{\Σ}}\underset{t\∈N_s}{\mathrm{\Σ}}I+\underset{s}{\mathrm{\Σ}}\log f$

其中，p表示用于判定分割类别的联合概率分布；Z表示归一化因子；exp表示求指数操作；Σ表示求和操作；s表示像素点在待分割合成孔径雷达SAR图像中的位置坐标；A表示用于判定分割类别的联合势能；t表示在待分割合成孔径雷达SAR图像中位置坐标为s的像素点的邻域系统内的像素点；N_s表示待分割合成孔径雷达SAR图像中坐标位置为s的像素点的邻域系统；I表示用于判定分割类别的邻域内像素对的相互势能；f表示用于判定分割类别的似然分布。

步骤4，获得后验边缘概率分布。

对用于判定分割类别的联合概率分布，利用吉布斯Gibbs采样方法，得到用于判定分割类别的后验边缘概率分布。

步骤5，分割SAR图像。

利用贝叶斯最大后验边缘概率准则，确定每个像素点新的标号，逐点更新待分割合成孔径雷达SAR图像的标号，得到分割图像，是按下式计算得到：

$v=\frac{1}{N_d\left(L\right)}\underset{a=1}{\overset{N_d\left(L\right)}{\mathrm{\Σ}}}{(y_s-y_{s_a})}^2$

其中，v表示待分割的合成孔径雷达SAR图像灰度值的半方差特征；N_d(L)表示在待分割的合成孔径雷达SAR图像中所选取滑动窗口内d方向上间隔距离为L的象素对的个数；d表示待分割的合成孔径雷达SAR图像中滑动窗口内像素对坐标的相对方向；L表示待分割的合成孔径雷达SAR图像中滑动窗口内像素对之间的距离；Σ表示求和操作；a表示待分割的合成孔径雷达SAR图像中滑动窗口内像素对的序号；y_s表示待分割的合成孔径雷达SAR图像中坐标为s的象素点的灰度值；s表示像素点在待分割合成孔径雷达SAR图像中的位置坐标；表示待分割的合成孔径雷达SAR图像中坐标为s_a的象素点的灰度值；s_a表示在待分割合成孔径雷达SAR图像中所选取滑动窗口内,与在待分割合成孔径雷达SAR图像中坐标为s的像素点,在d方向上距离为L的点的坐标位置。

步骤6，估计模型参数。

在待分割合成孔径雷达SAR图像的分割图像上，利用对比发散方法，估计用于判定分割类别的联合势能的参数，按如下步骤进行。

第一步，将对比发散方法的初始迭代次数设为1，将对比发散方法的最大迭代次数设为20，将用于判定分割类别的联合势能的参数的更新值设为0。

第二步，按照下式，计算用于获得联合势能参数的更新值的能量：

$E=1+\underset{n}{\mathrm{\Σ}}\underset{m}{\mathrm{\Σ}}\exp \left({\mathrm{\ω}}_{\mathrm{nm}}^{T}h\right)$

其中，E表示用于获得联合势能参数的更新值的能量；Σ表示求和操作；n表示待分割合成孔径雷达SAR图像的分割值；m表示待分割合成孔径雷达SAR图像的平稳态；exp表示求指数操作；ω_nm表示在分割值为n、平稳态为m时的联合势能参数；T表示转置操作；h表示待分割合成孔径雷达SAR图像的特征向量；

第三步，按照下式计算用于判定分割类别的联合势能的参数的更新值：

$\mathrm{\θ}=\mathrm{\τ}\underset{s\∈S}{\mathrm{\Σ}}\left\{\right[\mathrm{\δ}(x_s,n)\mathrm{\δ}(u_s,m)-q\left]h(1-\frac{\exp \left({\mathrm{\ω}}^{T}h\right)}{E})\right)\left]\right\}$

其中，θ表示用于判定分割类别的联合势能参数的更新值；Σ表示求和操作；τ表示调整用于判定分割类别的联合势能参数的更新值的因子；δ(x_s,n)表示一个常数，当分割图像中坐标位置为s的像素点和待分割合成孔径雷达SAR图像的分割值n相同时，δ(x_s,n)＝1，当分割图像中坐标位置为s的像素点和待分割合成孔径雷达SAR图像的分割值n不同时，δ(x_s,n)＝0；x_s表示分割图像s坐标位置的像素点；s表示像素点在待分割合成孔径雷达SAR图像中的坐标位置；n表示待分割合成孔径雷达SAR图像的分割值；δ(u_s,m)表示一个常数，当待分割合成孔径雷达SAR图像中坐标位置为s的像素点的平稳态和待分割合成孔径雷达SAR图像的平稳态m相同时，δ(u_s,m)＝1，当待分割合成孔径雷达SAR图像中坐标位置为s的像素点的平稳态和待分割合成孔径雷达SAR图像的平稳态m不同时，δ(u_s,m)＝0；u_s表示待分割合成孔径雷达SAR图像中s坐标位置的像素点的平稳态；m表示待分割合成孔径雷达SAR图像的平稳态；q表示对用于判定分割类别的联合概率，进行吉布斯Gibbs采样后的分割值和平稳态；h表示待分割合成孔径雷达SAR图像的特征向量；exp表示求指数操作；ω表示用于判定分割类别的联合势能的参数；T表示转置操作；E表示用于获得联合势能参数的更新值的能量；

第四步，按照下式更新用于判定分割类别的联合势能的参数ω：

ω^(γ)＝ω^(γ-1)+ρθ^(γ-1)+(1-ρ)θ^(γ)

其中，ω^(γ)表示第γ次参数估计时，用于判定分割类别的联合势能参数；γ表示当前的对比发散方法的迭代次数；ω^(γ-1)表示第γ-1次参数估计时，用于判定分割类别的联合势能的参数；ρ表示调整用于判定分割类别的联合势能参数的更新值的因子；θ^(γ-1)表示第γ-1次参数估计时，用于判定分割类别的联合势能的参数的更新值；θ^(γ)表示第γ次参数估计时，用于判定分割类别的联合势能的参数的更新值；

第五步，更新迭代次数γ：

将当前的迭代次数加1作为新的迭代次数，判断新的迭代次数是否超过对比发散方法所设定的最大迭代次数20，如果是，则得到用于判定分割类别的联合势能参数；否则，返回第二步。

利用最小二乘法，根据待分割合成孔径雷达SAR图像的分割图像，估计用于判定分割类别的邻域内像素对的相互势能参数。

步骤7，判断是否达到最大迭代次数。

将当前的迭代次数加1作为新的迭代次数，判定新的迭代次数是否达到图像分割所设定的最大迭代次数20，如果没有达到，返回步骤(3)；否则，执行下一步骤。

步骤8，输出分割结果。

下面结合图2对本发明的效果做进一步的描述。

图2为本发明的仿真图。图2(a)、图2(b)和图2(c)分别表示三幅原始待分割的合成孔径雷达SAR图像，图2(d)、图2(e)和图2(f)分别表示三幅原始待分割的合成孔径雷达SAR图像，用已有的三重马尔可夫随机场TMRF方法分割的仿真图，图2(g)、图2(h)和图2(i)分别表示三幅原始待分割的合成孔径雷达SAR图像，用本发明提出的三重判别随机场方法分割的仿真图。由图2可以看出，本发明提出的三重判别随机场方法在分割精度、分割区域一致性以及边界定位上明显优于已有的三重马尔可夫随机场TMRF方法。在本发明的分割结果中，误分割减少，分割区域更为平滑，边界定位更为准确，图像的结构信息保持更为精确。为了进一步证明三重判别随机场方法相较已有的三重马尔可夫随机场TMRF方法，有较好的分割准确性，分别用比值图像的对数归一化似然比|D|以及方差RI_var两个客观评价指标，比较两种方法的分割效果。其比较结果如表1所示：

表1三重判别随机场方法与TMRF方法对比

由表1可以看出本发明在对数归一化似然比D以及方差RI_var2个客观评价指标上均优于已有的三重马尔可夫随机场TMRF方法。这是由于本发明在对非平稳合成孔径雷达SAR图像进行统计建模时，利用合成孔径雷达SAR图像观测值的统计相关特性，来提高模型精度。基于图像观测数据特征构建相互势能函数的方法，使得所描述的空域相互作用关系更为精确，从而在边界定位这一性能上展现出本发明的优势，而且本发明综合考虑了合成孔径雷达SAR图像的统计分布特性以及统计纹理特征，实现了对非平稳SAR图像更为精确的统计建模，从而使得本发明取得了高的分割精度。

Claims (4)

1.一种基于三重判别随机场的SAR图像分割方法，包括如下步骤：

(1)输入待分割SAR图像；

(2)提取特征和初始化：

(2a)将半径为7的窗口在待分割的合成孔径雷达SAR图像内滑动，采用统计方法，计算滑窗内待分割的合成孔径雷达SAR图像像素点灰度值的均值、方差和熵三个图像统计特征，利用半方差方法，提取滑窗内待分割的合成孔径雷达SAR图像的四个半方差特征；

(2b)将半径为15的窗口在待分割的合成孔径雷达SAR图像内滑动，采用几何矩估计方法，计算滑窗内待分割合成孔径雷达SAR图像的非平稳各向异性高斯核参数；

(2c)采用直方图阈值方法，对待分割的合成孔径雷达SAR图像进行初始分割，得到待分割的合成孔径雷达SAR图像的初始分割图像；

(2d)采用模糊C均值聚类方法，对待分割合成孔径雷达SAR图像的平稳态进行初始分类，得到待分割的合成孔径雷达SAR图像的初始平稳态分类；

(2e)利用最小二乘法，根据待分割的合成孔径雷达SAR图像的初始分割图像和待分割的合成孔径雷达SAR图像的初始平稳态分类，估计用于判定分割类别的邻域内像素对的相互势能参数；

(2f)将用于判定分割类别的联合势能函数的参数全部初始化为1；

(2g)将当前的图像分割迭代次数设为1，将图像分割的最大迭代次数设为20；

(3)获得联合概率分布：

(3a)按照下式，构建用于判定分割类别的联合势能：

其中，A表示用于判定分割类别的联合势能；Σ表示求和操作；n表示待分割合成孔径雷达SAR图像的分割值；m表示待分割合成孔径雷达SAR图像的平稳态；δ(x_s,n)表示一个常数，当分割图像中s坐标位置的像素点和待分割合成孔 径雷达SAR图像的分割值n相同时，δ(x_s,n)＝1，当分割图像中s坐标位置的像素点和待分割合成孔径雷达SAR图像的分割值n不同时，δ(x_s,n)＝0；x_s表示分割图像s坐标位置的像素点；s表示像素点在待分割合成孔径雷达SAR图像中的坐标位置；δ(u_s,m)表示一个常数，当待分割合成孔径雷达SAR图像中s坐标位置的像素点的平稳态和待分割合成孔径雷达SAR图像的平稳态m相同时，δ(u_s,m)＝1，当待分割合成孔径雷达SAR图像中s坐标位置的像素点的平稳态和待分割合成孔径雷达SAR图像的平稳态m不同时，δ(u_s,m)＝0；u_s表示待分割合成孔径雷达SAR图像中s坐标位置的像素点的平稳态；log表示求对数操作；G表示待分割图像中s坐标位置的像素点为分割值n，且待分割合成孔径雷达SAR图像中s坐标位置的像素点的平稳态为平稳态m的概率；

(3b)按照下式，构建用于判定分割类别的邻域内像素对的相互势能：

其中，I表示用于判定分割类别的邻域内像素对的相互势能；s表示像素点在待分割合成孔径雷达SAR图像中的坐标位置；Σ表示求和操作；t表示在待分割合成孔径雷达SAR图像中坐标位置为s的像素点的邻域系统内的像素点；N_s表示待分割合成孔径雷达SAR图像中s坐标位置的像素点的邻域系统；δ(x_s,x_t)表示一个常数，当分割图像中s坐标位置的像素点和分割图像中t坐标位置的像素点相同时，δ(x_s,x_t)＝1，当分割图像中s坐标位置的像素点和分割图像中t坐标位置的像素点不同时，δ(x_s,x_t)＝1；x_s表示分割图像s坐标位置的像素点；x_t表示分割图像t坐标位置的像素点；α表示平稳态不同时，待分割合成孔径雷达SAR图像中s坐标位置的像素点和待分割合成孔径雷达SAR图像中t坐标位置的像素点之间相关性的参数；δ(u_s,u_t)表示一个常数，当待分割合成孔径雷达SAR图像中s坐标位置的像素点的平稳态和待分割合成孔径雷达SAR图像中t坐标位置的像素点的平稳态相同时，δ(u_s,u_t)＝1，当待分割合成孔径雷达SAR 图像中s坐标位置的像素点的平稳态和待分割合成孔径雷达SAR图像中t坐标位置的像素点的平稳态不同时，δ(u_s,u_t)＝0；u_s表示待分割合成孔径雷达SAR图像中s坐标位置的像素点的平稳态；u_t表示待分割合成孔径雷达SAR图像中t坐标位置的像素点的平稳态；m表示待分割合成孔径雷达SAR图像的平稳态；α_m表示在平稳态为m时，待分割合成孔径雷达SAR图像中s坐标位置的像素点和待分割合成孔径雷达SAR图像中t坐标位置的像素点的相关性的参数；δ^*(u_s,u_t,m)表示一个常数，当待分割合成孔径雷达SAR图像中s坐标位置的像素点的平稳态、待分割合成孔径雷达SAR图像中t坐标位置的像素点的平稳态和平稳态m均相同时，δ^*(u_s,u_t,m)＝1，当待分割合成孔径雷达SAR图像中s坐标位置的像素点的平稳态、待分割合成孔径雷达SAR图像中t坐标位置的像素点的平稳态和平稳态m不同时，δ^*(u_s,u_t,m)＝0；μ表示待分割合成孔径雷达SAR图像中s坐标位置的像素点的特征和待分割合成孔径雷达SAR图像中t坐标位置的像素点的特征，组成的特征向量；

(3c)按照下式，构建用于判定分割类别的似然分布：

其中，f表示用于判定分割类别的似然分布；Π表示求乘积操作；s表示像素点在待分割合成孔径雷达SAR图像中的坐标位置；β表示广义伽马Gamma分布的形状参数；y_s表示待分割合成孔径雷达SAR图像中s坐标位置的像素点；λ表示广义伽马Gamma分布的指示形状参数；α表示广义伽马Gamma分布的尺度参数；Γ表示伽马Gamma函数；

(3d)按照下式，构建用于判定分割类别的联合概率分布：

其中，p表示用于判定分割类别的联合概率分布；Z表示归一化因子；exp表 示求指数操作；Σ表示求和操作；s表示像素点在待分割合成孔径雷达SAR图像中的坐标位置；A表示用于判定分割类别的联合势能；t表示在待分割合成孔径雷达SAR图像中s坐标位置的像素点的邻域系统内像素点的坐标位置；N_s表示待分割合成孔径雷达SAR图像中s坐标位置的像素点的邻域系统；I表示用于判定分割类别的邻域内像素对的相互势能；f表示用于判定分割类别的似然分布；

(4)获得后验边缘概率分布：

对用于判定分割类别的联合概率分布，利用吉布斯Gibbs采样方法，得到用于判定分割类别的后验边缘概率分布；

(5)分割SAR图像：

利用贝叶斯最大后验边缘概率准则，确定每个像素点新的标号，逐点更新待分割合成孔径雷达SAR图像的标号，得到分割图像；

(6)估计模型参数：

(6a)在待分割合成孔径雷达SAR图像的分割图像上，利用对比发散方法，估计用于判定分割类别的联合势能参数；

(6b)利用最小二乘法，根据待分割合成孔径雷达SAR图像的分割图像，估计用于判定分割类别的邻域内像素对的相互势能参数；

(7)判断是否达到最大迭代次数：

将当前的分割迭代次数加1作为新的分割迭代次数，判定新的迭代次数是否达到图像分割所设定的分割最大迭代次数20，如果没有达到，返回步骤(3)；否则，执行步骤(8)；

(8)输出分割图像。

2.根据权利要求1所述的基于三重判别随机场的SAR图像分割方法，其特征在于，步骤(2a)所述的半方差方法，是按下式计算得到SAR图像灰度值的半方差特征：

其中，v表示待分割的合成孔径雷达SAR图像灰度值的半方差特征；N_d(L) 表示在待分割的合成孔径雷达SAR图像中所选取的滑动窗口内，d方向上间隔距离为L的像素对的个数；d表示待分割的合成孔径雷达SAR图像中滑动窗口内，像素对坐标的相对方向；L表示待分割的合成孔径雷达SAR图像中滑动窗口内，像素对之间的距离；Σ表示求和操作；表示待分割的合成孔径雷达SAR图像中滑动窗口内像素对的序号；y_s表示待分割的合成孔径雷达SAR图像中s坐标位置的像素点；s表示像素点在待分割合成孔径雷达SAR图像中的坐标位置；表示待分割的合成孔径雷达SAR图像中坐标位置的像素点；表示在待分割合成孔径雷达SAR图像中所选取滑动窗口内，第个像素对所对应的坐标。

3.根据权利要求1所述的基于三重判别随机场的SAR图像分割方法，其特征在于，步骤(5)所述的贝叶斯最大后验边缘概率准则公式如下：

其中，x_s表示分割图像s坐标位置的像素点；s表示像素点在待分割合成孔径雷达SAR图像中的坐标位置；arg max表示求最大值操作；u_s表示待分割合成孔径雷达SAR图像中s坐标位置的像素点的平稳态；Σ表示求和操作；r_s表示合成孔径雷达SAR图像中s坐标位置的像素点的用于判定分割类别的后验边缘概率。

4.根据权利要求1所述基于三重判别随机场的SAR图像分割方法，其特征在于，步骤(6a)所述的对比发散方法是按照下步骤获得用于判定分割类别的联合势能参数：

第一步，将对比发散方法的初始迭代次数设为1，将对比发散方法的最大迭代次数设为20，将用于判定分割类别的联合势能参数的更新值设为0；

第二步，按照下式，计算用于获得联合势能参数的更新值的能量：

其中，E表示用于获得联合势能参数的更新值的能量；Σ表示求和操作；n表示待分割合成孔径雷达SAR图像的分割值；m表示待分割合成孔径雷达SAR图像的平稳态；exp表示求指数操作；ω_nm表示在分割值为n、平稳态为m时的联合势能参数；T表示转置操作；h表示待分割合成孔径雷达SAR图像的特征向量；

第三步，按照下式计算用于判定分割类别的联合势能参数的更新值：

其中，θ表示用于判定分割类别的联合势能参数的更新值；Σ表示求和操作；τ表示调整用于判定分割类别的联合势能参数的更新值的因子；δ(x_s,n)表示一个常数，当分割图像中坐标位置为s的像素点和待分割合成孔径雷达SAR图像的分割值n相同时，δ(x_s,n)＝1，当分割图像中坐标位置为s的像素点和待分割合成孔径雷达SAR图像的分割值n不同时，δ(x_s,n)＝0；x_s表示分割图像s坐标位置的像素点；s表示像素点在待分割合成孔径雷达SAR图像中的坐标位置；n表示待分割合成孔径雷达SAR图像的分割值；δ(u_s,m)表示一个常数，当待分割合成孔径雷达SAR图像中坐标位置为s的像素点的平稳态和待分割合成孔径雷达SAR图像的平稳态m相同时，δ(u_s,m)＝1，当待分割合成孔径雷达SAR图像中坐标位置为s的像素点的平稳态和待分割合成孔径雷达SAR图像的平稳态m不同时，δ(u_s,m)＝0；u_s表示待分割合成孔径雷达SAR图像中s坐标位置的像素点的平稳态；m表示待分割合成孔径雷达SAR图像的平稳态；q表示对用于判定分割类别的联合概率，进行吉布斯Gibbs采样后的分割值和平稳态；h表示待分割合成孔径雷达SAR图像的特征向量；exp表示求指数操作；ω表示用于判定分割类别的联合势能的参数；T表示转置操作；E表示用于获得联合势能参数的更新值的能量；

第四步，按照下式更新用于判定分割类别的联合势能的参数ω：

ω^(γ)＝ω^(γ-1)+ρθ^(γ-1)+(1-ρ)θ^(γ)

其中，ω^(γ)表示第γ次参数估计时，用于判定分割类别的联合势能参数；γ表示当前的对比发散方法的迭代次数；ω^(γ-1)表示第γ-1次参数估计时，用于判定分割类别的联合势能的参数；ρ表示调整用于判定分割类别的联合势能参数的更新值的因子；θ^(γ-1)表示第γ-1次参数估计时，用于判定分割类别的联合势能的参数的更新值；θ^(γ)表示第γ次参数估计时，用于判定分割类别的联合势能的参数的更新值；

第五步，更新迭代次数γ：

将当前的迭代次数加1作为新的迭代次数，判断新的迭代次数是否超过对比发散方法所设定的最大迭代次数20，如果是，则得到用于判定分割类别的联合势能参数；否则，返回第二步。