無資料地區(qū)水文模擬及相似性分析

姚 成，章玉霞，李致家，孫 龍

(1.河海大學(xué)水文水資源學(xué)院，江蘇南京 210098;2.水利部水文局，北京 100053)

無資料地區(qū)水文預(yù)報(bào)研究涉及水文科學(xué)的基礎(chǔ)理論和方法，是水文科學(xué)研究中方興未艾的研究熱點(diǎn)和難點(diǎn)。為了推進(jìn)無資料地區(qū)水文科學(xué)理論的發(fā)展，國際水文科學(xué)協(xié)會于2003年啟動(dòng)了PUB(predictions in ungauged basins)國際水文計(jì)劃，意在用10 a時(shí)間集中研究無資料地區(qū)的水文預(yù)報(bào)問題，探索解決該問題的新途徑和新方法，實(shí)現(xiàn)水文理論的重大突破，促進(jìn)水文科學(xué)在解決實(shí)際問題中的應(yīng)用［1－2］。PUB計(jì)劃的啟動(dòng)極大地激發(fā)了學(xué)者們開展無資料地區(qū)水文預(yù)報(bào)研究的熱情［3－7］。

水文模型作為一種研究流域內(nèi)復(fù)雜水文現(xiàn)象的有效工具，一直在PUB研究領(lǐng)域中發(fā)揮著重要作用。由于在無資料地區(qū)沒有或缺乏歷史觀測資料，不能進(jìn)行模型參數(shù)的率定與檢驗(yàn)，因此，解決參數(shù)的確定問題是水文模型在無資料地區(qū)應(yīng)用的關(guān)鍵所在。目前，常用的模型參數(shù)確定方法包括參數(shù)移植、參數(shù)估計(jì)以及參數(shù)的區(qū)域回歸方法等。隨著對PUB問題研究的深入，參數(shù)移植法在實(shí)際預(yù)報(bào)工作中的應(yīng)用日趨廣泛［8］。

筆者基于新安江模型，以位于皖南山區(qū)的嵌套式屯溪流域?yàn)槔?，采用SCE－UA算法率定屯溪流域的模型參數(shù)，并假設(shè)其內(nèi)部子流域?yàn)闊o資料流域，對模型參數(shù)進(jìn)行移植，在不經(jīng)重新率定參數(shù)的情況下客觀評價(jià)新安江模型模擬嵌套站點(diǎn)出流過程的能力，并對影響參數(shù)移植精度的主要因素以及流域間的水文相似性進(jìn)行分析，以期為進(jìn)一步提高無資料地區(qū)水文預(yù)報(bào)精度奠定基礎(chǔ)。

1 新安江模型及其參數(shù)率定

1.1 模型簡介

新安江模型作為一個(gè)概念性水文模型，自建立以來已得到大量的應(yīng)用驗(yàn)證，取得良好的驗(yàn)證效果，是目前在中國的實(shí)時(shí)洪水預(yù)報(bào)中應(yīng)用最廣泛的水文模型之一。該模型先根據(jù)降水和下墊面特征將流域劃分為若干個(gè)單元，然后對每個(gè)單元分別進(jìn)行產(chǎn)匯流計(jì)算，得到各個(gè)單元流域的出流過程，最后將其演算至流域出口并進(jìn)行疊加，即可得到整個(gè)流域的出流過程。新安江模型由4個(gè)模塊組成，包括蒸散發(fā)模塊、產(chǎn)流模塊、分水源模塊以及匯流模塊，每個(gè)模塊分別對應(yīng)不同的模型參數(shù)［9］。在進(jìn)行模型的實(shí)際應(yīng)用時(shí)，如何確定合理可靠的模型參數(shù)往往是需要解決的首要問題［10］。

1.2 參數(shù)率定方法

