北京师范大学地理学与遥感科学学院.ppt

1、北京师范大学地理学与遥感科学学院,第一次全国水利普查技术细则,提 纲,一、水土流失普查遥感影像处理技术细则 二、水土流失普查专题信息提取技术细则 三、水土流失普查强度判别与成果处理技术细则,一、水土流失普查遥感影像处理技术细则,1.基础数据2.遥感数据预处理3.质量控制4.成果要求,1. 基础数据,1.1 HJ-1-A 、B数据HJ-1-A和HJ-1-B卫星数据对地刈宽为700公里、地面像元分辨率为30米、4个谱段（蓝、绿、红、近红外），重访周期为4天。覆盖全国需约80景HJ-1-A、HJ-1-B星CCD数据。,HJ-1-A卫星,洞庭湖-HJ星数据,HJ-1-A和HJ-1-B卫星CCD数据用于

2、计算植被覆盖，修正土地利用。数据获取基于中国资源卫星应用中心的数据查询系统（快视）进行，具体要求： （1）2009年至2010年季度数据； （2）不同侵蚀类型调查区图像的时相一致或相近； （3）图像清晰，地物层次分明，色调均一； （4）图像没有坏行、缺带，没有条带、斑点噪声和耀斑； （5）云层覆盖少的图像（即以晴空图像为优）； （6）数据覆盖全国。,1.2 MODIS数据,MODIS数据波段范围广，有36个波段，数据空间分辨率包括250米、500米和1000米三个尺度。,普查使用MODIS NDVI数据和MODIS Landcover数据，空间分辨率1km，用于计算全国植被盖度年内变化曲线，具

3、体要求： （1）2000年至2010年的数据，时间周期16天； （2）为实现HJ-1-A、HJ-1-B星CCD数据产品与MODIS数据 产品匹配，选取时间要尽可能一致； （3）对于不同侵蚀类型调查区图像的时相一致或相近； （4）图像没有坏行、缺带； （5）尽量挑选云层覆盖少的图像（即以晴空图像为优）； （6）数据覆盖全国,1.3 Aqua/AMSR-E数据,Aqua/AMSR-E亮温数据是辐射计数据，来自美国冰雪数据中心，数据在网上免费下载（https:/wist.echo.nasa.gov/api/），像元大小为0.25，用于冻融侵蚀温度相关因子反演。 具体要求： （1）2002年至2010

4、年数据，时间周期2天； （2）图像没有坏行、缺带； （3）数据覆盖冻融区。,2.1 几何精校正,2. 遥感数据预处理,2.1.1 HJ-1-A、HJ-1-B数据,几何精校正流程图,(1) 数据挑选数据源选取需能有效的覆盖普查范围，保证数据质量和方便色调处理。具体标准有：1）相邻区域时相一致或相近的图像；2）图幅方正，图像清晰，地物层次分明，色调均一；3）图像没有坏行缺带，没有条带、斑点噪声和耀斑；4）尽量挑选云层覆盖少的图像（即以晴空图像为优）；5）为了保证图像色调均一，尽量选用同一季节的图像。 (2) 几何精纠正预处理系统提供的标准数据产品是经辐射校正和系统几何校正的2级产品，只进行了单景数

5、据波段间的配准、纵横向随机条纹的基本滤除和CCD影像色调的平衡归一化校正。需对获取的2级标准产品影像进行以下处理：1）选取其中的1、2、3和4波段，采用遥感图像处理软件如ERDAS IMAGINE、ENVI等，处理采用ERDAS IMAGINE软件实现波段合成，转换成IMG格式文件。2）对有噪声的图像进行去噪声处理以及对模糊图像进行拉伸增强处理，使图像清晰化，调整影像的亮度、对比度和色阶比，达到图像最佳效果。,(3) 图像几何精校正图像的几何精纠正以1:50000和1:100000地形图或（同等或更高）分辨率卫星影像作为参考，人工选取控制点进行几何精纠正。纠正结果文件采用正轴等面积割圆锥投影，

6、又称亚尔勃斯（Albers）投影，以25N和47N两条纬线为标准纬线相割；椭球体采用国际通用的WGS-84；数据格式为Geotiff。几何精校正用ERDAS IMAGINE软件完成。纠正算法采用多项式校正，每幅影像的控制点均匀分布，选取20个以上控制点；平原地区可采用2次多项式进行几何校正，山区需采用3次多项式进行几何校正。几何精校正后的影像空间采样分辨率：HJ-1-A和HJ-1-B卫星CCD数据产品为30米。几何精校正产品格式：Geotiff；几何精校正的命名规则：Path-Row-卫星标识-获取日期-ref。,(4) 数据地形图分幅 数据分幅工作包括镶嵌及编辑处理、地形图分幅处理、质量控制

