郑州“7·20”暴雨移植子牙河流域方法及成果

□杨丰源

2021年汛期，郑州发生了历史罕见的“7·20”特大暴雨，造成了重大人员伤亡和财产损失。许多地区、部门以期通过郑州暴雨在本地区重演，来检验防范措施是否有效，应急预案是否完备可行，达到补齐短板、降低损失的目的。现以子牙河流域为例，从天气形势上分析目标流域发生郑州“7·20”暴雨的可能性，然后分析暴雨中心位置，进而介绍移植技术、移植结果，为子牙河流域防范郑州“7·20”级别特大暴雨提供参考。

1.移植可能性分析

历 史上1939年、1956年、1963年、1996年暴雨，均在子牙河系出现大洪水，天气形势与此次相近。从形成原因来看，1939年夏天赤道附近系统活跃，副高位置偏西北，台风北上次数多，海河流域汛期暴雨次数多；1956年受台风外围及副高影响，输送充足水汽；1963年暴雨成因为多次低压槽和低涡共同影响结果；1996年暴雨成因为9608号台风北上减弱为低气压与西风槽和冷空气结合的影响。“7·20”暴雨为集合南亚高压东伸、副热带高压西伸北抬、环流形势阻塞、低涡缓慢西移、强台风远程输送水汽、低空急流发展等多种形势和因素共同影响。因此，从天气形势上分析，子牙河系存在发生“7·20”暴雨的可能性；从地形特征来看，子牙河系西部太行山迎风坡抬升特性明显，“7·20”暴雨郑州市西北侧山地起到了相同的作用，因此将郑州“7·20”暴雨移植至子牙河系是可行的。

2.暴雨中心选择

结合典型历史暴雨洪水的暴雨中心位置，如1939年滏阳河一带，1956年和1963年暴雨中心分别发生在滹沱河上游井陉县和滏阳河獐犭厶，1996年暴雨中心在野沟门水库，2016年暴雨中心分别发生在邢台临城一带、石家庄赞皇、石家庄井陉一带，可以看出暴雨中心基本分布在太行山迎风坡，滹沱河主要以冶河中上游为中心，滏阳河主要以临城水库、朱庄水库等为中心。

综上分析，若将“7·20”郑州暴雨移植至子牙河流域，那么暴雨中心应在滏阳河朱庄、临城水库上游或者滹沱河冶河中上游，现以暴雨中心在冶河微水上游和泜河临城水库上游进行移植分析。

3.移植技术应用

基于数字高程，提取流域边界，使用反距离权重内插将平移后点雨量值内插至网格，再计算流域面雨量。

3.1 填洼处理

在依据数字高程圈画流域中，流域高低起伏。计算水流方向时，存在有些水不能流出洼地，使得水系产生误差或不连续。为了使水系在遇到洼地流向依然连续，需对DEM数据进行洼地处理。填洼前基于数字化流域电子水系，采用AGREE算法对原始DEM数据进行预处理，使DEM和输入的矢量图层相一致。先设置缓冲区、平滑、增益3个参数。河道及周围区域地形呈“V”形，以确保水最终流入河流网格。根据设置的增益参数，河道所在网格高程将再次降低相应值，以便根据矢量河道校正DEM。AGREE算法修正后流域数字高程模型（DEM）GTOPO30见图1。

图1 AGREE算法修正后流域数字高程模型（DEM）GTOPO30

填洼处理的目的是将洼地内部高程填至与出口高程一致，按照增加中心栅格高程方法填平洼地，使处理后的DEM能生成收敛连续的河道。

3.2 网格流向计算

在填洼后DEM基础上计算网格流向，使用D8单流向算法：根据公式（1）计算每个网格单元与周围网格的坡度，其中为中心网格与相邻网格中心之间的高程差。按最陡坡度原则设定该单元水流流向，流向用一系列数字代表。D8算法提取流向见图2。

图2 D8算法提取流向

3.3 集水面积确定

在确定流向后，可根据流向分析集水面积。若上游栅格的水量都汇集到某一栅格点，便可据此计算该栅格的累计水量。流域出口处，集水面积最大，流域内的水最终都流向它。根据图2中DEM确定的集水面积矩阵见图3。

图3 根据图2中DEM确定的集水面积矩阵

3.4 河道阈值确定及流域提取

划清河道及水系需搞清楚形成永久性水道所必须的集水面积，集水面积大于等于集水面积阈值的栅格被定义为河道，反之则为坡地。有些研究认为绝大多数侵蚀性斜坡上凸下凹，并且在两部分间有一个转折，转折出现的临界区域确定了河网的界限。在具有足够侵蚀介质的区域（例如具有良好土壤覆盖的流域），侵蚀能力（输沙率）f可表示为：

式中：

s—坡度；

q—流量；

k—侵蚀扩散率；

m、n—系数。

坡地的侵蚀主要由坡度决定，河道的侵蚀除坡度外也与流量关系密切。对于坡地n=0、m=1，对于河道n通常大于1、m大致为2。Smith和Bretherton进一步研究认为形成凹形坡的条件为：

