基于DEM的資料短缺地區設計洪水推求

張 婷，馮 平

(天津大學水利工程仿真與安全國家重點實驗室，天津 300072)

1 前言

資料短缺地區的水文研究從20世紀90年代開始逐步被全球的水文工作者所重視［1］。國際水文科學協會(International Association of Hydrological Sciences，IAHS)于2003年啟動了無資料(資料不足)流域(地區)的水文預報，即一個簡稱為PUB(prediction in ungauged basins)的國際水文計劃，意在未來10年全面開展資料短缺流域的水文計算研究［2］。

目前國內外學者已經對PUB展開了積極的工作。李紅霞等基于目前常用的幾種區域化方法，提出了“綜合相似法”來選擇參證流域［3］。柴曉玲等研究了IHACRES模型在資料短缺地區徑流模擬中的應用［4］。陳志明指出了地貌瞬時單位線采用地貌特征信息模擬資料短缺地區徑流時存在的一些問題［5］。在國外，典型的研究主要有用分解法進行資料短缺流域的水文預報，它實際上是我國地區綜合法的一個特例。

對于資料短缺地區的水文計算，可通過數字高程模型(digital elevation model，DEM)技術與水文模型耦合的方法，將詳細的地理信息(如高程、地理條件、土壤類型、植被、土地利用等)用于模型的參數估算［6］。這些地貌特征信息的輸入提高了資料短缺流域水文模擬模型的確定性。

2 流域概況

福州市位于福建省的東部、閩江中下游，總面積為12184 km2，海拔多為600～l 000 m，屬于典型河口盆地地貌。全區地貌類型以山地、丘陵為主，占全區總面積的72．7%。福州地區年平均氣溫為19．1 ～20．1 ℃，年平均降水量為1200 ～1700 mm［7］。

文章的研究流域為福州市北部山區，東起磨洋河，西至過溪水庫，面積約為44 km2，環繞福州江北中心城區。北部山區洪水通過過溪、新店溪、馬沙溪等11條河道排入江北中心城區水系，山區洪水與城區洪水相遇后一起排入閩江。因此對北部山區進行設計洪水的估算對江北城區的防洪排澇規劃管理至關重要。

根據水系分布特點將其劃分為11個子流域，利用DEM圖提取各子流域(見圖1)，其中過溪、新店溪、登云溪上游有過溪、八一、登云3座水庫對其水量進行調蓄。該流域有赤橋(1966—2008年)、葉洋(1970—2008年)、嶺頭(1970—2008年)等雨量站的降雨觀測資料(地理位置見圖1)和赤橋水文站的各時段最大降水量(1955—2008年)以及文山里水文站的蒸發資料(1979—2008年)可供使用，而只有赤橋水文站(1966—2008年)能提供流量資料，其位置在新店溪子流域北部，其他子流域均無流量資料?？梢?，福州市北部山區的設計洪水計算問題屬于資料短缺地區的設計洪水計算。

圖1 研究流域DEMFig．1 DEM of the research basin

3 方案設計

由于研究流域內徑流資料短缺，雨量資料相對充足，故采用通過暴雨資料推求設計洪水的方法。本流域地處南方濕潤地區，植被良好，雨量充沛，地下水豐富，符合蓄滿產流機制條件，故可采用新安江三水源模型計算產流。具體研究方案如下:

3．1 設計暴雨計算

根據短歷時最大降水量統計資料，采用Pearson-Ⅲ型曲線進行暴雨頻率分析［8］。給出不同標準下(頻率為2%、5%、10%和20%)的各時段暴雨設計值。然后選取典型暴雨，根據設計暴雨值和典型暴雨的時程分配推求不同頻率下的設計暴雨過程。

3．2 產流計算方案

采用新安江三水源蓄滿產流模型。由于赤橋流域與研究的各子流域地理位置相近，而且地形、地貌、植被、土壤、水文、地質條件等均非常相似，因此，模型參數可通過赤橋水文站的流量資料進行率定，再移植到其他子流域。

3．3 匯流計算方案

根據凈雨在流域上匯流的途徑不同，將匯流過程分成地面匯流、壤中流匯流和地下匯流3個部分。壤中流和地下徑流分別采用不同的線性調蓄水庫模擬其匯流過程。由于缺少徑流資料，文章通過DEM推求等流時線，然后采用等流時線法模擬流域地面匯流過程。

4 設計暴雨推求

采用Pearson-Ⅲ型曲線，由北部山區短歷時最大降水量統計資料，分別繪制其最大1 h、3 h、6 h、12 h和24 h降雨量頻率分布曲線，如圖2是最大24 h降雨量的頻率分布曲線。然后，給出各時段相應頻率為2%、5%、10%和20%的設計暴雨值(見表1)。

