CN103940974A - 基于gis的中尺度流域土壤侵蚀时空动态分析方法 - Google Patents

Abstract

本发明涉及一种基于GIS的中尺度流域土壤侵蚀时空动态分析方法，其具体过程为：把流域划分为若干栅格单元，作为计算单元；利用土壤侵蚀试验观测研究成果和认识，基于土壤侵蚀试验观测数据，提取关于土壤侵蚀的相关各因子；将上述各因子进行分析总结，建立侵蚀因子图集，完成对每个栅格单元土壤侵蚀量的估算；基于土壤侵蚀空间变化分析，完成对流域土壤侵蚀状况的评价和分析。本发明考虑了土壤侵蚀及其因子的空间不均匀性，因而将坡面模型与GIS集成，完成每个单元的运算，并适当考虑了物质传输和迁移，完成了对土壤流失的分布计算，输出了土壤侵蚀图。

Description

基于GIS的中尺度流域土壤侵蚀时空动态分析方法

技术领域

本发明涉及区域土壤侵蚀领域，尤其涉及一种基于GIS的中尺度流域土壤侵蚀时空动态分析方法。 

背景技术

土壤侵蚀是一个多时空尺度过程，土壤侵蚀的分布具有明显的时空变化特征，科学进行土壤侵蚀防治则需要研究其时空变化规律。土壤侵蚀模型包括经验模型，即从侵蚀产沙的基本成因出发，依据实际观测资料，采用数理统计分析的方法，建立坡面、流域或区域侵蚀产沙量与其主要影响因素之间的经验关系式。 

目前USLE模型具有较广的应用范围，其可以为土地利用规划者、水保专家估算不同土壤、坡长、坡度、作物管理条件下的土壤流失率，并确定可供选择的土地利用措施，以使土壤流失率限制在允许流失量以内。但是该模型的应用具有地域性，需要根据不同的地域特征进行修正和完善。 

现有技术中，已存在充分考虑地域差异的模型，其对全国进行分区建模，同时考虑径流、降雨、水土保持措施对侵蚀的影响，但该模型仍属于集总模型，无法反应区域内水土流失差异，也无法为水土保持措施配置提供参考，同时，模型基于一定的水土保持水平下统计分析结果，当水土保持措施发生变化时，需要重新测定。 

鉴于上述缺陷，本发明创作者经过长时间的研究和实践终于获得了本创作。 

发明内容

本发明的目的在于提供一种基于GIS的中尺度流域土壤侵蚀时空动态分析方法，用以克服上述技术缺陷。 

为实现上述目的，本发明提供一种基于GIS的中尺度流域土壤侵蚀时空动态分析方法，其具体过程为：其具体过程为： 

步骤1，把流域划分为若干栅格单元，作为计算单元； 

步骤2，利用土壤侵蚀试验观测研究成果和认识，基于土壤侵蚀试验观测数据，提取关于土壤侵蚀的相关各因子； 

步骤3，将上述各因子进行分析总结，建立侵蚀因子图集，完成对每个栅格单元土壤侵蚀量的估算； 

步骤4，基于土壤侵蚀空间变化分析，完成对流域土壤侵蚀状况的评价和分析； 

所述步骤4中采用的评价模型的基本形式如下述公式(1-1)、(1-2)、(1-3)、(1-4)、(1-5)所示，其为： 

A_c=f(R，K，LS,BET,g，w)                          (1-1) 

w=f(Vegcov,wind)                         (1-2) 

$A_t=\left(\underset{c=1}{\overset{n}{\mathrm{\Σ}}}(A_c\×s_c)\right)/\underset{c=1}{\overset{n}{\mathrm{\Σ}}}{s}_c---(1-3)$

$E_t=A_t\×S_t=\underset{c=1}{\overset{n}{\mathrm{\Σ}}}{A}_c\×s_c---(1-4)$

$E_r=100\×\frac{|E_t-S_o|}{S_o}---(1-5)$

式中：A_C为每个计算单元的年土壤侵蚀模数(t／km2·a)；R、K、LS、BET、g分别为降雨侵蚀因子、土壤可蚀性因子、坡度坡长因子、水土保持措施因子和沟蚀因子，w表示每个计算单元的风蚀强度等级；Vegcov、wind分别为植被盖度与风蚀范围；面蚀的计算将以CSLE为基础；At为计算流域的年平均土壤侵蚀模数(t／km2·a)，s_c为计算单元面积，S_t为流域总面积(km2)，Et为流域侵蚀产沙总量(104t)；S_o为水文观测的流域年输沙总量，Er为模型的计算相对误差(％)。 

