城市小流域暴雨洪水響應本構關系

關鍵詞：GSSHA水文模型；隨機暴雨移置；降雨時空特征；暴雨洪水

城市化是全世界范圍內(nèi)的普遍現(xiàn)象?？焖俪鞘袛U張通過增加不透水面積、改變河網(wǎng)水系連通性等顯著改變了流域下墊面水文水力學屬性[1-2]，使城市地區(qū)極易出現(xiàn)“峰高、量大和歷時短”的洪水過程，洪澇災害頻發(fā)，嚴重威脅社會經(jīng)濟可持續(xù)發(fā)展[3-5]。從流域尺度理解暴雨洪水響應規(guī)律，是城市水文研究的重點和難點，也是洪澇預報預警等城市防洪減災措施實施的科學依據(jù)。

水文模擬是研究流域尺度暴雨洪水響應規(guī)律的重要方法。自Freeze等[6]提出的“FH69藍本”開始，水文模型經(jīng)歷了從集總式模型向分布式物理模型的發(fā)展。分布式物理模型中的數(shù)學物理方程（如Richards方程、Saint-Venant方程）通常僅適用于描述微觀或中觀尺度的水文過程[7-8]，難以刻畫流域尺度暴雨洪水響應特征[9]，流域尺度暴雨洪水響應本構關系尚未得到充分認識[10-11]。本構關系廣義上可理解為“反映物質(zhì)宏觀性質(zhì)的數(shù)學模型”[7]，最早應用在力學領域，綜合反映材料或結構的力學性能。水文本構關系的研究對象是空間異質(zhì)性強的流域及水文過程[7]，主要目的是描述流域尺度產(chǎn)匯流過程需要滿足的質(zhì)量、動量、能量守恒方程及對應流域的幾何關系。前人對此開展了較為豐富的研究，例如，Tian等[12]采用蒙特卡羅法分析發(fā)現(xiàn)土壤水分運動規(guī)律在宏觀和微觀尺度上有明顯差異；郭偉建等[13]基于運動波方程，通過引入模擬尺度內(nèi)水位、流速及地形異質(zhì)性影響，在模擬空間對水文模型進行積分，實現(xiàn)模型尺度與模擬尺度的匹配。國內(nèi)外已有研究主要集中在點尺度土壤水運動、山坡尺度降雨產(chǎn)流機制[7]，缺少對流域尺度暴雨洪水響應本構關系的研究。

暴雨洪水響應受到流域下墊面水文水力學屬性與降雨時空特征的共同影響。已有研究主要關注下墊面水文水力學屬性對洪水過程的影響[5，14]，但降雨時空特征（暴雨移動速度、方向等）對洪水的影響規(guī)律也不可忽視[15]。例如，Yang等[16]研究發(fā)現(xiàn)隨著降雨量級的增大，降雨空間異質(zhì)性對城市小流域洪水量級的影響減弱；李冰雪等[17]通過分析多種臺風降雨移動情景下的城市內(nèi)澇過程，識別了降雨時空特征對城市內(nèi)澇積水深度的影響程度；Zhou等[18]分析了降雨時空變化對城市小流域洪水頻率計算結果的影響。相關研究證實了降雨時空特征對洪水響應過程的重要性，但尚未定量揭示其在流域尺度暴雨洪水響應本構關系中的貢獻[19]。

本文以南京市九鄉(xiāng)河流域為研究區(qū)，使用分布式水文模型GSSHA（GriddedSurface/SubsurfaceHydrologicAnalysis）、隨機暴雨移置和可解釋機器學習方法，識別控制流域洪水響應的致洪暴雨關鍵特征，在流域尺度上構建考慮暴雨時空特征的暴雨洪水響應宏觀本構關系，定量揭示長三角地區(qū)城市小流域暴雨洪水響應規(guī)律，為城市地區(qū)暴雨洪澇預報預警提供科學支撐。

1研究區(qū)概況

