基于貝葉斯理論及MCMC-MH算法推演地基土材料阻尼比的概率分布模型

曹艷梅， 李東偉， 張玉玉， 楊 林

(北京交通大學 土木建筑工程學院，北京 100044)

近年來，交通系統(tǒng)對周圍結(jié)構(gòu)物及其內(nèi)部儀器、設(shè)備等引發(fā)的微振動問題越來越受到重視。由于交通荷載引起的振動在土體中以波動形式進行傳播和衰減，因此振動的預(yù)測和評估需要首先確定場地土的動力特性[1]，而其中土體動參數(shù)的準確與否對振動預(yù)測和評估結(jié)果影響非常大[2]。根據(jù)既有文獻[3]，對地面振動預(yù)測影響最大的土體參數(shù)為土的剪切模量和材料阻尼比。Auersch[4]通過理論分析和試驗研究證明了場地振動時波動幅值隨距離的衰減指數(shù)主要受到材料阻尼和幾何阻尼的影響。土體的材料阻尼是由土骨架中固體顆粒之間的摩擦以及固態(tài)骨架與孔隙流體之間的相對運動而產(chǎn)生的，體現(xiàn)為波動能量的耗散，而幾何阻尼體現(xiàn)為波動的擴散，因此材料阻尼比的大小對振動能量的傳播與衰減起到至關(guān)重要的作用[5-6]。

目前，材料阻尼比主要是通過測試場地中表面波的空間衰減特性而反演得到，反演過程中場地土通常被模擬為各向同性均質(zhì)分層彈性半空間[7]。Lai等[8-9]開發(fā)了一種基于相位和幅值的回歸技術(shù)，同時確定土壤的剪切模量和材料阻尼比。該方法通過幾何衰減因子的迭代更新來計算幾何阻尼，再由衰減系數(shù)反演材料阻尼比分布。Rix等[10]通過解耦剪切波波速和材料阻尼比的方法簡化了反演地基土動參數(shù)的技術(shù)，但是仍然需要計算幾何衰減因子，后來采用與Lai等開發(fā)的技術(shù)類似的方式獲得衰減曲線，進而獲得了地基土材料阻尼比分布。Badsar等基于半功率帶寬法提出了一種確定材料阻尼比的新方法，不需要計算幾何衰減因子，直接將半功率帶寬法應(yīng)用于土壤響應(yīng)的波數(shù)域，有效地避開了不確定性剪切模量分布對幾何衰減因子計算準確性的影響。

然而，由于衰減曲線對小空間尺度和較大深度土壤性質(zhì)的變化不敏感，以及試驗條件的限制和試驗誤差等不確定性因素，地基土的動參數(shù)在一定范圍內(nèi)不是確定的，而是隨機分布的，如果采用確定的地基土參數(shù)對振動進行預(yù)測會出現(xiàn)一定的偏差。Schevenels等[11]研究了地基土動參數(shù)的不確定性對地面振動預(yù)測的影響，通過貝葉斯理論及與試驗頻散曲線構(gòu)造似然函數(shù)，引入馬爾科夫鏈得到了后驗地基土剪切模量概率分布模型。肖張波[12]結(jié)合地質(zhì)統(tǒng)計學和地球物理反演理論，在貝葉斯理論下基于馬爾科夫鏈的擾動對先驗空間隨機搜索，利用地震數(shù)據(jù)約束得到地球儲層物理參數(shù)的后驗概率模型。

為了能夠更加準確地預(yù)測交通系統(tǒng)引起的場地振動，本文建立地基土阻尼比的先驗概率分布模型，根據(jù)試驗衰減曲線構(gòu)建似然函數(shù)，并將其與正演算法獲得的理論衰減曲線進行比較，利用貝葉斯理論和MCMC-MH(Monte Carlo Markov chain-Metropolis)算法得到多組土體參數(shù)后驗樣本數(shù)據(jù)，進行獨立性和收斂性檢驗，得到有效的土體阻尼比后驗概率分布模型。進而應(yīng)用獨立的后驗樣本計算出自由場地的振動響應(yīng)，利用核密度估計得到一定置信度的置信區(qū)間，通過和試驗數(shù)據(jù)進行對比，驗證了土體阻尼比非確定性概率模型的合理性和可靠性。

