基于ENVI和ArcGis的云南省侵蚀模数图量算方法

刘 刚，刘阳容，杨林山，林 芸，房林东，田 帅，李绍平

(1.中国电建集团昆明勘测设计研究院有限公司，云南 昆明 650051；2.中国科学院西北生态环境资源研究院，甘肃 兰州 730000；3.成都市新都区城乡建设局，四川 成都 610500；4.湘西土家族苗族自治州水利水电勘测设计研究院，湖南 吉首 416000)

采用侵蚀模数图法来推求设计流域多年平均产沙量的方法已广泛应用于云南省中小流域水利水电工程泥沙计算中[1]。采用AutoCAD的传统云南省侵蚀模数图量算方法需以侵蚀模数图为底图，人工手绘侵蚀模数图图斑，并进行面积换算，再按照侵蚀类型取对应侵蚀模数，采用面积加权法计算流域平均侵蚀模数，进而推求设计流域多年平均产沙量。上述方法具有工作量大、耗时长、易出错、量算成果不能重复使用且与ArcGis水文分析成果不能有效衔接等缺点。

目前，在设计工作中使用的云南省侵蚀模数图为栅格数据(红、绿、蓝3个通道)，模数图坐标为大地坐标，有75幅之多。设计流域泥沙计算需要水文集水区划分成果(流域范围)作为计算边界，量算时需要单位换算，不便于直接量算。鉴于ArcGis水文集水区划分在设计中的广泛运用，且侵蚀模数图法在云南省使用广泛，迫切需要提出一种新的方法以实现云南省侵蚀模数图的快速量算，克服传统方法的诸多弊端，以提高设计人员的工作效率。

基于上述原因，本文提出了基于ENVI和ArcGis的云南省侵蚀模数图量算方法。本文利用ArcGis对云南省75幅侵蚀模数图进行处理，制作带投影信息的云南省侵蚀模数拼接栅格图。利用ENVI软件与ArcGis软件良好的兼容性，保留侵蚀模数拼接栅格图的投影信息，使用ENVI图像分类以及分类后处理功能对带投影信息的云南省侵蚀模数拼接栅格图进行矢量化[2]。侵蚀模数图矢量化后因带有投影信息，能够与基于ArcGis的水文集水区划分成果进行叠图分析，可推求流域(流域范围全在云南省)多年平均产沙量。

1 基础数据制作

基于ENVI和ArcGis的云南省侵蚀模数图量算方法，需要利用ArcGis与ENVI软件对设计中常用的侵蚀模数栅格图进行处理，得到云南省侵蚀模数矢量图。因设计项目往往分散，为方便工作避免因设计流域不同而重复基础数据制备环节，本文建议利用ArcGis对云南省75幅侵蚀模数图进行坐标校正、裁剪、拼接以及投影变化等处理，制作带投影信息的云南省侵蚀模数栅格拼接图，再利用ENVI图像分类、分类后处理等功能制作全云南省侵蚀模数矢量图。

a)数据制作步骤1。云南省侵蚀模数图共计75幅，每幅成图比例尺为1∶20 万。为保证制作的矢量图具有很高的空间位置精度，建议利用ArcGis软件的[Georeferencing]工具对每张侵蚀模数栅格图进行格网逐点校正，使用[Raster]-[Clip]工具裁去侵蚀模数栅格图多余边角，利用[Mosaic To New Raster]工具完成所有栅格数据的拼接。侵蚀模数栅格图坐标系为大地坐标，使用时需要转化为投影坐标，可使用[Project Raster]工具对拼接好的栅格数据进行投影转换，将大地坐标系转为云南省水文手册附图相同的投影坐标系[3]。

