邻域规则下的遥感图像分类后处理方法研究

王 巍，郑新奇，原智远，张路路

(中国地质大学(北京)信息工程学院，北京100083)

地物光谱的多样性和复杂性，导致传统分类方法得到的分类图像不能达到应用目的，必须进行分类后处理工作，从而提高遥感影像分类的精度。初始分类图像中常见的缺陷是分布于分类图像上的孤立点、孔洞等，即所谓的“椒盐现象”或“椒盐噪声”，不仅影响遥感图像分类精度，而且在栅格转矢量过程中大幅增加了计算量。商业遥感影像处理软件为用户提供了基于低通滤波、基于光谱特征或是基于形态学算子的聚类处理方法［1］，可以将临近的、同种类型的小块分区进行聚类并合并，但存在聚类结果不理想的情况。国内部分学者针对分类图像中的椒盐噪声开展了相关研究，主要解决有两条途径:一是针对中、低分辨率影像，多结合地统计学的相关理论，如武文波等提出基于数学形态学的遥感影像分类后处理方法，并在ENDAS平台上进行了试验［2］;靳志宾等在k-NN分类方法中加入地理权重，能一定程度减少椒盐噪声的产生［3］;姚娜等将地统计学理论引入遥感临域，利用临近像元所蕴含的空间依赖性来改善初始分类结果［4］;基于地统计学的方法一般都需要训练样本去构建分类后处理的模型，因此其效果必然会受到样本容量的影响。另外一种解决途径是采取面向对象分类，如殷亚秋、李家国等借助易康软件，通过面向对象分类的方式进行水体信息提取，从而在很大程度上抑制了椒盐噪声［5］，但面向对象方法更多适用于高分辨率影像，针对中、低分辨率影像则偏重于像素级别的处理。

遥感分类图像的本质是对地物信息的一种数值模拟，各像素间具备空间依赖性，即对于图像某一区域而言，分散在一定距离内的像素点在数值特征上通常会呈现空间依赖性［6］。因此，针对中分辨率遥感影像，笔者提出一种基于邻域规则的遥感分类后处理方法，该方法基于地理学第一定律的思想，借鉴元胞自动机中Moore(摩尔)邻域这一概念，构建由中心像素及其周围相邻8个像素组成的“像素元胞”，尝试对“元胞”邻域制定聚类规则来实现对分类图像的分析和调整，力求去除分类图像的椒盐噪声。

一、理论和方法

1．地理学第一定律

地理学第一定律来自于美籍瑞士地理学家Waldo Tobler的观点:任何事物都与其他事物相联系，但临近的事物比较远的事物联系更为紧密［7-8］。地理学第一定律表明地理事物或属性在空间分布上具有相关性，遥感分类图像是对现实世界的数值模拟。根据地理学第一定律，理论上可以对分类图像采用基于邻域规则的聚类方式进行分类后处理，其重点在于规则的制定。

2．摩尔邻域

元胞自动机(CA)是一种时间、空间、状态都离散，空间相互作用和时间因果关系都为局部的网格动力学模型［9］。摩尔邻域的概念来自于元胞自动机。标准CA中，摩尔邻域由中心元胞周围相邻8个元胞组成。如果说标准摩尔邻域的扩展半径为1，那么扩展型摩尔邻域则以中心元胞开始，往外层进行扩展，例如半径为2的扩展型摩尔邻域中会存在5×5的结构，即存在25个“元胞”(如图1所示)。将这种结构放到分类图像上，表现为一种3×3或5×5像素大小的一块区域，聚类规则便基于图像的摩尔邻域这种结构而制定。需要注意的是，以中心像素向外扩展过程中，在图像的边缘无法获得完整的摩尔邻域。例如标准摩尔邻域由中心像素向周围扩展一次，所获取的摩尔邻域可能不足9个像素，而是仅有4个或6个像素。

图1摩尔邻域(左)和扩展摩尔邻域(右)

3．聚类规则

由于分类图像的DN值代表地物类型，在研究中，对于分类的图像主要根据其DN值进行处理，可以忽略决定颜色显示的RGB值。同时，为了提高运算效率，在对中心像素DN值的检测过程中，检测区域设置为标准摩尔邻域。在进行分类图像调整时，主要遵从以下几条既定的邻域规则。

(1)众数规则

在一个摩尔邻域内，某一DN值的频数较高，则将此DN值设置为该摩尔邻域的众数。众数规则对中心像素DN值的影响表现为中心像素DN值倾向于等于该摩尔邻域的众数。

(2)最邻近规则

中心像素与其临近像素DN值相等的可能性高，而与距离较远的像素DN值相等的可能性低。最邻近规则对中心像素DN值的影响表现为中心像素DN值倾向于等于与其临近像素的DN值。

(3)像素权重规则

