姚 成，李致家，張 珂，朱躍龍，劉志雨，黃迎春，龔珺夫，張錦堂，童冰星

(1.河海大學(xué)水文水資源學(xué)院，江蘇 南京 210098；2.河海大學(xué)計(jì)算機(jī)與信息學(xué)院，江蘇 南京 210098；3.水利部信息中心，北京 100053；4.安徽省水文局，安徽 合肥 230022)

我國(guó)中小河流眾多，洪水頻發(fā)，災(zāi)害嚴(yán)重，是當(dāng)前防洪工作的重點(diǎn)薄弱環(huán)節(jié)。據(jù)統(tǒng)計(jì)，一般年份中小河流的水災(zāi)損失占全國(guó)水災(zāi)總損失的70%～80%，2010年以來水災(zāi)造成的人員死亡中有2/3以上發(fā)生在中小河流[1]。目前，我國(guó)中小河流的洪水預(yù)報(bào)預(yù)警整體水平不高，尚不能滿足防洪需求。水文實(shí)測(cè)、遙感遙測(cè)、DEM、土壤植被等組成的水文大數(shù)據(jù)日益豐富，研究水文大數(shù)據(jù)驅(qū)動(dòng)的精細(xì)化洪水預(yù)報(bào)已成為提高中小河流預(yù)報(bào)精度的必然途徑[2]。

精細(xì)預(yù)報(bào)要求預(yù)報(bào)模型或方法能夠精準(zhǔn)描述中小河流產(chǎn)匯流過程、考慮人類活動(dòng)影響，同時(shí)還能夠提供高分辨率、多要素預(yù)報(bào)產(chǎn)品，以滿足新時(shí)期影響預(yù)報(bào)與風(fēng)險(xiǎn)預(yù)警的業(yè)務(wù)需求[3]。大量研究表明，分布式水文模型是實(shí)現(xiàn)精細(xì)預(yù)報(bào)的重要支撐[4-11]，如，2015年美國(guó)國(guó)家海洋和大氣管理局聯(lián)合多個(gè)科研機(jī)構(gòu)與政府部門，以分布式模型WRF-Hydro為基礎(chǔ)，共同研發(fā)了美國(guó)國(guó)家水模型，美國(guó)天氣預(yù)報(bào)中心在2016年的業(yè)務(wù)預(yù)報(bào)中，對(duì)該模型進(jìn)行了初步運(yùn)行，實(shí)現(xiàn)了覆蓋美國(guó)大陸1 km網(wǎng)格上的連續(xù)計(jì)算，預(yù)報(bào)成果不僅包括美國(guó)270萬條河流斷面的流量與流速變化，也包括土壤含水量等變量的空間分布，進(jìn)一步驗(yàn)證了分布式模型在精細(xì)化洪水預(yù)報(bào)中的有效性[3]；2017年，李致家等[12]以正交網(wǎng)格分布式模型為基礎(chǔ)，提出了基于蓄滿產(chǎn)流與超滲產(chǎn)流空間組合框架的精細(xì)化降雨徑流模型，考慮參數(shù)的地理空間分布，建立了模型參數(shù)與下墊面特征間的定量關(guān)系，在濕潤(rùn)、半濕潤(rùn)半干旱流域均取得了良好的實(shí)際應(yīng)用效果。

本文以柵格型新安江模型(Grid-Xin’anjiang model，GXM)為基礎(chǔ)[13-17]，采用考慮中小河流下墊面條件復(fù)雜、人類活動(dòng)影響大等特點(diǎn)的精細(xì)化產(chǎn)匯流模擬方法，并對(duì)參數(shù)及其空間分布進(jìn)行客觀估計(jì)，旨在進(jìn)一步提高中小河流洪水預(yù)報(bào)精度，實(shí)現(xiàn)流域內(nèi)任意網(wǎng)格單元流量、土壤含水量等不同水文要素的精細(xì)預(yù)報(bào)，為中小河流影響預(yù)報(bào)與風(fēng)險(xiǎn)預(yù)警業(yè)務(wù)提供更豐富的預(yù)報(bào)產(chǎn)品。

1 模 型 原 理

