黃國(guó)新 謝小華 鄧 蕊

(1.江西省水文局，江西 南昌 330000；2.江西省吉安市水文局，江西 吉安 343000)

江西省位于長(zhǎng)江中游南岸，河網(wǎng)密布水系發(fā)達(dá)，獨(dú)特的地形氣候條件造成該地區(qū)洪澇災(zāi)害頻繁，洪水預(yù)報(bào)對(duì)水安全至關(guān)重要。但由于部分地區(qū)水文站網(wǎng)缺失，水文監(jiān)測(cè)資料缺乏，給洪水預(yù)報(bào)、水庫(kù)調(diào)度及安全度汛工作帶來(lái)了很大的挑戰(zhàn)。吉安市行政區(qū)劃范圍內(nèi)有大中型水庫(kù)47座，其中，無(wú)長(zhǎng)期水文觀測(cè)站點(diǎn)的水庫(kù)有36座，占全部大中型水庫(kù)數(shù)量的76.6%。為便于上述無(wú)資料地區(qū)進(jìn)行洪水預(yù)報(bào)，開(kāi)展地貌單位線適用性研究很有必要。

傳統(tǒng)的流域單位線從系統(tǒng)水文的角度進(jìn)行研究，把流域看作一個(gè)系統(tǒng)，降雨和出流分別視為輸入和輸出。因此，流域單位線成了黑箱模型的產(chǎn)物。流域地貌瞬時(shí)單位線是以流域地形地貌情況以及概率學(xué)隨機(jī)理論為基礎(chǔ)，具有物理意義的流域匯流隨機(jī)模型。流域地貌瞬時(shí)單位線以研究區(qū)獨(dú)特的地形地貌特點(diǎn)，揭示流域基本特性與單位線之間的關(guān)系；以研究區(qū)水文氣象特性來(lái)體現(xiàn)與單位線的非線性時(shí)變之間的關(guān)系，克服了神經(jīng)網(wǎng)絡(luò)等學(xué)習(xí)方法推求單位線的一些缺陷，與研究區(qū)現(xiàn)實(shí)情況更為貼切，更能反映區(qū)域特點(diǎn)。自1979年委內(nèi)瑞拉的水文學(xué)者提出該模型以來(lái)，近幾十年來(lái)得到了迅速發(fā)展[1]。尤其是隨著計(jì)算機(jī)技術(shù)和GIS技術(shù)的高速發(fā)展，同時(shí)伴隨數(shù)字高程信息化的誕生，為研究地貌瞬時(shí)單位線與研究區(qū)地形地貌特征之間的相互關(guān)系，提供了非常有效的技術(shù)支持和保障，更加快速地促進(jìn)了地貌瞬時(shí)單位線在區(qū)域洪水預(yù)報(bào)工作中的實(shí)踐與研究。地貌單位線采取數(shù)字高程模型(Digital Elevation Model)、地理信息系統(tǒng)(Geographic Information System)、遙感(Remote Sensing)等相關(guān)先進(jìn)科學(xué)手段，使我們可以實(shí)時(shí)獲得研究區(qū)的地形地貌，同時(shí)充分考慮了人類因素的影響導(dǎo)致的流域時(shí)空變化，讓研究更接近流域的實(shí)際情況。受流域面積及形狀的限制，地貌單位線一般適用于50～1000km2的流域范圍。本文選取吉安市南車水庫(kù)流域?yàn)檠芯繉?duì)象，基于DEM 數(shù)字高程信息，提取流域內(nèi)河網(wǎng)基本情況，獲取地貌參數(shù)，然后建立流域匯流模型并驗(yàn)證分析[2]。

1 研究區(qū)概況及數(shù)據(jù)來(lái)源

1.1 研究區(qū)及水庫(kù)概況

南車水庫(kù)位于贛江水系禾水一級(jí)支流牛吼江上，地處吉安市東南部，坐落于泰和縣橋頭鎮(zhèn)大壟坑村，興建于1992年9月，2003年6月建成蓄水，水庫(kù)壩上斷面以上集水面積459km2，流域范圍內(nèi)匯集有中朝水、少江水、六八河、高幣水等河流，覆蓋新江鄉(xiāng)、雙橋鄉(xiāng)和五斗江鄉(xiāng)等鄉(xiāng)鎮(zhèn)，于大壟坑村匯入南車水庫(kù)。研究區(qū)控制河長(zhǎng)57.9km，流域范圍見(jiàn)圖1。