隨著計(jì)算機(jī)技術(shù)與現(xiàn)代數(shù)學(xué)方法的發(fā)展，新安江模型參數(shù)的率定方法正逐漸由人工試錯(cuò)法向自動(dòng)優(yōu)化方法轉(zhuǎn)變。其中，比較具有代表性的是SCE－UA全局優(yōu)化算法。該算法將基于確定性的復(fù)合形搜索技術(shù)與自然界中的生物競爭進(jìn)化原理相結(jié)合，能夠根據(jù)目標(biāo)函數(shù)快速搜索到全局最優(yōu)解，是一種可以有效解決非線性約束最優(yōu)化問題的方法［11］。筆者即是采用SCE－UA算法對新安江模型的參數(shù)進(jìn)行率定和檢驗(yàn)。

當(dāng)采用SCE－UA算法進(jìn)行水文模型參數(shù)率定時(shí)，目標(biāo)函數(shù)的選取對于率定結(jié)果有著顯著影響，不同的目標(biāo)函數(shù)往往會得到不同的優(yōu)化結(jié)果。為了使SCE－UA算法中的目標(biāo)函數(shù)能夠較全面地描述由實(shí)測資料所反映的水文特征，根據(jù) GB/T 22482—2008《水文情報(bào)預(yù)報(bào)規(guī)范》中的精度評定指標(biāo)，采用綜合誤差系數(shù)(CEC)［12－13］的方法建立目標(biāo)函數(shù)(Fun):

式中:Ei——任意一個(gè)獨(dú)立的誤差或精度評定指標(biāo);n——評定指標(biāo)的個(gè)數(shù)。

對于洪水過程而言，以徑流深相對誤差(RRE)、洪峰相對誤差(RPE)、峰現(xiàn)時(shí)差(PTE)以及確定性系數(shù)(DC)作為評定指標(biāo)，則n=4;對于日徑流過程而言，可以只采用RPE與DC作為評定指標(biāo)，則n=2。CEC介于0～1之間，CEC值越大，目標(biāo)函數(shù)Fun值越小，表明模型計(jì)算結(jié)果越好;反之，CEC值越小，F(xiàn)un值越大，表明模型計(jì)算結(jié)果越差。

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

2.1 研究流域概況

選用位于皖南山區(qū)的屯溪流域作為研究流域，該流域鄰近中國東南沿海，屬亞熱帶季風(fēng)氣候區(qū)，四季分明，氣候溫和，多年平均氣溫約為17℃。屯溪流域面積為2692 km2，地勢西高東低，相對高差較大。該流域植被良好，雨量充沛，多年平均降雨量約為1 800 mm，降水在年內(nèi)、年際分配極不均勻，汛期內(nèi)的降雨量一般占年總降雨量的60%以上，為典型的濕潤流域。除屯溪水文站外，流域內(nèi)還有6個(gè)水文站(圖1)，均具有較豐富的日徑流資料與次洪過程資料。屯溪流域及其嵌套子流域的地貌特征值及其相關(guān)資料見表1。表1中的流域地貌特征值均是基于數(shù)字高程模型(DEM)［14］提取得到的。

2.2 應(yīng)用結(jié)果分析

基于SCE－UA算法，利用屯溪流域70%的實(shí)測資料率定新安江模型參數(shù)，剩余30%的資料用于模型參數(shù)檢驗(yàn)。對于月潭、萬安等嵌套子流域，模型參數(shù)直接采用屯溪流域的率定結(jié)果，不再進(jìn)行模型參數(shù)的重新率定。圖2、圖3為屯溪流域及其嵌套子流域的率定期與檢驗(yàn)期內(nèi)日徑流模擬結(jié)果與次洪模擬結(jié)果箱形圖，該圖描述了各個(gè)精度統(tǒng)計(jì)結(jié)果的最大值、上四分位數(shù)、中位數(shù)、下四分位數(shù)及其最小值，表明新安江模型在不同流域模擬精度的變化幅度。表2統(tǒng)計(jì)了率定期與檢驗(yàn)期內(nèi)模型模擬結(jié)果的合格率、平均相對誤差水平以及DC均值的整體效果。

圖1 屯溪流域水系及其嵌套水文站點(diǎn)分布Fig.1 Drainage network of Tunxi Basin and locations of nested hydrological stations

