花崗巖殘積土的土水特征曲線參數(shù)概率

羅小艷++劉偉平

摘要：土水特征曲線（SWCC）用于描述非飽和土中含水量與基質(zhì)吸力的關系，在非飽和土力學中具有重要的作用。在SWCC擬合模型中有多個擬合參數(shù)，各值通過實驗方法獲得且具有很大的不確定性。采用貝葉斯理論對擬合參數(shù)的不確定性進行分析，將SWCC擬合參數(shù)作為隨機變量，采用Van Genuchten模型擬合花崗巖殘積土的土水特征曲線試驗數(shù)據(jù)。以馬爾可夫鏈蒙特卡羅方法的延遲拒絕適應性算法獲得模型參數(shù)后驗分布抽樣，獲得了不同置信區(qū)間的SWCC及其對應參數(shù)值。

關鍵詞：貝葉斯理論；土水特征曲線；馬爾可夫鏈蒙特卡羅算法；花崗巖殘積土；置信區(qū)間

中圖分類號：TU411文獻標志碼：A文章編號：16744764（2017）01011206

收稿日期：20160916

基金項目：國家自然科學基金（51468041、 51268046）；江西省自然科學基金（20161BAB203078）；教育部博士點專項基金（20123601110001）

作者簡介：羅小艷（1978），女，主要從事巖土工程可靠度研究，（Email）luoxiaoyan2@126.com。

劉偉平（通信作者），男，副教授，（Email）wpliu@126.com。

Received：20160916

Foundation item：National Natural Science Foundation of China （No.51468041， 51268046）； Natural Science Foundation of Jiangxi （No. 20161BAB203078）； PhD Programs Foundation of Ministry of Education of China （No. 20123601110001）

Author brief：Luo Xiaoyan （1978）， main research interest： geotechnical engineering reliability， （Emial） luoxiaoyan2@126.com.

Liu Weiping （corresponding author）， associate professor， （Email） wpliu@126.com.Probabilistic analysis of the soilwater characteristic

curve of granite residual soil

Luo Xiaoyan1，2，3，Liu Weiping1

（1.School of Civil Engineering and Architecture， Nanchang University， Nanchang 330031， P. R. China；

2. School of Civil Engineering and Architecture， Jiangxi Science and Technology Normal University，

Nanchang 330013， P. R. China；3. The University of Newcastle， NSW 2308， Australia）

Abstract：The soilwater characteristic curve（SWCC）could be used to describe the relationship between the water content and the matrix suction in unsaturated soil and played important roles in unsaturated soil mechanics. Several parameters in fitting models and significant uncertainty of SWCC were obtained by experiment. The uncertainty of SWCC fitting parameters was analyzed using the Bayesian theory. The curvefitting parameters were regarded as the random vectors. And the Van Genuchten model was used to fit the experimental data of SWCC for the granite residual soils. The posterior distributions of model fitting parameters were obtained by the Markov Chain Monte Carlo simulation delayed rejection adaptive metropolis. Confidence intervals of uncertain model parameters of SWCC were obtained by proposed approach.

Keywords：bayesian theory； soilwater characteristic curve； markov chain monte carlo； granite residual soils； confidence interval

土水特征曲線（SWCC）描述土中含水量或飽和度與基質(zhì)吸力之間的關系，可以用來預測非飽和土中的滲透系數(shù)、抗剪強度、體積變化等[13]。由于每個測量點需要長時間去使基質(zhì)吸力平衡，土水特征曲線試驗耗時很長。通常只能在有限的數(shù)據(jù)的基礎上對土水特征曲線進行擬合。為了有效地預測土水特征曲線，不同的土水特征曲線的擬合方法和模型被提出。土水特征曲線廣泛應用到非飽和土工程的數(shù)值分析中，SWCC擬合參數(shù)輸入的可靠性直接影響著數(shù)值結(jié)果的正確性。

