馮曉雷， 黃 銘

(合肥工業(yè)大學(xué) 土木與水利工程學(xué)院，安徽 合肥 230009)

0 引 言

邊坡作為引水工程的重要組成部分，一旦出現(xiàn)滑坡對(duì)于引水工程的安全運(yùn)行有著重大的影響，且后期的維修耗時(shí)長(zhǎng)、成本高、難度大，因此需要對(duì)于重要邊坡進(jìn)行安全監(jiān)測(cè)、故障診斷和剩余使用壽命預(yù)測(cè)，以確保引水工程的安全運(yùn)行[1]。

邊坡的監(jiān)測(cè)位移是表征邊坡安全的重要信息，通過(guò)分析邊坡的位移監(jiān)測(cè)資料，來(lái)研究滑坡的預(yù)警方法，是目前的一個(gè)重點(diǎn)研究?jī)?nèi)容。徐衛(wèi)亞等[2]提出了以變形變化趨勢(shì)為邊坡破壞標(biāo)準(zhǔn)，以改進(jìn)的灰色-時(shí)序分析時(shí)變預(yù)測(cè)模型進(jìn)行預(yù)警預(yù)測(cè)；許強(qiáng)等[3]通過(guò)大量滑坡變形監(jiān)測(cè)數(shù)據(jù)，提出了加速度的滑坡預(yù)警方法與臨滑預(yù)警指標(biāo)。但是目前研究成果多集中于位移和位移速率及其隨時(shí)間變化規(guī)律方面，很難解釋引起邊坡蠕滑變形與失穩(wěn)的損傷機(jī)制與動(dòng)因，與邊坡穩(wěn)定系數(shù)相關(guān)的失穩(wěn)破壞位移判據(jù)也研究甚少[4]。

關(guān)于剩余壽命的預(yù)測(cè)方法主要包括基于失效機(jī)制模型、數(shù)據(jù)驅(qū)動(dòng)以及融合上述兩類(lèi)思想的預(yù)測(cè)方法[5]。其中數(shù)據(jù)驅(qū)動(dòng)方法的隨機(jī)過(guò)程常用來(lái)描述性能變化呈現(xiàn)出隨機(jī)特性的退化過(guò)程[6]，隨機(jī)過(guò)程模型中的維納過(guò)程由于具有優(yōu)良的統(tǒng)計(jì)特性是目前應(yīng)用較為廣泛的一種剩余壽命預(yù)測(cè)方法, 不同領(lǐng)域的學(xué)者利用單一性能指標(biāo)對(duì)目標(biāo)進(jìn)行了維納過(guò)程建模，獲得了較好的剩余壽命預(yù)測(cè)結(jié)果[7-9]。因此本文利用隨機(jī)過(guò)程中維納過(guò)程的優(yōu)點(diǎn)對(duì)于邊坡的退化過(guò)程進(jìn)行相關(guān)分析。

基于上述考慮，本文提出一種基于損傷力學(xué)的維納過(guò)程隨機(jī)變量模型對(duì)于邊坡的壽命預(yù)測(cè)方法。針對(duì)蠕滑邊坡的變形破壞特點(diǎn)，基于運(yùn)用損傷變量原理分析邊坡蠕滑變形特征與位移的演化規(guī)律，提出通過(guò)安全系數(shù)和邊坡的初始彈性變形量得到的位移預(yù)警值作為維納過(guò)程的閾值，利用樣本數(shù)據(jù)采用極大似然估計(jì)的方法確定模型參數(shù)的估計(jì)值，進(jìn)行邊坡剩余壽命預(yù)測(cè)；并將得到的估計(jì)值作為貝葉斯更新的初始值，利用貝葉斯更新隨機(jī)參數(shù)得到考慮后續(xù)邊坡位移值的邊坡剩余壽命快速準(zhǔn)確估計(jì)，并運(yùn)用于相關(guān)的工程。

