CN105513085A - 采用srgb-rmrf的sar海冰图像分割方法 - Google Patents

Abstract

本发明公开了一种针对SAR海冰图像的分割方法，主要解决SAR海冰图像受相干斑噪声影响严重，分割结果不能真实反映地表海冰分布情况的问题。其实现过程是：（1）采用双边滤波器对SAR海冰图像进行滤波；（2）采用改进的边缘权重计算模型和区域内部差异进行初始化分割；（3）基于初始化分割结构构件区域相邻图；（4）构建区域MRF模型，并在MRF中融入区域间强度差异惩罚函数；（5）擦除边缘像素点，得到最终分割结果。本发明显著提升了SAR海冰图像的分割性能。

Description

采用SRGB-RMRF的SAR海冰图像分割方法

技术领域

本发明涉及图像分割方法领域，具体是一种采用SRGB-RMRF的SAR海冰图像分割方法。

背景技术

在海冰覆盖的海域，海冰信息对于船舶的航线设计十分重要，随着全球气候的逐渐变暖使得极地冰川融化，海冰信息在气候研究邻域也具有重要的意义。合成孔径雷达(SyntheticApertureRadar,SAR)作为一种高分辨率微波遥感成像雷达，所具有的全天候、全天时、多波段、多极化、可变视角、穿透性强等特点，使得其被广泛的应用于海冰监测上。SAR海冰图像分割作为该领域的关键技术之一，越来越受到人们的关注。

经典的SAR海冰图像分割方法主要包括：阈值方法、结合边沿检测和区域合并混合方法、聚类方法、基于马尔可夫随机场(Markovrandomfield,MRF)的方法等。其中，基于MRF的分割方法根据贝叶斯理论，将MRF空间上下文模型与特征模型相结合，形成MAP-MRF(maximumaposterionMRF)框架，用更为合理的数学方式解决图像分割的难题。传统的基于像素的MRF模型以像素为单元进行分割，分割准确率低，并且，像素级的处理计算量大。因此，基于区域的MRF模型被相继提出并很好的应用到了SAR图像分割中。Deng和Clausi[6]通过在MRF模型和特征模型之间引入自适应权重，提高了算法的分割性能，但先验模型仅仅考虑了邻域标号的平滑性约束，忽视了区域间的差异信息。Yu和Clausi[7]以及Dawoud等[8]将边界强度融入到MLL模型中，使算法更适用于SAR图像分割。但是他们采用高斯模型对SAR图像的特征建模是不合理的，因为根据相干斑噪声模型，SAR图像服从的是Gamma分布，在此基础上利用边界强度进行分割并不能有效提升算法的分割性能。以上几种算法均利用分水岭方法对图像进行初始分割，该方法分割结果好坏完全依赖于图像的梯度估计。由于SAR独特的成像系统和机理，使得SAR图像中存在严重相干斑噪声和特征差异，分水岭算法得到的初始分割结果存在较大的区域数目，且在低对比度的边缘处，由于梯度较小往往导致欠分割现象。综上所述，如何提高初始分割质量，并提出更为合理的区域MRF模型，对SAR海冰图像分割有着重要的意义。

发明内容

本发明的目的是提供一种采用SRGB-RMRF的SAR海冰图像分割方法，以解决SAR海冰图像受相干斑噪声影响严重的问题，以相干斑抑制能力更强的SRGB方法为后续区域MRF提供更加准确的初始分割区域，并在MRF模型中引入区域间强度差异惩罚函数，以得到更为准确的分割结果。

为了达到上述目的，本发明所采用的技术方案为：

采用SRGB-RMRF的SAR海冰图像分割方法，其特征在于：包括以下步骤：

(1)、用SRGB方法对SAR海冰图像进行分割：

1a)对图像进行双边滤波；

给定输入图像y，则经过双边滤波输出图像h为：

