張 潔 陽 帥 譚澤穎 李 旭

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

0 引 言

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

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

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

1 土水特征曲線模型

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

(1)

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

(2)

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

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

圖1 土水特征曲線擬合圖

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

2.1 土水特征曲線及粒徑分布曲線參數(shù)的相關(guān)性分析

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

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

圖2 土水特征曲線

圖3 土樣粒徑分布曲線

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

表1 VG模型參數(shù)相關(guān)性分析

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

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

Y=b0+b1X1+b2X2+b3X3+ε

(3)

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

μY=b0+b1X1+b2X2+b3X3

(4)

σY=σ

(5)

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

(6)

式中：φ為一標(biāo)準(zhǔn)正態(tài)分布的概率密度函數(shù)。假設(shè)存在s組土樣，令δ={δ1，δ2，…，δs}。當(dāng)B已知的條件下，觀測到所有數(shù)據(jù)的概率為：

(7)

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

采用相同的方法，也可對lnn的預(yù)測方程進行標(biāo)定。

2.3 VG模型參數(shù)的概率分布預(yù)測

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

(8)

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

圖4 預(yù)測值與實際值對比

采用最大似然法對lnn進行標(biāo)定，可得其均值預(yù)測公式為：

(9)

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

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

(10)

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

圖5 模型參數(shù)殘差直方圖

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

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

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

(11)

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

(12)

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

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

3 模型驗證

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

圖7 土樣粒徑分布曲線

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

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

圖9 驗證土樣土水特征曲線及預(yù)測曲線

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

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

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

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

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

(13)

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

圖10 土水特征曲線置信區(qū)間

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

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

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

4 結(jié) 論

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