由于土水特征曲線受影響因素很多，土的種類、試驗方法、環(huán)境因素等都會使土水特征存在明顯的差異性。土水特征曲線的試驗數(shù)據(jù)與擬合參數(shù)也就具有明顯的不確定性。Phoon等 [4]采用四項Hermite展開式對砂土的土水特征曲線進行概率分析。譚曉慧等[5]、徐全等[6]分別對UNSODA數(shù)據(jù)庫的非飽和土、合肥非飽和土的土水特征曲線進行了概率和敏感性分析。為了分析模型參數(shù)不確定性，擬合參數(shù)可以作為隨機變量[7]。土水特征曲線的計算模型參數(shù)不易確定，該模型的預測精度嚴重影響邊坡穩(wěn)定性分析的可靠性。

目前，在巖土工程隨機分析中，貝葉斯方法主要用于土的抗剪強度[8]、滲流參數(shù)[9]、固結(jié)系數(shù)[10]等問題的分析中。利用實測土水特征曲線試驗數(shù)據(jù)進行擬合參數(shù)的貝葉斯更新分析還不足。本文基于貝葉斯理論分析花崗巖殘積土土水特征曲線擬合參數(shù)的不確定性，根據(jù)試驗測得的土水特征曲線數(shù)據(jù)獲得了Van Genuchten模型擬合參數(shù)的后驗分布。采用馬爾可夫鏈蒙特卡羅方法的延遲拒絕自適應算法進行參數(shù)后驗分布抽樣，分析土水特征曲線置信區(qū)間，對模型參數(shù)的不確定性進行分析。

1土水特征曲線

對于土水特征曲線的研究，學者們提出了多種不同的土水特征曲線模型用于預測基質(zhì)吸力與含水率或飽和度的關系，土中的土水相互作用較為復雜，且多為經(jīng)驗模型。本文采用Van Genuchten[11]提出的土水特征曲線模型進行分析，Van Genuchten（VG）模型表達式為θw=θr+（θs-θr）S（1）

S=1[1+（aψ）n]（1-n-1）（2）式中：θw為土體的體積含水率；θs為飽和含水率；θr為殘余含水率；S為飽和度；ψ為土體的基質(zhì)吸力；a、n為曲線擬合參數(shù)，a與進氣值相關參數(shù)；n為當基質(zhì)吸力超過進氣值時土中水流出率相關參數(shù)。

2貝葉斯方法

貝葉斯理論非常適合應用分析巖土工程問題，尤其在只能獲得有限數(shù)據(jù)的情況。該理論認為任一未知參數(shù)均可以看作成隨機變量，且具有不確定性，可以用概率分布來描述。貝葉斯公式可以表示為p（ζ|D）=p（D|ζ）p（ζ）p（D）（3）式中：p（ζ）為參數(shù)先驗概率；p（ζ|D）為似然函數(shù)，體現(xiàn)模型輸出與觀測數(shù)據(jù)的似然度；P（ζ|D）為參數(shù)后驗分布；p（D）為概率密度函數(shù)歸化常數(shù)；ζ是不確定輸入?yún)?shù)矢量。

根據(jù)選定的土水特征曲線擬合模型，擬合參數(shù)的不確定性可以通過貝葉斯理論估計，將擬合參數(shù)看成隨機變量ζ?？紤]到模型自身和參數(shù)的不確定性以及觀測誤差，真實值與模型預測值之間存在差值，其可以用一個誤差參數(shù)ε來定義，用ε來代表其不確定性。Sm=S（ζ）+ε（4）式中：Sm是測量飽和度；S（ζ）為計算模型的預測飽和度；ε代表模型誤差和測量誤差。

假設ε符合正態(tài)分布，則ε的概率密度函數(shù)為h（ε）=12πσεexp-（ε-με）22σ2ε（5）式中：με為ε均值，本文取με=0；σ2ε為ε方差。

