县域尺度水土流失抽样监测技术研究

马松增，何明月

(1.河南省水土保持监督监测总站，郑州 450008；2.河南省水利勘测有限公司，郑州 450003)

在水土流失动态监测中，美国主要采用抽样调查和模型法(USLE、WEQ)，澳大利亚及欧洲主要采用专家法和模型法(USLE、CEMS)[1,2]。我国从20世纪80年代开始，主要采用综合评判法对全国水土流失状况进行了3次遥感调查；2010-2012年全国第1次水利普查水蚀区水土流失专项调查中，采用了抽样调查与CSLE模型法，但该方法中部分参数因子的区域率定与计算仍需进一步研究[3-5]。目前的常规区域监测大多针对整个区域进行全面监测，国家层面的遥感调查周期基本上是5 a甚至更长时间，省级层面开展一次全面系统的水土流失调查大约需要2 a时间。受自然因素和人为因素的共同作用，区域水土流失分布呈现空间异质性。对于监测区域内人为干扰较弱、不易发生侵蚀、环境相对稳定不易发生变化的区域采取同样的监测方式，不但耗时长，而且需要动用大量人力物力。

本研究选择济源市作为试验研究范围，提出空间异质性区域划分，将受人为活动干扰、容易发生水土流失的区域，划定为水土流失易发区。在划定水土流失易发区与非易发区的基础上，进一步将易发区域内水土流失状况容易发生变化、人为干扰强烈、自然和社会环境易发生改变的区域划分为水土流失易变区。以区域划分成果和近年最新水土流失遥感监测成果为前提，通过3S技术和常规监测技术的有机结合，采用单元抽样等统计学方法，研究在不同抽样方案条件下，以抽样单元推算区域、以区域推算县域尺度，快速准确获取区域乃至县域水土流失情况的可行性及推算精度。

1 研究区域概况

济源市位于河南省西北部，黄河北岸，与山西省为邻，处于东经112°02′-112°45′、北纬34°54′-35°17′，面积为1 908 km2。地势呈西高东低、北高南低之势，全市东部海拔最低，为120 m，北部为太行山脉，鳌背山、斗顶峰海拔高度分别为1 930、1 955 m。平原、台地、丘陵、低山、中山地貌均有分布，其中山丘区面积占85.4%。暖温带大陆性季风气候，多年平均气温14.3 ℃，极端最低气温-20 ℃，极端最高气温43.4 ℃，大于10 ℃的活动积温4 539.6 ℃，无霜期平均223 d，多年平均日 2 375.4 h，多年平均降水量641.7 mm，降水季节分布不均。境内河流属黄河水系，黄河自晋入豫，沿济源市市南界流径57 km，汇纳逢石河、大峪河、砚瓦河、蟒河、沁河等15条一级支流。

研究区北部太行山区地质剖面底部为片麻岩、片岩与石英岩，中部为石灰岩、夹页岩及部分砂岩，上部为厚层的石灰岩；西部低山丘陵区主要为沙页岩、灰岩及泥岩，东南部黄土丘陵区大部分为新生界第四系黄土及黄土状物质；东部平面地面普遍为第四系黄土物质覆盖。土壤类型主要为褐土、粗骨土、石质土。植被属温带落叶阔叶林，除少部分原始森林，多为松柏林和灌木丛。

2 研究内容与方法

2.1 水土流失基础数据获取

(1)地形、地貌数据。基于ArcGIS平台，以DEM为基础，提取海拔、起伏度、坡度等信息，构建研究区地形地貌、坡度等本底信息数据。

(2)林草覆盖、土地利用本底信息数据。基于ArcGIS、ENVI平台，以1998年TM遥感影像为基础，计算归一化植被指数、植被覆盖度，构建研究区植被覆盖度本底信息数据。根据解译标志，以1998年TM遥感影像为基础，通过人机交互解译方式，提取土地利用信息，构建研究区8个一级类(耕地、园地、林地、草地、居民点及工矿用地、交通运输用地、水域及水利设施用地、其他土地)，10个二级类的土地利用本底信息数据。

(3)水土流失本底信息数据。基于ArcGIS平台，综合分析1998年研究区地貌、坡度、土地利用、植被覆盖度等因子，依据《土壤侵蚀分类分级标准》(SL190-2007)，采用综合评判三因子法，判定水土流失强度，构建研究区水土流失面积、强度和空间分布等本底信息数据。

(4)水土流失变化情况数据。以2009年TM遥感影像为基础，采用上述方法，获取2009年林草覆盖、土地利用、水土流失等状况数据；对比1998年本底信息数据，分析获取林草植被覆盖度、土地利用、水土流失变化趋势、数量及空间分布等特征。

