杜彥臻，劉珈伊，孫夢瑤，林洪孝，王 剛

(山東農(nóng)業(yè)大學(xué)水利土木工程學(xué)院，山東 泰安 271018)

準(zhǔn)確和可靠的徑流預(yù)報對于水資源的優(yōu)化管理十分重要。由于水文過程的非線性和非均質(zhì)性，通過水文模型校準(zhǔn)模擬時，主要在于水文模型參數(shù)的優(yōu)化選擇，通過優(yōu)選模型參數(shù)使校準(zhǔn)模擬系統(tǒng)與實際徑流特征相匹配。由于水文預(yù)報的降雨—徑流模型受到多種因素驅(qū)動，且預(yù)報模型的復(fù)雜性各不相同，大部分參數(shù)都需要估值，盡管近年來基于模型預(yù)報取得了進(jìn)展，但水文學(xué)家普遍認(rèn)為參數(shù)必須進(jìn)行一定程度的校準(zhǔn)才能較好獲得對觀測流量的可接受模擬。水文模型參數(shù)校準(zhǔn)主要有人工試錯法和自動優(yōu)化法，由于傳統(tǒng)手動校準(zhǔn)方法具有主觀性、費時性和不確定性，并且難以實現(xiàn)多個性能目標(biāo)(比如高流量、低流量、平均流量)的優(yōu)選。隨著先進(jìn)的水文預(yù)報系統(tǒng)和計算機(jī)技術(shù)的發(fā)展，大量的水資源自動校準(zhǔn)方法將其取代，如遺傳算法、SCE-UA算法[1]、單純形法、改進(jìn)算法[2]、Rosenbrock法[3]等。盡管傳統(tǒng)的單目標(biāo)自動校準(zhǔn)方法快速客觀，但單個標(biāo)準(zhǔn)進(jìn)行模擬不能全面考慮水文過程的綜合特征，從而未得到水文學(xué)家的廣泛接受。

考慮多目標(biāo)情況對水文模型能更準(zhǔn)確的校準(zhǔn)模擬，國內(nèi)外進(jìn)行了大量研究，主要有：①MOCOM-UA算法[4]是Yaop等(1998年)開發(fā)的將復(fù)雜混洗策略和單純形演化與Pareto優(yōu)勢概念相結(jié)合，解決多目標(biāo)校準(zhǔn)問題的復(fù)合演化全局優(yōu)化算法；②MOSCDE算法[5,6]是郭俊等(2013年)基于SCE-UA上利用種群進(jìn)化的相關(guān)經(jīng)驗，采用全局優(yōu)化和DE算法為核心提高算法自我調(diào)節(jié)能力，進(jìn)一步提高收斂速度和計算效率；③NSGA-II算法[7]是Srinivas和Deb(2000年)在NSGA基礎(chǔ)上提出，應(yīng)用快速非支配排序算法和擁擠度比較算子，復(fù)雜度比NSGA低，提高了優(yōu)化結(jié)果的精度；④MOPSO算法是基于群體智能遷移協(xié)作提出的搜索算法，能夠考慮多目標(biāo)對系統(tǒng)參數(shù)進(jìn)行有效優(yōu)化；⑤MOSCEM-UA算法是由阿姆斯特丹和亞利桑那大學(xué)合作開發(fā)的一種基于全局優(yōu)化的水文模型參數(shù)校準(zhǔn)方法。因此本文通過應(yīng)用MOPSO、MOSCEM-UA兩種算法對沁水河流域進(jìn)行模型參數(shù)多目標(biāo)優(yōu)化研究具有十分重要的意義。

1 研究方法

1.1 MOPSO算法

PSO算法是基于種群內(nèi)鳥類社會行為模擬的搜索算法。在大自然中，鳥類通過考慮個人經(jīng)驗和種群移動方向來尋求食物。受種群智能啟發(fā)促使肯尼迪和Eberhart[8]提出了PSO算法。PSO的簡單性和適用性增加了在解決大量工程和管理優(yōu)化問題方面的普及，可用于解決大規(guī)模、非線性、不可微和多峰值的復(fù)雜優(yōu)化問題，并發(fā)現(xiàn)在處理單目標(biāo)優(yōu)化方面非常有效。簡單高效的PSO算法促使研究人員在2002年以來將其應(yīng)用于解決MO問題，進(jìn)而研制出處理多目標(biāo)函數(shù)的MOPSO算法。

