大伙房流域降雨徑流模型

1、.大伙房流域降雨徑流模型大伙房流域降雨徑流預(yù)報(bào)模型又簡(jiǎn)稱“dhf”模型，該模型于1973 年由遼寧省大伙房水庫(kù)管理局劉愛(ài)杰、王本德等人提出，至今已使用30 余年，為水庫(kù)洪水調(diào)度做出了很大貢獻(xiàn)?！?dhf”模型是適用于我國(guó)濕潤(rùn)地區(qū)的超滲產(chǎn)流模型，目前已在遼寧省多個(gè)水庫(kù)的水情自動(dòng)測(cè)報(bào)系統(tǒng)中使用，效果較好。建立在“ dhf”模型基礎(chǔ)上的降雨徑流預(yù)報(bào)方案，在大伙房流域經(jīng)過(guò)調(diào)試和精度驗(yàn)證后進(jìn)行使用，在使用中平均精度令人滿意。尤其在“ 957”特大洪水調(diào)度中，發(fā)揮了顯著作用，準(zhǔn)確預(yù)報(bào)出了第一非常溢洪道溢流時(shí)間和水庫(kù)最高庫(kù)水位，為省防制定調(diào)度決策提供了科學(xué)依據(jù)，使水庫(kù)工程發(fā)揮了強(qiáng)大的調(diào)蓄作用，最大限度地配合

2、了下游搶險(xiǎn)，共減免下游直接經(jīng)濟(jì)損失74.89 億元?！?dhf”模型由兩部分組成，一是八參數(shù)超滲產(chǎn)流計(jì)算模型，引用雙層入滲曲線進(jìn)行扣損計(jì)算，并以拋物線描述表層蓄水量和下層滲率的分布狀況；二是八參數(shù)變強(qiáng)度、變速度的經(jīng)驗(yàn)單位線匯流計(jì)算模型，參數(shù)隨降雨分布而變，采用“前期影響凈雨”描述匯流速度的變化。這是一個(gè)集總的概念模型，模型的參數(shù)多半在滿足其物理意義的前提下確定，只有6 個(gè)需要優(yōu)選法選定或試錯(cuò)法確定。1 大伙房模型概化流程流域下墊面分為表層、下層和地下水蓄存三部分，計(jì)算流程如圖10-1 所示。2 大伙房模型產(chǎn)流計(jì)算產(chǎn)流模型將下墊面分為表層，下層和深層三部分。表層土壤中的張力水蓄量與植物截流、填洼

3、儲(chǔ)存合稱表層蓄水量sa ，其極值為表層蓄水容量 s0 ；下層土壤中的張力水蓄量稱為下層蓄水量u a ，其極值為下層蓄水容量 u 0 ；地下水儲(chǔ)水層的蓄水量以va 表示，其極值為地下水庫(kù)蓄水容量v0 。;.ed=kc*er/c雨強(qiáng) p蒸散發(fā)雨間蒸發(fā)pe =p-ed不透水 面積上pc=pe -y 0直 接 徑流蒸 散 發(fā)sa/s0y =g*p表 層 存儲(chǔ)+sa地 面 壤 中r=pc-s0rc=r-yu流蒸 散 發(fā)ua/u0 下層存儲(chǔ)rl =rc(u0-ua)蒸 散 發(fā)va/v0地下水存儲(chǔ)地下徑流y l=kw*r l深層地下水存儲(chǔ)總出流圖 10-1大伙房產(chǎn)流模型概化流程圖其中：s0- 表層蓄水容量；

4、u0- 下層蓄水容量；d0- 下層蓄水容量（ u0）與地下水庫(kù)蓄水容量（v0）之和；g- 不透水面積占全流域面積之比值；kc- 流域蒸散發(fā)能力與大水體蒸發(fā)量的比值；a- 表層蓄水容量曲線形狀參數(shù)；b- 下滲率拋物線分配曲線形狀參數(shù)；k- 下層下滲曲線曲率；2.1不透水面積上的徑流計(jì)算根據(jù)上面的流程圖，不透水面積上直接產(chǎn)生徑流，通過(guò)下式來(lái)計(jì)算：ycgpepepedpkc er / c（1）（2）式中 p 為時(shí)段降雨量，kc 是流域蒸散發(fā)能力與er 大水體蒸發(fā)量的比值，c 為日降雨觀測(cè)的時(shí)段數(shù)，g 為流域不透水面積占總流域面積的比值。;.2.2地面壤中流的計(jì)算地面壤中流通過(guò)下式計(jì)算yu rf（3）

