苏腾飞， 张圣微， 李洪玉

(内蒙古农业大学水利与土木建筑工程学院，呼和浩特 010018)

0 引言

从遥感影像中快速、准确地提取农田信息，是农业遥感的基础研究内容。近年来，很多学者都指出利用高空间分辨率遥感影像(high resolution remote sensing image，HRI)开展农田监测的必要性[1-3]，其中一个重要理由是： 很多地区，尤其是中国的一些区域，农田平均面积较小[4]； 如果利用中等空间分辨率的遥感影像，如空间分辨率为250 m的MODIS或30 m的Landsat数据，极易因像元混合导致错误的分割结果。因此，发展适用于HRI的农田信息提取算法是非常必要的。

HRI可以捕获众多地物的细节信息，使得很多地物在影像中以面目标的形式呈现，这就使HRI像素之间产生了较高的空间关联性。如何有效利用像素之间的空间信息，是发展HRI解译算法的关键。马尔科夫随机场(Markov random field，MRF)作为一种影像分割算法模型，在其Gibbs势能函数中直接将像素的光谱和空间信息相关联，以提升像素类别场在空间上的连续性，因而受到很多学者的青睐。

然而，基于MRF影像分割算法的求解通常具有较高的计算复杂度。主流的MRF模型求解都采用了Gibbs势能函数最小化的策略[5-9]，通过逐像素的迭代计算，试图求出势能函数最低时所对应的分割结果。为了提高MRF模型的求解效率，很多学者将超像素的图模型引入到MRF模型中，以最大限度地减少MRF的求解空间，从而提高计算效率。孙巍等[7]将简单线性迭代聚类(simple linear iterated clustering，SLIC)超像素生成的结果用于条件随机场的建模； 刘磊等[8]将高阶势能引入到超像素MRF算法中； Yu等[9]借助区域生长的策略，也发展了一套基于超像素的MRF分割算法； Qin等[10]将Yu的算法拓展到多波段影像的分割中； 苏腾飞等[11]将全局最优合并的策略引入到Yu的算法中。

虽然超像素MRF模型具有求解效率高的优势，但目前主流的超像素生成算法难以将超像素的边界与影像中地物的边界完全地吻合，而这类错误往往增大了MRF求解的不确定性。超像素算法的目的是将影像分割为大小近似相等、形状大致相似的若干区域(或超像素)[12-14]。目前主流的超像素生成算法，例如SLIC[12]，Mean-Shift[13]，TurboPixel[14]等，虽然都可以取得较好的超像素分割结果，但是也难以保证超像素边界与影像中地物的实际边界完全符合，这类现象会导致超像素MRF算法在交互势函数最小化的过程中发生错误。近年来虽然出现了较多基于超像素MRF算法的研究，但是都忽略了这一问题[7-11]。

在农田HRI的分割中，边界信息是非常重要的。很多相邻的农田地块在HRI中具有十分相似的光谱特征(例如种植同类作物的农田)，而它们之间间隔的通常是灌渠或道路，在HRI中对这类地物的解译主要依靠的是边界信息[15]。因此，将超像素MRF算法应用于农田HRI分割时，如何最大限度地利用边界信息，并且避免超像素生成算法引入的边界不完全符合的问题，对于分割精度的提高将是非常关键的因素。

针对农田地区的HRI，提出了一种基于超像素MRF的HRI分割算法。为了避免超像素边界不完全符合造成的误差对分割精度产生的负面影响，提出一种超像素的高阶邻域模型，并重新定义了超像素MRF的交互势函数，以全面利用边界信息来提高分割精度。在本文模型的基础上，还给出了一种交互势函数调节参数的自动估计方法，以进一步提高精度。

1 算法原理

1.1 超像素MRF模型

遥感影像可被看作一个矩形格阵：S={s(i，j)|1≤i≤H，1≤j≤W}，其中i和j分别代表像素的列、行号，W和H分别为影像的总列、行数。在S的基础上，定义2个随机场X=(Xs)s∈S与Y=(Ys)s∈S。X为标号场，其元素x的值取自整数集合L={l|1≤l≤L′}；Y为观测场，其元素y由像素光谱值表示。在基于超像素的MRF模型中，S被划分为在空间上互不重叠的若干区域(超像素)：R={rt|1≤t≤M}，rt={(i，j)|l(i，j)=t，(i，j)∈S}。其中每个超像素rt包含若干像素，并且所有rt的并集等于S，M为总区域数。

