矿坑涌水量预测模型参数灵敏度分析

王科奇, 刘佩贵, 陶月赞, 李 飞, 樊慧慧

(1.合肥工业大学 土木与水利工程学院,安徽 合肥 230009; 2.合肥工业大学 资源与环境工程学院,安徽 合肥 230009)

矿坑涌水量预测模型参数灵敏度分析

王科奇1, 刘佩贵1, 陶月赞1, 李 飞1, 樊慧慧2

(1.合肥工业大学 土木与水利工程学院,安徽 合肥 230009; 2.合肥工业大学 资源与环境工程学院,安徽 合肥 230009)

矿区水文地质条件一般较复杂,参数的空间变异性强,致使矿坑涌水量预测模型中的参数具有一定的不确定性,从而影响预测结果的可靠度。为定量确定影响预测结果的主要因子,文章以某金属矿区的矿坑涌水量计算模型为研究对象,运用扩展傅里叶幅度灵敏度检验法(extended Fourier amplitude sensitivity test,EFAST)对模型参数进行了全局灵敏度分析,得到相关参数的一阶灵敏度和总灵敏度,并分别与局部灵敏度分析法和Morris法的评价结果进行了详细对比分析。研究结果表明,采用EFAST法开展矿坑涌水量预测模型的参数灵敏度分析是切实可行的,也是较为有效的方法之一,矿坑涌水量预测结果对渗透系数最为敏感。

矿坑涌水量;渗透系数;灵敏度分析;扩展傅里叶幅度灵敏度检验法(EFAST);Morris法

矿坑涌水量的准确预测是矿区水文地质工作的重要内容之一,其预测方法主要有解析法、数值法、水文地质比拟法等,这些方法均要借助数学模型计算涌水量,而影响其计算精度的关键要素之一是参数的合理取值。然而受成矿条件等因素的影响,大多数矿区水文地质条件较为复杂,含水层空间变异性强,此外还受到勘察程度、认知水平等条件的限制,矿坑涌水量预测模型中的相关参数往往存在很大的不确定性,给参数的取值造成了一定的困难,也严重影响了涌水量预测结果的可靠度。研究预测结果对模型中重要参数的响应程度(即参数的灵敏度分析),是解决该问题的有效途径之一。

灵敏度分析方法可以定性或者定量分析模型输入参数的不确定性对模型结果产生的影响程度,合理确定主要与次要因子。灵敏度分析可分为局部灵敏度分析(local sensitivity analysis)和全局灵敏度分析[1](global sensitivity analysis)2类。其中,局部灵敏性分析因具有操作简便、易于理解等优点,在地下水数学模型中应用较广泛[2-4],但该方法未考虑参数间的耦合作用,其分析结果有一定的片面性。全局灵敏度分析则突破了该局限,且可以提供更加丰富的信息,分析方法主要有回归法[5]、Morris法[6]、Sobol’法[7]、傅里叶幅度灵敏度检验法[8](Fourier amplitude sensitivity test,FAST)以及扩展傅里叶幅度灵敏度检验法[9](extended FAST,EFAST)等。其中,Morris法具有计算简单、计算成本较低的优点,在地下水数学模型分析方面应用较多[10-12],但该方法是一种定性全局灵敏度分析方法,不能对参数的影响大小进行定量化,且参数的初始值和扰动值会对结果产生较大影响[11,13]。而EFAST法是一种基于方差分析的全局灵敏度分析方法,能够定量地获得参数的灵敏度,与其他方法相比可获得的信息多,相对较优。因此,本文以丘陵区某金属矿山的矿坑涌水量预测模型为例,选取模型中的相关参数,采用EFAST法分析水量预测结果对各个参数的敏感程度,以期为提高预测结果的可靠度提供科学依据。

1 傅里叶幅度灵敏度扩展检验法

EFAST法是一种基于方差分析的灵敏度分析方法,认为参数的灵敏度可以用模型结果的方差来反映,模型结果的方差由各参数及参数间的相互作用所导致[14]。该法既具有FAST方法计算参数一阶灵敏度的高效性,又与Sobol’法一样有计算参数的总灵敏度的能力。

