CN109579791B - 一种快速高效水土保持淤地坝调查与评估方法 - Google Patents

Abstract

本发明涉及一种快速高效水土保持淤地坝调查与评估方法，包括：初选淤地坝位置；获取高精度数字地形；求取各个栅格数据；确定需要进行地形特征判断的栅格周边的运算范围；根据地形特征判定淤地坝位置；判定淤地坝现状特征；计算坝后库容；淤积体积计算；淤地坝上游土壤侵蚀模数的计算。本发明先利用低分辨率的DEM数据确定可能出现淤地坝的位置，再利用航拍的方式构建高精度的DEM数字化地形，并利用数字化地形对淤地坝可能出现的区域进行分析，确定是否有淤地坝并利用数字化地形对确定有淤地坝的位置进行分析，最终判断出淤地坝的大小和淤积范围。在通过获取的淤积范围得到淤地坝的其他数据，并根据这些数据获取淤地坝上游水土保持的状态。

Description

一种快速高效水土保持淤地坝调查与评估方法

技术领域

本发明涉及一种快速高效水土保持淤地坝调查与评估方法，是一种水利设施的数字化调查方法。

背景技术

淤地坝是一种重要的水土保持工程措施，在世界范围内均有广泛的应用。据统计，仅在中国的黄河流域2003-2010年就建设了8180座淤地坝。这些淤地坝对河流流域的灌溉、防洪、水土保持、减少河流泥沙等方面做出了重要的贡献。然而，随着淤地坝使用时间的延伸，淤地坝及周围流域情况的不断变化，以及早期淤地坝建设的年代还没有成熟的数字技术，众多淤地坝的数字化信息或没有或残缺不全，难以实现数字化管理。如何补齐这些缺失的淤地坝数字化信息，需要进行缜密的调查，详细的了解各个淤地坝的状态，构建各个淤地坝的数字化模型，实现淤地坝的数字化管理。然而要实现这一详细的调查，现有的方式是采用人工调查的方式，对淤地坝进行现场测量，这需要耗费大量的人力物力。尽管有一些利用数字化地形图DEM进行淤地坝调查的方案，如对DEM进行三角化等方式，但这些方案都存在这样或那样的缺点，如分析过程复杂，没有充分的利用DEM数据，使用的计算机资源过多，使计算成本过高等。由于这些问题的存在，使利用计算机进行淤地坝调查停留在方案上，在实际中应用很少。

发明内容

为了克服现有技术的问题，本发明提出了一种快速高效水土保持淤地坝调查与评估方法。所述的方法利用河道的走向，快速构建数以千计的淤地坝数字模型，根据经年的记录，对淤地坝的淤积情况进行分析比较，最终得出淤地坝上游水土保持的状态。

本发明的目的是这样实现的：一种快速高效水土保持淤地坝调查与评估方法，其特征在于，所述方法的步骤如下：

步骤1，初选淤地坝位置：利用研究区域的低精度DEM数据，以汇水面积0.5km²为阈值提取流域数字河网，根据控制面积大小判断淤地坝可能的位置，运算的河网处的控制面积是否小于50km²且大于0.5km²，如果控制面积在此范围内，即为淤地坝的预选位置；

步骤2：获取高精度数字地形：利用流域数字河网和淤地坝预选位置制定无人机飞行路线，派出无人机按照飞行路线沿河网进行拍摄和航测，获取分辨率为5乘5米，高度相对误差小于1米的高精度DEM数据；

步骤3：求取各个栅格数据的步骤：根据高精度DEM数据，求取各个栅格DEM的坡向、坡度、曲率数据；

步骤4：确定需要进行地形特征判断的栅格周边的运算范围：根据淤地坝的坝长最大值Lmax和最小值Lmin，结合栅格尺寸l，以河网栅格的理论位置位于沟道正中的位置，考虑运行误差，同时提升淤地坝的识别精度，栅格运算的地形特征搜索判断半径取r=Lmax；

步骤5：根据地形特征判定淤地坝位置：对河网经过的栅格及其周边栅格进行扫描和坡向、坡度特征分析，分析范围按照栅格分辨率和坝长综合确定；

