基于LH-OAT的马斯京根模型参数全局敏感性分析

2017-03-21刘春烨

中国农村水利水电 2017年9期

刘春烨

(太原理工大学水利科学与工程学院，太原 030024)

马斯京根法是河道洪水流量演算中广泛应用的一种水文学方法，该法在建立马斯京根槽蓄方程的基础上，与水量平衡方程联立，推出上下游流量关系方程，进行河段洪水的演算[1]。马斯京根法应用的一个前提是模型参数的合理率定，通常对参数K、x进行率定，再由K、x率定值计算参数C0、C1、C2。参数K、x的率定方法有分析法和试算法，但这2种方法较为复杂且率定结果精确性较差。因此，部分学者采用了直接优化参数C0、C1、C2的方法，主要有POS算法[1]、遗传算法、最小二乘法、SCE-UA算法[2]等。

对模型参数进行敏感性分析是参数率定的基础。敏感性分析即分析参数对模型的影响程度，根据参数的影响程度，筛选出重要的影响参数，减少参数的率定次数和调参过程中的不确定性，提高率定效率。水文模型敏感性分析方法有局部敏感性分析法、全局抽样敏感性分析法。局部敏感性分析法原理简单但无法完全考虑参数间的相互影响，不能满足提高模型精确性的要求。全局敏感性分析法可综合分析多个参数的变化对模型运行结果的影响，且该方法已在各种水文模型中得到了应用，如RSA方法、GLUE方法[3]和LH-OAT方法[4]等。目前，国内外直接对马斯京根模型参数C0、C1、C2敏感性分析较少。本研究引入Latin-Hypercube抽样算法，并结合OAT 方法,采用LH-OAT全局敏感性分析方法直接对马斯京根模型参数C0、C1、C2进行全局敏感性分析。在敏感性分析的基础上，确定参数的敏感性大小顺序，采用假定最优参数法进行马斯京根模型参数率定，发现按敏感性顺序进行参数的率定，可提高马斯京根模型参数的率定效率。

1 LH-OAT全局敏感性分析方法

现阶段处理参数敏感性的方法有随机法、模糊法、区间分析法。前2种方法需要较多的数据，而区间估计法只需已知参数的上下界，方法较为简单。LH-OAT全局敏感性分析方法属于区间估计法，是采用无量纲参数来表示全局敏感性的一种方法，该法的适用条件为已知模型参数的取值范围。LH-OAT全局敏感性分析方法的理论依据分别为OAT敏感性分析法和LH抽样法。OAT敏感性分析法在Monte Carlo抽样方法的基础上进行重复试验，每次对其中一个参数进行微小的扰动，保持其他参数不变，带入敏感性公式，即可得到每个参数的敏感度，具有较高的精确性。但Monte Carlo抽样法在参数取值范围内随机抽样，会产生大量的样本，使OAT敏感性分析方法的运算量变大。LH抽样方法，假设抽样符合均匀分布，将每个参数区间均分为n层，分别从每个分层中随机抽样一次，生成一个LH抽样参数组，减少了采样频次，能提高敏感性分析的效率。LH-OAT全局敏感性分析法结合了LH抽样法的高效性和OAT算法的精确性。

LH-OAT全局敏感性分析法具体步骤如下：①采用LH抽样法，将整个参数空间均分为N层，分别从每个分层中随机抽样一次，组成一个抽样参数组。②假设有P个模型参数，按照OAT方法，对每个LH 抽样参数组进行P次参数的微小扰动并且每次只扰动其中一个参数，计算每次扰动前后目标函数的变化值。③根据LH-OAT方法定义，模型总共需运行N(P+1)次。每个参数的相对敏感度由以下公式计算得到:

(1)

式中：M为敏感性分析的目标函数；ei,k为第i个参数在第k层中抽样值；Δei,k为ei,k的扰动变化值；Si,k为参数ei,k在k层的相对敏感度。

由参数相对敏感度Si,k可计算参数的全局敏感度GSi,k：

(2)

式中：N为参数空间的总层数；GSi,k为全局敏感度。

2 马斯京根模型

现阶段河道洪水流量演算有水文学和水力学2种方法，在水文学中应用最多的为马斯京根法,其公式为：

Q2=C0I2+C1I1+C2Q1

(3)