假定模型为Y=f(X),模型的输入参数为X(x1,x2,…,xn),选择适当的搜索函数Gi将模型Y=f(x1,x2,…,xn)转化为Y=f(s),即实现了多维到一维的转换。取

xi(s)=G(sinwis),i=1,2，…,n

(1)

其中,s为标量,且-∞

通过对f(s)进行傅里叶级数展开可得到:

y=f(s)=∑|Ajcosjs+Bjcosjs|

(2)

其中, +∞

由参数xi的不确定性所引起的模型输出结果方差Vi为wi的整数倍振幅平方和[15]，即

(3)

其中,Z0=Z-{0}。

模型输出结果的总方差V为所有频率的振幅平方和:

(4)

参数xi的一阶全局灵敏度Si定义为:

Si=Vi/V

(5)

一阶全局灵敏度能够反映参数自身的不确定性对模型结果的直接影响程度。

模型输出结果的方差V分解[16]为:

(6)

对模型结果的各高阶方差进行归一化处理可以得到参数的各阶灵敏度,例如:

Sij=Vij/V,S12…n=V12…n/V

(7)

针对一个多参数的非线性模型,参数xi的总全局灵敏度为参数各阶灵敏度的总和[7],即

STi=Si+Sij+Sijm+…+S12…i…n

(8)

参数的总灵敏度STi不仅包含了参数xi自身的不确定性对模型输出结果方差的影响,同时也包含了xi与其余参数间的耦合作用对模型结果的影响。参数xi的总灵敏度STi与一阶灵敏度Si的差值STi-Si可以表示该参数与其他参数之间相互作用对模型结果的间接影响;若差值较小,则说明参数间的交互作用不明显[17]。

2 实例研究

2.1 研究区概况

研究区为低山丘陵地形,属亚热带季风气候区,多年平均降水量1 216.2 mm,年平均蒸发量1 497.5 mm。矿区含水岩组主要有松散岩类孔隙含水岩组,火山碎屑岩、熔岩类空洞-裂隙含水岩组,碎屑岩、碳酸盐岩类岩溶-裂隙含水岩组及岩浆岩类裂隙含水岩组。其中,碎屑岩、碳酸盐岩类岩溶-裂隙含水岩组为矿床的含矿层位,矿床内各含水岩组间无明显的隔水层存在,天然状态下有一定的水力联系,具有统一的水动力场。降雨入渗为该矿床的主要充水因素,由风化裂隙向深部补给。

2.2 参数灵敏度分析

矿体含水层为裂隙含水岩组,水力性质为承压水,并且无明显的侧向隔水介质,侧向水力联系较好;在开采过程中,随着排水疏干,地下水位会降到含水层顶板以下,水力性质转为承压-潜水水流。因此选用均质无限补给边界条件下承压转无压完整井的稳定流涌水量计算公式,即

R0=R+r0,S=H-h,

(9)

其中,Q为矿坑涌水量预测值;K为渗透系数;H为抽水前承压含水层从底板起算的水头值;M为承压含水层厚度;h为抽水后承压含水层从底板起算水头值;S为水位降深;R0为矿井开采区引用影响半径;r0为矿井开采区引用半径;R为单井影响半径;F为预计矿井开采面积。

根据勘察成果及矿区开展的相关水文地质工作,确定参数取值范围如下：K为0.42～0.65 m/d，M为419.56～546.43 m，H为557.86～632.39 m，r0为313.37～368.67 m。因资料有限,假定参数均服从均匀分布。

目前对于地下水数学模型,通常采用局部灵敏度法进行灵敏度分析。本文首先采用差分法对矿坑涌水量预测模型进行局部灵敏度分析,以各参数的平均值为基准值,通过逐一改变参数的变化幅度来对模型进行局部灵敏度分析,结果如图1所示。

图1 局部灵敏度系数法评价结果

由图1可知,通过局部灵敏度分析法得到的参数灵敏度排序依次为|H|>|K|>|r0|>|M|。

本文借助Simlab软件实现样本的抽样和结果的全局灵敏度分析,根据EFAST法的基本原理每个参数至少要采样65次,本次试验设置样本数为1 000。利用软件的内部函数编辑器输入该数学模型,直接调用抽样结果即可获得模型的输出结果。对模型结果进行灵敏度分析,可得各输入参数的一阶灵敏度和总灵敏度,见表1所列。

