一种基于DEM的动力学流域产沙量模拟方法与流程

一种基于dem的动力学流域产沙量模拟方法
技术领域
1.本发明涉及水土流失分析技术领域，特别涉及一种基于dem的动力学流域产沙量模拟方法。


背景技术：

2.近些年来，世界各国在集中人力和物力开发侵蚀预报经验模型的同时，注重侵蚀产沙物理过程及其物理模型的研究，并取得了大批创新性的研究成果，先后开发了通用土壤流失方程(usle)、修正的通用土壤流失方程(rusle)、产流产沙预报模型(wepp)、荷兰模型(lisem)、欧洲产流产沙预报模型(erosem)等。但由于国外物理过程模型(wepp，lisem，erosem)大都在平原地区开发和应用，缺乏我国观测数据的验证，对侵蚀物理过程描述相对简单，所以不能应用于我国黄土丘壑区陡峭多山的复杂地形领域，进行侵蚀产沙量的模拟计算。
3.国内流域侵蚀产沙模型的研究同欧美发达国家相比，起步较晚。20世纪80年代中期，以侵蚀产沙过程为基础的预报模型在国内得到了很快发展。现有技术针对黄土高原丘陵区小流域侵蚀地貌的特点，对梁坡和沟坡两大单位进行概化，利用河流推移质和悬移质输沙公式，建立了估算小流域产沙量的数学模型。现有技术还根据流域径流形成和侵蚀产沙机理，利用水文学和泥沙运动力学的基本理论，构建了流域产沙动力学模型。近几年以来，开发的流域侵蚀产沙预报模型在反映流域侵蚀产沙发生的空间分布方面取得了新的进展，能够利用黄土高原丘陵沟壑区安塞站径流小区的观测资料，建立由沟间地单元地块子模型和沟谷地单元地块子模型两部分组成的次降雨小流域地块侵蚀预报模型。
4.20世纪90年代以后，随着3s技术发展，分布式水沙模型研制得到了重视，研究人员以黑草河小流域为对象，建立了基于场次暴雨的小流域侵蚀产沙分布式模型，模型可计算出任意网格单元的产汇流和侵蚀产沙的时空分布过程，但模型没有针对丘壑区的地形进行研制，因此不能反映丘壑区中峁坡、沟坡和沟槽各自特有侵蚀特征，故此不能应用到黄土沟壑地区土壤侵蚀产沙预测评价中。
5.跟随社会发展，加大黄土丘陵沟壑区治理力度迫在眉睫，但关于黄土丘陵沟壑区的流域治理规划、设计及管理中仍然缺乏一种有效的、适用于黄土丘陵沟壑区的水土流失评价模型。因此，本技术提供了种一种基于dem的动力学流域产沙量模拟方法，采用坡面流运动波方程建立基于网格的分布式径流运动数值模型，计算得到黄土沟壑地区流域任意时刻每个网格径流深、流速、流量，根据力学概念推导出产流产沙量。


技术实现要素：

