泛第三极典型地区土地利用和水保措施解译精度分析

杨力华， 杨勤科， 朱梦阳， 王春梅， 庞国伟

(1.中国科学院 水利部 水土保持研究所， 陕西 杨陵 712100；2.西北大学 城市与环境学院， 陕西 西安 710127； 3.中国科学院大学， 北京 101400)

泛第三极地区是以青藏高原为主体的第三极向西扩展，包括青藏高原、帕米尔高原、兴都库什、天山、伊朗高原、高加索、喀尔巴阡等山脉，是全球环境与气候变化的敏感区和功能脆弱区[1]，土壤侵蚀是该区最主要的环境问题之一，且呈现明显的增加趋势[2-5]，但是泛第三极地区土壤侵蚀实测和制图数据匮乏。为了制定科学的土壤侵蚀防治对策，须开展土壤侵蚀抽样调查，其中土地利用和水保措施是调查的主要内容。土壤侵蚀抽样调查的目的是基于可免费下载、高分辨率的遥感影像，结合实地调查，快速、准确的获取抽样调查单元的土地利用和水土保持措施信息[6]，支持抽样调查单元土壤流失速率的精确计算，从而使大区域土壤侵蚀评价由目前的潜在评价(如Borrelli[7];Teng[3])逐渐逼近实际情况。为此，对土地利用解译结果做出精度评价十分重要。对于土地利用解译结果的评价，通常有两种方式： ①以野外调查或者更高分辨率影像解译结果作为真实土地利用信息参照，通过计算混淆矩阵进行评价[8]。如Chen[9]以Google Earth影像和无人机航拍影像解译结果为参照，用空间数据二级抽样检验方法，对30m分辨率土地覆盖数据(陈军等)[10]进行了精度评价；新近发布的全球10 m分辨率土地利用数据，通过Google Earth影像作为参考，基于OLI数据解译样本点，采用混淆矩阵进行了精度评价[11-12]；李骜等[13]基于土壤侵蚀野外抽样调查，利用土地利用变化转移矩阵对河北怀来县资源三号卫星遥感普查所得土地利用类型进行评价。 ②与已有土地利用数据集比较一致性，如Hansen[14]用平均单类分类精度比较了IGBP-DISCover数据集(1 km分辨率)和UMD数据集(1 km分辨率)的精度。Fritz[15]利用数据之间的叠合度评价了GLC-2000(1 km分辨率)和MODIS数据集(1 km分辨率)的一致性。Yang[16]通过计算总体一致性和单类分类一致性评价了全球7种土地覆盖数据集在中国的精度。上述精度评价所用参考影像多是高分辨率影像(如无人机航拍影像以及Google Earth影像)，其解译过程是人机交互目视解译进行的。然而，参考数据也会受到自身数据来源、分类方法和分类系统的限制，精度并不一定能满足验证需求。所以基于野外调查，探讨大范围区域(如泛第三极)、高分辨率遥感数据的土地利用解译结果的精度，是十分必要的。本研究中西藏样区以无人机航拍方式进行调查，泰国和巴基斯坦样区以现场调绘结果为参照，对Google Earth下载的高分辨率影像土地利用和水土保持措施解译结果做出精度评价，并通过与GLC 30 m[9]， GLC 10 m[12]的比较进行解译结果的区域代表性分析。本研究结果对于泛第三极抽样调查单元土地利用信息的精确提取，以及土地利用结构与变化特征分析具有重要意义。

1 研究区概况

