擴(kuò)散波與馬斯京根法在柵格匯流演算中的應(yīng)用比較

姚 成，章玉霞，李致家

(河海大學(xué)水文水資源學(xué)院，江蘇南京 210098)

水文模型作為一種研究流域內(nèi)復(fù)雜水文現(xiàn)象的有效工具，在水文學(xué)研究領(lǐng)域中一直發(fā)揮著重要作用。就目前的水文模型發(fā)展趨勢(shì)而言，分布式水文模型由于其能夠更好地考慮降雨和下墊面條件的空間變異性，能夠更好地利用GIS與RS技術(shù)等空間信息模擬流域的降雨徑流響應(yīng)，已成為水文模型研究的重要方向［1-2］。現(xiàn)有的分布式水文模型研究多以產(chǎn)流計(jì)算的空間變異性為出發(fā)點(diǎn)，在模型的構(gòu)建和應(yīng)用中更多地追求產(chǎn)流量的高精度模擬，對(duì)匯流過(guò)程卻缺乏足夠的重視。匯流過(guò)程作為洪水預(yù)報(bào)的一個(gè)重要環(huán)節(jié)，其模擬精度直接影響模型對(duì)流域出口斷面流量過(guò)程的預(yù)報(bào)精度，也是流域水文過(guò)程模擬的重點(diǎn)與難點(diǎn)。為了使分布式水文模型更好地用于洪水預(yù)報(bào)，必須加強(qiáng)對(duì)匯流過(guò)程的重視，需要對(duì)分布式水文模型中的匯流過(guò)程進(jìn)行更深一步的研究［3］。

常用的流域匯流演算方法有2種，即基于水量平衡與槽蓄方程的水文學(xué)方法和基于圣維南方程組及其簡(jiǎn)化算法的水力學(xué)方法［3-4］。比較而言，水文學(xué)方法(如馬斯京根法)計(jì)算簡(jiǎn)單，資料需求相對(duì)較低，但理論上具有近似性，使用條件會(huì)受到限制。而水力學(xué)方法可以使匯流演算更具物理性，能夠有效處理匯流演算時(shí)多方面的影響，但其計(jì)算較復(fù)雜，對(duì)資料的需求也相對(duì)較高。筆者以皖南山區(qū)呈村流域?yàn)槔?，采用新安江模型進(jìn)行柵格單元的產(chǎn)流計(jì)算，分別采用基于水力學(xué)方法的擴(kuò)散波模型與基于水文學(xué)方法的馬斯京根法進(jìn)行柵格單元間的匯流演算，將所構(gòu)建的柵格型新安江模型［5-7］用于呈村流域的洪水模擬，對(duì)2種柵格匯流演算方法進(jìn)行應(yīng)用比較，分析在分布式水文模型中采用擴(kuò)散波與馬斯京根法模擬匯流過(guò)程的時(shí)效性。

1 柵格匯流演算

柵格型新安江模型以流域內(nèi)每個(gè)DEM柵格作為計(jì)算單元，并假設(shè)在柵格單元內(nèi)的降雨及下墊面特征空間分布均勻，模型只考慮各個(gè)要素在不同柵格之間的變異性。模型以蓄滿產(chǎn)流原理為基礎(chǔ)，先計(jì)算出每個(gè)柵格單元的產(chǎn)流量，再采用自由水蓄水庫(kù)結(jié)構(gòu)將其劃分為地表徑流、壤中流以及地下徑流3種水源，最后根據(jù)柵格間的匯流演算次序依次將各種水源演算至流域出口(圖1)。對(duì)于柵格間的匯流過(guò)程模擬，筆者分別采用一維擴(kuò)散波模型與馬斯京根法。

1.1 一維擴(kuò)散波模型(DW)

擴(kuò)散波雖然是動(dòng)力波(圣維南方程組)的一種簡(jiǎn)化形式，但在許多實(shí)際情況下它都適用。而且，利用擴(kuò)散波模型進(jìn)行坡面匯流與河道匯流演算都能取得足夠高的精度，在很多時(shí)候其精度與動(dòng)力波匯流演算模型的精度相當(dāng)。本研究在采用擴(kuò)散波模型進(jìn)行匯流演算時(shí)假設(shè)任意柵格單元都是由坡地和河道組成，地下徑流與壤中流都直接匯入河道中，因此柵格間的匯流就由坡面匯流及河道匯流組成。

描述坡面水流運(yùn)動(dòng)的一維擴(kuò)散波方程組為

圖1 柵格匯流演算示意圖Fig.1 Schematic representation of cell-to-cell flow routing