GXM是以正交網(wǎng)格為計(jì)算單元，以新安江模型為基礎(chǔ)構(gòu)建的分布式水文模型，該模型由植被冠層截留模塊、河道降水模塊、三層蒸散發(fā)模塊、蓄滿產(chǎn)流模塊、分水源模塊、坡地匯流模塊與河道匯流模型組成。在產(chǎn)流與分水源計(jì)算時(shí)，新安江模型采用只有統(tǒng)計(jì)意義的概念性元件，即蓄水容量分布曲線，來考慮土壤含水量空間分布不均的問題。而GXM則是利用地形地貌以及土壤植被等下墊面分布特征，確定任意網(wǎng)格單元的張力水蓄水容量與自由水蓄水容量，進(jìn)而對(duì)土壤含水量的空間分布進(jìn)行具體描述。同時(shí)，為了更好地考慮中小河流下墊面條件復(fù)雜、人類活動(dòng)影響大等特點(diǎn)，GXM采用了更具針對(duì)性的產(chǎn)匯流模擬方法。

1.1 沿程水流再分配過程模擬

GXM在進(jìn)行逐網(wǎng)格匯流演算時(shí)，對(duì)沿程水流的再分配過程進(jìn)行了模擬。即對(duì)于任意網(wǎng)格單元而言，當(dāng)其土壤含水量未達(dá)到田間持水量之前，上游來水優(yōu)先補(bǔ)充該網(wǎng)格的土壤缺水量，直至其蓄滿為止。此外，若網(wǎng)格單元有河道存在，屬于河道網(wǎng)格時(shí)，則該網(wǎng)格地表徑流先按一定的比例匯入河道，然后再匯至下游網(wǎng)格。其中，網(wǎng)格單元上游來水的計(jì)算方法為

(1)

式中：Qsup,i——當(dāng)前網(wǎng)格i的上游網(wǎng)格入流量；Qsout,j——上游第j個(gè)網(wǎng)格地表徑流的出流量；fch,j——上游第j個(gè)網(wǎng)格地表徑流出流量匯入河道的比例，對(duì)于坡地網(wǎng)格，該比例取0；m——與網(wǎng)格i相鄰的上游網(wǎng)格總數(shù)。

1.2 庫(kù)塘壩蓄泄影響的定量化模擬

中小河流內(nèi)中小型水庫(kù)、水塘和水壩等蓄水工程種類多、數(shù)量大，且大都無閘門控制、缺乏實(shí)測(cè)水文資料，如何在作業(yè)預(yù)報(bào)中考慮庫(kù)塘壩蓄泄影響是中小河流洪水預(yù)報(bào)面臨的重要問題。GXM首先根據(jù)庫(kù)塘壩經(jīng)緯度坐標(biāo)，確定其所處的網(wǎng)格位置，建立庫(kù)塘壩的拓?fù)浣Y(jié)構(gòu)；對(duì)于有實(shí)測(cè)資料的庫(kù)塘壩，在模型匯流演算中，以庫(kù)塘壩所處的網(wǎng)格為阻斷節(jié)點(diǎn)，根據(jù)其實(shí)測(cè)資料進(jìn)行蓄泄，然后匯流至下游網(wǎng)格；對(duì)于無資料庫(kù)塘壩，引入蓄滿率θ及蓄滿率臨界狀態(tài)λC，通過設(shè)定蓄泄曲線的方式進(jìn)行蓄泄演算[18]。圖1為場(chǎng)次洪水過程中庫(kù)塘壩時(shí)段蓄泄率隨蓄滿率的變化過程示意圖。當(dāng)θ≤λC時(shí)，庫(kù)塘壩以攔蓄為主；反之，以泄流為主。對(duì)于一場(chǎng)洪水過程而言，庫(kù)塘壩的初始蓄滿率θ0低于λC，此時(shí)庫(kù)塘壩僅攔蓄；隨著降雨增加，θ增加至λC；若繼續(xù)降雨，θ>λC，此時(shí)庫(kù)塘壩表現(xiàn)為泄流。

圖1 無資料庫(kù)塘壩蓄泄關(guān)系示意圖Fig.1 Schematic of water storage and discharge curves of ungauged reservoirs and ponds

對(duì)于任意蓄水工程而言，模型假設(shè)其蓄泄曲線均可用冪函數(shù)表示：

(2)

