基于Kriging模型的重力壩非概率可靠性分析

(南昌大學(xué) 建筑工程學(xué)院，南昌 330031)

1 研究背景

重力壩特殊的結(jié)構(gòu)性質(zhì)及其復(fù)雜的工作環(huán)境，決定了對其進(jìn)行可靠性分析的復(fù)雜性。在我國，首次把結(jié)構(gòu)可靠性分析方法引入大壩可靠度計算的是河海大學(xué)“水工可靠度”科研組。武清璽等[1]提出了基于有限元法的重力壩可靠性分析，并進(jìn)一步計算了重力壩抗滑穩(wěn)定的可靠性。江勝華等[2]運(yùn)用加權(quán)響應(yīng)面法，計算了重力壩各斷層抗滑穩(wěn)定的可靠指標(biāo)，并與安全系數(shù)進(jìn)行了比較，表明其安全系數(shù)與可靠指標(biāo)并非存在完全的對等關(guān)系。概率可靠性模型雖能較好地描述參數(shù)的不確定性，但需不確定參數(shù)完整的概率分布信息，這對于絕大多數(shù)大壩工程通常難以獲得，其不確定參數(shù)相關(guān)信息通常均較缺乏；而且有些不確定參數(shù)(如受人為控制的庫水位)并非嚴(yán)格意義上的隨機(jī)變化。另一方面，已有研究表明，概率可靠性計算結(jié)果對概率模型參數(shù)極其敏感，而重力壩工程失事為小概率事件，不確定參數(shù)概率函數(shù)尾部的較小變化或計算方法本身誤差，都可能導(dǎo)致計算結(jié)果出現(xiàn)這些問題，導(dǎo)致概率可靠性理論應(yīng)用于實(shí)際工程受到限制[3-5]。

在大壩工程中，雖難以獲取不確定參數(shù)的準(zhǔn)確統(tǒng)計信息，但不確定參數(shù)的變化幅度較易確定[6-7]。近些年發(fā)展起來的N-PR分析方法，采用區(qū)間或橢球等凸集合模型[8](Convex model)描述這類不確定但有界的參數(shù)，僅要求知道不確定參數(shù)集合的界限，無需了解其概率分布情況，對于不確定參數(shù)信息量少、計算精度要求高的重力壩可靠性問題研究有著獨(dú)特優(yōu)勢和較強(qiáng)的工程適用性。然而重力壩因其結(jié)構(gòu)受荷載、環(huán)境、材料等多種因素協(xié)同影響，其極限狀態(tài)方程常表現(xiàn)為高度非線性且較難用顯式的數(shù)學(xué)方程給予刻畫，導(dǎo)致需借助數(shù)值模擬的手段將其顯式化表達(dá)，在結(jié)構(gòu)工程中，主要有響應(yīng)面法[9](RSM)、神經(jīng)網(wǎng)絡(luò)法、支持向量機(jī)法和Kriging模型[10]等方法。工程中常用RSM以二項(xiàng)式模型來近似擬合其功能函數(shù)，但在高度非線性問題中，RSM需大量的計算工作，且對于高度非線性可靠性問題常出現(xiàn)無法收斂和精度不足的情況[11]。Gaspar等[12]指出，與多項(xiàng)式回歸模型相比，Kriging模型求解的失效概率更為準(zhǔn)確。因Kriging模型較其他方法具有更優(yōu)的預(yù)測能力及靈活性[13]，近10 a，Kriging模型被廣泛應(yīng)用于結(jié)構(gòu)可靠性分析。文獻(xiàn)[14]利用Kriging模型實(shí)現(xiàn)了對某航空結(jié)構(gòu)極限狀態(tài)方程的精確擬合，且其可靠性計算效率和精度均明顯優(yōu)于RSM計算結(jié)果。

據(jù)此，考慮到重力壩不確定參數(shù)統(tǒng)計信息較難獲取，而其區(qū)間界限信息較易確定現(xiàn)象，本文基于凸集模型的區(qū)間變量來描述重力壩強(qiáng)度和穩(wěn)定N-PR。在此基礎(chǔ)上，借鑒Kriging模型強(qiáng)預(yù)測及擬合能力，并結(jié)合有限元ABAQUS和MatLab平臺優(yōu)勢，研發(fā)相應(yīng)實(shí)施程序，研究建立重力壩壩體單元強(qiáng)度破壞和沿壩基面抗滑失穩(wěn)的N-PR計算模型。

2 重力壩N-PR指標(biāo)計算模型

2.1 基于區(qū)間變量的結(jié)構(gòu)N-PR指標(biāo)