进一步，在上述步骤2中，提取关于土壤侵蚀的相关各因子的过程为： 

步骤21，计算降雨侵蚀力因子R； 

步骤22，计算土壤可蚀性因子K； 

步骤23，计算地形因素，包括坡度因子S、坡长因子L、坡度坡长因子LS； 

步骤24，计算水土保持措施因子，包括：植物措施因子B；工程措施因子E；耕作措施因子T； 

步骤25，计算流域沟蚀因子g； 

步骤26，计算流域风蚀量。 

进一步，在上述步骤21中计算降雨侵蚀力因子R，利用区域气象观测资料，按照下述公式(2-1)、(2-2)和(2-3)计算， 

$R={\mathrm{\α}}_2{F}^{{\mathrm{\β}}_2}---(2-2)$

$F=\left(\underset{i=1}{\overset{12}{\mathrm{\Σ}}}{P}_i^2\right)P^{-1}---(2-3)$

式中，P为年平均降水量(mm)，Pi为第i月的平均降雨量(mm)；R为多年平均降雨侵蚀力(MJ·mm·hm^-2·h^-1·a^-1)；α₂、β₂为模型参数，α₂取值为0.3589，β₂取值为1.9462，F指数mm大小与年平均雨量P的季节分布有关，取值范围在P·12-1～P之间。 

进一步，在上述步骤22，计算土壤可蚀性因子K时，根据得到的土壤类型数据图为数据依据，按照下述公式(2-4)进行计算， 

$K=\{1.2+1.3\exp [1.256\mathrm{SAN}(1-\mathrm{SIL}/100\left]\right\}\×{\left(\frac{\mathrm{SIL}}{\mathrm{CLA}+\mathrm{SIL}}\right)}^{0.3}\×$

$[1.0-\frac{0.25C}{C+\exp (3.72-2.95C)}]\×[1.0-\frac{0.7\mathrm{SN}1}{\mathrm{SN}1+\exp (-5.51+22.9\mathrm{SN}1)}]---(2-4)$

式中，SAN、SIL和CLA是砂粒、粘粒和有机碳含量(％)，SN1=1-SAN／100；K值单位均为美制单位，t·acre·h／(100acre·ft·tonf·in)即吨·英亩·小时／(百英亩·英尺·吨力·英寸)。 

进一步，坡度因子算法如公式(2-5)， 

S=10.8sinθ+0.03 

S=16.8sinθ-0.05θ<5° 

S=21.91sinθ-0.965°≤θ<14°(2-5) 

坡长因子算法如下，θ≥14° 

m=0.2θ≤1° 

L=(λ／22.1)mm=0.31°<θ≤3°(2-6) 

m=0.43°<θ≤5° 

m=0.5θ>5° 

式中，S为坡度因子；θ为坡度；L为坡长因子；入为由DEM提取 的波长；m为坡长指数，22.1为22.1m标准小区坡长；α为坡度坡长指数，黄土高原取值0.5。 

进一步，植物措施因子B根据土地利用图和植物盖度图，利用NDVI根据下述公式(2-7)计算， 

$F=100\&times;\frac{\mathrm{float}\left(\mathrm{NDVI}\right)-{\mathrm{NDVI}}_{\min }}{{\mathrm{NDVI}}_{\max }-{\mathrm{NDVI}}_{\min }}---(2-8)$

式中，F为植被盖度，NDVI为归一化植被指数，NDVImin和NDVImax分别为无植被地区的NDVI值和植被良好覆盖地区的NDVI值；在现场测试基础上，结合遥感图像典型地类(沙地、高覆盖林地等)的采样，确定最高和最低植被盖度的NDVI阈值；将采样的结果与对应的NDVI值进行统计比较，当NDVI值≥NDVImax时，盖度可视为100％，当NDVI值≤NDVImin时，盖度可视为0％。 

进一步，工程措施因子E根据下述公式(2-9)计算， 

$E=1-\frac{S_t}{S}\&times;\mathrm{\&alpha;}-\frac{S_d}{S}\&times;\mathrm{\&beta;}---(2-9)$

式中，S_t为梯田面积，S_d为淤地坝控制面积，S为土地总面积，α，β分别为梯田和淤地坝的减沙系数，α，β分别确定为为0.764和1。 

进一步， 

计算流域沟蚀因子g，在无植被覆盖的黄土坡条件下，沟蚀因子g主要受降雨因素和地形即坡度的影响，参见公式(2-10)所示， 

G=1+1.60sin(θ-15)                          (2-10)。 

进一步，根据上述各因子的计算结果，通过图形计算代数方法将上述各因子图层相乘，完成CSLE模型，得到土壤侵蚀模数图， 

在获取土壤侵蚀模数后，对计算结果进行精度验证，采用确定性系数法进行精度验证，所述确定性系数法，用来描述计算值对目标值拟合精度的无量纲统计参数，一般取值在0-1之间确定性系数E的计算公式(2-11)所示， 

$E=1-\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(x_i-y_i)}^2}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(x_i-\stackrel{\&OverBar;}{x})}^2}---(2-11)$

