CN103293285A - 一种在流域或区域尺度上的土壤侵蚀测定方法 - Google Patents

Abstract

本发明公开了一种在流域或区域尺度上的流域土壤侵蚀测定方法，包括以下步骤：S1、基础数据的归约化，通过野外GPS测量、遥感数据、数字地形图构建侵蚀沟道图和数字高程数据；S2、侵蚀地形的表达，获取水文地貌关系正确的DEM；S3、侵蚀坡长单元的提取，根据获取的部分区域的数字高程数据，计算坡度和坡长，分割单元坡长；S4、侵蚀学中L因子和S因子的提取，根据提取的坡度和坡长计算出L因子和S因子；S5、算法推广及应用，设计相应算法。本发明创造性的将水流分配策略的思想应用到单元坡长的提取与计算中，将单元坡长视为基础计算单元，对坡长按照水流分配策略进行分配，计算坡长、坡度，进而提取L因子和S因子，使得能够方便快捷的应用到较大区域或流域范围进行因子提取，进而测定土壤侵蚀，测量精度大大提高，且易推广使用。

Description

一种在流域或区域尺度上的土壤侵蚀测定方法

技术领域

本发明涉及土壤侵蚀和计算机技术领域，尤其是涉及一种在流域或区域尺度上进行坡长、坡度提取，进而测定流域尺度土壤侵蚀的方法。

背景技术

土壤侵蚀是指土壤或成土母质在外力(水、风)作用下被破坏剥蚀、搬运和沉积的过程。我国是世界土壤侵蚀最严重的国家之一，其范围遍及全国各地。土壤侵蚀的成因复杂，危害严重，主要侵蚀类型有水力侵蚀、风力侵蚀、重力侵蚀、冻融侵蚀和冰川侵蚀等。根据水利部遥感中心1990年调查统计，全国土壤侵蚀面积达492万km²，占国土面积的51％，其中水力侵蚀面积179万km2，风力侵蚀面积188万km²，冻融侵蚀面积125万km²。目前我国仍有水土流失面积356.92万平方公里，水土流失不但造成土壤退化，还引起江河淤积，水体污染，洪涝灾害加剧，导致生态恶化和区域贫困。因此水土流失已经成为我国实现可持续发展的严重障碍，国家对水土保持科技的需要比以往任何时候都迫切。重力侵蚀则不同程度地分别分布在以上3类土壤侵蚀分布范围内。现有技术中区域或流域土壤侵蚀测定方法多采用汇水面积、单位径流长度等值代替，而采用汇水面积、单位径流长度等值代替真值，误差较大，算法提取效率低，难以在大尺度范围进行快捷应用，不宜推广和应用。在坡面尺度上土壤侵蚀学的研究中，L因子和S因子是通过实测来获取的，长期以来，如何获得区域或流域尺度上的两个因子一直以来是侵蚀学研究者较困难解决的问题。

发明内容

本发明所要解决的技术问题是：提供一种在流域或区域尺度上的土壤侵蚀测定方法，创造性的将水流分配策略的思想应用到单元坡长的提取与计算中，将单元坡长视为基础计算单元，对坡长进行流向的分配，计算坡长、坡度，进而提取L因子和S因子，使得能够方便快捷的应用到较大区域或流域范围进行地形因子提取，进而基于统计模型(如：RUSLE)，计算土壤侵蚀。本方法对侵蚀计算精度大大提高，且易推广使用。

为解决上述技术问题，本发明的技术方案是：一种在流域或区域尺度上快速提取坡面长度，进而测定土壤侵蚀的方法，能够基于数字高程数据，快速高效的测定侵蚀地形因子，包括以下步骤：

S1、基础数据的归约化，通过野外GPS测量、遥感数据、数字地形图构建得到侵蚀沟道图和数字高程数据；

S2、侵蚀地形的表达，对步骤S1中获得的侵蚀沟道图和数字高程数据通过ANUDEM软件构建得到水文地貌关系正确的地形数据；

S3、侵蚀坡长单元的提取，根据步骤S1中获取的部分区域的数字高程数据，研究坡长的分配策略，计算坡度和坡长，采用图像分割方法，分割单元坡长；