假設擬合參數(shù)ζ=[a，n]滿足已知的先驗分布p（ζ），根據(jù)貝葉斯公式，則在數(shù)組D中的擬合參數(shù)更新后驗分布可以表示為[1213]p（ζ|D）=

c0p（ζ）（2π）-N2σ-Nεexp-N2σ2εJg（ζ|D）（6）式中：c0為歸化常數(shù)；N為測量記錄點數(shù)；p（ζ）為模型參數(shù)的先驗分布，代表主觀判斷；Jg（ζ|D）可以表示為擬合優(yōu)度。Jg（ζ|D）=1NNn=1[Sm（i）-S（ψ；ζ）]2 （7）式中：Sm（i）第i個測量點的飽和度；S（ψ；ζ）相對應模型參數(shù)所得到飽和度輸出值。

土水特征曲線擬合參數(shù)的最優(yōu)值可以通過函數(shù)擬合優(yōu)度最小的方法獲得，也就是后驗密度函數(shù)最大時，最優(yōu)值所對應的誤差方差等于擬合優(yōu)度最小值，即2ε=minJg（ζ|D）=Jg（|D） （8）由于土水特征曲線擬合參數(shù)值具有非負性，選取擬合參數(shù)ζ的先驗分布為對數(shù)分布[14]。假定擬合參數(shù)之間非相關，則參數(shù)的先驗分布的表達式為p（ζ）=∏mi=112πζiσlnζiexp-12ln（ζi）-μlnζiσlnζi2（9）式中：m為隨機變量的個數(shù)； μlnζi和 σlnζi分別為隨機變量對數(shù)均值和標準差。σlnζi=ln1+σζiμζi2（10）

μlnζi=lnμζi-12σ2lnζi（11）3土水特征曲線試驗數(shù)據(jù)

采用GEOExperts應力相關的土水特征曲線壓力儀測得花崗巖殘積土土樣的土水特征曲線數(shù)據(jù)。該設備通過增加或減小基質(zhì)吸力，實現(xiàn)土樣的脫濕和吸濕過程，利用高進氣值陶土板能夠在空氣和水之間起隔膜作用的特點，可測量不同吸力范圍的土水特征曲線。試驗土樣樣本采集于江西省贛州市于都縣金橋村附近的崩崗區(qū)，其物理性質(zhì)指標為：比重2.73、液限40.0%、塑限27.5%、塑性指數(shù)12.5、最大干密度1.76 g/cm3。將通過2 mm孔徑篩的土樣制成干密度為1.45、1.50、1.55、165 g/cm3的試樣；將分別通過0.5、1、2 mm孔徑篩的土樣制成干密度為1.60 g/cm3試樣，對不同干密度和不同顆粒粒徑大小的試樣進行土水特征曲線試驗。本文使用數(shù)據(jù)來源于上述七個土樣的土水特征曲線試驗數(shù)據(jù)。試驗均采用重塑真空抽氣飽和土樣，所有土樣進行一次脫濕試驗，基質(zhì)吸力施加值為0、5、25、50、100、200、300、400 kPa。具體試驗過程限于文章篇幅不再展開，試驗所得的各土樣飽和度與基質(zhì)吸力關系數(shù)據(jù)如圖1中數(shù)據(jù)點所示。

圖1最優(yōu)土水特征曲線與實測數(shù)據(jù)

Fig.1SWCC of optimal parameters and laboratory data4不確定性估計

基于Van Genuchten擬合模型和前述方法，選取上述花崗巖殘積土土水特征試驗所測數(shù)據(jù)點進行貝葉斯分析，所得的土水特征曲線模型擬合參數(shù)的最優(yōu)值及相應的其他值如表1所示。最優(yōu)值所對應的土水特征曲線及土水特征曲線測量數(shù)據(jù)點如圖1所示。圖2為土水特征曲線模型擬合參數(shù)a、n的后驗分布，分布的峰值點對應值（=0.016 6，=1.296 7）即為參數(shù)a、n的最優(yōu)值。表1模型擬合參數(shù)的后驗分布統(tǒng)計值