$\left\{\begin{array}{c}h\left(x\right)=\frac{1}{Z}\∫y\left(\mathrm{\ϵ}\right)c\left(\mathrm{\ϵ},x\right)g\left(y\left(\mathrm{\ϵ}\right),y\left(x\right)\right)d\mathrm{\ϵ}\\ Z=\∫c\left(\mathrm{\ϵ},x\right)g\left(y\left(\mathrm{\ϵ}\right),y\left(x\right)\right)d\mathrm{\ϵ}\end{array}\right.,$

其中，Z为归一化系数，x表示中心像素，ε为中心像素邻域中的像素点，y(ε)为ε点的像素值，σ_d和σ_r分别控制空间与灰度相似度的衰减程度，c(ε,x)＝exp(-(||ε-x||/σ_d)²/2)为像素的空间相似度权重，灰度相似度权重为g(y(ε),y(x))＝exp(-(||y(ε)-y(x)||/σ_r)²/2)；

1b)计算各点的边缘权重w(e)，按升序对边排列得到e₁,e₂,e₃......e_N；

$w\left(e\right)=\left|\frac{I_i-I_j}{I_i+I_j}\right|,$

其中，I_i，I_j为对应的强度值；

1c)对最小边缘权重的边e_n进行合并判断，假设其所连接区域为C₁和C₂，若C₁和C₂不属于同一区域且Dif(C₁,C₂)≤MInt(C₁,C₂)，则合并两个区域并继续执行第1d)步。

$Dif(C_1,C_2)\underset{v_i\∈C_1,v_j\∈C_2}{\min }w\left(\right(v_i,v_j\left)\right),$

其中，v_i和v_j为相邻像素点，Dif(C₁,C₂)为区域C₁和C₂的差异值，用区域间最小边缘权重来定义；

int(C)＝1/N*∑_e∈MST(C,E)w(e)，

Int(C)＝int(C)+τ(C)，

MInt(C₁,C₂)＝min(Int(C₁),Int(C₂))，

将区域C的内部差异值int(C)定义为C中边缘权重的均值，τ(C)＝k/|C|，k为可调节参数，k值越大，界定的可以区分两个区域的界限就越明显，MInt(C₁,C₂)为Int(C₁)和Int(C₂)中的最小值；

1d)更新int和区域标号，如果n≤N，则按照顺序，选择下一条边返回1c)执行，否则输出初始分割结果；

(2)、进行结合区域间强度差异的区域MRF分割；

2a)用EM算法估计特征模型参数；

EM算法通过迭代E步骤和M步骤估算得到各类的均值、方差，其中E步骤如下：

$p(y_s/u_i,{\mathrm{\σ}}_i)=\frac{1}{{\mathrm{\σ}}_i\sqrt{2\mathrm{\π}}}\exp (-\frac{{(y_s-u_i)}^2}{2{{\mathrm{\σ}}_i}^2}),$

$w_{si}=\frac{p(y_s/u_i,{\mathrm{\σ}}_i){\mathrm{\π}}_i}{{\mathrm{\Σ}}_{i=1}^{n}p(y_s/u_i,{\mathrm{\σ}}_i){\mathrm{\π}}_i},$

M步骤如下：

$u_i=\frac{{\mathrm{\Σ}}_{s\∈S}{w}_{si}{y}_s}{{\mathrm{\Σ}}_{s\∈S}{w}_{si}},$

${\mathrm{\σ}}_i^2=\frac{{\mathrm{\Σ}}_{s\∈S}{w}_{si}{(y_s-u_i)}^2}{{\mathrm{\Σ}}_{s\∈S}{w}_{si}},$

${\mathrm{\π}}_i=\frac{{\mathrm{\Σ}}_{s\∈S}{w}_{si}}{{\mathrm{\Σ}}_i^n{\mathrm{\Σ}}_{s\∈S}{w}_{si}},$

其中，i代表类别，n表示类别总数，p(y_s/u_i,σ_i)表示区域s属于类i的概率，y_s为区域s的强度值，u_i为类i的强度均值，σ_i为类i的强度方差，w_si为区域s属于i类的概率，π_i为类i的先验概率；

2b)计算各相邻区域间的强度差异，构造惩罚函数D；

$D\left(d_{RT}\right)=e^{-{(d_{RT}/K)}^2},$

$d_{RT}=\left|\frac{u_R-u_T}{u_R+u_T}\right|,$

K(t+1)＝1.1.K(t)，

其中，d_RT为区域间强度差异，区域R和T的强度均值是u_R和u_T，参数K控制惩罚衰减速度，初始K(0)值为0.01，K(t+1)＝1.1.K(t)，随着迭代次数t的增加，K值增大，不同相邻区域之间的惩罚差异减小；