多目標(biāo)粒子群優(yōu)化算法(MOPSO)有兩個屬性：①在MO問題Pareto前沿生成高質(zhì)量的全局非支配解；②尋找Pareto前沿解中適應(yīng)多樣性的單個粒子最佳解。MOPSO算法在慣性權(quán)重w和加速度系數(shù)c1和c2方面有較合理的適應(yīng)性，使其能在搜索空間的探索和利用之間取得良好的平衡，能夠較合理且高效的解決多目標(biāo)問題。但因其受速度和動態(tài)調(diào)節(jié)影響，易陷入局部最優(yōu)，導(dǎo)致收斂精度低，從而導(dǎo)致其生成Pareto解的聚類性和均勻分布性較差。

1.2 MOSCEM-UA算法

多目標(biāo)混合復(fù)雜演化大都會(MOSCEM-UA)算法[9,10]是Vrugt等在SCE-UA算法的基礎(chǔ)上結(jié)合馬爾可夫鏈Metropolis強(qiáng)度提出的較新的適應(yīng)性分配方法進(jìn)行模型參數(shù)校準(zhǔn)，是一種穩(wěn)健和有效的全局多目標(biāo)優(yōu)化方法。算法能夠：①客觀的采樣整個參數(shù)空間而不是離散預(yù)定一組參數(shù)值；②針對多個標(biāo)準(zhǔn)進(jìn)行優(yōu)化，從而提供對結(jié)果的優(yōu)選，即多目標(biāo)優(yōu)化識別可犧牲一個標(biāo)準(zhǔn)來改進(jìn)另一標(biāo)準(zhǔn)的情況，該算法迭代更新一組模型輸入?yún)?shù)的同時最小化若干優(yōu)化標(biāo)準(zhǔn)。由所有參數(shù)值組成的折中表面被稱為Pareto前沿，以突出解決方案非唯一性。

MOSCEM-UA算法使用改進(jìn)的適應(yīng)性分配方法，保留種群多樣性，同時其隨機(jī)性防止算法陷入某個單一 “最佳”參數(shù)集的相對較小區(qū)域，從而進(jìn)一步保留了采樣總體的多樣性并生成相當(dāng)均勻的Pareto前沿，進(jìn)行合理高效的搜索參數(shù)空間。但因其在國外流域模型應(yīng)用較多，對國內(nèi)垂向混合產(chǎn)流模型的適應(yīng)性及精度評價等方面有待進(jìn)一步探究應(yīng)用。

1.3 算法參數(shù)設(shè)置

MOPSO算法參數(shù)主要為：種群規(guī)模100，參數(shù)個數(shù)n=15，目標(biāo)函數(shù)3個，加速度系數(shù)c1+c2≤4，慣性權(quán)重w設(shè)置以0.9～0.01隨迭代次數(shù)而減少，最大迭代次數(shù)I=1 000。

MOSCEM-UA算法參數(shù)主要為：參數(shù)個數(shù)n=15，目標(biāo)函數(shù)3個，復(fù)合形個數(shù)p=8，每個復(fù)合形樣本點數(shù)m=2×15+1=31，種群大小s=8 31=248，每個復(fù)合形進(jìn)化代數(shù)L=m/4，后代接受概率比例因子β=1/2，最大迭代次數(shù)I=1 000。

2 多目標(biāo)優(yōu)化模型

2.1 沁水河流域概況