其中，坡度计算公式： $\mathrm{\θ}=\left\{\begin{array}{cc}\arctan \left[\frac{\mathrm{Max}({\mathrm{DEM}}_c-{\mathrm{DEM}}_i)}{\mathrm{Length}}\right]& {\mathrm{DEM}}_c\≠{\mathrm{DEM}}_i\\ 0.1& {\mathrm{DEM}}_c={\mathrm{DEM}}_i\end{array}\right.$

θ是坡度；

Length是单元格长度，当水平或垂直方向时为单元格长度，对角方向时是

单元格长度；

DEM_i是单元格高程；

其中，坡长计算公式： ${\mathrm{\λ}}_{i,j}=\underset{x=0,y=0}{\overset{x=i,y=j}{\mathrm{\Σ}}}\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\λ}}_{x,y}$

λ_x，y是栅格内的单元坡长；

S4、侵蚀学中L因子和S因子的提取，根据提取的坡度和坡长计算出S因子和L因子；

其中，LS因子计算公式：

LS＝L·S

L＝(λ/22.13)^m

m＝β/(1+β)

β＝(sinθ)/[3·(sinθ)^0.8+0.56]

S＝10.8·sinθ+0.03   θ＜9％

S＝16.8·sinθ-0.5   θ≥9％

λ是坡长；m是坡长指数；β是坡度指数；θ是坡度。

采用了上述技术方案，本发明的有益效果为：本发明创造性的将水流分配策略的思想应用到单元坡长的提取与计算中，将单元坡长视为基础计算单元，对坡长进行流向的分配，计算坡长、坡度，进而提取L因子和S因子，将LS因子作为地形因子，用于经验模型RUSLE中，使得能够方便快捷的应用到较大区域或流域范围进行侵蚀测算，精度大大提高，且易推广使用。

附图说明

图1是本发明实施例的流程示意图；

图2是一现有算法的采样点统计图；

图3是另一现有算法的采样点统计图；

图4是本发明实施例的算法的采样点统计图；

图5是本发明实施例的沟谷提取图；

具体实施方式

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

本发明基于地理信息系统技术和数字地形分析技术，基于数字高程模型，应用水流分配策略的思想应用到坡长分配中，提取坡长、坡度，进而提取L因子和S因子。

如图1所示，一种在流域或区域尺度上的土壤侵蚀测定方法，包括以下步骤：

S1、基础数据的归约化，通过野外GPS测量、遥感数据、数字地形图构建得到侵蚀沟道图和数字高程数据。利用大比例尺航空摄影测量方法，结合高分辨率卫星遥感影像，高精度GPS地形测量数据，在1-2km²范围的微小流域内，获取高精度、高分辨率的DEM数据。

S2、侵蚀地形的表达，对步骤S1中获得的侵蚀沟道图和数字高程数据通过ANUDEM软件构建得到水文地貌关系正确的地形数据。这种DEM是澳大利亚国立大学的Hutchinson教授于20世纪80年代提出的，能合理表现地表高程的连续与突变，强调坡向转折，保证流水线的连续性和流域边界的准确，在国际上得到广泛应用。

S3、侵蚀坡长单元的提取，根据步骤S1中获取的部分区域的数字高程数据，研究坡长的分配策略，计算坡度和坡长。水流分配策略可以分为2种，一种是单一流向策略，其公式如下所示：

$S=\left\{\begin{array}{cc}\mathrm{Max}\left(\arctan \left(\frac{{\mathrm{DEM}}_c-{\mathrm{DEM}}_i}{\mathrm{Length}}\right)\right)& {\mathrm{DEM}}_c\≠{\mathrm{DEM}}_i\\ 0.1& {\mathrm{DEM}}_c={\mathrm{DEM}}_i\end{array}\right.$

i＝1，2，4，8，16，32，64，128

S为方向，i为方向，其中方向在计算过程中有相应编码代替，编码所对应的方向如下表1和表2所示，DEM_i为各个方向所代表的的高程，Length为单元坡长；

表1    表2

一种是多流向策略。单一流向策略即坡长累积沿着单一方向累积，但与实际地形不甚符合。多流向策略将坡长按比例进行多个方向分配，与实际地形吻合度较高。多流向计算公式如下：

$F_i=\frac{L_i\tan {\mathrm{\β}}_i}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{L}_i\tan {\mathrm{\β}}_i}$