九鄉(xiāng)河發(fā)源于南京市江寧區(qū)湯山街道龍尚湖，流經(jīng)江寧和棲霞2個行政區(qū)，在棲霞山的山腳匯入長江，流域面積為133km2，平均海拔為46m，河道總長為21.7km，平均河寬為10m，見圖1。

本文基于1985—2019中國年度土地覆蓋數(shù)據(jù)集（CLCD）[20]和世界土壤數(shù)據(jù)庫（HWSD）[21]分析九鄉(xiāng)河流域下墊面特征。九鄉(xiāng)河上游面積約39km2，以農(nóng)田（占比45.7%）、森林（占比42.2%）為主，不透水面約占10%；下游流經(jīng)南京仙林大學城，河道順直，渠化明顯，不透水面約占32%（圖1（c））。九鄉(xiāng)河流域具有典型亞熱帶季風氣候，年均降水量為1106mm，降水年際變化大；土壤類型以壤土和砂壤土為主，有少量黏土和砂土（圖1（d）），土層較薄，以蓄滿產(chǎn)流機制為主。

2009年以來，九鄉(xiāng)河流域經(jīng)歷了快速城市化進程，暴雨洪澇災害頻發(fā)，是南京重點防洪區(qū)域之一[22]。2021年起，南京大學和南京市江寧區(qū)水務局聯(lián)合在九鄉(xiāng)河流域開展綜合水文試驗，在上游設置控制斷面，通過時差法超聲波測流系統(tǒng)進行流量監(jiān)測，并布設3個翻斗式雨量計進行降雨監(jiān)測。此外，南京市水文局在流域周邊布設有6個雨量站（圖1（b））。

2研究方法

本文在九鄉(xiāng)河流域建立分布式水文模型GSSHA，模擬場次洪水過程；采用隨機暴雨移置方法對雷達降雨資料進行分析，隨機生成大量降雨事件，用于驅(qū)動GSSHA模型，得到對應洪水過程；借助隨機森林模型和可解釋機器學習方法識別控制洪水響應的致洪暴雨關鍵特征，分析致洪暴雨特征與洪峰量級之間的關系，建立流域尺度暴雨洪水響應本構關系模型，并對洪峰量級進行預測。

2.1GSSHA水文模型構建

GSSHA模型是具有明確物理意義的分布式水文模型，由Downer等[23]在CASC2D模型的基礎上改進和開發(fā)。該模型以數(shù)字高程模型為基礎將流域劃分為規(guī)則的網(wǎng)格計算單元，使用有限差分的擴散波形式模擬二維坡面匯流[24]，采用一維運動波模擬河道匯流。該模型已在不同國家和地區(qū)得到較好應用，如沙特阿拉伯麥加流域、中國石馬河流域和黃家河流域等[25-27]。

九鄉(xiāng)河流域上游共有10場暴雨洪水記錄，選擇其中7場為率定期，其余3場為驗證期；選擇納什效率系數(shù)（ENS）為評價指標；采用Levenberg-Marquardt自動率定算法[23]對飽和水力傳導度、土壤孔隙度、初始土壤含水量、河道匯流糙率和坡面匯流糙率共5個敏感參數(shù)進行率定，將最優(yōu)參數(shù)用于全流域暴雨洪水過程模擬。

江蘇省氣象局提供雷達降雨資料（時間分辨率為6min，空間分辨率為1km），用平均偏差校正法[28]對雷達降雨誤差進行校正，校正后的降雨數(shù)據(jù)用于驅(qū)動GSSHA模型。參數(shù)率定結果顯示，場次洪水過程的ENS介于0.7～0.8之間，表明GSSHA模型能夠準確捕捉九鄉(xiāng)河流域暴雨洪水響應特征。

本文設置3種模擬情景：控制情景（S1）為流域當前下墊面狀況；天然情景（S2）將城市土地利用格點全部替換為農(nóng)田，表示城市化之前的下墊面條件；完全城市化情景（S3）將非城市土地利用格點全部替換為城市，表示完全城市化之后的下墊面條件。通過對比不同模擬情景，分析暴雨洪水響應特征對流域下墊面水文水力學屬性的敏感程度。

