结合MST划分和RHMRF-FCM算法的高分辨率遥感图像分割

林文杰，李 玉，赵泉华

辽宁工程技术大学测绘与地理科学学院遥感科学与应用研究所，辽宁 阜新 123000

遥感图像分割是将图像域划分成具有不同土地利用和土地覆盖(land use and land cover,LULC)类型区域的过程，是实现遥感图像信息自动提取的关键步骤，具有重要意义[1]。相对于传统中、低分辨率遥感图像，高空间分辨遥感图像(下文简称高分遥感图像)蕴含更加丰富的地物几何、纹理和光谱信息，呈现更多的地物细节信息[2]。高分遥感图像这些特点在地物目标提取方面极具潜力，但也给传统的遥感图像分割方法带来一些新的挑战[3]：①在更精细的空间尺度下，高分遥感图像中同类地物的光谱信息呈现更强的异质性，使得传统分割算法结果存在大量分割噪声；②城市遥感图像中呈现大量的琐碎地物类型，这些地物在分割结果中表现为几何噪声，传统分割算法难以处理这些几何噪声问题；③高分遥感图像中存在大量的“同谱异物”和“同物异谱”现象，分别造成同类地物光谱测度的聚类性降低和同类地物的可区分性降低，使得传统分割算法效果不佳；④数据量呈几何级数增长，呈现明显海量性的特点，使传统分割算法难以处理。

针对高分遥感图像的新特点，大部分学者在构建遥感图像分割模型过程中，将空间信息作为一种重要的约束信息，衍生出大量顾及空间约束的高分遥感图像分割算法，并取得丰硕的研究成果[4-6]。最小生成树(minimum spanning tree,MST)作为图论中的一种经典模型，因其能完备、高效地表达遥感图像中蕴含的空间关联和光谱信息而受到广泛关注[7-14]。MST是无向邻接图(下文简称为图)在全局权和最小约束条件下的最优表达，是一种典型的全局最优化问题。图像的MST表达结果中同质区域表现为空间上的集聚特性，而同质区域间的边界表现为不连通或连通节点间权值较大。视觉上，图像MST呈现边缘不连通的特性，使其对同质区域的不规则边界表达效果较为理想。但是，MST易受图像噪声影响，如在噪声处的权值较大，导致基于阈值的MST划分方法的结果中存在严重过分割现象。

而对于高分遥感图像海量性数据的特点，扩展现有稳健快速的分割方法，是解决海量遥感图像处理的有效途径之一。模糊c均值(fuzzy c-means,FCM)聚类算法作为遥感图像分割的经典算法之一，具有稳健性强、快速、稳定和易于扩展等特点，在遥感图像分割中应用广泛。该类方法包含以区域和像素为基本分割单元的两种方式[15-16]，其中基于区域的分割方法因其具备一定的抗噪性和对目标边界不规则性的良好表达能力，且以区域作为决策单元，大大降低了分割决策的复杂性，在图像分割中应用广泛[8]。在基于区域和统计的遥感图像分割方法中，根据区域划分与分割过程之间的关系，可分为动态和静态划分两类。动态划分方法是指在分割过程中根据迭代求解的参数动态调整几何划分参数，该类方法不仅能达到区域划分逼近目标不规则边界的目的，而且能得到分布参数的最优估计。常用的几何划分方法有规则划分[5]、Voronoi划分[17]和泊松划分[18]等。典型的如文献[16]采用Voronoi划分作为基本的几何划分方法，在此基础上构建区域隐马尔可夫随机场(regional hidden Markov random field,RHMRF)，并定义基于区域的模糊目标函数，对全色高分辨率遥感图像取得较高的分割精度。虽然该方法对彩色图像或纹理图像中的简单几何边界目标取得较好结果，但它存在迭代收敛慢、计算量大以及对复杂非规则目标边界精细结构表达效果不佳等缺点。另外，动态划分方法需迭代划分图像域，导致该类方法计算时间较长，在一定程度限制了其推广应用。

不同于动态划分方法，静态划分方法只在算法运行初期根据相似性准则对图像域进行一次划分，并在该划分得到的子区域基础上构建图像分割模型。由于在模型求解过程中只需估计模型参数，从而避免了求解过程对图像域的反复划分过程，因此极大地节省了算法运行时间。但该类方法存在分割结果依赖于初始划分结果的问题，因此建立有效表达地物目标边界非规则性的划分模型是该类方法成败的关键。常用的静态划分方法有超像素划分[19]、Mean-Shift划分[20]和MST划分[21-22]等，这些算法大都利用图像的光谱信息作为划分的主要判据，虽然具有一定的抗噪性且能获得局部最优结果，但是对于高分辨率遥感中广泛存在的几何噪声，易造成严重的过分割现象。

