基于DEM建立小型水庫流域R-V地貌瞬時單位線之初探

劉 芝

(四川省資陽市雁江區水務局，四川 資陽，641300)

1 概述

在水庫洪水預報中，傳統的預報方法是繪制水庫洪水預報綜合查算圖[1]，但該查算圖需要根據本流域大量的實測降雨流量資料進行配線，可移用性較差，必須做到一庫一圖，而對于實測資料缺乏的小型水庫，只能借用鄰近流域特性參數或地區綜合經驗參數繪制查算圖，預報結果的可信度不高。隨著社會經濟的發展和社會管理水平的提高，尤其在“智慧地球”、“智慧城市”、“智慧水利”的發展框架下，國家在逐漸筑牢大江大河防洪體系的同時，對水庫洪水預報能力的要求也在不斷提高。研究和揭示水庫洪水與流域地形地貌的關系，建立揭示流域水文響應與地形地貌關系的地貌瞬時單位線(GIUH)及具有物理基礎的分布式水文模型，是破解小型水庫洪水預報難題的必然出路。隨著計算機技術和地理信息技術的發展，尤其是數字高程模型(DEM)的廣泛應用，提取流域河網及地形地貌信息成為可便捷實現的過程。地貌瞬時單位線(GIUH)和分布式水文模型[2]在小型水庫洪水預報中的推廣和應用，成為未來的必然趨勢。

2 R-V GIUH基本理論

Rodriguez-Iturbe和Valdes于1979年首先提出了地貌瞬時單位線(GIUH)的概念，并將瞬時均勻降落在流域上的凈雨看作由無數的水質點組成，用概率論的方法來研究流域上水質點的運動，當組成凈雨的水質點間呈弱相互作用時，流域瞬時單位線等價于水質點持留時間的概率密度函數，這就是經典的R-V GIUH理論。

2.1 河流分級

R-V GIUH理論是基于Horton-Strahler河流分級法提出的。Horton-Strahler河流分級的具體方法為：將所有河網弧段中沒有支流的河網弧段定為第1級；兩個1級河網弧段匯成的河網弧段定為第2級；兩個2級河網弧段匯成的河網弧段定為第3級；如此下去，當且僅當同級別的兩條河網弧段匯成的河網弧段級別增加一級，直到河網出水口。

2.2 Horton地貌參數

R-V GIUH理論將瞬時單位線參數同流域地貌參數結合，提供了流域響應函數的地貌學解釋。自然界的水系在發育過程中遵循一定的規律。由于受地質構造和自然環境的影響，河網形態在平面上表現為有規律的排列并在不同尺度上具有一定的自相似性。Horton通過大量研究發現[3]，任何一個流域內，各級河道的數目、平均長度、平均集水面積、平均縱比降與河道的Horton-Strahler河流分級級別之間存在幾何級數關系，并據此提出了Horton地貌定律，即河數定律、河長定律、面積定律、比降定律。而Horton地貌定律中的河數率RB、河長率RL、面積率RA是推求R-V地貌瞬時單位線時重要的地貌參數。

Horton地貌參數(河數率RB、河長率RL、面積率RA)的計算公式如下：

RB=Ni-1/Nii=2，3，…，Ω

一般來說，不同級別的河流的河數率、河長率、面積率并不相同，應用時可取各級河流的平均值。且大量研究成果顯示：天然河網的Horton地貌參數存在一個較窄的變化范圍，河數率范圍為3～5，河長率范圍為1.5～3.5，面積率范圍為3～6。

2.3 R-V GIUH基本公式推導

基于Horton-Strahler河流分級方法，把水質點到達的某一級河流或坡面稱作狀態xi，把水質點從注入點到出口斷面順序經過的各狀態稱作路徑si，則水質點沿路徑si流至出口斷面的匯流時間為Tsi=Tx1+Tx2+…+TxΩ。根據流域瞬時單位線等價于水質點持留時間的概率密度函數，可導出流域地貌瞬時單位線基本表達式[4]：

μ(t)=fB(t)=∑siSfx1(t)×fx2(t)×…×fxΩ(t)·p(si)

(1)

