水文預報課程設計

1、水文預報課程設計報告 學 院：_水利與環(huán)境學院_專 業(yè)：_水文與水資源工程_班 級： 200905201 姓 名：_馬天玉_學 號：_20090520115_ 指導教師：_胡彩虹_ 第一章 基本任務1.1蒸發(fā)折算系數(shù)Kc的優(yōu)選根據(jù)已給數(shù)據(jù)資料及參數(shù)（本報告采用89-92年的歷史數(shù)據(jù)），將流域作為整體：（1）進行日模型產(chǎn)流量計算；（2）比較計算年徑流與實測年徑流；（3）通過誤差分析，優(yōu)選蒸發(fā)折算系數(shù)Kc；（4）8990年的歷時數(shù)據(jù)作為率定參數(shù)，9192年的數(shù)據(jù)作為模型檢驗。1.2暴雨預報根據(jù)已給的設計暴雨資料和任務一率定的Kc，將流域作為整體進行如下計算：（1）次洪產(chǎn)流量計算，劃分水源；（2）直

2、接徑流匯流，地下徑流匯流的計算。（3）采用2004年暴雨數(shù)據(jù)進行預報。 根據(jù)已給的資料、參數(shù)及做過的習題，自己編寫程序，將流域作為整體進行產(chǎn)流量計算；將計算年徑流與實測年徑流進行比較。第二章 基本資料2.1流域概況白盆珠水庫位于廣東省東江一級支流西枝江的上游，壩址以上集雨面積856 km2。流域地處粵東沿海的西部，海洋性氣候顯著，氣候溫和，雨量豐沛。暴雨成因主要是鋒面雨和臺風雨，常受熱帶風暴影響。降雨年際間變化大，年內(nèi)分配不均，多年平均降雨量為1800mm，實測年最大降雨量為3417mm，汛期49月降雨量占年降雨量的81左右：徑流系數(shù)0.50.7。流域內(nèi)地勢平緩，土壤主要有黃壤和砂壤，具有明顯

3、的腐殖層，淀積層和母質(zhì)土等層次結(jié)構(gòu)，透水性好。臺地、丘陵多生長松、杉、樟等高大喬木；平原則以種植農(nóng)作物和經(jīng)濟作物為主，植被良好。流域上游有寶口水文站，流域面積553km2，占白盆珠水庫壩址以上集雨面積的64.6%。白盆珠水庫有10年逐日入庫流量資料、逐日蒸發(fā)資料和時段入庫流量資料：流域內(nèi)有7個雨量站，其中寶口以上有4個。雨量站分布較均勻有10年逐日降水資料和時段降水資料；寶口水文站具有10年以上水位、流量資料；流域?qū)偕絽^(qū)性小流域且受到地形、地貌等下墊面條件影響，洪水陡漲緩落，匯流時間一般23h，有時更短；一次洪水總歷時25d。圖2-1 該地區(qū)水文站分布圖2.2數(shù)據(jù)資料(1) 計算流域面積為55

4、3km2。(2) 流域內(nèi)有四個雨量站，權(quán)重系數(shù)分別為0.33、0.14、0.33、0.20。(3) 日產(chǎn)流模型數(shù)據(jù)，具體見資料文件名：87-92data.xls，數(shù)據(jù)格式為：T(i) Q(i) E(i) P1(i) P2(i) P3(i) P4(i)(4)暴雨預報的數(shù)據(jù)，見表2-1表2-1 2004年暴雨過程數(shù)據(jù)表時間蒸散發(fā)(mm)降雨量（mm）禾多布馬山高潭寶口T（i）E（i）P1（i）P2（i）P3（i）P4（i）2004-9-23 12:001.36.29.921.617.32004-9-23 15:001.37.61620.612.62004-9-23 18:001.36.26.414

5、.915.92004-9-23 21:001.38.817.229.418.52004-9-24 0:001.22534.835.324.62004-9-24 3:000.929.929.243.937.82004-9-24 6:000.938.624.846.9332004-9-24 9:000.96.97.56.112.32004-9-24 12:000.928.329.934.228.52004-9-24 15:000.925.642.739.875.42004-9-24 18:000.993.9137.612413.22004-9-24 21:000.985.390.88575.920

