基于子像元交换算法和地形数据的土地覆盖亚像元制图研究

虞舟鲁， 王文超， 戎 奕， 沈掌泉

(1.浙江大学环境与资源学院,杭州 310058； 2.杭州学联土地规划设计咨询有限公司,杭州 310030)

基于子像元交换算法和地形数据的土地覆盖亚像元制图研究

虞舟鲁1， 王文超2， 戎 奕2， 沈掌泉1

(1.浙江大学环境与资源学院,杭州 310058； 2.杭州学联土地规划设计咨询有限公司,杭州 310030)

在将遥感影像应用于土地覆盖制图的过程中，混合像元会产生负面影响，尤其是在空间分辨率较低的情况下。软分类和超分辨率(亚像元/子像元)制图技术可以解决上述问题。子像元交换算法是一种简单而有效的亚像元制图技术，但也存在计算效率不够高和在超分辨率制图因子(scale factor,S)较大时制图精度较差的问题，其原因可能是亚像元制图只是从软分类结果中获得信息。因此，将数字高程模型(digital elevation model,DEM)及其衍生的数据作为子像元交换算法的辅助信息，研究提高其制图精度的有效性。研究结果表明： ①在DEM信息的辅助下，亚像元制图的精度得到了改善，即使在S较大时也是有效的； ②同时使用高程和坡度信息时的效果要好于单独应用高程或坡度信息； ③在S较大的情况下，制图精度对邻域范围大小的敏感性要低于S较小时； ④在DEM信息的辅助下，计算效率得到提高。实验结果表明，将DEM作为辅助信息进行亚像元制图是有效和可行的。

亚像元制图； 超分辨率制图； 子像元交换算法； 遥感； 数字高程模型(DEM)； 土地覆盖制图

0 引言

土地覆盖是生态系统中生物化学循环、水文变化及地表与大气之间相互作用的体现[1]，其信息可用于理解和管理生态环境。因此，获得准确和现势的土地覆盖信息在土地管理、规划和监测中是非常重要的。由于遥感技术具有大范围获取信息的能力，因此遥感影像是获得土地覆盖信息的有效来源[2]。但是，卫星遥感影像普遍存在混合像元的问题，特别是在空间分辨率较低的情况下[3]。因此，在利用遥感数据提取土地覆盖信息时，混合像元常成为一个不容忽视的问题。在处理混合像元时，与硬分类技术相比，软分类技术可以避免信息的损失，因为其可以获得各土地覆盖类型在像元中所占的比例，并形成一系列对应各土地覆盖类型的比例图像[4-5]。可是，软分类形成的比例图像并没有明确各类型在像元内的空间位置，而在很多应用中，确定各类型在像元内的空间位置并获得详细的土地覆盖图是非常必要和重要的。Atkinson[6]提出的亚像元制图(亦称子像元制图、超分辨率制图或锐化)的概念认为，根据软分类所获得的比例图像形成锐化的土地覆盖图，可以综合发挥软分类和硬分类的优势。目前，学者们已经提出了一系列的亚像元制图技术，诸如子像元交换算法[7-8]、图像锐化技术[9-10]、基于知识的分析技术[11]、Hopfield神经网络[1,12-13]、反卷积滤波器[14]、线性优化技术[15]、遗传算法[4]、前馈神经网络[16]、基于马尔可夫随机场的技术[2]、基于子像元/像元空间吸引模型的技术[17]、基于子像元/像元空间吸引模型的改进子像元交换算法[18]、子像元/像元空间吸引模型与单亲遗传算法相结合的技术[19]等。子像元交换算法具有算法简单和效率较高的特点，已成功地应用于马来西亚的海岸线制图[20]和英国Dorset的Christchurch地区的农村土地覆盖制图[21]； 但是当超分辨率因子(scale factor, S)较大时(即在空间分辨率较低的情况下)，其计算效率和制图精度就会下降，其原因可能是只依靠软分类结果中各类型的比例信息和子像元类型的随机初始化所导致的大量子像元交换。另外，子像元类型的随机初始化也会对最终的制图结果产生影响。因此，通过辅助信息来改进该算法的初始化过程，有可能提高其计算效率和制图精度。

