PA-DDS算法在HBV模型參數(shù)優(yōu)化中的應(yīng)用

代 旭, 陳元芳

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

0 引 言

水文模型是模擬降雨徑流關(guān)系的工具，這些模型在結(jié)合流域產(chǎn)匯流機(jī)制的基礎(chǔ)上定量分析了出流過程線[1]。模型模擬的精度跟許多因素有關(guān)，比如資料的誤差，模型結(jié)構(gòu)，參數(shù)的相關(guān)性及其準(zhǔn)確性等。模型參數(shù)優(yōu)化在確定了模型的結(jié)構(gòu)以后對(duì)模型的模擬結(jié)果起重大作用。

在水文模型參數(shù)率定的兩種方法中，試錯(cuò)法是人工率定的主要方法，簡(jiǎn)單易行，但需要大量的計(jì)算時(shí)間[2]。自動(dòng)率定方法是目前發(fā)展較快且應(yīng)用廣泛的參數(shù)率定方法。自動(dòng)率定法包括局部尋優(yōu)和全局尋優(yōu)兩種方法。針對(duì)單峰函數(shù)，局部尋優(yōu)法可以有效快速地找到最優(yōu)解，而全局尋優(yōu)法通過全局尋優(yōu)避免最優(yōu)解限制于局部空間內(nèi)。

水文模型參數(shù)率定方法分為單目標(biāo)[3]和多目標(biāo)率定兩種[4]。多目標(biāo)優(yōu)化算法作為全局優(yōu)化的重要組成部分，其目的在于搜尋可行的參數(shù)空間并找出其Pareto最優(yōu)解集[5]。目前多目標(biāo)優(yōu)化算法有NSGA-Ⅱ[6]法、MOGA法[7]、PAES法[8]、AMALGAM法[9]等。其中NSGA-Ⅱ算法作為其他多目標(biāo)尋優(yōu)算法性能檢驗(yàn)的標(biāo)準(zhǔn)，具有尋優(yōu)速度快，解集分布廣泛且均勻的優(yōu)點(diǎn)。AMALGAM算法通過對(duì)四種相關(guān)算法分配權(quán)重從而實(shí)現(xiàn)信息交換同時(shí)尋優(yōu)，解的收斂性能較好。本文將PA-DDS算法與NSGA-Ⅱ算法和AMALGAM算法在收斂性能方面進(jìn)行了對(duì)比，比較了PA-DDS算法與AMALGAM算法的Pareto前沿，并將其用于HBV參數(shù)優(yōu)化中，結(jié)合巴格瑪?shù)俸恿饔蛳嚓P(guān)資料驗(yàn)證了該方法的可行性和適用性。

1 PA-DDS算法

Tolson和Shoemaker于2007年提出了動(dòng)態(tài)維度搜索(Dynamically Dimensioned Search, DDS)算法[10]，并與混合競(jìng)爭(zhēng)優(yōu)化(Shuffled Complex Evolutions, SCE[11])算法進(jìn)行詳細(xì)比較，結(jié)果表明：當(dāng)?shù)螖?shù)有限且維度不小于10時(shí)，DDS算法明顯優(yōu)于SCE算法。Tolson和Asadzadeh在Pareto存檔進(jìn)化策略(Pareto Archived Evolution Strategy, PAES)的啟發(fā)下，使用DDS搜索機(jī)制，并在搜索過程中存檔所有非劣解，提出了處理多目標(biāo)問題的Pareto存檔動(dòng)態(tài)維度搜索(Pareto-Archived Dynamically Dimensioned Search, PA-DDS)算法[12]，應(yīng)用實(shí)例表明PA-DDS算法具有較高的計(jì)算效率。

PA-DDS方法的主要步驟(圖1)為：

步驟一，定義具有n個(gè)目標(biāo)函數(shù)的多目標(biāo)優(yōu)化問題。

步驟二，利用DDS算法實(shí)行n個(gè)優(yōu)化路徑從而最小化每一個(gè)目標(biāo)函數(shù)，同時(shí)存檔外部解集中的所有非劣解。

