方圣恩 張秋虎 林友勤 張笑華

(1.福州大學土木工程學院，福州 350116)(2.合肥工業(yè)大學土木與水利學院，合肥 230009)

引言

數(shù)值和數(shù)學模型已廣泛應用于復雜工程問題的模擬和計算，以進行結(jié)構(gòu)設計和優(yōu)化、靜動力分析、缺陷或損傷診斷、安全性評估等［1］.但現(xiàn)實結(jié)構(gòu)中總存在著不確定性，在分析大量結(jié)構(gòu)參數(shù)時容易導致分析結(jié)果的失真和矛盾，同時又需耗費大量的計算成本［2］.因此，靈敏度分析(sensitivity analysis，縮寫SA)［3］對確定模型參數(shù)來說是十分必要的，其通過量化參數(shù)不確定性對響應變異性(比如方差)的貢獻來識別參數(shù)重要性［4］.具體分為參數(shù)主效應(main effects)分析和相互效應(interaction effects)分析.

目前已有的SA方法大多基于數(shù)學或統(tǒng)計學方法，各有其適用范圍［5］.可以一次單獨分析一個參數(shù)，也可以采用隨機抽樣的方式同時分析多個參數(shù)［6］.然而對同一個問題，不同的SA方法可能給出不同的判斷結(jié)果.Saltelli等［3］將SA方法分為散點圖(scatterplot)法、導數(shù)(derivatives)法、回歸系數(shù)(regression coefficients)法及方差(variance-based)法等四大類.各種方法分別對應于線性或非線性模型的基本假設，有著各自的優(yōu)缺點.散點圖可以將參數(shù)和響應間關(guān)系通過可視化方式體現(xiàn)出來，非常直觀，但是判斷上往往帶有主觀性且無法自動獲取相關(guān)靈敏度指標［7］.通過偏微分求導的方式則要求參數(shù)—響應間關(guān)系可以函數(shù)化，該方法容易理解且計算效率高，但僅能對參數(shù)設計點附近小范圍內(nèi)的不確定性進行分析，無法提供靈敏度的全局信息［8］.此外，通過對參數(shù)—響應樣本進行線性或非線性回歸，可以搜索整個參數(shù)不確定性空間，得到的標準回歸系數(shù)往往就是靈敏度指標［9］.但為了得到足夠的分析精度，回歸法往往需要大量的樣本，導致對復雜模型來說計算量大為增加，實用性不強.最后，方差法可以很好地了解模型的靈敏度組成模式，尤其適用于模型的線性、單調(diào)性和可疊加性未知的情況［10，11］.方差法基于條件方差的計算，分解出各參數(shù)對響應方差的獨立影響，同時還能估計高階相互效應的影響.但這種方法實現(xiàn)過程比較復雜，計算量大，不利于實際應用.因此，可以在SA過程引入meta-model(如Kriging模型)來快速進行參數(shù)—響應間的不確定性傳播，以大幅提高計算效率［12，13］.但傳統(tǒng)meta-model對參數(shù)總體效應和高階相互效應的估計精度往往不足.由此可見，一個好的SA方法要具備全局性、計算高效性和實用性強等特點.

有鑒于此，本文結(jié)合隨機響應面(stochastic response surface，縮寫SRS)，利用偏導方式求解得到參數(shù)靈敏度系數(shù)，以量化參數(shù)對響應變異性(單位響應變化)的貢獻率.所提出方法通過數(shù)值梁進行驗證，并與方差分析(analysis of variance，縮寫ANOVA)法進行對比，驗證了本文方法的正確性.

1 隨機響應面理論

隨機響應面可以看作是對確定性響應面理論的擴展，前者通過基于Hermite多項式的多項式混沌展開式(polynomial chaos expansion)來建立不確定性參數(shù)x和響應y之間的聯(lián)系［14］:

式中ξ={ξi}表示服從標準正態(tài)分布N(0，1)的標準隨機變量，是原參數(shù)向量x的變換;ai為表達式各項的系數(shù)，通過回歸方式求解;Γp(·)表示多維p階Hermite多項式［14］:

