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

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

3、構，透水性好。臺地、丘陵多生長松、杉、樟等高大喬木；平原則以種植農作物和經濟作物為主，植被良好。流域上游有寶口水文站，流域面積553km2,占白盆珠水庫壩址以上集雨面積的64.6%。白盆珠水庫有10年逐日入庫流量資料、逐日蒸發資料和時段入庫流量資料：流域內有7個雨量站，其中寶口以上有4個。雨量站分布較均勻.有10年逐日降水資料和時段降水資料；寶口水文站具有10年以上水位、流量資料；流域屬山區性小流域且受到地形、地貌等下墊面條件影響，洪水陡漲緩落，匯流時間一般2-3h,有時更短；一次洪水總歷時25d。圖2-1該地區水文站分布圖2.2 數據資料計算流域面積為553km20流域內有四個雨量站，權重系

4、數分別為0.33、0.14、0.33、0.20。(3)日產流模型數據，具體見資料文件名：87-92data.xls,數據格式為:T(i)Q(i)E(i)P1(i)P2(i)P3(i)P4(i)(4)暴雨預報的數據，見表2-1表2-12004年暴雨過程數據表時間蒸散發(mm)降雨量(mm禾多布馬山高潭寶口T(i)E(i)P1(i)P2(i)P3(i)P4(i)2004-9-2312:001.36.29.921.617.32004-9-2315:001.37.61620.612.62004-9-2318:001.36.26.414.915.92004-9-2321:001.38.817.229.4

5、18.52004-9-240:001.22534.835.324.62004-9-243:000.929.929.243.937.82004-9-246:000.938.624.846.9332004-9-249:000.96.97.56.112.32004-9-2412:000.928.329.934.228.52004-9-2415:000.925.642.739.875.42004-9-2418:000.993.9137.612413.22004-9-2421:000.985.390.88575.92004-9-250:000.851.547.749.238.52004-9-253:00

6、1.139.870.342.197.72004-9-256:001.143.247.361.545.92004-9-259:001.120.513.315.813.12004-9-2512:001.110.581.83.32004-9-2515:001.17.48.47.610.92004-9-2518:001.11.82.82.14.62004-9-2521:001.10.200.302004-9-260:001.200002004-9-263:002.100002004-9-266:002.100002004-9-269:002.100002004-9-2612:002.100002004

7、-9-2615:002.100002004-9-2618:002.100002004-9-2621:0020000(5)計算參數數據，見表2-2表2-2計算參數表計算年份參數初始張力水蓄量1989-1990WmUmLmDmWWUWLWD140206060110104060BCFcIM0.20.16220.001(6)流域單位線單位線過程(m3/s)為：0,40,80,130,100,80,48,20,10,5,0(7)地下徑流匯流3Cg=0.978,Qg=55.3m3/s第三章計算公式該流域海洋性氣候顯著、氣候溫和、雨量豐沛，多年平均降雨量為1800mm徑流系數0.5-0.7,土壤主要有黃壤和

8、砂壤，層次結構明顯，透水性好，植被覆蓋度高，地勢平坦，由此可初步判定該流域的產流機制為蓄滿產流模式。3.1 產流計算3.1.1 蒸散發計算根據流域蓄滿產流特點，蒸散發計算采用的是三層蒸散發計算模式。三層蒸發模式的具體計算如下：1）當WU+PEpEU=Ep,EL=0,ED=0;2)當WU+PErWL>CXWLM,01.Q 。EU=WU+P,EL=(Ep-EU)XWL/WLM,ED=0;3)當WU+PErC(EP-EU)&WL<CXWLM,EU=WU+P,EL=C(EyEU),ED=0;4)當WU+PEP,WKC(EyEU),EU=WU+P,EL=WL,ED=C(E>EU

9、)-EL.式中：WUfe上層土壤蓄水量，WL*下層土圖2-1包氣帶蓄水容量曲線壤蓄水量，EU為上層土壤蒸發量，EL為下層土壤蒸發量，ED為深層土壤蒸發量，P為流域平均降雨量，Ep為流域平均蒸發能力，C為深層蒸散發擴散系數，WLM為下層張力水蓄水容量。3.1.2 產流量計算根據流域特點，產流量計算系根據蓄滿產流理論得出的。蓄滿產流，即任一地點上，土壤含水量達田間持水量前，降雨量全部補充土壤含水量，不產流；當土壤蓄滿后，其后續降雨量全部產生徑流。流域內各點包氣帶的蓄水容量是不同的，將各點包氣帶蓄水容量從小到大排列，以包氣帶達到田間持水量時的土壤含水量WM為縱坐標，以流域內小于等于該WM的面積占全流

