金 鑫，宋 穎，潘斌林

（桂林理工大學(xué) 環(huán)境科學(xué)與工程學(xué)院，廣西 桂林541004）

傳統(tǒng)的分布式水文模型在進(jìn)行匯流計算時，將計算單元沿水流方向的寬度作為計算采用的匯流寬度。對于坡面匯流，該做法將坡面流從始至終作為片狀薄層水流處理，而在實際中，坡面片狀水流只是存在于匯流開始的階段，由于水流的沖刷和侵蝕作用，片狀水流很快就匯集成細(xì)溝流。根據(jù)研究，坡面細(xì)溝形成后，在細(xì)溝間距不變、雨強(qiáng)相同和坡度相等條件下，細(xì)溝不斷增寬，溝內(nèi)斷面平均流速會隨之減小，細(xì)溝寬度發(fā)展到一定程度后，細(xì)溝流流速與細(xì)溝泥沙起動流速相等，從而形成穩(wěn)定的細(xì)溝寬度［1］。在通常的匯流模擬中，由于將以細(xì)溝流為主的坡面水流作為片狀水流對待，坡面匯流計算的匯流寬度遠(yuǎn)遠(yuǎn)大于實際細(xì)溝流寬度的總和。在進(jìn)行河道匯流計算時，傳統(tǒng)的做法仍是將數(shù)字河道單元的寬度直接作為河道水流寬度對待，在進(jìn)行單純的產(chǎn)匯流模擬研究時，通過對匯流過程的一定技術(shù)處理，可以得到滿足精度要求的水流過程，但據(jù)此得到的水流流速與實際流速有很大不同，模型無法進(jìn)行較為準(zhǔn)確的分布式輸出，更談不上分布式驗證。而且在進(jìn)行分布式土壤侵蝕模擬或其它與流域產(chǎn)匯流相關(guān)的研究時，由于計算的水流流速與實際情況差別較大，極大地影響相關(guān)模擬計算結(jié)果。為此，根據(jù)作者在黃土地區(qū)分布式土壤侵蝕模型研究中的體會，提出了考慮流域匯流寬度隨集水面積變化的產(chǎn)匯流模型。

1 產(chǎn)匯流模型

1.1 流域概況及流域特征提取

研究區(qū)為位于伊洛河的東灣流域，該流域為以超滲產(chǎn)流為主要產(chǎn)流方式的黃土丘陵區(qū)。流域除出口水文站東灣站外，其上游還有欒川水文站，集水面積分別為2 657和239 k m2。流域特征提取以100 m柵格型DEM為數(shù)據(jù)源，流向判斷采用D 8法［2］，洼地采用傳統(tǒng)的填平法處理［3］。經(jīng)過洼地處理之后，研究區(qū)域內(nèi)所有網(wǎng)格的水流都可以經(jīng)出口流出，形成連續(xù)的匯流網(wǎng)絡(luò)。按照通常的做法，匯流網(wǎng)絡(luò)可以確定研究區(qū)域內(nèi)各柵格的上游累積網(wǎng)格數(shù)，當(dāng)某柵格的上游累積網(wǎng)格數(shù)達(dá)到一定的集水面積閾值時，就認(rèn)為該柵格為河道柵格。依次對研究區(qū)域內(nèi)的所有柵格進(jìn)行判斷，就可以確定河道柵格，提取出數(shù)字水系。

1.2 超滲產(chǎn)流計算

產(chǎn)生徑流的有效降雨等于實際降雨扣除降雨損失的部分。有效降雨由流域上的實際降雨以及蒸發(fā)蒸騰、林冠截留、填洼、土壤入滲等降雨損失而決定。由于研究區(qū)的降雨損失主要為林冠截留和土壤入滲，蒸發(fā)蒸騰和填洼損失量很小，可以忽略。因此，在產(chǎn)流的計算中，只考慮林冠截留和土壤入滲的降雨損失。

1.2.1 林冠截留計算 Dickinsoll［4］認(rèn)為覆被最大截留量等于葉面積指數(shù)（LAI）乘以一個特定的存儲值。降雨過程中冠層截留量主要和覆被的種類和生長時段有關(guān)，這一因素常常通過葉面積指數(shù)LAI來表達(dá)：

式中：Iv——冠層截留能力（mm）；Kc——冠層截留系數(shù)，與覆被類型有關(guān)；dc——覆被覆蓋度，反映了覆被空間分布情況；LAI——葉面指數(shù)；t——計算時段。下同。