式中：θt——t時(shí)段庫(kù)塘壩的蓄滿率；ωs,t——t時(shí)段庫(kù)塘壩的蓄水率；ms——攔蓄曲線形狀參數(shù)，反映庫(kù)塘壩的蓄水能力，一般取值為0～1；ωd,t——t時(shí)段庫(kù)塘壩的泄流率；md——泄流曲線形狀參數(shù)，反映庫(kù)塘壩的泄流能力，一般取值為0～3。θ0可通過模型的連續(xù)計(jì)算，利用土壤含水量和前期累積降雨量進(jìn)行計(jì)算[18]；也可以通過Landsat、高分一號(hào)等衛(wèi)星影像數(shù)據(jù)，采用最大類間方差迭代提取方法，對(duì)庫(kù)塘壩水體面積進(jìn)行提取[19]，再根據(jù)水體面積-蓄量關(guān)系進(jìn)行計(jì)算。

若令泄流為正，則t時(shí)段的蓄泄量ΔQt為

(3)

式中：QSUM,t——t時(shí)段庫(kù)塘壩所處網(wǎng)格對(duì)應(yīng)的上游地表徑流、壤中流、地下徑流及河道徑流的累積流量。

1.3 考慮河道特征的數(shù)字河道匯流演算方法

GXM在建立之初，采用傳統(tǒng)的Muskingum法進(jìn)行逐網(wǎng)格河道匯流演算。為了進(jìn)一步適應(yīng)不同河道特征，并考慮中小河流河道水力要素的時(shí)空變化，后又引入基于擴(kuò)散波(diffusion wave，DW)與Muskingum-Cunge-Todini(MCT)[20]的數(shù)字河道匯流演算方法[21]。與Muskingum法相比，DW與MCT方法可以根據(jù)河道斷面特征，更好地考慮中小河流河道流量、流速、水深等要素的時(shí)空變異性，進(jìn)而更真實(shí)地模擬中小河流河道匯流過程。

GXM采用基于地形特征的中小河流河寬估算方法[22]對(duì)河道進(jìn)行數(shù)字化，以提高河道網(wǎng)格單元的河寬估算精度，計(jì)算方法為

B=δfte+B0

(4)

式中：B——河道網(wǎng)格i的河寬；fte——河道網(wǎng)格單元地形因子；δ——河寬比例系數(shù)；B0——基礎(chǔ)河寬。其中，δ與B0均可通過遙感影像圖量測(cè)或斷面實(shí)測(cè)資料進(jìn)行確定，fte可通過河道網(wǎng)格單元對(duì)應(yīng)的上游累積匯水面積與坡度原點(diǎn)矩進(jìn)行計(jì)算[22]。

由于GXM采用了蓄滿產(chǎn)流原理，因此模型理論上僅適用于濕潤(rùn)半濕潤(rùn)地區(qū)。若將其用于干旱半干旱地區(qū)，可采用蓄滿產(chǎn)流與超滲產(chǎn)流的時(shí)空動(dòng)態(tài)組合框架[12]。若將模型用于地下水超采地區(qū)時(shí)，如海河流域中小河流[23-24]，可在模型中增加基于“帶門檻水庫(kù)”的地下攔蓄模塊[25]。

2 模型參數(shù)空間分布估計(jì)

分布式模型網(wǎng)格單元多達(dá)數(shù)千個(gè)，其參數(shù)在空間上呈不均勻分布，若只依賴于出口斷面實(shí)測(cè)水文資料，難以保證參數(shù)確定和預(yù)報(bào)結(jié)果的合理性。解決該問題的主要途徑是以參數(shù)物理意義為基礎(chǔ)，采用理論推導(dǎo)與數(shù)據(jù)挖掘等技術(shù)，建立參數(shù)與地形、土壤、植被等下墊面特征的定量關(guān)系，對(duì)模型參數(shù)的空間分布進(jìn)行客觀估計(jì)，以降低參數(shù)對(duì)實(shí)測(cè)資料的依賴性，保證參數(shù)估計(jì)的合理性。