在遥感影像分割问题中，已知的是影像的像素值，求解的是各个像素的类别。MRF将该问题建模，为求取P(x)P(y|x)最大化的问题，一般将其进行负对数变换，即

U(x|y)=-ln(P(x|y))=-ln(P(x))-ln(P(y|x))=U(x)+U(y|x)，

(1)

进而转化为势能函数最小化的问题，即

(2)

式中：P(x|y)是X在Y条件下的条件概率分布，也是MRF模型的目标函数；P(x)为X的先验分布；P(y|x)是Y在X条件下的条件概率分布；U(x)为一阶势函数，反映了X与Y在概率分布上的关联性；U(y|x)为交互势函数，表示X与Y的元素在空间分布上的关系； β是2个势能函数项的调节参数，需要根据具体影像来设置或估算。

一阶势函数一般采用高斯分布模型来建模[9-11]，其计算方法为： 对于某个区域i，其属于类别c的势能可由下式估计，即

(3)

在MRF迭代求解的初期，马氏距离常因协方差矩阵的估计偏差而产生较大的错误。因此本文在具体实现中采用了范数距离取对数的方式，即

(4)

式中： ||·||表示2阶范数；Z2同样为归一化常数，与Z1类似。取自然对数的目的是为了减少范数距离的数值变化范围。

1.2 高阶邻域交互势函数建模

传统的交互势函数主要利用的是边界强度等影像的空间背景信息，引导求解过程向正确的方向收敛，并保证在空间上相邻的像素或超像素具有一致的类别，因此该项又被称为类边界惩罚函数。基于超像素MRF的交互势函数的一般形式为

(5)

式中：rs与rt表示2个相邻的区域；g(·)表示2个区域公共边界的边界强度均值；xs与xt分别为2个区域的类别标号。

一般而言，传统计算方法仅考虑了2个相邻区域的公共边界的像素[7-10]。假如超像素过分割的结果与实际地物的边界完全符合(如图1(a))，则传统计算方法可以准确地反映2个区域之间的边界强度信息。但在很多实际情况中，超像素生成算法需要兼顾超像素内部灰度上的一致性和超像素之间在形状上的相似性，这使得超像素之间的边界像素难以完美地与实际地物边界符合(如图1(b))，此时传统计算方法难以准确反映2个区域的边界强度信息，图1(a)和(b)中的颜色反映的是地物的光谱值。

(a) 与实际地物边界完全符合(b) 与实际地物边界有误差 (c)r1和r2的高阶邻域范围

图1超像素边界示意图

Fig.1Illustrationofsuper-pixels’boundary

为了解决以上问题，本文提出了一种高阶邻域的交互势函数模型。该模型考虑了超像素之间高阶邻域像素的边界强度值，以更准确地描述地物边界信息，从而提高MRF的分割精度。该模型与式(5)类似，但增加了一个高阶邻域范围参数hn，即

(6)

与式(5)不同的是，式(6)中计算方法考虑了2个超像素边界像素的高阶邻域，如图1(c)，浅灰色表明了hn为3时高阶邻域范围内的像素，在此基础上计算式可表示为

(7)

式中：D(r1，r2)是超像素r1与r2公共边界的像素集合；N(p，hn)是像素p以hn为范围的高阶邻域集合像素；p和m分别表示属于集合D与N的某一像素；qm表示像素m处的边界强度值，其范围在(0，1)之间；d(r1，r2)表示超像素r1与r2公共边界的长度，以像素为单位； 系数3为经验常数，是根据多次实验确定的。本文的边界强度计算方法采用了文献[10]的矢量梯度方法，该方法可以从多波段影像中提取单波段的边界强度影像。

1.3 基于边界惩罚的β参数估计方法