式中:hs——坡面水流的水深，m;us——坡面水流的平均流速，m/s;qs——單位時(shí)間內(nèi)所計(jì)算的坡面徑流深，m/s;t——時(shí)間，s;x——流徑長(zhǎng)度，m;Soh——沿出流方向的地表坡度;Sfh——沿出流方向的地表摩阻比降。

在進(jìn)行柵格間匯流演算時(shí)，式(1)需要在每個(gè)柵格單元上進(jìn)行離散，其中的連續(xù)性方程可表示為

式中:Agc——柵格單元的面積，m2;Qsin——柵格單元的地表徑流入流量，m3/s，Qsin=Qsup+Qs;Qsup——上游柵格地表入流流量，m3/s;Qs——當(dāng)前柵格的地表徑流量，m3/s;Qsout——柵格單元的地表徑流出流量，m3/s。

本研究采用基于兩步(預(yù)測(cè)步—校正步)MacCormack算法的二階顯式有限差分格式進(jìn)行式(1)、式(2)的求解［8-10］。該差分格式為顯式差分，其穩(wěn)定性條件為

式中:g——重力加速度，m/s2;Δt——計(jì)算時(shí)間步長(zhǎng)，s;Δx——柵格單元間匯流路徑的長(zhǎng)度(正向出流時(shí)Δx=dc，斜向出流時(shí) Δx =dc)，m;dc——柵格單元邊長(zhǎng)，m。

描述河道水流運(yùn)動(dòng)的一維擴(kuò)散波方程組為

式中:Ach——河道的過(guò)水?dāng)嗝婷娣e，m2;Qch——河道流量，m3/s;ql——單寬旁側(cè)入流，m2/s;hch——河道水深，m;Soc——河道坡度;Sfc——河道摩阻比降。

對(duì)于式(4)的求解，同樣采用基于MacCormack算法的顯式差分格式。

1.2 基于柵格的馬斯京根法(MK)

根據(jù)柵格間匯流演算次序，依次將地表徑流、壤中流、地下徑流以及河道水流分別演算至流域出口。以Qs為例(圖2)，a，b，c 柵格的流量分別為 Qa，Qb，Qc，出流量Qa'，Qb'，Qc'可以通過(guò)馬斯京根法計(jì)算得到:

在t時(shí)刻，若柵格d的產(chǎn)流量為 Qs，t，則其出流可表示為

圖2 基于柵格的馬斯京根法示意圖Fig.2 Schematic representation of cell-to-cell Muskingum method

其中

式中:xe——流量比重因子;ke——蓄量常數(shù);fch——地表徑流出流量匯入河道的比例。

1.3 柵格間匯流演算次序

柵格匯流演算次序即分布式水文模型在進(jìn)行匯流演算時(shí)所需用到的柵格間演算優(yōu)先次序。計(jì)算步驟如下:(a)對(duì)原始DEM矩陣進(jìn)行預(yù)處理，生成研究流域的流向矩陣與累計(jì)匯水面積矩陣。(b)初始化柵格演算次序矩陣，將研究流域內(nèi)的柵格單元全部賦值為零，在累計(jì)匯水面積矩陣的基礎(chǔ)上，對(duì)于矩陣內(nèi)屬于研究流域的且柵格值為1的柵格演算次序賦值n=1(n代表柵格演算次序)。(c)疊加n=n+1，逐個(gè)掃描n=0的柵格，假設(shè)當(dāng)前柵格為i，根據(jù)流向矩陣判斷如下:若流向i的相鄰柵格中演算次序均非零，則i的演算次序?yàn)閚;否則為零。(d)重復(fù)步驟c，直至研究流域的出口柵格演算次序已被賦值。

2 參數(shù)確定

對(duì)于柵格型新安江模型參數(shù)的確定，先根據(jù)參數(shù)的物理意義，利用地貌特征、植被覆蓋、土壤類型等下墊面空間信息，結(jié)合已有的應(yīng)用經(jīng)驗(yàn)對(duì)參數(shù)取值做出前期估計(jì)，然后再根據(jù)參數(shù)估計(jì)的應(yīng)用結(jié)果甄選出敏感參數(shù)，最后再利用實(shí)測(cè)的出口流量過(guò)程對(duì)敏感參數(shù)進(jìn)行率定和檢驗(yàn)［10-11］。