沁水河流域?qū)倌z東半島低山丘陵區(qū)，地勢為南高北地，流域發(fā)源于玉林店鎮(zhèn)大尖崮南山峰及水道鎮(zhèn)西直格莊一帶，流經(jīng)玉林店、大窯、文化、寧海和水道等5處鎮(zhèn)(街)匯入黃海，上游比降為0.2%～0.5%，下游比降為0.2%，流域集水面積為168 km2。沁水河為典型的山溪性雨源型河道，河道流量與降水量變化規(guī)律相一致，且年內(nèi)、年際變化大。枯季流量較小，洪水主要集中在汛期。流域?qū)倥瘻貛О霛駶櫞箨懶约撅L(fēng)氣候，四季變化分明。流域內(nèi)較大支流有3條：張家河長度約10.83 km、屯圈河長度約10.2 km、石頭河長度約10.1 km。該流域有5個雨量站分別是牟平、玉林店、十六里頭、西流莊和徐家疃，1個流量站牟平水文站，1個蒸發(fā)站門樓水庫。沁水河流域內(nèi)水利工程建設(shè)共計有9座小型水庫、6座閘壩和7座塘壩。

2.2 流域垂向混合產(chǎn)流模型

沁水河位于北方半濕潤地區(qū)，單一產(chǎn)流模型不能準(zhǔn)確模擬研究區(qū)產(chǎn)流量，因此模型考慮沁水河特殊水文特征及人類活動對水文循環(huán)過程要素的影響，提出適合該流域的混合產(chǎn)流模型見圖1所示。模型算法中，先假設(shè)較大的參數(shù)不確定性范圍，通過迭代運算不斷減小參數(shù)范圍，使模擬值不斷與實測值接近。該垂向混合產(chǎn)流模型參數(shù)[11]共有15個參數(shù)，參數(shù)意義、取值范圍及初始值見表1所示。

模型結(jié)構(gòu)主要分為4個層次[12]：①產(chǎn)流部分用蓄滿—超滲垂向混合產(chǎn)流模式；②蒸散發(fā)部分采用三層蒸發(fā)模型；③水源部分用EX次拋物線型自由水蓄水曲線劃分地表、壤中和地下徑流3種水源；④匯流部分地面匯流用單位線法、地下匯流用線性水庫法。降水到達(dá)地面，通過空間分布下滲曲線，劃分地面徑流RS和下滲水量FA，落在地面的部分為蓄滿產(chǎn)流，剩余部分為超滲產(chǎn)流或不產(chǎn)流。地面徑流RS根據(jù)改進(jìn)后的格林—安普特[13]下滲公式計算，計算公式為：

(1)

(2)

RS=P-FA

(3)

式中：FM為平均下滲能力，mm；fc為穩(wěn)定下滲率，mm；KF為滲透系數(shù)；WM為平均蓄水容量，mm；W為實際土壤含水量，mm；

圖1 沁水河流域垂向混合產(chǎn)流計算模型Fig.1 Vertical mixed flow calculation model in the Qinshui River Basin

序號符號意義范圍初始值1K蒸散發(fā)折算系數(shù)0.5～1.212WUM流域上層蓄水容量5～20203WLM流域下層蓄水容量60～90604C深層蒸散發(fā)系數(shù)0.01～0.50.125WM流域平均蓄水容量30～3001306B流域蓄水容量分布曲線指數(shù)0.1～0.60.27BF下滲率分布曲線指數(shù)0～21.628KF土壤對下滲率影響系數(shù)0～203.959fc流域?qū)嶋H穩(wěn)定下滲率1～201510KI自由水箱壤中流出流系數(shù)0.1～0.80.511KG自由水箱地下水出流系數(shù)0.1～0.40.412CS地面水線性水庫匯流系數(shù)0.10～0.990.1513CI壤中流線性水庫匯流系數(shù)0.10～0.990.5014CG地下水線性水庫匯流系數(shù)0.10～0.990.9915L滯時0～21

BF為下滲能力空間分布特征參數(shù)；P為扣除蒸發(fā)的降雨量，mm。

地下徑流RR采用蓄滿產(chǎn)流結(jié)構(gòu)計算，計算公式為：

(4)

總產(chǎn)流量計算公式為：