5、2.2.1 表層扣損 r的計(jì)算當(dāng)smpcas0時(shí)，r pcsa s0s0 (1 (smpc ) /(as0 ) a（4）當(dāng) smpcas0 時(shí)， rpc (s0sa )（ 5）式中 sa 是表層流域平均蓄水量，sm是相應(yīng)于 sa 的點(diǎn)蓄水量； a 為曲線形狀參數(shù)。sm和 sa 由下式計(jì)算得到：smas0 1 (1 sa / s0 ) 1/ a（6）sasm / a 1 (1 sm / sm )（7）式中 sm 為表層點(diǎn)最大蓄水量，當(dāng)smsm 時(shí)，則 sa s0 ，并有smas0（8）其中 pc 通過(guò)下式來(lái)計(jì)算：pcp edyc（9）2.2.2下滲計(jì)算時(shí)段平均總下滲率f與da的關(guān)系 :fz1(

6、d0da )(10)時(shí)段平均下層下滲率fu與da的關(guān)系 :f uz2(u0ua )(11)其中 z1 1 e k2 vtu0 / d0 ; z21e k2 vt， k1 和 k2 的關(guān)系如下 :k1k2u0/ d0(12)k1 和 k2 分別是反映總下滲和下層下滲曲線曲率系數(shù)；da 是 ua 的函數(shù)，二者的關(guān)系式如下 :dad01 (1 ua / u0 )u0 / d0(13)d0 是下層蓄水容量u0 與地下水庫(kù)蓄水容量v0 之和。2.2.3徑流深計(jì)算公式;.總徑流 y的計(jì)算公式：z2un rby r z2 (uau0 ) z2u01z2 bu0(14)當(dāng) rz2unz2 bu0 時(shí)，得y r

7、 z2 (u0ua ) r fu(15)地面壤中流yu 的計(jì)算公式：當(dāng) z1dnrz1bd0 時(shí)，yur z1 d0 (1 ua / u0 )u0 / d0z1d1 (1z1 dnr ) bz1bd0(16)當(dāng)z1dnrz1bd0 時(shí)，yur z1 d0(1 ua / u0 )u0 / d0(17)d n 和 un 通過(guò)下列兩式計(jì)算dnbd0 1 (1 ua / u0 )u0 /(bd 0 )(18)unbu01 (1 ua / u0 )1/ b(19)地下徑流 yl 的計(jì)算：ylkw rl(20)其中 kw 為 yl 與 rl 之比值，在干旱地區(qū)rl 全部滲入深層，不形成地下徑流，則kw0

8、；在濕潤(rùn)地區(qū) rl 全部形成地下徑流流入河道，所以kw1。 rl 為地下水庫(kù)下滲強(qiáng)度，為總徑流與壤中流之差，即：rly yu(21)2.3大伙房模型蒸散發(fā)計(jì)算（ 1）當(dāng)降雨量超過(guò)流域蒸發(fā)能力時(shí)，即凈滲雨強(qiáng)pc0 ，產(chǎn)生下滲，土壤蓄水量增加，此時(shí)：表層蓄水量計(jì)算公式sa(t 1)s0 1 (1 (smt pct ) / as0 )a(22)式中 t 為單位日，以8:00 為日分界；sm 以式（ 6）計(jì)算當(dāng) smtpct as0 時(shí)，上式可寫成;.sa (t 1)s0(23)依照流程圖的下滲程序，逐時(shí)段計(jì)算式，sa (i 1)saipciri(24)式中 i 為時(shí)段序；pc 、 r 分別用（ 9）

9、、（ 4）、（ 5）式進(jìn)行計(jì)算。下層土壤蓄水計(jì)算公式據(jù)流程圖，逐時(shí)段計(jì)算ua 的公式ua (i 1) uai riyi(25)若逐日計(jì)算，可用下式計(jì)算若 z2 untrtz2 bu0 時(shí)ua(t 1)uatz2(u0uat )z2u0 (1 z2untrt )bz2 bu0(26)若 z2untrtz2 bu0 時(shí)ua (t1)uatz2 (u0 uat )(27)式中的un與z2分別用式（19）及z21e k2 vt求出。（ 2）當(dāng)日流域平均降雨量小于流域蒸發(fā)能力，即ped 時(shí)，土壤蓄水量將減少進(jìn)行蒸散發(fā)計(jì)算表層蒸散發(fā)公式假定： a. 表層蒸發(fā)等于流域平均蒸散發(fā)能力；b. 每次下滲后的初次蒸

10、散發(fā)皆從地表開始。eu 為表層可蒸發(fā)量，計(jì)算公式為：es (1eb ,t 1 )a(1eb ,t ) au ,t0as0as0(28)式中：ec, teu, tptt1eb ,t 1ecneaeb,teb ,t 1ec,t(29)因此， t1日表層蓄水量計(jì)算式為sa ,t 1sa ,teu, t(30);.下層蓄水量 ua 計(jì)算公式假定： a. 表層蓄水量蒸發(fā)盡后，該點(diǎn)的下層開始蒸發(fā)；b. 下層蒸發(fā)量 el 與下層蓄水量 ua 成比例關(guān)系。因此，el(eceu )ua/ u0(31)當(dāng) ( s eu )0 時(shí)，則el( ecsa )ua / u0(32)這樣， t1日下層土壤蓄水量以下式計(jì)算u

11、a,t 1ua,tel(33)2.4大伙房模型匯流計(jì)算大伙房匯流模型是根據(jù)多年預(yù)報(bào)經(jīng)驗(yàn)，總結(jié)出的8參數(shù)變強(qiáng)度、變匯流速度的經(jīng)驗(yàn)單位線，其線型公式為cctk 3 e aa( ?t / n?tm )ddsintntmntm式中t匯流曲線，是tm 與 t 的函數(shù)；tm匯流曲線底寬；t 匯流時(shí)間（時(shí)刻）；k3 匯流曲線的比例系數(shù)；n 地下徑流與地面壤中流匯流曲線底寬的比例系數(shù)；aa、 dd、 cc反映匯流曲線形狀的參數(shù)；地表徑流與地下徑流具有不同的匯流曲線形狀，可用形狀參數(shù)區(qū)分；不同的降雨分布地表徑流也具有不同的匯流曲線，其參數(shù)也不同。tm 是時(shí)段凈雨（徑流）與前期影響凈雨（已產(chǎn)生的徑流）的函數(shù)，模型

12、中假定：tmlb ( yayt )k 0ya, t 1k ( ya ,tyt ,t )式中l(wèi)b 特征河長(zhǎng)lbl / b0 ；l 流域最大河長(zhǎng)；;.k 0 擬合線的斜率；ya 前期影響凈雨量，表征河槽蓄水狀況，制約匯流速度等；k 前期凈雨影響程度衰減系數(shù)，反映已產(chǎn)生的徑流在河槽中消退速度；2.5大伙房模型產(chǎn)流參數(shù)優(yōu)選大伙房模型產(chǎn)流參數(shù)的最優(yōu)化問(wèn)題是多目標(biāo)的有約束的非線性規(guī)劃問(wèn)題。目標(biāo)函數(shù)取為n次洪水產(chǎn)流的擬合精度和n次洪水計(jì)算產(chǎn)流與實(shí)際產(chǎn)流誤差總和的線性加權(quán)，各取0.5 的權(quán)重進(jìn)行求解。將問(wèn)題變?yōu)椋簃ax w1* acu 1 w 2 * tytc s.t. xi minxixi max式中acu

