基于損傷力學的維納過程引水工程邊坡剩余壽命預測

馮曉雷， 黃 銘

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

0 引 言

邊坡作為引水工程的重要組成部分，一旦出現滑坡對于引水工程的安全運行有著重大的影響，且后期的維修耗時長、成本高、難度大，因此需要對于重要邊坡進行安全監測、故障診斷和剩余使用壽命預測，以確保引水工程的安全運行[1]。

邊坡的監測位移是表征邊坡安全的重要信息，通過分析邊坡的位移監測資料，來研究滑坡的預警方法，是目前的一個重點研究內容。徐衛亞等[2]提出了以變形變化趨勢為邊坡破壞標準，以改進的灰色-時序分析時變預測模型進行預警預測；許強等[3]通過大量滑坡變形監測數據，提出了加速度的滑坡預警方法與臨滑預警指標。但是目前研究成果多集中于位移和位移速率及其隨時間變化規律方面，很難解釋引起邊坡蠕滑變形與失穩的損傷機制與動因，與邊坡穩定系數相關的失穩破壞位移判據也研究甚少[4]。

關于剩余壽命的預測方法主要包括基于失效機制模型、數據驅動以及融合上述兩類思想的預測方法[5]。其中數據驅動方法的隨機過程常用來描述性能變化呈現出隨機特性的退化過程[6]，隨機過程模型中的維納過程由于具有優良的統計特性是目前應用較為廣泛的一種剩余壽命預測方法, 不同領域的學者利用單一性能指標對目標進行了維納過程建模，獲得了較好的剩余壽命預測結果[7-9]。因此本文利用隨機過程中維納過程的優點對于邊坡的退化過程進行相關分析。

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

1 基于維納過程的建模

1.1 模型描述

基于隨機過程的建模方法能夠反映邊坡的真實退化過程，其描述為：

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

(1)

式中：X(t)為邊坡在t時刻的性能退化指標退化量；X(0)為性能退化指標初始值；λ為退化指標的漂移系數，表征邊坡的性能退化速度；B(t)～N(0,t)為標準維納過程，表征退化動態特性；σ為退化指標的擴散系數。本文將漂移系數λ看成隨機變量來表征隨機因素對邊坡的影響[10]，即λ～N(μλ,σ2λ)，式(1)可轉化為：

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

(2)

1.2 模型剩余壽命及可靠度函數

若退化量呈現遞增趨勢，其失效閾值為w，采用首達時間定義T為性能退化過程的壽命，即：

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

(3)

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

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表示關于實數z的Dawson積分。可以看出當σ2λ=0，即漂移系數λ為確定值，上述特性與一般維納過程完全相同，可見確定性參數模型是隨機變量模型的特殊情況。

對于當前時刻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)可以看出，當t=0時剩余壽命分布與壽命分布完全相同，即壽命分布是剩余壽命分布的特殊情況。

1.3 閾值的確定

1.3.1 蠕滑型邊坡損傷變量與其位移變形量分析

根據損傷力學基本原理，邊坡的蠕滑與塑性變形損傷過程可以用損傷變量Dt來定量刻畫與評價，其損傷變量[11]為：

Dt=1-EtE0

(10)

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

根據損傷力學中損傷變量的定義和邊坡壓縮滑移變形受力特點，可確定任意時刻t對應的邊坡損傷變量與其位移量的定量關系：

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

(11)

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

1.3.2 蠕滑型邊坡穩定系數與損傷變量的定量關系

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

K=DlimDc r

(12)

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

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

(13)

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

因此，將式(11)帶入式(13)，位移實時監測直接確定蠕滑型邊坡任意時刻t的動態穩定系數F(t)為：

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

(14)

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

1.3.3 基于安全系數的邊坡穩定性位移監測預警判據

根據引水工程邊坡具體情況，依據水利水電工程邊坡設計規范，可確定邊坡的安全系數。根據式(14)，將安全系數Fc r作為判定邊坡是否穩定的判據準則，由此可確定邊坡穩定性位移監測預警判據Sc r[6]：

Fc r=Sc rSc r-S0

(15)

即

Sc r=Fc rS0Fc r-1

(16)

1.4 參數估計

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

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

(17)

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

t1t2…t2

????

t1t2…tf]。

則關于參數Θ的對數似然函數為:

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

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

(18)

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

1.5 隨機參數的在線更新

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

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)}

根據λ服從正態分布的性質，可以得到：

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

(19)

因此，tm時刻對隨機參數λ更新后的結果為：

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

(20)

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

(21)

2 實例工程運用

2.1 工程概況

某引水工程戈壁明渠段位于荒漠平原地區，戈壁明渠穿越的地質條件復雜，存在砂巖段、泥巖砂巖互層段、強膨脹性泥巖段。最大挖深42 m，挖方段地質條件復雜，全斷面為泥巖，發生滑坡破壞的風險較大，故對該明渠高邊坡進行了相應的監測。

2.2 維納過程模型建立

該明渠左岸邊坡某測點累計水平位移數據具有典型破壞特征，以此作為代表，將該位移作為邊坡位移性能指標進行剩余使用壽命預測，其實測位移性能變化趨勢如圖1所示。對邊坡的位移進行維納過程剩余使用壽命建模分析如下。

圖1 邊坡實測位移變化曲線

2.2.1 確定閾值

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

2.2.2 維納過程建模與預測

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

2.2.3 貝葉斯更新

對于上述維納過程隨機參數模型，當獲得新的實測位移值時，需重復上述求解過程，計算復雜。本文采用貝葉斯更新的方式，以上述的參數估計值為初始值，有效利用獲得的新實測位移值，實現快速準確地貝葉斯參數更新，代入本文模型進行邊坡的剩余壽命計算。

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

以類似的方式，在獲得后續位移值時，可以進行快速便捷的參數更新和剩余壽命預測。本文依此方法，基于后續獲得的實測位移值，對各時刻對應的剩余壽命進行了預測，預測結果如圖2所示。

圖2 貝葉斯更新維納過程和常規維納過程預測壽命結果

另外，以2.2.2節重復求解的方式，也進行了相應的剩余壽命預測，結果如圖2所示，這兩種方法的誤差曲線如圖3所示。

圖3 貝葉斯更新維納過程和常規維納過程相對誤差圖

從圖2、圖3中可以看出，本文基于損傷變量的維納過程隨機變量模型的預測效果良好，同時貝葉斯更新參數的方法簡單方便，本例中預測效果略優于常規維納過程重復計算預測的結果。

3 結 論

(1)運用損傷力學的基本原理，根據邊坡損傷變量及其穩定性系數的關系，依據邊坡安全系數得到基于安全系數和初始彈性變形量S0的邊坡位移預警判據Sc r。

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

(3)在貝葉斯框架下，實現對維納過程參數的更新，計算結果表明貝葉斯更新方式有效利用了新的實測位移值，能夠便捷準確地實現參數估計和剩余壽命預測。