6、04-9-25 0:000.851.547.749.238.52004-9-25 3:001.139.870.342.197.72004-9-25 6:001.143.247.361.545.92004-9-25 9:001.120.513.315.813.12004-9-25 12:001.110.581.83.32004-9-25 15:001.17.48.47.610.92004-9-25 18:001.11.82.82.14.62004-9-25 21:001.10.200.302004-9-26 0:001.200002004-9-26 3:002.100002004-9-26 6:

7、002.100002004-9-26 9:002.100002004-9-26 12:002.100002004-9-26 15:002.100002004-9-26 18:002.100002004-9-26 21:0020000（5）計算參數(shù)數(shù)據(jù)，見表2-2 表2-2 計算參數(shù)表計算年份參數(shù)初始張力水蓄量19891990WmUmLmDmWWUWLWD140206060110104060BCFcIM0.20.16220.001（6）流域單位線單位線過程（m3/s）為：0，40，80，130，100，80，48，20，10，5，0（7）地下徑流匯流Cg=0.978，Qg=55.3m3/s第三章

8、 計算公式 該流域海洋性氣候顯著、氣候溫和、雨量豐沛，多年平均降雨量為1800mm，徑流系數(shù)0.5-0.7，土壤主要有黃壤和砂壤，層次結(jié)構(gòu)明顯，透水性好，植被覆蓋度高，地勢平坦，由此可初步判定該流域的產(chǎn)流機制為蓄滿產(chǎn)流模式。3.1產(chǎn)流計算 蒸散發(fā)計算根據(jù)流域蓄滿產(chǎn)流特點，蒸散發(fā)計算采用的是三層蒸散發(fā)計算模式。三層蒸發(fā)模式的具體計算如下：1)當WU+PEP, EU=Ep,EL=0,ED=0;2)當WU+P<EP, WLC×WLM, EU=WU+P,EL=(EP-EU)×WL/WLM,ED=0;3)當WU+P<EP, C(EP-EU) WL<C×WL

9、M, EU=WU+P,EL=C(EP-EU),ED=0;4)當WU+P<EP，WL<C(EP-EU)， EU=WU+P,EL=WL,ED=C(EP-EU)-EL.RW 圖2-1 包氣帶蓄水容量曲線式中：WU為上層土壤蓄水量，WL為下層土壤蓄水量，EU為上層土壤蒸發(fā)量，EL為下層土壤蒸發(fā)量，ED為深層土壤蒸發(fā)量，P為流域平均降雨量，Ep為流域平均蒸發(fā)能力，C為深層蒸散發(fā)擴散系數(shù)，WLM為下層張力水蓄水容量。產(chǎn)流量計算 根據(jù)流域特點，產(chǎn)流量計算系根據(jù)蓄滿產(chǎn)流理論得出的。蓄滿產(chǎn)流，即任一地點上，土壤含水量達田間持水量前，降雨量全部補充土壤含水量，不產(chǎn)流；當土壤蓄滿后，其后續(xù)降雨量全部產(chǎn)生

10、徑流。流域內(nèi)各點包氣帶的蓄水容量是不同的，將各點包氣帶蓄水容量從小到大排列，以包氣帶達到田間持水量時的土壤含水量WM為縱坐標，以流域內(nèi)小于等于該WM的面積占全流域的面積比為橫坐標，所繪的曲線稱為流域蓄水容量曲線。a=WMM× (1-(1-W/WM) 1/(b+1)PE>0，則產(chǎn)流；否則不產(chǎn)流。產(chǎn)流時： 1)當PE+aWMM： R=PE+W-WM+WM×(1-(PE+a)/WMM)b+12)當PE+a>WMM： R=PE+W-WM式中：PE為扣除蒸發(fā)量后的降雨量，a為土壤含水量W對應的土壤水深，WM為流域平均蓄水容量，WMM為流域各地點包氣帶蓄水容量的最大值，b為

11、流域包氣帶蓄水容量分布的不均勻指數(shù)，R為流域產(chǎn)流量。3.1.3二水源劃分流域坡地上的降雨產(chǎn)流量因產(chǎn)流過程的條件和運動路徑不同，受流域的調(diào)蓄作用不同，各徑流成分在流量過程線上的反應是不一樣的。在實際工作中，常需按各種徑流成分分別計算或模擬，因為要對產(chǎn)流量進行水源劃分。直接徑流和地下徑流水源劃分如下：1）當PE<=FC時：RS=0.0 RG=R2）當PEFC時：RG=FC*R/PE RS=R-RG式中：FC為穩(wěn)定下滲率，RS為直接徑流，RG為地下徑流。3.1.4各層蓄水量計算降雨補充土壤含水量，由前一天的土壤含水量推求第二天的土壤含水量，補充來源為降雨減去蒸散發(fā)減去徑流量，順序為上、下、深層