M=g(x)=g(x1,x2,…,xn) 。

(1)

式中：M為區(qū)間變量；g(x)為xi的連續(xù)函數(shù)。

令區(qū)間變量M的均值和離差分別為Mc和Mr，則結(jié)構(gòu)的N-PR指標(biāo)η為

η=Mc/Mr。

(2)

圖1 線性和非線性功能函數(shù)Fig.1 Linear and nonlinear functions

對xi進(jìn)行如下標(biāo)準(zhǔn)化變換，即

(3)

M=g(x1,x2,…,xn)=G(δ1,δ2,…,δn)=0 。(4)

對于常見功能函數(shù)為多區(qū)間變量的線性函數(shù)情況為

(5)

根據(jù)式(3)對ri和sj進(jìn)行如下變換

(6)

則式(5)功能方程可標(biāo)準(zhǔn)化為

(7)

則結(jié)構(gòu)的N-PR指標(biāo)η為

(8)

2.2 重力壩單元N-PR指標(biāo)η計算模型

以2.1節(jié)所述結(jié)構(gòu)N-PR指標(biāo)η計算原理為基礎(chǔ)，結(jié)合重力壩區(qū)間參數(shù)選取和功能函數(shù)選擇，構(gòu)建重力壩單元強(qiáng)度破壞和沿壩基面抗滑失穩(wěn)的N-PR指標(biāo)計算模型。

2.2.1 重力壩單元強(qiáng)度N-PR指標(biāo)η計算模型

考慮重力壩單元強(qiáng)度破壞，由于重力壩結(jié)構(gòu)受力狀態(tài)復(fù)雜，根據(jù)單元的第一、第二和第三主應(yīng)力，采用愈茂宏雙剪強(qiáng)度理論[17]，三維狀態(tài)下單元強(qiáng)度破壞的響應(yīng)值函數(shù)可表示為

(9)

式中：g(x)為重力壩單元強(qiáng)度的功能函數(shù)(響應(yīng)值)；σ1,σ2和σ3分別為重力壩單元的第一、第二和第三主應(yīng)力；α=ft/fc,ft和fc為重力壩單元的抗拉強(qiáng)度和抗壓強(qiáng)度。

2.2.2 重力壩抗滑穩(wěn)定性N-PR指標(biāo)計算模型

考慮重力壩抗滑穩(wěn)定性，其失穩(wěn)破壞的響應(yīng)值函數(shù)可表示為

(10)

式中：k為滑移面上單元總數(shù)；f和c為分別為材料的摩擦系數(shù)和黏聚力；σyi和τxyi分別為單元i的正應(yīng)力和剪應(yīng)力；bi為單元i沿滑動面的邊長。

3 基于Kriging模型的重力壩N-PR指標(biāo)η計算

Kriging模型求解重力壩N-PR指標(biāo)η，其核心是通過有限的數(shù)值仿真計算成果，擬合得到一個表達(dá)式來近似代替復(fù)雜的功能函數(shù)M=g(x)。

3.1 Kriging模型的基本原理

Kriging模型是一種統(tǒng)計近似模型，對于任意輸入變量x∈Rn與響應(yīng)值y有如下關(guān)系，即

式中：f(x)T為回歸模型，f(x)T=[f1(x),f2(x),…,fp(x)]，fi(x)在MatLab軟件DACE工具箱中有零階(常數(shù))、一階和二階模型可供選擇，通常采用二次多項(xiàng)式函數(shù)形式；β為回歸模型系數(shù)，β=[β1,β2,…,βp]T；z(x)為隨機(jī)誤差，其均值為0、方差為σ2，則z(w)和z(x)協(xié)方差為

cov[z(w),z(x)]=σ2R(θ;w,x) 。

(12)

式中：R(θ;w,x)是以θ為參數(shù)的相關(guān)模型，其中w，x∈Rn。

相關(guān)模型選用Guass函數(shù)形式，則

假定隨機(jī)抽樣m個計算樣本點(diǎn)組成的矩陣為S=[s1,s2,…,sm]T，si∈Rn(i=1,2,…,m)，其每個樣本點(diǎn)的組成的響應(yīng)值矩陣為Y=[y1,y2,…,ym]T。定義R=(Rij)m×m為S中點(diǎn)坐標(biāo)的相關(guān)矩陣，其中，Rij=R(θ;si,sj)，(i,j=1,2,…,m)。