1.2.2 產(chǎn)流計算 超滲產(chǎn)流模型是應(yīng)用下滲公式計算土壤下滲率f，然后依據(jù)霍頓下滲公式進(jìn)行產(chǎn)流計算，計算方法為：

式中：f——下滲能力（mm/min）；f0——最大下滲能力 （mm/min）；fc——穩(wěn) 定 下 滲 率 （mm/min）；k——下滲系數(shù)；I——雨強(qiáng)（mm/min）；R——徑流（mm）；ΔW——表示單位時間內(nèi)土壤含水量的增量（mm）。

因為下滲能力是土壤含水率的函數(shù)，土壤含水量的增量與下滲能力有如下關(guān)系：

根據(jù)霍頓公式，可得：

將（6）式帶入（5）式，得：

在實際計算中利用土壤含水量求下滲能力更為方便，故將霍頓下滲公式（2）與式（6）聯(lián)立，消除時間項，得到下滲能力與土壤含水率的函數(shù)關(guān)系：

可根據(jù)下滲能力與土壤含水量的關(guān)系通過積分求得實際下滲量，進(jìn)而求得產(chǎn)流量和時段末土壤含水量。

1.3 匯流計算

根據(jù)研究的對象和需要的不同，匯流過程可以用一維或二維的圣維南方程組來描述。對于平原區(qū)的匯流由于坡度較小，需要用二維模型來描述和模擬［5］；而對于山區(qū)性流域，一維圣維南方程組就能夠滿足模擬精度的需要。本文的研究區(qū)屬山區(qū)性流域，其匯流過程可用一維圣維南方程組來描述。

對于寬度為B，長度為L，凈雨強(qiáng)為Ie（x，t）的矩形坡面，用運動波來描述洪水波，用曼寧公式計算流速，圣維南方程組中的連續(xù)方程可寫為如下形式：

由于坡面上的水流屬寬淺水流，可用水深h近似替代水力半徑R，上式就成為關(guān)于水深的偏微分方程，利用Preissmann差分格式［6］對上式進(jìn)行離散，得：

其中：

根據(jù)初始及邊界條件，可得到一組關(guān)于水深的一元非線性方程，利用二分法及迭代計算可依次求得不同時刻各空間節(jié)點上的水深，由曼寧公式求得各節(jié)點上的流量、流速，得到該坡面的出流過程［7］。

對于有旁側(cè)入流的河道單元，圣維南方程組中的連續(xù)方程可寫為如下形式：

由于可根據(jù)河道形狀及匯流寬度將水力半徑表示為水深的函數(shù)，因此，上式也是關(guān)于水深的偏微分方程，對上式進(jìn)行離散有：

式中：a——水力半徑用水深表示后，連續(xù)方程中水深的指數(shù)。

河道匯流采用與坡面匯流同樣的方法得到一組一元非線性方程，最后得到河道單元的出流過程。

2 柵格單元匯流寬度的確定

在分布式水文模型的匯流計算中，常把實際水面面積等同于單元面積，即用單元尺寸代替實際的匯流寬度。這種假定只適用于自然條件下的2種特殊情況：一是坡面頂部產(chǎn)流時會有短暫的“片流”存在；二是流域中下游主河槽寬度剛好等于單元尺寸。

流速的大小不僅取決于河道的坡度，在流域的中下游更多地受到水位高差的影響。在分布式匯流計算中，一般只計算坡度對流速的影響，再加上DEM資料均化的影響，中下游計算出的流速往往偏小。為了擬合出口處的徑流過程線必然要人為增大流速計算系數(shù)，從而導(dǎo)致流域上游的流量過程線集中且提前。

分布式匯流的主要目的是尋找坡面流速與坡度、植被等下墊面特性之間的關(guān)系，更細(xì)致的還涉及到水量的累積，即從坡頂?shù)狡履_，水深與流量不斷增加，但水深的增量不與流量的增量成正比，因為流速也同時增加了，有研究表明坡面上水深的增量與流量增量的3/2次方成正比［8］。

