有限數據條件下土性參數波動范圍計算方法的有效性分析

田 密， 盛小濤

(1. 湖北工業大學 土木建筑與環境學院， 湖北 武漢 430068； 2. 武漢大學 水工巖石力學教育部重點實驗室， 湖北 武漢 430072；3. 長江科學院 水利部巖土力學與工程重點實驗室， 湖北 武漢 430010)

自然界土體因受復雜地質成因影響呈現出顯著的空間變異性[1]。20世紀70年代，Vanmarcke提出土體剖面的隨機場模型，提出用波動范圍反映空間兩點間土性參數的相關性。在波動范圍內，土性參數是相關的，大于該范圍可認為基本不相關[2]。波動范圍是聯系點變異性和空間平均特性的重要紐帶，也是隨機場模型應用于實際工程的一個重要參數[2～4]。因此，如何準確地計算波動范圍是應用隨機場模型恰當描述土性參數空間變異性的關鍵。

目前，在求解波動范圍方面國內外學者做了大量工作，現有波動范圍計算方法主要有平均零跨法[1]、遞推空間法[1]、相關函數法[5]、擬合方差折減函數法[6]以及回歸模擬法[7]等。近些年也有學者提出了確定土性參數統計特征值(如波動范圍)的貝葉斯方法[8～10]。雖然現有波動范圍計算方法較多，但不同方法的理論根據、可靠性和應用難易程度各不相同。如我國學者閆澍旺等[11]利用各種方法對大量實測的靜力觸探錐尖阻力和側摩阻力數據進行分析，得到遞推空間法、相關函數法及平均零跨法計算結果相差不多，計算結果可信程度比較高。程強等[12]通過數百個土層靜力觸探資料的計算得到遞推空間法和相關函數法是穩定可靠的計算方法，平均零跨法和回歸模擬法使用有其局限性，計算結果常有較大偏差。李小勇和謝康和[8]利用各種方法對太原工程場地粉質粘土層比貫入阻力進行了波動范圍的計算分析，統計模擬法計算結果明顯偏大，且變異性大，計算結果的穩定性差，而遞推空間法、相關函數法以及平均零跨法等方法計算結果相差不多，結果可信度比較高。李鏡培等[13]采用遞推空間法和相關函數法分析了靜力觸探比貫入阻力實測數據，得出遞推空間法與相關函數法的求解結果相當接近的結論，建議在大量工程場地的波動范圍計算時采用遞推空間法。譚曉慧等[14]根據安徽省某工程場地的粉質粘土層6個靜探孔數據資料對求波動范圍的相關函數法和遞推空間法進行了對比研究。

現有波動范圍計算方法準確性與可靠性的對比研究主要基于工程實測數據。然而，巖土工程勘察中所獲得的資料往往十分有限。基于有限數據采用不同方法計算的土性參數波動范圍會存在顯著差異，難以做出比較[15]。其次，天然土體中土性參數真實的波動范圍往往是未知的[16,17]，依據工程場地勘察資料采用各種波動范圍求解方法的可靠性尚待考究。因此，有必要研究有限數據條件下各種波動范圍求解方法的有效性，為合理選擇土性參數波動范圍計算方法提供參考依據，從而為應用隨機場模型準確描述土性參數空間變異性提供所必需的輸入信息。

本文模擬生成靜力觸探試驗數據，通過對有限模擬數據的分析對比了求解波動范圍的平均零跨法、相關函數法、擬合理論方差折減函數法、擬合簡化方差折減函數法及貝葉斯方法的有效性。首先假定靜力觸探試驗錐尖阻力真實的波動范圍，生成錐尖阻力模擬數據。然后，基于有限模擬數據分別采用不同方法計算其波動范圍，并與其真實值對比，以此說明各種計算方法的準確性，探討各種方法的適用性。

1 波動范圍計算方法

波動范圍是描述土性參數空間變異性的重要指標，計算方法主要有平均零跨法、遞推空間法、回歸模擬法、相關函數法、擬合方差折減函數法及貝葉斯方法等。本文主要對比分析平均零跨法、相關函數法、擬合理論方差折減函數法、擬合簡化方差折減函數法及貝葉斯方法的有效性，下面分別介紹這些方法。

(1)平均零跨法

(1)

(2)相關函數法

相關函數法是通過不同類型理論的相關函數ρ(τ)擬合樣本的相關函數，選擇擬合度最優的相關函數積分求解波動范圍[19～21]。應用相關函數法時，其計算結果受所選相關函數類型的影響很大。由于樣本數量的限制，根據測量數據求得的相關系數ρ(τ)在后半段會出現劇烈振蕩(見圖1)，給選取最優相關函數類型帶來了困難。即便采用最優的相關函數，由于ρ(τ)值在后半段劇烈振蕩，估計的波動范圍偏差也比較大。為了提高相關函數法估計波動范圍的可靠性，Uzielli等[22]提出利用ρ(τ)超過Bartlett值的初始部分數據擬合相關函數，從而估計波動范圍。

圖1 相關函數法

(3)擬合理論或簡化方差折減函數法

