基于流域GIS空间建模的雨洪模拟方法

2022-03-11徐建刚

地理空间信息 2022年2期

徐晗，徐建刚

（1.南京大学建筑与城市规划学院，江苏南京 210093；2.上海营邑城市规划设计股份有限公司，上海 200030）

通过近十年来的不断探索，中国海绵城市规划工作开始从流域入手，从流域整体排洪防涝来进行规划设计。为了更科学地确定海绵城市大中型海绵体的选址、滞水容量与淹没范围，需要引入水文模型对河道断面在暴雨情景下通过的流量进行模拟与预测。伴随着空间信息科学与技术的进步，基于流域栅格化技术的雨洪时空模拟方法逐渐发展起来。运用数学与物理方程对水文循环中的各个过程进行解释与描述，通过参数与变量来表达其空间分布上的差异性[1]。这种模拟现实物理法则的物理模型能够更准确地描述水文循环过程，并且具有很强的适应性[2]。因而更适用于对空间精度要求较高，空间差异性敏感的城市规划领域。本文在林蔚[3]等提出的暴雨汇流水文模型集成分析的基础上，结合海绵城市规划的实际设计工作需求，面向海绵城市的应用进行优化。通过对水文学、水力学相关模型的研究，基于多层次GIS空间分析模型技术，系统地建立了从流域雨洪过程相关影响因子栅格化处理、到产汇流过程模型的集成化模拟方法，首次实现了对流域任意栅格单元内河道断面的径流过程线的可视化测算。同时，还对于模型中单元水流长度与河道单元流速计算方法进行了优化改进。并以闽西国家历史文化名城长汀县城汀江观音桥水文站上游流域为实验区，使用水文站实测数据对模型进行对比验证，实现了通过对多个情景的模型参数可视化调整来提升模拟精度的技术应用。

1 分时段产流过程模拟

在流域径流形成过程中，流域内的降水经过植被截留、下渗、蒸发等损失而产生净雨过程即为产流过程。目前国内主要采用的产流模拟大多基于GIS平台对流域空间数据进行栅格化表达，进而构建了美国农业部水土保持局在上世纪50年代提出的SCS-CN流域水文模型，通过对暴雨过程中各个时段流域空间产流量的累积计算，以获得流域分时段产流量栅格数据。SCS-CN模型综合考虑了降雨量、土壤透水性、土地利用类型、前期土壤湿度等因子与流域产流的关系，

通过径流曲线数（CN）来推求最大可能滞留量（S），CN是一个经验性无因次参数，是对降雨前流域下垫面土壤类型、土地利用方式、前期土壤含水量所造成影响的一个综合反映。

公式如下：

式中，S为最大可能滞留量（mm）；CN为径流曲线数。

美国农业部土壤保持局在分析了大量长期的实验结果基础上，经过一系列推导得到以下计算产流量公式：

式中，R为实际产流量（mm）；P为降雨量（mm）；S为最大可能滞留量（mm）。

基于上述模型原理与国内相关研究成果，本研究结合实际海绵城市规划项目中所能获得的相关资料，对产流量的影响因子进行梳理，并运用空间插值、地图代数等GIS技术实现了参数的全部栅格化，最终形成产流过程模拟计算的系统性多层次概念模型（图1）。

图1 产流过程模拟概念模型

1.1 降雨量（P）栅格化

根据图1流程所示，降雨量指标值的空间分布式计算是模拟的基础。降雨资料通常为散布于流域内的多个雨量监测站点位数据，可利用IDW插值、泰森多边形等方法栅格化为覆盖全流域的降雨量栅格数据。

在实际应用中，雨量监测站通常记录的是各单位时段内的降雨量（mm）。由于SCS-CN模型中对产流量计算公式是基于一场完整降雨，公式（2）中的初损量0.2 s不会随着降雨时长变化。这使得在计算分时段降雨所形成的径流时需要将每个小时降雨都看作一场完整降雨，并分别运用公式（2）进行计算。则有，对于历时m小时的降雨场次α，其中第n小时（1≤n≤m）的累计降雨量（Pn）为n小时降雨量之和，公式如下：式中，Pn为第n小时降雨的累计降雨量（mm）；pi为第i小时降雨的降雨量（mm）。

1.2 径流曲线数（CN）栅格化

径流曲线数CN值在流域空间的栅格化推算是模拟的关键，本研究从流域下垫面三大因子特征量化、坡降影响修正和土壤前期湿度分级三方面探析，建立了适用于实验区的参数及其科学计算方法。首先，通过查阅美国农业部在网络发布的说明文件《National Engineering Handbook Hydrology Chapters》中的表格，可以获取到在不同土地利用类型、土壤透水度与植被覆盖度组合下的CN值。土地利用类型因子主要为流域土地利用现状图；土壤透水性因子则为流域土壤的HSG分级的空间分布数据，其包括了A/B/C/D 4个级别的指标，透水性由A至D依次降低；植被覆盖率基于多波段卫星影像计算获得的归一化植被指数（NDVI）栅格数据，根据杨胜天[5]等的研究内容计算获得，并依据流域所在地的植被特征，重分类为好、中、差3个等级的栅格数据。最后将三大因子的空间数据运用叠加标识、栅格重分类等GIS技术，进而生成平均土壤湿度状态下的流域CN值栅格数据。

