CN110136143A - 基于admm算法的马氏场下多分辨率遥感图像分割方法 - Google Patents

Abstract

本发明提出了一种基于ADMM算法的马氏场下多分辨率遥感图像分割方法，其步骤如下：对读取的像素级图像进行初始化过分割，得到对象级过分割区域邻接图、对象级图像的邻域系统、观测特征场和分割标记场；利用混合高斯模型对观测特征场进行概率建模得到特征场模型，利用Gibbs模型对分割标记场进行概率建模得到标记场模型，根据Bayes准则联合特征场模型和标记场模型得到分割标记场的后验分布模型；根据最大后验概率准则运用ADMM算法进行迭代更新求解，获取观测特征场的标记结果作为遥感图像的分割结果。本发明将特征场和标记场能量分开考虑，可更好地求解遥感影像分割标记，分割效率远远高于传统手工分割水平，提高了影像分割精度。

Description

基于ADMM算法的马氏场下多分辨率遥感图像分割方法

技术领域

本发明涉及图像分割的技术领域，尤其涉及一种基于ADMM算法的马氏场下多分辨率遥感图像分割方法。

背景技术

随着卫星技术的发展，卫星的空间分辨率越来越高，对地观测数据量也越来越大，但是我国对于遥感数据的处理速度，却不能够满足对于发展高分辨率遥感系统的迫切需求。因此，针对高分辨率图像的数据处理，挖掘更多遥感影像信息，成为当前研究的重要科学问题。在对高分辨率图像的数据处理分析过程中，图像分割是一个较为重要的基础步骤。

关于遥感影像的分割方法，其中阈值法是最简单的图像分割算法，一般单波段影像是由不同灰度级的像元所组成的，这一点可以从图像的直方图上明确看出，如图2所示，不同的灰度级都有对应的像元数量，选择合理的阈值就可以将图像分割为若干对象过分割区域。基于边缘检测的分割方法通常是寻找影像中灰度值突变的地方来确定局部过分割区域的边缘。基于过分割区域的分割方法主要是依据像素之间的相似性来形成局部过分割区域，从而获取分割结果。对于遥感图像的分割来说，每幅遥感图像可以看作是由随机变量组成的，并且每幅遥感图像所包含的数据量是极大的，因此基于概率统计的建模方法备受关注，其中，又以马尔科夫随机场方法(Markov Random Field，MRF)应用较为广泛。

在一幅数字图像中，大部分像素点仅与其周围数个相邻像素点有相互关系，这种先验信息对于图像分析有着很大帮助。而Markov随机场理论是分析空间依赖性的重要概率统计手段，它能够对上述先验信息进行很好的概率建模。因此，在图像分割中，Markov随机场方法可以通过标记之间的概率建模来有效地描述图像的空间信息。在MRF方法中，对观测场建模得到似然函数，对标记场建模得到先验函数，利用MAP准则对模型进行求解，常用的求解方法有生成式模型、判别式模型、过分割区域生长法和对偶法等。生成式模型是对后验概率建模，从统计的角度表示数据的分布情况，能够反映同类数据本身相似度的一种模型，如经典的像素级条件迭代法ICM。判别式模型是一种能反映一类数据之间差异的方法，如条件随机场方法，CRF。过分割区域生长法是将成组的像素或过分割区域发展成更大过分割区域的过程，从种子点的集合开始，从种子点的过分割区域增长是通过将与每个种子点有相似属性像强度、灰度级、纹理颜色等的相邻像素合并到此过分割区域。如IGRS，一种对象级的MRF算法。以上方法其中一个共同点就是将似然函数和先验函数一起考虑，同时综合似然信息和先验信息，这样做可能会模糊特征场或标记场中的重要信息，造成误分现象。

发明内容

针对现有MRF方法对图像分割时，会模糊特征场或标记场中的重要信息，造成误分现象的技术问题，本发明提出一种基于ADMM算法的马氏场下多分辨率遥感图像分割方法，将似然函数和先验函数分开考虑，采用ADMM算法来对模型进行求解，提高影像分割精度。

为了达到上述目的，本发明的技术方案是这样实现的：一种基于ADMM算法的马氏场下多分辨率遥感图像分割方法，其步骤如下：

步骤一：输入并读取一幅遥感图像，对读取的像素级图像进行初始化过分割，得到由过分割区域组成的对象级图像和对象级过分割区域邻接图，根据过分割区域邻接图定义对象级图像的邻域系统、观测特征场和分割标记场；

步骤二：利用混合高斯模型对观测特征场进行概率建模得到特征场模型，利用Gibbs模型对分割标记场进行概率建模得到标记场模型，根据Bayes准则联合特征场模型和标记场模型得到分割标记场的后验分布模型；

步骤三：根据最大后验概率准则运用ADMM算法对步骤二得到的分割标记场的后验分布模型进行迭代更新求解，获取观测特征场和分割标记场的标记结果，将观测特征场的最终标记结果作为遥感图像的分割结果。

所述步骤一中运用mean-shift方法对像素级图像进行过分割处理并获得对象级过分割区域邻接图的方法为：

