余勝男，張雪，邵年華，呂海強

(1.深圳市廣匯源環(huán)境水務(wù)有限公司，廣東深圳518020；2.深圳市西麗水庫管理處，廣東深圳518050)

隨著計算機技術(shù)的發(fā)展，洪水預(yù)報和水資源管理等眾多領(lǐng)域開始逐漸借助流域水文模型進行預(yù)報和管理[1]。流域水文模擬旨在應(yīng)用物理數(shù)學和水文學知識，在流域尺度范圍內(nèi)，對降雨徑流過程進行局部或綜合模擬，從而達到確定流域水文響應(yīng)的目的[2]。目前國內(nèi)外運用比較廣泛的新安江模型，是20世紀60年代中國著名水文學家趙人俊教授提出的流域水文模型，該模型以蓄滿產(chǎn)流為機制，適用于濕潤地區(qū)或半濕潤地區(qū)[3]。水文模型結(jié)構(gòu)比較復(fù)雜，在參數(shù)優(yōu)化過程中解集的選取存在計算效率和精度的矛盾[4]，參數(shù)優(yōu)化方法主要包括試錯法和自動率定法[5]，后者又包括局部尋優(yōu)和全局尋優(yōu)兩種方法。局部尋優(yōu)法可滿足單峰函數(shù)解的要求，而對于復(fù)峰函數(shù)，則需通過全局尋優(yōu)避免最優(yōu)解限制于局部空間內(nèi)[6]，所以全局優(yōu)化算法更多的應(yīng)用于多參數(shù)的水文模型參數(shù)率定中。

多目標優(yōu)化算法[7]作為全局優(yōu)化的重要組成部分，一般來說，多目標優(yōu)化問題并不存在一個最優(yōu)解，所有可能的解都成為非劣解，也稱Pareto解，多目標優(yōu)化目的在于搜尋可行的參數(shù)空間并找出其Pareto最優(yōu)解集。目前廣泛應(yīng)用的多目標優(yōu)化算法有NSGA-Ⅱ[8]法、AMALGAM法等。其中NSGA-Ⅱ算法因其解集分布廣泛且均勻，尋優(yōu)速度快的優(yōu)點已經(jīng)成為檢驗其他多目標尋優(yōu)算法性能的參照標準。AMALGAM算法解的收斂性能較好，近年國內(nèi)外學者也將其運用于數(shù)據(jù)大，模型結(jié)構(gòu)復(fù)雜的參數(shù)優(yōu)化中。PA-DDS算法作為新近提出的多目標優(yōu)化算法，在求解過程中動態(tài)存檔非劣解防止解的丟失，因此解集的廣度及精度較高。國內(nèi)對PA-DDS算法的應(yīng)用研究[9-10]尚處于起步階段，本文首次將其應(yīng)用于多參數(shù)新安江模型參數(shù)優(yōu)化中，且與NSGA-Ⅱ算法和AMALGAM算法在收斂性能和非劣解方面進行了對比，結(jié)合深圳西麗水庫入庫流量過程資料驗證了該方法的優(yōu)越性。

1 西麗水庫概況

西麗水庫(圖1)位于深圳市南山區(qū)西麗街道境內(nèi)，壩址位于珠江三角洲水系的大沙河上游河段右岸支流西麗水上。西麗水庫控制流域面積29.0 km2，總庫容4 000萬m3，多年平均降雨量1 565 mm。水庫于1959年動工興建，原設(shè)計以防洪、灌溉為主，結(jié)合發(fā)電。隨著經(jīng)濟的發(fā)展，西麗水庫成為以防洪、城市供水為主[12]，兼具原水調(diào)蓄功能的中型水庫,是供水網(wǎng)絡(luò)干線(東部供水水源工程)末端的調(diào)節(jié)水庫之一。水庫擔任著南山、寶安500萬人口的供水任務(wù)，水庫一旦失事，將影響居民的飲用水安全[11]。

西麗水庫設(shè)有降雨、水位、蒸發(fā)自動觀測站點。水位計于1997年建設(shè)，位于大壩右壩肩DN2200閘室旁，每隔15 min自動采集一個水位數(shù)據(jù)；雨量筒于2012年建設(shè)，位于西麗水庫管理處辦公樓樓頂，每隔15 min自動采集一個雨量數(shù)據(jù)；蒸發(fā)皿于2014年建設(shè)，位于大壩左壩肩背水坡，每日采集一個蒸發(fā)數(shù)據(jù)。

2 PA-DDS算法簡介

