高斯过程模型的全局灵敏度分析的参数选择及采样方法

2015-01-07万华平任伟王宁波

振动工程学报 2015年5期

关键词：高斯全局灵敏度

万华平，任伟新，王宁波

（1.合肥工业大学土木与水利工程学院，安徽合肥230009；2.中南大学土木工程学院，湖南长沙410004）

高斯过程模型的全局灵敏度分析的参数选择及采样方法

万华平1，2，任伟新1，2，王宁波2

（1.合肥工业大学土木与水利工程学院，安徽合肥230009；2.中南大学土木工程学院，湖南长沙410004）

提出了基于全局灵敏度分析的有限元模型修正参数选择方法，考虑了参数整个变化空间的作用及参数间的相互作用，具有适用于参数不确定性大和无模型限制等优势。全局灵敏度分析常采用蒙托卡罗方法计算灵敏度指标，因此足够大的采样次数是获得可靠灵敏度指标的前提，但是同时会造成计算成本的增加。为此，采用高斯过程模型取代耗时的有限元模型用于降低计算成本，同时探讨了拉丁超立方抽样、Halton序列和Sobol序列3种空间采样方法用于全局灵敏度分析的计算效率，旨在选择一种高效的采样方法。最后，一桁架人行桥实例验证了有限元模型修正参数选择和采样选择方法。

参数选择；有限元模型修正；全局灵敏度分析；高斯过程模型；拉丁超立方抽样

引言

在工程领域，有限元模型被广泛地应用于基于模型的研究工作，一个能够准确表征结构行为的高精度有限元模型至关重要，而有限元模型修正是实现高精度有限元模型的方法之一。有限元模型修正包括确定目标函数、参数选择和数学优化估计修正参数这3个重要步骤［1］。复杂土木工程结构包含的参数较多，将所有参数用于模型修正是不现实的。如果不敏感的参数未被排除在修正参数之外，修正结果很可能会失去物理意义，因为自动的数学优化过程将会导致不敏感参数的大幅度修正。此外，剔除非敏感性参数会带来模型降阶，从而大大减轻了计算负担。因此，参数选择是结构有限元模型修正的关键问题之一。

目前广泛用于有限元模型修正参数筛选的方法为局部灵敏度分析［2］，它采用直接求导法或者有限差分法计算参数的局部灵敏度。局部灵敏度分析反映的是单个参数在名义值处的局部梯度信息，不能度量参数在整个取值范围的作用和参数间的共同作用。因此，局部灵敏度分析不适用于参数变化范围较大的情况。全局灵敏度分析克服了局部灵敏度分析的不足，而且该方法不受模型限制，也就是说它适用于任何模型，比如高度非线性模型等。全局灵敏度分析已广泛的应用于各种研究领域，比如结构动力分析［3］，管路系统设计［4］和结构可靠度分析［5］等。

但是，全局灵敏度分析计算花费高，因为蒙托卡罗方法需要大量采样以确保灵敏度指标的收敛。土木结构复杂，其通常模拟为高精度有限元模型，这样使得蛮力蒙托卡罗采样法（brute-force MCS，即直接基于有限元模型）不切实际。因此，降低全局灵敏度分析应用于复杂土木结构的计算花费显得尤为重要。为此，本文采用高斯过程模型取代耗时的高精度有限元模型以降低全局灵敏度分析的计算成本。高斯过程模型具有评估预报的不确定性，能很好地处理高维、小样本和强非线性等复杂问题［6］。近年来，高斯过程模型已经在工程领域得到广泛的应用［5，7-11］。另一方面，选择一种高效的采样方法也是降低全局灵敏度分析计算花费的有效途径。空间采样方法因其具有样本良好的分布均匀性而受到研究者的青睐。本文对比了拉丁超立方抽样、Halton序列和Sobol序列3种常用的空间采样方法在全局灵敏度分析中的计算效率，旨在为选择一种高效的采样方法提供参考。

1 高斯过程模型