13、1-n次洪水產(chǎn)流的擬和精度tytc-n次洪水計(jì)算產(chǎn)流與實(shí)際產(chǎn)流誤差總和acu 1n合格 / ntytc(totaltytc) / total其中 n為洪水的次數(shù)，total為 n次洪水的實(shí)際凈雨的總和，tytc為 n次洪水的預(yù)報(bào)凈雨的總和，n 合格為擬合的合格洪水的次數(shù)。若某次洪水滿足以下條件之一即為合格：計(jì)算產(chǎn)流與實(shí)際產(chǎn)流的絕對(duì)誤差小于3mm；計(jì)算產(chǎn)流與實(shí)際產(chǎn)流的相對(duì)誤差小于20%且絕對(duì)誤差小于20mm；計(jì)算時(shí)判斷合格的絕對(duì)誤差和相對(duì)誤差的標(biāo)準(zhǔn)可以根據(jù)具體情況來(lái)設(shè)定；本程序采用可變?nèi)莶罘ㄇ蠼鈫?wèn)題，把多個(gè)約束求極小值問(wèn)題變?yōu)橐粋€(gè)單約束求極小值問(wèn)題。設(shè)所求解的非線性規(guī)劃問(wèn)題為：minf ( x

14、)xenhi ( x )0(i1 , 2,.,m )s .t .0(jm1 ,., p )gi ( x )本算法將上述兩問(wèn)題化為單約束問(wèn)題進(jìn)行求解:minf ( x )xe ns . t .( k)t( x )0式中 t(x) 是約束破壞的估計(jì)量，表達(dá)式為：mp1/ 2t (x)hi 2 ( x)u j g j 2 (x)i 1j m 1其中， uj 為 heaviside 算子。( k)是第 k步搜索中給出的關(guān)于可行性的可變?nèi)莶顪?zhǔn)則值，表達(dá)式為：;.(k )min(k 1) , m1x(k )x(k )r1ir 2(0)2(m 1)t其中t-初始多面體的大小 ;m-等式約束數(shù)目 ;r=n-m

