基于RS和GIS的大南川水库溃坝灾害模拟方法研究

韩丽蓉，张 涛

(1.青海大学地质系，青海西宁810016;2.青藏高原北缘新生代资源环境重点实验室，青海西宁810016;3.浙江省电力建设有限公司，浙江宁波315000)

一、引 言

我国已经建设的一部分水库大坝存在着安全系数低的问题，大坝一旦溃决，势必将对下游造成毁灭性的灾难［1-2］。GIS技术运用于洪水灾害的评估始于20世纪90年代，但是很多研究和技术方法都比较粗略，很少考虑利用其空间查询分析功能提高计算分析的速度与准确性。刘仁义、刘南提出了有源淹没和无源淹没两种模型;刘小生、黄玉生提出了利用“体积法”计算洪水淹没面积［3］;周品、李勇等利用DEM数据结合“体积法”的概念公式，得到淹没面积［4］;常燕卿、张福浩提出将洪水在特定河段内的DEM简化为一斜平面，把淹没范围归结为研究区域原始DEM被一个斜平面切割［5］。洪水演进模型主要是基于种子蔓延算法［1-7］，其不足之处在于必须编写程序比较复杂。

二、洪水演进模型

本文探讨了一种无需编程便可快速计算水库溃坝后下游淹没区的方法，适用于溃坝源头与淹没区落差较大地区。该方法结合本流域数字高程模型(DEM)，采用分段平面模拟方法。分段平面模拟法的思想是以本流域沿途乡镇作若干过水点特征横断面，将本流域划分为若干个淹没小段，把洪水在相邻横断面内的DEM简化为一个斜平面，并根据水动力演进模型模拟得到各断面淹没水深;然后，根据DEM解析出的山谷线高程和横断面的淹没水深，计算得到河道两端横断面的洪水高程;最后，通过将洪水淹没面DEM与研究区域原始DEM叠加，并进行填挖方处理，得到洪水淹没区情况。其中，填方区域表示洪水淹没区，填方区与挖方区交线即为淹没区边界，将填挖方DEM转换成面要素后删除挖方部分，保留的填方部分为淹没区范围。

具体处理时对水库及下游作如下简化:设水流的滑动摩擦因数为μ;水库坝高为x0，坝长为y0;水库距下游目标研究淹没区平距为D，高差为h。如图1所示。

图1 淹没演进模型

假设水库大坝由于地震等影响遭受剧烈破坏，造成1/2溃决，则可计算得出溃口面积大小

由于坝址处被已溃坝体阻挡，水流不会瞬间落至1/2坝体处，假设水流瞬间从坝顶作自由落体至3/4坝高，水库内水的原始速度为0，根据势能转换关系得到3/4坝高处的水流速度

由以上两式可得水库溃坝后库容放空时间

式中，K为水库的库容量;L为洪水溃口处流量;M0为洪水溃口面积。

水流出库后将沿着地表顺流而下，地表近似为一倾斜面，则水流在地表流动的加速度为

式中，θ为斜面坡度;μ为地面对水流的滑动摩擦因数。

设水流到达淹没区的时间为t，则根据运动方程有

经过计算得到水流到达淹没区特征过水点横断面的时刻t与流速vt

考虑水库溃坝后洪峰经过时沿途淹没高度，设在t时刻洪水淹没区的横截面面积为Mt=xtyt;洪水流量恒定不变，为Lt=Mtvt;溃口处流量为L0=Lt;即有xtytvt=M0v0，从而得到t时刻洪峰到达的淹没区的淹没高度

三、大南川水库溃坝灾害分析

1.建立研究区域DEM

大南川水库建于1974年，位于青海省湟中县，集雨区面积407 km2，水库库容1300万m3，坝顶高程 2 640.5 m，坝长460 m，最大坝高 46.5 m，属于湟水支流南川河水系。利用研究区域的地形图，在ArcMap中线性内插出高程数据，建立三角网TIN;然后把TIN格式转换成GRID格式DEM。

2.溃坝模拟

(1)提取研究区域山谷线

从DEM提取山谷线，山谷线与断面交点为该断面最低点，进而根据山谷线，并以洪水高度为依据确定淹没水位。首先利用ArcMap的栅格计算器得到与淹没前DEM地形相反的反地形;再按山脊线的提取方法得到山谷线［8］。

(2)提取特征断面

本文选择大南川水库流域的8个人口密集区为代表，即大南川河水库上游1个及下游7个横断面，绘制横断面图。断面方向与南川河法线(山谷线的法线)基本重合，考虑到大南川水库的有限库容量及洪水流经的极限流量，经反复试验，提取的断面最高点与最低点的高差基本保持在45 m，保证了在洪水深度的计算中能获得较高图解精度的断面长度。

本文选择绘制7个分段纵断面图。(3)计算洪水到达时间与洪水水深由大坝数据及洪水演进模型分析可知，溃口大小为

此时水流速度为

水流流量大小为

水库放空时间为