圖1 研究區(qū)范圍

南車水庫(kù)設(shè)計(jì)總庫(kù)容15380萬(wàn)m3，興利庫(kù)容9520萬(wàn)m3，防洪庫(kù)容415萬(wàn)m3，死庫(kù)容2810萬(wàn)m3，是一座以防洪、灌溉、發(fā)電等功能為主的混凝土面板碾壓堆石壩大(2)型水庫(kù)。

1.2 數(shù)據(jù)與研究方法

1.2.1 DEM數(shù)據(jù)

數(shù)字高程圖是數(shù)據(jù)分析的基礎(chǔ)，本文數(shù)字高程數(shù)據(jù)來(lái)源為地理空間數(shù)據(jù)云GDEM 30m分辨率數(shù)字高程數(shù)據(jù)。數(shù)據(jù)獲取后由ArcGIS提取處理[3]。

1.2.2 降雨徑流資料

本研究選取南車水庫(kù)流域2017—2018年6場(chǎng)降雨及徑流資料進(jìn)行模型分析及驗(yàn)證，經(jīng)過(guò)分析和計(jì)算提取相關(guān)指標(biāo)，見(jiàn)表1。

表1 降雨及徑流資料

降雨資料來(lái)源于流域范圍內(nèi)測(cè)站上報(bào)的實(shí)時(shí)水雨情數(shù)據(jù)庫(kù)，流量、峰現(xiàn)時(shí)間值根據(jù)庫(kù)容變化量按水量平衡原理倒推入庫(kù)流量得到。

徑流系數(shù)計(jì)算公式如下：

(1)

式中：W為場(chǎng)次洪水實(shí)測(cè)總量，m3；F為流域面積，km2；P為流域面降雨量，mm。

通過(guò)查閱《吉安市水文手冊(cè)》，南車水庫(kù)流域范圍內(nèi)降雨歷時(shí)一般在35h以內(nèi)，徑流系數(shù)在0.5～0.8。由表1可知，6場(chǎng)次流域降雨歷時(shí)均值為22h，徑流系數(shù)均值為0.6，所選場(chǎng)次洪水能夠反映流域的實(shí)際情況。

1.2.3 研究方法

地貌瞬時(shí)單位線是既有物理意義又簡(jiǎn)便易行的一種匯流模型，將其引進(jìn)水文模型中，直接參與水文計(jì)算。其基本思路是：設(shè)有多個(gè)(n個(gè))體積、質(zhì)量相等且互相沒(méi)有聯(lián)系、作用的水滴，在很短的時(shí)間匯到研究區(qū)，水滴在運(yùn)動(dòng)軌跡內(nèi)要經(jīng)過(guò)相應(yīng)的運(yùn)移時(shí)間才移動(dòng)到流域出口斷面位置，每個(gè)水滴都有自己相應(yīng)的運(yùn)移時(shí)間TB，但它們卻擁有一致的分布函數(shù)。基于區(qū)域內(nèi)的特定時(shí)間注入的水量與流出的水量之差等于水量變化量(水量動(dòng)態(tài)平衡原理)及大數(shù)法則，初步判斷流域瞬時(shí)單位線(IUH)相當(dāng)于水滴在原來(lái)位置到流域出口斷面運(yùn)移時(shí)間的概率密度函數(shù)。在天然下墊面條件下，無(wú)數(shù)下落到流域內(nèi)的水滴順著自己的軌跡，通過(guò)填坑、補(bǔ)洼后匯集到達(dá)流域出口斷面。每個(gè)水滴途經(jīng)的線路都有無(wú)數(shù)的方案，這些結(jié)果的概率函數(shù)可以采用Strahler以及Horton方法求得。計(jì)算全部結(jié)果的運(yùn)行軌跡的隨機(jī)運(yùn)行時(shí)間的概率密度函數(shù)同此運(yùn)動(dòng)軌跡的函數(shù)概率相乘，然后求和，即為區(qū)域內(nèi)水滴滯留時(shí)間的概率函數(shù)，就是本次應(yīng)用的洪水預(yù)報(bào)方法——地貌瞬時(shí)單位線(GIUH)[4]。

