王 維，馮忠倫，楊 偉，林洪孝，王 剛，刁艷芳

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

0 引 言

洪水是自然界對(duì)人類的主要威脅之一。洪水預(yù)報(bào)就是在流域上發(fā)生的暴雨及河流上游的水流補(bǔ)給，經(jīng)過水情信息采集、傳輸、存貯、處理等一系列的科學(xué)計(jì)算，最終預(yù)估在流域出口斷面或河流下游水文站即將發(fā)生的實(shí)際洪水過程。做好洪水預(yù)報(bào)，關(guān)鍵是選擇適應(yīng)流域情況的水文模型以及確保模擬結(jié)果精度。前者主要表現(xiàn)在產(chǎn)匯流結(jié)構(gòu)，后者則主要反映在所采用的參數(shù)組上[1]。

模型參數(shù)的確定方法總體可分為人工試錯(cuò)率定和自動(dòng)優(yōu)化率定兩種。人工率定得來的參數(shù)經(jīng)驗(yàn)性較強(qiáng)，具有極大的主觀性且極為耗時(shí)。故在生產(chǎn)實(shí)踐中人們開始尋求參數(shù)求解的程序化實(shí)現(xiàn)。而針對(duì)參數(shù)的不確定性，理解也從最初的局部尋優(yōu)過渡到尋找全局最優(yōu)解[2]，此過程涌現(xiàn)了諸如遺傳算法、粒子群優(yōu)化算法、SCE-UA等優(yōu)秀算法。其中SCE-UA算法是一種結(jié)合了基因算法、單純形算法優(yōu)點(diǎn)，以信息共享和生物演化規(guī)律為基礎(chǔ)的非線性混合算法。能在不需要顯式目標(biāo)函數(shù)或目標(biāo)函數(shù)的偏導(dǎo)數(shù)情況下有效解決高維參數(shù)全局優(yōu)化問題，在應(yīng)用過程中被證明更適應(yīng)各類水文模型求參[3,4]。

本文針對(duì)山東境內(nèi)一般山丘地形的田莊水庫流域，選取了三水源新安江與垂向混合產(chǎn)流兩種模型結(jié)構(gòu)，將SCE-UA優(yōu)化后的參數(shù)帶入，對(duì)典型洪水進(jìn)行模擬與檢驗(yàn)。

1 SCE-UA算法

SCE-UA(Shuffled Complex Evolution Algorithm)是美國亞利桑那大學(xué)的Duan博士等人在原有單純性算法基礎(chǔ)上結(jié)合生物與基因選擇原理后研究出的一種能夠有效解決非線性約束最優(yōu)化問題的算法，為高效篩選多參數(shù)組合的全局最優(yōu)解提供了一種新的思路。其算法流程見圖1。算法中：n為所選參數(shù)組參數(shù)個(gè)數(shù)，m代表所選每個(gè)復(fù)合形頂點(diǎn)個(gè)數(shù)，p為復(fù)合形個(gè)數(shù)，q為每個(gè)子復(fù)合形頂點(diǎn)數(shù)，s代表種群大小，α、β是父代產(chǎn)生的子代個(gè)數(shù)及代數(shù)。SCE-UA方法的主要步驟為[5,6]：

圖1 SCE-UA算法流程圖Fig.1 Flowchartof the SCE-UA algorithm

(1)初始化，選取參與進(jìn)化的復(fù)合形個(gè)數(shù)p(p≥1)，每個(gè)復(fù)合形包含定點(diǎn)數(shù)m(m≥n+1)與樣本數(shù)s；

(2)產(chǎn)生樣本，在可行空間隨機(jī)產(chǎn)生s個(gè)點(diǎn)x1,x2,…,xs并計(jì)算每個(gè)點(diǎn)xi處函數(shù)值fi=f(xi),=1,2,…,x；

(3)樣本點(diǎn)排序，將s個(gè)點(diǎn)以升序排列，并將它們貯存到數(shù)組中D={xi,fi,i=1,…,s}此時(shí)i=1代表了目標(biāo)最小的函數(shù)點(diǎn)；

(4)復(fù)合形群體劃分，將D劃分p個(gè)復(fù)合形分別為A1,…,Ap，每個(gè)復(fù)合形包含m個(gè)點(diǎn)，組成數(shù)組A，使：Ak={xk,j,fk,j|xk,j=xk+p,j-1,fk,j=fk+p,j-1,j=1,2,…,m}；