構(gòu)造m維權(quán)系數(shù)向量c=[c1,c2,…,cm]T，根據(jù)響應(yīng)值yi(i=1,2,…,m)的線性加權(quán)疊加來預(yù)測待測點(diǎn)x的響應(yīng)值，可得

(14)

(15)

根據(jù)以上內(nèi)容，最優(yōu)Kriging模型的構(gòu)建轉(zhuǎn)化為尋找最優(yōu)解θ*的優(yōu)化問題。

3.2 Kriging模型的檢驗(yàn)

Kriging模型作為一近似擬合模型，必然與真實(shí)模型存在一定誤差，為檢驗(yàn)Kriging模型擬合效果，特用復(fù)相關(guān)系數(shù)給以判斷，可用下式計算，即：

(16)

(17)

可知R2介于0到1之間，R2越趨于1表明擬合效果越佳，其檢驗(yàn)步驟如下：

(1)根據(jù)各參數(shù)區(qū)間范圍，采用拉丁超立方體獲取足夠樣本點(diǎn)，得到樣本點(diǎn)集U，再從中隨機(jī)選取樣本點(diǎn)組成子樣本集C1，C1∈U，并根據(jù)樣本子集C1建立Kriging模型。

(3)根據(jù)復(fù)相關(guān)系數(shù)R2判斷所建Kriging模型是否滿足預(yù)期要求，滿足則選用該模型，否則返回步驟1重新選擇樣本子集C1，重復(fù)以上步驟，直至滿足精度要求。

3.3 重力壩N-PR指標(biāo)η的計算

在運(yùn)用Kriging模型進(jìn)行N-PR指標(biāo)η計算時，需將擬合好的近似函數(shù)M=y(x)代替功能函數(shù)M=g(x)，并運(yùn)用數(shù)學(xué)規(guī)劃方法求解M=y(x)的最大值和最小值，因此可將N-PR指標(biāo)η的求解視為二次規(guī)化問題的極值求解。

二次規(guī)劃問題的數(shù)學(xué)模型可表述為

(18)

式中：x為重力壩參數(shù)變量集合,x={x1,x2,…,xn}。

由式(18)可得到功能函數(shù)的最小值Mmin，然后將目標(biāo)功能函數(shù)M=y(x)改為M=-y(x)便可得到最大值Mmax，進(jìn)而可求得N-PR指標(biāo)η。

基于Kriging模型的重力壩N-PR指標(biāo)η的計算步驟如下：

(1)根據(jù)重力壩基本資料確定不確定參數(shù)界限，抽樣1 000個樣本點(diǎn)形成總樣本點(diǎn)集U。

(2)隨機(jī)選取50個樣本點(diǎn)形成樣本子集C1，利用有限元ABAQUS建立重力壩模型進(jìn)行結(jié)構(gòu)分析，獲得各樣本點(diǎn)處的結(jié)構(gòu)響應(yīng)。

(3)依據(jù)應(yīng)力計算成果，根據(jù)式(9)計算重力壩壩體單元強(qiáng)度功能函數(shù)值，根據(jù)式(10)計算壩基面抗滑穩(wěn)定功能函數(shù)值。

(4)根據(jù)樣本點(diǎn)和功能函數(shù)值，利用MatLab調(diào)用Kriging模型對其進(jìn)行擬合得到近似功能函數(shù)M=y(x)，并根據(jù)3.2節(jié)所述檢驗(yàn)方法對其進(jìn)行判定，滿足要求則采用此M=y(x)作為近似功能函數(shù)，否則返回步驟2，直至得到符合精度要求的Kriging模型。

(5)基于擬合所得近似功能函數(shù)M=y(x)，根據(jù)式(18)執(zhí)行二次規(guī)劃求解，求得功能函數(shù)的最小值Mmin和最大值Mmax。

(6)根據(jù)功能函數(shù)的最小值和最大值，計算其均值Mc和離差Mr，進(jìn)而運(yùn)用式(8)計算壩體單元強(qiáng)度和沿壩基面抗滑穩(wěn)定N-PR指標(biāo)η。

4 工程應(yīng)用

某碾壓混凝土重力壩位于福建省永定縣境內(nèi)，本工程屬Ⅰ等樞紐工程。從左岸至右岸分別為1#—6#壩段，本文選用2#非溢流壩段作為研究對象，運(yùn)用本文方法分析其單元強(qiáng)度及沿壩基面抗滑穩(wěn)定可靠性。2#壩段全長50 m，壩頂高程179.0 m，最大壩高96 m，最小壩高56 m，大壩上游正常蓄水位173.0 m。

4.1 重力壩的有限元模型