假設(shè)某一時(shí)間降落到研究區(qū)內(nèi)分布均勻的降水量是由無(wú)數(shù)微小的水滴構(gòu)成，且水滴之間僅為很弱的相互影響和作用，則每個(gè)水質(zhì)點(diǎn)的運(yùn)動(dòng)都可被看作是一個(gè)由轉(zhuǎn)移概率控制的從某級(jí)別河道向同級(jí)或高級(jí)河道運(yùn)移的Markov過(guò)程，那么該流域的瞬時(shí)單位線——u(0,t) 同水滴匯集到流域出口時(shí)間的分布密度函數(shù)fB(t)相等，得到下式：

u(0,t)=fB(t)

(2)

據(jù)此，可以判定該流域匯流的S曲線S(t)與水滴匯集到流域出口時(shí)間的概率分布函數(shù)FB(t)相等，得到下式：

S(t)=FB(t)

(3)

式(2)、式(3)表明，只要求得水滴匯集到流域出口時(shí)間的概率分布密度或概率分布函數(shù)二者中的其一，就能夠確定流域的瞬時(shí)單位線或S曲線，與傳統(tǒng)流域單位線的分析方式不同的是，地貌單位線主要基于區(qū)域地形地貌要素等角度與流域水文之間的相互影響，通過(guò)以上方式求得[1]。

基于水滴匯流時(shí)間T等于其流路長(zhǎng)度L與其速度V之商的基本關(guān)系式，即T=L/V，根據(jù)概率論的理論推導(dǎo)得出下式：

(4)

(5)

式中：g(l)與φ(v)分別為L(zhǎng)與V的分布密度；G(l)與ψ(v)分別為L(zhǎng)與V的分布函數(shù)；vmax為流域中水滴的最大匯集速度。

式(4)、式(5)將確定流域地貌瞬時(shí)單位線或S曲線的問(wèn)題，轉(zhuǎn)化為推求水滴向流域出口斷面匯集的流路長(zhǎng)度和速度的分布函數(shù)或分布密度問(wèn)題。水滴匯流的流路長(zhǎng)度可通過(guò)DEM直接提取得到，而由于受地形和地貌的影響，匯流速度在空間分布上具有不均勻性，很難直接確定。根據(jù)《美國(guó)水土保持局工程手冊(cè)》中提出的方法，參照Manning公式的形式，根據(jù)DEM單元的坡度計(jì)算流速：

V=asb

(6)

式中：s為單元的平均坡度；a和b為經(jīng)驗(yàn)參數(shù)。

因此通過(guò)提取坡度分布律的方式來(lái)確定水滴匯流速度的分布律。以此為基礎(chǔ)，將式(5)改為離散形式，可得

(7)

式中：m為離散的坡度值的總數(shù)；Θ(s)為坡度的分布函數(shù)。

求得FB(t)后，就能夠得到某一Δt時(shí)段的地貌單位線：

U(Δt，t)=S(t)-S(t-Δt)=FB(t)-FB(t-Δt)

(8)

式中：U(Δt，t)為計(jì)算時(shí)段為Δt的地貌單位線。

2 基于DEM的流域地貌單位線獲取

2.1 南車水庫(kù)流域DEM的生成

提取得到的流域分水線(山脊線)是水源的起點(diǎn)，所以其匯流累積量為0。因此，在Arcmap中使用水文模塊對(duì)區(qū)域的匯流累積量的0值進(jìn)行提取，并將提取的數(shù)據(jù)與生成的正地形求交，得到消除了存在于負(fù)地形區(qū)域中的錯(cuò)誤的山脊線，獲得流域的分水線。

在Arcgis中創(chuàng)建一個(gè)線圖層，以生成的分水線作為參考圖層，結(jié)合南車水庫(kù)流域及其嵌套流域地貌圖，勾勒出流域邊界，并將邊界圖層轉(zhuǎn)化為面域圖對(duì)DEM圖進(jìn)行裁剪，得到南車水庫(kù)流域的DEM圖，見(jiàn)圖2。

