湯 嶺,王海軍,李致家,黃迎春,盛奕華

(1.河海大學(xué)水文水資源學(xué)院,江蘇 南京 210098; 2.山東省水文中心,山東 濟(jì)南 250013)

隨著計(jì)算機(jī)技術(shù)的發(fā)展,水文模型發(fā)展迅速,從單一模型不斷優(yōu)化改進(jìn)為綜合模型[1-3],再到如今廣泛使用的基于物理基礎(chǔ)的分布式水文模型[4-7]。分布式水文模型大多以DEM網(wǎng)格為計(jì)算單位,參數(shù)均具有明確的物理意義,每個(gè)網(wǎng)格代表的實(shí)際下墊面情況均可以通過(guò)參數(shù)反映[8],水流運(yùn)動(dòng)一般通過(guò)連續(xù)方程和動(dòng)力方程求解,能較全面地描述水文過(guò)程。但由于其計(jì)算量大、參數(shù)數(shù)量多[9-10],模型需要較長(zhǎng)時(shí)間的模擬,且相比單目標(biāo)優(yōu)化[11-13],多目標(biāo)優(yōu)化更是需要長(zhǎng)時(shí)間運(yùn)行模型才能找到全局最優(yōu)解。因此,分布式水文模型參數(shù)優(yōu)化的時(shí)效問(wèn)題成為了熱點(diǎn)和難點(diǎn)問(wèn)題[14-16],制約了模型的實(shí)用性。

基于代理建模的優(yōu)化方法是一種減少參數(shù)優(yōu)化計(jì)算量的有效方法,已廣泛應(yīng)用于水文建模和水資源管理中[17]。近年來(lái),出現(xiàn)了許多基于代理建模的多目標(biāo)優(yōu)化方法,其中一些是基于經(jīng)典的多目標(biāo)非支配排序,使用廉價(jià)代理模型替代昂貴的動(dòng)態(tài)模型。例如:Syberfeldt等[18]提出了多目標(biāo)并行代理輔助進(jìn)化算法(multi-objective parallel substitute-assisted evolution algorithm, MOPSA-EA),該算法使用一個(gè)ANN作為代理,考慮了代理計(jì)算的不確定性和并行性;宋曉猛等[19]使用響應(yīng)曲面對(duì)新安江模型進(jìn)行參數(shù)優(yōu)化,認(rèn)為代理模型可有效降低計(jì)算效率且穩(wěn)定性較好;Kourakos等[20]開(kāi)發(fā)了一種基于模塊化神經(jīng)網(wǎng)絡(luò)代理的多目標(biāo)優(yōu)化算法MOSA (multi-objective surrogate assisted algorithm ),在希臘圣托里尼的含水層管理中應(yīng)用效果很好。

為解決復(fù)雜水文模型參數(shù)優(yōu)化的時(shí)效問(wèn)題,本文以沂河下游臨沂水文站以上流域(臨沂流域)為研究區(qū)進(jìn)行洪水模擬,將多目標(biāo)自適應(yīng)代理模型優(yōu)化(multi-objective adaptive surrogate modeling-based optimization,MO-ASMO)算法應(yīng)用在TOPKAPI(topographic kinematic approximation and integration)模型的參數(shù)率定中,并與傳統(tǒng)多目標(biāo)優(yōu)化算法NSGA-Ⅱ(non-dominated sorting genetic algorithm Ⅱ)、NSGA-Ⅲ進(jìn)行比較,以期提高TOPKAPI模型在濕潤(rùn)半濕潤(rùn)地區(qū)的模擬精度。

1 研究方法

1.1 TOPKAPI模型

TOPKAPI模型由Todini[21]基于運(yùn)動(dòng)學(xué)和流域地形學(xué)結(jié)合的理念構(gòu)建,該模型將研究區(qū)域細(xì)分為多個(gè)網(wǎng)格,并假設(shè)當(dāng)前計(jì)算網(wǎng)格只與下游唯一的網(wǎng)格相連接,通過(guò)網(wǎng)格之間的高程差計(jì)算徑流路線(xiàn)。空間上采用非線(xiàn)性運(yùn)動(dòng)波方程模擬水流運(yùn)動(dòng),以3個(gè)結(jié)構(gòu)上相似的非線(xiàn)性水庫(kù)方程分別描述土壤、地表和河道中的水流運(yùn)動(dòng)規(guī)律,從而得到整個(gè)流域的水文過(guò)程。

1.2 NSGA-Ⅱ、NSGA-Ⅲ

NSGA-Ⅱ是Deb等[22]于2002年提出的多目標(biāo)優(yōu)化算法,由NSGA改進(jìn)而來(lái),降低了計(jì)算復(fù)雜度,并能保持種群的多樣性;引入了精英策略,擴(kuò)大了采樣空間,防止最佳個(gè)體的丟失。具體流程如圖1所示。

