来源：第三维度    作者：杨军,贾鹏,周廷刚,张敏,刘微    单位：三峡库区生态环境教育部重点实验          西南大学地理科学学院          中国科学院深圳先进技术研究院          重庆市地理空间信息工程技术研究中心
    摘要:以数字高程模型(DEM)为数据基础,分析了给定洪水水位与洪水水量两种洪水淹没范围计算方法,并利用ArcEngine开发组件实现了淹没区的虚拟现实表达.结果表明:两种计算方法可快速、高效地模拟洪水淹没范围,结合地理信息系统可直观地展示在三维可视化平台上,为防洪指挥调度及洪灾灾害的评估提供更加科学、直观的决策服务,对现代水利的数字化和信息化建设有一定的现实意义。
    我国是一个自然灾害频繁发生的国家,其中洪涝灾害尤为严重,准确的洪水预报和洪水分析可以为防洪应急预案的制定和指挥调度提供重要的科学依据.因此,准确、科学地模拟、预测洪水淹没范围和水深分布对防洪减灾具有重要的现实意义。
    自上世纪90年代以来,在水文水动力模型基础上利用GIS强大的空间分析功能和可视化功能,模拟洪水淹没区,一直是防洪减灾研究热点之一.如刘仁义等人提出了“无源淹没”与“有源淹没”的理论[3];丁志雄等人提出了基于三角形格网模型和任意多边形格网模型的洪水淹没计算方法[4].
    但现有的研究大多数是从理论角度或是基于矢量、三角网等数据形式的二维平面空间技术进行分析讨论,很少对洪水淹没进行分析与计算.随着流域信息化研究的深入和“数字流域”工程的逐步展开,洪水淹没必须上升到模拟仿真的高度,这已经成为洪水淹没分析的趋势之一[5].本文主要以高分辨率数字高程模型(DigitalElevationModel,DEM)为基础,采用高效计算模式,根据不同的洪流产生条件,模拟给定洪水水位和给定洪水水量两种洪水淹没方式,并对这两种淹没方式的算法思想进行阐述,最后应用VisualStudio.NET和ArcEngine平台对这两种洪水淹没计算结果进行了虚拟现实可视化及其三维仿真.
    1 洪水淹没分析方法
    洪水淹没是一个很复杂的动态变化过程,从流域内某一个特定的防洪区域来看,洪水造成的最终淹没范围可能有两种主要的形式,一种是由于降水引起流域水位上升或者积水超过堤防的高程,流入下游区域形成大面积淹没范围,即漫堤式淹没;另一种是洪水冲破堤防从决口处流入下游区域,形成淹没区,即决堤式淹没或堤防溃决[4].根据洪水发生的主要原因及洪水淹没的形态,可以把洪水淹没范围分析概化为以下两种不同的模型进行研究.
    1.1 给定洪水水位条件下的淹没分析
    给定洪水水位就是基于某一水位(H)下的淹没范围和水深的分布,适合漫堤式淹没[4].给定的洪水水位,可以是现状的水位,也可以是水文水力模型模拟预测的结果.此方法虽然没有考虑水动力学模拟水面形状变化等问题,但是这种给定水位近似“水平面”的方法,能便捷、快速地得到淹没面积和水深分布等核心数据,在防洪抗洪中有较大的实用价值[3].
    基于高分辨率DEM求给定水位条件下的洪水淹没,可以分为两种情况:一是把凡是高程值低于给定水位(H)的栅格点都计入淹没范围,没有考虑地域的连通性.这种情况就相当于整个区域大面积均匀降水,只考虑降水而引起的水位上升,没有地表径流的汇入;二是在前一种情形基础之上,考虑“连通”淹没的情况,洪水从源点漫堤后沿着地表流动,水流受到地表起伏特征的影响,即使某一低洼处高程低于给定水位,也可能由于地形的阻挡而不会淹没[1,3,6-7].这两种情况分别称为“无源淹没”和“有源淹没”[3].从计算机算法实现角度来看,无源淹没相对简单,只需找出低于给定水位的栅格即可.本文主要着眼于高效能计算,对有源淹没进行分析与计算.
    1.2 给定水量条件下的淹没分析
    在应对指挥调度等洪水灾情时,往往需要获得调度或者泄洪的水量(Q),或是决堤淹没发生后,洪水冲破堤坝从决口流入下游区域,可以获得大致的库容流出水量.基于该估算的洪水水量,根据分洪区地形特征计算洪水淹没的范围,及时为下游区域进行预报,从而进一步提高防洪的科学性.
    对给定洪水水量(Q)条件下的洪水淹没计算,是在给定洪水水位采用的有源淹没计算基础上进行,可考虑采用“体积法”计算[8].其主要思路是通过泄洪、溃坝、漫堤等流进下游区域的水量(Q)与淹没范围内总的水体体积相等的原理,来解算洪水淹没范围.
    2 基于DEM淹没分析思想与算法
    数字高程模型(DEM)是地形表面形态属性信息的数字表达,是带有空间位置特征和地形属性特征的数字描述[9].其中最基本规则格网(GRID)形式是由一系列地面点Z(x,y)及其所具有的高程Z组成,在计算机实现中则是一个二维数组,可以很容易地用计算机进行处理和运算[6,10].其结构模型一般用数字表达式:
    DEM={Z(x,y)}x=1,2,…,m;y=1,2,…,n     (1)
    2.1 给定洪水水位下的有源淹没分析思想与算法
    有源淹没是从淹没源点(一般指水库泄洪口、漫坝处等水源输出地)开始向低于水位线并且连通着的洼地淹没,所以它具有两个明显的特征,一是淹没点的高程值低于水位值,二是有低于水位值的连通路径连接到淹没点.根据规则格网(GRID)模型,每一个栅格Z(x,,y)四周有8个领域的栅格Z(i,j).对于给定的洪水水位(H),可以从淹没源点Z(x,y)开始充分利用多线程技术进行8邻域搜索的方法获得洪水淹没范围.其流程如图1.