针对以上问题，本文试图将MST作为静态划分方法对图像域进行划分，在此基础上构建区域HMRF模型和结合KL散度信息(Kullback-Leibler divergence information)与信息熵正则化项的FCM目标函数；然后采用求偏导的方式对各模型参数进行求解，从而得到全局最优分割结果。其中，MST静态划分包括图像的MST表达和在其基础上的顾及区域形状信息的同质区域划分方法，试图通过引入形状信息抑制划分的几何噪声。同时，本文方法将隐马尔可夫随机场[8](hidden Markov random field,HMRF)从建模像素间邻域关系扩展到建模同质区域间的邻域关系，更方便有效地获取模型参数估计和图像最优分割。

1 算法描述

1.1 图像MST表达与划分

令Z={Zs=Z(us,vs);(us,vs)∈P,s∈N}为定义在栅格P(图像域)上的离散随机场，其中，(us,vs)为像素s的平面坐标，s为像素索引，N={1,2,…,n}为像素索引集，n为像素个数，Zs={Zsb;b=1,2,…,B}为表征第s个像素光谱测度矢量的随机变量，即定义在P上的连续随机函数Z在(us,vs)上的采样，b为波段索引，B为波段个数。给定遥感图像z={zs=z(us,vs);(us,vs)∈P,s∈N}可视为离散随机场Z的一个实现，Z的所有可能实现构成的光谱矢量特征空间记为ΩZ。将像素映射为节点，并以节点8邻域作为节点间的邻接关系，构建图像的图G=(V,E,W)。其中V={vs=(us,vs,zs);s∈N}为节点集，E={ec=(vs,vs′);c=1,2，…,C,s∈N,s′∈Np(s)}为边集，C为边的条数，Np(s)为节点s的8邻域，W={wc=w(vs,vs′);c=1,2,…,C}为边集对应的权值矢量，wc定义了边ec连接的两个节点间光谱矢量的非相似性测度

(1)

式中，dss′和θss′分别为节点vs和vs′之间光谱矢量欧氏距离和夹角；t1和t2∈(0,1)分别用于调节dss′和θss′在权函数中的作用强度，本文中t1和t2分别取0.1和0.8。

图G的生成树T′=(V,E′)是G中具有树状结构的子图，即E′⊂E，且所有节点具有连通性。MST(T=(V,E′))是T中权重和最小的一棵，即图像的MST表达T去除了图像的图表达中的冗余连接边，且是在全局连接权和最小约束下的全局最优简化。另外，T有如下特性：边的个数|E′|=|V|-1，即边的个数比节点个数少1；T是无环图。因此，图像的MST表达具有如下优点：基于T上的操作时间复杂度为O(|V|)，即具有线性复杂度；T的每条边都是割边，即删除或者恢复一条连接边会增加或者减少一个划分区域，极大简化了区域划分操作的复杂度，且每个子区域都对应一个子树；目标的边界表现为沿边生成的特性，有利于表达目标的复杂不规则边界。目前常用的经典MST生成算法有两种：Prim和Kruskal，二者均基于贪心策略[23]。其中Prim算法是一种点驱动的MST生长算法，基于斐波那契堆(Fibonacci heap)和邻接的数据结构下的该算法时间复杂度为O(|E|+|V|log|V|)，适用于稠密图；Kruskal算法为边驱动算法，其时间复杂度为O(|E|log|E|)，适用于稀疏图。由图像8邻域构建的图G可认为是稀疏图，因此本文基于Kruskal算法构建图像的MST。

为解决高分辨率遥感中大量存在的几何噪声和同质区域内异质性增强的问题，借助最小异质区域划分的思想，在图像MST表达的基础上，将光谱测度相似性和划分区域的形状差异作为区域合并的准则。首先，假设每个节点均为一个同质区域，T的每条边ec∈E′都是割边，因此将其作为区域合并边。然后定义区域合并的准则，即最小异质区域准则：对于ec中连接的两个区域Pi、Pi′，根据判据h(Pi,Pi′)≤S决定是否合并这两个区域。最后，依次遍历每条割边，直至完成所有的区域合并。其中，S为合并阈值，定义了最小均质区域的尺度，合并判据函数h(Pi,Pi′)为区域间光谱测度相似性和形状参数的线性组合，定义为

hPi,Pi′=αhcolor+(1-α)hshape

(2)

式中，α∈[0,1]为光谱测度相似性权重；光谱测度相似性hcolor定义为

