李致家，臧帥宏 ，劉志雨，徐 杰，董鵬飛，馬亞楠

(1.河海大學(xué)水文水資源學(xué)院，江蘇 南京 210098； 2.水利部水文局，北京 100053；3.煙臺市福山區(qū)水利局，山東 煙臺 265500)

新安江模型是由趙人俊[1-2]提出的概念性水文模型，在中國濕潤、半濕潤地區(qū)的降雨徑流模擬及預(yù)報中得到廣泛應(yīng)用，并取得了良好的效果[3]。目前新安江模型在河道匯流計算中采用的是Muskingum分段連續(xù)演算法[4]，一般適用于地勢較陡的地區(qū)，另外，其中的參數(shù)(K、x)沒有考慮時間和空間上的變異性，不能隨著河道特征以及水力要素的變化而變化，因此有必要對新安江模型河道匯流進(jìn)行研究，增加新的匯流方案，進(jìn)一步增強模型的適用范圍。

河道流量演算以圣維南方程組[5]為理論基礎(chǔ)，利用上斷面流量過程演算下斷面流量過程。由于圣維南方程的復(fù)雜性，目前無法求得其精確解析解，因此在實踐中常采用近似的計算方法。目前常用的河道演算方法可分為2種：一種是基于河段水量平衡方程和槽蓄方程的水文學(xué)方法，如Muskingum法[6-9]；1938年MCCarthy提出一種原始的“水文”洪水演算方法，即Muskingum法，其中所有的地形特征和河段的水力學(xué)特性歸為模型的一系列參數(shù)。1969年Cunge在Muskingum法的基礎(chǔ)上提出了Muskingum-Cunge(MC)演算法[10]，將Muskingum演算法擴(kuò)展到時間可變參數(shù)。1978年P(guān)ouce和Yevjevich在Cunge的基礎(chǔ)上，考慮到波速的變動問題，建議采用變動參數(shù)表示動力波速和水面寬度，提出了變動參數(shù)MC方法[11]。但在變動參數(shù)的實際應(yīng)用中存在質(zhì)量不守恒問題，由此，Todini對變動參數(shù)MC模型引入修正因子對模型進(jìn)行修改，提出Muskingum-Cunge-Todini(MCT)方法[12]。另一種是直接求解完全圣維南方程組或其簡化方程組的方法，即水力學(xué)法，為水文學(xué)方法提供了理論學(xué)基礎(chǔ)。隨著計算機(jī)計算速度提高，求解圣維南方程組的速度有了飛速提高，所以基于圣維南方程的各種簡化形式也得到了廣泛應(yīng)用[13-15]，例如運動波、擴(kuò)散波方法。運動波不考慮動力方程中的慣性項和附加比降項，適用于坡度較陡的山區(qū)，而擴(kuò)散波不考慮慣性項，能反映回水影響，可用于平原河網(wǎng)地區(qū)[16]，所以擴(kuò)散波比運動波適用范圍更廣。

擴(kuò)散波以及MCT這2種方法可以較好地反映洪水波在河道中的運動特點，2種方法均能夠根據(jù)河道斷面特征及河道的下墊面特征來推求河道斷面流量，與現(xiàn)行求解完全圣維南方程組的各種方法相比，這2種演算法解法相對容易，所需資料相對較少。因此，擴(kuò)散波演算法和MCT方法既能保證計算精度，又便于求解與編程。本文采用這2種方法，同時考慮旁側(cè)入流建立具有物理基礎(chǔ)的數(shù)字河道匯流模型，將其與原來的新安江模型進(jìn)行了耦合，對新安江模型進(jìn)行改進(jìn)，改進(jìn)后模型命名為XAJ-DCH模型(Xin’anjiang Digital Channel Model)，并將XAJ-DCH模型應(yīng)用到呈村流域來驗證模型的合理性和適用性。

1 研究流域和數(shù)據(jù)

圖1 研究流域示意圖Fig.1 Map of study catchments