7、三步。1）镶嵌及编辑处理经几何精纠正后，以较高的精度把各景影像拼接起来，形成覆盖全国的整幅影像。2）地形图分幅处理采用基于网格新裁切方案，具体按1:250000比例尺要求，生成覆盖某个区域（如省为单位）的网格矢量线，利用网格线对该区域图像进行自动裁切，生成一系列该比例尺的影像地图。,Aqua/AMSR-E数据校正，利用头文件信息直接将地理坐标信息赋给图像，在ENVI软件完成。,2.1.2 MODIS数据,MODIS数据校正，利用头文件信息直接将地理坐标信息赋给图像，采用ENVI软件完成。,2.1.3 Aqua/AMSR-E亮温数据,2.2 大气效应纠正,利用HJ-1-A和HJ-1-B星的蓝光(

8、Alt+1)、绿光（Alt+2）、红光（Alt+3）和近红外（Alt+4）共4个通道的遥感影像，以及与其相对应的头文件（*.xml），根据灰度和RGB不同波段组合方式浏览HJ-1-A和HJ-1-B星影像，并判断是否存在云（或云阴影）、水体和浓密植被。,2.2.1 图像信息获取,2.2.2 图像预处理,利用公式（1）将DN值转换为表观辐亮度,（式 1）,是原始图像像元灰度值； 是绝对辐射定标系数，可从头文件或中国资源卫星应用中心网站中获取。,公式（2）将表观辐亮度转换 为表观反射率,（式 2）,是表观反射率；,是日地距离纠正因子；,是大气外太阳光谱辐照度；,是太阳天顶角。,上述参数中有些可从原始

9、数据头文件获得， 或http:/网站中获取。,图像的暗目标自动提取依据比值植被指数 (RVI)，土壤调整植被指数(SAVI)和归一化水体指数 (NDWI)的综合分析法实现，各指数的计算公式如下：,（1）最小反射率法最小反射率法利用浓密植被和水体分别在遥感影像的蓝光、近红外通道具有非常小的反射率来自动提取暗目标并获取大气参数。技术流程:依据比值植被指数(RVI)，土壤调整植被指数(SAVI)和归一化水体指数(NDWI)的综合分析法，通过决策树方法逐步实现对清洁水体和蓝光波段浓密植被作为图像暗目标的自动提取。,2.2.3 大气效应纠正,利用自主开发的大气订正软件完成。HJ-1-A和HJ-1-B卫星

10、数据的大气纠正主要采用了最小反射率法、MODIS辅助法与气象台站数据法三种方法。,图像暗目标自动提取技术流程图,最小反射率法大气纠正技术流程图,操作步骤: 1)将HJ-1-A、HJ-1-B的CCD数据产品转换为大气顶的表观辐亮度和表观反射率： 2)软件建立了以气溶胶光学厚度（AOD）和太阳天顶角（ ）为索引的查找表（Look-Up-Table），气溶胶类型主要是大陆乡村型。 3)在实现图像暗目标自动提取后，依据查找表获取气溶胶光学厚度； 4)依据气溶胶和太阳天顶角，通过查找表获取其它大气参数，进行大气纠正。,（2）MODIS辅助法技术流程:利用已经进行大气纠正的现成遥感图像数据产品对卫星影像进

11、行大气纠正得到地表反射率。,1）采用MODIS星上定标产品得到MODIS图像获取时太阳在大气顶的辐照度，将定标好的MODIS L1B辐亮度产品MOD02转换为大气顶的表观反射率； 2) 将得到的MODIS大气顶反射率和地表反射率产品MOD09配合解得大气中气溶胶光学厚度和水汽含量等参数； 3) 假设MODIS与HJ过境时大气状况没有发生变化，再利用简单的大气辐射传输模型得到HJ-1星角度上的大气反射率、大气透过率，进行大气纠正,最终得到HJ-1星的地表反射率。,大气效应纠正技术流程图,操作步骤: 1) 定标好的MODIS L1B辐亮度产品MOD02转换为大气顶的表观反射率。 2）MODIS大气

12、顶反射率和地表反射率产品配合解得大气参数。 3）HJ-1-A、HJ-1-B卫星数据大气效应纠正。,（3）气象台站数据法当获取与HJ-1-A、HJ-1-B卫星数据相应的气象资料数据时，可以用于提出大气纠正所需的大气参数。气象数据应包括：蓝光(440nm)、绿光(550nm)、红光(660nm)和近红外(880nm)四个波段任意两波段的气溶胶数据和大气柱水汽含量(单位：g/cm2)，输入的气象数据中，AOD和CWV不能小于0，或可直接输入气象数据进行大气纠正。,2.3 角度效应纠正,将HJ-1-A、HJ-1-B星数据方向性地表反射率产品（30m）进行角度效应纠正，得到垂直向下观测的归一化植被指数N

13、DVI。 技术流程 1) HJ-1-A、HJ-1-B星数据方向性地表反射率根据NDVI的定义得到带有方向性特征的NDVI数据。 2) HJ-1-A、HJ-1-B星数据方向性NDVI数据在HJ-1-A、HJ-1-B星30m像元为均匀植被的假设前提下使用简单的余弦纠正初步得到垂直观测的NDVI数据。,影像角度效应纠正技术流程图,操作步骤:HJ-1-A、HJ-1-B星数据方向性地表反射率根据NDVI的定义得到带有方向性特征的NDVI数据。,3.1 质量控制内容,质量控制是卫星数据产品质量保证的重要途径，具体内容包括HJ-1-A、HJ-1-B卫星数据产品质量检查、几何精校正产品质量精度检查、大气纠正产