(3)

hshape=γhcompt+1-γhsmooth

(4)

(5)

(6)

式中，γ∈[0,1]为区域形状紧致度权重；lmerge、li和li′分别为合并后区域、区域i和i′的轮廓线长度；bmerge、bi和bi′分别为合并后区域、区域i和i′的外包矩形周长。MST静态划分过程如图1所示。其中，图1(a)为模拟图像，该图像由3个同质区域(区域1、2和3)构成，其中区域1包含区域4(1个像素)，用该区域模拟几何噪声；图1(b)为图1(a)的无向带权图，邻接关系为8邻域，且为了能清楚地体现边权特性，只标注边权大于20的边，由此图可见在图表达中，边缘像素的连接边具有较大边权值；图1(c)为生成的图像MST，可见图像MST能清楚地表达边界；图1(d)为采用MST静态划分得到的划分树，每个子树对应一个均质区域。所以，静态MST划分在较好地表达目标的不规则边界的同时也能有效处理几何噪声。

图1 MST均质区域划分Fig.1 Homogeneous region tessellation by MST

1.2 区域隐马尔可夫场模型

上述MST划分将图像域划分为若干个子区域，但是这些子区域的划分是局部最优的，而非全局最优的。所以，分割算法一般都需要采用额外聚类算法，以得到全局最优分割结果。为此，本文借助像素HMRF-FCM的图像分割思想，在MST划分的基础上构建RHMRF-FCM算法，以获得全局最优分割结果。

(7)

式中，β为标号场邻域作用强度

(8)

在光谱测度空间的随机场p(z)和子区域标号随机场p(x)的基础上，参考基于像素的HMRF模型定义，RMRF定义为

p(z,x)=pzxp(x)

(9)

式中，p(z|x)采用似然函数方式进行计算，即

(10)

根据HMRF定义，p(z)为光谱测度观测量的随机场，而p(x)为非观测状态量(隐变量)的MRF，故将观测量随机场与隐变量MRF构建的模型称为HMRF。假设同标号数据服从同一多值高斯分布，则

(11)

式中，μxi和Σxi分别为第xi类多值高斯分布的均值和协方差矩阵；θxi=μxi,Σxi；θ={θ1,θ2,…,θq}。

1.3 模糊目标函数

文献[24]提出一种基于像素的改进HMRF-FCM算法。基于该算法思想，定义如下基于区域的HMRF-FCM目标函数

(12)

式中，第2项为KL散度信息正则化项，用于抑制过分割；λ为模糊度参数；R={rij;i=1,2,…,m，j=1,2,…,q}为模糊关系矩阵，rij为模糊关系函数，刻画了第i个区域数据隶属于第j个聚类间的程度；η={ηij;i=1,2,…,m，j=1,2,…,q}为标号HMRF先验分布，ηij为区域标号HMRF模型的先验分布，定义为

(13)

Dij为第i个区域数据关于第j类的信息量和，定义为

(14)

基于区域划分的全局最优分割结果即为目标函数式(12)的最小化解。对模糊关系矩阵R的反模糊化就可以得到每个区域的明晰类属性。反模糊化目前一般采用最大化模糊关系的方式进行求解，即

(15)

1.4 模型参数求解

目标函数式(12)中的未知参数包括μj、Σj和rij，其中参数μj和Σj为关于目标函数的显式参数，因此分别通过求目标函数关于这两个参数的偏导数进行求解。结果为

(16)

(17)

(18)

1.5 算法流程

综合以上讨论，最优分割算法的流程可总结如下：

(1) 设置循环索引d=0、聚类数q、模糊度参数λ、迭代停止次数T。

(7) 如果达到停止迭代次数，则退出循环；否则，d=d+1并返回步骤(3)继续迭代。

2 试验与讨论

为验证本文图像分割算法的可行性和有效性，在Intel Xeon CPU E5-2670 2.6 GHz/8 G内存/Matlab环境下，对WorldView3 Pan-Sharpen 0.5 m分辨率遥感图像进行了分割试验。如图2所示，图2(a)为512×512的待分割图像；图2(b)为人工判读的分割结果图，包含5个同质区域，分别由灰度值为28(Ⅰ)、75(Ⅱ)、112(Ⅲ)、140(Ⅳ)和255(Ⅴ)5种灰度标识，将其作为标准分割结果，并用于定量评价。首先，为了验证本文算法的有效性和可靠性，设计两个对比试验：①定量分析静态划分算法中尺度参数S、光谱测度参数α和区域形状紧致度权重γ对最终分割结果的影响；②定量分析本文算法和eCognition软件中的多分辨率分割方法以及分水岭算法。其次，为验证本文算法的实用性，对整景WorldView 3遥感图像进行试验。另外，提出本文算法各个参数统一设置如下：FCM停止迭代次数T=30；模糊度参数λ=0.1；标号场邻域作用强度β=0.1。