β参数的作用是调节一阶势函数与交互势函数对MRF总势能函数的贡献。β可由用户根据实际情况来设置，但此方法受主观因素影响较大，且往往需要多次实验，费时费力。β还可根据影像中的统计特征量来估计[9-11]，但该方法需要较多的类别样本，而且对于一景影像的分割，该方法的β是固定的。根据式(5)和(6)，可知交互势函数主要惩罚的是类别不同区域的边界，而不同类别区域的边界具有不同程度的边界强度信息，因此在类别边界惩罚中应考虑到具体类别的统计信息。基于该思路，本文提出了一种自动的β估计方法，即

β(xs,xt)=ln(‖μs,t‖)，

(8)

式中μs,t为类别xs与xt的光谱均值向量相减所得的向量。需要注意的是，其计算结果与式(4)中的对数范数距离类似，因此，式(8)在考虑到具体类别统计信息的同时，也在数值上使交互势函数和一阶势函数更好地相互影响，从而更为精确地引导MRF分割算法的求解过程。

1.4 求解方法与算法流程

很多主流的MRF模型求解算法采用的是迭代条件模式(iterated conditional mode，ICM)[5-6]。该方法是在最大先验概率(maximum a priori，MAP)的理论上建立的。ICM通过逐像素的势能函数最小化迭代计算，试图求出Gibbs势能函数最低时所对应的分割结果。在本文超像素MRF模型的基础上，给出基于ICM的求解方法及具体流程。

输入包括原始HRI、单波段的边界强度影像、超像素生成算法的参数、类别数目参数c、高阶邻域范围参数hn和最大迭代次数Tmax。

具体流程如下：

1)利用超像素生成算法，对输入的HRI进行过分割。

2)利用最近邻域法，初始化各个超像素的类别。

3)设置迭代次数T=0，并进行迭代计算。①估计各个类别的光谱均值向量μ； ②根据式(8)计算出不同类别边界的参数β； ③利用式(2)计算各个超像素的类别； ④T加1。

4)若T≤Tmax，返回步骤3。

5)结束并输出结果。

由算法流程可知，在算法初始阶段，需对各个超像素的类别进行初始化。该步骤可用非监督聚类与大多数投票的方法来实现[11]，即先利用K-Means等算法给出粗略的聚类，再将每一个超像素类别初始化为其中大多数像素所属的类别。然而在实验中，该方法易受到K-Means初始化的影响而使后续步骤产生较差的分割效果。

因此本文利用了一种简单、高效的半监督策略来完成超像素类别的初始化，该方法包含3个步骤： ①用户在原始影像中选取样本，每一个类别仅选取5～10个具有代表性的像素作为样本； ②利用样本计算各个类别的光谱均值向量，作为每个类别的中心； ③根据各个类别的中心向量，利用最近邻域法对各个超像素进行类别初始化。实验表明，该方法可在用户干预较少的前提下保证较高的精度。

步骤3和4实际上是利用ICM实现的。值得一提的是，步骤3中第③步，在估计超像素类别的计算时，将式(4)，(6)和(8)带入式(2)，以估计本文超像素MRF模型的势能函数最小时所对应的超像素类别。

2 实验验证

2.1 实验数据

为了验证本文算法的性能，共选取了2景农田地区的高空间分辨率多光谱遥感影像，以开展分割实验。

本文进行实验的计算机配置为： Windows7，4 G内存，CPU： Intel I5 4200 m(2.5 GHz)。2景影像分别由中国的资源三号卫星和美国的OrbView-3卫星获取，以下对其分别简称为S1和S2，其基本信息见表1。

表1 实验影像数据基本信息Tab.1 Basic information of the image data used for experiment

2景影像均包含近红外、红光、绿光和蓝光4个波段。由于原始影像幅宽巨大，而本文实验的主要目的是验证算法对于农田地区的分割精度，因此在S1和S2中选取了具有代表性的子影像开展实验，子影像大小均为400像素×400像素，采用近红外(R)，红光(G)，绿光(B)波段假彩色合成影像，如图2所示。

(a) S1子影像(b) S2子影像

图2实验用的子影像

Fig.2Subsetsforexperiment

