基于R-V GIUH的臨界雨量計算研究

王新宏，吳 巍，唐永鵬,2，龔立堯，寧利中

(1西安理工大學 省部共建西北旱區(qū)生態(tài)水利國家重點實驗室，陜西 西安 710048; 2水利部新疆維吾爾自治區(qū)水利水電勘測設計研究院，新疆 烏魯木齊 830000)

臨界雨量是山洪災害預警預報和防治體系的重要參數(shù)，臨界雨量的計算精度對于減少人員傷亡和財產(chǎn)損失有非常關(guān)鍵的作用[1]。如果計算的臨界雨量偏大，可能導致災害已經(jīng)發(fā)生而人員還未轉(zhuǎn)移或來不及轉(zhuǎn)移；反之，又會導致人員過早轉(zhuǎn)移，使得預警時間過長而失去預警的意義。

臨界雨量的影響因素主要有3個，包括降雨、土壤含水量及下墊面。在計算中，不可能孤立的只考慮某一個因素，否則既不符合實際情況也不滿足精度要求。近年來，國內(nèi)外學者在臨界雨量的主要研究中，均是在同時考慮3種因素的基礎上建立相關(guān)的數(shù)學模型來計算臨界雨量指標。王鑫等[2]基于二維淺水方程，建立了暴雨山洪的水動力學模型，對流域的臨界雨量計算進行了詳細研究。李昌志等[3]利用分布式水文模型,采用試算的辦法對臨界雨量進行確定，并結(jié)合土壤含水量變化，相應計算了臨界雨量的動態(tài)變化過程。葉金印等[4]以新安江模型為基礎，提出在考慮土壤含水量飽和度的基礎上對山洪災害進行動態(tài)預警。陳桂亞等[5]利用實測雨量統(tǒng)計法，借助大量實測資料統(tǒng)計了區(qū)域臨界雨量初值，進而確定出區(qū)域臨界雨量的變幅，以此給出預警對象的臨界雨量。江錦紅等[6]著眼于河道的安全泄洪流量，提出了最小臨界雨量和臨界雨力的新概念和計算方法，并以此構(gòu)建暴雨臨界曲線作為山洪災害預警標準。陳文輝[7]利用比擬法，根據(jù)臨界雨量的計算初值,結(jié)合內(nèi)插法和山洪災害實例調(diào)查法綜合確定流域的臨界雨量。

臨界雨量計算的另一個重要環(huán)節(jié)是匯流計算，而流域的匯流過程往往比較復雜，涉及的因素也較多，因此國內(nèi)外學者從不同角度對匯流進行了相關(guān)研究，隨著理論的創(chuàng)新與技術(shù)的發(fā)展，不斷建立了相關(guān)的匯流模型來模擬流域匯流。王建金等[8]利用與馬斯京根匯流模型耦合的BP神經(jīng)網(wǎng)絡來計算匯流，有效提高了山洪災害預警預報的精度。岳麗麗等[9]將時變線性匯流模型應用于半濕潤半干旱地區(qū)，很好模擬了流域的洪水過程。姚蕾等[10]將Nash匯流模型應用到無資料地區(qū)的匯流計算中，表現(xiàn)出很好的適用性。

上述幾種臨界雨量的計算方法各具特點，但應用于不同地區(qū)時也有很大的局限性，比如不同方法適用的流域面積范圍各不相同，并且有的方法對資料的要求較高，模型的建立必須基于詳細的水文和地形地貌資料，無法與我國水文資料匱乏的基本情況相結(jié)合。此外，有些方法(如實測雨量統(tǒng)計法)缺乏相應的理論基礎，僅通過統(tǒng)計實測資料獲得臨界雨量，可靠性無法得到保障。鑒于此，本研究以陳家河流域為例，在不需要詳細水文資料的基礎上，考慮下墊面因素建立R-V GIUH(Geomorphologic instantaneous unit hydrograph)匯流模型，利用水位流量反推法[11-12]計算臨界雨量，以期為提高臨界雨量的計算精度提供支持。

1 研究區(qū)域概況

圖1 隴縣陳家河流域的地理位置Fig.1 Geographical location of Chenjia River watershed in Longxian

2 R-V GIUH匯流模型的建立

2.1 模型理論