呈村流域位于安徽省黃山市，屬于錢塘江流域，鄰近中國東南沿海，位于亞熱帶季風(fēng)氣候區(qū)，年平均溫度17 ℃，冬季盛行西北風(fēng)，天氣晴冷干燥；夏季多東南風(fēng)，氣溫高，光照強，空氣濕潤；春秋兩季氣旋活動頻繁，冷暖變化大。春季及初夏多鋒面雨，夏秋之際多臺風(fēng)，季風(fēng)環(huán)流的方向與主要山脈走向基本正交，山脈起著阻滯北方寒流和臺風(fēng)的作用。年平均降水量1 600 mm，其中4—6月多雨，占50%，易發(fā)生洪澇災(zāi)害；7—9月，占20%，旱災(zāi)頻繁。河川徑流年內(nèi)、年際變化較大，是典型的濕潤地區(qū)。呈村流域洪水具有漲落陡、洪峰高、歷時短等特點。呈村流域總面積為298.024 km2，如圖1所示，共含有7個雨量站，分別是樟源口、汪村、棣甸、左龍、田里、馮村、大連，1個水文站為呈村站。數(shù)字高程資料(DEM)來自美國地質(zhì)調(diào)查局(USGS)免費提供的全球90 m×90 m的原始DEM數(shù)據(jù)，主要用于填洼、確定流向、水系提取等處理。 高程變化范圍為262～1 614 m，地勢南高北低。

2 模型構(gòu)建及應(yīng)用

2.1 河道演算順序及河道斷面概化

圖2 河道計算順序編碼圖Fig.2 Sequentially coding diagram of river channel computation

首先在數(shù)字高程模型的基礎(chǔ)上，利用Arcgis軟件提取流域的填洼、流向、匯流累積、水系等，采用自然子流域法劃分子流域。首先找到每個子流域的出口點位置，然后根據(jù)流向文件確定每個子流域出口點到整個流域出口點之間的河道，這部分河道用于河道演算；將提取的整個流域的出口點設(shè)為編號1，按照流向文件，根據(jù)D8法原則，從編號為1的柵格開始向上游搜索入流柵格，依次搜索設(shè)置編號2，3，…直至各支流河鏈的起始點，如果在搜索中發(fā)現(xiàn)有2個入流柵格，則這2個入流柵格編號相同。如圖2所示，編號1為流域出口點，紅色星標(biāo)代表著上游子流域的出口點，在進(jìn)行河道演算時，從上游柵格最大編號開始計算，尋找每個編號出現(xiàn)的次數(shù)，對相同編號的河段依次進(jìn)行河道演算。本文假設(shè)計算河道是矩形河道，同時考慮旁側(cè)入流的影響。河道中每個斷面的河寬計算采用Topkapi模型中河寬計算公式[17]在進(jìn)行河道演算時，根據(jù)流向，假設(shè)上下游2個柵格的中心點代表河道斷面，則兩者之間的直線距離就代表該河段之間的距離。

2.2 模擬方法

2.2.1 Muskingum-Cunge-Todini(MCT)方法

MCT方法是在變參數(shù)Muskingum-Cunge方法基礎(chǔ)上引入校正因子β、校正后的庫朗數(shù)C*和雷諾數(shù)Re來保證計算過程中的質(zhì)量守恒[10]，計算公式如下：

(1)

(2)

(3)

(4)

O2=C1I2+C2I1+C3O1+C4

(5)

2.2.2 擴(kuò)散波方法

忽略動力方程中的慣性項，假定河道演算為紊流形式，流量計算采用曼寧公式：

(6)

(7)

式中：Sf——河道摩阻比降；S0——河底比降;h——河道水深,m；Q——河道總流量，m3/s；n——河道的曼寧糙率系數(shù)；A——流斷面面積，m2；R——水力半徑,m。

網(wǎng)格單元k點的河底比降、摩擦比降向前一個網(wǎng)格l點進(jìn)行向前差分作近似計算，計算公式如下：

(8)

(9)

式中：S0k——河底比降；Sfk——k點的摩擦比降；E——河道底部高程，m；Δx——k、l這2個網(wǎng)格點之間的距離，m。

