基于MOSCEM-UA的水文模型多目標(biāo)參數(shù)優(yōu)化研究

杜彥臻，劉紅利，趙天宇，林洪孝，王 剛

(1.山東農(nóng)業(yè)大學(xué) 水利土木工程學(xué)院，山東 泰安 271018； 2.南水北調(diào)東線山東干線有限責(zé)任公司，山東 濟(jì)南 250000)

0 引 言

國(guó)內(nèi)外對(duì)流域水文模型的研究不斷深入，通過對(duì)水文現(xiàn)象的合理概化，流域水文過程及各水文要素計(jì)算構(gòu)建流域水文模型。隨著計(jì)算機(jī)技術(shù)的發(fā)展，水文模型在水文預(yù)報(bào)、水庫(kù)防洪以及水資源管理等方面應(yīng)用廣泛。世界氣象組織以模型的精度為標(biāo)準(zhǔn)進(jìn)行模型對(duì)比，發(fā)現(xiàn)水文模型在濕潤(rùn)地區(qū)模擬精度較好。其中國(guó)內(nèi)應(yīng)用較多的新安江模型，也主要適應(yīng)于以蓄滿產(chǎn)流為主的濕潤(rùn)地區(qū)，在典型的半濕潤(rùn)、半干旱地區(qū)常會(huì)有蓄滿和超滲產(chǎn)流兩種產(chǎn)流機(jī)制并存的流域。因此，研制適用于半濕潤(rùn)、半干旱地區(qū)的混合產(chǎn)流模型十分必要。水文模型預(yù)報(bào)精度與模型的參數(shù)優(yōu)化率定選擇密切相關(guān)，模型參數(shù)是用來評(píng)估實(shí)測(cè)與模擬數(shù)據(jù)的接近程度。國(guó)內(nèi)外水文學(xué)者對(duì)自動(dòng)校準(zhǔn)方法進(jìn)行大量研究，自動(dòng)優(yōu)化法利用計(jì)算機(jī)的速度和能力，客觀且相對(duì)容易操作，從而優(yōu)選出理想?yún)?shù)。

模型參數(shù)的優(yōu)化方法主要有單目標(biāo)和多目標(biāo)優(yōu)化法，目前國(guó)內(nèi)應(yīng)用較多的是遺傳算法、羅森布洛克法及SCE-UA等單目標(biāo)算法[1]，對(duì)多目標(biāo)優(yōu)化算法應(yīng)用相對(duì)較少。單目標(biāo)算法雖簡(jiǎn)單快捷，但由于水文模型結(jié)構(gòu)不足，單目標(biāo)只能捕捉單一方面的水文特征，不能較全面的對(duì)水文特征進(jìn)行綜合考慮，從而沒有得到水文學(xué)家的廣泛認(rèn)可，要使模擬結(jié)果的準(zhǔn)確性進(jìn)一步提高，需進(jìn)行多目標(biāo)優(yōu)化。水文模型參數(shù)校準(zhǔn)實(shí)踐經(jīng)驗(yàn)表明，預(yù)報(bào)過程中多個(gè)目標(biāo)之間相互制約，任何單目標(biāo)函數(shù)都不足以準(zhǔn)確的評(píng)價(jià)觀測(cè)數(shù)據(jù)的所有特征。因此為提高模型預(yù)報(bào)精度采用多目標(biāo)優(yōu)化方法勢(shì)在必行，國(guó)外應(yīng)用較多的多目標(biāo)方法如：多目標(biāo)粒子群算法[2]、MOSCDE算法[3]、NSGA-Ⅱ[4]、MOCOM[5]以及MOSCEM-UA算法等。

以沁水河流域?yàn)槔?，采用MOSCEM-UA算法基于垂向混合產(chǎn)流模型的基礎(chǔ)上，對(duì)流域的降雨、徑流、蒸發(fā)、洪水資料進(jìn)行次洪模型多目標(biāo)參數(shù)優(yōu)化。

1 研究區(qū)垂向混合產(chǎn)流模型