10、域的面積比a為橫坐標，所繪的曲線稱為流域蓄水容量曲線。1/(b+1)a=WMM(1-(1-W/WM)1PE>Q則產流；否則不產流。產流時：1)當PE+aCWMMb+1R=PE+W-WM+WM(PE+a)/WMM)2) 當PE+a>WM：MR=PE+W-WM式中：PE為扣除蒸發量后的降雨量，a為土壤含水量W對應的土壤水深，WM為流域平均蓄水容量，WMMJ流域各地點包氣帶蓄水容量的最大值，b為流域包氣帶蓄水容量分布的不均勻指數，R為流域產流量。3.1.3 二水源劃分流域坡地上的降雨產流量因產流過程的條件和運動路徑不同，受流域的調蓄作用不同，各徑流成分在流量過程線上的反應是不一樣的。在實

11、際工作中，常需按各種徑流成分分別計算或模擬，因為要對產流量進行水源劃分。直接徑流和地下徑流水源劃分如下：1)當PE<=FC寸：RS=0.0RG=R2)當PE>FC時：RG=FC*R/PERS=R-RG式中：FC為穩定下滲率，RS為直接徑流，RG為地下徑流。3.1.4 各層蓄水量計算降雨補充土壤含水量，由前一天的土壤含水量推求第二天的土壤含水量，補充來源為降雨減去蒸散發減去徑流量，順序為上、下、深層依此補充。三層蓄水量變化的具體計算如下：1) WUi+Pi-EUi-Ri<=UM，WUi+1=WUi+Pi-EUi-Ri;WLi+1=WLi-ELi;WDi+1=WDi-EDi;2)

12、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<=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-E

13、Li+(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匯流計算根據流域凈雨和流域徑流單位線，采用卷積的差分形式算出流域出口的流量過程。直接徑流匯流可根據該流域的時段單位線推求，地下徑流匯流由線性水庫演算法推求。1)直接徑流匯流計算公式：QS(i尸RS(i)XUH;式中：UH為該流域的單位線2)地下徑流匯流計算公式：QG(i尸CGXQG(i-1)+(1-CG)xRG(i)xuQ(i

14、尸QS(i)+QG(i)/行一流域面積F(km2)式中：U為單包轉換系數，U=3.6*：t(h)3.3模型參數1) Kc：蒸散發能力折算系數，它主要反映流域平均高程與蒸發站高程之間差別的影響和蒸發皿散發與路面蒸散發問差別的影響；2) WM流域平均張力水容量，它表示流域蓄滿的標準；3) WUMh層張力水蓄水容量，它包括了植物截留量；4) WLM下層張力水蓄水容量；5) b:流域包氣帶蓄水容量分布的不均勻指數，在一般情況下其取值與單元流域面積有關；6) C:深層蒸散發擴散系數，它主要取決于流域內深根植物的覆蓋范圍。7) IM:不透水面積占全流域面積的比例，它的值可由大比例尺的地形圖，通過地理信息系

15、統現代技術量測出來，也可用歷史上干旱期小洪水資料來分析。第四章模型檢驗、結果評定及分析水文預報是一項直接服務于國家安全和國民經濟建設的不可或缺的重要基本工作，是幫助人類有效地預防洪水、減少洪災損失，有效利用水資源的非工程措施之一。隨著經濟、社會發展及其全球化進程的需要，水文預報的服務面進一步拓展，對水文預報提出了更高的要求。水文預報結果的準確率與可信程度是衡量服務質量的前提，為了更好地為國家安全和國民經濟建設服務，必須對水文預報結果的可靠性和有效性進行評定和檢驗。4.1 產流模式的檢驗定性分析該流域集雨面積856km2。流域地處南方，海洋性氣候顯著，氣候溫和，雨量豐沛。暴雨成因主要是鋒面雨和臺