Table 1Summary of statistics of posterior distributionsμaσa μn σn 2ε0.016 60.016 78×10-71.296 41.296 71.0×10-40.002 6圖2參數(shù)的后驗分布

Fig.2Posterior probability density5馬爾可夫鏈蒙特卡羅方法

當給定先驗分布時，參數(shù)可以通過貝葉斯理論進行更新，對于輸入?yún)?shù)的后驗分布通常不能取得理論值，這時，后驗分布就需通過隨機的方法產(chǎn)生隨機樣本。馬爾可夫鏈蒙特卡羅抽樣方法可以很好地解決這個問題，實現(xiàn)復雜模型貝葉斯公式的數(shù)值計算。MCMC方法的基本思想是建立一個進行抽樣的馬爾可夫鏈，并最終趨于一個穩(wěn)定分析，通過馬爾可夫鏈獲得后驗分布的隨機樣本，再利用蒙特卡羅法求相關值，如期望值、最大后驗分布概率、置信區(qū)間等。這個方法目前常用于后驗分布抽樣，也可應用于非線性、參數(shù)高維、參數(shù)高度相關的模型分析，并可以獲得足夠多的樣本[10，16]，還并能有效處理大量的隨機參數(shù)，且不依賴于先驗分布。MCMC采用的算法有適應性Metropolis算法（Adapive metropolis）、延遲拒絕算法（Delayed rejection）、延遲拒絕適應性算法（delayed rejection adaptive metropolis algorithm， DRAM）等。本文采用延遲拒絕適應性算法（DRAM）[15]，其穩(wěn)定收斂速度快，適應用性強，即使目標分布偏離高斯分布很大時，DRAM也能夠獲得好的抽樣結(jié)果，目前，該算法在巖土工程領域應用較少。DRAM 算法結(jié)合了延遲拒絕算法和適應性Mertopolis算法的優(yōu)點，在抽樣過程中根據(jù)前面所得的抽樣結(jié)果不斷更新提議協(xié)方差，使其不斷逼近目標分布函數(shù)。 Ci=C0i≤i0

sdcov（ζ0，…，ζi-1）+sdδIdi>i0 （12）式中：C0為初始協(xié)方差；sd為取決于空間維數(shù)的縮放因子，sd=2.42d[17]； cov（ζ0，…，ζi-1）為樣本協(xié)方差矩陣；Id代表d維單位矩陣；δ代表小常數(shù)，保證Ci為非奇異矩陣。

協(xié)方差矩陣表達式為cov（ζ0，…，ζk）=1k（ki=1ζiζTi-（k+1）kTk）（13）式中：k=1k+1ki=0ζi。

對于i>i0，將式（13）代入式（12），協(xié)方差滿足遞歸方程，即Cn+1=n-1nCn+sdn（nn-1Tn-1-

（n+1）nTn+ζnζTn+εId）（14）提出初始參數(shù)值ζ0，如果ζ（1）i被接受，此次接受概率公式為α1（ζi-1，ζ（1）i）=min1，p（θ（1）i|D）q1（ζ（1）i，ζi-1）p（ζi-1|D）q1（ζi-1，ζ（1）i）（15）式中：q1為第1次形成的分布。

如果ζ（1）i被拒絕，根據(jù)當前值第2次迭代繼續(xù)，提議第2次分布q2，其接受概率為α2（ζi-1，ζ（2）i）=min1，p（ζ（2）i|D）q1（ζ（2）i，ζ（1）i）q2（ζ（2）i，ζ（1）i，ζi-1）p（ζi-1|D）q1（ζi-1，ζ（1）i）q2（ζi-1，ζ（1）i，ζ（2）i）

×[1-α1（ζ（2）i，ζ（1）i）][1-α1（ζi-1，ζ（1）i）]（16）式中：q2是第2次形成的分布。