i＝1，2，4，8，16，32，64，128

F_i为该方向上分配比例，L_i为该方向的单元坡长，β_i为该方向的水平夹角。

单一分配策略的算法实现较容易，且在沟道提取中的效果较好，因此采用单一分配策略对坡长进行提取。

在进行因子计算与提取前，需要设置有效的沟道截断阈值，可采用图像分割方法(C-V主动轮廓发)，将分割出来的小块单元，分析其坡长结果，根据不同小块单元的坡长变化特征，进而得到沟道截断阈值；也可以由领域专家来进行设置；

S4、侵蚀学中L因子和S因子的提取，根据提取的坡度和坡长计算出S因子和L因子。

LS因子计算公式如下：

LS＝L·S

L＝(λ/22.13)^m

m＝β/(1+β)

β＝(sinθ)/[3·(sinθ)^0.8+0.56]

S＝10.8·sinθ+0.03   θ＜9％

S＝16.8·sinθ-0.5   θ≥9％

λ是坡长；

m是坡长指数；

β是坡度指数；

θ是坡度。

以国内县南沟流域为例，图2和图3是两个现有算法的采样点统计图，图4是本发明实施例的算法的采样点统计图，可明显看出本算法的效果比现有算法的精度更高。

如图5所示，本实施例的沟谷提取图与地形图基本吻合。

综上所述，能够方便快捷的应用到较大区域或流域范围进行因子提取，测量精度大大提高，且易推广使用。

本发明不局限于上述具体的实施方式，本领域的普通技术人员从上述构思出发，不经过创造性的劳动，所作出的种种变换，均落在本发明的保护范围之内。

Claims (1)

1.一种在流域或区域尺度上快速提取坡面长度，进而测定土壤侵蚀的方法，其特征在于能够基于数字高程数据，快速高效的测定侵蚀地形因子，包括以下步骤：

S1、基础数据的归约化，通过野外GPS测量、遥感数据、数字地形图构建得到侵蚀沟道图和数字高程数据；

S2、侵蚀地形的表达，对步骤S1中获得的侵蚀沟道图和数字高程数据通过ANUDEM软件构建得到水文地貌关系正确的地形数据；

S3、侵蚀坡长单元的提取，根据步骤S1中获取的部分区域的数字高程数据，研究坡长的分配策略，计算坡度和坡长，采用图像分割方法，分割单元坡长；

其中，坡度计算公式： $\mathrm{\θ}=\left\{\begin{array}{cc}\arctan \left[\frac{\mathrm{Max}({\mathrm{DEM}}_c-{\mathrm{DEM}}_i)}{\mathrm{Length}}\right]& {\mathrm{DEM}}_c\≠{\mathrm{DEM}}_i\\ 0.1& {\mathrm{DEM}}_c={\mathrm{DEM}}_i\end{array}\right.$

θ是坡度；

Length是单元格长度，当水平或垂直方向时为单元格长度，对角方向时是

单元格长度；

DEM_i是单元格高程；

其中，坡长计算公式： ${\mathrm{\λ}}_{i,j}=\underset{x=0,y=0}{\overset{x=i,y=j}{\mathrm{\Σ}}}\underset{k=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\λ}}_{x,y}$

λ_x，y是栅格内的单元坡长；

S4、侵蚀学中L因子和S因子的提取，根据提取的坡度和坡长计算出S因子和L因子；

其中，LS因子计算公式：

LS＝L·S

L＝(λ/22.13)^m

m＝β/(1+β)

β＝(sinθ)/[3·(sinθ)^0.8+0.56]

S＝10.8·sinθ+0.03   θ＜9％

S＝16.8·sinθ-0.5   θ≥9％

λ是坡长；m是坡长指数；β是坡度指数；θ是坡度。

Images