地形是基本的地球物理形态，包含了一个地区地壳变化、地质构造和气候变化等丰富的信息，在土地覆盖分类中经常作为辅助的信息[22]。详细和准确的数字高程模型(digital elevation model,DEM)信息可通过SPOT及ASTER等遥感数据获取，也可通过机载激光扫描系统(light detecting and ranging,LiDAR)获得，美国国家航空航天局(National Aeronautics and Space Administration,NASA)则通过航天飞机雷达地形测绘任务(shuttle Radar topography mission,SRTM)获得了全球大部分陆地的地形信息。Nguyen等[22]提出将LiDAR获得的DEM作为辅助信息通过Hopfield神经网络(Hopfield neural network,HNN)进行超分辨率制图，其研究结果表明，采用LiDAR获得的DEM数据结合光学遥感数据，HNN可以在子像元尺度上准确地预测土地覆盖类型，使制图精度明显提升，尤其是建筑物制图。其应用的最大问题是在研究区需要有可用的LiDAR数据，否则无法将其应用于通过软分类获得的比例图像结果中。基于上述情况，本文通过统计分析地形数据与土地覆盖类型之间的关系，并将地形数据作为辅助信息用于子像元交换算法的初始化过程，来代替原算法中的随机初始化过程，以研究将地形作为辅助信息进行基于子像元交换算法的超分辨率制图的有效性和可行性。

1 研究区概况与数据源

研究区位于浙江省兰溪市东北部，研究区土地覆盖、高程和坡度如图1所示。

(a) 土地覆盖(b) 高程 (c) 坡度

图1研究区土地覆盖、高程和坡度

Fig.1Landcover,elevationandslopeofstudyarea

研究区面积约为420 hm2，覆盖研究区的航空遥感影像包含512像元×512像元，空间分辨率为4 m。土地覆盖类型图来自于对航空遥感影像的目视解译。土地覆盖类型主要有耕地、林地、建设用地和水域。通过对1∶5万比例尺地形图数字化后建立不规则三角网(triangulated irregular network,TIN)获得DEM数据，然后利用DEM计算得到坡度数据(以上运算均在ArcGIS9.3的spatial analysis扩展模块中完成)。

在本文中，将由子像元制图的参考图像生成比例图像的过程称为退化，各像元的土地覆盖类型的比例是以参考图像为依据、按照S计算得到[4,17]。计算得到的比例图像，作为子像元交换算法的输入。以参考的土地覆盖类型图为基础，通过退化获得S=2,4,8,16，32这5组比例图像。硬分类是传统分类的结果，以比例图像为基础，将比例最大的土地覆盖类型作为相应像元所有子像元的类型。图2与图3分别为研究区不同土地覆盖类型的比例图像和相应硬分类结果的示例，图上所示为S=8时的结果。

(a) 耕地(b) 林地 (c) 建设用地 (d) 水域

图2比例图像示例

Fig.2Exampleoffractionimages

图3 硬分类结果示例Fig.3 Example of hard classifcation result

2 实验算法

2.1 子像元交换算法

子像元交换算法的基本思想是通过交换像元内子像元的空间位置使相邻子像元之间的空间相关性最大。在算法的开始，按照各类型的比例，随机赋予像元中各子像元的土地覆盖类型。在初始化完成后，各像元中子像元的空间位置可以改变，而土地覆盖类型的比例保持不变； 按照邻近子像元之间及像元的空间相关性最大化的原则，改变像元中子像元的空间位置。该算法包含3个步骤：

1)对于每个像元，按照与距离成反比的规则计算每个子像元的空间吸引力。假设pi,j为像元Pa,b的一个子像元，pk为pi,j相邻的子像元，pk属于Pa,b或其相邻的像元，则Pa,b的空间吸引力OPa,b为

(1)

式中：S为超分辨率制图因子，每个像元包含S2个子像元；Opi,j为第i行、第j列子像元的空间吸引力，以其邻域中的子像元为基础，可得

(2)

