混凝土徐變的Gamma模型參數(shù)估計研究

陳夢成，楊 超，方 葦，溫清清，許開成，黃 宏

(1.華東交通大學(xué) 土木建筑學(xué)院，江西 南昌 330013；2.省部共建軌道交通基礎(chǔ)設(shè)施性能監(jiān)測與保障國家重點實驗室，江西 南昌 330013)

混凝土結(jié)構(gòu)在土木工程和基礎(chǔ)設(shè)施中應(yīng)用廣泛，服役期間混凝土在長期荷載作用下會產(chǎn)生徐變變形，并且隨時間不斷增長。徐變會引起橋梁結(jié)構(gòu)性能退化，導(dǎo)致其服役時間未到設(shè)計壽命，提前失效。1996年9月26日，帕勞共和國的Koror-Babeldaob橋建成不到19年突然發(fā)生垮塌。專家分析，其倒塌的主要原因是混凝土徐變引起的跨中下?lián)线^大。因此，Ba?ant等[1]在國際混凝土結(jié)構(gòu)協(xié)會會議上呼吁人們要高度重視混凝土徐變及徐變效應(yīng)問題。

目前，國內(nèi)外對于混凝土徐變研究較多[2-7]，但主要集中在確定性徐變問題上。事實上，影響混凝土徐變的因素有內(nèi)在因素，如材料組分和材料性能等，也有外在因素，如環(huán)境與荷載等，這些因素大多是不確定的。另外，混凝土徐變也是隨時間增長的。因此，混凝土徐變發(fā)展過程是一個單調(diào)非負(fù)、隨時間增長的隨機(jī)過程，傳統(tǒng)的確定性計算方法已經(jīng)不適用于混凝土徐變預(yù)測要求。趙人達(dá)等[8]、徐騰飛等[9]、Ma等[10]和余志武等[11]基于隨機(jī)有限元和統(tǒng)計方法對混凝土徐變引起的橋梁變形進(jìn)行了研究，均得到了非常有益的結(jié)果，推動了混凝土徐變研究進(jìn)程，但這些研究成果尚未考慮徐變隨時間的變化。1965年，Benjamin等[12]首次將混凝土徐變看作隨機(jī)過程。1975年，Abdel-Hameed[13]提出采用隨機(jī)Gamma過程建立結(jié)構(gòu)性能隨時間隨機(jī)退化的模型。常用來描述隨機(jī)過程的有預(yù)定義狀態(tài)類別的Markov過程，其屬無自相關(guān)隨機(jī)過程類。Nootwijk[14]闡釋了離散馬爾可夫過程(如Markov鏈)和連續(xù)Markov過程(如Brownian運(yùn)動)與Poisson、Levy和Gamma過程之間的差異性，并指出Gamma隨機(jī)過程是一個單調(diào)、獨立、非負(fù)退化增量的隨機(jī)過程，適合模擬與時間相關(guān)的、漸變的結(jié)構(gòu)性能退化過程，如裂紋、疲勞損傷、磨耗、腐蝕、收縮徐變等引起的結(jié)構(gòu)性能退化過程。Pandy等[15]也進(jìn)一步說明Gamma過程適用于捕捉這類過程。目前，國內(nèi)外使用隨機(jī)Gamma過程模擬隨時間單調(diào)、緩慢增長過程的研究很多，如混凝土收縮和徐變[16-17]、疲勞裂紋擴(kuò)展[18-20]、鋼材銹蝕[21-25]和產(chǎn)品性能老化[19-20,26]等。但是，模擬混凝土徐變發(fā)展過程的論文屈指可數(shù)[16-17]，且尚不夠成熟。

本文首先采用Gamma隨機(jī)過程對混凝土徐變發(fā)展過程進(jìn)行建模，并將模型參數(shù)作為隨機(jī)變量，根據(jù)它們的先驗分布，導(dǎo)出其后驗分布和后驗估計；其次，基于尺度參數(shù)β后驗估計值計算得到形狀參數(shù)冪函數(shù)表達(dá)式中系數(shù)c的后驗矩估計和極大似然估計值；最后運(yùn)用累計誤差平方和最小準(zhǔn)則，結(jié)合應(yīng)用仿真算例，對這幾種估計方法進(jìn)行對比。

