一种流域侵蚀产沙指数模拟方法和系统

1.本发明涉及模拟流域侵蚀产沙量技术领域，更具体地，涉及一种流域侵蚀产沙指数模拟方法和系统。


背景技术：

2.流域侵蚀产沙采用单位面积产沙量syi进行表征，是指在特定时段内，流域产沙总量与流域面积的比值，目前流域侵蚀产沙量模拟的方法主要有两类，分别是：土壤侵蚀模型(rusle)和水蚀预测模型(wepp)。
3.rusle土壤侵蚀模型的表达式为a＝rklspc，式中，a为单位面积年平均土壤侵蚀量，t/hm2；r为降雨侵蚀力因子，mj
·
mm/(hm2·h·
a)；k为土壤可蚀性因子，t
·
h/(mj
·
mm)；l为坡长因子；s为坡度因子；p为水土保持措施因子；c为植被覆盖因子。相关研究表明，rusle模型已取得了很好应用。该模型考虑的影响因子全面，但是也存在一定的局限性。应用中需要考虑标准小区的概念(长22.13m、宽1.85m、坡度为9％的地块)，数据获取比较繁琐，且数据的计算处理量大；同时，它只能适用于坡度不太大的坡耕地的预测，且下垫面仅考虑水土保持和植被覆盖等因素，未能充分考虑土地利用变化及前期降雨等因素对流域侵蚀产沙的影响。
4.水蚀预测模型(wepp)是美国农业部于1985年主持开发的一种基于侵蚀过程的模型，较rusle模型更为先进，具有侵蚀模块、水文过程模块、植物生长模块、天气随机生成模块、残留物分解模块、土壤模块、冬季过程模块、灌溉模块、地表径流模块等9个功能模块。wepp模型作为土壤侵蚀预测的技术工具，对土壤侵蚀科学研究提供了许多借鉴，已在国内外得到了快速的发展和应用。但模型中包含的子模型和参数较多，需要大量资料对模型参数进行修正和验证。
5.现有技术的耦合不同时空尺度模型的流域侵蚀产沙量预测方法，根据流域的降雨侵蚀力因子rl、土壤可蚀性因子kl、坡长因子lsl、植被因子cl、工程措施因子el、耕作措施因子tl和沟蚀系数因子gl，计算wl的年侵蚀产沙量，系统参数多，数据获取繁琐，数据计算处理量大。


技术实现要素：