式中：N为邻域的子像元数；Z(pi,j,pk)为二值函数，若子像元pi,j与其邻域中的子像元pk的土地覆盖类型相同，则其值为1，否则为0；λk为权重系数，与距离有关，即

(3)

式中：a为基于对数的非线性系数；h(pi,j,pk)为子像元pi,j与其邻域中的子像元pk之间的距离，即

(4)

2)以当前像元的子像元空间位置为基础，逐像元进行评估。

3)在同一像元中，选择2个具有不同土地覆盖类型且空间吸引力最小的子像元，如果交换其位置后该像元的空间吸引力增加，则交换其位置； 否则不进行交换。

对上述3个步骤进行迭代，直到子像元的空间吸引力不再增加时，运算过程停止。

Atkinson[7-8]提出的子像元交换算法是针对2个土地覆盖类型的情况，而Thornton等[21]则将其扩展到可用于多个土地覆盖类型的超分辨率制图。

2.2 基于DEM的初始化过程

将从DEM获得的高程和坡度数据用于像元中子像元的初始化，以代替原算法中的随机初始化过程。在本研究中，还分析了只用高程数据、只用坡度数据和同时使用高程与坡度数据3种情况。初始化过程包括4个步骤：

1)根据DEM数据和对应的土地覆盖类型统计二者的关系。由于高程和坡度为连续值，因此需要将其转换为类型值，在ArcMap的spatial analysis扩展模块中应用natural break将高程和坡度分为15个子区间，各子区间的范围如表1所示。

表1 各子区间高程和坡度数值范围Tab.1 Data ranges of elevation and slope for partitioning

对分区后的高程和坡度数据与土地覆盖类型之间的关系进行统计，假设土地覆盖类型为i，高程或坡度的子区间为j，则它们之间的关系Cij为

(5)

式中nij为土地覆盖类型i落在第j个高程或坡度子区间的子像元个数。因此，c是一个二维或三维矩阵，一个维度是高程和(或)坡度，另一个维度为土地覆盖类型。图4示出通过参考图像数据和高程/坡度统计得到的cij。

(a) 土地覆盖类型与高程的关系 (b) 土地覆盖类型与坡度的关系

(c)各土地覆盖类型与高程和坡度的关系图4 土地覆盖类型与DEM信息的相互关系Fig.4 Interrelationships between land cover types and DEM data

从图4可以看出，不同的土地覆盖类型与高程、坡度子区间之间存在明显的差异。

（2）当时，为通项递减的正项级数．因为，所以当时，；当时，．由定理2可知，当时，级数收敛；当时，级数发散．

2)以高程或坡度数据为基础，获得像元中子像元各土地覆盖类型的置信值。假定pi,j是像元Pa,b中的一个子像元，根据pi,j的高程或坡度与土地覆盖类型之间的关系cij，可以知道对于土地覆盖类型l的置信值是cl,pi,j。

3)将置信值进行归一化，即

(6)

以DEM数据为基础，逐像元地进行初始化。

2.3 子像元交换算法参数

(a)N=1(b)N=2(c)N=4 (d)N=6(e)N=8

(S为4； 中间黑色点为子像元pi,j； 其他灰色点为其邻域的子像元； 粗线条为像元范围)

图5邻域大小定义的示例

Fig.5Illustrationofdefinitionofneighborhoodsize

2.4 算法比较与评价标准

分别采用调整的Kappa值、所需CPU时间、迭代次数和总交换次数4个指标来比较和评价算法的表现。

调整的Kappa值通过计算结果与参考图像的比较获得，用于评价超分辨率制图的精度。与Kappa值不同的是只统计混合像元，因此对亚像元制图精度的描述更敏感和更合理。所需CPU时间用于比较算法的计算效率，在相同计算环境下获得。迭代次数和总交换次数也用于衡量算法的表现和效率，如在相同参数下，执行的迭代次数多，则所需的CPU时间就长； 同样，如果总交换次数少而制图精度高，则说明初始化过程的效果好、计算效率高。

在本研究中，各算法是由MATLAB的脚本m语言实现的，为了使计算结果更为可靠和更具可比性，对每组参数均重复5次，然后进行平均。

3 实验结果与分析