14、品质量和精度检查。利用典型地物光谱库数据或实测地物反射率数据HJ-1-A、HJ-1-B星大气纠正产品进行检查，分析并估算大气纠正产品的精度。,3. 质量控制,3.2 质量控制方法,（1）过程抽检成立项目质量检查组，在项目执行过程中，定期开展质量抽查，发现问题及时解决，严格控制把关各个作业环节的质量，保证优质高效地完成项目。 （2）阶段成果的“二检”在项目实施过程中，针对阶段成果，严格执行作业员自检、作业员互检的“二检”制度。作业员自检：作业人员自己进行全面的检查工作，检查比例可根据本身的作业水平决定。作业员互检：作业员之间的交换互检，检查比例为100%，并作详细记录。 （3）最终成果预检组织专

15、家会同项目技术组组成专门检查小组，在上述抽检和二检的基础上，根据质量检查组检查报告，对项目的最终成果进行全面的预检，填写相应的记录表。,（1）辐射精度中等反射率地物，HJ-1-A、HJ-1-B星大气纠正产品的误差在10%左右； （2）空间精度HJ-1-A、HJ-1-B星数据几何精纠正精度：平原地区小于1.5个像元,山区小于3个像元。分幅接边精度：平原地区2.5个像元，山区3.5个像素。,3.3 精度要求,4.1 数据投影,成果数据投影采用正轴等面积割圆锥投影，即亚尔勃斯（Albers）投影，空间数据参数如下：,4. 成果要求,HJ-1-A、HJ-1-B卫星、MODIS、Aqua/AMSR-E校

16、正数据产品格式为Geotiff格式。依据国家基础地理信息1：25万数据标准分幅存储。影像地图标注县、乡镇地名及省级公路，按县域分幅成图。,4.2 数据格式与分幅,二、水土流失普查专题信息提取技术细则,（一）水力侵蚀专题信息提取技术细则 （二）风力侵蚀专题信息提取技术细则 （三）冻融侵蚀专题信息提取技术细则,（一）水力侵蚀专题信息提取技术细则,（1）降雨侵蚀力因子：是指降雨导致土壤侵蚀发生的潜在能力，用一次降雨总动能E与该次降雨最大30min雨强I30的乘积EI30表示。反映了雨滴对土壤颗粒的击溅分离以及降雨形成径流对土壤冲刷的综合作用。 （2）降雨侵蚀力季节分布：是指一年中某时段降雨侵蚀力占全

17、年降雨侵蚀力的百分比，用作权重因子计算水土保持生物措施因子值。 （3）降雨侵蚀力等值线图：空间上多年平均年降雨侵蚀力相等点的连线称为等侵蚀力线，由等侵蚀力线构成的空间等值线分布图称为降雨侵蚀力等值线图，反映了多年平均降雨侵蚀力的空间变化特征。,1. 降雨侵蚀力因子,1.1 术语,相关数据包括：各县建立的所属气象站1981-2010年逐日雨量电子数据（.dat格式或登记表电子文档），及其说明文件；装订成册的各县“气象数据登记表”。 按降雨侵蚀力因子计算公式要求，进行以下数据处理： （1）非侵蚀性降雨剔除。如果日雨量小于（不含）12mm，则将该日雨量设为0。 （2）处理后的数据按原格式（.dat）

18、重新存入相同目录下，命名为“侵蚀性降雨”。,1.2 数据处理,1.3 因子计算,质量控制包括两个方面：一是对降雨资料进行质量控制。每省最少抽查1个县2-3年的记录，核对气象数据登记表数据及其电子数据。二是对降雨侵蚀力因子计算结果进行抽查。在水蚀区最少抽查一个县的计算结果，与利用已有分钟降水资料计算的降雨侵蚀力因子进行对比，精度75%。,1.4 质量控制,1.5 成果构成,（1）全国各县多年平均年降雨侵蚀力因子值。 （2）全国降雨侵蚀力等值线图，及其30m网格的栅格数据。 （3）全国各县多年平均24个半月降雨侵蚀力占年降雨侵蚀力百分比。 （4）全国多年24个半月降雨侵蚀力占年降雨侵蚀力百分比等值

19、线图，及 其30m网格的栅格数据。,数字高程模型：表现某高程基准下地面高程空间分布的有序数字阵列。数字高程模型可以基于地形图、地面或遥感测量等方式获取的高程数据，经内插建立。 坡度因子：指CSLE的坡度因子，定义为某一坡度土壤流失量与坡度为5.13，其它条件都一致的坡面产生的土壤流失量之比率。对于面上的普查，采用每个栅格的坡度因子值参与土壤流失量的计算，对于抽样调查，将取每个评价单元坡度的中位数参与土壤流失量的计算。 坡长因子：指CSLE的坡长因子，定义为某一坡面土壤流失量与坡长为22.13m、其它条件都一致的坡面产生的土壤流失量之比率。用法同坡度因子值。,2. 坡度坡长因子,2.1 术语,数