1 基于維納過(guò)程的建模

1.1 模型描述

基于隨機(jī)過(guò)程的建模方法能夠反映邊坡的真實(shí)退化過(guò)程，其描述為：

X(t)=X(0)+λt+σB(t)

(1)

式中：X(t)為邊坡在t時(shí)刻的性能退化指標(biāo)退化量；X(0)為性能退化指標(biāo)初始值；λ為退化指標(biāo)的漂移系數(shù)，表征邊坡的性能退化速度；B(t)～N(0,t)為標(biāo)準(zhǔn)維納過(guò)程，表征退化動(dòng)態(tài)特性；σ為退化指標(biāo)的擴(kuò)散系數(shù)。本文將漂移系數(shù)λ看成隨機(jī)變量來(lái)表征隨機(jī)因素對(duì)邊坡的影響[10]，即λ～N(μλ,σ2λ)，式(1)可轉(zhuǎn)化為：

X(t)=X(0)+λ(μλ,σ2λ)t+σB(t)

(2)

1.2 模型剩余壽命及可靠度函數(shù)

若退化量呈現(xiàn)遞增趨勢(shì)，其失效閾值為w，采用首達(dá)時(shí)間定義T為性能退化過(guò)程的壽命，即：

T=inf{t:X(t)≥w}

(3)

考慮漂移系數(shù)λ為隨機(jī)變量的情況，剩余壽命T的概率密度函數(shù)、可靠度函數(shù)和期望分別為：

fT(t)=w(2πt3(σ2λt+σ2)exp[-(w-μλt)22t(σ2λt+σ2]

(4)

RT(t)=φ(w-μλtσλ2t2+σ2t)-exp(2wμλσ2+2σ2λw2σ4)×

φ(-2σ2λwt+σ2(w+μλt)σ2σ2λt2+σ2t)

(5)

E(T)=E(E(w│λ))=E(wλ)=

wσ2λexp(-μ2λ2σ2λ) ∫μλ0exp(x22σ2λ)dx=2(w)σλD(μλ2σλ)

(6)

式中：D(z)=exp(-z2)∫z0exp(x2)dx表示關(guān)于實(shí)數(shù)z的Dawson積分。可以看出當(dāng)σ2λ=0，即漂移系數(shù)λ為確定值，上述特性與一般維納過(guò)程完全相同，可見(jiàn)確定性參數(shù)模型是隨機(jī)變量模型的特殊情況。

對(duì)于當(dāng)前時(shí)刻t，剩余壽命可以定義為St=inf{st:X(t+st)≥w)│X(t)

fSt(st)=fst(t+st)

(7)

RSt(st)=1-∫stfst(n)dn

(8)

E(St)=2 [w-X(t)]σλD(μλ2σλ)

(9)

從式(7)～(9)可以看出，當(dāng)t=0時(shí)剩余壽命分布與壽命分布完全相同，即壽命分布是剩余壽命分布的特殊情況。

1.3 閾值的確定

1.3.1 蠕滑型邊坡?lián)p傷變量與其位移變形量分析

根據(jù)損傷力學(xué)基本原理，邊坡的蠕滑與塑性變形損傷過(guò)程可以用損傷變量Dt來(lái)定量刻畫(huà)與評(píng)價(jià)，其損傷變量[11]為：

Dt=1-EtE0

(10)

式中：Dt為邊坡巖土體任意時(shí)刻t對(duì)應(yīng)的損傷變量;E0為初始彈性模量;Et為發(fā)生損傷后的變形模量。材料未損傷時(shí)，Et=E0，Dt=0；材料完全損傷破壞時(shí)，Et=0，Dt=1。

根據(jù)損傷力學(xué)中損傷變量的定義和邊坡壓縮滑移變形受力特點(diǎn)，可確定任意時(shí)刻t對(duì)應(yīng)的邊坡?lián)p傷變量與其位移量的定量關(guān)系：