表2示出本文提出的初始化算法与原算法中随机初始化算法的超分辨率制图精度与计算所需的CPU时间。

表2 各初始化算法的调整Kappa值和计算所需CPU时间Tab.2 Adjusted Kappa values derived by different initializing algorithms and their consuming CPU time

①HC为硬分类； ②R为随机初始化； ③H为基于高程的初始化； ④S为基于坡度的初始化； ⑤H&S为基于高程和坡度的初始化。

从表2可知，在S≤8的情况下，本文算法制图精度要低于硬分类； 而在S>8时，其精度高于硬分类。除S=2外，本文算法的制图精度均好于随机初始化算法； 在S>8时，优势更为明显。在3种DEM辅助初始化算法中，S>8时高程数据辅助初始化算法的精度最高。与随机初始化算法相比，DEM辅助初始化所需的处理时间要长，而且随着S的加大，所需时间也明显增加。

图6是不同初始化算法的结果，图中线条表示参考的土地覆盖类型图中各土地覆盖类型边界。

(a) 随机初始化

(b) 基于高程辅助的初始化

图6-1 不同初始化算法5种S的制图结果Fig.6-1 Sub-pixel mapping results with 5 degraded S generated by different initializing algorithms

(c) 基于坡度辅助的初始化

(d) 基于高程和坡度辅助的初始化

图6-2 不同初始化算法5种S的制图结果Fig.6-2 Sub-pixel mapping results with 5 degraded S generated by different initializing algorithms

通过目视比较可以看出，基于DEM辅助初始化算法的制图结果要好于随机初始化算法。因此，应用基于DEM辅助的初始化算法可以为后续的子像元交换提供更好的基础。图7示出基于子像元交换算法的制图精度与硬分类和原算法制图精度的对比。在S>4时，本文算法能取得更好的制图精度。在DEM辅助的3种算法中，应用高程辅助算法的精度最好，而仅基于坡度辅助算法的精度最低。基于DEM信息辅助算法的不足在于，当S≤8时，对邻域大小这个参数的变化比较敏感。

图7 不同算法的子像元制图精度Fig.7 Sub-pixel mapping accuracy for different algorithms

图8示出不同参数情况下各算法子像元制图所需的CPU时间。

图8不同算法子像元制图所需的CPU时间

Fig.8CPUtimeinsub-pixelmappingcomputationwithdifferentalgorithms

尽管初始化过程所需的计算时间较多，但基于DEM辅助初始化的算法过程所需CPU时间要远少于原算法； 在基于DEM辅助的3种算法中，基于高程辅助算法所需的CPU时间最少，而基于坡度辅助算法所需的CPU时间最多。其原因是本文提出的算法无论是迭代次数还是总交换次数均明显少于原算法(图9和图10)。说明本文算法可以提高制图精度，同时提高运算效率。

图9不同算法子像元制图运算过程中的迭代次数

Fig.9Numberofiterationinsub-pixelmappingcomputationwithdifferentalgorithms

图10 不同算法子像元制图运算过程中的总交换次数Fig.10 Total swapping numbers in sub-pixel mapping computation with different algorithms

图11示出硬分类和4种子像元交换算法的超分辨率制图的最终结果。

(a) 硬分类

(b) 基于随机初始化的子像元交换算法

(c) 基于高程辅助初始化的子像元交换算法

(d) 基于坡度辅助初始化的子像元交换算法

图11-1不同算法子像元制图最终结果

Fig.11-1Sub-pixelmappingresultswith5degradedscalefactorsderivedbydifferentalgorithms

(e) 基于高程和坡度辅助初始化的子像元交换算法

图11-2不同算法子像元制图最终结果

Fig.11-2Sub-pixelmappingresultswith5degradedscalefactorsderivedbydifferentalgorithms

从图11可以看出，在S=32时，所有子像元交换算法的结果均好于硬分类，但本文算法获得的结果更加合理； 在大部分S值下，特别是当S较大时，本文提出的算法均好于原算法。

4 结论

