基于聚類分析的核密度估計重要抽樣方法

摘"要：對于結(jié)構(gòu)體系失效概率的計算，當(dāng)其失效概率較小時，馬爾科夫蒙特卡洛方法難以獲取足夠的失效樣本點，且無法覆蓋全部失效區(qū)域，導(dǎo)致采用核密度估計重要抽樣方法計算失效概率精度較低.因此，提出一種基于聚類分析的核密度估計重要抽樣方法.該方法首先利用拉丁超立方抽樣得出第1次馬爾科夫蒙特卡洛方法抽樣的初始樣本點，并通過抽樣獲得失效樣本；然后將獲得的失效樣本進行聚類分析，獲得結(jié)構(gòu)失效模式的個數(shù)，并選取出具有代表性的失效樣本作為初始失效樣本點，進行第2次馬爾科夫蒙特卡洛方法抽樣;最后基于第2次抽樣獲得的失效樣本點，采用核密度估計的重要抽樣方法進行失效概率計算.算例分析表明，該方法針對小失效概率的體系可靠度問題具有良好的計算精度.

關(guān)鍵詞：結(jié)構(gòu)體系失效概率；馬爾科夫蒙特卡洛方法；核密度估計；聚類分析；重要抽樣方法

中圖分類號：TU311.2

文獻標(biāo)志碼：A

0"引"言

受施工設(shè)計等條件的影響，結(jié)構(gòu)參數(shù)及作用呈現(xiàn)隨機性的特征，所以結(jié)構(gòu)的失效也是隨機的，從不確定性的角度評估結(jié)構(gòu)的失效概率則十分重要.現(xiàn)代結(jié)構(gòu)體系復(fù)雜，往往有多個失效模式，當(dāng)結(jié)構(gòu)破壞時，結(jié)構(gòu)失效模式的出現(xiàn)是隨機的，因此復(fù)雜結(jié)構(gòu)的失效概率求解需要考慮結(jié)構(gòu)體系失效概率.蒙特卡洛模擬（Monte Carlo simulation，MCS）方法概念明確、方法簡單，且可以普遍適用于各種復(fù)雜的結(jié)構(gòu)體系失效概率的計算分析.但當(dāng)結(jié)構(gòu)失效概率較小時，運用MCS方法計算需較多樣本，單一樣本的獲取需要調(diào)用功能函數(shù)進行計算，所需的計算機時耗費巨大，效率十分低下.

為解決MCS方法計算量過大的問題，眾多學(xué)者提出了基于MCS方法改進的結(jié)構(gòu)體系失效概率的計算方法，例如，基于驗算點的重要抽樣方法^［1-2］、方向重要抽樣方法^［3-4］及基于核密度估計的重要抽樣方法^［5-6］等.但對于結(jié)構(gòu)體系失效概率問題，結(jié)構(gòu)的極限狀態(tài)方程往往存在多個驗算點，而基于驗算點的重要抽樣方法和方向重要抽樣方法都需要事先確定極限狀態(tài)方程的多個驗算點或者多個重要抽樣方向，在實際求解過程中，驗算點及重要抽樣方向均難以獲得，導(dǎo)致采用上述2種方法求解結(jié)構(gòu)體系失效概率比較困難.而基于核密度估計的重要抽樣方法不需要事先獲得驗算點，進一步降低了失效概率的求解難度.Ang等^［5］提出了核密度估計的重要抽樣方法，該方法首先需獲取一定量的失效區(qū)域內(nèi)的樣本，然后根據(jù)這些樣本采用核密度估計方法對失效域的概率密度函數(shù)進行估計，將估計得到的概率密度函數(shù)作為重要抽樣函數(shù)進行失效概率求解，但是當(dāng)失效概率極小時，獲取失效區(qū)域內(nèi)的點變得十分困難.Au等^［6］根據(jù)上述思路，采用馬爾科夫鏈的蒙特卡洛方法（MCMC）來獲取失效區(qū)域內(nèi)的點，極大程度提高了失效區(qū)域樣本點的獲取效率.對于結(jié)構(gòu)動力失效概率計算領(lǐng)域，也可將結(jié)構(gòu)動力失效概率問題轉(zhuǎn)換為體系失效概率進行求解，因此在地震易損性領(lǐng)域，基于密度估計的易損性分析也得到了較多的應(yīng)用^［7-9］.