本研究在土地利用研究比较少，解译难度比较大的西藏、泛第三极暖湿区的代表性区域泰国以及缺乏土壤侵蚀观测数据的巴基斯坦选择典型样区进行调查。西藏样区的调查在海拔相对较低、土地利用较为丰富，土壤水蚀和风蚀均比较严重的“一江两河”(西藏东南部)地区进行。该区面积6.65×104km2，海拔3 200～4 900 m，气候类型为高原季风温带半干旱性气候，年降水量约300 mm(西部)到500 mm(东部)[17]。该区地形较为破碎，主要以高原、山地、河谷为主。研究区耕地面积相对较多，自然植被以林灌草混交为主。“一江两河”区域河谷地带沟蚀剧烈，宽谷地带覆沙坡面的片状侵蚀发育明显[2]。泰国位于亚洲南部，地处热带季风气候区域。本次调查地区为泰国北部清莱地区，面积约1.20×104km2，是全国地势最高的山区，平均海拔约1 600 m。该区土地利用类型比较复杂，水田、旱地和园地均有分布。此外，泰国北部地区可利用耕地面积有限，毁林开荒式轮作和不合理的耕作措施使得该区成为泰国土壤侵蚀最剧烈的地区[4]。巴基斯坦位于南亚次大陆西北部，属于热带气候，降雨空间分布不均，地形多样。本次调查区波特瓦尔高原位于巴基斯坦北部，主要由平原区和沟壑区组成，坡度范围为8%～40%，面积约220 km2，是巴基斯坦重要的粮食生产地[18]。土地利用类型比较单一，多为耕地、草地和林地。年降雨量为500～1 400 mm，约50%～80%的降雨分布在5—10月，降雨侵蚀力较大，侵蚀多为片蚀和沟蚀[5]，是巴基斯坦水土流失最严重的地区。

2 材料与方法

2.1 数据及处理

本次调查目的是对室内解译的土地利用和水保措施信息进行现场复查和精度评价，后续相关分析中将对调查区域土地利用和水保措施的空间变异做出评估。因此调查单元沿着环境梯度变化比较大的断面布设。

2.1.1 抽样单元的确定 土壤侵蚀抽样调查的单元，是一个具有固定位置，有一定面积，可表现土壤侵蚀基本特征，在土壤侵蚀因子(特别是水土保持措施)方面表现出空间异质性的微小地理区域[6]。在丘陵区和山区调查单元为面积约0.2～3 km2的小流域，地形平坦区域则为 1 km×1 km的矩形区域。对于抽样单元的布设(图1)有3个方面的考虑： ①在研究区内选择沿某环境要素梯度变化大的方向，并考虑可通达性布设断面，在西藏主要考虑了海拔和植被类型的变异，在泰国北部主要考虑了地形的变异，在巴基斯坦主要考虑了降水和地形的变异； ②根据野外工作的难易程度布设样点，西藏相对便于纵深作业所以样点间距50 km，但是每个点作业难度比较大，因此布设了11个样点，完成了9个小流域的无人机航拍，泰国北部地域狭窄且泰国政府有一些限制，因此样点之间间距大约10 km，但是每个点上工作难度较小，所以完成了18个样点的野外调查，巴基斯坦样区沿伊斯兰堡东北—西南方向每隔10 km布设一个样点，共布设15个； ③样区布设考虑其可通达性。

图1 研究样区位置与样点布设示意

2.1.2 抽样单元室内土地利用解译 从Google Earth下载42个抽样单元的近期高分辨率(0.46～0.56 m)遥感影像，影像主要来源于2017—2018年Worldview卫星和航拍数据的整合，其中Worldview卫星影像是全色波段与多光谱波段的融合。根据“泛第三极土壤侵蚀抽样调查单元遥感解译工作大纲”的统一规定，在ArcGIS中调用已准备好的解译模板，利用人机交互目视解译方法，完成对土地利用和水保措施的提取，结果保存为42个shape文件。土地利用和水保措施分类采用第一次全国水利普查水土保持情况普查中应用的分类体系，参考2017年发布的《土地利用现状分类(GB/T21010-2007)》修改，将12个一级类简化为8个一级类：耕地、园地、林地、草地、建设用地、交通运输用地、水域及水利设施用地、其它土地，并对部分二级类进行归并和调整。关于解译技术流程，详见参考文献[6]。

2.2 野外调查

2.2.1 西藏东南部野外调查 野外调查于2018年8月12—27日进行，考察路线沿藏东地区的“一江两河”及雅鲁藏布江下游的主要支流尼洋曲。经过对当地水保局以及水文局的访问，初步了解西藏样区的土地利用、土壤侵蚀及其治理情况。然后利用微型无人机进行低空航拍摄影，同步对调查区代表性土地利用和水保措施进行现场识别，并用奥维互动地图和GPS工具箱拍照，整理现场代表性土地利用类型。无人机飞行航高均控制在300 m以内，使用Pix4dmapper进行单张航片拼接，获取的数字正射影像(DOM)分辨率均在6 cm以内。根据野外实地同步调查的数据为参考，使用ArcGIS软件解译无人机获取的数字正射影像(DOM)，以其解译结果作为参照数据。