2c)计算相同标号的相邻区域合并前后的能量差值δE_ij，并找出最小值δE_min，合并δE_min为负值的相邻区域；

$\begin{array}{c}{\mathrm{\&delta;E}}_{ij}=E_h-E_i-E_j\\ =l\left(\underset{s\&Element;R_h}{\mathrm{\&Sigma;}}{\mathrm{ln\&mu;}}_{L_h}-\underset{s\&Element;R_i}{\mathrm{\&Sigma;}}{\mathrm{ln\&mu;}}_{L_i}-\underset{s\&Element;R_j}{\mathrm{\&Sigma;}}{\mathrm{ln\&mu;}}_{L_j}\right)-{\mathrm{\&beta;}}_0\underset{s\&Element;R_i,t\&Element;R_j}{\underset{<s,t>\&Element;c}{\mathrm{\&Sigma;}}}D\left(d_{st}\right)\end{array}$

其中，R_i和R_j为相邻区域，R_h为合并后的新区域，E_h、E_i、E_j分别为区域h、i、j的能量函数值，若δE_ij值小于0，说明R_i和R_j的合并有利于能量最小化；

2d)判断迭代次数，若已达到最大次数，输出分割结果；若未达到，增大系数K，返回2a)继续迭代。

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

1)本发明针对经典的Graph-Based方法做出改进，得到分割新方法SRGB，使其在自身优势的基础上，具有较强的相干斑抑制能力，实现边缘的准确定位，并减少过分割现象。

2)本发明提出的SRGB算法根据相干斑噪声特性，采用双边滤波器和归一化强度差异来度量像素间的边缘权重，并在此基础上将梯度公式扩展到曼哈顿距离内计算，进一步提升算法的相干斑抑制能力。

3)本发明构建了更为合理的区域MRF模型，采用Gamma分布对MRF模型建模，并将区域间强度差异作为惩罚函数融入MRF模型中，有效的减少误分割现象。

附图说明

图1是本发明的采用SRGB-RMRF的SAR海冰图像分割方法流程图。

图2是本发明中分割方法得到的合成SAR海冰图像分割结果，其中：图2a为区域级MRF(RMRF)算法分割结果图，图2b为基于语义迭代区域生长(IRGS)算法分割结果图，图2c为基于Graph-Based方法的区域级MRF算法分割结果图，图2d为本发明方法分割结果图。

具体实施方式

如图1所示，采用SRGB-RMRF的SAR海冰图像分割方法，包括以下步骤：

(1)、用SRGB方法对SAR海冰图像进行分割：

1a)对图像进行双边滤波；

给定输入图像y，则经过双边滤波输出图像h为：

$\left\{\begin{array}{c}h\left(x\right)=\frac{1}{Z}\&Integral;y\left(\mathrm{\&epsiv;}\right)c\left(\mathrm{\&epsiv;},x\right)g\left(y\left(\mathrm{\&epsiv;}\right),y\left(x\right)\right)d\mathrm{\&epsiv;}\\ Z=\&Integral;c\left(\mathrm{\&epsiv;},x\right)g\left(y\left(\mathrm{\&epsiv;}\right),y\left(x\right)\right)d\mathrm{\&epsiv;}\end{array}\right.,$

其中，Z为归一化系数，x表示中心像素，ε为中心像素邻域中的像素点，y(ε)为ε点的像素值，σ_d和σ_r分别控制空间与灰度相似度的衰减程度,c(ε,x)＝exp(-(||ε-x||/σ_d)²/2)为像素的空间相似度权重，灰度相似度权重为g(y(ε),y(x))＝exp(-(||y(ε)-y(x)||/σ_r)²/2)；

1b)计算各点的边缘权重w(e)，按升序对边排列得到e₁,e₂,e₃......e_N；

$w\left(e\right)=\left|\frac{I_i-I_j}{I_i+I_j}\right|,$

其中，I_i，I_j为对应的强度值；

1c)对最小边缘权重的边e_n进行合并判断，假设其所连接区域为C₁和C₂，若C₁和C₂不属于同一区域且Dif(C₁,C₂)≤MInt(C₁,C₂)，则合并两个区域并继续执行第1d)步。

$Dif(C_1,C_2)\underset{v_i\&Element;C_1,v_j\&Element;C_2}{\min }w\left(\right(v_i,v_j\left)\right),$