在忽略次網(wǎng)格匯流的情況下，單元上的過境水流從坡頂?shù)钠髦饾u匯聚成股，進(jìn)入細(xì)溝、淺溝、切溝、小河直到大河，水流的集中程度逐漸增加。如果用單元內(nèi)所用過水通道的寬度之和除以單元尺寸作為縱坐標(biāo)，用集水面積或集水網(wǎng)格數(shù)量作為橫坐標(biāo)，就會得到一條曲線（見圖1），即：在坡頂片流階段過水寬度占滿全部網(wǎng)格，對應(yīng)曲線的“A”點；隨著水流迅速匯聚成股過水寬度比例迅速減少，對應(yīng)曲線的“B”點；當(dāng)水流進(jìn)入小河、大河時水流的集中程度仍然在增加，但由于集水面積的增加，水量也隨著增加，需要更大的過水通道，在曲線中表現(xiàn)為從“C”點拐向“D”點；隨著水量的進(jìn)一步增加，河道寬度也隨之增加，直到再次占滿整個單元，即曲線中的“E”點；當(dāng)集水面積很大時，相鄰的幾個網(wǎng)格都屬于河道，此時對應(yīng)著曲線中的“F”點。這條曲線可稱為“匯流寬度曲線”。

圖1 設(shè)想的柵格匯流寬度隨集水面積變化的曲線

2.1 匯流寬度曲線的建立

對于匯流寬度曲線，目前可以確定的是，該曲線一定存在一個等同于坡面單元寬度的起點A，一個水面寬度最小的轉(zhuǎn)折點C，以及流域出口點F，但要確定A—C與C—F點之間曲線的確切走勢，尚需一定的實測資料支持。但該曲線在A—C段逐漸下降、在C—F逐漸上升的趨勢是存在的，也是理論上可以解釋的。在A—C和C—F這2段的確切走勢無法確定前，可暫時概化為直線。

對于只有一個出口站的流域，根據(jù)多年不同流量下實測的河寬資料，可以確定出過水通道寬度曲線的終點F，起點A的寬度即為坡面柵格的寬度；對于水面寬度最小的轉(zhuǎn)折點C，其橫坐標(biāo)可以將其定為河道的閾值，而縱坐標(biāo)就需要通過率定來確定。匯流寬度曲線對于提高匯流模擬以及產(chǎn)輸沙模擬的精度是有實際意義的。

根據(jù)實測河寬資料得到東灣、欒川站的代表河寬分別為65和35 m。結(jié)合兩站的集水面積2 684和338 k m2，就得到了過水通道寬度曲線中的2個點；研究區(qū)的坡面－河道閾值為10 k m2，這樣得到一個東灣流域簡化的匯流寬度曲線（該曲線反映的是進(jìn)行匯流計算時所采用的水流寬度，并不是真正的水流寬度，見圖2）。

圖2 東灣站（含欒川）的匯流寬度簡化曲線

得到圖2的匯流寬度曲線之后，分布式匯流計算中可以根據(jù)單元水量計算出相應(yīng)的水深，根據(jù)曼寧公式計算出流速。相對于以前不考慮匯流寬度沿程變化的方法，該方法在理論上更合理一些。

2.2 匯流寬度曲線的意義

（1）統(tǒng)一坡面匯流參數(shù)與河道匯流參數(shù)。為了描述水流集中、水量累積等因素對匯流的影響，在分布式匯流中常通過兩套不同的參數(shù)對坡面匯流與河道匯流分別計算，二者通過集水面積閾值來區(qū)分。水流從坡面流到河道水流是一個連續(xù)的過程，用匯流寬度曲線來描述更加合理，而且不再需要兩套不同的參數(shù)。

（2）便于確定合適的DEM分辨率。當(dāng)區(qū)域的集水面積在“E”點左側(cè)時，河道匯流與坡面匯流的計算方法都符合流域特征提取時最陡坡度法中唯一出口的假定；當(dāng)集水面積在“E”點右側(cè)時，就會同時有若干相鄰網(wǎng)格有水流通過，這需要同時改變流域特征提取方法和河道匯流計算方法。因此，可以把“E”點作為確定DEM分辨率的指標(biāo)：只要能夠保證計算區(qū)域出口處的集水面積在“E”點左側(cè)，就可以取更小的分辨率。

（3）計算水位、矯正流速后計算輸沙。分布式水文模型的水流流速并不是計算真正的流速，它包含了影響出口處匯流結(jié)果的所有誤差；而在輸沙計算中，流速誤差將以幾何級數(shù)傳遞到攜沙能力的誤差中［9］。匯流寬度曲線的使用，使得不同網(wǎng)格的流速計算相對于原來更趨合理，有理由進(jìn)一步降低輸沙計算的誤差。