研究中对分类图像的所有像素采用循环遍历的方式进行处理。在循环过程中，一部分像素已经经过分析处理，而另一部分像素则还未进行分析处理。因此，像素权重规则在这里主要指进行规则聚类时，经过分析处理的像素对中心像素的影响程度要高于未经过分析处理的像素。在摩尔邻域内，此规则体现出两个方面的含义:一是前一像素DN值权重高于后一像素;二是上一行像素DN值权重高于下一行像素。像素权重规则对中心像素DN值的影响表现为中心像素DN值倾向于等于与其临近的已经过分析处理的像素DN值。

(4)规则的差异性

前面提到过，图像边缘无法获取完整摩尔邻域。因此，在对分类图像进行循环遍历时，根据像素在图像上的不同部位，上述规则需要做一定程度的修改。规则的差异性主要有3类:一是在图像的4个角点处，此时摩尔邻域仅有4个像素;二是在图像边缘处，此时摩尔邻域有6个像素;三是循环至图像的内部，才能出现完整摩尔邻域，即9个像素。3种摩尔邻域情况如图2所示，图中黑色部分表示中心像素，灰色部分表示摩尔邻域。

图2 图像上不同位置的摩尔邻域示意图

二、数据与试验

1．数据简介

试验数据为一景成像于2013年10月的北京地区Landsat8中分辨率遥感影像，在ENVI5．1平台下对Landsat8影像进行数据预处理，过程包括辐射标定和大气校正;经过预处理的影像数据占用较大的硬盘空间，处理时间较长，需要进行试验区的挑选和影像裁剪。通过浏览北京市影像发现，密云地区地类种类较全，包含了植被、水体、土壤以及不透水面4类主要地类，适用于本次试验。通过矩形框对经过预处理的Landsat8影像进行裁剪，矩形裁剪框的范围是北纬40°30'14″至40°19'46″，东经116°40'50″至116°59'27″。研究区在北京市的地理位置如图3所示。

图3 研究区在北京市的地理位置

2．分类图像的获取

城市遥感中，许多学者借助一种“软”分类的思想，提出一种V-I-S分类模型用于城市景观的分析［10-11］。V-I-S指的是植被(vegetation)、不透水面(impervious surface)以及土壤(soil)。通过对试验区影像解译可以发现，水体在影像中占有一定的比例。因此，研究中对V-I-S模型进行了扩展，添加水体类型，构建V-I-S-W分类模型，并在ENVI平台中进行监督分类，分类方法选择最大似然法。获取的分类图像包含4类地物信息:植被、土壤、水体和不透水面。

借助ENVI平台实现初始分类图像的获取，主要步骤如下:

(1)训练样本的选择

使用ENVI提供的ROI工具进行训练样本的选取。表1提供的数据表明，样本可区分度均在1．80以上，满足软件计算要求。

表1 训练样本的可区分度

(2)基于最大似然法的分类

使用ENVI中Maximum Likelihood工具进行最大似然法分类，设定保存路径和文件名，其他参数采取默认设置。运行工具后，软件自动进行计算，并将初始分类图像保存到设置路径中。

3．规则的实现及应用

本文在IDL平台下实现文中提出的聚类规则。选择IDL的原因，是因为IDL不仅提供了功能强大图像处理函数库，而且与ENVI软件高度集成，使其功能得到进一步加强。通过IDL实现的程序包含下列功能:

(1)分类数据的输入和输出

ENVI标准格式的分类图像，其DN值、地物类型及显示颜色三者间的关联信息都写在图像的头文件中，而分析处理重点是仅包含DN值的二维图像(矩阵)，可以在IDL批处理模式中，调用ENVI内置函数进行相关数据的读取操作;写出数据时，ENVI和IDL也提供了多种写出方式，可以将数据文件和头文件分开输出到本地储存，也可以通过ENVI函数将数据文件和头文件一并写入本地存储。表2中列举了IDL代码中用到的数据读写函数及其功能。

表2 IDL批处理模式下用到的数据读写函数

(2)众数的求取

分类图像DN值的数值类型为长整形。IDL提供的Histogram()函数来实现直方图统计的功能，根据直方图就可以计算数组的众数。

在IDL编译器中，求取众数的实现方式如下面代码所示:

FUNCTION CAL_MODE，data

ht=HISTOGRAM(data，OMIN=omin)

in=WHERE(ht NE 0)

maxv=MAX(ht［in］，index)

RETURN，omin+WHERE(ht EQ maxv)

END

(3)构建摩尔邻域

对IDL二维数组中元素的访问方式除了行列序号外，借助IDL数组的位置索引机制也是一种非常便利的访问手段。IDL中规定:索引从0开始计数，由数组左上角至右下角，按列逐渐递增，因此二维数组最后一位元素的索引为二维数组元素个数减1。构建摩尔邻域时，借助索引机制，可以方便地获取摩尔邻域内像素的位置和DN值信息。

下面以图像的内部的摩尔邻域为例说明(如图4所示)。首先为摩尔邻域进行编码，序号从1—8，9号位是中心像素，用编码表示以中心像素向外进行摩尔扩展时可获取的像素位置信息;然后根据图像矩阵的行列号设定转换规则，将位置编码映射到图像矩阵像素的位置索引上;最后通过位置索引信息便可以获取摩尔邻域内所有像素的DN值。通过这种映射方式，只需获取位置编码、索引或DN值三者其中之一，就能查询其他两者的信息。