式中，x_i和y_i分别为实测值和模型计算值，n为年份，为实测侵蚀量年平均值。 

进一步，在上述步骤23中，计算地形因素的过程为： 

步骤31，首先建立孤山川流域的DEM，采用原始数据是1：5万地形图，生成10m分辨率的H_C-DEM； 

步骤32，提取坡度和坡长，本发明对坡度的提取采用D8算法，若得到值为0，则将栅格值赋值为0.1°，以确保栅格之间径流连续，生成坡度图； 

坡长计算首先通过DEM进行填挖运算得到填挖后的DEM图，再计算其流向，通过运行流向图来计算流域的栅格坡长，同时结合坡度图确定坡面的径流原点及径流终点； 

在具体计算时，径流终点计算的标志用2.75°的坡度来实现，再根据计算的栅格坡长结合径流原点及径流终点来设置坡长，再次结合坡度图计算坡长指数，由坡长指数与流向图进行迭代计算得到流域的累计坡长，并最终获得流域坡长图； 

步骤33，根据上述式(2-5)、(2-6)计算坡度因子S，坡长因子L和坡度坡长因子LS。 

与现有技术相比较本发明的有益效果在于：本发明充分利用了GIS空间分析功能，其输出的结果(土壤侵蚀强度图)可以给出流域内各部位的侵蚀强度，因而既可满足沟道水利水保工程设计的要求，也可满足水土保持措施布设与规划的需求；同时，这种方法建立的模型较好的吸收了坡面土壤侵蚀模型的成果，特别是USLE、CSLE等模型，比较容易推广。至于单元的划分，可以是不规则的多边形(称为地块)，也可以是规则的栅格。 

本发明考虑了土壤侵蚀及其因子的空间不均匀性，因而将坡面模型与GIS集成，完成每个单元的运算，并适当考虑了物质传输和迁移，完成了对土壤流失的分布计算，输出了土壤侵蚀图。从图形的空间结构和各种统计特征看，基本可以反应该区域土壤侵蚀的时间和空间变化特征，反应土壤侵蚀与环境要素的关系，模型精度基本达到70％。 

附图说明

图1为本发明基于GIS的中尺度流域土壤侵蚀时空动态分析方法的流程图； 

图2为本发明获取相关因子的流程图； 

图3-1为本发明1975-2006年孤山川流域各站年平均降雨量； 

图3-2为本发明1975-2006年孤山川流域年降水量频率曲线图； 

图3-3为本发明1975年孤山川流域的年降雨侵蚀力图； 

图3-4为本发明1986年孤山川流域的年降雨侵蚀力图； 

图3-5为本发明1997年孤山川流域的年降雨侵蚀力图； 

图3-6为本发明2006年孤山川流域的年降雨侵蚀力图； 

图4-1为本发明为孤山川流域土壤类型图； 

图4-2为本发明孤山川流域K值图； 

图5-1为本发明孤山川流域S因子图； 

图5-2为本发明孤山川流域L因子图； 

图5-3为本发明孤山川流域LS因子图； 

图6-1为本发明孤山川流域1975的B值图； 

图6-2为本发明孤山川流域1986的B值图； 

图6-3为本发明孤山川流域1997的B值图； 

图6-4为本发明孤山川流域2006的B值图； 

图7为本发明孤山川流域不同土壤侵蚀模数所占面积百分比柱状图。 

具体实施方式

以下结合附图，对本发明上述的和另外的技术特征和优点作更详细的说明。 

本发明基于GIS的中尺度流域土壤侵蚀时空动态分析方法，依据分布式的经验模型，实现对大中流域土壤侵蚀的定量评价。 