20、据处理过程包括地形图数据处理、流水线和分水线提取以及DEM生成。 （1）地形图数据处理：需要1:50000地形图、1:250000地形图和100m分辨率DEM。对国家测绘局基础地理信息中心提供的数字化地形图或自行数字化的地形图，经过必要修改编辑和质量控制，以确保每条等高线和每个高程点上标注有合理的高程值、河流由高向低流动、河流和等高线位置关系正确。 （2）流水线和分水线提取：利用GIS水文地貌分析功能，通过编程自动提取流水线和分水线。以100m分辨率DEM为基础，划分出各大流域的2级支流。划分出的流域单元，面积控制在1*1042*104km2。平原区流域划分和河流提取结果，须与比较大比例尺地形

21、图对照，并做出必要的修编。,2.2 数据处理,（3）DEM生成：坡度坡长因子的提取，以1:10000和1:50000数字地形图建立的DEM为数据基础。其中基于1:10000地形图的DEM适用于野外调查单元对坡度坡长因子的计算，基于1:50000地形图的DEM适用于全国坡度坡长因子的提取与计算。利用经过质量控制的数字地形图（包括等高线、高程点、河流、湖泊和水库等要素）和专业软件ANUDEM，经过插值生成DEM。1:10000和1:50000地形图插值，在丘陵地区和山区，分辨率设置为10m和25m，平原地区和东北漫岗丘陵区分别设置为5m和10m。,2.3 因子计算,以DEM为基础提取坡度坡长因子，

22、根据坡度坡长因子的定义和坡度坡长因子值的算法（式1-11、式1-12和式1-13），利用LS计算专用程序（基于C+语言开发）完成流域坡度坡长因子值提取与计算 。,LS计算专用程序界面,地形因子的质量控制采用分层次、分区控制法。首先对DEM建立、坡度和坡长的计算分别控制；其次在主要水土流失类型区各选3-5个样点，通过实测坡度和坡长，完成对质量的评价，要求精度达到75%。,2.4 质量控制,2.5 成果构成,（1）全国DEM、坡度、坡长、坡度坡长因子栅格文件，精度为30m网格。 （2）以1:250000地形图图幅为基础，对流域坡度坡长因子值进行重新组织（拼接和切割），完成图幅坡度坡长因子值。该因子

23、值可以1:250000地形图编号来命名。,土壤可蚀性因子：表征土壤被冲被蚀的难易程度，反映土壤对侵蚀外营力剥蚀和搬运的敏感性，是影响土壤侵蚀的内在因素。国际上常用K表示。K值的大小是由土壤性质本身所决定的，在众多的土壤性质中土壤颗粒组成和有机质对K值影响最大，这也是利用土壤性质在宏观尺度上估算土壤可蚀性的理论基础。理论上的获取方法，是指标准小区上单位降雨侵蚀力引起的土壤流失量，单位为thm2 h/（hm2 MJmm），K值是进行土壤侵蚀和水土流失定量评价的重要依据。,3. 土壤可蚀性因子,3.1 术语,完成全国水蚀区土壤可蚀性K值的计算及其分布规律，一方面需要土壤数据库，包括土壤亚类和土属空间

24、数据和属性数据；此外，为满足本次土壤侵蚀普查的精度要求，K值的计算需要精确到全国909个土属，而且计算K值的土壤属性数据必须对应到中比例尺的水蚀区土壤类型图上。具体的数据获取和处理方法说明如下： （1）土壤基本属性和空间数据获取 收集全国34个省市的土种志和1:500000土壤类型图，对照中国土壤发生分类系统，查阅各土种的理化性质和机械组成。主要包括土种的名称、土种所在的地点、分布面积、该土种所属的亚类和土属、表层的有机质含量(%)、粗砂2-0.2mm(%)、细砂0.2-0.02mm(%)、粉砂0.02-0.002mm(%)和粘粒0.002mm(%)。将这些数据进行预处理，建立土壤理化性质的属

25、性数据库。同时，对于纸质版的土壤类型图进行几何纠正、配准和数字化，建立全国1：50万土壤类型图图形数据库，用于土壤可蚀性的计算。,3.2 数据处理,（2）试点县的野外采样过程 结合土壤侵蚀野外抽样调查工作，完成单元内土壤样品的采集工作。在试点省份的试点县内，将野外调查布点图与土壤分布图进行叠加分析，取该点土壤样品，并集中到省级普查中心，进行土壤理化性质。土壤样品采集方法如下： 1）记录样点信息。土壤样品采样点用GPS进行空间定位，记录经纬度坐标，在信息表上记录采样点基本信息，包括地点（县、市、乡、村）、土壤名称、位置坐标、土地利用类型、地貌部位（坡度、坡向）、水土保持措施、植被盖度等。 2）取

