劉 春 燁

(太原理工大學(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)參過(guò)程中的不確定性，提高率定效率。水文模型敏感性分析方法有局部敏感性分析法、全局抽樣敏感性分析法。局部敏感性分析法原理簡(jiǎn)單但無(wú)法完全考慮參數(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ù)的率定效率。

1 LH-OAT全局敏感性分析方法

現(xiàn)階段處理參數(shù)敏感性的方法有隨機(jī)法、模糊法、區(qū)間分析法。前2種方法需要較多的數(shù)據(jù)，而區(qū)間估計(jì)法只需已知參數(shù)的上下界，方法較為簡(jiǎn)單。LH-OAT全局敏感性分析方法屬于區(qū)間估計(jì)法，是采用無(wú)量綱參數(shù)來(lái)表示全局敏感性的一種方法，該法的適用條件為已知模型參數(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為全局敏感度。

2 馬斯京根模型

現(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)。

馬斯京根法模擬河道流量過(guò)程中不可避免地存在著不確定性，對(duì)模型參數(shù)進(jìn)行全局敏感性分析是研究參數(shù)不確定性的分析方法之一。在傳統(tǒng)的馬斯京根模型參數(shù)分析過(guò)程中，是通過(guò)公式(4)、(5)、(6)來(lái)計(jì)算C0、C1、C2的值，本文提出直接對(duì)C0、C1、C2進(jìn)行敏感性分析，避免了由K、x取值不準(zhǔn)確對(duì)參數(shù)C0、C1、C2造成的影響。

3 實(shí)例分析

3.1 模型參數(shù)敏感性計(jì)算

為定性了解參數(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敏感性的無(wú)量綱數(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

3.2 參數(shù)敏感性結(jié)果分析

以徑流量作為敏感性評(píng)價(jià)目標(biāo)，從模型參數(shù)上來(lái)看，1～4時(shí)段(即河道漲水階段)馬斯京根模型參數(shù)針對(duì)不同流量會(huì)產(chǎn)生相應(yīng)波動(dòng)，但參數(shù)敏感性總趨勢(shì)為C0>C1>C2；5-8時(shí)段(即河道落水階段)參數(shù)敏感性趨勢(shì)為C0VC1>VC2。

圖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ù)的敏感性相差不大，無(wú)主要敏感性參數(shù)。分析敏感性變化的原因，是由于洪峰流量較河道普通流量大，當(dāng)流量過(guò)大時(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ù)的偏離程度、率定效率來(lái)評(píng)價(jià)2種方法的優(yōu)劣，結(jié)果見表4、圖3。方法1率定11次，精度為99.8%，方法2率定16次，精度為99.6%。因此，按敏感性由大到小進(jìn)行參數(shù)率定可提高69%的率定效率。以洪峰流量為率定目標(biāo)時(shí)，由于沒(méi)有主要敏感性參數(shù)，參數(shù)的敏感性差別不大，率定效率會(huì)有所下降。

表4 率定結(jié)果 %

圖3 率定結(jié)果比較Fig.3 Comparison of the calibration results

4 結(jié) 語(yǔ)

采用LH-OAT全局敏感性分析方法分析馬斯京根模型參數(shù)的敏感性，并通過(guò)實(shí)例進(jìn)行演算，在分析演算過(guò)程中得出以下結(jié)論。

(1)直接采用LH-OAT敏感性分析方法對(duì)C0、C1、C2進(jìn)行敏感性分析，確定了各參數(shù)對(duì)馬斯京根模型的影響程度，對(duì)河道資料較少、無(wú)法獲得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ù)敏感性差別較小，無(wú)主要敏感性參數(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.