【正文】 時(shí)間步長內(nèi)到達(dá)地面。 (516) (HarGreaves 1994)為宇宙射輻射MJ/(m2day)；為日平均氣溫（日最大氣溫與最小氣溫的均值）；為日氣溫變化 (日最大氣溫與最小氣溫的差值)。常用的方法有如下：；(如能量平衡法,空氣動力學(xué)法,能量平衡空氣動力學(xué)法和渦度相關(guān)技術(shù)等)；，利用衛(wèi)星遙感并結(jié)合模式模擬研究非均勻陸面上的蒸散發(fā)是一個(gè)新的趨勢,該項(xiàng)技術(shù)在我國剛剛起步。由上幾個(gè)式可知：=. (512)因此，臨界降水量完全由最大葉面持水深度和葉面指數(shù)兩個(gè)參數(shù)確定，于是可得植被截留量公式： (513)由于本文沒有葉面指以及最大葉面持水深度資料，而葉面指數(shù)和最大葉面持水深度如上面提到都有范圍，結(jié)合實(shí)地情況對其初始估計(jì)值如表(52)： 表(52) 葉面指數(shù)、最大葉面持水深度初始估計(jì)值　葉面指數(shù)葉面持水深度(mm)森林、灌木林7草地2作物耕地3取=，王志安（2005，Wang et al 2005） 蒸散發(fā)蒸散發(fā)包括蒸發(fā)和散發(fā)，蒸發(fā)包括水面蒸發(fā)和土壤蒸發(fā)，散發(fā)指植物散發(fā)，也就是由植物莖葉向空氣中散發(fā)水分的過程。最大葉面持水深度(mm),～。圖54 植被截留量與降水量P的關(guān)系植被蓋度(%)：指植被(包括葉、莖、枝)在單位面積內(nèi)植被的垂直投影面積所占百分比。其計(jì)算方法可用如下公式[36]:()= (59)式中是截留系數(shù)，它是降水量P的函數(shù)。，柵格單元的氣溫按下面的公式計(jì)算： (56) (57) (58) 式中：高程為Z的柵格單元的氣溫()；、回歸系數(shù)；為測站個(gè)數(shù)；測站氣溫()；測站高程(m)；個(gè)測站的平均高程(m)；個(gè)測站的平均氣溫()。，柵格單元的氣溫按下面的公式計(jì)算： (54) 式中：高程為Z的柵格單元的氣溫()；Z高程(m)；測站氣溫()；測站高程(m)。在無法取得降雨量空間分布函數(shù)的情況下，泰森多邊形法為許多水文計(jì)算選用，本文選擇泰森多邊法計(jì)算燕山水庫流域降雨分布（如圖53）。此法是將相鄰雨量站用直線連接而成為若干個(gè)三角形，然后對三角形各邊作垂直平分線，連接這些垂直平分線的交點(diǎn)，得若干個(gè)多邊形，每個(gè)多邊形內(nèi)各有一個(gè)雨量站，該多邊形即是該雨量站的控制面積。，這種方法只能用在雨量站稠密的流域，否則，這種方法不太準(zhǔn)確；，再繪制等雨量，再求得相鄰雨量線之間的面積，并計(jì)算相鄰降雨量線的平均降雨深，再計(jì)算平均降雨量，這種方法同樣需要足夠的雨量站；,降水分布模型可以通過三維空間插值法。但是雨量站網(wǎng)觀測的降雨量在空間上是非規(guī)則離散分布的，并不能完全反映實(shí)際降雨在空間上的連續(xù)分布，因此，在研究中常常要從雨量站網(wǎng)觀測資料來合理提取降雨空間分布的特征值，如流域平均雨量、面雨量或流域降雨空間分布等值線等。目前，獲取降雨資料的手段有兩種方法：雨量站觀測和雷達(dá)遙感法。圖51 柵格地表、地下水輸出輸入 步長內(nèi)柵格地表水量平衡方程 (51) (52) 步長內(nèi)柵格地下水量平衡方程 (53)計(jì)算時(shí)間步長；步長內(nèi)總降雨量；步長內(nèi)植被截流；步長內(nèi)總蒸散發(fā)量；步長內(nèi)鄰近柵格地表入流量； 步長內(nèi)柵格地表出流量； 步長內(nèi)下滲量；S步長內(nèi)柵格地表蓄水量；步長內(nèi)地下土層含水量；步長內(nèi)到達(dá)地面的雨量；步長內(nèi)鄰近柵格地下入流量；步長內(nèi)柵格地下出流量； 步長內(nèi)柵格地下蓄水量；步長內(nèi)深層下滲量； 步長內(nèi)柵格土壤蒸發(fā)量。對柵格地表、地下層的匯流計(jì)算是分別進(jìn)行產(chǎn)匯流演算的，二層之間相互聯(lián)系、又相互獨(dú)立。 柵格地表水、地下水輸出輸入原理 柵格產(chǎn)流不僅受降雨、植被截流、下滲、蒸散發(fā)等影響，同時(shí)也受周圍柵格匯流的影響，對于分布式水文模型，一般假定產(chǎn)流由地表徑流和地下徑流(主要為壤中流)組成，當(dāng)包氣帶達(dá)到田間持水量，或降雨強(qiáng)度大于下滲能力，地表就開始積水，并在重力的作用下順坡面流動，形成的坡面薄層水流就是地表徑流；同時(shí)，地下土壤的含水量由初始含水量逐漸增大，當(dāng)土壤包氣帶達(dá)到田間持水量后不但有地表徑流產(chǎn)生還會產(chǎn)生地下徑流。在傳統(tǒng)的集總式水文模型中，它能夠大大地簡化計(jì)算量，提高模型的適應(yīng)能力。目前，建立這樣復(fù)雜的水文數(shù)學(xué)物理方程還比較困難，也難應(yīng)用到復(fù)雜的大流域。，用數(shù)學(xué)物理方法對水文現(xiàn)象進(jìn)行模擬，依據(jù)物理學(xué)質(zhì)量、動量與能量守恒定律以及產(chǎn)匯流的特性，推導(dǎo)出描述地表徑流和地下徑流的微分方程組。 圖417 燕山水庫水位庫容曲線 第五章 柵格單元水文模型 單元水文模型構(gòu)建方法單元水文模型是分布式水文模型的核心部分，其內(nèi)部結(jié)構(gòu)由地表水模型、土壤水模型等構(gòu)成，所涉及到的水文過程有蒸散發(fā)、下滲、地表徑流、壤中流等。庫容曲線是水庫水位、水量擬合曲線。本文的庫容值是用用GIS軟件ArcView計(jì)算的，其原理是就是求解DEM體積，即在TIN的基礎(chǔ)上求解下表面為TIN上表面為水平面組成的空間體的體積，通過計(jì)算三棱柱體積進(jìn)行累加得到，每一個(gè)三棱柱體積計(jì)算公式為： ；；式中zi，xi，yi為TIN中一個(gè)三角形三個(gè)頂點(diǎn)的坐標(biāo)值；H為水位高程；S是三棱柱底面積也即TIN中的一個(gè)三角形的面積。但是，由于地形地勢的復(fù)雜性，隨著水位不同其值是變化的，傳統(tǒng)水庫庫容計(jì)算方法成果生產(chǎn)周期長，精度低且數(shù)據(jù)更新困難。測定和計(jì)算各特征水位下的相應(yīng)庫容對于水庫管理中實(shí)行有計(jì)劃地蓄水、泄水及編制水庫控制運(yùn)用計(jì)劃都是十分重要的。圖416就是采用以上方法提取的燕山水庫流域的流水累積量矩陣，由柵格顏色深與淺表示流水積累量的多與少，由圖可以看出流水累積量矩陣的圖和河網(wǎng)一致，流域流水積累量最大值與流域的柵格總數(shù)一致。權(quán)值矩陣是一個(gè)連續(xù)的數(shù)字矩陣，它表示了在一次暴雨中的平均降雨量，可用來計(jì)算一個(gè)流域內(nèi)的流走的降雨量的多少,并且假定所有的降雨量沒有被地表水截留、蒸發(fā)或損失，輸出的流水累積量矩陣表示了流經(jīng)每個(gè)柵格的降雨量。其基本方法是：以規(guī)則柵格表示的數(shù)字地面高程模型每點(diǎn)處有一個(gè)單位的水量，按照自然水流從高處流往低處的自然規(guī)律，根據(jù)區(qū)域地形的水流方向數(shù)字矩陣計(jì)算每點(diǎn)處所流過的水量數(shù)值，沒有另外柵格水量注入的，其流水累積量為0；有另外一個(gè)柵格的流水注入的其流水累積量為1；有另外2個(gè)柵格的流水注入的其流水累積量為2；以此類推……，便可得到該區(qū)域水流累積數(shù)字矩陣，如圖415便是由一個(gè)簡單的水流方向矩陣產(chǎn)生的區(qū)域水流累積數(shù)字矩陣，為了方便，我們也把河網(wǎng)示意圖放在中間。 圖414 燕山水庫流域及子流域 流水累積量在流域的虛擬河網(wǎng)中，柵格水流流入到一個(gè)柵格的柵格個(gè)數(shù)的多少，就叫這一個(gè)柵格的流水累積量。一旦得到子流域號矩陣，就可以得出每個(gè)子流域的范圍，然后找出所以子流域的邊界點(diǎn)，于是就可以得到流域圖。如果不是交點(diǎn)，則該點(diǎn)賦值為當(dāng)前的流域號；如果是交點(diǎn)，則判斷流經(jīng)該點(diǎn)的鄰近的河流柵格的集水面積值是否大于閾值，如果其中任意一個(gè)大于閾值，則停止該條流路的追索，否則，繼續(xù)追索下去，直到當(dāng)前子流域不能再追索下去，然后增加子流域號。如圖413就是一個(gè)簡單的46水流方向矩陣形成的一個(gè)小流域，紅色箭頭表示水流方向流向流域外。圖412 燕山水庫流域河網(wǎng)圖 流域及子流域的劃分流域是指一個(gè)面積范圍，在這個(gè)范圍內(nèi)的積水都流向同一個(gè)出口。 圖411 由柵格水流方向到河網(wǎng)形示意圖本文運(yùn)用D8法提取的燕山水庫流域的虛擬河網(wǎng)如圖(412)。但是，這種方法的缺點(diǎn)有兩個(gè)：①會漏掉面積較大的山谷。識別谷點(diǎn)法就是通過比較網(wǎng)格的高程與相鄰網(wǎng)格高程以識別谷點(diǎn)。 圖49 多流向法示意圖如果中心柵格出流水量為，且坡度為、 、的柵格的高程比中心柵格低，則中心柵格出流的水量分配如下：本文采用D8法提取的整個(gè)流域內(nèi)的柵格的流向如(圖410)：圖410 燕山水庫流域柵格流向 河網(wǎng)的提取提取河網(wǎng)方法一般有兩種：谷點(diǎn)提取法[30]和網(wǎng)格流向法[31]。如果相鄰柵格的水流方向值為l，2，4，…，128，且此方向沒有指向中心柵格，則以此柵格的方向值作為中心柵格的方向值；，直至沒有任何柵格能被賦以方向值；對方向值不為1，2，4，…，128的柵格賦以負(fù)值(這種情況在經(jīng)洼地填充處理的DEM中不會出現(xiàn))。在極端情況下，如果8個(gè)相鄰柵格的高程值都與中心柵格高程值相同，則中心柵格方向值賦以255。這里假定計(jì)算區(qū)域是另一更大數(shù)據(jù)區(qū)域的一部分；，計(jì)算其對8個(gè)鄰域柵格的距離權(quán)落差值；，執(zhí)行以下步驟：①如果最大落差值小于0，則賦以負(fù)值以表明此柵格方向未定(這種情況在經(jīng)洼地填充處理的DEM中不會出現(xiàn))；②如果最大落差值大于或等于0，且最大值只有一個(gè)，則將對應(yīng)此最大值的方向值作為中心柵格處的方向值； ③如果最大落差值大于0，且有一個(gè)以上的最大值，則在邏輯上以查表方式確定水流方向。距離權(quán)落差是指中心柵格與相鄰柵格的高程差除以兩柵格間的距離，柵格間的距離與方向有關(guān)，如果鄰域柵格對中心柵格的方向值為3128，則柵格間的距離為2的開平方根，否則距離為1。輸出的方向值以2的冪值指定是因?yàn)榇嬖跂鸥袼鞣较虿荒艽_定的情況，此時(shí)須將數(shù)個(gè)方向值相加，這樣在后續(xù)處理中從相加結(jié)果便可以確定相加時(shí)中心柵格的鄰域柵格狀況。該方法認(rèn)為網(wǎng)格的產(chǎn)流點(diǎn)是網(wǎng)格中心點(diǎn),河道用一維的線描述,水流方向概括為8種可能流向,故稱之為D8法。目前，應(yīng)用最廣泛的單流向法是D8法。圖45是燕山水庫流域預(yù)處理后的柵格圖： 圖45 預(yù)處理后燕山水庫流域柵格圖 流向判定流向判斷是建立在33的DEM網(wǎng)格上,其判斷方法有單流向法和多流向法之分。，將帶有小平地標(biāo)記的柵格點(diǎn)的高程值增加一個(gè)小的高程增量(如DEM垂直分辨率的十分之一)，返回到步驟1，重復(fù)步驟1；。對于凹地，把洼地集水區(qū)域內(nèi)所有低于出流點(diǎn)的柵格高程升高至出流點(diǎn)高程，這樣凹地就成為一個(gè)平地，將平地內(nèi)的所有柵格點(diǎn)標(biāo)記。如果沒有潛在的出流點(diǎn)，或者存在任何洼地集水區(qū)域的邊界柵格，它的高程低于最低的潛在出流點(diǎn)，那么擴(kuò)大窗口，重復(fù)步驟2；。從洼地集水區(qū)域中找出那些潛在的出流點(diǎn)。 洼地標(biāo)定及抬升算法洼地柵格指相鄰8個(gè)柵格高程都不低于本柵格高程的柵格。在ArcView內(nèi)直接用菜單[Surface]下的命令[Derive Slope]，本文提取的燕山水庫流域Slope圖 (如圖44) 圖44 燕山水庫區(qū)域柵格坡度圖 DEM的預(yù)處理DEM預(yù)處理[26]是為了將DEM中的洼地和小平地改造成斜坡，也就是使DEM數(shù)據(jù)反映的地形特征均由斜坡構(gòu)成，這樣處理的目的就是使水流能夠根據(jù)地表徑流、漫流模型流出流域邊界：同時(shí)也為了流域信息提取的后續(xù)步驟的正確進(jìn)行，必須預(yù)處理洼地和小平地。三階差分計(jì)算方法：不同的常見的商用GIS軟件中采用的不同算法。 圖42 坡度和坡向圖 圖43 33窗口網(wǎng)格圖坡度的計(jì)算如下：Slope = 上式的計(jì)算通常有下面幾種算法：采用二階差分計(jì)算方法： 式中，為該網(wǎng)格中心的高程值，為網(wǎng)格間距（下同）。經(jīng)研究和試驗(yàn)發(fā)現(xiàn)擬合曲面法是解求坡度的最佳方法，擬合曲面法一般采用二次曲面，即如圖33的窗口，每個(gè)窗口中心為一個(gè)高程點(diǎn)，窗口在DEM高程矩陣中連續(xù)移動后，完成整幅圖的計(jì)算。 圖41 500m*500m區(qū)域柵格高程圖 地形坡度分析與計(jì)算坡度的計(jì)算方法有很多種，其中主要的計(jì)算地面平均坡度的方法共有三種，即四塊法、空間矢量分析法和擬合平面法[25]。本文采用網(wǎng)格的化分方法，將燕山水庫流域化分為500m*500m的柵格，提取的柵格高程如圖(41)，由圖上數(shù)據(jù)可以看出，,柵格柵格高程的變化區(qū)間不大。當(dāng)然，子流域還可以根據(jù)需要進(jìn)行第二級的劃分。首先，根據(jù)DEM進(jìn)行河網(wǎng)和子流域的提取，然后，基于等流時(shí)線的概念，將子流域分為若干條匯流網(wǎng)帶，在每一個(gè)匯流網(wǎng)帶上，圍繞河道劃分出若干個(gè)矩形坡面，在每個(gè)矩形坡面上，根據(jù)山坡水文學(xué)原理建立單元水文模型，進(jìn)行坡面產(chǎn)匯流計(jì)算,最后，進(jìn)行河網(wǎng)匯流演算[24]?；跂鸥竦膭澐址椒ㄒ曆芯繀^(qū)的不同,又可細(xì)分為兩類:一類是對于較小的實(shí)驗(yàn)場或小流域(幾百平方公里以內(nèi)),直接用DEM網(wǎng)格劃分,每個(gè)網(wǎng)格的大小一般為30m30m或50m50m等,該類方法在一些小尺度分布式水文物理模型(如SHE模型等)中比較流行；另一類是針對幾十萬到幾百萬平方公里的大流域,如一些大尺度的分布式水文模型,通常將研究區(qū)劃分為1km1km或更大的網(wǎng)格，每個(gè)網(wǎng)格單元根據(jù)DEM分辨率和模型精度要求,又可分為更小的網(wǎng)格。目前,流域離散的方法主要有三種[23]：單元網(wǎng)格(Grid)、山坡(Hillslop)和自然子流域(Subwatershed)。因?yàn)榈匦卧瓟?shù)據(jù)里有等高線和高程