26、样。选择面积约20cm20cm的小区域，清除土层表面枯枝落叶，挖一个小垂直断面，用铁锹垂直取1个土柱，土柱深约20cm，重量大约500-800g土样即可，装入密封袋。在山区，土壤样点主要依据地形选择：3个样点一般分别分布在坡面上部、中部和底部，同时还要兼顾阴阳坡和土地利用类型。在平原区，可在调查单元内随机选择3个土壤样点，要避免3个样点之间的距离太近，最后混合装袋成一个样品。 3）填写标签与封袋。用铅笔填写标签（包括采集日期、地点、土壤名称、采集人等）。一式两份。一个放入袋中，一个贴于土样袋外侧。将密封袋中的空气慢慢挤出，封好密封袋，采样完毕。 4）采样点处理。用铁锹将取样时从样点挖掘出的土壤

27、回填到样点，可用脚踩实，以枯枝落叶覆盖样点。 5）将野外采集的土壤样品带回实验室，进行理化分析，分析试验指标包括土壤有机质含量和土壤机械组成，机械组成采用美国制标准，分5个粒级，即：2-0.1mm、0.1-0.05mm、0.05-0.02mm、0.02-0.002mm和0.002mm。,（3）共享分析测试中心的土壤理化分析数据 通过共享中国科学院南京土壤研究所分析测试中心自1980年以来所有的土壤理化分析资料，计算近年来不同土壤类型的K值，用于更新部分土壤类型的K值。主要共享全国范围内的土壤样品的采集地点、土壤类型、地理坐标、机械组成、有机质含量等数据，用于计算土壤可蚀性K值。（4）整理已经发

28、表的文献资料 利用CSCD文库，查阅自1980年以来所有文献上的相关土壤类型的土壤理化分析结果，以及相关的土壤可蚀性资料，用于补充和修正利用第二次土壤普查资料获取的K值。,3.3 因子计算,以全国数字化土壤类型图为蓝本，以全国第二次土壤普查的样点数据（土种）为基础，先计算每个土种的土壤可蚀性K值。再利用中国土壤发生分类系统中土属与土种的对应关系，将土种的K值利用面积加权平均法归并到土属上来。再将土属的K值分级后链接到全国土壤类型图上来，得到全国土壤可蚀性K值分布图，具体的流程图如下：,关于土壤可蚀性K值的计算，在获取了上述土壤理化分析资料后，可以根据理化性质的具体情况，分别选用Williams

29、模型和Wischmeier模型。,（1）Williams模型Williams等人在EPIC（ErosionProductivity Impact Calculator）模型中发展了土壤可蚀性因子K 值的估算方法，只需要土壤有机碳和颗粒组成资料，,（2） Wischmeier模型,如果利用第二次土壤普查资料或者获取的土壤分析资料不是美国制，在计算土壤可蚀性K值前就必须进行土壤颗粒分析结果的转化，即由国际制转化为美国制。 因为第二次土壤普查中土壤颗粒分析标准是国际制。在国际制向美国制的转换中（见表1-4），一般采用具体用图解法。即在半对数纸上先画出国际制的土壤颗粒大小分配曲线，然后查出美国制某一粒

30、径的百分数含量。,此方法是在计算机还没有普及的情况下进行使用，根据这个思路，我们进行了拟合曲线法，得到了很好的拟合效果。即应用Y=aX+b和Y=aX2+bX+c进行拟合, 得出了每个土壤类型的转换方程，其中两个方程中X=ln(P), P为粒径大小(mm)，Y是小于P粒径的累计颗粒含量百分数(%)。,数据获取和处理过程中肯定存在一定的误差，都可能影响到最终成果的准确性。为了尽可能提高K值的计算精度，需要在计算流程的多个关键环节进行质量控制： （1）为保证土壤属性数据的代表性和全面性，在将土种的理化性质向土属归并时，需保证该土属内的土种数占该土属内所有土种总数的85%以上。 （2）在进行土壤样品采

31、集时，采样点的布设需保证新采集的土壤样品数能够覆盖水蚀区试点省份80%以上的亚类。 （3）在进行土壤样品分析时，标样和样品之间的误差不能大于5%。 （4）在进行1:500000土壤类型图数字化和属性赋值过程中，要保证土壤图斑空间位置和属性信息的准确性大于95%。 （5）在K值计算中存在国际制向美制单位转换的问题，需要借助于转换方程完成，在计算过程中要保证转换方程的确定性系数（R2）不小于90%。,3.4 质量控制,（1）建立全国1:500000土壤分布图与第二次土壤普查资料的土壤属性数据之间的关系，构建全国土壤基础属性和空间信息数据库。（2）通过采样和更新土壤理化性质，获取全国909个土属的最