存檔動態(tài)維度搜索算法PA-DDS算法由動態(tài)維度搜索算法DDS發(fā)展而來，動態(tài)維度搜索(Dynamically Dimensioned Search, DDS[12])算法起始于全局搜索，即在全搜索域上產(chǎn)生初始解，隨著迭代的進行，算法通過以服從正態(tài)分布的概率動態(tài)地減少解的變化維度實現(xiàn)搜索空間由全局向局部的轉(zhuǎn)化過程，逐漸局限在一個局部空間內(nèi)，最終求得目標函數(shù)最優(yōu)解，但這種尋優(yōu)過程容易出現(xiàn)將解丟失的情況。在動態(tài)維度搜索DDS算法的啟發(fā)下且為避免該算法的弱點，Tolson和Asadzadeh結(jié)合Pareto存檔進化策略(Pareto Archived Evolution Strategy, PAES)，提出了處理多目標問題的Pareto存檔動態(tài)維度搜索(Pareto-Archived Dynamically Dimensioned Search, PA-DDS[13])算法，該算法能在搜索過程中動態(tài)存檔所有非劣解，防止解的丟失，最終實現(xiàn)模型參數(shù)全局尋優(yōu)。

PA-DDS方法的主要步驟(圖2)為：①確定待解決的優(yōu)化問題及其包含的目標函數(shù)，本研究選擇的目標函數(shù)是相對誤差及確定性系數(shù)；②采用DDS算法初始化種群，通過迭代尋優(yōu)求解最優(yōu)參數(shù)及最優(yōu)函數(shù)值，并存檔非劣解，形成Pareto前沿；③計算擁擠半徑并據(jù)此找出步驟二中的Pareto前沿的最優(yōu)解集；④設(shè)置擾動參數(shù)并對當前最優(yōu)解集產(chǎn)生一定范圍的擾動，重復(fù)步驟二，產(chǎn)生新的解集，若新解集為最優(yōu)非劣解，則替代原解集，否則繼續(xù)重復(fù)步驟二，直到找到最優(yōu)非劣解且計算量達到設(shè)置的閾值為止。

3 優(yōu)化過程

3.1 確定參數(shù)尋優(yōu)空間

參數(shù)尋優(yōu)空間由參數(shù)的上邊界和下邊界組成，可以根據(jù)人工經(jīng)驗法確定，也可以采用推薦的參數(shù)范圍[14]。表1、2分別列出了日模型和次洪模型優(yōu)化參數(shù)的上下邊界。各參數(shù)意義如下：K為流域蒸散發(fā)折算系數(shù)；B為張力水蓄水容量曲線方次；C為深層蒸散發(fā)折算系數(shù)；Wm為流域平均張力水容量；Wum為上層張力水容量；Wlm為下層張力水容量；Im為不透水面積占全流域面積的比例；Ex為表層自由水蓄水容量曲線方次；Sm為表層自由水蓄水容量；Kg為表層自由水蓄水庫對地下水的日出流系數(shù)；Ki為表層自由水對壤中流的日出流系數(shù)；Cg為地下水消退系數(shù)；Ci為壤中流消退系數(shù)；Cs為河網(wǎng)蓄水消退系數(shù)；X為馬斯京根演算參數(shù)；KE為馬斯京根法演算參數(shù)；L為滯時。

表1 日模型參數(shù)上下邊界

注：KE=ΔT1，取時段間隔，L=0。

表2 次洪模型參數(shù)上下邊界

3.2 目標函數(shù)

本文選擇確定性系數(shù)和平均各場次洪峰相對誤差作為目標函數(shù)，對應(yīng)的方程分別是：

(1)

(2)

式中yoi、yci——實測流量和模擬流量序列；n——實測資料的長度。

4 應(yīng)用實例

本文將PA-DDS算法收斂性能與NSGA-Ⅱ算法和AMALGAM算法進行對比，并將非劣解分布的均勻性及解的相似性方面與AMALGAM進行比較，選取深圳市西麗水庫入庫流量過程數(shù)據(jù)進行分析。

4.1 入庫洪水還原計算

西麗水庫無實測流量資料，實測的入庫洪水過程可以采用水量平衡方程[15]進行計算：

(3)

上式可以寫為

(4)

根據(jù)西麗水庫逐日/時水位、降雨監(jiān)測數(shù)據(jù)及引調(diào)水資料，根據(jù)式4可還原計算出等式左邊數(shù)值即為日/次洪模型入庫流量作為新安江模型的真值系列。

4.2算法收斂性能分析

基于超體積(Hyper-Volume)在評價Pareto前端收斂性、寬廣性和均勻性時表現(xiàn)出來的客觀性，常被用來評價多目標問題的搜索結(jié)果，本文選擇超體積作為指標來判斷PA-DDS算法與NSGA-Ⅱ算法和AMALGAM算法的求解質(zhì)量，評價各種算法的收斂性能。其定義如下：

(5)

式中NPF——最終Pareto前沿上解的個數(shù)；νi——Pareto前沿上第i個解與參照點圍成的體積。

