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

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

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

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

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

1 混凝土徐變模型構建

1.1 Gamma過程

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

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

(1)

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

(1)當概率為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ù)Noortwijk等[27]研究成果，Gamma隨機過程模型中比例參數(shù)β為常數(shù)，形狀參數(shù)α(t)可表示為冪函數(shù)形式，即

α(t)=c·tb

(4)

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

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

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

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

1.2.1 矩估計法

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

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

(5)

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

(6)

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

(7)

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

(8)

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

(9)

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

(10)

(11)

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

(12)

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

(13)

(14)

(15)

式(14)、式(15)表明，用矩估計法簡單，而且容易計算。但是當參數(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ù)的乘積構成的，即

(16)

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

(17)

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

(18)

(19)

(20)

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

(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ù)學習方法。記ti時刻的混凝土徐變增量為隨機變量Di，δi為來自Di的樣本數(shù)據(jù)或檢測數(shù)據(jù)。由前面討論可知，它們?yōu)楠毩⑼植迹腉amma分布，相對應的有未知形狀參數(shù)為b、c和尺度參數(shù)β，其中參數(shù)b一般是已知的常數(shù)(與特定的混凝土徐變過程有關)。在Bayesian理論中，c和β不再是未知參數(shù)，而是隨機變量，可以根據(jù)經驗和已完成的試驗或檢測歷史數(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)更新先驗分布而得到后驗分布。下面討論當c為定值時尺度參數(shù)β的后驗分布和估計的解析求解方法。

當參數(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為未知的隨機變量。

當參數(shù)c未知且為定值時，依據(jù)式(14)或式(22)，尺度參數(shù)β依賴于c，因此，β的形狀參數(shù)和尺度參數(shù)可以分別寫為：a(c)和λ(c)。當選擇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時混凝土徐變量預測值可表示為

(31)

1.3.2 后驗矩估計

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

1.3.3 后驗極大似然估計

2 實例應用與討論

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估計法、后驗矩估計法和后驗極大似然估計法進行參數(shù)估計。

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

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

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

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

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

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

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

2.1.2 b未知

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

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

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

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

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

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

2.2 鋼管混凝土徐變

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

在模擬鋼管混凝土徐變發(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種估計算法的鋼管混凝土徐變隨時間變化的模擬結果見圖5。圖6為3種參數(shù)估計算法下鋼管混凝土徐變模擬結果的誤差平方比較。

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

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

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

3 結論

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

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

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

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