扫描过程为：以河道栅格（i，j）为中心，以n为搜索半径，n=int（r/l），由内向外检测各个栅格的坡向、坡度数据；i、j为河道中心网格坐标，以1为起始的自然数；

特征分析过程为：以坡向为主要判断条件，以坡度为辅助条件，计算比较（i-n，j-n）~（i+n，j+n）范围内栅格的地形变化特征；包括如下子步骤：

子步骤1：对（i，j）栅格周围的八个栅格的坡向数据进行比较，出现下述两种情况：

情况1：如果出现两个对边或者两个对角线边的坡向差值表现为相差在0附近，与之垂直或近似垂直方向两侧栅格的坡向差值相差在180附近，且前坡向差值表现为相差约180的一对栅格的坡度差值的绝对值大于河道平均比降；

情况2：如果出现有几个同侧栅格的坡向表现为0，另一侧的几个栅格的坡向相同；

如果不能确定出现上述特征，则以（i，j）栅格为中心逐步扩大搜索范围，直至达到搜索半径，如果还不能确定上述特征则停止在（i，j）栅格的搜索更换其他栅格进行搜索；

子步骤2：对（i，j）栅格周围八个栅格外的十六个栅格的坡向数据进行比较，出现下述两种情况：

情况1：坡向差值在0附近的栅格排列是否表现为直线或近似直线，其直线两侧的栅格坡向差值是否也表现在180附近；

情况2：坡向差值在0附近的栅格数是否在在增加，坡向相同的栅格是否沿一条直线延伸；

如果“否”或者不能确定上述特征，则以（i，j）栅格为中心继续逐步扩大搜索范围，直至达到搜索半径，如果还不能确定则停止（i，j）栅格为中心的搜索；如果“是”则进入下一步骤；

子步骤3：对所述直线两端的栅格的周围栅格进行坡向和坡度判断，是否同样表现出其沿直线方向栅格的特征；如果“否”或者不能确定则以正在搜索的栅格为中心逐步扩大搜索范围，直至达到搜索半径，如果还不能确定则停止搜索，如果“是”则继续沿直线扩大判断范围，直至沿直线方向栅格的特征消失，则所述直线所在栅格的位置为淤地坝的位置，直线的长度为淤地坝的长度；

其中，表现为情况1的是还未淤满的淤地坝，表现为情况2的已经淤满的淤地坝；

步骤6，判定淤地坝现状特征：根据所确定的淤地坝位置，提取淤地坝位置两侧栅格的地表曲率数据，则两侧曲率较大的位置即为沟道折线，折线延长的间距约等于淤地坝的坝长，以此利用曲率信息，检测并确认淤地坝的走向和长度；

步骤7，计算坝后库容：利用库区地形图，确定量算淤地面积的起点高程和等高距；分别量出各层等高线所包围的面积，即为该高程的淤地面积；用两层面积的平均值乘其等高距，得出层间库容；累加各层库容，得出各高程下的总库容，绘制水位—库容和水位—淤地面积关系曲线；

步骤8，淤积体积计算：通过水位-淤地面积、水位-库容曲线提取淤积深度，淤积深度和淤积面积相乘，得到淤积体积；以淤地坝为点，进行盆域分析，得到淤地坝控制面积；将有效库容减淤积体积，得到剩余有效库容；

步骤9，淤地坝上游土壤侵蚀模数的计算：土壤侵蚀模数（t/（km²·ya））=淤积体积/淤地坝控制面积/淤积年限。

本发明产生的有益效果是：本发明先利用低分辨率的DEM数据确定可能出现淤地坝的位置，再利用航拍的方式构建高精度的DEM数字化地形，并利用数字化地形对淤地坝可能出现的区域进行分析，确定是否有淤地坝并利用数字化地形对确定有淤地坝的位置进行分析，最终判断出淤地坝的大小和淤积范围。在通过获取的淤积范围得到淤地坝的其他数据，并根据这些数据获取淤地坝上游水土保持的状态。通过航拍可以节省大量人力物力，能够以最低的成本建立数以千计的淤地坝的数字化模型，并能够快速的确定淤地坝的位置，以及淤积状态和未来的变化，实现淤地坝的数字化管理。

附图说明

下面结合附图和实施例对本发明作进一步说明。

图1是本发明的实施例所述系统的结构示意图；

