亚像元定位之邻域像元选取问题研究

张 建 ，李宗南，谢 静

（1.华中农业大学 资源与环境学院，武汉 430070；2.农业部 长江中下游耕地保育重点实验室，武汉430070 3.华中农业大学 理学院，武汉 430070）

遥感卫星传感器成像机制所导致的混合像元普遍存在于各类遥感影像中，混合像元现象的存在导致了混合像元地类属性的不确定性、混合像元内部不同地类所占比率（丰度）的不确定性及其空间分布的不确定性等问题［1－2］.前两个问题归结起来是混合像元分解需要解决的内容，而混合像元分解从上世纪70年代混合像元分解问题开始备受关注［3］，到目前为止针对不同的混合像元成因已经建立了较为完善的像元分解模型［4］，通过在地物与光谱特征之间的定量解算，可以获得不同类型地物的比率.

混合像元分解虽然在一定程度上提高了遥感影像定性和定量解译精度，但对于同样甚至更为重要的混合像元内部各类地物空间分布信息却无法获取，而其对于以获取地物空间范围或分布特征等为目的的研究或应用又是极其重要的.相对于混合像元分解，亚像元定位问题的研究起步较晚，亚像元定位最早被用于在亚像元尺度下确定同质农业用地的精确边［5－6］，在此之后亚像元定位的基本概念由Atkinson于1997年正式提出［7］，通过引入美国地理学家W.R Toblerz在1970年提出来地理学第一定律即：“空间相关性理论”这个假设条件［8］，定义亚像元定位的准则及目标就是使地物分布空间相关性最大［9］.

空间相关性是目前大多数亚像元定位方法的理论基础，而距离是主要决定因素.距离决定了在亚像元邻域范围内参与其所属类别判定的像元位置和个数，从而也决定亚像元空间排布结果.因此针对亚像元定位过程中邻域像元选取问题，本文就常见的邻域选取方法进行介绍和对比试验，在此基础上展开分析和探讨.

1 理论基础和方法

1.1 理论基础

作为空间相关性理论最直接的应用，反距离加权方法是基于该理论亚像元定位方法中被用来具体制定空间化规则最早的一种方法，同时也具有明确的物理意义，在亚像元定位方法中具有代表性，本文以反距离模型为理论基础，比较不同邻域选取方法的效果.反距离模型基于相近相似的原理：即两个物体离得近，它们的性质就越相似，反之，离得越远则相似性就越小，它以预测点与样本点间的距离为权重进行加权平均，预测点越近的样本点赋予的权重就越大［10］.

式（1）中，r幂值是距离权重因子，用来控制样点距离对插值结果的影响程度.当r取较大值，则最近处样点对插值结果的影响加强，反之则减弱.通过对亚像元j成为第k类地物的概率即（k0）值进行比较，选取其中最大所属类别的概率来确定该亚像元的所属类别，进行亚像元空间化.式（2）中，x，y分别是样点的横坐标、纵坐标.

1.2 邻域像元选取方法

距离是空间相关性的重要条件，在以正方形为采样单元的遥感影像数据中，每个亚像元的邻域像元最多只有8个，但亚像元到这8个邻域像元的距离是不一样的，一些研究者认为只有距离最近的若干个邻域像元才会对中心像元所包含的亚像元位置有显著影响，而不是所有的8个邻域像元.因此在计算亚像元所属类别时所选取的邻域像元选择方式对于亚像元定位处理会产生关键影响.

下面就常用的邻域像元选取方式进行介绍.首先通过简单的图例来说明亚像元与邻域像元的关系，如图1所示，在一个3×3栅格单元中，中央的像元分成若干个亚像元，这里列出了2个不同尺度的情况：当尺度因子S＝2时，假设像元里面亚像元分布在四个象限，新的栅格变为原栅格大小的1／4如图1a所示；当尺度因子S＝3 时，假设像元里面亚像元分布在9个不同区域，新的栅格变为原栅格大小的1／9如图1b所示.

1.2.1 固定邻域像元选取 亚像元所属类别概率计算中的距离参数是通过计算它与其它相邻像元之间的欧式距离来判断的，距离越近，相关性越大，概率贡献值越大.通过计算图1中，分别计算亚像元A 的中心点到它的八个邻域像元中心点的欧氏距离L，并将进行排序，得到的结果是：

