基于ArcPy技术的特大型水库库容量算研究<br/>——以丹江口水库为例

基于ArcPy技术的特大型水库库容量算研究
——以丹江口水库为例

2017-07-24夏定华孙慧敏

水利水电快报 2017年6期

关键词：丹江口水库库容分块

杨松夏定华孙慧敏刘沛黄军

(长江水利委员会水文局汉江水文水资源勘测局，湖北襄阳 441022)

基于ArcPy技术的特大型水库库容量算研究
——以丹江口水库为例

杨松夏定华孙慧敏刘沛黄军

(长江水利委员会水文局汉江水文水资源勘测局，湖北襄阳 441022)

针对特大型水库特点进行精确、快速的库容量算对编制水库调度方案具有至关重要的作用。介绍了常用的DEM法库容量算原理和方法，并对GIS平台下的ArcPy技术展开了研究。研究得出，采用ArcPy技术可快速实现特大型丹江口水库的库容分层分块复核量算，提高了特大型水库库容量算效率。

库容量算；DEM；ArcPy；特大型水库；丹江口水库

目前，水库库容量算方法主要包括等高线面积法、断面法、测点构网法、数字高程模型(DEM)法等。等高线面积法是假定两等高线之间体积变化呈线性，无法精确反映两高程之间的地形起伏，特别对地势较为平坦地区，计算结果误差较大；断面法主要是在库区设定相互平行的断面，要求断面间的地形变化较为均匀，呈线性，而对于面积较大、支汊较多的水库，计算精度难以保证；测点构网法和数字高程模型(DEM)法，是目前较为普遍的方法，但是对于特大型不规则水库的库容量算，采用现有的方法则计算效率较低。因此，针对丹江口特大型水库支汊多、面积大的特点，研究基于ArcPy技术实现不规则特大型水库库容量算，计算效率得到了提高，库容量算成果得到了复核评价，满足水库调度运行需要。

1 DEM法库容量算原理和方法

DEM法计算水域库容的基本思路是用不规则的三角形面(TIN)或者栅格单元面(栅格DEM)尽可能地与实际、真实的水下地形曲面紧密吻合，建立水下地形曲面模型，水下地形曲面模型和水域某一高程平面相交则形成水域库容的实体模型[1]。根据数据结构类型DEM法主要分为测点构网法和栅格DEM法。

1.1 测点构网法

根据实测的水道地形图建立TIN，累积计算TIN在每个小区域上的容积，即为水库的总容积。

如图1中，设空间三角形的3个顶点为A、B、C，顶点三维坐标为(xa，ya，za)、(xb，yb，zb)、(xc，yc，zc)，且za≥zb≥zc，可以通过排序得到这样的假设。

图1 三角形区域上的容积计算示意

假设计算高程面为z，三角形A′B′C内角A′、C的正弦值分别为sinA′、sinC，A′B′、B′C、CA′边长分别为c、a、b，则容积(v)、接触表面积(r)的计算公式推算如下：

如果z≤zc，则a=0，v=0；

如果zc≤z≤zb，则

(1)

(2)

如果zb≤z≤za，则

(3)

(4)

如果z≥za，则

(5)

(6)

图2 丹江口水库库区示意

如果数字高程模型为TIN(不规则三角网)模型，则直接采用上式计算各三角形区域上延的体积，然后累加即为容积。

1.2 栅格DEM法

水库容积定义为某一水位下或两水位之间的蓄水容积。DEM法采用数字高程模型将库区水体微分为若干个底面为正方形的柱状体，通过正方体体积积分得到水库库容，其数学模型为：

(7)

式中，V为库容，Ps为单个栅格的面积，H为计算水位值，hi为处于计算水位H以下的栅格单元的高程值，n为处于计算水位H以下的栅格单元个数。

2 ArcPy技术在丹江口水库中的应用

2.1 库区概况

丹江口水库为丹(丹江)、汉(汉江)两江并联式特大型水库，枢纽大坝位于湖北省丹江口市城区，在丹江入汇汉江的口门下游0.8 km处。丹江口水库汉江库区(简称汉库)自上而下行政隶属分别为陕西省白河县、湖北省郧西县、郧县、丹江口市；丹江库区(简称丹库)自上而下行政隶属分别为河南省淅川县、湖北省丹江口市。