圖2 南車水庫(kù)流域邊界及流域DEM圖

2.2 無(wú)洼地DEM生成

考慮到內(nèi)插函數(shù)存在一定的局限性和研究區(qū)內(nèi)某些非常規(guī)地形的約束，通常會(huì)造成人為內(nèi)插出的DEM高程面有一些凹陷地區(qū)，在對(duì)這些區(qū)域進(jìn)行水文分析計(jì)算和水滴運(yùn)移軌跡推算時(shí)，容易生成不合理或者錯(cuò)誤的水力條件，進(jìn)而獲取錯(cuò)誤的水流方向信息和匯流時(shí)間。所以，在對(duì)流域內(nèi)各水滴進(jìn)行水力分析和水文計(jì)算之前，要優(yōu)先對(duì)提取的原始DEM高程數(shù)據(jù)進(jìn)行處理，使得各洼地得以填充，得到處理后的研究區(qū)的DEM高程圖[5]。

2.3 河網(wǎng)提取

利用地理信息系統(tǒng)Arcgis10.2，利用基于流向的河網(wǎng)提取方法，得到柵格河網(wǎng)數(shù)據(jù)，見(jiàn)圖3。

圖3 南車水庫(kù)流域河網(wǎng)

2.4 地貌瞬時(shí)單位線計(jì)算

地貌單位線預(yù)報(bào)方案中河流級(jí)別的劃分采用斯特拉勒法，它是基于修正的霍頓河流級(jí)別劃分方式而得出的。斯特拉勒分級(jí)法定義從河源出發(fā)的河流為1級(jí)河流，同級(jí)的兩條河流交匯所形成的河流級(jí)數(shù)增加1級(jí)，不同級(jí)別的兩條河流交匯所形成的河流級(jí)別為二者中較高者。

選取地理信息系統(tǒng)操作軟件Arcgis10.2，利用已經(jīng)提取的河網(wǎng)信息和參數(shù)，運(yùn)用斯特拉勒分級(jí)法對(duì)河網(wǎng)進(jìn)行分級(jí)運(yùn)算，獲取具備所研究流域相應(yīng)級(jí)別的河網(wǎng)屬性數(shù)據(jù)，以此為基礎(chǔ)可獲取流域各級(jí)別河段的數(shù)量、流程、集水面積和其他與地貌單位線推求相關(guān)的各項(xiàng)地貌特征值，同時(shí)運(yùn)用霍頓河系定律進(jìn)行計(jì)算就能夠獲取流域內(nèi)的河數(shù)率、河長(zhǎng)率和面積率，其表達(dá)式分別為

河數(shù)率：

(9)

河長(zhǎng)率：

(10)

面積率：

(11)

各河道具體數(shù)據(jù)見(jiàn)表2，劃分的集水區(qū)域矢量圖見(jiàn)圖4。

表2 各河道屬性數(shù)據(jù)

圖4 集水區(qū)域矢量圖

根據(jù)地貌單位線通用公式，進(jìn)行瞬時(shí)單位線計(jì)算。為簡(jiǎn)化工作量，本文采用編程方式，將上述計(jì)算過(guò)程封裝為計(jì)算工具。把河道三參數(shù)輸入工具，計(jì)算得到本流域的地貌瞬時(shí)單位線，見(jiàn)圖5，進(jìn)而得到流域1h地貌瞬時(shí)單位線[6]。

圖5 地貌瞬時(shí)單位線計(jì)算

3 基于地貌單位線的洪水預(yù)報(bào)分析

為便于洪水預(yù)報(bào)分析，將地貌單位線導(dǎo)入到江西省吉安市水庫(kù)防汛調(diào)度預(yù)報(bào)系統(tǒng)中，該系統(tǒng)能夠自動(dòng)讀取歷史、實(shí)時(shí)數(shù)據(jù)，輔助完成對(duì)預(yù)報(bào)方案的率定，預(yù)報(bào)期可以達(dá)到72h。該系統(tǒng)現(xiàn)已經(jīng)在實(shí)際洪水預(yù)報(bào)工作中應(yīng)用，目前運(yùn)行穩(wěn)定可靠。