12、依此補充。三層蓄水量變化的具體計算如下：1）WUi+Pi-EUi-Ri<=UM，WUi+1=WUi+Pi-EUi-Ri;WLi+1=WLi-ELi;WDi+1=WDi-EDi;2）WUi+Pi-EUi-Ri>UM，WLi-ELi+(WUi+Pi-EUi-Ri-UM)<=LM， WUi+1=UM;WLi+1=WLi-ELi+(WUi+Pi-EUi-Ri-UM);WDi+1=WDi-EDi;3)WUi+Pi-EUi-Ri>UM，WLi-ELi+(WUi+Pi-EUi-Ri-UM)>LM， WDi-EDi+WLi-ELi+(WUi+Pi-EUi-Ri-UM)-LM<

13、;=DM時， WUi+1=UM; WLi+1=LM; WDi+1=WDi-EDi+WLi-ELi+(WUi+Pi-EUi-Ri-UM)-LM; 4)WUi+Pi-EUi-Ri>UM，WLi-ELi+(WUi+Pi-EUi-Ri-UM)>LM， WDi-EDi+WLi-ELi+(WUi+Pi-EUi-Ri-UM)-LM>DM時， WUi+1=UM; WLi+1=LM; WDi+1=DM; Wi+1=WUi+1+WLi+1+WDi+1；式中：i表示第i天。3.2匯流計算根據(jù)流域凈雨和流域徑流單位線，采用卷積的差分形式算出流域出口的流量過程。直接徑流匯流可根據(jù)該流域的時段單位線推求

14、，地下徑流匯流由線性水庫演算法推求。1）直接徑流匯流計算公式：QS(i)=RS(i)×UH；式中：UH為該流域的單位線2）地下徑流匯流計算公式：QG(i)=CG×QG(i-1)+(1-CG)×RG(i)×UQ(i)=QS(i)+QG(i)式中：U為單位轉(zhuǎn)換系數(shù)，3.3模型參數(shù)1）Kc：蒸散發(fā)能力折算系數(shù)，它主要反映流域平均高程與蒸發(fā)站高程之間差別的影響和蒸發(fā)皿散發(fā)與路面蒸散發(fā)間差別的影響；2）WM:流域平均張力水容量，它表示流域蓄滿的標準；3）WUM:上層張力水蓄水容量，它包括了植物截留量；4）WLM：下層張力水蓄水容量；5）b：流域包氣帶蓄水容量分布的

15、不均勻指數(shù)，在一般情況下其取值與單元流域面積有關；6）C:深層蒸散發(fā)擴散系數(shù)，它主要取決于流域內(nèi)深根植物的覆蓋范圍。7）IM:不透水面積占全流域面積的比例，它的值可由大比例尺的地形圖，通過地理信息系統(tǒng)現(xiàn)代技術(shù)量測出來，也可用歷史上干旱期小洪水資料來分析。第四章 模型檢驗、結(jié)果評定及分析水文預報是一項直接服務于國家安全和國民經(jīng)濟建設的不可或缺的重要基本工作，是幫助人類有效地預防洪水、減少洪災損失，有效利用水資源的非工程措施之一。隨著經(jīng)濟、社會發(fā)展及其全球化進程的需要，水文預報的服務面進一步拓展，對水文預報提出了更高的要求。水文預報結(jié)果的準確率與可信程度是衡量服務質(zhì)量的前提，為了更好地為國家安全和