高斯过程模型是基于贝叶斯思想构建的：高斯先验用来定义模型输出，通过对训练样本集的最大似然估计得到预报值的后验高斯分布，高斯过程模型完全由它的均值函数和协方差函数完全决定。对于均值函数一般采用零均值函数，因为缺乏对潜在函数的总体趋势先验知识［12］，并且零均值函数可以简化高斯过程模型的推导；对于协方差函数使用平方指数函数，该协方差函数具有使拟合的函数光滑且无限可导的优势

假设有n个观测值的训练样本集D=（X，T），其中X=［XT1，XT2，…，XTn］T，T=［t1，t2，…，tn］T。非观测点x*的预报值t*亦为高斯，其均值和方差如下

式中（·）T表示转置运算，C*=［C（X*，X1），C（X*，X2），…，C（X*，Xn）］T，C=C（X，X）和α=C-1T。

高斯过程模型的超参数Θ通过边缘释然函数最大化求得，即负对数边缘释然函数最小化

式中L（Θ）为负对数边缘释然函数，它的表达式及其对超参数的梯度如下［13］：

式中 |·|和tr（·）分别表示求行列式运算和对矩阵求迹运算。

式（4）的优化问题可以通过共轭梯度优化算法求解，但容易陷入局部最优，为此多初始值策略（Multi-Starting Point Strategy）用于防止优化过程的局部最优问题［11］。在决定使用构建的高斯过程模型代替有限元模型之前，留一法交叉验证（Leave-One-Out Cross-Validation）用来验证高斯过程模型的精度。

2 全局灵敏度分析

全局灵敏度分析的基本思想是将模型输出的总方差分解为由各参数及其相互作用贡献的子方差之和。具有d个独立随机参数的模型y（x）可分解为按维数递增形式的主效应和交互效应之和［14］

E（y）代表模型y（x）的均值，zi（xi）表示参数xi的主效应，zi，j（xi，j）为参数xi和xj间的相互作用，以此类推可得高阶项。

式（7）右边除E（y）之外的2d-1项互相独立［14］，因此可以获得如下表达式［15］

其中，

灵敏度指标定义于子方差与总方差的比值，如下

评价参数xi对模型输出作用大小的指标为一阶灵敏度指标Si和总灵敏度指标STi，STi由下式求得

由式（8）～（10）可知，求出Vi，V-i和V便可获得参数xi的一阶灵敏度指标Si和总灵敏度指标STi。方差项Vi，V-i和V可通过蒙托卡罗方法求取：

式中X（M1）和X（M2）为两个独立的N×d维抽样矩阵，（X（M2）-im，X（M1）im）表示用数组X（M1）的第i列向量替代数组X（M2）的第i列向量后形成的数组，同理（X（M1）-im，X（M2）im）表示用数组X（M2）的第i列向量替代数组X（M1）的第i列向量后形成的数组。

3 采样方法

3.1 拉丁超立方抽样

拉丁超立方抽样是由Mc Kay等人［16］提出的，它是一种基于分层抽样思想的方法。拉丁超立方抽样将变量的概率分布函数分成N（N为抽样次数）个互不重叠的等概率子区间，然后在每个子区间内分别进行独立随机采样，以确保在每个子区间精确地采样一次。拉丁超立方抽样的样本可由下式产生

式中π为1到N的独立随机排列数，U为［0，1］之间独立于π的随机数。

3.2 Halton序列

Halton序列［17］是vander Corput序列的扩展。它是通过将一系列整数表示成素数基的数位的形式，然后将这些数位反序排列再在前面加小数点得到的值。

设任意整数n和任意素数R，则整数n均有唯一的R进制表示

式中m=［lnn/lnR］，［·］表示取整运算。然后，将数位反序排列再在前面加小数点得到新的值如下

则d维的Halton序列可表示为

式中（R1，R2，…，Rd）为前d个素数。

3.3 Sobol序列

Sobol序列［18］是对Halton序列的重新排列，它是基于一组直接数而构造的。Halton序列以一组连续素数为基，而Sobol序列均以最小素数2为基，这样可以避免Halton方法产生高维序列需要计算大的素数的麻烦。

任意整数n以2为基表示如下

式中m=［lnn/lnR］，［·］表示取整运算，n0，…，nm取0或1。

Sobol序列的生成需借助不可约的本原多项式