2.3 子流域與河道匹配問題

實現(xiàn)子流域與計算河道之間的耦合是一個關(guān)鍵問題。在時間尺度上，新安江模型產(chǎn)流計算和坡面匯流計算是以1 h為時間步長，而在用擴(kuò)散波方法或者M(jìn)CT方法進(jìn)行河道演算時，為了保證計算的穩(wěn)定性條件，是以秒為時間步長，因此涉及一個時間尺度轉(zhuǎn)換的問題，在這里采用的是將1 h的產(chǎn)流量均分。

圖3 用擴(kuò)散波方法和MCT方法進(jìn)行河道演算的河道部分Fig.3 Extracted channel routed by diffusion wave and MCT method

在空間連接上，主要把子流域進(jìn)行2種劃分，第一種是位于計算河道上游的子流域，如圖3所示，左龍、馮村、田里、大連這幾個子流域，在進(jìn)行三水源劃分之后，用滯后演算法進(jìn)行子流域河網(wǎng)匯流演算，得到每個子流域出口點的流量，作為起始河段的上游斷面的輸入流量；第二種是計算河道所在的子流域，如圖3所示的棣甸、樟源口、汪村、呈村，這些子流域進(jìn)行三水源劃分之后，根據(jù)每個子流域中河道的柵格數(shù)量，將水量進(jìn)行均分，平攤到每個河段中，作為側(cè)向入流。

2.4 模型率定

為檢驗XAJ-DCH模型的適用性，選取呈村流域的歷史洪水進(jìn)行模擬分析。在選擇洪水時保證大、中、小洪水全面選擇，并且各種峰形及歷時長短的洪水盡量覆蓋。為了保證率定參數(shù)的準(zhǔn)確性，將大約2/3的場次洪水用于參數(shù)率定，然后將剩下的進(jìn)行檢驗。共選擇了1986—1998年共23場洪水進(jìn)行模擬，前15場用于率定期，后8場用于驗證。

3 結(jié)果及討論

3.1 參數(shù)率定結(jié)果

為保證計算穩(wěn)定性，擴(kuò)散波方法河道演算以2 s為時間步長，MCT方法以600 s為時間步長。在使用擴(kuò)散波方法時，對于一場洪水的模擬基本控制在20 min以內(nèi)，使用MCT方法時基本控制在2 min以內(nèi)，所以XAJ-DCH模型可用于實時預(yù)報。由表1可以看出，新安江模型中共含有17個參數(shù)，而XAJ-DCH模型在此基礎(chǔ)上增加了一個河道曼寧糙率系數(shù)n。因為2個模型只是匯流模塊不一致，所以2個模型除了匯流參數(shù)，其他的參數(shù)率定值保持一致。河網(wǎng)蓄水消退系數(shù)Cs[18]是子流域河網(wǎng)匯流計算中一個極其敏感的參數(shù)，反映流域退水速率的快慢。XAJ-DCH模型中擴(kuò)散波方法的Cs率定值和MCT方法Cs率定值分別為0.896、0.890，兩者的Cs率定值相近，而新安江模型中Cs率定值為0.780，比XAJ-DCH模型中所采用的2種方法的Cs率定值偏小。Cs通過影響子流域的河網(wǎng)調(diào)蓄作用控制洪峰的形狀，n通過影響河道調(diào)蓄作用對洪峰形狀起控制作用。這2個參數(shù)共同對洪峰起控制作用，兩者處于相互協(xié)調(diào)、相互制約的關(guān)系。因為采用不同的河道匯流方法對河道調(diào)蓄作用也存在著差異，為了達(dá)到相同的洪水模擬效果，所以Cs值也會存在差異。

表1 XAJ和XAJ-DCH模型參數(shù)率定值Table 1 Calibrated parameters of XAJ model and XAJ-DCH model