1 混凝土徐變模型構(gòu)建

1.1 Gamma過程

混凝土徐變是混凝土在恒定荷載作用下隨時間遞增而單調(diào)增長的變形。將混凝土在t時刻的徐變量定義為X(t)，并且有X(0)=0，其隨時間的增加而逐漸單調(diào)增長。因此，對任意時刻ti、tj，如果tj>ti，則必有X(tj)-X(ti)>0，符合隨機(jī)Gamma過程的特點。這里使用Gamma分布來描述混凝土徐變量的隨機(jī)性，使用隨機(jī)Gamma過程來描述混凝土徐變量隨時間變化的過程。

混凝土徐變量X的概率密度函數(shù)可用Gamma分布表示，具體為

(1)

假定形狀參數(shù)α(t)為非負(fù)、右連續(xù)的實值函數(shù)，且t>0時，有α(0)=0。當(dāng)Gamma過程形狀參數(shù)α(t)>0和尺度參數(shù)β>0時，混凝土徐變隨機(jī)過程{X(t):t≥0}可以看作是一個時間連續(xù)隨機(jī)Gamma過程，且具有以下特征：

(1)當(dāng)概率為1時，有X(0)=0。

(2)對所有τ>t≥0時，有[X(τ)-X(t)]～Ga(α(τ)-α(t),β)，即X(τ)-X(t)服從Gamma分布。

(3)X(t)的增量是獨立的，即對任意時間節(jié)點t12，各增量X(t2)-X(t1)，…，X(tn)-X(tn-1)相互獨立。

因此，t時刻的混凝土徐變量X(t)的概率密度函數(shù)可以表示為

fX(t)(x)=Ga(x|α(t),β)=

(2)

由式(2)可得到混凝土徐變量的期望值和方差為

(3)

由于混凝土徐變?yōu)榉蔷€性過程，因此，它是一個非穩(wěn)態(tài)Gamma隨機(jī)過程。根據(jù)Noortwijk等[27]研究成果，Gamma隨機(jī)過程模型中比例參數(shù)β為常數(shù)，形狀參數(shù)α(t)可表示為冪函數(shù)形式，即

α(t)=c·tb

(4)

式中：c、b為參數(shù)

將式(4)代入式(2)，Gamma模型中未知參數(shù)變?yōu)閏、b和β三個。在通常情況下，Noortwijk[14]建議參數(shù)b取為常數(shù)，它取決于特定的退化過程。當(dāng)b=1時，隨機(jī)Gamma過程為穩(wěn)態(tài)的，單調(diào)過程與時間為線性關(guān)系；若b≠1，隨機(jī)Gamma過程為非穩(wěn)態(tài)的，單調(diào)過程與時間為非線性關(guān)系。

1.2 模型參數(shù)的經(jīng)典估計方法

為使Gamma模型能解決徐變問題，需要運(yùn)用統(tǒng)計學(xué)方法對模型中三個未知參數(shù)進(jìn)行估計[28-29]。構(gòu)建一組數(shù)據(jù)集，它由檢測時間節(jié)點ti和相應(yīng)混凝土徐變量的檢測值xi組成(i=1,2,…,n)，并且0=t0

1.2.1 矩估計法

依據(jù)大數(shù)定理，只要樣本容量足夠大，便可認(rèn)為樣本矩?zé)o限接近總體矩。因此可用樣本一階原點矩估計總體平均值，用樣本二階中心矩估計總體方差。

根據(jù)式(3)，t時混凝土徐變量的期望值和方差可寫為

(5)

當(dāng)參數(shù)b已知時，通過變換，將日歷時間t轉(zhuǎn)化為換算時間z(t)，即z(t)=tb后，非平穩(wěn)隨機(jī)Gamma過程可以很容易轉(zhuǎn)化為平穩(wěn)隨機(jī)Gamma過程。對于所有時間t，如果混凝土徐變量的增量X(t+h)-X(t)的概率分布僅與h(h>0)有關(guān)，那么混凝土徐變量的隨機(jī)變化過程就具有穩(wěn)定的增量。將式(5)中日歷時間t進(jìn)行逆轉(zhuǎn)換，即令t(z)=z1/b，則式(5)變?yōu)?/p>

