刘凌雪，敖天其,2，黎小东，胡 正，胡富昶，李孟芮

(1.四川大学水利水电学院, 成都 610065; 2.四川大学水力学与山区河流开发保护国家重点实验室, 成都 610065)

随着水文气象资料和下垫面信息的愈加丰富及计算机科学的发展，分布式水文物理模型成为当今水文科学研究的前沿方向之一[1]。与传统的集总式模型相比，分布式水文模型在拥有众多优点的同时也面临着因参数过多而难以率定的问题。对于足够复杂且具有高维度和高计算需求的分布式水文模型而言，伴随着参数相互作用的过度参数化导致模型参数不是唯一可识别的，Beven[2]将此现象称为“异参同效”现象，即不同的参数集会导致相同或相似的模拟结果。避免过度参数化或不可识别性的一种可能的方式是减少参数个数[3]。

模型敏感性分析(SA)是过滤模型关键参数的有效工具[4]，它旨在从不敏感的参数中识别出最敏感的参数以减小参数维数[5]，降低参数率定的不确定性，提高模型的参数优化效率。目前，SA方法已广泛用于许多领域的实际研究中：Collins和Avissar[6]采用傅立叶振幅灵敏度测试(FAST)方法评估LAID模型中参数对显热和潜热的影响情况；Werkhoven等[7]将全局SA方法作为筛选工具来降低参数维数，进行参数的多目标优化率定，并将其成功应用于SAC-SMA模型；Sun等[8]考虑了多目标敏感性分析和多目标条件的参数优化，将Morris筛选法与非支配排序差分进化(NSDE)算法耦合并用于MIKE/NAM这一降雨径流模型的参数敏感性分析和率定中。

BTOPMC模型[9]属分布式水文物理模型，有5个可以定量反映流域物理特征的参数，主要涉及每个网格单元的土壤类型，植被和土地利用(土地覆盖)。尽管该模型的参数具有明确的物理意义，但由于测量困难，参数值往往通过参数率定来确定。基于BTOPMC模型参数的物理意义及以往的模拟经验，在不同的研究区，由于流域物理特征及类型的不同，模型参数需要分为不同类型的土壤、植被、地形以及地表覆盖情况来率定[10]，这就使得BTOPMC的参数率定过程通常需要较长时间，且由于“异参同效”现象使得参数最优值也存在较大的不确定性。本文分别采用定性分析方法MOAT[11]和定量分析方法DGSM[12]对BTOPMC模型的参数进行敏感性分析，筛选较敏感参数以减少需要被率定的模型参数个数，从而达到减少参数率定值的不确定性、提高参数率定效率的目的。

1 BTOPMC模型

BTOPMC是由日本山梨大学[9]在TOPMODEL的概念基础上加入了洼地消除、数字河网生成、马斯京根-康奇法等子模型而开发的一个分布式水文物理模型。BTOPMC模型的结构组成见图1[13]。目前，该模型在河网技术、自然子流域划分、产汇流模型参数以网格为单位的空间分布给定方法、汇流精度等方面取得了明显进步，涉及的主要内容有[10]：①发现了按数字高程模型(DEM)中的洼地的空间位置计算高程微增量，从而消除洼地的自动解析新技术；②研究了Pfafstertter流域自动分割法用于分布式水文模型的问题并提出改进措施；③提出了用Pfafstertter法分割流域时的稳定分割程度的判别方法；④开发了解决马斯京跟汇流方法中的负出流问题的近似处理方法和时空步长实时调节方法；⑤提出了将TOPMODEL用作大流域的产流模型时，以各网格(而不是由网格组成的子流域)中土壤、植被等流域物理特性给定个模型参数的新模式；⑥提出了用Pareto解析法减少模型参数优化过程及建立参数转移函数过程中的不确定性的必要性。

图1 BTOPMC模型结构组成Fig.1 The structure of BTOPMC model