6.本发明为克服上述技术问题，提供数据处理量小，模型参数较少，考虑土地利用变化的一种流域侵蚀产沙指数模拟方法和系统。
7.本发明技术方案如下：
8.一种流域侵蚀产沙指数模拟方法，用于计算流域一段时间内的侵蚀产沙指数，包括步骤：
9.s1、获取流域控制面上的第一cn值和对应的单位面积产沙量syi；
10.s2、将cn值作为自变量，单位面积产沙量syi作为因变量，建立线性回归模型；
11.s3、利用第一cn值和对应的单位面积产沙量syi对线性回归模型参数进行率定和
验证，得到syi计算模型；
12.s4、获取第二cn值，所述第二cn值用于计算待测的单位面积产沙量syi，将第二cn值输入syi计算模型，计算得到单位面积产沙量syi。
13.本技术方案提出了一种流域侵蚀产沙指数模拟方法，将径流曲线数cn值作为自变量，单位面积产沙量syi作为因变量，建立线性回归模型，系统参数少，简化了计算过程，cn值通过查表得到，数据获取简便，无需处理大量数据，并且考虑土地利用变化以及前期土壤湿润度等级数据用于侵蚀产沙量的模拟和预测，提高了预测精度。
14.进一步地，得到所述土地利用类型数据的方法为，获取流域地形图和tm/etm影像图，通过流域地形图和tm/etm影像图得到土地利用/覆盖图，利用土地利用/覆盖图得到土地利用类型数据；
15.前期土壤湿润度等级数据的计算方法为，先获取流域控制面上各个雨量站前五天降雨总量，查询前期土壤湿润度分类表得到前期土壤湿润度等级数据。
16.上述技术方案中，数据获取方便，目前可以直接从资源环境数据云平台(http://www.resdc.cn)下载所需年份的土地利用遥感监测数据作为土地类型数据。
17.进一步地，步骤s2所述线性回归模型公式为：
18.syi＝a0*cn β019.其中，a0为模型的回归系数，β0为模型的常量。
20.进一步地，步骤s3利用最小二乘法对线性回归模型参数进行率定，确定模型的回归系数a0和常量β0的取值，得到syi计算模型。
21.进一步地，步骤s3进行率定时，将第一cn值和单位面积产沙量syi按照时间尺度分为逐日、逐半月、逐年三组数据，按照每一组时间尺度的第一cn值和单位面积产沙量syi，分别得到该时间尺度对应的回归系数a0和常量β0的取值；
22.其中，各个时间尺度的第一cn值获取方法是：
23.获取流域控制面上各个雨量站的日cn值，运用arcgis空间插值技术得到流域控制面上的逐日时间尺度的第一cn值；
24.获取流域控制面上各个雨量站的半月cn均值，运用arcgis空间插值技术得到流域控制面上的逐半月时间尺度的第一cn值；
25.获取流域控制面上各个雨量站的年cn均值，运用arcgis空间插值技术得到流域控制面上的逐年时间尺度的第一cn值；
26.步骤s4所述第二cn值为待测的单位面积产沙量syi对应的cn值，根据待测的单位面积产沙量syi的时段获取第二cn值，若待测的单位面积产沙量syi为逐日时段，则从雨量站获取逐日时间尺度的cn值作为第二cn值；若待测的单位面积产沙量syi为半月时段，则运用arcgis空间插值技术得到逐半月时间尺度的第二cn值；若待测的单位面积产沙量syi为年时段，则运用arcgis空间插值技术得到逐年时间尺度的第二cn值。
27.进一步地，步骤s2所述线性回归模型还包括参数：降雨侵蚀力r，线性回归模型公式为：
28.syi＝a1*cn b1*r β129.其中，cn值和降雨侵蚀力r作为自变量，0≤cn≤100，单位面积产沙量syi作为因变量，a1和b1为模型的回归系数，β1为模型的常量，r为降雨侵蚀力，通过半月时段降雨侵蚀力
模型利用日雨量资料计算降雨侵蚀力r。
30.进一步地，步骤s3利用最小二乘法对线性回归模型参数进行率定，确定模型的回归系数a1、b1，模型的常量β1，得到syi计算模型。
31.进一步地，步骤s3进行率定时，将降雨侵蚀力r、第一cn值和单位面积产沙量syi按照时间尺度分为逐半月、逐年两组数据，按照每一组时间尺度的降雨侵蚀力r、第一cn值和单位面积产沙量syi，分别得到该时间尺度对应的回归系数a1和b1，以及常量β1的取值；
32.其中，各个时间尺度的第一cn值获取方法是：
33.获取流域控制面上各个雨量站的半月cn均值，运用arcgis空间插值技术得到流域控制面上的逐半月时间尺度的第一cn值；
34.获取流域控制面上各个雨量站的年cn均值，运用arcgis空间插值技术得到流域控制面上的逐年时间尺度的第一cn值；
35.各个时间尺度的降雨侵蚀力r获取方法是：
36.半月时段降雨侵蚀力模型将每个月前15天作为第1个时段，剩下天数作为第2个时段，因此每年划分为24个时段；具体公式如下：
[0037][0038]
式中：ri表示第i个半月的侵蚀力值，mj
·
mm/(hm2·
h)，i＝1，2，3，
…
，24；k表示该半月时段内的降雨天数；pj表示半月时段内日雨量≥12mm的第j天日雨量，《12mm的日雨量不纳入计算；反映当地降雨特性的模型参数α和β由以下公式确定：
[0039][0040]
α＝21.586β-7.1891
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)
[0041]
式中：p
d12
表示日雨量≥12mm的日平均雨量，p
y12
表示日雨量≥12mm的年平均雨量；利用式(1)～式(3)计算流域控制面上各个雨量站逐半月的降雨侵蚀力值，再累加1年的24个半月时段得到各个雨量站的年降雨侵蚀力值。
[0042]
一种流域侵蚀产沙指数模拟系统，包括：数据获取单元、模型构建单元、模型率定单元、syi计算单元；
[0043]
数据获取单元获取流域控制面上的第一cn值和对应的单位面积产沙量syi；模型构建单元将cn值作为自变量，单位面积产沙量syi作为因变量，建立线性回归模型；模型率定单元利用第一cn值和对应的单位面积产沙量syi对线性回归模型参数进行率定和验证，得到syi计算模型；syi计算单元获取第二cn值，所述第二cn值用于计算待测的单位面积产沙量syi，将第二cn值输入syi计算模型，计算得到单位面积产沙量syi；
[0044]
其中，所述第一cn值是通过土地利用类型数据、土壤类型数据及前期土壤湿润度等级数据，查询cn值表得到。
[0045]
进一步地，得到所述土地利用类型数据的方法为，获取流域地形图和tm/etm影像图，通过流域地形图和tm/etm影像图得到土地利用/覆盖图，利用土地利用/覆盖图得到土地利用类型数据；
[0046]
前期土壤湿润度等级数据的计算方法为，先获取流域控制面上各个雨量站前五天降雨总量，查询前期土壤湿润度分类表得到前期土壤湿润度等级数据。
[0047]
本发明提出了一种流域侵蚀产沙指数模拟方法和系统，与现有技术相比，本发明技术方案的有益效果是：本技术方案将径流曲线数cn值和降雨侵蚀力r作为自变量，单位面积产沙量syi作为因变量，建立线性回归模型，系统参数少，简化了计算过程，cn值通过查表得到，数据获取简便，无需处理大量数据，并且考虑土地利用变化以及前期土壤湿润度等级数据用于侵蚀产沙量的模拟和预测，提高了预测精度。
附图说明
[0048]
图1为实施例1流域侵蚀产沙指数模拟方法步骤示意图；
[0049]
图2为本发明流域侵蚀产沙指数模拟方法流程图；
[0050]
图3为龙川子流域逐日cn均值箱线图；
[0051]
图4为博罗子流域逐日cn均值箱线图；
[0052]
图5为实施例3流域侵蚀产沙指数模拟方法步骤示意图；
[0053]
图6为年平均降雨侵蚀力空间分布图；
[0054]
图7为土地利用/覆盖示意图；
[0055]
图8为土壤类型示意图。
具体实施方式
[0056]
附图仅用于示例性说明，不能理解为对本专利的限制；为了更好说明本实施例，附图某些部件会有省略、放大或缩小，并不代表实际产品的尺寸；对于本领域技术人员来说，附图中某些公知结构及其说明可能省略是可以理解的。下面结合附图和实施例对本发明的技术方案做进一步的说明。
[0057]
实施例1
[0058]
一种流域侵蚀产沙指数模拟方法如图1所示，包括步骤：
[0059]
s1、获取流域控制面上的第一cn值和对应的单位面积产沙量syi；
[0060]
s2、将cn值作为自变量，单位面积产沙量syi作为因变量，建立线性回归模型；
[0061]
s3、利用第一cn值和对应的单位面积产沙量syi对线性回归模型参数进行率定和验证，得到syi计算模型；
[0062]
s4、获取待测的单位面积产沙量syi对应第二cn值，将第二cn值输入syi计算模型，计算得到单位面积产沙量syi。
[0063]
本实施例将径流曲线数cn值作为自变量，单位面积产沙量syi作为因变量，建立线性回归模型，系统参数少，简化了计算过程，cn值通过查表得到，数据获取简便，无需处理大量数据，并且考虑土地利用变化以及前期土壤湿润度等级数据用于侵蚀产沙量的模拟和预测，提高了预测精度。
[0064]
实施例2
[0065]
本实施例提供一种流域侵蚀产沙指数模拟方法，流程图如图2所示，在实施例1的基础上，得到所述土地利用类型数据的方法为，获取流域地形图和tm/etm影像图，通过流域地形图和tm/etm影像图结合制作的方法得到土地利用/覆盖图，利用土地利用/覆盖图得到土地利用类型数据；
[0066]
前期土壤湿润度等级数据的计算方法为，先获取流域控制面上各个雨量站前五天
降雨总量，查询前期土壤湿润度分类表得到前期土壤湿润度等级数据。
[0067]
步骤s2所述线性回归模型公式为：
[0068]
syi＝a0*cn β0[0069]
其中，a0为模型的回归系数，β0为模型的常量。
[0070]
步骤s3利用最小二乘法对线性回归模型参数进行率定，确定模型的回归系数a0和常量β0的取值，得到syi计算模型。
[0071]
步骤s3进行率定时，将第一cn值和单位面积产沙量syi按照时间尺度分为逐日、逐半月、逐年三组数据，按照每一组时间尺度的第一cn值和单位面积产沙量syi，分别得到该时间尺度对应的回归系数a0和常量β0的取值；
[0072]
其中，各个时间尺度的第一cn值获取方法是：
[0073]
获取流域控制面上各个雨量站的日cn值，运用arcgis空间插值技术得到流域控制面上的逐日时间尺度的第一cn值；
[0074]
获取流域控制面上各个雨量站的半月cn均值，运用arcgis空间插值技术得到流域控制面上的逐半月时间尺度的第一cn值；
[0075]
获取流域控制面上各个雨量站的年cn均值，运用arcgis空间插值技术得到流域控制面上的逐年时间尺度的第一cn值；
[0076]
步骤s4所述第二cn值为待测的单位面积产沙量syi对应的cn值，根据待测的单位面积产沙量syi的时段获取第二cn值，若待测的单位面积产沙量syi为逐日时段，则从雨量站获取逐日时间尺度的cn值作为第二cn值；若待测的单位面积产沙量syi为半月时段，则运用arcgis空间插值技术得到逐半月时间尺度的第二cn值；若待测的单位面积产沙量syi为年时段，则运用arcgis空间插值技术得到逐年时间尺度的第二cn值。
[0077]
利用统计分析软件spss绘出逐日cn值的年际分布箱线图，逐日cn值的年际分布箱线图如图3和图4所示。
[0078]
实施例3
[0079]
流域侵蚀产沙指数模拟方法步骤示意图如图5所示，本实施例在实施例1的基础上，步骤s2所述线性回归模型还包括参数：降雨侵蚀力r，线性回归模型公式为：
[0080]
syi＝a1*cn b1*r β1[0081]
其中，cn值和降雨侵蚀力r作为自变量，0≤cn≤100，单位面积产沙量syi作为因变量，a1和b1为模型的回归系数，β1为模型的常量，r为降雨侵蚀力，通过半月时段降雨侵蚀力模型利用日雨量资料计算降雨侵蚀力r。
[0082]
步骤s3利用最小二乘法对线性回归模型参数进行率定，确定模型的回归系数a1、b1，模型的常量β1，得到syi计算模型。
[0083]
步骤s3进行率定时，将降雨侵蚀力r、第一cn值和单位面积产沙量syi按照时间尺度分为逐半月、逐年两组数据，按照每一组时间尺度的降雨侵蚀力r、第一cn值和单位面积产沙量syi，分别得到该时间尺度对应的回归系数a1和b1，以及常量β1的取值；
[0084]
其中，各个时间尺度的第一cn值获取方法是：
[0085]
获取流域控制面上各个雨量站的半月cn均值，运用arcgis空间插值技术得到流域控制面上的逐半月时间尺度的第一cn值；
[0086]
获取流域控制面上各个雨量站的年cn均值，运用arcgis空间插值技术得到流域控
制面上的逐年时间尺度的第一cn值；
[0087]
各个时间尺度的降雨侵蚀力r获取方法是：
[0088]
半月时段降雨侵蚀力模型将每个月前15天作为第1个时段，剩下天数作为第2个时段，因此每年划分为24个时段；具体公式如下：
[0089][0090]
式中：ri表示第i个半月的侵蚀力值，mj
·
mm/(hm2·
h)，i＝1，2，3，
…
，24；k表示该半月时段内的降雨天数；pj表示半月时段内日雨量≥12mm的第j天日雨量，《12mm的日雨量不纳入计算；反映当地降雨特性的模型参数α和β由以下公式确定：
[0091][0092]
α＝21.586β-7.1891
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)
[0093]
式中：p
d12
表示日雨量≥12mm的日平均雨量，p
y12
表示日雨量≥12mm的年平均雨量；利用式(1)～式(3)计算流域控制面上各个雨量站逐半月的降雨侵蚀力值，再累加1年的24个半月时段就可以得到各个雨量站的年降雨侵蚀力值。
[0094]
本实施例采用线性回归模型，系统参数少，简化了计算过程，查表得到cn值，通过降雨侵蚀力模型得到r值，无需处理大量数据，并且考虑土地利用变化以及前期土壤湿润度等级数据用于侵蚀产沙量的模拟和预测，提高了预测精度。
[0095]
利用克里金法对降雨侵蚀力r进行插值处理，得到平均降雨侵蚀力空间分布图，如图6所示。
[0096]
实施例4
[0097]
本实施例以广东省东江流域为例，提供一种流域侵蚀产沙指数模拟方法，包括步骤：
[0098]
s1、数据获取，计算cn值、降雨侵蚀力r、单位面积产沙量syi；
[0099]
通过实地考察和相关部门调研等收集补充整理流域气象、水文、流域地形、tm/etm影像图、土地利用/覆盖图。根据流域地貌特征及水系特点，结合流域dem，运用arcgis空间分析技术划分子流域，龙川、博罗是东江流域的2个主要水文站。本实施例分别以龙川、博罗2个水文站控制的子流域为研究单元，将东江流域划分为龙川子流域、博罗子流域进行比较研究；结合遥感技术，从现有的地形图和tm/etm影像上获取土地利用/覆盖示意图如图7所示，建立项目研究数据库。
[0100]
cn值的确定：cn值依据流域土壤类型、土地利用/覆盖措施及前期土壤湿润度(amc)综合确定。根据土壤水文性质或土壤最小下渗率，可以将土壤分为a、b、c、d四类，土地类型示意图如图8所示。东江流域20世纪80年代的主要用地面积分别为耕地(18％)、林地(73％)、草地(4％)、水域(3％)。美国土壤保持局根据前5d降雨总量把前期土壤湿润度划分为amcⅰ(干旱)、amcⅱ(正常)、amcⅲ(湿润)三个等级，如表1所示。根据流域上各站点的土地利用类型、土壤类型及前期土壤湿润度等级，通过查表确定两个子流域控制面上各雨量站位置逐日时段的cn值。运用arcgis空间插值技术，插值得到两个子流域控制面上逐日时段的cn均值，用统计分析软件spss绘出这两组数据的年际分布箱线图，如图3和图4所示。逐半月和逐年时段的cn值以同样的插值方法得到。
[0101][0102]
表1前期土壤湿润度分类
[0103]
降雨侵蚀力r的计算：采用半月时段降雨侵蚀力模型，利用日雨量资料计算降雨侵蚀力。将各雨量站逐半月时段的降雨侵蚀力值进行累加得到各雨量站的年降雨侵蚀力值。由年降雨侵蚀力值计算可得到各站点多年平均降雨侵蚀力值，运用arcgis空间插值技术克里金法插值得到流域上的多年平均降雨侵蚀力空间分布图，如图6所示。
[0104]
模型输出参数：单位面积产沙量syi根据逐日平均悬移质输沙率计算得到。
[0105]
s2、构建经验产沙指数模型；
[0106]
将龙川、博罗两个子流域分别作为一个产沙系统，以cn作为自变量，syi作为因变量。假设cn与syi呈线性关系，建立二者的一元线性回归模型syi＝a0*cn β0。为了提高模型的模拟精度，考虑在模型中进一步增加降雨侵蚀力r参数作为自变量。同理，构建多元线性回归模型syi＝a1*cn b1*r β1。模型的回归系数及常量由最小二乘法确定。
[0107]
s3、模型的率定与验证。
[0108]
运用实测数据对两个模型进行不同时间尺度的率定和验证，结果见表2。
[0109][0110]
表2模型率定与验证结果
[0111]
采用绝对误差(ape)、平方误差(ise)、效率系数(ce)及相关系数(cc)模型性能指
标验证模型精度，ape值小于0.35、ise值小于0.15、ce值大于0.60、cc值大于0.75的模型被认为是拟合较好的模型。模型精度验证结果见表3。
[0112][0113]
表3模型性能评价结果
[0114]
本实施例用过去某一时段的观测资料建立模型，并用另一时段的观测资料对建立的模型的模拟效果进行验证，确定该模型能较准确的用于预测后，就可以对流域内过去某一时段中缺少观测记录的那部分数据进行预测。
[0115]
并且，本实施例中，通过回归模型结果比较得出：半月时间尺度下，同时包含cn和r两个参数的模型能较准确的用于东江流域的侵蚀产沙量预测。
[0116]
对于缺乏产沙量观测资料的时段，借助该模型进行预测时，只需确定该时段下的
cn和r两个参数的值，代入模型的回归方程中就能推求出相应缺乏实测资料时段的预测syi值。
[0117]
从龙川和博罗两个子流域不同时间尺度模型的模型性能评价结果(表3)中可以看出，只包含cn参数的模型在逐日时间尺度和逐年时间尺度上的模型模拟精度不高，而在半月时间尺度上的模型模拟精度较好，龙川子流域和博罗子流域在半月时间尺度上的效率系数分别为0.61和0.62。在模型中增加了r参数后，龙川子流域和博罗子流域在逐年时间尺度上的模拟精度都满足了模型模拟精度的要求，效率系数分别达到0.62和0.81；同时，龙川子流域和博罗子流域在半月时间尺度上的效率系数也都得到了提升，平均提高了16个百分点。
[0118]
实施例5
[0119]
一种流域侵蚀产沙指数模拟系统，包括：数据获取单元、模型构建单元、模型率定单元、syi计算单元；
[0120]
数据获取单元获取流域控制面上的第一cn值和对应的单位面积产沙量syi；模型构建单元将cn值作为自变量，单位面积产沙量syi作为因变量，建立线性回归模型；模型率定单元利用第一cn值和对应的单位面积产沙量syi对线性回归模型参数进行率定和验证，得到syi计算模型；syi计算单元获取第二cn值，所述第二cn值用于计算待测的单位面积产沙量syi，将第二cn值输入syi计算模型，计算得到单位面积产沙量syi；
[0121]
其中，所述第一cn值的时段为包括逐日时段、逐半月时段和逐年时段，逐日时段的第一cn值是通过土地利用类型数据、土壤类型数据及前期土壤湿润度等级数据，查询cn值表得到。得到所述土地利用类型数据的方法为，获取流域地形图和tm/etm影像图，通过流域地形图和tm/etm影像图得到土地利用/覆盖图，利用土地利用/覆盖图得到土地利用类型数据；在本发明其他实施例中，还可以直接从资源环境数据云平台(http://www.resdc.cn)下载所需年份的土地利用遥感监测数据作为土地类型数据。
[0122]
前期土壤湿润度等级数据的计算方法为，先获取流域控制面上各个雨量站前五天降雨总量，查询前期土壤湿润度分类表得到前期土壤湿润度等级数据。