當(dāng)采用DW方法進(jìn)行柵格匯流演算時(shí)，匯流參數(shù)包括Soh，nh(地表曼寧糙率系數(shù))，Soc，nc(河道曼寧糙率系數(shù))以及α(柵格單元的河道寬度指數(shù))。而當(dāng)采用MK方法時(shí)，匯流參數(shù)包括xe與ke。由于模型在匯流演算時(shí)考慮了河道排水網(wǎng)絡(luò)的影響，因此2種匯流方法的參數(shù)還包括fch。根據(jù)已有的研究成果，任意柵格單元的Soh與Soc可以直接通過(guò)DEM分析得到，fch也是在流域數(shù)字化的基礎(chǔ)上利用面積比例法進(jìn)行計(jì)算，α則是通過(guò)分析流域內(nèi)已有的實(shí)測(cè)斷面資料獲?。?，11］。nh與nc在擴(kuò)散波匯流演算中屬于敏感參數(shù)，nh的確定可以根據(jù)柵格單元的植被類型，參考應(yīng)用經(jīng)驗(yàn)賦予對(duì)應(yīng)的估值。在此基礎(chǔ)上，本文定義了一個(gè)空間分布均勻的糙率比例因子μnh，并利用觀測(cè)資料對(duì)該因子進(jìn)行率定，即任意柵格單元的坡面匯流糙率nh'=μnhnh。nc的計(jì)算公式為［12］

式中:Ad——柵格單元的上游匯水面積(可利用累積匯水面積矩陣確定)，km2;n0——計(jì)算系數(shù)，可以由流域出口點(diǎn)實(shí)測(cè)的水位-流量關(guān)系反算得到，當(dāng)出口點(diǎn)無(wú)實(shí)測(cè)資料時(shí)需要對(duì)n0進(jìn)行率定。

對(duì)于xe與ke的確定，暫沒(méi)有考慮其空間分布不均的問(wèn)題，而是采用流域內(nèi)柵格單元統(tǒng)一賦值的方法。xe相對(duì)不敏感，按一般經(jīng)驗(yàn)定值即可。ke屬于敏感參數(shù)，反映的是不同徑流成分在柵格單元的匯流時(shí)間，ke越小表明匯流速度越快，則對(duì)應(yīng)洪水過(guò)程的峰現(xiàn)時(shí)間將有所提前，峰值也會(huì)相對(duì)增大。由于河道水流、地表徑流、壤中流與地下徑流在匯流速度上存在差異，因此在率定ke時(shí)需對(duì)其進(jìn)行區(qū)分。

3 應(yīng)用實(shí)例

3.1 研究流域概況

選擇位于皖南山區(qū)的呈村流域作為研究流域，該流域鄰近東南沿海，四季分明，氣候溫和，屬亞熱帶季風(fēng)氣候區(qū)。呈村流域控制面積為290km2，地勢(shì)南高北低，相對(duì)高差較大，平均海拔高程為583 m，流域河道平均坡度為0.95%，最大匯流路徑長(zhǎng)度為36 km。流域內(nèi)植被良好，雨量充沛，多年平均降雨量約為2 100 mm，流域降水在年內(nèi)、年際分配極不均勻，為典型的濕潤(rùn)流域。該流域植被類型主要包括常綠針葉林、落葉闊葉林、混合林、森林地、林地草原、牧草地與作物地，土壤類型主要為黏壤土。流域內(nèi)有10個(gè)雨量站，具有1986—1999年32場(chǎng)洪水過(guò)程的時(shí)段資料，其中前22場(chǎng)洪水過(guò)程資料用于模型參數(shù)率定，其余資料用于模型參數(shù)檢驗(yàn)。呈村流域水系及其雨量站分布見(jiàn)圖3。

3.2 結(jié)果分析

在應(yīng)用模型進(jìn)行呈村流域的洪水模擬時(shí)，柵格單元的分辨率為1 km×1 km(30"×30")，率定的參數(shù)主要包括μnh，n0，xe及ke，其余參數(shù)則是在已有經(jīng)驗(yàn)的基礎(chǔ)上利用下墊面信息確定。例如，當(dāng)根據(jù)植被類型確定nh時(shí)，對(duì)于常綠針葉林、落葉闊葉林、混合林及森林地，nh的取值均為0.10，而對(duì)于林地草原、牧草地與作物地，nh的取值分別為0.30，0.17 及0.35［11］。參數(shù) μnh，n0，xe的率定結(jié)果分別為1.50，0.06，0.30;河道水流、地表徑流、壤中流和地下徑流對(duì)應(yīng)的ke率定結(jié)果分別為0.06，0.10，0.60和1.00。