将分辨率为m×n的像素级图像I(R,G,B)过分割为l个最小面积为S_min的过分割区域，每个过分割区域赋予标号得到标号矩阵L_s＝{l_s|s∈S}，其中，元素l_s∈{1,2,…,l},s∈S；位置索引集合S＝{s＝(x,y)|1≤x≤m,1≤y≤n}，s表示像素级图像的位置索引，位置索引s的数量有m×n个，(x,y)为像素级图像中像素点的位置坐标，m为像素级图像的长度，n为像素级图像的宽度；且像素级的观测特征集Y⁰＝{y_s＝(r_s,g_s,b_s)}，其中，y_s表示位置索引s处像素点的观测特征向量，r_s、g_s、b_s分别为位置索引s处像素级图像的R、G、B分量的值；

得到对象级图像的位置索引集合R＝{r₁,r₂,…,r_l}，其中，过分割区域r_i＝{s|l_s＝i,i＝1,2,…,l}，i表示过分割区域的标号值，i＝1,2,…,l；l表示过分割区域个数；

定义过分割区域邻接图G＝(V,E)，其中，V＝{V₁,V₂,…,V_l}是位置集合，每个位置V_i代表对应的过分割区域r_i，E＝{e_ij|1≤i,j≤l}表示邻接关系，元素e_ij表示过分割区域r_i中与过分割区域r_j相邻的像素个数，e_ij≠0当且仅当过分割区域r_i和过分割区域r_j是相邻的。

所述对象级图像的邻域系统、观测特征场和分割标记场在过分割区域邻接图G上定义的实现方法为：在过分割区域邻接图G上定义对象级的观测特征场Y＝{Y_i|1≤i≤l}和对象级分割标记场X＝{X_i|1≤i≤l}，其中，表示过分割区域r_i的观测特征，y_s表示像素级图像位置索引s处的特征观测向量，|r_i|表示过分割区域r_i内的像素点个数；随机变量X_i∈K＝{1,2,…,k}，K为分割类别集合，k为预先给定的分割类别数；

过分割区域邻接图G＝(V,E)确定邻域系统：N＝{N_i|1≤i≤l}，其中，N_i＝{r_j|e_ij≠0}。

所述步骤二中利用混合高斯模型对观测特征场进行概率建模的方法为：在给定分割标记场X后，观测特征场Y中各个对象的特征向量之间相互独立，即：

其中，m×n表示像素级图像的像素个数，s表示像素级图像的像素的位置索引，P(Y＝y|X＝x)是给定观测数据Y＝y的条件下，对象级分割标记场X＝x的后验概率；x_s是位置索引s处的标记值，y_s是位置索引s处的影像特征的取值；

每一个对象所包含的像素的特征向量互相独立同分布服从高斯分布，且对象特征向量的似然函数等价于对象包含的像素的特征向量的似然函数乘积，即：

其中，x是分割标记场X的一个实现，p是特征向量的维数，μ_h表示分割类别h的特征均值，Σ_h表示分割类别h的特征协方差矩阵，使用极大似然法对各个类别的特征均值μ_h和特征协方差矩阵Σ_h的估计结果为：

k为分割类别数。

所述步骤二中利用Gibbs模型对分割标记场进行概率建模的方法为：根据Harmmersley-Clifford定理可知，具备Markov性的分割标记场的联合概率P(x)服从Gibbs分布，即：

其中，是一个归一化常量，x为分割标记场X的一个实现，Ω为所有实现x的集合，U(x)为能量函数：U(x)＝∑_c∈CV_c(x)，V_c(x)是集团的势；T为温度常量；

采用仅考虑双点基团的MLL方式去近似求解概率分布P(x)，能量函数U(x)可以通过下式来计算：

其中，V₂(x_i,x_j)表示过分割区域r_i和过分割区域r_j的势函数的取值，即：

β_c是势团参数，x_i,x_j分别表示过分割区域r_i对应的随机变量X_i和过分割区域r_j对应的随机变量X_j的实现；

由Markov性，在仅考虑双点的情况下得到局部概率的表达式：

其中，X是一个随机场，X_i代表过分割区域r_i标记取值的随机变量；N_i表示位置i的邻域，N_i＝{r_j|e_ij≠0}，U(x_i)表示标记场X_i的取值为x_i时的能量函数；

所述步骤二中根据Bayes准则得到分割标记场的后验分布模型的方法为：

由Bayes公式：

及最大后验概率准则得：

其中，y表示观测特征场Y的观测值，x表示标记场X的一个实现，y_i表示观测特征场Y_i的观测值，x_i表示标记场X_i的实现；

求分割标记的最佳结果转化为求分割标记场X后验分布最大化的问题，即：

其中，表示最终得到的分割标记场X_i的最优实现。

所述步骤三中根据最大后验概率准则运用ADMM算法对步骤二得到的分割标记场的后验分布模型进行迭代更新求解的方法为：

