基于粒徑分布曲線的非飽和砂土土水特征曲線概率預測模型*

張 潔 陽 帥 譚澤穎 李 旭

(①同濟大學地下建筑與工程系，上海 200092，中國)(②北京交通大學土木建筑工程學院，北京 100044，中國)

0 引 言

土水特征曲線是研究非飽和土學的重要手段，定義了非飽和土的基質吸力和含水量之間的關系(李志清等，2007)，反映了土的持水能力特性。土水特征曲線廣泛應用于巖土工程，環境工程以及農業工程。土水特征曲線可用來推斷非飽和土的滲透系數、抗剪強度(Vanapalli et al.，1996; 高游等，2017)等性質，對涉及到非飽和土的工程具有重要意義(包承綱，2004)，目前對于土水特征曲線的影響因素開展了許多研究(石振明等，2018；李同錄等，2019)。為方便工程應用，很多學者提出不同的擬合方程來描述土水特征曲線，如VG模型(van Genuchten，1980)、FX模型(Fredlund et al.，1994)、Gardner模型(Gardner et al.，1958)等。其中：VG模型由于參數較少、形式簡單，是目前使用最廣泛的土水特征曲線模型之一。

非飽和土土水特征曲線可通過室內試驗獲得(Li et al.，2014)，但相關試驗一般比較費時且費用高昂，對非飽和土力學的工程應用帶來了挑戰。另一方面，由于土壤粒徑分布曲線和土水特征曲線形狀相似，一些學者發現可基于土體粒徑分布曲線來預測土體的土水特征曲線。其中：Arya和Paris在1982年(Arya et al.，1982)首次嘗試通過粒徑分布曲線參數預測土水特征曲線(A & P模型)，將顆粒累積分數轉化為含水率，將顆粒半徑通過經驗系數α轉化為孔隙半徑，并通過毛細方程計算吸力。一些研究者還在A & P模型的基礎上提出一些物理或概念模型對土水特征曲線進行預測(如Mohammadi et al.，2011)。這些研究結果表明，可以基于一些易于獲得的物理指標，如粒徑分布，通過一定的轉換方法來估算土的土水特征曲線，從而實現用較低的代價來獲得非飽和土的主要工程性質。上述方法雖然為快速預測土體土水特征曲線提供了一個很好的思路，但本質上是一種經驗方法，不可避免地存在預測誤差。目前，現有研究中極少對上述誤差進行估算，導致不少工程人員直接將由粒徑分布曲線估算獲得的土水特征曲線作為真實的土水特征曲線，無法考慮土水特征曲線估算誤差的影響。

本文將以VG模型為例，提出一種基于土體粒徑分布曲線的土水特征曲線概率預測方法。與傳統方法相比，本文提出的方法可以全面考慮土水特征曲線的估算誤差以及由此帶來的不確定性。本文主要分為以下幾個部分：首先建立非飽和砂土的數據庫，利用極大似然估計法擬合出VG模型的參數；然后，利用粒徑分布曲線參數和孔隙率等參數建立VG模型參數的線性方程；進而，通過回歸分析建立土水特征曲線的經驗公式；最后，對建立的概率分布模型進行驗證。本文提出的方法對非飽和土土水特征曲線的預測和非飽和土可靠度分析均具有重要的意義。

1 土水特征曲線模型

土水特征曲線定義了非飽和土的基質吸力和含水量之間的關系，反映了土的持水能力特性。在VG模型(van Genuchten，1980)中，其表達式為：

(1)

式中：S為土的體積飽和度；ψ為吸力；a、m、n均為擬合參數。一般而言，a與土的進氣狀態有關，n與土的孔徑分布有關，參數m與土體特征曲線的整體對稱性有關。為簡便計算，m可通過參數n基于如下公式進行估算(van Genuchten，1980)：

(2)

基于式(2)，VG模型的模型參數可進一步減少為2個，即a和n。土水特征曲線模型參數可通過最小二乘法等方法來擬合。本文將利用最小二乘法對參數進行擬合。

圖1給出了某土樣實測土水特征曲線值?；谏鲜鲎钚《朔?，其最優模型參數值為a=0.032、n=2.15、m=0.46，模型的擬合優度為0.99。圖中也同時給出了基于上述參數獲得的預測土水特征曲線。由圖可知，實測曲線與由VG模型獲得的擬合曲線十分接近，由此根據擬合參數a、n、m確定土水特征曲線的參數相比于傳統圖解法能更精確(潘登麗等，2020)。