式中：Q1为计算时段初始时刻的出流量；Q2为时段末的河段出流量；I1、I2分别为计算时段始、末的河段入流量；C0、C1、C2均为流量演算系数，计算公式如下：

(4)

(5)

(6)

式中：K为蓄量常数，即蓄量流量关系曲线的坡度；x为流量比重系数；Δt为计算时段长。

马斯京根法模拟河道流量过程中不可避免地存在着不确定性，对模型参数进行全局敏感性分析是研究参数不确定性的分析方法之一。在传统的马斯京根模型参数分析过程中，是通过公式(4)、(5)、(6)来计算C0、C1、C2的值，本文提出直接对C0、C1、C2进行敏感性分析，避免了由K、x取值不准确对参数C0、C1、C2造成的影响。

3 实例分析

3.1 模型参数敏感性计算

为定性了解参数C0、C1、C2对马斯京根模型模拟洪水流量的影响程度，提出采用LH-OAT全局敏感性分析方法对模型参数进行全局敏感性分析。以文献[5]中数据为例，分析C0、C1、C2的敏感性，河段上断面入流、下断面出流情况见表1，所选河道稳定，计算时段长为 Δt=6 h。其中时段1-4为连续时段，时段5-8为连续时段。由公式(4)、(5)、(6)利用多元函数求极值的方法得到C0、C1、C2参数的取值范围[6]分别为(-9/11，3/13)，(1/4，19/21)，(1/21,1)。

表1 不同时段流量 m3/s

首先，将马斯京根模型的3个参数C0、C1、C2在其取值范围内均分成10份，共11层；然后利用Latin-Hypercube抽样法进行取样组合，可确定11种组合，且C0、C1、C2应满足公式：C0+C1+C2=1，以避免组合参数的不合理；按照OAT分析计算方法，对每种组合参数微小扰动10%(保持2个参数不变，对另一个参数进行微小扰动)后利用公式(1)进行计算，每组数据需计算11×(3+1)=44次，4个连续时段可以确定3组计算数据，有2组不同的连续时段，共需计算44×3×2=264次。利用公式(2)对相对敏感性计算结果进行均值处理，最终得到6组反映3个参数C0、C1、C2敏感性的无量纲数值，具体计算结果见表2。从表1中随机抽取一组数据，分别把参数C0、C1、C2的扰动幅度变为20%、30%，每改变扰动幅度一次需计算44次，变2次扰动幅度，故需计算88次，分析变幅后的敏感性，计算结果见表3。

表2 参数的全局敏感度Tab.2 Global sensitivity of parameters

表3 不同扰动幅度参数敏感度Tab.3 Sensitivity of different disturbance amplitude parameters

3.2 参数敏感性结果分析

以径流量作为敏感性评价目标，从模型参数上来看，1～4时段(即河道涨水阶段)马斯京根模型参数针对不同流量会产生相应波动，但参数敏感性总趋势为C0>C1>C2；5-8时段(即河道落水阶段)参数敏感性趋势为C0VC1>VC2。

图1 不同时段参数敏感性比较Fig.1 Comparison of parameter sensitivity

图2 不同扰动幅度参数敏感度趋势Fig.2 Sensitivity trend of different disturbance amplitude parameters

以洪峰流量作为敏感性评价目标，当出现洪峰时(3～4时段)，C0、C1、C2虽遵循涨水阶段的规律，但C0、C1的敏感性相对前一时段的值减小，C2的敏感性相对增大，3个参数的敏感性相差不大，无主要敏感性参数。分析敏感性变化的原因，是由于洪峰流量较河道普通流量大，当流量过大时，河道的稳定性减弱，影响因素相应发生变化且影响因素增加。

采用假定最优参数法分析参数敏感性对率定效率和率定结果的影响，首先选定一组参数作为最优参数，分别采用2种方法进行参数率定，第1种方法是根据参数的敏感性大小顺序进行参数率定，第2种方法是根据扰乱参数敏感性顺序进行率定。以率定结果与最优参数的偏离程度、率定效率来评价2种方法的优劣，结果见表4、图3。方法1率定11次，精度为99.8%，方法2率定16次，精度为99.6%。因此，按敏感性由大到小进行参数率定可提高69%的率定效率。以洪峰流量为率定目标时，由于没有主要敏感性参数，参数的敏感性差别不大，率定效率会有所下降。