1)利用ICM算法实现像素级MRF方法得到每个像素点的分割所属类别：x^P＝{x_s|s∈S}，得到初始迭代的对象分割标记场实现X⁽⁰⁾＝{X₁ ⁽⁰⁾,X₂ ⁽⁰⁾,…,X_l ⁽⁰⁾}，其中，X_i ⁽⁰⁾＝mode(x_s|s∈r_i)，mode()为众数函数，x_s表示位置索引s的标记，X_i ⁽⁰⁾表示过分割区域r_i的标记；由对象级分割标记场在第t步的实现X^o(t)＝{X₁ ^o(t),X₂ ^o(t),…,X_l ^o(t)}得到对应各个类别的特征均值μ^(t)＝{μ₁ ^(t),μ₂ ^(t),…,μ_k ^(t)}和特征协方差Σ^(t)＝{Σ₁ ^(t),Σ₂ ^(t),…,Σ_k ^(t)}：

其中，X_i ^(t)表示第t次迭代时过分割区域r_i的标记、μ_h ^(t)、Σ_h ^(t)分别表示第t次迭代时第h类的均值向量和协方差矩阵；

2)对对象特征向量的似然函数即概率密度函数取对数，即计算每一像素位置对应的特征能量，每一个过分割区域r_i的特征能量为过分割区域内每一像素位置对应的特征能量的总和，即特征场能量为：

其中，μ_h为第h类的均值向量，Σ_h为第h类的协方差矩阵，y_s为位置索引s处的对象级的特征向量，()’为矩阵的转置；

分割标记场建模的局部概率取对数并去掉常数项得到过分割区域r_i的标记场能量为：

其中，X_i为过分割区域r_i的标记，X_i＝1,2,…,k，不同的类别取值对应不同的标记场能量；

将特征场能量和标记场能量分开考虑，将特征场标记记为h_y，标记场标记记为h_x，则特征场能量为：

标记场能量为：

EL(h_x)＝U(X_i)＝U(x)＝∑_cV₂(X_i,X_j)；

3)增加约束条件h_y-h_x＝0，步骤二中分割标记场的后验分布模型转化为ADMM算法的模型为：

其增广拉格朗日函数为：

4)增广拉格朗日函数L_ρ(h_y,h_x,)的ADMM算法的求解方法为：

h_y ^t+1：＝arg min_lL_ρ(h_y,h_x ^t,λt)，

h_x ^t+1：＝arg min_sL_ρ(h_y ^t+1,h_x,λ^t)，

λ^t+1：＝λ^t+ρ(h_y ^t+1-h_x ^t+1)，

其中，h_y、h_x、λ分别表示由观测特征场得到的标记结果、由分割标记场得到的标记结果以及对偶变量，h_x ^t、λ^t、h_y ^t+1、h_x ^t+1、λ^t+1分别表示h_y、h_x、λ在第t次和第t+1次的迭代结果。

本发明的有益效果：对似然函数和先验函数进行交替求解，充分考虑似然信息和先验信息，并增加相应的等式约束条件，通过对变量进行交替迭代来得到最优解；将特征场和标记场能量分开考虑，可更好地求解遥感影像分割标记，分割效率远远高于传统手工分割水平，提高了影像分割精度。本发明对影像分割后的实验图做定量分析可知，分割结果获得了最优或较优的kappa和OA值；因此，无论是定性分析还是定量分析，都验证了本发明在遥感影像分割中的有效性。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明的流程图。

图2为本发明一幅真实遥感图像及其直方图，其中，(a)表示真实遥感图像，(b)表示(a)的直方图。

图3为本发明初始化处理的示例图，其中，(a)为过分割结果图，(b)为(a)的局部图，(c)为(b)对应的过分割区域邻接图。

图4为本发明实验一影像图像处理的对比图，其中，(a1)为原彩色影像的灰度图，(a2)为真实手工分割结果图，(a3)为ICM方法分割结果图，(a4)为MRMRF方法分割结果图，(a5)为IRGS方法分割结果图，(a6)为OMRF方法分割结果图，(a7)为本发明方法分割结果图。

图5为本发明试验二影像处理对比图，其中，(b1)为原彩色影像的灰度图，(b2)为真实手工分割结果图，(b3)为ICM方法分割结果图，(b4)为MRMRF方法分割结果图，(b5)为IRGS方法分割结果图，(b6)为OMRF方法分割结果图，(b7)为本发明方法分割结果图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有付出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1所示，一种基于ADMM算法的马氏场下多分辨率遥感图像分割方法，运用一种新的方法即ADMM方法对模型进行求解，提高对遥感影像的分割精度。以下是具体的建模步骤：

步骤一：输入并读取一幅遥感图像，对读取的像素级图像进行初始化过分割，得到由过分割区域组成的对象级图像和对应的对象级过分割区域邻接图RAG，根据过分割区域邻接图RAG定义对象级图像的邻域系统N、对象级观测特征场Y和对象级分割标记场X。

步骤一的具体实施过程如下：

1)对输入的高空间分辨率三通道图像I(R,G,B)进行位置索引集合和像素级观测特征集合的定义，假设图像I(R,G,B)分辨率为m×n，得到位置索引集合S＝{s＝(x,y)|1≤x≤m,1≤y≤n}，像素级观测特征集Y⁰＝{y_s＝(r_s,g_s,b_s)}，其中，y_s表示位置索引s处像素点的观测特征值，也称为特征向量，r_s、g_s、b_s分别为位置索引s处图像的R、G、B分量的值，m为图像的长度，n为图像的宽度，(x,y)为图像中像素点的位置坐标，其位置索引为s，位置索引s的数量有m×n个。