步驟三，計(jì)算外部解集中解的擁擠距離，并基于此選擇非劣解。

步驟四，通過擾動(dòng)參數(shù)產(chǎn)生新解，并將新解與外部解集中的解進(jìn)行比較，如果未達(dá)到計(jì)算量，且新解是非劣解，那么新解即為當(dāng)前最優(yōu)解，否則重新對(duì)當(dāng)前解擾動(dòng)產(chǎn)生新解；若達(dá)到計(jì)算量，則停止迭代，外部解集中的解即為折中解的估計(jì)值。

2 HBV模型

HBV(Hydrologiska Fyrans Vattenbalans Modell)主要用于河流流量預(yù)測(cè)和河流污染物傳播，由瑞典水利氣象研究中心開發(fā)[13]。張建新[14]等將HBV模型應(yīng)用于中國(guó)東北多冰雪地區(qū)，建立了考慮融雪的洪水預(yù)報(bào)系統(tǒng)；張小琳[15]等對(duì)HBV的實(shí)際蒸散發(fā)模型進(jìn)行了評(píng)估；李海川[16]等將HBV模型應(yīng)用到湯旺河流域進(jìn)行參數(shù)敏感性分析，并在此基礎(chǔ)上率定參數(shù)。

2.1 模型特征

HBV模型是半分布式的概念模型，包括氣象插值、積雪累積融化、蒸散發(fā)估計(jì)子程序、土壤水分計(jì)算程序、徑流產(chǎn)生程序、路徑選擇程序。模型輸入數(shù)據(jù)包括觀測(cè)雨量、氣溫以及可能蒸散發(fā)估計(jì)量。氣溫用來計(jì)算積雪累積和融化，在無積雪地區(qū)溫度項(xiàng)可以忽略。

2.2 模型參數(shù)

模型參數(shù)包括8個(gè)初始給定參數(shù)，分別是：SP為雪蓋厚度，單位為mm，WC為雪的等效水量，單位為mm，SM為土壤水分，單位為mm，UZ為上層土壤需水量，單位為mm，LZ為下層土壤蓄水量，單位為mm，WHC為雪的滯水能力因子，CFR為再冰凍因子，INOUT表示地下水入流量/出流量，Inflow>0，Outflow<0。8個(gè)初始給定參數(shù)的數(shù)值見表1。

表1 HBV模型初始給定參數(shù)及其數(shù)值Tab.1 Initial parameters and their numerical values of HBV model

圖1 PA-DOS流程圖Fig.1 Flow chart of PA-DOS

模型還包括9個(gè)待優(yōu)化可變參數(shù)，分別是：FC為土壤蓄水量閾值，單位為mm，LP為潛在蒸散發(fā)能力閾值，單位為mm，ALFA為描述隨著響應(yīng)水庫(kù)含水量的減少?gòu)搅鳟a(chǎn)生非線性衰減的參數(shù)，BETA為控制響應(yīng)函數(shù)的參數(shù)或每毫米降雨積雪融化所得到的土壤蓄水量的增加值，K為上層響應(yīng)水庫(kù)的消退系數(shù)，K4為下層響應(yīng)水庫(kù)的消退系數(shù)，PERC為從上層相應(yīng)水庫(kù)到下層響應(yīng)水庫(kù)的滲漏，CFLUX為毛管含水量的閾值，MAXBAS為轉(zhuǎn)移函數(shù)參數(shù)。

3 優(yōu)化過程

利用PA-DDS優(yōu)化算法對(duì)HBV模型9個(gè)待優(yōu)化可變參數(shù)進(jìn)行優(yōu)化。

3.1 參數(shù)初始值及尋優(yōu)空間

基于研究流域的特點(diǎn)及參數(shù)的物理意義，給出9個(gè)參數(shù)的初始值并通過人工率定法給出參數(shù)的閾值(表2)。

表2 HBV模型待優(yōu)化參數(shù)的初始值及閾值Tab.2 The initial value and threshold value ofthe optimized parameters

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

