敏感性分析是模型验证, 模型优化和模型诊断的一种很有效的工具, 其广泛应用于经济、生态、工程、化学等领域大型非线性模型中.敏感性分析(Sensitivity Analysis), 就是评价不确定的模型输入参数对模型输出不确定性的贡献率^[1-2].敏感性分析方法可以分为局部敏感性分析和全局敏感性分析两大类, 局部敏感性分析只单独分析单个参数对模型输出的影响, 全局敏感性分析则可以同时分析多个参数对模型输出的总影响, 从而不会忽略多个参数共同作用对输出的影响.常见的全局敏感性分析方法包括Morris法^[3]、FAST法^[4]、GLUE法^[5]、Sobol法^[6]、Extend FAST法^[7]等.其中Sobol方法基于方差分解的原理, 可用于非线性、非单调的数学模型、结果稳健可靠、并且能够给出参数的各阶敏感性系数, 应用广泛.

Sobol方法的基础是蒙特卡洛方法, 通过蒙特卡洛模拟随机抽样计算敏感性分析系数.蒙特卡洛方法^[8]是以概率统计理论为基础的一种数值计算方法, 它的关键是随机数的生成.随机数的产生可以分为三大类:查表方法, 物理方法, 数学方法.在计算机上使用数学方法产生随机数是目前发展最快, 使用最广的方式, 已有位移法, 线性同余法, 反馈位移寄存法等等.伪随机数^[9]便是由上述类似的确定的算法生成的.目前的伪随机数已经被证明有很多缺陷, 例如样本点之间的"间隙"减少了它在分布上的均匀性.因此, 一些其他方法被提出来以克服这些缺点, 如准随机序列, 拉丁超立方体抽样, 它们相较于伪随机数具有更好的均匀性.均匀性会影响蒙特卡洛方法的收敛速度与精确度^[10], 使用伪随机序列蒙特卡洛方法的收敛速度为$O(N^{-1/2})$, 而使用准随机序列的蒙特卡洛方法的收敛速度为$O((\log N)^kN^{-1})$, 在文献[10-11]中均有证明.因此能产生更均匀的随机数的抽样器, 对于加速蒙特卡洛方法计算, 并得到更为精确的结果, 从而影响Sobol方法所得敏感性系数的精度至关重要.因此我们对比上述三种随机抽样器所产生随机数的性质(主要是均匀性)以及其所得敏感性系数, 以得到更好的抽样器代替伪随机抽样计算敏感性系数, 便于我们更好的应用敏感性系数进行一系列研究.

本文首先介绍准随机序列和拉丁超立方体抽样的基本算法, 以及基于蒙特卡洛方法的敏感性分析Sobol方法; 接着比较伪随机抽样器与准随机和拉丁超立方体抽样器的抽样结果, 收敛性与均匀性的比较; 然后比较它们对于最终敏感性系数计算结果的影响.发现准随机抽样相对于另外两种抽样方法具有明显优势.

准随机数的生成有很多方法, 如Van Der Corput (VDC)生成器, Hault's生成器, Sobol's生成器等等.它们遵循类似的策略, 使得产生的数尽可能的避开其他数, 以防止样本点聚集和留下"孔洞".本文中使用基于VDC序列的Halton序列生成准随机数.

Halton序列是一维Van Der Corput (VDC)序列一般化到多维序列^[9].每一维Halton序列是一个VDC序列, 可以表示为

$ x(k), ~~k=1, 2, 3, \cdots . $

(2.1)

接着VDC使用反函数方法来构造准随机数$x(k)$.分成三步, 如下

1) 指定一个素数作为VDC的基础, $b=2, 3, 5, 7, \cdots$;

2) 对每个指标$k$, 基于预先定义的宽度$w$给出按位表示:$k={(d_w d_{w-1}\cdots d_1 d_0)}_b$;

3) 使用反函数: inverse$(k)=\frac{d_0}{b}+\frac{d_1}{b^2}+\cdots +\frac{d_w}{b^{w+1}}$, 这就是$x(k)$的值, 在形式上$x(k)={(0.d_0 d_1\cdots d_w)}_b$.

然而, 当Halton序列扩展到高维时, 样本点在高维子空间的投影失去了它们的均匀分布的特性.所以在我们的实现过程中, 我们更改了VDC的步骤3, 使用折叠基的反函数^[12]