本发明实验数据基于孤山川流域的资料和计算结果，分析器土壤侵蚀的时空变化特征，探讨其空间变化特征原因。 

该模型的计算过程为： 

步骤1，把流域划分为若干栅格单元； 

步骤2，利用土壤侵蚀试验观测研究成果和认识，基于土壤侵蚀试验观测数据，提取关于土壤侵蚀的各因子； 

步骤3，将上述各因子进行分析总结，建立侵蚀因子图集，完成对每个栅格单元土壤侵蚀量的估算； 

步骤4，基于土壤侵蚀空间变化分析，完成对流域土壤侵蚀状况的评价和分析。 

本发明步骤4中采用的评价模型的基本形式如下述公式(1-1)、(1-2)、(1-3)、(1-4)、(1-5)所示，其为： 

A_c=f(R，K，LS,BET,g，w)                          (1-1) 

w=f(Vegcov,wind)                          (1-2) 

$A_t=\left(\underset{c=1}{\overset{n}{\mathrm{\&Sigma;}}}(A_c\&times;s_c)\right)/\underset{c=1}{\overset{n}{\mathrm{\&Sigma;}}}{s}_c---(1-3)$

$E_t=A_t\&times;S_t=\underset{c=1}{\overset{n}{\mathrm{\&Sigma;}}}{A}_c\&times;s_c---(1-4)$

$E_r=100\&times;\frac{|E_t-S_o|}{S_o}---(1-5)$

式中：A_C为每个计算单元的年土壤侵蚀模数(t／km2·a)；R、K、LS、BET、g分别为降雨侵蚀力、土壤可蚀性、坡度坡长、水土保持措施因子和沟蚀因子，w表示每个计算单元的风蚀强度等级；Vegcov、wind分别为植被盖度与风蚀范围(地表形态)；面蚀的计算将以CSLE为基础；At为计算流域的年平均土壤侵蚀模数(t／km2·a)，s_c为计算单元面积，S_t为流域总面积(km2)，Et为流域侵蚀产沙总量(104t)；S_o为水文观测的流域年输沙总量，Er为模型的计算相对误差(％)。 

如上所述，在上述步骤2中，需要考虑降雨因素、土壤因素、地形因素、水土保持因素、沟道措施因素和风蚀因素；提取的因子包括、K、LS、BET、g分别为降雨侵蚀力因子R、土壤可蚀性因子、坡度坡长因子、水土保持措施因子和沟蚀因子。 

具体过程为： 

步骤21，计算降雨侵蚀力因子R；利用孤山川流域气象观测资料，按照下述公式(2-1)、(2-2)和(2-3)计算， 

$R={\mathrm{\&alpha;}}_2{F}^{{\mathrm{\&beta;}}_2}---(2-2)$

$F=\left(\underset{i=1}{\overset{12}{\mathrm{\&Sigma;}}}{P}_i^2\right)P^{-1}---(2-3)$

式中，P为年平均降水量(mm)，Pi为第i月的平均降雨量(mm)；R为多年平均降雨侵蚀力(MJ·mm·hm^-2·h^-1·a^-1)；α₂、β₂为模型参数，α₂取值为0.3589，β₂取值为1.9462，F指数mm大小与年平均雨量P的季节分布有关，取值范围在P·12-1～P之间。 

步骤22，计算土壤可蚀性因子K；根据得到的土壤类型数据图为数据依据，按照下述公式(2-4)进行计算， 

