安航永，郭亞偉，張 浩

（1.廣州禺山水務(wù)勘測(cè)設(shè)計(jì)股份有限公司,廣東 廣州 511400；2.江蘇省灌溉總渠管理處,江蘇 淮安 223200；3.濮陽黃河河務(wù)局范縣黃河河務(wù)局,河南 濮陽 457506）

0 引言

庫容曲線是水庫規(guī)劃設(shè)計(jì)和管理調(diào)度的重要依據(jù)。庫容曲線的不準(zhǔn)確將造成調(diào)洪演算中的水量不平衡,進(jìn)而影響水庫的發(fā)電、防洪、灌溉等效益。常用庫容曲線計(jì)算方法有三類：圖上量測(cè)方法、基于DEM（數(shù)字高程模型）和GIS 技術(shù)的計(jì)算方法、基于遙感數(shù)據(jù)的計(jì)算方法。圖上量測(cè)的方法主要有“格網(wǎng)法”“等高線法”“橫斷面法”等[1]?！案窬W(wǎng)法”“等高線法”和“橫斷面法”通常在精度比較低的紙質(zhì)地形圖或電子地圖上量算,計(jì)算工作量大,耗時(shí)較長,計(jì)算結(jié)果精度較低?；贒EM 和GIS 技術(shù)的計(jì)算方法以航測(cè)獲得庫區(qū)DEM為基礎(chǔ),利用ArcGIS 進(jìn)行分析計(jì)算,計(jì)算精度較高[2]。但是,該方法要求工作人員具有較高的ArcGIS 技術(shù)功底和一定的專業(yè)知識(shí),且當(dāng)計(jì)算水位間隔較小時(shí),人工重復(fù)操作較多,工作量明顯增加,耗時(shí)較長[3]?；谶b感數(shù)據(jù)的計(jì)算方法以庫區(qū)遙感影像數(shù)據(jù)為基礎(chǔ),提取水域面積,然后查尋影像拍攝時(shí)的水位,進(jìn)而推求庫容曲線。但是,該方法僅適用于水庫建成后的情況,且由于無法捕捉死水位以下及高水位（設(shè)計(jì)洪水或校核洪水）時(shí)的遙感影像,該方法對(duì)死水位以下庫容及高水位相應(yīng)庫容計(jì)算誤差較大[4]。

本文基于DEM 和“填充法”理念提出庫容曲線計(jì)算方法,并編寫MATLAB 計(jì)算程序,一次性得出庫容曲線全部數(shù)值。

1 基本原理

本文庫容曲線計(jì)算的基本原理來源于“填充法”的理念：將庫容計(jì)算視為一個(gè)盛水的容器,水庫庫容即為填充水量,見圖1。

圖1 “填充法”庫容計(jì)算基本原理

庫容曲線計(jì)算過程：

（1）將計(jì)算范圍內(nèi)現(xiàn)狀地形細(xì)分成多個(gè)邊長相等的正方形計(jì)算單元,認(rèn)為各單元表面高程相同,相鄰單元表面高程可以不同。

（2）根據(jù)實(shí)測(cè)地形數(shù)據(jù),采用MATLAB 三維插值函數(shù)插值得各計(jì)算單元表面高程。

（3）每個(gè)計(jì)算單元填充水體體積為計(jì)算水位與地面高程之間的立方體體積,即V=dh×s,dh 為計(jì)算水位線與計(jì)算單元地面高程的差值,s 為計(jì)算單元表面積。

（4）將計(jì)算范圍內(nèi)V 值小于0 的單元體積賦值為0（即計(jì)算單元地面高程高于計(jì)算水位線）。

（5）累加所有計(jì)算單元水體體積即為對(duì)應(yīng)于計(jì)算水位的水庫容積。

（6）重復(fù)（1）～（5）計(jì)算不同水位所對(duì)應(yīng)的水庫容積,即可繪制水位～庫容曲線。

水庫容積計(jì)算公式：

式中：V 為計(jì)算水位對(duì)應(yīng)的水庫容積；h 為計(jì)算水位；hi為計(jì)算范圍內(nèi)某一單元地面高程；s 為計(jì)算單元面積；n 為計(jì)算范圍內(nèi)計(jì)算單元個(gè)數(shù)。

2 程序框架與設(shè)計(jì)

2.1 程序框架設(shè)計(jì)

