盛奕華，衣學(xué)軍，黃迎春，李京兵，姚 成，李 昊

(1.河海大學(xué)水文水資源學(xué)院，南京 210098；2.煙臺(tái)市水文局，山東 煙臺(tái) 264000；3.安徽省水文局，合肥 230022；4.煙臺(tái)市城市水源工程運(yùn)行維護(hù)中心，山東 煙臺(tái) 264000)

流域水文模型作為一種研究流域內(nèi)復(fù)雜水文現(xiàn)象的有效工具，在洪水預(yù)報(bào)、水資源規(guī)劃和管理等工作中發(fā)揮著重要作用[1]。水文模型的參數(shù)主要依賴于水文氣象觀測(cè)資料，通過率定的方式來確定，因此無資料地區(qū)水文模型參數(shù)的確定問題是模型應(yīng)用的主要難點(diǎn)。2003年國(guó)際水文科學(xué)協(xié)會(huì)(IAHS)就此提出了無資料流域水文預(yù)測(cè)的PUB計(jì)劃(Prediction in Ungauged Basins)[2-4]，該計(jì)劃通過10多年的研究取得了豐碩成果，形成了無資料地區(qū)水文模擬和洪水預(yù)報(bào)的方法集。目前國(guó)內(nèi)外常用的無資料地區(qū)參數(shù)移植方法主要基于區(qū)域化方法，常用的方法有：多元回歸法[5]、空間臨近法[5-7]、流域物理特性相似法[5,8,9]以及全局平均法[5,10]。

水文降雨-徑流過程容易受到流域氣象條件的影響，研究表明同一流域在干旱、濕潤(rùn)等不同氣象條件下，水文響應(yīng)可能存在一定的差異性。水文模型參數(shù)的率定依賴于所選取的歷史數(shù)據(jù)，若率定期的數(shù)據(jù)代表性不足，則容易出現(xiàn)所率定的模型參數(shù)與率定期特殊的氣象條件[11]“過度擬合”的現(xiàn)象。這種氣象條件的不確定性增加了參數(shù)的不確定性[12]，使得在無資料地區(qū)模型參數(shù)移植的穩(wěn)定性不高。目前的區(qū)域化方法若只考慮兩個(gè)流域之間的參數(shù)移植，難以有效解決由于氣象因素導(dǎo)致的水文模型“過度擬合”問題，而全局平均法[13]雖然考慮了多個(gè)流域，但該方法只是將模型參數(shù)的平均值直接進(jìn)行移植，取得的模擬效果較差。

基于以上問題，本文開展基于多目標(biāo)優(yōu)化函數(shù)[14-16]的多流域HBV模型同步率定方法的研究，以降低上述參數(shù)不確定性的影響，通過獲取多個(gè)流域能夠同時(shí)適用的共享參數(shù)，提高共享參數(shù)在無資料地區(qū)的可移植性。本文以概念性水文模型HBV模型[17]為例，對(duì)安徽屯溪流域進(jìn)行模型率定和參數(shù)移植，通過與單個(gè)流域之間的模型參數(shù)移植效果進(jìn)行比較，分析基于多流域同步率定方法的參數(shù)移植效果。

1 研究區(qū)域與數(shù)據(jù)

選取安徽省南部山區(qū)的屯溪流域作為本文的研究區(qū)，該流域臨近中國(guó)東南沿海，居新安江上游，處于亞熱帶季風(fēng)氣候區(qū)，季節(jié)分明，氣候溫和，多年平均氣溫約17 ℃。其多年平均降雨量約為1 800 mm，降雨主要集中在6至9月的汛期。流域面積為2 692 km2，地勢(shì)西北高東南低，在屯溪水文站上游還嵌套有6個(gè)水文站，分別是榆村水文站、流口水文站、新亭水文站、呈村水文站、萬安水文站和月潭水文站(圖1)。

選取屯溪流域及其嵌套子流域的徑流資料用于HBV模型參數(shù)的率定和驗(yàn)證，各流域的地貌特征及資料系列見表1。用于日徑流模擬的資料序列較長(zhǎng)，其中豐、枯水年份兼具，具有一定的代表性。表中流域地貌特征等數(shù)據(jù)是基于數(shù)字高程模型(DEM)提取得到的。