(6)

i=1,2，…,n。為了數(shù)學(xué)處理上的便利，同時令

(7)

對所有i=1,2,…,n，式(7)中混凝土徐變量增量Di滿足形狀參數(shù)為cwi和尺度參數(shù)為β的Gamma分布，并且D1,D2，…,Dn獨立。注意到Xi、Di、Yi為隨機(jī)變量，xi、δi、yi為相對應(yīng)的檢測值，那么對于每個時間節(jié)點i，Yi的一階和二階矩分別為

(8)

為便于理解，引入總換算時間內(nèi)混凝土徐變量增量平均值定義，即

(9)

聯(lián)立式(8)和式(9)，可得

(10)

(11)

由于E(Yi)=0，那么由式(7)右邊第二項求期望值可得

(12)

根據(jù)式(10)～式(12)，并對等式(11)兩邊求期望值，有

(13)

(14)

(15)

式(14)、式(15)表明，用矩估計法簡單，而且容易計算。但是當(dāng)參數(shù)b未知時，矩估計法將失去效用。

1.2.2 極大似然估計法

現(xiàn)在討論模型形狀參數(shù)c、b和尺度參數(shù)β的極大似然估計方法。首先建立混凝土徐變增量的似然函數(shù)，然后再對其兩邊取自然對數(shù)，并使自然對數(shù)似然函數(shù)取最大值。在混凝土徐變增量服從Gamma分布的總體中，檢測(或抽樣)到一組數(shù)據(jù)樣本δi=xi-xi-1，i=1,2,…,n，獨立且同分布。因此，混凝土徐變增量的似然函數(shù)是由其獨立增量Gamma密度函數(shù)的乘積構(gòu)成的，即

(16)

對式(16)兩邊取自然對數(shù)，有

(17)

為了使式(17)取得極大值，令

(18)

(19)

(20)

式中：ψ(a)為Digamma函數(shù)，由伽馬函數(shù)對數(shù)求導(dǎo)得

(21)

其中，a>0。ψ(a)可使用Bernardo[30]提出的算法精確計算。

(22)

(23)

(24)

(25)

(26)

1.3 模型參數(shù)的后驗估計方法

1.3.1 尺度參數(shù)β的Bayesian后驗分布和估計

若考慮未知參數(shù)的不確定性，則需要使用Bayesian定理[14]。Bayesian定理是為解決概率問題而提供的一種數(shù)據(jù)學(xué)習(xí)方法。記ti時刻的混凝土徐變增量為隨機(jī)變量Di，δi為來自Di的樣本數(shù)據(jù)或檢測數(shù)據(jù)。由前面討論可知，它們?yōu)楠毩⑼植?，均服從Gamma分布，相對應(yīng)的有未知形狀參數(shù)為b、c和尺度參數(shù)β，其中參數(shù)b一般是已知的常數(shù)(與特定的混凝土徐變過程有關(guān))。在Bayesian理論中，c和β不再是未知參數(shù)，而是隨機(jī)變量，可以根據(jù)經(jīng)驗和已完成的試驗或檢測歷史數(shù)據(jù)確定，有明確的概率分布，是先驗分布，可以是有信息分布或無信息分布。假定參數(shù)c和尺度參數(shù)β是相互獨立的，那么根據(jù)Bayesian定理，有

(27)

式中：L(δ1,…,δn|c,β)為檢測數(shù)據(jù)δ1,…,δn的似然函數(shù)，分布形式已由式(16)給出；π(c,β)為獲得新檢測數(shù)據(jù)δi(i=1,2,…,n)前參數(shù)(c,β)的先驗概率密度(先驗分布)；π(c,β|δ1,…,δn)為獲得新檢測數(shù)據(jù)δi(i=1,2,…,n)后(c,β)的后驗概率密度(后驗分布)。

