0 引言

敏感度分析(Sensitivity Analysis)用于定性或定量地评价模型参数不确定性对模型结果产生的影响，包括模型输出和模型自身性能的影响.敏感度分析具有诊断和预测能力，在工程应用领域中被看作是科学建模和模型分析的先决条件.敏感度分析可以分为局部敏感度分析和全局敏感度分析^[1].局部敏感度分析只考察各个模型不确定输入在局部范围内的小变化对模型结果的影响程度，只有当模型输入和模型输出是线性或近似线性依赖关系时，局部敏感度分析才能较好地把握模型参数不确定性对整体模型性能的影响.全局敏感度分析则研究的是各个模型不确定输入在其允许的变化范围内对模型结果的全局影响，它并不要求模型是线性系统，因此全局敏感度分析还包括考察模型参数之间相互作用对模型结果的影响^[2].

全局敏感度分析方法主要有多元回归法^[3-5]、Morri法^[6]、傅里叶幅度敏感度检验法(FAST)^[7]以及基于方差分解的Sobol法^[8]等.Sobol敏感度分析方法是一种基于方差分解的蒙特卡罗法.基于方差分解的Sobol方法通过衡量随机输入参数对其输出响应量的方差贡献程度来表征各个参数对模型系统的影响程度，其优点是考虑了随机输入参数在整个取值范围内对输出响应的贡献以及随机输入参数的交互作用，并结合蒙特卡罗法定量地计算Sobol全局敏感度指标.Sobol敏感度分析方法作为一种蒙特卡罗法，根据随机输入概率分布采用随机数产生器生成大量的输入样本，对每一个输入样本进行一次数值模拟得到相应的输出响应样本，然后对这些大量的响应样本进行统计推断得到Sobol全局敏感度指标.基于蒙特卡罗的Sobol方法所需的样本数量非常巨大，对于实际的复杂工程问题，数值模拟的每一次计算开销往往非常大，很难实现蒙特卡罗方法所需庞大样本的计算.

为了避免对实际问题的庞大样本模拟，可以先构造高阶或高分辨率响应面模型(如Kriging响应面模型^[9-10])来替代真实物理模型，然后对响应面模型采用基于蒙特卡罗的Sobol方法计算Sobol全局敏感度指标，从而实现基于方差分解的全局敏感度分析.一般地，这种响应面替代方法会涉及两种数值误差，一种是替代响应面模型与真实物理模型的误差，另一种是蒙特卡罗误差(样本数无穷大时，误差才会收敛到零).两种误差的叠加将会对Sobol全局灵敏度指标计算结果的置信度产生重大影响.相反，如果可以直接避免蒙特卡罗误差，将极大地改善全局敏感度指标计算的精度.近年来，多项式混沌方法在不确定量化领域取得了巨大成功^[11-12]，主要原因是多项式混沌方法直接避免了蒙特卡罗误差.如果采用多项式混沌展开作为高阶替代模型，则由展开式系数不仅可以直接计算输出响应的均值和方差等概率统计量，也可以直接计算Sobol全局敏感度指标.

实际上，全局敏感度分析与不确定量化研究密切相关，Sudret^[13]最先将多项式混沌方法推广运用到基于方差分解的全局敏感度分析中，替代传统的蒙特卡罗法(即Sobol方法)，实现对Sobol全局敏感度指标的计算.多项式混沌方法的一大重要优势在于：如果响应面函数是光滑的，所需计算展开式系数的数量较少，具有谱方法的收敛速度，相比传统蒙特卡罗方法所需计算时间较少.但是，多项式混沌方法需要求解展开式系数，对于高维随机空间或多项式混沌高阶展开，多项式混沌方法将会面临严重的“维数灾难”问题(curse of dimension).为了克服这一困难，Blatman和Sudret^[14-15]提出了稀疏多项式混沌方法用于基于方差分解的敏感度分析中，他们采用自适应最小二乘回归技术^[14]或最小角回归技术^[15]来保留多项式混沌展开中的主要贡献项，由此避免求解所有多项式混沌展开项系数.稀疏多项式混沌展开方法还可以采用近年来流行的压缩感知技术(compressive sensing)来获得重要的稀疏多项式混沌展开系数^[16-17]，在样本数远小于多项式混沌展开项数的情况下，压缩感知技术将l₀范数极小化问题松弛为l₁范数极小化问题(简称为l₁优化问题)，从而可以采用各种高效的l₁优化方法来计算稀疏多项式混沌展开系数.

1 多项式混沌方法

多项式混沌方法最早起源于Wiener各向同性混沌理论中的随机变量谱展开^[18]，用于湍流问题研究.上世纪九十年代，Ghanem和Spanos将这种随机变量的谱展开用来刻画固体力学有限元分析中的不确定因素^[19]，成功应用于工程结构力学中的各种问题.本世纪初，Xiu和Karniadakis将早期基于高斯随机变量的Hermite多项式混沌推广到Askey多项式混沌^[20-21]，适用于更一般的随机变量概率分布(均匀分布、Gamma指数分布等)，这种广义多项式混沌方法被他们首先应用于计算流体力学(CFD)不确定量化分析中.而Sudret^[13]将多项式混沌方法应用于基于方差分解的全局敏感度分析中，可以便捷地实现对Sobol全局敏感度指标的计算.下面就多项式混沌方法进行基本的介绍.