表1 屯溪流域及其嵌套子流域概況Tab.1 Characteristics of Tunxi catchment and the nested sub-catchment

2 模型與方法

2.1 HBV模型

HBV(Hydrologiska By?rns Vattenbalansavdelning)模型是由瑞典氣象水文研究所(SMHI)開發(fā)的概念性水文模型，該模型由于結(jié)構(gòu)簡(jiǎn)單，需要通過歷史資料率定的參數(shù)較少且應(yīng)用簡(jiǎn)便，在瑞典、美國(guó)等國(guó)家得到了廣泛應(yīng)用，取得了良好的應(yīng)用效果[17]。該模型依據(jù)水量平衡原理進(jìn)行建模，首先通過度日因子法計(jì)算融積雪量，并通過實(shí)際土壤蓄水量函數(shù)進(jìn)行地下水補(bǔ)給量和流域?qū)嶋H蒸散發(fā)量的計(jì)算，通過線性水庫(kù)方程將徑流分為地表徑流、壤中流和地下徑流3種，最后基于單位線方法的三角權(quán)重函數(shù)進(jìn)行出口流量過程的演算[18,19]。HBV模型結(jié)構(gòu)主要包含融雪計(jì)算、土壤水分計(jì)算和流域產(chǎn)匯流3個(gè)部分[20]。由于本文研究的區(qū)域位于亞熱帶地區(qū)，降雪比重低，因此在模型中不考慮融雪計(jì)算部分。

本研究所構(gòu)建的HBV模型需要率定的參數(shù)有7個(gè)：FC、Beta、L、K0、K1、KD和K2，各個(gè)參數(shù)的物理意義[21,22]及取值范圍如表2所示。

表2 參數(shù)取值范圍Tab.2 Description and threshold values of the calibrated parameters

2.2 基于ROPE方法的模型參數(shù)優(yōu)選

ROPE(Robust Parameter Estimation algorithm)[14,15]是由Bardossy和Singh提出的一種高效穩(wěn)健的模型參數(shù)自動(dòng)優(yōu)化算法。該方法基于深度函數(shù)的概念，把模型參數(shù)組作為一個(gè)多維度空間，通過蒙特卡洛采樣不斷迭代生成高維空間內(nèi)部的新參數(shù)組，從而獲得預(yù)定數(shù)量的、具有穩(wěn)定模擬效果的參數(shù)組。ROPE算法采用的是半空間深度函數(shù)判斷新參數(shù)組的深度值。半空間深度函數(shù)可計(jì)算數(shù)據(jù)點(diǎn)相對(duì)于數(shù)據(jù)集的半空間深度(HD)，其定義如下：

設(shè)x是d維空間中的一點(diǎn)，F(xiàn)是d維空間上的一個(gè)分布函數(shù)，那么x相對(duì)于F的半空間深度就是包含x的半空間上的最小概率。

HD(x,F)=inf{F(H)，x∈H}，x∈Rd

(1)

式中：H是包含x的封閉半空間；inf為取下界。

根據(jù)半空間深度函數(shù)的性質(zhì)，深度值大于1的點(diǎn)是位于d維空間的分布函數(shù)內(nèi)部的，因此通過隨機(jī)采樣和深度值計(jì)算，可以生成d維空間內(nèi)部的新點(diǎn)距。

ROPE方法的主要步驟如下：

(1)選取需要進(jìn)行率定的模型參數(shù)，根據(jù)參考文獻(xiàn)等資料[2,14-16]，確定模型參數(shù)的取值范圍；

(2)根據(jù)參數(shù)取值范圍，通過蒙特卡洛隨機(jī)采樣方法生成n組參數(shù)，記為Xn(本研究中n=10 000)；

(3)確定模型的模擬效率評(píng)判準(zhǔn)則，對(duì)所有的參數(shù)組Xn進(jìn)行水文模擬；