圖2 北部山區最大24 h降雨量頻率分布曲線Fig．2 The frequency distribution curve of the biggest 24 h rainfall in northern mountain area

選取1966年9月3日10:00至9月4日10:00的24 h暴雨為典型暴雨。選擇原因為:該場降雨雨量大，強度也大，暴雨核心部分出現在后期，形成的洪水主峰出現較遲。根據設計暴雨值和典型暴雨的時程分配推求出不同頻率下的設計暴雨過程。

5 模型計算

5．1 產流計算

新安江三水源產流模型主要由三部分組成:蓄滿產流計算、流域水源劃分和流域蒸散發計算［9］。

5．1．1 蓄滿產流

設R是凈雨量，降雨量為P，蒸發量為E，蓄水容量為Wm，降雨開始時的實際蓄水量為W，對于流域中某點而言，蓄滿前，R=0;蓄滿后:

表1 不同頻率下的設計暴雨值Table 1 Design rainstorm in different frequencies

5．1．2 劃分三水源

設RS、RSS、RG分別為地面徑流、壤中流、地下徑流，則有:

5．1．3 蒸散發計算

設E為流域蒸發量，EU、EL、ED分別為上層、下層、深層蒸發量，則有:

5．2 匯流計算

5．2．1 地面徑流匯流

采用等流時線法。隨著流域的地形地貌逐漸可用DEM來表示［10］，于是便可以通過流域的DEM而不是地圖來自動生成等流時線分布圖［11］。具體提取的方法如下［12］。

1)首先對原始DEM數據進行洼地填充，得到無洼地的 DEM［13］。

2)通過Arcgis軟件或者Swat軟件提取各個子流域。

3)以柵格為單位，通過Arcgis軟件的Flow Direction工具提取水流方向，用Slope工具提取坡度。

4)通過經驗公式來計算流經每個柵格的速度vi，如:

式(4)中，i為柵格序號;Si為某一水流出流方向上的坡度;a為參數，具有速度的量綱;b為冪指數，反映坡度大小對流速的影響。在利用式(4)來計算流經柵格的速度時，必須先合理地確定公式中的參數a和b。方法有兩種，一種是通過水力學方法來直接測定，另外一種是通過水文學中的單位線理論來間接率定［14］。文章采用第二種方法，利用赤橋水文站資料率定得到參數 a=1．55、b=0．3，然后移用到其他子流域。

5)根據式(4)計算出速度圖層后，利用柵格計算器計算該圖層的倒數，即1/vi。各個柵格的匯流時間可以通過提取水流長度(Flow Length工具)，在“輸入權重光柵”(Input weight raster)一欄中選擇1/vi圖層作為權重得到。各個子流域的匯流時間為15～40 min，按照5 min的間距reclass進行分類，通過查詢屬性表中每個時段的柵格個數即可得到各個子流域的等流時面積。圖3為第5子流域的等流時面積分布。

圖3 第5子流域等流時面積分布Fig．3 The isochrones-divided area of the fifth sub basin

等流時線匯流具體計算方法:假設流域被劃分為N塊等流時面積，記為fi(i=1，2，…，N)，一場空間分布均勻的凈雨有M個時段，記為hj(j=1，2，…，M)，計算時段及等流時線時均距為Δt，此處Δt=5 min，出口斷面有流量的時段數 T=M+N－1，每個時段的流域出口斷面流量(地面徑流)記為QS(k)(k=1，2，…，T)，計算公式為:

5．2．2 壤中流和地下徑流匯流

壤中流和地下徑流分別采用線性水庫調蓄模型和線性水庫蓄泄模型模擬其匯流過程，設QSS(k)、RSS(k)分別為第k個時段的壤中流流量和平均壤中流凈雨深，設Qg(k)、Rg(k)分別為第k個時段的地下徑流流量和平均地下凈雨深，KKSS、KKG分別為壤中流的消退系數和地下徑流的消退系數，f是子流域面積，則有:

5．2．3 總徑流量

第k個時段總徑流量Q(k)的計算式為:

5．3 模型參數率定

選取赤橋水文站1979—2008年中較典型的15場洪水進行參數的率定和驗證(前10場用于率定，后5場用于驗證)，得到的模型參數見表2。圖4為其中900731次洪水過程的模擬結果。圖5為第5子流域不同標準下的設計洪水。表3是赤橋站模擬與實測峰量對比情況。

表2 模型參數率定值Table 2 Adjusted parameters of the model

續表

圖4 900731次洪水過程模擬Fig．4 Simulation of the 900731 flood process

圖5 第5子流域不同標準下的設計洪水Fig．5 The design flood in different frequencies of the fifth sub basin

表3 實測與模擬峰量對比表Table 3 The peak flow and quantity of floods obtained by measurement and simulation