32、新可蚀性K值，构建全国土壤可蚀性因子数据库，再将K值链接到土壤类型图上，进行分级后可得到全国土壤可蚀性因子图，精度为30m网格。,3.5 成果构成,（1）生物措施因子：一定条件下耕作农地上的土壤流失量与同等条件下连续休闲对照裸地上的土壤流失量之比，为一无量纲数，其值大小变化于01之间。 （2）植被盖度：是指植被（包括叶、茎、枝）在单位面积内植被的垂直投影面积所占百分比。 （3）归一化植被指数（NDVI）： NDVI=(NIR-R)/ (NIR+R)，常用遥感观测的近红外波段反射率（NIR）和红光波段反射率（R）计算，可用于检测植被生长状态、植被盖度和消除部分辐射误差等。,4. 生物措施因子,4

33、.1 术语,在获取时间序列高时空分辨率NDVI数据的过程中，用到了以下数据： （1）30m左右的高空间分辨率的数据：HJ-1多光谱； （2）全国2008年1:100000土地利用图； （3）低空间分辨率时间序列的NDVI数据：MODIS 反射率产品（MCD43B4，Nadir BRDF-Adjusted Reflectance, NBAR），空间分辨率为1km，时间分辨率为16天，时间序列5年，用于生成可靠的NDVI 时间序列； （4）低空间分辨率的MODIS分类产品：MOD12Q1，空间分辨率为1km，时间2004年。,4.2 植被盖度数据处理,4.2.1 数据说明,时间序列高分辨率NDVI

34、计算的具体流程如下图,4.2.2 时间序列NDVI计算,（1）HJ-1数据和MODIS 数据的配准MODIS产品的投影方式与HJ-1数据的UTM投影不同，将它们投影转为一致。,（2）HJ-1 NDVI分布图在对HJ-1多光谱数据大气纠正和角度订正的基础上，利用NDVI的定义（如公式所示），生成不同时相HJ-1 NDVI分布图。,（3）各个不同地类MODIS NDVI时间序列数据的提取MODIS分类产品（MOD12Q1）中包含五大分类体系，根据应用目的选取了的植被功能分类（PFT）产品，其包含地物类别分别为：常绿针叶树、常绿阔叶树、落叶针叶树、落叶阔叶树、灌木、草地、谷类作物、阔叶作物、城镇、水

35、体等。,选取多个混合像元后，通过线性模型分解，利用最小二乘的方法便可以得到亚像元上各地类的NDVI。 （4）高空间、高时间分辨率NDVI产品的生成融合MODIS NDVI 和HJ-1 NDVI数据得到HJ-1空间尺度上全年半月尺度的各地物类NDVI数据产品。,4.2.3 植被盖度计算,从NDVI提取植被盖度产品的技术流程如下图,（1）NDVI到植被盖度转换系数确定以MODIS NDVI产品为辅助判断标准，同时判断每一种气候地理区划里每一地类一年中MODIS NDVI的最大和最小值。在每一种气候地理区划里找到一年中每一30m分辨率不同植被类型的直方图中最大的2-5%左右选取为NDVImax， N

36、DVImin的选取固定不同植被类型周边小范围空间区域内裸土地类的NDVI平均值。,（2）NDVI到植被覆盖度转换,结合地面调查单元的植被覆盖度实测数据，通过以上公式计算得到时间分辨率15天空间分辨率30m的植被盖度。（3）植被盖度转换系数的确定和检验通过空间定位得到与地面调查单元测量点对应的遥感NDVI图像像元，选用部分地面实测盖度和遥感像元NDVI数据拟合NDVImax， NDVImin与k。使用另外部分数据对转换系数进行精度验证。,进行植被盖度季节分布监测，获取不同地区主要植被类型的植被盖度季节分布标准曲线，这是计算生物措施B因子的一个重要基础，可以对宏观遥感调查计算生物措施B因子中，对遥

37、感反演的植被盖度进行验证。（1）植被盖度动态监测准备依托各地水土保持监测站的试验小流域开展工作，依据全国水土流失类型区，选择30个小流域为样点进行植被盖度动态监测（表1-5）。,4.2.4植被盖度季节分布标准曲线,各小流域工作人员调查收集本流域的土地利用状况，确定本地区的土地利用类型，并填写右表（表1-6）。,（2）野外选择调查地块要求在正式进行野外植被盖度照相测量前，确定野外地块的类型、数量和位置。每个小流域（样点）按农地、人工草地、天然草地、人工林地、果园等5种类型共选择12个地块进行观测。具体方法如下：1、地块的选择应该代表该流域某种土地利用类型的主要生长状况，同时兼顾组织实施的方便性，

38、如靠近交通线；2、每个地块的植被覆盖均匀，面积尽量在6060m2以上，如果达不到可适当降低标准，但地块面积应尽可能相对大一些；3、各土地利用类型地块的选择和数量分配：农地结合该地区主要耕作制度进行，选择2种代表性作物,每种作物选择2个地块；人工草地、天然草地、人工林地、果园等每种类型各选择2个地块，果园应该选择当地主要的果树类型；,如果设定的某种植被类型在本流域没有，则可以在其它植被类型上增添地块的数量，使得地块的总数量保持控制为12个。a)根据手持GPS，确定每个地块的位置；b)生成流域（样点）地块布设表和空间分布图。流域地块布设表记录了每个地块编号、植被类型、相对位置、经纬度坐标，给出农地