將x轉(zhuǎn)換為ξ可以保證Γp(·)的正交性，即E [ Γp·Γq]=0(p≠q).對SA來說，這種轉(zhuǎn)換能保證分析過程同等對待(無量綱化)不同類型的參數(shù).然后通過概率配點法［14］(probabilistic collocation method)計算擬合SRS所需的樣本，再利用回歸分析求得ai，最后得到SRS的顯式多項式表達式.相關(guān)理論的詳細介紹可以參閱文獻［14］，本文不再具體闡述.

2 靈敏度分析的隨機響應面法

本節(jié)中基于公式(1)推導參數(shù)的靈敏度系數(shù)矩陣.首先對ξ(而非x)求偏導數(shù)，以此同等對待不同類型的參數(shù).同時考慮到偏導數(shù)問題往往難以求解，需要結(jié)合數(shù)值分析方法，使得求解過程復雜化.而SRS表達式中的參數(shù)本身就是隨機參數(shù)，且多項式分解后對ξ求偏導容易實現(xiàn).此外，求解過程考慮了參數(shù)不確定性的完整空間，是一種全局方法.這里以常用的二階SRS為例:

響應y對ξi偏導的期望值為:

式中隨機標準變量的均值為0，即E(ξi)=0.因此，向量y=(yk)(k=1，2，…，m)的偏導矩陣θ可以表示為:

所對應的期望值矩陣為:

式中矩陣元素aki即式(1)中一階項的系數(shù)，其代表了參數(shù)對響應的主效應.該結(jié)論可以通過3參數(shù)的2階SRS來證明，其中xi=μi+σiξi(i=1，2，3).即:

由此求得:

同時，y對原參數(shù)(x1，x2，x3)的偏導數(shù)為:

式中σi(i=1，2，3)是參數(shù)xi的標準差.因此容易求得:

因為xi的概率分布是已知的，E［?y/?xi］實際上是對應于xi均值的常數(shù)(假設為bi)，所以可以進一步得到:

由上式可見，系數(shù)ai實際上綜合考慮了xi的均值和標準差，其中均值可以看作是參數(shù)xi的確定性部分，而標準差則代表了不確定性部分，相當于傳統(tǒng)靈敏度分析中的“對參數(shù)擾動”，而這種“擾動”必須基于參數(shù)的某個設計點，即本文中的均值.若是標準差脫離了均值進行分析，那就失去了SA的意義.最后值得一提的是，系數(shù)ai是一種穩(wěn)定估計，即更高階(如3階)SRS表達式中一階項所對應的ai幾乎不會有變化［14］，這一特點對本文提出的SA方法來說是有利的.因為在構(gòu)建metamodel時，可能無法確定模型的階數(shù)，若是基于高階和低階模型得到的SA結(jié)果是不同的，那么靈敏度分析方法也就失去了可靠性.但是采用SRS則不存在這個問題，因此其具有更優(yōu)的適用性.

3 靈敏度分析的方差分析法

為了進行比較，本文同時采用了ANOVA法估計參數(shù)不確定性對響應變異性的影響.ANOVA通過估計不同樣本組的組間與組內(nèi)均方(mean of square)來研究樣本不確定性或誤差的來源［15，16］，其中樣本來源于服從正態(tài)分布的總體.ANOVA法可以分離各參數(shù)的效應，以估計參數(shù)對模型總體變異性的獨立影響，或不同參數(shù)組合相互效應對模型總體變異性的影響.具體實施上，本文采用了2k因子設計［15，16］(2kfactorial design，縮寫2kFD)來實現(xiàn)參數(shù)的SA.該設計基于ANOVA，每個參數(shù)僅有2個水平(可看作上下界值)，通過實驗設計方法得到不同的參數(shù)和水平組合，在通過試驗或數(shù)值計算求得樣本;然后通過回歸分析擬合樣本得到一個線性模型，模型表達式的各項系數(shù)即相當于參數(shù)不確定性的效應系數(shù).