其次，中藥價值鏈方面因素。價值鏈的上、中、下游分別為研發(fā)、生產(chǎn)和流通等環(huán)節(jié)。中藥價值鏈?zhǔn)菄@中藥產(chǎn)品形成的涉及多個行業(yè)的鏈?zhǔn)浇Y(jié)構(gòu)，它是以中藥制造為核心，由研發(fā)、制造、市場銷售及服務(wù)等活動共同構(gòu)成，這些活動經(jīng)過不同階段，將產(chǎn)品和服務(wù)傳遞給最終消費者。當(dāng)前中藥產(chǎn)業(yè)價值鏈條的向上（研發(fā)）、向下（服務(wù)、貿(mào)易、物流）縱向延伸能力差。具體表現(xiàn)為價值鏈上游優(yōu)秀中醫(yī)藥資源的保護(hù)、研發(fā)不足，即繼承與創(chuàng)新不足，在中醫(yī)藥理論創(chuàng)新、標(biāo)準(zhǔn)創(chuàng)新上亟待加強(qiáng)[10]；價值鏈下游中藥流通環(huán)節(jié)，信息化程度較低、可追溯體系尚未建立、品牌建設(shè)需要加強(qiáng)、服務(wù)體系亟待完善等問題。

R=RS+RR

(5)

2.3 多目標(biāo)函數(shù)的優(yōu)化

多目標(biāo)與單目標(biāo)優(yōu)化是通過一定的優(yōu)化算法獲得目標(biāo)函數(shù)的最優(yōu)解，當(dāng)目標(biāo)函數(shù)有兩個或兩個以上時稱為多目標(biāo)優(yōu)化(Multi-objective Optimization Problem, MOP)，其不同于單目標(biāo)為有限解，多目標(biāo)優(yōu)化的解通常是一組均衡解。水文模擬常使用統(tǒng)計理論指標(biāo)，如均方根誤差作為模型校準(zhǔn)的主要目標(biāo)，但統(tǒng)計指標(biāo)不能確保流域整體水文特征的代表性。本文選用三個指標(biāo)：強(qiáng)調(diào)高、低流量的統(tǒng)計指標(biāo)和強(qiáng)調(diào)水平衡的水文指標(biāo)?？傊?，指標(biāo)涵蓋了三個不同且重要的部分，同時允許在統(tǒng)計和水文指標(biāo)之間進(jìn)行比較。

第一目標(biāo)是均方根誤差(RMSE)，強(qiáng)調(diào)擬合水文圖的高流量部分。第二目標(biāo)是變換均方根誤差(TRMSE)，將實測和模擬流量時間序列先通過Box-Cox變換，轉(zhuǎn)換值λ=0.3，再計算變換流量的RMSE，加大小流量權(quán)重以獲得強(qiáng)調(diào)水文圖的低流量部分。第三目標(biāo)是徑流系數(shù)誤差(ROCE)，強(qiáng)調(diào)水平衡的整體精度。計算公式[14]為：

(6)

(7)

(8)

(9)

3 結(jié)果分析

3.1 參數(shù)敏感度分析

研究選用流域1982～1990年共9a的降雨、徑流和蒸散發(fā)量等資料用于流域垂向混合產(chǎn)流模型參數(shù)優(yōu)化率定，1982～1986年為率定期，1987～1990年為驗證期。在校準(zhǔn)周期內(nèi)通過目標(biāo)函數(shù)對模型的15個參數(shù)進(jìn)行優(yōu)化，參數(shù)初始值作為算法的起點。模型參數(shù)眾多且有的參數(shù)細(xì)微變化會導(dǎo)致模擬結(jié)果有較大變化，但有的參數(shù)對模型影響較小，因此要對參數(shù)進(jìn)行敏感性分析。Morris篩選法是一種局部靈敏度分析法，本文采用修正的摩爾斯篩選法對參數(shù)進(jìn)行靈敏度分析，采用自變量固定步長計算取其均值，計算公式[15]為：

(10)

式中：S為靈敏度判別因子；Yi、Yi+1分別為模型i、i+1次輸出值；Y0為參數(shù)率定后結(jié)果初始值；Pi、Pi+1分別為i、i+1次模型運算參數(shù)值對率定參數(shù)后參數(shù)值變化的百分率；n為運行次數(shù)。