考虑一个概率空间(Ω,F,P)，其中Ω为样本空间,F为定义其上的一个σ代数,P是定义在可测空间(Ω,F)的概率测度.考虑一组定义在概率空间(Ω,F,P)上的N维随机变量x=(x₁,x₂,…,x_N)^T，每个边缘概率密度函数为f_i(x_i).假设这组N维随机变量是相互独立的，则联合概率密度函数为f(x)=f₁(x₁)f₂(x₂)…f_N(x_N).二阶随机响应面函数可以展开成广义多项式混沌的级数形式^[20]

$\begin{align} & Y\left( x \right)={{c}_{0}}{{\psi }_{0}}+\sum\limits_{{{i}_{1}}=1}^{N}{{{c}_{{{i}_{1}}}}{{\psi }_{1}}\left( {{\xi }_{{{i}_{1}}}}\left( x \right) \right)}+ \\ & \sum\limits_{{{i}_{1}}=1}^{N}{\sum\limits_{{{i}_{2}}=1}^{N}{{{c}_{{{i}_{1}},{{i}_{2}}}}{{\psi }_{2}}{{\xi }_{{{i}_{1}}}}\left( \left( x \right),{{\xi }_{{{i}_{2}}}}\left( x \right) \right)}}+\cdots \\ \end{align}$

(1)

可以写成紧凑形式

$Y\left( x \right)=\sum\limits_{i=0}^{\infty }{{{c}_{i}}{{\Phi }_{i}}\left( \xi \right)}.$

(2)

广义多项式混沌由单变量正交多项式的张量积构成，其满足正交关系

$E[{{\Phi }_{i}}{{\Phi }_{j}}]={{\delta }_{ij}}E[{{\Phi }_{i}}^{2}],$

(3)

这里，δ_ij为Kronecker算符，E[·]为数学期望.表 1给出了不同概率分布(连续或离散)下对应的多项式混沌展开基函数，比如高斯分布对于Hermite正交多项式，均匀分布对应Legendre正交多项式.对于任意概率分布的情形，还可以采用Stieltjes方法或改进Chebyshev方法^[22]来构造基于任意概率测度(连续或离散测度)的多项式混沌基函数.

多项式混沌展开式的系数可由正交关系确定为

${{c}_{i}}=E[Y\left( x \right){{\Phi }_{i}}\left] /E \right[{{\Phi }_{i}}^{2}].$

(4)

如果采用不高于p阶的多项式混沌进行截断近似，则有

$\tilde{Y}\left( x \right)=\sum\limits_{i=0}^{P-1}{{{c}_{i}}{{\Phi }_{i}}\left( \xi \right)},$

(5)

截断近似保留的多项式混沌有P项，可由下式确定

$P=\left( \begin{matrix} N+p \\ p \\ \end{matrix} \right)=\frac{\left( N+p \right)!}{N!p!}.$

(6)

多项式混沌展开最大的好处就是可以由展开式系数直接计算输出响应的均值和方差等概率统计量，形式如下

$E\left[ \tilde{Y}\left( x \right) \right]={{c}_{0}},$

(7)

$Var\left[ \tilde{Y}\left( x \right) \right]=\sum\limits_{i=1}^{P-1}{{{c}_{i}}^{2}E[{{\Phi }_{i}}^{2}]}.$

(8)

多项式混沌方法可以分成嵌入式(intrusive)和非嵌入式(non-intrusive)两种，嵌入式方法主要基于Galerkin投影的思想，将多项式混沌展开式代入物理控制方程，利用多项式混沌的正交性特征，得到以多项式混沌展开系数为待求函数的控制方程，然后采用数值方法计算求解新的耦合控制方程，显然针对新控制方程需要重新编写求解器.而非嵌入式方法可以采用已有的求解器，基于配置法的思想，在随机变量空间选取配置点，将配置点参数作为求解器的输入参数计算目标函数值，然后计算目标函数的多项式混沌展开系数.本文主要考虑非嵌入式多项式混沌方法，相关的计算方法可以分成两大类：投影方法和回归方法.投影方法主要包括蒙特卡罗积分法和高斯积分法.

1.1 蒙特卡罗积分法

蒙特卡罗积分法主要取决于随机变量的抽样，这包括蒙特卡罗抽样(MC)^[23]、拟蒙特卡罗抽样(QMC)^[24]、拉丁方抽样(LH)^[25]. 蒙特卡罗积分法不依赖于随机变量的维数，但是它要获得足够的精度需要成千上万的样本数，计算开销巨大.

1.2 高斯积分法

高斯数值积分是一种有效的高精度数值计算方法，可以针对不同概率密度函数进行高斯积分计算展开式系数(4) ，形式如下

$\begin{align} & E[Y\left( x \right){{\Phi }_{i}}]=\sum\limits_{{{i}_{1}}=1}^{K}{\cdots }\sum\limits_{{{i}_{N}}=1}^{K}{{{\omega }_{{{i}_{1}}}}},\cdots , \\ & {{\omega }_{{{i}_{N}}}}Y(x({{\xi }_{{{i}_{1}}}},\cdots ,{{\xi }_{{{i}_{N}}}}))\times {{\Phi }_{i}}(x({{\xi }_{{{i}_{1}}}},\cdots ,{{\xi }_{{{i}_{N}}}})), \\ \end{align}$

(9)