Dt=1-EtE0=1-σεtσε0=εt-ε0εt=S(t)-S0S(t)

(11)

式中：E0=σε0，為邊坡彈性模量；Et=σεt，為任意時(shí)刻t對(duì)應(yīng)的彈性模量；σ為邊坡滑移面上的應(yīng)力；ε0，εt分別為起始點(diǎn)、任意時(shí)刻t對(duì)應(yīng)的滑移方向的應(yīng)變；S0為邊坡的起始彈性位移量；S(t)為任意時(shí)刻的位移量。

1.3.2 蠕滑型邊坡穩(wěn)定系數(shù)與損傷變量的定量關(guān)系

根據(jù)張道兵等對(duì)于邊坡安全系數(shù)的定義，邊坡在逐漸損傷直到失穩(wěn)破壞過(guò)程中，損傷變量Dt逐漸趨近于容許最大損傷變量Dcr，因此安全系數(shù)被定義為極限損傷變量Dlim與容許最大損傷變量Dc r之比，即：

K=DlimDc r

(12)

由此可知，邊坡穩(wěn)定性系數(shù)F(t)則為邊坡極限損傷變量Dlim與任意時(shí)刻損傷變量Dt之比[6]，即:

F(t)=Dlim/Dc r=1Dt

(13)

式中：Dt為t時(shí)刻的損傷變量；Dlim=1為極限損傷變量。

因此，將式(11)帶入式(13)，位移實(shí)時(shí)監(jiān)測(cè)直接確定蠕滑型邊坡任意時(shí)刻t的動(dòng)態(tài)穩(wěn)定系數(shù)F(t)為：

F(t)=1Dt=S(t)S(t)-S0

(14)

式中：S0為邊坡的起始彈性變形量;S(t)為邊坡任意時(shí)刻t的位移變形量。

1.3.3 基于安全系數(shù)的邊坡穩(wěn)定性位移監(jiān)測(cè)預(yù)警判據(jù)

根據(jù)引水工程邊坡具體情況，依據(jù)水利水電工程邊坡設(shè)計(jì)規(guī)范，可確定邊坡的安全系數(shù)。根據(jù)式(14)，將安全系數(shù)Fc r作為判定邊坡是否穩(wěn)定的判據(jù)準(zhǔn)則，由此可確定邊坡穩(wěn)定性位移監(jiān)測(cè)預(yù)警判據(jù)Sc r[6]：

Fc r=Sc rSc r-S0

(15)

即

Sc r=Fc rS0Fc r-1

(16)

1.4 參數(shù)估計(jì)

假定邊坡在時(shí)刻tf獲得的退化數(shù)據(jù)為{X(ti),1≤ti≤tf}，對(duì)參數(shù)Θ=(μλ0,σλ02,σ2)＇進(jìn)行估計(jì)，其中(μλ0,σλ02)表示漂移系數(shù)λ均值與方差的離線估計(jì)值，令t=(t1,…tf)′，x=(x(t1),…,x(tf))′，那么x服從多變量高斯分布，其均值和協(xié)方差分別為：

μ=μλ0t,Σ=σ2λ0tt＇+Ω

(17)

其中，Ω=σ2，=[t1t1…t1

t1t2…t2

????

t1t2…tf]。

則關(guān)于參數(shù)Θ的對(duì)數(shù)似然函數(shù)為:

(Θ│x)=-f2ln(2π)-12ln|Σ|-

12(x-μ)′Σ-1(x-μ)

(18)

極大化似然函數(shù)(18)就可以得到參數(shù)Θ的最優(yōu)解(μλ0,σ2λ0,σ2)。

1.5 隨機(jī)參數(shù)的在線更新

