概率积分法模型参数全局灵敏度分析*

王 宁，朱 健

(江苏省电力设计院，江苏 南京 211102)

0 引言

中国是个煤炭大国，煤炭资源的开采在现在乃至未来的几十年仍然在国民经济发展中占据主要的地位，且作用极其重大。煤炭资源的开采势必引起岩层的破坏波及地表，对地面建(构)筑物等方面都产生不同程度的影响。国内学者刘宝琛、廖国华从随机介质理论中发展出的概率积分法模型［1］预测的地表变形完全可以满足工程应用，用该模型预测采煤引起的下沉值相对中误差可以控制在10%以内，而且模型参数物理意义明确，因此该模型是国内煤炭行业应用最为广泛而且最为成熟的模型。采用此模型来研究开采沉陷引起的地表移动破坏规律已经取得了丰富的成果。在其发展过程中，很多学者研究了模型中各个参数的取值方法以及与地质采矿条件的数学关系，然而模型本身各个参数的灵敏度情况却鲜有研究，周大伟［2］利用正交设计统计分析过概率积分法参数的敏感性。

灵敏度分析目的在于确定模型中哪些参数最容易在模型描述过程中引入不确定性［3，4］。通过灵敏度分析可以确定模型各参数对输出结果影响的大小。本文采用全局灵敏度分析方法，分析概率积分法模型中各参数的灵敏性，为模型参数的选取提供参考数学层面的依据。对于灵敏度高的参数，在选取其取值时，应更加贴近实际情况。

1 概率积分法模型

根据开采沉陷学［1］，如图1所示。有限开采时地表移动盆地倾向主断面上的下沉预计公式为:

式中:m为采厚;q为下沉系数;α为煤层倾角;θ0为开采影响传播角;L为倾向工作面计算长度;s1、s2为下、上拐点偏移距。

图1 有限开采时地表倾向主断面移动和变形计算原理图Fig.1 Calculate diagram of movement and deformation on surface tend towards main section for finite mining conditions

2 研究区概况

以淮南矿区1212(3)工作面为例，该工作面走向长D=525 m，倾向宽L=122 m，煤层倾角10°，上边界采深H2=300 m，下边界采深H1=353 m，平均采深H0=326 m。模拟计算时，概率积分法参数的取值见表1，各参数的取值服从均匀分布。

表1 概率积分法模型参数Tab.1 Parameters of probability integral method

3 研究方法

局部灵敏度分析方法仅能反映单个参数变化对模型输出的影响。全局灵敏度分析是一种新的分析手段，对模型中所有的参数在其定义的变化范围内，通过改变参数的输入值，分析模型全部的可能输出结果。全局灵敏度分析的方法有多元逐步回归法，Morris法，Sobol法，傅里叶幅度灵敏度检验法(Fourier Amplitude Sensitivity Test)等［5－10］。其中，多元逐步回归法由于计算相对简单，易于操作而被广泛应用，本文采用该方法进行概率积分法模型参数的灵敏度分析。

多元逐步回归分析的目的是为了分析每个自变量对因变量的影响程度，运用回归分析原理采用双检验原则，逐步引入和剔除自变量而建立最优回归方程的优选方法。具体含义是:

1)每步有二个过程，即引进变量和剔除变量，且引进变量和剔除变量均需作F检验后方可继续进行，故又称为双重检验回归分析法。

2)引入变量，引入变量的原则是未引进变量中偏回归平方和最大者并经F显著性检验，若显著则引进，否则终止。

3)剔除变量，剔除原则是在引进的自变量中偏回归平方和最小者，并经F检验不显著，则剔除。

4)终止条件，即最优条件，再无显著自变量引进，也没有不显著自变量可以剔除，这也是最优回归方程的实质。

在多元逐步回归分析中，每个输入变量计算所得的系数可用于度量输出对于相应的输入分布的灵敏度。具体含义如下:参数被选入回归模型中顺序定性的表明了参数灵敏度的高低，标准回归系数(SRC)的绝对值是参数灵敏度的定量指标，SRC的符号表明了该参数对模型输出正的影响或负的影响;没有被选入回归模型的参数，说明该参数对模型输出的贡献很小，可忽略。

4 结果分析

