一种基于CAESAR-Lisflood模型的流域泥沙溯源模拟方法与流程

一种基于caesar
‑
lisflood模型的流域泥沙溯源模拟方法
技术领域
1.本发明涉及地表(地貌)数值模拟技术领域，特别是涉及一种基于caesar
‑
lisflood模型的流域泥沙溯源模拟方法。


背景技术：

2.目前国内外运用各类数值模拟及物理模型对地球表层(地貌演变)系统的模拟研究较多，对流域物质运移及河道变化的模拟研究已经发展为新的研究热点。例如运用土壤侵蚀预报模型来研究水土流失，如通用土壤流失方程(usle)和weep土壤侵蚀模型。随着地理信息系统(gis)以及遥感技术的发展，更多研究尝试将数字高程模型(dem)应用在土壤侵蚀类型与过程研究中，应用gis系统建立流域降雨条件下土壤侵蚀的定量计算方法，用以模拟土壤侵蚀的空间分布格局及变化规律。
3.目前结合坡面失稳模块的地表(地貌)演变模型较少，且相关模型缺乏对流域不同区域或事件的产沙及其后续运移过程的追踪技术。例如，常见的地表(地貌)演变模型结果大都只能输出流域总的产沙量，无法单独对滑坡产沙及其后续运移过程进行追踪模拟，即无法将滑坡产沙和非滑坡产沙区分开来。另外，现有的地表(地貌)演变模型虽然能模拟不同情景下不同时空尺度流域演变情况，能对流域特定情景下总产沙量及地表高程变化进行预测，但对极端灾害情景(例如滑坡)下的物质产生及运移过程模拟还存在欠缺。
4.因此，需要新的技术方法，以至少部分解决现有技术中存在的问题。


技术实现要素：

