引言

全球性城市化过程导致城市人口增多、城市面积迅猛扩张，致使城市成为一个巨大的承灾体。这个巨大的承灾体在面对自然灾害的冲击时，往往不堪一击。其中极端暴雨导致的城市积涝已经并将持续地给全球城市带来冲击及影响，给人类带来巨大的人员及财产损失(扈海波等, 2013；Zhou et al, 2013；Shepherd et al, 2011)。

这促使众多的气象水文专家将研究的重点放在城市洪涝及积涝的风险识别、预警及预报上(Sharif et al，2006；He et al，2013；Moreno et al，2013)。而其中最为主流及最有应用前景的预报预警模式即为结合短时临近预报系统及城市水文模型的开发应用(Ogden et al，2000；Smith et al，2005a；2005b；Vivoni et al，2006；Sharif et al，2006；Ntelekos et al，2006；Wright et al，2012；Barthold et al，2015)。其经典的流程方法是用雷达反演雨量强迫到城市水文模型，实时模拟城市地区未来的积水状况(积水的深度及持续时间)来实施城市暴雨积涝风险预警及监控(赵琳娜等，2012；Smith et al，2005a；2005b；Silvestro et al，2017；Davolio et al，2017；Poletti et al，2019)，甚至是洪涝灾害损失评估及预评估(Meyer et al，2009)。

这种预警模式能否成功应用的关键依赖于城市水文模型模拟输出结果的有效性。此有效性包括积水模拟是否正确以及预警时间提前量能否满足城市应急的需求(Sharif et al，2006)。实质上模型模拟的有效性根本依从于预报雨量的准确度(Silvestro et al，2017；Thorndahl，2017；Poletti et al，2019)。而预报雨量与应用效果之间的矛盾是预报提前量越大，其不确定性也越大，如强迫到水文模型中，效果不会太好(Kober et al，2012；Silvestro et al，2017；Thorndahl et al，2107)。只是预报提前量越大越满足预警时效及应急准备的需求。以短时临近预报系统的定量降水预报(QPFs)为例，此类数据在落区及雨量大小上普遍存有偏差，将其强迫到水文模型中，会产生让现有的基于确定性模式的预报及预警作业无法承受的误差(Kober et al，2012)。但QPF又有较理想的时间提前量，它对预警很关键。用QPF等外推及预报雨量做水文模型强迫的应用不是特别多(Silvestro et al，2017；Poletti et al，2019)。在特定个例情况下，如预报结果较好也会有不错的应用效果，比如用融合了数值预报模式(NWP)雨量场的短时临近预报做强迫的应用(Davolio et al，2017；Nerini et al，2017；Poletti et al，2019)。尽管提高QPF等雨量预报产品的质量对城市积涝预警很重要，但现有的短时临近数值预报产品的质量很难突破到能一如既往地给第三方模型提供满意强迫值的程度(Liechti et al，2013)。

尽管如此，短时临近预报系统中较为精准的雷达降水估计(QPEs)，尤其是经过雨量站校准后，其雨量及落区准确度确实能满足城市水文模型强迫的需求(Anagnostou and Krajewski, 1998)，而模型输出结果也会产生一定的时间提前量。Sharif et al(2006)认为城市水文模型与短时临近预报产品的QPE结合进行城市积涝灾害的预警也是可行的，因为积水过程的延迟效应会导致积水峰值明显落后于降水峰值。根据Sharif et al(2006)的研究，一般的降水过程大概有20~40 min的积水延迟，这就是说即便用实况雨量来做积涝风险预警，也能达到一定的时间提前量。而20~40 min的提前量对一些应急能力较强的城市来说，还是能减缓不少积涝灾害的影响及损失。尽管雨量实时强迫到水文模型的积涝预警模式具有此类优点，但其最终的模型输出结果经由短时临近系统处理，然后强迫到水文模型中，两个模型的加工会加大输出的不确定性(Amponsah et al，2016)。这将制约其在确定性预报及预警模式中的应用效果，尤其是在预报员无法把握模拟结果的情况下，甚至繁琐的数据处理及对模型的不熟练操作，都会干扰到预报员的决策判断而造成模型被严重闲置。Barthold et al(2015)及Herman and Schumacher(2018)等曾强烈呼吁气象界专业人士需要与水文界专业人士合作，共同解决暴雨洪涝及城市积涝预报及预警难题。Barthold et al(2015)甚至认为正是由于积涝预警在技术应用层面上遭遇到水文学与气象学上的难题及挑战，从而使得积涝能为所欲为地成为最为致命的气象灾害，每年导致的人员伤亡甚至超过了台风及龙卷。

除了用雨量实时强迫到城市水文模型作模拟的预警模式外，另一种脱机模式更值得应用及关注。美国国家天气局(NWS)的积涝指导系统(FFGs及GFFG)就采取脱机运行方式，用水文模型反推及验算可造成流域内暴雨积涝的1、3及6 h临界雨量及临界产流值(Thresh-R)并用做预警判断(Ntelekos et al，2006；Barthold et al，2015)。FFG的临界雨量的降雨持续时间分别对应雷达外推雨量的1、3及6 h降雨时长，用于表述该时段流域可导致洪涝或积涝的临界雨量(Georgakakos，1986)。NWS的天气预报办公室这样做的目的是可对照雷达外推雨量做出预警及预报判定，便于实施极端暴雨洪涝及积涝灾害的应急防控(Gourley and Vieux, 2005)。FFGs临界雨量值主要通过计算临界产流来判定，此临界产流由块状地表水文反应模型“Sacramento土壤湿度核算模型”推算(The Sacramento Soil Moisture Accounting model，SAC-SMA)(Carpenter et al，1999)，此模型以流域为核算单元。FFG建立了覆盖全美的数千个不同汇水面积的流域临界雨量(Gourley and Vieux, 2005)。Schmidt et al(2007)提出用格点化的方法替换FFG使其成为GFFG(gridded FFG)。随后NWS的河流预报中心(River Forecast Center，RFC)开发了GFFG并投入运行。GFFG使用了5年回归期雨量在美国自然资源保护局(Natural Resources Conservation Service，NRCS)三角形网格单元的水文模型上计算得出。Barthold et al(2015)认为FFG及GFFG需要进一步被改进使其更加适合积涝应用，尤其要改变其块状水文模型的局限使其能基于QPF做对流尺度的积涝预报。至此，WRF的水文模块WRF-Hydro着手发展用对流模型输出及其集合预报雨量强迫到分布式水文模型的积涝预报应用(Gochis et al，2015)。