本程序框架主要分為3 部分：首先是數(shù)據(jù)預(yù)處理階段,該階段主要實(shí)現(xiàn)數(shù)據(jù)的預(yù)處理和導(dǎo)入,其中數(shù)據(jù)預(yù)處理通過CAD 完成。其次是參數(shù)設(shè)置階段,該階段主要設(shè)置計(jì)算單元邊長及高程計(jì)算步長。最后是計(jì)算及結(jié)果輸出階段,該階段根據(jù)導(dǎo)入數(shù)據(jù)及計(jì)算參數(shù),完成庫容曲線計(jì)算工作,并將計(jì)算結(jié)果導(dǎo)出到Result.xls 中。

圖2 庫容曲線計(jì)算程序框架

2.2 程序詳細(xì)設(shè)計(jì)

本程序以實(shí)測(cè)地形數(shù)據(jù)（txt 文件）為計(jì)算起點(diǎn),以生成水位～庫容曲線為最終目標(biāo),程序功能具有嚴(yán)密的流程,計(jì)算數(shù)據(jù)具有明顯的時(shí)序性。因此,本程序采用流程式架構(gòu),即根據(jù)計(jì)算數(shù)據(jù)運(yùn)算的先后順序劃分程序模塊。程序總體上分為4 個(gè)模塊：

（1）數(shù)據(jù)導(dǎo)入模塊：讀取經(jīng)過預(yù)處理的txt 文件,采用load（）函數(shù)將高程點(diǎn)數(shù)據(jù)賦值給data 數(shù)組,將計(jì)算范圍邊界坐標(biāo)數(shù)據(jù)賦值給bj 數(shù)組,水庫實(shí)測(cè)地形見圖3,預(yù)處理后高程點(diǎn)及計(jì)算范圍邊界示意圖見圖4。

圖3 水庫實(shí)測(cè)地形圖

圖4 預(yù)處理后高程點(diǎn)及計(jì)算范圍示意圖

（2）參數(shù)設(shè)置模塊：分別設(shè)置計(jì)算單元邊長d 及高程計(jì)算步長b。

張清元干活很賣力，院里種有幾塊菜地，挖田、起溝、挑糞的活都是他干，他整的田垅平平整整，土塊也打得細(xì)。黎院長常常獎(jiǎng)勵(lì)他，比如給他炕一個(gè)紅苕，或是給一條黃瓜吃。張清元在院里過得也算開心。

（3）數(shù)據(jù)計(jì)算模塊：①根據(jù)計(jì)算范圍邊界坐標(biāo)數(shù)據(jù)及計(jì)算單元邊長,采用meshgrid（）函數(shù)對(duì)計(jì)算范圍進(jìn)行網(wǎng)格劃分,見圖5；②采用griddata（）函數(shù)進(jìn)行曲面差值,生成DEM 數(shù)組,見圖6；③根據(jù)預(yù)設(shè)數(shù)學(xué)算法計(jì)算水庫容積。

圖5 計(jì)算范圍網(wǎng)格劃分示意圖

圖6 插值生成DEM 示意圖

（4）數(shù)據(jù)導(dǎo)出模塊：計(jì)算完成后,采用xlswrite（）函數(shù)將計(jì)算成果導(dǎo)出到Result.xls 文件中。

2.3 計(jì)算程序

本程序庫容曲線計(jì)算代碼如下：

2.4 使用方法

（1）數(shù)據(jù)預(yù)處理

在實(shí)際工作中,地形資料通常以CAD 文件格式給出。為了方便程序計(jì)算,需要對(duì)原始測(cè)量數(shù)據(jù)進(jìn)行預(yù)處理。

高程點(diǎn)數(shù)據(jù)提?。豪肅AD 軟件中“數(shù)據(jù)提取”命令,提取實(shí)測(cè)地形文件高程點(diǎn)的坐標(biāo)及高程信息,并保存為txt 文件。txt 文件存儲(chǔ)格式為“X Y Z”（其中,X 列為橫坐標(biāo)值,Y列為縱坐標(biāo)值,Z 列為高程值）,數(shù)據(jù)中間以空格斷開。

計(jì)算范圍數(shù)據(jù)提?。菏紫?根據(jù)實(shí)測(cè)地形及壩頂高程等信息,利用CAD 軟件中“多段線”命令框定計(jì)算范圍；然后,選中計(jì)算范圍線,利用“l(fā)ist”命令提取計(jì)算范圍線折點(diǎn)坐標(biāo)值,并保存為txt 文件。txt 文件存儲(chǔ)格式為“X Y”（其中,X列為橫坐標(biāo)值,Y 列為縱坐標(biāo)值）,數(shù)據(jù)中間以空格斷開。由于本程序只計(jì)算計(jì)算范圍內(nèi)水庫容積,為保證計(jì)算精度和計(jì)算速度,應(yīng)根據(jù)實(shí)際情況綜合確定計(jì)算范圍。

