鞏祥瑞， 呂震宙， 孫天宇， 張雷雷， 封雷

(1. 西安現(xiàn)代控制技術(shù)研究所, 西安 710065； 2. 西北工業(yè)大學(xué) 航空學(xué)院, 西安 710072)

靈敏度分析方法主要研究結(jié)構(gòu)系統(tǒng)的輸入變量對輸出響應(yīng)的影響程度[1]，已廣泛應(yīng)用于可靠性工程[2]、危險分析[3-4]、環(huán)境科學(xué)[5]、計算物理化學(xué)[6-7]等領(lǐng)域。靈敏度分析一般分為局部靈敏度分析方法和全局靈敏度分析方法。局部靈敏度分析方法研究模型輸入變量的分布參數(shù)對結(jié)構(gòu)系統(tǒng)輸出響應(yīng)的影響程度[8]，只能反映輸入變量在某一名義值處對輸出響應(yīng)的影響，具有很大的局限性[9]。

全局靈敏度分析，也稱為重要性測度分析，研究一個或多個輸入變量對模型輸出響應(yīng)的影響程度，能夠反映輸入變量對結(jié)構(gòu)系統(tǒng)輸出響應(yīng)的平均影響程度，因此，該分析方法在實際工程中得到了更廣泛的應(yīng)用。重要性測度分析方法的核心思想是對影響結(jié)構(gòu)系統(tǒng)輸出響應(yīng)的輸入變量進(jìn)行重要性排序，從而得到對結(jié)構(gòu)系統(tǒng)輸出響應(yīng)影響較大的輸入變量，為結(jié)構(gòu)系統(tǒng)的設(shè)計改進(jìn)提供理論依據(jù)。近年來，研究人員已經(jīng)提出了很多重要性測度分析方法。Saltelli 和 Helton提出了無參的方法，但缺乏模型的獨(dú)立性[9-10]。Sobol[11]和 Iman[12]等提出了基于方差的重要性測度分析方法和相應(yīng)的求解算法，該方法隱含地假設(shè)了方差可以充分描述結(jié)構(gòu)系統(tǒng)輸出響應(yīng)的不確定性。然而，相關(guān)學(xué)者指出隨機(jī)變量的任何一階矩只能反映其部分信息，必然會丟失一些信息[13]。為了分析輸入變量在其整個分布區(qū)域上對結(jié)構(gòu)系統(tǒng)輸出響應(yīng)的影響程度，Borgonovo[14]、Liu[15]和 Li[16]等又分別提出了矩獨(dú)立重要性測度分析(Moment-independent Importance Measures analysis，M-IM)方法。

目前，Borgonovo[14]所提出的矩獨(dú)立重要性測度分析方法應(yīng)用最為廣泛。但在可靠性分析中，結(jié)構(gòu)系統(tǒng)的失效概率通常是工程技術(shù)人員所關(guān)注的焦點。Li等[16]又提出了基于失效概率的矩獨(dú)立重要性測度分析方法，該方法能夠計算輸入變量在其分布區(qū)域的任一固定點時，結(jié)構(gòu)系統(tǒng)的無條件失效概率和條件失效概率的平均偏差。然而，在實際工程中將輸入變量固定于某一點，這是很難實現(xiàn)的，只能使輸入變量在其分布區(qū)域的某一區(qū)間上縮減變化。

為了解決這一問題，Wei等[17]提出了一種新的輸入變量縮減區(qū)間的獲取方法，該方法的基本思想是：首先產(chǎn)生2個隨機(jī)數(shù)q1和q2，q1和q2服從均勻分布，q1～U(0,1)，q2～U(q1,1)，然后根據(jù)輸入變量累積分布函數(shù)的逆函數(shù)F-1(·)，得到輸入變量的縮減區(qū)間[F-1(q1),F-1(q2)]。該方法的優(yōu)點是很容易獲取輸入變量的縮減區(qū)間，然而它并沒有考慮輸入變量縮減區(qū)間獲取的任意性和等可能性。

