地貌单位线在卡洛特水电站洪水预报中的运用

Application of Geomorphic Unit Line in Flood Forecasting of Karot Hydropower Station

摘要: 为了解决在卡洛特水电站洪水预报中面临的大流域无或缺水文资料的难题，寻找一种能确保预报精度的有效预报方法。通过采用90 m分辨率的DEM数据，用ArcGIS软件提取地貌参数，构建基于统计物理学的地貌单位线。选用2017~2018年流域内发生的2场全流域洪水进行汇流模拟，结果表明：地貌单位线在集水面积大于2000 km²的流域也可取得较好的汇流模拟效果。地貌单位线可成为解决该流域无或缺资料区域洪水预报问题的有效方法。

Abstract: In order to solve the problem of without or lacking hydrological data in large watershed in flood fore-casting of Karot Hydropower Station, an effective forecasting method to ensure the accuracy of forecasting is sought. By using DEM data with 90 m resolution and ArcGIS software to extract geomorphic parameters, geomorphic unit line based on statistical physics is constructed. Two floods occurred in the whole basin from 2017 to 2018 were selected to carry out the confluence calculation. Geomorphological unit line can also achieve good simulation results in watersheds with catchment area larger than 2000 km². Geomorphic unit line can be an effective way to solve the problem of flood forecasting in areas without or lacking data in the basin.

1. 引言

随着“一带一路”倡议的提出，中国水电站建设事业也快速向周边国家发展。伴随水电建设事业的不断推进，建设区域内水文资料缺乏愈发普遍，国内外目前使用较成熟的水文预报方法均局限于水文资料的缺乏而无法达到理想精度 [1]。目前常用的水文比拟法、经验公式等方法，应用在流域面积较大的国外水电项目上，水文情势跟可参照流域存在较大差异，具有一定的任意性，预报精度无法得到保障 [2]；率定合理参数需要长时间水雨情资料，这跟电站建设周期和对洪水预报需求矛盾。因此探索在集水面积较大且无资料区域的水文预报方法就显得尤为重要和迫切。

Rodriguez-Itrube和Gupta [3] 等人设想降落在流域上的雨水是由许多弱相互作用的雨滴所组成，应用统计物理学的方法，提出了地貌瞬时单位线理论。在此基础上，国内水文研究者进行了进一步研究和应用，特别是随着信息技术的发展，GIS和数字高程模型(DEM)在水文领域广泛应用。1988年文康、李琪等 [4] 推导出地貌单位线通用公式，并应用在集水面积64~333 km²的12个流域内且拟合成果较好；1988年谢平 [5] 用集水面积为797 km²的密赛流域对地貌单位线预报精度进行验证；2008年林靓靓用DEM数据对集水面积17.8 km²冯家圪垛流域进行汇流模拟取得不错效果 [6]；2012年孙龙、石鹏等将基于DEM的地貌单位线应用在集水面积96 km²榆村流域 [7]；2018年纪小敏、陈颖冰等对江苏省无资料山丘区洛阳河流域(集水面积150.26 km²)进行径流模拟，模拟精度较高 [8]。上述研究认为地貌单位线汇流模拟效果甚好，但研究流域都为集水面积小于1000 km²小流域，地貌单位线在较大集水面积的参数率定和使用研究文章较少。结合巴基斯坦卡洛特水电站洪水预报实际情况，开展了在较大集水面积条件下地貌单位线汇流模拟精度的研究，对类似项目解决流域内无或缺资料区域洪水预报提供了一个有效途径。

2. 地貌瞬时单位线原理

地貌单位线主要由发生场，即初始状态概率 $Q_i\left(0\right)$ 和传扩场，即时段t内的状态转移概率 ${\theta }_{ij}\left(t\right)$ 组成，亦即单位水体流经出口断面的时间分布函数 [9]。基于地貌单位线理论的汇流模型需要输入霍顿(Horton)河数率、河长率和面积率等流域信息。

地貌瞬时单位线(R-V GIUH)利用斯特拉勒分级方案 [10] 和HORDUN地貌参数推求流域汇流过程，地貌单位线通用公式如下 [4]。

$u\left(t\right)=-{\displaystyle \underset{j=1}{\overset{\Omega }{\sum }}{\lambda }_j\left[{\displaystyle \underset{i=1}{\overset{j}{\sum }}{\theta }_{i,\Omega }\left(0\right)A_{ij}}\right]}\exp \left(-{\lambda }_{j}t\right)$ (1)

$A_{ij}$ 为关于 ${\lambda }_j$ 和 $P_{ij}$ 函数的概化系数：