$K=\{1.2+1.3\exp [1.256\mathrm{SAN}(1-\mathrm{SIL}/100\left]\right\}\&times;{\left(\frac{\mathrm{SIL}}{\mathrm{CLA}+\mathrm{SIL}}\right)}^{0.3}\&times;$

$[1.0-\frac{0.25C}{C+\exp (3.72-2.95C)}]\&times;[1.0-\frac{0.7\mathrm{SN}1}{\mathrm{SN}1+\exp (-5.51+22.9\mathrm{SN}1)}]---(2-4)$

式中，SAN、SIL和CLA是砂粒、粘粒和有机碳含量(％)，SN1=1-SAN／100。K值单位均为美制单位，t·acre·h／(100acre·ft·tonf·in)即吨·英亩·小时／(百英亩·英尺·吨力·英寸)。 

步骤23，计算坡度因子S、坡长因子L、坡度坡长因子LS，根据生成的DEM，结合下述公式(2-5)、(2-6)计算， 

坡度因子算法如公式(2-5) 

S=10.8sinθ+0.03 

S=16.8sinθ-0.05θ<5° 

S=21.91sinθ-0.965°≤θ<14°(2-5) 

θ≥14° 

坡长因子算法如下， 

m=0.2θ≤1° 

L=(λ／22.1)m m=0.31°<θ≤3°(2-6) 

m=0.43°<θ≤5° 

m=0.5θ>5° 

式中，S为坡度因子；θ为坡度；L为坡长因子；入为由DEM提取的波长；m为坡长指数，22.1为22.1m标准小区坡长；α为坡度坡长指数，黄土高原取值0.5。 

步骤24，计算水土保持措施因子，本发明中水土保持措施因子包括：植物措施因子B；工程措施因子E；耕作措施因子T； 

其中，植物措施因子B根据土地利用图和植物盖度图，利用N DVI根据下述公式(2-7)计算， 

$F=100\&times;\frac{\mathrm{float}\left(\mathrm{NDVI}\right)-{\mathrm{NDVI}}_{\min }}{{\mathrm{NDVI}}_{\max }-{\mathrm{NDVI}}_{\min }}---(2-7)$

式中，F为植被盖度，NDVI为归一化植被指数，NDVImin和NDVImax分别为无植被地区的NDVI值和植被良好覆盖地区的NDVI值；在现场测试基础上，结合遥感图像典型地类(沙地、高覆盖林地等)的采样，确定最高和最低植被盖度的NDVI阈值。将采样的结果与对应的NDVI值进行统计比较，当NDVI值≥NDVImax时，盖度可视为100％，当NDVI 值≤NDVImin时，盖度可视为0％。 

NDVI数据计算的过程为： 

对于较小的流域(面积<=1000km2)，可利用TM数据计算NDVI，请参见公式(2-8)， 

$\mathrm{NDVI}=\frac{\mathrm{TM}4-\mathrm{TM}3}{\mathrm{Tm}4+\mathrm{TM}3}---(2-8)$

式中，NDVI为所求，归一化植被指数；TM3、TM4分别为遥感影像的近红外和红外波段。 

对于较大流域(面积>1000km2)，可直接获取。下载MODIS或SPOTVEGETATION的NDVI数据产品。 

工程措施因子E根据下述公式(2-9)计算， 

$E=1-\frac{S_t}{S}\&times;\mathrm{\&alpha;}-\frac{S_d}{S}\&times;\mathrm{\&beta;}---(2-9)$

式中，S_t为梯田面积，S_d为淤地坝控制面积，S为土地总面积，α，β分别为梯田和淤地坝的减沙系数，α，β分别确定为为0.764和1。 

计算耕作措施因子T：耕作措施因子值T主要是根据不同坡度条件下等高耕作减少土壤流失来确定，请参阅表1所示。 

表1孤山川流域不同坡度下耕作措施因子值 

步骤25，计算流域沟蚀因子g，在无植被覆盖的黄土坡条件下，沟蚀因子g主要受降雨因素和地形即坡度的影响，参见公式(2-10)所示， 

G=1+1.60sin(θ-15)                  (2-10)。 

步骤26，计算流域风蚀量； 

表2风力侵蚀的强度等级 

步骤27，土壤侵蚀评价；根据上述各因子的计算结果，通过图形计算代数方法将上述各因子图层相乘，完成CSLE模型，得到土壤侵蚀模数图。 

在上述步骤27之后，获取土壤侵蚀模数后，需对计算结果进行精度验证，可采用回归分析法和确定性系数法，其中，回归分析方法，分析模型计算值与实测值之间的相关系数，判断模型对实测值的拟合程度； 

所述确定性系数法，用来描述计算值对目标值拟合精度的无量纲统计参数，一般取值在0-1之间，确定性系数E的计算公式(2-11)所示， 

$E=1-\frac{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(x_i-y_i)}^2}{\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(x_i-\stackrel{\&OverBar;}{x})}^2}---(2-11)$

式中，x_i和y_i分别为实测值和模型计算值，n为年份，为实测侵蚀量年平均值。确定性系数的评价标准参见表3所示。 

表3确定性系数评定标准 

Table.2-3A.ssessment standards of certainty factor 

本实施例，对孤山川流域内各站点1975、1986、1997和2006年的各因子数值进行时空变化分析，现描述如下： 