16、國民經(jīng)濟建設服務，必須對水文預報結(jié)果的可靠性和有效性進行評定和檢驗。4.1產(chǎn)流模式的檢驗 定性分析該流域集雨面積856km2。流域地處南方，海洋性氣候顯著，氣候溫和，雨量豐沛。暴雨成因主要是鋒面雨和臺風雨，常受熱帶風暴影響。降雨年際間變化大，年內(nèi)分配不均，多年平均降雨量為1800mm，實測年最大降雨量為3417mm，汛期49月降雨量占年降雨量的81左右：徑流系數(shù)0.50.7。流域內(nèi)地勢平緩，土壤主要有黃壤和砂壤，具有明顯的腐殖層，淀積層和母質(zhì)土等層次結(jié)構(gòu)，透水性好。臺地、丘陵多生長松、杉、樟等高大喬木；平原則以種植農(nóng)作物和經(jīng)濟作物為主，植被良好。流域上游有一水文站，控制流域面積553km2，占

17、流域集雨面積的64.6%。該水文站以上有4個雨量站。雨量站分布較均勻，有10年逐日降水資料和時段降水資料；該水文站具有10年以上水位、流量資料；流域?qū)偕絽^(qū)性小流域且受到地形、地貌等下墊面條件影響，洪水陡漲緩落，匯流時間一般23h，有時更短；一次洪水總歷時25d。由流域概況可以看出，該地區(qū)屬于濕潤多雨地區(qū)，雨量集中，地勢平坦，土壤層容易蓄滿，而且有場系列的降雨資料和水位流量資料，綜合從氣候條件、土壤狀況、植被組成以及洪水機制看，新安江模型適用于該地區(qū)的水文預報。4.2率定期優(yōu)選蒸發(fā)折算系數(shù)Kc4.2.1計算年徑流與實測年徑流的比較結(jié)果及誤差分析采用試算法，設置Kc取值在0.9-4.9之間，步長為

18、0.001，分別用1987-1992年的資料數(shù)據(jù)進行計算，得到各年全年徑流量理論計算值和實測值的相對誤差值隨著Kc的變化過程如下圖1：圖1 各年全年徑流量理論計算值和實測值的相對誤差值隨著Kc的變化自左至右分別為1991、1988、1992、1989、1987、1990年，各年相對誤差最小時，Kc取值如下表：表1 各年Kc最佳取值年份Kc最佳取值相對誤差19871.4830.0001333619881.2270.00029593519891.4235.32E-0519902.5650.00013631819911.0820.00010660919921.4180.000149347可以看出，1

19、990年的資料不可取，因此舍棄不用。4.2.2 蒸發(fā)折算系數(shù)Kc值的優(yōu)選方法在流域整體的計算徑流量與實測徑流量的相對誤差滿足5%以內(nèi)的前提下，盡可能使1987-1989年連續(xù)三年的相對誤差均在10%以內(nèi)，并使得各年的相對誤差差別盡可能小。根據(jù)相對誤差規(guī)律得以下三種方法，同時說明：通過人為觀察、比較，人工選取所給Kc的最優(yōu)值；計算1987-1989年三年的相對誤差累計值，相對誤差累計值在一定程度上反映了Kc對計算徑流量造成的影響，相對誤差累計值越小越好；計算1987-1989年三年相對誤差的均方差，它在一定程度上說明了各相對誤差偏離平均相對誤差的程度，也反映了Kc取值造成誤差的穩(wěn)定程度，相對誤差

20、的均方差越小越好。確定Kc的取值：縮小Kc的取值范圍，取Kc=0.9-1.3，步長為0.001進行計算，得到各年全年徑流量理論計算值和實測值的相對誤差值隨著Kc的變化過程如下圖2.圖2 各年全年徑流量理論計算值和實測值的相對誤差值隨著Kc的變化表2 各年Kc最佳取值年份Kc最佳取值相對誤差19871.30.07413011719881.2270.00029593519891.30.06879280519911.0820.00010660919921.30.049825263因為1990年數(shù)據(jù)舍棄，原定4年率定，改為采用1987-1989年資料進行率定，得到KC取1.263-1.430之間數(shù)據(jù)時

21、，誤差均在5%以內(nèi)，其中，Kc取值為1.3450時，誤差最小。4.3 模型檢驗及評價4.3.1 模型檢驗結(jié)果通過對Kc取值在1.263-1.430之間進行變化繪制日徑流量理論計算值和實測值的對比，得到Kc取值為1.300時，擬合較好，此時各年相對誤差如下表3.表3 Kc=1.300時的各年份相對誤差值年份878889909192相對誤差0.0741 -0.0169 0.0983 0.5440 -0.0754 0.0715 4.3.2 1987-1992年各年計算徑流與實測徑流的擬合結(jié)果Kc取1.300時，各年逐日徑流理論計算值和實測值對比圖見下圖3-8. 圖3 1987年逐日徑流理論計算值和實