FFG及GFFG的应用现状及发展趋势给用城市水文模型脱机运行方式建立城市地区积涝预警模式提供了借鉴和参考。尽管二者在模型使用、数据基础及预报空间尺度上表现不一，但脱机模式的优势同样能体现在城市积涝预警中。第一，临界雨量阈值的概念和方法更适用于对城市积涝易发点(城市立交桥、十字路口、地势低洼地带等)的积涝监测及预警，要点是通过水文模型的脱机模拟分析可得出这些易发点可能发生各种程度的积水及积涝的雨量阈值。第二，通过脱机方式得出城市地区易发点的对应雷达外推1、3及6 h的临界雨量，便于预报员做出预警选择及判断。除此之外，预报员还可依托不同的预报产品(NWP及融合NWP的预报雨量)，而不仅仅限于短时临近预报产品，灵活地做出城市积涝的预警及预报(赵琳娜等，2012；Poletti et al，2019)。第三，利用城市水文模型的脱机运行方式可不受模型运行时效性及时间快捷性的限制及要求，甚至可不惜耗费巨大的计算时间及资源来搜索临界雨量，因此可配置更高或超高分辨率的水文模型，使用更精致的空间尺度及更短的时间步长来模拟推导临界雨量(Seo et al，2013)，最终满足城市暴雨积涝预警及预报的要求。而且这种运行方式可避免城市水文模型实时运行带来的时间延迟及数据解读过程，有助于提高预报员的判断及预警反应能力。

鉴于上文所论，本研究开发了高时空分辨率城市分布式水文模型，在研究区的城市空间面上模拟1、3及6 h的降水，并强迫到该模型。根据模型模拟出的易发积涝点的积水深度及状况，给出不同预警级别的雨量阈值留做预警指导。在应用中发现其技术路线可行，并对模型模拟及实际应用效果做了初步的检验及分析。当然，这套城市暴雨积涝的预警及预报机制仍然需要在应用及实践中进一步被验证及发展，而从众多的文献中可觑见该模型方法是城市暴雨积涝预警及预报的发展及应用趋势(Ntelekos et al，2006；Gourley and Vieux, 2005；Barthold et al，2015；Herman and Schumacher, 2018)，值得进一步被开发及应用。

1 数据资料 1.1 雷达反演资料

雷达QPE/QPF资料为水文模型强迫的常用雨量数据资源，因其与雨量站数据相比有更好的时空分辨率及空间连续覆盖的优势(van de Beek et al，2010)。研究用北京市气象局瑞图系统(RAMPS-in)的短时临近预报资料做模型强迫。RAMPS-in提供511×581网格的1 km×1 km分辨率的1、3及6 h的QPE/QPF，时间分辨率为10 min。其雨量数据反演自京津冀地区的北京S波段(BJRS)、天津S波段(TJRS)、石家庄S波段(SJZRS)、张北C波段(ZBRC)及承德C波段(CDRC)雷达传送的基数据(陈明轩等，2013)，该数值产品从空间上完整覆盖北京区域。

雷达反演雨量由雷达反射率因子与降水估计之间的Z-R关系换算得到，并非直接观测值。除雷达信号因受遮挡、虚假回波及回波噪音等的干扰，不同降水条件下的雨滴谱分布差异使得雷达降水估计的Z-R关系发生改变，其反演雨量需用雨量站数据做进一步的订正(Anagnostou and Krajewski, 1998)。因此自2012年，北京地区的雷达反演雨量用了近2 000个经数据质量控制的自动站雨量做订正。同时，1、3及6 h的QPFs与高分辨率的非静态数值预报模式(NWP)雨量场作融合(Kober et al，2012)，融合结果减少了原本数据不确定性，提高了这类雨量预报数值产品的质量及可用性(Poletti et al，2019)。

这些QPEs/QPFs(包括融合产品)雨量数据用反距离插值(IDW)从1 km分辨率的Cartesian网格降尺度到30 m的城市水文模型网格，该30 m网格匹配水文模型采用DEM栅格数据，即可完成雨量数据对水文模型的对接或强迫。

城市水文模型模拟分析临界雨量时，是用自定义雨量值为强迫，以遍历方式模拟查证城市积涝易发点的临界雨量值。原理上，暂不需要观测雨量资料作强迫。但是基于实时监测的预报及预警模式仍需用雷达反演及外推雨量做强迫。

1.2 土地利用与土地覆盖

水文模型使用的土地利用和土地覆盖分类数据多从遥感资料中获取(Velásquez et al，2020；Sy et al，2020)。遥感资料无疑是城市水文模型的重要数据源。对于现用的30 m网格高精度土地利用和土地覆盖数据似乎不能完全取自遥感资料。本研究用国家测绘中心提供的2010年30 m×30 m分辨率的全球土地覆盖产品(GlobeLand30)辅助提取城市不透水地表及透水地表类型。GlobeLand30产品主要源自Landsat TM/ETM+影像数据。2010年的GlobeLand产品还用中国环境及灾害卫星(HJ-1)的影像数据辅助生成(Chen J et al，2015)。另外，其他类别的遥感影像也用来验证及补充该分类产品，包括MODIS NDVI、全球测绘数据及DEM数据，以及线上的高分辨率影像，比如GoogleEarth影像、Bing和OpenStreet影像数据。GlobLand30有10种土地分类结果：耕地、森林、草地、灌木、湿地、水体、人工地面及裸地等，可用于提取透水及不透水地表组成。而用于城市水文模型参数输入的30 m格点土地利用及土地覆盖数据集除了必须有透水及不透水地表组成外，还应有地表水动力及与水文参数属性相关的土地利用及土地覆盖分类结果。遥感数据并不单独提供或者完整覆盖及替换此类数据。因此，北京市测绘局提供的2010年的1∶2 000高分辨率GIS基础地理信息数据被用来补充提取所需要的土地利用及土地覆盖数据(图 1)。

