山洪预报模型历史数据影响及解决方法探析

2022-12-22苏晓林秦根泉黄少辉林智伟丁呈涛

防灾减灾学报 2022年4期

苏晓林，秦根泉，黄少辉，林智伟，丁呈涛，3

（1.福州大学土木工程学院，福建福州 350000；2.四创科技有限公司，福建福州 350000；3.福建省防灾减灾信息应用工程技术研究中心，福建福州 350000）

0 引言

山区小流域洪水常在暴雨后发生，具有汇流时间短、水位陡涨陡落、成灾时间快的特点，因而具有很强的突发性和破坏性，是我国防汛工作的难点和薄弱环节。随着山洪灾害防御体系在全国各地项目建设中陆续建立，我国山洪防治取得了许多成果。据统计，1961年以来，中国极端强降水事件呈增多趋势，而2011年至2019年山洪灾害年平均死亡人数近333人，较项目实施前的2000年到2010年的1179人大幅减少七成[1]。但是由于我国山区地形复杂、降雨时空分布不均，加上城镇化的发展，导致山洪风险依旧严重。据统计，绝大部分已建的山洪监测预警系统还处于初级阶段，随着我国山洪灾害防治非工程措施建设的深入，山洪监测预警系统建设正面临从“感知型”向“智能型”的转变，即从山洪灾害调查评价得到的静态风险区成果进行经验性预报，向以山洪模型模拟成果进行预警转变，以期得到更加精准和实时的山洪灾害预测预警来指导避险[2]。

鉴于历史水文资料的缺失错漏是制约山洪预报模型成功搭建的重要因素，本文以海南白沙河流域为例，探索了一种通过选取典型流域进行洪水预报模型参数率定，再采用水文模型参数物理特性相似法将得到的参数方案移植到其它子流域的方法，以期解决山洪小流域洪水预报模型当前普遍存在的数据缺失问题。

1 研究现状

1.1 山洪预报模型研究成果

当前我国的山洪预警预报方法总体上可分为两种：一种是以统计归纳所得的临界预警指标如临界雨量阈值进行预报的经验性方法；另一种是以模拟降雨产流过程的洪水预报模型进行预报的机理性方法[3]。随着研究的发展，也出现了结合两种方法的动态临界预警指标法等。一般洪水预报模型从实现原理上可分为水文模型、水动力模型和人工智能预测模型。其中由于新兴的人工智能预测模型的预测精度和输入的历史数据量有关，所以难以在普遍资料不足的山洪流域满足精度要求。同样的，严格基于物理过程的水动力模型具有更高的准确性，但因其对河道地形等资料的高要求，也不适合应用在无资料的山洪流域，而半经验半物理过程的水文模型更加简单、高效，在山区暴雨洪水模拟上应用更为广泛[4]。

当前常用的洪水模拟水文模型方法主要有：新安江模型、HEC-HMS（Hydrologic Engineering Center-Hydrologic Modeling System）模型、TOPMODEL（Topgraphy Based Hydro-logical Model）模型、FFMS（Flash Flood Modul Simulation System）模型、NAM模型等[5-7]。新安江模型是赵人俊团队研发的集总式水文模型，在较湿润地区有较好的适用性，研究应用成果较多；HEC-HMS模型是美国陆军工程兵团开发的半分布式次降雨径流模型，在长时间尺度的模拟应用上较少，主要应用于湿润和半湿润地区；TOPMODEL模型是基于地形的半分布式流域水文模型，它考虑了地形因素，但忽略了土壤、植被等要素，主要适用于小流域尺度和半干旱地区；FFMS模型是刘昌军团队研发的模块化小流域分布式水文模型，它是基于超渗、蓄满时空动态组合的原理，主要运用于小流域暴雨洪水的模拟；DHI（Danish Hydraulic Institute）开发的集总式概念性的NAM模型，具有较简单的结构和明确物理意义的参数，对国内外许多流域的降雨径流过程都有较好的模拟效果，此外由于NAM模型对数据量的要求相对较少，因此在资料条件较差的区域也具有一定的可行性[8-11]。

1.2 模型数据存在的问题

依据相关规范在洪水预报中，必须使用至少10年的水文气象资料，且还应包括各种代表性年份和一定场次洪水资料，以此作为精度评定，确定预报等级，其中对于骤发性洪水（如短历时、高强度暴雨形成的小范围山洪）的预报可采用气象与水文相结合的方法，可不进行精度评定[12]。山洪预报虽然可以不进行精度评定，但在建立水文模型过程中，仍然需要较好的水文资料进行参数率定，才能得到反映地区情况较准确的水文模型参数方案。并且随着监测手段和山洪治理建设的推进，对防洪安全和模型预报精度的要求也将有所提高。