22、測值對比圖圖4 1988年逐日徑流理論計算值和實測值對比圖圖5 1989年逐日徑流理論計算值和實測值對比圖圖6 1990年逐日徑流理論計算值和實測值對比圖圖7 1991年逐日徑流理論計算值和實測值對比圖圖8 1992年逐日徑流理論計算值和實測值對比圖以上各圖是在日徑流理論計算時未計算基流的情況下得到的，可見其整個起伏趨勢還是相當契合的，但局部仍存在差異，基本可以滿足常次預報的要求。因此確定Kc取1.300。4.4誤差來源設計的蓄滿產(chǎn)流模型結(jié)構(gòu)與流域的實際產(chǎn)流過程和規(guī)律不完全相符， 出現(xiàn)的問題以及可能誤差影響因素包括：（1）1990年數(shù)據(jù)計算出的結(jié)果誤差過大，可能是因為1990年數(shù)據(jù)存在問題。（

23、2）各年先對誤差均已經(jīng)很小，但不能全部達到5%的范圍之內(nèi)，可能是因為調(diào)試不夠，或者是因為模型假設與市級的情況不盡相同。（3）逐日徑流計算理論值和實測值之間存在較大的相對誤差，可能是因為未對基流進行計算，模型本身精度也有限，導致這種誤差.。（4）由于流域地理、氣候、氣象、水文條件上與模型假設條件存在一定程度上的差異，導致計算結(jié)果存在誤差。4.5模型的應用-暴雨預報對2004年暴雨過程進行洪水預報。運用單位線法和出流系數(shù)法分別計算直接徑流出流量和地下徑流出流量，兩者之和即為總的流量，其中Kc取1.300。具體結(jié)果見下表1。表4 2004年暴雨過程時間月日時RRsRg地下徑流Qg直接徑流Qs總徑流Q

24、923123.70 2.73 0.97 55.30 0.00 55.30 154.30 3.18 1.12 55.35 10.92 66.27 183.64 2.35 1.28 55.57 34.58 90.15 218.44 6.87 1.57 56.11 70.38 126.49 2424.95 22.21 2.74 57.96 115.01 172.98 24334.24 31.24 3.00 60.07 228.10 288.17 636.53 33.53 3.00 62.13 454.13 516.26 96.04 3.04 3.00 64.14 781.11 845.25 1228

25、.75 25.75 3.00 66.11 984.05 1050.16 1540.88 37.88 3.00 68.03 1095.58 1163.62 1892.05 89.05 3.00 69.91 1106.62 1176.53 2182.33 79.33 3.00 71.76 1495.12 1566.88 2446.05 43.05 3.00 73.56 2053.13 2126.69 25354.26 51.26 3.00 75.32 2673.36 2748.68 648.20 45.20 3.00 77.04 2953.15 3030.19 914.31 11.31 3.00

26、78.72 2909.37 2988.10 123.69 0.69 3.00 80.37 2668.85 2749.23 156.16 3.16 3.00 81.98 2147.52 2229.50 180.45 0.00 0.45 80.69 1500.58 1581.27 210.00 0.00 0.00 78.91 964.96 1043.87 240.00 0.00 0.00 77.17 540.63 617.81 2630.00 0.00 0.00 75.48 254.57 330.05 60.00 0.00 0.00 73.82 122.02 195.84 90.00 0.00 0

27、.00 72.19 50.44 122.64 120.00 0.00 0.00 70.60 12.66 83.26 150.00 0.00 0.00 69.05 3.50 72.55 180.00 0.00 0.00 67.53 1.58 69.11 210.00 0.00 0.00 66.05 0.00 66.05 2464.59 0.00 64.59 27363.17 0.00 63.17 661.78 0.00 61.78 960.42 0.00 60.42 1259.09 0.00 59.09 1557.79 0.00 57.79 1856.52 0.00 56.52 2155.28