式中pi（x）为第i维的本原多项式，Si为本原多项式pi（x）的度数，系数a1，i，a2，i，…，aSi-1，i取0或1。采用递推关系定义一组正数｛m1，i，m2，i，…，mSi，i｝，如下

式中 ⊕表示二进制按位异或运算，即1⊕0=0⊕1=0，1⊕1=0⊕0=1；mk，i（k=1，2，…，Si）为小于2i的奇数。

定义直接数νk，i=（k=1，2，…，Si），则基于直接数的Sobol序列可由下式求得

Antonov等人［19］提出了一种基于格雷码（Gray code）快速生成Sobol序列的方法。利用格雷码的性质，Sobol序列可由如下递归公式生成［20］

式中下标ck表示整数k二进制表示的最右零位。比如，对于0=（0000）2，1=（0001）2，2=（0010）23= （0011）2，可得到c0=1，c1=2，c2=1，c3=3。

4 实桥算例

4.1 桥梁描述及有限元模拟

风荷桥是位于湖南省常德市白马湖公园内的一座钢桁架人行桥（图1）。该桥宽3.5 m，长115 m （20 m+3×25 m+20 m），桥型布置如图2所示。该桁架人行桥上部结构由主桁、腹杆、横纵联、桥面系和栏杆扶手。主桁、腹杆、横纵联均由空心钢管构成，其尺寸分别为Φ273 mm×16 mm，Φ140 mm× 12 mm，Φ114 mm×3.5 mm。构成桥面系的实木尺寸为3.5 m（长度）×0.2 m（宽度）×0.15 m（厚度），实木之间间隙为0.02 m。栏杆立柱均匀焊接在主桁钢管上，立柱间间距为2.2 m。

图1 风荷桥Fig.1 Feheng Bridge

该桁架桥采用ANSYS进行有限元建模，桥面板采用质量单元（MASS21），其余构件均采用梁单元（BEAM44）。需要指出的是，通常被认作为附属结构的栏杆扶手在有限元模拟时采用了梁单元而非质量单元是因为该桥的栏杆扶手足够粗壮，其刚度贡献不应忽略。建立的有限元模型总共有3 035个节点和4 962个单元（4 832个梁单元，130个质量单元），如图3所示。

图2 桥型布置图（单位：mm）Fig.2 Conguration of Fenghe Bridge（Unit：mm）

图3 有限元模型Fig.3 Three-dimensional finite element model

4.2 参数选择及采样方法比较

该桥实测的前五阶模态频率为4.132，5.096，6.044，6.409，7.248 Hz；对应的初始有限元模型计算频率为4.352，5.369，6.365，6.721，7.634 Hz。假设桥梁模态频率的降低是由结构的刚度折减而引起。另外钢桥在实际服役过程中，钢材弹性模量易受工作环境影响（比如温度等）而发生改变。因此，本文选择弹性模量作为修正参数，根据构件所处的位置和作用划分为6个子块（如表1所示）。

表1 参数列表Tab.1 The list of parameters

采用基于频率残差的无量纲目标函数

式中fiGP和fiTest分别表示高斯过程预测频率和实测频率。由文献［21］可知，即使实测频率和计算频率相差不大，修正参数的幅度也可能很大。文献［21］中弹性模量修正量最大为30%，所以参数先验范围设定为名义值210 GPa的±30%改变量。

对前面叙述的3种空间采样方法，均采用足够大样本数5×106计算各参数的一阶灵敏度指标Si和总灵敏度指标STi，结果如图4所示（图中数字仅给出了有限位精度）。结果表明：

图4 参数对目标函数的灵敏度指标Fig.4 Sensitivity indices of parameters with respect to objective function

·参数E2，E3，E6与其他参数的相互作用明显，因为它们各自的总灵敏度指标明显大于一阶灵敏度指标；

·参数E2对目标函数起着主导作用，对目标函数的总效应占82.3%，这是因为主桁构成了该桥的主骨架，起着重要作用，因此参数E2应该确定为修正参数；

·参数E3，E6对目标函数影响也比较显著，所以它们均应该选择为修正参数；