蓄滿和超滲是兩種典型的產(chǎn)流模型。濕潤(rùn)地區(qū)降雨產(chǎn)流以蓄滿產(chǎn)流為主，干旱地區(qū)以超滲產(chǎn)流為主，即使在典型的濕潤(rùn)或干旱地區(qū)也常會(huì)出現(xiàn)兩種產(chǎn)流機(jī)制并存的情況。包為民[6]教授首次提出垂向混合產(chǎn)流模型[7]，該模型在縱向上將蓄滿—超滲產(chǎn)流進(jìn)行混合。沁水河流域位于北方的半濕潤(rùn)地區(qū)，單一的產(chǎn)流模型不能準(zhǔn)確的模擬研究區(qū)產(chǎn)流量，因此提出適合該流域的混合產(chǎn)流計(jì)算方法十分重要。

模型結(jié)構(gòu)主要分為4部分[8]：①產(chǎn)流部分用蓄滿—超滲垂向混合產(chǎn)流概念；②蒸散發(fā)部分按上、下、深三層蒸發(fā)模式；③水源部分用EX次拋物線形自由水蓄水曲線劃分地表、壤中和地下徑流3種水源；④匯流部分地面匯流用單位線法、地下匯流用線性水庫(kù)法。此外，模型依據(jù)沁水河流域特殊的水文特征，考慮人類活動(dòng)對(duì)水文循環(huán)要素的影響，即加入河道內(nèi)農(nóng)業(yè)灌溉、生活取水和水利工程等影響因素，構(gòu)成了完整的適應(yīng)該流域的垂向混合產(chǎn)流模型。

混合產(chǎn)流下滲計(jì)算，雨量P到達(dá)地面，通過空間分布下滲曲線，劃分地面徑流RS和下滲水量FA，F(xiàn)A向下運(yùn)動(dòng)，土壤缺水量大時(shí)補(bǔ)充含水量，不產(chǎn)流；土壤缺水量小時(shí)補(bǔ)足土壤缺水后，產(chǎn)生地面以下徑流RR。依此劃分為地面徑流和地面以下徑流(壤中流和地下徑流)。地面徑流RS根據(jù)改進(jìn)后的格林—安普特下滲公式[9]計(jì)算，計(jì)算公式為：

(1)

RS=P-FA

(3)

式中：FM為流域平均下滲能力，mm；fc為穩(wěn)定下滲率，mm；KF為滲透系數(shù)；WM為流域平均蓄水容量；W為流域?qū)嶋H土壤含水量，mm；BF為反映下滲能力空間分布特征的參數(shù)；P為扣除雨間蒸發(fā)的降雨量，mm。

地下徑流RR采用蓄滿產(chǎn)流結(jié)構(gòu)計(jì)算，計(jì)算公式為：

(4)

式中：Wmm為蓄水量最大值。

因此，產(chǎn)流總量為R=RS+RR。

2 MOSCEM-UA算法

2.1 MOSCEM-UA算法簡(jiǎn)介

多目標(biāo)混合復(fù)雜演化MOSCEM-UA(Multi-objective Shuffled Complex Evolution Metropolis)算法[10]是Vrugt等[11]在2003年提出了一種新的多目標(biāo)算法，由阿姆斯特丹大學(xué)和亞利桑那大學(xué)合作開發(fā)，是在Duan等提出的SCE-UA算法[12]的擴(kuò)展，它結(jié)合了SCE-UA與馬爾可夫鏈Metropolis強(qiáng)度的一種新的適應(yīng)性分配方法。MOSCEM-UA可被視為建立的水文模型校準(zhǔn)基準(zhǔn)，并成功地進(jìn)行模型參數(shù)校準(zhǔn)，被證明為一種穩(wěn)健和有效的全局優(yōu)化方法。

MOSCEM-UA是對(duì)混合復(fù)雜演化都市(SCEM-UA)全局優(yōu)化算法的一種改進(jìn)，它使用Pareto主導(dǎo)的概念，而不是直接進(jìn)行單目標(biāo)函數(shù)評(píng)估，將最初的點(diǎn)群發(fā)展為一組解決方案穩(wěn)定分布Pareto解集，該算法迭代效率高，逼近速度快。其算法的流程見圖1、圖2。