16、風雨，常受熱帶風暴影響。降雨年際間變化大，年內分配不均，多年平均降雨量為1800mmi實測年最大降雨量為3417mmi汛期49月降雨量占年降雨量的81%左右：徑流系數0.50.7。流域內地勢平緩，土壤主要有黃壤和砂壤，具有明顯的腐殖層，淀積層和母質土等層次結構，透水性好。臺地、丘陵多生長松、杉、樟等高大喬木；平原則以種植農作物和經濟作物為主，植被良好。流域上游有一水文立控制流域面積553km2,占流域集雨面積的64.6%。該水文站以上有4個雨量站。雨量站分布較均勻，有10年逐日降水資料和時段降水資料；該水文站具有10年以上水位、流量資料；流域屬山區性小流域且受到地形、地貌等下墊面條件影響，洪水

17、陡漲緩落，匯流時間一般2-3h,有時更短；一次洪水總歷時25d。由流域概況可以看出，該地區屬于濕潤多雨地區，雨量集中，地勢平坦，土壤層容易蓄滿，而且有場系列的降雨資料和水位流量資料，綜合從氣候條件、土壤狀況、植被組成以及洪水機制看，新安江模型適用于該地區的水文預報。4.2 率定期優選蒸發折算系數Kc4.2.1 計算年徑流與實測年徑流的比較結果及誤差分析采用試算法，設置Kc取值在0.9-4.9之間，步長為0.001，分別用1987-1992年的資料數據進行計算,得到各年全年徑流量理論計算值和實測值的相對誤差值隨著Kc的變化過程如下圖1:圖1各年全年徑流量理論計算值和實測值的相對誤差值隨著Kc的變

18、化自左至右分別為1991、1988、1992、1989、1987、1990年，各年相對誤差最小時，Kc取值如下表：表1各年Kc最佳取值年份Kc最佳取值相對誤差19871.4830.0001333619881.2270.00029593519891.4235.32E-0519902.5650.00013631819911.0820.00010660919921.4180.000149347可以看出，1990年的資料不可取，因此舍棄不用4.2.2 蒸發折算系數Kc值的優選方法在流域整體的計算徑流量與實測徑流量的相對誤差滿足5%Z內的前提下，盡可能使1987-1989年連續三年的相對誤差均在10%Z

19、內，并使得各年的相對誤差差別盡可能小。根據相對誤差規律得以下三種方法，同時說明：通過人為觀察、比較，人工選取所給Kc的最優值；計算1987-1989年三年的相對誤差累計值，相對誤差累計值在一定程度上反映了Kc對計算徑流量造成的影響，相對誤差累計值越小越好；計算1987-1989年三年相對誤差的均方差，它在一定程度上說明了各相對誤差偏離平均相對誤差的程度，也反映了Kc取值造成誤差的穩定程度，相對誤差的均方差越小越好。4.2.3 確定Kc的取化縮小Kc的取值范圍，取Kc=0.9-1.3,步長為0.001進行計算，得到各年全年徑流量理論計算值和實測值的相對誤差值隨著Kc的變化過程如下圖2.圖2各年全

20、年徑流量理論計算值和實測值的相對誤差值隨著Kc的變化表2各年Kc最佳取值年份Kc最佳取值相對誤差19871.30.07413011719881.2270.00029593519891.30.06879280519911.0820.00010660919921.30.049825263因為1990年數據舍棄，原定4年率定，改為采用1987-1989年資料進行率定，得到KC取1.263-1.430之間數據時，誤差均在5犯內，其中，Kc取值為1.3450時，誤差最小。4.3 模型檢驗及評價4.3.1 模型檢驗結果通過對Kc取值在1.263-1.430之間進行變化繪制日徑流量理論計算值和實測值的對比，

21、得到Kc取值為1.300時，擬合較好，此時各年相對誤差如下表3.表3Kc=1.300時的各年份相對誤差值年份878889909192相對誤差0.0741-0.01690.09830.5440-0.07540.07154.3.21987-1992年各年計算徑流與實測徑流的擬合結果Kc取1.300時，各年逐日徑流理論計算值和實測值對比圖見下圖3-8.圖31987年逐日徑流理論計算值和實測值對比圖圖41988年逐日徑流理論計算值和實測值對比圖圖51989年逐日徑流理論計算值和實測值對比圖圖61990年逐日徑流理論計算值和實測值對比圖圖71991年逐日徑流理論計算值和實測值對比圖圖81992年逐日徑流