这里，ξ_ij(i_j=1,…,K)在每个随机变量上高斯积分点，而ω_ij(i_j=1,…,K)为相应的权重因子.张量积形式的高斯积分法所需配置点个数为K^N(K~p)，这将造成“维数灾难”问题，随机变量的维数较大时会超过蒙特卡罗方法所需的样本数.为了减轻“维数灾难”问题，可以采用Smolyak稀疏网格积分点^[26]，使得所需配置点数显著减少，而积分精度不受显著影响.Smolyak稀疏网格方法将在附录中作简要介绍.

1.3 回归方法

采用普通最小二乘回归方法(OLS,Ordinary Least Square)来估计响应面系数(即回归系数)，这时回归所需样本点可以采用蒙特卡罗随机抽样获得，响应面系数近似估计由下式得到

$c={{({{\Phi }^{T}}\Phi )}^{-1}}{{\Phi }^{T}}Y.$

(10)

这里，Φ=Φ_j(ξⁱ),i=1,…,M; j=0,1,…,P-1，Φ^TΦ为信息矩阵,c为响应面系数构成的列向量,Y为随机抽样点对应真实的响应面值所构成的列向量.信息矩阵求逆可以采用奇异值分解方法，而回归方法所需抽样个数必须满足M>P以保证信息矩阵可逆.Sudret^[13]给出了关于抽样个数的一个经验最优公式为

$M=\left( N-1 \right)P.$

(11)

2 稀疏多项式混沌方法

对于随机响应面函数的多项式混沌展开而言，随机空间维数和展开阶数的增加，多项式混沌展开式项数将成几何级数急剧增大，出现“维数灾难”现象.表 2给出了在不同随机空间维数和展开阶数下的多项式混沌展开项数，从表中可以清楚地看出这一现象.因此，随机输入变量较多时不能进行高阶的多项式混沌展开，进行高阶的多项式混沌展开时只能有较少的随机输入变量个数，这就限制了多项式混沌方法的适应范围.

随机空间维数和展开阶数的增加，对投影高斯积分法而言，还会引起高斯积分样本点数目的急剧增大；而对最小二乘回归方法而言，所需样本数也远大于多项式混沌展开项数(由上节所述).可见，高维高阶问题要求提供的样本数将非常巨大，而实际工程计算中每个响应样本的获得都需要花费不菲的机时开销，根本无法提供.相反，如果可以提供数目巨大的响应样本，直接采用维数无关的蒙特卡罗积分法计算就可以了.在可提供样本数目较小的情况下，不可能计算出所有多项式混沌展开项；然而实际工程计算中从输入到输出的响应面函数并不一定需要所有的多项式混沌展开项，大量的多项式混沌展开项系数值可能非常小，甚至为零，完全可以舍去.因此，可以构造稀疏的多项式混沌展开来近似响应面函数.至此，我们所面临的问题就是如何构造稀疏多项式混沌展开式.一种笨拙的方法就是计算出所有的展开项系数，然后排序筛除掉系数值小的展开项，得到稀疏多项式混沌展开；这种方法显然没有克服或减轻“维数灾难”问题.另一种方法就是直接定位和计算贡献大(系数值大)的展开项系数，而这种直接定位也不需要太大的样本数，这种方法可以直接借鉴当前流行的“压缩感知”稀疏重构技术.

2.1 基于压缩感知的稀疏重构

如果多项式混沌展开是稀疏的，如何求解定位和求解稀疏多项式混沌展开系数，就构成了稀疏重构问题，可以表示成如下的P₀规划问题

${{P}_{0}}:~\underset{c\in {{R}^{P}}}{\mathop{\min }}\,\|c{{\|}_{0}}~s.t.\text{ }Y=\Phi c.$

(12)

这里，‖c‖₀表示c中非零元素的个数，Φ为测量矩阵,Y为样本观测数据(样本响应面值).由此，求解P₀规划问题就可以得到稀疏多项式混沌展开式.但是，P₀规划问题就是l₀优化问题，它是一个NP-hard问题.为此，可以将这个l₀优化问题松弛为可求解的l₁优化问题，即如下的P₁规划问题

${{P}_{1}}:~\underset{c\in {{R}^{P}}}{\mathop{\min }}\,\|c{{\|}_{1}}~s.t.\text{ }Y=\Phi c.$

(13)

l₁优化问题(13) 也称为BP(basis pursuit)问题^[27]，如果考虑测量噪声或截断误差，可以推广为BPDN(basis pursuit denoising)问题

${{P}_{1,\delta }}:\underset{c\in {{R}^{P}}}{\mathop{\min }}\,\|c{{\|}_{1}}~s.t.\text{ }\|\Phi c-Y{{\|}_{2}}＜\varepsilon .$

(14)

Candes,Romberg和Tao等人^[28]建立的RIP理论(压缩感知的奠基性理论)将P₀和P₁规划问题统一了起来：只要测量矩阵满足所谓的RIP性质，并且RIC常数小于某个指定值(RIP条件)时，可以证明P₀和P₁规划问题是等价的.

矩阵RIP性质(Restricted Isometry Property)的定义：对任意‖x‖₀≤s的稀疏向量x，如果存在常数δ_s∈[0,1) 使得

$(1-{{\delta }_{s}})\|x{{\|}^{2}}\le \|\Phi x{{\|}^{2}}\le (1+{{\delta }_{s}})\|x{{\|}^{2}},$