28、0.00 55.28 2454.06 0.00 54.06 28352.87 0.00 52.87 將直接徑流、地下徑流出流過程和總流量過程繪制出洪水流量過程線，見下圖9.圖9 洪水流量過程線第五章 總結(jié)和心得此次課程設計，做了很久，期間碰到諸多棘手的問題，終于一一解決，得以完成，雖然還有些不盡人意，但總體上還是較好得完成了此次課程設計的各項要求。通過此次課程設計，收獲良多。首先是對于水文預報這門課程，為了完成這次課設，對課本進行了深入系統(tǒng)的復習，尤其是與新安江模型相關的內(nèi)容，使得對這門課程的掌握更加扎實牢固，理解也更加深入；其二，本次課設采用的是matlab軟件進行數(shù)據(jù)處理，這也是基于數(shù)據(jù)較

29、多，計算分析比較復雜的緣故，利用原有的一點兒編程基礎，在這期間進一步不斷地學習，對該軟件的使用能力有了很大程度的增強，這對以后的學習工作都是大有裨益的。這次課設的意義就是在于不斷逼迫自己去學習更多的新東西，并把以前學過的東西進行整合貫通，達到提升自身水平的效果；其三，在做課程設計的過程中，老師對我們進行了多次耐心認真的輔導，同學們之間也不斷進行互相的交流，也得到了很多的樂趣?？偟膩碚f，這次課程設計受益匪淺。這是畢業(yè)設計之前的最后一次課設，也是由書本上知識運用到實踐中的一次嘗試，讓我們意識到自身知識的匱乏和有限，長嘆“書到用時方恨少”，同時這也是一個有趣的過程。試想，通過自己的不懈努力終于把它攻

30、克，把成果展現(xiàn)在自己眼前時的喜悅，又有什么能比得上？感謝老師，感謝此次課設。附件：1.1987-1992年連續(xù)六年，前三年率定期和后兩年檢驗期的計算徑流量與實測徑流量絕對誤差和相對誤差計算程序：%*% 1987-1989年計算結(jié)果 % %*%clcclear ZL=load('1987到1989資料.txt');%ZL表示導入的數(shù)據(jù)WM=140;UM=20;LM=60;DM=60;B=0.3;C=0.16; IM=0.002;%流域平均張力水容量WM（mm),上層張力水容量UM（mm),下土層張力水容量LM（mm),深層張力水容量DM(mm),張力水蓄水容量曲線方次B，深層蒸散

31、發(fā)折算系數(shù)C，%不透水面積占全流域面積的比例IMQ=ZL(:,3); %日徑流量實測值Q（m3/s）E0=ZL(:,4); %日蒸發(fā)量蒸發(fā)皿實測值（mm）P1=ZL(:,5);P2=ZL(:,6);P3=ZL(:,7);P4=ZL(:,8); %P1,P2,P3,P4分別表示流域四個地區(qū)的日降雨量實測值（mm）sumQ=0; %sumQ表示年徑流量實測值（mm）sumR=zeros(4000，1); h=length(ZL);%求出數(shù)據(jù)的天數(shù)R=zeros(h,1) ; %R表示日徑流量的理論計算值EU=zeros(1,h); for i = 1:h %流域平均降雨量計算P(i) = 0.33

