基于小区域精确淹没算法的地道桥淹没分析与动态模拟

张文胜， 解 骞， 谢进通， 宋金辉， 张春娈， 肖洪生

(1.石家庄铁道大学交通运输学院，河北石家庄 050043;2.河北省交通安全与控制重点实验室，河北石家庄 050043;3.衡水龙翔公路工程勘察设计咨询有限公司，河北衡水 053000)

地道桥是城市道路下穿其它交通设施形成的一种凹型路段。由于地道桥最低处常比周围低4～5 m，形成一大片地势低洼区域，暴雨时极易发生积水。行人、车辆在不知水深的情况下冒然通过，容易造成车辆熄火，对车辆和人员造成威胁。因此，精确模拟地道桥暴雨淹没过程，对防洪减灾和交通应急有重要的意义。20世纪90年代开始发展起来利用GIS技术进行基于地形的洪水淹没分析，但该技术还只是停留在二维阶段，21世纪初开始以数字高程模型(DEM)为基础分析淹没区域。刘仁义等利用种子蔓延法对给定水位的条件下进行有源淹没和无源淹没分析［1］，开辟了国内使用DEM分析淹没范围的先河;葛小平等根据水力演进模型进行洪水的模拟预测，采用分区平面模拟方法模拟淹没范围的提取及显示［2］;张书亮从管网的排水能力角度计算城市汇水区域，为雨水管网设计提供技术参考［3］;丁志雄等将DEM转成三角形和任意多边形网格，在给定水位或水量两种情况下，利用网格的联通性算出淹没区域，研究出计算淹没区域的新思路［4］;Aiman Abdel-Lattif运用GIS和遥感技术分析出径流网络，研究出某地区发生洪水时的风险区域和危险等级［5］;姜仁贵等将DEM转成Grid格网，根据每个格网平均高程与水位高的差计算积水量，并给出三维仿真显示［6］;国外知名软件暴雨洪水管理模型(SWMM)主要是对城市降雨的水量和水质的模拟，可分析每个管道和河水的水量、水质情况，但没有对淹没区综合的分析和显示功能。目前已有的淹没分析都是大区域、低精度、静态的汇水分析，针对某一具体小区域进行动态的淹没算法分析与模拟的研究较少。本文针对地道桥这一小区域的暴雨过程，利用淹没算法对积水区域及深度进行分析与动态模拟。以GIS的空间分析为基础，采用种子蔓延法和流域D8算法结合确定汇水区域，根据地道桥的长宽及车道数等动态生成三维模型，在给出降雨强度的情况下动态模拟水淹过程并显示水深。

1 小区域精确淹没算法

暴雨刚开始时雨水会被干燥的地面吸收，一段时间后地面吸水饱和，雨水开始在地面形成径流，高处的水会往低处流。地道桥区周围一片区域地面水会往桥内汇，在桥底部形成积水，随着降雨的进行，积水会不断的上涨，积水区域也会扩大，因此研究小区域的汇水范围及汇水面积与深度的关系是小区域精确淹没算法的关键。

1.1 汇水区域范围确定及算法改进

地道桥处地势低洼，暴雨时周围区域雨水往桥内汇集形成汇水区域，汇水区域的确定对于暴雨淹没的分析与计算至关重要，而确定汇水区域是一个复杂的过程，有以下几种情况:桥建在洼地，四周的水都往桥内汇;桥建在坡上，只有一面区域水会往桥内汇;桥建在地势高处，只有地道桥区域的水会往里汇，周边区域的降水不会汇入。

有源淹没是指洪水从某点开始淹没，然后根据DEM推算出淹没区域。地道桥的汇水区域可以用洪水淹没中有源淹没的逆向过程去推算。目前一般采用种子蔓延法［7］来推算出淹没区域。种子蔓延法是指在洪水的源头“种下”种子，种子开始向四周蔓延，只要高程比源头低的区域就会蔓延到，反之，则停止向这个方向蔓延。如果淹没区域是一片洼地，这样可以确定种子只在洼地内区域蔓延，不会蔓延到洼地外比洪水源头高程低的区域。种子蔓延法是一种无序的蔓延，若仅仅把地道桥汇水考虑为比桥高的区域降水都会流入桥区，这仅仅是桥建在洼地的最低端是成立的，其它地方都不成立。比如桥建在一面斜坡上，将会造成很多比桥区高的桥两旁区域都会被计入汇水区域。所以必须把种子蔓延定为有向的蔓延，种子蔓延的方向为降水水流的逆方向。

水流方向算法中最经典的是D8算法，如图1所示。该算法简单，易于实现，且执行效率高［8］。D8算法是一种单流向算法，假设单个栅格中的水只从一个方向流入周围8个方向的栅格中任意一个，它用最陡坡度法确定水流方向，即这8个方向中比该栅格高程小，且落差除以栅格中心距得数最大的方向即为水流出方向。

图1 D8算法示意图

