基于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 和“填充法”理念的計算程序,方便、快速地完成了七盞燈水庫庫容曲線的計算工作。該方法不但操作簡單、計算快速,而且計算精度也較高,可作為水利工程技術人員完成庫容曲線計算工作的一個手段。