表1 屯溪流域及其嵌套子流域概況Table 1 Characteristics of Tunxi Basin and nested sub-basins

圖2 新安江日模型模擬結(jié)果統(tǒng)計(jì)箱形圖Fig.2 Boxplots of statistics for the daily Xin’anjiang model simulations

圖3 新安江次洪模型模擬結(jié)果統(tǒng)計(jì)箱形圖Fig.3 Boxplots of statistics for hourly Xin’anjiang model simulations

由圖2、圖3及表2可知，采用SCE－UA優(yōu)化算法，新安江模型能夠很好地用于屯溪流域的日徑流與洪水過程模擬，可以取得較高的模擬精度。對于嵌套子流域而言，采用參數(shù)移植法的模型模擬徑流深與過程線的精度相對較高，對RRE與DC的控制相對較好;而模擬的洪峰與峰現(xiàn)時(shí)間精度相對較低。若以DC作為精度等級評定標(biāo)準(zhǔn)，參數(shù)移植精度等級均可達(dá)乙級以上;若以RPE合格率作為等級評定標(biāo)準(zhǔn)，只有月潭、萬安2個(gè)流域的參數(shù)移植精度可達(dá)乙級，其他4個(gè)流域的模擬精度甚至達(dá)不到丙級。比較而言，子流域面積越大，洪峰模擬精度相對越高。該結(jié)果表明，在嵌套式流域采用參數(shù)移植法進(jìn)行新安江模型實(shí)際應(yīng)用時(shí)，徑流深計(jì)算結(jié)果一般能夠滿足精度需求，獲得的洪水過程比較合理，但洪峰誤差相對較大，在實(shí)際應(yīng)用時(shí)需加以重視。

在用新安江模型進(jìn)行次洪過程模擬時(shí)，流域蓄水消退系數(shù)(CS)的取值對模擬結(jié)果影響相當(dāng)大［15］。由參數(shù)移植結(jié)果可知，即使在嵌套式流域內(nèi)進(jìn)行CS的移植，洪峰模擬精度仍然不高。分析其原因，可能主要是流域的坡度(SA)差異導(dǎo)致。一般而言，SA越小，流域的匯流速率越慢，此時(shí)受到的流域調(diào)蓄作用相對越強(qiáng)，CS取值也應(yīng)該較大;反之，SA越大，匯流越集中，流域的匯流速率越快，受到的流域調(diào)蓄作用越弱，CS取值也應(yīng)該較小。由于屯溪流域及嵌套子流域間的坡度值存在較大的差異，直接采用屯溪流域的CS模擬子流域的洪水過程，其精度勢必會受到影響。因此，在新安江次洪模型的參數(shù)移植中，有必要對地貌特征加強(qiáng)考慮。

表2 新安江次洪模型模擬結(jié)果誤差統(tǒng)計(jì)Table 2 Accuracy statistics of hourly Xin’anjiang model simulation results

3 水文相似性分析

通過以上分析可以發(fā)現(xiàn)，在嵌套式流域進(jìn)行新安江模型的參數(shù)移植，能夠取得較好的徑流深模擬結(jié)果，但洪峰模擬精度不高，表明研究流域在產(chǎn)流過程上具有一定的相似性，而流域匯流過程差異較大。一般而言，流域水文響應(yīng)是降雨和下墊面條件綜合作用的產(chǎn)物，當(dāng)降雨條件一定時(shí)，流域的地形地貌特征對流域水文響應(yīng)有著重要的影響［16］。因此，筆者以研究流域的地貌特征為基礎(chǔ)，開展流域的水文相似性分析，以探尋提高無資料地區(qū)參數(shù)移植精度的方法。