b)数据制作步骤2。ENVI软件与ArcGis软件具有良好的兼容性，可保留侵蚀模数拼接栅格图的投影信息。在ArcGis软件中，将云南省侵蚀模数拼接栅格图导出为ENVI可识别的栅格数据形式，并利用ENVI遥感软件对已经拼接好的云南省侵蚀模数图进行监督分类[4]。依据侵蚀模数图类别(无明显侵蚀、轻度侵蚀、中度侵蚀、强度侵蚀、极强度侵蚀、剧烈侵蚀、冻融侵蚀、坝子、水体、城市)分别建立训练样本。使用最小距离分类方法进行分类[5]。其中，训练样本的分离性可采用ENVI中ROI Separability Calculation工具进行计算。分离系数范围为0～2，大于1.9表明样本之间分离性好；小于1.8需要重新选择样本；低于1需合并分类样本。本次云南省侵蚀模数图样本分离系数中，冻融侵蚀与城市类型分离系数为0，两者区分度较低。可结合冻融侵蚀分布特点，后期在ArcGis软件中人工区分。

云南省侵蚀模数栅格图经过监督分类后，分类结果会不可避免的出现细小图斑，可利用ENVI软件中Majority/Minority分析、聚类处理(Clump)及过滤处理(Sieve)等功能进行处理[6-7]。经上述处理后，利用ENVI软件将分类结果文件转为ArcGis软件支持的Shapefile(.shp)矢量数据格式，即可得到云南省侵蚀模数矢量(图1)。

图1 云南省侵蚀模数矢量图制作流程示意

2 数据精度评估

为评价上述方法制作的云南省侵蚀模数矢量图的分类精度，本文借用Kappa系数进行评价。Kappa系数是指一种对遥感图像分类精度和误差矩阵进行评价的多元离散方法[8]。Kappa系数计算基于混淆矩阵，为创建混淆矩阵，综合考虑评估工作量大小，本文利用ArcGis [Create Random Points]工具随机生成200 个点，并利用云南省侵蚀模数栅格图人工赋值真实侵蚀类型，使用[Spatial Join]工具对随机点添加本次制作的云南省侵蚀模数矢量图侵蚀类型。云南省侵蚀模数矢量图精度评估点位置分布示意见图2。依据随机点属性表，通过简单的统计即可构建混淆矩阵(表1)[9]。

图2 云南省侵蚀模数矢量图精度评估点位置分布示意

一般Kappa系数计算值在0.00～1.00之间。其中，0.00～0.20表示分类质量差；0.21～0.40表示分类质量一般；0.41～0.60表示分类质量好；0.61～0.80表示分类质量很好；0.81～1.00表示分类质量极好[10]。计算得Kappa系数为1，表明本次云南省侵蚀模数矢量图分类成果质量极好，与原侵蚀模数栅格图一致，能够用于设计工作中。

表1 云南省侵蚀模数矢量图精度评估混淆矩阵

3 量算方法

基于ArcGis的水文集水区划分能够提供流域泥沙计算边界，计算边界因带有空间参考信息，可与步骤2制作的云南省侵蚀模数矢量图进行叠图分析，进而推求流域多年平均产沙量。结合水文集水区划分成果，可使用ArcGis中的 [Clip]工具对步骤2中云南省侵蚀模数矢量图进行裁剪，可利用[Summary Statistics]工具对裁剪结果文件按侵蚀类型分类汇总，统计各侵蚀类型面积。裁剪前设计流域与云南省侵蚀模数矢量图应具有相同坐标系，且裁剪后需要更新裁剪成果属性表的面积字段。

依据SL 190—2007《土壤侵蚀分类标准》，土壤水力侵蚀强度共分为6级，分别为微度、轻度、中度、强烈、极强烈及剧烈。在设计计算中，云南省侵蚀模数图中城市、水体、坝子无明显侵蚀类型归为微度级别，轻度侵蚀类型归为轻度级别，中度侵蚀类型归为中度级别，强度侵蚀类型归为强度级别，极强烈侵蚀类型归为极强级别，剧烈侵蚀类型归为剧烈级别。在缺少实测及调查侵蚀模数资料时，不同级别侵蚀模数可参考SL 190—2007《土壤侵蚀分类标准》进行取值。设计中使用的云南省侵蚀模数图不同侵蚀类型对应侵蚀模数取值见表2。

表2 云南省侵蚀模数图不同侵蚀类型对应侵蚀模数取值