22、理論計算值和實測值對比圖以上各圖是在日徑流理論計算時未計算基流的情況下得到的，可見其整個起伏趨勢還是相當契合的，但局部仍存在差異，基本可以滿足常次預報的要求。因止匕確定Kc取1.300。4.4 誤差來源設計的蓄滿產流模型結構與流域的實際產流過程和規律不完全相符，出現的問題以及可能誤差影響因素包括:(1)1990年數據計算出的結果誤差過大，可能是因為1990年數據存在問題。(2)各年先對誤差均已經很小，但不能全部達到5%勺范圍之內，可能是因為調試不夠，或者是因為模型假設與市級的情況不盡相同。(3)逐日徑流計算理論值和實測值之間存在較大的相對誤差，可能是因為未對基流進行計算，模型本身精度也有限，導

23、致這種誤差.o(4)由于流域地理、氣候、氣象、水文條件上與模型假設條件存在一定程度上的差異，導致計算結果存在誤差。4.5模型的應用-暴雨預報對2004年暴雨過程進行洪水預報。運用單位線法和出流系數法分別計算直接徑流出流量和地下徑流出流量，兩者之和即為總的流量，其中Kc取1.300。具體結果見下表1。表42004年暴雨過程時間月日時RRsRg地下徑流Qg直接徑流Qs總徑流Q923123.702.730.9755.300.0055.30154.303.181.1255.3510.9266.27183.642.351.2855.5734.5890.15218.446.871.5756.1170.38

24、126.492424.9522.212.7457.96115.01172.9824334.2431.243.0060.07228.10288.17636.5333.533.0062.13454.13516.2696.043.043.0064.14781.11845.251228.7525.753.0066.11984.051050.161540.8837.883.0068.031095.581163.621892.0589.053.0069.911106.621176.532182.3379.333.0071.761495.121566.882446.0543.053.0073.562053.

25、132126.6925354.2651.263.0075.322673.362748.68648.2045.203.0077.042953.153030.19914.3111.313.0078.722909.372988.10123.690.693.0080.372668.852749.23156.163.163.0081.982147.522229.50180.450.000.4580.691500.581581.27210.000.000.0078.91964.961043.87240.000.000.0077.17540.63617.812630.000.000.0075.48254.5

26、7330.0560.000.000.0073.82122.02195.8490.000.000.0072.1950.44122.64120.000.000.0070.6012.6683.26150.000.000.0069.053.5072.55180.000.000.0067.531.5869.11210.000.000.0066.050.0066.052464.590.0064.5927363.170.0063.17661.780.0061.78960.420.0060.421259.090.0059.091557.790.0057.791856.520.0056.522155.280.0

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

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

29、的喜悅，又有什么能比得上？感謝老師，感謝此次課設。附件：1.1987-1992年連續六年，前三年率定期和后兩年檢驗期的計算徑流量與實測徑流量絕對誤差和相對誤差計算程序：%*%1987-1989年計算結果%*%clcclearZL=load（'1987到1989資料.txt'）;%ZL表示導入的數據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,深層蒸散發折算系數C，%不透水面積占全流域面積

30、的比例IMQ=ZL(:,3);%日徑流量實測值Q(m3/s)E0=ZL(:,4);%日蒸發量蒸發皿實測值(mm)P1=ZL(:,5);P2=ZL(:,6);P3=ZL(:,7);P4=ZL(:,8);%P1,P2,P3,P4分別表示流域四個地區的日降雨量實測值(mm)sumQ=0;%sumQ表示年徑流量實測值(mm)sumR=zeros(4000，1);h=length(ZL);%求出數據的天數R=zeros(h,1);%R表示日徑流量的理論計算值EU=zeros(1,h);fori=1:h%流域平均降雨量計算P(i)=0.33*P1(i)+0.14*P2(i)+0.33*P3(i)+0.2*