對于結(jié)構(gòu)體系失效概率的計算，往往存在多個失效區(qū)域，當(dāng)失效概率較小時，MCMC也難以跨越邊界，獲取另一區(qū)域足夠多的失效樣本點，從而導(dǎo)致重要抽樣函數(shù)估計精度不夠，失效概率計算誤差較大.針對上述結(jié)構(gòu)體系失效概率中存在的問題，在小失效概率的計算問題上，本研究提出了一種基于聚類分析核密度估計的結(jié)構(gòu)體系失效概率的重要抽樣方法，克服了文獻［6］中的失效樣本點難以在各個失效區(qū)域中全面獲得的缺陷.該方法首先采用MCMC，獲取一定量的失效樣本點，然后將采用聚類分析獲得可能失效區(qū)域內(nèi)最靠近極限狀態(tài)方程的失效點，再采用MCMC對樣本點附近區(qū)域進行重點抽樣，從而獲得較多的失效樣本點，最后基于失效樣本點根據(jù)核密度估計重要抽樣方法構(gòu)造的核密度重要抽樣函數(shù)，計算結(jié)構(gòu)體系失效概率.

1"基于核密度估計的重要抽樣方法

核密度估計的重要抽樣方法的基本思路是通過失效區(qū)域樣本點構(gòu)造合理的重要抽樣函數(shù)，進而進行失效概率的計算^［5］，可不用獲得驗算點，降低了結(jié)構(gòu)失效概率的求解難度.

1.1"獲取失效域內(nèi)的樣本點

1）定義馬爾科夫鏈在失效域中的極限分布，可按式（1）進行構(gòu)造，

q_F（x）=I_f（x）f_x（x）P_f（1）

式中，q_F（x）表示馬爾科夫鏈在失效域中的極限分布，f_x（x）表示隨機變量的聯(lián)合概率密度函數(shù)，P_f為失效概率，I_f（x）為指示函數(shù).

2）選擇合理的建議分布，建議分布一般有2種形式可供選取，多維均勻分布與多維正態(tài)分布.本研究主要介紹多維均勻分布形式，其分布函數(shù)式為，

p^*（θ|x）=1∏nk=1l_k""ifx_k－θ_k≤l_k2，k=1，…，n0"其他 （2）

式中，x_k為抽取的樣本，l_k為樣本偏離上一個樣本的最大距離，θ_k為向量θ的第k個分量，n為變量空間的維數(shù).

3）選取馬爾科夫鏈?zhǔn)в虺跏紭颖?，一般可通過工程經(jīng)驗選取.

4）生成馬爾科夫鏈第j個狀態(tài)向量，馬爾科夫鏈的當(dāng)前狀態(tài)主要由前一個狀態(tài)通過建議分布和Metropolis-Hastings準(zhǔn)則確定，首先由建議分布和前一狀態(tài)的狀態(tài)向量計算得到當(dāng)前狀態(tài)向量，并將向量代入極限分布，則兩者比值為，

r=I_f（ε）f_x（ε）I_f（x_j－1）f_x（x_j－1）（3）

式中，ε為備選向量，根據(jù)Metropolis-Hastings準(zhǔn)則，可按照式（4）來判斷是否將ε作為下一狀態(tài)向量，

x_j=ε"min（1，r）gt;rand（0，1）x_j－1"min（1，r）≤rand（0，1） （4）

式中，rand（0，1）表示［0，1］范圍內(nèi)的均勻分布隨機變量.

5）重復(fù)步驟（4），生成馬爾科夫鏈樣本，直至失效樣本數(shù)達到指定的樣本數(shù)，則抽樣完成.

1.2"利用核密度估計構(gòu)造重要抽樣函數(shù)