其次，坡降作为水文学中影响径流的重要因子，在地形变化较大的山地丘陵地区应用SCS-CN模型时还应考虑坡度影响并对径流的影响。本研究采用陈正维[6]的研究成果，对CN值进行修正：

式中，CN为径流曲线数，slp为坡降值。

最后，由于土壤前期湿度会显著影响土壤对降雨的初损吸收量，因此，SCS-CN模型以当次降雨前5 d的雨量为基准来判断土壤前期湿度，其将土壤前期湿度分为AMCI，AMCII，AMCIII三级指标，分别代表干、平均、湿3种状态。本研究采用罗鹏[7]的研究成果对AMC进行修正，以上述经过坡度修正的CN值作为平均状态，其他2种状态的情况则根据下面2个公式进行转换：

式中，CNI为干状态下的CN值；CNII为平均状态下的CN值；CNIII为湿状态下的CN值。

最终获得3种状态下的流域CN值栅格数据，根据暴雨情景前5 d雨量选择使用。

1.3 分时段产流量计算

基于地图代数GIS技术，将前述所获得的Pn栅格化成果及CN值栅格代入SCS-CN模型公式（1）、（2）中作为参数，计算得出降雨n小时后流域的累计产流量（Qn），而降雨在第n小时的产流量（qn）则为第n场降雨累计产流量（Qn）与第n-1场累计产流量（Qn-1）之差：

式中，qn为第n小时降雨的分时段产流量（mm）；Qn为第n小时降雨的累计产流量（mm）。计算获得流域一场降雨中的分时段产流量栅格。

2 汇流过程模拟

汇流过程，即流域净雨流入河道中的聚集过程。流域汇流过程模拟采用空间分布式旅行法模型[8]，模型通过追踪各单元径流汇流至流域出口所经过路径，并应用曼宁公式、连续方程等水力学公式计算获得径流通过各单元所花费的时间成本，进而累加计算各流域单元到达流域出口所要花费的汇流时间，生成流域等流时面栅格，即所产生径流同时达到出口断面的栅格单元数目，定义为该时刻的无因次流量。最终结合产流模型所生成的分时段流域产流量栅格，汇总统计生成流域出口断面流量过程线。这一汇流过程计算模式的概念层次模型如下（图2）。

图2 汇流过程模拟概念模型

2.1 流域等流时面栅格计算

降雨径流在流域进行汇流过程时存在坡面漫流与河道流2种不同的汇流模式，需要将流域划分为坡面单元与河道单元，计算流域单元水流长度并分别计算各单元的汇流成本，进而生成等流时面。因此，计算流域各单元的径流通过时间成本需要划分为5个计算步骤。

2.1.1 坡面单元与河道单元划分计算

本研究引入河道提取阈值来界定坡面与河道两类单元。该域值被定义为形成河道汇流所需要的最小集水面积占流域总面积的比值。基于由流域DEM生成的流域流向栅格，可以通过向流动方向累加的方法生成的无因次流量栅格（累计数为栅格像元个数，图3），则河道单元即是像元值大于“河道提取阈值*流域栅格总数”的栅格像元，其余像元便是坡面单元。由于降雨量的不同会导致最小集水面积的变化，这使得河道提取阈值需要通过实验对比实现河流情况获得。

2.1.2 单元水流长度计算

传统单元水流长度的确定，通常依据栅格像元大小计算，当流向为正南北向或正东西向时，取栅格像元边长大小；当流向为斜对角线方向时，取栅格像元斜对角线长度（栅格像元边长×20.5）。这种方法存在2个问题：①在现实中水流通常沿地形成曲线汇流，这种方法会导致分析用栅格像元越大，所获取值的径流长度相比现实径流长度越小。②这种方法所获取到的仅是路径像元的平面长度，而实际水流是在路径单元的坡面进行汇流，忽略了路径坡度的影响，坡度越大计算结果相较现实径流长度越小。为解决以上2个问题，研究提出以下优化方法。

1）地图精度纠正。如图4所示，不同比例尺地图下同一条汇水线的形态与长度并不相同，受比例尺精度影响，比例尺越小，误差越大。根据国内学者对数字地图曲线长度精度的相关研究[9-10]，本研究采用公式（8）进行径流长度纠正。

图4 基于不同比例尺地形图所生成的汇水线差异

式中，ld为修正后径流长度；li为纠正前径流长度；d为地图精度纠正系数。

通过对河流进行不同栅格单元大小（设单元边长为r）的栅格化处理，并统计获得各个栅格的栅格数N，将测量结果做log(N)～log(r)双对数散点图，拟合得一负斜率直线，该斜率的绝对值为分维值D。作者通过对研究流域内多条主要河流进行上述计算，获得的分维值D约为1.05～1.11，综合分析实验与前人研究的结果，本研究取1.1为实证流域单位水流长度的地图精度纠正系数。

