基于Sobol方法的MIKE NAM模型参数敏感性分析

王复生，田 娟

(安徽省水利水电勘测设计研究总院有限公司,安徽 合肥 230088)

降雨径流模拟是防洪排涝计算、流域洪水演算以及水资源保护规划的重要研究方法。随着数字化技术的发展，基于模型对流域内的降雨径流过程进行模拟计算已经成为了一种主要趋势。具有物理机制的降雨径流水文模型可以综合考虑降雨、土壤特性、下垫面条件等多方面因素，反演不同边界条件下的降雨产汇流规律及空间分布特征[1- 2]，对洪水演进理论发展和工程应用均具有重要价值。

MIKE NAM模型是一种集总式水文模型，可以用于一个流域或由多个子分区组成的较大河流及由河流、渠道构成的复杂河系的降雨径流模拟过程，在国内外得到了广泛的研究与应用，并取得了良好的模拟效果[3- 5]。NAM模型的参数较多，模拟过程中存在模型参数的不确定性与敏感性，参数和边界条件之间的相互影响等问题[6- 7]。参数取值的合理性以及对参数含义的准确理解对模型精度会产生很大影响，因此对模型参数进行全局敏感性分析，筛选出高敏感性参数进行重点率定十分有必要。

本文以巢湖流域桃溪水文站控制子流域为研究对象，建立MIKE NAM模型，综合考虑NAM模型中各类参数的取值情况，运用Sobol方法对模型参数的敏感性进行定量分析，评估各参数对模型输出结果方差的影响程度，将高敏感性参数作为主要率定参数，以期提高模型率定效率和模拟精度。

1 模型与方法

1.1 MIKE NAM模型

MIKE是由丹麦水资源与水环境研究所(DHI)开发的专业工程软件，适用于河流、灌溉渠道及其它水体的一二维水动力、水质、泥沙运输模拟等[8- 9]。MIKE NAM模型作为一种新兴的水文模拟工具，在流域降雨产汇流模拟中具有良好的适用性，被广泛应用到世界各地不同气象水文条件的流域，在国内长江、黄河、松花江流域均有应用案例，是一个经过大量工程实践验证的模型工具。NAM模型将土壤含水层分成积雪储水层(Snow Storage)、地表储水层(Surface Storage)、浅层储水层(Lower Zone Storage)和地下水储水层(Ground Water Storage)4个部分，模型分别进行连续计算以模拟流域中各种相应的水文过程。

1.2 Sobol参数敏感性分析方法

Sobol方法是一种基于方差分解的定量全局敏感性分析算法，相比于定性敏感性分析方法，该方法能够通过敏感性量化指标的计算直接给出模型参数的敏感性大小，其核心是将目标函数的总方差分解为单个参数的方差和多个参数间相互作用的方差，目前已被广泛应用于经济、环境等多个领域[10]。

Sobol方法假设模型结构可以表示为函数形式u=f(x)，其中模型参数x=x1，x2，…，xn可视为n维离散点，u为标量输出。

假设函数f(x)可积，xi在[0，1]上服从均匀分布，且f(x)满足：

(1)

式中，k=i1，…，is。

则函数f(x)的方差分解表达式可以表示为：

(2)

式中，1≤i1<…is≤n(1≤s≤n)，则被加数共有2n项。

类似的，模型的总方差也可以分解为单个参数和多个参数相互作用的组合：

(3)

式中，D—模型的总方差；Di—参数xi产生的方差；Dij—参数xi和xj相互作用产生的方差，D1，2，…，n—n个参数共同作用产生的方差。

将上式归一化即可得到各参数和参数之间相互作用的敏感性：

(4)

则模型的敏感性指数可以表示为：

式中，Si—参数xi单独作用时的敏感度；Sij—参数xi和xj相互作用的敏感度；STi—参数xi和其它参数共同作用的敏感度；D～i—除参数xi外其它参数共同作用产生的方差。

2 研究区域与模型构建

2.1 研究区域概况

本研究对巢湖流域上游丰乐河桃溪水文站控制流域进行建模分析。

丰乐河发源于大别山余脉，流经双河、桃溪、丰乐至三河镇下大潭湾与杭埠河汇合后流入巢湖，河道全长117.5km，流域面积2124km2，是杭埠河最大的支流。丰乐河干流上的桃溪水文站设立于1951年7月，控制集水面积1510km2，约占丰乐河流域总面积的75%，有长系列的实测流量资料，且精度可靠，是巢湖流域水文分析计算重要的参证站点。丰乐河桃溪站控制流域地形图如图1所示。

图1 桃溪站控制流域地形图

根据研究区域内实测降雨资料及雨量站的分布情况，选取9个代表站点进行降雨径流计算，利用泰森多边形法分析得到各个雨量站的计算权重值见表1。

