基于DEM 和“填充法”的库容曲线计算方法

安航永，郭亚伟，张 浩

（1.广州禺山水务勘测设计股份有限公司,广东 广州 511400；2.江苏省灌溉总渠管理处,江苏 淮安 223200；3.濮阳黄河河务局范县黄河河务局,河南 濮阳 457506）

0 引言

库容曲线是水库规划设计和管理调度的重要依据。库容曲线的不准确将造成调洪演算中的水量不平衡,进而影响水库的发电、防洪、灌溉等效益。常用库容曲线计算方法有三类：图上量测方法、基于DEM（数字高程模型）和GIS 技术的计算方法、基于遥感数据的计算方法。图上量测的方法主要有“格网法”“等高线法”“横断面法”等[1]。“格网法”“等高线法”和“横断面法”通常在精度比较低的纸质地形图或电子地图上量算,计算工作量大,耗时较长,计算结果精度较低。基于DEM 和GIS 技术的计算方法以航测获得库区DEM为基础,利用ArcGIS 进行分析计算,计算精度较高[2]。但是,该方法要求工作人员具有较高的ArcGIS 技术功底和一定的专业知识,且当计算水位间隔较小时,人工重复操作较多,工作量明显增加,耗时较长[3]。基于遥感数据的计算方法以库区遥感影像数据为基础,提取水域面积,然后查寻影像拍摄时的水位,进而推求库容曲线。但是,该方法仅适用于水库建成后的情况,且由于无法捕捉死水位以下及高水位（设计洪水或校核洪水）时的遥感影像,该方法对死水位以下库容及高水位相应库容计算误差较大[4]。

本文基于DEM 和“填充法”理念提出库容曲线计算方法,并编写MATLAB 计算程序,一次性得出库容曲线全部数值。

1 基本原理

本文库容曲线计算的基本原理来源于“填充法”的理念：将库容计算视为一个盛水的容器,水库库容即为填充水量,见图1。

图1 “填充法”库容计算基本原理

库容曲线计算过程：

（1）将计算范围内现状地形细分成多个边长相等的正方形计算单元,认为各单元表面高程相同,相邻单元表面高程可以不同。

（2）根据实测地形数据,采用MATLAB 三维插值函数插值得各计算单元表面高程。

（3）每个计算单元填充水体体积为计算水位与地面高程之间的立方体体积,即V=dh×s,dh 为计算水位线与计算单元地面高程的差值,s 为计算单元表面积。

（4）将计算范围内V 值小于0 的单元体积赋值为0（即计算单元地面高程高于计算水位线）。

（5）累加所有计算单元水体体积即为对应于计算水位的水库容积。

（6）重复（1）～（5）计算不同水位所对应的水库容积,即可绘制水位～库容曲线。

水库容积计算公式：

式中：V 为计算水位对应的水库容积；h 为计算水位；hi为计算范围内某一单元地面高程；s 为计算单元面积；n 为计算范围内计算单元个数。

2 程序框架与设计

2.1 程序框架设计

本程序框架主要分为3 部分：首先是数据预处理阶段,该阶段主要实现数据的预处理和导入,其中数据预处理通过CAD 完成。其次是参数设置阶段,该阶段主要设置计算单元边长及高程计算步长。最后是计算及结果输出阶段,该阶段根据导入数据及计算参数,完成库容曲线计算工作,并将计算结果导出到Result.xls 中。

图2 库容曲线计算程序框架

2.2 程序详细设计

本程序以实测地形数据（txt 文件）为计算起点,以生成水位～库容曲线为最终目标,程序功能具有严密的流程,计算数据具有明显的时序性。因此,本程序采用流程式架构,即根据计算数据运算的先后顺序划分程序模块。程序总体上分为4 个模块：

（1）数据导入模块：读取经过预处理的txt 文件,采用load（）函数将高程点数据赋值给data 数组,将计算范围边界坐标数据赋值给bj 数组,水库实测地形见图3,预处理后高程点及计算范围边界示意图见图4。

图3 水库实测地形图

图4 预处理后高程点及计算范围示意图

（2）参数设置模块：分别设置计算单元边长d 及高程计算步长b。

张清元干活很卖力，院里种有几块菜地，挖田、起沟、挑粪的活都是他干，他整的田垅平平整整，土块也打得细。黎院长常常奖励他，比如给他炕一个红苕，或是给一条黄瓜吃。张清元在院里过得也算开心。