根据空间相关性在距离上的定义以及亚像元A 的中心点到它的8 个邻域像元中心点的欧氏距离L的排序结果，一些研究者通过距离的排序作为指标选择最近的5 个和选择3 个邻域像元，作为亚像元定位处理中影响亚像元A 的邻域像元［11］.

图1 基于空间距离的邻域选取Fig.1 Neighboring pixels choosing based on the distance

图2 基于可变邻域距离的邻域像元选取方法示意图Fig.2 Variable neighboring pixels choosing method（S＝4）

1.2.3 象限邻域像元选取 利用象限的概念Koen还提出一种象限邻域选取方法来确定亚像元的邻域像元［13］.象限邻域像元选取方法中将亚像元所在的像元中心作为象限中心，选取参与亚像元定位计算的邻域像元时只考虑和它在同一象限的.图3中选取了3种不同的尺度因子（S＝2，3，4）分别进行说明，亚像元和相邻像元用浅灰色、深灰色和黑色标识.根据邻域像元的定义，和亚像元同种颜色的邻域像元作为该亚像元的象限邻域像元.当尺度因子为偶数时，例如S＝4的时候，浅灰色亚像元的象限邻域像元是左上、左中和中上的邻域像元.而深灰色亚像元的象限邻域像元则是中上、右上和右中3个邻域像元.其中邻域像元中上既是浅灰色亚像元又是深灰色亚像元的象限邻域像元，所以用两种颜色间隔表示.

当尺度因子为奇数时，位于坐标轴上和坐标轴中心点的亚像元的邻域像元选取需要特别处理.Koen在文中没有具体定义尺度因子为奇数时的象限邻域选取方法［13］，本文基于象限邻域的思想，对尺度因子为奇数时的象限邻域选取方法进行了改进.如图3所示，当S＝3时，位于坐标轴上的亚像元，例如亚像元A 所对应的象限邻域像元时，将坐标轴整体顺时针旋转45°，定义右边对应的3个邻域像元作为邻域象限像元.当亚像元位于坐标轴中心，例如亚像元B 的时候，定义包含该亚像元所在像元的8 个邻域像元作为象限邻域像元.

图3 象限邻域像元选取示意图Fig.3 Quadrant neighboring pixels choosing method

2 实验设计与精度评定

2.1 模拟数据实验

基于上文所提到的邻域选取方法，分别进行实验和结果分析.实验数据采用具有不同方向和形状特点的4幅二值模拟数据，如图4所示.每幅模拟数据影像的大小是256×256个像元，其灰度值分别为255和0.

首先对原始影像进行局部均值重采样，这里使用的重采样尺度因子为S＝4，即4×4个像元合成为一个像元，得到对应的重采样结果，合成过程中灰度值255和0两种类型的像元比例被计入对应的混合像元中，用于模拟混合像元分解之后得到的各类地物的丰度数据.然后利用最大似然法对重样影像进行分类，得到重采样影像的硬分类结果.此外基于模拟数据的重采样结果，采用本文前面阐述的8邻域、5邻域、3邻域、可变邻域和象限邻域5种亚像元邻域像元选取方法分别进行实验.在此需要特别说明的是，在利用不同邻域选取方法对亚像元每种地物类别的所属概率进行计算后，选用的是同一排序方法对其进行比较来最终判定亚像元的所属类别.

图4 模拟图Fig.4 Four type synthetic imagery

图5 模拟图实验结果Fig.5 Results of four type synthetic imagery

2.2 实例数据实验

实例数据使用ETM＋影像，在影像中截取一块大小为800×760的影像作为分析区域，如图6所示，ETM＋多光谱影像分辨率为30 m，该影像覆盖武汉市大部分城区和一部分长江、汉江区域.

因为基于反距离模型的亚像元定位方法主要适用于遥感影像上大于像元采样尺度的影像目标，因此在实例实验中通过提取图6中长江、汉江区域来作为实验对象.将提取结果进行重采样到与MODIS 影像相近的分辨率（尺度因子S＝8，240m），使用该重采样结果在MODIS影像尺度上对长江、汉江区域进行亚像元定位，试图获取与ETM＋影像相同分辨率的亚像元定位结果.图7A是长江、汉江范围是对ETM＋影像进行分类后提取的结果，图7B是其重采样结果（S＝8），图7C 是硬分类结果，余下的是基于反距离模型，在不同邻域选取方法下的处理结果.