39、植被的播种时间、出苗时间和收割时间，和地块景观照片编号，并在表后附上每个地块的景观照片；空间分布图以主要交通线、流域界线（或行政界线）和标志地名作为底图,表示该样点各地块的空间分布。填表时，对于农地作物的播种、出苗和收获时间，如现阶段尚不能明确，可以在以后的植被盖度动态监测的过程中补充填写。c)确定全部地块的测量顺序。,（3）植被盖度季节分布监测实施植被盖度季节分布监测从2010年开始共实施两年，每半个月对所选择的全部地块进行一次植被盖度测量，要求两天之内完成全部地块测量，相当于一个测量周期。每年测量24次，每次测量要求至少2人，为保证全年测量的顺利进行，中途最好不更换测量人员。每次测量安排在

40、每个半月的后半段（即上半月的10号15号期间，下半月的2531号期间）进行，采用照相法测量植被盖度，同时测量植被高度，并填写野外照相法测量植被盖度调查表。最后及时上传照相测量记录表和数码植被图像，根据专用软件及时计算植被盖度。,（4）植被盖度计算采用植被盖度自动计算系统（PCOVER）自动计算植被盖度。植被盖度自动计算系统（PCOVER）（计算机软件著作权登记号：2008SR12421）是专为便携型植被盖度摄影仪所开发的配套软件，运行环境为Windows 2000或Windows XP，系统安装要求约30M硬盘空间。（5）获取植被盖度季节分布标准曲线完成一年植被盖度动态监测后，对计算的植被盖度

41、信息进行整理，完成各小流域主要植被类型的植被盖度季节分布标准曲线的绘制。,4.3 生物措施因子计算,以时间分辨率15天、空间分辨率30m的植被盖度影像系列为基础，结合时间分辨率15天、空间分辨率30m的的降雨侵蚀力比例图和分辨率30m的土地利用图，进行生物措施B因子的计算。具体步骤：,（1）结合分辨率30m的土地利用图，按农、林、草等3大类型，根据时间分辨率15天、空间分辨率30m的植被盖度影像系列分别计算每半月时段的土壤流失比率SLR（水体和建筑用地忽略），得到各半月时段的土壤流失比率SLR图。 （2）根据各半月时段的土壤流失比例SLR和各半月时段降雨侵蚀力占全年的比率，计算得到生物措施B因

42、子图：,采取两种方式对基于遥感的生物措施B因子进行质量控制。 （1）遥感反演的植被盖度季节分布进行验证和质量控制 全国选择30个样点，每个样点按农地、人工草地、天然草地、人工林地、果园等选择2个地块，确定每个地块的经纬度位置。 监测获取每个地块的植被盖度季节分布数据。 从时间分辨率15天、空间分辨率30m的植被盖度影像系列上，获取30个验证样点每个地块的植被盖度季节分布数据。 针对30个样点的每个地块，对实测的地表植被盖度季节分布和遥感反演的植被盖度季节分布进行对比分析，精度75%。 （2）生物措施B因子进行验证和质量控制 在每个水蚀类型区，统计整理该类型区的径流小区的侵蚀资料和生物措施B因子

43、计算成果，精度75%。 根据历史径流小区计算的生物措施B因子结果与基于遥感的生物措施B因子计算进行对比分析。,4.4 质量控制,（1）全国30个小流域主要植被类型的植被盖度季节分布标准曲线 （2）全国半月尺度植被覆盖度图，精度为30m网格 （3）全国半月B因子图和年B因子图，精度为30m网格,4.5 成果构成,工程措施：是指通过改变小地形（如坡改梯等平整土地的措施），拦蓄地表径流，增加土壤降雨入渗，改善农业生产条件，充分利用光、温、水土资源，建立良性生态环境，减少或防止土壤侵蚀，合理开发、利用水土资源而采取的措施。 耕作措施：是以保水保土保肥为主要目的，以提高农业生产为宗旨，以犁、锄、耙等为耕

44、（整）地农具所采取的措施。 地块：是指空间上连续的、具有同一土地利用类型、同一覆盖度、同一水土保持措施的区域。,5. 工程措施和耕作措施因子,5.1 术语,（1）野外调查数据工程措施和耕作措施因子依据野外调查单元确定及数据采集技术规定中相关要求完成。 （2）高分辨率影像相结合的ET专题图处理工程措施因子E和耕作措施因子T因子主要依托土壤侵蚀野外调查单元获得的工程措施和耕作措施数据，并结合高分辨率的CBERS全色影像，采用人机交互式解译的方法获取典型区内的水土保持工程措施和耕作措施数据，最后以土地利用为依据，按照比例平均分配的方法获取。专题图以行政县为单位，根据该县1:100000的土地利用图，