地貌瞬時單位線由委內(nèi)瑞拉的水文學者羅德里古茲-伊特布與瓦爾德斯在1979年提出[13]，他們經(jīng)過長期的研究認為，將流域的河網(wǎng)分級后，雨滴降落在流域內(nèi)每級河網(wǎng)都要服從一定的概率，降落在每級河網(wǎng)上的雨滴匯集的過程也服從相應的概率分布。因此從這個角度將數(shù)學中的概率理論運用到流域的匯流過程中。雨滴降落在流域之后不可能很快就會到達流域的出口斷面，而是會在流域上滯留一段時間，而且每個雨滴在流域上滯留的時間都是不同的，但是雨滴的滯留時間服從相同的概率分布函數(shù)，基于上述理論建立了地貌瞬時單位線模型。流域的地貌瞬時單位線是在一定的匯流路徑上建立的，因此必須根據(jù)陳家河流域數(shù)字高程模型(DEM)的數(shù)據(jù)，提取河網(wǎng)后得到水流路徑，進而計算出流域的地貌參數(shù)，便可以建立模型。下面主要介紹模型的表達式以及相關(guān)參數(shù)的計算方法與過程。

2.2 河網(wǎng)的提取

運用ArcGIS10.0軟件對包含陳家河流域的DEM進行處理，得到陳家河流域DEM如圖2所示。經(jīng)過洼地填充[14]、流向計算以及匯流累積量計算等水文分析后，得到用于河網(wǎng)提取的數(shù)據(jù)。閾值的確定對于河網(wǎng)提取的合理性有著決定性作用[15]。因此，對于閾值的優(yōu)化確定都需要進行深入的研究和分析。本研究按照流域的自然發(fā)育規(guī)律，結(jié)合流域臨界支撐面積的原理[14]來確定流域的面積閾值，具體的方法如下：首先選取不同的面積閾值，然后在每一個面積閾值下計算出河網(wǎng)總長度、河網(wǎng)平均比降，然后根據(jù)計算結(jié)果繪制面積閾值～河網(wǎng)總長度關(guān)系曲線和面積閾值～河網(wǎng)平均比降關(guān)系曲線，流域的臨界面積閾值一般出現(xiàn)在上述這2種關(guān)系曲線的轉(zhuǎn)折點處[16-18]。因此，將上述兩曲線擬合為函數(shù)關(guān)系并進行求導，得到其變化率曲線關(guān)系，在變化率趨于0處即為臨界支撐面積。根據(jù)上述方法得到支撐陳家河流域的臨界面積即為0.65 km2，在此基礎上提取河網(wǎng)如圖3所示。

圖2 陳家河流域DEM圖Fig.2 DEM map of the Chenjia River watershed

2.3 地貌參數(shù)的計算

地貌瞬時單位線中地貌參數(shù)的提取包括河數(shù)率、河長率和面積率。要得到地貌參數(shù)，首先需對河網(wǎng)進行分級運算[19](斯特拉勒分級法)，從而獲得包括不同級別河流的河網(wǎng)(圖4)，統(tǒng)計出不同級別河流數(shù)量及不同級別河流總長度。對陳家河流域進行子流域劃分(圖5)，得到每個子流域的面積，根據(jù)霍頓河系定律計算得到各個地貌參數(shù)見表1。表1中河數(shù)率、河長率和面積率的計算公式如下：

(1)

圖4 陳家河流域河流分級圖 Fig.4 Classification of the extracted channelnetwork in Chenjia River watershed

河流級別Stream order河流數(shù)目Stream number總河長/kmRiver length平均河長/kmAverage length總集水面積/km2Catchment area平均面積/km2Average area河數(shù)率Bifurcation ratio河長率Length ratio面積率Area ratio15654.890.9869.131.2321123.212.1116.231.483220.9810.4921.1810.594.202.663.11418.998.9910.0110.01總計 Total70108.0622.56116.5423.31

2.4 模型的建立與驗證

R-V GIUH理論的匯流模型是基于地貌參數(shù)和流域水動力參數(shù)而建立的。本研究利用文康等[20]的通用公式來計算流域的地貌瞬時單位線，即：

(2)

