劉 春 燁
(太原理工大學(xué)水利科學(xué)與工程學(xué)院,太原 030024)
馬斯京根法是河道洪水流量演算中廣泛應(yīng)用的一種水文學(xué)方法,該法在建立馬斯京根槽蓄方程的基礎(chǔ)上,與水量平衡方程聯(lián)立,推出上下游流量關(guān)系方程,進(jìn)行河段洪水的演算[1]。馬斯京根法應(yīng)用的一個(gè)前提是模型參數(shù)的合理率定,通常對(duì)參數(shù)K、x進(jìn)行率定,再由K、x率定值計(jì)算參數(shù)C0、C1、C2。參數(shù)K、x的率定方法有分析法和試算法,但這2種方法較為復(fù)雜且率定結(jié)果精確性較差。因此,部分學(xué)者采用了直接優(yōu)化參數(shù)C0、C1、C2的方法,主要有POS算法[1]、遺傳算法、最小二乘法、SCE-UA算法[2]等。
對(duì)模型參數(shù)進(jìn)行敏感性分析是參數(shù)率定的基礎(chǔ)。敏感性分析即分析參數(shù)對(duì)模型的影響程度,根據(jù)參數(shù)的影響程度,篩選出重要的影響參數(shù),減少參數(shù)的率定次數(shù)和調(diào)參過程中的不確定性,提高率定效率。水文模型敏感性分析方法有局部敏感性分析法、全局抽樣敏感性分析法。局部敏感性分析法原理簡(jiǎn)單但無法完全考慮參數(shù)間的相互影響,不能滿足提高模型精確性的要求。全局敏感性分析法可綜合分析多個(gè)參數(shù)的變化對(duì)模型運(yùn)行結(jié)果的影響,且該方法已在各種水文模型中得到了應(yīng)用,如RSA方法、GLUE方法[3]和LH-OAT方法[4]等。目前,國(guó)內(nèi)外直接對(duì)馬斯京根模型參數(shù)C0、C1、C2敏感性分析較少。本研究引入Latin-Hypercube抽樣算法,并結(jié)合OAT 方法,采用LH-OAT全局敏感性分析方法直接對(duì)馬斯京根模型參數(shù)C0、C1、C2進(jìn)行全局敏感性分析。在敏感性分析的基礎(chǔ)上,確定參數(shù)的敏感性大小順序,采用假定最優(yōu)參數(shù)法進(jìn)行馬斯京根模型參數(shù)率定,發(fā)現(xiàn)按敏感性順序進(jìn)行參數(shù)的率定,可提高馬斯京根模型參數(shù)的率定效率。
現(xiàn)階段處理參數(shù)敏感性的方法有隨機(jī)法、模糊法、區(qū)間分析法。前2種方法需要較多的數(shù)據(jù),而區(qū)間估計(jì)法只需已知參數(shù)的上下界,方法較為簡(jiǎn)單。LH-OAT全局敏感性分析方法屬于區(qū)間估計(jì)法,是采用無量綱參數(shù)來表示全局敏感性的一種方法,該法的適用條件為已知模型參數(shù)的取值范圍。LH-OAT全局敏感性分析方法的理論依據(jù)分別為OAT敏感性分析法和LH抽樣法。OAT敏感性分析法在Monte Carlo抽樣方法的基礎(chǔ)上進(jìn)行重復(fù)試驗(yàn),每次對(duì)其中一個(gè)參數(shù)進(jìn)行微小的擾動(dòng),保持其他參數(shù)不變,帶入敏感性公式,即可得到每個(gè)參數(shù)的敏感度,具有較高的精確性。但Monte Carlo抽樣法在參數(shù)取值范圍內(nèi)隨機(jī)抽樣,會(huì)產(chǎn)生大量的樣本,使OAT敏感性分析方法的運(yùn)算量變大。LH抽樣方法,假設(shè)抽樣符合均勻分布,將每個(gè)參數(shù)區(qū)間均分為n層,分別從每個(gè)分層中隨機(jī)抽樣一次,生成一個(gè)LH抽樣參數(shù)組,減少了采樣頻次,能提高敏感性分析的效率。LH-OAT全局敏感性分析法結(jié)合了LH抽樣法的高效性和OAT算法的精確性。
LH-OAT全局敏感性分析法具體步驟如下:①采用LH抽樣法,將整個(gè)參數(shù)空間均分為N層,分別從每個(gè)分層中隨機(jī)抽樣一次,組成一個(gè)抽樣參數(shù)組。②假設(shè)有P個(gè)模型參數(shù),按照OAT方法,對(duì)每個(gè)LH 抽樣參數(shù)組進(jìn)行P次參數(shù)的微小擾動(dòng)并且每次只擾動(dòng)其中一個(gè)參數(shù),計(jì)算每次擾動(dòng)前后目標(biāo)函數(shù)的變化值。③根據(jù)LH-OAT方法定義,模型總共需運(yùn)行N(P+1)次。每個(gè)參數(shù)的相對(duì)敏感度由以下公式計(jì)算得到:
(1)
式中:M為敏感性分析的目標(biāo)函數(shù);ei,k為第i個(gè)參數(shù)在第k層中抽樣值;Δei,k為ei,k的擾動(dòng)變化值;Si,k為參數(shù)ei,k在k層的相對(duì)敏感度。
由參數(shù)相對(duì)敏感度Si,k可計(jì)算參數(shù)的全局敏感度GSi,k:
(2)
式中:N為參數(shù)空間的總層數(shù);GSi,k為全局敏感度。
現(xiàn)階段河道洪水流量演算有水文學(xué)和水力學(xué)2種方法,在水文學(xué)中應(yīng)用最多的為馬斯京根法,其公式為:
Q2=C0I2+C1I1+C2Q1
(3)
式中:Q1為計(jì)算時(shí)段初始時(shí)刻的出流量;Q2為時(shí)段末的河段出流量;I1、I2分別為計(jì)算時(shí)段始、末的河段入流量;C0、C1、C2均為流量演算系數(shù),計(jì)算公式如下:
(4)
(5)
(6)
式中:K為蓄量常數(shù),即蓄量流量關(guān)系曲線的坡度;x為流量比重系數(shù);Δt為計(jì)算時(shí)段長(zhǎng)。
馬斯京根法模擬河道流量過程中不可避免地存在著不確定性,對(duì)模型參數(shù)進(jìn)行全局敏感性分析是研究參數(shù)不確定性的分析方法之一。在傳統(tǒng)的馬斯京根模型參數(shù)分析過程中,是通過公式(4)、(5)、(6)來計(jì)算C0、C1、C2的值,本文提出直接對(duì)C0、C1、C2進(jìn)行敏感性分析,避免了由K、x取值不準(zhǔn)確對(duì)參數(shù)C0、C1、C2造成的影響。
為定性了解參數(shù)C0、C1、C2對(duì)馬斯京根模型模擬洪水流量的影響程度,提出采用LH-OAT全局敏感性分析方法對(duì)模型參數(shù)進(jìn)行全局敏感性分析。以文獻(xiàn)[5]中數(shù)據(jù)為例,分析C0、C1、C2的敏感性,河段上斷面入流、 下斷面出流情況見表1,所選河道穩(wěn)定,計(jì)算時(shí)段長(zhǎng)為 Δt=6 h。其中時(shí)段1-4為連續(xù)時(shí)段,時(shí)段5-8為連續(xù)時(shí)段。由公式(4)、(5)、(6)利用多元函數(shù)求極值的方法得到C0、C1、C2參數(shù)的取值范圍[6]分別為(-9/11,3/13),(1/4,19/21),(1/21,1)。
表1 不同時(shí)段流量 m3/s
首先,將馬斯京根模型的3個(gè)參數(shù)C0、C1、C2在其取值范圍內(nèi)均分成10份,共11層;然后利用Latin-Hypercube抽樣法進(jìn)行取樣組合,可確定11種組合,且C0、C1、C2應(yīng)滿足公式:C0+C1+C2=1,以避免組合參數(shù)的不合理;按照OAT分析計(jì)算方法,對(duì)每種組合參數(shù)微小擾動(dòng)10%(保持2個(gè)參數(shù)不變,對(duì)另一個(gè)參數(shù)進(jìn)行微小擾動(dòng))后利用公式(1)進(jìn)行計(jì)算,每組數(shù)據(jù)需計(jì)算11×(3+1)=44次,4個(gè)連續(xù)時(shí)段可以確定3組計(jì)算數(shù)據(jù),有2組不同的連續(xù)時(shí)段,共需計(jì)算44×3×2=264次。利用公式(2)對(duì)相對(duì)敏感性計(jì)算結(jié)果進(jìn)行均值處理,最終得到6組反映3個(gè)參數(shù)C0、C1、C2敏感性的無量綱數(shù)值,具體計(jì)算結(jié)果見表2。從表1中隨機(jī)抽取一組數(shù)據(jù),分別把參數(shù)C0、C1、C2的擾動(dòng)幅度變?yōu)?0%、30%,每改變擾動(dòng)幅度一次需計(jì)算44次,變2次擾動(dòng)幅度,故需計(jì)算88次,分析變幅后的敏感性,計(jì)算結(jié)果見表3。
表2 參數(shù)的全局敏感度Tab.2 Global sensitivity of parameters
表3 不同擾動(dòng)幅度參數(shù)敏感度Tab.3 Sensitivity of different disturbance amplitude parameters
以徑流量作為敏感性評(píng)價(jià)目標(biāo),從模型參數(shù)上來看,1~4時(shí)段(即河道漲水階段)馬斯京根模型參數(shù)針對(duì)不同流量會(huì)產(chǎn)生相應(yīng)波動(dòng),但參數(shù)敏感性總趨勢(shì)為C0>C1>C2;5-8時(shí)段(即河道落水階段)參數(shù)敏感性趨勢(shì)為C0
圖1 不同時(shí)段參數(shù)敏感性比較Fig.1 Comparison of parameter sensitivity
圖2 不同擾動(dòng)幅度參數(shù)敏感度趨勢(shì)Fig.2 Sensitivity trend of different disturbance amplitude parameters
以洪峰流量作為敏感性評(píng)價(jià)目標(biāo),當(dāng)出現(xiàn)洪峰時(shí)(3~4時(shí)段),C0、C1、C2雖遵循漲水階段的規(guī)律,但C0、C1的敏感性相對(duì)前一時(shí)段的值減小,C2的敏感性相對(duì)增大,3個(gè)參數(shù)的敏感性相差不大,無主要敏感性參數(shù)。分析敏感性變化的原因,是由于洪峰流量較河道普通流量大,當(dāng)流量過大時(shí),河道的穩(wěn)定性減弱,影響因素相應(yīng)發(fā)生變化且影響因素增加。
采用假定最優(yōu)參數(shù)法分析參數(shù)敏感性對(duì)率定效率和率定結(jié)果的影響,首先選定一組參數(shù)作為最優(yōu)參數(shù),分別采用2種方法進(jìn)行參數(shù)率定,第1種方法是根據(jù)參數(shù)的敏感性大小順序進(jìn)行參數(shù)率定,第2種方法是根據(jù)擾亂參數(shù)敏感性順序進(jìn)行率定。以率定結(jié)果與最優(yōu)參數(shù)的偏離程度、率定效率來評(píng)價(jià)2種方法的優(yōu)劣,結(jié)果見表4、圖3。方法1率定11次,精度為99.8%,方法2率定16次,精度為99.6%。因此,按敏感性由大到小進(jìn)行參數(shù)率定可提高69%的率定效率。以洪峰流量為率定目標(biāo)時(shí),由于沒有主要敏感性參數(shù),參數(shù)的敏感性差別不大,率定效率會(huì)有所下降。
表4 率定結(jié)果 %
圖3 率定結(jié)果比較Fig.3 Comparison of the calibration results
采用LH-OAT全局敏感性分析方法分析馬斯京根模型參數(shù)的敏感性,并通過實(shí)例進(jìn)行演算,在分析演算過程中得出以下結(jié)論。
(1)直接采用LH-OAT敏感性分析方法對(duì)C0、C1、C2進(jìn)行敏感性分析,確定了各參數(shù)對(duì)馬斯京根模型的影響程度,對(duì)河道資料較少、無法獲得K、x值的地區(qū)的河道演算具有較好的適用性。
(2)以徑流量為評(píng)價(jià)目標(biāo),當(dāng)河道漲水時(shí),參數(shù)C0較敏感,C1、C2一般敏感;當(dāng)河道落水時(shí),C0、C1一般敏感,C2較為敏感,因此,建議在漲水和落水階段分別分析參數(shù)的敏感性。以洪峰流量作為評(píng)價(jià)目標(biāo)時(shí),C0、C1的敏感性相對(duì)減弱,C2的敏感性相對(duì)增加,但3個(gè)參數(shù)敏感性差別較小,無主要敏感性參數(shù)。
(3)按敏感性大小順序進(jìn)行參數(shù)率定,可以大幅度提高馬斯京根模型參數(shù)的率定效率。在實(shí)際工程計(jì)算時(shí)需針對(duì)研究河道實(shí)際情況研究參數(shù)C0、C1、C2的敏感性,以期最接近河道流量的真實(shí)情況,敏感性分析成果為馬斯京根模型參數(shù)的合理率定奠定了基礎(chǔ)。
□
[1] 李 匡,胡宇豐,梁犁麗,等.PSO 算法在馬斯京根法參數(shù)率定中的應(yīng)用[J].水電站機(jī)電技術(shù),2016,39(8):88-91.
[2] 李鴻雁,趙 娟,王玉新,等.擴(kuò)域搜索遺傳算法優(yōu)化馬斯京根參數(shù)及其應(yīng)用[J].吉林大學(xué)學(xué)報(bào),2011,41(3):861-865.
[3] 李 丹,張 翔,張 揚(yáng).水文模型參數(shù)敏感性的區(qū)間分析[J].水利水電科技進(jìn)展,2011,31(1):29-32,41.
[4] 劉慧如,吳建華,孫 毅,等.基于LH-OAT的雙曲正切模型參數(shù)敏感性分析[J].水電能源科學(xué),2016,34(9):84-86.
[5] 林三益.水文預(yù)報(bào)[M]. 北京:中國(guó)水利水電出版社,2001.
[6] 李 匡,付 力,胡宇豐,等.馬斯京根法參數(shù)C0、C1、C2取值范圍的確定[J].南水北調(diào)與水利科技,2012,10(5):43-45,55.