2）坡度纠正。在计算路径像元长度时，可乘以该像元坡度的正割值，进行坡度纠正。可转化为与坡降值（即像元坡度的正切值）相关的公式：

式中，lp为坡度修正后路径像元径流长度；li为路径像元径流长度；slp为路径像元坡降值

2.1.3 坡面单元通过成本计算

坡面单元通过成本时间由稳态运动波近似值及曼宁公式计算获得,计算公式如下[11]：

式中，to为坡面单元时间成本（s）；L为单元水流长度（m）；n为曼宁糙率；ie为平均净雨强度（m/s）；slp为坡降值。

其中平均净雨强度栅格ie表示由累计产流量栅格除以降雨时长获得：

式中，ie为平均净雨强度（m/s）；Rn为第n时段降雨的累计产流量（m）；T为降雨时长（s）。

坡降值栅格由流域DEM栅格，为每个像元计算从该像元到与其相邻的像元方向上的最大变化率，并输出像元值为坡降值的坡降栅格。这样计算获得的非常小的甚至为0的坡降值，这使得部分坡面单元的流速将接近于0，与现实径流情况不符。研究通过设定最小坡降值，去除栅格中过小的坡降值，通常设值为0.01。

曼宁糙率栅格基于流域土地利用栅格生成，通过查阅相关文献标准，赋予对应土地利用糙率数值。

2.1.4 河道单元通过成本计算

河道单元流速公式由宽河道的稳态连续方程结合曼宁公式转化而来，计算公式如下[11]

式中，tc为河道单元时间成本（m/s）；B为河道的有效宽度（m）；n为曼宁糙率；L为单元水流长度（m）；Qe为平均累计汇流量（m3/s）；sl为坡降值。

其中对于河道有效宽度的计算，由于现实中河道宽度与其通过流量成正相关，则对于无因次流量栅格中的河道像元而言，其汇积值越大，该像元的河道有效宽度越大。从而就可以通过计算流域出口断面的河道有效宽度，运用河道的无因次流量栅格推求河道有效宽度栅格。公式如下：

式中，Bi为当前像元河道有效宽度（m）；Bm为流域出口断面河道有效宽度（m）；Qi为当前像元无因次流量；Qm为流域出口无因次流量。

河道平均累计汇流量栅格计算类似无因次流量，采用累计流量法，并采用平均净雨强度栅格作为权重栅格计算获得。

对于河道曼宁糙率值，由于地形图及土地利用图等资料难以真实表达河道底部情况，故在计算时均采用固定数值。河道坡降取值与坡面坡降一致，取坡降栅格值。

2.1.5 汇流时间栅格生成

基于流域流向栅格可追踪获得各栅格像元流向流域出口的唯一汇流路径，结合汇流成本栅格，就可以通过累加计算得到该像元产生径流流至流域出口所花费的总时间，将时间作为该像元的新栅格值，便生成了汇流时间栅格。公式如下：

式中，Tm为流域内某一栅格像元汇流至流域出口的时间；m为汇流路径上的路径像元数量；τi为径流通过某一路径像元所花费的时间。

在计算生成流域汇流时间栅格后，进而通过对时间栅格值按单位时段长度进行分类，通常与降雨数据的分时段长度一致，从而获得流域等流时面栅格（图5），表达在同一时段内汇流至流域出口流域栅格像元。

图5 左：汇流时间栅格；右：等流时面栅格

2.2 流域出口断面流量过程线计算

以产流模型所计算获得的分时段产流量栅格数据作为权重，对流域等流时面栅格进行分区统计，统计各时段等流时面内对应产流量像元值的总和，进而可生成该时段产流在流域出口所形成的流量过程线。最后应用传统单位线法的叠加假定，即出口断面的流量过程线等于各个时段流量过程错开时段叠加之和，对分时段流量过程线错开时段进行叠加，最终生成流域出口断面流量过程线。

3 模拟实证

本研究以长汀县观音桥水文站控制断面以上的汀江上游流域为模拟实证应用示范区，流域面积约为382 km2。以林蔚[3]研究成果为对照组进行方法优化设计与结果精度验证，选择观音桥水文站20 150 519场次洪水，对应重现期十年一遇。该场次为2012－2015年间最大洪水场次，且林蔚[3]研究中的模型效率系数最低、峰现误差时间最高。结果如表1所示,实验设定在其他参数及数据一致的条件下进行。其中对照组不进行任何优化，实验组01仅对单元水流长度进行优化计算；实验组02仅对河道有效宽度进行优化计算；实验组03对二者均进行优化计算，各实验组形成流域出口断面过程线对比图见图6。

表1 模型精度对照表

图6 实验流域出口断面流量过程线对比

根据实验结果可以看出，2种优化措施效果均较为明显，采用优化后的模拟方法得到的结果在模型精度验证的4项指标中均有提升，模型等级由丙级（0.70≥E＞0.50）提升为乙级（0.90≥E＞0.70）。按照城市总体规划一般设定的20～30年期限，该结果可以作为海绵城市规划中的大型海绵体（湿地公园等）控制性指标依据。