式中：i、j為河流分級數(shù)(i≠j)；λi、λj為第i、j級河流的平均滯留時間；Ω為河網(wǎng)級別，本研究經(jīng)計算為4；θi,Ω(0)為河網(wǎng)級別為Ω的第i級河流初始狀態(tài)概率；Aij是與平均滯留時間λi和狀態(tài)轉(zhuǎn)移概率Pij相關(guān)的系數(shù)，文中狀態(tài)轉(zhuǎn)移概率Pij為i級河流中流入j級河流的河流數(shù)占i級河流總數(shù)的比率；t為時間。

流域水動力參數(shù)包括流域面積AΩ、河道比降S0、曼寧粗糙系數(shù)n和平均河寬B。流域面積AΩ已經(jīng)在河網(wǎng)提取中獲得，為116.54 km2，面積閾值為0.65 km2，與此對應的河道比降為30.42‰；曼寧粗糙系數(shù)n的取值為0.036。河道平均寬度B往往隨著凈雨強度的變化而變化，其確定具有相當?shù)母呕?。在本研究中首先根?jù)ArcGIS10.0計算的結(jié)果初步取為1.328 m，在此基礎上初步建立匯流模型，然后通過模擬3場實測的洪水過程對B進行率定，取不同的平均河寬值B，使得模擬的時段單位線不斷與實測的洪水過程線接近，直到最接近時的平均河寬值即為凈雨強度ir對應的平均河寬值，得到3組數(shù)值后，建立B關(guān)于ir的關(guān)系曲線，并對曲線進行擬合，可以得到其關(guān)系式為B=0.407lnir+1.006。

研究選擇的菌種均為實驗室常規(guī)典型菌種，與引起淡水魚腐敗變質(zhì)的菌種（假單胞菌屬、無色桿菌屬、黃桿菌屬、摩氏桿菌屬等[4]）存在一定差異。在研究保鮮劑的抑菌效果時，不同保鮮劑對不同菌種的抑殺菌效果會有所不同，但研究方法應該是一樣的。試驗中1.30 g/100 mL檸檬酸抑菌效果最好，作為一種有機酸，不僅改變環(huán)境pH值，影響微生物抗熱性，還賦予鮮切鯽魚一定的酸味，是一種理想的水產(chǎn)品保鮮劑，但到底是單一使用的效果好還是與其他生物防腐劑復配效果好還需進一步探討。

但是在基于地貌瞬時單位線理論的匯流模型中B關(guān)于ir的關(guān)系并未經(jīng)過驗證，因此利用另外2場實測洪水過程對模型進行驗證。模型驗證既是對陳家河流域地貌瞬時單位線匯流模型是否符合當?shù)貙嶋H情況的檢驗，也是對于臨界雨量計算精度的保證。利用B關(guān)于ir的關(guān)系，計算2場實測洪水過程的平均河寬B，進而得到其地貌瞬時單位線，然后利用實測過程對2場洪水進行驗證，結(jié)果如圖6所示。從圖6可以看到，洪水過程的模擬值和實測值在峰值上擬合較好，而且趨勢上也與實測洪水吻合很好，洪峰出現(xiàn)的時差和洪峰流量的誤差也都比較小，平均誤差僅為-4.0%(表2)。在漲洪階段，模擬洪水過程和實測值基本吻合，誤差很小，而在退洪階段，發(fā)現(xiàn)實測洪水過程退洪很快，模擬值較實測值顯得相對滯后一些，有一定誤差。綜合來看，基于地形地貌的地貌瞬時單位線匯流模型精度較高。臨界雨量和準備轉(zhuǎn)移雨量的計算是基于洪峰流量和漲洪階段。因此，可以利用本研究建立的陳家河流域地貌瞬時單位線匯流模型來模擬洪水過程及計算預警指標。

表2 洪水計算值與實測值的比較Table 2 Comparison between calculated and measured values of flood

圖6 洪水過程的模擬值和實測值比較Fig.6 Comparison of calculated and measured data of flood process

2.5 匯流計算

從模型驗證的結(jié)果來看，可以利用模型對陳家河流域的匯流進行計算。因此，為了計算預警對象的臨界雨量，利用驗證后的模型對每個典型頻率的設計降雨進行匯流計算，結(jié)果見表3。