2.2 水土流失空间异质性区域划分

(1)水土流失易发区划分。水土流失易发区划分结合水土流失发生发展的特点，根据空间连续性、综合性、地带性、因地制宜性等原则，利用ArcGIS的空间数据的转换与处理，矢量数量、栅格数据的空间分析等功能，对获取的植被覆盖度、坡度、地貌类型、土地利用等水土流失影响因子进行叠置分析，并根据遥感影像的人为干扰迹象等因素进行人工修正，最终将受人为活动干扰、容易发生水土流失的区域，划定为水土流失易发区[6，7]。

(2)水土流失易变区划分。水土流失易变区划分原则与易发区划分原则相同，主要根据水土流失易发区内2期对比数据的土地利用、植被覆盖度、土壤侵蚀状况变化叠置分析，结合野外调查验证和当地土地利用规划、经济发展需求等进行综合分析进行划定。将易发区内水土流失状况容易发生变化、人为干扰强烈、自然和社会环境易发生改变的区域进一步划分水土流失易变区。

(3)异质性最小单元构建与优化。在异质性最小单元内，土地利用、坡度等级、植被覆盖度等级具有均一性，是本研究中的抽样最小单元。以1998年土壤侵蚀图斑数据为本底数据，根据土地利用、坡度和植被盖度分级数据，结合遥感影像、矢量化地形图，初步构建异质性最小单元，即初始异质性最小单元[8]。

对初始异质性最小单元，进行优化处理，采用逐级拆分或合并的方式，将大的抽样单元逐级拆分，小的单元逐级合并，最终形成面积分布集中于30～80 hm2的优化单元，即优化异质性最小单元。

2.3 水土流失监测抽样方案

本研究通过基础数据空间分析，划分水土流失易发区、水土流失易变区，在区域内，以初始异质性最小单元和优化异质性最小单元作为抽样样本单元，采用统计学样本容量计算公式计算抽样数量，确定抽样位置，通过转移概率矩阵外推获取水土流失易发区、易变区的区域水土流失数据，并与2009年判读数据进行对比，评价其区域尺度的精度。将区域水土流失数据、非参与抽样区域及本底水土流失数据叠加，外推获取县域尺度监测数据，并再次进行县域尺度的精度评价。

根据研究目的和要求[9-11]，共设定4种水土流失监测抽样方法及推算的对比研究方案，包括：基于易发区初始异质性最小单元抽样、基于易发区优化异质性最小单元抽样、基于易发区优化异质性最小单元不等概抽样、基于易变区优化异质性最小单元抽样。为保证抽样单元的随机性与代表性，每种方案均设置2～3组平行的抽样单元。

3 结果与分析

3.1 水土流失本底数据构建

水土流失本底数据构建内容包括地形地貌、林草覆盖、土地利用和水土流失本底信息等数据。

(1)地形地貌本底数据。济源市地貌类型主要有平原、台地、丘陵、低山、中山5种，呈现丘陵>低山>平原>中山>台地的状态，丘陵和低山为主。其中平原233.73 km2，占12.25%；台地44.96 km2，占2.36%；丘陵929.24 km2，占48.69%；低山区615.83 km2，占32.27%；中山区84.49 km2，占4.43%。地貌类型见图1。

图1 地貌类型分布Fig.1 Geomorphologic type map

济源市坡度整体坡度较缓，集中分布在0°～25°范围内，占县域面积的80.38%，但仍有部分地势较陡，大于25°的陡坡、急陡坡所占面积为374.52 km2，占19.62%。以平缓坡为主，面积为701.58 km2，占36.77%；其次为15°～25°陡坡和8°～15°斜坡，分别占20.61%和16.73%。坡度分级图见图2。

图2 坡度分级分布Fig.2 Slope Degree classified Map

(2)林草覆盖本底数据。 与1998年相比， 2009年济源市林草面积增加53.55 km2。其中，中高覆盖明显增加，从1998年的451.31 km2增加到2009年的615.91 km2，增加了164.60 km2，所占林草面积比例增加了12.90%；其次为高覆盖，增加了38.43 km2，变化幅度为27.28%；中覆盖与中低覆盖面积明显减少，分别从1998年358.20 km2和71.85 km2减少到2009年的242.29 km2和37.65 km2，减少幅度分别为32.36%、47.60%；低覆盖变化不显著，1998-2009年仅存在0.63 km2的差异。2期林草覆盖见图3、图4。

图3 1998年林草覆盖度分布Fig.3 Map of forest cover in 1998

图4 2009年林草覆盖度分布Fig.4 Map of forest cover in 2009