本文提出了一種新的基于失效概率的矩獨(dú)立重要性測度分析方法，并給出了一種更加合理的輸入變量縮減區(qū)間的獲取方法。同時，引入自適應(yīng)超球重要抽樣(Adaptive Radial-Based Importance Sampling，ARBIS)方法[18-19]，構(gòu)建所提新指標(biāo)的高效求解算法。應(yīng)用數(shù)值算例和工程算例，驗證本文所提新指標(biāo)的合理性和求解算法的高效性。

1 新的矩獨(dú)立重要性測度指標(biāo)

1.1 新指標(biāo)的定義

EQ(E(IF)-E(IF|Xi∈Ui(Q)))2

(1)

式中：Pf為非條件失效概率；Q為分位數(shù)矩陣，用來生成輸入變量Xi的縮減區(qū)間；當(dāng)Xi的分布區(qū)間縮減到其所有可能的子區(qū)間Ui(Q)時，就能夠得到條件失效概率Pf|Xi∈Ui(Q)；EQ(·)為對Q求期望；IF為指示函數(shù)，當(dāng)X位于失效域(即g(X)≤0)時，IF(x)=1，否則IF(x)=0。

1.2 輸入變量縮減區(qū)間獲取方法

為了等概率地獲取輸入變量所有可能的縮減區(qū)間，首先定義一個分位數(shù)矩陣Hk(3×2)：

(2)

式中：hk(1)和hk(2)為[0,1]區(qū)間上的2個隨機(jī)數(shù)，滿足0≤hk(1)

Hk(k=1,2,…,N)組成一個(3N×2)的分位數(shù)矩陣Q：

(3)

式中：Hk(k=1,2,…,N) 是被等概率獲取的，例如P{Hk}=1/N。

為了更好地理解式(3)，可表示為

(4)

式中：q3k-2(1)=0，q3k-2(2)=q3k-1(1)=hk(1)，q3k-1(2)=q3k(1)=hk(2)，q3k(2)=1 (k=1,2,…,N)。

因此，根據(jù)分位數(shù)矩陣Q就可以定義縮減區(qū)間矩陣U為

(5)

對于輸入變量Xi，其縮減區(qū)間可表示為

(6)

上述是一種新的區(qū)間劃分技術(shù)來獲取輸入變量在其分布區(qū)域上所有可能的縮減區(qū)間。本文所提出的這種新的區(qū)間獲取方法具有2個優(yōu)點：①輸入變量所有可能的縮減區(qū)間能夠被等可能地獲取，從而就能夠更加合理地計算當(dāng)輸入變量在其分布區(qū)域的縮減區(qū)間上變化時對結(jié)構(gòu)系統(tǒng)失效概率的平均影響程度；②全期望公式在這種區(qū)間劃分技術(shù)中成立，基于全期望公式，第2節(jié)中就能夠?qū)⒈疚奶岢龅男碌幕谑Ц怕实木鬲?dú)立重要性測度指標(biāo)轉(zhuǎn)化為更加便于計算的基于方差的重要性測度指標(biāo)。

2 新指標(biāo)與基于方差的指標(biāo)的關(guān)系

由于1.2節(jié)中提出的輸入變量縮減區(qū)間獲取方法中全期望公式成立，就能夠推導(dǎo)出新的矩獨(dú)立重要性測度指標(biāo)與基于方差的重要性測度指標(biāo)間的關(guān)系，同時得到一種新的基于方差的矩獨(dú)立重要性測度指標(biāo)表示方式。

輸入變量在縮減區(qū)間上的全期望公式表示為

E(IF)=EQ(E(IF|Xi∈U(Q)))

(7)

證明過程詳見附錄A。

VQ(E(IF|Xi∈Ui(Q)))

(8)

從式(8)中可以看出，新的矩獨(dú)立重要性測度指標(biāo)能夠很容易地轉(zhuǎn)換為基于方差的重要性測度指標(biāo)，從而新的矩獨(dú)立重要性測度指標(biāo)就能夠采用基于方差的重要性測度分析方法來進(jìn)行求解。

3 新指標(biāo)的高效求解算法

新的矩獨(dú)立重要性測度指標(biāo)傳統(tǒng)的計算方法是雙層重復(fù)抽樣蒙特卡羅(Double-Loop-Repeat-Set Monte Carlo，DLRS MC)[17]方法，該方法先用N個樣本點計算非條件失效概率，然后在每一個子區(qū)間上，再產(chǎn)生N個新的樣本點來計算條件失效概率，最后計算非條件失效概率和條件失效概率差異的平均值。