2.2隨機暴雨洪水事件

隨機暴雨移置是降雨資料時空延展的主流方法之一[29]，采用“以空間換時間”思想，以目標區(qū)域為中心劃定暴雨移置區(qū)，在移置區(qū)內(nèi)選擇一系列強降雨事件組成“暴雨目錄”，結合概率重采樣和地理位置移置的方法延長暴雨序列[18]。本文使用校正的雷達降水數(shù)據(jù)，通過隨機暴雨移置生成大量暴雨過程，驅(qū)動GSSHA模型生成相應洪水過程。具體步驟如下：

（1）確定移置區(qū)。結合研究區(qū)地理位置和氣候條件特征，通過對歷史雨量資料進行分析，選定降雨空間分布最均勻的一個區(qū)域作為本文暴雨移置區(qū)T0，面積約7200km2（圖1（a））。將T0內(nèi)與九鄉(xiāng)河流域大小、形狀和走向基本一致的任意區(qū)域稱為一個“子移置區(qū)”。為簡化計算，本文采用規(guī)則矩形將九鄉(xiāng)河流域進行概化，將該矩形在移置區(qū)內(nèi)遍歷移動，獲得子移置區(qū)集T1。

（2）生成暴雨目錄。計算子移置區(qū)集T1內(nèi)滑動Xh窗口內(nèi)的累計降水量（Ptot），統(tǒng)計每個子移置區(qū)的最大Xh累計降水量（Pmax），篩選Ptot超過Pmax的80%分位數(shù)的Xh降雨過程，組成暴雨目錄S和對應子移置區(qū)集T2。對不同降雨歷時Ptot統(tǒng)計分析發(fā)現(xiàn)，當X取6h時，Ptot變異性最大，表明暴雨目錄的代表性較強。

（3）篩選暴雨目錄。對暴雨目錄和子移置區(qū)集進行篩選，剔除重復的降雨過程。具體做法為：對Ptot降序排列，得到子移置區(qū)集T2'。從T2'（1）開始，逐一檢驗其與剩余子移置區(qū)T2'（2）～T2'（n）是否有重疊，若有重疊，保留T2'（1），排除與之重疊的子移置區(qū)。以此類推，直至T2'中任何2個子移置區(qū)都不相互重疊，得到最終暴雨目錄S1和子移置區(qū)集T3。

（4）提取降雨過程。提取T3中Pmax對應的6h降雨過程，共得到187場暴雨過程，組成暴雨目錄S2。

（5）生成洪水過程。將S2平移到九鄉(xiāng)河流域所在位置，驅(qū)動GSSHA模型，得到187場洪水過程。開展場次洪水過程模擬時，流域前期土壤含水量保持不變。

2.3致洪暴雨特征因子

本文選擇3組指標描述致洪暴雨特征，分別為：①降雨基本特征：降雨累積量（Rtot）、強降雨覆蓋率（A10）、最大降雨強度（Rmax）；②降雨相對流域空間位置特征：標準化降雨權重流距（DRW）、強降雨權重流距（DRWr）；③降雨相對流域運動特征：降雨強度梯度（dR/dt）、降雨權重流距梯度（dDRW/dt）、強降雨權重流距梯度（dDRWr/dt）。

降雨累積量、強降雨覆蓋率的計算公式如下：

式中：DRW取值為0～1，數(shù)值越接近于0，代表降雨空間分布越靠近流域出口。當使用強降雨強度作為權重時，得到標準化強降雨權重流距。本文所有與降雨權重流距有關的特征因子均由標準化降雨權重流距計算所得。通過計算所有降雨事件的平均權重，得到九鄉(xiāng)河流域平均降雨權重流距為0.47km。

降雨強度梯度表示雨強隨時間變化的快慢程度，降雨權重流距梯度表示暴雨相對流域的移動速度和方向，正值表示向流域上游移動，負值則表示向下游移動，絕對值越大表明移動速度越快。二者分別用降雨過程的雨強、DRW隨時間的變化率表示，采用最小二乘法擬合法對斜率進行計算[18]。