但是暴雨时，不是桥周围所有比地道桥高的区域都会往桥里汇水，这种情况只是桥建在地势最低处，桥周围地势都比桥区高时成立。当桥建在坡上或在地势高处并不成立，所以不能单纯的使用种子蔓延法。为此本文对其进行了改进，改进后的种子蔓延算法提高了汇水区域计算准确度，更能符合实际地形汇水特征。具体改进方法是，在蔓延的路径上加上一个限定的条件，根据已建出的地道桥模型，在其边界线上等距离的取一系列点，放入最初始种子点集seedlist中，从seedlist中取出一个种子，将DEM抽象成网格状，如图2所示，这粒种子在网格的中心，然后向周围8个方向寻找种子的蔓延方向。以先向E方向寻找为例，用D8算法判断E格水的流向，若E格的水向Seed格流，种子向E格蔓延，将E格中种下种子，即把E点放入seedlist中。若8个方向中没有一个方向的水流入Seed格中，则将Seed放入边界点集boundarylist中，确定这点为汇流区域的边界点。收敛条件是当seedlist中再没有种子，则蔓延结束。将边界点集中的点按一定顺序连起来形成汇流区域。

图2 种子蔓延示例

1.2 汇水量计算方法

城市地表大多由沥青路面、水泥路面、屋顶这些不透水区域和少量绿化草地构成，容易形成积水。汇水量可以考虑为降雨积累量扣除初期损失和累计下渗量得出［9-10］，即

图3 地道桥纵断面汇水量及淹没深度变量图

式中，R为某一时刻汇水总量;P为某一时刻的累计雨量;Ss为初期损失，以地表滞蓄量和植物截留量为主;Fm为累计下渗量。

城市地表的产流是指降雨经过植物截留、下渗、蒸发后，产生地表径流的部分降雨，地表的产流是个动态的过程，在降雨刚开始时滞留、下渗量较大，随着降雨的持续，下渗量越来越小［11］。地表的产流采用变径流系数法计算

式中，Ψ为某一时刻的径流系数;Ψe为最终径流系数;Ψ0为初始径流系数;P为累计雨量;c为常数。

确定最终汇水量公式为

式中，S为汇水区域的面积;i为降雨强度;t为降雨时间。

1.3 淹没深度的计算方法

已知积水体积的情况下，要确定淹没深度，需找出地道桥的积水体积与淹没深度之间的关系。可以先根据地道桥的引道长度、道路的等级模拟出地道桥纵断面的曲线方程［12］。地道桥的纵断面可抽象为由3段圆曲线和2段直线表示，根据道路等级、停车视距计算出最底部的凹曲线和两边的凸曲线的曲线半径及投影长度，直线的斜率和长度则可由两边曲线确定，假设地道桥两边对称，曲线方程为

同样可计算出图3中纵断面淹没面积S与淹没面积投影长度x0的关系式

这样可知在x已知的情况下，地道桥积水量为

式中，S1为机动车道的纵断面淹没面积;w1为机动车道的宽度;S2为非机动车道的纵断面淹没面积;w2为非机动车道的宽度。

通过动态的输入降雨强度，根据公式(3)可计算出地道桥周围的汇水量R，由汇水量得出淹没深度是一个不断迭代的过程。若直接令积水量V等于汇水量R，由V反解出x0，最后将x0带入式(1)中计算出水面高度，这样看似简单，其实实现起来很困难，因为积水量V是分段复杂函数，求复杂函数的反函数十分繁琐。为此，采取以下算法予实现，如图4所示，设f(x0)=V－R，在R确定的情况下，f(x0)为单增函数且必有零点，可以采用二分法逼近零点。x0属于区间(0，L)，零点最初在区间(0，L)内，如果 f(x0)≈0，则 L/2就是零点;如果f(x0) ＜0，零点在区间(L/2，L)内，将函数区间改为(L/2，L);如果f(L/2)＞0，零点在区间(0，L/2)内，将函数区间改为(0，L/2)，取区间的中间值继续进行迭代。通过每次迭代将f(x0)所在的区间缩小一半，使区间的端点逼近函数的零点，求出零点的近似值。最后将得出零点值的xp带入式(4)中计算出水面高度，并实时更新水面高度。通过上述算法，简化了函数的计算过程，使计算效率得到了较大的提高。

图4 积水深度计算流程

2 基于三维GIS的地道桥淹没动态模拟及算法实现

根据本文提供的淹没算法及建模方法，采用C#编程语言和SQL Server数据库平台，二次开发三维GIS软件skyline，实现地道桥淹没分析与动态模拟。skyline提供了由坐标可在三维场景中自由画线、画面和构造简单体的功能，但它不能直接构造像地道桥这类特殊形体的建筑。三维空间建模技术可分为基于面模型的准三维建模和基于体模型的真三维建模［13］。面模型是用点、线、面表示实体的表面，如地形表面，地质层面;体模型是用体元分割和体元模型比如四面体、棱柱体进行实体表达。地道桥属于地面的交通设施，将它与地形一体化建模可与地面更好的融合，因此可以采用构造地形的方法构造地道桥模型。地道桥地形简单规整，本文采用基于面模型构模方法中的Grid模型，用规则的格网构造地道桥，如图5所示。格网的划分可通过将地道桥从引道头到引道尾分成等距段，每段的分成中间机动车道和两边非机动车道三个格网。调用skyline二次开发接口IobjectManager5中Create2Dpolygon画面方法，由上述的曲线方程(4)可得格网上每点的高程，根据每点的三维坐标画出方格，最后由方格拼出整个地道桥。由构造出的地道桥坐标及地面高程，通过上述改进的水文学算法确定汇水区域范围。根据实时输入的降雨强度和控制模拟速率，以及上述的淹没深度计算方法得到实时的淹没深度和淹没面积，用时间组件Timer实现动态调整积水水位和积水深度值的显示，最终淹没模拟如图6所示。算法实现过程具体如下:

图5 基于Grid的构造示意图

图6 基于动态算法的地道桥淹没模拟

(1)指定对某一地道桥进行暴雨模拟。通过给出的地道桥初始信息构造出地道桥三维模型，并由上述的汇水区域范围确定算法计算出该地道桥的汇水区域。

(2)动态输入降雨强度。通过动态输入的降雨强度计算出累计降雨量，由上述的汇水量计算方法用降雨量和降雨时间计算下渗量，降雨量减去下渗量可得实时的汇水量。

(3)由汇水量和地道桥立体体积计算淹没深度。通过地道桥的初始信息可知地道桥立体的体积，采用二分法逼近汇水量从而算出淹没深度。

本文对石家庄某地道桥进行暴雨淹没模拟，该桥桥长560 m，最大深度5 m，有4条机动车道和2条非机动车道，由于地势低洼，暴雨时经常被淹。将某一历时3 h暴雨的降雨强度数据动态输入程序中模拟地道桥淹没过程，结果如图7、图8所示。结果表明暴雨刚发生时，由于地面干燥，降雨量与下渗量基本相等，地面积水极少，当降雨进行一定时间后，地面含水量达到饱和，下渗量开始慢慢减少，地面开始形成积水，20 min后桥底积水深度已达到0.467 m，机动车不能通行。随着降雨的继续进行，积水越来越深，当降雨2 h后，积水已漫过非机动车道，淹没深度达到了1.967 m，非机动车和行人不能通过。通过和现实暴雨时地道桥被淹情况的对比，模拟过程基本与现实相符。

图7 降雨强度历时数据

图8 实时淹深度

3 结论

为了预防城市暴雨造成地道桥积水给行人和交通带来事故灾害，改进经典水文学算法结合GIS的空间分析功能研究出小区域精确淹没分析算法，对城市暴雨地道桥进行淹没分析及动态仿真模拟。着重考虑了根据DEM确定地道桥汇水区域的方法，解决了小区域汇流区域面积不准确的问题。并运用迭代逼近算法得出汇水量，动态给出淹没区域及积水深，为小区域精确淹没分析提供了步骤及方法。在暴雨发生时，道路上的雨蓖常常被垃圾、树叶堵塞，导致排水设施基本失效，所以在本文中将排水量考虑在下渗量之中，更能贴近实际。

［1］刘仁义，刘南.基于GIS复杂地形洪水淹没区计算方法［J］.地理学报，2001，56(1):1-6.

［2］葛小平，许有鹏.GIS支持下的洪水淹没范围模拟［J］.水科学进展，2002，13(4):456-460.

［3］张书亮，曾巧玲，姜永发，等.GIS支持下的城市暴雨积水计算的可视化［J］.水利学报，2004，12:92-98.

［4］丁志雄，李纪人，李琳.基于 GIS格网模型的洪水淹没分析方法［J］.水利学报，2004，6:56-60.

［5］Aiman Abdel-Lattif，Sherif Y.Morphmetric analysis and flash floods of Wadi Sudr and Wadi Wardan，Gulf of Suez，Egypt:using digital elevation modal［J］.Arab J Geosci，2012，5:181-195.

［6］姜仁贵，解建仓，李建勋，等.基于数字地球的洪水淹没分析及仿真研究［J］.计算机工程与应用，2011，47(13):219-222.

［7］赵秀英，王耀强.基于DEM的有源淹没算法设计与实现——以种子蔓延法为例［J］.科技导报，2012，30(8):61-64.

［8］邬伦，汪大明，张毅.基于DEM的水流方向算法研究［J］.中国图象图形学报，2006，11(7):998-1003.

［9］杨戈，吴升.城市暴雨积水模拟方法分析及研究［J］.测绘信息及工程，2009，34(1):35-37.

［10］Yang Liu，Purius Soonthornnonda.Stormwater runoff characterized by GIS determined source areas and runoff volumes［J］.Environment Management，2011，47:201-217.

［11］任伯帜.城市设计暴雨及雨水径流计算模型研究［D］.重庆:重庆大学，2004.

［12］赵静雅.地道桥引桥纵断面设计浅析［J］.北方交通，2010，7:14-17.

［13］吴立新，史文中.地理信息系统原理与算法［M］.北京:科学出版社，2003:48-51.