2.3 目標(biāo)函數(shù)

2.3.1 單一目標(biāo)優(yōu)化

對(duì)于流域自身的參數(shù)率定，采用單一目標(biāo)優(yōu)化的方法，選取納什效率系數(shù)(NS)為評(píng)價(jià)指標(biāo)，本研究中將編號(hào)為i的流域中參數(shù)組θ的目標(biāo)函數(shù)表示為Oi(θ)，其計(jì)算公式如下：

(2)

NS用于描述洪水過程模擬序列與實(shí)測(cè)序列之間偏差的大小，NS值越接近1，表明模擬序列偏離實(shí)測(cè)序列的程度越小，模擬精度越高。

2.3.2 多目標(biāo)優(yōu)化

對(duì)于多個(gè)流域同步率定，需要綜合考慮所有參與率定流域的模擬效果。本文通過構(gòu)建基于折衷規(guī)劃算法的多目標(biāo)優(yōu)化函數(shù)來實(shí)現(xiàn)多個(gè)流域水文模型的同時(shí)率定，獲取能夠適用于所有參與率定流域的參數(shù)組。由于每個(gè)流域都有自身目標(biāo)函數(shù)的理想解(即單一目標(biāo)優(yōu)化時(shí)的最優(yōu)解)，而目標(biāo)函數(shù)之間存在一定的沖突，所以考慮折衷規(guī)劃的思路，尋求與所有理想解距離最小的折衷解。由于本研究中的目標(biāo)函數(shù)具有相同的量度，于是簡(jiǎn)化為如下目標(biāo)函數(shù)：

(3)

3 結(jié)果與分析

3.1 單個(gè)流域參數(shù)率定與移植

基于上述ROPE參數(shù)率定優(yōu)選算法，選擇單一目標(biāo)函數(shù)，采用HBV模型對(duì)7個(gè)研究流域分別進(jìn)行參數(shù)率定，并將單個(gè)流域的率定參數(shù)移用到其他6個(gè)流域，從而分析單個(gè)流域參數(shù)率定和移植的效果。

通過ROPE方法迭代，每個(gè)流域均得到了10 000組模擬結(jié)果非常接近的最佳參數(shù)組，由于這10 000組參數(shù)的模擬結(jié)果非常接近，為簡(jiǎn)單起見，

采用10 000組參數(shù)的模擬均值作為每個(gè)流域的模擬結(jié)果。圖2展示了HBV模型在不同流域率定期與驗(yàn)證期的10 000組參數(shù)NS效率系數(shù)的平均值。可以看出模擬結(jié)果NS平均值均高于0.8，說明基于ROPE方法進(jìn)行參數(shù)優(yōu)選的HBV模型能很好地用于屯溪嵌套流域的洪水模擬，精度較高。

圖2 單個(gè)流域率定參數(shù)的模擬NS均值Fig.2 Model performance for single catchment calibration and validation

圖3顯示了單個(gè)流域率定的10 000組參數(shù)移植到其他6個(gè)流域的模擬結(jié)果平均值。由圖3可以看出，單個(gè)流域率定參數(shù)結(jié)果均優(yōu)于移用其他流域參數(shù)的模擬結(jié)果。此外，不同流域的參數(shù)移植效果不同，大部分流域通過移植其他流域參數(shù)，能夠獲取理想的模擬結(jié)果(WA，CC，TX，LK，YT)。同時(shí)，部分移植參數(shù)的模擬結(jié)果較差(XT，YC)，表明這些流域率定的參數(shù)不適合移植到其他流域。

圖3 移植參數(shù)的模擬效果NS平均值Fig.3 Model performance by transferring parameters from other catchments