表3 陳家河流域洪峰流量計算結(jié)果Table 3 Calculated results of flood peak flow in Chenjia River watershed

將地貌瞬時單位線計算的洪水過程線結(jié)果與文獻[11]中傳統(tǒng)匯流計算方法——推理公式法[21]計算的洪水過程線結(jié)果相比較(以百年一遇的洪水過程線為例)，結(jié)果見圖7。由圖7可以發(fā)現(xiàn)，相對于地貌瞬時單位線法，推理公式法的洪峰流量偏大，并且洪峰出現(xiàn)時間偏差也較大，相對滯后。因此在預警時，會出現(xiàn)來不及轉(zhuǎn)移的情況，造成險情。推理公式法計算的洪水過程線采用五點概化，在漲洪和退洪階段都有一定的差距。從地貌瞬時單位線的模型驗證結(jié)果來看，在計算的過程中考慮了下墊面因素，精度相對較高，從而保證了臨界雨量的計算精度。

圖7 兩種方法計算洪水過程線結(jié)果的對比Fig.7 Comparison of two flood hydrographs

3 臨界雨量的計算

本研究中，臨界雨量的計算方法用水位流量反推法。該方法假定降雨和洪水頻率相同，因此只需要預警對象陳家河流域所在河道控制斷面的地形資料，就可計算預警對象的臨界雨量。對控制斷面進行水力計算，得到該流域的成災水位和成災流量分別為1 137.2 m和364 m3/s。依據(jù)陳家河流域不同洪水頻率的洪峰流量計算結(jié)果(表3)繪制洪水頻率-洪峰流量(P-Qmp)關(guān)系曲線(圖8)，然后在該曲線上由成災流量查出成災頻率為1.75%(相當于57年一遇)，推算出該成災頻率下不同降雨歷時的降雨量即為臨界雨量。表4中給出了土壤含水量(Pa)為20，50和80 mm時的臨界雨量計算結(jié)果。

圖8 洪水頻率與洪峰流量關(guān)系曲線Fig.8 Relationship of flood frequency and peak flow

mm

準備轉(zhuǎn)移雨量[22]是指在山洪發(fā)生但尚未成災時，在合理的時間內(nèi)組織人員轉(zhuǎn)移時的雨量。準備轉(zhuǎn)移雨量的計算是基于洪水過程線，而且與洪峰流量密切相關(guān)。根據(jù)陳家河流域的匯流情況，組織人員轉(zhuǎn)移的時間取為30 min，準備轉(zhuǎn)移雨量的計算過程主要包括以下步驟：(1)首先將流域的成災流量與不同頻率的洪峰流量值做一個對比，并找出大于成災流量但最接近的洪峰流量與相應的洪水頻率；(2)在該洪水頻率的洪水過程線上找到漲洪階段與洪峰流量相差30 min的流量,即為準備轉(zhuǎn)移雨量；(3)根據(jù)準備轉(zhuǎn)移流量，同樣的運用計算臨界雨量的方法便可得到不同時段的準備轉(zhuǎn)移雨量；(4)對準備轉(zhuǎn)移雨量值進行綜合分析，一般準備轉(zhuǎn)移雨量和臨界雨量的頻率關(guān)系以2～3倍為宜，如果頻率差別較大，則需要調(diào)整轉(zhuǎn)移準備時間再次計算[23]。

基于臨界雨量確定準備轉(zhuǎn)移雨量時，在洪水過程線上，按照成災水位流量出現(xiàn)前30 min左右對應的流量，反算相應的時段雨量，將該雨量作為準備轉(zhuǎn)移雨量。

根據(jù)以上方法，便可得到陳家河流域的準備轉(zhuǎn)移流量為1%的洪水過程線上與洪峰流量相差30 min的流量，在P-Qmp關(guān)系曲線上推求出準備轉(zhuǎn)移頻率[24]，設計出該頻率下不同降雨歷時的雨量即為準備轉(zhuǎn)移雨量，結(jié)果見表5。

表5 陳家河流域不同土壤含水量時的準備轉(zhuǎn)移雨量Table 5 Preparing transfer rainfall of Chenjia River watershed at different soil moisture contents mm