图6 原始假彩色合成影像（波段5，4，3）Fig.6 The raw remote sensing imagery

图7 实例数据处理结果Fig.7 Results of remote sensing imagery

2.3 实验结果精度评价方法

本文通过PCC（Percent Correctly Classified）、Kappa一致性系数和Adjusted Kappa（AK）系数3个精度评价指标对实验结果进行精度评定.其中，AK 系数是针对亚像元定位结果精度评定在Kappa系数的基础上改变而来的，它仅对混合像元所在位置的亚像元定位结果进行精度分析［14］.采用AK 系数的目的是排除纯像元参与结果精度评定，因为它的参与只会单纯提升Kappa值，掩盖了真实情况，选择AK 系数针对混合像元的亚像元分解结果进行精度评价，可以更加科学和直观分析实验结果.

3 实验结果与分析

3.1 模拟实验数据分析

从视觉效果上来看，原始模拟数据经过重采样后变得模糊，对该结果进行硬分类，可以看硬分类结果的边缘信息丢失的很严重，呈现严重的锯齿.对重采样后的模拟数据进行亚像元定位处理后，可发现使用各种邻域选取方法进行亚像元定位处理后的效果都要优于硬分类的结果.

通过表1、表2、表3 和表4，可以进一步发现无论选择何种邻域选取方法进行亚像元定位后的精度结果均比硬分类的结果要好，视觉感觉在精度评定结果中得到了印证.其中亚像元定位的方法在PCC总分类精度评定指标上一般都提高了1%左右，Kappa系数评定指标上一般都提高了5%左右，而去除了纯像元部分的AK 系数评定指标更加清楚的反映了，亚像元定位方法优于硬分类的结论，其中AK 系数最高提高了58%（字母模拟数据），最低也提高了40%（多边形模拟数据），而且可以发现，图形越是复杂，亚像元定位方法相对于硬分类方法的精度提高的越多.

进一步分析4个精度评定表，可以发现除了在多边形模拟数据结果中，采用可变邻域方法的结果精度上超过了8邻域方法，其它3种模拟数据结果都表明，8邻域方法在几种邻域选取方法中，精度最高.精度结果从高到低依次是八邻域、可变邻域、象限邻域、3邻域、5邻域.8邻域考虑了待定位亚像元所在的混合像元周围的全部邻域像元，最全面的考虑到了各个方向上的影响，所以结果最为客观可信.可变邻域、3邻域和5邻域方法利用比较待定位亚像元到邻域像元的距离，选取最近若干个像元作为参与亚像元定位的邻域像元.当对象比较简单的时候，舍弃较远的几个像元，对提高精度是有好处的，例如在多边形模拟数据的结果中可变邻域的精度结果就比8邻域的高.但是问题的关键是，在不知道影像分布规律的情况下，选择8邻域方法作为基于空间引力的亚像元定位的邻域选取方法最为合适和稳健.

表1 线性要素亚像元定位结果精度评价Tab.1 Accuracy results of the line synthetic imagery

表2 椭圆形要素亚像元定位结果精度评价Tab.2 Accuracy results of the ellipse synthetic imagery

表3 多边形要素亚像元定位结果精度评价Tab.3 Accuracy results of the polygon synthetic imagery

表4 字母要素亚像元定位结果精度评价Tab.4 Accuracy results of the character synthetic imagery

3.2 实例实验数据分析

实例实验中8邻域的邻域选择方法也获得了最好的精度，其中AK 系数达到了86.76%，如表5.对比亚像元定位结果和长江、汉江区域原始提取影像，也可以发现，除了一部分比较复杂的位置，出现了一些误差之外，大部分区域的亚像元定位效果较好.其中，汉江区域亚像元定位效果较差的原因是因为，汉江比较窄，在进行尺度为8的重采样之后，汉江区域的空间结构信息丢失非常严重，已经超出了基于反距离权重模型亚像元定位方法的适用范围，所以结果不理想.

表5 实例数据亚像元定位实验结果精度评定表Tab.5 Accuracy results of the remote sensing imagery