2.4暴雨洪水響應本構模型構建

SHAP（ShapleyAddictiveexPlanations）是一種可解釋機器學習方法，由Lundberg等[31]基于博弈論提出，已被應用于分析城市洪水模擬、洪水預警管理等[32-33]影響因子的重要性。本文采用隨機森林模型[34]和SHAP定量識別致洪暴雨特征對洪水響應的貢獻。從暴雨目錄中隨機選擇60%為訓練集，剩余40%為驗證集，建立暴雨特征因子與洪峰量級之間的隨機森林模型，剔除方差膨脹因子大于5的特征因子來消除自變量之間的多重共線性。設置隨機森林決策數(shù)目為100，最小葉子數(shù)為5，以確定性系數(shù)（R2）、平均絕對誤差（EMA）、均方根誤差（ERMS）為模型檢驗指標，多次循環(huán)保證SHAP結果可靠性。

通過分析洪峰與各特征因子之間的關系，構建本構關系模型，隨機選擇80%洪水過程進行參數(shù)率定，其余20%用于模型驗證。采用經(jīng)驗頻率公式對洪峰量級進行重現(xiàn)期計算，以5年一遇為閾值將事件分為低量級和高量級洪水事件，計算擬合精度R2、EMA、ERMS，驗證本構關系模型對不同量級洪水預測的適用性。3種下墊面情景的模擬采用同樣流程處理。

3結果與分析

3.1隨機致洪暴雨事件特征

結果表明，約70%降雨過程的強降雨覆蓋整個九鄉(xiāng)河流域，表明該地區(qū)暴雨空間尺度較大，容易發(fā)生全流域洪水響應。對最大雨強進行頻率分析表明，降雨事件重現(xiàn)期主要集中在5～50a，最極端事件重現(xiàn)期超過200a。

從空間分布上看，暴雨目錄與站點觀測的平均降雨累積量在空間上具有高度一致性，整體呈現(xiàn)下游偏高的分布特征。從量級上看，暴雨目錄結果高于站點觀測，表明通過隨機暴雨移置可以擴展目標流域歷史致洪暴雨的表征（圖2）。65%降雨過程的DRW小于流域平均值，表明大部分過程的暴雨中心靠近流域中下游，與累積降雨量在流域下游偏多的分布特征一致。

3.2場次暴雨洪水響應特征分析

挑選2組典型降雨過程，分析洪水響應特征，結果如圖3所示。A組最大小時雨強相似，洪峰流量（Qp）量級差異較大，其中事件A1的洪峰量級是事件A2的2倍，事件A2雨峰峰現(xiàn)時間比事件A1出現(xiàn)晚；B組洪峰流量量級相似，最大小時雨強相差較大，其中事件B2最大小時雨強為事件B1的2.5倍，事件B1洪水過程線比事件B2更為平緩。

從降雨時程分布特征看，歷時短且時程分布集中的暴雨能產(chǎn)生更高洪峰量級，且洪水過程線形狀更加尖瘦；歷時長且時程分布分散的暴雨所產(chǎn)生的洪水過程線相對平緩（圖3（c）—圖3（d））。從降雨空間運動特征看，事件A1和事件A2的最大雨強相同，但事件A1的dDRW/dt為負值，表明降雨向流域下游移動，通過與洪水波演進發(fā)生“疊加效應”，放大了洪峰量級。事件A1和事件B12場致洪暴雨均向流域下游移動，但事件B1降雨時程分布更均勻，表明穩(wěn)定且持續(xù)的降雨有利于土壤下滲，導致產(chǎn)流和洪水響應較為緩慢，因此事件B1的洪峰量級小于事件A1。