45、在GIS软件的支持下，修改图层的属性表，添加工程措施因子E和生物措施因子T两个属性字段，将调查单元上获得的各土地利用类型的措施因子对应赋值到土地利用的属性表上，根据工程措施因子E和生物措施因子T两个属性分别生成工程措施因子E和生物措施因子T的专题图，从而实现措施因子的空间离散。,5.2 数据处理,ET因子获取技术流程图,5.3 因子计算,（1）措施数量统计在野外调查单元获取的工程措施和耕作措施的数据基础上，以各县的土地利用类型为基本计算单元，统计该县内所有野外调查单元内采取了各种工程措施和耕作措施的土地总面积，并计算其占野外调查单元内该土地利用类型的比例。 （2）ET因子计算按照工程措施、耕作

46、措施面积比例计算工程措施因子和耕作措施因子。某种土地利用类型的ET因子即为该土地利用类型中采取的每种措施的ET因子（基于野外调查单元）与该措施在该土地利用中的比例乘积的累积 ，,注：水土保持措施未涉及地区的工程措施（E）因子值尚在收集当中。部分地区工程措施（E）措施因子值将结合调查人员收集的径流小区资料进行计算、更新。,注：水土保持措施未涉及地区的耕作措施（T）因子值尚在收集当中。部分地区耕作措施（T）因子值将结合调查人员收集的径流小区资料进行更新。,工程措施和耕作措施因子的质量控制采用分层次、分区控制法。在主要水土流失类型区各选3-5个样点，通过实地调查检验，完成对质量的评价，要求精度达到7

47、5%。,5.4 质量控制,5.5 成果构成,（1）野外调查单元内获取的工程措施和耕作措施数据 （2）全国ET因子专题图，精度为30m网格,（二）风力侵蚀专题信息提取技术细则,风力因子：是指不同等级风速和不同风向对土壤风蚀发生的潜在能力，单位m/s。风速和风向采用气象台站观测数据。 临界侵蚀风速：是指风力作用于地表，产生土壤风蚀的最小风速。大于和等于临界侵蚀风速的各等级风速，均为侵蚀风速。考虑到流沙地表的临界侵蚀风速约5m/s（在各类地表状况中的临界侵蚀风速最小），因而在统计风力因子时，自5m/s的风速开始统计。 各等级风速累积时间分布栅格图：全年大于和等于临界侵蚀风速的各等级风速的累积时间（m

48、in）在空间上的分布，反映了各等级风速累积时间的空间变化特征，采用30m分辨率的栅格数据表达各等级风速累积时间的空间变化。,1. 风力因子,1.1 术语,在实际操作中，不同气象站和不同时段内，可能既有逐日整点风速风向数据，又有逐日4次观测的风速风向数据。为了统一数据格式，便于计算，提出以下2条技术途径。 （1）若收集到的风速和风向资料为逐日24小时整点风速风向，则按1m/s间隔直接统计全年5m/s的各等级风速的累积时间。统计步骤如下：1）按照给定的“逐日24小时整点风速风向统计表”完成逐日24小时整点风速风向电子版本（Excel格式）。2）统计全年5.05.9m/s风速发生的频数，每次代表1小

49、时，按下式换算为该等级风速的累积时间。累积时间=频数1hr/次60min/hr.为便于计算，介于5.05.9m/s的风速等级取平均值U5.5m/s，6.06.9m/s的风速等级取平均值U6.5m/s，依此类推。3）按步骤（2）方法，依次统计5.05.9m/s、6.06.9m/s、7.07.9m/s、各等级风速的累积时间，最高风速统计到该气象台站记录到的最大风速（Excel格式）。,1.2 数据处理,（2）有的气象台站如果仅有每天4次（观测时间为2时、8时、14时、20时）风速风向记录，则采用线性插值的方法，获取24小时整点风速估计值，然后按1m/s间隔统计全年5m/s的各等级风速的累积时间。统

50、计步骤如下：1） 按照给定的“逐日24小时整点风速风向统计表”完成逐日4次风速风向电子版本（Excel格式）。逐日4次的风速风向填写在对应的时刻（2时、8时、14时、20时），待完成步骤2和步骤3后，将各整点风速插值填入对应时刻。2） 按照线性插值方法，计算相邻两个整点时刻的风速，,3）完成逐日24小时整点风速风向统计表（Excel格式）。 4）按照“方法一”中“3）”的相同方法，统计全年各等级风速发生的频数。,1.3 因子计算,考虑到风力因子在空间上的变化是渐进的，即风速和风向在空间上是逐渐变化的，相邻气象台站之间的风速和风向变化呈逐渐过渡关系。因此，风力因子计算可以在GIS软件和自主开发程序的支持下实现自动运算。 （1） 逐年统计19812010年15月和1012月的逐日逐时不同风速等级和相应风向。 （2）将每个月分为上半月和下半月（每月115日、1630日或31日）两个时段，统计每半月时段内各风速等级的累积时间（单位：min）和16方位风向出现频率。例如，某气象台站某一等级的风速累积时间和16方位风向出现频率为：累积时间（min）频数1hr/次60min/hr.式中的各风力等级按照5.05.9m/s取其中值5.5m/s，6.06.9m/s取其中值6.5m/s，7.07.9m/s取其中值7.5m/s，。,