根據(jù)現(xiàn)有針對(duì)GXM參數(shù)及其空間分布估計(jì)的研究方法[17,26-27]，植被葉面指數(shù)、植被高度與坡面匯流曼寧糙率系數(shù)等參數(shù)，可以根據(jù)網(wǎng)格單元植被類型，利用已有研究結(jié)論和文獻(xiàn)，直接進(jìn)行賦值；網(wǎng)格單元張力水蓄水容量、自由水蓄水容量、壤中流與地下徑流出流系數(shù)等參數(shù)，均可以通過其與地形指數(shù)、土壤類型、土層厚度等特征的定量關(guān)系進(jìn)行估算；河道網(wǎng)格單元對(duì)應(yīng)的河道寬度、河道長(zhǎng)度等參數(shù)，可直接利用DEM、遙感影像等進(jìn)行計(jì)算或測(cè)量；地表徑流匯入河道比例參數(shù)可采用基于流向的面積比例法進(jìn)行估算。在此基礎(chǔ)上，再利用實(shí)測(cè)水文資料，對(duì)蓄滿率、蓄滿率臨界狀態(tài)等剩余參數(shù)進(jìn)行率定或?qū)浪銋?shù)進(jìn)行微調(diào)，即可獲得每個(gè)網(wǎng)格單元的模型參數(shù)組。

3 降雨場(chǎng)與下墊面空間分布提取

在GXM應(yīng)用中，所需的降雨場(chǎng)主要包括實(shí)況降雨場(chǎng)和預(yù)報(bào)降雨場(chǎng)。其中，實(shí)況降雨場(chǎng)主要通過兩種途徑獲?。阂环N是直接利用地面雨量站實(shí)時(shí)數(shù)據(jù)，采用泰森多邊形、反距離平方等方法進(jìn)行空間插值；另一種則是利用雷達(dá)數(shù)據(jù)、衛(wèi)星數(shù)據(jù)與地面雨量站數(shù)據(jù)，采用混合地理加權(quán)回歸和基于貝葉斯與機(jī)器學(xué)習(xí)的多源信息融合等方法[2]，以進(jìn)一步提高實(shí)況降雨場(chǎng)的時(shí)空分辨率。對(duì)于預(yù)報(bào)降雨場(chǎng)而言，模型現(xiàn)階段主要采用WRF模式以及EC、NCEP、GRAPES等數(shù)值預(yù)報(bào)產(chǎn)品[28]。

在進(jìn)行下墊面空間分布提取時(shí)，主要采用DEM高程、植被覆蓋/利用、土壤類型與水利普查等數(shù)據(jù)，數(shù)據(jù)源包括SRTM數(shù)據(jù)庫(kù)、UMD-LULC數(shù)據(jù)庫(kù)、HWSD數(shù)據(jù)庫(kù)等。同時(shí)，利用DEM高程數(shù)據(jù)，可進(jìn)一步提取得到流向、水系、坡度、坡向、流徑長(zhǎng)度和地形指數(shù)等地形地貌特征。在此基礎(chǔ)上，也可以利用每個(gè)網(wǎng)格單元的流向和上游累積匯水面積，獲得逐網(wǎng)格匯流演算方法所需的網(wǎng)格單元演算次序[15]。

4 實(shí) 例 應(yīng) 用

GXM自建立以來，已成功應(yīng)用于江蘇、浙江、安徽、陜西等不同省份多個(gè)濕潤(rùn)半濕潤(rùn)中小河流的洪水模擬及業(yè)務(wù)預(yù)報(bào)，均取得了良好的應(yīng)用效果[2,13-17]。本文以新安江源頭橫江為例，開展GXM的應(yīng)用分析。

4.1 研究區(qū)域概況

橫江為典型的濕潤(rùn)區(qū)中小河流，位于皖南山區(qū)，其主要控制站屯溪站以上流域面積為2 670 km2，流域內(nèi)植被豐富，森林覆蓋率達(dá)70%以上。屯溪流域降水主要集中在汛期，汛期內(nèi)降水約占全年總降水量的60%以上，暴雨類型以鋒面型暴雨和低壓型暴雨為主。屯溪流域?qū)偕絽^(qū)性河流，地勢(shì)西高東低，高差變化懸殊，河道坡度大，洪水陡漲陡落，流速大、歷時(shí)短。流域內(nèi)中型水庫(kù)1座，庫(kù)容為2 900萬m3。流域上游已建小型水庫(kù)總庫(kù)容2 200萬m3，塘壩總庫(kù)容2 555萬m3，各小型水庫(kù)與塘壩基本無調(diào)度運(yùn)用方案。屯溪流域水系及中小水庫(kù)塘壩分布見圖2。

圖2 屯溪流域水系及中小水庫(kù)與塘壩分布Fig.2 Water system and spatial distribution of small and medium-sized reservoirs and ponds in Tunxi Basin