图4 摩尔邻域构建及映射关系

(4)情景推导

程序载入分类图像后，通过对图像所有像素循环遍历的方式来调整图像，根据既定的规则逐一进行处理，每次仅调整一个像素，笔者将这一过程命名为情景推导。表3是通过伪代码形式展示情景推导过程，其书写方式基于IDL语法规则，因篇幅原因，更加复杂的推导过程以Rule_Set替代。部分推导过程中，需要将待检测中心像素周围的像素作为新的中心像素进行规则判别，这一过程类似于函数递归调用，因此以RECURSION替代。另外，IDL中没有In运算符，查询功能在IDL代码中可以通过Where()函数实现;而or在IDL中是位运算符，实现条件选择功能可以通过if-else组合来完成。

表3 情景推导过程

三、试验结果

为了对比分析，研究中对初始分类图像应用ENVI软件的Clump Classes工具进行了聚类处理，形态学算子大小设置为3×3，与摩尔邻域大小相同。图5显示了影像不同区域上原始图像、初始分类图、Clump Classes工具处理结果以及本研究提出方法处理结果的对比，初始分类图像经过基于邻域规则的聚类处理后，椒盐噪声的抑制效果较为显著，且视觉上较Clump Classes工具的处理结果而言更加接近于原始影像呈现的状态。

分类精度方面，因缺少地表真实信息，在统计分类精度时采用人工解译获取地表真实ROI的方式，将ROI作为地表真实信息。用于精度验证的ROI数量分别是:植被668像素、土壤653像素、水体666像素以及不透水面693像素。所有操作均基于ENVI平台，精度统计结果见表4。通过结果可以看出，本次试验中，基于邻域规则的聚类方法处理的分类图像精度最高，达到95．63%，Kappa系数为0．941 8。

表4 分类精度统计

图5 部分区域影像及分类结果对比

四、结束语

通过试验结果可以看出，应用基于邻域规则的聚类方法对初始分类图像进行分类后处理，可以在很大程度上抑制椒盐噪声，提高分类精度，满足了实际应用的需求。同时，使用地表真实ROI的精度验证结果表明，基于邻域规则的聚类结果较之ENVI软件提供的Clump Classes工具的处理结果而言精度要更高，形态上更加接近于影像所呈现的状态;但是在地块边缘处理上，Clump Classes工具给出的结果更加平滑，能更好消除地块的锯齿状边缘。条件所限，两者都缺乏真实地物信息的检验，因而结果还有待进一步验证。

基于邻域规则的聚类方法关键点在于对规则的制定，规则越合理、推导越深入、复杂度越高，处理的结果应该更加精确;如果加入对地块边缘的处理规则，也能在一定程度上消除锯齿状边缘。在本研究中，规则设计是建立在对图像所有像素进行循环遍历处理的基础之上，随着规则复杂度的提高、计算量的增大，处理效率必然会受到影响。因此，如何设计行之有效、效率更高、适用性更强的规则将是未来的研究重点。

［1］ 邓书斌．ENVI遥感影像处理方法［M］．北京:科学出版社，2011．

［2］ 武文波，杨贵军，陈步尚．基于数学形态学遥感影像分类后优化处理［J］．辽宁工程技术大学学报(自然科学版)，2003，22(2):172-175．

［3］ 靳志宾，蒲英霞，陈刚，等．基于地理加权的k-NN高分辨率遥感影像分类算法改进［J］．遥感技术与应用，2013(1):97-102．

［4］ 姚娜，林宗坚，张锦雄，等．遥感影像分类后处理的地统计方法［J］．武汉大学学报(信息科学版)，2013，38(1):15-18．

［5］ 殷亚秋，李家国，余涛，等．基于高分辨率遥感影像的面向对象水体提取方法研究［J］．测绘通报，2015(1):81-85．

［6］ 张锦雄，GOODCHILD M F．野外空间采样的渐进式策略［J］．武汉大学学报(信息科学版)，2008，33(5):441-445．

［7］ TOBLER W．On the First Law of Geography:A Reply［J］．Annals of the Association of American Geographers，2011，94(2):304-310．

［8］ 孙俊，潘玉君，和瑞芳，等．地理学第一定律之争及其对地理学理论建设的启示［J］．地理研究，2012(10):1749-1763．

［9］ 黎夏，叶嘉安，刘小平，等．地理模拟系统:元胞自动机与多智能体［M］．北京:科学出版社，2007．

［10］ SETIAWANA H，MATHIEU R，et al．Assessing the Applicability of the V-I-SModel to Map Urban Land Use in the Developing World:Case Study of Yogyakarta，Indonesia［J］．Computers，Environment and Urban Systems，2006，30(4):503-522．

［11］ WENG Qihao，LU Dengsheng．Landscape as a Continuum:an Examination of the Urban Landscape Structures and Dynamics of Indianapolis City，1991—2000，by Using Satellite Images［J］．International Journal of Remote Sensing，2009，30(10):2547-2577．