基于数学重采样的遥感土地利用调查精度验证方法

基于数学重采样的遥感土地利用调查精度验证方法

1遥感土地利用覆盖的应用通过研究遥感土地利用，可以为生态环境动态监测系统提供强有力的数据，但验证数据非常重要，被忽视。很多情况下,由遥感解译得到的数据,往往只经过简单的野外实地验证,进行一些定性的评价,就交付使用或入库,缺乏较严格认真的精度评估。我国土地资源遥感动态监测正处于试点阶段,对生态环境基础数据需求大而迫切。应当看到,现今信息技术和网络通讯高度发达,保证生态环境基础数据的质量,比过去任何时候都重要,广泛共享的、高质量的基础数据,是信息时代生态环境动态监测和建设的需要,直接关系到生态恢复建设的质量,意义重大。评价遥感土地利用成果的精度仍然是整个土地资源遥感动态监测工作的核心。目前,在遥感土地利用(覆盖)调查中应用最为广泛的是TM(ETM+)数据。土地利用调查成果数据的几何精度评价,通常可采用理论分析的方法,即采用误差传播理论进行分析。用该方法进行几何精度评价必须获得数据取得过程中各环节所产生的误差。另一种方法是用遥感图像与调查数据叠合并计算叠合的吻合度。这种方法随着GIS与遥感集成化程度和水平的提高将会受到重视。第三种方法是先从微观角度入手,计算图幅的中误差,最后推出总面积的中误差;最后从宏观角度出发,计算评价数据的整体几何精度、总面积的相对精度等。本文结合“九五”科技攻关“西部大开发”专项“西北地区遥感本底调查及动态监测系统开发”的研究,以定西遥感土地利用调查为例,介绍一种利用GIS叠置的方法,对项目区进行再次解译,对成果进行叠置分层分析计算,可生成完整的数据库精度评价结果或分类系统不同级别的精度评价结果。2遥感本底调查及动态监测系统开发甘肃省生态环境遥感本底调查及动态监测系统开发是国家“九五”科技攻关“西部大开发”专项“西北地区遥感本底调查及动态监测系统开发”的子专题之一。数据源为TM图像(分辨率为30m),内容包括生态环境遥感本底调查、生态环境评价、生态环境治理规划、生态环境数据库的建立等。2.1地貌单元及其古地貌试验区位于甘肃省中部黄土高原、甘南高原、陇南山地的交汇地带,介于东经103°44′～105°38′,北纬34°50′～36°02′之间,包括定西、陇西、渭源和通渭四县。海拔大部分在1900～3000m之间。全区大的地貌单元属于黄土高原丘陵山地类型,是在新第三系、第四系早期古地形基础上,覆盖了厚度不等的风蚀黄土,经水流侵蚀切割而成。地貌差异大,梁峁起伏,沟壑纵横,地形破碎是区内最显著的特点。2.2治理生物措施相结合,进行小流域综合治理近几年来,区内各县以小流域为单元,以“水保立县”为发展战略,山、水、田、林统一规划,蓄、引、灌、排相结合,工程措施和生物措施相结合,大力发展小流域综合治理。对坡度大于25度的坡耕地进行退耕还林、还草,发展种草种树,进行生态环境重建;发展以“121雨水集流”工程为主的集雨节灌,进行土豆、中药材以及经济林种植,全面发展特色农业;同时,进行了旧城镇改造、小城镇和农村居民地建设;另外,宝兰二线、搀柳高速公路、国道省道改造和工矿区的建设等都大大地改变了定西地区的土地利用现状。3定性描述方向发展和定量方向发展方向遥感技术的发展带动了环境科学也不断向深度和广度方向发展,并且逐步由定性描述向定量方向发展。Landsat-5TM影像能客观、直接地提供相关的地理环境信息,因此调查中选择TM作为主要信息源。3.1土地利用分类系统本专题所采用的总的土地利用分类系统为国家“863-308”西部金睛行动《中国生态环境本底调查分类系统》(土地利用/土地覆盖),在没有改变总的分类系统的同时,结合试验区的特点进行了增减,形成了既统一又具有本地区特点的新的土地利用分类系统(见表1)。3.2辅助数据评价屏幕目视解译是解译人员以自身对该地区地域特征的了解为基础,根据影像的纹理、色调、形态及其在该地区的分布规律,对应一定的分类系统及技术要求,对地物判读的过程;解译过程是一个反复的过程,首先确定分类系统的第一级单元;然后进一步细分一级单元。在这一步需要对辅助数据进行认真的评价。细分大的土地覆盖单元要反复进行直到超过最小制图单位的所有土地覆盖单元被描述和确定为止。然后,对这些初步确定的单元进行野外调查,解决剩余的解译问题和解译综合验证。原始数据收集和验证的土地覆盖类型在一个点是确定的。与收集原始数据相比,解译是反向进行的。解译采样点的类型别时,不考虑周围环境的影响。忽视最小制图单元的要求将导致完全不同的解译结果。显示时应把比例尺固定在小于1∶40000的某一个值。当缩小图像时,能够看到图像的全貌。没有采样点的区域,可以显示原先的解译结果,这样有助于解译者熟悉已经采用的解译方法。在源数据收集时,详细记录野外调查结果,以节约经费。然后检查验证数据的拓扑错误。记录解译过程中的错误。最后建立一个所有土地覆盖单元的汇总表。3.3不同产土类型土地利用结构及其变化试验区四县土地总面积为1103014hm2,其中耕地面积为562936.94hm2,占土地总面积的54.4%;林地(包括落叶林、常绿林、防护林、灌木林等)66344.17hm2,占总面积的6.01%;草地108647.85hm2,占总面积的9.85%;公路、铁路及厂矿用地为282.12hm2,占0.025%;城镇、农村建筑用地为17208.15hm2,占总面积的1.56%;荒地为185981.37hm2,占16.86%;园地为118.8hm2,占0.01%;水体面积为488.37hm2,占0.04%;难利用地和其它用地为123324.56hm2,占总面积的11.18%。土地利用结构之比例为:农∶林∶牧∶水∶其它=54∶6∶10∶0.03∶29.97。图1为试验区1990年土地利用结构比例图。经过十年的发展变化,项目区土地利用和土地覆盖发生了很大改变。其中耕地面积为600120hm2,占土地总面积的58.86%;林地(包括落叶林、常绿林、防护林、灌木林等)73797.63hm2,占总面积的6.69%;草地83778.10hm2,占总面积的7.6%;公路、铁路及厂矿用地为341.57hm2,占0.03%;城镇、农村建筑用地为19766.66hm2,占总面积的1.79%;荒地为149357.70hm2,占13.54%;园地为231.68hm2,占0.02%;水体面积为400.67hm2,占0.03%;难利用地和其它用地为125965.94hm2,占总面积的11.42%。土地利用结构之比例为:农∶林∶牧∶水∶其它=59∶7∶8∶0.03∶25.97。图2为试验区2000年土地利用结构比例图。5采样点的确定和验证精度验证是在对成果进行定性评价的基础上的定量分析,是确定土地覆盖单元类别和面积划定和描述的精度。验证时,采用相同的数据源按照同样的解译方法对采样点再次解译,然后对两次解译结果进行比较。比较结果提供一个混淆矩阵(图3所示),用以区分特定的土地覆盖类别的错分程度1。4.1采样点地理位置第一步,选定层并确定选定层的选择点数。第二步,确定采样点,如确定每个土地覆盖单元内采样点的地理位置。验证时将确定解译错的面积占每层总面积的比例。4.1.1解译误差的解译算法确定h层需要的采样点数按下列公式计算:nh=Ρh(1-Ρh)σ2hnh=Ph(1−Ph)σ2h(1)其中:nh是采样点数,Ph是估计的解译误差率,σh是h层可接受的解译误差的绝对标准差。通过对选定的单元进行两次解译可估计解译误差率。利用GIS的叠置方法,统计计算两次解译的代码相同的单元和不同的单元的总面积,后者除以样点解译区的总面积得解译误差率。在得到验证结果、知道错误率的基础上就可以计算标准差。必须估计每层的错误率和确定每层可接受的标准差。然后,就可计算分层的采样点数。若每层的面积比较小,应增加点的密度。但考虑到实际操作的可能性,每平方公里的点数须限制于一定的数量。因此,面积相对小的层的采样点的总数将减少,这将导致比较高的标准差。4.1.2采样点的选择和网格宽度的确定第一步,步面积是根据每层的总面积除以每层的采样点数计算的。土地覆盖单元的绘制和每个选定的土地覆盖单元的采样点数的计算是按累积方法进行的。每层土地覆盖单元的面积逐渐累加,然后,从随机采样的起点按步面积系统选择土地覆盖单元,同时计算每个土地覆盖单元的采样点数。图4说明了单元选择和采样点的确定。其中Ah表示层h的总面积;Ahi表示层h第i类的土地覆盖类别的面积;nh表示层h的随机采样点数;AhnhAhnh表示步面积。土地覆盖单元面积逐渐累计,如,单元面积累计从图幅的西北角开始,到西南角结束。起始点选择:随机选择0到步面积之间的随机值,垂直箭头表示随机采样点,每个单元当被采样过程选中时,有一个采样点。图4表明单元1和3有一个采样点,单元5有两个采样点,单元2和4没有采样点。这种方法保证了随机选择在空间上的均匀分布。当单元面积大于步面积时,总会被选中。采样点数取决于单元的大小,因而选择是与面积成正比的。第二步,确定单元内随机采样点的位置。结合局部随机采样方法确定单元内的采样点。这种方法保证采样点不会集中在单元内的一个部位,而是在单元内最优分配。下面是采样时的3种不同情况:第一种情况如果给定的单元只有一个点,采用纯粹的随机采样。根据单元的坐标范围,确定能完全包含单元的最小矩形(边界矩形)。如果随机选中的点在单元内,那它就是有效的,否则,采样一直进行到一个选中的点有效为止。第二种情况如果有多个点,系统采样和局部采样要结合起来。第三种情况如果有n个点,边界矩形就要按平行于其最短边的直线分为n个同样大小的小矩形。每个小矩形内采一个点。如果没有一个点有效,这个过程就要再次进行。如果至少有一个点有效但不是所有的点有效,剩余的点为r(0<r<n)。为分散这些点,随机选择n个矩形中的r个矩形,然后确定一个样点。如果r个随机选择点都是有效的,就确定了所有的点。否则将重复以上步骤直到剩余r′个点(0<r′<r)。如果采样多于5个(包括5个),给定单元的矩形就要用一个指向北的规则的正方形网格覆盖。随机选择的网格交点决定网格的位置。网格的宽度由以下确定:一是单元上绘制的采样单元的数目n;二是从北到南和从东到西单元的最大范围。即矩形边的长度和宽度。为确定网格宽度d,n应该大于或等于覆盖单元的网格的最大数目,这是d的最小值。采样方法的影响随d值增大而增大,因此,d值应该尽可能地接近最小值。如果R和H分别代表矩形的长和宽,则:s=([ΔR]d+2)([ΔΗ]d+2)s=([ΔR]d+2)([ΔH]d+2)(2)d未知时上面方程不能解。网格宽度可由另一个方程确定,即:n=([ΔR]d+2)([ΔΗ]d+2)n=([ΔR]d+2)([ΔH]d+2)(3)然而多数情况下,实际得到的宽度值太大。因此,上面的方程要用以下的方程替代。n=([ΔR]d+1)([ΔΗ]d+1)n=([ΔR]d+1)([ΔH]d+1)(4)从而得到以下结果:d=√ΔRΔΗn-1+c2+c和c=ΔR+ΔΗ2(n-1)如果m是覆盖单元的矩形数目,那么在多数情况下,不等式n<m能满足。m个矩形确定之后,以后的步骤与第二种情况相同。n=a×m+r(r<m),定义了一个总数。m个矩形中的每一个只能采一个点。由于无效的采样过程,这将剩余r′个点。如果r′<m,这个过程以r′代替n重复进行。如果r′>m,剩余的r′个点按照第二种情况的方法采样。对于一个给定的单元n<m,公式(4)的使用是可能的,很明显开始的时候a=0,r′=n,即:开始时总样点数等于剩余的样点数。完成随机采样后,就确定了每个点的地理位置。现在,就能解译点的土地覆盖类型。4.2熟悉解译方法的人完成对随机采样点的解译,由没有参与收集原始数据但熟悉解译方法的人完成。使用相同的数据,对随机采样点周围1km的范围进行解译。解译员不知道随机采样点的位置,但同样参考一些专题图件数据。4.3确定主观上的确定估计值完成验证数据的解译后,对验证结果与采样点的原先解译结果进行叠置分析,建立混淆矩阵。正确分类的点数除以每层的总点数得到每层分类精度的估计值。ˆΡh˜h=Xh˜hnh(5)ˆΡh为样点的估计值,在数据收集过程和验证过程中,属于h层且是分类系统中的同一类的比值;Xh为h层正确分类的点数;nh为h层的总点数。因为样点数与面积成正比,每层正确解译区域的面积等于正确分类的点数乘以步面积就可得到。可根据分类系统的类进行分层。在这种情况下,分类系统中每类和总面积的精度就能够得到估计。多层精度的估计值与面积成比例,即ˆΡ=∑hAhAˆΡh˜h(6)ˆΡh为h层正确分类面积比例的估计值;Ah为h层的总面积;A为调查区域的总面积。根据分类系统进行分层时,汇总分类系统中不同级别的验证结果是很有意义的。更高级的正确分类面积的估计值按两步进行计算。第一步,将每层中的更高级别的正确解译结果点的估计值加起来,即:ˆΡh˜Η=∑˜h∈˜ΗXh˜hnh(7)其中:ˆΡh为在汇总的˜Η级h层正确分类面积的估计值;Xh˜h为在汇总的˜Η级h层正确分类点数的估计值。然后,对H级的样方不同层的估计值进行汇总,这要根据每层面积的大小按比例进行,即:ˆΡΗ˜Η=∑h∈ΗAhAΗˆΡh˜Η(8)ˆΡH为汇总的H级正确分类面积的估计值;AH为汇总的H级的所有层的总面积。最后,计算每一个正确分类区域估计值的标准差。每层的标准差可根据二项式分布计算。估计标准差的公式没有考虑选择面积单元的技术和在单元内确定多个样点的技术。因为这些技术通常增加了估计的精度,从而保证标准差比计算值要小。公式(5)、(6)和(8)可用于估计误差方差:Vˆar(ˆΡh˜h)=ˆΡh˜h(1-ˆΡh˜h)nh-1(9)这是一层的正确分类点比例估计值的方差。而每层面积总精度的方差为:Vˆar(ˆΡ)=∑h(A2hAˆΡh˜h(1-ˆΡh˜h)nh-1(10)汇总的级别的方差估计值按同样的方法计算:Vˆar(ˆΡΗ˜Η)=∑h∈Η(A2hAΗˆΡh˜Η(1-ˆΡh˜