图1 给定洪水水位(H)淹没范围计算主要流程
    1) 洪水水位高程的确定.特定淹没区的洪水水位高程的确定主要来自实测数据或者堤坝高程和溃坝处高程,也可能是汛期中一些水位高程记录点或者是大洪水发生时航测照片的数据,也可以是有水文水力学模型演算的数字[12].
    2) 建立可以存储栅格淹没区的二维数组和布尔型数组.建立一个ArryList型的二维数组是为了用来存储从DEM中读取的Z值低于水位高程(H)栅格点,同时建立一个Bool型的数组来存储该栅格点是否被访问过,如果访问过则设为“True”,否则为“False”.
    3) 获得DEM相关数据.主要通过ArcEngine开发组件接口IRasterCursor和IPixelBlock3来获得DEM的坐标及其Value值.通过IRasterCursor获得整个区域DEM的栅格游标,然后调用IPixelBlock3获得二维数组等有效信息.
    4) 循环8邻域栅格扩散.首先把淹没源点加入ArryList型的二维数组,并把相应Bool型的数组赋予True.然后从源点开始搜索其周围8邻域Z(i,j)栅格,把满足Z(i,j)≤H点则加入ArryList型的二维数组中,同时移除其中第一个栅格点,数组中第一个索引点则为邻域中的一个点,再次以第一个索引点为淹没源点进行8邻域搜索循环,直到ArryList型的二维数组为空则终止循环.Bool型的数组的标记就完整记录下洪水淹没范围,所有后续处理都可以应用这个标记去做各种相应的分析,获取感兴趣的数据.
    下面是实现该算法的关于C#语言核心伪码:
    Pre增加淹没源点到ArryList型的二维数组中.
    Post判断ArryList型的二维数组是否为空,不为空则循环执行下列代码.
    Return标记了是否淹没的Bool型的数组SearchFlood(i,j)获得淹没ArryList型的二维数组中索引号为0(即淹没源点)坐标把相应坐标的Bool型的数组赋予True从ArryList型的二维数组中移除索引号为0元素
     For(i=x-1;i≤x+1;i++)
         For(j=y-1;j≤y+1;j++) 
             If(i≠x&j≠y&Z(i,j)≤H)
    把该点加入ArryList型的二维数组中判断ArryList型的二维数组是否为空,不为空则循环SearchFlood(i,j)2.2给定洪水水量(Q)条件分析思想与算法如前所述,给定洪水水量(Q)条件下的洪水淹没计算是在有源淹没计算的基础采用“体积法”的原理来计算.在有源淹没条件下,规则格网(GRID)形式DEM每个栅格点Z(i,j)的大小是相等并且固定,它所代表的相应地表的面积也是确定的,只要通过栅格淹没点的个数就可以很容易完成洪水淹没面积(S)和体积(V)的计算.
                   (2)

    式中:S淹没为洪水淹没总面积;Si为每个栅格点代表的实际地表面积;n为计算到的有源淹没点数目.
    有源淹没体积计算实现方法为:根据计算的栅格的水深和相应的栅格实际面积即可求解相应的淹没体积:
       （3）
    式中:V淹没为淹没区的体积;Si为每个栅格点代表的实际地表面积;H为给定洪水水位高程;Zi为第i个栅格点高程值.
    在特定的洪水淹没计算区域,对于给定的洪水量(Q)淹没到下游区域所形成的淹没高度和范围,主要通过与给定水位高程条件下淹没范围内总的水体体积相等的原理来计算.从计算区域最低处(必须与源点具有连通性)开始不断改变给定洪水水位(H)的条件,利用公式(3)求出相应的体积V淹没与给定的洪水量(Q)比较,其表达式可以描述为:
   (4)
    由公式(4)看出,当f(H)→0时,在Q水量下的洪水淹没范围及高程就是由此时H值所决定的,利用该H值来模拟在Q洪水量下的淹没范围就得到很好地解算.
    1) 获得DEM中的Z值最小点.最低点要求能与源点连通,即有低于水位值的连通路径连接到淹没源点.调用ArcEngine组件库中DataSourcesRaster命名空间的IRasterBandCollection和IRasterBand两个接口获得DEM相关信息,利用IRasterStatistics接口下的Minimum属性即可,最后遍历ArryList型的二维数组中的元素,取得最小点的下标(坐标).
    2) 计算给定水量(Q)下的水位.首先把获得高程最小点之一作为有源淹没的源点,把最低高程值定义为最小Hmin,定义一个比那好步长△H,这里设为1M,循环执行计算有源淹没,并返回洪水体积V,每执行一次加△H,利用公式(4),直到当F(H)≥0或者F(H+△h)≤0的H值为止,H+△h/2即为Q洪水量下的近似淹没水位,利用这个数据进行洪水淹没分析及其虚拟现实表达三维可视化仿真,就能得到防洪决策最关心的核心数据.
    3 算法实现与虚拟现实表达
    虚拟现实技术使人们能够在三维虚拟环境中更好地对洪水淹没的现象和规律进行观察、控制和相关分析,因此它是开展洪水演进研究工作的最佳选择之一[1].根据以上算法思想与实际需要,进行了洪水淹没虚拟现实表达实现.该测试系统的开发语言为C#,在ESRI公司推出的ArcGISEngine产品基础上进行.
    ArcGISEngine是ArcGIS系列产品中的二次开发组件库,是一组跨平台的嵌入式ArcObjects,用来构建定制GIS和桌面制图应用程序,也可以开发出合乎用户需求的软件系统[13].
    洪水淹没现实表达主要是在VisualStudio.NET2008开发平台基础上利用ArcEngine9.3开发组件库进行代码实现,其中主要运用了Scene Control控件进行三维可视化仿真,调用I3D Properties与ILayerExtensions等接口进行了三维空间的建模,调用了IRaster Stretch Color RampR enderer和IAlgorithmic ColorRamp等接口进行了地形的渲染;充分利用了多线程技术和多核计算机的硬件优势,提高了程序的执行效率.实验数据DEM分辨率为5m×5m,区域地形主要通过786×767个像元表达出来,区域内最低高程为438m,最高为4361m.分别虚拟了给定洪水水位条件和给定洪水水量两种条件下的洪水淹没分析系统,可以比较快速、准确地计算出淹没范围和水深分布,并在三维空间里虚拟表达(图2,图3).