图2是本发明的实施例所述方法的流程图；

图3是本发明的实施例所述的9栅格分析示意图；

图4是本发明的实施例所述的25栅格分析示意图。

具体实施方式

实施例：

本实施例是一种快速高效水土保持淤地坝调查与评估方法。为实现完全自动化快捷的淤地坝调查，尽量减少人工干预，本实施例的基本思路是：首先利用免费公开的低精度DEM数据，进行河网分析，排除大量无用的信息，只对有河网的地域进行分析，同时根据淤地坝的特点进一步排除没有可能建立淤地坝的区域，确定可能出现淤地坝的位置，进一步减小研究区域。然后对已经限定的研究区域进行航拍，得到高精度的DEM数据。之后利用获取的高精度DEM数据对可能出现淤地坝的位置进行分析，最后确认淤地坝的位置和走向，以及淤地坝的淤积面积等信息。分析完全由计算机完成，这些计算机系统可以设置在地面，也可以直接安装在无人机中，在航拍的过程中直接分析完成淤地坝的识别，即在无人机飞行过程中传回的数据不是DEM数据，而是淤地坝的数据，这样可以减小无线通讯数据量，减少数据丢失。至于其DEM数据可以在无人机落地后再进行数据传输，获取全面完整的数据。

本实施例所述方法所使用的系统包括，如图1所示：航拍子系统和处理子系统，所述的航拍子系统包括：搭载摄影装置、激光测距仪、飞行卫星定位器，以及图像和信息实时传输装置的无人机，所述的处理子系统包括：图像和信息实时接收系统，所述的图像实时接收系统与数字地形辨别装置连接，所述的数字地形连接装置与数字河网分析装置连接，所述的数字河网分析装置与DEM数据求取装置连接，所述DEM数据求取装置与河道和坝址分析装置连接，所述的河道和坝址分析装置与淤地坝走向和长度分析装置连接。

所述的航拍系统是以无人机为主体的系统，在无人机上安装了各种设备并依据事先确定的航拍路线对淤地坝存在的河流流域进行航拍。

所述的摄影装置可以是光学数字拍摄装置，如照相机摄像机等，也可以雷达拍摄装置，如多普勒地形雷达，或者红外拍摄装置，或者同时安装几种拍摄装置，同时对地面进行光学、电磁波和红外拍摄，以获取精确的数字化地形。

在拍摄的过程中，同时利用激光测距仪对无人机的拍摄地面的实际高度进行测量，或直接获得地面的高程变化，或印证多普勒雷达的高度和地形数据。

卫星定位器可以是GPS定位，也可以是北斗定位，以及其他卫星定位器。通过卫星定位器可以获得准确的地形坐标，并与航拍图像结合，得到精确的数字化地形。

所述的处理子系统包括图像和信息实时接收系统，通过计算机进行相邻影像匹配，量测每个像素的视差，利用空中三角测量，自计算、识别同名像点，求得其像方坐标，经解算获得被摄体的物方空间坐标并生成DEM，最终将视差转换到相对或是绝对地面高程，再使用ARCMAP水文分析模块，对获取的DEM进行填洼处理，以获取无洼地DEM。

数字地形辨别装置对接收到的DEM数据进行流向分析，获取DEM流向数字信息。

数字河网分析装置通过DEM提取河网特征，根据这一特征找到淤地坝坝址。数字河网分析装置通过栅格转ASCII将得到的河网的DEM转换为.txt格式，命名为河网。进行坡向和坡度分析，分别得到DEM的坡度和坡向信息。再通过栅格转ASCII将坡度和坡向DEM分别转换为.txt格式，分别命名为坡度和坡向，因为河网、坡度、坡向三者中的栅格位置互相严格对应，所以得到的各单元格位置信息也是严格对应的。