對于目標函數(shù)最小的雙目標優(yōu)化問題，超體積值越大，表明Pareto前沿解集廣度越大，因此可以采用超體積作為衡量算法所得非劣解集優(yōu)劣性的判別標準。取西麗水庫2014年入庫洪水過程的新安江模型采用3種優(yōu)化算法求目標函數(shù)最優(yōu)解，并將求解過程中的迭代次數(shù)及計算超體積特征值繪制在圖3中。從圖3中可以看出，由于PA-DDS算法初始種群數(shù)量小于AMALGAM算法與NSGA-Ⅱ算法，所以在優(yōu)化過程初期，PA-DDS算法超體積值略小，但PA-DDS算法達到收斂時的超體積值大于AMALGAM算法和NSGA-Ⅱ算法，表明其最終得到的解集廣度更廣且更加接近理論最優(yōu)Pareto前沿，因此PA-DDS算法收斂性能較另外2種算法好。

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

4.3.1日模型優(yōu)化結(jié)果

從圖3中可以看出，PA-DDS算法和AMALGAM算法在達到收斂時兩者對應(yīng)的超體積值相對接近，所以將兩種算法求解的目標函數(shù)值進行對比分析。選取2003—2013年共10 a的西麗水庫實測入庫流量過程作為日模型參數(shù)率定序列，將兩種算法的Pareto前沿繪制在圖4中。

由于PA-DDS算法在求解過程中會動態(tài)存儲所有非劣解，所以其Pareto前沿的廣度和均勻性要優(yōu)于AMALGAM算法。同時在目標函數(shù)同時達到最小值的區(qū)域附近，PA-DDS算法非劣解分布更加密集且均勻，從圖中可以看出來Pareto曲線為下凸曲線，有明顯的拐點，從而滿足獲得最優(yōu)解的條件，其算法表現(xiàn)相對更優(yōu)。

將PA-DDS算法所得Pareto前沿下凸拐點時的參數(shù)組作為最優(yōu)參數(shù)組，結(jié)果見表3。采用2014—2016年共3 a的資料用于參數(shù)檢驗。新安江模型日模平均相對誤差為9.96%，平均確定性系數(shù)為0.75(表4)。本文羅列部分日模型實測及模擬過程線，見圖5。

表3 日模型優(yōu)化參數(shù)

表4 日模型模擬效果

4.3.2次洪模型優(yōu)化結(jié)果

選取40場洪水進行參數(shù)率定，10場洪水進行參數(shù)檢驗。將日模型最優(yōu)參數(shù)組合中的K、B、C、Wm、Wum、Wlm、Im、Ex等8個參數(shù)直接運用在次洪模型中，次洪模型只需優(yōu)化計算剩余的7個參數(shù)。將PA-DDS和AMALGAM兩種算法的Pareto前沿繪制在圖6中，從圖中可以看出，2種算法的Pareto分布特點和圖4相類似。將PA-DDS算法所得Pareto前沿下凸拐點時的參數(shù)組作為最優(yōu)參數(shù)組，結(jié)果見表5。率定期和檢驗期洪水擬合精度統(tǒng)計情況見表6。從圖7中可以看出，次洪模型模擬流量和實測流量擬合較好，平均確定性系數(shù)達到近0.84，部分確定性系數(shù)達到0.97。本文羅列部分次洪實測及模擬過程線，見圖7。

表5 次洪模型優(yōu)化參數(shù)

表6 次洪模型模擬效果

5 結(jié)論

本文將PA-DDS 、NSGA-Ⅱ、AMALGAM 3種優(yōu)化算法進行比較優(yōu)選，優(yōu)選結(jié)果表明，PA-DDS多目標優(yōu)化算法在收斂性能方面優(yōu)于NSGA-Ⅱ算法和AMALGAM算法，PA-DDS多目標優(yōu)化算法在非劣解分布的均勻性及解的相似性方面優(yōu)于AMALGAM。

本文首次利用多目標優(yōu)化方法，采用PA-DDS算法對三水源新安江模型參數(shù)進行優(yōu)化，以深圳市西麗水庫入庫流量過程資料進行分析，結(jié)果表明不論日模型和次洪模型，PA-DDS算法都表現(xiàn)出尋優(yōu)速度快，非劣解分布范圍廣并且穩(wěn)定等特點，表明PA-DDS針對新安江模型洪水預(yù)報這種多目標多參數(shù)優(yōu)化中優(yōu)勢明顯，并且模型模擬效果較好，可為西麗水庫洪水預(yù)報模型提供技術(shù)支撐，為防洪與供水調(diào)度提供理論支撐，通過本文的應(yīng)用實踐，說明PA-DDS方法能夠提高新安江模型的參數(shù)率定效率，為其參數(shù)優(yōu)化提供新的思路方法，具有很好的借鑒意義，可以在其它流域推廣應(yīng)用。