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

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

0 引 言

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

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

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

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

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

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

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

(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)計算，計算公式為：

(4)

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

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

2 MOSCEM-UA算法

2.1 MOSCEM-UA算法簡介

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

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

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

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

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

SCE-UA算法參數(shù)設(shè)置具體為：m=2n+1，q=n+1，s=pm，α=1，β=2n+1。MOSCEM-UA算法有四個參數(shù)由必須具體問題確定：種群大小s，復(fù)合形個數(shù)p，進(jìn)而確定每個復(fù)合形樣本點數(shù)m=s/p，每個復(fù)合形進(jìn)化步數(shù)L，后代接受概率比例因子β，本算法使用L=m/4和β=1/2。因此，MOSCEM算法需要指定變量是s和p。本文所選的垂向混合產(chǎn)流模型中參數(shù)個數(shù)n=15，復(fù)合形個數(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非劣解原理簡介

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

Pareto非劣解集定義了參數(shù)中的最小不確定性，在不犧牲另一個參數(shù)的情況下主觀相對優(yōu)選使F(θ)的特定分量最小化，其目的是同時最小化兩個參數(shù)的目標(biāo)函數(shù)F(θ1)、F(θ2)。Pareto解決方案集見圖3，點A和點B表示每個單獨的目標(biāo)函數(shù)F(θ1)、F(θ2)最小化的解決方案，連接A和B的線段為Pareto解集， 在多目標(biāo)意義上γ是解集Δδ中一個比任何點都具優(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)行次洪模型模擬時以Δt為時段長，模型共有15個參數(shù)見表1，其中不敏感參數(shù)直接取其中間值，其中模型模擬較為敏感5個參數(shù)為BF、KF、fc、CS、L。

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

3.3 目標(biāo)函數(shù)及精度評價指標(biāo)

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

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

(6)

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

(7)

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

(8)

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

(9)

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

(10)

(11)

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

4 實例應(yīng)用

4.1 研究區(qū)概況

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

4.2 結(jié)果分析

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

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

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

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

圖6 Pareto解二維散點圖Fig.6 Pareto solution to two-dimensional scatter plot

預(yù)報方案的效果評定和有效性檢驗，一般將良好代表性的系列資料分為率定期和驗證期，通過模擬與實測水文要素間的比較來檢驗預(yù)報方案的效果與有效性?！端那閳箢A(yù)報規(guī)范》(SL250-2000)規(guī)定：評價和檢驗方法采用統(tǒng)一的許可誤差和有效性標(biāo)準(zhǔn)對預(yù)報方案進(jìn)行評定和校驗。依據(jù)規(guī)范在預(yù)報過程當(dāng)中，洪量、洪峰相對誤差取預(yù)測期內(nèi)實測值的20%作為允許誤差。由表2看出，洪量指標(biāo)RE在四個目標(biāo)函數(shù)方案中的率定期合格率達(dá)到75%，驗證期達(dá)到100%，其中三目標(biāo)同時考慮時誤差最小，其次是洪量目標(biāo)F2；洪峰指標(biāo)REPF在F3中最??；NSE值的精度等級大多數(shù)在乙級及以上，且三目標(biāo)優(yōu)于單目標(biāo)F1、F2及F3，預(yù)報結(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ù)報精度評價結(jié)果Tab2 Vertical mixed flow model forecast accuracy evaluation results

圖7，8 分別是選取率定期和驗證期以協(xié)調(diào)解參數(shù)模擬過程的效果圖，基于洪量目標(biāo)、均方根目標(biāo)、洪峰目標(biāo)及三目標(biāo)的四種組合下進(jìn)行洪水模擬， 19810701場次洪水可以看出，三目標(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)行檢驗期19970820場洪水模擬，且效果較好。結(jié)果表明，多目標(biāo)優(yōu)化在整體流量過程線上吻合趨勢較好，且率定期的效果優(yōu)于驗證期的效果，進(jìn)一步驗證了多目標(biāo)參數(shù)率定的有效性和優(yōu)越性。

圖7 率定期19810701場洪水流量模擬過程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 驗證期19970820場洪水驗證期流量模擬過程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)方法的不同之處在于要最小化包含多個目標(biāo)函數(shù)的向量。因此，優(yōu)化方法的任務(wù)是映射出一個折中表面也稱為pareto前沿，而不是單目標(biāo)情況下尋找單個標(biāo)量值最優(yōu)，本文采用MOSCEM-UA算法對多個目標(biāo)函數(shù)產(chǎn)生Pareto非劣解，并在目標(biāo)函數(shù)空間中定義Pareto前沿來實現(xiàn)。

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

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