DEM数据求取装置，即输入河网数据，其中有河网经过的单元格数字为1，没有河网经过的单元格数字为-9999，利用软件将河网单元格中数字为1的确定为有用单元格，并查到其相应的位置。输入坡向、坡度文件，确定坡向所需单元格信息，根据河网数据确定的有用单元格位置坐标，从第一个有用单元格开始，使用if语句开始判定，以河网栅格有用单元格为中心，以n=int（r/l）为半径，以坡向中两个等距单元格数值的差值绝对值是否在170到180之间为判定条件，如果差值在170到180之间，则进行坡度数据的判定。坡度判定相同位置单元格数值的差值，如果两侧栅格坡度值差值的绝对值大于河道平均比降，则说明此处为淤地坝，位置信息保存到文件中，再进行下一组坡向和坡度的判定，最终将确定为淤地坝的位置信息汇总到文件中，将输出的单元格进行标号，且定义一个和河网表同行同列的表格，将输出的单元格放入表格对应位置，定义表格前六行和初始河道表格前六行信息相同，其他非淤地坝位置定义为-9999。将所得文件导入gis，并且加载真实地形影像检验所得淤地坝位置结果。到此，淤地坝位置确定。

淤地坝走向和长度分析装置：根据上述计算确定的淤地坝位置，直接提取该位置两侧栅格的地表曲率数据，则两侧曲率较大的位置即为沟道折线，折线间的间距约等于淤地坝的坝长。

所述的分析子系统可以设置在地面，也可以安装在无人机上，使获取的DEM数据在无人机上即进行处理，直接得到淤地坝的位置信息和其他参数。其数据传输的方式可以是：实时传输淤地坝的数据和分析中的关键参数，但将分析前使用的DEM数据储存在无人机中，当无人机回到地面后再将所有数据通过线传至地面的系统中。这样可以节约大量的无线数据传输资源，以利用无人机的飞行控制和关键数据传输。

本实施例所述方法的具体步骤如下，流程图见图2所示：

步骤1，初选淤地坝位置：利用研究区域的低精度DEM数据，以汇水面积0.5km²为阈值提取流域数字河网，根据控制面积大小判断淤地坝可能的位置，运算的河网处的控制面积是否小于50km²且大于0.5km²，如果控制面积在此范围内，即为淤地坝的预选位置。

本步骤的作用是缩小航拍范围，为航拍路线建立数据基础。首先利用低精度的DEM数据进行河网分析。在河网分析过程中将河流流域以外的地区排出，这样就可以排出一大部分地域数据，减少了计算量。

在河网分析过程中，通过公开的DEM数据能够很容易的得到河道的数字模型和所有信息，包括：深度、流向、比降等，不论在河道上是否构建了淤地坝，其河道数据可以显示出河道的这些特征。获得这些河道信息后就可以根据淤地坝的特点排出不可能建立淤地坝的区域。淤地坝只会建立在河道支流上，不可能出现在河道干流，因此可以将河道干流排除。淤地坝的控制面积一般小于50km²，如果大于50km²就是水库了，因此，可以将控制面积大于50km²的水域排除。至于控制面积最小应当是多少可以排除，主要根据调查工作要求而定，例如调查的出资人是当地政府，即调查的淤地坝主要是当地政府出资兴建的淤地坝，这些淤地坝的控制面积一般都大于0.5km²，因此可以将控制面积的下限定位0.5km²，小于这个控制面积的淤地坝通常由当地农民或其他出资人建立，不在调查范围内，因此可以排除。

所述的低精度DEM数据是指可以免费获得的DEM数据。这种免费的EDM数据的分辨率和精度具有一定的水平，能够达到粗略确定淤地坝位置的目的。例如：2015年发布的ASTER GDEM V2数据，其空间分辨率达到30乘30米，高度误差小于20m，具有95%置信度

通过免费的DEM数据，可以自动的获得研究区域的河网，以及河网的形态和周围环境，因此，通过DEM数据的处理，就能够在计算机系统中分辨研究区域的流域、主干流和支流，可能出现淤地坝的控制范围等数据，通过这些数据，在计算机系统中就能够确定淤地坝的预选位置。

步骤2：获取高精度数字地形：利用流域数字河网和淤地坝预选位置制定无人机飞行路线，派出无人机按照飞行路线沿河网进行拍摄和航测，获取分辨率为5乘5米，高度相对误差小于1米的高精度DEM数据。