图2 试验数据和标准分割结果Fig.2 Experiment image and standard segment

2.1 静态划分参数对最终分割结果的影响

首先，为了定性和定量评价分割尺度对分割

结果的影响，分别对S=40、100、160和220，α=0.9，γ=0.1进行分割试验，试验结果分别如图3(a1)—(a2)、(b1)—(b2)、(c1)—(c2)和(d1)—(d2)所示。其中，图3(a1)—(d1)为各静态划分参数下的划分结果。图3(a2)—(d2)为对应的RHMRF-FCM算法最优分割结果。对比图3(a1)—(d1)的静态划分结果，随着S的增大，静态划分结果中的过划分现象逐渐减少；对比图3(a2)—(d2)的最优分割结果，随着尺度参数的增大，分割结果中包含的几何噪声逐渐减少。从视觉上看，图3(a2)—(d2)中的最优分割结果的复杂地物边界表达较为完整、准确，如图3黑框部分所示，从而验证了RHMRF-FCM算法的有效性。

为了定量评价4个尺度参数S=40、100、160和220的分割精度，采用混淆矩阵[25]作为定量分割精度评价模型，评价结果如表1所示。由表1可知，在S为40时，同质区域Ⅲ的用户精度低于0.7。理论上，这是由于在该尺度参数下，初始静态分割结果过分割较为严重，对最优分割过程造成一定程度的扰动；其他尺度参数的各精度指标均在0.7以上，且总精度和Kappa值均优于0.8，说明当尺度参数取[100,220]时，对分割结果影响较小。另外，在S=100下的各项精度指标总体最优。因此，将S=100作为最优划分尺度，进一步探讨在该尺度参数下，α和γ对静态划分结果和最优分割结果的影响。

图3 不同尺度下的分割结果Fig.3 Segmentation result from different S

表1 不同S分割结果精度评价

其次，为评价不同α对分割结果的影响，分别对α=0.9、0.8和0.7(S=100，γ=0.1)进行分割试验，试验结果分别如图3(b1)—(b2)、图4(a1)—(a2)和(b1)—(b2)所示。其中，图3(b1)和图4(a1)—(b1)为各静态划分参数下的划分结果；图3(b2)和图4(a2)—(b2)为对应的RHMRF-FCM算法最优分割结果。视觉上，对比各静态划分结果，随着α的减小，静态划分结果的过划分现象略微增加，最优分割结果中的复杂地物边界表达较为完整、准确，如图3(b1)—(b2)和图4中的黑框所示。

为了定量评价取上述α值所得的分割精度，得到如表2所示的评价结果。由表2可知，各α参数的精度指标均在0.8以上，且相互间的差距较小，这说明当α∈[0.7,0.9]时，本文算法的静态划分结果对最优分割过程影响较小，从而验证算法的稳健性。另外，在α=0.8时的各项精度指标最优。因此，将α=0.8作为最优光谱测度参数，进一步探讨在该参数下，γ对静态划分结果和最优分割结果的影响。

最后，为评价不同区域形状紧致度参数γ对分割结果的影响，分别对γ=0.1、0.2、0.3(S=100，α=0.8)进行分割试验，试验结果如图4(a1)—(a2)、图5(a1)—(a2)和(b1)—(b2)所示。其中，图4(a1)、图5(a1)—(b1)分别为各静态划分参数下的划分结果。图4(a2)和图5(a2)—(b2)分别为对应的RHMRF-FCM算法最优分割结果。视觉上，静态划分结果的过划分现象随着光谱测度参数γ的增大而轻微增大，但对最终的分割结果影响不是很明显。图5的定量评价结果如表3所示，表中列出各γ参数的各项精度指标。在γ=0.1时的各项精度指标最优。

图4 S=100、γ=0.1下不同α的分割结果Fig.4 Segmentation result from different α with S=100 and γ=0.1

表2 不同α下的分割结果精度评价

综上所述，初始分割参数中的尺度参数对初始划分结果影响最为明显，而α与γ对初始划分结果影响较小。试验结果表明，S、α和γ分别取[100,220]、[0.7,0.9]和[0.1,0.3]时，尺度参数对最终的分割结果影响较小，从而证明本文算法的稳健性。

