孔祥意，劉志雨，，周國良，李致家

(1.水利部信息中心(水利部水文水資源監(jiān)測預報中心)，北京 100053； 2.河海大學水文水資源學院，江蘇 南京 210098)

流域內的水文過程具有四維性，即三維空間和時間屬性，具有物理基礎的分布式流域水文模型[1]正是為描述水文過程的四維性而提出的。該類模型的根本特點為模型的輸入(包括參數、初始條件、邊界條件等)和輸出均具有空間分布性[2-3]。目前，具有物理基礎的分布式流域水文模型大都以DEM柵格作為模型的計算單元[4-5]。該類模型認為，在每一個計算單元網格內的氣象、下墊面及地形特征分布均勻[6]。

計算單元網格越小，越能表征流域的地形地貌特征，保留更多的流域地形地貌信息，真實體現流域的水文響應特性。但在流域面積一定的情況下，采用小尺度網格離散流域會使計算單元網格數量急劇增加，隨之而來的是繁重的計算量，影響分布式流域水文模型在實際應用中參數的率定和實時計算的時效性。因此，在較大流域中，大都采用較低分辨率DEM(如1 000 m ×1 000 m)構建具有物理基礎的分布式流域水文模型[7-8]。

目前普遍認為，具有物理基礎的分布式流域水文模型的最顯著優(yōu)勢是：其中的參數可以根據下墊面條件及相關文獻資料(如公開出版發(fā)行的手冊等，下同)中提供的參數值與下墊面條件對照關系直接獲得，不需要利用實測水文資料率定，使得該類模型可以用于無資料流域。在有資料流域的應用結果表明，絕大多數情況下直接采用文獻資料中的參數，尤其是糙率(流速參數)，模型的模擬效果不好。要想得到理想的模擬結果，必須重新率定，且不同大小計算單元的模型率定后其糙率值差別明顯[9-10]。分析認為，文獻中列出的參數與下墊面條件的對應關系，是基于真實的地形特征獲得的；而分布式水文模型采用DEM柵格為計算單元，對流域地表進行了離散處理，每一個柵格內的地形特征(如坡度)被均化，已非真實地表特征。計算單元大小不同，對地表特征均化程度不同。因此，在具有物理基礎的分布式流域水文模型的實際應用過程中，必須對由采用柵格作為計算單元導致的均化作用進行矯正。

目前的矯正方法主要分為兩類：第一類為依據網格尺度直接矯正流速(或水流傳播時間)[11]，第二類為矯正地表糙率或矯正地形坡度[12]。目前，在有實測資料流域采用各種優(yōu)化方法由水文氣象資料率定模型參數，其實就是基于提取的非真實地形坡度對地表糙率進行矯正。王曉燕等[13]建立了流域平均坡度與DEM分辨率間的相關關系，并得到線性回歸方程，用以矯正地形坡度，但所用DEM柵格尺度小于100 m，研究結論是否適用于更低分辨率的DEM有待驗證。

本文所用模型為TOKASIDE(topgraphic kinematic approximation and saturation-infiltration double excess grid-based distributed model)[14]，是根據Todini 等[15-16]提出的利用非線性水庫方法模擬水文過程的理論構建的基于地形、運動波方法以及蓄滿超滲雙產流機制的網格化分布式模型。該模型為具有物理基礎的分布式流域水文模型，以DEM柵格為計算單元，采用非線性運動波方程模擬水流運動[17-18]。模型將每一個計算單元中的水文過程概化為4個結構上相似的非線性水庫方程，分別描述淺層土壤與深層土壤中的排水以及坡地地表徑流和河道徑流。模型采用了飽和及超滲產流機制，在整個降雨過程中，隨著土壤含水量與降雨強度的變化，超滲與蓄滿產流機制可能在每一個計算單元網格內交替發(fā)生[16]。

本文基于TOKASIDE模型，從固定參數下計算單元邊長對模擬結果的影響機理及不同計算單元邊長模型間參數移用誤差矯正方法兩個方面進行研究。

1 研究區(qū)概況