圖1 土水特征曲線擬合圖

2 基于粒徑分布曲線的土水特征曲線參數概率預測方法

2.1 土水特征曲線及粒徑分布曲線參數的相關性分析

已有研究表明，土體土水特征曲線參數與粒徑分布曲線參數具有高度的相關性(Cosby et al.，1984；Saxton et al.，1986)。Li et al.(2014)指出土體土水特征曲線與土體粒徑分布曲線的d10、d30、d60、不均勻系數Cu、曲率系數Cc等特征參數具有密切關系。

根據《巖土工程勘察規范》(中華人民共和國國家標準編寫組，2009)定義，砂土指粒徑大于2 mm的顆粒質量不超過總質量的50%且粒徑大于0.075 mm的顆粒質量超過總質量的50%的土。為研究砂土土水特征曲線與土體粒徑分布曲線之間的關系，本文首先從美國農業安全部UNSODA數據庫(Leij，1996)中選取了100組砂土數據的脫濕土水特征曲線及土壤物理性質參數。圖2和圖3給出了這100組砂土數據的土水特征曲線及土粒徑分布曲線。

圖2 土水特征曲線

圖3 土樣粒徑分布曲線

為建立砂土的土水特征曲線經驗預測模型，論文首先采用最小二乘法對上述100組土水特征曲線進行了擬合。表1給出了100組砂土的土體粒徑分布曲線擬合參數與VG模型參數的相關性分析。

表1 VG模型參數相關性分析

2.2 基于最大似然法的回歸分析

首先介紹lna的回歸方程標定方法。為方便表述，令X={lnd10，lnd30，lnd60}，令Y=lna。假設Y和X1，X2，X3的關系可用如下方程表示：

Y=b0+b1X1+b2X2+b3X3+ε

(3)

式中：ε為一均值為0、標準差為σ的正態分布隨機變量。上述方程中，未知參數包括B={b0，b1，b2，b3，σ}。當B已知的條件下，基于式(3)，Y的均值和標準差分別為：

μY=b0+b1X1+b2X2+b3X3

(4)

σY=σ

(5)

令X1i，X2i，X3i代表第i組土樣的lnd10，lnd30，lnd60值，令δi代表第i個土樣的Y值。給定B條件下，觀測到δi的概率為：

(6)

式中：φ為一標準正態分布的概率密度函數。假設存在s組土樣，令δ={δ1，δ2，…，δs}。當B已知的條件下，觀測到所有數據的概率為：

(7)

上述表達式即為回歸方程參數B的似然函數。根據最大似然性原理(趙軍圣等，2010)，使得式(5)取值最大的B即為模型參數的最優值。實際優化過程中，式(5)中由于連乘的存在會導致似然性函數的數值絕對值很小，在優化中容易受到計算誤差的影響而導致算法難以收斂。此時，可對似然函數取對數，通過對似然函數對數值取最大值而獲得模型參數的最優值。

采用相同的方法，也可對lnn的預測方程進行標定。

2.3 VG模型參數的概率分布預測

下面將采用上述最大似然估計法對lna、lnn進行預測。為方便表述，令θ={lna，lnn}。令μ1、μ2分別代表lna、lnn的均值，令σ1、σ2分別代表lna、lnn的標準差。首先采用最大似然法基于式(7)對lna的回歸方程進行標定，可得B={0.32、-0.27、1.5、-0.19、1.49}。將上述參數帶入方程(4)，可得μ1的計算公式為：

(8)

基于式(5)，可得σ1=1.49。μ1代表了lna的最可能取值，或最優預測值，σ1衡量了預測lna過程中的變異性。圖4a給出了采用式(8)對數據庫中100組土樣預測lna與實際lna的對比圖。通過計算，兩者的相關系數為0.79，表示計算值與預測值具有較好的相關性。

圖4 預測值與實際值對比

采用最大似然法對lnn進行標定，可得其均值預測公式為：

(9)

lnn的標準差為σ2=0.43。圖4b給出了基于式(9)獲得的lnn的預測值與實測值的對比圖。由圖中可知，兩者之間的相關系數為0.58。相比而言，lnn預測值與實測值之間的相關性要低一些，表明lnn更難預測。