$A_{ij}=\left[B_{ij}\left(\Omega \right)/{\displaystyle \underset{a=1}{\overset{\Omega }{\prod }}\left[-{\lambda }_j\left({\lambda }_a-{\lambda }_j\right)\right]}\right]\left(a\ne j\right)$ (2)

${\theta }_{1,\Omega }\left(0\right)$ 为 $\Omega$ 级流域j级河流水质点的初始概率；

${\theta }_{1,\Omega }\left(0\right)=\frac{R_B^{\Omega -j}}{R_A^{\Omega -1}}\left[R_A^{j-1}-{\displaystyle \underset{i=1}{\overset{j}{\sum }}{R}_A^{i-1}{R}_B^{j-1}{P}_{ij}}\right]$ (3)

状态转移概率 $P_{ij}$，表示雨滴由i级河流流入j级河流概率。

$p_{ij}=\left[\left(R_B-2\right)R_B^{\Omega -j-1}\frac{{\displaystyle \underset{k=1}{\overset{j-i-1}{\prod }}\left(R_B^{\Omega -j-k}-1\right)}}{{\displaystyle \underset{k=1}{\overset{j-i}{\prod }}\left(2R_B^{\Omega -j-k}-1\right)}}\right]+\frac{2}{R_B}{\delta }_{i+1,j}$ (4)