圖1 SCE-UA算法流程圖Fig.1 Flowchart of the SCE-UA algorithm

圖2 MOSCEM-UA大都會(huì)演化算法流程圖Fig.2 Flowchart of the Sequence Evolution employed in the MOSCEM‐UA algorithm.

2.2 算法參數(shù)設(shè)置

SCE-UA算法是在單純形法基礎(chǔ)上根據(jù)生物競(jìng)爭(zhēng)進(jìn)化和基因算法的原理綜合而成，可一致、快速、合理地搜索到水文模型參數(shù)的全局最優(yōu)。雖然算法參數(shù)較多，但絕大部分都可取經(jīng)驗(yàn)值，只有復(fù)合形個(gè)數(shù)p根據(jù)具體問題確定。p主要決定收斂精度及快慢，在取值范圍內(nèi)p值越大越適應(yīng)于高階非線性問題，根據(jù)沁水河流域水文要素特征本文選取p=8時(shí)對(duì)參數(shù)全局最優(yōu)及參數(shù)間的關(guān)聯(lián)性都相對(duì)較好，且能夠?qū)?shù)結(jié)果進(jìn)行較合理的分析優(yōu)選。

SCE-UA算法參數(shù)設(shè)置具體為：m=2n+1，q=n+1，s=pm，α=1，β=2n+1。MOSCEM-UA算法有四個(gè)參數(shù)由必須具體問題確定：種群大小s，復(fù)合形個(gè)數(shù)p，進(jìn)而確定每個(gè)復(fù)合形樣本點(diǎn)數(shù)m=s/p，每個(gè)復(fù)合形進(jìn)化步數(shù)L，后代接受概率比例因子β，本算法使用L=m/4和β=1/2。因此，MOSCEM算法需要指定變量是s和p。本文所選的垂向混合產(chǎn)流模型中參數(shù)個(gè)數(shù)n=15，復(fù)合形個(gè)數(shù)p=8，種群大小s=8×31=248，最大迭代次數(shù)I=1 000，算法終止條件為達(dá)到設(shè)定最大迭代次數(shù)或經(jīng)過多次循環(huán)后參數(shù)值和目標(biāo)函數(shù)均無明顯提高。

3 模型參數(shù)的多目標(biāo)優(yōu)化

3.1 Pareto非劣解原理簡(jiǎn)介

對(duì)于單目標(biāo)問題，確定目標(biāo)函數(shù)的最優(yōu)非劣解時(shí)常采用權(quán)重法對(duì)各目標(biāo)函數(shù)賦予權(quán)重系數(shù)來確定。但在多目標(biāo)參數(shù)優(yōu)化時(shí)，權(quán)重法會(huì)導(dǎo)致“均化效應(yīng)”，難以準(zhǔn)確評(píng)價(jià)觀測(cè)數(shù)據(jù)的多樣性特征。因此，基于多目標(biāo)函數(shù)參數(shù)優(yōu)化不再是單一的“最佳”參數(shù)集，而是包括在可行參數(shù)空間中對(duì)多目標(biāo)之間折中解決方案的Pareto集。

Pareto非劣解集定義了參數(shù)中的最小不確定性，在不犧牲另一個(gè)參數(shù)的情況下主觀相對(duì)優(yōu)選使F(θ)的特定分量最小化，其目的是同時(shí)最小化兩個(gè)參數(shù)的目標(biāo)函數(shù)F(θ1)、F(θ2)。Pareto解決方案集見圖3，點(diǎn)A和點(diǎn)B表示每個(gè)單獨(dú)的目標(biāo)函數(shù)F(θ1)、F(θ2)最小化的解決方案，連接A和B的線段為Pareto解集， 在多目標(biāo)意義上γ是解集Δδ中一個(gè)比任何點(diǎn)都具優(yōu)越性的元素。數(shù)學(xué)上定義Pareto最優(yōu)集為不劣于參數(shù)空間中任一參數(shù)的解的集合，即參數(shù)θ1劣于參數(shù)θ2，當(dāng)且僅當(dāng)滿足不存在i使