尽管子像元交换算法是一种简单、有效的亚像元制图算法，但在超分辨率制图因子(S)较大时，仍然存在计算量大和制图精度较低的问题。本研究基于地形信息和土地覆盖类型之间的内在联系，提出了一种新的算法，将原算法中的随机初始化过程替换为基于DEM辅助的初始化过程。研究结果表明，本文提出的算法可以获得更好的制图结果，尤其是在S较大的情况下，优势更加明显； 而且，因为减少了迭代次数和总交换次数，其计算效率更高，所需的计算时间更少。在以地形作为辅助信息时，无论是制图精度还是计算效率，基于高程辅助初始化的效果都最好。

本文所提出算法的不足是，在应用本算法前，需要先获得地形信息与土地覆盖类型之间的关系，因此初始化过程和制图精度可能受其影响，对此需要进一步研究和评估。

志谢： 感谢美国密歇根州立大学全球变化与对地观测中心Qi Jiaguo教授在研究过程中给予的指导和帮助。

[1] Tatem A J,Lewis H G,Atkinson P M,et al.Super-resolution land cover pattern prediction using a Hopfield neural network[J].Remote Sensing of Environment,2002,79(1):1-14.

[2] Kasetkasem T,Arora M K,Varshney P K.Super-resolution land cover mapping using a Markov random field based approach[J].Remote Sensing of Environment,2005,96(3/4):302-314.

[3] Foody G M.Remote Sensing Image Analysis:Including the Spatial Domain[M].Norwell,MA:Kluwer,2004:37-49.

[4] 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.

[5] 任 武,葛 咏.遥感影像亚像元制图方法研究进展综述[J].遥感技术与应用,2011,26(1):33-44.

Ren W,Ge Y.Progress on sub-pixel mapping methods for remotely sensed images[J].Remote Sensing Technology and Application,2011,26(1):33-44.

[6] Atkinson P M.Mapping sub-pixel boundaries from remotely sensed images[M]//Kemp Z.Innovations in GIS IV.London:Taylor and Francis,1997:166-180.

[7] Atkinson P M.Super-resolution target mapping from soft-classified remotely sensed imagery[C]//Proceedings of the 5th International Conference on GeoComputation.Leeds,UK:University of Leeds,2001.

[8] Atkinson P M.Sub-pixel target mapping from soft-classified,remotely sensed imagery[J].Photogrammetric Engineering and Remote Sensing,2005,71(7):839-846.

[9] Foody G M.Sharpening fuzzy classification output to refine the representation of sub-pixel land cover distribution[J].International Journal of Remote Sensing,1998,19(13):2593-2599.

[10] Grossa H N,Schotta J R.Application of spectral mixture analysis and image fusion techniques for image sharpening[J].Remote Sensing of Environment,1998,63(2):85-94.

[11] Schneider W.Land use mapping with subpixel accuracy from Landsat TM image data[C]//Proceedings of 25th International Symposium,Remote Sensing and Global Environmental Change.Austria,Graz:[s.n.],1993,2:155-161.

[12] Tatem A J,Lewis H G,Atkinson P M,et al.Super-resolution target identification from remotely sensed images using a Hopfield neural network[J].IEEE Transaction on Geoscience and Remote Sensing,2001,39(4):781-796.

[13] Nguyen M Q,Atkinson P M,Lewis H G.Superresolution mapping using a Hopfield neural network with fused images[J].IEEE Transactions on Geoscience and Remote Sensing,2006,44(3):736-749.

[14] Ruiz C P,López F J A.Restoring SPOT images using PSF-derived deconvolution filters[J].International Journal of Remote Sensing,2002,23(12):2379-2391.

[15] Verhoeye J.De Wulf R.Land cover mapping at sub-pixel scales using linear optimization techniques[J].Remote Sensing of Environment,2002,79(1):96-104.

[16] Mertens K C,Verbeke L P C,Westra T,et al.Sub-pixel mapping and sub-pixel sharpening using neural network predicted wavelet coefficients[J].Remote Sensing of Environment,2004,91(2):225-236.

[17] Mertens K C,De Baets B,Verbeke L P C,et al.A sub-pixel mapping algorithm based on sub-pixel/pixel spatial attraction models[J].International Journal of Remote Sensing,2006,27(15):3293-3310.