這個程序可以不斷地迭代進行，直至獲得足夠的抽樣樣本。

6置信區(qū)間分析

基于MCMC的土水特征曲線置信區(qū)間分析計算具體步驟：

1）確定土水特征曲線擬合參數(shù)向量ζ的先驗分布p（ζ），假定隨機變量的先驗分布為獨立對數(shù)分布，a均值和標準差分別為0.05和0.05；n均值和標準差分別為1.5和0.3；

2）根據(jù)土水特征曲線實測數(shù)據(jù)點確定似然函數(shù)p（D|ζ）；

3）以后驗概率密度函數(shù)P（ζ|D）為目標函數(shù)，采用延遲拒絕適應性算法獲得隨機樣本；

4）刪除收斂前的參數(shù)樣本，將收斂后的樣本作為后驗概率密度函數(shù)的樣本；

5）計算后驗分布的置信區(qū)間。

考慮土水特征曲線擬合參數(shù)的不確定性，不同置信區(qū)間所對應的水土特征曲線和相應擬合參數(shù)可以用前述方法獲得。通過MCMC方法，在分布空間取樣，對樣本點建立迭代條件，抽取收斂于后驗分布函數(shù)的樣本，獲得符合Van Genuchten模型的20 000個土水特征曲線樣本。圖3為產(chǎn)生參數(shù)a、n隨機變量樣本的時程變化圖。

圖3隨機樣本的時程變化圖

Fig.3Evolution processes of random samples去掉收斂前的3 000個樣本，將剩余的17 000個樣本作為后驗分布樣本，即生成了17 000組模型參數(shù)[a，n]。所以，對于一個給定的基質(zhì)吸力，有17 000個相應的飽和度。通過對這17 000個飽和度重新排列，可以獲得不同置信度的土水特征曲線和相關擬合參數(shù)。各置信度對應的擬合參數(shù)值如表2所示。土水特征曲線的擬合參數(shù)a的置信區(qū)間范圍[0.009 3，0.036 1]，n的置信區(qū)間范圍[1.191 9，1447 2]，所得的參數(shù)可用于分析非飽和土工程的可靠度分析。百分點可用來表征擬合參數(shù)的置信區(qū)間，如50%置信區(qū)間為[25百分點，75百分點]，75%置信區(qū)間為[12.5百分點，87.5百分點]。不同置信區(qū)間對應的土水特征曲線如圖4所示，圖中給定了不同置信區(qū)間的上下限。表2不同置信度百分點的擬合參數(shù)值

Table 2Fitting parameter for different

percentiles of SWCC samples百分點an2.50.036 11.447 250.032 11.416 712.50.027 11.375 9250.022 61.334 7750.013 81.247 187.50.011 71.223 5950.010 31.203 097.50.009 31.191 9圖4土水特征曲線的置信區(qū)間

Fig. 4Confidence intervals of SWCC7結(jié)論

根據(jù)貝葉斯理論框架，采用最大似然法，獲得花崗巖殘積土土水特征曲線的最優(yōu)值。建立了基于延遲拒絕自適應（DRAM）算法的土水特征曲線參數(shù)估計方法，該算法具有很好的魯棒性和靈活性。通過對模型參數(shù)后驗分布抽樣，采用后驗分布的穩(wěn)態(tài)隨機樣本可以對土水特征曲線的某一置信水平的區(qū)間進行預測。融合馬爾可夫鏈蒙特卡羅方法為巖土工程中各類模型參數(shù)優(yōu)選和不確定性分析提供有效手段。（致謝：感謝澳大利亞紐卡斯爾大學黃勁松教授對本文的指導和幫助。）

參考文獻：

[1] ASSOULINE S. A model for soil relative hydraulic conductivity based on the water retention characteristic curve [J]. Water Resources Research， 2001， 37（2）：265271.