表1顯示了致洪暴雨特征因子與洪峰量級的相關系數(shù)。3種情景中，S1中洪峰量級與Rmax相關系數(shù)最大，其次為Rtot，與DRWr和dDRWr/dt相關性較弱；S3與S1相似；S2中洪峰量級與Rmax、Rtot的相關性減弱，與DRWr、dDRWr/dt的相關性增強。洪峰量級與dDRWr/dt存在明顯負相關關系，表明當暴雨云團移動和洪水演進形成“疊加效應”時，容易產(chǎn)生更高洪峰量級，二者相關性在S2中最強，表明天然條件下，流域產(chǎn)匯流特征的空間異質(zhì)性較強，洪水響應受到致洪暴雨時空結構特征影響較大。

隨機森林模型能夠較為準確地預測洪峰量級，3種情景中訓練集和驗證集的平均R2達0.7～0.9，EMA為15～30m3/s，ERMS為30～60m3/s。使用SHAP法對特征因子的貢獻度排序，發(fā)現(xiàn)Rtot、Rmax、dDRWr/dt的SHAP值明顯高于其他特征因子，表明對隨機森林模擬精度貢獻較大（圖4）。

對于控制情景，高值Rtot和Rmax比DRWr的貢獻度大，表明Rtot和Rmax越大，洪峰量級越大；當降雨中心位于流域上游時，對洪峰量級起到增強作用，但相比于Rtot、Rmax和dDRWr/dt對洪峰的影響，這種影響較小。dDRWr/dt較小時，貢獻度較高，表明DRWr隨時間基本不發(fā)生變化時，洪峰量級越大。對于天然情景，dDRWr/dt對洪峰的貢獻最大，表明洪水響應受降雨空間運動影響變大，可能與自然流域下墊面屬性異質(zhì)性較強有關，城市擴張減弱了下墊面空間異質(zhì)性，導致洪水響應對降雨空間運動的敏感程度變?nèi)酢τ谕耆鞘谢榫?，暴雨洪水響應關系較為簡單，洪峰量級主要由Rtot和Rmax決定，與降雨時空特征關系較弱（圖4（c））。

3.3暴雨洪水本構關系構建及適用性

圖5展示了洪峰量級與單個暴雨特征因子的統(tǒng)計關系。Rmax和Rtot與洪峰量級表現(xiàn)為冪函數(shù)關系，dDRWr/dt與洪峰在形狀上呈現(xiàn)非單調(diào)函數(shù)關系?；诮y(tǒng)計關系提出的流域暴雨洪水響應本構關系模型形式為

基于構建的本構關系模型，對洪峰量級進行預測，結果如圖6所示，預測值與模型模擬值的R2為0.7～0.9，EMR為0.154。天然情景的模擬效果略低于控制情景，完全城市化情景的模擬效果優(yōu)于控制情景。同一情景下，模型對小量級洪水的預測精度低于大洪水，這可能與構建模型時小量級洪水樣本較少有關。

4結論

本研究借助隨機暴雨移置和分布式水文模型GSSHA，獲得九鄉(xiāng)河流域187場洪水過程，通過解析致洪暴雨特征與洪水響應的關系，構建了流域尺度暴雨洪水響應本構關系模型。主要結論如下：

（1）分布式水文模型GSSHA可以準確刻畫九鄉(xiāng)河流域的暴雨洪水響應特征。基于隨機暴雨移置和分布式水文模擬方法生成了致洪暴雨和流域洪水響應過程的大樣本，為缺資料小流域暴雨洪水響應研究提供了新范式。

（2）致洪暴雨時空結構特征是控制流域洪水響應的關鍵因素，城市擴張對流域下墊面屬性的改變在一定程度上弱化了洪水響應對致洪暴雨特征的敏感程度。當流域完全城市化時，暴雨洪水響應關系趨于簡單，洪峰量級主要由累積降水量和雨強控制。

（3）構建了考慮致洪暴雨時空特征的暴雨洪水響應本構關系模型，可以用于預測洪峰量級，預測值與模擬值的平均相對誤差約為0.154，預測精度較高。未來需要在本構關系模型中增加對流域下墊面水文水力學屬性（如不透水面分布、土壤屬性等）的表征，提高本構關系模型的普適性。