依据各侵蚀级别面积，结合表1中不同侵蚀级别对应的平均侵蚀模数，可依据式(1)计算设计流域平均侵蚀模数。依据计算得到的流域平均侵蚀模数采用式(2)计算得到设计流域多年平均产沙量。

Mave=∑Mi·Ai

(1)

SA=Mave·F

(2)

式中Mave——流域平均侵蚀模数，t/(km2·a)；Mi——不同侵蚀类型侵蚀模数，t/(km2·a)；Ai——不同侵蚀类型面积，km2；F——设计流域面积，km2；SA——设计流域多年平均产沙量，t/a。

4 实例应用

4.1项目区概况

金龙河位于云南省中西部，属澜沧江流域，是澜沧江一级支流黑惠江的上游段，在汇入剑湖之前称为金龙河，出剑湖口之后始称黑惠江。金龙河发源于丽江市古城区九河乡白汉场罗凤山，先后流经九河、三河、梅园、金龙、金星等地，于金华镇柳营村附近汇入剑湖。金龙河流域面积674.1 km2，河长40.6 km，河道比降3.7‰。下文将以此流域多年平均产沙量计算为例说明本方法的实现过程。

4.2量算过程

本文以金龙河流域范围ASTER GDEM 2 DEM为水文分析地形数据[11]，利用ArcGis水文分析工具进行水文集水区划分，初步确定金龙河流域边界。在地势平坦地区，利用ASTER GDEM 2数据进行水文分析，其集水区划分成果可能存在较大偏差。本文采用该流域附近1∶5万扫描地形图为底图，对水文集水区划分成果进行校正。

使用校正后的金龙河流域边界，依本文前述量算方法对云南省侵蚀模数矢量图进行裁剪，分类统计各侵蚀类型面积。金龙河流域侵蚀模数图类型有坝子、城市、水体、无明显侵蚀、轻度侵蚀、中度侵蚀、强度侵蚀。对比金龙河流域侵蚀模数矢量图斑(图3)与原侵蚀模数栅格图可知，金龙河流域侵蚀模数矢量图图斑边界与栅格图图斑边界几乎完全一致，表明基于ENVI和ArcGis的云南省侵蚀模数图量算方法是可行的。本次金龙河流域共涉及140个图斑，采用本法能够大大缩短侵蚀模数图量算时间，提高设计人员的工作效率。

a)金龙河流域侵蚀模数矢量

b)金龙河流域侵蚀模数栅格矢量边界对比图3 金龙河流域侵蚀模数矢量、栅格对比情况

按照表1侵蚀级别合并侵蚀类型，经统计，微度级别面积占比72.6%，轻度级别面积占比19.3%，中度级面积别占比8.0%，强度级别面积占比为0.1%。依据侵蚀级别面积汇总成果，按式(1)计算流域平均侵蚀模数。在实际工作中，通常拟定3个计算方案(平均侵蚀模数取高值、中值及低值)，并结合相关资料及工程经验选取合适的计算方案。本次考虑到金龙河流域主要以微度侵蚀为主，且随着流域内小流域综合治理工作的深入开展，流域水土流失将减少。因此，本次金龙河不同侵蚀级别侵蚀模数取低值方案，即微度级别取200 t/(km2·a)，轻度级别取500 t/(km2·a)，中度级别取2 500 t/(km2·a)，强度级别取5 000 t/(km2·a)。依据计算，金龙河流域多年平均侵蚀模数为448 t/(km2·a)，流域多年平均产沙量为30.2 万t。

5 结语

本方法对于单幅侵蚀模数图量算同样适用，只是数据制作步骤1更加简单。从方便工作角度出发，建议制作全云南省侵蚀模数矢量图。本文采用的侵蚀模数图量算方法具有量算速度快、精度高、计算简单的优点。可有效减少设计人员工作量，对设计人员实际工作有较大帮助。同时，在使用设计流域裁剪侵蚀模数矢量图后，可将裁剪结果与侵蚀模数栅格图进行对比，根据项目需要对局部不一致的地方进行修正。