NSGA-Ⅱ在求解目標(biāo)維數(shù)較低的問(wèn)題時(shí)較為有效,但面對(duì)三維及以上的高維目標(biāo)優(yōu)化問(wèn)題效果不理想[23]。NSGA-Ⅲ的整體流程與NSGA-Ⅱ一致,皆通過(guò)交叉變異和非支配排序得到Pareto前沿,兩者的區(qū)別在保持種群多樣性方面,NSGA-Ⅱ是利用擁擠度比較操作,NSGA-Ⅲ則通過(guò)參考點(diǎn)機(jī)制。NSGA-Ⅲ具體流程如圖2所示。

1.3 MO-ASMO算法

MO-ASMO算法可以滿(mǎn)足復(fù)雜的基于物理基礎(chǔ)的分布式水文模型的參數(shù)優(yōu)化[24],主要有3個(gè)步驟:

a.選擇最合適的初始采樣和代理建模算法。用好格子點(diǎn)法[25]作為初始采樣方法,高斯回歸(GPR)[26]作為代理建模算法。GPR是一種具有多種協(xié)方差函數(shù)的靈活回歸方法,可以適應(yīng)不同的插值和回歸任務(wù)。此外,超參數(shù)值還可以控制GPR的行為。例如,如果噪聲項(xiàng)設(shè)置為0,則GPR可視為多維插值方法。如果特征長(zhǎng)度較大,則GPR是平滑的,如果較小,則GPR對(duì)輸出的微小變化很敏感。為了選擇合適的超參數(shù)值,使用SCE-UA算法最大化邊緣似然函數(shù)。

b.使用最具代表性的Pareto子集進(jìn)行模擬。這種自適應(yīng)抽樣策略可以簡(jiǎn)單而有效地平衡收斂性和多樣性。

c.通過(guò)設(shè)置MO-ASMO算法參數(shù)調(diào)整模型運(yùn)行終止條件,本文設(shè)置算法參數(shù)為運(yùn)行迭代次數(shù)。

MO-ASMO算法具體流程如下:①采用GLP抽樣生成初始待率定參數(shù)矩陣,放入TOPKAPI模型生成模擬流量過(guò)程線(xiàn)并計(jì)算目標(biāo)函數(shù)矩陣。②將參數(shù)矩陣和目標(biāo)函數(shù)矩陣輸入GPR代理模型,得到相關(guān)關(guān)系。③原MO-ASMO算法用NSGA-Ⅱ?qū)Υ砟Ｐ瓦M(jìn)行多目標(biāo)優(yōu)化,本文考慮到有3個(gè)目標(biāo)函數(shù),將其改進(jìn)使用NSGA-Ⅲ對(duì)代理模型優(yōu)化得到Pareto解集,取20%的Pareto最優(yōu)解放入TOPKAPI模型運(yùn)行并生成目標(biāo)函數(shù)矩陣。④利用新生成參數(shù)矩陣和目標(biāo)函數(shù)矩陣更新訓(xùn)練數(shù)據(jù),再次構(gòu)建代理模型,在滿(mǎn)足回歸迭代次數(shù)條件前,重復(fù)步驟②～④。⑤輸出最終Pareto前沿對(duì)應(yīng)的參數(shù)樣本和目標(biāo)函數(shù)值。

1.4 Morris法

Morris法是Morris[27]于1991年提出的一種全局敏感性分析法, 通過(guò)獨(dú)立隨機(jī)的“一次變化法”抽樣,再評(píng)價(jià)每個(gè)參數(shù)值的改變對(duì)目標(biāo)函數(shù)的影響,從而可以定性分析出敏感參數(shù)。

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

2.1 研究區(qū)概況

臨沂流域面積為10316km2,呈扇形,屬溫帶大陸性季風(fēng)氣候。該流域雨量分布時(shí)空不均,夏秋多,冬春少,年降水量約800mm;最高氣溫36.5℃,最低氣溫-11℃,全年平均氣溫14.1℃,沂河上興建了大量的水利攔河工程,有跋山水庫(kù)、唐村水庫(kù)、岸堤水庫(kù)、許家崖水庫(kù)4個(gè)大型水庫(kù),4個(gè)水庫(kù)控制區(qū)域產(chǎn)生的徑流作為入流處理。

2.2 基礎(chǔ)數(shù)據(jù)及預(yù)處理