由图2可见，S1和S2分别展现了不同地貌特征的农田地区。S1覆盖了内蒙古锡林郭勒盟草原地区的农田，主要作物包括玉米和小麦，且当地的农业主要依赖雨水和河水灌溉。该区的农田平均面积较大，但农田内部的光谱变化也较为明显。另外，该区的草场由于牧草收割活动而呈现出明显的草场区域边界。S2为宁夏中卫市的郊区区域。当地属半干旱气候，依赖黄河水灌溉，盛产小麦和葵花等经济作物。相比S1，S2中农田的平均面积较小，但其农田内部的光谱均匀性更高。值得一提的是，S2中存在较多的休耕田地，且其内部光谱由于杂草的存在而呈现出较大变化。

为了开展算法精度的定量评价，以专家手动提取的分割结果作为基准，计算不同MRF算法分割结果的总精度(overall accuracy，OA)与Kappa系数进行定量对比分析。

2.2 参数hn对算法性能的影响

为了进一步分析本文算法性能(分割精度与计算速度)与高阶邻域参数hn的关系。研究了S1和S2分割实验的运算时间、OA和Kappa随hn的变化规律，如图3所示。hn的范围被设置为[1，15]。值得一提的是，当hn=1时，高阶邻域交互势函数与普通交互势函数的计算相同。

(a) S1子影像OA与hn(b) S2子影像OA与hn

(c) S1子影像Kappa与hn(d) S2子影像Kappa与hn

(e) S1子影像运算时间与hn(f) S2子影像运算时间与hn

图3本文算法性能与hn的关系

Fig.3Relationshipbetweentheperformanceoftheproposedalgorithmandhn

在分割精度方面，S1和S2子影像的最高精度均是在hn>1时得到的，这较为有力地证明了本文提出的高阶邻域交互势函数对于提升分割精度的作用。对于S1子影像，hn=3时精度达到最高，而S2子影像最高精度所对应的hn为6。这表明，对于不同的影像，其最佳hn的取值可能存在差异。另外一个较为明显的规律是，随着hn的增加，虽然S1和S2子影像的精度起初增加并达到最高值，但随后精度呈下降趋势，且在hn较大时S1和S2子影像的精度均低于hn=1时的精度。这在一定程度上说明，超像素边界像素的高阶邻域在一定范围内对提高分割精度具有积极作用，但是过大的高阶邻域范围反而会误导MRF的求解过程，以致精度降低。需要说明的是，本文高阶邻域交互势函数提出的目的，是为了降低超像素与实际地物边界不完全符合造成的误差对MRF求解的负面影响。在此前提下，hn应近似等于超像素分割与实际地物边界的误差。因此，过高或过低的hn均不利于边界强度等空间背景信息在MRF求解过程中起到积极作用。

在计算速度方面，S1和S2子影像均显示出了同一规律： 运算时间随hn的增大而变长。有趣的是，虽然S1和S2子影像的像素数目相同，超像素的数目也相近，但S2的运算时间明显长于S1。这主要是由于： 根据式(6)，只有在相邻超像素类别不同时，高阶邻域交互势函数的计算才会进行； S1中的地物平均面积较大，因此相邻超像素之间类别不同的情况较少，而S2中很多地物的面积较小，这使得相邻超像素之间类别不同的情况更多，因此交互势函数的计算次数更多，导致耗时更久。

2.3 算法对比分析

为了全面、客观地验证本文算法的性能，选取了另外3种超像素MRF算法，以开展对比实验。算法分别为： ①无β参数自动估计的本文算法(以下简称M1)，该算法也采用了高阶邻域交互势函数，且其hn的设置与本文算法相同； ②一种高阶MRF影像分割算法[8](以下简称M2)； ③一种融合了MRF与区域生长的多波段影像分割算法[10](以下简称M3)。本文算法简称为M0。进行M0与M1对比的目的是验证本文提出的β参数自动估计方法对分割精度的影响。

为了消除初始化对分割结果的影响，4种算法均采用SLIC生成的超像素，其结果见图4(a)和(b)。另外，4种算法求解过程中所需的边界强度影像，也都由同一方法计算[10]，见图4(c)和(d)。此外，4种算法在分割S1和S2子影像时，类别数目参数c均分别设置为5和4； 分割S1和S2子影像时的hn分别被设置为3和6；Tmax均设置为50，以平衡分割精度与收敛速度。