本文通過繪制流域參數(shù)率定和移植的結(jié)果矩陣圖(圖4)，以進(jìn)一步分析流域之間的參數(shù)移植效果。從矩陣圖中可以明顯看出，來自不同流域的移植參數(shù)的模擬效果值NS出現(xiàn)較大的變化。同時(shí)可以看出，該參數(shù)移植效果矩陣是不對(duì)稱的，即存在兩個(gè)流域A和B，當(dāng)A流域率定的參數(shù)用于B流域時(shí)，可以取得理想的模擬結(jié)果，而B流域率定的參數(shù)卻不適用于A流域的模擬。以YT流域和YC流域的NS指標(biāo)為例，使用YC流域率定得到的參數(shù)在YT流域模擬，模擬效果NS平均值比YT流域自身率定NS值低0.2，移植效果差；而使用YT流域率定得到的參數(shù)在YC流域模擬，模擬效果NS平均值比YC流域自身率定NS值低0.07，移植效果合理。這種不對(duì)稱現(xiàn)象，進(jìn)一步印證了“過度擬合”現(xiàn)象的存在。即率定期內(nèi)某些主要的徑流過程可能在一個(gè)流域內(nèi)出現(xiàn)頻繁，參數(shù)在率定過程中體現(xiàn)了該過程的影響，使得參數(shù)的可移植性增強(qiáng)；而在另一個(gè)流域則可能歷史洪水信息不夠完整，參數(shù)對(duì)氣象條件“過度擬合”，從而導(dǎo)致參數(shù)從一個(gè)流域到另一個(gè)流域的移植性較弱。另外考慮到由于屯溪流域及其嵌套子流域間的坡度值存在較大差異，對(duì)于匯流過程影響較大，造成流域的調(diào)蓄作用強(qiáng)弱不一，從而參數(shù)移植效果良莠不一。

圖4 參數(shù)移植的NS平均值矩陣Fig.4 Color code matrix for parameter transfer

3.2 多流域同步率定

獲取多個(gè)流域能夠共同適用的模型參數(shù)，是降低參數(shù)不確定性、提高參數(shù)的可移植性的關(guān)鍵。本文通過折衷規(guī)劃方法構(gòu)建多個(gè)流域同時(shí)率定的多目標(biāo)優(yōu)化函數(shù)，并利用ROPE參數(shù)優(yōu)選方法，對(duì)所有7個(gè)流域的模型參數(shù)進(jìn)行同步率定。通過ROPE方法迭代，優(yōu)選出了在7個(gè)流域均具有較好模擬效果的10 000 組參數(shù)，這10 000 組參數(shù)在7個(gè)流域模擬結(jié)果的NS最大值與最小值的波動(dòng)都很小。7個(gè)流域通過同時(shí)率定獲得的NS效率系數(shù)平均值分別為0.89、0.87、0.87、0.82、0.78、0.87、0.83。圖5顯示了7個(gè)研究流域分別在率定期與驗(yàn)證期進(jìn)行單個(gè)流域率定的NS平均值，和7個(gè)流域的同步率定的參數(shù)模擬NS平均值，以及單個(gè)流域在率定期之間參數(shù)移植的模擬結(jié)果。由圖5可以看出，對(duì)于率定期來說，同步率定的參數(shù)模擬效果略差于流域自身率定的效果，但二者差距并不大，在部分流域(LK)二者的模擬效果相當(dāng)。而將同步率定的參數(shù)模擬效果，與單個(gè)流域率定參數(shù)的移植效果對(duì)比，可以明顯發(fā)現(xiàn)，多個(gè)流域同步率定獲取的參數(shù)組模擬效果良好，優(yōu)于大多數(shù)單個(gè)流域率定參數(shù)的移植效果。驗(yàn)證期的模擬結(jié)果驗(yàn)證了這一結(jié)論，在絕大多數(shù)流域移植參數(shù)模擬結(jié)果確定性系數(shù)與自身接近，只在YC流域下降略顯明顯，分析其原因可能是YC流域面積較其他流域明顯為小，流域的調(diào)蓄作用較弱，造成移植參數(shù)模擬效果不好。綜合以上說明同步率定獲得的參數(shù)組較直接移植更為穩(wěn)定，在一定程度上能夠降低“過度擬合”的影響。

圖5 同步率定的模擬效果NS平均值Fig.5 Model performance for simultaneous calibration and validation

3.3 基于留一法的參數(shù)移植