注：Ke為流域蒸散發(fā)折算系數(shù)；Wum為上層張力水蓄水容量；Wlm為下層張力水蓄水容量；C為深層蒸散發(fā)折算系數(shù)；B為張力水蓄水容量曲線次方；Wm為流域平均張力水蓄水容量；Im為不透水面積占全流域面積的比例；Sm為表層自由水蓄水容量；Ex為表層自由水蓄水容量曲線的指數(shù)；Kg為表層自由水蓄水庫對地下水的日出流系數(shù)；Ki為表層自由水蓄水庫對壤中流的日出流系數(shù)；Cg為地下水消退系數(shù)；Ci為壤中流消退系數(shù)；Lag為滯后時間；k為Muskingum法洪水波在河段內(nèi)傳播時間；x為Muskingum法權(quán)重系數(shù)。

3.2 模型模擬結(jié)果分析

根據(jù)GB/T 22482—2008《水文情報預(yù)報規(guī)范》，選擇洪量合格率、洪峰合格率、確定性系數(shù)以及峰現(xiàn)時間誤差4個指標(biāo)作為評價標(biāo)準(zhǔn)，洪量相對誤差及洪峰相對誤差不超過20%為合格，確定性系數(shù)大于等于0.7為合格，峰現(xiàn)時間誤差小于等于3為合格[19]。圖4中C1、C2、C3分別代表率定期Muskingum法、MCT方法、擴(kuò)散波方法的模擬情況，V1、V2、V3分別代表驗證期使用原新安江模型中Muskingum法、XAJ-DCH模型MCT方法和擴(kuò)散波方法的模擬情況。表2列出了采用3種演算方法得到的率定期及驗證期所有次洪模擬特征值。圖5為呈村1986061909號洪水用以上3種演算方法得到的洪水流量過程線。

圖4 呈村流域洪峰相對誤差、洪量相對誤差、確定性系數(shù)箱線圖Fig.4 Box plots of the relative peak discharge error, relative runoff volume error and Nash-Sutcliffe efficiency statistics of all hourly simulations in Chengcun catchment

圖5 呈村1986061909號洪水過程線Fig.5 Observed and simulated hydrographs of flood No. 1986061909 by using different channel routing methods

表2 次洪模擬結(jié)果特征值統(tǒng)計Table 2 Characteristic statistics of the results of flood event simulation %

由圖4可以看出，呈村流域率定期和驗證期3種演算方法的洪峰相對誤差處于-20%～20%之間，洪量相對誤差處于-10%～10%之間，確定性系數(shù)基本上都大于0.7，并且3種方法的洪峰相對誤差、洪量相對誤差、確定性系數(shù)模擬的結(jié)果相似。總體來說，3種方法的模擬結(jié)果均在可以接受的范圍。

4 結(jié)論及展望

4.1 結(jié)論

a. 新安江模型中的馬斯京根法計算簡單，計算效率高，只要有充足的實測資料，通過率定參數(shù)可獲得良好的模擬效果，但新安江模型沒有考慮河道參數(shù)的時空變異性。而XAJ-DCH模型中采用的MCT方法和擴(kuò)散波方法考慮了河道的時空變異性，如河寬、河道坡降、初始水深等，并且XAJ-DCH模型可同時計算出各河道斷面的水位和流量。

b. XAJ-DCH模型簡單，運行時間短，可用于實時洪水預(yù)報。

c. XAJ-DCH模型中加入的擴(kuò)散波方法考慮了回水的影響，可用于平原地區(qū)，適用范圍比馬斯京根法廣。

4.2 展望

對于新安江模型河道匯流的改進(jìn)還有諸多因素考慮不周，許多問題有待下一步研究。

a. 目前模型只概化為矩形河道，對于其他河道形狀未給予考慮，在接下來要做的是考慮建立梯形河道、三角形河道并且考慮河漫灘情況，使用戶可以根據(jù)實際的河道斷面幾何形狀對模型中的河道概化形狀進(jìn)行選擇。

b. 對于計算出的水位，由于沒有實測資料進(jìn)行對比，無法得知模擬的準(zhǔn)確性。

c. 因為MCT方法和擴(kuò)散波方法的差分不同，所以兩者需要不同的河道演算時間步長，對于2種方法的穩(wěn)定條件有待于進(jìn)一步研究。