[2] RAO B H， SINGH D N. Establishing soilwater characteristic curve of a finegrained soil from electrical measurements [J]. Journal of Geotechnical and Geoenvironmental Engineering， 2010， 136（5）： 751754.

[3] FREDLUND M D， FREDLUND D G， WILSON G W. Prediction of the soilwater characteristic curve from grainsize distribution and volumemass properties [C]// 3rd Brazilian Symposium on Unsaturated Soils， Rio de Janeiro， Brazil， 1997： 1323.

[4] PHOON K K， SANTOSO A， CHENG Y G. Probabilistic analysis of soil water characteristic curves from sandy clay loam [C]// Proc.，GeoCongress 2008： The Challenge of Sustainability in the Geoenvironment， Geotechnical Special Publication， 2008，178： 917925.

[5] 譚曉慧，李丹，沈夢芬，等. 土水特征曲線參數(shù)的概率統(tǒng)計及敏感性分析[J]. 土木建筑與環(huán)境工程，2012，34（6）：97103.

TAN X H， LI D， SHEN M F， et al. Probability statistics and sensitivity analysis of parameters of soilwater characteristic curves [J]. Journal of Civil， Architectural & Environmental Engineering， 2012，34（6）：97103.（in Chinese）

[6] 徐全，譚曉慧，辛志宇，等. 土水特征曲線的概率分析[J].水文地質(zhì)工程地質(zhì)，2015，42（3）：7985.

XU Q， TAN X H， XIN Z Y， et al. Probabilistic analysis of the soilwater characteristic curve [J]. Hydrogeology and Engineering Geology， 2015， 42（3）：7985. （in Chinese）

[7] PHOON K K， SANTOSO A， QUEK S T. Probabilistic analysis of soilwater characteristic curves [J]. Journal of Geotechnical and Geoenvironmental Engineering， 2010， 136（3）： 445455.

[8] WANG Y， AU S K， CAO Z J. Bayesian approach for probabilistic characterization of sand friction angles [J]. Engineering Geology， 2010， 114： 354363.

[9] 左自波，張璐璐，程 演，等. 基于MCMC 法的非飽和土滲流參數(shù)隨機反分析 [J]. 巖土力學，2013，34（8）：23932400.

ZUO Z B， ZHANG L L， CHENG Y， et al. Probabilistic back analysis of unsaturated soil seepage parameters based on Markov Chain Monte Carlo method [J].Rock and Soil Mechanics， 2013，34（8）：23932400.（in Chinese）

[10] KELLY R， HUANG J. Bayesian updating for onedimensional consolidation measurements [J].Canadian Geotechnical Journal， 2015， 52： 13181330.

[11] VAN GENUCHTEN M T. A closedform equation for predicting the hydraulic conductivity of unsaturated soils [J].Soil Science Society of America Journal， 1980， 44：892898.

[12] YAN W M， YUEN K V， YOON G L. Bayesian probabilistic approach for the correlations of compression index for marine clays [J].Journal of Geotechnical and Geoenvironmental Engineering， 2009， 135： 19321940.

[13] YUEN K V. Recent developments of Bayesian model class selection and applications in civil engineering [J]. Structural Safety， 2010， 32 ：338346.

[14] HUANG J， GRIFFITHS D V， FENTON G A. System reliability of slopes by RFEM [J].Soils and Foundations， 2010， 50（3）： 343353.

[15] HAARIO H， LAINE M， MIRA A， et al. DRAM： efficient adaptive MCMC [J]. Statistics and Computing， 2006， 16：339354.

[16] ZHANG L L， ZHANG J， ZHANG L M， et al. Back analysis of slope failure with Markov chain Monte Carlo simulation [J]. Computers and Geotechnics， 2010， 37： 905912.

[17] HAARIO H， SAKSMAN E， TAMMINEN J. An adaptive metropolis algorithm [J]. Bernoulli， 2001， 7： 223242.