前面尋找失效域樣本點的目的就是利用這些樣本點構(gòu)造重要抽樣函數(shù)，并結(jié)合重要抽樣函數(shù)進行可靠度計算.文獻［6］中采用了自適應(yīng)寬核密度估計的方法來構(gòu)造重要抽樣函數(shù)，是一種非參數(shù)概率密度估計方法，其基本原理是在固定寬核密度函數(shù)估計的基礎(chǔ)上，通過修正窗口寬度參數(shù)來實現(xiàn)的，其具體形式為，

k（x）=1M∑Mj=11（wλ_j）ⁿKx－x^jwλ_j（5）

式中，k（·）為核概率密度函數(shù)，其形式一般可以按照文獻［6］進行選取；λ_j為局部寬度因子;w和λ_j的求解公式見式（6）和式（7），

λ_j=∏Mk=1f（x^k）^1Mf（x^j）^α（6）

w=M^－1n+4_d（7）

式中，α為靈敏因子，0≤α≤1;M_d為不同樣本的個數(shù)，M_d≤M.

1.3"基于重要抽樣函數(shù)進行樣本抽樣

由于核函數(shù)為M項累加的形式，需采用如下過程進行隨機變量抽樣：首先在［1，M］區(qū)間上生成1個均勻的整數(shù)隨機變量u，如果u=j，則可選擇第j個分量的核密度抽樣函數(shù)產(chǎn)生隨機變量樣本，第j分量核密度抽樣函數(shù)為，

k_j（x）=1WⁿKx－x^jW（8）

式中，當(dāng)采用自適應(yīng)寬核密度函數(shù)估計時，W可取wλ_j，λ_j和w分別按照式（6）和式（7）計算.

1.4"求解失效概率

利用重要抽樣概率密度函數(shù)產(chǎn)生一定量的樣本后，就可以采用重要抽樣方法計算失效概率.

2"基于聚類分析的改進核密度估計的重要抽樣方法

基于核密度估計的重要抽樣方法的關(guān)鍵在于采用MCMC獲取失效域的樣本，樣本沿極限狀態(tài)方程分布越均勻，樣本越多，則采用核密度估計的重要抽樣方法計算結(jié)果越精確，而失效點選取質(zhì)量的好壞主要取決于l_k，l_k越大，則失效點覆蓋的失效區(qū)域越大，但樣本點比較少；l_k越小，失效點覆蓋的失效區(qū)域越小，但樣本點比較多.

仔細(xì)考察式（4）所提出的Metropolis-Hastings接受與拒絕準(zhǔn)則，可以發(fā)現(xiàn)，從概率角度上面來看，r大于rand（0，1）的概率為，

p{rgt;rand（0，1）}=∫^r₀11－0dx=r（9）

通過式（9）可以看出，r越大，則新生成樣本越容易被接受，而r值是根據(jù)式（3）進行計算的，從上述分析中可以確定，當(dāng)新生成樣本的聯(lián)合概率密度值較大時，相應(yīng)的r值也會較大，同時也增加了新樣本被接受的機會.換言之，新樣本會向聯(lián)合概率密度值較大的區(qū)域聚集，當(dāng)l_k越小時，這個特點越明顯.如圖1所示，當(dāng)失效概率較小時，要獲取A點附近區(qū)域樣本點的概率極小，導(dǎo)致MCMC獲得的樣本一直在區(qū)域1，則此時只能增加l_k或者增大MCMC樣本點個數(shù)才能越過A點附近區(qū)域，進而獲取區(qū)域2的樣本點.當(dāng)l_k較小時，需要抽取巨大數(shù)量的樣本點才能達到期望的計算精度，但l_k取較大值時，同樣也會導(dǎo)致計算精度低的問題.

為了失效樣本點的抽取達到比較理想的效果，本研究提出采用2次MCMC尋找核密度估計的重要抽樣方法所需要的失效點.所采取策略的基本思想是第1次MCMC取用較大的l_k，以便第1次的失效樣本點能夠盡可能覆蓋整個極限狀態(tài)方程所有的失效區(qū)域；引入聚類分析獲取各失效區(qū)域的較優(yōu)樣本點，并采用較小的l_k進行MCMC抽樣，以便獲得的失效樣本點能盡可能多，通過2次MCMC抽樣，使得獲取的失效域樣本能覆蓋所有失效區(qū)域，并且具有一定的數(shù)量；最后采用第2次獲得的失效樣本點，基于核密度估計的重要抽樣方法計算失效概率.具體過程如下：