式中：fx1(t)，fx2(t)，…，fxΩ(t)表示水質點在當前狀態下的等待時間概率密度函數；p(si)表示水質點取路徑si的概率；S表示所有路徑si的集合，有S={s1，s2，……}；Ω表示河流最高級別；符號×表示卷積相乘。

式(1)中的路徑概率p(si)有：

p(si)=θx1px1x2px2x3…pxΩ-1xΩ

(2)

式中：θxi為水質點處于初始狀態的概率，簡稱初始概率；pxixj(i

(3)

(4)

式(1)中等待時間概率密度函數fxi(t)的確定，實質上就是確定水質點流速的概率密度函數。常用下列單參數指數型函數來進行擬合：

(5)

式中：Ki為水質點在狀態xi的平均持留時間，等于狀態的空間尺度與水質點在該狀態的平均流速的比值，可用地形地貌學方法確定。

將R-V GIUH基本表達式導出為數學計算公式的過程比較復雜，且當河流最高級別Ω不同時，導出的計算公式形式亦不相同，使用時需要分河網級別予以導出。R-V GIUH理論原作者給出了三級河網的R-V GIUH計算公式，由文康[6]等于1988年給出了四級河網和五級河網的R-V GIUH計算公式。

2.4 流域平均流速的推求

(6)

式中：v為流域平均流速；L為主河道長度；τ為流域匯流時間，可用下列公式推求：

(7)

(8)

式中：ψ為洪峰徑流系數；F為流域集水面積；L為主河道長度；J為河道平均比降；S為暴雨雨力；n為暴雨公式指數；m為匯流參數。以上參數均能利用《四川省水文手冊》進行查算。

3 基于DEM提取流域河網和Horton地貌參數

基于DEM提取流域河網和地貌參數可在ArcGIS軟件中實現。ArcGIS軟件是美國環境系統研究所ESRI開發的一款無縫擴展GIS軟件，具有強大的地理信息分析處理能力。ArcGIS軟件中的Hydrology工具集能夠方便地提取河網結構等大量水文信息。

3.1 DEM的獲取

數字高程模型(DEM)是表示地面高程隨空間變化的數學模型，是對地形地貌的一種離散的數字表達。DEM的獲取方式主要有三種：一是通過互聯網上既有的DEM數據庫直接下載，常用的DEM產品有ASTERGDEM、SRTM、GLS2005、TANDEM等，分辨率多為90m和30m；二是通過現有的地形等高線圖生成DEM，即將現有地形等高線圖數字矢量化后，進行TIN建模，進而生成DEM；三是根據實測數據源生成DEM，DEM數據源的獲取方式包括地面測量、數字攝影、激光雷達、立體遙感、GPS等。普通的研究和應用采用實測數據源生成DEM的可能性不高，通常采用第一、二種DEM獲取方式，且以互聯網下載方式最為便捷。

3.2 DEM流域位置定位

互聯網下載獲得的DEM是分區塊的，在區塊DEM數據上定位到具體研究區域的位置，需要借助于行政邊界矢量數據。在ArcGIS軟件中打開DEM數據，并將下載到的行政邊界矢量數據加載到DEM上(需注意坐標系統的一致性)，然后使用ArcToolbox中的Extract by Mask工具進行剪切，即得到行政邊界區域內的DEM數據。將對應行政區域的水系圖通過空間校正疊加到DEM數據上，可定位出具體流域的位置，再沿集水區域邊界剪切，便可得到具體流域范圍的DEM數據。

3.3 DEM洼地填充處理

受DEM分辨率的限制，以及DEM生產過程中的系統誤差，原始DEM表面存在著一些凹陷或尖峰的區域，這些是DEM數據的缺陷，會導致不合理的甚至是錯誤的水流方向計算結果。但并非所有的洼地都是數據誤差造成的，也有可能是真實存在的地形。故需要進行洼地計算，計算出洼地的位置、區域、深度，以此判斷哪些洼地是數據誤差，哪些洼地是真實地形。具體計算步驟如下：①用Hydrology工具集下的FlowDirection工具計算出水流方向；②依據水流方向，用Sink工具計算出洼地位置；③用Watershed工具計算出洼地貢獻區域；④用Zonal工具集下的ZonalStatistic工具計算出洼地貢獻區域最低高程；⑤用Zonal Fill工具計算出洼地貢獻區域出口最低高程；⑥用MapAlgebra工具集下的Raster Calculator工具計算出洼地深度，洼地深度=洼地貢獻區域出口最低高程-洼地貢獻區域最低高程。依據計算出的洼地深度設置合理的洼地填充閾值，用Hydrology工具集下的Fill工具對洼地進行填平處理，得到接近于真實地形的無洼地DEM。

3.4 河網提取

河網是水在地形地貌上的產物，同時水的匯集與流動又受到河網的約束。河網特征反映一個流域的綜合水文特征。常用的河網提取方法為連續河網生成法，即采用D8法按最陡坡度確定DEM中每個網格單元的水流方向，根據水流方向計算出每個網格單元的上游集水面積(即匯流累積量)，設置一個集水面積閾值，以該閾值作為河道認定的門檻，由匯流累積量大于該閾值的網格單元組成河網。河網提取具體步驟如下：①用Hydrology工具集下的Flow Direction工具計算出無洼地DEM的水流方向；②依據無洼地DEM水流方向，用Flow Accumulation工具計算出匯流累積量；③設置適當的集水面積閾值，用Map Algebra工具集下的RasterCalculator工具生成柵格河網；④用Hydrology工具集下的Stream to Feature工具將柵格河網轉換為矢量河網。

3.5 流域劃分

流域又稱集水區域，是指由分水線所包圍的河流集水區，流經其中的水流從一個公共的出水口流出。流域的邊界可依據DEM的數字地形分析確定。流域劃分具體步驟如下：①依據提取出的流域河網，用Hydrology工具集下的Stream Link工具計算出河網結點；②依據無洼地DEM水流方向，用Basin工具劃分出流域盆地；③依據河網結點計算結果，用Snap Pour Point工具計算出匯水區出水口；④用Watershed工具生成集水流域；⑤用Conversion Tools工具集下的Raster to Polygon工具對柵格集水流域進行矢量化，得到矢量集水流域。

3.6 Horton地貌參數的提取

Horton地貌定律是以Horton-Strahler河流分級為基礎的。河網的分級包含重要的水文信息，可用以研究水流的運行和匯流模式等。河網分級具體步驟為：①根據提取的柵格河網，用Hydrology工具集下的Stream Order工具，選擇Strahler分級法，完成河網的分級；②用Stream to Feature工具對柵格分級河網進行矢量化，得到矢量分級河網。

Horton地貌參數的提取步驟具體如下：①打開矢量分級河網的屬性列表，添加列，用Calculate Geometry工具計算出河段長度，導出含河段長度的屬性列表，便可獲得每一級別河段的條數和長度等信息；②打開矢量集水流域的屬性列表，添加列，用Calculate Geometry工具計算出集水流域面積，導出含集水流域面積的屬性列表，便可獲得每一級河流的集水面積等信息；③依據Horton地貌參數計算公式，計算得到河數率RB、河長率RL、面積率RA等地貌參數值。

4 應用實例

4.1 研究區域概況

雁江區位于四川盆地中部，全區共有小型水庫78座，水庫流域面積介于0.1km2～450km2之間，普遍沒有降雨和流量觀測站，受預報方法和管理水平的限制，其水庫洪水預報工作基本上是空白的。本次研究以高板橋水庫流域作為研究對象，進行流域河網和Horton地貌參數的提取，并建立水庫流域的R-V地貌瞬時單位線。高板橋水庫位于雁江區東峰鎮境內，屬沱江水系，地理坐標為E104°53′18.3″、N30°3′50.3″，總庫容402萬m3，壩址以上集雨面積33.77km2，主河道長13.355km，河道平均坡比降1.317‰，是一座以防洪為主、兼顧農業灌溉等綜合效益的小(1)型水庫。

4.2 流域河網和地貌參數的提取

研究使用的DEM為從地理空間數據云免費下載的ASTER GDEM V2 30m分辨率DEM，雁江區的行政邊界矢量數據為從全國地理信息資源目錄服務系統免費下載獲得，流域定位所用的雁江區水系圖比例尺為1∶85000。通過圖層疊加和剪切，可分別獲得雁江區和高板橋水庫流域的DEM，見圖1。

(a)雁江區DEM (b)高板橋水庫流域DEM

按照洼地計算方法對原始DEM進行洼地計算，得到洼地的位置、區域和深度。經計算，高板橋水庫流域共有洼地21處，洼地深度最大值25m，最小值1m。由于本研究使用的原始DEM分辨率只有30m，故結合實際地形考慮，認為以上洼地均是由于數據誤差產生的偽洼地，故將所有洼地予以填平處理，得到研究區域的無洼地DEM。

設置適當的集水面積閾值，利用ArcGIS軟件中的Hydrology工具集進行流域河網提取，得到高板橋水庫流域的流域河網，見圖2。

(a)雁江區流域河網 (b)高板橋水庫流域河網

對提取出的高板橋水庫流域河網進行Strahler分級運算，獲得具有河流分級屬性的流域河網數據，并以此計算出各級河流的條數、長度、集水面積。根據Horton地貌參數計算公式，可計算得到高板橋水庫流域的河數率、河長率、面積率等地貌參數，見表1。

表1 高板橋水庫流域Horton地貌參數

可見，基于DEM提取出的高板橋水庫流域河網共分成4級，總集水面積42.0783km2，與從萬分之一地形圖上量算到的水庫壩址以上集雨面積33.77km2存在一定差異。高板橋水庫流域的河數率為4.3667、河長率為2.3821、面積率為4.9665，均在自然河流Horton地貌參數的取值范圍內，可用于進行流域R-V地貌瞬時單位線的推求。

4.3 R-V地貌瞬時單位線推求

根據《四川省水文手冊》用推理公式求得高板橋水庫流域的平均匯流時間約為8h，進而求得高板橋水庫流域平均流速約1.67m/s。

將提取的高板橋水庫流域地貌參數和計算的流域平均流速帶入四級河網流域地貌瞬時單位線公式進行計算，發現初始概率θx4<0，這是不合理的，此處采用文康[6]等推薦的處理方式，即當不滿足RB/RA<0.8這個約束條件，初始概率出現負值時，直接以不包括低一級河道集水面積在內的各級河網集水面積與流域總面積之比表示初始概率。

最后，可得出高板橋水庫流域的地貌瞬時單位線為：

μ(t)=-0.1671e-9.2492t-20.5336e-4.5271t+6.5354e-0.9442t+17.2233e-3.058t

將上述地貌瞬時單位線與適當的產流模型結合，就能實現高板橋水庫流域的水庫洪水預報。

5 結論和討論

(1)本文以雁江區高板橋水庫流域為實例，基于DEM提取出流域河網和Horton地貌參數，并成功建立了水庫流域的R-V地貌瞬時單位線，為進一步探索小型水庫流域分布式水文預報模型奠定了基礎，對推動小型水庫洪水預報方法研究具有積極意義。

(2)本研究采用中小流域暴雨洪水推理公式計算流域平均流速，計算公式中不僅包含下墊面特性參數，還充分考慮了降雨強度對流域平均流速的影響，且其參數值均可通過《四川省水文手冊》查算，具有一定的實用和推廣價值。

(3)本研究所使用的DEM為從互聯網免費下載獲得，分辨率僅為30m，且經過源數據插值處理，數據精度欠佳，以致提取出的流域參數與實際流域有一定出入。在實際的水庫洪水預報應用中，為更準確的獲取流域地形地貌信息，尚需采用精度更高的DEM產品。

(4)研究中發現，柵格河網生成時，設置不同的集水面積閾值，可得到不同密度的河網，主觀隨意性比較大。可嘗試選取不同地區不同尺度的多個流域做專題研究，探尋河網生成閾值的固有特性，尋求不同區域不同尺度下的河網生成閾值常值。

(5)R-V GIUH理論未考慮河網的水動力特性空間變異化，即認為水質點在流域中的傳播速度是固定不變的，然而實踐證明水流在河網中的傳播是一個非線性過程，且水流在從上游向下游傳播的過程中，其流速緩慢增加。同時，R-V GIUH理論幾乎完全忽略了坡面匯流過程，將坡面匯流時間忽略不計，這在坡面匯流滯時比重較大的小流域上，是無法真實模擬流域匯流過程的。故此，R-V地貌瞬時單位線用于小型水庫流域匯流計算，是存在一定理論缺陷的，尚需探索相關方面的改進。