(5)復(fù)合群體進(jìn)化，根據(jù)競(jìng)爭(zhēng)復(fù)合形演化算法對(duì)每一個(gè)復(fù)合形(Ak,k=1,2,…,p)進(jìn)行演化；

(6)混合復(fù)合形，將A1,…,Ap代替到D中，有新的D={Ak,k=1,…,p}，再次對(duì)D按目標(biāo)函數(shù)的升序進(jìn)行排列；

(7)收斂性判斷，如果滿足收斂條件則停止，否則返回步驟(4)。

2 所選模型結(jié)構(gòu)與參數(shù)

2.1 三水源新安江模型

新安江模型是以蓄滿產(chǎn)流為主的水文預(yù)報(bào)模型[7]。在濕潤與半濕潤地區(qū)有較為廣泛的應(yīng)用[8]。本文所選新安江模型，產(chǎn)流部分用蓄滿產(chǎn)流計(jì)算徑流量，蒸散發(fā)部分按上層、下層和深層三層計(jì)算，水源部分采用自由水蓄水庫劃分為地表、壤中和地下徑流3種，地面匯流采用單位線法，壤中流與地下徑流用線性水庫計(jì)算，河網(wǎng)匯流采用滯后演算法。因所選田莊水庫流域，以水庫以上面積為一個(gè)完整單元進(jìn)行模擬預(yù)報(bào)，前期對(duì)流域洪水由距法分析推得瞬時(shí)單位線參數(shù)n=1.90，k=4.41。

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

混合產(chǎn)流模型是為解決蓄滿與超滲均無法忽略地區(qū)產(chǎn)流問題而提出的。混合產(chǎn)流至今還沒有形成一套獨(dú)立的產(chǎn)流理論[9]，其研究通常以蓄滿產(chǎn)流和超滲產(chǎn)流兩種典型理論為基礎(chǔ)。目前常用的混合產(chǎn)流計(jì)算方法是垂向混合結(jié)構(gòu)法[10]。此方法是把流域產(chǎn)流過程中地面徑流部分按超滲產(chǎn)流，地面以下按蓄滿產(chǎn)流進(jìn)行計(jì)算，具體結(jié)構(gòu)見圖2。

圖2 垂向結(jié)合法結(jié)構(gòu)圖Fig.2 Vertical binding assay structure

按此方法計(jì)算的流域蓄滿與超滲面積比例反應(yīng)在前期土壤含水量與實(shí)際下滲量隨時(shí)間的變化上，可用下式表示：

(1)

式中：γ為蓄滿產(chǎn)流的面積比例系數(shù)；FA為實(shí)際下滲量；Wmm代表流域最大蓄水量；B代表流域蓄水容量分布曲線指數(shù)；a相應(yīng)于初始土壤平均含水量為W時(shí)的縱坐標(biāo)值。

超滲產(chǎn)流計(jì)算部分即地面徑流RS，用格林-安普特下滲曲線[11]進(jìn)行計(jì)算，計(jì)算公式為：

(3)

RS=PE-FA

(4)

式中：FM為流域平均下滲能力；FC為穩(wěn)定下滲率；WM為流域平均蓄水容量；W為流域?qū)嶋H土壤含水量；KF為土壤對(duì)下滲率影響系數(shù)；BF為下滲率分布曲線指數(shù)；PE為扣除雨間蒸發(fā)的降雨量。

蓄滿產(chǎn)流計(jì)算部分即地面以下的徑流(壤中流與地下徑流之和)RR，計(jì)算公式為：

(6)

綜上所述，產(chǎn)流量R=RS+RR，其余部分與新安江模型描述相同。

3 優(yōu)化過程

3.1 確定模型參數(shù)邊界

參數(shù)邊界由模型應(yīng)用情況人為給定，三水源新安江模型應(yīng)用較多，參數(shù)范圍參照參數(shù)物理意義及此方面已有研究成果[8,12]給出；對(duì)于垂向混合產(chǎn)流模型研究成果有限，參照瞿思敏教授在沐河上游青峰嶺水庫研究結(jié)果[10]，對(duì)參數(shù)范圍進(jìn)行調(diào)整。模型參數(shù)邊界選擇見表1。

表1 模型參數(shù)范圍Tab.1 The range of mode parameters

3.2 選取SCE-UA算法參數(shù)

算法內(nèi)部參數(shù)取值，一般取m=2n+1，q=n+1，α=1，β=m，模擬效果較好[13]。p的取值反映收斂的快慢及收斂精度。過小會(huì)降低參數(shù)模擬精度，過大則會(huì)延誤收斂速度，應(yīng)視具體問題進(jìn)行分析。一般地，p值越大，越適宜于高階非線性問題，本文所選三水源新安江模型變化參數(shù)n1=13，垂向混合產(chǎn)流模型n2=16，所選流域統(tǒng)一取p=8[14]。