在上一節(jié)的同步率定中，所有流域都參與了模型率定，在本節(jié)中，采用留一法以進(jìn)一步驗(yàn)證同時(shí)率定情況下參數(shù)移植效果。即對(duì)于所選的7個(gè)流域，每次保留一個(gè)流域假定為無資料流域，利用剩余6個(gè)流域進(jìn)行同時(shí)率定，并將同時(shí)率定獲得的參數(shù)移植到該無資料流域。圖6為其在率定期與驗(yàn)證期的模擬結(jié)果NS平均值。

圖6 留一法同步率定的模擬效果NS平均值Fig.6 Model performance for leave one out calibration and validation

從圖6中可以看出，對(duì)于率定期來說，與上一節(jié)7個(gè)流域的同步率定的模擬結(jié)果NS平均值相比，留一法多流域同步率定的模擬結(jié)果NS平均值因未使用目標(biāo)流域的數(shù)據(jù)，模擬結(jié)果略差于采用所有流域進(jìn)行同步率定的結(jié)果，但與單個(gè)流域之間的參數(shù)移植結(jié)果NS平均值相比，留一法多流域同步率定的參數(shù)移植到TX、WA、CC和LK流域的結(jié)果具有明顯優(yōu)勢(shì)。同時(shí)注意到，留一法多流域同步率定的模擬結(jié)果NS平均值與流域自身率定的NS平均值相差不大，甚至在部分流域(LK)二者模擬效果相當(dāng)。驗(yàn)證期的模擬結(jié)果同樣滿足該規(guī)律，注意到與上節(jié)同步率定相比，留一法移植后的參數(shù)模擬結(jié)果NS值下降非常微小，移植參數(shù)具有較高的確定性。該結(jié)果表明，在無資料地區(qū)洪水模擬中，可用多流域同步率定法進(jìn)行參數(shù)移植，以便獲得較為穩(wěn)定的移植效果，該方法在一定程度上降低了模型參數(shù)的不確定性，這將有助于提高無資料地區(qū)洪水預(yù)報(bào)的精度。

4 結(jié) 論

(1)采用基于深度函數(shù)的ROPE參數(shù)優(yōu)選方法，HBV模型能很好地用于屯溪嵌套流域的洪水過程模擬，可以取得較好的模擬精度。該參數(shù)優(yōu)選方法能夠獲取具有良好穩(wěn)定性的模型參數(shù)，是水文預(yù)報(bào)參數(shù)選擇的一種新途徑。

(2)對(duì)于屯溪及其嵌套的6個(gè)子流域，單個(gè)流域之間的HBV模型參數(shù)移植結(jié)果參差不一，應(yīng)考慮到觀測(cè)資料對(duì)于參數(shù)不確定性的影響，不宜直接進(jìn)行參數(shù)移植。

(3)多流域同步率定能有效降低HBV模型參數(shù)不確定性的影響，并獲得較優(yōu)的移植效果，提高無資料地區(qū)參數(shù)移植穩(wěn)定性和洪水預(yù)報(bào)的參數(shù)確定性。水文模型不確定性來源眾多，而本文僅討論了針對(duì)模型參數(shù)的不確定性，對(duì)于自然不確定性、數(shù)據(jù)不確定性以及模型結(jié)構(gòu)的不確定性還需加以研究，多流域同步率定能否普遍適用于水文模型的參數(shù)移植尚需進(jìn)一步驗(yàn)證。同時(shí)，在多流域同步率定方法中還應(yīng)進(jìn)一步加強(qiáng)對(duì)水文相似性、流域個(gè)數(shù)等因素的考慮。本文所進(jìn)行的是逐日資料的模擬，因此在模擬中主要關(guān)注了長(zhǎng)期的水量平衡，所以未對(duì)峰現(xiàn)時(shí)間和洪峰模擬結(jié)果進(jìn)行比對(duì)。在中小河流洪水預(yù)報(bào)中，小時(shí)尺度洪水的模擬更為重要，同步率定方法在次洪模型中的應(yīng)用及洪水主要特征值的分析是我們下一步研究的重點(diǎn)。

□