[18] Shen Z Q,Qi J G,Wang K.Modification of pixel-swapping algorithm with initialization from a sub-pixel/pixel spatial attraction model[J].Photogrammetric Engineering and Remote Sensing,2009,75(5):557-567.

[19] 沈掌泉,虞舟鲁.基于单亲遗传算法和子像元/像元空间吸引模型的亚像元制图研究[J].遥感技术与应用,2015,30(1):129-134.

Shen Z Q,Yu Z L.Sub-pixel mapping with partheno-genetic algorithm and sub-pixel/pixel spatial attraction model[J].Remote Sensing Technology and Application,2015,30(1):129-134.

[20] Muslim A M,Foody G M,Atkinson P M.Localized soft classification for super-resolution mapping of the shoreline[J].International Journal of Remote Sensing,2006,27(11):2271-2285.

[21] Thornton M W,Atkinson P M,Holland D A.Sub-pixel mapping of rural land cover objects from fine spatial resolution satellite sensor imagery using super-resolution pixel-swapping[J].International Journal of Remote Sensing,2006,27(3):473-491.

[22] Nguyen M Q,Atkinson P M,Lewis H G.Superresolution mapping using a Hopfield neural network with LiDAR data[J].IEEE Geoscience and Remote Sensing Letters,2005,2(3):366-370.

Sub-pixelmappingoflandcoverusingsub-pixelswappingalgorithmandtopographicdata

YU Zhoulu1， WANG Wenchao2， RONG Yi2， SHEN Zhangquan1

(1.CollegeofEnvironmentalandResourceSciences，ZhejiangUniversity，Hangzhou310058，China；2.HangzhouFederationofLandPlanningandDesignConsultingCo.Ltd.，Hangzhou310030，China)

When remote sensing images are used to provide information for land cover mapping, it is negatively affected by the occurrence of mixed pixels in the remote sensing images, particularly in the case of coarse spatial resolutions. Soft classification and super-resolution(sub-pixel) mapping techniques can solve this kind of problems. The pixel-swapping(PS)algorithm is a simple and efficient technique for sub-pixel mapping. However, its computation is inefficient and yields poor mapping accuracy when the super-resolution scale factor (S)is large. A possible reason for this is that it only relies on the information from the fraction images. In this study, the digital elevation model(DEM) and their derivative data are employed as supplementary information for the PS algorithm so as to improve super-resolution mapping accuracy. Some conclusions have been reached: ① The sub-pixel mapping accuracy could be improved with the assistance of the DEM even if the scale factor is large； ② The mapping accuracy by incorporating both elevation and slope information is better than that of using elevation or slope data alone; ③ Mapping accuracy is less sensitive to the number of neighbors when scale factor is large；④ The computing efficiency is improved when incorporating DEM in pixel-swapping. Thus, it is feasible to use DEM as supplemental information for sub-pixel mapping.

sub-pixel mapping； super-resolution mapping； sub- pixel swapping algorithm； remote sensing； digital elevation model(DEM)； land cover mapping

10.6046/gtzyyg.2017.04.14

虞舟鲁,王文超,戎奕,等.基于子像元交换算法和地形数据的土地覆盖亚像元制图研究[J].国土资源遥感,2017,29(4)：88-97.(Yu Z L,Wang W C,Rong Y,et al.Sub-pixel mapping of land cover using sub-pixel swapping algorithm and topographic data[J].Remote Sensing for Land and Resources,2017,29(4)：88-97.)

TP 751.1

A

1001-070X(2017)04-0088-10

2016-05-05；

2016-06-28

国家科技支撑计划项目“旱区多遥感平台农田信息精准获取技术集成与服务”(编号： 2012BAH29B04)资助。

虞舟鲁(1986-),男,研究实习员,主要从事农业遥感与信息技术应用、土地利用等方面的研究。Email： yuzl@zju.edu.cn。

沈掌泉(1969-),男,副教授,主要从事农业遥感与信息技术、计算机应用等方面的研究。Email： zhqshen@zju.edu.cn。

(责任编辑：张仙)