(15)

则矩阵Φ∈R^n×N满足s阶RIP性质，而δ_s为RIC常数.Candes等人^［28]在RIP理论之初给出的RIP条件较为粗糙，很快Candes^[29]将RIP条件改进为δ_2s≤$\sqrt{2}$-1，Foucart和Lai^[30]又改进为δ_2s≤0.453 1，Cai^[31]则进一步改进为δ_2s≤0.472.特别是，Mo和Li^[32]将RIP条件改进为δ_2s≤0.493 1.满足RIP条件的矩阵就可以简称为RIP矩阵.

基于RIP理论可以用来指导设计测量矩阵，不少学者提出了一些满足RIP条件的确定性测量矩阵，如Chirp测量矩阵、Alltop序列形成的测量矩阵、随机卷积形成的测量矩阵.但是，实际应用时往往采用随机矩阵来作为测量矩阵，而满足RIP条件的随机矩阵主要有Gaussian随机矩阵与Bernoulli随机矩阵.所谓Gaussian随机矩阵,即指矩阵中的元素是独立的随机变量且服从期望为0，方差为1/n的Gaussian分布.所谓Bernoulli矩阵，即指矩阵中的元素以相同的概率取1/$\sqrt{n}$或-1/$\sqrt{n}$.Gaussian随机矩阵与Bernoulli随机矩阵是依据概率满足RIP条件的，有如下定理^[33].

定理1 假定n×N的矩阵Φ是一个Gaussian或者Bernoulli随机矩阵，当

$s\le {{c}_{1}}n/log\left( N/s \right),$

(16)

矩阵Φ是一个s阶RIP矩阵的概率不小于

$1-\exp (-{{c}_{2}}n),$

(17)

这里，常数c₁,c₂仅仅依赖于RIC常数δ_s.

针对稀疏多项式混沌展开，需要考虑的是由多项式混沌基函数(高维正交多项式)构成的测量矩阵是否满足RIP性质，由此保证l₁优化方法可以构造稀疏多项式混沌展开.为此，Rauhut^[34]证明了以有界正交系统(Bounded Orthonormal Systems)为基底的随机测量矩阵依据概率满足RIP性质.显然，Legendre正交多项式就是最常见的有界正交系统，对高维随机输入系统，可以给出如下重要定理：

定理2^{[17, 34]} 设{Φ_j}_0≤j≤P-1为N维独立均匀分布随机变量下的Legendre多项式混沌基函数，并且它们的阶数小于或等于p.从N维均匀分布随机空间[-1,1]^N独立地抽取M个随机样本{ξⁱ}_1≤i≤M，构造随机测量矩阵Φ=Φ_j(ξⁱ)；如果当N>p时，

$M\ge C{{3}^{p}}{{\delta }^{-2}}s{{\log }^{3}}\left( s \right)\log \left( P \right),$

(18)

则在不小于1-P^-γlog3(s)的概率下，矩阵Φ/$\sqrt{M}$的RIC常数δ_s满足δ_s≤δ.其中,C和γ为不依赖于N,p,M的正数.

2.2 基于l₁优化的稀疏多项式混沌方法

由于Legendre正交多项式混沌为基底的随机测量矩阵具有RIP性质(见定理2) ，因此均匀分布的多项式混沌展开的稀疏重构在理论上是有保障的，因此本文只考虑均匀分布的稀疏多项式混沌方法.可以将基于l₁优化的稀疏多项式混沌方法分成如下几个步骤

1) 在均匀分布随机空间采用蒙特卡罗抽样(MC)、拟蒙特卡罗抽样(QMC)、拉丁方抽样(LH)选取样本点{ξⁱ}_i=1^M；

2) 计算得到抽样样本观测数据Y=y(ξⁱ),i=1,…,M，将样本点的随机空间坐标代入多项式混沌展开，得到M×P的测量矩阵Φ=Φ_j(ξⁱ),i=1,…,M; j=0,1,…,P-1；

3) 利用l₁优化算法计算BP问题或BPDN问题；如果考虑BPDN问题，可以采用交叉验证技术(cross-validation)对ε进行推算.(后续算例只考虑BP问题)

4) 最后得到稀疏向量c，由此得到稀疏多项式混沌展开用于全局敏感度分析.

求解P₁规划问题(即BP问题)，最简单的方法就是定义新的待求向量c⁺,c^-.