$ {\rm{inverse}}(k) = \frac{{({d_0} + 0)\bmod b}}{b} + \frac{{({d_1} + 1)\bmod b}}{{{b^2}}} + \cdots + \frac{{({d_w} + w)\bmod b}}{{{b^{w + 1}}}}, $

(2.2)

使用这个方法, 我们可以将Halton序列扩展到最大150维.

我们模仿伪随机抽样器得到准随机抽样器.所得准随机抽样相较于伪随机抽样的不同点是它有两种抽样方式.第一种顺序法采用样本的大小作为它的自变量, 产生一个样本(一个向量参数).第二种方法平行法是专门针对准随机抽样的, 因为准随机数之间有很强的相关性, 因此连续产生的样本之间可能有从属性, 这将降低蒙特卡洛模拟的性能.例如, 在Sobol敏感性分析方法中, 它每次需要两个样本(抽样和重抽样).因此, 第二种方法以一种平行的方式产生样本.例:我们要对$m$个参数抽样, 样本数量为$n$.我们若需要$k$个这样的样本, 在顺序法中, 我们需要抽样$k$次, 但是在平行的方式中, 我们只需要1次平行抽样方法就能同时得到$k$个样本.事实上, 平行抽样的准随机生成器将使用一个$m\times{k}$维的序列代替在顺序法中的一个$m$维序列.尽管顺序抽样也包含随机成分尽可能保证样本的独立性, 但是实验结果显示平行抽样仍然比顺序抽样更好.平行抽样的缺点是它很容易超过准随机数生成器支持的最大维度(Halton是150).

我们使用上述方式得到所需的准随机抽样器产生所需准随机数.

我们根据拉丁超立方体抽样定义^[13]的经典算法实现拉丁超立方体抽样器.首先创造一个拉丁超立方体, 并均匀分段.接着, 通过在每一维中两个样本点不会有相同的坐标的方式, 在超立方体中选择样本点.因此它保证了样本点均匀地展开在"立方体"中.然后在每个选择的分段区间挑选一个数(这可以被看作是生成一个均匀分布的伪随机数, 它的均匀分布界限设置为分段区间).

拉丁超立方体抽样的缺点是耗时, 因为它的时间复杂度为$O(N^2)$, 其中$N$是样本数.在一些应用的代码中, 这个方法也消耗内存, 因为它将产生临时矩阵存储来计算结果.我们在实现中使用STL向量容器来客服它, 这将不会产生任何临时矩阵, 有效地加快算法.

Sobol敏感性分析方法是基于方差分解的^[6].假设$f(x)$是定义在$I^n$上的可积函数.于是, 它能够表示为如下形式

$ f(x)=f_0+\sum\limits_{s=1}^{n}\sum\limits_{i_1 <\cdots <i_s}^{n}f_{i_1\cdots {i_s}}(x_{i_1}\cdots {x_{i_s}}), $

(2.3)

该表示具有如下特性

$ \int_{0}^{1}f_{i_1\cdots {i_s}}(x_{i_1}\cdots {x_{i_s}}){\rm d}x_k=0, \;\;\;k=i_1, \cdots, i_s . $

(2.4)

接着就能推断出

$ \int_{0}^{1}f(x){\rm d}x=f_0, {\rm d}x=d_1d_1\cdots {d_n}. $

(2.5)

$ \int_{0}^{1}f^2(x){\rm d}x=f_0^2+\sum\limits_{s=1}^{n}\sum\limits_{i_1<\cdots <i_s}^{n}f_{i_1\cdots {i_s}}^2(x_{i_1}\cdots {x_{i_s}}) . $

(2.6)

从(2.3)式, 能够推断出$f(x)$的方差分解.

$ V_{i_1\cdots {i_s}}=\int_{0}^{1}f_{i_1\cdots {i_s}}^2{\rm d}x_{i_1}\cdots {{\rm d}x_{i_s}} . $

(2.7)

$ V=\int_{0}^{1}f^2(x){\rm d}x-f_0^2=\sum\limits_{s=1}^{n}\sum\limits_{i_1<\cdots <i_s}^{n}V_{i_1\cdots {i_s}} . $

(2.8)

在Sobol的策略中, 输入参数的敏感性的评价来自它对输出参数总方差的贡献量.因此, 我们有

$ S_{i_1\cdots {i_s}}=\frac{V_{i_1\cdots {i_s}}}{V}. $

(2.9)

常用的敏感性系数有:一阶敏感性系数和总一阶敏感性系数.其中一阶敏感性系数$(S_i)$由下式计算

$ S_i=\frac{V_i}{V}. $