以淮河流域息縣水文站以上部分(簡稱淮河流域)為研究區(qū)，地形情況如圖1所示。所用DEM數據來源于地理空間數據云(http：//www.gscloud.cn)提供的ASTER GDEM數據，柵格分辨率為30 m×30 m，所用其他分辨率DEM均為由ARCGIS工具再采樣生成得到。由30 m×30 m分辨率DEM提取的流域面積為10 050 km2。為研究需要進行了子流域劃分，劃分結果如圖2所示，子流域1～8的面積分別為2 166 km2、1 934 km2、694 km2、604 km2、268 km2、791 km2、2 012 km2和1 579 km2。

圖1 淮河流域地形Fig.1 Topography of Huaihe River Basin

圖2 子流域劃分結果Fig.2 Sub watershed Division result

2 計算單元邊長對TOKASIDE模擬結果的影響

2.1 對模擬徑流過程的影響

由于研究流域面積較大，如果采用高分辨率DEM構建模型，則參數率定及洪水模擬的時間均會很長，因此采用基于2 000 m×2 000 m分辨率DEM的模型進行參數率定。另外又構建了基于300 m×300 m、500 m×500 m、800 m×800 m、1 000 m×1 000 m以及1 500 m×1 500 m分辨率DEM的5個模型，將率定的模型參數直接用于該5個模型，對201602號洪水的模擬結果如圖3所示。由圖3可以看出：(a)隨著計算單元邊長的減小，洪峰值逐漸增大，峰現時間逐漸提前；(b)洪峰值及峰現時間隨單元尺度呈非均勻變化，分別在300 m和800 m處出現突變；(c)在某一單元尺度范圍(如500～800 m，1 000～1 500 m)，模擬結果相近。

圖3 不同計算單元邊長模型的模擬結果Fig.3 Simulation results of different cell scale models

2.2 對水量平衡要素的影響

選擇模擬的流量過程退水段同時刻為結束點，基于不同計算單元邊長模型模擬出的201602號洪水徑流成分變化統計結果見表1。由表1可以看出：隨著計算單元邊長的增大，流域內的土壤含水量增量、地表水增量以及河道內水量(由于在模擬初期設定了15 m3/s的河道基流，表中負值代表洪水模擬結束后河道內總水量減少)呈增加趨勢，而流出流域出口的徑流量明顯減少。這些變化均是由于坡地及河道內坡度或比降的減小所致。

表1 201602號洪水徑流成分統計

3 計算單元邊長對模擬結果影響機理分析

3.1 對提取的流域坡度的影響機理

利用ARCGIS工具提取13個分辨率DEM對應的各子流域及全流域的平均坡度，如表2和圖4所示。

由表2和圖4可以看出：(a)山區(qū)性子流域，隨著網格尺度的增大，平均坡度減小明顯，且坡度減小的梯度逐漸減??；(b)處于平原內的子流域(如4號和 5號)，DEM 柵格大小對坡度提取結果影響很小。

圖4 流域坡度與DEM柵格尺度的關系Fig.4 Relationship between watershed slope and DEM grid scale

為了進一步研究流域坡度受 DEM 分辨率影響的變化規(guī)律，對表2進行歸一化處理，即以 30 m×30 m 分辨率 DEM 對應的坡度為基數，用其他分辨率 DEM 對應的坡度與該基數之比，表示DEM分辨率降低后形成的流域坡度是30 m×30 m分辨率DEM對應坡度(認為是流域的真實坡度)的多少倍，即歸一化處理后得到的結果為坡度變化與 DEM 柵格尺度之間的關系。歸一化處理后的結果見表3和圖5。

表2 由不同分辨率DEM提取的流域坡度

由表3和圖5可以看出，歸一化處理后不同子流域(含總流域)間坡度隨DEM柵格尺度的變化規(guī)律高度相似，可以取平均線作為共同的變化趨勢或變化規(guī)律，同時也說明不同流域間可以采用統一的表示方法，便于向無高分辨率DEM流域移用。從圖5可見，提取的淮河流域平均坡度隨DEM柵格尺度的增大呈非均勻減小時，在400 m和900 m處分別出現一次突變。

表3 不同分辨率DEM的流域坡度變化值

圖5 流域坡度比與DEM柵格尺度的關系Fig.5 Relationship between watershed slope ratio and DEM grid scale

3.2 對提取的河段特征的影響機理

圖6為由30 m×30 m分辨率DEM 提取的淮河流域水系(采用500 km2面積閾值)，圖7為4號、5號和8號3個河段河長度隨DEM柵格尺度的變化。由圖7可以看到河段長度整體呈現隨柵格尺度增大而減小的變化趨勢。分析認為，DEM分辨率對不同類型河段(內鏈和外鏈)的影響不同。對外鏈河段的影響主要表現在3個方面：(a)下游節(jié)點位置的變化；(b)由于分水線的變化，在同樣集水面積閾值情況下，上游節(jié)點位置的變化(如8號子流域在DEM柵格邊長大于800 m以后)；(c)河段彎曲程度的變化。對內鏈河段的影響主要表現在兩個方面：(a)上、下節(jié)點位置的變化；(b)河段彎曲程度的變化。

圖6 河段及編碼Fig.6 River reach and code

圖7 河段長度隨DEM柵格尺度的變化Fig.7 Variation of river reach scale with DEM grid scale

3.3 計算單元邊長的影響分析

根據對提取的流域特征受DEM分辨率影響的分析可以看出，由于DEM柵格尺度的增大對地面點的高程具有均化作用，造成絕大部分柵格的坡度減小，進而減小了流域的平均坡度。盡管各子流域平均坡度可能相差很大，但平均坡度隨著DEM柵格尺度的變化規(guī)律(相對變化)卻非常一致。隨著DEM柵格尺度的增大，流域坡度并非以均勻的梯度減小，最大減小梯度發(fā)生在200 m范圍內。

柵格尺度的增大，改變了柵格的高程，從而改變了相鄰網格間徑流輸入與輸出關系，影響到每一個柵格內徑流向流域出口匯集的路徑，最終影響了整個流域的排水網絡空間分布結構，表現為各河段的走勢、長短的變化。DEM柵格尺度增大對排水網絡的影響，主要表現在減小了河段長度，對比降影響不明顯。

在其他條件完全相同的情況下，河段長度的減小，將縮短流域內各點向流域出口匯流的時間，導致模擬的洪峰出現時間提前，峰值增大。DEM柵格坡度及流域平均坡度的減小，又會增大流域各點的匯流時間，導致峰現時間推遲，峰值減小?？梢?，河段長度變化與坡度變化對徑流過程的影響有相互抵消作用，可認為該抵消作用導致了圖3中計算單元邊長為500 m和800 m以及1 000 m和1 500 m的模擬結果相似。根據計算單元邊長對模擬過程的影響結果(圖3)可以看出，基于不同分辨率DEM模型的模擬結果之間之所以存在差異，主控因素為柵格坡度的不同(外部直接表現為提取的流域平均坡度的不同)。圖3模擬的洪水過程分別在計算單元邊長為300 ～500 m之間和800～1 000 m出現的突變與圖5(b)中坡度隨DEM柵格尺度變化的兩個突變點完全吻合，證明了坡度對模擬結果影響的主控作用。

從模型結構及計算方法來看，計算單元坡度變化通過影響模型中計算的流速，進而影響單元間徑流的轉換強度，最終影響了流域出口的流量過程。

4 計算單元邊長對模型影響的矯正方法

4.1 矯正方法

國內外廣泛采用的坡面及河道流速計算公式中的下墊面糙率都是針對坡地與河道的真實坡度或比降(對應高分辨率DEM)，因此，對于基于某一低分辨率DEM的模型，在直接采用提取的坡度或比降值計算坡面或河道內流速時，必須增加一個反映坡度影響的矯正因子，對非真實坡度或比降進行矯正，計算公式為

(1)

式中：Vc——計算單元內水體積；rc——計算單元內時段降水量；Qc——計算步長內計算單元的周邊來水；s0——計算單元地表坡度；γ——邊坡角度；nc——計算單元地表曼寧糙率系數；X——計算單元邊長；J——坡度矯正系數，用以矯正因DEM空間分辨率變化引起的坡度提取誤差，J的大小與流域平均坡度提取誤差(用基于兩個DEM提取的流域坡度的相對比值表示)有關。在實際應用時，必須有J與流域坡度比值間的關系曲線，該曲線需要由實測水文氣象資料的流域分析得到。分析過程為：首先率定基于不同分辨率DEM的水文模型的J，然后由DEM分辨率根據歸一化的流域坡度與DEM柵格尺度間的關系線得到流域坡度比，最后得到J值與流域坡度比之間的關系曲線。

可以看出，使用提出的矯正方法時需要兩條關系曲線，第一條曲線是歸一化的流域坡度與DEM柵格尺度間的關系線，第二條曲線是矯正系數與流域坡度比之間的關系線。

4.2 實例驗證

將由計算單元邊長為2 000 m模型率定的參數經過矯正后用于計算單元邊長為400 m模型，以檢驗本文提出的矯正方法的有效性。采用4場洪水過程來分析J與流域坡度比值間的關系。由于率定參數時使用的是計算單元邊長為2 000 m的模型，所以分析出的是相對單元尺度為2 000 m模型的矯正系數。計算結果如表4和圖8所示。由表4和圖8可以看出，不同大小計算單元的模型之間的矯正系數差別明顯。

表4 矯正系數與計算單元邊長及流域坡度比值關系

圖8 2 000 m基準流域坡度比與 矯正系數關系Fig.8 Relationship between slope ratio and correction coefficient in 2000 m watershed

由圖5可得400 m分辨率與30 m分辨率的平均坡度間的比值為0.30，2 000 m分辨率與30 m分辨率的平均坡度間的比值為0.11，則400 m分辨率與2 000 m分辨率的平均坡度間的比值為2.7。根據圖8可得計算單元邊長為400 m的模型對于計算單元邊長為2 000 m模型的坡度矯正系數為0.222。

將由計算單元邊長為2 000 m的模型率定所得參數以及矯正系數0.222用于計算單元邊長為400 m的模型，模擬結果如圖9所示。由圖9可以看出，矯正后的模擬結果無論是流量過程線、洪峰流量還是峰現時間，均更加接近實測值，顯著改善了模型的模擬效果。

圖9 矯正前后模擬洪水過程對比Fig.9 Comparison of simulated flood process before and after correction

4.3 矯正方法的使用

基于提出的矯正方法，在無資料大流域中可以按照以下步驟構建分布式流域水文模型：

a.目前較高分辨率DEM(如30 m×30 m)可以免費獲得，更高分辨率的DEM也可以得到。借助GIS工具，經過再采樣可以得到不同分辨率情況下的研究流域，進而得到歸一化的流域坡度與DEM柵格尺度之間的關系曲線。

b.假定高分辨率DEM(如30 m×30 m或更小尺度)代表流域真實的地表特征，對應的分布式流域水文模型的糙率可直接取用相關文獻中的數值。

c.借用有實測水文資料的流域，分析計算單元邊長為30 m模型的矯正系數與坡度比值之間的關系曲線。

d.為減少計算量，增加模型運行的時效性，選擇較低分辨率DEM構建分布式流域水文模型，根據模型計算單元邊長確定矯正系數值，同時根據流域的下墊面條件從相關文獻中查出對應的糙率及其他參數值，最后將矯正系數及參數值代入模型中，完成具有物理基礎的分布式流域水文模型的構建。

5 結 語

本文研究了相同參數情況下計算單元邊長對具有物理基礎的分布式流域水文模型模擬結果的影響，根據DEM分辨率與提取的流域特征之間關系，分析了計算單元邊長對模擬結果的影響機理。由于不同流域的真實地形坡度不同，使得隨著DEM分辨率的變化，提取的流域坡度變化絕對值差異明顯，但經過本文提出的歸一化處理后，由DEM提取的流域坡度隨DEM分辨率變化的相對值非常接近?？梢姡煌饔蜷g流域坡度受DEM分辨率的影響具有極強的規(guī)律性。對于同一流域，相同地面糙率時分布式流域水文模型的模擬效果受計算單元邊長大小的影響明顯，主要是由計算單元對流域地形坡度均化作用的不同引起，且通常情況下不同計算單元邊長的分布式流域水文模型之間糙率系數不宜直接移用。。

在對流域坡度與DEM分辨率關系研究的基礎上，提出了可以消除計算單元邊長對分布式流域水文模型影響的矯正方法，該方法可以在很大程度上消除由柵格均化作用導致的分布式流域水文模型的模擬誤差，因此，能夠有效地解決具有物理基礎的分布式流域水文模型在無資料流域的應用問題。