在此进一步借鉴Muslim 在其对海岸线进行亚像元定位结果的基于位置精度评定方法［15］，对上面提取的长江边界进行空间位置上的精度评定.评定方法主要是计算采用了8邻域方式亚像元定位处理生成的长江边界和原始影像上对应位置边界的偏差值.精度评定选取实验影像中所提取的长江区域的右岸边界进行分析.图8中表明了长江右岸边界中每一个边界点上亚像元定位结果与原始边界的偏差.

图8 边界位置偏差图Fig.8 The deviation of the edge

如图9，对边界位置偏差值进行统计，可以发现边界位置偏差值为零的占到了74%，偏差绝对值为1的占到了15%，其它的占到了11%.因为是模拟从MODIS影像尺度（250 m）到ETM＋影像尺度（30m）进行亚像元定位，尺度相隔较大，加上长江边界比较复杂，所以该精度还是可以接受的.

图9 边界位置偏差值分布统计图ig.9 The deviation distribution of the edge

4 结论与展望

本文针对遥感影像亚像元定位问题中邻域像元选取方法进行了介绍和对比分析，通过模拟数据和实例数据实验结果发现在大多数情况下，特别是影像内容较为复杂的时候，选择八邻域像元信息参与亚像元定位的计算，结果稳健，可以取得最好的效果.此外，在确定邻域选取方法之后对亚像元每种地物类别的所属概率进行计算还存在不同的概率排序方法，其对结果也存在较大的影响，值得进一步探讨.

［1］承继承，郭华东，史文中.遥感数据的不确定性问题［M］.北京：科学出版社，2004.

［2］凌 峰，吴胜军，肖 飞，等.遥感影像亚像元定位研究综述［J］.中国图形图像学报，2011，16（8）：1335－1345.

［3］Hallum C R.On a model for optimal proportions estimation for category mixture［C］∥Proceedings of Eighth Int’1 Symp，Remote Sensing of Environment，Ann arbor，Mi，1972：951－958.

［4］Charles Ichoku，Arnon Karnieli.A review of mixture modeling techniques for sub－pixel land cover estimation［J］.Remote Sensing Reviews，1996，13：161－186.

［5］Schneider W.Land use mapping with sub－pixel accuracy from LANDSAT－TM image data［C］∥Proceedings of 25th international symposium，remote sensing and global environmental change，Graz，Austria，1993.

［6］Flack J，Gahegan M，West G.The use of sub－pixel measures to improve the classification of remotely－sensed imagery of agricultural land［C］∥Proceedings of the 7th Australasian Remote Sensing Conference，Melbourne，Australia，1994：531－541.

［7］Atkinson P M.Mapping sub－pixel boundaries from remotely sensed images［M］.Innovations in GIS 4，London，UK：Taylor and Francis，1997：166－180.

［8］Tobler W.A computer movie simulating urban growth in the Detroit region［J］.Economic Geography，1970，46 （2）：234－240.

［9］李小文，曹春香，常超一.地理学第一定律与时空邻近度的提出［J］.自然杂志，2007，29（2）：69－72.

［10］Zhan Qingming，Molenaar Martien，Lucieer Arko.Pixel Unmixing at the Sub－pixel Scale Based on Land Cover Class Probabilities：Application to Urban Areas［J］.Uncertainty in Remote Sensing and GIS，2002：59－76.

［11］吴 柯.基于神经网络的混合像元分解与亚像元定位研究［D］.武汉：武汉大学博士学位论文，2008.

［12］Mertens K C，Baets B D，Verbeke L P C，Wulf Robert R De.A sub－pixel mapping algorithm based on sub－pixel spatial attraction models［J］.International Journal of Remote Sensing，2006，27（15）：3293－3310.

［13］Mertens K C，De Baets B，Verbeke L P C，et al.Direct sub－pixel mapping exploiting spatial dependence［C］∥Geoscience and Remote Sensing Symposium，2004.IGARSS＇04.Proceedings.IEEE International，20－24Sept.2004，5：3046－3049.

［14］Mertens K C，Verbeke L P C，Ducheyne E I，et al.Using Genetic Algorithms in Sub－pixel Mapping［J］.International Journal of Remote Sensing，2003，24（21）：4241－4247.

［15］Muslim A M，Foody Giles M，Atkinson Peter M.Localized soft classification for super－resolution mapping of the shoreline［J］.International Journal of Remote Sensing，2006，27（11）：2271－2285.