基于DEM的空間流速場構建方法及應用

謝志剛，石 朋,2，紀小敏，趙蘭蘭，瞿思敏，崔彥萍

(1.河海大學水文水資源學院，南京 210098；2.河海大學 水文水資源與水利工程科學國家重點實驗室，南京 210098；3.江蘇省水文水資源勘測局，南京 210029；4.水利部信息中心，北京 100053)

已有研究表明[1,2]，流域匯流可被視為凈雨在地貌擴散和水動力擴散綜合作用下的匯集過程。Redriguez-Iturbe和Gupta等[3,4]提出假定速度項表現為指數分布，則其線性反應函數為時間指數分布，由此流域地貌瞬時單位線可以解釋為一滴雨抵達流域上任何處，其匯流時間的概率密度函數。芮孝芳等[5]研究發現，確定等待時間概率密度函數的實質就是確定水質點流速的概率密度函數，故把推求地貌瞬時單位線問題轉化為確定流速的概率密度函數問題。

計算流速的方法有很多，主要的流速計算公式包括曼寧公式、梅德門特公式、謝才公式等[6]。1993年，梅德門特教授[7]基于柵格假定了一個時不變的空間流速場；1995年，Muzik[8]在梅德門特教授的基礎上，結合曼寧公式與連續方程進行流速計算，并利用GIS確定流速場；2003年，Assefa M.Melesse等[9]人總結出結合曼寧公式與運動波近似方程來估算坡面流速，聯合連續方程與曼寧公式求解河道流速，以此構建流速場；2006年，石朋等[10,11]采用曼寧公式計算柵格單元流速并以此建立了一個在空間上變化而時間上不變的空間流速場，以沿渡河流域為研究對象，取得了較好的計算效果?？追舱艿萚12]區分坡地單元與河道網格單元，以SCS公式計算坡地流速，以梅德門特公式計算河道流速，由此計算出整個流域的流速分布?？紤]到在一次降雨徑流過程中，雨強對流速的影響是不容忽視的[13]，本文擬通過將雨強納入考慮范圍并改進原有的曼寧公式，計算流域各網格內流速，并應用于定安河流域，以提高洪水模擬過程的精度。

1 研究思路

1.1 流速計算方法

流域下墊面條件存在空間差異性，不同地形地貌條件下流域各處的流速大小與方向各不相同。流速大小的計算是本文的核心所在。本文主要選用了兩種流速計算方法，分別是：傳統曼寧公式及考慮雨強的曼寧公式。

1.1.1 曼寧公式

傳統的曼寧公式如下：

(1)

式中：n為地表粗糙系數，可以通過土地利用情況獲得；S為地表坡度，可以通過DEM獲得；R為水力半徑，可以近似以水深代替。

由已有研究可知[14]，水深的計算可以建立水深與匯水面積的函數關系：

H=φpAψp

(2)

式中：H為具有超越概率p的柵格平均水深；A為柵格上游流域匯水面積；φp為網格常數，與洪水頻率有關；ψp為幾何常數，與洪水頻率有關。

1.1.2 考慮雨強的曼寧公式

實際洪水過程中存在降雨過程，不同的降雨強度對匯流速度也存在影響，因此采用考慮雨強的曼寧公式進行流速計算，可提高計算精度。計算公式如下：

(3)

式中：i為雨強大小的無量綱因子；b為經驗系數。

雨強作為影響流速大小的影響因子納入曼寧公式中，如公式(3)所示。在GIS中利用克里金插值法得出流域雨強柵格圖，即可知每個柵格內的雨強大小，再由已知的糙率系數、坡度和水力半徑通過公式(3)可求出受雨強影響的柵格流速。由于整個洪水過程中存在降雨空間分布不均勻的情況[15]，不同地區雨強大小不同，并且在同一時段可能只有部分地區降雨，因此要對雨強柵格圖進行0處理，即利用柵格計算器把所有0柵格替換成1，在無降雨地區i取值1，等價于直接用曼寧公式計算。

1.2 匯流時間計算

通過上述計算所得柵格流速，依據柵格大小可以計算出水流在柵格內的滯留時間。計算公式如下：

Δτ=L/v

(4)

(5)

式中：Δτ為柵格內滯留時間；L為柵格的邊長。

公式(4)表示的是水流方向同柵格邊線平行的情況，公式(5)表示的是水流方向同柵格對角線平行的情況。在GIS中通過D8算法可以確定水流方向，即確定每個柵格內水流方向，通過公式(4)和公式(5)計算各柵格內滯留時間，在GIS中沿水流方向累積各柵格到達流域出口的時間即可確定各柵格匯流累積時間，計算公式如下：

(6)

式中：τ為某柵格到達流域出口的累積匯流時間；m為某柵格徑流路徑上的柵格數。

1.3 單位線的確定