虽然我国对山洪的治理工作已经取得初步成效，但由于山洪高发地区通常在较为偏远的山区，限于较差的实际条件，通常难以得到足够多的符合质量要求的历史序列数据资料，使得山洪模型在构建和参数率定等方面存在许多困难。当前监测数据主要存在的问题有：①水文监测种类较少，河道通常只有水位数据缺乏相应断面的水位流量关系；②水文资料的代表性不强，不能很好地代表流域情况；③水文资料存在不连续、不稳定的问题，数据质量不佳，降低了数据可靠性，也影响数据的筛选使用；④水文资料时间序列过少，数据量不足，难以满足预报要求。

分析现有的数据问题，可以发现影响数据质量的原因主要有：①山区洪水流速快且含沙量大难以准确监测，监测设备不够齐全，监测手段较为落后；②监测站点位置不够合理，例如雨量监测站的密度不够，不能准确代表整个流域的情况；③监测站点级别较低，多为自建站点，监测仪器的稳定性不足，存在缺测、误测等问题；④山洪区域监测站点建设时间晚，监测时间短，缺乏多年长序列实测资料。

1.3 数据缺陷解决思路

在缺乏资料或资料质量欠佳的山洪高发地区，水文模型参数受限于现状条件而难以率定，对此在实际工程中通常是采用水文模型参数区域化的方法，将其它有资料地区的水文模型参数参照运用于本地区中。当前水文模型参数区域化方法主要可以分为两类：一类是参数移植法，主要有空间邻近法和物理特性相似法等；另一类是参数回归法，主要有多元回归法和全局平均法等[13]。空间邻近法是借用邻近地区的模型参数的方法，它假定相近地区具有相似的水文条件，以空间距离作为主要判定依据，对于地形地貌变化不大的地区有较好的准确性。物理特性相似法是将地理条件、气象条件等流域特征作为参数移用的判断依据，与特征参数的选择有较密切的关系。多元回归法是建立模型参数与选取的地形特征参数的回归方程，对特征参数选取和回归方式有较高的要求，且需要较多的流域样本。全局平均法是将流域内划分的有资料的小流域所得参数进行算术平均，用以代表整个流域特征情况，一般适用于不大的流域。

由于参数回归法对数据条件要求较高，因此当前应用更为广泛的是参数移植法，而山洪预报机理性方法是基于水循环过程对降雨径流过程的模拟，各基本假定和参数具有一定的物理意义，采用物理特性相似法更为合理。当前对流域相似性的判断主要考虑气候条件和下垫面特征[14]。对于处在同一流域中特别是同一小流域中的子流域来说，其气候条件和下垫面特征通常都不会有太大的变化，具有较大的相似程度，在参数的移植借鉴上有着较高的可信度。NAM模型的参数在物理意义上主要与下垫面特征有关，对于具有相似下垫面特征的山洪区域，模型参数具有较小的差异性，在参数移植上有较好的可行性。

2 参数确定方法

2.1 研究区域概况

本文选取海南省的乐东县白沙河流域开展研究。白沙河位于海南省的西南面，河道比降9.95‰。干流全长26km，流域集雨面积170km2。流域处亚热带海洋性季风气候区，雨水十分充沛。白沙河流域中的水文站点数量较少，水位站点仅有两个，一是流域出口处的丹村水位站（已损坏且缺乏历史水位资料），另一个是支流上的牛尼头水位站；而雨量站的数量也较少，且多分布于白沙河上支流区域。因此，对于整个白沙河流域洪水预报模型的搭建，站点分布及数据均不能满足要求。子流域划分采用ArcGIS水文分析工具处理白沙河流域DEM数据，以重要断面（如水位站点）和支流汇入断面作为流域划分依据，将白沙河流域划分为10个子流域。选取白沙河流域内历史数据情况相对较好的牛尼头水位站点控制的子流域5作为典型区域，子流域5内有尖峰岭和牛尼头两个雨量站，根据子流域形状和站点分布，由泰森多边形作图法确定其权重分别为0.9和0.1（图1）。

图1 白沙河流域划分及站点分布图Fig.1 Map of Baisha River Basin and hydrological stations

通过NAM模型对子流域5进行山洪模拟和参数率定，再通过对子流域典型特征的比对将参数移植到其它子流域模型。牛尼头水位站建站年月为2013年1月，其上游控制子流域面积约为11km2。牛尼头水位站的实测水位数据序列为2015年6月1日—2020年12月8日，有5.5年，本次基于水位数据变化和降雨数据变化相匹配的原则，对数据准确性和完整性进行判断，筛选可用洪水场次，结果详见表1。

表1 牛尼头水位站水位监测数据情况表