假设洪水的过水横断面近似为矩形，地面对水流的滑动摩擦因数μ=0.015，重力加速度g=10 m/s2，利用各乡镇横断面与纵断面，并结合淹没模型，可得到各横断面长度yt、各纵断面坡度及坡角、洪水到达各横断面的时间、洪水到达时的流速及洪水深度等。

(4)计算淹没后洪水DEM

根据以上计算得出的洪水深度，在ArcMap中模拟绘制出洪水淹没面DEM。由于两横断面间距较大，因此采用插值方法加绘洪水淹没面等高线，等高线的走向与山谷线的法线基本重合;然后利用ArcMap的3D分析工具获取等高线与山谷线交点处的高程属性，并将交点处的山谷线高程值加上洪水深度，即得到该等高线所在位置的洪水面高程;最后把得到的洪水面高程赋予绘制的等高线，生成TIN格式的DEM，并将其转换为分辨率为5 m的GRID格式的DEM。因为研究区域落差较大，最高点与最低点之差达400 m，因此如果以某特定高程面为洪水淹没面会带来很大的误差，甚至不能得出结果，导致模拟失败。

(5)确定淹没区范围

将上面得到的洪水DEM与研究区域原始DEM叠加，并进行填挖方处理，得到填挖方范围。其中填方区域表示洪水淹没区，填方区与挖方区交线即为淹没区边界。挖方区在本研究中无意义，应将其删除。但由于在ArcGIS中不易对DEM进行裁剪处理，因此需要先将填挖方DEM转换成面要素，然后删除挖方部分，保留填方部分，即得到洪水淹没区，如图2所示。

图2 洪水淹没区范围

(6)淹没区域三维可视化

在ArcScene中以前面讨论得到的研究区域DEM的值为高程值，在其表面叠加显示遥感影像和淹没区范围［9］。由于ArcScene中无法对行政区地名进行标注，因此利用ArcScene的3D编辑器制作三维地名及洪水到达时间的三维符号，得到淹没区的三维景观，如图3所示。

图3 淹没区域三维场景

3.淹没区损失分析

(1)淹没区遥感影像处理

首先，利用ArcGIS的数据格式转换功能将面要素格式的淹没区范围转换成栅格格式(像素大小为5 m)，再转换为Erdas能用的有地理参考的IMG格式(其存在背景噪声);其次，在ERDAS视窗中对淹没范围进行边缘检测，去除背景噪声，不关闭视窗，在Vector菜单中选择相应功能将淹没范围转换成矢量格式，再转换为栅格格式;再次，在ERDAS中以淹没区域为裁剪区域，对研究区域遥感影像进行掩模裁剪，得到淹没区范围的遥感影像;最后，对淹没区遥感影像进行非监督分类［1］，将淹没区地物分成5类:房屋、耕地、林地、水面、空闲地，进而得到各类的面积。

(2)损失情况分析统计

淹没区域主要覆盖西宁市城南新区和城中区，面积合计151 000 000 m2，人口合计29.7万。本次洪水淹没面积为30 726 900 m2，估计受灾人口为30 726 900÷151 000 000×29.7=6.0万人。

四、总结与展望

本文综合运用ERDSA和ArcGIS等手段，以西宁市湟中县的大南川水库为例，讨论了洪水演进模型，提出了一种快速简便的适用于落差较大地区洪水淹没分析的洪水演进模型。通过该模型可以计算出下游乡镇遭遇洪水的时间、洪水淹没的深度、洪水的流速等重要洪水参数，进而能确定洪水淹没范围。模拟演示了流域的三维场景，并对淹没范围的遥感影像进行非监督分类，估计得出溃坝洪水给下游造成的损失，效果比较理想。

本研究中还存在不少需要改进之处:①所建立的洪水演进模型只考虑了洪峰演进情况，没有考虑库容大小、上游来水和集雨区降水量，会给研究结果造成一定的误差;②所使用的图件成图时间较早，分辨率偏低，影响研究结果的准确性。洪水灾害的预测是一个非常复杂的课题，其受到很多因素的影响，还有待进一步研究。

［1］王韶玉.基于DEM的坝堤溃决洪水淹没评价模型与方法研究［D］.武汉:华中科技大学，2010:56-58.

［2］周远方.大南川水库溃坝的数值模拟研究［D］.长沙:长沙理工大学，2010:63-64.

［3］刘小生，黄玉生.基于Arc/Info的洪水淹没面积的计算方法［J］.测绘通报，2003(6):46-48.

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

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

［6］何宗宜，黄春雄，韩用顺.水库下游洪水淹没灾害评估系统的研究［J］.测绘通报，2003(2):24-27.

［7］冯平，纪恩福，卢永兰.水库下游洪灾淹没损失的估算［J］.灾害学，1995，10(1):5-9.

［8］汤国安，杨昕.地理信息系统空间分析实验教程［M］.北京:科学出版社，2006:85-86.

［9］张成才，刘丹丹，于欣，等.基于ArcGIS的东平湖洪水淹没场景三维可视化［J］.郑州大学学报:工学版，2008，29(1):88-90.

［10］党安荣，贾海峰，陈晓峰，等.遥感图像处理教程［M］.北京:清华大学出版社，2010:187-192.