基于Morris—Sobol的临界雨量参数敏感性分析

原文林 高倩雨 张晓蕾 郝鹏

摘要：山洪灾害临界雨量计算参数众多，部分取值受主观因素影响较大。对此，应用敏感性分析理论，建立基于Morris-Sobol的临界雨量计算参数敏感性分析二层模型。该模型首先利用修正的Morris筛选法对参数进行定性分析，以定性分析结果为基础，选取中等灵敏度以上参数作为定量分析的参数，再通过Sobol指数法对选取参数进行定量敏感性分析。通过定性分析与定量计算相结合的分层递进，分析计算各参数对临界雨量的影响程度。以河南省安阳县都里乡都里前街小流域为研究对象，对其临界雨量计算过程中的参数敏感性进行分析，结果表明：最大th点雨量均值H1的敏感性最大，说明在临界雨量计算中，短历时暴雨对计算结果的影响较大。不论是利用设计暴雨还是实测降雨进行临界雨量计算，均应注意短歷时暴雨的资料精度。

关键词：山洪灾害；临界雨量；参数敏感性分析；Moms-Sobol二层模型：水位流量反推法

中图分类号：TV121+.1

文献标志码：A

doi： 10. 3969/j.issn.1000-1379.2018.07.008

1引言

山洪灾害是具有成灾迅速、破坏力强、危害严重等特点的一种自然灾害，国内外一般选取临界雨量作为山洪灾害预警预报的重要指标。临界雨量是指导致一个流域或区域发生山溪洪水可能致灾时，降雨量达到或超过的量级和强度

。

目前，国外最具代表性的临界雨量计算方法为美国的Flash Flood Guidance（ FFG）方法。国内山洪灾害临界雨量计算方法主要有数据驱动法和水文水力学法。我国山丘区面积占国土面积的近三分之二，保证山丘区小流域预警预报质量才能有效减少灾害损失，而山丘区小流域雨量、水文等监测设施不完善，缺乏实测资料，临界雨量计算通常采用水文水力学法。在应用水文水力学法时，部分计算参数取值受主观因素影响较大，需通过参数敏感性及其相互关联性分析，辨识参数的敏感程度，确定各参数对临界雨量计算结果的影响，为临界雨量计算及参数率定提供参考依据。

敏感性分析一般分为局部分析和全局分析，局部敏感性分析仅针对单个参数对计算结果的影响，全局敏感性分析考虑整个参数空间对计算结果的影响8。全局敏感性分析方法又可分为定性和定量两类方法，定性分析法包括Morris筛选法、多元回归法和LH -OAT方法等，定量分析法包括Sobol指数法、ExtendFAST法和GLUE方法等。修正的Morris筛选法利用微分计算各参数的敏感性，具有计算量小且能够处理多个输人参数的优点，但缺乏参数间交互效应的分析：Sobol指数法基于方差分解的原理，具有对单个参数的主效应、全效应及多个参数的交互效应进行分析的优点，在很多领域得到应用，但该方法计算量较大。

山洪灾害临界雨量计算涉及参数众多，为了研究各参数对临界雨量的影响程度以及各参数间的交互效应对临界雨量的影响，本文利用修正的Moms筛选法计算量小的优点对各参数进行敏感性定性分析，结合Sobol指数法的定量计算，建立基于Morris-Sobol的临界雨量计算参数敏感性分析二层模型。以河南省安阳县都里乡都里前街小流域的临界雨量计算为例，结合临界雨量计算模型中的参数，辨识计算参数对临界雨量的影响程度及相互关联性，为临界雨量计算参数取值提供参考依据。

2临界雨量计算参数分析

水位流量反推法是水文水力学法中最典型的方法，本文主要针对水位流量反推法计算临界雨量中涉及的参数进行敏感性分析与讨论。水位流量反推法假定降雨与洪水同频率，根据河道控制断面成灾水位，由水位流量关系计算对应的成灾流量，由流量频率曲线确定成灾洪水频率，最后由降雨频率曲线确定临界雨量。水位流量反推法计算临界雨量主要包括设计暴雨计算、设计洪水计算和水位流量关系推求三部分。

（1）设计暴雨计算参数。山丘区小流域设计暴雨主要是基于当地暴雨图集等资料进行计算，即暴雨图集法。根据研究对象所处地理位置，利用暴雨图集确定相应参数取值并进行设计暴雨计算。涉及的参数主要有最大24h点雨量均值H24、最大1h点雨量均值H1、变差系数Cv和暴雨递减指数n1、n2、n3。