3.3 建立目標(biāo)函數(shù)

目標(biāo)函數(shù)的選取對(duì)優(yōu)化的最終結(jié)果有很大影響，針對(duì)田莊水庫流域洪水實(shí)際情況，在考慮洪量與洪水過程后，將洪峰峰值作為一個(gè)單獨(dú)的指標(biāo)納入目標(biāo)函數(shù)。對(duì)洪水模擬，選取以下3個(gè)方程之和作為最終目標(biāo)函數(shù)FObj。有：

(10)

式中：QS來表示實(shí)測(cè)洪水?dāng)?shù)值；Qj表示模擬洪水?dāng)?shù)值。方程fObj1用來控制洪量，保證水量平衡；fObj2反應(yīng)洪水過程，用來保證單位時(shí)段內(nèi)模擬流量與實(shí)測(cè)流量變化趨于一致；fObj3用來控制模擬洪峰與實(shí)測(cè)洪峰誤差。

3.4 設(shè)定迭代終止條件

為防止參數(shù)尋優(yōu)過程中出現(xiàn)死循環(huán)，為程序運(yùn)行過程設(shè)置終止條件，設(shè)計(jì)滿足以下兩個(gè)條件中任意一個(gè)，程序停止迭代，輸出參數(shù)組最優(yōu)解：①目標(biāo)函數(shù)在5次循環(huán)后仍無法提高0.01%的精度情況下；②5次循環(huán)后算法仍無法顯著改變參數(shù)值并且無法改善目標(biāo)函數(shù)。

4 案例分析

田莊水庫地處沂河上游，屬暖溫帶半濕潤季風(fēng)氣候，集水面積417 km2。流域均為山區(qū)，多年平均降水量850 mm，多年平均徑流深320 mm。模型輸入部分從流域內(nèi)雙廟、徐家莊等8個(gè)雨量站收集點(diǎn)雨資料轉(zhuǎn)換為面雨量，蒸發(fā)資料移用流域上游鎮(zhèn)后水文站；流量資料選擇田莊水庫流域內(nèi)配套參證水文站朱家莊的數(shù)據(jù)與模擬結(jié)果比較。

表2 模型參數(shù)優(yōu)化結(jié)果Tab.2 Results of parameter optimization

模型參數(shù)按3.1節(jié)所給邊界范圍內(nèi)隨機(jī)產(chǎn)生后納入流程循環(huán)體中優(yōu)選，選取田莊水庫流域1985～2011年共27 a的18場(chǎng)典型洪水進(jìn)行模擬，其中前12場(chǎng)洪水用于參數(shù)率定，后6場(chǎng)洪水用于參數(shù)檢驗(yàn)。參數(shù)優(yōu)化結(jié)果見表2?？梢钥闯?，優(yōu)化后同一參數(shù)在兩種模型內(nèi)取值略有不同，主要體現(xiàn)在垂向混合產(chǎn)流模型的蓄水容量面積曲線指數(shù)B與深層蒸散發(fā)系數(shù)C比三水源新安江模型偏大，這是由于產(chǎn)流部分不盡相同，參數(shù)取值基本能反應(yīng)所代表的物理實(shí)質(zhì)。洪水模擬成果見表3，其中率定期與驗(yàn)證期各選兩場(chǎng)，典型洪水模擬過程見圖3。

從洪峰洪量角度，由表3可以看出：兩種模型模擬的洪峰，洪量誤差大都能夠控制在容許范圍內(nèi)。從誤差的區(qū)間分布上來看，正負(fù)分布相對(duì)平衡，說明模擬的系統(tǒng)誤差較小。按洪峰洪量不超過實(shí)測(cè)20%來評(píng)判，模擬的18場(chǎng)洪水中，三水源新安江模型洪量相對(duì)誤差有3場(chǎng)不合格，全部集中在率定期，洪量合格率83.3%，洪峰相對(duì)誤差在率定期有1場(chǎng)不合格，洪峰合格率94.4%；垂向混合產(chǎn)流模型洪量相對(duì)誤差率定期2場(chǎng)、驗(yàn)證期1場(chǎng)，洪量合格率83.3%，洪峰相對(duì)誤差驗(yàn)證期有1場(chǎng)不合格，洪峰合格率94.4%。