圖2 hΓ2(h)-h曲線

擬合方差折減函數法是把曲線hΓ2(h)-h上最大值(hΓ2(h))max對于h作為平穩點，即為hmax。然后根據數值(h,hΓ2(h)) (0≤h≤hmax)用理論的方差折減函數或簡化的方差折減函數進行曲線擬合[1](如圖3所示)。

圖3 擬合方差折減函數法

(4)貝葉斯方法

貝葉斯方法充分融合多源信息(包括有限勘探數據和先驗信息)，將未知參數看作隨機變量，認為在獲得測量數據前就已存在一個概率分布，稱之為先驗分布，在測量數據下未知參數的條件分布稱為后驗分布，后驗分布是對未知參數進行統計推斷的依據[23]。令X為所關心的未知參數，基于貝葉斯理論，未知參數X的后驗分布可以表示成[23]

P(X|Data)=K-1P(Data|X)P(X)

(2)

式中：P(Data|X)為給定參數X情況下測量數據的概率密度函數，稱為似然函數；X為未知參數；Data為測量數據；P(X)為參數X的先驗分布，反映了獲取測量數據前對X的認識；K=P(Data)，是與X無關的常數；P(X|Data)為參數X的后驗分布，反映了在先驗信息和測量數據下參數X的更新信息。

2 靜力觸探試驗數據模擬

靜力觸探試驗(Cone Penetration Test,CPT)是原位試驗中最常用的一種測試技術，它可以采集連續的樣本數據，取樣間距較小。最小取樣間距為0.02 m，測量精度非常高，測量誤差可以忽略不計。并且靜力觸探還可以以不同取樣間距采集數據來適應不同成因土性參數波動范圍的要求，廣泛地用于研究土性參數的波動范圍[11～14]。本文以靜力觸探試驗的錐尖阻力為研究對象，模擬錐尖阻力的測量數據。通過模擬數據分析，比較各種波動范圍計算方法的有效性。

ρ(τ)=exp(-2|τ|/λ)

(3)

式中：ρ(τ)為垂直方向上兩點的相關系數；τ為垂直方向上兩點的相對距離。

ξ=lnqN的空間變異性可以定量表征為：

(4)

文獻資料[22]中，CPT錐尖阻力qN的均值與標準差分別取為固定值μ=300，σ=120。取CPT模擬數據的樣本間距ΔD=0.02 m，lnqN波動范圍λ分別取 0.1，0.4，0.8，1.2 m，樣本長度分別取I=D/λ=5，10，20，30，40，50。對不同波動范圍(如λ=0.1，0.4，0.8，1.2 m)，根據式(4)分別在不同樣本長度I條件下重復模擬20套lnqN數據。例如已知qN的均值與標準差，lnqN相關函數為指數型，CPT樣本間距為0.02 m。當lnqN波動范圍λ取 0.1 m時，CPT樣本長度I為5時，由方程(4)可以模擬一套lnqN數據(如圖4)，重復生成20套lnqN隨機模擬數據。

圖4 CPT錐尖阻力模擬數據(λ=0.1 m，樣本長度I=5)

分別采用平均零跨法、相關函數法、擬合理論方差折減函數法、擬合簡化方差折減函數法以及貝葉斯方法計算錐尖阻力波動范圍值。采用貝葉斯方法時，錐尖阻力波動范圍的先驗分布根據其經驗取值范圍取為均勻分布。由文獻[22]和[24]可知，波動范圍經驗取值范圍為[0.1 m, 1.2 m]。基于先驗信息和模擬的CPT 數據，即可得到波動范圍的后驗分布，即方程(2)，具體可參考文獻[10]。本文采用馬爾科夫鏈蒙特卡洛模擬(Markov Chain Monte Carlo Simulation，MCMCS)的Metropolis-Hastings 算法[10]求解后驗分布，生成50000個波動范圍的MCMCS 樣本，然后基于波動范圍的MCMCS隨機樣本估計其后驗均值。

3 波動范圍計算方法有效性分析

3.1 波動范圍估計值比較

分別采用不同方法計算錐尖阻力波動范圍值，并統計20套重復模擬數據計算的波動范圍平均值λm及變異系數Covλm。圖5給出了五種方法計算結果。

圖5 不同波動范圍計算方法估計結果

由圖5還可以看出，不同波動范圍計算方法的適用性受CPT樣本數據量的影響。如圖5a，當錐尖阻力真實波動范圍取λ=0.1 m，CPT樣本長度I=5，10時，相關函數法無法估計波動范圍值。當真實波動范圍λ=0.1 m，CPT樣本長度I=5時，擬合理論方差折減函數法和擬合簡化方差折減函數法無法計算波動范圍。原因在于當樣本長度非常小時，相關函數法中相關系數超過Bartlett值的數據十分有限，當數據量恰好等于1時則無法計算兩點間相關系數，從而導致無法估計波動范圍。擬合理論方差折減函數法和擬合簡化方差折減函數法因難以計算實際的方差折減函數，所以無法計算波動范圍值。然而，當真實波動范圍取λ=0.1 m，CPT樣本長度I為5，10時，貝葉斯方法和平均零跨法可以估計波動范圍值。由此可見，當錐尖阻力真實波動范圍值較小，且靜力觸探試驗樣本數據量較少時，相關函數法、擬合理論方差折減函數法以及擬合簡化方差折減函數法無法計算錐尖阻力波動范圍，而貝葉斯方法和平均零跨法可用于估計巖土參數波動范圍值。但需要注意的是，平均零跨法適用于巖土參數相關結構為高斯型，對于本文錐尖阻力真實相關函數為指數型，平均零跨法計算結果偏差較大。