（2）设计洪水计算参数。小流域设计洪水的计算方法主要有推理公式法、等流时线法和瞬时单位线法等，本文采用推理公式法进行设计洪水计算，其中涉及的参数为降雨损失参数μ。

（3）水位流量关系参数。水位流量关系推求常用方法为曼宁公式法。曼宁公式中涉及的参数主要有糙率n和水面比降i。

3基于Morris-Sobol的临界雨量计算参数敏

感性分析二层模型

在进行参数敏感性分析时，为了提高计算效率，笔者提出了Morris筛选法定性分析和Sobol指数法定量计算的临界雨量计算参数敏感性分析二层模型，其基本流程见图1。

3.1定性分析

山洪灾害临界雨量计算中的参数敏感性定性分析采用修正的Morris筛选法，其基本原理为选取山洪灾害临界雨量计算模型中的某一参数xi，以固定步长C进行扰动变化，最大变幅为M，计算其灵敏度判别因子S，其计算公式为式中：Yo为参数xi初始取值对应的临界雨量；Yi、Yi+1分别为参数xi第i次、第i+1次扰动变化后对应的临界雨量；Pi、Pi+1分别为参数xi第i次、第i+1次扰动变化后的取值相对于初始值的变化率，%：Z为扰动变化总次数，即试验次数，其值由固定步长C和最大变幅M决定。

根据各参数的敏感性判别因子进行参数敏感性分级，敏感性分级标准见表1。

根据各参数敏感性定性分析结果，即敏感性分级情况，去除不灵敏参数，并对剩余参数进行敏感性定量分析。

3.2定量分析

山洪灾害临界雨量计算中的参数敏感性定量分析采用Sobol指数法，其基本原理是将临界雨量计算系统分解成单个参数及参数间相互组合的函数，计算单个参数和参数组合的方差，并计算分析各部分方差对总方差的影响，即分解方差，以此来分析参数的重要程度和参数间的相互影响程度

。

假设Y=f（x）=f（x1，x2，…，xi，…，xn），其中：Y 为临界雨量，xi为敏感性定量分析选取的参数。根据方差分解理论有

主效应指标Sxi反映单个参数xi单独对临界雨量的影响，主效应指标越大，说明该参数对临界雨量的影响越大。全效应指标ST包含了参数xi的主效应和该参数与其他参数的交互作用，若参数的全效应指标较大，说明该参数不仅对临界雨量影响较大，而且与其他参数的交互效应也很大。

采用蒙特卡洛法模拟得到参数的一阶、二阶及更高阶次的敏感度。在利用修正的Morris筛选法进行定性分析的基础上，根据灵敏度判别因子S选取山洪灾害临界雨量计算中的k个参数进行定量分析，在k个参数的取值范围内利用拉丁超立方抽样t次，抽取A和B两个矩阵：

矩阵的每一行表示k个参数的组合，每一列表示t次随机抽取参数xi的值，将式（5）中矩阵A的第i列换成矩阵B的第i列，其余列保持不变记为Ci；再将矩阵B的第i列换成矩阵A的第i列，其余列保持不变，得到矩阵记为C-i。

利用A、B、Ci、C-i矩阵中的每一列参数分别进行临界雨量计算，并按照式（7）至式（12）进行相应的方差估计和参数敏感性指标计算。式中：f为矩阵A和B均值的估计值；V（Y）为临界雨量的方差估计值；Ui为矩阵A和Ci均值的估计值；U-i为A和C-i均值的估计值；Sxi为主效应指标的估计值；ST为全效应指标的估计值。

4实例分析

4.1研究区概况

以河南省安阳县都里乡都里前街小流域为研究对象，分析山洪灾害临界雨量计算中的参数敏感性。安阳县地处河南省北部，北纬35°35-36°21'，东经113°35'-114°45，地势由西到东呈阶梯降低，有山区、丘陵、平原、洼地多种地表形态。安阳县属暖温带季风气候区，全年气温变化较大，四季分明，因受季风影响，故降雨时空分布极不均匀，多年平均降雨量510mm。研究小流域形状见图2，小流域面积为69.11km^2，河长为21.07km。研究对象控制断面成灾水位为220.15m。