就流域產(chǎn)流而言，可將地形指數(shù)作為一個(gè)相似性判別標(biāo)準(zhǔn)，因?yàn)樗从沉送寥廊彼颗c產(chǎn)流過程的空間分布［17－18］。地形指數(shù)的數(shù)學(xué)表達(dá)為，其中α為單寬集水面積，tanβ為局部地表坡度，描述了徑流在流域中任意一點(diǎn)的累積趨勢以及在重力作用下順坡移動(dòng)的趨勢［19］。一般而言，在濕潤地區(qū)地形指數(shù)大的地方往往位于河道附近或河谷低洼處，而這些區(qū)域的地下水埋深較淺，包氣帶相對較薄，易于產(chǎn)流;相反，地形指數(shù)小的地方大多位于流域的上游山坡，遠(yuǎn)離河道，包氣帶相對較厚，不易產(chǎn)流，這與實(shí)際的產(chǎn)流現(xiàn)象基本相符［12］。因此，可以認(rèn)為流域上地形指數(shù)相同的點(diǎn)具有相同的降雨－徑流響應(yīng)，而通過流域的地形指數(shù)分布曲線即可模擬變動(dòng)的產(chǎn)流面積。對于不同的研究流域，即可認(rèn)為是地形指數(shù)分布曲線相同或相似的流域，在降雨－徑流過程上具有一定的相似性。

根據(jù)地形指數(shù)的計(jì)算方法，利用DEM信息，對屯溪流域及其嵌套子流域的地形指數(shù)進(jìn)行計(jì)算，提取各個(gè)流域的地形指數(shù)分布曲線如圖4(a)所示?？梢钥闯?，屯溪流域與其各子流域的地形指數(shù)分布曲線很相似，結(jié)合新安江模型的產(chǎn)流模擬結(jié)果(圖2(a)、圖3(a)與表2中RRE的統(tǒng)計(jì)結(jié)果)，進(jìn)一步表明，在地形指數(shù)分布曲線相似的流域進(jìn)行模型產(chǎn)流參數(shù)移植，可以取得較高的產(chǎn)流計(jì)算精度，地形指數(shù)分布曲線的相似性可作為判別新安江模型產(chǎn)流參數(shù)能否進(jìn)行移植的標(biāo)準(zhǔn)。

就流域匯流而言，它是地貌擴(kuò)散與水動(dòng)力擴(kuò)散共同作用的結(jié)果，基于統(tǒng)計(jì)物理學(xué)方法的地貌單位線理論使流域匯流同流域的地形地貌特征有機(jī)地結(jié)合起來［20］，因此也可將地貌單位線作為一個(gè)相似性判別標(biāo)準(zhǔn)。文獻(xiàn)［21］采用雨滴匯流時(shí)間等于其流路長度與其速度之商這一基本關(guān)系式，基于DEM信息，提出由流路長度分布律和坡度分布律確定地貌單位線的方法，取得了良好的驗(yàn)證效果。筆者采用該方法提取了屯溪流域與其各子流域的地貌單位線，如圖4(b)所示?？梢钥闯?，各個(gè)流域的地貌單位線存在較大的差異性，表明流域間匯流過程的相似性較小，這與前述的分析結(jié)論相吻合。

4 基于地貌單位線的參數(shù)移植分析

通過新安江模型參數(shù)移植的應(yīng)用結(jié)果以及相似性分析可以看出，若流域間地形指數(shù)分布曲線相似度較高，則可以直接移用新安江模型的產(chǎn)流參數(shù)，且能夠獲得較好的徑流深計(jì)算結(jié)果;但若在地貌單位線相似度不高的情況下直接移用模型的匯流參數(shù)，則不易獲得高精度的洪峰計(jì)算結(jié)果。新安江模型中流域匯流演算采用滯后演算法，其中敏感性參數(shù)CS受流域地貌特征的調(diào)蓄作用，在不同流域間難以直接移用。考慮到地貌單位線理論能夠更好地反映地貌特征對流域匯流的影響，筆者嘗試將該理論用于新安江模型中，即采用新安江模型進(jìn)行流域的產(chǎn)流計(jì)算，采用地貌單位線進(jìn)行流域的匯流計(jì)算，進(jìn)而分析此種情況下的參數(shù)移植精度。

在確定地貌單位線時(shí)所用的坡面流速計(jì)算公式為［19］