將tm時(shí)刻所獲得的退化數(shù)據(jù)，即xj=X(tj),j=0,1,…,m,m>f，其對(duì)應(yīng)的退化數(shù)據(jù)集為X=(x0,x1,…,xm)′，Δxj=xj-xj-1對(duì)應(yīng)于從時(shí)刻tj-1到時(shí)刻tj的增量，Δti為時(shí)間間隔，通過(guò)前文參數(shù)估計(jì)確定未知參數(shù)值，當(dāng)漂移系數(shù)λ服從均值為μλ方差為σ2λ的正態(tài)分布時(shí)，在已知X和Θ的條件下，λ的后驗(yàn)估計(jì)仍然服從正態(tài)分布，其均值為μ(m)λ且方差為σ2(m)λ。在貝葉斯框架下，λ的后驗(yàn)分布可以很容易通過(guò)貝葉斯公式進(jìn)行更新，如下所示：

p(λ│X,Θ)∝p(X│λ,Θ)p(λ│Θ)

∝exp[-∑mj=1(Δxj-λΔtj)22σ2Δtj]exp[-(λ-μλ0)22σ2λ0]∝exp

{-12[∑mj=1λ2Δtjσ2-∑mj=1(2λΔxj)σ2+λ2-2λ2μλ0+μλ02σ2λ0]}

∝exp{-12[(tmσ2-1σ2λ0)λ2-2(xmσ2-μ2λ0σ2λ0)λ]}

∝exp{-[λ-(xmσ2λ0+μλ0σ2)(tmσ2λ0+σ2)]22σ2σ2λ0(tmσλ02+σ2)}

根據(jù)λ服從正態(tài)分布的性質(zhì)，可以得到：