15、 ， f ( x) 的自由度 ;xi ( k )n- e 多面體中第 k次搜索時(shí)的第 i 個(gè)頂點(diǎn)值 ;xr 2(k )xh-多面體（除）外的形心頂點(diǎn) ;由于各參數(shù)的初值相差較大，而且變化范圍也各不相同，所以將參數(shù)歸一化，即初值都取為1.0 。根據(jù)物理意義將e6、 e7、 e8、 e9的變化范圍取為 0.5-10.0，其它參數(shù)取為 0.5-5.0 。2.6大伙房模型匯流參數(shù)優(yōu)選大伙房匯流模型需要優(yōu)選的參數(shù)有b0、 k0、 k、n、 dd、cc，其中 dd和cc根據(jù)降雨的分布分為上游、中游和下游三套參數(shù)，所以需要優(yōu)選的參數(shù)一共有10個(gè)。 k3和 aa可以通過(guò)優(yōu)選出的參數(shù)帶入下列兩式計(jì)算求得：ntm

16、cct )k 3( ntm0 ) /e aa( t / ntm )ddsin(0ntmaacc / dd ( t0 ) /(ntm ) dd1 tan ( t0 ) / ntm大伙房匯流參數(shù)優(yōu)選的基礎(chǔ)是多次歷史洪水，其最優(yōu)化問(wèn)題是多目標(biāo)的有約束的非線性規(guī)劃問(wèn)題。目標(biāo)函數(shù)取為 n次洪水過(guò)程的擬合精度、n次洪水計(jì)算洪峰與實(shí)際洪峰誤差總和以及n次洪水的計(jì)算峰現(xiàn)時(shí)間和實(shí)際的峰現(xiàn)時(shí)間誤差總和的線性加權(quán)，各取0.4 、 0.3 、 0.3 的權(quán)重進(jìn)行求解。將問(wèn)題變?yōu)椋簃inw1qe 1w2 qe 2w3 qe 3s.t .xi minxixi max式中qe1-n次洪水匯流的擬和精度;qe2-n次洪水計(jì)算

17、峰值與實(shí)際峰值誤差總和;qe3-n次洪水計(jì)算峰現(xiàn)時(shí)間與實(shí)際峰現(xiàn)時(shí)間誤差總和;本程序采用可變?nèi)莶罘ㄇ蠼鈫?wèn)題，把多個(gè)約束求極小值問(wèn)題變?yōu)橐粋€(gè)單約束求極小值問(wèn)題. 設(shè)所求解的非線性規(guī)劃問(wèn)題為:minf ( x )xnehi ( x )0(i1, 2,., m )s .t .0(jm 1,., p )gi ( x )本算法將上述兩問(wèn)題化為下式進(jìn)行求解:;.minf ( x ) xe ns.t .( k )t ( x )0式中 t(x) 是約束破壞的估計(jì)量，表達(dá)式為:m1/ 2pt (x)hi 2 ( x)u j g j 2 (x)i 1j m 1其中， uj 為 heaviside 算子。( k)是第

18、 k步搜索中給出的關(guān)于可行性的可變?nèi)莶顪?zhǔn)則值，表達(dá)式為：( k )min( k1) , m1xi ( k ) xr 2 ( k )(0)r12( m1)t其中t-初始多面體的大小;m-等式約束數(shù)目 ;r=n-m ， f ( x) 的自由度 ;xi ( k ) -en 多面體中第 k次搜索時(shí)的第 i 個(gè)頂點(diǎn)值 ;(k )xr 2- 多面體（除xh ）外的形心頂點(diǎn) ;2.7大伙房產(chǎn)流模型參數(shù)為提高產(chǎn)流預(yù)報(bào)精度，按降雨開始時(shí)的起始流量大小進(jìn)行分類。分類預(yù)報(bào)參數(shù)如表10-1 表 10-2 所示。表 10-1降雨起始流量100m3/s 的大伙房產(chǎn)流模型參數(shù)gabks0u0d00.054.1291.1810.2845100130續(xù)表 10-1流域蒸發(fā)能力 ed月份6789日雨量 5mm3.74.34.33.75mm 日雨量 1.91.91.91.930mm日雨量 30mm1.01.01.01