將基于R-V GIUH的臨界雨量計算結(jié)果與文獻[11]中基于傳統(tǒng)匯流計算方法，即推理公式法計算的臨界雨量結(jié)果相比較，結(jié)果見圖9。由圖9可以看出：基于推理公式法的臨界雨量計算結(jié)果相對于R-V GIUH的臨界雨量計算結(jié)果偏小，同一時段臨界雨量值小5～8 mm，主要是由于推理公式法計算的洪峰流量偏大，這與許多設計人員認為的推理公式法在計算結(jié)果上偏于保守相符合，同時從側(cè)面反映了基于R-V GIUH的臨界雨量計算結(jié)果是合理的，而且這種方法是建立在小流域DEM上，因此其更符合小流域自身情況，也更加準確。

準備轉(zhuǎn)移雨量的計算基于洪水過程線，相對于推理公式法的洪水過程只用五點概化，地貌瞬時單位線具有更為詳細且準確的洪水過程。因此，在準備轉(zhuǎn)移雨量的計算過程中推理公式法會將誤差進一步放大。從圖10可見，基于R-V GIUH和推理公式法計算的同一時段準備轉(zhuǎn)移雨量相差4～10 mm。

從總體來看，由于地貌瞬時單位線法匯流模型能夠很好模擬陳家河流域的洪水過程，所以臨界雨量的計算結(jié)果較為準確；而推理公式法計算的臨界雨量偏于保守，在預警過程中勢必會導致人員過早轉(zhuǎn)移，失去預警的意義。

圖9 基于推理公式法和R-V GIUH的臨界雨量計算結(jié)果對比Fig.9 Comparison of calculated critical rainfalls based on formula derivation and the R-V GIUH model

4 結(jié) 論

本研究建立了基于地形地貌的地貌瞬時單位線匯流模型，并對模型進行參數(shù)率定，使之成為符合當?shù)貙嶋H情況的匯流模型，以此來提高匯流計算的精度，從而提高臨界雨量的計算精度。本研究詳細描述了地貌參數(shù)和水動力參數(shù)的獲取方法，并利用提取的上述參數(shù)，根據(jù)R-V GIUH基本理論，建立了陳家河流域的地貌瞬時單位線匯流模型。為了從數(shù)據(jù)上和洪水過程線上反映地貌瞬時單位線匯流模型相對于傳統(tǒng)方法——推理公式法有更好的精度和適用性，將兩種匯流方式下的臨界雨量結(jié)果進行了對比,具體的結(jié)論如下。

1)運用ArcGIS10.0提取陳家河流域的DEM數(shù)據(jù)，對其進行填洼處理、流向計算、匯流累積量計算等水文分析計算，對流域的臨界支撐面積進行計算，得到面積閾值為0.65 km2的河網(wǎng)符合陳家河流域地貌的發(fā)育規(guī)律，避免了河網(wǎng)提取的隨意性。對提取的河網(wǎng)進行斯特拉勒分級法分級運算，計算得到霍頓(Horton)河數(shù)率為4.20，河長率為2.66，面積率為3.11，均符合自然流域河數(shù)率為3～5，河長率為2.5～3.5，面積率為3～6的范圍要求。

2)基于陳家河流域DEM，并利用地貌參數(shù)、水動力參數(shù)，建立了陳家河流域的地貌瞬時單位線匯流模型，并利用5場洪水對模型進行了參數(shù)率定和驗證，結(jié)果表明地貌瞬時單位線匯流模型能夠較好地模擬陳家河流域的洪水過程，并且精度也比較高，平均誤差為-4.0%。在此匯流模型下，利用設計的各個典型頻率的平均凈雨強度，計算得到1%，2%，5%，10%，20%頻率對應的洪峰流量為436.7，340.3，233.0，163.2和99.5 m3/s。

3)根據(jù)臨界雨量的計算結(jié)果，當土壤含水量Pa=50 mm，在降雨歷時1，3，6，12，24 h內(nèi)降雨量達到44，78，91，110，134 mm時應立刻組織人員轉(zhuǎn)移，否則將出現(xiàn)險情。在預警時，流域不同降雨歷時臨界雨量的大小與土壤含水量關(guān)系密切，在預警過程中要視土壤含水量的大小來確定臨界雨量和準備轉(zhuǎn)移雨量的大小。