1.3 城市地区管网排水能力

城市排水管网的数据对城市水文模型的作用至关重要。GIS城市路网数据可用于提取网格化城市管网排水参数。提取的先决假定是城市排水管网常与路网平行修建，而管网的排水标准可参照GB50318—2017(中华人民共和国住房和城乡建设部，2017), 此标准定义了不同道路等级的管网排水能力(表 1)。依此假定，能根据路网等级推算对应的管网排水能力。

管网排水能力用设计重现期暴雨强度来度量(Zhu et al，2019)。Seo et al(2013)认为汇水接入城市管网系统的城市汇水区的汇流水量曲线变化直接影响排水管网的径流量。城市排水管线系统的排水能力，原则上覆盖环绕路网系统的GIS线型缓冲区域。因此水文模型网格单元的排水能力D_cell用网格单元几何体与管线缓冲区做GIS叠加计算得到, 即

${D_{{\rm{cell}}}} = \frac{{\sum\limits_{i = 1}^n {Inters[B\left(i \right), {O_{{\rm{cell}}}}] \times {D_{{\rm{buffer}}}}\left(i \right)} }}{{{A_{{\rm{cell}}}}}}$

(1)

式中：Inters[B(i)，O_cell]是GIS叠合操作用于计算第i个排水管线缓冲区B(i)与目标网格单元(O_cell)的叠加面积(单位：m²)；n是网格单元内的排水管线缓冲区数量(单位：个)；D_buffer(i)为第i个排水管线的排水能力，常用暴雨重现期(单位：a)来表示该管线所能应对的地表产流量，具体计算参见2.2.1节；A_cell为网格单元面积(单位：m²)。

上述方法依据路网的排水设计标准来提取管网的排水能力，间接获得城市地区网格化管网排水能力。排水能力的网格化参数方案便于二维水文模型的演算，便于模型的功能实现及节省计算量，但在某种程度上会增大模型计算结果的不确定性。尤其是在有人为目标设定的情况下，比如主要积涝点在布设水泵的条件下，会导致模型模拟的结果存在其不能与实际积水状况相吻合的情况。因此，在此建议如果在掌握了城市完整的排水管网及排水设施分布资料的前提下，可依据此类资料提取城市地区网格化管网排水能力指标值。如果只掌握部分“关键性”排水管网资料，建议可在路网资料提取的基础上做适当的修正，减少数据的不确定性。图 2即为提取出的北京城区及近郊区的管网排水能力分布，也基本反映了市中心地区基础设施较完善，管网排水能力稍强的基本状况。

2 城市水文模型

城市水文模型主要包含三个子模型：(1)地表水反应模型用于推算产流量。产流分别在不透水地表及透水地表上生成。不透水地表上的产流生成可简化为降水量减去地表洼蓄和蒸发量而得。透水地表产流则主要为Hortonian渗出(Velásquez et al，2020)。(2)径流路径搜索。(3)水动力模型。水动力模型以2D浅水方程为基础模拟计算积水及积涝状况。

2.1 城市地表产流计算

地表洼蓄主要包含地表吸管蓄水及重力填洼蓄水等(Velásquez et al，2020)，可依据地表土地利用与覆盖类型给出经验值。暖季降水的产流量会因蒸发出现轻微的变化，包括对土壤湿度及洼蓄的影响(Moreno et al，2013)，因此模型对产流的估算考虑地表水汽蒸发。蒸发量(E)的计算引用通用陆面模式(common land model，CLM)的子模块来计算(Meng et al，2009)，就此将水汽蒸散纳入产流的计算。在北京“7·21”特大暴雨个例研究中，CLM估算的降水蒸发量对不透水地表产流大致产生0.24 mm·(24 h)^-1的削减量，而在有植被覆盖的透水地表上的影响则较大，尤其在植被的蒸腾作用下可产生0.5 mm·h^-1的蒸发量。

不透水地表产流R_p(单位：mm)可看作是降水量R(单位：mm)减去洼蓄B_g(单位：mm)及地表蒸发量E(单位：mm)：

