不確定性參數(shù)靈敏度分析的隨機(jī)響應(yīng)面法*

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

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

引言

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

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

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

1 隨機(jī)響應(yīng)面理論

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

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

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

2 靈敏度分析的隨機(jī)響應(yīng)面法

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

響應(yīng)y對ξi偏導(dǎo)的期望值為:

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

所對應(yīng)的期望值矩陣為:

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

由此求得:

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

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

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

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

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

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

4 算例

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

圖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ù)與材料特性均假設(shè)為名義值(均值).梁的有限元模型被劃分為20個(gè)相同梁單元，并假設(shè)單元10包含了4個(gè)不確定性參數(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

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

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

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

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

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

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

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

最后要說明的是，表1、2僅給出了參數(shù)的主效應(yīng)，因?yàn)楸舅憷袇?shù)的相互效應(yīng)影響非常小(EI、ED和ID對5階頻率的總體貢獻(xiàn)率僅分別為0.40%、0.28%、0.02%、0.51%和0.31%)，所以未考慮.由此可以認(rèn)為頻率的總體變異性是各參數(shù)不確定性獨(dú)立影響的疊加.

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

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

4.3 方法適用性討論

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

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

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

5 結(jié)論

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