1）確定第1次MCMC初始失效樣本.Au等^［6］提出初始失效樣本點可以通過數(shù)值方法或者工程經(jīng)驗來確定，通過工程經(jīng)驗的方法確定失效點往往距離極限狀態(tài)方程比較遠(yuǎn)，本研究提出采用拉丁超立方抽樣（LHS）方法來確定初始失效點，一般LHS方法將隨機變量區(qū)域按照概率等分，由于離均值點較遠(yuǎn)區(qū)域樣本點的概率較小，因此采用這種方法獲得離均值較遠(yuǎn)的樣本較少，但一般離均值較遠(yuǎn)的樣本更加可能為失效樣本.基于以上原因，本研究將隨機變量區(qū)域采用等距離進行分割，使得離均值較遠(yuǎn)的樣本點增多，從而達到獲得失效樣本數(shù)的概率增加的目的.

2）第1次MCMC抽樣獲取失效樣本.第1次MCMC抽樣可盡量使失效點樣本能夠覆蓋更多的失效區(qū)域.因此，通過增大l_k或者增大樣本點抽樣次數(shù)的方法均可使得備選樣本點達到更廣的區(qū)域，而增大樣本點勢必會增加計算成本，所以本研究建議采用增大l_k的方法來獲取分布更加廣泛的失效樣本點，第1次MCMC抽樣l_k可取8倍的隨機變量標(biāo)準(zhǔn)差.

3）確定第2次MCMC初始失效樣本.根據(jù)前面對Metropolis-Hastings準(zhǔn)則的分析表明，MCMC失效點一般分布在隨機變量聯(lián)合概率密度較大的區(qū)域，樣本點表現(xiàn)出高度區(qū)域聚集性，即圖1所示的區(qū)域1和區(qū)域2.根據(jù)上述特點，可采用聚類分析獲取第1次MCMC獲取的失效域樣本的聚類數(shù)與對應(yīng)的聚類樣本點，在每個聚類中選取最靠近極限狀態(tài)方程的失效點作為第2次MCMC抽樣的初始失效樣本，以保證MCMC抽樣能夠獲取各區(qū)域內(nèi)的失效樣本點.

本研究采用的K-means聚類分析方法進行聚類分析^［10］，其基本步驟如下：

a）確定聚類數(shù)k及需要進行聚類分析的數(shù)據(jù).

b）初始化k個聚類中心，首先在所有數(shù)據(jù)中隨機選取1個點作為第1個聚類中心，在數(shù)據(jù)中選取離第1個聚類中心最遠(yuǎn)的點作為第2個聚類中心，計算前面2個聚類中心的中心點，選取數(shù)據(jù)中離這個中心點最遠(yuǎn)的點作為第3個聚類中心，依次選出k個初始聚類中心X^k_c.

c）對于數(shù)據(jù)集中的某一數(shù)據(jù)X_i，根據(jù)d^k_i=X_i－X^k_c計算該數(shù)據(jù)與每個聚類中心的距離，并將該數(shù)據(jù)歸于距離最小的聚類中心X^k_c所在的類.

d）根據(jù)第k個聚類，重新計算新聚類中心，聚類中心的計算公式為，

X^k_c=1n∑ni=1X^k_i（10）

e）計算測度函數(shù)E為，

E=∑ki=1∑nj=1Xⁱ_j－Xⁱ_c²（11）

f）根據(jù)E_i－E_i－1lt;ε計算2次測度函數(shù)之間的誤差是否小于預(yù)設(shè)值，如果滿足則迭代完畢，如果不滿足，則進入步驟c）重新迭代計算.

最佳聚類數(shù)的基本計算步驟如下^［10-11］：

a）確定聚類數(shù)的搜索范圍為k∈［k_min，k_max］.

b）對于每一聚類數(shù)，采用K-means聚類分析方法進行聚類分析.

c）計算聚類數(shù)對應(yīng)的評價指標(biāo)為，

I=1k∑kj=11n∑ni=1X_i－X^j_c²min1≤i≤j≤kXⁱ_c－X^j_c²（12）