${c_j}^ + = \left\{ \matrix{ {c_j},if\quad {c_j} > 0, \hfill \cr 0,if\quad {c_j} \le 0, \hfill \cr} \right.\quad {c_j}^ - = \left\{ \matrix{ 0,if\quad {c_j} > 0, \hfill \cr - {c_j},if\quad {c_j} \le 0. \hfill \cr} \right.$

(19)

根据c=c⁺-c^-和c=c⁺+c^-，将P₁规划问题转化为线性规划问题，形式如

$\underset{{{c}^{+}},{{c}^{-}}\in {{R}^{P}}}{\mathop{\text{minimize}}}\,\sum\limits_{j=1}^{P}{({{c}_{j}}^{+}+{{c}_{j}}^{-})}s.t.Y=\left[ \Phi |-\Phi \right]\left[ \begin{matrix} {{c}^{+}} \\ {{c}^{-}} \\ \end{matrix} \right].$

(20)

线性规划问题可以采用传统的单纯形算法和内点算法，存在可靠和便于使用的标准软件.但是，这种转化成线性规划问题的处理方式，扩大了计算维数空间，计算效率不高，实际计算可以采用专门针对l₁优化的高效算法(包括改进的内点算法).有效集算法(Active Set)^[34-36]和投影梯度法(Projected Gradient)^[37-39]是普遍采用的两类l₁优化方法.实际上，Blatman和Sudret^[15]提出的基于最小角回归的稀疏多项式混沌方法就属于一种有效集算法.本文计算采用的l₁优化算法是van den Berg和Friedlander^[39]提出的谱投影梯度算法(Spectral Projected Gradient algorithm (SPGL1) )，该方法适用于大尺度稀疏信号的重构计算.

3 全局敏感度分析

进行全局敏感度分析首先需要定义全局敏感度指标，用以表征不同随机输入对模型响应输出的影响程度.目前基于方差分解的全局敏感度指标是最普遍应用的一个指标形式，通过计算单个输入以及多个输入对输出方差的贡献，评估单个输入和多个输入交互效应的敏感度.它最先由Sobol提出^[40]用作敏感度指标，因此也称为Sobol敏感度指标.Sobol敏感度指标可以直接从多项式混沌展开式得到，这为基于方差分解的全局敏感度分析带来了极大的方便^[13].

3.1 方差分解(ANOVA分解)与Sobol敏感度指标

随机输出响应面函数可以唯一地分解成

$\begin{align} & Y\left( x \right):=\left( x \right)={{M}_{0}}+\sum\limits_{i=1}^{N}{{{M}_{i}}({{x}_{i}})}+ \\ & \sum\limits_{1\le i＜j\le N}{{{M}_{ij}}\left( {{x}_{i}},{{x}_{j}} \right)}+\cdots +{{M}_{1,2,\ldots ,N}}\left( x \right), \\ \end{align}$

(21)

各个分解项之间满足正交关系

$\begin{align} & E[{{M}_{{{i}_{1}},\cdots ,{{i}_{s}}}}({{x}_{{{i}_{1}}}},\cdots ,{{x}_{{{i}_{s}}}}){{M}_{{{j}_{1}},\cdots ,{{j}_{t}}}}({{x}_{{{j}_{1}}}},\cdots ,{{x}_{{{j}_{t}}}})]=0, \\ & s,t=1,\cdots ,N,\{{{i}_{1}},\cdots ,{{i}_{s}}\ne {{j}_{1}},\cdots ,{{j}_{t}}\}, \\ \end{align}$

(22)

对(21) 等式两边取方差，得到ANOVA分解

$VarM\left( x \right):=D=\sum\limits_{i=1}^{N}{{{D}_{i}}}+\sum\limits_{1\le i＜j\le N}{{{D}_{ij}}}+\cdots +{{D}_{1,2,\cdots ,N}}.$

(23)

由上式可以看出分解项表征了不同单个输入及其相互作用对输出响应方差的贡献，因此可以定义Sobol敏感度指标为

${{S}_{{{i}_{1}},\cdots ,{{i}_{s}}}}={{D}_{{{i}_{1}},\cdots ,{{i}_{s}}}}/D,$

(24)

其满足

$\sum\limits_{i=1}^{N}{{{S}_{i}}}+\sum\limits_{1\le i＜j\le N}{{{S}_{i,j}}}+\cdots +{{S}_{1,2,\cdots ,N}}=1.$

(25)

其中的S_i被称为主效应敏感度指标，它表征单个输入对响应方差的贡献.进一步地，还可以定义总效应敏感度指标

$S_{i}^{T}={{S}_{i}}+\sum\limits_{j＜i}{{{S}_{j,i}}}+\sum\limits_{1＜k＜i}{{{S}_{j,k,i}}}+\cdots +{{S}_{1,2,\ldots ,N}}.$

(26)

它表征单个输入及其与其它输入相互作用对响应方差的总贡献.

3.2 Sobol敏感度指标的多项式混沌计算方法

Sobol敏感度指标早期是采用蒙特卡罗方法计算的^{[8, 41]}，广义多项式混沌展开法提出后发现它可以直接由展开式系数得到.首先将多项式混沌截断展开式(5) 改写成

$\begin{align} & \tilde{Y}\left( x \right)={{c}_{0}}+\sum\limits_{i=1}^{N}{\sum\limits_{\alpha \in {{I}_{i}}}{{{c}_{\alpha }}{{\psi }_{\alpha }}({{x}_{i}})}}+ \\ & \sum\limits_{1\le {{i}_{1}}＜{{i}_{2}}\le N}{\sum\limits_{\alpha \in {{I}_{i}},{{i}_{2}}}{{{c}_{\alpha }}{{\psi }_{\alpha }}\left( {{x}_{{{i}_{1}}}},{{x}_{{{i}_{2}}}} \right)}}+\cdots + \\ & \sum\limits_{1\le {{i}_{1}}＜\cdots s\le N}{\sum\limits_{\alpha \in {{I}_{{{i}_{1}},\cdots ,{{i}_{s}}}}}{{{c}_{\alpha }}{{\psi }_{\alpha }}\left( {{x}_{{{i}_{1}}}},\cdots ,{{x}_{{{i}_{s}}}} \right)}}+\cdots + \\ & \sum\limits_{\alpha \in {{I}_{1,\cdots ,N}}}{{{c}_{\alpha }}{{\psi }_{\alpha }}({{x}_{1}},\cdots ,{{x}_{N}})}, \\ \end{align}$

(27)

这里，I_i1,…,is={α∈(α₁,α₂,…,α_N): α_k=0 ⇔ k∉(i₁,…,i_s),∀k=1,…,N}.由此得到

$M({{x}_{{{i}_{1}}}},\cdots ,{{x}_{{{i}_{s}}}})=\sum\limits_{\alpha \in {{I}_{{{i}_{1}},\cdots ,{{i}_{s}}}}}{{{c}_{\alpha }}{{\psi }_{\alpha }}({{x}_{{{i}_{1}}}},\cdots ,{{x}_{{{i}_{s}}}})},$

(28)

然后根据多项式混沌的正交性，可以得到

${{D}_{{{s}_{1}},\cdots ,{{i}_{s}}}}=\sum\limits_{\alpha \in {{I}_{{{i}_{1}},\cdots ,{{i}_{s}}}}}{{{c}_{\alpha }}^{2}},$

(29)

由此,根据定义(24) 就可以计算Sobol敏感度指标了.

4 基本算例

多项式混沌法可以方便有效地进行基于方差分解的全局敏感度分析，下面选取几个典型的随机空间维数较多或多项式混沌展开所需阶数较高的标准模型，计算它们的Sobol敏感度指标，并与解析解或蒙特卡罗数值解进行细致对比，以此展示多项式混沌方法在全局敏感度分析中的成功应用.

4.1 Ishigami函数

考虑Ishigami函数^[42]，

$Y\left( x \right)=\sin {{x}_{1}}+a{{\sin }^{2}}{{x}_{2}}+b{{x}_{3}}^{4}\sin {{x}_{1}}.$

(30)

其中，x₁,x₂,x₃是在[-π,π]上均匀分布的三个独立输入随机变量.Ishigami函数是一个非线性、非单调函数，它是文献中校核全局敏感度方法的一个常用标准模型(benchmark model).Ishigami函数的方差和Sobol敏感度指标可以解析地得到

$\begin{align} & D=\frac{{{a}^{2}}}{8}+\frac{b{{\pi }^{4}}}{5}+\frac{{{b}^{2}}{{\pi }^{8}}}{18}+\frac{1}{2}, \\ & {{D}_{1}}=\frac{b{{\pi }^{4}}}{5}+\frac{{{b}^{2}}{{\pi }^{8}}}{50}+\frac{1}{2},{{D}_{2}}=\frac{{{a}^{2}}}{8},{{D}_{3}}=0, \\ & {{D}_{12}}={{D}_{23}}=0,{{D}_{13}}=\frac{8{{b}^{2}}{{\pi }^{8}}}{225},{{D}_{123}}=0. \\ \end{align}$

在全局敏感度方法校核中，通常采用a=7,b=0.1.首先采用基于张量积高斯积分(Gauss-Legendre积分配置点)的多项式混沌方法计算了Ishigami函数的Sobol主效应和总效应敏感度指标，并与精确解进行比较，见表 3.从表中可以发现多项式混沌展开到9阶，Sobol敏感度指标才得到显著收敛；而张量积高斯积分所需配置点数目为1 331.然后再采用基于稀疏网格(由Gauss-Patterson积分配置点构造)积分的多项式混沌方法计算Ishigami函数的Sobol主效应和总效应敏感度指标，见表 4.从表中可以看出多项式混沌展开到9阶，稀疏网格计算结果具有同样的收敛性，而所需配置点数目(等于495) 远小于张量积高斯积分方法.Sudret^[13]已经采用最小二乘回归方法对Ishigami函数计算了不同展开阶数下的Sobol灵敏度指标，在多项式混沌展开为9阶时，只需要291个拟蒙特卡罗抽样点就可以达到相当的精度，虽然精度不及高斯积分法.最后，采用基于l₁优化的稀疏多项式混沌方法计算Sobol灵敏度指标(见表 5)，只需要120个拟蒙特卡罗抽样点就得到了与回归法相当精度的结果，而所需样本点数目不到回归法的一半；将样本点数增大到160时，精度没有进一步地改进；将样本点数减小到80时，精度稍稍变差，但结果仍优于Blatman和Sudret的采用逐步回归技术的稀疏多项式混沌方法^[14](采用的样本数为84) .

4.2 Sobol函数

考虑Sobol函数^[43]，

$Y\left( x \right)=\prod\limits_{i=1}^{q}{\frac{\left| 4{{x}_{i}}-2 \right|+{{a}_{i}}}{1+{{a}_{i}}}}.$

(31)

其中，x_i(i=1,…,q)是在[0,1]上均匀分布的独立输入随机变量，输出响应面函数Y是一个间断函数，参数a_i为非负数.Sobol函数的方差及其Sobol敏感度指标也可以解析地得到

$\begin{align} & D=\prod\limits_{i=1}^{q}{({{D}_{i}}+1)}-1,{{D}_{i}}=\frac{1}{3{{(1+{{a}_{i}})}^{2}}}, \\ & {{S}_{{{i}_{1}},\cdots ,{{i}_{s}}}}=\frac{1}{D}\prod\limits_{j=1}^{s}{{{D}_{{{i}_{j}}}}}. \\ \end{align}$

(32)

显然，可以看出随着参数a_i的增大，相应输入x_i对输出响应的影响重要性减小.

在校核全局敏感度方法时，通常采用q=8以及a=[1,2,5,10,20,50,100,500].为了克服“维数灾难”，Sudret^[13]最早先考察低阶多项式混沌展开(p=2) 来计算精度较差的Sobol敏感度值，从中判断出前4个随机输入变量对输出方差有重要影响(这种过程属于定性敏感度分析)；然后再固定影响不大的后4个随机输入变量，重新采用高阶多项式混沌展开来研究前4个随机输入变量对输出方差的敏感度效应.Blatman和Sudret^[14]则采用基于逐步回归技术的稀疏多项式混沌方法，重新考察了Sobol函数的全局敏感度，在保持较高展开阶的同时，所需样本也大大地减小了.我们首先采用了最小二乘回归法计算了Sobol函数的全局敏感度指标(见表 6)，发现当阶数p=4时敏感度指标的计算值与精确值符合得较好，而此时需要的样本点个数为3 466.在此阶数下如果采用高斯张量积分或稀疏网格积分来计算全尺度的多项式混沌，所需的样本数是非常庞大的.为此，采用基于l₁优化的稀疏多项式混沌方法计算Sobol敏感度指标(见表 7)，抽样方式为拟蒙特卡罗抽样，最高阶为6时样本数为290，最高阶为10时样本数为660.从表中可以看出，精度大大降低了，这在Blatman和Sudret^[14]的相同算例中也可以看出，这与Sobol函数是一个间断函数有关：多项式混沌方法作为全局方法，不能很好地处理局部间断问题.尽管如此，计算表明对稀疏多项式混沌展开来说，Sobol敏感度排序是有保证的.

4.3 Corner peak函数

考虑Genz提出^[44]的Corner peak函数

$Y\left( x \right)=1+\sum\limits_{i=1}^{n}{{{a}_{i}}{{x}_{i}}^{-(n+1)}},$

(33)

其中，x_i(i=1,…,n)是在[0,1]上均匀分布的独立输入随机变量.随着参数a_i的增大，相应输入x_i的影响重要性增大，也即响应的Sobol主效应和总效应敏感度指标增大.

考察独立随机输入个数n=10以及参数a₁=0.01，a₂=0.02，… ,a₁₀=0.1，首先采用拟蒙特卡罗抽样(QMC)实现的蒙特卡罗方法计算Corner peak函数的Sobol主效应和总效应敏感度指标，用以作为参考值.然后，采用基于最小二乘回归和稀疏网格积分的多项式混沌展开法，计算多项式混沌展开为3阶时Corner peak函数的Sobol主效应和总效应敏感度指标，得到了与蒙特卡罗方法非常接近的计算结果(见表 8).最后，采用基于l₁优化的稀疏多项式混沌方法计算Sobol敏感度指标(见表 9)，使用了远少于最小二乘回归和稀疏网格积分的抽样样本数(p=3,M=200) ，就达到了相当的精度，这充分说明了稀疏多项式混沌方法在全局敏感度分析中的巨大效益.此外，展开阶数的增大，稀疏多项式混沌方法计算的结果越精确，说明了高阶稀疏多项式展开可以对光滑输出响应函数进行有效逼近.

4.4 Morris函数

考虑一个高维随机输入的Morris函数^[2]，其定义为

$\begin{align} & Y\left( x \right)={{\beta }_{0}}+\sum\limits_{i=1}^{20}{{{\beta }_{i}}{{w}_{i}}}+\sum\limits_{i＜j}^{20}{{{\beta }_{ij}}{{w}_{i}}{{w}_{j}}}+ \\ & \sum\limits_{i＜j＜l}^{20}{{{\beta }_{ijl}}{{w}_{i}}{{w}_{j}}{{w}_{l}}}+\sum\limits_{i＜j＜l＜s}^{20}{{{\beta }_{ijls}}{{w}_{i}}{{w}_{j}}{{w}_{l}}{{w}_{s}}}, \\ \end{align}$

(34)

${{w}_{i}}=\left\{ \begin{matrix} 2[1.1{{x}_{i}}/({{x}_{i}}+0.1)-0.5], & if~i=3,5,7, \\ 2({{x}_{i}}-0.5), & 其它. \\ \end{matrix} \right.$

(35)

这里，x_i(i=1,…,20)是在[0,1]上均匀分布的独立输入随机变量，式(34) 中的系数赋值如下

$\left\{ \begin{matrix} {{\beta }_{i}}=20, & for~i=1,\ldots ,10; \\ {{\beta }_{i}}={{\left( -1 \right)}^{i}}, & 其它; \\ {{\beta }_{ij}}=-15, & for~i,j=1,\ldots ,6; \\ {{\beta }_{ij}}={{\left( -1 \right)}^{i+j}}, & 其它; \\ {{\beta }_{ijl}}=-10, & for~i,j,l=1,\ldots ,5; \\ {{\beta }_{ijls}}=5, & for~i,j,l,s=1,\cdots ,4, \\ \end{matrix} \right.$

(36)

其它未赋值系数为零.显然，Morris函数没有关于均值、方差以及Sobol敏感度指标的解析表达式，需要采用蒙特卡罗方法对它们进行计算.Blatman和Sudret^[14]采用440 000个样本点计算了Morris函数的Sobol敏感度指标，给出了前10个随机输入的总效应灵敏度的置信区间，并以此为参考值来验证他们提出的稀疏多项式混沌方法的有效性.基于全尺度多项式混沌方法在高维随机输入下所面临的“维数灾难”问题，我们采用基于l₁优化的稀疏多项式混沌方法计算了Sobol敏感度指标，并将结果与拟蒙特卡罗抽样(QMC)实现的蒙特卡罗方法的结果进行对比，见表 10.表中给出了前10个随机输入的主效应和总效应敏感度指标，不仅蒙特卡罗法计算的总效应敏感度指标是在Blatman和Sudret给出的总效应敏感度置信区间内的，而且稀疏多项式混沌方法计算的结果也是在这个置信区间内.由此说明，在这种高维情况下，在回归法和稀疏网格积分法已经失效的情况下，基于l₁优化的稀疏多项式混沌方法能够有效地进行全局敏感度分析，而所需的配置样本点也远少于多项式混沌展开项数.

5 结论与展望

多项式混沌方法是一种最为有效地计算Sobol敏感度指标的方法，是进行基于方差分解的全局敏感度分析的重要手段.多项式混沌方法非常适用于全局光滑函数，在进行不确定度量化和全局敏感度分析时，具有谱方法的收敛速度，相比传统蒙特卡罗方法所需计算时间较少.但是，在处理高维数随机空间和进行高阶多项式混沌展开时将面临“维数灾难”问题.本文主要针对“维数灾难”问题，系统回顾了最小二乘回归法和稀疏网格积分法，以及介绍了基于l₁优化的稀疏重构技术(即压缩感知技术)，用以减少非嵌入式多项式混沌方法所需的样本配置点数目，采用几个典型的解析响应面模型来验证基于多项式混沌的全局敏感度分析方法的有效性.

多项式混沌展开法近年来在全局敏感度研究领域得到了长足的应用.多项式混沌方法不仅可以处理许多常见的概率分布类型作为随机输入变量，如正态分布、均匀分布、Gamma分布以及Beta分布，还可以通过针对任意概率密度函数(连续或离散测度)构造多项式混沌基函数，以此来处理任意概率分布的随机输入问题.此外，多项式混沌方法不仅可以用于基于方差分解的全局敏感度分析，还可以非常方便地推广到基于导数敏感度测度的全局敏感度研究^[45]，甚至还被作为高精度响应面模型来实现矩无关敏感度指标的计算^[46-47]，用于解决具有高偏斜性(high-skewed)和多峰性(multi-modal)随机输出响应的全局敏感度研究.

附录

● Smolyak稀疏网格方法

首先定义一维单变量随机函数的高斯积分，形式如下

${{l}^{(1)}}f={{\int }_{a}}^{b}f\left( \xi \right)p\left( \xi \right)d\xi \approx {{Q}^{(1)}}f\equiv \sum ni=1f({{\xi }^{(i)}}){{\omega }^{(i)}},$

(A.1)

这里，ξ⁽ⁱ⁾为高斯积分配置点，ω⁽ⁱ⁾为高斯权重值。积分配置点可以是嵌套类型，也可以是非嵌套类型。嵌套类型的积分配置点可以复用前一层级的积分点，可以减少高维随机函数积分所需配置点数目。图 1给出了嵌套类型的不同层级(level)上的Chebyshev积分配置点。

高维随机函数的张量积高斯积分可以表示成

$lf\approx {{Q}^{(N)}}f=\left( {{Q}_{l}}^{(1)}\otimes \cdots \otimes {{Q}_{l}}^{(1)} \right)f,$

(A.2)

这里,l表示积分层级，不同层级对于不同数目的配置点。例如，在不同层级上Chebyshev积分配置点数目是

${{n}_{l}}=({{2}^{l}})\pm 1.$

(A.3)

为了写出高维随机函数的Smolyak稀疏网格积分形式，定义

${{\Delta }_{l}}^{(1)}\equiv \left( {{Q}_{l}}^{(1)}-{{Q}_{l-1}}^{(1)} \right)f,{{Q}_{0}}^{(1)}f\equiv 0.$

(A.4)

由此，高维随机函数的Smolyak稀疏网格积分表达式为

${{Q}_{l}}^{(N)}f\equiv \sum\limits_{\left| l \right|\le l+N-1}{({{\Delta }_{{{l}_{1}}^{(1)}}}\otimes \cdots \otimes {{\Delta }_{{{l}_{N}}^{(1)}}})f}.$

(A.5)

这里，

$\left| l \right|\equiv \sum\limits_{i=1}^{N}{{}}$

图 2给出了在张量积和Smolyak稀疏网格下的二维Chebyshev积分配置点，显然可以发现Smolyak稀疏网格的所需配置点要少很多。

Chebyshev积分配置点下的Smolyak稀疏网格并不适用于一般概率分布的计算。均匀概率分布可以采用Gauss-Legendre积分配置点构造的稀疏网格，但是这种配置点却不是嵌套类型的；而适用于均匀概率分布的嵌套稀疏网格可以采用Gauss-Patterson积分配置点^[48]。需要指出的是，对于一维积分,n个Gauss-Legendre积分配置点可以达到2n-1阶的积分精度，而n个Gauss-Patterson积分配置点可以达到1.5n+0.5阶的积分精度。类似地，对于高斯概率分布可以采用Genz-Keister积分配置点^[49]构造的嵌套稀疏网格。