BTOPMC模型具有5个需要被率定的参数。基于BTOPMC模型参数的物理意义及以往的模拟经验，在不同的研究区，由于流域物理特征及类型的不同，模型参数需要分为不同类型的土壤、植被、地形以及地表覆盖情况来率定[10]。根据所使用的数据源[14]，本研究中模型参数的具体分类见表1，其中“标识”列中的Xi为对应于不同参数的在本研究中需要被率定的变量。

表1 BTOPMC模型参数表[10]Tab.1 Model parameters of BTOPMC[10]

BTOPMC模型中用于模型精度评价的指标为Nash-Sutcliffe效率(NSE)[14]。NSE的数值越大，模拟结果越好，其计算公式如下：

NSE=

(1)

式中：N为总的时间步长；Qsim和Qobs分别为时间步长t所对应的模拟流量和实测流量；Qav是整个时期观测流量的平均值。

2 研究区域及资料

富士川流域位于日本本州岛中央，是日本三大险峻河流之一，属一级河川，研究区范围示意图见图2。河流主干流长度约128 km，流域控制面积约786.2 km2。全流域大概88%被森林覆盖，流域内最高点是著名的富士山。全流域平均降雨量约1 550 mm，年平均气温在4～18 ℃之间。年内，4月的降雨较3月有明显增多；6月，气压配置成为“梅雨型”，直到7月上、中旬梅雨结束。当梅雨气压停滞日本上空时，若有南方热带低气压出现，舌状的暖气流就会进入日本上空，“湿舌”与“梅雨”相遇便会导致流域内出现暴雨天气。8月末、9月上旬至9月末、10月上旬是易降雨的季节，而且这一时期也是台风多发期，降大雨的可能性很大。

图2 研究区范围示意图Fig.2 Schematic diagram of the study area

本文选取富士川流域1991年9月的一场降雨过程，7个测站的降雨资料均来自气象资料自动集取系统(AMeDAS)的每小时降雨观测记录，基本无数据丢失。年潜在蒸散发量采用Penman-Monteith方法估计为1 800 mm。对于洪水的模拟，潜在蒸散发量被假定为5.5 mm/h。流域的植被覆盖图来自IGBP第二版(USGS)，土壤图来自FAO，且植被和土壤的分布情况是对每一个网格给出的[14]。

3 研究方法

近年来，基于实验设计(DoE)的全局SA方法十分受欢迎，原因在于它们在保持计算效率的同时还可定量表征全局敏感度。典型的基于DoE的SA过程主要有两个步骤：首先，使用所选择的取样方法在可行参数空间内生成一组样本参数；然后，利用SA方法获得不同参数变化对模型输出的影响情况。目前，有许多常用的DoE方法，如蒙特卡洛(MC)、拉丁超立方体(LH)[15]、对称拉丁超立方体(SLH)等。

3.1 取样方法

本文初步选取了4种DoE方法：MC是一种最常用的DoE方法，它依赖于重复随机采样来获得未知概率实体的指定分布函数的近似数值，其敛散性依赖于独立的随机参数的个数；准蒙特卡洛(QMC)是确定性的MC方法，也被称为低偏差序列，它是一组填充样本区域的“有效”点[16]；LH是一种根据多维分布生成合理的参数集合分布的统计方法，它适合于任意维数的抽样设计；SLH与LH均具有“充满空间”性质，理论上来讲，LH的试验点带有随机性，抽样表现不稳定，SLH则能够产生分布性更好的试验点。

对于BTOPMC，我们使用300个采样点来评估这4种方法的空间填充属性，分别重复运行MC、LH和SLH 6次，因QMC是一个确定性序列所以只运行一次，计算空间填充性能好坏的度量数值：修正的L2偏差(MD2)，中心化L2偏差(CD2)，对称化L2偏差(SD2)和可卷的L2偏差(WD2)。对于这些度量，数值越小，DoE算法的空间填充性能越好[17]。4种DoE算法的比较结果见图3。对比分析之后，我们选择LH作为后续BTOPMC敏感性分析的DoE方法。