航拍路线的制定十分关键，如果采用大型的无人机并搭载高档的仪器设备，可以在较高的高度进行大范围航拍，能够快速的大面积拍摄，在很短的时间内完成拍摄。但高空拍摄无人机使用成本较高，并且受天气影响较大。可以使用小型无人机进行低空拍摄，甚至可以使用多架无人机进行星座飞行，各架无人机分别携带不同的仪器设备，相互协作完成拍摄。这种方式虽然设备成本较低，但由于是低空拍摄，需要较长的飞行路线才能囊括较大的区域。因此，拍摄路线的制定除上述因素之外，还要根据项目的大小、资金情况等因素确定。

步骤3：求取各个栅格数据的步骤：根据高精度DEM数据，求取各个栅格DEM的坡向、坡度、曲率数据。

本步骤可以通过软件系统完成，无需人工干预。

步骤4：确定需要进行地形特征判断的栅格周边的运算范围：根据淤地坝的坝长最大值Lmax和最小值Lmin，结合栅格尺寸l，以河网栅格的理论位置位于沟道正中的位置，考虑运行误差，同时提升淤地坝的识别精度，栅格运算的地形特征搜索判断半径取r=Lmax。

现有的淤地坝最长的在1千米以上，但很少有超过2千米的淤地坝。至于最短的淤地坝如前所述根据投资方的要求而确定，例如，投资方是地方政府，一般不会修筑10米以下的淤地坝。

特征判断半径的选取目的在于限制判断的范围，以节约成本。较小半径的精细判断有助于得到精细的数据，但成本过高，而较大半径的粗略判断虽然能够节约成本，但数据的误差可能影响整体判断的正确性，因此把握判断半径的大小也是十分重要的。

步骤5：根据地形特征判定淤地坝位置：对河网经过的栅格及其周边栅格进行扫描和坡向、坡度特征分析，分析范围按照栅格分辨率和坝长综合确定。

本步骤的思路是：在最可能出现淤地坝的河道中心位置选择一个栅格（i，j）为中心，逐个比较周围栅格与栅格（i，j）的坡向和坡度。这里有一个前提，由于淤地坝都是土石的重力坝，水平投影均为直线，不会是弯曲的或者是折线形，同时，坝顶的宽度不会超过10米。根据这一特点，本步骤判断的过程是：如果周围栅格出现几个栅格有坡面直线分界线变化的倾向，在经过判断分界线两侧坡度的判断，就确定淤地坝的位置。

判断的可以先从九宫判断开始，如图2所示，再进一步扩大为25宫，如图3所示，再扩大为36宫等，直至n´n宫。

坡面直线分界线变化倾向有两种，一种是坝后尚未淤满，另一种是坝后已经淤满。前者的坡向相差180，而后者的直线一侧坡面均表现为零。这两个特点在判断中都会表现得十分明显，在经过河道比降的比较，就能够很方便快捷的得到淤地坝的准确位置。

扫描过程为：以河道栅格（i，j）为中心，以n为搜索半径，n=int（r/l），由内向外检测各个栅格的坡向、坡度数据；i、j为河道中心网格坐标，以1为起始的自然数。

根据DEM数据能够很容易的得到河道的数据，通过河道的数据，就可以找到河道的中心栅格，并对中心栅格周围的栅格进行扫描分析。

特征分析过程为：以坡向为主要判断条件，以坡度为辅助条件，计算比较（i-n，j-n）~（i+n，j+n）范围内栅格的地形变化特征；包括如下子步骤：

子步骤1：对（i，j）栅格周围的八个栅格的坡向数据进行比较，出现下述两种情况：

情况1：如果出现两个对边或者两个对角线边的坡向差值表现为相差在0附近，与之垂直或近似垂直方向两侧栅格的坡向差值相差在180附近，且前坡向差值表现为相差约180的一对栅格的坡度差值的绝对值大于河道平均比降。

情况2：如果出现有几个同侧栅格的坡向表现为0，另一侧的几个栅格的坡向相同。

如果不能确定出现上述特征，则以（i，j）栅格为中心逐步扩大搜索范围，直至达到搜索半径，如果还不能确定上述特征则停止在（i，j）栅格的搜索更换其他格栅进行搜索。

子步骤2：对（i，j）栅格周围八个栅格外的十六个栅格的坡向数据进行比较，出现下述两种情况：

情况1：坡向差值在0附近的栅格排列是否表现为直线或近似直线，其直线两侧的格栅坡向差值是否也表现出在180附近；