對(duì)南車水庫(kù)流域確定的6場(chǎng)洪水過(guò)程進(jìn)行模擬統(tǒng)計(jì)，采用地貌單位線進(jìn)行流域匯流計(jì)算。模擬結(jié)果洪量、洪峰、峰現(xiàn)時(shí)間的平均相對(duì)誤差分別為17%、19%和3h，場(chǎng)次洪水綜合預(yù)報(bào)合格率50%。圖6～圖11 為實(shí)測(cè)與模擬洪水過(guò)程線比較結(jié)果。

圖6 20170628洪水預(yù)報(bào)過(guò)程注 圖中綠色線條代表實(shí)測(cè)流量過(guò)程，紅色為預(yù)報(bào)流量過(guò)程。

圖7 20170711洪水預(yù)報(bào)過(guò)程注 圖中綠色線條代表實(shí)測(cè)流量過(guò)程，紅色為預(yù)報(bào)流量過(guò)程。

圖8 20180606洪水預(yù)報(bào)過(guò)程注 圖中綠色線條代表實(shí)測(cè)流量過(guò)程，紅色為預(yù)報(bào)流量過(guò)程。

圖9 20180612洪水預(yù)報(bào)過(guò)程注 圖中綠色線條代表實(shí)測(cè)流量過(guò)程，紅色為預(yù)報(bào)流量過(guò)程。

圖10 20180621洪水預(yù)報(bào)過(guò)程注 圖中綠色線條代表實(shí)測(cè)流量過(guò)程，紅色為預(yù)報(bào)流量過(guò)程。

圖11 20180622洪水預(yù)報(bào)過(guò)程注 圖中綠色線條代表實(shí)測(cè)流量過(guò)程，紅色為預(yù)報(bào)流量過(guò)程。

場(chǎng)次洪水評(píng)定參照水文情報(bào)預(yù)報(bào)規(guī)范，洪峰、洪量預(yù)報(bào)以預(yù)見(jiàn)期內(nèi)實(shí)測(cè)變幅的20%作為許可誤差，峰現(xiàn)時(shí)間按3h計(jì)，預(yù)報(bào)結(jié)果見(jiàn)表3[7]。

表3 洪水預(yù)報(bào)結(jié)果統(tǒng)計(jì)

續(xù)表

4 結(jié) 語(yǔ)

本文根據(jù)童冰星等[8]提出的地貌單位線確定方法，以吉安市南車水庫(kù)流域?yàn)槔?，利用DEM信息提取了流域的地貌單位線，并將之用于該流域的洪水模擬，通過(guò)對(duì)6場(chǎng)次洪水預(yù)報(bào)與實(shí)測(cè)值比較，峰現(xiàn)時(shí)間合格率為83.3%，洪峰流量合格率為66.7%，洪水總量合格率為66.7%，場(chǎng)次洪水綜合預(yù)報(bào)合格率為50%，取得了較好的模擬結(jié)果，進(jìn)一步驗(yàn)證了通過(guò)流路長(zhǎng)度分布，確定地貌單位線的有效性，并將此地貌單位線法應(yīng)用到吉安市31個(gè)無(wú)資料水庫(kù)洪水預(yù)報(bào)工作中。

通過(guò)研究分析預(yù)報(bào)存在誤差的原因：一方面水庫(kù)庫(kù)容采用4段制上報(bào)，入庫(kù)流量采用水量平衡原理倒推入庫(kù)，人為地對(duì)實(shí)測(cè)洪峰值進(jìn)行了均化處理，預(yù)報(bào)數(shù)值計(jì)算步長(zhǎng)為1h，影響了洪水預(yù)報(bào)精度；另一方面，預(yù)報(bào)流域的基流大小對(duì)預(yù)報(bào)結(jié)果影響較大，而對(duì)無(wú)資料地區(qū)的基流獲取尚且存在困難。隨著實(shí)際應(yīng)用和實(shí)踐，對(duì)于無(wú)資料地區(qū)的水文預(yù)報(bào)問(wèn)題，有望展開(kāi)更加深入的研究并解決上述問(wèn)題。