1 地基土材料阻尼比先驗概率模型的建立

1.1 地基土材料阻尼比非高斯概率分布模型

基于Badsar等研究中的多通道表面波分析法(MASW)現(xiàn)場測試及其提出的半功率帶寬法，可知場地土振動傳遞函數(shù)的頻率-波數(shù)譜可表示為

(1)

(2)

圖1 某場地土表面波衰減曲線

圖2 反演的地基土材料阻尼比分布

從圖2可看出，土層材料阻尼比在不同深度具有隨機分布的特性，通過地震錐穿透原位測試也顯示出材料阻尼比隨深度而變化的特點。本文忽略土體參數(shù)在水平方向上的空間變異性，僅對土體材料阻尼比隨深度變化的不確定性進行研究。由于表面波的振動能量主要集中在地表，因此文中取地表以下6 m的深度范圍來考慮材料阻尼比分布的隨機性，而6 m以下半空間土體的阻尼比參數(shù)與6 m深度處的取值相同。計算時，假定除阻尼比以外的其它土參數(shù)均已確定：土體密度在整個地層中均為1 800 kg·m-3，泊松比υ為1/3；剪切波波速Cs從地表至1.5 m深度地層為128 m/s，1.5～4 m深度為175 m/s，4～6 m及半空間范圍內(nèi)為355 m/s。

由于阻尼比隨深度具有一定的分散性，且所有土層材料的阻尼比參數(shù)不能均用簡單的一個隨機變量來表示，為了更加合理地模擬土參數(shù)分布，本文采用非高斯隨機過程X(h,θ)來模擬[13]，其中θ是隨機維度中的坐標，h是表征地基土深度的豎向坐標。文中采用24個均勻薄層+底部半空間對土體進行建模，每一薄層的厚度取為lc=0.25 m。

基于隨機過程模擬理論[14]，非高斯隨機過程X(h,θ)可以離散化為概率密度函數(shù)pX(θ)和協(xié)方差函數(shù)CX(h)。由圖2所反演出的材料阻尼比可看出，阻尼比的取值范圍在0.03～0.05。由于嚴格的先驗隨機土體模型會引起后驗土體模型的偏差，且較廣的先驗概率模型能包含更多可能的土層剖面，所以本文先驗概率分布模型的概率密度函數(shù)pX(θ)在深度0～6 m內(nèi)采用0.01～0.08的均勻分布，以便能更好地模擬土層參數(shù)隨深度的變化。協(xié)方差函數(shù)與兩個取樣點間的距離有關(guān)，如圖3所示，它能較好地模擬土層參數(shù)隨深度的變化，在地質(zhì)空間統(tǒng)計中有廣泛應(yīng)用。

圖3 非高斯過程的協(xié)方差函數(shù)

協(xié)方差函數(shù)CX(h)可以采用平穩(wěn)的馬特恩協(xié)方差函數(shù)來表示[15]

(3)

式中：KX為第二類修正貝塞爾函數(shù)；Γ為伽馬函數(shù)；lc和v分別為確定隨機過程X(h,θ)變化范圍和平滑度的嚴格正參數(shù),文中取lc=0.25 m和v=1；σ為概率密度函數(shù)pX(θ)的均方差。

1.2 地基土材料阻尼比先驗概率分布模型的生成

對非高斯隨機過程X(h,θ)直接進行Karhunen-Loeve分解，得到的隨機變量ξk(θ)為互不相關(guān)的非高斯隨機變量。由于概率分布較難確定不便于MATLAB軟件程序生成實現(xiàn)，故文中首先將非高斯隨機過程基于Nataf變換理論[16-17]和Gauss-Hermite多項式變換[3]轉(zhuǎn)換為標準高斯隨機過程η(h,θ)，其相關(guān)函數(shù)Rη如圖4所示。

圖4 標準高斯隨機過程的相關(guān)函數(shù)

根據(jù)Karhunen-Loeve分解理論，標準高斯概率分布可以表示為具有隨機系數(shù)的確定性函數(shù)的乘積和的形式，即

(4)