(1)降雨侵蚀力，根据1975、1986、1997和2006年的月降雨量和降雨量侵蚀力计算模型分年度计算各站点的降雨侵蚀力，将得到的各站点的年降雨侵蚀力在ArcGIS下用IDW进行插值后，得到降雨侵蚀力插值表面，用孤山川流域的边界进行裁剪，得到孤山川流域的降雨侵蚀力空间插值表面，从中分析孤山川流域降雨侵蚀力的空间差异。请参见表4所示。 

表4孤山川流域降雨量侵蚀力特征统计表 

图3-1所示，为1975-2006年孤山川流域各站年平均降雨量； 

图3-2所示，为1975-2006年孤山川流域年降水量频率曲线图。 

图3-3、3-4、3-5、3-6分别为1975、1986、1997和2006年孤山川流域的年降雨侵蚀力图。 

(2)土壤可蚀性，根据上述(2-4)计算土壤可蚀性因子K，首先对收集到的府谷、准旗土壤类型图形进行数字化编辑及合并，得到孤山川流域类型图，再根据中国土种志、中国土壤、伊克昭盟土壤和榆林土壤等有关资料分析土壤类别，得到粘粒、粉粒、砂粒、有机质、有机含碳量等土壤理化性质，根据上述公式(2-4计算)。 

请参阅图4-1所示，其为孤山川流域土壤类型图； 

请参阅图4-2所示，其为孤山川流域K值图。 

(3)地形因素，步骤31，首先建立孤山川流域的DEM，采用原始数据是1：5万地形图，生成10m分辨率的H_C-DEM。 

步骤32，提取坡度和坡长，本发明对坡度的提取采用D8算法，若得到值为0，则将栅格值赋值为0.1°，以确保栅格之间径流连续，生成坡度图。 

坡长计算首先通过DEM进行填挖运算得到填挖后的DEM图，再计算其流向，通过运行流向图来计算流域的栅格坡长，同时结合坡度图确定坡面的径流原点及径流终点。在具体计算时，径流终点计算的标志用2.75°的坡度来实现，再根据计算的栅格坡长结合径流原点及径流终点来设置坡长，再次结合坡度图计算坡长指数，由坡长指数与流向图进行迭代计算得到流域的累计坡长，并最终获得流域坡长图。 

步骤33，根据上述式(2-5)、(2-6)计算坡度因子S，坡长因子L和坡度坡长因子LS。 

请参阅图5-1，5-2和5-3，分别为孤山川流域S、L和LS因子图。 

表5孤山川流域地形因子特征统计表 

(4)植物措施因子，本发明将农作物归为植物措施。根据上述式(2-7)计算植物措施因子。 

首先根据区域各年度的遥感影像图，根据式(2-7)计算出植被覆盖图，再结合区域植被覆盖图与土地利用图，得到植物措施因子措施图。 

请参阅图6-1、6-2、6-3、6-4所示，其为孤山川流域1975、1986、1997和2006年B值图。 

(5)工程措施因子，黄土地区工程措施主要有淤地坝、梯田、拦泥坝、涝池陡塘等。根据上述(2-9)计算工程措施因子。 

(6)耕作措施因子，本发明主要考虑等高耕作。 

(7)沟道措施因子，根据公式(2-10)计算沟蚀因子g。 

表6孤山川流域沟蚀因子值统计特征表 

(8)风力侵蚀，根据流域的风蚀范围图和植被覆盖图计算流域的风力侵蚀值，其中，风蚀范围图是将土地利用图中的沙地作为风蚀地段，将面积稍大点的几个图斑作为大片流动沙丘。 

表7风蚀评价标准 

在上述步骤4中，本发明对区域内的水蚀评价依据中国土壤流失方程进行，基于各侵蚀因子的计算结果，并通过地图代数计算方法完成模型的计算，得到区域水蚀模数图。 

请参阅图7所示，其为本发明孤山川流域不同土壤侵蚀模数所占面积百分比柱状图，将土壤侵蚀图中度以上侵蚀部分提取出来，再分别与地形、降雨侵蚀力图进行叠加，可发现土壤侵蚀与地形有较好的垂直对应关系。 

从流域的土地利用图可得出孤山川流域水土流失时空格局的土地利用关系，将区域中的年侵蚀图中度以上的侵蚀图与土地利用图叠加，可研究土壤侵蚀与土地利用的关系。 

本发明引入沟蚀系数，提高评价精度。 