3.2 SA方法

Morris One At A Time(MOAT)是一种简单、高效的参数筛选方法，通过少量的模型计算就可以获得模型参数的定性排序，尤其适合参数较多的复杂模型[11]。该方法虽然可能将不重要参数判断为重要参数，但不会出现将重要参数判为不重要参数的错误[18]。每个参数的总体效果可以利用所采取样本的对应参数梯度的均值μ和标准差δ来近似估计(梯度是指目标函数值的差值与样本参数值的差值之比)。均值μ的大小表征参数的敏感度，而标准差δ表征参数之间相互作用的程度。Campolongo等[19]提出了修正的均值μ*，以解决部分参数对均值μ存在负效应的问题，公式如下：

(2)

di=[f(X1,…,Xi-1,Xi+Δ,Xi+1,…,Xn)-f(X)]/Δ

(3)

其中X=(X1,X2,…,Xn)是一个试验范围内的任意选定值，是一个n维p级的网格，Xi可以从{0,1/(p-1),2/(p-1),…,1}中取值；Δ是1/(p-1)预定的倍数，当p为偶数时，Δ通常取p/[2×(p-1)]。

对于MOAT方法而言，μ*越大，参数越敏感。本文中，我们使用修正的平均值μ*作为MOAT方法的敏感性度量。此外，MOAT[11]有特定的取样方法―Morris，所取样本大小通常设置为n+1的倍数，其中n是参数个数。由表1可得，本研究中n取19。

图3 四种DoE方法对应的度量值比较Fig.3 Comparison of metrics for the four DoE methods

DGSM[12]是一种新的基于导数的全局SA方法，该方法隶属于微分法，它是建立在求解响应函数对输入变量在多个点处的偏导数的基础上，可以看作是局部敏感性在全局范围内的扩展。对于涉及参数较多的高维模型而言，DGSM进行数值评估所需的计算时间比广泛应用的Sobol敏感度指数的计算时间低许多个数量级。在本研究中，我们将计算所得的各参数敏感性指数占所有参数敏感性指数的比例作为该方法的敏感性度量。各参数敏感性指数的确定公式如下：

STi=vi/π2D

(4)

4 结果与讨论

4.1 定性敏感性分析

根据第3.2节中关于Morris取样方法的介绍，取300个参数样本对BTOPMC进行定性参数筛选，结果见图4(a)，可以发现X8，X9，XC，X7，X6，X4，XB，X3，X5，XF，XE，XA的μ*值(每个Xi的箱线图中红色线的对应值)明显大于0，可以确定为敏感变量，涉及的参数类型有m(m)，Srmax(m)，Sbar0(m)和n0。此外，LH抽样也可用于MOAT进行定性SA分析，样本数量同样取300，分析结果见图4(b)，我们可以观察到的敏感变量有X9，X8，X6，XE，X7，X3，XF，涉及到参数的类型有m(m)，Srmax(m)和n0。综合分析，两种取样方法所对应的MOAT结果是具有一定程度上的一致性的，由于采样不确定性，它们之间所存在的差异是合乎情理的。

根据MOAT方法虽然可能将不重要参数判断为重要参数，但不会出现将重要参数判为不重要参数的错误[18]的这一分析特征，将m(m)，Srmax(m)和n0确定为BTOPMC的敏感参数类型，其中Srmax(m)中所包含的X6，X7，X8，X9和m(m)中所包含的X3，X4，X5均可被确定为需要被率定的敏感变量，n0中XE和XF的敏感性较为突出，同样被认定为敏感变量。

图4 BTOPMC模型参数敏感性的定性分析结果Fig.4 Qualitative sensitively analysis results of the parameters in the BTOPMC model

4.2 定量敏感性分析