Fi(θ1)Fi(θ2)i=1,…,m

(5)

圖3 Pareto非劣解集圖Fig.3 Pareto noninferior solutions set

3.2 優(yōu)選參數(shù)設(shè)置

垂向混合產(chǎn)流模型進(jìn)行次洪模型模擬時(shí)以Δt為時(shí)段長(zhǎng)，模型共有15個(gè)參數(shù)見表1，其中不敏感參數(shù)直接取其中間值，其中模型模擬較為敏感5個(gè)參數(shù)為BF、KF、fc、CS、L。

表1 模型參數(shù)上下限Tab.1 Model parameters upper and lower limits

3.3 目標(biāo)函數(shù)及精度評(píng)價(jià)指標(biāo)

不同的目標(biāo)函數(shù)用來評(píng)價(jià)水文過程的不同特征，多目標(biāo)優(yōu)選時(shí)，多個(gè)目標(biāo)之間往往存在沖突，當(dāng)某一目標(biāo)函數(shù)較優(yōu)時(shí)必定會(huì)犧牲其他目標(biāo)函數(shù)，故目標(biāo)函數(shù)的選擇對(duì)水文模型的預(yù)報(bào)精度有較大的影響，本文采用3個(gè)目標(biāo)函數(shù)作為垂向混合產(chǎn)流模型次洪模擬標(biāo)準(zhǔn)，公式如下。

(1)均方誤差RMSE：表示實(shí)測(cè)與模擬值的吻合程度。

(6)

(2)洪量目標(biāo)控制：表示整體洪量擬合以保證水量平衡。

(7)

(3)洪峰目標(biāo)控制：表示洪峰相對(duì)誤差側(cè)重于洪峰流量擬合。

(8)

水文模型精度誤差評(píng)定指標(biāo)主要以場(chǎng)次洪量相對(duì)誤差(RE)、洪峰相對(duì)誤差(REPF)以及Nash-Sutcliffe效率系數(shù)(NSE)為標(biāo)準(zhǔn)，具體公式為：

(9)

REPF=|Qo,i-Qs,i|/Qo,i×100%

(10)

(11)

式中：Qo,i為實(shí)測(cè)流量，m3/s；Qs,i為模擬流量，m3/s；m為資料序列長(zhǎng)度；N為場(chǎng)次洪水?dāng)?shù)。

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

4.1 研究區(qū)概況

沁水河流域?qū)倌z東半島低山丘陵區(qū)，地勢(shì)為南高北低，流域?qū)倥瘻貛О霛駶?rùn)大陸性季風(fēng)氣候，多年平均降水量為651.9 mm，年平均氣溫12.7 ℃，四季變化分明。沁水河為典型的山溪性雨源型河道，河道流量與降水量變化規(guī)律相一致，且年內(nèi)、年際變化大，枯季流量較小，洪水主要集中在汛期。沁水河上游河床受水流沖刷較大，中游側(cè)蝕嚴(yán)重，下游沉積，形成流水不暢現(xiàn)象。中下游建立攔河閘壩后，中低水時(shí)，攔蓄上游來水，汛期，受上游閘壩攔蓄影響，洪水過程發(fā)生形變，洪峰比天然河道滯后，洪峰洪量也往往比天然河道偏大。該流域測(cè)驗(yàn)斷面上游共計(jì)有9座小型水庫(kù)，測(cè)驗(yàn)斷面及以上有4座攔河閘壩，測(cè)驗(yàn)斷面以下有2座攔河閘壩，7座塘壩。流域水文站監(jiān)測(cè)斷面以上控制流域面積為168 km2，共有5個(gè)雨量(分別是牟平、十六里頭、玉林店、徐家疃、西柳莊)、1個(gè)流量站(牟平水文站)和1個(gè)蒸發(fā)站(門口水庫(kù))。本研究應(yīng)用該流域1981-2000年20 a的實(shí)測(cè)日降雨、蒸發(fā)、流量資料用于垂向混合產(chǎn)流模型參數(shù)優(yōu)化率定，其中1981-1992年為率定期，1993-2000年為驗(yàn)證期。