為保證計算結(jié)果的準(zhǔn)確性，從壩踵向上取1.5倍最大壩高，壩址向下游取1.5倍最小壩高，壩基深度取83 m。本文結(jié)合有限元ABAQUS建模獲得各壩體單元應(yīng)力，其2#壩段三維有限元模型如圖2所示。其中，該模型總節(jié)點(diǎn)數(shù)28 851個、總單元數(shù) 24 711個。其中取垂直水流方向、順?biāo)鞣较蚝豌U直方向分別為建模坐標(biāo)x,y和z方向。

圖2 2#壩段網(wǎng)格劃分Fig.2 Meshing of dam segment 2#

4.2 材料參數(shù)和荷載

考慮壩體、壩基和混凝土的物理與力學(xué)參數(shù)變異性對重力壩強(qiáng)度和穩(wěn)定計算結(jié)果的影響，影響較小的參數(shù)視為定值考慮，如壩體上游畜水位因受人為控制可視為定值處理，及壩體和壩基泊松比等參數(shù)對其強(qiáng)度和穩(wěn)定影響較??；而影響較大的參數(shù)視為區(qū)間變量，如壩體和壩基彈性模量對其強(qiáng)度影響較為明顯，抗剪斷參數(shù)f和c對壩體穩(wěn)定影響較大。根據(jù)水庫管理局和工程局提供參數(shù)指標(biāo)，其定值參數(shù)統(tǒng)計見表1，不確定參數(shù)界限見表2。材料力學(xué)模型采用Drucker-Prager準(zhǔn)則。計算施加的荷載組合為：自重+上游水壓力+揚(yáng)壓力。

表1 定值參數(shù)統(tǒng)計Table 1 Statistics of deterministic parameters

表2 不確定參數(shù)界限Table 2 Ranges of uncertain parameters

注：上述材料參數(shù)界限均根據(jù)實(shí)測值確定

4.3 重力壩N-PR指標(biāo)η計算結(jié)果

根據(jù)3.3節(jié)計算步驟求解該重力壩壩段N-PR指標(biāo)η，其壩體單元N-PR指標(biāo)η分布如圖3(a)所示，巖基單元N-PR指標(biāo)η分布如圖3(b)所示，滑移面單元的N-PR指標(biāo)η分布及統(tǒng)計如圖4所示。

圖3 壩體單元和巖基單元的N-PR指標(biāo)η分布Fig.3 Distribution of N-PR index η of dam body androck foundation elements

圖4 滑移面單元N-PR指標(biāo)η分布及統(tǒng)計Fig.4 Distribution and statistics of N-PR index η ofsliding surface element

4.4 結(jié)果分析

圖3(a)給出了壩體單元強(qiáng)度N-PR指標(biāo)η分布，主要在壩踵附近出現(xiàn)少量η<1的單元，大量單元η>1，處于安全狀態(tài)。其中η最小值為0.56，位于7 482號單元。圖3(b)給出了巖基單元N-PR指標(biāo)η分布，在壩基面靠近壩踵處出現(xiàn)3個單元η<1。由于不考慮下游水壓力作用，下游壩基單元η較上游壩基單元η大。圖4給出了滑移面單元統(tǒng)計信息，超過99%巖基單元處于安全狀態(tài)。由圖3(a)可知，重力壩的可能失效單元主要分布在壩踵附近，故在重力壩建設(shè)過程中需針對壩踵部位實(shí)施有效的加固措施。根據(jù)式(10)得到抗滑穩(wěn)定N-PR指標(biāo)η為5.28，表明壩體整體處于穩(wěn)定狀態(tài)。

5 結(jié) 論

(1)考慮到重力壩失事概率小且不易獲取到不確定參數(shù)概率分布信息等的工程特點(diǎn)，構(gòu)建了一種基于不確定參數(shù)界限的重力壩單元N-PR計算模型，有效規(guī)避了經(jīng)典概率可靠性分析中要求不確定參數(shù)隨機(jī)變化且計算結(jié)果對參數(shù)高度敏感的局限。

(2)本文研究了一種Kriging模型通過合理選取和擬合功能函數(shù)以代替重力壩功能函數(shù)，較好適應(yīng)了重力壩功能函數(shù)高度非線性甚至難以顯式表達(dá)的特點(diǎn)，可實(shí)現(xiàn)重力壩N-PR指標(biāo)η的高效計算。

(3)以凸集來描述影響大壩安全的不確定因素，以一種非概率的安全指標(biāo)來度量大壩的安全程度，工程適用性和實(shí)用性較強(qiáng)，可深入分析重力壩安全性能。