曹艷梅, 李東偉, 張玉玉, 楊 林
(北京交通大學(xué) 土木建筑工程學(xué)院,北京 100044)
近年來(lái),交通系統(tǒng)對(duì)周圍結(jié)構(gòu)物及其內(nèi)部?jī)x器、設(shè)備等引發(fā)的微振動(dòng)問(wèn)題越來(lái)越受到重視。由于交通荷載引起的振動(dòng)在土體中以波動(dòng)形式進(jìn)行傳播和衰減,因此振動(dòng)的預(yù)測(cè)和評(píng)估需要首先確定場(chǎng)地土的動(dòng)力特性[1],而其中土體動(dòng)參數(shù)的準(zhǔn)確與否對(duì)振動(dòng)預(yù)測(cè)和評(píng)估結(jié)果影響非常大[2]。根據(jù)既有文獻(xiàn)[3],對(duì)地面振動(dòng)預(yù)測(cè)影響最大的土體參數(shù)為土的剪切模量和材料阻尼比。Auersch[4]通過(guò)理論分析和試驗(yàn)研究證明了場(chǎng)地振動(dòng)時(shí)波動(dòng)幅值隨距離的衰減指數(shù)主要受到材料阻尼和幾何阻尼的影響。土體的材料阻尼是由土骨架中固體顆粒之間的摩擦以及固態(tài)骨架與孔隙流體之間的相對(duì)運(yùn)動(dòng)而產(chǎn)生的,體現(xiàn)為波動(dòng)能量的耗散,而幾何阻尼體現(xiàn)為波動(dòng)的擴(kuò)散,因此材料阻尼比的大小對(duì)振動(dòng)能量的傳播與衰減起到至關(guān)重要的作用[5-6]。
目前,材料阻尼比主要是通過(guò)測(cè)試場(chǎng)地中表面波的空間衰減特性而反演得到,反演過(guò)程中場(chǎng)地土通常被模擬為各向同性均質(zhì)分層彈性半空間[7]。Lai等[8-9]開(kāi)發(fā)了一種基于相位和幅值的回歸技術(shù),同時(shí)確定土壤的剪切模量和材料阻尼比。該方法通過(guò)幾何衰減因子的迭代更新來(lái)計(jì)算幾何阻尼,再由衰減系數(shù)反演材料阻尼比分布。Rix等[10]通過(guò)解耦剪切波波速和材料阻尼比的方法簡(jiǎn)化了反演地基土動(dòng)參數(shù)的技術(shù),但是仍然需要計(jì)算幾何衰減因子,后來(lái)采用與Lai等開(kāi)發(fā)的技術(shù)類似的方式獲得衰減曲線,進(jìn)而獲得了地基土材料阻尼比分布。Badsar等基于半功率帶寬法提出了一種確定材料阻尼比的新方法,不需要計(jì)算幾何衰減因子,直接將半功率帶寬法應(yīng)用于土壤響應(yīng)的波數(shù)域,有效地避開(kāi)了不確定性剪切模量分布對(duì)幾何衰減因子計(jì)算準(zhǔn)確性的影響。
然而,由于衰減曲線對(duì)小空間尺度和較大深度土壤性質(zhì)的變化不敏感,以及試驗(yàn)條件的限制和試驗(yàn)誤差等不確定性因素,地基土的動(dòng)參數(shù)在一定范圍內(nèi)不是確定的,而是隨機(jī)分布的,如果采用確定的地基土參數(shù)對(duì)振動(dòng)進(jìn)行預(yù)測(cè)會(huì)出現(xiàn)一定的偏差。Schevenels等[11]研究了地基土動(dòng)參數(shù)的不確定性對(duì)地面振動(dòng)預(yù)測(cè)的影響,通過(guò)貝葉斯理論及與試驗(yàn)頻散曲線構(gòu)造似然函數(shù),引入馬爾科夫鏈得到了后驗(yàn)地基土剪切模量概率分布模型。肖張波[12]結(jié)合地質(zhì)統(tǒng)計(jì)學(xué)和地球物理反演理論,在貝葉斯理論下基于馬爾科夫鏈的擾動(dòng)對(duì)先驗(yàn)空間隨機(jī)搜索,利用地震數(shù)據(jù)約束得到地球儲(chǔ)層物理參數(shù)的后驗(yàn)概率模型。
為了能夠更加準(zhǔn)確地預(yù)測(cè)交通系統(tǒng)引起的場(chǎng)地振動(dòng),本文建立地基土阻尼比的先驗(yàn)概率分布模型,根據(jù)試驗(yàn)衰減曲線構(gòu)建似然函數(shù),并將其與正演算法獲得的理論衰減曲線進(jìn)行比較,利用貝葉斯理論和MCMC-MH(Monte Carlo Markov chain-Metropolis)算法得到多組土體參數(shù)后驗(yàn)樣本數(shù)據(jù),進(jìn)行獨(dú)立性和收斂性檢驗(yàn),得到有效的土體阻尼比后驗(yàn)概率分布模型。進(jìn)而應(yīng)用獨(dú)立的后驗(yàn)樣本計(jì)算出自由場(chǎng)地的振動(dòng)響應(yīng),利用核密度估計(jì)得到一定置信度的置信區(qū)間,通過(guò)和試驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比,驗(yàn)證了土體阻尼比非確定性概率模型的合理性和可靠性。
基于Badsar等研究中的多通道表面波分析法(MASW)現(xiàn)場(chǎng)測(cè)試及其提出的半功率帶寬法,可知場(chǎng)地土振動(dòng)傳遞函數(shù)的頻率-波數(shù)譜可表示為
(1)
(2)
圖1 某場(chǎng)地土表面波衰減曲線
圖2 反演的地基土材料阻尼比分布
從圖2可看出,土層材料阻尼比在不同深度具有隨機(jī)分布的特性,通過(guò)地震錐穿透原位測(cè)試也顯示出材料阻尼比隨深度而變化的特點(diǎn)。本文忽略土體參數(shù)在水平方向上的空間變異性,僅對(duì)土體材料阻尼比隨深度變化的不確定性進(jìn)行研究。由于表面波的振動(dòng)能量主要集中在地表,因此文中取地表以下6 m的深度范圍來(lái)考慮材料阻尼比分布的隨機(jī)性,而6 m以下半空間土體的阻尼比參數(shù)與6 m深度處的取值相同。計(jì)算時(shí),假定除阻尼比以外的其它土參數(shù)均已確定:土體密度在整個(gè)地層中均為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ù)不能均用簡(jiǎn)單的一個(gè)隨機(jī)變量來(lái)表示,為了更加合理地模擬土參數(shù)分布,本文采用非高斯隨機(jī)過(guò)程X(h,θ)來(lái)模擬[13],其中θ是隨機(jī)維度中的坐標(biāo),h是表征地基土深度的豎向坐標(biāo)。文中采用24個(gè)均勻薄層+底部半空間對(duì)土體進(jìn)行建模,每一薄層的厚度取為lc=0.25 m。
基于隨機(jī)過(guò)程模擬理論[14],非高斯隨機(jī)過(guò)程X(h,θ)可以離散化為概率密度函數(shù)pX(θ)和協(xié)方差函數(shù)CX(h)。由圖2所反演出的材料阻尼比可看出,阻尼比的取值范圍在0.03~0.05。由于嚴(yán)格的先驗(yàn)隨機(jī)土體模型會(huì)引起后驗(yàn)土體模型的偏差,且較廣的先驗(yàn)概率模型能包含更多可能的土層剖面,所以本文先驗(yàn)概率分布模型的概率密度函數(shù)pX(θ)在深度0~6 m內(nèi)采用0.01~0.08的均勻分布,以便能更好地模擬土層參數(shù)隨深度的變化。協(xié)方差函數(shù)與兩個(gè)取樣點(diǎn)間的距離有關(guān),如圖3所示,它能較好地模擬土層參數(shù)隨深度的變化,在地質(zhì)空間統(tǒng)計(jì)中有廣泛應(yīng)用。
圖3 非高斯過(guò)程的協(xié)方差函數(shù)
協(xié)方差函數(shù)CX(h)可以采用平穩(wěn)的馬特恩協(xié)方差函數(shù)來(lái)表示[15]
(3)
式中:KX為第二類修正貝塞爾函數(shù);Γ為伽馬函數(shù);lc和v分別為確定隨機(jī)過(guò)程X(h,θ)變化范圍和平滑度的嚴(yán)格正參數(shù),文中取lc=0.25 m和v=1;σ為概率密度函數(shù)pX(θ)的均方差。
對(duì)非高斯隨機(jī)過(guò)程X(h,θ)直接進(jìn)行Karhunen-Loeve分解,得到的隨機(jī)變量ξk(θ)為互不相關(guān)的非高斯隨機(jī)變量。由于概率分布較難確定不便于MATLAB軟件程序生成實(shí)現(xiàn),故文中首先將非高斯隨機(jī)過(guò)程基于Nataf變換理論[16-17]和Gauss-Hermite多項(xiàng)式變換[3]轉(zhuǎn)換為標(biāo)準(zhǔn)高斯隨機(jī)過(guò)程η(h,θ),其相關(guān)函數(shù)Rη如圖4所示。
圖4 標(biāo)準(zhǔn)高斯隨機(jī)過(guò)程的相關(guān)函數(shù)
根據(jù)Karhunen-Loeve分解理論,標(biāo)準(zhǔn)高斯概率分布可以表示為具有隨機(jī)系數(shù)的確定性函數(shù)的乘積和的形式,即
(4)
式中:ξk(θ)為標(biāo)準(zhǔn)高斯隨機(jī)變量;λk和fk(h)分別為相關(guān)函數(shù)Rη矩陣分解的特征值和特征向量,如圖5所示。
圖5 相關(guān)函數(shù)Rη矩陣分解的特征值及特征向量
式(4)中,隨著階數(shù)的增高,相應(yīng)的能量貢獻(xiàn)不斷降低,為了降低先驗(yàn)?zāi)P娃D(zhuǎn)為后驗(yàn)?zāi)P瓦^(guò)程中的計(jì)算量,本文先驗(yàn)?zāi)P椭腥〗財(cái)嚯A數(shù)M=18,這些多元隨機(jī)變量的概率密度函數(shù)為
(5)
通過(guò)MATLAB隨機(jī)生成M個(gè)服從標(biāo)準(zhǔn)高斯分布的隨機(jī)數(shù)θi,經(jīng)式(4)及等概率轉(zhuǎn)換原則[18]獲得地基土材料阻尼比的非高斯概率分布模型,即先驗(yàn)概率模型。從先驗(yàn)概率模型中提取4組樣本數(shù)據(jù),并根據(jù)TML-PLM理論及頻率波數(shù)域-半功率帶寬法正演出樣本數(shù)據(jù)對(duì)應(yīng)的4組表面波衰減曲線,如圖6所示。因?yàn)檫€未考慮試驗(yàn)數(shù)據(jù)的約束,從圖6中可以看出先驗(yàn)材料阻尼比樣本以及對(duì)應(yīng)的先驗(yàn)衰減曲線與圖1、圖2的試驗(yàn)數(shù)據(jù)都具有較大的差異性。
圖6 土體材料阻尼比先驗(yàn)概率模型的4組樣本數(shù)據(jù)及先驗(yàn)衰減曲線
根據(jù)貝葉斯理論[19-20],土體材料阻尼比的后驗(yàn)概率分布模型可以表示為
(6)
式中:π(ξ|E)為受試驗(yàn)數(shù)據(jù)約束的后驗(yàn)概率分布模型;p(ξ)為不受試驗(yàn)數(shù)據(jù)約束的先驗(yàn)概率分布模型;L(E|ξ)為似然函數(shù),它是建立在土體模型參數(shù)與試驗(yàn)數(shù)據(jù)的匹配度之上的?;谙闰?yàn)?zāi)P椭械睦碚撍p曲線(見(jiàn)圖6(b))與試驗(yàn)測(cè)得的多組衰減曲線下界的相關(guān)性,構(gòu)建似然函數(shù)L(E|ξ)
圖7 基于試驗(yàn)數(shù)據(jù)構(gòu)建的頻散曲線范圍
由于似然函數(shù)L(E|ξ)是一個(gè)與試驗(yàn)數(shù)據(jù)匹配的函數(shù),因此式(6)不能直接計(jì)算后驗(yàn)概率分布模型,本文聯(lián)合應(yīng)用馬爾科夫鏈蒙特卡洛方法和Metropolis-Hastings算法[21]來(lái)獲得符合后驗(yàn)概率π(ξ|E)的土體剖面。根據(jù)MCMC方法的原理,首先將圖2所示的土體材料阻尼比分布通過(guò)NATAF變換生成馬爾科夫鏈的初始狀態(tài)ξ1,其次通過(guò)概率轉(zhuǎn)移函數(shù)q(ξj+1|ξj)產(chǎn)生新的狀態(tài)ξj+1,隨后生成均勻分布隨機(jī)變量U~Uniform[0,1],若滿足式(8)的條件則接受轉(zhuǎn)移狀態(tài)ξj+1,否則拒絕該轉(zhuǎn)移狀態(tài)。定義轉(zhuǎn)移函數(shù)q(ξj+1|ξj)為
(8)
(9)
式中,均方差σ用來(lái)確定馬爾可夫鏈的步長(zhǎng),σ越大馬爾可夫鏈的擾動(dòng)越大,收斂越快,但狀態(tài)轉(zhuǎn)移中候選狀態(tài)的接受率會(huì)越低。本文取σ=0.12。
對(duì)于式(9),由于π(ξ)無(wú)明確的概率計(jì)算公式,故在計(jì)算時(shí)分兩個(gè)階段來(lái)判斷是否接受候選狀態(tài)ξj+1。首先使用先驗(yàn)概率p(ξ)來(lái)考慮擾動(dòng)后的新?tīng)顟B(tài)ξj+1是否通過(guò)第一階段,定義
(10)
若第一階段U<α1再進(jìn)行第二階段,即使用似然函數(shù)L(E|ξ)來(lái)考慮是否接受擾動(dòng)后的新?tīng)顟B(tài)ξj+1,此時(shí)定義
(11)
兩階段過(guò)程中,第一階段控制馬爾科夫鏈的擾動(dòng)達(dá)到平穩(wěn)收斂,第二階α2段等于0或1進(jìn)行試驗(yàn)數(shù)據(jù)的約束得到后驗(yàn)概率分布模型。將以上的MCMC-MH算法采用MATLAB編程實(shí)現(xiàn),由此計(jì)算出的馬爾科夫鏈前十組的后驗(yàn)材料阻尼比分布,如圖8所示。
圖8 材料阻尼比的10組后驗(yàn)樣本數(shù)據(jù)
在MCMC-MH算法下,分別從變量ξi對(duì)應(yīng)不同樣本組數(shù)時(shí)的均值和均方差兩個(gè)方面來(lái)判斷推演所得的后驗(yàn)概率分布模型的馬爾可夫鏈?zhǔn)欠袷諗俊?/p>
(12)
(13)
式中,n為狀態(tài)ξi的組數(shù),每組ξi包含M個(gè)相關(guān)變量(M=18)。根據(jù)式(12)和式(13)計(jì)算出馬爾科夫鏈的均值mξ和均方差σξ,分別如圖9(a)和圖9(b)所示。從圖中可以看出,隨著樣本組數(shù)n的增多,均值mξ和均方差σξ隨著馬爾科夫鏈的擾動(dòng)不斷變小,最終逐漸趨于平穩(wěn)。
平穩(wěn)隨機(jī)過(guò)程的馬爾科夫鏈可通過(guò)相關(guān)系數(shù)ρξ來(lái)描述樣本之間的相關(guān)性,相關(guān)系數(shù)ρξ=0的樣本即認(rèn)為是相互獨(dú)立的后驗(yàn)概率模型。相關(guān)系數(shù)和相關(guān)函數(shù)的定義為
(14)
(15)
式中:τ為馬爾科夫鏈上兩組樣本的間隔步數(shù);Rξ為兩組樣本的相關(guān)函數(shù);N為樣本總組數(shù)。
圖9 變量ξi不同樣本組數(shù)的均值、均方差以及不同間隔步長(zhǎng)的相關(guān)系數(shù)
根據(jù)式(14)和式(15)計(jì)算得到相關(guān)系數(shù),如圖9(c)所示??梢钥闯?,對(duì)于本文的算例,樣本間隔步數(shù)大約在500步時(shí)相關(guān)系數(shù)接近0,此時(shí)即認(rèn)為樣本之間相互獨(dú)立。取馬爾科夫鏈間隔步數(shù)大于500步的狀態(tài)為相互獨(dú)立的后驗(yàn)概率模型,其中10組相互獨(dú)立的土體阻尼比后驗(yàn)概率分布如圖10所示。
對(duì)于N組獨(dú)立的土體動(dòng)參數(shù)的后驗(yàn)樣本,通過(guò)TLM-PML理論可以得到N組地基土自由場(chǎng)豎向位移的傳遞函數(shù)H(x,ω),即土體表面的豎直位移u(x,ω)與錘擊力p(ω)的比值。本文將距振源不同距離D和不同頻率f的振動(dòng)位移傳遞函數(shù)通過(guò)核密度估計(jì)[22]表示為具有概率密度函數(shù)的經(jīng)驗(yàn)分布,進(jìn)而用具有一定置信度水平的置信區(qū)間來(lái)表示。
圖10 10組相互獨(dú)立的后驗(yàn)材料阻尼比樣本
為了驗(yàn)證較少數(shù)量獨(dú)立樣本的核密度估計(jì)是否準(zhǔn)確,下面采用非獨(dú)立的馬爾科夫鏈上的所有后驗(yàn)樣本與獨(dú)立的100組后驗(yàn)樣本的核密度估計(jì)進(jìn)行對(duì)比,如圖11所示,其中實(shí)線為獨(dú)立的100組后驗(yàn)樣本數(shù)據(jù)計(jì)算所得,虛線為非獨(dú)立的25 000組后驗(yàn)樣本數(shù)據(jù)計(jì)算所得。從圖11中可以發(fā)現(xiàn),分別用收斂的馬爾科夫鏈中提取的100組獨(dú)立樣本與已經(jīng)收斂的非獨(dú)立的馬爾科夫鏈后驗(yàn)樣本所得的概率密度,數(shù)值的大小和形狀趨勢(shì)基本是一致。由于獨(dú)立樣本數(shù)量少,因而在兩種圖形的峰值附近存在較小的差異性,但是這種差異性對(duì)較大置信度的雙側(cè)置信區(qū)間基本無(wú)影響。
圖11 自由場(chǎng)表面豎向位移傳遞函數(shù)H(x,ω)的核密度估計(jì)
根據(jù)概率密度函數(shù)計(jì)算出距振源不同距離D的傳遞函數(shù)H(x,ω)在95%置信度的雙側(cè)置信區(qū)間,并與Badsar等研究中的試驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比,如圖12所示。
從圖12中可以看出,在后驗(yàn)概率模型約束的頻率范圍15~70 Hz內(nèi),95%置信度的置信區(qū)間基本上涵蓋了所有的試驗(yàn)數(shù)據(jù)。個(gè)別頻率的不一致可能是由于試驗(yàn)中存在的不利影響因素和試驗(yàn)數(shù)據(jù)的處理誤差等引起的。因此,通過(guò)隨機(jī)的后驗(yàn)概率土體模型可以在一定頻率范圍內(nèi)以一定置信度的置信區(qū)間預(yù)測(cè)地基土的振動(dòng)傳遞,可以有效地解決地基土反演問(wèn)題中存在的非唯一性問(wèn)題,避免交通荷載環(huán)境振動(dòng)預(yù)測(cè)中的各種不確定性因素的影響。
從圖12中還可以看出,在較低頻率下(約低于15 Hz),預(yù)測(cè)的置信區(qū)間和試驗(yàn)數(shù)據(jù)差異較大,這是由于低頻時(shí),土體中振動(dòng)波的波長(zhǎng)較大,表面波測(cè)試無(wú)法測(cè)得較深的土層,而理論模型也未考慮土層深度大于6 m時(shí)材料阻尼比的空間變化,因此可以通過(guò)增加表面波測(cè)試的空間分辨率獲得較低頻率情況下的衰減曲線,進(jìn)而得到低頻范圍內(nèi)可靠的后驗(yàn)傳遞函數(shù)置信區(qū)間。
當(dāng)距振源的距離D較遠(yuǎn)時(shí),如圖12(c)(D=100 m),試驗(yàn)數(shù)據(jù)的高頻階段高于預(yù)測(cè)的置信區(qū)間,而距離D=10 m和D=50 m情況下均未發(fā)現(xiàn)高頻階段的差異性,因此可以判定這種遠(yuǎn)距離處的高頻階段的差異性主要是由周圍環(huán)境的振動(dòng)或噪聲等干擾引起的,進(jìn)一步說(shuō)明了環(huán)境振動(dòng)遠(yuǎn)場(chǎng)響應(yīng)的預(yù)測(cè)可能會(huì)存在一定的誤差??梢酝ㄟ^(guò)優(yōu)化表面波測(cè)試或使用其他的測(cè)試方法來(lái)得到更高精度的衰減曲線,可以獲得更大的預(yù)測(cè)頻率范圍,縮小一定置信度下的置信區(qū)間,進(jìn)一步來(lái)減少軌道交通引起的振動(dòng)預(yù)測(cè)中的不確定性。
圖12 后驗(yàn)概率模型傳遞函數(shù)的95%置信區(qū)間(灰色陰影)及試驗(yàn)數(shù)據(jù)(黑色實(shí)線)
本文基于非高斯隨機(jī)過(guò)程建立了地基土材料阻尼比的先驗(yàn)概率模型,通過(guò)貝葉斯理論和似然函數(shù)以MCMC-MH算法得到土體動(dòng)參數(shù)的后驗(yàn)概率模型,進(jìn)而以一定置信度的置信區(qū)間預(yù)測(cè)了地基土的振動(dòng)傳遞,并與試驗(yàn)數(shù)據(jù)進(jìn)行比較,驗(yàn)證了本文提出的土體阻尼比非確定性概率模型的合理性和可靠性。通過(guò)本文的研究,得到如下主要結(jié)論:
(1)地基土材料阻尼比的隨機(jī)分布特征對(duì)場(chǎng)地振動(dòng)的傳播與衰減具有較大的影響,因此建議在場(chǎng)地振動(dòng)的研究中考慮進(jìn)地基土參數(shù)取值的非確定性。
(2)本文提出的地基土材料阻尼比后驗(yàn)概率模型的推演方法是一種非常有效的隨機(jī)振動(dòng)分析方法,可以有效地解決地基土反演問(wèn)題中存在的非唯一性問(wèn)題,能夠應(yīng)用于交通環(huán)境場(chǎng)地振動(dòng)的預(yù)測(cè)、評(píng)估、振動(dòng)傳遞和分布的隨機(jī)性分析和研究中。
(3)基于多通道表面波分析的半功率帶寬法可以避免不確定性的剪切模量概率分布對(duì)求解衰減系數(shù)的影響,能夠準(zhǔn)確有效地應(yīng)用于地基土材料阻尼比的反演和隨機(jī)分析中。
(4)可以通過(guò)擾動(dòng)系數(shù)來(lái)控制馬爾可夫鏈的步長(zhǎng)及收斂速度,但由于馬爾科夫鏈的收斂性需要大量的后驗(yàn)樣本數(shù)據(jù)進(jìn)行檢驗(yàn),所以建議采用平行馬爾科夫鏈并行運(yùn)算方法來(lái)減少計(jì)算機(jī)的運(yùn)行時(shí)間,提高計(jì)算效率。
(5)后驗(yàn)材料阻尼比土體參數(shù)具有很大的隨機(jī)性,基于多組相互獨(dú)立的后驗(yàn)樣本數(shù)據(jù),從概率的角度以一定置信度的置信區(qū)間預(yù)測(cè)地基土的振動(dòng)傳遞,避免了不確定因素的影響,對(duì)環(huán)境振動(dòng)預(yù)測(cè)及控制具有重要意義。