2)根据设置的最小面积运用mean-shift方法对像素级图像进行过分割处理：将图像I(R,G,B)过分割为l个最小面积为S_min的过分割区域，每个过分割区域赋予标号，得到标号矩阵L_s＝{l_s|s∈S}，其中，元素l_s∈{1,2,…,l},s∈S；由此得到对象级图像的位置索引集合R＝{r₁,r₂,…,r_l}，其中，过分割区域r_i＝{s|l_s＝i,i＝1,2,…,l}，i表示过分割区域的标号值，i＝1,2,…,l；s表示像素级图像的位置索引，l表示过分割区域个数。

3)定义过分割区域邻接图(RAG)G＝(V,E)，其中，V＝{V₁,V₂,…,V_l}是位置集合，每个位置V_i代表对应的过分割区域r_i，E＝{e_ij|1≤i,j≤l}表示邻接关系，元素e_ij表示过分割区域r_i中与过分割区域r_j相邻的像素个数，e_ij≠0当且仅当过分割区域r_i和过分割区域r_j是相邻的。

如图3中(a)所示，首先用MATLAB程序得到过分割结果图，为了更好的说明区域邻接图的定义选取过分割结果图的一部分，如图3(b)所示，根据上述区域邻接图的定义，构造了图(b)的区域邻接图(c)，简单说明了图(b)中各过分割区域的邻接关系。

4)在过分割区域邻接图G上定义对象级的观测特征场Y＝{Y_i|1≤i≤l}和对象级分割标记场X＝{X_i|1≤i≤l}，其中，表示过分割区域r_i的观测特征，其中，y_s表示像素级图像位置索引s处的观测值，|r_i|表示过分割区域r_i内的像素点个数；X是一个随机场，X_i是一个随机变量，X_i∈K＝{1,2,…,k}，其中，K为分割类别集合，k为预先给定的分割类别数。

5)依据对象级的过分割区域邻接图G＝(V,E)确定邻域系统：N＝{N_i|1≤i≤l}，其中，N_i＝{r_j|e_ij≠0}。

步骤二：利用混合高斯模型对观测特征场Y进行概率建模得到特征场模型，利用Gibbs模型对分割标记场X进行概率建模得到标记场模型，根据Bayes准则联合特征场模型和标记场模型得到分割标记场X的后验分布模型。

步骤二的具体实施过程如下：

1)特征场建模

对特征场建模我们采用混合高斯(mixed-gaussion)模型，在过分割区域邻接图中，假设每一对象的观测特征在给定所属类别的条件下相互独立，对于对象级的观测特征场Y，大部分对象级马氏场方法均对其有两个假设：

假设一：在给定分割标记场X后，特征场中各个对象的特征向量之间相互独立，即：

其中，m×n表示输入图像的像素个数，s表示像素的位置索引，P(Y＝y|X＝x)是给定观测数据Y＝y的条件下，对象级的分割标记场X＝x的后验概率；x_s是位置索引s处的标记值，y_s是位置索引s处的影像特征的取值。

假设二：每一个对象所包含的像素的特征向量互相独立同分布服从高斯分布，且对象特征向量的似然函数等价于对象包含的像素的特征向量的似然函数乘积，即：

其中，s代表像素级图像中的位置索引，y_s代表位置索引s处的特征，x是分割标记场X的一个实现，p是特征向量的维数，μ_h表示分割类别h的特征均值，Σ_h表示分割类别h的特征协方差矩阵。在上述公式中，需要估计出特征均值μ_h和协方差矩阵Σ_h，使用极大似然法对这两个参数进行估计。各个类别的特征均值μ_h和特征协方差矩阵Σ_h的估计结果如下：

2)标记场建模

对标记场建模，采用的是Gibbs模型。在MRF模型中，假设标记场具备Markov性，根据Harmmersley-Clifford定理(已知格网的位置集合S具有邻域系统N，如果位置集合S上的随机场X是一个MRF模型。那么随机场X也是一个GRF)可知，标记场的联合概率P(x)服从Gibbs分布，即：

其中，是一个归一化常量，x为分割标记场X的一个实现，Ω为所有实现x的集合，U(x)为能量函数：U(x)＝∑_c∈CV_c(x)，V_c(x)是集团的势。T为温度常量，除非特殊需要一般都设为1，为了计算简便将其设为1。能量函数U(x)是势团集合C中所有势团的总和。对于各向同性的GRF来说，为了保证计算精度的同时，尽量使计算方便，采用仅考虑双点基团的MLL方式去近似求解概率分布P(x)，能量函数U(x)可以通过下式来计算：

其中，V₂(x_i,x_j)表示过分割区域r_i和过分割区域r_j的势函数的取值，即：

β_c是对应的势团参数，取值由人工给定，x_i,x_j分别表示过分割区域r_i和过分割区域r_j对应的随机变量X_i和X_j的一个实现。