5.凯撒二维水动力地表景观演变模型(caesar
‑
lisflood)模型是基于gis分布式的地表景观演变模型caesar以及二维水动力模型lisflood而开发的。该模型具备一定的河流二维水动力学的物理基础，在特定降雨情景下可以对流域物质运移过程进行模拟，并生成流域物质侵蚀堆积空间分布格局。其核心算法使用元胞自动机，每一时刻的状态受相邻的单元的影响，体现了时间和空间上相邻区域相互影响的思想。该模型可以根据不同降雨情景，计算相应地表径流的流速和深度。其最突出的功能是能够借助二维水动力模型较好的描述河流动力学特征和泥沙的运输过程，进而估计河床侵蚀和河道变化。
6.模型的输入主要包括地形、降雨、土壤粒径分布、植被、土壤蠕移速率、边坡失稳临界值等参数。模型输出包括洪水淹没水深及范围、流域出口产沙量以及地形高程变化等。在模型中，研究流域被划分为一系列大小一致的网格单元，每一个网格单元都存储着该对应位置的高程、植被条件、粒径分布、流量、水深、流速等数值。在模型运行过程中，每个网格存储的数值根据周围相邻的网格单元通过一系列规则和算法实时迭代更新，进而对流域不同位置各地表要素进行演算。这些运算规则大致可以分为四个模块：(1)水文模块；(2)河流动力模块；(3)泥沙运移模块；(4)坡面模块。
7.研究发现，caesar
‑
lisflood模型虽然也包含坡面模块，能对流域坡面的失稳和物质蠕移进行一定程度的模拟。但模型过于简化坡面失稳过程，不能有效的体现滑坡触发模
拟的空间异质性。且模型输出结果中，流域出口产沙量不能区分滑坡产沙和非滑坡产沙，因此模型不能实现对滑坡产沙及其运移进行追踪模拟，从而无法量化滑坡对流域总产沙量的贡献及定量影响。
8.本发明基于caesar
‑
lisflood二维地表景观演变模型，通过对模型坡面模块进行改进，并向模型添加地表泥沙(物质)运移追踪方法，能够实现滑坡触发模拟及其产生物质的运移过程追踪，并由此实现了本发明。
9.根据本发明的一方面，提供一种基于caesar
‑
lisflood模型的流域泥沙溯源模拟方法，该方法包括：
10.s1，获取流域内的地质，地貌以及水文要素，并且基于这些要素来确定流域的滑坡失稳坡度角阈值；
11.s2，利用滑坡失稳坡度角阈值，建立流域滑坡失稳坡度角阈值空间分布文件，并且利用空间索引将caesar
‑
lisflood模型与该空间分布文件建立联系；
12.s3，建立二维物质源索引矩阵，对不同物质源区对应的空间范围进行初始定义，并且将caesar
‑
lisflood模型与二维物质源索引矩阵建立联系；
13.s4，基于caesar
‑
lisflood模型，建立代表地表活动地层系统的四维数组变量strata[xyindex；z；n，t]以及代表泥沙粒径组成的三维数组变量grain[xyindex；n，t]，其中，xyindex代表位置索引，z存储地层索引，取值范围为1
‑
10，代表10个活动地层；n存储对应泥沙粒径范围索引，取值范围为1
‑
9，代表9个不同泥沙粒径范围，t代表追踪索引值，取值范围1
‑
m，m为预定义的追踪索引值个数，大于1的自然数，代表m个不同物质源区；以及
[0014]
s5，利用caesar
‑
lisflood模型模拟计算各物质源区物质运移并溯源不同物质源区物质，其中包括各物质源区在流域出口逐日产沙量以及各物质源物质运移空间分布。
[0015]
根据本发明的实施方案，s1中，地质地貌以及水文要素包括土壤覆盖层特征，地质岩组类型，岩体结构，地面坡度以及地下水分布。
[0016]
根据本发明的实施方案，s2中，所述空间分布文件为与流域数字高程模型(dem)文件区域范围和空间分辨率相一致的二维矩阵文件，每个矩阵元素存储着对应位置坡面失稳的坡度角阈值。
[0017]
根据本发明的实施方案，s3中，建立二维物质源索引矩阵包括：
[0018]
预定义物质源追踪索引个数m，分别代表m个不同物质源区；以及
[0019]
建立与数字高程模型(dem)文件区域范围和空间分辨率相一致的二维矩阵文件，矩阵索引值个数与预定义物质源追踪索引个数m一致。
[0020]
根据本发明的实施方案，s5包括：
[0021]
caesar
‑
lisflood模型模型读取流域滑坡失稳坡度角阈值空间分布文件，计算滑坡失稳区域范围并触发滑坡，计算由于滑坡发生而导致的物质运移；以及读取所述二维物质源索引矩阵，追踪不同物质源区物质运移过程。
[0022]
本发明基于caesar
‑
lisflood二维地表景观演变模型，通过对模型坡面模块滑坡触发模拟进行改进，并向模型添加地表泥沙(物质)运移追踪方法，实现滑坡触发模拟及其产生物质的运移过程追踪。改进后的坡面模块可以通过输入一个二维矩阵存储流域不同位置对应的坡面失稳角度阈值，实现对空间不同区域失稳坡度角阈值的滑坡触发进行模拟。同时，模型新增的流域泥沙(物质)追踪方法通过输入一个二维索引矩阵存储不同物质源区
对应的索引值来标记对应空间位置的泥沙(物质)，从而实现对不同泥沙(物质)源区的泥沙(物质)运移过程进行追踪。改进后的模型能更加有效的模拟不同区域的滑坡触发机制，且能实现对滑坡物质运移过程的追踪，例如能够定量区分模拟结果中滑坡产沙及非滑坡产沙，为研究滑坡对流域总产沙量及流域泥沙(物质)运移时空规律的影响提供有效的模拟工具及技术思路。
附图说明
[0023]
图1是根据本发明一个实施方案的基于caesar
‑
lisflood模型的流域泥沙溯源模拟方法的工作流程图；
[0024]
图2是根据本发明一个实施方案的流域滑坡失稳阈值空间分布的示意图；
[0025]
图3是根据本发明一个实施方案的基于caesar
‑
lisflood模型的流域泥沙溯源模拟方法的流域侵蚀堆积分布图(滑坡物质分布)模拟结果示意图；
[0026]
图4是根据本发明一个实施方案的物质源区追踪指数值预定义空间分布示意图；
[0027]
图5是根据本发明一个实施方案的基于caesar
‑
lisflood模型的流域泥沙溯源模拟方法的滑坡产沙(物质)追踪模拟结果示意图。
具体实施方式
[0028]
为清楚的说明本发明中的方案，下面给出优选的实施例并结合附图详细说明。以下的说明本质上仅仅是示例性的而并不是为了限制本公开的应用或用途。应当理解的是，在全部的附图中，对应的附图标记表示相同或对应的部件和特征。
[0029]
图1是根据本发明一个实施方案的基于caesar
‑
lisflood模型的流域泥沙溯源模拟方法的工作流程图。参考图1，该实施例的基于caesar
‑
lisflood模型的流域泥沙溯源模拟方法可以包括：
[0030]
s1，获取流域内的地质，地貌以及水文要素，并且基于这些要素来确定流域的滑坡失稳坡度角阈值。例如地质地貌以及水文要素可以包括土壤覆盖层特征，地质岩组类型，岩体结构，地面坡度以及地下水分布等等。可以通过现场勘探、现有历史资料等来获取这些资料。然后，例如可以基于历史滑坡数据库的坡度角阈值统计法、基于安全指数法的阈值坡度计算等已知的方法来确定滑坡失稳坡度角阈值，本领域技术人员可以根据不同流域的具体情况来确定。
[0031]
s2，利用滑坡失稳坡度角阈值，建立流域滑坡失稳坡度角阈值空间分布文件，并且利用空间索引将caesar
‑
lisflood模型与该空间分布文件建立联系。
[0032]
更具体地，坡度角阈值空间分布文件，也即，二维坡角阈值索引文件为与流域dem(数字高程模型)文件区域范围和空间分辨率相一致的二维矩阵文件，每个矩阵元素存储着对应位置坡面失稳的坡度角阈值。通过在模型运行窗口坡面模块添加一个文本输入框，用于输入带有文件后缀的索引文件名。在模型读取文件并初始化时，将索引文件写入模型并转换成二维矩阵，并通过空间索引将坡面失稳坡度角阈值与模型其他数据(例如地形数据，用来计算流域坡度角)建立联系。
[0033]
s3，建立二维物质源索引矩阵，对不同物质源区对应的空间范围进行初始定义，并且将caesar
‑
lisflood模型与二维物质源索引矩阵建立联系。
[0034]
更具体地，首先可以预定义物质源追踪索引个数m，分别代表m个不同物质源区。然后建立与dem(数字高程模型)文件区域范围和空间分辨率相一致的二维矩阵文件(索引个数可以自定义；如果物质源只分滑坡区域和非滑坡区域，则索引个数为2，标记索引值为1，2；如果定义三个物质源区，则索引个数为3，标记索引值为1，2，3；依此类推)。通过在模型运行窗口文件输入模块添加两个文本输入框，一个用于输入预定义索引个数，一个用于输入带有文件后缀的追踪索引文件名。在模型读取文件并初始化时，将追踪索引文件写入模型并转换成二维矩阵，矩阵每个值对应该位置的物质源标记索引值，也即，矩阵索引值个数与预定义物质源追踪索引个数m一致。
[0035]
s4，基于caesar
‑
lisflood模型，建立代表地表活动地层系统的四维数组变量strata[xyindex；z；n，t]以及代表泥沙粒径组成的三维数组变量grain[xyindex；n，t]，其中，xyindex代表位置索引，z存储地层索引，取值范围为1
‑
10，代表10个活动地层；n存储对应泥沙粒径范围索引，取值范围为1
‑
9，代表9个不同泥沙粒径范围，t代表追踪索引值，取值范围1
‑
m，m为预定义的追踪索引值个数，大于1的自然数，代表m个不同物质源区。
[0036]
更具体地，在caesar
‑
lisflood模型中，模型通过建立一个活动地层系统(strata)来代表地表泥沙(物质)组成。活动地层系统包括最上面的一个顶部活动地层，多个下埋活动地层和一个基岩层。模型能对每个活动地层定义多达9个不同粒径大小的泥沙颗粒组成，技术人员可以根据具体情况来自行设定粒径范围个数和数值大小。同时模型中用泥沙组成变量(grain)来存储九种不同粒径大小泥沙体积占相应地层体积的比例，用来量化流域不同位置不同粒径泥沙量及后续运移量。模型通过迭代演算地层系统以及对应泥沙组成模拟流域不同位置泥沙(物质)运移结果，进而计算流域地形变化量及产沙量。地层系统及泥沙组成具体定义如下：
[0037]
定义grain[xyindex；n，t]变量为三维数组。其中，xyindex代表位置索引，为对应模型任意位置[x，y]的唯一值。n代表对应泥沙粒径范围索引，取值范围为1
‑
9，代表9个不同泥沙粒径范围；t代表追踪索引值，取值范围1
‑
m，m为预定义的追踪索引值个数，代表m个不同泥沙(物质)追踪源区；
[0038]
定义strata[xyindex；z；n，t]变量为四维数组。其中，xyindex代表位置索引，z存储地层索引，取值范围为1
‑
10，代表10个活动地层；n存储对应泥沙粒径范围索引，取值范围为1
‑
9，代表9个不同泥沙粒径范围；t代表追踪索引值，取值范围1
‑
m，m为预定义的追踪索引值个数，代表m个不同泥沙(物质)追踪源区。
[0039]
s5，利用caesar
‑
lisflood模型模拟计算各物质源区物质运移并溯源不同物质源区物质，其中包括各物质源区在流域出口逐日产沙量以及各物质源物质运移空间分布。
[0040]
更具体地，caesar
‑
lisflood模型读取流域滑坡失稳坡度角阈值空间分布文件，计算滑坡失稳区域范围并触发滑坡，计算由于滑坡发生而导致的物质运移。基于读取的二维索引矩阵，获得每个位置对应的坡面失稳坡度角阈值。在每个模型运行步长内基于dem(数字高程模型)计算每个位置坡度角后，与对应位置的坡面失稳坡度角阈值比较。如果该位置坡度角小于失稳坡度角阈值，则该处坡体处于稳定状态，不发生滑坡。如果该位置坡度角大于失稳坡度角阈值，则触发滑坡。该过程会重复迭代，直至流域各区域的坡度角小于或等于失稳坡度角阈值。然后，基于caesar
‑
lisflood模型的坡面过程和泥沙运移过程，并读取所述二维物质源索引矩阵，来模拟并追踪不同物质源区物质运移过程，并输出例如每个物质
源区生成一个日产沙量数据文件；各物质源区对应泥沙(物质)运移空间分布文件，每个泥沙(物质)源区对应一个相应的空间分布输出文件等等。
[0041]
下面结合实例，对本发明的具体实施方式作进一步详细描述。以下实例用于说明本发明，但不用来限制本发明的适用范围。
[0042]
本研究选取四川省某一流域为研究区，流域面积约179平方千米。2008年发生地震期间，流域共发生滑坡188起，滑坡产生的大量碎屑物质堆积在山间汇水山谷，在降雨的作用下滑坡碎屑物质顺着河谷向河道运移，导致河床淤升，河道发生改变。
[0043]
本研究通过统计研究区地震期间滑坡失稳坡度角阈值，建立流域滑坡失稳坡度角阈值空间分布文件。并基于不同滑坡空间范围定义的追踪索引值，建立基于不同滑坡位置的物质源区追踪索引值的空间分布文件。基于改进后的caesar
‑
lisflood模型，还原当时滑坡触发情景，并对滑坡物质运移过程进行追踪模拟，探索滑坡发生后(2008年
‑
2100年)对研究流域出口产沙量以及流域地表(地貌)演变的影响。
[0044]
本研究模型主要输入数据包括：初始地形数据(dem)，降雨数据(2008
‑
2100年gcm全球模式预测数据)，泥沙粒径组成比例数据(也即，泥沙级配数据，通过在研究流域不同地点实地采样土样，然后在实验室通过筛分法计算得到，具体见下表1)，流域坡面失稳坡角阈值空间分布文件以及流域物质源区追踪索引文件；caesar
‑
lisflood模型本身为本领域技术人员所熟知，例如运行所需数据及参数等等，因此，本文在此不作详细阐述，而是着重对本发明改进进行详细阐述。
[0045]
表1：泥沙粒径级配数据
[0046][0047]
本研究通过统计研究区地震期间滑坡失稳阈值，建立流域滑坡失稳坡度角阈值空间分布文件。此处滑坡失稳坡度角阈值取值范围为10
°
至65
°
，其中非滑坡区域坡度角阈值设置为65
°
(参见图2)。
[0048]
基于生成的滑坡失稳坡度角阈值空间分布文件，获得每个位置对应的坡面失稳坡度角阈值。计算滑坡失稳区域范围并触发滑坡。
[0049]
流域在滑坡影响下的侵蚀主要分布在滑坡影响区，沉积主要分布在滑坡底部和主要下游河道。其中滑坡侵蚀最大深度超过4米，滑坡堆积最大深度可达10米。应该注意的是，滑坡扰动集中在下游地区，对上游的高程变化影响较小。遭受大面积侵蚀和沉积的区域主要集中在研究流域的下游，这与滑坡发生的位置一致。山体滑坡的发生对流域地形演变有重大的影响(参见图3，正值(红色)代表侵蚀，负值(蓝色)代表堆积)。
[0050]
模型读取二维物质追踪索引文件，并预定义物质源追踪索引个数。基于2008年该地区地震引发的滑坡情景，本案例在研究流域内一共定义21个追踪指标，分别代表21个物质追踪源区(参见图4)。除三个子流域的追踪源区(追踪指数：1
‑
3)外，滑坡所在区域也被设置为追踪源区(追踪指数：4
‑
21)。其中，在前16个最大的滑坡中，每个滑坡都标有一个追踪指数(追踪指数：4
‑
19)，以观察单个滑坡的泥沙输送情况。其余较小的滑坡分为两组，分别位于子流域1和子流域2，用另外两个追踪指数标记(追踪指数：20
‑
21)。
[0051]
在原始模型基础上，定义新的地层系统和泥沙(物质)组成等数组变量。结合模型原有运行框架，对滑坡等各物质源区的物质运移进行演算模拟，模拟时间尺度为93年(2008
‑
2100)，在模拟结束时保存各个追踪物质源区对应的输出文件。
[0052]
滑坡产沙(物质)追踪结果统计图(图5)展示了滑坡的空间分布以及相应泥沙(物质)运移状况。每个滑坡的饼状图显示了山体滑坡产生泥沙(物质)的不同组成比例，包括就地沉积停留在原地的泥沙、输送到河道的泥沙和输出流域出口泥沙的比例。显然，不同的滑坡体之间的总产沙量以及不同运移状态的泥沙所占的比例都显示出很大的差异。与河道网络相连或相交的滑坡有较高比例的泥沙被输送到主河道和输出流域出口。而位于远离河网的山坡上的滑坡有较少的泥沙输出流域出口。值得注意的是，1号滑坡和2号滑坡是面积最大的两个滑坡，其大小规模相近。而2号滑坡虽然离流域出口更远，但滑坡产生的泥沙的输送比例更高。这是因为2号滑坡与径流相交，滑坡物质运移能力大大增强。而1号滑坡只有很少一部分滑坡体有径流经过，导致滑坡体物质运移速率较低。上述模拟结果与地震后实地勘测结果基本一致，反映了本发明方法的科学可靠。
[0053]
本案例基于河流域滑坡失稳角统计值对流域坡面失稳临界角度预定义，对流域内2008年地震同震型滑坡进行触发还原。同时在全球变化背景下，采用基于全球气候模式(gcm)的2008年
‑
2100年未来降雨数据作为模型水文输入。运用改进的caesar
‑
lisflood模型对流域长时间尺度(百年尺度)地表演变以及滑坡物质追踪进行模拟。结果显示，改进后的模型输出能较好的还原由于滑坡造成的地形变化以及长时间尺度下流域地表演变。同时，改进后的模型能对用户自定义的特定滑坡产生的泥沙(物质)进行追踪模拟，解决了实际运用研究中无法识别河道以及流域出口泥沙来源的难题。为流域泥沙管理以及滑坡风险管控提供了新的工具以及技术手段。
[0054]
综上所述，以上所述内容仅为本发明的实施例，仅用于说明本发明的原理，并非用于限定本发明的保护范围。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。