（3）数据计算模块：①根据计算范围边界坐标数据及计算单元边长,采用meshgrid（）函数对计算范围进行网格划分,见图5；②采用griddata（）函数进行曲面差值,生成DEM 数组,见图6；③根据预设数学算法计算水库容积。

图5 计算范围网格划分示意图

图6 插值生成DEM 示意图

（4）数据导出模块：计算完成后,采用xlswrite（）函数将计算成果导出到Result.xls 文件中。

2.3 计算程序

本程序库容曲线计算代码如下：

2.4 使用方法

（1）数据预处理

在实际工作中,地形资料通常以CAD 文件格式给出。为了方便程序计算,需要对原始测量数据进行预处理。

高程点数据提取：利用CAD 软件中“数据提取”命令,提取实测地形文件高程点的坐标及高程信息,并保存为txt 文件。txt 文件存储格式为“X Y Z”（其中,X 列为横坐标值,Y列为纵坐标值,Z 列为高程值）,数据中间以空格断开。

计算范围数据提取：首先,根据实测地形及坝顶高程等信息,利用CAD 软件中“多段线”命令框定计算范围；然后,选中计算范围线,利用“list”命令提取计算范围线折点坐标值,并保存为txt 文件。txt 文件存储格式为“X Y”（其中,X列为横坐标值,Y 列为纵坐标值）,数据中间以空格断开。由于本程序只计算计算范围内水库容积,为保证计算精度和计算速度,应根据实际情况综合确定计算范围。

（2）操作方法

将高程点数据文件（XYZ.txt）及计算范围线数据文件（BJ.txt）存放在本程序.m 文件所在文件夹内。同时,在该文件夹内创建一个名称为“Result”的excel 文件,并将该文件中的一个工作表命名为“Result”,用来存放计算结果。然后,在MATLAB 平台中打开本程序所在.m 文件,点击运行即可开始计算。

（3）结果导出

本程序计算完成后将通过“xlswrite”函数把计算结果写入Result.xls 文件内Result 工作表中。计算结果存储格式为“A B”（其中,A 列为水位数值,单位为m；B 列为库容数值,单位为万m3）,水位间隔即为上述设定高程计算步长b。

3 应用实例

本文以七盏灯水库为例,采用上述程序计算七盏灯水库库容曲线。七盏灯水库地处丘陵地带,位于广州市番禺区沙头街西北部大夫山森林公园内,库长约400 m,库水面宽130 m～155 m,汇水区内最高点为大坝南西方向的大乌岗,高程226 m左右；坝址区为丘陵前的冲积洼地,两岸低丘高程多在19.00 m～54.00 m,左岸低丘最高点高程110.00 m,右岸最高点高程88.00 m,地形坡度5°～38°；大坝下游河谷出口处地形较平缓,地面高程15.00 m～18.00 m,现为西大夫山森林公园。

七盏灯水库于1958 年建成,集雨面积0.72 km2,为小（2）型水库,工程等别为Ⅴ等,主坝级别为5 级。水库防洪标准为20 年一遇,校核标准为50 年一遇。水库大坝现状为斜墙土坝,坝长180 m,坝高8.40 m,坝顶高程25.00 m。大坝溢洪道为单孔涵式,进口段为5.00 m（宽）×1.95 m（高）,总长度103.20 m。大坝放水洞位于溢洪道下方,进口暗涵内尺寸1.50 m×1.50 m。设有一扇1.70 m×1.70 m 平板钢闸门,配备手动螺杆启闭机。

七盏灯水库曲线计算采用1∶500 实测地形。首先,根据实测地形及坝顶高程框定计算范围。然后,利用CAD 软件中“list”命令,提取计算范围线坐标,将坐标数据保存到BJ.txt文件中。接着,利用CAD 软件中“数据提取”命令,提取高程点坐标及高程数据,将高程点数据保存到XYZ.txt 文件中。在地形数据预处理完成后,将BJ.txt、XYZ.txt 文件复制到本程序.m 文件所在文件夹。在MATLAB 平台中打开本程序所在.m 文件,设置计算单元边长d=1,高程计算步长b=0.1,并点击运行。

软件运行结束后,计算结果自动存入Result.xls 文件内Result 工作表中,计算得七盏灯水库库容曲线见图7、表1。

图7 七盏灯水库库容曲线

表1 七盏灯水库库容曲线

4 结语

本文以七盏灯水库为例,采用基于DEM 和“填充法”理念的计算程序,方便、快速地完成了七盏灯水库库容曲线的计算工作。该方法不但操作简单、计算快速,而且计算精度也较高,可作为水利工程技术人员完成库容曲线计算工作的一个手段。