DLRS MC方法可以獲得高精度的結(jié)果值，但其計算量很大，計算效率較低。本節(jié)提出了一種新的求解算法——ARBIS方法，旨在對新的矩獨(dú)立重要性測度指標(biāo)進(jìn)行高效求解。

3.1 自適應(yīng)超球重要抽樣方法

ARBIS方法最初是用來計算失效概率的，而計算新的基于失效概率的矩獨(dú)立重要性測度指標(biāo)的核心就是計算非條件失效概率和條件失效概率。因此，本節(jié)基于ARBIS方法構(gòu)建新指標(biāo)的求解算法。

ARBIS方法[18-19]的基本思想是自適應(yīng)地尋找到一個落在功能函數(shù)安全域內(nèi)的最大的超球。當(dāng)抽樣點落入超球內(nèi)部時，就將該抽樣點歸類為安全點，不必再代入結(jié)構(gòu)系統(tǒng)功能函數(shù)進(jìn)行求解。因此，該方法在計算結(jié)構(gòu)系統(tǒng)輸出響應(yīng)的失效概率時，計算量將大幅度降低。同時， ARBIS方法采用自適應(yīng)的方法來獲取滿足條件的最大超球半徑，相比于傳統(tǒng)的梯度搜索算法，效率更高，穩(wěn)健性也更好。

任何一個隨機(jī)變量通過合適的變換方法，都可以轉(zhuǎn)換到標(biāo)準(zhǔn)正態(tài)空間中[20-21]。因此，以下步驟都假定在標(biāo)準(zhǔn)正態(tài)空間中進(jìn)行。為了更好地理解ARBIS方法，給出該方法的一些關(guān)鍵步驟。

步驟1設(shè)定初始化超球半徑β0。 為了確保超球與失效域相交，初始半徑β0可通過式(9)來確定：

(9)

步驟2根據(jù)隨機(jī)變量的分布類型，通過轉(zhuǎn)化得到標(biāo)準(zhǔn)正態(tài)空間的樣本點uk(k=1,2,…,N)。

(10)

圖1 自適應(yīng)策略獲取最優(yōu)化半徑Fig.1 Adaptive strategy for obtaining optimal radius

步驟6重復(fù)步驟4和步驟5，直到收斂條件被滿足[22]，即可得到最優(yōu)的超球半徑βopt。

3.2 采用自適應(yīng)超球重要抽樣方法求解新指標(biāo)

當(dāng)輸入變量Xi縮減到區(qū)間Ui(Q)時，Xi就是一個截斷變量，則所有輸入變量的聯(lián)合概率密度函數(shù)可表示為

(11)

式中：K為截斷系數(shù)，可表示為

(12)

對于式(12)的估計，并不需要計算結(jié)構(gòu)系統(tǒng)的功能函數(shù)值，其計算量可以忽略。因此，條件失效概率Pf|Xi∈Ui(Q)可重新表示為

(13)

從式(13)中可以看出，條件失效概率Pf|Xi∈Ui(Q)的計算已轉(zhuǎn)化為一個不含截斷變量的多模態(tài)并聯(lián)系統(tǒng)問題。經(jīng)過這種轉(zhuǎn)化，構(gòu)建ARBIS方法計算新的矩獨(dú)立重要性測度指標(biāo)的算法如下：

步驟1自適應(yīng)搜索獲取超球半徑βopt，并計算無條件失效概率Pf；同時得到半徑為βopt的超球外失效點組成的樣本點矩陣B。

步驟2當(dāng)輸入變量Xi縮減到區(qū)間Ui(Q)時，在樣本點矩陣B中搜索輸入變量Xi落在縮減區(qū)間Ui(Q)的樣本點，記樣本點個數(shù)為m。

步驟3計算條件失效概率：

(14)

式中：N為采用ARBIS方法的初始樣本點個數(shù)。

4 算例分析

4.1 數(shù)值算例

某一結(jié)構(gòu)系統(tǒng)的功能函數(shù)可表示為

(15)

式中：X1、X2和X3為相互獨(dú)立的輸入隨機(jī)變量，均服從標(biāo)準(zhǔn)正態(tài)分布，即Xi～N(0,1),i=1,2,3。