式中:v——坡面流速;S——坡面坡度;a——反映坡面糙率的經(jīng)驗(yàn)系數(shù)。

由于S可以根據(jù)DEM計(jì)算得到，在確定地貌單位線時(shí)需要率定的參數(shù)只有a。利用式(3)提取不同流域地貌單位線時(shí)，若流域的坡度存在差異，即使采用相同的a值，所獲得的結(jié)果亦不相同，仍然能夠在一定程度上反映流域間匯流過程的差異性。

圖4 基于DEM提取的地形指數(shù)分布曲線和1 h單位線g.4 Distribution functions of topographic index and 1-h unit hydrograph u(1，t)obtained from DEM

根據(jù)屯溪流域的次洪過程資料，通過率定得a=1.76 m/s，將該參數(shù)值移用于月潭、萬安等6個(gè)嵌套子流域，所得到的洪水模擬結(jié)果見表3。可以看出，在新安江模型中采用地貌單位線法進(jìn)行匯流計(jì)算也能取得較高的計(jì)算精度。對比表2和表3后發(fā)現(xiàn)，對于除月潭以外的5個(gè)子流域而言，移用地貌單位線參數(shù)a時(shí)的洪峰模擬精度明顯優(yōu)于移用滯后演算法參數(shù)CS時(shí)的洪峰模擬精度。該結(jié)果表明，在無資料地區(qū)洪水模擬預(yù)測中，可用地貌單位線法計(jì)算流域的匯流過程，以便更好地考慮地貌特征的影響，這將有助于提高無資料地區(qū)的參數(shù)移植精度。

表3 基于地貌單位線的新安江次洪模型模擬結(jié)果誤差統(tǒng)計(jì)Table 3 Accuracy statistics of hourly Xin’anjiang model simulation results based on geomorphologic unit hydrograph

5 結(jié) 語

a.采用基于綜合誤差系數(shù)的SCE－UA優(yōu)化算法，新安江模型能夠很好地用于屯溪流域的日徑流與洪水過程模擬，可以取得甲級模擬精度。

b.對于嵌套的6個(gè)子流域而言，進(jìn)行新安江模型的參數(shù)移植，徑流深模擬結(jié)果較好，獲得的出流過程比較合理，但洪峰誤差和峰現(xiàn)時(shí)差較大。在次洪模擬時(shí)應(yīng)對地貌特征加強(qiáng)考慮，建立模型匯流參數(shù)與地貌特征的定量關(guān)系，尋求參數(shù)的地貌規(guī)律，降低其對觀測資料的依賴性。

c.地形指數(shù)分布曲線的相似性可以作為判別新安江模型產(chǎn)流參數(shù)能否進(jìn)行移植的標(biāo)準(zhǔn)。若流域間地形指數(shù)分布曲線相似度較高，可以考慮直接移用模型產(chǎn)流參數(shù)。對于地貌單位線相似性較小的流域，不建議直接移用模型的滯后演算法參數(shù)。在進(jìn)行無資料地區(qū)匯流過程模擬時(shí)，采用地貌單位線法有助于提高參數(shù)移植的精度與適應(yīng)性。地貌單位線的相似性能否直接作為水文模型匯流參數(shù)移植的判別標(biāo)準(zhǔn)尚需進(jìn)一步驗(yàn)證。同時(shí)，在進(jìn)行無資料地區(qū)參數(shù)移植時(shí)還應(yīng)該進(jìn)一步加強(qiáng)對土壤、植被等因素的考慮。

［1］SIVAPALAN M，TAKEUCHI K，F(xiàn)RANKS S，et al.IAHS decade on predictions in ungauged basins(PUB)，2003－2012:shaping an exciting future for the hydrological sciences［J］.Hydrological Sciences Journal，2003，48(6):857－880.

［2］賈仰文，王浩，倪廣恒，等.分布式流域水文模型原理與實(shí)踐［M］.北京:中國水利水電出版社，2005.