由Markov性，可以在仅考虑双点的情况下得到局部概率的表达式：

其中，X是一个随机场，X_i是一个代表过分割区域r_i标记取值的随机变量；N_i表示位置i的邻域，N_i＝{r_j|e_ij≠0}，元素e_ij表示过分割区域r_i中与过分割区域r_j相邻的像素个数，e_ij≠0当且仅当过分割区域r_i和过分割区域r_hj是相邻的；U(x_i)表示标记场X_i的取值为x_i时的能量函数。

3)由Bayes公式：

及最大后验概率准则(MAP)得

其中，y表示观测特征场Y的观测值，x表示标记场X的一个实现，y_i表示观测特征场Y_i的观测值，x_i表示标记场X_i的实现。

所以，求分割标记的最佳结果就转化为求分割标记场X后验分布最大化的问题，即：

其中，表示最终得到的分割标记场X_i的最优实现。

通过循环迭代，更新分割标记，最终得到分割结果。

步骤三：根据最大后验概率准则运用ADMM算法对步骤二得到的分割标记场X的后验分布模型进行迭代更新求解，获取观测特征场和分割标记场的标记结果，将观测特征场的最终标记结果作为遥感图像的分割结果。

步骤三的具体循环迭代的具体过程如下：

1)首先由经典的ICM算法实现像素级MRF方法，得到每个像素点的分割所属类别，即像素级分割场结果：x^P＝{x_s|s∈S}，进而得到初始迭代的对象分割标记场实现X⁽⁰⁾＝{X₁ ⁽⁰⁾,₂ ⁽⁰⁾,…,X_l ⁽⁰⁾}，其中，X_i ⁽⁰⁾＝mode(x_s|s∈r_i)，mode为众数函数，x_s表示位置索引s的标记，X_i ⁽⁰⁾表示过分割区域r_i的标记。由对象级分割标记场在第t步的实现X^o(t)＝{X₁ ^o(t),X₂ ^{o (t)},…,X_l ^o(t)}，根据下式得到对应各个类别的特征均值μ^(t)＝{μ₁ ^(t),μ₂ ^(t),…,μ_k ^(t)}和特征协方差Σ^(t)＝{Σ₁ ^(t),Σ₂ ^(t),…,Σ_k ^(t)}：

其中，X_i ^(t)表示第t次迭代时过分割区域r_i的标记、μ_h ^(t)、Σ_h ^(t)分别表示第t次迭代时第h类的均值向量和协方差矩阵。

2)在计算概率密度时，为了避免出现过小的密度值，实际编程运算时常对概率密度函数取对数，即计算每一像素位置对应的特征能量，对于每一个过分割区域R_i的特征能量来说，为过分割区域内每一像素位置对应的特征能量的总和。假设第h类的均值向量为μ_h，协方差矩阵为Σ_h，位置索引s处的对象级的特征向量为y_s，则位置s处的能量值为：

Energy_s＝(y_s-μ_h)^-1Σ_h ^-1(y_s-μ_h)′+log|Σ_h|。

每一位置处都要进行相同的计算，从而计算出所有像素位置属于第h类的特征能量，对于每一个过分割区域r_i，其特征能量为：

由于均值向量μ_l和协方差矩阵Σ_h是第h类的均值向量和协方差矩阵，会随着类别标记h的变化而变化，也就是说，特征能量EO(h)的大小会随着类别标记的变化而变化，将h的值取遍所有的类别标记值后，能使特征能量EO(h)的值达到最小的标记值即为根据特征场能量得到的过分割区域r_i的标记。因此，可以将特征能量EO(h)看成是关于类别标记h的函数。

3)同样的，对分割标记场建模后在计算时取对数并去掉常数项得到过分割区域r_i的标记场能量为：

其中，X_i为过分割区域r_i的标记，X_i＝1,2,…,k,不同的类别取值对应不同的标记场能量，能使标记场能量EL(h)的值达到最小的标记值即为根据标记场能量得到的过分割区域r_i的标记。

4)原来对模型求解的方法都是将特征场能量和标记场能量加起来后对整体能量求最小，即求解arg min_hEO(h)+EL(h))问题，通过循环迭代，能使总能量达到最小的标记值h^*即为最终标记，现在引进ADMM算法，对优化问题：

其中，a与b是优化变量，f(a)+g(b)是待最小化的目标函数，它由与优化变量a相关的f(a)和与优化变量b相关的g(b)这两部分构成，Aa+Bb-C＝0是等式约束条件。其增广拉格朗日函数为：