表1 桃溪站控制流域面雨量计算各雨量站权重

2.2 NAM模型相关参数

NAM模型所需数据资料主要为流域降雨及蒸发数据，模型的初始条件包括开始时刻流域地表储水层和根区储水层的土壤含水率，以及坡面流、壤中流和基流的初始值[11]。NAM模型的输出结果主要包括地表坡面径流、壤中流和基流过程等。模型主要参数及取值范围详见表2。作为集总式水文模型，NAM模型中每个子流域的参数均为本区域的参数平均值。

表2 NAM模型主要参数

3 计算结果与分析

3.1 参数敏感性分析

由于Umax与Lmax相关性较高，可视为一个参数，本文以Umax、CQOF、CKIF、CK1、2、TOF、TIF、TG、CKBF共8个模型参数作为自变量函数，用蒙特卡罗法(Monte Carlo Simulation)对自变量进行抽样，每个参数得到10组样本值，根据Sobol参数敏感性分析方法产生10×(8+2)=100个参数组合，将这100组参数分别带入桃溪站控制流域NAM模型中，对2020年6月8日—8月8日的区域降雨进行模拟。分别以洪峰流量、峰现时间和最大7d洪量作为NAM模型输出结果的目标函数，用一阶敏感性指数Si和全局敏感性指数STi来量化输出结果对各个参数的敏感度高低。模型参数敏感性分析成果如图2所示。

图2 参数敏感性分析成果图

由敏感性分析成果可知：

以洪峰流量作为目标函数时，坡面流汇流系数CQOF和坡面流汇流时间常数CK1、2的一阶敏感性指数、全局敏感性指数远高于其它参数，说明这两个参数对洪峰流量大小起决定性作用。CQOF的一阶敏感性指数大于0，CK1、2的一阶敏感性指数小于0，说明坡面流汇流系数与流量呈正相关，坡面流汇流时间常数与流量呈负相关，这与参数代表的物理意义相符。

以峰现时间作为目标函数时，坡面流汇流时间常数CK1、2的敏感性指数高于其它参数，说明该参数对洪峰出现的时间有较大影响。

以最大7d洪量作为目标函数时，坡面流汇流系数CQOF的敏感性指数高于其它参数，说明该参数对洪量有较大影响。

3.2 模型率定

根据敏感性分析结论，NAM模型模拟结果对不同参数变化的敏感性是有差异的，参数CQOF和CK1、2对NAM模型精度的影响比其它参数大，率定时需重点考虑。因此本研究在对2020年6月8日—8月8日降雨进行模拟时，先采用NAM模型自带的SCE全局寻优算法对8个参数进行自动率定，迭代次数设置为1000次，在此基础上对CQOF和CK1、2两个重点参数进行人工率定，最终得到的参数率定成果如表3所示。桃溪站模拟洪水与实测洪水拟合效果见图3，洪量的拟合效果见图4。

表3 桃溪站参数率定成果

图3 桃溪站模拟洪水与实测洪水拟合效果

图4 桃溪站累计洪量拟合效果

2020年桃溪站实测洪峰流量为1239m3/s，模型计算得到的洪峰流量为1304m3/s，差值为5.2%，原因在于2020年桃溪站上游柏林圩等圩口发生了溃破，对丰乐河洪峰流量起到了削峰作用。2020年桃溪站实测最大7d洪量为4.26亿m3/s，模型计算最大7d洪量为4.41亿m3，误差为3.5%，桃溪站实测洪峰与模型模拟洪峰的出现时间相差10h，误差在合理范围内。考虑破圩还原后，误差进一步减小，计算误差结果见表4。

表4 计算误差统计表

4 结论

(1)基于Sobol方法对MIKE NAM模型中的主要参数进行全局敏感性分析，分析结果表明：参数CQOF和CK1、2为高敏感度参数，对模型精度影响较大，建议作为重点率定参数。

(2)对CQOF和CK1、2两个参数采用人工率定结果，其它参数采用模型自动率定结果，以丰乐河桃溪水文站控制流域为研究区，模拟了2020年汛期降雨作用下的区域降雨径流过程，并根据实测流量数据对模拟效果进行拟合检验，结果表明模型计算误差在合理范围内。

(3)本文构建的NAM模型中未考虑流域内圩区的调蓄以及流域下游水位实时变化的作用，对模拟精度有一定的影响，并且参数敏感性分析成果缺少不同区域不同场次降雨的验证。后续可以将NAM模型与MIKE11水动力模型进行耦合，针对不同降雨强度、不同流域特征进行建模和率定，进而实现巢湖全流域的降雨径流模拟。本研究可为巢湖流域水文水动力模型的构建和验证提供技术支撑和应用基础。