代理模型參數率定方法在TOPKAPI模型中的應用

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

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

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

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

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

1 研究方法

1.1 TOPKAPI模型

TOPKAPI模型由Todini[21]基于運動學和流域地形學結合的理念構建,該模型將研究區域細分為多個網格,并假設當前計算網格只與下游唯一的網格相連接,通過網格之間的高程差計算徑流路線?？臻g上采用非線性運動波方程模擬水流運動,以3個結構上相似的非線性水庫方程分別描述土壤、地表和河道中的水流運動規律,從而得到整個流域的水文過程。

1.2 NSGA-Ⅱ、NSGA-Ⅲ

NSGA-Ⅱ是Deb等[22]于2002年提出的多目標優化算法,由NSGA改進而來,降低了計算復雜度,并能保持種群的多樣性;引入了精英策略,擴大了采樣空間,防止最佳個體的丟失。具體流程如圖1所示。

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

1.3 MO-ASMO算法

MO-ASMO算法可以滿足復雜的基于物理基礎的分布式水文模型的參數優化[24],主要有3個步驟:

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

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

c.通過設置MO-ASMO算法參數調整模型運行終止條件,本文設置算法參數為運行迭代次數。

MO-ASMO算法具體流程如下:①采用GLP抽樣生成初始待率定參數矩陣,放入TOPKAPI模型生成模擬流量過程線并計算目標函數矩陣。②將參數矩陣和目標函數矩陣輸入GPR代理模型,得到相關關系。③原MO-ASMO算法用NSGA-Ⅱ對代理模型進行多目標優化,本文考慮到有3個目標函數,將其改進使用NSGA-Ⅲ對代理模型優化得到Pareto解集,取20%的Pareto最優解放入TOPKAPI模型運行并生成目標函數矩陣。④利用新生成參數矩陣和目標函數矩陣更新訓練數據,再次構建代理模型,在滿足回歸迭代次數條件前,重復步驟②～④。⑤輸出最終Pareto前沿對應的參數樣本和目標函數值。

1.4 Morris法

Morris法是Morris[27]于1991年提出的一種全局敏感性分析法, 通過獨立隨機的“一次變化法”抽樣,再評價每個參數值的改變對目標函數的影響,從而可以定性分析出敏感參數。

2 實例應用

2.1 研究區概況

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

2.2 基礎數據及預處理

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

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

2.3 評價指標

在半濕潤半干旱的淮河流域,洪水過程常常歷時短、漲水退水過程迅速,因此選擇洪峰流量、峰現時間和納什效率系數為主要評價指標。

洪峰流量誤差REp表示單場洪水過程實測洪峰與模擬洪峰的相對誤差絕對值,越小表明預報洪峰流量偏離實測洪峰流量的幅度越小。峰現時間誤差TEp表示單場洪水過程模擬的峰現時間與實測峰現時間的絕對誤差;納什效率系數NSE表示單場洪水過程模擬系列與實測系列之間擬合程度,越接近1,表明模擬系列與實測系列整體擬合越好。設MREp、MTEp、MNSE分別為多場洪水的平均洪峰流量誤差、峰現時間誤差和納什效率系數。

2.4 Morris法分析結果

TOPKAPI模型參數眾多[28],優化不敏感參數會降低率定效率。因此本文采用Morris法對23個TOPKAPI模型參數抽樣100次,對得到的2300組樣本進行敏感性分析,選取率定期10場洪水的各參數敏感性平均值作為該參數的整體敏感性,敏感性分析結果見表1。

表1 敏感性分析結果Table 1 Sensitivity analysis results

可以看到影響3個目標函數的敏感參數具有高度相似性,土壤厚度為最敏感的一類參數,其次為匯流參數中的河道曼寧系數,主要是四級和五級河道的曼寧系數影響較大,土壤橫、縱向水力傳導度對目標函數也有一定的影響,地表曼寧系數中存在部分類別對峰現時間誤差有一定的影響。綜合上述結果,選取影響最大的12個參數進行參數率定,分別為L1、L2、ksv1、ksv2、n_c1、n_c2、n_c3、n_c4、n_c5、n_o2、n_o3、n_o5。

2.5 參數率定結果

從臨沂流域2011—2021年資料中選取15場洪水進行模擬分析,其中前10場洪水為率定期,后5場為驗證期?？紤]到計算機的運行效率,3種多目標優化方法均驅動原始模型(TOPKAPI模型)5400次,不同參數優化方法得到Pareto前沿如圖5所示。

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

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

圖6 率定期3種算法目標函數比較Fig.6 Comparison of objective functions of three algorithms in calibration period

圖7 驗證期3種算法目標函數比較Fig.7 Comparison of objective functions of three algorithms in verification period

表2 3種算法目標函數對比Table 2 Comparison of objective functions of three algorithms

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

圖8 部分洪水實測與模擬結果對比Fig.8 Comparison of partial measured floods and simulated results

3 結 語

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