影響模型模擬的15個參數(shù)靈敏度計算結(jié)果見表2所示，參數(shù)WUM、WLM為中等靈敏度參數(shù)，WM、L靈敏參數(shù)，BF、KF、fc、CS高靈敏參數(shù)，其余為不靈敏參數(shù)。在參數(shù)優(yōu)化之前進(jìn)行敏感性分析以評估各參數(shù)的相對重要性，對于低靈敏度參數(shù)，可以根據(jù)經(jīng)驗取值且相對固定，故需要對最有影響的參數(shù)進(jìn)行校準(zhǔn)以進(jìn)行更好的模擬。

表2 模型參數(shù)變化靈敏度分析結(jié)果Tab.2 Model parameter change sensitivity analysis result

注：0 ≤|S|<0.05 為不靈敏參數(shù)，0.05 ≤|S|<0.2為中等靈敏參數(shù)，0.2≤|S|<1為靈敏參數(shù)，|S|≥1為高靈敏參數(shù)。

3.2 多目標(biāo)校準(zhǔn)結(jié)果分析

圖2 MOPSO與MOSCEM算法的非劣前沿Fig.2 Pareto Frontiers of MOPSO and MOSCEM Algorithms

圖3為1982-1990年兩種算法模擬過程的月流量偏差，可以看出MOPSO與MOSCEM算法在汛期的月流量百分比偏差(Bias)基本在10%以內(nèi)，達(dá)到較合理的范圍(理想情況Bias<1%)，其中7、8月偏差基本在5%以內(nèi)，擬合較好。從整體月偏差統(tǒng)計來看，MOSCEM算法在大多數(shù)月份的偏差比MOPSO校準(zhǔn)明顯偏小，這種差異在非汛期尤為顯著，對比得出MOSCEM算法旨在更好地估計高流量及匹配峰值流量，同時減少低流量百分比偏差，整體統(tǒng)計更符合模擬過程。

圖3 兩種算法模擬過程的月流量偏差(1982-1990年)Fig.3 Monthly flow deviations of the two simulation processes from 1982 to 1990

圖4為沁水河流域在率定期1982年兩種算法校準(zhǔn)模擬結(jié)果，MOPSO和MOSCEM算法優(yōu)選參數(shù)模擬實測流量的水文圖見圖4(a)、4(c)，模擬過程中實測和模擬流量之間的差異即校準(zhǔn)殘差見圖4(b)、4(d)。由圖看出，兩種算法模擬都能捕獲大部分高流量事件，其中MOPSO在高流量區(qū)比低流量區(qū)的模擬結(jié)果較好，但對于模擬高流量事件上MOSCEM比MOPSO校準(zhǔn)結(jié)果更優(yōu)，且MOSCEM在模擬衰退流量方面有較好的估計。通過表3校準(zhǔn)結(jié)果回歸統(tǒng)計分析進(jìn)一步看出，兩種算法模擬情況下：①MOSCEM算法模擬的MultipleR(相關(guān)關(guān)系)為0.90,比MOPSO的0.76大，說明MOSCEM模擬結(jié)果較好，相關(guān)性較強(qiáng)；②MOSCEM的標(biāo)準(zhǔn)誤差為0.60， MOPSO為1.04，因而MOSCEM算法的擬合效果更優(yōu)。

圖4 率定期1982年模擬結(jié)果Fig.4 Results of simulations from 1982 in calibration period

圖5為驗證期1989年校準(zhǔn)模擬水文圖及殘差分析，圖5(a)、5(c)看出，對于降雨—徑流模型MOPSO算法模擬結(jié)果與MOSCEM相差較小，且都有明顯的提高，兩種算法都略微提高對大流量事件的模擬，同時較好地適應(yīng)了小流量衰退。通過校準(zhǔn)殘差圖5(b)、5(d)可以看出MOSCEM校準(zhǔn)比MOPSO稍好一些，且在匹配整個水平年的大部分實測流量上整體表現(xiàn)較優(yōu)。原因可能是MOSCEM算法能夠?qū)⒏叩土髁康哪繕?biāo)函數(shù)賦與合理權(quán)重，采樣整個空間對汛期與非汛期綜合拆分進(jìn)行模擬，更適應(yīng)沁水河流域的整個流量擬合過程。

