基于MOSCEM-UA的水文模型多目標參數優化研究

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

(1.山東農業大學 水利土木工程學院，山東 泰安 271018； 2.南水北調東線山東干線有限責任公司，山東 濟南 250000)

0 引 言

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

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

以沁水河流域為例，采用MOSCEM-UA算法基于垂向混合產流模型的基礎上，對流域的降雨、徑流、蒸發、洪水資料進行次洪模型多目標參數優化。

1 研究區垂向混合產流模型

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

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

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

(1)

RS=P-FA

(3)

式中：FM為流域平均下滲能力，mm；fc為穩定下滲率，mm；KF為滲透系數；WM為流域平均蓄水容量；W為流域實際土壤含水量，mm；BF為反映下滲能力空間分布特征的參數；P為扣除雨間蒸發的降雨量，mm。

地下徑流RR采用蓄滿產流結構計算，計算公式為：

(4)

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

因此，產流總量為R=RS+RR。

2 MOSCEM-UA算法

2.1 MOSCEM-UA算法簡介

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

MOSCEM-UA是對混合復雜演化都市(SCEM-UA)全局優化算法的一種改進，它使用Pareto主導的概念，而不是直接進行單目標函數評估，將最初的點群發展為一組解決方案穩定分布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 算法參數設置

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

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

3 模型參數的多目標優化

3.1 Pareto非劣解原理簡介

對于單目標問題，確定目標函數的最優非劣解時常采用權重法對各目標函數賦予權重系數來確定。但在多目標參數優化時，權重法會導致“均化效應”，難以準確評價觀測數據的多樣性特征。因此，基于多目標函數參數優化不再是單一的“最佳”參數集，而是包括在可行參數空間中對多目標之間折中解決方案的Pareto集。

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

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

(5)

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

3.2 優選參數設置

垂向混合產流模型進行次洪模型模擬時以Δt為時段長，模型共有15個參數見表1，其中不敏感參數直接取其中間值，其中模型模擬較為敏感5個參數為BF、KF、fc、CS、L。

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

3.3 目標函數及精度評價指標

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

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

(6)

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

(7)

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

(8)

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

(9)

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

(10)

(11)

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

4 實例應用

4.1 研究區概況

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

4.2 結果分析

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

圖4 Pareto解對應參數歸一化空間Fig.4 Pareto solution corresponding parameter normalized space

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

圖5 三目標參數優化Pareto解散點圖Fig.5 The three-dimensional plots of the Pareto solutions

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

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

表2 垂向混合產流模型預報精度評價結果Tab2 Vertical mixed flow model forecast accuracy evaluation results

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

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

5 結 論

本文應用MOSCEM-UA算法基于改進的垂向混合產流模型上，考慮水資源開發利用及水利工程的影響，在沁水河流域參數優化研究中應用。主要結論如下：

圖8 驗證期19970820場洪水驗證期流量模擬過程Fig.8 In validation period 19970820 flood rate flow simulation process in validation period

表3 垂向混合產流模型參數多目標優化結果Tab 3 Multi-objective optimization results of vertical mixed flow model parameters

(1)在垂向混合產流模型的基礎上應用水文模型多目標優化參數MOSCEM-UA算法，能夠快速高效的生成Pareto非劣解集，提高整體的模擬精度。

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

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

(4)在優化過程中，MOSCEM-UA算法對于模型的適應性在國內應用較少，對pareto非劣解參數優選結果的準確性有待提高，以及進一步改善目標參數優化過程中多個目標間的相互約束及模型校準的異參同效現象。