余勝男，張雪，邵年華，呂海強(qiáng)

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

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

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

1 西麗水庫(kù)概況

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

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

2 PA-DDS算法簡(jiǎn)介

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

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

3 優(yōu)化過(guò)程

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

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

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

注：KE=ΔT1，取時(shí)段間隔，L=0。

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

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

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

(1)

(2)

式中yoi、yci——實(shí)測(cè)流量和模擬流量序列；n——實(shí)測(cè)資料的長(zhǎng)度。

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

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

4.1 入庫(kù)洪水還原計(jì)算

西麗水庫(kù)無(wú)實(shí)測(cè)流量資料，實(shí)測(cè)的入庫(kù)洪水過(guò)程可以采用水量平衡方程[15]進(jìn)行計(jì)算：

(3)

上式可以寫(xiě)為

(4)

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

4.2算法收斂性能分析

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

(5)

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

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

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

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

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

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

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

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

表4 日模型模擬效果

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

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

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

表6 次洪模型模擬效果

5 結(jié)論

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

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