从表1可知，该站点的实测水位序列不足，质量也欠佳，中间存在许多缺测和误测数据，而且雨量站存在部分降雨场次的缺失，综合考虑历史资料的完整性和代表性后选取3场较为可信的降雨和水位数据供模型参数率定和验证。

2.2 参数率定方法

2.2.1 模型原理

NAM模型将降雨径流过程分为四个相互影响的储水层的计算，其具体结构见图2[15]。

图2 NAM模型计算结构图Fig.2 NAM model calculation structure diagram

模型的参数是对流域平均情况的代表取值，一般无法通过实测获得，只能通过率定得到适合该流域的参数方案，参数物理意义及取值范围见表2[16]。

表2 NAM模型主要参数意义及取值范围

2.2.2 模型率定过程

水文产汇流模型参数率定需要将实测水位数据转换成对应流量数据。根据实测数据，牛尼头监测断面属于宽浅型河道，且断面附近河道较为顺直，断面尺寸变化不大，可以采用如式（1）的曼宁公式计算水位流量关系。

式中，Q为流量，单位m3·s-1；A为断面过流面积，单位m2；v为过流断面平均流速，单位m·s-1；R为水力半径；n为糙率，取值0.040；J为洪水水面线比降，由于缺乏资料，实际采用河道比降，取值8.8‰。

模型参数的确定，采用自动率定和人工率定相结合的方式。

自动率定采用SCE-UA随机搜索参数优化算法，该算法自1992年由Duan提出以来在水文模型参数率定中得到广泛应用和验证[17]。自动率定的目标为：①总水量平衡，即模拟的径流量与实测值相等；②洪水过程线的形状相近，变化规律类似；③洪水流量的峰值接近，主要考虑洪峰出现时间、洪峰流量大小以及水量等条件；④低流量吻合，即基流情况吻合[15]。

在进行自动率定之后，再进行人工率定，反复进行参数调整步骤，最终得到最符合各场次洪水的总体情况的参数方案。人工调整步骤如下：①通过调整参数Lmax和Umax来平衡水量；②通过调整参数CK1、CK2和CQOF来调整洪峰流量的大小；③通过调整参数CKBF来调整基流[15]。使用Nash-Sutcliffe系数评价模拟结果与实测资料的符合程度，计算公式见公式2。

式中，R2为Nash-Sutcliffe系数，越接近1表示模拟结果和实测资料的差异性越小；Qobs（i）为实测流量，单位m3·s-1；Qsim（i）为模拟流量，单位m3·s-1；为实测流量平均值，单位m3·s-1；N为实测流量资料序列的长度；wi为权重系数。

3 参数结果与应用

3.1 参数率定结果分析

根据站点现有数据情况，尽量多选取可用于模型计算的山洪场次进行模型率定和验证。经分析，选取2016年的2场山洪进行模型参数率定，2017年的1场山洪作为验证，参数率定的模型模拟结果见表3和图3-4。基于同一模型参数方案，不同场次降雨间模拟效果差异的原因还可为山洪预报对降雨分布的敏感性，不同场次降雨雨强中心和站点断面的距离存在差异，这就导致雨量站的实际分配效果和由泰森多边形计算出的雨量分配权重方案有所差异，且将直接影响到汇流时间的长短，从而引起模型参数之间的差异。同时降雨强度的不同会影响径流量和下渗量的分配比例，也是影响模型参数的原因之一。

从表3中可以看出，洪水率定期选取的2个场次洪水的模拟洪峰相对误差均在洪峰许可误差20%之内，符合洪峰模拟精度要求。从图3、4的模拟结果中可以看出，模拟洪水过程和实测情况整体吻合性较好，基本上能够正确地反映出洪水过程，出现洪峰的时间也较为接近。一般情况下洪峰出现时间要晚于降雨峰值出现的时间，这一雨洪滞后性也在模拟结果中得到了验证。

表3 牛尼头站断面NAM模型场次山洪模拟结果

图3 20160924场次洪水实测和模拟结果对比图Fig.3 Comparison of measured and simulated results of 20160924 floods

图4 20161018场次洪水实测和模拟结果对比图Fig.4 Comparison of measured and simulated results of 20161018 floods

但是，其中20160924场次洪水模拟的Nash-Sutcliffe系数不佳，洪水模拟过程也不甚理想。分析原因是在率定参数的过程中，为了使2场洪水总体模拟效果较好，在自动率定之后进行手动参数调整。考虑到20161018场次洪峰流量更大，风险程度更大，所以在参数调整过程中，更加侧重该场次的模拟效果，在保证洪峰误差和峰现时间误差满足要求的前提下，适度牺牲20160924场次模拟精度。