2.2.2 泰国北部野外调查 野外调查于2018年11月22—29日进行，考察地区主要在泰国清莱地区。出发前将所选择的调查单元室内初步解译结果打印，并将每个单元导出成KML文件、再导入手机版奥维互动地图；在现场一方面直接在纸质图纸上修正地块边界和相应的属性代码，同时用奥维互动地图定点拍照记录代表性地物；参考现场调绘结果，使用ArcGIS软件对室内初步解译的底图进行修正，以其作为参照数据对初步解译结果进行评价。

2.2.3 巴基斯坦野外调查 野外调查于2019年4月15—24日进行，考察团沿着穆里—伊斯兰堡—拉瓦尔品第—古杰尔汗—杰格瓦尔—保恩路线进行调查。野外调查前对所选小流域的解译结果进行打印出图。野外过程中，一组借助奥维互动地图软件，实时定位并进行拍照记录；另一组根据小流域解译纸质图进行现场勾绘修改。野外校核结束后，总结室内解译存在的问题，使用ArcGIS软件对15个小流域的原始解译结果进行修改整理并作为参照数据评价室内初步解译结果。

2.3 质量评价方法

地理数据的质量评价，通常从空间准确性、专题准确性和时间准确性3个方面考虑[19]。空间准确性指图斑定位，专题准确性指土地利用和水保措施判读的定性是否准确。时间方面，由于本次室内解译使用了2017—2018年的影像，野外验证在2018年夏季和2019年春季，可以认为时相是一致的，所以时间准确性方面暂不做考虑。本研究具体从以下三方面进行高分影像解译结果精度评价。

(1) 空间准确性—混淆矩阵分析。将室内解译矢量图和参照图(修改后的、准确的图)均转换为栅格图层，用相对应的两个图层的所有栅格数作为样本点建立混淆矩阵，统计各调查单元的总体分类精度和kappa系数。

总体分类精度是正确分类的样本数与所有样本数的比值[8]，kappa系数[20]用来描述影像分类结果的一致性，其值越大精度越高(表1)。

表1 kappa系数分类评价标准

(2) 专题准确性—土地利用面积分类误差。遥感解译完成后，计算每一个图斑的面积，统计每一类水土保持措施的面积。以基于野外调查的解译结果作为正确值，将2种方式解译的土地利用图在ArcGIS中进行交集制表(tabulate intersection)统计，定量分析单类土地利用正确分类面积以及易错分的类型和面积。

(3) 土地利用结果区域代表性：本研究采用直方图相似度指数来表示室内解译结果与区域性土地覆盖数据产品(GLC30[9]和GLC10[12])在土地利用结构方面的相似性，用以表明解译结果的区域代表性。直方图相似度指数包括距离直方图相似度指数和相交直方图相似度指数，本文采用相交直方图相似度指数，即直方图的交(HI)[21]进行相似度计算。HI取值范围为0～1，HI值越大，说明两套数据集的相似度越高，因而说明抽样单元遥感解译的土地利用更具有区域代表性。

土地利用结构计算公式为：

(1)

式中：L为某种土地利用结构比例；Li为Uj中第i种土地利用类型的面积；i是3套数据集分类系统统一后的类别，为耕地、林地、草地、人造地表、水域、裸地6个一级类之一；Uj为第j个抽样单元的面积；k为抽样单元数量，k值分别取9，18，15。

3 结果与分析

3.1 西藏东南部样区的解译精度

3.1.1 土地利用分类精度分析

(1) 根据混淆矩阵统计西藏东南部9个抽样单元解译精度(表2)。总体分类精度平均值为83.48%，kappa系数平均值为0.74，8个抽样单元分类结果一致性较好，1个抽样单元(5号单元)解译精度欠佳。

表2 西藏东南部总体解译精度

(2) 统计西藏东南部土地利用面积转移比例(表3)表明，基于Google Earth影像解译的草地和耕地精度较高；水域解译正确率为80.42%，水域边界的不明显影响其解译正确率；裸地分类正确率较低，主要因为高寒地区裸地与草地在影像上难以区分，裸地易被错分成草地或者灌木林地；林地解译正确率偏低，主要因为藏东林芝地区林草混交，表现在林地中的二级类灌木林地与草地混分严重(如5号抽样单元)，小流域内种植大片沙蒿(小灌木林)，室内解译错将沙蒿分为草地(图2)；交通运输用地解译正确率最低，因为土地利用变化较快地区(如3,4号抽样单元)，在影像上未能反映。