4.2 結(jié)果分析

MOSCEM-UA算法優(yōu)化得到Pareto集合的參數(shù)分布見圖4所示，對(duì)15個(gè)參數(shù)值均基于上下限(0， 1)之間進(jìn)行歸一化處理，連接參數(shù)的每條線表示一個(gè)Pareto集。本例選取5個(gè)敏感參數(shù)BF、KF、fc、CS、L作為優(yōu)選對(duì)象對(duì)其取值范圍的邊界進(jìn)行歸一化處理。由圖看出，Pareto集中最優(yōu)參數(shù)的變化幅度較大，取值范圍較小的參數(shù)容易確定，取值范圍較大的參數(shù)不易確定。因此，算法得到Pareto解能夠捕獲大多數(shù)參數(shù)的不確定性信息，且參數(shù)最優(yōu)分布基本在取值范圍之內(nèi)沒有越界趨勢(shì)。

圖4 Pareto解對(duì)應(yīng)參數(shù)歸一化空間Fig.4 Pareto solution corresponding parameter normalized space

采用MOSCEM-UA算法產(chǎn)生Pareto集合對(duì)多目標(biāo)(RMSE目標(biāo)F1、洪量目標(biāo)F2、洪峰目標(biāo)F3)進(jìn)行計(jì)算得到50個(gè)Pareto散點(diǎn)解，三維目標(biāo)參數(shù)優(yōu)化Pareto散點(diǎn)解見圖5所示。二維Pareto解見圖6所示，左邊的黑點(diǎn)表示目標(biāo)函數(shù)空間上的Pareto前沿，Pareto前沿點(diǎn)是最好的單一目標(biāo)函數(shù)解決方案?？梢钥闯?，Pareto前沿解目標(biāo)函數(shù)間存在明顯的制約關(guān)系，兩者不能同時(shí)最優(yōu)，必須犧牲一方來取得彼此的最優(yōu)值，因此MOSCEM-UA算法對(duì)多目標(biāo)優(yōu)化計(jì)算可以得到一組多個(gè)目標(biāo)都較優(yōu)的Pareto非劣解集。

圖5 三目標(biāo)參數(shù)優(yōu)化Pareto解散點(diǎn)圖Fig.5 The three-dimensional plots of the Pareto solutions

圖6 Pareto解二維散點(diǎn)圖Fig.6 Pareto solution to two-dimensional scatter plot

預(yù)報(bào)方案的效果評(píng)定和有效性檢驗(yàn)，一般將良好代表性的系列資料分為率定期和驗(yàn)證期，通過模擬與實(shí)測(cè)水文要素間的比較來檢驗(yàn)預(yù)報(bào)方案的效果與有效性。《水文情報(bào)預(yù)報(bào)規(guī)范》(SL250-2000)規(guī)定：評(píng)價(jià)和檢驗(yàn)方法采用統(tǒng)一的許可誤差和有效性標(biāo)準(zhǔn)對(duì)預(yù)報(bào)方案進(jìn)行評(píng)定和校驗(yàn)。依據(jù)規(guī)范在預(yù)報(bào)過程當(dāng)中，洪量、洪峰相對(duì)誤差取預(yù)測(cè)期內(nèi)實(shí)測(cè)值的20%作為允許誤差。由表2看出，洪量指標(biāo)RE在四個(gè)目標(biāo)函數(shù)方案中的率定期合格率達(dá)到75%，驗(yàn)證期達(dá)到100%，其中三目標(biāo)同時(shí)考慮時(shí)誤差最小，其次是洪量目標(biāo)F2；洪峰指標(biāo)REPF在F3中最?。籒SE值的精度等級(jí)大多數(shù)在乙級(jí)及以上，且三目標(biāo)優(yōu)于單目標(biāo)F1、F2及F3，預(yù)報(bào)結(jié)果較好。由表可以看出，在考慮單目標(biāo)最優(yōu)情況下，最終的優(yōu)選結(jié)果偏向于該目標(biāo)函數(shù)所體現(xiàn)的水文特征，其他目標(biāo)函數(shù)未必能達(dá)到令人滿意的結(jié)果。因此單一目標(biāo)優(yōu)化不能準(zhǔn)確全面的校準(zhǔn)水文模型及水文過程的性能特征，需要通過多目標(biāo)優(yōu)化來完成。