表3結果表明徑流的合格率較為理想，其中10場洪水峰值相對誤差在±10%以內，12場洪水在±20%以內，平均相對誤差為10%。另外，其中有7場洪水洪量相對誤差在±10%以內。多數場次洪水的確定性系數在0．65以上?？傮w來說，模擬效果比較好，其中洪峰的模擬比洪量要好，這更有利于小流域設計洪水的計算。

6 設計洪水推求

根據計算得到的不同頻率下的設計暴雨過程，利用產流參數在各子流域上進行產流計算，然后利用筆者建議的基于DEM方法進行匯流計算，就可以給出各子流域的設計洪水過程，其洪峰和洪量結果如表4所示。

表4 各子流域設計洪峰洪量統計表Table 4 The design peak flow and quantity of floods in the sub basins

7 對比驗證

為了驗證表4給出的設計洪水的合理性，也采用水文比擬法計算了各子流域的設計洪峰值。根據赤橋水文站資料分析，流域洪峰流量隨面積的遞減指數n=0．835。這樣，選擇赤橋站為參證站，利用式(9)就可以推求不同頻率的設計洪峰流量(見表5)。

表5 水文比擬法與DEM方法計算洪峰流量結果對比(單位:m3/s)Table 5 The result of flood peak flow calculated by hydrological analogy approach and DEM approach(unit:m3/s)

對比二者計算結果:當設計標準較低時，基于DEM方法的結果偏大;當設計標準較高時，水文比擬法結果偏大，但是二者基本相近。水文比擬法因未考慮地形、地質、人類活動等的影響，存在一定缺陷。因此，對于資料短缺地區的設計洪水計算，基于DEM的方法可能更合理一些。

8 結語

以缺少徑流資料的福州市北部山區為研究流域，根據水系特點將其劃分為11個子流域，研究了資料短缺地區的設計洪水推求方法。根據流域資料情況，建議通過暴雨資料推求設計洪水，主要做了以下工作:

1)結合研究流域的實際情況，通過DEM技術與水文模型相結合的方法，探討資料短缺地區設計洪水的計算方法。

2)通過參證站水文模型參數的移用，解決流域內其他子流域的產匯流參數的估算問題。特別是在匯流計算中，將近年來發展非常迅速的地理信息系統工具運用到研究當中，通過DEM推求等流時線，避免了在地圖上人工勾繪等流時線的繁瑣工作。隨著全國各地數字高程模型的建立，將其運用到水文計算當中有重要的現實意義。

文章采用的這種推求設計洪水的方法為其他資料短缺地區的設計洪水推求提供了一個范例。另外，通過DEM推求等流時線的方法也可與其他水文模型相結合，用于其他流域的徑流模擬和設計洪水的推求。

［1］談 戈，夏 軍，李 新．資料短缺地區水文預報研究的方法與出路［J］．冰川凍土，2004，26(2):192－196．

［2］劉蘇峽，夏 軍，莫興國．資料短缺流域水文預報(PUB計劃)研究進展［J］．水利水電技術，2005，36(2):9 －12．

［3］李紅霞，張永強，敖天其，等．無資料地區徑流預報方法比較與改進［J］．長江科學院院報，2010，27(2):11－15．

［4］柴曉玲，郭生練，彭定志，等．IHACRES模型在無資料地區徑流模擬中的應用研究［J］．水文，2006，26(2):30－33．

［5］陳志明．流域地貌瞬時單位線法剖析［J］．水電能源科學，1993，11(6):105 －111．

［6］李紅霞．無徑流資料流域的水文預報研究［D］．大連:大連理工大學，2009．

［7］戴志忠．福州地區水文特性分析［J］．水利科技，2004(3):9－11．

［8］詹道江，葉守澤．工程水文學［M］．北京:中國水利水電出版社，2007．

［9］趙人?。饔蛩哪M［M］．北京:水利電力出版社，1984．

［10］任立良，劉新仁．數字高程模型在流域水系拓撲結構計算中的應用［J］．水科學進展，1999，10(2):129 －134．

［11］Maidment D R．Developing a spatially distributed unit hydrograph by using GIS［C］//Kovar K，Nachtnebel H P．Hydro GIS 93，Application of GIS in Hydrology and Water Resources:Proceedings of the Vienna Conference．IAHS，1993(211):181 －192．

［12］熊立華，彭定志．基于數字高程模型的等流時線推求與應用［J］．武漢大學學報(工學版)，2003，36(3):1 －3．

［13］湯國安，楊 昕．Arcgis地理信息系統空間分析實驗教程［M］．北京:科學出版社，2006．

［14］熊立華，郭生練．分布式流域水文模型［M］．北京:中國水利水電出版社，2004．