式中：ξk(θ)為標準高斯隨機變量；λk和fk(h)分別為相關(guān)函數(shù)Rη矩陣分解的特征值和特征向量，如圖5所示。

圖5 相關(guān)函數(shù)Rη矩陣分解的特征值及特征向量

式(4)中，隨著階數(shù)的增高，相應(yīng)的能量貢獻不斷降低，為了降低先驗?zāi)Ｐ娃D(zhuǎn)為后驗?zāi)Ｐ瓦^程中的計算量，本文先驗?zāi)Ｐ椭腥〗財嚯A數(shù)M=18，這些多元隨機變量的概率密度函數(shù)為

(5)

通過MATLAB隨機生成M個服從標準高斯分布的隨機數(shù)θi，經(jīng)式(4)及等概率轉(zhuǎn)換原則[18]獲得地基土材料阻尼比的非高斯概率分布模型，即先驗概率模型。從先驗概率模型中提取4組樣本數(shù)據(jù)，并根據(jù)TML-PLM理論及頻率波數(shù)域-半功率帶寬法正演出樣本數(shù)據(jù)對應(yīng)的4組表面波衰減曲線，如圖6所示。因為還未考慮試驗數(shù)據(jù)的約束，從圖6中可以看出先驗材料阻尼比樣本以及對應(yīng)的先驗衰減曲線與圖1、圖2的試驗數(shù)據(jù)都具有較大的差異性。

圖6 土體材料阻尼比先驗概率模型的4組樣本數(shù)據(jù)及先驗衰減曲線

2 材料阻尼比后驗概率模型的建立

2.1 試驗數(shù)據(jù)約束的似然函數(shù)

根據(jù)貝葉斯理論[19-20]，土體材料阻尼比的后驗概率分布模型可以表示為

(6)

式中：π(ξ|E)為受試驗數(shù)據(jù)約束的后驗概率分布模型；p(ξ)為不受試驗數(shù)據(jù)約束的先驗概率分布模型；L(E|ξ)為似然函數(shù)，它是建立在土體模型參數(shù)與試驗數(shù)據(jù)的匹配度之上的?；谙闰?zāi)Ｐ椭械睦碚撍p曲線(見圖6(b))與試驗測得的多組衰減曲線下界的相關(guān)性，構(gòu)建似然函數(shù)L(E|ξ)

圖7 基于試驗數(shù)據(jù)構(gòu)建的頻散曲線范圍

2.2 MCMC-MH算法生成后驗概率分布模型

由于似然函數(shù)L(E|ξ)是一個與試驗數(shù)據(jù)匹配的函數(shù)，因此式(6)不能直接計算后驗概率分布模型，本文聯(lián)合應(yīng)用馬爾科夫鏈蒙特卡洛方法和Metropolis-Hastings算法[21]來獲得符合后驗概率π(ξ|E)的土體剖面。根據(jù)MCMC方法的原理，首先將圖2所示的土體材料阻尼比分布通過NATAF變換生成馬爾科夫鏈的初始狀態(tài)ξ1，其次通過概率轉(zhuǎn)移函數(shù)q(ξj+1|ξj)產(chǎn)生新的狀態(tài)ξj+1,隨后生成均勻分布隨機變量U～Uniform[0,1]，若滿足式(8)的條件則接受轉(zhuǎn)移狀態(tài)ξj+1，否則拒絕該轉(zhuǎn)移狀態(tài)。定義轉(zhuǎn)移函數(shù)q(ξj+1|ξj)為

(8)

(9)

式中，均方差σ用來確定馬爾可夫鏈的步長，σ越大馬爾可夫鏈的擾動越大，收斂越快，但狀態(tài)轉(zhuǎn)移中候選狀態(tài)的接受率會越低。本文取σ=0.12。

對于式(9)，由于π(ξ)無明確的概率計算公式，故在計算時分兩個階段來判斷是否接受候選狀態(tài)ξj+1。首先使用先驗概率p(ξ)來考慮擾動后的新狀態(tài)ξj+1是否通過第一階段，定義

(10)

若第一階段U<α1再進行第二階段，即使用似然函數(shù)L(E|ξ)來考慮是否接受擾動后的新狀態(tài)ξj+1，此時定義