（2）操作方法

將高程點(diǎn)數(shù)據(jù)文件（XYZ.txt）及計(jì)算范圍線數(shù)據(jù)文件（BJ.txt）存放在本程序.m 文件所在文件夾內(nèi)。同時(shí),在該文件夾內(nèi)創(chuàng)建一個(gè)名稱為“Result”的excel 文件,并將該文件中的一個(gè)工作表命名為“Result”,用來存放計(jì)算結(jié)果。然后,在MATLAB 平臺(tái)中打開本程序所在.m 文件,點(diǎn)擊運(yùn)行即可開始計(jì)算。

（3）結(jié)果導(dǎo)出

本程序計(jì)算完成后將通過“xlswrite”函數(shù)把計(jì)算結(jié)果寫入Result.xls 文件內(nèi)Result 工作表中。計(jì)算結(jié)果存儲(chǔ)格式為“A B”（其中,A 列為水位數(shù)值,單位為m；B 列為庫容數(shù)值,單位為萬m3）,水位間隔即為上述設(shè)定高程計(jì)算步長b。

3 應(yīng)用實(shí)例

本文以七盞燈水庫為例,采用上述程序計(jì)算七盞燈水庫庫容曲線。七盞燈水庫地處丘陵地帶,位于廣州市番禺區(qū)沙頭街西北部大夫山森林公園內(nèi),庫長約400 m,庫水面寬130 m～155 m,匯水區(qū)內(nèi)最高點(diǎn)為大壩南西方向的大烏崗,高程226 m左右；壩址區(qū)為丘陵前的沖積洼地,兩岸低丘高程多在19.00 m～54.00 m,左岸低丘最高點(diǎn)高程110.00 m,右岸最高點(diǎn)高程88.00 m,地形坡度5°～38°；大壩下游河谷出口處地形較平緩,地面高程15.00 m～18.00 m,現(xiàn)為西大夫山森林公園。

七盞燈水庫于1958 年建成,集雨面積0.72 km2,為?。?）型水庫,工程等別為Ⅴ等,主壩級(jí)別為5 級(jí)。水庫防洪標(biāo)準(zhǔn)為20 年一遇,校核標(biāo)準(zhǔn)為50 年一遇。水庫大壩現(xiàn)狀為斜墻土壩,壩長180 m,壩高8.40 m,壩頂高程25.00 m。大壩溢洪道為單孔涵式,進(jìn)口段為5.00 m（寬）×1.95 m（高）,總長度103.20 m。大壩放水洞位于溢洪道下方,進(jìn)口暗涵內(nèi)尺寸1.50 m×1.50 m。設(shè)有一扇1.70 m×1.70 m 平板鋼閘門,配備手動(dòng)螺桿啟閉機(jī)。

七盞燈水庫曲線計(jì)算采用1∶500 實(shí)測(cè)地形。首先,根據(jù)實(shí)測(cè)地形及壩頂高程框定計(jì)算范圍。然后,利用CAD 軟件中“l(fā)ist”命令,提取計(jì)算范圍線坐標(biāo),將坐標(biāo)數(shù)據(jù)保存到BJ.txt文件中。接著,利用CAD 軟件中“數(shù)據(jù)提取”命令,提取高程點(diǎn)坐標(biāo)及高程數(shù)據(jù),將高程點(diǎn)數(shù)據(jù)保存到XYZ.txt 文件中。在地形數(shù)據(jù)預(yù)處理完成后,將BJ.txt、XYZ.txt 文件復(fù)制到本程序.m 文件所在文件夾。在MATLAB 平臺(tái)中打開本程序所在.m 文件,設(shè)置計(jì)算單元邊長d=1,高程計(jì)算步長b=0.1,并點(diǎn)擊運(yùn)行。

軟件運(yùn)行結(jié)束后,計(jì)算結(jié)果自動(dòng)存入Result.xls 文件內(nèi)Result 工作表中,計(jì)算得七盞燈水庫庫容曲線見圖7、表1。

圖7 七盞燈水庫庫容曲線

表1 七盞燈水庫庫容曲線

4 結(jié)語

本文以七盞燈水庫為例,采用基于DEM 和“填充法”理念的計(jì)算程序,方便、快速地完成了七盞燈水庫庫容曲線的計(jì)算工作。該方法不但操作簡(jiǎn)單、計(jì)算快速,而且計(jì)算精度也較高,可作為水利工程技術(shù)人員完成庫容曲線計(jì)算工作的一個(gè)手段。