（4）劃分坡面單元細(xì)溝和細(xì)溝間的比例、溝道單元溝道和坡面的比例，避免在單元柵格內(nèi)只進(jìn)行單一方式侵蝕計算帶來的誤差，提高侵蝕產(chǎn)沙計算的精度。

（5）匯流特征寬度及其空間分布的引入，除了細(xì)化水量－流速關(guān)系外，還為坡面流再入滲的描述、直接產(chǎn)流面積（水面面積）時空變化的描述提供了可能。

3 模型應(yīng)用

選擇東灣和欒川2站有同步資料的場次洪水19810714和19831003，使用本文所建立的產(chǎn)匯流模型，對東灣站的洪水過程進(jìn)行模擬，把東灣站的洪水過程線擬合好之后，在對應(yīng)欒川站位置的單元上可以同時得到欒川站的過程線，并用同一套參數(shù)直接模擬欒川站的洪水過程，模擬結(jié)果見表1。另外，對上述2場洪水，按照傳統(tǒng)方法，不考慮匯流寬度隨集水面積的變化，直接用柵格寬度代替匯流寬度進(jìn)行匯流計算，兩站的模擬結(jié)果見表2。

表1 考慮匯流寬度變化的模擬結(jié)果

表2 不考慮匯流寬度變化的模擬結(jié)果

從東灣站的模擬結(jié)果來看，兩種方法的模擬精度大致相當(dāng)，而對欒川站的模擬，兩種方法所模擬的精度有較大差別。從欒川站的模擬結(jié)果可以看出，考慮匯流寬度隨集水面積變化的模擬方法，其精度雖不十分令人滿意，但比傳統(tǒng)方法的精度有較大提高，且計算的洪水過程與實測更為接近。考慮匯流寬度變化直接進(jìn)行欒川站的洪水模擬過程見圖3—4，不考慮匯流寬度變化直接進(jìn)行欒川站的洪水模擬過程見圖5—6。

圖3 考慮匯流寬度變化的欒川站19810714洪水過程

圖5 不考慮匯流寬度變化的欒川站19810714洪水過程

圖4 考慮匯流寬度變化的欒川站19831003洪水過程

圖6 不考慮匯流寬度變化的欒川站19831003洪水過程

4 結(jié)論

匯流寬度隨集水面積變化的分布式產(chǎn)匯流模型與傳統(tǒng)分布水水文模型的不同，在于使用了匯流寬度曲線，通過該曲線對不同集水面積的柵格匯流寬度進(jìn)行確定，改變了傳統(tǒng)方法將柵格寬度直接作為匯流寬度的做法。雖然匯流寬度曲線變化的確切趨勢尚需大量的實測資料進(jìn)行驗證，但從東灣流域的應(yīng)用情況來看，匯流寬度曲線的引入，使坡面匯流和河道匯流可以使用統(tǒng)一的匯流參數(shù)，并且在一定程度上提高了分布式輸出的精度。

［1] 楊具瑞，曹叔尤，劉興年，等.黃土坡面細(xì)溝侵蝕穩(wěn)定寬度的動力學(xué)研究［J］.昆明理工大學(xué)學(xué)報：理工版，2004，29（4）：159－163.

［2] O＇Callaghan J F，Mark D M.The extraction of drainage net works from digital elevation data.［J］.Computer Vision，Graphics，and Image Pr ocessing，1984，28：323－344.

［3] Jenson S K，Do mingue J O.Extraction topographic str ucture fro m digital elevation data f or geographic information system analysis［J］.Photogrammetric Engineering and Remote Sensing.1988，54（11）：1593－1600.

［4] 王加虎.分布式水文模型原理與方法研究［D］.南京：河海大學(xué)，2006.

［5] Giammarco P D，Todini E，Lamberti P.A conservative finite elements approach to overland flow：The control volu me finite element for mulation［J］.Journal of Hydrology，1996，175：267－291.

［6] 張建云.水文學(xué)手冊［M］.李紀(jì)生，譯.北京：科學(xué)出版社，2002.

［7] 李麗.分布式水文模型的匯流演算研究［D］.南京：河海大學(xué)，2006.

［8] Wood E F，Hebson C S.On Hydrologic si milarity 1：Derivation of t he di mensionless flood frequency cur ve［J］.Water Resour Res.，1986，22（11）：1549－1554.

［9] 中國水利學(xué)會泥沙專業(yè)委員會.泥沙手冊［M］.北京：中國環(huán)境科學(xué)出版社，1992.