d）根據(jù)步驟c）中計算指標(biāo)I的最小值確定最佳聚類數(shù).

4）第2次MCMC抽樣獲取失效樣本.設(shè)聚類數(shù)為n，則通過第3）步獲得n個失效樣本，將獲取的n個失效樣本點作為初始失效樣本，分別進行n次MCMC抽樣，第2次MCMC抽樣時l_k可按照較小值進行計算，使得每個區(qū)域內(nèi)獲得的失效點較多.設(shè)每次抽樣獲得失效域樣本數(shù)為m，則總共可獲得n×m個失效域樣本.

5）計算失效概率.將得到的失效域樣本采用核密度估計的重要抽樣方法進行失效概率計算.

3"算例

設(shè)有1串聯(lián)體系，其極限狀態(tài)方程為，

G（x，y）=x+y－a=0x－y－a=0－x－y－a=0y－x－a=0 （13）

式中，x和y為獨立的標(biāo)準(zhǔn)正態(tài)分布隨機變量，其均值為0，標(biāo)準(zhǔn)差σ為1，a為變量.根據(jù)概率論的基本原理，a可取標(biāo)準(zhǔn)正態(tài)分布隨機變量標(biāo)準(zhǔn)差σ的整數(shù)倍.當(dāng)a取3時（即a=3σ），失效區(qū)域樣本點獲取概率為2.699E-3;當(dāng)a取4時（即a=4σ），失效區(qū)域樣本點獲取概率為6334E-5;當(dāng)a取6時（即a=6σ），失效區(qū)域樣本點獲取概率為1.973E-9.因此，a取值6σ時，從檢驗本文方法在小失效概率的計算有效性上，可完全滿足要求.綜上，為對比本文中改進的方法和文獻［6］中的方法，將a分別取值為3、4和6，然后進行失效域樣本抽樣和失效概率的計算.

圖2～圖4給出了a取值不同時的MCMC抽樣獲得的失效域樣本點.從圖中可以看出，當(dāng)a取值較小時，本文方法和文獻［6］的方法均能獲得較好的效果，但隨著a值的增大，本文方法計算的失效樣本在整個失效區(qū)域仍然分布比較均勻，而文獻［6］的方法得出的失效樣本只能集中在部分區(qū)域，出現(xiàn)這種情況的原因主要在于Metropolis-Hastings準(zhǔn)則的

限定，根據(jù)該準(zhǔn)則可知，失效區(qū)域概率密度較小時備選樣本很難被接受，而對于a值取值比較大的情況下，4個角點的聯(lián)合概率密度均較小，由于文獻［6］中的l_k較小，難以跨過角點，導(dǎo)致失效樣本過分集中在某個區(qū)域，而本文提出的方法因為引入了聚類分析，采取了2次MCMC抽樣的方法，使得l_k的取值可根據(jù)抽樣的需要進行調(diào)整，因此在獲取失效樣本時則不存失效樣本不均勻的問題.

1給出了失效概率計算結(jié)果，從表中可以看出，隨著a值的增大，失效概率隨之減小，出現(xiàn)這種情況的主要原因是當(dāng)a值越大時，根據(jù)正態(tài)分布累計概率計算的基本原理，落在失效區(qū)域外（圖4正方形區(qū)域外）的樣本點急劇減少，從隨機抽樣的角度看，式（13）所示的串聯(lián)體系的體系失效概率則會變小.此外，本文方法在獲取失效點基礎(chǔ)上，使用核密度估計的重要抽樣方法獲得了較為精確的結(jié)果，但是文獻［6］中的方法由于失效樣本無法在整個失效域上分布均勻，導(dǎo)致當(dāng)a=4和6時，失效概率計算結(jié)果過小，產(chǎn)生了較大誤差，僅當(dāng)a=3時能獲得較好的結(jié)果.

4"結(jié)"論

本研究在核密度估計重要抽樣方法的基礎(chǔ)上，引入聚類分析，采用2次MCMC抽樣的策略，提出一種基于聚類分析的核密度估計重要抽樣方法，得到了以下結(jié)論：