其中，v_i和v_j为相邻像素点，Dif(C₁,C₂)为区域C₁和C₂的差异值，用区域间最小边缘权重来定义；

int(C)＝1/N*∑_e∈MST(C,E)w(e)，

Int(C)＝int(C)+τ(C)，

MInt(C₁,C₂)＝min(Int(C₁),Int(C₂))，

将区域C的内部差异值int(C)定义为C中边缘权重的均值，τ(C)＝k/|C|，k为可调节参数，k值越大，界定的可以区分两个区域的界限就越明显，MInt(C₁,C₂)为Int(C₁)和Int(C₂)中的最小值；

1d)更新int和区域标号，如果n≤N，则按照顺序，选择下一条边返回1c)执行，否则输出初始分割结果；

(2)、进行结合区域间强度差异的区域MRF分割；

2a)用EM算法估计特征模型参数；

EM算法通过迭代E步骤和M步骤估算得到各类的均值、方差，其中E步骤如下：

$p(y_s/u_i,{\mathrm{\&sigma;}}_i)=\frac{1}{{\mathrm{\&sigma;}}_i\sqrt{2\mathrm{\&pi;}}}\exp (-\frac{{(y_s-u_i)}^2}{2{{\mathrm{\&sigma;}}_i}^2}),$

$w_{si}=\frac{p(y_s/u_i,{\mathrm{\&sigma;}}_i){\mathrm{\&pi;}}_i}{{\mathrm{\&Sigma;}}_{i=1}^{n}p(y_s/u_i,{\mathrm{\&sigma;}}_i){\mathrm{\&pi;}}_i},$

M步骤如下：

$u_i=\frac{{\mathrm{\&Sigma;}}_{s\&Element;S}{w}_{si}{y}_s}{{\mathrm{\&Sigma;}}_{s\&Element;S}{w}_{si}},$

${\mathrm{\&sigma;}}_i^2=\frac{{\mathrm{\&Sigma;}}_{s\&Element;S}{w}_{si}{(y_s-u_i)}^2}{{\mathrm{\&Sigma;}}_{s\&Element;S}{w}_{si}},$

${\mathrm{\&pi;}}_i=\frac{{\mathrm{\&Sigma;}}_{s\&Element;S}{w}_{si}}{{\mathrm{\&Sigma;}}_i^n{\mathrm{\&Sigma;}}_{s\&Element;S}{w}_{si}},$

其中，i代表类别，n表示类别总数，p(y_s/u_i,σ_i)表示区域s属于类i的概率，y_s为区域s的强度值，u_i为类i的强度均值，σ_i为类i的强度方差，w_si为区域s属于i类的概率，π_i为类i的先验概率；

2b)计算各相邻区域间的强度差异，构造惩罚函数D；

$D\left(d_{RT}\right)=e^{-{(d_{RT}/K)}^2},$

$d_{RT}=\left|\frac{u_R-u_T}{u_R+u_T}\right|,$

K(t+1)＝1.1.K(t)，

其中，d_RT为区域间强度差异，区域R和T的强度均值是u_R和u_T，参数K控制惩罚衰减速度，初始K(0)值为0.01，K(t+1)＝1.1.K(t)，随着迭代次数t的增加，K值增大，不同相邻区域之间的惩罚差异减小；

2c)计算相同标号的相邻区域合并前后的能量差值δE_ij，并找出最小值δE_min，合并δE_min为负值的相邻区域；

$\begin{array}{c}{\mathrm{\&delta;E}}_{ij}=E_h-E_i-E_j\\ =l\left(\underset{s\&Element;R_h}{\mathrm{\&Sigma;}}{\mathrm{ln\&mu;}}_{L_h}-\underset{s\&Element;R_i}{\mathrm{\&Sigma;}}{\mathrm{ln\&mu;}}_{L_i}-\underset{s\&Element;R_j}{\mathrm{\&Sigma;}}{\mathrm{ln\&mu;}}_{L_j}\right)-{\mathrm{\&beta;}}_0\underset{s\&Element;R_i,t\&Element;R_j}{\underset{<s,t>\&Element;c}{\mathrm{\&Sigma;}}}D\left(d_{st}\right)\end{array}$