情况2：坡向差值在0附近的栅格数是否在在增加，坡向相同的栅格是否沿一条直线延伸；

如果“否”或者不能确定上述特征，则以（i，j）栅格为中心继续逐步扩大搜索范围，直至达到搜索半径，如果还不能确定则停止（i，j）栅格为中心的搜索；如果“是”则进入下一步骤；

子步骤3：对所述直线两端的栅格的周围栅格进行坡面和坡度判断，是否同样表现出其沿直线方向格栅的特征；如果“否”或者不能确定则以正在搜索的栅格为中心逐步扩大搜索范围，直至达到搜索半径，如果还不能确定则停止搜索，如果“是”则继续沿直线扩大判断范围，直至沿直线方向栅格的特征消失，则所述直线所在栅格的位置为淤地坝的位置，直线的长度为淤地坝的长度。

本步骤1中的情况1为淤地坝没有淤满的情况，而情况2则是淤地坝已经淤满的情况。

步骤6，判定淤地坝现状特征：根据所确定的淤地坝位置，提取淤地坝位置两侧栅格的地表曲率数据，则两侧曲率较大的位置即为沟道折线，折线延长的间距约等于淤地坝的坝长以此利用曲率信息，检测并确认淤地坝的走向和长度。

提取曲率数据能够找到淤地坝的淤积体与坝体之间的分界线，以此确定淤地坝后的淤积程度，如果已经淤满，这一分界线则没有了，但通过前述的判断已经得知淤地坝是否已经淤满，因此如果是情况2，就不需要进行本步骤了。

步骤7，计算坝后库容：利用库区地形图，确定量算淤地面积的起点高程和等高距；分别量出各层等高线所包围的面积，即为该高程的淤地面积；用两层面积的平均值乘其等高距，得出层间库容；累加各层库容，得出各高程下的总库容，绘制水位—库容和水位—淤地面积关系曲线。

利用等高线法绘制坝高—淤地坝面积曲线与坝高—库容曲线。利用库区地形图，确定量算淤地面积的起点高程和等高距（等高距一般选1m~3m）。累加各层库容，得出各高程下的总库容，绘制水位—库容和水位—淤地面积关系曲线。

步骤8，淤积体积计算：通过水位-淤地面积、水位-库容曲线提取淤积深度，淤积深度和淤积面积相乘，得到淤积体积；以淤地坝为点，进行盆域分析，得到淤地坝控制面积；将有效库容减淤积体积，得到剩余有效库容。

淤积体积：通过水位-面积、水位-库容曲线提取淤积深度，淤积深度和淤积面积相乘，得到淤积体积。

计算过程为：分别量出各层等高线所包围的面积，即为该高程的淤地面积。用两层面积的平均值乘其等高距，得出层间库容。相邻水层之间的体积可按下式进行计算：

式中：

—两水层间的体积，单位为立方米（m³）；

—两水层间的高距，单位为米（m）；

、

—分别为相邻上下两水层的面积，单位为平方米（m²）

步骤9，淤地坝上游土壤侵蚀模数的计算：土壤侵蚀模数（t/（km²·ya））=淤积体积/淤地坝控制面积/淤积年限。

淤地坝的上游水土保持状态可以根据淤地坝的淤积情况进行判断，这是淤地坝使用过程中对上游水土保持情况的有力证据，体现出上游的水土保持措施作用。

最后应说明的是，以上仅用以说明本发明的技术方案而非限制，尽管参照较佳布置方案对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案（比如调查的具体方式、各种淤地坝的形式、步骤的先后顺序等）进行修改或者等同替换，而不脱离本发明技术方案的精神和范围。

Claims (1)

1.一种快速高效水土保持淤地坝调查与评估方法，其特征在于，所述方法的步骤如下：

步骤1，初选淤地坝位置：利用研究区域的低精度DEM数据，以汇水面积0.5km²为阈值提取流域数字河网，根据控制面积大小判断淤地坝可能的位置，运算的河网处的控制面积是否小于50km²且大于0.5km²，如果控制面积在此范围内，即为淤地坝的预选位置；

步骤2：获取高精度数字地形：利用流域数字河网和淤地坝预选位置制定无人机飞行路线，派出无人机按照飞行路线沿河网进行拍摄和航测，获取分辨率为5乘5米，高度相对误差小于1米的高精度DEM数据；