(3)土地利用本底数据。与1998年相比，2009年土地利用结构仍以耕地、林地为主，新增园地、旱梯田。土地利用一级类中，耕地面积减少，其他土地利用类型所占面积均呈增加趋势。2期土地利用见图5、图6。

图5 1998年土地利用分布Fig.5 Land use map in 1998

图6 2009年土地利用分布Fig.6 Land use map in 2009

(4)水土流失变化数据。1998-2009年，县域总体土壤侵蚀呈下降态势，侵蚀强度降低，水土流失面积减少。水土流失减少面积达201.37 km2，占全县面积10.55%。其中：中度侵蚀面积变化最大，减少260.01 km2；轻度侵蚀变化次之，增加72.09 km2；强烈侵蚀减少13.46 km2，占其原面积的54.47%。

全县有456.53 km2的土壤侵蚀强度等级发生变化，占全县面积的23.92 km2，其中445.78 km2侵蚀强度降低，10.75 km2侵蚀强度升高。土壤侵蚀变化最为显著的为中度侵蚀变化为轻度侵蚀，集中分布在西南丘陵区，尤其是小浪底水库水利工程周边，耕地退耕，林地面积增加，林地植被覆盖增加的区域较多；在北部丘陵与中低山地貌过渡带，耕地减少，林地增加的现象也很显著。其次为轻度侵蚀变化为微度侵蚀，主要分布在城镇周边的耕地，坡耕地整治为梯田，侵蚀强度降低，涉及邵原镇、王屋乡、轵城乡、克井乡；部分变化在中低山区呈零散分布，主要是由于植被覆盖度变化；同时受小浪底水利工程影响，部分耕地、林地变成水域及水利设施用地，侵蚀强度降低。中度侵蚀变为微度侵蚀，集中分布在小浪底水库周边，同样是由于土地利用发生变化导致。其他侵蚀强度等级间的转换较小，多小于5 km2，呈点状分布。2期水土流失见图7、图8。

图7 1998年土壤侵蚀强度分布Fig.7 Soil erosion intensity map in 1998

图8 2009年土壤侵蚀强度分布Fig.8 Soil erosion intensity map in 2009

3.2 水土流失空间异质性区域划分

(1)结合济源市实际状况，划分水土流失易发区面积1 167.2 km2，占全县面积比为61.2%。分布高程范围200～600 m，地貌以丘陵为主，此外还包括平原向丘陵的过渡带、丘陵向低山的过渡地带。土地利用以旱地、林地为主，林地植被覆盖度以中覆盖为主。主要涉及邵原镇、王屋乡、下冶乡、大峪乡、坡头乡、轵城乡、承留乡、思礼乡、克井乡、五龙口等10个乡镇。划分水土流失易变区面积为622.3 km2，占水土流失易发区面积的53.3%。水土流失易发区、易变区分布见图9、图10。

图9 水土流失易发区分布Fig.9 Map of soil erosion prone area

图10 水土流失易变区分布Fig.10 Map of soil erosion volatility area

(2)水土流失易发区构建初始异质性最小单元6 863 个。其中耕地2 435 个，林地4 364 个，草地27 个，园地37 个。最大单元面积为87 hm2，最小单元为3 hm2，平均单元面积为25 hm2。异质性最小单元面积分布集中在10～40 hm2，以20～30 hm2面积区间内数量最大，所占比例为37.2%。初始异质性最小单元分布见图11。

图11 初始异质性最小单元分布Fig.11 Map of initial heterogeneity minimum element

(3)在初始最小异质性单元基础上，分别对面积大于100、90和80 hm2的单元进行逐级拆分优化，得到优化异质性最小单元3 862 个，最大单元面积为85.9 hm2，最小单元面积为6.3 hm2，平均面积为44.6 hm2。面积分布在30～80 hm2范围内的异质性最小单元所占比例为80.6%。优化异质性最小单元分布见图12。

图12 优化异质性最小单元分布Fig.12 Map of Optimal heterogeneity minimum element

3.3 基于易发区初始异质性最小单元抽样方案及其推算

以水土流失易发区作为抽样监测对象，初始异质性最小单元为抽样单元，计算抽样数量，确定抽样位置，通过抽样单元内2期侵蚀变化，建立转移矩阵，通过面积外推法获取易发区外推数据，并与判读数据对比，进行区域尺度的精度评价。将易发区外推数据、非参与抽样数据与非易发区本底数据叠加分析，获取县域尺度监测数据，并进行县域尺度的精度评价。