表3 西藏东南部遥感解译对应野外调查的土地利用转移比例

图2 西藏样区5号抽样单元室内解译与野外调查对比

3.1.2 水土保持措施分类精度分析 野外调查(表4)表明，基于Google Earth影像解译和无人机影像解译的人工乔木林措施、水平阶措施面积一致性较好，四旁林措施面积相差较大。西藏东南部处于高寒地区，谷歌影像颜色多呈灰白色，建设用地周边零碎的四旁树较难划分，导致该种生物措施的面积解译出现较大偏差。总体上看，水土保持措施的解译尚比较理想。

表4 西藏东南部样区水土保持措施解译面积比较 km2

3.2 泰国北部样区的解译精度

3.2.1 土地利用分类精度分析

(1) 泰国北部抽样单元解译精度评价结果(表5)表明，18个单元的平均分类精度为85.37%，平均Kappa值为0.75，17个抽样单元分类结果达合格标准，1个抽样单元(1号单元)解译不合格。

表5 泰国北部总体解译精度

(2) 泰国北部土地利用面积转移比例统计(表6)表明，耕地解译正确率最高，达92.61%；其次较好的是建设用地和水域，由于土地利用转移，出现小部分错分情况；园地和林地解译精度均达80%以上，其中园林间种、套种等立体种植模式降低了其解译精度，如图3中15号抽样单元野外调查为021园地(咖啡)和031林地(橡树)混种，室内解译结果统分为021(园地)；分类精度最低的是草地，主要由于泰国北部地区地物类型比较复杂，种植面积不集中，斑块细碎化导致草地的错分(如1,2,3,18抽样单元)

表6 泰国北部遥感解译对应野外调查的土地利用转移比例

图3 泰国样区15号抽样单元室内解译与野外调查对比

3.2.2 水土保持措施分类精度分析 泰国北部多采用人工种植经果林，兼顾经济效益和水土保持效益。表7表明，基于Google Earth影像和野外调查解译的四旁林措施和土坎水平梯田面积一致性相对较高，人工乔木林和经果林解译面积一致性较低，主要由于土地利用类型解译错误，造成对应生物措施的面积不一致(如林地错分为园地，导致室内解译的经果林面积偏大)。

表7 泰国北部样区水土保持措施解译面积比较 km2

3.3 巴基斯坦样区的解译精度

3.3.1 土地利用分类精度分析

(1) 巴基斯坦北部抽样单元解译精度评价结果(表8)表明，15个单元的平均分类精度为90.22%，平均kappa值为0.82。14个抽样单元分类结果一致性较好，1个抽样单元(12号单元)解译精度欠佳。

(2) 巴基斯坦北部土地利用面积转移比例统计(表9)表明，基于Google Earth影像解译的园地和水域定性和定位均正确；其次较好的是耕地和林地，解译正确率均为90%以上；建设用地的解译正确率为49.10%，主要原因是该区域受人为因素干扰大，存在许多开垦的却未投入利用的土地易被错分为耕地，如2，7，9(图4)，10号单元；草地的解译正确率偏低，因为该区域存在多处开垦的新造地裸地(如1号单元)，易被误判为草地，且居民地旁边的零星草地未被解译(如14号抽样单元)；与西藏样区相似，对于土地利用类型转移较快的交通运输用地，解译正确率最低。

表8 巴基斯坦总体解译精度

表9 巴基斯坦遥感解译对应野外调查的土地利用转移比例

图4 巴基斯坦9号抽样单元室内解译与野外调查对比

3.3.2 水土保持措施分类精度分析 巴基斯坦北部样区地物类型相对而言比较简单，水土保持措施判读精度总体较好。表10表明，两种方法解译的四旁林、经果林面积一致性很高，土坎水平梯田解译精度欠佳，主要因为波特瓦尔高原地形多变，土坎水平梯田较多且较为分散，由于影像拍摄角度等，一些沟底耕地水保工程措施难以识别。此外，该区草本植物植株较低，人工种草措施在影像上纹理不清楚，出现室内解译漏判情况。