其中，R_i和R_j为相邻区域，R_h为合并后的新区域，E_h、E_i、E_j分别为区域h、i、j的能量函数值，若δE_ij值小于0，说明R_i和R_j的合并有利于能量最小化；

2d)判断迭代次数，若已达到最大次数，输出分割结果；若未达到，增大系数K，返回2a)继续迭代。

至此，SRGB-RMRF的的SAR海冰图像分割步骤基本完成。以下通过合成SAR海冰图像分割结果对本发明的有效性进行验证。

1.测试数据：

为无噪的合成SAR海冰图像添加视数L为1、2、3、4、6、8、10、16的相干斑噪声进行测试，视数越大，代表相干斑噪声程度越小。

2.对比方法：

将本发明中提及分割方法与区域级MRF(RMRF)分割算法、基于语义迭代区域生长(IRGS)分割算法以及基于Graph-Based方法的区域级MRF算法(GB-RMRF)进行比较。

3.评价指标：

衡量分割算法的准确性一般使用分割准确率和Kappa系数。分割准确率是指像素被正确标记的概率，Kappa系数用来评价分割图像的精度问题，主要用于精确性评价和图像的一致性判断，其公式为：

$Kappa=\frac{P_0-P_c}{1-P_c},$

式中，P₀为整体分类精度，P_c为期望分类精度，Kappa值越大，说明分割效果越好。

4.测试结果(分割准确率/Kappa系数)如表1所示：

表1分割结果比较(整体准确率％/kappa系数)

IRGS算法在视数为1时，分割性能相对较差；RMRF算法在视数小于3时，分割性能劣于GB-RMRF算法；相比而言，本发明算法优于其他三种算法，分割精度和Kappa系数得到大幅度提高，噪声越大，分割优势越明显。由于SAR独特的成像机理，真实的SAR海冰图像含有较强的相干斑噪声，故本发明算法在SAR海冰图像分割中更具实用性。为进一步说明本文算法的有效性，图2给出L＝2时各算法的分割结果图。本发明算法分割结果明显优于其他3种算法，在高噪声情况下，基本完成海冰和水的准确分割，很好的保护了海冰的边界。

Claims (1)

1.采用SRGB-RMRF的SAR海冰图像分割方法，其特征在于：包括以下步骤：

(1)、用SRGB方法对SAR海冰图像进行分割：

1a)对图像进行双边滤波；

给定输入图像y，则经过双边滤波输出图像h为：

$\left\{\begin{array}{c}h\left(x\right)=\frac{1}{Z}\&Integral;y\left(\mathrm{\&epsiv;}\right)c(\mathrm{\&epsiv;},x)g\left(y\right(\mathrm{\&epsiv;}),y(x\left)\right)d\mathrm{\&epsiv;}\\ Z=\&Integral;c(\mathrm{\&epsiv;},x)g\left(y\right(\mathrm{\&epsiv;}),y(x\left)\right)d\mathrm{\&epsiv;}\end{array}\right.,$

其中，Z为归一化系数，x表示中心像素，ε为中心像素邻域中的像素点，y(ε)为ε点的像素值，σ_d和σ_r分别控制空间与灰度相似度的衰减程度，c(ε,x)＝exp(-(||ε-x||/σ_d)²/2)为像素的空间相似度权重，灰度相似度权重为g(y(ε),y(x))＝exp(-(||y(ε)-y(x)||/σ_r)²/2)；

1b)计算各点的边缘权重w(e)，按升序对边排列得到e₁,e₂,e₃......e_N；

$w\left(e\right)=\left|\frac{I_i-I_j}{I_i+I_j}\right|,$

其中，I_i，I_j为对应的强度值；

1c)对最小边缘权重的边e_n进行合并判断，假设其所连接区域为C₁和C₂，若C₁和C₂不属于同一区域且Dif(C₁,C₂)≤MInt(C₁,C₂)，则合并两个区域并继续执行第1d)步。

$Dif(C_1,C_2)=\underset{v_i\&Element;C_1,v_j\&Element;C_2}{\min }w\left(\right(v_i,v_j\left)\right),$

其中，v_i和v_j为相邻像素点，Dif(C₁,C₂)为区域C₁和C₂的差异值，用区域间最小边缘权重来定义；

int(C)＝1/N*∑_e∈MST(C,E)w(e)，