${\left\{ {\begin{array}{*{20}{l}} {{R_{\rm{p}}} = R - {B_{\rm{g}}} - E}&{R \ge {B_{\rm{g}}} + E}\\ {{R_{\rm{p}}} = 0}&{R < {B_{\rm{g}}} + E} \end{array}} \right.}$

(2)

对于透水地表的产流R_i(单位：mm)则由Horton模型算出地表下渗量f(单位：mm；Dooge，1992)，然后再减去蒸发及洼蓄：

${\left\{ {\begin{array}{*{20}{l}} {{R_{\rm{i}}} = R - {B_{\rm{g}}} - E - f}&{R \ge {B_{\rm{g}}} + E + f}\\ {{R_{\rm{i}}} = 0}&{R < {B_{\rm{g}}} + E + f} \end{array}} \right.}$

(3)

需要注意的是上述的地表产流为未经由城市排水系统排泄的自然产流。地表产流减去城市排水后，产生地表径流，在排水不良的地势低洼处形成积水及积涝。

2.2 城市排水估计 2.2.1 雨强估计

雨强是估算城市排水能力的一个关键核算指标，一般用Chicago雨量图方法(Chicago hyetograph mthod, CHM)推算回归期内的暴雨雨强值序列而得出(Chen Y et al，2015)。局地雨强值由一定重现期的降雨率q(单位：mm·min^-1) 来表示，用暴雨雨强公式估算(扈海波，2016；Zhu et al，2019)：

$q = \frac{{{A_1}(1 + c\lg P)}}{{{{\left({t + B} \right)}^n}}}$

(4)

式中：t是降雨持续时间(单位：min)，P为暴雨重现期(单位：a)，A₁、B、c及n为区域性暴雨强度特征相关的控制参数(均为无量纲参数，以下如无特别指出的均为无量纲值参数或系数)。A₁为降雨参数(单位：mm)；c是降雨变化参数，b是降雨持续时间纠正参数，n是暴雨消减指数(与回归期相关)。表 2列出北京地区暴雨雨强估算的参数取值(A₁、b、c和n)(扈海波，2016；Zhu et al, 2019), 而全国其他省(自治区、直辖市)城市规划管理部门多根据本地区的暴雨强度特征颁布了不同标准(Zhu et al，2019)。

2.2.2 排水能力推算方法

排水能力指排水管网能有效应付一定重现期暴雨雨强的能力，重现期雨强值取决于地方性或区域性气候特征，其计算公式定义为

$Qs = q' \psi F$

(5)

式中：Qs为城市管网系统的设计排水流量(单位：L·s^-1)，q′为对应的设计雨强(单位：L·s^-1·hm^-2)，ψ是产流系数，F为汇水面积(单位：hm²)。ψ控制地表产流量，其在不同土地利用和土地覆盖类型上是不同的，取值可参见GB 50318-2017(中华人民共和国住房和城乡建设部，2017)。

用式(4)和式(5)可推导出计算管网排水能力的设计雨量的公式：

$q\prime = \frac{{{A_1}(1 + c\lg P)}}{{\psi {{\left({t + B} \right)}^n}}}$

(6)

这样，城市地区单位时间产流r_d(单位：mm·min^-1)等于其透水地表及不透水地表上单位时间产流r_s(单位：mm·min^-1，其中r_s=R_p+R_i)减去管网设计排水量q′，公式表述为：

${r_{\rm{d}}} = {r_{\rm{s}}} - \frac{{{A_1}(1 + c\lg P)}}{{\psi {{\left({t + B} \right)}^n}}}$

(7)

2.3 城市地表径流路径选向

地表径流路径的搜索常用基于滑动窗口的D8算法(图 3a)，该算法在DEM中以单个高程格点八个方向的最大下降坡度作为径流流向，以格点中心点的连线作为径流路径(Jenson and Domingue, 1988)。Gallant and Hutchinson(2011)认为当径流紊流和分流效应导致了径流分散，坡度线和径流线在水文地形形态上是不一样的。在此情况下，多重径流流向的寻径才是描述径流分散的关键及首选(Seibert and McGlynn, 2007)。在算法实现上，多重径流流向的判断仍然可依据格点与周边相邻四个格点的坡度来确定，即水流从四个方向流向水位较低的格点(图 3b)。众多文献提到的LISFLOOD-FP模型就用这种四朝向分流算法(Bates and De Roo，2000；Zhu et al，2019)。此算法简便，格点间产流的流出和流入依据与周边格点的坡度来判断。而径流交换加快，格点间的水位变化会加快，这就要求水动力模型Manning水流的时间分辨率较高，不然容易产生径流不连线性及偏差过大的情况(Hodges，2019)。D8算法和四流向算法的坡度线和对应的径流网络均可直接在DEM等高程数据的等高线上搜索得到(Moretti and Orlandini, 2008)。除此之外，在水动力模型中加入回水机制，同样可达到多重径流的效果。Orlandini et al(2012)指出多重径流算法所产生的径流分流更多地受格点单元大小的限制，在使用该算法搜索径流分散时需要反复对其检验及验证。因此，Orlandini et al(2012)提出D8的改进算法(D8-LTD)，该算法通过径流路径规整及规划，完善及优化了径流路径的生成，这样D8-LTD算法除了保有路径明确(确定及唯一性)、无分散、结果可靠及计算效率较高的优点外，还能近似在任何空间尺度上保持坡度线和径流路径的一致(Moretti and Orlandini, 2008)。

尽管如此，本文认为在城市环境中，即便在地势起伏不大而城市建筑等阻碍物较多的情况下，城市地表径流的分散及紊流情况也会比较明显及复杂。因此四朝向分流算法(图 3b)及2.4节介绍的具有回水计算过程的水动力模型(图 3c)比单路径的D8算法更适用于城市水文模型的地表径流路径搜寻及水动力模型的计算。当然，为保证径流的连续性特征，需要规划水动力模型的时间步长。

2.4 城市水动力模型

城市水动力模型是城市水文模型的核心。此模型建于积水深度均一的二维浅水方程上。积水过程的浅水假设会降低模拟的准确性，增加模拟结果的不确定性和可接受程度。但是城市水文模型的模拟精度更多受数据不确定性的影响(比如地形参数和地表粗糙度等边界条件的制约)，而动力过程的模型简化对此的影响有限。这种非惯性及扩散波动模型可加入回水效应，实现多径流寻向的效果。城市积涝的水流模式其实是介于分散波动和全动力方程之间。Bates and De Roo(2000)也认为城市积水的水动力过程可简化为消除了局地加速、对流加速及对流压力项的一维St.Venant方程，甚至假定摩擦与重力驱动是平衡的。

用浅水方程作为城市水动力模型基础，其表述式如下：

连续方程为：

$\frac{{\partial h}}{{\partial t}} + d\left({\frac{{\partial u}}{{\partial x}} + \frac{{\partial v}}{{\partial y}}} \right) = q$

(8)

动力平衡方程为(x方向和y方向)：

$\frac{{\partial u}}{{\partial t}} + g\frac{{\partial h}}{{\partial x}} + g{S_{{\rm{f}}x}} = 0$

(9)

$\frac{{\partial v}}{{\partial t}} + g\frac{{\partial h}}{{\partial y}} + g{S_{{\rm{f}}y}} = 0$

(10)

式中：h是水位(水深加相对高程，单位：mm)；d是径流深度(单位：mm)；t是径流时间(单位：s)；u和v分别为在x方向和y方向上的径流速度(单位：m·s^-1)；S_fx和S_fy分别为在x方向及y方向上的摩擦坡度，g为重力加速度。摩擦项通过Manning公式计算，即：

${S_{{\rm{f}}x}} = \frac{{{n^2}u\sqrt {{u^2} + {v^2}} }}{{{d^{\frac{4}{3}}}}}$

(11)

${S_{{\rm{f}}y}} = \frac{{{n^2}v\sqrt {{u^2} + {v^2}} }}{{{d^{\frac{4}{3}}}}}$

(12)

式中：n为Manning系数。

在式(8)中q为潜层水流项，在城市水动力建模中假定为0。Seo et al(2013)认为在城市环境中透水地面被周边不透水地面(水泥陆面及建筑物)阻断的情况下，其部分下渗水流仍然会潜流到城市排水系统中，影响排水系统中汇流曲线10%的流量，因此提出了在城市排水模型系统中针对城市潜层水流的计算方法。但是对城市潜层水的估算非常复杂，即便出流到排水系统中，也对地表产流、积水及水淹状况的影响不大。因此，此处仍然设定q值为0。

需要指出的是基于D8算法及四流向算法的水动模型，以单一格点为基本计算单元来核算其进水及出水量，而且假定摩擦与重力驱动是平衡的，再用Manning公式直接计算格点的进水及出水的流速(Bates and De Roo，2000)，那样确实是简化了计算。但在计算漫流过程中，比如用四流向算法，在单步执行时，随着时间的推移，就会出现径流交换加快，而只能限制模型的时间步长来减少计算偏差。当然时间步长越短，偏差就越少。因此，这里还是认为用有限差分的方法来计算水动力方程可暂时解决这个难题。

将式(8)改写为以下离散方程：

$\begin{array}{l} h_{i, j}^{n + 1} = \Delta t\left[ {\left({\frac{{{d_{i - 1, j}} + {d_{i, j}}}}{{2\Delta x}}} \right)u_{i - 1, j}^{n + 1} - } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\left({\frac{{{d_{i, j}} + {d_{i + 1, j}}}}{{2\Delta x}}} \right)u_{i, j}^{n + 1} + \left({\frac{{{d_{i, j - 1}} + {d_{i, j}}}}{{2\Delta y}}} \right)v_{i - 1, j}^{n + 1} - \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\left({\frac{{{d_{i, j}} + {d_{i, j + 1}}}}{{2\Delta y}}} \right)v_{i, j}^{n + 1}} \right] + q\Delta t + h_{i, j}^n \end{array}$

(13)

x和y方向上的水动力平衡方程[式(9)、式(10)]分别可改写为

$u_{i, j}^{n + 1} = - g\frac{{\Delta t}}{{\Delta x}}\left({h_{i + 1, j}^{n + 1} - h_{i, j}^{n + 1}} \right) - u_{i, j}^n\left({1 - {S_x}} \right)$

(14)

$v_{i, j}^{n + 1} = - g\frac{{\Delta t}}{{\Delta y}}\left({h_{i + 1, j}^{n + 1} - h_{i, j}^{n + 1}} \right) - v_{i, j}^n\left({1 - {S_y}} \right)$

(15)

将式(14)和式(15)简化：gΔtS_fx=uS_x，gΔtS_fy=vS_y。

然后，将式(14)、式(15)分别代入式(13)整理后，可以得到类似下面的简写式

$ah_{i, j}^{n + 1} + bh_{i + 1, j}^{n + 1} + ch_{i, j + 1}^{n + 1} + dh_{i - 1, j}^{n + 1} + eh_{i, j - 1}^{n + 1} = fh_{i, j}^n$

(16)

式中：h_i，jⁿ⁺¹为格点(i, j)处第n+1时间步的水位值，h_i，jⁿ为对应格点第n时间步的水位值，为已知量，其余水位值变量的定义与此类似。a、b、c、d、e和f系数均为已知值，可由式(13)~(15)的系数变量推导得出，具体可自行推导，此处就不赘述。

由此可见式(16)中，有5个未知量(h_i-1，jⁿ⁺¹, h_i，j-1ⁿ⁺¹, h_i，jⁿ⁺¹, h_i+1，jⁿ⁺¹, h_i，j+1ⁿ⁺¹)，而直接求解式(16)的难度很大，需要在水文模型(I × J)网格上建立大小为(I +1)(J +1)的二维矩阵，去解类似Ax= b的矩阵方程组，而不管是用矩阵的LU分解或Cholesky分解, 其计算及存储量都超级大。因此Wang and Chen(2003)提出了用变向隐式法(alternating-direction-implicit，ADI)(图 3c)来计算类似于式(16)的差分方程式，也就是半步差分法。ADI将二(三)维计算难题简化为二(三)个连续的一维计算过程，相关算法参见Wang and Chen(2003)，该文献还给出了Peaceman-Rachford算法及Douglas-Gunn算法的两套计算方案。附录为ADI的计算结果。需要提及的是这些算法最后均用Thomas算法来求解二维对角矩阵(https://www.sciencedirect.com/topics/computer-science/thomas-algorithm)，而Thomas算法的时间执行度仅为O(n)。

2.5 水动力模型时间步长的限定

用有限差分方法来计算格点的水位值，在浅水方程中水位的调整在x和y方向上进行，这种调整能够体现径流的回水效应，因此多重径流方向的选取与四流向算法一样。但四流向算法是直接用Manning公式来计算格点的进出水情况，就如同之前讨论的那样，必须限定水动模型的时间步长，而且由时间延迟产生的计算偏差随时存在，延迟愈大偏差就越大。而有限差分方法虽然不会因时间步长的存在而产生偏差，但如果时间步长过大，会导致迭代计算失效，也就是这种ADI会导致计算不稳定(Caviglia and Dragani, 1996)。因此，水动力模型的时间步长，需由式(17)作限定条件

$\left({\frac{{|u| + \sqrt {gd} }}{{\Delta x}} + \frac{{|v| + \sqrt {gd} }}{{\Delta y}}} \right)\Delta t\leqslant 1$

(17)

式中：Δt为模型时间步长(单位：s)，Δx、Δy为在x和y方向上的格距，d为水深。

3 基于城市水文模型的暴雨风险预警实例应用及研究分析

城市水文模型的实例应用分为两部分。一是模拟雷达反演雨量强迫到水文模型的在线运行方式。该方式依据实况运行得到的积水深度及水淹面积来做出积涝发生及发展的判断及预警。二是用水文模型的模拟来搜索城市地区积涝易发点的积涝风险阈值或一定预警等级的临界雨量。

3.1 以北京“7·21”特大暴雨为个例的城市水文模型应用及效果检验

城市水文模型网格格点数为1 521×1 908，分辨率为30 m，网格左上角经纬度坐标为40.104558°N、116.554678°E。其图幅范围覆盖北京主要城市区域。模型积分的时间步长为2 s。在30 m网格大小情况下，完全能满足式(17)对水动力模型时间步长的限定。

城市水文模型以2012年7月21日北京地区特大暴雨为个例，用QPE做强迫，模拟地表产流、管网排水量、径流量、积水淹没深度及范围等参量。模型从“7·21”当天的11时模拟到22时。当模拟到21时时，结果显示整个北京城市区域呈现大面积城市积涝情景(图 4)。此情景与“7·21”当天的实际发生情况几乎一致。而且积水严重地点多集中在道路主干道、道路交叉口及立交桥等地势低洼处(对照图 4中的积水图斑)。检查12个积涝易发点[广安门、公主坟及西客站隧道等(图略)]，这些积涝点的平均最大积水深度可达2.22 m左右，与当天的积涝情景吻合。从叠加上地形图的水淹分布图(图略)上也可看到积水区域主要分布在地势低洼处。这里仅列出广安门、莲花桥和六里桥的积水模拟结果及与从网上搜集到的“7·21”的实况灾情图片做对照组合展示，结果显示水文模型的模拟结果和实际情况比较贴近(图 4)。

3.2 基于城市水文模型模拟的北京城市积涝易发点暴雨积涝风险阈值搜索及检验分析

为找到城市暴雨积涝风险预警的临界雨量，用城市水文模型脱线模拟的方式搜索临界雨量。方法用初始的1、3及6 h雨量[比如10 mm·(1 h)^-1、10 mm·(3 h)^-1及10 mm·(6 h)^-1雨量]强迫到城市水文模型，然后强迫雨量值按一定增量(1 mm)递增，循环模拟及搜索在不同雨量强迫下，城市地区有可能出现0.2、0.5、0.8、1和1.2 m等不同积水深度情况下的临界雨量值。由于以1 mm为增量递增，因此雨量分辨率为1 mm。研究统计了不同雨量强迫条件下0.2、0.5、0.8和1.2 m等积水深度的格点面积占总格点面积的占比，其1、3及6 h的模拟结果情况见图 5。研究发现在1 h雨量不大于12 mm，3 h雨量不大于19 mm及6 h雨量不大于31 mm的情况下，研究区图幅范围内的积水都没有大于0.2 m深度的格点。如果以0.2、0.5、0.8和1.2 m积水深度的格点面积占总格点面积比例刚好达到0.5%作为城市地区积涝灾害的蓝、黄、橙和红色预警划分标准，那么1 h预警等级临界雨量分别是24、37、45和52 mm，3 h临界雨量是30、40、48和62 mm，6 h临界雨量则是39、48、55和69 mm。

对照Gourley and Vieux(2005)对GFFG的统计分析可发现红色预警的1 h临界雨量超过了美国二千多个流域的1 h的暴雨洪涝风险雨量阈值的众数(55 mm左右)，而橙色预警的1 h临界雨量大致在GFFG的中等偏下位置；3 h临界雨量情况与1 h的类似；而6 h红色预警临界雨量值略等于GFFG的众数值。由于选取的北京城区范围的不透水地表面积比例较大而且集中，地势也较为平坦，属于积涝较敏感地区(扈海波等，2013；尹志聪等，2015)，因此得到这个结果也是可解释的。

根据2007年北京市防汛办提供的资料显示，北京地区积涝易发点共计49个(图 6)，从图中可发现这些积涝易发点也主要分布在城市交通主干道、道路十字路口及立交桥下。研究在模型模拟结果的图幅范围内搜索这些积涝易发点附近0~300 m范围内，当积水深度超过0.2、0.5、0.8和1.2 m等的临界雨量值。搜索得到的积涝易发点临界雨量结果见图 7。这里仍然以积水深度分别超过0.2、0.5、0.8和1.2 m的临界雨量作为易发点“蓝、黄、橙、红”暴雨积涝风险预警等级指标值。当然，预警指标明确地按分级划分临界雨量虽然迎合了确定性预报的要求。但此类指标的不确定性依然存在，例如土地利用及土地覆盖类型改变、城市排水能力的改进(排水泵的增加及管线的改造)及周边地形地物的变化等都会导致风险阈值发生改变。因此，风险阈值也需要在水文模型基本参数及输入更新的情况下，随时更新及修补。另外，基于临界雨量的风险预警还应适当运用不确定性的预警及预报方法。Barthold et al(2015)认为用风暴边缘区域的75%和90%的FFG依然能够有效识别及预知暴雨积涝。因此在实际应用中还应灵活地参照临界雨量，考虑不确定性因素，适量滑动临界雨量取值。

城市积涝易发点临界雨量值开发应用的一个难点是如何有效地检验和验证临界雨量的有效性及合理性。包括城市水文模型模拟结果的检验和验证一直是个国际难题(Welles et al, 2007；Amponsah et al, 2016；Sy et al, 2020)。Amponsah et al(2016)认为城市积涝事件及其水文过程的对照观测非常稀少，因为积涝事件的时空尺度很小，在受到此类事件影响的地区缺少相对应的气象水文观测，尤其没有产流强度及水动力过程的观测，因此针对此类事件的气象水文观测验证是几乎不可能实现的。Sy et al(2020)则提出通过实地收集和考察来补充及验证历史积涝灾情信息，认为城市居民对灾情的反馈是积涝分析的重要信息源，甚至把由此可获得的分析结果定义为“市民科学”。因此关于临界雨量结果的正确性及合理性除了需要与国内外的文献结果做下纵向的对比外，还需要在实际应用中不断检验和验证。尤其是“市民科学”的积极参与及投入。

本研究在完成城市积涝易发点不同风险预警等级下的临界雨量值划分后，开展对应的应用个例检验。这里以2020年7月12日发生在北京南部城区的局地暴雨为例来说明个例在临界雨量中的检验和应用情况(图 8)。该暴雨积涝事件发生在北京市丰台区南沙窝桥及岳各庄桥一线，从21时左右的该地路面的积涝实情(图略)来看积水已比较严重，积水深度至少在0.5 m左右或以上。对照易发点的暴雨风险预警临界雨量(图 8a)，20时的RMAPS-in的(图 8b)1 h的QPE达到25 mm左右，这个1 h临界雨量值对应的预警等级应达到橙色预警程度。如前文所述，即便用QPE做暴雨积涝风险的预警仍然有一定的时间提前量。而从19时起报的QPF来看其1 h雨量为15~20 mm(图 8c)，这个雨量值能对应到蓝色及黄色预警等级。可见QPF用来做预警的效果大多受限制于预报结果的准确性，包括落区及雨量准确性。但QPF在预警服务中仍受欢迎，因其具有一定的时间提前量，这个在预警中也是一个关键指标，具无法替代的优势和作用(Silvestro et al, 2017；Thorndahl et al, 2017；Poletti et al, 2019)。

4 结论与讨论

城市暴雨积涝风险预警及预报需要结合气象及水文两个领域的相关模型进行交叉应用及开发。本文重点论及城市水文模型基本原理及模型方法实现，以及雷达反演雨量强迫到城市水文模型的积水及积涝情景模拟。然后给出基于城市水文模型的城市暴雨积涝风险预警及预报技术应用。

从模型应用的情况来看，雷达反演雨量的定量降水估计QPE在雨量值及落区准确性上基本能满足城市水文模型强迫的需求。由于地表产流及径流生成的延迟效应，一般径流峰值会落后降水峰值近1 h左右，因此可充分利用这个时间差来发挥QPE在积涝预警中的作用。也就是说结合QPE和临界雨量值可做暴雨积涝风险预警，其时间提前量可控制在1 h范围内。而定量降水预报(QPF)的雨量值及落区有较大偏差，尤其是3 h及6 h的预报。此类数据直接强迫到城市水文模型会产生一定的不确定性模拟结果。但是其具相当的时间提前量，同样是值得积涝预报使用的重要数据，而其1 h预报产品从初步使用的情况来看，基本上还是能根据其雨量场结果做出有效的预警判断。今后工作中，更应注重结合QPFs及临界雨量判断的暴雨积涝预警应用。当然，也更加期待QPFs的1、3及6 h数据准确性和有效性能进一步完善，真正达到城市暴雨积涝预报的要求。

土地利用及土地覆盖是城市水文模型的重要参数，该参数控制及贯穿应用于城市水文模型的地表水反应过程及水动力过程，因此其数据的准确性及分辨率等均影响模型模拟结果的可靠性。除不透水地表组成是地表产流分析的重要基础外，地表土地类型属性及分类更是控制水动力模型参数的基础，比如水动力模型中的地表粗糙度及产流系数等关键参数。此数据的不确定性很大程度制约了水文模型的模拟精度。城市水文模型的高空间分辨率(30 m)特性对土地利用及土地覆盖的分类结果及精度要求也较高。这类资料最好从高精度的基础地理信息数据中提取，再辅以高时空分辨率的遥测资料做分类订正。高分辨率土地利用及土地覆盖等参数在城市水文模型的模拟分析中的优势体现在有助于获得较为准确的地表产流(流量及流速)及汇流，从而能提高模型模拟的整体精确度。

城市水文模型的地表水反应过程的模拟主要在不透水地表及透水地表上展开。在产流估算上，不透水地表的产流损失分量只有少量的洼蓄及蒸发，而透水地表的产流还要持续性减去Hotanian下渗量，因此城市不透水地表面积越大，地表产流越多，积涝风险就越大。另外，用路网GIS资料间接推算的城市管网的排水能力可参数化到城市水文模型中，其导致模拟的不确定性及有效性需要进一步评测，但这种参数化方案在计算资源分配及可参数化应用上有一定的优势。提取到的城市排水能力参数分布表明城区的管网排水能力大于近郊区，但不透水面积比例较大，因此反而更容易发生积涝事件。

需要特别提及的是在径流路径寻径上，城市地区的地表建筑物等相关阻碍物较多，地表径流的紊流、分流及散流情况也多，因此建议城市水文模型的径流路径寻向用多流向算法或通过与多流向算法效果类似的回水机制来实现。

城市水动力模型可建立在基于积水深度均匀的二维浅水方程上。浅水假定消除了局地加速、对流加速及对流压力项，但其对模拟精度的影响不大，甚至远低于地形参数和地表粗糙度等边界条件不确定性对模拟精度的影响。浅水下的低速缓流状态完全符合城市地表径流及积水汇流的水动力情况。浅水假定还便于城市水动力方程的求解。本文提出借用ADI分两步求解水动力模型的差分方程。而且水动力模型差分方程的求解包含了回水机制，同时对模型的时间步长做了规约，确保模型格点单元的水量连续性及均衡性(Hodges, 2019)。需要注意的是尽管Thomas算法只是个小的计算技巧，该算法却是求解差分方程的关键。

本文论及的暴雨积涝风险预警及预报模式主要有在线及脱线运行两种方式。从实例应用及分析来看，北京“7·21”特大暴雨个例的实况雷达反演雨量强迫到水文模型，大体能模拟到当天的积涝情景。尽管用城市水文模型实时模拟积水状况并用模拟结果实施暴雨积涝风险预警在技术上可行，但水文模型在线方式在气象业务运行过程中会因模型模拟参数读取及结果解读的繁杂处理流程、较高的数据处理及计算量、结果的不确定性、模型运用的熟练度及运用技巧限制等缺陷而使得其应用前景并非看上去那么美好。因此，本文提出参照NWS的FGGs的脱线模式不失为一种更为合理的城市暴雨积涝预警及预报模式。基于城市水文模型的脱线运行方式把积涝易发点对应一定积水深度的临界雨量模拟输出，此临界雨量可作为预报员在做暴雨预警时的重要参考，方便预报员直观做出合理的预警判断及发布。研究以0.2、0.5、0.8和1.2 m积水深度的格点面积占总格点面积比例达到0.5%作为城市地区积涝灾害的蓝、黄、橙及红色预警划分标准，得出北京地区1 h的对应预警等级的临界雨量分别为24、37、45和52 mm；3 h临界雨量则分别为30、40、48和62 mm，6 h分别为39、48、55和69 mm。然后，用城市水文模型模拟出了北京地区至少49个积涝易发点在积水深度达到0.2、0.5、0.8和1.2 m时的1、3和6 h的临界雨量，并以这些临界雨量作为积涝易发点暴雨积涝的“蓝、黄、橙、红”预警等级判断及实时监控标准。通过与国内外文献资料的对比分析及个别积涝个例的初步检验，基本可认为这些临界雨量及对应的预警等级标准有一定的可信度及合理性。

附录：

算法需要将从时间步长n到n+1步分解成两个子步长来完成：从n到$ n + \frac{1}{2}$步及从$ n + \frac{1}{2}$到n+1步。第一步中，在x方向上的隐式计算(implicit)式如下，x方向的计算式可组合为一个线性方程：

${a_x}h_{i - 1, j}^{n + \frac{1}{2}} + {b_x}h_{i, j}^{n + \frac{1}{2}} + {c_x}h_{i + 1, j}^{n + \frac{1}{2}} = {d_x}$

(1)

其中：

$ \begin{array}{l} \;\;\;\;\;\;\;\;{a_x} = - \frac{{g\Delta {t^2}}}{{2\Delta {x^2}}}\left({{d_{i - 1, j}} + {d_{i, j}}} \right)\\ {b_x} = g\Delta {t^2}\left[ {\left({\frac{{{d_{i - 1, j}} + 2{d_{i, j}} + {d_{i + 1, j}}}}{{2\Delta {x^2}}}} \right) + } \right.\\ \;\;\;\;\;\;\;\;\left. {\left({\frac{{{d_{i, j - 1}} + 2{d_{i, j}} + {d_{i, j + 1}}}}{{2\Delta {y^2}}}} \right)} \right] + 1\\ \;\;\;\;\;\;\;\;{c_x} = - \frac{{g\Delta {t^2}}}{{2\Delta {x^2}}}\left({{d_{i, j}} + {d_{i + 1, j}}} \right)\\ {d_x} = \Delta t\left[ {\left({\frac{{{d_{i - 1, j}} + {d_{i, j}}}}{{2\Delta x}}} \right)u_{i - 1, j}^n\left({1 - {S_x}} \right) - } \right.\\ \;\;\;\;\;\;\;\;\left({\frac{{{d_{i, j}} + {d_{i + 1, j}}}}{{2\Delta x}}} \right)u_{i, j}^n\left({1 - {S_x}} \right) + \\ \;\;\;\;\;\;\;\;\frac{{g\Delta t}}{{2\Delta {y^2}}}\left({{d_{i, j}} + {d_{i, j - 1}}} \right)h_{i, j - 1}^n + \\ \;\;\;\;\;\;\;\;\left({\frac{{{d_{i, j}} + {d_{i, j - 1}}}}{{2\Delta y}}} \right)v_{i, j - 1}^n\left({1 - {S_y}} \right) + \\ \;\;\;\;\;\;\;\;\frac{{g\Delta t}}{{2\Delta {y^2}}}\left({{d_{i, j + 1}} + {d_{i, j}}} \right)h_{i, j + 1}^n - \\ \;\;\;\;\;\;\;\;\left. {\left({\frac{{{d_{i, j + 1}} + {d_{i, j}}}}{{2\Delta y}}} \right)v_{i, j}^n\left({1 - {S_y}} \right)} \right] + \\ \;\;\;\;\;\;\;\;h_{i, j}^n + q\Delta t \end{array} $

第二步在y方向上的显式(explicit)计算式如下, 同样可组合为线性方程组:

${a_y}h_{i,j - 1}^{n + 1} + {b_y}h_{i,j}^{n + 1} + {c_y}h_{i,j + 1}^{n + 1} = {d_y}$

(2)

$ \begin{array}{l} \;\;\;\;\;\;\;\;{a_y} = - \frac{{g\Delta {t^2}}}{{2\Delta {y^2}}}\left( {{d_{i,j - 1}} + {d_{i,j}}} \right)\\ {b_y} = g\Delta {t^2}\left[ {\left( {\frac{{{d_{i - 1,j}} + 2{d_{i,j}} + {d_{i + 1,j}}}}{{2\Delta {x^2}}}} \right) + } \right.\\ \left. {\;\;\;\;\;\;\;\;\left( {\frac{{{d_{i,j - 1}} - 2{d_{i,j}} + {d_{i,j + 1}}}}{{2\Delta {y^2}}}} \right)} \right] + 1\\ \;\;\;\;\;\;\;\;{c_y} = - \frac{{g\Delta {t^2}}}{{2\Delta {y^2}}}\left( {{d_{i,j + 1}} + {d_{i,j}}} \right)\\ {d_y} = \Delta t\left[ {\frac{{g\Delta t}}{{2\Delta {x^2}}}\left( {{d_{i - 1,j}} + {d_{i,j}}} \right)h_{i - 1,j}^{n + \frac{1}{2}} + } \right.\\ \;\;\;\;\;\;\;\;\left( {\frac{{{d_{i - 1,j}} + {d_{i,j}}}}{{2\Delta x}}} \right)u_{i - 1,j}^n\left( {1 - {S_x}} \right) + \end{array} $

$ \begin{array}{l} \frac{{g\Delta t}}{{2\Delta {x^2}}}\left({{d_{i, j}} + {d_{i + 1, j}}} \right)h_{i + {1},j}^{n + {\frac{1}{2}}} - \\ \left({\frac{{{d_{i, j}} + {d_{i + 1, j}}}}{{2\Delta x}}} \right)u_{i, j}^n\left({1 - {S_x}} \right) + \\ \left({\frac{{{d_{i, j}} + {d_{i, j - 1}}}}{{2\Delta y}}} \right)v_{i, j - 1}^n\left({1 - {S_y}} \right) - \\ \left. {\left({\frac{{{d_{i, j + 1}} + {d_{i, j}}}}{{2\Delta y}}} \right)v_{i, j}^n\left({1 - {S_y}} \right)} \right] + \\ h_{i, j}^{n + \frac{1}{2}} + q \end{array} $

在式(1)和式(2)中，每一个子式中均含有$ \left( {h_{i - 1,j}^{n + \frac{1}{2}},h_{i,j}^{n + \frac{1}{2}},h_{i + 1,j}^{n + \frac{1}{2}}} \right)$以及(h_i-1，jⁿ⁺¹，h_i，jⁿ⁺¹, h_i+1，jⁿ⁺¹)两组未知值，该未知值的计算可由两组已知系数值(a, b, c)，分别在x方向及y方向上，形成一个对角三角形矩阵方程组。