此外，由圖5還可以看到，隨著CPT樣本長度增加，貝葉斯方法、相關函數法與擬合理論方差折減函數法估計的λm/λ逐漸趨近于1，即λm逐漸逼近真實值λ，而且偏差越來越小。擬合簡化方差折減函數法與平均零跨法雖然也有逼近真實值的趨勢，但總體上估計結果仍偏小，如圖5a所示，當樣本長度I=50時，此時樣本容量n=Iλ/ΔD=250，擬合簡化方差折減函數法與平均零跨法計算的λm/λ仍然偏離1，而其他三種方法則接近于1。由圖5還可知，當靜力觸探試驗樣本長度I≥20時，貝葉斯方法、相關函數法與擬合理論方差折減函數法估計的波動范圍值非常接近。說明當靜力觸探試驗樣本數據充分時，均可采用三種方法計算波動范圍值。

3.2 波動范圍估計結果的不確定性分析

為了進一步說明不同波動范圍計算方法的可靠性，圖6對比分析了貝葉斯方法、相關函數法以及擬合理論方差折減函數法計算波動范圍值的不確定性。因平均零跨法與擬合簡化方差折減函數法估計波動范圍值存在顯著偏差，此節僅對比其他三種方法估計結果的不確定性。

圖6 不同方法估計波動范圍值的不確定性對比

由圖6可知，在相同CPT數據量條件下貝葉斯方法計算的波動范圍的變異系數Covλm最小，相關函數法與擬合理論方差折減函數法計算的變異系數Covλm較為接近。說明貝葉斯方法求解波動范圍值的不確定性較小，計算結果較為準確可靠，相關函數法與擬合理論方差折減函數法計算結果的不確定性較為接近。這是因為相關函數法和擬合理論方差折減函數法則是基于部分測量數據擬合理論相關函數或理論方差折減函數計算波動范圍，所以兩種方法的計算結果相近。而貝葉斯方法則充分融合了多源信息(包括有限觀測數據和先驗信息)，所以計算結果的不確定性較小。綜合波動范圍估計值以及計算結果不確定性的對比可知，在這些方法中貝葉斯方法計算準確性最高、適用性最強，它既適用于數據量較多的情況，也適用于數據量較少的情況。

4 結 論

本文基于模擬的靜力觸探試驗數據對比了計算土性參數波動范圍的平均零跨法、相關函數法、擬合理論方差折減函數法、擬合簡化方差折減函數法以及貝葉斯方法的有效性。通過不同求解方法計算波動范圍值與變異系數的對比，主要得出以下結論：

(1)通過真實相關函數為指數型相關函數的靜力觸探試驗模擬數據分析，不同計算方法確定波動范圍的準確性存在一定差異。貝葉斯方法確定的波動范圍最接近于真實值，其次是相關函數法與擬合理論方差折減函數法，平均零跨法因其是在假定相關函數為高斯型的前提下提出來的，計算結果存在較大誤差，使用時應當注意其假定條件，擬合簡化方差折減函數法因近似的方差折減函數不一定與實際方差折減函數相同，計算結果存在一定偏差。

(2)CPT樣本數據量影響波動范圍計算方法的適用性。當錐尖阻力真實波動范圍較小，且靜力觸探試驗樣本數據量非常有限時，貝葉斯方法和平均零跨法可以估計波動范圍值，但相關函數法、擬合理論方差折減函數法和擬合簡化方差折減函數法無法計算波動范圍。如本文錐尖阻力真實波動范圍取λ=0.1 m，CPT樣本長度I=5，10時，相關函數法無法估計波動范圍。錐尖阻力真實波動范圍λ=0.1 m，且CPT樣本長度I=5時，擬合理論方差折減函數法和擬合簡化方差折減函數法不適用于確定波動范圍。

(3)當靜力觸探試驗數據量較多時，如本文CPT樣本長度I≥20時，貝葉斯方法、相關函數法與擬合理論方差折減函數法估計的波動范圍非常接近，此時均可采用這三種方法計算波動范圍值。

(4)在靜力觸探試驗數據量相同的條件下，相關函數法與擬合理論方差折減函數法計算的不確定性接近，貝葉斯方法因其充分融合了多源信息(包括有限觀測數據和先驗信息)，計算波動范圍值的不確定性最小。綜合波動范圍估計值以及計算結果不確定性的對比可知，貝葉斯方法準確性最高、適用性最強，既適用于數據量較多的情況，也適用于數據量較少的情況。