p(λ│X,Θ)=12πσ2(m)λexp[(λ-μ(m)λ22σ2(m)λ]

(19)

因此，tm時(shí)刻對(duì)隨機(jī)參數(shù)λ更新后的結(jié)果為：

μ(m)λ=(xmσ2λ0+μλ0σ2)(tmσ2λ0+σ2)

(20)

σ2(m)λ=2σ2σ2λ0(tmσ2λ0+σ2)

(21)

2 實(shí)例工程運(yùn)用

2.1 工程概況

某引水工程戈壁明渠段位于荒漠平原地區(qū)，戈壁明渠穿越的地質(zhì)條件復(fù)雜，存在砂巖段、泥巖砂巖互層段、強(qiáng)膨脹性泥巖段。最大挖深42 m，挖方段地質(zhì)條件復(fù)雜，全斷面為泥巖，發(fā)生滑坡破壞的風(fēng)險(xiǎn)較大，故對(duì)該明渠高邊坡進(jìn)行了相應(yīng)的監(jiān)測(cè)。

2.2 維納過(guò)程模型建立

該明渠左岸邊坡某測(cè)點(diǎn)累計(jì)水平位移數(shù)據(jù)具有典型破壞特征，以此作為代表，將該位移作為邊坡位移性能指標(biāo)進(jìn)行剩余使用壽命預(yù)測(cè)，其實(shí)測(cè)位移性能變化趨勢(shì)如圖1所示。對(duì)邊坡的位移進(jìn)行維納過(guò)程剩余使用壽命建模分析如下。

圖1 邊坡實(shí)測(cè)位移變化曲線

2.2.1 確定閾值

由該邊坡的S(t)曲線可知，其起始位移變化量S0=28.45 mm；依據(jù)《水利水電工程邊坡設(shè)計(jì)規(guī)范》(SL386-2007)，綜合確定該邊坡的臨界安全系數(shù)Fc r=1.1，根據(jù)式(18)確定邊坡穩(wěn)定性位移預(yù)警判據(jù)為Sc r=w=312.95 mm。

2.2.2 維納過(guò)程建模與預(yù)測(cè)

以上述所得到的位移預(yù)警判據(jù)作為本文模型的閾值。以30 d為間隔，將2005年8月至2011年6月該邊坡的退化數(shù)據(jù)作為樣本，采用極大似然估計(jì)得到邊坡維納過(guò)程隨機(jī)參數(shù)模型參數(shù)為Θ=(μλ0,σ2λ0,σ2)′=(0.120 8,0.004 3,0.185 9)′，并進(jìn)行邊坡維納過(guò)程隨機(jī)參數(shù)模型的剩余壽命預(yù)測(cè)，得到其剩余壽命為508.31 d，與實(shí)際維修時(shí)間相比相對(duì)誤差是5.90%，效果良好。

2.2.3 貝葉斯更新

對(duì)于上述維納過(guò)程隨機(jī)參數(shù)模型，當(dāng)獲得新的實(shí)測(cè)位移值時(shí)，需重復(fù)上述求解過(guò)程，計(jì)算復(fù)雜。本文采用貝葉斯更新的方式，以上述的參數(shù)估計(jì)值為初始值，有效利用獲得的新實(shí)測(cè)位移值，實(shí)現(xiàn)快速準(zhǔn)確地貝葉斯參數(shù)更新，代入本文模型進(jìn)行邊坡的剩余壽命計(jì)算。

以2011年7月為例，實(shí)測(cè)獲得位移值為256.63 mm，采用公式(20)、(21)進(jìn)行參數(shù)更新，得到更新后的參數(shù)為(0.120 5，0.000 07，0.185 9)，將此更新后的參數(shù)代入本文模型，預(yù)測(cè)得到該時(shí)刻的剩余壽命為466.14 d，與實(shí)際維修時(shí)間相比的相對(duì)誤差為3.59%，效果理想。

以類(lèi)似的方式，在獲得后續(xù)位移值時(shí)，可以進(jìn)行快速便捷的參數(shù)更新和剩余壽命預(yù)測(cè)。本文依此方法，基于后續(xù)獲得的實(shí)測(cè)位移值，對(duì)各時(shí)刻對(duì)應(yīng)的剩余壽命進(jìn)行了預(yù)測(cè)，預(yù)測(cè)結(jié)果如圖2所示。

圖2 貝葉斯更新維納過(guò)程和常規(guī)維納過(guò)程預(yù)測(cè)壽命結(jié)果

另外，以2.2.2節(jié)重復(fù)求解的方式，也進(jìn)行了相應(yīng)的剩余壽命預(yù)測(cè)，結(jié)果如圖2所示，這兩種方法的誤差曲線如圖3所示。

圖3 貝葉斯更新維納過(guò)程和常規(guī)維納過(guò)程相對(duì)誤差圖

從圖2、圖3中可以看出，本文基于損傷變量的維納過(guò)程隨機(jī)變量模型的預(yù)測(cè)效果良好，同時(shí)貝葉斯更新參數(shù)的方法簡(jiǎn)單方便，本例中預(yù)測(cè)效果略優(yōu)于常規(guī)維納過(guò)程重復(fù)計(jì)算預(yù)測(cè)的結(jié)果。

3 結(jié) 論

(1)運(yùn)用損傷力學(xué)的基本原理，根據(jù)邊坡?lián)p傷變量及其穩(wěn)定性系數(shù)的關(guān)系，依據(jù)邊坡安全系數(shù)得到基于安全系數(shù)和初始彈性變形量S0的邊坡位移預(yù)警判據(jù)Sc r。

(2)以位移預(yù)警判據(jù)Sc r為閾值，代入維納過(guò)程，構(gòu)建了基于損傷力學(xué)的維納過(guò)程隨機(jī)變量模型，所建模型計(jì)算效果理想。

(3)在貝葉斯框架下，實(shí)現(xiàn)對(duì)維納過(guò)程參數(shù)的更新，計(jì)算結(jié)果表明貝葉斯更新方式有效利用了新的實(shí)測(cè)位移值，能夠便捷準(zhǔn)確地實(shí)現(xiàn)參數(shù)估計(jì)和剩余壽命預(yù)測(cè)。