图2 给定洪水水位800m 和1000m 淹没范围模拟




图3给定洪水量50×108m3和200×108m3淹没范围模拟
    4 优势与特点
    1) 模拟了洪水淹没的最终水面.洪水淹没的机理是水源区和淹没区有连通性和存在水位差,洪水淹没最终的结果应该是水位达到平衡状态,这个时候的淹没区就应该是最终淹没区.基于水动力模型可以将洪水淹没过程模拟出来,对分析洪水淹没过程是非常有用,但对于最终的淹没范围结果与淹没概化分析模型没有多大的区别[4].同时以上两种方法能够在GIS显示平台上以最快的速度显示这些淹没结果.
    2) 淹没范围计算的准确性在很大程度上依赖于淹没区DEM的分辨率.给定洪水水位与水量概化的这两种分析模型,都是以DEM为数据基础,利用8邻域搜索的原理来获得淹没面积的.由于DEM的每一个像元代表了地面上一定的区域,且都拥有自己的高程值和精度.依据像元的高程值精度,选择恰当的阈值作为判断像元归宿的条件,该值称为种子蔓延探测分辨率.设空间数据精度为m,探测分辨率为f,经分析和计算实验,其关系式表达为[3,8]:
         z(5)

    式(5)中n为倍频系数,通常取2,k为平衡参数,用以平衡探测分辨率和运算效率间的矛盾.
    3) 快速、高效地将洪水淹没结果虚拟表达出来.不同的洪水淹没模型的计算结果都应该包括某个时刻洪水淹没的范围、淹没高度、水量等淹没结果参数.充分利用这些参数快速、直观地将淹没结果以可视化的方式展示出来.虽然这种近似“水平面”的洪水面没有考虑水面形状的变化问题,但是在防洪救灾工作中确有着较大的实际意义[4].又能在可视化平台叠加各种专题数据进行拓扑叠加和空间分析,为提出快速、准确、直观和有效的灾害预防方案打下快速决策的基础.
    5 结论
    实验表明,在高分辨率DEM基础之上结合相关地理信息技术、高效的计算机技术,以及三维可视化技术,可以快速、直观地将洪水淹没范围用三维虚拟环境方式模拟出来,为洪水灾害的客观预测分析和灾害评估提供了重要基础和依据.淹没范围的计算精度和准确性主要取决于空间数据(DEM)精度的优劣.值得说明的是,虽然计算淹没范围算法在一定程度上能够满足洪水淹没的实现需要,也存在很多问题需要去解决和完善,特别是在构建水动力模型方面去模拟洪水淹没过程和以及在处理大范围数据方面需要改进和添加算法设计,这也是以后继续研究学习的方向和重点.
    参考文献:
    [1]张尚弘,赵登峰,姜晓明.洪水淹没预报的虚拟仿真表达[J].清华大学学报,2010,50(6):861-864.
    [2]葛小平,许有鹏,张琪,等.GIS支持下的洪水淹没范围模拟[J].水科学进展,2002,13(4):456-460.
    [3]刘仁义,刘南.基于GIS复杂地形洪水淹没区计算方法[J].地理学报,2001,56(1):1-6.

    [4]丁志雄,李纪人,李琳.基于GIS格网模型的洪水淹没分析方法[J].水利学报,2006(6):56-59.
    [5]张成才,常静,孙喜梅,等.基于GIS的洪水淹没场景三维可视化研究[J].北京师范大学学报,2010,46(3):329-332.
    [6]刘坤,吕学强,陈丹,等.基于区域生长法的洪水淹没分析[J].北京机械工业学院学报,2009,24(1):1-3.
    [7]祝红英,顾华奇,桂新,等.基于ArcGIS的洪水淹没分析模拟及可视化[J].测绘通报,2009(5):66-68.
    [8]刘小生,黄玉生.“体积法”洪水淹没范围模拟计算[J].测绘通报,2004(12):47-51.
    [9]李志林,朱庆.数字高程模型[M].武汉:武汉测绘科技大学出版社,2000.
    [10]邬伦,张晶,赵伟,等.地理信息系统[M].北京:电子工业出版社,2002.
    [11]周品,李勇,谭建军,等.基于DEM的洪水淹没计算机算法优化研究[J].微计算机信息,2007,23(1-3):196-198.
    [12]甘郝新,邓抒豪,郑斌,等.基于GIS的洪水淹没范围计算[J].人民珠江,2007(6):98-100.
    [13]AL-SABHANW,MULLIGANM,BLACKBURNGA.AReal-TimeHydrologica lModel for FloodPrediction Using GISandtheWWW[J].Computers,Environmentand Urban Systems,2003,27(1):9.