【正文】 w to reservoir in the two floodwaters, isochrones process of the watershed is simulated, the results fit well.Key Words: DEM。摘 要隨著現(xiàn)代社會(huì)數(shù)字化進(jìn)程的加快，基于數(shù)字高程模型(Digital Elevation Model, DEM)的分布式水文模型(Distributed Hydrology Model)在水利行業(yè)的主導(dǎo)地位也日益明顯，它避開(kāi)了傳統(tǒng)水文學(xué)依靠多年的水文觀測(cè)資料的計(jì)算方法，而依據(jù)流域的相關(guān)參數(shù)和地理信息進(jìn)行水文計(jì)算，給水利行業(yè)注入了新的活力。 Distributed Hydrology Model。分布式水文模型是在20世紀(jì)70年代提出的，經(jīng)歷了一個(gè)循序漸進(jìn)的發(fā)展歷程：先是在一個(gè)小場(chǎng)地實(shí)地測(cè)量每個(gè)塊的高程、坡度，后來(lái)隨著計(jì)算機(jī)的應(yīng)用和計(jì)算機(jī)圖形處理的發(fā)展，特別是在出現(xiàn)GIS后，DEM成為提取地形參數(shù)的有效平臺(tái)，在這個(gè)方面許多學(xué)者做出大量研究，目前，在柵格高程、柵格坡流、柵格流向方面取得一致，在虛擬河網(wǎng)的提取方面目前有多種方法，哪種方法更行之有效，還沒(méi)有定論，但是，分布式水文模型建立方法還是被學(xué)者們認(rèn)可。概念性模型(conceptual model)和物理模型(physicallybased model)；從反映水流運(yùn)動(dòng)空間變化能力的角度出發(fā)，水文模型又可以分為集總式模型(lumped model)和分布式模型(distributed model)，還有一種水文模型介于集總式模型和分布式模型之間，稱為“半分布式”水文模型[2]。 國(guó)內(nèi)研究現(xiàn)狀國(guó)內(nèi)的分布式水文模型的研究起步較早，20世紀(jì)60年代,趙人俊等[9]提出旨在考慮包氣帶持水能力的空間分布對(duì)蓄滿產(chǎn)流影響的新安江模型就是一個(gè)經(jīng)典的半分布式水文模型。本文第一章是對(duì)國(guó)內(nèi)外和本人對(duì)所研究課題的探討；第二章是對(duì)流域的下墊信息、水文參數(shù)提??；第三章是在流域地形圖基礎(chǔ)上生成流域的DEM(數(shù)字高程模型)；第四章是由流域DEM提取分辨率為500m500m柵格的高程、坡度、河網(wǎng)，流量積累等地形參數(shù)；第五章是建立柵格單元水文模型，第六章對(duì)地表徑流模擬計(jì)算；最后一章是對(duì)模型的驗(yàn)證。燕山水庫(kù)大壩屬土石混合壩型，等別和規(guī)模屬二等大(Ⅱ)型，地震設(shè)防烈度為6度。 地形和土壤燕山水庫(kù)流域海拔除幾個(gè)山頂外基本在80m～250m之間，落差較小、地勢(shì)相對(duì)平坦、坡度變化不大，流域的北、西南、東南均有丘陵，因此，燕山水庫(kù)流域處于一個(gè)狹長(zhǎng)斷陷盆地之中，也形成一個(gè)天然的閉合流域。影響曼寧系數(shù)的因素很多，其中包括河網(wǎng)形狀、地形、植被覆蓋、地表侵蝕特征等，Engman(1986)給出了適于分布式水文模型的坡面徑流和河網(wǎng)匯流演算的有效曼寧系數(shù)[15]，對(duì)于本流域的土地利用實(shí)際情況，選取的曼寧系數(shù)初值如下：表(23) 曼寧系數(shù)取值 覆蓋類型曼寧系數(shù)初始取值曼寧系數(shù)合理范圍農(nóng)田森林、灌木林荒草地作物耕地城鎮(zhèn)用地第三章 燕山水庫(kù)流域數(shù)字高程模型的構(gòu)建 數(shù)字高程模型概念及表示方法 數(shù)字高程模型的概念及應(yīng)用數(shù)字高程模型DEM(Digital Elevation Model),其目的是用攝影測(cè)量或其他技術(shù)手段獲得地形數(shù)據(jù)給地理空間定位，它是利用坐標(biāo)場(chǎng)中大量的已知（X，Y，Z）坐標(biāo)點(diǎn)及屬性值對(duì)連續(xù)地面的數(shù)字表示，是一個(gè)描述地球表面形態(tài)多種信息空間分布的有序數(shù)值陣列，DEM的實(shí)質(zhì)就是對(duì)地球表面地形地貌的一種離散的數(shù)學(xué)表達(dá)[16]。圖31 DEM表示方法 DEM的三種主要模型數(shù)字高程模型構(gòu)建有以下幾種方法：；；； ：通過(guò)遙感立體像對(duì)，根據(jù)視差模型自動(dòng)選配左右影像，用攝影測(cè)量的方法建立空間高程立體模型，量取密集高程數(shù)據(jù)，建立數(shù)字高程模型；：根據(jù)有限個(gè)離散點(diǎn)的高程、采用多項(xiàng)式或樣條函數(shù)求得擬合公式（即曲面函數(shù)方程），從而可以逐一計(jì)算各點(diǎn)的高程，可得到擬合的柵格DEM。等高線模型指根據(jù)已知的高程值的集合來(lái)表示整個(gè)區(qū)域的高程狀況，因?yàn)槊恳粭l等高線對(duì)應(yīng)一個(gè)已知的高程值，所以一系列等高線集合和它們的高程值一起就構(gòu)成了一種地面高程模型。所以TIN是一個(gè)三維空間的分段線性模型，在整個(gè)區(qū)域內(nèi)連續(xù)但不可微。 數(shù)字高程模型的構(gòu)建方法 數(shù)字高程模型構(gòu)建有以下幾種方法：：在地形圖上蒙上一層透明的方格膜片，逐格讀取中心或角點(diǎn)的高程值，寫(xiě)成數(shù)字高程矩陣的形式，構(gòu)成數(shù)字高程模型。過(guò)去100多年來(lái)，攝影測(cè)量學(xué)經(jīng)歷了模擬法和解析法這兩個(gè)階段，現(xiàn)已進(jìn)人了它的第三個(gè)階段即數(shù)字?jǐn)z影測(cè)量階段[18]。思路為：首先，在數(shù)據(jù)點(diǎn)集中取其中任一點(diǎn)P1，在其余各點(diǎn)中尋找與此點(diǎn)距離最近的點(diǎn)P2，連接P1和P2構(gòu)成第一邊，然后在其余所有點(diǎn)中尋找與這條邊最近的點(diǎn)，找到后即構(gòu)成第一個(gè)三角形；再以這個(gè)三角形新生成的兩邊為底邊分別尋找距它們最近的點(diǎn)構(gòu)成第二個(gè)三角形、第三個(gè)三角形……，依此類推，直到把所有的點(diǎn)全部連入三角網(wǎng)中，這是最簡(jiǎn)單的生成TIN的方法。這種三角形稱作Delaunay三角形，它與Voronoi多邊形存在對(duì)偶關(guān)系，Delaunay三角形的外接圓圓心是與三角形相關(guān)的Voronoi多邊形的一個(gè)頂點(diǎn)。合并算法、逐點(diǎn)插入法及三角網(wǎng)生長(zhǎng)算法[20] [21]等?；跂鸥竦膭澐址椒ㄒ曆芯繀^(qū)的不同,又可細(xì)分為兩類:一類是對(duì)于較小的實(shí)驗(yàn)場(chǎng)或小流域(幾百平方公里以內(nèi)),直接用DEM網(wǎng)格劃分,每個(gè)網(wǎng)格的大小一般為30m30m或50m50m等,該類方法在一些小尺度分布式水文物理模型(如SHE模型等)中比較流行；另一類是針對(duì)幾十萬(wàn)到幾百萬(wàn)平方公里的大流域,如一些大尺度的分布式水文模型,通常將研究區(qū)劃分為1km1km或更大的網(wǎng)格，每個(gè)網(wǎng)格單元根據(jù)DEM分辨率和模型精度要求,又可分為更小的網(wǎng)格。 圖41 500m*500m區(qū)域柵格高程圖 地形坡度分析與計(jì)算坡度的計(jì)算方法有很多種，其中主要的計(jì)算地面平均坡度的方法共有三種，即四塊法、空間矢量分析法和擬合平面法[25]。在ArcView內(nèi)直接用菜單[Surface]下的命令[Derive Slope]，本文提取的燕山水庫(kù)流域Slope圖 (如圖44) 圖44 燕山水庫(kù)區(qū)域柵格坡度圖 DEM的預(yù)處理DEM預(yù)處理[26]是為了將DEM中的洼地和小平地改造成斜坡，也就是使DEM數(shù)據(jù)反映的地形特征均由斜坡構(gòu)成，這樣處理的目的就是使水流能夠根據(jù)地表徑流、漫流模型流出流域邊界：同時(shí)也為了流域信息提取的后續(xù)步驟的正確進(jìn)行，必須預(yù)處理洼地和小平地。對(duì)于凹地，把洼地集水區(qū)域內(nèi)所有低于出流點(diǎn)的柵格高程升高至出流點(diǎn)高程，這樣凹地就成為一個(gè)平地，將平地內(nèi)的所有柵格點(diǎn)標(biāo)記。該方法認(rèn)為網(wǎng)格的產(chǎn)流點(diǎn)是網(wǎng)格中心點(diǎn),河道用一維的線描述,水流方向概括為8種可能流向,故稱之為D8法。在極端情況下，如果8個(gè)相鄰柵格的高程值都與中心柵格高程值相同，則中心柵格方向值賦以255。但是，這種方法的缺點(diǎn)有兩個(gè)：①會(huì)漏掉面積較大的山谷。如果不是交點(diǎn)，則該點(diǎn)賦值為當(dāng)前的流域號(hào)；如果是交點(diǎn)，則判斷流經(jīng)該點(diǎn)的鄰近的河流柵格的集水面積值是否大于閾值，如果其中任意一個(gè)大于閾值，則停止該條流路的追索，否則，繼續(xù)追索下去，直到當(dāng)前子流域不能再追索下去，然后增加子流域號(hào)。權(quán)值矩陣是一個(gè)連續(xù)的數(shù)字矩陣，它表示了在一次暴雨中的平均降雨量，可用來(lái)計(jì)算一個(gè)流域內(nèi)的流走的降雨量的多少,并且假定所有的降雨量沒(méi)有被地表水截留、蒸發(fā)或損失，輸出的流水累積量矩陣表示了流經(jīng)每個(gè)柵格的降雨量。本文的庫(kù)容值是用用GIS軟件ArcView計(jì)算的，其原理是就是求解DEM體積，即在TIN的基礎(chǔ)上求解下表面為T(mén)IN上表面為水平面組成的空間體的體積，通過(guò)計(jì)算三棱柱體積進(jìn)行累加得到，每一個(gè)三棱柱體積計(jì)算公式為： ；；式中zi，xi，yi為T(mén)IN中一個(gè)三角形三個(gè)頂點(diǎn)的坐標(biāo)值；H為水位高程；S是三棱柱底面積也即TIN中的一個(gè)三角形的面積。目前，建立這樣復(fù)雜的水文數(shù)學(xué)物理方程還比較困難，也難應(yīng)用到復(fù)雜的大流域。圖51 柵格地表、地下水輸出輸入 步長(zhǎng)內(nèi)柵格地表水量平衡方程 (51) (52) 步長(zhǎng)內(nèi)柵格地下水量平衡方程 (53)計(jì)算時(shí)間步長(zhǎng)；步長(zhǎng)內(nèi)總降雨量；步長(zhǎng)內(nèi)植被截流；步長(zhǎng)內(nèi)總蒸散發(fā)量；步長(zhǎng)內(nèi)鄰近柵格地表入流量； 步長(zhǎng)內(nèi)柵格地表出流量； 步長(zhǎng)內(nèi)下滲量；S步長(zhǎng)內(nèi)柵格地表蓄水量；步長(zhǎng)內(nèi)地下土層含水量；步長(zhǎng)內(nèi)到達(dá)地面的雨量；步長(zhǎng)內(nèi)鄰近柵格地下入流量；步長(zhǎng)內(nèi)柵格地下出流量； 步長(zhǎng)內(nèi)柵格地下蓄水量；步長(zhǎng)內(nèi)深層下滲量； 步長(zhǎng)內(nèi)柵格土壤蒸發(fā)量。此法是將相鄰雨量站用直線連接而成為若干個(gè)三角形，然后對(duì)三角形各邊作垂直平分線，連接這些垂直平分線的交點(diǎn)，得若干個(gè)多邊形，每個(gè)多邊形內(nèi)各有一個(gè)雨量站，該多邊形即是該雨量站的控制面積。其計(jì)算方法可用如下公式[36]:()= (59)式中是截留系數(shù)，它是降水量P的函數(shù)。常用的方法有如下：；(如能量平衡法,空氣動(dòng)力學(xué)法,能量平衡空氣動(dòng)力學(xué)法和渦度相關(guān)技術(shù)等)；，利用衛(wèi)星遙感并結(jié)合模式模擬研究非均勻陸面上的蒸散發(fā)是一個(gè)新的趨勢(shì),該項(xiàng)技術(shù)在我國(guó)剛剛起步。 (516) (HarGreaves 1994)為宇宙射輻射MJ/(m2day)；為日平均氣溫（日最大氣溫與最小氣溫的均值）；為日氣溫變化 (日最大氣溫與最小氣溫的差值)。圖54 植被截留量與降水量P的關(guān)系植被蓋度(%)：指植被(包括葉、莖、枝)在單位面積內(nèi)植被的垂直投影面積所占百分比。在無(wú)法取得降雨量空間分布函數(shù)的情況下，泰森多邊形法為許多水文計(jì)算選用，本文選擇泰森多邊法計(jì)算燕山水庫(kù)流域降雨分布（如圖53）。目前，獲取降雨資料的手段有兩種方法：雨量站觀測(cè)和雷達(dá)遙感法。在傳統(tǒng)的集總式水文模型中，它能夠大大地簡(jiǎn)化計(jì)算量，提高模型的適應(yīng)能力。庫(kù)容曲線是水庫(kù)水位、水量擬合曲線。圖416就是采用以上方法提取的燕山水庫(kù)流域的流水累積量矩陣，由柵格顏色深與淺表示流水積累量的多與少，由圖可以看出流水累積量矩陣的圖和河網(wǎng)一致，流域流水積累量最大值與流域的柵格總數(shù)一致。一旦得到子流域號(hào)矩陣，就可以得出每個(gè)子流域的范圍，然后找出所以子流域的邊界點(diǎn)，于是就可以得到流域圖。 圖411 由柵格水流方向到河網(wǎng)形示意圖本文運(yùn)用D8法提取的燕山水庫(kù)流域的虛擬河網(wǎng)如圖(412)。如果相鄰柵格的水流方向值為l，2，4，…，128，且此方向沒(méi)有指向中心柵格，則以此柵格的方向值作為中心柵格的方向值；，直至沒(méi)有任何柵格能被賦以方向值；對(duì)方向值不為1，2，4，…，128的柵格賦以負(fù)值(這種情況在經(jīng)洼地填充處理的DEM中不會(huì)出現(xiàn))。輸出的方向值以2的冪值指定是因?yàn)榇嬖跂鸥袼鞣较虿荒艽_定的情況，此時(shí)須將數(shù)個(gè)方向值相加，這樣在后續(xù)處理中從相加結(jié)果便可以確定相加時(shí)中心柵格的鄰域柵格狀況。，將帶有小平地標(biāo)記的柵格點(diǎn)的高程值增加一個(gè)小的高程增量(如DEM垂直分辨率的十分之一)，返回到步驟1，重復(fù)步驟1；。 洼地標(biāo)定及抬升算法洼地柵格指相鄰8個(gè)柵格高程都不低于本柵格高程的柵格。經(jīng)研究和試驗(yàn)發(fā)現(xiàn)擬合曲面法是解求坡度的最佳方法，擬合曲面法一般采用二次曲面，即如圖33的窗口，每個(gè)窗口中心為一個(gè)高程點(diǎn)，窗口在DEM高程矩陣中連續(xù)移動(dòng)后，完成整幅圖的計(jì)算。首先，根據(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]。對(duì)于凸殼的生成采用格雷厄姆算法，該算法是求解平面點(diǎn)集凸殼問(wèn)題的最佳算法，然后求出包含新插入點(diǎn)p的外接圓的三角形，并按LOP算法優(yōu)化三角網(wǎng)，迭代計(jì)算完成整個(gè)區(qū)域內(nèi)的三角網(wǎng)的建立。這樣，在生成初步三角剖分網(wǎng)的同時(shí)，也實(shí)現(xiàn)了圖形優(yōu)化的處理。常用Delaunay三角剖分方法。利用全數(shù)字?jǐn)z影測(cè)量系統(tǒng)進(jìn)行數(shù)據(jù)采集必須對(duì)用高精度掃描儀將航片或衛(wèi)片掃描成數(shù)字影像或直接獲取遙感數(shù)據(jù)進(jìn)行內(nèi)定向，從而確定數(shù)字影像的掃描坐標(biāo)系和影像坐標(biāo)系之間的相互轉(zhuǎn)換關(guān)系。：對(duì)有限個(gè)離散點(diǎn)或等高線，按一定優(yōu)化原則將每三個(gè)鄰近點(diǎn)聯(lián)接成三角形，每個(gè)三角形代表一個(gè)局部平面，再根據(jù)每個(gè)平面方程，可計(jì)算各柵格點(diǎn)高程，生成DEM。圖(34)TIN模型。如圖可知，等高線線模型中不同的等高線是離散的，只表達(dá)了區(qū)域的部分高程值，無(wú)論等高線之間距離多么小，都不可能表達(dá)整個(gè)地形。如圖32所示。由于用數(shù)字形式表達(dá)地形表面,DEM具有如下顯著特點(diǎn)：，產(chǎn)生多種比例尺的地形圖、縱橫斷面圖和立體圖；；。流域內(nèi)土壤土層較厚，約2米左右，根據(jù)河南省土壤分布圖、南陽(yáng)土壤分布圖以及中國(guó)土壤數(shù)據(jù)庫(kù), 流域內(nèi)表層土壤為典型的淋溶褐土系列，主要以黃褐土為主,其次為粗骨土，黃棕壤、石質(zhì)土、灰潮土，土壤分布平面圖見(jiàn)圖(22)。 流域水文氣象燕山水庫(kù)流域?qū)儆谂瘻貛Ъ撅L(fēng)氣候，氣候變化完全受季風(fēng)支配：在夏季季風(fēng)活躍，南方暖濕汽流不斷襲來(lái)，造成夏季潮濕多雨，氣溫日差大；極端最高氣溫可達(dá)400C；九月以后，季風(fēng)逐漸南移，北方冷空氣不斷南下，冷高壓氣團(tuán)控制，形成冬季寒冷，降雨量少。本文采用了河南省水利廳設(shè)計(jì)院提供1：5000的燕山水庫(kù)流域電子地形圖以及河南省水文局提供的相關(guān)水文地質(zhì)資料、1：10000河南省地形圖數(shù)據(jù)和因特網(wǎng)上提供的河南省土地利用圖，全國(guó)土壤分布圖以及相關(guān)參數(shù)建立時(shí)空關(guān)系，來(lái)反映地形、植被和土壤