(a) S1子影像超像素分割(b) S2子影像超像素分割(c) S1子影像边界提取结果 (d) S2子影像边界提取结果

图4S1和S2子影像超像素分割与边界强度提取结果

Fig.4Super-pixelsegmentationresultsandedgestrengthextractionresultsofS1andS2subsets

S1子影像4种不同MRF分割算法的分割结果如图5所示。

(a) M0 (b) M1 (c) M2 (d) M3

图5S1子影像4种MRF算法的分割结果

Fig.5S1subsetsegmentationresultsof4differentMRFalgorithms

观察图5(a)可以发现，本文算法M0得到了较好的分割结果，其中大部分地物内部的类别具有较好的类别一致性。虽然个别超像素被分为错误的类别，例如上方中央位置的农田(这主要是由于该处较大的光谱变化增加了分割难度)。M1的分割结果与M0相近，但M1未能精确分割S1子影像下方中央位置牧草场和收割后草场的交界区域(图5(b))，这主要是由于2种地物的光谱较为接近，再加上该处超像素与实际地物边界未能完全符合，导致交互势函数难以正确引导分割结果，这也在一定程度上显示了本文β自动估计方法的优越性。通过观察图5(c)和(d)，可以发现M2和M3将很多河流周边的区域错分为河流，导致精度较低。这主要是由于M2和M3的交互势函数未能充分利用空间背景信息以及不同类别之间的β参数差异，导致分割效果欠佳。另外，M3的分割结果中很多玉米被错分为小麦，除了其交互势函数模型的欠缺外，其超像素合并的求解策略容易加大类间区分的难度，也是影响其精度的一个重要原因。

S2子影像4种算法的分割结果如图6所示。

(a) M0 (b) M1 (c) M2 (d) M3

图6S2子影像4种MRF算法的分割结果

Fig.6S2subsetsegmentationresultsof4differentMRFalgorithms

对比图6(a)和(b)可以发现，M0与M1的分割结果在总体上较为接近，其主要差异在于玉米和葵花这2个类别的区分上。由图2(b)可见，S2中的玉米和葵花具有较为相似的光谱特征，且这2个类别在空间上交错分布，增加了这2个类别区分的难度。通过仔细对比M0与M1的结果，可以观察到M0对于玉米和葵花的区分效果更佳，这一现象在S2子影像右下方的农田较为明显，M1将其中的很多葵花错分为玉米，而M0在该区的分割更为成功。尽管整体上M2的分割精度与M0十分接近，但图6(c)显示M2将一些玉米错分为葵花。这在一定程度上说明，M2虽然可以较为有效地区分易混淆的类别，但在S2的分割上其效果依旧次于本文算法。对于M3，图6(d)显示其分割结果中很多休耕地和农田被错分为城镇，分割效果较差。

为了进一步定量评价S1与S2子影像4种算法的分割精度，各算法的OA和Kappa如表2所示。

表2 S1与S2子影像4种算法的分割精度定量评价Tab.2 Quantitative evaluation of segmentationaccuracy by 4 algorithms for S1 and S2 subsets

从表2中可以看出，S1子影像M0的分割精度最高，M3的精度最差，M0和M1的分割精度较为接近，这与图5的观察是一致的； S2子影像虽然整体精度都较高，但同样是M0的分割精度最高，M3的精度最差，也定量地说明本文算法具有最佳的分割精度。以上实验结果和分析都较为有力地证明了本文算法的优势。对于农田区域的HRI，本文的超像素MRF模型可以更为准确地利用边界强度等空间背景信息，从而引导MRF求解过程得到更为精确的分割结果。

3 结论

本文发展了一种超像素MRF的影像分割算法，以提升高空间分辨率遥感数据在农田地区的信息提取精度。本文算法将高阶邻域模型引入到MRF的交互势函数中，并采用了逐类别的β参数估计方法。通过实验验证，得出结论如下：

