計(jì)算單元邊長對分布式流域水文模型模擬結(jié)果的影響

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

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

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

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

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

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

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

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

1 研究區(qū)概況

以淮河流域息縣水文站以上部分(簡稱淮河流域)為研究區(qū)，地形情況如圖1所示。所用DEM數(shù)據(jù)來源于地理空間數(shù)據(jù)云(http：//www.gscloud.cn)提供的ASTER GDEM數(shù)據(jù)，柵格分辨率為30 m×30 m，所用其他分辨率DEM均為由ARCGIS工具再采樣生成得到。由30 m×30 m分辨率DEM提取的流域面積為10 050 km2。為研究需要進(jìn)行了子流域劃分，劃分結(jié)果如圖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 子流域劃分結(jié)果Fig.2 Sub watershed Division result

2 計(jì)算單元邊長對TOKASIDE模擬結(jié)果的影響

2.1 對模擬徑流過程的影響

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

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

2.2 對水量平衡要素的影響

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

表1 201602號洪水徑流成分統(tǒng)計(jì)

3 計(jì)算單元邊長對模擬結(jié)果影響機(jī)理分析

3.1 對提取的流域坡度的影響機(jī)理

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

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

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

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

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

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

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

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

3.2 對提取的河段特征的影響機(jī)理

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

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

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

3.3 計(jì)算單元邊長的影響分析

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

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

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

從模型結(jié)構(gòu)及計(jì)算方法來看，計(jì)算單元坡度變化通過影響模型中計(jì)算的流速，進(jìn)而影響單元間徑流的轉(zhuǎn)換強(qiáng)度，最終影響了流域出口的流量過程。

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

4.1 矯正方法

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

(1)

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

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

4.2 實(shí)例驗(yàn)證

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

表4 矯正系數(shù)與計(jì)算單元邊長及流域坡度比值關(guān)系

圖8 2 000 m基準(zhǔn)流域坡度比與 矯正系數(shù)關(guān)系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。根據(jù)圖8可得計(jì)算單元邊長為400 m的模型對于計(jì)算單元邊長為2 000 m模型的坡度矯正系數(shù)為0.222。

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

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

4.3 矯正方法的使用

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

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

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

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

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

5 結(jié) 語

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

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