其中，λ是对偶变量(或称为拉格朗日乘子)，ρ＞0是惩罚参数，L_ρ(a,b,λ)的增广是加入了二次惩罚项(ρ/2‖Ax+By-C||₂ ²。则该问题第t+1次的ADMM迭代求解方法为：

a^t+1：＝arg min_aL_ρ(a,b^t,λ^t)，

b^t+1：＝arg min_bL_ρ(a^t+1,b,λ^t)，

λ^t+1：＝λ^t+ρ(Aa^t+1+Bb^t+1-C)；

其中，b^t、λ^t分别表示待优化变量b和对偶变量λ的第t次迭代结果，a^t+1、b^t+1、λ^t+1分别表示待优化变量a、b及对偶变量λ的第t+1次迭代结果。

现在将原来放在一起考虑的特征场能量和标记场能量分开考虑，将特征场标记记为h_y，标记场标记记为h_x，则特征场能量为

标记场能量为：

EL(h)＝U(X_i)＝U(x)＝∑_cV₂(X_i,X_j)；

增加约束条件h_y-h_x＝0，则原有的问题可以转化为：

其增广拉格朗日函数为：

则上述问题的ADMM求解方法为：

h_y ^t+1：＝arg min_lL_ρ(h_y,h_x ^t,λ^t)，

h_x ^t+1：＝arg min_sL_ρ(h_y ^t+1,h_x,λ^t)，

λ^t+1：＝λ^t+ρ(h_y ^t+1-h_x ^t+1)，

其中，h_y、h_x、λ分别表示由观测特征场得到的标记结果、由分割标记场得到的标记结果以及对偶变量，h_x ^t、λ^t、h_y ^t+1、h_x ^t+1、λ^t+1分别表示h_y、h_x、λ在第t次和第t+1次的迭代结果。

关于评价标准，对于图像的分割问题，评价可以分为两种——定性评价和定量评价。定性评价就是视觉直观判断，看分割结果是否有较好的过分割区域性，是否具有较为准确的分区边缘，是否与人工判读一致。定量评价则是利用分割结果的具体数据算出各类具体的评价指标。其中常用的有Kappa系数，全局精度和错误率。本发明对于算法的定量评价指标为Kappa系数和全局精度。下面简单说明一下上述指标。

(1)Kappa系数

Kappa系数的指标主要用于一致性检验以及衡量分类精度。计算该指标首先需要计算混淆矩阵。以图像分割为M类为例，混淆矩阵可以表示为表1形式。在表1中，N_ij表示数据的真实标记是第i类，但算法将其分割为第j类。N_+j表示测试分割为第j类的总数，N_i+表示图像中真实分割应属于第i类的总数。基于混淆矩阵，可以得到Kappa系数的表达式：

其中，N表示影像的像素总个数。

表1混淆矩阵

(2)全局精度

该指标是最常用到的评价指标，就是被分类正确的数据的数目占所有数据总量的比例，即全局精度OA系数为：

为了验证本发明的效果，与以下四种分割方法做了对比。ICM(iterativecondition model)方法是一种像素级的“贪婪算法”，是一个基于局部条件概率的确定性算法，采用生成式模型，通过逐点更新影像标记完成影像分割。MRMRF(Multi-ResolutionMarkov Random Field)方法是一种与小波结合的算法，通过小波的金字塔结构，利用不同的尺度描述地物，扩大空间描述范围，且更有效的提取特征进行分析。IRGS(iterativeregion growing using semantics)方法是一种对象级的MRF算法，IRGS方法在构建空间上下文模型的过程中，引入了边缘强度信息，采用迭代过分割区域生长技术，以过分割区域为基元，得到较为准确的分割结果。OMRF(Object Markov Random Field)方法是另一种经典的对象级马氏场模型，与像素级不同的是，该方法将一个过分割区域视为一个基元来考虑更多的空间邻域关系，并采用生成式的迭代方法来求解模型。

上述四种方法有一个共同特点：标记场和特征场是一起考虑的。由于本发明将特征场与标记场分开求解，因此通过与这四种方法的对比，可以展示出本发明的优势。为了验证本发明的效果，用一幅拼接遥感图和一幅泰州过分割区域的航拍影像作为实验数据进行测试，其中实验一影像是一幅大小为256×256×3的遥感图，含有城镇、河流、森林和农田四种类别，如图4所示；实验二影像是一幅大小为512×512×3的遥感图，含有城镇、河流和农田三种类别，如图5所示。

本发明的运行平台是：core(TM)i5-4590@3.30GHz，RAM：4G，64位win7系统，2016a版matlab。为了验证本发明的有效性，对一幅拼接遥感图和一幅航拍遥感图分别做了一组实验，并且不同方法中相同参数设置是相同的。实验影像一的彩色影像的灰度图如图4(a1)所示，真实手工分割如图4(a2)所示。对实验影像一用ICM方法，使β＝10，得到的分割结果的灰度影像如图4(a3)所示。对实验影像一用MRMRF方法，使potential＝10，level-number＝3，wavelebase＝’haar’，得到的分割结果的灰度影像如图4(a4)所示。对实验影像一用IRGS方法，得到的分割结果的灰度影像如图4(a5)所示。对实验影像一用OMRF方法，并且使β＝10，MinRA＝50，得到的分割结果的灰度影像如图4(a6)所示。对实验影像一用本发明方法，并且使用β＝1、MinRA＝50、λ＝100、ρ＝0.4，得到的分割结果的灰度影像如图4(a7)所示。遥感影像实验二的彩色影像如图5(b1)所示，真实手工分割如图5(b2)所示。对实验影像二用ICM方法，使β＝10，得到的分割结果的灰度影像如图5(b3)所示。对实验影像二用MRMRF方法，使potential＝10，level-number＝3，wavelebase＝’haar’，得到的分割结果的灰度影像如图5(b4)所示。对实验影像二用IRGS方法，beta＝10，yi＝Icm(y,3,0.5)，得到的分割结果的灰度影像如图5(b5)所示。对实验影像二用OMRF方法，并且使β＝10，MinRA＝600，得到的分割结果的灰度影像如图5(b6)所示。对实验影像二用本发明方法，并且使用beta＝10，MinRA＝600,lambda＝100,sigma＝10，得到的分割结果的灰度影像如图5(b7)所示。实验影像一和实验影像二分割结果的Kappa系数如表2所示，分割结果的总精确度OA如表3所示。

表2 Kappa系数

表3 OA系数

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (7)

1.一种基于ADMM算法的马氏场下多分辨率遥感图像分割方法，其特征在于，其步骤如下：

步骤一：输入并读取一幅遥感图像，对读取的像素级图像进行初始化过分割，得到由过分割区域组成的对象级图像和对象级过分割区域邻接图，根据过分割区域邻接图定义对象级图像的邻域系统、观测特征场和分割标记场；

步骤二：利用混合高斯模型对观测特征场进行概率建模得到特征场模型，利用Gibbs模型对分割标记场进行概率建模得到标记场模型，根据Bayes准则联合特征场模型和标记场模型得到分割标记场的后验分布模型；

步骤三：根据最大后验概率准则运用ADMM算法对步骤二得到的分割标记场的后验分布模型进行迭代更新求解，获取观测特征场和分割标记场的标记结果，将观测特征场的最终标记结果作为遥感图像的分割结果。

2.根据权利要求1所述的基于ADMM算法的马氏场下多分辨率遥感图像分割方法，其特征在于，所述步骤一中运用mean-shift方法对像素级图像进行过分割处理并获得对象级过分割区域邻接图的方法为：

将分辨率为m×n的像素级图像I(R，G，B)过分割为l个最小面积为S_min的过分割区域，每个过分割区域赋予标号得到标号矩阵L_s＝{l_s|s∈S}，其中，元素l_s∈{1，2，...，l}，s∈S；位置索引集合S＝{s＝(x，y)|1≤x≤m，1≤y≤n}，s表示像素级图像的位置索引，位置索引s的数量有m×n个，(x，y)为像素级图像中像素点的位置坐标，m为像素级图像的长度，n为像素级图像的宽度；且像素级的观测特征集Y⁰＝{y_s＝(r_s，g_s，b_s)}，其中，y_s表示位置索引s处像素点的观测特征向量，r_s、g_s、b_s分别为位置索引s处像素级图像的R、G、B分量的值；

得到对象级图像的位置索引集合R＝{r₁，r₂，...，r_l}，其中，过分割区域r_i＝{s|l_s＝i，i＝1，2，...，l}，i表示过分割区域的标号值，i＝1，2，...，l；l表示过分割区域个数；

定义过分割区域邻接图G＝(V，E)，其中，V＝{V₁，V₂，...，V_l}是位置集合，每个位置V_i代表对应的过分割区域r_i，E＝{e_ij|1≤i，j≤l}表示邻接关系，元素e_ij表示过分割区域r_i中与过分割区域r_j相邻的像素个数，e_ij≠0当且仅当过分割区域r_i和过分割区域r_j是相邻的。

3.根据权利要求2所述的基于ADMM算法的马氏场下多分辨率遥感图像分割方法，其特征在于，所述对象级图像的邻域系统、观测特征场和分割标记场在过分割区域邻接图G上定义的实现方法为：在过分割区域邻接图G上定义对象级的观测特征场Y＝{Y_i|1≤i≤l}和对象级分割标记场X＝{X_i|1≤i≤l}，其中，表示过分割区域r_i的观测特征，y_s表示像素级图像位置索引s处的特征观测向量，|r_i|表示过分割区域r_i内的像素点个数；随机变量X_i∈K＝{1，2，...，k}，K为分割类别集合，k为预先给定的分割类别数；

过分割区域邻接图G＝(V，E)确定邻域系统：N＝{N_i|1≤i≤l}，其中，N_i＝{r_j|e_ij≠0}。

4.根据权利要求1或3所述的基于ADMM算法的马氏场下多分辨率遥感图像分割方法，其特征在于，所述步骤二中利用混合高斯模型对观测特征场进行概率建模的方法为：在给定分割标记场X后，观测特征场Y中各个对象的特征向量之间相互独立，即：

其中，m×n表示像素级图像的像素个数，s表示像素级图像的像素的位置索引，P(Y＝y|X＝x)是给定观测数据Y＝y的条件下，对象级分割标记场X＝x的后验概率；x_s是位置索引s处的标记值，y_s是位置索引s处的影像特征的取值；

每一个对象所包含的像素的特征向量互相独立同分布服从高斯分布，且对象特征向量的似然函数等价于对象包含的像素的特征向量的似然函数乘积，即：

其中，x是分割标记场X的一个实现，p是特征向量的维数，μ_h表示分割类别h的特征均值，∑_h表示分割类别h的特征协方差矩阵，使用极大似然法对各个类别的特征均值μ_h和特征协方差矩阵∑_h的估计结果为：

k为分割类别数。

5.根据权利要求4所述的基于ADMM算法的马氏场下多分辨率遥感图像分割方法，其特征在于，所述步骤二中利用Gibbs模型对分割标记场进行概率建模的方法为：根据Harmmersley-Clifford定理可知，具备Markov性的分割标记场的联合概率P(x)服从Gibbs分布，即：

其中，是一个归一化常量，x为分割标记场X的一个实现，Ω为所有实现x的集合，U(x)为能量函数：U(x)＝∑_c∈CV_c(x)，V_c(x)是集团的势；T为温度常量；

采用仅考虑双点基团的MLL方式去近似求解概率分布P(x)，能量函数U(x)可以通过下式来计算：

其中，V₂(x_i，x_j)表示过分割区域r_i和过分割区域r_j的势函数的取值，即：

β_c是势团参数，x_i，x_j分别表示过分割区域r_i对应的随机变量X_i和过分割区域r_j对应的随机变量X_j的实现；

由Markov性，在仅考虑双点的情况下得到局部概率的表达式：

其中，X是一个随机场，X_i代表过分割区域r_i标记取值的随机变量；N_i表示位置i的邻域，N_i＝{r_j|e_ij≠0}，U(x_i)表示标记场X_i的取值为x_i时的能量函数。

6.根据权利要求5所述的基于ADMM算法的马氏场下多分辨率遥感图像分割方法，其特征在于，所述步骤二中根据Bayes准则得到分割标记场的后验分布模型的方法为：

由Bayes公式：

及最大后验概率准则得：

其中，y表示观测特征场Y的观测值，x表示标记场X的一个实现，y_i表示观测特征场Y_i的观测值，x_i表示标记场X_i的实现；

求分割标记的最佳结果转化为求分割标记场X后验分布最大化的问题，即：

其中，表示最终得到的分割标记场X_i的最优实现。

7.根据权利要求1或6所述的基于ADMM算法的马氏场下多分辨率遥感图像分割方法，其特征在于，所述步骤三中根据最大后验概率准则运用ADMM算法对步骤二得到的分割标记场的后验分布模型进行迭代更新求解的方法为：

1)利用ICM算法实现像素级MRF方法得到每个像素点的分割所属类别：x^P＝{x_s|s∈S}，得到初始迭代的对象分割标记场实现X⁽⁰⁾＝{X₁ ⁽⁰⁾，X₂ ⁽⁰⁾，...，X_l ⁽⁰⁾}其中，X_i ⁽⁰⁾＝mode(x_s|s∈r_i)，mode()为众数函数，x_s表示位置索引s的标记，X_i ⁽⁰⁾表示过分割区域r_i的标记；由对象级分割标记场在第t步的实现X^o(t)＝{X₁ ^o(t)，X₂ ^o(t)，...，X_l ^o(t)}得到对应各个类别的特征均值μ^(t)＝{μ₁ ^(t)，μ₂ ^(t)，...，μ_k ^(t)}和特征协方差∑^(t)＝{∑₁ ^(t)，∑₂ ^(t)，...，∑_k ^(t)}：