·参数E1，E4，E5应该从修正参数组移除，因为它们对目标函数的影响极其微小。

接下来可考察3种空间采样方法用于计算灵敏度指标的效率。采样次数依次为103，5×103，10× 103，15×103，…，100×103，共21种采样。引入绝对误差（Absolute Error，AE）评价灵敏度指标计算的收敛效果如下

式中SMaxi（SMaxTi）为采样次数为5×106的一阶（总）灵敏度指标设为指标，SSami（SSamTi）表示采样次数为103，5×103，…，100×103的一阶（总）灵敏度指标，d为参数个数。3种空间采样方法效果对比如图5所示。

图5 采样方法计算收敛效果比较Fig.5 Convergence comparison of sampling methods

由图5可知，Sobol序列用于全局灵敏度分析计算效率最高，Halton序列次之，最后为拉丁超立方抽样。它们计算效率的高低可以从其均匀性优劣得到解释：拉丁超立方保证了一维的均匀性，但未考虑维度间的均匀性；Sobol序列和Halton序列为2种低差异性（Low Discrepancy）的拟蒙托卡罗采样（Quasi-Monte Carlo）方法，它们在维度间呈现出良好的均匀性；Sobol序列避免了Halton序列高维抽样时的集聚效应，前者的均匀性优于后者［22］。

5 结论

本文提出了基于全局灵敏度分析的有限元模型修正参数选择方法。相比传统的局部灵敏度分析方法，全局灵敏度分析方法具有量化参数整个变化范围作用，考虑了参数间的相互作用等优点。但是，全局灵敏度分析应用于复杂的土木结构计算花费高。为此，一方面，本文采用高斯过程模型用于降低全局灵敏度分析方法计算成本高；另一方，比较了拉丁超立方抽样、Halton序列、Sobol序列3种空间采样方法用于全局灵敏度分析的计算效率，对寻找一种高效的采样方法用于减少全局灵敏度分析计算花费具有指导意义。

［1］ Ren W X，Chen H B.Finite element model updating in structural dynamics by using the response surface method［J］.Engineering Structures，2010，32（8）：2 455—2 465.

［2］ Brownjohn J M W，Xia P Q，Hao H，et al.Civil structure condition assessment by fe model updating：Methodology and case studies［J］.Finite Elements in Analysis and Design，2001，37（10）：761—775.

［3］于德介，李睿.Sobol’法在非线性隔振系统灵敏度分析中的应用研究［J］.振动工程学报，2004，17（2）：210—213. Yu D J，Li R.Application of Sobol′method to sensitivity analysis of a non-linear passive vibration isolators ［J］.Journal of Vibration Engineering，2004，17（2）：210—213.

［4］陈刚，汪玉，毛为民，等.冲击载荷作用下舰艇管路系统全局参数灵敏度分析［J］.振动与冲击，2007，26 （3）：45—48. Cheng G，Wang Y，Mao Y M，et al.Global parameter sensitivity analysis of shipboard piping systems under shock loads［J］.Journal of Vibration and Shock，2006，26（3）：45—48.

［5］ Wang P，Lu Z，Tang Z.An application of the Kriging method in global sensitivity analysis with parameter uncertainty［J］.Applied Mathematical Modelling，2013，37（9）：6 543—6 555.

［6］ Simpson T W，Poplinski J D，Koch P N，et al.Metamodels for computer-based engineering design：Survey and recommendations［J］.Engineering with Computers，2001，17（2）：129—150.

［7］刘开云，刘保国，徐冲.基于遗传-组合核函数高斯过程回归算法的边坡非线性变形时序分析智能模型［J］.岩石力学与工程学报，2009，28（10）：2 128—2 134. Liu K Y，Liu B G，Xu C.Intelligent analysis model of slope nonlinear displacement time series based on genetic-Gaussian process regression algorithm of combined kernel function［J］.Chinese Journal of Rock Mechanics and Engineering，2009，28（10）：2 128—2 134.

［8］刘信恩，肖世富，莫军.高斯过程响应面法研究［J］.应用力学学报，2010，27（1）：190—195. Liu X E，Xiao S F，Mo J.A new response surface method based on Gaussian process［J］.Chinese Journal of Applied Mechanics，2010，27（1）：190—195.