Madsen[17]和Mertens[18]提出了對(duì)于優(yōu)化算法以均方根誤差(RMSE)為目標(biāo)函數(shù)，因此本文目標(biāo)函數(shù)選作高低流量均方根誤差，對(duì)應(yīng)的方程分別是：

(2)

式中：Qobs、Pobs分別為高低流量觀測(cè)值；Qcal、Pcal分別為高低流量模擬值；n為實(shí)測(cè)資的長(zhǎng)度。

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

巴格瑪?shù)俸游挥谀岵礌柤拥聺M都谷地，發(fā)源于肖普里山，匯入恒河，降水是主要補(bǔ)給源，流域面積2 900 km2。

本文將PA-DDS算法收斂性能與NSGA-Ⅱ算法和AMALGAM算法進(jìn)行對(duì)比(圖2)，并將非劣解分布的均勻性及解的相似性方面與AMALGAM進(jìn)行比較(圖3)，選取巴格瑪?shù)俸?005-2011年期間洪水日徑流過程數(shù)據(jù)進(jìn)行分析，多目標(biāo)同時(shí)考慮高流量與低流量的均方根誤差，優(yōu)化過程采用PA-DDS算法進(jìn)行計(jì)算，分析生成Pareto最優(yōu)解，從圖中可以看出來Pareto曲線為下凸曲線，有明顯的拐點(diǎn)，從而滿足獲得最優(yōu)解的條件。

4.1 算法收斂性能分析

基于超體積(Hyper-Volume)在評(píng)價(jià)Pareto前端收斂性、寬廣性和均勻性時(shí)表現(xiàn)出來的客觀性[19]，本文選擇超體積作為指標(biāo)來判斷PA-DDS算法與NSGA-Ⅱ算法和AMALGAM算法的求解質(zhì)量。其定義如下：

Hv=volume(∪NPFi=1vi)

(3)

式中：NPF為最終Pareto前沿上解的個(gè)數(shù)；vi為Pareto前沿上第i個(gè)解與參照點(diǎn)圍成的體積。

圖2 不同算法超體積收斂曲線Fig.2 Hyper-Volume convergence curve of different algorithms

從圖2中可以看出，由于AMALGAM算法與NSGA-Ⅱ算法初始種群數(shù)都要大于PA-DDS算法，所以在優(yōu)化過程初期，前者超體積大于后者。對(duì)于目標(biāo)函數(shù)最小化的雙目標(biāo)優(yōu)化問題，超體積值越大，表明Pareto前沿解的分布廣度越大。從圖2中也可以看出，當(dāng)?shù)螖?shù)大于200次時(shí)，PA-DDS非劣解的廣度都要優(yōu)于其他兩種方法，同時(shí)AMALGAM算法、NSGA-Ⅱ算法在迭代次數(shù)分別達(dá)到100次和200次時(shí)，超體積值穩(wěn)定，算法收斂；而PADDS算法當(dāng)?shù)螖?shù)達(dá)到270次時(shí)才趨于穩(wěn)定，表明AMALGAM算法和NSGA-Ⅱ算法在HBV參數(shù)優(yōu)化過程中相比于PA-DDS算法更快趨向于收斂。

4.2 非劣解分布分析

從圖2中可以看出，PA-DDS算法和AMALGAM算法在達(dá)到收斂時(shí)兩者對(duì)應(yīng)的超體積值相對(duì)接近，所以將兩種算法求解的目標(biāo)函數(shù)值進(jìn)行對(duì)比分析。將兩種算法的Pareto前沿繪制在圖3中。

圖3 不同算法得到的Pareto前沿Fig.3 The Pareto front obtained by different algorithms