通過上述分析獲得了lna和lnn的邊緣概率分布。為獲得lna和lnn之間的聯合概率分布，還需獲得lna和lnn之間的相關性。為獲取上述相關性，可對100個土樣lna、lnn預測值與實測值之間的殘差進行分析。令δ1i、δ2i分別代表第i個土樣的lna、lnn的實測值，令e1i、e2i分別代表第i個土樣關于lna和lnn的預測誤差。第i個土樣關于lna和lnn的殘差可按下式計算：

(10)

圖5a、圖5b分別給出了100組土樣lna、lnn的預測殘差直方圖。圖中同時給出了根據直方圖擬合獲得的正態分布曲線。Kolmogorov-Smirnov檢驗(侯澍旻等，2007)表明，在0.05顯著性水平下，lna、lnn的預測殘差均服從正態分布。

圖5 模型參數殘差直方圖

圖6給出了lna、lnn殘差的散點圖。由圖可知，lna、lnn殘差之間呈正相關關系，兩者之間的相關系數為0.41，即ρ12=0.41。

圖6 ln a-ln n殘差散點圖

令μθ代表θ的均值，即令Cθ代表θ的協方差。Cθ可按下式計算：

(11)

由于lna、lnn分別服從正態分布，θ={lna，lnn}服從多元正態分布?；诙嘣龖B分布特性，θ的聯合概率密度函數為：

(12)

當θ={lna，lnn}服從多元正態分布時，{a，n}服從多元對數正態分布。在實際應用中，可先根據式(12)生成θ={lna，lnn}的樣本，再通過指數變換獲得{a，n}的樣本，由此獲得對應的土水特征曲線。

在式(12)中，μθ={μlna，μlnn}代表了θ的最可能取值，由此生成的土水特征曲線為最可能的土水特征曲線；Cθ給出了θ的估算誤差，界定了θ可能的取值范圍。因此，本文提出的方法不但可提供土水特征曲線的最可能位置，還能提供土水特征曲線可能的取值范圍。

3 模型驗證

為進一步驗證方法的可靠性，本文搜集了文獻外30組砂土的脫濕土水特征曲線及其粒徑分布曲線，采用本文方法對其土水特征曲線進行了預測。圖7給出了其中4組土樣的粒徑分布曲線，圖8給出了對應土樣的實測土水特征曲線。

圖7 土樣粒徑分布曲線

圖8 土樣實測土水特征曲線

以土樣1為例，根據該土樣粒徑分布曲線，其d10、d30、d60，Cu、Cc參數分別為0.09 mm、0.18 mm、0.24 mm、2.72、0.70。將上述參數帶入式(8)、式(9)，可得μθ={-1.36，0.47}。該分布的協方差為式(11)，Cθ={2.231，0.266；0.266，0.188}。根據土水特征曲線參數的均值，圖9a中給出了土樣1的最可能土水特征曲線。由圖可知，最可能土水特征曲線與實際土水特征曲線具有相似性，但不完全一致，給定吸力水平下，最可能土水特征曲線與實測特征曲線飽和度最大差別為0.091。

圖9 驗證土樣土水特征曲線及預測曲線

將μθ、Cθ代入式(12)，即可獲得土樣土水特征曲線參數 {lna，lnn}的分布。根據該分布，獲得了 {lna，lnn}的50個樣本，由此生成了50條土水特征曲線。作為比較，圖中也給出了上述50條土水特征曲線。由圖可知，由于lna、lnn預測誤差的存在，土水特征曲線也存在多種可能性，根據粒徑分布曲線無法唯一確定土水特征曲線。

另一方面，與圖2中100條砂土曲線的土水特征曲線相比，圖9a中土水特征曲線的離散性有顯著降低，且實測土水特征曲線位于預測土水特征曲線范圍之內，說明粒徑分布曲線顯著降低了土水特征曲線的不確定性。

對于圖9a中的每一個吸力，50條預測曲線分別對應50個飽和度。根據上述50個飽和度數值，可分別獲得該吸力水平下90%的置信區間。改變不同的吸力值，由此可獲得不同吸力水平下飽和度的90%置信區間。圖9a中比較了實測的土水特征曲線以及不同吸力水平下飽和度的90%置信區間。由圖可知，實測土水特征曲線位于預測土水特征曲線的90%置信區間內。