表2 垂向混合產(chǎn)流模型預(yù)報(bào)精度評(píng)價(jià)結(jié)果Tab2 Vertical mixed flow model forecast accuracy evaluation results

圖7，8 分別是選取率定期和驗(yàn)證期以協(xié)調(diào)解參數(shù)模擬過程的效果圖，基于洪量目標(biāo)、均方根目標(biāo)、洪峰目標(biāo)及三目標(biāo)的四種組合下進(jìn)行洪水模擬， 19810701場(chǎng)次洪水可以看出，三目標(biāo)優(yōu)化率定統(tǒng)籌了洪量過程及洪峰峰值等特性，整體結(jié)果較優(yōu)。模型參數(shù)多目標(biāo)優(yōu)化結(jié)果通過最終優(yōu)化得到的pareto參數(shù)協(xié)調(diào)解見表3，由此參數(shù)進(jìn)行檢驗(yàn)期19970820場(chǎng)洪水模擬，且效果較好。結(jié)果表明，多目標(biāo)優(yōu)化在整體流量過程線上吻合趨勢(shì)較好，且率定期的效果優(yōu)于驗(yàn)證期的效果，進(jìn)一步驗(yàn)證了多目標(biāo)參數(shù)率定的有效性和優(yōu)越性。

圖7 率定期19810701場(chǎng)洪水流量模擬過程Fig.7 19810701 flood rate flow simulation process in calibration period

5 結(jié) 論

本文應(yīng)用MOSCEM-UA算法基于改進(jìn)的垂向混合產(chǎn)流模型上，考慮水資源開發(fā)利用及水利工程的影響，在沁水河流域參數(shù)優(yōu)化研究中應(yīng)用。主要結(jié)論如下：

圖8 驗(yàn)證期19970820場(chǎng)洪水驗(yàn)證期流量模擬過程Fig.8 In validation period 19970820 flood rate flow simulation process in validation period

表3 垂向混合產(chǎn)流模型參數(shù)多目標(biāo)優(yōu)化結(jié)果Tab 3 Multi-objective optimization results of vertical mixed flow model parameters

(1)在垂向混合產(chǎn)流模型的基礎(chǔ)上應(yīng)用水文模型多目標(biāo)優(yōu)化參數(shù)MOSCEM-UA算法，能夠快速高效的生成Pareto非劣解集，提高整體的模擬精度。

(2)多目標(biāo)與單目標(biāo)方法的不同之處在于要最小化包含多個(gè)目標(biāo)函數(shù)的向量。因此，優(yōu)化方法的任務(wù)是映射出一個(gè)折中表面也稱為pareto前沿，而不是單目標(biāo)情況下尋找單個(gè)標(biāo)量值最優(yōu)，本文采用MOSCEM-UA算法對(duì)多個(gè)目標(biāo)函數(shù)產(chǎn)生Pareto非劣解，并在目標(biāo)函數(shù)空間中定義Pareto前沿來實(shí)現(xiàn)。

(3)MOSCEM-UA算法結(jié)合了SCE-UA算法復(fù)雜混洗的優(yōu)勢(shì)與Metropolis算法基于概率協(xié)方差的搜索方法，Zitzler和Thiele提出改進(jìn)的適應(yīng)性分配概念，以提高Pareto解集的高效性和一致性。

(4)在優(yōu)化過程中，MOSCEM-UA算法對(duì)于模型的適應(yīng)性在國(guó)內(nèi)應(yīng)用較少，對(duì)pareto非劣解參數(shù)優(yōu)選結(jié)果的準(zhǔn)確性有待提高，以及進(jìn)一步改善目標(biāo)參數(shù)優(yōu)化過程中多個(gè)目標(biāo)間的相互約束及模型校準(zhǔn)的異參同效現(xiàn)象。