(2.10)

这只包括某一输入参数在输出方差上的影响, 但是忽略了它与其他输入参数间的相关性(二阶, 或更高阶).

总一阶敏感性系数$(S_{T_i})$定义如下

$ S_{T_i}=\frac{V-V_{\sim{i}}}{V}, $

(2.11)

其中$V_{\sim{i}}$表示除了$i$以外的所有输入参数带来的方差.因此$(S_{T_i})$意味着来自输入参数$i$的方差的总贡献, 包括它的一阶和与其他输入参数相关的更高阶方差.

而为了计算Sobol敏感性系数, 使用了蒙特卡洛方法.假设$y$是所有输入参数的子集$K$中的一个元素, $z$是$K$外面的另一个参数, 则有

$ V_y=\int{f(x)f(y, z'){\rm d}x{\rm d}z'}-f_0^2. $

(2.12)

有如下证明

$ \int{f(x)f(y, z') {\rm d}x{\rm d}z'}=\int{\rm d}y \int{f(y, z){\rm d}z}\int{f(y, z'){\rm d}z'} =\int{{\rm d}y\Big[\int{f(y, z){\rm d}z}\Big]^2}, $

(2.13)

$ \int{f(y, z){\rm d}z}=f_0+\sum\limits_{s=1}^{m}\sum\limits_{(i_1 <\cdots <i_s)\in{K}}f_{i_1\cdots {i_s}}(x_{i_1}, \cdots, x_{i_s}), $

(2.14)

$ \int{f(x)f(y, z'){\rm d}x{\rm d}z'}=f_0^2+\sum\limits_{s=1}^{m}\sum\limits_{(i_1<\cdots <i_s)\in{K}}V_{i_1\cdots {i_s}}=f_0^2+V_y. $

(2.15)

对于$(S_{T_i})$的计算, 有

$ V_z=\int{f(x)f(y', z){\rm d}x{\rm d}y'}-f_0^2 $

(2.16)

和

$ S_{T_i}=1-S_z. $

(2.17)

因此只需要利用蒙特卡洛方法计算如下4个值

$ \int{f(x){\rm d}x}, $

(2.18)

$ \int{f^2(x){\rm d}x}, $

(2.19)

$ \int{f(x)f(y, z'){\rm d}x{\rm d}z'}, $

(2.20)

$ \int{f(x)f(y', z){\rm d}x{\rm d}y'}, $

(2.21)

这里$x$是抽样点, $(y, z, y', z')$是重抽样点, 它们需要相互独立.如果$y$或$z$只有一个元素, 我们仅仅用一个抽样矩阵$A$和一个重抽样矩阵$B$, 以及$A_B^i$来表示除了第$i$列来自$B$, 所有列都来自$A$的矩阵(矩阵$B_A^i$, 与$A_B^i$对应).那么上面的式子可以写为

$ \frac{1}{N}\sum\limits_{j=1}^{N}f(A)_j, $

(2.22)

$ \frac{1}{N}\sum\limits_{j=1}^{N}f^2(A)_j, $

(2.23)

$ \frac{1}{N}\sum\limits_{j=1}^{N}f(A)_jf(B_A^i)_j, $

(2.24)

$ \frac{1}{N}\sum\limits_{j=1}^{N}f(A)_jf(A_B^i)_j . $

(2.25)

在文献[14]中, 作者提出了一个方法用来改进Sobol经典方法中的一阶和总一阶系数的计算.在该版本中, 在计算上述4个值时对抽样和重抽样点求均值.在我们的实验中, 采用该方法来计算Sobol系数.

$ \frac{1}{2N}\sum\limits_{j=1}^{N}(f(A)_j+f(B)_j), $

(2.26)

$ \frac{1}{2N}\sum\limits_{j=1}^{N}(f^2(A)_j+f^2(B)_j), $

(2.27)

$ \frac{1}{2N}\sum\limits_{j=1}^{N}(f(A)_jf(B_A^i)_j+f(B)_jf(A_B^i)_j), $

(2.28)

$ \frac{1}{2N}\sum\limits_{j=1}^{N}(f(A)_jf(A_B^i)_j+f(B)_jf(B_A^i)_j) . $

(2.29)

在解决很多缺少或没有有效的解析方法的复杂问题时, Monte Carlo(MC)方法是一个有效方法. MC的第一步就是得到输入参数的大量样本.假设每个输入参数对它的可能值有一个范围, 那么每个样本(输入参数的一组可能值)能够由多维向量空间中的一个点来表示.一个好的抽样器产生的样本点应该均匀分布在抽样空间.最普通的抽样方法就是用伪随机数.然而, 伪随机抽样已经被证明有很多缺陷, 例如样本点之间的"间隙"减少了它在分布上的均匀性.所以, 其他一些抽样方法被提出以克服这些缺点.两个候选方法便是第2节中提到的:准随机序列和拉丁超立方体抽样.

准随机序列(也称为低偏差序列)与伪随机数的不同点体现在两方面: 1)即使样本的大小很小, 它也能产生均匀分布的样本点; 2)从统计学的角度看, 序列中的数字并不是相互独立的.由于缺少了"随机"性, 准随机序列并不能替代伪随机数, 但是在蒙特卡洛模拟中, 它在分布上的超均匀性将改善MC方法的精确度.

拉丁超立方体抽样采用另一种方法改进伪随机抽样.它使用的策略是分层抽样, 就是将每一个变量值的区间分为$N$个间隔($N$经常等于抽样的大在$M\times N$超立方体中($M$是输入参数的维度), 每一个样本点有不同的坐标.在每一个选择的间隔中, 生成伪随机数.因为每个可能的间隔都被用在样本中, 它在分布上比纯伪随机数有更好的均匀性.

在图 1中展示了伪随机抽样器的二维抽样, 可以很清楚的看到样本不是完全均匀分布的, 样本点之间有许多"间隙"或"孔洞".当样本数量增加时, 这些"间隙"将慢慢被填满, 但仍会留下一些很小的"孔洞"见图 2.这个缺陷将会降低蒙特卡洛模拟的性能, 这也是为什么提出准随机抽样和拉丁超立方体抽样来克服它.

接下来我们比较准随机抽样, 拉丁超立方体抽样和伪随机抽样的结果.为了比较它们的特性, 我们首先使用均匀分布的数作为测试目标.假设随机数分布在0, 1之间, 则样本的理论均值为0.5, 标准偏差是0.288675.在我们的测试中, 我们计算理论值与实际值之差的绝对值, 分析它们的收敛速度.

从图 3, 图 4中可以看出准随机抽样比伪随机抽样的收敛速度更快.对于标准偏差的收敛, 准随机抽样在样本数量100左右有一个转折点.在样本数量小于100之前, 标准偏差的收敛速度是减小的, 当样本数量足够大时它开始变得正常.这个现象可能是由构造准随机数的机制导致的.在刚开始, 数字猛烈地"跳跃"来充满整个大空间, 这使得它的基于均值的标准偏差增加.但是随后由于大的"间隙"逐渐被充满, 它只需要"跳跃"一个小区域便可以了.

同时, 我们也可以发现拉丁超立方体抽样器相对于伪随机抽样器有重大的改进.当样本数量超过100时, 它和准随机抽样器有一样好的特性.

另一种评价样本特性的方式是样本分布的均匀性.我们分别产生100, 200, 400, 600, 800, 1000的伪随机, 准随机和拉丁超立方体随机样本, 同时采用卡方检验来确定样本是否服从均匀分布, 规定: 0假设为样本服从均匀分布, 1假设为样本不服从均匀分布.采用$P$值$(P\in{[0, 1]})$衡量, $P$值越趋近于0, 表示越有理由拒绝0假设, 即样本不服从均匀分布; $P$值越趋近于1, 表示越有理由接受0假设, 即样本服从均匀分布.我们进行多次实验(1000次), 其中一次结果如表 1所示.我们发现准随机抽样和拉丁超立方体抽样的$P$值基本都接近1, 而伪随机抽样的$P$值随样本数量没有显著趋势, 比较随机且$P$值离1都较远, 在多次实验中发现伪随机抽样$P$值不稳定, 且都远离1.可以看出准随机抽样和拉丁超立方体抽样比伪随机抽样有更好的均匀性.

目前已有很多不同的敏感性分析方法.其中有些研究给定解决方案的敏感性, 称作"局部敏感性"分析; 另一些方法分析整个模型的敏感性, 称作"全局敏感性"分析.在我们的文章中, 我们使用由Sobol发展的全局敏感性计算方法.该算法是基于蒙特卡洛模拟的(参见表 2和表 3).

在实验中, 我们需要一些基准模型进行评价.我们选择两个敏感性分析经典解析模型: Ishigami函数和$G$函数, 它们的一阶和总一阶系数都有理论值.

Ishigami函数是由Ishigami和Homma(1990)提出的, 参见文献[15].它有3个输入参数

$ f(X_1, X_2, X_3)=\sin X_1+a\sin^2X_2+bX_3^4\sin X_1, $

(4.1)

它的总方差$V$和局部方差$V_i$如下所示

$ V=\frac{a^2}{8}+\frac{b\pi^4}{5}+\frac{b^2\pi^8}{18}+\frac{1}{2}, $

(4.2)

$ V_1=\frac{b\pi^4}{5}+\frac{b^2\pi^8}{50}+\frac{1}{2}, $

(4.3)

$ V_2=\frac{a^2}{8}, $

(4.4)

$ V_3=0 . $

(4.5)

$G$函数是另一个由Sobol提出的解析模型, 它是一个非单调和非可加函数^[16].它有$k$个输入参数$x_1, x_2, \cdots, x_k$.每个输入参数都分布在$[0, 1]$区间.

$ g(x_1, x_2, \cdots, x_k)=\prod\limits_{i=1}^{k}g_i(x_i), $

(4.6)

其中

$ g_i(x_i)=\frac{\mid4x_i-2\mid+a_i}{1+a_i}. $

(4.7)

因此我们有$g_i(x_i)$的方差和总方差

$ V_i=\frac{1}{3(1+a_i)^2}, $

(4.8)

$ V=\prod\limits_{i=1}^{k}(1+V_i)-1, $

(4.9)

$ V_{T_i}=V_i\prod\limits_{j\neq{i}}(1+V_j), $

(4.10)

在我们接下来的实验中, 我们取$k=6, a_i=\{0, 0.5, 3, 9, 99, 99\}$.

正如前面所讨论的, Sobol方法计算敏感性系数是基于蒙塔卡洛模拟的.而随机抽样是蒙特卡洛模拟的第一步, 因此, 抽样的质量将影响计算的敏感性系数的精确度.这里, 我们研究伪随机, 准随机和拉丁超立方体随机抽样的性能.为了增加对比的客观性, 我们分别对Ishigami和$G$函数模型进行独立实验, 而且我们同时评价一阶系数和总一阶系数.

在图 5, 图 6中, 我们可以看到准随机抽样的$S_i$比其他两种抽样方式有显著优势.一方面, 准随机抽样的$S_i$相较于另两种方式有更精确的均值(表 2, 表 3).另一方面, 它的置信区间比伪随机抽样和拉丁超立方体抽样都要窄很多, 其结果精确度更高.

总一阶系数$S_{T_i}$是另一个经常使用的系数.在图 7和图 8中, 我们仍旧只能够观察到准随机抽样的$S_{T_i}$胜过其他的抽样. 表 4和表 5也显示了准随机抽样比其他两中抽样方式更加精确.

上面实验的抽样量都是1000, 拉丁超立方体不能展示出它相对于伪随机抽样的优势.这个现象是很奇怪的, 因为拉丁超立方体抽样比伪随机抽样有更好的收敛性和均匀性.然而, 如果我们将样本数从1000降低到100, 拉丁超立方体抽样就展示出相对于伪随机抽样显著的改进.在Ishigami模型中分别计算$S_i$, $S_{T_i}$, 有如表 6和表 7中的结果, 可以清楚的看到在低样本量时, 拉丁超立方体抽样相对于伪随机抽样的的结果更为精确.

我们的主要工作是实现不同的随机抽样器.并通过一系列对Sobol敏感性系数计算的实验, 我们发现Sobol方法基于准随机抽样比伪随机抽样和拉丁超立方体抽样具有显著优势.这个优势来自于准随机抽样的样本点有更好的分布均匀性的事实, 这改善了蒙特卡洛模拟的性能.对于Sobol方法, 它展示在两方面

1) 敏感性系数更精确;

2) 敏感性系数的置信区间更窄.

但是我们也观察到一些现象需要我们在未来的工作中进一步研究.

在实验中, 我们发现准随机抽样的顺序法比平行法的性能差.这可能是因为准随机抽样随机性的缺乏, 这导致连续产生的样本之间有轻微的依赖性.应该提出一个随机性更好的方法克服它.因为顺序法对使用者来说是更普遍的, 且它能减少维度的要求.

拉丁超立方体抽样的性能并不如预期的那么好.在小样本量时, 其结果比伪随机更精确, 而在大样本量时, 它相较于伪随机抽样并没有改进.因此, 应该探究它的缺点的原因, 并提出一种更好的实现方法.