步骤3：求取各个栅格数据的步骤：根据高精度DEM数据，求取各个栅格DEM的坡向、坡度、曲率数据；

步骤4：确定需要进行地形特征判断的栅格周边的运算范围：根据淤地坝的坝长最大值Lmax和最小值Lmin，结合栅格尺寸l，以河网栅格的理论位置位于沟道正中的位置，考虑运行误差，同时提升淤地坝的识别精度，栅格运算的地形特征搜索判断半径取r=L_max；

步骤5：根据地形特征判定淤地坝位置：对河网经过的栅格及其周边栅格进行扫描和坡向、坡度特征分析，分析范围按照栅格分辨率和坝长综合确定；

扫描过程为：以河道栅格（i，j）为中心，以n为搜索半径，n=int（r/l），由内向外检测各个栅格的坡向、坡度数据；i、j为河道中心网格坐标，以1为起始的自然数；

特征分析过程为：以坡向为主要判断条件，以坡度为辅助条件，计算比较（i-n，j-n）~（i+n，j+n）范围内栅格的地形变化特征；包括如下子步骤：

子步骤1：对（i，j）栅格周围的八个栅格的坡向数据进行比较，出现下述两种情况：

情况1：如果出现两个对边或者两个对角线边的坡向差值表现为相差在0附近，与之垂直或近似垂直方向两侧格栅的坡向差值相差在180附近，且前坡向差值表现为相差约180的一对格栅的坡度差值的绝对值大于河道平均比降；

情况2：如果出现有几个同侧栅格的坡向表现为0，另一侧的几个格栅的坡向相同；

如果不能确定出现上述特征，则以（i，j）栅格为中心逐步扩大搜索范围，直至达到搜索半径，如果还不能确定上述特征则停止在（i，j）栅格的搜索更换其他栅格进行搜索；

子步骤2：对（i，j）栅格周围八个栅格外的十六个栅格的坡向数据进行比较，出现下述两种情况：

情况3：坡向差值在0附近的栅格排列是否表现为直线或近似直线，其直线两侧的栅格坡向差值是否也表现出在180附近；

情况4：坡向差值在0附近的栅格数是否在在增加，坡向相同的栅格是否沿一条直线延伸；

如果“否”或者不能确定上述特征，则以（i，j）栅格为中心继续逐步扩大搜索范围，直至达到搜索半径，如果还不能确定则停止（i，j）格栅为中心的搜索；如果“是”则进入下一步骤；

子步骤3：对所述直线两端的栅格的周围栅格进行坡面和坡度判断，是否同样表现出其沿直线方向格栅的特征；如果“否”或者不能确定则以正在搜索的栅格为中心逐步扩大搜索范围，直至达到搜索半径，如果还不能确定则停止搜索，如果“是”则继续沿直线扩大判断范围，直至沿直线方向栅格的特征消失，则所述直线所在栅格的位置为淤地坝的位置，直线的长度为淤地坝的长度；

其中，表现为情况1、3的是还未淤满的淤地坝，表现为情况2、4的是已经淤满的淤地坝；

步骤6，判定淤地坝现状特征：根据所确定的淤地坝位置，提取淤地坝位置两侧栅格的地表曲率数据，则两侧曲率较大的位置即为沟道折线，折线延长的间距约等于淤地坝的坝长以此利用曲率信息，检测并确认淤地坝的走向和长度；

步骤7，计算坝后库容：利用库区地形图，确定量算淤地面积的起点高程和等高距；分别量出各层等高线所包围的面积，即为该高程的淤地面积；用两层面积的平均值乘其等高距，得出层间库容；累加各层库容，得出各高程下的总库容，绘制水位—库容和水位—淤地面积关系曲线；

步骤8，淤积体积计算：通过水位-淤地面积、水位-库容曲线提取淤积深度，淤积深度和淤积面积相乘，得到淤积体积；以淤地坝为点，进行盆域分析，得到淤地坝控制面积；将有效库容减淤积体积，得到剩余有效库容；

步骤9，淤地坝上游土壤侵蚀模数的计算：土壤侵蚀模数（t/（km²·ya））=淤积体积/淤地坝控制面积/淤积年限。

Images