［3］劉蘇峽，劉昌明，趙衛(wèi)民.無測站流域水文預(yù)測(PUB)的研究方法［J］.地理科學(xué)進(jìn)展，2010，29(11):1333－1339.(LIU Suxia，LIU Changming，ZHAO Weimin.Towards the methodology for predictions in ungauged basins［J］.Process in Geography，2010，29(11):1333－1339.(in Chinese))

［4］BáRDOSSY A，SINGH S.Robust estimation of hydrological model parameters［J］.Hydrology and Earth System Sciences，2008，5:1641－1675.

［5］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.

［6］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.

［7］姚成，紀(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))

［8］WAGENER T，WHEATER H，GUPTA V.Rainfall－runoff modelling in gauged and ungauged catchments［M］.London:Imperial College Press，2004.

［9］趙人俊.流域水文模擬［M］.北京:水利電力出版社，1984.

［10］芮孝芳，凌哲，劉寧寧，等.新安江模型的起源及對其進(jìn)一步發(fā)展的建議［J］.水利水電科技進(jìn)展，2012，32(4):1－5.(RUI Xiaofang，LING Zhe，LIU Ningning，et al.Origin of Xin'anjiang model and its further development［J］.Advances in Science and Technology of Water Resources，2012，32(4):1－5.(in Chinese))

［11］DUAN Q，SOROOSHIAN S，GUPTA V.Effective and efficient global optimization for conceptual rainfall－runoff models［J］.Water Resource Research，1992，28(4):1015－1031.

［12］姚成.基于柵格的新安江(Grid－Xin'anjiang)模型研究［D］.南京:河海大學(xué)，2009.

［13］葉金印，姚成，李京兵，等.綜合誤差系數(shù)在新安江模型參數(shù)全局優(yōu)化中的應(yīng)用［J］.河海大學(xué)學(xué)報(bào):自然科學(xué)版，2013，41(1):1－5.(YE Jinyin，YAOCheng，LI Jingbing，et al.Global parameter optimization for the Xin'anjiang model based on the comprehensive error coefficient［J］.Journal of Hohai University:Natural Sciences，2013，41(1):1－5.(in Chinese))

［14］中國科學(xué)院計(jì)算機(jī)網(wǎng)絡(luò)信息中心國際科學(xué)數(shù)據(jù)服務(wù)平臺.SRTM 90米分辨率數(shù)字高程原始數(shù)據(jù)［EB/OL］.［2009］.http://datamirror.csdb.cn.

［15］ZHAO Renjun，LIU Xinren.The Xin'anjiang model［C］//SINGH V.Computer Models of Watershed Hydrology.Colorado:Water Resources Publications，1995:215－232.

［16］李致家，于莎莎，李巧玲，等.降雨－徑流關(guān)系的區(qū)域規(guī)律［J］.河海大學(xué)學(xué)報(bào):自然科學(xué)版，2012，40(6):597－604.(LI Zhijia，YU Shasha，LI Qiaoling，et al.Regional pattern of rainfall－runoff relationship［J］.Journal of Hohai University:Natural Sciences，2012，40(6):597－604.(in Chinese))

［17］BEVEN K，KIRKBY M.A physically based variable contributing area model of basin hydrology［J］.Hydrological Sciences Bulletin，1979，24(1):43－69.

［18］孔凡哲，芮孝芳.基于地形特征的流域水文相似性［J］.地理研究，2003，22(6):709－715.(KONG Fanzhe，RUI Xiaofang.Hydrological similarity of catchments based on topography［J］.Geographical Research，2003，22(6):709－715.(in Chinese))

［19］QUINN P，BEVEN K，LAMB R.The ln(α/tanβ)index:how to calculate it and how to use it in the TOPMODEL framework［J］.Hydrological Processes，1994，9:161－185.

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

［21］芮孝芳.由流路長度分布律和坡度分布律確定地貌單位線［J］.水科學(xué)進(jìn)展，2003，14(5):602－606.(RUI Xiaofang.Study of determining geomorphologic unit hydrograph by means of probability density functions of path length and slope［J］.Advances in Water Science，2003，14(5):602－606.(in Chinese))