采用類似的方法，圖9b、圖9c和圖9d分別對其他3種土樣的實測土水特征曲線與預測土水特征曲線進行了比較。由圖可知，另外3種土樣的實測土水特征曲線雖與預測最可能土水特征曲線不完全一致，但都位于預測土水特征曲線90%置信區間范圍內。

為對本文提出的方法進行進一步的驗證，采用本文方法對驗證數據庫中30個土樣的土水特征曲線進行了預測。為方便比較，對每條實測土水特征曲線按其90%的上下界進行了歸一化，具體方法介紹如下。令φi代表某一具體的吸力水平，令Si，SiU和SiL分別代表其實測飽和度、實測飽和度90%置信區間上界、實測飽和度90%置信區間下界。與φi對應的歸一化飽和度為：

(13)

采用類似的方法，也可對該吸力水平下飽和度的90%置信區間上下界進行標準化。由式(13)可以看出，歸一化后飽和度90%置信區間上界的值為1、下界的值為0。按上述方法，圖10給出了上述30組土樣的歸一化的土水特征曲線。由圖可知，所有土水特征曲線的歸一化吸力值均位于0～1之間，表明所有實測土水特征曲線均位于預測土水特征曲線90%置信區間內。

圖10 土水特征曲線置信區間

總體而言，由于土體的孔隙結構對土水特征曲線的形態具有絕對性的影響，土的持水能力是其微觀孔隙分布的宏觀表現(胡冉等，2013)，因此土的粒徑分布為土水特征曲線的預測提供有用的信息，可用來減少土水特征曲線的不確定性。另一方面，即使級配相同的土體，其在不同的應力、干濕循環歷史作用下，會形同不同的孔隙結構和儲水狀態，對土水特征曲線也存在一定的影響(張雪東等，2010)，導致土水特征曲線無法由粒徑分布曲線唯一確定，采用粒徑分布曲線對土水特征曲線進行預測時，不可避免地存在模型預測誤差。

與通過非飽和土試驗獲得土水特征曲線相比，基于土體粒徑分布曲線估算土水特征曲線效率較高，但其伴隨的預測誤差不可忽視。如前所述，文獻中已提出了不少基于粒徑分布曲線的土水特征曲線預測方程，但極少對預測方程的誤差進行研究，導致基于粒徑分布曲線的土水特征曲線預測方程存在被濫用的危險。與已有研究相比，本文提出的方法不但能提供土水特征曲線參數的最優取值，還能給出土水特征曲線的變異范圍，為分析土水特征曲線預測誤差對非飽和巖土體性能評價的影響、評定預測誤差是否位于可接受范圍內、以及是否應采納其他更為精確的土水特征曲線測試和預測方法提供了基礎。例如，在邊坡穩定性分析中，如土水特征曲線變異范圍內邊坡始終處于穩定或不穩定狀態，表明土水特征曲線的預測誤差對邊坡是否需要進行加固的決策影響不大，該預測誤差可視為可接受；如土水特征曲線在其預測誤差范圍內變化時，邊坡是否需要加固的結論發生顯著變化，說明土水特征曲線預測誤差對工程決策具有重要影響，此時可考慮采用更精確的方法對土水特征曲線進行預測或測量。

土水特征曲線變異性對非飽和土工程系統的影響也可采用可靠度理論來進行定量分析。如可靠度分析表明土水特征曲線的不確定性不是巖土體性能預測的主導不確定性因素，則預測誤差可視為可接受；當土水特征曲線的不確定性對非飽和土系統的性能預測具有重大影響時，可考慮采用更精確的土水特征曲線測試或預測方法。

4 結 論

本文基于100組砂土的土水特征曲線和粒徑分布曲線試驗數據，分析了土水特征曲線VG模型參數與粒徑分布曲線參數的相關性，采用線性回歸方法獲得了基于粒徑分布曲線對VG模型參數的預測方程，通過殘差分析建立了土水特征曲線參數的概率預測模型，并基于實測數據對提出的模型進行了驗證。研究表明，基于粒徑分布曲線無法唯一確定土體的土水特征曲線，采用粒徑分布曲線對土水特征曲線進行預測時不可避免地存在模型誤差。與已有研究相比，本文提出的方法不但能提供土水特征曲線參數的最優取值，還能給出土水特征曲線的變異范圍，為明晰土水特征曲線預測誤差對非飽和土工程力學性能分析與評價的影響奠定了基礎。