基于柵格型新安江模型的中小河流精細化洪水預報

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

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

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

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

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

1 模 型 原 理

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

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

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

(1)

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

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

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

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

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

(2)

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

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

(3)

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

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

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

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

B=δfte+B0

(4)

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

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

2 模型參數空間分布估計

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

根據現有針對GXM參數及其空間分布估計的研究方法[17,26-27]，植被葉面指數、植被高度與坡面匯流曼寧糙率系數等參數，可以根據網格單元植被類型，利用已有研究結論和文獻，直接進行賦值；網格單元張力水蓄水容量、自由水蓄水容量、壤中流與地下徑流出流系數等參數，均可以通過其與地形指數、土壤類型、土層厚度等特征的定量關系進行估算；河道網格單元對應的河道寬度、河道長度等參數，可直接利用DEM、遙感影像等進行計算或測量；地表徑流匯入河道比例參數可采用基于流向的面積比例法進行估算。在此基礎上，再利用實測水文資料，對蓄滿率、蓄滿率臨界狀態等剩余參數進行率定或對估算參數進行微調，即可獲得每個網格單元的模型參數組。

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

在GXM應用中，所需的降雨場主要包括實況降雨場和預報降雨場。其中，實況降雨場主要通過兩種途徑獲取：一種是直接利用地面雨量站實時數據，采用泰森多邊形、反距離平方等方法進行空間插值；另一種則是利用雷達數據、衛星數據與地面雨量站數據，采用混合地理加權回歸和基于貝葉斯與機器學習的多源信息融合等方法[2]，以進一步提高實況降雨場的時空分辨率。對于預報降雨場而言，模型現階段主要采用WRF模式以及EC、NCEP、GRAPES等數值預報產品[28]。

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

4 實 例 應 用

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

4.1 研究區域概況

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

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

4.2 結果分析

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

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

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

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

5 結 語

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