表1 DLRS MC方法和ARBIS方法計算的新的矩獨(dú)立重要性測度指標(biāo)值Table 1 New moment-independent importance measure indices of numerical example computed by DLRS MC and ARBIS methods

圖2 DLRS MC方法和ARBIS方法計算的新的矩獨(dú)立重要性測度指標(biāo)的收斂曲線Fig.2 Convergence curves of new moment-independent importance measure indices of numerical example computed by DLRS MC and ARBIS methods

4.2 工程算例

在汽車結(jié)構(gòu)中，車橋承載著大部分汽車重量，其通過懸臂與車架相接，將來自車輪的牽引力和制動力，還有側(cè)向力經(jīng)過懸架傳遞給車架，起主要承載作用的就是汽車前軸[23]。由于工字型截斷梁能夠提高抗彎強(qiáng)度，因此前軸通常采用工字結(jié)構(gòu)梁。

圖3為前軸的結(jié)構(gòu)示意圖。危險截面處的最大正應(yīng)力為σ=M/Wx，最大切應(yīng)力為τ=T/Wρ，M為彎矩，T為扭矩，Wx為截面系數(shù)，Wρ為極截面系數(shù)，且有

(16)

Wρ=0.8bt2+0.4[a3(h-2t)/t]

(17)

考慮前軸結(jié)構(gòu)靜強(qiáng)度失效，有以下極限狀態(tài)方程：

(18)

式中：σS為靜強(qiáng)度屈服極限，根據(jù)前軸材料特性有σS=460 MPa;g為裕度值。

前軸結(jié)構(gòu)尺寸參數(shù)和承受外載看作獨(dú)立正態(tài)隨機(jī)變量，分布參數(shù)如表2所示。

采用DLRS MC方法和ARBIS方法計算的新的矩獨(dú)立重要性測度指標(biāo)結(jié)果如表3所示，同時也給出了計算功能函數(shù)的總次數(shù)?？梢钥闯觯?種算法計算結(jié)果相近。圖4給出了DLRS MC方法和ARBIS方法計算結(jié)果標(biāo)準(zhǔn)差的收斂曲線?？梢?，估計值的標(biāo)準(zhǔn)差都比較小，因此計算結(jié)果可靠穩(wěn)健。DLRS MC方法的總次數(shù)為446 000，而ARBIS方法計算總次數(shù)僅為4 475，因此ARBIS方法的計算效率很高。

2種算法計算的新的矩獨(dú)立重要性測度指標(biāo)的排序均相同：t>T>b>a>h>M。這表明當(dāng)分別固定6個輸入隨機(jī)變量在其各自的分布區(qū)域內(nèi)縮減變化時，t對汽車前軸結(jié)構(gòu)的失效概率影響程度最大，而M對其失效概率的影響幾乎可以忽略不計。

圖3 汽車前軸結(jié)構(gòu)示意圖Fig.3 Schematic of automobile front axle structure

表2 汽車前軸結(jié)構(gòu)各輸入變量分布參數(shù)Table 2 Distribution parameters of input variables of automobile front axle structure

表3 DLRS MC方法和ARBIS方法計算汽車前軸指標(biāo)結(jié)果Table 3 Results of indices of automobile front axle computed by DLRS MC and ARBIS methods

圖4 DLRS MC方法和ARBIS方法計算汽車前軸的收斂曲線Fig.4 Convergence curves of automobile front axle computed by DLRS MC and ARBIS methods

5 結(jié) 論

本文提出了一種新的基于失效概率的矩獨(dú)立重要性測度分析方法來計算輸入隨機(jī)變量對結(jié)構(gòu)系統(tǒng)輸出響應(yīng)失效概率的平均影響程度。

1) 建立了一種新的區(qū)間劃分技術(shù)來等可能地獲得輸入隨機(jī)變量所有可能的縮減區(qū)間，并給出了相應(yīng)證明。

2) 引入了自適應(yīng)超球重要抽樣方法來進(jìn)行新指標(biāo)的求解，提高了計算效率。

3) 給出了一個數(shù)值算例和一個工程算例，說明了本文所提新指標(biāo)的意義，同時也驗證了新算法的高效性。