由此可見，一旦獲得新的檢測數(shù)據(jù)樣本δi，就可利用Bayesian定理式(27)更新先驗分布而得到后驗分布。下面討論當(dāng)c為定值時尺度參數(shù)β的后驗分布和估計的解析求解方法。

當(dāng)參數(shù)c已知時，尺度參數(shù)β的先驗分布考慮采用共軛先驗形式，即π(β|c)服從Gamma分布，記為β～Ga(a,λ)，其中a為形狀參數(shù)，λ為尺度參數(shù)。由式(27)有

π(β|c,δ1,…,δn)∝L(δ1,…,δn|c,β)·π(β|c)=

L(δ1,…,δn|c,β)·π(β)=

(28)

(29)

式中：B為未知的隨機(jī)變量。

當(dāng)參數(shù)c未知且為定值時，依據(jù)式(14)或式(22)，尺度參數(shù)β依賴于c，因此，β的形狀參數(shù)和尺度參數(shù)可以分別寫為：a(c)和λ(c)。當(dāng)選擇a(c)=cτb(τ>0)和λ(c)=λ時，并且c取定值時，尺度參數(shù)β的后驗估計值可以表示為

(30)

由式(5)可知，t時混凝土徐變量X(t)的期望值E(X(t))=ctb/β，在檢測到一組新數(shù)據(jù)δi(i=1,2,…,n)后，t時混凝土徐變量預(yù)測值可表示為

(31)

1.3.2 后驗矩估計

若b未知時，由于此時還有2個未知量c和β，矩估計不再適用。

1.3.3 后驗極大似然估計

2 實例應(yīng)用與討論

2.1 素混凝土徐變

Ba?ant收集了全球大量的徐變試驗數(shù)據(jù)并建立徐變數(shù)據(jù)庫，本文挑選其中一根加載齡期為2 d的試件，試驗測得混凝土徐變觀測數(shù)據(jù)xi(i=1,2,…,15)，并且x1=0[31]。

2.1.1 b已知

考慮形狀參數(shù)為時間t的冪函數(shù)形式，即α(t)=c·tb。根據(jù)Ba?ant獲得的徐變試驗數(shù)據(jù)xi，可計算得到徐變增量δi，并且假定b是已知的，b=0.22。分別采用矩估計法(式(14)、式(15))、極大似然估計法(式(22)、式(23))、Bayesian估計法、后驗矩估計法和后驗極大似然估計法進(jìn)行參數(shù)估計。

(1)根據(jù)b=0.22，使用矩估計法、Newton-Raphson迭代法，求解得到參數(shù)c和尺度參數(shù)β的經(jīng)典估計值，并將其代入式(25)，計算得到混凝土徐變隨時間變化數(shù)據(jù)。

(2)根據(jù)b=0.22，用極大似然估計法、Newton-Raphson迭代法，求解得到參數(shù)c和尺度參數(shù)β的經(jīng)典估計值，并將其代入式(25)，計算得到混凝土徐變隨時間變化數(shù)據(jù)。

5種算法下素混凝土徐變模擬結(jié)果比較見圖1。從圖1可知，5種估計方法模擬結(jié)果均與試驗觀測數(shù)據(jù)吻合較好，其中，因式(14)和式(22)的緣故，經(jīng)典矩估計、極大似然估計方法和后驗矩估計、極大似然估計得到的曲線重合。5種算法下damma模型參數(shù)估計見表1。

圖1 5種算法下素混凝土徐變模擬結(jié)果比較

表1 5種算法下Gamma模型參數(shù)估計值

圖2為5種估計算法與試驗觀測值誤差平方的比較，由圖2可知，用經(jīng)典矩估計和極大似然估計與用后驗估計和極大似然估計得到的徐變計算值隨徐變時間增長波動較大，而用Bayesian估計得到的徐變值波動較小，波動相對比較平緩，且Bayesian估計的累計誤差平方和最小(見表1)，因此，用Bayesian估計的效果最好，優(yōu)于其他4種方法。

圖2 5種算法下素混凝土徐變誤差平方比較

2.1.2 b未知