表3 1982年校準(zhǔn)結(jié)果回歸統(tǒng)計分析Tab 3 1982 calibration results regression statistical analysis

注：MultipleR：x和y的相關(guān)系數(shù)r，一般在-1～1之間，絕對值越靠近1相關(guān)性越強(qiáng)，越靠近0相關(guān)性越弱；標(biāo)準(zhǔn)誤差：衡量擬合程度的大小，值越小，說明擬合程度越好。

圖5 驗證期1989年模擬結(jié)果Fig.5 Results of simulations from 1989 in validation period

使用兩種算法對前5年數(shù)據(jù)進(jìn)行初始校準(zhǔn)后，通過運行對后4年數(shù)據(jù)周期(1987-1990)進(jìn)行校準(zhǔn)參數(shù)驗證。表4為流域的總體RMSE、Bias和相關(guān)系數(shù)(MultipleR)對模型性能分析對比：①均方根誤差方面：9年校準(zhǔn)周期內(nèi)總體RMSE基本在允許誤差范圍內(nèi)；②百分比偏差：雖在率定期部分年份1982、1985年的MOPSO擬合稍差一些，但總體上理想偏差合格率從率定期60%提高到驗證期87.5%，其中MOPSO算法合格率由40%提高到75%，MOSCEM由80%提高到99%，偏差基本在理想情況內(nèi)；③相關(guān)系數(shù)：MOPSO較MOSCEM算法的相關(guān)性差，其中MOSCEM算法的相關(guān)系數(shù)78%都達(dá)到0.90以上，擬合結(jié)果相關(guān)性較強(qiáng)?？傮w上MOSCEM算法校準(zhǔn)更傾向于匹配整個流量過程，能夠達(dá)到合理質(zhì)量校準(zhǔn)模擬。

表4 MOPSO與MOSCEM算法模型性能對比結(jié)果Tab.4 Model performance comparison results of MOPSO and MOSCEM Algorithm

4 結(jié) 語

本文考慮水利工程對降雨—徑流模型的影響，建立了適應(yīng)于沁水河流域的垂向混合產(chǎn)流模型，同時考慮了多目標(biāo)對模型參數(shù)校準(zhǔn)優(yōu)化，采用兩種廣泛應(yīng)用的多目標(biāo)優(yōu)化算法對模型參數(shù)進(jìn)行優(yōu)選對比。從模擬結(jié)果可以看出。

(1)垂向混合產(chǎn)流模型在沁水河流域具有一定的適應(yīng)性，能夠模擬該流域降雨徑流過程。

(2)兩種算法所得到的模擬效果均較好，其中MOSCEM算法綜合考慮高、低流量模擬，百分比偏差合格率達(dá)到87.5%，且相關(guān)性系數(shù)基本在0.90以上，與實測更接近，說明此算法更適合本模型模擬。

(3)徑流模擬結(jié)果中，算法在汛期模擬誤差偏大，非汛期偏小，考慮主要原因為汛期受水利工程影響較大，導(dǎo)致模擬誤差偏大，而非汛期降雨徑流較小，擬合較為簡易。

(4)由于多目標(biāo)優(yōu)化算法在國內(nèi)水文模型方面應(yīng)用較少，且僅考慮流量部分作為多目標(biāo)，考慮目標(biāo)情況較為單一，沒有綜合土壤濕度、氣候等方面影響，從而多目標(biāo)參數(shù)校準(zhǔn)方面有待進(jìn)一步提高。因此在今后水文模型參數(shù)多目標(biāo)優(yōu)化研究方面，應(yīng)選用適應(yīng)的水文模型，對多目標(biāo)算法進(jìn)一步研究，以便對水資源模擬進(jìn)行合理的評價，提高校準(zhǔn)過程效率。

□