表8考虑沟蚀后的精度分析表 

表9考虑风蚀、沟蚀后的精度分析表 

表9为考虑风蚀、沟蚀后的精度分析表，进一步提高分析精度。 

本发明充分利用了GIS空间分析功能，其输出的结果(土壤侵蚀强度 图)可以给出流域内各部位的侵蚀强度，因而既可满足沟道水利水保工程设计的要求，也可满足水土保持措施布设与规划的需求；同时，这种方法建立的模型较好的吸收了坡面土壤侵蚀模型的成果，特别是USLE、CSLE等模型，比较容易推广。至于单元的划分，可以是不规则的多边形(称为地块)，也可以是规则的栅格。 

本发明考虑了土壤侵蚀及其因子的空间不均匀性，因而将坡面模型与GIS集成，完成每个单元的运算，并适当考虑了物质传输和迁移，完成了对土壤流失的分布计算，输出了土壤侵蚀图。从图形的空间结构和各种统计特征看，基本可以反应该区域土壤侵蚀的时间和空间变化特征，反应土壤侵蚀与环境要素的关系，模型精度基本达到70％。 

以上所述仅为本发明的较佳实施例，对发明而言仅仅是说明性的，而非限制性的。本专业技术人员理解，在发明权利要求所限定的精神和范围内可对其进行许多改变，修改，甚至等效，但都将落入本发明的保护范围内。 

Claims (10)

1.一种基于GIS的中尺度流域土壤侵蚀时空动态分析方法，其特征在于，其具体过程为： 

步骤1，把流域划分为若干栅格单元，作为计算单元； 

步骤2，利用土壤侵蚀试验观测研究成果和认识，基于土壤侵蚀试验观测数据，提取关于土壤侵蚀的相关各因子； 

步骤3，将上述各因子进行分析总结，建立侵蚀因子图集，完成对每个栅格单元土壤侵蚀量的估算； 

步骤4，基于土壤侵蚀空间变化分析，完成对流域土壤侵蚀状况的评价和分析； 

所述步骤4中采用的评价模型的基本形式如下述公式(1-1)、(1-2)、(1-3)、(1-4)、(1-5)所示，其为： 

A_c=f(R，K，LS,BET,g，w)                         (1-1) 

w=f(Vegcov，wind)                         (1-2) 

式中：A_C为每个计算单元的年土壤侵蚀模数(t／km2·a)；R、K、LS、BET、g分别为降雨侵蚀因子、土壤可蚀性因子、坡度坡长因子、水土保持措施因子和沟蚀因子，w表示每个计算单元的风蚀强度等级；Vegcov、wind分别为植被盖度与风蚀范围；面蚀的计算将以CSLE为基础；At为计算流域的年平均土壤侵蚀模数(t／km2·a)，s_c为计算单元面积，S_t为流域总面积(km2)，Et为流域侵蚀产沙总量(104t)；S_o为水文观测的流域年输沙总量，Er为模型的计算相对误差(％)。 

2.根据权利要求1所述的基于GIS的中尺度流域土壤侵蚀时空动态分析方法，其特征在于，在上述步骤2中，提取关于土壤侵蚀的相关各因子的过程为： 

步骤21，计算降雨侵蚀力因子R； 

步骤22，计算土壤可蚀性因子K； 

步骤23，计算地形因素，包括坡度因子S、坡长因子L、坡度坡长因子LS； 

步骤24，计算水土保持措施因子，包括：植物措施因子B；工程措施因子E；耕作措施因子T； 

步骤25，计算流域沟蚀因子g； 

步骤26，计算流域风蚀量。 

3.根据权利要求2所述的基于GIS的中尺度流域土壤侵蚀时空动态分析方法，其特征在于，在上述步骤21中计算降雨侵蚀力因子R，利用区域气象观测资料，按照下述公式(2-1)、(2-2)和(2-3)计算， 

式中，P为年平均降水量(mm)，Pi为第i月的平均降雨量(mm)；R为多年平均降雨侵蚀力(MJ·mm·hm^-2·h^-1·a^-1)；α₂、β₂为模型参数，α₂取值为0.3589，β₂取值为1.9462，F指数mm大小与年平均雨量P的季节分布有关，取值范围在P·12-1～P之间。 

4.根据权利要求2或3所述的基于GIS的中尺度流域土壤侵蚀时空动态分析方法，其特征在于，在上述步骤22，计算土壤可蚀性因子K时，根据得到的土壤类型数据图为数据依据，按照下述公式(2-4)进行计算， 

