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

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

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

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

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

1 模 型 原 理

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

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

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

(1)

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

1.2 庫塘壩蓄泄影響的定量化模擬

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

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

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

(2)

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

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

(3)

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

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

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

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

B=δfte+B0

(4)

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

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

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

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

根據(jù)現(xiàn)有針對GXM參數(shù)及其空間分布估計的研究方法[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)格單元對應(yīng)的河道寬度、河道長度等參數(shù)，可直接利用DEM、遙感影像等進(jìn)行計算或測量；地表徑流匯入河道比例參數(shù)可采用基于流向的面積比例法進(jìn)行估算。在此基礎(chǔ)上，再利用實測水文資料，對蓄滿率、蓄滿率臨界狀態(tài)等剩余參數(shù)進(jìn)行率定或?qū)浪銋?shù)進(jìn)行微調(diào)，即可獲得每個網(wǎng)格單元的模型參數(shù)組。

3 降雨場與下墊面空間分布提取

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

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

4 實 例 應(yīng) 用

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

4.1 研究區(qū)域概況

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

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

4.2 結(jié)果分析

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

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

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

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