其中，X_i ^(t)表示第t次迭代时过分割区域r_i的标记、μ_h ^(t)、∑_h ^(t)分别表示第t次迭代时第h类的均值向量和协方差矩阵；

2)对对象特征向量的似然函数即概率密度函数取对数，即计算每一像素位置对应的特征能量，每一个过分割区域r_i的特征能量为过分割区域内每一像素位置对应的特征能量的总和，即特征场能量为：

其中，μ_h为第h类的均值向量，∑_h为第h类的协方差矩阵，y_s为位置索引s处的对象级的特征向量，()’为矩阵的转置；

分割标记场建模的局部概率取对数并去掉常数项得到过分割区域r_i的标记场能量为：

其中，X_i为过分割区域r_i的标记，X_i＝1，2，...，k，不同的类别取值对应不同的标记场能量；

将特征场能量和标记场能量分开考虑，将特征场标记记为h_y，标记场标记记为h_x，则特征场能量为：

标记场能量为：

EL(h_x)＝U(X_i)＝U(x)＝∑_cV₂(X_i，X_j)；

3)增加约束条件h_y-h_x＝0，步骤二中分割标记场的后验分布模型转化为ADMM算法的模型为：

其增广拉格朗日函数为：

4)增广拉格朗日函数L_ρ(h_y，h_x，λ)的ADMM算法的求解方法为：

h_y ^t+1：＝arg min_l L_ρ(h_y，h_x ^t，λ^t)，

h_x ^t+1：＝arg min_s L_ρ(h_y ^t+1，h_x，λ^t)，

λ^t+1：＝λ^t+ρ(h_y ^t+1-h_x ^t+1)，

其中，h_y、h_x、λ分别表示由观测特征场得到的标记结果、由分割标记场得到的标记结果以及对偶变量，h_x ^t、λ^t、h_y ^t+1、h_x ^t+1、λ^t+1分别表示h_y、h_x、λ在第t次和第t+1次的迭代结果。