图5 S=100、α=0.8下不同γ的分割结果Fig.5 Segmentation result from different γ with S=100 and α= 0.8

2.2 对比试验

为验证本文算法的有效性，将本文算法与eCognition 9.0软件中的多分辨率分割方法、分水岭算法结果做定性、定量分析。其中，eCognition软件中的多分辨率分割方法和分水岭算法都是目

前高分辨图像分割中常用的方法，但是这两种算法都只包含静态划分部分，并未包含最优分割部分。因此，本试验只用eCognition软件处理得到多分辨率分割方法和分水岭算法的静态划分结果，并导入Matlab中，用本文提出的RHMRF-FCM算法处理得到最优分割结果。算法结果如图6所示，其中图6(a1)—(a2)为本文算法在S=100、α=0.8和γ=0.1下处理得到的分割结果；图6(b1)—(b2)为eCognition软件多分辨率分割算法在相同参数下的分割结果；图6(c1)—(c2)为eCognition软件分水岭算法在长度因子15下的分割结果。

表3 不同γ下的分割结果精度评价

图6 对比算法分割结果Fig. 6 Segmentation results from comparing algorithms

从静态划分结果上看，3个算法都能较好地拟合地物的不规则边界；从划分对象的个数上看，多分辨分割算法划分的子区域个数最少，本文算法个数偏多，而分水岭算法的个数太多。但是，对比3个算法的划分结果，依然可以从中发现一些细微的差别。如，本文算法对图像中的地物边缘特征保持得更好，而多分辨率分割方法对于一些弱对比度的边缘保持的效果较差，这是由于本文算法借助图像的MST表达提高了对地物不规则边缘的辨识。而本文算法和多分辨率分割方法均采用顾及光谱测度相似性和形状参数的划分准则，都能较为理想地划分出城市场景中复杂建筑物边缘，而分水岭算法则未考虑形状参数，也在一定程度上使得该算法对大的复杂几何形体划分上呈现过分割的特性。另外，对比3个算法的最优分割结果(图6(a2)—(c2))，视觉上看复杂地物边界表达较为完整，这也说明了本文提出的最优分割算法的有效性。

3种对比算法的精度评价结果如表4所示，对比发现本文算法的各项精度指标均优于对比算法，且本文算法的总精度相对于多分辨率分割方法提高了18%，相对于分水岭算法提高了15%。

表4 对比算法分割结果精度评价

2.3 整景图像分割试验

为了验证本文算法的实用性，对1景8192×8192的WorldView-3 0.5 m高分遥感图像进行分割试验。分割参数设置为：S=100、α=0.8、γ=0.1和q=6，算法结果如图7所示。其中，图7(a1)为原图，黑框区域为局部放大图位置，如图7(b1)所示；图7(a2)为原图与划分结果边界叠加图，包含128 587个均质区域；图7(a3)为最优分割结果；图7(b2)—(b3)分别为图7(b1)的分割结果。视觉上看，静态划分结果能较好地划分均质区域，且对复杂地物边界表达较为完整，最优分割结果中各同质区域划分较为合理。算法总体运行时间大约为3 h，对于6千万以上节点规模的MST求解来说，时间上可以接受，从而证明了本文算法的实用性。

图7 整景图像分割结果Fig.7 Segmentation results of whole scene image

3 结 论

高分辨遥感图像中包含大量的几何噪声，具有地物目标的光谱测度异质性增强等特点，采用传统以像素为基本决策单元的分割方法，容易造成分割结果中包含大量噪声的问题，严重影响分割结果的高层分析应用。为此，本文提出一种结合静态MST划分和区域HMRF-FCM算法的高分辨率遥感图像分割方法。该方法在图像MST表达的基础上，将形状信息引入同质区域的划分中，在有效抑制几何噪声的同时能较好地表达地物目标的不规则边界。在此基础上结合区域HMRF模型构建的模糊目标函数，能够有效克服地物目标光谱测度异质性高的问题。为了验证本文算法的有效性和可行性，以WorldView-3高分遥感图像为试验数据，探讨了静态划分算法中尺度参数、光谱测度参数、区域形状紧致度权重对最终分割结果的影响，以及对比分析了本文算法和eCognition软件中的多分辨率分割方法、分水岭分割方法的分割。试验结果表明，本文算法能有效抑制分割噪声，且能较好地拟合地物的不规则边界，从而自适应、快速有效地分割遥感图像，相对于对比算法，分割精度提高了10%以上。由于该算法较为高效，因此在今后的工作中将该方法应用于大尺度高分辨遥感图像的分割，并设计相应的并行算法，以进一步提高算法的实用性。