基于PA-DDS算法的動(dòng)態(tài)存檔機(jī)制，其在求解過程中不會(huì)剔除所有非劣解，所以從圖3中可以看出，PA-DDS算法所得到的Pareto前沿的廣度和數(shù)量要優(yōu)于AMALGAM算法。同時(shí)在目標(biāo)函數(shù)同時(shí)達(dá)到最小值的區(qū)域附近，兩種算法Pareto解集比較接近，說明適用于巴格瑪?shù)俸恿饔虻腍BV模型參數(shù)優(yōu)化問題中兩種算法的尋優(yōu)性能相似，但PA-DDS算法非劣解分布更加密集且均勻，其算法表現(xiàn)相對(duì)更優(yōu)。

4.3 尋優(yōu)速度及穩(wěn)定性分析

為分析PA-DDS算法和AMALGAM算法的尋優(yōu)速度和穩(wěn)定性，記錄PA-DDS算法和AMALGAM算法各自獨(dú)立運(yùn)行15次的時(shí)間，并將兩者得到的Pareto前端繪于圖4。

統(tǒng)計(jì)15次實(shí)驗(yàn)運(yùn)行時(shí)間PA-DDS算法平均耗時(shí)135 s，AMALGAM算法平均耗時(shí)144 s，所以前者尋優(yōu)效率優(yōu)于后者。從圖4中可以看出，兩種算法在目標(biāo)函數(shù)同時(shí)達(dá)到最小值的附近區(qū)域，解集分布都相對(duì)集中；在其中一個(gè)目標(biāo)函數(shù)趨于較大值，另一個(gè)目標(biāo)函數(shù)趨于較小值的區(qū)域附近，解集分布都相對(duì)分散，表明兩種算法能得到質(zhì)量較高的非劣解。同時(shí)可以看出，無論在哪個(gè)區(qū)域，PA-DDS算法所得非劣解集分布更為集中，表明該算法在HBV模型參數(shù)優(yōu)化問題中具有良好的穩(wěn)定性。

圖4 不同算法15次獨(dú)立實(shí)驗(yàn)得到的Pareto前沿Fig.4 The Pareto front obtained by 15 independent experiments with different algorithms

4.4 優(yōu)化結(jié)果

將PA-DDS算法所得Pareto前沿中目標(biāo)函數(shù)同時(shí)達(dá)到最小值時(shí)的參數(shù)組作為最優(yōu)參數(shù)組，結(jié)果見表3。選取2013年作為驗(yàn)證期來檢驗(yàn)PA-DDS算法參數(shù)優(yōu)化結(jié)果(圖5、表4)。

表3 最優(yōu)參數(shù)組數(shù)值Tab.3 The value of optimum parameter array

圖5 2013年流量擬合效果圖Fig.5 Discharge Fitting effect in 2013

表4 PA-DDS算法與HBV結(jié)合的結(jié)果統(tǒng)計(jì)表Tab.4 PA-DDS algorithm combined with HBV results statistics

由表4可以看出，用PA-DDS多目標(biāo)優(yōu)化算法對(duì)HBV模型參數(shù)優(yōu)化結(jié)果，擬合5場(chǎng)洪水確定性系數(shù)平均值為0.86，其中0.85以上的有3場(chǎng)，接近0.9的有2場(chǎng)。以上結(jié)果表明，用PA-DDS多目標(biāo)優(yōu)化算法對(duì)HBV模型參數(shù)優(yōu)化的結(jié)果令人滿意。

5 結(jié) 語

(1)本文將PA-DDS多目標(biāo)優(yōu)化算法在收斂性能方面對(duì)比NSGA-Ⅱ算法和AMALGAM算法，并將非劣解分布的均勻性及解的相似性方面與AMALGAM進(jìn)行比較，利用尼泊爾巴格瑪?shù)俸恿饔蛸Y料，采用PA-DDS算法對(duì)HBV模型參數(shù)進(jìn)行優(yōu)化，并選取五場(chǎng)歷史洪水對(duì)優(yōu)化結(jié)果進(jìn)行了檢驗(yàn)，發(fā)現(xiàn)PA-DDS算法尋優(yōu)速度快，非劣解分布范圍廣并且穩(wěn)定，以確定性系數(shù)為指標(biāo)，模型擬合效果較好。