對(duì)于呈村流域率定期內(nèi)的22場(chǎng)洪水過(guò)程而言，DW方法模擬的徑流深相對(duì)誤差(RRE)介于 -12.4% ～15.4%之間，洪峰相對(duì)誤差(RPE)介于 -19.7% ～29.5%之間，峰現(xiàn)時(shí)差(PTE)介于-3～4 h之間，確定性系數(shù)(DC)介于0.68～0.99之間;MK方法模擬的RRE，RPE，PTE誤差范圍分別為-13.4% ～19.7%，-21.6% ～27.8%，-3～6h，DC介于0.80～0.98之間。對(duì)于檢驗(yàn)期內(nèi)的10場(chǎng)洪水過(guò)程而言，DW 方法模擬的RRE，RPE，PTE誤差范圍分別為 -12.7% ～9.6%，-18.6% ～16.0%，-10～0 h，DC介于0.93～0.97之間;MK方法模擬的 RRE，RPE，PTE誤差范圍分別為 -15.0% ～4.2%，-23.4% ～17.7%，0～9 h，DC介于0.92～0.98之間。在計(jì)算效率方面，MK方法明顯優(yōu)于采用顯示差分格式的DW方法。平均而言，在30"分辨率情況下，MK方法對(duì)于呈村流域每次洪水過(guò)程(平均約7 d)的計(jì)算時(shí)間約為1.5 s，而DW方法則需要近22 s。分辨率越高，兩者之間的差別越大。表1為采用2種匯流方法模擬呈村流域洪水過(guò)程的合格率、平均相對(duì)誤差水平以及DC均值的統(tǒng)計(jì)結(jié)果。圖4為摘錄的2次實(shí)測(cè)與模擬洪水過(guò)程線的比較。應(yīng)用結(jié)果表明，DW和MK方法對(duì)于呈村流域洪水過(guò)程的模擬精度基本相當(dāng)，均能達(dá)到甲級(jí)，采用這2種方法進(jìn)行柵格間的匯流演算均能夠取得較好的應(yīng)用效果。

圖3 呈村流域水系及雨量站分布Fig.3 Drainage network and locations of gauge stations in Chengcun watershed

表1 采用DW法與MK法的模型模擬結(jié)果誤差統(tǒng)計(jì)Table 1 Accuracy statistics of model simulation results using cell-to-cell diffusion wave and Muskingum routing methods

雖然2種匯流方法在模擬精度上較為相近，但在實(shí)際應(yīng)用時(shí)仍然存在一定的差異。DW方法參數(shù)的物理意義較明確，參數(shù)的空間分布能夠通過(guò)地貌特征、植被類型等信息獲取，能夠較好地反映下墊面特征的空間變異性，不僅能輸出流域出口斷面的水位、流量過(guò)程，而且能輸出各個(gè)計(jì)算時(shí)段對(duì)應(yīng)的水深、流速等變量的空間分布。且DW方法考慮了水流運(yùn)動(dòng)中的壓力項(xiàng)，因而對(duì)坡度平緩地區(qū)的適用性較好。此外，DW方法對(duì)于無(wú)資料地區(qū)的洪水演算也具有一定的適用性。以呈村流域?yàn)槔词沽瞀蘮h=1，即不利用觀測(cè)資料對(duì)其率定，只根據(jù)植被類型確定任意柵格單元的nh值，也能取得一定的應(yīng)用效果。DW方法的缺陷主要是對(duì)于輸入信息的需求較高，計(jì)算相對(duì)復(fù)雜，運(yùn)算效率較低。當(dāng)采用顯示差分格式對(duì)擴(kuò)散波方程進(jìn)行求解時(shí)，為了保證差分格式的穩(wěn)定，計(jì)算時(shí)間步長(zhǎng)不宜過(guò)大，柵格單元分辨率越高則該時(shí)間步長(zhǎng)越小，需要的模型運(yùn)算時(shí)間也就越長(zhǎng)。與DW方法不同的是，MK方法是由水量平衡方程與槽蓄方程推導(dǎo)而來(lái)，屬于水文學(xué)方法，對(duì)輸入信息的要求相對(duì)較低，計(jì)算簡(jiǎn)便，運(yùn)算效率較高，實(shí)用性較好。但該方法對(duì)于參數(shù)及變量空間分布的描述相對(duì)不足，不能充分考慮植被、土壤等下墊面特征對(duì)匯流過(guò)程的影響，對(duì)降雨徑流觀測(cè)資料的依賴性較高，在無(wú)資料地區(qū)的適用性仍需進(jìn)一步研究。

圖4 采用DW法與MK方法模擬的洪水過(guò)程線比較Comparison of flood hydrographs simulated by cell-to-cell diffusion wave and Muskingum routing methods