表1 参数的一阶灵敏度和总灵敏度

对比上述2种方法的分析结果可知,EFAST法得到的参数灵敏度排序与局部灵敏度分析法并不一致。这是因为局部灵敏度分析方法是在线性模型的基础上发展而来的,对于非线性模型，该方法只能分析因子对模型局部的影响。全局灵敏度分析方法能够分析参数的变化范围对模型结果的影响,与局部灵敏度分析方法相比,对于非线性模型其结果更趋于实际,这也正说明在非线性模型灵敏度分析中采用全局方法的重要性。

为了与EFAST法作对照，本文采用在地下水数学模型中应用较为广泛的全局灵敏度分析法Morris法[10],利用Simlab软件对该数学模型进行灵敏度分析,获得各参数的灵敏度分析结果,如图2所示。图2中,σ为各参数单位效应的标准差;μ*为各参数单位效应的均值。

图2 基于Morris法的灵敏度分析结果

2.3 结果分析与讨论

在局部灵敏度分析法中，某一特定参数k的灵敏度系数S是通过只改变参数k的值,使其由ak变化为ak+Δak,相应模型结果由yi(ak)变化为yi(ak+Δak),由(10)式获得S的近似值,即

(10)

由(10)式可知局部灵敏度分析法能够在线性模型中获得正确的结果,而对于非线性模型,参数在基准值处变动过大时就不能提供有效分析结果,这正是局部灵敏度分析法获得的参数灵敏度排序与EFAST法不同的原因。

利用EFAST法能够定量获得各参数的一阶灵敏度及总灵敏度,由表1数据可知,参数的总灵敏度与一阶灵敏度之间的差值较小,说明参数之间的耦合作用对结果的影响程度较小;渗透系数K、水头高度H及含水层厚度M这3个参数的灵敏度相对较大,矿坑涌水量预测值对这3个参数较为敏感;其中K的不确定性对矿坑涌水量影响最大,一阶灵敏度高达0.772 1。

Morris法能够获得各参数单位效应的均值μ*和标准差σ,其中μ*可以表征参数灵敏度的相对大小,值越大,说明参数的敏感性越强;σ表征参数间的相互作用程度[18]。由图2可得参数的灵敏度大小依次为|K|>|H|>|M|>|r0|。EFAST法获得的一阶灵敏度的排序与Morris法的排序结果是一致的,间接证明了EFAST法分析结果的可靠性,也表明应用该方法在矿坑涌水量预测模型参数灵敏度分析方面的可行性。

3 结 论

本文基于EFAST法,以某矿区的矿坑涌水量计算模型为例,进行了参数的灵敏度分析,分析结果表明矿区含水层的渗透系数对水量预测结果影响最大,即表明渗透系数最为敏感。为了减小矿坑涌水量预测的不确定性进而获得更高精度的的预测值,在勘察工作中应尽量提高渗透系数的精度。

对于非线性模型,全局灵敏度分析结果比局部灵敏度的分析结果更趋于实际,EFAST法是定量刻画全局灵敏度的分析方法之一,不仅克服了局部灵敏度分析方法的片面性,而且能将参数的敏感程度定量化,而不是像Morris法仅能获得参数的灵敏性排序,提供了更综合全面的信息。

本文选取的是解析模型,后续研究将进行多参数数值模型的定量全局灵敏度分析。

[1] 雷晓辉,田雨,贾仰文,等.WEP 模型全局参数敏感性分析及其在汉江上游流域的应用[J].水文,2010,30(6):14-18,41.

[2] 束龙仓,朱元生，孙庆义，等.地下水允许开采量确定的风险分析[J].水利学报,2000,31(3):77-81.

[3] 束龙仓,朱元生,孙庆义,等.地下水资源评价结果的可靠性探讨 [J].水科学进展,2000,11(1):21-24.

[4] 李森,陈家军,叶慧海,等.地下水流数值模拟中随机因素的灵敏度分析[J].水利学报,2006,37(8):977-984.

[5] MCKAY M D,BECKMAN R J,CONOVER W J.Comparison of three methods for selecting values of input variables in the analysis of output from a computer code[J].Technometrics,1979,21(2):239-245.