TOPKAPI模型需要輸入的數(shù)據(jù)包括DEM、土壤類(lèi)型、氣象數(shù)據(jù)等。本文DEM (圖3)來(lái)自地理空間數(shù)據(jù)云SRTM數(shù)據(jù),分辨率為1 km;土壤類(lèi)型數(shù)據(jù)(圖4(a))來(lái)自世界土壤數(shù)據(jù)庫(kù),分辨率為1 km;土地利用數(shù)據(jù)(圖4(b))來(lái)自中國(guó)科學(xué)院資源環(huán)境科學(xué)數(shù)據(jù)中心,采用中國(guó)土地利用現(xiàn)狀遙感監(jiān)測(cè)數(shù)據(jù)庫(kù)數(shù)據(jù)集,分辨率為1 km;通過(guò)Strahler方法(圖4(c))對(duì)河道分級(jí);氣溫?cái)?shù)據(jù)來(lái)源于中國(guó)氣象數(shù)據(jù)網(wǎng)(http://data.cma.cn),采用中國(guó)地面氣候資料日值數(shù)據(jù)集(V3.0);降水量、實(shí)測(cè)流量和實(shí)測(cè)入流等數(shù)據(jù)由淮河水文局提供,時(shí)間跨度為2011—2021年,時(shí)段間隔為2h。

圖4 臨沂流域土壤類(lèi)型、土地利用分布和Strahler分級(jí)Fig.4 Soil type, land use distribution and Strahler grading chart in the Linyi Basin

2.3 評(píng)價(jià)指標(biāo)

在半濕潤(rùn)半干旱的淮河流域,洪水過(guò)程常常歷時(shí)短、漲水退水過(guò)程迅速,因此選擇洪峰流量、峰現(xiàn)時(shí)間和納什效率系數(shù)為主要評(píng)價(jià)指標(biāo)。

洪峰流量誤差REp表示單場(chǎng)洪水過(guò)程實(shí)測(cè)洪峰與模擬洪峰的相對(duì)誤差絕對(duì)值,越小表明預(yù)報(bào)洪峰流量偏離實(shí)測(cè)洪峰流量的幅度越小。峰現(xiàn)時(shí)間誤差TEp表示單場(chǎng)洪水過(guò)程模擬的峰現(xiàn)時(shí)間與實(shí)測(cè)峰現(xiàn)時(shí)間的絕對(duì)誤差;納什效率系數(shù)NSE表示單場(chǎng)洪水過(guò)程模擬系列與實(shí)測(cè)系列之間擬合程度,越接近1,表明模擬系列與實(shí)測(cè)系列整體擬合越好。設(shè)MREp、MTEp、MNSE分別為多場(chǎng)洪水的平均洪峰流量誤差、峰現(xiàn)時(shí)間誤差和納什效率系數(shù)。

2.4 Morris法分析結(jié)果

TOPKAPI模型參數(shù)眾多[28],優(yōu)化不敏感參數(shù)會(huì)降低率定效率。因此本文采用Morris法對(duì)23個(gè)TOPKAPI模型參數(shù)抽樣100次,對(duì)得到的2300組樣本進(jìn)行敏感性分析,選取率定期10場(chǎng)洪水的各參數(shù)敏感性平均值作為該參數(shù)的整體敏感性,敏感性分析結(jié)果見(jiàn)表1。

表1 敏感性分析結(jié)果Table 1 Sensitivity analysis results

可以看到影響3個(gè)目標(biāo)函數(shù)的敏感參數(shù)具有高度相似性,土壤厚度為最敏感的一類(lèi)參數(shù),其次為匯流參數(shù)中的河道曼寧系數(shù),主要是四級(jí)和五級(jí)河道的曼寧系數(shù)影響較大,土壤橫、縱向水力傳導(dǎo)度對(duì)目標(biāo)函數(shù)也有一定的影響,地表曼寧系數(shù)中存在部分類(lèi)別對(duì)峰現(xiàn)時(shí)間誤差有一定的影響。綜合上述結(jié)果,選取影響最大的12個(gè)參數(shù)進(jìn)行參數(shù)率定,分別為L(zhǎng)1、L2、ksv1、ksv2、n_c1、n_c2、n_c3、n_c4、n_c5、n_o2、n_o3、n_o5。

2.5 參數(shù)率定結(jié)果

從臨沂流域2011—2021年資料中選取15場(chǎng)洪水進(jìn)行模擬分析,其中前10場(chǎng)洪水為率定期,后5場(chǎng)為驗(yàn)證期?？紤]到計(jì)算機(jī)的運(yùn)行效率,3種多目標(biāo)優(yōu)化方法均驅(qū)動(dòng)原始模型(TOPKAPI模型)5400次,不同參數(shù)優(yōu)化方法得到Pareto前沿如圖5所示。

圖5 NSGA-Ⅱ、NSGA-Ⅲ和MO-ASMO算法的Pareto前沿Fig.5 Pareto frontiers of NSGA-Ⅱ, NSGA-Ⅲ and MO-ASMO