反之，则容易形成一个凸斜率。对于流量Q和流域A，通常满足Q∝Ax，指标x与流量回归周期有关。当考虑较大的回归周期(如年流量)时，x=1。Krikby即以单位等高线的集水面积a取代式（2）和（3）中q，得到：

式中：

F—输沙量；

β、u、w—系数。

(5)式建立时，容易形成凹坡，反之则容易形成凸坡。经过进一步的研究，Tarboton得出结论，在当地条件下有一个固定的输沙速率。在此假设下，输沙量只与流域面积有关:

U为常数，（6）式两边分别对a求导并都乘以a得到公式（7）进而得到公式（8）：

（8）式右边与（5）式形式相近，因而公式（9）成立时易形成凹形坡，a(∂F/∂s)为正值，故凹形坡满足ds/da<0；反之，凸形坡则为ds/da>0。

基于上述理论，阈值通常由集水区面积和河道平均坡度之间的关系得出。首先给出集水区的取值范围，计算每个数字排水系统平均坡度，绘制平均坡度与流域面积关系曲线。平均坡度与流域面积关系曲线拐点对应的流域面积，可视为流域内河流地貌发育的关键支撑区域，即流域面积阈值。

通过计算分析，结合数字化矢量水系，推求出流域集水面积阈值为50栅格。确定集水面积阈值后，集水面积大于等于50的栅格被定义为河道，其余均为坡地。海河流域数字水系见图4。

图4 海河流域数字水系

3.5 流域确定

确定分水岭边界的算法如下：确定阈值后，搜索DEM矩阵以找到分水岭面积大于或等于阈值的网格，网格是流域的出口（DEM矩阵可以包含多个流域），按顺序编号（n=1,2,3，…）。对于任何编号的出口网格，需根据流量矩阵表向上搜索，任何流向它的网格都分配给这个一编号的网格，并继续循环搜索，直到再也没有流向该网格为止。水系矩阵（a）与流域矩阵（b）见图5。

图5 水系矩阵（a）与流域矩阵（b）

3.6 控制站分区

控制站分区一般按水文站控制，先划分最下游站控制的集水区，再依次划分上游的集水区，这样上游区会覆盖在原来的下游区中。各区域在生成过程中，以流量控制站的栅格为该集水区成员，标上区号，不断循环直到所有栅格被标上区号。

3.7 降水空间分布及控制站分区面雨量计算

空间分布，也称为空间插值。根据是否考虑高程影响，降水的空间分布可分为仅考虑平面位置影响的二维空间插值、仅考虑高程影响的降雨—高程线性回归和同时考虑平面位置和高程影响的三维空间插值。考虑到数据的局限性，此文采用二维空间插值，只考虑插值。

泰森多边形网格和泰森多边形雨量测量站之间的距离是流域内所有雨量测量站中最短的。然而，根据泰森多边形划分降水量将导致降水量不均匀。在泰森多边形的内部，降雨是一致的，泰森多边形之间的降雨可能会存在较大差异。

与泰森多边形法不同，反距离加权平均法在降雨连续性方面更好。思路就是网格降雨与雨量站所在的网格降雨成正比，与从网格到每个雨量测量站所在网格的距离成反比。计算公式如下：

式(10)（11）中：

P—期望格网的降雨值；

N—雨测站个数；

P（si）—第i个雨测站的降雨值

λ1—第i个雨测站的权重

d1—第i个雨测站与期望格网的距离；

b—权 重指标，b=0时，式(10)演化为算术平均算法;b=1时，式(10)为倒数距离比法;b=2时，式(10)为距离平方的倒数比法，此文使用了广泛使用的距离平方倒数比法。

待控制站分区内逐栅格雨量插值后，计算控制站分区内所有栅格雨量时间序列的算数平均值，可得到控制站分区的面雨量时间序列。距离反比权重插值见图6。

图6 距离反比权重插值

4.移植结果分析

尖岗水库7月18日—22日5天累积雨量950.8mm，其中7月20日15～16时最大1h雨量147mm。暴雨中心移植至不同位置时主要区间面雨量统计见表1。移植降雨量如下：方案一移植至冶河微水位置，子牙河流域面雨量184.3mm，其中滹沱河岗南水库以上241.6mm，岗黄区间448mm，滏阳河山区产流区214.5mm，暴雨中心位于冶河中游区域。方案二移植至滏阳河临城水库，子牙河流域面雨量202.9mm，其中滹沱河岗南水库以上110.8mm，岗黄区间382.5mm，滏阳河山区产流区339.9mm，暴雨中心位于滏阳河山区水库上游。通过对比，两种方案仅滹沱河降雨较“63·8”小，方案一冶河降水略大于“96·8”和“16·7”降水，滏阳河山区降水小于“96·8”和“16·7”降水。方案二冶河降水小于“96·8”和“16·7”降水，滏阳河山区降水接近“96·8”降水，大于“16·7”降水。总体上，移植降水与“96·8”降水量级大致相当。□

表1 暴雨中心移植至不同位置时主要区间面雨量统计