3.4 遥感抽样法解译土地利用区域代表性分析

3.4.1 典型地区土地利用结构相似性 图5表明，西藏样区、泰国样区、巴基斯坦样区室内解译与GLC 30 m数据集[9]土地利用结构相比，HI值分别为66.05%，80.41%和79.44%；与GLC 10 m数据集[12]相比，HI值分别为71.78%，80.42%和92.96%。室内遥感解译结果与参考土地覆盖数据集具有很高的相似性，且更符合GLC 10 m数据集的分类，可见抽样解译的土地利用具有区域代表性。

表10 巴基斯坦样区水土保持措施解译面积比较 km2

图5 室内遥感解译和GLC 30 m，GLC 10 m土地利用结构图

3.4.2 典型地区土地利用结构差异性 对于西藏样区，室内解译结果与GLC 30 m[9]，GLC 10 m[12]参考数据集差异主要在耕地、林地、裸地，特别是裸地，相差20%以上，参考数据集所给裸地面积偏大，但是通过分析Google Earth影像以及其历史影像，参考数据集显示的许多土地并非裸地，且经过野外校验，大多裸地是没有明显绿色作物生长的耕地或者草地。对于泰国样区，差异主要在耕地和林地，统一后的分类系统耕地包括种植农作物的土地(水田、旱地、水浇地)和种植经济作物的土地(园地)，参考数据集耕地分类面积偏多、林地分类面积偏少，经过野外校验，泰国北部地区斑块细碎化严重影响了园地、林地的分类，造成参考数据集与室内解译在园林土地利用结构上有着不同程度的差异。对于巴基斯坦样区，差异主要在耕地和草地，地块破碎化以及人类大面积开垦土地均造成了室内解译结果和GLC 30 m[9]，GLC 10 m[12]土地利用结构的差异。

4 讨论与结论

4.1 讨 论

(1) 土地利用与水保措施判读精度。基于Google Earth影像解译西藏样区、泰国样区、巴基斯坦样区的土地利用总体分类精度分别为83.46%，85.37%和90.22%，其中西藏样区精度比中国第一次全国水利普查水土流失动态监测解译的精度(64%～91%[13])略高。但是一些类型如西藏样区、巴基斯坦样区灌木林地和泰国北部园地解译正确率相对较低，这一情况也存在于其他相关数据中[9,12](如GLC 30 m，GLC 10 m)。受到土地利用解译难点问题的影响，水土保持措施(特别是林草措施和比较小尺度的工程措施)的判读也存在精度偏低的问题。这些情况只能通过使用更为现势性的、更高分辨率的数据，同时辅助必要的野外工作来解决。

(2) 遥感解译的数据源。本文所用的Google Earth影像分辨率在0.46～0.56 m之间，与正在开展的全国水土流失动态监测与公告项目中采用的高分1号和资源3号数据相比[13]，其分辨率较高，对于水土保持措施的分辨能力较好。但是该数据上有些影像特征和土地利用类型依然较难辨识，且部分地区可免费下载的数据现势性较差或不一定有高分辨率影像(0.5 m左右)。所以遥感解译需与野外调查、微型无人机航拍相结合，并不断丰富解译标志库，才能快速精确的进行大区域土地利用调查。

4.2 结 论

(1) 基于Google Earth高分遥感影像(0.46～0.56 m)完成的抽样单元土地利用和水土保持措施的解译，可满足土壤侵蚀评价的要求，3个典型地区土地利用类型总体精度达80%以上，平均kappa系数在0.7以上。通过对比数据(GLC 30 m，GLC 10 m)比较不同数据集土地利用结构一致性，说明抽样单元解译结果在土地利用结构方面具有良好的区域代表性，可支持区域土壤侵蚀调查和制图。

(2) 基于Google Earth高分遥感影像(0.46～0.56 m)在抽样单元上提取空间尺度比较小的水土保持措施是可行的。本研究通过对比野外调查和基于高分辨率遥感影像解译水土保持措施面积，表明在西藏、泰国和巴基斯坦样区遥感解译的生物措施和工程措施总体与野外调查相符，可为土壤流失速率计算中的生物措施因子和工程措施因子提供数据支撑。