(2)分析PA-DDS算法對(duì)于其他水文模型的適應(yīng)性是今后研究的重點(diǎn)。

□

[1] 唐運(yùn)憶，欒承梅. SCE-UA算法在新安江模型及TOPMODEL參數(shù)優(yōu)化應(yīng)用中的研究[J]. 水文，2007，(6)：33-35,74.

[2] 欒承梅. 流域水文模型參數(shù)優(yōu)化問題研究[D]. 南京：河海大學(xué)，2005.

[3] 江 燕，劉昌明，武夏寧. 水文模型參數(shù)優(yōu)選的改進(jìn)粒子群算法參數(shù)分析[J]. 水電能源科學(xué),2009,(1):24-27.

[4] 王宇暉，雷曉輝，蔣云鐘，等. Hymod模型參數(shù)敏感性分析和多目標(biāo)優(yōu)化[J]. 水電能源科學(xué)，2010，(11):15-17,122.

[5] 趙香山. 基坑工程變形計(jì)算的土層參數(shù)確定方法分析[D]. 上海：上海交通大學(xué)，2015.

[6] Bekele E G，Ncklow J W. Multi-obijective Automatic Calibration of SWAT Using NSGA-Ⅱ[J]. Journal of Hydrology，2007，341(3-4)：165-176.

[7] 劉 彤. 解多目標(biāo)優(yōu)化問題的進(jìn)化算法[D].西安：西安電子科技大學(xué),2010.

[8] F.Luna，A.J.Nebro，E.Alba. Observations in using Grid-enabled technologies for solving multi-objective optimization problems[J]. Parallel Computing. 2006，32：377-393.

[9] Jasper A Vrugt，Bruce A Robinson. Improved Evolutionary Optimization from Genetically Adaptive Multimethod Search[J]. Proceedings of the National Academy of Sciences of the United States of America.2007,104(3)：708-711.

[10] Tolson B A，Shoemaker C A. Dynamically dimensioned search algorithm for computationally effcient watershed model calibration[J]. Water Resources Research，2007，43(1)：W01413.1-W01413.16.

[11] Duan，Q.，V.K.Gupta，and S.Sorooshian. Shuffled complex evolution approach for effective and effcient global minimization[J]. J.Optim.Theory Appl，1993，76:501-521.

[12] Asadzadeh M，Tolson B A. A new multi-objective algorithm，Pareto archived DDS[C]∥Proceedings of the 11th Annual Conference Companion on Genetic and Evolutionary Computation Conference：Late Breaking Papers. ACM，2013.

[13] Bergstr?m. The HBV model[M]∥ Singh V P. Computer Models of Watershed Hydrology. Little Town，Colorado，USA：Water Ressources Publications，1995：443-476.

[14] 張建新，趙孟芹，章樹安，等. HBV模型在中國(guó)東北多冰雪地區(qū)的應(yīng)用研究[J]. 水文，2007，(4)：31-34.

[15] 張小琳，王衛(wèi)光，陳曙光，等. 基于HBV模型的實(shí)際蒸散發(fā)估算模型評(píng)估[J]. 水電能源科學(xué)，2015，(2)：15-18,41.

[16] 李海川，郝振純，尼瑪?shù)┰觯? HBV模型在湯旺河流域的應(yīng)用研究[J]. 三峽大學(xué)學(xué)報(bào)(自然科學(xué)版)，2016，(5)：24-28,35.

[17] Madsen H. Parameter estimation in distributed hydrological catchment modelling using automatic cali-bration with multiple objectives[J]. Advances in Water Resources，2003，26(2)：205-216.

[18] Mertens J，Madsen H，F(xiàn)eyen L. Including prior information in the estimation of effective soil para-meters in unsaturated zone modelling[J]. Journal of Hydrology，2004，294(4)：251-269.

[19] Konwles J， Corne D. Properties of an adaptive archiving algorithm for storing nondominatedvectors[J]. Evolutionary Computation，IEEE Transactions on，2003，7(2)：100-116.