1)将高阶邻域模型引入到超像素MRF影像分割算法中，可以抑制超像素与实际地物边界不完全符合造成的误差对MRF求解的负面影响，使MRF更充分地利用空间背景信息，从而提高分割精度。

2)提出β自动估计方法是基于范数距离计算的，相比于传统的马氏距离，该方法能够更为精确地估计不同类别信息对β参数调节的贡献，有利于提高分割精度。

通过对高阶邻域参数hn的分析表明，不同影像的最佳hn可能存在差异。因此在未来的工作中，发展最佳hn的自动估计方法是十分必要的。另外，本文仅考虑了农田地区的光学遥感影像； 如何将本文算法拓展到SAR影像或非农田区域的遥感影像，也是值得探索的方向。

[1] Löw F,Conrad C,Michel U.Decision fusion and non-parametric classifiers for land use mapping using multi-temporal RapidEye data[J].ISPRS Journal of Photogrammetry and Remote Sensing,2015,108:191-204.

[2] Kim H O,Yeom J M.Effect of red-edge and texture features for object-based paddy rice crop classification using RapidEye multi-spectral satellite image data[J].International Journal of Remote Sensing,2014,35(19):7046-7068.

[3] 杨闫君,占玉林,田庆久,等.基于GF-1/WFV NDVI时间序列数据的作物分类[J].农业工程学报,2015,31(24):155-161.

Yang Y J,Zhan Y L,Tian Q J,et al.Crop classification based on GF-1/WFV NDVI time series[J].Transactions of the Chinese Society of Agricultural Engineering,2015,31(24):155-161.

[4] Liu M W,Ozdogan M,Zhu X J.Crop type classification by simultaneous use of satellite images of different resolutions[J].IEEE Transactions on Geoscience and Remote Sensing,2014,52(6):3637-3649.

[5] Besag J.On the statistical analysis of dirty pictures[J].Journal of the Royal Statistical Society,Series B Methodological,1986,48(3):259-302.

[6] Kumar S,Hebert M.Discriminative random fields[J].International Journal of Computer Vision,2006,68(2):179-201.

[7] 孙 巍,郭 敏.基于SLIC与条件随机场的图像分割算法[J].计算机应用研究,2015,32(12):3817-3820,3824.

Sun W,Guo M.Image segmentation based on SLIC and conditional random field[J].Application Research of Computers,2015,32(12):3817-3820,3824.

[8] 刘 磊,石志国,宿浩茹,等.基于高阶马尔可夫随机场的图像分割[J].计算机研究与发展,2013,50(9):1933-1942.

Liu L,Shi Z G,Su H R,et al.Image segmentation based on higher order Markov random field[J].Journal of Computer Research and Development,2013,50(9):1933-1942.

[9] Yu Q Y,Clausi D A.IRGS:Image segmentation using edge penalties and region growing[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,2008,30(12):2126-2139.

[10] Qin A K,Clausi D A.Multivariate image segmentation using semantic region growing with adaptive edge penalty[J].IEEE Transactions on Image Processing,2010,19(8):2157-2170.

[11] 苏腾飞,李洪玉.一种结合GBM的MRF遥感图像分割算法[J].内蒙古农业大学学报(自然科学版),2015,36(1):143-149.

Su T F,Li H Y.GBM-combined MRF method for remote sensing image segmentation[J].Journal of Inner Mongolia Agricultural University(Natural Science Edition),2015,36(1):143-149.

[12] Achanta R,Shaji A,Smith K,et al.SLIC superpixels compared to state-of-the-art superpixel methods[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,2012,34(11):2274-2282.

[13] Comaniciu D,Meer P.Mean shift:A robust approach toward feature space analysis[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,2002,24(5):603-619.

[14] Levinshtein A,Stere A,Kutulakos K N,et al.TurboPixels:Fast superpixels using geometric flows[J].IEEE Transactions on Pattern Analysis and Machine Intelligence,2009,31(12):2290-2297.

[15] Su T F,Li H Y,Zhang S W,et al.Image segmentation using mean shift for extracting croplands from high-resolution remote sensing imagery[J].Remote Sensing Letters,2015,6(12):952-961.