${\delta }_{i+1},j=\{\begin{array}{l}1\left(j=i+1\right)\\ 0\left(j\ne i+1\right)\end{array}\left(i,j=1,2,\cdots ,\Omega \right)$ (5)

式中 ${\lambda }_a$ 、 ${\lambda }_j$ 分别为a、j级河流的平均汇流时间，其由流域平均流速确定； $R_A$ 为面积率， $R_B$ 为河数率， $B_{ij}\left(\Omega \right)$ 为 $S={\lambda }_j$ 代入 $P{\left(S\right)}_{j,\Omega +1}$ 后的计算结果。本文对研究流域提取三级河网。

根据公式(1)至(5)计算三级河网地貌瞬时单位线见式(6)。将式(6)对时间t进行积分，得到S(t)曲线，最后根据S(t)曲线原理，转换成时段单位线。

$\begin{array}{c}\mu \left(t\right)={\theta }_1\left(0\right)[\frac{{\lambda }_1{\lambda }_2\left({\lambda }_2-{\lambda }_1{p}_{13}\right)}{\left(\lambda {}_2-{\lambda }_1\right)\left({\lambda }_3-{\lambda }_1\right)}\exp \left(-{\lambda }_{1}t\right)+\frac{{\lambda }_1{\lambda }_2{\lambda }_3{p}_{12}}{\left({\lambda }_1-{\lambda }_2\right)\left({\lambda }_3-{\lambda }_2\right)}\exp \left(-{\lambda }_{2}t\right)+\frac{{\lambda }_3\left({\lambda }_1{\lambda }_2-{\lambda }_1{\lambda }_3{p}_{13}\right)}{\left({\lambda }_1-\lambda {}_3\right)\left({\lambda }_2-{\lambda }_3\right)}\exp \left(-{\lambda }_{3}t\right)]\\ +{\theta }_2\left(0\right)[\frac{{\lambda }_2{\lambda }_3}{{\lambda }_3-{\lambda }_2}\exp \left(-{\lambda }_{2}t\right)+\frac{{\lambda }_2{\lambda }_3}{{\lambda }_2-{\lambda }_3}\exp \left(-{\lambda }_{3}t\right)]+{\theta }_3\left(0\right){\lambda }_3\exp \left(-{\lambda }_{3}t\right)\end{array}$ (6)

$S\left(t\right)={\displaystyle {\int }_0^{t}u\left(t\right)\text{d}t}$ (7)

$u\left(\Delta t,t\right)=\frac{A}{\Delta t}\left[S\left(t\right)-S\left(t-\Delta t\right)\right]$ (8)

式中，A为流域面积， $\Delta t$ 为时段长，其它符号意义同上。

3. 应用实例

3.1. 工程和流域简介

巴基斯坦卡洛特水电站建设项目被列为“一带一路”重点项目和“中巴经济走廊”优先实施能源合作项目。卡洛特水电站所处吉拉姆河流域，绝大部分地区处于克什米尔地区，境内安全环境、河道特征相对较差。吉拉姆河流域多为高山峡谷，洪水陡涨陡落，为典型的山区性河流。径流以融雪水和季节性降雨补给为主；流域全年降雨比较集中，干湿季节明显，呈双雨季特征。基于地理环境、洪水预见期等因素，本文研究卡拉斯站和卡洛特(专用)站区间流域，面积约2039.4 km²，流域水系见图1。

3.2. 水文资料情况介绍

卡拉斯站和卡洛特(专用)站区间，历史水文资料由巴基斯坦水电发展署(The pakistan water and power development authority, WAPDA)负责收集管理。水文要素观测主要以人工观测为主，水位基面为不固定的测站基面，年际间水位起算点不完全一致；流量采用流速仪法施测，主要依托断面附近桥梁开展渡河流量测验；雨量资料绝大部分仅有日雨量。WAPDA刊布的成果基本以日资料为主，难以全面掌握河流水文特性，编制水文预报方案可利用的资料有限，这极大的制约了常用洪水预报方法的使用。为了满足电站洪水预报的需要，流域站网于2016年10月新建，建有10个雨量站，卡拉斯和卡洛特(专用)站2个水文站，站网布置见图1。

Figure 1. Distribution maps of basin & stations

图1. 流域水系与测站分布图

3.3. 地貌参数计算

通过GIS水文分析工具对90 *90 m DEM经过填洼、水流方向提取、汇流累积量计算、河网分级和分水岭计算等处理 [11] [12] [13]。提取的地貌参数能否有效地反应流域的实际情况关键在集水面积阈值的确定，不同大小的集水面积阈值提取出来的模拟河网信息不同。通过设置一系列集水面积阈值提取出不同密度的河网 [14] [15] [16]，分析二者关系可以发现，当集水面积阈值增大时，最初的河网密度减小速率很快，然后在集水面积为40.5 km²时减小速率明显变缓慢。故认为该流域的地貌发育的集水面积阈值选40.5 km²比较合适(见图2)。卡拉斯至卡洛特(专用)站区间提取并计算的地貌参数见表1。

Figure 2. Relation between river network density and catchment area threshold

图2. 河网密度与集水面积阈值的关系

Table 1. Geomorphological parameters of watershed from kallas to Karot (special) station

表1. 卡拉斯至卡洛特(专用)站区间流域地貌参数

3.4. 地貌单位线计算

${\lambda }_i$ 为i级河流的平均汇流时间，是单位线动力因子的反映。本文研究流域坡度陡峻(平均比降约在3‰)且比较均一，可认为流速沿河长基本不变，平均流速对洪水洪峰影响有限 [5]。所以，流域平均流速采用卡洛特(专用)水文站2017年4月6日实测断面平均流速3.38 m/s计算 ${\lambda }_i$。

通过将地貌参数代入公式(6)计算出流域地貌瞬时单位线：

$\begin{array}{l}u\left(t\right)=0.00712\exp \left(-0.361t\right)+0.02764\exp \left(-0.126t\right)\\ \text{+}0.065273\exp \left(-0.044t\right)\end{array}$ (9)

将(9)式通过(7)式积分得到S(t)曲线：

$\begin{array}{l}S\left(t\right)=-0.01972\exp \left(-0.361t\right)-0.21943\exp \left(-0.126t\right)\\ -1.48348\exp \left(-0.044t\right)+1.72263\end{array}$ (10)

根据该流域洪水特性，时段长选用1 h，根据式(8)得到流域时段单位线，见图3。

Figure 3. 1 mm unit line

图3. 1 mm时段单位线

3.5. 地貌单位线洪水计算

根据流域内2017~2018年实测洪水资料，选取仅有2017年4月6日和2018年4月21日两场全流域降雨场次洪水进行模拟分析。采用蓄满产流模型 [2] [17] [18] 计算净雨，地貌单位线进行汇流计算，洪水流量过程见图4和图5，误差统计见表2。由表2可见，地貌瞬时单位线模拟效果较好，洪峰与峰现时间误差分别小于3%和1 h。但在退水阶段，计算洪水小于实测洪水，这可能是由于退水阶段忽视了壤中流等汇流过程。

Figure 4. Comparison of flood discharge processes on April 6, 2017

图4. 2017年4月6日洪水流量过程对比

Figure 5. Comparison of flood discharge processes on April 21, 2018

图5. 2018年4月21日洪水流量过程对比

Table 2. Error statistics of peak discharge and the time of peak

表2. 洪峰流量与峰现时间误差统计

4. 结论

1) 地貌单位线能够密切结合流域地貌气侯特征作流域汇流反应函数的计算，摆脱了从输入输出为已知量反推黑箱内反应函数的手段 [19]，利用地貌单位线进行汇流模拟计算，物理意义更加明确。

2) 对面积为2039.4 km²的卡洛特电站洪水预报流域进行了河网提取，确定40.5 km²为该流域地貌特征的集水面积阈值，计算和率定地貌单位线。通过对该流域的2场洪水进行汇流模拟计算结果令人满意。

3) 地貌单位线在集水面积大于2000 km²的吉拉姆河流域也可较好的运用，表明地貌单位线在大集水面积情况下汇流模拟效果依然较好，为该流域内无或缺资料区域洪水预报提供了有效方法。