1）通過2次MCMC抽樣的策略，使l_k可采用不同的取值，首次MCMC抽樣時l_k取較大值，令失效樣本可分布較廣的區(qū)域，但樣本數(shù)量少；基于第1次獲得失效樣本，采用聚類分析獲得最優(yōu)MCMC失效樣本點，并采用較小的l_k進行MCMC抽樣.該策略可有效解決小失效概率時的失效樣本抽樣效率低的問題.

2）基于本研究提出的失效樣本點的抽樣方法，提出了改進的核密度估計的重要抽樣方法，通過算例研究表明，本研究方法能夠有效解決小失效概率下的系統(tǒng)可靠度問題，為復(fù)雜系統(tǒng)可靠度問題提供了一種新思路.

參考文獻：

［1］Bucher C G.Adaptive sampling-an iterative fast Monte Carlo procedure ［J］.Struct Saf，1988，5（2）：119-126.

［2］Mori Y，Kato T.Multi-normal integrals by importance sampling for series system reliability ［J］.Struct Saf，2003，25（4）：363-378.

［3］Ditlevsen O，Olesen R，Mohr G.Solution of a class of load combination problems by directional simulation［J］.Struct Saf，1986，4（2）：95-109

［4］Jinsuo N，Ellingwood B R.Directional methods for structural reliability analysis［J］.Struct Saf，2000，22（3）：233-249.

［5］Ang G L，Ang A S，Tang W H.Optimal importance sampling density estimator［J］.J Eng Mech，1992，118（6）：1146-63.

［6］Au S K，Beck J L.A new adaptive importance sampling scheme for reliability Calculations［J］.Struct Saf，1999，（21）：135-158.

［7］江輝，李辰，馮夢瑤，等.基于核密度估計的干接縫裝配式橋墩概率性地震損傷特性分析［J］.東南大學(xué)學(xué)報（自然科學(xué)版），2021，51（4）：566-574.

［8］董俊，曾永平，單德山.核密度估計的橋梁結(jié)構(gòu)地震易損性分析［J］.哈爾濱工業(yè)大學(xué)學(xué)報，2018，50（3）：109-117.

［9］單德山，張二華，董俊，等.基于核密度估計的鐵路橋梁構(gòu)件地震易損性分析［J］.鐵道學(xué)報，2019，41（8）：108-116.

［10］翟東海，魚江，高飛，等.最大距離法選取初始簇中心的K-means文本聚類算法的研究［J］.計算機應(yīng)用研究，2014，31（3）：713-719.

［11］周世兵，徐振源，唐旭清.K-means算法最佳聚類數(shù)確定方法［J］.計算機應(yīng)用，2010，30（8）：1995-1998.

（實習(xí)編輯：羅"媛）

An Importance Sampling Method for Kernel Density Estimation Based on Cluster Analysis

SHAO Junhu¹，F(xiàn)ENG Miao¹，CHEN Xiaoping¹，ZHOU Yang¹，SONG Shuai²，WANG Fangfang¹

（1.School of Architecture and Civil Engineering，Chengdu University，Chengdu 610106，China;

2.College of Engineering and Technology，Chengdu University of Technology，Leshan 614000，China）

Abstract：

For the calculation of the failure probability of structural system，when the failure probability is low，it is difficult to obtain sufficient failure sample points by using Markov Monte Carlo method，and covering all failure areas is also impossible.Therefore，this leads to a low accuracy in calculating failure probability by using importance sampling method for kernel density estimation.The proposed method uses Latin hypercube sampling to obtain the initial sample points，which is the same as the one obtained by using Markov Monte Carlo method and the failure samples are also obtained through sampling.Then，the number of structural failure modes are obtained by the cluster analysis of the failure samples，and representative failure samples are selected as the initial failure sample points for the second Markov Monte Carlo method sampling.Finally，based on the failure sample points obtained from the second sampling，the importance sampling method for kernel density estimation is used to calculate the failure probability.The analysis of numerical examples shows that the proposed method has good accuracy for the system reliability with small failure probability.

Key words：

failure probability of structural systems;Markov Monte Carlo method;kernel density estimation;cluster analysis;an importance sampling method