在b未知的情況下，有3個未知量c、b和β，用矩估計法無法求解。此時，4種估算法，只剩下3種，即經(jīng)典極大似然估計法、Bayesian估計法和后驗極大似然估計法。

表2 3種算法下Gamma模型參數(shù)估計值

圖3 3種算法下素混凝土徐變模擬結(jié)果比較

3種參數(shù)估計算法下素混凝土徐變模擬結(jié)果的誤差平方比較見圖4。

圖4 3種算法下素混凝土徐變誤差平方比較

由圖3可以看出，3種估計算法均和試驗觀測值接近，但由圖4和表2可知，經(jīng)典極大似然估計誤差最大，后驗極大似然估計次之，Bayesian估計最小。因此，Bayesian估計和后驗極大似然估計優(yōu)于經(jīng)典極大似然估計，其中Bayesian估計最優(yōu)。

2.2 鋼管混凝土徐變

2006年，陳政清等[32]進(jìn)行一批圓鋼管高強(qiáng)混凝土徐變試驗，選取其中8號鋼管混凝土試件進(jìn)行徐變發(fā)展過程模擬。試件尺寸及材料等參數(shù)參考文獻(xiàn)[32]。在本實例中，仍然假定鋼管混凝土徐變量X(t)～Ga[α(t),β]；徐變量變化的Gamma隨機(jī)過程為時間冪函數(shù)形式，即α(t)=c·tb，其中：c和b均為待確定的確定性參數(shù)；尺度參數(shù)β為隨機(jī)變量。

在模擬鋼管混凝土徐變發(fā)展過程中，采用類似2.1.2節(jié)的算法，用極大似然估計法得到鋼管混凝土徐變Gamma模型的初始參數(shù)c和b的估計值，見表3。假定尺度參數(shù)β的先驗分布仍然服從Ga(a(c),λ(c))分布，形狀參數(shù)a(c)=cτb，并取τ=2>0，尺度參數(shù)λ(c)=3。

表3 3種算法下Gamma模型參數(shù)估計值

3種估計算法的鋼管混凝土徐變隨時間變化的模擬結(jié)果見圖5。圖6為3種參數(shù)估計算法下鋼管混凝土徐變模擬結(jié)果的誤差平方比較。

圖5 3種算法下鋼管混凝土徐變模擬結(jié)果比較

圖6 3種算法下鋼管混凝土徐變誤差平方比較

由圖5可以看出，3種估計算法的模擬結(jié)果均接近試驗觀測值，但由圖6和表3可知，Bayesian估計和后驗極大似然估計的模擬結(jié)果優(yōu)于經(jīng)典極大似然估計，其中尤以Bayesian估計為最佳。

3 結(jié)論

本文根據(jù)混凝土徐變發(fā)展過程特點，建立了基于非穩(wěn)態(tài)Gamma隨機(jī)過程的混凝土徐變隨機(jī)模型，Gamma隨機(jī)過程反映了混凝土徐變量是時間的冪函數(shù)形式。應(yīng)用矩估計法、極大似然估計法和Bayesian估計法，融合試驗觀測數(shù)據(jù)，對Gamma隨機(jī)過程參數(shù)進(jìn)行估計，從而實現(xiàn)了對混凝土徐變過程的模擬。在此基礎(chǔ)上，本文還應(yīng)用Bayesian理論得到了后驗的參數(shù)估計方法，并將經(jīng)典參數(shù)估計法的模擬效果和后驗參數(shù)估計法的效果進(jìn)行比較，得到如下結(jié)論：

(1)提出的五種參數(shù)估計法的模擬效果均較好，均能和試驗觀測數(shù)據(jù)相吻合。

(2)當(dāng)參數(shù)b已知時，后驗的矩估計和極大似然估計的混凝土徐變效果與經(jīng)典的矩估計和極大似然估計的效果完全相同，Bayesian估計的效果優(yōu)于矩估計和極大似然估計的效果。

(3)當(dāng)參數(shù)b為未知時，矩估計法不再適用；Bayesian估計法和后驗極大似然估計法的模擬結(jié)果均優(yōu)于經(jīng)典極大似然估計，其中以Bayesian估計法的模擬結(jié)果為最佳。