(11)

兩階段過程中，第一階段控制馬爾科夫鏈的擾動達到平穩(wěn)收斂，第二階α2段等于0或1進行試驗數(shù)據(jù)的約束得到后驗概率分布模型。將以上的MCMC-MH算法采用MATLAB編程實現(xiàn)，由此計算出的馬爾科夫鏈前十組的后驗材料阻尼比分布，如圖8所示。

圖8 材料阻尼比的10組后驗樣本數(shù)據(jù)

2.3 收斂性和獨立性檢驗

在MCMC-MH算法下，分別從變量ξi對應(yīng)不同樣本組數(shù)時的均值和均方差兩個方面來判斷推演所得的后驗概率分布模型的馬爾可夫鏈是否收斂。

(12)

(13)

式中，n為狀態(tài)ξi的組數(shù)，每組ξi包含M個相關(guān)變量(M=18)。根據(jù)式(12)和式(13)計算出馬爾科夫鏈的均值mξ和均方差σξ，分別如圖9(a)和圖9(b)所示。從圖中可以看出，隨著樣本組數(shù)n的增多，均值mξ和均方差σξ隨著馬爾科夫鏈的擾動不斷變小，最終逐漸趨于平穩(wěn)。

平穩(wěn)隨機過程的馬爾科夫鏈可通過相關(guān)系數(shù)ρξ來描述樣本之間的相關(guān)性，相關(guān)系數(shù)ρξ=0的樣本即認為是相互獨立的后驗概率模型。相關(guān)系數(shù)和相關(guān)函數(shù)的定義為

(14)

(15)

式中：τ為馬爾科夫鏈上兩組樣本的間隔步數(shù)；Rξ為兩組樣本的相關(guān)函數(shù)；N為樣本總組數(shù)。

圖9 變量ξi不同樣本組數(shù)的均值、均方差以及不同間隔步長的相關(guān)系數(shù)

根據(jù)式(14)和式(15)計算得到相關(guān)系數(shù)，如圖9(c)所示。可以看出，對于本文的算例，樣本間隔步數(shù)大約在500步時相關(guān)系數(shù)接近0，此時即認為樣本之間相互獨立。取馬爾科夫鏈間隔步數(shù)大于500步的狀態(tài)為相互獨立的后驗概率模型，其中10組相互獨立的土體阻尼比后驗概率分布如圖10所示。

3 地基土材料阻尼比后驗概率模型的應(yīng)用

對于N組獨立的土體動參數(shù)的后驗樣本，通過TLM-PML理論可以得到N組地基土自由場豎向位移的傳遞函數(shù)H(x,ω)，即土體表面的豎直位移u(x,ω)與錘擊力p(ω)的比值。本文將距振源不同距離D和不同頻率f的振動位移傳遞函數(shù)通過核密度估計[22]表示為具有概率密度函數(shù)的經(jīng)驗分布，進而用具有一定置信度水平的置信區(qū)間來表示。

圖10 10組相互獨立的后驗材料阻尼比樣本

為了驗證較少數(shù)量獨立樣本的核密度估計是否準確，下面采用非獨立的馬爾科夫鏈上的所有后驗樣本與獨立的100組后驗樣本的核密度估計進行對比，如圖11所示，其中實線為獨立的100組后驗樣本數(shù)據(jù)計算所得，虛線為非獨立的25 000組后驗樣本數(shù)據(jù)計算所得。從圖11中可以發(fā)現(xiàn)，分別用收斂的馬爾科夫鏈中提取的100組獨立樣本與已經(jīng)收斂的非獨立的馬爾科夫鏈后驗樣本所得的概率密度，數(shù)值的大小和形狀趨勢基本是一致。由于獨立樣本數(shù)量少，因而在兩種圖形的峰值附近存在較小的差異性，但是這種差異性對較大置信度的雙側(cè)置信區(qū)間基本無影響。

圖11 自由場表面豎向位移傳遞函數(shù)H(x,ω)的核密度估計

根據(jù)概率密度函數(shù)計算出距振源不同距離D的傳遞函數(shù)H(x,ω)在95%置信度的雙側(cè)置信區(qū)間，并與Badsar等研究中的試驗數(shù)據(jù)進行對比,如圖12所示。