Int(C)＝int(C)+τ(C)，

MInt(C₁,C₂)＝min(Int(C₁),Int(C₂))，

将区域C的内部差异值int(C)定义为C中边缘权重的均值，τ(C)＝k/|C|，k为可调节参数，k值越大，界定的可以区分两个区域的界限就越明显，MInt(C₁,C₂)为Int(C₁)和Int(C₂)中的最小值；

1d)更新int和区域标号，如果n≤N，则按照顺序，选择下一条边返回1c)执行，否则输出初始分割结果；

(2)、进行结合区域间强度差异的区域MRF分割；

2a)用EM算法估计特征模型参数；

EM算法通过迭代E步骤和M步骤估算得到各类的均值、方差，其中E步骤如下：

$p(y_s/u_i,{\mathrm{\&sigma;}}_i)=\frac{1}{{\mathrm{\&sigma;}}_i\sqrt{2\mathrm{\&pi;}}}\exp (-\frac{{(y_s-u_i)}^2}{2{{\mathrm{\&sigma;}}_i}^2}),$

$w_{si}=\frac{p(y_s/u_i,{\mathrm{\&sigma;}}_i){\mathrm{\&pi;}}_i}{{\mathrm{\&Sigma;}}_{i=1}^{n}p(y_s/u_i,{\mathrm{\&sigma;}}_i){\mathrm{\&pi;}}_i},$

M步骤如下：

$u_i=\frac{{\mathrm{\&Sigma;}}_{s\&Element;S}{w}_{si}{y}_s}{{\mathrm{\&Sigma;}}_{s\&Element;S}{w}_{si}},$

${\mathrm{\&sigma;}}_i^2=\frac{{\mathrm{\&Sigma;}}_{s\&Element;S}{w}_{si}{(y_s-u_i)}^2}{{\mathrm{\&Sigma;}}_{s\&Element;S}{w}_{si}},$

${\mathrm{\&pi;}}_i=\frac{{\mathrm{\&Sigma;}}_{s\&Element;S}{w}_{si}}{{\mathrm{\&Sigma;}}_i^n{\mathrm{\&Sigma;}}_{s\&Element;S}{w}_{si}},$

其中，i代表类别，n表示类别总数，p(y_s/u_i,σ_i)表示区域s属于类i的概率，y_s为区域s的强度值，u_i为类i的强度均值，σ_i为类i的强度方差，w_si为区域s属于i类的概率，π_i为类i的先验概率；

2b)计算各相邻区域间的强度差异，构造惩罚函数D；

$D\left(d_{RT}\right)=e^{-{(d_{RT}/K)}^2},$

$d_{RT}=\left|\frac{u_R-u_T}{u_R+u_T}\right|,$

K(t+1)＝1.1.K(t)，

其中，d_RT为区域间强度差异，区域R和T的强度均值是u_R和u_T，参数K控制惩罚衰减速度，初始K(0)值为0.01，K(t+1)＝1.1.K(t)，随着迭代次数t的增加，K值增大，不同相邻区域之间的惩罚差异减小；

2c)计算相同标号的相邻区域合并前后的能量差值δE_ij，并找出最小值δE_min，合并δE_min为负值的相邻区域；

$\begin{array}{c}{\mathrm{\&delta;E}}_{ij}=E_h-E_i-E_j\\ =l(\underset{s\&Element;R_h}{\&Sigma;}{\mathrm{ln\&mu;}}_{L_h}-\underset{s\&Element;R_i}{\&Sigma;}{\mathrm{ln\&mu;}}_{L_i}-\underset{s\&Element;R_j}{\&Sigma;}{\mathrm{ln\&mu;}}_{L_j})-{\mathrm{\&beta;}}_0\underset{s\&Element;R_i,t\&Element;R_j}{\underset{<s,t>\&Element;c}{\mathrm{\&Sigma;}}}D\left(d_{st}\right)\end{array}$

其中，R_i和R_j为相邻区域，R_h为合并后的新区域，E_h、E_i、E_j分别为区域h、i、j的能量函数值，若δE_ij值小于0，说明R_i和R_j的合并有利于能量最小化；

2d)判断迭代次数，若已达到最大次数，输出分割结果；若未达到，增大系数K，返回2a)继续迭代。