4.2 結(jié)果分析

GXM在屯溪流域應(yīng)用時(shí)，網(wǎng)格大小采用1 km×1 km，即可以更好地與土壤、植被等數(shù)據(jù)的網(wǎng)格大小統(tǒng)一，也可以在保證應(yīng)用精度的前提下，提高模型運(yùn)算效率[29]。本文利用屯溪站1982—2020年日資料及同期66場(chǎng)洪水過程資料，對(duì)模型參數(shù)及其空間分布進(jìn)行了估計(jì)、率定與檢驗(yàn)。其中，1982—2008年資料用于參數(shù)的估計(jì)與率定，剩余資料用于參數(shù)的檢驗(yàn)。結(jié)果表明，GXM在屯溪流域應(yīng)用的徑流深合格率為90.9%，洪峰合格率為83.3%，峰現(xiàn)時(shí)差合格率為90.9%，確定性系數(shù)均值為0.92。圖3為屯溪站及休寧站2019年、2020年洪水過程預(yù)報(bào)結(jié)果，休寧站為屯溪站的嵌套站?？梢钥闯觯珿XM不僅可以實(shí)現(xiàn)對(duì)屯溪斷面出流過程的高精度模擬預(yù)報(bào)，同時(shí)可以在不修正參數(shù)及其空間分布情況下，直接給出上游嵌套斷面休寧的高精度結(jié)果，也進(jìn)一步驗(yàn)證了參數(shù)空間分布估算的合理性與準(zhǔn)確性。

圖3 屯溪站及休寧站2019年、2020年洪水過程預(yù)報(bào)結(jié)果Fig.3 Forecast results of flood hydrographs in 2019 and 2020 for Tunxi and Xiuning stations

除了預(yù)報(bào)斷面的洪水過程，GXM還可以提供流量、水位、流速等不同水文要素的時(shí)空分布。圖4為流量、缺水量、庫(kù)塘壩蓄量變化量以及流速空間分布的模擬結(jié)果，驗(yàn)證了GXM對(duì)于不同水文要素時(shí)空分布預(yù)報(bào)的能力，相關(guān)預(yù)報(bào)產(chǎn)品可為進(jìn)一步開展中小河流影響預(yù)報(bào)與風(fēng)險(xiǎn)預(yù)警業(yè)務(wù)提供支撐。

圖4 流量、缺水量、庫(kù)塘壩蓄量變化量以及流速的空間分布模擬結(jié)果Fig.4 Simulation results of spatial distribution of discharge,water shortage,storage change of reservoir and pond,and flow velocity

5 結(jié) 語

本文以GXM為基礎(chǔ)，針對(duì)中小河流下墊面條件復(fù)雜、中小水庫(kù)與塘壩數(shù)量大等特點(diǎn)，采用沿程水流再分配方法處理網(wǎng)格單元之間以及網(wǎng)格與河道之間的水流交換，利用庫(kù)塘壩蓄泄影響的定量化模擬方法反映庫(kù)塘壩蓄泄對(duì)產(chǎn)匯流過程的影響，采用數(shù)字河道匯流演算方法模擬中小河流河道特征，并為不同條件下進(jìn)行逐網(wǎng)格匯流演算提供了Muskingum、DW、MCT等多種演算方法。此外，針對(duì)分布式水文模型參數(shù)空間分布估算難題，GXM以參數(shù)與地形地貌、土壤植被等下墊面特征的定量關(guān)系挖掘?yàn)榛A(chǔ)，對(duì)模型參數(shù)的空間分布實(shí)現(xiàn)客觀估計(jì)，降低參數(shù)對(duì)實(shí)測(cè)水文資料的依賴性。以安徽橫江屯溪流域?yàn)槔?，?duì)GXM精細(xì)預(yù)報(bào)效果進(jìn)行驗(yàn)證，結(jié)果表明，GXM不僅可以實(shí)現(xiàn)對(duì)屯溪出口斷面洪水過程的高精度預(yù)報(bào)，也可以實(shí)現(xiàn)對(duì)流域內(nèi)部嵌套斷面洪水過程的高精度預(yù)報(bào)；同時(shí)還可以獲取流量、土壤含水量、流速、水深等不同要素的時(shí)空分布，可以為中小河流洪水防控提供豐富的預(yù)報(bào)產(chǎn)品。