本文采用@RISK软件进行计算，@RISK采用模拟(有时称为“蒙特卡洛模拟”)方法进行灵敏度分析。这种模拟是用一种方法让计算机反复重新计算自己的工作表，并生成可能的结果分布，每次以随机方式选择不同的值组合，获得输出值和其概率分布。事实上，计算机是尝试所有有效的输入变量组合，以模拟所有可能的结果。

图2的结果显示了模型输出变量对于输入分布的灵敏度，由此可确定模型中最“关键”的输入变量。可以看出下沉系数是影响下沉量最敏感的参数，相关系数值0.85。最不敏感的参数为开采影响传播角，相关系数值为0.004。主要影响角正切和拐点偏移距的相关系数值分别为0.09和0.03。

这个排序是符合开采沉陷规律的，最大下沉值的计算公式W0=mqcosα，下沉系数q直接影响下沉值的大小;主要影响角正切tanβ决定了下沉曲线的开口大小，进而影响下沉值，因为在采煤空间一定的情形下，曲线开口越大，最大下沉值就越小;拐点偏移距s/H0影响着工作面的有效计算尺寸，拐点偏移距越大，有效计算尺寸越小，下沉值就越小;开采影响传播角θ0主要决定了下沉曲线的拐点的分布形态，对下沉值影响极小。

下沉系数的相关系数值0.85的计算方法如下:在模拟过程中，采用拉丁超立方抽样所得的下沉系数和模拟过程中每次迭代计算所得的下沉值之间的关系可用图3所示的散点图表示。该图详细描述了作为输入项的下沉系数和模拟所得到的下沉值之间的关系。由该散点图可以算出下沉系数与下沉值的相关系数为0.85。

图2 各参数敏感指数Fig.2 Sensitivity indexs of parameters

图3 下沉系数－下沉值散点图Fig.3 q－ W scatter diagram

遍历下沉系数的特定输入时，如表2所示，运行每一次模拟都会给出下沉值的平均值和标准差等。输出的下沉值与每次模拟输入所使用的下沉系数值绘图，即可生成图4(a)所示的线形图，该图描述了输入的每一个所选的(下沉系数)步长值运行模拟时，被跟踪模拟的下沉值是如何变化的。图4中的(b)、(c)、(d)分别描述了开采影响传播角、主要影响角正切和拐点偏移距s/H0与下沉值的关系。

表2 下沉系数的特定输入Tab.2 Particular inputs of subsidence coefficient

图4 下沉值与各输入项的线形图Fig.4 Alignment of subsidence values and inputs

5 结论

利用多元逐步回归的方法对概率积分法模型进行参数的全局敏感性分析，结果表明:下沉系数q是影响模型的最关键因素，依次为主要影响角正切tanβ，拐点偏移距s/H0，开采影响传播角θ0。因此，在对地表进行开采沉陷预计时，应尽可能准确的选择q和tanβ，这一结论与工程实践相一致。本文从数学分析的层面为此提供了依据。

［1］何国清，杨伦，凌赓娣.矿山开采沉陷学［M］.徐州:中国矿业大学出版社，1991.

［2］周大伟，吴侃，刁新鹏，等.测量误差对概率积分法参数精度影响［J］.辽宁工程技术大学学报:自然科学版，2013，32(3):316－320.

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

［4］赵春虎.高承压含水层水文地质参数局部和全局灵敏度分析［J］.煤炭科学技术，2013，41(8):110 －113.

［5］姜志伟，陈仲新，周清波，等.CERES－Wheat作物模型参数全局敏感性分析［J］.农业工程学报，2011，27(1):236 －242.

［6］孙艳伟，把多铎，王文川，等.SWMM模型径流参数全局灵敏度分析［J］.农业机械学报，2012，43(7):42 －48.

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

［8］刘兴坡.基于径流系数的城市降雨径流模型参数校准方法［J］.给水排水，2009，35(11):213 －217.

［9］吴志伟，宋汉周.基于全局灵敏度分析的大坝温度场影响因子探讨［J］.水利学报，2011，42(6):737 －742.

［10］徐崇刚，胡远满，常禹，等.生态模型的灵敏度分析［J］.应用生态学报，2004，15(6):1 056 －1 062.