6.本发明的目的在于提供一种基于dem的动力学流域产沙量模拟方法，采用坡面流运动波方程建立基于网格的分布式径流运动数值模型，计算得到黄土沟壑地区流域任意时刻每个网格径流深、流速、流量，根据力学概念推导出产流产沙量。
7.本发明提供了一种基于dem的动力学流域产沙量模拟方法，包括以下步骤：
8.获取黄土沟壑地区的地形数据和降雨数据；
9.基于黄土沟壑地区的地形数据，获取黄土沟壑地区地表的倾斜角度，按照地表的倾斜角度对黄土沟壑地区进行区域划分；
10.根据黄土沟壑地区的地形数据，分别获取区域划分后各区域的透水面积数据、不透水面积数据以及坡面面积数据；
11.以dem网格分割黄土沟壑地区，分别建立dem网格单元的产汇流模型；
12.分别将各区域的降雨数据、透水面积数据、不透水面积数据以及坡面面积数据，输入对应dem网格单元的产汇流模型中，获得区域划分后各区域透水面积区域的产流、不透水面积区域的产流以及坡面区域的汇流；
13.建立各区域的水力侵蚀模型；
14.分别将各区域透水面积区域的产流、不透水面积区域的产流、坡面区域的汇流和降雨数据输入对应区域的水力侵蚀模型，得到对应区域的产沙量；
15.叠加各区域的产沙量，获得黄土沟壑地区的产沙量。
16.进一步地，所述基于黄土沟壑地区的地形数据，获取黄土沟壑地区地表的倾斜角度，按照地表的倾斜角度对黄土沟壑地区进行区域划分，将黄土沟壑地区划分为梁峁坡区域、沟谷坡区域和沟槽区域。
17.进一步地，所述分别将各区域的降雨数据、透水面积数据、不透水面积数据以及坡面面积数据，输入对应dem网格单元的产汇流模型中，获得区域划分后各区域透水面积区域的产流的步骤，包括：
18.获取各区域透水面积区域的降雨数据，提取透水面积区域的降雨强度数据；
19.根据黄土沟壑地区的降雨数据，采用horton下渗方程获得下渗强度数据；
20.进行降雨强度数据和下渗强度数据的数值大小对比，若降雨强度数据大于下渗强度数据，则输出透水面积区域的产流为降雨强度数据与下渗强度数据的差值；若降雨强度数据不大于下渗强度数据，则输出透水面积区域的产流为0。
21.进一步地，所述分别将各区域的降雨数据、透水面积数据、不透水面积数据以及坡面面积数据，输入对应dem网格单元的产汇流模型中，获得区域划分后各区域不透水面积区域的产流的方法包括：
22.r(t)＝p(t)
‑
e(t)
23.其中，r(t)表示降雨径流深，p(t)表示每个dem栅格上的降雨过程,e(t)表示对应的蒸发过程。
24.进一步地，所述分别将各区域的降雨数据、透水面积数据、不透水面积数据以及坡面面积数据，输入对应dem网格单元的产汇流模型中，获得区域划分后各区域坡面区域的汇流的步骤，包括：
25.获取各区域的坡面区域的降雨数据；
26.根据坡面区域的降雨数据，建立圣维南方程组描述坡面区域的水流运动；
27.求解圣维南方程组，得到坡面区域的汇流结果。
28.进一步地，所述求解圣维南方程组，得到坡面区域的汇流结果的步骤，包括：
29.根据达西定律和曼宁公式，将圣维南方程组转化为一阶拟线性坡面流偏微分方程；
30.采用preissmann四点隐式差分求解一阶拟线性坡面流偏微分方程；
31.采用牛顿迭代法获得一阶拟线性坡面流偏微分方程中的未知量，得到圣维南方程组的结果；
32.输出圣维南方程组的结果作为坡面区域的汇流结果。
33.进一步地，所述基于黄土沟壑地区的地形数据，获取黄土沟壑地区地表的倾斜角度，按照地表的倾斜角度对黄土沟壑地区进行区域划分，将黄土沟壑地区划分为梁峁坡区域、沟谷坡区域和沟槽区域。
34.进一步地，所述分别将各区域透水面积区域的产流、不透水面积区域的产流、坡面区域的汇流和降雨数据输入对应区域的水力侵蚀模型，得到对应区域的产沙量，其中获取梁峁坡区域产沙量的步骤，包括：
35.建立梁峁坡区域的水力侵蚀模型，
[0036][0037]
其中，q
c
、q
h
分别为清水流量和浑水流量；j为沟道以千分之表示的比降；j0为坡面比降；
[0038]
根据梁峁坡区域的水力侵蚀模型，得出梁峁坡区域的含沙量：
[0039][0040]
其中，s为含沙量；r为清水容重。
[0041]
进一步地，所述分别将各区域透水面积区域的产流、不透水面积区域的产流、坡面区域的汇流和降雨数据输入对应区域的水力侵蚀模型，得到对应区域的产沙量，其中获取沟谷坡区域产沙量的步骤，包括：
[0042]
建立沟谷坡区域的水力侵蚀模型，包括：
[0043]
设土壤侵蚀的能量系数为e2，则
[0044][0045]
其中，e
g
为沟谷坡单宽土壤侵蚀率；a
r
＝a'/tgα1，由实测资料适线率定；r
s
、r
m
分别为泥沙密实干容重和浑水容重；τ0‑
τ
′
c
是沟谷坡区域的有效应力；v为沟谷坡平均流速，得到沟谷坡区域的土壤侵蚀率；
[0046]
令a
g
＝e2a
r
，则：
[0047][0048]
得到沟谷坡区域的有效应力：
[0049]
τ0‑
τ
′
c
＝γ
m
h2j2 (γ
s
‑
γ
m
)d sinα2‑
f(γ
s
‑
γ
m
)d cosα2[0050]
式中，h2为沟谷坡水深；j2为沟谷坡比降；α2为沟谷坡坡度；
[0051]
沟谷坡单宽土壤侵蚀率e
g
乘以对应的产流量和过水宽度，得到沟谷坡区域的产沙量。
[0052]
进一步地，所述分别将各区域透水面积区域的产流、不透水面积区域的产流、坡面
区域的汇流和降雨数据输入对应区域的水力侵蚀模型，得到对应区域的产沙量，其中获取沟谷坡区域产沙量的步骤，包括：
[0053]
建立沟槽区域的水力侵蚀模型：
[0054][0055]
式中，e
c
为侵蚀率，u
*
为摩阻流速，h3为沟槽水深；j3为沟槽比降；κ为卡门常数；e3沟槽发生冲刷的能量系数；c为系数，由实测资料适线率定，得到沟槽区域的土壤侵蚀率e
c
；
[0056]
沟槽区域的土壤侵蚀率e
c
乘以对应的产流量和过水宽度，得到沟槽区域的产沙量。
[0057]
与现有技术相比，本发明具有如下显著优点：
[0058]
本发明提出的一种基于dem的动力学流域产沙量模拟方法，通过获取黄土沟壑地区的地形数据和降雨数据，将黄土沟壑地区按照地形划分为梁峁坡区域、沟谷坡区域和沟槽区域三种垂直分布单元，采用坡面流运动波方程建立基于网格的分布式径流运动数值模型，计算得到黄土沟壑地区流域任意时刻每个网格径流深、流速、流量，再根据力学概念推导得出产流产沙量，完成动力学流域产沙量的模拟工作。
附图说明
[0059]
图1为本发明实施例提供的黄土沟壑地区侵蚀垂直分带性示意图；
[0060]
图2为本发明实施例提供的岔巴沟流域曹坪站计算与实测四场径流过程验证结果对比图；
[0061]
图3为本发明实施例提供的岔巴沟流域曹坪站计算与实测四场输沙过程验证结果对比图。
具体实施方式
[0062]
下面结合本发明中的附图，对本发明实施例的技术方案进行清楚、完整的描述，显然，所描述的实施例是本发明的一部分实施例，而不是全部实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都应属于本发明保护的范围。
[0063]
黄土地区属干旱少雨的大陆性气候，多年平均降水量在500mm左右，中小天然流域由于植被稀少，土层厚，包气带缺水量大，表土下渗能力小。暴雨历时短，强度大，形成的洪水径流中，很少有地下径流，产流方式以超渗产流为主。
[0064]
参照图1
‑
3，本发明提供了一种基于dem的动力学流域产沙量模拟方法，包括以下步骤：
[0065]
获取黄土沟壑地区的地形数据和降雨数据；
[0066]
基于黄土沟壑地区的地形数据，获取黄土沟壑地区地表的倾斜角度，按照地表的倾斜角度对黄土沟壑地区进行区域划分，将黄土沟壑地区划分为梁峁坡区域、沟谷坡区域和沟槽区域。
[0067]
划分为梁峁坡区域、沟谷坡区域和沟槽区域三种垂直分布单元的划分规则为：
[0068]
一、梁峁坡区域顶端坡度多在5
°
以下，坡长10
‑
20m，侵蚀以溅蚀和面蚀为主。梁峁坡上部，坡度在20
°
以下，坡长20
‑
30m，以细沟侵蚀为主，间或有浅沟侵蚀发生。梁峁坡中下部地形比较复杂，坡度在20
°‑
30
°
之间，坡长15
‑
20m，细沟侵蚀进一步发育，以浅沟侵蚀为主，间或有坡面切沟和陷穴侵蚀发生。
[0069]
二、沟谷坡区域指峁边线以下至坡脚线以上的地带。是以切沟侵蚀、重力侵蚀及洞穴侵蚀为主的剧烈侵蚀地带。谷坡侵蚀受梁峁坡来水影响很大，若控制坡面来水，则谷坡侵蚀量将会大大减少。沟谷坡坡度一般在30
°‑
45
°
之间。
[0070]
三、沟槽区域指坡脚线以下的地段，包括沟条地和沟床地两部分。沟床包括石质沟床与土质沟床，后者以沟岸扩张、沟底下切和溯源侵蚀为主。
[0071]
根据黄土沟壑地区的降雨数据，建立基于网格单元的产汇流模型，分别计算梁峁坡区域、沟谷坡区域和沟槽区域的透水面积区域的产流、不透水面积区域的产流和坡面区域的汇流；
[0072]
黄土沟壑区坡面具有明显的垂直分带特征。通常，可以根据其地形、地貌特征和侵蚀形式分为梁峁坡、沟谷坡和沟槽(或沟谷底)三种垂直分带单元，如图1所示，其中：ⅰ代表梁峁坡；ⅱ代表沟谷坡；ⅲ代表沟槽或沟谷底。建立网格单元的产汇流模型。
[0073]
计算梁峁坡区域、沟谷坡区域和沟槽区域的透水面积区域的产流的步骤包括：
[0074]
获取梁峁坡区域、沟谷坡区域和沟槽区域透水面积区域的降雨数据，提取透水面积区域的降雨强度数据i；
[0075]
根据梁峁坡区域、沟谷坡区域和沟槽区域透水面积区域的降雨数据，采用horton下渗方程获得下渗强度数据f，f＝f
c
(f0‑
f
c
)e
‑
kt
，其中，i是降雨强度数据；r是产流量；f是下渗能力；f0是最大下渗能力；f
c
是稳定下渗能力；k是随土质而变的指数；t是时间；
[0076]
进行降雨强度数据和下渗强度数据的数值大小对比，若降雨强度数据大于下渗强度数据，则输出透水面积区域的产流为降雨强度数据与下渗强度数据的差值；若降雨强度数据不大于下渗强度数据，则输出透水面积区域的产流为0，数值大小对比公式为：
[0077]
计算梁峁坡区域、沟谷坡区域和沟槽区域的不透水面积区域的产流的方法为：
[0078]
r(t)＝p(t)
‑
e(t)
[0079]
其中，r(t)表示降雨径流深，p(t)表示每个dem栅格上的降雨过程,e(t)表示对应的蒸发过程。
[0080]
计算梁峁坡区域、沟谷坡区域和沟槽区域的坡面区域的汇流的步骤包括：
[0081]
获取梁峁坡区域、沟谷坡区域和沟槽区域的坡面区域的降雨数据；
[0082]
根据坡面区域的降雨数据，建立圣维南方程组描述坡面区域的水流运动，坡面水流连续方程表示为：式中，q为单宽流量；h为坡面流水深；r
e
(t)为净雨过程即产流过程；圣维南方程组的运动波近似特别适用于表面粗糙、坡度陡、旁侧来流少的水流运动，坡面流都可以用运动波方程来描述，即：
[0083][0084]
求解圣维南方程组，根据达西定律和曼宁公式将圣维南方程组转化为一阶拟线性坡面流偏微分方程
[0085]
一阶拟线性坡面流偏微分方程的初始条件和边界条件为：
[0086][0087]
采用preissmann四点隐式差分求解一阶拟线性坡面流偏微分方程：
[0088][0089]
其中，θ为权重系数，0≤θ≤1，从计算格式稳定性需要出发，θ宜大于0.5，最好取在0.6≤θ≤1，用牛顿迭代法直接解得其中唯一的未知量据此，能推求出任意时空不均匀降雨的坡面流量过程；
[0090]
采用牛顿迭代法获得一阶拟线性坡面流偏微分方程中的未知量，得到圣维南方程组的结果；
[0091]
输出圣维南方程组的结果作为坡面区域的汇流结果。
[0092]
建立梁峁坡区域的水力侵蚀模型、沟谷坡区域的水力侵蚀模型和沟槽区域的水力侵蚀模型；
[0093]
将梁峁坡区域透水面积区域的产流、不透水面积区域的产流以及坡面区域的汇流的数据和降雨数据输入梁峁坡区域的水力侵蚀模型，得到梁峁坡区域的产沙量、土壤侵蚀率和有效应力的步骤包括：
[0094]
建立梁峁坡区域的水力侵蚀模型，包括：
[0095][0096]
其中，q
c
、q
h
分别为清水流量和浑水流量；j为沟道以千分之表示的比降；j0为坡面比降；
[0097]
根据梁峁坡区域的水力侵蚀模型，得出梁峁坡区域的含沙量：
[0098]
[0099]
其中，s为含沙量；r为清水容重；
[0100]
根据梁峁坡区域的含沙量，得到梁峁坡区域的土壤侵蚀率：
[0101][0102]
其中，r
s
、r
m
分别为泥沙密实干容重和浑水容重；
[0103][0104]
式中：e
r
为梁峁坡单宽土壤侵蚀率；；τ0‑
τ
c
梁峁坡区域的有效应力；v为梁峁坡平均流速；；a
r
＝a'/tgα1，由实测资料适线率定；
[0105]
所述梁峁坡区域的地形数据中水深为h1，比降为j1，坡角为α1，得到梁峁坡区域的有效应力：
[0106]
τ0‑
τ
c
＝r
m
h1j1 (r
s
‑
r
m
)d sinα1‑
f(r
s
‑
r
m
)d cosα1[0107]
式中：d为泥沙粒径；f为摩擦系数，取0.047。
[0108]
将沟谷坡区域透水面积区域的产流、不透水面积区域的产流以及坡面区域的汇流的数据和降雨数据输入沟谷坡区域的水力侵蚀模型，得到沟谷坡区域的产沙量、土壤侵蚀率和有效应力的步骤包括：
[0109]
建立沟谷坡区域的水力侵蚀模型，包括：
[0110]
设土壤侵蚀的能量系数为e2，则
[0111][0112]
其中，e
g
为沟谷坡单宽土壤侵蚀率；a
r
＝a'/tgα1，由实测资料适线率定；r
s
、r
m
分别为泥沙密实干容重和浑水容重；τ0‑
τ
′
c
是沟谷坡区域的有效应力；v为沟谷坡平均流速，得到沟谷坡区域的土壤侵蚀率；
[0113]
令a
g
＝e2a
r
，则：
[0114][0115]
得到沟谷坡区域的有效应力：
[0116]
τ0‑
τ
′
c
＝γ
m
h2j2 (γ
s
‑
γ
m
)d sinα2‑
f(γ
s
‑
γ
m
)d cosα2[0117]
式中，h2为沟谷坡水深；j2为沟谷坡比降；α2为沟谷坡坡度；
[0118]
沟谷坡单宽土壤侵蚀率e
g
乘以对应的产流量和过水宽度，得到沟谷坡区域的产沙量。
[0119]
将沟槽区域透水面积区域的产流、不透水面积区域的产流以及坡面区域的汇流的数据和降雨数据输入沟槽区域的水力侵蚀模型，得到沟槽区域的产沙量、土壤侵蚀率和有效应力的步骤包括：
[0120]
建立沟槽区域的水力侵蚀模型，包括：
[0121][0122]
式中，e
c
为侵蚀率，u
*
为摩阻流速，h3为沟槽水深；j3为沟槽比降；κ为卡门常数；e3沟槽发生冲刷的能量系数；c为系数，由实测资料适线率定，得到沟槽区域的土壤侵蚀率e
c
；
[0123]
沟槽区域的土壤侵蚀率e
c
乘以对应的产流量和过水宽度，得到沟槽区域的产沙量。
[0124]
实施例1
[0125]
黄土沟壑地区主沟壑密度大，地形破碎，植被稀少，土质松，暴雨多，侵蚀严重。黄土沟壑地区主要地貌可概括为沟谷地和沟间地两种。由于该区坡面切割强烈，地形破碎、沟壑纵横，加之长期以来不合理的土地利用，植被稀少，暴雨集中并且强度较大，水土流失极为严重，土壤侵蚀模数大多在1万t/(km2
·
a)以上，部分地区竞高达3万t/(km2
·
a)，是黄河粗泥沙的主要来源区。本技术选择该区观测资料比较完善的典型小流域岔巴沟流域进行侵蚀产沙模型研究。黄土沟壑地区坡度很陡，使得洪水波传播速度很快，沿程坦化小，具有运动波的传播特征，所以在本技术中采用坡面流运动波方程来建立基于网格的分布式径流运动数值模型，利用该模型可以计算得到流域任意时刻每个网格径流深、流速、流量，而且重要的是，本模型是具有一定物理基础的小流域分布式径流模拟数值模型，参数明确，具有明确的物理意义。
[0126]
岔巴沟流域位于东经109
°
47'、北纬37
°
21'，自然地理区划属于黄土丘陵沟壑区第一副区，在无定河流域的西南部与无定河的一级支流大理河相汇，流域面积为187km2，沟道长24.1km。该流域地貌的基本特征是土壤侵蚀严重，沟谷发育剧烈，全流域被大小沟道切割成支离破碎、沟壑纵横的典型黄土地貌景观。地面坡度变化复杂，且不连续，同一峁、梁的各个方向，及同一方向的上、中、下坡面坡度变化都很急剧。沟谷坡面多为陡峭坡面，坡度一般大于60
°
；主沟上游及较大支流坡度在54
°‑
60
°
之间；在沟头及支沟上部减至30
°‑
45
°
之间；峁梁坡
‑
梁顶部坡度平缓，约在5
°‑
10
°
以内；梁的两侧较陡，峁腰上部较陡，下部较缓，变化范围在15
°‑
30
°
之间。
[0127]
坡面及沟道的侵蚀物质必须输移到出口断面才能成为流域的产沙量，侵蚀物质的输送受到水流条件、泥沙条件和沟道边界条件等多种因素的制约。研究认为，在黄土高原绝大部分地区的沟道或河道内都是基岩出露，或者虽然不是基岩河床，但沉积物都较粗属砾石或粗沙，其泥沙输移比接近于1.0。泥沙输移比的大小还与流域面积有关，流域越小，泥沙输移比越大。本技术研究的流域属小流域，因此，假定泥沙输移比为1.0，能够满足计算精度要求，并且可使输移计算大大简化。所以，流域出口断面的产沙量可由各单元侵蚀产沙量错开若干个传播时段迭加求得，可参考汇流时间进而确定传播时段。
[0128]
模型参数率定依据岔巴沟流域实测水文泥沙资料，从1970
‑
2000年的十七场次洪资料洪水中选取了前9场水沙过程进行了模型参数率定，模型率定结果如表1所示。
[0129]
表1水沙侵蚀模型参数率定结果
[0130][0131]
取各组参数的平均值作为最后参数对1970
‑
2000年后8场水沙过程进行验证，水文模型和侵蚀产沙模型验证结果特征值统计见下表。
[0132]
表2水沙侵蚀模型参数验证结果
[0133][0134]
本次模型验证的八场产沙过程中，峰现时间把握的比较好，只有990720次滞后1个小时、000704和010818次均滞后0.5个小时，其他场次均无时差。960731次模拟产沙过程中模拟产沙量误差最大(100％)，而沙峰误差较小(40.5％)；950902次模拟产沙过程中模拟产沙量误差次之(63.4％)，沙峰误差为最大(162％)；940804为验证资料中实测总产沙量最大的场次，模拟产沙量误差为41.6％，沙峰误差为36.7％，效果尚可。010818次模拟产沙过程中模拟产沙量误差次之(33.3％)，沙峰误差较大(67.8％)；其他场次产沙量和沙峰误差均较小，最大分别为25.5％和54.9％。侵蚀产沙模型和产汇流模型相比，计算精度较低。但与前人的成果相比，由于采用基于栅格的分布式模型和力学概念推导出产流产沙公式，对不同的下垫面分别采取了不同的参数值，公式参数具有清晰明确的物理概念，模拟效果比整个流域采取单一的参数值更加合理和接近实际。
[0135]
以上公开的仅为本发明的几个具体实施例，但是，本发明实施例并非局限于此，任何本领域的技术人员能思之的变化都应落入本发明的保护范围。