(1)抽样数量。在耕地、林地、草地3种土地利用类型中按各级侵蚀强度进行抽样。济源市易发区抽样单元总计245个，其中耕地108个，林地76个，草地61个。草地抽样比例最大，为26.6%，耕地与林地的抽样比例分别为14.8%、7.6%。

3组抽样单元总面积分别为145.6、144.0、142.8 km2，抽样单元总面积比较接近，组间相对稳定。经计算分析，3组抽样精度分别为98.9%、94.4%、84.6%，拟合精度分别为98.8%、97.9%、97.3%，表明3组抽样单元均能很好地代表样本总体。

(2)区域尺度精度评价。以微度侵蚀相对误差作为抽样外推精度评价指标，以各级侵蚀强度相对误差面积加权平均作为拟合度指标。由表1可知，3组抽样所对应的抽样外推精度为97.3%、73.6%、85.7%；拟合精度为94.8%、87.6%、92.4%。第1次抽样外推精度较高，与第2次、第3次抽样结果有一定差异，产生这种差异的原因是抽样单元具有均一性和2期对比水土流失数据空间分布的规律性。

表1 基于易发区初始单元外推区域数据Tab.1 Extrapolation of regional data based on the initial element of the prone area

(3)县域尺度精度评价。将水土流失易发区数据与居民点水域及工矿用地、交通运输用地、水域及水利设施用地等非易发区本底数据进行叠加分析，获取县域监测成果，并与2009年判读数据对比分析进行精度分析。由表2可以看出，3次抽样的相对误差均小于5%，推算精度大于95%，绝对误差小于2%，各级侵蚀强度加权平均精度大于95%。表明该水土流失监测方法进行水土流失估算是可行的，监测成果的精度可以满足要求。

表2 县域尺度成果数据Tab.2 County scale data

3.4 基于易发区优化异质性最小单元抽样方案及其推算

该方法与基于易发区初始异质性单元抽样法相似，但抽样单元采用优化异质性最小单元，减少了抽样单元面积不均一所造成的影响。

(1)抽样数量。设置2组抽样，分别对面积大于90、80 hm2的初始单元进行拆分优化，优化处理后再进行基于易发区的优化异质性最小单元抽样。2组抽样单元数分别为245、241个，耕地抽样比例分别为13.3%、11.8%，林地抽样比例分别为7.0%、6.4%，草地抽样比例分别为24.7%、21.7%。

(2)区域尺度精度评价。由表3可以看出，第1组推算精度大于第2组，即对面积大于90 hm2的初始单元拆分后抽样的方法优于对80 hm2的初始单元进行拆分抽样的方法。分析原因，进一步拆分图斑引起单元数量的增加，但由于抽样单元的数量不会同步增加，进而导致抽样比例下降，非变化单元的抽中概率增加，因此反映在抽样结果上就是样本精度及推算精度均降低。

(3)县域尺度精度评价。根据县域尺度成果数据精度成果表4，第1组推算县域尺度数据与2009年判读数据相比，相对误差3.2%，外推精度为96.8%，绝对误差为1.7%；各级侵蚀强度加权平均误差3.5%，精度96.5%；水土流失误差3.5%。第2组数据相对误差4.7%，外推精度为95.3%，绝对误差为2.4%；各级侵蚀强度加权平均误差4.9%，精度95.1%；水土流失误差5.1%。研究表明，抽样精度的稳定性受异质性单元拆分面积的均一性的影响，拆分面积均一性高，抽样精度更加稳定；最优拆分面积在40～50 hm2，单元集中分布在30～80 hm2。

表3 基于易发区优化单元外推区域尺度数据Tab.3 Extrapolation of regional scale data based on theoptimization unit of the prone area

表4 县域尺度成果数据Tab.4 County scale data

3.5 基于易发区优化异质性最小单元不等概抽样方案及其推算

在水土流失易发区中的易变区和非易变区同时抽样，为使抽样单元尽可能多的落入侵蚀易变区，在抽样过程中加入限制条件。分别对易变区、非易变区抽样结果进行外推计算，并将2者外推结果进行求和处理，获得该抽样方法的外推结果，并与判读数据进行对比分析。

易变区抽样数目138个，占易变区总数的6.4%。其中耕地86个，抽样比例7.4%，林地52个，抽样比例5.1%。非易变区抽样数目偏少，仅34个，占非易变区总数的2.0%。

3组抽样研究结果表明，区域尺度精度分别为86.8%、87.9%和97.1%，县域尺度精度分别为79.9%、81.9%、98.8%。在3组抽样中，仅有1组抽样满足精度要求，整体水土流失面积判定误差较大，说明偏向于易变区的抽样方法已经干扰了侵蚀变化概率，不能保证抽样单元空间分布的随机性与均匀性。同时该方法在外推过程中，各级强度外推数据与判读数据误差较大，不能完全保证县域尺度外推数据的准确性。

