我要投稿
投稿指南
RSS订阅近年，我国成功研制出全球首套30m分辨率地表覆盖数据产品（GlobalLand30），受到国际社会的广泛关注和高度评价[1]。为了推动其广泛应用和细化完善，须对其进行精度评价和全球验证。严格的高分辨率制图精度评估的过程包括抽取样本、获取样本的真值和精度评价[2]。因此，GlobalLand30数据抽样是对该数据验证的基础环节，一组科学合理的样本能够用最少的成本得出GlobalLand30数据的真正精度，而不合理的抽样往往会导致GlobalLand30的真正精度不符（或高或低）。所以，根据GlobalLand30数据特点设计科学合理的抽样方法是其验证的重要基础。
GlobalLand30抽样作为空间抽样的一种特殊形式，它与常规空间抽样有着共同的目标，即用较少的样本量获得较高精度的统计推断[3]。主要解决一下问题：（1）确定样本量大小；（2）样本如何布设 [4]。
然而，GlobalLand30作为地表覆盖数据，又区别于其他空间数据，它是将地球表面各种物质类型及其自然属性与特征总体归纳为特有的分类体系，共计10大地类：耕地、林地、草地、灌木、湿地、水体、苔原、人造覆盖、裸地、永久冰雪[5]。各个地类以图斑的形式相互连接，整张地图呈现为不同地类的图斑之间的拼接。不同的地区不同的地类之间，图斑的空间分布不同，呈现出来极高的空间异质性，例如在利比亚占地一半以上的都是大片的均质分布的裸地，只有在北部沿海等少数地区是异质程度高的草地等地类；在西班牙虽然耕地占有面积也在50%，但是除了在阿拉贡、莱昂、拉曼达、安达卢西亚几个地区内是以耕地为主的较为均质的空间分布，其他地区都是以林、灌、耕等为主的破碎度极高的空间分布形式。这一重大差异也决定了GlobalLand30数据抽样在确定样本量与布设样本时需要分析各地区各地类之间的分布差异，参考其他地表覆盖数据抽样[6]，类似于大面积沙漠这样的连续均质区域需要的样本是大大少于多种地类马赛克式分布的异质程度高的地区。而且，GlobalLand30数据在不同的地区中存在大量的稀有地类，这些稀有类以小图斑的形式分布，在抽样时需要适量考虑这些稀有地类的样本分布是否充足有效。
以西班牙为例，图1（a）所示西班牙境内在2010年的地表覆盖分布，表1所示各地类的占有面积。境内以耕地与林地为主（占69%），而水体和湿地的占有面积特别稀少（分别是0.7%和0.03%）。通常按照各地类的面积作为计算样本量的标准，那么，湿地与水体的样本数就会过于稀少。例如在总量3000时，湿地与水体分配的样本才为1和22。一般情况下，通过增加总样本量，或由专家人为调节的方式来增加湿地与水体的样本数量，这从某种程度上来说，是对成本的增加。此外，不同的地类也呈现不同的分布方式，例如林地、人造覆盖是以破碎图斑的形式分布，而耕地则在不同的地区呈现不同的分本方式。在西北部，耕地呈现极度破碎的状态，以图1（b）为例，在这一地区，耕地面积只占1%，但却存在耕地总图斑数的25%；而西班牙中东部（阿拉贡、莱昂、拉曼达、安达卢西亚）的地区集中了10%的耕地，呈现大面积连续的均质状态分布。在布设耕地样本时，应兼顾均质大面积地区（中东部）与异质破碎地区（西北部）之间的平衡。
|
|
耕 |
林 |
草 |
灌 |
湿地 |
水 |
人造 |
裸 |
|
面积比 |
51% |
18% |
5% |
20% |
0.03% |
0.7% |
2% |
2% |
至今全球共生产了5套全球地表覆盖产品，其中除了UMD产品，其他数据都进行了严格的精度评价[6-8]。多数采用了分层抽样，而其中DISCover、MOD12、GlobCOVER产品采用了完全随机的抽样方法，USGS产品采用两级整群抽样，都没有考虑各地区的空间异质分布。GLC2000因重点考虑耕地、林地等地类，抽样时引入了shannon多样性指数，初步考虑了地表覆盖各类型分布的空间异质性，但这是针对1km分辨率的地表覆盖制图的抽样方法。而常规空间抽样中[4]，“三明治”空间抽样模型将空间异质性作为布设因素之一，但因需要专家对每一地区都进行分区或分类，这对于全球地表覆盖数据抽样而言存在困难。
针对这一问题，本文针对地表覆盖数据空间异质问题，提出了引入景观形状指数（LSI）对地表覆盖数据进行抽样，可以从计算样本量与布设样本两个方面进行阐述。
2 引入景观形状指数
在上述地表覆盖数据抽样中，GLC2000抽样时引入的shannon多样性指数是景观生态学中的概念，是一种分析景观格局的方法，称为景观指数法。景观指数是指用来量化分析景观格局的指标，度量单元的类型、数量、形状、空间分布、复杂程度等[9-11]。近年来，景观指数的应用分析越来越广，主要应用在两个领域，一是用于分析评价土地利用/地表覆盖的格局和变化分析，另一个就是景观生态功能分析。常用指数较多[12]，可以分为描述形状的景观指数（形状指数）与描述斑块镶嵌的景观指数（丰富度、优势度等）。根据地表覆盖数据验证的目标，即验证图斑边界是否正确，验证边界内的分类是否正确，笔者认为地表覆盖验证与图斑的内部、边界相关，因此，选择景观形状指数进行地表覆盖数据格局分析。
景观形状指数是通过计算形状与相同面积的圆或正方形之间的偏离程度来测算其复杂程度，其表达式如下[10]：
式（1）
式（1）中，E为所有边界的总长度，A为总面积。在栅格数据中，E为处于某一地类图斑边界的像元总数，A为某一地类的像元总数。
笔者认为景观形状指数（LSI）能够综合反映地类的面积与异质程度。如下图2所示的地区，各地类的面积与总面积之比与LSI值统计如表2所示。根据LSI与面积占分析，耕地、林地、人造覆盖占有面积高的地类，LSI数值也高；稀有地类草地、湿地、水体、裸地的面积比较小，LSI数值相对低。此外，相比较耕地与林地，虽然耕地的占有面积较少，但是由于整体破碎度高，空间结构更复杂，所以LSI又高于林地。LSI值能够综合量化分析各地类的面积与整体分布复杂程度。因此，若将各地类LSI值作为各层样本量的分配标准，则可以降低分布广、连续性大的均质地类的样本量；又可以适度提升复杂的稀有地类的样本量。这样可以减少根据面积分配样本的人工干预。
|
|
面积比% |
LSI |
|
耕地 |
25.48 |
34.66 |
|
林地 |
48.89 |
25.11 |
|
草地 |
0.93 |
8.97 |
|
灌木 |
7.43 |
16.55 |
|
湿地 |
0.03 |
2.20 |
|
水体 |
0.18 |
7.84 |
|
人造覆盖 |
16.20 |
23.89 |
|
裸地 |
0.81 |
6.22 |
图3 2010年GlobalLand30在巴塞罗那及邻近地区人造覆盖格网及LSI值
3 引用LSI的地表覆盖数据抽样方法
引用LSI的地表覆盖数据抽样方法，采用分层抽样方法。首先根据某一区域的大小与空间异质程度计算该地区的样本总量，其次根据每一层的LSI值将样本总量分配到各层中。然后在每一层布设格网，根据格网内的LSI值选取有代表性的格网，再在每一格网内随机布设样本。
3.1 样本量计算方法
常规的样本量的计算是根据概率统计学推断的，例如TRASP[13]，只考虑了需要抽样的区域大小、均值等情况，没有考虑空间数据的异质分布。引用LSI的地表覆盖数据抽样方法在计算样本量时，首先计算目标区域内的各层LSI值（
）与LSI总值（
），然后根据区域内的面积与LSI总值计算样本总量，再根据每一层的LSI值分配这些样本。
计算某一层的LSI值（
）时，考虑到目标区域的大小，有时无法在一个图幅内全部包含，因此引入了虚拟图幅。假设该区域多个图幅覆盖，图幅边框与目标区域多边形存在两种拓扑关系，即“包含”与“相交”。 为“包含”关系的图幅数量为m，这些图幅在该层的LSI值（
）可以直接相加；假设为“相交”关系的图幅数量为n，这些图幅在该层的LSI值需与该区域相交部分的面积与图幅总面积之比相乘。LSI总值（
）则将各层在目标区域内的LSI总值相加得到，计算方法如式（2）所示，
式（2）
式（2）中，
为第k层在目标区域内的LSI总值，
为第k层在第i个图幅内的LSI值，m是与目标区域多边形呈“包含”关系的图幅数量，n是与目标区域多边形呈“相交”关系的图幅数量，
是与目标区域多边形“相交”的第j个图幅中相交部分的面积，
是与目标区域多边形“相交”的第j个图幅中的总面积。
在计算样本总量之前，选择一参照地区，求得该地区的LSI总值为
，总面积（或总像元数）为
，该地区的总样本量为
（计算方法可参照TRASP（童，2011）等方法或自行定义）。当计算目标区域的样本总量时，计算得目标区域的LSI总值为
，面积（或总像元数）为
，将目标区域与参照地区作对比，得出总样本量公式，如式（3）所示，
式（3）
求第k层的样本量
时，根据该层的LSI值进行计算，可直接这一层的LSI与总LSI值之比作为分层权重，公式为式（4）所示，
式（4）
3.2 基于格网的LSI值布设样本方法
得出每一层的样本量
，就在每一层布设样本了。布设方法如下：
1在目标区域内划分格网，要求格网数远远大于
，一般格网大小的选择会考虑数据的空间分辨率，GlobalLand30格网选择为1km或2km，格网的大小对结果影响不大；
2在该层内，求取每一格网内的LSI值；
3筛选LSI值，对于矢量数据而言，LSI均大于等于1，由于栅格数据造成边界像元数等于像元总数的情况，则存在LSI值小于1的情况。因LSI小于1的格网内，目标层的数量极少，不具有代表性，因此删除LSI小于某一阈值的格网。一般可选取0.8或0.9；
4将格网按照LSI值大小排序，将异质度高的格网与异质度低格网区分开；
5将全部格网从大到小等间隔划分为
个阶段，若存在余数，可归为第一阶段和最后一阶段；
6在每一阶段随机或等序选取1个格网，在每一格网内随机选取一个样本。
总体计算流程如下图4所示：
图4引用LSI的地表覆盖数据抽样方法总流程
4 实验分析
实验的目的是验证根据引用LSI的地表覆盖数据抽样方法选取的样本是否科学合理，即对于全体数据是否具有代表性。因此，实验思路是选取某一地区，将该地区的GlobalLand30数据与另一地表覆盖数据进行全部像素的对比。再用引用LSI的地表覆盖数据抽样方法对这一地区进行抽样，将样本与另一地表覆盖数据进行对比，得出样本对比结果。比较分析样本对比结果与全像素对比结果是否一致。若结果一致说明数据具有较高的代表性，反之说明结果较差。
实验区域选择西班牙，西班牙境内的地表覆盖分布情况较为复杂，以林地、耕地、人造覆盖为主，草地、湿地、水体等分布稀少，且在大部分地区呈现异质程度较高的情况。实验数据是Globalland30，与之对比的另一地表覆盖数据选取为Corine Land Cover。
在西班牙境内，Globalland30数据共计像元600,000,000多个，将其全部与Corine Land Cover对比，得到对比结果结果如下表3所示。
表3 西班牙境内GlobalLand30与Corine Land Cover全像素对比结果
|
|
耕 |
林 |
草 |
灌 |
水 |
湿 |
人造 |
裸 |
|
对比 结果 |
0.9297 |
0.9024 |
0.8031 |
0.8266 |
0.7211 |
0.6054 |
0.5045 |
0.8988 |
根据引用LSI的地表覆盖数据抽样方法，参考区域选取的是空间异质程度为西班牙平均水平的地区，选取参考区域内的像素总和
= 12101382，LSI总值
= 60，参照TRASP计算总样本量
，得出样本共计3044个，各层样本数量分布如下表4所示。观察可知，景观指数抽样修正了每一地类的样本量，大大降低了耕地这一呈现连续均质分布的类别的样本量，大大提升了草地、水、湿地、人造覆盖等异质程度高的样本量。
表4 根据引用LSI的地表覆盖数据抽样方法在西班牙计算所得样本量
|
|
耕 |
林 |
草 |
灌 |
水 |
湿 |
人造 |
裸 |
|
样本量 |
533 |
570 |
403 |
658 |
34 |
250 |
354 |
186 |
而将样本与Corine Land Cover对比之后，得出全像素对比与样本对比结果如下图5所示。对比全像素对比曲线（绿色）与引用LSI的地表覆盖数据抽样方法所得样本对比曲线（红色），观察可知，两曲线基本相吻合，说明样本与全像素的代表性较高。只有在人造覆盖、湿地这两类异质程度极高的地类存在部分差异。而观察蓝色曲线（使用随机抽样的方法所得样本对比结果），发现与全像素对比结果相差甚远，可知引用LSI的地表覆盖数据抽样方法远远优于随机抽样的方法。
图5 西班牙抽样对比、全像素对比结果比较
5 结论
本文针对Globalland30数据抽样时需要考虑空间分布差异的特殊性，提出了引用LSI的地表覆盖数据抽样方法。该方法引入景观形状指数LSI对每一层的数据进行量化分析，综合考虑该层的面积与空间异质程度，对样本量计算进行优化。在布设样本时，将格网与景观形状指数LSI相结合，将目标区域根据空间异质程度不同划分出不同的等级，在每一等级中抽取相应的样本，使样本分布兼顾不同区域的空间异质程度之间的差异。且选择西班牙地区进行方法验证实验，实验采用全像素对比与样本对比结果比较的方法，验证了引用LSI的地表覆盖数据抽样方法的适用性。实验结果分析表明，该方法提高了空间异质程度高的地区与地类的样本量，降低了连续均质分布的地区的样本量；且根据LSI值抽取的样本更能代表西班牙Globalland30的数据。
本文下一步工作将这一方法改编为基于网络的地表覆盖抽样，再结合网络多源异构的地表覆盖参考数据的自适应基础，发展Globalland30基于网络的验证方法。
参考文献
1. 陈军, 陈晋,廖安平等. 全球30m地表覆盖遥感制图的总体技术[J].测绘学报, 2014, 43(6), pp: 551-556.
2. CONGALTON, R.G. GREEN, K.. Assessing the Accuracy of Remotely Sensed Data: Principles and Practices[M] .Boca Raton, FL: CRC Press. 1999
3. 王劲峰,姜成晟等. 空间抽样与统计推断[M]. 北京.科学出版社. 2009.
4. 姜成晟,王劲峰等. 地理空间抽样理论研究综述[J].地理学报, 2009, 64(3), pp: 368-380.
5. 陈军, 陈晋, 宫鹏等. 全球地表覆盖高分辨率遥感制图[J]. 地理信息世界, 2011,4 (2), pp: 12-14.
6. Olofsson P. A global land-cover validation data set, part I: fundamental design principles[J]. International Journal of Remote Sensing. 2012.33(18), pp. 5768-5788.
7. SCEPAN, J. Thematic validation of Global Land Cover Data Sets –Procedures and Interpretation Methods[C].1999, Conference: Geoscience and Remote Sensing Symposium, 2001. IGARSS '01. IEEE 2001 International, Volume: 3
8. Philippe Mayaux. Validation of the Global Land Cover 2000 Map[J]. IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, 2006,44(7): 1728-1739
9. Evelyn Uuemaa, ülo Mander, Riho Marja，Trends in the use of landscape spatial metrics as landscape indicators: A review[J]. Ecological Indicators . 2013,28 : 100–106
10. 邬建国. 景观生态学-格局、过程、尺度与等级. 北京.高等教育出版社. 2001.
11. Virginia H. Dale, Keith L. Kline. Issues in using landscape indicators to assess land changes[J]. Ecological Indicators, 2013,28 :91–99
12. Forman, R.T.T., 1995. Land Mosaics: The Ecology of Landscapes and Regions[M]. Cambridge University Press.
13. Tong X., 2011. Designing a two-rank acceptance sampling plan for quality inspection of geospatial data products. Computers & Geosciences, 37(10), pp. 1570-1583.