從圖12中可以看出,在后驗概率模型約束的頻率范圍15～70 Hz內(nèi)，95%置信度的置信區(qū)間基本上涵蓋了所有的試驗數(shù)據(jù)。個別頻率的不一致可能是由于試驗中存在的不利影響因素和試驗數(shù)據(jù)的處理誤差等引起的。因此，通過隨機的后驗概率土體模型可以在一定頻率范圍內(nèi)以一定置信度的置信區(qū)間預(yù)測地基土的振動傳遞，可以有效地解決地基土反演問題中存在的非唯一性問題，避免交通荷載環(huán)境振動預(yù)測中的各種不確定性因素的影響。

從圖12中還可以看出，在較低頻率下(約低于15 Hz)，預(yù)測的置信區(qū)間和試驗數(shù)據(jù)差異較大，這是由于低頻時，土體中振動波的波長較大，表面波測試無法測得較深的土層，而理論模型也未考慮土層深度大于6 m時材料阻尼比的空間變化，因此可以通過增加表面波測試的空間分辨率獲得較低頻率情況下的衰減曲線，進而得到低頻范圍內(nèi)可靠的后驗傳遞函數(shù)置信區(qū)間。

當距振源的距離D較遠時，如圖12(c)(D=100 m)，試驗數(shù)據(jù)的高頻階段高于預(yù)測的置信區(qū)間，而距離D=10 m和D=50 m情況下均未發(fā)現(xiàn)高頻階段的差異性，因此可以判定這種遠距離處的高頻階段的差異性主要是由周圍環(huán)境的振動或噪聲等干擾引起的，進一步說明了環(huán)境振動遠場響應(yīng)的預(yù)測可能會存在一定的誤差?？梢酝ㄟ^優(yōu)化表面波測試或使用其他的測試方法來得到更高精度的衰減曲線，可以獲得更大的預(yù)測頻率范圍，縮小一定置信度下的置信區(qū)間，進一步來減少軌道交通引起的振動預(yù)測中的不確定性。

圖12 后驗概率模型傳遞函數(shù)的95%置信區(qū)間(灰色陰影)及試驗數(shù)據(jù)(黑色實線)

4 結(jié) 論

本文基于非高斯隨機過程建立了地基土材料阻尼比的先驗概率模型，通過貝葉斯理論和似然函數(shù)以MCMC-MH算法得到土體動參數(shù)的后驗概率模型，進而以一定置信度的置信區(qū)間預(yù)測了地基土的振動傳遞，并與試驗數(shù)據(jù)進行比較，驗證了本文提出的土體阻尼比非確定性概率模型的合理性和可靠性。通過本文的研究，得到如下主要結(jié)論：

(1)地基土材料阻尼比的隨機分布特征對場地振動的傳播與衰減具有較大的影響，因此建議在場地振動的研究中考慮進地基土參數(shù)取值的非確定性。

(2)本文提出的地基土材料阻尼比后驗概率模型的推演方法是一種非常有效的隨機振動分析方法，可以有效地解決地基土反演問題中存在的非唯一性問題，能夠應(yīng)用于交通環(huán)境場地振動的預(yù)測、評估、振動傳遞和分布的隨機性分析和研究中。

(3)基于多通道表面波分析的半功率帶寬法可以避免不確定性的剪切模量概率分布對求解衰減系數(shù)的影響，能夠準確有效地應(yīng)用于地基土材料阻尼比的反演和隨機分析中。

(4)可以通過擾動系數(shù)來控制馬爾可夫鏈的步長及收斂速度，但由于馬爾科夫鏈的收斂性需要大量的后驗樣本數(shù)據(jù)進行檢驗，所以建議采用平行馬爾科夫鏈并行運算方法來減少計算機的運行時間，提高計算效率。

(5)后驗材料阻尼比土體參數(shù)具有很大的隨機性，基于多組相互獨立的后驗樣本數(shù)據(jù)，從概率的角度以一定置信度的置信區(qū)間預(yù)測地基土的振動傳遞，避免了不確定因素的影響，對環(huán)境振動預(yù)測及控制具有重要意義。