受限于较少的可用洪水场次，且考虑到偶然性的影响，难以给出多种工况情景下参数方案，同时也为了便于模型的部署应用，综合考虑现有数据资料率定结果，最终确定一套参数方案，如表4所示，各参数均在合理范围之内，将被作为该地区的参考方案。

表4 山洪模拟NAM模型主要参数取值方案

3.2 模型模拟验证分析

选取2017年的一场山洪进行参数验证，模拟结果见表5和图5。从表5可以看出，验证期的洪水场次模拟洪峰误差和峰现时间误差基本满足要求，Nash-Sutcliffe系数满足要求，此外从图5的模拟结果也可以看出洪水过程的模拟情况同实测较为吻合，说明采用的参数方案具备一定的准确性。

图5 20170826场次洪水实测和模拟结果对比图Fig.5 Comparison of measured and simulated results of 20170826 floods

表5 牛尼头站断面NAM模型场次山洪模拟验证结果

3.3 参数移植方法

3.3.1 参数移植应用思路

由于不同山洪流域的土壤和地形可能存在较大差异，因而不能直接使用代表性区域的参数方案，本次研究根据物理特性相似法制定NAM模型在山洪流域的参数移植方法，参数移植思路分两种：一是分析模型参数的敏感性，筛选对模拟结果影响较大的模型参数和相关的流域物理特征，并建立模型参数和对应流域物理特征之间的数学关系，以流域物理特征计算模型参数取值；二是选取足够数量的具有不同流域物理特征的研究区分别进行参数率定，进而组合成一套适用多种类型流域的参数方案，建模时通过流域物理特征选取参数方案。本研究区域流域物理特征如表6所示。

表6 白沙河牛尼头小流域物理特征表

要建立整个白沙河流域的降雨径流模型，需要先建立各个子流域的水文产汇流模型，然后建立各个子流域汇流出口到流域总出口断面的河道水动力模型。由于在白沙河流域内，除牛尼头子流域外的其它子流域都缺乏站点和资料，采用上述方法二，以牛尼头子流域的参数方案为基础，对模型参数进行调整，移植运用到其它子流域，从而搭建出整个白沙河流域的降雨径流模型。

3.3.2 参数敏感性与相关流域物理特征分析

模型参数的敏感性指的是参数变化对模拟结果影响的大小，在现有的NAM模型的参数敏感性分析研究中，姚璇[18]采用基于扰动分析法的相对灵敏度方法对NAM模型9个主要参数进行了局部敏感度分析，结果表明参数Lmax、Umax、CQOF和TG（参数含义参见表2）具有较高的敏感性，而CK1和CK2在日径流模拟结果影响较小。陈福容等[19]则认为参数Lmax、CK1和CK2是较为敏感的参数。此外李磊等[20]研究发现参数Lmax、Umax、CQOF对径流总量影响明显，CK1、CK2对峰值的影响明显。研究普遍认为参数Lmax、Umax、CQOF是对水量平衡有较大影响的参数，而影响汇流时间的CK1、CK2同样具有较大的敏感性，且在本次手动率定过程中同样发现了类似的规律，因此对模拟结果影响较大的这5个参数，需要在参数移植中重点考虑。根据各参数的物理意义，总结其主要相关的流域特征参数，如表7所示。

表7 NAM模型部分参数及主要相关流域特征参数表

参数Lmax、Umax分别是浅层蓄水层最大含水量和地表蓄水层最大含水量，通常前者约为后者的10倍，主要与土壤有关，且强烈受到植被的影响。与洪峰流量有显著相关的CK1、CK2和CQOF这三个参数和流域特征之间一般具有如下关系：当土壤透水率越大时，CK1、CK2和CQOF越小；当坡度越大时，CK1、CK2越小，而CQOF越大；当植被覆盖率越大时，CK1、CK2也越大，而CQOF越小；当集水面积越大而流域形状系数越小时，CK1、CK2越大，而CQOF变化较小。

3.3.3 白沙河子流域参数调整取值

白沙河各子流域的流域特征情况基本类似，因此，可以考虑在牛尼头子流域模型参数方案的基础上调整Lmax、Umax、CQOF、CK1和CK2参数，而保持其它参数不变。白沙河各个子流域的模型参数取值方案如表8所示。

表8 白沙河各子流域物理特征和模型参数取值表

采用上述参数方案建立白沙河各子流域的NAM模型，得到各子流域汇流的径流量，再通过水动力模型将其汇流到流域总出口断面（即子流域9汇流断面处），得到流域总出口断面的径流量过程。图6为干流沿程各断面位置处的流量模拟结果。从图6可以看出，从河源到流域出口，沿程河道断面的洪水流量不断增加，且峰现时间也逐渐延后，符合洪水演进的正常规律，后期将结合流域出口新建的水情站点监测数据，对模型参数进行优化，逐步提高模型预报精度。