根据《河南省中小流域设计暴雨洪水图集》，结合研究区自然地理、地形地貌、水文气象等情况及其所在水文分区，确定各分析参数取值范围及初始值，见表2。

4.2定性分析结果

根据修正的Morris筛选法，以步长C=0.5%对某一参数值进行扰动，扰动最大变幅M=±20%，且保证其余参数固定不变。计算各参数的灵敏度判别因子，结果见表3，敏感参数分布见图3。

由表3可知，山洪灾害临界雨量计算中不灵敏参数有3个，分别是变差系数Cv，、暴雨递减指数n1、暴雨递减指数n2；中等灵敏参数有1个，为降雨损失参数μ；灵敏参数有5个，分别是最大24 h点雨量均值H24、最大1h点雨量均值H1、暴雨递减指数n2、糙率n和水面比降i。

在模型第一层定性分析结果的基础上，去除不灵敏参数，减少定量分析輸入参数，提高定量分析计算效率，同时对研究系统中的参数进行充分的分析计算，初步选取中等灵敏参数及以上等级共6个参数作为定量分析的参数，即降雨损失参数μ、最大24h点雨量均值H24、最大1 h点雨量均值H，、暴雨递减指数n2、糙率n和水面比降i。

4.3定量分析结果

根据选定的6个敏感性定量分析参数及其取值范围，进行拉丁超立方抽样，得到2000个样本，采用Sobol指数法计算6个参数的主效应指标和全效应指标估计值，结果见表4。

由表4可见，6个参数的敏感性由强到弱依次为最大1h点雨量均值H、糙率n、水面比降i、暴雨递减指数n2降雨损失参数μ、最大24h点雨量均值H24。

最大1h点雨量均值H，不仅自身的变化对临界雨量计算结果的影响远大于其余参数，且其与其他参数的交互效应对临界雨量的影响程度也最大，这说明在临界雨量计算中，短历时暴雨对计算结果的影响较大。因此，不论是利用设计暴雨还是实测降雨进行临界雨量计算，均应注意短历时暴雨的数据精度。

糙率n对临界雨量的影响程度也较大。在实际工作中，其取值通常根据河道特征参照天然河道糙率表进行确定。为了提高临界雨量计算精度，在确定山洪灾害临界雨量计算参数值的时候，应重点分析糙率n的取值。当实测水文资料充足时，应利用实测水文资料进行推算：当实测水文资料缺乏时，应综合考虑河流的沟道特征、植被生长和床面粗糙情况等多方面因素，选取一个最佳值。

在确定山洪灾害临界雨量计算参数取值时，按照参数的敏感性大小排序情况，加强对敏感性较强参数的研究，适当减少对敏感性较弱参数的研究，提高临界雨量计算精度，保证山洪灾害预謦预报精准度。

此外，通过基于Morris-Sobol的临界雨量计算参数敏感性分析二层模型对计算中涉及的设计暴雨参数、设计洪水参数和水位流量关系参数进行敏感性分析，综合分析结果可知，水位流量关系参数糙率n和水面比降i的敏感性分别排在第2位和第3位，对临界雨量的影响均较大。在确定山洪灾害临界雨量计算参数取值时，应加强水位流量关系参数的研究，提高山洪灾害预警预报质量。

5结论

利用Morris筛选法和Sobol指数法建立敏感性分析模型，通过定性分析与定量计算相结合的分层递进，分析计算参数对临界雨量的影响程度，可为山洪灾害临界雨量计算时参数取值提供参考依据。以河南省安阳县都里前街小流域为例的分析，结果表明：在山洪灾害临界雨量计算参数中，最大1h点雨量均值H，对山洪灾害临界雨量的影响程度最大，因此在山洪灾害临界雨量计算时应注意短历时暴雨的资料精度。

此外，通过对各参数的敏感性进行综合分析，水位流量关系参数的敏感性较强，因此应对水位流量关系参数进行重点研究，尤其是糙率n对山洪灾害临界雨量的影响程度较大，应利用实测水文资料进行推算，或综合考虑河流的沟道特征、植被生长和床面粗糙情况等多方面因素，选取最佳值。

本研究仅针对最为常用的水位流量反推法进行了参数敏感性讨论，在今后的分析研究中可针对多种方法分别研究不同方法下的敏感性参数，以便实际应用选取不同方法时作为参考。