［9］万华平，任伟新，魏锦辉.基于高斯过程响应面的结构有限元模型修正方法［J］.振动与冲击，2012，31 （24）：82—87. Wan H P，Ren W X，Wei J H.Strucutral finite element model updating based on Gaussian process response surface methodology［J］.Journal of Vibration and Shock，2012，31（24）：82—87.

［10］张冬冬，郭勤涛.Kriging响应面代理模型在有限元模型确认中的应用［J］.振动与冲击，2013，32（9）：187—191. Zhang D D，Guo Q T.Application of Kriging response surface in finite element model validation［J］.Journal of Vibration and Shock，2013，32（9）：187—191.

［11］Wan H P，Mao Z，Todd M D，et al.Analytical uncertainty quantification for modal frequencies with structural parameter uncertainty using a Gaussian process metamodel［J］.Engineering Structures，2014，75：577—589.

［12］Neal R M.Regression and classification using Gaussian process priors.In：Bernardo J，Berger J，Dawid A，et al.Bayesian Statistics 6［M］.Oxford University Press，1998，475—501.

［13］Rasmussen C E，Williams C K I.Gaussian Processes for machine learning［M］.MIT Press，2006.

［14］Efron B，Stein C.The jackknife estimate of variance ［J］.The Annals of Statistics，1981：586—596.

［15］Sobol I M.Sensitivity estimates for nonlinear mathematical models［J］.Mathematical Modeling and Computational Experiment，1993，1（4）：407—414.

［16］Mc Kay M D，Beckman R J，Conover W J.A comparison of three methods for selecting values of input variables in the analysis of output from a computer code ［J］.Technometrics，1979，21：239—245.

［17］Halton J H.On the efficiency of certain quasi-random sequences of points in evaluating multi-dimensional integrals［J］.Numerische Mathematik，1960，2（1）：84—90.

［18］Sobol I M.On the distribution of points in a cube and the approximate evaluation of integrals［J］.USSR Computational Mathematics and Mathematical Physics，1967，7：86—112.

［19］Antonov I A，Saleev，V M.An economic method of computing LPτ-sequences［J］.USSR Computational Mathematics and Mathematical Physics，1979，19（1），252—256.

［20］Bratley P，Fox B L.Algorithm 659：Implementing Sobol＇s quasirandom sequence generator［J］.ACM Transactions on Mathematical Software（TOMS），1988，14（1），88—100.

［21］Jaishi B，Ren W X.Structural finite element model updating using ambient vibration test results［J］.Journal of Structural Engineering，2005，131（4）：617—628.

［22］Cheng J，Druzdzel M J.Computational investigation of low-discrepancy sequences in simulation algorithms for bayesian networks［A］.Proceedings of the Sixteenth Conference on Uncertainty in Artificial Intelligence ［C］.Morgan Kaufmann Publishers，Stanford，CA，2000，72—81.

WAN Hua-ping1，2，REN Wei-xin1，2，WANG Ning-bo2
（1.School of Civil and Hydraulic Engineering，Hefei University of Technology，Hefei 230009，China 2.School of Civil Engineering，Central South University，Changsha 410004，China）

This paper proposes a global sensitivity analysis（GSA）approach for parameter selection in finite element modal updating.The GSA has the capability to quantify the effects of individual parameters over their entire space and interaction effects among parameters.This approach has a couple of advantages，such as capability to high degree of parameter uncertainty and model independence.However，it is still computationally intensive because a large number of model evaluations for Monte Carlo simulation（MCS）are needed to obtain a confident estimate of the sensitivity indices（SIs）.Therefore，the fast-running Gaussian process model is adopted to replace the time-consuming finite element model to alleviate the computational burden.In addition，this study explores the efficiency of three space-filling sampling methods，i.e.，Latin hypercube sampling，Halton sequence and Sobol sequence，in the SIs evaluation，which has great value for choosing an efficient sampling method.At last，a real-world truss pedestrian bridge is used to demonstrate the proposed framework.

parameter selection；finite element modal updating；global sensitivity analysis；Gaussian process model；Latin hypercube sampling