32、 * P1(i) + 0.14 * P2(i) + 0.33 * P3(i) + 0.2 * P4(i); %逐日降雨深的計算sumQ = sumQ + Q(i) * 24 * 3.6 / 553; %實際測得年徑流量的計算end%求逐日降雨深和實測全年徑流量W(1) = 110; WU(1) = 10; WL(1) = 40; WD(1) = 60; %初始的土壤總，上，下，深層土壤含水量%流域三層蒸發(fā)計算WMM = WM * (1 + B); %初始土壤含水量的確定a(1) = WMM * (1 - (1 - (W(1) / WM) (1 / (1 + B); for j=1:4000 K

33、c(j)=0.9+0.001*j; %設定Kc值在0.9到4.9的范圍內(nèi)變動，步長取0.001 for i = 1:h EP(i) = E0(i) * Kc(j); %流域逐日蒸散發(fā)深if WU(i) + P(i) >= EP(i) EU(i) = EP(i);EL(i) = 0;ED(i) = 0;endif WU(i) + P(i) < EP(i) if WL(i) >= C * LM EU(i) = WU(i) + P(i); EL(i) = (EP(i) - EU(i) * WL(i) / LM; ED(i) = 0; elseif WL(i) < C * LM

34、 & WL(i) >= C * (EP(i) - EU(i) EU(i) = WU(i) + P(i); EL(i) = (EP(i) - EU(i) * C; ED(i) = 0; elseif WL(i) < C * (EP(i) - EU(i) EU(i) = WU(i) + P(i); EL(i) = WL(i); ED(i) = (EP(i) - EU(i) * C - EL(i); end endE(i) = EU(i) + EL(i) + ED(i);%求得各層的蒸發(fā)量 PE(i) = P(i) - E(i); %流域產(chǎn)流計算凈降雨量PE%三層蒸散發(fā)計算，求總

35、蒸發(fā)量和流域凈降水量if PE(i) > 0 %當產(chǎn)流時 if PE(i) + a(i) < WMM R(i) = PE(i) + W(i) - WM + WM * (1 - (PE(i) + a(i) / WMM) (B + 1); W(i + 1) = W(i) + PE(i) - R(i); a(i + 1) = PE(i) + a(i); elseif PE(i) + a(i) >= WMM R(i) = PE(i) + W(i) - WM; W(i + 1) = WM; a(i + 1) = WMM; end end %產(chǎn)流計算完畢if WU(i) + P(i) -

36、 EU(i) - R(i) <= UM WU(i + 1) = WU(i) + P(i) - EU(i) - R(i); WL(i + 1) = WL(i) - EL(i); WD(i + 1) = WD(i) - ED(i); else WU(i + 1) = UM; if WL(i) - EL(i) + (WU(i) + P(i) - EU(i) - R(i) - UM) <= LM WL(i + 1) = WL(i) - EL(i) + (WU(i) + P(i) - EU(i) - R(i) - UM); WD(i + 1) = WD(i) - ED(i); else WL

37、(i + 1) = LM; if WD(i) - ED(i) + WL(i) - EL(i) + (WU(i) + P(i) - EU(i) - R(i) - UM) - LM <= DM WD(i + 1) = WD(i) - ED(i) + WL(i) - EL(i) + (WU(i) + P(i) - EU(i) - R(i) - UM) - LM; else WD(i + 1) = DM; end endend%此處計算的是降雨補充土壤含水量，由前一天的土壤含水量推求第二天的土壤含水量，補充來源為降雨減去蒸散發(fā)減去徑流量，順序為上、下、深依此補充。if PE(i) <= 0

38、 %當不產(chǎn)流時 R(i) = 0; W(i + 1) = W(i) + PE(i); a(i + 1) = WMM * (1 - (1 - W(i + 1) / WM) (1 / (1 + B);end sumR(j) = sumR(j) + R(i); end end%選出誤差最小的年份和差額cha=zeros(4000,1);for j=1:4000kc(j)=0.9+0.001*j;cha(j)=abs(sumR(j)-sumQ)/sumQ); %cha表示年徑流量理論值和實測值的相對誤差end plot(kc,cha); %計算1989年和1990年相對誤差隨著Kc值變化圖xlabel

39、('Kc值');ylabel('相對誤差值'); %坐標軸表示對象標簽%1987-1989cha=zeros(4000,1); %直接篩選出1987年誤差最小的年份和差額i=1;for j=1:4000cha(j)=abs(sumR(j)-sumQ)/sumQ); if cha(j)<=0.05KC(i)=kc(j);i=i+1 ; endend %由此得到KC取1.263-1.430之間數(shù)據(jù)時，誤差均在5%以內(nèi)y,l=min(cha); %y是誤差值，l是誤差最小時j的取值Kc8=0.9+l*0.001; %Kc的最優(yōu)取值為1.34502.1987-19

40、92年六年日模產(chǎn)流量計算程序：此處僅以1987年為例%*% 1987年日徑流量理論值和實測值的檢驗 % %*%計算步驟%首先進行三層蒸發(fā)計算，得到凈降雨量，然后進行徑流劃分，然后計算次日到初始土壤含水量%然后根據(jù)劃分出的各時段的地表和地下徑流，求取洪水過程clcclear %清空數(shù)據(jù) hs=load('1987年資料.txt');WM=140;UM=20;LM=60;DM=60;B=0.3;C=0.16; IM=0.002;Fc=24;Cg=0.978; %導入常數(shù)值;RR=hs(:,3);E0=hs(:,4);P1=hs(:,5);P2=hs(:,6);P3=hs(:,7);

41、P4=hs(:,8);Kc=1.440; %確定Kc取值m=length(hs); %確定暴雨持續(xù)時間R=zeros(1,m); EU=zeros(1,m);for i = 1:m %流域平均降雨量計算P(i) = 0.33 * P1(i) + 0.14 * P2(i) + 0.33 * P3(i) + 0.2 * P4(i); %逐日降雨深的計算end%求逐日降雨深 W(1) = 110; WU(1) = 10; WL(1) = 40; WD(1) = 60; %流域三層蒸發(fā)計算，先給出初始 上、下、深層土壤含水量WMM = WM * (1 + B);%初始土壤含水量的確定a(1) = WM

42、M * (1 - (1 - (W(1) / WM) (1 / (1 + B); for i = 1:m EP(i) = E0(i) * Kc;%逐日蒸散發(fā)深if WU(i) + P(i) >= EP(i) EU(i) = EP(i);EL(i) = 0;ED(i) = 0;endif WU(i) + P(i) < EP(i) if WL(i) >= C * LM EU(i) = WU(i) + P(i); EL(i) = (EP(i) - EU(i) * WL(i) / LM; ED(i) = 0; elseif WL(i) < C * LM & WL(i) &

43、gt;= C * (EP(i) - EU(i) EU(i) = WU(i) + P(i); EL(i) = (EP(i) - EU(i) * C; ED(i) = 0; elseif WL(i) < C * (EP(i) - EU(i) EU(i) = WU(i) + P(i); EL(i) = WL(i); ED(i) = (EP(i) - EU(i) * C - EL(i); end endE(i) = EU(i) + EL(i) + ED(i);%求得各層的蒸發(fā)量 PE(i) = P(i) - E(i); %流域產(chǎn)流計算凈降雨量PE%三層蒸散發(fā)計算，求總蒸發(fā)量和流域凈雨量if PE

44、(i) > 0 %當產(chǎn)流時 if PE(i) + a(i) < WMM R(i) = PE(i) + W(i) - WM + WM * (1 - (PE(i) + a(i) / WMM) (B + 1); W(i + 1) = W(i) + PE(i) - R(i); a(i + 1) = PE(i) + a(i); elseif PE(i) + a(i) >= WMM R(i) = PE(i) + W(i) - WM; W(i + 1) = WM; a(i + 1) = WMM; end end %產(chǎn)流計算完畢if WU(i) + P(i) - EU(i) - R(i) &

45、lt;= UM WU(i + 1) = WU(i) + P(i) - EU(i) - R(i); WL(i + 1) = WL(i) - EL(i); WD(i + 1) = WD(i) - ED(i); else WU(i + 1) = UM; if WL(i) - EL(i) + (WU(i) + P(i) - EU(i) - R(i) - UM) <= LM WL(i + 1) = WL(i) - EL(i) + (WU(i) + P(i) - EU(i) - R(i) - UM); WD(i + 1) = WD(i) - ED(i); else WL(i + 1) = LM; i

46、f WD(i) - ED(i) + WL(i) - EL(i) + (WU(i) + P(i) - EU(i) - R(i) - UM) - LM <= DM WD(i + 1) = WD(i) - ED(i) + WL(i) - EL(i) + (WU(i) + P(i) - EU(i) - R(i) - UM) - LM; else WD(i + 1) = DM; end endend%此處計算的是降雨補充土壤含水量，由前一天的土壤含水量推求第二天的土壤含水量，補充來源為降雨減去蒸散發(fā)減去徑流量，順序為上、下、深依此補充。if PE(i) <= 0 %當不產(chǎn)流時 R(i) = 0; W(i + 1) = W(i) + PE(i); a(i + 1) = WMM * (1 - (1 - W(i + 1) / WM) (1 / (1 + B);end end%以上產(chǎn)流計算完畢f(xié)or i=1:m s(i)=i;endplot(s ,RR,s,R); %計算1989年和1990年相對誤差隨著Kc值變化圖xlabel('天數(shù)');ylabel('流量值'); %坐標軸表示對象標簽Rr=0;for i=1:mRr=Rr+RR(i)/m; %Rr是實測值的均值endX=0;Y=0;for i=1:mX=X+(R(i)-RR(i)2