4 結(jié) 語(yǔ)

考慮到分布式水文模型自身的特點(diǎn)，以及無(wú)資料地區(qū)洪水預(yù)報(bào)工作的實(shí)際需要，采用何種方法進(jìn)行計(jì)算單元間的匯流演算已成為水文科學(xué)研究中的熱點(diǎn)與難點(diǎn)問(wèn)題之一。本文以皖南山區(qū)呈村流域?yàn)槔?，針?duì)DW與MK這2種柵格匯流演算方法開(kāi)展應(yīng)用比較，結(jié)果表明，它們均能較好地用于呈村流域的洪水模擬，能夠取得較高的模擬精度。通過(guò)對(duì)2種方法進(jìn)行對(duì)比分析，發(fā)現(xiàn)DW方法有較好的物理基礎(chǔ)，能夠利用地貌特征、植被覆蓋、土壤類型等下墊面信息對(duì)參數(shù)及其空間分布作出估計(jì)，對(duì)于無(wú)資料地區(qū)具有一定的適用性，但其運(yùn)算效率相對(duì)較低，對(duì)于大尺度流域的實(shí)用性尚需探討。而MK方法卻計(jì)算簡(jiǎn)便，實(shí)用性較好，但其對(duì)于觀測(cè)資料有較高的依賴性，不能很好地考慮參數(shù)及變量的空間變異性。如何進(jìn)一步挖掘2種方法的物理內(nèi)涵，使其更好地用于分布式水文模型，還有待更深入的研究。

［1］芮孝芳，黃國(guó)如.分布式水文模型的現(xiàn)狀與未來(lái)［J］.水利水電科技進(jìn)展，2004，24(2):55-58.(RUI Xiaofang，HUANG Guoru.The present and future of the distributed hydrological models［J］.Advances in Science and Technology of Water Resources，2004，24(2):55-58.(in Chinese))

［2］熊立華，郭生練.分布式流域水文模型［M］.北京:中國(guó)水利水電出版社，2004.

［3］芮孝芳.水文學(xué)研究進(jìn)展［M］.南京:河海大學(xué)出版社，2007.

［4］李致家，孔凡哲，王棟，等.現(xiàn)代水文模擬與預(yù)報(bào)技術(shù)［M］.南京:河海大學(xué)出版社，2010.

［5］李致家，姚成，汪中華.基于柵格的新安江模型的構(gòu)建和應(yīng)用［J］.河海大學(xué)學(xué)報(bào):自然科學(xué)版，2007，35(2):131-134.(LI Zhijia，YAO Cheng，WANG Zhonghua.Development and application of grid-based Xin’anjiang model［J］.Journal of Hohai University:Nature Sciences，2007，35(2):131-134.(in Chinese))

［6］李致家，姚成，章玉霞，等.柵格型新安江模型的研究［J］.水力發(fā)電學(xué)報(bào)，2009，28(2):25-33.(LI Zhijia，YAO Cheng，ZHANG Yuxia，et al.Study on gird-based Xin’anjiang model［J］.Journal of Hydroelectric Engineering，2009，28(2):25-33.(in Chinese))

［7］YAOCheng，LI Zhijia，BAOHongjun，et al.Application of a developed Grid-Xin’anjiang model to Chinese watersheds for flood forecasting purpose［J］.Journal of Hydrologic Engineering，2009，14(9):923-934.

［8］WANG M，HJELMFELT A.DEM based overland flow routing model［J］.Journal of Hydrologic Engineering，1998，3(1):1-8.

［9］JAIN M，SINGH V.DEM-based modelling of surface runoff using diffusion wave equation［J］.Journal of Hydrology，2005，302:107-126.

［10］YAO Cheng，LI Zhijia，YU Zhongbo，et al.A priori parameter estimates for a distributed，grid-based Xin’anjiang model using geographically based information［J］.Journal of Hydrology，2012，468/469:47-62.

［11］姚成，紀(jì)益秋，李致家，等.柵格型新安江模型的參數(shù)估計(jì)及應(yīng)用［J］.河海大學(xué)學(xué)報(bào):自然科學(xué)版，2012，40(1):42-47.(YAO Cheng，JI Yiqiu，LI Zhijia，et al.Parameter estimation and application of grid-based Xin’anjiang model［J］.Journal of Hohai University:Natural Sciences，2012，40(1):42-47.(in Chinese))

［12］KOREN V，REED S，SMITH M，et al.Hydrology laboratory research modeling system(HL-RMS)of the national weather service［J］.Journal of Hydrology，2004，291:297-318.