31、P4(i);%逐日降雨深的計算sumQ=sumQ+Q(i)*24*3.6/553;%實際測得年徑流量的計算end%求逐日降雨深和實測全年徑流量W(1)=110;WU(1)=10;WL(1)=40;WD(1)=60;%初始的土壤總，上，下，深層土壤含水量%流域三層蒸發計算WMM=WM*(1+B);%初始土壤含水量的確定a(1)=WMM*(1-(1-(W(1)/WM)A(1/(1+B);forj=1:4000Kc(j)=0.9+0.001*j;%設定Kc值在0.9到4.9的范圍內變動，步長取0.001fori=1:hEP(i)=E0(i)*Kc(j);%流域逐日蒸散發深ifWU(i)+P(i)&g

32、t;=EP(i)EU(i)=EP(i);EL(i)=0;ED(i)=0;endifWU(i)+P(i)<EP(i)ifWL(i)>=C*LMEU(i)=WU(i)+P(i);EL(i)=(EP(i)-EU(i)*WL(i)/LM;ED(i)=0;elseifWL(i)<C*LM&WL(i)>=C*(EP(i)-EU(i)EU(i)=WU(i)+P(i);EL(i)=(EP(i)-EU(i)*C;ED(i)=0;elseifWL(i)<C*(EP(i)-EU(i)EU(i)=WU(i)+P(i);EL(i)=WL(i);ED(i)=(EP(i)-EU(i)*

33、C-EL(i);endendE(i)=EU(i)+EL(i)+ED(i);%求得各層的蒸發量PE(i)=P(i)-E(i);%流域產流計算凈降雨量PE%三層蒸散發計算，求總蒸發量和流域凈降水量ifPE(i)>0%當產流時ifPE(i)+a(i)<WMMR(i)=PE(i)+W(i)-WM+WM*(1-(PE(i)+a(i)/WMM)A(B+1);W(i+1)=W(i)+PE(i)-R(i);a(i+1)=PE(i)+a(i);elseifPE(i)+a(i)>=WMMR(i)=PE(i)+W(i)-WM;W(i+1)=WM;a(i+1)=WMM;endend%產流計算完畢if

34、WU(i)+P(i)-EU(i)-R(i)<=UMWU(i+1)=WU(i)+P(i)-EU(i)-R(i);WL(i+1)=WL(i)-EL(i);WD(i+1)=WD(i)-ED(i);elseWU(i+1)=UM;ifWL(i)-EL(i)+(WU(i)+P(i)-EU(i)-R(i)-UM)<=LMWL(i+1)=WL(i)-EL(i)+(WU(i)+P(i)-EU(i)-R(i)-UM);WD(i+1)=WD(i)-ED(i);elseWL(i+1)=LM;ifWD(i)-ED(i)+WL(i)-EL(i)+(WU(i)+P(i)-EU(i)-R(i)-UM)-LM<

35、;=DMWD(i+1)=WD(i)-ED(i)+WL(i)-EL(i)+(WU(i)+P(i)-EU(i)-R(i)-UM)-LM;elseWD(i+1)=DM;endendend%此處計算的是降雨補充土壤含水量，由前一天的土壤含水量推求第二天的土壤含水量，補充來源為降雨減去蒸散發減去徑流量，順序為上、下、深依此補充。ifPE(i)<=0%當不產流時R(i)=0;W(i+1)=W(i)+PE(i);a(i+1)=WMM*(1-(1-W(i+1)/WM)A(1/(1+B);endsumR(j)=sumR(j)+R(i);endend%選出誤差最小的年份和差額cha=zeros(4000,1

36、);forj=1:4000kc(j)=0.9+0.001*j;cha(j)=abs(sumR(j)-sumQ)/sumQ);%cha表示年徑流量理論值和實測值的相對誤差endplot(kc,cha);%計算1989年和1990年相對誤差隨著Kc值變化圖xlabel('Kc值');ylabel('相對誤差值');%坐標軸表示對象標簽%1987-1989cha=zeros(4000,1);%直接篩選出1987年誤差最小的年份和差額i=1;forj=1:4000cha(j)=abs(sumR(j)-sumQ)/sumQ);ifcha(j)<=0.05KC(i)=

37、kc(j);i=i+1;endend%由此得到KC取1.263-1.430之間數據時，誤差均在5%以內y,l=min(cha);%y是誤差值，l是誤差最小時j的取值Kc8=0.9+l*0.001;%Kc的最優取值為1.34502.1987-1992年六年日模產流量計算程序：此處僅以1987年為例%*%1987年日徑流量理論值和實測值的檢驗%*%計算步驟%首先進行三層蒸發計算，得到凈降雨量，然后進行徑流劃分，然后計算次日到初始土壤含水量%然后根據劃分出的各時段的地表和地下徑流，求取洪水過程clcclear%清空數據hs=load('1987年資料.txt');WM=140;UM=

38、20;LM=60;DM=60;B=0.3;C=0.16;IM=0.002;Fc=24;Cg=0.978;%導入常數值;RR=hs(:,3);E0=hs(:,4);P1=hs(:,5);P2=hs(:,6);P3=hs(:,7);P4=hs(:,8);Kc=1.440;%確定Kc取值m=length(hs);%確定暴雨持續時間R=zeros(1,m);EU=zeros(1,m);fori=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)=4

39、0;WD(1)=60;%流域三層蒸發計算，先給出初始上、下、深層土壤含水量WMM=WM*(1+B)；%初始土壤含水量的確定a(1)=WMM*(1-(1-(W(1)/WM)A(1/(1+B);fori=1:mEP(i)=E0(i)*Kc;%逐日蒸散發深ifWU(i)+P(i)>=EP(i)EU(i)=EP(i);EL(i)=0;ED(i)=0;endifWU(i)+P(i)<EP(i)ifWL(i)>=C*LMEU(i)=WU(i)+P(i);EL(i)=(EP(i)-EU(i)*WL(i)/LM;ED(i)=0;elseifWL(i)<C*LM&WL(i)>

40、;=C*(EP(i)-EU(i)EU(i)=WU(i)+P(i);EL(i)=(EP(i)-EU(i)*C;ED(i)=0;elseifWL(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);endendE(i)=EU(i)+EL(i)+ED(i);%求得各層的蒸發量PE(i)=P(i)-E(i);%流域產流計算凈降雨量PE%三層蒸散發計算，求總蒸發量和流域凈雨量ifPE(i)>0%當產流時ifPE(i)+a(i)<WMMR(i)=PE(i)+W(i)-WM+WM*(1-(PE(i

41、)+a(i)/WMM)A(B+1);W(i+1)=W(i)+PE(i)-R(i);a(i+1)=PE(i)+a(i);elseifPE(i)+a(i)>=WMMR(i)=PE(i)+W(i)-WM;W(i+1)=WM;a(i+1)=WMM;endend%產流計算完畢ifWU(i)+P(i)-EU(i)-R(i)<=UMWU(i+1)=WU(i)+P(i)-EU(i)-R(i);WL(i+1)=WL(i)-EL(i);WD(i+1)=WD(i)-ED(i);elseWU(i+1)=UM;ifWL(i)-EL(i)+(WU(i)+P(i)-EU(i)-R(i)-UM)<=LMWL

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

43、ifPE(i)<=0%當不產流時R(i)=0;W(i+1)=W(i)+PE(i);a(i+1)=WMM*(1-(1-W(i+1)/WM)A(1/(1+B);endend%以上產流計算完畢fori=1:ms(i)=i;endplot(s,RR,s,R);%計算1989年和1990年相對誤差隨著Kc值變化圖xlabel('天數');ylabel('流量值');%坐標軸表示對象標簽Rr=0;fori=1:mRr=Rr+RR(i)/m;%Rr是實測值的均值endX=0;Y=0;fori=1:mX=X+(R(i)-RR(i)A2;Y=Y+(RR(i)-Rr)A2;D

44、C=1-X/Y;%DC是確定性系數end3.暴雨預報計算程序：%計算步驟%首先進行三層蒸發計算，得到凈降雨量，然后進行徑流劃分，然后計算次日到初始土壤含水量%然后根據劃分出的各時段的地表和地下徑流，求取洪水過程%*%2004年場次暴雨的洪水預報計算%*%clcclear%清空數據hs=load('2004暴雨.txt');WM=140;UM=20;LM=60;DM=60;B=0.3;C=0.16;IM=0.002;Fc=24;Cg=0.978;%導入常數值;E0=hs(:,3);P1=hs(:,4);P2=hs(:,5);P3=hs(:,6);P4=hs(:,7);Kc=1.9

45、540;%確定Kc取值m=length(hs);%確定暴雨持續時間R=zeros(1,m);fori=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;%流域三層蒸發計算，先給出初始上、下、深層土壤含水量WMM=WM*(1+B);%初始土壤含水量的確定a(1)=WMM*(1-(1-(W(1)/WM)A(1/(1+B);fori=1:mEP(i)=E0(i)*Kc;%逐日蒸散發深ifWU(i)+P(i)>=EP(i)EU(i)=EP(i);EL(i)=0;ED(i)=0;endifWU(i)+P(i)<EP(i)ifWL(i)>=C*LMEU(i)=WU(i)+P(i);EL(i)=(EP(i)-EU(i)*WL(i)/LM;ED(i)=0;elseifWL(i)<C*LM&WL(i)>=C*(EP(i)-EU(i)EU(i)=WU(i)+P(i);EL(i)=(EP(i)-EU(i)*C;ED(i)=0;elseifWL(i)<C*(EP(i