3.6 基于易变区优化异质性最小单元抽样方案及其推算

基于易变区优化异质性最小单元抽样仅在水土流失易变区进行抽样及外推计算，在非易变区不进行抽样，最终将侵蚀易变区侵蚀抽样外推结果与非易变区侵蚀本底数据累加，构成水土流失易发区推算结果。该结果进而与非易发区数据进行汇总获得县域水土流失监测成果。

(1)抽样数量。由于草地单元数量不足以构成抽样样本，因此将草地单元与林地单元合并抽样。单元总数为1 129个，其中耕地517个，林草地612个。本次抽样共抽取单元125块，其中耕地71块，抽样比例为13.7%；林草地54块，抽样比例为8.8%。

(2)区域尺度精度评价。3组基于易变区的区域尺度外推精度分别为87.5%、87.0%、87.4%，平均外推精度87.3%。3组拟合精度分别为91.9%、93.0%、92.7%，平均拟合精度92.6%，拟合精度高于外推精度。总体而言，区域尺度精度满足要求，可以与非易发区数据进行汇总进行县域尺度数据推算。

(3)县域尺度精度评价。研究数据表5表明，第2组抽样数据外推精度最高，相对误差3.3%，外推精度为96.7%，绝对误差为1.9%。3组抽样相对误差均小于5%，外推精度大于95%，绝对误差小于2%，各级侵蚀强度加权平均精度大于90%。结果表明，该水土流失监测方法进行水土流失估算是可行的，监测成果的精度可以满足要求。

表5 县域尺度成果数据Tab.5 County scale data

4 结论与讨论

(1)县域尺度中最优异质性最小单元构建是本研究水土流失监测方法的关键。优化异质性最小单元的面积分布集中于30～80 hm2范围内为最优情况。

(2)基于易变区优化异质性最小单元抽样既考虑了抽样单元面积的均一性，又兼顾了人为干扰异质性。该方法各类精度、拟合精度均高于其他2种抽样方案，本研究推荐使用基于易变区优化异质性最小单元。

(3)基于易发区优化异质性最小单元不等概抽样，抽样数据与本底数据吻合程度偏低，干扰了侵蚀变化概率，外推过程中各级侵蚀强度外推数据与判读数据误差较大，不能完全保证外推数据的准确性，不建议采用该方法进行水土流失监测。

(4)易变区分布规律、易变区划分、变化概率的掌握均是根据不同时期内土壤侵蚀图谱的叠置统计分析所得，本研究基于1998年、2009年2期节点数据探索的结果，虽然验证了抽样方法的可行性，获得了水土流失易变区的分布规律，但仅靠2期节点数据还是有限的，因此下一步研究将增加1～2个节点或者更多，探索分布规律、变化趋势与特征。

(5)通过抽样单元监测数据推算县域尺度水土流失数据的方法，在河南省西部区域是可行的，但若将该监测方法推广至省内或省外其他区域，是否适用，有待进一步研究。

□

[1] 谢 云,赵 莹,张玉平.美国土壤侵蚀调查的历史与现状[J].中国水土保持,2013,(10):53-60.

[2] Morgan R P C, Quinton J N, Smith R E, et al. The Europeansoil erosion model ( EUROSEM) : a dynamic approach for predicting sediment transport from fields and small catchments[J].Earth Surface Processes and Landforms, 1998,(23):527-544.

[3] 郭索彦.水土保持监测理论与方法[M].北京：中国水利水电出版社,2010.

[4] 李智广,符素华,刘宝元.我国水力侵蚀抽样调查方法[J].中国水土保持科学,2012,10(1):77-81.

[5] 张新玉,鲁胜力,王 莹.我国水土保持监测工作现状及探讨[J].中国水土保持,2014,10(1):6-9.

[6] 陆欢欢.水土流失易发区划定研究[D]. 南京：南京农业大学,2013.

[7] 李 万. 关于自然地理区划原则的探讨[J]. 地理科学，1985,5(3):268-272.

[8] 罗志东.我国水土保持基础管理单元“水保斑”的认识与探索[J].中国水土保持科学,2015,13(8):127-131.

[9] 邵志强.抽样调查中样本容量的确定方法[J].统计与决策,2012,(22):12-14.

[10] 田 浩.水土保持监测抽样方法探讨[J].改革与开放, 2010,(5):97.

[11] 张国友.关于抽样调查中样本容量的确定[J].安徽理工大学学报,2003,5(1):30-33.