通過以上兩種流速計算方法在GIS中可利用柵格計算器計算出每個柵格的流速大小，以此可得出兩種不同的空間流速場，利用空間流速場以公式(4)、(5)可以計算得到每個柵格的滯留時間，通過公式(6)沿水流方向得出各柵格匯流累積時間，在GIS中利用重分類模塊統計時段內通過流域出口的柵格數，即得到時段內出流面積。以匯流時間為橫坐標，以時段內出流面積與總面積比值為縱坐標，可得到流域匯流時間的概率密度分布。根據徑流過程形成的“粒子學說”[2-4]，即水質點在弱相互作用情況下的匯流時間概率密度分布函數等價于地貌瞬時單位線，可知以上所得匯流時間概率密度分布即為地貌瞬時單位線，將此應用于實例當中，對比徑流模擬計算結果，分析兩種流速計算方法存在的問題與優點。

2 實例應用

2.1 流域資料處理

本文所選流域是萬泉河水系內的定安河流域，流域面積1 333 km2，內設長征、大墩、合羅、加報、羅擔、木色、瓊中、石古、思河、烏坡和烏石11個雨量站，具有完善的降雨資料。流域地處熱帶季風氣候區，常年雨量充沛，多年平均降雨量約1 639 mm。流域是萬泉河水系的一級支流，受水利工程影響較小，具有多年連續的降雨徑流資料，適合本文模擬研究計算。

本文基于流域DEM(90×90)進行空間流速場的構建，推導出相應的地貌瞬時單位線，并應用于降雨徑流模擬分析。流域坡度分布如圖1所示。

圖1 流域坡度分布Fig.1 Slope distribution of River Basin

根據流域土地利用情況確定了糙率系數n的柵格數據圖；利用GIS中水文分析模塊可計算各柵格上游匯水面積，以此根據公式(2)可求出流域水深柵格數據圖；通過克里金插值法得出流域時段雨強柵格數據圖。利用GIS中柵格計算器，以上述得出的柵格數據通過公式(1)～(3)計算每個柵格的流速，分別構建出以曼寧公式計算出的空間流速場和以考慮雨強的曼寧公式計算出的空間流速場，雨強曼寧公式以20010825次洪為例。圖2(a)為曼寧公式構建出空間上變化而時不變流速場，圖2(b)為2001年8月25日次洪降雨第2個時段雨強的空間流速場柵格圖。

圖2 計算所得空間流速場Fig.2 Calculated spatial velocity field

根據匯流時間計算公式(4)～(6)可以計算出相應的匯流累積時間柵格分布圖，圖3(a)為曼寧公式計算出的累積時間柵格圖，圖3(b)為2001年8月25日次洪降雨第2個時段雨強的計算出的累積時間柵格圖。

圖3為沿水流方向到達流域出口的匯流累積時間分布圖，距離匯流出口越遠所需匯流時間越長，與實際情況相符合，據此統計匯流時間的概率密度分布函數[16]，得出地貌瞬時單位線。雨強曼寧公式由于各時段雨強不同所得出的單位線也不同，故采用時段出流方式進行匯流計算。

圖3 匯流累積時間分布Fig.3 Accumulation time distribution of confluence

2.2 降雨徑流模擬計算

利用上述所得地貌瞬時單位線轉換成1 h時段單位線，應用于定安河流域進行降雨徑流模擬分析，選取該流域2000-2013年中的7場洪水進行模擬計算，流域每年洪澇災害多發生于5-11月之間，強降雨主要集中在7-10月，考慮雨強對流速的影響，故選取此期間7場洪水進行模擬計算，所得計算結果如表1所示。經過對比分析，由曼寧公式計算值相比較于考慮雨強的曼寧公式計算值，7場洪水平均徑流相對誤差由11.52%下降到5.93%，平均洪峰相對誤差由6.32%下降到3.85%，平均確定性系數由0.761上升到0.840。由此可知考慮雨強大小影響的曼寧公式計算所得結果更符合實際情況，其徑流深、洪峰流量、確定性系數模擬效果更優[17]。圖4為20010825和20071011兩場洪水由曼寧公式計算的洪水過程線和考慮雨強的曼寧公式計算的洪水過程性與實測洪水過程線的對比圖。

表1 定安河徑流計算結果Tab.1 Runoff calculation results of the Ding An River

注：此平均值為絕對值平均值。

圖4 洪水模擬計算對比Fig.4 Comparison of flood simulation

3 結 語

本文核心是在流域流速計算方法上，而傳統地貌瞬時單位線的確定是建立在流速空間分布均勻的情況下，實際流域中流速大小受地形地貌影響以及水力條件影響較大，故本文提出了考慮雨強影響的曼寧公式計算流域流速，建立了整個流域的空間流速場。曼寧公式充分考慮了土地利用情況、坡度以及水深對流速的影響，由此可以確定一個空間上變化而時不變的空間流速場，但是實際洪水過程中存在降雨空間分布不均勻性，不同區域雨強不同亦對流速產生重大影響，故此加入雨強這一變量因素，以此改進曼寧公式求得流域空間流速場，應用于定安河流域，取得了比較好的模擬效果，故此可推廣應用于缺資料地區的匯流計算過程中。

本文所取雨強為1 h時段雨強，未考慮不同時段大小雨強對流速的影響，后期重點研究不同時段大小雨強對流速計算精度的影響，提高徑流模擬效果。本文所取柵格大小是90 m×90 m，未考慮不同分辨率大小情況下坡度值精度不同，由此產生流速計算誤差，故在后續研究過程中會進一步考慮分辨率對流速計算產生的影響。