丹江口水利枢纽于1958年9月开工兴建，1968年蓄水运用，初期工程于1973年竣工。大坝坝顶高程为162 m，当正常蓄水位157 m时，汉江库区回水长177.4 km、丹江库区回水长85.4 km，形成了汉库、丹库两个容积与面积基本相当的并联式水库(见图2)。丹江口大坝加高工程完成后，水库大坝坝顶高程为176.6 m，正常蓄水位调整为170 m。

丹江口水库库区四周多为山地，沿程宽、窄相间。正常蓄水位170 m时，汉库内最宽处约5 km、最窄处约400 m；丹库内最宽处(俗称“小太平洋”)约17 km、最窄处约250 m。在水库变动回水区，河床以卵石、砾石、中粗沙为主；在水库常年回水区以中细沙、淤泥絮泥为主。库区气候年内温差较大，盛夏最高气温超过40℃，冬天最低气温达-10℃。主汛期一般在7～9月份，水量较丰沛。

2008年丹江口水库库容量算时考虑整个丹江口库区地形数据量庞大，当时计算采用的程序软件计算数据量有限，因此，库容量算方法采用了分块计算，将整个丹江口水库共划分为25块分别进行计算。2015年对库容量算结果进行复核计算仍然按照2008年量算分块的相同范围进行，并对如何提高计算效率做了进一步研究[2]。

2.2 ArcPy简介[3]

空间数据的处理常常复杂而费时，而且其过程具有很强的重复性，基于窗口的应用程序极大程度上限制了处理数据能力，有必要进行自动化处理。任何可以支持COM的脚本语言都可执行ArcGIS的地理处理工具，如Python、Jscript和VBScript等。这些脚本语言都是公开的，而且易学易用，脚本可以将大批量的处理通过流程化的方式，动态地传入、传出参数，达到批量、自动处理数据的目的。

Python是一种面向对象的、跨平台的开源脚本编程语言，其功能非常强大而且较为简单，易于学习[4-5]。ESRI公司自ArcGIS 9.0版中开始引入此语言，并已成为其首选的脚本语言，能够用于数据分析、数据转换、数据管理和地图自动化等，可以很好地根据用户需求定制专门的工具，提高工作效率。

ESRI公司将ArcGIS中的地理处理工具以及一些扩展模块封装成了ArcPy站点包，其中还有多种可用于处理和询问GIS数据的函数和类。用户可以借助Python和ArcPy调用工具，开发出各种实用程序处理地理数据。

利用ArcPy开展地理处理批处理脚本程序时，只需要通过引入地理处理模块，创建地理处理对象，就可以开展程序开发。主要步骤分为两步：

(1) 第1步，导入地理处理模块：

import arcpy

(2)第2步，调用地理处理模块的属性和方法，以缓冲区分析函数为例：

arcpy.env.workspace = "C:/data"

arcpy.Buffer_analysis("roads", "C:/output/majorrdsBuffered", "100 Feet", "FULL", "ROUND", "LIST", "Distance")

2.3 ArcPy实现批量库容量算批处理

(1) 根据ArcGIS处理单个计算区域库容量算原理和方法，可以根据TIN法计算库容，基于同样的原理和方法，结合ArcPy实现多个计算区域库容量算的批处理，见图3。

图3 Arcpy批量库容量算处理流程

(2) 数据文件组织方式，以丹江口水库为例，可将25个分块的TIN文件放于同一文件夹路径，根据处理流程，给定计算水位值的最小高程值和最大高程值，循环遍历所有文件，并设置计算水位值间隔，以此计算所有TIN文件的库容曲线。以丹江口水库库容分块量算为例，数据组织见图4。

图4 ArcPy计算库容的数据组织方式(TIN数据)

(3) 创建大水体库容分块量算的脚本文件(见图5)。

#输入ArcPy站点包

import arcpy

from arcpy import env

arcpy.CheckOutExtension("3D")

files =arcpy.ListFiles()

for fc in files:

for i in range(n,m):

#循环计算各水位值，间隔为1m

dl=fc+str(i)+".txt"

arcpy.SurfaceVolume_3d(fc,dl,"BELOW",mgs,1,0)

图5 新建GIS脚本工具界面

图6 大水体库容分块量算GIS插件界面

(4) 创建大水体库容分块量算GIS工具流程(见图6)。首先，在Arccatalog中选择New->Toolbox，并设置文件名；其次，进入新建的Toolbox，空白处右键Add->Script…，并按照向导完成(属性可以留空，以后使用工具属性进行设置)；最后，导入计算功能的.py文件，并设置输入参数，输入参数类型与工具的参数类型一定要一致。