式中，SAN、SIL和CLA是砂粒、粘粒和有机碳含量(％)，SN1=1-SAN／100；K值单位均为美制单位，t·acre·h／(100acre·ft·tonf·in)即吨·英亩·小时／(百英亩·英尺·吨力·英寸)。 

5.根据权利要求4所述的基于GIS的中尺度流域土壤侵蚀时空动态分析方法，其特征在于，坡度因子算法如公式(2-5)， 

S=10.8sinθ+0.03 

S=16.8sinθ-0.05θ<5° 

S=21.91sinθ-0.965°≤θ<14°(2-5) 

坡长因子算法如下，θ≥14° 

m=0.2θ≤1° 

L=(λ／22.1)^mm=0.31°<θ≤3°(2-6) 

m=0.43°<θ≤5° 

m=0.5θ>5° 

式中，S为坡度因子；θ为坡度；L为坡长因子；入为由DEM提取的波长；m为坡长指数，22.1为22.1m标准小区坡长；α为坡度坡长指数，黄土高原取值0.5。 

6.根据权利要求5所述的基于GIS的中尺度流域土壤侵蚀时空动态分析方法，其特征在于，植物措施因子B根据土地利用图和植物盖度图，利用N DVI根据下述公式(2-7)计算， 

式中，F为植被盖度，NDVI为归一化植被指数，NDVImin和NDVImax分别为无植被地区的NDVI值和植被良好覆盖地区的NDVI值；在现场测试基础上，结合遥感图像典型地类(沙地、高覆盖林地等)的采样，确定最高和最低植被盖度的NDVI阈值；将采样的结果与对应的NDVI值进行统计比较，当NDVI值≥NDVImax时，盖度可视为100％，当NDVI值≤NDVImin时，盖度可视为0％。 

7.根据权利要求6所述的基于GIS的中尺度流域土壤侵蚀时空动态分析方法，其特征在于，工程措施因子E根据下述公式(2-9)计算， 

式中，S_t为梯田面积，S_d为淤地坝控制面积，S为土地总面积，α，β分别为梯田和淤地坝的减沙系数，α，β分别确定为为0.764和1。 

8.根据权利要求2所述的基于GIS的中尺度流域土壤侵蚀时空动态分析方法，其特征在于， 

计算流域沟蚀因子g，在无植被覆盖的黄土坡条件下，沟蚀因子g主要受降雨因素和地形即坡度的影响，参见公式(2-10)所示， 

G=1+1.60sin(θ-15)                         (2-10)。 

9.根据权利要求2所述的基于GIS的中尺度流域土壤侵蚀时空动态 分析方法，其特征在于，根据上述各因子的计算结果，通过图形计算代数方法将上述各因子图层相乘，完成CSLE模型，得到土壤侵蚀模数图， 

在获取土壤侵蚀模数后，对计算结果进行精度验证，采用确定性系数法进行精度验证，所述确定性系数法，用来描述计算值对目标值拟合精度的无量纲统计参数，一般取值在0-1之间确定性系数E的计算公式2-11)所示， 

式中，x_i和y_i分别为实测值和模型计算值，n为年份，为实测侵蚀量年平均值。 

10.根据权利要求2或3所述的基于GIS的中尺度流域土壤侵蚀时空动态分析方法，其特征在于，在上述步骤23中，计算地形因素的过程为： 

步骤31，首先建立孤山川流域的DEM，采用原始数据是1：5万地形图，生成10m分辨率的H_C-DEM； 

步骤32，提取坡度和坡长，本发明对坡度的提取采用D8算法，若得到值为0，则将栅格值赋值为0.1°，以确保栅格之间径流连续，生成坡度图； 

坡长计算首先通过DEM进行填挖运算得到填挖后的DEM图，再计算其流向，通过运行流向图来计算流域的栅格坡长，同时结合坡度图确定坡面的径流原点及径流终点； 

在具体计算时，径流终点计算的标志用2.75°的坡度来实现，再根据计算的栅格坡长结合径流原点及径流终点来设置坡长，再次结合坡度图计算坡长指数，由坡长指数与流向图进行迭代计算得到流域的累计坡长，并最终获得流域坡长图； 

步骤33，根据上述式(2-5)、(2-6)计算坡度因子S，坡长因子L和坡度坡长因子LS。