由圖5～7和表2知,驗(yàn)證期MNSE效果普遍比率定期好,原因可能是率定期的10場(chǎng)洪水中出現(xiàn)個(gè)別模擬較差的場(chǎng)次,從而拉低了率定期的MNSE,驗(yàn)證期MREp與MTEp比率定期稍差。驗(yàn)證期MO-ASMO算法和NSGA-Ⅲ表現(xiàn)各有優(yōu)劣,MTEp評(píng)價(jià)指標(biāo)方面,NSGA-Ⅲ得到更好更集中的Pareto前沿,而MNSE和MREp方面,MO-ASMO算法表現(xiàn)更好。綜合來(lái)看,在TOPKAPI模型中,MO-ASMO算法比NSGA-Ⅱ和NSGA-Ⅲ能取得更好效果。NSGA-Ⅲ和MO-ASMO算法各項(xiàng)指標(biāo)均優(yōu)于NSGA-Ⅱ,原因是NSGA-Ⅲ增加了參考點(diǎn)機(jī)制,從而導(dǎo)致樣本點(diǎn)分布更分散,自變量多樣性增加,從而更容易得到全局最優(yōu)解。

圖6 率定期3種算法目標(biāo)函數(shù)比較Fig.6 Comparison of objective functions of three algorithms in calibration period

圖7 驗(yàn)證期3種算法目標(biāo)函數(shù)比較Fig.7 Comparison of objective functions of three algorithms in verification period

表2 3種算法目標(biāo)函數(shù)對(duì)比Table 2 Comparison of objective functions of three algorithms

考慮到3種多目標(biāo)優(yōu)化算法得到Pareto解集個(gè)數(shù)不一致,本文對(duì)每個(gè)Pareto解集進(jìn)行歸一化,再求其最小歐幾里得距離得到一組參數(shù)作為各個(gè)方法的最終參數(shù),運(yùn)行TOPKAPI模型得到不同方法每場(chǎng)洪水的目標(biāo)函數(shù)值(表2)。率定期10場(chǎng)洪水,NSGA-Ⅱ、NSGA-Ⅲ和MO-ASMO算法的洪峰合格率分別為50%、70%和60%,峰現(xiàn)時(shí)間合格率為60%,70%和70%;驗(yàn)證期5場(chǎng)洪峰合格率分別為40%,40%和60%,峰現(xiàn)時(shí)間合格率為40%,20%和60%,MO-ASMO算法在率定期和驗(yàn)證期對(duì)洪峰誤差的控制較好,峰現(xiàn)時(shí)間誤差也合格,綜合看比其他2種方法更優(yōu)。圖8為部分洪水實(shí)測(cè)與模擬結(jié)果對(duì)比,3種不同率定方法下各場(chǎng)洪水模擬與實(shí)測(cè)過(guò)程線(xiàn)在低水階段較吻合,高水處出現(xiàn)不同程度偏差,總體來(lái)說(shuō)MO-ASMO算法結(jié)果較好,其中2012070700號(hào)洪水第一個(gè)峰3種方法模擬均偏大,原因?yàn)門(mén)OPKAPI模型退水時(shí)間較長(zhǎng),洪水后期的降雨直接轉(zhuǎn)化為徑流,從而導(dǎo)致峰值模擬偏大。

圖8 部分洪水實(shí)測(cè)與模擬結(jié)果對(duì)比Fig.8 Comparison of partial measured floods and simulated results

3 結(jié) 語(yǔ)

為提高復(fù)雜水文模型參數(shù)優(yōu)化效率,通過(guò)Morris敏感性分析方法成功篩選出12個(gè)敏感性參數(shù),減少了參與率定的參數(shù)個(gè)數(shù)。在多目標(biāo)參數(shù)優(yōu)化方面,將多目標(biāo)自適應(yīng)代理模型優(yōu)化(MO-ASMO)算法應(yīng)用在TOPKAPI模型的參數(shù)率定中,對(duì)比NSGA-Ⅱ和NSGA-Ⅲ算法,MO-ASMO算法在相同模型運(yùn)行次數(shù)下產(chǎn)生更好的Pareto前沿。通過(guò)分析相對(duì)最優(yōu)解的模擬結(jié)果,MO-ASMO算法在各項(xiàng)指標(biāo)均取得較高的合格率,有效提高了模型參數(shù)優(yōu)化的效率,驗(yàn)證了該方法的適用性。但本文僅從模型參數(shù)率定方面探究其對(duì)模擬結(jié)果的影響,實(shí)際流域的情況比模擬復(fù)雜許多,下一步的研究中,需要考慮流域的實(shí)時(shí)變化,深入探討代理模型在水文模型中參數(shù)實(shí)時(shí)校正等問(wèn)題。