4 算例

本文基于數(shù)值懸臂梁(圖1)來驗證所提出方法的可行性和可靠性.采用數(shù)值算例可以避免試驗數(shù)據(jù)中混雜了其他的不確定性來源(比如環(huán)境噪聲、人為誤差等)，這些不確定性可能導致分析過程無法辨別響應變異性是否源于參數(shù)的不確定性，導致分析目標不明確.

圖1 數(shù)值懸臂梁示意圖Fig.1 Schematic diagram of the numerical cantilever beam

所采用的懸臂梁長2 m，截面尺寸0.2 m×0.2 m;材料參數(shù)為楊氏模量(E)30 GPa，密度(D)2400 kg/m3，泊松比(P)0.2.上述幾何參數(shù)與材料特性均假設為名義值(均值).梁的有限元模型被劃分為20個相同梁單元，并假設單元10包含了4個不確定性參數(shù)E、I(截面慣性矩)、D、P，以此分析參數(shù)不確定性對梁前5階振動頻率的影響.此外，由圖2可見，梁的5階振動模態(tài)中包括了4階彎曲模態(tài)和1階軸向變形模態(tài).

圖2 懸臂梁振動模態(tài)Fig.2 Vibrational modes of the cantilever beam

本算例設計了兩種情形:1)參數(shù)不確定性水平相同，即每個參數(shù)的名義值作為均值，而表示不確定性的標準差設為1%;2)參數(shù)不確定性水平不同，即假設I的標準差為2%，其他3個參數(shù)的標準差仍為1%.第一種情況可以在同等條件下估計參數(shù)靈敏度，而第二種情況則考慮到現(xiàn)實結(jié)構(gòu)中不同參數(shù)的不確定性水平可能不同，因此需要區(qū)別對待.

4.1 參數(shù)不確定性水平相同

首先建立包含4個參數(shù)的2階SRS，即通過概率配點法得到20個樣本，然后基于回歸分析求得SRS的表達式系數(shù).由第2節(jié)可知，表達式中1階項的系數(shù)反映了各參數(shù)對頻率的主效應.與此同時，為了進行比較，還利用2kFD進行SA，所需的樣本數(shù)為24=16.兩種方法的SA結(jié)果分別列于表1和2，其中為了便于比較，表中數(shù)值為當頻率發(fā)生單位變化時，各參數(shù)對頻率變異性的貢獻率.

表1 參數(shù)不確定性對頻率變異性的貢獻率(情形1:SRS)Table 1 Parameter uncertainty percentage contributions to frequency variability(scenario 1:SRS)

表2 參數(shù)不確定性對頻率變異性的貢獻率(情形1:ANOVA)Table 2 Parameter uncertainty percentage contributions to frequency variability(scenario 1:ANOVA)

由表可見，兩種方法給出了完全一致的分析結(jié)果，驗證了本文方法的正確性.對第1階頻率而言，E和I的貢獻率相同，且比D高10%左右，說明聯(lián)系單元10截面抗彎剛度的參數(shù)對頻率的影響比密度大;而在第2、5階頻率上，E、I、D三者的貢獻率基本相同.同時對第3階軸向變形模態(tài)來說，截面慣性矩I此時不起作用，因為其僅和抗彎剛度相關(guān)，與截面軸向剛度無關(guān);而且在第4階模態(tài)上，由于單元10處于模態(tài)節(jié)點處，此時I的影響也幾乎為零.值得注意是，P對所有5階頻率都沒有任何影響.上述觀察結(jié)果可以通過懸臂梁振動頻率的理論公式［17］來解釋:

式中L和A分別表示梁長和截面積.由上式可見，對彎曲頻率來說，E、I、D的影響應該是相同的;而P對頻率則無任何影響.算例分析中，第1階頻率中D的貢獻率可能受到單元10在1階模態(tài)振型中位置的影響，此時E、I的影響會變大一點.這點推測可以在第4階頻率的分析中得到證實:該階模態(tài)振型中單元10正好處于模態(tài)節(jié)點處，使得此單元的抗彎剛度EI不發(fā)揮作用(主要是截面慣性矩I不起作用)，因此I的影響幾乎為0.而同樣作為材料參數(shù)，此時D的貢獻率卻是E的2.5倍，說明E的影響也被大大削弱了.最后，對第3階軸向振動模態(tài)來說，僅E和D對頻率的變異性有貢獻，這點和理論公式一致.由上述分析可見，SA可以有效地對參數(shù)不確定性的影響進行量化，從而發(fā)現(xiàn)對分析模型重要的參數(shù).理論公式雖然能一定程度上反映問題，但容易造成誤判(因為理論公式中的參數(shù)是對整根梁而言的，而不是針對某個局部單元)，特別是對復雜工程結(jié)構(gòu)來說，通常無法得到頻率的理論公式，此時SA方法更能體現(xiàn)其價值.

最后要說明的是，表1、2僅給出了參數(shù)的主效應，因為本算例中參數(shù)的相互效應影響非常小(EI、ED和ID對5階頻率的總體貢獻率僅分別為0.40%、0.28%、0.02%、0.51%和0.31%)，所以未考慮.由此可以認為頻率的總體變異性是各參數(shù)不確定性獨立影響的疊加.

4.2 參數(shù)不確定性水平不同

本小節(jié)中根據(jù)E、I、D的不同不確定性水平，重新構(gòu)建了SRS和2kFD進行分析，結(jié)果列于表3和4.基于4.1節(jié)的分析結(jié)果，本分析中不再考慮P.由表可見，在I增加了1倍的不確定性后，第1、2、5階頻率中I的貢獻率也變成了E和D的2倍.該結(jié)果是合理的，因為在頻率理論公式中3個參數(shù)的權(quán)重是相同的.同時對剩下的兩階頻率而言，此時I依然沒有影響，即便其不確定性水平提高了1倍，原因還是I對這兩階模態(tài)不相關(guān)，這從另一方面也說明了上節(jié)的分析結(jié)果是正確的.最后，對第4階頻率來說，D的貢獻率仍然是E的2.5倍，說明其沒有受到I不確定性水平改變的影響，同時也證明了上節(jié)的分析結(jié)果是正確的.

4.3 方法適用性討論

由算例結(jié)果可知，SRS法給出的SA結(jié)果是可靠和準確的.但也存在一個問題，即既然ANOVA法可以給出準確的分析結(jié)果，那么提出SRS法是否必要?為此從以下幾點進行討論.首先，SRS和ANOVA方法都要求參數(shù)不確定性服從正態(tài)分布，從這點上說，二者皆屬于概率SA法.但要注意的是，SRS可以直接給出不確定參數(shù)和響應間關(guān)系的隨機表達式，應用上更簡便，同時可以表示參數(shù)和響應的非線性關(guān)系;而2kFD建立的是線性模型，要求參數(shù)和響應間不存在強非線性.同時，2kFD分析時僅考慮參數(shù)的上下界，而對界限內(nèi)的概率分布情況是一種弱要求，這對ANOVA分析的結(jié)果是有一定的影響的;而SRS理論則是完全基于參數(shù)概率分布假設的，分析過程更嚴格.最為重要的是，對擁有較多參數(shù)的復雜結(jié)構(gòu)來說，2kFD所需的樣本數(shù)呈指數(shù)級增長，使得計算成本激增.比如10個參數(shù)情況，2kFD至少需要210=1024個樣本;而相同參數(shù)數(shù)目下，SRS僅需要132個樣本，計算量上要小得多.由上述分析可見，兩種方法各有其自身的適用范圍，在某些情況下，SRS法更具優(yōu)勢.

表3 參數(shù)不確定性對頻率變異性的貢獻率(情形2:SRS)Table 3 Parameter uncertainty percentage contributions to frequency variability(scenario 2:SRS)