從確定性系數(shù)角度，由表3可以看出：SCE-UA優(yōu)化后三水源新安江模型擬合洪水中確定性系數(shù)大于0.9的有3場(chǎng)，0.8～0.9之間的有5場(chǎng)，0.7～0.8之間的有5場(chǎng)，0.7以下的有5場(chǎng)；SCE-UA優(yōu)化后三水源新安江模型擬合洪水中確定性系數(shù)大于0.9的有2場(chǎng)，0.8～0.9之間的有9場(chǎng)，0.7～0.8之間的有4場(chǎng)，0.7以下的有3場(chǎng)?？梢?，優(yōu)化后的垂向混合產(chǎn)流模型模擬的洪水過程更符合實(shí)際。

圖3 典型洪水模擬過程線Fig.3 Simulation hydrograph for typical flood

5 結(jié) 論

本文將SCE-UA算法應(yīng)用到新安江模型和垂向混合產(chǎn)流模型參數(shù)優(yōu)化過程中，從田莊流域的模擬可以看出：①SCE-UA算法適于兩種模型的參數(shù)優(yōu)化，并能取得較好的模擬效果；②對(duì)于田莊水庫流域，優(yōu)化后的垂向混合產(chǎn)流模型模擬結(jié)果要略優(yōu)于三水源新安江模型；③SCE-UA優(yōu)化后的垂向混合產(chǎn)流模型比三水源新安江模型更適合模擬相關(guān)流域的洪水過程；④確定性系數(shù)不是很高，考慮可能因?yàn)樗x較流域處于北方半濕潤地區(qū)，流域面積小，洪水漲落迅速且流域洪水成因復(fù)雜；⑤兩種模擬洪水過程線退水段與實(shí)際有一定的差別，地下產(chǎn)流及匯流結(jié)構(gòu)有改善的空間。

□

[1] 孔凡哲，宋曉猛，占車生，等. 水文模型參數(shù)敏感性快速定量評(píng)估的RSMSobol方法[J]. 地理學(xué)報(bào)，2011，66(9):1 270-1 280.

[2] 郭 俊，周建中，周 超，等. 概念性流域水文模型參數(shù)多目標(biāo)優(yōu)化率定[J]. 水科學(xué)進(jìn)展，2012，23(4):447-456.

[3] 宋星原，舒全英，王海波，等. SCE-UA、遺傳算法和單純形優(yōu)化算法的應(yīng)用[J]. 武漢大學(xué)學(xué)報(bào)(工學(xué)版)，2009，42(1):6-9，15.

[4] Kuczera G. Efficient subspace probabilistic parameter optimization for catchment models[J]． Water Resources Research， 1997，33(1): 177-185．

[5] Sorooshian S， Duan Q Y， Gupta V K. Optimal use of the SCE-UA global optimization method for calibrating watershed models[J]. Journal of Hydrology， 1994，158(3-4):265-284.

[6] Eckhardt K， Arnold J G. Automatic calibration of a distributed catchment model [J]. Journal of Hydrology，2001，251(1-2):103-109.

[7] 董潔平，李致家，戴健男. 基于SCE-UA算法的新安江模型參數(shù)優(yōu)化及應(yīng)用[J]. 河海大學(xué)學(xué)報(bào)(自然科學(xué)版)，2012，40(5):485-490.

[8] 趙人俊. 流域水文模擬—新安江模型與陜北模型[M]. 北京:水利電力出版社，1984.

[9] 包為民，王從良. 垂向混合產(chǎn)流模型及應(yīng)用[J]. 水文，1997，(3):19-22.

[10] 瞿思敏，包為民，張 明，等. 新安江模型與垂向混合產(chǎn)流模型的比較[J]. 河海大學(xué)學(xué)報(bào)(自然科學(xué)版)，2003，31(4):374-377.

[11] 包為民. 格林—安普特下滲曲線的改進(jìn)和應(yīng)用[J]. 人民黃河，1993，(9):1-3.

[12] 舒 暢，劉蘇峽，莫興國，等. 新安江模型參數(shù)的不確定性分析[J]. 地理研究，2008，27(2):343-352.

[13] Duan Q Y， Gupta V K， S Sorooshian. Shuffled complex evolution approach for effective and efficient global minimization[J]. Journal of Optimization Theory and Applications， 1993， 76(3): 501-521.

[14] 李向陽，程春田，武新宇，等. 水文模型模糊多目標(biāo)SCE-UA參數(shù)優(yōu)選方法研究[J]. 中國工程科學(xué)，2007，9(3):52-57.