[6] MORRIS M D.Factorial sampling plans for preliminary computational experiments[J].Technometrics,1991,33(2):161-174.

[7] SOBOL’ I M.Sensitivity analysis for nonlinear mathematical models[J].Mathematical Modeling and Computational Experiment,1993,1(4):407-414.

[8] CUKIER R I,FORTUIN C M,SHULER K E,et al.Study of the sensitivity of coupled reaction systems to uncertainties in rate coefficients：Ⅰ. theory[J].The Journal of Chemical Physics,1973,59(8):3873-3878.

[9] SALTELLI A,TARANTOLA S,CHAN K P S.A quantitative model-independent method for global sensitivity analysis of model output[J].Technometrics,1999,41(1):39-56.

[10] 束龙仓,王茂枚,刘瑞国,等.地下水数值模拟中的参数灵敏度分析[J].河海大学学报(自然科学版),2007,35(5):491-495.

[11] 鲁程鹏,束龙仓,刘丽红,等.基于灵敏度分析的地下水数值模拟精度适应性评价[J].河海大学学报(自然科学版),2010,38(1):26-30.

[12] 刘瑶林,刘国东,徐涛,等.多层含水层地下水数值模型参数灵敏度分析[J].环境科学与技术,2014,37(12):33-37,44.

[13] 张丽华,苏小四,孟祥菲,等.地下水流数值模拟参数全局灵敏度分析[J].中国农村水利水电,2014(8):92-97.

[14] 吴锦,余福水,陈仲新,等.基于 EPIC 模型的冬小麦生长模拟参数全局敏感性分析[J].农业工程学报,2009,25(7):136-142.

[15] 王建栋,郭维栋,李红祺.拓展傅里叶幅度敏感性检验 (EFAST) 在陆面过程模式中参数敏感性分析的应用探索[J].物理学报,2013,62(5):35-41.

[16] SALTELLI A,BOLADO R.An alternative way to compute Fourier amplitude sensitivity test (FAST)[J].Computational Statistics & Data Analysis,1998,26(4):445-460.

[17] 余衍然,李成,姚林泉,等.基于傅里叶幅值检验扩展法的轨道车辆垂向模型全局灵敏度分析[J].振动与冲击,2014,33(6):77-81.

[18] 宋明丹,冯浩,李正鹏,等.基于 Morris 和 EFAST 的 CERES-Wheat 模型敏感性分析[J].农业机械学报,2014,45(10):124-131,166.

(责任编辑 张淑艳)

Sensitivity analysis for forecasting model parameters of mine water inflow

WANG Keqi1, LIU Peigui1, TAO Yuezan1, LI Fei1, FAN Huihui2

(1.School of Civil and Hydraulic Engineering, Hefei University of Technology, Hefei 230009, China; 2.School of Resources and Environmental Engineering, Hefei University of Technology, Hefei 230009, China)

The uncertainties of mine water inflow forecasting model parameters have often been witnessed due to the complex hydrogeological conditions and high spatial variation of parameters of mining area, thus affecting the reliability of predictions. In order to quantitatively determine the main factors influencing the size of the water flow, taking the calculation model of pit water inflow of a metal mining area as the research object, the global sensitivity analysis on model parameters with the application of extended Fourier amplitude sensitivity test(EFAST) method is utilized. The first-order sensitivity of relevant parameters and the total sensitivity are obtained. In addition, the results are also compared with the parameter sensitivity sequencing obtained by the Morris method and local sensitivity analysis. The results show that the EFAST method is practical and feasible for the parameter sensitivity analysis of mine water inflow prediction model and it is one of the most effective methods, and the prediction result is most sensitive to the hydraulic conductivity.

mine water inflow; hydraulic conductivity; sensitivity analysis; extended Fourier amplitude sensitivity test(EFAST) method; Morris method

2015-09-25；

2016-02-16

国家自然科学基金资助项目(51309071；51509064)

王科奇(1991-),男,河南南阳人,合肥工业大学硕士生; 陶月赞(1964-),男,安徽巢湖人,博士,合肥工业大学教授,博士生导师.

10.3969/j.issn.1003-5060.2017.04.014

P641.41

A

1003-5060(2017)04-0502-04