表4 參數(shù)不確定性對頻率變異性的貢獻率(情形2:ANOVA)Table 4 Parameter uncertainty percentage contributions to frequency variability(scenario 2:ANOVA)

5 結(jié)論

本文針對工程中不確定性參數(shù)的靈敏度分析問題，提出了一種隨機響應面法，即基于對SRS表達式的偏導求解，推導出參數(shù)靈敏度系數(shù)，實現(xiàn)對參數(shù)不確定性影響的量化.文中基于一根包含不確定單元參數(shù)的數(shù)值懸臂梁來驗證所提出的方法，并與ANOVA法進行對比.分析結(jié)果表明，本文方法可以準確地估計參數(shù)靈敏度，并給出各參數(shù)對頻率響應變異性的貢獻率，使得靈敏度分析結(jié)果更客觀.同時，SRS的構(gòu)建過程無需大量樣本，在計算成本上有著一定的優(yōu)勢.最后，所提出的方法可以進一步應用于復雜結(jié)構(gòu)的參數(shù)靈敏度分析上.

1 Saltelli A，Tarantola S，Campolongo F，Ratto M.Sensitivity analysis in practice:a guide to assessing scientific models.Chichester:John Wiley＆Sons Ltd，2004

2 Hornberger G M，Spear R C.An approach to the preliminary analysis of environmental systems.Journal of Environmental Management，1981，12:7～18

3 Saltelli A，Ratto M，Andres T，Campolongo F，Cariboni J，Gatelli D，Saisana M，Tarantola S.Global sensitivity analysis:the primer.Chichester:John Wiley＆Sons Ltd，2008

4 Chan K，Saltelli A，Tarantola S.Sensitivity analysis of model output:variance-based methods make the difference.In:Proceedings of the 1997 Winter Simulation Conference，Atlanta，USA，1997

5 Iman R L，Helton JC.An investigation of uncertainty and sensitivity analysis techniques for computer models.Risk Analysis，1988，8(1):71～90

6 Hamby D M.A review of techniques for parameter sensitivity analysis of environmental models.Environmental Monitoring and Assessment，1994，32(2):135～154

7 Chan Y H，Correa CD，Ma K L.The generalized sensitivity scatterplot.IEEE Transactions on Visualization and Computer Graphics，2013，19(10):1768～1781

8 Punzo V，Ciuffo B.Sensitivity analysis of car-following models:methodology and application.In:Proceedings of the 90th Transportation Research Board Annual Meeting，Washington，2011

9 Ratto M，Pagano A，Young P.State dependent parameter metamodelling and sensitivity analysis.Computer Physics Communications，2007，177(11):863～876

10 Cukier R I，Levine H B，Shuler K E.Nonlinear sensitivity analysis of multiparameter model systems.Journal of Computational Physics，1978，26(1):1～42

11 Saltelli A.Making best use of model evaluations to compute sensitivity indices.Computer Physics Communications，2002，145(2):280～297

12 Storlie C，Helton J.Multiple predictor smoothing methods for sensitivity analysis:description of techniques.Reliability Engineering and System Safety，2008，93:28～54

13 Ciuffo B，Punzo V，Qualietta E.Kriging meta-modelling to verify traffic micro-simulation calibration methods.In:Proceedings of the 90th Transportation Research Board Annual Meeting，Washington，2011

14 Isukapalli S S.Uncertainty analysis of transport-transformation models［PhD Thesis］.New Brunswick Rutgers:The State University of New Jersey，1999

15 Montgomery D C.Design and analysis of experiments(6th edition).New York:J.Wiley＆Sons，2004

16 Myers R H，Montgomery D C，Anderson-Cook C M.Response surface methodology:process and product optimization using designed experiments(3rd edition).New York:John Wiley＆Sons，2009

17 Chopra A K.Dynamics of structures:theory and applications to earthquake engineering(4th Edition).Upper Saddle River:Prentice Hall，2011