为了验证MOAT方法对敏感参数的筛选结果，并量化每个敏感性参数的贡献率，我们利用DGSM这一种定量SA方法对BTOPMC的参数敏感性进行量化分析。本研究将LH作为DGSM的取样方法。由于定量SA方法通常比定性SA方法需要更多的样本，所以本文将样本数量取为500，分析结果分别见图5。可明显看出，X9，X8，X7，X6的敏感性较强，此外，X3，X4.X5,XE，XF的总效应也很突出，对应到参数类别，我们可以将Srmax(m)，m(m)和n0作为BTOPMC模型的敏感参数类型。此外，虽然XA，XB，XC，XD对于SA结果有一定程度上的影响，但相比较而言影响较小，在此不将其作为敏感性参数考虑。

比较由MOAT方法所得的定性分析结果和由DGSM方法所得的定量分析结果，可以发现所筛选的BTOPMC的敏感参数基本一致。由于采样方法及分析方法的不同，参数敏感程度的大小排序存在差异也是合乎情理的[20]。

图5 BTOPMC模型参数敏感性的定量分析结果Fig.5 Quantitative sensitively analysis results of the parameters in the BTOPMC model

4.3 SA结果的物理解释

BTOPMC中，T0和m是反映流域土壤类型影响的参数类型。作为T0的衰减因子，m的变化将导致T0明显变化[10]，所以，当m和T0发生相同的变化时，由m所引起的模型性能的变化是更加明显的。也就是说，在BTOPMC中m比T0更敏感。从表1可知，m中的X3，X4和X5分别是T0中X0，X1和X2的衰减因子，相应地，X3，X4和X5更为敏感。

Srmax主要反映了植被和土地利用对模拟结果的影响情况。在根区的储水值达到Srmax之前，水将在重力的作用下渗入不饱和区域。富士川流域面积的88%被植被覆盖，Srmax变化对BTOPMC在该流域的模拟结果影响较大是合理的。此外，与深根植被区相比，浅根植被区的Srmax对地表径流的影响更大，而灌溉区和不透水区的Srmax对模型模拟结果的影响最大，对应到相应参数，X9和X8比X7更敏感，X6的敏感性则相对较小。

n0即河道对水流的摩擦阻力，是BTOPMC汇流模型中唯一需要被率定的参数，受多种因素影响，如河岸或河床材料、河道结构、河岸不规则程度、植被覆盖情况等。因此，汇流模型通常对n0的变化比较敏感。此外，黏土、砂土和粉土的n0基本相同，且都小于河网和建筑区，因此XE和XF的变化对模型性能的影响要大于XG，XH和XI。

Sbar0是反映流域地形对模型模拟结果影响情况的参数类型，需要针对每个子流域进行率定。与针对每个网格进行率定的前4种参数类型相比，大面积均质化可以解释其在BTOPMC中的相对不敏感性。

当然，由于SA方法和样本选取方法的局限性、观测误差、模型本身的结构不完善等原因，使得并非所有筛选结果都可以有精确的物理解释[21]。

5 结 论

在本次研究中，采用MOAT这一定性SA方法和DGSM这一定量SA方法分别对BTOPMC模型的参数敏感性进行定性和定量分析，将Srmax(m)，m(m)和n0确定为BTOPMC模型用于富士川流域时的敏感参数。

本文的敏感参数筛选结果与模型参数的物理解释一致，这也证明了采用全局SA方法寻找复杂模型中重要参数的可行性。

虽然在富士川流域通过三种方法筛选出了敏感参数，但还需要在此基础上对该模型参数的“异参同效”现象进行分析，并进行参数率定以进一步明确该研究结果对BTOPMC参数优化效率的影响情况。此外，通过减少参数个数来降低BTOPMC模型参数率定结果的不确定性和提高参数率定效率的普遍性需要在更多具有不同特征的流域上得到进一步研究。